From 787bf99efd18c050b5e9f172d8bcb38f17f8830e Mon Sep 17 00:00:00 2001 From: Eli Date: Thu, 16 Jul 2026 00:43:59 -0300 Subject: [PATCH 1/4] Implement GMRES(m) restarted solver --- examples/krylov/gmres.rs | 57 +++ examples/krylov/main.rs | 11 +- src/krylov/gmres.rs | 642 +++++++++++++++++++++++++++++++++ src/krylov/mod.rs | 4 + tests/property/krylov/gmres.rs | 93 +++++ tests/property/krylov/mod.rs | 1 + 6 files changed, 805 insertions(+), 3 deletions(-) create mode 100644 examples/krylov/gmres.rs create mode 100644 src/krylov/gmres.rs create mode 100644 tests/property/krylov/gmres.rs diff --git a/examples/krylov/gmres.rs b/examples/krylov/gmres.rs new file mode 100644 index 0000000..5c72ecb --- /dev/null +++ b/examples/krylov/gmres.rs @@ -0,0 +1,57 @@ +use rustebra::krylov::gmres; +use rustebra::sparse::CsrMatrix; +use rustebra::storage::Basis; + +pub(crate) fn run() { + println!("\n== GMRES(m) =="); + + // Non-symmetric, non-SPD system, unlike Conjugate Gradient's requirements: + // [[4, 1], [2, 3]] x = [1, 2]. Solution: x = [0.1, 0.6]. + let a = CsrMatrix::new( + 2, + 2, + vec![0, 2, 4], + vec![0, 1, 0, 1], + vec![4.0_f64, 1.0, 2.0, 3.0], + ) + .expect("valid CSR"); + let b = [1.0, 2.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 4]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).expect("converges"); + println!("full-basis solve: x = {out_x:?}"); + + // Restart size smaller than the problem dimension (M = 1 < n = 3): the solution still + // emerges, just carried forward across several restart cycles instead of one. + let a3 = CsrMatrix::new( + 3, + 3, + vec![0, 2, 5, 7], + vec![0, 1, 0, 1, 2, 1, 2], + vec![5.0_f64, 1.0, 1.0, 4.0, 1.0, 1.0, 3.0], + ) + .expect("valid CSR"); + let b3 = [6.0, 6.0, 4.0]; + let x0_3 = [0.0, 0.0, 0.0]; + let mut out_x3 = [0.0; 3]; + let mut buffer3 = [0.0; 3]; + let mut basis3 = Basis::::new(&mut buffer3, 3).unwrap(); + let mut scratch3 = [0.0; 3]; + + gmres( + &a3, + &b3, + &x0_3, + 500, + 1e-10, + &mut out_x3, + &mut basis3, + &mut scratch3, + ) + .expect("converges across restarts"); + println!("restarted GMRES(1) solve: x = {out_x3:?}"); +} diff --git a/examples/krylov/main.rs b/examples/krylov/main.rs index 02fb8bf..80b3812 100644 --- a/examples/krylov/main.rs +++ b/examples/krylov/main.rs @@ -1,14 +1,19 @@ //! Tour `rustebra::krylov`'s subspace methods: Lanczos iteration, which builds an orthonormal //! basis of a Krylov subspace and the symmetric tridiagonal matrix a symmetric operator -//! projects onto within it, and Arnoldi iteration, its non-symmetric counterpart producing an -//! upper Hessenberg projection. +//! projects onto within it, Arnoldi iteration, its non-symmetric counterpart producing an +//! upper Hessenberg projection, and GMRES(m), the restarted linear solver built on top of it. //! -//! Run with: `cargo run --example krylov` +//! Run with: `cargo run --example krylov` (add `--features alloc` for the GMRES section, since +//! it takes its operator as a `SparseLinearOp`). mod arnoldi; +#[cfg(feature = "alloc")] +mod gmres; mod lanczos; fn main() { lanczos::run(); arnoldi::run(); + #[cfg(feature = "alloc")] + gmres::run(); } diff --git a/src/krylov/gmres.rs b/src/krylov/gmres.rs new file mode 100644 index 0000000..7cefb04 --- /dev/null +++ b/src/krylov/gmres.rs @@ -0,0 +1,642 @@ +use super::ConvergenceError; +use super::power_iteration::{Slice, normalize}; +use crate::algorithm::vector::{dot, norm}; +use crate::scalar::Scalar; +use crate::sparse::SparseLinearOp; +use crate::storage::Basis; + +/// `y -= coefficient * x`, the orthogonalization step's only vector update. +fn subtract_scaled(y: &mut [T], coefficient: T, x: &[T]) { + for (slot, &x_i) in y.iter_mut().zip(x.iter()) { + *slot = slot.sub(coefficient.mul(x_i)); + } +} + +/// Computes `r = b - a * x` into `scratch` and returns `‖r‖`, or `NonFinite` if that norm is +/// neither positive nor zero (the only way a norm can fail both comparisons). +fn residual_norm( + a: &impl SparseLinearOp, + b: &[T], + x: &[T], + scratch: &mut [T], +) -> Result { + a.apply(x, scratch) + .map_err(|_| ConvergenceError::DimensionMismatch)?; + for (slot, &b_i) in scratch.iter_mut().zip(b.iter()) { + *slot = b_i.sub(*slot); + } + let r_norm = norm(&Slice { data: &*scratch }); + if r_norm > T::zero() || r_norm == T::zero() { + Ok(r_norm) + } else { + Err(ConvergenceError::NonFinite) + } +} + +/// Solves the general (possibly non-symmetric) linear system `A x = b` via restarted GMRES, +/// GMRES(`M`): Arnoldi iteration builds an `M`-dimensional Krylov basis from the current +/// residual, the resulting least-squares problem over that basis is solved via Givens +/// rotations, and the cycle restarts from the improved iterate until either the residual +/// meets `tol` or `max_restarts` cycles are exhausted. +/// +/// `A` is supplied as a [`SparseLinearOp`] rather than a dense matrix: applying it never +/// allocates, so restart cycles reuse the same workspace (`out_x`, `basis`, `scratch`) +/// throughout. +/// +/// # Algorithm +/// +/// Each restart cycle: +/// +/// 1. Computes the residual `r = b - A x` and its norm `β = ‖r‖`, returning `Ok` immediately +/// if `β <= tol`. +/// 2. Runs Arnoldi iteration from `q_0 = r / β`, building an orthonormal basis `Q` of up to +/// `M` vectors and the upper Hessenberg projection `H`, stopping early (before `M` steps) +/// on breakdown — an invariant subspace found before the basis filled up, the same +/// "success, not failure" case documented on [`super::arnoldi`]. +/// 3. Solves `min_y ‖β e_1 - H y‖` (a `(reached + 1) x reached` least-squares problem) via +/// incremental Givens rotations, then updates `x <- x + Q y`. +/// +/// # Convergence +/// +/// GMRES's residual norm decreases monotonically within a cycle (each additional basis +/// vector can only improve the least-squares fit) and never increases across a restart, +/// because restarting recomputes the same residual the next cycle continues from. It is not +/// guaranteed to decrease *strictly* every cycle, though: a starting vector aligned with an +/// invariant subspace the operator doesn't expand (breakdown on the very first Arnoldi step) +/// leaves `x` unchanged, and the iteration stagnates. Restarting also discards the larger +/// Krylov subspace full (non-restarted) GMRES would have kept building, so GMRES(`M`) can +/// converge slower, or stagnate on problems full GMRES would resolve — the restart budget +/// `max_restarts` bounds the cost of that risk rather than eliminating it. +/// +/// # Errors +/// +/// - [`ConvergenceError::DimensionMismatch`] if `a.rows() != a.cols()`, `b`, `x0`, `out_x`, +/// or a `basis` vector doesn't have exactly `a.rows()` elements, or `M > a.rows()`. +/// - [`ConvergenceError::NonFinite`] if a residual or Arnoldi iterate goes `NaN` or infinite. +/// - [`ConvergenceError::Breakdown`] if the small least-squares system built from `H` has a +/// zero pivot that Arnoldi's own breakdown test didn't already catch (a coincidental exact +/// singularity in the projected system). +/// - [`ConvergenceError::MaxIterationsExceeded`] if the residual hasn't met `tol` after +/// `max_restarts` cycles. +/// +/// # Examples +/// +/// ``` +/// use rustebra::krylov::gmres; +/// use rustebra::sparse::CsrMatrix; +/// use rustebra::storage::Basis; +/// +/// // [[4, 1], [2, 3]] x = [1, 2]. Solution: x = [0.1, 0.6]. +/// let a = CsrMatrix::new(2, 2, vec![0, 2, 4], vec![0, 1, 0, 1], vec![4.0_f64, 1.0, 2.0, 3.0]) +/// .unwrap(); +/// let b = [1.0, 2.0]; +/// let x0 = [0.0, 0.0]; +/// let mut out_x = [0.0; 2]; +/// let mut buffer = [0.0; 4]; +/// let mut basis = Basis::::new(&mut buffer, 2).unwrap(); +/// let mut scratch = [0.0; 2]; +/// +/// gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).unwrap(); +/// +/// assert!((out_x[0] - 0.1).abs() < 1e-8); +/// assert!((out_x[1] - 0.6).abs() < 1e-8); +/// ``` +#[allow(clippy::too_many_arguments)] +pub fn gmres( + a: &impl SparseLinearOp, + b: &[T], + x0: &[T], + max_restarts: usize, + tol: T, + out_x: &mut [T], + basis: &mut Basis<'_, T, M>, + scratch: &mut [T], +) -> Result<(), ConvergenceError> +where + T: Scalar + PartialOrd, +{ + let n = a.rows(); + if a.cols() != n + || b.len() != n + || x0.len() != n + || out_x.len() != n + || basis.vector_len() != n + || scratch.len() != n + || M > n + { + return Err(ConvergenceError::DimensionMismatch); + } + + out_x.copy_from_slice(x0); + + let mut beta = residual_norm(a, b, out_x, scratch)?; + if beta <= tol { + return Ok(()); + } + + for _ in 0..max_restarts { + if M == 0 { + // No basis vectors can be built; the residual can't be reduced this cycle. + beta = residual_norm(a, b, out_x, scratch)?; + if beta <= tol { + return Ok(()); + } + continue; + } + + // q_0 = r / beta, where `r` (unnormalized) is already sitting in `scratch` from the + // residual computation above. + normalize(scratch)?; + { + let q_0 = basis + .vector_mut(0) + .ok_or(ConvergenceError::DimensionMismatch)?; + q_0.copy_from_slice(scratch); + } + + let mut h = [[T::zero(); M]; M]; + let mut h_sub = [T::zero(); M]; + let mut reached = 0usize; + + for j in 0..M { + let norm_aq = { + let q_j = basis.vector(j).ok_or(ConvergenceError::DimensionMismatch)?; + a.apply(q_j, scratch) + .map_err(|_| ConvergenceError::DimensionMismatch)?; + norm(&Slice { data: &*scratch }) + }; + + for (i, row) in h.iter_mut().enumerate().take(j + 1) { + let q_i = basis.vector(i).ok_or(ConvergenceError::DimensionMismatch)?; + let h_ij = dot(&Slice { data: q_i }, &Slice { data: &*scratch }) + .map_err(|_| ConvergenceError::DimensionMismatch)?; + row[j] = h_ij; + subtract_scaled(scratch, h_ij, q_i); + } + + let h_next = norm(&Slice { data: &*scratch }); + // `x - x` is `0` for every finite `x` and `NaN` for `NaN`/±infinity — the only + // values unequal to themselves. + let probe = h_next.sub(h_next); + #[allow(clippy::eq_op)] + let non_finite = probe != probe; + if non_finite { + return Err(ConvergenceError::NonFinite); + } + + reached = j + 1; + // Breakdown: the Krylov subspace built so far is already invariant. Stop + // extending the basis; the least-squares solve below uses what was built. + if h_next <= tol.mul(norm_aq) { + break; + } + h_sub[j] = h_next; + + if j + 1 < M { + let q_next = basis + .vector_mut(j + 1) + .ok_or(ConvergenceError::DimensionMismatch)?; + let inv = T::one().div(h_next); + for (slot, &w_i) in q_next.iter_mut().zip(scratch.iter()) { + *slot = w_i.mul(inv); + } + } + } + + // Reduce the (reached + 1) x reached Hessenberg block `[h; h_sub]` to upper + // triangular form via incremental Givens rotations, tracking the same rotations' + // effect on the right-hand side `g` (initialized to `beta * e_1`). + let mut g = [T::zero(); M]; + g[0] = beta; + let mut cs = [T::zero(); M]; + let mut sn = [T::zero(); M]; + + for j in 0..reached { + let mut col = [T::zero(); M]; + for (i, slot) in col.iter_mut().enumerate().take(j + 1) { + *slot = h[i][j]; + } + let sub = h_sub[j]; + + for i in 0..j { + let old_i = col[i]; + let old_i1 = col[i + 1]; + col[i] = cs[i].mul(old_i).add(sn[i].mul(old_i1)); + col[i + 1] = cs[i].mul(old_i1).sub(sn[i].mul(old_i)); + } + + let r = Scalar::sqrt(col[j].mul(col[j]).add(sub.mul(sub))); + let (c, s) = if r == T::zero() { + (T::one(), T::zero()) + } else { + (col[j].div(r), sub.div(r)) + }; + cs[j] = c; + sn[j] = s; + col[j] = r; + + let g_j = g[j]; + let g_next = if j + 1 < M { g[j + 1] } else { T::zero() }; + g[j] = c.mul(g_j).add(s.mul(g_next)); + if j + 1 < M { + g[j + 1] = c.mul(g_next).sub(s.mul(g_j)); + } + + for (i, &value) in col.iter().enumerate().take(j + 1) { + h[i][j] = value; + } + } + + // Back substitution: R y = g, with R the (now upper triangular) leading block of `h`. + let mut y = [T::zero(); M]; + for i in (0..reached).rev() { + let mut sum = g[i]; + for (k, &y_k) in y.iter().enumerate().take(reached).skip(i + 1) { + sum = sum.sub(h[i][k].mul(y_k)); + } + let diag = h[i][i]; + if diag == T::zero() { + return Err(ConvergenceError::Breakdown); + } + y[i] = sum.div(diag); + } + + // x <- x + Q y + for (k, &y_k) in y.iter().enumerate().take(reached) { + let q_k = basis.vector(k).ok_or(ConvergenceError::DimensionMismatch)?; + for (slot, &q_ki) in out_x.iter_mut().zip(q_k.iter()) { + *slot = slot.add(y_k.mul(q_ki)); + } + } + + beta = residual_norm(a, b, out_x, scratch)?; + if beta <= tol { + return Ok(()); + } + } + + Err(ConvergenceError::MaxIterationsExceeded) +} + +#[cfg(test)] +mod tests { + use super::gmres; + use crate::krylov::ConvergenceError; + use crate::sparse::CsrMatrix; + use crate::storage::Basis; + + fn assert_close(actual: f64, expected: f64, tol: f64) { + assert!( + (actual - expected).abs() < tol, + "expected {expected}, got {actual}" + ); + } + + fn csr_from_dense(a: &[f64], n: usize) -> CsrMatrix { + let mut row_ptr = vec![0_u32]; + let mut col_indices = vec![]; + let mut values = vec![]; + for r in 0..n { + for c in 0..n { + let v = a[r * n + c]; + if v != 0.0 { + col_indices.push(c as u32); + values.push(v); + } + } + row_ptr.push(col_indices.len() as u32); + } + CsrMatrix::new(n, n, row_ptr, col_indices, values).unwrap() + } + + fn residual_norm(a: &[f64], n: usize, x: &[f64], b: &[f64]) -> f64 { + let mut r_sq = 0.0; + for row in 0..n { + let mut ax = 0.0; + for col in 0..n { + ax += a[row * n + col] * x[col]; + } + let r = b[row] - ax; + r_sq += r * r; + } + r_sq.sqrt() + } + + #[test] + fn solves_a_small_nonsymmetric_system() { + // [[4, 1], [2, 3]] x = [1, 2]. Solution: x = [0.1, 0.6]. + let a = [4.0, 1.0, 2.0, 3.0]; + let m = csr_from_dense(&a, 2); + let b = [1.0, 2.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 4]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + gmres(&m, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).unwrap(); + + assert_close(out_x[0], 0.1, 1e-8); + assert_close(out_x[1], 0.6, 1e-8); + } + + #[test] + fn solves_a_diagonal_system_in_one_restart() { + // diag(2, 4, 5) x = [4, 8, 10] => x = [2, 2, 2]. Full-basis GMRES solves an + // n-dimensional system exactly within one restart (n steps of Arnoldi). + let a = [2.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 5.0]; + let m = csr_from_dense(&a, 3); + let b = [4.0, 8.0, 10.0]; + let x0 = [0.0, 0.0, 0.0]; + let mut out_x = [0.0; 3]; + let mut buffer = [0.0; 9]; + let mut basis = Basis::::new(&mut buffer, 3).unwrap(); + let mut scratch = [0.0; 3]; + + gmres(&m, &b, &x0, 1, 1e-10, &mut out_x, &mut basis, &mut scratch).unwrap(); + + assert_close(out_x[0], 2.0, 1e-8); + assert_close(out_x[1], 2.0, 1e-8); + assert_close(out_x[2], 2.0, 1e-8); + } + + #[test] + fn restarts_accumulate_progress_toward_the_solution() { + // A well-conditioned 3x3 system, solved with a restart size (M = 1) too small to + // reach the solution in a single cycle: correctness relies on restarts carrying the + // residual forward, not on any single cycle's Krylov subspace being big enough. + let a = [5.0, 1.0, 0.0, 1.0, 4.0, 1.0, 0.0, 1.0, 3.0]; + let m = csr_from_dense(&a, 3); + let b = [6.0, 6.0, 4.0]; + let x0 = [0.0, 0.0, 0.0]; + let mut out_x = [0.0; 3]; + let mut buffer = [0.0; 3]; + let mut basis = Basis::::new(&mut buffer, 3).unwrap(); + let mut scratch = [0.0; 3]; + + gmres( + &m, + &b, + &x0, + 500, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch, + ) + .unwrap(); + + assert!(residual_norm(&a, 3, &out_x, &b) < 1e-8); + } + + #[test] + fn m_equals_one_degenerates_toward_steepest_descent_like_behavior() { + // A symmetric positive-definite system: with M == 1, every restart cycle can only + // move along the current residual direction, the same search direction gradient + // descent would take. It still converges, just gradually, given enough restarts. + let a = [3.0, 1.0, 1.0, 2.0]; + let m = csr_from_dense(&a, 2); + let b = [4.0, 3.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 2]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + gmres( + &m, + &b, + &x0, + 200, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch, + ) + .unwrap(); + + assert!(residual_norm(&a, 2, &out_x, &b) < 1e-8); + } + + #[test] + fn zero_restart_budget_only_accepts_an_already_converged_guess() { + let a = [2.0, 0.0, 0.0, 2.0]; + let m = csr_from_dense(&a, 2); + let b = [2.0, 2.0]; + + // x0 is already the solution: converges without spending any restart. + let x0_exact = [1.0, 1.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 2]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + gmres( + &m, + &b, + &x0_exact, + 0, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch, + ) + .unwrap(); + assert_close(out_x[0], 1.0, 1e-10); + assert_close(out_x[1], 1.0, 1e-10); + + // x0 is not the solution and no restarts are budgeted: no progress can be made. + let x0_wrong = [0.0, 0.0]; + let result = gmres( + &m, + &b, + &x0_wrong, + 0, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch, + ); + assert_eq!(result, Err(ConvergenceError::MaxIterationsExceeded)); + } + + #[test] + fn stagnation_across_restarts_reports_exhaustion_not_a_hang() { + // A 90-degree rotation: for *any* direction q, `A * q` is orthogonal to `q`. With + // M == 1, GMRES minimizes `‖r - y * (A * q_0)‖` over the scalar `y`, and since + // `A * q_0 ⟂ r` (r is a multiple of q_0), that minimum sits exactly at `y == 0` — + // every restart cycle picks the step that changes nothing. A tight budget must fail + // fast with `MaxIterationsExceeded`, not loop forever making zero progress. + let a = [0.0, -1.0, 1.0, 0.0]; + let m = csr_from_dense(&a, 2); + let b = [1.0, 0.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 2]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + let result = gmres(&m, &b, &x0, 20, 1e-10, &mut out_x, &mut basis, &mut scratch); + + assert_eq!(result, Err(ConvergenceError::MaxIterationsExceeded)); + // x never moved off the stagnation point. + assert_close(out_x[0], 0.0, 1e-10); + assert_close(out_x[1], 0.0, 1e-10); + } + + #[test] + fn m_zero_never_reduces_a_nonzero_residual() { + let a = [2.0, 0.0, 0.0, 2.0]; + let m = csr_from_dense(&a, 2); + let b = [2.0, 2.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer: [f64; 0] = []; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + let result = gmres(&m, &b, &x0, 5, 1e-10, &mut out_x, &mut basis, &mut scratch); + + assert_eq!(result, Err(ConvergenceError::MaxIterationsExceeded)); + } + + #[test] + fn already_converged_initial_guess_returns_immediately() { + let a = [2.0, 0.0, 0.0, 2.0]; + let m = csr_from_dense(&a, 2); + let b = [2.0, 4.0]; + let x0 = [1.0, 2.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 2]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + gmres(&m, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).unwrap(); + + assert_close(out_x[0], 1.0, 1e-10); + assert_close(out_x[1], 2.0, 1e-10); + } + + #[test] + fn mismatched_dimensions_are_an_error_not_a_panic() { + let a = csr_from_dense(&[2.0, 0.0, 0.0, 1.0], 2); + let b = [1.0, 0.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 2]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + // `b` too short. + assert_eq!( + gmres( + &a, + &[1.0], + &x0, + 10, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch + ), + Err(ConvergenceError::DimensionMismatch) + ); + + // `x0` too short. + assert_eq!( + gmres( + &a, + &b, + &[1.0], + 10, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch + ), + Err(ConvergenceError::DimensionMismatch) + ); + + // `out_x` too short. + let mut out_x_short = [0.0; 1]; + assert_eq!( + gmres( + &a, + &b, + &x0, + 10, + 1e-10, + &mut out_x_short, + &mut basis, + &mut scratch + ), + Err(ConvergenceError::DimensionMismatch) + ); + + // `scratch` too short. + let mut scratch_short = [0.0; 1]; + assert_eq!( + gmres( + &a, + &b, + &x0, + 10, + 1e-10, + &mut out_x, + &mut basis, + &mut scratch_short + ), + Err(ConvergenceError::DimensionMismatch) + ); + + // `basis` vectors of the wrong length. + let mut buffer_wrong = [0.0; 1]; + let mut basis_wrong = Basis::::new(&mut buffer_wrong, 1).unwrap(); + assert_eq!( + gmres( + &a, + &b, + &x0, + 10, + 1e-10, + &mut out_x, + &mut basis_wrong, + &mut scratch + ), + Err(ConvergenceError::DimensionMismatch) + ); + + // M > n: a 2-dimensional space has no 3 orthonormal directions. + let mut buffer_deep = [0.0; 6]; + let mut basis_deep = Basis::::new(&mut buffer_deep, 2).unwrap(); + assert_eq!( + gmres( + &a, + &b, + &x0, + 10, + 1e-10, + &mut out_x, + &mut basis_deep, + &mut scratch + ), + Err(ConvergenceError::DimensionMismatch) + ); + } + + #[test] + fn non_square_operator_is_a_dimension_mismatch() { + // 2x3, not square: GMRES requires a square operator. + let a = + CsrMatrix::new(2, 3, vec![0, 2, 3], vec![0, 1, 2], vec![1.0_f64, 1.0, 1.0]).unwrap(); + let b = [1.0, 1.0]; + let x0 = [0.0, 0.0, 0.0]; + let mut out_x = [0.0; 3]; + let mut buffer = [0.0; 6]; + let mut basis = Basis::::new(&mut buffer, 3).unwrap(); + let mut scratch = [0.0; 3]; + + let result = gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch); + + assert_eq!(result, Err(ConvergenceError::DimensionMismatch)); + } +} diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index f45d7dc..ae657f8 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -1,5 +1,7 @@ mod arnoldi; mod conjugate_gradient; +#[cfg(feature = "alloc")] +mod gmres; mod hessenberg; mod inverse_power_iteration; mod lanczos; @@ -8,6 +10,8 @@ mod tridiagonal; pub use self::arnoldi::arnoldi; pub use self::conjugate_gradient::conjugate_gradient; +#[cfg(feature = "alloc")] +pub use self::gmres::gmres; pub use self::hessenberg::HessenbergMatrix; pub use self::inverse_power_iteration::inverse_power_iteration; pub use self::lanczos::lanczos; diff --git a/tests/property/krylov/gmres.rs b/tests/property/krylov/gmres.rs new file mode 100644 index 0000000..5458db5 --- /dev/null +++ b/tests/property/krylov/gmres.rs @@ -0,0 +1,93 @@ +#![cfg(feature = "alloc")] + +//! Property test for `gmres`: on diagonally dominant systems (which guarantee a well-posed, +//! well-conditioned solve), the final residual `‖b - A x‖` must fall below the requested +//! tolerance, checked against an independently computed matrix-vector product rather than +//! anything the solver itself produced internally. + +use proptest::prelude::*; +use rustebra::krylov::gmres; +use rustebra::sparse::CsrMatrix; +use rustebra::storage::Basis; + +const N: usize = 4; +const TOL: f64 = 1e-9; + +/// Dense row-major `a` as a `CsrMatrix`, storing every entry (including exact zeros) so the +/// generator's sparsity pattern never has to be tracked separately. +fn csr_from_dense(a: &[f64; N * N]) -> CsrMatrix { + let mut row_ptr = vec![0_u32]; + let mut col_indices = vec![]; + let mut values = vec![]; + for r in 0..N { + for c in 0..N { + col_indices.push(c as u32); + values.push(a[r * N + c]); + } + row_ptr.push(col_indices.len() as u32); + } + CsrMatrix::new(N, N, row_ptr, col_indices, values).unwrap() +} + +fn residual_norm(a: &[f64; N * N], x: &[f64; N], b: &[f64; N]) -> f64 { + let mut sq = 0.0; + for r in 0..N { + let mut ax = 0.0; + for c in 0..N { + ax += a[r * N + c] * x[c]; + } + let ri = b[r] - ax; + sq += ri * ri; + } + sq.sqrt() +} + +prop_compose! { + /// A row-major, strictly diagonally dominant `N x N` matrix: each diagonal entry's + /// magnitude exceeds the sum of the magnitudes of the rest of its row, which guarantees + /// non-singularity and a well-conditioned solve for GMRES to converge on within a small + /// restart budget. + fn diagonally_dominant_matrix()( + off_diagonal in prop::array::uniform16(-1.0..1.0f64), + diagonal_boost in prop::array::uniform4(5.0..10.0f64), + ) -> [f64; N * N] { + let mut a = off_diagonal; + for r in 0..N { + a[r * N + r] = 0.0; + } + let mut row_sums = [0.0; N]; + for r in 0..N { + row_sums[r] = (0..N).map(|c| a[r * N + c].abs()).sum(); + } + for r in 0..N { + a[r * N + r] = row_sums[r] + diagonal_boost[r]; + } + a + } +} + +proptest! { + /// GMRES(N) — a full-dimension restart size — solves a diagonally dominant system to + /// within `TOL` in a small restart budget, verified against an independently computed + /// residual rather than any value the solver reports about itself. + #[test] + fn residual_is_small_after_convergence( + a in diagonally_dominant_matrix(), + b in prop::array::uniform4(-10.0..10.0f64), + ) { + let m = csr_from_dense(&a); + let x0 = [0.0; N]; + let mut out_x = [0.0; N]; + let mut buffer = [0.0; N * N]; + let mut basis = Basis::::new(&mut buffer, N).unwrap(); + let mut scratch = [0.0; N]; + + gmres(&m, &b, &x0, 10, TOL, &mut out_x, &mut basis, &mut scratch).unwrap(); + + let residual = residual_norm(&a, &out_x, &b); + prop_assert!( + residual < 1e-6, + "residual {residual} too large for a converged solve", + ); + } +} diff --git a/tests/property/krylov/mod.rs b/tests/property/krylov/mod.rs index a131948..88509e3 100644 --- a/tests/property/krylov/mod.rs +++ b/tests/property/krylov/mod.rs @@ -3,6 +3,7 @@ mod arnoldi; #[allow(dead_code)] mod common; +mod gmres; mod inverse_power_iteration; mod lanczos; mod power_iteration; From 0023393d20629e88befcae386177c6969ee26878 Mon Sep 17 00:00:00 2001 From: Eli Date: Thu, 16 Jul 2026 01:13:44 -0300 Subject: [PATCH 2/4] - Added gmres property tests - Added krylov edge cases --- src/krylov/gmres.rs | 13 +++-- src/krylov/power_iteration.rs | 22 +++++---- tests/edge_cases/krylov.rs | 60 ++++++++++++++++++++++- tests/property/diff/gmres.rs | 89 +++++++++++++++++++++++++++++++++++ tests/property/diff/mod.rs | 1 + 5 files changed, 171 insertions(+), 14 deletions(-) create mode 100644 tests/property/diff/gmres.rs diff --git a/src/krylov/gmres.rs b/src/krylov/gmres.rs index 7cefb04..622e4b2 100644 --- a/src/krylov/gmres.rs +++ b/src/krylov/gmres.rs @@ -13,7 +13,7 @@ fn subtract_scaled(y: &mut [T], coefficient: T, x: &[T]) { } /// Computes `r = b - a * x` into `scratch` and returns `‖r‖`, or `NonFinite` if that norm is -/// neither positive nor zero (the only way a norm can fail both comparisons). +/// `NaN` or infinite. fn residual_norm( a: &impl SparseLinearOp, b: &[T], @@ -26,10 +26,15 @@ fn residual_norm( *slot = b_i.sub(*slot); } let r_norm = norm(&Slice { data: &*scratch }); - if r_norm > T::zero() || r_norm == T::zero() { - Ok(r_norm) - } else { + // `x - x` is `0` for every finite `x` and `NaN` for `NaN`/±infinity — the only values + // unequal to themselves. `r_norm > T::zero()` alone would accept `+Infinity`. + let probe = r_norm.sub(r_norm); + #[allow(clippy::eq_op)] + let non_finite = probe != probe; + if non_finite { Err(ConvergenceError::NonFinite) + } else { + Ok(r_norm) } } diff --git a/src/krylov/power_iteration.rs b/src/krylov/power_iteration.rs index edaefb7..6eb999c 100644 --- a/src/krylov/power_iteration.rs +++ b/src/krylov/power_iteration.rs @@ -189,22 +189,26 @@ where /// Scales `v` in place to unit Euclidean length. /// /// Errors with `ZeroVector` when `‖v‖` is exactly zero — that covers the zero vector and the -/// empty (`n == 0`) vector. Errors with `NonFinite` when `‖v‖` is neither strictly positive -/// nor exactly zero, which is only possible when it's `NaN`: a norm poisoned by non-finite -/// input compares false against every ordinary value, including `0`, so it fails both checks -/// rather than being divided by. +/// empty (`n == 0`) vector. Errors with `NonFinite` when `‖v‖` is `NaN` or infinite: `x - x` is +/// `0` for every finite `x` and `NaN` for `NaN`/±infinity, the only values unequal to +/// themselves, so that self-subtraction is checked before the ordinary comparisons below — +/// `+Infinity > 0` is `true`, so a naive `length > T::zero()` check would silently divide by +/// an infinite length and zero out an overflowed (not actually zero) vector instead of erroring. pub(super) fn normalize(v: &mut [T]) -> Result<(), ConvergenceError> { let length = norm(&Slice { data: &*v }); - if length > T::zero() { + let probe = length.sub(length); + #[allow(clippy::eq_op)] + let non_finite = probe != probe; + if non_finite { + Err(ConvergenceError::NonFinite) + } else if length == T::zero() { + Err(ConvergenceError::ZeroVector) + } else { let inv = T::one().div(length); for slot in v.iter_mut() { *slot = slot.mul(inv); } Ok(()) - } else if length == T::zero() { - Err(ConvergenceError::ZeroVector) - } else { - Err(ConvergenceError::NonFinite) } } diff --git a/tests/edge_cases/krylov.rs b/tests/edge_cases/krylov.rs index d7932e6..0fff9cb 100644 --- a/tests/edge_cases/krylov.rs +++ b/tests/edge_cases/krylov.rs @@ -1,11 +1,15 @@ //! Curated, fixed-matrix edge cases for `power_iteration`, `inverse_power_iteration`, -//! `lanczos`, and `arnoldi`: degenerate spectra, singular shifts, non-finite inputs, +//! `lanczos`, `arnoldi`, and `gmres`: degenerate spectra, singular shifts, non-finite inputs, //! dimension mismatches, and the zero vector — cases the property harness deliberately never //! generates. use rustebra::krylov::{ ConvergenceError, arnoldi, inverse_power_iteration, lanczos, power_iteration, }; +#[cfg(feature = "alloc")] +use rustebra::krylov::gmres; +#[cfg(feature = "alloc")] +use rustebra::sparse::CsrMatrix; use rustebra::storage::{Basis, StaticStorage}; use crate::common::{ @@ -509,6 +513,60 @@ fn arnoldi_zero_initial_vector_is_a_zero_vector_error() { assert_eq!(result, Err(ConvergenceError::ZeroVector)); } +/// Dense row-major `a` as a `CsrMatrix`, storing every entry (including exact zeros). +#[cfg(feature = "alloc")] +fn csr_from_dense(a: &[f64], n: usize) -> CsrMatrix { + let mut row_ptr = vec![0_u32]; + let mut col_indices = vec![]; + let mut values = vec![]; + for r in 0..n { + for c in 0..n { + col_indices.push(c as u32); + values.push(a[r * n + c]); + } + row_ptr.push(col_indices.len() as u32); + } + CsrMatrix::new(n, n, row_ptr, col_indices, values).unwrap() +} + +#[test] +#[cfg(feature = "alloc")] +fn gmres_overflowing_residual_is_non_finite_not_a_spurious_breakdown() { + // b's entries are individually finite, but `‖r‖²` (an unscaled sum of squares) overflows + // f64 before the sqrt: a naive `r_norm > 0.0` check accepts `+Infinity`, silently + // normalizes it down to a zero vector, and mislabels the resulting degenerate Arnoldi + // step as `Breakdown` instead of reporting the real cause, `NonFinite`. + let a = csr_from_dense(&[1.0, 0.0, 0.0, 1.0], 2); + let b = [1e200, 1e200]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 4]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + let result = gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch); + + assert_eq!(result, Err(ConvergenceError::NonFinite)); +} + +#[test] +#[cfg(feature = "alloc")] +fn gmres_non_finite_matrix_entry_is_an_error_not_a_panic() { + for poison in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + let a = csr_from_dense(&[2.0, poison, 0.0, 1.0], 2); + let b = [1.0, 1.0]; + let x0 = [0.0, 0.0]; + let mut out_x = [0.0; 2]; + let mut buffer = [0.0; 4]; + let mut basis = Basis::::new(&mut buffer, 2).unwrap(); + let mut scratch = [0.0; 2]; + + let result = gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch); + + assert!(result.is_err(), "gmres accepted a {poison} entry: {result:?}"); + } +} + #[test] fn arnoldi_dimension_mismatch_is_an_error() { // 4 entries can't be a 3x3 matrix. diff --git a/tests/property/diff/gmres.rs b/tests/property/diff/gmres.rs new file mode 100644 index 0000000..b8182dc --- /dev/null +++ b/tests/property/diff/gmres.rs @@ -0,0 +1,89 @@ +#![cfg(feature = "alloc")] + +//! Differential property test: `gmres` against `nalgebra`'s direct LU solve on random +//! diagonally dominant systems, comparing the solution vector itself rather than only the +//! residual the solver reports about its own answer. + +use super::approx_eq; +use nalgebra::{DMatrix, DVector}; +use proptest::prelude::*; +use rustebra::krylov::gmres; +use rustebra::sparse::CsrMatrix; +use rustebra::storage::Basis; + +const N: usize = 4; +const TOL: f64 = 1e-6; + +/// Dense row-major `a` as a `CsrMatrix`, storing every entry (including exact zeros). +fn csr_from_dense(a: &[f64; N * N]) -> CsrMatrix { + let mut row_ptr = vec![0_u32]; + let mut col_indices = vec![]; + let mut values = vec![]; + for r in 0..N { + for c in 0..N { + col_indices.push(c as u32); + values.push(a[r * N + c]); + } + row_ptr.push(col_indices.len() as u32); + } + CsrMatrix::new(N, N, row_ptr, col_indices, values).unwrap() +} + +prop_compose! { + /// A row-major, strictly diagonally dominant `N x N` matrix, guaranteeing a well-posed + /// system nalgebra's direct solve and GMRES(`N`) both converge on. + fn diagonally_dominant_matrix()( + off_diagonal in prop::array::uniform16(-1.0..1.0f64), + diagonal_boost in prop::array::uniform4(5.0..10.0f64), + ) -> [f64; N * N] { + let mut a = off_diagonal; + for r in 0..N { + a[r * N + r] = 0.0; + } + let mut row_sums = [0.0; N]; + for r in 0..N { + row_sums[r] = (0..N).map(|c| a[r * N + c].abs()).sum(); + } + for r in 0..N { + a[r * N + r] = row_sums[r] + diagonal_boost[r]; + } + a + } +} + +proptest! { + /// GMRES(N) — a full-dimension restart size, so no restart is ever needed for a + /// diagonally dominant system — must land on the same solution nalgebra's direct LU solve + /// produces for `A x = b`. + #[test] + fn solution_matches_nalgebra_direct_solve( + a in diagonally_dominant_matrix(), + b in prop::array::uniform4(-10.0..10.0f64), + ) { + let m = csr_from_dense(&a); + let x0 = [0.0; N]; + let mut out_x = [0.0; N]; + let mut buffer = [0.0; N * N]; + let mut basis = Basis::::new(&mut buffer, N).unwrap(); + let mut scratch = [0.0; N]; + + gmres(&m, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).unwrap(); + + let a_na = DMatrix::from_row_slice(N, N, &a); + let b_na = DVector::from_row_slice(&b); + let x_na = a_na + .lu() + .solve(&b_na) + .expect("diagonally dominant matrices are always invertible"); + + for i in 0..N { + prop_assert!( + approx_eq(out_x[i], x_na[i], TOL), + "x[{}]: ours={} vs nalgebra={}", + i, + out_x[i], + x_na[i] + ); + } + } +} diff --git a/tests/property/diff/mod.rs b/tests/property/diff/mod.rs index df07f1a..cdb51bd 100644 --- a/tests/property/diff/mod.rs +++ b/tests/property/diff/mod.rs @@ -4,6 +4,7 @@ mod arnoldi; mod cholesky; +mod gmres; mod lanczos; mod lu; mod qr; From 432ae7421d907f2741715c449c1620b29c01e353 Mon Sep 17 00:00:00 2001 From: Eli Date: Thu, 16 Jul 2026 22:56:56 -0300 Subject: [PATCH 3/4] Add Arnoldi and GMRES(m) book pages Document both under Krylov Methods, linking them from the section README and SUMMARY nav. Co-Authored-By: Claude Sonnet 5 --- docs/SUMMARY.md | 2 ++ docs/book/08-krylov/README.md | 5 +++ docs/book/08-krylov/arnoldi.md | 48 ++++++++++++++++++++++++++++ docs/book/08-krylov/gmres.md | 57 ++++++++++++++++++++++++++++++++++ 4 files changed, 112 insertions(+) create mode 100644 docs/book/08-krylov/arnoldi.md create mode 100644 docs/book/08-krylov/gmres.md diff --git a/docs/SUMMARY.md b/docs/SUMMARY.md index d2082be..66fae45 100644 --- a/docs/SUMMARY.md +++ b/docs/SUMMARY.md @@ -32,6 +32,8 @@ - [Power Iteration](book/08-krylov/power-iteration.md) - [Conjugate Gradient](book/08-krylov/conjugate-gradient.md) - [Lanczos Iteration](book/08-krylov/lanczos.md) + - [Arnoldi Iteration](book/08-krylov/arnoldi.md) + - [GMRES(m)](book/08-krylov/gmres.md) - [Eigenvalues & Eigenvectors](book/09-eigenvalues.md) - [Numerical Stability](book/12-numerical-stability/README.md) - [NaN/Inf Policy](book/12-numerical-stability/nan-inf-policy.md) diff --git a/docs/book/08-krylov/README.md b/docs/book/08-krylov/README.md index db0fddd..5bc93d8 100644 --- a/docs/book/08-krylov/README.md +++ b/docs/book/08-krylov/README.md @@ -7,6 +7,11 @@ - **Conjugate Gradient (CG)** for solving symmetric positive-definite linear systems - **Lanczos iteration**, which builds an orthonormal basis of a Krylov subspace and the symmetric tridiagonal matrix a symmetric operator projects onto within it +- **Arnoldi iteration**, the non-symmetric counterpart to Lanczos, which builds an orthonormal + basis of a Krylov subspace and the upper Hessenberg matrix a general operator projects onto + within it +- **GMRES(m)**, a restarted iterative solver for general (possibly non-symmetric) linear + systems, built on top of Arnoldi iteration Unlike the direct decompositions in [Decompositions](../06-decompositions/README.md), these refine an estimate (or a basis) over many iterations and can fail to converge — or, for diff --git a/docs/book/08-krylov/arnoldi.md b/docs/book/08-krylov/arnoldi.md new file mode 100644 index 0000000..e018b4b --- /dev/null +++ b/docs/book/08-krylov/arnoldi.md @@ -0,0 +1,48 @@ +# Arnoldi Iteration + +`arnoldi` reduces a general, row-major `n x n` matrix `a` to upper Hessenberg form over a +`K`-dimensional Krylov subspace: starting from a normalized `v0`, it builds an orthonormal +basis `Q` of `span{v0, a*v0, ..., a^{K-1}*v0}` and returns the projection `H = Qᵗ * a * Q`, +which is upper Hessenberg. Unlike [Lanczos Iteration](lanczos.md), `a` need not be symmetric — +there is no three-term recurrence to exploit, so every step orthogonalizes the candidate +vector against the *entire* basis built so far, by modified Gram-Schmidt, rather than just the +two previous vectors. + +```rust +{{#include ../../../examples/krylov/arnoldi.rs}} +``` + +## Orthogonalization + +Modified Gram-Schmidt (subtracting each projection immediately, rather than computing all +projections against the original candidate and subtracting them at the end) is the standard +trade for Arnoldi: markedly more stable than classical Gram-Schmidt at the same `O(K * n)` +per-step cost, though still less stable than Householder Arnoldi, which trades that extra +stability for losing the explicit basis vectors the Krylov projection needs. No +reorthogonalization pass is added on top, unlike Lanczos's full reorthogonalization — Arnoldi's +every-vector orthogonalization does not erode as quickly as Lanczos's three-term recurrence +does, so a second pass doesn't earn its keep the same way. + +## Breakdown is success, not failure + +When the candidate vector's norm falls to (numerically) zero relative to `‖a * q_j‖` after +orthogonalization, `q_0, ..., q_j` already span an invariant subspace of `a`: there's no new +direction to extend the basis with, but the vectors and the leading block of `H` already built +are exact and useful. This is reported as `Ok((h, reached))` with `reached < K`, not an error — +unlike Lanczos, which reports the analogous condition as `ConvergenceError::Breakdown`. The +difference is what the caller does next: a Lanczos caller that requested `K` vectors and got +fewer has nothing it can use without changing its request, while GMRES, built on top of +Arnoldi, can solve directly in the smaller subspace `Ok` reports — the exact subspace is often +exactly where the true solution already lives. Callers that do need the full `K` vectors +distinguish this case from a complete run by checking `reached < K`. + +## Gotchas + +- `K` is a `const` generic on the caller's `Basis` buffer, not a runtime parameter — see + [Krylov Basis-Size Const-Generic Convention](../../specs/krylov-basis-size-const-generics.md). + `K > n` is a `DimensionMismatch`. +- Both `ConvergenceError::ZeroVector` and `ConvergenceError::NonFinite` on `v0` are checked up + front, even when `K == 0` means no basis vector is ever written. +- `tol` has no auto-computed default — see + [Krylov Tolerance and Convergence Criteria](../../specs/krylov-tolerance-and-convergence.md). + A `tol` of `0` detects only exact breakdown. diff --git a/docs/book/08-krylov/gmres.md b/docs/book/08-krylov/gmres.md new file mode 100644 index 0000000..f90820a --- /dev/null +++ b/docs/book/08-krylov/gmres.md @@ -0,0 +1,57 @@ +# GMRES(m) + +`gmres` solves the general (possibly non-symmetric) linear system `A x = b` via restarted +GMRES, GMRES(`M`): unlike [Conjugate Gradient](conjugate-gradient.md), `A` need not be +symmetric positive-definite. Each restart cycle runs [Arnoldi Iteration](arnoldi.md) from the +current residual to build an `M`-dimensional Krylov basis, solves the resulting small +least-squares problem via Givens rotations, and updates `x` — restarting from the improved +iterate until either the residual meets `tol` or `max_restarts` cycles are exhausted. + +`A` is supplied as a sparse linear operator rather than a dense matrix: applying it never +allocates, so restart cycles reuse the same workspace (`out_x`, `basis`, `scratch`) throughout. + +```rust +{{#include ../../../examples/krylov/gmres.rs}} +``` + +## Algorithm + +Each restart cycle: + +1. Computes the residual `r = b - A x` and its norm `β = ‖r‖`, returning `Ok` immediately if + `β <= tol`. +2. Runs Arnoldi iteration from `q_0 = r / β`, building an orthonormal basis `Q` of up to `M` + vectors and the upper Hessenberg projection `H`, stopping early (before `M` steps) on + breakdown — an invariant subspace found before the basis filled up, the same + "success, not failure" case documented on [Arnoldi Iteration](arnoldi.md). +3. Solves `min_y ‖β e_1 - H y‖` via incremental Givens rotations, then updates `x <- x + Q y`. + +## Convergence + +GMRES's residual norm decreases monotonically within a cycle (each additional basis vector can +only improve the least-squares fit) and never increases across a restart, because restarting +recomputes the same residual the next cycle continues from. It is not guaranteed to decrease +*strictly* every cycle, though: a starting vector aligned with an invariant subspace the +operator doesn't expand (breakdown on the very first Arnoldi step) leaves `x` unchanged, and the +iteration stagnates. Restarting also discards the larger Krylov subspace full (non-restarted) +GMRES would have kept building, so GMRES(`M`) can converge slower, or stagnate on problems full +GMRES would resolve — the restart budget `max_restarts` bounds the cost of that risk rather than +eliminating it. + +## Gotchas + +- The restart size `M` is a `const` generic on the caller's `Basis` buffer, not a runtime + parameter — see + [Krylov Basis-Size Const-Generic Convention](../../specs/krylov-basis-size-const-generics.md). + `M > n` (the operator's dimension) is a `DimensionMismatch`. +- `M == 0` never reduces a nonzero residual: no basis vector can be built, so every restart + cycle is a no-op check against `tol`, and a nonzero residual exhausts `max_restarts` without + ever moving `x`. +- `tol` has no auto-computed default — see + [Krylov Tolerance and Convergence Criteria](../../specs/krylov-tolerance-and-convergence.md). +- Non-finite residuals or Arnoldi iterates (`NaN` or infinite, including values that overflow + `f64` mid-computation) return `ConvergenceError::NonFinite` rather than silently producing a + wrong answer. +- A zero pivot in the small least-squares system built from `H` — a coincidental exact + singularity in the projected system that Arnoldi's own breakdown test didn't already catch — + returns `ConvergenceError::Breakdown`. From 8cf8afd56f8a62b5a9e51edd3f794338596cf45d Mon Sep 17 00:00:00 2001 From: Eli Date: Fri, 17 Jul 2026 00:27:50 -0300 Subject: [PATCH 4/4] Fix import formatting after rebase conflict resolution Co-Authored-By: Claude Sonnet 5 --- tests/edge_cases/krylov.rs | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/edge_cases/krylov.rs b/tests/edge_cases/krylov.rs index 0fff9cb..c3ef5e8 100644 --- a/tests/edge_cases/krylov.rs +++ b/tests/edge_cases/krylov.rs @@ -3,12 +3,12 @@ //! dimension mismatches, and the zero vector — cases the property harness deliberately never //! generates. +#[cfg(feature = "alloc")] +use rustebra::krylov::gmres; use rustebra::krylov::{ ConvergenceError, arnoldi, inverse_power_iteration, lanczos, power_iteration, }; #[cfg(feature = "alloc")] -use rustebra::krylov::gmres; -#[cfg(feature = "alloc")] use rustebra::sparse::CsrMatrix; use rustebra::storage::{Basis, StaticStorage}; @@ -563,7 +563,10 @@ fn gmres_non_finite_matrix_entry_is_an_error_not_a_panic() { let result = gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch); - assert!(result.is_err(), "gmres accepted a {poison} entry: {result:?}"); + assert!( + result.is_err(), + "gmres accepted a {poison} entry: {result:?}" + ); } }