Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
add comment

3 Answers 3

up vote 11 down vote accepted
+50
  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$

share|improve this answer
add comment

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.

share|improve this answer
    
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
show 6 more comments

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

Edit: Here is an implementation.

share|improve this answer
5  
Funny ^^ but hardly an answer to his question. –  Olivier Bégassat Oct 25 '11 at 17:07
1  
or get octave for free=) –  Kabumbus Oct 25 '11 at 17:38
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.