Computing SVD and pseudoinverse

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

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 I in SVD is analogous to D in diagonalization. I is diagonal, though it may not be square. The matrices on either side of I 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 I 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 = MM^* = 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 I⁺ is formed from I by taking the reciprocal of all the non-zero elements, leaving all the zeros alone, and making the matrix the right shape: if I is an m by n matrix, then I⁺ 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 I 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.

Related posts

• Pseudoinverse is not an adjoint
• Don't invert that matrix
• Applied linear algebra