evd: nonsymmetric EVD fixes — ?GEEV input scaling, AED sort termination, ?LALN2 2x2 solves (#112, #144, #145)

docs/lapack-validation.md

@@ -130,11 +130,18 @@ fixtures and with the same residual formula:
   grows like an accumulated rounding error and still catches a residual that grows with n:
   the clamp defect fixed in #121 gave 62.5 at n = 63, which a linear `20 * n / 16` (78.75)
   would have let pass. With it the f32, f64 and c64 groups and every parallel group pass.
-  What remains is c32 only: kinds 6 and 7 (scaled to `sqrt(MAX) / 16` and
-  `sqrt(MIN) * 16`) at n = 63 to 65, ratios 41.7 to 59.4 against 39.7 to 40.3.
-  Single-precision LAPACK `cgeev` (through scipy, OpenBLAS) gives 29 to 40 on the same
-  matrices, itself 39.77 against 39.69 on kind 6 at n = 63 (left vectors). That test stays
-  in `scripts/lapack-known-failures.json` under #112.
+  What remained was c32 only: kinds 6 and 7 (scaled to `sqrt(MAX) / 16` and
+  `sqrt(MIN) * 16`) at n = 63 to 65, ratios 41.7 to 59.4 against 39.7 to 40.3, where
+  single-precision LAPACK `cgeev` (through scipy, OpenBLAS) gives 29 to 40. A stage-by-stage
+  comparison with `cgeev` placed the difference in the complex Schur iteration (backward
+  error 226 to 258 against 67 to 93 on the same matrix scaled to unit size): `?GEEV` scales
+  its input into `[smlnum, bignum]` (`smlnum = sqrt(safmin) / eps`) before the Hessenberg
+  reduction, and faer did not. `evd_imp` now does the same, by a power of two, and scales the
+  eigenvalues back; the c32 group passes and the known-failures entry is gone. A sweep over
+  f32/f64/c32/c64, n = 16 to 257 and scales up to the range limits showed that the missing
+  scaling was far more serious below: on random, Grcar and near-Jordan matrices the ratio grew
+  from about `2^-90` (32-bit), and at `2^-100` (32-bit) and `2^-960` (64-bit) and below it was
+  `1e4` to `2e15` where LAPACK gives 2 to 70. Regression test: `geev_scaling_tests`.
 
 ### Phase 4 — Balancing and back-transformation
 

faer/src/linalg/evd/mod.rs

@@ -901,14 +901,19 @@ fn solve_real_shifted_upper_quasi_triangular_system<T: RealField>(
 				);
 				x[i - 1] -= dot0;
 				x[i] -= dot1;
-				let ref a = &A[(i, i)] - p;
-				let ref b = A[(i - 1, i)];
-				let ref c = A[(i, i - 1)];
-				let ref r0 = x[i - 1];
-				let ref r1 = x[i];
-				let ref inv_det = (a.abs2() - b * c).recip();
-				(x[i - 1], x[i]) =
-					((a * r0 - b * r1) * inv_det, (a * r1 - c * r0) * inv_det);
+				// The same `?LALN2` solve as for a complex shift, with a zero
+				// imaginary part. The inverse through `a^2 - b c` underflowed to a
+				// zero determinant for a block whose entries are below the square
+				// root of the smallest positive number, and returned `inf`.
+				let ((y0, _), (y1, _)) = solve_shifted_2x2_cplx(
+					(&A[(i, i)] - p, zero()),
+					A[(i - 1, i)].copy(),
+					A[(i, i - 1)].copy(),
+					(x[i - 1].copy(), zero()),
+					(x[i].copy(), zero()),
+					&eps_norm,
+				);
+				(x[i - 1], x[i]) = (y0, y1);
 				let a0 = x[i - 1].abs();
 				let a1 = x[i].abs();
 				let m = if a0 > a1 { a0 } else { a1 };
@@ -966,6 +971,70 @@ fn solve_real_shifted_upper_quasi_triangular_system<T: RealField>(
 		f0 * f1
 	}
 }
