Skip to content

Commit

Permalink
Adding rtol and atol for pinv and nullspace (#29998)
Browse files Browse the repository at this point in the history 
* Add code for rtol and atol

* add tests

* resolve comment

* fix typo

* fix tests

* add news.md item

* Not deprecated yet.

* Tweak docs slightly

* typo from diff

[skip ci]
  • Loading branch information
@sam0410 @simonbyrne
sam0410 authored and simonbyrne committed Dec 11, 2018
1 parent 217d330 commit 5b2e3e7
Show file tree
Hide file tree
Showing 4 changed files with 46 additions and 23 deletions.
1 change: 1 addition & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ Standard library changes
* `mul!`, `rmul!` and `lmul!` methods for `UniformScaling` ([#29506]).
* `Symmetric` and `Hermitian` matrices now preserve the wrapper when scaled with a number ([#29469]).
* Exponentiation operator `^` now supports raising an `Irrational` to an `AbstractMatrix` power ([#29782]).
* Added keyword arguments `rtol`, `atol` to `pinv`, `nullspace` and `rank` ([#29998], [#29926]).

#### Random
* `randperm` and `randcycle` now use the type of their argument to determine the element type of
Expand Down
57 changes: 34 additions & 23 deletions stdlib/LinearAlgebra/src/dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1237,19 +1237,21 @@ factorize(A::Transpose) = transpose(factorize(parent(A)))
## Moore-Penrose pseudoinverse

"""
pinv(M[, rtol::Real])
pinv(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
pinv(M, rtol::Real) = pinv(M; rtol=rtol) # to be deprecated in Julia 2.0
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*σ₁)` where `σ₁` is the largest singular value of `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
`eps(real(float(one(eltype(M)))))*minimum(size(M))`, which is essentially machine epsilon
for the real part of a matrix element multiplied by the larger matrix dimension. For
inverting dense ill-conditioned matrices in a least-squares sense,
The optimal choice of absolute (`atol`) and relative tolerance (`rtol`) varies
both with the value of `M` and the intended application of the pseudoinverse.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest
dimension of `M`, and `ϵ` is the [`eps`](@ref) of the element type of `M`.
For inverting dense ill-conditioned matrices in a least-squares sense,
`rtol = sqrt(eps(real(float(one(eltype(M))))))` is recommended.
For more information, see [^issue8859], [^B96], [^S84], [^KY88].
Expand Down Expand Up @@ -1280,7 +1282,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
Expand All @@ -1289,9 +1291,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
Expand All @@ -1302,17 +1305,14 @@ 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
function pinv(x::Number)
xi = inv(x)
return ifelse(isfinite(xi), xi, zero(xi))
Expand All @@ -1321,13 +1321,16 @@ end
## Basis for null space

"""
nullspace(M[, rtol::Real])
nullspace(M; atol::Real=0, rtol::Rea=atol>0 ? 0 : n*ϵ)
nullspace(M, rtol::Real) = nullspace(M; rtol=rtol) # to be deprecated in Julia 2.0
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
`rtol` is the smallest dimension of `A` multiplied by the [`eps`](@ref)
of the [`eltype`](@ref) of `A`.
vectors of A whose singular have magnitude are greater than `max(atol, rtol*σ₁)`,
where `σ₁` is `M`'s largest singularvalue.
By default, the relative tolerance `rtol` is `n*ϵ`, where `n`
is the size of the smallest dimension of `M`, and `ϵ` is the [`eps`](@ref) of
the element type of `M`.
# Examples
```jldoctest
Expand All @@ -1343,21 +1346,29 @@ julia> nullspace(M)
0.0
1.0
julia> nullspace(M, 2)
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::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)

"""
cond(M, p::Real=2)
Expand Down
3 changes: 3 additions & 0 deletions stdlib/LinearAlgebra/src/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,6 @@

# To be deprecated in 2.0
rank(A::AbstractMatrix, tol::Real) = rank(A,rtol=tol)
nullspace(A::AbstractVector, tol::Real) = nullspace(reshape(A, length(A), 1), rtol= tol)
nullspace(A::AbstractMatrix, tol::Real) = nullspace(A, rtol=tol)
pinv(A::AbstractMatrix{T}, tol::Real) where T = pinv(A, rtol=tol)
8 changes: 8 additions & 0 deletions stdlib/LinearAlgebra/test/dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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=0.001), 2) == 0
@test size(nullspace(b, atol=100*εb), 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)
Expand All @@ -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])
Expand Down

0 comments on commit 5b2e3e7

Please sign in to comment.