主成分分析（PCA）学习笔记

学生而已

已于 2023-10-18 19:44:54 修改

阅读量516

点赞数 4

分类专栏：主成分分析（PCA）学习文章标签：机器学习人工智能

于 2023-06-30 16:53:58 首次发布

本文链接：https://blog.csdn.net/vv0610_/article/details/131472696

版权

主成分分析（PCA）学习专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、前言

二、简介

三、第一个问题，降维的过程是怎么实现的？

四、第二个问题，怎么知道 d 个成分就是主要成分？

十、求解PCA（一）：统计学原理（拉格朗日乘子法求解）

十一、求解PCA（二）：几何学原理（拉格朗日乘子法求解）

十二、求解PCA（三）：梯度下降法（或上升法）求解

十三、后记

十四、推荐阅读

一、前言

历经了多日的折磨，PCA总算是学明白了一些。我一直坚信学懂的评判标准是能清晰地把知识点讲出来，故作此文。另一方面，也希望保留一些笔记，方便自己查阅。

在开始读这篇文章之前，我需要提醒你，这篇文章中会用到高等数学、线性代数、概率论与统计的多个知识点，如果有实在看不懂的地方，请先查阅相关资料。

最后，仅此一篇文章是不会覆盖到PCA的方方面面的，但此文应该可以帮助小白快速入门。

二、简介

主成分分析（Principal Component Analysis，PCA），我并不想在一开始就聊一些太复杂的概念，这样对于初学者来说是地狱一般的体验。所以，我们举个形象点的例子：

假如我们有一组多元随机变量（Multivariate Random Variable）：

$\boldsymbol{x} = \begin{Bmatrix} x_{1} & ... & x_{D} \end{Bmatrix}^{T}\in \mathbb{R}^{D}$

请注意，x 是一个 D 维列向量，后文我们将用加粗的 x 表示向量（Vector），不加粗的 x 表示标量（Scalar）。ps. 我真要吐槽 csdn 的文本编辑器，公式和文字在一行就会对不齐。

很显然，x 中有 D 个标量，我们把这些 xi 叫做 x 的成分，这可以表示很多含义，例如 x 是某个学生的考试成绩，其中 x1 是数学成绩，x2 是英语成绩，以此类推。假如说，我们要从 x 的 D 个成分中，提取出 d 个成分（其中 d << D），我们该如何操作呢？假如我们提取出的 d 个成分构成一个新的多元随机变量 y，那么我们可以这样表示它：

$\boldsymbol{y} = \begin{Bmatrix} y_{1} & ... & y_{d} \end{Bmatrix}^{T}\in \mathbb{R}^{d}$

于是我们会发现，从 x 到 y，数据的维度从 D 变成了 d，维度减少了。那么如果我们提取出的 d 个成分是我们想要的 x 的 d 个主要成分（目前暂时可以理解成对数据影响比较大的关键成分），那么这一过程就是主成分分析。所以不难理解，PCA 是一个降维的过程。

那么接下来的问题就是，如何实现从 x 到 y 的这一过程，而在这一过程中，我们会涉及到很多要解决的问题。

三、第一个问题，降维的过程是怎么实现的？

思考一下，前面的数据从 D 维变到了 d 维，这其中有什么玄机？没错，上文提到的 x 其实只是单点数据，或者叫一个数据。我们还拿上面的成绩单来举个例子，x 只是其中某个人的各科成绩。而这一个 x，就是在 D 维空间中的一个“点”，如果要画图的话，其各维度分量其实就是对应的坐标（比如有三门成绩语文数学英语，可以画出三个坐标轴分别代表语文数学英语，每门课的分数就是对应坐标轴上的坐标）。其实单独考虑一个数据点是没有太大意义的，在实际 PCA 的应用过程中，我们有成千上万个这样的数据点，这些 x 共同构成了一个数据矩阵（数据集） X，假设一共有 N 个数据，我们就可以这么写：

$X = \begin{Bmatrix} \mathbf{X}_{1} & ... & \mathbf{X}_{N} \end{Bmatrix}\in \mathbb{R}^{D\times N}$

数据矩阵 X 中的每一个 X_i，其实都是一个 D 维的列向量。那么这时再选取主成分之后会变成什么样呢？见下：

$Y = \begin{Bmatrix} \mathbf{Y}_{1} & ... & \mathbf{Y}_{N} \end{Bmatrix}\in \mathbb{R}^{d\times N}$

不难发现，选取主成分并没有改变数据的个数，而是降低了数据的维度。这样的话，我们就需要一个统一的标准，来对数据矩阵 X 中的每个数据 X_i 进行同样的降维操作（为了方便讲解，我们后面还是用 x 来表示某个不具体的 X_i）。而对于这一操作，我们是这样进行的：

$\boldsymbol{y} = U^{T}\boldsymbol{x} ~~~~~~~~~~ s.t.~~~~~ U^{T}U=I_{d},~~U\in \mathbb{R}^{D\times d}$

跳跃有点大？没关系，我们来分析一下其中奥妙。来看一下其中各个字母表示的含义：首先 s.t. 是 subject to 的缩写，也就是前面的式子服从于后面的条件。其次，第一个约束条件中的 I 是 d 阶单位矩阵（Identity Matrix）。然后，从第一个约束条件可以看出 U 是一个标准正交矩阵（Orthonormal Matrix）。最后，U^T 是 U 的转置（Transpose）。好了，到此我们可以开始分析。U^T 毫无疑问是一个 dxD 的矩阵，其中每一个列向量的维度都是 d，它施加在 x 上会有什么效果呢？其实说来简单，dxD 的矩阵乘以 Dx1 的列向量，结果是 dx1，你会发现维度降低了……我知道这可能没有说服力，但是现阶段这是最好的理解方式了。可能有人会问，为什么要用转置乘以 x 呢？直接用 U 乘以 x 不行吗？答案是可以，这时候你只需把 U 的维度改成 dxD 使得矩阵乘法公式成立，并且在后面的一系列推导中用 U 代替 U^T 即可。

