return SVD from pinv(::SVD) (#1398)

As suggested by @andreasnoack in
#1387 (comment),
this changes the `pinv(::SVD)` function to instead return the SVD of the
pseudo-inverse, rather than an explicit matrix, so that it can be
applied stably.

(No backwards-compatibility issue since the `pinv(::SVD)` method was
introduced in #1387.)

Author: stevengj

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.

src/svd.jl

@@ -308,17 +308,22 @@ 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 (F.S[1:k] .\ F.Vt[1:k, :])' * 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)
     checksquare(F)
     @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
 
+# 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))
 

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
@@ -239,14 +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) ≈ 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 ≈ pinv(svd(A), rtol=1e-3)                        rtol=1e-13
-    @test pinv_3 ≈ 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 ≈ 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_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