PCA 入门 从方差方向到重构误差

DATA NOTEBOOK阅读地图

6个主要章节

1段代码示例

6组图片与图表

先确认变量、单位与假设，再核对手算、图表和解释范围。

章节速览6 节 · 展开选择

1. 01先中心化，再讨论方向
2. 02最大方差为什么对应特征向量
3. 03四个点的完整手算
4. 04舍弃的方差，怎样变成重构误差
5. 05尺度、常量列与训练集边界
6. 06完整 NumPy 实现与合成核验

NumPy 统计能分别描述每一列的均值和方差。当两列一起变化时，还需要回答一个几何问题：把点投影到哪个方向，才能保留最多的变化？

主成分分析（PCA）把这个问题写成正交投影。本文用四个合成点走通中心化、协方差、方向、投影和重构，再说明“解释了 80% 方差”究竟对应多少误差。这里讨论实数、有限、完整的矩阵，不做缺失值填补、白化或非线性降维。

先中心化，再讨论方向

约定数据矩阵 X 有 n 行样本、p 列特征，n≥2。每一列减去训练集的均值，得到中心化矩阵 C；本文采用分母为 n−1 的样本协方差 S：

μ = 每列在训练集上的均值
C = X − μ
S = CᵀC / (n−1)

中心化把坐标原点移到数据均值。否则 XᵀX 同时包含均值偏移与围绕均值的变化；它所找到的方向可能主要指向远离原点的均值。PCA 重构时要把 μ 加回去，所以最终拟合的是经过均值的仿射子空间。

scikit-learn 的 PCA 文档采用中心化、默认不逐列缩放的约定，解释方差也使用 n−1 分母。若全篇一致改用 n，所有协方差特征值会乘同一个因子，方向和解释方差比例保持不变；涉及绝对误差的公式则必须跟随分母改变。

最大方差为什么对应特征向量

取单位列向量 v，中心化数据在这个方向上的得分为 z=Cv。因为各列已经中心化，z 的均值也为零，样本方差是：

Var(z) = zᵀz / (n−1) = vᵀSv，约束 vᵀv=1

S 是实对称半正定矩阵，可以取一组正交单位特征向量 v₁,…,vₚ，特征值按 λ₁≥…≥λₚ≥0 排列。把任意单位方向展开为 v=Σ_j a_j v_j，有 Σ_j a_j²=1，于是：

vᵀSv = Σ_j λ_j a_j² ≤ λ₁

v₁ 达到上界，就是第一主成分方向。在与 v₁ 正交的方向中重复这个论证，得到第二方向，依此类推。scikit-learn 的分解指南将 PCA 表述为依次解释最大方差的正交成分；上面的展开直接给出有限维情况下的理由。

方向与得分要分开。 v_j 是 p 维特征空间里的方向；Cv_j 才是 n 个样本在该方向上的坐标。把 v_j 和全部对应得分同时取反，投影与重构不变。若多个特征值相等，该特征空间内还可以旋转基底；若截断边界恰好切过一组相等特征值，最优子空间本身也可能不唯一。

不同主成分得分在这份中心化样本上的协方差为零，因为 v_iᵀSv_j=0；这不等于一般意义下的统计独立，也不能从主成分方向直接推出因果关系。

不同主成分得分在这份中心化样本上的协方差为零，因为 v_iᵀSv_j=0；这不等于一般意义下的统计独立，也不能从主成分方向直接推出因果关系。

方差比较沿方向进行。 原创几何图，使用下文四个合成点；箭头只标示正交方向，不把箭头长度当作方差刻度，也不代表真实数据或实测效果。查看原图

四个点的完整手算

原始四行数据取 (13,21)、(11,23)、(9,17)、(7,19)，均值为 μ=(10,20)。中心化后为 (3,1)、(1,3)、(−1,−3)、(−3,−1)，从而：

CᵀC = [[20,12],[12,20]]
S = [[20/3,4],[4,20/3]]
v₁ = (1,1)/√2，λ₁ = 32/3
v₂ = (1,−1)/√2，λ₂ = 8/3

保留 k=1 时，第一方向的解释方差比例为 (32/3)/(40/3)=0.8。每行的重构是 μ+(c·v₁)v₁；下面把每行两个坐标的平方误差相加。

样本	原始 x	中心化 c	第一得分 c·v₁	第二得分 c·v₂	一维重构 x̂	行平方误差
1	(13,21)	(3,1)	2√2	√2	(12,22)	2
2	(11,23)	(1,3)	2√2	−√2	(12,22)	2
3	(9,17)	(−1,−3)	−2√2	√2	(8,18)	2
4	(7,19)	(−3,−1)	−2√2	−√2	(8,18)	2

