Skip to content

Commit d7942fb

Browse files
committed
added stevengj's suggestion to qr. reversed changes in svd
1 parent 7764601 commit d7942fb

File tree

3 files changed

+3
-46
lines changed

3 files changed

+3
-46
lines changed

src/qr.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -538,10 +538,10 @@ end
538538
function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol))
539539
m = min(size(A)...)
540540
m == 0 && return 0
541-
rdiag = diag(getfield(A, :factors))
542-
tol = max(atol, rtol*abs(rdiag[1]))
541+
factors = getfield(A, :factors)
542+
tol = max(atol, rtol*abs(factors[1,1]))
543543

544-
return something(findfirst(abs.(rdiag) .<= tol), m+1) - 1
544+
return something(findfirst(i -> abs(factors[i,i]) <= tol, 1:m), m+1) - 1
545545
end
546546

547547
# Julia implementation similar to xgelsy

src/svd.jl

Lines changed: 0 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -245,37 +245,6 @@ svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))]
245245
svdvals(x::Number) = abs(x)
246246
svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T}
247247

248-
"""
249-
rank(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}
250-
251-
Compute the numerical rank of a `n × m` matrix by counting how many singular values are greater
252-
than `max(atol, rtol*σ₁)` where `σ₁` is `A`'s largest calculated singular value.
253-
`atol` and `rtol` are the absolute and relative tolerances, respectively.
254-
The default relative tolerance is `n*ϵ`,
255-
where `n` is the size of the smallest dimension of `A`,
256-
and `ϵ` is the [`eps`](@ref) of the element type of `A`.
257-
!!! note
258-
Numerical rank can be a sensitive and imprecise characterization of
259-
ill-conditioned matrices with singular values that are close to the threshold
260-
tolerance `max(atol, rtol*σ₁)`. In such cases, slight perturbations to the
261-
singular-value computation or to the matrix can change the result of `rank`
262-
by pushing one or more singular values across the threshold. These variations
263-
can even occur due to changes in floating-point errors between different Julia
264-
versions, architectures, compilers, or operating systems.
265-
266-
!!! compat "Julia 1.1"
267-
The `atol` and `rtol` keyword arguments requires at least Julia 1.1.
268-
In Julia 1.0 `rtol` is available as a positional argument, but this
269-
will be deprecated in Julia 2.0.
270-
271-
"""
272-
273-
function rank(S::SVD{<:Any,T}; atol::Real = 0.0, rtol::Real = (min(size(S)...)*eps(real(float(one(eltype(S))))))) where {T}
274-
svals = getfield(S, :S)
275-
tol = max(atol, rtol*svals[1])
276-
count(>(tol), svals)
277-
end
278-
279248
### SVD least squares ###
280249
function ldiv!(A::SVD{T}, B::AbstractVecOrMat) where T
281250
m, n = size(A)

test/svd.jl

Lines changed: 0 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -294,16 +294,4 @@ end
294294
@test F.S F32.S
295295
end
296296

297-
@testset "rank svd" begin
298-
# Test that the rank of an svd is computed correctly
299-
@test rank(svd([1.0 0.0; 0.0 1.0])) == 2
300-
@test rank(svd([1.0 0.0; 0.0 0.9]), rtol=0.95) == 1
301-
@test rank(svd([1.0 0.0; 0.0 0.9]), atol=0.95) == 1
302-
@test rank(svd([1.0 0.0; 0.0 1.0]), rtol=1.01) == 0
303-
@test rank(svd([1.0 0.0; 0.0 1.0]), atol=1.01) == 0
304-
305-
@test rank(svd([1.0 2.0; 2.0 4.0])) == 1
306-
@test rank(svd([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0])) == 2
307-
end
308-
309297
end # module TestSVD

0 commit comments

Comments
 (0)