PCA分析公式推导

2021-12-29

字数统计: 7.3k | 阅读时长≈ 34 分钟

吴老师的CS229课程里讲到了PCA分析，但是感觉有些地方不清晰，查阅了一些资料，自己整理了一下。

特征值分解与奇异值分解

首先我们需要了解一些特征值分解与奇异值分解的一些背景知识。

特征值分解

特征值和特征向量

矩阵 $A$ 的特征值和特征向量满足关系： \[ A x=\lambda x, \quad x \neq 0 \] 我们注意到对于任何特征向量 $x \in \mathbb{C}^{n}$ ，和特征值 $t \in \mathbb{C}$ ，$A(c x)=c A x=c \lambda x=\lambda(c x)$ ，因此 $c x$ 一样是 $A$ 的特征向量。因此，当我们提到与特征值 $\lambda$ 的特征向量 $x$ 时，我们一般假定特征向量 $x$ 是标准化的，即 $\|x\| = 1$ 。但是，我们注意到，$x$ 和 $-x$ 都是特征向量，因此我们无法保证特征向量的符号。

我们可以将上式改写为： \[ |(\lambda I-A)|=0 \] 方阵 $A \in \mathbb{R}^{n \times n}$ 一定有 $n$ 个特征值（可能是复数）。

特征值和特征向量有以下性质：

方阵 $A$ 的迹等于所有特征值之和。

\[ \operatorname{tr} A=\sum_{i=1}^{n} \lambda_{i} \]

证明如下： \[ \begin{aligned} \operatorname{tr}(A) &=\operatorname{tr}\left(Q \Lambda Q^{-1}\right) \quad \because \text { 对协方差矩阵进行特征值分解 } \\ &=\operatorname{tr}\left(Q^{-1} Q \Lambda\right) \\ &=\operatorname{tr}(\Lambda) \\ &=\sum_{i=1}^{n} \lambda_{i} \end{aligned} \]

方阵 $A$ 的行列式等于其所有特征值的乘积（缺证明）

\[ |A|=\prod_{i=1}^{n} \lambda_{i} \]

方阵 $A$ 的秩等于其非零特征值的数目（缺证明）
如果 $A$ 非奇异，那么 $1/\lambda$ 和 $x$ 为 $A^{-1}$ 的特征值和特征向量。

证明，对 $Ax = \lambda x$ 左右两项均左乘一个 $A^{-1}$ ，得到 \[ A^{-1}A x = A^{-1} \lambda x \\ (1/\lambda) x = A^{-1}x \]
对角矩阵$D=\operatorname{diag}\left(d_{1}, \ldots d_{n}\right)$ 的特征值就是其对角线元素

证明：因为，$Dx = [d_{1}x_{1} \quad d_{2}x_{2} \quad \cdots \quad d_{n}x_{n} ]^{T}$ ， $d_{i}x = [d_{i}x_{1} \quad d_{i}x_{2} \quad \cdots \quad d_{i}x_{n} ]^{T}$ ，那么我们只需要设 $x_{i}$ 为任一非零实数，而其他 $x_{j}$ ($j \neq i$ ) 均为0，那么就满足 $Dx = d_{i}x$ ，即对角线元素 $d_{i}$ 为 $D$ 的特征值。

我们可以将特征方程写为 \[ A X=X \Lambda \] 其中 $X \in \mathbb{R}^{n \times n}$ ，每一列表示一个特征向量，而 $\Lambda$ 为对角矩阵，对角线为特征值，即： \[ X \in \mathbb{R}^{n \times n}=\left[\begin{array}{cccc} \mid & \mid & & \mid \\ x_{1} & x_{2} & \cdots & x_{n} \\ \mid & \mid & & \mid \end{array}\right], \Lambda=\operatorname{diag}\left(\lambda_{1}, \ldots, \lambda_{n}\right) \] 证明如下： \[ \begin{aligned}A X &=A\left[x_{1}, x_{2}, \ldots, x_{n}\right] \\&=\left[\lambda_{1} x_{1}, \lambda_{2} x_{2}, \ldots, \lambda_{n} x_{n}\right] \\&=\left[x_{1}, x_{2}, \ldots, x_{n}\right] \Lambda \\&=X \Lambda\end{aligned} \] 如果 $A$ 的特征向量之间线性无关，那么 $X$ 可逆，因此存在下式 \[ A=X \Lambda X^{-1} \]

最后，我们的约束条件是 $A$ 有 $n$个线性无关的特征向量，而不是矩阵 $A$ 本身是否满秩，矩阵 $A$ 无论是否满秩都有可能有 $n$ 个线性无关的特征向量。

对称矩阵的特征值和特征向量

对称矩阵 $A$ 的特征值和特征向量有两个重要的性质，首先，对称矩阵所有的特征值均是实数；第二，对称矩阵的特征向量彼此正交（缺证明）。因此，此时的 $X$ 矩阵为正交矩阵，此时我们定义特征向量的矩阵为 $U$ ，即我们表示对称矩阵 $A$ 为下式 \[ A=Q \Lambda Q^{-1}=Q \Lambda Q^{T} \]

