The Householder Reduction

The Householder rotation P has the following form:

$$\boldsymbol{P} = \boldsymbol{1} - 2 \boldsymbol{w}\otimes\boldsymbol{w}$$ (3.40)

where 1 is the Kronecker delta $\delta_{ij}$ and $\boldsymbol{w}\otimes\boldsymbol{w}$ is the tensor product, i.e., $\boldsymbol{w}\otimes\boldsymbol{w} \equiv w_i w_j$.

Of vector w we demand only that $\boldsymbol{w}\cdot\boldsymbol{w} = 1$.

Matrix P is clearly symmetric, hence P^T = P.

Furthermore matrix P is its own inverse, and this can be seen as follows:

$$\begin{eqnarray*}\boldsymbol{P} \cdot \boldsymbol{P} &=& \left( \boldsymbol{1} ... ... \left(\boldsymbol{w}\otimes\boldsymbol{w}\right) \\ &=& \ldots \end{eqnarray*}$$

Now, observe that

$$\left(\boldsymbol{w}\otimes\boldsymbol{w}\right) \cdot \left(... ...j w_j w_k = w_i w_k \equiv \boldsymbol{w}\otimes\boldsymbol{w}$$

and, of course, $\boldsymbol{1}\cdot\boldsymbol{1} = \boldsymbol{1}$, hence

$$\ldots = \boldsymbol{1} - 4 \boldsymbol{w}\otimes\boldsymbol{w} + 4 \boldsymbol{w}\otimes\boldsymbol{w} = \boldsymbol{1}$$

Consequently: P = P^T = P^-1, so P is orthogonal.

We can use any vector uin place of w if we normalize it at the same time:

$$\boldsymbol{P} = \boldsymbol{1} - \frac{2\boldsymbol{u}\otime... ...\boldsymbol{1} - \frac{\boldsymbol{u}\otimes\boldsymbol{u}}{H}$$ (3.41)

where

$$H = \frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}$$ (3.42)

So far P isn't specific enough to deserve a proud name of a Housholder matrix. But now we choose u to be

$$\boldsymbol{u} = \boldsymbol{x} \mp \vert\boldsymbol{x}\vert\boldsymbol{e}_1$$ (3.43)

where e₁ is the first basis vector.

The Householder operator with the u vector so defined rotates vector x onto e₁, and this can be seen as follows:

 

$$\boldsymbol{P}\cdot\boldsymbol{x} \left( \boldsymbol{1} - \frac{\boldsymbol{u}\otimes\boldsymbol{u}}{H} \right) \cdot \boldsymbol{x}$$ =

$$\left( \boldsymbol{1} - \frac{2\boldsymbol{u}\otimes \left(\bolds... ...l{e}_1\right)} {\boldsymbol{u}\cdot\boldsymbol{u}} \right) \cdot \boldsymbol{x}$$ =

What's $\boldsymbol{u}\cdot\boldsymbol{u}$?

 

$$\boldsymbol{u}\cdot\boldsymbol{u} \left(\boldsymbol{x}\mp\vert\boldsymbol{x}\vert\boldsymbol{e}_1\right) \cdot \left(\boldsymbol{x}\mp\vert\boldsymbol{x}\vert\boldsymbol{e}_1\right)$$ =

$$\boldsymbol{x}\cdot\boldsymbol{x} \mp 2 \vert\boldsymbol{x}\vert\... ...oldsymbol{x} + \vert\boldsymbol{x}\vert^2 \boldsymbol{e}_1\cdot\boldsymbol{e}_1$$ =

$$\vert\boldsymbol{x}\vert^2 \mp 2 \vert\boldsymbol{x}\vert x_1 + \vert\boldsymbol{x}\vert^2$$ =

$$2 \vert\boldsymbol{x}\vert^2 \mp 2 \vert\boldsymbol{x}\vert x_1$$ = (3.44)

Substituting equation (3.45) into (3.44) yields:

$$\boldsymbol{P}\cdot\boldsymbol{x} \left( \boldsymbol{1} - \frac{2\boldsymbol{u}\otimes \left(\bolds... ...mbol{x}\vert^2 \mp 2 \vert\boldsymbol{x}\vert x_1} \right) \cdot \boldsymbol{x}$$ =

$$\boldsymbol{x} - \frac{2\boldsymbol{u} \left( \vert\boldsymbol{x}... ...t x_1 \right)}{2 \vert\boldsymbol{x}\vert^2 \mp 2 \vert\boldsymbol{x}\vert x_1}$$ =

= x - u

$$\boldsymbol{x} - \boldsymbol{x} \pm \vert\boldsymbol{x}\vert\boldsymbol{e}_1$$ =

$$\pm \vert\boldsymbol{x}\vert\boldsymbol{e}_1$$ = (3.45)

Remember that P is not a projection operator. Projection operators are not orthogonal. P is a rotation, which rotates vector x onto the e₁ direction. In the process the length of the vector does not change. Whereas it would have changed for a projection operator.

The procedure for transforming matrix A into a tridiagonal form works as follows.

The first Householder operator, P₁ is selected to rotate the sub-column of the first column:

$$\left( \begin{array}{c} a_{21} \\ a_{21} \\ \vdots \\ a_{n1} \end{array}\right)$$

onto

$$\left( \begin{array}{c} a'_{21} \\ 0 \\ \vdots \\ 0 \end{array}\right)$$

And, to accomplish that, the operator has to have the following form:

$$\boldsymbol{P}_1 = \left( \begin{array}{ccccc} 1 & 0 & 0 & ... ...}^{(n-1)}\boldsymbol{P}_1 & \\ 0 & & & & \end{array}\right)$$

Multiplying A by P₁ from the left:

$$\boldsymbol{P}_1 \cdot \boldsymbol{A} = \left( \begin{array}... ...& \\ \vdots & & & & \\ a_{n1} & & & & \end{array}\right)$$

attacks the first sub-column with ^(n-1)P₁and if that sub-operator has been chosen to rotate the first sub-column onto its first dimension then the resulting matrix A' will look as follows:

$$\left( \begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \cdo... ... & & & \\ \vdots & & & & \\ 0 & & & & \end{array}\right)$$

Observe that the first row will not be changed by that operation at all, but the lower right corner, of course, will get changed. But now if we apply operator P₁ to A' from the right the same thing will happen to the first row (remember that P is symmetric), so that the resulting matrix A'' will look as follows:

$$\boldsymbol{P}_1 \cdot \boldsymbol{A} \cdot \boldsymbol{P}_1 ... ... & & & \\ \vdots & & & & \\ 0 & & & & \end{array}\right)$$

The second Householder matrix is going to look as follows:

$$\boldsymbol{P}_2 = \left( \begin{array}{ccccc} 1 & 0 & 0 & ... ...(n-2)}\boldsymbol{P}_2 & \\ 0 & 0 & & & \end{array}\right)$$

and, as you should understand by now, it will do to

$$\left( \begin{array}{c} a_{32} \\ a_{42} \\ \vdots \\ a_{n2} \end{array}\right)$$

the same as P₁ did to

$$\left( \begin{array}{c} a_{21} \\ a_{31} \\ a_{41} \\ \vdots \\ a_{n1} \end{array}\right)$$

In other words it will rotate it onto its first direction, thus leaving only one term under the diagonal. The identity block in the upper left corner ensures that tridiagonalization achieved so far will not be spoiled during successive rotations.

It is now obvious that in n-2 such steps the whole matrix must become triadiagonalized.

Because $\boldsymbol{P} = \boldsymbol{1} - \frac{\boldsymbol{u}\otimes\boldsymbol{u}}{H}$ we can write down our operations on A in more detail. First

$$\boldsymbol{A}\cdot\boldsymbol{P} = \boldsymbol{A}\cdot\left... ...rac{\boldsymbol{A}\cdot\boldsymbol{u}}{H}\otimes\boldsymbol{u}$$ (3.46)

Introducing a new vector p:

$$\boldsymbol{p} = \frac{\boldsymbol{A}\cdot\boldsymbol{u}}{H}$$ (3.47)

we get

$$\boldsymbol{A}\cdot\boldsymbol{P} = \boldsymbol{A} - \boldsymbol{p} \otimes\boldsymbol{u}$$ (3.48)

So now:

$$\boldsymbol{P}\cdot\boldsymbol{A}\cdot\boldsymbol{P} \left(\boldsymbol{1} - \frac{\boldsymbol{u}\otimes\boldsymbol{u}}... ...ght) \cdot \left( \boldsymbol{A} - \boldsymbol{p} \otimes\boldsymbol{u} \right)$$ =

$$\boldsymbol{A} - \boldsymbol{p} \otimes \boldsymbol{u} - \boldsym... ...{H} \boldsymbol{u}\otimes\boldsymbol{u}\cdot\boldsymbol{p}\otimes\boldsymbol{u}$$ =

$$\boldsymbol{A} - \boldsymbol{p} \otimes \boldsymbol{u} - \boldsym... ... \frac{\boldsymbol{u}\cdot\boldsymbol{p}}{H}\boldsymbol{u}\otimes\boldsymbol{u}$$ =

$$\boldsymbol{A} - \boldsymbol{p} \otimes \boldsymbol{u} - \boldsymbol{u}\otimes\boldsymbol{p} + 2 K \boldsymbol{u}\otimes\boldsymbol{u}$$ = (3.49)

where

$$K = \frac{\boldsymbol{u}\cdot\boldsymbol{p}}{2H}$$ (3.50)

Our expression for $\boldsymbol{P}\cdot\boldsymbol{A}\cdot\boldsymbol{P}$ can be further simplified by introducing vector q:

q = p - Ku (3.51)

and observing that:

$$\boldsymbol{P}\cdot\boldsymbol{A}\cdot\boldsymbol{P} = \bolds... ...{q}\otimes\boldsymbol{u} - \boldsymbol{u}\otimes\boldsymbol{q}$$ (3.52)

Let us summarise the flow of computation now. For a particular square submatrix of A we'll have

$$\begin{eqnarray*}\sigma &=& a_{21}^2 + a_{31}^2 + a_{41}^2 + \cdots + a_{n1}^2 \... ...ol{q}\otimes\boldsymbol{u} -\boldsymbol{u}\otimes\boldsymbol{q} \end{eqnarray*}$$

and the accumulated transform, which we are going to need for the eigenvectors is:

$$\boldsymbol{Q} = \boldsymbol{P}_1\cdot\boldsymbol{P}_2 \cdot\ldots\cdot\boldsymbol{P}_{n-1}$$ (3.53)

Of course, you now appreciate that if eigenvectors are needed, this is going to cost us additional operations, because we'll have to compute $\boldsymbol{P} = \boldsymbol{1} - \boldsymbol{u}\otimes\boldsymbol{u}/ H$ for every rotation, whereas without eigenvectors we can get away with just the $\sigma$, u, H, p, K, q, A' sequence.