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
@@ -715,10 +715,9 @@ end
rank(A[, tol::Real])
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 `rtol*σ₁ + atol` where `σ₁` is
`A`'s largest singular values, atol and rtol are the absolute and relative
This conversation was marked as resolved by sam0410

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 10, 2018

Member

largest singular value

Also, put backticks around atol and rtol

tolerance respectively.
This conversation was marked as resolved by sam0410

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 10, 2018

Member

comma before "respectively"

# Examples
```jldoctest
@@ -728,16 +727,19 @@ 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=0.00001, rtol=0.00001)
This conversation was marked as resolved by sam0410

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 10, 2018

Member

Instead of this, may be better to have an example that actually uses the atol, e.g.

julia> rank(diagm(0 => [1, 0.001, 2]), atol=1.5)
1
3
```
"""
function rank(A::AbstractMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
function rank(A::AbstractMatrix; atol=0.0, rtol=0.0)
This conversation was marked as resolved by sam0410

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 8, 2018

Member

This seems like a breaking change? We need to keep a rank(A, tol) = rank(A, rtol=tol) method in Julia 1.x.

Why is there no nonzero default rtol? (Also, similar to #22742, we should probably set the default rtol to zero if atol > 0 was specified.)

This comment has been minimized.

Copy link
@StefanKarpinski

StefanKarpinski Nov 8, 2018

Member

Good point; leave a TODO comment about deprecating tol in Julia 2.0, however.

s = svdvals(A)
count(x -> x > tol*s[1], s)
count(x -> x > (atol + rtol * s[1]), s)
This conversation was marked as resolved by sam0410

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 8, 2018

Member

This should be made consistent with the current behavior of isapprox, which was changed in #22742 to use max: test x > max(atol, rtol*s[1]).

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 8, 2018

Member

(I don't know if the compiler can hoist a max(atol, rtol*s[1]) call out of the count loop, but to be safe we should probably do it manually? That is, set tol = max(atol, rtol*s[1]) and test x > tol.)

This comment has been minimized.

Copy link
@sam0410

sam0410 Nov 8, 2018

Author Contributor

Hi @stevengj , if I add tol = max(atol, rtol*s[1]), this test case throws an error stating ERROR: BoundsError: attempt to access 0-element Array{Float64,1} at index [1]. I feel, we can remove this test-case. What are your opinions on this?

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 9, 2018

Member

A 0x0 matrix is well defined, and its rank should be zero (since its column space is 0-dimensional). Similarly for any 0xN or Mx0 matrix.

So, I would add a line

isempty(A) && return 0 # 0-dimensional case

at the top of the function (before the call to svdvals).

This comment has been minimized.

Copy link
@stevengj

stevengj Nov 9, 2018

Member

Or you could just leave it as count(x -> x > max(atol, rtol*s[1]), s). It's not like the performance matters for the count loop, since in any practical case it will be swamped by the cost of svdvals.

end
rank(x::Number) = x == 0 ? 0 : 1

@@ -194,7 +194,9 @@ end
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)==0
@test qr(big.([0 1; 0 0])).R == [0 1; 0 0]

@test norm([2.4e-322, 4.4e-323]) ≈ 2.47e-322
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.