faster implementation of rank(::QRPivoted) fixes #1128

src/qr.jl

@@ -538,8 +538,10 @@ end
 function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol))
     m = min(size(A)...)
     m == 0 && return 0
-    tol = max(atol, rtol*abs(A.R[1,1]))
-    return something(findfirst(i -> abs(A.R[i,i]) <= tol, 1:m), m+1) - 1
+    rdiag = diag(getfield(A, :factors))
+    tol = max(atol, rtol*abs(rdiag[1]))
+
+    return something(findfirst(abs.(rdiag) .<= tol), m+1) - 1
 end
 
 # Julia implementation similar to xgelsy

src/svd.jl

@@ -245,6 +245,37 @@ svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))]
 svdvals(x::Number) = abs(x)
 svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T}
 
+"""
+    rank(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}
+
+Compute the numerical rank of a `n × m` matrix by counting how many singular values are greater 
+than `max(atol, rtol*σ₁)` where `σ₁` is `A`'s largest calculated 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`.
+!!! note
+    Numerical rank can be a sensitive and imprecise characterization of
+    ill-conditioned matrices with singular values that are close to the threshold
+    tolerance `max(atol, rtol*σ₁)`. In such cases, slight perturbations to the
+    singular-value computation or to the matrix can change the result of `rank`
+    by pushing one or more singular values across the threshold. These variations
+    can even occur due to changes in floating-point errors between different Julia
+    versions, architectures, compilers, or operating systems.
+
+!!! compat "Julia 1.1"
+    The `atol` and `rtol` keyword arguments requires at least Julia 1.1.
+    In Julia 1.0 `rtol` is available as a positional argument, but this
+    will be deprecated in Julia 2.0.
+
+"""
+
+function rank(S::SVD{<:Any,T}; atol::Real = 0.0, rtol::Real = (min(size(S)...)*eps(real(float(one(eltype(S))))))) where {T}
+    svals = getfield(S, :S)
+    tol = max(atol, rtol*svals[1])
+    count(>(tol), svals)
+end
+
 ### SVD least squares ###
 function ldiv!(A::SVD{T}, B::AbstractVecOrMat) where T
     m, n = size(A)

test/svd.jl

@@ -294,4 +294,16 @@ end
     @test F.S ≈ F32.S
 end
 
+@testset "rank svd" begin
+    # Test that the rank of an svd is computed correctly
+    @test rank(svd([1.0 0.0; 0.0 1.0])) == 2
+    @test rank(svd([1.0 0.0; 0.0 0.9]), rtol=0.95) == 1
+    @test rank(svd([1.0 0.0; 0.0 0.9]), atol=0.95) == 1
+    @test rank(svd([1.0 0.0; 0.0 1.0]), rtol=1.01) == 0
+    @test rank(svd([1.0 0.0; 0.0 1.0]), atol=1.01) == 0
+
+    @test rank(svd([1.0 2.0; 2.0 4.0])) == 1
+    @test rank(svd([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0])) == 2
+end
+
 end # module TestSVD