Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
np.random.seed(15)
n = 5
A = np.random.randn(n, n)
Now compute the eigenvalues and eigenvectors of $A^TA$ as eigvals
and eigvecs
using la.eig
or la.eigh
(symmetric):
eigvals, eigvecs = la.eigh(A.T @ A)
eigvals
array([0.08637178, 0.457892 , 2.04177547, 2.34383161, 8.37000184])
Eigenvalues are real and positive. Coincidence?
eigvecs.shape
(5, 5)
Check that those are in fact eigenvectors and eigenvalues:
B = A.T @ A
B - eigvecs @ np.diag(eigvals) @ la.inv(eigvecs)
array([[ 1.77635684e-15, 2.22044605e-15, 3.10862447e-15, 8.88178420e-16, 8.88178420e-16], [ 2.22044605e-16, -8.88178420e-16, 0.00000000e+00, -8.88178420e-16, 2.63677968e-16], [-8.88178420e-16, 4.44089210e-16, 8.88178420e-16, -2.22044605e-16, 6.66133815e-16], [ 0.00000000e+00, -4.44089210e-16, 8.88178420e-16, 0.00000000e+00, 7.77156117e-16], [ 0.00000000e+00, -9.02056208e-17, 3.88578059e-16, 0.00000000e+00, 3.33066907e-16]])
eigvecs
are orthonormal! (Why?)
Check:
la.norm(eigvecs.T @ eigvecs - np.eye(n))
2.094158390904296e-15
Now piece together the SVD:
Sigma = np.diag(np.sqrt(eigvals))
V = eigvecs
U = A @ V @ la.inv(Sigma)
Check orthogonality of U
:
U @ U.T - np.eye(n)
array([[-1.11022302e-15, -7.67324533e-16, -6.67804308e-16, 1.47963782e-15, 1.44717415e-16], [-7.67324533e-16, 6.66133815e-16, -5.43672425e-16, 6.51481016e-16, -2.64998805e-16], [-6.67804308e-16, -5.43672425e-16, 2.66453526e-15, -1.42306617e-15, 1.55464267e-15], [ 1.47963782e-15, 6.51481016e-16, -1.42306617e-15, 3.10862447e-15, -1.13201575e-15], [ 1.44717415e-16, -2.64998805e-16, 1.55464267e-15, -1.13201575e-15, 4.44089210e-16]])