证明： $Q^{-1}=Q^{T}$ (其实就是正交矩阵的基本性质) \[ \begin{aligned} &Q^{T} Q=\left[\begin{array}{c} - q_{1} - \\ - q_{2} - \\ \cdots \\ - q_{n} - \end{array}\right]\left[\begin{array}{ccc} | & | & \vdots & | \\ q_{1} & q_{2} & \vdots & q_{n} \\ | & | & \vdots & | \end{array}\right]\\ &=\left[\begin{array}{cccc} q_{1}^{T} q_{1} & q_{1}^{T} q_{2} & \cdots & q_{1}^{T} q_{n} \\ q_{2}^{T} q_{1} & q_{2}^{T} q_{2} & \cdots & q_{2}^{T} q_{n} \\ \cdots & \cdots & \cdots & \cdots \\ q_{n}^{T} q_{1} & q_{n}^{T} q_{2} & \cdots & q_{n}^{T} q_{n} \end{array}\right]\\ \\ &=I \quad \because 对角线为模长度，非对角线根据特征向量两两正交的性质均为0\ \end{aligned} \]

此时，我们可以表示二次型为 \[ x^{T} A x=x^{T} U \Lambda U^{T} x=y^{T} \Lambda y=\sum_{i=1}^{n} \lambda_{i} y_{i}^{2} \] 其中 $y=U^{T} x$ 。

因为 $y_{i}^{2} $，因此，二次型的符号完全由特征值 $\lambda_{i}$ 决定。如果所有的特征值均满足 $\lambda_{i} > 0$ ，那么这个矩阵便是正定矩阵；如果所有的特征值 $\lambda_{i} \geq 0$ （也就是存在特征值为0），那么这个矩阵就是半正定矩阵。相似地，如果所有的特征值满足 $\lambda_{i}<0$ 或 $\lambda_{i} \leq 0$ ，此时矩阵 $A$ 就是一个负定或半负定矩阵。如果矩阵 $A$ 同时有正数和负数的特征值，说明 $A$ 为不定矩阵

奇异值分解

下面这段话全部来自于别人的文章[3]

当给定一个大小为 $m \times n$ 的矩阵 $A$ ，虽然矩阵 $A$ 不一定是方阵，但大小为 $m \times m$ 的 $A A^{T}$ 和 $n \times n$ 的 $A^{T} A$ 却是对称矩阵，可以进行特征值分解。若 $A A^{T}=P \Lambda_{1} P^{T}$ ， $A^{T} A=Q \Lambda_{2} Q^{T}$ 则矩阵 $A$ 的奇异值分解为 \[ A=P \Sigma Q^{T} \] 满足 \[ P^{T} = P^{-1}\\ Q^{T} = Q^{-1}\\ \] 其中，矩阵 $P=\left(\vec{p}_{1}, \vec{p}_{2}, \ldots, \vec{p}_{m}\right)$ 的大小为 $m \times m$ ，列向量 $\vec{p}_{1}, \vec{p}_{2}, \ldots, \vec{p}_{m}$ 是 $A A^{T}$ 的特征向量，也被称为矩阵 $A$ 的左奇异向量 (left singular vector) ; 矩阵 $Q=\left(\vec{q}_{1}, \vec{q}_{2}, \ldots, \vec{q}_{n}\right)$ 的大小为 $n \times n$ ，列向量 $\vec{q}_{1}, \vec{q}_{2}, \ldots, \vec{q}_{n}$ 是 $A^{T} A$ 的特征向量，也被称为矩阵 $A$ 的右奇异向量 (right singular vector) ; 矩阵 $\Lambda_{1}$ 大小为 $m \times m$ ，矩阵 $\Lambda_{2}$ 大小为 $n \times n$ ，两个矩阵对角线上的非零元素相同 (即矩阵 $A A^{T}$ 和矩阵 $A^{T} A$ 的非零特征值相同，推导过程见下) ; 矩阵 $\Sigma$ 的大小为 $m \times n$ ，位于对角线上的元素被称为奇异值 (singular value)。

