280 MiB 的矩阵做 SVD 用了 1.5 GiB

主线位置：直接计算主线，TypedTensor::svd(&mut engine) 落到 tlinalg 和 faer。线性代数怎么落到 faer 见 线性代数；这块内存逐份花在哪里、本站 705 × 52 120 的真实数据上的完整对比，见 用 tenferro 做潜在语义分析。

flowchart TB
  a["A 280 MiB"] --> svd["svd"]
  svd --> work["faer 工作区 V 和 QR"]
  svd --> res["tlinalg 残差矩阵"]
  svd --> vt["Vt 输出"]
  work --> peak["峰值约 5.4 倍 A"]
  res --> peak
  vt --> peak
  a --> g["G = A A^T"]
  g --> eigh["eigh"]
  eigh --> small["峰值约 1.3 倍 A"]

实验 crate exp-tenferro8/pit2-exp，bin p24_svd_peak，release，CpuEngine::new 一个 8 线程的池，f64，M2 Max、macOS 27、rustc 1.99.0。每种做法一个进程，套在 /usr/bin/time -l 下读「peak memory footprint」，各跑 3 次取最快的那次的时间；开始时负载 2.7。矩阵是确定性的伪随机稠密矩阵，不是本站的数据。

现象

一个 705 × 52 120 的 f64 矩阵（列主序 280.3 MiB），只是想要它的前几个奇异值和左奇异向量：

let mut engine = CpuEngine::new(Pool::builder().threads(8).build()?);
let a = TypedTensor::<f64>::from_vec_col_major(vec![705, 52_120], data)?;
let (u, s, vt) = a.svd(&mut engine)?;

input 705 x 52120 f64 = 280.3 MiB
svd: u [705, 705] s [705] vt [705, 52120], s[0] = 73.497600, 1333 ms
        1.43 real         8.88 user         0.16 sys
          1613070336  maximum resident set size
...
          1610647280  peak memory footprint

进程峰值 1 610 647 280 字节，即 1536 MiB；只建矩阵、不做分解的进程是 282.3 MiB。分解本身多用了约 1254 MiB，是输入的 4.5 倍。三次运行的峰值相差不到 0.3%。

同一个矩阵的另外两种做法：

eigh 给出的最大奇异值 sqrt(max w) = 73.497600，和 svd 的 s[0] 相同。

出处

svd_impl 从 engine 的缓冲池取 U、S、Vt 三个输出和一块工作区，交给 tlinalg 的 svd_into。同时活着的大块有：输入 A；faer 对扁长矩阵先做 QR 时复制进工作区的一份；faer 在工作区里填的 52 120 × 705 的 V；tlinalg 逐元素转置写出的 Vt；较小一维不小于 128 时 tlinalg 为了残差检查 A − U diag(S) Vᴴ 分配的 705 × 52 120 残差矩阵。每份都是 280 MiB，五份加起来就是实测的 1.5 GiB。逐项的源码位置和为什么要做残差检查（faer 的分治 SVD 在某些输入上会给出不准确的因子而不报错），见 用 tenferro 做潜在语义分析 的「截断 SVD：tenferro 能给什么」。

两点让它躲不开：

• tenferro-linalg 没有截断 SVD，也没有「只要 U」的入口。svd 一定返回完整的 Vt，而 Vt 恰好是 52 120 那一侧。
• traced 程序的输出剪枝只认「Svd 只用到 S」这一种，改写成 SvdVals；只用 U 和 S、丢掉 Vt，不在其中（同一页的实测）。

svdvals 一律走 QR 算法，没有残差检查、不要 V，峰值降到 612 MiB；比只建矩阵多出的 329 MiB 约是一份矩阵加 faer 的工作区。

修法

只要左奇异向量和奇异值时，算 Gram 矩阵 G = A Aᵀ（705 × 705，3.8 MiB），对它做 eigh：G 的特征值是 σ²，特征向量就是 U，52 120 那一侧根本不出现。用 dot_general 直接收缩两边的第 1 维，不要先 transpose 出一份 Aᵀ：

let cfg = DotGeneralConfig {
  lhs_contracting_dims: [1].as_slice().into(),
  rhs_contracting_dims: [1].as_slice().into(),
  lhs_batch_dims: [].as_slice().into(),
  rhs_batch_dims: [].as_slice().into(),
};
let g = a.dot_general(&a, cfg, &mut engine)?;
let (w, v) = g.eigh(&mut engine)?;   // w 升序，σ = sqrt(w)，v 的列是 U

input 705 x 52120 f64 = 280.3 MiB
eigh(A A^T): w [705] v [705, 705], sqrt(max w) = 73.497600, 201 ms
        0.29 real         1.83 user         0.03 sys
...
           386286456  peak memory footprint

• 比 svd 快 6.6 倍，分解部分多用的内存从 1254 MiB 降到 86 MiB。先 a.transpose(&[1, 0], &mut engine) 再 matmul 的写法多一份 280 MiB 的 Aᵀ（第一次试验里峰值 649 MiB）。
• 代价是精度：小奇异值经过平方，相对误差放大。只要前 k 个、谱的动态范围不大时（LSA、PCA 这类用法）没有问题；要很小的奇异值就不能这样做。
• 两侧都要（U 和 V）时只能用 svd，这时按 1.5 GiB 准备内存；或把矩阵按行分块、自己写随机化 SVD。

怎么确认没中招

• 对大矩阵调 svd 之前算一下：m × n × 8 字节乘 5，就是要准备的峰值（较小一维小于 128 时没有残差矩阵，少一份）。
• 只用到 U 和 S 的代码，看 vt 是不是一直被丢掉（let (u, s, _) = ...）；是，就换成 Gram 加 eigh。
• 在 /usr/bin/time -l（macOS）或 /usr/bin/time -v（Linux）下跑一次，看峰值。

小结

• svd 总是返回完整的薄 U、S、Vt；705 × 52 120 的 f64 矩阵（280 MiB）一次 svd 让进程峰值到 1536 MiB，分解多用的是输入的 4.5 倍。
• 来源是 faer 的 QR 复制和 V、tlinalg 的残差矩阵、转置出来的 Vt 输出，各 280 MiB；没有截断或只要 U 的入口。
• 只要 U 和 S 时用 G = A Aᵀ（dot_general，不转置）加 eigh：201 ms、峰值 368 MiB，最大奇异值和 svd 的一致；代价是小奇异值的精度。