# About QR decomposition - linear algebra

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.

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.

3 Likes

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.

It’s a nice write-up regardless

2 Likes

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

1 Like

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

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.

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

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.

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

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).

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

A Jupyter notebook markdown cell works pretty well, too.

I like VS code for this: Markdown editing with Visual Studio Code. This works out of the box. For math in markdown (which discourse supports), you can install this extension: Markdown+Math - Visual Studio Marketplace.

You could probably try atom editor and markdown preview package

1 Like