接下来，我们来看看矩阵 $\Sigma$ 与矩阵 $A A^{T}$ 和矩阵 $A^{T} A$ 的关系。令常数 $k$ 是矩阵 $A$ 的秩，则 $k \leq \min (m, n)$ ，当 $m \neq n$ 时，很明显，矩阵 $\Lambda_{1}$ 和矩阵 $\Lambda_{2}$ 的大小不同，但矩阵 $\Lambda_{1}$ 和矩阵 $\Lambda_{2}$ 对角线上的非零元素却是相同的，若将矩阵 $\Lambda_{1}$ (或矩阵 $\Lambda_{2}$ ) 对角线上的非零元素分别为 $\lambda_{1}, \lambda_{2}, \ldots, \lambda_{k}$ （因为 $\operatorname{rank}\left(\boldsymbol{A} \boldsymbol{A}^{\mathrm{T}}\right)=\operatorname{rank}\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}\right)=\operatorname{rank}(\boldsymbol{A})$ ），其中，这些特征值也都是非负的（需证明 $AA^T$ 和 $A^TA$ 是半正定的，见下），再令矩阵 $\ Sigma$ 对角线上的非零元素分别为 $\sigma_{1}, \sigma_{2}, \ldots, \sigma_{k}$ ，则 \[ \sigma_{1}=\sqrt{\lambda_{1}}, \sigma_{2}=\sqrt{\lambda_{2}}, \ldots, \sigma_{k}=\sqrt{\lambda_{k}} \] 即非零奇异值的平方对应着矩阵 $\Lambda_{1}$ (或矩阵 $\Lambda_{2}$ ) 的非零特征值，到这里，我们就不难看出奇异值分解与对称对角化分解（特征分解）的关系了，即我们可以由对称对角化分解得到我们想要的奇异值分解。为了便于理解，在这里，给定一个大小为 $2 \times 2$ 的矩阵 $A=\left[\begin{array}{cc}4 & 4 \\ -3 & 3\end{array}\right]$ ，虽然这个矩阵是方阵，但却不是对称矩阵，我们来看看它的奇异值分解是怎样的。由 $A A^{T}=\left[\begin{array}{cc}32 & 0 \\ 0 & 18\end{array}\right]$ 进行对称对角化分解，得到特征值为 $\lambda_{1}=32 ， \lambda_{2}=18$ ，相应地，特征向量为 $\vec{p}_{1}=(1,0)^{T} ， \vec{p}_{2}=(0,1)^{T}$; 由 $A^{T} A=\left[\begin{array}{cc}25 & 7 \\ 7 & 25\end{array}\right]$ 进行对称对角化分解，得到特征值为 $\lambda_{1}=32 ， \lambda_{2}=18 ，$ 相应地，特征向量为 $\vec{q}_{1}=\left(\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}\right)^{T} ， \vec{q}_{2}=\left(-\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}\right)^{T}$ 。取 $\Sigma=\left[\begin{array}{cc}4 \sqrt{2} & 0 \\ 0 & 3 \sqrt{2}\end{array}\right]$ ，则矩阵 $A$ 的奇异值分解为 \[ \begin{aligned} &A=P \Sigma Q^{T}=\left(\vec{p}_{1}, \vec{p}_{2}\right) \Sigma\left(\vec{q}_{1}, \vec{q}_{2}\right)^{T} \\ &=\left[\begin{array}{ll} 1 & 0 \\ 0 & 1 \end{array}\right]\left[\begin{array}{cc} 4 \sqrt{2} & 0 \\ 0 & 3 \sqrt{2} \end{array}\right]\left[\begin{array}{cc} \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\ -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \end{array}\right]=\left[\begin{array}{cc} 4 & 4 \\ -3 & 3 \end{array}\right] \end{aligned} \]

如果 $A$ 是一个对称方阵，则 $AA^T=A^TA$ ，此时左奇异向量和右奇异向量构成的矩阵也是相等的，即 $P=Q$ 。更为神奇的是，对称方阵的奇异值分解与特征分解结果相同，证明见下。

证明：矩阵 $A A^{T}$ 和矩阵 $A^{T} A$ 的非零特征值相同

下面这段话来自网络[4]

For any $m \times n$ and $n \times m$ matrices $A$ and $B$, the nonzero eigenvalues of $A B$ and $B A$ are the same. Namely, if $A B v=\lambda v$ with $\lambda \neq 0$ and $v \neq 0$, then $B v \neq 0$ （因为 $A B v=\lambda v \neq 0$ ，如果 $B v = 0$，则 $ABv = 0$，与原条件相悖，因此 $B v \neq 0$ ）and $B A(B v)=B(A B v)=\lambda B v$ ，即 $\lambda$ 同样为 $BA$ 的特征值，因此 $AB$ 和 $BA$ 的非零特征值相同。

证明：矩阵 $A A^{T}$ 和矩阵 $A^{T} A$ 半正定

首先，证明 $A^{T} A$ 半正定，证明如下： \[ x^{T}\left(A^{T} A\right) x=\left(x^{T} A^{T}\right)(A x)=(A x)^{T}(A x) \geq 0 \] 同理，可证 $A A^{T}$ 半正定。

证明：对称方阵的奇异值分解与特征分解结果相同

当 $A$ 为对称方阵时，$AA^T=A^TA = A^2$ 。

假设 $\lambda$ 和 $u$ 为 $A$ 的一个特征值和相应的特征方程，此时存在： \[ \begin{aligned} &A^{2} u \\ &=A(Au) \\ &=A(\lambda u) \\ &=\lambda A u \\ &=\lambda^2 u \\ \end{aligned} \] 因此， $\lambda^2$ 和 $u$ 为 $A^2$ 的特征值和特征向量。

因此，此时，$P, Q$ 为特征分解的特征向量矩阵，而 $\Sigma$ 为特征值组成的对角矩阵。也就是说，对称方阵的奇异值分解与特征分解结果相同

PCA 预处理

在使用PCA分析前，我们需要先进行特征缩放。

Let $\mu=\frac{1}{m} \sum_{i=1}^{m} x^{(i)}$.

Replace each $x^{(i)}$ with $x^{(i)}-\mu$.

