机器学习: PCA

前言：

此博客仅仅为了自己复习方便，表述很不严谨，读者要了解PCA，建议直接移步参考文献。

---

一、基本介绍

主成分分析(Principal components analysis)，以下简称PCA, 广泛地运用在数据压缩和噪声消除中，是一种很重要的无监督学习算法。

这一方法使用正交变换, 把由线性相关变量表示的观测数据转换为少数几个由线性无关变量表示的数据，线性无关的变量称为主成分。

PCA算法包括协方差矩阵的特征值分解和数据矩阵的奇异值分解方法。

PCA的目标是:尽可能保留信息地对数据实现降维。

二、方法

基于PCA的目的, 总体来说可以从两个方向进行优化:

(1): 能够最大可能地重构原始数据。参考[1]的小节2。
(2): 使得数据在低维空间差异尽可能大。差异程度用各坐标轴上的方差之和度量。参考[1]的小节3。

两种思路的优化结果完全等价的: 最后的优化目标:
$$arg\underset{W}{\max }tr(W^{T}X{X}^{T}W)s.t.W^{T}W=I$$
其中$W\in R^{K\times K^{\prime }},X\in R^{K\times N}$, 其中$K，K^{\prime }$分别为原空间和降维空间的维度。

值得注意的是，若从第(2)个思路出发, $Z=W^{T}X\in R^{K^{\prime }\times N}$是数据$X$在降维空间里的投影。很明显降维空间样本在各方向的方差和正比于$tr(Z{Z}^T)=tr(W^{T}X{X}^{T}W)$. 因此凭直觉就可以直接得出优化目标。

最后得到的最优$W\in R^{K\times K^{\prime }}$由样本协方差矩阵$S=X{X}^T$的前$K^{\prime }$大的特征值对应的特征向量组成。

准确来说$S=\frac{X{X}^T}{N-1}$, 但是一个常数因子不影响PCA结果，为了简便，下方认为$S=X{X}^T$

三、两个步骤

PCA降维实际上可以分为两个步骤。(1) 通过一个正交变换对原坐标轴进行旋转，(2) 然后去除掉旋转后数据坐标中的一些维度。这在[3]中有比较好的动画展示。

因此我们本质上来说是要寻找一个正交矩阵$W^{\prime }\in R^{K\times K}$， 对数据空间进行旋转操作，求出数据在新空间中的坐标: $W^{\prime T}X\in R^{K\times N}$, 最后舍弃一些维度。将两个步骤合起来，即进行操作: $W^{T}X\in R^{K^{\prime }\times N}$，从而实现降维。

---

对$S$的特征值进行降序排列: $({\lambda }_1,{\lambda }_2,...,{\lambda }_K)$

[2]中P301页定理16.1指出: 样本第$k$主成分${\mathbf{y}}_{\mathbf{k}}$ 对应着 $S$第$k$个特征向量，且$var({\mathbf{y}}_{\mathbf{k}})={\lambda }_k$, 因此这也揭示了特征值大小的物理意义: 衡量了数据在其对应特征向量方向上的信息量大小。

注：要满足$var({\mathbf{y}}_{\mathbf{k}})={\lambda }_k$ 需要$S=\frac{X{X}^T}{N-1}$

无论坐标轴如何旋转，由于个体到坐标轴原点的距离不变，归一化后的数据均值为0， 因此各坐标轴方差之和始终等于样本到原点距离平方之和，要使得方差最大，则选择对应方差大的主成分作为降维后的坐标轴。

四、两个实现方法

1、基于协方差矩阵

由目标推导, 基于样本协方差矩阵$S=X{X}^T$实现PCA算法过程如下:

---

输入:

• 样本集 $X\in R^{K\times N}$
• 低秩空间维度$K^{\prime }$

过程：

