pca算法学习笔记

写在开头

最近在研究生导师的引导下入了机器学习的坑，开始漫长啃论文作码农的生活。本篇文章主要是记录下我学习pca算法的笔记，更多是对论文的翻译（A TUTORIAL ON PRINCIPAL COMPONENT ANALYSIS，论文链接放在文章最后，很经典的一篇文章），尽量做到逻辑通顺，有错误恳请大家指正，轻喷。
PCA(Principal component analysis 主成分分析)是广泛使用的降维算法，基于线性变换的线性降维方法，即符合PX = Y（P为方阵，描述一种线性变换）。为什么需要降维，当数据样本维度很高时，无论是计算数据间距离，还是做其他计算操作，计算机开销都不小，这种情形被称为“维度灾难”（西瓜书里看的）。而在保留大部分数据信息的情况下，减少数据特征类型，缓解维度灾难。

基础

目标：找到最有‘意义’的基底来重新表示嘈杂、充满混乱的数据集。
The Goal: Principal component analysis computes the most meaningful basis to re-express a noisy, garbled data set.

框架-基变换（原理）

前文中提过，PCA是线性降维方法，满足PX = Y。设X为nxm的矩阵，此时，Y矩阵的列向量可由P矩阵的列向量线性表出，就是说P是一个空间的基底，而X（数据集）是对应的坐标值，Y是X在P空间中的映射。（$p_1$为P矩阵列向量，$y_1$为Y矩阵列向量）
$$\left[\begin{array}{ccc}{p}_1& \cdots & p_n\end{array}\right]\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1m}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nm}\end{array}\right]=\left[\begin{array}{ccc}{y}_1& \cdots & y_m\end{array}\right]\text{\hspace{1em}(1)}$$

