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!
 A: 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. $$
A: 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.
