fix restart

Author: Yanyu000

highs/pdlp/hipdlp/pdhg.cc

@@ -26,6 +26,15 @@ static constexpr double kDivergentMovement = 1e10;
 
 using namespace std;
 
+void vecPrint(const std::vector<double>& vec, const char* name) {
+  std::cout << name << ": [";
+  for (size_t i = 0; i < vec.size(); ++i) {
+    std::cout << vec[i];
+    if (i < vec.size() - 1) std::cout << ", ";
+  }
+  std::cout << "]" << std::endl;
+}
+
 void PDLPSolver::preprocessLp() {
   logger_.info(
       "Preprocessing LP using cupdlp formulation (slack variables for "
@@ -412,9 +421,9 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
 
     // --- 3. Convergence and Restart Check (BEFORE iterate update) ---
     bool bool_checking = (iter < 10) ||
-                         (iter == (params_.max_iterations - 1)) ||
-                         (iter > 0 && (iter % PDHG_CHECK_INTERVAL == 0));
-
+                         (iter == (params_.max_iterations - 1));
+    
+    bool_checking = (bool_checking || iter % PDHG_CHECK_INTERVAL == 0);
     if (bool_checking) {
       ComputeAverageIterate(Ax_avg, ATy_avg);
       // Reset the average iterate accumulation
@@ -467,36 +476,35 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
           restart_scheme_.Check(iter, current_results, average_results);
 
       if (restart_info.should_restart) {
-        std::vector<double> restart_x, restart_y;
         if (restart_info.restart_to_average) {
-          // Restart from the average iterate
-          restart_x = x_avg_;
-          restart_y = y_avg_;
+          restart_scheme_.primal_feas_last_restart_ = average_results.primal_feasibility;
+          restart_scheme_.dual_feas_last_restart_ = average_results.dual_feasibility;
+          restart_scheme_.duality_gap_last_restart_ = average_results.duality_gap;
+
+          x_current_ = x_avg_;
+          y_current_ = y_avg_;
 
           Ax_cache_ = Ax_avg;
           ATy_cache_ = ATy_avg;
         } else {
-          // Restart from the current iterate
-          restart_x = x_current_;
-          restart_y = y_current_;
+          restart_scheme_.primal_feas_last_restart_ = current_results.primal_feasibility;
+          restart_scheme_.dual_feas_last_restart_ = current_results.dual_feasibility;
+          restart_scheme_.duality_gap_last_restart_ = current_results.duality_gap;
         }
 
         // Perform the primal weight update using z^{n,0} and z^{n-1,0}
         PDHG_Compute_Step_Size_Ratio(working_params);
         current_eta_ = working_params.eta;
         restart_scheme_.passParams(&working_params);
-        restart_scheme_.UpdateBeta(working_params.omega * working_params.omega);
 
-        x_at_last_restart_ = restart_x;  // Current becomes the new last
-        y_at_last_restart_ = restart_y;
+        x_at_last_restart_ = x_current_;  // Current becomes the new last
+        y_at_last_restart_ = y_current_;
 
         std::fill(x_sum_.begin(), x_sum_.end(), 0.0);
         std::fill(y_sum_.begin(), y_sum_.end(), 0.0);
         sum_weights_ = 0.0;
 
-        // Set x_current_ and y_current_ for the next iteration's starting point
-        x_current_ = restart_x;
-        y_current_ = restart_y;
+        restart_scheme_.last_restart_iter_ = iter;
 
         // Recompute Ax and ATy for the restarted iterates
         linalg::Ax(lp, x_current_, Ax_cache_);
@@ -607,15 +615,15 @@ void PDLPSolver::PDHG_Compute_Step_Size_Ratio(
   // iterates.
   double primal_diff_norm = linalg::diffTwoNorm(x_at_last_restart_, x_current_);
   double dual_diff_norm = linalg::diffTwoNorm(y_at_last_restart_, y_current_);
-
+ 
   double dMeanStepSize = std::sqrt(stepsize_.primal_step *
                                        stepsize_.dual_step);
 
 
   // 2. Update the primal weight (beta = omega^2) if movements are significant.
   if (std::min(primal_diff_norm, dual_diff_norm) > 1e-10) {
     double beta_update_ratio = dual_diff_norm / primal_diff_norm;
-    double old_beta = working_params.omega * working_params.omega;
+    double old_beta = stepsize_.beta;
 
     double dLogBetaUpdate = 0.5 * std::log(beta_update_ratio) +
                             0.5 * std::log(std::sqrt(old_beta));
@@ -626,6 +634,7 @@ void PDLPSolver::PDHG_Compute_Step_Size_Ratio(
   stepsize_.dual_step = stepsize_.primal_step * stepsize_.beta;
   working_params.eta = std::sqrt(stepsize_.primal_step * stepsize_.dual_step);
   working_params.omega = std::sqrt(stepsize_.beta);
+  restart_scheme_.UpdateBeta(stepsize_.beta);
 }
 
 void PDLPSolver::UpdateAverageIterates(const std::vector<double>& x,
@@ -1365,3 +1374,4 @@ double PDLPSolver::ComputeNonlinearity(const std::vector<double>& delta_primal,
   }
   return std::abs(nonlinearity);
 }
+

highs/pdlp/hipdlp/restart.cc

@@ -156,7 +156,7 @@ RestartInfo RestartScheme::Check(int current_iter,
   debugPdlpRestartLog(debug_pdlp_log_file_, current_iter, mu_current,
                           mu_average);
   double candidate_score;
-  bool restart_to_average;
+  bool restart_to_average = false;
   if (mu_current < mu_average) {
     candidate_score = mu_current;
     restart_to_average = false;
@@ -168,11 +168,12 @@ RestartInfo RestartScheme::Check(int current_iter,
   // 2. Check the three restart conditions in order
   bool should_restart = false;
   PDHG_restart_choice restart_choice_for_logic = PDHG_NO_RESTART;
-
+ 
   if ((current_iter - last_restart_iter_) >= 0.36 * current_iter) {
     // Condition 1: Artificial Restart
     should_restart = true;
   } else {
+
     double mu_last_restart = compute_score(beta_, primal_feas_last_restart_,
                                            dual_feas_last_restart_, duality_gap_last_restart_);
 
@@ -201,9 +202,13 @@ RestartInfo RestartScheme::Check(int current_iter,
     duality_gap_last_candidate_ = current_results.duality_gap;
   }
 
-  // NOTE: Your main solver loop is responsible for calling `SetLastRestartIter(current_iter)`
-  // when it receives `should_restart = true`. This function does not modify that state directly,
-  // it only determines if it *should* happen.
+  if (should_restart ) {
+    if (restart_to_average){
+      std::cout << "Last restart was iter " << last_restart_iter_ << ": average\n";
+    } else {
+      std::cout << "Last restart was iter " << last_restart_iter_ << ": current\n";
+    }
+  }
 
   return RestartInfo(should_restart, restart_to_average);
 }

highs/pdlp/hipdlp/restart.hpp

@@ -42,6 +42,7 @@ class RestartScheme {
                     const SolverResults& average_results);
 
   int GetLastRestartIter() const { return last_restart_iter_; }
+  int last_restart_iter_ = 0;
   void passParams(PrimalDualParams* params) { params_ = params; };
   void passLogOptions(const HighsLogOptions* log_options) {
     log_options_ = log_options;
@@ -59,6 +60,14 @@ class RestartScheme {
   };
   void SetLastRestartIter(int iter) { last_restart_iter_ = iter; };
 
+    // State for adaptive restart
+  //dPrimalFeasLastRestart = primal_feas_last_restart_
+  //dDualFeasLastRestart = dual_feas_last_restart_;
+  //dDualityGapLastRestart = duality_gap_last_restart_;
+  double primal_feas_last_restart_ = 0.0;
+  double dual_feas_last_restart_ = 0.0;
+  double duality_gap_last_restart_ = 0.0;
+
  private:
   PrimalDualParams* params_;
   const HighsLogOptions* log_options_;
@@ -67,13 +76,7 @@ class RestartScheme {
 
   RestartStrategy strategy_ = RestartStrategy::NO_RESTART;
   int fixed_restart_interval_ = 100;
-  int last_restart_iter_ = 0;
   double beta_;
-
-  // State for adaptive restart
-  double primal_feas_last_restart_ = 0.0;
-  double dual_feas_last_restart_ = 0.0;
-  double duality_gap_last_restart_ = 0.0;
 
   double primal_feas_last_candidate_ = 0.0;
   double dual_feas_last_candidate_ = 0.0;