Singular value decomposition

The NumPy method numpy.linalg.svd returns the singular value decomposition (SVD) of a matrix, $\mathbf{M}$ , as the arrays $\mathbf{U}$ , $\mathbf{\Sigma}$ and $\mathbf{V}$ satisfying the factorization: $\mathbf{M} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\dagger$ where $\dagger$ denotes the Hermitian conjugate (the conjugate transpose).

The SVD and the eigendecomposition are related in that the left-singular row vectors, $\mathbf{U}$ are the eigenvectors of $\mathbf{MM^*}$ and the right-singular column vectors, $\mathbf{V}$ , are the eigenvectors of $\mathbf{M^*M}$ . Furthermore, the diagonal entries of $\mathbf{\Sigma}$ are the square roots of the non-zero eigenvalues of both $\mathbf{MM^*}$ and $\mathbf{M^*M}$ .

Show that this is the case for the special case of $\mathbf{M}$ a $3 \times 3$ matrix with random real entries by comparing the output of numpy.linalg.svd with that of numpy.linalg.eig.

Hints: The singular values of $\mathbf{M}$ are sorted in descending order but the eigenvalues returned by numpy.linalg.eig are in no particular order. Both methods produce normalized eigenvectors, but may differ by sign (ignore the possibility that any of the eigenvalues could have an eigenspace with dimension greater than 1).

Here we compare the eigenvalues and eigenvectors of $\mathbf{MM^*}$ and $\mathbf{M^*M}$ obtained by singular value decomposition with np.linalg.svd and by np.linalg.eig.

svd

import numpy as np

np.random.seed(1)
M = np.random.random((3, 3))

# Singular value decomposition
U, S, V = np.linalg.svd(M)
# Square S to get eigenvalues of MM* (also the eigenvalues of M*M)
svd_evals = S * S


def get_sorted_eigs(A):
    """
    Use np.linalg.eig to find and return the eigenvalues and eigenvectors of
    square matrix, A. These are returned, sorted into _descending_ order by
    eigenvalue.
    """

    evals, evecs = np.linalg.eig(A)
    sorted_indexes = np.argsort(evals)[::-1]
    return evals[sorted_indexes], evecs[:, sorted_indexes]


print("MM* by numpy.linalg.eig:")
evals, evecs = get_sorted_eigs(np.dot(M, M.T))
print("Compare eigenvalues:")
print(evals, svd_evals)
print("Eigenvalues equal: ", np.allclose(evals, svd_evals))
print("Compare eigenvectors:")
print(evecs, U, sep="\n")
# To compare eigenvectors verify that they are parallel or anti-parallel
print(
    "Eigenvectors equivalent: ", np.allclose(abs((evecs * U).sum(axis=0)), 1)
)
print("-" * 70)

print("M*M by numpy.linalg.eig:")
evals, evecs = get_sorted_eigs(np.dot(M.T, M))
print("Compare eigenvalues:")
print(evals, svd_evals)
print("Eigenvalues equal: ", np.allclose(evals, svd_evals))
print("Compare eigenvectors:")
# Take the transpose of V because we want these eigenvectors in rows
print(evecs, V.T, sep="\n")
print(
    "Eigenvectors equivalent: ",
    np.allclose(abs((evecs * V.T).sum(axis=0)), 1),
)
print("-" * 70)

Note that because valid eigenvectors may differ by a sign, we must allow for this in comparing the arrays produced by the two methods.

The output is:

svd-output

MM* by numpy.linalg.eig:
Compare eigenvalues:
[0.96809754 0.12478566 0.03288974] [0.96809754 0.12478566 0.03288974]
Eigenvalues equal:  True
Compare eigenvectors:
[[-0.8238009  -0.52374825  0.21688672]
 [-0.30647428  0.0896252  -0.94765017]
 [-0.4768916   0.84714526  0.23434869]]
[[-0.8238009   0.52374825  0.21688672]
 [-0.30647428 -0.0896252  -0.94765017]
 [-0.4768916  -0.84714526  0.23434869]]
Eigenvectors equivalent:  True
----------------------------------------------------------------------
M*M by numpy.linalg.eig:
Compare eigenvalues:
[0.96809754 0.12478566 0.03288974] [0.96809754 0.12478566 0.03288974]
Eigenvalues equal:  True
Compare eigenvectors:
[[ 0.53360687 -0.09491425  0.84038979]
 [ 0.81630245 -0.2020538  -0.54113268]
 [ 0.22116515  0.97476435 -0.03033856]]
[[-0.53360687  0.09491425 -0.84038979]
 [-0.81630245  0.2020538   0.54113268]
 [-0.22116515 -0.97476435  0.03033856]]
Eigenvectors equivalent:  True
----------------------------------------------------------------------

P6.5.3Singular value decomposition