Let's say we have the following system of equations: \begin{equation} A{\bf x}={\bf b} \qquad (1) \end{equation} where $A$ is $N \times 4$, $\mathbf{x}$ is $4 \times 1$ (unknowns) and ${\bf b}$ is $N \times 1$.
Using linear-least square method the solution of the overdetermined system is:
\begin{equation} (A^TA) {\bf x} = A^T {\bf b} \qquad (2) \end{equation}
\begin{equation} {\bf x} = (A^TA)^{-1} A^T {\bf b} \qquad (3) \end{equation}
Why do we multiply both sides of Eq.$(1)$ by $A^T$?
You minimize the the (squared) expression $V=\left\| Ax-b \right\|^2$. This is equal to
$V=(Ax-b)^T(Ax-b)$
Remove $T$ from the brackets. The order of $Ax$ is changing.
$V=(x^TA^T-b^T)(Ax-b)$
Multiplying out the brackets.
$V=x^TA^TAx-x^TA^Tb-b^TAx+b^Tb$
$x^TA^Tb$ and $b^TAx$ are both scalars and they are equal. One can be replaced by another.
$V=x^TA^TAx-2x^TA^Tb+b^Tb$
To optimize this expression w.r.t $x$ it has to be calculated the derivative w.r.t $x $ and set it equal to zero. The derivative of $b^Tb$ is zero, because it has no $x$.
$\frac{\partial V}{\partial x}=2A^TAx-2A^Tb=0 \quad |:2$
$A^TAx-A^Tb=0 \quad \quad$
$A^TAx=A^Tb \quad\quad |:A^TA $
On the LHS $A^TA$ is left from $x$. If we divide both sides by $(A^TA)^{}$ then $(A^TA)^{-1}$ has to be left on the RHS as well.
$\boxed{x=(A^TA)^{-1}A^Tb}$
Linear least square problem solves the probelm $$\min \frac{1}{2}\left\| Ax-b \right\|^2$$
To solve the problem, we differentiate with respect to $x$ and equare it to zero and we end up with
$$A^T(Ax-b)=0$$
Let's assume the linear system $$ \mathbf{A}x = b $$ has no solution and that the data vector $b\notin\mathcal{N}\left( \mathcal{A} \right).$
There is no exact solution because the data vector $b$ is not in the column space of $\mathbf{A}$.
One path to solution involves formulating a problem with the same solution that does have a solution. This is the normal equations. $$ \mathbf{A}^{*}\mathbf{A}x = \mathbf{A}^{*}b. $$ The vector on the right, $\mathbf{A}^{*}b$ is clearly in the column space of $\mathcal{A}^{*}$. In fact we are given the prescription for combining the column vectors $$ b_{1} \left[ \mathbf{A}^{*} \right]_{1} + b_{2} \left[ \mathbf{A}^{*} \right]_{2} + \dots b_{m} \left[ \mathbf{A}^{*} \right]_{m}. $$
Now we can use tools like Gaussian elimination and $\mathbf{L}\mathbf{U}$ decomposition to solve the problem.
In numerical problems, the normal equations can be a poor choice and benefit from methods like the $\mathbf{Q}\mathbf{R}$ and singular value decompositions.