至此，我们就有了一个简单的“降维标准”，在这里我们称之为投影会准确得多。

四、第二个问题，怎么知道 d 个成分就是主要成分？

很合理的问题。你说你选 d 个成分就是主要成分了？凭什么？

这里我们就不得不提到 PCA 的两个要求：最大可分性和最近重构性。

要解释最大可分性比较容易一些，毕竟降维是将数据点的维度减小了，所以很有可能出现多个数据点靠的很近甚至重合的情况，这时，我们就希望数据投影后的离散程度越大越好。而表示离散程度的方法，就是方差（Variance），某个多元随机变量 x 的方差记为 Var(x)。所以，最大可分性就是要求：

$\max \mathrm{Var}(\boldsymbol{x})~~~~~(4.1)$

再说说最近重构性，这个就比较麻烦一些了。首先我们要明白空间中的“近”是怎么判断的？这里我们就要提到大名鼎鼎的欧氏距离（Euclidean Distance）ps. 是欧几里得不是欧拉… 欧氏距离怎么说的来着？如果在平面上有两点 a = (x1, y1) 和 b = (x2, y2) 那么它们之间的距离定义为：

$||\boldsymbol{a}-\boldsymbol{b}||_{2}=\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}$

其实欧氏距离也是经典的 2-范数（norm），对于一般情况下，多元随机变量的 p-范数我们定义如下：

$||\boldsymbol{x}||_{p}=\sqrt[p]{|x_{1}|^{p}+|x_{2}|^{p}+...+|x_{n}|^{p}}=(\sum_{i=1}^{n}|x_{i}|^{p})^{\frac{1}{p}}$

这里稍微多说一些哈，这样后面就能省点力气了~我们会发现 p-范数是针对多元随机变量来说的，也就是一个列向量，或者 n 维空间中的点。那么对于矩阵情况又如何呢？

在矩阵范数中，有一个佛罗贝尼乌斯范数（Frobenius Norm），简称叫 F-范数，其定义式是长这样的：

$||A||_{F}= \sqrt{\sum_{i}^{m}\sum_{j}^{n}|a_{ij}|^{2}}$

其中 A 是一个 mxn 的矩阵，其实 F-范数就是求矩阵中所有元素的平方和然后开根号，这和向量的 2-范数是类似的。方便起见，在后文中，如果我们不加说明，||A|| 默认表示矩阵的 F-范数。

至此我们明白了空间中的距离是如何度量的，那么最近重构性就可以进一步解释了。“最近”也就是要求欧氏距离最近，而重构则是一个新的概念。还记得这个式子吗？

$\boldsymbol{y} = U^{T}\boldsymbol{x} ~~~~~~~~~~~~~ s.t.~~~ U^{T}U=I_{d},~~~U\in \mathbb{R}^{D\times d}$

我们前面说过，这个过程叫做投影。那重构是什么呢？我们再来写一下：

$\boldsymbol{x^{*}} = UU^{T}\boldsymbol{x} ~~~~~~~~~~~~~ s.t.~~~ U^{T}U=I_{d},~~~U\in \mathbb{R}^{D\times d}$

众所周知，如果 U 是一个标准正交矩阵，那么 U^T = U^(-1) （前提，如果有的话，U 的转置和逆矩阵（Inverse Matrix）是一样的。所以如果把 U^T 视作某种投影函数的话，U 就是这个过程的逆过程。即：U^T 对 x 进行降维得到 y，U 对 y 升维得到 x*。这便是重构的过程。至此不免有些疑问，既然 U^T = U^(-1)，那么 UU^T 便等于 I（单位阵），显然 x*=x 呀。话虽如此，对于一个数据点我们当然可以这么理解，但是如果有多个数据点，降维的过程中可能会发生信息的损失（这个说来话长了，简单理解为不同数据点重叠，没法完全恢复即可），那么重构的过程必然会伴随着误差，定义为 ||x-x*||。这便是欧氏距离的定义了。没错，最近重构性指的就是重构之后的数据点距离原数据点最近，即：

$\min ||\boldsymbol{x}-UU^{T}\boldsymbol{x}||~~~~~(4.2)$

现在让我们来回答最初的问题，如何说明 d 个成分就是主要成分呢？注意到 (4.1) 式和 (4.2) 式中只有 U 是变量。所以，如果能找到一个 U，同时满足 (4.1) 和 (4.2) 式，那么 U 便是我们想要的投影矩阵，经过 U^T 的作用我们得到 y，而 y 中的分量便是 x 的 d 个主要成分。这里我们可以直接丢出最终的结论了：

$U^{*}=\arg \max _{U} \mathrm{Var}(U^{T}\boldsymbol{x})=\arg \min _{U} ||\boldsymbol{x}-UU^{T}\boldsymbol{x}||$

我们距离求解这个问题还很远很远，也许你还会有很多问题，不要急躁，我们一点点来。

五、看一个例子：二维世界

这个例子一定是简单到不能再简单的：在二维世界中，我们每组数据的维度只有两维，也就是 x=[x1, x2]^T。假设现在有四个点，x1,x2,x3,x4。如下图：

显然他们的均值为0（这里是指，x_11+x_21+x_31+x_41=0 且 x_12+x_22+x_32+x_42=0），过原点画一条直线（想象一下二维空间任意一个向量都可以表示成过原点的有向线段，这个就是方向向量，然后过某一方向向量画一条直线）：

