# 2.10. 主成分分析
主成分分析（PCA）在特征空间里找出一组正交的方向，也就是主轴。第一条主轴捕获数据中最大的方差。第二条主轴与第一条正交，并在此前提下捕获剩余方差中最大的一份。后面每一条主轴都遵循同样的规律。
把数据投影到前 `k` 条主轴上，就得到保留方差最多的 `k` 维线性子空间；这个子空间同时也是让重建平方误差最小的那一个。保留方差最多的子空间，恰好也是丢弃信息最少的子空间。正是这条对偶性质，让 PCA 成为压缩、去噪、给高维数据画图时最常用的第一步。
RustyML 的 `PCA` 对标 `sklearn.decomposition.PCA`：在特征矩阵上调用 fit，用 `transform` 得到得分，用 `inverse_transform` 把得分映射回特征空间。
RustyML 的 `PCA` 和 scikit-learn 的 `PCA` 有 3 处不同。第一，它只做中心化，从不做缩放。你喂给它的数值，比在一条会自动标准化的流水线里更加要紧。
第二，求解器这个设置是一个名为 `SVDSolver` 的枚举，提供 3 种具体策略。该选哪一种，取决于特征数量和你想要的主成分数量。
第三，每一步分解都跑在手写的纯 Rust 代码上。既没有 LAPACK，也没有 BLAS。这让各求解器之间的取舍是实打实的，而不是纸上谈兵。
## 2.10.1. 你要打交道的接口
`PCA` 和 `SVDSolver` 都位于 `rustyml::machine_learning::decomposition` 之下。这两个类型也都能通过 [Prelude](../Chapter-01/1.5._Prelude与模块导入.md) 拿到。这个模型是无监督的：它接收一个 `f64` 特征矩阵，行是样本、列是特征，不使用任何标签。
```rust,ignore
// 构造。n_components 必须 > 0，这里会检查。
PCA::new(n_components: usize) -> Result<PCA, Error>
PCA::default() // n_components = 2，Full 求解器
pca.with_svd_solver(solver: SVDSolver) -> PCA // builder，消费 self 并返回
// 学习与投影。fit 接收 &mut self。
pca.fit(&x) -> Result<&mut PCA, Error>
pca.transform(&x) -> Result<Array2<f64>, Error> // (n_samples, n_components)
pca.fit_transform(&x) -> Result<Array2<f64>, Error>
pca.inverse_transform(&scores) -> Result<Array2<f64>, Error> // (n_samples, n_features)
// 拟合后的状态。mean、components、方差以及样本/特征数量在 fit 前为 None，
// fit 后为 Some(&...)。n_components 和 svd_solver 则随时可用。
pca.get_mean() -> Option<&Array1<f64>> // 每个特征的中心化均值
pca.get_components() -> Option<&Array2<f64>> // (n_components, n_features)，行即主轴
pca.get_explained_variance() -> Option<&Array1<f64>>
pca.get_explained_variance_ratio() -> Option<&Array1<f64>>
pca.get_singular_values() -> Option<&Array1<f64>>
pca.get_n_components() -> usize
pca.get_svd_solver() -> SVDSolver
pca.get_n_samples() -> Option<usize>
pca.get_n_features() -> Option<usize>
```
构造函数只检查 `n_components > 0`。更严的那条界限，`n_components <= min(n_samples, n_features)`，要等到 `fit` 时才检查，因为它取决于数据本身。用这些 getter 来判断模型是否已经拟合：在成功的 `fit` 调用把值填进去之前，每一个都返回 `None`。
`get_components` 把主轴作为 `(n_components, n_features)` 矩阵的**行**返回，这与 scikit-learn 的 `components_` 布局一致。所以 `transform` 就是中心化数据乘 `components.T`，`inverse_transform` 就是得分乘 `components`。
下面这个例子在一个小的二维数据集上跑完整流程，只保留 1 个主成分：
```rust
use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;
fn main() {
let x = array![
[2.5, 2.4],
[0.5, 0.7],
[2.2, 2.9],
[1.9, 2.2],
[3.1, 3.0],
[2.3, 2.7],
];
// fit_transform 需要 &mut self；它就是先 fit 再对同一份数据 transform。
let mut pca = PCA::new(1).unwrap();
let scores = pca.fit_transform(&x).unwrap();
println!("scores shape: {:?}", scores.shape()); // [6, 1]
println!("mean: {:?}", pca.get_mean().unwrap());
println!("component: {:?}", pca.get_components().unwrap());
println!("variance: {:?}", pca.get_explained_variance().unwrap());
println!("ratio: {:?}", pca.get_explained_variance_ratio().unwrap());
println!("singular: {:?}", pca.get_singular_values().unwrap());
}
```
`fit_transform(&x)` 不是一条会给出不同结果的融合捷径，它就是对同一个矩阵先调 `fit` 再调 `transform`，和手动分两步做完全一样。真正需要用两步写法的场景，是把*另一个*矩阵（比如新样本）投影到一个已经拟合好的模型上。由于 `fit` 返回 `&mut Self`，你也可以直接链式调用一个 getter，不过更常见的做法还是通过 `get_*` 方法把状态读回来。
## 2.10.2. 只做中心化，不做缩放
在分解数据之前，`fit` 会先减去每个特征的均值，并把它存进 `get_mean()`。这是 `fit` 唯一做的预处理。它不会除以标准差，不会做白化，也不会碰各列的量纲。
`transform` 复用训练时存下的均值来中心化新样本，`inverse_transform` 再把均值加回去。均值是拟合模型的一部分，不是每次调用时临时算出来的。
跳过这一步，是把 PCA 用歪的最常见原因。方差不具备尺度不变性：把一个特征从米换算成毫米，数值会乘以 1000，方差却要乘以 1000^2，也就是一百万倍。于是 PCA 会把第一主成分整个分配给那一列。
解法是在 fit 之前先标准化特征。标准化把每一列都拉到单位方差，让 PCA 站在同一起跑线上比较它们。用 [`standardize`](../Chapter-04/4.2._标准化与归一化.md) 并传入 `StandardizationAxis::Column`。下面的例子展示了前后的差别：
```rust
use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;
use rustyml::utils::standardize::{standardize, StandardizationAxis};
fn main() {
// 特征 0 的量纲极小，特征 1 的量纲极大。信息量相同，
// 方差却天差地别。
let x = array![
[0.1, 1000.0],
[0.2, 1005.0],
[0.3, 1002.0],
[0.4, 998.0],
[0.5, 1010.0],
];
let mut raw = PCA::new(2).unwrap();
raw.fit(&x).unwrap();
println!("raw ratios: {:?}", raw.get_explained_variance_ratio().unwrap());
let xs = standardize(&x, StandardizationAxis::Column).unwrap();
let mut scaled = PCA::new(2).unwrap();
scaled.fit(&xs).unwrap();
println!("standardized ratios: {:?}", scaled.get_explained_variance_ratio().unwrap());
}
```
在原始数据上，第一个比值接近 `1.0`：第 1 列的数值跨度淹没了第 0 列的信号，PC1 也几乎正好指向第 1 列。按列标准化之后，两个特征的贡献变得大致相当，比值也就反映出两者真实的相关结构。
要不要标准化是一个建模决策。如果特征本来就处在同一个有意义的量纲下，比如像素强度，那么只做中心化就够了。请务必有意识地做这个决定，因为 rustyml 不会替你做主。统计学上把这个选择描述为“对协方差矩阵做 PCA”还是“对相关矩阵做 PCA”。先标准化再做 PCA，等价于对相关矩阵做 PCA。
## 2.10.3. 选择求解器，以及背后的纯 Rust 机制
`SVDSolver` 决定 PCA 怎么计算主成分。这个名字不算准确：3 个变体里只有 1 个真的会构造 SVD。但接口在这 3 者之间保持一致：无论选哪一个，你拿到的主成分行、可解释方差、奇异值都一样，只是每个求解器走了不同的路线，成本和精度也各不相同。
```rust,ignore
pub enum SVDSolver {
Full, // 默认：对协方差矩阵做精确的特征分解
Randomized(u64), // 随机化范围查找，由 u64 播种
PowerIteration, // 用带收缩的幂迭代提取前 k 个特征对
}
```
| `Full`（默认） | 构建 `d x d` 协方差 `X^T X / (n - 1)`，并对其做精确分解（先 Householder 三对角化，再做隐式位移 QL）。 | 中小规模数据，样本和特征数大致都在 10,000 以下。精度接近机器精度，是稳妥的默认选择。 |
| `Randomized(seed)` | 把 `X` 草绘到一个 `k` 宽的随机子空间，跑几轮子空间迭代，再在降维后的空间里做一次小规模 SVD。全程不构建完整的 `d x d` 协方差。 | 又大又宽的数据（10,000+ 特征），看重速度、也能接受一点近似误差。`seed` 让随机草图可复现。 |
| `PowerIteration` | 先构建协方差，再用带 Hotelling 收缩的幂迭代只提取前 `k` 个特征对，而不做完整的特征分解。 | 只需要少数几个主成分（`k` 远小于 `d`），想跳过完整的 `O(d^3)` 特征求解。 |
每种求解器的成本遵循一个清晰的规律。`Full` 和 `PowerIteration` 都要构建稠密的 `d x d` 协方差，两者都要付出 `O(n * d^2)` 的构建开销和 `O(d^2)` 的存储开销。分岔发生在下一步：`Full` 要做一次完整的 `O(d^3)` 特征分解。
`PowerIteration` 则逐个提取 `k` 个特征对，每提取一个只用有限次矩阵-向量乘法，并在找下一个之前把当前这个收缩掉。当 `k` 远小于 `d` 时，这种方式在计算量上更划算。
`Randomized` 则完全不构建 `d x d` 矩阵：它的工作集是 `n x k` 的草图和 `k x d` 的降维投影，因此在 `d` 很大时是最省内存的选择。在小规模问题上，3 种求解器能吻合到好几位有效数字；只有当矩阵变大，近似求解器的速度优势或数值漂移才会显现出来。
```rust
use ndarray::array;
use rustyml::machine_learning::decomposition::{PCA, SVDSolver};
fn main() {
let x = array![
[2.5, 2.4],
[0.5, 0.7],
[2.2, 2.9],
[1.9, 2.2],
[3.1, 3.0],
[2.3, 2.7],
[2.0, 1.6],
[1.0, 1.1],
];
for solver in [SVDSolver::Full, SVDSolver::Randomized(42), SVDSolver::PowerIteration] {
let mut pca = PCA::new(1).unwrap().with_svd_solver(solver);
pca.fit(&x).unwrap();
let sigma = pca.get_singular_values().unwrap();
println!("{:?}: sigma_1 = {:.6}", solver, sigma[0]);
}
}
```
`PowerIteration` 和 `Randomized` 只暴露内部设置中很小的一部分。`PowerIteration` 没有公开的迭代次数或容差设置：它内部每个主成分最多跑 1000 次迭代，特征值容差为 `1e-6`，种子固定。这些值都无法通过公开 API 调节。
`PowerIteration` 还可能收敛失败。当数据的有效秩低于你请求的主成分数量时，就可能发生这种情况。举例来说，重复的列，或者特征之间存在精确的线性相关，都可能引发这个问题。一旦收缩步骤把可提取的方差耗尽，`fit` 就会返回 `Error::NotConverged`。`Full` 和 `Randomized` 都不会以这种方式失败，因为它们都不依赖逐个主成分的收敛检查。
`Randomized` 只暴露种子这一个参数，过采样量和子空间迭代次数都是写死的。所以 PCA 的全部配置面就是 `n_components` 加上 `with_svd_solver`，这是刻意做小的。
这 3 种求解器都跑在 crate 自己的 `machine_learning::linalg` 模块上。这个模块直接在 `ndarray` 数组上重新实现了对称特征分解、SVD（单边 Jacobi）和 QR（改进的 Gram-Schmidt）。因此没有 LAPACK 依赖需要安装、链接或打包。这与 crate 其余部分[纯 Rust 数值工具](../Chapter-06/6.0._数学工具.md)的做法一脉相承。底层的协方差与投影乘法用的是并行的 [`gemmkit`](../Chapter-06/6.2._矩阵乘法.md) 后端。
## 2.10.4. 解读方差，选定 `n_components`
fit 之后，有 3 个向量描述方差如何分布在保留的各条主轴上。`get_singular_values` 按降序返回中心化数据的奇异值。`get_explained_variance` 返回每个奇异值的平方，再除以 `(n - 1)`：这就是每个主成分方向上的真实方差，单位与原始数据一致。
`get_explained_variance_ratio` 把上面这些值各自除以中心化数据的总方差。这个总方差是完整的迹，是对*所有*特征求和，而不只是你保留的那些。所以比值里的每一项，都是该主成分解释了总体方差中多大的比例。
这个分母的取法很关键：如果保留的主成分比特征少，这些比值加起来会小于 1，差的那部分正是你丢弃的方差。
正是这条性质，让比值成为选 `n_components` 的合适工具。标准做法是碎石图加累积法则：先做一次满秩分解，读出累积比值，再保留能越过某个方差目标（比如 90%、95%、99%）所需的最少主轴数。
```rust
use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;
fn main() {
// 8 个样本，3 个特征。特征 2 几乎是特征 0 的线性回声，
// 所以数据实际上是低秩的。
let x = array![
[1.0, 0.2, 1.19],
[2.0, 0.1, 2.09],
[3.0, 0.4, 3.38],
[4.0, 0.2, 4.19],
[5.0, 0.5, 5.48],
[6.0, 0.3, 6.29],
[7.0, 0.1, 7.09],
[8.0, 0.6, 8.58],
];
// 先做满秩分解（n_components = min(n_samples, n_features) = 3），
// 再决定实际保留几条主轴。
let mut pca = PCA::new(3).unwrap();
pca.fit(&x).unwrap();
let ratio = pca.get_explained_variance_ratio().unwrap();
let mut cumulative = 0.0;
for (i, r) in ratio.iter().enumerate() {
cumulative += r;
println!("PC{}: individual = {:.4}, cumulative = {:.4}", i + 1, r, cumulative);
}
// 累积比值越过 95% 的最小 k。
let mut acc = 0.0;
let mut k = ratio.len();
for (i, r) in ratio.iter().enumerate() {
acc += r;
if acc >= 0.95 {
k = i + 1;
break;
}
}
println!("keep {} component(s) to retain >= 95% variance", k);
}
```
有个细节值得留意：你满秩 fit 一次是为了读出谱，如果想要更小的投影，再用选定的 `k` 重新 fit 一次。重新 fit 在这里很便宜，能给你一个 `transform` 恰好输出 `k` 列的模型。
如果在你的数据上重新 fit 代价很高，可以走一条捷径：满秩 fit 的前 `k` 个主成分行，和 `k` 成分 fit 得到的是同一批主轴，唯一的差别是下一节要讲的符号约定。所以你可以对满秩结果切片，而不必重新 fit。不要只拿单个孤立的比值去定 `k`，而要看累积曲线，它才能告诉你收益在哪里趋平。
## 2.10.5. 重建与 `inverse_transform`
`inverse_transform` 用公式 `reconstructed = scores * components + mean` 把得分映射回特征空间。当你保留全部 `min(n_samples, n_features)` 个主成分时，这次往返是无损的，至多有浮点误差。保留得更少时，重建就是数据在保留子空间上的正交投影。残差就是活在被丢弃主轴上的方差。度量这份残差，是给压缩代价标上一个数字的直接办法。
```rust
use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;
fn main() {
let x = array![
[1.0, 0.2, 1.19],
[2.0, 0.1, 2.09],
[3.0, 0.4, 3.38],
[4.0, 0.2, 4.19],
[5.0, 0.5, 5.48],
[6.0, 0.3, 6.29],
];
// 只保留一条主轴，再往返回到原始的三维特征空间。
let mut pca = PCA::new(1).unwrap();
pca.fit(&x).unwrap();
let scores = pca.transform(&x).unwrap(); // (6, 1)
let reconstructed = pca.inverse_transform(&scores).unwrap(); // (6, 3)
// 残差的 Frobenius 范数 = 随 PC2、PC3 一起被丢掉的方差。
let err: f64 = (&reconstructed - &x).iter().map(|d| d * d).sum::<f64>().sqrt();
println!("scores shape: {:?}", scores.shape());
println!("reconstruction shape: {:?}", reconstructed.shape());
println!("reconstruction error: {:.6}", err);
}
```
留意这里的形状约定，因为它和 `transform` 正好相反。`transform` 接收 `n_features` 列，返回 `n_components` 列。`inverse_transform` 接收 `n_components` 列，返回 `n_features` 列。传给它一个列数对不上 `n_components` 的矩阵，你会得到 `Error::DimensionMismatch`，而不是悄无声息的广播。
下一节要讲的符号修正，会把一整条主轴连同它的得分一起翻转。这让重建结果对符号翻转不敏感：乘积 `scores * components` 不会变，不论某条主轴有没有被取反。符号翻转永远不会破坏一次往返。
## 2.10.6. 确定性：符号约定与随机种子
特征向量和奇异向量只在相差一个符号的意义下唯一：`-v` 和 `v` 张成同一条主轴。scikit-learn 用户对此并不陌生，这正是 PC 符号有时会在不同运行或不同库版本间翻转的原因。
RustyML 在分解之后把符号钉死，让结果保持确定性：必要时对每个主成分行取反，使其绝对值最大的载荷变为非负。结果就是，3 种求解器在同一份数据上，对每一条主轴的朝向都能达成一致。同一个求解器的重复运行也会达成一致。符号保持稳定，`fit` 无需任何额外步骤就能复现。
这带来 2 个推论。其一，这套符号约定是 rustyml 自己特有的：它取决于主成分向量本身，而不是像 scikit-learn 的 `svd_flip` 那样取决于 `U` 因子。所以同一份数据下，某个主成分的符号可能和 scikit-learn 打印出来的不同。这个差异纯属外观：主轴、它解释的方差、每一次重建都完全相同。
其二，运行间唯一真正的变动来源是 `SVDSolver::Randomized`：它的随机草图由你传入的 `u64` 播种。同一个种子给出逐比特一致的输出，不同种子给出略有差异的近似子空间。`Full` 和 `PowerIteration` 都不带外部随机性，`PowerIteration` 用一个固定值在内部给自己的起始向量播种。
如果你在意整条流水线的精确可复现性，参见[可复现性与随机种子](../Chapter-07/7.1._可复现性与随机种子.md)。这里唯一的杠杆就是 `Randomized` 的种子。
拟合好的 `PCA` 用 `save_to_path` 和 `load_from_path` 序列化。这两个方法把整个模型（均值、主成分、方差、奇异值）写入或读出为紧凑的 [postcard](https://docs.rs/postcard) 二进制格式。文件扩展名无关紧要：不管你叫它什么，字节内容都是 postcard 格式。重新加载的模型，transform 出来的结果与原模型完全一致：
```rust,ignore
pca.save_to_path("pca_model.bin")?;
let loaded = PCA::load_from_path("pca_model.bin")?;
let scores = loaded.transform(&x_new)?; // 与保存前的模型完全一致
```
其中的机制，包括版本兼容方面的考量、二进制格式在什么情况下适合长期持久化，参见[深入模型持久化](../Chapter-07/7.2._深入模型持久化.md)。
## 2.10.7. 错误与边界情况
PCA 会校验输入，并返回带类型的[错误](../Chapter-01/1.6._错误处理.md)，而不是直接 panic。这些示例里的 `.unwrap()` 只是为了让代码简短。不要在生产代码里这样使用 `.unwrap()`。下表列出了实际会撞上的错误，外加一条 `PowerIteration` 求解器特有的错误：
| `PCA::new(0)` | `Error::InvalidParameter` |
| `fit` 时 `n_components > min(n_samples, n_features)` | `Error::InvalidParameter` |
| `fit` 时特征矩阵为空 | `Error::EmptyInput` |
| `fit` 时样本少于 2 个 | `Error::InvalidInput` |
| `fit` 或 `transform` 时输入含 `NaN` 或 `Inf` | `Error::NonFinite` |
| `fit` 之前调用 `transform` 或 `inverse_transform` | `Error::NotFitted` |
| `transform` 传入的特征数不对 | `Error::DimensionMismatch` |
| `inverse_transform` 传入的得分列数不对 | `Error::DimensionMismatch` |
| `SVDSolver::PowerIteration` 收敛失败，例如数据的有效秩低于 `n_components` 时 | `Error::NotConverged` |
`NotConverged` 这种情况是 `PowerIteration` 特有的。`Full` 和 `Randomized` 都不会抛出它，因为它们都不依赖严格的逐个主成分收敛检查。
2 个样本的下限并非随意定的：方差需要 `n - 1` 作分母，单独一行没有什么可分解的。`n_components` 的上限 `min(n_samples, n_features)` 是秩的界，你没法提取出比数据所张成的还要多的正交方向。
要得更多会引发参数错误，而不是被悄悄截断。rustyml 是对着数据检查这条界限的，而不是只对着模型本身。所以同一个 `PCA::new(5)` 实例，在 100 x 20 矩阵上能成功，在 3 x 4 矩阵上就会失败。
PCA 生来就是线性的：它只能找出输入特征的线性组合方向。盘踞在弯曲流形上的结构，对它而言始终是隐形的。当碎石图怎么都不肯趋平、重建在每个 `k` 上都很差时，这通常就是该转向非线性方法的信号。
[核主成分分析](./2.11._核主成分分析.md)把同一套方差最大化的思路搬进核特征空间。[t-SNE](./2.12._t-SNE.md)是把非线性邻域结构做二维可视化的工具。如果目标是分开带标签的类别，而不是捕捉原始方差，就用它的有监督对应物[线性判别分析](./2.6._线性判别分析.md)：它把数据朝类别可分的方向投影，而不是朝总体散布的方向。