Rework symmetric generalized eigen/eigvals (#49673)

Author: aravindh-krishnamoorthy

NEWS.md

@@ -92,6 +92,11 @@ Standard library changes
   (real symmetric) part of a matrix ([#31836]).
 * The `norm` of the adjoint or transpose of an `AbstractMatrix` now returns the norm of the
   parent matrix by default, matching the current behaviour for `AbstractVector`s ([#49020]).
+* `eigen(A, B)` and `eigvals(A, B)`, where one of `A` or `B` is symmetric or Hermitian,
+  are now fully supported ([#49533])
+* `eigvals/eigen(A, cholesky(B))` now computes the generalized eigenvalues (`eigen`: and eigenvectors)
+  of `A` and `B` via Cholesky decomposition for positive definite `B`. Note: The second argument is
+  the output of `cholesky`.
 
 #### Printf
 * Format specifiers now support dynamic width and precision, e.g. `%*s` and `%*.*g` ([#40105]).

stdlib/LinearAlgebra/src/diagonal.jl

@@ -796,12 +796,11 @@ function eigen(A::AbstractMatrix, D::Diagonal; sortby::Union{Function,Nothing}=n
     end
     if size(A, 1) == size(A, 2) && isdiag(A)
         return eigen(Diagonal(A), D; sortby)
-    elseif ishermitian(A)
+    elseif all(isposdef, D.diag)
         S = promote_type(eigtype(eltype(A)), eltype(D))
-        return eigen!(eigencopy_oftype(Hermitian(A), S), Diagonal{S}(D); sortby)
+        return eigen(A, cholesky(Diagonal{S}(D)); sortby)
     else
-        S = promote_type(eigtype(eltype(A)), eltype(D))
-        return eigen!(eigencopy_oftype(A, S), Diagonal{S}(D); sortby)
+        return eigen!(D \ A; sortby)
     end
 end
 

stdlib/LinearAlgebra/src/eigen.jl

@@ -524,7 +524,7 @@ true
 """
 function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA), TB)
-    eigen!(eigencopy_oftype(A, S), eigencopy_oftype(B, S); kws...)
+    eigen!(copy_similar(A, S), copy_similar(B, S); kws...)
 end
 eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
 
@@ -619,7 +619,7 @@ julia> eigvals(A,B)
 """
 function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA), TB)
-    return eigvals!(eigencopy_oftype(A, S), eigencopy_oftype(B, S); kws...)
+    return eigvals!(copy_similar(A, S), copy_similar(B, S); kws...)
 end
 
 """

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -156,6 +156,11 @@ end
 eigmax(A::RealHermSymComplexHerm{<:Real}) = eigvals(A, size(A, 1):size(A, 1))[1]
 eigmin(A::RealHermSymComplexHerm{<:Real}) = eigvals(A, 1:1)[1]
 
+function eigen(A::HermOrSym{TA}, B::HermOrSym{TB}; kws...) where {TA,TB}
+    S = promote_type(eigtype(TA), TB)
+    return eigen!(eigencopy_oftype{S}(A), eigencopy_oftype(B, S); kws...)
+end
+
 function eigen!(A::HermOrSym{T,S}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
     vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
@@ -164,26 +169,32 @@ function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,Not
     vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
-function eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    return _choleigen!(A, B, sortby)
-end
-function eigen!(A::StridedMatrix{T}, B::Union{RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    return _choleigen!(A, B, sortby)
+
+function eigen(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=nothing)
+    if ishermitian(A)
+        eigen!(eigencopy_oftype(Hermitian(A), eigtype(eltype(A))), C; sortby)
+    else
+        eigen!(copy_similar(A, eigtype(eltype(A))), C; sortby)
+    end
 end
-function _choleigen!(A, B, sortby)
-    U = cholesky(B).U
-    vals, w = eigen!(UtiAUi!(A, U))
-    vecs = U \ w
+function eigen!(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=nothing)
+    # Cholesky decomposition based eigenvalues and eigenvectors
+    vals, w = eigen!(UtiAUi!(A, C.U))
+    vecs = C.U \ w
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
 
 # Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal}
-UtiAUi!(A::StridedMatrix, U) = _UtiAUi!(A, U)
+UtiAUi!(A, U) = _UtiAUi!(A, U)
 UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo))
 UtiAUi!(A::Hermitian, U) = Hermitian(_UtiAUi!(copytri!(parent(A), A.uplo, true), U), sym_uplo(A.uplo))
-
 _UtiAUi!(A, U) = rdiv!(ldiv!(U', A), U)
 
+function eigvals(A::HermOrSym{TA}, B::HermOrSym{TB}; kws...) where {TA,TB}
+    S = promote_type(eigtype(TA), TB)
+    return eigen!(eigencopy_oftype{S}(A), eigencopy_oftype(B, S); kws...)
+end
+
 function eigvals!(A::HermOrSym{T,S}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
     vals = LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))[1]
     isnothing(sortby) || sort!(vals, by=sortby)
@@ -195,3 +206,15 @@ function eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,N
     return vals
 end
 eigvecs(A::HermOrSym) = eigvecs(eigen(A))
+
+function eigvals(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=nothing)
+    if ishermitian(A)
+        eigvals!(eigencopy_oftype(Hermitian(A), eigtype(eltype(A))), C; sortby)
+    else
+        eigvals!(copy_similar(A, eigtype(eltype(A))), C; sortby)
+    end
+end
+function eigvals!(A::AbstractMatrix{T}, C::Cholesky{T, <:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    # Cholesky decomposition based eigenvalues
+    return eigvals!(UtiAUi!(A, C.U); sortby)
+end

stdlib/LinearAlgebra/test/symmetriceigen.jl

@@ -0,0 +1,72 @@
+# This file is a part of Julia. License is MIT: https://julialang.org/license
+
+module TestSymmetricEigen
+
+using Test, LinearAlgebra
+
+@testset "chol-eigen-eigvals" begin
+    ## Cholesky decomposition based
+
+    # eigenvalue sorting
+    sf = x->(real(x),imag(x))
+
+    ## Real valued
+    A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
+    H = (A+A')/2
+    B = Float64[2 1 4 3; 0 3 1 3; 3 1 0 0; 0 1 3 1]
+    BH = (B+B')/2
+    # PD matrix
+    BPD = B*B'
+    # eigen
+    C = cholesky(BPD)
+    e,v = eigen(A, C; sortby=sf)
+    @test A*v ≈ BPD*v*Diagonal(e)
+    # eigvals
+    @test eigvals(A, BPD; sortby=sf) ≈ eigvals(A, C; sortby=sf)
+
+    ## Complex valued
+    A =  [1.0+im 1.0+1.0im 0 0; 1.0+1.0im 2.0+3.0im 1.0+1.0im 0; 0 1.0+2.0im 3.0+4.0im 1.0+5.0im; 0 0 1.0+1.0im 4.0+4.0im]
+    AH = (A+A')/2
+    B =  [2.0+2.0im 1.0+1.0im 4.0+4.0im 3.0+3.0im; 0 3.0+2.0im 1.0+1.0im 3.0+4.0im; 3.0+3.0im 1.0+4.0im 0 0; 0 1.0+2.0im 3.0+1.0im 1.0+1.0im]
+    BH = (B+B')/2
+    # PD matrix
+    BPD = B*B'
+    # eigen
+    C = cholesky(BPD)
+    e,v = eigen(A, C; sortby=sf)
+    @test A*v ≈ BPD*v*Diagonal(e)
+    # eigvals
+    @test eigvals(A, BPD; sortby=sf) ≈ eigvals(A, C; sortby=sf)
+end
+
+@testset "issue #49533" begin
+    ## Real valued
+    A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
+    B = Matrix(Diagonal(Float64[1:4;]))
+    # eigen
+    e0,v0 = eigen(A, B)
+    e1,v1 = eigen(A, Symmetric(B))
+    e2,v2 = eigen(Symmetric(A), B)
+    @test e0 ≈ e1 && v0 ≈ v1
+    @test e0 ≈ e2 && v0 ≈ v2
+    # eigvals
+    @test eigvals(A, B) ≈ eigvals(A, Symmetric(B))
+    @test eigvals(A, B) ≈ eigvals(Symmetric(A), B)
+
+    ## Complex valued
+    A =  [1.0+im 1.0+1.0im 0 0; 1.0+1.0im 2.0+3.0im 1.0+1.0im 0; 0 1.0+2.0im 3.0+4.0im 1.0+5.0im; 0 0 1.0+1.0im 4.0+4.0im]
+    AH = (A+A')/2
+    B =  [2.0+2.0im 1.0+1.0im 4.0+4.0im 3.0+3.0im; 0 3.0+2.0im 1.0+1.0im 3.0+4.0im; 3.0+3.0im 1.0+4.0im 0 0; 0 1.0+2.0im 3.0+1.0im 1.0+1.0im]
+    BH = (B+B')/2
+    # eigen
+    sf = x->(real(x),imag(x))
+    e1,v1 = eigen(A, Hermitian(BH))
+    e2,v2 = eigen(Hermitian(AH), B)
+    @test A*v1 ≈ Hermitian(BH)*v1*Diagonal(e1)
+    @test Hermitian(AH)*v2 ≈ B*v2*Diagonal(e2)
+    # eigvals
+    @test eigvals(A, BH; sortby=sf) ≈ eigvals(A, Hermitian(BH); sortby=sf)
+    @test eigvals(AH, B; sortby=sf) ≈ eigvals(Hermitian(AH), B; sortby=sf)
+end
+
+end # module TestSymmetricEigen