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

-

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$

-

First step is to buy Matlab. Second step is to type pinv.

Edit: Here is an implementation.

-
Funny ^^ but hardly an answer to his question. –  Olivier Bégassat Oct 25 '11 at 17:07
or get octave for free=) –  Kabumbus Oct 25 '11 at 17:38

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.

-
But the question is about the Moore-Penrose pseudoinverse, not the left inverse. If left inverses exist, the Moore-Penrose pseudoinverse is one of them, isn't it? And it can be computed using the algorithm in J.M.'s answer. So I'm not sure what you mean by "there is no direct formula (or algorithm)". –  Rahul Feb 28 '13 at 20:17
@Rahul The algorithm in J.M.'s answer is iterative. What I am saying is that this algorithm computes the particular left inverse which is indeed the Moore-Penrose pseudoinverse, and does so in exact arithmetic, not by convergence. –  adam W Feb 28 '13 at 20:29
I liked your idea much, so I tried to implement it to see, whether I've got it correctly. Here is a webpage where I just put my procedere together:go.helms-net.de/math/divers/… . Has that taken your proposal correctly? –  Gottfried Helms Feb 28 '13 at 21:45
@Gottfried Well, I do get the same numbers up to accuracy. For instance, $-0.012390122873089741$ is what I get for the lower right element. The way you describe "putting A into upper triangular shape" seems like you are doing regular LU or QR though...so not too certain. –  adam W Feb 28 '13 at 23:02
Is this faster than SVD? It's $\Omega(mn^2)$, right? –  user7530 Mar 1 '13 at 2:05