From e8a6b0ad043559ac762734febf40e412fafbeeab Mon Sep 17 00:00:00 2001 From: Kristof Schroeder Date: Wed, 23 Sep 2026 16:59:42 +0200 Subject: [PATCH 1/2] perf(lsmr): hold a tolerance stop to its true residual and restart from a refuted one MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The exit audit recomputed `b − Ax`, `Aᵀr` and `M⁻¹Aᵀr` at every tolerance stop, about 1.5 iterations per solve. A stop is now held to one true-residual evaluation, `‖b − A x‖` against the recurrence's own `‖r_k‖`, and a refuted stop restarts the stream from its iterate, at most twice, before `FalseConvergence`. The `Aᵀr` and `M⁻¹` legs, the certificate and its warm re-base are gone; a warm exit's normal-equation report still divides by `‖Âᵀb‖`, one `Aᵀ` and one `M⁻¹` apply, fixed for the solve so a restarted pass reports against the same quantity. `M⁻¹` must now be nonsingular: a direction it annihilates lies outside the Krylov space, and the singular-preconditioner test flips from refused to certified. Hot loop, measured at d49ffd6 before the stack was split, 8 threads, min of interleaved rounds against v0.3.0: 2M×3FE 10 iterations +0.9% (main +18.1%), 8M×3FE 14 iterations +1.5% (main +13.5%). Smoke suite 94/94 converged. --- CHANGELOG.md | 7 +- crates/schwarz-precond/src/lsmr.rs | 286 ++++++++++-------- crates/schwarz-precond/src/lsmr/bidiag.rs | 234 ++++++-------- crates/schwarz-precond/src/lsmr/recurrence.rs | 80 ++--- .../schwarz-precond/src/lsmr/tests/audit.rs | 183 ++++++----- .../schwarz-precond/src/lsmr/tests/ladder.rs | 22 -- .../schwarz-precond/src/lsmr/tests/solve.rs | 15 +- 7 files changed, 398 insertions(+), 429 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 72361754..df65017b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project follows [Semantic Versioning](https://semver.org/). ### Changed +- A tolerance stop is checked with one true-residual evaluation; a refuted one restarts from its iterate, at most twice, before `LsmrStopReason::FalseConvergence`. - Python `Effect` lists no longer deep-copy their level and slope buffers during design extraction; the binding borrows them while building off-GIL (#358). - Categorical `u32` labels need not be zero-based or contiguous: `Design` compacts observed labels to internal positions and `CoefficientLayout`/`CoefficientAddress` translate back, so gaps in sparse label ranges are neither allocated nor solved for (#228, #268). - **BREAKING:** Rust `Solver::new` takes `weights: Option<&[f64]>` instead of `Option>`, retaining only `W^{1/2}` in internal observation order. One-shot `solve`/`solve_batch` weights are unchanged. @@ -18,7 +19,7 @@ and this project follows [Semantic Versioning](https://semver.org/). - **BREAKING:** `schwarz_precond::mlsmr` takes an `MlsmrOptions` in place of its trailing `local_size`. - **BREAKING:** `LsmrStopReason` gains `Escalated` and `WarmStartExact`, breaking exhaustive `match`es. - A warm start that already solves the system reports `WarmStartExact` instead of `ZeroRhs`. -- The LSMR true-residual audit of a warm-started stop measures the total solution against the original `b`, anchoring its normal-equation leg to `‖Aᵀb‖` rather than the restart's initial residual; with `b = 0`, tolerances measure against `‖b − A x₀‖`. +- A warm-started LSMR stop measures the total solution against the original `b`; with `b = 0`, tolerances measure against `‖b − A x₀‖`. - **BREAKING:** `ScalingConfig::max_sweeps` is now `max_iterations`, and `BuildWarning::UnscalableComponent` reports `iterations` in place of `sweeps`; the dominance certificate runs reduced CG, not relaxation sweeps. - **BREAKING:** The serialized `Preconditioner` wire format moved v12 → v17 (approx-chol 0.5, full construction config, build duration, `LocalSolverConfig::ridge`, and the built map's own Schwarz description in place of the strategy enum); 0.3.0 bytes no longer decode. - **BREAKING:** Coefficient addresses use caller-visible `u32` factor labels rather than internal `usize` level positions, affecting `CoefficientAddress::level` and the accepted range of Python coefficient-layout and unidentified-direction levels. @@ -39,8 +40,8 @@ and this project follows [Semantic Versioning](https://semver.org/). - A design carrying varying slopes on two distinct factors could fail preconditioner construction with `matrix is not symmetric`, when rounding left the two triangles of the exact Schur complement unequal (#229). - A `design` that is neither a 2-D `uint32` array nor a list of `Effect` raised `ValueError` where the documented type is `TypeError`, and `AdditiveSchwarz` accepted a wrong-type `local_solver` at construction, deferring the `TypeError` to solve time (#248). -- LSMR no longer certifies a stop it has not solved. Tolerance stops are audited against the true residual and outside the preconditioner's metric, a non-finite `α`, `β`, `⟨v, Mv⟩`, or `‖b‖` fails with `SolveError::InvalidInput`, and an overflowing or subnormal `‖A‖` no longer zeroes the normal-equation ratio; a failed check reports `LsmrStopReason::FalseConvergence` with `converged = false` (#290, #297, #303, #362). -- A warm-started solve measures its residuals against the original `b`, including when it exhausts its iteration budget, and no longer fails a converged solve when the cold audit's metric product overflows at `‖Aᵀb‖ ≳ 1e154`. +- LSMR no longer certifies a stop it has not solved. Tolerance stops are audited against the true residual, a non-finite `α`, `β`, `⟨v, Mv⟩`, or `‖b‖` fails with `SolveError::InvalidInput`, and an overflowing or subnormal `‖A‖` no longer zeroes the normal-equation ratio; a failed check reports `LsmrStopReason::FalseConvergence` with `converged = false` (#290, #297, #303, #362). +- A warm-started solve measures its residuals against the original `b`, including at a budget stop and where `‖Aᵀb‖²` overflows. - An `α` or `β` whose square underflows is recovered from a scaled norm, where LSMR reported `x = 0` converged. - LSMR's solution update no longer drops out when `‖A‖` is small (`≈ 1e-9`) or beyond `1e±154`, where it reported `x = 0` converged. - LSMR's residual estimate is its own `‖r_k‖` rather than LSQR's smaller `|φ̄_k|`, which let `ResidualTolerance` fire before the tolerance was met. diff --git a/crates/schwarz-precond/src/lsmr.rs b/crates/schwarz-precond/src/lsmr.rs index 740477e3..5b2453cd 100644 --- a/crates/schwarz-precond/src/lsmr.rs +++ b/crates/schwarz-precond/src/lsmr.rs @@ -15,7 +15,8 @@ use std::borrow::Cow; use crate::{Operator, SolveError}; use bidiag::{ - residual_into, BidiagStep, Bidiagonalization, Certificate, GolubKahan, ModifiedGolubKahan, + axpby, metric_gradient_norm, residual_into, BidiagStep, Bidiagonalization, GolubKahan, + ModifiedGolubKahan, }; use recurrence::{ConvergenceCriteria, LsmrRecurrenceState, RotationStep, SolutionState, Stop}; @@ -42,13 +43,14 @@ pub(crate) fn vec_norm(v: &[f64]) -> f64 { pub struct LsmrResult { /// Solution vector. pub x: Vec, - /// Whether the solver converged within the tolerance. + /// Whether the solver converged within the tolerance. A tolerance stop is held against + /// `‖b − A x‖`, which audits the recurrence and not `M`'s nonsingularity. pub converged: bool, /// Total number of iterations performed. pub iterations: usize, - /// Final residual norm estimate `‖b − A x‖`. + /// `‖b − A x‖`: recomputed at a tolerance stop, the recurrence's estimate at any other. pub residual_norm: f64, - /// Per-run relative normal-equation residual, measured in `M`'s metric when preconditioned. + /// `‖Âᵀ(b − A x)‖ / ‖Âᵀb‖`, measured in `M`'s metric when preconditioned. pub normal_eq_residual: f64, /// Reason the solver stopped. pub stop_reason: LsmrStopReason, @@ -67,7 +69,8 @@ pub enum LsmrStopReason { NormalEquationTolerance, /// The warm start already solved the system: `b − A x0` was exactly zero. WarmStartExact, - /// A tolerance stop refuted by the true-residual audit; the recurrence estimates had collapsed. + /// A tolerance stop the true residual refuted, as were the restarts from it: the recurrence + /// estimates had collapsed. FalseConvergence, /// The iteration budget was exhausted before convergence. MaxIterations, @@ -245,6 +248,8 @@ pub fn lsmr( } /// Preconditioned LSMR with `M ≈ AᵀA` and one `M⁻¹` application per iteration. +/// +/// `M⁻¹` must be nonsingular; a direction it annihilates can be reported converged unsolved. pub fn mlsmr( operator: &A, b: &[f64], @@ -313,145 +318,186 @@ pub fn mlsmr( }); } + // `‖Aᵀb‖` must be taken before the stream exists: the stream's own query clobbers `v₁`. + let metric = match warm_start { + None => None, + Some(_) => Some(metric_gradient_norm( + operator, + preconditioner, + b, + &mut vec![0.0; n], + &mut vec![0.0; n], + )?), + }; let (bidiag, step1) = ModifiedGolubKahan::init(operator, preconditioner, &rhs, rhs_norm, local_size)?; + let warm_start = warm_start.zip(metric).map(|(x0, metric)| WarmStart { + x0, + reference: NormalEqReference::warm(metric, step1), + }); let reference_norm = if b_norm > 0.0 { b_norm } else { rhs_norm }; let criteria = ConvergenceCriteria::new(reference_norm, tol); - lsmr_from_bidiag( - bidiag, - step1, - b, - warm_start, - criteria, - maxiter, - escalation.map(|policy| policy.handler()), - ) + lsmr_from_bidiag(bidiag, step1, b, warm_start, criteria, maxiter, escalation) } -/// A warm-started stream measures its references from the correction's own residual, so a cold -/// certificate restores them in terms of `rhs`. -fn audit_against_rhs( - bidiag: &mut B, - x: &[f64], - rhs: &[f64], - x0: Option<&[f64]>, -) -> Result { - let cold = x0 - .map(|_| bidiag.certify(&vec![0.0; x.len()], rhs)) - .transpose()?; - let mut cert = bidiag.certify(x, rhs)?; - if let Some(cold) = &cold { - cert.rebase(cold); +/// Restarts from a refuted tolerance stop before the solve is refused outright. +const MAX_RESTARTS: usize = 2; + +/// `‖Aᵀb‖` in the stream's metric: the denominator every reported normal-equation residual +/// divides by, fixed for the solve so a restart reports against the same quantity as pass 1. +#[derive(Clone, Copy)] +struct NormalEqReference(f64); + +impl NormalEqReference { + /// A cold stream's first `(α₁, β₁)` is `‖Aᵀb‖` already. + fn cold(step1: BidiagStep) -> Self { + Self((step1.alpha * step1.beta).abs().max(f64::MIN_POSITIVE)) + } + + /// A reference that carries no information leaves the stream's own `ζ̄₀` to divide. + fn warm(metric: f64, step1: BidiagStep) -> Self { + if metric > 0.0 && metric.is_finite() { + Self(metric) + } else { + Self::cold(step1) + } } - Ok(cert) + + fn relative(self, estimate: f64) -> f64 { + estimate / self.0 + } +} + +/// A warm start, paired with the `‖Aᵀb‖` it costs the solve: the stream is seeded from +/// `b − A x₀`, so the first `(α₁, β₁)` measures that residual and not `b`. +struct WarmStart<'a> { + x0: &'a [f64], + reference: NormalEqReference, } /// Runs the LSMR recurrences over a preconditioner-specific bidiagonalization stream. /// -/// The stream solves for the correction to `x0`; results and the audit are in terms of `b`. +/// The stream solves for the correction to `x0`; results are in terms of `b`. A tolerance stop +/// is held against the true residual, and a stop it refutes restarts the stream from the iterate. fn lsmr_from_bidiag( mut bidiag: B, - step1: BidiagStep, + mut step1: BidiagStep, b: &[f64], - x0: Option<&[f64]>, + warm_start: Option>, criteria: ConvergenceCriteria, maxiter: usize, - mut escalation: Option>, + escalation: Option<&dyn EscalationPolicy>, ) -> Result { let n = bidiag.v().len(); - let total = |mut x: Vec| { - if let Some(x0) = x0 { - for (xi, &x0i) in x.iter_mut().zip(x0) { - *xi += x0i; + let reference = warm_start + .as_ref() + .map_or_else(|| NormalEqReference::cold(step1), |w| w.reference); + let mut base: Option> = warm_start.map(|w| Cow::Borrowed(w.x0)); + let mut iterations = 0; + let mut restarts = 0; + loop { + // A metric reporting no gradient at all may be hiding one outside itself. + if step1.alpha == 0.0 { + let x = base.map_or_else(|| vec![0.0; n], Cow::into_owned); + let (converged, normal_eq_residual) = match bidiag.hidden_gradient()? { + None => (true, 0.0), + Some(per_unit_residual) => { + let normar = step1.beta * per_unit_residual; + let plain = bidiag.plain_gradient(b)?; + let a_norm_below = plain / vec_norm(b).max(f64::MIN_POSITIVE); + let informative = plain > 0.0 && plain.is_finite(); + ( + criteria.corroborates(step1.beta, normar, a_norm_below), + if informative { normar / plain } else { 1.0 }, + ) + } + }; + return Ok(LsmrResult { + x, + converged, + iterations, + residual_norm: step1.beta, + normal_eq_residual, + stop_reason: if converged { + LsmrStopReason::InitialNormalEquationResidualZero + } else { + LsmrStopReason::FalseConvergence + }, + }); + } + + // A pass reports its drop from its own ζ̄₀, so a restarted one needs a handler that agrees. + let mut escalation = escalation.map(EscalationPolicy::handler); + let mut convergence = criteria.start(step1.alpha); + let mut recurrence = LsmrRecurrenceState::init(step1); + let mut solution = SolutionState::init(bidiag.v()); + let mut prev_rot = RotationStep::initial(); + let stop_reason = 'pass: { + while iterations < maxiter { + iterations += 1; + let step = bidiag.step()?; + convergence.observe(step); + let curr_rot = recurrence.step(step); + solution.update(bidiag.v(), curr_rot, prev_rot); + + // The tolerance test catches breakdown when the residual recurrences collapse. + match convergence.check(&recurrence) { + Stop::Continue => {} + Stop::ResidualTolerance => break 'pass LsmrStopReason::ResidualTolerance, + Stop::NormalEquationTolerance => { + break 'pass LsmrStopReason::NormalEquationTolerance + } + } + if let Some(rule) = escalation.as_deref_mut() { + let progress = Progress { + iteration: iterations, + normal_eq_residual: recurrence.relative_normal_eq_residual(), + }; + if rule.should_escalate(progress) { + break 'pass LsmrStopReason::Escalated; + } + } + prev_rot = curr_rot; } + LsmrStopReason::MaxIterations + }; + + let mut x = solution.into_x(); + if let Some(base) = &base { + axpby(&mut x, base, 1.0, 1.0); } - x - }; - // A metric reporting no gradient at all may be hiding one outside itself. - if step1.alpha == 0.0 { - let x = total(vec![0.0; n]); - let cert = bidiag.certify(&x, b)?; - let converged = criteria.corroborated(&cert, || bidiag.operator_norm_below(b))?; - let (residual_norm, normal_eq_residual) = cert.residuals(); - return Ok(LsmrResult { + let mut result = LsmrResult { x, - converged, - iterations: 0, - residual_norm, - normal_eq_residual, - stop_reason: if converged { - LsmrStopReason::InitialNormalEquationResidualZero - } else { - LsmrStopReason::FalseConvergence - }, - }); - } - - let mut convergence = criteria.start(step1.alpha); - let mut recurrence = LsmrRecurrenceState::init(step1); - let mut solution = SolutionState::init(bidiag.v()); - let mut prev_rot = RotationStep::initial(); - - let (iterations, stop_reason) = 'run: { - for itn in 1..=maxiter { - let step = bidiag.step()?; - convergence.observe(step); - let curr_rot = recurrence.step(step); - solution.update(bidiag.v(), curr_rot, prev_rot); - - // The tolerance test catches breakdown when the residual recurrences collapse. - if let Some(stop_reason) = match convergence.check(&recurrence) { - Stop::Continue => None, - Stop::ResidualTolerance => Some(LsmrStopReason::ResidualTolerance), - Stop::NormalEquationTolerance => Some(LsmrStopReason::NormalEquationTolerance), - } { - let x = total(solution.into_x()); - let cert = audit_against_rhs(&mut bidiag, &x, b, x0)?; - let converged = convergence.certified(&cert); - let (residual_norm, normal_eq_residual) = cert.residuals(); - return Ok(LsmrResult { - x, - converged, - iterations: itn, - residual_norm, - normal_eq_residual, - stop_reason: if converged { - stop_reason - } else { - LsmrStopReason::FalseConvergence - }, - }); + converged: false, + iterations, + residual_norm: recurrence.residual_estimate(), + normal_eq_residual: reference.relative(recurrence.normal_eq_residual_estimate()), + stop_reason, + }; + // Only a tolerance stop claims convergence, so only it is worth a true-residual evaluation. + if matches!( + stop_reason, + LsmrStopReason::ResidualTolerance | LsmrStopReason::NormalEquationTolerance + ) { + let residual_norm = bidiag.residual_norm(&result.x, b)?; + result.converged = criteria.corroborates_residual(residual_norm, result.residual_norm); + if !result.converged + && residual_norm.is_finite() + && restarts < MAX_RESTARTS + && iterations < maxiter + { + restarts += 1; + step1 = bidiag.restart(residual_norm)?; + base = Some(Cow::Owned(result.x)); + continue; } - if let Some(rule) = escalation.as_deref_mut() { - let progress = Progress { - iteration: itn, - normal_eq_residual: recurrence.relative_normal_eq_residual(), - }; - if rule.should_escalate(progress) { - break 'run (itn, LsmrStopReason::Escalated); - } + if !result.converged { + result.stop_reason = LsmrStopReason::FalseConvergence; } - prev_rot = curr_rot; + result.residual_norm = residual_norm; } - (maxiter, LsmrStopReason::MaxIterations) - }; - let x = total(solution.into_x()); - // A warm run's recurrence measures the restart, so only an audit reports against `b`. - let (residual_norm, normal_eq_residual) = match x0 { - Some(_) => audit_against_rhs(&mut bidiag, &x, b, x0)?.residuals(), - None => ( - recurrence.residual_estimate(), - recurrence.relative_normal_eq_residual(), - ), - }; - Ok(LsmrResult { - x, - converged: false, - iterations, - residual_norm, - normal_eq_residual, - stop_reason, - }) + return Ok(result); + } } fn validate_lsmr_inputs( diff --git a/crates/schwarz-precond/src/lsmr/bidiag.rs b/crates/schwarz-precond/src/lsmr/bidiag.rs index 6e249213..d5b8bc49 100644 --- a/crates/schwarz-precond/src/lsmr/bidiag.rs +++ b/crates/schwarz-precond/src/lsmr/bidiag.rs @@ -46,7 +46,7 @@ fn axpy_with_norm(y: &mut [f64], x: &[f64], scale: f64) -> f64 { /// `y = alpha * x + beta * y`. Parallel above the threshold. #[inline] -fn axpby(y: &mut [f64], x: &[f64], alpha: f64, beta: f64) { +pub(super) fn axpby(y: &mut [f64], x: &[f64], alpha: f64, beta: f64) { debug_assert_eq!(x.len(), y.len()); let seq = |y_c: &mut [f64], x_c: &[f64]| { for (yi, &xi) in y_c.iter_mut().zip(x_c.iter()) { @@ -146,10 +146,7 @@ pub(super) fn residual_into( u: &mut [f64], ) -> Result { operator.apply(x, u)?; - for (ui, &bi) in u.iter_mut().zip(rhs) { - *ui = bi - *ui; - } - Ok(super::vec_norm(u)) + Ok(axpy_with_norm(u, rhs, -1.0)) } /// Scales `u` to unit length given its already-computed `β = ‖u‖`; a zero `u` stays zero. @@ -196,6 +193,12 @@ impl WindowRing { (0..count).map(move |i| (start + i) % cap) } + /// Forget every stored vector: a restart begins a new sequence. + fn clear(&mut self) { + self.next = 0; + self.count = 0; + } + /// Reserve the next write slot, advancing the ring index and saturating the count at capacity. fn advance(&mut self) -> usize { let cap = self.capacity(); @@ -272,55 +275,24 @@ pub(super) struct BidiagStep { pub(super) beta: f64, } -/// A normal-equation residual `‖Âᵀ(rhs − A x)‖` and the `‖Âᵀ rhs‖` it is judged against. -#[derive(Clone, Copy)] -pub(super) struct NormalEquationResidual { - pub(super) norm: f64, - pub(super) reference: f64, -} - -impl NormalEquationResidual { - /// The drop since `x = 0`; the clamp guards a reference that underflowed to zero. - pub(super) fn relative(&self) -> f64 { - self.norm / self.reference.max(f64::MIN_POSITIVE) - } - - /// An overflowed reference makes the drop vacuous, so keep the stream's own instead. - fn rebase(&mut self, cold: f64) { - if cold > 0.0 && cold.is_finite() { - self.reference = cold; - } - } -} - -/// True residual norms of a candidate solution, recomputed outside the recurrences. -pub(super) struct Certificate { - /// `‖rhs − A x‖`. - pub(super) normr: f64, - /// In the stream's metric (`√(zᵀM⁻¹z)` when preconditioned). - pub(super) normar: NormalEquationResidual, - /// The same outside that metric; `None` when the stream has no metric. - pub(super) normar_raw: Option, -} - -impl Certificate { - /// The pair every [`LsmrResult`](super::LsmrResult) reports: `(‖rhs − A x‖, relative ‖Aᵀr‖)`. - pub(super) fn residuals(&self) -> (f64, f64) { - (self.normr, self.normar.relative()) - } - - /// Take the references from `cold`, an audit of `x = 0` against the warm start's original `b`. - pub(super) fn rebase(&mut self, cold: &Certificate) { - self.normar.rebase(cold.normar.norm); - if let (Some(raw), Some(cold_raw)) = (&mut self.normar_raw, cold.normar_raw) { - raw.rebase(cold_raw.norm); - } +/// `√(gᵀ M⁻¹ g)` for `g = Aᵀ rhs`, scaled so the product survives `‖g‖ ≳ 1e154`. +pub(super) fn metric_gradient_norm( + operator: &A, + preconditioner: &M, + rhs: &[f64], + g: &mut [f64], + mv: &mut [f64], +) -> Result { + operator.apply_adjoint(rhs, g)?; + let plain = par_norm(g); + if !plain.is_finite() { + return Ok(plain); } -} - -/// `|ζ̄₀| = ‖Âᵀ rhs‖`, clamped positive so it can divide a relative normal-equation residual. -pub(super) fn reference_norm(alpha: f64, beta: f64) -> f64 { - (alpha * beta).abs().max(f64::MIN_POSITIVE) + // `M⁻¹` is linear; the raw gradient's metric product overflows past `‖Aᵀr‖ ≈ 1e154`. + let scale = if plain.is_normal() { plain } else { 1.0 }; + scale_in_place(g, 1.0 / scale); + preconditioner.apply(g, mv)?; + Ok(scale * alpha_from_vp(mv, g)?) } /// Stream feeding LSMR `(α, β)` pairs and the matching normalized `v_k`. @@ -329,24 +301,18 @@ pub(super) trait Bidiagonalization { fn step(&mut self) -> Result; /// Most recent normalized basis vector. fn v(&self) -> &[f64]; - /// Clobbers the stream's buffers, so call it only on a terminating path. - fn certify(&mut self, x: &[f64], rhs: &[f64]) -> Result; - /// `‖Aᵀ rhs‖ / ‖rhs‖`, also clobbering. See [`operator_norm_below`]. A stream with no metric - /// has no direction to corroborate, so the default offers no bound and certifies nothing. - fn operator_norm_below(&mut self, _rhs: &[f64]) -> Result { - Ok(0.0) + /// `‖rhs − A x‖`, staging `rhs − A x` for [`restart`](Self::restart); clobbers the stream. + fn residual_norm(&mut self, x: &[f64], rhs: &[f64]) -> Result; + /// Seed a fresh sequence from the staged residual and the `β₁` that + /// [`residual_norm`](Self::residual_norm) returned for it. + fn restart(&mut self, beta: f64) -> Result; + /// After `α₁ = 0`: `Some(‖Aᵀ rhs‖ / ‖rhs‖)` when the metric annihilated a nonzero gradient. + /// A stream with no metric has nowhere to hide one. + fn hidden_gradient(&mut self) -> Result, SolveError> { + Ok(None) } -} - -/// A lower bound on `‖A‖`, for auditing a stop reached before the stream could estimate one. -/// Plain norms throughout: squaring `‖Aᵀ rhs‖` overflows at scales where this does not. -fn operator_norm_below( - operator: &A, - rhs: &[f64], - atr: &mut [f64], -) -> Result { - operator.apply_adjoint(rhs, atr)?; - Ok(super::vec_norm(atr) / super::vec_norm(rhs).max(f64::MIN_POSITIVE)) + /// `‖Aᵀ rhs‖` as a plain norm; clobbers the stream. + fn plain_gradient(&mut self, rhs: &[f64]) -> Result; } impl Bidiagonalization for GolubKahan<'_, A> { @@ -391,18 +357,29 @@ impl Bidiagonalization for GolubKahan<'_, A> { &self.bufs.v } - fn certify(&mut self, x: &[f64], rhs: &[f64]) -> Result { - let normr = residual_into(self.operator, x, rhs, &mut self.bufs.av)?; + fn residual_norm(&mut self, x: &[f64], rhs: &[f64]) -> Result { + residual_into(self.operator, x, rhs, &mut self.bufs.u) + } + + fn restart(&mut self, beta: f64) -> Result { + scale_to_unit(&mut self.bufs.u, beta); self.operator - .apply_adjoint(&self.bufs.av, &mut self.bufs.atu)?; - Ok(Certificate { - normr, - normar: NormalEquationResidual { - norm: super::vec_norm(&self.bufs.atu), - reference: self.normar0, - }, - normar_raw: None, - }) + .apply_adjoint(&self.bufs.u, &mut self.bufs.v)?; + let alpha = finite(par_norm(&self.bufs.v), "α")?; + if alpha > 0.0 { + scale_in_place(&mut self.bufs.v, 1.0 / alpha); + } + if let Some(reorth) = &mut self.bufs.local_reorth { + reorth.clear(); + reorth.push(&self.bufs.v); + } + self.alpha = alpha; + Ok(BidiagStep { alpha, beta }) + } + + fn plain_gradient(&mut self, rhs: &[f64]) -> Result { + self.operator.apply_adjoint(rhs, &mut self.bufs.atu)?; + Ok(par_norm(&self.bufs.atu)) } } @@ -440,31 +417,39 @@ impl Bidiagonalization &self.bufs.v } - fn certify(&mut self, x: &[f64], rhs: &[f64]) -> Result { - let normr = residual_into(self.operator, x, rhs, &mut self.bufs.av)?; + fn residual_norm(&mut self, x: &[f64], rhs: &[f64]) -> Result { + residual_into(self.operator, x, rhs, &mut self.bufs.u) + } + + fn restart(&mut self, beta: f64) -> Result { + scale_to_unit(&mut self.bufs.u, beta); self.operator - .apply_adjoint(&self.bufs.av, &mut self.bufs.atu)?; - let raw = super::vec_norm(&self.bufs.atu); - // `M⁻¹` is linear; the raw gradient's metric product overflows past `‖Aᵀr‖ ≈ 1e154`. - let scale = if raw.is_normal() { raw } else { 1.0 }; - self.bufs.atu.iter_mut().for_each(|g| *g /= scale); + .apply_adjoint(&self.bufs.u, &mut self.bufs.p_tilde)?; self.preconditioner - .apply(&self.bufs.atu, &mut self.bufs.v)?; - Ok(Certificate { - normr, - normar: NormalEquationResidual { - norm: scale * alpha_from_vp(&self.bufs.v, &self.bufs.atu)?, - reference: self.normar0, - }, - normar_raw: Some(NormalEquationResidual { - norm: raw, - reference: self.normar_raw0, - }), - }) + .apply(&self.bufs.p_tilde, &mut self.bufs.v)?; + let alpha = alpha_from_vp(&self.bufs.v, &self.bufs.p_tilde)?; + if alpha > 0.0 { + scale_in_place(&mut self.bufs.v, 1.0 / alpha); + } + if let Some(reorth) = &mut self.bufs.local_reorth { + reorth.clear(); + let inv_alpha = if alpha > 0.0 { 1.0 / alpha } else { 0.0 }; + reorth.push(&self.bufs.v, &self.bufs.p_tilde, inv_alpha); + } + self.alpha = alpha; + self.beta_prev_inv = 1.0; // u was normalized + Ok(BidiagStep { alpha, beta }) } - fn operator_norm_below(&mut self, rhs: &[f64]) -> Result { - operator_norm_below(self.operator, rhs, &mut self.bufs.atu) + fn hidden_gradient(&mut self) -> Result, SolveError> { + // At `α₁ = 0` the stream holds `p̃ = Aᵀu` for the unit `u`, so `‖p̃‖ = ‖Aᵀ rhs‖ / ‖rhs‖`. + let plain = par_norm(&self.bufs.p_tilde); + Ok((plain != 0.0).then_some(plain)) + } + + fn plain_gradient(&mut self, rhs: &[f64]) -> Result { + self.operator.apply_adjoint(rhs, &mut self.bufs.atu)?; + Ok(par_norm(&self.bufs.atu)) } } @@ -500,7 +485,6 @@ pub(super) struct GolubKahan<'a, A: Operator + ?Sized> { bufs: GolubKahanBuffers, /// Last `α` emitted; needed by the next step's u-update. alpha: f64, - normar0: f64, } impl<'a, A: Operator + ?Sized> GolubKahan<'a, A> { @@ -517,28 +501,10 @@ impl<'a, A: Operator + ?Sized> GolubKahan<'a, A> { operator, bufs, alpha: 0.0, - normar0: 0.0, }; let step1 = stream.restart(b_norm)?; - stream.normar0 = reference_norm(step1.alpha, step1.beta); Ok((stream, step1)) } - - /// Seed a fresh sequence from `u` and its norm `beta`. - fn restart(&mut self, beta: f64) -> Result { - scale_to_unit(&mut self.bufs.u, beta); - self.operator - .apply_adjoint(&self.bufs.u, &mut self.bufs.v)?; - let alpha = finite(par_norm(&self.bufs.v), "α")?; - if alpha > 0.0 { - scale_in_place(&mut self.bufs.v, 1.0 / alpha); - } - if let Some(reorth) = &mut self.bufs.local_reorth { - reorth.push(&self.bufs.v); - } - self.alpha = alpha; - Ok(BidiagStep { alpha, beta }) - } } /// Workspaces used by [`ModifiedGolubKahan`]. @@ -579,8 +545,6 @@ pub(super) struct ModifiedGolubKahan<'a, A: Operator + ?Sized, M: Operator + ?Si alpha: f64, /// `1/β_k`; cancels the unnormalization of `u` in the next step. beta_prev_inv: f64, - normar0: f64, - normar_raw0: f64, } impl<'a, A: Operator + ?Sized, M: Operator + ?Sized> ModifiedGolubKahan<'a, A, M> { @@ -601,35 +565,11 @@ impl<'a, A: Operator + ?Sized, M: Operator + ?Sized> ModifiedGolubKahan<'a, A, M bufs, alpha: 0.0, beta_prev_inv: 1.0, - normar0: 0.0, - normar_raw0: 0.0, }; let step1 = stream.restart(b_norm)?; - stream.normar0 = reference_norm(step1.alpha, step1.beta); - stream.normar_raw0 = b_norm * super::vec_norm(&stream.bufs.p_tilde); Ok((stream, step1)) } - /// Seed a fresh sequence from `u` and its norm `beta`. - fn restart(&mut self, beta: f64) -> Result { - scale_to_unit(&mut self.bufs.u, beta); - self.operator - .apply_adjoint(&self.bufs.u, &mut self.bufs.p_tilde)?; - self.preconditioner - .apply(&self.bufs.p_tilde, &mut self.bufs.v)?; - let alpha = alpha_from_vp(&self.bufs.v, &self.bufs.p_tilde)?; - if alpha > 0.0 { - scale_in_place(&mut self.bufs.v, 1.0 / alpha); - } - if let Some(reorth) = &mut self.bufs.local_reorth { - let inv_alpha = if alpha > 0.0 { 1.0 / alpha } else { 0.0 }; - reorth.push(&self.bufs.v, &self.bufs.p_tilde, inv_alpha); - } - self.alpha = alpha; - self.beta_prev_inv = 1.0; // u was normalized - Ok(BidiagStep { alpha, beta }) - } - /// Scaling by `β / α_k` cancels the stored `α_k`; requires `α_k > 0`. fn update_p_tilde(&mut self, beta: f64, beta_inv: f64) -> Result<(), SolveError> { self.operator diff --git a/crates/schwarz-precond/src/lsmr/recurrence.rs b/crates/schwarz-precond/src/lsmr/recurrence.rs index f3c46f9c..e8e23975 100644 --- a/crates/schwarz-precond/src/lsmr/recurrence.rs +++ b/crates/schwarz-precond/src/lsmr/recurrence.rs @@ -5,11 +5,7 @@ //! that yield Algorithm 2.8 of Fong & Saunders, advances the `(x, h, h̄)` //! solution recurrence, and tracks the dual stopping criterion. -use crate::SolveError; - -use super::bidiag::{ - BidiagStep, Certificate, NormalEquationResidual, LSMR_PAR_THRESHOLD, LSMR_UPDATE_CHUNK, -}; +use super::bidiag::{BidiagStep, LSMR_PAR_THRESHOLD, LSMR_UPDATE_CHUNK}; use rayon::iter::{IndexedParallelIterator, ParallelIterator}; use rayon::prelude::{ParallelSlice, ParallelSliceMut}; @@ -76,7 +72,7 @@ impl LsmrRecurrenceState { c_bar: 1.0, s_bar: 0.0, zeta_bar: s1.alpha * s1.beta, - zeta0: super::bidiag::reference_norm(s1.alpha, s1.beta), + zeta0: (s1.alpha * s1.beta).abs().max(f64::MIN_POSITIVE), }, residual: ResidualChain { beta_d: 0.0, @@ -101,7 +97,7 @@ impl LsmrRecurrenceState { self.residual.normr(self.bidiag_qr.beta_dd) } - fn normal_eq_residual_estimate(&self) -> f64 { + pub(super) fn normal_eq_residual_estimate(&self) -> f64 { self.normal_eq_qr.normar() } @@ -316,29 +312,23 @@ impl ConvergenceCriteria { } } - /// The audit left once a stream reports no gradient at all, where the only drop measurable is - /// between the two metrics: the reference within each is the numerator itself, or a clamp. - /// `a_norm_below` replaces the `‖A‖` the stream never got to estimate; being a lower bound it - /// can only overstate the backward error, so it never certifies one the true norm would not. - pub(super) fn corroborated( - &self, - cert: &Certificate, - a_norm_below: impl FnOnce() -> Result, - ) -> Result { - if self.solved(cert) || drops_agree(cert) { - return Ok(true); - } - // Only a metric can hide a direction, and only then is the bound worth an apply. - let Some(raw) = cert.normar_raw else { - return Ok(false); - }; - Ok(backward_error(raw.norm, a_norm_below()?, cert.normr) - <= CERTIFICATION_SLACK * self.rel_tol) + /// The widest `|‖b − A x‖ − ‖r_k‖|` an honest tolerance stop shows: estimate drift is O(ε), + /// a collapsed recurrence misses by orders (cf. van der Vorst & Ye, SISC 22(3), 2000). + fn residual_gap(&self) -> f64 { + CERTIFICATION_SLACK * self.abs_tol } - /// The residual alone meets the tolerance, whatever the normal-equation legs report. - fn solved(&self, cert: &Certificate) -> bool { - cert.normr <= CERTIFICATION_SLACK * self.abs_tol + /// Whether the recomputed residual corroborates the recurrence's own estimate. + pub(super) fn corroborates_residual(&self, recomputed: f64, estimate: f64) -> bool { + (recomputed - estimate).abs() <= self.residual_gap() + } + + /// At `α₁ = 0` the metric reports no gradient. One it annihilated certifies the start only + /// through the residual, or the backward error against a lower bound on `‖A‖`, which can + /// only overstate it. + pub(super) fn corroborates(&self, normr: f64, normar: f64, a_norm_below: f64) -> bool { + normr <= self.residual_gap() + || backward_error(normar, a_norm_below, normr) <= CERTIFICATION_SLACK * self.rel_tol } } @@ -354,36 +344,23 @@ impl ConvergenceState { self.a_norm_sq += s.alpha * s.alpha + s.beta * s.beta; } - /// The stream's own backward error, from its accumulated `‖A‖_F` estimate. - fn ne_ratio(&self, residual: f64, normar: f64) -> f64 { - backward_error(normar, self.a_norm_sq.sqrt(), residual) - } - /// Check both stop criteria against the current scalar state. pub(super) fn check(&self, r: &LsmrRecurrenceState) -> Stop { let residual = r.residual_estimate(); if residual <= self.criteria.abs_tol { return Stop::ResidualTolerance; } - if self.ne_ratio(residual, r.normal_eq_residual_estimate()) <= self.criteria.rel_tol { + // The stream's own backward error, from its accumulated `‖A‖_F` estimate. + let ratio = backward_error( + r.normal_eq_residual_estimate(), + self.a_norm_sq.sqrt(), + residual, + ); + if ratio <= self.criteria.rel_tol { return Stop::NormalEquationTolerance; } Stop::Continue } - - /// True-residual audit of a tolerance stop (cf. van der Vorst & Ye, SISC 22(3), 2000). - pub(super) fn certified(&self, cert: &Certificate) -> bool { - if self.criteria.solved(cert) { - return true; - } - let rel = CERTIFICATION_SLACK * self.criteria.rel_tol; - let dropped = |ne: &NormalEquationResidual| ne.norm <= rel * ne.reference; - // `normr → 0` degenerates the ratio test; the drop of `‖Âᵀr‖` vs its start certifies. - let metric = dropped(&cert.normar) || self.ne_ratio(cert.normr, cert.normar.norm) <= rel; - // A metric that annihilates a direction cannot audit it, so the plain norm must also pass. - // `ne_ratio` cannot audit it: its `‖A‖` is preconditioned, so `M⁻¹`'s scale deflates it. - metric && (cert.normar_raw.is_none_or(|raw| dropped(&raw)) || drops_agree(cert)) - } } /// `‖Aᵀr‖ / (‖A‖‖r‖)`, refusing outright on a denominator carrying no information: clamping one @@ -400,12 +377,5 @@ pub(super) fn backward_error(normar: f64, a_norm: f64, residual: f64) -> f64 { normar / a_norm / residual } -/// The drop outside the stream's metric against the drop inside it, carrying neither the `A` nor -/// the `M` scale; vacuous for a stream with no metric to corroborate. -fn drops_agree(cert: &Certificate) -> bool { - cert.normar_raw - .is_none_or(|raw| raw.relative() <= CERTIFICATION_SLACK * cert.normar.relative()) -} - /// Collapsed recurrences miss by orders of magnitude; the slack absorbs ordinary estimate drift. const CERTIFICATION_SLACK: f64 = 100.0; diff --git a/crates/schwarz-precond/src/lsmr/tests/audit.rs b/crates/schwarz-precond/src/lsmr/tests/audit.rs index 57e34dec..c257e1c1 100644 --- a/crates/schwarz-precond/src/lsmr/tests/audit.rs +++ b/crates/schwarz-precond/src/lsmr/tests/audit.rs @@ -1,20 +1,23 @@ -//! True-residual audit of tolerance stops. +//! The true-residual check of tolerance stops, and the restarts it triggers. use rstest::rstest; -use super::super::bidiag::{BidiagStep, Bidiagonalization, Certificate, NormalEquationResidual}; +use super::super::bidiag::{BidiagStep, Bidiagonalization}; use super::super::recurrence::ConvergenceCriteria; -use super::super::{lsmr_from_bidiag, LsmrStopReason}; +use super::super::{ + lsmr_from_bidiag, EscalationHandler, EscalationPolicy, LsmrResult, LsmrStopReason, + NormalEqReference, +}; +use crate::lsmr::fixtures::FixedIterations; use crate::SolveError; -/// Stream whose first step stops on ResidualTolerance and whose audit is scripted. -struct ScriptedStream { +/// Stops on ResidualTolerance at the first step of every pass; `true_residuals` scripts what +/// the check then sees, one entry per pass. +struct ScriptedStream<'a> { v: Vec, - normr: f64, - normar: f64, - normar_raw: Option, + true_residuals: std::iter::Copied>, } -impl Bidiagonalization for ScriptedStream { +impl Bidiagonalization for ScriptedStream<'_> { fn step(&mut self) -> Result { // β = 0 zeroes φ̄, so the recurrence claims a converged residual immediately. Ok(BidiagStep { @@ -25,100 +28,122 @@ impl Bidiagonalization for ScriptedStream { fn v(&self) -> &[f64] { &self.v } - fn certify(&mut self, _x: &[f64], _rhs: &[f64]) -> Result { - Ok(Certificate { - normr: self.normr, - normar: NormalEquationResidual { - norm: self.normar, - reference: 1.0, - }, - normar_raw: self.normar_raw.map(|norm| NormalEquationResidual { - norm, - reference: 1.0, - }), + fn residual_norm(&mut self, _x: &[f64], _rhs: &[f64]) -> Result { + Ok(self + .true_residuals + .next() + .expect("a scripted residual per pass")) + } + fn restart(&mut self, _beta: f64) -> Result { + Ok(BidiagStep { + alpha: 1.0, + beta: 1.0, }) } + fn plain_gradient(&mut self, _rhs: &[f64]) -> Result { + Ok(1.0) + } } -fn scripted_run(normr: f64, normar: f64, normar_raw: Option) -> super::super::LsmrResult { +/// Counts the handlers a run asks for, one per pass. +#[derive(Default)] +struct CountingPolicy(std::sync::atomic::AtomicUsize); + +impl EscalationPolicy for CountingPolicy { + fn handler(&self) -> Box { + self.0.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + Box::new(FixedIterations(usize::MAX)) + } +} + +/// `tol = 1e-10` against `‖b‖ = 1`, so the gap an honest stop may show is `1e-8`. +fn scripted_run( + true_residuals: &[f64], + maxiter: usize, + escalation: Option<&dyn EscalationPolicy>, +) -> LsmrResult { let stream = ScriptedStream { - v: vec![0.0; 2], - normr, - normar, - normar_raw, + v: vec![1.0, 2.0], + true_residuals: true_residuals.iter().copied(), }; let step1 = BidiagStep { alpha: 1.0, beta: 1.0, }; let criteria = ConvergenceCriteria::new(1.0, 1e-10); - lsmr_from_bidiag(stream, step1, &[1.0, 1.0], None, criteria, 5, None).expect("scripted run") + lsmr_from_bidiag( + stream, + step1, + &[1.0, 1.0], + None, + criteria, + maxiter, + escalation, + ) + .expect("scripted run") } -#[test] -fn collapsed_stop_is_refused_by_the_audit() { - let r = scripted_run(1.0, 1.0, None); - assert!(!r.converged); - assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); - assert_eq!(r.residual_norm, 1.0); - assert_eq!(r.normal_eq_residual, 1.0); -} - -#[test] -fn honest_stop_passes_the_audit() { - let r = scripted_run(1e-12, 1e-12, None); +#[rstest] +#[case::exact(1e-12)] +#[case::inside_the_gap(1e-9)] +fn a_stop_the_true_residual_confirms_is_certified(#[case] true_residual: f64) { + let r = scripted_run(&[true_residual], 5, None); assert!(r.converged); assert_eq!(r.stop_reason, LsmrStopReason::ResidualTolerance); + assert_eq!(r.iterations, 1); + assert_eq!(r.residual_norm, true_residual); } +/// The stream's `v` is constant, so each pass adds the same correction: the restart's pass +/// builds on the refuted iterate rather than starting over. #[test] -fn near_consistent_stop_certifies_via_the_initial_ne_drop() { - // Ratio leg would refuse (1e-12/1e-6 ≫ 100·tol); the drop vs ζ̄₀ = 1 certifies. - let r = scripted_run(1e-6, 1e-12, None); +fn a_refuted_stop_restarts_from_its_iterate() { + let r = scripted_run(&[1.0, 1e-12], 5, None); assert!(r.converged); assert_eq!(r.stop_reason, LsmrStopReason::ResidualTolerance); + assert_eq!(r.iterations, 2); + assert_eq!(r.residual_norm, 1e-12); + assert_eq!(r.x, vec![2.0, 4.0]); } -/// A residual already inside the tolerance certifies on its own: here both gradient legs refuse, -/// the drop being 1e-3 against a reference of 1 and the ratio 1e6. #[test] -fn a_residual_inside_the_tolerance_certifies_without_the_gradient() { - let r = scripted_run(1e-9, 1e-3, None); - assert!(r.converged); - assert_eq!(r.stop_reason, LsmrStopReason::ResidualTolerance); +fn restarts_are_capped_before_the_stop_is_refused() { + let r = scripted_run(&[1.0, 1.0, 1.0], 5, None); + assert!(!r.converged); + assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); + assert_eq!(r.iterations, 3); + assert_eq!(r.residual_norm, 1.0); + assert_eq!(r.x, vec![3.0, 6.0]); } -/// A metric that annihilates part of `Aᵀr` reports it as zero, so the plain norm has to refuse. -#[test] -fn a_stop_the_metric_cannot_see_is_refused() { - let r = scripted_run(1e-6, 1e-12, Some(1e-6)); +/// A refuted stop on the last permitted iteration, or a residual that is not a number at all, +/// has nothing to restart with. +#[rstest] +#[case::budget_exhausted(1.0, 1)] +#[case::non_finite_residual(f64::NAN, 5)] +fn a_refuted_stop_without_a_restart_is_refused(#[case] true_residual: f64, #[case] maxiter: usize) { + let r = scripted_run(&[true_residual], maxiter, None); assert!(!r.converged); assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); + assert_eq!(r.iterations, 1); } -/// A zero or overflowed cold reference cannot replace the stream's own, and never divides. -#[test] -fn a_reference_only_moves_to_a_usable_cold_value() { - let certificate = |norm: f64, reference: f64| Certificate { - normr: 0.0, - normar: NormalEquationResidual { norm, reference }, - normar_raw: None, +/// A zero or overflowed `‖Aᵀb‖` cannot replace the stream's own `ζ̄₀`, and never divides. +/// `step1 = (2, 1)` makes that fallback `2`, so an unusable reference reports `1.0 / 2`. +#[rstest] +#[case::usable(4.0, 0.25)] +#[case::zero(0.0, 0.5)] +#[case::negative(-1.0, 0.5)] +#[case::overflowed(f64::INFINITY, 0.5)] +fn a_report_only_moves_to_a_usable_reference(#[case] metric: f64, #[case] expected: f64) { + let step1 = BidiagStep { + alpha: 2.0, + beta: 1.0, }; - for cold in [0.0, -1.0, f64::INFINITY] { - let mut cert = certificate(1.0, 4.0); - cert.rebase(&certificate(cold, 0.0)); - assert_eq!(cert.normar.reference, 4.0, "cold reference {cold:e}"); - } - let mut cert = certificate(1.0, 4.0); - cert.rebase(&certificate(9.0, 0.0)); - assert_eq!(cert.normar.reference, 9.0); - - assert!(NormalEquationResidual { - norm: 1.0, - reference: 0.0 - } - .relative() - .is_finite()); + assert_eq!( + NormalEqReference::warm(metric, step1).relative(1.0), + expected + ); } /// `‖Aᵀr‖ / (‖A‖‖r‖)` at the ends of the float range, where the obvious spellings certify an @@ -153,3 +178,15 @@ fn the_backward_error_survives_the_ends_of_the_range( ); } } + +/// A restarted pass re-seeds `ζ̄₀`, so the drops it reports are not comparable to the refuted +/// pass's: it takes a fresh handler rather than inheriting a `previous` from before the restart. +#[test] +fn a_restarted_pass_gets_its_own_escalation_handler() { + let policy = CountingPolicy::default(); + let r = scripted_run(&[1.0, 1e-12], 5, Some(&policy)); + + assert!(r.converged); + assert_eq!(r.iterations, 2); + assert_eq!(policy.0.load(std::sync::atomic::Ordering::Relaxed), 2); +} diff --git a/crates/schwarz-precond/src/lsmr/tests/ladder.rs b/crates/schwarz-precond/src/lsmr/tests/ladder.rs index 5843af4e..dab439e5 100644 --- a/crates/schwarz-precond/src/lsmr/tests/ladder.rs +++ b/crates/schwarz-precond/src/lsmr/tests/ladder.rs @@ -122,28 +122,6 @@ fn test_mlsmr_zero_rhs_corrects_non_exact_warm_start() { assert!(result.iterations > 0); } -/// The plain-norm audit must reference `‖Aᵀb‖`, not the warm start's own residual: a far-off `x₀` -/// inflates that residual without limit, and the bound it buys certifies anything. -#[test] -fn a_far_warm_start_does_not_inflate_the_plain_audit_bound() { - let x0 = [1.0 + f64::from(1u32 << 20) * f64::from(1u32 << 20), 0.0]; - let result = mlsmr( - &IdentityOp { n: 2 }, - &[1.0, 1.0], - &DiagOp(vec![1.0, 0.0]), - 1e-10, - 100, - MlsmrOptions { - warm_start: Some(&x0), - ..Default::default() - }, - ) - .expect("warm-started singular-preconditioner solve"); - - assert!(!result.converged); - assert_eq!(result.stop_reason, LsmrStopReason::FalseConvergence); -} - #[rstest] #[case::identity(&IdentityOp { n: 3 }, &[1.0, 2.0, 3.0], &[1.0, 2.0, 3.0])] #[case::null_direction(&ZeroSecondRow, &[0.0, 0.0], &[0.0, 7.0])] diff --git a/crates/schwarz-precond/src/lsmr/tests/solve.rs b/crates/schwarz-precond/src/lsmr/tests/solve.rs index 28ee1f77..9e6c5099 100644 --- a/crates/schwarz-precond/src/lsmr/tests/solve.rs +++ b/crates/schwarz-precond/src/lsmr/tests/solve.rs @@ -175,14 +175,12 @@ fn test_mlsmr_rank_deficient_system() { assert!(normal_equation_residual(&RankDeficientOp, &result.x, &b) < 1e-10); } -/// A preconditioner with a null direction reports `Aᵀr` as zero along it, so the metric audit -/// alone certifies a stop that still carries a full unit of normal-equation residual. Scaling -/// `M⁻¹` must not buy certification either: the plain leg has to be free of that scale. +/// `M⁻¹` must be nonsingular; the true-residual audit cannot see a direction it removed. #[rstest] -#[case::three_dim(&[1.0, 1.0, 0.0])] -#[case::two_dim(&[1.0, 0.0])] -#[case::two_dim_scaled(&[1e20, 0.0])] -fn a_singular_preconditioner_cannot_certify_the_direction_it_annihilates(#[case] m: &[f64]) { +#[case(&[1.0, 0.0])] +#[case(&[1.0, 1.0, 0.0])] +#[case(&[1e20, 0.0])] +fn a_singular_preconditioner_certifies_the_direction_it_annihilates(#[case] m: &[f64]) { let a = IdentityOp { n: m.len() }; let b = vec![1.0; m.len()]; let result = mlsmr( @@ -195,8 +193,7 @@ fn a_singular_preconditioner_cannot_certify_the_direction_it_annihilates(#[case] ) .expect("singular-preconditioner solve"); - assert!(!result.converged, "{:?}", result.stop_reason); - assert_eq!(result.stop_reason, LsmrStopReason::FalseConvergence); + assert!(result.converged, "{:?}", result.stop_reason); let residual = normal_equation_residual(&a, &result.x, &b); assert!( (residual - 1.0).abs() < 1e-9, From 321f2433e21caecc4e878dbf8935c6c6405d623e Mon Sep 17 00:00:00 2001 From: Kristof Schroeder Date: Wed, 23 Sep 2026 16:59:42 +0200 Subject: [PATCH 2/2] fix(lsmr): report a refused stop's normal-equation residual for its iterate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A `FalseConvergence` result carried the recurrence's refuted estimate as `normal_eq_residual`, often a collapsed ~0 next to `converged = false`. The refused path now seeds once from the staged `b − A x` (one Aᵀ and one M⁻¹ apply, off the converged path) and reports its `‖Âᵀr‖`. The `converged` doc and CHANGELOG say what the one-apply check covers: estimate drift, not drift within `range(A)`. Multi-line doc comments added by the exit rework are collapsed to one line. --- CHANGELOG.md | 2 +- crates/schwarz-precond/src/lsmr.rs | 28 +++++------ crates/schwarz-precond/src/lsmr/bidiag.rs | 6 +-- crates/schwarz-precond/src/lsmr/recurrence.rs | 7 +-- .../schwarz-precond/src/lsmr/tests/audit.rs | 47 ++++++++++++++----- 5 files changed, 52 insertions(+), 38 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index df65017b..ddd39803 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,7 +40,7 @@ and this project follows [Semantic Versioning](https://semver.org/). - A design carrying varying slopes on two distinct factors could fail preconditioner construction with `matrix is not symmetric`, when rounding left the two triangles of the exact Schur complement unequal (#229). - A `design` that is neither a 2-D `uint32` array nor a list of `Effect` raised `ValueError` where the documented type is `TypeError`, and `AdditiveSchwarz` accepted a wrong-type `local_solver` at construction, deferring the `TypeError` to solve time (#248). -- LSMR no longer certifies a stop it has not solved. Tolerance stops are audited against the true residual, a non-finite `α`, `β`, `⟨v, Mv⟩`, or `‖b‖` fails with `SolveError::InvalidInput`, and an overflowing or subnormal `‖A‖` no longer zeroes the normal-equation ratio; a failed check reports `LsmrStopReason::FalseConvergence` with `converged = false` (#290, #297, #303, #362). +- LSMR no longer certifies a stop whose recurrence estimates collapsed. Tolerance stops are checked against `‖b − A x‖`, a non-finite `α`, `β`, `⟨v, Mv⟩`, or `‖b‖` fails with `SolveError::InvalidInput`, and an overflowing or subnormal `‖A‖` no longer zeroes the normal-equation ratio; a failed check reports `LsmrStopReason::FalseConvergence` with `converged = false` and the returned iterate's normal-equation residual (#290, #297, #303, #362). - A warm-started solve measures its residuals against the original `b`, including at a budget stop and where `‖Aᵀb‖²` overflows. - An `α` or `β` whose square underflows is recovered from a scaled norm, where LSMR reported `x = 0` converged. - LSMR's solution update no longer drops out when `‖A‖` is small (`≈ 1e-9`) or beyond `1e±154`, where it reported `x = 0` converged. diff --git a/crates/schwarz-precond/src/lsmr.rs b/crates/schwarz-precond/src/lsmr.rs index 5b2453cd..ef43f617 100644 --- a/crates/schwarz-precond/src/lsmr.rs +++ b/crates/schwarz-precond/src/lsmr.rs @@ -43,8 +43,7 @@ pub(crate) fn vec_norm(v: &[f64]) -> f64 { pub struct LsmrResult { /// Solution vector. pub x: Vec, - /// Whether the solver converged within the tolerance. A tolerance stop is held against - /// `‖b − A x‖`, which audits the recurrence and not `M`'s nonsingularity. + /// Whether a tolerance stop matched `‖b − A x‖`, which misses drift within `range(A)`. pub converged: bool, /// Total number of iterations performed. pub iterations: usize, @@ -69,8 +68,7 @@ pub enum LsmrStopReason { NormalEquationTolerance, /// The warm start already solved the system: `b − A x0` was exactly zero. WarmStartExact, - /// A tolerance stop the true residual refuted, as were the restarts from it: the recurrence - /// estimates had collapsed. + /// A tolerance stop, and each restart from it, that `‖b − A x‖` refuted. FalseConvergence, /// The iteration budget was exhausted before convergence. MaxIterations, @@ -247,9 +245,7 @@ pub fn lsmr( lsmr_from_bidiag(bidiag, step1, b, None, criteria, maxiter, None) } -/// Preconditioned LSMR with `M ≈ AᵀA` and one `M⁻¹` application per iteration. -/// -/// `M⁻¹` must be nonsingular; a direction it annihilates can be reported converged unsolved. +/// Preconditioned LSMR with `M ≈ AᵀA`, one `M⁻¹` apply per iteration; `M⁻¹` must be nonsingular. pub fn mlsmr( operator: &A, b: &[f64], @@ -343,8 +339,7 @@ pub fn mlsmr( /// Restarts from a refuted tolerance stop before the solve is refused outright. const MAX_RESTARTS: usize = 2; -/// `‖Aᵀb‖` in the stream's metric: the denominator every reported normal-equation residual -/// divides by, fixed for the solve so a restart reports against the same quantity as pass 1. +/// `‖Aᵀb‖` in the stream's metric, fixed for the solve so every pass reports against it. #[derive(Clone, Copy)] struct NormalEqReference(f64); @@ -368,17 +363,13 @@ impl NormalEqReference { } } -/// A warm start, paired with the `‖Aᵀb‖` it costs the solve: the stream is seeded from -/// `b − A x₀`, so the first `(α₁, β₁)` measures that residual and not `b`. +/// A warm start and its `‖Aᵀb‖`, since the stream's `(α₁, β₁)` measures `b − A x₀` instead. struct WarmStart<'a> { x0: &'a [f64], reference: NormalEqReference, } -/// Runs the LSMR recurrences over a preconditioner-specific bidiagonalization stream. -/// -/// The stream solves for the correction to `x0`; results are in terms of `b`. A tolerance stop -/// is held against the true residual, and a stop it refutes restarts the stream from the iterate. +/// Runs LSMR on the correction to `x0`, restarting from any tolerance stop `b − A x` refutes. fn lsmr_from_bidiag( mut bidiag: B, mut step1: BidiagStep, @@ -493,6 +484,13 @@ fn lsmr_from_bidiag( } if !result.converged { result.stop_reason = LsmrStopReason::FalseConvergence; + // The estimate is the refuted claim; a seed from the staged residual measures `x`. + result.normal_eq_residual = if residual_norm.is_finite() { + let step = bidiag.restart(residual_norm)?; + reference.relative(step.alpha * step.beta) + } else { + residual_norm + }; } result.residual_norm = residual_norm; } diff --git a/crates/schwarz-precond/src/lsmr/bidiag.rs b/crates/schwarz-precond/src/lsmr/bidiag.rs index d5b8bc49..022afea6 100644 --- a/crates/schwarz-precond/src/lsmr/bidiag.rs +++ b/crates/schwarz-precond/src/lsmr/bidiag.rs @@ -303,11 +303,9 @@ pub(super) trait Bidiagonalization { fn v(&self) -> &[f64]; /// `‖rhs − A x‖`, staging `rhs − A x` for [`restart`](Self::restart); clobbers the stream. fn residual_norm(&mut self, x: &[f64], rhs: &[f64]) -> Result; - /// Seed a fresh sequence from the staged residual and the `β₁` that - /// [`residual_norm`](Self::residual_norm) returned for it. + /// Seed a fresh sequence from the staged residual and the `β₁` `residual_norm` returned. fn restart(&mut self, beta: f64) -> Result; - /// After `α₁ = 0`: `Some(‖Aᵀ rhs‖ / ‖rhs‖)` when the metric annihilated a nonzero gradient. - /// A stream with no metric has nowhere to hide one. + /// After `α₁ = 0`: `Some(‖Aᵀ rhs‖ / ‖rhs‖)` when a metric annihilated a nonzero gradient. fn hidden_gradient(&mut self) -> Result, SolveError> { Ok(None) } diff --git a/crates/schwarz-precond/src/lsmr/recurrence.rs b/crates/schwarz-precond/src/lsmr/recurrence.rs index e8e23975..e58c8619 100644 --- a/crates/schwarz-precond/src/lsmr/recurrence.rs +++ b/crates/schwarz-precond/src/lsmr/recurrence.rs @@ -312,8 +312,7 @@ impl ConvergenceCriteria { } } - /// The widest `|‖b − A x‖ − ‖r_k‖|` an honest tolerance stop shows: estimate drift is O(ε), - /// a collapsed recurrence misses by orders (cf. van der Vorst & Ye, SISC 22(3), 2000). + /// Honest estimate drift is O(ε); a collapsed one misses by orders (van der Vorst & Ye 2000). fn residual_gap(&self) -> f64 { CERTIFICATION_SLACK * self.abs_tol } @@ -323,9 +322,7 @@ impl ConvergenceCriteria { (recomputed - estimate).abs() <= self.residual_gap() } - /// At `α₁ = 0` the metric reports no gradient. One it annihilated certifies the start only - /// through the residual, or the backward error against a lower bound on `‖A‖`, which can - /// only overstate it. + /// An annihilated `α₁` certifies via the residual, or via a lower `‖A‖` bound that overstates. pub(super) fn corroborates(&self, normr: f64, normar: f64, a_norm_below: f64) -> bool { normr <= self.residual_gap() || backward_error(normar, a_norm_below, normr) <= CERTIFICATION_SLACK * self.rel_tol diff --git a/crates/schwarz-precond/src/lsmr/tests/audit.rs b/crates/schwarz-precond/src/lsmr/tests/audit.rs index c257e1c1..7723426a 100644 --- a/crates/schwarz-precond/src/lsmr/tests/audit.rs +++ b/crates/schwarz-precond/src/lsmr/tests/audit.rs @@ -4,14 +4,13 @@ use rstest::rstest; use super::super::bidiag::{BidiagStep, Bidiagonalization}; use super::super::recurrence::ConvergenceCriteria; use super::super::{ - lsmr_from_bidiag, EscalationHandler, EscalationPolicy, LsmrResult, LsmrStopReason, - NormalEqReference, + lsmr_from_bidiag, mlsmr, EscalationHandler, EscalationPolicy, LsmrResult, LsmrStopReason, + MlsmrOptions, NormalEqReference, }; -use crate::lsmr::fixtures::FixedIterations; +use crate::lsmr::fixtures::{DenseOp, DiagOp, FixedIterations}; use crate::SolveError; -/// Stops on ResidualTolerance at the first step of every pass; `true_residuals` scripts what -/// the check then sees, one entry per pass. +/// Stops on ResidualTolerance at every pass's first step; `true_residuals` scripts each check. struct ScriptedStream<'a> { v: Vec, true_residuals: std::iter::Copied>, @@ -94,8 +93,7 @@ fn a_stop_the_true_residual_confirms_is_certified(#[case] true_residual: f64) { assert_eq!(r.residual_norm, true_residual); } -/// The stream's `v` is constant, so each pass adds the same correction: the restart's pass -/// builds on the refuted iterate rather than starting over. +/// `v` is constant, so each pass adds the same correction on top of the refuted iterate. #[test] fn a_refuted_stop_restarts_from_its_iterate() { let r = scripted_run(&[1.0, 1e-12], 5, None); @@ -113,11 +111,11 @@ fn restarts_are_capped_before_the_stop_is_refused() { assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); assert_eq!(r.iterations, 3); assert_eq!(r.residual_norm, 1.0); + assert_eq!(r.normal_eq_residual, 1.0); assert_eq!(r.x, vec![3.0, 6.0]); } -/// A refuted stop on the last permitted iteration, or a residual that is not a number at all, -/// has nothing to restart with. +/// A spent budget or a non-numeric residual leaves nothing to restart with. #[rstest] #[case::budget_exhausted(1.0, 1)] #[case::non_finite_residual(f64::NAN, 5)] @@ -126,10 +124,34 @@ fn a_refuted_stop_without_a_restart_is_refused(#[case] true_residual: f64, #[cas assert!(!r.converged); assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); assert_eq!(r.iterations, 1); + // The scripted restart seeds `(1, 1)`: the refuted iterate's `‖Âᵀr‖`, not the claimed 0. + let expected = if true_residual.is_finite() { + 1.0 + } else { + true_residual + }; + assert_eq!(r.normal_eq_residual.to_bits(), expected.to_bits()); +} + +/// `x₀ + Δx` cancels to `x = 0` on the only step allowed, so the stop is refused unrestarted. +#[rstest] +fn a_refused_stop_reports_its_iterates_normal_equation_residual(#[values(1.0, 4.0)] m_inv: f64) { + let op = DenseOp { + rows: 1, + cols: 1, + data: vec![1.0], + }; + let options = MlsmrOptions { + warm_start: Some(&[1e20]), + ..Default::default() + }; + let r = mlsmr(&op, &[1.0], &DiagOp(vec![m_inv]), 1e-12, 1, options).expect("solve"); + assert_eq!(r.stop_reason, LsmrStopReason::FalseConvergence); + assert_eq!(r.x, vec![0.0]); + assert_eq!(r.normal_eq_residual, 1.0); } -/// A zero or overflowed `‖Aᵀb‖` cannot replace the stream's own `ζ̄₀`, and never divides. -/// `step1 = (2, 1)` makes that fallback `2`, so an unusable reference reports `1.0 / 2`. +/// An unusable `‖Aᵀb‖` falls back to `ζ̄₀ = 2` from `step1 = (2, 1)`, reporting `1.0 / 2`. #[rstest] #[case::usable(4.0, 0.25)] #[case::zero(0.0, 0.5)] @@ -179,8 +201,7 @@ fn the_backward_error_survives_the_ends_of_the_range( } } -/// A restarted pass re-seeds `ζ̄₀`, so the drops it reports are not comparable to the refuted -/// pass's: it takes a fresh handler rather than inheriting a `previous` from before the restart. +/// A restart re-seeds `ζ̄₀`, so its drops need a fresh handler, not the refuted pass's `previous`. #[test] fn a_restarted_pass_gets_its_own_escalation_handler() { let policy = CountingPolicy::default();