Any rectangular <math>m \times n</math> matrix <math>\mathit A</math> can be decomposed to a product of an orthogonal matrix <math>\mathit Q</math> and an upper (right) triangular matrix <math>\mathit R</math>, as described in [[wp:QR decomposition|QR decomposition]].

'''Task'''

Demonstrate the QR decomposition on the example matrix from the [[wp:QR_decomposition#Example_2|Wikipedia article]]:

::<math>A = \begin{pmatrix}
12 & -51 & 4 \\
6 & 167 & -68 \\
-4 & 24 & -41 \end{pmatrix}</math>

and the usage for linear least squares problems on the example from [[Polynomial regression]]. The method of [[wp: Householder transformation|Householder reflections]] should be used:

'''Method'''

Multiplying a given vector <math>\mathit a</math>, for example the first column of matrix <math>\mathit A</math>, with the Householder matrix <math>\mathit H</math>, which is given as

::<math>H = I - \frac {2} {u^T u} u u^T</math>

reflects <math>\mathit a</math> about a plane given by its normal vector <math>\mathit u</math>. When the normal vector of the plane <math>\mathit u</math> is given as

::<math>u = a - \|a\|_2 \; e_1</math>

then the transformation reflects <math>\mathit a</math> onto the first standard basis vector

::<math>e_1 = [1 \; 0 \; 0 \; ...]^T</math>

which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of <math>a_1</math>:

::<math>u = a + \textrm{sign}(a_1)\|a\|_2 \; e_1</math>

and normalize with respect to the first element:

::<math>v = \frac{u}{u_1}</math>

The equation for <math>H</math> thus becomes:

::<math>H = I - \frac {2} {v^T v} v v^T</math>

or, in another form

::<math>H = I - \beta v v^T</math>

with
::<math>\beta = \frac {2} {v^T v}</math>

Applying <math>\mathit H</math> on <math>\mathit a</math> then gives

::<math>H \; a = -\textrm{sign}(a_1) \; \|a\|_2 \; e_1</math>

and applying <math>\mathit H</math> on the matrix <math>\mathit A</math> zeroes all subdiagonal elements of the first column:

::<math>H_1 \; A = \begin{pmatrix}
r_{11} & r_{12} & r_{13} \\
0    & *    & * \\
0    & *    & * \end{pmatrix}</math>

In the second step, the second column of <math>\mathit A</math>, we want to zero all elements but the first two, which means that we have to calculate <math>\mathit H</math> with the first column of the ''submatrix'' (denoted *), not on the whole second column of <math>\mathit A</math>.

To get <math>H_2</math>, we then embed the new <math>\mathit H</math> into an <math>m \times n</math> identity:

::<math>H_2 = \begin{pmatrix}
1 & 0 & 0 \\
0 & H & \\
0 &   & \end{pmatrix}</math>

This is how we can, column by column, remove all subdiagonal elements of <math>\mathit A</math> and thus transform it into <math>\mathit R</math>.

::<math>H_n \; ... \; H_3 H_2 H_1 A = R</math>

The product of all the Householder matrices <math>\mathit H</math>, for every column, in reverse order, will then yield the orthogonal matrix <math>\mathit Q</math>.

::<math>H_1 H_2 H_3 \; ... \; H_n = Q</math>

The QR decomposition should then be used to solve linear least squares ([[Multiple regression]]) problems <math>\mathit A x = b</math> by solving

::<math>R \; x = Q^T \; b</math>

When <math>\mathit R</math> is not square, i.e. <math>m > n</math> we have to cut off the <math>\mathit m - n</math> zero padded bottom rows.

::<math>R =
\begin{pmatrix}
R_1 \\
0 \end{pmatrix}</math>

and the same for the RHS:

::<math>Q^T \; b =
\begin{pmatrix}
q_1 \\
q_2 \end{pmatrix}</math>

Finally, solve the square upper triangular system by back substitution:

::<math>R_1 \; x = q_1</math>