主成分分析（PCA）原理总结与实例

微信公众号：幼儿园的学霸

最近在点云处理中,需要获取目标点云的最小包围盒(OBB),在网上看到很多利用PCA求解包围盒的代码,但是代码大多都比较简洁,属于PCA的一个应用,而没有原理的具体描述.因此,抽时间对PCA进行了一些梳理,以方便自己理解.特此记录.

目录

介绍

主成分分析（Principal components analysis，PCA）是最重要的降维方法之一.PCA的主要思想是寻找到数据的主轴方向，由主轴构成一个新的坐标系，新坐标系的维数可以比原维数低，然后数据由原坐标系向新的坐标系投影，这个投影的过程就可以是降维的过程。

PCA将一系列可能相关联的高维变量减少为一系列被称为主成分的低维线性不相关变量。这些低维数据会尽可能地保留原始数据的方差。比如我们有如下的数据分布：

图中数据从原点到右上角呈现散点分布，我们可以通过x轴和y轴来描述整个数据分布的构型。但是如果我们只想使用一个描述符来描述整个数据怎么办？可以预想这个描述符能够最大化解释整个数据集。比如下图：

在图中，我们发现沿着红色虚线的方向能够最大程度将数据分开，或者说在该方向上数据集的方差达到最大；而如果考虑绿色虚线的方向，此时就很难将数据分开。实际上，在这里红色虚线是第一个主成分，绿色虚线是第二个主成分，它们俩必须正交。如果要选择一个主成分那必须是红色虚线的方向。

所以主成分分析事实上就是重新依次寻找方差最大的特征向量作为一个主成分，每个后面的主成分会保留剩余方差的最大值，唯一的限制是它必须和其他的主成分正交。

再比如描述一个人信息时会用体重、身高、发型、爱好、收入、职业等信息，有时根据一个人的体重、身高、发型基本可以确定其性别，从众多属性中选取一两个，而无需其他属性作为参考就确定了一个分类，PCA就是这样一个处理数据常用手段，即利用较少的属性对一组数据分类，PCA是一个降维的数据处理手段,现实中的数据在计算机中可以用一个m行n列矩阵表示，n列可以是体重、身高、发型、爱好、收入、职业等属性信息，每一行就代表一个具体人。利用PCA降维的过程是利用线性代数处理数据的过程,先介绍与PCA有关的线性代数知识。

向量的表示与基变换

内积与投影

两个向量的内积(点积,dot product)被定义为:
$$\vec{A}\cdot \vec{B}=(a_1,a_2,\dots ,a_n{)}^T\cdot (b_1,b_2,\dots ,b_n{)}^T=a_1{b}_1+a_2{b}_2+\cdots +a_n{b}_n$$
内积的结果为一个数值.其几何意义呢?
为方便理解,假设两个二维向量$\vec{A}=(x_1,y_1),\vec{B}=(x_2,y_2)$,在二维平面上的表示如下:

在二维空间中,将向量的内积表示为我们熟悉的形式:
$$\vec{A}\cdot \vec{B}=|\vec{A}||\vec{B}|cos\alpha$$
即:A与B的内积等于A到B的投影长度乘以B的模.再进一步,假设B的模为1,即$|\vec{B}|=1$,则:
$$\vec{A}\cdot \vec{B}=|\vec{A}|cos\alpha$$
也就是说,设向量B的模为1,则A与B的内积值等于向量A向向量B所在直线投影的矢量长度.这个概念后面将会反复用到!

注意,此处用了矢量长度,以和标量长度进行区分.矢量长度可能为负值.

基/基向量

下面继续在二维空间内讨论向量.上文说过,一个二维向量可以对应二维笛卡尔直角坐标系中从原点出发的一个有向线段.例如下面这个向量:

在代数表示方面，我们经常用线段终点的点坐标表示向量，例如上面的向量可以表示为(3,2)，这是我们再熟悉不过的向量表示。
不过我们常常忽略，只有一个(3,2)本身是不能够精确表示一个向量的。我们仔细看一下，这里的3实际表示的是向量在x轴上的投影值是3，在y轴上的投影值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。那么一个向量(3,2)实际是说在x轴投影为3而y轴的投影为2。注意投影是一个矢量，所以可以为负。
更正式的说，向量(x,y)实际上表示线性组合：
$$x(1,0{)}^T+y(0,1{)}^T$$
不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基(or基向量)。

所以，要准确描述向量，首先要确定一组基(基向量)，然后给出在基(基向量)所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。

我们之所以默认选择(1,0)和(0,1)为基，当然是比较方便，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。

