computing principal components and singular values

pca.py

import numpy as np
import numpy.linalg as la

np.random.seed(42)
X = np.random.randn(100, 5)

assert la.matrix_rank(X) == 5


def pca_svd(X):
    # svd, randomized svd
    _, s, V = la.svd(X - X.mean(axis=0))
    return V.T, s


def pca_nipals(X, n=2, eps=1e-10):
    E = X - X.mean(axis=0)
    ps = []
    n = min(X.shape[1], n)
    for i in range(n):
        t = E[:, 0]
        while True:
            # Might not be orthogonal if X doesn't have full rank.
            p = E.T.dot(t)
            p /= la.norm(p, ord=2)
            t_new = E.dot(p)
            if np.abs(t.dot(t) - t_new.dot(t_new)) < eps:
                E = E - t_new[:, None].dot(p[None, :])
                ps.append(p)
                break
            t = t_new

    ps = np.vstack(ps).T
    return ps, X.dot(ps * np.sqrt(X.shape[0])).std(axis=0)


def pca_eig(X):
    """
    cov(X) = X_.dot(X_.T)/(X_.shape[1]-1)
    """
    u, V = la.eig(np.cov((X - X.mean(axis=0)).T, ddof=0) * X.shape[0])
    s = np.sqrt(u)  # Assumes X has full rank.
    sort = np.argsort(s)[::-1]
    return V[:, sort], s[sort]


def pca_eig2(X):
    X_ = X - X.mean(axis=0)
    u, V = la.eig(X_.T.dot(X_))
    s = np.sqrt(u)  # Assumes X has full rank.
    sort = np.argsort(s)[::-1]
    return V[:, sort], s[sort]


pc, s = pca_svd(X)
pc2, s2 = pca_nipals(X, n=5)
pc3, s3 = pca_eig(X)
pc4, s4 = pca_eig2(X)

print(s)
print(s2)
print(s3)
print(s4)

print(pc)
print(pc2)
print(pc3)
print(pc4)