diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 3656791a98..2dabf07d3c 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/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" diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh b/cpp/include/cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh index 6aee5213a7..46ebcc815e 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh +++ b/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; diff --git a/cpp/src/grpc/codegen/field_registry.yaml b/cpp/src/grpc/codegen/field_registry.yaml index 718b79a052..3b98703f34 100644 --- a/cpp/src/grpc/codegen/field_registry.yaml +++ b/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.) + # 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 # ───────────────────────────────────────────────────────────────────────────── diff --git a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto b/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto index cf00d8c760..a143877a49 100644 --- a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto +++ b/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; } diff --git a/cpp/src/grpc/codegen/generated/generated_pdlp_settings_to_proto.inc b/cpp/src/grpc/codegen/generated/generated_pdlp_settings_to_proto.inc index 07d3d59bbe..ca9f136ac2 100644 --- a/cpp/src/grpc/codegen/generated/generated_pdlp_settings_to_proto.inc +++ b/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(settings.pdlp_precision)); + pb_settings->set_do_curtis_reid_scaling(settings.hyper_params.do_curtis_reid_scaling); diff --git a/cpp/src/grpc/codegen/generated/generated_proto_to_pdlp_settings.inc b/cpp/src/grpc/codegen/generated/generated_proto_to_pdlp_settings.inc index 4f749cb301..f34e1cf94f 100644 --- a/cpp/src/grpc/codegen/generated/generated_proto_to_pdlp_settings.inc +++ b/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(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(); + } diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 5a4ab72c32..6aba0e8d0a 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -260,6 +260,8 @@ solver_settings_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 = { diff --git a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu b/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu index 28f8b47257..68e485694f 100644 --- a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu +++ b/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu @@ -9,6 +9,8 @@ #include +#include + #include #include #include @@ -138,6 +140,14 @@ void pdlp_initial_scaling_strategy_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 +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(raft::max(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 +__global__ void curtis_reid_row_kernel(const typename mip::problem_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{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(raft::abs(op_problem.coefficients[row_offset + j] * row_scale * + cummulative_variable_scaling[col]), + std::numeric_limits::min()); + accumulated_value += -raft::log(abs_val) - col_log_scale[col]; + } + + accumulated_value = deterministic_block_reduce(shared_span, accumulated_value); + + if (threadIdx.x == 0) { + row_log_scale[row] = nnz_in_row > 0 ? accumulated_value / static_cast(nnz_in_row) : f_t(0); + } +} + +// Column analogue of curtis_reid_row_kernel, over the transposed matrix (mirrors +// pock_chambolle_scaling_kernel_col). +template +__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{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( + raft::abs(A_T[col_offset + j] * col_scale * cummulative_constraint_matrix_scaling[row]), + std::numeric_limits::min()); + accumulated_value += -raft::log(abs_val) - row_log_scale[row]; + } + + accumulated_value = deterministic_block_reduce(shared_span, accumulated_value); + + if (threadIdx.x == 0) { + col_log_scale[col] = nnz_in_col > 0 ? accumulated_value / static_cast(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 +void pdlp_initial_scaling_strategy_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 + <<>>( + 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 + <<>>( + 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(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(clamp_bound), + stream_view_.get()); +} + template void pdlp_initial_scaling_strategy_t::pock_chambolle_scaling(f_t alpha) { @@ -1080,7 +1253,24 @@ pdlp_initial_scaling_strategy_t::view() const typename pdlp_initial_scaling_strategy_t::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( \ + const typename mip::problem_t::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 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) diff --git a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cuh b/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cuh index 96d7f0629c..1f8e24ad26 100644 --- a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cuh +++ b/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). diff --git a/cpp/tests/linear_programming/grpc/grpc_client_test.cpp b/cpp/tests/linear_programming/grpc/grpc_client_test.cpp index 506c7fa62b..9a71a3b0aa 100644 --- a/cpp/tests/linear_programming/grpc/grpc_client_test.cpp +++ b/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 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 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); diff --git a/cpp/tests/linear_programming/pdlp_distributed_test.cu b/cpp/tests/linear_programming/pdlp_distributed_test.cu index d76f775c08..aa7f05fd2d 100644 --- a/cpp/tests/linear_programming/pdlp_distributed_test.cu +++ b/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 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(&handle, problem); auto base = solve_lp(base_op, base_settings); diff --git a/cpp/tests/linear_programming/pdlp_test.cu b/cpp/tests/linear_programming/pdlp_test.cu index 715d43b7a4..311adefe61 100644 --- a/cpp/tests/linear_programming/pdlp_test.cu +++ b/cpp/tests/linear_programming/pdlp_test.cu @@ -511,8 +511,14 @@ TEST(pdlp_class, run_sub_mittleman) } } -constexpr double initial_step_size_afiro = 1.4893; -constexpr double initial_primal_weight_afiro = 0.0141652; +// Golden reference values for afiro's initial step size/primal weight, computed under +// cuOpt's current default scaling (Curtis-Reid -> Ruiz -> Pock-Chambolle). This test only +// cares about *whether* update_step_size_on_initial_solution/ +// update_primal_weight_on_initial_solution change these from their as-computed defaults, +// not their specific values, so these need re-baselining whenever the default scaling +// pipeline changes. +constexpr double initial_step_size_afiro = 1.402293; +constexpr double initial_primal_weight_afiro = 0.02019181; constexpr double factor_tolerance = 1e-4f; // Should be added to google test @@ -817,6 +823,93 @@ TEST(pdlp_class, initial_solution_test) } } +// Explicitly-stored zero coefficients (present in the CSR with value 0, not omitted) +// must not produce non-finite Curtis-Reid scale factors in float. The row/col kernels +// floor |a_ij| at std::numeric_limits::min() before taking a log so raft::log +// never sees 0. +TEST(pdlp_class, curtis_reid_scaling_explicit_zero_coefficient_float) +{ + const raft::handle_t handle_{}; + + cuopt::mathematical_optimization::optimization_problem_t op_problem(&handle_); + op_problem.set_maximize(false); + + // 1 constraint, 2 variables. Variable 1's coefficient is an *explicit* stored zero + // (present in the CSR with value 0.0, not simply omitted). + std::vector A_values = {1.0f, 0.0f}; + std::vector A_indices = {0, 1}; + std::vector A_offsets = {0, 2}; + op_problem.set_csr_constraint_matrix(A_values.data(), + static_cast(A_values.size()), + A_indices.data(), + static_cast(A_indices.size()), + A_offsets.data(), + static_cast(A_offsets.size())); + + std::vector constraint_lower = {0.0f}; + std::vector constraint_upper = {10.0f}; + op_problem.set_constraint_lower_bounds(constraint_lower.data(), + static_cast(constraint_lower.size())); + op_problem.set_constraint_upper_bounds(constraint_upper.data(), + static_cast(constraint_upper.size())); + + std::vector objective = {1.0f, 1.0f}; + op_problem.set_objective_coefficients(objective.data(), static_cast(objective.size())); + + std::vector var_lower = {0.0f, 0.0f}; + std::vector var_upper = {5.0f, 5.0f}; + op_problem.set_variable_lower_bounds(var_lower.data(), static_cast(var_lower.size())); + op_problem.set_variable_upper_bounds(var_upper.data(), static_cast(var_upper.size())); + + cuopt::mathematical_optimization::mip::problem_t problem(op_problem); + + auto solver_settings = pdlp_solver_settings_t{}; + // We only care about the scaling computed at construction time, not an actual solve. + solver_settings.iteration_limit = 0; + solver_settings.method = cuopt::mathematical_optimization::method_t::PDLP; + solver_settings.hyper_params.do_curtis_reid_scaling = true; + // Isolate Curtis-Reid: its own exp+clamp fold into the cumulative scale (clamp_bound = + // 30) would otherwise silently absorb a -inf/NaN log-domain value before it reaches the + // final scale factors, masking the bug this test targets. Checking the pre-fold + // log-domain values directly (via get_iteration_*_scaling() below) needs Ruiz/ + // Pock-Chambolle disabled so they don't overwrite those scratch buffers afterward. + solver_settings.hyper_params.do_ruiz_scaling = false; + solver_settings.hyper_params.do_pock_chambolle_scaling = false; + + // pdlp_solver_t's constructor builds a real pdhg_solver_t and wires it into the initial + // scaling strategy (running_mip=false, skip_ruiz_pock_compute=false), which runs + // compute_scaling_vectors() -- and therefore curtis_reid_scaling() -- right here. Going + // through pdlp_solver_t (rather than constructing pdlp_initial_scaling_strategy_t + // directly) avoids passing a null pdhg_solver_ptr, which trips its "PDHG solver pointer + // is null" assertion when running_mip is false. + cuopt::mathematical_optimization::pdlp::pdlp_solver_t solver(problem, + solver_settings); + auto& scaling = solver.get_initial_scaling_strategy(); + + // Pre-fold log-domain row/col scale (curtis_reid_scaling()'s direct output, before the + // exp+clamp that turns it into a multiplicative factor) -- this is what actually goes + // non-finite if raft::log() sees an unfloored zero. + auto row_log_scale = + host_copy(scaling.get_iteration_constraint_matrix_scaling(), handle_.get_stream()); + auto col_log_scale = host_copy(scaling.get_iteration_variable_scaling(), handle_.get_stream()); + for (float v : row_log_scale) { + EXPECT_TRUE(std::isfinite(v)) << "row log-scale is not finite: " << v; + } + for (float v : col_log_scale) { + EXPECT_TRUE(std::isfinite(v)) << "col log-scale is not finite: " << v; + } + + // Final (post exp+clamp) cumulative scale factors should also be finite. + auto row_scale = host_copy(scaling.get_constraint_matrix_scaling_vector(), handle_.get_stream()); + auto col_scale = host_copy(scaling.get_variable_scaling_vector(), handle_.get_stream()); + for (float v : row_scale) { + EXPECT_TRUE(std::isfinite(v)) << "row scale factor is not finite: " << v; + } + for (float v : col_scale) { + EXPECT_TRUE(std::isfinite(v)) << "col scale factor is not finite: " << v; + } +} + TEST(pdlp_class, initial_primal_weight_step_size_test) { const raft::handle_t handle_{}; @@ -1256,6 +1349,8 @@ TEST(pdlp_class, first_primal_feasible_stable3) solver_settings.pdlp_solver_mode = pdlp_solver_mode_t::Stable3; solver_settings.method = cuopt::mathematical_optimization::method_t::PDLP; solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; cuopt::mathematical_optimization::io::mps_data_model_t op_problem = cuopt::mathematical_optimization::io::read_mps(path); @@ -1470,6 +1565,8 @@ TEST(pdlp_class, first_primal_feasible_and_per_constraint_residual_stable3) solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; solver_settings.method = cuopt::mathematical_optimization::method_t::PDLP; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; cuopt::mathematical_optimization::io::mps_data_model_t op_problem = cuopt::mathematical_optimization::io::read_mps(path); @@ -1506,6 +1603,8 @@ TEST(pdlp_class, first_primal_feasible_and_per_constraint_residual_batch_stable3 constexpr double kOptimalityTolerance = 1e-2; solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; constexpr int batch_size = 2; @@ -1553,6 +1652,8 @@ TEST(pdlp_class, first_primal_feasible_and_per_constraint_residual_batch_differe constexpr double kOptimalityTolerance = 1e-2; solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; constexpr int batch_size = 2; @@ -1620,6 +1721,8 @@ TEST(pdlp_class, all_primal_feasible_and_per_constraint_residual_batch_different constexpr double kOptimalityTolerance = 1e-2; solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; constexpr int batch_size = 2; @@ -1687,6 +1790,8 @@ TEST(pdlp_class, all_primal_feasible_and_per_constraint_residual_batch_many_diff constexpr double kOptimalityTolerance = 1e-2; solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; const auto& original_lb = op_problem.get_constraint_lower_bounds(); const auto& original_ub = op_problem.get_constraint_upper_bounds(); @@ -1798,6 +1903,8 @@ TEST(pdlp_class, all_primal_feasible_and_per_constraint_residual_batch_many_diff constexpr double kOptimalityTolerance = 1e-2; solver_settings.set_optimality_tolerance(kOptimalityTolerance); solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; const auto& original_lb = op_problem.get_constraint_lower_bounds(); const auto& original_ub = op_problem.get_constraint_upper_bounds(); @@ -1951,6 +2058,8 @@ TEST(pdlp_class, warm_start) solver_settings.detect_infeasibility = false; solver_settings.method = cuopt::mathematical_optimization::method_t::PDLP; solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling changing iteration counts (woodlands09 only). + solver_settings.hyper_params.do_curtis_reid_scaling = (instance_name != "woodlands09"); cuopt::mathematical_optimization::io::mps_data_model_t mps_data_model = cuopt::mathematical_optimization::io::read_mps(path); @@ -2250,6 +2359,8 @@ TEST(pdlp_class, simple_batch_different_bounds) auto solver_settings = pdlp_solver_settings_t{}; solver_settings.method = cuopt::mathematical_optimization::method_t::PDLP; solver_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + solver_settings.hyper_params.do_curtis_reid_scaling = false; const std::vector& variable_lower_bounds = op_problem.get_variable_lower_bounds(); const std::vector& variable_upper_bounds = op_problem.get_variable_upper_bounds(); diff --git a/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py b/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py index ae58752ffb..66757e6f3e 100644 --- a/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py +++ b/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py @@ -24,6 +24,7 @@ CUOPT_ITERATION_LIMIT, CUOPT_METHOD, CUOPT_MIP_HEURISTICS_ONLY, + CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING, CUOPT_PDLP_SOLVER_MODE, CUOPT_PRIMAL_INFEASIBLE_TOLERANCE, CUOPT_RELATIVE_DUAL_TOLERANCE, @@ -530,6 +531,8 @@ def test_parse_var_names(): settings.set_parameter(CUOPT_METHOD, SolverMethod.PDLP) settings.set_parameter(CUOPT_PDLP_SOLVER_MODE, PDLPSolverMode.Stable2) settings.set_parameter(CUOPT_PRESOLVE, 0) + # Expected primal values below were recorded prior to implementing Curtis-Reid scaling. + settings.set_parameter(CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING, False) solution = solver.Solve(data_model_obj, settings) expected_dict = {