User:Fephisto/Householder QR Decomposition

QR decomposition

Householder transformations can be used to calculate a QR decomposition. Consider a matrix tridiangularized up to column $i$, then our goal is to construct such Householder matrices that act upon the principal submatrices of a given matrix

${\begin{bmatrix}a_{11}&a_{12}&\cdots &&&a_{1n}\\0&a_{22}&\cdots &&&a_{1n}\\\vdots &&\ddots &&&\vdots \\0&\cdots &0&x_{1}=a_{ii}&\cdots &a_{in}\\0&\cdots &0&\vdots &&\vdots \\0&\cdots &0&x_{n}=a_{ni}&\cdots &a_{nn}\end{bmatrix}}$

via the matrix

${\begin{bmatrix}I_{i-1}&0\\0&P_{v}\end{bmatrix}}$.

(note that we already established before that Householder transformations are unitary matrices, and since the multiplication of unitary matrices is itself a unitary matrix, this gives us the unitary matrix of the QR decomposition)

If we can find a ${\vec {u}}$ so that

$P_{v}{\vec {x}}={\vec {e}}_{1}$

we could accomplish this. Thinking geometrically speaking, we are looking for a plane so that the reflection of ${\overline {v}}$ about the plane happens to land directly on the basis vector. In other words,

${\vec {x}}-2(<{\vec {x}},{\vec {v}}>){\vec {v}}=\alpha {\vec {e}}_{1}$ 1

for some constant $\alpha$. However, this also shows that

${\vec {v}}\propto {\vec {x}}-\alpha {\vec {e}}_{1}$.

And since ${\vec {v}}$ is a unit vector, this means that

${\vec {v}}=\pm {\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\vert \vert {\vec {x}}-\alpha {\vec {e}}_{1}\vert \vert _{2}}}$ 2

Now if we apply equation (2) back into equation (1).

${\vec {x}}-\alpha {\vec {e}}_{1}=2(<{\vec {x}},{\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\vert \vert {\vec {x}}-\alpha {\vec {e}}_{1}\vert \vert _{2}}}>){\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\vert \vert {\vec {x}}-\alpha {\vec {e}}_{1}\vert \vert _{2}}}$

Or, in other words, by comparing the scalars in front of the vector ${\vec {x}}-\alpha {\vec {e}}_{1}$ we must have

$\vert \vert {\vec {x}}-\alpha {\vec {e}}_{1}\vert \vert _{2}^{2}=2(<{\vec {x}},{\vec {x}}-\alpha e_{1}>)$.

Or

$2(\vert \vert {\vec {x}}\vert \vert _{1}^{2}-\alpha x_{1})=\vert \vert {\vec {x}}\vert \vert _{2}^{2}-2\alpha x_{1}+\alpha ^{2}$

Which means that we can solve for $\alpha$ as

$\alpha =\pm \vert \vert {\vec {x}}\vert \vert _{2}^{2}$

This completes the construction; however, in order to avoid catastrophic cancellation in equation (2) we choose the sign of $\alpha$ so that

$\alpha =sign(Re(x_{1}))\vert \vert {\vec {x}}\vert \vert _{2}^{2}$