6 ms·

A software engineer's circuitous journey to calculate eigenvalues

109 points by reikonomusha 4y ago

ogogmad 4y ago
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.
- adamnemecek 4y ago
  Yeah I have been realizing this lately as well. I have written up a bit on this https://github.com/adamnemecek/adjoint/ 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.
  - ogogmad 4y ago
    What I'm suggesting doesn't appear to be mentioned there.
    - adamnemecek 4y ago
      The thing that makes it tick is adjoint and norm.
  - elcritch 4y ago
    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.
  - JadeNB 4y ago
    > 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.
    - haskellandchill 4y ago
      Eigenvectors for eigenvalue 0 characterize the null space. The eigenvalue 0 can have a multiplicity of eigenvectors, these eigenvectors span the nullspace.
- bigbillheck 4y ago
  You might be interested in GraphBLAS: https://graphblas.org https://graphblas.org
- throwawaymaths 4y ago
  Is this not basically what julialang does?
- reikonomusha 4y ago
  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?
aaaaaaaaaaab 4y ago
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...
- ogogmad 4y ago
  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.
  - reikonomusha 4y ago
    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.
- reikonomusha 4y ago
  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.
deleted 4y ago
[deleted]
chaxor 4y ago
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).
- reikonomusha 4y ago
  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/ https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/ [2] http://www.sbcl.org/manual/index.html#sb_002dsimd http://www.sbcl.org/manual/index.html#sb_002dsimd
  - adgjlsfhk1 4y ago
    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.
cycomanic 4y ago
> 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?
- deleted 4y ago
  [deleted]
- reikonomusha 4y ago
  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.
math-dev 4y ago
Wonderful article, even if it’s not the best way to do things (per your conclusion and the discussion in this thread).

109 points by reikonomusha 4y ago

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.

adamnemecek 4y ago
Yeah I have been realizing this lately as well. I have written up a bit on this https://github.com/adamnemecek/adjoint/ 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.
- ogogmad 4y ago
  What I'm suggesting doesn't appear to be mentioned there.
  - adamnemecek 4y ago
    The thing that makes it tick is adjoint and norm.
- elcritch 4y ago
  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.
- JadeNB 4y ago
  > 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.
  - haskellandchill 4y ago
    Eigenvectors for eigenvalue 0 characterize the null space. The eigenvalue 0 can have a multiplicity of eigenvectors, these eigenvectors span the nullspace.
bigbillheck 4y ago
You might be interested in GraphBLAS: https://graphblas.org https://graphblas.org
throwawaymaths 4y ago
Is this not basically what julialang does?
reikonomusha 4y ago
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?

aaaaaaaaaaab 4y ago

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

ogogmad 4y ago
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.
- reikonomusha 4y ago
  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.
reikonomusha 4y ago
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.

deleted 4y ago

[deleted]

chaxor 4y ago

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

reikonomusha 4y ago
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/ https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/ [2] http://www.sbcl.org/manual/index.html#sb_002dsimd http://www.sbcl.org/manual/index.html#sb_002dsimd
- adgjlsfhk1 4y ago
  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.

cycomanic 4y ago

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

deleted 4y ago
[deleted]
reikonomusha 4y ago
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.

math-dev 4y ago

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