Let $\sigma_{j}^{2}=\frac{1}{m} \sum_{i}\left(x_{j}^{(i)}\right)^{2}$

Replace each $x_{j}^{(i)}$ with $x_{j}^{(i)} / \sigma_{j}$.

前两步将均值改为 0，后两步将方差均调整为1，使得不同特征具有相同的取值范围 ( “scale”) 。如果你的特征本身就已经是范围一致，那就不用进行后两步，举个例子，假如你的输入特征是图像像素的灰度，那么你的特征取值都是 $\{0,1, \ldots, 255\}$ ，那么就不需要对范围做处理。

PCA 直观解释

吴老师说 PCA 分析结果可以有八九种解释的方式，但是我知道的只有两种。

第一种解释方式是找到一个超平面，使得样本点在这个超平面上尽可能分开，也就是投影点的方差尽可能大。如下图，叉形点是二维空间的样本点，我们想找到某个单位向量 $u$ ，使得数据投影到 $u$ 方向的点的方差最大化。

第二种解释是我们希望所有的样本到这个超平面的距离足够近（所有样本点到超平面距离的平方和最小），对应着下图的虚线。

推导第一主成分的方向

这里就是说，如果我们就像上面图中一样，找一个单位向量 $u$ 构成的一条直线方向，使得数据投影到 $u$ 方向的点的方差最大化。这里设样本数为 $m$ ，特征数为 $n$，我们知道给定一个单位向量 $u$ 和某个点 $x$ ，$x$ 投射到 $u$ 上投影点的长度为 $x^{T}u$ （内积的性质）。我们上面提到了最大化投影点的方差，因此，我们需要最大化下式[5]： \[ \begin{aligned} &\max\frac{1}{m} \sum_{i=1}^{m}\left(x^{(i)^{T}} u\right)^{2} \\ &=\frac{1}{m} \sum_{i=1}^{m} (u^{T} x^{(i)}) (x^{(i)^{T}} u) \\ &=\frac{1}{m} \sum_{i=1}^{m} u^{T} x^{(i)} x^{(i)^{T}} u \\ &=u^{T}\left(\frac{1}{m} \sum_{i=1}^{m} x^{(i)} x^{(i)^{T}}\right) u \\ \\ & s.t. \quad \|\mu\|_{2} = 1 \end{aligned} \] 这里，一个隐含的性质是，投影的点的均值也为 $\boldsymbol{0}$ ，证明如下： \[ \begin{aligned} &\frac{1}{m} \sum_{i=1}^{m}\left(x^{(i)^{T}} u\right) \\ &=\frac{1}{m} \sum_{i=1}^{m}\left(x^{(i)^{T}} \right) u \\ &= 0u \quad \because x进行了中心化，均值为 0\\ &=0 \end{aligned} \] 我们发现这个解为 $x$ 的协方差矩阵 $\Sigma=\frac{1}{m} \sum_{i=1}^{m} x^{(i)}{x}^{(i)^{T}}$ 的 主特征值(principal eigenvector) （特征值中模最大的特征值称为主特征值，相应的特征向量称为主特征向量）。

证明：我们可以通过构建拉格朗日乘子式来求解上面的问题： \[ \begin{aligned} L(u, \lambda) &=u^{T} \Sigma u-\lambda\left(u^{T} u-1\right) \\ &=u^{T} \Sigma u-\lambda\left(u^{T} I u-1\right) \\ \nabla_{u} L &=\left(\Sigma+\Sigma^{T}\right) u-\lambda\left(I+I^{T}\right) u \quad \because \frac{\partial\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{X}\right)}{\partial \boldsymbol{X}}=\left(\boldsymbol{A}+\boldsymbol{A}^{\mathrm{T}}\right) \boldsymbol{X} \\ &=2(\Sigma u-\lambda u) \\ \end{aligned} \] 因此，下式成立，即 $u $ 和 $\lambda$ 分别为 $\Sigma$ 的特征向量和特征值。 \[ \Sigma u=\lambda u \] 代入原式，得到我们要找的特征向量 $u$ 就是主特征向量。 \[ \begin{aligned} &\operatorname{max}\left(u^{T} \Sigma u\right) & \\ =& \max \left(u^{T} \lambda u\right) \\ =& \max (\lambda) \end{aligned} \] 于是，我们得到投影数据最佳的一维子空间，也就是第一主成分。

下面，我们将其一般化，推导 PCA 前 $n'$个主成分。

PCA一般推导

推导过程主要来自于西瓜书和南瓜书[6][7]。

最小投影距离

假设 $\mathrm{m}$ 个 $\mathrm{n}$ 维数据 $\{x^{(1)}, x^{(2)}, \ldots, x^{(m)}\}$ 都已经进行了中心化，即 $\sum_{i=1}^{m} x^{(i)}=0$ 。经过投影变换后得到的新坐标系为 $\left\{w_{1}, w_{2}, \ldots, w_{n}\right\}$, 其中 $w_{i}$ 是标准正交基向量 ($n \times 1$ 维向量)，即 $\|w_{i}\|_{2}=1, w_{i}^{T} w_{j}=0$ 。

