wip

Author: Yanyu000

highs/pdlp/hipdlp/pdhg.cc

@@ -645,7 +645,6 @@ if (DEBUG_MODE){
     }
 
     CUDA_CHECK(cudaGraphLaunch(graphExec, gpu_stream_));
-    CUDA_CHECK(cudaStreamSynchronize(gpu_stream_));
 #else
     for (int i = 2; i <= PDHG_CHECK_INTERVAL - 1; i++) {
       performHalpernPdhgStep(false, i);
@@ -1362,7 +1361,8 @@ double PDLPSolver::computePrimalFeasibility(
     }
   }
 
-  return linalg::vectorNorm(primal_residual);
+  return linalg::vectorNorm(primal_residual) /
+         scaling_.getConstraintBoundScaling();
 }
 
 void PDLPSolver::computeDualSlacks(const std::vector<double>& dualResidual,
@@ -1448,7 +1448,8 @@ double PDLPSolver::computeDualFeasibility(const std::vector<double>& ATy_vector,
     }
   }
 
-  double dual_feasibility = linalg::vectorNorm(dual_residual);
+  double dual_feasibility = linalg::vectorNorm(dual_residual) /
+                            scaling_.getObjectiveVectorScaling();
 
   return dual_feasibility;
 }
@@ -1493,7 +1494,7 @@ PDLPSolver::computeDualityGap(const std::vector<double>& x,
 double PDLPSolver::computeDualObjective(const std::vector<double>& y,
                                         const std::vector<double>& dSlackPos,
                                         const std::vector<double>& dSlackNeg) {
-  double dual_obj = lp_.offset_;
+  double dual_obj = 0.0;
 
   // Compute b'y (or rhs'y in cuPDLP notation)
   for (HighsInt i = 0; i < lp_.num_row_; ++i) {
@@ -1514,7 +1515,10 @@ double PDLPSolver::computeDualObjective(const std::vector<double>& y,
     }
   }
 
-  return dual_obj;
+  return dual_obj /
+             (scaling_.getConstraintBoundScaling() *
+              scaling_.getObjectiveVectorScaling()) +
+         lp_.offset_;
 }
 
 bool PDLPSolver::checkConvergence(
@@ -1537,27 +1541,31 @@ bool PDLPSolver::checkConvergence(
   results.dual_feasibility = dual_feasibility;
 
   // Compute objectives
-  double primal_obj = lp_.offset_;
+  double primal_obj = 0.0;
   for (HighsInt i = 0; i < lp_.num_col_; ++i) {
     primal_obj += lp_.col_cost_[i] * x[i];
   }
-  results.primal_obj = primal_obj;
+  primal_obj /=
+      scaling_.getConstraintBoundScaling() *
+      scaling_.getObjectiveVectorScaling();
+  results.primal_obj = primal_obj + lp_.offset_;
 
   // Pass the now-populated slack vectors to computeDualObjective
   double dual_obj = computeDualObjective(y, dSlackPos, dSlackNeg);
   results.dual_obj = dual_obj;
 
   // Compute duality gap
-  double duality_gap = primal_obj - dual_obj;
+  double duality_gap = results.primal_obj - dual_obj;
   results.duality_gap = std::abs(duality_gap);
 
   // Compute relative gap (matching cuPDLP formula)
   double relative_obj_gap =
-      std::abs(duality_gap) / (1.0 + std::abs(primal_obj) + std::abs(dual_obj));
+      std::abs(duality_gap) /
+      (1.0 + std::abs(results.primal_obj) + std::abs(dual_obj));
   results.relative_obj_gap = relative_obj_gap;
 
 #if PDLP_DEBUG_LOG
-  debugPdlpFeasOptLog(debug_pdlp_log_file_, iter, primal_obj, dual_obj,
+  debugPdlpFeasOptLog(debug_pdlp_log_file_, iter, results.primal_obj, dual_obj,
                       relative_obj_gap,
                       primal_feasibility / (1.0 + unscaled_rhs_norm_),
                       dual_feasibility / (1.0 + unscaled_c_norm_), type);
@@ -1936,10 +1944,11 @@ void PDLPSolver::unscaleSolution(std::vector<double>& x,
   y_current_ = y;
 
   const std::vector<double>& col_scale = scaling_.getColScaling();
+  const double objective_scale = scaling_.getObjectiveVectorScaling();
   if (!dSlackPos_.empty() && col_scale.size() == dSlackPos_.size()) {
     for (size_t i = 0; i < dSlackPos_.size(); ++i) {
-      dSlackPos_[i] *= col_scale[i];
-      dSlackNeg_[i] *= col_scale[i];
+      dSlackPos_[i] *= col_scale[i] / objective_scale;
+      dSlackNeg_[i] *= col_scale[i] / objective_scale;
     }
   }
 }
@@ -1994,7 +2003,7 @@ void PDLPSolver::initializeStepSizes() {
   best_primal_weight_ = primal_weight_;
   stepsize_.beta = primal_weight_ * primal_weight_;
   params_.omega = primal_weight_;
-  
+
   if (params_.step_size_strategy != StepSizeStrategy::FIXED &&
       params_.step_size_strategy != StepSizeStrategy::PID) {
     logger_.info(
@@ -2828,13 +2837,16 @@ bool PDLPSolver::checkConvergenceGpu(const HighsInt iter, const double* d_x,
                              cudaMemcpyDeviceToHost, gpu_stream_));
   CUDA_CHECK(cudaStreamSynchronize(gpu_stream_));
 
+  const double constraint_scale = scaling_.getConstraintBoundScaling();
+  const double objective_scale = scaling_.getObjectiveVectorScaling();
+  const double combined_scale = constraint_scale * objective_scale;
   double primal_feas_sq = h_results[0];
   double dual_feas_sq = h_results[1];
-  double primal_obj = h_results[2] + lp_.offset_;
-  double dual_obj = h_results[3] + lp_.offset_;
+  double primal_obj = h_results[2] / combined_scale + lp_.offset_;
+  double dual_obj = h_results[3] / combined_scale + lp_.offset_;
 
-  results.primal_feasibility = std::sqrt(primal_feas_sq);
-  results.dual_feasibility = std::sqrt(dual_feas_sq);
+  results.primal_feasibility = std::sqrt(primal_feas_sq) / constraint_scale;
+  results.dual_feasibility = std::sqrt(dual_feas_sq) / objective_scale;
   results.primal_obj = primal_obj;
   results.dual_obj = dual_obj;
 

highs/pdlp/hipdlp/scaling.cc

@@ -22,6 +22,8 @@ void Scaling::initialize(const HighsLp& lp) {
   col_scale_.assign(lp.num_col_, 1.0);
   row_scale_.assign(lp.num_row_, 1.0);
   is_scaled_ = false;
+  constraint_bound_scale_ = 1.0;
+  objective_vector_scale_ = 1.0;
 
   // use linalg to compute norms
   norm_cost_ = linalg::computeCostNorm(lp, 2.0);
@@ -51,6 +53,10 @@ void Scaling::scaleProblem() {
     applyL2Scaling();
     is_scaled_ = true;
   }
+
+  highsLogDev(params_->log_options_, HighsLogType::kInfo,
+              "Applying bound-objective scaling...\n");
+  applyBoundObjectiveScaling();
 }
 
 void Scaling::applyRuizScaling() {
@@ -228,6 +234,45 @@ void Scaling::applyL2Scaling() {
   }
 }
 
+void Scaling::applyBoundObjectiveScaling() {
+  double rhs_norm_sq = 0.0;
+  for (HighsInt i = 0; i < lp_->num_row_; ++i) {
+    const double lower = lp_->row_lower_[i];
+    const double upper = lp_->row_upper_[i];
+    if (std::isfinite(lower) && lower != upper) rhs_norm_sq += lower * lower;
+    if (std::isfinite(upper)) rhs_norm_sq += upper * upper;
+  }
+
+  const double rhs_norm = std::sqrt(rhs_norm_sq);
+  const double obj_norm = computeNorm(lp_->col_cost_.data(), lp_->num_col_, 2.0);
+
+  constraint_bound_scale_ = 1.0 / (1.0 + rhs_norm);
+  objective_vector_scale_ = 1.0 / (1.0 + obj_norm);
+
+  for (HighsInt i = 0; i < lp_->num_row_; ++i) {
+    if (lp_->row_lower_[i] > -kHighsInf) {
+      lp_->row_lower_[i] *= constraint_bound_scale_;
+    }
+    if (lp_->row_upper_[i] < kHighsInf) {
+      lp_->row_upper_[i] *= constraint_bound_scale_;
+    }
+  }
+
+  for (HighsInt i = 0; i < lp_->num_col_; ++i) {
+    lp_->col_cost_[i] *= objective_vector_scale_;
+    if (lp_->col_lower_[i] > -kHighsInf) {
+      lp_->col_lower_[i] *= constraint_bound_scale_;
+    }
+    if (lp_->col_upper_[i] < kHighsInf) {
+      lp_->col_upper_[i] *= constraint_bound_scale_;
+    }
+  }
+
+  if (constraint_bound_scale_ != 1.0 || objective_vector_scale_ != 1.0) {
+    is_scaled_ = true;
+  }
+}
+
 void Scaling::applyScaling(const std::vector<double>& col_scaling,
                            const std::vector<double>& row_scaling) {
   // Scale cost vector: c_scaled = c / col_scaling
@@ -273,12 +318,12 @@ void Scaling::unscaleSolution(std::vector<double>& x,
 
   // Unscale primal variables: x_original = x_scaled / col_scale
   for (size_t i = 0; i < x.size(); ++i) {
-    x[i] /= col_scale_[i];
+    x[i] /= (col_scale_[i] * constraint_bound_scale_);
   }
 
   // Unscale dual variables: y_original = y_scaled / row_scale
   for (size_t i = 0; i < y.size(); ++i) {
-    y[i] /= row_scale_[i];
+    y[i] /= (row_scale_[i] * objective_vector_scale_);
   }
 }
 

highs/pdlp/hipdlp/scaling.hpp

@@ -41,6 +41,8 @@ class Scaling {
   bool isScaled() const { return is_scaled_; }
   const std::vector<double>& getColScaling() const { return col_scale_; }
   const std::vector<double>& getRowScaling() const { return row_scale_; }
+  double getConstraintBoundScaling() const { return constraint_bound_scale_; }
+  double getObjectiveVectorScaling() const { return objective_vector_scale_; }
 
   double getNormCost() const { return norm_cost_; }
   double getNormRhs() const { return norm_rhs_; }
@@ -51,6 +53,8 @@ class Scaling {
   std::vector<double> col_scale_;
   std::vector<double> row_scale_;
   bool is_scaled_ = false;
+  double constraint_bound_scale_ = 1.0;
+  double objective_vector_scale_ = 1.0;
 
   double norm_cost_;
   double norm_rhs_;
@@ -59,6 +63,7 @@ class Scaling {
   void applyRuizScaling();
   void applyPockChambolleScaling();
   void applyL2Scaling();
+  void applyBoundObjectiveScaling();
 
   // Helper function to apply scaling factors to the problem
   void applyScaling(const std::vector<double>& col_scaling,