What is step by step logic of pinv (pseudoinverse)? So we have a matrix $A$ size of $M \times N$ with elements $a_{i,j}$. What is a step by step algorithm that returns the Moore-Penrose inverse $A^+$ for a given $A$ (on level of manipulations/operations with $a_{i,j}$ elements, not vectors)? 

How to we use the elements of $\mathbf{A}$ to construct the pseudoinverse $\mathbf{A}^{\dagger}$?
Anatomy of the SVD
The target matrix has $m$ rows, $n$ columns, and has rank $\rho$
$$
\mathbf{A} \in \mathbb{C}^{m\times n}_{\rho}
$$
Define the singular value decomposition:
$$
\begin{align}
  \mathbf{A} &=
  \mathbf{U} \, \Sigma \, \mathbf{V}^{*} \\
%
 &=
% U 
  \left[ \begin{array}{cc}
     \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}}
  \end{array} \right]  
% Sigma
  \left[ \begin{array}{cccc|cc}
     \sigma_{1} & 0 & \dots &  &   & \dots &  0 \\
     0 & \sigma_{2}  \\
     \vdots && \ddots \\
       & & & \sigma_{\rho} \\\hline
       & & & & 0 & \\
     \vdots &&&&&\ddots \\
     0 & & &   &   &  & 0 \\
  \end{array} \right]
% V 
  \left[ \begin{array}{c}
     \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \\ 
     \color{red}{\mathbf{V}_{\mathcal{N}}}^{*}
  \end{array} \right]  \\
%
  & =
% U
   \left[ \begin{array}{cccccccc}
    \color{blue}{u_{1}} & \dots & \color{blue}{u_{\rho}} & \color{red}{u_{\rho+1}} & \dots & \color{red}{u_{n}}
  \end{array} \right]
% Sigma
  \left[ \begin{array}{cc}
     \mathbf{S}_{\rho\times \rho} & \mathbf{0} \\
     \mathbf{0} & \mathbf{0} 
  \end{array} \right]
% V
   \left[ \begin{array}{c}
    \color{blue}{v_{1}^{*}} \\ 
    \vdots \\
    \color{blue}{v_{\rho}^{*}} \\
    \color{red}{v_{\rho+1}^{*}} \\
    \vdots \\ 
    \color{red}{v_{n}^{*}}
  \end{array} \right]
%
\end{align}
$$
Coloring separates $\color{blue}{range}$ space entities from $\color{red}{null}$ space entities.
Using the thin, or economical, or nullspace free variants, we have
$$
  \mathbf{A} = \color{blue}{\mathbf{U}_{\mathcal{R}}} \, \mathbf{S} \, \color{blue}{\mathbf{V}^{*}_{\mathcal{R}}}
\qquad \Rightarrow \qquad
\mathbf{A}^{\dagger} = 
\color{blue}{
\mathbf{V}_{\mathcal{R}}
} \, \mathbf{S}^{-1} \, \color{blue}%{\mathbf{U}^{*}_{\mathcal{R}}}.
$$
Construction of the SVD


*

*Form the product matrix $\mathbf{A}^{*}\mathbf{A}$.

*Solve for the eigenvalues $\lambda \left( \mathbf{A}^{*}\mathbf{A}\right)$. This is a messy affair because it demands finding the roots of the characteristic equation which is the numerical challenge of solving for the zeros of an arbitrary polynomial.

*Order the nonzero values to assemble the matrix of singular values $\mathbf{S}=\text{diagonal}\left( \sigma \right)$ where the singular value spectrum is
$$
  \sigma = \sqrt{\lambda \left( \mathbf{A}^{*}\mathbf{A} \right)}
$$
with
$$
  \sigma_{1} \ge \sigma_{2} \ge \dots \sigma_{\rho} > 0.
$$
A problem immediately arises over the value of $\rho$. Are small numbers numerical artifacts? The issue is vital because the contribution to the pseudoinverse is weighted by the reciprocal of the singular value. Very small number make very large contributions. 