假设这个黄色直线的方向向量是 u_1。现在，将四个点投影到黄色直线上。

可以发现，四条蓝色线段（代表投影的方向）其实都是沿着同一方向相同或相反的，我们可以把它们都平移到原点，并且将方向向量命名为 u_2 （正负不重要）。

很显然我们可以知道 u_1⊥u_2，这就意味着 u_1^Tu_2 = 0，而 u_1 和 u_2 又都是单位向量。所以 u_1^Tu_1 = u_2^Tu_2 = 1。我们现在单独来看 x_1。

我们来回忆一下前面说过的最大可分性和最近重构性。最大可分性指的是数据距离样本中心（均值）的距离最小，而可分性在图中就是投影点到原点的距离，正好可以对应红线。最近重构性指的是投影点距离原数据点的距离最小，而重构性在途中就是投影点到原数据点的距离，正好对应蓝线。紫线就是数据点自己本身决定的，也就是数据点到原点的欧氏距离。我们会发现，蓝色线长度^2+红色线长度^2 = 紫色线长度^2。而紫色线长度是固定的，所以红色线长度越大，蓝色线长度越小。翻译一下就是：投影点方差越大，投影点距离原数据点距离越小。所以，最大可分性和最近重构性是等价的，解决了一个，就解决了另一个。这也为我们之后要说的，PCA的三种形式有着密切关系。

六、需要的工具：期望与方差

跨度确实有点大，从线性代数跑到概率论与统计来了（也许是简单随机抽样）。

我们先明确一组概念。我们可以用样本均值（Sample Mean）来估计总体均值（Population Mean），而总体均值和期望（Expected Value）是一样的。上述一串话可以这么表示（假设样本有N个）：

$E[\boldsymbol{x}]=\frac{1}{N}\sum_{i=1}^{N}\boldsymbol{x_{i}}=\bar{\boldsymbol{x}}$

注意区别，我们用 E[x] 或 μ 表示期望，用 x-bar 表示均值。

有了期望之后，总体方差（Population Variance）样本方差（Sample Variance）也就有了：

$\mathrm{Var}(\boldsymbol{x})=E[(\boldsymbol{x}-E(\boldsymbol{x}))^{2}]=E[(\boldsymbol{x}-E(\boldsymbol{x}))(\boldsymbol{x}-E(\boldsymbol{x}))^T]$

$\sigma ^{2}=\frac{1}{N}\sum_{i=1}^{N}(x_{i}-\mu )^{2}$

$S^{2}=\frac{1}{N-1}\sum_{i=1}^{N}(x_{i}-\bar{\boldsymbol{x}})$

需要注意，求期望内部的平方展开后是外积（请参考期望定义）。另外，我们用 Var(x) 或 σ^2 表示总体方差，用 S^2 表示样本方差。

你可能会问为什么 μ 变成 x-bar 后分母变成了 N-1。好问题，请参阅 Bessel's Correction

关于正儿八经地解读，在这里我只能浅薄地引用一下，总体均值已知（μ 已知）时，我们用 N 做分母。但是当总体均值未知，而我们用样本均值去估计总体均值，会导致自由度减一（原本有 N 个独立样本，现在每个都混进了一个 x-bar，自由度变少了）。但是如果换种说话是不是就能理解了呢？当分母使用 N-1 时，S^2 是 σ^2 的无偏估计量（Unbiased Estimator）。N-1 就是用来修正用 N 做分母时的偏移的。而实际上这么修正之后 S^2 的 MSE 会变大，因为开方操作是个凹函数。当然，这不是今天要讨论的问题。

这和我们之前提到的有什么关系呢？请放心！我们很快就会用到了！

七、需要的工具：协方差矩阵

我们还是先从协方差讲起吧。假如有两个多元随机变量 x 和 y。它们的协方差（Covariance）定义为：

