Principal Component Analysis.

主成分分析(Principal Component Analysis,PCA)是一种线性数据降维方法。该算法在降维时去除数据的相关性，且最大限度保持原始数据的方差信息。

主成分分析的步骤是对归一化数据的协方差矩阵进行特征值分解(或直接对归一化的数据矩阵进行奇异值分解)。

本文目录：

1. 几何解释
2. 线性变换的解释
3. 最大投影方差角度
4. 最小重构代价角度
5. 奇异值分解角度
6. 主坐标分析
7. 概率主成分分析
8. PCA的实现

1. 几何解释

如上图所示的二维数据，如果将这些数据投影到一维空间中，选择数据方差最大的方向进行投影，才能最大化数据的差异性，保留更多的原始数据信息。

从概率的角度理解，投影后数据方差越大，表示数据分布越分散，则数据分布概率越分散，每个数据点的概率取值就比较小；此时该分布的信息熵更大，包含更多信息。

由上图可以看出，这$N$个样本点沿着$f_1$轴方向有最大的离散性，$f_1$的方向称为第一个主成分。

为了去掉相关性，第二个主成分应该正交于第一个主成分。

主成分分析试图在力保数据信息丢失最少的原则下，去除数据的相关性，对高维空间的数据降维处理。

2. 线性变换的解释

若记$N$个样本点\(\{x_1,x_2,...,x_N\},x_n \in \Bbb{R}^d\)，降维后得到\(\{z_1,z_2,...,z_N\},z_n \in \Bbb{R}^k\)，通常降维后的数据维度低于原维度，即$k<d$；

值得一提的是，PCA会对原始数据进行归一化，即对每个样本减去所有样本的平均值:

线性降维可以用一个线性变换矩阵$W \in \Bbb{R}^{k×d}$来表示:

对降维后的样本进行重构，为减小参数量，选择重构矩阵为$W^T$:

目标函数为最小化重构误差：

注意到矩阵$W^TW$是对称矩阵，可以进行相似对角化：

其中：

• $U^T$是正交矩阵，表示对一个向量进行旋转，不改变其长度；
• $Σ$是对角矩阵，由于$W$的秩不超过$k$,$Σ$的秩也不超过$k$,最多有$k$个对角元素不为零；表示对旋转后向量的最多$k$个维度进行长度的缩放，其余维度置零；
• $U$是正交矩阵，表示对缩放后的向量旋转回原来的方向。

则优化的目标函数进一步表示为：

先求最优的$Σ$：最优的$Σ$使$I-Σ$接近零矩阵，由于$Σ$最多有$k$个对角元素不为零，则最优的$Σ$可以表示为\(Σ= \begin{bmatrix} I_k & 0 \\ 0 & 0 \end{bmatrix}\)，

则优化的目标函数进一步表示为：

由于$U$是正交矩阵，可以把最小化问题转换成最大化问题：

先考虑$k=1$的情况，此时优化问题为：

上述问题等价于求瑞利商的极值。使用拉格朗日乘数法解上述约束优化问题，定义拉格朗日函数：

上式对$u$求偏导并置$0$，得：

解得：

即$u$和$λ$分别是矩阵$XX^T \in \Bbb{R}^{d×d}$的特征向量和特征值，此时优化问题化简为：

$λ$是矩阵$XX^T$最大的特征值。

根据数学归纳法，该结论可推广至$k>1$的情况。当$k>1$时，$λ$是矩阵$XX^T$最大的前$k$个特征值，$u$是与之对应的特征向量。

从重构角度可以选择合适的降维维度$k$。人为设定重构阈值$t$，然后选取使得下式成立的最小$k$值：

3. 最大投影方差角度

若记$N$个样本点\(\{x_1,x_2,...,x_N\},x_n \in \Bbb{R}^d\)，降维后得到\(\{z_1,z_2,...,z_N\},z_n \in \Bbb{R}^k\)，通常降维后的数据维度低于原维度，即$k<d$；

