12-降维-PCA-最大投影方差/最小重构代价

1.PCA(principal component analysis)主成分分析定义

PCA主要思想是将N维特征映射到K维上，这里的K维是全新的正交特征，是相互垂直的线性无关向量，也称为主成分，K维是在原有的N维特征的基础上重新构造出来的K维特征。PCA的工作是从原始的空间中顺序的找到一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。

2.PCA中心思想：一个中心，两个基本点

• 一个中心：原始特征空间的重构；
  将特征数据由线性相关转换成线性无关$(线性相关\mapsto 线性无关)$
• 两个基本点：最大投影方差和最小重构距离
  两种对原始特征空间的重构特征方式，我们通过要求转换后的数据拥有最大投影方差和最小重构距离来实现最好的线性相关到线性无关的转换。
  
  $注：当数据向u_1方向投影的时候，可以看出数据的投影距离比较扩散，相当于方差较大，所以u_1方向好\surd$
  $注：当数据向u_2方向投影的时候，可以看出数据的投影距离比较集中，相当于方差较小，所以u_2方向差\times$
  因为数据在$u_1$方向的重构距离更小，投影方差最大，所以满足最小重构距离和最大投影方差这两个特性

3.PCA数学建模-最大投影方差

3.1数据中心化处理-最大投影方差

将数据移动到原点，使得期望为零，方便计算
$$X_i-\bar{X};中心化处理\text{\hspace{1em}(1)}$$
$定义L_i为向量X_i-\bar{X}在u_1方向的投影,为方便计算，我们定义||u_1||=||u_1^T{u}_1||=1$
$$L_i=(X_i-\bar{X}{)}^T{u}_1;st||u_1^T{u}_1||=1\text{\hspace{1em}(2)}$$
我们的目标是满足：最大投影方差,所以我们定义目标函数J
$由于L_i为投影，并且数据已经中心化了，所以L_i^2为投影方差$
$$J=argmax\sum _{i=1}^N\frac{1}{N}{L}_i^2=argmax\frac{1}{N}\sum _{i=1}^N((X_i-\bar{X}{)}^T{u}_1{)}^2;st:||u_1^T{u}_1||=1\text{\hspace{1em}(3)}$$
$$J=argmax\frac{1}{N}\sum _{i=1}^N{u}_1^T\cdot (X_i-\bar{X})\cdot (X_i-\bar{X}{)}^T\cdot u_{1}st:||u_1^T{u}_1||=1\text{\hspace{1em}(4)}$$
$$J=argmax{u}_1^T[\frac{1}{N}\sum _{i=1}^N(X_i-\bar{X})\cdot (X_i-\bar{X}{)}^T]u_{1}st:||u_1^T{u}_1||=1\text{\hspace{1em}(4)}$$
$$J=argmax(u_1^T\cdot S\cdot u_1)st:||u_1^T{u}_1||=1\text{\hspace{1em}(5)}$$
我们的目标是为了找到方向$u_1$,所以我们用拉格朗日乘子法求解带约束问题：
$$\hat{u}=argmax(u_1^T\cdot S\cdot u_1);st||u_1^T{u}_1||=1\text{\hspace{1em}(6)}$$
转换成：
$$L(u_1,\lambda )=u_1^T\cdot S\cdot u_1+\lambda (1-u_1^T{u}_1)\text{\hspace{1em}(7)}$$
对上式$L(u_1,\lambda )求导：$
$$\frac{\partial L(u_1,\lambda )}{\partial u_1}=2S{u}_1-2\lambda u_1=0\text{\hspace{1em}(8)}$$

3.2 结论-最大投影方差

$$结论：S\widehat{u_1}=\lambda \widehat{u_1}\text{\hspace{1em}(9)}$$
$注：\widehat{u_1}=方差矩阵S所对应的特征向量；\lambda =特征向量对应的特征值$
$我们要求的最大投影方差所在的向量为样本方差S所在的特征向量\widehat{u_1}$

4.PCA数学建模-最小重构距离

4.1 核心思想

最小重构距离的含义是：首先我们将原始空间的特征向量进行重构，转换成由一组线性无关的向量组成，最小重构就是我们要求转换后的数据$\hat{X}与原来的数据X之间的损失最小,即我们用新的坐标轴去恢复数据的时候损失最小。$
$$J=\frac{1}{N}\sum _{i=1}^N||X_i-\widehat{X_i}|{|}^2\text{\hspace{1em}(10)}$$

4.2 坐标转换

