17) Gram-Schmidt and QR

Last time¶

• Linear Algebra

• Polynomial evaluation

• Orthogonality

Today¶

1. Revisit orthogonality

2. Constructing orthogonal bases
   2.1 Orthogonal matrices
   2.2 Gram-Schmidt orthogonalization

3. QR factorization

1. Revisit orthogonality¶

Polynomials can be orthogonal too!¶

• Which inner product will be zero?

  • Which functions are even and odd?

Pairwise inner product¶

The pairwise inner products between two sets of vectors can be expressed by collecting the sets as columns in matrices and writing $A = X^T Y$ where $A_{i,j} = x_i^T y_j$. It follows from this definition that

Normalization and orthogonalization¶

2. Orthogonality¶

2.1 Orthogonal matrices¶

If two vectors $x$ and $y$ are such that $x^T y = 0$ then we say $x$ and $y$ are orthogonal (or “$x$ is orthogonal to $y$”).

A vector is said to be normalized if $\lVert x \rVert = 1$. If $x$ is orthogonal to $y$ and $\lVert x \rVert = \lVert y \rVert = 1$ then we say $x$ and $y$ are orthonormal.

A square matrix with orthonormal columns is said to be an orthogonal matrix (or orthonormal matrix).

We typically use $Q$ or $U$ and $V$ for matrices that are known/constructed to be orthogonal.

Orthogonal matrices are always full rank -- the columns are linearly independent.

The inverse of an orthogonal matrix is its transpose:

Orthogonal matrices are a powerful building block for robust numerical algorithms.

2.2 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.

A naive algorithm¶

What do orthogonal polynomials look like?¶

Computations with orthogonal polynomials form the basis of spectral methods, one of the most powerful techniques for the numerical solution of partial differential equations (PDE).

What happens if we use more than 50 values of $x$?

Theorem:¶

Every full-rank $m\times n$ matrix ($m \ge n$) has a unique reduced $Q R$ factorization with $R_{j,j} > 0$¶

The algorithm we’re using generates this positive matrix due to the line:

        R[j,j] = norm(v)

Solving equations using $QR = A$¶

If $A x = b$ then $Rx = Q^T b$ (because Q is orthogonal, and it is easy to solve with $R$ since is triangular).

How accurate is it?¶

A variant with more parallelism¶

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

Here $I-q_iq_i^T$ removes the component of a vector in the direction $q_i$, assuming $q_i$ is a unit vector. From right to left, we are first applying $I-q_1q_1^T$, which removes the $q_1$ component. Then we apply $I-q_2q_2^T$, which removes the $q_2$ component. In the above, the cross term vanishe as $q_1 \perp q_2$.

This operation 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.

The order of operations matter!

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}}$