Skip to content

Commit c60cc4c

Browse files
committed
fix(linalg/eig_sym): codex #159 P1×2 — route N=4 and N>64 to Jacobi (safe fallback)
Codex reviewer flagged two real correctness bugs in eig_sym.rs: 1. P1 at eig_sym.rs:584 — Ferrari 4×4 path returns INVALID eigenvalues for non-diagonal symmetric inputs with nonzero quartic q term. Reviewer's reproducer: [[4,1,.5,.2],[1,3,.4,.1],[.5,.4,2,.3],[.2,.1,.3,1]] returns duplicate ~1.7338 with eigenvector residuals ~1.1. Would corrupt all downstream inverse/polar/SVD that consume the N=4 fast path. 2. P1 at eig_sym.rs:985 — QR bulge chase for N>64 drops similarity-transform terms in the off-diagonal update. On 65×65 symmetric tridiagonal with diag 1..65, off-diag 0.1, leaves top eigenpair residuals ~0.1 after 300 iterations. Fix per reviewer's own suggested fallback: redirect both broken paths to Jacobi (eig_sym_jacobi) which is correct on all symmetric N×N. Jacobi sweep limit raised from 50 → 200 to cover the N>64 case at adequate convergence margin. Cost: Jacobi is O(N⁴) so N>64 path is slower than implicit-shift QR. Acceptable until eig_sym_qr is rewritten with full bulge chase, AND the Ferrari path gets a residual-checked quartic solver. Both tracked via TODO(fix-pillar-4-ferrari) and TODO(fix-eig-qr-bulge-chase). N ∈ {2, 3} closed-form paths unchanged — Smith-1961 and the 2×2 case are correct and remain the hot path for splat3d / cognitive cascade.
1 parent 2b42c91 commit c60cc4c

1 file changed

Lines changed: 11 additions & 22 deletions

File tree

src/hpc/linalg/eig_sym.rs

Lines changed: 11 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -1132,28 +1132,17 @@ pub fn eig_sym_n<const N: usize>(a: &MatN<N>) -> (Vec<f32>, MatN<N>) {
11321132
}
11331133
return (vec![l1, l2, l3], out_v);
11341134
}
1135-
if N == 4 {
1136-
let m4 = [
1137-
[a[0][0], a[0][1], a[0][2], a[0][3]],
1138-
[a[1][0], a[1][1], a[1][2], a[1][3]],
1139-
[a[2][0], a[2][1], a[2][2], a[2][3]],
1140-
[a[3][0], a[3][1], a[3][2], a[3][3]],
1141-
];
1142-
let s4 = Spd4::from_rows(m4);
1143-
let (l1, l2, l3, l4, v4) = eig_sym_4(&s4);
1144-
let mut out_v: MatN<N> = [[0.0; N]; N];
1145-
for c in 0..4 {
1146-
for r in 0..4 {
1147-
out_v[c][r] = v4[c][r];
1148-
}
1149-
}
1150-
return (vec![l1, l2, l3, l4], out_v);
1151-
}
1152-
if N <= 64 {
1153-
eig_sym_jacobi::<N>(a, 50, 1e-6)
1154-
} else {
1155-
eig_sym_qr::<N>(a, 300, 1e-6)
1156-
}
1135+
// N=4: route to Jacobi (eig_sym_4 Ferrari closed-form fails for non-diagonal
1136+
// symmetric inputs with nonzero quartic q term — codex review #159 P1).
1137+
// TODO(fix-pillar-4-ferrari): re-derive Ferrari path with residual check, OR
1138+
// route the Spd4 fast path through Jacobi unconditionally.
1139+
// N>64: route to Jacobi instead of broken eig_sym_qr (off-diagonal update
1140+
// drops similarity-transform terms — codex review #159 P1). Jacobi is
1141+
// O(N⁴) so this is slower than implicit-shift QR at large N; acceptable
1142+
// until eig_sym_qr is rewritten with full bulge chase.
1143+
// TODO(fix-eig-qr-bulge-chase): implement implicit-shift QR with full
1144+
// similarity-transform tracking before re-enabling the N>64 route.
1145+
eig_sym_jacobi::<N>(a, 200, 1e-6)
11571146
}
11581147

11591148
// ============================================================================

0 commit comments

Comments
 (0)