记原数据空间中的一个单位向量$u$，对原始数据沿该方向进行投影，希望投影差异最大。则每一个样本点$x_n$在该方向上的投影为$u^Tx_n$，这样计算的是“绝对”差异，即所有从原点出发的数据向量投影计算的差异。我们并不关心数据本身的大小，而是希望最大化数据与其数据平均中心的差异，需要对原始数据进行预处理：

将预处理后的数据投影到$u$方向上，便可以得到在该方向上的投影差异(注意投影结果是标量)：

对投影差异的解释如下。若向量$x_n - \overline{x}$与向量$u$之间的夹角为$\theta$，则将向量$x_n - \overline{x}$直接投影到向量$u$上的长度为$|x_n - \overline{x}|\cdot cos \theta$。直接计算两向量内积$(x_n - \overline{x}) \cdot u=|x_n - \overline{x}|\cdot |u|\cdot cos \theta$；当$|u|=1$时两者等价，因此用内积表示投影差异。

定义投影方差为所有数据投影差异的平方平均值：

这是一个约束优化问题：

上述问题等价于求瑞利商的极值。使用拉格朗日乘数法解上述约束优化问题，定义拉格朗日函数：

上式对$u$求偏导并置$0$，得：

解得：

记样本的协方差矩阵$S=\frac{1}{N}\sum_{n=1}^{N} {(x_n - \overline{x})(x_n - \overline{x})^T} \in \Bbb{R}^{d×d}$，则上式表示为：

即$u$和$λ$分别是矩阵$S$的特征向量和特征值，此时优化问题化简为：

$λ$是矩阵$S$最大的特征值。

值得一提的是，有时也用样本的散布矩阵(scatter matrix) $S=\sum_{n=1}^{N} {(x_n - \overline{x})(x_n - \overline{x})^T} \in \Bbb{R}^{d×d}$进行计算。由于散布矩阵和协方差矩阵仅相差一个系数$\frac{1}{N}$，不影响结果，因此以下不做区分。

4. 最小重构代价角度

若记$N$个样本点\(\{x_1,x_2,...,x_N\},x_n \in \Bbb{R}^d\)，降维后得到\(\{z_1,z_2,...,z_N\},z_n \in \Bbb{R}^k\)，通常降维后的数据维度低于原维度，即$k<d$；

在原向量空间取一组完备的单位正交基\(\{u_1,...,u_d\}\)，其中前$k$维代表降维后的方向。

原向量空间中的数据$x_n$可以表示为（任意向量由基向量线性表示）:

注意到此处$u_i^Tx_n$表示将向量$x_n$投影到单位向量$u_i$上的投影长度。

而降维后的数据$\hat{x}_n$表示为：

定义所有样本的重构代价为：

其中$S=\frac{1}{N}\sum_{n=1}^{N} {(x_n - \overline{x})(x_n - \overline{x})^T} \in \Bbb{R}^{d×d}$是协方差矩阵。

这是一个约束优化问题：

上述问题等价于求瑞利商的极值。使用拉格朗日乘数法解上述约束优化问题，可以得到$u_{k+1}$到$u_{d}$是$S$最小的$d-k$个特征值，即$u_{1}$到$u_{k}$是$S$最大的$k$个特征值。

5. 奇异值分解角度

由上述分析可知，主成分分析可以通过对样本的散布矩阵$S$(或协方差矩阵)进行特征值分解，其最大的前$k$个特征值对应的特征向量就是所求的$k$个主方向。

回顾样本的协方差矩阵$S$计算为：

记样本矩阵$X\in \Bbb{R}^{n \times d}$，单位矩阵$I_n \in \Bbb{R}^{n \times n}$和全$1$矩阵$\Bbb{1}_n \in \Bbb{R}^{n \times 1}$，进行化简：

记中心矩阵(centering matrix,该矩阵没找到合适的中文描述,主要用于为散布矩阵提供一个简洁的表达方法)$H=I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T$，$H$矩阵是对称的幂等矩阵，其性质如下：

• $H^T=H$
  1. 证明:$H^T=(I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T)^T = I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T =H$
• $H^2=H$
  1. 证明:

注意到矩阵$H$相当于对样本矩阵$X$进行预处理，即$\overline{X}=HX$。通过引入矩阵$H$，散布矩阵$S$可以被表示为：

