[HN Gopher] A software engineer's circuitous journey to calculat... ___________________________________________________________________ A software engineer's circuitous journey to calculate eigenvalues Author : reikonomusha Score : 28 points Date : 2022-08-11 16:34 UTC (2 days ago) (HTM) web link (www.stylewarning.com) (TXT) w3m dump (www.stylewarning.com) | ogogmad wrote: | 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 wrote: | 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. | ogogmad wrote: | What I'm suggesting doesn't appear to be mentioned there. | adamnemecek wrote: | The thing that makes it tick is adjoint and norm. | aaaaaaaaaaab wrote: | 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 wrote: | 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 wrote: | 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 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!). | | Lastly, "antiquated" refers to a LAPACK 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 95+, with voluminous | bug fixes in each release. As such, this cannot be translated | with this tool, unless somebody writes a full-blown F95 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. | amelius wrote: | But the lapack functions have been thoroughly tested, over | decades. | | Don't roll your own lapack functions is generally good | advice. | reikonomusha wrote: | 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 release April 2022. ___________________________________________________________________ (page generated 2022-08-13 23:00 UTC)