eigvals doesn't accept a (Matrix, Symmetric) pencil, although eigen does

[jishnub]

julia> using LinearAlgebra

julia> A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
4×4 Matrix{Float64}:
 1.0  1.0  0.0  0.0
 1.0  2.0  1.0  0.0
 0.0  1.0  3.0  1.0
 0.0  0.0  1.0  4.0

julia> B = Matrix(Diagonal(Float64[1:4;]))
4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0

julia> eigen(A, Symmetric(B))
GeneralizedEigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
4-element Vector{Float64}:
 0.16959260913516316
 0.7541879474129647
 1.245812052587038
 1.830407390864838
vectors:
4×4 Matrix{Float64}:
 -0.59105     0.38815    -0.38815    0.59105
  0.490812   -0.0954118  -0.0954118  0.490812
 -0.224098   -0.341243    0.341243   0.224098
  0.0674664   0.347057    0.347057   0.0674664

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

Closest candidates are:
  eigvals!(::StridedMatrix{var"#s124"} where var"#s124"<:Union{Float32, Float64}; permute, scale, sortby)
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:306
  eigvals!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64}
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:578
  eigvals!(::Union{Hermitian{T, S}, Symmetric{T, S}}, ::Union{Hermitian{T, S}, Symmetric{T, S}}; sortby) where {T<:Union{Float32, Float64}, S<:(StridedMatrix{T} where T)}
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/symmetriceigen.jl:187