在线性代数里，矢量空间的一组元素中，若没有矢量可用有限个其他矢量的线性组合所表示，则称为线性无关或线性独立(linearly independent)，反之称为线性相关(linearly dependent)。
线性相关/线性相关的具体解释可以搜索资料进一步查看.

例如，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}),(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$

现在，我们想获得(3,2)在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，我们只要分别计算(3,2)和两个基的内积，不难得到新的坐标为$(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}})$。下图给出了新的基以及(3,2)在新基上坐标值的示意图：

另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。

基变换的矩阵表示

下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换：

太漂亮了！其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标。可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个2行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：
$$\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ 1\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{ccc}1& 2& 3\\ 1& 2& 3\end{array}\right)=\left(\begin{array}{ccc}\frac{2}{\sqrt{2}}& \frac{4}{\sqrt{2}}& \frac{6}{\sqrt{2}}\\ 0& 0& 0\end{array}\right)$$
于是一组向量的基变换被干净的表示为矩阵的相乘。

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果.

数学表示为:

$$\left(\begin{array}{c}{p}_1\\ p_2\\ ⋮\\ p_R\end{array}\right)\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_M\end{array}\right)=\left(\begin{array}{cccc}{p}_1{a}_1& p_1{a}_2& \cdots & p_1{a}_M\\ p_2{a}_1& p_2{a}_2& \cdots & p_2{a}_M\\ ⋮& ⋮& \ddots & ⋮\\ p_R{a}_1& p_R{a}_2& \cdots & p_R{a}_M\end{array}\right)$$
其中pi是一个行向量，表示第i个基，aj是一个列向量，表示第j个原始数据记录。

特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换。

最后，上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说，一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪，但是如果明白了矩阵相乘的物理意义，其合理性就一目了然了。

协方差矩阵及优化目标

上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？

要完全数学化这个问题非常繁杂，这里我们用一种非形式化的直观方法来看这个问题。

为了避免过于抽象的讨论，我们仍以一个具体的例子展开。假设我们的数据由五条记录组成，将它们表示成矩阵形式：

$$\left(\begin{array}{ccccc}1& 1& 2& 4& 2\\ 1& 3& 3& 4& 4\end{array}\right)$$
其中每一列为一条数据记录，而一行为一个字段。为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0（这样做的道理和好处后面会看到）。

我们看上面的数据，第一个字段均值为2，第二个字段均值为3，所以变换后：

$$\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)$$
我们可以看下五条数据在平面直角坐标系内的样子：

现在问题来了：如果我们必须使用一维来表示这些数据，又希望尽量保留原始的信息，你要如何选择？

通过上一节对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

以上图为例，可以看出如果向x轴投影，那么最左边的两个点会重叠在一起，中间的两个点也会重叠在一起，于是本身四个各不相同的二维点投影后只剩下两个不同的值了，这是一种严重的信息丢失，同理，如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。

下面，我们用数学方法表述这个问题。

方差

上文说到，我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：

$$Var(a)=\frac{1}{m}\sum _{i=1}^m(a_i-\mu {)}^2$$
由于上面我们已经将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示：
$$Var(a)=\frac{1}{m}\sum _{i=1}^m{a}_i^2$$
于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

协方差

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

$$Cov(a,b)=\frac{1}{m}\sum _{i=1}^m{a}_i{b}_i$$
可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。

当协方差为0时，表示两个字段完全独立表示两个字段不相关。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。

至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）。

协方差矩阵

上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。

我们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。于是我们来了灵感：

假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：

$$X=\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_m\\ b_1& b_2& \cdots & b_m\end{array}\right)$$

然后我们用X乘以X的转置，并乘上系数1/m：

$$\frac{1}{m}X{X}^{\mathsf{T}}=\left(\begin{array}{cc}\frac{1}{m}{\sum }_{i=1}^m{a}_i^2& \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i\\ \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i& \frac{1}{m}{\sum }_{i=1}^m{b}_i^2\end{array}\right)$$
奇迹出现了！这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。

根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：
设我们有m个n维数据记录，将其按列排成nxm的矩阵X，设$C=\frac{1}{m}X{X}^{\mathsf{T}}$，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差。

协方差矩阵对角化

根据上述推导，我们发现要达到优化目前，等价于将协方差矩阵对角化：即除对角线外的其它元素化为0，并且在对角线上将元素按大小从上到下排列，这样我们就达到了优化目的。这样说可能还不是很明晰，我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系：

设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵，设Y=PX，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：