如果我们将数据从 $n$ 维降到 $n^{\prime}$ 维，即丟弃新坐标系中的部分坐标，则新的坐标系为 $\left\{w_{1}, w_{2}, \ldots, w_{n^{\prime}}\right\}$。样本点 $x^{(i)}$ 在 $n^{\prime}$ 维坐标系中的投影为: $z^{(i)}=\left(z_{1}^{(i)}, z_{2}^{(i)}, \ldots, z_{n^{\prime}}^{(i)}\right)^{T}$. 其中， $z_{j}^{(i)}=w_{j}^{T} x^{(i)}$ 是 $x^{(i)}$ 在低维坐标系里第 $j$ 维的坐标。

如果我们用 $z^{(i)}$ 来恢复原始数据 $x^{(i)}$, 则得到的恢复数据 $\hat{x}^{(i)}=\sum_{j=1}^{n^{\prime}} z_{j}^{(i)} w_{j}=W z^{(i)}$, 其中， $W$ 为标准正交基组成的 $n \times n'$ 矩阵。 \[ W = \left[\begin{array}{ccc} | & | & \vdots & | \\ w_{1} & w_{2} & \vdots & w_{n'} \\ | & | & \vdots & | \end{array}\right]\\ \] 带入上式，易知，$z^{(i)}=W^{T} x^{(i)}$ 。并且存在，$W^TW=I$ ，证明如下：

\[ \begin{aligned} &W^{T} W=\left[\begin{array}{c} - w_{1} - \\ - w_{2} - \\ \cdots \\ - w_{n'} - \end{array}\right]\left[\begin{array}{ccc} | & | & \vdots & | \\ w_{1} & w_{2} & \vdots & w_{n'} \\ | & | & \vdots & | \end{array}\right]\\ &=\left[\begin{array}{cccc} w_{1} \cdot w_{1} & w_{1} \cdot w_{2} & \cdots & w_{1} \cdot w_{n'} \\ w_{2} \cdot w_{1} & w_{2} \cdot w_{2} & \cdots & w_{2} \cdot w_{n'} \\ \cdots & \cdots & \cdots & \cdots \\ w_{n'} \cdot w_{1} & w_{n'}\cdot w_{2} & \cdots & w_{n'} \cdot w_{n'} \end{array}\right]\\ \\ &=I \quad \because 对角线为模长度，非对角线根据特征向量两两正交的性质均为0\ \end{aligned} \] 现在我们考虑整个样本集，我们莃望所有的样本到这个超平面的距离足够近，即最小化下式: \[ \sum_{i=1}^{m}\left\|\hat{x}^{(i)}-x^{(i)}\right\|_{2}^{2} \] 将这个式子进行整理，可以得到: \[ \begin{aligned} \sum_{i=1}^{m}\left\|\hat{x}^{(i)}-x^{(i)}\right\|_{2}^{2} &=\sum_{i=1}^{m}\left\|W z^{(i)}-x^{(i)}\right\|_{2}^{2} \\ &=\sum_{i=1}^{m}\left(\mathbf{W} \boldsymbol{z}^{(i)}-\boldsymbol{x}^{(i)}\right)^{\mathrm{T}}\left(\mathbf{W} \boldsymbol{z}^{(i)}-\boldsymbol{x}^{(i)}\right) \\ &=\sum_{i=1}^{m}\left(\boldsymbol{z}^{(i)T} \mathbf{W}^{\mathrm{T}} \mathbf{W} \boldsymbol{z}^{(i)}-\boldsymbol{z}^{(i)T} \mathbf{W}^{\mathrm{T}} \boldsymbol{x}^{(i)}-\boldsymbol{x}^{(i)T} \mathbf{W} \boldsymbol{z}^{(i)}+\boldsymbol{x}^{(i)T} \boldsymbol{x}^{(i)}\right) \\ &=\sum_{i=1}^{m}\left(\boldsymbol{z}^{(i)T} \boldsymbol{z}^{(i)}-2 \boldsymbol{z}^{(i)T} \mathbf{W}^{\mathrm{T}} \boldsymbol{x}^{(i)}+\boldsymbol{x}^{(i)T} \boldsymbol{x}^{(i)}\right) \quad \because \boldsymbol{z}^{(i)T} \mathbf{W}^{\mathrm{T}} \boldsymbol{x}^{(i)} \text{ 和 } \boldsymbol{x}^{(i)T} \mathbf{W} \boldsymbol{z}^{(i)} 互为转置且为标量，因此二者相等\\ &=\sum_{i=1}^{m} z^{(i) T} z^{(i)}-2 \sum_{i=1}^{m} z^{(i) T} z^{(i)}+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \\ &=-\sum_{i=1}^{m} z^{(i) T} z^{(i)}+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \\ &=-\sum_{i=1}^{m} \operatorname{tr} \left(z^{(i)} z^{(i) T}\right)+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \quad \because \operatorname{tr} \left(z^{(i)} z^{(i) T}\right) = z^{(i) T} z^{(i)}\\ &=- \operatorname{tr} \left(\sum_{i=1}^{m} z^{(i)} z^{(i) T}\right)+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \\ &=- \operatorname{tr} \left(\sum_{i=1}^{m} W^{T} x^{(i)} x^{(i)T} W \right)+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \\ &=-\operatorname{tr}\left(W^{T}\left(\sum_{i=1}^{m} x^{(i)} x^{(i) T}\right) W\right)+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \\ &=-\operatorname{tr}\left(W^{T} X^{T} X W\right)+\sum_{i=1}^{m} x^{(i) T} x^{(i)} \end{aligned} \] 注意到 $\sum_{i=1}^{m} x^{(i)} x^{(i) T}$ 是数据集的协方差矩阵（从协方差矩阵准取定义看，这个公式少了 $1/m$ ，但是对我们这里推导无影响，下同），可以写成设计矩阵的格式（行是样本，列是特征，因此 $X$ 是一个 $m \times n$ 的矩阵）。而 $\sum_{i=1}^{m} x^{(i) T} x^{(i)}$ 是一个常量。

