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