结果有三个不同口径：总平方误差 SSE=8；平均每行平方误差 SSE/n=2；平均每坐标平方误差 MSE=SSE/(np)=1。报告误差时，应写出分母。80% 解释率也不表示“80% 样本被正确分类”。

舍弃的方差，怎样变成重构误差

把前 k 个方向按列放进 V_k，满足 V_kᵀV_k=I。中心化重构为 Ĉ=CV_kV_kᵀ，残差与保留子空间正交。每行沿完整正交基展开后，平方长度等于各方向坐标的平方和；因此：

SSE_k = ‖C − CV_kV_kᵀ‖_F²
      = Σ_{j>k} ‖Cv_j‖²
      = (n−1) Σ_{j>k} λ_j

这个等式针对参与拟合的中心化数据、普通欧氏距离与正交投影。用同一模型重构测试集时仍能直接计算残差，但训练谱的尾和不再等于测试集误差。

为什么前 k 个方向也是最小误差方案？设任意 p×k 正交列矩阵为 Q，它保留的平方长度是 (n−1)tr(QᵀSQ)。在 S 的特征基中，这个迹可写成 Σ_j λ_j b_j，其中 0≤b_j≤1、Σ_j b_j=k；最大值为前 k 个特征值之和。总平方长度固定，保留最多就使残差最小。这个结论只比较维数为 k 的线性正交投影，并不覆盖任意非线性编码器。

损失等于被舍弃方向的平方坐标之和。 原创合成手算图；保留 32 与舍弃 8 使用同一平方长度单位，图中的误差口径来自四点算术，不是模型评测结果。查看原图

尺度、常量列与训练集边界

改变某一列的单位会改变欧氏距离和协方差，因此可能改变 PCA 方向。中心化保留原始单位下的变动；逐列标准化则把“相对自身标准差的变化”作为比较尺度。两者对应不同问题，应根据特征含义选择。标准化后的重构误差也在标准化坐标中计量，不能直接当成原始单位误差。

情况	应怎样理解
某列是常量	中心化后为零；若另做标准化，要单独处理零标准差
全部行完全相同	总方差为零，解释比例的分母为零；本文程序约定返回全零比例
p 大于 n	中心化矩阵秩至多 n−1，协方差里必有零特征值
特征值相等或非常接近	具体方向可能改变；比较投影矩阵、子空间和误差更合适
低方差方向含有预测信号	PCA 没有使用标签，不能保证保留它，也不能保证提升预测效果
存在异常点、缺失值或极端量级	先处理数据质量与数值范围；本文不提供稳健 PCA 或缺失值算法

均值、可选的标准差和 PCA 方向都只能在训练部分拟合。测试集减去保存的训练均值，再使用保存的方向，不重新估计测试集均值；交叉验证时，每一折都要在自己的训练部分重新拟合。scikit-learn 的数据泄漏说明明确把 PCA 列为需要遵守这一规则的预处理步骤。

自检：给每列都加上常数，PCA 方向会变吗？

如果重新在平移后的训练集拟合，中心化矩阵保持不变，协方差、方向子空间和训练重构误差也保持不变；均值与最终重构一起平移。对已有模型的新输入单独平移，却不调整模型均值，则不是同一情形。

自检：保留一维解释了 80%，测试集误差一定也是总量的 20% 吗？

不一定。80% 来自当前训练样本的谱分解。新样本可能沿舍弃方向变化更多；要使用训练均值与方向，直接计算测试集重构残差。若评价下游预测任务，还需要独立评价其目标指标。

完整 NumPy 实现与合成核验

下面需要安装 NumPy，本文实际运行环境为 Python 3 与 NumPy 1.24.4。程序直接对中心化矩阵做精简 SVD：若 C=U diag(s) Vᵀ，则 λ_j=s_j²/(n−1)，Vᵀ 的行就是方向。NumPy SVD 文档给出分解的形状、奇异值降序与正交性质；没有必要先显式形成协方差矩阵再用于正文实现。

此实现允许 0≤k≤min(n,p)，k=0 代表只用训练均值重构。精简 SVD 未列出的 p−min(n,p) 个协方差特征值均为零。模型的 components 只保存前 k 个方向，variance 与 ratio 则保留精简 SVD 的完整谱，便于观察舍弃部分；这与 scikit-learn 同名概念属性只保留所选成分的接口不同。程序统一拒绝布尔（包括混入数值列表的布尔元素）、复数、文本、非有限值、非法形状与产生非有限中间结果的输入；它仍是浮点教学实现，不提供任意精度或极端病态矩阵保证。

