2.1 奇异值分解与主成分

动机

数据表虽然有很多列，独立的变化模式却可能少得多。两项测量可以一起变化；位于高维空间中的观测，也可以大致集中在一个平面附近。降维要用较少坐标描述这些结构。要把这一要求写成数学问题，我们必须说清楚：保留什么信息，以及怎样度量重建误差。

主成分分析（PCA）选取相互正交的方向，尽量保留数据的变异。奇异值分解（SVD）提供这些坐标，也精确记录舍弃的变异。同一计算还可以理解为：用一个低维仿射子空间拟合观测。这种联系解释了为什么要先中心化，以及为什么方差公式使用奇异值的平方。

PCA 不像回归那样，指定某一列为因变量；各项特征共同决定几何。因此，测量单位会影响答案。把米改成厘米，如果不相应调整预处理，就会改变普通欧氏距离。选择哪些特征、使用什么尺度，都是数学模型的一部分。

定义

设 $X\in\mathbb R^{n\times d}$ 包含 $n\ge2$ 个观测。第 $i$ 行是 $x_i^T$，其中 $x_i\in\mathbb R^d$ 是特征的列向量。以下始终采用“行代表观测、列代表特征”的约定：方向属于特征空间 $\mathbb R^d$，全部观测的得分则组成 $\mathbb R^n$ 中的向量。

定义

中心化与样本协方差

平均特征向量与中心化矩阵定义为

$$\mu=\frac1nX^T\mathbf1,\qquad H=I_n-\frac1n\mathbf1\mathbf1^T,\qquad A=HX=X-\mathbf1\mu^T.$$

这里 $\mathbf1\in\mathbb R^n$ 的每一项都是一。样本协方差矩阵是

$$C=\frac1{n-1}A^TA\in\mathbb R^{d\times d}.$$

由 $\mathbf1^TH=0$，可知 $A$ 每列的总和都是零。另外 $H^T=H$、$H^2=H$，所以再次中心化不会改变结果。常数向量位于 $H$ 的零空间，故 $\operatorname{rank}(A)\le\min(n-1,d)$。这解释了为什么原表即使具有满行秩，中心化后也可能出现零奇异值。减去均值移除的是整体位置，而不是围绕该位置的变化方向。

给定单位方向 $v\in\mathbb R^d$，观测 $i$ 的得分是 $(x_i-\mu)^Tv$。得分向量为 $Av$，其均值是零，所以样本方差为

$$\operatorname{Var}_{\mathrm{sample}}(Av) =\frac{\|Av\|_2^2}{n-1}=v^TCv.$$

定义

奇异值、载荷方向与得分

把约化 SVD 写成

$$A=U\Sigma V^T,\qquad p=\min(n,d),$$

其中 $U\in\mathbb R^{n\times p}$、$V\in\mathbb R^{d\times p}$ 的列向量分别正交归一，而 $\Sigma=\operatorname{diag}(\sigma_1,\ldots,\sigma_p)$，满足 $\sigma_1\ge\cdots\ge\sigma_p\ge0$。正奇异值的个数是 $r=\operatorname{rank}(A)$。右奇异向量 $v_j$ 是 PCA 的载荷方向，$Av_j=\sigma_ju_j$ 是对应的得分向量。

完整 SVD 的 $U$、$V$ 都是方形正交矩阵，中间则是矩形对角矩阵；它按需要补入零方向。对任何 $k\le r$，两种分解都给出相同的前 $k$ 个载荷组成的矩阵 $V_k=[v_1\ \cdots\ v_k]$，以及降维坐标 $Z=AV_k=U_k\Sigma_k$。载荷是用来量取坐标的方向，得分是某个观测在该方向的坐标；两者的维数也不同。

定理／命题

定理

主方向使保留方差最大

设 $A$ 已中心化，$C=A^TA/(n-1)$。如有需要，把奇异值序列补零至第 $d$ 项。协方差特征值为

$$\lambda_j=\frac{\sigma_j^2}{n-1}.$$

在所有单位方向中，$v_1$ 使得分方差最大。更一般地，若 $W\in\mathbb R^{d\times k}$ 满足 $W^TW=I_k$，则 $\operatorname{tr}(W^TCW)$ 的最大值为 $\sum_{j=1}^k\lambda_j$，由前 $k$ 个右奇异向量达到。

