# Calculating the Moore-Penrose pseudoinverse

I have a problem with a project requiring me to calculate the Moore-Penrose pseudo inverse. I've also posted about this on StackOverflow, where you can see my progress.

From what I understand from Planet Math you can simply compute the pseudoinverse only the first formula which I can understand, but it also says that this is for general cases, and you have to do SVD (singular value decomposition) and the formula becomes much complicated (the second formula) which I don't understand... I mean,

1. What is V? What is S? I'm confused.
2. How can I calculate SVD?
3. Can you please help me building the code/algorithm, or just a simple advice?
4. Why there are two pseudo inverse formulas?

Left pseudo inverse formula $$A_\text{left} = (A^TA)^{-1}A^T$$ Right pseudo inverse formula $$A_\text{right}=A^T(AA^T)^{-1}$$

Thank you very much, Daniel.

These formulas are for different matrix formats of the rectangular matrix $A$.

The matrix to be (pseudo-)inverted should have full rank. (added:) If $A\in I\!\!R^{m\times n}$ is a tall matrix, $m>n$, then this means $rank(A)=n$, that is, the columns have to be linearly independent, or $A$ as a linear map has to be injective. If $A$ is a wide matrix, $m<n$, then the rows of the matrix have to be independent to give full rank. (/edit)

If full rank is a given, then you are better off simplifying these formulas using a QR decomposition for $A$ resp. $A^T$. There the R factor is square and $Q$ is a narrow tall matrix with the same format as $A$ or $A^T$,

If $A$ is tall, then $A=QR$ and $A^{\oplus}_{left}=R^{-1}Q^T$

If $A$ is wide, then $A^T=QR$, $A=R^TQ^T$, and $A^{\oplus}_{right}=QR^{-T}$.

You only need an SVD if $A$ is suspected to not have the maximal rank for its format. Then a reliable rank estimation is only possible comparing the magnitudes of the singular values of $A$. The difference is $A^{\oplus}$ having a very large number or a zero as a singular value where $A$ has a very small singular value.

Added, since wikipedia is curiosly silent about this: Numerically, you first compute or let a library compute the SVD $A=U\Sigma V^T$ where $Σ=diag(σ_1,σ_2,\dots,σ_r)$ is the diagonal matrix of singular values, ordered in decreasing size $σ_1\ge σ_2\ge\dots\ge σ_r$.

Then you estimate the effective rank by looking for the smallest $k$ with for instance $σ_{k+1}<10^{-8}σ_1$ or as another strategy, $σ_{k+1}<10^{-2}σ_k$, or a combination of both. The factors defining what is "small enough" are a matter of taste and experience.

With this estimated effective rank $k$ you compute $$Σ^⊕=diag(σ_1^{-1},σ_2^{-1},\dots,σ_k^{-1},0,\dots,0)$$ and $$A^⊕=VΣ^⊕U^T.$$

Note how the singular values in $Σ^⊕$ and thus $A^⊕$ are increasing in this form, that is, truncating at the effective rank is a very sensitive operation, differences in this estimation lead to wildly varying results for the pseudo-inverse.

• 1-> what do you mean by different matrix formats ?? 2-> What if the matrix dosen't have full rank? What happens, i have to use QR decomposition as you sugested? 3-> ok, but don't SVD help me calculating the pseudo inverse? at least for every case 4-> can't i just use the first formula from A-left formula? And if not, why? 5-> sorry, i am not a math guru, i understand the baiscs of matrix, i also understand how to calculate the inverse, but this is a little bit tricky Jan 1, 2014 at 22:41
• 1) more rows than columns->tall, more columns than rows->wide -- 2) If you know the structural reason for the rank deficit, you better treat it by reducing redundancy in the problem formulation. If the rank deficit happens dynamically, you better use the better analytic power of the SVD. Jan 1, 2014 at 23:01
• -- 3) Yes, but SVD is still an iterative process with time $O(n^3|\log ε|)$ to reach a certain precision $ε$. The initial reduction to bidiagonal form corresponds to 2 QR decompositions. Meaning, if QR works, it is much faster. -- 4) You can, but for large matrices, the matrix product makes the condition quadratically worse. -- 5) try out Householder reflectors and Givens rotations, the mathematics is not that bad, keeping track of the indices for the positions to modify however... Anyway, the proposed library has all these fancy decompositions, so you only need to put the results together. Jan 1, 2014 at 23:02
• ** 3-> ** so, correct me if i understand, First formula A-left > SVD > QR decomposition ** 4-> ** so First formula A-left might work for 3x3, 4x4, but the higher you go, you might find bigger errors, is that right? ** 5-> ** i know they are better algorithms, put for starters let me try First formula A-left to see if it works, i can test it here calculator-fx.com/calculator/linear-algebra/… ** extra-> ** thank you very much for answering, i had never been into mathematics this deep and coding before ... Jan 2, 2014 at 0:03
• * 3) The computational effort for the 'first formula' is less than that for the SVD, numerical stability is worse. And, as said, if it works, QR is faster than SVD, the numerical results will be about the same. * 4) Yes, it will also work for 10x10, but for 1000x1000 I would expect problems. * 5) By all means, try out all methods. Jan 2, 2014 at 7:22

