PCA主成分降维原理简单推导 + Python代码实现

前言

​ PCA主成分降维是降维的常用手段。其原理简单易理解。本文将从原理简单地进行推导求解，过程并不严谨。主要是在其代码实现。

原理推导

​ 首先，何为降维？

​ 在日常生活中，我们的数据往往维度非常的高，导致我们的运算过于复杂，也会导致过拟合，因此，为了解决这个问题。我们除了直接去除一些不必要的维度以外，还可以进行数据的降维（不直接舍弃某些维度，而是采用数学求解）。

​ 在PCA中，我们常说的降维实际上是投影。啥是投影？我们来看看下面这张图

​ 在坐标系中，存在向量u,v。

投影

​ 有一说法，向量u在向量v的投影。如下图的向量AD就是投影（相当于从u点用手电筒向着v的方向照射得出来的影子）

​ 那么投影的公式又是什么呢？在此处我们不作表达，因为我们真正用到的只是投影的长度AD。而投影长度的公式很直观，当燃是$|u|cos(a)$，其中$a$是向量u，v的夹角。$|u|$是u的模长。

公式转换

​ 定义向量u，v
$$u=\left(\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_n\end{array}\right);v=\left(\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_n\end{array}\right)$$

​ 上过高中的都知道，我们有计算向量内积的公式
$$u\cdot v=|u||v|cos(a)$$
​ 假如我们令$|v|$恒等于1，那么公式就可以变成$u\cdot v=|u|cos(a)$，发现了吗，后面这一段正是投影长度的公式。

​ 又因为在数据中，我们只有数据点，没有夹角，所以按照我们另一个计算向量内积的公式。
$$u\cdot v=u_1{v}_1+u_2{v}_2+\cdots u_n{v}_n$$
​ 其中，$u_1{v}_1$分别代表u，v的第一个维度分量。以此类推

​ 所以
$$u_1{v}_1+u_2{v}_2+\cdots u_n{v}_n=|u|cos(a)$$
​ 转化成矩阵相乘。就变成了
$$u^{T}v=|u|cos(a)$$

数学推导

函数损失

​ 对于坐标系中的点A-E，我们希望其投影到向量u之后方差尽可能的大。那么很显然可以设定损失函数即为方差。

​ 方差的计算我们都很熟悉，但是在此之前我们不妨做个定义，假设所有的点都已经中心化了，也就是对于原始坐标轴中的每个点，他们在各个坐标方向上的点都是均值为0，实际上就是每一个点都减去该点的均值。最后再进行投影。那么投影长度就是
$$(x-\bar{x}{)}^{T}u$$
​ 因此原始向量的均值为0，所以，投影后的向量均值也是0，那么函数的损失就可以写成
$$\begin{array}{c}\begin{array}{rl}f(u)& =\frac{1}{n}\sum _{i=1}^n{\left[(x_i-\bar{x}{)}^{T}u\right]}^2\\ & s.t.u^{T}u=1\end{array}\end{array}$$
​ 其中，$u^{T}u=1$代表$|u|=1$，只是转化成了矩阵相乘的方式罢了。

​ 对函数求平方进行展开，又因为式子中的$(x_i-\bar{x}{)}^{T}u$是一个常数，所以我们对其进行转置是不会改变最终的结果
$$\begin{array}{c}\begin{array}{rl}f(u)& =\frac{1}{n}\sum _{i=1}^n\left[(x_i-\bar{x}{)}^{T}u\ast (x-\bar{x}{)}^{T}u\right]\\ & =\frac{1}{n}\sum _{i=1}^n\left[u^T(x_i-\bar{x})(x_i-x{)}^{T}u\right]\\ & =u^T\left(\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x})(x_i-\bar{x}{)}^T\right)u\\ & s.t.u^{T}u=1\end{array}\end{array}$$
​ 可以看出，括号之中的不就是协方差矩阵的公式吗（不熟悉的请自行网上搜索查阅，此处不做解释）？为了表达简便，我们用$S=\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x})(x_i-\bar{x})$。

​ 因此，我们的目标就是
$$\begin{array}{c}\begin{array}{rl}& \underset{u}{\max }f(u)\\ =& \underset{u}{\max }\left(u^{T}Su\right)\\ & s.t.u^{T}u=1\end{array}\end{array}$$

求解（拉格朗日乘数法）

​ 面对带约束的极值问题，最让人容易想到的，当燃是使用拉格朗日乘数法（不熟悉的请自行网上搜索查阅，此处不做解释）求解。

​ 构造拉格朗日函数。
$$L(u)=u^{T}Su+\lambda (1-u^{T}u)$$
​ 对$u$求导
$$\frac{\partial L(u)}{\partial u}=2Su-2\lambda u=0$$
​ 移项得
$$Su=\lambda u$$
​ 啥意思呢？很容易看出u其实是S的一个特征向量，而λ是对应的特征值。将式子回代入原函数。
$$f(u)=u^{T}Su=u^T\lambda u=\lambda u^{T}u=\lambda$$
​ 因此，可以得出，投影长度方差最终就等于特征值，要最大化投影长度方差，就相当于要取协方差矩阵的最大特征值，而求最大次投影长度方差，也对应其次最大特征值，依此类推。

实现流程

​ ①对数据进行中心化。

​ ②求出其对应的协方差矩阵。

​ ③选定k个维度。

​ ⑤计算协方差矩阵的特征值。

​ ⑥对特征值进行从大到小排序，并保留前k个，并计算对应的特征向量。

​ ⑦把得到的特征向量按列拼接，用原始数据相乘，实现降维。

代码实现

import numpy as np
from scipy import stats
class PCA():
    def __init__(self):
        pass
    def train(self,data,k):
        '''
        :param data: 数据
        :param k: 保留前k维度
        :return: 返回降维后的维度
        '''
        data_mean=data-np.mean(data,axis=0) #中心化
        S=1/(data.shape[0]-1)*data_mean.T@data_mean #协方差矩阵
        eigenvalues,eigenvectors=np.linalg.eig(S) #获取特征值和特征向量
        indexs=np.argsort(-eigenvalues)[:k] #取前k个最大的特征值对应的索引
        vectors=eigenvectors[:,indexs]#取出前k个最大特征值对应的特征向量
        result=data_mean@vectors #数据乘以特征向量完成降维
        return result
if __name__ == '__main__':
    x=stats.norm.rvs(0,10,(100,2)) #生成数据
    pca=PCA() #初始化
    result=pca.train(x,1) #降维

结束

​ 过程并不严谨，有一些地方没有给予说明，由于过程对于我来说过于复杂，因此本文并没有提出，但大体方向应当是对的。有何错误，遗漏的地方，烦请指正。阿里嘎多