Difference between orthogonal projection and least squares solution When you find the least squares solution you solve $$A^TA = A\vec b$$ but to find the orthogonal projection into the "subspace" A, you multiply this result (the least squares solution) with the original matrix. Why is this? 
If you use the analogy with the light shining orthogonally on to the subspace and the orthogonal projection is the shadow in the subspace, isn't this shadow also the least squares solution?
 A: You want to minimize $||Ax-b||^2$. If $A$ is a full rank matrix, the unique solution is given by
$$x=(A^TA)^{-1}A^Tb$$
So $Ax=A(A^TA)^{-1}A^Tb=Hb$, with $H=A(A^TA)^{-1}A^T$ the projection matrix on the subspace formed by the columns of $A$.
A: Given the linear system
$$
\mathbf{A} x = b
$$
we look for solutions of the form 
$$
\mathbf{A} x - b = 0.
$$
If the data vector $b$ is in the column space of $\mathbf{A}$, the the above equation has an exact solution,
$$
 x = \mathbf{A}^{-1}b.
$$
If $b$ has a nullspace component, we can't satisfy the difference equation, so we relax the requirement and ask the $\mathbf{A}x - b$ be as small as possible. This leads to the definition of the least squares problem: 
$$
x_{LS} = \left\{ x\in\mathbb{C}^{n} \colon \lVert \mathbf{A} x_{LS} - b \rVert_{2}^{2} \text{ is minimized} \right\}
$$
The general least squares solution is
$$
 x_{LS} = \mathbf{A}^{\dagger} b + \left( \mathbf{I}_{n} - \mathbf{A}^{\dagger}\mathbf{A} \right)y, \quad y \in\mathbb{C}^{n}.
$$
The are multiple avenues for solution. For example, the normal equations which you allude to:
$$
\mathbf{A}^{*} \mathbf{A} x = \mathbf{A}^{*} b 
$$
which offers the solution
$$
x_{LS} = \left( \mathbf{A}^{*}\mathbf{A} \right)^{-1} \mathbf{A}^{*} b
$$
(Construction of the normal equations: Overdetermined System Ax=b)
The fundamental projectors are defined using the SVD in Least squares solutions and the orthogonal projector onto the column space. You are asking about the projector onto the range space of $\mathbf{A}$:
$$
\mathbf{P}_\color{blue}{\mathcal{R}\left( \mathbf{A} \right)} = \mathbf{A}\mathbf{A}^{\dagger}
$$
The operator projects $m$ vectors on the range space (column space, image).
If the classic inverse exists, then it is also the pseudoinverse Pseudo inverse of a singular value decomposition SVD is equal to its “real” inverse for a square matrix?
For your normal equations example, the inverse is also the pseudoinverse. Existence of a normal equations solution implies we have at least as many columns as rows $m\ge n$, and the matrix has full column rank $m=\rho$. The singular value decomposition is
$$
 \mathbf{A} = \mathbf{U}\, \Sigma \, \mathbf{V}^{*} 
=
\underbrace{\left[ \begin{array}{cc}
  \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}}
\end{array} \right]}_{m\times m}
%
\underbrace{\left[ \begin{array}{c}
  \mathbf{S}_{\rho \times \rho} \\\mathbf{0}
\end{array} \right]}_{m\times n}
%
\underbrace{\left[ \begin{array}{c}
  \color{blue}{\mathbf{V}_{\mathcal{R}}}
\end{array} \right]^{*}}_{n\times n}
.
$$
This implies the pseudoinverse is
$$
 \mathbf{A}^{\dagger} = \mathbf{V}\, \Sigma^{\dagger} \, \mathbf{U}^{*} 
=
%
\underbrace{\left[ \begin{array}{c}
  \color{blue}{\mathbf{V}_{\mathcal{R}}}
\end{array} \right]}_{n\times n}
%
\underbrace{\left[ \begin{array}{cc}
  \mathbf{S}_{\rho \times \rho}^{-1} & \mathbf{0}
\end{array} \right]}_{n\times m}
%
\underbrace{\left[ \begin{array}{cc}
  \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}}
\end{array} \right]^{*}}_{m\times m}
$$
The normal equations solution is
$$
\begin{align}
  \left( \mathbf{A}^{*}\mathbf{A} \right)^{-1} \mathbf{A}^{*} 
%
  &=
% 
  \left(  
  \left( \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{*} \right) 
  \left( \mathbf{U}\, \Sigma           \, \mathbf{V}^{*} \right) 
  \right)^{-1} 
    \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{*} \\
%
  &=
%
  \left(  
   \mathbf{V}\, \mathbf{S}^{2} \mathbf{V}^{*} 
  \right)^{-1}
    \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{*} \\
%
  &=
%
  \left(  
   \mathbf{V}\, \mathbf{S}^{-2} \mathbf{V}^{*} 
  \right)
    \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{*} \\
%
  &=
%  \left(  
   \mathbf{V}\, \Sigma^{\dagger} \mathbf{U}^{*} \\
%
 &=
%
 \mathbf{A}^{\dagger}
%
\end{align}
$$
A: A least squares solution is not the shadow you refer to in the shining light analogy.  This shadow is the orthogonal projection of $b$ onto the column space of $A$, and it is unique.  Call this projection $p$.  A least squares solution of $Ax = b$ is a vector $x$ such that $Ax = p$. The vector $x$ need not be unique.
Consider the matrix
$$A = \begin{bmatrix}
4 &8 \\
6 &12 
\end{bmatrix}$$
and the vector
$$b = \begin{bmatrix}
5 \\
1
\end{bmatrix}
$$
which is not in $C(A) = \textrm{span} \left( \begin{bmatrix} 4 \\ 6 \end{bmatrix} \right)$.  The orthogonal projection of $b$ onto $C(A)$ is given by
$$ p = \frac{\begin{bmatrix} 4 \\ 6 \end{bmatrix} \cdot \begin{bmatrix} 5 \\ 1 \end{bmatrix}}{\begin{bmatrix} 4 \\ 6 \end{bmatrix} \cdot \begin{bmatrix} 4 \\ 6 \end{bmatrix}} \begin{bmatrix} 4 \\ 6 \end{bmatrix} = \begin{bmatrix} 2 \\ 3 \end{bmatrix}$$
A least squares solution of $Ax = b$ is a vector $x$ such that
$$Ax = p$$
This system has infinitely many solutions.  The solution set is
$$x = \left\{ t \begin{bmatrix} -2 \\ 1 \end{bmatrix} + \begin{bmatrix} \frac{1}{2} \\ 0 \end{bmatrix}, t \in \mathbb{R} \right\} $$
Therefore, both $x_1 = \begin{bmatrix} \frac{1}{2} \\ 0 \end{bmatrix}$ and $x_2 = \begin{bmatrix} -\frac{7}{2} \\ 2 \end{bmatrix}$, for instance, are least squares solutions, because both $Ax_1 = p$ and $Ax_2 = p$.  But neither of these solutions is the "shadow" you refer to in the shining light analogy.  Rather, $p$ is the shadow, and $x_1$ and $x_2$ are simply vectors you could multiply $A$ by to get $p$.