(我们定义数据是实现样本X中心化后:$X_i-\bar{X},u_1为此方向的单位向量；大小：L_i=(X_i-\bar{X}{)}^T{u}_1；方向：u_1$
假设样本数据$X_i$在原来的坐标系$x_1,x_2上的坐标为(x_1,x_2),我们将其转换至新的坐标系u_1,u_2上，可得：$
$$X_{u1}=((X_i-\bar{X}{)}^T{u}_1)u_1\text{\hspace{1em}(11)}$$
同理对于第k个方向，我们可以得到$u_k$方向的新坐标：
$$X_{uk}=((X_i-\bar{X}{)}^T{u}_k)u_k\text{\hspace{1em}(12)}$$
那么我们定义原来的数据X在新的坐标系下的表示(p维)：
$$X_i=\sum _{k=1}^p((X_i-\bar{X}{)}^T{u}_k)u_k\text{\hspace{1em}(13)}$$
那么我们定义降维后的数据$\hat{X}$在新的坐标系下的表示(q维)：
$$\widehat{X_i}=\sum _{k=1}^q((X_i-\bar{X}{)}^T{u}_k)u_k\text{\hspace{1em}(14)}$$
那么我们可以将J表示如下：
$$J=\frac{1}{N}\sum _{i=1}^N||X_i-\widehat{X_i}|{|}^2\text{\hspace{1em}(15)}$$
$$J=\frac{1}{N}\sum _{i=1}^N||\sum _{k=1}^p((X_i-\bar{X}{)}^T{u}_k)u_k-\sum _{k=1}^q((X_i-\bar{X}{)}^T{u}_k)u_k|{|}^2\text{\hspace{1em}(16)}$$
$$J=\frac{1}{N}\sum _{i=1}^N||\sum _{k=q+1}^p((X_i-\bar{X}{)}^T{u}_k)u_k|{|}^2\text{\hspace{1em}(17)}$$
$$J=\frac{1}{N}\sum _{i=1}^N\sum _{k=q+1}^p[(X_i-\bar{X}{)}^T{u}_k{]}^2\text{\hspace{1em}(18)}$$
$$J=\sum _{k=q+1}^p\sum _{i=1}^N\frac{1}{N}[(X_i-\bar{X}{)}^T{u}_k{]}^2\text{\hspace{1em}(19)}$$
$由于{\sum }_{i=1}^N\frac{1}{N}[(X_i-\bar{X}{)}^T{u}_k{]}^2=u_k^{T}S{u}_k;且，u_k^T{u}_k=1$
$$J=\sum _{k=q+1}^p{u}_k^{T}S{u}_k;st:u_k^T{u}_k=1;\text{\hspace{1em}(20)}$$
同样我们用拉格朗日乘子法：
$$J=argmi{m}_{(u_{q+1},\cdots ,u_p,\lambda )}\sum _{k=q+1}^p{u}_k^{T}S{u}_k+\lambda (1-u_k^T{u}_k)\text{\hspace{1em}(21)}$$
损失函数最小取在S特征向量对应的特征值中剩下最小的几个值之和。

5. PCA 和最小二乘法的关系

5.1 PCA主成分分析

• PCA 步骤如下：
  – step1: 将样本 X标准化
  $$X_{\sigma }=\frac{X-\mu }{\sigma }\text{\hspace{1em}(22)}$$
  – step2: 求协方差矩阵$X_{\sigma }{X}_{\sigma }^T$：本质上是构建投影向量$u_0$,投影方程如下：
  $$L=u_0^T{X}_{\sigma }{X}_{\sigma }^T{u}_0,st:u_0^T{u}_0=1\text{\hspace{1em}(23)}$$
  – step3:应用拉格朗日乘子法可得：
  $$L=u_0^T{X}_{\sigma }{X}_{\sigma }^T{u}_0+\lambda (1-u_0^T{u}_0)\text{\hspace{1em}(24)}$$
  – step4: 求偏导：
  $$\frac{\partial L}{\partial u_0}=2X_{\sigma }{X}_{\sigma }^T{u}_0-2\lambda u_0=0\to X_{\sigma }{X}_{\sigma }^T{u}_0=\lambda u_0\text{\hspace{1em}(25)}$$
  – step5: 也就是说协方差矩阵$X_{\sigma }{X}_{\sigma }^T$的特征值和特征向量就是使得投影尺寸最大或者最小，如果我们选择特征值大的k部分和对应的特征向量，那么我们就能得到最大投影尺寸，也就是最小重构尺寸。
  – step6: 对特征值从大到小排序，选择其中最大的k个。然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵,将协方差矩阵$X_{\sigma }{X}_{\sigma }^T$的最大值k和对应的特征向量构成新的矩阵$X_k$
  – step7: 将数据转换到k个特征向量构建的新空间中

5.2 最小二乘法

最小二乘法步骤：
– step1: 构建损失函数：
$$L=(Ax-b{)}^T(Ax-b)=(x^T{A}^T-b^T)(Ax-b)=x^T{A}^{T}Ax-2x^T{A}^{T}b+b^{T}b\text{\hspace{1em}(26)}$$

• step2: 求偏导：
  $$\frac{\partial L}{\partial u_0}=2A^{T}Ax-2A^{T}b=0\to A^{T}Ax=A^{T}b\text{\hspace{1em}(25)}$$
• step3: 如果$A^{T}A$可逆，那么可得：
  $$\frac{\partial L}{\partial u_0}=2A^{T}Ax-2A^{T}b=0\to A^{T}Ax=A^{T}b\text{\hspace{1em}(25)}$$
• step4: 可得解如下：
  $$\hat{x}=(A^{T}A{)}^{-1}{A}^{T}b\text{\hspace{1em}(25)}$$

5.3 关系如图：

5.3 思考

• 思考：$Av=\sigma u$ ，那么可以根据四个子空间得到如下图：
  
  假设我们有一群人，
• 最小二乘法求解： 因为v属于行空间，可以理解为主要特征，
  $$A=U\Sigma V^T\to A^{T}A=V{\Sigma }^T\Sigma V^T$$

1：我们可以按照身高，年龄，性别等主要特征区分，就是按行空间V区分,。每一个v无法组成单独的个体，是每个人的特征===>>> 最小二乘法求解,最小二乘法的目的是为了找到一条直线，这个直线可以近似拟合所有的点，所以最小二乘法是为了找到样本中的点的特征

• PCA 主成分分析： 因为u 属于列空间，就是找几个典型样本代表整体
  $$A=U\Sigma V^T\to A{A}^T=U\Sigma {\Sigma }^T{U}^T$$

2：我们可以按每一个样本人来区分就是按照列空间U区分，每一个u都是单独的个体，比如可以按照中国人，日本人区分 ===>>> PCA主成分分析。PCA的目的是为了用小部分的样本来近似替代总体的样本，所以pca主成分分析是为了找到主要的群体，这部分群体占比比较大，可以将这部分人近似的代表总数人。

6. python 代码[有点问题，后续更新]

• 代码：

import numpy as np

np.set_printoptions(suppress=True, precision=3)
from sklearn.decomposition import PCA


class PcaTest(object):
    def __init__(self, matrix, k):
        self.matrix = matrix
        self.k = k
        self.standard_matrix = np.zeros_like(self.matrix)
        # self.pca_result =
        self.xxt_matrix = np.zeros((self.matrix.shape[0], self.matrix.shape[0]))
        self.pca_step()

    def pca_step(self):
        self.get_standard_matrix()
        self.get_xxt_matrix()
        self.get_xk_matrix()

    def get_standard_matrix(self):
        matrix_mu = np.mean(self.matrix, axis=0)
        # matrix_sigma = np.std(self.matrix, axis=0)
        matrix_sigma = np.std(self.matrix, axis=0, ddof=1)
        self.standard_matrix = (self.matrix - matrix_mu) / matrix_sigma
        print(f"matrix=\n{self.matrix}")
        print(f"standard_matrix=\n{self.standard_matrix}")
        return self.standard_matrix

    def get_xxt_matrix(self):
        # self.xxt_matrix = self.standard_matrix @ self.standard_matrix.T /(self.matrix.shape[0]-1)
        self.xxt_matrix = np.cov(self.standard_matrix, rowvar=False)
        print(f"xxt_matrix=\n{self.xxt_matrix}")
        return self.xxt_matrix

    def get_xk_matrix(self):
        x_k = min(self.k, self.matrix.shape[1])
        pca_xxt_matrix = self.get_xxt_matrix()
        xxt_evalue, xxt_evector = np.linalg.eig(pca_xxt_matrix)
        sorted_index = np.argsort(xxt_evalue)
        sorted_evector = xxt_evector[:, sorted_index[:-self.k - 1:-1]]
        print(f"xxt_evalue={xxt_evalue}")
        print(f"xxt_evector={xxt_evector}")
        pca_p = sorted_evector
        pca_x = self.standard_matrix
        pca_y = pca_x@pca_p
        print(f"pca_y=\n{pca_y}")

        return xxt_evector


if __name__ == "__main__":
    my_matrix = np.arange(20).reshape(4, 5)
    k = 3
    pca_matrix = PcaTest(my_matrix, k)
    sklearn_pca = PCA(n_components=k)
    sklearn_pca.fit(my_matrix)
    sklearn_y = sklearn_pca.transform(my_matrix)
    print(f"sklearn_y=\n{sklearn_y}")

• 结果：

matrix=
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]]
standard_matrix=
[[-1.162 -1.162 -1.162 -1.162 -1.162]
 [-0.387 -0.387 -0.387 -0.387 -0.387]
 [ 0.387  0.387  0.387  0.387  0.387]
 [ 1.162  1.162  1.162  1.162  1.162]]
xxt_matrix=
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
xxt_matrix=
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
xxt_evalue=[ 5.  0.  0.  0. -0.]
xxt_evector=[[-0.447 -0.444 -0.    -0.     0.   ]
 [-0.447 -0.641 -0.866  0.    -0.215]
 [-0.447  0.362  0.289 -0.    -0.719]
 [-0.447  0.362  0.289 -0.707  0.467]
 [-0.447  0.362  0.289  0.707  0.467]]
pca_y=
[[ 2.598  0.     0.   ]
 [ 0.866  0.     0.   ]
 [-0.866 -0.    -0.   ]
 [-2.598 -0.    -0.   ]]
sklearn_y=
[[ 16.771   0.      0.   ]
 [  5.59    0.      0.   ]
 [ -5.59   -0.     -0.   ]
 [-16.771  -0.     -0.   ]]