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.

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

Add relative and absolute tolerance for rank. #29926

Merged
merged 16 commits into from
Dec 6, 2018
25 changes: 17 additions & 8 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -712,13 +712,13 @@ end
###########################################################################################

"""
rank(A[, tol::Real])
rank(A::AbstractArray, atol::Real, rtol::Real)
rank(A[, tol::Real]) = rank(A, rtol=tol)
sam0410 marked this conversation as resolved.
Show resolved Hide resolved

Compute the rank of a matrix by counting how many singular
values of `A` have magnitude greater than `tol*σ₁` where `σ₁` is
`A`'s largest singular values. By default, the value of `tol` is the smallest
dimension of `A` multiplied by the [`eps`](@ref)
of the [`eltype`](@ref) of `A`.
values of `A` have magnitude greater than `max(atol, rtol*σ₁)` where `σ₁` is
`A`'s largest singular value, `atol` and `rtol` are the absolute and relative
tolerance, respectively.
sam0410 marked this conversation as resolved.
Show resolved Hide resolved

# Examples
```jldoctest
Expand All @@ -731,14 +731,23 @@ julia> rank(diagm(0 => [1, 0, 2]))
julia> rank(diagm(0 => [1, 0.001, 2]), 0.1)
2

julia> rank(diagm(0 => [1, 0.001, 2]), 0.00001)
julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.1)
2

julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.00001)
3

julia> rank(diagm(0 => [1, 0.001, 2]), atol=1.5)
1
```
"""
function rank(A::AbstractMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
function rank(A::AbstractMatrix; atol::Real = 0.0, rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(atol))
isempty(A) && return 0 # 0-dimensional case
s = svdvals(A)
count(x -> x > tol*s[1], s)
tol = max(atol, rtol*s[1])
count(x -> x > tol, s)
end
rank(A::AbstractMatrix, tol::Real) = rank(A,rtol=tol) # TODO: deprecate tol in 2.0
sam0410 marked this conversation as resolved.
Show resolved Hide resolved
rank(x::Number) = x == 0 ? 0 : 1

"""
Expand Down
3 changes: 3 additions & 0 deletions stdlib/LinearAlgebra/test/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,9 @@ end

@test rank(fill(0, 0, 0)) == 0
@test rank([1.0 0.0; 0.0 0.9],0.95) == 1
@test rank([1.0 0.0; 0.0 0.9],rtol=0.95) == 1
@test rank([1.0 0.0; 0.0 0.9],atol=0.95) == 1
@test rank([1.0 0.0; 0.0 0.9],atol=0.95,rtol=0.95)==1
@test qr(big.([0 1; 0 0])).R == [0 1; 0 0]

@test norm([2.4e-322, 4.4e-323]) ≈ 2.47e-322
Expand Down