from dataclasses import dataclass
import numpy as np


def matrix(values, name, min_rows=1, width=None):
    try:
        elements = np.asarray(values, dtype=object)
        if any(isinstance(v, (bool, np.bool_)) for v in elements.flat):
            raise ValueError("boolean element")
        raw = np.asarray(values)
        if raw.dtype.kind not in "iuf":
            raise ValueError("expected real numeric values")
        x = np.array(raw, dtype=np.float64, copy=True)
    except (TypeError, ValueError, OverflowError) as exc:
        raise ValueError(name + ": invalid numeric matrix") from exc
    if x.ndim != 2 or x.shape[0] < min_rows:
        raise ValueError(name + ": invalid shape")
    if (width is None and x.shape[1] == 0) or (
        width is not None and x.shape[1] != width
    ):
        raise ValueError(name + ": invalid feature count")
    if not np.isfinite(x).all():
        raise ValueError(name + ": non-finite value")
    return x


def readonly(values):
    x = np.array(values, dtype=np.float64, copy=True)
    x.setflags(write=False)
    return x


@dataclass(frozen=True)
class PCA:
    mean: np.ndarray
    components: np.ndarray
    variance: np.ndarray
    ratio: np.ndarray
    training_sse: float

    def transform(self, values):
        x = matrix(values, "X", width=self.mean.size)
        try:
            with np.errstate(over="raise", invalid="raise"):
                z = (x - self.mean) @ self.components.T
            if not np.isfinite(z).all():
                raise ValueError("non-finite projection")
            return z
        except FloatingPointError as exc:
            raise ValueError("projection overflow") from exc

    def inverse_transform(self, values):
        z = matrix(values, "Z", width=self.components.shape[0])
        try:
            with np.errstate(over="raise", invalid="raise"):
                x = z @ self.components + self.mean
            if not np.isfinite(x).all():
                raise ValueError("non-finite reconstruction")
            return x
        except FloatingPointError as exc:
            raise ValueError("reconstruction overflow") from exc


def fit_pca(values, k):
    x = matrix(values, "X", min_rows=2)
    if isinstance(k, (bool, np.bool_)) or not isinstance(k, (int, np.integer)):
        raise ValueError("k must be an integer")
    if not 0 <= k <= min(x.shape):
        raise ValueError("k out of range")
    try:
        with np.errstate(over="raise", invalid="raise", divide="raise"):
            mean = x.mean(axis=0)
            centered = x - mean
            _, singular, vt = np.linalg.svd(centered, full_matrices=False)
            squared = singular * singular
            variance = squared / (x.shape[0] - 1)
            total = squared.sum()
            ratio = squared / total if total > 0 else np.zeros_like(squared)
            sse = float(squared[k:].sum())
        if not all(np.isfinite(a).all() for a in (mean, vt, variance, ratio)):
            raise ValueError("non-finite decomposition")
        if not np.isfinite(sse):
            raise ValueError("non-finite loss")
    except (FloatingPointError, np.linalg.LinAlgError) as exc:
        raise ValueError("decomposition failed or overflowed") from exc
    # 固定每行最大绝对分量的符号；不消除重特征值子空间的旋转自由度。
    for row in vt:
        pivot = int(np.argmax(np.abs(row)))
        if row[pivot] < 0:
            row *= -1
    return PCA(readonly(mean), readonly(vt[:k]), readonly(variance),
               readonly(ratio), sse)


if __name__ == "__main__":
    x = np.array([[13, 21], [11, 23], [9, 17], [7, 19]], dtype=float)
    model = fit_pca(x, 1)
    restored = model.inverse_transform(model.transform(x))
    print("mean:", model.mean)
    print("variance:", model.variance)
    print("ratio:", model.ratio)
    print("reconstructed:", np.round(restored, 8))
    sse = float(np.sum((x - restored) ** 2))
    print("SSE / mean row loss / coordinate MSE:",
          round(sse, 8), round(sse / len(x), 8), round(sse / x.size, 8))

合成例输出的均值为 [10,20]，谱约为 [10.66666667,2.66666667]，比例为 [0.8,0.2]；一维重构四行为 [12,22]、[12,22]、[8,18]、[8,18]，最后三个误差为 8、2、1。

内部核验位于 tests/verify-pca-article.cjs 与 tests/fixtures/pca-checks.py：提取并运行正文代码，使用独立协方差特征分解、二维闭式谱和随机正交旋转检查投影、损失及边界。验证只涉及这些合成矩阵与数学关系，没有运行真实任务，也没有把训练重构误差当成预测收益。