Document status and goals of high-level library (was: What is the status of the high-level interface?)

As LLA is now an abandonedware, I am looking into Magicl as a potential back-end of numcl. I would like to be briefed about the library in order to make myself 100% sure about the decision, especially the state of its high-level interface.

I hope there is not much overlap between numcl and magicl yet --- numcl is still focusing on the high-level interface, and I believe magicl is still more on the low-level BLAS/LAPACK side. I only see einsum and couple of other things in the src/high-level/ directory, and hope not many libraries are depending on it yet.
Could you brief me about the datastructure mainly passed around by magicl functions. numcl works on the natural CL simple 1D arrays in order to make future migration to C API easy.
Handling of element-types. I know LLA handles element-type nicely, dispatches to the right functions and returns correct arrays. Is it handled by magicl already?
Requirement to the pointer pinning and word alignment. Some Lisp allows to pin a pointer so that GC does not move them for some CFFI applications. Does magicl have any particular requirement regarding the raw pointers?

Hey @guicho271828:

As LLA is now an abandonedware, I am looking into Magicl as a potential back-end of numcl. I would like to be briefed about the library in order to make myself 100% sure about the decision, especially the state of its high-level interface.
* I hope there is not much overlap between numcl and magicl yet --- numcl is still focusing on the high-level interface, and I believe magicl is still more on the low-level BLAS/LAPACK side. I only see einsum and couple of other things in the `src/high-level/` directory, and hope not many libraries are depending on it yet.

I doubt much depends on MAGICL, only because it's not such a friendly or generic library to use.

The high-level interface was supposed to be a one-stop shop for numerical linear algebra, with nice CLOS/structure objects for each kind of array, handling polymorphism and "do what I mean" semantics. It was meant to support dense tensors, sparse tensors, and the specialties of BLAS/LAPACK, like banded data, strided data, etc.

The unfortunate reality, however, is that we got the library "good enough" to do efficient matrix arithmetic over (complex double-float)s, we got the special functions we wanted (matrix exponentiation and logarithm, QR, SVD, zuncsd, etc.), and that allowed us to move forward with our work in quilc and qvm (other projects we work on). We haven't, unfortunately, had the time to further invest in the high-level interface, although it has always been on the backburner.

We have, in the back of our minds, hoped that some very talented and interested individual would look at the MAGICL interface, see that the high-level is a wreck, and say "I can do better." But alas, it hasn't happened. (Though I give special mention to @kilimanjaro who provided some cool macros for doing einsum in Lisp, which you've seen.) Writing a linear algebra library is thankless and difficult work, and is mostly an exercise in good UI design, plumbing, etc.

We currently have no plans during "work hours" to improve MAGICL, though things we've wanted in a high-level interface are:

Support for single/double-precision real and complex numbers in all functions. Many previous libraries have utterly lacked the ability to work with complex numbers.
Easy to use, easy to optimize, easy to access unsafe internals.
Support only linear algebra, not scientific computing more broadly.
Support for different storage layouts, like dense/sparse/banded/diagonal/strided/etc.

We have not had a desire to be compatible with CL arrays. The CL standard library for N-dimension arrays is lacking anyway.

In short, we have been happy to throw away the high-level interface and start from scratch, provided we can continue to service quilc and qvm.

* Could you brief me about the datastructure mainly passed around by magicl functions. numcl works on the natural CL simple 1D arrays in order to make future migration to C API easy.

In the high-level interface, there is a 2D matrix structure which is backed by a 1D specialized simple array. All of the low-level functions take column major data (as required by Fortran!). This is in direct conflict with CL's usual row-major array storage. No special data structure is used at the low level; we just use these specialized vectors, finding their underlying storage pointer, pinning that storage, and calling into Fortran. They must be one of these specialized arrays:

(simple-array single-float (*))
(simple-array double-float (*))
(simple-array (complex single-float) (*))
(simple-array (complex double-float) (*))

Note that some Lisp implementations, like LispWorks and (still in development!) ECL have not had support for specialized complex arrays.

* Handling of element-types. I know LLA handles element-type nicely, dispatches to the right functions and returns correct arrays. Is it handled by magicl already?

No dispatching exists in MAGICL anywhere. There aren't even good ways to construct and operate on arrays polymorphically at all. There's no generic arithmetic, or anything. You can see damning functions like multiply-complex-matrices, which works for one and only one type.

