Least Squares with QR, linear algebra I lesson

Can the least-squares line for the points $(0, 1)$ , $(1, 2)$ and $(2, 2)$ be found without forming the matrix $X^{\top}X$ of the normal equations? A reduced QR factorization of the design matrix gives the model’s column space perpendicular unit directions, and the least-squares problem then shrinks to one triangular system solved from the bottom up. This lesson derives that system, carries out a complete fit with it, and keeps track of the part of the data the triangular system never sees, the residual that the model cannot explain.

Key terms

Least squares by QR factorization

Definition Let $X$ be a real $m \times n$ matrix with linearly independent columns, so $m \ge n$ , and let $\mathbf{y}$ be a vector in $\mathbb{R}^m$ . The reduced QR factorization $X = QR$ has an $m \times n$ matrix $Q$ with orthonormal columns, so $Q^{\top}Q = I_n$ , whose columns span $\operatorname{Col} X$ , and an $n \times n$ upper triangular matrix $R$ with positive diagonal entries, which is therefore invertible. The method has three steps. Compute $Q$ and $R$ , for example by Gram-Schmidt. Compute the vector $Q^{\top}\mathbf{y}$ . Solve the triangular system $R\mathbf{c} = Q^{\top}\mathbf{y}$ by back substitution. Its solution is the unique least-squares solution $\hat{\mathbf{c}}$ of $X\mathbf{c} = \mathbf{y}$ . On this page a tuple written for a vector stands for a column.

In plain words QR replaces the possibly slanted directions of the design columns by perpendicular unit directions that span the same subspace. The data and the model stay the same, and only the coordinates used inside the column space change. In exact arithmetic this route gives the same answer as the normal equations. In floating point, forming $X^{\top}X$ can make rounding errors more harmful, and a QR route avoids forming it, although a small exact example cannot show which route is more accurate.

Example Fit the line $a + bt$ to the points $(0, 1)$ , $(1, 2)$ and $(2, 2)$ . The design columns are $(1, 1, 1)$ and $(0, 1, 2)$ . The second column has dot product $3/\sqrt{3} = \sqrt{3}$ with the first unit vector, its projection is $(1, 1, 1)$ , and the remainder $(-1, 0, 1)$ has length $\sqrt{2}$ . So Gram-Schmidt gives

\mathbf{q}_1 = \frac{1}{\sqrt{3}}\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}, \qquad \mathbf{q}_2 = \frac{1}{\sqrt{2}}\begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix},

and with $Q = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix}$ ,

R = \begin{bmatrix} \sqrt{3} & \sqrt{3} \\ 0 & \sqrt{2} \end{bmatrix}.

As a check, $\sqrt{3}\,\mathbf{q}_1 = (1, 1, 1)$ and $\sqrt{3}\,\mathbf{q}_1 + \sqrt{2}\,\mathbf{q}_2 = (0, 1, 2)$ , the two columns of $X$ .

Orthonormal coordinates of the data

Definition Let the orthonormal columns $\mathbf{q}_1, \dots, \mathbf{q}_n$ of $Q$ span $\operatorname{Col} X$ . The orthonormal coordinates of the data are the entries of

Q^{\top}\mathbf{y} = \begin{bmatrix} \mathbf{q}_1 \cdot \mathbf{y} \\ \vdots \\ \mathbf{q}_n \cdot \mathbf{y} \end{bmatrix}.

They are the coordinates of the fitted vector $\mathbf{p}$ , the orthogonal projection of $\mathbf{y}$ onto $\operatorname{Col} X$ , in the orthonormal basis $(\mathbf{q}_1, \dots, \mathbf{q}_n)$ , so

\begin{aligned} \mathbf{p} &= QQ^{\top}\mathbf{y} \\[0.4em] &= (\mathbf{q}_1 \cdot \mathbf{y})\,\mathbf{q}_1 + \cdots + (\mathbf{q}_n \cdot \mathbf{y})\,\mathbf{q}_n . \end{aligned}

They are coordinates along the columns of $Q$ , not along the design columns, so in general they differ from the model coefficients. Since $\mathbf{p} = X\hat{\mathbf{c}} = Q(R\hat{\mathbf{c}})$ and the columns of $Q$ are independent, the two are related by $Q^{\top}\mathbf{y} = R\hat{\mathbf{c}}$ .

In plain words Each entry of $Q^{\top}\mathbf{y}$ measures how far the data reach along one perpendicular unit direction of the model. Together the entries locate the best prediction inside the column space, but in those new directions, not in terms of intercept and slope.

Example For the three points, $\mathbf{y} = (1, 2, 2)$ and

\begin{gathered} \mathbf{q}_1 \cdot \mathbf{y} = \frac{1 + 2 + 2}{\sqrt{3}} = \frac{5}{\sqrt{3}}, \\[0.4em] \mathbf{q}_2 \cdot \mathbf{y} = \frac{-1 + 0 + 2}{\sqrt{2}} = \frac{1}{\sqrt{2}}. \end{gathered}

