From cec4c3210986ba045f0d8ecf6396647337361d8e Mon Sep 17 00:00:00 2001 From: Alan Edelman Date: Tue, 7 Jan 2020 06:30:22 -0500 Subject: [PATCH] Update svd.jl (#30239) * Update svd.jl * Update svd.jl Make doc more useful to the reader. * Fix grammar and one trailing whitespace * Update svd.jl * Fix doctests * Avoid duplicating documentation between mutating and non-mutating versions of svd functions. * Reorganize generalized svd docstring and doctests * Break some long lines Co-authored-by: Viral B. Shah Co-authored-by: Andreas Noack --- stdlib/LinearAlgebra/src/svd.jl | 189 +++++++------------------------- 1 file changed, 38 insertions(+), 151 deletions(-) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index 843235a615d3b..ccf025be062e4 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -87,32 +87,7 @@ default_svd_alg(A) = DivideAndConquer() svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD `svd!` is the same as [`svd`](@ref), but saves space by -overwriting the input `A`, instead of creating a copy. - -# 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> F = svd!(A); - -julia> F.U * Diagonal(F.S) * F.Vt -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> A -4×5 Array{Float64,2}: - -2.23607 0.0 0.0 0.0 0.618034 - 0.0 -3.0 1.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 -2.0 0.0 0.0 +overwriting the input `A`, instead of creating a copy. See documentation of [`svd`](@ref) for details. ``` """ function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where T<:BlasFloat @@ -161,25 +136,21 @@ Another (typically slower but more accurate) option is `alg = QRIteration()`. # 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> A = rand(4,3); -julia> F = svd(A); +julia> F = svd(A); # Store the Factorization Object -julia> F.U * Diagonal(F.S) * F.Vt -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> A ≈ F.U * Diagonal(F.S) * F.Vt +true -julia> u, s, v = F; # destructuring via iteration +julia> U, S, V = F; # destructuring via iteration -julia> u == F.U && s == F.S && v == F.V +julia> A ≈ U * Diagonal(S) * V' +true + +julia> Uonly, = svd(A); # Store U only + +julia> Uonly == U true ``` """ @@ -217,29 +188,6 @@ Base.propertynames(F::SVD, private::Bool=false) = Return the singular values of `A`, saving space by overwriting the input. See also [`svdvals`](@ref) and [`svd`](@ref). - -# 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> svdvals!(A) -4-element Array{Float64,1}: - 3.0 - 2.23606797749979 - 2.0 - 0.0 - -julia> A -4×5 Array{Float64,2}: - -2.23607 0.0 0.0 0.0 0.618034 - 0.0 -3.0 1.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 -2.0 0.0 0.0 ``` """ svdvals!(A::StridedMatrix{T}) where {T<:BlasFloat} = isempty(A) ? zeros(real(T), 0) : LAPACK.gesdd!('N', A)[2] @@ -409,41 +357,7 @@ Base.iterate(S::GeneralizedSVD, ::Val{:done}) = nothing svd!(A, B) -> GeneralizedSVD `svd!` is the same as [`svd`](@ref), but modifies the arguments -`A` and `B` in-place, instead of making copies. - -# 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> F = svd!(A, B); - -julia> F.U*F.D1*F.R0*F.Q' -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 -1.0 - -julia> F.V*F.D2*F.R0*F.Q' -2×2 Array{Float64,2}: - 0.0 1.0 - 1.0 0.0 - -julia> A -2×2 Array{Float64,2}: - 1.41421 0.0 - 0.0 -1.41421 - -julia> B -2×2 Array{Float64,2}: - 1.0 -0.0 - 0.0 -1.0 +`A` and `B` in-place, instead of making copies. See documentation of [`svd`](@ref) for details. ``` """ function svd!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasFloat @@ -458,12 +372,11 @@ end svd(A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat} = svd!(copy(A),copy(B)) """ + svd(A, B) -> GeneralizedSVD Compute the generalized SVD of `A` and `B`, returning a `GeneralizedSVD` factorization -object `F`, such that `A = F.U*F.D1*F.R0*F.Q'` and `B = F.V*F.D2*F.R0*F.Q'`. - -For an M-by-N matrix `A` and P-by-N matrix `B`, +object `F` such that `[A;B] = [F.U * F.D1; F.V * F.D2] * F.R0 * F.Q'` - `U` is a M-by-M orthogonal matrix, - `V` is a P-by-P orthogonal matrix, @@ -477,35 +390,36 @@ For an M-by-N matrix `A` and P-by-N matrix `B`, Iterating the decomposition produces the components `U`, `V`, `Q`, `D1`, `D2`, and `R0`. -The entries of `F.D1` and `F.D2` are related, as explained in the LAPACK -documentation for the -[generalized SVD](http://www.netlib.org/lapack/lug/node36.html) and the -[xGGSVD3](http://www.netlib.org/lapack/explore-html/d6/db3/dggsvd3_8f.html) -routine which is called underneath (in LAPACK 3.6.0 and newer). +The generalized SVD is used in applications such as when one wants to compare how much belongs +to `A` vs. how much belongs to `B`, as in human vs yeast genome, or signal vs noise, or between +clusters vs within clusters. (See Edelman and Wang for discussion: https://arxiv.org/abs/1901.00485) + +It decomposes `[A; B]` into `[UC; VS]H`, where `[UC; VS]` is a natural orthogonal basis for the +column space of `[A; B]`, and `H = RQ'` is a natural non-orthogonal basis for the rowspace of `[A;B]`, +where the top rows are most closely attributed to the `A` matrix, and the bottom to the `B` matrix. +The multi-cosine/sine matrices `C` and `S` provide a multi-measure of how much `A` vs how much `B`, +and `U` and `V` provide directions in which these are measured. # 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> A = randn(3,2); B=randn(4,2); julia> F = svd(A, B); -julia> F.U*F.D1*F.R0*F.Q' -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 -1.0 +julia> U,V,Q,C,S,R = F; -julia> F.V*F.D2*F.R0*F.Q' -2×2 Array{Float64,2}: - 0.0 1.0 - 1.0 0.0 +julia> H = R*Q'; + +julia> [A; B] ≈ [U*C; V*S]*H +true + +julia> [A; B] ≈ [F.U*F.D1; F.V*F.D2]*F.R0*F.Q' +true + +julia> Uonly, = svd(A,B); + +julia> U == Uonly +true ``` """ function svd(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} @@ -582,33 +496,6 @@ end Return the generalized singular values from the generalized singular value decomposition of `A` and `B`, saving space by overwriting `A` and `B`. See also [`svd`](@ref) and [`svdvals`](@ref). - -# 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> svdvals!(A, B) -2-element Array{Float64,1}: - 1.0 - 1.0 - -julia> A -2×2 Array{Float64,2}: - 1.41421 0.0 - 0.0 -1.41421 - -julia> B -2×2 Array{Float64,2}: - 1.0 -0.0 - 0.0 -1.0 ``` """ function svdvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasFloat