Add nullspace and cond functions for QRPivoted matrices #54519
Conversation
Co-authored-by: Daniel Karrasch <[email protected]>
Co-authored-by: Daniel Karrasch <[email protected]>
| 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 |
There was a problem hiding this comment.
Since both Q and the permutation don't change norms, isn't this simply
| 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.
There was a problem hiding this comment.
There was a problem hiding this comment.
Q doesn't change 2-norms, but it changes 1 norms and infinity norms.
There was a problem hiding this comment.
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).
| 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]) |
There was a problem hiding this comment.
| 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.
There was a problem hiding this comment.
This algorithm seems incorrect. If cond(A) == cond(R), but cond(R) is not the ratio of diagonal entries of R.
| R = UpperTriangular(A.R) | ||
| try | ||
| Rinv = inv(R) | ||
| return opnorm(R, p)*opnorm(Rinv, p) |
There was a problem hiding this comment.
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.
This PR closes JuliaLang/LinearAlgebra.jl#1068 and also closes JuliaLang/LinearAlgebra.jl#855.