• (1) 对数据进行归一化：$x_{ij}=\frac{x_{ij}-{\bar{x}}_i}{\sqrt{s_i}}$. ($x_{ij}$表示第$j$个数据的第$i$个特征)
• (2) 计算样本协方差矩阵$S=X{X}^T$
• (3) 对$S$进行特征值分解
• (4) 取前$K^{\prime }$大特征值对应的特征向量、并单位化，得到$W=(w_1,w_2,...,w_{K^{\prime }})$
• (5) 计算主成分变量${\mathbf{y}}_{\mathbf{k}}$和原变量${\mathbf{x}}_i$的相关关系: $\rho ({\mathbf{y}}_{\mathbf{k}},{\mathbf{x}}_{\mathbf{i}})=\frac{\sqrt{{\lambda }_k}{w}_{ik}}{\sqrt{s_i}}$
• (6) 计算$k$个主成分变量对原始变量${\mathbf{x}}_i$的贡献率:$v_i={\rho }^2({\mathbf{x}}_i,({\mathbf{y}}_1,...,{\mathbf{y}}_{K^{\prime }}))={\sum }_{k=1}^{K^{\prime }}\frac{{\lambda }_k{w}_{ik}^2}{s_i}$. （${\mathbf{y}}_{\mathbf{i}}、{\mathbf{y}}_{\mathbf{j}}$不相关）
• (7) 计算降维后的数据: $Y=W^{T}X$

输出:

• 降维后的数据$Y\in R^{K^{\prime }\times N}$
• 主成分变量和原变量的相关性矩阵$\mathrm{P}\in R^{K^{\prime }\times K}$
• $k$个主成分变量对原始变量的贡献率向量: $\mathbf{v}\in R^K$

---

2、基于SVD分解

SVD分解原理课参考[4]. 对于样本集$X\in R^{K\times N}$, 我们可以对其进行SVD分解:
$$X=U\Sigma V^T$$
其中$U\in R^{K\times K}$, $\Sigma \in R^{K\times N}$, $V\in R^{N\times N}$. 且$U^{T}U=I$, $V^{T}V=I$.

协方差矩阵：
$$S=X{X}^T=U\Sigma V^{T}V{\Sigma }^T{U}^T=U{\Sigma }^2{U}^T$$
进而：
$$X^{T}XU=U{\Sigma }^2$$

因此$U\in R^{K\times K}$的列向量即是$S$的特征向量。特征值${\lambda }_k={\sigma }_{kk}^2$. 基于此，有基于SVD的PCA算法:

输入:

• 样本集 $X\in R^{K\times N}$
• 低秩空间维度$K^{\prime }$

过程：

• (1) 对数据进行归一化：$x_{ij}=\frac{x_{ij}-{\bar{x}}_i}{\sqrt{s_i}}$. ($x_{ij}$表示第$j$个数据的第$i$个特征)
• (2) 对$X$进行SVD分解，得到$X=U\Sigma V^T$
• (3) 按${\sigma }_{kk}$从大到小顺序从$U$中选取$K^{\prime }$个列向量.
• (4) 对选取的向量进行单位化，得到$W=(u_1,u_2,...,u_{K^{\prime }})$
• (5) 计算主成分变量${\mathbf{y}}_{\mathbf{k}}$和原变量${\mathbf{x}}_i$的相关关系: $\rho ({\mathbf{y}}_{\mathbf{k}},{\mathbf{x}}_{\mathbf{i}})=\frac{{\sigma }_{kk}{w}_{ik}}{\sqrt{s_i}}$
• (6) 计算$k$个主成分变量对原始变量${\mathbf{x}}_i$的贡献率:$v_i={\rho }^2({\mathbf{x}}_i,({\mathbf{y}}_1,...,{\mathbf{y}}_{K^{\prime }}))={\sum }_{k=1}^{K^{\prime }}\frac{{\sigma }_{kk}^2{w}_{ik}^2}{s_i}$
• (7) 计算降维后的数据: $Y=W^{T}X$

输出:

• 降维后的数据$Y\in R^{K^{\prime }\times N}$
• 主成分变量和原变量的相关性矩阵$\mathrm{P}\in R^{K^{\prime }\times K}$
• $k$个主成分变量对原始变量的贡献率向量: $\mathbf{v}\in R^K$

