Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding rtol and atol for pinv and nullspace #29998

Merged
merged 12 commits into from
Dec 11, 2018
Merged

Adding rtol and atol for pinv and nullspace #29998

Changes from 1 commit
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
759afc8
Add code for rtol and atol
Nov 10, 2018
62c9adf
add tests
Nov 11, 2018
278321e
Merge branch 'master' of https://github.com/JuliaLang/julia into pinv…
Nov 20, 2018
82afdf2
resolve comment
Nov 21, 2018
46d3cd6
fix typo
Nov 22, 2018
c3d8fc6
fix tests
Nov 23, 2018
20346b5
Merge branch 'master' of https://github.com/JuliaLang/julia into pinv…
Dec 7, 2018
3e28425
add news.md item
Dec 7, 2018
6dbd52a
Not deprecated yet.
simonbyrne Dec 7, 2018
d90190d
Tweak docs slightly
simonbyrne Dec 7, 2018
6418489
Merge branch 'master' into pinvNull
sam0410 Dec 10, 2018
367c63b
typo from diff
simonbyrne Dec 11, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Add code for rtol and atol
  • Loading branch information
sam0410 committed Nov 10, 2018
commit 759afc8b6e6fc39716062d9755f32aaa97f4739e
48 changes: 33 additions & 15 deletions stdlib/LinearAlgebra/src/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1201,13 +1201,14 @@ factorize(A::Transpose) = transpose(factorize(parent(A)))
## Moore-Penrose pseudoinverse

"""
pinv(M[, rtol::Real])
pinv(M, atol::Real, rtol::Real)
sam0410 marked this conversation as resolved.
Show resolved Hide resolved
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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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))
Expand All @@ -1285,11 +1285,12 @@ end
## Basis for null space

"""
nullspace(M[, rtol::Real])
nullspace(M, atol::Real, rtol::Real)
sam0410 marked this conversation as resolved.
Show resolved Hide resolved
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`.

Expand All @@ -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)
Expand Down