linalg: Eigenvalues and Eigenvectors #816

perazz · 2024-05-16T08:28:27Z

Computing eigenvalues and eigenvectors of a square matrix: $A \cdot \bar{v} - \lambda \cdot \bar{v} = 0$ .
Based on LAPACK General (*GEEV) or real symmetric (*SYEV) / complex Hermitian (*HEEV) operations.

General matrix

Prior art

Numpy: eigenvalues, eigenvectors = eig(a)
Numpy: eigenvalues = eigvals(a)
Scipy: eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False, check_finite=True, homogeneous_eigvals=False)

Proposed implementation

Eigenvalues only: two function interfaces lambda = eigvals(A [, err])
Eigendecomposition: subroutine interface call eig(a,lambda [,right] [,left] [,overwrite_a] [,err])
- Eigenvalues of real matrices can be complex-conjugate. To avoid confusion, right, left are postprocessed and always returned as complex arrays. LAPACK instead stores (re,im) of conjugates as consecutive real values [j, j+1] in a real array, this would be very confusing and hard to track down.
- The eigendecomposition can be used either for eigenvalues only, or to retrieve left or right eigenvectors.

Real Symmetric / Complex Hermitian

Prior art

Numpy: eigenvalues, eigenvectors = eigh(a, uplo='L')
Numpy: eigenvalues = eigvalsh(a)
Scipy: eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=<object object>, eigvals=<object object>, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None)

Proposed implementation

Eigenvalues only: two function interfaces lambda = eigvalsh(A [, err])
Eigendecomposition: subroutine interface call eigh(a,lambda [,vectors] [,upper_a] [,overwrite_a] [,err] )
- option to access lower/upper triangle part only
- Avoid allocation when possible e.g. when providing vectors that can be used as temporary storage for a.

Note: the LAPACK backends are not pure due to functions with intent(out) arguments. The current proposed interfaces are ready to be made pure (e.g., function interface with all intent(in) arguments) to make it easy when the backend will be pure.

I believe this is ready for consideration, thank you @jvdp1 @jalvesz @fortran-lang/stdlib

jvdp1

Thank you. Here are some comments regarding the specs. I'll focus on the code later

doc/specs/stdlib_linalg.md

example/linalg/example_eig.f90

perazz · 2024-06-04T17:06:55Z

From Fortran Monthly call:

For general matrix, option to also return real array of singular values (it is complex by default). But, stop on errors if the matrix had complex eigenvalue pairs.

For example, NumPy returns real eigenvalues when all imaginary parts are zero:

    if not isComplexType(t) and all(w.imag == 0.0):
        w = w.real
        vt = vt.real
        result_t = _realType(result_t)
    else:
        result_t = _complexType(result_t)

Thank you @jvdp1 for the reviews. As pointed out during our last Monthly call with @everythingfunctional @jalvesz, I've added the option to return the real part of the eigenvalues only. In this latest case, we must check that there are no meaningful imaginary components, and raise an error if so:

stdlib/src/stdlib_linalg_eigenvalues.fypp

Lines 346 to 350 in e6aa0f5

    
           ! Check that no eigenvalues have meaningful imaginary part 
        
           if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then  
        
              err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & 
        
                                 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) 
        
           endif

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

jalvesz · 2024-06-08T13:11:39Z