Repository navigation
fix(gevd): rotate the columns of Q when the complex QZ chases a zero down B (#324) - #338
Conversation
…down B (#324) hessenberg_to_qz_unblocked (the complex QZ path for n < 15) chases a negligible T[(j, j)] at j == ilo, or below a negligible H[(j, j - 1)], to the bottom with rotations on rows jch, jch + 1 of H and T. Q has to take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path already do (LAPACK zhgeqz rotates Q(1,JCH) and Q(1,JCH+1) with C, DCONJG(S)). It rotated rows jch, jch + 1 of Q with s. Q stayed unitary but A = Q S Z^H no longer held, so every left eigenvector of gevd_cplx was wrong, while the eigenvalues and the right eigenvectors, which only use S, T and Z, were right. Whether a pencil reaches the branch depends on rounding: the issue's 12 x 12 c32 pencil does on x86-64 with AVX2/FMA (left residual 1.962e-1 -> 2.079e-7) and does not on aarch64. Tests: a unit test hands the unblocked iteration an exact zero on B's diagonal (j == ilo and below a zero subdiagonal, c32/c64) and checks the reconstruction; an integration test embeds the issue's bit patterns and checks left and right residuals at c32 and c64. The same lines are in upstream faer (Codeberg main 5833f4e). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
GPT-6.1-sol adversarial reviewCRITICAL: none Critical (must not merge)none Needs changes
Unverified claims
Checked and finenone Not checked: the full local suite or the original cross-architecture diagnostic runs; the mutation probe establishes a test-oracle failure, not a production failure of the corrected rotation. Reviewed head: |
… fix the size claim `residuals()` in `faer/tests/gevd_infinite_eigenvalue.rs` folded each column into its maxima with `f64::max`, which returns the other operand when one is `NaN`. An all-zero left eigenvector column gives `0 / 0`, so the `NaN` was dropped and the column read as perfectly accurate: with `ul` forced to zeros the test passed. Each column now asserts, before it is folded in, that both residual numerators and `|beta| / (|alpha| + |beta|)` are finite and that `||u||`, `||v||` and the residual denominator are non-zero, reporting the column index and values. The same injection now fails at `column 0: left eigenvector has norm 0e0`. The `n < 15` size claim was wrong. `hessenberg_to_qz_blocked` calls `hessenberg_to_qz_unblocked` and returns when `n < nmin`, and calls it again at the end on the whole range when it does not, so the branch is reachable at every `n`. `nmin` is `Ord::max(15, blocking_threshold)` and the default `blocking_threshold` is 75, so with the default parameters the whole iteration is unblocked for `n < 75`. `CHANGELOG.md`, `PATCHSET.md` and the test's module comment now say that. No change to the rotation logic.
GPT-6.1-sol adversarial reviewCRITICAL: none Critical (must not merge)none Needs changes
Unverified claims
Checked and fine
Not checked: the full test suite or local x86-64/AVX2 execution; production source files were not modified. Reviewed head: |
`left` and `right` are accumulated with `f64::max`, which returns the other operand when one is `NaN`; a `NaN` therefore leaves the maximum where it was and reads as a perfect residual. Each column now computes its two residuals, asserts both finite with the numerator, norm and denominator in the message, and only then folds them in. `ul = Mat::<T>::zeros(n, n)` after the `gevd_cplx` call turns the left assert red and the mirror mutation on `ur` turns the right one red; both were green before.
The residual was `numerator / (||v|| * denominator)`. For a column whose eigenvector norm is large enough the product overflows to infinity and the quotient is 0, which passes `is_finite` and reads as a perfect residual. Replacing a column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, left the test green with `nu * den = inf` and the residual recorded as `0e0` against a staged `2.978e-1`. Two divisions `(numerator / ||v||) / denominator` make the same mutation fail at `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5`.
GPT-6.1-sol adversarial reviewCRITICAL: none Critical (must not merge)none Needs changes
Unverified claims
Checked and fine
Reviewed head: |
`residuals()` built each column's numerator out of the eigenvector and the homogeneous pair as they come. The residual is invariant under both scalings, so the value was right whenever it survived, but two columns made it read as a perfect residual instead: - `alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0`: the denominator `|beta| ||A|| + |alpha| ||B||` overflows to `inf`, and `x / inf` is `0`; - `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`: the numerator underflows to `0` before the division. Both passed every guard and both are wrong by `O(1)` once normalized (`2.977833e-1` left, `2.941382e-1` right). Dividing afterwards cannot recover information the numerator never had, so the eigenvector is now divided by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)` before the residual is formed, the two norms, that scale and the denominator are required to be finite and positive, and the residual is one division again because `||u|| = 1` and `|alpha|, |beta| <= 1` make the numerator at most the denominator. The two columns are now the test `residual_oracle_rejects_unnormalized_eigenpairs`, which the unnormalized oracle fails with `left residual 0e0 < 1.43e-5` and `right residual 0e0 < 1.43e-5` for both. The columns of the issue's pencil read the same as before (`c32` 1.772195055005881e-7 / 1.3466182542302985e-7, `c64` 3.838923294839156e-16 / 2.9362584762522423e-16 against `3.864850603564751e-16` / `2.9009055001556393e-16` unnormalized; tolerance 2.7e-14).
GPT-6.1-sol adversarial reviewCRITICAL: none Critical (must not merge)none Needs changesnone Unverified claims
Checked and fine
Not checked: local Rust tests and mutations; Reviewed head: |
Fixes #324.
hessenberg_to_qz_unblockedhas a branch for a negligibleT[(j, j)]atj == iloor below a negligibleH[(j, j - 1)], which is an infinite eigenvalue. It chases that zero to the bottom with rotations on rowsjch, jch + 1ofHandT. BecauseA = Q H Z^H, the rotationGon the rows has to be undone asQ <- Q G^H, on columnsjch, jch + 1ofQwithconj(s). The QZ sweep (qz_cplx/mod.rs:371), theother zero-chase branch (
:169,trot(c, -s)), the blocked path (:1347,:1375) and LAPACKzhgeqz(
ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))) all do this. Onlyqz_cplx/mod.rs:123-127rotated rowsjch, jch + 1ofQwiths.Qstayed unitary, butA = Q S Z^Hfailed byO(1), so every left eigenvectorQ ywas wrong. The eigenvalues andthe right eigenvectors were still right, because they use only
S,TandZ. The real QZ (qz_real) rotatescolumns in the same branch and is not affected.
Which sizes reach the branch
hessenberg_to_qz_unblockedis not a "small matrix" path.hessenberg_to_qz_blockedruns it and returns whenn < nmin, and also calls it at the end on the wholeilo..ihirange when it does not, so the branch is reachable atevery
n.nminisOrd::max(15, params.blocking_threshold)(qz_cplx/mod.rs:1225), and the defaultblocking_thresholdis75(linalg/evd/schur/mod.rs:31, propagated throughGeneralizedSchurParamsinlinalg/gevd/mod.rs:112). With the default parameters, then, the whole iteration takes the unblocked path forn < max(15, blocking_threshold), i.e.n < 75, and above that the blocked loop runs first and the same rotationsrun in the final unblocked call. An earlier revision of this description, of
CHANGELOG.md, of thePATCHSET.mdrow and of the test's module comment stated a hard 15-element cutoff; that was wrong.
Measured rather than read, with a temporary
eprintln!at both call sites inhessenberg_to_qz_blockedand thedefault
params(a 20 x 20 and a 100 x 100 pencil throughgevd_cplx):So a 20 x 20 call is unblocked by default, not only under a 15-element cutoff, and a 100 x 100 call still ends in the same code. The
eprintln!is not part of this branch.Why only c32, and only on some machines
Whether a pencil reaches the branch depends on rounding in the reduction and the iteration. Dumps of
A,B,Qand
Zat each stage ofcompute_gevd_generic, for the issue's pencil, residuals relative to||A0||and||B0||:c32)Q^H A0 Z - AQ^H B0 Z - BQunitaryBhessenberg_to_qzA least-squares fit shows that
Zis consistent with(S, T)(a unitaryQ'withA0 Z = Q' S,B0 Z = Q' T,These tables come from local instrumented runs on the targets and revisions named above and below; the temporary
eprintln!and dumps behind them are not in the branch, and their raw logs are not attached to this PR. What a reviewer can rerun from the branch alone is the two tests in the next section: they run on every target and thec32one is red before the fix on x86-64 AVX2.residual 2.9e-7). A temporary trace in the branch showed it runs once (
j = 0,ilast = 11,ilazro) forc32onx86-64 with AVX2/FMA. It does not run for
c64, nor forc32on aarch64 (M5 native) or on x86-64 without AVX(Rosetta 2 exposes no AVX), and those all pass. The layout dependence in the issue was not reproduced: on x86-64 AVX2
the result is bitwise the same for a padded
Mat(column stride 16), a column stride of 12 at 64-byte offsets0/8/32/40, and after another
gevd_cplxcall in the same process.Verification
Before, on fork
maine756f13and ond374ec0(the issue's revision), and after, onae68797. Maximumrelative residual over the columns, computed in f64, with
cargo run --releaseandfeatures = ["std", "linalg"]:c32leftc32rightc64leftc64rightTests:
linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columnsgives the unblocked iteration an exact zeroon the diagonal of
B, atj == iloand below a zero subdiagonal, forc32andc64withn = 2, 3, 6, 12. Itchecks that
QandZare unitary, thatA = Q S Z^HandB = Q T Z^Hhold to10 n eps(LAPACK's QZ testsaccept up to 30 for the ratio divided by
n * ulp), thatSandTare upper triangular, and that exactly onebetais zero. This runs the branch on every target. Before the fix it fails on aarch64:(&Q * &S * Z.adjoint() - &A).norm_max() = 2.509against7.0e-15.faer/tests/gevd_infinite_eigenvalue.rsembeds the issue's bit patterns and checks the left and right residuals atc32andc64against10 n eps, plus a near-infinite eigenvalue. Before the fix it fails on x86-64 AVX2 withc32 left residual 1.962e-1 >= 1.43e-5. On aarch64 it passes before and after, because that path does not reachthe branch.
The regression test's own oracle was vacuous, and is now guarded
residuals()in that test folded each column in withf64::max, which returns the other operand when one isNaN.An all-zero left eigenvector column gives
0 / 0for the left residual, theNaNwas dropped bymax, and thecolumn read as perfectly accurate. The right residual is accumulated the same way, so it had the same hole. Injecting
zeros into the eigenvector outputs left the test green on
d8c7353's version of the file. That is a hole in the testthis PR adds, so it is closed here: per column, before anything is folded in, the test now asserts that the left and
right residuals, both numerators and
|beta| / (|alpha| + |beta|)are finite, and that||u||,||v||and theresidual denominator are non-zero. The message carries the column index and the values.
Measured, mutation
let ul = Mat::<T>::zeros(n, n);inserted right aftergevd_cplx(...).unwrap():On
d8c7353(the reviewedae68797rebased, before these asserts) — green:On this head — red, at
faer/tests/gevd_infinite_eigenvalue.rs:220:The mirror mutation on the right eigenvectors,
let ur = Mat::<T>::zeros(n, n);, is red too:A third hole of the same family, from the next review round: the residual was computed as
numerator / (||v|| * denominator). The product overflows to infinity for a column with a large enough eigenvectornorm, and the quotient is then
0— a value that passesis_finite()and reads as a perfect residual. Splitting itinto two divisions closed that one, and the next review round found two more of the same family, both of which every
guard in the file passed:
alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0: the denominator|beta| ||A|| + |alpha| ||B||overflows toinf, andx / infis0(numerator8.877e7, denominatorinf);alpha = 1e-150, beta = 0, u = v = 1e-200 * e0: the numerator underflows to0before any division runs.Both are wrong by
O(1)once normalized:2.977833357222273e-1left and2.941381743992797e-1right, about 20000xthe
c32tolerance. Dividing afterwards cannot recover information the numerator never had, so the oracle now dividesthe eigenvector by its norm and
(alpha, beta)bymax(|alpha|, |beta|)before forming the residual — the residualis invariant under both scalings, so the value is unchanged wherever it used to survive — requires the two norms, that
scale and the denominator to be finite and positive, and divides once instead of twice, because
||u|| = 1and|alpha|, |beta| <= 1make the numerator at most the denominator.Those two columns are now the test
residual_oracle_rejects_unnormalized_eigenpairs, which hands them to the oracledirectly and requires a residual at or above the tolerance. Red first, with that test written against the old
residuals()body (extracted unchanged intocolumn_residual, so what fails is the oracle's value and not a compileerror; the line numbers below are that working state's, not this head's):
Green on this head, same command:
The normalization is not a relaxation: the issue's pencil reads the same as before. On the M5,
c321.772195055005881e-7 / 1.3466182542302985e-7against1.7721950550479728e-7 / 1.3466182541980509e-7unnormalized, and
c643.838923294839156e-16 / 2.9362584762522423e-16against3.864850603564751e-16 / 2.9009055001556393e-16, with a2.7e-14tolerance; the differences are the reordering ofone division, not a different quantity.
Neither mutation is part of this branch. The test file in the pushed commits has no mutation (
grep -c 'MUTATION' faer/tests/gevd_infinite_eigenvalue.rsis 0, and the same for the temporaryeprintln!inqz_cplx/mod.rs) andpasses:
Note on the command:
cargo test -p t4a-faer gevd_infinite_eigenvalue(the filter form) runs zero tests, becausethe filter matches test names and neither test in that binary is called that. Its real output ends:
--test gevd_infinite_eigenvalueis the form that runs it. The library's gevd tests:The runs quoted for the two counterexamples and for the clean file are on the M5 (
arm64), the aarch64 row of theresiduals table, where the
c32column reads1.772e-7before and after the QZ fix because this host does not takethe branch; the earlier runs quoted in that table were made on x86-64 with
avx2andfmain/proc/cpuinfo, whichdoes.
No performance claim: the same number of rotations, on a column pair instead of a row pair.
CHANGELOG.md(Unreleased) and the
PATCHSET.mdrow are in the same commit. The same lines are in upstream Codebergmain5833f4e, from57e25ac(I read the source and did not run it there). Following the deferred-feedback rule inAGENTS.md, this PR adds noupstream-reports/entry.🤖 Generated with Claude Code