Nice properties about SVD

svd.py

import numpy as np
import matplotlib.pyplot as plt

# Relation of SVD to PCA and eigen-problems
# A = USV'
# A'A = VSU'USV' = VS^2V'
# A'AV = VS^2V'V
# A'AV = VS^2 
# which is an eigenvector problem. Means V are the eigenvectors of A'A.
# A similar argument leads to U being the eigenvectors AA'.

def sorted_eig(X):
    w,e = np.linalg.eig(X)
    idx = w.argsort()[::-1]
    return w[idx], e[:,idx]

def svd(X):
    u,s,v = np.linalg.svd(A, full_matrices=0)
    return u,s,v.T

def draw_frame(ax, o, v, s, ec='red'):
    ss = s / np.sum(s)
    plt.arrow(o[0], o[1], ss[0]*v[0,0], ss[0]*v[1,0], head_width=0.1, head_length=0.1, ec=ec)
    plt.arrow(o[0], o[1], ss[1]*v[0,1], ss[1]*v[1,1], head_width=0.1, head_length=0.1, ec=ec)

A = np.random.multivariate_normal([10,10], cov=[[1, -1.5],[-1.5, 3]], size=200)
mu = np.mean(A, axis=0)
A = A - mu[None, :]

u,s,v = svd(A)
w,e = sorted_eig(A.T @ A)
w = np.sqrt(w)

plt.scatter(A[:,0], A[:, 1])
draw_frame(plt, [0., 0.], v, s, 'red')
draw_frame(plt, [0., 0.], e, w, 'green')

plt.axes().set_aspect('equal', 'datalim')
plt.show()