+/// Solves `[[a, b], [c, a]] [x0, x1]^T = [r0, r1]^T` for complex `a` (as `(re, im)`), real
+/// `b`, `c` and complex right-hand sides, as LAPACK's `?LALN2` does for a complex shift of a
+/// 2 x 2 block: when every entry is below `smin` in `|re| + |im|`, the block is taken as
+/// `smin * I`; otherwise Gaussian elimination with complete pivoting, with the second pivot
+/// raised to `smin` when it is smaller. The caller rescales a large solution.
+fn solve_shifted_2x2_cplx<T: RealField>(
+	a: (T, T),
+	b: T,
+	c: T,
+	r0: (T, T),
+	r1: (T, T),
+	smin: &T,
+) -> ((T, T), (T, T)) {
+	type C<T> = (T, T);
+	fn abs1<T: RealField>(z: &C<T>) -> T {
+		z.0.abs() + z.1.abs()
+	}
+	fn mul<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		(&x.0 * &y.0 - &x.1 * &y.1, &x.0 * &y.1 + &x.1 * &y.0)
+	}
+	fn sub<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		(&x.0 - &y.0, &x.1 - &y.1)
+	}
+	fn div<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		// scaled by the larger component of `y`, so `|y|^2` neither overflows nor underflows
+		let s = y.0.abs().fmax(y.1.abs());
+		let (yr, yi) = (&y.0 / &s, &y.1 / &s);
+		let (xr, xi) = (&x.0 / &s, &x.1 / &s);
+		let d = &yr * &yr + &yi * &yi;
+		((&xr * &yr + &xi * &yi) / &d, (&xi * &yr - &xr * &yi) / &d)
+	}
+	let zero = zero::<T>;
+	let m = [[(a.0.copy(), a.1.copy()), (b, zero())], [(c, zero()), a]];
+	let (mut p, mut q) = (0usize, 0usize);
+	let mut cmax = abs1(&m[0][0]);
+	for (i, j) in [(0, 1), (1, 0)] {
+		let v = abs1(&m[i][j]);
+		if v > cmax {
+			cmax = v;
+			(p, q) = (i, j);
+		}
+	}
+	if cmax < *smin {
+		let inv = smin.recip();
+		return (
+			(&r0.0 * &inv, &r0.1 * &inv),
+			(&r1.0 * &inv, &r1.1 * &inv),
+		);
+	}
+	// pivot row `p`, pivot column `q`; the other row and column are `1 - p`, `1 - q`
+	let (rp, ro) = if p == 0 { (r0, r1) } else { (r1, r0) };
+	let u11 = (m[p][q].0.copy(), m[p][q].1.copy());
+	let u12 = (m[p][1 - q].0.copy(), m[p][1 - q].1.copy());
+	let l21 = div(&m[1 - p][q], &u11);
+	let mut u22 = sub(&m[1 - p][1 - q], &mul(&l21, &u12));
+	if abs1(&u22) < *smin {
+		u22 = (smin.copy(), zero());
+	}
+	let y2 = sub(&ro, &mul(&l21, &rp));
+	let x2 = div(&y2, &u22);
+	let x1 = div(&sub(&rp, &mul(&u12, &x2)), &u11);
+	// unknown `q` is `x1`, the other is `x2`
+	if q == 0 { (x1, x2) } else { (x2, x1) }
+}
 fn solve_cplx_shifted_upper_quasi_triangular_system<T: RealField>(
 	A: MatRef<'_, T>,
 	p: T,
@@ -1043,25 +1112,24 @@ fn solve_cplx_shifted_upper_quasi_triangular_system<T: RealField>(
 				let ref r0_im = x[(i - 1, 1)];
 				let ref r1_re = x[(i, 0)];
 				let ref r1_im = x[(i, 1)];
-				let mut z_re = a_re.abs2() - a_im.abs2() - b * c;
-				let mut z_im = (a_re * a_im).mul_pow2(from_f64::<T>(2.0));
-				let mut z = z_re.hypot(&z_im);
-				if z < eps_norm {
-					z_re = eps_norm.copy();
-					z_im = zero();
-					z = z_re.copy();
-				}
-				let ref z = z;
-				let ref z_re = (z_re / z) / z;
-				let ref z_im = (-z_im / z) / z;
-				let ref x0_re = (a_re * r0_re - a_im * r0_im) - b * r1_re;
-				let ref x0_im = (a_re * r0_im + a_im * r0_re) - b * r1_im;
-				let ref x1_re = (a_re * r1_re - a_im * r1_im) - c * r0_re;
-				let ref x1_im = (a_re * r1_im + a_im * r1_re) - c * r0_im;
-				x[(i - 1, 0)] = x0_re * z_re - x0_im * z_im;
-				x[(i - 1, 1)] = x0_re * z_im + x0_im * z_re;
-				x[(i, 0)] = x1_re * z_re - x1_im * z_im;
-				x[(i, 1)] = x1_re * z_im + x1_im * z_re;
+				// `?LALN2`'s complex 2 x 2 solve (#145): complete pivoting, the
+				// second pivot kept at least `eps_norm`, and `eps_norm * I` for a
+				// block whose entries are all below it. Inverting the block through
+				// its determinant compared that determinant (a product of two
+				// entries) with the entry-sized `eps_norm`, and clamping only the
+				// determinant does not give a perturbed system that is solved.
+				let ((y0_re, y0_im), (y1_re, y1_im)) = solve_shifted_2x2_cplx(
+					(a_re.copy(), a_im.copy()),
+					b.copy(),
+					c.copy(),
+					(r0_re.copy(), r0_im.copy()),
+					(r1_re.copy(), r1_im.copy()),
+					&eps_norm,
+				);
+				x[(i - 1, 0)] = y0_re;
+				x[(i - 1, 1)] = y0_im;
+				x[(i, 0)] = y1_re;
+				x[(i, 1)] = y1_im;
 				let mut m = x[(i - 1, 0)].abs();
 				let v = x[(i - 1, 1)].abs();
 				if v > m {
@@ -1833,6 +1901,47 @@ pub fn eigenvectors_from_schur_cplx<T: ComplexField>(
 		Err(EvdError::NoConvergence)
 	}
 }
+/// `?GEEV`'s scaling of its input: when the largest entry is below
+/// `smlnum = sqrt(min_positive) / eps` or above `bignum = 1 / smlnum`, the power of two that
+/// brings it to just inside that bound, else `None`. LAPACK scales by `smlnum / anrm` (or
+/// `bignum / anrm`); a power of two lands within a factor of two of the same bound and does not
+/// round (an entry far below the largest can still underflow when scaling down, as in LAPACK).
+/// The largest entry is the largest `|re|` or `|im|`, within a factor `sqrt(2)` of LAPACK's
+/// modulus: faer's complex modulus is itself inf for `2^1023 + 0i` and zero for `2^-1074 + 0i`,
+/// which would select a zero factor or none. Without the scaling, the Hessenberg reduction and
+/// the Schur iteration of a matrix whose entries are near the bottom of the range work in
+/// subnormals: on Grcar, random and near-Jordan matrices the GET22 ratio grew from about `2^-90`
+/// (32-bit) and reached `1 / eps` by `2^-100` (32-bit) and `2^-960` (64-bit); near the top, the
+/// complex Schur iteration lost a factor of three against `cgeev` (#112).
+fn geev_scaling<T: ComplexField>(A: MatRef<'_, T>) -> Option<T::Real> {
+	let n = A.nrows();
+	let mut anrm = zero::<T::Real>();
+	for j in 0..n {
+		for i in 0..n {
+			anrm = anrm.fmax(A[(i, j)].real().abs()).fmax(A[(i, j)].imag().abs());
+		}
+	}
+	let smlnum = min_positive::<T::Real>().sqrt() / eps::<T::Real>();
+	let bignum = smlnum.recip();
+	let two = from_f64::<T::Real>(2.0);
+	let half = from_f64::<T::Real>(0.5);
+	if anrm == zero::<T::Real>() || (anrm >= smlnum && anrm <= bignum) {
+		return None;
+	}
+	let mut f = one::<T::Real>();
+	if anrm < smlnum {
+		// `anrm * f` in `[smlnum, 2 smlnum)`
+		while &anrm * &f < smlnum {
+			f = f * &two;
+		}
+	} else {
+		// `anrm * f` in `(bignum / 2, bignum]`
+		while &anrm * &f > bignum {
+			f = f * &half;
+		}
+	}
+	Some(f)
+}
 fn evd_imp<T: ComplexField>(
 	A: MatRef<'_, T>,
 	s: ColMut<'_, T>,
@@ -1879,6 +1988,11 @@ fn evd_imp<T: ComplexField>(
 		None
 	};
 	H.copy_from(A);
+	// the eigenvectors of `f A` are those of `A`; its eigenvalues are scaled back at the end
+	let scaling = geev_scaling(A);
+	if let Some(f) = &scaling {
+		z!(H.rb_mut()).for_each(|uz!(h)| *h = h.mul_real(f));
+	}
 	{
 		let (mut householder, stack) =
 			unsafe { temp_mat_uninit::<T, _, _>(bs, n - 1, stack) };
@@ -2013,6 +2127,13 @@ fn evd_imp<T: ComplexField>(
 			return Err(EvdError::NoConvergence);
 		}
 	}
+	if let Some(f) = &scaling {
+		let f = f.recip();
+		z!(s.rb_mut()).for_each(|uz!(x)| *x = x.mul_real(&f));
+		if let Some(s_im) = s_im.rb_mut() {
+			z!(s_im).for_each(|uz!(x)| *x = x.mul_real(&f));
+		}
+	}
 	Ok(())
 }
 /// computes the matrix $A$'s eigendecomposition
@@ -3101,6 +3222,134 @@ mod tridiagonal_range_tests {
 	}
 }
 
+#[cfg(test)]
+mod geev_scaling_tests {
+	//! Regression tests for the input scaling of `evd_imp` (#112). A matrix whose largest entry is
+	//! near the bottom of the range went through the Hessenberg reduction and the Schur iteration
+	//! in subnormals: a Grcar or random matrix scaled by `2^-120` (`f32`, `c32`) or `2^-1010`
+	//! (`f64`, `c64`) had a GET22 ratio of about `1 / eps` (`4e6` and `2e15`). Near the top the
+	//! complex Schur iteration lost a factor of three against `cgeev`. Residuals are measured in
+	//! `f64` on the dimensionless matrix `a / scale`, exact for a power of two. The left side is
+	//! checked as the right eigenvectors of `a^H`; the `u_left` branch of `evd_imp` itself is
+	//! exercised by the LAPACK nonsymmetric EVD tests.
+	use crate::{assert, c32, c64, Mat, Par};
+	use num_complex::Complex64 as Z;
+
+	/// GET22's `||a u - u diag(w)||_1 / (||a||_1 ||u||_1 eps)` for one side
+	fn get22(a: &[Vec<Z>], u: &[Vec<Z>], w: &[Z], eps: f64) -> f64 {
+		let n = a.len();
+		let (mut an, mut un, mut err) = (0.0f64, 0.0f64, 0.0f64);
+		for j in 0..n {
+			let (mut ac, mut uc, mut ec) = (0.0, 0.0, 0.0);
+			for i in 0..n {
+				ac += a[i][j].norm();
+				uc += u[i][j].norm();
+				let mut r = -u[i][j] * w[j];
+				for k in 0..n {
+					r += a[i][k] * u[k][j];
+				}
+				ec += r.norm();
+			}
+			an = an.max(ac);
+			un = un.max(uc);
+			err = err.max(ec);
+		}
+		err / an / un / eps
+	}
+
+	/// the Grcar matrix (ones on the diagonal and the three superdiagonals, minus one on the
+	/// subdiagonal) and a pseudo-random matrix, both with largest entry one
+	fn unit(n: usize, cplx: bool) -> [Vec<Vec<Z>>; 2] {
+		let grcar = (0..n)
+			.map(|i| {
+				(0..n)
+					.map(|j| match j as isize - i as isize {
+						-1 => Z::new(-1.0, 0.0),
+						0..=3 => Z::new(1.0, 0.0),
+						_ => Z::new(0.0, 0.0),
+					})
+					.collect()
+			})
+			.collect();
+		let mut state = 0x2545_f491_4f6c_dd1du64;
+		let mut next = move || {
+			state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
+			// a multiple of `2^-10` in `[-1, 1)`, exact in every precision
+			((state >> 54) as f64 - 512.0) / 512.0
+		};
+		let rand = (0..n)
+			.map(|_| (0..n).map(|_| Z::new(next(), if cplx { next() } else { 0.0 })).collect())
+			.collect();
+		[grcar, rand]
+	}
+
+	macro_rules! check {
+		($t:ty, $mk:expr, $z:expr, $eps:expr, $scales:expr, $cplx:expr) => {{
+			let n = 32;
+			// the GET22 threshold of the LAPACK nonsymmetric EVD tests
+			let thresh = 20.0 * f64::max(1.0, (n as f64 / 16.0).sqrt());
+			for unit in unit(n, $cplx) {
+				for e in $scales {
+					let scale = 2.0f64.powi(e);
+					let a = Mat::<$t>::from_fn(n, n, |i, j| $mk(unit[i][j] * scale));
+					for (side, m) in [("right", a.clone()), ("left", a.adjoint().to_owned())] {
+						let dimless: Vec<Vec<Z>> = (0..n).map(|i| (0..n).map(|j| $z(m[(i, j)]) / scale).collect()).collect();
+						let evd = m.eigen(Par::Seq).unwrap();
+						let w: Vec<Z> = (0..n).map(|i| Z::new(evd.S()[i].re as f64, evd.S()[i].im as f64) / scale).collect();
+						let u: Vec<Vec<Z>> = (0..n).map(|i| (0..n).map(|j| Z::new(evd.U()[(i, j)].re as f64, evd.U()[(i, j)].im as f64)).collect()).collect();
+						// `f64::max` drops a NaN, so a nonfinite result is rejected first
+						assert!(w.iter().all(|w| w.is_finite()) && u.iter().flatten().all(|u| u.is_finite()));
+						let r = get22(&dimless, &u, &w, $eps);
+						assert!(r <= thresh, "{} scale 2^{e} {side}: GET22 {r}", stringify!($t));
+					}
+				}
+			}
+		}};
+	}
+
+	/// The ends of the range for a complex matrix, where faer's complex modulus overflows
+	/// (`2^1023 + 0i`, `2^127 + 0i`) or underflows (`2^-1074`, `2^-149`): the eigenvalues come
+	/// back exactly.
+	#[test]
+	fn evd_scaling_at_the_complex_range_ends() {
+		let top = c64::new(2.0f64.powi(1023), 0.0);
+		let evd = Mat::from_fn(1, 1, |_, _| top).eigen(Par::Seq).unwrap();
+		assert!(evd.S()[0] == top);
+		let top = c32::new(2.0f32.powi(127), 0.0);
+		let evd = Mat::from_fn(1, 1, |_, _| top).eigen(Par::Seq).unwrap();
+		assert!(evd.S()[0] == top);
+		// a rotation by the smallest subnormal: eigenvalues `+-i m`
+		let m = f64::from_bits(1);
+		let a = Mat::from_fn(2, 2, |i, j| match (i, j) {
+			(0, 1) => c64::new(-m, 0.0),
+			(1, 0) => c64::new(m, 0.0),
+			_ => c64::new(0.0, 0.0),
+		});
+		let evd = a.eigen(Par::Seq).unwrap();
+		let mut w = [evd.S()[0], evd.S()[1]];
+		w.sort_by(|a, b| a.im.partial_cmp(&b.im).unwrap());
+		assert!(w == [c64::new(0.0, -m), c64::new(0.0, m)]);
+		let m = f32::from_bits(1);
+		let a = Mat::from_fn(2, 2, |i, j| match (i, j) {
+			(0, 1) => c32::new(-m, 0.0),
+			(1, 0) => c32::new(m, 0.0),
+			_ => c32::new(0.0, 0.0),
+		});
+		let evd = a.eigen(Par::Seq).unwrap();
+		let mut w = [evd.S()[0], evd.S()[1]];
+		w.sort_by(|a, b| a.im.partial_cmp(&b.im).unwrap());
+		assert!(w == [c32::new(0.0, -m), c32::new(0.0, m)]);
+	}
+
+	#[test]
+	fn evd_scales_its_input_into_range() {
+		check!(f32, |z: Z| z.re as f32, |x: f32| Z::new(x as f64, 0.0), f32::EPSILON as f64, [-120, -100, 0, 100, 120], false);
+		check!(c32, |z: Z| c32::new(z.re as f32, z.im as f32), |x: c32| Z::new(x.re as f64, x.im as f64), f32::EPSILON as f64, [-120, -100, 0, 100, 120], true);
+		check!(f64, |z: Z| z.re, |x: f64| Z::new(x, 0.0), f64::EPSILON, [-1010, -1000, -960, 0, 960, 1015], false);
+		check!(c64, |z: Z| c64::new(z.re, z.im), |x: c64| Z::new(x.re, x.im), f64::EPSILON, [-1010, -1000, -960, 0, 960, 1015], true);
+	}
+}
+
 #[cfg(test)]
 #[path = "tests/mod.rs"]
 mod lapack;