Generalized eigenvalue decomposition

[daanhb]

The generalized eigenvalue decomposition is not implemented in GenericLinearAlgebra:

julia> using GenericLinearAlgebra, DoubleFloats, LinearAlgebra

julia> A = rand(Double64, 5, 5); B = rand(Double64, 5, 5);

julia> eigvals(A, B)
ERROR: MethodError: no method matching eigvals!(::Matrix{Double64}, ::Matrix{Double64})

I am working on a pull request to implement this. It is based on code kindly provided to me by @thijssteel, who did his PhD on the topic. It uses the single shift or double shift QZ algorithm. The PR currently implements eigvals(A, B):

julia> using GenericLinearAlgebra

julia> eigvals(A, B)
5-element Vector{Complex{Double64}}:
    -7.55847194974678743120673992727391194 + 0.0im
 -7.6691169090262111372117306819298531e-01 + 0.0im
 9.77704965303924447091419817528601559e-03 - 6.82491740678931004263096986923651618e-01im
 9.77704965303924447091419817528601559e-03 + 6.82491740678931004263096986923651618e-01im
     2.67081937010893995823217172351462554 + 0.0im

The speed and accuracy of the implementation seem somewhat comparable to that of LinearAlgebra for Float64:

julia> n = 150; A = rand(n, n); B = rand(n, n);

julia> @time v1 = eigvals(A, B);
  0.014203 seconds (16 allocations: 411.875 KiB)

julia> @time v2 = GenericLinearAlgebra.generalized_eigvals(A, B);
  0.017426 seconds (5.10 k allocations: 754.531 KiB)

julia> norm(abs.(v1)-abs.(v2))   # abs to account for different ordering of complex conjugate pairs
8.884694390317172e-13

Further improvements may be possible, as the code does not include any of the optimizations that @thijssteel worked on (and which are made available in recent versions of LAPACK).

I have a use case for generalized eigenvalues as part of the AAA algorithm for computing rational approximations. It would be nice if it works with GenericLinearAlgebra.

[daanhb]

I haven't followed up on the PR for generalized eigenvalues, apologies, but I do use that code from time to time. In order to make it available (to make experiments in papers reproducible) and to speed up development, I have moved the code to a separate package: https://github.com/daanhb/GeneralizedEigenvalues.jl. It is really just the qz.jl file of this pull request.

My intention is that this is temporary. I'll add tests in the new repo and once I'm confident about the quality I'll make a new and more polished PR. So for now I'll close the PR, but I'd like to leave the issue open.

[NAThompson]

Note that the generalized eigenvalue problem $Ax = \lambda Bx$ is solved in this package, at least for $B \succ 0$. Example code:

using GenericLinearAlgebra: eigen, cholesky, Diagonal, Symmetric


A = fill(big(NaN), 2, 2)
A[1,1] = 6
A[1,2] = 2
A[2,1] = 2
A[2,2] = 3

B = fill(big(NaN), 2,2)
B[1,1] = 3
B[1,2] = 1
B[2,1] = 1
B[2,2] = 1

A = Symmetric(A)
B = Symmetric(B)
F = eigen(A, cholesky(B))
println(F)
R = A*F.vectors - B*F.vectors*Diagonal(F.values)
println("Residual = $R")

Output:

LinearAlgebra.GeneralizedEigen{BigFloat, BigFloat, Matrix{BigFloat}, Vector{BigFloat}}(BigFloat[2.0, 3.5], BigFloat[0.5773502691896257645091487805019574556476017512701268760186023264839776723029325 -0.4082482904638630163662140124509818986609912467761116880721154278751600629095498; 0.0 1.224744871391589049098642037352945695982973740328335064216346283625480188728649])
Residual = BigFloat[0.0 0.0; 0.0 0.0]

Note that unlike the eigen function defined in LinearAlgebra, we must explicitly pass in the Cholesky decomposition of B into the eigen call. (Question to the maintainer: Would a PR making the eigen call just call cholesky on B be welcome?)

It appears this is undocumented, so if I'm doing anything wrong here I'd appreciate corrections. But to the best of my understanding now, this "just works".