注: $A=A^T\in R^{N\times N}$ $\Rightarrow$ $A$可以对角化 $\iff$ $A$有$N$个线性无关特征向量

---

五、代码

博客完整代码可以参看Github

1、基于协方差矩阵

def PCA_with_CM(X, k):
    """
    功能：使用协方差矩阵实现PCA算法
    输入:
        X: Tensor, (K,N)         # 数据矩阵
        k: int, (1)              # 目标维度
    输出:
        Y: Tensor, (k,N)         # 降维后数据
        R: Tensor, (k,K)         # Rij表示yi和xi的相关性
        v: Tensor, (K)           # vi表示所有选择的主成分对xi的贡献率
    """
    
    #(1) 归一化
    mean = torch.mean(X, 1).unsqueeze(1)  #(K, 1)
    var = torch.var(X, 1).unsqueeze(1)    #(K, 1)
    X = (X - mean) / torch.sqrt(var)
    
    
    #(2) 计算协方差矩阵
    S = X.matmul(X.t())               #(K, K)
    
    #(3) 特征分解
    lam, W = torch.eig(S, True)
    lam = lam[:, 0]
    
    #(4) 单位化， 并取
    M = torch.sqrt(torch.sum(W*W, 0)).unsqueeze(0)
    W = W / M
    top_k = torch.argsort(-lam)[:k]
    W = W[:,top_k]                           # (K, k)
    lam = lam[top_k]
    
    #(5) 计算相关系数
    s = torch.diag(S).unsqueeze(0)
    R = torch.sqrt(lam).unsqueeze(1) * W.t()/ torch.sqrt(s)   #(k, K)
    
    #(6) 计算贡献率
    v = torch.sum(R * R, 0)                #(K)
    
    #(7) 计算目标数据
    Y = W.t().matmul(X)
    
    return Y, R, v

2、基于SVD分解

def PCA_with_SVD(X, k):
    """
    功能：使用SVD分解实现PCA算法
    输入:
        X: Tensor, (K,N)         # 数据矩阵
        k: int, (1)              # 目标维度
    输出:
        Y: Tensor, (k,N)         # 降维后数据
        R: Tensor, (k,K)         # Rij表示yi和xi的相关性
        v: Tensor, (K)           # vi表示所有选择的主成分对xi的贡献率
    """
    #(1) 归一化
    mean = torch.mean(X, 1).unsqueeze(1)  #(K, 1)
    var = torch.var(X, 1).unsqueeze(1)    #(K, 1)
    X = (X - mean) / torch.sqrt(var)
    
    #(2) SVD分解
    U, Sigma, Vt = torch.svd(X, some=False)
    top_k = torch.argsort(-Sigma)[:k]
    
    #(3) 选择topk
    W = U[:,top_k]                           # (K, k)
    lam = Sigma[top_k] * Sigma[top_k]
    
    #(4) 计算相关系数
    S = X.matmul(X.t())
    s = torch.diag(S).unsqueeze(0)
    R = torch.sqrt(lam).unsqueeze(1) * W.t()/torch.sqrt(s)   #(k, K)
    
    #(5) 计算贡献率
    v = torch.sum(R * R, 0)                  #(K)
    
    #(6) 计算目标数据
    Y = W.t().matmul(X)
    
    return Y, R, v

---

3、基于sklearn库

[5] 中比较详细的描述了使用sklearn进行PCA降维的使用方式。

不过遗憾的是sklearn官方似乎没有提供相关系数矩阵$R$, 和贡献率向量$\mathbf{v}$的计算方式。这里补充他们的计算代码:

from sklearn.decomposition import PCA

pca = PCA(n_components=2, whiten=True)
pca.fit(X)

