Fall back to the best iterate when terminating with a numerical error

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

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.

src/solver/implementations/default/info.rs

@@ -51,6 +51,18 @@ pub struct DefaultInfo<T> {
     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<T>,
+    ) {
+        // 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<T>,
+    ) {
+        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>) -> 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<f64>,
+        vars: &mut DefaultVariables<f64>,
+        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::<f64>::default();
+        let mut info = DefaultInfo::<f64>::new();
+        let mut vars = DefaultVariables::<f64>::new(1, 1);
+        let mut best = DefaultVariables::<f64>::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::<f64>::default();
+        let mut info = DefaultInfo::<f64>::new();
+        let mut vars = DefaultVariables::<f64>::new(1, 1);
+        let mut best = DefaultVariables::<f64>::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::<f64>::default();
+        let mut info = DefaultInfo::<f64>::new();
+        let mut vars = DefaultVariables::<f64>::new(1, 1);
+        let mut best = DefaultVariables::<f64>::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::<f64>::default();
+        let mut info = DefaultInfo::<f64>::new();
+        let mut vars = DefaultVariables::<f64>::new(1, 1);
+        let mut best = DefaultVariables::<f64>::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::<f64>::default();
+        let mut info = DefaultInfo::<f64>::new();
+        let mut vars = DefaultVariables::<f64>::new(1, 1);
+        let mut best = DefaultVariables::<f64>::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());
+    }
+}

src/solver/implementations/default/solver.rs

@@ -105,12 +105,13 @@ where
         let step_rhs  = DefaultVariables::<T>::new(data.n,data.m);
         let step_lhs  = DefaultVariables::<T>::new(data.n,data.m);
         let prev_vars = DefaultVariables::<T>::new(data.n,data.m);
+        let best_vars = DefaultVariables::<T>::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(),