Add nullspace and cond functions for QRPivoted matrices

[danilonumeroso]

This PR closes #54354 and also closes #40970.

[dkarrasch]

Since both Q and the permutation don't change norms, isn't this simply

Suggested change

	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
	cond(Q::QRPivoted, p::Real=2) = cond(Q.R, p)

? If we add this, we could also add the corresponding method for Q::QR.

[dkarrasch]

Indeed: #53214 (comment).

[stevengj]

Q doesn't change 2-norms, but it changes 1 norms and infinity norms.

[dkarrasch]

Right. So, together with your comment below, I agree that we should forward this to cond(Q.R, 2) in the p=2 case and otherwise throw and recommend using the original matrix, or, cond(Matrix(Q), p).

[stevengj]

If I recall correctly, we don't currently have efficient methods for slicing Q matrices (which are internally indirectly represented by linear operators). See https://discourse.julialang.org/t/what-is-the-most-efficient-way-of-obtaining-the-orthogonal-column-space-of-a-matrix/72212/13

So you should probably need to do something like:

d = n + 1 - indstart # dimension of nullspace
I_slice = zeros(eltype(A), m, d) # construct last d columns of I
for i = 1:d; I_slice[i, indstart + d - 1] = 1; end
return A.Q * I_slice

[dkarrasch]

We do:

julia/stdlib/LinearAlgebra/src/abstractq.jl

Lines 90 to 96 in 98f3947

	@inline function _getindex(Q::AbstractQ, ::Colon, J::AbstractVector{<:Integer})
	Y = zeros(eltype(Q), size(Q, 2), length(J))
	@inbounds for (i,j) in enumerate(J)
	Y[j,i] = oneunit(eltype(Q))
	end
	lmul!(Q, Y)
	end

[stevengj]

Suggested change

	return iszero(maxv) ? oftype(real(maxv), Inf) : maxv / abs(A.R[m,m])
	return iszero(maxv) ? oftype(maxv, Inf) : maxv / abs(A.R[m,m])

real seems redundant here.

[stevengj]

This algorithm seems incorrect. If $A = QR$, then cond(A) == cond(R), but cond(R) is not the ratio of diagonal entries of R.

[stevengj]

This algorithm also seems wrong?

I'm not there that there is any good way to compute the p=1 or Inf condition numbers from the QR factorization? Might better to just throw an error in this case, since the user is better off calling cond on the original matrix AFAICT.