A fully pivoted QR-factorization of a (m,n)-matrix A of floating point real or complex numbers has the purpose to provide a de-composition
(1) Pr * A * Pc = Q * R
where Pr and Pc are permutations of rows and columns of A, Q is unitary, and R is a quasi triangular matrix of the form.
(2) R = [R₁ R₂ ; Z₁ Z₂] # size(R₁) = (k,k), size(Z₂) = (m-k,n-k), Z₁, Z₂ are zero
(3) Q = [Q₁ Q₂] # size(Q₁) = (n,k), size(Q₂) = (n,m-k)
(4) Q * R = [Q₁*R₁ Q₁*R₂]
Here R₁ is upper square triangular with all diagonal elements nonzero (or positive) with a size of (k,k) where k = rank(A), size(R₂) = (k,n-k), and Z₁, Z₂ are zero matrices of appropriate sizes.
The decomposition is not unique. For numerical reasons the non-zeroness of the diagonal of R₁ has to be approximated by allowing a tolerance.
About R₂ there is not much to say; a column permutation strategy would try to keep the norms of the columns in R₂ smaller than those of R₂.
A “pseudo” solution of the system
(5) A * x ≈ b # we cannot expect "=" in general but minimize ∥A*x-b∥ in the 2-norm
could be defined by
(6) x₁ = inv(R₁) * ( Q₁' * Pr' * b - R₂ * x₂ )
with
(7) x = Pc * [x₁; x₂]
Here n-k of the x-components can be set to arbitrary values, while k x-componets are determined by those and R₁, Q₁, b. If k=n there is no ambiguity.
A typical implementation would set x₂ = 0 and calculate x₁ according to (6).
The “residue” b-A*x is given by Q₂ * Q₂' * Pr * b and its 2-norm ∥Q₂' * Pr * b∥ cannot be reduced by any choice of x. If k=m the residue is zero.
If no pivoting is done, the permutation matrices are both identity and the diagonal elements of R¹ may be zero or small. There is no numerical instability in the calculation of the components, but the further usability of R may be restricted. The size of R¹ may be greater than rank(A), if there are zero diagonal elements. R¹ may have zero diagonal elements, even if rank(A)=m. So working with pivot=false should be done only if rank(A)=n can be assumed. That is usually the case for least-squares fitting problems whith m≥n.
The implementation in Julia is done by the methods of qrfact. In Julia0.7 there are at least 5 implementations according to the values and types of arguments.
a. qrfact(A [, Val(false)]) when A is a dense matrix of BlasFloat element type
b. qrfact(A, Val(true)) when A is a dense matrix of BlasFloat element type
c. qrfact(A [, Val(false)]) when A is dense with other types
d. qrfact(A, Val(true)) when A is dense with other types
e. qrfact(A [, tol]) when A is sparse with element types Float64 and Complex{Float64}
a. and b. are implemented by the LAPACK C-library, c.and d. in the LinearAlgebra package of Julia, and e. by the SPQR C-library.
Notes.
The dense variants allow to select pivoting true, when only column pivoting is performed, while the sparse variant always performs column- and row-pivoting.
Only the sparse variant need to use a tolerance value to check if a pivot element is acceptable or not, because the other variants either accept zeros (if pivot=false) or select the maximal column. For the sparse case the criteria for pivot selection is the sparsity of the resulting matrix rather than the size of the pivot elements.
BlasFloat types are Float32/64 and Complex{Float32/64}.
The generic cases c. and d. support all types T, which can evaluate at least zero(T)/sqrt(abs2(one(T))).
The return types of qrfact are different subtypes of Factorization
C | typeof(F) | typeof(F.Q) | typeof(F.R)
- | --------- | --------- | -----------
a. | LinearAlgebra.QRCompactWY{Float64} | LinearAlgebra.QRCompactWYQ{Float64} | Array{Float64,2}
b. | QRPivoted{Float64} | LinearAlgebra.QRPackedQ{Float64} | Array{Float64,2}
c. | QR{BigFloat} | LinearAlgebra.QRPackedQ{BigFloat} | Array{Float64,2}
d. | QRPivoted{Float16} | LinearAlgebra.QRPackedQ{Float16} | Array{Float64,2}
e. | SuiteSparse.SPQR.QRSparse{Float64,Int64} | SuiteSparse.SPQR.QRSparseQ{Float64,Int64} |SparseMatrixCSC{Float64,Int64}
The following code prints information about the sizes of F.Q and F.R:
using LinearAlgebra
using SparseArrays
types = (("a", Matrix, Float64),
("b", Matrix, Float64, Val(true)),
("c", Matrix, BigFloat),
("d", Matrix, Float16, Val(true)),
("e", sparse, Float64) )
function sizelist(A)
println("sz | F | Q | Matrix(Q)| R | Q'*A | Q*R | R'*Q' | Q'\\A | F\\A")
println("-- | ----- | ------ | ------ | ------ | ------ | ------ | ------ | ------ | ----")
sizm(X, Y) = try size(X * Y) catch; "ERROR." end
sizd(X, Y) = try size(X \ Y) catch; "ERROR." end
sizf(X) = try size(Matrix(X)) catch; "ERROR." end
for t in types
(case, AT, ET) = t[1:3]
F = qrfact(AT(ET.(A)), t[4:end]...)
println("$case. | $(size(F)) | $(size(F.Q)) | $(sizf(F.Q)) | $(size(F.R)) | $(sizm(F.Q', A)) | $(sizm(F.Q, F.R)) | $(sizm(F.R', F.Q')) | $(sizd(F.Q', A)) | $(sizd(F, A))")
end
end;
A = [1 2 1 1 1 1 0; 1 2 0 1 1 1 0; 1 2 0 1 1 1 0; 0 0 0 0 0 0 0];
size(A) = (4,7) and rank(A) = 2.
Testing case m<n
julia> sizelist(A)
| sz | F | Q | Matrix(Q) | R | Q’*A | Q*R | R’*Q’ | Q’\A | F\A |
|---|---|---|---|---|---|---|---|---|---|
| a. | (4, 7) | (4, 4) | (4, 4) | (4, 7) | (4, 7) | (4, 7) | (7, 4) | (4, 7) | ERROR. |
| b. | (4, 7) | (4, 4) | (4, 4) | (4, 7) | (4, 7) | (4, 7) | (7, 4) | (4, 7) | (7, 7) |
| c. | (4, 7) | (4, 4) | (4, 4) | (4, 7) | (4, 7) | (4, 7) | (7, 4) | (4, 7) | ERROR. |
| d. | (4, 7) | (4, 4) | (4, 4) | (4, 7) | (4, 7) | (4, 7) | (7, 4) | (4, 7) | ERROR. |
| e. | (4, 7) | (4, 4) | (4, 1) | (2, 7) | (4, 7) | ERROR. | ERROR. | (4, 7) | (7, 7) |
This table has inconsistencies in some of the sizes of row e and exceptions in row e and column s(F\A).
The size(Q) == size(Matrix(Q)) should also hold for the sparse case. The deviation is due to a bug #26367 in Matrix().
The size(F.R) = (2,6) in the sparse case deviates from the other cases. It is not clear, what should be the correct value. See #26368.
The size of R is (m,n) in the dense cases, and (k,n) in the sparse case (where k = rank(A)). A PR is on the way to change the size in the sparse case to (m,n) as well.
The errors in row e are a consequence of the outcome of size(F.R). But it could be argued, that multiplication with Q-matrices should also allow right-hand sides with row count identical to the rank of A: they could be padded with zeros, which is a common use case when working with Q-R factorizations.
The error in row a (pivot=false, LAPACK), throws DimensionMismatch("matrix is not square: dimensions are (4, 6)" in ldiv! .../qr.jl:764 that does not look as the right exception.
The error in row c (pivot=false, generic), throws SingularException(4) in qr..jl:833 . That looks OK, because R has small diagonal elements R[i,i] with i <= rank(A).
The error in row d (pivot=true, generic), throws SingularException(4) in qr..jl:833.
That looks not optimal, because R has no small diagonal elements R[i,i] with i <= rank(A) and rows b and e (pivot=true, LAPACK and sparse) prove, that a better solution is possible.
(Solutions in b and e are not (approximately) the same, but for both A * ( F \ A) ≈ A).
Now for m>n, the situation is quite different:
julia> sizelist(A')
| sz | F | Q | Matrix(Q) | R | Q’*A | Q*R | R’*Q’ | Q’\A | F\A |
|---|---|---|---|---|---|---|---|---|---|
| a. | (7, 4) | (7, 7) | (7, 4) | (4, 4) | (7, 4) | (7, 4) | ERROR. | (7, 4) | ERROR. |
| b. | (7, 4) | (7, 7) | (7, 4) | (4, 4) | (7, 4) | (7, 4) | ERROR. | (7, 4) | (4, 4) |
| c. | (7, 4) | (7, 7) | (7, 4) | (4, 4) | (7, 4) | (7, 4) | ERROR. | (7, 4) | ERROR. |
| d. | (7, 4) | (7, 7) | (7, 4) | (4, 4) | (7, 4) | (7, 4) | ERROR. | (7, 4) | ERROR. |
| e. | (7, 4) | (7, 7) | (7, 2) | (2, 4) | (7, 4) | ERROR. | ERROR. | (7, 4) | ERROR. |
In this case, there are additional inconsitencies:
Case e. the size of Matrix(Q) deviates from cases a.-d. due to the same bug.
This time size(Q) != size(Matrix(Q)) in all cases. That is unexpected and yet undocumented.
Q'*A' returns a (7,4) matrix in all cases, while Matrix(Q)'*A would return a (4,4)-matrix, what is a consequence of the previous item and probably unexpected yet undocumented.
size(R) = (4,4) in cases a.-d. and (2,4) in case e.. That is the same deviation, as in the first table and should be made consistent.
Q*R is possible in cases a.-d. despite the contradicting sizes of Q and R. That makes sense due to a common use case, and because we want to support Q * R = A[:,F.p].
If Q*R is possible (and it should!), the also R'*Q' should succeed, because R'*Q' = (Q*R)'. That it throws an exeption is a bug in my opinion.
Column F\A is even worse than in the first table, because here also the sparse case fails. A proper implementation should deliver useful results at least in cases b., d., e., which use pivoting.
In my opionion, ther should be the following consequences.
-
Decide, if the row count of
Rusesmin(m,n)as is in the dense cases today orrank(A)as is in the sparse cases. -
Store this value as a new field in all types of
AbstractQ. -
Make
size(Matrix(Q)) == (m, size(R,1)), whatever the decision in 1. was. -
Allow multiplications
Q*BandB'*Q'whensize(B,1) in (m, size(R,1), min(m,n)). -
Let
F\Breturn useful results and not exceptions in all pivoted cases. -
Provide a public method
rank(F)which returns useful values in all pivoted cases.
