From 235967782a3211810159e241c00fede2b6edbae1 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 13 May 2019 13:33:40 +0200 Subject: [PATCH 1/6] add svd for realsym/hermitian, and tests --- stdlib/LinearAlgebra/src/svd.jl | 9 +++++++++ stdlib/LinearAlgebra/test/svd.jl | 20 ++++++++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index 76d9fcb934633..d3db8d2c81b7d 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -119,6 +119,15 @@ function svd(A::Transpose; full::Bool = false) s = svd(A.parent, full = full) return SVD(transpose(s.Vt), s.S, transpose(s.U)) end +function svd(A::Union{Hermitian,Symmetric{<:Real}}, full::Bool=false) + s = eigen(A) + vals = s.values + σ = map(x -> x >= 0 ? 1 : -1, vals) + vals .*= σ + I = sortperm(vals; rev=true) + vecs = s.vectors[:,I] + return SVD(vecs, vals[I], (vecs*Diagonal(σ[I]))') +end function getproperty(F::SVD, d::Symbol) if d == :V diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index 71f9a41caeef1..857fb9335e18a 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -92,6 +92,26 @@ a2img = randn(n,n)/2 b = rand(eltya,n) @test usv\b ≈ transpose(a)\b end + usv = svd(Symmetric(real(asym))) + @testset "singular value decomposition of real Symmetric" begin + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ real(asym) + @test convert(Array, usv) ≈ real(asym) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ real(asym)\b + end + usv = svd(Hermitian(asym)) + @testset "singular value decomposition of Hermitian" begin + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ Hermitian(asym) + @test convert(Array, usv) ≈ Hermitian(asym) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ Hermitian(asym)\b + end @testset "Generalized svd" begin a_svd = a[1:n1, :] gsvd = svd(a,a_svd) From d0446d2f4ba1bb26303ff0989d6f714fe1eab9b7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 13 May 2019 17:56:34 +0200 Subject: [PATCH 2/6] shift improved code from svd.jl to symmetric.jl --- stdlib/LinearAlgebra/src/svd.jl | 9 --------- stdlib/LinearAlgebra/src/symmetric.jl | 11 +++++++++++ 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index d3db8d2c81b7d..76d9fcb934633 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -119,15 +119,6 @@ function svd(A::Transpose; full::Bool = false) s = svd(A.parent, full = full) return SVD(transpose(s.Vt), s.S, transpose(s.U)) end -function svd(A::Union{Hermitian,Symmetric{<:Real}}, full::Bool=false) - s = eigen(A) - vals = s.values - σ = map(x -> x >= 0 ? 1 : -1, vals) - vals .*= σ - I = sortperm(vals; rev=true) - vecs = s.vectors[:,I] - return SVD(vecs, vals[I], (vecs*Diagonal(σ[I]))') -end function getproperty(F::SVD, d::Symbol) if d == :V diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 14bae24d8c020..e1fe95a5c3d32 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -677,6 +677,17 @@ eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}) where {T<:BlasComplex,S<:StridedM eigvecs(A::HermOrSym) = eigvecs(eigen(A)) +function svd(A::RealHermSymComplexHerm, full::Bool=false) + vals, vecs = eigen(A) + σ = map(x -> x >= 0 ? 1 : -1, vals) + vals .= abs.(vals) + I = sortperm(vals; rev=true) + permute!(vals, I) + Σ = Diagonal(σ[I]) + Base.permutecols!!(vecs, I) + return SVD(vecs, vals, (vecs*Σ)') +end + function svdvals!(A::RealHermSymComplexHerm) vals = eigvals!(A) for i = 1:length(vals) From b951f311b92a19cddb509cc80fd54281515b0d46 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 14 May 2019 09:22:35 +0200 Subject: [PATCH 3/6] minimize allocations --- stdlib/LinearAlgebra/src/symmetric.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index e1fe95a5c3d32..2220ef250102c 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -679,13 +679,18 @@ eigvecs(A::HermOrSym) = eigvecs(eigen(A)) function svd(A::RealHermSymComplexHerm, full::Bool=false) vals, vecs = eigen(A) - σ = map(x -> x >= 0 ? 1 : -1, vals) - vals .= abs.(vals) - I = sortperm(vals; rev=true) + I = sortperm(vals; by=abs, rev=true) permute!(vals, I) - Σ = Diagonal(σ[I]) - Base.permutecols!!(vecs, I) - return SVD(vecs, vals, (vecs*Σ)') + Base.permutecols!!(vecs, I) # left-singular vectors + V = copy(vecs) # right-singular vectors + # shifting -1 from singular values to right-singular vectors + @inbounds for i = 1:length(vals) + if vals[i] < 0 + vals[i] = -vals[i] + for j = 1:size(V,1); V[j,i] = -V[j,i]; end + end + end + return SVD(vecs, vals, V') end function svdvals!(A::RealHermSymComplexHerm) From c7b9ab3a96ba9fecb15e7c94060cf1ad1cd30ec8 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 May 2019 14:40:05 +0200 Subject: [PATCH 4/6] combine tests --- stdlib/LinearAlgebra/test/svd.jl | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index 857fb9335e18a..c21475f50b5c1 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -92,25 +92,17 @@ a2img = randn(n,n)/2 b = rand(eltya,n) @test usv\b ≈ transpose(a)\b end - usv = svd(Symmetric(real(asym))) - @testset "singular value decomposition of real Symmetric" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ real(asym) - @test convert(Array, usv) ≈ real(asym) - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ real(asym)\b - end - usv = svd(Hermitian(asym)) - @testset "singular value decomposition of Hermitian" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ Hermitian(asym) - @test convert(Array, usv) ≈ Hermitian(asym) - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ Hermitian(asym)\b + @testset "singular value decomposition of Hermitian/real-Symmetric" begin + for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,)) + usv = svd(T(asym)) + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ T(asym) + @test convert(Array, usv) ≈ T(asym) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ T(asym)\b + end end @testset "Generalized svd" begin a_svd = a[1:n1, :] From 829469df26c0dc8a0f63d9d2c24e39be53e40eb8 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 May 2019 15:31:27 +0200 Subject: [PATCH 5/6] move test out of loop --- stdlib/LinearAlgebra/test/svd.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index c21475f50b5c1..ed445cf7eb7fc 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -92,18 +92,6 @@ a2img = randn(n,n)/2 b = rand(eltya,n) @test usv\b ≈ transpose(a)\b end - @testset "singular value decomposition of Hermitian/real-Symmetric" begin - for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,)) - usv = svd(T(asym)) - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ T(asym) - @test convert(Array, usv) ≈ T(asym) - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ T(asym)\b - end - end @testset "Generalized svd" begin a_svd = a[1:n1, :] gsvd = svd(a,a_svd) @@ -133,6 +121,18 @@ a2img = randn(n,n)/2 @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ c end end + @testset "singular value decomposition of Hermitian/real-Symmetric" begin + for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,)) + usv = svd(T(asym)) + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ T(asym) + @test convert(Array, usv) ≈ T(asym) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ T(asym)\b + end + end if eltya <: LinearAlgebra.BlasReal @testset "Number input" begin x, y = randn(eltya, 2) From bfc2044ff7e490827b947c40b32909e5727fe50a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 May 2019 17:12:47 +0200 Subject: [PATCH 6/6] clean up svd tests --- stdlib/LinearAlgebra/test/svd.jl | 44 ++++++++++---------------------- 1 file changed, 13 insertions(+), 31 deletions(-) diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index ed445cf7eb7fc..7a29d86974a61 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -35,25 +35,15 @@ end n = 10 -# Split n into 2 parts for tests needing two matrices -n1 = div(n, 2) -n2 = 2*n1 - Random.seed!(1234321) areal = randn(n,n)/2 aimg = randn(n,n)/2 -a2real = randn(n,n)/2 -a2img = randn(n,n)/2 @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) - aa2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) asym = aa' + aa # symmetric indefinite - apd = aa' * aa # symmetric positive-definite - for (a, a2) in ((aa, aa2), (view(aa, 1:n, 1:n), view(aa2, 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) - + for a in (aa, view(aa, 1:n, 1:n)) usv = svd(a) @testset "singular value decomposition" begin @test usv.S === svdvals(usv) @@ -72,28 +62,20 @@ a2img = randn(n,n)/2 @test svdz.Vt ≈ Matrix{eltya}(I, 0, 0) end end - usv = svd(a') - @testset "singular value decomposition of adjoint" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ a' - @test convert(Array, usv) ≈ a' - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ a'\b - end - usv = svd(transpose(a)) - @testset "singular value decomposition of transpose" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ transpose(a) - @test convert(Array, usv) ≈ transpose(a) - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ transpose(a)\b + @testset "singular value decomposition of adjoint/transpose" begin + for transform in (adjoint, transpose) + usv = svd(transform(a)) + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ transform(a) + @test convert(Array, usv) ≈ transform(a) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ transform(a)\b + end end @testset "Generalized svd" begin - a_svd = a[1:n1, :] + a_svd = a[1:div(n, 2), :] gsvd = svd(a,a_svd) @test gsvd.U*gsvd.D1*gsvd.R*gsvd.Q' ≈ a @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ a_svd