A Course on Numerical Linear Algebra¶

Instructor: Deniz Bilman

Textbook: Numerical Linear Algebra, by Trefethen and Bau.

Lecture 16: Stability of Householder Triangularization¶

In [5]:
using LinearAlgebra

# Define the size of the matrix
m = 50

# Generate a random matrix and convert to upper triangular
R = triu(rand(m, m))

Q, X = qr(randn(m, m))       # Compute QR factorization of some other random matrix to obtain Q

# Set A to be the product QR up to rounding errors (coming just from the product).

A = Q*R

# Compute QR factorization of A by Householder triangularization.
# In Julia, the qr function from the LinearAlgebra package performs QR factorization using Householder reflections by default. 
# This is the standard approach for dense matrices.

Q2, R2 = qr(A)

LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}
Q factor: 50×50 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
R factor:
50×50 Matrix{Float64}:
 0.31951  0.657193  0.0932518  0.955495  …  0.951841   0.460561  0.215868
 0.0      0.180917  0.844836   0.687551     0.839704   0.504348  0.728806
 0.0      0.0       0.0901526  0.573472     0.282417   0.934994  0.175477
 0.0      0.0       0.0        0.187741     0.17888    0.661017  0.632674
 0.0      0.0       0.0        0.0          0.649704   0.376047  0.183774
 0.0      0.0       0.0        0.0       …  0.789325   0.407612  0.289946
 0.0      0.0       0.0        0.0          0.737896   0.32322   0.454282
 0.0      0.0       0.0        0.0          0.783517   0.81968   0.59265
 0.0      0.0       0.0        0.0          0.270473   0.568733  0.05098
 0.0      0.0       0.0        0.0          0.283912   0.42023   0.919315
 ⋮                                       ⋱                       
 0.0      0.0       0.0        0.0          0.0938852  0.755723  0.386912
 0.0      0.0       0.0        0.0          0.962039   0.952543  0.502407
 0.0      0.0       0.0        0.0          0.844864   0.381722  0.990536
 0.0      0.0       0.0        0.0          0.86463    0.844068  0.205038
 0.0      0.0       0.0        0.0       …  0.958373   0.156297  0.190568
 0.0      0.0       0.0        0.0          0.0463459  0.342806  0.566779
 0.0      0.0       0.0        0.0          0.394683   0.770051  0.551017
 0.0      0.0       0.0        0.0          0.0        0.120719  0.346685
 0.0      0.0       0.0        0.0          0.0        0.0       0.479349
In [8]:
# How accurate is Q2?

norm(Matrix(Q2)-Matrix(Q))

0.08974853446442728

This is a disaster! We are working with double precision arithmetic --- $\epsilon_{\text{M}}\approx 2.22\times 10^{-16}$.

How about the relative error in computing $R$?

In [9]:
norm(Matrix(R2)-Matrix(R))/norm(Matrix(R))

0.0014230423742656228

Our computations were performed using 16 digits of accuracy, but the results seem to be accurate to only 2 or 3 digits. So, the individual rounding errors have been amplified by factors $\approx 10^{13}$...

However....

In [10]:
norm(A - Matrix(Q2)*Matrix(R2))/norm(A)

7.410078779986838e-16

An astonishing thing happens when we multiply the two inaccurate matrices Q2 and R2. The product $Q_2 R_2$ is accurate to a full 16 digits!

So, the huge errors in $Q_2$ and $R_2$ must be correlated as they cancel out when one takes the product.

Note that arbitrary such errors do not cancel out when we take a product.

Just take small random perturbations of $Q$ and $R$ (for which $A=QR$) and see what happens:

In [13]:
Q3 = Matrix(Q) + (1e-4)*rand(m,m)
R3 = Matrix(R) + (1e-4)*rand(m,m)

norm(A-Q3*R3)/norm(A)

0.0018034039691867726

Here $Q_3$ is no better or worse than $Q_2$, when compared to $Q$ and $R_3$ is no better or worse than $R_2$, when compared to $R$. But the product $Q_3 R_3$ is not accurate at all.

So, what's going on?

The errors in $Q_2$ and $R_2$ are forward errors.

In general, a large forward error can be the result of an ill-conditioned problem or an unstable algorithm (Theorem 15.1).

In our experiment, they are due to the former: an ill-conditiond problem. As a rule, the sequences of column spaces of a random triangular matrix are exceedingly ill-conditioned as a function of the entries of the matrix.

The error in $Q_2 R_2$ is the backward error or residual. The smallness of this error suggests that Householder triangularization is backward stable.

Recall the axioms:

A:

For all $x \in \mathbb{R}$, there exists $\epsilon$ with $|\epsilon| \leq \frac{1}{2}\epsilon_{\text{M}}$ such that $\operatorname{fl}(x)=x(1+\epsilon)$.

B:

For all $x, y \in \mathbb{F}$, there exists $\epsilon$ with $|\epsilon| \leq \frac{1}{2}\epsilon_{\text{M}}$ such that

$$ x \circledast y=(x * y)(1+\epsilon) . $$

Householder triangularization is backward stable for all matrices $A$ and all computers satisfying axioms A and B.

In a statement of the form $A=\tilde{Q}\tilde{R} + \delta A$, we let $\tilde{R}$ denote the computed $R$ matrix in the QR factorization with Householder triangularization.

The definition of $\tilde{Q}$ in the above statement is more subtle: we mean a certain matrix that is exactly unitary.

Recall that the exact $Q$ is equal to the product $Q_1 Q_2 \cdots Q_n$ where each $Q_k$ is the Householder reflector defined by the vector $v_k$, which is determined at the $k$th step of the algorithm.

In the floating point computation, we obtain instead a sequence of vectors $\tilde{v}_k$. We let $\tilde{Q}_k$ denote the exactly unitary reflector defined - mathematically, not on the computer - by the floating point vector $\tilde{v}_k$. Then define

$$ \tilde{Q}=\tilde{Q}_1 \tilde{Q}_2 \cdots \tilde{Q}_n $$

This exactly unitary matrix $\tilde{Q}$ will take the role of our "computed $Q$." This definition may seem odd at first, but it is the natural one.

In applications, as discussed in Lecture 10, the matrix $Q$ is generally not formed explicitly anyway, so it would not be useful to define a "computed $Q$ " of the more obvious kind.

It is the vectors $\tilde{v}_k$ that are formed explicitly, and these are what enter into (16.2).

Theorem:¶

Let the $Q R$ factorization $A=Q R$ of a matrix $A \in \mathbb{C}^{m \times n}$ be computed by Householder triangularization (Algorithm 10.1) on a computer satisfying the axioms A (13.5) and B (13.7), and let the computed factors $\tilde{Q}$ and $\tilde{R}$ as introduced above. Then we have

$$ \tilde{Q} \tilde{R}=A+\delta A, \qquad \frac{\|\delta A\|}{\|A\|}=O\left(\epsilon_{\text {M}}\right) $$

for some $\delta A\in\mathbb{C}^{m\times n}$.

Solving $Ax = b$ by QR Factorization (employing Householder Triangularization)¶

We have seen that Householder triangularization is backward stable but not always accurate in the forward sense.

The QR factorization in general is not an end itself --- it is means for computing other quantities: e.g., solving a system of linear equations, solving a least squares problem, computing eigenvalues, ...

Question: Is the backward stability of QR factorization enough to make it a satisfactory component of a larger algorithm?

We want to solve a nonsingular $m\times m$ system $Ax=b$. Recall the algorithm employing QR factorization:

Algorithm 16.1:

  • Compute $Q$, $R$ : $A = QR$.
  • Set $y= Q^* b$.
  • Solve the triangular system $Rx = y$ by backward substitution.

On the computer, the first step produces $\tilde{Q}$ and $\tilde{R}$ with $A \approx \tilde{Q} \tilde{R}$.

The next step produces $\tilde{y} \approx \tilde{Q}^* b$. This multiplication on the right-hand side introduces roundoff errors, so the computed vector on the computer is not exactly $\tilde{Q}^* b$, but some vector $\tilde{y}$.

But this vector $\tilde{y}$ satisfies the following backward stability estimate (since multiplication is backward stable):

$$ (\tilde{Q}+\delta Q) \tilde{y}=b, \quad\text{for some}\quad\|\delta Q\|=O(\epsilon_{\text {M}}). $$

(Note that the $2$-norm of $\tilde{Q}$ is $1$ since we are building the exactly unitary matrix from computed vectors.)

This exact equality tells us that the result of applying the Householder reflectors in floating point arithmetic is exactly equivalent to multiplying $b$ by a slightly perturbed matrix $(\tilde{Q}+\delta Q)^{-1}$.

In the final step, with the $\tilde{R}$ and $\tilde{y}$ available at hand (the exact $R$ and $y$ are long forgotten at this point), we solve $\tilde{R} x = \tilde{y}$ by backward substitution and compute a $\tilde{x}$. This step will also introduce roundoff errors, but it turns out that the computation is backward stable as well. We have

$$ (\tilde{R}+\delta R) \tilde{x}=\tilde{y}, \qquad \frac{\|\delta R\|}{\|\tilde{R}\|}=O(\epsilon_{\text{M}}). $$

As always the equality on the left is edxact. We will revisit this last stability result in the next lecture.

What happens up to this point is a typical backward stability result in numerical linear algebra.

Theorem 16.2¶

Algorithm 16.1 is bacward stable, satisfying

$$ (A+\Delta A) \tilde{x}=b, \quad \frac{\|\Delta A\|}{\|A\|}=O\left(\epsilon_{\text {M}}\right) $$