$cov(\boldsymbol{x},\boldsymbol{y})=E[(\boldsymbol{x}-E(\boldsymbol{x})(\boldsymbol{y}-E(\boldsymbol{y})]$

不难发现，如果尝试求 x 自己的协方差其实得到的就是 x 的方差：

$cov(\boldsymbol{x},\boldsymbol{x})=E[(\boldsymbol{x}-E(\boldsymbol{x})^{2}]=\mathrm{Var}(\boldsymbol{x})$

那么协方差矩阵长什么样呢？它有如下定义：

$K_{ij}=cov(x_i,x_j)=E[(x_i-E(\boldsymbol{x})(x_j-E(\boldsymbol{x})]$

也就是说协方差矩阵 K 的 (i,j) 位置上的数是 xi 和 xj 的协方差。

我们把前面的例子搬过来说，假如 x 是这样的：

$\boldsymbol{x} = \begin{Bmatrix} x_{1} & ... & x_{D} \end{Bmatrix}^{T}\in \mathbb{R}^{D}$

这个时候请允许我再做一个假设，即 x-bar 也就是均值为 0。实际上这个假设并不是什么太难达成的事情，如果我原本的 x 均值不为 0，只需要令其中的每个分量减去 x-bar，就得到了一个新的均值为 0 （zero-mean）的 x。如果 x 此时是对于某个总体抽样了 D 次得出的一组样本，并且均值为 0，那么关于 x_i 和 x_j 的协方差我们可以很轻松地写出来（注意是标量）：

$cov(x_{i},x_{j})=E[x_{i}x_{j}]$

从这里你能看出协方差矩阵的大小吗？cov(x_1,x_1), cov(x_1,x_2) ... cov(x1,x_D)，我想这个矩阵应该是DxD的，我们尝试画一下：

$K=\begin{bmatrix} E[x_1x_1]&E[x_1x_2] & ... & E[x_1x_D]\\ E[x_2x_1]&E[x_2x_2] & ... & E[x_2x_D] \\ ... & ... &... &... \\ E[x_Dx_1]&E[x_Dx_2] & ... & E[x_Dx_D] \end{bmatrix}$

其实通过前面的介绍，如果把 E 看成一个函数，会发现它其实是线性的，只是给括号里的东西乘了个系数，所以可以提出去：

$K=E[\begin{bmatrix} x_1x_1&x_1x_2 & ... & x_1x_D\\ x_2x_1&x_2x_2 & ... & x_2x_D \\ ... & ... &... &... \\ x_Dx_1&x_Dx_2 & ... & x_Dx_D \end{bmatrix}]$

我们仔细观察内部，不难发现每一列都可以提出一个 x，每一行都可以提出一个 x^T。然后，让我给 K 改个名：

$\sum _x = E[\boldsymbol{xx}^T]$

Σx代表 x 的协方差矩阵，而右边，就是上面推断而来的结果！这个公式十分重要，请牢记！

前面说过 E 其实只是个系数，那这个系数是多少呢？统计学上我们一般用 N-1，但是也有用 N 的。我没法确定，但是我可以提供个建议，当 N 为大样本的时候，其实没差啦。

我们刚才讨论的 x 是一个多元随机变量，是个列向量。现在我们把数据矩阵 X 搬过来，看看会发生什么。直接套公式~

$XX^T=\begin{bmatrix} \boldsymbol{X_1} & ... & \boldsymbol{X_N} \end{bmatrix}\begin{bmatrix} \boldsymbol{X_1}^T\\ ...\\ \boldsymbol{X_N}^T \end{bmatrix}=\boldsymbol{X_1X_1}^T+...+\boldsymbol{X_NX_N}^T$

我们会发现，所有数据的协方差矩阵只是将所有单组数据的协方差矩阵加在一块了。我们尝试验证一下，比如第 i 行和第 j 行的协方差应该是多少？只需要找到 XX^T 的 (i,j) 位置即可，它是将所有 X_iX_i^T_(i,j) 加在一块了，结果仍然是 i 行和 j 行的协方差。

顺带一提，我们会发现协方差矩阵对角线元素上都是方差。

以防有人混淆，我这里稍微说点。前面我们用到的 x 和 X_i 是一个概念，即一列数据。所以，x_i 是一个标量，而 X_i 是一个列向量。那么 x_i 应该对应谁呢？应该对应某一数据列中的第 i 行（或者叫第 i 个属性）。

八、需要的工具：矩阵的迹

矩阵的迹（Trace）指的是矩阵主对角线上所有元素之和，用公式来表示即：

$tr(A)=\sum_{i=1}^{n}a_{ii}=a_{11}+a_{22}+...+a_{nn}$

这里主要说几个性质就行：

第一个：

$tr(AB)=tr(BA)$

这个性质很好证，回想一下对角线上元素怎么来的即可，这里就不证了。

第二个：

$\boldsymbol{x}^TA\boldsymbol{x}=tr(A\boldsymbol{xx}^T)$

首先我们要知道 x^TAx 是一个常数，所以 x^TAx = tr(x^TAx)，再运用第一个性质，就得来了这个式子。

第三个：

$\frac{\partial }{\partial \boldsymbol{X}}tr(\boldsymbol{X}\boldsymbol{X}^T\boldsymbol{B})=\boldsymbol{B}\boldsymbol{X}+\boldsymbol{B}^T\boldsymbol{X}$

关于矩阵求导公式，来自于 The Matrix Cookbook 一书，不做深究。

九、正式开始：PCA三种形式的等价性

前面的一切准备就是为了这里，首先明确一下我们想要求解的式子是什么？前面说过最大可分性和最小重构性，分别对应投影点方差最大和投影点距离最近。而我们要求的则是这个最大值和最小值之中的矩阵 U。我们把第四部分中最后一个式子修改一下，把多元随机变量 x 修改为数据矩阵 X。

$U^{*}=\arg \max _{U} \sum_{i}^{}\mathrm{Var}(U^{T}\boldsymbol{X_i})=\arg \min _{U}\sum_{i}^{} ||\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i}|| \\s.t. U^TU=I_d$

前面我们用例子分析了第二个等号左右两边是等价的，但是现在还需要证一下：

$\arg \min _{U} \sum_{i}^{}||\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i}||^2\\ =\arg \min _{U}\sum_{i}^{} (\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i})^T(\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i})$

这是第一步，矩阵所有元素的平方和等于矩阵自己和自己的内积。

这是第二步，转置放进括号，需要改变矩阵乘法运算顺序，但是对减法没有影响。

这是第三步，需要注意 X^T 是 D 维的，提出来之后的单位矩阵是 D 阶。

第四步，这里用到了约束条件 U^TU=Id。

$\arg \min _{U} \sum_{i}^{}||\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i}||^2\\ =\arg \min _{U}\sum_{i}^{} (\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i})^T(\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i})\\ =\arg \min _{U}\sum_{i}^{}(\boldsymbol{X_i}^T-\boldsymbol{X_i}^TUU^T)(\boldsymbol{X_i}-UU^{T}\boldsymbol{X_i})\\ =\arg \min _{U}\sum_{i}^{}\boldsymbol{X_i}^T(I_D-UU^T)(I_D-UU^T)\boldsymbol{X_i}\\ =\arg \min _{U} \sum_{i}^{}\boldsymbol{X_i}^T(I_D-2UU^T+UU^TUU^T)\boldsymbol{X_i}\\ =\arg \min _{U} \sum_{i}^{}\boldsymbol{X_i}^T(I_D-UU^T)\boldsymbol{X_i}\\ =\arg \min _{U} \sum_{i}^{}(\boldsymbol{X_i}^T\boldsymbol{X_i}-\boldsymbol{X_i}^TUU^T\boldsymbol{X_i})\\ =\arg \max _{U} \sum_{i}^{}\boldsymbol{X_i}^TUU^T\boldsymbol{X_i}$

