Computing the SVD¶

import numpy as np
import numpy.linalg as la

1) For a square matrix¶

m = 4
n = m
A = np.random.randn(m, n)
print(A)

[[ 1.42306309 -1.24203693 -0.93723067 -2.0891819 ]
 [-0.75566914  2.22185875 -0.72213972  1.72900854]
 [-0.53615049 -1.67399197 -0.94643061  0.32664051]
 [ 0.59359804  1.58629795  1.5037636  -0.98707451]]

Using numpy.linalg.svd¶

U, S, Vt = la.svd(A)

print(U)
print(U.shape)

[[-0.66398367 -0.2123579  -0.71524382  0.04955897]
 [ 0.68138956  0.1960955  -0.67713964  0.19681649]
 [-0.25888667  0.56114452  0.12748125  0.77578544]
 [ 0.16676826 -0.77560783  0.11686194  0.59746475]]
(4, 4)

print(Vt)
print(Vt.T.shape)

[[-0.30110577  0.74822109  0.15426018  0.5707051 ]
 [-0.39693586 -0.4816553  -0.53726629  0.5672698 ]
 [-0.33998961 -0.43358196  0.81740312  0.16812463]
 [-0.79761163  0.14205653 -0.13928701 -0.56941616]]
(4, 4)

print(S)
print(S.shape)

[4.05837523 3.05248106 1.48570324 0.17488016]
(4,)

Using eigen-decomposition¶

Now compute the eigenvalues and eigenvectors of $A^TA$ as eigvals and eigvecs

eigvals, eigvecs = la.eig(A.T.dot(A))

Eigenvalues are real and positive. Coincidence?

eigvals

array([16.47040952,  9.31764063,  0.03058307,  2.20731411])

eigvecs are orthonormal! Check:

eigvecs.T @ eigvecs

array([[ 1.00000000e+00,  3.88578059e-16,  0.00000000e+00,
        -1.38777878e-16],
       [ 3.88578059e-16,  1.00000000e+00, -2.77555756e-17,
        -3.33066907e-16],
       [ 0.00000000e+00, -2.77555756e-17,  1.00000000e+00,
         4.16333634e-16],
       [-1.38777878e-16, -3.33066907e-16,  4.16333634e-16,
         1.00000000e+00]])

Now piece together the SVD:

S2 = np.sqrt(eigvals)
V2 = eigvecs
U2 = A @ V2 @ la.inv(np.diag(S2))

2) For a nnon-square square matrix¶

m = 3
n = 5
A = np.random.randn(m, n)
print(A)

[[ 1.18363425 -0.34818317 -1.89959274  0.55446767 -1.59147727]
 [ 0.19794113 -1.95160646  0.27558341  0.93666225 -0.38931728]
 [ 1.45876852 -0.14058534  0.72291651  1.43042735  0.99827354]]

U, S, Vt = la.svd(A,full_matrices=True)

print(U)
print(U.shape)

print(Vt)
print(Vt.T.shape)

print(S)
print(S.shape)

[[-0.89303791  0.34321089  0.29101474]
 [-0.43033767 -0.46240657 -0.77523523]
 [-0.13150204 -0.81754905  0.56064317]]
(3, 3)
[[-0.45359036  0.3975683   0.50415236 -0.36937045  0.49557149]
 [-0.33549725  0.34312616 -0.5237104  -0.53970428 -0.45183265]
 [ 0.56937104  0.75220407 -0.2038258   0.13386068  0.22481552]
 [-0.5954303   0.25135987 -0.22752045  0.72785242  0.02731793]
 [-0.05455968 -0.30856302 -0.61501813 -0.15683409  0.70637658]]
(5, 5)
[2.941074   2.61673664 1.77187214]
(3,)

eigvals, eigvecs = la.eig(A.T.dot(A))
print(eigvals)

[ 3.13953087e+00  8.64991627e+00  6.84731066e+00  2.37949537e-16
 -7.70846804e-17]