最小化上式等价于: \[ \underbrace{\arg \min }_{W}-\operatorname{tr}\left(W^{T} X^{T} X W\right) \\ \text { s.t. } W^{T} W=I \]

构建拉格朗日函数如下[8][9] ，这里拉格朗日乘子矩阵 $\Theta$ 是 $n' \times n'$ 的矩阵，与约束条件维度相同 \[ \begin{aligned} L(\mathbf{W}, \Theta) &=-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\left\langle\Theta, \mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right\rangle \quad \because 这里是两个矩阵的内积，等于两个矩阵相同位置元素乘积之和，可用迹函数表示\\ &=-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\operatorname{tr}\left(\Theta^{\mathrm{T}}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right)\right) \end{aligned} \]

若此时仅考虑约束 $w_{i}^{\mathrm{T}} w_{i}=1$ ，即只考虑 $W^TW$ 的对角线元素的约束，此时拉格朗日乘子矩阵 $\Theta$ 为对角矩阵，令新的拉格朗日乘子矩阵为 $\Lambda=\operatorname{diag}\left(\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n^{\prime}}\right)$ , 此时的拉格朗日函数为 \[ L(\mathbf{W}, \Lambda)=-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\operatorname{tr}\left(\Lambda^{\mathrm{T}}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right)\right) \] 对拉格朗日函数关于 $\mathbf{W}$ 求导可得 \[ \begin{aligned} \frac{\partial L(\mathbf{W}, \Lambda)}{\partial \mathbf{W}} &=\frac{\partial}{\partial \mathbf{W}}\left[-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\operatorname{tr}\left(\Lambda^{\mathrm{T}}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right)\right)\right] \\ &=-\frac{\partial}{\partial \mathbf{W}} \operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\frac{\partial}{\partial \mathbf{W}} \operatorname{tr}\left(\Lambda^{\mathrm{T}}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right)\right) \end{aligned} \] 由矩阵微分公式 $\frac{\partial}{\partial \mathbf{X}} \operatorname{tr}\left(\mathbf{X}^{\mathrm{T}} \mathbf{B} \mathbf{X}\right)=\mathbf{B} \mathbf{X}+\mathbf{B}^{\mathrm{T}} \mathbf{X}, \frac{\partial}{\partial \mathbf{X}} \operatorname{tr}\left(\mathbf{B X}^{\mathrm{T}} \mathbf{X}\right)=\mathbf{X B}^{\mathrm{T}}+\mathbf{X B}$ 可得 \[ \begin{aligned} \frac{\partial L(\mathbf{W}, \Lambda)}{\partial \mathbf{W}} &=-2 \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}+\mathbf{W} \Lambda+\mathbf{W} \Lambda^{\mathrm{T}} \\ &=-2 \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}+\mathbf{W}\left(\Lambda+\Lambda^{\mathrm{T}}\right) \\ &=-2 \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}+2 \mathbf{W} \Lambda \end{aligned} \] 令 $\frac{\partial L(\mathbf{W}, \Lambda)}{\partial \mathbf{W}}=0$ 可得 \[ \begin{aligned} -2 \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}+2 \mathbf{W} \Lambda &=\mathbf{0} \\ \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W} &=\mathbf{W} \Lambda \end{aligned} \] 将 $\mathrm{W}$ 和 $\Lambda$ 展开可得 \[ \mathbf{X}^{\mathrm{T}} \mathbf{X}\boldsymbol{w}_{i}=\lambda_{i} \boldsymbol{w}_{i}, \quad i=1,2, \ldots, n^{\prime} \]

