diff --git a/src/solver/core/solver.rs b/src/solver/core/solver.rs index 6af492c5..93ed73e8 100644 --- a/src/solver/core/solver.rs +++ b/src/solver/core/solver.rs @@ -133,6 +133,7 @@ where pub step_lhs: V, pub step_rhs: V, pub prev_vars: V, + pub best_vars: V, pub info: I, pub solution: SO, pub(crate) settings: SE, // not public to avoid unchecked modifications @@ -304,6 +305,11 @@ where self.info.print_status(&self.settings).unwrap(); }} + // remember the best iterate seen so far, so that a later numerical + // failure can fall back to it rather than a degraded final iterate + self.info + .save_best_iterate(&self.variables, &mut self.best_vars, &self.settings); + // termination checks // -------------- @@ -447,6 +453,22 @@ where } timeit! {timers => "post-process"; { + // on numerical failure the final iterate is often severely degraded + // (e.g. a blown-up primal residual): fall back to the best iterate + // seen so that the "almost" convergence check below applies to it. + // Iteration and time limits are included because `post_process` + // runs that same check for them, and their final iterate is no more + // likely to be the best one seen than a failed solve's is. + if self.info.get_status().is_errored() + || matches!( + self.info.get_status(), + SolverStatus::MaxIterations | SolverStatus::MaxTime + ) + { + self.info + .reset_to_best_iterate(&mut self.variables, &self.best_vars, &self.settings); + } + //check for "almost" convergence case and then extract solution self.info.post_process(&self.residuals, &self.settings); self.solution diff --git a/src/solver/core/traits.rs b/src/solver/core/traits.rs index cc2453ed..da2f406a 100644 --- a/src/solver/core/traits.rs +++ b/src/solver/core/traits.rs @@ -212,6 +212,21 @@ where /// restore a prior iterate fn reset_to_prev_iterate(&mut self, variables: &mut Self::V, prev_variables: &Self::V); + /// save the current iterate as the best seen so far if it improves on it + fn save_best_iterate( + &mut self, + variables: &Self::V, + best_variables: &mut Self::V, + settings: &Self::SE, + ); + /// restore the best iterate seen if it improves on the current one + fn reset_to_best_iterate( + &mut self, + variables: &mut Self::V, + best_variables: &Self::V, + settings: &Self::SE, + ); + /// Record some of the top level solver's choice of various /// scalars. `μ = ` normalized gap. `α = ` computed step length. /// `σ = ` multiplier for the updated centering parameter. diff --git a/src/solver/implementations/default/info.rs b/src/solver/implementations/default/info.rs index 0f777a54..f221a348 100644 --- a/src/solver/implementations/default/info.rs +++ b/src/solver/implementations/default/info.rs @@ -51,6 +51,18 @@ pub struct DefaultInfo { pub(crate) prev_gap_abs: T, /// relative duality gap from previous iteration pub(crate) prev_gap_rel: T, + + // best iterate seen so far (by worst-case residual/gap merit), used as a + // fallback when the solver terminates with a numerical error on a + // degraded final iterate + pub(crate) best_cost_primal: T, + pub(crate) best_cost_dual: T, + pub(crate) best_res_primal: T, + pub(crate) best_res_dual: T, + pub(crate) best_gap_abs: T, + pub(crate) best_gap_rel: T, + pub(crate) best_ktratio: T, + pub(crate) best_merit: T, /// solve time pub solve_time: f64, /// solver status @@ -88,6 +100,7 @@ where self.status = SolverStatus::Unsolved; self.iterations = 0; self.solve_time = 0f64; + self.best_merit = T::infinity(); timers.reset_timer("solve"); } @@ -252,6 +265,67 @@ where variables.copy_from(prev_variables); } + fn save_best_iterate( + &mut self, + variables: &Self::V, + best_variables: &mut Self::V, + settings: &DefaultSettings, + ) { + // Iterates on an infeasibility path (κ/τ > 1) are never candidates. + if self.ktratio > T::one() { + return; + } + + // Strict improvement, tested so that a non-finite merit is rejected. Every + // comparison against NaN is false, so `merit >= best_merit` would treat a NaN + // iterate as an improvement, overwrite a good one with it, and then disable the + // fallback altogether, since the restore below requires a finite best_merit -- + // and a solve whose iterates have gone non-finite is exactly the case this is + // meant to rescue. The same test rejects an infinite merit, which is what a + // zero reduced tolerance would produce. + let merit = self.termination_merit(settings); + if !(merit < self.best_merit) { + return; + } + + self.best_merit = merit; + self.best_cost_primal = self.cost_primal; + self.best_cost_dual = self.cost_dual; + self.best_res_primal = self.res_primal; + self.best_res_dual = self.res_dual; + self.best_gap_abs = self.gap_abs; + self.best_gap_rel = self.gap_rel; + self.best_ktratio = self.ktratio; + + best_variables.copy_from(variables); + } + + fn reset_to_best_iterate( + &mut self, + variables: &mut Self::V, + best_variables: &Self::V, + settings: &DefaultSettings, + ) { + if !self.best_merit.is_finite() { + return; + } + + let merit = self.termination_merit(settings); + if self.ktratio <= T::one() && merit <= self.best_merit { + return; + } + + self.cost_primal = self.best_cost_primal; + self.cost_dual = self.best_cost_dual; + self.res_primal = self.best_res_primal; + self.res_dual = self.best_res_dual; + self.gap_abs = self.best_gap_abs; + self.gap_rel = self.best_gap_rel; + self.ktratio = self.best_ktratio; + + variables.copy_from(best_variables); + } + fn save_scalars(&mut self, μ: T, α: T, σ: T, iter: u32) { self.mu = μ; self.step_length = α; @@ -362,6 +436,32 @@ where } } + // How far the current iterate is from satisfying `is_solved` at the *reduced* + // tolerances: each quantity divided by the tolerance that will judge it, combined + // exactly as in that test, so the duality gap enters through whichever of its + // absolute/relative forms is closer to passing. A merit below 1 passes. + // + // The reduced tolerances are the right yardstick because this merit only ever ranks + // candidates for an unsuccessful exit, where `check_convergence_almost` -- which uses + // exactly these tolerances -- decides whether the restored iterate can still be + // reported as AlmostSolved. Ranking by them makes the selection safe: the chosen + // iterate has the smallest reduced merit of all candidates, so if any candidate would + // have passed that check, the chosen one passes it too. Comparing the raw quantities + // instead can prefer an iterate that fails the check over one that passes, because the + // tolerances differ from each other -- by default `reduced_tol_feas` is 1e-4 while + // `reduced_tol_gap_rel` is 5e-5. + fn termination_merit(&self, settings: &DefaultSettings) -> T { + let gap = T::min( + self.gap_abs / settings.reduced_tol_gap_abs, + self.gap_rel / settings.reduced_tol_gap_rel, + ); + let feas = T::max( + self.res_primal / settings.reduced_tol_feas, + self.res_dual / settings.reduced_tol_feas, + ); + T::max(gap, feas) + } + fn is_solved(&self, tol_gap_abs: T, tol_gap_rel: T, tol_feas: T) -> bool { ((self.gap_abs < tol_gap_abs) || (self.gap_rel < tol_gap_rel)) && (self.res_primal < tol_feas) @@ -388,3 +488,175 @@ where && (self.res_dual_inf < -tol_infeas_rel * residuals.dot_qx) } } + +#[cfg(test)] +mod tests { + use super::*; + use crate::solver::core::traits::Info; + + // fabricate an iterate whose merit-relevant fields are set directly, with x[0] + // tagging which iterate it is so we can see which one survives + fn set_iterate( + info: &mut DefaultInfo, + vars: &mut DefaultVariables, + res: f64, + ktratio: f64, + tag: f64, + ) { + info.res_primal = res; + info.res_dual = res; + info.gap_rel = res; + info.gap_abs = res; + info.ktratio = ktratio; + vars.x[0] = tag; + } + + #[test] + fn best_iterate_survives_a_degraded_final_iterate() { + let settings = DefaultSettings::::default(); + let mut info = DefaultInfo::::new(); + let mut vars = DefaultVariables::::new(1, 1); + let mut best = DefaultVariables::::new(1, 1); + info.reset(&mut Default::default()); + + // a good iterate, then a better one, then the blowup that ends the solve + set_iterate(&mut info, &mut vars, 1e-6, 0.5, 1.0); + info.save_best_iterate(&vars, &mut best, &settings); + set_iterate(&mut info, &mut vars, 1e-10, 0.5, 2.0); + info.save_best_iterate(&vars, &mut best, &settings); + set_iterate(&mut info, &mut vars, 3e-5, 0.5, 3.0); + info.save_best_iterate(&vars, &mut best, &settings); + + info.reset_to_best_iterate(&mut vars, &best, &settings); + + assert_eq!(vars.x[0], 2.0, "the 1e-10 iterate should be restored"); + assert_eq!(info.res_primal, 1e-10); + assert_eq!(info.gap_rel, 1e-10); + } + + #[test] + fn a_better_final_iterate_is_left_alone() { + let settings = DefaultSettings::::default(); + let mut info = DefaultInfo::::new(); + let mut vars = DefaultVariables::::new(1, 1); + let mut best = DefaultVariables::::new(1, 1); + info.reset(&mut Default::default()); + + set_iterate(&mut info, &mut vars, 1e-6, 0.5, 1.0); + info.save_best_iterate(&vars, &mut best, &settings); + set_iterate(&mut info, &mut vars, 1e-9, 0.5, 2.0); + + info.reset_to_best_iterate(&mut vars, &best, &settings); + + assert_eq!( + vars.x[0], 2.0, + "the current iterate is the best one; keep it" + ); + assert_eq!(info.res_primal, 1e-9); + } + + // The acceptance check applied on an unsuccessful exit uses a different tolerance for + // each quantity -- by default reduced_tol_feas is 1e-4 but reduced_tol_gap_rel is 5e-5 -- + // so ranking candidates by the raw worst quantity can prefer one that fails the check + // over one that passes it. Here the first candidate has the smaller raw maximum (8e-5 + // against 9e-5) yet its gap is outside tolerance, while the second is inside on every + // count; the second must be the one restored. + #[test] + fn merit_ranks_by_distance_from_the_acceptance_check() { + let settings = DefaultSettings::::default(); + let mut info = DefaultInfo::::new(); + let mut vars = DefaultVariables::::new(1, 1); + let mut best = DefaultVariables::::new(1, 1); + info.reset(&mut Default::default()); + + // candidate 1: gap outside the reduced tolerance, residuals well inside + info.gap_abs = 8e-5; + info.gap_rel = 8e-5; + info.res_primal = 1e-5; + info.res_dual = 1e-5; + info.ktratio = 0.5; + vars.x[0] = 1.0; + info.save_best_iterate(&vars, &mut best, &settings); + assert!(!info.is_solved( + settings.reduced_tol_gap_abs, + settings.reduced_tol_gap_rel, + settings.reduced_tol_feas + )); + + // candidate 2: larger raw maximum, but inside tolerance on every count + info.gap_abs = 2e-5; + info.gap_rel = 2e-5; + info.res_primal = 9e-5; + info.res_dual = 9e-5; + vars.x[0] = 2.0; + info.save_best_iterate(&vars, &mut best, &settings); + assert!(info.is_solved( + settings.reduced_tol_gap_abs, + settings.reduced_tol_gap_rel, + settings.reduced_tol_feas + )); + assert_eq!( + best.x[0], 2.0, + "the iterate that passes the check must be preferred" + ); + + // and it is the one a degraded final iterate is replaced by + info.gap_abs = 1e-2; + info.gap_rel = 1e-2; + info.res_primal = 1e-2; + info.res_dual = 1e-2; + vars.x[0] = 3.0; + info.reset_to_best_iterate(&mut vars, &best, &settings); + assert_eq!(vars.x[0], 2.0); + assert_eq!(info.res_primal, 9e-5); + } + + // A solve whose iterates go non-finite is the case this fallback exists for, so a NaN + // iterate must not be able to displace the good one that preceded it. Comparisons + // against NaN are all false, which makes this easy to get wrong. + #[test] + fn a_nan_iterate_cannot_displace_the_best_one() { + let settings = DefaultSettings::::default(); + let mut info = DefaultInfo::::new(); + let mut vars = DefaultVariables::::new(1, 1); + let mut best = DefaultVariables::::new(1, 1); + info.reset(&mut Default::default()); + + set_iterate(&mut info, &mut vars, 1e-10, 0.5, 1.0); + info.save_best_iterate(&vars, &mut best, &settings); + + // the solve blows up: residuals, gap and κ/τ are all NaN + set_iterate(&mut info, &mut vars, f64::NAN, f64::NAN, 2.0); + info.save_best_iterate(&vars, &mut best, &settings); + assert_eq!(best.x[0], 1.0, "the NaN iterate must not be saved"); + assert!( + info.best_merit.is_finite(), + "the fallback must remain armed" + ); + + info.reset_to_best_iterate(&mut vars, &best, &settings); + assert_eq!(vars.x[0], 1.0, "the good iterate must be restored"); + assert_eq!(info.res_primal, 1e-10); + } + + #[test] + fn iterates_on_the_infeasibility_path_are_not_candidates() { + let settings = DefaultSettings::::default(); + let mut info = DefaultInfo::::new(); + let mut vars = DefaultVariables::::new(1, 1); + let mut best = DefaultVariables::::new(1, 1); + info.reset(&mut Default::default()); + + // tiny residuals but ktratio > 1: this is an infeasibility certificate path + set_iterate(&mut info, &mut vars, 1e-12, 5.0, 1.0); + info.save_best_iterate(&vars, &mut best, &settings); + set_iterate(&mut info, &mut vars, 3e-5, 5.0, 2.0); + info.reset_to_best_iterate(&mut vars, &best, &settings); + + assert_eq!( + vars.x[0], 2.0, + "nothing was ever saved, so nothing is restored" + ); + assert!(!info.best_merit.is_finite()); + } +} diff --git a/src/solver/implementations/default/solver.rs b/src/solver/implementations/default/solver.rs index 515141e7..ab2ce9e8 100644 --- a/src/solver/implementations/default/solver.rs +++ b/src/solver/implementations/default/solver.rs @@ -105,12 +105,13 @@ where let step_rhs = DefaultVariables::::new(data.n,data.m); let step_lhs = DefaultVariables::::new(data.n,data.m); let prev_vars = DefaultVariables::::new(data.n,data.m); + let best_vars = DefaultVariables::::new(data.n,data.m); // configure empty user callbacks output = Self{ data,variables,residuals,kktsystem, - step_lhs,step_rhs,prev_vars,info, + step_lhs,step_rhs,prev_vars,best_vars,info, solution,cones,settings, timers: None, callbacks: SolverCallbacks::default(),