18) Householder QR

Last time¶

• Gram-Schmidt orthogonalization

• Classical vs Modified Gram-Schmidt

• QR factorization

Today¶

1. Recap from last time

2. Householder QR

1. Recap Gram-Schmidt orthogonalization¶

For many applications, we find ourselves interested in the column spaces of a matrix $A$:

The idea of QR factorization is the construction of a sequence of orthonormal vectors, $q_1, q_2, \ldots$ that span these successive spaces.

Thus, suppose we want to find an orthogonal basis for the span of the columns of $A$:

Given $a_1, a_2, \dots$, we can construct vectors $q_1, q_2, \ldots$ and entries $r_{ij}$, by an iterative process of successive orthogonalization.

Gram-Schmidt with more parallelism¶

We are performing the following (from right to left):

which can be factored in blocks to expoloit a bit of parallelism. Let’s call this classical Gram-Schmidt.

Classical Vs Modified Gram-Schmidt¶

Why does the order of operations matter?¶

We are performing the following (from right to left):

which can be factored in blocks and is not exact in finite arithmetic.

Let’s look at this a bit more closely. In the classical Gram-Schmidt (CGS), we take each vector, one at a time, and make it orthogonal to all previous vectors. If an error is made in computing $q_2$ in CGS, so that $q^T_1q_2=\delta$ is small, but non-zero in finite arithmetic, this will not be corrected for in any of the computations that follow, and it will actually be propagated. In this case, $v_3$ will not be orthogonal to $q_1$ or $q_2$.

Because of rounding errors, in classical Gram-Schmidt, $Q_{k−1}$ does not have truly orthogonal columns. In modified Gram–Schmidt (MGS) we compute the length of the projection of $w = v_k$ onto $q_1$ and subtract that projection (and the rounding errors) from $w$. Next, we compute the length of the projection of the computed $w$ onto $q_2$ and subtract that projection (and the rounding errors) from $w$, and so on, but always orthogonalizing against the computed version of $w$.

We can look at the size of what’s left over¶

We project out the components of our vectors in the directions of each $q_j$.

The next vector is almost linearly dependent¶

We can see that $v$ has substantial components in some $Q$-directions, while it is nearly orthogonal to the other directions.

Right-looking modified Gram-Schmidt¶

Each outer step of the modified Gram-Schmidt algorithm can be interpreted as a right-multiplication by a square upper-triangular matrix.

Classical versus modified?¶

• Classical

  • Really unstable, orthogonality error of size $1 \gg \epsilon_{\text{machine}}$

  • Don’t need to know all the vectors in advance

• Modified

  • Needs to be right-looking for efficiency

  • Less unstable, but orthogonality error $10^{-9} \gg \epsilon_{\text{machine}}$

2. Householder triangularization¶

Householder triangularization is numerically more stable than Gram-Schmidt orthogonaliztion, though it lacks the latter’s applicability as a basis for iterative methods. The Householder algorithm is a process of “orthogonal triangularization”, making a matrix triangular by a sequence of unitary matrix operations.

Gram-Schmidt constructed a triangular matrix $R$ to orthogonalize $A$ into $Q$. Each step was an orthogonal projector, which is a rank-deficient operation.

Image credits: Pictures taken from the reference book: Trefethen and Bau (1997), Numerical Linear Algebra (available at the library).

Householder uses orthogonal transformations (reflectors) to triangularize.

The reflection, as depicted above by Trefethen and Bau (1997) can be written $F = I - 2 \frac{v v^T}{v^T v}$.

The structure of the algorithm is

Constructing the $Q_j$¶

Each of our $Q_j$ will have the form

where $F$ is a “reflection” that achieves

where $x$ is the column of $R$ from the diagonal down. This transformation is a reflection across a plane with normal $v = Fx - x = \lVert x \rVert e_1 - x$.

The two methods can be summarized as follows:

• Gram-Schmidt: a triangular orthogonalization

• Householder: an orthogonal triangularization

Adventures in reflection¶

Householder: A naive algorithm¶

How to interpret $V$ as $Q$?¶

The following two programs compute a matrix containing vectors that generate the Householder reflectors whose product is Q.

Great, but we can still break it¶

We had the lines

    v = copy(R[j:end, j])
    v[1] -= norm(v)
    v = normalize(v)

What happens when R is already upper triangular?

In this case

In general, there are two possible reflections. For numerical stability, it is important to choose the one that moves $x$ the largest distance, so that it is not too close to $x$ itself.

In fact, suppose that in the above figure, the angle of $H^{+}$ and the $e_1$ axis is much smaller than $x$ or $|| x || e_1$. Thus, the calculation of $v$ represents a subtraction of nearby quantities and will tend to suffer from cancellation errors. We can pick the sign of $v$ so that we can avoid such effects by ensuring that $||v||$ is never smaller than $||x||$.

Householder: An improved algorithm¶

Householder is backward stable¶

Orthogonality is preserved¶

Summary:¶

• Classic Gram-Schmidt: Usually very poor orthogonality.

• Modified Gram-Schmidt: Depends upon condition of $A$. Fails completely when $A$ is singular.

• Householder triangularization: Always good orthogonality and backward stable.