总变异非零时，第 $j$ 个主成分的解释方差比例，以及前 $k$ 个主成分的累计比例，分别是

$$\rho_j=\frac{\sigma_j^2}{\sum_{\ell=1}^r\sigma_\ell^2},\qquad R_k=\frac{\sum_{j=1}^k\sigma_j^2}{\sum_{\ell=1}^r\sigma_\ell^2}.$$

比例中的 $n-1$ 被约去。若以 $n$ 作为总体协方差的分母，特征值会改变，主方向及上述比例则不变。若全部观测相同，就有 $A=0$；所有得分都是零，没有正的总变异可用来定义百分比。

定理

投影误差与截断 SVD

对正交归一的 $W$，中心化坐标中的重建为 $AWW^T$，而

$$\|A-AWW^T\|_F^2 =\|A\|_F^2-\|AW\|_F^2.$$

因此 PCA 也使总平方重建误差最小。其解为

$$A_k=AV_kV_k^T=U_k\Sigma_kV_k^T,\qquad \|A-A_k\|_F^2=\sum_{j=k+1}^r\sigma_j^2.$$

此外，在所有秩不超过 $k$ 的矩阵 $B$ 中，$A_k$ 使 $\|A-B\|_F$ 最小。返回原坐标时，重建为 $\widehat X=\mathbf1\mu^T+A_k$。

Frobenius 范数的平方等于矩阵各项的平方和，所以上式合计全部观测、全部特征的偏差，不是每个观测的方差，也不是每项数据的平均误差。解释秩的选择时，要同时说明这个尺度。保留八成变异是否足够，取决于用途；一个方差较小的特征仍可能很重要。最优化结论精确说明 PCA 所解的是哪一种平方距离问题。

降维坐标保留什么

得分表除了各列均值为零，列与列也相互正交。事实上，

$$\frac1{n-1}Z^TZ =\frac1{n-1}\Sigma_k^2 =\operatorname{diag}(\lambda_1,\ldots,\lambda_k).$$

因此不同得分的样本协方差为零，方差按大小排列，可从对角线直接读出每个保留坐标的贡献。零协方差描述的是这个得分表的二阶矩，仅凭它不能推断原来测量在统计上独立。PCA 构造正交线性坐标，使总变异最大，并没有为观测指定概率模型。

重建时沿相同载荷反向变换：把得分行乘以载荷矩阵的转置，就返回特征空间。重建行位于所选子空间，残差垂直于每个保留方向。这也提供计算检查：残差表乘以载荷矩阵，结果应为零（数值计算允许舍入误差）。投影会把不同原始观测映射为同一得分，所以仅靠保留得分不能恢复被舍弃的部分。

PCA 拟合直线时度量的是指定特征度量下最短的垂直距离；普通回归则把解释变量视为给定，沿因变量方向量取误差。因此，即使只有两项特征和一条直线，两种目标仍可能选择不同直线。使用算法前，要辨明误差模型允许哪些测量变化，以及最终需要的是预测，还是数据的几何表示。

证明思路

证明

从谱定理到主成分

$A^TA$ 对称且半正定，因为 $v^TA^TAv=\|Av\|^2\ge0$。按特征值从大到小选择一组正交归一特征向量。对正特征值 $\sigma_j^2$，令 $u_j=Av_j/\sigma_j$。于是 $u_i^Tu_j=v_i^TA^TAv_j/(\sigma_i\sigma_j)$，在 $i=j$ 时为一，其余为零。零特征值则给出 $\|Av_j\|^2=0$，所以 $Av_j=0$，不用除以零。这 $r$ 个 $u_j$ 张成列空间，补成 $\mathbb R^n$ 的正交归一基得到方形 $U$，并保留完整特征基 $V$。取 $\Sigma\in\mathbb R^{n\times d}$，前 $r$ 个对角项为正 $\sigma_j$，其余为零。逐列有 $AV=U\Sigma$，所以 $A=U\Sigma V^T$；矩形或秩不足矩阵也成立。

