QR factorization and least squares

Section 4.4 \(QR\) factorization and least squares

Subsection From Gram–Schmidt to \(QR\)

The Gram–Schmidt algorithm 4.2.17 processes the independent columns of a matrix from left to right. It constructs orthonormal directions \(\mathbf{q}_1,\ldots,\mathbf{q}_n\) and records coefficients \(r_{ij}\text{.}\) These same quantities form a matrix factorization.

🔗

Let

\begin{equation*} A = \begin{bmatrix} \mathbf{a}_1\amp\mathbf{a}_2\amp\cdots\amp\mathbf{a}_n \end{bmatrix}. \end{equation*}

🔗

At column \(j\text{,}\) the Gram–Schmidt algorithm 4.2.17 computes

\begin{equation*} \mathbf{p}_j = \sum_{i\lt j} r_{ij}\mathbf{q}_i, \qquad \mathbf{v}_j = \mathbf{a}_j-\mathbf{p}_j, \end{equation*}

\begin{equation*} r_{jj} = \|\mathbf{v}_j\|, \qquad \mathbf{q}_j = \frac{\mathbf{v}_j}{r_{jj}}. \end{equation*}

🔗

Since

\begin{equation*} \mathbf{v}_j = r_{jj}\mathbf{q}_j, \end{equation*}

we can reconstruct the original column:

\begin{equation*} \mathbf{a}_j = \mathbf{p}_j+\mathbf{v}_j = \sum_{i\lt j} r_{ij}\mathbf{q}_i + r_{jj}\mathbf{q}_j = \sum_{i\leq j} r_{ij}\mathbf{q}_i. \end{equation*}

🔗

Thus column \(j\) uses only \(\mathbf{q}_1,\ldots,\mathbf{q}_j\text{.}\)

🔗

Collect the orthonormal directions in

\begin{equation*} Q = \begin{bmatrix} \mathbf{q}_1\amp\mathbf{q}_2\amp\cdots\amp\mathbf{q}_n \end{bmatrix} \end{equation*}

and the coefficients in

\begin{equation*} R = \begin{bmatrix} r_{11}\amp r_{12}\amp\cdots\amp r_{1n}\\ 0\amp r_{22}\amp\cdots\amp r_{2n}\\ \vdots\amp\ddots\amp\ddots\amp\vdots\\ 0\amp\cdots\amp0\amp r_{nn} \end{bmatrix}. \end{equation*}

🔗

The column equations above become the single matrix equation

\begin{equation*} A=QR. \end{equation*}

🔗

The zeros below the diagonal of \(R\) record the order of the algorithm: column \(j\) uses only directions constructed by stage \(j\text{.}\)

🔗

Definition 4.4.1. Reduced \(QR\) factorization.

Let \(A\) be an \(m\times n\) matrix with linearly independent columns. A reduced \(QR\) factorization of \(A\) is a factorization

\begin{equation*} A=QR, \end{equation*}

where

🔗

\(Q\) is an \(m\times n\) matrix with orthonormal columns;
🔗

🔗
\(R\) is an invertible upper-triangular \(n\times n\) matrix.
🔗

🔗

For the Gram–Schmidt convention used in this unit, the diagonal entries of \(R\) are positive:

\begin{equation*} r_{jj}=\|\mathbf{v}_j\|>0. \end{equation*}

🔗

Activity 4.4.2. Reading \(Q^T Q\) (U4-LO1).

Let

\begin{equation*} Q = \begin{bmatrix} 1/\sqrt{3}\amp1/\sqrt{6}\\ 1/\sqrt{3}\amp-2/\sqrt{6}\\ 1/\sqrt{3}\amp1/\sqrt{6} \end{bmatrix} = \begin{bmatrix} \mathbf{q}_1\amp\mathbf{q}_2 \end{bmatrix}. \end{equation*}

🔗

Compute \(\|\mathbf{q}_1\|^2\) and \(\|\mathbf{q}_2\|^2\text{.}\)
🔗

🔗
Compute \(\mathbf{q}_1^T\mathbf{q}_2\text{.}\)
🔗

🔗
Compute \(Q^TQ\text{.}\)
🔗

