The most straightforward way I know is to pass through the normal equations:

$$A^T A x = A^T b$$

and substitute in the$QR$ decomposition of$A$ (with the convention$Q \in \mathbb{R}^{m \times n},R \in \mathbb{R}^{n \times n}$). Thus you get

$$R^T Q^T Q R x = R^T Q^T b.$$

But$Q^T Q=I_n$. (Note that in this convention$Q$ isn't an orthogonal matrix, so$Q Q^T \neq I_m$, but this doesn't matter here.) Thus:

$$R^T R x = R^T Q^T b.$$

If$A$ has linearly independent columns (as is usually the case with overdetermined systems), then$R^T$ is injective, so by multiplying both sides by the left inverse of$R^T$ you get

$$Rx=Q^T b.$$

This system is now easy to solve numerically.

For numerical purposes it's important that the removal of$Q^T Q$ and$R^T$ from the problem is done analytically, and in particular$A^T A$ is never constructed numerically.

• 1
  $\begingroup$Is there any reason to make this so convoluted? From $Ax=b$ you have $QRx=b$, multiply by $Q^T$ on the left.$\endgroup$

  Martin Argerami

  – Martin Argerami

  2019-04-12 15:45:35 +00:00

  CommentedApr 12, 2019 at 15:45
• $\begingroup$@MartinArgerami Because actually the least squares solution usually does not satisfy $Ax=b$. This simple perspective only shows you that this approach gives you a solution when a solution exists. Now you could argue directly that multiplying both sides by $Q^T$ furnishes an equation whose solution is the least squares solution. (Such an argument would resemble the usual geometric argument for deriving the normal equations.) This would make a good alternative answer to mine.$\endgroup$

  Ian

  – Ian

  2019-04-12 15:46:31 +00:00

  CommentedApr 12, 2019 at 15:46

Note that$Rx$ has the form$$Rx = \begin{bmatrix} y_1 \\ y_2 \\ 0\end{bmatrix} $$, so if$$ Q^{-1}b = \begin{bmatrix} z_1 \\ z_2 \\ z_3\end{bmatrix}$$then$|| Rx - Q^{-1}b||$ will be minimal for$y_1 = z_1$,$y_2=z_2$. This set of equation is no longer overdetermined.

Using matrix notation, if tou write$R = \begin{bmatrix} R_1 \\ 0\end{bmatrix}$ and intoduce$P=\begin{bmatrix}1 & 0 & 0 \\ 0 & 1& 0\end{bmatrix}$, then you have$$ R_1x = PQ^{-1}b$$$$ x = (R_1)^{-1}PQ^{-1}b$$

• 1
  $\begingroup$The key trick in this answer is that by the orthogonality, $\| Ax - b \| = \| Rx - Q^T b \|$.$\endgroup$

  Ian

  – Ian

  2019-04-12 20:34:30 +00:00

  CommentedApr 12, 2019 at 20:34
• $\begingroup$@Ian, That's something that OP has alredy obtained on his own (since $Q$ is orthogonal, $Q^{-1}=Q^T$).$\endgroup$

  Adam Latosiński

  – Adam Latosiński

  2019-04-13 11:14:27 +00:00

  CommentedApr 13, 2019 at 11:14

• linear-algebra
• numerical-methods
• numerical-linear-algebra