把单位向量展开为 $v=\sum_j a_jv_j$，则 $v^TCv=\sum_j\lambda_ja_j^2\le\lambda_1$，因为 $\sum_ja_j^2=1$。若选取多个方向，令 $q_j=\|W^Tv_j\|^2$。正交投影给出 $0\le q_j\le1$、$\sum_jq_j=k$，故 $\operatorname{tr}(W^TCW)=\sum_j\lambda_jq_j$ 不超过最大 $k$ 个特征值的总和。选取它们的特征向量就达到上界。

投影等式可逐个观测证明：每个中心化行向量分成子空间内的投影，以及与该空间垂直的余项；勾股定理把两部分的平方长度相加。再对所有行求和，就得到 Frobenius 等式。在 SVD 坐标中，被舍弃的各个秩一项也在 Frobenius 内积下正交，所以每项对误差的贡献恰是奇异值的平方。

为什么任意秩 $k$ 矩阵也不能给出更好的结果？其所有行都位于某个维数不超过 $k$ 的子空间 $S$。把 $A$ 每行正交投影到 $S$，误差不会大于选择该候选矩阵的对应行。而在所有这些子空间中，上述方差论证已经找出最佳者。这同时证明了低秩近似和几何 PCA 的结论。

截断对谱范数 $\|E\|_2=\max_{\|w\|_2=1}\|Ew\|_2$ 也最优。对 $0\le k<r$，

$$\min_{\operatorname{rank}(B)\le k}\|A-B\|_2 =\|A-A_k\|_2=\sigma_{k+1}.$$

把秩不超过 $k$ 的候选 $B$ 限制在 $\operatorname{span}\{v_1,\ldots,v_{k+1}\}$，因定义域维数为 $k+1$，其零空间含单位向量 $w$。所以 $Bw=0$，而 $\|(A-B)w\|_2=\|Aw\|_2\ge\sigma_{k+1}$。截断余项的最大奇异值恰为首个舍弃值，达到下界。谱误差度量单位方向的最坏放大，Frobenius 误差的平方则合计全部舍弃方向的平方损失。

例题

例题

1. 完整计算两项特征的 PCA

考虑四个观测：

$$X=\begin{bmatrix}12&5\\8&5\\10&6\\10&4\end{bmatrix},\qquad \mu=\begin{bmatrix}10\\5\end{bmatrix},\qquad A=\begin{bmatrix}2&0\\-2&0\\0&1\\0&-1\end{bmatrix}.$$

$A$ 的两列正交，故

$$A^TA=\begin{bmatrix}8&0\\0&2\end{bmatrix},\qquad C=\begin{bmatrix}8/3&0\\0&2/3\end{bmatrix}.$$

取 $v_1=(1,0)^T$、$v_2=(0,1)^T$，则 $\sigma_1=2\sqrt2$、$\sigma_2=\sqrt2$，而

$$U=\begin{bmatrix} 1/\sqrt2&0\\-1/\sqrt2&0\\0&1/\sqrt2\\0&-1/\sqrt2 \end{bmatrix},\qquad V=I_2.$$

相乘 $U\operatorname{diag}(2\sqrt2,\sqrt2)V^T$ 可重新得到 $A$。第一个得分向量是 $(2,-2,0,0)^T$，样本方差为 $8/3$，保留的变异比例是 $8/(8+2)=4/5$。只保留此坐标，再加回均值，得到

$$\widehat X= \begin{bmatrix}12&5\\8&5\\10&5\\10&5\end{bmatrix}.$$

最后两个观测各有大小为一的误差，总平方误差是 $2$，等于被舍弃的 $\sigma_2^2$。除以八项数据后，均方误差是 $1/4$，逐项的均方根误差是 $1/2$。

这个例子区分了四个常被混用的量：奇异值、奇异值平方、样本方差、解释比例。它们描述同一方向，却回答不同问题。奇异值是得分空间的长度；把它平方再除以 $n-1$ 才是样本方差。比例也要使用平方，因为正交分解中可相加的是平方能量，不是未平方的长度。

例题

2. 中心化改变拟合的几何

同一批原始观测满足

$$X^TX=A^TA+n\mu\mu^T =\begin{bmatrix}408&200\\200&102\end{bmatrix}.$$

非对角项变成 $200$，所以 $(1,0)^T$ 不是 $X^TX$ 的特征向量。直接分解 $X$，会选到另一个受整体位置影响的首方向。若只把原始观测投影到通过原点的横轴，整个第二坐标都被舍弃，平方误差为 $25+25+36+16=102$。