W = pca.components_                                     # 组成主成分的特征向量(K', K)
lam_sqrt = np.expand_dims(pca.singular_values_, 1)      # 对应主成分特征值的二范数，即对应主成分的奇异值，由于广播机制，因此增加维度
X -= np.mean(X, axis=0)                                 # 求协方差矩阵需要减去均值
s = np.expand_dims(np.diag(X.transpose().dot(X)), 0)    # 求方差
R = lam_sqrt * W / np.sqrt(s)                           # Rij, 第i个主成分yi与xj的相关系数
v = np.sum(R*R, 0)                                      # vj: 所选K' 主成分对xj的贡献率
X_new = pca.transform(X)                                # 降维后数据

注: sklearn输入的$X\in R^{N\times K}$

六、非线性PCA

待理解

七、概率PCA

待理解

参考资料

• [1] 刘建平Pinard: 主成分分析（PCA）原理总结
• [2] 李航统计学习方法第二版P298.
• [3] Principal Component Analysis
• [4] 刘建平Pinard：SVD分解原理
• [5] 刘建平Pinard: 用scikit-learn学习主成分分析(PCA)

附录

参考资料[1]中推导的疑惑之处:

• 疑惑一: 小节2中，使用$W{z}^i$得到恢复数据${\bar{x}}^{(i)}$
• 解答： 从过渡矩阵角度理解。新空间的一组基按列向量排列组成$W$， W是新空间到原始空间的过渡矩阵。
• 疑惑二: 小节2中式(5)->式(6)的推导过程:
• 解答 : 当$B^T$和$A$同形时，$tr(AB)=tr(BA)$
• 疑惑三: 小节3中，指出样本$x^{(i)}$在新坐标系中的投影方差为$W^T{x}^{(i)}{x}^{(i)T}W$.
• 解答: 这句话表述错了，原文中的样本应该表述为个体，方差应该是$tr(W^T{x}^{(i)}{x}^{(i)T}W)$或者$x^{(i)T}W{W}^T{x}^{(i)}$, 更严谨的，不可以说一个个体的方差，应该说样本的方差，但是最后的优化目标确实没错。
• 疑惑四：小节4中，作者将样本的协方差矩阵写为$X{X}^T$
• 解答：严格来说，应该是$\frac{X{X}^T}{n-1}$, 不过他们特征向量一样，只是特征值有成倍的关系，也不影响结果。

---

参考资料[2]中推导的疑惑之处:

• 疑惑一: P299, $\Sigma =cov(\text{x},\text{x})=E[(\text{x}-\text{u})(\text{x}-\text{u}{)}^T]$
• 解答: 对随机变量$X$, $cov(X,X)=var(X)=E[(X-u{)}^2]$. $cov(X,Y)=E[(X-u_x)(Y-u_y)]$. 书中，$\text{u}$是确定向量，$\text{x}$是$m$维随机变量($\text{x}$中的每一个维度都相当于一个$X$或$Y$)，$E$是对矩阵中的每个元素取均值。
• 疑惑二: P300，式(16.2): $E(y_i)={\alpha }_i^{T}u$
• 解答: $E(X+Y)=E(X)+E(Y)$
• 疑惑三： P300，式(16.3): $var(y_i)={\alpha }_i^T\Sigma {\alpha }_i$
• 解答:
  • $var(X+Y)=var(X)+var(Y)+2cov(X,Y)$
  • $var(aX)=a^{2}var(X)$
  • $Cov(aX,bY)=abCov(X,Y)$
  • ${\text{x}}^{T}A\text{x}={\sum }_i{\sum }_j{x}_i{x}_j{A}_{ij}$
• 疑惑四: P300，式(16.4): $cov(y_i,y_j)={\alpha }_i^T\Sigma {\alpha }_j$
• 解答:
  • $Cov(X+Y,Z)=Cov(X,Z)+Cov(Y,Z)$
  • ${\text{x}}^{T}A\text{x}={\sum }_i{\sum }_j{x}_i{x}_j{A}_{ij}$
• 疑惑五: P308页，计算$k$个主成分对第$i$个变量${\mathbf{x}}_i$的贡献率:${\rho }^2({\mathbf{x}}_i,({\mathbf{y}}_1,...,{\mathbf{y}}_k))$.
• 解答: 这是多重相关系数，计算公式可参考这里