* Requirement to the pointer pinning and word alignment. Some Lisp allows to pin a pointer so that GC does not move them for some CFFI applications. Does magicl have any particular requirement regarding the raw pointers?

Yes, the data must be pinned and the data must be laid out as usual (sequential 32/64 bit single/double floats). Specifically, the underlying storage vector—a Lisp object—must be pinned. Once that is pinned, you can make the call into Fortran, then un-pin it. That code lives here.

I hope this helps, and we are definitely happy to see the Lisp numerical ecosystem improve.

Thank you for a lot of clarifications! By the way, these clarifications I believe can be put in somewhere in the repository as part of the documentation (with some format/style cleanup, of course.)

I'm going to keep this issue open as a reminder as to documenting the goals of the lib.

I just came across this issue today!

I have been intermittently working on a few projects this past 1-1.5 year that might address some of the issues mentioned above. But

I lack in experience in both industry terms as well as with magicl or linear algebra in general. My experience is primarily from taking a few machine learning courses online and during undergrad studies.
I might have (hopefully only slightly!) differing high level goals.

Support for different storage layouts, like dense/sparse/banded/diagonal/strided/etc.

I specifically wanted numpy-like multidimensional strided arrays for copy-free slicing - the result is dense-arrays. But then, I felt this might be abstracted out into a separate abstract-arrays system, which people can then use to implement sparse-arrays and more.

Even here, I am yet to take a look at abstract-classes. But the custom-arrays specializing on element-types and ranks using cl:satisfies type seems to be the limit of CL type system. Beyond this, I suspect one needs a custom type system, or coalton.

No dispatching exists in MAGICL anywhere.

I also invested some time in polymorphic-functions primarily for dispatching on specialized array types - CL's generic functions not providing for this. Examples of dispatching is at two-arg-fn-float.lisp and two-arg-fn-all.lisp both in numericals. However, this too is beyond ANSI CL.

Beyond this, I suspect a limiting factor would indeed be how much of this can be done without breaking anything in quilc, qvm, or most-anything plenty of other users are using this library for.

A lot has changed about MAGICL since 2019. The general sentiment of "we got it to a point where it works for us" is still true, and thus it's still not a one-stop shop for all things linear algebra, however, the interface has matured.

I'd also say that compared to my comment in 2019, there are now more people who are regularly working with and contributing back to MAGICL. Nobody is hacking on it full-time, but it is receiving more attention. We even have meetings every 2 weeks on the Quil-Lang Discord, among other things.

As described in the README, there is a dispatching mechanism now for backend function implementations, as well as dispatching on tensor type. We decided to forego polymorphic types altogether, and just encode the float type in the class name of the tensor (e.g. MATRIX/DOUBLE-FLOAT). There has also been work to not require BLAS/LAPACK, since it's a heavy dependency. There's also been porting effort: MAGICL now is known to work on SBCL, Allegro, and LispWorks.

The high-level interface, by and large, attempts to be as plain as possible. The dispatching/backend stuff felt like a necessary evil, but we're not trying to be any more creative than that, despite being called "MAGICL". At the end of the day, most MAGICL users want a reliable, efficient, predictable, and direct way to do computational linear algebra without too much typing or voodoo.

We've inched forward a bit on the documentation front. The README is a little better, and there's some information in the doc directory. But still, I think there's plenty of work to do here.

We've inched forward a bit on the documentation front. The README is a little better, and there's some information in the doc directory. But still, I think there's plenty of work to do here.

Playing around with MAGICL, this definitely seems to be the case. There is plenty of magic here; but more spells seem to be needed to be added into the spell-book.

~~Also, I attempted a rudimentary backend for dense-arrays using magicl: https://github.com/digikar99/dense-arrays/commit/77bc4b5eb6204e2e22a721b580f04a644922dbbe~~ May 29, 2022, EDIT: After noting magicl::storage is a cl:vector, I no longer need a separate class, I can simply work with standard-dense-array. What this also means is that the functionality is dense-numericals can now be used alongside magicl!

A conversion step is certainly required, but once that is done, the performance isn't actually bad. ~~Admittedly, I still need to work on the layout (besides bug fixes themselves), without which the magicl backend might not be of much use, but so, we already have another frontend.~~ Layout should also be ready now as of this commit; bugs may still persist though!

