Computing the pseudoinverse of a 3x3 matrix

Question

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.

Answer 1

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

Answer 2

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
!=======================================================================
   function adj(a)
!
!  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)
   inv = adj(a)*detr
   end function

I use double precision for reals:

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