【主成分分析】PCA降维算法及Matlab代码实现
前言
机器学习中经常会遇到高维变量的大数据集，并且大数据集的很多高维变量之间并不是独立的，它们之间往往存在相关关系。这些变量一方面为机器学习提供了大量的信息，另一方面由于信息冗余也增加了数据处理的复杂度和计算开支。
我们经常需要通过一系列特征属性来对研究对象进行分析。特征属性越多越有利于刻画和分析事物，但属性越多带来的数据计算量也会越大，这就需要在两者之间取得一个平衡点。既能够减少研究对象的特征维度从而降低计算复杂度，又能够尽量减少因为降维导致的信息损失，就是数据降维的关键和首要目标。
主成分分析（Principal Component Analysis PCA）就是机器学习中一种常用且有效的数据降维方法。
PCA降维的核心思想
数据降维是因为高维数据中部分维度之间相关性较高造成了数据冗余，但数据降维的过程中又不可避免地造成部分信息的损失。如何做到既能够通过数据降维降低原始数据的冗余程度，又能够尽可能多地保留原始数据的信息呢？这其实就是PCA的核心思想。
核心思想一：基变化向量投影
要实现数据降维，一个常见的方法就是，高维向量向低维空间的投影。
1.向量投影与向量内积
向量内积的代数定义和几何定义表明：一个向量在另一个单位向量上的投影值，就是这两个向量的内积数值。如：aaa在模长为111的向量bbb投影为∣a∣cosθ|a|cosθ∣a∣cosθ,即a⋅ba·ba⋅b。
（1）向量内积的代数定义
两个向量内积的运算规则是，参与向量内积的两个向量必须维度相等，向量内积运算时将两个向量对应位置上的元素分别相乘之后求和即可得到向量内积的结果。向量内积结果是一个标量。例如，假设有两个维度相同的向量：
a=[a1a2a3⋮an]和b=[b1b2b3⋮bn]\begin{gathered} a= \begin{bmatrix} a_1\\a_2\\a_3\\\vdots\\a_n\end{bmatrix} \end{gathered}\begin{gathered} 和b= \begin{bmatrix} b_1\\b_2\\b_3\\\vdots\\b_n\end{bmatrix} \end{gathered}a=⎣⎢⎢⎢⎢⎢⎡a1a2a3⋮an⎦⎥⎥⎥⎥⎥⎤和b=⎣⎢⎢⎢⎢⎢⎡b1b2b3⋮bn⎦⎥⎥⎥⎥⎥⎤
即这两个向量的内积为：
a⋅b=a1b1+a2b2+a3b3⋯anbna·b=a_1b_1+a_2b_2+a_3b_3\cdots a_nb_na⋅b=a1b1+a2b2+a3b3⋯anbn
（2）向量内积的几何定义
向量内积的几何定义用来表征某个向量aaa在另一个向量bbb方向上的投影长度乘以向量b的模长，即a⋅b=∣a∣∣b∣cosθa·b=|a||b|cosθa⋅b=∣a∣∣b∣cosθ。如果向量bbb是单位向量（即模长为111），那么向量aaa与向量bbb的内积结果就等于向量a在向量b方向上的投影长度∣a∣cosθ|a|cosθ∣a∣cosθ这就是向量内积的几何定义。
2.向量投影与矩阵向量乘法
前文已经讲过，向量的表示依赖于基向量的选择。同一个向量，选择的标准正交基不同，所得到的向量坐标表示也不同。
（1）二维空间中
二维空间中一个向量实际上位于以两个标准正交向量与向量为[42]\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]实际上位于两个标准正交向量[01]\begin{gathered} \begin{bmatrix} 0\\1\end{bmatrix}\end{gathered}[01]与向量[10]\begin{gathered}\begin{bmatrix}1\\0\end{bmatrix}\end{gathered}[10]为基底所张成的二维空间，并且向[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]的坐标值是分别在标准正交基向量上的投影值。例如，向量[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在标准正交基向量[10]\begin{gathered}\begin{bmatrix} 1\\0\end{bmatrix}\end{gathered}[10]的投影值，就是两个向量的内积[10]⋅[42]\begin{gathered}\begin{bmatrix} 1\\0\end{bmatrix}\end{gathered}·\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[10]⋅[42]=4。这个投影值就是向量[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在x轴上的坐标值。向量[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在标准正交基向量[01]\begin{gathered}\begin{bmatrix} 0\\1\end{bmatrix}\end{gathered}[01]的投影值，就是两个向量的内积[01]⋅[42]\begin{gathered}\begin{bmatrix} 0\\1\end{bmatrix}\end{gathered}·\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[01]⋅[42]=2。这个投影值就是向量[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在y轴上的坐标值。上述向量可以写成矩阵向量乘积[1001]\begin{gathered}\begin{bmatrix} 1&0\\0&1\end{bmatrix}\end{gathered}[1001][42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]=[[10][01]⋅[42]⋅[42]]\begin{gathered}\begin{bmatrix} \begin{gathered}\begin{bmatrix} 1\\0\end{bmatrix}\\\begin{gathered} \begin{bmatrix} 0\\1\end{bmatrix}\end{gathered}\end{gathered}\begin{gathered}· \begin{bmatrix} 4\\2\end{bmatrix}\\\begin{gathered}· \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}\end{gathered} \end{bmatrix}\end{gathered}⎣⎢⎢⎢⎡[10][01]⋅[42]⋅[42]⎦⎥⎥⎥⎤=[42]\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]。
对于二维空间中的一个向量[42]\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]，如果选择两个标准正交向量[1212]\begin{gathered} \begin{bmatrix} \cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}⎣⎢⎢⎡2121⎦⎥⎥⎤和[−1212]\begin{gathered} \begin{bmatrix} -\cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}⎣⎢⎢⎡−2121⎦⎥⎥⎤作为基底，则向量[42]\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在标准正交基向量[1212]\begin{gathered} \begin{bmatrix} \cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}⎣⎢⎢⎡2121⎦⎥⎥⎤上的投影值为[1212]⋅[42]\begin{gathered} \begin{bmatrix} \cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}\cdot\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}⎣⎢⎢⎡2121⎦⎥⎥⎤⋅[42]=323\sqrt 232，向量[42]\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]在标准正交基向量[−1212]\begin{gathered} \begin{bmatrix} -\cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}⎣⎢⎢⎡−2121⎦⎥⎥⎤上的投影值为[−1212]⋅\begin{gathered} \begin{bmatrix} -\cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}\cdot⎣⎢⎢⎡−2121⎦⎥⎥⎤⋅[42]=−2\begin{gathered} \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}=-\sqrt 2[42]=−2。上述向量可以写成矩阵向量乘积[12−121212]\begin{gathered} \begin{bmatrix} \cfrac{1}{\sqrt 2}&-\cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}&\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}⎣⎢⎢⎡2121−2121⎦⎥⎥⎤[42]\begin{gathered}\begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}[42]=[[1212][−1212]⋅[42]⋅[42]]\begin{gathered}\begin{bmatrix} \begin{gathered}\begin{bmatrix} \cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\\\begin{gathered} \begin{bmatrix} -\cfrac{1}{\sqrt 2}\\\cfrac{1}{\sqrt 2}\end{bmatrix}\end{gathered}\end{gathered}\begin{gathered}· \begin{bmatrix} 4\\2\end{bmatrix}\\\begin{gathered}· \begin{bmatrix} 4\\2\end{bmatrix}\end{gathered}\end{gathered} \end{bmatrix}\end{gathered}⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡⎣⎢⎢⎡2121⎦⎥⎥⎤⎣⎢⎢⎡−2121⎦⎥⎥⎤⋅[42]⋅[42]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=[32−2]\begin{gathered} \begin{bmatrix} 3\sqrt 2\\-\sqrt 2\end{bmatrix}\end{gathered}[32−2]。
（2）推广到m维空间
mmm维空间向量[x1x2x3⋮xm]\begin{gathered}\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xm⎦⎥⎥⎥⎥⎥⎤实际上以mmm个标准正交基向量[100⋮0]\begin{gathered}\begin{bmatrix}1\\0\\0\\\vdots\\0\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡100⋮0⎦⎥⎥⎥⎥⎥⎤，[010⋮0]\begin{gathered}\begin{bmatrix}0\\1\\0\\\vdots\\0\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡010⋮0⎦⎥⎥⎥⎥⎥⎤，[001⋮0]\begin{gathered}\begin{bmatrix}0\\0\\1\\\vdots\\0\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡001⋮0⎦⎥⎥⎥⎥⎥⎤，⋯\cdots⋯,[000⋮1]\begin{gathered}\begin{bmatrix}0\\0\\0\\\vdots\\1\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡000⋮1⎦⎥⎥⎥⎥⎥⎤为基底所张成mmm维空间中，并且坐标值是向量[x1x2x3⋮xm]\begin{gathered}\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xm⎦⎥⎥⎥⎥⎥⎤分别在标准正交基向量上的投影值。上述向量投影可写成矩阵向量乘积形式，即标准正交基向量作为行向量所构成的m×mm\times mm×m矩阵与向量相乘，具体如下。
[100⋯0010⋯0001⋯0⋮⋮⋮⋱⋮000⋯1][x1x2x3⋮xm]=[x1x2x3⋮xm]\begin{bmatrix} 1 & 0 &0 & \cdots & 0 \\ 0 & 1 &0 & \cdots & 0 \\ 0 & 0 &1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots &\vdots \\ 0 & 0 &0 & \cdots & 1 \\ \end{bmatrix} \begin{gathered}\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}= \begin{gathered}\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered} ⎣⎢⎢⎢⎢⎢⎡100⋮0010⋮0001⋮0⋯⋯⋯⋱⋯000⋮1⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xm⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xm⎦⎥⎥⎥⎥⎥⎤
由于选择的标准正交基向量不同，向量的坐标表示也不同。对于某个mmm维向量[x1x2x3⋮xm]\begin{gathered}\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xm⎦⎥⎥⎥⎥⎥⎤，想将其变换到由kkk个mmm维标准正交基向量表示的新空间中，那么首先将kkk个mmm维标准正交基向量按行组成矩阵QQQ，最终矩阵向量乘积Q⋅xQ\cdot xQ⋅x就是变换结果。
推广来看，如果有nnn个mmm维列向量aja_jaj，要将其变换到由kkk个mmm维标准正交列向量qiq_iqi作为基底表示的新空间中，首先需要将kkk个mmm维标准正交基向量按行组成矩阵QQQ，然后将nnn个mmm维向量按列组成矩阵AAA，那么两个矩阵的乘积Q⋅AQ\cdot AQ⋅A就是变换结果，如下所示：
[q1Tq2Tq3T⋮qkT][a1a2a3⋯an]=[q1Ta1q1Ta2q1Ta3⋯q1Tanq2Ta1q2Ta2q2Ta3⋯q2Tanq3Ta1q3Ta2q3Ta3⋯q3Tan⋮⋮⋮⋱⋮qkTa1qkTa2qkTa3⋯qkTan]\begin{gathered}\begin{bmatrix}q_1^T\\q_2^T\\q_3^T\\\vdots\\q_k^T\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}\begin{gathered} \begin{bmatrix}a_1&a_2&a_3&\cdots&a_n\end{bmatrix}\end{gathered}\begin{gathered}\end{gathered}=\begin{bmatrix} q_1^Ta_1 & q_1^Ta_2 &q_1^Ta_3 & \cdots & q_1^Ta_n \\ q_2^Ta_1 & q_2^Ta_2 &q_2^Ta_3 & \cdots & q_2^Ta_n \\ q_3^Ta_1 & q_3^Ta_2 &q_3^Ta_3 & \cdots & q_3^Ta_n \\ \vdots & \vdots & \vdots & \ddots &\vdots \\ q_k^Ta_1 & q_k^Ta_2 &q_k^Ta_3 & \cdots & q_k^Ta_n \\ \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎡q1Tq2Tq3T⋮qkT⎦⎥⎥⎥⎥⎥⎤[a1a2a3⋯an]=⎣⎢⎢⎢⎢⎢⎡q1Ta1q2Ta1q3Ta1⋮qkTa1q1Ta2q2Ta2q3Ta2⋮qkTa2q1Ta3q2Ta3q3Ta3⋮qkTa3⋯⋯⋯⋱⋯q1Tanq2Tanq3Tan⋮qkTan⎦⎥⎥⎥⎥⎥⎤
其中qiq_iqi是mmm维标准正交列向量，因此qiTq_i^TqiT是一个行向量，表示第i个基向量；aja_jaj是一个mmm维列向量，表示第jjj个原始数据记录。
这里的k(k≤m)k(k\le m)k(k≤m)决定了变换后数据的维度，也就是说可以将mmm维数据aja_jaj变换到更低维度的空间中，变换后的维度取决于基向量的数量kkk，因此矩阵相乘也可以表示降维变换。
上述矩阵相乘的几何意义就是，两个矩阵相乘的结果是将右边矩阵中的每一个列向量变换到以左边矩阵中每一个行向量为基底所表示的空间中。也就是前文所说的，一个矩阵代表着一种线性变换。因此，要实现对矩阵AAA的线性变换可以考虑左乘一个以标准正交基为行向量的矩阵QQQ。
核心思想二：协方差归零投影
既然数据降维的起因是高维数据的维度之间存在较高的相关性导致数据信息存在冗余，那么数据降维的一个核心思想自然就是，数据降维后的维度之间尽可能相对独立，也就是降维之后的数据维度之间的协方差为0。
此处我们回顾一下方差与协方差的基础知识。
1.方差
方差和标准差是最常用的度量一组数据分散程度的指标。对于一组含有nnn个样本的集合，我们容易知道以下公式。
均值：
xˉ=∑i=1nxin\bar x= \frac {\sum_{i=1}^n x_i}{n}xˉ=n∑i=1nxi
方差：
S2=∑i=1n(xi−xˉ)2n−1S^2= \frac {\sum_{i=1}^n (x_i-\bar x)^2}{n-1}S2=n−1∑i=1n(xi−xˉ)2
标准差：
S=∑i=1n(xi−xˉ)2n−1S=\sqrt\frac {\sum_{i=1}^n (x_i-\bar x)^2}{n-1}S=n−1∑i=1n(xi−xˉ)2
例如两组数据为[1,4,8,13,24][1,4,8,13,24][1,4,8,13,24]和[8,9,10,11,12][8,9,10,11,12][8,9,10,11,12]且它们的均值都是101010，但这两组数据的分布情况差别很大。前者标准差为9.039.039.03，后者标准差为1.581.581.58。也就是说，虽然两组数据均值相等，但是后者相对前者数据分布上更为“拥挤”和“集中”。方差和标准差描述的就是这种“散布度”。
一组数据如果越“拥挤”，数据点就越接近，每个数据点与均值的差距(xi−xˉ)2(x_i-\bar x)^2(xi−xˉ)2就越小，对应的方差和标准差就越小。所以方差和标准差能够较好地反映一组数据的分散程度。
2.协方差
协方差度量的是维度和维度之间的关系。例如我们收集到某年级所有学生的语文、历史、地理、英语、数学、物理、化学、生物等科目的考试成绩，对于这样的数据集我们可以对每一维数据也就是每个科目数据进行独立的方差计算，但通常我们还想了解这些科目成绩之间的关系。于是，我们就需要用到协方差这个概念。
假设两组数据分别是xxx和yyy，那么这两组数据的协方差为：
cov(x,y)=∑i=1n(xi−xˉ)(yi−yˉ)n−1cov(x,y)= \frac {\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)}{n-1}cov(x,y)=n−1∑i=1n(xi−xˉ)(yi−yˉ)
上述协方差公式的通俗理解就是，协方差表达了两个变量在变化过程中变化的方向一致性和变化大小的程度。例如两个变量是同方向变化，还是反方向变化？两个变量同向或反向变化的程度如何？
xxx变大，同时yyy也变大，说明两个变量xxx和yyy是同向变化的，这时协方差就为正。协方差数值越大，两个变量同向变化的程度也就越高。
xxx变大，但同时yyy变小，说明两个变量xxx和yyy是反向变化的，这时协方差就为负。协方差数值（绝对值）越大，两个变量反向变化的程度也就越高。
xxx变化趋势和yyy变化趋势相互独立的时候，协方差就为000。
总的来说，协方差是度量各个维度偏离其均值程度的一个指标。协方差为正说明两者是正相关的，协方差为负说明两者是负相关的，协方差为000说明两者的关系就是统计上说的“相互独立”。
3.协方差矩阵
方差和标准差主要用来处理一维数据，协方差只能处理二维数据，那么多维数据怎么办呢？多维数据就需要多次计算协方差，也就是将多维数据中的维度数据两两计算协方差，例如一个101010维数据就需要计算45(C102=45)45(C_{10}^2=45)45(C102=45)个不同维度间的协方差。
如此多的协方差该如何处理呢？这个时候矩阵就派上用场了，我们使用矩阵来组织这些协方差，即构建协方差矩阵。协方差矩阵就是度量维度和维度之间关系的矩阵。
假设向量X=(X1,X2,⋯，Xn)TX=(X_1,X_2,\cdots，X_n)^TX=(X1,X2,⋯，Xn)T为nnn维随机变量，则称矩阵：
C=(Cij)n×n=[C11C12C13⋯C1nC21C22C23⋯C2nC31C32C33⋯C3n⋮⋮⋮⋱⋮Cm1Cm2Cm3⋯Cmn]C=(C_{ij})_{n\times n}=\begin{bmatrix} C_{11} & C_{12} &C_{13} & \cdots & C_{1n} \\ C_{21} & C_{22} &C_{23} & \cdots & C_{2n} \\ C_{31} & C_{32} &C_{33} & \cdots & C_{3n} \\ \vdots & \vdots & \vdots & \ddots &\vdots \\ C_{m1} & C_{m2} &C_{m3} & \cdots & C_{mn} \\ \end{bmatrix}C=(Cij)n×n=⎣⎢⎢⎢⎢⎢⎡C11C21C31⋮Cm1C12C22C32⋮Cm2C13C23C33⋮Cm3⋯⋯⋯⋱⋯C1nC2nC3n⋮Cmn⎦⎥⎥⎥⎥⎥⎤
为nnn维随机变量XXX的协方差矩阵，其中Cij=cov(Xi,Xj),(i,j=1,2,⋯,n)C_{ij}=cov(X_i,X_j),(i,j=1,2, \cdots ,n )Cij=cov(Xi,Xj),(i,j=1,2,⋯,n)为XXX的分量XiX_iXi和XjX_jXj的协方差。
例如，某年级学生成绩单上有语文、数学、物理333门课程的成绩，记作x=x=x=语文、y=y=y=数学、z=z=z=物理。那么这333门课程的协方差矩阵如下：
C=[cov(x,x)cov(x,y)cov(x,z)cov(y,x)cov(y,y)cov(y,z)cov(z,x)cov(z,y)cov(z,z)]C=\begin{bmatrix} cov(x,x)& cov(x,y) &cov(x,z) \\ cov(y,x) & cov(y,y) &cov(y,z) \\ cov(z,x) & cov(z,y) &cov(z,z) \\ \end{bmatrix}C=⎣⎡cov(x,x)cov(y,x)cov(z,x)cov(x,y)cov(y,y)cov(z,y)cov(x,z)cov(y,z)cov(z,z)⎦⎤
而且不难得知，协方差矩阵是一个对称矩阵，主对角线元素就是各个维度的方差，非主对角线元素就是不同维度之间的协方差。
所以，我们降维之后希望各个维度之间相互独立，也就是希望降维之后不同维度之间的协方差为000，同样也就是希望上面的协方差矩阵除了主对角线之外的部分都为000。
核心思想三：最大方差投影
我们知道二维空间数据降维的方法是将数据点投影到某条直线上，那么应该选择哪条直线才合理呢？只要进行数据降维，将二维空间中的数据点投影到一维空间的某条直线上，就一定会带来信息的损失，这是无法避免的。所以我们考虑选择哪条直线来投影比较合理其实就是考虑选择哪条直线来投影带来的信息损失最小。
那二维空间中数据点最重要的信息是什么呢？我们有理由认为二维空间中数据点最重要的信息就是这些数据点之间的分布规律，这种分布规律必然要通过数据点的“差异”来体现。同时，我们要注意从高维空间向低维空间降维的过程中，数据点一定会变得更加“拥挤”。例如二维空间的数据点[21]\begin{bmatrix}2\\1\end{bmatrix}[21]和[41]\begin{bmatrix}4\\1\end{bmatrix}[41]如果向yyy轴进行投影，那么得到的投影数值都是111。原来二维空间的222个数据点都变变成了一维空间（直线）上的111个数值点，变得更加“拥挤”了。投影效果展示如图所示：
既然降维过程中天然存在将高维空间中的数据点在低维空间中变得更加“拥挤”从而造成信息损失的倾向，那么合理、科学的降维投影方式就应该是使投影之后的数据点尽可能分散，尽可能降低这种降维带来的“拥挤”程度。
这种投影数据点的分散实际上就是要求原始数据矩阵降维处理之后的新矩阵的维度的方差尽可能大，也就是降维之后矩阵的协方差矩阵的对角线元素尽可能大。我们可以将这种降维投影的要求称为“最大方差投影”。
降维关键：协方差矩阵对角化
PCA降维本质上是通过对矩阵AAA进行线性变换来实现的，但是这种降维并不是随意的，而是要求降维之后的矩阵YYY能够最大程度地保持原有矩阵的性质，也就是原始数据的失真程度应尽可能小。
如何才能保证降维过程中数据失真尽可能小呢？我们知道协方差矩阵度量的是维度和维度之间的关系，协方差矩阵主对角线上的元素就是各个维度上的方差，非主对角线上的元素就是各维度之间的相关性（协方差）。一个合理的降维过程应该满足 “协方差归零投影”和“最大方差投影”的要求，也就是，降维之后新矩阵YYY的协方差矩阵CyC_yCy的非主对角线元素尽可能为000，而主对角线元素尽可能大。满足上述要求的矩阵是一个对角矩阵，所以降维的实质就是要求降维之后的新矩阵YYY的协方差矩阵CyC_yCy是对角矩阵。
详解PCA降维
计算矩阵YYY的协方差矩阵CyC_yCy
1.为何计算协方差矩阵
PCA降维的原理推导过程中，我们通过对原矩阵AAA进行线性变换得到。紧接着，我们就开始计算矩阵YYY的协方差矩阵。那么，我们为什么会考虑计算矩阵YYY的协方差矩阵CyC_yCy呢？
这是因为矩阵YYY的协方差矩阵CyC_yCy的主对角线元素是降维后新维度的方差，非主对角线元素是降维后各新维度的协方差。而PCA降维的核心思想就是“协方差归零投影”和“最大方差投影”，也就是要求降维之后的矩阵Y=QAY=QAY=QA的协方差为000，而维度方差尽可能大。
矩阵YYY的协方差矩阵CyC_yCy的主对角线元素就是新维度的方差，而非主对角线元素就是新维度之间的协方差。所以，一个合理的降维过程对应的新矩阵YYY的协方差矩阵应该是一个对角矩阵。
2.详解协方差矩阵的表达式
（1）Cy=1n−1YYTC_y=\frac{1}{n-1}YY^TCy=n−11YYT
设矩阵X=[a1a2a3a4a5b1b2b3b4b5]\begin{gathered}X=\begin{bmatrix} a_1&a_2&a_3&a_4&a_5\\b_1&b_2&b_3&b_4&b_5\end{bmatrix} \end{gathered}X=[a1b1a2b2a3b3a4b4a5b5]，XXX的转置矩阵XT=[a1b1a2b2a3b3a4b4a5b5]\begin{gathered}X^T=\begin{bmatrix} a_1&b_1\\a_2&b_2\\a_3&b_3\\a_4&b_4\\a_5&b_5\end{bmatrix} \end{gathered}XT=⎣⎢⎢⎢⎢⎡a1a2a3a4a5b1b2b3b4b5⎦⎥⎥⎥⎥⎤ 。XXTXX^TXXT为：
[a12+a22+a32+a42+a52a1b1+a2b2+a3b3+a4b4+a5b5a1b1+a2b2+a3b3+a4b4+a5b5b12+b22+b32+b42+b52]\begin{bmatrix} a_1^2+a_2^2+a_3^2+a_4^2+a_5^2 & a_1b_1+a_2b_2+a_3b_3+a_4b_4+a_5b_5 \\ a_1b_1+a_2b_2+a_3b_3+a_4b_4+a_5b_5 & b_1^2+b_2^2+b_3^2+b_4^2+b_5^2 \end{bmatrix}[a12+a22+a32+a42+a52a1b1+a2b2+a3b3+a4b4+a5b5a1b1+a2b2+a3b3+a4b4+a5b5b12+b22+b32+b42+b52]
观察上面矩阵中的各个元素，不难发现以下结论：
第一，根据方差公式：S2=∑i=1n(xi−xˉ)2n−1S^2= \frac {\sum_{i=1}^n (x_i-\bar x)^2}{n-1}S2=n−1∑i=1n(xi−xˉ)2 当对矩阵XXX进行零均值化处理后，方差公式就简化为：S2=∑i=1nxi2n−1S^2= \frac {\sum_{i=1}^n x_i^2}{n-1}S2=n−1∑i=1nxi2
于是，上述XXTXX^TXXT结果矩阵的主对角线元素正好就是矩阵X “维度”的方差的(n-1)倍。
第二，根据协方差公式：
cov(x,y)=∑i=1n(xi−xˉ)(yi−yˉ)n−1cov(x,y)= \frac {\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)}{n-1}cov(x,y)=n−1∑i=1n(xi−xˉ)(yi−yˉ)
当对矩阵XXX进行零均值化处理后，协方差公式就简化为：
cov(x,y)=∑i=1nxiyin−1cov(x,y)= \frac {\sum_{i=1}^n x_iy_i}{n-1}cov(x,y)=n−1∑i=1nxiyi
于是，上述XXTXX^TXXT结果矩阵的非主对角线元素正好就是矩阵 XXX“维度”的协方差的(n−1)(n-1)(n−1)倍。
总结，我们将上述XXTXX^TXXT结果矩阵乘以1/（n−1）1/（n-1）1/（n−1）正好就是矩阵XXX的协方差矩阵CxC_xCx。也就是说，一个对行数据进行归零化处理的矩阵XXX，其协方差矩阵为：
Cx=1n−1XXTC_x=\frac{1}{n-1}XX^TCx=n−11XXT
由于原始数据的矩阵AAA是经过行数据归零化处理的，因此有：C=1n−1AATC=\frac{1}{n-1}AA^TC=n−11AAT
那么，矩阵Y=QAY=QAY=QA是否也经过行数据归零化处理了呢？假设行数据归零化后的矩阵A=[a1a2a3⋯am]A=[a_1a_2 a_3\cdots a_m]A=[a1a2a3⋯am]，则a1+a2+a3+⋯am=0a_1+a_2+a_3+\cdots a_m=0a1+a2+a3+⋯am=0。于是Y=QA=[Qa1Qa2Qa3⋯Qam]Y=QA=[Qa_1 Qa_2Qa_3 \cdots Qa_m]Y=QA=[Qa1Qa2Qa3⋯Qam]，则也Qa1+Qa2+Qa3+⋯+Qam=0Qa_1+Qa_2+Qa_3+ \cdots +Qa_m=0Qa1+Qa2+Qa3+⋯+Qam=0。也就是说，矩阵Y=QAY=QAY=QA的行数据也是符合归零化条件的。因此，可以知道矩阵YYY的协方差矩阵：Cy=1n−1YYTC_y=\frac{1}{n-1}YY^TCy=n−11YYT
（2）Cy=1n−1YYT=QCQTC_y=\frac{1}{n-1}YY^T=QCQ^TCy=n−11YYT=QCQT
由于Y=QAY=QAY=QA，因此Cy=1n−1YYT=1n−1QA(QA)TC_y=\frac{1}{n-1}YY^T=\frac{1}{n-1}QA(QA)^TCy=n−11YYT=n−11QA(QA)T。由转置矩阵性质可知(QA)T=ATQT(QA)^T=A^TQ^T(QA)T=ATQT，所以上式可写作：
Cy=1n−1YYT=1n−1QAATQTC_y=\frac{1}{n-1}YY^T=\frac{1}{n-1}QAA^TQ^TCy=n−11YYT=n−11QAATQT
观察发现其中1n−1AAT\frac{1}{n-1}AA^Tn−11AAT就是矩阵A的协方差矩阵CCC，因此上式可以写作为：
Cy=1n−1YYT=QCQTC_y=\frac{1}{n-1}YY^T=QCQ^TCy=n−11YYT=QCQT
矩阵YYY的协方差矩阵CyC_yCy对角化
数据降维后新矩阵YYY的各个新维度要线性无关，也就是要求矩阵YYY的协方差矩阵Cy=QCQTC_y=QCQ^TCy=QCQT能够实现对角化。
1.什么样的矩阵QQQ能够对角化CyC_yCy
根据实对称矩阵正交对角化的定理可以得到：
Λ=UTCUΛ=U^TCUΛ=UTCU
所以，要使矩阵YYY的协方差矩阵Cy=QCQTC_y=QCQ^TCy=QCQT是对角矩阵，只需要Q=UTQ=U^TQ=UT即可。这个结果表明，矩阵QQQ如果是由原矩阵AAA的协方差矩阵CCC的特征向量构成的矩阵，矩阵AAA经过矩阵QQQ线性变换之后的矩阵YYY的协方差矩阵CyC_yCy就为对角矩阵。
2.实对称矩阵对角化性质
实对称矩阵有一个非常好的性质，那就是其可以转化为对角矩阵。原矩阵的协方差矩阵CCC满足实对称矩阵的条件，所以可以通过线性变换将CCC转化为对角矩阵ΛΛΛ，具体来说就是：
A=UTCUA=U^TCUA=UTCU
举个例子，假设有实对称矩阵：
A=[1−20−22−20−23]A=\begin{bmatrix} 1& -2 &0 \\ -2 & 2 &-2 \\ 0 & -2 &3 \\ \end{bmatrix}A=⎣⎡1−20−22−20−23⎦⎤
我们求解得到它的特征值为555、222和−1-1−1。单位正交化后的特征向量分别为[0.333−0.6670.667]\begin{bmatrix} 0.333 \\ -0.667\\ 0.667 \\ \end{bmatrix}⎣⎡0.333−0.6670.667⎦⎤、[−0.6670.3330.667]\begin{bmatrix} -0.667 \\ 0.333 \\ 0.667 \\ \end{bmatrix}⎣⎡−0.6670.3330.667⎦⎤和[0.6670.3330.667]\begin{bmatrix} 0.667 \\ 0.333 \\ 0.667 \\ \end{bmatrix}⎣⎡0.6670.3330.667⎦⎤。由这些单位正交的特征向量组成的矩阵就是使成立的矩阵：
U=[0.333−0.667−0.667−0.6670.333−0.6670.6670.667−0.333]U=\begin{bmatrix} 0.333& -0.667 & -0.667\\ -0.667& 0.333 & -0.667\\ 0.667& 0.667 & -0.333\\ \end{bmatrix}U=⎣⎡0.333−0.6670.667−0.6670.3330.667−0.667−0.667−0.333⎦⎤
就是使Λ=UTAUΛ=U^TAUΛ=UTAU成立的矩阵UUU我们可以验证一下 ：
UTAU=[50002000−1]U^TAU=\begin{bmatrix} 5& 0 & 0\\ 0& 2 & 0\\ 0& 0 & -1\\ \end{bmatrix}UTAU=⎣⎡50002000−1⎦⎤
观察上式可以发现，由实对称矩阵AAA的单位正交的特征向量组成的矩阵就是使矩阵AAA对角化的矩阵UUU，并且矩阵AAA对角化后的主对角线元素就是矩阵AAA的特征值。
总结来说，原矩阵AAA的协方差矩阵CCC是一个实对称矩阵，因此一定可以找到单位正交矩阵UUU使得Λ=UTCUΛ=U^TCUΛ=UTCU。而我们通过PCA降维的原理推导过程已经知道矩阵YYY的协方差矩阵Cy=QCQTC_y=QCQ^TCy=QCQT。
我们降维的目标是使Cy=QCQTC_y=QCQ^TCy=QCQT是对角矩阵ΛΛΛ，那么如何使Cy=QCQTC_y=QCQ^TCy=QCQT是对角矩阵ΛΛΛ呢？通过对比，我们发现只需要Q=UTQ=U^TQ=UT就可以实现对角化。也就是说，我们寻找的矩阵QQQ就是矩阵UTU^TUT，也就是由协方差矩阵CCC的单位正交的特征向量组成的矩阵的转置矩阵。
3.求解降维矩阵PPP
通过前文的讲解，我们知道PCA推导的逻辑链条如下：
(1)对原矩阵A(m×n)A(m×n)A(m×n)进行降维就是对矩阵AAA进行某种线性变换，也就等效于寻找某个合适的矩阵QQQ来乘以矩阵AAA，从而得到降维后的新矩阵Y=QAY=QAY=QA。
(2)这个合适的矩阵QQQ应该满足一些条件，具体来说就是“协方差归零投影”和“最大方差投影”。其中“协方差归零投影”也就是希望矩阵YYY的协方差矩阵CyC_yCy是对角矩阵。
(3)但是“最大方差投影”体现在哪里呢？我们知道原矩阵A(m×n)A(m×n)A(m×n)是mmm维，而降维矩阵P(k×m)P(k×m)P(k×m)是kkk维，kkk取值不同对应的降维程度就不同。“最大方差投影”要求我们降维的时候，首先选取协方差矩阵CyC_yCy对角化矩阵主对角线元素中最大值（也就是协方差矩阵CCC的特征值的最大值）所对应的特征向量，依此类推。这才是我们最后对应某个kkk值的降维矩阵PPP。
| PCA降维的步骤 |
|---|
| 计算原矩阵AAA的协方差矩阵CCC。 |
| 计算协方差矩阵CCC的单位正交的特征向量与对应的特征值。 |
| 根据降维要求，确定kkk值大小。将CCC的特征值从大到小排列，选取前kkk个特征值所对应的特征向量。 |
| 将这些特征向量作为行向量，求解出降维矩阵PPP。 |
| 将降维矩阵P乘以原矩阵AAA即可降维，得到Y=PAY=PAY=PA。 |
Matlab代码实现
load fisheriris; %导入数据集
X = meas; % n = 150, m = 4
meanX = ones(size(X,1), 1) * mean(X); % 中心化处理centredX = X - meanX;C = cov(centredX); % 直接调用cov直接计算协方差矩阵即可[W, Lambda] = eig(C); % W是特征向量组成的矩阵（4×4），Lambda是特征值组成的对角矩阵
ev = (diag(Lambda))'; % 提取特征值
ev = ev(:, end:-1:1); % eig计算出的特征值是升序的，这里手动倒序（W同理）
W = W(:, end:-1:1);
sum(W.*W, 1) % 可以验证每个特征向量各元素的平方和均为Wr = W(:, 1:2); % 提取前两个主成分的特征向量
Tr = centredX * Wr; % 新坐标空间的数据点
% 作图
figure;stairs(cumsum(ev)/sum(ev), 'LineWidth',1.5);axis([1 4 0 1]);xlabel('$ k $', 'Interpreter', 'latex');ylabel('$ f(k)=\frac{\sum _{i=1}^i \lambda_k}{\sum_{i=1}^m \lambda_i} $',...'Interpreter', 'latex');hold on;plot([1 4], [0.95 0.95], '--'); % 从图中可以看出，r为方差贡献率，取r = 2
figure;scatter(Tr(:,1), Tr(:,2), 130, categorical(species), '.');colormap(winter);xlabel('Principal Component 1');ylabel('Principal Component 2');
[U, Sigma, V] = svd(X); % 可以检验，W和V完全相同（向量的正负号不影响）
Vr = V(:, 1:2); % 提取前两个主成分的特征向量
Tr = X * Vr; % 新坐标空间的数据点
% 画图部分同上
[loadings, scores] = pca(X, 'NumComponents', r);
[Wr, Tr, ev] = pca(X, 'NumComponents',2); % 画图部分
【主成分分析】PCA降维算法及Matlab代码实现相关推荐
- PCA降维算法原理及代码实现（python和matlab）
常见的数据降维算法有:奇异值分解(SVD).主成分分析(PCA).因子分析(FA).独立成分分析(ICA). PCA降维的基本思想:通过计算数据矩阵的协方差矩阵,然后得到协方差矩阵的特征值.特征向量. ...
- PCA 降维算法详解 以及代码示例
1. 前言 PCA : principal component analysis ( 主成分分析) 最近发现我的一篇关于PCA算法总结以及个人理解的博客的访问量比较高, 刚好目前又重新学习了一下PCA ...
- 基于PCA主成分分析的BP神经网络回归预测MATLAB代码
基于PCA主成分分析的BP神经网络回归预测MATLAB代码 代码注释清楚. 先对数据集进行主成分分析,自主根据贡献率选择主成分:同时计算KMO验证值:用PCA以后数据进行BP神经网络回归预测. 可以读 ...
- 【视频】主成分分析PCA降维方法和R语言分析葡萄酒可视化实例|数据分享
最近我们被客户要求撰写关于主成分分析PCA的研究报告,包括一些图形和统计输出.降维技术之一是主成分分析 (PCA) 算法,该算法将可能相关变量的一组观察值转换为一组线性不相关变量.在本文中,我们将讨论 ...
- PSO-LSSVM算法及其MATLAB代码
挺完整的一篇博客,这里转载记录一下. 原文链接:PSO-LSSVM算法及其MATLAB代码 一.PSO 1.概念 粒子群优化算法(PSO:Particle swarm optimization)是一种 ...
- 哈里斯鹰优化(HHO)算法(含MATLAB代码)
先做一个声明:文章是由我的个人公众号中的推送直接复制粘贴而来,因此对智能优化算法感兴趣的朋友,可关注我的个人公众号:启发式算法讨论.我会不定期在公众号里分享不同的智能优化算法,经典的,或者是近几年提出 ...
- 直流电机，传递函数，模糊控制pid算法，matlab代码
直流电机,传递函数,模糊控制pid算法,matlab代码 ID:1466674483342194夹竹丶桃
- 最短路径 Dijkstra算法的Matlab代码实现
为了搞清楚最短路径的算法过程,自己编写代码实现dijkstra算法寻找路径 % 文件名:dijkstra.m % 时间:2020年9月12日 % 来源:https://blog.csdn.net/li ...
- 下料问题的列生成算法，matlab代码，简单易懂
下料问题的列生成算法,matlab代码,简单易懂 ID:18300606016240721绿兔叽
最新文章
- 道歉无用！被 Linux “拉黑”的明尼苏达大学还在努力重获信任
- leetcode算法题--多米诺与托米诺平铺★
- Keras + Ubuntu环境搭建
- visualSVN仓库创建及操作
- pandas删除数据库 python_Python常见的科学计算库
- Linux Vim三种工作模式（命令模式、输入模式和编辑模式）详解
- 单片机测量代码运行时间方法-STM32
- 辨异 —— 逻辑之辨、人文社科观念
- 班尼路信息化系统基础选型的简单分析
- 二级python和office哪个难_对于操作office来说，python能与vbs相比吗？谁强谁弱呢？...
- 快逸报表数据库密码加密解决方案
- 写给电脑小白的电脑科普
- 宏观经济学——GDP
- 热插拔技术详解（上）
- 用MATLAB沉降观测实验,观测沉降论文,关于建筑物沉降观测数据处理相关参考文献资料-免费论文范文...
- 分门别类刷leetcode——贪心算法（C++实现）
- widows升级nodejs版本
- android自定义多选框 带图片,Android自定义单选多选下拉列表的实例代码
- JavaScript前端判断文件是否存在（案例详解）
- 关键词推广怎么做比较好？抖音宣传做关键词推广有哪些好的方法