DENSE-ARRAYS-PLUS-LITE> (let* ((a (make-array '(100 100) :element-type 'double-float
                                                         :initial-element 0.0d0))
                               (b 0.0d0))
                          (declare (optimize speed)
                                   (type (array double-float) a)
                                   (type double-float b))
                          (time (loop repeat 1000
                                      do (do-arrays ((a-elt a double-float))
                                           (declare (type double-float b a-elt))
                                           (incf b a-elt))))
                          b)
Evaluation took:
  0.552 seconds of real time
  0.551673 seconds of total run time (0.551673 user, 0.000000 system)
  100.00% CPU
  1,218,525,560 processor cycles
  320,045,040 bytes consed

0.0d0
DENSE-ARRAYS-PLUS-LITE> (let* ((a (magicl:zeros '(100 100) :type 'double-float))
                               (b 0.0d0))
                          (declare (optimize speed)
                                   (type double-float b))
                          (time (loop repeat 1000
                                      do (loop :for i :below 100
                                               :do (loop :for j :below 100
                                                         :do (incf b (magicl:tref a i j))))))
                          b)
Evaluation took:
  1.200 seconds of real time
  1.200351 seconds of total run time (1.200351 user, 0.000000 system)
  100.00% CPU
  2,650,436,110 processor cycles
  320,077,824 bytes consed