🔗
What does the \((i,j)\)-entry of \(Q^TQ\) measure?
🔗

🔗
Explain why

\begin{equation*} Q^TQ=I_2 \end{equation*}

is the matrix check that the columns of \(Q\) are orthonormal.

🔗

🔗

Solution.

We have

\begin{equation*} \|\mathbf{q}_1\|^2 = \frac13+\frac13+\frac13 = 1 \end{equation*}

and

\begin{equation*} \|\mathbf{q}_2\|^2 = \frac16+\frac46+\frac16 = 1. \end{equation*}

🔗

Also,

\begin{equation*} \mathbf{q}_1^T\mathbf{q}_2 = \frac{1}{\sqrt{18}} - \frac{2}{\sqrt{18}} + \frac{1}{\sqrt{18}} = 0. \end{equation*}

🔗

Therefore

\begin{equation*} Q^TQ = \begin{bmatrix} \mathbf{q}_1^T\mathbf{q}_1\amp \mathbf{q}_1^T\mathbf{q}_2\\ \mathbf{q}_2^T\mathbf{q}_1\amp \mathbf{q}_2^T\mathbf{q}_2 \end{bmatrix} = \begin{bmatrix} 1\amp0\\ 0\amp1 \end{bmatrix} = I_2. \end{equation*}

🔗

The \((i,j)\)-entry of \(Q^TQ\) is the dot product \(\mathbf{q}_i^T\mathbf{q}_j\text{.}\) The diagonal entries check unit length, and the off-diagonal entries check orthogonality.

🔗

Theorem 4.4.3. Gram–Schmidt produces \(QR\).

Every \(m\times n\) matrix \(A\) with linearly independent columns has a reduced \(QR\) factorization

\begin{equation*} A=QR \end{equation*}

in which the columns of \(Q\) are orthonormal and \(R\) is upper triangular with positive diagonal entries.

🔗

Why is this true?.

Apply the Gram–Schmidt algorithm 4.2.17 to the columns

\begin{equation*} \mathbf{a}_1,\ldots,\mathbf{a}_n. \end{equation*}

🔗

At stage \(j\text{,}\) the algorithm gives

\begin{equation*} \mathbf{a}_j = \sum_{i\leq j} r_{ij}\mathbf{q}_i. \end{equation*}

🔗

Collecting these equations column by column gives

\begin{equation*} A=QR. \end{equation*}

🔗

The columns of \(Q\) are orthonormal. The condition \(i\leq j\) makes \(R\) upper triangular, and

\begin{equation*} r_{jj} = \|\mathbf{v}_j\| > 0 \end{equation*}

because the input columns are linearly independent.

🔗

Corollary 4.4.4. Square case.

Every invertible square matrix has a \(QR\) factorization

\begin{equation*} A=QR, \end{equation*}

where \(Q\) is an orthogonal matrix and \(R\) is upper triangular with positive diagonal entries.

🔗

Why is this true?.

An invertible square matrix has linearly independent columns. Apply “Gram–Schmidt produces QR 4.4.3.” Since \(Q\) is square and has orthonormal columns,

\begin{equation*} Q^TQ=QQ^T=I. \end{equation*}

Thus \(Q\) is an orthogonal matrix.

🔗

Activity 4.4.5. Build \(Q\) and \(R\) from the Gram–Schmidt ledger (U4-LO2, U4-LO6).

Let

\begin{equation*} A = \begin{bmatrix} 3\amp1\\ 4\amp2 \end{bmatrix} = \begin{bmatrix} \mathbf{a}_1\amp\mathbf{a}_2 \end{bmatrix}. \end{equation*}

🔗

Compute

\begin{equation*} r_{11}=\|\mathbf{a}_1\| \end{equation*}

and

\begin{equation*} \mathbf{q}_1 = \frac{\mathbf{a}_1}{r_{11}}. \end{equation*}

🔗

🔗
Compute the projection coefficient

\begin{equation*} r_{12} = \mathbf{q}_1^T\mathbf{a}_2. \end{equation*}

🔗

🔗
Compute