第五步，把 X_i^T 和 X_i 乘进去之后，由于 X_i^TX_i 是数据本身决定的，不会随着 U 变动，因此把求最小值转换成了求最大值的问题，最后一个式子是不是有点眼熟？我们再写一步：

最后三行你需要知道，首先，因为样本 X 均值为 0，所以：

$\mathrm{Var}(\boldsymbol{x})=E[(\boldsymbol{x}-E(\boldsymbol{x}))^{2}]=E[\boldsymbol{x}^2]$

再者，用到一些坐标变换的知识：

$\boldsymbol{y} = U^{T}\boldsymbol{x}$

||y||=||x||，因为 U^T 只将 x 旋转了，并没有改变 x 到原点的距离（也就是 x 的模长）。

最后，E 是线性的，还记得吗？所以如果 x* 能令 f(x) 取最大值，同样也可以令 cf(x)+b 取最大值。

结合以上三点，便有了最后三个等于号。至此，我们证明出了最大方差形式和最近重建形式的等价性！我们前面学到的知识都用上了~

还有第三个形式，从第五步开始：

至此，第三形式就出来了，这个形式求解 U 的话会比较方便一些。

十、求解PCA（一）：统计学原理（拉格朗日乘子法求解）

历史上，PCA首先被用于统计学中估计一组多元随机变量 x 的主要成分。如果给定一个均值为 0 的多元随机变量 x∈R^D 以及一个整数 d < D。x 的 d 个主要成分 y∈R^d 是 x 的 d 个不相关的线性成分（译自 Generalized Principal Component Analysis 一书）。

在文章开头我们有这样一个式子，其实就是对上面一段话的描述：

$\boldsymbol{y} = U^{T}\boldsymbol{x} ~~~~~~~~~~ s.t.~~~~~ U^{T}U=I_{d},~~U\in \mathbb{R}^{D\times d}$

但是这个式子处理起来略微有些复杂，所以我们把它拆开来看。U 毫无疑问是一个 Dxd 的矩阵，其中每一列我们用 u_i 来表示，于是我们可以把上式改写一下：

$y_i = \boldsymbol{u_{i}}^{T}\boldsymbol{x}\in \mathbb{R},~~~~~\boldsymbol{u_{i}}\in \mathbb{R}^D,~~~~~i=1,2,...,d$

