“深入浅出数学建模”之主成分分析

主成分分析（PCA）

好多小伙伴在学习数学建模的时候会遇到这么一个算法—主成分分析，顾名思义，是用来分析主要成分的。简单点说，就是对评价指标集进行分析，得到新的若干个成分，来代替原始的评价指标（听不懂没关系！继续往下看，保证可以听得懂！）。

它是一种普遍的降维算法，在机器学习、生物分析等各个领域都有着广泛的应用。接下来，跟着我的思路，带你掌握这个算法，并学会如何去应用！

1 例题引入

大家设想一下这个例子，为了衡量一个人的体育综合能力，我们设有如下评价指标集{BMI（身高体重），肺活量，400米跑成绩，800米跑成绩，跳高，跳远}。接着我们分析问题的时候，觉得这些指标太多了（似乎并不多，但是假设你觉得多），对我们分析问题比较麻烦，这个时候我需要对指标进行处理，这个时候我就进行主成分分析。很多同学对这有个误解，认为主成分分析出来的是从原来指标集中直接挑出来的，绝对不是这样！大家请看下图：

在这个图中，左边是我们刚刚设定的变量，这里我用1-6来表示，右边是我们求出来的主成分（PC），为什么会有6个主成分呢？实际上我们用matlab软件求出来的就是多少个指标对应着多少个主成分。有人会问，不是说好的降维吗，为啥还是六个主成分呢？别急，我们先保留这个问题。

大家还需要注意的一点就是，每一个主成分和原来的指标都是存在一个线性关系的，我这里图中只画出来PC1，其他的PC也是这样的，也就说，我们可以把PC1用下面公式表示：
$$P{C}_1=a{X}_1+b{X}_2+...+f{X}_6$$
什么意思呢？也就是说，每一个主成分都和原来的变量有关系，并不是说主成分就是从原来变量中直接挑出来的，但是他是有一定意义的。比如说，假设这个主成分PC1在X2,X3,X4(对应着肺活量，400米跑，800米跑)前面的系数比较大且为正值，那我们可以分析得到这个主成分主要侧重于人的耐力方面（肺活量，400米，800米这三个指标本身就体现了一个人的耐力）；同理如果某个主成分PC2，在X5,X6,对着的前面系数比较大，说明我们分析得到的这个主成分是衡量人的爆发力方面。但是注意，并不是每一个主成分都可以分析出它们的意义。

我们可以看到前面两个主成分是彼此独立的，衡量人体育成绩不同方面的，这也是符合我们寻找主成分的定义的。接下来我们看一下，主成分分析的具体原理又是怎么实现的呢？

2 基本思想

在展开PCA数学原理之前，我们先整体把握一下其思想：

PCA算法是一种降维算法，降维就不可避免的涉及到数据的损失，就好比在三维世界中你看到的一个球，在二维世界中你只能看到一个面。怎么让我们的数据经过降维之后，尽可能的保留原始信息，这是我们需要解决的。

我们期望找到主成分彼此之间尽可能的独立，就比如前面例题中找到的主成分，分别来侧重说明某一方面，相互独立。如果不是相互独立，说明我们找的主成分不具有代表性，是不符合要求的。

总结一下：

• 降维后数据尽可能损失小
• 得到主成分之间尽可能独立

以上就是我们PCA算法的基本思想，想深入继续理解的同学可以继续往下看，不想深入了解的，直接想调用代码的可以直接略过下小节。

3 数学原理

在这里，我们不多赘述复杂的过程与推导，我们只关心理解算法，了解算法，知道算法是怎么来的，怎么解决的，有大体的了解就ok了。

我们首先有数据集如下图所示：

这里每一行代表一个特征（N），每一列代表一个样本（M），同时我们选择N个R维的正交基，如下图所示：

接下来我们会有：

大家对这一步不必过于纠结，这一步就是进行一个线性的变化，进行降维的过程，可以看到经过一个矩阵P进行处理后，我们得到的矩阵X’的维度发生了变化，每一行仍然代表是特征，但是此时的特征数变成R（R<N）,样本的数目没有变，达到了我们降维的目的。同时对于变化后的每一个变量，有：
$$p_i{x}_j=p_{i1}{x}_{1j}+p_{i2}{x}_{2j}+...+p_{iN}{x}_{NN}$$

是线性的关系，也就是为什么我们说PCA是一种线性降维。当然我们也有非线性降维，后续我们还会继续讨论这个问题。

我们回头再看我们的问题，我们可以找任何一个P矩阵对我们原始矩阵进行降维，但是怎么才能找到我们最好的，最优的P矩阵呢？此时问题转化为求P矩阵的最优问题，优化条件则是我们上文讲过的两个问题：一是降维后数据尽可能的损失小；二是降维得到的主成分尽可能的独立。