for some $\Delta A \in \mathbb{C}^{m\times n}$.

Proof:¶

We know that the computed solution $\tilde{x}$ satisfies

$$ (\tilde{R} + \delta R) \tilde{x} = \tilde{Q}^* \otimes b $$

for some $\delta R$ such that $\frac{\|\delta R\|}{\|\tilde{R}\|}=O(\epsilon_{\text{M}})$. The RHS satisfies $\tilde{Q}^* \otimes b = (\tilde{Q}^* + \delta Q)b= (\mathbb{I} + (\delta Q)\tilde{Q})\tilde{Q}^* b$ for some $\delta Q$. Thus, we obtain, for a different $\delta Q$ of the same size $\|\delta Q\| = O(\epsilon_\text{M})$,

$$ \tilde{Q}(\mathbb{I} + \delta Q)(\tilde{R}+\delta R) \tilde{x} = b. $$

Expanding yields

$$ \left[Q \tilde{R} + \tilde{Q} (\delta Q)\tilde{R} + \tilde{Q} (\delta R) + \tilde{Q} (\delta Q)(\delta R)\right]\tilde{x} = b, $$

which is

$$ \left[A + \delta A + \tilde{Q} (\delta Q)\tilde{R} + \tilde{Q} (\delta R) + \tilde{Q} (\delta Q)(\delta R)\right]\tilde{x} = b, $$

for some $\delta A$ satisfying $\| \delta A\|/\|A\| = O(\epsilon_\text{M})$. This equation is indeed of the form $(A+\Delta A )\tilde{x} = b$.

We need to show that the rest of the terms in $\Delta A$ also satisfy the same estimate:

$$ \frac{\|\tilde{Q} (\delta Q)\tilde{R} + \tilde{Q} (\delta R) + \tilde{Q} (\delta Q)(\delta R)\|}{\|A\|} = O(\epsilon_\text{M}). $$

First, since $\tilde{Q}$ is unitary, we have

$$ \|\tilde{R}\|=\|\tilde{Q}^*(A+\delta A)\| \leq\|\tilde{Q}^*\|(\|A\|+\|\delta A\|) \leq \|\tilde{Q}^*\| \|A\|(1 + O(\epsilon_\text{M})). $$

Then,

$$ \|\tilde{Q} (\delta Q) \tilde{R}\| \leq \|\tilde{Q}\| \| \delta Q \| \|\tilde{Q}^*\|\|A\|(1 + O(\epsilon_\text{M})) = \| A \| O(\epsilon_\text{M}) $$

since $\| \delta Q \|=O(\epsilon_\text{M})$.

Next, we observe:

$$ \| \tilde{Q} (\delta R)\| \leq \| \tilde{Q}\| \|\delta R\| \leq \|\tilde{Q}\| \| \tilde{R} \| O(\epsilon_\text{M}) = \|\tilde{Q}\| \| \tilde{Q}^*(A+\delta A) \| O(\epsilon_\text{M}) \leq \|\tilde{Q}\| \| \tilde{Q}^*\| \| A+\delta A \| O(\epsilon_\text{M}) \leq \|A\| O(\epsilon_\text{M}). $$

Finally,

$$ \| \tilde{Q} (\delta Q) (\delta R) \| \leq \| \tilde{Q} \| \| \delta Q\| \| \delta R \| \leq \| \tilde{Q} \| O(\epsilon_\text{M}) O(\epsilon_\text{M})\|\tilde{R}\| \leq \| \tilde{Q} \| O(\epsilon_\text{M}) O(\epsilon_\text{M})\|\tilde{Q}^*\|\|A+\delta A\| \leq O(\epsilon_\text{M}^2) \|A\| + O(\epsilon_\text{M}^3) \|A\| = O(\epsilon_\text{M}^2) \|A\|. $$

Combining these estimates yields

$$ \frac{\|\Delta A\|}{\|A\|} \leq \frac{\|\delta A\|}{\|A\|} + \frac{\|\tilde{Q} (\delta Q) \tilde{R}\|}{\|A\|} + \frac{\| \tilde{Q} (\delta R)\|}{\|A\|} + \frac{\| \tilde{Q} (\delta Q) (\delta R) \|}{\| A\|} = O(\epsilon_\text{M}), $$

which establishes bacward stability.

Now, let's recall the condition number for the problem of computing $A^{-1} b$ where the input is $A$. The problem is $f(A)=A^{-1} b$. Recall that the (relative) condition number for this problem's instance is

$$ \kappa(A) = \|A\| \|A^{-1}\|. $$

Thus, we have the following accuracy result.

Theorem:¶

The solution $\tilde{x}$ of $Ax=b$ computed by Algorithm 16.1 satisfies

$$ \frac{\|\tilde{x}-x\|}{\|x\|}=O(\kappa(A) \epsilon_{\text {M}}). $$