Rework symmetric generalized eigen/eigvals

NEWS.md

@@ -92,6 +92,12 @@ 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))` and `eigvals/eigen(A, bunchkaufman(B))` now compute the generalized
+  eigenvalues (`eigen`: and eigenvectors) of `A` and `B` via Cholesky decomposition for positive definite `B`
+  and Bunchkaufmann (LDLT) decomposition for symmetric/Hermitian `B`, respectively. Note: The second argument is
+  the output of `cholesky` or `bunchkaufman`.
 
 #### Printf
 * Format specifiers now support dynamic width and precision, e.g. `%*s` and `%*.*g` ([#40105]).

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -164,26 +164,44 @@ 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)
+function eigen(A::RealHermSymComplexHerm{T,<:StridedMatrix}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    return eigen!(Matrix{T}(A), eigencopy_oftype(B, T); 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::StridedMatrix{T}, B::Union{RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    return eigen!(eigencopy_oftype(A, T), Matrix{T}(B); sortby)
 end
-function _choleigen!(A, B, sortby)
-    U = cholesky(B).U
-    vals, w = eigen!(UtiAUi!(A, U))
-    vecs = U \ w
+
+function eigen(A::AbstractMatrix{T}, C::Cholesky{T, <:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    # Cholesky decomposition based eigenvalues and eigenvectors
+    vals, w = eigen!(UtiAUi!(copy(A), C.U))
+    vecs = C.U \ w
+    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
+end
+
+function eigen(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    # Bunchkaufman decomposition based eigenvalues and eigenvectors
+    if B.uplo == 'U'
+        vals, w = eigen!(ldiv(lu!(copy(B.D)), UiAUti!(A[B.p,B.p], B.U)); sortby)
+        vecs = (B.U' \ w)[invperm(B.p),:]
+    else # B.uplo == 'L'
+        vals, w = eigen!(ldiv(lu!(copy(B.D)), UiAUti!(A[B.p,B.p], B.L)); sortby)
+        vecs = (B.L' \ w)[invperm(B.p),:]
+    end
     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)
 
+# Perform U \ A / U' in-place, where U::Union{UpperTriangular,Diagonal}
+UiAUti!(A, U) = _UiAUti!(A, U)
+UiAUti!(A::Symmetric, U) = Symmetric(_UiAUti!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo))
+UiAUti!(A::Hermitian, U) = Hermitian(_UiAUti!(copytri!(parent(A), A.uplo, true), U), sym_uplo(A.uplo))
+_UiAUti!(A, U) = rdiv!(ldiv!(U, A), U')
+
 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 +213,26 @@ function eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,N
     return vals
 end
 eigvecs(A::HermOrSym) = eigvecs(eigen(A))
+
+# Note: No specilized LAPACK routines for Matrix+Symmetric and Symmetric+Matrix exist. Hence, the calls are forwarded to conventional Matrix eigvals!
+function eigvals(A::StridedMatrix{T}, B::Union{RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    return eigvals!(eigencopy_oftype(A, T), Matrix{T}(B); sortby)
+end
+
+function eigvals(A::RealHermSymComplexHerm{T,<:StridedMatrix}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    return eigvals!(Matrix{T}(A), eigencopy_oftype(B, T); sortby)
+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!(eigencopy_oftype(A, T), C.U); sortby)
+end
+
+# Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues
+function eigvals(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    if B.uplo == 'U'
+        return eigvals!(ldiv(lu!(copy(B.D)), UiAUti!(A[B.p,B.p], B.U)); sortby)
+    else # B.uplo == 'L'
+        return eigvals!(ldiv(lu!(copy(B.D)), UiAUti!(A[B.p,B.p], B.L)); sortby)
+    end
+end

stdlib/LinearAlgebra/test/symmetriceigen.jl

@@ -0,0 +1,120 @@
+# This file is a part of Julia. License is MIT: https://julialang.org/license
+
+module TestSymmetricEigen
+
+using Test, LinearAlgebra, Random
+
+Random.seed!(555)
+
+@testset "bk-eigen-eigals" begin
+    # Bunchkaufman decomposition based
+
+    # Real-valued random matrix
+    N = 10
+    A = randn(N,N)
+    B = randn(N,N)
+    BH = (B+B')/2
+    sf = x->(real(x),imag(x))
+    # eigen
+    e0 = eigvals(A,BH; sortby=sf)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    # eigvals
+    e0 = eigvals(A,BH; sortby=sf)
+    el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ el
+    @test e0 ≈ eu
+
+    # Complex-valued random matrix
+    N = 10
+    A = complex.(randn(N,N),randn(N,N))
+    B = complex.(randn(N,N),randn(N,N))
+    BH = (B+B')/2
+    sf = x->(real(x),imag(x))
+    # eigen
+    e0 = eigvals(A,BH; sortby=sf)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    # eigvals
+    e0 = eigvals(A,BH; sortby=sf)
+    el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ el
+    @test e0 ≈ eu
+end
+
+@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