对于问题一，我们可以自然而然想到，损失尽可能的小，即让数据在降维后的点尽可能的分散，这比较好理解，比如原来三维空间中的点彼此分离，如果将它们投影到二维空间那么必定会有一些点相互重叠，造成数据的损失，而我们要做的即让经过P矩阵处理后，点尽可能的分散。大家注意，这里的点指的是样本（也就是每列，在这里强调是要和后面区分），我们采用的方法是用方差（刻画样本分散程度的统计量）来进行刻画：
$$\mathrm{var}(x_j)=\frac{1}{M}\sum _{k=1}^M{(x_{jk}-\mu )}^2$$
对于问题二，主成分之间尽可能独立，我们采用的是协方差来刻画：
$$Cov(x_i,x_j)=E[x_i{x}_j]=\frac{1}{M}\sum _{k=1}^M{x}_{ik}{x}_{jk}$$
在这个公式中x和前面公式表述并不是一个意思，这里x代表我们的主成分，也就是降维后特征，我们还没有进行选择，但是要求特征两两之间尽可能的独立，用来刻画事务的某一方面。

那么如何去构造一个矩阵，能够包括以上的信息呢？首先我们对原始数据矩阵进行0均值处理得到矩阵A（方便计算，对我们来说无所谓，因为我们交给matlab计算，在这里处理方便观察解释公式）：

对矩阵A进行如下构造：

整理后可得结果：

到此为止，聪明的同学可以看到，对角线上的都是方差，除了对角线上的都是协方差，我们的优化目标是期望方差尽可能的大，协方差为零。到这里我们的主成分思想全部讲完了，它的精髓可能还需要你多看几遍才能明白，我们要求的是矩阵P，但是优化目标在后面的构造的矩阵中。

4 MATLAB 求解步骤

上面是我总结的主成分分析的步骤，实际上如果你仅仅是为了降维的话，第4步就已经结束了，接下来的5，6是进一步的扩展，是利用主成分分析进行一个综合评价的过程，它的思想有点类似于熵权法求综合评价值，后面我也会讨论这个专题。

最重要的环节！上代码！！

x=[1.0000   40.4000   24.7000    7.2000    6.1000    8.3000    8.7000    2.4420   20.0000;
    2.0000   25.0000   12.7000   11.2000   11.0000   12.9000   20.2000    3.5420    9.1000;
    3.0000   13.2000    3.3000    3.9000    4.3000    4.4000    5.5000    0.5780    3.6000;
    4.0000   22.3000    6.7000    5.6000    3.7000    6.0000    7.4000    0.1760    7.3000;
    5.0000   34.3000   11.8000    7.1000    7.1000    8.0000    8.9000    1.7260   27.5000;
    6.0000   35.6000   12.5000   16.4000   16.7000   22.8000   29.3000    3.0170   26.6000;
    7.0000   22.0000    7.8000    9.9000   10.2000   12.6000   17.6000    0.8470   10.6000;
    8.0000   48.4000   13.4000   10.9000    9.9000   10.9000   13.9000    1.7720   17.8000;
    9.0000   40.6000   19.1000   19.8000   19.0000   29.7000   39.6000    2.4490   35.8000;
   10.0000   24.8000    8.0000    9.8000    8.9000   11.9000   16.2000    0.7890   13.7000;
   11.0000   12.5000    9.7000    4.2000    4.2000    4.6000    6.5000    0.8740    3.9000;
   12.0000    1.8000    0.6000    0.7000    0.7000    0.8000    1.1000    0.0560    1.0000;
   13.0000   32.3000   13.9000    9.4000    8.3000    9.8000   13.3000    2.1260   17.1000;
   14.0000   38.5000    9.1000   11.3000    9.5000   12.2000   16.4000    1.3270   11.6000;
   15.0000   26.2000   10.1000    5.6000   15.6000    7.7000   30.1000    0.1260   25.9000];
g=zscore(x);
r=corrcoef(g);
[x,y,z]=pcacov(r);
num=3;%这个可以进行修改，根据计算得到的Z值进行判断
df=g*x(:,1:num);%行代表公司，所以全部选择
tf=df*z(1:num)/100;
[stf,ind]=sort(tf,'descend'); %stf是综合得分，ind是公司标签
stf=stf', ind=ind'

还剩一个小尾巴没给大家说，大家尝试运行代码后发现得到的主成分个数和你的变量数是一样的，在开头的时候给大家埋下了伏笔，怎么去选择主成分呢？

在pcacov这个函数返回三个值，x每一列代表每一个主成分PC与原始变量之间的系数，y返回的是特征值，z返回的是贡献率（其实是根据某一个特征值在所有特征值的占比得到的，CCA算法也有涉及），主成分都是按照贡献率大小排序的，我们选择累计贡献率大于85%的主成分即可，毕竟我们降维不可能避免数据损失，按照上面的例题，我们三个主成分贡献率达89%，也就说能够解释原始数据的89%，剩下的数据在降维过程丢失了（其实是人为舍去了）。

奈何本人实力有限，必有疏漏，欢迎批评指正！后续推出各种算法，欢迎大家订阅！