From 2764a180757d8434937d2bfb9e14b804e6a1c5fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Petr=20Balv=C3=ADn?= Date: Sun, 27 Sep 2026 11:58:07 +0200 Subject: [PATCH] fix(linalg): converge the Hermitian Jacobi sweep on near-diagonal matrices Assisted-by: GLM 5.3 Flash --- CHANGELOG.md | 4 ++++ linalg/decomp3.go | 13 ++++++++++++- linalg/decomp3_test.go | 33 +++++++++++++++++++++++++++++++++ 3 files changed, 49 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e34e3e..6ce8330 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 update refuse a matrix or vector holding a NaN or infinity instead of returning an all-NaN factor or spectrum with a nil error, the refusal the sparse solvers already make. +- `EigenComplex` converges on a Hermitian matrix whose off-diagonal + entries all sit just under the per-entry skip level: the sweep no + longer exhausts its passes on a matrix it should have declared + converged. ## [1.0.0] - 2026-09-03 diff --git a/linalg/decomp3.go b/linalg/decomp3.go index 3714860..61fa043 100644 --- a/linalg/decomp3.go +++ b/linalg/decomp3.go @@ -122,6 +122,17 @@ func hermitianJacobi(h []complex128, v []complex128, n int) error { // the floor already diagonal and return its raw diagonal, identity // eigenvectors included. threshold := 1e-13 * norm + // The per-entry skip must sit where a matrix whose every + // off-diagonal entry is skippable also passes the aggregate + // convergence test: with P pairs the off-norm of all-skipped entries + // reaches skip·√P, so a skip at the plain threshold leaves a matrix + // whose entries are uniformly just under it frozen above the test, + // and the sweep exhausts its passes without turning a single + // rotation. + skip := threshold + if pairs := float64(n * (n - 1) / 2); pairs > 1 { + skip = threshold / math.Sqrt(pairs) + } offMass := func() float64 { off := 0.0 for p := range n { @@ -139,7 +150,7 @@ func hermitianJacobi(h []complex128, v []complex128, n int) error { for p := range n { for q := p + 1; q < n; q++ { z := h[p*n+q] - if base.AbsComplex(z) <= threshold { + if base.AbsComplex(z) <= skip { continue } // Step 1: a diagonal phase turn makes the pair entry diff --git a/linalg/decomp3_test.go b/linalg/decomp3_test.go index e9e4ed6..27dd88a 100644 --- a/linalg/decomp3_test.go +++ b/linalg/decomp3_test.go @@ -294,6 +294,39 @@ func subComplex(a, b []complex128) []complex128 { return out } +// TestEigenComplexNearDiagonalConverges pins the Jacobi sweep against a +// matrix whose every off-diagonal entry sits just under the per-entry +// skip level while the aggregate off-norm stays above the convergence +// threshold: the skip must not freeze the sweep above its own +// convergence test, which used to exhaust the passes and report the +// nearly diagonal matrix as unconverged. +func TestEigenComplexNearDiagonalConverges(t *testing.T) { + const n = 30 + cv := make([]complex128, n*n) + for i := range n { + cv[i*n+i] = 1 + } + for i := range n { + for j := i + 1; j < n; j++ { + cv[i*n+j] = complex(0.9e-13, 0) + cv[j*n+i] = complex(0.9e-13, 0) + } + } + a, err := core.FromComplexes(cv, n, n) + if err != nil { + t.Fatalf("FromComplexes: %v", err) + } + vals, _, err := EigenComplex(a) + if err != nil { + t.Fatalf("EigenComplex: %v", err) + } + for i := range n { + if math.Abs(vals.FloatAt(i)-1) > 1e-11 { + t.Fatalf("eigenvalue %d is %.12g, want 1", i, vals.FloatAt(i)) + } + } +} + // TestEigenComplexRefusesNonFinite pins that a poisoned matrix never // reads as Hermitian: the mirror comparison cannot see a NaN // difference, so the entry is refused outright, the way the sparse