Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update svd.jl #30239

Merged
merged 9 commits into from
Jan 7, 2020
Merged

Update svd.jl #30239

Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
2f8e42e
Update svd.jl
alanedelman Dec 2, 2018
bbbc8dd
Update svd.jl
alanedelman Dec 17, 2018
817fab3
Fix grammar and one trailing whitespace
ViralBShah Dec 20, 2018
3ecd192
Update svd.jl
alanedelman Jan 3, 2019
e980f79
Fix doctests
andreasnoack Jan 4, 2019
a5b1586
Avoid duplicating documentation between mutating and non-mutating
andreasnoack Jan 4, 2019
6606597
Reorganize generalized svd docstring and doctests
andreasnoack Jan 4, 2019
7b3f77e
Merge branch 'master' into patch-1
andreasnoack Jan 7, 2020
3826d18
Break some long lines
andreasnoack Jan 7, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
189 changes: 38 additions & 151 deletions stdlib/LinearAlgebra/src/svd.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
```
"""
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand All @@ -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}
Expand Down Expand Up @@ -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
Expand Down