lisp-stat / linear-algebra Goto Github PK

One of the things we're going to need is a good pattern for dispatching to various backends. Here's one possible pattern that I used for a scatterplot smoother:

 ;; What about the case where both lla and linear-algebra are loaded?
		    beta #+lla (lla:solve A b)
			 #+linear-algebra (linear-algebra:solve A b)
		         #-(or lla linear-algebra) (solve-matrix (aops:stack-cols A b))

As I start to use it, I think it's more of an anti-pattern because:

• It doesn't easily handle the case where I might want multiple back-ends in one application (is this an edge case?)
• It is only evaluated at compile-time, so if I want to change I need to recompile
• It sort of abuses *features* in a non-specific way I don't like

I know that MGL-MAT has a way to automatically dispatch methods depending on whether BLAS or CUDA are available. Generally his architecture seems similar to that of linear-algebra, with various 'kernels' (lisp, CUDA, BLAS) that are operated on by lisp code. It might be worth studying.

linear-algebra has a dependency on floating-point, where float-equal has an overlap in functionality with numerical-utilities num=.

Looking at the code,%float-equal is much like num= with *epsilon* replaced by *tolerance*, with the exception that num= works on vectors and arrays. However the definitions of error are different. In num=:

(<= (abs (- a b)) (* (max 1 (abs a) (abs b)) tolerance))

and in %float-equal:

(defun %relative-error (exact approximate)
  "Return the relative error of the numbers."
  (abs
   (if (or (zerop exact) (zerop approximate))
       (- exact approximate)
       (/ (- exact approximate) exact))))

Perhaps someone with more recent experience here can weigh in on whether or not these are actually equivalent? By my reading they are for all practical purposes.

Advantages, as I see them for num=:

• Works on rational types as well as vectors and arrays
• Tested rather extensively during the development of the special-functions system
• Requires less rework; no one is currently using floating-point

Advantages for floating-point:

• Documented reference Numerical Algorithms with C

Finally, looking at the exports for float-point:

(defpackage :floating-point
  (:use :common-lisp)
  (:nicknames :fp)
  ;; Control parameters
  (:export :*epsilon* :*significant-figures*)
  ;; Constants
  (:export :pi/2)
  ;; Error analysis
  (:export :default-epsilon
           :relative-error)
  ;; Floating point predicates
  (:export :float-equal
	   :sigfig-equal))

I can't see anything there that num= doesn't give us, so I'm inclined to suggest we archive floating-point and use num= as a 'standard' for linear-algebra. There's a few bits and pieces that would be nice to move to numerical-utilities, like constants.lisp and the epsilon stuff, but none of that is critical. If it turns out something is missing, then we can move it to num= right away.

Thoughts? Comments? Did I miss something?

unary-operations implements the norm functions, but in a rather clumsy way. Refactor this using existing utilities.

When doing a fresh compile, 7 style warnings are generated:

 caught STYLE-WARNING:
;   The binding of SB-C::Y is not a REAL:
;    NIL
;   See also:
;     The SBCL Manual, Node "Handling of Types"

They don't affect the test results, but it would be nice to eliminate them.

lisp-unit is abandoned (?) and we should consider moving it to something like cl-unit2, both as an exercise to learn linear-algebra better, but also to bring all the dependencies to supported versions and, finally, to 'standardise' on a single test framework throughout the ecosystem (the Lisp-Stat ecosystem).

As far as I can tell this is a search-and-replace operation, as the systems (lisp-unit and cl-unit) are rather similar, and at the same time standardising on the numeric equality predicate (see issue #3).

Commit 8c9271f re-enabled tests for triangular-matrix types, and exposed a bug in make-matrix. Here is one of the failing tests:

  (let ((matrix (make-matrix 10 10 :matrix-type 'upper-triangular-matrix :initial-element 1.0)))
    (assert-true (matrixp matrix))
    (assert-true (typep matrix 'upper-triangular-matrix))
    (assert-num= matrix (upper-triangular-array))) ; failure is here

Looking at the compared values:

LS-USER> (linear-algebra-test::upper-triangular-array)
#2A((1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0)
    (0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0))

That is a correct upper-triangular-matrix.

LS-USER> (linear-algebra:make-matrix 10 10 :matrix-type 'linear-algebra:upper-triangular-matrix :initial-element 1.0)
#<LINEAR-ALGEBRA:UPPER-TRIANGULAR-MATRIX {102D7167C3}>
LS-USER> (linear-algebra::contents *)
#2A((1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0)
    (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0))

but make-matrix returns an array of all ones.

@beberman, is this something you can look into?

This consists of:

1. Developing a high level generic API based on numpy's linear algebra operations. Best location is probably array-operations
2. Implement them in lisp by refactoring or replacing the existing functions

Then implement the same in LLA so the API's match.

If we are going to refactor to use the builtin vector method then I would consider doing a full refactor of all the underlying classes in this system.

The system should support arbitrary tensors of any number of dimensions if you want to make this useful for more general computations.

The system should also support a wider class of specific matrix types - in particular diagonal, tridiagonal, upper tridiagonal etc.

My suggestion would be that the underlying type be based on vector for storage, a permutation operator that can be bound, and a ref that takes an index list.

So
(defclass tensor ()
((data - as a vector)
(dimensions - list)
(permutation - a list of integers or nil)))

Then a core reference function
(defmethod ref ((index ref) (tensor tensor))
which checks if permutation is nil/ non nil. If nil it just computes the required offset into the data given the dimensions and the index

If non-nil it first applies the permutation to the index then computes the offset.

For specialized matrixes we override this function and return 0 if the index is outside the range.

So then matrix is just
(defclass matrix (tensor)
)

with
(defgeneric make-matrix (&key &otherkeys)
(:method (&key row column) building a default)
(:method (&key row column type) building a specific type))

(defmethod transpose ((m matrix))
; this just sets the permutation to nil if it exists otherwise it sets it to '(1 0))
)

Then the generic matrix multiplies will work with all types if it uses ref underneath and specialized ones can be written to minimize the lookups as needed.

The system includes several specialised matrix types:

• hermitian
• square
• symmetric
• triangular

There are also existing versions of these matrix types in numerical-utilities. IMO, the implementation in numerical-utilities far outclass those in linear-algebra. The only thing missing from them is a lisp-based solver, which linear-algebra does have. There are LAPACK solvers in LLA, for those specialised matrices though.

This issue is to track a review of both of these and (probably) figure out a way to implement the solve function on the matrix implementations in numerical-utilities. An alternative may be to say that we're not going to implement a lisp based solver. I have used the gauss solver in an implementation of lowess, and find the speed very slow. This might be the reason that a lisp-based solver was never implemented for the specialised types in numerical-utilities.