\begin{equation*} \mathbf{p}_2 = r_{12}\mathbf{q}_1 \end{equation*}

and

\begin{equation*} \mathbf{v}_2 = \mathbf{a}_2-\mathbf{p}_2. \end{equation*}

🔗

🔗
Compute

\begin{equation*} r_{22} = \|\mathbf{v}_2\| \end{equation*}

and

\begin{equation*} \mathbf{q}_2 = \frac{\mathbf{v}_2}{r_{22}}. \end{equation*}

🔗

🔗
Assemble

\begin{equation*} Q = \begin{bmatrix} \mathbf{q}_1\amp\mathbf{q}_2 \end{bmatrix} \end{equation*}

and

\begin{equation*} R = \begin{bmatrix} r_{11}\amp r_{12}\\ 0\amp r_{22} \end{bmatrix}. \end{equation*}

🔗

🔗
Check

\begin{equation*} Q^TQ=I_2 \end{equation*}

and

\begin{equation*} QR=A. \end{equation*}

🔗

🔗

Solution.

The first column gives

\begin{equation*} r_{11} = \|\mathbf{a}_1\| = 5, \qquad \mathbf{q}_1 = \frac15 \begin{bmatrix} 3\\ 4 \end{bmatrix}. \end{equation*}

🔗

The projection coefficient is

\begin{equation*} r_{12} = \mathbf{q}_1^T\mathbf{a}_2 = \frac35(1)+\frac45(2) = \frac{11}{5}. \end{equation*}

🔗

Therefore

\begin{equation*} \mathbf{p}_2 = r_{12}\mathbf{q}_1 = \frac{11}{5} \begin{bmatrix} 3/5\\ 4/5 \end{bmatrix} = \begin{bmatrix} 33/25\\ 44/25 \end{bmatrix}, \end{equation*}

and

\begin{equation*} \mathbf{v}_2 = \mathbf{a}_2-\mathbf{p}_2 = \begin{bmatrix} 1\\ 2 \end{bmatrix} - \begin{bmatrix} 33/25\\ 44/25 \end{bmatrix} = \begin{bmatrix} -8/25\\ 6/25 \end{bmatrix}. \end{equation*}

🔗

Its length is

\begin{equation*} r_{22} = \|\mathbf{v}_2\| = \frac{2}{5}, \end{equation*}

\begin{equation*} \mathbf{q}_2 = \frac{\mathbf{v}_2}{r_{22}} = \begin{bmatrix} -4/5\\ 3/5 \end{bmatrix}. \end{equation*}

🔗

Thus

\begin{equation*} Q = \begin{bmatrix} 3/5\amp-4/5\\ 4/5\amp3/5 \end{bmatrix}, \qquad R = \begin{bmatrix} 5\amp11/5\\ 0\amp2/5 \end{bmatrix}. \end{equation*}

🔗

The orthonormality check gives

\begin{equation*} Q^TQ = I_2, \end{equation*}

and direct multiplication gives

\begin{equation*} QR = \begin{bmatrix} 3\amp1\\ 4\amp2 \end{bmatrix} = A. \end{equation*}

🔗

Activity 4.4.6. An upper-triangular matrix already has a \(QR\) factorization (U4-LO6).

Let

\begin{equation*} A = \begin{bmatrix} 6\amp2\\ 0\amp3 \end{bmatrix}. \end{equation*}

🔗

Explain why

\begin{equation*} Q=I_2, \qquad R=A \end{equation*}

is a \(QR\) factorization. Check each condition in the definition.

🔗

Solution.

The columns of \(I_2\) are orthonormal, so

\begin{equation*} Q^TQ=I_2. \end{equation*}

🔗

The matrix

\begin{equation*} R = \begin{bmatrix} 6\amp2\\ 0\amp3 \end{bmatrix} \end{equation*}

is upper triangular and has positive diagonal entries. Finally,

\begin{equation*} QR = I_2A = A. \end{equation*}

🔗

Thus \(A=QR\) is a \(QR\) factorization.

🔗

Warning 4.4.7. Rectangular \(Q\).