While the SVD yields a "clean" way to construct the pseudoinverse, it is sometimes an "overkill" in terms of efficiency.

The Moore-Penrose pseudoinverse can be seen as follows: Let $\ell:\mathbb R^n\rightarrow\mathbb R^m$ be a linear map. Then $\ell$ induces an isomorphism $\ell':{\rm Ker}(\ell)^\perp\rightarrow {\rm Im}(\ell)$. Then the Moore-Penrose pseudoinverse $\ell^+:\mathbb R^m\rightarrow \mathbb R^n$ can be described as follows.

$$\ell^+(x)=\ell'^{-1}(\Pi(x)),$$ where $\Pi$ is the orthogonal projection of $x$ on ${\rm Im}(\ell)$.

In other words, what you need is to compute orthonormal bases of ${\rm Im}(\ell)$ and of ${\rm Ker}(\ell)^\perp$ to contruct the Moore-Penrose pseudoinverse.

For an algorithm, you may be interested by the iterative method here

edit: roughly speaking, one way to see why the SVD might be an "overkill" is that if $A$ is a matrix with rational coefficients, then $A^+$ also have rational coefficients (see e.g. this paper), while the entries of the SVD are algebraic numbers.

• 1-> what do you mean about "overkill"? like is process consuming? first i just want to see this working, then i will think about optimisations 2-> I don't quite understand this formula ℓ+(x)=ℓ′−1(Π(x)) i'm used to work with formulas like uppwards 3-> ok, but don't SVD help me calculating the pseudo inverse? at least for every case 4-> sorry, i am not a math guru, i understand the baiscs of matrix, i also understand how to calculate the inverse, but this is a little bit tricky 5-> can't i just use the first formula from A-left formula? And if not, why? Jan 1, 2014 at 22:49
• 1-> By "overkill", I wanted to say that the SVD is a very powerful tool, maybe too powerful for computing the Moore-Penrose pseudoinverse, but there is indeed a very nice relation, see en.wikipedia.org/wiki/… 2->It is just a notation: a matrix actually represents a linear map (the map $u\mapsto A\cdot u$) and reasoning in terms of maps is often convenient. 3->Yes, but the SVD involves computing with algebraic numbers, and often makes exact computations impossible.
– emeu
Jan 2, 2014 at 4:07
• 5-> A-left only works if the linear map $u\mapsto A\cdot u$ is surjective. A-right only works if the map is injective. (note that if $A$ is a square invertible matrix, then both formulas give the same result: $A^{-1}$)
– emeu
Jan 2, 2014 at 4:08
• 2-> for me, a matrix is just a table that i can referr to its elements like A[i][j], thats the way i'm seeing it, because i worked a lot with arrays 5-> yes, i observed that, A-left gives same result as the inverse, and i cannot understand a something, given an array (1,2,3 - 4,5,6 - 7 8 9) that has the determinant 0 gives me a NULL matrix, see here (calculator-fx.com/calculator/linear-algebra/matrix-determinant) ... why is that? calculator-fx.com/calculator/linear-algebra/… gives me the right answer. Jan 2, 2014 at 18:27
• 2-> Yes, seeing a matrix as a table is also a very good way to think of it, especially for computational purposes. However, it is often a good idea (if you want to study this topic deeper) to have both representations in mind: as a table and as a linear map. 5-> A-left only works if the associated linear map is surjective, which is not the case if the determinant of a square matrix $A$ is zero. Consequently, in that case, the A-left formula does not work (but the Moore-Penrose pseudoinverse is still well defined and there are other ways to compute it)
– emeu
Jan 3, 2014 at 7:32