next up previous
Next: Linear Models in R Up: stat6341 Previous: Practical Guidelines for Monte

Least Squares and the QR Decomposition

Suppose we would like to solve the linear system of equations,

\begin{displaymath}
Ax = b,\ \ A\in \Re^{n \times p},
\end{displaymath}

where n>p. Since there are more equations than unknowns, an exact solution does not exist. In this case we look for an approximate solution by minimizing a norm of the difference,

\begin{displaymath}
\Vert Ax-b\Vert
\end{displaymath}

When the 2-norm is used here, this is referred to as the least squares problem. This section considers the case in which A is full rank. Less than full rank problems are discussed in later sections.

An algebraic expression for the solution can be derived from the normal equations,

\begin{displaymath}
A^TAx_{ls} = A^Tb.
\end{displaymath}

These equations can be solved using the Cholesky decomposition of $A^TA$:
  1. Compute $S = A^TA$.
  2. Obtain Cholesky decomposition of S, $S=GG^T$, where G is lower triangular.
  3. Let $d=A^Tb$, solve the lower triangular system $Gy=d$ via forward substitution.
  4. Solve the upper triangular system, $G^Tx_{ls}=y$ via backward substitution.
The biggest problem associated with this approach is the increased floating point error introduced in the first step.

A more numerically stable approach is provided by the QR decomposition. Let

\begin{displaymath}
A=QR,
\end{displaymath}

where Q is an $n\times n$ orthogonal matrix and R is an $n\times p$ upper triangular matrix. If A is full rank, then all of the diagonal elements of R are non-zero. This decomposition always exists but it is not necessarily unique. Multiplication of the ith column of Q and the ith row of R by -1 does not change the product. Therefore, the QRD is unique if a requirement is added that diagonal elements of R must be positive.

Properties of the QRD. Let

\begin{displaymath}
Q = [Q_1\ Q_2],\ Q_1\in \Re^{n\times p},\ Q_2\in \Re^{n\times (n-p)}.
\end{displaymath}

  1. The columns of $Q_1$ form an orthonormal basis for $range(A)$ and the columns of $Q_2$ form an orthonormal basis for $range(A)^\perp$. More generally, let

    \begin{eqnarray*}
A &=& [a_1\ \cdots \ a_k],\ 1\le k\le p\\
Q &=& [q_1\ \cdots \ q_k],\ 1\le k\le p
\end{eqnarray*}

    denote submatrices of A,Q, respectively. Then

    \begin{displaymath}
span\{a_1,\cdots,a_k\} = span\{q_1,\cdots,q_k\},\ 1\le k\le p.
\end{displaymath}

  2. $A = Q_1R_1$, where $R_1$ contains the first p rows of R. Note that the remaining rows of R are 0 and that $R_1$ is upper triangular. This is referred to as the reduced QRD.

  3. \begin{displaymath}
A^TA = (Q_1R_1)^T(Q_1R_1) = R_1^TQ_1^TQ_1R_1 = R_1^TR_1.
\end{displaymath}

    Hence, $R_1^T$ is the Cholesky factor of the symmetric matrix $A^TA$.
  4. $Q_1Q_1^T$ is the projection matrix onto $range(A)$.
  5. $I=Q_1Q_1^T$ is the projection matrix onto $range(A)^\perp$.

To show how these properties generate solutions to least squares problems, let $A=QR$ be the full QRD of A. Use the column partition of Q defined above to obtain

\begin{displaymath}
Q^Tb = \left[
\begin{array}{l}
Q_1^Tb\\
Q_2^Tb
\end{array}\right] = \left[
\begin{array}{l}
c\\
d
\end{array}\right]
\end{displaymath}

$Q^T$ is an orthogonal matrix and so

\begin{eqnarray*}
\Vert Ax-b\Vert _2^2 &=& \Vert Q^T(Ax-b)\Vert _2^2\\
&=& \Ve...
...\Vert _2^2\\
&=& \Vert R_1x - c\Vert _2^2 + \Vert d\Vert _2^2.
\end{eqnarray*}

Since d does not involve x, this is minimized by solving the square, upper triangular system of equations

\begin{displaymath}
R_1x = Q_1^T.
\end{displaymath}

Note that just the reduced QRD is needed to obtain the LS solution.

Fitted values are given by

\begin{displaymath}
\hat{b} = Ax_{ls} = Q_1R_1x_{ls} = Q_1Q_1^Tb,
\end{displaymath}

so fitted values are the projection of b onto $range(A)$. Residuals are given by

\begin{displaymath}
e = b - \hat{b} = b - Q_1Q_1^Tb = (I-Q_1Q_1^T)b,
\end{displaymath}

which is the projection of b onto $range(A)^\perp$.


next up previous
Next: Linear Models in R Up: stat6341 Previous: Practical Guidelines for Monte
Larry Ammann
2017-11-01