A software engineer's circuitous journey to calculate eigenvalues
stylewarning.com
stylewarning.com
Assume you've got a routine for unitarily diagonalising a self-adjoint matrix with entries in the complex numbers. Call the resulting decomposition the spectral decomposition.
Now imagine keeping that same routine, but changing the complex numbers to other *-algebras. That means redefining what the operations {plus,times,minus,div,conjugate} do, while keeping the algorithm the same.
The observation is that for one choice of *-algebra, the spectral decomposition changes to the SVD, and an SVD-algorithm is obtained directly. For another *-algebra, the spectral decomposition changes to the Takagi decomposition, and a Takagi-algorithm is obtained directly.
The above uses operator overloading to turn one matrix decomposition into another.
https://github.com/adamnemecek/adjoint/
Check the raw source there’s a bunch of links.
Fixed points, diagonalizations and eigenshit are all the same thing.
This is true in a moral sense, but not literally. Fixed points are eigenvectors for a particular eigenvalue \lambda = 1 (together with the 0 vector, which is always fixed but not an eigenvector, by fiat). (Over a field,) eigenvectors corresponding to non-0 eigenvalues are not fixed points, although they lie on fixed lines; and they can be made fixed points by re-scaling the original matrix, if that is convenient. I don't know an obvious way to think of eigenvectors with eigenvalue \lambda = 0 as fixed vectors.
A matrix that is diagonalizable has sufficiently many eigenvectors, but a matrix may have eigenvectors even if it is not diagonalizable.
How do you think this sort of generalized approach would work if you want things to behave nicely with, say, poorly conditioned matrices, or matrices with certain structure that might get clobbered when changing to another algebra? Are these things even necessary to consider in such a design?
Or would this direction be mostly for bootstrapping a functional—but perhaps not performant or numerically stable—library?
Why? LAPACK is probably the most well tested bit of numerical code out there and must be a piece of the most well tested software in general. What is the technical justification except for some NIH?
The above quote is suggesting to minimize dependence on a mechanical translation of 15 year old code that's solely used as a fallback to the native (albeit foreign) routines. There is already a way to load a full, native LAPACK library in MAGICL (including netlib's reference, but also e.g. Apple's Accelerate framework), and that's always available precisely because of the quality of modern LAPACK distributions.
To be clear, the article discusses how ~10 lines of code (+ 3000 words of mathematical justification) can allow an existing LAPACK function that works on one data type to become useful for another data type. That hardly seems like NIH to me.
A more nuanced discussion is elsewhere in this thread.
If they wanted to minimize the use of "antiquated" Fortran code, that's fine I guess. But why didn't they just use ZGEEV for both the real and the complex case? Instead, they went on great lengths to transform complex matrices into real matrices, and recover the eigenvalues of the original complex matrix somehow...
It's educational, sure. But it also seems unnecessary.
As such, ZGEEV would require
(1) getting COMPLEX*16 values and arrays to translate appropriately into Common Lisp (where the representation is as a boxed, heap-allocated object, typically), and
(2) introducing approximately double the amount of code (totaling well over 100,000 lines of mechanically translated code... just to get one function!).
To (1), I think F2CL supports this, but it's evident from the longest-maintained distribution of the software that it hasn't been as thoroughly exercised as the real-type code. So, if LOC doesn't cost anything, this might be a fine starting point.
Lastly, "antiquated" specifically refers to a 15+-year-old LAPACK distribution written in FORTRAN 77, for which F2CL works—sort of. F2CL is actually more strict than F77, which makes translation of real-world F77 code difficult. For instance, in F77, you can pass an entire string to a function that declares itself as accepting only a single character. However, after translation, Lisp will raise a type error when that happens. So, our choice is to fix 30-year-old F77-to-Lisp translation code, or manually fix those Fortran type errors in the output. The latter is what was done to make the real code work.
Reference LAPACK is now written in Fortran 90+, with many numerical bug fixes in each release. As such, this cannot be translated with this tool, unless somebody writes a full-blown F90 parser and refactors the translator.
If one wants to use modern LAPACK directly in MAGICL, there's already a MAGICL/EXT-LAPACK system one can load, but it requires your system to have the binary library installed.
Don't roll your own lapack functions is generally good advice.
If your constrained to using Common Lisp, you can
- load and distribute LAPACK binaries (already well supported in principle, though Windows users typically have difficulties)
- fix the F77 translator and produce millions of lines of Common Lisp code for all LAPACK types (or, at worst, >100k for ZGEEV and dependencies)
- write a new F90+ translator and produce millions of lines of Common Lisp code for all LAPACK types
- write a numerically stable and accurate implementation from scratch that works for all manner of conditioned matrices (really hard!)
- translate a small subset of reliable F77 routines that clock in at a few 10k LOC and prove mathematically they can do more than what's written on the tin.
The first and last options are what have been chosen, the first allowing the most modern implementation of numerical functions for all data types, and the last allowing expedient independence from foreign binaries (i.e., it's pure Common Lisp).
Which option are you suggesting is optimal if you're seeking to relieve dependence—as a fallback for when it's unavailable—on a few shared libraries (/ foreign code)?
(It's probably out of the scope of this thread to get into the weeds of how Common Lisp is compiled as a standard language, and how specific compiler implementations like SBCL work. It suffices to say though that it's almost always not a simple matter of "get .o files and link them to your compiled Lisp", in no small part due to the popularity of image-based development.)
It's perfectly reasonable to consider different deployment strategies, but "statically link Fortran code (essentially into the Lisp runtime)", to a user of MAGICL, is not an ordinary thing to do.
Quelle surprise at the Windows issues, I'd probably look at using conda as a distribution tool to solve that (there are probably other options, but conda has the mindshare).
In Common Lisp specifically, you can effectively integrate with libraries (C ABI), but this heavily pivots on your definition of "effectively". One person's "effective" is another person's "impossible". For instance, for me, it is effective to
$ brew install lapack
$ sbcl
* (ql:quickload :magicl)
We can ask ourselves many questions about improving packaging and deployment, with different trade-offs. For instance, maybe we put in a lot of effort in MAGICL's build system to fetch hosted dependencies if they don't exist. But then some users wouldn't want that, so we might add a configuration option. But then some users might not discover that, so we should add a build-time error message that points the user to the configuration.
Alternatively, maybe we get Common Lisp's dependency manager (ASDF and Quicklisp) to integrate with something like Conda. That's possible if you take the time to implement it and persuade the respective maintainers it's a good idea. (It doesn't help if your customer doesn't use or doesn't want Conda, but c'est la vie.)
What used to be a "installing the dependencies is your problem; I'm going to keep my own library build simple" has turned into a complicated morass of configuration and communication. Some ecosystems, like Python's, have put in person-decades to make this easy, and they're still full of problems and confusion.
In this specific case, MAGICL is sometimes constrained to using portable Common Lisp because a customer does not want to (1) install conda on Windows and (2) wants to build the software from source. In fact, the customer isn't even using the linear algebra library directly, they're using an application which uses the library, which is why "just use numpy" isn't an option.
Imagine the pain once you find out some customers deploy this software onto air-gapped networks requiring all code that goes onto the network to be in text form and auditable. :) Now we are truly in for a doozy.
If your software is actually used, especially industrially, the "real world" sometimes just has really annoying corner cases, or users that have odd setups which you don't really empathize with. While seemingly obtuse (just install the damn dependency! it's one file!), using a little math to get something "for free" to get a larger program running without dependencies doesn't seem like such a bad deal.
Presumably if people want to build everything from scratch, then they know how to set up a correct build environment, and know that they need trustworthy libraries which are numerically correct, rather than a library which hasn't been as rigorously vetted which in a easier-to-install language. If not, I'd be curious what they would do with security/cryptographic libraries, would they be happy with someone's custom SHA256 library that they'd knocked up to avoid needing to install openssl for example.
I guess I expect people to treat numerical libraries in the same way as security/cryptographic libraries, as there's just as much subtlety in how those libraries need to be implemented.
> In addition to (and as part of) its support for multi-dimensional arrays, Julia provides native implementations of many common and useful linear algebra operations which can be loaded with using LinearAlgebra.
> [...]
> Linear algebra functions in Julia are largely implemented by calling functions from LAPACK.
While the performance of numerical code in Common Lisp is quite excellent due to very good support for fast floating point math (especially with SBCL's recent addition of SB-SIMD [2]), the library support isn't quite there like it is for Julia or Python. That's slowly changing, with MAGICL, numcl, Numericals, and Lisp-Stat.
A nice thing about Common Lisp compared to Julia or Python though is true AoT compilation without sacrificing a modicum of interactivity.