显然, 此式为矩阵特征值和特征向量的定义式, 其中 $\lambda_{i}, w_{i}$ 分别表示矩阵 $\mathbf{X}^{\mathrm{T}} \mathbf{X}$ 的特征值和单位特征向量。由于以上是仅考虑约束 $w_{i}^{\mathrm{T}} w_{i}=1$ 所求得的结果, 而 $w_{i}$ 还需满足约束 $w_{i}^{\mathrm{T}} \boldsymbol{w}_{j}=0(i \neq j)$ 。观察 $\mathbf{X}^{\mathrm{T}} \mathbf{X}$ 的定义可知, $\mathbf{X}^{\mathrm{T}} \mathbf{X}$ 是一个实对称矩阵, 实对称矩阵的不同特征值所对应的特征向量之间相互正交, 同一特征值的不同特征向量可以通过施密特正交化使其变得正交, 所以通过上式求得的 $w_{i}$ 可以同时满足约束 $\boldsymbol{w}_{i}^{\mathrm{T}} \boldsymbol{w}_{i}=1, \boldsymbol{w}_{i}^{\mathrm{T}} \boldsymbol{w}_{j}=0(i \neq j)$ 。根据拉格朗日乘子法的原理可知, 此时求得的结果仅是最优解的必要条件, 而且 $\mathbf{X}^{\mathrm{T}} \mathbf{X}$ 有 $n$ 个相互正交的单位特征向量, 所以还需要从这 $n$ 个特征向量里找出 $n^{\prime}$ 个能使得目标函数达到最优值的特征向量作为最优解。将 $\mathbf{X}^{\mathrm{T}} \mathbf{X} W= W \Lambda$ 代人目标函数可得[9]

\[ \begin{aligned} \min _{\mathrm{W}}-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right) &=\max _{\mathbf{W}} \operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right) \\ &=\max _{\mathbf{W}} \operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W} \Lambda \right) \\ &=\max _{\mathbf{W}} \operatorname{tr}\left( \Lambda \right) \quad \because \mathbf{W}^{\mathrm{T}} \mathbf{W} = \mathbf{I} \\ &=\max _{\mathbf{W}} \sum_{i=1}^{n^{\prime}} \lambda_{i} \end{aligned} \]

显然, 此时只需要令 $\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n^{\prime}}$ 和 $w_{1}, w_{2}, \ldots, w_{n^{\prime}}$ 分别为矩阵 $\mathrm{XX}^{\mathrm{T}}$ 的前 $n^{\prime}$ 个最大的特征值和单位特征向量就能使得目标函数达到最优值。

最大投影方差

我们知道样本点在超平面上的投影为 $z^{(i)}=W^{T} x^{(i)}$ ，首先我们证明投影点的均值同样为 $\mathbf{0}$ 。 \[ \begin{aligned} & \frac{1}{m} \sum_{i=1}^{m} z^{(i)} \\ =& \frac{1}{m} \sum_{i=1}^{m} W^{T} x^{(i)} \\ =& W^{T}\left(\frac{1}{m} \sum_{i=1}^{m} x^{\left(i\right)}\right) \\ =& W^{T} \boldsymbol{0} \\ =& \boldsymbol{0} \end{aligned} \]

因此，投影点的协方差矩阵为：

\[ \begin{aligned} & \frac{1}{m} \sum_{i=1}^{m} z^{(i)} z^{(i) T} \\ =& \frac{1}{m} \sum_{i=1}^{m} W^{T} x^{(i)} x^{(i) T} W \\ =& \frac{1}{m} W^{T} \sum_{i=1}^{m}\left( x^{(i)} x^{(i) T} \right) W \\ =& \frac{1}{m} W^{T} \sum_{i=1}^{m}\left( x^{(i)} x^{(i) T} \right) W \\ =& \frac{1}{m} W^{T} \mathbf{X}^{\mathrm{T}} \mathbf{X} W \\ \end{aligned} \]

投影点的方差为其协方差矩阵的迹，即为 $\frac{1}{m} \operatorname{tr}\left(W^{T} X^{T} X W\right)$

于是，此时我们的优化目标为

\[ \begin{array}{ll} \underbrace{\arg \max }_{W} & \operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right) \\ \text { s.t. } & \mathbf{W}^{\mathrm{T}} \mathbf{W}=\mathbf{I} \end{array} \]

这和上面的最小投影距离的优化目标一致。

选择合适的 $n'$

当你对协方差矩阵 $\frac{1}{m} \mathbf{X}^{\mathrm{T}} \mathbf{X}$ 进行奇异值分解或特征值分解后，你可以得到总共 $n$ 个特征值。这时你可以指定一个阈值 $\mathrm{t}$, 这个阈值 $\mathrm{t}$ 在 $(0,1)$ 之间，表示保留的方差比例，一般惯例可选 0.90，0.95 和 0.99。假如我们的 $\mathrm{n}$ 个特征值为 $\lambda_{1} \geq \lambda_{2} \geq \ldots \geq \lambda_{n}$ ，则n’可以通过下式得到:

\[ \frac{\sum_{i=1}^{n^{\prime}} \lambda_{i}}{\sum_{i=1}^{n} \lambda_{i}} \geq t \]

吴老师提到过这个公式，很多人的博客也提到了这个公式，但是没有告诉你这个公式怎么来的。这个公式在我看来没有这么一目了然，应该给证明过程的。

我思考了半天，才弄懂原理。这个公式蕴含的一个前提是，协方差矩阵所有特征值之和等于样本点的方差。

