Why does this hacky derivation for least-squares regression work?

Question

Consider the regression problem where we have $m$ measurements of the dependent variable and a model with $n$ degrees of freedom, where $m>n$. We can write the dependent variable measurements in a vector $\mathbf{y}$ (size $m\times 1$) and the model parameters in a vector $\mathbf{x}$ (size $n\times 1$), and arrange the independent measurements appropriately in a matrix $\mathbf{A}$ (size $m\times n$). The problem is now to choose $\mathbf{x}$ such that $\mathbf{y}$ is represented as closely as possible by $\mathbf{A}\mathbf{x}$. The most common way of quantifying "as closely as possible" is in the least-squares sense. We write: $$ E = \left(\mathbf{y} - \mathbf{A}\mathbf{x}\right)^T\left(\mathbf{y} - \mathbf{A}\mathbf{x}\right) \tag{1} $$

By taking the derivative of $E$ wrt $\mathbf{x}$ and setting it to zero, we end up with the least-squares solution:

$$ \mathbf{x}=\left(\mathbf{A}^T \mathbf{A}\right)^{-1} \mathbf{A} \mathbf{y} \tag{2} $$

---

I once had a professor show me a hacky way to have to neither remember this formula, nor do the tedious derivation in an exam situation. He was very clear that it was a hack and that I should never use it as a serious derivation. It goes as follows - start out by writing:

$$ \mathbf{A}\mathbf{x} = \mathbf{y} \tag{3} $$

Observe that we cannot solve this equation by taking the inverse of $\mathbf{A}$ because this is not a square matrix. No problem! - multiply each side by $A^T$ (size $n\times m$) to get: $$ \left(\mathbf{A}^T\mathbf{A}\right) \mathbf{x} = \mathbf{A}^T \mathbf{y} \tag{4} $$ Now $A^T A$ is a square matrix (size $n\times n$), which means it's (potentially) invertible. Multiply each side by $\left(\mathbf{A}^T \mathbf{A}\right)^{-1}$ to get: $$ \mathbf{x} = \left(\mathbf{A}^T \mathbf{A}\right)^{-1} \mathbf{A} \mathbf{y} \tag{5} $$

We get the right result, in the least-squares sense! Of course the hack is that, in general, Eq. (3) is not true to begin with; it is an inconsistent equation with no solution.

---

My question is: Is it just pure coincidence that this hack works in this particular case? Or is there perhaps a deeper reason? Maybe an insight as to why it leads to the least-squares solution as opposed to other one? Thank you!

Answer 1

If $\|y - Ax\|$ is minimized, then $Ax$ is the closest point to $y$ in the image space $\operatorname{im}(A)$ of $A$; therefore the residual $r = y - Ax$ is orthogonal to $\operatorname{im}(A)$. Now $\operatorname{im}(A)^\perp = \ker(A^T)$ in general, so $r \in \ker(A^T)$ and $$ A^T y = A^T (Ax + r) = A^T Ax + A^T r = A^TAx. $$

Answer 2

Why this works is explained in Overdetermined System Ax=b.

How does this connect to least squares? Start with the sequence of measurements $\left\{ x_{k}, y_{k} \right\}_{k=1}^{m}$, and the trial function $$ y(x) = \sum_{k=1}^{n} c_{k} f_{k}(x) $$ where the functions $f$ are linearly independent. Use least squares to find the best amplitudes $c$.

The linear system is $$ \begin{align} \mathbf{A} c &= Y \\ % \left[ \begin{array}{cccc} f_{1}\left( x_{1} \right) & f_{2}\left( x_{1} \right) & \dots & f_{n}\left( x_{1} \right) \\ % f_{1}\left( x_{2} \right) & f_{2}\left( x_{2} \right) & \dots & f_{n}\left( x_{2} \right) \\ % \vdots & \vdots & & \vdots \\ % f_{1}\left( x_{m} \right) & f_{2}\left( x_{m} \right) & \dots & f_{n}\left( x_{m} \right) \\ % \end{array} % \right] \left[ \begin{array}{c} c_{1} \\ c_{2} \\ \vdots \\ c_{n} \end{array} \right] % &= % \right] \left[ \begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{array} \right] % \end{align} $$

The normal equations are $$ \begin{align} \mathbf{A}^{*}\mathbf{A} c &= \mathbf{A}^{*}Y \\ % \left[ \begin{array}{ccc} \sum_{k=1}^{m} f_{1}\left( x_{k} \right)f_{1}\left( x_{k} \right) & \sum_{k=1}^{m} f_{1}\left( x_{k} \right)f_{2}\left( x_{k} \right) & \dots & \sum_{k=1}^{m} f_{1}\left( x_{k} \right)f_{n}\left( x_{k} \right) \\ % \sum_{k=1}^{m} f_{2}\left( x_{k} \right)f_{1}\left( x_{k} \right) & \sum_{k=1}^{m} f_{2}\left( x_{k} \right)f_{2}\left( x_{k} \right) & \dots & \sum_{k=1}^{m} f_{2}\left( x_{k} \right)f_{n}\left( x_{k} \right) \\ % \vdots & \vdots && \vdots \\ % \sum_{k=1}^{m} f_{n}\left( x_{k} \right)f_{1}\left( x_{k} \right) & \sum_{k=1}^{m} f_{n}\left( x_{k} \right)f_{2}\left( x_{k} \right) & \dots & \sum_{k=1}^{m} f_{n}\left( x_{k} \right)f_{n}\left( x_{k} \right) \\ % \end{array} \right] \left[ \begin{array}{c} c_{1} \\ c_{2} \\ \vdots \\ c_{n} \end{array} \right] % &= \left[ \begin{array}{c} \sum_{k=1}^{m} f_{1}\left( x_{k} \right) y\left( x_{k} \right) \\ \sum_{k=1}^{m} f_{2}\left( x_{k} \right) y\left( x_{k} \right) \\ \vdots \\ \sum_{k=1}^{m} f_{n}\left( x_{k} \right) y\left( x_{k} \right) \end{array} \right] % \end{align} $$

We recover the normal equations when we attack the minimization problem with the calculus. The goal is the minimize the vector function: $$ r^{2}(c) = \sum_{k=1}^{m} \left( y_{k} - c_{1}f_{1}\left( x_{k} \right) - c_{2}f_{2}\left( x_{k} \right) - \dots - c_{n}f_{n}\left( x_{k} \right) \right)^{2} $$ The method is to demand that the $n$ first derivates simultaneously be 0. $$ \frac{\partial} {\partial c_{j}} r^{2}(c) = -2 \sum_{k=1}^{m} \left( y_{k} - c_{1}f_{1}\left( x_{k} \right) - c_{2}f_{2}\left( x_{k} \right) - \dots - c_{n}f_{n}\left( x_{k} \right) \right) f_{j}(x_{k}) = 0 $$ Reorder this equation as $$ c_{1} \sum_{k=1}^{m} f_{1}\left( x_{k} \right) f_{j}(x_{k}) + c_{2} \sum_{k=1}^{m} f_{2}\left( x_{k} \right) f_{j}(x_{k}) + \dots + c_{n} \sum_{k=1}^{m} f_{n}\left( x_{k} \right) f_{j}(x_{k}) = \sum_{k=1}^{m} f_{j}(x_{k}) y\left( x_{k} \right) $$ to see that it is row $j$ in the matrix equation for the normal equations.