From 8fedd692c334528d5d7a88507295ba1050499aeb Mon Sep 17 00:00:00 2001 From: Sacha Verweij Date: Tue, 15 May 2018 15:21:44 -0700 Subject: [PATCH] Deprecate svd(...) in favor of svdfact(...) and factorization destructuring. --- NEWS.md | 40 ++++++ doc/src/manual/noteworthy-differences.md | 2 +- doc/src/manual/parallel-computing.md | 4 +- doc/src/manual/performance-tips.md | 4 +- stdlib/IterativeEigensolvers/test/runtests.jl | 4 +- stdlib/LinearAlgebra/docs/src/index.md | 21 ++- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 - stdlib/LinearAlgebra/src/bitarray.jl | 2 +- stdlib/LinearAlgebra/src/deprecated.jl | 73 ++++++++++ stdlib/LinearAlgebra/src/diagonal.jl | 8 +- stdlib/LinearAlgebra/src/svd.jl | 132 +++++------------- stdlib/LinearAlgebra/src/triangular.jl | 2 +- stdlib/LinearAlgebra/test/bidiag.jl | 6 +- stdlib/LinearAlgebra/test/diagonal.jl | 2 +- stdlib/LinearAlgebra/test/lapack.jl | 2 +- stdlib/LinearAlgebra/test/svd.jl | 4 +- stdlib/LinearAlgebra/test/triangular.jl | 1 - 17 files changed, 176 insertions(+), 132 deletions(-) diff --git a/NEWS.md b/NEWS.md index 11b4ac06c4c16f..f6b1764f23b26f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1040,6 +1040,46 @@ Deprecated or removed (`S, T, Q, Z, α, β = schurfact(A, B)`) or as a `GeneralizedSchur` object (`schurf = schurfact(A, B)`) ([#26997]). + * `svd(A::Abstractarray; thin=true)` has been deprecated in favor of + `svdfact(A; full=false)`. Note that the `thin` keyword and its replacement + `full` have opposite meanings. Additionally, whereas `svd` returns a + tuple of arrays `(U, S, V)` such that `A ≈ U*Diagonal(S)*V'`, `svdfact` returns + an `SVD` object that nominally provides `U, S, Vt` such that `A ≈ U*Diagonal(S)*Vt`. + So for a direct replacement, + use `(U, S, Vt = svdfact(A[; full=...]); (U, S, copy(Vt')))`. But going forward, + consider using the direct result of `svdfact(A[; full=...])` instead, + either destructured into its components (`U, S, Vt = svdfact(A[; full=...])`) + or as an `SVD` object (`svdf = svdfact(A[; full=...])`) ([#26997]). + + * `svd(x::Number; thin=true)` has been deprecated in favor of + `svdfact(x; full=false)`. Note that the `thin` keyword and its replacement + `full` have opposite meanings. Additionally, whereas `svd(x::Number[; thin=...])` + returns a tuple of numbers `(u, s, v)` such that `x ≈ u*s*conj(v)`, + `svdfact(x::Number[; full=...])` returns + an `SVD` object that nominally provides `u, s, vt` such that `x ≈ u*Diagonal(s)*conj(vt)`. + So for a direct replacement, + use `(u, s, vt = first.((svdfact(x[; full=...]...,)); (u, s, conj(vt)))`. + But going forward, + consider using the direct result of `svdfact(x[; full=...])` instead, + either destructured into its components (`U, S, Vt = svdfact(A[; full=...])`) + or as an `SVD` object (`svdf = svdfact(x[; full=...])`) ([#26997]). + + * `svd(A::Abstractarray, B::AbstractArray)` has been deprecated in favor of + `svdfact(A, B)`. Whereas the former returns a tuple of arrays, + the latter returns a `GeneralizedSVD` object. So for a direct replacement, + use `(svdfact(A, B)...,)`. But going forward, + consider using the direct result of `svdfact(A, B)` instead, + either destructured into its components (`U, V, Q, D1, D2, R0 = svdfact(A, B)`) + or as a `GeneralizedSVD` object (`gsvdf = svdfact(A, B)`) ([#26997]). + + * `svd(x::Number, y::Number)` has been deprecated in favor of + `svdfact(x, y)`. Whereas the former returns a tuple of numbers, + the latter returns a `GeneralizedSVD` object. So for a direct replacement, + use `first.((svdfact(x, y)...,))`. But going forward, + consider using the direct result of `svdfact(x, ys)` instead, + either destructured into its components (`U, V, Q, D1, D2, R0 = svdfact(x, y)`) + or as a `GeneralizedSVD` object (`gsvdf = svdfact(x, y)`) ([#26997]). + * The timing functions `tic`, `toc`, and `toq` are deprecated in favor of `@time` and `@elapsed` ([#17046]). diff --git a/doc/src/manual/noteworthy-differences.md b/doc/src/manual/noteworthy-differences.md index 20ba324c098e32..80941d9465739d 100644 --- a/doc/src/manual/noteworthy-differences.md +++ b/doc/src/manual/noteworthy-differences.md @@ -66,7 +66,7 @@ may trip up Julia users accustomed to MATLAB: parentheses may be required (e.g., to select elements of `A` equal to 1 or 2 use `(A .== 1) .| (A .== 2)`). * In Julia, the elements of a collection can be passed as arguments to a function using the splat operator `...`, as in `xs=[1,2]; f(xs...)`. - * Julia's [`svd`](@ref) returns singular values as a vector instead of as a dense diagonal matrix. + * Julia's [`svdfact`](@ref) returns singular values as a vector instead of as a dense diagonal matrix. * In Julia, `...` is not used to continue lines of code. Instead, incomplete expressions automatically continue onto the next line. * In both Julia and MATLAB, the variable `ans` is set to the value of the last expression issued diff --git a/doc/src/manual/parallel-computing.md b/doc/src/manual/parallel-computing.md index 76bb41af49448c..bcfb77374d8847 100644 --- a/doc/src/manual/parallel-computing.md +++ b/doc/src/manual/parallel-computing.md @@ -457,7 +457,7 @@ we could compute the singular values of several large random matrices in paralle ```julia-repl julia> M = Matrix{Float64}[rand(1000,1000) for i = 1:10]; -julia> pmap(svd, M); +julia> pmap(svdvals, M); ``` Julia's [`pmap`](@ref) is designed for the case where each function call does a large amount @@ -486,7 +486,7 @@ As an example, consider computing the singular values of matrices of different s ```julia-repl julia> M = Matrix{Float64}[rand(800,800), rand(600,600), rand(800,800), rand(600,600)]; -julia> pmap(svd, M); +julia> pmap(svdvals, M); ``` If one process handles both 800×800 matrices and another handles both 600×600 matrices, we will diff --git a/doc/src/manual/performance-tips.md b/doc/src/manual/performance-tips.md index ff8ecd570ef485..1ecac2d6baadee 100644 --- a/doc/src/manual/performance-tips.md +++ b/doc/src/manual/performance-tips.md @@ -544,7 +544,7 @@ function norm(A) if isa(A, Vector) return sqrt(real(dot(A,A))) elseif isa(A, Matrix) - return maximum(svd(A)[2]) + return maximum(svdvals(A)) else error("norm: invalid argument") end @@ -555,7 +555,7 @@ This can be written more concisely and efficiently as: ```julia norm(x::Vector) = sqrt(real(dot(x,x))) -norm(A::Matrix) = maximum(svd(A)[2]) +norm(A::Matrix) = maximum(svdvals(A)) ``` ## Write "type-stable" functions diff --git a/stdlib/IterativeEigensolvers/test/runtests.jl b/stdlib/IterativeEigensolvers/test/runtests.jl index 6cbad717dec5d5..f02e464a227224 100644 --- a/stdlib/IterativeEigensolvers/test/runtests.jl +++ b/stdlib/IterativeEigensolvers/test/runtests.jl @@ -186,7 +186,7 @@ end @testset "real svds" begin A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], [2.0, -1.0, 6.1, 7.0, 1.5]) S1 = svds(A, nsv = 2) - S2 = svd(Array(A)) + S2 = (F = svdfact(Array(A)); (F.U, F.S, F.V)) ## singular values match: @test S1[1].S ≈ S2[2][1:2] @@ -248,7 +248,7 @@ end @testset "complex svds" begin A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], exp.(im*[2.0:2:10;]), 5, 4) S1 = svds(A, nsv = 2) - S2 = svd(Array(A)) + S2 = (F = svdfact(Array(A)); (F.U, F.S, F.V)) ## singular values match: @test S1[1].S ≈ S2[2][1:2] diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 7ab9c2fb63fc44..16f5423a188ced 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 | [`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 | | | | +| Matrix type | LAPACK | [`eigfact`](@ref) | [`eigvals`](@ref) | [`eigvecs`](@ref) | [`svdfact`](@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: @@ -349,7 +349,6 @@ LinearAlgebra.ordschur LinearAlgebra.ordschur! LinearAlgebra.svdfact LinearAlgebra.svdfact! -LinearAlgebra.svd LinearAlgebra.svdvals LinearAlgebra.svdvals! LinearAlgebra.Givens diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 156d20f4e94419..985103e6c68471 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -133,7 +133,6 @@ export rdiv!, schurfact!, schurfact, - svd, svdfact!, svdfact, svdvals!, diff --git a/stdlib/LinearAlgebra/src/bitarray.jl b/stdlib/LinearAlgebra/src/bitarray.jl index ea127ddfeebf4a..995c7ce73b087e 100644 --- a/stdlib/LinearAlgebra/src/bitarray.jl +++ b/stdlib/LinearAlgebra/src/bitarray.jl @@ -87,7 +87,7 @@ end ## norm and rank -svd(A::BitMatrix) = svd(float(A)) +svdfact(A::BitMatrix) = svdfact(float(A)) qr(A::BitMatrix) = qr(float(A)) ## kron diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 7ed0abcf42aa2f..b67aaab244b40a 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -1345,3 +1345,76 @@ function schur(A::StridedMatrix, B::StridedMatrix) " or as a `GeneralizedSchur` object (`schurf = schurfact(A, B)`)."), :schur) return (schurfact(A, B)...,) end + +# deprecate svd(...) in favor of svdfact(...) and factorization destructuring +export svd +function svd(A::AbstractArray; full::Bool = false, thin::Union{Bool,Nothing} = nothing) + depwarn(string("`svd(A::Abstractarray; thin=true)` has been deprecated ", + "in favor of `svdfact(A; full=false)`. Note that the `thin` keyword ", + "and its replacement `full` have opposite meanings. Additionally, whereas ", + "`svd` returns a tuple of arrays `(U, S, V)` such that `A ≈ U*Diagonal(S)*V'`, ", + "`svdfact` returns an `SVD` object that nominally provides `U, S, Vt` ", + "such that `A ≈ U*Diagonal(S)*Vt`. So for a direct replacement, use ", + "`(U, S, Vt = svdfact(A[; full=...]); (U, S, copy(Vt')))`. But going forward, ", + "consider using the direct result of `svdfact(A[; full=...])` instead, ", + "either destructured into its components (`U, S, Vt = svdfact(A[; full=...])`) ", + "or as an `SVD` object (`svdf = svdfact(A[; full=...])`)."), :svd) + F = svdfact(A, full = (thin != nothing ? !thin : full)) + return F.U, F.S, copy(F.Vt') +end +function svd(x::Number; full::Bool = false, thin::Union{Bool,Nothing} = nothing) + depwarn(string("`svd(x::Number; thin=true)` has been deprecated ", + "in favor of `svdfact(x; full=false)`. Note that the `thin` keyword ", + "and its replacement `full` have opposite meanings. Additionally, whereas ", + "`svd(x::Number[; thin=...])` returns a tuple of numbers `(u, s, v)` ", + " such that `x ≈ u*s*conj(v)`, `svdfact(x::Number[; full=...])` returns ", + "an `SVD` object that nominally provides `u, s, vt` such that ", + "`x ≈ u*s*vt`. So for a direct replacement, use ", + "`(u, s, vt = first.((svdfact(A[; full=...])...,)); (u, s, conj(vt)))`. ", + "But going forward, ", + "consider using the direct result of `svdfact(A[; full=...])` instead, ", + "either destructured into its components (`u, s, vt = svdfact(A[; full=...])`) ", + "or as an `SVD` object (`svdf = svdfact(A[; full=...])`)."), :svd) + u, s, v = first.((svdfact(x)...,)) + return u, s, conj(vt) +end + +function svd(A::AbstractMatrix, B::AbstractMatrix) + depwarn(string("`svd(A::Abstractarray, B::AbstractArray)` has been deprecated ", + "in favor of `svdfact(A, B)`. Whereas the former returns a tuple of arrays, ", + "the latter returns a `GeneralizedSVD` object. So for a direct replacement, ", + "use `(svdfact(A, B)...,)`. But going forward, ", + "consider using the direct result of `svdfact(A, B)` instead, ", + "either destructured into its components ", + "(`U, V, Q, D1, D2, R0 = svdfact(A, B)`) ", + "or as a `GeneralizedSVD` object (`gsvdf = svdfact(A, B)`)."), :svd) + return (svdfact(A, B)...,) +end +function svd(x::Number, y::Number) + depwarn(string("`svd(x::Number, y::Number)` has been deprecated ", + "in favor of `svdfact(x, y)`. Whereas the former returns a tuple of numbers, ", + "the latter returns a `GeneralizedSVD` object. So for a direct replacement, ", + "use `first.((svdfact(x, y)...,))`. But going forward, ", + "consider using the direct result of `svdfact(x, y)` instead, ", + "either destructured into its components ", + "(`U, V, Q, D1, D2, R0 = svdfact(x, y)`) ", + "or as a `GeneralizedSVD` object (`gsvdf = svdfact(x, y)`)."), :svd) + return first.((svdfact(x, y)...,)) +end + +@inline function _simpledepsvd(A) + depwarn(string("`svd(A)` has been deprecated ", + "in favor of `svdfact(A)`. Whereas `svd` ", + "returns a tuple of arrays `(U, S, V)` such that `A ≈ U*Diagonal(S)*V'`, ", + "`svdfact` returns an `SVD` object that nominally provides `U, S, Vt` ", + "such that `A ≈ U*Diagonal(S)*Vt`. So for a direct replacement, use ", + "`(U, S, Vt = svdfact(A); (U, S, copy(Vt')))`. But going forward, ", + "consider using the direct result of `svdfact(A)` instead, ", + "either destructured into its components (`U, S, Vt = svdfact(A)`) ", + "or as an `SVD` object (`svdf = svdfact(A)`)."), :svd) + U, S, Vt = svdfact(A) + return U, S, copy(Vt') +end +svd(A::BitMatrix) = _simpledepsvd(float(A)) +svd(D::Diagonal{<:Number}) = _simpledepsvd(D) +svd(A::AbstractTriangular) = _simpledepsvd(copyto!(similar(parent(A)), A)) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 95406701e30482..bfcb0aeb6c86bb 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -455,18 +455,18 @@ end #Singular system svdvals(D::Diagonal{<:Number}) = sort!(abs.(D.diag), rev = true) svdvals(D::Diagonal) = [svdvals(v) for v in D.diag] -function svd(D::Diagonal{<:Number}) +function svdfact(D::Diagonal{<:Number}) S = abs.(D.diag) piv = sortperm(S, rev = true) U = Diagonal(D.diag ./ S) Up = hcat([U[:,i] for i = 1:length(D.diag)][piv]...) V = Diagonal(fill!(similar(D.diag), one(eltype(D.diag)))) Vp = hcat([V[:,i] for i = 1:length(D.diag)][piv]...) - return (Up, S[piv], Vp) + return SVD(Up, S[piv], copy(Vp')) end function svdfact(D::Diagonal) - U, s, V = svd(D) - SVD(U, s, copy(V')) + U, s, Vt = svdfact(D) + return SVD(U, s, Vt) end # dismabiguation methods: * of Diagonal and Adj/Trans AbsVec diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index 0a2070633b5c38..d8a26e93481a30 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -10,6 +10,14 @@ struct SVD{T,Tr,M<:AbstractArray} <: Factorization{T} end SVD(U::AbstractArray{T}, S::Vector{Tr}, Vt::AbstractArray{T}) where {T,Tr} = SVD{T,Tr,typeof(U)}(U, S, Vt) +# iteration for destructuring into factors +Base.start(::SVD) = Val(:U) +Base.next(F::SVD, ::Val{:U}) = (F.U, Val(:S)) +Base.next(F::SVD, ::Val{:S}) = (F.S, Val(:Vt)) +Base.next(F::SVD, ::Val{:Vt}) = (F.Vt, Val(:done)) +Base.done(F::SVD, ::Val{:done}) = true +Base.done(F::SVD, ::Any) = false + """ svdfact!(A; full::Bool = false) -> SVD @@ -92,6 +100,11 @@ julia> F.U * Diagonal(F.S) * F.Vt 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 + +julia> U, S, Vt = svdfact(A); # destructuring via iteration + +julia> U*Diagonal(S)*Vt ≈ A +true ``` """ function svdfact(A::StridedVecOrMat{T}; full::Bool = false, thin::Union{Bool,Nothing} = nothing) where T @@ -125,62 +138,6 @@ function svdfact(x::Integer; full::Bool = false, thin::Union{Bool,Nothing} = not return svdfact(float(x), full = full) end -""" - svd(A; full::Bool = false) -> U, S, V - -Computes the SVD of `A`, returning `U`, vector `S`, and `V` such that -`A == U * Diagonal(S) * V'`. The singular values in `S` are sorted in descending order. - -If `full = false` (default), a "thin" SVD is returned. For a ``M -\\times N`` matrix `A`, in the full factorization `U` is `M \\times M` -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. - -`svd` is a wrapper around [`svdfact`](@ref), extracting all parts -of the `SVD` factorization to a tuple. Direct use of `svdfact` is therefore more -efficient. - -# Examples -```jldoctest -julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.] -4×5 Array{Float64,2}: - 1.0 0.0 0.0 0.0 2.0 - 0.0 0.0 3.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 2.0 0.0 0.0 0.0 - -julia> U, S, V = svd(A); - -julia> U * Diagonal(S) * V' -4×5 Array{Float64,2}: - 1.0 0.0 0.0 0.0 2.0 - 0.0 0.0 3.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 2.0 0.0 0.0 0.0 -``` -""" -function svd(A::AbstractArray; full::Bool = false, thin::Union{Bool,Nothing} = nothing) - # DEPRECATION TODO: remove deprecated thin argument and associated logic after 0.7 - if thin != nothing - Base.depwarn(string("the `thin` keyword argument in `svd(A; thin = $(thin))` has ", - "been deprecated in favor of `full`, which has the opposite meaning, ", - "e.g `svd(A; full = $(!thin))`."), :svd) - full::Bool = !thin - end - F = svdfact(A, full = full) - F.U, F.S, copy(F.Vt') -end -function svd(x::Number; full::Bool = false, thin::Union{Bool,Nothing} = nothing) - # DEPRECATION TODO: remove deprecated thin argument and associated logic after 0.7 - if thin != nothing - Base.depwarn(string("the `thin` keyword argument in `svd(A; thin = $(thin))` has ", - "been deprecated in favor of `full`, which has the opposite meaning, ", - "e.g. `svd(A; full = $(!thin))`."), :svd) - full::Bool = !thin - end - return first.(svd(fill(x, 1, 1))) -end function getproperty(F::SVD, d::Symbol) if d == :V @@ -278,6 +235,17 @@ function GeneralizedSVD(U::AbstractMatrix{T}, V::AbstractMatrix{T}, Q::AbstractM GeneralizedSVD{T,typeof(U)}(U, V, Q, a, b, k, l, R) end +# iteration for destructuring into factors +Base.start(::GeneralizedSVD) = Val(:U) +Base.next(F::GeneralizedSVD, ::Val{:U}) = (F.U, Val(:V)) +Base.next(F::GeneralizedSVD, ::Val{:V}) = (F.V, Val(:Q)) +Base.next(F::GeneralizedSVD, ::Val{:Q}) = (F.Q, Val(:D1)) +Base.next(F::GeneralizedSVD, ::Val{:D1}) = (F.D1, Val(:D2)) +Base.next(F::GeneralizedSVD, ::Val{:D2}) = (F.D2, Val(:R0)) +Base.next(F::GeneralizedSVD, ::Val{:R0}) = (F.R0, Val(:done)) +Base.done(F::GeneralizedSVD, ::Val{:done}) = true +Base.done(F::GeneralizedSVD, ::Any) = false + """ svdfact!(A, B) -> GeneralizedSVD @@ -341,10 +309,10 @@ For an M-by-N matrix `A` and P-by-N matrix `B`, - `U` is a M-by-M orthogonal matrix, - `V` is a P-by-P orthogonal matrix, - `Q` is a N-by-N orthogonal matrix, -- `R0` is a (K+L)-by-N matrix whose rightmost (K+L)-by-(K+L) block is - nonsingular upper block triangular, - `D1` is a M-by-(K+L) diagonal matrix with 1s in the first K entries, - `D2` is a P-by-(K+L) matrix whose top right L-by-L block is diagonal, +- `R0` is a (K+L)-by-N matrix whose rightmost (K+L)-by-(K+L) block is + nonsingular upper block triangular, `K+L` is the effective numerical rank of the matrix `[A; B]`. @@ -377,6 +345,14 @@ julia> F.V*F.D2*F.R0*F.Q' 2×2 Array{Float64,2}: 0.0 1.0 1.0 0.0 + +julia> U, V, Q, D1, D2, R0 = svdfact(A, B); # destructuring via iteration + +julia> U*D1*R0*Q' ≈ A +true + +julia> V*D2*R0*Q' ≈ B +true ``` """ function svdfact(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} @@ -388,46 +364,6 @@ end # version svdfact(x::Number, y::Number) = svdfact(fill(x, 1, 1), fill(y, 1, 1)) -""" - svd(A, B) -> U, V, Q, D1, D2, R0 - -Wrapper around [`svdfact`](@ref) extracting all parts of the -factorization to a tuple. Direct use of -`svdfact` is therefore generally more efficient. The function returns the generalized SVD of -`A` and `B`, returning `U`, `V`, `Q`, `D1`, `D2`, and `R0` such that `A = U*D1*R0*Q'` and `B = -V*D2*R0*Q'`. - -# Examples -```jldoctest -julia> A = [1. 0.; 0. -1.] -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 -1.0 - -julia> B = [0. 1.; 1. 0.] -2×2 Array{Float64,2}: - 0.0 1.0 - 1.0 0.0 - -julia> U, V, Q, D1, D2, R0 = svd(A, B); - -julia> U*D1*R0*Q' -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 -1.0 - -julia> V*D2*R0*Q' -2×2 Array{Float64,2}: - 0.0 1.0 - 1.0 0.0 -``` -""" -function svd(A::AbstractMatrix, B::AbstractMatrix) - F = svdfact(A, B) - F.U, F.V, F.Q, F.D1, F.D2, F.R0 -end -svd(x::Number, y::Number) = first.(svd(fill(x, 1, 1), fill(y, 1, 1))) - @inline function getproperty(F::GeneralizedSVD{T}, d::Symbol) where T Fa = getfield(F, :a) Fb = getfield(F, :b) diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 2bdabd65a85b6b..6bfc450c0e5eb0 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -2435,7 +2435,7 @@ end eigfact(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A)) # Generic singular systems -for func in (:svd, :svdfact, :svdfact!, :svdvals) +for func in (:svdfact, :svdfact!, :svdvals) @eval begin ($func)(A::AbstractTriangular) = ($func)(copyto!(similar(parent(A)), A)) end diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index 6d83933ba248b9..c48d52c12afdf9 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -244,8 +244,8 @@ srand(1) if (elty <: BlasReal) @test AbstractArray(svdfact(T)) ≈ AbstractArray(svdfact!(copy(Tfull))) @test svdvals(Tfull) ≈ svdvals(T) - u1, d1, v1 = svd(Tfull) - u2, d2, v2 = svd(T) + u1, d1, v1 = (F = svdfact(Tfull); (F.U, F.S, copy(F.Vt'))) + u2, d2, v2 = (F = svdfact(T); (F.U, F.S, copy(F.Vt'))) @test d1 ≈ d2 if elty <: Real test_approx_eq_modphase(u1, u2) @@ -253,7 +253,7 @@ srand(1) end @test 0 ≈ vecnorm(u2*Diagonal(d2)*v2'-Tfull) atol=n*max(n^2*eps(relty),vecnorm(u1*Diagonal(d1)*v1'-Tfull)) @inferred svdvals(T) - @inferred svd(T) + @inferred svdfact(T) end end diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 06a5d78ba2c84d..f44ec2e24d092b 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -262,7 +262,7 @@ srand(1) end @testset "svd (#11120/#11247)" begin - U, s, V = svd(D) + U, s, V = (F = svdfact(D); (F.U, F.S, F.V)) @test (U*Diagonal(s))*V' ≈ D @test svdvals(D) == s @test svdfact(D).V == V diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl index b27ab8112edc4d..ad80e1b322b991 100644 --- a/stdlib/LinearAlgebra/test/lapack.jl +++ b/stdlib/LinearAlgebra/test/lapack.jl @@ -165,7 +165,7 @@ end @testset "gesvd, ggsvd" begin @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) A = rand(elty,10,5) - U,S,V = svd(A) + U,S,V = (F = svdfact(A); (F.U, F.S, F.V)) lU,lS,lVt = LAPACK.gesvd!('S','S',A) @test U ≈ lU @test S ≈ lS diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index 3af0a4f2e9e040..a719cef86d1c86 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -82,7 +82,7 @@ a2img = randn(n,n)/2 β = svdfact(α) @test β.S == [abs(α)] @test svdvals(α) == abs(α) - u,v,q,d1,d2,r0 = svd(a,a_svd) + u,v,q,d1,d2,r0 = svdfact(a,a_svd) @test u ≈ gsvd.U @test v ≈ gsvd.V @test d1 ≈ gsvd.D1 @@ -103,10 +103,8 @@ a2img = randn(n,n)/2 x, y = randn(eltya, 2) @test svdfact(x) == svdfact(fill(x, 1, 1)) @test svdvals(x) == first(svdvals(fill(x, 1, 1))) - @test svd(x) == first.(svd(fill(x, 1, 1))) @test svdfact(x, y) == svdfact(fill(x, 1, 1), fill(y, 1, 1)) @test svdvals(x, y) ≈ first(svdvals(fill(x, 1, 1), fill(y, 1, 1))) - @test svd(x, y) == first.(svd(fill(x, 1, 1), fill(y, 1, 1))) end end if eltya != Int diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 49c8a13f1f726e..ea77322e304d71 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -264,7 +264,6 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo end if !(elty1 in (BigFloat, Complex{BigFloat})) # Not implemented yet - svd(A1) svdfact(A1) elty1 <: BlasFloat && svdfact!(copy(A1)) svdvals(A1)