SVD & PCA

本文讲解SVD和PCA

Singular Value Decomposition

abbr. SVD

The Singular Value Decomposition is a highlight of linear algebra. $A$ is any $m$ by $n$ matrix, square or rectangular. Its rank is $r$. We will diagonalize this $A$, but not by $S^{-1} AS$. The eigenvectors in $S$ have three big problems: They are usually not orthogonal, there are not always enough eigenvectors, and $Ax = \lambda x$ requires $A$ to be square. The singular vectors of $A$ solve all those problems in a perfect way.

” $A$ is diagonalized”

$Av_1=\sigma_1u_1\qquad Av_2=\sigma_2u_2\qquad \dots \qquad Av_r=\sigma_ru_r$

The singular vectors $v_1 \dots v_r$ are in the row space of A. The singular values $\sigma_1 \dots \sigma_r$ are all positives number. When the $v's$ and $u's$ go into the columns of $V$ and $U$, orthogonality gives $V^TV=I$ and $U^TU=I$. The $\sigma's$ go into diagonal matrix $\Sigma$.

proof:
$v$是$A^TA$的eigen vector
那么$v$在$A$的row space上
$A^{T}Av_i=\sigma_i^2v_i$
$A^{T}A$一定是半正定的($X^TA^TAX\geq 0$)
所以存在$r$个eigenvalue 大于0 (r为rank)

row space $\perp$ nullspace
(proof: row space是$A^Ty$的集合与$nullspace$中的$x$的dot product
$X^T(A^Ty)=(AX)^Ty=0$)

$V$在$R^n$上，可以找出$r$个正交基($v's$是singular vector)，可以在nullspace上找出$n-r$个正交基，可以构造$V_{n\times n}$，显然$V^T=V$

目标是$AV = U\Sigma$
$V^T=V^{-1}$
$A=U\Sigma V^T$
$A^TAv_i=\sigma_i^2 v_i$
$v_i^TA^TAv_i=\sigma_i^2 v_i^Tv_i$
$\|Av_i\|^2=\sigma^2_i$ so that $\|Av_i\|=\sigma_i$
$AA^TAv_i=\sigma^2_iAv_i$ gives $u_i=Av_i/\sigma_i$ as a unit eigenvector of $AA^T$

$u_i$在A的column space上，所以存在$r$个正交基，与$v's$对应，同理可以在left null space上找到m-r个正交基组成$U$。

$A=U\Sigma V^T=u_1\sigma_1v_1^T+\dots+u_r\sigma_rv_r^T$

将SVD用于降维：
$A_{m\times n}V_{n\times r}=U_{m\times r}\Sigma_{r\times r}$
就是将$A$投影到$V$所以定义的space上，从而达到了降维。可以按照singular value的大小进行排列，如果前$k$个singular value的大小占据了90%（只是个例子），那么相应的$v_i$就可以丢掉，从而继续降维。eigenvalue相当于把数据投影到eigenvector上后，所有点的方差（因为要构造unit vector）。

Principal Component Analysis

abbr. PCA
SVD是PCA的一种方法，还有一种常用的方法就是将SVD降维中所用到的$V$换成了$A$的协方差矩阵的eigenvalue所组成的投影矩阵。
PCA的目的是使投影后的矩阵各向量之间的协方差为0，既去除了各个向量之间的线性关系，并且使得投影后的矩阵的在各个方向上的方差达到最大化。可以参考主成分分分析这篇文章所讲的最大方差理论。简单的对于分类来说，方差越大那么区分度越高。

算法步骤：

1. 求矩阵$A_{m\times n}$中的每一列的均值，可以得到一个$n\times 1$的vector $m$。其中A的每一行代表一个数据点，而每一列就是一个feature；
2. 对$A$的每一行都减去$m$；
3. 求$A$的转置的协方差矩阵$C$；
4. 求$C$的eigenvalue和eigenvector；
5. 根据eigenvalue对eigenvector进行排序，可根据需求舍弃较小的eigenvalue所对应的eigenvector；
6. 将eigenvector所组成的矩阵进行转置($r\times n$)与转置的$A$($n\times m$)进行dot product，将结果转置既为所求结果。

因为用到了协方差矩阵，协方差矩阵中每一项度量的是两个随机变量的相关性，$\Sigma_{ij}=cov(X_i,Xj)=E[(X_i-u_i)(X_j-u_j)^T]$，协方差定义中的$X_i$是行向量，对应$A_{m\times n}$中的一列，既对一个random variable的一组观测值。所以需要对$A$进行转置，使得A的行向量成为feature。$cov(A^T)$的对角线是每一个feature的方差，其他的则是两个feature的协方差。
假设投影后的矩阵为$Y$，那么$Y^T=P^TA^T$，这里的表述方法是为了方便表达协方差矩阵。并且$P^T$是一个rotation matrix ，因为$P^T=P^{-1}$。
那么$cov(Y^T)=E((Y^T)(Y^T)^T)$
$=E((P^TA^T)(P^TA^T)^T)$
$=E(P^T(A^T(A^T)^T)P)$
$=P^Tcov(A^T)P$
由于$cov(A^T)$是一个对称矩阵，所以可以将其对角化，而对角化所用的矩阵就是$cov(A^T)的$eigenvetor所组成的矩阵(spectral theorem)。所以$p_i$就是$cov(A^T)$的eigenvector。在推导的过程中，要注意A是已经减去了均值的，要不不能这么推，这个过程反推也是有意义的，也能说明第一步为什么不需要减去均值。
接着就可以按照eigenvalue对eigenvector进行排序，根据需要cover多少variance来选几个eigenvector了。为什么这么选，可以想一下coviance matrix对角线上的值的意义。

import numpy as np

def pca(data,normalise=1):

    # centre data
    m = np.mean(data, axis=0)
    data -= m

    # Covariance matrix
    C = np.cov(np.transpose(data))

    # Compute eigenvalues and sort into descending order
    evals, evecs = np.linalg.eig(C)
    indices = np.argsort(evals)
    indices = indices[::-1]
    evecs = evecs[:, indices]
    evals = evals[indices]


    evalsNorm = evals / np.sum(evals)
    cumEvals = np.cumsum(evalsNorm)
    evecs = evecs[:,:(np.where(cumEvals > 0.9)[0][0] + 1)]

    if normalise:
       for i in range(np.shape(evecs)[1]):
           evecs[:, i] = evecs[:, i] /  np.linalg.norm(evecs[:, i])

    # Produce the new data matrix
    # x = np.transpose(np.dot(np.transpose(evecs), np.transpose(data)))
    x = np.dot(data, evecs)

    # Compute the original data again
    y = np.transpose(np.dot(evecs, np.transpose(x))) + m

    return x,y,evals,evecs

iris = np.genfromtxt(r"/home/real/pca/iris.data", delimiter=",", dtype="str")
iris = iris[:, range(np.shape(iris[0])[0] - 1)]
iris = iris.astype(np.float)
x,y,evals,evecs = pca(iris)
np.savetxt("result.csv", x, delimiter=",")
np.savetxt("resulty.csv", y, delimiter=",")