The QR factorization - Fundamentals of Numerical Computation

Sets of vectors satisfying a certain property are useful both theoretically and computationally.

In two and three dimensions, orthogonality is the same as perpendicularity.

Orthogonal vectors simplify inner products. For example, if $\mathbf{q}_1$ and $\mathbf{q}_2$ are orthogonal, then

\| \mathbf{q}_1 - \mathbf{q}_2 \|_2^2 = (\mathbf{q}_1-\mathbf{q}_2)^T(\mathbf{q}_1-\mathbf{q}_2) = \mathbf{q}_1^T\mathbf{q}_1 - 2 \mathbf{q}_1^T\mathbf{q}_2 + \mathbf{q}_2^T\mathbf{q}_2 = \|\mathbf{q}_1\|_2^2 + \|\mathbf{q}_2\|_2^2.

(3.3.2)

As in the rest of this chapter, we will be using the 2-norm exclusively.

Equation (3.3.2) is the key to the computational attractiveness of orthogonality. Figure 3.3.1 shows how nonorthogonal vectors can allow a multidimensional version of subtractive cancellation, in which $\|\mathbf{x}-\mathbf{y}\|$ is much smaller than $\|\mathbf{x}\|$ and $\|\mathbf{y}\|$ . As the figure illustrates, orthogonal vectors do not allow this phenomenon. By (3.3.2), the magnitude of a vector difference or sum is larger than the magnitudes of the original vectors.

Figure 3.3.1:Nonorthogonal vectors can cause cancellation when subtracted, but orthogonal vectors never do.

3.3.1Orthogonal and ONC matrices¶

Statements about orthogonal vectors are often made more easily in matrix form. Let $\mathbf{Q}$ be an $n\times k$ matrix whose columns $\mathbf{q}_1, \ldots, \mathbf{q}_k$ are orthogonal vectors. The orthogonality conditions (3.3.1) become simply that $\mathbf{Q}^T\mathbf{Q}$ is a diagonal matrix, since

\mathbf{Q}^T \mathbf{Q} = \begin{bmatrix} \mathbf{q}_1^T \\[1mm] \mathbf{q}_2^T \\ \vdots \\ \mathbf{q}_k^T \end{bmatrix} \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_k \end{bmatrix} = \begin{bmatrix} \mathbf{q}_1^T\mathbf{q}_1 & \mathbf{q}_1^T\mathbf{q}_2 & \cdots & \mathbf{q}_1^T\mathbf{q}_k \\[1mm] \mathbf{q}_2^T\mathbf{q}_1 & \mathbf{q}_2^T\mathbf{q}_2 & \cdots & \mathbf{q}_2^T\mathbf{q}_k \\ \vdots & \vdots & & \vdots \\ \mathbf{q}_k^T\mathbf{q}_1 & \mathbf{q}_k^T\mathbf{q}_2 & \cdots & \mathbf{q}_k^T\mathbf{q}_k \end{bmatrix}.

(3.3.3)

If the columns of $\mathbf{Q}$ are orthonormal, then $\mathbf{Q}^T\mathbf{Q}$ is the $k\times k$ identity matrix. This is such an important property that we will break with common practice here and give this type of matrix a name.

Of particular interest is a square ONC matrix.^[1]

Orthogonal matrices have properties beyond Theorem 3.3.1.

3.3.2Orthogonal factorization¶

Now we come to another important way to factor a matrix: the QR factorization. As we will show below, the QR factorization plays a role in linear least squares analogous to the role of LU factorization in linear systems.

In most introductory books on linear algebra, the QR factorization is derived through a process known as Gram–Schmidt orthogonalization. However, while it is an important tool for theoretical work, the Gram–Schmidt process is numerically unstable. We will consider an alternative construction in Computing QR factorizations.

When $m$ is much larger than $n$ , which is often the case, there is a compressed form of the factorization that is more efficient. In the product

\mathbf{A} = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_m \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ 0 & r_{22} & \cdots & r_{2n} \\ \vdots & & \ddots & \vdots\\ 0 & 0 & \cdots & r_{nn} \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix},

(3.3.5)

the vectors $\mathbf{q}_{n+1},\ldots,\mathbf{q}_m$ always get multiplied by zero. Nothing about $\mathbf{A}$ is lost if we delete them and reduce the factorization to the equivalent form

\mathbf{A} = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_n \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ 0 & r_{22} & \cdots & r_{2n} \\ \vdots & & \ddots & \vdots\\ 0 & 0 & \cdots & r_{nn} \end{bmatrix} = \hat{\mathbf{Q}} \hat{\mathbf{R}}.

(3.3.6)

Example 3.3.1 (QR factorization)

Julia

MATLAB

Python

Example 3.3.1

Julia provides access to both the thin and full forms of the QR factorization.

