主成分分析的统计学视角

主成分分析是将高维空间中的数据集拟合成一个低维子空间的方法，到目前为止它已成功应用于数学建模、数据压缩、数据可视化等领域。

主成分分析是将高维空间的数据集 { x i ∈ R D ∣ i = 1 , 2 , ⋯   , n } \{\boldsymbol{x}_i\in\mathbb{R}^D\vert i=1,2,\cdots,n\} {xi​∈RD∣i=1,2,⋯,n} 拟合到一个低维放射子空间 $S$ 中，且其维数 $d\ll D$。该问题可视为统计问题或者代数几何问题。

PCA 的统计学视角

多维随机变量 $\mathbf{x}\in {\mathbb{R}}^D$ 满足 $\mathbb{E}[\mathbf{x}]=0$，可寻找 $d(\ll D)$ 个主元 $y_i(i=1,2,\cdots ,d)$，使 $\mathbf{y}=[y_1,y_2,\cdots ,y_d{]}^{\top }$ 可表示为 $\mathbf{x}$ 的 $d$ 个不线性相关的成分
$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_d\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\mathbf{x}\\ {\mathbf{u}}_2^{\top }\mathbf{x}\\ ⋮\\ {\mathbf{u}}_d^{\top }\mathbf{x}\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\\ {\mathbf{u}}_2^{\top }\\ ⋮\\ {\mathbf{u}}_d^{\top }\end{array}\right]\mathbf{x}=U^{\top }\mathbf{x}$$
或
$$y_i={\mathbf{u}}_i^{\top }\mathbf{x},i=1,2,\cdots ,d$$
满足 ${\mathbf{u}}_i^{\top }{\mathbf{u}}_i=1,{\mathbf{u}}_i^{\top }{\mathbf{u}}_j=0$ 且 $\text{Var}[y_1]\ge \text{Var}[y_2]\ge \cdots \ge \text{Var}[y_d]$，其中，$y_1,y_2,\cdots ,y_d$ 分别称为 $\mathbf{x}$ 的第1、第2、$\cdots$、第 $d$ 个主成分.

1. 寻找第一个主成分

以第一主成分为例，我们试图寻找向量 ${\mathbf{u}}_1^{\ast }$ 使得
$$\begin{array}{r}\underset{{\mathbf{u}}_1^{\ast }\in {\mathbb{R}}^D}{\max }\text{Var}[{\mathbf{u}}_1^{\top }\mathbf{x}]\\ s.t.{\mathbf{u}}_1^{\top }{\mathbf{u}}_1=1\end{array}$$
定理：（随机变量的主成分）

对于随机变量 $\mathbf{x}\in {\mathbb{R}}^D$ 且满足 $\mathbb{E}[\mathbf{x}]=0$，协方差矩阵为 ${\Sigma }_{\mathbf{x}}=\mathbb{E}[\mathbf{x}{\mathbf{x}}^{\top }]$，假设 $\text{rank}({\Sigma }_{\mathbf{x}})\ge d$，则多维随机变量 $\mathbf{x}$ 的第 $i$ 个主成分 $y_i$ 可表示为
$$y_i={\mathbf{u}}_i^{\top }\mathbf{x}$$
其中， { u i } i = 1 d \{\boldsymbol{u}_i\}_{i=1}^d {ui​}i=1d​ 是协方差矩阵 ${\Sigma }_{\mathbf{x}}$ 的第 $i$ 个最大特征值对应的特征向量（相互正交），且 ${\mathbf{\lambda }}_i=\text{Var}[{\mathbf{y}}_i]$.