\[ \begin{aligned} \operatorname{Var}(x) &=\frac{1}{m} \sum_{i=1}^{m}\left\|x^{(i)}\right\|^{2} \\ &=\operatorname{tr}\left(\operatorname{Cov}(x)\right) \\ &=\operatorname{tr}\left(\frac{1}{m} X^{T} X\right) \\ &=\operatorname{tr}\left(Q \Lambda Q^{-1}\right) \quad \because 对协方差矩阵进行特征值分解\\ &=\operatorname{tr}\left(Q^{-1} Q \Lambda\right) \\ &=\operatorname{tr}(\Lambda) \\ &=\sum_{i=1}^{n} \lambda_{i} \end{aligned} \]

因此，当你需要计算方差比例的时候，总方差可以通过对所有特征值求和得到，也可以通过对协方差矩阵对角线所有元素求和得到。

疑问

为什么对称矩阵就一定可以进行特征值分解（存在 $n$ 个线性无关的特征向量），并且其特征向量是一个正交矩阵？

我现在知道这个性质来自于谱定理(Spectral theorem)，吴老师说这可能是线性代数中最重要的一个定理，以后有机会再看吧。
为什么用 $z^{(i)}$ 来恢复原始数据 $x^{(i)}$ 的公式为 $\hat{x}^{(i)}=\sum_{j=1}^{n^{\prime}} z_{j}^{(i)} w_{j}=W z^{(i)}$ ？

这种转换其实就是基变换，基变换的内容可以看我的博客矩阵乘法的内涵和本质。$\mathbb{R}^{n}$ 的标准基为 $\{e_1, e_2, \cdots, e_n\}$ ，$n$ 个特征向量 $\left\{w_{1}, w_{2}, \ldots, w_{n}\right\}$ 也是$\mathbb{R}^{n}$ 的一组基（为方便称呼，我下面称呼其为”新基“），从”新基“转换到标准基的转移矩阵就是 $\mathbf{W}$ ，因此如果我们已知某个样本点在”新基“中的坐标为 $z^{(i)}$ ，其中我们只用前 $n'$ 维的信息，等同于我们设 $n'$ 维以后的维度的坐标均为 0，因此存在 \[ \hat{x}^{(i)}=W z^{(i)}=\sum_{j=1}^{n^{\prime}} z_{j}^{(i)} w_{j} + \sum_{j=n'+1}^{n} 0 \times w_{j} =\sum_{j=1}^{n^{\prime}} z_{j}^{(i)} w_{j} \]
不是很清楚拉格朗日函数的矩阵形式

其实压根就没有什么所谓的矩阵形式的拉格朗日函数，我们采用矩阵的目标是为了方便书写和运算，比如一个线性方程组你直接写出多个等式，你也可以用矩阵表达为 $Ax=b$ ，二者的内容是一致的。

这里第一次构建拉格朗日乘子式，限制条件是 $W^{T} W-I=\mathbf{0}$ ，如果我们把它写成多个等式的形式，其实就是 $n \times n = n^{2}$ 个等式的约束条件，即可以将 $W^{T} W-I=\mathbf{0}$ 写为下面这种形式： \[ \sum_{k=1}^{n} w_{ik}w_{jk} -1 = 0, \quad \text{ if } i=j \\ \sum_{k=1}^{n} w_{ik}w_{jk} = 0, \quad \text{ if } i \neq j \\ i = 1,\cdots, n; \quad j = 1,\cdots, n \] 此时构建的拉格朗日乘子式即为下式，后面的一项就是 $n^{2}$ 个等式的约束条件，只是表示为迹的形式。 \[ \begin{aligned} L(\mathbf{W}, \Theta) &=-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\left\langle\Theta, \mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right\rangle \quad \because 这里是两个矩阵的元素乘积，等于两个矩阵相同位置元素乘积之和，可用迹函数表示\\ &=-\operatorname{tr}\left(\mathbf{W}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \mathbf{W}\right)+\operatorname{tr}\left(\Theta^{\mathrm{T}}\left(\mathbf{W}^{\mathrm{T}} \mathbf{W}-\mathbf{I}\right)\right) \end{aligned} \] 但是我认为，这里其实没有 $n^{2}$ 个等式的约束条件，因为 $W^{T} W-I$ 是一个对称矩阵，实际上只有 $\frac{1}{2}n(n+1)$ 个等式约束条件（上三角的元素数目），因此这里应该要需要限制 $\Theta$ 为一个对称矩阵，从而将约束条件数目降为 $\frac{1}{2}n(n+1)$ 。

第二次构建拉格朗日函数时，我们只考虑 $W^TW$ 的对角线元素的约束，此时只有 $n$ 个约束条件，我们将 $\Theta$ 改为一个对角矩阵 $\Lambda$ 。

参考文献

吴恩达老师CS229课程讲义
https://www.cnblogs.com/pinard/p/6251584.html
https://zhuanlan.zhihu.com/p/26306568
https://math.stackexchange.com/questions/1249497/largest-eigenvalues-of-aa-equals-to-aa
https://www.bilibili.com/video/BV1EW411R7g6?p=14
《机器学习》- 周志华
https://www.cnblogs.com/pinard/p/6239403.html
https://math.stackexchange.com/questions/1104376/how-to-set-up-lagrangian-optimization-with-matrix-constrains
《南瓜书》

版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！