分解的梯度
主线位置：TRG 主线「一步 SVD 截断加收缩」之前，先确认 svd、eigh、qr、solve 这几个分解在两边求导的结果一致。
版本：JAX 0.12.0.1（CPU，x64）、tenferro-rs 0.8.0（CPU，native，单线程）
flowchart TB a["矩阵 A"] --> j["jnp.linalg"] a --> t["TracedTensorLinalgExt"] j --> lj["与符号无关的损失"] t --> lt["与符号无关的损失"] lj --> gj["jax.grad"] lt --> gt["AdContext grad"] gj --> fd["对照中心差分"] gt --> fd
两边的分解都是薄分解，求导规则也是同一套公式，tenferro 把系数 F 展宽了，只在奇异值相等或极接近时结果和 JAX 不同（见 简并奇异值）；它的 linearize_svd 连 JAX 的复数规范修正项都照搬了（规则本身见 线性代数 的「SVD 的反向规则」）。所以在普通矩阵上，梯度两边相同到 1e-15 量级。要注意的是另外几件事：分解得到的向量符号不同；测试矩阵别选秩亏的；solve 的右端项形状。
tenferro 一侧的求导要用带 tenferro_linalg::semantic_ad_rules() 的 AdContext，TracedTensorAdExt::grad 的规则集是空的，经过 SVD 会报错，见 grad、vjp、jvp 的「AdContext 和 TracedTensorAdExt」和常见坑 unsupported vjp AD rule。
实验代码：JAX 一侧是 grads.py，tenferro 一侧是 bin grads（release，CpuEngine::serial()），对照脚本 compare_grads.py。M2 Max、macOS 27、rustc 1.99.0。
函数对照
| JAX | tenferro（TracedTensorLinalgExt） | 说明 |
|---|---|---|
jnp.linalg.svd(a, full_matrices=False) → u, s, vh | a.svd()? → (u, s, vt) | 都是薄 SVD，s 降序；jnp.linalg.svd(a) 默认是满的，见 SVD 的输出 |
jnp.linalg.eigh(a) → w, v | a.eigh()? → (w, v) | 特征值升序，v 的列是特征向量 |
jnp.linalg.qr(a) → q, r | a.qr()? → (q, r) | 默认都是 reduced：6×4 的 a 得到 6×4 的 q、4×4 的 r |
jnp.linalg.solve(a, b) | a.solve(&b)? | tenferro 的 b 必须是二维，见下文 |
eps 一类的参数没有 | svd_with_options(SvdOptions)、eigh_with_options、qr_with_options | 求导用的 derivative_eps 和符号规范 |
四个损失，两边对差分
损失和试验阶段的 linalg_grads 一样：对每个分解构造一个只依赖「和符号无关」的量的标量，比如奇异值、投影矩阵 U₂U₂ᵀ、diag(R)²。矩阵元素是 sin(1 + 0.7i² + 2.3j + 0.5ij²)，权重是 cos(i − 2j)。
def loss_svd(a):
u, s, vh = jnp.linalg.svd(a, full_matrices=False)
u2 = u[:, :2]
return jnp.sum(jnp.array([1.0, 0.5, -0.3, 0.2]) * s) + jnp.sum(wmat(6, 6) * (u2 @ u2.T))
g = jax.jit(jax.grad(loss_svd))(a) # 再和步长 1e-6 的中心差分比
fn loss_svd(a: &TracedTensor) -> R<TracedTensor> {
let (u, s, _vt) = a.svd()?;
let u2 = (&u * &colmask(6, 4, &[0, 1])?)?; // u[:, :2] 写成掩码，原因见「SVD 的输出」
Ok((&wsum(&s, &vecc(&[1.0, 0.5, -0.3, 0.2])?)? + &wsum(&outer(&u2)?, &wmat(6, 6)?)?)?)
}
let g = ad.grad(&loss_svd(&x)?, &x)?; // ad = AdContext + linalg 规则；编译、准备、执行后读回
eigh 的损失是 A + Aᵀ 的特征值加权和加最大特征向量的投影；qr 是 QQᵀ 和 diag(R)² 的加权和；solve 是 (A + 3I)x = b 的解的加权和。
# JAX
svd loss=+3.115374357904 max|grad|=7.877e-01 max|grad-fd|=1.556e-09 nan=0
eigh loss=+2.571532979194 max|grad|=1.480e+00 max|grad-fd|=1.222e-09 nan=0
qr loss=+3.339249699594 max|grad|=4.612e+00 max|grad-fd|=1.279e-09 nan=0
solve loss=-0.500244841286 max|grad|=1.517e-01 max|grad-fd|=1.526e-10 nan=0
# tenferro
svd loss=+3.115374357904 max|grad|=7.877e-1 max|grad-fd|=2.082e-9 nan=0
eigh loss=+2.571532979194 max|grad|=1.480e0 max|grad-fd|=1.295e-9 nan=0
qr loss=+3.339249699594 max|grad|=4.612e0 max|grad-fd|=1.388e-9 nan=0
solve loss=-0.500244841286 max|grad|=1.517e-1 max|grad-fd|=1.970e-10 nan=0
# 两边的梯度逐元素比
grad svd max|jax - tenferro| = 2.359e-15 max|grad| = 0.788
grad eigh max|jax - tenferro| = 5.204e-16 max|grad| = 1.480
grad qr max|jax - tenferro| = 2.665e-15 max|grad| = 4.612
grad solve max|jax - tenferro| = 2.776e-17 max|grad| = 0.152
- 损失值 12 位相同，对差分的误差都在 1e-9 量级。这个量级来自差分本身：步长 1e-6 的中心差分，截断误差和舍入误差加起来就是这么大，两边差不多。
- 两边的梯度相差 1e-15 量级，相对梯度的大小（0.15 到 4.6）是几个 ulp。分解由不同的库算（JAX 走 LAPACK，tenferro 走 faer），反向里的矩阵乘顺序也不同，所以不逐位相同。
- 奇异值、特征值、
R两边差 1.8e-15、1.8e-15、3.3e-16。
奇异向量的符号不同
同一个矩阵的分解因子直接相减：
svd_s max|jax - tenferro| = 1.776e-15
eigh_w max|jax - tenferro| = 1.776e-15
qr_r max|jax - tenferro| = 3.331e-16 sign per vector = [1, 1, 1, 1] after sign fix 3.331e-16
svd_u max|jax - tenferro| = 1.253e+00 sign per vector = [1, 1, -1, 1] after sign fix 4.441e-16
svd_vt max|jax - tenferro| = 1.586e+00 sign per vector = [1, 1, -1, 1] after sign fix 6.661e-16
qr_q max|jax - tenferro| = 3.331e-16 sign per vector = [1, 1, 1, 1] after sign fix 3.331e-16
eigh_v max|jax - tenferro| = 1.580e+00 sign per vector = [-1, -1, -1, -1, -1] after sign fix 3.469e-16
sign per vector 是两边对应向量内积的符号。SVD 的第 3 个奇异向量方向相反（U 的列和 Vᵀ 的行一起反，U S Vᵀ 不变）；eigh 的 5 个特征向量全反；QR 两边一样，diag(R) 都是 [-1.742, 1.692, -1.585, -1.010]。把符号对齐后，差都回到 1e-16 量级。
奇异向量和特征向量本来就只确定到一个符号，哪个符号由算法的细节决定，LAPACK 和 faer 各有各的。这对搬代码有两个后果：
- 损失要和符号无关，比如投影
uuᵀ、奇异值、|u|²。JAX 的测试里写u[:, 0]的加权和，在 tenferro 里梯度可能整个反号，差分也一样跟着反，两边各自是「对的」，但不能互相比。 - 把 JAX 算出的向量当参考值存下来、再和 tenferro 比，要先对齐符号。
tenferro 有一个 JAX 没有的选项：SvdOptions::gauge(SvdGauge::CanonicalPivot)，把 U 每列绝对值最大的元素变成正数，并把同一个符号乘到 Vᵀ 的对应行上（SvdGauge）。EighOptions、QrOptions 也有类似的选项。JAX 一侧只能手工做同样的规范：
let (uc, _sc, vtc) = a.svd_with_options(SvdOptions::default().gauge(SvdGauge::CanonicalPivot))?;
piv = np.sign(u[np.argmax(np.abs(u), axis=0), range(4)]) # JAX 的 U 每列主元的符号
u_canon, vt_canon = u * piv, vh * piv[:, None]
svd_u CanonicalPivot vs JAX U with pivot made positive: 4.441e-16 pivot signs of JAX U [1, -1, -1, 1]
svd_vt CanonicalPivot vs JAX Vt with the same row signs: 6.661e-16
规范之后两边的 U、Vᵀ 差 1e-16 量级。默认（SvdGauge::Raw）不做这一步。
PR #1338 引入 SvdOptions 和 SvdGauge::CanonicalPivot，同一天的 PR #1339（「Add opt-in linalg gauge conventions」）补上 EighGauge、QrGauge 和 qr_with_options，两者都首次发布于 v0.3.0。之前三个分解只有后端给的原始符号。#1339 随附的工作记录写明了取舍：「Kept gauge fixing opt-in. svd, eigh, and qr still use backend raw signs/phases by default.」以及「Did not add new gauge-aware AD rules. The existing AD formulas still operate on the chosen forward outputs」。规范只是前向输出的后处理，求导规则不变；所以上面四个损失的梯度不受这个选项影响，而一个依赖符号的损失在 Raw 和 CanonicalPivot 下会得到不同的值。
秩亏的测试矩阵
试验阶段记下一条：秩亏的测试矩阵会让 JAX 的 svd、qr 梯度和差分对不上。重跑时两边一起测。把 6×4 矩阵的第 4 列换成第 1、2 列之和（浮点相加，秩在舍入意义上是 3）：
dup = base.copy()
dup[:, 3] = dup[:, 0] + dup[:, 1]
# JAX
dup svd s = [2.54062902e+00 2.17416489e+00 1.46539222e+00 1.21167608e-16] diag(R) = [-1.74204082e+00 1.69209801e+00 -1.58541566e+00 -3.81934349e-16]
svd dup loss=+3.524811646420 max|grad|=1.209e+00 max|grad-fd|=6.858e-02 nan=0
qr dup loss=+2.823253820524 max|grad|=1.620e+15 max|grad-fd|=1.620e+15 nan=0
# tenferro
dup svd s = ["2.5406290203799005e0", "2.17416488595396e0", "1.4653922170352183e0", "2.3602035230609313e-16"] diag(R) = ["-1.7420408237631937e0", "1.6920980060104915e0", "-1.5854156582908705e0", "-2.775555224048112e-16"]
svd dup loss=+3.524811646420 max|grad|=1.209e0 max|grad-fd|=6.860e-2 nan=0
qr dup loss=+2.474531276122 max|grad|=3.476e15 max|grad-fd|=3.476e15 nan=0
两边都对不上，而且对不上的方式一样。问题出在测试点，不在求导：
- svd：损失里有
0.2 · s₄，s₄在 0 附近像|x|，在 0 处不可导。差分取的是两侧斜率的平均，梯度取的是其中一侧，于是差 7e-2。两边的梯度差同样多，对差分的误差 6.858e-2 和 6.860e-2。 - qr：
Q的第 4 列在秩亏处没有定义，算出来是什么由舍入决定。连损失值都不一样：LAPACK 给 2.823，faer 给 2.475。梯度含R⁻¹，r₄₄是 1e-16 量级，梯度就到 1e15。
所以测试分解的梯度，矩阵要随机生成或者像上面那样用解析式构造，并检查最小奇异值离 0 和彼此之间都足够远；想测秩亏或简并的行为，就单独测，并且损失只用定义良好的量（见 简并奇异值）。
s₄ 就恰好是 0。svd 损失的 24 个梯度分量里，JAX 和 tenferro 各有几个 NaN？提示：6×4 是 m > n。tenferro 卷的常见坑 SVD 的梯度全是 NaN 讲了 tenferro 一侧的机制，那里把它记作上游 bug。JAX 会不会不一样？
solve 的右端项
jnp.linalg.solve(a, b) 接受一维的 b，当成一列。tenferro 不接受：
let b1 = TracedTensor::from_vec_col_major(vec![5], (0..5).map(|i| (i as f64).cos()).collect())?;
a5p.solve(&b1)
solve(a, b[5]): Err(tenferro-linalg.lu_factor_solve: rank mismatch: expected 2, actual 1)
错误在建图那一行就返回。右端项写成 [5, 1]，解也是 [5, 1]；上面 solve 的损失就是这样写的，权重也是 [5, 1]。
动脑筋答案
# JAX
zero svd s = [2.3266923 1.83134488 1.09677589 0. ] diag(R) = [-1.74204082 1.69209801 -1.58541566 0. ]
svd zero loss=+2.308468906516 max|grad|=nan max|grad-fd|=nan nan=24
qr zero loss=+1.120156703074 max|grad|=nan max|grad-fd|=nan nan=18
# tenferro
zero svd s = ["2.3266923038930245e0", "1.8313448767848501e0", "1.0967758874392364e0", "0e0"] diag(R) = ["-1.7420408237631937e0", "1.6920980060104915e0", "-1.5854156582908705e0", "0e0"]
svd zero loss=+2.308468906516 max|grad|=0.000e0 max|grad-fd|=0.000e0 nan=24
qr zero loss=+1.120156703074 max|grad|=inf max|grad-fd|=inf nan=18
两边都是 24 个里 24 个 NaN。tenferro 一行的
max|grad|=0.000e0是打印的问题：fold(0.0, f64::max)遇到 NaN 时取另一边，NaN 全被跳过了，要看nan=24。原因两边一样。JAX 0.12.0.1 的
_svd_jvp_rule（jax/_src/lax/linalg.py）对 m > n 加的投影修正项直接除以s：if m > n: dAV = dA @ V dU = dU + (dAV - U @ (Ut @ dAV)) / s_dim.astype(A.dtype)同一个规则里，
1/s的对角项用s_inv = 1 / (s + s_zeros) - s_zeros挡住了s = 0，这一行没有挡。tenferro 的 这一行 也直接除以原始的s。s₄ = 0时这一列是0/0，再经矩阵乘扩散到全部 24 个分量。- 所以 tenferro 卷记下的这个坑，JAX 0.12.0.1 也有。绕法也通用：补零成方阵，或者把恰好为 0 的奇异值扰动开（见那一页）。
qr 两边都是 18 个 NaN，tenferro 剩下的分量里还有
inf。损失用了整个QQᵀ，r₄₄ = 0时这个函数本身不可导。
小结
- svd、eigh、qr、solve 的梯度两边对中心差分都在 1e-9 量级，彼此差 1e-15 量级；两边用的是同一套公式，不同的库和运算顺序造成末几位不同。
- 奇异向量、特征向量的符号两边不同（这个矩阵上 SVD 反了 1 个，eigh 反了 5 个，QR 相同）。损失写成和符号无关的形式；要比较向量本身，先对齐符号，或者用
SvdGauge::CanonicalPivot再在 JAX 一侧手工做同样的规范。 - 秩亏的测试矩阵两边都对不上差分：最小奇异值在 0 处是尖点，
Q的多余列没有定义。恰好为 0 的奇异值在 m > n 时两边都给全 NaN，JAX 的规则和 tenferro 一样在投影修正项里除以s。 - tenferro 的
solve要求b是二维，一维的b在建图时报rank mismatch: expected 2, actual 1。