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 <viral@juliacomputing.com>
Co-authored-by: Andreas Noack <andreas@noack.dk>

Author: 3 people

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