中心化 PCA 拟合的是通过 $(10,5)$ 的水平线，误差只有 $2$。两者所解的拟合问题不同；最佳原始秩一方向也不一定水平。上述比较说明，原始横轴投影不等于中心化后的重建。

对新观测 $(11,7)$，先减去训练数据的均值得到 $(1,2)$。第一得分为 $1$，用一个主成分重建后是 $(11,5)$。处理后续观测时，要沿用相同的均值和载荷，使所有得分使用同一坐标系。

仿射观点在此很清楚：保留的中心化子空间通过零，加回均值后才移回观测所在的位置。因此，维度可以用来描述数据周围的变化，不用花一个方向描述整体离开原点的位移。若要从得分表解码回原坐标，除了载荷，还需要保存平均向量。

例题

3. 特征单位与重复特征值

把第二项特征乘以三。中心化后的行是 $(2,0)$、$(-2,0)$、$(0,3)$、$(0,-3)$，散布矩阵为 $\operatorname{diag}(8,18)$。首载荷变成 $(0,1)^T$，保留比例为 $18/26=9/13$，秩一平方误差为 $8$。

另一个选择是把原来每列除以其样本标准差，即 $\sqrt{8/3}$ 和 $\sqrt{2/3}$。标准化后的散布矩阵是 $3I_2$，协方差矩阵是 $I_2$。每个单位方向的得分方差都为一，没有特别优先的单一首主成分；任意正交归一的一对方向都张成同一个二维空间。

标准化选择了另一种距离度量，让不同特征使用可比较的尺度。当单位较任意时，这可能有用；但如果变异大小本身有意义，标准化也会消除这项差异。解释保留比例前，先说清楚误差在原单位还是标准化单位中度量。零方差的列没有正标准差可除，它对中心化变异的贡献也是零。

幂迭代与伪逆

从实特征基到幂迭代

先备概念见 特征值、实对角化及可对角化判据。设实方阵 $B=VDV^{-1}$ 有实特征基，各列 $v_j$ 长度取一，但未必正交。一组模最大条件是

$$|\lambda_1|>\max_{j>1}|\lambda_j|,\qquad c=V^{-1}q_0,\qquad c_1\ne0.$$

从单位初值 $q_0$ 开始，每次相乘后归一化；在固定特征基中展开得到

$$q_{k+1}=\frac{Bq_k}{\|Bq_k\|_2},\qquad B^kq_0=c_1\lambda_1^k \left(v_1+\sum_{j>1}\frac{c_j}{c_1} \left(\frac{\lambda_j}{\lambda_1}\right)^k v_j\right).$$

其余系数相对首项几何衰减，方向就趋近首特征线；首系数非零也保证归一化分母不为零。这些比例属于特征基坐标，欧氏范数界还取决于基的几何。下文对称散布矩阵才可使用正交坐标的较强估计。

$c_1$ 由解 $Vc=q_0$ 得到，不能用 $v_1^Tq_0$ 代替：非正交基中，点积为零仍可有首系数，点积非零也未必有。缺少的首系数不会在相乘中出现；严格模间隔则排除另一个同量级模式持续竞争。

若 $\lambda_1<0$，归一化向量可正负交替，收敛的是直线。Rayleigh 商 $\rho(q)=q^TBq/(q^Tq)$ 在该直线两个单位代表附近连续，且不受符号影响，所以趋近 $\lambda_1$。非对称时没有对称矩阵的方差极值解释，中途商值也可超出实特征值范围。

例题

负主特征值与非正交基

取

$$B=\begin{bmatrix}-3&4\\0&1\end{bmatrix},\qquad V=\begin{bmatrix}1&1\\0&1\end{bmatrix},\qquad D=\operatorname{diag}(-3,1).$$

这里的基列尚未归一化。$q_0=(0,1)^T$ 的坐标为 $(-1,1)^T$，含首系数；把基归一化只会缩放系数，不改变零与非零。原始迭代为 $(4,1)^T$、$(-8,1)^T$，而 $q_2=(-8,1)^T/\sqrt{65}$ 有 $\rho(q_2)=-223/65<-3$。方向趋近横线，符号交替。若沿 $(1,1)^T$ 开始，则一直留在特征值一的方向，虽与横向特征向量点积非零。

