diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 8caf001d9c2887..2720dfd06ef972 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -198,16 +198,16 @@ Legend: ### Matrix factorizations -| Matrix type | LAPACK | [`eig`](@ref) | [`eigvals`](@ref) | [`eigvecs`](@ref) | [`svd`](@ref) | [`svdvals`](@ref) | -|:------------------------- |:------ |:------------- |:----------------- |:----------------- |:------------- |:----------------- | -| [`Symmetric`](@ref) | SY | | ARI | | | | -| [`Hermitian`](@ref) | HE | | ARI | | | | -| [`UpperTriangular`](@ref) | TR | A | A | A | | | -| [`LowerTriangular`](@ref) | TR | A | A | A | | | -| [`SymTridiagonal`](@ref) | ST | A | ARI | AV | | | -| [`Tridiagonal`](@ref) | GT | | | | | | -| [`Bidiagonal`](@ref) | BD | | | | A | A | -| [`Diagonal`](@ref) | DI | | A | | | | +| Matrix type | LAPACK | [`eigfact`](@ref) | [`eigvals`](@ref) | [`eigvecs`](@ref) | [`svd`](@ref) | [`svdvals`](@ref) | +|:------------------------- |:------ |:----------------- |:----------------- |:----------------- |:------------- |:----------------- | +| [`Symmetric`](@ref) | SY | | ARI | | | | +| [`Hermitian`](@ref) | HE | | ARI | | | | +| [`UpperTriangular`](@ref) | TR | A | A | A | | | +| [`LowerTriangular`](@ref) | TR | A | A | A | | | +| [`SymTridiagonal`](@ref) | ST | A | ARI | AV | | | +| [`Tridiagonal`](@ref) | GT | | | | | | +| [`Bidiagonal`](@ref) | BD | | | | A | A | +| [`Diagonal`](@ref) | DI | | A | | | | Legend: @@ -334,7 +334,6 @@ LinearAlgebra.lqfact LinearAlgebra.lq LinearAlgebra.bkfact LinearAlgebra.bkfact! -LinearAlgebra.eig LinearAlgebra.eigvals LinearAlgebra.eigvals! LinearAlgebra.eigmax diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index f01f87176dd178..cbb223243e101a 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -80,7 +80,6 @@ export diagind, diagm, dot, - eig, eigfact, eigfact!, eigmax, diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index a455acbe8a625d..103a775054170e 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -1282,3 +1282,43 @@ function lu(x::Number) "`LU` object (`lup = lufact(x)`)."), :lu) return first.((lufact(x)...,)) end + +# deprecate eig(...) in favor of eigfact and factorization destructuring +export eig +function eig(A::Union{Number, StridedMatrix}; permute::Bool=true, scale::Bool=true) + depwarn(string("`eig(A[, permute, scale])` has been deprecated in favor of ", + "`eigfact(A[, permute, scale])`. Whereas `eig(A[, permute, scale])` ", + "returns a tuple of arrays, `eigfact(A[, permute, scale])` returns ", + "an `Eigen` object. So for a direct replacement, use ", + "`(eigfact(A[, permute, scale])...,)`. But going forward, consider ", + "using the direct result of `eigfact(A[, permute, scale])` instead, ", + "either destructured into its components ", + "(`vals, vecs = eigfact(A[, permute, scale])`) ", + "or as an `Eigen` object (`eigf = eigfact(A[, permute, scale])`)."), :eig) + return (eigfact(A, permute=permute, scale=scale)...,) +end +function eig(A::AbstractMatrix, args...) + depwarn(string("`eig(A, args...)` has been deprecated in favor of ", + "`eigfact(A, args...)`. Whereas `eig(A, args....)` ", + "returns a tuple of arrays, `eigfact(A, args...)` returns ", + "an `Eigen` object. So for a direct replacement, use ", + "`(eigfact(A, args...)...,)`. But going forward, consider ", + "using the direct result of `eigfact(A, args...)` instead, ", + "either destructured into its components ", + "(`vals, vecs = eigfact(A, args...)`) ", + "or as an `Eigen` object (`eigf = eigfact(A, args...)`)."), :eig) + return (eigfact(A, args...)...,) +end +eig(A::AbstractMatrix, B::AbstractMatrix) = _geneig(A, B) +eig(A::Number, B::Number) = _geneig(A, B) +function _geneig(A, B) + depwarn(string("`eig(A::AbstractMatrix, B::AbstractMatrix)` and ", + "`eig(A::Number, B::Number)` have been deprecated in favor of ", + "`eigfact(A, B)`. Whereas the former each return a tuple of arrays, ", + "the latter returns a `GeneralizedEigen` object. So for a direct ", + "replacement, use `(eigfact(A, B)...,)`. But going forward, consider ", + "using the direct result of `eigfact(A, B)` instead, either ", + "destructured into its components (`vals, vecs = eigfact(A, B)`), ", + "or as a `GeneralizedEigen` object (`eigf = eigfact(A, B)`)."), :eig) + return (eigfact(A, B)...,) +end diff --git a/stdlib/LinearAlgebra/src/eigen.jl b/stdlib/LinearAlgebra/src/eigen.jl index 56890060f13f2a..560a9884f77533 100644 --- a/stdlib/LinearAlgebra/src/eigen.jl +++ b/stdlib/LinearAlgebra/src/eigen.jl @@ -20,6 +20,13 @@ end GeneralizedEigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V} = GeneralizedEigen{T,V,typeof(vectors),typeof(values)}(values, vectors) +# iteration for destructuring into factors +Base.start(::Union{Eigen,GeneralizedEigen}) = Val(:values) +Base.next(F::Union{Eigen,GeneralizedEigen}, ::Val{:values}) = (F.values, Val(:vectors)) +Base.next(F::Union{Eigen,GeneralizedEigen}, ::Val{:vectors}) = (F.vectors, Val(:done)) +Base.done(F::Union{Eigen,GeneralizedEigen}, ::Val{:done}) = true +Base.done(F::Union{Eigen,GeneralizedEigen}, ::Any) = false + isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values) """ @@ -94,6 +101,20 @@ julia> F.values 18.0 julia> F.vectors +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.0 + +julia> vals, vecs = F; # destructuring via iteration + +julia> vals +3-element Array{Float64,1}: + 1.0 + 3.0 + 18.0 + +julia> vecs 3×3 Array{Float64,2}: 1.0 0.0 0.0 0.0 1.0 0.0 @@ -107,38 +128,6 @@ function eigfact(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) wher end eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1)) -function eig(A::Union{Number, StridedMatrix}; permute::Bool=true, scale::Bool=true) - F = eigfact(A, permute=permute, scale=scale) - F.values, F.vectors -end - -""" - eig(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> D, V - eig(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> D, V - eig(A, permute::Bool=true, scale::Bool=true) -> D, V - -Computes eigenvalues (`D`) and eigenvectors (`V`) of `A`. -See [`eigfact`](@ref) for details on the -`irange`, `vl`, and `vu` arguments -(for [`SymTridiagonal`](@ref), [`Hermitian`](@ref), and -[`Symmetric`](@ref) matrices) -and the `permute` and `scale` keyword arguments. -The eigenvectors are returned columnwise. - -# Examples -```jldoctest -julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]) -([1.0, 3.0, 18.0], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]) -``` - -`eig` is a wrapper around [`eigfact`](@ref), extracting all parts of the -factorization to a tuple; where possible, using [`eigfact`](@ref) is recommended. -""" -function eig(A::AbstractMatrix, args...) - F = eigfact(A, args...) - F.values, F.vectors -end - """ eigvecs(A; permute::Bool=true, scale::Bool=true) -> Matrix @@ -378,6 +367,18 @@ julia> F.values 0.0 - 1.0im julia> F.vectors +2×2 Array{Complex{Float64},2}: + 0.0-1.0im 0.0+1.0im + -1.0-0.0im -1.0+0.0im + +julia> vals, vecs = F; # destructuring via iterationw + +julia> vals +2-element Array{Complex{Float64},1}: + 0.0 + 1.0im + 0.0 - 1.0im + +julia> vecs 2×2 Array{Complex{Float64},2}: 0.0-1.0im 0.0+1.0im -1.0-0.0im -1.0+0.0im @@ -390,38 +391,6 @@ end eigfact(A::Number, B::Number) = eigfact(fill(A,1,1), fill(B,1,1)) -""" - eig(A, B) -> D, V - -Computes generalized eigenvalues (`D`) and vectors (`V`) of `A` with respect to `B`. - -`eig` is a wrapper around [`eigfact`](@ref), extracting all parts of the -factorization to a tuple; where possible, using [`eigfact`](@ref) is recommended. - -# Examples -```jldoctest -julia> A = [1 0; 0 -1] -2×2 Array{Int64,2}: - 1 0 - 0 -1 - -julia> B = [0 1; 1 0] -2×2 Array{Int64,2}: - 0 1 - 1 0 - -julia> eig(A, B) -(Complex{Float64}[0.0+1.0im, 0.0-1.0im], Complex{Float64}[0.0-1.0im 0.0+1.0im; -1.0-0.0im -1.0+0.0im]) -``` -""" -function eig(A::AbstractMatrix, B::AbstractMatrix) - F = eigfact(A,B) - F.values, F.vectors -end -function eig(A::Number, B::Number) - F = eigfact(A,B) - F.values, F.vectors -end """ eigvals!(A, B) -> values diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index c6ef222cacc7b7..2bdabd65a85b6b 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -2156,7 +2156,7 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat R[i,i+1] = i / sqrt((2 * i)^2 - 1) R[i+1,i] = R[i,i+1] end - x,V = eig(R) + x,V = eigfact(R) w = Vector{Float64}(undef, m) for i = 1:m x[i] = (x[i] + 1) / 2 diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index a096ffd43c27f1..6d83933ba248b9 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -231,8 +231,8 @@ srand(1) @testset "Eigensystems" begin if relty <: AbstractFloat - d1, v1 = eig(T) - d2, v2 = eig(map(elty<:Complex ? ComplexF64 : Float64,Tfull)) + d1, v1 = eigfact(T) + d2, v2 = eigfact(map(elty<:Complex ? ComplexF64 : Float64,Tfull)) @test (uplo == :U ? d1 : reverse(d1)) ≈ d2 if elty <: Real test_approx_eq_modphase(v1, uplo == :U ? v2 : v2[:,n:-1:1]) diff --git a/stdlib/LinearAlgebra/test/eigen.jl b/stdlib/LinearAlgebra/test/eigen.jl index eec694ed7f475e..15f11eda728331 100644 --- a/stdlib/LinearAlgebra/test/eigen.jl +++ b/stdlib/LinearAlgebra/test/eigen.jl @@ -28,12 +28,12 @@ aimg = randn(n,n)/2 α = rand(eltya) β = rand(eltya) - eab = eig(α,β) + eab = (eigfact(α,β)...,) @test eab[1] == eigvals(fill(α,1,1),fill(β,1,1)) @test eab[2] == eigvecs(fill(α,1,1),fill(β,1,1)) @testset "non-symmetric eigen decomposition" begin - d, v = eig(a) + d, v = eigfact(a) for i in 1:size(a,2) @test a*v[:,i] ≈ d[i]*v[:,i] end @@ -70,7 +70,7 @@ aimg = randn(n,n)/2 @test eigvecs(f) === f.vectors @test_throws ErrorException f.Z - d,v = eig(asym_sg, a_sg'a_sg) + d,v = eigfact(asym_sg, a_sg'a_sg) @test d == f.values @test v == f.vectors end @@ -89,7 +89,7 @@ aimg = randn(n,n)/2 @test eigvecs(a1_nsg, a2_nsg) == f.vectors @test_throws ErrorException f.Z - d,v = eig(a1_nsg, a2_nsg) + d,v = eigfact(a1_nsg, a2_nsg) @test d == f.values @test v == f.vectors end @@ -98,11 +98,11 @@ end @testset "eigenvalue computations with NaNs" begin for eltya in (NaN16, NaN32, NaN) - @test_throws(ArgumentError, eig(fill(eltya, 1, 1))) - @test_throws(ArgumentError, eig(fill(eltya, 2, 2))) + @test_throws(ArgumentError, eigfact(fill(eltya, 1, 1))) + @test_throws(ArgumentError, eigfact(fill(eltya, 2, 2))) test_matrix = rand(typeof(eltya),3,3) test_matrix[2,2] = eltya - @test_throws(ArgumentError, eig(test_matrix)) + @test_throws(ArgumentError, eigfact(test_matrix)) end end diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 0a1c42a103e817..1cf83b466e531a 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -51,7 +51,7 @@ dimg = randn(n)/2 -eps(real(one(eltya)))/4 eps(real(one(eltya)))/2 -1.0 0; -0.5 -0.5 0.1 1.0]) F = eigfact(A, permute=false, scale=false) - eig(A, permute=false, scale=false) + (eigfact(A, permute=false, scale=false)...,) @test F.vectors*Diagonal(F.values)/F.vectors ≈ A F = eigfact(A) # @test norm(F.vectors*Diagonal(F.values)/F.vectors - A) > 0.01 diff --git a/stdlib/LinearAlgebra/test/schur.jl b/stdlib/LinearAlgebra/test/schur.jl index 52ba3957fb5052..2dac81d3b201c2 100644 --- a/stdlib/LinearAlgebra/test/schur.jl +++ b/stdlib/LinearAlgebra/test/schur.jl @@ -26,7 +26,7 @@ aimg = randn(n,n)/2 view(apd, 1:n, 1:n))) ε = εa = eps(abs(float(one(eltya)))) - d,v = eig(a) + d,v = eigfact(a) f = schurfact(a) @test f.vectors*f.Schur*f.vectors' ≈ a @test sort(real(f.values)) ≈ sort(real(d)) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 2a6a279a397a09..dfadea2a33bc7c 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -218,14 +218,14 @@ end @testset "symmetric eigendecomposition" begin if eltya <: Real # the eigenvalues are only real and ordered for Hermitian matrices - d, v = eig(asym) + d, v = eigfact(asym) @test asym*v[:,1] ≈ d[1]*v[:,1] @test v*Diagonal(d)*transpose(v) ≈ asym @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1])) @test abs.(eigfact(Symmetric(asym), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2) - eig(Symmetric(asym), 1:2) # same result, but checks that method works + (eigfact(Symmetric(asym), 1:2)...,) # same result, but checks that method works @test abs.(eigfact(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2) - eig(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2) # same result, but checks that method works + (eigfact(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2)...,) # same result, but checks that method works @test eigvals(Symmetric(asym), 1:2) ≈ d[1:2] @test eigvals(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2] # eigfact doesn't support Symmetric{Complex} @@ -233,14 +233,14 @@ end @test eigvecs(Symmetric(asym)) ≈ eigvecs(asym) end - d, v = eig(aherm) + d, v = eigfact(aherm) @test aherm*v[:,1] ≈ d[1]*v[:,1] @test v*Diagonal(d)*v' ≈ aherm @test isequal(eigvals(aherm[1]), eigvals(aherm[1:1,1:1])) @test abs.(eigfact(Hermitian(aherm), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2) - eig(Hermitian(aherm), 1:2) # same result, but checks that method works + (eigfact(Hermitian(aherm), 1:2)...,) # same result, but checks that method works @test abs.(eigfact(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2) - eig(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2) # same result, but checks that method works + (eigfact(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2)...,) # same result, but checks that method works @test eigvals(Hermitian(aherm), 1:2) ≈ d[1:2] @test eigvals(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2] @test Matrix(eigfact(aherm)) ≈ aherm @@ -365,7 +365,7 @@ end end @testset "Issues #8057 and #8058. f=$f, A=$A" for f in - (eigfact, eigvals, eig), + (eigfact, eigvals, eigfact), A in (Symmetric([0 1; 1 0]), Hermitian([0 im; -im 0])) @test_throws ArgumentError f(A, 3, 2) @test_throws ArgumentError f(A, 1:4) diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 50b94cc6f77647..49c8a13f1f726e 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -249,7 +249,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo # eigenproblems if !(elty1 in (BigFloat, Complex{BigFloat})) # Not handled yet - vals, vecs = eig(A1) + vals, vecs = eigfact(A1) if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues. @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(eps(float(real(one(vals[1])))))*(norm(A1,Inf)*n)^2 end diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index 679ed2dfd20a08..639e5ced030f24 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -243,7 +243,7 @@ end w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) evecs = LAPACK.stein!(b, a, w) - (e, v) = eig(SymTridiagonal(b, a)) + (e, v) = eigfact(SymTridiagonal(b, a)) @test e ≈ w test_approx_eq_vecs(v, evecs) end @@ -266,10 +266,10 @@ end end @testset "eigenvalues/eigenvectors of symmetric tridiagonal" begin if elty === Float32 || elty === Float64 - DT, VT = @inferred eig(A) - @inferred eig(A, 2:4) - @inferred eig(A, 1.0, 2.0) - D, Vecs = eig(fA) + DT, VT = @inferred eigfact(A) + @inferred eigfact(A, 2:4) + @inferred eigfact(A, 1.0, 2.0) + D, Vecs = eigfact(fA) @test DT ≈ D @test abs.(VT'Vecs) ≈ Matrix(elty(1)I, n, n) test_approx_eq_modphase(eigvecs(A), eigvecs(fA)) diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index 23a33b39985214..0f82d111a618ee 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -1009,7 +1009,6 @@ function factorize(A::LinearAlgebra.RealHermSymComplexHerm{Float64,<:SparseMatri end chol(A::SparseMatrixCSC) = error("Use cholfact() instead of chol() for sparse matrices.") -eig(A::SparseMatrixCSC) = error("Use IterativeEigensolvers.eigs() instead of eig() for sparse matrices.") function Base.cov(X::SparseMatrixCSC; dims::Int=1, corrected::Bool=true) vardim = dims diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index d2a48d8a32b8da..f9a86449fbc0de 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -1780,7 +1780,6 @@ end C, b = A[:, 1:4], fill(1., size(A, 1)) @test !Base.USE_GPL_LIBS || factorize(C)\b ≈ Array(C)\b @test_throws ErrorException chol(A) - @test_throws ErrorException eig(A) @test_throws ErrorException inv(A) end