证明： 为简单起见，假定 ${\Sigma }_{\mathbf{x}}$ 无重复特征值。由 ${\Sigma }_{\mathbf{x}}{\mathbf{u}}_j={\lambda }_j{\mathbf{u}}_j$ 或 ${\mathbf{u}}_j^{\top }{\Sigma }_{\mathbf{x}}={\lambda }_j{\mathbf{u}}_j^{\top }$ 知
$$\begin{gathered} {\mathbf{u}}_i^{\top }\underbrace{{\Sigma }_{\mathbf{x}}{\mathbf{u}}_j}={\lambda }_j{\mathbf{u}}_i^{\top }{\mathbf{u}}_j \\ \underbrace{{\mathbf{u}}_i^{\top }{\Sigma }_{\mathbf{x}}}{\mathbf{u}}_j={\lambda }_i{\mathbf{u}}_i^{\top }{\mathbf{u}}_j \end{gathered}$$
即 $({\mathbf{\lambda }}_i-{\mathbf{\lambda }}_j){\mathbf{u}}_i^{\top }{\mathbf{u}}_j=0$，又由于 ${\mathbf{\lambda }}_i\ne {\mathbf{\lambda }}_j$，可知 ${\mathbf{u}}_i^{\top }{\mathbf{u}}_j=0$

由于
$$\begin{array}{rl}\mathrm{Var}\left[{\mathbf{y}}_i\right]& =\mathrm{Var}\left[{\mathbf{u}}_{\mathbf{i}}^{\top }\mathbf{x}\right]=E\left[{\left(u_i^{\top }x\right)}^2\right]\\ & =E\left[{\mathbf{u}}_{\mathbf{i}}^{\top }\mathbf{x}{\mathbf{x}}^{\top }{\mathbf{u}}_i\right]={\mathbf{u}}_{\mathbf{i}}^{\top }E\left[\mathbf{x}{\mathbf{x}}^{\top }\right]u_i={\mathbf{u}}_i{\Sigma }_x{\mathbf{u}}_i\end{array}$$
则优化问题$\max \mathrm{Var}\left[{\mathbf{y}}_1\right]$可建模为
$$\{\begin{array}{l}{\max }_{{\mathbf{u}}_1\in {\mathbb{R}}^D}{\mathbf{u}}_1^{\top }{\Sigma }_x{\mathbf{u}}_1\\ \text{ s.t. }{\mathbf{u}}_1^{\top }{\mathbf{u}}_1=1\end{array}$$

构造拉格朗日函数，将约束优化化成无约束优化
$$\mathcal{L}\left({\mathbf{u}}_1\right)={\mathbf{u}}_1^{\top }{\Sigma }_{\mathbf{x}}{\mathbf{u}}_1+\mathbf{\lambda }\left(1-{\mathbf{u}}_1^{\top }{\mathbf{u}}_1\right)$$
偏导数值为零
$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_1\right)}{\partial {\mathbf{u}}_1}=2{\Sigma }_x{\mathbf{u}}_1-2\mathbf{\lambda }{\mathbf{u}}_1=2\left({\Sigma }_x{\mathbf{u}}_1-\mathbf{\lambda }{\mathbf{u}}_1\right)=0$$
即
$${\Sigma }_x{\mathbf{u}}_1=\mathbf{\lambda }{\mathbf{u}}_1$$
可知$u_1$是协方差矩阵${\sum }_x$的特征值$\mathbf{\lambda }$对应的特征向量，最优值${\mathbf{u}}_1^{\top }{\sum }_x{\mathbf{u}}_1=\mathbf{\lambda }{\mathbf{u}}_1{\mathbf{u}}_1^{\top }={\mathbf{\lambda }}_1>0$。

2. 获取第二个主成分

