Add nullspace and cond functions for QRPivoted matrices

NEWS.md

@@ -87,7 +87,9 @@ Standard library changes
 
 #### LinearAlgebra
 
-* `rank` can now take a `QRPivoted` matrix to allow rank estimation via QR factorization ([#54283]).
+* `rank` can now take a `QRPivoted` matrix to allow rank estimation via QR factorization ([#54283])
+* `nullspace` can now take a `QRPivoted` matrix to allow kernel computation via QR factorization ([#54519]).
+* `cond` can now take a `QRPivoted` matrix to allow condition number estimation via QR factorization ([#54519]).
 
 #### Logging
 

stdlib/LinearAlgebra/src/qr.jl

@@ -542,6 +542,36 @@ function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(
     return something(findfirst(i -> abs(A.R[i,i]) <= tol, 1:m), m+1) - 1
 end
 
+function nullspace(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol))
+    m, n = size(A, 1), size(A, 2)
+    (m == 0 || n == 0) && return Matrix{eigtype(eltype(A.Q))}(I, n, n)
+
+    indstart = rank(A) + 1
+    return A.Q[:, indstart:end]
+end
+
+function cond(A::QRPivoted, p::Real=2)
+    m = min(size(A.R)...)
+    if p == 2
+        maxv = abs(A.R[1,1])
+        return iszero(maxv) ? oftype(real(maxv), Inf) : maxv / abs(A.R[m,m])
+    elseif p == 1 || p == Inf
+        checksquare(A.R)
+        R = UpperTriangular(A.R)
+        try
+            Rinv = inv(R)
+            return opnorm(R, p)*opnorm(Rinv, p)
+        catch e
+            if isa(e, LAPACKException) || isa(e, SingularException)
+                return convert(float(real(eltype(A))), Inf)
+            else
+                rethrow()
+            end
+        end
+    end
+    throw(ArgumentError(lazy"p-norm must be 1, 2 or Inf, got $p"))
+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)

stdlib/LinearAlgebra/test/qr.jl

@@ -214,6 +214,15 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = Matrix(Q)
                 @test_throws DimensionMismatch mul!(Vector{eltya}(undef, n+1), q, b)
             end
         end
+
+        @testset "Test nullspace (issue #54354)" begin
+            a15null = nullspace(qr(areal[:,1:n1], ColumnNorm()))
+            @test rank([areal[:,1:n1] a15null]) == 10
+            @test norm(areal[:,1:n1]'a15null,Inf) ≈ zero(eltya) atol=300ε
+            @test norm(a15null'areal[:,1:n1],Inf) ≈ zero(eltya) atol=400ε
+            @test size(nullspace(qr(a', ColumnNorm())), 2) == 0
+            @test nullspace(qr(zeros(eltya,n, n), ColumnNorm())) == Matrix(I, n, n)
+        end
     end
 end
 
@@ -540,4 +549,24 @@ end
     @test rank(qr([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0], ColumnNorm())) == 2
 end
 
+@testset "Test matrix condition number (issue #54354)" begin
+    ainit = rand(n, n)
+    @testset "for $elty" for elty in (Float32, Float64, ComplexF32, ComplexF64)
+        ainit = convert(Matrix{elty}, ainit)
+        for a in (copy(ainit), view(ainit, 1:n, 1:n))
+            F = qr(a, ColumnNorm())
+            Rinv = inv(F.R)
+            @test cond(F, 1)   == opnorm(F.R, 1) * opnorm(Rinv, 1)
+            @test cond(F, Inf) == opnorm(F.R, Inf) * opnorm(Rinv, Inf)
+            @test_throws ArgumentError cond(F,3)
+        end
+    end
+    @testset "Singular matrices" for p in (1, 2, Inf)
+        @test cond(qr(zeros(Int, 2, 2), ColumnNorm()), p) == Inf
+        @test cond(qr(zeros(2, 2), ColumnNorm()), p)      == Inf
+        @test cond(qr([0 0; 1 1], ColumnNorm()), p)       == Inf
+        @test cond(qr([0. 0.; 1. 1.], ColumnNorm()), p)   == Inf
+    end
+end
+
 end # module TestQR