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
3 changes: 2 additions & 1 deletion stdlib/LinearAlgebra/src/deprecated.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base: @deprecate, depwarn
# To be deprecated in 2.0
rank(A::AbstractMatrix, tol::Real) = rank(A,rtol=tol)
25 changes: 16 additions & 9 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -712,13 +712,15 @@ end
###########################################################################################

"""
rank(A[, tol::Real])
rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
rank(A::AbstractMatrix, rtol::Real) = rank(A; rtol=rtol) # to be deprecated in Julia 2.0

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
tolerances, respectively. The default relative tolerance is `n*ϵ`, where `n`
is the size of the smallest dimension of `A`, and `ϵ` is the [`eps`](@ref) of
the element type of `A`.

# Examples
```jldoctest
Expand All @@ -728,16 +730,21 @@ julia> rank(Matrix(I, 3, 3))
julia> rank(diagm(0 => [1, 0, 2]))
2

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

julia> rank(diagm(0 => [1, 0.001, 2]), 0.00001)
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(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