In a reduced \(QR\) factorization of an \(m\times n\) matrix with \(m\gt n\text{,}\) the matrix \(Q\) has orthonormal columns, so

\begin{equation*} Q^TQ=I_n. \end{equation*}

🔗

However, \(Q\) is not square, and generally

\begin{equation*} QQ^T\neq I_m. \end{equation*}

🔗

Instead,

\begin{equation*} QQ^T\mathbf{b} = \operatorname{proj}_{\operatorname{col}(Q)}(\mathbf{b}) = \operatorname{proj}_{\operatorname{col}(A)}(\mathbf{b}). \end{equation*}

🔗

Warning 4.4.8. Signs in a \(QR\) factorization.

The positive-diagonal convention fixes the signs of the Gram–Schmidt factorization. Without that convention, one may replace a column \(\mathbf{q}_j\) by \(-\mathbf{q}_j\) and replace row \(j\) of \(R\) by its negative. The product \(QR\) does not change.

🔗

Numerical software may therefore return \(Q\) and \(R\) with signs different from a hand calculation. Check

\begin{equation*} Q^TQ=I \qquad\text{and}\qquad QR=A \end{equation*}

rather than requiring every entry to have the same sign.

🔗

Subsection Least squares with \(QR\)

The normal equations 4.3.3 give a direct route from a \(QR\) factorization to a least-squares solution.

🔗

Suppose

\begin{equation*} A=QR \end{equation*}

is a reduced \(QR\) factorization. The normal equations 4.3.3 are

\begin{equation*} A^TA\widehat{\mathbf{x}} = A^T\mathbf{b}. \end{equation*}

🔗

Substituting \(A=QR\) gives

\begin{equation*} R^TQ^TQR\widehat{\mathbf{x}} = R^TQ^T\mathbf{b}. \end{equation*}

🔗

Since

\begin{equation*} Q^TQ=I, \end{equation*}

this becomes

\begin{equation*} R^TR\widehat{\mathbf{x}} = R^TQ^T\mathbf{b}. \end{equation*}

🔗

The matrix \(R\) is invertible, so \(R^T\) is invertible. Canceling \(R^T\) gives

\begin{equation*} R\widehat{\mathbf{x}} = Q^T\mathbf{b}. \end{equation*}

🔗

Algorithm 4.4.9. Least squares with \(QR\).

Input. An \(m\times n\) matrix \(A\) with linearly independent columns and a vector \(\mathbf{b}\in\mathbb{R}^m\text{.}\)

🔗

Compute a reduced \(QR\) factorization

\begin{equation*} A=QR. \end{equation*}

🔗

🔗
Compute

\begin{equation*} \mathbf{y}=Q^T\mathbf{b}. \end{equation*}

🔗

🔗
Solve the upper-triangular system

\begin{equation*} R\widehat{\mathbf{x}}=\mathbf{y} \end{equation*}

by back-substitution.

🔗

🔗

Output. The vector \(\widehat{\mathbf{x}}\) is the least-squares coefficient vector. The fitted vector is

\begin{equation*} A\widehat{\mathbf{x}} = QR\widehat{\mathbf{x}} = QQ^T\mathbf{b}. \end{equation*}

🔗

Activity 4.4.10. Back-substitution (U2-LO2).

Solve

\begin{equation*} R\mathbf{x} = \mathbf{y}, \end{equation*}

where

\begin{equation*} R = \begin{bmatrix} 4\amp2\amp7\\ 0\amp3\amp5\\ 0\amp0\amp2 \end{bmatrix}, \qquad \mathbf{y} = \begin{bmatrix} 29\\ 21\\ 6 \end{bmatrix}. \end{equation*}

🔗

Work from the last equation upward.

🔗

Solution.

The last equation gives

\begin{equation*} 2x_3=6, \qquad x_3=3. \end{equation*}

🔗

The second equation gives

\begin{equation*} 3x_2+5x_3=21, \end{equation*}

\begin{equation*} 3x_2+15=21, \qquad x_2=2. \end{equation*}

🔗

The first equation gives

\begin{equation*} 4x_1+2x_2+7x_3=29, \end{equation*}

\begin{equation*} 4x_1+4+21=29, \qquad x_1=1. \end{equation*}

