Implement Curtis Reid Scaling

cpp/include/cuopt/mathematical_optimization/constants.h

@@ -158,6 +158,9 @@
 /* @brief QCQP (barrier) scaling hyper-parameters */
 #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration"
 
+/* @brief PDLP scaling hyper-parameter: Curtis-Reid prescaling toggle */
+#define CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING "pdlp_hyper_enable_curtis_reid_scaling"
+
 /* @brief Barrier initial point safeguard */
 #define CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD "barrier_initial_point_safeguard"
 

cpp/include/cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh

@@ -14,6 +14,8 @@ struct pdlp_hyper_params_t {
   int default_l_inf_ruiz_iterations                 = 10;
   bool do_pock_chambolle_scaling                    = true;
   bool do_ruiz_scaling                              = true;
+  bool do_curtis_reid_scaling                       = true;
+  int number_of_curtis_reid_iterations              = 10;
   double default_alpha_pock_chambolle_rescaling     = 1.0;
   double default_artificial_restart_threshold       = 0.36;
   bool compute_initial_step_size_before_scaling     = false;

cpp/src/grpc/codegen/field_registry.yaml

@@ -790,6 +790,19 @@ pdlp_settings:
       from_proto_cast: "pdlp_precision_t"
       optional: true
 
+    # PDLP scaling hyper-parameters (nested: settings.hyper_params.<field>)
+    # Advanced/internal tuning knob; see pdlp_hyper_params_t in
+    # cpp/include/cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh.
+  - hyper_params:
+    - do_curtis_reid_scaling:
+        description: >-
+          Whether Curtis-Reid prescaling runs before Ruiz/Pock-Chambolle
+          matrix scaling.
+        default: "true"
+        field_num: 35
+        type: bool
+        optional: true
+
 # ─────────────────────────────────────────────────────────────────────────────
 # MIP Solver Settings
 # ─────────────────────────────────────────────────────────────────────────────

cpp/src/grpc/codegen/generated/cuopt_remote_data.proto

@@ -314,6 +314,9 @@ message PDLPSolverSettings {
   // later iterations: -1 automatic (uses the built-in heuristic), 0 disabled,
   // 1 enabled. (default: -1 (automatic))
   optional int32 barrier_adaptive_regularization = 34;
+  // Whether Curtis-Reid prescaling runs before Ruiz/Pock-Chambolle matrix
+  // scaling. (default: true)
+  optional bool do_curtis_reid_scaling = 35;
   PDLPWarmStartData warm_start_data = 50;
 }
 

cpp/src/grpc/codegen/generated/generated_pdlp_settings_to_proto.inc

@@ -40,3 +40,4 @@
   pb_settings->set_save_best_primal_so_far(settings.save_best_primal_so_far);
   pb_settings->set_first_primal_feasible(settings.first_primal_feasible);
   pb_settings->set_pdlp_precision(static_cast<int32_t>(settings.pdlp_precision));
+  pb_settings->set_do_curtis_reid_scaling(settings.hyper_params.do_curtis_reid_scaling);

cpp/src/grpc/codegen/generated/generated_proto_to_pdlp_settings.inc

@@ -90,3 +90,6 @@
   if (pb_settings.has_pdlp_precision()) {
     settings.pdlp_precision = static_cast<pdlp_precision_t>(pb_settings.pdlp_precision());
   }
+  if (pb_settings.has_do_curtis_reid_scaling()) {
+    settings.hyper_params.do_curtis_reid_scaling = pb_settings.do_curtis_reid_scaling();
+  }

cpp/src/math_optimization/solver_settings.cu

@@ -260,6 +260,8 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_")
     {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"},
     {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"},
+    // PDLP scaling hyper-parameter (hidden from default --help: name contains "hyper_")
+    {CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING, &pdlp_settings.hyper_params.do_curtis_reid_scaling, true, "Curtis-Reid prescaling, run before Ruiz/Pock-Chambolle scaling"},
   };
   // String parameters
   string_parameters = {

cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu

@@ -9,6 +9,8 @@
 
 #include <utilities/copy_helpers.hpp>
 
+#include <limits>
+
 #include <cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh>
 #include <cuopt/mathematical_optimization/utilities/segmented_sum_handler.cuh>
 #include <mip_heuristics/mip_constants.hpp>
@@ -138,6 +140,14 @@ void pdlp_initial_scaling_strategy_t<i_t, f_t>::compute_scaling_vectors(
   // master pdlp_solver_t from a shape-0 placeholder)
   if (primal_size_h_ == 0 || dual_size_h_ == 0) return;
 
+  // Curtis-Reid runs first as a prescale: its global log-domain fit removes broad
+  // magnitude skew, leaving Ruiz and Pock-Chambolle to equilibrate locally.
+  //
+  // Skipped under MIP: MIP resets scaling on integer columns. Enabling it needs more
+  // benchmarking.
+  if (hyper_params_.do_curtis_reid_scaling && !running_mip_) {
+    curtis_reid_scaling(hyper_params_.number_of_curtis_reid_iterations);
+  }
   if (hyper_params_.do_ruiz_scaling) { ruiz_inf_scaling(number_of_ruiz_iterations); }
   if (hyper_params_.do_pock_chambolle_scaling) { pock_chambolle_scaling(alpha); }
 }
@@ -383,6 +393,169 @@ __global__ void pock_chambolle_scaling_kernel_col(
   if (threadIdx.x == 0) initial_scaling_view.iteration_variable_scaling[col] = accumulated_value;
 }
 
+template <typename f_t>
+struct a_times_exp_clamped_b {
+  a_times_exp_clamped_b(f_t clamp_bound) : clamp_bound_(clamp_bound) {}
+  HDI f_t operator()(f_t a, f_t b)
+  {
+    f_t clamped = raft::min<f_t>(raft::max<f_t>(b, -clamp_bound_), clamp_bound_);
+    return a * raft::exp(clamped);
+  }
+  f_t clamp_bound_;
+};
+
+// One row of Curtis-Reid's log-domain least-squares fit (see curtis_reid_scaling() below
+// for the full description and references): row_log_scale[row] = mean over row's nonzeros
+// of (-log|a_ij| - col_log_scale[col]). One block per row (like
+// pock_chambolle_scaling_kernel_row), deterministic block-reduce, no atomics. Reads the
+// coefficient as-if-already-scaled by whatever cumulative row/column scale exists so far
+// (matching inf_norm_row_kernel/pock_chambolle_scaling_kernel_row's convention) -- although
+// Curtis-Reid always runs first in the current fixed sequence (cummulative_* still all-1
+// at that point), this keeps the fit correct if that ever changes.
+template <typename i_t, typename f_t, int BLOCK_SIZE>
+__global__ void curtis_reid_row_kernel(const typename mip::problem_t<i_t, f_t>::view_t op_problem,
+                                       const f_t* cummulative_constraint_matrix_scaling,
+                                       const f_t* cummulative_variable_scaling,
+                                       const f_t* col_log_scale,
+                                       f_t* row_log_scale)
+{
+  __shared__ f_t shared[BLOCK_SIZE / raft::WarpSize];
+  auto shared_span      = raft::device_span<f_t>{shared, BLOCK_SIZE / raft::WarpSize};
+  f_t accumulated_value = f_t(0);
+
+  int row        = blockIdx.x;
+  i_t row_offset = op_problem.offsets[row];
+  i_t nnz_in_row = op_problem.offsets[row + 1] - row_offset;
+  f_t row_scale  = cummulative_constraint_matrix_scaling[row];
+
+  for (int j = threadIdx.x; j < nnz_in_row; j += blockDim.x) {
+    i_t col     = op_problem.variables[row_offset + j];
+    f_t abs_val = raft::max<f_t>(raft::abs(op_problem.coefficients[row_offset + j] * row_scale *
+                                           cummulative_variable_scaling[col]),
+                                 std::numeric_limits<f_t>::min());
+    accumulated_value += -raft::log(abs_val) - col_log_scale[col];
+  }
+
+  accumulated_value = deterministic_block_reduce<f_t, BLOCK_SIZE>(shared_span, accumulated_value);
+
+  if (threadIdx.x == 0) {
+    row_log_scale[row] = nnz_in_row > 0 ? accumulated_value / static_cast<f_t>(nnz_in_row) : f_t(0);
+  }
+}
+
+// Column analogue of curtis_reid_row_kernel, over the transposed matrix (mirrors
+// pock_chambolle_scaling_kernel_col).
+template <typename i_t, typename f_t, int BLOCK_SIZE>
+__global__ void curtis_reid_col_kernel(i_t n_variables,
+                                       const f_t* A_T,
+                                       const i_t* A_T_offsets,
+                                       const i_t* A_T_indices,
+                                       const f_t* cummulative_constraint_matrix_scaling,
+                                       const f_t* cummulative_variable_scaling,
+                                       const f_t* row_log_scale,
+                                       f_t* col_log_scale)
+{
+  __shared__ f_t shared[BLOCK_SIZE / raft::WarpSize];
+  auto shared_span      = raft::device_span<f_t>{shared, BLOCK_SIZE / raft::WarpSize};
+  f_t accumulated_value = f_t(0);
+
+  int col        = blockIdx.x;
+  i_t col_offset = A_T_offsets[col];
+  i_t nnz_in_col = A_T_offsets[col + 1] - col_offset;
+  f_t col_scale  = cummulative_variable_scaling[col];
+
+  for (int j = threadIdx.x; j < nnz_in_col; j += blockDim.x) {
+    i_t row     = A_T_indices[col_offset + j];
+    f_t abs_val = raft::max<f_t>(
+      raft::abs(A_T[col_offset + j] * col_scale * cummulative_constraint_matrix_scaling[row]),
+      std::numeric_limits<f_t>::min());
+    accumulated_value += -raft::log(abs_val) - row_log_scale[row];
+  }
+
+  accumulated_value = deterministic_block_reduce<f_t, BLOCK_SIZE>(shared_span, accumulated_value);
+
+  if (threadIdx.x == 0) {
+    col_log_scale[col] = nnz_in_col > 0 ? accumulated_value / static_cast<f_t>(nnz_in_col) : f_t(0);
+  }
+}
+
+// Curtis-Reid prescaling (A. R. Curtis, J. K. Reid, "On the Automatic Scaling of Matrices
+// for Gaussian Elimination", IMA J. Applied Mathematics, 1972; also IIASA Collaborative
+// Paper CP-81-037, https://pure.iiasa.ac.at/id/eprint/1766/7/CP-81-037.pdf): a log-domain
+// least-squares fit run *before* Ruiz/Pock-Chambolle, minimizing
+// sum((log|a_ij| - row_log_scale[i] - col_log_scale[j])^2) via alternating per-row/
+// per-column log-mean fixed-point iteration. This port's sequence and defaults are
+// inspired by the HPR-LP-C codebase (https://github.com/PolyU-IOR/HPR-LP-C).
+template <typename i_t, typename f_t>
+void pdlp_initial_scaling_strategy_t<i_t, f_t>::curtis_reid_scaling(
+  i_t number_of_curtis_reid_iterations)
+{
+#ifdef PDLP_DEBUG_MODE
+  RAFT_CUDA_TRY(cudaDeviceSynchronize());
+  std::cout << "Doing curtis_reid_scaling" << std::endl;
+#endif
+  // Reuse the iteration_* scratch buffers as this phase's row/col log-scale vectors --
+  // same size, same "this phase's working value" role Ruiz/Pock-Chambolle give them.
+  // curtis_reid_row_kernel/curtis_reid_col_kernel read op_problem_scaled_'s coefficients
+  // as-if-already-scaled by the current cummulative_* factors (like Ruiz/Pock-Chambolle's
+  // own kernels do); Curtis-Reid always runs first in compute_scaling_vectors(), so
+  // cummulative_* is still all-1 here in practice.
+  auto& row_log_scale = iteration_constraint_matrix_scaling_;
+  auto& col_log_scale = iteration_variable_scaling_;
+  RAFT_CUDA_TRY(
+    cudaMemsetAsync(row_log_scale.data(), 0, sizeof(f_t) * dual_size_h_, stream_view_.get()));
+  RAFT_CUDA_TRY(
+    cudaMemsetAsync(col_log_scale.data(), 0, sizeof(f_t) * primal_size_h_, stream_view_.get()));
+
+  constexpr i_t number_of_threads = 128;
+  for (i_t iter = 0; iter < number_of_curtis_reid_iterations; ++iter) {
+    curtis_reid_row_kernel<i_t, f_t, number_of_threads>
+      <<<dual_size_h_, number_of_threads, 0, stream_view_.get()>>>(
+        op_problem_scaled_.view(),
+        cummulative_constraint_matrix_scaling_.data(),
+        cummulative_variable_scaling_.data(),
+        col_log_scale.data(),
+        row_log_scale.data());
+    RAFT_CUDA_TRY(cudaPeekAtLastError());
+
+    curtis_reid_col_kernel<i_t, f_t, number_of_threads>
+      <<<primal_size_h_, number_of_threads, 0, stream_view_.get()>>>(
+        primal_size_h_,
+        A_T_.data(),
+        A_T_offsets_.data(),
+        A_T_indices_.data(),
+        cummulative_constraint_matrix_scaling_.data(),
+        cummulative_variable_scaling_.data(),
+        row_log_scale.data(),
+        col_log_scale.data());
+    RAFT_CUDA_TRY(cudaPeekAtLastError());
+  }
+
+  if (running_mip_) { reset_integer_variables(); }
+
+  // Fold the converged log-domain fit into the cumulative scale (exp + clamp, see
+  // a_times_exp_clamped_b): cummulative *= exp(clamp(log_scale)). Unlike Ruiz/
+  // Pock-Chambolle's a_divides_sqrt_b_bounded fold (which incorporates *this iteration's*
+  // norm into a running cumulative), Curtis-Reid already produces the final multiplicative
+  // scale factor directly, so a straight multiply is correct here.
+  //
+  // clamp_bound = 30 (exp(+-30) ~ [9.4e-14, 1.07e13]) bounds the scale-factor range a
+  // pathological log-domain fit could produce.
+  constexpr f_t clamp_bound = f_t(30);
+  raft::linalg::binaryOp(cummulative_constraint_matrix_scaling_.data(),
+                         cummulative_constraint_matrix_scaling_.data(),
+                         row_log_scale.data(),
+                         dual_size_h_,
+                         a_times_exp_clamped_b<f_t>(clamp_bound),
+                         stream_view_.get());
+  raft::linalg::binaryOp(cummulative_variable_scaling_.data(),
+                         cummulative_variable_scaling_.data(),
+                         col_log_scale.data(),
+                         primal_size_h_,
+                         a_times_exp_clamped_b<f_t>(clamp_bound),
+                         stream_view_.get());
+}
+
 template <typename i_t, typename f_t>
 void pdlp_initial_scaling_strategy_t<i_t, f_t>::pock_chambolle_scaling(f_t alpha)
 {
@@ -1080,7 +1253,24 @@ pdlp_initial_scaling_strategy_t<i_t, f_t>::view()
     const typename pdlp_initial_scaling_strategy_t<int, F_TYPE>::view_t initial_scaling_view, \
     F_TYPE* A_T,                                                                              \
     int* A_T_offsets,                                                                         \
-    int* A_T_indices);
+    int* A_T_indices);                                                                        \
+                                                                                              \
+  template __global__ void curtis_reid_row_kernel<int, F_TYPE, 128>(                          \
+    const typename mip::problem_t<int, F_TYPE>::view_t op_problem,                            \
+    const F_TYPE* cummulative_constraint_matrix_scaling,                                      \
+    const F_TYPE* cummulative_variable_scaling,                                               \
+    const F_TYPE* col_log_scale,                                                              \
+    F_TYPE* row_log_scale);                                                                   \
+                                                                                              \
+  template __global__ void curtis_reid_col_kernel<int, F_TYPE, 128>(                          \
+    int n_variables,                                                                          \
+    const F_TYPE* A_T,                                                                        \
+    const int* A_T_offsets,                                                                   \
+    const int* A_T_indices,                                                                   \
+    const F_TYPE* cummulative_constraint_matrix_scaling,                                      \
+    const F_TYPE* cummulative_variable_scaling,                                               \
+    const F_TYPE* row_log_scale,                                                              \
+    F_TYPE* col_log_scale);
 
 #if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT
 INSTANTIATE(float)

cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cuh

@@ -129,6 +129,10 @@ class pdlp_initial_scaling_strategy_t {
   void ruiz_iter_local();
   // Shard-local end-to-end Pock-Chambolle pass. Exposed for distributed PDLP:
   void pock_chambolle_scaling(f_t alpha);
+  // Curtis-Reid prescaling pass -- see the implementation in initial_scaling.cu for
+  // details and references. Not exposed to distributed PDLP yet (no cross-shard-coherent
+  // version written).
+  void curtis_reid_scaling(i_t number_of_curtis_reid_iterations);
   // Iteration_* scratch buffers used by ruiz_iter_local /
   // pock_chambolle_scaling. Exposed mutably so distributed PDLP can grow
   // them back to full size after the ctor's release (see distributed_scaling).

cpp/tests/linear_programming/grpc/grpc_client_test.cpp

@@ -2277,6 +2277,8 @@ TEST(MapperRoundtrip, PDLPSettingsAllFields)
   orig.pdlp_precision               = pdlp_precision_t::MixedPrecision;
   orig.save_best_primal_so_far      = true;
   orig.first_primal_feasible        = true;
+  orig.hyper_params.do_curtis_reid_scaling =
+    false;  // not the default true, to detect overwrite-on-decode
 
   cuopt::remote::PDLPSolverSettings pb;
   map_pdlp_settings_to_proto(orig, &pb);
@@ -2319,6 +2321,7 @@ TEST(MapperRoundtrip, PDLPSettingsAllFields)
   EXPECT_EQ(restored.pdlp_precision, pdlp_precision_t::MixedPrecision);
   EXPECT_EQ(restored.save_best_primal_so_far, true);
   EXPECT_EQ(restored.first_primal_feasible, true);
+  EXPECT_EQ(restored.hyper_params.do_curtis_reid_scaling, false);
 }
 
 TEST(MapperRoundtrip, PDLPSettingsIterationLimitSentinel)
@@ -2418,6 +2421,28 @@ TEST(MapperRoundtrip, PDLPSettingsBarrierIterativeRefinementExplicitFalseRoundtr
   EXPECT_FALSE(restored.barrier_iterative_refinement);
 }
 
+TEST(MapperRoundtrip, PDLPSettingsCurtisReidScalingOmittedPreservesDefault)
+{
+  cuopt::remote::PDLPSolverSettings pb;
+
+  pdlp_solver_settings_t<int32_t, double> fresh;
+  ASSERT_TRUE(fresh.hyper_params.do_curtis_reid_scaling);
+  map_proto_to_pdlp_settings(pb, fresh);
+  EXPECT_TRUE(fresh.hyper_params.do_curtis_reid_scaling)
+    << "Omitted optional bool must preserve the C++ default `true`";
+}
+
+TEST(MapperRoundtrip, PDLPSettingsCurtisReidScalingExplicitFalseRoundtrips)
+{
+  cuopt::remote::PDLPSolverSettings pb;
+  pb.set_do_curtis_reid_scaling(false);
+  ASSERT_TRUE(pb.has_do_curtis_reid_scaling());
+
+  pdlp_solver_settings_t<int32_t, double> restored;
+  map_proto_to_pdlp_settings(pb, restored);
+  EXPECT_FALSE(restored.hyper_params.do_curtis_reid_scaling);
+}
+
 // Wide-coverage sanity: a default-constructed proto (no fields touched on the
 // wire) must, after the mapper, leave every C++ scalar settings field at its
 // in-class default. Spot-checks a representative cross-section of the fields
@@ -2457,6 +2482,7 @@ TEST(MapperRoundtrip, PDLPSettingsDefaultProtoPreservesAllCppDefaults)
   EXPECT_EQ(after.dual_postsolve, fresh.dual_postsolve);
   EXPECT_EQ(after.eliminate_dense_columns, fresh.eliminate_dense_columns);
   EXPECT_EQ(after.barrier_iterative_refinement, fresh.barrier_iterative_refinement);
+  EXPECT_EQ(after.hyper_params.do_curtis_reid_scaling, fresh.hyper_params.do_curtis_reid_scaling);
   // Numeric defaults != 0.
   EXPECT_EQ(after.num_gpus, fresh.num_gpus);
   EXPECT_EQ(after.folding, fresh.folding);

cpp/tests/linear_programming/pdlp_distributed_test.cu

@@ -47,6 +47,8 @@ static void expect_distributed_matches_base(raft::handle_t const& handle,
 
   pdlp_solver_settings_t<int, double> base_settings{};
   base_settings.method = method_t::PDLP;
+  // Curtis-Reid scaling is not supported yet for multi-GPU.
+  base_settings.hyper_params.do_curtis_reid_scaling = false;
 
   auto base_op = mps_data_model_to_optimization_problem<int, double>(&handle, problem);
   auto base    = solve_lp(base_op, base_settings);