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

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!

• The solution to the normal equations $$A^T A x = A^T b$$ is, in general, not a solution to the original linear system $A x = b$. Note that the normal equations always have a solution. To show that, use the SVD. Commented Oct 11, 2016 at 18:49
• It's not coincidence. You have the orthogonal projection $P$ onto the range of $A$. Write $y = Py + (y - Py)$. What you look for is an $x$ with $Ax = Py$. Now $A^T(y-Py) = 0$ since $y-Py \in (\operatorname{range} A)^{\perp}$, so $A^Ty = A^TPy = A^TAx$. Commented Oct 11, 2016 at 18:49
• As $x$ varies, $Ax$ ranges over all vectors in the column space of $A$. When $x$ is chosen so that $Ax$ is as close as possible to $b$, you can see visually that the residual $b - Ax$ is orthogonal to the column space of $A$. Equivalently, $b - Ax$ is orthogonal to each column of $A$. In other words, $A^T(b - Ax) = 0$. Commented Mar 8, 2017 at 23:41
• adding to @littleO 's answer: first draw the simplest case, just one column, $b$ and $A$ in the plane. Then for the multi-column case, if $b - Ax$ isn't orthogonal to say the first column $A_1$ of $A$, you could change $x_1$ to get closer to $b$ . Commented Nov 25, 2017 at 16:57

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.$$

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.