fix(gevd): rotate the columns of Q when the complex QZ chases a zero down B (#324)

CHANGELOG.md

@@ -22,6 +22,12 @@ Changes since `t4a-faer` 0.32.0 (tag `faer-v0.32.0`). Numbers are tensor4all/fae
 
 ## Correctness
 
+- `gevd_cplx` returned wrong left eigenvectors for some pencils with an infinite eigenvalue (a zero on the diagonal of `B` during the QZ
+  iteration's unblocked path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path
+  also ends in), while the eigenvalues and the right eigenvectors were right. The rotations that chase that zero to the bottom were
+  applied to the rows of `Q` instead of its columns. Whether a pencil reaches that branch depends on rounding, so the same input could fail on
+  x86-64 with AVX2 and pass on aarch64 (#324).
+
 - `stats::col_mean` and `stats::col_varm` wrote their results in reverse order when the input had a negative row stride (for example
   `mat.reverse_rows()`), and `row_mean`/`row_varm` did the same for a negative column stride (`mat.reverse_cols()`). Both `NanHandling` modes
   were affected, and `varm` also subtracted the mean of the mirrored row. The input rows were reversed to get a positive stride, but `out` and the

PATCHSET.md

@@ -194,6 +194,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. |
 | Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::<I>(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B, which is the supernodal symbolic workspace rather than the COLAMD one. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus the process-wide allocation counters in `no_step_alloc` and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which pass standalone and flake only inside the whole-suite parallel run. |
 | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex<T>` SIMD path, and the scalar fallbacks of `ComplexImpl<f32>` / `ComplexImpl<f64>` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex<f64>` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. |
+| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it normalizes the eigenvector by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)`, which leaves the residual unchanged, and requires the two eigenvector norms, that scale and the residual denominator to be finite and positive; a numerator formed from the values as they come loses the residual instead, and no later division brings it back: `u = v = 1e-200 * e0` with `alpha = 1e-150` underflows it to `0`, and `alpha = f64::MAX / 2` with `beta = 0` sends the denominator to `inf`, whose quotient is `0` as well. It then asserts both residuals finite, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::<T>::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it; and `residual_oracle_rejects_unnormalized_eigenpairs` hands the oracle the review's two columns directly (`alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0` and `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`): on the oracle without normalization both read `left residual 0e0 < 1.430511474609375e-5` and `right residual 0e0 < 1.430511474609375e-5`, and with it `2.977833357222273e-1` and `2.941381743992797e-1`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). |
 
 ## Ground rules
 

faer/src/linalg/gevd/qz_cplx/mod.rs

@@ -120,11 +120,18 @@ fn hessenberg_to_qz_unblocked<T: ComplexField>(
 											.get_mut(.., jch + 1..ilastm + 1)
 											.two_rows_mut(jch, jch + 1);
 										rot(c.copy(), s.copy(), x, y);
-										if let Some(mut Q) = Q.rb_mut() {
-											let (x, y) = Q
-												.rb_mut()
-												.two_rows_mut(jch, jch + 1);
-											rot(c.copy(), s.copy(), x, y);
+										// the rotation acts on rows `jch, jch + 1` of `H`
+										// and `T`, so `Q` takes its adjoint on the right:
+										// columns `jch, jch + 1`, with `conj(s)` (faer#324)
+										if let Some(Q) = Q.rb_mut() {
+											let (x, y) =
+												Q.two_cols_mut(jch, jch + 1);
+											rot(
+												c.copy(),
+												s.conj(),
+												x.transpose_mut(),
+												y.transpose_mut(),
+											);
 										}
 										if ilazr2 {
 											H[(jch, jch - 1)] =
@@ -1574,6 +1581,91 @@ mod tests {
 			}
 		}
 	}
+	/// faer#324: a zero on the diagonal of `B` at `j`, with `j == ilo` or a
+	/// negligible `A[(j, j - 1)]`, is chased to the bottom by rotations applied
+	/// from the left. those rotations must update the columns of `Q`; they used
+	/// to rotate its rows, so `Q` stayed unitary but no longer satisfied
+	/// `A = Q S Z^H`, and every left eigenvector built from it was wrong.
+	fn check_unblocked_zero_chase<T: ComplexField>(
+		A: Mat<T>,
+		B: Mat<T>,
+		zero_at: usize,
+	) {
+		let n = A.nrows();
+		let mut A = A;
+		let mut B = B;
+		for j in 0..n {
+			for i in j + 2..n {
+				A[(i, j)] = zero();
+			}
+			for i in j + 1..n {
+				B[(i, j)] = zero();
+			}
+		}
+		B[(zero_at, zero_at)] = zero();
+		if zero_at > 0 {
+			A[(zero_at, zero_at - 1)] = zero();
+		}
+		let mut Q = Mat::<T>::identity(n, n);
+		let mut Z = Mat::<T>::identity(n, n);
+		let mut S = A.clone();
+		let mut P = B.clone();
+		let mut alpha = Col::<T>::zeros(n);
+		let mut beta = Col::<T>::zeros(n);
+		hessenberg_to_qz_unblocked(
+			0,
+			n - 1,
+			S.as_mut(),
+			P.as_mut(),
+			Some(Q.as_mut()),
+			Some(Z.as_mut()),
+			alpha.as_mut(),
+			beta.as_mut(),
+			false,
+		);
+		// LAPACK's QZ tests divide by `n * ulp` and accept ratios up to 30
+		let ref tol = from_f64::<T::Real>(10.0 * n as f64) * eps::<T::Real>();
+		let I = Mat::<T>::identity(n, n);
+		assert!((Q.adjoint() * &Q - &I).norm_max() < *tol);
+		assert!((Z.adjoint() * &Z - &I).norm_max() < *tol);
+		assert!((&Q * &S * Z.adjoint() - &A).norm_max() < tol * A.norm_max());
+		assert!((&Q * &P * Z.adjoint() - &B).norm_max() < tol * B.norm_max());
+		for j in 0..n {
+			for i in j + 1..n {
+				::core::assert!(
+					S[(i, j)] == zero::<T>() && P[(i, j)] == zero::<T>()
+				);
+			}
+		}
+		::core::assert_eq!(
+			beta.iter().filter(|b| **b == zero::<T>()).count(),
+			1
+		);
+	}
+	#[test]
+	fn unblocked_zero_chase_updates_q_columns() {
+		let rng = &mut StdRng::seed_from_u64(324);
+		for n in [2, 3, 6, 12] {
+			let rand = stats::CwiseMatDistribution {
+				nrows: n,
+				ncols: n,
+				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
+			};
+			let A: Mat<c64> = rand.sample(rng);
+			let B: Mat<c64> = rand.sample(rng);
+			let to_c32 = |M: &Mat<c64>| {
+				Mat::from_fn(n, n, |i, j| {
+					c32::new(M[(i, j)].re as f32, M[(i, j)].im as f32)
+				})
+			};
+			// `zero_at == 0` takes the `j == ilo` path, `zero_at == n / 2` the
+			// negligible subdiagonal path
+			for zero_at in [0, n / 2] {
+				check_unblocked_zero_chase(A.clone(), B.clone(), zero_at);
+				check_unblocked_zero_chase(to_c32(&A), to_c32(&B), zero_at);
+			}
+		}
+	}
 	#[test]
 	fn test_qz_cplx_blocked() {
 		let rng = &mut StdRng::seed_from_u64(0);