$$\begin{array}{lll}D& =& \frac{1}{m}Y{Y}^{\mathsf{T}}\\ & =& \frac{1}{m}(PX)(PX{)}^{\mathsf{T}}\\ & =& \frac{1}{m}PX{X}^{\mathsf{T}}{P}^{\mathsf{T}}\\ & =& P(\frac{1}{m}X{X}^{\mathsf{T}})P^{\mathsf{T}}\\ & =& PC{P}^{\mathsf{T}}\end{array}$$

现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足$PC{P}^{\mathsf{T}}$是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。
由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：
1）实对称矩阵不同特征值对应的特征向量必然正交,其点乘为0。
2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化.
3）数据点在特征向量上投影的方差，为对应的特征值，选择特征值大的特征向量，就是选择点投影方差大的方向，即是具有高信息量的主成分；次佳投影方向位于最佳投影方向的正交空间，是第二大特征值对应的特征向量，以此类推.

特征值和特征向量的相关内容建议进行额外了解

由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为$e_1,e_2,\dots ,e_n$，我们将其按列组成矩阵：
$$E=\left(\begin{array}{cccc}{e}_1& e_2& \cdots & e_n\end{array}\right)$$
则对协方差矩阵C有如下结论：
$$E^{\mathsf{T}}CE=\Lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$$
其中Λ为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。
以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。
到这里，我们发现我们已经找到了需要的矩阵P：
$$P=E^{\mathsf{T}}$$
P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是C的一个特征向量。如果设P按照Λ中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。

至此我们完成了整个PCA的数学原理讨论。

PCA算法步骤

算法步骤

设有m条n维数据,进行PCA算法降维度的过程如下:
1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

The aim of this step is to standardize the range of the continuous initial variables so that each one of them contributes equally to the analysis.

3）求出协方差矩阵$C=\frac{1}{m}X{X}^T$

也有这样的求法:$C=\frac{1}{m-1}X{X}^T$,这可能是标准的协方差矩阵，但这对最后找到的特征向量没有影响，对特征值之间的大小关系也没有影响
m或者m-1 对应的分别对应方差的有偏估计和无偏估计,建议查阅相关资料

4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）$Y=PX$即为降维到k维后的数据

数学实例

数学过程

这里以上文提到的
$$X=\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)$$
为例，我们用PCA方法将这组二维数据其降到一维。

因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：
$$C=\frac{1}{m}X{X}^T=\frac{1}{5}\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)\left(\begin{array}{cc}-1& -2\\ -1& 0\\ 0& 0\\ 2& 1\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}\frac{6}{5}& \frac{4}{5}\\ \frac{4}{5}& \frac{6}{5}\end{array}\right)$$
然后求其特征值和特征向量.
(1)根据矩阵C的特征方程,求特征值:
$$\left|\begin{array}{c}C-\lambda E\end{array}\right|=\left|\begin{array}{cc}\frac{6}{5}-\lambda & \frac{4}{5}\\ \frac{4}{5}& \frac{6}{5}-\lambda \end{array}\right|=(\frac{2}{5}-\lambda )(2-\lambda )=0$$
因此特征值为:
$${\lambda }_1=2,{\lambda }_2=2/5$$
2.求解特征向量.
将每个特征值λ代入线性方程组$(C-\lambda E)x=0$,求基础解系.
2.1 当λ=2时,解解线性方程组:$(C-2E)x=0$
$$\begin{gathered} (C-2E)x=\left(\begin{array}{cc}-\frac{4}{5}& \frac{4}{5}\\ \frac{4}{5}& -\frac{4}{5}\end{array}\right)x⟶ \\ \left(\begin{array}{cc}-1& 1\\ 0& 0\end{array}\right)x=0⟶ \\ \{\begin{array}{rl}-x_1+x_2& =0\\ 0\ast x_1+0\ast x_2& =0\end{array}⟶ \\ \{\begin{array}{rl}{x}_1& =1\\ x_2& =1\end{array}⟶得基础解系 \\ {P}_1=\left(\begin{array}{c}1\\ 1\end{array}\right) \end{gathered}$$
同理,
$$P_2=\left(\begin{array}{c}-1\\ 1\end{array}\right)$$

因此,其对应的特征向量分别是：
$$c_1\left(\begin{array}{c}1\\ 1\end{array}\right),c_2\left(\begin{array}{c}-1\\ 1\end{array}\right)$$
其中对应的特征向量分别是一个通解，c1和c2可取任意实数。那么标准化后的特征向量为：
$$\left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right),\left(\begin{array}{c}-1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)$$
因此我们的矩阵P是：
$$P=E^T={\left(\begin{array}{cc}{P}_1& P_2\end{array}\right)}^T={\left(\begin{array}{cc}1/\sqrt{2}& -1/\sqrt{2}\\ 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)}^T=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ -1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)$$
可以验证协方差矩阵C的对角化：
$$PC{P}^{\mathsf{T}}=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ -1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{cc}6/5& 4/5\\ 4/5& 6/5\end{array}\right)\left(\begin{array}{cc}1/\sqrt{2}& -1/\sqrt{2}\\ 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)=\left(\begin{array}{cc}2& 0\\ 0& 2/5\end{array}\right)$$
最后我们用P的第一行乘以数据矩阵，就得到了降维后的表示：
$$Y=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)=\left(\begin{array}{ccccc}-3/\sqrt{2}& -1/\sqrt{2}& 0& 3/\sqrt{2}& 1/\sqrt{2}\end{array}\right)$$
降维投影结果如下图：

