diff --git a/NEWS.md b/NEWS.md index 168a2d03293d4..6a3ad9246d1a1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -87,6 +87,8 @@ Standard library changes #### LinearAlgebra +* `rank` can now take a `QRPivoted` matrix to allow rank estimation via QR factorization ([#54283]). + #### Logging #### Printf diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 06c2fba2932f5..0f81a07e12b08 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -535,6 +535,13 @@ function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T} return B 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 +end + # Julia implementation similar to xgelsy function ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractMatrix{T}, rcond::Real) where {T<:BlasFloat} require_one_based_indexing(B) diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index e339706598a8a..aef2e23740487 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -528,4 +528,16 @@ end end end +@testset "issue #53214" begin + # Test that the rank of a QRPivoted matrix is computed correctly + @test rank(qr([1.0 0.0; 0.0 1.0], ColumnNorm())) == 2 + @test rank(qr([1.0 0.0; 0.0 0.9], ColumnNorm()), rtol=0.95) == 1 + @test rank(qr([1.0 0.0; 0.0 0.9], ColumnNorm()), atol=0.95) == 1 + @test rank(qr([1.0 0.0; 0.0 1.0], ColumnNorm()), rtol=1.01) == 0 + @test rank(qr([1.0 0.0; 0.0 1.0], ColumnNorm()), atol=1.01) == 0 + + @test rank(qr([1.0 2.0; 2.0 4.0], ColumnNorm())) == 1 + @test rank(qr([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0], ColumnNorm())) == 2 +end + end # module TestQR