Hacker News new | past | comments | ask | show | jobs | submit login
A software engineer's circuitous journey to calculate eigenvalues (stylewarning.com)
109 points by reikonomusha on Aug 13, 2022 | hide | past | favorite | 27 comments



I've recently had the following idea that might make making new linear algebra libraries like this one easier. There might be some overlap with the ideas in the article, but generalised.

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.


A lot of the trouble of implementing numerical linear algebra is getting decent performance and numerical stability using finite-precision operations with bounded numerical representations.

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?


Yeah I have been realizing this lately as well. I have written up a bit on this

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.


> Fixed points, diagonalizations and eigen… 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.


Eigenvectors for eigenvalue 0 characterize the null space. The eigenvalue 0 can have a multiplicity of eigenvectors, these eigenvectors span the nullspace.


Interesting, I learned recently the reason that Galois fields are so important in error correction codes (and encryption) is that they enable invertable matrix operations using integers which makes them easy compute in dedicated hardware. Apparently development of the error correction codes used on CDs wouldn't have been feasible with the computing power of the time without being able to build custom silicon for it.


What I'm suggesting doesn't appear to be mentioned there.


The thing that makes it tick is adjoint and norm.


You might be interested in GraphBLAS: https://graphblas.org


Is this not basically what julialang does?


> Due to the age of the routines, I personally think it’s prudent to minimize its usage.

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?


It's netlib's FORTRAN 77 LAPACK distribution (over 15 years old), mechanically translated into Common Lisp. LAPACK isn't just one thing; there are multiple distributions of it, and the reference distribution is versioned, frequently updated, and no longer written in FORTRAN 77. So it's a bit erroneous to call LAPACK "well tested"; it's roughly like saying Python (2.5) is super well tested because half the programming world uses Python (3.x), Cython, Unladen Swallow, and/or Numba now.

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.


I don't get the whole premise of this article.

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...


Getting DGEEV to work requires 50,000 lines of mechanically translated code to work from HOMPACK, BLAS, and LAPACK. At no point is the "COMPLEX*16" data type used in the public interface to DGEEV because the routine demands real and imaginary parts are split into separate columns of the eigenvalue array and eigenvector matrix, i.e., actual complex numbers are never constructed. So things are "easy".

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.


But the lapack functions have been thoroughly tested, over decades.

Don't roll your own lapack functions is generally good advice.


LAPACK has been written in Fortran 90+ since 2008 in version 3.2. That's nearly 15 years ago. The last release of the LAPACK reference was 3.10.1 released April 2022.

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)?


Why would you be constrained to using Common Lisp? If you can't integrate with (effectively) non-reimplementable libraries (or would propose to rewrite crypto-code in CL?), than that to me says you shouldn't be using that language for that task?

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).


I'm not sure any language ecosystem, or operating system, has built a one-size-fits-all way of specifying dependencies.

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)
but it is not effective for a Windows user with a Visual C++ and commercially licensed Lisp implementation installed.

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.


Ah, I was arguing BLAS et. al. are effectively non-reimplementable (I'm aware of ATLAS and OpenBLAS, but I'd argue they demonstrate that it's a rare thing to be able to do). I would expect if people aren't willing to use existing tools like conda, then they understand they're going to need to put in the hard yards setting up a proper build environment.

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.


Use a Fortran compiler to compile the specific parts that you need, and link statically if that is what you want?


This is essentially the first option I listed, which is already supported by MAGICL by loading MAGICL/EXT-LAPACK [1]. Since MAGICL is a library and not a compiled application, using LAPACK in the way you suggest is still a foreign dependency (now a C compiler, Fortran compiler, Lisp compiler, and the transitive dependencies of the *GEEV Fortran functions), even if at the end of the day it's basically packaged up in an end-user's application binary.

(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.

[1] https://github.com/quil-lang/magicl#extensions


The eigenvalues and eigenvectors are going to be complex some of the time anyway. So I agree, why not just copy-and-paste the real algorithm (or use operator overloading) while remembering to insert complex-conjugates?

It's educational, sure. But it also seems unnecessary.


I think it's a very non-trivial task to "just" copy and paste a three-decade-old program that has tens of thousands of lines of dependencies, making sure to remember to add conjugates as necessary. :) Eigenvalues and eigenvectors are notoriously difficult programs to get right for general matrices with no known properties a priori.


Wonderful article, even if it’s not the best way to do things (per your conclusion and the discussion in this thread).


I don't know anything about common lisp, this seems to be a problem where Julia would shine marvelously. Given that most Julia users tend to be excellent mathematicians and physicists, I would imagine any advances or performance improvements in these types of calculations from the literature could be developed in Julia in a faction of the time that Fortran or C would take. In this case (numerics), the performance of Julia actually would probably be just as good as Fortran (or potentially faster due to being able to implement new ideas where performance enhancements quickly).


Julia indeed has great support for built-in linear algebra. They do lean on LAPACK (like just about any sane person/group would), but also re-implemented a lot of basic/common stuff natively in Julia. [1]

> 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.

[1] https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/

[2] http://www.sbcl.org/manual/index.html#sb_002dsimd


I'm hopeful that within the next year, StaticCompiler.jl gets good enough to aot compile lapack type code. we're close, but there's still a lot of rough edges.




Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact

Search: