return SVD from pinv(::SVD)

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

@@ -306,17 +306,18 @@ 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)
     # TODO: 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
 
 size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim)

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