perform pca via eigendecomposition

pca_eigh.py

def dual_pca(X, max_num_pcs):
    """
    assuming N << D (X.shape[0] << X.shape[1]), 
    we first compute the kernel matrix K = XX^T, perform an
    eigendecomposition and subsequently reconstruct
    truncated principle components
    """
    
    
    K = X.dot(X.T)
    
    S, U = np.linalg.eigh(K)
    S = S[::-1]
    U = U[:,::-1]
    
    # reduce dimensionality
    U = U[:,0:max_num_pcs]
    S = S[0:max_num_pcs]

    # reconstruct principle components
    PCs = - U * np.sqrt(S)

    # compare to scikits (which uses SVD)
    from sklearn.decomposition import PCA
    pca = PCA(n_components=max_num_pcs)
    PC2 = pca.fit_transform(G)
    
    # PCs and PC2 will be the same except for differences in sign