A = rand(1.:9., 6, 4)
@show m,n = size(A);

(m, n) = size(A) = (6, 4)

Here is a standard call:

Q,R = qr(A);
Q

6×6 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}

4×4 Matrix{Float64}:
 -17.0587  -12.838    -11.9001   -12.1346
   0.0       8.89863    2.16065    0.36146
   0.0       0.0        4.44071   -2.5183
   0.0       0.0        0.0       -2.29784

If you look carefully, you see that we seemingly got a full $\mathbf{Q}$ but a thin $\mathbf{R}$ . However, the $\mathbf{Q}$ above is not a standard matrix type. If you convert it to a true matrix, then it reverts to the thin form.

Q̂ = Matrix(Q)

6×4 Matrix{Float64}:
 -0.468968   0.334814   -0.0684959  -0.442051
 -0.351726   0.0544507   0.382097    0.576839
 -0.527589  -0.311643   -0.361428   -0.348326
 -0.117242   0.617494    0.51132    -0.279296
 -0.527589  -0.42402     0.368817    0.139263
 -0.293105   0.476154   -0.56675     0.503106

We can test that $\mathbf{Q}$ is an orthogonal matrix:

opnorm(Q' * Q - I)

8.331709950177064e-16

The thin $\hat{\mathbf{Q}}$ cannot be an orthogonal matrix, because it is not square, but it is still ONC:

Q̂' * Q̂ - I

4×4 Matrix{Float64}:
 -1.11022e-16   9.51871e-17  -2.86165e-17  5.50694e-17
  9.51871e-17   0.0          -2.55659e-17  3.04621e-17
 -2.86165e-17  -2.55659e-17   0.0          1.13916e-16
  5.50694e-17   3.04621e-17   1.13916e-16  2.22045e-16

Example 3.3.1

MATLAB provides access to both the thin and full forms of the QR factorization.

A = magic(5);
A = A(:, 1:4);
[m, n] = size(A)

Here is the full form:

[Q, R] = qr(A);
szQ = size(Q), szR = size(R)

We can test that $\mathbf{Q}$ is an orthogonal matrix:

QTQ = Q' * Q
norm(QTQ - eye(m))

With a second input argument given to qr, the thin form is returned. (This is usually the one we want in practice.)

[Q_hat, R_hat] = qr(A, 0);
szQ_hat = size(Q_hat), szR_hat = size(R_hat)

Now $\hat{\mathbf{Q}}$ cannot be an orthogonal matrix, because it is not square, but it is still ONC. Mathematically, $\hat{\mathbf{Q}}^T \hat{\mathbf{Q}}$ is a $4\times 4$ identity matrix.

Q_hat' * Q_hat - eye(n)

Example 3.3.1

MATLAB provides access to both the thin and full forms of the QR factorization.

A = 1.0 + floor(9 * random.rand(6,4))
A.shape

(6, 4)

Here is the full form:

from numpy.linalg import qr
Q, R = qr(A, "complete")
print(f"size of Q is {Q.shape}")
print("R:")
print(R)

size of Q is (6, 6)
R:
[[-14.76482306 -10.83656738 -12.39432395 -13.88435196]
 [  0.         -10.12762595  -5.20241107  -4.89169769]
 [  0.           0.          -7.63646862  -5.03691218]
 [  0.           0.           0.          -3.99068666]
 [  0.           0.           0.           0.        ]
 [  0.           0.           0.           0.        ]]

We can test that $\mathbf{Q}$ is an orthogonal matrix:

print(f"norm of (Q^T Q - I) is {norm(Q.T @ Q - eye(6)):.3e}")

norm of (Q^T Q - I) is 7.294e-16

The default for qr, and the one you usually want, is the thin form.

Q_hat, R_hat = qr(A)
print(f"size of Q_hat is {Q_hat.shape}")
print("R_hat:")
print(R_hat)

size of Q_hat is (6, 4)
R_hat:
[[-14.76482306 -10.83656738 -12.39432395 -13.88435196]
 [  0.         -10.12762595  -5.20241107  -4.89169769]
 [  0.           0.          -7.63646862  -5.03691218]
 [  0.           0.           0.          -3.99068666]]

Now $\hat{\mathbf{Q}}$ cannot be an orthogonal matrix, because it is not square, but it is still ONC. Mathematically, $\hat{\mathbf{Q}}^T \hat{\mathbf{Q}}$ is a $4\times 4$ identity matrix.

print(f"norm of (Q_hat^T Q_hat - I) is {norm(Q_hat.T @ Q_hat - eye(4)):.3e}")

norm of (Q_hat^T Q_hat - I) is 6.135e-16

3.3.3Least squares and QR¶

If we substitute the thin factorization (3.3.6) into the normal equations (3.2.3), we can simplify expressions a great deal.

\begin{split} \mathbf{A}^T\mathbf{A} \mathbf{x} &= \mathbf{A}^T \mathbf{b}, \\ \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \hat{\mathbf{Q}} \hat{\mathbf{R}} \mathbf{x} &= \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \mathbf{b}, \\ \hat{\mathbf{R}}^T \hat{\mathbf{R}} \mathbf{x}& = \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \mathbf{b}. \end{split}

(3.3.7)

In order to have the normal equations be well posed, we require that $\mathbf{A}$ is not rank-deficient (as proved in Theorem 3.2.2). This is enough to guarantee that $\hat{\mathbf{R}}$ is nonsingular (see Exercise 3.3.4). Therefore, its transpose is nonsingular as well, and we arrive at

\hat{\mathbf{R}} \mathbf{x}=\hat{\mathbf{Q}}^T \mathbf{b}.

(3.3.8)

This algorithm is implemented in Function 3.3.2. It is essentially the algorithm used internally by Julia when A\b is called. The execution time is dominated by the factorization, the most common method for which is described in Computing QR factorizations.

Algorithm 3.3.2 (lsqrfact)

Julia

MATLAB

Python

Solution of least squares by QR factorization

lsqrfact.jl

1
2
3
4
5
6
7
8
9
10
11
12
"""
    lsqrfact(A, b)

Solve a linear least-squares problem by QR factorization. Returns
the minimizer of ||b-Ax||.
"""
function lsqrfact(A, b)
    Q, R = qr(A)
    c = Q' * b
    x = backsub(R, c)
    return x
end

Solution of least squares by QR factorization

lsqrfact.m

1
2
3
4
5
6
7
8
9
10
11
12
function x = lsqrfact(A, b)
% LSQRFACT   Solve linear least squares by QR factorization.
% Input: 
%   A     coefficient matrix (m by n, m > n)
%   b     right-hand side (m by 1)
% Output:
%   x     minimizer of || b - Ax ||

    [Q, R] = qr(A, 0);    % thin factorization
    c = Q' * b;
    x = backsub(R, c);
end

Solution of least squares by QR factorization

lsqrfact.py

1
2
3
4
5
6
7
8
9
10
11
def lsqrfact(A, b):
    """
    lsqrfact(A, b)
    
    Solve a linear least squares problem by QR factorization. Returns the
    minimizer of ||b-Ax||.
    """
    Q, R = np.linalg.qr(A)
    c = Q.T @ b
    x = backsub(R, c)
    return x

The solution of least-squares problems via QR factorization is more stable than when the normal equations are formulated and solved directly.

Example 3.3.2 (Stability of least-squares via QR)

Julia

MATLAB

Python

Example 3.3.2

We’ll repeat the experiment of Example 3.2.1, which exposed instability in the normal equations.

t = range(0, 3, 400)
f = [ x -> sin(x)^2, x -> cos((1 + 1e-7) * x)^2, x -> 1. ]
A = [ f(t) for t in t, f in f ]
x = [1., 2, 1]
b = A * x;

The error in the solution by Function 3.3.2 is similar to the bound predicted by the condition number.

observed_error = norm(FNC.lsqrfact(A, b) - x) / norm(x);
@show observed_error;
@show error_bound = cond(A) * eps();

observed_error = 4.665273501889628e-9


error_bound = cond(A) * eps() = 4.053030228488391e-9

Example 3.3.2

We’ll repeat the experiment of Example 3.2.1, which exposed instability in the normal equations.

t = linspace(0, 3, 400)';
A = [ sin(t).^2, cos((1+1e-7)*t).^2, t.^0 ];
x = [1; 2; 1];
b = A * x;

The error in the solution by Function 3.3.2 is similar to the bound predicted by the condition number.

observed_error = norm(lsqrfact(A, b) - x) / norm(x)
error_bound = cond(A) * eps

Unrecognized function or variable 'lsqrfact'.

Example 3.3.2

We’ll repeat the experiment of Example 3.2.1, which exposed instability in the normal equations.

t = linspace(0, 3, 400)
A = array([ [sin(t)**2, cos((1+1e-7)*t)**2, 1] for t in t ])
x = array([1, 2, 1])
b = A @ x

The error in the solution by Function 3.3.2 is similar to the bound predicted by the condition number.

print(f"observed error: {norm(FNC.lsqrfact(A, b) - x) / norm(x):.3e}")
print(f"conditioning bound: {cond(A) * finfo(float).eps:.3e}")

observed error: 5.186e-09
conditioning bound: 4.053e-09

3.3.4Exercises¶

Footnotes¶

Confusingly, a square matrix whose columns are orthogonal is not necessarily an orthogonal matrix; the columns must be orthonormal, which is a stricter condition.
↩

Preface

The normal equations

Preface

Computing QR factorizations