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