# About QR decomposition - linear algebra

**URL:** <https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606>\
**Category:** Numerics\
**Created:** [March 9, 2018, 6:44am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606 "2018-03-09T06:44:28Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [March 9, 2018, 6:44am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/1 "2018-03-09T06:44:29Z")

</div>

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

```julia
(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.

```julia
(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

```julia
(5) A * x ≈ b # we cannot expect "=" in general but minimize ∥A*x-b∥ in the 2-norm

```

could be defined by

```julia
(6) x₁ = inv(R₁) * ( Q₁' * Pr' * b - R₂ * x₂ )

```

with

```julia
(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.

```julia
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`:

```julia
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
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](https://github.com/JuliaLang/julia/issues/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](https://github.com/JuliaLang/julia/issues/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
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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 9, 2018, 10:43am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/2 "2018-03-09T10:43:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 9, 2018, 10:49am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/3 "2018-03-09T10:49:34Z")

</div>

It’s a nice write-up regardless 🙂

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [March 9, 2018, 10:53am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/4 "2018-03-09T10:53:58Z")

</div>

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](https://github.com/JuliaLang/julia/issues/26368)

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 9, 2018, 11:00am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/5 "2018-03-09T11:00:20Z")

</div>

> [@klacru](#):
>
> The decomposition is not unique.

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.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [March 9, 2018, 11:05am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/6 "2018-03-09T11:05:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 9, 2018, 11:12am UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/7 "2018-03-09T11:12:28Z")

</div>

> [@klacru](#):
>
> clarification of the types and sizes of the results of QR factorizations

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

> [@klacru](#):
>
> but let me finish my text

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.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [March 9, 2018, 1:51pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/8 "2018-03-09T13:51:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 9, 2018, 1:54pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/9 "2018-03-09T13:54:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [March 9, 2018, 1:55pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/10 "2018-03-09T13:55:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jmert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmert/32/3161_2.png) [@jmert](https://discourse.julialang.org/u/jmert)\
**Post date:** [March 9, 2018, 3:08pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/11 "2018-03-09T15:08:36Z")

</div>

A Jupyter notebook markdown cell works pretty well, too.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [March 9, 2018, 3:24pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/12 "2018-03-09T15:24:27Z")

</div>

I like VS code for this: [Markdown editing with Visual Studio Code](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: [Markdown+Math - Visual Studio Marketplace](https://marketplace.visualstudio.com/items?itemName=goessner.mdmath).

---

<div class="post-metadata">

**Author:** ![Liso](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@Liso](https://discourse.julialang.org/u/Liso)\
**Post date:** [March 9, 2018, 4:58pm UTC](https://discourse.julialang.org/t/about-qr-decomposition-linear-algebra/9606/13 "2018-03-09T16:58:45Z")

</div>

You could probably try atom editor and [markdown preview package](https://atom.io/packages/markdown-preview)