*Compute the eigenvectors $\hat{u}_{k}$, $k=1,\rho$ for each eigenvalue:
$$
 \mathbf{A}^{*}\mathbf{A}\, \hat{u}_{k} = \lambda_{k}\,\hat{u}_{k} 
$$
This is the origin of the sign ambiguity in the SVD.

*Normalize these vectors and make them the column vectors of the codomain matrix $\color{blue}{\mathbf{U}^{*}_{\mathcal{R}}}$:
$$
  \color{blue}{\mathbf{U}^{*}_{\mathcal{R}}} = \left[ \begin{array}{cccc}
  \color{blue}{\frac{ \hat{u}_{1} } {\lVert \hat{u}_{1} \rVert}} &
  \color{blue}{\frac{ \hat{u}_{2} } {\lVert \hat{u}_{2} \rVert}} & \dots &
  \color{blue}{\frac{ \hat{u}_{\rho} } {\lVert \hat{u}_{\rho} \rVert}}
\end{array} \right]
$$

*Optional: Want $\color{red}{\mathbf{U}_{\mathcal{N}}}$? Apply Gram-Schmidt orthogonalization to the column vectors. The orientation is a free choice and a reason while different programs create different nullspace spans. If you choose this option, you need to pack $\mathbf{S}$ into a sabot matrix $\Sigma$ and compute the companion nullspace for $\mathbf{V}$.

*With $\mathbf{S}$ and $\color{blue}{\mathbf{U}^{*}_{\mathcal{R}}}$ in hand, construct the column vectors of $\color{blue}{\mathbf{V}^{*}_{\mathcal{R}}}$ using
$$
  \left[ \color{blue}{\mathbf{V}_{\mathcal{R}}} \right]_{k} = \sigma_{k}^{-1} \left[ 
\color{blue}{\mathbf{U}^{*}_{\mathcal{R}}} \mathbf{A} 
\right]_{k}
$$
Notice the reciprocal of the singular value.

*Optional: Want $\color{red}{\mathbf{V}_{\mathcal{N}}}$? Apply Gram-Schmidt orthogonalization to the column vectors of $\color{red}{\mathbf{V}_{\mathcal{N}}}$.

*Invert the diagonal matrix $\mathbf{S}$:
$$
 \mathbf{S} = 
\left[ \begin{array}{ccc}
  \sigma_{1} \\ & \sigma_{2} \\ && \ddots \\ &&& \sigma_{\rho}
\end{array} \right]
%
\qquad \Rightarrow \qquad
%
 \mathbf{S}^{-1} = 
\left[ \begin{array}{ccc}
  \frac{1}{\sigma_{1}} \\ & \frac{1}{\sigma_{2}} \\ && \ddots \\ &&& \frac{1}{\sigma_{\rho}}
\end{array} \right]
$$
Again, the reciprocal of the singular values.

*Form the Hermitian conjugate by taking the transpose of the conjugate or the conjugate of the transpose:
$$
  \color{blue}{\mathbf{V}^{*}_{\mathcal{R}}} = 
\color{blue}{\overline{\mathbf{V}}_{\mathcal{R}}^{\mathrm{T}}} =
\color{blue}{\overline{\mathbf{V}^{\mathrm{T}}_{\mathcal{R}}}}
$$
Construction of the pseudoinverse


*Assemble
$$
  \mathbf{A}^{\dagger} = 
\color{blue}{\mathbf{V}_{\mathcal{R}}}
\, \mathbf{S}^{-1} \,
\color{blue}{\mathbf{U}^{*}_{\mathcal{R}}} =
  \left[ \begin{array}{cc}
     \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}}
  \end{array} \right]  
% Sigma
  \left[ \begin{array}{cccccc}
    \mathbf{S}^{-1} &  \mathbf{0} \\
     \mathbf{0} &  \mathbf{0}
  \end{array} \right]
% V 
  \left[ \begin{array}{c}
     \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \\ 
     \color{red}{\mathbf{V}_{\mathcal{N}}}^{*}
  \end{array} \right]  
$$

