Adding rtol and atol for pinv and nullspace

stdlib/LinearAlgebra/src/dense.jl

@@ -1201,13 +1201,14 @@ factorize(A::Transpose) = transpose(factorize(parent(A)))
 ## Moore-Penrose pseudoinverse
 
 """
-    pinv(M[, rtol::Real])
+    pinv(M, atol::Real, rtol::Real)
+    pinv(M[, tol::Real]) = pinv(M, rtol=tol)
 
 Computes the Moore-Penrose pseudoinverse.
 
 For matrices `M` with floating point elements, it is convenient to compute
 the pseudoinverse by inverting only singular values greater than
-`rtol * maximum(svdvals(M))`.
+`max(atol, rtol * maximum(svdvals(M)))`.
 
 The optimal choice of `rtol` varies both with the value of `M` and the intended application
 of the pseudoinverse. The default value of `rtol` is
@@ -1244,7 +1245,7 @@ julia> M * N
 
 [^KY88]: Konstantinos Konstantinides and Kung Yao, "Statistical analysis of effective singular values in matrix rank determination", IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. [doi:10.1109/29.1585](https://doi.org/10.1109/29.1585)
 """
-function pinv(A::AbstractMatrix{T}, rtol::Real) where T
+function pinv(A::AbstractMatrix{T}; atol::Real = 0.0, rtol::Real = (eps(real(float(one(T))))*min(size(A)...))*iszero(atol)) where T
     m, n = size(A)
     Tout = typeof(zero(T)/sqrt(one(T) + one(T)))
     if m == 0 || n == 0
@@ -1253,9 +1254,10 @@ function pinv(A::AbstractMatrix{T}, rtol::Real) where T
     if istril(A)
         if istriu(A)
             maxabsA = maximum(abs.(diag(A)))
+            tol = max(rtol*maxabsA, atol)
             B = zeros(Tout, n, m)
             for i = 1:min(m, n)
-                if abs(A[i,i]) > rtol*maxabsA
+                if abs(A[i,i]) > tol
                     Aii = inv(A[i,i])
                     if isfinite(Aii)
                         B[i,i] = Aii
@@ -1266,17 +1268,15 @@ function pinv(A::AbstractMatrix{T}, rtol::Real) where T
         end
     end
     SVD         = svd(A, full = false)
+    tol         = max(rtol*maximum(SVD.S), atol)
     Stype       = eltype(SVD.S)
     Sinv        = zeros(Stype, length(SVD.S))
-    index       = SVD.S .> rtol*maximum(SVD.S)
+    index       = SVD.S .> tol
     Sinv[index] = one(Stype) ./ SVD.S[index]
     Sinv[findall(.!isfinite.(Sinv))] .= zero(Stype)
     return SVD.Vt' * (Diagonal(Sinv) * SVD.U')
 end
-function pinv(A::AbstractMatrix{T}) where T
-    rtol = eps(real(float(one(T))))*min(size(A)...)
-    return pinv(A, rtol)
-end
+pinv(A::AbstractMatrix{T}, tol::Real) where T = nullspace(A, rtol=tol) # TODO: deprecate tol in 2.0
 function pinv(x::Number)
     xi = inv(x)
     return ifelse(isfinite(xi), xi, zero(xi))
@@ -1285,11 +1285,12 @@ end
 ## Basis for null space
 
 """
-    nullspace(M[, rtol::Real])
+    nullspace(M, atol::Real, rtol::Real)
+    nullspace(M[, tol::Real]) = nullspace(M, rtol=tol)
 
 Computes a basis for the nullspace of `M` by including the singular
-vectors of A whose singular have magnitude are greater than `rtol*σ₁`,
-where `σ₁` is `A`'s largest singular values. By default, the value of
+vectors of A whose singular have magnitude are greater than `max(atol, rtol*σ₁)`,
+where `σ₁` is `A`'s largest singular value. By default, the value of
 `rtol` is the smallest dimension of `A` multiplied by the [`eps`](@ref)
 of the [`eltype`](@ref) of `A`.
 
@@ -1312,16 +1313,33 @@ julia> nullspace(M, 2)
  0.0  1.0  0.0
  1.0  0.0  0.0
  0.0  0.0  1.0
+
+julia> nullspace(M, rtol=3)
+3×3 Array{Float64,2}:
+ 0.0  1.0  0.0
+ 1.0  0.0  0.0
+ 0.0  0.0  1.0
+
+julia> nullspace(M, atol=0.95)
+3×1 Array{Float64,2}:
+ 0.0
+ 0.0
+ 1.0
 ```
 """
-function nullspace(A::AbstractMatrix, rtol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
+function nullspace(A::AbstractMatrix; atol::Real = 0.0, rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(atol))
     m, n = size(A)
     (m == 0 || n == 0) && return Matrix{eltype(A)}(I, n, n)
     SVD = svd(A, full=true)
-    indstart = sum(s -> s .> SVD.S[1]*rtol, SVD.S) + 1
+    tol = max(atol, SVD.S[1]*rtol)
+    indstart = sum(s -> s .> tol, SVD.S) + 1
     return copy(SVD.Vt[indstart:end,:]')
 end
-nullspace(a::AbstractVector, rtol::Real = min(size(a)...)*eps(real(float(one(eltype(a)))))) = nullspace(reshape(a, length(a), 1), rtol)
+nullspace(A::AbstractMatrix, tol::Real) = nullspace(A, rtol=tol) # TODO: deprecate tol in 2.0
+
+nullspace(A::AbstractVector; atol::Real = 0.0, rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(atol)) = nullspace(reshape(a, length(a), 1), rtol= rtol, atol= atol)
+
+nullspace(A::AbstractVector, tol::Real) = nullspace(reshape(a, length(a), 1), rtol= tol) # TODO: deprecate tol in 2.0
 
 """
     cond(M, p::Real=2)

stdlib/LinearAlgebra/test/dense.jl

@@ -72,6 +72,8 @@ bimg  = randn(n,2)/2
                 @test norm(a[:,1:n1]'a15null,Inf) ≈ zero(eltya) atol=300ε
                 @test norm(a15null'a[:,1:n1],Inf) ≈ zero(eltya) atol=400ε
                 @test size(nullspace(b), 2) == 0
+                @test size(nullspace(b, rtol=2), 2) != 0
+                @test size(nullspace(b, atol=2), 2) == 0
                 @test size(nullspace(b, 100*εb), 2) == 0
                 @test nullspace(zeros(eltya,n)) == Matrix(I, 1, 1)
                 @test nullspace(zeros(eltya,n), 0.1) == Matrix(I, 1, 1)
@@ -82,6 +84,12 @@ bimg  = randn(n,2)/2
         end
     end # for eltyb
 
+@testset "Test pinv (rtol, atol)" begin
+    M = [1 0 0; 0 1 0; 0 0 0]
+    @test pinv(M,atol=1)== zeros(3,3)
+    @test pinv(M,rtol=0.5)== M
+end
+
     for (a, a2) in ((copy(ainit), copy(ainit2)), (view(ainit, 1:n, 1:n), view(ainit2, 1:n, 1:n)))
         @testset "Test pinv" begin
             pinva15 = pinv(a[:,1:n1])