What is step by step logic of pinv (pseudoinverse)?

Question

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

Answer 1

1. Perform a singular value decomposition $\mathbf A=\mathbf U\mathbf \Sigma\mathbf V^\top$

2. Check for "tiny" singular values in $\mathbf \Sigma$. A common criterion is that anything less than $\|\mathbf A\|\epsilon$ is "tiny". Set those singular values to zero.

3. Form $\mathbf \Sigma^+=\mathrm{diag}(1/\sigma_1,1/\sigma_2,\dots, 0)$. That is, reciprocate the nonzero singular values, and leave the zero ones untouched.

4. $\mathbf A^+=\mathbf V\mathbf \Sigma^+\mathbf U^\top$

Answer 2

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.