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)? 

As far as I know and from what I have read, there is no direct formula (or algorithm) for the (left) psuedo-inverse of $A$, other than the "famous formula"
$$(A^\top A)^{-1}A^\top$$
which is computationally expensive. I will describe here an algorithm that does indeed calculate it efficiently. As I said, I believe this may be previously
un-published, so cite me appropriately, or if you know of a reference please state it. Here is my algorithm.

Set $B_0 = A$ and for each iteration step, take a column of $B_i$ and orthogonalize against the columns of $A$. Here is the algorithm ($A$ has $n$ independent columns):

 1. Initialize $B_0=A$. 
 2. For $j \in 1\ldots n$ do
\begin{align}
  \mathbf{t} &= B_i^\top A_{\star j} \\
  R\mathbf{t} &= e_j & \text{Find such $R$, an elementary row operation matrix}\\
  B_{i+1}^\top &= R B_i^\top \\
\end{align}
Notice that each step does one vector/matrix multiplication and one elementary matrix row operation. This will require
much less computation than the direct evaluation of $(A^\top A)^{-1}A^\top$. Notice also that there is a nice opportunity here for parallel computation. (One more thing to notice--this may calculate a regular inverse, starting with $B_0=A$, or with $B_0=I$.)
The step that calculates $R$ may partially be orthogonal
 rather than elementary (orthonormal/unitary -thus better
in the sense of error propagation)
if desired, but should not make use of the previously orthonormalized rows of $B$, since they must remain orthogonalized in the process. If this is done, the process becomes a "cousin" to the QR algorithm, whereas before it would be considered a "cousin" to the LU factorization. "Cousin" because the matrix operations of those are self referential, and this algorithm references the matrix $A$.
The following is a hopefully enlightening
description.

See first edits for a more wordy description. But I think the following is more concise and to the point.

Consider the conformable matrix $B^\top$ such that $B$ has the same dimension as $A$ and
\begin{align}
 B^\top A &= I \tag{1}\\
 B &= A B_A \tag{2}\\
\end{align}
Here $B_A$ is arbitrary and exists only to ensure that the matix $B$ shares the same column space as $A$. (1) states that $B^\top$ is a (not necessarily unique) left inverse for $A$. (2) states that $B$ shares the same column space as $A$.
Claim: $B^\top$ is the Moore-Penrose pseudoinverse for $A$.
Proof:
The Moore-Penrose pseudoinverse for $A$ is $$A^\dagger = \left(A^\top A\right)^{-1}A^\top$$
Given (1) and substituting (2) we have
\begin{align}
 \left(AB_A\right)^\top A &= I \\
 B_A^\top A^\top A &= I \\
 B_A^\top &= \left( A^\top A\right)^{-1} \\
 B_A &= \left( A^\top A\right)^{-1} \\
\end{align}
Now solving for $B$ from (2):
\begin{align}
  B &= A\left(B_A\right) \\
  B &= A\left( A^\top A\right)^{-1} \\
 \Rightarrow B^\top &=\left( A^\top A\right)^{-1}A^\top \\
\end{align}
Q.E.D.