Then $\mathbf{p} = \tfrac53(1, 1, 1) + \tfrac12(-1, 0, 1) = (\tfrac76, \tfrac53, \tfrac{13}{6})$ , the same fitted vector that the normal equations give.

Triangular least-squares equation

Definition With $X = QR$ as above, the triangular least-squares equation is

R\mathbf{c} = Q^{\top}\mathbf{y}.

A vector $\mathbf{c}$ is a least-squares solution of $X\mathbf{c} = \mathbf{y}$ if and only if it solves this equation, and since $R$ is invertible there is exactly one. To derive it, note that the columns of $Q$ span $\operatorname{Col} X$ , so the residual is orthogonal to $\operatorname{Col} X$ exactly when $Q^{\top}(\mathbf{y} - QR\mathbf{c}) = \mathbf{0}$ , and $Q^{\top}Q = I_n$ reduces this to $Q^{\top}\mathbf{y} - R\mathbf{c} = \mathbf{0}$ . The normal equations lead to the same place, because $X^{\top}X = R^{\top}R$ and $X^{\top}\mathbf{y} = R^{\top}Q^{\top}\mathbf{y}$ , and the invertible factor $R^{\top}$ can be canceled. The equation does not claim $QR\mathbf{c} = \mathbf{y}$ , which has no solution when the data are inconsistent.

In plain words Multiplying by $Q^{\top}$ keeps only the part of the data that lies along the model’s column space. What remains is a small square triangular system, solved from the last row upward.

Example For the three points the equation reads

\begin{aligned} \sqrt{3}\,a + \sqrt{3}\,b &= \frac{5}{\sqrt{3}}, \\[0.4em] \sqrt{2}\,b &= \frac{1}{\sqrt{2}}. \end{aligned}

The last row gives $b = \tfrac12$ . Dividing the first row by $\sqrt{3}$ gives $a + b = \tfrac53$ , so $a = \tfrac53 - \tfrac12 = \tfrac76$ . This is the line $y = \tfrac76 + \tfrac12 t$ , the same answer the normal equations give, with fitted values $\tfrac76$ , $\tfrac53$ and $\tfrac{13}{6}$ .

Residual outside the model space

Definition The residual outside the model space is

\mathbf{y} - QQ^{\top}\mathbf{y} = (I - QQ^{\top})\mathbf{y},

the component of $\mathbf{y}$ orthogonal to $\operatorname{Col} X$ . It equals the least-squares residual $\mathbf{r} = \mathbf{y} - X\hat{\mathbf{c}}$ . Since $Q^{\top}(I - QQ^{\top})\mathbf{y} = \mathbf{0}$ , neither $Q^{\top}\mathbf{y}$ nor the triangular equation sees it. In particular, adding to $\mathbf{y}$ any vector $\mathbf{w}$ orthogonal to $\operatorname{Col} X$ leaves $\hat{\mathbf{c}}$ unchanged and changes the residual to $\mathbf{r} + \mathbf{w}$ .

In plain words Applying $Q^{\top}$ sets aside the part of the data that is perpendicular to every prediction the model can make. That part does not disappear. It is the residual, and solving the triangular system exactly does not make inconsistent data consistent.

Example For the three points the residual is $(-\tfrac16, \tfrac13, -\tfrac16)$ . The vector $\mathbf{w} = (1, -2, 1)$ is orthogonal to both design columns, since $1 - 2 + 1 = 0$ and $0 - 2 + 2 = 0$ . Subtracting it from the data gives $(0, 4, 1)$ . This leaves $Q^{\top}\mathbf{y}$ and the fitted line $y = \tfrac76 + \tfrac12 t$ unchanged, but the residual becomes $\mathbf{r} - \mathbf{w} = (-\tfrac76, \tfrac73, -\tfrac76)$ , and the sum of squared errors grows from $\tfrac16$ to $\tfrac{49}{6}$ .

Split of the squared error

Definition With $X = QR$ as above, for every $\mathbf{c}$ in $\mathbb{R}^n$ ,

\begin{aligned} \lVert X\mathbf{c} - \mathbf{y} \rVert^2 &= \lVert R\mathbf{c} - Q^{\top}\mathbf{y} \rVert^2 \\[0.4em] &\quad + \lVert (I - QQ^{\top})\mathbf{y} \rVert^2 . \end{aligned}

The reason is the decomposition

\begin{aligned} X\mathbf{c} - \mathbf{y} &= Q(R\mathbf{c} - Q^{\top}\mathbf{y}) \\[0.4em] &\quad - (I - QQ^{\top})\mathbf{y}, \end{aligned}

in which the first term lies in $\operatorname{Col} X$ and the second is orthogonal to it, together with $\lVert Q\mathbf{z} \rVert = \lVert \mathbf{z} \rVert$ for every $\mathbf{z}$ in $\mathbb{R}^n$ , which holds because $Q^{\top}Q = I_n$ . The first term is zero exactly when $R\mathbf{c} = Q^{\top}\mathbf{y}$ , and the second does not depend on $\mathbf{c}$ . So the smallest sum of squared errors is

