【五万字总结PCA】【李航|PRML】统计学习方法--16.主成分分析（详细推导）
作者：互联网
本文略微有点长，请大家耐心观看，你一定会有收获
文章目录
PCA数学原理
首先，给出不太严谨的数学推导PCA
PCA（Principal Component Analysis）是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。网上关于PCA的文章有很多，但是大多数只描述了PCA的分析过程，而没有讲述其中的原理。这篇文章的目的是介绍PCA的基本数学原理，帮助读者了解PCA的工作机制是什么。
当然我并不打算把文章写成纯数学文章，而是希望用直观和易懂的方式叙述PCA的数学原理，所以整个文章不会引入严格的数学推导。希望读者在看完这篇文章后能更好的明白PCA的工作原理。
数据的向量表示及降维问题
一般情况下，在数据挖掘和机器学习中，数据被表示为向量。例如某个淘宝店2012年全年的流量及交易情况可以看成一组记录的集合，其中每一天的数据是一条记录，格式如下：
(日期, 浏览量, 访客数, 下单数, 成交数, 成交金额)
其中“日期”是一个记录标志而非度量值，而数据挖掘关心的大多是度量值，因此如果我们忽略日期这个字段后，我们得到一组记录，每条记录可以被表示为一个五维向量，其中一条看起来大约是这个样子：
( 500 , 240 , 25 , 13 , 2312.15 ) ⊤ (500,240,25,13,2312.15)^{\top} (500,240,25,13,2312.15)⊤
注意这里我用了转置，因为习惯上使用列向量表示一条记录（后面会看到原因），本文后面也会遵循这个准则。不过为了方便有时我会省略转置符号，但我们说到向量默认都是指列向量。
我们当然可以对这一组五维向量进行分析和挖掘，不过我们知道，很多机器学习算法的复杂度和数据的维数有着密切关系，甚至与维数呈指数级关联。当然，这里区区五维的数据，也许还无所谓，但是实际机器学习中处理成千上万甚至几十万维的情况也并不罕见，在这种情况下，机器学习的资源消耗是不可接受的，因此我们必须对数据进行降维。
降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。
举个例子，假如某学籍数据有两列M和F，其中M列的取值是如何此学生为男性取值1，为女性取值0；而F列是学生为女性取值1，男性取值0。此时如果我们统计全部学籍数据，会发现对于任何一条记录来说，当M为1时F必定为0，反之当M为0时F必定为1。在这种情况下，我们将M或F去掉实际上没有任何信息的损失，因为只要保留一列就可以完全还原另一列。
当然上面是一个极端的情况，在现实中也许不会出现，不过类似的情况还是很常见的。例如上面淘宝店铺的数据，从经验我们可以知道，“浏览量”和“访客数”往往具有较强的相关关系，而“下单数”和“成交数”也具有较强的相关关系。这里我们非正式的使用“相关关系”这个词，可以直观理解为“当某一天这个店铺的浏览量较高（或较低）时，我们应该很大程度上认为这天的访客数也较高（或较低）”。后面的章节中我们会给出相关性的严格数学定义。
这种情况表明，如果我们删除浏览量或访客数其中一个指标，我们应该期待并不会丢失太多信息。因此我们可以删除一个，以降低机器学习算法的复杂度。
上面给出的是降维的朴素思想描述，可以有助于直观理解降维的动机和可行性，但并不具有操作指导意义。例如，我们到底删除哪一列损失的信息才最小？亦或根本不是单纯删除几列，而是通过某些变换将原始数据变为更少的列但又使得丢失的信息最小？到底如何度量丢失信息的多少？如何根据原始数据决定具体的降维操作步骤？
要回答上面的问题，就要对降维问题进行数学化和形式化的讨论。而PCA是一种具有严格数学基础并且已被广泛采用的降维方法。下面我不会直接描述PCA，而是通过逐步分析问题，让我们一起重新“发明”一遍PCA。
向量的表示及基变换
既然我们面对的数据被抽象为一组向量，那么下面有必要研究一些向量的数学性质。而这些数学性质将成为后续导出PCA的理论基础。
内积与投影
下面先来看一个高中就学过的向量运算：内积。两个维数相同的向量的内积被定义为：
( a 1 , a 2 , ⋯ , a n ) ⊤ ⋅ ( b 1 , b 2 , ⋯ , b n ) ⊤ = a 1 b 1 + a 2 b 2 + ⋯ + a n b n \left(a_{1}, a_{2}, \cdots, a_{n}\right)^{\top} \cdot\left(b_{1}, b_{2}, \cdots, b_{n}\right)^{\top}=a_{1} b_{1}+a_{2} b_{2}+\cdots+a_{n} b_{n} (a1,a2,⋯,an)⊤⋅(b1,b2,⋯,bn)⊤=a1b1+a2b2+⋯+anbn
内积运算将两个向量映射为一个实数。其计算方式非常容易理解，但是其意义并不明显。下面我们分析内积的几何意义。假设A和B是两个n维向量，我们知道n维向量可以等价表示为n维空间中的一条从原点发射的有向线段，为了简单起见我们假设A和B均为二维向量，则A=(x1,y1)，B=(x2,y2)。则在二维平面上A和B可以用两条发自原点的有向线段表示，见下图：
好，现在我们从A点向B所在直线引一条垂线。我们知道垂线与B的交点叫做A在B上的投影，再设A与B的夹角是a，则投影的矢量长度为 ∣ A ∣ cos ( a ) |A| \cos (a) ∣A∣cos(a) ，其中 ∣ A ∣ = x 1 2 + y 1 2 |A|=\sqrt{x_{1}^{2}+y_{1}^{2}} ∣A∣=x12+y12 的是向量A的模，也就是A线段的标量长度。
注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于0，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。
到这里还是看不出内积和这东西有什么关系，不过如果我们将内积表示为另一种我们熟悉的形式：
A ⋅ B = ∣ A ∣ ∣ B ∣ cos ( a ) A \cdot B=|A||B| \cos (a) A⋅B=∣A∣∣B∣cos(a)
现在事情似乎是有点眉目了：A与B的内积等于A到B的投影长度乘以B的模。再进一步，如果我们假设B的模为1，即让|B|=1，那么就变成了：
A ⋅ B = ∣ A ∣ cos ( a ) A \cdot B=|A| \cos (a) A⋅B=∣A∣cos(a)
也就是说， 设向量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 ) ⊤ + y ( 0 , 1 ) ⊤ x(1,0)^{\top}+y(0,1)^{\top} x(1,0)⊤+y(0,1)⊤
不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基。
所以， 要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了 。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。
我们之所以默认选择(1,0)和(0,1)为基，当然是比较方便，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。
例如，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为 ( 1 2 , 1 2 ) \left(\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right) (2 1,2 1) 和 ( − 1 2 , 1 2 ) \left(-\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right) (−2 1,2 1) 。
现在，我们想获得(3,2)在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，我们只要分别计算(3,2)和两个基的内积，不难得到新的坐标为 ( 5 2 , − 1 2 ) \left(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}}\right) (2 5,−2 1)。下图给出了新的基以及(3,2)在新基上坐标值的示意图：
另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。
基变换的矩阵表示
下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换：
( 1 / 2 1 / 2 − 1 / 2 1 / 2 ) ( 3 2 ) = ( 5 / 2 − 1 / 2 ) \left(\begin{array}{cc} 1 / \sqrt{2} & 1 / \sqrt{2} \\ -1 / \sqrt{2} & 1 / \sqrt{2} \end{array}\right)\left(\begin{array}{l} 3 \\ 2 \end{array}\right)=\left(\begin{array}{c} 5 / \sqrt{2} \\ -1 / \sqrt{2} \end{array}\right) (1/2 −1/2 1/2 1/2 )(32)=(5/2 −1/2 )
太漂亮了！其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标。可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：
( 1 / 2 1 / 2 − 1 / 2 1 / 2 ) ( 1 2 3 1 2 3 ) = ( 2 / 2 4 / 2 6 / 2 0 0 0 ) \left(\begin{array}{cc} 1 / \sqrt{2} & 1 / \sqrt{2} \\ -1 / \sqrt{2} & 1 / \sqrt{2} \end{array}\right)\left(\begin{array}{lll} 1 & 2 & 3 \\ 1 & 2 & 3 \end{array}\right)=\left(\begin{array}{ccc} 2 / \sqrt{2} & 4 / \sqrt{2} & 6 / \sqrt{2} \\ 0 & 0 & 0 \end{array}\right) (1/2 −1/2 1/2 1/2 )(112233)=(2/2 04/2 06/2 0)
于是一组向量的基变换被干净的表示为矩阵的相乘。
一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果 。
数学表示为：
( p 1 p 2 ⋮ p R ) ( a 1 a 2 ⋯ a M ) = ( p 1 a 1 p 1 a 2 ⋯ p 1 a M p 2 a 1 p 2 a 2 ⋯ p 2 a M ⋮ ⋮ ⋱ ⋮ p R a 1 p R a 2 ⋯ p R a M ) \left(\begin{array}{c} p_{1} \\ p_{2} \\ \vdots \\ p_{R} \end{array}\right)\left(\begin{array}{llll} 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} \\ \vdots & \vdots & \ddots & \vdots \\ p_{R} a_{1} & p_{R} a_{2} & \cdots & p_{R} a_{M} \end{array}\right) ⎝⎜⎜⎜⎛p1p2⋮pR⎠⎟⎟⎟⎞(a1a2⋯aM)=⎝⎜⎜⎜⎛p1a1p2a1⋮pRa1p1a2p2a2⋮pRa2⋯⋯⋱⋯p1aMp2aM⋮pRaM⎠⎟⎟⎟⎞
其中 p i p_i pi是一个行向量，表示第i个基， a j a_j aj是一个列向量，表示第j个原始数据记录。
特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换。
最后，上述分析同时给矩阵相乘找到了一种物理解释： 两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去 。更抽象的说，一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪，但是如果明白了矩阵相乘的物理意义，其合理性就一目了然了。
协方差矩阵及优化目标
上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？
要完全数学化这个问题非常繁杂，这里我们用一种非形式化的直观方法来看这个问题。
为了避免过于抽象的讨论，我们仍以一个具体的例子展开。假设我们的数据由五条记录组成，将它们表示成矩阵形式：
( 1 1 2 4 2 1 3 3 4 4 ) \left(\begin{array}{lllll} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end{array}\right) (1113234424)
其中每一列为一条数据记录，而一行为一个字段。为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0（这样做的道理和好处后面会看到）。
我们看上面的数据，第一个字段均值为2，第二个字段均值为3，所以变换后：
( − 1 − 1 0 2 0 − 2 0 0 1 1 ) \left(\begin{array}{ccccc} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{array}\right) (−1−2−10002101)
我们可以看下五条数据在平面直角坐标系内的样子：
现在问题来了：如果我们必须使用一维来表示这些数据，又希望尽量保留原始的信息，你要如何选择？
通过上一节对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。
那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。
以上图为例，可以看出如果向x轴投影，那么最左边的两个点会重叠在一起，中间的两个点也会重叠在一起，于是本身四个各不相同的二维点投影后只剩下两个不同的值了，这是一种严重的信息丢失，同理，如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。
下面，我们用数学方法表述这个问题。
方差
上文说到，我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：
Var ( a ) = 1 m ∑ i = 1 m ( a i − μ ) 2 \operatorname{Var}(a)=\frac{1}{m} \sum_{i=1}^{m}\left(a_{i}-\mu\right)^{2} Var(a)=m1i=1∑m(ai−μ)2
由于上面我们已经将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示：
Var ( a ) = 1 m ∑ i = 1 m a i 2 \operatorname{Var}(a)=\frac{1}{m} \sum_{i=1}^{m} a_{i}^{2} Var(a)=m1i=1∑mai2
于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。
协方差
对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。
如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。
数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：
Cov ( a , b ) = 1 m ∑ i = 1 m a i b i \operatorname{Cov}(a, b)=\frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} Cov(a,b)=m1i=1∑maibi
可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。
当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。
至此，我们得到了降维问题的优化目标： 将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差） 。
协方差矩阵
上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。
我们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。于是我们来了灵感：
假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：
X = ( a 1 a 2 ⋯ a m b 1 b 2 ⋯ b m ) X=\left(\begin{array}{llll} a_{1} & a_{2} & \cdots & a_{m} \\ b_{1} & b_{2} & \cdots & b_{m} \end{array}\right) X=(a1b1a2b2⋯⋯ambm)
然后我们用X乘以X的转置，并乘上系数1/m：
1 m X X ⊤ = ( 1 m ∑ i = 1 m a i 2 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m b i 2 ) \frac{1}{m} X X^{\top}=\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) m1XX⊤=(m1∑i=1mai2m1∑i=1maibim1∑i=1maibim1∑i=1mbi2)
奇迹出现了！这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。
根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：
设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设 C = 1 m X X T C=\frac{1}{m} X X^{\mathrm{T}} C=m1XXT，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差 。
协方差矩阵对角化
根据上述推导，我们发现要达到优化目前，等价于将协方差矩阵对角化：即除对角线外的其它元素化为0，并且在对角线上将元素按大小从上到下排列，这样我们就达到了优化目的。这样说可能还不是很明晰，我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系：
设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵，设Y=PX，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：
D = 1 m Y Y ⊤ = 1 m ( P X ) ( P X ) ⊤ = 1 m P X X ⊤ P ⊤ = P ( 1 m X X ⊤ ) P ⊤ = P C P ⊤ \begin{aligned} D &=\frac{1}{m} Y Y^{\top} \\ &=\frac{1}{m}(P X)(P X)^{\top} \\ &=\frac{1}{m} P X X^{\top} P^{\top} \\ &=P\left(\frac{1}{m} X X^{\top}\right) P^{\top} \\ &=P C P^{\top} \end{aligned} D=m1YY⊤=m1(PX)(PX)⊤=m1PXX⊤P⊤=P(m1XX⊤)P⊤=PCP⊤
现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了 寻找一个矩阵P，满足 P C P ⊤ P C P^{\top} PCP⊤是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件 。
至此，我们离“发明”PCA还有仅一步之遥！
现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。
由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：
1）实对称矩阵不同特征值对应的特征向量必然正交。
2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化。
由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为 e 1 , e 2 , ⋯ , e n e_{1}, e_{2}, \cdots, e_{n} e1,e2,⋯,en，我们将其按列组成矩阵：
E = ( e 1 e 2 ⋯ e n ) E=\left(\begin{array}{llll} e_{1} & e_{2} & \cdots & e_{n} \end{array}\right) E=(e1e2⋯en)
则对协方差矩阵C有如下结论：
E ⊤ C E = Λ = ( λ 1 λ 2 ⋱ λ n ) E^{\top} C E=\Lambda=\left(\begin{array}{llll} \lambda_{1} & & & \\ & \lambda_{2} & & \\ & & \ddots & \\ & & & \lambda_{n} \end{array}\right) E⊤CE=Λ=⎝⎜⎜⎛λ1λ2⋱λn⎠⎟⎟⎞
其中Λ为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。
以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。
到这里，我们发现我们已经找到了需要的矩阵P：
P = E T P=E^{\mathrm{T}} P=ET
P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是C的一个特征向量。如果设P按照Λ中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。
至此我们完成了整个PCA的数学原理讨论。在下面的一节，我们将给出PCA的一个实例。
算法
为了巩固上面的理论，我们在这一节给出一个具体的PCA实例。
总结一下PCA的算法步骤：
设有m条n维数据。
1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵 = C = 1 m X X T =C=\frac{1}{m} X X^{\mathrm{T}} =C=m1XXT
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）Y=PX即为降维到k维后的数据
进一步讨论
根据上面对PCA的数学原理的解释，我们可以了解到一些PCA的能力和限制。PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。
因此，PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑Kernel PCA，通过Kernel函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。
最后需要说明的是，PCA是一种无参数技术，也就是说面对同样的数据，如果不考虑清洗，谁来做结果都一样，没有主观参数的介入，所以PCA便于通用实现，但是本身无法个性化的优化。
希望这篇文章能帮助朋友们了解PCA的数学理论基础和实现原理，借此了解PCA的适用场景和限制，从而更好的使用这个算法。
有了上面的基础，那么再来看李航老师书上的推导（其中补充了推导不详细的部分），他这个新的部分在于SVD来求特征向量。
主成分分析（李航）
引出降维，二维三维上的结论在高维就不work了
- 主成分分析(principal component analysis, PCA)是一种常用的无监督学习方法,这一方法利用正交变换把由线性相关变量表示的观测数据转换为少数几个由线性无关变量表示的数据,线性无关的变量称为主成分。
- 主成分的个数通常小于原始变量的个数,所以主成分分析属于降维方法。
- 主成分分析主要用于发现数据中的基本结构,即数据中变量之间的关系
16.1 总体主成分分析
16.1.1 基本想法
主成分分析中,首先对给定数据进行规范化,使得数据每一变量的平均值为0,方差为1。之后对数据进行正交变换,原来由线性相关变量表示的数据,通过正交变换变成由若干个线性无关的新变量表示的数据。新变量是可能的正交变换中变量的方差的和(信息保存)最大的,方差表示在新变量上信息的大小。将新变量依次称为第一主成分、第二主成分等。这就是主成分分析的基本思想。
通过主成分分析,可以利用主成分近似地表示原始数据,这可理解为发现数据的“基本结构” ;也可以把数据由少数主成分表示,这可理解为对数据降维。
下面再看方差最大的解释。
假设有两个变量 x 1 x_{1} x1 和 x 2 x_{2} x2, 三个样本点 A 、 B , C A 、 B, C A、B,C, 样 本分布在由 x 1 x_{1} x1 和 x 2 x_{2} x2 轴组成的坐标系中, 如下图所示。对坐标系进行旋转变换, 得到新的坐标轴 y 1 y_{1} y1, 表示新的变量 y 1 y_{1} y1 。样本点 A 、 B 、 C A 、 B 、 C A、B、C 在 y 1 y_{1} y1 轴上投影, 得到 y 1 y_{1} y1 轴的 坐标值 A ′ 、 B ′ 、 C ′ 。 A^{\prime}、 B^{\prime} 、 C^{\prime} 。 A′、B′、C′。 坐标值的平方和 O A ′ 2 + O B ′ 2 + O C ′ 2 O A^{\prime 2}+O B^{\prime 2}+O C^{\prime 2} OA′2+OB′2+OC′2 表示样本在变量 y 1 y_{1} y1 上的方 差和。主成分分析旨在选取正交变换中方差最大的变量, 作为第一主成分, 也就是旋转变换中坐标值的平方和最大的轴。==注意到旋转变换中样本点到原点的距离的平方和 O A 2 + O B 2 + O C 2 O A^{2}+O B^{2}+O C^{2} OA2+OB2+OC2 保持不变, 根据勾股定理, 坐标值的平方和 O A ′ 2 + O B ′ 2 + O C ′ 2 O A^{\prime 2}+O B^{\prime 2}+O C^{\prime 2} OA′2+OB′2+OC′2 最大等价于样本点到 y 1 y_{1} y1 轴的距离的平方和 A A ′ 2 + B B ′ 2 + C C ′ 2 A A^{\prime 2}+B B^{\prime 2}+C C^{\prime 2} AA′2+BB′2+CC′2 最小。==所以, 等价地, 主成分分析在旋转变换中选取离样本点的距离平方和最小的轴, 作为第一主成分。第二主成分等的选取, 在保证与已选坐标轴正交的条件下, 类似地进行。
16.1.2 定义和导出
假设 x = ( x 1 , x 2 , ⋯ , x m ) T \boldsymbol{x}=\left(x_{1}, x_{2}, \cdots, x_{m}\right)^{\mathrm{T}} x=(x1,x2,⋯,xm)T 是 m m m 维随机变量, 其均值向量是 μ \boldsymbol{\mu} μ
μ = E ( x ) = ( μ 1 , μ 2 , ⋯ , μ m ) T \boldsymbol{\mu}=E(\boldsymbol{x})=\left(\mu_{1}, \mu_{2}, \cdots, \mu_{m}\right)^{\mathrm{T}} μ=E(x)=(μ1,μ2,⋯,μm)T
协方差矩阵是 Σ \Sigma Σ
Σ = cov ( x , x ) = E [ ( x − μ ) ( x − μ ) T ] \Sigma=\operatorname{cov}(\boldsymbol{x}, \boldsymbol{x})=E\left[(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\mathrm{T}}\right] Σ=cov(x,x)=E[(x−μ)(x−μ)T]
考虑由 m m m 维随机变量 x x x 到 m m m 维随机变量 y = ( y 1 , y 2 , ⋯ , y m ) T \boldsymbol{y}=\left(y_{1}, y_{2}, \cdots, y_{m}\right)^{\mathrm{T}} y=(y1,y2,⋯,ym)T 的线性变换
y i = α i T x = α 1 i x 1 + α 2 i x 2 + ⋯ + α m i x m y_{i}=\alpha_{i}^{\mathrm{T}} \boldsymbol{x}=\alpha_{1 i} x_{1}+\alpha_{2 i} x_{2}+\cdots+\alpha_{m i} x_{m} yi=αiTx=α1ix1+α2ix2+⋯+αmixm
其中
α
i
T
=
(
α
1
i
,
α
2
i
,
⋯
,
α
m
i
)
,
i
=
1
,
2
,
⋯
,
m
\alpha_{i}^{\mathrm{T}}=\left(\alpha_{1 i}, \alpha_{2 i}, \cdots, \alpha_{m i}\right), i=1,2, \cdots, m
αiT=(α1i,α2i,⋯,αmi),i=1,2,⋯,m 。
由随机变量的性质可知,
E ( y i ) = α i T μ , i = 1 , 2 , ⋯ , m var ( y i ) = α i T Σ α i , i = 1 , 2 , ⋯ , m cov ( y i , y j ) = α i T Σ α j , i = 1 , 2 , ⋯ , m ; j = 1 , 2 , ⋯ , m \begin{aligned} &E\left(y_{i}\right)=\alpha_{i}^{\mathrm{T}} \mu, \quad i=1,2, \cdots, m \\ &\operatorname{var}\left(y_{i}\right)=\alpha_{i}^{\mathrm{T}} \Sigma \alpha_{i}, \quad i=1,2, \cdots, m \\ &\operatorname{cov}\left(y_{i}, y_{j}\right)=\alpha_{i}^{\mathrm{T}} \Sigma \alpha_{j}, \quad i=1,2, \cdots, m ; \quad j=1,2, \cdots, m \end{aligned} E(yi)=αiTμ,i=1,2,⋯,mvar(yi)=αiTΣαi,i=1,2,⋯,mcov(yi,yj)=αiTΣαj,i=1,2,⋯,m;j=1,2,⋯,m
详细推导：
E ( y i ) = E ( α i T x ) = α i T E ( x ) = α i T μ E\left(y_{i}\right)=E(\alpha_{i}^{\mathrm{T}} \boldsymbol{x})=\alpha_{i}^{\mathrm{T}}E(\boldsymbol{x})=\alpha_{i}^{\mathrm{T}} \mu E(yi)=E(αiTx)=αiTE(x)=αiTμ
var ( y i ) = E ( ( y i − E ( y i ) ) 2 ) = E ( ( α i T x − α i T μ ) ) = E ( α i T ( x − μ ) ( α i T ( x − μ ) T ) = E ( α i T ( x − μ ) ( x − μ ) T α i ) = α i T E [ ( x − μ ) ( x − μ ) T ] α i = α i T Σ α i \operatorname{var}\left(y_{i}\right)=E((y_i-E(y_i))^2)=E((\alpha_{i}^{\mathrm{T}} \boldsymbol{x}-\alpha_{i}^{\mathrm{T}} \mu))=E(\alpha_i^T(\boldsymbol{x}-\mu)(\alpha_i^T(\boldsymbol{x}-\mu)^T)=E(\alpha_i^T(\boldsymbol{x}-\mu)(\boldsymbol{x}-\mu)^T\alpha_i)=\alpha_i^TE\left[(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\mathrm{T}}\right]\alpha_i=\alpha_{i}^{\mathrm{T}} \Sigma \alpha_{i} var(yi)=E((yi−E(yi))2)=E((αiTx−αiTμ))=E(αiT(x−μ)(αiT(x−μ)T)=E(αiT(x−μ)(x−μ)Tαi)=αiTE[(x−μ)(x−μ)T]αi=αiTΣαi
cov ( y i , y j ) = E [ [ y i − E ( y i ) ] [ y j − E ( y j ) ] ] = E [ ( y i − α i T μ ) ( y j − α j T μ ) ] = E [ ( α i T x − α i T μ ) ( α j T x − α j T μ ) ] = E [ ( α i T x − α i T μ ) ( ( α j T x − α j T μ ) ) T ] = E [ α i T ( x − μ ) ( x − μ ) T α j ] = α i T Σ α j \begin{aligned}\operatorname{cov}\left(y_{i}, y_{j}\right)&=E[[y_i-E(y_i)][y_j-E(y_j)]]\\&=E[(y_i-\alpha_i^T\mu)(y_j-\alpha_j^T\mu)]\\&=E[(\alpha_i^Tx-\alpha_i^T\mu)(\alpha_j^Tx-\alpha_j^T\mu)]\\&=E[(\alpha_i^Tx-\alpha_i^T\mu)((\alpha_j^Tx-\alpha_j^T\mu))^T]\\&=E[\alpha_i^T(x-\mu)(x-\mu)^T\alpha_j]\\&=\alpha_i^T\Sigma \alpha_{j}\end{aligned} cov(yi,yj)=E[[yi−E(yi)][yj−E(yj)]]=E[(yi−αiTμ)(yj−αjTμ)]=E[(αiTx−αiTμ)(αjTx−αjTμ)]=E[(αiTx−αiTμ)((αjTx−αjTμ))T]=E[αiT(x−μ)(x−μ)Tαj]=αiTΣαj
总体主成分
定义
16.1
16.1
16.1 (总体主成分） 给定一个如式
y
i
=
α
i
T
x
=
α
1
i
x
1
+
α
2
i
x
2
+
⋯
+
α
m
i
x
m
y_{i}=\alpha_{i}^{\mathrm{T}} \boldsymbol{x}=\alpha_{1 i} x_{1}+\alpha_{2 i} x_{2}+\cdots+\alpha_{m i} x_{m}
yi=αiTx=α1ix1+α2ix2+⋯+αmixm所示的线性变换, 如果它们满足下列条件:
(1) 系数向量
α
i
T
\alpha_{i}^{\mathrm{T}}
αiT 是单位向量, 即
α
i
T
α
i
=
1
,
i
=
1
,
2
,
⋯
,
m
\alpha_{i}^{\mathrm{T}} \alpha_{i}=1, i=1,2, \cdots, m
αiTαi=1,i=1,2,⋯,m;
(2) 变量
y
i
y_{i}
yi 与
y
j
y_{j}
yj 互不相关, 即
cov
(
y
i
,
y
j
)
=
0
(
i
≠
j
)
\operatorname{cov}\left(y_{i}, y_{j}\right)=0(i \neq j)
cov(yi,yj)=0(i=j);
(3) 变量
y
1
y_{1}
y1 是
x
x
x 的所有线性变换中方差最大的;
y
2
y_{2}
y2 是与
y
1
y_{1}
y1 不相关的
x
x
x 的所有 线性变换中方差最大的;一般地,
y
i
y_{i}
yi 是与
y
1
,
y
2
,
⋯
,
y
i
−
1
(
i
=
1
,
2
,
⋯
,
m
)
y_{1}, y_{2}, \cdots, y_{i-1}(i=1,2, \cdots, m)
y1,y2,⋯,yi−1(i=1,2,⋯,m) 都不相关 的
x
x
x 的所有线性变换中方差最大的; 这时分别称
y
1
,
y
2
,
⋯
,
y
m
y_{1}, y_{2}, \cdots, y_{m}
y1,y2,⋯,ym 为
x
x
x 的第一主成分、第二主成分、… 第
m
m
m 主成分。
定义中的条件
(
1
)
(1)
(1) 表明线性变换是正交变换,
α
1
,
α
2
,
⋯
,
α
m
\alpha_{1}, \alpha_{2}, \cdots, \alpha_{m}
α1,α2,⋯,αm 是其一组标准正交基,
α i T α j = { 1 , i = j 0 , i ≠ j \alpha_{i}^{\mathrm{T}} \alpha_{j}= \begin{cases}1, & i=j \\ 0, & i \neq j\end{cases} αiTαj={1,0,i=ji=j
条件 ( 2 ) ( 3 ) (2)(3) (2)(3) 给出了一个求主成分的方法: 第一步, 在 x x x 的所有线性变换
α 1 T x = ∑ i = 1 m α i 1 x i \alpha_{1}^{\mathrm{T}} \boldsymbol{x}=\sum_{i=1}^{m} \alpha_{i 1} x_{i} α1Tx=i=1∑mαi1xi
中, 在 α 1 T α 1 = 1 \alpha_{1}^{\mathrm{T}} \alpha_{1}=1 α1Tα1=1 条件下, 求方差最大的, 得到 x \boldsymbol{x} x 的第一主成分; 第二步, 在与 α 1 T x \alpha_{1}^{\mathrm{T}} \boldsymbol{x} α1Tx 不相关的 x \boldsymbol{x} x 的所有线性变换
α 2 T x = ∑ i = 1 m α i 2 x i \alpha_{2}^{\mathrm{T}} \boldsymbol{x}=\sum_{i=1}^{m} \alpha_{i 2} x_{i} α2Tx=i=1∑mαi2xi
中, 在 α 2 T α 2 = 1 \alpha_{2}^{\mathrm{T}} \alpha_{2}=1 α2Tα2=1 条件下, 求方差最大的, 得到 x \boldsymbol{x} x 的第二主成分; 第 k k k 步, 在与 α 1 T x , α 2 T x , ⋯ , α k − 1 T x \alpha_{1}^{\mathrm{T}} \boldsymbol{x}, \alpha_{2}^{\mathrm{T}} \boldsymbol{x}, \cdots, \alpha_{k-1}^{\mathrm{T}} \boldsymbol{x} α1Tx,α2Tx,⋯,αk−1Tx 不相关的 x \boldsymbol{x} x 的所有线性变换
α k T x = ∑ i = 1 m α i k x i \alpha_{k}^{\mathrm{T}} \boldsymbol{x}=\sum_{i=1}^{m} \alpha_{i k} x_{i} αkTx=i=1∑mαikxi
中, 在 α k T α k = 1 \alpha_{k}^{\mathrm{T}} \alpha_{k}=1 αkTαk=1 条件下, 求方差最大的, 得到 x \boldsymbol{x} x 的第 k k k 主成分; 如此继续下去, 直到 得到 x x x 的第 m m m 主成分。
16.1.3 主要性质
定理 16.1 16.1 16.1 设 x x x 是 m m m 维随机变量, Σ \Sigma Σ 是 x x x 的协方差矩阵, Σ \Sigma Σ 的特征值分别是 λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ λ m ⩾ 0 \lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{m} \geqslant 0 λ1⩾λ2⩾⋯⩾λm⩾0, 特征值对应的单位特征向量分别是 α 1 , α 2 , ⋯ , α m \alpha_{1}, \alpha_{2}, \cdots, \alpha_{m} α1,α2,⋯,αm, 则 x x x 的 第 k k k 主成分是
y k = α k T x = α 1 k x 1 + α 2 k x 2 + ⋯ + α m k x m , k = 1 , 2 , ⋯ , m y_{k}=\alpha_{k}^{\mathrm{T}} \boldsymbol{x}=\alpha_{1 k} x_{1}+\alpha_{2 k} x_{2}+\cdots+\alpha_{m k} x_{m}, \quad k=1,2, \cdots, m yk=αkTx=α1kx1+α2kx2+⋯+αmkxm,k=1,2,⋯,m
x x x 的第 k k k 主成分的方差是
var ( y k ) = α k T Σ α k = λ k , k = 1 , 2 , ⋯ , m \operatorname{var}\left(y_{k}\right)=\alpha_{k}^{\mathrm{T}} \Sigma \alpha_{k}=\lambda_{k}, \quad k=1,2, \cdots, m var(yk)=αkTΣαk=λk,k=1,2,⋯,m
即协方差矩阵 Σ \Sigma Σ 的第 k k k 个特征值。
证明 采用拉格朗日乘子法求出主成分。
首先求 x x x 的第一主成分 y 1 = α 1 T x y_{1}=\alpha_{1}^{\mathrm{T}} \boldsymbol{x} y1=α1Tx, 即求系数向量 α 1 \alpha_{1} α1 。由定义 16.1 16.1 16.1 知, 第一主成 分的 α 1 \alpha_{1} α1 是在 α 1 T α 1 = 1 \alpha_{1}^{\mathrm{T}} \alpha_{1}=1 α1Tα1=1 条件下, x x x 的所有线性变换中使方差var ( α 1 T x ) = α 1 T Σ α 1 \operatorname{var}\left(\alpha_{1}^{\mathrm{T}} \boldsymbol{x}\right)=\alpha_{1}^{\mathrm{T}} \Sigma \alpha_{1} var(α1Tx)=α1TΣα1
达到最大的。
求第一主成分就是求解约束最优化问题:max α 1 α 1 T Σ α 1 s.t. α 1 T α 1 = 1 \begin{array}{ll} \max _{\alpha_{1}} & \alpha_{1}^{\mathrm{T}} \Sigma \alpha_{1} \\ \text { s.t. } & \alpha_{1}^{\mathrm{T}} \alpha_{1}=1 \end{array} maxα1 s.t. α1TΣα1α1Tα1=1
定义拉格朗日函数
a r g m a x α L = α 1 T Σ α 1 − λ ( α 1 T α 1 − 1 ) argmax_\alpha L=\alpha_{1}^{\mathrm{T}} \Sigma \alpha_{1}-\lambda\left(\alpha_{1}^{\mathrm{T}} \alpha_{1}-1\right) argmaxαL=α1TΣα1−λ(α1Tα1−1)
其中 λ \lambda λ 是拉格朗日乘子。将拉格朗日函数L对 α 1 \alpha_{1} α1 求导, 并令其为 0 , 得
Σ α 1 − λ α 1 = 0 \Sigma \alpha_{1}-\lambda \alpha_{1}=0 Σα1−λα1=0
因此, λ \lambda λ 是 Σ \Sigma Σ 的特征值, α 1 \alpha_{1} α1 是对应的单位特征向量。于是, 目标函数
α 1 T Σ α 1 = α 1 T λ α 1 = λ α 1 T α 1 = λ \alpha_{1}^{\mathrm{T}} \Sigma \alpha_{1}=\alpha_{1}^{\mathrm{T}} \lambda \alpha_{1}=\lambda \alpha_{1}^{\mathrm{T}} \alpha_{1}=\lambda α1TΣα1=α1Tλα1=λα1Tα1=λ
假设 α 1 \alpha_{1} α1 是 Σ \Sigma Σ 的最大特征值 λ 1 \lambda_{1} λ1 对应的单位特征向量, 显然 α 1 \alpha_{1} α1 与 λ 1 \lambda_{1} λ1 是最优化问题的解。所以, α 1 T x \alpha_{1}^{\mathrm{T}} x α1Tx 构成第一主成分, 其方差等于协方差矩阵的最大特征值
var ( α 1 T x ) = α 1 T Σ α 1 = λ 1 \operatorname{var}\left(\alpha_{1}^{\mathrm{T}} \boldsymbol{x}\right)=\alpha_{1}^{\mathrm{T}} \Sigma \alpha_{1}=\lambda_{1} var(α1Tx)=α1TΣα1=λ1
接着求 x \boldsymbol{x} x 的第二主成分 y 2 = α 2 T x ∘ y_{2}=\alpha_{2}^{\mathrm{T}} \boldsymbol{x}_{\circ} y2=α2Tx∘ 第二主成分的 α 2 \alpha_{2} α2 是在 α 2 T α 2 = 1 \alpha_{2}^{\mathrm{T}} \alpha_{2}=1 α2Tα2=1, 且 α 2 T x \alpha_{2}^{\mathrm{T}} \boldsymbol{x} α2Tx 与 α 1 T x \alpha_{1}^{\mathrm{T}} \boldsymbol{x} α1Tx 不相关的条件下, x \boldsymbol{x} x 的所有线性变换中使方差
var ( α 2 T x ) = α 2 T Σ α 2 \operatorname{var}\left(\alpha_{2}^{\mathrm{T}} \boldsymbol{x}\right)=\alpha_{2}^{\mathrm{T}} \Sigma \alpha_{2} var(α2Tx)=α2TΣα2
达到最大的。
求第二主成分需要求解约束最优化问题max α 2 α 2 T Σ α 2 s.t. α 1 T Σ α 2 = 0 , α 2 T Σ α 1 = 0 α 2 T α 2 = 1 \begin{array}{ll} \max _{\alpha_{2}} & \alpha_{2}^{\mathrm{T}} \Sigma \alpha_{2} \\ \text { s.t. } & \alpha_{1}^{\mathrm{T}} \Sigma \alpha_{2}=0, \quad \alpha_{2}^{\mathrm{T}} \Sigma \alpha_{1}=0 \\ & \alpha_{2}^{\mathrm{T}} \alpha_{2}=1 \end{array} maxα2 s.t. α2TΣα2α1TΣα2=0,α2TΣα1=0α2Tα2=1
注意到
α 1 T Σ α 2 = α 2 T Σ α 1 = α 2 T λ 1 α 1 = λ 1 α 2 T α 1 = λ 1 α 1 T α 2 \alpha_{1}^{\mathrm{T}} \Sigma \alpha_{2}=\alpha_{2}^{\mathrm{T}} \Sigma \alpha_{1}=\alpha_{2}^{\mathrm{T}} \lambda_{1} \alpha_{1}=\lambda_{1} \alpha_{2}^{\mathrm{T}} \alpha_{1}=\lambda_{1} \alpha_{1}^{\mathrm{T}} \alpha_{2} α1TΣα2=α2TΣα1=α2Tλ1α1=λ1α2Tα1=λ1α1Tα2
以及
α 1 T α 2 = 0 , α 2 T α 1 = 0 \alpha_{1}^{\mathrm{T}} \alpha_{2}=0, \quad \alpha_{2}^{\mathrm{T}} \alpha_{1}=0 α1Tα2=0,α2Tα1=0
定义拉格朗日函数
α 2 T Σ α 2 − λ ( α 2 T α 2 − 1 ) − ϕ α 2 T α 1 \alpha_{2}^{\mathrm{T}} \Sigma \alpha_{2}-\lambda\left(\alpha_{2}^{\mathrm{T}} \alpha_{2}-1\right)-\phi \alpha_{2}^{\mathrm{T}} \alpha_{1} α2TΣα2−λ(α2Tα2−1)−ϕα2Tα1
其中 λ , ϕ \lambda, \phi λ,ϕ 是拉格朗日乘子。对 α 2 \alpha_{2} α2 求导, 并令其为 0 , 得
2 Σ α 2 − 2 λ α 2 − ϕ α 1 = 0 2 \Sigma \alpha_{2}-2 \lambda \alpha_{2}-\phi \alpha_{1}=0 2Σα2−2λα2−ϕα1=0
将方程左乘以 α 1 T \alpha_{1}^{\mathrm{T}} α1T 有
2 α 1 T Σ α 2 − 2 λ α 1 T α 2 − ϕ α 1 T α 1 = 0 2 \alpha_{1}^{\mathrm{T}} \Sigma \alpha_{2}-2 \lambda \alpha_{1}^{\mathrm{T}} \alpha_{2}-\phi \alpha_{1}^{\mathrm{T}} \alpha_{1}=0 2α1TΣα2−2λα1Tα2−ϕα1Tα1=0
α 1 T α 2 \alpha_1^T\alpha_2 α1Tα2正交=> α 1 T α 2 = 0 \alpha_1^T\alpha_2=0 α1Tα2=0
此式前两项为 0 , 且 α 1 T α 1 = 1 \alpha_{1}^{\mathrm{T}} \alpha_{1}=1 α1Tα1=1, 导出 ϕ = 0 \phi=0 ϕ=0, 由此推出
Σ α 2 − λ α 2 = 0 \Sigma \alpha_{2}-\lambda \alpha_{2}=0 Σα2−λα2=0
由此, λ \lambda λ 是 Σ \Sigma Σ 的特征值, α 2 \alpha_{2} α2 是对应的单位特征向量。于是, 目标函数
α 2 T Σ α 2 = α 2 T λ α 2 = λ α 2 T α 2 = λ \alpha_{2}^{\mathrm{T}} \Sigma \alpha_{2}=\alpha_{2}^{\mathrm{T}} \lambda \alpha_{2}=\lambda \alpha_{2}^{\mathrm{T}} \alpha_{2}=\lambda α2TΣα2=α2Tλα2=λα2Tα2=λ
假设 α 2 \alpha_{2} α2 是 Σ \Sigma Σ 的第二大特征值 λ 2 \lambda_{2} λ2 对应的单位特征向量, 显然 α 2 \alpha_{2} α2 与 λ 2 \lambda_{2} λ2 是以上 最优化问题的解 E \mathbb{E} E 。于是 α 2 T x \alpha_{2}^{\mathrm{T}} \boldsymbol{x} α2Tx 构成第二主成分, 其方差等于协方差矩阵的第二大特征值，
var ( α 2 T x ) = α 2 T Σ α 2 = λ 2 \operatorname{var}\left(\alpha_{2}^{\mathrm{T}} \boldsymbol{x}\right)=\alpha_{2}^{\mathrm{T}} \Sigma \alpha_{2}=\lambda_{2} var(α2Tx)=α2TΣα2=λ2
一般地, x \boldsymbol{x} x 的第 k k k 主成分是 α k T x \alpha_{k}^{\mathrm{T}} \boldsymbol{x} αkTx, 并且 var ( α k T x ) = λ k \operatorname{var}\left(\alpha_{k}^{\mathrm{T}} \boldsymbol{x}\right)=\lambda_{k} var(αkTx)=λk, 这里 λ k \lambda_{k} λk 是 Σ \Sigma Σ 的第 k k k 个 特征值并且 α k \alpha_{k} αk 是对应的单位特征向量。可以从个第 k − 1 k-1 k−1 主成分出发递推证明第 k k k 个主成分的情况。
按照上述方法求得第一、第二、直到第 m m m 主成分, 其系数向量 α 1 , α 2 , ⋯ , α m \alpha_{1}, \alpha_{2}, \cdots, \alpha_{m} α1,α2,⋯,αm 分别是 Σ \Sigma Σ 的第一个、第二个、直到第 m m m 个单位特征向量, λ 1 , λ 2 , ⋯ , λ m \lambda_{1}, \lambda_{2}, \cdots, \lambda_{m} λ1,λ2,⋯,λm 分别是对应 的特征值。并且, 第 k k k 主成分的方差等于 Σ \Sigma Σ 的第 k k k 个特征值，var ( α k T x ) = α k T Σ α k = λ k , k = 1 , 2 , ⋯ , m \operatorname{var}\left(\alpha_{k}^{\mathrm{T}} x\right)=\alpha_{k}^{\mathrm{T}} \Sigma \alpha_{k}=\lambda_{k}, \quad k=1,2, \cdots, m var(αkTx)=αkTΣαk=λk,k=1,2,⋯,m
定理证毕。
推论
16.1
m
16.1 \quad m
16.1m 维随机变量
y
=
(
y
1
,
y
2
,
⋯
,
y
m
)
T
\boldsymbol{y}=\left(y_{1}, y_{2}, \cdots, y_{m}\right)^{\mathrm{T}}
y=(y1,y2,⋯,ym)T 的分量依次是
x
\boldsymbol{x}
x 的第一主成分 到第
m
m
m 主成分的充要条件是:
(1)
y
=
A
T
x
,
A
\boldsymbol{y}=A^{\mathrm{T}} \boldsymbol{x}, A
y=ATx,A 为正交矩阵
A = [ α 11 α 12 ⋯ α 1 m α 21 α 22 ⋯ α 2 m ⋮ ⋮ ⋮ α m 1 α m 2 ⋯ α m m ] A=\left[\begin{array}{cccc} \alpha_{11} & \alpha_{12} & \cdots & \alpha_{1 m} \\ \alpha_{21} & \alpha_{22} & \cdots & \alpha_{2 m} \\ \vdots & \vdots & & \vdots \\ \alpha_{m 1} & \alpha_{m 2} & \cdots & \alpha_{m m} \end{array}\right] A=⎣⎢⎢⎢⎡α11α21⋮αm1α12α22⋮αm2⋯⋯⋯α1mα2m⋮αmm⎦⎥⎥⎥⎤
(1) 为了叙述方便，这里将变量和其最优值用同一符号表示。
(2)
y
\boldsymbol{y}
y 的协方差矩阵为对角矩阵
cov ( y ) = diag ( λ 1 , λ 2 , ⋯ , λ m ) λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ λ m \begin{aligned} &\operatorname{cov}(\boldsymbol{y})=\operatorname{diag}\left(\lambda_{1}, \lambda_{2}, \cdots, \lambda_{m}\right) \\ &\lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{m} \end{aligned} cov(y)=diag(λ1,λ2,⋯,λm)λ1⩾λ2⩾⋯⩾λm
其中 λ k \lambda_{k} λk 是 Σ \Sigma Σ 的第 k k k 个特征值, α k \alpha_{k} αk 是对应的单位特征向量, k = 1 , 2 , ⋯ , m k=1,2, \cdots, m k=1,2,⋯,m. 以上证明中, λ k \lambda_{k} λk 是 Σ \Sigma Σ 的第 k k k 个特征值, α k \alpha_{k} αk 是对应的单位特征向量, 即
Σ α k = λ k α k , k = 1 , 2 , ⋯ , m \Sigma \alpha_{k}=\lambda_{k} \alpha_{k}, \quad k=1,2, \cdots, m Σαk=λkαk,k=1,2,⋯,m
用矩阵表示即为
Σ A = A Λ \Sigma A=A \Lambda ΣA=AΛ
这里 A = [ α i j ] m × m , Λ A=\left[\alpha_{i j}\right]_{m \times m}, \Lambda A=[αij]m×m,Λ 是对角矩阵, 其第 k k k 个对角元素是 λ k \lambda_{k} λk. 因为 A A A 是正交矩阵, 即 A T A = A A T = I A^{\mathrm{T}} A=A A^{\mathrm{T}}=I ATA=AAT=I, 由式 ( 16.14 ) (16.14) (16.14) 得到两个公式
A T Σ A = Λ A^{\mathrm{T}} \Sigma A=\Lambda ATΣA=Λ
和
Σ = A Λ A T \Sigma=A \Lambda A^{\mathrm{T}} Σ=AΛAT
总体主成分的性质：
(1)总体主成分
y
\boldsymbol{y}
y 的协方差矩阵是对角矩阵
cov ( y ) = Λ = diag ( λ 1 , λ 2 , ⋯ , λ m ) \operatorname{cov}(\boldsymbol{y})=\Lambda=\operatorname{diag}\left(\lambda_{1}, \lambda_{2}, \cdots, \lambda_{m}\right) cov(y)=Λ=diag(λ1,λ2,⋯,λm)
(2)总体主成分 y \boldsymbol{y} y 的方差之和等于随机变量 x \boldsymbol{x} x 的方差之和, 即
∑ i = 1 m λ i = ∑ i = 1 m σ i i \sum_{i=1}^{m} \lambda_{i}=\sum_{i=1}^{m} \sigma_{i i} i=1∑mλi=i=1∑mσii
其中 σ i i \sigma_{i i} σii 是随机变量 x i x_{i} xi 的方差, 即协方差矩阵 Σ \Sigma Σ 的对角元素。事实上, 利用式 Σ = A Λ A T \Sigma=A \Lambda A^{\mathrm{T}} Σ=AΛAT 及矩阵的迹 (trace) 的性质, 可知
t r ( A B ) = t r ( B A ) tr(AB)=tr(BA) tr(AB)=tr(BA)
A A A为正交矩阵==> A T A = I A^TA=I ATA=I
∑ i = 1 m var ( x i ) = tr ( Σ T ) = tr ( A Λ A T ) = tr ( A T A Λ ) = tr ( A T Λ A ) = tr ( Λ ) = ∑ i = 1 m λ i = ∑ i = 1 m var ( y i ) \begin{aligned} \sum_{i=1}^{m} \operatorname{var}\left(x_{i}\right) &=\operatorname{tr}\left(\Sigma^{\mathrm{T}}\right)=\operatorname{tr}\left(A \Lambda A^{\mathrm{T}}\right)=\operatorname{tr}\left(A^{\mathrm{T}}A \Lambda \right)=\operatorname{tr}\left(A^{\mathrm{T}} \Lambda A\right) \\ &=\operatorname{tr}(\Lambda)=\sum_{i=1}^{m} \lambda_{i}=\sum_{i=1}^{m} \operatorname{var}\left(y_{i}\right) \end{aligned} i=1∑mvar(xi)=tr(ΣT)=tr(AΛAT)=tr(ATAΛ)=tr(ATΛA)=tr(Λ)=i=1∑mλi=i=1∑mvar(yi)
(3) 第 k k k 个主成分 y k y_{k} yk 与变量 x i x_{i} xi 的相关系数 ρ ( y k , x i ) \rho\left(y_{k}, x_{i}\right) ρ(yk,xi) 称为因子负荷量 (factor loading）, 它表示第 k k k 个主成分 y k y_{k} yk 与变量 x i x_{i} xi 的相关关系。计算公式是
ρ ( y k , x i ) = λ k α i k σ i i , k , i = 1 , 2 , ⋯ , m ① \rho\left(y_{k}, x_{i}\right)=\frac{\sqrt{\lambda_{k}} \alpha_{i k}}{\sqrt{\sigma_{i i}}}, \quad k, i=1,2, \cdots, m① ρ(yk,xi)=σii λk αik,k,i=1,2,⋯,m①
因为
ρ ( y k , x i ) = cov ( y k , x i ) var ( y k ) var ( x i ) = cov ( α k T x , e i T x ) λ k σ i i \rho\left(y_{k}, x_{i}\right)=\frac{\operatorname{cov}\left(y_{k}, x_{i}\right)}{\sqrt{\operatorname{var}\left(y_{k}\right) \operatorname{var}\left(x_{i}\right)}}=\frac{\operatorname{cov}\left(\alpha_{k}^{\mathrm{T}} \boldsymbol{x}, e_{i}^{\mathrm{T}} \boldsymbol{x}\right)}{\sqrt{\lambda_{k}} \sqrt{\sigma_{i i}}} ρ(yk,xi)=var(yk)var(xi) cov(yk,xi)=λk σii cov(αkTx,eiTx)
其中 e i e_{i} ei 为基本单位向量, 其第 i i i 个分量为 1, 其余为 0 。再由协方差的性质
cov ( α k T x , e i T x ) = α k T Σ e i = e i T Σ α k = λ k e i T α k = λ k α i k \operatorname{cov}\left(\alpha_{k}^{\mathrm{T}} \boldsymbol{x}, e_{i}^{\mathrm{T}} \boldsymbol{x}\right)=\alpha_{k}^{\mathrm{T}} \Sigma e_{i}=e_{i}^{\mathrm{T}} \Sigma \alpha_{k}=\lambda_{k} e_{i}^{\mathrm{T}} \alpha_{k}=\lambda_{k} \alpha_{i k} cov(αkTx,eiTx)=αkTΣei=eiTΣαk=λkeiTαk=λkαik
故得式 ① ① ① 。
(4) 第 k k k 个主成分 y k y_{k} yk 与 m m m 个变量的因子负荷量满足
∑ i = 1 m σ i i ρ 2 ( y k , x i ) = λ k \sum_{i=1}^{m} \sigma_{i i} \rho^{2}\left(y_{k}, x_{i}\right)=\lambda_{k} i=1∑mσiiρ2(yk,xi)=λk
由式 ① ① ① 有
∑ i = 1 m σ i i ρ 2 ( y k , x i ) = ∑ i = 1 m λ k α i k 2 = λ k α k T α k = λ k \sum_{i=1}^{m} \sigma_{i i} \rho^{2}\left(y_{k}, x_{i}\right)=\sum_{i=1}^{m} \lambda_{k} \alpha_{i k}^{2}=\lambda_{k} \alpha_{k}^{\mathrm{T}} \alpha_{k}=\lambda_{k} i=1∑mσiiρ2(yk,xi)=i=1∑mλkαik2=λkαkTαk=λk
（5） m m m 个主成分与第 i i i 个变量 x i x_{i} xi 的因子负荷量满足
∑ k = 1 m ρ 2 ( y k , x i ) = 1 ② \sum_{k=1}^{m} \rho^{2}\left(y_{k}, x_{i}\right)=1② k=1∑mρ2(yk,xi)=1②
由于 y 1 , y 2 , ⋯ , y m y_{1}, y_{2}, \cdots, y_{m} y1,y2,⋯,ym 互不相关, 故
ρ 2 ( x i , ( y 1 , y 2 , ⋯ , y m ) ) = ∑ k = 1 m ρ 2 ( y k , x i ) \rho^{2}\left(x_{i},\left(y_{1}, y_{2}, \cdots, y_{m}\right)\right)=\sum_{k=1}^{m} \rho^{2}\left(y_{k}, x_{i}\right) ρ2(xi,(y1,y2,⋯,ym))=k=1∑mρ2(yk,xi)
又因 x i x_{i} xi 可以表为 y 1 , y 2 , ⋯ , y m y_{1}, y_{2}, \cdots, y_{m} y1,y2,⋯,ym 的线性组合, 所以 x i x_{i} xi 与 y 1 , y 2 , ⋯ , y m y_{1}, y_{2}, \cdots, y_{m} y1,y2,⋯,ym 的相关系数的 平方为 1 , 即
ρ 2 ( x i , ( y 1 , y 2 , ⋯ , y m ) ) = 1 \rho^{2}\left(x_{i},\left(y_{1}, y_{2}, \cdots, y_{m}\right)\right)=1 ρ2(xi,(y1,y2,⋯,ym))=1
故得式 ② ② ② 。
16.1.4 主成分的个数
定理 16.2 16.2 16.2 对任意正整数 q , 1 ⩽ q ⩽ m q, 1 \leqslant q \leqslant m q,1⩽q⩽m, 考虑正交线性变换
y = B T x y=B^{\mathrm{T}} \boldsymbol{x} y=BTx
其中 y y y 是 q q q 维向量, B T B^{\mathrm{T}} BT 是 q × m q \times m q×m 矩阵, 令 y y y 的协方差矩阵为
Σ y = B T Σ B \Sigma_{\boldsymbol{y}}=B^{\mathrm{T}} \Sigma B Σy=BTΣB
则 Σ y \Sigma_{y} Σy 的迹 tr ( Σ y ) \operatorname{tr}\left(\Sigma_{y}\right) tr(Σy) 在 B = A q B=A_{q} B=Aq 时取得最大值, 其中矩阵 A q A_{q} Aq 由正交矩阵 A A A 的前 q q q 列组成。
证明 令 β k \beta_{k} βk 是 B B B 的第 k k k 列, 由于正交矩阵 A A A 的列构成 m m m 维空间的基, 所以 β k \beta_{k} βk 可以由 A A A 的列表示, 即
β k = ∑ j = 1 m c j k α j , k = 1 , 2 , ⋯ , q \beta_{k}=\sum_{j=1}^{m} c_{j k} \alpha_{j}, \quad k=1,2, \cdots, q βk=j=1∑mcjkαj,k=1,2,⋯,q
等价地
B = A C B=A C B=AC
其中 C C C 是 m × q m \times q m×q 矩阵, 其第 j j j 行第 k k k 列元素为 c j k c_{j k} cjk 。
首先,B T Σ B = C T A T Σ A C = C T Λ C = ∑ j = 1 m λ j c j c j T B^{\mathrm{T}} \Sigma B=C^{\mathrm{T}} A^{\mathrm{T}} \Sigma A C=C^{\mathrm{T}} \Lambda C=\sum_{j=1}^{m} \lambda_{j} c_{j} c_{j}^{\mathrm{T}} BTΣB=CTATΣAC=CTΛC=j=1∑mλjcjcjT
其中 c j T c_{j}^{\mathrm{T}} cjT 是 C C C 的第 j j j 行。因此
tr ( B T Σ B ) = ∑ j = 1 m λ j tr ( c j c j T ) = ∑ j = 1 m λ j tr ( c j T c j ) = ∑ j = 1 m λ j c j T c j = ∑ j = 1 m ∑ k = 1 q λ j c j k 2 \begin{aligned} \operatorname{tr}\left(B^{\mathrm{T}} \Sigma B\right) &=\sum_{j=1}^{m} \lambda_{j} \operatorname{tr}\left(c_{j} c_{j}^{\mathrm{T}}\right) \\ &=\sum_{j=1}^{m} \lambda_{j} \operatorname{tr}\left(c_{j}^{\mathrm{T}} c_{j}\right) \\ &=\sum_{j=1}^{m} \lambda_{j} c_{j}^{\mathrm{T}} c_{j} \\ &=\sum_{j=1}^{m} \sum_{k=1}^{q} \lambda_{j} c_{j k}^{2} \end{aligned} tr(BTΣB)=j=1∑mλjtr(cjcjT)=j=1∑mλjtr(cjTcj)=j=1∑mλjcjTcj=j=1∑mk=1∑qλjcjk2
其次, 由式 B = A C B=A C B=AC 及 A A A 的正交性知
C = A T B C=A^{\mathrm{T}} B C=ATB
由于 A A A 是正交的, B B B 的列是正交的, 所以
C T C = B T A A T B = B T B = I q C^{\mathrm{T}} C=B^{\mathrm{T}} A A^{\mathrm{T}} B=B^{\mathrm{T}} B=I_{q} CTC=BTAATB=BTB=Iq
即 C C C 的列也是正交的。于是
tr ( C T C ) = tr ( I q ) = q ∑ j = 1 m ∑ k = 1 q c j k 2 = q \begin{aligned} &\operatorname{tr}\left(C^{\mathrm{T}} C\right)=\operatorname{tr}\left(I_{q}\right)=q \\ &\sum_{j=1}^{m} \sum_{k=1}^{q} c_{j k}^{2}=q \end{aligned} tr(CTC)=tr(Iq)=qj=1∑mk=1∑qcjk2=q
这样, 矩阵 C C C 可以认为是某个 m m m 阶正交矩阵 D D D 的前 q q q 列。正交矩阵 D D D 的行也 正交, 所以满足
d j T d j = 1 , j = 1 , 2 , ⋯ , m d_{j}^{\mathrm{T}} d_{j}=1, \quad j=1,2, \cdots, m djTdj=1,j=1,2,⋯,m
其中 d j T d_{j}^{\mathrm{T}} djT 是 D D D 的第 j j j 行。由于矩阵 D D D 的行包括矩阵 C C C 的行的前 q q q 个元素, 所以
c j T c j ⩽ 1 , j = 1 , 2 , ⋯ , m c_{j}^{\mathrm{T}} c_{j} \leqslant 1, \quad j=1,2, \cdots, m cjTcj⩽1,j=1,2,⋯,m
即
∑ k = 1 q c j k 2 ⩽ 1 , j = 1 , 2 , ⋯ , m \sum_{k=1}^{q} c_{j k}^{2} \leqslant 1, \quad j=1,2, \cdots, m k=1∑qcjk2⩽1,j=1,2,⋯,m
注意到在式 ( 16.26 ) (16.26) (16.26) 中 ∑ k = 1 q c j k 2 \sum_{k=1}^{q} c_{j k}^{2} ∑k=1qcjk2 是 λ j \lambda_{j} λj 的系数, 由式 ∑ j = 1 m ∑ k = 1 q c j k 2 = q \sum_{j=1}^{m} \sum_{k=1}^{q} c_{j k}^{2}=q ∑j=1m∑k=1qcjk2=q 这些系数之和是 q q q, 且 由式 ∑ k = 1 q c j k 2 ⩽ 1 , j = 1 , 2 , ⋯ , m \sum_{k=1}^{q} c_{j k}^{2} \leqslant 1, \quad j=1,2, \cdots, m ∑k=1qcjk2⩽1,j=1,2,⋯,m 知这些系数小于等于 1 。因为 λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ λ q ⩾ ⋯ ⩾ λ m \lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{q} \geqslant \cdots \geqslant \lambda_{m} λ1⩾λ2⩾⋯⩾λq⩾⋯⩾λm, 显然, 当能找到 c j k c_{j k} cjk 使得
∑ k = 1 q c j k 2 = { 1 , j = 1 , ⋯ , q 0 , j = q + 1 , ⋯ , m ① \sum_{k=1}^{q} c_{j k}^{2}= \begin{cases}1, & j=1, \cdots, q \\ 0, & j=q+1, \cdots, m\end{cases}① k=1∑qcjk2={1,0,j=1,⋯,qj=q+1,⋯,m①
时， ∑ j = 1 m ( ∑ k = 1 q c j k 2 ) λ j \sum_{j=1}^{m}\left(\sum_{k=1}^{q} c_{j k}^{2}\right) \lambda_{j} ∑j=1m(∑k=1qcjk2)λj 最大。而当 B = A q B=A_{q} B=Aq 时, 有
c j k = { 1 , 1 ⩽ j = k ⩽ q 0 , 其他 c_{j k}= \begin{cases}1, & 1 \leqslant j=k \leqslant q \\ 0, & \text { 其他 }\end{cases} cjk={1,0,1⩽j=k⩽q 其他
满足式 ① ① ① 。所以，当 B = A q B=A_{q} B=Aq 时, tr ( Σ y ) \operatorname{tr}\left(\Sigma_{y}\right) tr(Σy) 达到最大值。
定理 16.2 16.2 16.2 表明, 当 x x x 的线性变换 y y y 在 B = A q B=A_{q} B=Aq 时, 其协方差矩阵 Σ y \Sigma_{y} Σy 的迹 tr ( Σ y ) \operatorname{tr}\left(\Sigma_{y}\right) tr(Σy) 取得最大值, 这就是说, 当取 A A A 的前 q q q 列取 x x x 的前 q q q 个主成分时, 能够最大限度地保留原有变量方差的信息。
定理 16.3 16.3 16.3 考虑正交变换
y = B T x \boldsymbol{y}=B^{\mathrm{T}} \boldsymbol{x} y=BTx
这里
B
T
B^{\mathrm{T}}
BT 是
p
×
m
p \times m
p×m 矩阵,
A
A
A 和
Σ
y
\Sigma_{y}
Σy 的定义与定理
16.2
16.2
16.2 相同, 则
tr
(
Σ
y
)
\operatorname{tr}\left(\Sigma_{y}\right)
tr(Σy) 在
B
=
A
p
B=A_{p}
B=Ap 时 取得最小值, 其中矩阵
A
p
A_{p}
Ap 由
A
A
A 的后
p
p
p 列组成。
定理
16.3
16.3
16.3 可以理解为, 当 舍弃
A
A
A 的后
p
p
p 列, 即舍弃变量
x
x
x 的后
p
p
p 个主成分时, 原有变量的方差的信息损失最少。
以上两个定理作为选择
k
k
k 个主成分的理论依据。
具体选择
k
k
k 的方法, 通常利用方差贡献率。
定义
16.2
16.2
16.2 第
k
k
k 主成分
y
k
y_{k}
yk 的方差贡献率定义为
y
k
y_{k}
yk 的方差与所有方差之和的比, 记作
η
k
\eta_{k}
ηk
η k = λ k ∑ i = 1 m λ i \eta_{k}=\frac{\lambda_{k}}{\sum_{i=1}^{m} \lambda_{i}} ηk=∑i=1mλiλk
k k k 个主成分 y 1 , y 2 , ⋯ , y k y_{1}, y_{2}, \cdots, y_{k} y1,y2,⋯,yk 的累计方差贡献率定义为 k k k 个方差之和与所有方差之和的比
∑ i = 1 k η i = ∑ i = 1 k λ i ∑ i = 1 m λ i \sum_{i=1}^{k} \eta_{i}=\frac{\sum_{i=1}^{k} \lambda_{i}}{\sum_{i=1}^{m} \lambda_{i}} i=1∑kηi=∑i=1mλi∑i=1kλi
通常取
k
k
k 使得累计方差贡献率达到规定的百分比以上, 例如
70
%
∼
80
%
70 \% \sim 80 \%
70%∼80% 以上。累计方差贡献率反映了主成分保留信息的比例, 但它不能反映对某个原有变量
x
i
x_{i}
xi 保留信息的比例, 这时通常利用
k
k
k 个主成分
y
1
,
y
2
,
⋯
,
y
k
y_{1}, y_{2}, \cdots, y_{k}
y1,y2,⋯,yk 对原有变量
x
i
x_{i}
xi 的贡献率。
定义
16.3
k
16.3 \quad k
16.3k 个主成分
y
1
,
y
2
,
⋯
,
y
k
y_{1}, y_{2}, \cdots, y_{k}
y1,y2,⋯,yk 对原有变量
x
i
x_{i}
xi 的贡献率定义为
x
i
x_{i}
xi 与
(
y
1
,
y
2
,
⋯
,
y
k
)
\left(y_{1}, y_{2}, \cdots, y_{k}\right)
(y1,y2,⋯,yk) 的相关系数的平方, 记作
ν
i
\nu_{i}
νi
ν i = ρ 2 ( x i , ( y 1 , y 2 , ⋯ , y k ) ) \nu_{i}=\rho^{2}\left(x_{i},\left(y_{1}, y_{2}, \cdots, y_{k}\right)\right) νi=ρ2(xi,(y1,y2,⋯,yk))
计算公式如下:
ν i = ρ 2 ( x i , ( y 1 , y 2 , ⋯ , y k ) ) = ∑ j = 1 k ρ 2 ( x i , y j ) = ∑ j = 1 k λ j α i j 2 σ i i \nu_{i}=\rho^{2}\left(x_{i},\left(y_{1}, y_{2}, \cdots, y_{k}\right)\right)=\sum_{j=1}^{k} \rho^{2}\left(x_{i}, y_{j}\right)=\sum_{j=1}^{k} \frac{\lambda_{j} \alpha_{i j}^{2}}{\sigma_{i i}} νi=ρ2(xi,(y1,y2,⋯,yk))=j=1∑kρ2(xi,yj)=j=1∑kσiiλjαij2
16.1 .5 规范化变量的总体主成分
在实际问题中，直接求主成分有时会产生不合理的结果。为了消除这个影响, 常常对各个随机变量实施规范化, 使其均值为 0 , 方差为 1 。
设
x
=
(
x
1
,
x
2
,
⋯
,
x
m
)
T
x=\left(x_{1}, x_{2}, \cdots, x_{m}\right)^{\mathrm{T}}
x=(x1,x2,⋯,xm)T 为
m
m
m 维随机变量,
x
i
x_{i}
xi 为第
i
i
i 个随机变量,
i
=
i=
i=
1
,
2
,
⋯
,
m
1,2, \cdots, m
1,2,⋯,m, 令
x i ∗ = x i − E ( x i ) var ( x i ) , i = 1 , 2 , ⋯ , m x_{i}^{*}=\frac{x_{i}-E\left(x_{i}\right)}{\sqrt{\operatorname{var}\left(x_{i}\right)}}, \quad i=1,2, \cdots, m xi∗=var(xi) xi−E(xi),i=1,2,⋯,m
其中
E
(
x
i
)
,
var
(
x
i
)
E\left(x_{i}\right), \operatorname{var}\left(x_{i}\right)
E(xi),var(xi) 分别是随机变量
x
i
x_{i}
xi 的均值和方差, 这时
x
i
∗
x_{i}^{*}
xi∗ 就是
x
i
x_{i}
xi 的规范化随机 变量。
显然, 规范化随机变量的协方差矩阵就是相关矩阵
R
∘
R_{\circ}
R∘ 主成分分析通常在规范化 随机变量的协方差矩阵即相关矩阵上进行。
对照总体主成分的性质可知, 规范化随机变量的总体主成分有以下性质：
（1）规范化变量主成分的协方差矩阵是
Λ ∗ = diag ( λ 1 ∗ , λ 2 ∗ , ⋯ , λ m ∗ ) \Lambda^{*}=\operatorname{diag}\left(\lambda_{1}^{*}, \lambda_{2}^{*}, \cdots, \lambda_{m}^{*}\right) Λ∗=diag(λ1∗,λ2∗,⋯,λm∗)
其中
λ
1
∗
⩾
λ
2
∗
⩾
⋯
⩾
λ
m
∗
⩾
0
\lambda_{1}^{*} \geqslant \lambda_{2}^{*} \geqslant \cdots \geqslant \lambda_{m}^{*} \geqslant 0
λ1∗⩾λ2∗⩾⋯⩾λm∗⩾0 为相关矩阵
R
R
R 的特征值。
（2）协方差矩阵的特征值之和为
m
m
m
∑ k = 1 m λ k ∗ = m \sum_{k=1}^{m} \lambda_{k}^{*}=m k=1∑mλk∗=m
（3）规范化随机变量 x i ∗ x_{i}^{*} xi∗ 与主成分 y k ∗ y_{k}^{*} yk∗ 的相关系数 ( ( ( 因子负荷量 ) ) ) 为
ρ ( y k ∗ , x i ∗ ) = λ k ∗ e i k ∗ , k , i = 1 , 2 , ⋯ , m \rho\left(y_{k}^{*}, x_{i}^{*}\right)=\sqrt{\lambda_{k}^{*}} e_{i k}^{*}, \quad k, i=1,2, \cdots, m ρ(yk∗,xi∗)=λk∗ eik∗,k,i=1,2,⋯,m
其中
e
k
∗
=
(
e
1
k
∗
,
e
2
k
∗
,
⋯
,
e
m
k
∗
)
T
e_{k}^{*}=\left(e_{1 k}^{*}, e_{2 k}^{*}, \cdots, e_{m k}^{*}\right)^{\mathrm{T}}
ek∗=(e1k∗,e2k∗,⋯,emk∗)T 为矩阵
R
R
R 对应于特征值
λ
k
∗
\lambda_{k}^{*}
λk∗ 的单位特征向量。
（4）所有规范化随机变量
x
i
∗
x_{i}^{*}
xi∗ 与主成分
y
k
∗
y_{k}^{*}
yk∗ 的相关系数的平方和等于
λ
k
∗
\lambda_{k}^{*}
λk∗
∑ i = 1 m ρ 2 ( y k ∗ , x i ∗ ) = ∑ i = 1 m λ k ∗ e i k ∗ 2 = λ k ∗ , k = 1 , 2 , ⋯ , m \sum_{i=1}^{m} \rho^{2}\left(y_{k}^{*}, x_{i}^{*}\right)=\sum_{i=1}^{m} \lambda_{k}^{*} e_{i k}^{* 2}=\lambda_{k}^{*}, \quad k=1,2, \cdots, m i=1∑mρ2(yk∗,xi∗)=i=1∑mλk∗eik∗2=λk∗,k=1,2,⋯,m
（5）规范化随机变量 x i ∗ x_{i}^{*} xi∗ 与所有主成分 y k ∗ y_{k}^{*} yk∗ 的相关系数的平方和等于 1
∑ k = 1 m ρ 2 ( y k ∗ , x i ∗ ) = ∑ k = 1 m λ k ∗ e i k ∗ 2 = 1 , i = 1 , 2 , ⋯ , m \sum_{k=1}^{m} \rho^{2}\left(y_{k}^{*}, x_{i}^{*}\right)=\sum_{k=1}^{m} \lambda_{k}^{*} e_{i k}^{* 2}=1, \quad i=1,2, \cdots, m k=1∑mρ2(yk∗,xi∗)=k=1∑mλk∗eik∗2=1,i=1,2,⋯,m
16.2 样本主成分分析
预备知识
X = ( x 1 x 2 x 3 … x n ) n × p = ( x 1 T x 2 T … x N T ) = ( x 11 x 12 … x 1 p x 21 x 22 … x 2 p x N 1 x N 2 … x N p ) n × p X=(x_1 \quad x_2 \quad x_3 \quad … \quad x_n)_{n\times p}=\begin{pmatrix} x_1^T \\ x_2^T\\ …\\ x_N^T \end{pmatrix}= {\begin{pmatrix} x_{11} & x_{12} &…& x_{1p} \\ x_{21} & x_{22} &…& x_{2p} \\ x_{N1} & x_{N2} &…& x_{Np} \end{pmatrix}}_{n\times p} X=(x1x2x3…xn)n×p=⎝⎜⎜⎛x1Tx2T…xNT⎠⎟⎟⎞=⎝⎛x11x21xN1x12x22xN2………x1px2pxNp⎠⎞n×p
x i ∈ R p i = 1 , 2 , … , N x_i\in R^p \quad i=1,2,…,N xi∈Rpi=1,2,…,N
期望： x ˉ = 1 N ∑ i = 1 N x i \bar{x}=\frac{1}{N} \sum_{i=1}^{N} x_{i} xˉ=N1∑i=1Nxi
协方差： S = 1 N ∑ i = 1 N ( x i − x ˉ ) ( x i − x ˉ ) T S=\frac{1}{N} \sum_{i=1}^{N}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T} S=N1∑i=1N(xi−xˉ)(xi−xˉ)T
S = 1 N ∑ i = 1 N ( x i − x ˉ ) ( x i − x ˉ ) T = 1 N ( x 1 − x ˉ , x 2 − x ˉ , ⋯ , x N − x ˉ ) ( x 1 − x ˉ , x 2 − x ˉ , ⋯ , x N − x ˉ ) T = 1 N ( X T − 1 N X T I N 1 I N 1 T ) ( X T − 1 N X T I N 1 I N 1 T ) T = 1 N X T ( E N − 1 N I N 1 I 1 N ) ( E N − 1 N I N 1 I 1 N ) T X = 1 N X T H N H N T X = 1 N X T H N H N X = 1 N X T H X \begin{aligned} S &=\frac{1}{N} \sum_{i=1}^{N}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T} \\ &=\frac{1}{N}\left(x_{1}-\bar{x}, x_{2}-\bar{x}, \cdots, x_{N}-\bar{x}\right)\left(x_{1}-\bar{x}, x_{2}-\bar{x}, \cdots, x_{N}-\bar{x}\right)^{T} \\ &=\frac{1}{N}\left(X^{T}-\frac{1}{N} X^{T} \mathbb{I}_{N 1} \mathbb{I}_{N 1}^{T}\right)\left(X^{T}-\frac{1}{N} X^{T} \mathbb{I}_{N 1} \mathbb{I}_{N 1}^{T}\right)^{T} \\ &=\frac{1}{N} X^{T}\left(E_{N}-\frac{1}{N} \mathbb{I}_{N 1} \mathbb{I}_{1 N}\right)\left(E_{N}-\frac{1}{N} \mathbb{I}_{N 1} \mathbb{I}_{1 N}\right)^{T} X \\ &=\frac{1}{N} X^{T} H_{N} H_{N}^{T} X \\ &=\frac{1}{N} X^{T} H_{N} H_{N} X=\frac{1}{N} X^{T} H X \end{aligned} S=N1i=1∑N(xi−xˉ)(xi−xˉ)T=N1(x1−xˉ,x2−xˉ,⋯,xN−xˉ)(x1−xˉ,x2−xˉ,⋯,xN−xˉ)T=N1(XT−N1XTIN1IN1T)(XT−N1XTIN1IN1T)T=N1XT(EN−N1IN1I1N)(EN−N1IN1I1N)TX=N1XTHNHNTX=N1XTHNHNX=N1XTHX
中心矩阵 H H H
16.2.1 样本主成分的定义与性质
假设对 m m m 维随机变量 x = ( x 1 , x 2 , ⋯ , x m ) T x=\left(x_{1}, x_{2}, \cdots, x_{m}\right)^{\mathrm{T}} x=(x1,x2,⋯,xm)T 进行 n n n 次独立观测, x 1 , x 2 , ⋯ , x n \boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \cdots, \boldsymbol{x}_{n} x1,x2,⋯,xn 表示观测样本, 其中 x j = ( x 1 j , x 2 j , ⋯ , x m j ) T x_{j}=\left(x_{1 j}, x_{2 j}, \cdots, x_{m j}\right)^{\mathrm{T}} xj=(x1j,x2j,⋯,xmj)T 表示第 j j j 个观测样本, x i j x_{i j} xij 表示第 j j j 个 观测样本的第 i i i 个变量, j = 1 , 2 , ⋯ , n j=1,2, \cdots, n j=1,2,⋯,n 。 观测数据用样本矩阵 X \boldsymbol{X} X 表示, 记作
X = [ x 1 x 2 ⋯ x n ] = [ x 11 x 12 ⋯ x 1 n x 21 x 22 ⋯ x 2 n ⋮ ⋮ ⋮ x m 1 x m 2 ⋯ x m n ] X=\left[\begin{array}{llll} x_{1} & x_{2} & \cdots & x_{n} \end{array}\right]=\left[\begin{array}{cccc} x_{11} & x_{12} & \cdots & x_{1 n} \\ x_{21} & x_{22} & \cdots & x_{2 n} \\ \vdots & \vdots & & \vdots \\ x_{m 1} & x_{m 2} & \cdots & x_{m n} \end{array}\right] X=[x1x2⋯xn]=⎣⎢⎢⎢⎡x11x21⋮xm1x12x22⋮xm2⋯⋯⋯x1nx2n⋮xmn⎦⎥⎥⎥⎤
给定样本矩阵 X X X, 可以估计样本均值, 以及样本协方差。样本均值向量 x ˉ \bar{x} xˉ 为
x ˉ = 1 n ∑ j = 1 n x j \bar{x}=\frac{1}{n} \sum_{j=1}^{n} \boldsymbol{x}_{j} xˉ=n1j=1∑nxj
样本协方差矩阵 S S S 为
S = [ s i j ] m × m s i j = 1 n − 1 ∑ k = 1 n ( x i k − x ˉ i ) ( x j k − x ˉ j ) , i , j = 1 , 2 , ⋯ , m \begin{aligned} &S=\left[s_{i j}\right]_{m \times m} \\ &s_{i j}=\frac{1}{n-1} \sum_{k=1}^{n}\left(x_{i k}-\bar{x}_{i}\right)\left(x_{j k}-\bar{x}_{j}\right), \quad i, j=1,2, \cdots, m \end{aligned} S=[sij]m×msij=n−11k=1∑n(xik−xˉi)(xjk−xˉj),i,j=1,2,⋯,m
其中
x
ˉ
i
=
1
n
∑
k
=
1
n
x
i
k
\bar{x}_{i}=\frac{1}{n} \sum_{k=1}^{n} x_{i k}
xˉi=n1∑k=1nxik 为第
i
i
i 个变量的样本均值,
x
ˉ
j
=
1
n
∑
k
=
1
n
x
j
k
\bar{x}_{j}=\frac{1}{n} \sum_{k=1}^{n} x_{j k}
xˉj=n1∑k=1nxjk 为第
j
j
j 个变量的样本 均值。
样本相关矩阵
R
R
R 为
R = [ r i j ] m × m , r i j = s i j s i i s j j , i , j = 1 , 2 , ⋯ , m R=\left[r_{i j}\right]_{m \times m}, \quad r_{i j}=\frac{s_{i j}}{\sqrt{s_{i i} s_{j j}}}, \quad i, j=1,2, \cdots, m R=[rij]m×m,rij=siisjj sij,i,j=1,2,⋯,m
定义 m m m 维向量 x = ( x 1 , x 2 , ⋯ , x m ) T \boldsymbol{x}=\left(x_{1}, x_{2}, \cdots, x_{m}\right)^{\mathrm{T}} x=(x1,x2,⋯,xm)T 到 m m m 维向量 y = ( y 1 , y 2 , ⋯ , y m ) T \boldsymbol{y}=\left(y_{1}, y_{2}, \cdots, y_{m}\right)^{\mathrm{T}} y=(y1,y2,⋯,ym)T 的线性 变换
y = A T x y=A^{\mathrm{T}} \boldsymbol{x} y=ATx
其中
A = [ a 1 a 2 ⋯ a m ] = [ a 11 a 12 ⋯ a 1 m a 21 a 22 ⋯ a 2 m ⋮ ⋮ ⋮ a m 1 a m 2 ⋯ a m m ] a i = ( a 1 i , a 2 i , ⋯ , a m i ) T , i = 1 , 2 , ⋯ , m \begin{aligned} &A=\left[\begin{array}{llll} a_{1} & a_{2} & \cdots & a_{m} \end{array}\right]=\left[\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 m} \\ a_{21} & a_{22} & \cdots & a_{2 m} \\ \vdots & \vdots & & \vdots \\ a_{m 1} & a_{m 2} & \cdots & a_{m m} \end{array}\right] \\ &a_{i}=\left(a_{1 i}, a_{2 i}, \cdots, a_{m i}\right)^{\mathrm{T}}, \quad i=1,2, \cdots, m \end{aligned} A=[a1a2⋯am]=⎣⎢⎢⎢⎡a11a21⋮am1a12a22⋮am2⋯⋯⋯a1ma2m⋮amm⎦⎥⎥⎥⎤ai=(a1i,a2i,⋯,ami)T,i=1,2,⋯,m
考虑式 y = A T x y=A^{\mathrm{T}} \boldsymbol{x} y=ATx 的任意一个线性变换
y i = a i T x = a 1 i x 1 + a 2 i x 2 + ⋯ + a m i x m , i = 1 , 2 , ⋯ , m y_{i}=a_{i}^{\mathrm{T}} x=a_{1 i} x_{1}+a_{2 i} x_{2}+\cdots+a_{m i} x_{m}, \quad i=1,2, \cdots, m yi=aiTx=a1ix1+a2ix2+⋯+amixm,i=1,2,⋯,m
其中 y i y_{i} yi 是 m m m 维向量 y y y 的第 i i i 个变量, 相应于容量为 n n n 的样本 x 1 , x 2 , ⋯ , x n , y i x_{1}, x_{2}, \cdots, x_{n}, y_{i} x1,x2,⋯,xn,yi 的 样本均值 y ˉ i \bar{y}_{i} yˉi 为
y ˉ i = 1 n ∑ j = 1 n a i T x j = a i T x ‾ \bar{y}_{i}=\frac{1}{n} \sum_{j=1}^{n} a_{i}^{\mathrm{T}} \boldsymbol{x}_{j}=a_{i}^{\mathrm{T}} \overline{\boldsymbol{x}} yˉi=n1j=1∑naiTxj=aiTx
其中 x ˉ \bar{x} xˉ 是随机向量 x x x 的样本均值
x ˉ = 1 n ∑ j = 1 n x j \bar{x}=\frac{1}{n} \sum_{j=1}^{n} x_{j} xˉ=n1j=1∑nxj
y i y_{i} yi 的样本方差 var ( y i ) \operatorname{var}\left(y_{i}\right) var(yi) 为
var ( y i ) = 1 n − 1 ∑ j = 1 n ( a i T x j − a i T x ‾ ) 2 = a i T [ 1 n − 1 ∑ j = 1 n ( x j − x ˉ ) ( x j − x ‾ ) T ] a i = a i T S a i \begin{aligned} \operatorname{var}\left(y_{i}\right) &=\frac{1}{n-1} \sum_{j=1}^{n}\left(a_{i}^{\mathrm{T}} \boldsymbol{x}_{j}-a_{i}^{\mathrm{T}} \overline{\boldsymbol{x}}\right)^{2} \\ &=a_{i}^{\mathrm{T}}\left[\frac{1}{n-1} \sum_{j=1}^{n}\left(\boldsymbol{x}_{j}-\bar{x}\right)\left(\boldsymbol{x}_{j}-\overline{\boldsymbol{x}}\right)^{\mathrm{T}}\right] a_{i}=a_{i}^{\mathrm{T}} S a_{i} \end{aligned} var(yi)=n−11j=1∑n(aiTxj−aiTx)2=aiT[n−11j=1∑n(xj−xˉ)(xj−x)T]ai=aiTSai
对任意两个线性变换 y i = α i T x , y k = α k T x y_{i}=\alpha_{i}^{\mathrm{T}} \boldsymbol{x}, y_{k}=\alpha_{k}^{\mathrm{T}} \boldsymbol{x} yi=αiTx,yk=αkTx, 相应于容量为 n n n 的样本 x 1 , x 2 , ⋯ , x n \boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \cdots, \boldsymbol{x}_{n} x1,x2,⋯,xn, y i , y k y_{i}, y_{k} yi,yk 的样本协方差为
cov ( y i , y k ) = a i T S a k \operatorname{cov}\left(y_{i}, y_{k}\right)=a_{i}^{\mathrm{T}} S a_{k} cov(yi,yk)=aiTSak
样本主成分的定义：
定义
16.4
16.4
16.4 (样本主成分)
\quad
给定样本矩阵
X
X
X. 样本第一主成分
y
1
=
a
1
T
x
y_{1}=a_{1}^{\mathrm{T}} \boldsymbol{x}
y1=a1Tx 是在
a
1
T
a
1
=
1
a_{1}^{\mathrm{T}} a_{1}=1
a1Ta1=1 条件下, 使得
a
1
T
x
j
(
j
=
1
,
2
,
⋯
,
n
)
a_{1}^{\mathrm{T}} x_{j}(j=1,2, \cdots, n)
a1Txj(j=1,2,⋯,n) 的样本方差
a
1
T
S
a
1
a_{1}^{\mathrm{T}} S a_{1}
a1TSa1 最大的
x
x
x 的线性 变换; 样本第二主成分
y
2
=
a
2
T
x
y_{2}=a_{2}^{\mathrm{T}} \boldsymbol{x}
y2=a2Tx 是在
a
2
T
a
2
=
1
a_{2}^{\mathrm{T}} a_{2}=1
a2Ta2=1 和
a
2
T
x
j
a_{2}^{\mathrm{T}} x_{j}
a2Txj 与
a
1
T
x
j
(
j
=
1
,
2
,
⋯
,
n
)
a_{1}^{\mathrm{T}} \boldsymbol{x}_{j}(j=1,2, \cdots, n)
a1Txj(j=1,2,⋯,n) 的样本协方差
a
1
T
S
a
2
=
0
a_{1}^{\mathrm{T}} S a_{2}=0
a1TSa2=0 条件下, 使得
a
2
T
x
j
(
j
=
1
,
2
,
⋯
,
n
)
a_{2}^{\mathrm{T}} \boldsymbol{x}_{j}(j=1,2, \cdots, n)
a2Txj(j=1,2,⋯,n) 的样本方差
a
2
T
S
a
2
a_{2}^{\mathrm{T}} S a_{2}
a2TSa2 最大的
x
x
x 的线性变换;一般地, 样本第
i
i
i 主成分
y
i
=
a
i
T
x
y_{i}=a_{i}^{\mathrm{T}} \boldsymbol{x}
yi=aiTx 是在
a
i
T
a
i
=
1
a_{i}^{\mathrm{T}} a_{i}=1
aiTai=1 和
a
i
T
x
j
a_{i}^{\mathrm{T}} \boldsymbol{x}_{j}
aiTxj 与
a
k
T
x
j
(
k
<
i
,
j
=
1
,
2
,
⋯
,
n
)
a_{k}^{\mathrm{T}} \boldsymbol{x}_{j}(k<i, j=1,2, \cdots, n)
akTxj(k<i,j=1,2,⋯,n) 的样本协方差
a
k
T
S
a
i
=
0
a_{k}^{\mathrm{T}} S a_{i}=0
akTSai=0 条件下, 使得
a
i
T
x
j
(
j
=
1
,
2
,
⋯
,
n
)
a_{i}^{\mathrm{T}} x_{j}(j=1,2, \cdots, n)
aiTxj(j=1,2,⋯,n) 的样本方差
a
i
T
S
a
i
a_{i}^{\mathrm{T}} S a_{i}
aiTSai 最大的
x
x
x 的线性变换。
样本主成分与总体主成分具有同样的性质。 这从样本主成分的定义容易看出。只 要以样本协方差矩阵
S
S
S 代替总体协方差矩阵
Σ
\Sigma
Σ 即可。总体主成分的定理
16.2
16.2
16.2 及定理
16.3
16.3
16.3 对样本主成分依然成立。
在使用样本主成分时, 一般假设样本数据是规范化的, 即对样本矩阵作如下变换:
x i j ∗ = x i j − x ˉ i s i i , i = 1 , 2 , ⋯ , m ; j = 1 , 2 , ⋯ , n x_{i j}^{*}=\frac{x_{i j}-\bar{x}_{i}}{\sqrt{s_{i i}}}, \quad i=1,2, \cdots, m ; \quad j=1,2, \cdots, n xij∗=sii xij−xˉi,i=1,2,⋯,m;j=1,2,⋯,n
其中
x ˉ i = 1 n ∑ j = 1 n x i j , i = 1 , 2 , ⋯ , m s i i = 1 n − 1 ∑ j = 1 n ( x i j − x ˉ i ) 2 , i = 1 , 2 , ⋯ , m \begin{aligned} &\bar{x}_{i}=\frac{1}{n} \sum_{j=1}^{n} x_{i j}, \quad i=1,2, \cdots, m \\ &s_{i i}=\frac{1}{n-1} \sum_{j=1}^{n}\left(x_{i j}-\bar{x}_{i}\right)^{2}, \quad i=1,2, \cdots, m \end{aligned} xˉi=n1j=1∑nxij,i=1,2,⋯,msii=n−11j=1∑n(xij−xˉi)2,i=1,2,⋯,m
将规范化变量 x i j ∗ x_{i j}^{*} xij∗ 仍记作 x i j x_{i j} xij, 规范化的样本矩阵仍记作 X X X 。 这时, 样本协方差矩阵 S S S 就是样本相关矩阵 R R R
R = c o v ( x i , x j ) v a r ( x i ) v a r ( x j ) = c o v ( x i , x j ) = 1 n − 1 X X T R=\frac{cov(x_i,x_j)}{\sqrt{var(x_i)var(x_j)}}=cov(x_i,x_j)=\frac{1}{n-1} X X^{\mathrm{T}} R=var(xi)var(xj) cov(xi,xj)=cov(xi,xj)=n−11XXT
样本协方差矩阵 S S S 是总体协方差矩阵 Σ \Sigma Σ 的无偏估计, 样本相关矩阵 R R R 是总体相关矩阵的无偏估计, S S S 的特征值和特征向量是 Σ \Sigma Σ 的特征值和特征向量的极大似然估计。
16.2.2 相关矩阵的特征值分解算法
数据的协方差矩阵或相关矩阵的特征值分解方法。
给定样本矩阵
X
X
X, 利用数据的样本协方差矩阵或者样本相关矩阵的特征值分解进行主成分分析。具体步骤如下:
(1) 对观测数据按式
x
i
j
∗
=
x
i
j
−
x
ˉ
i
s
i
i
,
i
=
1
,
2
,
⋯
,
m
;
j
=
1
,
2
,
⋯
,
n
x_{i j}^{*}=\frac{x_{i j}-\bar{x}_{i}}{\sqrt{s_{i i}}}, \quad i=1,2, \cdots, m ; \quad j=1,2, \cdots, n
xij∗=sii
xij−xˉi,i=1,2,⋯,m;j=1,2,⋯,n 进行规范化处理, 得到规范化数据矩阵, 仍以
X
X
X 表示。
(2) 依据规范化数据矩阵, 计算样本相关矩阵
R
R
R
R = [ r i j ] m × m = 1 n − 1 X X T R=\left[r_{i j}\right]_{m \times m}=\frac{1}{n-1} X X^{\mathrm{T}} R=[rij]m×m=n−11XXT
其中
r i j = 1 n − 1 ∑ l = 1 n x i l x l j , i , j = 1 , 2 , ⋯ , m r_{i j}=\frac{1}{n-1} \sum_{l=1}^{n} x_{i l} x_{l j}, \quad i, j=1,2, \cdots, m rij=n−11l=1∑nxilxlj,i,j=1,2,⋯,m
(3) 求样本相关矩阵
R
R
R 的
k
k
k 个特征值和对应的
k
k
k 个单位特征向量。
求解
R
R
R 的特征方程
∣ R − λ I ∣ = 0 |R-\lambda I|=0 ∣R−λI∣=0
得 R R R 的 m m m 个特征值
λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ λ m \lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{m} λ1⩾λ2⩾⋯⩾λm
求方差贡献率
∑
i
=
1
k
η
i
\sum_{i=1}^{k} \eta_{i}
∑i=1kηi 达到预定值的主成分个数
k
。
k_{\text {。 }}
k。
求前
k
k
k 个特征值对应的单位特征向量
a i = ( a 1 i , a 2 i , ⋯ , a m i ) T , i = 1 , 2 , ⋯ , k a_{i}=\left(a_{1 i}, a_{2 i}, \cdots, a_{m i}\right)^{\mathrm{T}}, \quad i=1,2, \cdots, k ai=(a1i,a2i,⋯,ami)T,i=1,2,⋯,k
(4) 求
k
k
k 个样本主成分
以
k
k
k 个单位特征向量为系数进行线性变换, 求出
k
k
k 个样本主成分
y i = a i T x , i = 1 , 2 , ⋯ , k y_{i}=a_{i}^{\mathrm{T}} \boldsymbol{x}, \quad i=1,2, \cdots, k yi=aiTx,i=1,2,⋯,k
(5) 计算
k
k
k 个主成分
y
j
y_{j}
yj 与原变量
x
i
x_{i}
xi 的相关系数
ρ
(
x
i
,
y
j
)
\rho\left(x_{i}, y_{j}\right)
ρ(xi,yj), 以及
k
k
k 个主成分对原 变量
x
i
x_{i}
xi 的贡献率
ν
i
\nu_{i}
νi 。
(6) 计算
n
n
n 个样本的
k
k
k 个主成分值
将规范化样本数据代入
k
k
k 个主成分式
y
i
=
a
i
T
x
,
i
=
1
,
2
,
⋯
,
k
y_{i}=a_{i}^{\mathrm{T}} \boldsymbol{x}, \quad i=1,2, \cdots, k
yi=aiTx,i=1,2,⋯,k, 得到
n
n
n 个样本的主成分值。第
j
j
j 个 样本
x
j
=
(
x
1
j
,
x
2
j
,
⋯
,
x
m
j
)
T
\boldsymbol{x}_{j}=\left(x_{1 j}, x_{2 j}, \cdots, x_{m j}\right)^{\mathrm{T}}
xj=(x1j,x2j,⋯,xmj)T 的第
i
i
i 主成分值是
y i j = ( a 1 i , a 2 i , ⋯ , a m i ) ( x 1 j , x 2 j , ⋯ , x m j ) T = ∑ l = 1 m a l i x l j i = 1 , 2 , ⋯ , m , j = 1 , 2 , ⋯ , n \begin{gathered} y_{i j}=\left(a_{1 i}, a_{2 i}, \cdots, a_{m i}\right)\left(x_{1 j}, x_{2 j}, \cdots, x_{m j}\right)^{\mathrm{T}}=\sum_{l=1}^{m} a_{l i} x_{l j} \\ i=1,2, \cdots, m, \quad j=1,2, \cdots, n \end{gathered} yij=(a1i,a2i,⋯,ami)(x1j,x2j,⋯,xmj)T=l=1∑malixlji=1,2,⋯,m,j=1,2,⋯,n
主成分分析得到的结果可以用于其他机器学习方法的输入。比如，将样本点投影到以主成分为坐标轴的空间中，然后应用聚类算法, 就可以对样本点进行聚类。
16.2.3 数据矩阵的奇异值分解算法
给定样本矩阵
X
X
X, 利用数据矩阵奇异值分解进行主成分分析。具体过程如下。 这 里假设有
k
k
k 个主成分。
参照式
A
≈
U
k
Σ
k
V
k
T
A \approx U_{k} \Sigma_{k} V_{k}^{\mathrm{T}}
A≈UkΣkVkT, 对于
m
×
n
m \times n
m×n 实矩阵
A
A
A, 假设其秩为
r
,
0
<
k
<
r
r, 0<k<r
r,0<k<r, 则可以将矩阵
A
A
A 进行截断奇异值分解
A ≈ U k Σ k V k T A \approx U_{k} \Sigma_{k} V_{k}^{\mathrm{T}} A≈UkΣkVkT
式中
U
k
U_{k}
Uk 是
m
×
k
m \times k
m×k 矩阵,
V
k
V_{k}
Vk 是
n
×
k
n \times k
n×k 矩阵,
Σ
k
\Sigma_{k}
Σk 是
k
k
k 阶对角矩阵;
U
k
,
V
k
U_{k}, V_{k}
Uk,Vk 分别由取
A
A
A 的完全奇异值分解的矩阵
U
,
V
U, V
U,V 的前
k
k
k 列,
Σ
k
\Sigma_{k}
Σk 由取
A
A
A 的完全奇异值分解的矩阵
Σ
\Sigma
Σ 的 前
k
k
k 个对角线元素得到。
定义一个新的
n
×
m
n \times m
n×m 矩阵
X
′
X^{\prime}
X′
X ′ = 1 n − 1 X T X^{\prime}=\frac{1}{\sqrt{n-1}} X^{\mathrm{T}} X′=n−1 1XT
X ′ X^{\prime} X′ 的每一列均值为零。不难得知,
X ′ T X ′ = ( 1 n − 1 X T ) T ( 1 n − 1 X T ) = 1 n − 1 X X T \begin{aligned} X^{\prime \mathrm{T}} X^{\prime} &=\left(\frac{1}{\sqrt{n-1}} X^{\mathrm{T}}\right)^{\mathrm{T}}\left(\frac{1}{\sqrt{n-1}} X^{\mathrm{T}}\right) \\ &=\frac{1}{n-1} X X^{\mathrm{T}} \end{aligned} X′TX′=(n−1 1XT)T(n−1 1XT)=n−11XXT
即 X ′ T X ′ X^{\prime \mathrm{T}} X^{\prime} X′TX′ 等于 X X X 的协方差矩阵 S X ∘ S_{X^{\circ}} SX∘
S X = X ′ T X ′ S_{X}=X^{\prime \mathbf{T}} X^{\prime} SX=X′TX′
主成分分析归结于求协方差矩阵
S
X
S_{X}
SX 的特征值和对应的单位特征向量, 所以问题转化 为求矩阵
X
′
T
X
′
X^{\prime \mathrm{T}} X^{\prime}
X′TX′ 的特征值和对应的单位特征向量。
假设
X
′
X^{\prime}
X′ 的截断奇异值分解为
X
′
=
U
Σ
V
T
X^{\prime}=U \Sigma V^{\mathrm{T}}
X′=UΣVT, 那么
V
V
V 的列向量就是
S
X
=
X
′
T
X
′
S_{X}=X^{\prime \mathbf{T}} X^{\prime}
SX=X′TX′ 的单位特征向量。因此,
V
V
V 的列向量就是
X
X
X 的主成分。于是, 求
X
X
X 主成分可以通过 求
X
′
X^{\prime}
X′ 的奇异值分解来实现。
S x = X ′ T X ′ = ( U Σ V T ) T U Σ V T = V Σ T U T U Σ V T V T S x V = V T V Σ 2 V T V = = > V T S x V = Σ 2 \begin{aligned} S_x&=X^{\prime \mathbf{T}} X^{\prime}\\ &=(U \Sigma V^{\mathrm{T}})^TU \Sigma V^{\mathrm{T}}\\ &=V \Sigma^TU^TU\Sigma V^T\\ &V^TS_xV=V^TV\Sigma^2V^TV==>V^TS_xV=\Sigma^2 \end{aligned} Sx=X′TX′=(UΣVT)TUΣVT=VΣTUTUΣVTVTSxV=VTVΣ2VTV==>VTSxV=Σ2
具体算法如下。
算法
16.1
16.1
16.1 (主成分分析算法)
输入:
m
×
n
m \times n
m×n 样本矩阵
X
X
X, 其每一行元素的均值为零;
输出:
k
×
n
k \times n
k×n 样本主成分矩阵
Y
Y
Y 。
参数: 主成分个数
k
k
k
(1) 构造新的
n
×
m
n \times m
n×m 矩阵
X ′ = 1 n − 1 X T X^{\prime}=\frac{1}{\sqrt{n-1}} X^{\mathrm{T}} X′=n−1 1XT
X
′
X^{\prime}
X′ 每一列的均值为零。
(2) 对矩阵
X
′
X^{\prime}
X′ 进行截断奇异值分解, 得到
X ′ = U Σ V T X^{\prime}=U \Sigma V^{\mathrm{T}} X′=UΣVT
有
k
k
k 个奇异值、奇异向量。矩阵
V
V
V 的前
k
k
k 列构成
k
k
k 个样本主成分。
(3) 求
k
×
n
k \times n
k×n 样本主成分矩阵
Y = V T X Y=V^{\mathrm{T}} X Y=VTX
至此，PCA的讲解还不是太透彻，我们再来看PCA中从两个角度来推导的PCA
PRML主成分分析
有两种经常使用的PCA的定义, 它们会给出同样的算法。PCA可以被定义为数据在低维线性 空间上的正交投影, 这个线性空间被称为主子空间（principal subspace）, 使得投影数据的方差 被最大化（Hotelling, 1933）。等价地, 它也可以被定义为使得平均投影代价最小的线性投影。 平均投影代价是指数据点和它们的投影之间的平均平方距离（Pearson, 1901）。正交投影的过 程如下图所示。我们依次讨论这些定义。
图1 : 主成分分析寻找⼀个低维空间，被称为主⼦平⾯，⽤紫⾊的线表⽰，使得数据点（红点）在⼦空
间上的正交投影能够最⼤化投影点（绿点）的⽅差。 PCA的另⼀个定义基于的是投影误差的平⽅和的最
⼩值，⽤蓝线表⽰。
最⼤⽅差形式
考虑一组观测数据集
{
x
n
}
\left\{\boldsymbol{x}_{n}\right\}
{xn}, 其中
n
=
1
,
…
,
N
n=1, \ldots, N
n=1,…,N, 因此
x
n
\boldsymbol{x}_{n}
xn 是一个
D
D
D 维欧几里得空间中的变量。 我们的目标是将数据投影到维度
M
<
D
M<D
M<D 的空间中, 同时最大化投影数据的方差。现阶段, 我们 假设
M
M
M 的值是给定的。稍后在本章中, 我们会研究从数据中确定合适的
M
M
M 值的方法。
首先, 考虑在一维空间
(
M
=
1
)
(M=1)
(M=1) 上的投影。我们可以使用
D
D
D 维向量
u
1
\boldsymbol{u}_{1}
u1 定义这个空间的方 向。为了方便（并且不失一般性）, 我们假定选择一个单位向量, 从而
u
1
T
u
1
=
1
\boldsymbol{u}_{1}^{T} \boldsymbol{u}_{1}=1
u1Tu1=1 (注意, 我们 只对
u
1
\boldsymbol{u}_{1}
u1 的方向感兴趣, 而对
u
1
\boldsymbol{u}_{1}
u1 本身的大小不感兴趣
)
)
) 。这样, 每个数据点
x
n
\boldsymbol{x}_{n}
xn 被投影到一个标量 值
u
1
T
x
n
\boldsymbol{u}_{1}^{T} \boldsymbol{x}_{n}
u1Txn 上。投影数据的均值是
u
1
T
x
‾
\boldsymbol{u}_{1}^{T} \overline{\boldsymbol{x}}
u1Tx, 其中,
x
‾
\overline{\boldsymbol{x}}
x 是样本集合的均值, 形式为
x ‾ = 1 N ∑ n = 1 N x n \overline{\boldsymbol{x}}=\frac{1}{N} \sum_{n=1}^{N} \boldsymbol{x}_{n} x=N1n=1∑Nxn
投影数据的方差为
1 N ∑ n = 1 N { u 1 T x n − u 1 T x ‾ } 2 = u 1 T S u 1 \frac{1}{N} \sum_{n=1}^{N}\left\{\boldsymbol{u}_{1}^{T} \boldsymbol{x}_{n}-\boldsymbol{u}_{1}^{T} \overline{\boldsymbol{x}}\right\}^{2}=\boldsymbol{u}_{1}^{T} \boldsymbol{S} \boldsymbol{u}_{1} N1n=1∑N{u1Txn−u1Tx}2=u1TSu1
其中 S S S 是数据的协方差矩阵, 定义为
S = 1 N ∑ n = 1 N ( x n − x ‾ ) ( x n − x ‾ ) T \boldsymbol{S}=\frac{1}{N} \sum_{n=1}^{N}\left(\boldsymbol{x}_{n}-\overline{\boldsymbol{x}}\right)\left(\boldsymbol{x}_{n}-\overline{\boldsymbol{x}}\right)^{T} S=N1n=1∑N(xn−x)(xn−x)T
我们现在关于 u 1 \boldsymbol{u}_{1} u1 最大化投影方差 u 1 T S u 1 \boldsymbol{u}_{1}^{T} \boldsymbol{S} \boldsymbol{u}_{1} u1TSu1 。很明显, 最大化的过程必须满足一定的限制来防 止 ∥ u 1 ∥ → ∞ \left\|\boldsymbol{u}_{1}\right\| \rightarrow \infty ∥u1∥→∞ 。恰当的限制来自归一化条件 u 1 T u 1 = 1 \boldsymbol{u}_{1}^{T} \boldsymbol{u}_{1}=1 u1Tu1=1 。为了强制满足这个限制, 我们引入拉格 朗日乘数, 记作 λ 1 \lambda_{1} λ1, 然后对下式进行一个无限制的最大化
u 1 T S u 1 + λ 1 ( 1 − u 1 T u 1 ) \boldsymbol{u}_{1}^{T} \boldsymbol{S} \boldsymbol{u}_{1}+\lambda_{1}\left(1-\boldsymbol{u}_{1}^{T} \boldsymbol{u}_{1}\right) u1TSu1+λ1(1−u1Tu1)
通过令它关于 u 1 \boldsymbol{u}_{1} u1 的导数等于零, 我们看到驻点满足
S u 1 = λ 1 u 1 \boldsymbol{S} \boldsymbol{u}_{1}=\lambda_{1} \boldsymbol{u}_{1} Su1=λ1u1
这表明 u 1 \boldsymbol{u}_{1} u1 一定是 S \boldsymbol{S} S 的一个特征向量。如果我们左乘 u 1 T \boldsymbol{u}_{1}^{T} u1T, 使用 u 1 T u 1 = 1 \boldsymbol{u}_{1}^{T} \boldsymbol{u}_{1}=1 u1Tu1=1, 我们看到方差为
u 1 T S u 1 = λ 1 \boldsymbol{u}_{1}^{T} \boldsymbol{S} \boldsymbol{u}_{1}=\lambda_{1} u1TSu1=λ1
因此当我们将 u 1 \boldsymbol{u}_{1} u1 设置为与具有最大的特征值 λ 1 \lambda_{1} λ1 的特征向量相等时, 方差会达到最大值。这个特 征向量被称为第一主成分。
我们可以用一种增量的方式定义额外的主成分, 方法为: 在所有与那些已经考虑过的方向正 交的所有可能的方向中, 将新的方向选择为最大化投影方差的方向。如果我们考虑
M
M
M 维投影 空间的一般情形, 那么最大化投影数据方差的最优线性投影由数据协方差矩阵
S
S
S 的
M
M
M 个特征 向量
u
1
,
…
,
u
M
\boldsymbol{u}_{1}, \ldots, \boldsymbol{u}_{M}
u1,…,uM 定义, 对应于
M
M
M 个最大的特征值
λ
1
,
…
,
λ
M
\lambda_{1}, \ldots, \lambda_{M}
λ1,…,λM 。可以通过归纳法很容易地证明出 来。
总结一下, 主成分分析涉及到计算数据集的均值
x
ˉ
\bar{x}
xˉ 和协方差矩阵
S
S
S, 然后寻找
S
S
S 的对应 于
M
M
M 个最大特征值的
M
M
M 个特征向量。寻找特征值和特征向量的算法以及与特征向量分解相关的 定理, 可以参考Golub and Van Loan (1996)。注意, 计算一个
D
×
D
D \times D
D×D 矩阵的完整的特征向量 分解的代价为
O
(
D
3
)
O\left(D^{3}\right)
O(D3) 。如果我们计划将我们的数据投影到前
M
M
M 个主成分中, 那么我们只需寻 找前
M
M
M 个特征值和特征向量。这可以使用更高效的方法得到, 例如幂方法（power method） (Golub and Van Loan, 1996)，它的时间复杂度为
O
(
M
D
2
)
O\left(M D^{2}\right)
O(MD2), 或者我们也可以使用EM算法。
最小误差形式（最小重构距离）
我们现在讨论 P C A \mathrm{PCA} PCA 的另一种形式, 基于误差最小化的投影。为了完成这一点, 我们引入 D D D 维 基向量的一个完整的单位正交集合 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui}, 其中 i = 1 , … , D i=1, \ldots, D i=1,…,D, 且满足
u i T u j = δ i j \boldsymbol{u}_{i}^{T} \boldsymbol{u}_{j}=\delta_{i j} uiTuj=δij
由于基是完整的, 因此每个数据点可以精确地表示为基向量的一个线性组合, 即
x n = ∑ i = 1 D α n i u i \boldsymbol{x}_{n}=\sum_{i=1}^{D} \alpha_{n i} \boldsymbol{u}_{i} xn=i=1∑Dαniui
其中, 系数 α n i \alpha_{n i} αni 对于不同的数据点来说是不同的。这对应于将坐标系旋转到了一个由 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui} 定义的 新坐标系, 原始的 D D D 个分量 { x n 1 , … , x n D } \left\{x_{n 1}, \ldots, x_{n D}\right\} {xn1,…,xnD} 被替换为一个等价的集合 { α n 1 , … , α n D } \left\{\alpha_{n 1}, \ldots, \alpha_{n D}\right\} {αn1,…,αnD} 。与 u j \boldsymbol{u}_{j} uj 做内 积, 然后使用单位正交性质, 我们有 α n j = x n T u j \alpha_{n j}=\boldsymbol{x}_{n}^{T} \boldsymbol{u}_{j} αnj=xnTuj, 因此不失一般性, 我们有
x n = ∑ i = 1 D ( x n T u i ) u i \boldsymbol{x}_{n}=\sum_{i=1}^{D}\left(\boldsymbol{x}_{n}^{T} \boldsymbol{u}_{i}\right) \boldsymbol{u}_{i} xn=i=1∑D(xnTui)ui
然而, 我们的目标是使用限定数量 M < D M<D M<D 个变量的一种表示方法来近似数据点, 这对应于 在低维子空间上的一个投影。不失一般性, M M M 维线性子空间可以用前 M M M 个基向量表示, 因此我 们可以用下式来近似每个数据点 x n \boldsymbol{x}_{n} xn
x ~ n = ∑ i = 1 M z n i u i + ∑ i = M + 1 D b i u i \tilde{\boldsymbol{x}}_{n}=\sum_{i=1}^{M} z_{n i} \boldsymbol{u}_{i}+\sum_{i=M+1}^{D} b_{i} \boldsymbol{u}_{i} x~n=i=1∑Mzniui+i=M+1∑Dbiui
其中 { z n i } \left\{z_{n i}\right\} {zni} 依赖于特定的数据点, 而 { b i } \left\{b_{i}\right\} {bi} 是常数, 对于所有数据点都相同。我们可以任意选 择 { u i } , { z n i } \left\{\boldsymbol{u}_{i}\right\},\left\{z_{n i}\right\} {ui},{zni} 和 { b i } \left\{b_{i}\right\} {bi}, 从而最小化由维度降低所引入的失真。作为失真的度量, 我们使用原始数 据点与它的近似点 x ~ n \tilde{x}_{n} x~n 之间的平方距离, 在数据集上取平均。因此我们的目标是最小化
J = 1 N ∑ n = 1 N ∥ x n − x ~ n ∥ 2 J=\frac{1}{N} \sum_{n=1}^{N}\left\|\boldsymbol{x}_{n}-\tilde{\boldsymbol{x}}_{n}\right\|^{2} J=N1n=1∑N∥xn−x~n∥2
首先考虑关于 { z n i } \left\{z_{n i}\right\} {zni} 的最小化。消去 x ~ n \tilde{\boldsymbol{x}}_{n} x~n, 令它关于 z n j z_{n j} znj 的导数为零, 然后使用单位正交的条 件, 我们有
z n j = x n T u j z_{n j}=\boldsymbol{x}_{n}^{T} \boldsymbol{u}_{j} znj=xnTuj
其中 j = 1 , … , M j=1, \ldots, M j=1,…,M 。类似地, 令 J J J 关于 b i b_{i} bi 的导数等于零, 再次使用单位正交的关系, 我们有
b j = x ‾ T u j b_{j}=\overline{\boldsymbol{x}}^{T} \boldsymbol{u}_{j} bj=xTuj
其中 j = M + 1 , … , D j=M+1, \ldots, D j=M+1,…,D 。如果我们消去 x ~ n = ∑ i = 1 M z n i u i + ∑ i = M + 1 D b i u i \tilde{\boldsymbol{x}}_{n}=\sum_{i=1}^{M} z_{n i} \boldsymbol{u}_{i}+\sum_{i=M+1}^{D} b_{i} \boldsymbol{u}_{i} x~n=∑i=1Mzniui+∑i=M+1Dbiui中的 z n i z_{n i} zni 和 b i b_{i} bi, 使用一般的展开式 x n = ∑ i = 1 D ( x n T u i ) u i \boldsymbol{x}_{n}=\sum_{i=1}^{D}\left(\boldsymbol{x}_{n}^{T} \boldsymbol{u}_{i}\right) \boldsymbol{u}_{i} xn=∑i=1D(xnTui)ui, 我们 有
x n − x ~ n = ∑ i = M + 1 D { ( x n − x ‾ ) T u i } u i \boldsymbol{x}_{n}-\tilde{\boldsymbol{x}}_{n}=\sum_{i=M+1}^{D}\left\{\left(\boldsymbol{x}_{n}-\overline{\boldsymbol{x}}\right)^{T} \boldsymbol{u}_{i}\right\} \boldsymbol{u}_{i} xn−x~n=i=M+1∑D{(xn−x)Tui}ui
从中我们看到, 从 x n \boldsymbol{x}_{n} xn 到 x ~ n \tilde{\boldsymbol{x}}_{n} x~n 的位移向量位于与主子空间垂直的空间中, 因为它是 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui} 的线性组 合, 其中 i = M + 1 , … , D i=M+1, \ldots, D i=M+1,…,D, 如图 1所示。这与预期相符, 因为投影点 x ~ n \tilde{x}_{n} x~n 一定位于主子空间 内, 但是我们可以在那个子空间内自由移动投影点, 因此最小的误差由正交投影给出。 于是, 我们得到了失真度量 J J J 的表达式, 它是一个纯粹的关于 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui} 的函数, 形式为
J = 1 N ∑ n = 1 N ∑ i = M + 1 D ( x n T u i − x ‾ T u i ) 2 = ∑ i = M + 1 D u i T S u i J=\frac{1}{N} \sum_{n=1}^{N} \sum_{i=M+1}^{D}\left(\boldsymbol{x}_{n}^{T} \boldsymbol{u}_{i}-\overline{\boldsymbol{x}}^{T} \boldsymbol{u}_{i}\right)^{2}=\sum_{i=M+1}^{D} \boldsymbol{u}_{i}^{T} \boldsymbol{S} \boldsymbol{u}_{i} J=N1n=1∑Ni=M+1∑D(xnTui−xTui)2=i=M+1∑DuiTSui
剩下的任务是关于 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui} 对 J J J 进行最小化, 这必须是具有限制条件的最小化, 因为如果不这 样, 我们会得到 u i = 0 \boldsymbol{u}_{i}=0 ui=0 这一没有意义的结果。限制来自于单位正交条件, 并且正如我们将看到 的那样, 解可以表示为协方差矩阵的特征向量展开式。在考虑一个形式化的解之前, 让我们试 着直观地考察一下这个结果。考虑二维数据空间 D = 2 D=2 D=2 以及一维主子空间 M = 1 M=1 M=1 的情形。我们必 须选择一个方向 u 2 \boldsymbol{u}_{2} u2 来最小化 J = u 2 T S u 2 J=\boldsymbol{u}_{2}^{T} \boldsymbol{S} \boldsymbol{u}_{2} J=u2TSu2, 同时满足限制条件 u 2 T u 2 = 1 \boldsymbol{u}_{2}^{T} \boldsymbol{u}_{2}=1 u2Tu2=1 。使用拉格朗日乘数 λ 2 \lambda_{2} λ2 来 强制满足这个限制, 我们考虑最小化
J ~ = u 2 T S u 2 + λ 2 ( 1 − u 2 T u 2 ) \tilde{J}=\boldsymbol{u}_{2}^{T} \boldsymbol{S} \boldsymbol{u}_{2}+\lambda_{2}\left(1-\boldsymbol{u}_{2}^{T} \boldsymbol{u}_{2}\right) J~=u2TSu2+λ2(1−u2Tu2)
令关于
u
2
\boldsymbol{u}_{2}
u2 的导数等于零, 我们有
S
u
2
=
λ
2
u
2
\boldsymbol{S} \boldsymbol{u}_{2}=\lambda_{2} \boldsymbol{u}_{2}
Su2=λ2u2, 从而
u
2
\boldsymbol{u}_{2}
u2 是
S
\boldsymbol{S}
S 的一个特征向量, 且特征值为
λ
2
\lambda_{2}
λ2 。因 此任何特征向量都会定义失真度量的一个驻点。为了找到
J
J
J 在最小值点处的值, 我们将
u
2
\boldsymbol{u}_{2}
u2 的解 代回到失真度量中, 得到
J
=
λ
2
J=\lambda_{2}
J=λ2 。于是, 我们通过将
u
2
\boldsymbol{u}_{2}
u2 选择为对应于两个特征值中较小的那个 特征值的特征向量, 可以得到
J
J
J 的最小值。因此, 我们应该将主子空间与具有较大的特征值的 特征向量对齐。这个结果与我们的直觉相符, 即为了最小化平均平方投影距离, 我们应该将主 成分子空间选为穿过数据点的均值并且与最大方差的方向对齐。对于特征值相等的情形, 任何 主方向的选择都会得到同样的
J
J
J 值。
对于任意的
D
D
D 和任意的
M
<
D
M<D
M<D, 最小化
J
J
J 的一般解都可以通过将
{
u
i
}
\left\{\boldsymbol{u}_{i}\right\}
{ui} 选择为协方差矩阵的特 征向量的方式得到, 即
S u i = λ i u i \boldsymbol{S} \boldsymbol{u}_{i}=\lambda_{i} \boldsymbol{u}_{i} Sui=λiui
其中 i = 1 , … , D i=1, \ldots, D i=1,…,D, 并且与平常一样, 特征向量 { u i } \left\{\boldsymbol{u}_{i}\right\} {ui} 被选为单位正交的。失真度量的对应的值为
J = ∑ i = M + 1 D λ i J=\sum_{i=M+1}^{D} \lambda_{i} J=i=M+1∑Dλi
这就是与主子空间正交的特征值的加和。于是, 我们可以通过将这些特征向量选择成 D − M D-M D−M 个 最小的特征值对应的特征向量, 来得到 J J J 的最小值, 因此定义了主子空间的特征向量是对应 于 M M M 个最大特征值的特征向量。
虽然我们已经考虑了
M
<
D
M<D
M<D 的情形, 但是
P
C
A
\mathrm{PCA}
PCA 对于
M
=
D
M=D
M=D 的情形仍然成立, 这种情况下没 有维度的降低, 仅仅是将坐标轴旋转, 与主成分对齐即可。
最后, 值得注意的时, 存在一个与此密切相关的线性维度降低的方法, 被称为典型相关分析 (canonical correlation analysis)，或者CCA（Hotelling, 1936; Bach and Jordan, 2002）。PCA操 作的对象是一个随机变量, 而CCA考虑两个（或者更多）的变量, 并且试图找到具有较高的交 叉相关性的线性子空间对, 从而在一个子空间中的每个分量都与另一个子空间的一个分量具有 相关性。它的解可以表示为一般的特征向量问题。
PCA的应用
-
数据压缩
-
数据预处理
- 对数据的白化或者球形化
-
数据可视化
主成分分析的另一个常见应用是数据可视化。这里, 每个数据点被投影到二维 ( M = 2 ) (M=2) (M=2) 的 主子空间中, 从而数据点 x n \boldsymbol{x}_{n} xn 被画在了一个笛卡尔坐标系中, 坐标系由 x n T u 1 \boldsymbol{x}_{n}^{T} \boldsymbol{u}_{1} xnTu1 和 x n T u 2 \boldsymbol{x}_{n}^{T} \boldsymbol{u}_{2} xnTu2 定义, 其 中 u 1 u_{1} u1 和 u 2 u_{2} u2 是特征向量, 对应于最大的和第二大的特征值。对于石油流数据集, 这种图的一个例子如下图所示。
标签：李航,right,16,--,boldsymbol,cdots,alpha,left,mathrm 来源： https://blog.csdn.net/wl1780852311/article/details/121184742