🔗

Therefore

\begin{equation*} \mathbf{x} = \begin{bmatrix} 1\\ 2\\ 3 \end{bmatrix}. \end{equation*}

🔗

Activity 4.4.11. Solving least squares with \(QR\) in code (U4-LO6, U4-LO7).

The following code solves the regression problem from “Fit a line to three data points 4.3.8” using \(QR\text{.}\)

🔗

import numpy as np

t = np.array([0.0, 1.0, 2.0])
y = np.array([1.0, 2.0, 2.0])

X = np.column_stack([np.ones_like(t), t])

Q, R = np.linalg.qr(X, mode="reduced")

c_qr = np.linalg.solve(R, Q.T @ y)
c_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]

(
    Q.shape,
    R.shape,
    np.allclose(Q.T @ Q, np.eye(Q.shape[1])),
    np.allclose(Q @ R, X),
    c_qr,
    c_lstsq,
    np.allclose(c_qr, c_lstsq),
)

What does the assignment
🔗
```
Q, R = np.linalg.qr(X, mode="reduced")
```
store in Q and R?
🔗

🔗
Why does the code use mode="reduced"?
🔗

🔗
What should Q.shape and R.shape be?
🔗

🔗
What does
🔗
```
Q.T @ Q
```
check?
🔗

🔗
What does
🔗
```
Q @ R
```
check?
🔗

🔗
What array is constructed by
🔗
```
np.eye(Q.shape[1])
```
in this example?
🔗

🔗
What mathematical vector is computed by
🔗
```
Q.T @ y
```
?
🔗

🔗
Why does the code solve a system with R?
🔗

🔗
Why is
🔗
```
np.linalg.solve(R, Q.T @ y)
```
preferable to explicitly computing an inverse of R?
🔗

🔗
Why should c_qr and c_lstsq agree?
🔗

🔗
Does their agreement mean that

\begin{equation*} X\mathbf{c}=\mathbf{y} \end{equation*}

has an exact solution?

🔗

🔗
Why may the entries of Q and R have signs different from a hand Gram–Schmidt calculation?
🔗

🔗

Solution.

The call

🔗

Q, R = np.linalg.qr(X, mode="reduced")

computes a reduced \(QR\) factorization and stores the two returned arrays in Q and R.

🔗

Since X has shape (3, 2), reduced \(QR\) gives

🔗

Q.shape == (3, 2)
R.shape == (2, 2)

The expression

🔗

Q.T @ Q

checks whether the columns of Q are orthonormal. The expression

🔗

Q @ R

checks whether the factors reconstruct X.

🔗

Here,

🔗

np.eye(Q.shape[1])

constructs the \(2\times2\) identity matrix.

🔗

The vector

🔗

Q.T @ y

is the right-hand side \(Q^T\mathbf{y}\) in the reduced system

\begin{equation*} R\widehat{\mathbf{c}} = Q^T\mathbf{y}. \end{equation*}

🔗

The matrix R is square, invertible, and upper triangular, so the code solves this system with np.linalg.solve. It does not need to construct an inverse matrix.

🔗

Both methods compute the same least-squares coefficient vector, so

🔗

c_qr

and

🔗

c_lstsq

should both be close to

🔗

array([1.16666667, 0.5])

and the final np.allclose check should return True.

🔗

This agreement does not mean that

\begin{equation*} X\mathbf{c}=\mathbf{y} \end{equation*}

has an exact solution. The residual from the regression activity 4.3.8 is nonzero.

🔗

A valid \(QR\) factorization can use simultaneous sign changes in one column of \(Q\) and the corresponding row of \(R\text{.}\) Therefore numerical output may have different signs from a hand calculation while still satisfying

\begin{equation*} Q^TQ=I \qquad\text{and}\qquad QR=X. \end{equation*}

🔗

The normal equations 4.3.3 and \(QR\) compute the same least-squares solution. The normal equations expose residual orthogonality directly. \(QR\) reuses the orthonormal directions produced by Gram–Schmidt 4.2.17 and reduces the computation to an upper-triangular solve.

🔗

Prev Top Next