solver: return the best iterate seen when terminating unsuccessfully

src/solver/core/solver.rs

@@ -132,7 +132,7 @@ where
     pub cones: C,
     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
@@ -426,8 +426,9 @@ where
                 StrategyCheckpoint::Fail => {α = T::zero(); break}
             }
 
-            // Copy previous iterate in case the next one is a dud
-            self.info.save_prev_iterate(&self.variables,&mut self.prev_vars);
+            // Record progress scalars, and checkpoint this iterate if it
+            // is the best seen so far (restored should the solve fail)
+            self.info.checkpoint_iterate(&self.variables,&mut self.best_vars,&self.settings);
 
             self.variables.add_step(&self.step_lhs, α);
 
@@ -447,6 +448,18 @@ where
         }
 
         timeit! {timers => "post-process"; {
+            // if the solver failed, report the best iterate seen rather
+            // than the point at which it gave up (unless the final iterate
+            // is trending towards an infeasibility certificate).   The
+            // "almost" convergence check below is then evaluated against
+            // the restored iterate.
+            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
@@ -592,10 +605,13 @@ mod internal {
                 // there is no problem, so nothing to do
                 output = StrategyCheckpoint::NoUpdate;
             } else {
-                // recover old iterate since "insufficient progress" often
-                // involves actual degradation of results
-                self.info
-                    .reset_to_prev_iterate(&mut self.variables, &self.prev_vars);
+                // recover the best iterate since "insufficient progress"
+                // often involves actual degradation of results
+                self.info.reset_to_best_iterate(
+                    &mut self.variables,
+                    &self.best_vars,
+                    &self.settings,
+                );
 
                 // If problem is asymmetric, we can try to continue with the dual-only strategy
                 if !self.cones.is_symmetric() && (scaling == ScalingStrategy::PrimalDual) {

src/solver/core/traits.rs

@@ -207,10 +207,26 @@ where
     /// Return `true` if termination conditions have been reached.
     fn check_termination(&mut self, residuals: &Self::R, settings: &Self::SE, iter: u32) -> bool;
 
-    /// save a prior iterate
-    fn save_prev_iterate(&mut self, variables: &Self::V, prev_variables: &mut Self::V);
-    /// restore a prior iterate
-    fn reset_to_prev_iterate(&mut self, variables: &mut Self::V, prev_variables: &Self::V);
+    /// Record the current iterate's scalars for the next iteration's
+    /// progress checks, and checkpoint the iterate itself (into
+    /// `best_variables`) if it is the best seen so far by the
+    /// termination-criteria merit.
+    fn checkpoint_iterate(
+        &mut self,
+        variables: &Self::V,
+        best_variables: &mut Self::V,
+        settings: &Self::SE,
+    );
+    /// Restore the best checkpointed iterate, if one exists, is better
+    /// than the iterate currently held, and the current iterate is not
+    /// trending towards an infeasibility certificate.  Returns true if
+    /// the restore was performed.
+    fn reset_to_best_iterate(
+        &mut self,
+        variables: &mut Self::V,
+        best_variables: &Self::V,
+        settings: &Self::SE,
+    ) -> bool;
 
     /// Record some of the top level solver's choice of various
     /// scalars. `μ = ` normalized gap.  `α = ` computed step length.

src/solver/implementations/default/info.rs

@@ -51,6 +51,27 @@ 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 termination-criteria merit.
+    // Restored on solver failure so that the reported solution is
+    // the best one encountered rather than the point at which the
+    // solver gave up.   `best_merit` is ∞ until a checkpoint is taken.
+    /// merit of the best iterate (∞ if none yet checkpointed)
+    pub(crate) best_merit: T,
+    /// primal objective value at the best iterate
+    pub(crate) best_cost_primal: T,
+    /// dual objective value at the best iterate
+    pub(crate) best_cost_dual: T,
+    /// primal residual at the best iterate
+    pub(crate) best_res_primal: T,
+    /// dual residual at the best iterate
+    pub(crate) best_res_dual: T,
+    /// absolute duality gap at the best iterate
+    pub(crate) best_gap_abs: T,
+    /// relative duality gap at the best iterate
+    pub(crate) best_gap_rel: T,
+    /// κ/τ ratio at the best iterate
+    pub(crate) best_ktratio: T,
     /// solve time
     pub solve_time: f64,
     /// solver status
@@ -88,6 +109,7 @@ where
         self.status = SolverStatus::Unsolved;
         self.iterations = 0;
         self.solve_time = 0f64;
+        self.best_merit = T::infinity();
 
         timers.reset_timer("solve");
     }
@@ -230,26 +252,76 @@ where
         self.status != SolverStatus::Unsolved
     }
 
-    fn save_prev_iterate(&mut self, variables: &Self::V, prev_variables: &mut Self::V) {
+    fn checkpoint_iterate(
+        &mut self,
+        variables: &Self::V,
+        best_variables: &mut Self::V,
+        settings: &DefaultSettings<T>,
+    ) {
+        // scalars from the previous iteration, used by the
+        // poor-progress tests in check_termination
         self.prev_cost_primal = self.cost_primal;
         self.prev_cost_dual = self.cost_dual;
         self.prev_res_primal = self.res_primal;
         self.prev_res_dual = self.res_dual;
         self.prev_gap_abs = self.gap_abs;
         self.prev_gap_rel = self.gap_rel;
 
-        prev_variables.copy_from(variables);
+        // Additionally checkpoint the iterate itself if it is the best
+        // seen so far, so that a failed solve can still report the best
+        // point encountered rather than the one at which it gave up.
+        // Only iterates on the optimality branch (κ/τ ≤ 1) are candidates:
+        // for κ/τ > 1 the iterate is trending towards an infeasibility
+        // certificate and its cost/gap values are not meaningful.
+        if self.ktratio <= T::one() {
+            let merit = self.termination_merit(settings);
+            if merit.is_finite() && merit < self.best_merit {
+                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_prev_iterate(&mut self, variables: &mut Self::V, prev_variables: &Self::V) {
-        self.cost_primal = self.prev_cost_primal;
-        self.cost_dual = self.prev_cost_dual;
-        self.res_primal = self.prev_res_primal;
-        self.res_dual = self.prev_res_dual;
-        self.gap_abs = self.prev_gap_abs;
-        self.gap_rel = self.prev_gap_rel;
+    fn reset_to_best_iterate(
+        &mut self,
+        variables: &mut Self::V,
+        best_variables: &Self::V,
+        settings: &DefaultSettings<T>,
+    ) -> bool {
+        // nothing to restore if no iterate was ever checkpointed.   If the
+        // current iterate has κ/τ > 1 it is trending towards an
+        // infeasibility certificate, which a restore would mask, so keep it.
+        if !self.best_merit.is_finite() || self.ktratio > T::one() {
+            return false;
+        }
+
+        // Never trade the current iterate for a worse one.   The iteration
+        // loop can exit *before* checkpointing -- an undersized step or a
+        // failed KKT solve breaks out of the loop body ahead of the
+        // checkpoint -- so the iterate held at exit is frequently better
+        // than anything stored, and restoring unconditionally would discard
+        // the best answer instead of recovering it.
+        if self.termination_merit(settings) <= self.best_merit {
+            return false;
+        }
+
+        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(prev_variables);
+        variables.copy_from(best_variables);
+        true
     }
 
     fn save_scalars(&mut self, μ: T, α: T, σ: T, iter: u32) {
@@ -368,6 +440,34 @@ where
             && (self.res_dual < tol_feas)
     }
 
+    // Distance of the current iterate from satisfying the `is_solved`
+    // test at the *reduced* tolerances: each termination quantity
+    // normalized by its tolerance, combined exactly as in that test (the
+    // duality gap counts via whichever of its absolute/relative forms is
+    // closer to passing).  An iterate with merit < 1 satisfies the test.
+    //
+    // The reduced tolerances are the right yardstick because this merit
+    // only ever selects among iterates 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 gives a guarantee: the selected
+    // iterate has the smallest reduced merit among all candidates, so if
+    // any candidate would have passed that check, the selected one passes
+    // it too.  Normalizing by the full tolerances instead can rank an
+    // iterate that fails the reduced check above one that passes it,
+    // turning a reduced-accuracy success into a reported failure.
+    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_primal_infeasible(
         &self,
         residuals: &DefaultResiduals<T>,
@@ -388,3 +488,121 @@ where
             && (self.res_dual_inf < -tol_infeas_rel * residuals.dot_qx)
     }
 }
+
+// ---------------------------------------------------------------------------
+// best-iterate checkpointing tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::solver::core::traits::Info;
+
+    // gives an info with the stated termination quantities and κ/τ,
+    // with best_merit initialized as at the start of a solve
+    fn info_with(res: f64, gap: f64, ktratio: f64) -> DefaultInfo<f64> {
+        let mut info = DefaultInfo::<f64>::default();
+        info.best_merit = f64::INFINITY;
+        info.res_primal = res;
+        info.res_dual = res;
+        info.gap_abs = gap;
+        info.gap_rel = gap;
+        info.ktratio = ktratio;
+        info
+    }
+
+    fn set_current(info: &mut DefaultInfo<f64>, res: f64, gap: f64, ktratio: f64) {
+        info.res_primal = res;
+        info.res_dual = res;
+        info.gap_abs = gap;
+        info.gap_rel = gap;
+        info.ktratio = ktratio;
+    }
+
+    #[test]
+    fn test_best_iterate_checkpoint_and_restore() {
+        let settings = DefaultSettings::<f64>::default();
+        let mut info = info_with(1e-3, 1e-3, 0.5);
+        let mut best = DefaultVariables::<f64>::new(2, 1);
+        let mut vars = DefaultVariables::<f64>::new(2, 1);
+
+        // first candidate is checkpointed
+        vars.x[0] = 1.0;
+        info.cost_primal = 10.0;
+        info.checkpoint_iterate(&vars, &mut best, &settings);
+        assert_eq!(best.x[0], 1.0);
+
+        // a worse iterate is not
+        set_current(&mut info, 1e-2, 1e-2, 0.5);
+        vars.x[0] = 2.0;
+        info.checkpoint_iterate(&vars, &mut best, &settings);
+        assert_eq!(best.x[0], 1.0);
+
+        // a better one is
+        set_current(&mut info, 1e-6, 1e-6, 0.5);
+        vars.x[0] = 3.0;
+        info.cost_primal = 30.0;
+        info.checkpoint_iterate(&vars, &mut best, &settings);
+        assert_eq!(best.x[0], 3.0);
+
+        // an infeasibility-trending iterate is never a candidate,
+        // however good its (meaningless) residuals look
+        set_current(&mut info, 1e-9, 1e-9, 2.0);
+        vars.x[0] = 4.0;
+        info.checkpoint_iterate(&vars, &mut best, &settings);
+        assert_eq!(best.x[0], 3.0);
+
+        // restore declines while the current iterate trends infeasible...
+        assert!(!info.reset_to_best_iterate(&mut vars, &best, &settings));
+        assert_eq!(vars.x[0], 4.0);
+
+        // ...and otherwise restores the checkpointed variables and scalars
+        set_current(&mut info, 1e-1, 1e-1, 0.5);
+        info.cost_primal = 99.0;
+        assert!(info.reset_to_best_iterate(&mut vars, &best, &settings));
+        assert_eq!(vars.x[0], 3.0);
+        assert_eq!(info.cost_primal, 30.0);
+        assert_eq!(info.res_primal, 1e-6);
+    }
+
+    // The iteration loop can break before checkpointing (undersized step,
+    // failed KKT solve), so the iterate held at exit is often better than
+    // any checkpoint.   Restoring then would discard the best answer, so a
+    // restore must only happen when it is an improvement.
+    #[test]
+    fn test_best_iterate_never_restores_something_worse() {
+        let settings = DefaultSettings::<f64>::default();
+        let mut info = info_with(1e-3, 1e-3, 0.5);
+        let mut best = DefaultVariables::<f64>::new(2, 1);
+        let mut vars = DefaultVariables::<f64>::new(2, 1);
+
+        // checkpoint a mediocre iterate
+        vars.x[0] = 1.0;
+        info.checkpoint_iterate(&vars, &mut best, &settings);
+        assert_eq!(best.x[0], 1.0);
+
+        // now hold a *better* iterate that was never checkpointed
+        set_current(&mut info, 1e-9, 1e-9, 0.5);
+        vars.x[0] = 2.0;
+        assert!(!info.reset_to_best_iterate(&mut vars, &best, &settings));
+        assert_eq!(vars.x[0], 2.0); // kept, not overwritten
+
+        // but a genuinely worse current iterate is replaced
+        set_current(&mut info, 1e-1, 1e-1, 0.5);
+        vars.x[0] = 3.0;
+        assert!(info.reset_to_best_iterate(&mut vars, &best, &settings));
+        assert_eq!(vars.x[0], 1.0);
+    }
+
+    #[test]
+    fn test_best_iterate_no_checkpoint_no_restore() {
+        let settings = DefaultSettings::<f64>::default();
+        let mut info = info_with(1e-3, 1e-3, 0.5);
+        info.best_merit = f64::INFINITY; // nothing checkpointed
+        let best = DefaultVariables::<f64>::new(2, 1);
+        let mut vars = DefaultVariables::<f64>::new(2, 1);
+        vars.x[0] = 7.0;
+        assert!(!info.reset_to_best_iterate(&mut vars, &best, &settings));
+        assert_eq!(vars.x[0], 7.0);
+    }
+}

src/solver/implementations/default/solver.rs

@@ -104,13 +104,13 @@ where
         // work variables for assembling step direction LHS/RHS
         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,best_vars,info,
             solution,cones,settings,
             timers: None,
             callbacks: SolverCallbacks::default(),