Add absolute and relative tolerances for rank

[lsorber]

Related julia-dev discussion: https://groups.google.com/d/msg/julia-dev/GAdcYzmibyo/iOLpyVQc8YIJ

The function rank should allow the user to set both absolute and relative tolerances using keyword arguments:

rank(A,reltol=1e-6)
rank(A,abstol=1e-6)

where reltol is a relative tolerance with respect to the largest singular value of A.

As a bonus, it would also be nice to have a mlrank (multilinear rank), which is one way to generalize rank to tensors (the only computationally tractable way). The multilinear rank of a tensor T is defined as a vector/tuple of ranks R, one for each mode of the tensor (in the matrix case, the column rank equals the row rank). R[n] is defined as the dimension of the space spanned by mode-n vectors of T. To get R[1], you take the svd of a matrix containing all column vectors of T and look at its dimensionality. I.e., R[1] = sum(svdvals(T[:,:]) .> abstol). For R[2], you do the same but for the row vectors of T, and so on.

[kshyatt]

This would be a relatively easy first PR if someone wants to tackle it.

[miguelraz]

@kshyatt Thanks a lot for the heads up - this is my very first issue correction attempt. Any comments are more than appreciated!

I have at least 2 problems with my submission so far:

1. I don't know how to easily parse keywords arguments into different cases.
2. I don't know how to make the 'else' case the default without throwing a useful error there.

Any comments for a noobie are more than appreciated.This is my attempt at modifying the code in linalg/generics.jl at about ~ line 559.

"""
    rank(M[; abstol::Real])

Compute the rank of a matrix by counting how many singular
values of `M` have magnitude greater than `abstol`.
By default, the value of `abstol` is the largest
dimension of `M` multiplied by the [`eps`](:func:`eps`)
of the [`eltype`](:func:`eltype`) of `M`, and `reltol` is the
relative tolerance.
"""

function rank(A::AbstractMatrix; kwargs...)
    if kwargs == abstol
        rank(A::AbstractMatrix; kwargs...) = sum(svdvals(A) .> kwargs)        
    elseif kwargs == reltol 
         rank(A::AbstractMatrix; kwargs...) = sum(div(svdvals(A),maximum(A)) .> kwargs)
    else
        throw(ArgumentError("Did not get a matrix or appropriate keyword arguments"))
    end
end

[jw3126]

I would try something like this:

function rank(A::AbstractMatrix; abstol=0., reltol=0.)
    thresh = abstol + reltol * maximum(svdvals(A))
    sum(svdvals(A) .> thresh )        
end

[miguelraz]

Thanks a lot @jw3126 ! I just tried to submit it as a first PR.
If this goes through I will try attacking the mlrank(A) problem.

The PR submission is JuliaLang/julia#19014 and I appreciate any and all criticism.
Thanks!

[miguelraz]

Stubborn Learning attempt JuliaLang/julia#3 this new PR.

[sam0410]

Hi @miguelraz, If you are not working on this, Can I take this up?

[StefanKarpinski]

Sure, no need to ask—just go for it and make a PR!

[miguelraz]

Go for it @sam0410 !

[simonbyrne]

Fixed by JuliaLang/julia#29926.