不难发现，y_i 首先是一个标量，其次会有 d 个这样的标量，而 y_i 正是 x 的 d 个主要成分。另外，y_1 是第一主成分，y_2是第二主成分，以此类推。而我们区分第 i 主成分的关键性证据就是，Var(y_i) >= Var(y_i+1），我们希望越靠前的主成分上数据的离散程度越高。也因此，对上式我们是有约束条件的：

$\boldsymbol{u_{i}}^T\boldsymbol{u_{i}}=1~~~~~and~~~~~\mathrm{Var}(y_1)\geq \mathrm{Var}(y_2)\geq...\geq\mathrm{Var}(y_d)> 0$

如果我们将 Var(y_1) 视为目标函数来求解 u_1，其应当在约束条件下取得最大值，如下：

$\boldsymbol{u_{1}}^{*}=\arg \max \mathrm{Var}(\boldsymbol{u_{1}}^{T}\boldsymbol{x})~~~s.t.~~~\boldsymbol{u_{1}}^{T}\boldsymbol{u_{1}}=1$

这个式子需要稍微变一下才好求解，这里需要回忆一下第七部分中推出的公式：

$\sum _x = E[\boldsymbol{xx}^T]$

以及求期望是一个线性的函数。那么，不失一般性，对任何一个 u（其实是 u_i）有：

$\mathrm{Var}(\boldsymbol{u}^T\boldsymbol{x})=E[(\boldsymbol{u}^T\boldsymbol{x})^2]=E[\boldsymbol{u}^T\boldsymbol{x}\boldsymbol{x}^T\boldsymbol{u}]=\boldsymbol{u}^T\sum_{x}^{}\boldsymbol{u}$

于是前面的目标函数就等价于：

$\boldsymbol{u_{1}}^{*}=\arg \max \boldsymbol{u_1}^T\sum_{x}^{}\boldsymbol{u_1}~~~s.t.~~~\boldsymbol{u_{1}}^{T}\boldsymbol{u_{1}}=1$

然后，我们便用到了拉格朗日乘子法（Lagrange Multiplier），这个方法在高数里面讲过，我们这里就不详细展开说了。总之，构造的拉格朗日函数为：

$L(\boldsymbol{u_1},\lambda _1) = \boldsymbol{u_1}^T\sum_{x}^{}\boldsymbol{u_1}+\lambda _1(1-\boldsymbol{u_1}^T\boldsymbol{u_1})$

其中 λ_1 便是拉格朗日乘子（这里需要先假设 Σ_x 的特征值各不相同），接下来我们需要使用 L 分别对 u_1 和 λ_1 求偏导。求导结果为：

$\frac{\partial L}{\partial \boldsymbol{u_1}}=2(\sum_{x}^{}\boldsymbol{u_1}-\lambda _1\boldsymbol{u_1})=0\\\\ \frac{\partial L}{\partial \lambda _1}=1-\boldsymbol{u_1}^T\boldsymbol{u_1}=0$

解得：

$\sum_{x}^{}\boldsymbol{u_1}=\lambda _1\boldsymbol{u_1}~~~and~~~\boldsymbol{u_1}^T\boldsymbol{u_1}=1$

这个式子说明，u_1 是 Σ_x，也即 x 协方差矩阵的特征向量，对应的特征值为 λ_1。而第一主成分方差：

$\mathrm{Var}(y_1)=\boldsymbol{u_1}^T\sum_{x}^{}\boldsymbol{u_1}=\lambda_1\boldsymbol{u_1}^T\boldsymbol{u_1}=\lambda_1$

还记得第一主成分方差是最大的，也因此第一主成分对应的特征值也是最大的。由此我们可以想象一下每个主成分 y_i 都对应着一个 u_i 和 λ_i。但还需证明，从求解 u_2 开始。

首先，要想求解 u_2，需要保证 y_1 和 y_2 不相关（这是主成分最一开始的要求），这里的不相关就是概率论中的不相关（Uncorrelatedness），即 y_1 和 y_2 的协方差为 0：

$cov(y_1,y_2)=E[y_1y_2^T]-E[y_1]E[y_2]=0$

我们在最一开始假设了 x 的均值为 0，所以 y_i 的期望也为 0。那么，我们只需让 y_1 和 y_2 的外积的均值为 0 就好。

$0=E[y_1y_2^T]=E[\boldsymbol{u_1}^T\boldsymbol{x}\boldsymbol{x}^T\boldsymbol{u_2}]=\boldsymbol{u_1}^T\sum_{x}\boldsymbol{u_2}=\lambda_1\boldsymbol{u_1}^T\boldsymbol{u_2}$

这个式子的头和尾，单独拿出来看：

$\boldsymbol{u_1}^T\boldsymbol{u_2}=0$

它告诉我们 u_1 和 u_2 是正交的。现在我们可以写出 u_2 的目标函数和约束了：

$\boldsymbol{u_{2}}^{*}=\arg \max \boldsymbol{u_2}^T\sum_{x}^{}\boldsymbol{u_2}~~~s.t.~~~\boldsymbol{u_{2}}^{T}\boldsymbol{u_{2}}=1~~~and~~~\boldsymbol{u_1}^T\boldsymbol{u_2}=0$

再次构造拉格朗日函数：

$L(\boldsymbol{u_2},\lambda_2,\gamma )=\boldsymbol{u_2}^T\sum_{x}^{}\boldsymbol{u_2}+\lambda_2(1-\boldsymbol{u_{2}}^{T}\boldsymbol{u_{2}})+\gamma\boldsymbol{u_1}^T\boldsymbol{u_2}$

对 u_2、λ_2、γ 求偏导，可以得到：

$\sum_x\boldsymbol{u_2}+\frac{\gamma}{2}\boldsymbol{u_1}=\lambda_2\boldsymbol{u_2},~~~\boldsymbol{u_{2}}^{T}\boldsymbol{u_{2}}=1~~~and~~~\boldsymbol{u_1}^T\boldsymbol{u_2}=0$

对于上式，可以两边同时左乘 u_1^T：

$\boldsymbol{u_1}^T\sum_x\boldsymbol{u_2}+\frac{\gamma}{2}\boldsymbol{u_1}^T\boldsymbol{u_1}=\lambda_1\boldsymbol{u_1}^T\boldsymbol{u_2}+\frac{\gamma}{2}= \lambda_2\boldsymbol{u_1}^T\boldsymbol{u_2}=0$

不难发现 γ=0。所以式子简化成了：

$\sum_x\boldsymbol{u_2}=\lambda_2\boldsymbol{u_2},~~~\boldsymbol{u_{2}}^{T}\boldsymbol{u_{2}}=1~~~and~~~\boldsymbol{u_1}^T\boldsymbol{u_2}=0$

和前面的分析一样，这个式子说明，u_2 是 Σ_x，也即 x 协方差矩阵的特征向量，对应的特征值为 λ_2。而第二主成分方差：

$\mathrm{Var}(y_2)=\boldsymbol{u_2}^T\sum_{x}^{}\boldsymbol{u_2}=\lambda_2\boldsymbol{u_2}^T\boldsymbol{u_2}=\lambda_2$

事实上，当 Σ_x 的特征值都是单值的时候，以上的分析方法对于任意 i ≠ j 的 u_i 和 u_j 都是成立的，各位感兴趣可以自己试一下。整个求解过程的核心是：① 找出一个 u_i，② u_j 和 u_i 不相关。

现在，我们需要考虑 Σ_x 有重复特征值的情况。首先需要明确，相同特征值对应的不同特征向量是互相正交的，其次，这些正交的特征向量组成的特征向量空间和其他特征值的特征向量是正交的。由此可以推断，也可以证明，若 λ_1 是重根，则 u_1 可以选择 λ_1 对应的任意一个特征向量，而 u_2 可以选择 λ_1 对应的特征向量空间中与 u_1 正交的任意其他特征向量，以此类推。所以，我们做一个小总结：

若想取得 x 的 d 个主要成分，则我们可以把 x 的协方差矩阵 Σ_x 的特征值从大到小排列，取前 d 个最大的特征值所对应的特征向量共同组成一个矩阵 U，则 y_i = u_i^T x，y = U^Tx，这便是 x 的 d 个主成分。

可能你还会有疑问，前面的分析基于的假设是 x 是均值为 0 的多元随机变量，那么均值不为 0 呢？我们在实际运用PCA的过程中，中心化（均值归零）的操作是很容易实现的，若想深究，可以参考 Generalized Principal Component Analysis 一书。

十一、求解PCA（二）：几何学原理（拉格朗日乘子法求解）

其实几何视角看PCA，正是我们前面说过的投影、降维、重建等一系列问题，这里就不再细说了，不太理解的读者可以向前翻翻。这里直接列出我们要求解的式子（也是从前面搬过来的）：

上式的约束条件是：

$U^TU=I_d$

构造拉格朗日函数：

$L=tr(UU^T\boldsymbol{X}\boldsymbol{X}^T)+tr((I_d-U^TU)\Lambda )$

其中 Λ 是实对称矩阵。上式中，添加了很多处理手法。首先目标函数是 tr()，也就是一个标量。而如果直接引入 I_d - U^TU 做为约束项的话，就会出现标量+矩阵的问题，这是不利于求解的。而将 λ 变为 Λ 的原因是 λ 并不只有一个，Λ = diag{λ_1, λ_2, ... , λ_d}，而 I_d - U^TU 也是一个对角阵，那么 Λ(I_d - U^TU) 和 (I_d - U^TU)Λ 是相等的。

于是，我们用到前文第八部分提到的求导公式，对 L 求偏导后可以得到：

$XX^TU=U\Lambda$

因此：

$U^TXX^TU=\Lambda$

在这一情况下，稍微改写一下目标函数便可以得到，其最大值为 trace(Λ) 的最大值：

$\arg \max _{U} tr(UU^T\boldsymbol{X}\boldsymbol{X}^T)\\=\arg \max _{U} tr(U^T\boldsymbol{X}\boldsymbol{X}^TU)\\ =\arg \max _{U} tr(\Lambda)$

也许会有人问为什么 Λ 是对角矩阵呢？如果不是对角矩阵会如何？这其实很好处理，我们在上式中，用 U' 替换 U，用 Λ 代表一个非对角矩阵（但仍然是实对称，对上式两边取转置就可以解释）。而我们注意到 U 是正交的，它仅仅起到了旋转的功能，也因此，我们可以令 U' = UR。则上式变为：

$U'^TXX^TU'=R^TU^TXX^TUR=R^T\Lambda R=\Lambda'$

显然，因为 Λ 是实对称矩阵，其必然可以相似对角化。若 R 选择为 Λ 的特征向量的时候，Λ' 就是对角矩阵了。所以，不失一般性，我们可以直接令 Λ 为对角矩阵。也即：

$U^TXX^TU=\Lambda$

现在，因为 Λ 是 dxd 的对角阵。要想得到 tr(Λ) 的最大值，则需要将 U 取为 XX^T 的前 d 个最大的特征值对应的特征向量。换言之，也就是 X = U_X Σ_X V_X^T 的前 d 个奇异向量（用到了奇异值分解，Singular Value Decomposition），或者说 U_X 的前 d 列（左奇异向量）。这里稍微展开一下，SVD也即奇异值分解，是这样写的：

$X=U_X\sum_XV_X^T$

其中 U_X 和 V_X 分别称为 X 的左奇异矩阵和右奇异矩阵，U_X 和 V_X 都是正交矩阵，而 Σ_X 则是 X 的奇异值矩阵，其是一个对角矩阵，Σ_X = diag{σ_1, σ_2, ... , σ_D}，σ_i 就是矩阵 X 的奇异值，并且规定 σ_1>=σ_2>=....>=σ_D>0。

接下来可以计算两个矩阵：

$XX^T=U_X\sum_XV_X^TV_X\sum_XU_X^T=U_X\sum_X^2U_X^T\\\\ X^TX=V_X\sum_XU_X^TU_X\sum_XV_X^T=V_X\sum_X^2V_X^T$

第二个式子暂时用不到，第一个式子再稍微变化一下：

$\Lambda_X=U_X^TXX^TU_X=U_X^TU_X\sum_X^2U_X^TU_X=\sum_X^2\\\\ \Lambda~~=U^TXX^TU~~=U^TU_X\sum_X^2U_X^TU~~=\sum^2$

我们就可以发现，σ_i^2 = λi 。

由此，目标函数的最大值（迹形式）为：

$\sum_{i=1}^{d}\sigma_i^2$

其中 σi 为矩阵 X 的第 i 个奇异值。同时我们注意到，在前面的化简过程中有这样一步：

这个式子中，Σ符号右边的地方是一个标量（1xD@DxD@Dx1，@表示矩阵乘法）。这意味着，重建的 X_i 到原始的 X_i 之前的最小二乘（Least Squares）误差为：

$LSE = \sum_{i}^{}\boldsymbol{X_i}^T(I_D-UU^T)\boldsymbol{X_i}\\ =\sum_{i}^{}tr(\boldsymbol{X_i}^T(I_D-UU^T)\boldsymbol{X_i})\\ =\sum_{i}^{}tr(\boldsymbol{X_i}^T\boldsymbol{X_i})-tr(\boldsymbol{X_i}^TUU^T\boldsymbol{X_i})\\ =tr(\boldsymbol{X}^T\boldsymbol{X})-tr(\boldsymbol{X}^TUU^T\boldsymbol{X})\\=tr(\boldsymbol{X}\boldsymbol{X}^T)-tr(\boldsymbol{U}^T\boldsymbol{X}\boldsymbol{X}^T\boldsymbol{U})\\ =tr(\sum_X^2)-tr(\sum^2)$

其中第四行到第五行用到了迹的交换性质。而上式，PCA的最小二乘误差最终等于：

$\sum_{i=d+1}^D\sigma_i^2$

通过上述分析，我们了解到 SVD 给出了 PCA 的最佳优化方案。显然，上述分析过程也是基于 0 均值进行的，若不是零均值该如何处理？请参阅其他书籍。

十二、求解PCA（三）：梯度下降法（或上升法）求解

这里我提一下思路，首先梯度下降法的原理就不说了。这里的目标函数是前面构造的拉格朗日函数，必须要将约束项算进去，否则很难收敛（实际上带上约束项也很难收敛）。而不管是目标函数取了最大值还是最小值，在最值点处的梯度一定为 0。这便和拉格朗日求解过程中令偏导为 0 有着异曲同工之妙。原理分析只能到这，因为梯度下降法是个计算方法，没法理论论证最终结果。各位可以写代码试一试。如果觉得矩阵求导比较困难，可以参考 The Matrix Cookbook 一书，里面的矩阵求导公式非常详尽。

@20231018 终于用 MATLAB 实现了人脸 PCA 的梯度下降法。

前半部分降维方法参考 Matthew Turk, Alex Pentland; Eigenfaces for Recognition. J Cogn Neurosci 1991; 3 (1): 71–86. doi: Eigenfaces for Recognition | Journal of Cognitive Neuroscience | MIT Press

%导入图片
clc;
clear;
pathOfImages = "E:\Projects\Eigenface\dataset";
imds = imageDatastore(pathOfImages, ...
    "ReadSize", 1)

%读取图片
numOfFaces = 16;

greyFaces = zeros(224, 224, numOfFaces);
flattenGreyFaces = zeros(numOfFaces, 224 * 224);

for i = 1: numOfFaces
    face = readimage(imds, i);
    face = imresize(face, [224, 224]);

    greyFaces(:, :, i) = face(:, :, 1);
    flattenGreyFaces(i, :) = reshape(face(:, :, 1), 1, []);
end

size(greyFaces), size(flattenGreyFaces)
greyFaces;
flattenGreyFaces;

montage(uint8(greyFaces));

%计算平均脸和差异脸
averageGreyFace = mean(greyFaces, 3);
averageFlattenGreyFace = mean(flattenGreyFaces, 1);

size(averageGreyFace), size(averageFlattenGreyFace)

biasFlattenGreyFace = flattenGreyFaces - averageFlattenGreyFace;

imshow(uint8(averageGreyFace));

%GD
% rng(5) % 0.02
rng(15) % 0.009
V = randn(numOfFaces);
lr = 0.00001;
numOfEpoochs = 2000;

calMSE(V, biasFlattenGreyFace)

optimizedV = gradientDescent(V, biasFlattenGreyFace, lr, numOfEpoochs);

%评价
calMSE(optimizedV, biasFlattenGreyFace)

[eigV, ~] = eig(biasFlattenGreyFace * biasFlattenGreyFace');
sum ((optimizedV - eigV) .^ 2, "all") / (16 * 16)

%构建特征脸
%% 单通道
eigenFlattenGreyFaces = optimizedV' * biasFlattenGreyFace;
size(eigenFlattenGreyFaces)
eigenGreyFaces = permute(reshape(eigenFlattenGreyFaces, numOfFaces, 224, 224), [2, 3, 1]);
montage(uint8(eigenGreyFaces))

%重建
reconstructedFlattenGreyFaces = optimizedV * eigenFlattenGreyFaces + averageFlattenGreyFace;
reconstructedGreyFaces = permute(reshape(reconstructedFlattenGreyFaces, numOfFaces, 224, 224), [2, 3, 1]);

montage(uint8(reconstructedGreyFaces))

%函数定义
function loss = calMSE(V, A)
    loss = sum((V * V' * A - A) .^ 2, "all") / (size(A, 1) * size(A, 2));
end

function grad = calGradient(V, A)
    grad = (V * V' * A - A) * A' * V / (size(A, 1) * size(A, 2)) * 4;
end

function optimizedParam = gradientDescent(V, A, lr, numOfEpochs)
    optimizedParam = V;
    for i = 1: numOfEpochs
        grad = calGradient(optimizedParam, A);
        optimizedParam = optimizedParam - lr * grad;
        loss = calMSE(optimizedParam, A);
        fprintf('Epoch %d, Loss: %f\n', i, loss);
    end
end

十三、后记

写了很多，必然有疏漏的地方，望各位批评指正。

写这篇文章的初衷是我自己在学习的过程中查阅资料很费劲，缺少一个领路人一样的存在吧，所以我希望这篇文章可以给各位指明一些方向。

十四、推荐阅读

很多东西我已经把链接放在文中了，另外一些资料列在下方（给各位道歉，我用到的资料都是英文版，并不清楚国内有没有译注）：

程序员的数学3 线性代数，人民邮电出版社, 图灵程序设计丛书·程序员的数学, 2016，（日）平冈和幸，堀玄著（注：很诙谐，但是我带入不了）
The Matrix Cookbook，November 15, 2012, 2012，Kaare Brandt Petersen, Michael Syskind Pedersen（注：工具书）
Introduction to Linear Algebra，Wellesley-Cambridge Press, 5th, 2016，Gilbert Strang（注：有很多宝贵的线性代数思想，有MIT公开课，B站可查）
Linear Algebra and Learning from Data，Wellesley-Cambridge Press, 1, 2019，Gilbert Strang（注：线代和深度学习的综合，有MIT公开课，B站可查）
Generalized Principal Component Analysis，Springer, Interdisciplinary applied mathematics 40, 1st ed., 2016，René Vidal, Yi Ma, S.S. Sastry (auth.)（注：所有PCA相关的都在这）

此致，献给所有热爱知识的人。