In a nutshell, given the singular decomposition of a matrix *A*,

the Moore-Penrose pseudoinverse is given by

This post will explain what the terms above mean, and how to compute them in Python and in Mathematica.

## Singular Value Decomposition (SVD)

The singular value decomposition of a matrix is a sort of change of coordinates that makes the matrix simple, a generalization of diagonalization.

### Matrix diagonalization

If a square matrix *A* is diagonalizable, then there is a matrix *P* such that

where the matrix *D* is diagonal. You could think of *P* as a change of coordinates that makes the action of *A* as simple as possible. The elements on the diagonal of *D* are the eigenvalues of *A* and the columns of *P* are the corresponding eigenvectors.

Unfortunately not all matrices can be diagonalized. Singular value decomposition is a way to do something like diagonalization for *any* matrix, even non-square matrices.

### Generalization to SVD

Singular value decomposition generalizes diagonalization. The matrix Σ in SVD is analogous to *D* in diagonalization. Σ is diagonal, though it may not be square. The matrices on either side of Σ are analogous to the matrix *P* in diagonalization, though now there are two different matrices, and they are not necessarily inverses of each other. The matrices *U* and *V* are square, but not necessarily of the same dimension.

The elements along the diagonal of Σ are not necessarily eigenvalues but **singular values**, which are a generalization of eigenvalues. Similarly the columns of *U* and *V* are not necessarily eigenvectors but **left singular vectors** and **right singular vectors** respectively.

The star superscript indicates **conjugate transpose**. If a matrix has all real components, then the conjugate transpose is just the transpose. But if the matrix has complex entries, you take the conjugate and transpose each entry.

The matrices *U* and *V* are **unitary**. A matrix *M* is unitary if its inverse is its conjugate transpose, i.e. *M*^{*} *M* = *M**M*^{*} = *I*.

## Pseudoinverse and SVD

The (Moore-Penrose) **pseudoinverse** of a matrix generalizes the notion of an inverse, somewhat like the way SVD generalized diagonalization. Not every matrix has an inverse, but every matrix has a pseudoinverse, even non-square matrices.

Computing the pseudoinverse from the SVD is simple.

If

then

where Σ^{+} is formed from Σ by taking the reciprocal of all the non-zero elements, leaving all the zeros alone, and making the matrix the right shape: if Σ is an *m* by *n* matrix, then Σ^{+} must be an *n* by *m* matrix.

We’ll give examples below in Mathematica and Python.

### Computing SVD in Mathematica

Let’s start with the matrix *A* below.

We can find the SVD of *A* with the following Mathematica commands.

A = {{2, -1, 0}, {4, 3, -2}} {U, S, V} = SingularValueDecomposition[A]

From this we learn that the singular value decomposition of *A* is

Note that the last matrix is not *V* but the transpose of *V*. Mathematica returns *V* itself, not its transpose.

If we multiply the matrices back together we can verify that we get *A* back.

U . S. Transpose[V]

This returns

{{2, -1, 0}, {4, 3, -2}}

as we’d expect.

### Computing pseudoinverse in Mathematica

The Mathematica command for computing the pseudoinverse is simply `PseudoInverse`

. (The best thing about Mathematica is it’s consistent, predictable naming.)

PseudoInverse[A]

This returns

{{19/60, 1/12}, {-(11/30), 1/6}, {1/12, -(1/12)}}

And we can confirm that computing the pseudoinverse via the SVD

Sp = {{1/Sqrt[30], 0}, {0, 1/2}, {0, 0}} V . Sp . Transpose[U]

gives the same result.

### Computing SVD in Python

Next we compute the singular value decomposition in Python (NumPy).

>>> a = np.matrix([[2, -1, 0],[4,3,-2]]) >>> u, s, vt = np.linalg.svd(a, full_matrices=True)

Note that `np.linalg.svd`

returns the *transpose* of *V*, not the *V* in the definition of singular value decomposition.

Also, the object `s`

is not the diagonal matrix Σ but a vector containing only the diagonal elements, i.e. just the singular values. This can save a lot of space if the matrix is large. The NumPy method `svd`

has other efficiency-related options that I won’t go into here.

We can verify that the SVD is correct by turning `s`

back into a matrix and multiply the components together.

>>> ss = np.matrix([[s[0], 0, 0], [0, s[1], 0]]) >>> u*ss*vt

This returns the matrix *A*, within floating point accuracy. Since Python is doing floating point computations, not symbolic calculation like Mathematica, the zero in *A* turns into `-3.8e-16`

.

Note that the singular value decompositions as computed by Mathematica and Python differ in a few signs here and there; the SVD is not unique.

### Computing pseudoinverse in Python

The pseudoinverse can be computed in NumPy with `np.linalg.pinv`

.

>>> np.linalg.pinv(a) matrix([[ 0.31666667, 0.08333333], [-0.36666667, 0.16666667], [ 0.08333333, -0.08333333]])

This returns the same result as Mathematica above, up to floating point precision.