Stacktrace:
 [1] eigvals(A::Matrix{Float64}, B::Symmetric{Float64, Matrix{Float64}}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:622
 [2] eigvals(A::Matrix{Float64}, B::Symmetric{Float64, Matrix{Float64}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:620
 [3] top-level scope
   @ REPL[5]:1

julia> VERSION
v"1.10.0-DEV.1127"

It seems like this might benefit from the Cholesky factorization as well

[aravindh-krishnamoorthy]

The case eigvals!(A::StridedMatrix{T}, B::Union{LinearAlgebra.RealHermSymComplexHerm{T},Diagonal{T}} is indeed not defined. Here, we have a choice between these two definitions:

1. We can either map it to the generic version function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) in eigen.jl, which uses LAPACK's GGEV3.
2. Or, we could define the following function, based on cholesky, to handle it:

function LinearAlgebra.eigvals!(A::StridedMatrix{T}, B::Union{LinearAlgebra.RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
    return _choleigvals!(A, B, sortby)
end

function _choleigvals!(A, B, sortby)
    U = cholesky(B).U
    vals = eigvals!(LinearAlgebra.UtiAUi!(A, U))
    return vals
end

In the following, the benchmark results on my PC for the two methods above. First, for option number 1 based on the generic function in eigen.jl:

julia> @benchmark eigvals(A,B)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.219 μs … 169.958 μs  ┊ GC (min … max): 0.00% … 97.47%
 Time  (median):     1.280 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.346 μs ±   2.285 μs  ┊ GC (mean ± σ):  2.32% ±  1.38%

     ██▄▂
  ▄▆████▅▄▃▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▁▂▁▂▂▂▂▂▁▁▁▂▁▂▂▁▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂ ▃
  1.22 μs         Histogram: frequency by time        2.21 μs <

 Memory estimate: 1.00 KiB, allocs estimate: 9.

Next, for option number 2 with cholesky:

julia> @benchmark eigvals(A,Symmetric(B))
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.788 μs … 111.861 μs  ┊ GC (min … max): 0.00% … 95.43%
 Time  (median):     2.003 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.144 μs ±   2.383 μs  ┊ GC (mean ± σ):  2.34% ±  2.13%

      ▅█▃
  ▃▆████▆▅▄▄▄▄▅▇█▇██▆▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂ ▃
  1.79 μs         Histogram: frequency by time        3.37 μs <

 Memory estimate: 2.56 KiB, allocs estimate: 14.

It seems like we're better off using the generic function from eigen.jl for eigvals(A,Symmetric(B))!

Note: For completeness, Symmetric(B) takes about 100 ns and has negligible impact.

[aravindh-krishnamoorthy]

@jishnub An example implementation, based on the above discussion, is here: aravindh-krishnamoorthy/julia@efbeb98

Could you please check if it works for you now?

[jishnub]

I find that this trend changes for larger matrices (the break-even is somewhere around 1000x1000).

julia> using LinearAlgebra, BenchmarkTools

julia> import LinearAlgebra: BlasReal, BlasComplex, HermOrSym

julia> A = Matrix(Tridiagonal(ones(2999), collect(Float64.(1:3000)), ones(2999)));

julia> B = Symmetric(Matrix(Diagonal(Float64[1:3000;])));

julia> function _choleigvals!(A, B, sortby)
           U = cholesky(B).U
           vals = eigvals!(LinearAlgebra.UtiAUi!(A, U))
           return vals
       end
_choleigvals! (generic function with 1 method)

julia> function eigvalsc!(A::StridedMatrix{T}, B::Union{LinearAlgebra.RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
           return _choleigvals!(A, B, sortby)
       end
eigvalsc! (generic function with 1 method)

julia> function eigvalsm!(A::StridedMatrix{T}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
           return eigvals!(A, Matrix{T}(B); sortby)
       end
eigvalsm! (generic function with 1 method)

julia> @btime eigvalsm!($(copy(A)), $(copy(B)));
  31.037 s (20 allocations: 69.76 MiB)

julia> @btime eigvalsc!($(copy(A)), $(copy(B)));
  16.086 s (23 allocations: 69.63 MiB)

This would suggest using the Cholesky approach, as evaluation times are relevant mainly for large matrices.

[aravindh-krishnamoorthy]

Thank you for this info. Indeed, I can also reproduce this on my end and note that for 1000x1000 matrices, Cholesky decomposition based method is better.

However, there is still an issue with this approach, which I did not mention last time, and I still prefer to delegate the calls to GGEV3:

The function eigen, which is based on cholesky, claims to take Symmetric input, but Cholesky decomposition can only be computed for positive (semi) definite matrices. So, the following is a future issue, which will force us to change the implementation anyway:

julia> using LinearAlgebra

julia> A = [1 2; 2 1]
2×2 Matrix{Int64}:
 1  2
 2  1

julia> B = [2 1; 1 -2]
2×2 Matrix{Int64}:
 2   1
 1  -2

julia> eigen(A, Symmetric(B))
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
[...]

Perhaps this reveals to us that we need a positive definite matrices class as will in stdlib, which can really improve several other algorithms, such as inverse, as well.

Nevertheless, I have updated the comments in the commit. The relevant commits are as follows: aravindh-krishnamoorthy/julia@efbeb98 + aravindh-krishnamoorthy/julia@86e72bb

I will let this lie around for a bit to gather further comments and opinions...

[jishnub]

It should be possible to check if the Cholesky decomposition fails, and dispatch to ggev3 in that case. This should address all cases, I think? Admittedly, this comes at the expense of an extra attempt at a Cholesky decomposition in the non-positive-definite case, but this would be fast compared to ggev3, I would imagine.

[aravindh-krishnamoorthy]

@jishnub Thank you for your reply. But, I still feel that attempting cholesky on a Symmetric matrix is wrong, as cholesky is only defined for positive (semi) definite matrices. I would suggest to fix this issue (non-existent function for Matrix, Symmetric) with the fix above which uses the generic function (since LAPACK does no have a Matrix+Symmetric function), and open a new ticket to decide how to go ahead with positive (semi) definite matrices.

Would this be acceptable?

[jishnub]

The way to check if a matrix is positive-definite in Julia is to perform a Cholesky and see if it fails, so the implementation is actually in the other way. I understand your concerns, but given that a Cholesky decomposition provides a significant boost to a common real-life case, perhaps it's better to attempt it anyway? The cholesky function in LinearAlgebra accepts a check = false argument that may be used to check if this succeeds instead of throwing an error.

Are generalized eigenvalue problems with a non-positive-definite B commonly encountered? I'm not really familiar with them.

Wonder what @dkarrasch has to say about this?

[aravindh-krishnamoorthy]

Sure, we could wait for the comment of others. But, I do not think that attempting Cholesky decomposition on symmetric matrices is a good idea. I'd rather have a new positive (semi) definite (PSD) matrix type for such functions. That way, we could have the performance boost, and, at the same time, be consistent.

[jishnub]

Such types already exist in packages, e.g. https://github.com/JuliaStats/PDMats.jl, so this doesn't need to be added to LinearAlgebra. Using such a matrix type would automatically bypass this whole issue.

[dkarrasch]

Fixed by JuliaLang/julia#49673.