# Computing the pseudoinverse of a 3x3 matrix

I need to compute the (Moore-Penrose) pseudoinverse of fixed-size 3x3 matrices. I would prefer to have a simple method without bringing in the industrial strength machinery of Lapack. Are there any simple methods to do this?

Edit: I should add two notes: first, I am prioritizing robustness and simplicity over speed. Second, I really only need the action of the pseudo-inverse on a 3x1 vector.

-

The book Generalized Inverses: Theory and Applications, Second Edition, by Ben-Israel and Greville states the following (Chapter 1, Theorem 5, page 48):

Theorem: If $A$ is an $m$ by $n$ complex matrix of rank $r > 0$, and has the full rank factorization

\begin{align} A = XY, \end{align}

then its Moore-Penrose ($\{1,2,3,4\}$-)generalized inverse, $A^{\dagger}$ is given by

\begin{align} A^{\dagger} = Y^{*}(X^{*}AY^{*})^{-1}X^{*}, \end{align}

where the asterisk denotes Hermitian transpose.

Like many theoretical results, this result is amenable to proof, but not immediately useful for computation.

In numerical analysis, there's no routine for "full rank factorization". Rather, you get that information using the SVD or a QR factorization. Consequently, you are much better off using LAPACK to calculate one of those two well-known factorizations.

Suppose that the SVD of $A$ is $A = U\Sigma V^{*}$. Suppose that $\Sigma = \mathrm{diag}(\sigma_{1}, \ldots, \sigma_{n})$. Define $\Sigma_{\varepsilon}^{-} = \mathrm{diag}(\sigma_{1}^{-}, \ldots, \sigma_{n}^{-})$, where

\begin{align} \sigma_{i}^{-} = \left\{\begin{array}{cc}\sigma_{i}^{-1}, & \textrm{if $|\sigma_{i}| \geq \varepsilon$,} \\ 0, & \textrm{if $|\sigma_{i}| < \varepsilon$,}\end{array}\right. \end{align}

where singular values with magnitude less than $\varepsilon > 0$ are treated as "numerically zero".

Then

\begin{align} A^{\dagger} = V\Sigma_{\varepsilon}^{-}U^{*}. \end{align}

This method is not simple, but it will be more robust, and it will be correct.

Edit: Two common formulas for the Moore-Penrose generalized inverse are given on Wikipedia.

1. If $A$ is an $m$ by $n$ matrix with linearly independent columns (so $m \geq n$), then \begin{align} A^{\dagger} = (A^{*}A)^{-1}A^{*}. \end{align}
2. If $A$ is an $m$ by $n$ matrix with linearly independent rows (so $m \leq n$), then

\begin{align} A^{\dagger} = A^{*}(AA^{*})^{-1}. \end{align}

A problem with both of these formula is the requirement of linear independence: in the square case, if all columns are linearly independent, then all rows are also linearly independent, and the matrix is invertible. The inverse of a matrix is also its Moore-Penrose generalized inverse, and if the matrix inverse is truly needed for a computation (often times -- but not always -- it is not, and can be replaced by solving an appropriate linear system) it should be calculated directly using LU, QR, or SVD, rather than using one of the two formulas for Moore-Penrose generalized inverses.

Another problem with both of the above formulas in the square case, assuming that the conditions are satisfied, is that they will be less accurate numerically than an SVD- or QR-based approach, because the condition number of $AA^{*}$ or $A^{*}A$ will be the square of the condition number of $A$ (all condition numbers mentioned in this explanation are for linear systems). In the case where $A$ is not square, but one of the two linear independence conditions is satisfied, the product of $A$ with its Hermitian transpose (in the appropriate order) will have a condition number equal to the square of the ratio of the largest singular value of $A$ to the smallest nonzero singular value of $A$. Again, this condition number could be large. Such a matrix product is best avoided unless there are compelling reasons for using it in an application (speed, for instance).

-
Since I only need the action of the MP-inverse, can I just use Lapack's ZGELSY least norm solution (essentially, rank revealing QR, followed by rank determination, RZ factorization, and then QR-solve) on the vector I'm trying to apply it to? –  Victor Liu Feb 20 '13 at 20:03
It sounds plausible; Moore-Penrose generalized inverses are used in overdetermined least squares problems. I'd have to double-check for sure. The rank-revealing QR is definitely on the right track; the basic problem becomes "how do I handle a singular $R$ matrix?", and as long as that's done sensibly, then you should be okay. –  Geoff Oxberry Feb 21 '13 at 4:20

For 3x3 matrices this should be easy. Here's a least-squares routine in Fortran that forms Moore-Penrose pseudoinverse in the process of solution, it might help:

   function solve_leastsq(A,b,m,n) result(x)
!
!  Solve system with m x n system matrix in least square sense (minimizing Euclidean norm).
!  System is overdetermined so we solve A'A * x = A'b, where A' is transpose of A.
!  The solution is given by x = (A'A)^(-1)A'*b
!

!  Result
real(dp), dimension(n) :: x

!  Input
real(dp), dimension(m,n), intent(in) :: A
real(dp), dimension(m), intent(in) :: b

!  Locals
real(dp), dimension(n,m) ::  At
real(dp), dimension(n,n) ::  AtA, invAtA
real(dp), dimension(n) ::  Atb

At = transpose(A)      ! transpose of A : A'
AtA = matmul(At,A)     ! left-hand side multiplication of A by A': A'A
Atb = matmul(At,b)     ! A'b
invAtA = inv(AtA)      ! an inverse of A'A
x = matmul(invAtA,Atb) ! x = (A'A)^(-1) * A'b

end function


For this to work you need these routines to find an inverse of a 3x3 matrix:

   function det_3x3(a)
!  Result
real(dp) :: det_3x3
!  Input
real(dp), dimension(3,3), intent(in) :: a

det_3x3 =  a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2) &
-a(1,3)*a(2,2)*a(3,1)-a(1,2)*a(2,1)*a(3,3)-a(1,1)*a(2,3)*a(3,2)

end function
!=======================================================================
!
!  Adjugate of a 3x3 matrix (the transpose of the cofactor matrix).
!

!  Result
real(dp), dimension(3,3) :: adj
!  Input
real(dp), dimension(3,3), intent(in) :: a

adj(1,1) = a(2,2)*a(3,3) - a(2,3)*a(3,2)
adj(1,2) = a(1,3)*a(3,2) - a(1,2)*a(3,3)
adj(1,3) = a(1,2)*a(2,3) - a(1,3)*a(2,2)
adj(2,1) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
adj(2,2) = a(1,1)*a(3,3) - a(1,3)*a(3,1)
adj(2,3) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
adj(3,1) = a(2,1)*a(3,2) - a(2,2)*a(3,1)
adj(3,2) = a(1,2)*a(3,1) - a(1,1)*a(3,2)
adj(3,3) = a(1,1)*a(2,2) - a(1,2)*a(2,1)

end function
!=======================================================================
function inv(a)
!
!  An inverse of a 3x3 matrix.
!

!  Result
real(dp), dimension(3,3) :: inv
!  Input
real(dp), dimension(3,3), intent(in) :: a
!  Locals
real(dp) :: detr
detr = 1./det_3x3(a)

   integer, parameter :: dp = kind(1.0d0)

Your program calculates $A^{\dagger}$ using $A^{\dagger} = (A^{*}A)^{-1}A$ and leaves out a critical hypothesis for the validity of the formula: linear independence of the columns. In the event that the hypothesis holds for the 3 by 3 case, the matrix inverse can be used directly instead. In the event that the hypothesis does not hold, the formula will return an incorrect result. –  Geoff Oxberry Feb 20 '13 at 17:15