对单位右特征向量，减去 $\lambda_1v_1v_1^T$ 仍把首向量映成零，但通常改变其他原特征向量。要保留其余特征对，可取左特征向量 $w_1^TB=\lambda_1w_1^T$、$w_1^Tv_1=1$，减去 $\lambda_1v_1w_1^T$。本例 $v_1=(1,0)^T$、$w_1=(1,-1)^T$ 给出

$$B+3v_1v_1^T=\begin{bmatrix}0&4\\0&1\end{bmatrix},\qquad B+3v_1w_1^T=\begin{bmatrix}0&1\\0&1\end{bmatrix}.$$

前者把原向量 $(1,1)^T$ 映成 $(4,1)^T$，配对更新则保留它。对称矩阵可取相同的左右单位特征向量，就恢复正交消去公式。

近似首 PCA 方向

只需要首主方向时，可用幂迭代代替完整分解。设 $S=A^TA$ 的特征值为 $\lambda_1>\lambda_2\ge\cdots\ge\lambda_d\ge0$，正交归一特征向量为 $v_j$。取单位初值 $q_0=\sum_jc_jv_j$，要求 $c_1\ne0$，反复计算

$$q_{k+1}=\frac{Sq_k}{\|Sq_k\|_2}.$$

迭代方向与 $S^kq_0$ 相同。提出 $\lambda_1^k$ 后，其余方向的系数带有 $(\lambda_j/\lambda_1)^k$，所以垂直部分与首分量之比不超过

$$\left(\frac{\lambda_2}{\lambda_1}\right)^k \frac{\sqrt{\sum_{j=2}^d c_j^2}}{|c_1|}.$$

方向因此以几何速率趋近 $v_1$，允许符号相反。归一化控制长度，不改变比例；谱间隔越小，收敛越慢。初值若垂直首方向，就永远没有该分量。单位向量的 $q_k^TSq_k$ 近似首特征值；已知精确首特征对后，$S-\lambda_1v_1v_1^T$ 可移除该模式。大型数据可直接计算 $Sq=A^T(Aq)$，不用显式建立散布矩阵。

例题

4. 跟踪主要散布方向

取 $S=\operatorname{diag}(4,1)$、$q_0=(1,1)^T/\sqrt2$。两次迭代给出

$$q_1=\frac{(4,1)^T}{\sqrt{17}},\qquad q_2=\frac{(16,1)^T}{\sqrt{257}}.$$

第二坐标与第一坐标之比依次为 $1$、$1/4$、$1/16$，展示谱比 $\lambda_2/\lambda_1$ 的作用；方向趋近首载荷 $(1,0)^T$。

选择最小范数最小二乘系数

对实矩阵 $M\in\mathbb R^{m\times p}$，用正奇异值写成 $M=U_r\Sigma_rV_r^T$，定义

$$M^+=V_r\Sigma_r^{-1}U_r^T,\qquad z_*=M^+y.$$

直接相乘 SVD，可核对四条 Penrose 等式：

$$MM^+M=M,\qquad M^+MM^+=M^+,\qquad (MM^+)^T=MM^+,\qquad(M^+M)^T=M^+M.$$

其中 $MM^+=U_rU_r^T$、$M^+M=V_rV_r^T$ 分别正交投影到列空间和行空间，均对称且幂等。满行秩时前者是 $I_m$，满列秩时后者是 $I_p$；方形可逆时 $M^+=M^{-1}$。秩不足时，两个投影辨明数据决定哪些响应与系数方向。

对任意 $y$，它是唯一最小范数最小二乘系数。因 $Mz_*=U_rU_r^Ty$ 是 $y$ 在列空间的正交投影，残差最小。其他最小二乘系数都是 $z_*+h$，其中 $h\in\ker M$。而 $z_*$ 属于 $V_r$ 张成空间，垂直零空间，所以 $\|z_*+h\|_2^2=\|z_*\|_2^2+\|h\|_2^2$；只有 $h=0$ 才同范数。

方程相容时，这给出最短精确解；不相容时，则在同样准确的拟合中选择最短系数。只对正奇异值取倒数，无须逆转奇异矩阵。同一 SVD 因而兼顾降维和系数选择。