对散布矩阵$S$进行特征值分解：

取$Q$的前$k$列$Q_k \in \Bbb{R}^{d \times k}$，将样本矩阵的特征维度从$d$降至$k$：

在上述流程中，对散布矩阵$S$进行特征值分解有时会有较大的计算量(当特征维度很高时)。若直接对预处理后的样本矩阵$HX$进行奇异值分解：$HX=U \Sigma V^T$。则散布矩阵$S$也可以表示为：

注意到$V=Q$,$\Sigma^2=\Lambda$。因此实践中不直接对散布矩阵$S$进行特征值分解，而是对预处理的样本矩阵$HX$进行奇异值分解，也可以实现线性降维。降维结果可以计算为：

6. 主坐标分析

对于样本矩阵$X \in \Bbb{R}^{n \times d}$，有时样本数量小于特征维度($n<d$)，此时对散布矩阵$S=X^TH^THX \in \Bbb{R}^{d \times d}$进行特征值分解会有较大的计算量。

考虑另一矩阵$T = HXX^TH \in \Bbb{R}^{n \times n}$。引入样本矩阵$HX$的奇异值分解$HX=U \Sigma V^T$，则$T$可以被表示为：

根据第$5$节的结论，降维结果可以表示为$X_{\text{dim-reduced}} = U \Sigma_k$。因此直接对矩阵$T$进行特征值分解也能得到降维结果。

上述方法被称作主坐标分析(Principal Coordinate Analysis,PCoA)。

在PCA中计算的协方差矩阵(散布矩阵) $S$衡量样本的不同特征维度之间的相关性，维度为$d \times d$；而PCoA中计算的矩阵$T$衡量不同样本之间的相关性，维度为$n \times n$。

7. 概率主成分分析

概率主成分分析(probabilistic PCA)把原始样本特征$x$看作观测变量(observed variable)，把降维后的特征$z$看作隐变量(latent variable)。

把特征$z$的先验分布建模成$k$维标准正态分布$N(0,I_k)$，使用线性高斯模型(linear Gaussian model)建模特征$z$与$d$维特征$x$的关系：

其中引入噪声$\epsilon$服从分布$N(0,σ^2I_d)$。

pPCA主要解决两个问题：

• 学习 learning：估计模型参数$w$,$\mu$,$σ^2$，通常用EM算法求解。
• 推断 inference：由样本特征$x$推断特征$z$的后验分布$P(z | x)$。

下面介绍求解后验分布$P(z | x)$的过程。

则数据特征$x$的分布可表示为$P(x) \text{~} N(\mu,w^Tw+σ^2I_d)$；

构造联合分布$P(\begin{bmatrix} x \ z \ \end{bmatrix})$:

其中协方差矩阵$Σ_{xz}$计算为：

根据联合分布$P(\begin{bmatrix} x \ z \ \end{bmatrix})$和特征$x$的分布$P(x)$可以求条件分布$P(z | x)$。

8. PCA的实现

① PCA from scratch

由上述介绍，PCA的过程是对归一化样本的散布矩阵$S$(或协方差矩阵)进行特征值分解，其最大的前$k$个特征值对应的特征向量就是所求的$k$个主方向。

def PCA(data, k):
    # 数据归一化
    data = data - np.mean(data, axis=0, keepdims=True) # (n, d)
    # 计算散布矩阵
    S = np.dot(data.T, data) # (d, d)
    # 计算特征值和特征向量
    eig_values, eig_vectors = np.linalg.eig(S)
    # 选择前k个最大的特征值标号
    index = np.argsort(-eig_values)[:k]
    # 选择对应的特征向量(主方向)
    PCA_vectors = eig_vectors[index, :] # (k, d)
    # 降维
    reduced_data = np.dot(data, PCA_vectors.T) # (n, k)
    return reduced_data

② PCA from sklearn

PCA也可以通过sklearn库快速实现：

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

digits = load_digits()
X_pca = PCA(n_components=2).fit_transform(digits.data)

plt.figure(figsize=(10, 5))
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=digits.target,label="PCA")
plt.legend()
plt.show()