# Computing the SVD¶

In [1]:
import numpy as np
import numpy.linalg as la

In [2]:
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):

In [11]:
eigvals, eigvecs = la.eigh(A.T @ A)

In [12]:
eigvals

Out[12]:
array([0.08637178, 0.457892  , 2.04177547, 2.34383161, 8.37000184])

Eigenvalues are real and positive. Coincidence?

In [13]:
eigvecs.shape

Out[13]:
(5, 5)

Check that those are in fact eigenvectors and eigenvalues:

In [14]:
B = A.T @ A
B - eigvecs @ np.diag(eigvals) @ la.inv(eigvecs)

Out[14]:
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:

In [15]:
la.norm(eigvecs.T @ eigvecs  - np.eye(n))

Out[15]:
2.094158390904296e-15

Now piece together the SVD:

In [16]:
Sigma = np.diag(np.sqrt(eigvals))

In [17]:
V = eigvecs

In [18]:
U = A @ V @ la.inv(Sigma)


Check orthogonality of U:

In [19]:
U @ U.T - np.eye(n)

Out[19]:
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]])