97 lines
3.5 KiB
Plaintext
97 lines
3.5 KiB
Plaintext
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>
|
|
|