About QR decomposition - linear algebra


#1

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 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 may be greater than rank(A), if there are zero diagonal elements. 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.

  1. Decide, if the row count of R uses min(m,n) as is in the dense cases today or rank(A) as is in the sparse cases.

  2. Store this value as a new field in all types of AbstractQ.

  3. Make size(Matrix(Q)) == (m, size(R,1)), whatever the decision in 1. was.

  4. Allow multiplications Q*B and B'*Q' when size(B,1) in (m, size(R,1), min(m,n)).

  5. Let F\B return useful results and not exceptions in all pivoted cases.

  6. Provide a public method rank(F) which returns useful values in all pivoted cases.


#2

I am following the incremental edits with interest, but it is not clear to me if this is a question, a suggestion, or something else.


#3

It’s a nice write-up regardless :slight_smile:


#4

It is thought to be a clarification of the types and sizes of the results of QR factorizations F and about the operations, which are possible with F, F.Q, F.R, …
Related to #26368


#5

If the A is invertible and the diagonal in R_1 is chosen to be positive, then it is unique.

Addition:

Q itself is divided into [Q_1,Q_2] where Q_1 is an n x k matrix and Q_2 is n x (n-k). The column space of Q_1 is the column space of A, and the columns of Q_2 are orthogonal to the column space of A, which makes them basis of the null space of A^T.

Also for sparse matrices, the Q matrix is not formed explicitly unless required, it is stored as a sequence of low rank operators, which can be multiplied by a vector, and can be transposed easily without forming the actual matrix.


#6

All true if m == n, but let me finish my text, which is not at first about linear algebra, but the qrfact implementations.


#7

OK, I get the context now, especially given the issue.

Composing in an editor first and posting something finished would be less confusing to readers, but I ain’t complaining, since this is a really nice writeup.


#8

can you recommend an (offline-) md-previewer for Linux?


#9

I use various Emacs extensions, but if you are not into that, retext is a simple solution (but please don’t get me wrong, I appreciate this writeup, just didn’t know what was going on).


#10

I’m using https://typora.io/ on Windows and it also seems to work on Linux.


#11

A Jupyter notebook markdown cell works pretty well, too.


#12

I like VS code for this: https://code.visualstudio.com/Docs/languages/markdown#_editor-and-preview-synchronization. This works out of the box. For math in markdown (which discourse supports), you can install this extension: https://marketplace.visualstudio.com/items?itemName=goessner.mdmath.


#13

You could probably try atom editor and markdown preview package