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]])