python实现

上面过程的求解可以利用python实现,代码如下:

##Python实现PCA
import numpy as np
def pca(X, k): # k is the components you want
    """

    :param X: 输入数据 shape=[N,M].N:每一行表示一个维度.M:每一列表示一个数据
    :param k: 保留的主成分数量/将输入数据X降为k维
    :return:
    """

    # 求每一维度的均值
    n_features, n_samples = X.shape
    mean = np.mean(X, axis=-1).reshape(n_features,1)
    # 去中心化.实际上是否去中心化对求到的协方差矩阵并无影响,只是方便后面进行降维
    norm_X = X - mean

    # 协方差矩阵(Covariance matrix).使用矩阵乘法的形式
    C = np.dot(norm_X,np.transpose(norm_X))/(n_samples-1)
    # 协方差矩阵.使用numpy函数
    # C = np.cov(norm_X,rowvar=True)#rowvar=True时,每一行为特征,而每一列表示一条记录. rowvar=False相反,表示数据的每一列代表一个feature
    # print(C)

    # 计算特征值和特征向量 每个特征值对应的是特征矩阵的每个列向量
    eig_val, eig_vec = np.linalg.eig(C)
    # 求前k大的特征向量,组成特征矩阵P
    sorted_indices = np.argsort(-eig_val)#将特征值按从大到小排序(通过给特征值取反),index保留的是对应原eig_val中的下标
    E = eig_vec[:,sorted_indices]
    P = np.transpose(E)[:k]
    print(P)

    # # 求取数据在新基下的结果
    # # Y=PX
    # Y = np.dot(P,X)
    # or,网络上常见用法为此种.
    # 但存疑:因为我们要求解的是输入数据在新基下的结果,而不是去中心化后的数据在新基下的结果
    Y = np.dot(P,norm_X)

    return Y

X = np.array([[-1,-2],[-1,0],[0,0],[2,1],[0,1]]).transpose()
# X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]]).transpose()
print(pca(X, 1))

运行结果如下:

#自定义实现
[[-2.12132034 -0.70710678  0.          2.12132034  0.70710678]]

#sklearn库实现
[[ 2.12132034]
 [ 0.70710678]
 [-0.        ]
 [-2.12132034]
 [-0.70710678]]

自定义实现结果和上面的数学求解结果一致.但是和sklearn库的结果在数值上正好相反,为什么呢?感兴趣可以思考下.

我的理解是在求解转换基P的过程中,在对特征向量进行归一化时,选择的方向不同,导致基的方向也不同,因此结果正好相反.如果将上面代码中的P = np.transpose(E)[:k]修改为P = -np.transpose(E)[:k] ,可以在观察下结果有何不同.

应用实例

(略)
可以查看OpenCV和PCL中关于PCA的应用

PCA算法总结

PCA算法的主要优点有：
1）仅仅需要以方差衡量信息量，不受数据集以外的因素影响。　
2）各主成分之间正交，可消除原始数据成分间的相互影响的因素。
3）计算方法简单，主要运算是特征值分解，易于实现。

PCA算法的主要缺点有：
1）主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。
2）方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。

一点点扩展:对二维点集进行 1)最小二乘拟合直线,得到直线的方向向量;2)进行PCA求解主方向.那么拟合的直线方向向量和求解的主方向是否一致呢? 感兴趣的可以试下,基于上面的程序,基本3行代码即可验证结果

参考资料

1.PCA的数学原理
2.Principal components analysis (PCA).
3.Principal Component Analysis (PCA) : An Overview

---

下面的是我的公众号二维码图片，按需关注。