第二个最优解 ${\mathbf{u}}_2$ 需要满足随机变量 $y_1={\mathbf{u}}_1^{\top }\mathbf{x}$ 与随机变量 $y_2={\mathbf{u}}_2^{\top }\mathbf{x}$ 不相关，即${\mathbf{u}}_1\perp {\mathbf{u}}_2$. 由于 $\mathbb{E}[\mathbf{x}]=0$，则 $\mathbb{E}[y_i]=\mathbb{E}[{\mathbf{u}}_i^{\top }\mathbf{x}]=0$. 两个随机变量的协方差可表示为
$$\begin{array}{l}\mathrm{Cov}\left(y_1,y_2\right)=\mathrm{Cov}\left({\mathbf{u}}_1^{\top }\mathbf{x},{\mathbf{u}}_2^{\top }\mathbf{x}\right)=E\left[\left({\mathbf{u}}_1^{\top }\mathbf{x}\right){\left({\mathbf{u}}_2^{\top }\mathbf{x}\right)}^{\top }\right]\\ =E\left[{\mathbf{u}}_1^{\top }\mathbf{x}{\mathbf{x}}^{\top }{\mathbf{u}}_2\right]={\mathbf{u}}_1^{\top }{\Sigma }_{\mathbf{x}}{\mathbf{u}}_2={\mathbf{\lambda }}_1{\mathbf{u}}_1^{\top }{\mathbf{u}}_2=0\end{array}$$
可知 ${\mathbf{u}}_1^{\top }{\mathbf{u}}_2=0$，

则优化模型为
$$\begin{array}{l}{\max }_{{\mathbf{u}}_2\in {\mathbb{R}}^D}\mathrm{Var}\left[y_2\right]={\mathbf{u}}_2^{\top }{\Sigma }_x{\mathbf{u}}_2\\ \text{ s.t. }{\mathbf{u}}_2^{\top }{\mathbf{u}}_2=1\\ {\mathbf{u}}_1^{\top }{\mathbf{u}}_2=0\end{array}$$
构造拉格朗日函数
$$\mathcal{L}\left({\mathbf{u}}_2,{\mathbf{\lambda }}_2,\mathbf{\gamma }\right)={\mathbf{u}}_2^{\top }{\Sigma }_x{\mathbf{u}}_2+{\mathbf{\lambda }}_2\left(1-{\mathbf{u}}_2^{\top }{\mathbf{u}}_2\right)+\mathbf{\gamma }{\mathbf{u}}_1^{\top }{\mathbf{u}}_2$$
置偏导数为0，得
$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_2,{\mathbf{\lambda }}_2,\mathbf{\gamma }\right)}{\partial {\mathbf{u}}_2}=2{\Sigma }_x{\mathbf{u}}_2-2{\mathbf{\lambda }}_2{\mathbf{u}}_2+\mathbf{\gamma }{\mathbf{u}}_1=0\text{\hspace{1em}(1)}$$

$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_2,{\mathbf{\lambda }}_2,\mathbf{\gamma }\right)}{\partial {\mathbf{\lambda }}_2}=1-{\mathbf{u}}_2^{\top }{\mathbf{u}}_2=0$$

$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_2,{\mathbf{\lambda }}_2,\mathbf{\gamma }\right)}{\partial \mathbf{\gamma }}={\mathbf{u}}_2^{\top }{\mathbf{u}}_2=0$$

