Singular value decomposition

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)