From e52ca47b820cd391384fd983b622cef2d80a76cc Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 2 Jul 2025 10:18:27 -0400 Subject: [PATCH 1/9] return SVD from pinv(::SVD) --- src/dense.jl | 3 ++- src/svd.jl | 7 ++++--- test/svd.jl | 7 ++++--- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index 7f996ec8..f6f1accb 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -1756,7 +1756,8 @@ SVD-based algorithm, it is better to employ the SVD directly via `svd(M; rtol, a or `ldiv!(svd(M), b; rtol, atol)`. One can also pass `M = svd(A)` as the argument to `pinv` in order to re-use -an existing [`SVD`](@ref) factorization. +an existing [`SVD`](@ref) factorization. In this case, `pinv` will return +the SVD of the pseudo-inverse, which can be applied accurately, instead of an explicit matrix. !!! compat "Julia 1.13" Passing an `SVD` object to `pinv` requires Julia 1.13 or later. diff --git a/src/svd.jl b/src/svd.jl index 00ec8793..809cddbe 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -306,9 +306,9 @@ function ldiv!(F::SVD{T}, B::AbstractVecOrMat; atol::Real=0, rtol::Real = (eps(r return B end -function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where T +function pinv(F::SVD{T,<:Any,M}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where {T,M} k = _count_svdvals(F.S, atol, rtol) - @views (F.S[1:k] .\ F.Vt[1:k, :])' * F.U[:,1:k]' + @views SVD(M(F.Vt[1:k, :]'), inv.(F.S[1:k]), M(F.U[:,1:k]')) end function inv(F::SVD) @@ -316,7 +316,8 @@ function inv(F::SVD) @inbounds for i in eachindex(F.S) iszero(F.S[i]) && throw(SingularException(i)) end - pinv(F; rtol=eps(real(eltype(F)))) + k = _count_svdvals(F.S, 0, eps(real(eltype(F)))) + return @views (F.S[1:k] .\ F.Vt[1:k, :])' * F.U[:,1:k]' end size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim) diff --git a/test/svd.jl b/test/svd.jl index fb4e416e..62a5c136 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -239,12 +239,13 @@ end @testset "SVD pinv and truncation" begin m, n = 10,5 A = randn(m,n) * [1/(i+j-1) for i = 1:n, j=1:n] # badly conditioned Hilbert matrix - @test pinv(A) ≈ pinv(svd(A)) rtol=1e-13 + @test pinv(A) ≈ Matrix(pinv(svd(A))) rtol=1e-13 pinv_3 = pinv(A, rtol=1e-3) - @test pinv_3 ≈ pinv(svd(A), rtol=1e-3) rtol=1e-13 - @test pinv_3 ≈ pinv(svd(A, rtol=1e-3)) rtol=1e-13 + @test pinv_3 ≈ Matrix(pinv(svd(A), rtol=1e-3)) rtol=1e-13 + @test pinv_3 ≈ Matrix(pinv(svd(A, rtol=1e-3))) rtol=1e-13 b = float([1:m;]) # arbitrary rhs @test pinv_3 * b ≈ svd(A, rtol=1e-3) \ b rtol=1e-13 + @test pinv_3 * b ≈ pinv(svd(A, rtol=1e-3)) * b rtol=1e-13 @test pinv_3 * b ≈ ldiv!(svd(A), copy(b), rtol=1e-3)[1:n] rtol=1e-13 @test pinv(A, atol=100) == pinv(svd(A), atol=100) == pinv(svd(A, atol=100)) == zeros(5,10) end From e9328ebf2ec407a82eb58b5079b5dc782bdf62a2 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 2 Jul 2025 15:12:07 -0400 Subject: [PATCH 2/9] Update svd.jl --- test/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/svd.jl b/test/svd.jl index 62a5c136..bdcc9e9c 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -247,7 +247,7 @@ end @test pinv_3 * b ≈ svd(A, rtol=1e-3) \ b rtol=1e-13 @test pinv_3 * b ≈ pinv(svd(A, rtol=1e-3)) * b rtol=1e-13 @test pinv_3 * b ≈ ldiv!(svd(A), copy(b), rtol=1e-3)[1:n] rtol=1e-13 - @test pinv(A, atol=100) == pinv(svd(A), atol=100) == pinv(svd(A, atol=100)) == zeros(5,10) + @test pinv(A, atol=100) == Matrix(pinv(svd(A), atol=100)) == Matrix(pinv(svd(A, atol=100))) == zeros(5,10) end @testset "Issue 40944. ldiv!(SVD) should update rhs" begin From 7549fea886b2d3e101db4c07a18c2c99b257331e Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 2 Jul 2025 21:20:38 -0400 Subject: [PATCH 3/9] add SVD * X --- src/svd.jl | 4 ++++ test/svd.jl | 18 +++++++++++------- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/svd.jl b/src/svd.jl index 809cddbe..c61c49c2 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -320,6 +320,10 @@ function inv(F::SVD) return @views (F.S[1:k] .\ F.Vt[1:k, :])' * F.U[:,1:k]' end +# multiplying SVD by matrix/vector, mainly useful for pinv(::SVD) output +(*)(F::SVD, A::AbstractVecOrMat{<:Number}) = F.U * (Diagonal(F.S) * (F.Vt * A)) +(*)(A::AbstractMatrix{<:Number}, F::SVD) = ((A*F.U) * Diagonal(F.S)) * F.Vt + size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim) size(A::SVD) = (size(A, 1), size(A, 2)) diff --git a/test/svd.jl b/test/svd.jl index bdcc9e9c..8afbbb6c 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -239,15 +239,19 @@ end @testset "SVD pinv and truncation" begin m, n = 10,5 A = randn(m,n) * [1/(i+j-1) for i = 1:n, j=1:n] # badly conditioned Hilbert matrix - @test pinv(A) ≈ Matrix(pinv(svd(A))) rtol=1e-13 + F = svd(A) + @test pinv(A) ≈ Matrix(pinv(F)) rtol=1e-13 pinv_3 = pinv(A, rtol=1e-3) - @test pinv_3 ≈ Matrix(pinv(svd(A), rtol=1e-3)) rtol=1e-13 - @test pinv_3 ≈ Matrix(pinv(svd(A, rtol=1e-3))) rtol=1e-13 + F_3 = svd(A, rtol=1e-3) + @test pinv_3 ≈ Matrix(pinv(F, rtol=1e-3)) rtol=1e-13 + @test pinv_3 ≈ Matrix(pinv(F_3)) rtol=1e-13 b = float([1:m;]) # arbitrary rhs - @test pinv_3 * b ≈ svd(A, rtol=1e-3) \ b rtol=1e-13 - @test pinv_3 * b ≈ pinv(svd(A, rtol=1e-3)) * b rtol=1e-13 - @test pinv_3 * b ≈ ldiv!(svd(A), copy(b), rtol=1e-3)[1:n] rtol=1e-13 - @test pinv(A, atol=100) == Matrix(pinv(svd(A), atol=100)) == Matrix(pinv(svd(A, atol=100))) == zeros(5,10) + @test pinv_3 * b ≈ F_3 \ b rtol=1e-13 + @test pinv_3 * b ≈ pinv(F_3) * b rtol=1e-13 + @test pinv_3 * b ≈ ldiv!(F, copy(b), rtol=1e-3)[1:n] rtol=1e-13 + c = float([1:n;]) # arbitrary rhs + @test c' * pinv_3 ≈ c' * pinv(F_3) rtol=1e-13 + @test pinv(A, atol=100) == Matrix(pinv(F, atol=100)) == Matrix(pinv(svd(A, atol=100))) == zeros(5,10) end @testset "Issue 40944. ldiv!(SVD) should update rhs" begin From 5cf400b95f081ba2ead9efaf1e97e95504a6bc74 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 3 Jul 2025 22:21:56 -0400 Subject: [PATCH 4/9] copy gives the right type --- src/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svd.jl b/src/svd.jl index c61c49c2..783aafce 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -308,7 +308,7 @@ end function pinv(F::SVD{T,<:Any,M}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where {T,M} k = _count_svdvals(F.S, atol, rtol) - @views SVD(M(F.Vt[1:k, :]'), inv.(F.S[1:k]), M(F.U[:,1:k]')) + @views SVD(copy(F.Vt[1:k, :]'), inv.(F.S[1:k]), copy(F.U[:,1:k]')) end function inv(F::SVD) From cecdb700fe993ea262a23bc1d22b8eb0620fcd73 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 3 Jul 2025 22:23:33 -0400 Subject: [PATCH 5/9] rm used type param --- src/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svd.jl b/src/svd.jl index 783aafce..ed42befb 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -306,7 +306,7 @@ function ldiv!(F::SVD{T}, B::AbstractVecOrMat; atol::Real=0, rtol::Real = (eps(r return B end -function pinv(F::SVD{T,<:Any,M}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where {T,M} +function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where {T} k = _count_svdvals(F.S, atol, rtol) @views SVD(copy(F.Vt[1:k, :]'), inv.(F.S[1:k]), copy(F.U[:,1:k]')) end From 8adb4560281c83301903b67e14d44435903e5718 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 3 Jul 2025 22:24:07 -0400 Subject: [PATCH 6/9] rm extraneous brackets --- src/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svd.jl b/src/svd.jl index ed42befb..85e3f894 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -306,7 +306,7 @@ function ldiv!(F::SVD{T}, B::AbstractVecOrMat; atol::Real=0, rtol::Real = (eps(r return B end -function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where {T} +function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where T k = _count_svdvals(F.S, atol, rtol) @views SVD(copy(F.Vt[1:k, :]'), inv.(F.S[1:k]), copy(F.U[:,1:k]')) end From 8626a91289ac7f612fd0085f9fce73578692bb56 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 3 Jul 2025 22:25:26 -0400 Subject: [PATCH 7/9] keep svdvals sorted --- src/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svd.jl b/src/svd.jl index 85e3f894..70fe20b9 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -308,7 +308,7 @@ end function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where T k = _count_svdvals(F.S, atol, rtol) - @views SVD(copy(F.Vt[1:k, :]'), inv.(F.S[1:k]), copy(F.U[:,1:k]')) + @views SVD(copy(F.Vt[1:k, :]'), reverse!(inv.(F.S[1:k])), copy(F.U[:,1:k]')) end function inv(F::SVD) From 52b2d58e0c76a39cb2fec5da6405b270872fad8b Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 3 Jul 2025 22:26:17 -0400 Subject: [PATCH 8/9] whoops --- src/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svd.jl b/src/svd.jl index 70fe20b9..04cb6e03 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -308,7 +308,7 @@ end function pinv(F::SVD{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(F)...))*iszero(atol)) where T k = _count_svdvals(F.S, atol, rtol) - @views SVD(copy(F.Vt[1:k, :]'), reverse!(inv.(F.S[1:k])), copy(F.U[:,1:k]')) + @views SVD(copy(F.Vt[k:-1:1, :]'), inv.(F.S[k:-1:1]), copy(F.U[:,k:-1:1]')) end function inv(F::SVD) From b3921af31e74233befab2c8160d0f56794cbceff Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 4 Jul 2025 09:12:17 -0400 Subject: [PATCH 9/9] fix test --- test/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/svd.jl b/test/svd.jl index 8afbbb6c..e9d039ea 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -56,7 +56,7 @@ using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted @test_throws DimensionMismatch inv(svd(Matrix(I, 3, 2))) @test inv(svd(Matrix(I, 2, 2))) ≈ I @test inv(svd([1 2; 3 4])) ≈ [-2.0 1.0; 1.5 -0.5] - @test pinv(svd([1 0 1; 0 1 0])) ≈ [0.5 0.0; 0.0 1.0; 0.5 0.0] + @test Matrix(pinv(svd([1 0 1; 0 1 0]))) ≈ [0.5 0.0; 0.0 1.0; 0.5 0.0] @test_throws SingularException inv(svd([0 0; 0 0])) @test inv(svd([1+2im 3+4im; 5+6im 7+8im])) ≈ [-0.5 + 0.4375im 0.25 - 0.1875im; 0.375 - 0.3125im -0.125 + 0.0625im] end