(1) 式两边同时左乘${\mathbf{u}}_1^{\top }$得
$$\begin{array}{l}2{\mathbf{u}}_1^{\top }{\Sigma }_x{\mathbf{u}}_2-2{\mathbf{\lambda }}_2{\mathbf{u}}_1^{\top }{\mathbf{u}}_2+\mathbf{\gamma }{\mathbf{u}}_1^{\top }{\mathbf{u}}_1=0\\ 2{\mathbf{\lambda }}_1{\mathbf{u}}_1^{\top }{\mathbf{u}}_2-2{\mathbf{\lambda }}_2{\mathbf{u}}_1^{\top }{\mathbf{u}}_2+\mathbf{\gamma }=0\end{array}$$
即
$$\mathbf{\gamma }=2\left({\mathbf{\lambda }}_2-{\mathbf{\lambda }}_1\right){\mathbf{u}}_1^{\top }{\mathbf{u}}_2=0$$
则 (1) 式可简化为
$${\Sigma }_{\mathbf{x}}{\mathbf{u}}_2={\mathbf{\lambda }}_2{\mathbf{u}}_2$$
说明最优解 ${\mathbf{u}}_2$ 为协方差矩阵${\Sigma }_{\mathbf{x}}$ 的第二大特征值${\mathbf{\lambda }}_2$对应的特征向量，此时的极值
$$\max {\mathbf{u}}_2^{\top }{\Sigma }_x{\mathbf{u}}_2={\mathbf{\lambda }}_2{\mathbf{u}}_2^{\top }{\mathbf{u}}_2={\mathbf{\lambda }}_2$$
对于其余的主元$y_i$与$y_i(i\ne j)$需满足$y_i={\mathbf{u}}_i^{\top }\mathbf{x}$ 与 $y_j={\mathbf{u}}_j^{\top }\mathbf{x}$ 不相关，即
$$\mathrm{Cov}\left(y_i,y_j\right)=E\left[{\mathbf{u}}_{\mathbf{i}}^{\top }\mathbf{x}{\mathbf{x}}^{\top }{\mathbf{u}}_j\right]={\mathbf{u}}_{\mathbf{i}}^{\top }{\Sigma }_x{\mathbf{u}}_j=0$$
假设${\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_{i-1}$为协方差矩阵${\Sigma }_x$的最大$i-1$个归一化的特征向量，而最优解${\mathbf{u}}_i$定义为第$i$个主元${\mathbf{y}}_i$对应的向量（未必为特征向量）。由前过程可知
$${\Sigma }_x{\mathbf{u}}_j={\mathbf{\lambda }}_j{\mathbf{u}}_{j}j=1,2,\cdots ,i-1$$
且满足
$${\mathbf{u}}_{\mathbf{i}}^{\top }{\Sigma }_x{\mathbf{u}}_j={\mathbf{\lambda }}_j{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_j=0j=1,2,\cdots ,i-1,{\lambda }_j>0$$
即
$${\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_j=0j=1,2,\cdots ,i-1$$
最优化模型为
$$\{\begin{array}{l}\max Var[y_i]={\mathbf{u}}_i^{\top }{\Sigma }_{\mathbf{x}}{\mathbf{u}}_i\\ \text{ s.t. }{\mathbf{u}}_i^{\top }{\mathbf{u}}_i=1\\ {\mathbf{u}}_i^{\top }{\mathbf{u}}_j=0j=1,2,\cdots ,i-1\end{array}$$
构造拉格朗日函数
$$\mathcal{L}\left({\mathbf{u}}_i,{\mathbf{\lambda }}_i,{\mathbf{\gamma }}_j\right)={\mathbf{u}}_{\mathbf{i}}^{\top }{\Sigma }_x{\mathbf{u}}_i+{\mathbf{\lambda }}_i\left(1-{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_i\right)+\sum _{j=1}^{i-1}{\mathbf{\gamma }}_j{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_j$$
置偏导数为0，得
$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_i,{\mathbf{\lambda }}_i,{\mathbf{\gamma }}_j\right)}{\partial {\mathbf{u}}_i}=2{\Sigma }_x{\mathbf{u}}_i-2{\mathbf{\lambda }}_i{\mathbf{u}}_i+\sum _{j=1}^{i-1}{\mathbf{\gamma }}_j{\mathbf{u}}_j=0\text{\hspace{1em}(2)}$$

$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_i,{\lambda }_i,{\mathbf{\gamma }}_j\right)}{\partial {\lambda }_i}=1-{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_i=0$$

$$\frac{\partial \mathcal{L}\left({\mathbf{u}}_i,{\mathbf{\lambda }}_i,{\mathbf{\gamma }}_j\right)}{\partial {\mathbf{\gamma }}_j}={\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_j=0j=1,2,\cdots ,i-1$$

(2)式两边同时左乘${\mathbf{u}}_{\mathbf{j}}^{\top }$,得
$$\begin{gathered} 2{\mathbf{u}}_{\mathbf{j}}^{\top }{\Sigma }_x{\mathbf{u}}_i-2{\lambda }_i{\mathbf{u}}_{\mathbf{j}}^{\top }{\mathbf{u}}_i+\sum _{j=1}^{i-1}{\mathbf{\gamma }}_j{\mathbf{u}}_j^{\top }{\mathbf{u}}_j=0 \\ 2{\lambda }_j{\mathbf{u}}_j^{\top }{\mathbf{u}}_i-2{\lambda }_j{\mathbf{u}}_{\mathbf{j}}^{\top }{\mathbf{u}}_i+\sum _{j=1}^{i-1}{\mathbf{\gamma }}_j=0 \end{gathered}$$
即
$$\sum _{j=1}^{i-1}{\mathbf{\gamma }}_j=2\left({\lambda }_j-{\lambda }_i\right){\mathbf{u}}_{\mathbf{j}}^{\top }{\mathbf{u}}_i=0$$

由拉格朗日乘子${\mathbf{\gamma }}_j$非负，则${\mathbf{\gamma }}_j=0j=1,\cdots ,i-1$。

(72) 式可简化为
$${\Sigma }_{\mathbf{x}}{\mathbf{u}}_i={\mathbf{\lambda }}_i{\mathbf{u}}_i$$
即最优解${\mathbf{u}}_i$为协方差矩阵${\Sigma }_{\mathbf{x}}$ 第 $i$ 个特征值${\lambda }_i$ 对应的特征向量，此时的极值为
$$\text{max}{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{\Sigma }}_x{\mathbf{u}}_i={\lambda }_i{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_i={\lambda }_i=Var[y_i]$$
对于${\Sigma }_{\mathbf{x}}$ 有重复特征根的情形亦如此，略。

由上述定理可知，随机变量 $\mathbf{x}$ 的 $d$ 个主元要优于一个主元，将所有的 $d$ 个主元表示成一个向量
$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_d\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\mathbf{x}\\ {\mathbf{u}}_2^{\top }\mathbf{x}\\ ⋮\\ {\mathbf{u}}_d^{\top }\mathbf{x}\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\\ {\mathbf{u}}_2^{\top }\\ ⋮\\ {\mathbf{u}}_d^{\top }\end{array}\right]\mathbf{x}=U^{\top }\mathbf{x}$$
其中$\mathbf{y}\in R^d,U\in R^{D\times d}$,此时 $\mathbf{y}$ 的协方差矩阵可表示为
$${\Sigma }_{\mathbf{y}}=E[\mathbf{y}{\mathbf{y}}^{\top }]=E[U^{\top }\mathbf{x}{\mathbf{x}}^{\top }U]=U^{\top }{\Sigma }_{\mathbf{x}}U$$
满足$U^{\top }U=I_d$。

​ 由线性代数知识可知，对于任意可对角化的矩阵$A$，则存在由$A$的特征向量组成的列表示的矩阵$V$，有$\mathbf{\Lambda }=V^{-1}AV$，而当矩阵$A$是实对称半正定矩阵时，其特征值${\mathbf{\lambda }}_i\ge 0$，特征向量互相正交，且$V^{-1}=V^{\top }$。因此，由于${\Sigma }_x$是实对称正定矩阵，则方程${\Sigma }_y=U^{\top }{\Sigma }_{x}U$中$U$的列是协方差矩阵${\Sigma }_x$的$d$个特征向量组成。而${\Sigma }_y$是一个对角矩阵，对角元为${\Sigma }_x$的$d$个特征值。因为我们的目标是极大化${\mathbf{y}}_i$的方差$Var[{\mathbf{y}}_i]={\lambda }_i$，所以我们的结论是协方差矩阵${\Sigma }_x$的前$d$个最大特征值对应的特征向量做为$U$的列，即为目标的最优解，其极值则为${\Sigma }_y$的对角元上$d$个特征值。

3. 非零均值随机变量的主元

当 $\mathbf{x}\in R^D$ 有非零均值，则$\mathbf{x}$ 的 $d$ 个不相关主元定义为
$$y_i={\mathbf{u}}_{\mathbf{i}}^{\top }\mathbf{x}+a_{i}i=1,2,\cdots ,d$$
满足
$${\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{u}}_i=1，Var({\mathbf{y}}_1)\ge Var({\mathbf{y}}_2)\ge \cdots \ge Var({\mathbf{y}}_d)>0$$
由于随机变量 $y_i$ 满足
$$\begin{gathered} \mathbb{E}[y_i]=0 \\ \text{cov}(y_i,y_j)=0 \\ \mathbb{E}[y_i]=\mathbb{E}[{\mathbf{u}}_{\mathbf{i}}^{\top }\mathbf{x}+a_i]={\mathbf{u}}_{\mathbf{i}}^{\top }\mathbb{E}[\mathbf{x}]+{\mathbf{a}}_i={\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{\mu }}_x+a_i=0i=1,2,\cdots ,d \end{gathered}$$
因此 $a_i=-{\mathbf{u}}_{\mathbf{i}}^{\top }{\mathbf{\mu }}_x$

则
$$\begin{gathered} Var[y_1]=Var[{\mathbf{u}}_1^{\top }\mathbf{x}+a_1]=Var[{\mathbf{u}}_1^{\top }\mathbf{x}-{\mathbf{u}}_1^{\top }{\mathbf{\mu }}_x]=Var[{\mathbf{u}}_1^{\top }\left(\mathbf{x}-{\mathbf{\mu }}_x\right)] \\ =E[{\mathbf{u}}_1^{\top }(\mathbf{x}-{\mathbf{\mu }}_x)(\mathbf{x}-{\mathbf{\mu }}_x{)}^{\top }{\mathbf{u}}_1]={\mathbf{u}}_1^{\top }E[(\mathbf{x}-{\mathbf{\mu }}_x)(\mathbf{x}-{\mathbf{\mu }}_x{)}^{\top }]{\mathbf{u}}_1={\mathbf{u}}_1^{\top }{\Sigma }_{\mathbf{x}}{\mathbf{u}}_1 \end{gathered}$$
则最优解 ${\mathbf{u}}_1$ 的计算可描述为 $\text{max}Var[y_1]$

即
$$\begin{gathered} \underset{{\mathbf{u}}_1}{\max }{\mathbf{u}}_1^{\top }{\Sigma }_x{\mathbf{u}}_1 \\ {\mathbf{u}}_1^{\top }{\mathbf{u}}_1=1 \end{gathered}$$
构造拉格朗日函数
$$\mathcal{L}({\mathbf{u}}_1)={\mathbf{u}}_1^{\top }{\Sigma }_{\mathbf{x}}{\mathbf{u}}_1+{\mathbf{\lambda }}_i\left(1-{\mathbf{u}}_1^{\top }{\mathbf{u}}_1\right)$$
置拉格朗日函数偏导数为0
$$\frac{\partial \mathcal{L}({\mathbf{u}}_1)}{\partial {\mathbf{u}}_1}=2{\Sigma }_x{\mathbf{u}}_1-2{\mathbf{\lambda }}_1{\mathbf{u}}_1=0$$
得
$${\Sigma }_x{\mathbf{u}}_1={\lambda }_1{\mathbf{u}}_1$$
由此可知${\lambda }_1$和${\mathbf{u}}_1$分别为协方差矩阵${\Sigma }_{\mathbf{x}}=(\mathbf{x}-\mathbf{\mu })(\mathbf{x}-\mathbf{\mu }{)}^{\top }$ 的最大特征值与其对应的特征向量。对于地$i$个最优解${\mathbf{u}}_i$的解与前面定理的证明完全一致。

4. 零均值随机变量的样本主元

在实际应用中，我们并不知道随机变量的协方差矩阵，只能由样本点进行估计，对于独立同分布且期望为0的样本 { x i } i = 1 N \left \{ \boldsymbol x_i \right \} _{i=1}^N {xi​}i=1N​，构造样本矩阵$\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N]$，其样本协方差为
$${\Sigma }_N=\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_i{\mathbf{x}}_{\mathbf{i}}^{\top }=\frac{1}{N}\mathbf{X}{\mathbf{X}}^{\top }$$
则$d$个样本主元为
$$y_i={\hat{\mathbf{u}}}_{\mathbf{i}}^{\top }\mathbf{x}i=1,2,\cdots ,d$$
其中 { u i } i = 1 d \left \{ \boldsymbol u_i \right \} _{i=1}^d {ui​}i=1d​为${\hat{\Sigma }}_N=\frac{1}{N}\mathbf{X}{\mathbf{X}}^{\top }$或$\mathbf{X}{\mathbf{X}}^{\top }$的前$d$个特征向量。

由于$\mathbf{X}{\mathbf{X}}^{\top }\in {\mathbf{R}}^{D\times D}$是一个非常大的矩阵，所以我们可以利用$\mathbf{X}$的奇异值获得最优解，即
$$X=U_x{\Sigma }_x{V}_x^{\top }$$

$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_d\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\mathbf{x}\\ {\mathbf{u}}_2^{\top }\mathbf{x}\\ ⋮\\ {\mathbf{u}}_d^{\top }\mathbf{x}\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\\ {\mathbf{u}}_2^{\top }\\ ⋮\\ {\mathbf{u}}_d^{\top }\end{array}\right]\mathbf{x}=U^{\top }\mathbf{x}$$

5. PCA 降维案例

We will first demonstrate PCA on a 13-dimensional dataset, by loading wine dataset from sklearn, see info here.

This dataset contains chemical analysis of N=178 different wines by three different cultivators.


The analysis contains the folowing measurements:

• Alcohol

• Malic acid

• Ash

• Alcalinity of ash

• Magnesium

• Total phenols

• Flavanoids

• Nonflavanoid phenols

• Proanthocyanins

• Colour intensity

• Hue

• OD280/OD315 of diluted wines

• Proline

So overall, we have N=178 data points, lying in ${\mathbb{R}}^D$, with D=13. We stack all points together into a matrix X_wine $\in {\mathbb{R}}^{D\times N}$.

We have labels 0,1, or 2 for each wine (clutivator). The true labels are given in L_wine.

We want to see whether PCA can be helpful in the unsupervised task of clustering the 178 wines.

We start by loading the dataset, and printing the first 5 data points, just to get a general impression.

# 主成分分析算法
# 输入：
#      X: 数据矩阵大小为 n*D，每一行为 D 维向量（样本点）
# 参数：
#      dims_remain: 降维后保留的维数
#      with_std: 是否进行标准化操作，默认为进行标准化
# 返回：
#      X_reduction: 降维后的数据,数据矩阵大小为 n*d，每一行为 d 维向量（样本点）

from numpy.linalg import svd
from sklearn.preprocessing import StandardScaler

class PCA_PC:
    def __init__(self,dims_remain=2,with_std=True):
        self.dims_remain = dims_remain
        self.with_std = with_std
        
    def fit_transform(self,X):
        if self.with_std:
            ss = StandardScaler() # 此对象针对的是模式矩阵
            ss.fit(X)
            XS = ss.transform(X)
            U,_,_ = svd(XS.T) # 特征值分解函数的输入是模式矩阵的转置，输出 U 的每一列为新坐标轴
            X_reduction = XS@U[:,0:self.dims_remain]
        else:
            U,_,_ = svd(X.T)
            X_reduction = X@U[:,0:self.dims_remain]
        
        return X_reduction

调用函数

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_wine

if __name__ == '__main__':     
    X_wine, L_wine = load_wine(return_X_y=True)
    np.set_printoptions(suppress=True)
    
    model1 = PCA_PC(dims_remain=2,with_std=False)
    X_reduct1 = model1.fit_transform(X_wine)
    plt.figure(figsize=(15,6))
    plt.subplot(121),plt.scatter(X_reduct1[:,0], X_reduct1[:,1], c=L_wine)
    plt.title('Unstandard Preprocessing')
    
    model2 = PCA_PC(dims_remain=2,with_std=True)
    X_reduct2 = model2.fit_transform(X_wine)
    plt.figure(figsize=(15,6))
    plt.subplot(121),plt.scatter(X_reduct2[:,0], X_reduct2[:,1], c=L_wine)
    plt.title('Standard Preprocessing')
    
    plt.show()