$$\{\begin{array}{llcc}{y}_1=& p_1{x}_{11}+p_2{x}_{21}& \cdots & +p_n{x}_{n1}\\ y_2=& p_1{x}_{12}+p_2{x}_{22}& \cdots & +p_n{x}_{n2}\\ ⋮& & & \\ y_m=& p_1{x}_{1m}+p_2{x}_{2m}& \cdots & +p_n{x}_{nm}\end{array}\text{\hspace{1em}(2)}$$

从式中可以得出，P矩阵的维度（列向量的维数）与Y向量相同，当P维度数小于m（原数据集X的维度）时，新的数据矩阵Y维度也会小于m，即实现了降维。
换句话说，PCA算法的本质是高维数据在低维空间中的投影。

协方差的应用

在一个线性系统中，会造成混乱的数据一般分成两类$noise、redundancy.$（噪声和冗余），下面分别进行讨论。

Noise

顾名思义，噪声是指发生错误的数据，也就是无意义的数据，会对数据分析产生影响。
噪声的衡量不存在绝对的尺度，但所有的噪声都是相对于测量值来进行测量的。一个常用的测量方法是信噪比，也就是信息数据方差与噪声方差的比值
There exists no absolute scale for noise but rather all noise is measured relative to the measurement. A common measure is the the signalto-noise ratio (SNR), or a ratio of variances ${\sigma }^2$
$$SNR=\text{σnoise2σsignal2}\text{\hspace{1em}(3)}$$

信噪比表示出数据的精度，高信噪比的数据精度高，低信噪比的数据说明噪声数据多。

Redundancy

冗余数据是更隐晦的问题。冗余数据可以理解为两种不同特征间存在某种联系，如$y=kx$，那么y数据对应的特征就是非必需的，也就是多余的。而想知道两种测量数据间是否存在某种关系，最简单的方式是将一种数据作为横轴，另一数据作为纵轴，画出平面图观察。
下图中是一系列任意两种测量种类（也就是数据特征，维度）间的图像关系。

从图中可以看出，（c）中表示两种数据存在明显的线性关系$r_2-k{r}_1$，故在数据测量时仅需收集其中一种数据，即可以使数据集更简洁同时减小数据的采集量，起到良好的数据降维作用。

协方差

只需计算方差便可得出$SNR$，而量化数据间的关系则用到了协方差。下面以$A,B$为例，设$A,B$中数据样本已进行中心化（即每个数据减去平均值，即${\sum }_i{x}_i=0$）

A = { a 1 , a 2 , ⋯   , a n } ,   B = { b 1 , b 2 , ⋯   , b n } \begin{gathered} A=\{{a_1,a_2,\cdots,a_n}\},\ B=\{{b_1,b_2,\cdots,b_n}\} \end{gathered} A={a1​,a2​,⋯,an​}, B={b1​,b2​,⋯,bn​}​

定义方差和协方差
$$\begin{array}{c}{\sigma }_A^2=\sum a_i{a}_i,{\sigma }_B^2=\sum b_i{b}_i\\ covarianceofAandB={\sigma }_{AB}^2=\sum a_i{b}_i\end{array}$$

从式中可推出两个重要结论：

$1.{\sigma }_{AB}^2=0$ 时，$A,B$间完全不相关

$2.{\sigma }_{AB}^2={\sigma }_A^2$ 时，$A=B$

将$A、B$数据作为列向量

$$\begin{array}{c}a=[a_1{a}_2\cdots a_n]\\ b=[b_1{b}_2\cdots b_n]\end{array}$$

最后可以导出协方差公式
${\sigma }_{AB}^2$矩阵斜对角线上的数据表示方差，而非对角线元素表示不同特征间的协方差。

$$\begin{array}{c}{\sigma }_{AB}^2=\text{n−11}a{b}^T\end{array}$$

为什么不是1/n？因为1/n是有偏估计的标准化，特别是在n比较小的时候，而1/(n-1)为无偏估计的标准化。（概率论没怎么学，俺不太懂，直接放结论）
The simplest possible normalization is 1/n. However, this provides a biased estimation of variance particularly for small n. It is beyond the scope of this paper to show that the proper normalization for an unbiased estimator is 1/n−1

特征向量

在求出协方差矩阵后，让我们回到对混乱数据的讨论；协方差为0时表示数据间毫无关系，所以在我们进行基变换后的新矩阵Y，我们希望它的协方差矩阵是一个对角矩阵（除对角线元素外，其他元素都为0的矩阵），这样可以大幅度减少数据间的冗余。
而协方差矩阵是实对称矩阵（这个可以自己在稿纸上验算一下），必正交相似于对角矩阵。写到这里我们对PCA降维有了大致思想。

数学公式推导

在推导开始前，先总结一下PCA算法的思路：我们的目标是找到一个更适合的基底，来更好地重新表示我们的数据集。更好的数据集意味着数据集不混乱，于是我们定义了两种会造成混乱的数据——噪声与冗余，而为了量化两种混乱数据，我们引入方差和协方差来衡量数据集中无意义数据的量。最后我们通过对协方差矩阵的优化，得出最合适的基。

接下第一步是求出协方差矩阵。
假定数据集$x$是$m$ x $n$的数据矩阵，其中m表示测量类型的数量（也就是特征，维度），n表示数据的个数。首先算出m行数据的平均值，再用原数据矩阵减去平均值矩阵。

$$x=x-\overline{x}\text{\hspace{1em}(赋值)}$$

计算协方差矩阵$S_x=\frac{1}{(n-1)}x{x}^T$
注：如果$x$数据集中，m是特征，而n是数据个数，算平均值是每个特征的所有数据求平均，算方差也是每个特征的所有数据做方差，所以这里是$x\ast x^T$（m x m的矩阵）

第二步推导Y的协方差矩阵$S_Y$
将$Y=PX$带入$S_Y=\frac{1}{(n-1)}Y{Y}^T$

令$A=x\ast x^T$，进行对角化可得$A=E\ast D\ast E^T$；E为A矩阵特征向量组成的正交矩阵，D是A矩阵特征值作对角线元素的对角矩阵。满足$E^{-1}=E^T$

然后我们让P矩阵的每一行都是A矩阵的特征向量，即$P=E^T$。可得$A=P^T\ast D\ast P$，于是

这样，我们最后得到了一个协方差矩阵是对角矩阵的新矩阵Y，同时得到了适合的基底P——协方差矩阵的特征向量。

但新的问题出现，协方差矩阵$A$为mxm的矩阵，求出的特征向量矩阵也是mxm的，这样Y=PX得到的P矩阵仍是mxn的矩阵，与原本的维度相同。
这时，我们可以分析特征值矩阵D，由于部分特征值非常小，我们可以舍去这部分特征值，达到降维效果。
假设我们选取$k$个特征值（k<m），得出的P矩阵为kxm（P矩阵的行向量是特征向量），得到的Y矩阵变成kxn，完成降维。

小结

· $X$原数据集矩阵的主成分是协方差矩阵$X\ast X^T$的特征向量。

· 为了达到降维的效果，我们计算了量化数据混乱程度的协方差矩阵，同时使用正交对角化优化协方差矩阵，并舍去部分特征值，最后得到新的数据矩阵$Y$。

和SVD的联系

$SVD$是奇异值分解的缩写，一般我们做PCA降维时都会使用SVD。SVD简单来说是对非方阵进行的特征值分解（只有方阵可以进行特征值分解，也就是写成特征值与特征向量表示的形式）。其他概念这里不在讲述，可以自行百度（太多了，懒）。

$$\begin{array}{c}X=U\Sigma V^T\end{array}$$

而奇异值分解中需要求出协方差矩阵的特征值和特征向量（应该是右奇异值），所以在我个人认为使用SVD做PCA降维和直接使用协方差矩阵是一样的。

部分代码

论文中的代码（matlab）：

使用SVD：

我自己写的（和sklearn库对比后有些bug，仅供参考）：

import numpy as np

def pca_try_nm(self, data, n_com):
		 """
        :param data: 数据集，为了方便转成了nxm，n为数据个数，m为维度
        :param n_com: 降维后数据保留的百分比
        :return: 
        """
        y = np.array(data)
        [N, M] = y.shape
        y_aver = y.mean(axis=0)
		# 平均值处理
        a = np.zeros(shape=y.shape)
        for i in range(M):
            a[:, i] = y[:, i] - y_aver[i]
        # SVD分解
        Y = a / np.sqrt(N - 1)
        U, sigma, VT = np.linalg.svd(Y)
        # 根据输入的特征保留比例，进行降维
        eige = 0
        re_sum = 0
        for i in sigma:
            eige = eige + i
            if eige >= sigma.sum() * n_com:
                break
            else:
                re_sum = re_sum + 1
        # 得到新的数据矩阵
        signal = np.matmul(a, VT.T[:, 0:(re_sum + 1)])

本人菜鸡一枚，写得不太行，大家随便看看，希望能给大家带来点收获(｢･ω･)｢

论文：https://www.cs.princeton.edu/picasso/mats/PCA-Tutorial-Intuition_jp.pdf