diff --git a/README.md b/README.md index ea14b47..4c065cc 100644 --- a/README.md +++ b/README.md @@ -87,7 +87,14 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--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` | -| `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | +| `--l_inf_ruiz_iter` | `int` | Iterations for L-inf Ruiz rescaling| `10` | +| `--no_pock_chambolle` | `flag` | Disable Pock-Chambolle rescaling | `enabled` | +| `--pock_chambolle_alpha` | `float` | Value for Pock-Chambolle alpha | `1.0` | +| `--no_bound_obj_rescaling` | `flag` | Disable bound objective rescaling | `enabled` | +| `--eval_freq` | `int` | Termination evaluation frequency | `200` | +| `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` | +| `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | +| `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | #### Output Files diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 23033bf..60aa5cf 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -85,6 +85,8 @@ extern "C" bool bound_objective_rescaling; bool verbose; int termination_evaluation_frequency; + int sv_max_iter; + double sv_tol; termination_criteria_t termination_criteria; restart_parameters_t restart_params; double reflection_coefficient; diff --git a/internal/solver.h b/internal/solver.h index ab88cb9..02b4eac 100644 --- a/internal/solver.h +++ b/internal/solver.h @@ -27,8 +27,6 @@ extern "C" const pdhg_parameters_t *params, const lp_problem_t *original_problem); - void set_default_parameters(pdhg_parameters_t *params); - #ifdef __cplusplus } #endif diff --git a/internal/utils.h b/internal/utils.h index 4998801..a6e36c4 100644 --- a/internal/utils.h +++ b/internal/utils.h @@ -140,6 +140,8 @@ extern "C" void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state); + void set_default_parameters(pdhg_parameters_t *params); + #ifdef __cplusplus } diff --git a/python/README.md b/python/README.md index 50c8fc9..9023910 100644 --- a/python/README.md +++ b/python/README.md @@ -155,6 +155,8 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `ReflectionCoeff` | `reflection_coefficient` | float | `1.0` | Reflection coefficient. | | `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.| | `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual. | +| `SVMaxIter` | `sv_max_iter` | int | 5000 | Maximum number of iterations for the power method | +| `SVTol`| `sv_tol` | float | `1e-4` | Termination tolerance for the power method | They can be set in multiple ways: diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index cc2733d..0d73b46 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -52,5 +52,8 @@ "ReflectionCoeff": "reflection_coefficient", # feasibility polishing "FeasibilityPolishing": "feasibility_polishing", - "FeasibilityPolishingTol": "eps_feas_polish_relative" + "FeasibilityPolishingTol": "eps_feas_polish_relative", + # singular value estimation (power method) + "SVMaxIter": "sv_max_iter", + "SVTol": "sv_tol", } \ No newline at end of file diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index 841b598..44e3e47 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -240,6 +240,10 @@ 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; + // power method for singular value estimation + d["sv_max_iter"] = p.sv_max_iter; + d["sv_tol"] = p.sv_tol; + return d; } @@ -288,6 +292,10 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p // Feasibility Polishing getb("feasibility_polishing", p->feasibility_polishing); getf("eps_feas_polish_relative", p->termination_criteria.eps_feas_polish_relative); + + // power method for singular value estimation + geti("sv_max_iter", p->sv_max_iter); + getf("sv_tol", p->sv_tol); } // view of matrix from Python diff --git a/src/cli.c b/src/cli.c index b9588ad..40855c3 100644 --- a/src/cli.c +++ b/src/cli.c @@ -161,24 +161,38 @@ void print_usage(const char *prog_name) fprintf(stderr, "Options:\n"); fprintf(stderr, " -h, --help Display this help message.\n"); - fprintf(stderr, " -v, --verbose Enable verbose " - "logging (default: false).\n"); - fprintf(stderr, " --time_limit Time limit in seconds " - "(default: 3600.0).\n"); + fprintf(stderr, " -v, --verbose " + "Enable verbose logging (default: false).\n"); + fprintf(stderr, " --time_limit " + "Time limit in seconds (default: 3600.0).\n"); fprintf( stderr, " --iter_limit Iteration limit (default: %d).\n", INT32_MAX); - fprintf(stderr, " --eps_opt Relative optimality " - "tolerance (default: 1e-4).\n"); - fprintf(stderr, " --eps_feas Relative feasibility " - "tolerance (default: 1e-4).\n"); - fprintf(stderr, " --eps_infeas_detect Infeasibility " - "detection tolerance (default: 1e-10).\n"); + fprintf(stderr, " --eps_opt " + "Relative optimality tolerance (default: 1e-4).\n"); + fprintf(stderr, " --eps_feas " + "Relative feasibility tolerance (default: 1e-4).\n"); + fprintf(stderr, " --eps_infeas_detect " + "Infeasibility detection tolerance (default: 1e-10).\n"); + fprintf(stderr, " --l_inf_ruiz_iter " + "Iterations for L-inf Ruiz rescaling (default: 10).\n"); + fprintf(stderr, " --no_pock_chambolle " + "Disable Pock-Chambolle rescaling (default: enabled).\n"); + fprintf(stderr, " --pock_chambolle_alpha " + "Value for Pock-Chambolle alpha (default: 1.0).\n"); + fprintf(stderr, " --no_bound_obj_rescaling " + "Disable bound objective rescaling (default: enabled).\n"); + fprintf(stderr, " --eval_freq " + "Termination evaluation frequency (default: 200).\n"); + fprintf(stderr, " --sv_max_iter " + "Max iterations for singular value estimation (default: 5000).\n"); + fprintf(stderr, " --sv_tol " + "Tolerance for singular value estimation (default: 1e-4).\n"); + fprintf(stderr, " -f --feasibility_polishing " + "Enable feasibility use feasibility polishing (default: false).\n"); fprintf(stderr, " --eps_feas_polish Relative feasibility " "polish tolerance (default: 1e-6).\n"); - fprintf(stderr, " -f --feasibility_polishing Enable feasibility " - "use feasibility polishing (default: false).\n"); } int main(int argc, char *argv[]) @@ -196,6 +210,13 @@ int main(int argc, char *argv[]) {"eps_infeas_detect", required_argument, 0, 1005}, {"eps_feas_polish", required_argument, 0, 1006}, {"feasibility_polishing", no_argument, 0, 'f'}, + {"l_inf_ruiz_iter", required_argument, 0, 1007}, + {"pock_chambolle_alpha", required_argument, 0, 1008}, + {"no_pock_chambolle", no_argument, 0, 1009}, + {"no_bound_obj_rescaling", no_argument, 0, 1010}, + {"sv_max_iter", required_argument, 0, 1011}, + {"sv_tol", required_argument, 0, 1012}, + {"eval_freq", required_argument, 0, 1013}, {0, 0, 0, 0}}; int opt; @@ -227,9 +248,30 @@ int main(int argc, char *argv[]) case 1006: // --eps_feas_polish_relative params.termination_criteria.eps_feas_polish_relative = atof(optarg); break; - case 'f': // --feasibility_polishing + case 'f': // --feasibility_polishing params.feasibility_polishing = true; break; + case 1007: // --l_inf_ruiz_iter + params.l_inf_ruiz_iterations = atoi(optarg); + break; + case 1008: // --pock_chambolle_alpha + params.pock_chambolle_alpha = atof(optarg); + break; + case 1009: // --no_pock_chambolle + params.has_pock_chambolle_alpha = false; + break; + case 1010: // --no_bound_obj_rescaling + params.bound_objective_rescaling = false; + break; + case 1011: // --sv_max_iter + params.sv_max_iter = atoi(optarg); + break; + case 1012: // --sv_tol + params.sv_tol = atof(optarg); + break; + case 1013: // --eval_freq + params.termination_evaluation_frequency = atoi(optarg); + break; case '?': // Unknown option return 1; } diff --git a/src/solver.cu b/src/solver.cu index a9e5f2b..f918d62 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -760,7 +760,7 @@ initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, { double max_sv = estimate_maximum_singular_value( state->sparse_handle, state->blas_handle, state->constraint_matrix, - state->constraint_matrix_t, 5000, 1e-4); + state->constraint_matrix_t, params->sv_max_iter, params->sv_tol); state->step_size = 0.998 / max_sv; if (params->bound_objective_rescaling) @@ -952,32 +952,6 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state) return results; } -void set_default_parameters(pdhg_parameters_t *params) -{ - params->l_inf_ruiz_iterations = 10; - params->has_pock_chambolle_alpha = true; - params->pock_chambolle_alpha = 1.0; - params->bound_objective_rescaling = true; - params->verbose = false; - params->termination_evaluation_frequency = 200; - params->feasibility_polishing = false; - params->reflection_coefficient = 1.0; - - params->termination_criteria.eps_optimal_relative = 1e-4; - params->termination_criteria.eps_feasible_relative = 1e-4; - params->termination_criteria.eps_infeasible = 1e-10; - params->termination_criteria.time_sec_limit = 3600.0; - params->termination_criteria.iteration_limit = INT32_MAX; - params->termination_criteria.eps_feas_polish_relative = 1e-6; - - params->restart_params.artificial_restart_threshold = 0.36; - params->restart_params.sufficient_reduction_for_restart = 0.2; - params->restart_params.necessary_reduction_for_restart = 0.5; - params->restart_params.k_p = 0.99; - params->restart_params.k_i = 0.01; - params->restart_params.k_d = 0.0; - params->restart_params.i_smooth = 0.3; -} //Feasibility Polishing void feasibility_polish(const pdhg_parameters_t *params, pdhg_solver_state_t *state) diff --git a/src/utils.cu b/src/utils.cu index 172e731..4135fe9 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -317,9 +317,62 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, return do_restart; } +void set_default_parameters(pdhg_parameters_t *params) +{ + params->l_inf_ruiz_iterations = 10; + params->has_pock_chambolle_alpha = true; + params->pock_chambolle_alpha = 1.0; + params->bound_objective_rescaling = true; + params->verbose = false; + params->termination_evaluation_frequency = 200; + params->feasibility_polishing = false; + params->reflection_coefficient = 1.0; + + params->sv_max_iter = 5000; + params->sv_tol = 1e-4; + + params->termination_criteria.eps_optimal_relative = 1e-4; + params->termination_criteria.eps_feasible_relative = 1e-4; + params->termination_criteria.eps_infeasible = 1e-10; + params->termination_criteria.time_sec_limit = 3600.0; + params->termination_criteria.iteration_limit = INT32_MAX; + params->termination_criteria.eps_feas_polish_relative = 1e-6; + + params->restart_params.artificial_restart_threshold = 0.36; + params->restart_params.sufficient_reduction_for_restart = 0.2; + params->restart_params.necessary_reduction_for_restart = 0.5; + params->restart_params.k_p = 0.99; + params->restart_params.k_i = 0.01; + params->restart_params.k_d = 0.0; + params->restart_params.i_smooth = 0.3; +} + +#define PRINT_DIFF_INT(name, current, default_val) \ + do { \ + if ((current) != (default_val)) { \ + printf(" %-18s : %d\n", name, current); \ + } \ + } while(0) + +#define PRINT_DIFF_DBL(name, current, default_val) \ + do { \ + if (fabs((current) - (default_val)) > 1e-9) { \ + printf(" %-18s : %.1e\n", name, (double)(current)); \ + } \ + } while(0) + +#define PRINT_DIFF_BOOL(name, current, default_val) \ + do { \ + if ((current) != (default_val)) { \ + printf(" %-18s : %s\n", name, (current) ? "on" : "off"); \ + } \ + } while(0) + void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem) { + pdhg_parameters_t default_params; + set_default_parameters(&default_params); if (!params->verbose) { return; @@ -339,7 +392,7 @@ void print_initial_info(const pdhg_parameters_t *params, printf("problem:\n"); printf(" variables : %d\n", problem->num_variables); printf(" constraints : %d\n", problem->num_constraints); - printf(" nnz(A) : %d\n", problem->constraint_matrix_num_nonzeros); + printf(" nonzeros(A) : %d\n", problem->constraint_matrix_num_nonzeros); printf("settings:\n"); printf(" iter_limit : %d\n", @@ -352,6 +405,33 @@ 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); + PRINT_DIFF_INT("l_inf_ruiz_iter", + params->l_inf_ruiz_iterations, + default_params.l_inf_ruiz_iterations); + PRINT_DIFF_DBL("pock_chambolle_alpha", + params->pock_chambolle_alpha, + default_params.pock_chambolle_alpha); + PRINT_DIFF_BOOL("has_pock_chambolle_alpha", + params->has_pock_chambolle_alpha, + default_params.has_pock_chambolle_alpha); + PRINT_DIFF_BOOL("bound_obj_rescaling", + params->bound_objective_rescaling, + default_params.bound_objective_rescaling); + PRINT_DIFF_INT("sv_max_iter", + params->sv_max_iter, + default_params.sv_max_iter); + PRINT_DIFF_DBL("sv_tol", + params->sv_tol, + default_params.sv_tol); + PRINT_DIFF_INT("evaluation_freq", + params->termination_evaluation_frequency, + default_params.termination_evaluation_frequency); + PRINT_DIFF_BOOL("feasibility_polishing", + params->feasibility_polishing, + default_params.feasibility_polishing); + PRINT_DIFF_DBL("eps_feas_polish_relative", + params->termination_criteria.eps_feas_polish_relative, + default_params.termination_criteria.eps_feas_polish_relative); printf("---------------------------------------------------------------------" "------------------\n"); @@ -364,6 +444,10 @@ void print_initial_info(const pdhg_parameters_t *params, "------------------\n"); } +#undef PRINT_DIFF_INT +#undef PRINT_DIFF_DBL +#undef PRINT_DIFF_BOOL + void pdhg_final_log(const pdhg_solver_state_t *state, bool verbose, termination_reason_t reason) {