From e89c21b3789ec5217e2b84a3ea4729b7b5c3d05b Mon Sep 17 00:00:00 2001 From: Jae-Mo Lihm Date: Wed, 22 May 2024 14:56:24 +0200 Subject: [PATCH] Add interface to use symmetric eigendecomposition with different lapack method (#49355) --- NEWS.md | 3 + stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 + stdlib/LinearAlgebra/src/lapack.jl | 5 -- stdlib/LinearAlgebra/src/svd.jl | 5 +- stdlib/LinearAlgebra/src/symmetriceigen.jl | 81 +++++++++++++++++++--- stdlib/LinearAlgebra/test/symmetric.jl | 16 +++++ 6 files changed, 94 insertions(+), 17 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6a3ad9246d1a1..b5c304714d878 100644 --- a/NEWS.md +++ b/NEWS.md @@ -88,6 +88,9 @@ Standard library changes #### LinearAlgebra * `rank` can now take a `QRPivoted` matrix to allow rank estimation via QR factorization ([#54283]). +* Added keyword argument `alg` to `eigen`, `eigen!`, `eigvals` and `eigvals!` for self-adjoint + matrix types (i.e., the type union `RealHermSymComplexHerm`) that allows one to switch + between different eigendecomposition algorithms ([#49355]). #### Logging diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 7c4a40c11ef24..34605700f73cd 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -190,6 +190,7 @@ end abstract type Algorithm end struct DivideAndConquer <: Algorithm end struct QRIteration <: Algorithm end +struct RobustRepresentations <: Algorithm end # Pivoting strategies for matrix factorization algorithms. abstract type PivotingStrategy end diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index fd485b096395e..b7c624dbb494d 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -5762,11 +5762,6 @@ syevr!(jobz::AbstractChar, range::AbstractChar, uplo::AbstractChar, A::AbstractM Finds the eigenvalues (`jobz = N`) or eigenvalues and eigenvectors (`jobz = V`) of a symmetric matrix `A`. If `uplo = U`, the upper triangle of `A` is used. If `uplo = L`, the lower triangle of `A` is used. - -Use the divide-and-conquer method, instead of the QR iteration used by -`syev!` or multiple relatively robust representations used by `syevr!`. -See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for -a comparison of the accuracy and performatce of different methods. """ syevd!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index dc8d717932be3..7a88c4a6e14c4 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -149,8 +149,9 @@ and `V` is ``N \\times N``, while in the thin factorization `U` is ``M \\times K`` and `V` is ``N \\times K``, where ``K = \\min(M,N)`` is the number of singular values. -If `alg = DivideAndConquer()` a divide-and-conquer algorithm is used to calculate the SVD. -Another (typically slower but more accurate) option is `alg = QRIteration()`. +`alg` specifies which algorithm and LAPACK method to use for SVD: +- `alg = DivideAndConquer()` (default): Calls `LAPACK.gesdd!`. +- `alg = QRIteration()`: Calls `LAPACK.gesvd!` (typically slower but more accurate) . !!! compat "Julia 1.3" The `alg` keyword argument requires Julia 1.3 or later. diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl index 8126df8ff9077..666b9a9bc81df 100644 --- a/stdlib/LinearAlgebra/src/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl @@ -5,13 +5,19 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) eigencopy_oftype(A::Symmetric, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) -# Eigensolvers for symmetric and Hermitian matrices -eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) = - Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...) +default_eigen_alg(A) = DivideAndConquer() -function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing) - S = eigtype(eltype(A)) - eigen!(eigencopy_oftype(A, S), sortby=sortby) +# Eigensolvers for symmetric and Hermitian matrices +function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) + if alg === DivideAndConquer() + Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...) + elseif alg === QRIteration() + Eigen(sorteig!(LAPACK.syev!('V', A.uplo, A.data)..., sortby)...) + elseif alg === RobustRepresentations() + Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...) + else + throw(ArgumentError("Unsupported value for `alg` keyword.")) + end end function eigen(A::RealHermSymComplexHerm{Float16}; sortby::Union{Function,Nothing}=nothing) S = eigtype(eltype(A)) @@ -21,6 +27,36 @@ function eigen(A::RealHermSymComplexHerm{Float16}; sortby::Union{Function,Nothin return Eigen(values, vectors) end +""" + eigen(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A)) -> Eigen + +Compute the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F` +which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the +matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.) + +Iterating the decomposition produces the components `F.values` and `F.vectors`. + +`alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition: +- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`. +- `alg = QRIteration()`: Calls `LAPACK.syev!`. +- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`. + +See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for +a comparison of the accuracy and performance of different algorithms. + +The default `alg` used may change in the future. + +!!! compat "Julia 1.12" + The `alg` keyword argument requires Julia 1.12 or later. + +The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref). +""" +function eigen(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) + S = eigtype(eltype(A)) + eigen!(eigencopy_oftype(A, S), alg; sortby) +end + + eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...) @@ -71,17 +107,42 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real) eigen!(eigencopy_oftype(A, S), vl, vh) end -function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) - vals = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1] + +function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) + vals::Vector{real(eltype(A))} = if alg === DivideAndConquer() + LAPACK.syevd!('N', A.uplo, A.data) + elseif alg === QRIteration() + LAPACK.syev!('N', A.uplo, A.data) + elseif alg === RobustRepresentations() + LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1] + else + throw(ArgumentError("Unsupported value for `alg` keyword.")) + end !isnothing(sortby) && sort!(vals, by=sortby) return vals end -function eigvals(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing) +""" + eigvals(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A))) -> values + +Return the eigenvalues of `A`. + +`alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition: +- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`. +- `alg = QRIteration()`: Calls `LAPACK.syev!`. +- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`. + +See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for +a comparison of the accuracy and performance of different methods. + +The default `alg` used may change in the future. +""" +function eigvals(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) S = eigtype(eltype(A)) - eigvals!(eigencopy_oftype(A, S), sortby=sortby) + eigvals!(eigencopy_oftype(A, S), alg; sortby) end + """ eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index efd8aa7b8d300..8f1e2b7420b22 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -720,6 +720,22 @@ end end end +@testset "eigendecomposition Algorithms" begin + using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations + for T in (Float64, ComplexF64, Float32, ComplexF32) + n = 4 + A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n)) + d, v = eigen(A) + for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations()) + @test (@inferred eigvals(A, alg)) ≈ d + d2, v2 = @inferred eigen(A, alg) + @test d2 ≈ d + @test A * v2 ≈ v2 * Diagonal(d2) + end + end +end + +const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) using .Main.ImmutableArrays