diff --git a/.gitignore b/.gitignore index cd9fb58..866c33f 100644 --- a/.gitignore +++ b/.gitignore @@ -41,6 +41,7 @@ *.su *.idb *.pdb +Testing/* # Kernel Module Compile Results *.mod* @@ -60,3 +61,5 @@ test/* /.vscode /.venv /_b +*.whl +*.pyc diff --git a/README.md b/README.md index 320a9eb..daf779b 100644 --- a/README.md +++ b/README.md @@ -84,6 +84,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `-v`, `--verbose` | `flag` | Enable verbose logging. | `false` | | `--time_limit` | `double` | Time limit in seconds. | `3600.0` | | `--iter_limit` | `int` | Iteration limit. | `2147483647` | +| `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` | | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` | | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` | | `--eps_infeas_detect` | `double` | Infeasibility detection tolerance. | `1e-10` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index b2cb101..c213b97 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -37,6 +37,12 @@ extern "C" TERMINATION_REASON_FEAS_POLISH_SUCCESS } termination_reason_t; + typedef enum + { + NORM_TYPE_L2 = 0, + NORM_TYPE_L_INF = 1 + } norm_type_t; + typedef struct { int num_variables; @@ -93,6 +99,7 @@ extern "C" restart_parameters_t restart_params; double reflection_coefficient; bool feasibility_polishing; + norm_type_t optimality_norm; bool presolve; } pdhg_parameters_t; diff --git a/internal/utils.h b/internal/utils.h index 7889869..1b6b7c9 100644 --- a/internal/utils.h +++ b/internal/utils.h @@ -107,7 +107,7 @@ extern "C" int get_print_frequency(int iter); - void compute_residual(pdhg_solver_state_t *state); + void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm); void compute_infeasibility_information(pdhg_solver_state_t *state); @@ -133,9 +133,9 @@ extern "C" void pdhg_feas_polish_final_log(const pdhg_solver_state_t *primal_state, const pdhg_solver_state_t *dual_state, bool verbose); - void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state); + void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, norm_type_t optimality_norm); - void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state); + void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, norm_type_t optimality_norm); void set_default_parameters(pdhg_parameters_t *params); diff --git a/python/README.md b/python/README.md index 6e8e042..c5c1342 100644 --- a/python/README.md +++ b/python/README.md @@ -141,6 +141,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `IterationLimit` | `iteration_limit` | int | `2147483647` | Maximum number of iterations. | | `OutputFlag`, `LogToConsole` | `verbose` | bool | `False` | Enable (`True`) or disable (`False`) console logging output. | | `TermCheckFreq` | `termination_evaluation_frequency` | int | `200` | Frequency (in iterations) at which termination conditions are evaluated. | +| `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. | | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. | | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. | | `InfeasibleTol` | `eps_infeasible` | float | `1e-10` | Threshold for declaring infeasibility. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 2a56350..0d2eec7 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -53,6 +53,8 @@ # feasibility polishing "FeasibilityPolishing": "feasibility_polishing", "FeasibilityPolishingTol": "eps_feas_polish_relative", + # termination criteria + "OptimalityNorm": "optimality_norm", # singular value estimation (power method) "SVMaxIter": "sv_max_iter", "SVTol": "sv_tol", diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index c22d66b..b18aba7 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -137,6 +137,21 @@ static const int32_t *get_index_ptr_i32(py::object obj, const char *name, throw std::invalid_argument(std::string(name) + " must be int32 or int64."); } +// helper function to convert norm string to enum +static norm_type_t parse_norm_string(const std::string& s) +{ + std::string lower = s; + std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower); + + if (lower == "l2") { + return NORM_TYPE_L2; + } else if (lower == "linf") { + return NORM_TYPE_L_INF; + } else { + throw std::invalid_argument("Unknown norm type: " + s + ". Use 'l2' or 'linf'."); + } +} + // ensure 1D array or None with expected length static void ensure_len_or_null(py::object obj, const char *name, int expect_len) { @@ -247,6 +262,8 @@ static py::dict get_default_params_py() d["feasibility_polishing"] = p.feasibility_polishing; d["eps_feas_polish_relative"] = p.termination_criteria.eps_feas_polish_relative; + // Termination criteria norm + d["optimality_norm"] = (p.optimality_norm == NORM_TYPE_L_INF) ? "linf" : "l2"; // power method for singular value estimation d["sv_max_iter"] = p.sv_max_iter; d["sv_tol"] = p.sv_tol; @@ -270,6 +287,19 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p { if (d.contains(k)) tgt = py::cast(d[k]); }; auto getb = [&](const char *k, bool &tgt) { if (d.contains(k)) tgt = py::cast(d[k]); }; + auto get_norm = [&](const char *k, norm_type_t &tgt) + { + if (d.contains(k)) { + py::object val = d[k]; + if (py::isinstance(val)) { + std::string sval = py::cast(val); + tgt = parse_norm_string(sval); + } + else { + throw std::invalid_argument("optimality_norm must be a string ('l2'/'linf')"); + } + } + }; // verbosity getb("verbose", p->verbose); @@ -303,6 +333,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getb("feasibility_polishing", p->feasibility_polishing); getf("eps_feas_polish_relative", p->termination_criteria.eps_feas_polish_relative); + // Termination criteria norm + get_norm("optimality_norm", p->optimality_norm); // power method for singular value estimation geti("sv_max_iter", p->sv_max_iter); getf("sv_tol", p->sv_tol); diff --git a/src/cli.c b/src/cli.c index 690a9ff..f0bfe4e 100644 --- a/src/cli.c +++ b/src/cli.c @@ -201,6 +201,8 @@ void print_usage(const char *prog_name) "Enable feasibility use feasibility polishing (default: false).\n"); fprintf(stderr, " --eps_feas_polish Relative feasibility " "polish tolerance (default: 1e-6).\n"); + fprintf(stderr, " --opt_norm " + "Norm for optimality criteria: l2 or linf (default: l2).\n"); fprintf(stderr, " --no_presolve " "Disable presolve (default: enabled).\n"); } @@ -227,7 +229,8 @@ int main(int argc, char *argv[]) {"sv_max_iter", required_argument, 0, 1011}, {"sv_tol", required_argument, 0, 1012}, {"eval_freq", required_argument, 0, 1013}, - {"no_presolve", no_argument, 0, 1014}, + {"opt_norm", required_argument, 0, 1014}, + {"no_presolve", no_argument, 0, 1015}, {0, 0, 0, 0}}; int opt; @@ -283,7 +286,20 @@ int main(int argc, char *argv[]) case 1013: // --eval_freq params.termination_evaluation_frequency = atoi(optarg); break; - case 1014: // --no_presolve + case 1014: // --opt_norm + { + const char *norm_str = optarg; + if (strcmp(norm_str, "l2") == 0) { + params.optimality_norm = NORM_TYPE_L2; + } else if (strcmp(norm_str, "linf") == 0) { + params.optimality_norm = NORM_TYPE_L_INF; + } else { + fprintf(stderr, "Error: opt_norm must be 'l2' or 'linf'\n"); + return 1; + } + } + break; + case 1015: // --no_presolve params.presolve = false; break; case '?': // Unknown option diff --git a/src/solver.cu b/src/solver.cu index 0e58dbf..3e61fad 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -125,7 +125,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, if ((state->is_this_major_iteration || state->total_count == 0) || (state->total_count % get_print_frequency(state->total_count) == 0)) { - compute_residual(state); + compute_residual(state, params->optimality_norm); if (state->is_this_major_iteration && state->total_count < 3 * params->termination_evaluation_frequency) { @@ -176,7 +176,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, if (state->termination_reason == TERMINATION_REASON_UNSPECIFIED) { state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; - compute_residual(state); + compute_residual(state, params->optimality_norm); display_iteration_stats(state, params->verbose); } @@ -408,32 +408,59 @@ initialize_solver_state(const pdhg_parameters_t *params, free(temp_host); double sum_of_squares = 0.0; + double max_val = 0.0; + double val = 0.0; for (int i = 0; i < n_vars; ++i) { - sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i]; + if (params->optimality_norm == NORM_TYPE_L_INF) { + val = fabs(working_problem->objective_vector[i]); + if (val > max_val) max_val = val; + } else { + sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i]; + } + } + + if (params->optimality_norm == NORM_TYPE_L_INF) { + state->objective_vector_norm = max_val; + } else { + state->objective_vector_norm = sqrt(sum_of_squares); } - state->objective_vector_norm = sqrt(sum_of_squares); sum_of_squares = 0.0; + max_val = 0.0; + val = 0.0; for (int i = 0; i < n_cons; ++i) { double lower = working_problem->constraint_lower_bound[i]; double upper = working_problem->constraint_upper_bound[i]; - if (isfinite(lower) && (lower != upper)) - { - sum_of_squares += lower * lower; + if (params->optimality_norm == NORM_TYPE_L_INF) { + if (isfinite(lower) && (lower != upper)) { + val = fabs(lower); + if (val > max_val) max_val = val; + } + if (isfinite(upper)) { + val = fabs(upper); + if (val > max_val) max_val = val; + } + } else { + if (isfinite(lower) && (lower != upper)) { + sum_of_squares += lower * lower; + } + if (isfinite(upper)) { + sum_of_squares += upper * upper; + } } + } - if (isfinite(upper)) - { - sum_of_squares += upper * upper; - } + if (params->optimality_norm == NORM_TYPE_L_INF) { + state->constraint_bound_norm = max_val; + } else { + state->constraint_bound_norm = sqrt(sum_of_squares); } - state->constraint_bound_norm = sqrt(sum_of_squares); state->num_blocks_primal = (state->num_variables + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK; state->num_blocks_dual = @@ -1129,7 +1156,7 @@ void primal_feasibility_polish(const pdhg_parameters_t *params, pdhg_solver_stat { if ((state->is_this_major_iteration || state->total_count == 0) || (state->total_count % get_print_frequency(state->total_count) == 0)) { - compute_primal_feas_polish_residual(state, ori_state); + compute_primal_feas_polish_residual(state, ori_state, params->optimality_norm); state->cumulative_time_sec = (double)(clock() - start_time) / CLOCKS_PER_SEC; @@ -1175,7 +1202,7 @@ void dual_feasibility_polish(const pdhg_parameters_t *params, pdhg_solver_state_ { if ((state->is_this_major_iteration || state->total_count == 0) || (state->total_count % get_print_frequency(state->total_count) == 0)) { - compute_dual_feas_polish_residual(state, ori_state); + compute_dual_feas_polish_residual(state, ori_state, params->optimality_norm); state->cumulative_time_sec = (double)(clock() - start_time) / CLOCKS_PER_SEC; diff --git a/src/utils.cu b/src/utils.cu index b86b307..2029de9 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -224,8 +224,8 @@ bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol) { return state->relative_dual_residual < rel_feas_tol && - state->relative_primal_residual < rel_feas_tol && - state->relative_objective_gap < rel_opt_tol; + state->relative_primal_residual < rel_feas_tol && + state->relative_objective_gap < rel_opt_tol; } bool primal_infeasibility_criteria_met(const pdhg_solver_state_t *state, @@ -349,6 +349,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->restart_params.k_d = 0.0; params->restart_params.i_smooth = 0.3; + params->optimality_norm = NORM_TYPE_L2; params->presolve = true; } @@ -407,6 +408,10 @@ void print_initial_info(const pdhg_parameters_t *params, params->termination_criteria.eps_feasible_relative); printf(" eps_infeas_detect : %.1e\n", params->termination_criteria.eps_infeasible); + if (params->optimality_norm != default_params.optimality_norm) { + printf(" optimality_norm : %s\n", + params->optimality_norm == NORM_TYPE_L_INF ? "L_inf" : "L2"); + } PRINT_DIFF_INT("l_inf_ruiz_iter", params->l_inf_ruiz_iterations, @@ -680,7 +685,7 @@ static double get_vector_sum(cublasHandle_t handle, int n, double *ones_d, return sum; } -void compute_residual(pdhg_solver_state_t *state) +void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm) { cusparseDnVecSetValues(state->vec_primal_sol, state->pdhg_primal_solution); cusparseDnVecSetValues(state->vec_dual_sol, state->pdhg_dual_solution); @@ -707,13 +712,26 @@ void compute_residual(pdhg_solver_state_t *state) state->constraint_upper_bound_finite_val, state->num_constraints, state->num_variables); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, - state->primal_residual, 1, - &state->absolute_primal_residual)); + if (optimality_norm == NORM_TYPE_L_INF) { + state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, + state->num_constraints, state->primal_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, + state->primal_residual, 1, + &state->absolute_primal_residual)); + } + state->absolute_primal_residual /= state->constraint_bound_rescaling; - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, - state->dual_residual, 1, - &state->absolute_dual_residual)); + + if (optimality_norm == NORM_TYPE_L_INF) { + state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, + state->num_variables, state->dual_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, + state->dual_residual, 1, + &state->absolute_dual_residual)); + } + state->absolute_dual_residual /= state->objective_vector_rescaling; CUBLAS_CHECK(cublasDdot( @@ -736,12 +754,15 @@ void compute_residual(pdhg_solver_state_t *state) state->objective_vector_rescaling) + state->objective_constant; - state->relative_primal_residual = + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); + state->objective_gap = fabs(state->primal_objective_value - state->dual_objective_value); + state->relative_objective_gap = state->objective_gap / (1.0 + fabs(state->primal_objective_value) + fabs(state->dual_objective_value)); @@ -1218,7 +1239,7 @@ __global__ void compute_dual_feas_polish_residual_kerenl( } } -void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state) +void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, norm_type_t optimality_norm) { cusparseDnVecSetValues(state->vec_primal_sol, state->pdhg_primal_solution); cusparseDnVecSetValues(state->vec_primal_prod, state->primal_product); @@ -1230,15 +1251,24 @@ void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_ state->constraint_upper_bound, state->constraint_rescaling, state->num_constraints); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->primal_residual, 1, &state->absolute_primal_residual)); + if (optimality_norm == NORM_TYPE_L_INF) { + state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, + state->num_constraints, state->primal_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, + state->primal_residual, 1, + &state->absolute_primal_residual)); + } + state->absolute_primal_residual /= state->constraint_bound_rescaling; + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); CUBLAS_CHECK(cublasDdot(state->blas_handle, state->num_variables, ori_state->objective_vector, 1, state->pdhg_primal_solution, 1, &state->primal_objective_value)); state->primal_objective_value = state->primal_objective_value / (state->constraint_bound_rescaling * state->objective_vector_rescaling) + state->objective_constant; } -void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state) +void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, norm_type_t optimality_norm) { cusparseDnVecSetValues(state->vec_dual_sol, state->pdhg_dual_solution); cusparseDnVecSetValues(state->vec_dual_prod, state->dual_product); @@ -1256,8 +1286,16 @@ void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_so ori_state->constraint_upper_bound_finite_val, state->num_variables, state->num_constraints); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->dual_residual, 1, &state->absolute_dual_residual)); + if (optimality_norm == NORM_TYPE_L_INF) { + state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, + state->num_variables, state->dual_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, + state->dual_residual, 1, &state->absolute_dual_residual)); + } + state->absolute_dual_residual /= state->objective_vector_rescaling; + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); double base_dual_objective; diff --git a/test/test.sh b/test/test.sh old mode 100644 new mode 100755