\lVert (I - QQ^{\top})\mathbf{y} \rVert^2 = \lVert \mathbf{y} \rVert^2 - \lVert Q^{\top}\mathbf{y} \rVert^2 .

In plain words The error of any coefficient choice has two perpendicular parts. One measures how far $R\mathbf{c}$ is from the orthonormal coordinates of the data, and solving the triangular equation removes it. The other is the residual outside the model space, which no coefficients can remove.

Example For the three points and the line $y = 1 + t$ , so $\mathbf{c} = (1, 1)$ ,

\begin{aligned} R\mathbf{c} - Q^{\top}\mathbf{y} &= \begin{bmatrix} 2\sqrt{3} - 5/\sqrt{3} \\ \sqrt{2} - 1/\sqrt{2} \end{bmatrix} \\[0.4em] &= \begin{bmatrix} 1/\sqrt{3} \\ 1/\sqrt{2} \end{bmatrix}, \end{aligned}

whose squared length is $\tfrac13 + \tfrac12 = \tfrac56$ . The fixed term is $\lVert \mathbf{y} \rVert^2 - \lVert Q^{\top}\mathbf{y} \rVert^2 = 9 - \tfrac{25}{3} - \tfrac12 = \tfrac16$ . Their sum is $1$ , which matches the direct computation from the residuals $0$ , $0$ and $-1$ of this line.

Check yourself

A reduced QR factorization of a design matrix has $R$ with rows $\begin{bmatrix} 3 & 1 \end{bmatrix}$ and $\begin{bmatrix} 0 & 2 \end{bmatrix}$ , and the data give $Q^{\top}\mathbf{y} = (7, 4)$ . Find the least-squares coefficients.

Show the answer

The last row says $2c_2 = 4$ , so $c_2 = 2$ . The first row says $3c_1 + c_2 = 7$ , so $3c_1 = 5$ and $c_1 = \tfrac53$ . The coefficients are $(\tfrac53, 2)$ . The numbers $7$ and $4$ are coordinates along the columns of $Q$ , not coefficients of the original design columns.

Use QR to fit the line $y = a + bt$ to the points $(-1, 1)$ , $(0, 0)$ and $(1, 2)$ , and find the smallest sum of squared errors.

Show the answer

The design columns $(1, 1, 1)$ and $(-1, 0, 1)$ are already orthogonal, so Gram-Schmidt only normalizes them. That gives $\mathbf{q}_1 = (1, 1, 1)/\sqrt{3}$ , $\mathbf{q}_2 = (-1, 0, 1)/\sqrt{2}$ and

R = \begin{bmatrix} \sqrt{3} & 0 \\ 0 & \sqrt{2} \end{bmatrix}.

With $\mathbf{y} = (1, 0, 2)$ , the orthonormal coordinates are $Q^{\top}\mathbf{y} = (3/\sqrt{3}, 1/\sqrt{2}) = (\sqrt{3}, 1/\sqrt{2})$ . The triangular system $\sqrt{3}\,a = \sqrt{3}$ , $\sqrt{2}\,b = 1/\sqrt{2}$ gives $a = 1$ and $b = \tfrac12$ . The smallest sum of squared errors is $\lVert \mathbf{y} \rVert^2 - \lVert Q^{\top}\mathbf{y} \rVert^2 = 5 - 3 - \tfrac12 = \tfrac32$ . As a check, the fitted values are $\tfrac12$ , $1$ , $\tfrac32$ , the residual is $(\tfrac12, -1, \tfrac12)$ , and $\tfrac14 + 1 + \tfrac14 = \tfrac32$ .

A model has the design columns $(1, 1, 1)$ and $(2, 2, 2)$ , and the data are $\mathbf{y} = (1, 2, 6)$ . What goes wrong with the QR method, and what is the least-squares fit?

Show the answer

Gram-Schmidt gives $\mathbf{q}_1 = (1, 1, 1)/\sqrt{3}$ . The projection of the second column onto $\mathbf{q}_1$ is $2\sqrt{3}\,\mathbf{q}_1 = (2, 2, 2)$ , the whole column, so the remainder is $\mathbf{0}$ and cannot be normalized. The columns are dependent, so there is no reduced QR factorization with two orthonormal columns and an invertible $R$ , and the triangular method does not apply as stated. Least squares still has an answer. The column space is the line spanned by $(1, 1, 1)$ , and the projection of $\mathbf{y}$ onto it is $(3, 3, 3)$ , since the mean of $1$ , $2$ and $6$ is $3$ . The fitted vector $(3, 3, 3)$ is unique, with residual $(-2, -1, 3)$ and sum of squared errors $14$ . The coefficients are not unique, because every coefficient vector $(c_1, c_2)$ with $c_1 + 2c_2 = 3$ gives that fitted vector.

Least Squares with QR

Key terms

Least squares by QR factorization

Orthonormal coordinates of the data

Triangular least-squares equation

Residual outside the model space

Split of the squared error

Check yourself

How this lesson connects

Builds on

Used later in