From 3c94a27f02a97296179fdc438f38896598f3f2be Mon Sep 17 00:00:00 2001 From: Sacha Verweij Date: Sun, 13 May 2018 15:01:10 -0700 Subject: [PATCH] Deprecate schur(...) in favor of schurfact(...) and factorization destructuring. --- NEWS.md | 15 +++ stdlib/LinearAlgebra/docs/src/index.md | 1 - stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 - stdlib/LinearAlgebra/src/dense.jl | 12 +- stdlib/LinearAlgebra/src/deprecated.jl | 23 ++++ stdlib/LinearAlgebra/src/schur.jl | 155 +++++++++++++++------- stdlib/LinearAlgebra/test/lapack.jl | 4 +- stdlib/LinearAlgebra/test/schur.jl | 14 +- 8 files changed, 157 insertions(+), 68 deletions(-) diff --git a/NEWS.md b/NEWS.md index f87772497c1e5..11b4ac06c4c16 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1025,6 +1025,21 @@ Deprecated or removed (`vals, vecs = eigfact(A, B)`), or as a `GeneralizedEigen` object (`eigf = eigfact(A, B)`) ([#26997]). + * `schur(A::AbstractMatrix)` has been deprecated in favor of `schurfact(A)`. + Whereas the former returns a tuple of arrays, the latter returns a `Schur` object. + So for a direct replacement, use `(schurfact(A)...,)`. But going forward, consider + using the direct result of `schurfact(A)` instead, either destructured into + its components (`T, Z, λ = schurfact(A)`) or as a `Schur` object + (`schurf = schurfact(A)`) ([#26997]). + + * `schur(A::StridedMatrix, B::StridedMatrix)` has been deprecated in favor of + `schurfact(A, B)`. Whereas the former returns a tuple of arrays, the latter + returns a `GeneralizedSchur` object. So for a direct replacement, use + `(schurfact(A, B)...,)`. But going forward, consider using the direct result + of `schurfact(A, B)` instead, either destructured into its components + (`S, T, Q, Z, α, β = schurfact(A, B)`) or as a `GeneralizedSchur` object + (`schurf = schurfact(A, B)`) ([#26997]). + * The timing functions `tic`, `toc`, and `toq` are deprecated in favor of `@time` and `@elapsed` ([#17046]). diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 2720dfd06ef97..7ab9c2fb63fc4 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -345,7 +345,6 @@ LinearAlgebra.hessfact LinearAlgebra.hessfact! LinearAlgebra.schurfact LinearAlgebra.schurfact! -LinearAlgebra.schur LinearAlgebra.ordschur LinearAlgebra.ordschur! LinearAlgebra.svdfact diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index cbb223243e101..156d20f4e9441 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -131,7 +131,6 @@ export lqfact, rank, rdiv!, - schur, schurfact!, schurfact, svd, diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 9e1bfa4ab0454..cca39e9bdd758 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -418,7 +418,7 @@ function schurpow(A::AbstractMatrix, p) retmat = retmat * powm!(UpperTriangular(float.(A)), real(p - floor(p))) end else - S,Q,d = schur(complex(A)) + S,Q,d = schurfact(complex(A)) # Integer part R = S ^ floor(p) # Real part @@ -605,7 +605,7 @@ matrix function is returned whenever possible. If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used, if `A` is triangular an improved version of the inverse scaling and squaring method is employed (see [^AH12] and [^AHR13]). For general matrices, the complex Schur form -([`schur`](@ref)) is computed and the triangular algorithm is used on the +([`schurfact`](@ref)) is computed and the triangular algorithm is used on the triangular factor. [^AH12]: Awad H. Al-Mohy and Nicholas J. Higham, "Improved inverse scaling and squaring algorithms for the matrix logarithm", SIAM Journal on Scientific Computing, 34(4), 2012, C153-C169. [doi:10.1137/110852553](https://doi.org/10.1137/110852553) @@ -662,7 +662,7 @@ that is the unique matrix ``X`` with eigenvalues having positive real part such If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to compute the square root. Otherwise, the square root is determined by means of the -Björck-Hammarling method [^BH83], which computes the complex Schur form ([`schur`](@ref)) +Björck-Hammarling method [^BH83], which computes the complex Schur form ([`schurfact`](@ref)) and then the complex square root of the triangular factor. [^BH83]: @@ -1409,8 +1409,8 @@ julia> A*X + X*B + C ``` """ function sylvester(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) where T<:BlasFloat - RA, QA = schur(A) - RB, QB = schur(B) + RA, QA = schurfact(A) + RB, QB = schurfact(B) D = -(adjoint(QA) * (C*QB)) Y, scale = LAPACK.trsyl!('N','N', RA, RB, D) @@ -1453,7 +1453,7 @@ julia> A*X + X*A' + B ``` """ function lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:BlasFloat} - R, Q = schur(A) + R, Q = schurfact(A) D = -(adjoint(Q) * (C*Q)) Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D) diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 103a775054170..7ed0abcf42aa2 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -1322,3 +1322,26 @@ function _geneig(A, B) "or as a `GeneralizedEigen` object (`eigf = eigfact(A, B)`)."), :eig) return (eigfact(A, B)...,) end + +# deprecate schur(...) in favor of schurfact(...) and factorization destructuring +export schur +function schur(A::Union{StridedMatrix,Symmetric,Hermitian,UpperTriangular,LowerTriangular,Tridiagonal}) + depwarn(string("`schur(A::AbstractMatrix)` has been deprecated in favor of ", + "`schurfact(A)`. Whereas the former returns a tuple of arrays, the ", + "latter returns a `Schur` object. So for a direct replacement, ", + "use `(schurfact(A)...,)`. But going forward, consider using ", + "the direct result of `schurfact(A)` instead, either destructured ", + "into its components (`T, Z, λ = schurfact(A)`) or as a `Schur` ", + "object (`schurf = schurfact(A)`)."), :schur) + return (schurfact(A)...,) +end +function schur(A::StridedMatrix, B::StridedMatrix) + depwarn(string("`schur(A::StridedMatrix, B::StridedMatrix)` has been ", + "deprecated in favor of `schurfact(A, B)`. Whereas the former returns ", + "a tuple of arrays, the latter returns a `GeneralizedSchur` object. ", + "So for a direct replacement, use `(schurfact(A, B)...,)`. But going ", + "forward, consider using the direct result of `schurfact(A, B)` instead, ", + "either destructured into its components (`S, T, Q, Z, α, β = schurfact(A, B)`) ", + " or as a `GeneralizedSchur` object (`schurf = schurfact(A, B)`)."), :schur) + return (schurfact(A, B)...,) +end diff --git a/stdlib/LinearAlgebra/src/schur.jl b/stdlib/LinearAlgebra/src/schur.jl index ec4595351bfe2..67c1c93d23e52 100644 --- a/stdlib/LinearAlgebra/src/schur.jl +++ b/stdlib/LinearAlgebra/src/schur.jl @@ -9,6 +9,14 @@ struct Schur{Ty,S<:AbstractMatrix} <: Factorization{Ty} end Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) where {Ty} = Schur{Ty, typeof(T)}(T, Z, values) +# iteration for destructuring into factors +Base.start(::Schur) = Val(:T) +Base.next(F::Schur, ::Val{:T}) = (F.T, Val(:Z)) +Base.next(F::Schur, ::Val{:Z}) = (F.Z, Val(:values)) +Base.next(F::Schur, ::Val{:values}) = (F.values, Val(:done)) +Base.done(F::Schur, ::Val{:done}) = true +Base.done(F::Schur, ::Any) = false + """ schurfact!(A::StridedMatrix) -> F::Schur @@ -78,10 +86,34 @@ julia> F.vectors * F.Schur * F.vectors' 2×2 Array{Float64,2}: 5.0 7.0 -2.0 -4.0 + +julia> T, Z, λ = schurfact(A); # destructuring via iteration + +julia> T +2×2 Array{Float64,2}: + 3.0 9.0 + 0.0 -2.0 + +julia> Z +2×2 Array{Float64,2}: + 0.961524 0.274721 + -0.274721 0.961524 + +julia> λ +2-element Array{Float64,1}: + 3.0 + -2.0 ``` """ schurfact(A::StridedMatrix{<:BlasFloat}) = schurfact!(copy(A)) schurfact(A::StridedMatrix{T}) where T = schurfact!(copy_oftype(A, eigtype(T))) +schurfact(A::Tridiagonal) = _schurfact_structured(A) +schurfact(A::Symmetric) = _schurfact_annotated(A) +schurfact(A::Hermitian) = _schurfact_annotated(A) +schurfact(A::UpperTriangular) = _schurfact_annotated(A) +schurfact(A::LowerTriangular) = _schurfact_annotated(A) +_schurfact_annotated(A) = schurfact(copyto!(similar(parent(A)), A)) +_schurfact_structured(A) = schurfact(Matrix(A)) function getproperty(F::Schur, d::Symbol) if d == :Schur @@ -106,47 +138,6 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Schur) show(io, mime, F.values) end -""" - schur(A::StridedMatrix) -> T::Matrix, Z::Matrix, λ::Vector - -Computes the Schur factorization of the matrix `A`. The methods return the (quasi) -triangular Schur factor `T` and the orthogonal/unitary Schur vectors `Z` such that -`A = Z * T * Z'`. The eigenvalues of `A` are returned in the vector `λ`. - -See [`schurfact`](@ref). - -# Examples -```jldoctest -julia> A = [5. 7.; -2. -4.] -2×2 Array{Float64,2}: - 5.0 7.0 - -2.0 -4.0 - -julia> T, Z, lambda = schur(A) -([3.0 9.0; 0.0 -2.0], [0.961524 0.274721; -0.274721 0.961524], [3.0, -2.0]) - -julia> Z * Z' -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 1.0 - -julia> Z * T * Z' -2×2 Array{Float64,2}: - 5.0 7.0 - -2.0 -4.0 -``` -""" -function schur(A::StridedMatrix) - SchurF = schurfact(A) - SchurF.T, SchurF.Z, SchurF.values -end -schur(A::Symmetric) = schur(copyto!(similar(parent(A)), A)) -schur(A::Hermitian) = schur(copyto!(similar(parent(A)), A)) -schur(A::UpperTriangular) = schur(copyto!(similar(parent(A)), A)) -schur(A::LowerTriangular) = schur(copyto!(similar(parent(A)), A)) -schur(A::Tridiagonal) = schur(Matrix(A)) - - """ ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur @@ -209,6 +200,17 @@ function GeneralizedSchur(S::AbstractMatrix{Ty}, T::AbstractMatrix{Ty}, alpha::V GeneralizedSchur{Ty, typeof(S)}(S, T, alpha, beta, Q, Z) end +# iteration for destructuring into factors +Base.start(::GeneralizedSchur) = Val(:S) +Base.next(F::GeneralizedSchur, ::Val{:S}) = (F.S, Val(:T)) +Base.next(F::GeneralizedSchur, ::Val{:T}) = (F.T, Val(:Q)) +Base.next(F::GeneralizedSchur, ::Val{:Q}) = (F.Q, Val(:Z)) +Base.next(F::GeneralizedSchur, ::Val{:Z}) = (F.Z, Val(:α)) +Base.next(F::GeneralizedSchur, ::Val{:α}) = (F.α, Val(:β)) +Base.next(F::GeneralizedSchur, ::Val{:β}) = (F.β, Val(:done)) +Base.done(F::GeneralizedSchur, ::Val{:done}) = true +Base.done(F::GeneralizedSchur, ::Any) = false + """ schurfact!(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur @@ -226,6 +228,67 @@ and `F.T`, the left unitary/orthogonal Schur vectors can be obtained with `F.lef `F.Q` and the right unitary/orthogonal Schur vectors can be obtained with `F.right` or `F.Z` such that `A=F.left*F.S*F.right'` and `B=F.left*F.T*F.right'`. The generalized eigenvalues of `A` and `B` can be obtained with `F.α./F.β`. + +# Examples +```jldoctest +julia> A = [5. 7.; -2. -4.] +2×2 Array{Float64,2}: + 5.0 7.0 + -2.0 -4.0 + +julia> B = [6. 8.; -3. -5.] +2×2 Array{Float64,2}: + 6.0 8.0 + -3.0 -5.0 + +julia> F = schurfact(A, B) +GeneralizedSchur{Float64,Array{Float64,2}} +S factor: +2×2 Array{Float64,2}: + 9.65794 0.323267 + -0.504339 0.604369 +T factor: +2×2 Array{Float64,2}: + 11.5642 0.0 + 0.0 0.518842 +Q factor: +2×2 Array{Float64,2}: + -0.864443 0.50273 + 0.50273 0.864443 +Z factor: +2×2 Array{Float64,2}: + -0.578929 0.815378 + -0.815378 -0.578929 +α: +2-element Array{Complex{Float64},1}: + 2.0000000000000044 + 2.738953629966988e-8im + 2.9999999999999907 - 4.10843044495046e-8im +β: +2-element Array{Float64,1}: + 2.0000000000000058 + 2.9999999999999925 + +julia> F.Q * F.S * F.Z' +2×2 Array{Float64,2}: + 5.0 7.0 + -2.0 -4.0 + +julia> F.Q * F.T * F.Z' +2×2 Array{Float64,2}: + 6.0 8.0 + -3.0 -5.0 + +julia> S, T, Q, Z, α, β = schurfact(A, B); # destructuring via iteration + +julia> Q * S * Z' ≈ A +true + +julia> Q * T * Z' ≈ B +true + +julia> α. / β ≈ eigvals(A, B) +true +``` """ schurfact(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:BlasFloat} = schurfact!(copy(A),copy(B)) function schurfact(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} @@ -300,16 +363,6 @@ end Base.propertynames(F::GeneralizedSchur) = (:values, :left, :right, fieldnames(typeof(F))...) -""" - schur(A::StridedMatrix, B::StridedMatrix) -> S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, α::Vector, β::Vector - -See [`schurfact`](@ref). -""" -function schur(A::StridedMatrix, B::StridedMatrix) - SchurF = schurfact(A, B) - SchurF.S, SchurF.T, SchurF.Q, SchurF.Z, SchurF.α, SchurF.β -end - function show(io::IO, mime::MIME{Symbol("text/plain")}, F::GeneralizedSchur) println(io, summary(F)) println(io, "S factor:") diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl index 2eefa52b69acc..b27ab8112edc4 100644 --- a/stdlib/LinearAlgebra/test/lapack.jl +++ b/stdlib/LinearAlgebra/test/lapack.jl @@ -590,7 +590,7 @@ end for job in ('N', 'E', 'V', 'B') for c in ('V', 'N') A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4]) - T,Q,d = schur(A) + T,Q,d = schurfact(A) s, sep = LinearAlgebra.LAPACK.trsen!(job,c,Array{LinearAlgebra.BlasInt}([0,1,0,0]),T,Q)[4:5] @test d[1] ≈ T[2,2] @test d[2] ≈ T[1,1] @@ -616,7 +616,7 @@ end @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) for c in ('V', 'N') A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4]) - T,Q,d = schur(A) + T,Q,d = schurfact(A) LinearAlgebra.LAPACK.trexc!(c,LinearAlgebra.BlasInt(1),LinearAlgebra.BlasInt(2),T,Q) @test d[1] ≈ T[2,2] @test d[2] ≈ T[1,1] diff --git a/stdlib/LinearAlgebra/test/schur.jl b/stdlib/LinearAlgebra/test/schur.jl index 2dac81d3b201c..beb6fdc71543c 100644 --- a/stdlib/LinearAlgebra/test/schur.jl +++ b/stdlib/LinearAlgebra/test/schur.jl @@ -35,15 +35,15 @@ aimg = randn(n,n)/2 @test convert(Array, f) ≈ a @test_throws ErrorException f.A - sch, vecs, vals = schur(UpperTriangular(triu(a))) + sch, vecs, vals = schurfact(UpperTriangular(triu(a))) @test vecs*sch*vecs' ≈ triu(a) - sch, vecs, vals = schur(LowerTriangular(tril(a))) + sch, vecs, vals = schurfact(LowerTriangular(tril(a))) @test vecs*sch*vecs' ≈ tril(a) - sch, vecs, vals = schur(Hermitian(asym)) + sch, vecs, vals = schurfact(Hermitian(asym)) @test vecs*sch*vecs' ≈ asym - sch, vecs, vals = schur(Symmetric(a + transpose(a))) + sch, vecs, vals = schurfact(Symmetric(a + transpose(a))) @test vecs*sch*vecs' ≈ a + transpose(a) - sch, vecs, vals = schur(Tridiagonal(a + transpose(a))) + sch, vecs, vals = schurfact(Tridiagonal(a + transpose(a))) @test vecs*sch*vecs' ≈ Tridiagonal(a + transpose(a)) tstring = sprint((t, s) -> show(t, "text/plain", s), f.T) @@ -111,7 +111,7 @@ aimg = randn(n,n)/2 @test S.Z ≈ SchurNew.Z @test S.alpha ≈ SchurNew.alpha @test S.beta ≈ SchurNew.beta - sS,sT,sQ,sZ = schur(a1_sf,a2_sf) + sS,sT,sQ,sZ = schurfact(a1_sf,a2_sf) @test NS.Q ≈ sQ @test NS.T ≈ sT @test NS.S ≈ sS @@ -119,7 +119,7 @@ aimg = randn(n,n)/2 end end @testset "0x0 matrix" for A in (zeros(eltya, 0, 0), view(rand(eltya, 2, 2), 1:0, 1:0)) - T, Z, λ = LinearAlgebra.schur(A) + T, Z, λ = LinearAlgebra.schurfact(A) @test T == A @test Z == A @test λ == zeros(0)