在左奇异坐标中，$v_j$ 的系数是 $(u_j^Ty)/\sigma_j$。$y$ 垂直列空间的部分无法拟合，留在残差；零右奇异方向则允许加入零空间向量。最小范数选择把自由坐标设为零，每个正奇异值只需一次除法。若只保留部分正奇异方向，解的就是该张成空间内的受限拟合，不是原来不受限的最小范数最小二乘问题。

例题

5. 解决欠定拟合

设 $M=[1\ \ 2]$、$y=5$，则 $M^+=(1,2)^T/5$、$z_*=(1,2)^T$。所有精确解是 $z=(5-2t,t)^T$，且

$$\|z\|_2^2=(5-2t)^2+t^2=5(t-2)^2+5.$$

在 $t=2$ 达到最小值，核对了 $z_*$ 及平方范数 $5$。与它的差为 $(t-2)(-2,1)^T$，正是垂直 $(1,2)^T$ 的零空间方向。

计算得分与重建

以下程序实现第一个例子。数值 SVD 返回的是 $V$ 的转置，所以 Vt 的各行存放载荷方向。

import numpy as np

X = np.array([[12., 5.], [8., 5.], [10., 6.], [10., 4.]])
mean = X.mean(axis=0)
A = X - mean
U, s, Vt = np.linalg.svd(A, full_matrices=False)
k = 1
scores = A @ Vt[:k].T
reconstruction = scores @ Vt[:k] + mean
fractions = s**2 / np.sum(s**2)
squared_error = np.sum((X - reconstruction)**2)
print(fractions)      # [0.8 0.2]
print(squared_error)  # 2.0, up to rounding

数值程序可以同时反转一个载荷和对应得分列的符号，两者相乘的重建不变。比较实现时，可检查子空间、解释方差和重建矩阵。若截断位置两侧的奇异值相同，最优保留子空间也可能不唯一；若该处存在严格间隔，子空间就确定，虽然基向量的符号仍可改变。

常见错误

常见错误

把奇异值直接当成方差比例

第一例中，未平方的比例 $2\sqrt2/(2\sqrt2+\sqrt2)=2/3$ 不是解释方差比例，正确值是 $8/(8+2)=4/5$。同时要区分：$V_k$ 的列描述特征方向，$AV_k$ 的行描述观测坐标。混淆两者通常首先造成矩阵维数不符。

总结

先中心化各列、选定距离度量，再对中心化矩阵作 SVD。右奇异向量提供载荷，乘上这些方向就得到得分。协方差特征值等于奇异值平方除以 $n-1$。前几个主成分同时使保留方差最大、正交平方重建误差最小；返回原特征坐标时，要加回均值。

练习

思考检查

1. 四个观测的例子只保留一个主成分。求解释比例、总平方重建误差，以及新观测 (11,7) 的重建。

把中心化、投影、加回均值分开计算。

思考检查

2. 某中心化矩阵的正奇异值是 6、3、1。要至少保留 90% 方差，最少需要保留多少秩？余下平方重建误差是多少？

先把奇异值平方，再计算累计总和。

思考检查

3. 某协方差矩阵的特征值是 5、5、0。第一载荷是否唯一？前二维主子空间是否唯一？

区分保留空间内的重复值，以及跨过截断位置的重复值。

解答

解答 · 1. 投影并加回均值

解释比例是 $4/5$，平方误差是 $2$。新观测的中心化坐标为 $(1,2)$，第一得分为 $1$。中心化重建是 $(1,0)$，加回 $(10,5)$ 后得到 $(11,5)$。

解答 · 2. 从平方谱选择秩

平方总和为 $36+9+1=46$。一个主成分保留 $36/46$，低于 $0.9$；两个保留 $45/46$，高于 $0.9$，故最少保留秩二。被舍弃的平方奇异值是 $1$，也就是余下总平方重建误差。

解答 · 3. 基可改变，子空间确定

特征值五的特征空间中，每个单位向量都使方差最大，所以第一载荷除符号外仍不唯一。完整二维特征空间则唯一，因为其特征值与零有间隔。任意正交归一基都保留总方差 $10$；更换基只会旋转得分，不改变重建。

相关笔记

• 最小二乘与残差
• 二次函数的梯度下降
• 课程读物