0.0d0
DENSE-ARRAYS-PLUS-LITE> (let* ((a (magicl:zeros (list 1000 1000) :type 'double-float))
                               (b (magicl:zeros (list 1000 1000) :type 'double-float)))
                          (declare (optimize speed))
                          (magicl.backends:with-backends (:blas)
                            (time (loop repeat 100 do (magicl:.+ a b a)))))
Evaluation took:
  1.035 seconds of real time
  11.865939 seconds of total run time (4.808226 user, 7.057713 system)
  [ Run times consist of 0.343 seconds GC time, and 11.523 seconds non-GC time. ]
  1146.47% CPU
  2,279,682,150 processor cycles
  800,001,600 bytes consed

NIL
DENSE-ARRAYS-PLUS-LITE> (let* ((a (zeros (list 1000 1000) :type 'double-float))
                               (b (zeros (list 1000 1000) :type 'double-float)))
                          (declare (optimize speed)
                                   (type (magicl-array double-float) a b))
                          (magicl.backends:with-backends (:blas)
                            (time (loop repeat 100 do (magicl-funcall #'magicl:.+ a b a)))))
Evaluation took:
  0.448 seconds of real time
  5.097049 seconds of total run time (1.922890 user, 3.174159 system)
  1137.72% CPU
  989,807,680 processor cycles
  131,040 bytes consed

NIL
DENSE-ARRAYS-PLUS-LITE> (let* ((a (magicl:zeros (list 1000) :type 'double-float))
                               (b (magicl:zeros (list 1000) :type 'double-float)))
                          (declare (optimize speed))
                          (magicl.backends:with-backends (:blas)
                            (time (loop repeat 100000 do (magicl:.+ a b a)))))
Evaluation took:
  0.688 seconds of real time
  0.687958 seconds of total run time (0.687958 user, 0.000000 system)
  100.00% CPU
  1,519,060,902 processor cycles
  812,798,608 bytes consed

NIL
DENSE-ARRAYS-PLUS-LITE> (let* ((a (zeros (list 1000) :type 'double-float))
                               (b (zeros (list 1000) :type 'double-float)))
                          (declare (optimize speed)
                                   (type (magicl-array double-float) a b))
                          (magicl.backends:with-backends (:blas)
                            (time (loop repeat 100000 do (magicl-funcall #'magicl:.+ a b a)))))
Evaluation took:
  1.688 seconds of real time
  1.689025 seconds of total run time (1.689025 user, 0.000000 system)
  100.06% CPU
  3,729,354,166 processor cycles
  102,400,000 bytes consed

NIL
DENSE-ARRAYS-PLUS-LITE> (let ((a (rand 3 3 :type 'single-float)))
                          (print a)
                          (magicl-funcall #'magicl:svd a))

#<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3 SINGLE-FLOAT
   (  0.102       0.590       0.073    )
   (  0.374       0.817       0.980    )
   (  0.916       0.891       0.975    )
 {105E99A073}> 
#<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3 SINGLE-FLOAT
   ( -0.218      -0.924      -0.314    )
   ( -0.619      -0.118       0.776    )
   ( -0.755       0.363      -0.547    )
 {105E9A1D03}>
#<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3 SINGLE-FLOAT
   (  2.108       0.000e+00   0.000e+00)
   (  0.000e+00   0.410       0.000e+00)
   (  0.000e+00   0.000e+00   0.320    )
 {105E9A1DD3}>
#<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3 SINGLE-FLOAT
   ( -0.448      -0.620      -0.644    )
   (  0.473      -0.776       0.418    )
   ( -0.759      -0.117       0.641    )
 {105E9A1EA3}>
DENSE-ARRAYS-PLUS-LITE> (let ((a (rand 3 3 :type 'single-float)))
                          (print a)
                          (multiple-value-bind (u s vt) (magicl-funcall #'magicl:svd a)
                            (magicl-funcall #'magicl:@ u s vt)))

#<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3 SINGLE-FLOAT
   (  0.293       0.588       0.419    )
   (  0.842       0.509       0.615    )
   (  0.552       0.412       0.359    )
 {10637DE403}> 
#<STANDARD-DENSE-ARRAY :COLUMN-MAJOR 3x3 SINGLE-FLOAT
   (  0.293       0.588       0.419    )
   (  0.842       0.509       0.615    )
   (  0.552       0.412       0.359    )
 {1063950B53}>
DENSE-ARRAYS-PLUS-LITE> (let ((a (aref (rand 3 3 3 :type 'single-float) nil nil 0)))
                          (print a)
                          (describe a)
                          (multiple-value-bind (u s vt) (magicl-funcall #'magicl:svd a)
                            (magicl-funcall #'magicl:@ u s vt)))

#<STANDARD-DENSE-ARRAY NIL 3x3 SINGLE-FLOAT
   (  0.480       0.506       0.480    )
   (  0.963       0.805       0.096    )
   (  0.979       0.775       0.738    )
 {1066F8E873}> 
#<STANDARD-DENSE-ARRAY NIL 3x3 SINGLE-FLOAT {1066F8E873}>
  [standard-object]

Slots with :INSTANCE allocation:
  STORAGE                        = #(0.4798323 0.080942154 0.24663782 0.50591195 0.3678155 0.009541988..
  DIMENSIONS                     = (3 3)
  ELEMENT-TYPE                   = SINGLE-FLOAT
  RANK                           = 2
  TOTAL-SIZE                     = 9
  STRIDES                        = (9 3)
  OFFSETS                        = (0 0)
  LAYOUT                         = NIL
  ROOT-ARRAY                     = #<STANDARD-DENSE-ARRAY :ROW-MAJOR 3x3x3 SINGLE-FLOAT {1066F8D253}>
#<STANDARD-DENSE-ARRAY :COLUMN-MAJOR 3x3 SINGLE-FLOAT
   (  0.480       0.506       0.480    )
   (  0.963       0.805       0.096    )
   (  0.979       0.775       0.738    )
 {1066F9C633}>

Any plan on CUDA integration?

Any plan on CUDA integration?

It would be cool but I don't have any hardware. So no plans from me at least.

With dense-arrays, there is a cl-cuda backend, and when I started out with dense-numericals, I also wanted to add support for cuda using cl-cuda. The former should be in a more working condition than the latter, while the latter needs to be updated and made a part of numericals. But the focus these days is more on getting more functions up to speed for the CPU-only backend.

And so, yes, I do have CUDA hardware, so if there are more users excited about CUDA backend for magicl, I'd want to try contributing.

PS: Is there a serious use case even for something like "distributed arrays" aka arrays that are spread over multiple machines over the network?

PS: Is there a serious use case even for something like "distributed arrays" aka arrays that are spread over multiple machines over the network?

The DQVM project (part of the QVM repo of this org) distributes a very large vector across machines. I couldn't imagine using a library for it though since it's such bespoke code, but it definitely is a distributed numerical object.

quil-lang / magicl

Document status and goals of high-level library (was: What is the status of the high-level interface?) #50