From e4ae4328c27e92810487ee62323ab45a4571f3cb Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Thu, 10 Sep 2026 09:24:08 -0400 Subject: [PATCH] feat: active-set stepsize boosting --- README.md | 14 + include/cupdlpx.h | 4 + include/cupdlpx_types.h | 17 ++ internal/active_set_boost.h | 51 ++++ internal/cuda_to_hip.h | 3 + internal/internal_types.h | 66 ++++ internal/utils.h | 22 +- python/README.md | 12 + python/cupdlpx/PDLP.py | 13 + python/cupdlpx/_core.py | 2 +- python/cupdlpx/model.py | 72 ++--- python_bindings/_core_bindings.cpp | 131 +++++--- src/active_set_boost.cu | 458 ++++++++++++++++++++++++++++ src/cli.c | 87 ++++++ src/cupdlpx.c | 83 +++++ src/feasibility_polish.cu | 12 +- src/solver.cu | 139 +++++---- src/utils.cu | 468 +++++++++++++++++++++-------- test/test_active_set_boost.py | 224 ++++++++++++++ test/test_api_surface.py | 6 + test/test_interface.c | 2 +- 21 files changed, 1610 insertions(+), 276 deletions(-) create mode 100644 internal/active_set_boost.h create mode 100644 src/active_set_boost.cu create mode 100644 test/test_active_set_boost.py diff --git a/README.md b/README.md index 517fdaa..b14f5dc 100644 --- a/README.md +++ b/README.md @@ -103,6 +103,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `-h`, `--help` | `flag` | Display the help message. | N/A | | `-v`, `--verbose` | `flag` | Verbose logging (enabled by default). | `true` | | `-q`, `--quiet` | `flag` | Disable verbose logging. | `false` | +| `--debug` | `flag` | Developer diagnostics (implies verbose). | `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` | @@ -119,6 +120,19 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--no_presolve` | `flag` | Disable presolve | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | +| `--no_active_set_boost` | `flag` | Disable the active-set stepsize boost | `enabled` | +| `--asb_activation_tol` | `double` | Residual threshold at which the boost activates | `1e-4` | +| `--asb_window_iter` | `int` | Number of recent iterations used to identify the active set | `10000` | +| `--asb_safety_factor` | `double` | Boosted step = factor / estimated singular value | `0.9` | +| `--asb_max_reverts` | `int` | Divergences tolerated before the boost turns off; a diverged step is always reverted | `2` | +| `--asb_min_raise_ratio` | `double` | Minimum ratio for a step increase | `1.1` | +| `--asb_reestimate_change_ratio` | `double` | Fraction of the active set that must change before re-estimating | `0.01` | +| `--asb_constraint_tol` | `double` | Tolerance for treating a constraint as binding | `1e-8` | +| `--asb_variable_tol` | `double` | Tolerance for treating a variable as at its bound | `1e-8` | +| `--asb_divergence_ceiling_ratio` | `double` | Step ceiling after a revert, relative to the diverged step | `0.7` | +| `--asb_divergence_margin` | `double` | Allowed fixed-point error increase before a revert | `0.05` | + +The active-set boost enlarges the stepsize late in the solve, based on the constraints and bounds identified as active. #### Output Files The solver generates three text files in the specified . The filenames are derived from the input file's basename. For an input `INSTANCE.mps.gz`, the output will be: diff --git a/include/cupdlpx.h b/include/cupdlpx.h index 05839d3..f556fc5 100644 --- a/include/cupdlpx.h +++ b/include/cupdlpx.h @@ -42,6 +42,10 @@ extern "C" // parameter void set_default_parameters(pdhg_parameters_t *params); + /* Return 0 if params are valid. Otherwise, return nonzero and write the first + error to error_message when provided. */ + int cupdlpx_validate_parameters(const pdhg_parameters_t *params, char *error_message, size_t error_message_size); + void cupdlpx_result_free(cupdlpx_result_t *results); void lp_problem_free(lp_problem_t *prob); diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 16d3f04..daf03a4 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -99,6 +99,7 @@ extern "C" double pock_chambolle_alpha; bool bound_objective_rescaling; bool verbose; + bool debug; int termination_evaluation_frequency; int sv_max_iter; double sv_tol; @@ -111,6 +112,17 @@ extern "C" double matrix_zero_tol; double infinite_bound; int geometric_mean_iterations; + bool active_set_boost; + double asb_activation_tol; + int asb_window_iter; + double asb_safety_factor; + int asb_max_reverts; + double asb_min_raise_ratio; + double asb_reestimate_change_ratio; + double asb_constraint_tol; + double asb_variable_tol; + double asb_divergence_ceiling_ratio; + double asb_divergence_margin; } pdhg_parameters_t; typedef struct @@ -149,6 +161,11 @@ extern "C" termination_reason_t termination_reason; double feasibility_polishing_time; int feasibility_iteration; + + /* active-set step controller statistics (zero when the controller is disabled) */ + int asb_raise_count; + int asb_revert_count; + int asb_pi_iterations; } cupdlpx_result_t; // matrix formats diff --git a/internal/active_set_boost.h b/internal/active_set_boost.h new file mode 100644 index 0000000..2a9a938 --- /dev/null +++ b/internal/active_set_boost.h @@ -0,0 +1,51 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#pragma once + +#include "internal_types.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + + typedef enum + { + ASB_PHASE_WAITING = 0, + ASB_PHASE_ACTIVE = 1, + ASB_PHASE_OFF = 2 + } asb_phase_t; + + typedef enum + { + ASB_ACTION_NONE = 0, + ASB_ACTION_REVERT = 1 + } asb_action_t; + + void active_set_boost_init(pdhg_solver_state_t *state); + + void active_set_boost_free(pdhg_solver_state_t *state); + + void active_set_boost_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params); + + asb_action_t active_set_boost_check(pdhg_solver_state_t *state, const pdhg_parameters_t *params); + + void active_set_boost_on_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params); + +#ifdef __cplusplus +} +#endif diff --git a/internal/cuda_to_hip.h b/internal/cuda_to_hip.h index fc9975e..ced2bf8 100644 --- a/internal/cuda_to_hip.h +++ b/internal/cuda_to_hip.h @@ -50,6 +50,7 @@ limitations under the License. #define cudaMemcpy hipMemcpy #define cudaMemcpyAsync hipMemcpyAsync #define cudaMemset hipMemset +#define cudaMemsetAsync hipMemsetAsync // Memory copy kinds #define cudaMemcpyHostToDevice hipMemcpyHostToDevice @@ -66,6 +67,7 @@ limitations under the License. #define cudaStream_t hipStream_t #define cudaStreamCreate hipStreamCreate #define cudaStreamDestroy hipStreamDestroy +#define cudaStreamSynchronize hipStreamSynchronize // Device synchronization #define cudaDeviceSynchronize hipDeviceSynchronize @@ -127,6 +129,7 @@ static inline const char *cublasGetStatusName(hipblasStatus_t status) #define cusparseCreate hipsparseCreate #define cusparseDestroy hipsparseDestroy #define cusparseSetStream hipsparseSetStream +#define cusparseGetStream hipsparseGetStream #define cusparseStatus_t hipsparseStatus_t #define CUSPARSE_STATUS_SUCCESS HIPSPARSE_STATUS_SUCCESS #define cusparseGetErrorName hipsparseGetErrorName diff --git a/internal/internal_types.h b/internal/internal_types.h index bde5a4e..04e8b44 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -37,6 +37,37 @@ typedef struct int *transpose_map; } cu_sparse_matrix_csr_t; +/* Power-iteration estimator of the maximum singular value. The context owns the + device workspace and SpMV plans for one (A, AT) pair and is defined in utils.cu; + after a non-degenerate run it warm-starts the next run from its final eigenvector. */ +typedef struct sv_estimator_ctx sv_estimator_ctx_t; + +/* Per-run options; zero-initialize and set what you need. */ +typedef struct +{ + int max_iterations; + double tolerance; + const bool *d_row_mask; /* NULL = unmasked */ + const bool *d_col_mask; /* NULL = unmasked */ + double abort_singular_value_threshold; /* > 0: stop once the running Rayleigh lower + bound reaches it; 0 = never */ +} sv_estimator_opts_t; + +typedef enum +{ + SV_ESTIMATOR_CONVERGED = 0, /* residual test passed */ + SV_ESTIMATOR_ABORTED = 1, /* abort threshold crossed; the estimate is a lower bound */ + SV_ESTIMATOR_MAX_ITER = 2, /* max_iterations exhausted without convergence */ + SV_ESTIMATOR_DEGENERATE = 3, /* zero or non-finite iterate; max_singular_value is 0.0 */ +} sv_estimator_status_t; + +typedef struct +{ + sv_estimator_status_t status; + double max_singular_value; + int iterations; /* power iterations actually run */ +} sv_estimator_result_t; + typedef struct { int num_variables; @@ -72,6 +103,7 @@ typedef struct double *reflected_dual_solution; double *primal_product; double step_size; + double base_step_size; double *d_primal_step_size; double *d_dual_step_size; double primal_weight; @@ -88,6 +120,7 @@ typedef struct double objective_vector_rescaling; double *primal_slack; double *dual_slack; + double *infeasibility_dual_scratch; /* n-vector scratch for the ray certificate; dual_slack must survive */ double rescaling_time_sec; clock_t start_time; double cumulative_time_sec; @@ -114,8 +147,41 @@ typedef struct double initial_fixed_point_error; double last_trial_fixed_point_error; int inner_count; + int restart_count; /* restarts performed, ASB reverts included */ + const char *last_restart_reason; /* criterion that fired for the most recent adaptive restart */ + int logged_restart_count; /* restart_count at the last printed iteration row */ int *d_inner_count; + /* active-set step controller */ + sv_estimator_ctx_t *asb_sv_ctx; /* cached power-iteration workspace + SpMV plans */ + int asb_phase; + double asb_sv; + double asb_step_ceiling; /* post-divergence cap on the target step (INFINITY = none) */ + double asb_anchor_primal_weight; + /* primal-weight PID state captured with the anchor; a revert restores it too */ + double asb_anchor_pw_error_sum; + double asb_anchor_pw_last_error; + double asb_anchor_best_pw; + double asb_anchor_best_pd_residual_gap; + int asb_revert_count; + int asb_raise_count; + int asb_pi_early_exit_count; + int asb_pi_iterations; + int asb_sv_failed_count; + bool asb_no_raise_certified; /* a NO_RAISE abort covers the unchanged union */ + long asb_changes_since_estimate; /* mask entries added or removed since the last sv estimate */ + long asb_free_variables; /* variables in the mask: not confidently clamped within the window */ + long asb_binding_constraints; /* constraints in the mask: binding within the window */ + int *d_asb_var_last_free; + int *d_asb_row_last_binding; + bool *d_asb_col_mask; + bool *d_asb_row_mask; + double *d_asb_primal_anchor; + double *d_asb_dual_anchor; + double *d_asb_dual_slack_anchor; + double *d_asb_dual_projection_input; + int *d_asb_count; + cusparseHandle_t sparse_handle; cublasHandle_t blas_handle; void *spmv_ctx; diff --git a/internal/utils.h b/internal/utils.h index dc04747..303e4db 100644 --- a/internal/utils.h +++ b/internal/utils.h @@ -71,12 +71,14 @@ extern "C" void *safe_realloc(void *ptr, size_t new_size); - double estimate_maximum_singular_value(cusparseHandle_t sparse_handle, - cublasHandle_t blas_handle, - const cu_sparse_matrix_csr_t *A, - const cu_sparse_matrix_csr_t *AT, - int max_iterations, - double tolerance); + sv_estimator_ctx_t *sv_estimator_create(cusparseHandle_t sparse_handle, + cublasHandle_t blas_handle, + const cu_sparse_matrix_csr_t *A, + const cu_sparse_matrix_csr_t *AT); + + sv_estimator_result_t sv_estimator_run(sv_estimator_ctx_t *ctx, const sv_estimator_opts_t *opts); + + void sv_estimator_free(sv_estimator_ctx_t *ctx); bool cupdlpx_use_spmvop_by_default(void); @@ -123,6 +125,8 @@ extern "C" const restart_parameters_t *restart_params, int termination_evaluation_frequency); + bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol); + void check_termination_criteria(pdhg_solver_state_t *solver_state, const termination_criteria_t *criteria); void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem); @@ -137,7 +141,9 @@ extern "C" void pdhg_final_log(const cupdlpx_result_t *result, const pdhg_parameters_t *params); - void display_iteration_stats(const pdhg_solver_state_t *solver_state, bool verbose); + void display_iteration_header(const pdhg_parameters_t *params); + + void display_iteration_stats(pdhg_solver_state_t *solver_state, const pdhg_parameters_t *params); const char *termination_reason_to_string(termination_reason_t reason); @@ -145,6 +151,8 @@ extern "C" void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm); + void sync_step_sizes_to_gpu(pdhg_solver_state_t *state); + void compute_infeasibility_information(pdhg_solver_state_t *state); void fill_or_copy(double **dest, int n, const double *src, double fill_value); diff --git a/python/README.md b/python/README.md index 0a84ed1..7486b2b 100644 --- a/python/README.md +++ b/python/README.md @@ -179,6 +179,18 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `Presolve`| `presolve` | bool | `True` | Whether to use presolve. | | `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.| | `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual. | +| `Debug` | `debug` | bool | `False` | Developer diagnostics (implies `OutputFlag`). | +| `ActiveSetBoost` | `active_set_boost` | bool | `True` | Enable the active-set stepsize boost. | +| `ASBActivationTol` | `asb_activation_tol` | float | `1e-4` | Residual threshold at which the boost activates. | +| `ASBWindowIter` | `asb_window_iter` | int | `10000` | Number of recent iterations used to identify the active set. | +| `ASBSafetyFactor` | `asb_safety_factor` | float | `0.9` | Boosted step = factor / estimated singular value. | +| `ASBMaxReverts` | `asb_max_reverts` | int | `2` | Divergences tolerated before the boost turns off; a diverged step is always reverted. | +| `ASBMinRaiseRatio` | `asb_min_raise_ratio` | float | `1.1` | Minimum ratio for a step increase. | +| `ASBReestimateChangeRatio` | `asb_reestimate_change_ratio` | float | `0.01` | Fraction of the active set that must change before re-estimating. | +| `ASBConstraintTol` | `asb_constraint_tol` | float | `1e-8` | Tolerance for treating a constraint as binding. | +| `ASBVariableTol` | `asb_variable_tol` | float | `1e-8` | Tolerance for treating a variable as at its bound. | +| `ASBDivergenceCeilingRatio` | `asb_divergence_ceiling_ratio` | float | `0.7` | Step ceiling after a revert, relative to the diverged step. | +| `ASBDivergenceMargin` | `asb_divergence_margin` | float | `0.05` | Allowed fixed-point error increase before a revert. | They can be set in multiple ways: diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index fb0cb2f..93a30b8 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -36,6 +36,7 @@ "IterationLimit": "iteration_limit", "OutputFlag": "verbose", "LogToConsole": "verbose", + "Debug": "debug", # termination evaluation cadence "TermCheckFreq": "termination_evaluation_frequency", # tolerances @@ -67,4 +68,16 @@ "Presolve": "presolve", "MatrixZeroTol": "matrix_zero_tol", "InfiniteBound": "infinite_bound", + # active-set step boost + "ActiveSetBoost": "active_set_boost", + "ASBActivationTol": "asb_activation_tol", + "ASBWindowIter": "asb_window_iter", + "ASBSafetyFactor": "asb_safety_factor", + "ASBMaxReverts": "asb_max_reverts", + "ASBMinRaiseRatio": "asb_min_raise_ratio", + "ASBReestimateChangeRatio": "asb_reestimate_change_ratio", + "ASBConstraintTol": "asb_constraint_tol", + "ASBVariableTol": "asb_variable_tol", + "ASBDivergenceCeilingRatio": "asb_divergence_ceiling_ratio", + "ASBDivergenceMargin": "asb_divergence_margin", } diff --git a/python/cupdlpx/_core.py b/python/cupdlpx/_core.py index cdd7c95..9d718f1 100644 --- a/python/cupdlpx/_core.py +++ b/python/cupdlpx/_core.py @@ -14,4 +14,4 @@ """Thin re-export of the compiled cuPDLPx core extension (_cupdlpx_core).""" -from ._cupdlpx_core import solve_once, get_default_params, read_mps \ No newline at end of file +from ._cupdlpx_core import solve_once, get_default_params, validate_params, read_mps \ No newline at end of file diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index dcb7f6c..1439b9f 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -21,7 +21,7 @@ import numpy as np import scipy.sparse as sp -from ._core import solve_once, get_default_params, read_mps +from ._core import solve_once, get_default_params, validate_params, read_mps from . import PDLP # array-like type @@ -33,10 +33,12 @@ _BOOL_PARAMS = frozenset( { "verbose", + "debug", "has_pock_chambolle_alpha", "bound_objective_rescaling", "feasibility_polishing", "presolve", + "active_set_boost", } ) _INT_PARAMS = frozenset( @@ -46,9 +48,10 @@ "geometric_mean_iterations", "l_inf_ruiz_iterations", "sv_max_iter", + "asb_window_iter", + "asb_max_reverts", } ) -_POSITIVE_INT_PARAMS = frozenset({"termination_evaluation_frequency", "sv_max_iter"}) _FLOAT_PARAMS = frozenset( { "eps_optimal_relative", @@ -65,32 +68,25 @@ "sv_tol", "matrix_zero_tol", "infinite_bound", + "asb_activation_tol", + "asb_safety_factor", + "asb_min_raise_ratio", + "asb_reestimate_change_ratio", + "asb_constraint_tol", + "asb_variable_tol", + "asb_divergence_ceiling_ratio", + "asb_divergence_margin", } ) -_POSITIVE_FLOAT_PARAMS = frozenset( - { - "eps_optimal_relative", - "eps_feasible_relative", - "eps_infeasible_relative", - "eps_feas_polish_relative", - "sv_tol", - "infinite_bound", - } -) -_NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol"}) _STRING_PARAMS = frozenset({"optimality_norm"}) - # int params are stored as C int32 on the backend _INT32_MAX = np.iinfo(np.int32).max -# every backend param must fall in one typed set, else it goes unvalidated +# every backend param must fall in one typed set, else it goes uncoerced; value +# ranges are not mirrored here: setParam hands the full dict to validate_params, +# the same C check optimize() runs, so range rules live in one place _CLASSIFIED_PARAMS = _BOOL_PARAMS | _INT_PARAMS | _FLOAT_PARAMS | _STRING_PARAMS -# refinement sets must be subsets of their base type sets -assert _POSITIVE_INT_PARAMS <= _INT_PARAMS -assert _POSITIVE_FLOAT_PARAMS <= _FLOAT_PARAMS -assert _NONNEGATIVE_FLOAT_PARAMS <= _FLOAT_PARAMS - def _as_dense_f64_c(a: ArrayLike) -> np.ndarray: """ Convert input to an owned C-contiguous numpy array of float64. @@ -522,7 +518,10 @@ def _resolve_param_key(self, name: str) -> str: raise KeyError(f"Unknown parameter '{name}'. Valid names: {valid}") return key - def _validate_param_value(self, key: str, value: Any) -> Any: + def _convert_param_value(self, key: str, value: Any) -> Any: + # type conversion only: turn what the user passed (numpy scalars, 0/1, integer-valued + # floats, mixed-case norm names) into the Python type the backend expects; value + # ranges are checked by validate_params, the same C rules optimize() applies if key in _BOOL_PARAMS: # accept a real bool, numpy bool, or an integer 0/1; store a Python bool if isinstance(value, (bool, np.bool_)): @@ -543,12 +542,8 @@ def _validate_param_value(self, key: str, value: Any) -> Any: value = int(value) else: raise TypeError(f"Parameter '{key}' must be an int.") - if key in _POSITIVE_INT_PARAMS and value <= 0: - raise ValueError(f"Parameter '{key}' must be positive.") - if key not in _POSITIVE_INT_PARAMS and value < 0: - raise ValueError(f"Parameter '{key}' must be nonnegative.") - if value > _INT32_MAX: - raise ValueError(f"Parameter '{key}' must not exceed {_INT32_MAX} (int32 range).") + if not -_INT32_MAX - 1 <= value <= _INT32_MAX: + raise ValueError(f"Parameter '{key}' must fit in an int32.") return value if key in _FLOAT_PARAMS: @@ -557,14 +552,7 @@ def _validate_param_value(self, key: str, value: Any) -> Any: raise TypeError(f"Parameter '{key}' must be a number.") if not isinstance(value, (int, float, np.integer, np.floating)): raise TypeError(f"Parameter '{key}' must be a number.") - value = float(value) - if not np.isfinite(value): - raise ValueError(f"Parameter '{key}' must be finite.") - if key in _POSITIVE_FLOAT_PARAMS and value <= 0.0: - raise ValueError(f"Parameter '{key}' must be positive.") - if key in _NONNEGATIVE_FLOAT_PARAMS and value < 0.0: - raise ValueError(f"Parameter '{key}' must be nonnegative.") - return value + return float(value) if key in _STRING_PARAMS: if not isinstance(value, str): @@ -580,9 +568,12 @@ def setParam(self, name: str, value: Any) -> None: """ Set the value of a solver parameter by name. """ - # resolve name and store + # resolve name, convert the type, validate the resulting parameter set, then store key = self._resolve_param_key(name) - self._params[key] = self._validate_param_value(key, value) + candidate = dict(self._params) + candidate[key] = self._convert_param_value(key, value) + validate_params(candidate) + self._params = candidate def getParam(self, name: str) -> Any: """ @@ -596,11 +587,12 @@ def setParams(self, /, **kwargs) -> None: """ Set multiple solver parameters by name. """ - updates = {} + candidate = dict(self._params) for k, v in kwargs.items(): key = self._resolve_param_key(k) - updates[key] = self._validate_param_value(key, v) - self._params.update(updates) + candidate[key] = self._convert_param_value(key, v) + validate_params(candidate) + self._params = candidate def resetParams(self) -> None: """ diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index ffc6816..b4d7708 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -14,8 +14,8 @@ See the License for the specific language governing permissions and limitations under the License. */ -#include #include +#include #include #include #include @@ -24,6 +24,7 @@ limitations under the License. #include #include #include +#include #include #include "cupdlpx.h" @@ -239,11 +240,9 @@ static void validate_result_dimensions(const cupdlpx_result_t *res, int expected { if (res->num_variables != expected_n || res->num_constraints != expected_m) { - throw std::runtime_error("solve_lp_problem returned result dimensions " + - std::to_string(res->num_variables) + "x" + - std::to_string(res->num_constraints) + - ", expected " + std::to_string(expected_n) + - "x" + std::to_string(expected_m) + "."); + throw std::runtime_error("solve_lp_problem returned result dimensions " + std::to_string(res->num_variables) + + "x" + std::to_string(res->num_constraints) + ", expected " + + std::to_string(expected_n) + "x" + std::to_string(expected_m) + "."); } if (expected_n > 0 && !res->primal_solution) { @@ -268,6 +267,7 @@ static py::dict get_default_params_py() // verbosity d["verbose"] = p.verbose; + d["debug"] = p.debug; d["termination_evaluation_frequency"] = p.termination_evaluation_frequency; // tolerances @@ -311,6 +311,19 @@ static py::dict get_default_params_py() d["matrix_zero_tol"] = p.matrix_zero_tol; d["infinite_bound"] = p.infinite_bound; + // active-set step boost + d["active_set_boost"] = p.active_set_boost; + d["asb_activation_tol"] = p.asb_activation_tol; + d["asb_window_iter"] = p.asb_window_iter; + d["asb_safety_factor"] = p.asb_safety_factor; + d["asb_max_reverts"] = p.asb_max_reverts; + d["asb_min_raise_ratio"] = p.asb_min_raise_ratio; + d["asb_reestimate_change_ratio"] = p.asb_reestimate_change_ratio; + d["asb_constraint_tol"] = p.asb_constraint_tol; + d["asb_variable_tol"] = p.asb_variable_tol; + d["asb_divergence_ceiling_ratio"] = p.asb_divergence_ceiling_ratio; + d["asb_divergence_margin"] = p.asb_divergence_margin; + return d; } @@ -321,8 +334,12 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p return; py::dict d = params_obj.cast(); + // every key the readers below look at; anything else in the dict is a typo + std::unordered_set known_keys; + auto getf = [&](const char *k, double &tgt) { + known_keys.insert(k); if (d.contains(k)) { py::object val = d[k]; @@ -330,7 +347,14 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p { throw std::invalid_argument(std::string(k) + " must be a number."); } - tgt = py::cast(val); + try + { + tgt = py::cast(val); + } + catch (const py::cast_error &) + { + throw std::invalid_argument(std::string(k) + " must be a number."); + } if (!std::isfinite(tgt)) { throw std::invalid_argument(std::string(k) + " must be finite."); @@ -339,6 +363,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p }; auto geti = [&](const char *k, int &tgt) { + known_keys.insert(k); if (d.contains(k)) { py::object val = d[k]; @@ -346,11 +371,19 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p { throw std::invalid_argument(std::string(k) + " must be an int."); } - tgt = py::cast(val); + try + { + tgt = py::cast(val); + } + catch (const py::cast_error &) + { + throw std::invalid_argument(std::string(k) + " must fit in an int32."); + } } }; auto getb = [&](const char *k, bool &tgt) { + known_keys.insert(k); if (d.contains(k)) { py::object val = d[k]; @@ -363,6 +396,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p }; auto get_norm = [&](const char *k, norm_type_t &tgt) { + known_keys.insert(k); if (d.contains(k)) { py::object val = d[k]; @@ -380,6 +414,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p // verbosity getb("verbose", p->verbose); + getb("debug", p->debug); geti("termination_evaluation_frequency", p->termination_evaluation_frequency); // tolerances @@ -423,32 +458,33 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getf("matrix_zero_tol", p->matrix_zero_tol); getf("infinite_bound", p->infinite_bound); - if (p->termination_evaluation_frequency <= 0) - throw std::invalid_argument("termination_evaluation_frequency must be positive."); - if (p->termination_criteria.iteration_limit < 0) - throw std::invalid_argument("iteration_limit must be nonnegative."); - if (p->geometric_mean_iterations < 0) - throw std::invalid_argument("geometric_mean_iterations must be nonnegative."); - if (p->l_inf_ruiz_iterations < 0) - throw std::invalid_argument("l_inf_ruiz_iterations must be nonnegative."); - if (p->sv_max_iter <= 0) - throw std::invalid_argument("sv_max_iter must be positive."); - if (p->termination_criteria.eps_optimal_relative <= 0.0) - throw std::invalid_argument("eps_optimal_relative must be positive."); - if (p->termination_criteria.eps_feasible_relative <= 0.0) - throw std::invalid_argument("eps_feasible_relative must be positive."); - if (p->termination_criteria.eps_infeasible_relative <= 0.0) - throw std::invalid_argument("eps_infeasible_relative must be positive."); - if (p->termination_criteria.eps_feas_polish_relative <= 0.0) - throw std::invalid_argument("eps_feas_polish_relative must be positive."); - if (p->sv_tol <= 0.0) - throw std::invalid_argument("sv_tol must be positive."); - if (p->termination_criteria.time_sec_limit < 0.0) - throw std::invalid_argument("time_sec_limit must be nonnegative."); - if (p->infinite_bound <= 0.0) - throw std::invalid_argument("infinite_bound must be positive."); - if (p->matrix_zero_tol < 0.0) - throw std::invalid_argument("matrix_zero_tol must be nonnegative."); + // active-set step boost + getb("active_set_boost", p->active_set_boost); + getf("asb_activation_tol", p->asb_activation_tol); + geti("asb_window_iter", p->asb_window_iter); + getf("asb_safety_factor", p->asb_safety_factor); + geti("asb_max_reverts", p->asb_max_reverts); + getf("asb_min_raise_ratio", p->asb_min_raise_ratio); + getf("asb_reestimate_change_ratio", p->asb_reestimate_change_ratio); + getf("asb_constraint_tol", p->asb_constraint_tol); + getf("asb_variable_tol", p->asb_variable_tol); + getf("asb_divergence_ceiling_ratio", p->asb_divergence_ceiling_ratio); + getf("asb_divergence_margin", p->asb_divergence_margin); + + for (auto item : d) + { + std::string key = py::str(item.first); + if (!known_keys.count(key)) + { + throw std::invalid_argument("Unknown parameter '" + key + "'."); + } + } + + char validation_error[256]; + if (cupdlpx_validate_parameters(p, validation_error, sizeof(validation_error)) != 0) + { + throw std::invalid_argument(validation_error); + } } // throw if a 1D array's length differs from the expected value @@ -461,8 +497,8 @@ static void expect_len(py::object obj, py::ssize_t expected, const char *name) } if (arr.size() != expected) { - throw std::invalid_argument(std::string(name) + " has wrong length: expected " + - std::to_string(expected) + ", got " + std::to_string((long long)arr.size())); + throw std::invalid_argument(std::string(name) + " has wrong length: expected " + std::to_string(expected) + + ", got " + std::to_string((long long)arr.size())); } } @@ -512,8 +548,8 @@ static void validate_bounds(const double *lower, const double *upper, int size, } // validate a compressed (CSR/CSC) index structure -static void validate_compressed(const int32_t *indptr, const int32_t *indices, int major, int minor, int nnz, - const char *fmt) +static void +validate_compressed(const int32_t *indptr, const int32_t *indices, int major, int minor, int nnz, const char *fmt) { if (indptr[0] != 0) { @@ -809,6 +845,11 @@ static py::dict solve_once(py::object A, info["PrimalRayLinObj"] = res->primal_ray_linear_objective; info["DualRayObj"] = res->dual_ray_objective; + // active-set step boost statistics + info["ASBRaiseCount"] = res->asb_raise_count; + info["ASBRevertCount"] = res->asb_revert_count; + info["ASBPowerIterations"] = res->asb_pi_iterations; + // res freed by res_guard on return return info; } @@ -870,12 +911,26 @@ static py::dict read_mps_py(const std::string &filename) } // module +// Validate a params dict the way solve_once will: defaults overlaid with the dict, +// then cupdlpx_validate_parameters. Raises ValueError on the first violation. +static void validate_params_py(py::object params_obj) +{ + pdhg_parameters_t p; + set_default_parameters(&p); + parse_params_from_python(params_obj, &p); +} + PYBIND11_MODULE(_cupdlpx_core, m) { m.doc() = "cupdlpx core bindings (auto-detect dense/CSR/CSC/COO; initialize default params here)"; m.def("get_default_params", &get_default_params_py, "Return default PDHG parameters as a dict"); + m.def("validate_params", + &validate_params_py, + py::arg("params"), + "Validate a params dict against the solver's parameter rules; raises ValueError"); + m.def("read_mps", &read_mps_py, py::arg("filename"), diff --git a/src/active_set_boost.cu b/src/active_set_boost.cu new file mode 100644 index 0000000..c5c810e --- /dev/null +++ b/src/active_set_boost.cu @@ -0,0 +1,458 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "active_set_boost.h" +#include "utils.h" +#include +#include + +static inline bool asb_raise_possible(double target, double step, double min_raise_ratio) +{ + return target > step && target >= min_raise_ratio * step; +} + +__global__ void asb_var_window_kernel(const double *__restrict__ lb, + const double *__restrict__ ub, + const double *__restrict__ dual_slack, + const double *__restrict__ objective, + double variable_tol, + int *__restrict__ last_free, + bool *__restrict__ mask, + int now, + int window_start, + int *__restrict__ delta_count, + int n) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) + return; + /* The sign of dual_slack = (xbar - projection_input) / step certifies clamping. */ + double tol = variable_tol * fmax(1.0, fabs(objective[i])); + bool fixed = lb[i] == ub[i]; + bool lower_projected = isfinite(lb[i]) && dual_slack[i] > tol; + bool upper_projected = isfinite(ub[i]) && dual_slack[i] < -tol; + bool confidently_clamped = fixed || lower_projected || upper_projected; + if (!confidently_clamped) + last_free[i] = now; + bool m = last_free[i] >= window_start; + if (m != mask[i]) + { + mask[i] = m; + atomicAdd(&delta_count[m ? 0 : 1], 1); + } +} + +__global__ void asb_row_window_kernel(const double *__restrict__ dual_projection_input, + const double *__restrict__ lb, + const double *__restrict__ ub, + double constraint_tol, + int *__restrict__ last_binding, + bool *__restrict__ mask, + int now, + int window_start, + int *__restrict__ delta_count, + int m_rows) +{ + int j = blockIdx.x * blockDim.x + threadIdx.x; + if (j >= m_rows) + return; + /* The row is inactive exactly when the dual projection input q lies safely + inside [-ub,-lb]; a margin retains boundary and uncertain rows. */ + double q = dual_projection_input[j]; + bool binding = lb[j] == ub[j]; + if (!binding) + { + bool safely_inside = true; + if (isfinite(ub[j])) + { + double margin = constraint_tol * fmax(1.0, fmax(fabs(q), fabs(ub[j]))); + safely_inside = safely_inside && (q + ub[j] > margin); + } + if (isfinite(lb[j])) + { + double margin = constraint_tol * fmax(1.0, fmax(fabs(q), fabs(lb[j]))); + safely_inside = safely_inside && (-lb[j] - q > margin); + } + binding = !safely_inside; + } + if (binding) + last_binding[j] = now; + bool m = last_binding[j] >= window_start; + if (m != mask[j]) + { + mask[j] = m; + atomicAdd(&delta_count[m ? 2 : 3], 1); + } +} + +void active_set_boost_init(pdhg_solver_state_t *state) +{ + state->asb_step_ceiling = INFINITY; + state->asb_anchor_best_pd_residual_gap = INFINITY; + size_t n = (size_t)state->num_variables; + size_t m = (size_t)state->num_constraints; + CUDA_CHECK(cudaMalloc(&state->d_asb_var_last_free, n * sizeof(int))); + CUDA_CHECK(cudaMalloc(&state->d_asb_row_last_binding, m * sizeof(int))); + CUDA_CHECK(cudaMalloc(&state->d_asb_col_mask, n * sizeof(bool))); + CUDA_CHECK(cudaMalloc(&state->d_asb_row_mask, m * sizeof(bool))); + CUDA_CHECK(cudaMalloc(&state->d_asb_primal_anchor, n * sizeof(double))); + CUDA_CHECK(cudaMalloc(&state->d_asb_dual_anchor, m * sizeof(double))); + CUDA_CHECK(cudaMalloc(&state->d_asb_dual_slack_anchor, n * sizeof(double))); + CUDA_CHECK(cudaMalloc(&state->d_asb_dual_projection_input, m * sizeof(double))); + CUDA_CHECK(cudaMalloc(&state->d_asb_count, 4 * sizeof(int))); + CUDA_CHECK(cudaMemset(state->d_asb_var_last_free, 0xff, n * sizeof(int))); /* -1 */ + CUDA_CHECK(cudaMemset(state->d_asb_row_last_binding, 0xff, m * sizeof(int))); /* -1 */ + CUDA_CHECK(cudaMemset(state->d_asb_col_mask, 0, n * sizeof(bool))); + CUDA_CHECK(cudaMemset(state->d_asb_row_mask, 0, m * sizeof(bool))); +} + +void active_set_boost_free(pdhg_solver_state_t *state) +{ + if (state->asb_sv_ctx) + { + sv_estimator_free(state->asb_sv_ctx); + state->asb_sv_ctx = NULL; + } + if (state->d_asb_var_last_free) + cudaFree(state->d_asb_var_last_free); + if (state->d_asb_row_last_binding) + cudaFree(state->d_asb_row_last_binding); + if (state->d_asb_col_mask) + cudaFree(state->d_asb_col_mask); + if (state->d_asb_row_mask) + cudaFree(state->d_asb_row_mask); + if (state->d_asb_primal_anchor) + cudaFree(state->d_asb_primal_anchor); + if (state->d_asb_dual_anchor) + cudaFree(state->d_asb_dual_anchor); + if (state->d_asb_dual_slack_anchor) + cudaFree(state->d_asb_dual_slack_anchor); + if (state->d_asb_dual_projection_input) + cudaFree(state->d_asb_dual_projection_input); + if (state->d_asb_count) + cudaFree(state->d_asb_count); +} + +static void asb_restore_anchor(pdhg_solver_state_t *state) +{ + double *p_dst[3] = {state->initial_primal_solution, state->current_primal_solution, state->pdhg_primal_solution}; + double *d_dst[3] = {state->initial_dual_solution, state->current_dual_solution, state->pdhg_dual_solution}; + for (int k = 0; k < 3; ++k) + { + CUDA_CHECK(cudaMemcpyAsync(p_dst[k], + state->d_asb_primal_anchor, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaMemcpyAsync(d_dst[k], + state->d_asb_dual_anchor, + state->num_constraints * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + } + CUDA_CHECK(cudaMemcpyAsync(state->dual_slack, + state->d_asb_dual_slack_anchor, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaStreamSynchronize(state->stream)); +} + +static void +asb_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params, int *delta_add, int *delta_remove) +{ + int now = state->total_count; + int window_start = now - params->asb_window_iter; + CUDA_CHECK(cudaMemsetAsync(state->d_asb_count, 0, 4 * sizeof(int), state->stream)); + asb_var_window_kernel<<num_blocks_primal, THREADS_PER_BLOCK, 0, state->stream>>>( + state->variable_lower_bound, + state->variable_upper_bound, + state->dual_slack, + state->objective_vector, + params->asb_variable_tol, + state->d_asb_var_last_free, + state->d_asb_col_mask, + now, + window_start, + state->d_asb_count, + state->num_variables); + asb_row_window_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( + state->d_asb_dual_projection_input, + state->constraint_lower_bound, + state->constraint_upper_bound, + params->asb_constraint_tol, + state->d_asb_row_last_binding, + state->d_asb_row_mask, + now, + window_start, + state->d_asb_count, + state->num_constraints); + int counts[4] = {0, 0, 0, 0}; + CUDA_CHECK(cudaMemcpyAsync(counts, state->d_asb_count, 4 * sizeof(int), cudaMemcpyDeviceToHost, state->stream)); + CUDA_CHECK(cudaStreamSynchronize(state->stream)); + state->asb_free_variables += counts[0] - counts[1]; + state->asb_binding_constraints += counts[2] - counts[3]; + *delta_add = counts[0] + counts[2]; + *delta_remove = counts[1] + counts[3]; +} + +typedef enum +{ + ASB_SV_OK = 0, /* converged; current estimate stored in state->asb_sv */ + ASB_SV_NO_RAISE = 1, /* early PI termination proved that no step increase is possible */ + ASB_SV_FAILED = 2, /* did not converge or produced an invalid result */ +} asb_sv_status_t; + +static double asb_pi_abort_threshold(const pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + return params->asb_safety_factor / (params->asb_min_raise_ratio * state->step_size); +} + +static asb_sv_status_t +asb_estimate_max_singular_value(pdhg_solver_state_t *state, const pdhg_parameters_t *params, int *out_pi_iterations) +{ + if (!state->asb_sv_ctx) + { + state->asb_sv_ctx = sv_estimator_create( + state->sparse_handle, state->blas_handle, state->constraint_matrix, state->constraint_matrix_t); + } + sv_estimator_opts_t opts = {}; + opts.max_iterations = params->sv_max_iter; + opts.tolerance = params->sv_tol; + opts.d_row_mask = state->d_asb_row_mask; + opts.d_col_mask = state->d_asb_col_mask; + opts.abort_singular_value_threshold = asb_pi_abort_threshold(state, params); + sv_estimator_result_t r = sv_estimator_run(state->asb_sv_ctx, &opts); + state->asb_pi_iterations += r.iterations; + *out_pi_iterations = r.iterations; + if (r.status == SV_ESTIMATOR_ABORTED) + { + state->asb_pi_early_exit_count++; + if (params->debug) + { + printf("[active-set boost] iter %d: sv estimate aborted early (PI %d iters, lower bound %.3e)\n", + state->total_count, + r.iterations, + r.max_singular_value); + } + return ASB_SV_NO_RAISE; + } + if (r.status != SV_ESTIMATOR_CONVERGED) + { + if (state->step_size <= state->base_step_size) + { + state->asb_sv = 0.0; + } + return ASB_SV_FAILED; + } + state->asb_sv = r.max_singular_value; + return ASB_SV_OK; +} + +static asb_action_t asb_revert(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + state->restart_count++; + double diverged_step = state->step_size; + state->asb_step_ceiling = params->asb_divergence_ceiling_ratio * diverged_step; + asb_restore_anchor(state); + compute_residual(state, params->optimality_norm); + /* The infeasibility rays describe the discarded iterate; the anchor has none. */ + state->max_primal_ray_infeasibility = 0.0; + state->max_dual_ray_infeasibility = 0.0; + state->primal_ray_linear_objective = 0.0; + state->dual_ray_objective = 0.0; + state->step_size = state->base_step_size; + state->primal_weight = state->asb_anchor_primal_weight; + state->primal_weight_error_sum = state->asb_anchor_pw_error_sum; + state->primal_weight_last_error = state->asb_anchor_pw_last_error; + state->best_primal_weight = state->asb_anchor_best_pw; + state->best_primal_dual_residual_gap = state->asb_anchor_best_pd_residual_gap; + state->inner_count = 0; + state->last_trial_fixed_point_error = INFINITY; + state->asb_sv = 0.0; + state->asb_no_raise_certified = false; + sync_step_sizes_to_gpu(state); + state->asb_revert_count++; + if (state->asb_revert_count >= params->asb_max_reverts) + { + state->asb_phase = ASB_PHASE_OFF; + } + if (params->debug) + { + printf("[active-set boost] iter %d: divergence #%d, stepsize reverted " + "(fixed-point error %.2e vs span initial %.2e, new stepsize ceiling %.3e)%s\n", + state->total_count, + state->asb_revert_count, + state->fixed_point_error, + state->initial_fixed_point_error, + state->asb_step_ceiling, + state->asb_revert_count >= params->asb_max_reverts ? ", controller off" : ""); + } + else if (params->verbose) + { + printf("[active-set boost] iter %d: divergence #%d, stepsize reverted%s\n", + state->total_count, + state->asb_revert_count, + state->asb_revert_count >= params->asb_max_reverts ? ", controller off" : ""); + } + return ASB_ACTION_REVERT; +} + +static double asb_worst_residual(const pdhg_solver_state_t *state) +{ + return fmax(state->relative_primal_residual, fmax(state->relative_dual_residual, state->relative_objective_gap)); +} + +static double asb_compute_target(const pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + double global_step = state->base_step_size; + double target = (state->asb_sv > 0.0) ? fmax(global_step, params->asb_safety_factor / state->asb_sv) : global_step; + return fmax(global_step, fmin(target, state->asb_step_ceiling)); +} + +static bool asb_change_ready(const pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + return state->asb_changes_since_estimate > 0 && + (double)state->asb_changes_since_estimate >= + params->asb_reestimate_change_ratio * (double)(state->asb_free_variables + state->asb_binding_constraints); +} + +static double asb_refresh_target(pdhg_solver_state_t *state, const pdhg_parameters_t *params, int *out_pi_iterations) +{ + if (!((state->asb_sv <= 0.0 && !state->asb_no_raise_certified) || asb_change_ready(state, params))) + return state->step_size; + /* The ceiling already precludes an increase: certify without running the power iteration. */ + if (!asb_raise_possible( + fmax(state->base_step_size, state->asb_step_ceiling), state->step_size, params->asb_min_raise_ratio)) + { + state->asb_no_raise_certified = true; + state->asb_changes_since_estimate = 0; + return state->step_size; + } + asb_sv_status_t status = asb_estimate_max_singular_value(state, params, out_pi_iterations); + if (status == ASB_SV_OK) + { + state->asb_changes_since_estimate = 0; + state->asb_no_raise_certified = false; + return asb_compute_target(state, params); + } + if (status == ASB_SV_FAILED) + { + state->asb_sv_failed_count++; + } + else + { + state->asb_no_raise_certified = true; + state->asb_changes_since_estimate = 0; + } + return state->step_size; +} + +void active_set_boost_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + if (state->asb_phase == ASB_PHASE_OFF) + return; + int delta_add = 0, delta_remove = 0; + asb_update_window(state, params, &delta_add, &delta_remove); + state->asb_changes_since_estimate += delta_add + delta_remove; +} + +asb_action_t active_set_boost_check(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + if (state->asb_phase == ASB_PHASE_OFF) + return ASB_ACTION_NONE; + /* A stop certified by the residuals keeps its iterate; only a limit stop may still revert. */ + if (state->termination_reason != TERMINATION_REASON_UNSPECIFIED && + state->termination_reason != TERMINATION_REASON_ITERATION_LIMIT && + state->termination_reason != TERMINATION_REASON_TIME_LIMIT) + return ASB_ACTION_NONE; + + bool residuals_finite = isfinite(state->relative_primal_residual) && isfinite(state->relative_dual_residual) && + isfinite(state->relative_objective_gap); + double res = asb_worst_residual(state); + + if (state->asb_phase == ASB_PHASE_WAITING) + { + if (!residuals_finite || res >= params->asb_activation_tol) + return ASB_ACTION_NONE; + if (state->asb_free_variables == state->num_variables && + state->asb_binding_constraints == state->num_constraints) + return ASB_ACTION_NONE; + state->asb_phase = ASB_PHASE_ACTIVE; + } + + bool boosted = state->step_size > state->base_step_size; + double fpe = state->fixed_point_error; + bool tripped = boosted && + (!residuals_finite || !isfinite(fpe) || + (state->initial_fixed_point_error > 0.0 && + fpe > (1.0 + params->asb_divergence_margin) * state->initial_fixed_point_error)); + if (tripped) + { + return asb_revert(state, params); + } + + return ASB_ACTION_NONE; +} + +static void asb_snapshot_anchor(pdhg_solver_state_t *state) +{ + CUDA_CHECK(cudaMemcpyAsync(state->d_asb_primal_anchor, + state->initial_primal_solution, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaMemcpyAsync(state->d_asb_dual_anchor, + state->initial_dual_solution, + state->num_constraints * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + /* dual_slack describes the iterate that became the anchor; compute_residual + needs it to reproduce the anchor's residuals after a revert */ + CUDA_CHECK(cudaMemcpyAsync(state->d_asb_dual_slack_anchor, + state->dual_slack, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + state->asb_anchor_primal_weight = state->primal_weight; + state->asb_anchor_pw_error_sum = state->primal_weight_error_sum; + state->asb_anchor_pw_last_error = state->primal_weight_last_error; + state->asb_anchor_best_pw = state->best_primal_weight; + state->asb_anchor_best_pd_residual_gap = state->best_primal_dual_residual_gap; +} + +void active_set_boost_on_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +{ + if (state->asb_phase != ASB_PHASE_ACTIVE) + return; + int pi_iterations = 0; + double target = asb_refresh_target(state, params, &pi_iterations); + if (asb_raise_possible(target, state->step_size, params->asb_min_raise_ratio)) + { + if (params->debug) + { + printf("[active-set boost] iter %d: stepsize %.3e -> %.3e (sv %.3e, PI iter %d)\n", + state->total_count, + state->step_size, + target, + state->asb_sv, + pi_iterations); + } + state->step_size = target; + state->asb_raise_count++; + asb_snapshot_anchor(state); + } +} diff --git a/src/cli.c b/src/cli.c index 999c812..cc85c93 100644 --- a/src/cli.c +++ b/src/cli.c @@ -125,6 +125,9 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir, fprintf(outfile, "Precondition time (sec): %e\n", result->rescaling_time_sec); fprintf(outfile, "Runtime (sec): %e\n", result->cumulative_time_sec); fprintf(outfile, "Iterations Count: %d\n", result->total_count); + fprintf(outfile, "ASB Step Raise Count: %d\n", result->asb_raise_count); + fprintf(outfile, "ASB Revert Count: %d\n", result->asb_revert_count); + fprintf(outfile, "ASB Power Iterations: %d\n", result->asb_pi_iterations); fprintf(outfile, "Primal Objective Value: %e\n", result->primal_objective_value); fprintf(outfile, "Dual Objective Value: %e\n", result->dual_objective_value); fprintf(outfile, "Relative Primal Residual: %e\n", result->relative_primal_residual); @@ -163,6 +166,9 @@ void print_usage(const char *prog_name) fprintf(stderr, " -v, --verbose " "Enable verbose logging (enabled by default; kept for compatibility).\n"); + fprintf(stderr, + " --debug " + "Developer diagnostics (implies verbose).\n"); fprintf(stderr, " -q, --quiet " "Disable verbose logging.\n"); @@ -218,6 +224,39 @@ void print_usage(const char *prog_name) fprintf(stderr, " --infinite_bound . " "Bounds at or beyond this are treated as infinite (default: 1e20).\n"); + fprintf(stderr, + " --no_active_set_boost " + "Disable the active-set stepsize boost (default: enabled).\n"); + fprintf(stderr, + " --asb_activation_tol " + "Residual threshold at which the boost activates (default: 1e-4).\n"); + fprintf(stderr, + " --asb_window_iter " + "Number of recent iterations used to identify the active set (default: 10000).\n"); + fprintf(stderr, + " --asb_safety_factor " + "Boosted step = factor / estimated singular value (default: 0.9).\n"); + fprintf(stderr, + " --asb_max_reverts " + "Divergences tolerated before the boost turns off; a diverged step is always reverted (default: 2).\n"); + fprintf(stderr, + " --asb_min_raise_ratio " + "Minimum ratio for a step increase (default: 1.1).\n"); + fprintf(stderr, + " --asb_reestimate_change_ratio " + "Fraction of the active set that must change before re-estimating (default: 0.01).\n"); + fprintf(stderr, + " --asb_constraint_tol " + "Tolerance for treating a constraint as binding (default: 1e-8).\n"); + fprintf(stderr, + " --asb_variable_tol " + "Tolerance for treating a variable as at its bound (default: 1e-8).\n"); + fprintf(stderr, + " --asb_divergence_ceiling_ratio " + "Step ceiling after a revert, relative to the diverged step (default: 0.7).\n"); + fprintf(stderr, + " --asb_divergence_margin " + "Allowed fixed-point error increase before a revert (default: 0.05).\n"); } int main(int argc, char *argv[]) @@ -246,6 +285,18 @@ int main(int argc, char *argv[]) {"no_presolve", no_argument, 0, 1015}, {"matrix_zero_tol", required_argument, 0, 1016}, {"infinite_bound", required_argument, 0, 1017}, + {"no_active_set_boost", no_argument, 0, 1019}, + {"asb_activation_tol", required_argument, 0, 1029}, + {"debug", no_argument, 0, 1030}, + {"asb_window_iter", required_argument, 0, 1020}, + {"asb_safety_factor", required_argument, 0, 1021}, + {"asb_max_reverts", required_argument, 0, 1022}, + {"asb_min_raise_ratio", required_argument, 0, 1023}, + {"asb_reestimate_change_ratio", required_argument, 0, 1024}, + {"asb_constraint_tol", required_argument, 0, 1025}, + {"asb_variable_tol", required_argument, 0, 1026}, + {"asb_divergence_ceiling_ratio", required_argument, 0, 1027}, + {"asb_divergence_margin", required_argument, 0, 1028}, {0, 0, 0, 0}}; int opt; @@ -331,6 +382,42 @@ int main(int argc, char *argv[]) case 1017: // --infinite_bound params.infinite_bound = atof(optarg); break; + case 1019: // --no_active_set_boost + params.active_set_boost = false; + break; + case 1020: // --asb_window_iter + params.asb_window_iter = atoi(optarg); + break; + case 1021: // --asb_safety_factor + params.asb_safety_factor = atof(optarg); + break; + case 1022: // --asb_max_reverts + params.asb_max_reverts = atoi(optarg); + break; + case 1023: // --asb_min_raise_ratio + params.asb_min_raise_ratio = atof(optarg); + break; + case 1024: // --asb_reestimate_change_ratio + params.asb_reestimate_change_ratio = atof(optarg); + break; + case 1025: // --asb_constraint_tol + params.asb_constraint_tol = atof(optarg); + break; + case 1026: // --asb_variable_tol + params.asb_variable_tol = atof(optarg); + break; + case 1027: // --asb_divergence_ceiling_ratio + params.asb_divergence_ceiling_ratio = atof(optarg); + break; + case 1028: // --asb_divergence_margin + params.asb_divergence_margin = atof(optarg); + break; + case 1029: // --asb_activation_tol + params.asb_activation_tol = atof(optarg); + break; + case 1030: // --debug + params.debug = true; + break; case '?': // Unknown option return 1; } diff --git a/src/cupdlpx.c b/src/cupdlpx.c index 159b291..10622e9 100644 --- a/src/cupdlpx.c +++ b/src/cupdlpx.c @@ -179,6 +179,85 @@ void set_start_values(lp_problem_t *prob, const double *primal, const double *du } } +/* The single authority on parameter ranges. Front ends only parse and convert + types; every semantic check lives here. */ +#define CUPDLPX_CHECK_PARAM(cond, ...) \ + do \ + { \ + if (!(cond)) \ + { \ + if (error_message && error_message_size > 0) \ + { \ + snprintf(error_message, error_message_size, __VA_ARGS__); \ + } \ + return 1; \ + } \ + } while (0) + +int cupdlpx_validate_parameters(const pdhg_parameters_t *p, char *error_message, size_t error_message_size) +{ + CUPDLPX_CHECK_PARAM(p != NULL, "params must not be NULL"); + CUPDLPX_CHECK_PARAM(p->termination_evaluation_frequency >= 3, + "termination_evaluation_frequency must be >= 3 (got %d)", + p->termination_evaluation_frequency); + CUPDLPX_CHECK_PARAM(p->termination_criteria.iteration_limit >= 0, + "iteration_limit must be nonnegative (got %d)", + p->termination_criteria.iteration_limit); + CUPDLPX_CHECK_PARAM(p->geometric_mean_iterations >= 0, + "geometric_mean_iterations must be nonnegative (got %d)", + p->geometric_mean_iterations); + CUPDLPX_CHECK_PARAM( + p->l_inf_ruiz_iterations >= 0, "l_inf_ruiz_iterations must be nonnegative (got %d)", p->l_inf_ruiz_iterations); + CUPDLPX_CHECK_PARAM(p->sv_max_iter > 0, "sv_max_iter must be positive (got %d)", p->sv_max_iter); + CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_optimal_relative > 0.0, + "eps_optimal_relative must be positive (got %g)", + p->termination_criteria.eps_optimal_relative); + CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_feasible_relative > 0.0, + "eps_feasible_relative must be positive (got %g)", + p->termination_criteria.eps_feasible_relative); + CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_infeasible_relative > 0.0, + "eps_infeasible_relative must be positive (got %g)", + p->termination_criteria.eps_infeasible_relative); + CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_feas_polish_relative > 0.0, + "eps_feas_polish_relative must be positive (got %g)", + p->termination_criteria.eps_feas_polish_relative); + CUPDLPX_CHECK_PARAM(p->sv_tol > 0.0, "sv_tol must be positive (got %g)", p->sv_tol); + CUPDLPX_CHECK_PARAM(p->termination_criteria.time_sec_limit >= 0.0, + "time_sec_limit must be nonnegative (got %g)", + p->termination_criteria.time_sec_limit); + CUPDLPX_CHECK_PARAM(p->infinite_bound > 0.0, "infinite_bound must be positive (got %g)", p->infinite_bound); + CUPDLPX_CHECK_PARAM(p->matrix_zero_tol >= 0.0, "matrix_zero_tol must be nonnegative (got %g)", p->matrix_zero_tol); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_activation_tol) && p->asb_activation_tol >= 0.0, + "asb_activation_tol must be finite and >= 0 (got %g)", + p->asb_activation_tol); + CUPDLPX_CHECK_PARAM(p->asb_window_iter > 0, "asb_window_iter must be positive (got %d)", p->asb_window_iter); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_safety_factor) && p->asb_safety_factor > 0.0, + "asb_safety_factor must be finite and positive (got %g)", + p->asb_safety_factor); + CUPDLPX_CHECK_PARAM(p->asb_max_reverts >= 0, "asb_max_reverts must be nonnegative (got %d)", p->asb_max_reverts); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_min_raise_ratio) && p->asb_min_raise_ratio >= 1.0, + "asb_min_raise_ratio must be finite and >= 1 (got %g)", + p->asb_min_raise_ratio); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_reestimate_change_ratio) && p->asb_reestimate_change_ratio >= 0.0, + "asb_reestimate_change_ratio must be finite and >= 0 (got %g)", + p->asb_reestimate_change_ratio); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_constraint_tol) && p->asb_constraint_tol >= 0.0, + "asb_constraint_tol must be finite and >= 0 (got %g)", + p->asb_constraint_tol); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_variable_tol) && p->asb_variable_tol >= 0.0, + "asb_variable_tol must be finite and >= 0 (got %g)", + p->asb_variable_tol); + CUPDLPX_CHECK_PARAM(p->asb_divergence_ceiling_ratio > 0.0 && p->asb_divergence_ceiling_ratio <= 1.0, + "asb_divergence_ceiling_ratio must be in (0, 1] (got %g)", + p->asb_divergence_ceiling_ratio); + CUPDLPX_CHECK_PARAM(isfinite(p->asb_divergence_margin) && p->asb_divergence_margin >= 0.0, + "asb_divergence_margin must be finite and >= 0 (got %g)", + p->asb_divergence_margin); + return 0; +} + +#undef CUPDLPX_CHECK_PARAM + cupdlpx_result_t *solve_lp_problem(lp_problem_t *prob, const pdhg_parameters_t *params) { // argument checks @@ -198,6 +277,10 @@ cupdlpx_result_t *solve_lp_problem(lp_problem_t *prob, const pdhg_parameters_t * { set_default_parameters(&local_params); } + if (local_params.debug) + { + local_params.verbose = true; + } // call optimizer cupdlpx_result_t *res = optimize(&local_params, prob); diff --git a/src/feasibility_polish.cu b/src/feasibility_polish.cu index 0584752..544da5c 100644 --- a/src/feasibility_polish.cu +++ b/src/feasibility_polish.cu @@ -175,10 +175,7 @@ void primal_feasibility_polish(const pdhg_parameters_t *params, state->total_count += params->termination_evaluation_frequency; check_feas_polishing_termination_criteria(state, ori_state, ¶ms->termination_criteria, true); - if (state->total_count % get_print_frequency(state->total_count) == 0) - { - display_feas_polish_iteration_stats(state, params->verbose, true); - } + display_feas_polish_iteration_stats(state, params->verbose, true); // Check Adaptive Restart do_restart = @@ -250,10 +247,7 @@ void dual_feasibility_polish(const pdhg_parameters_t *params, state->total_count += params->termination_evaluation_frequency; check_feas_polishing_termination_criteria(state, ori_state, ¶ms->termination_criteria, false); - if (state->total_count % get_print_frequency(state->total_count) == 0) - { - display_feas_polish_iteration_stats(state, params->verbose, false); - } + display_feas_polish_iteration_stats(state, params->verbose, false); // Check Adaptive Restart do_restart = @@ -276,6 +270,7 @@ static pdhg_solver_state_t *initialize_primal_feas_polish_state(const pdhg_solve { pdhg_solver_state_t *primal_state = (pdhg_solver_state_t *)safe_malloc(sizeof(pdhg_solver_state_t)); *primal_state = *original_state; + primal_state->d_asb_dual_projection_input = NULL; int num_var = original_state->num_variables; int num_cons = original_state->num_constraints; @@ -396,6 +391,7 @@ static pdhg_solver_state_t *initialize_dual_feas_polish_state(const pdhg_solver_ { pdhg_solver_state_t *dual_state = (pdhg_solver_state_t *)safe_malloc(sizeof(pdhg_solver_state_t)); *dual_state = *original_state; + dual_state->d_asb_dual_projection_input = NULL; int num_var = original_state->num_variables; int num_cons = original_state->num_constraints; diff --git a/src/solver.cu b/src/solver.cu index 636c7b5..8317eb7 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -14,6 +14,7 @@ See the License for the specific language governing permissions and limitations under the License. */ +#include "active_set_boost.h" #include "cupdlpx.h" #include "feasibility_polish.h" #include "internal_types.h" @@ -111,6 +112,7 @@ __global__ void compute_next_dual_solution_major_kernel(double *__restrict__ cur const double *__restrict__ const_ub, int n, const double *__restrict__ d_step_size, + double *__restrict__ dual_projection_input, const int *__restrict__ d_base_count, int k_offset, double reflection_coeff); @@ -122,7 +124,6 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, const int k_offset, const double reflection_coefficient, bool is_major); -static void sync_step_sizes_to_gpu(pdhg_solver_state_t *state); void sync_inner_count_to_gpu(pdhg_solver_state_t *state); static void check_params_validity(const pdhg_parameters_t *params); @@ -152,10 +153,16 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * } pdhg_solver_state_t *state = initialize_solver_state(working_problem, params, original_problem->objective_sense); - display_iteration_stats(state, params->verbose); - initialize_step_size_and_primal_weight(state, params); sync_step_sizes_to_gpu(state); + if (params->active_set_boost) + { + active_set_boost_init(state); + } + + display_iteration_header(params); + compute_residual(state, params->optimality_norm); + display_iteration_stats(state, params); state->start_time = clock(); bool do_restart = false; @@ -207,26 +214,31 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; - // Logging - if (state->total_count % get_print_frequency(state->total_count) == 0) + if (params->active_set_boost) { - display_iteration_stats(state, params->verbose); + active_set_boost_update_window(state, params); } // Check Termination check_termination_criteria(state, ¶ms->termination_criteria); - if (state->termination_reason != TERMINATION_REASON_UNSPECIFIED) + bool terminated = state->termination_reason != TERMINATION_REASON_UNSPECIFIED; + + if (params->active_set_boost && active_set_boost_check(state, params) == ASB_ACTION_REVERT) { - break; + do_restart = true; } - - // Check Adaptive Restart - do_restart = - should_do_adaptive_restart(state, ¶ms->restart_params, params->termination_evaluation_frequency); - if (do_restart) + else if (!terminated && + should_do_adaptive_restart(state, ¶ms->restart_params, params->termination_evaluation_frequency)) { perform_restart(state, params); - sync_step_sizes_to_gpu(state); + do_restart = true; + } + + display_iteration_stats(state, params); + + if (terminated) + { + break; } } @@ -239,12 +251,14 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * { state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; compute_residual(state, params->optimality_norm); - display_iteration_stats(state, params->verbose); + display_iteration_stats(state, params); } if (params->feasibility_polishing && state->termination_reason != TERMINATION_REASON_DUAL_INFEASIBLE && state->termination_reason != TERMINATION_REASON_PRIMAL_INFEASIBLE) { + state->step_size = state->base_step_size; + sync_step_sizes_to_gpu(state); feasibility_polish(params, state); } @@ -265,7 +279,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * return result; } -static void sync_step_sizes_to_gpu(pdhg_solver_state_t *state) +void sync_step_sizes_to_gpu(pdhg_solver_state_t *state) { double current_primal_step = state->step_size / state->primal_weight; double current_dual_step = state->step_size * state->primal_weight; @@ -284,11 +298,10 @@ void sync_inner_count_to_gpu(pdhg_solver_state_t *state) static void check_params_validity(const pdhg_parameters_t *params) { - if (params->termination_evaluation_frequency < 3) + char error_message[256]; + if (cupdlpx_validate_parameters(params, error_message, sizeof(error_message)) != 0) { - fprintf(stderr, - "Error: termination_evaluation_frequency must be >= 3 (got %d).\n", - params->termination_evaluation_frequency); + fprintf(stderr, "Error: %s.\n", error_message); exit(EXIT_FAILURE); } } @@ -452,6 +465,7 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ ALLOC_ZERO(state->reflected_primal_solution, var_bytes); ALLOC_ZERO(state->dual_product, var_bytes); ALLOC_ZERO(state->dual_slack, var_bytes); + ALLOC_ZERO(state->infeasibility_dual_scratch, var_bytes); ALLOC_ZERO(state->dual_residual, var_bytes); ALLOC_ZERO(state->delta_primal_solution, var_bytes); @@ -633,30 +647,6 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ state->dual_product); free(ones_dual_h); - if (params->verbose) - { - printf("---------------------------------------------------------------------" - "------------------\n"); - printf("%s | %s | %s | %s \n", - " runtime ", - " objective ", - " absolute residuals ", - " relative residuals "); - printf("%s %s | %s %s | %s %s %s | %s %s %s \n", - " iter", - " time ", - " pr obj ", - " du obj ", - " pr res", - " du res", - " gap ", - " pr res", - " du res", - " gap "); - printf("---------------------------------------------------------------------" - "------------------\n"); - } - return state; } @@ -816,6 +806,7 @@ __global__ void compute_next_dual_solution_major_kernel(double *__restrict__ cur const double *__restrict__ const_ub, int n, const double *__restrict__ d_step_size, + double *__restrict__ dual_projection_input, const int *__restrict__ d_base_count, int k_offset, double reflection_coeff) @@ -827,6 +818,8 @@ __global__ void compute_next_dual_solution_major_kernel(double *__restrict__ cur if (i < n) { double temp = current_dual[i] / step_size - primal_product[i]; + if (dual_projection_input) + dual_projection_input[i] = temp; double temp_proj = fmax(-const_ub[i], fmin(temp, -const_lb[i])); pdhg_dual[i] = (temp - temp_proj) * step_size; reflected_dual[i] = 2.0 * pdhg_dual[i] - current_dual[i]; @@ -943,6 +936,7 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, state->constraint_upper_bound, state->num_constraints, state->d_dual_step_size, + state->d_asb_dual_projection_input, state->d_inner_count, k_offset, reflection_coefficient); @@ -965,6 +959,11 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params) { + if (params->debug) + { + printf("[restart] iter %d: %s\n", state->total_count, state->last_restart_reason); + } + state->restart_count++; compute_delta_solution_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( state->initial_primal_solution, state->pdhg_primal_solution, @@ -1028,6 +1027,12 @@ static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t state->inner_count = 0; state->last_trial_fixed_point_error = INFINITY; + + if (params->active_set_boost) + { + active_set_boost_on_restart(state, params); + } + sync_step_sizes_to_gpu(state); } static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, const pdhg_parameters_t *params) @@ -1038,15 +1043,28 @@ static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, c } else { - double max_sv = estimate_maximum_singular_value(state->sparse_handle, - state->blas_handle, - state->constraint_matrix, - state->constraint_matrix_t, - params->sv_max_iter, - params->sv_tol); - state->step_size = 0.998 / max_sv; + sv_estimator_opts_t opts = {}; + opts.max_iterations = params->sv_max_iter; + opts.tolerance = params->sv_tol; + sv_estimator_ctx_t *estimator = sv_estimator_create( + state->sparse_handle, state->blas_handle, state->constraint_matrix, state->constraint_matrix_t); + sv_estimator_result_t sv = sv_estimator_run(estimator, &opts); + sv_estimator_free(estimator); + state->step_size = (sv.status == SV_ESTIMATOR_DEGENERATE) ? 1.0 : 0.998 / sv.max_singular_value; + if (params->verbose) + { + const char *note = sv.status == SV_ESTIMATOR_CONVERGED ? "converged" + : sv.status == SV_ESTIMATOR_DEGENERATE ? "degenerate" + : "not converged"; + printf("\nEstimating step size\n"); + printf(" Power iterations : %d (%s)\n", sv.iterations, note); + printf(" Max singular value : %.4e\n", sv.max_singular_value); + printf(" Step size : %.4e\n", state->step_size); + } } + state->base_step_size = state->step_size; + if (params->bound_objective_rescaling) { state->primal_weight = 1.0; @@ -1091,7 +1109,10 @@ static void compute_fixed_point_error(pdhg_solver_state_t *state) state->delta_primal_solution, 1, &cross_term)); - interaction = 2 * state->step_size * cross_term; + // measure in the M-norm of the smaller (global) step: M is indefinite for a + // controller-boosted step + double norm_step = fmin(state->step_size, state->base_step_size); + interaction = 2 * norm_step * cross_term; state->fixed_point_error = sqrt(movement + interaction); } @@ -1103,6 +1124,10 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) return; } + /* the cached estimator references the solver's handles and matrix storage: + destroy it first */ + active_set_boost_free(state); + if (state->spmv_ctx) cupdlpx_spmv_ctx_destroy(state->spmv_ctx); if (state->sparse_handle) @@ -1178,6 +1203,8 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) CUDA_CHECK(cudaFree(state->primal_slack)); if (state->dual_slack) CUDA_CHECK(cudaFree(state->dual_slack)); + if (state->infeasibility_dual_scratch) + CUDA_CHECK(cudaFree(state->infeasibility_dual_scratch)); if (state->primal_residual) CUDA_CHECK(cudaFree(state->primal_residual)); if (state->dual_residual) @@ -1278,13 +1305,11 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, co results->primal_ray_linear_objective = state->primal_ray_linear_objective; results->dual_ray_objective = state->dual_ray_objective; results->termination_reason = state->termination_reason; + results->asb_raise_count = state->asb_raise_count; + results->asb_revert_count = state->asb_revert_count; + results->asb_pi_iterations = state->asb_pi_iterations; results->feasibility_polishing_time = state->feasibility_polishing_time; results->feasibility_iteration = state->feasibility_iteration; - // if (presolve_stats != NULL) { - // results->presolve_stats = *presolve_stats; - // } else { - // memset(&(results->presolve_stats), 0, sizeof(PresolveStats)); - // } return results; } diff --git a/src/utils.cu b/src/utils.cu index 7ff69ce..413797f 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -17,6 +17,7 @@ limitations under the License. #include "utils.h" #include #include +#include #ifndef CUPDLPX_VERSION #define CUPDLPX_VERSION "unknown" @@ -63,35 +64,65 @@ void *safe_realloc(void *ptr, size_t new_size) return tmp; } -double estimate_maximum_singular_value(cusparseHandle_t sparse_handle, - cublasHandle_t blas_handle, - const cu_sparse_matrix_csr_t *A, - const cu_sparse_matrix_csr_t *AT, - int max_iterations, - double tolerance) +__global__ void elementwise_mask_kernel(double *__restrict__ v, const bool *__restrict__ mask, int n) { - const int m = A->num_rows; - const int n = A->num_cols; - double *eigenvector_d, *next_eigenvector_d, *dual_product_d; - - CUDA_CHECK(cudaMalloc(&eigenvector_d, m * sizeof(double))); - CUDA_CHECK(cudaMalloc(&next_eigenvector_d, m * sizeof(double))); - CUDA_CHECK(cudaMalloc(&dual_product_d, n * sizeof(double))); + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n) + v[i] = mask[i] ? v[i] : 0.0; +} - double *eigenvector_h = (double *)safe_malloc(m * sizeof(double)); - for (int i = 0; i < m; ++i) +__global__ void warm_start_fill_kernel(double *__restrict__ v, const bool *__restrict__ mask, double amplitude, int n) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) + return; + if (mask[i] && v[i] == 0.0) { - eigenvector_h[i] = dist(gen); + /* Deterministic per-index noise (splitmix64 hash to [-0.5, 0.5)) preserves + reproducibility without RNG state or a host-to-device copy. */ + unsigned long long z = (unsigned long long)i + 0x9E3779B97F4A7C15ULL; + z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL; + z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL; + z = z ^ (z >> 31); + double u = (double)(z >> 11) * (1.0 / 9007199254740992.0); /* [0, 1) */ + v[i] = amplitude * (u - 0.5); } +} - CUDA_CHECK(cudaMemcpy(eigenvector_d, eigenvector_h, m * sizeof(double), cudaMemcpyHostToDevice)); - free(eigenvector_h); - - double sigma_max_sq = 1.0; - const double one = 1.0; - +/* Masks are applied elementwise per run, so a cached context stays valid across mask + changes; it must be recreated only if the matrix itself changes. */ +struct sv_estimator_ctx +{ + cusparseHandle_t sparse_handle; + cublasHandle_t blas_handle; + int num_rows, num_cols; + double *eigenvector_d; + double *next_eigenvector_d; + double *dual_product_d; cusparseSpMatDescr_t matA, matAT; - CUSPARSE_CHECK(cusparseCreateCsr(&matA, + cusparseDnVecDescr_t vecEigen, vecNextEigen, vecDual; + void *descrAT, *descrA; + void *planAT, *planA; + void *dBufferAT, *dBufferA; + bool have_warm_start; /* eigenvector_d contains the final vector from the previous estimate */ +}; + +sv_estimator_ctx_t *sv_estimator_create(cusparseHandle_t sparse_handle, + cublasHandle_t blas_handle, + const cu_sparse_matrix_csr_t *A, + const cu_sparse_matrix_csr_t *AT) +{ + sv_estimator_ctx_t *ctx = (sv_estimator_ctx_t *)safe_calloc(1, sizeof(sv_estimator_ctx_t)); + ctx->sparse_handle = sparse_handle; + ctx->blas_handle = blas_handle; + const int m = A->num_rows; + const int n = A->num_cols; + ctx->num_rows = m; + ctx->num_cols = n; + CUDA_CHECK(cudaMalloc(&ctx->eigenvector_d, m * sizeof(double))); + CUDA_CHECK(cudaMalloc(&ctx->next_eigenvector_d, m * sizeof(double))); + CUDA_CHECK(cudaMalloc(&ctx->dual_product_d, n * sizeof(double))); + CUSPARSE_CHECK(cusparseCreateCsr(&ctx->matA, A->num_rows, A->num_cols, A->num_nonzeros, @@ -102,7 +133,7 @@ double estimate_maximum_singular_value(cusparseHandle_t sparse_handle, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F)); - CUSPARSE_CHECK(cusparseCreateCsr(&matAT, + CUSPARSE_CHECK(cusparseCreateCsr(&ctx->matAT, AT->num_rows, AT->num_cols, AT->num_nonzeros, @@ -113,68 +144,162 @@ double estimate_maximum_singular_value(cusparseHandle_t sparse_handle, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F)); - - cusparseDnVecDescr_t vecEigen, vecNextEigen, vecDual; - CUSPARSE_CHECK(cusparseCreateDnVec(&vecEigen, m, eigenvector_d, CUDA_R_64F)); - CUSPARSE_CHECK(cusparseCreateDnVec(&vecNextEigen, m, next_eigenvector_d, CUDA_R_64F)); - CUSPARSE_CHECK(cusparseCreateDnVec(&vecDual, n, dual_product_d, CUDA_R_64F)); - - void *descrAT = NULL; - void *descrA = NULL; - void *planAT = NULL; - void *planA = NULL; - - void *dBufferAT = NULL; - void *dBufferA = NULL; + CUSPARSE_CHECK(cusparseCreateDnVec(&ctx->vecEigen, m, ctx->eigenvector_d, CUDA_R_64F)); + CUSPARSE_CHECK(cusparseCreateDnVec(&ctx->vecNextEigen, m, ctx->next_eigenvector_d, CUDA_R_64F)); + CUSPARSE_CHECK(cusparseCreateDnVec(&ctx->vecDual, n, ctx->dual_product_d, CUDA_R_64F)); size_t bufferSizeAT = 0, bufferSizeA = 0; - cupdlpx_spmv_buffer_size(sparse_handle, matAT, vecNextEigen, vecDual, &bufferSizeAT); - cupdlpx_spmv_buffer_size(sparse_handle, matA, vecDual, vecEigen, &bufferSizeA); + cupdlpx_spmv_buffer_size(sparse_handle, ctx->matAT, ctx->vecNextEigen, ctx->vecDual, &bufferSizeAT); + cupdlpx_spmv_buffer_size(sparse_handle, ctx->matA, ctx->vecDual, ctx->vecEigen, &bufferSizeA); + CUDA_CHECK(cudaMalloc(&ctx->dBufferAT, bufferSizeAT)); + CUDA_CHECK(cudaMalloc(&ctx->dBufferA, bufferSizeA)); + cupdlpx_spmv_prepare( + sparse_handle, ctx->matAT, ctx->vecNextEigen, ctx->vecDual, ctx->dBufferAT, &ctx->descrAT, &ctx->planAT); + cupdlpx_spmv_prepare( + sparse_handle, ctx->matA, ctx->vecDual, ctx->vecEigen, ctx->dBufferA, &ctx->descrA, &ctx->planA); + return ctx; +} - CUDA_CHECK(cudaMalloc(&dBufferAT, bufferSizeAT)); - CUDA_CHECK(cudaMalloc(&dBufferA, bufferSizeA)); +void sv_estimator_free(sv_estimator_ctx_t *ctx) +{ + if (!ctx) + return; + cupdlpx_spmv_release(ctx->descrAT, ctx->planAT); + cupdlpx_spmv_release(ctx->descrA, ctx->planA); + CUDA_CHECK(cudaFree(ctx->dBufferAT)); + CUDA_CHECK(cudaFree(ctx->dBufferA)); + CUSPARSE_CHECK(cusparseDestroySpMat(ctx->matA)); + CUSPARSE_CHECK(cusparseDestroySpMat(ctx->matAT)); + CUSPARSE_CHECK(cusparseDestroyDnVec(ctx->vecEigen)); + CUSPARSE_CHECK(cusparseDestroyDnVec(ctx->vecNextEigen)); + CUSPARSE_CHECK(cusparseDestroyDnVec(ctx->vecDual)); + CUDA_CHECK(cudaFree(ctx->eigenvector_d)); + CUDA_CHECK(cudaFree(ctx->next_eigenvector_d)); + CUDA_CHECK(cudaFree(ctx->dual_product_d)); + free(ctx); +} + +sv_estimator_result_t sv_estimator_run(sv_estimator_ctx_t *ctx, const sv_estimator_opts_t *opts) +{ + cusparseHandle_t sparse_handle = ctx->sparse_handle; + cublasHandle_t blas_handle = ctx->blas_handle; + const int m = ctx->num_rows; + const int n = ctx->num_cols; + const bool *d_row_mask = opts->d_row_mask; + const bool *d_col_mask = opts->d_col_mask; + double *eigenvector_d = ctx->eigenvector_d; + double *next_eigenvector_d = ctx->next_eigenvector_d; + + sv_estimator_result_t result = {}; + result.status = SV_ESTIMATOR_MAX_ITER; + double max_singular_value_squared = 1.0; + const double one = 1.0; - cupdlpx_spmv_prepare(sparse_handle, matAT, vecNextEigen, vecDual, dBufferAT, &descrAT, &planAT); - cupdlpx_spmv_prepare(sparse_handle, matA, vecDual, vecEigen, dBufferA, &descrA, &planA); + cudaStream_t handle_stream = 0; + CUSPARSE_CHECK(cusparseGetStream(sparse_handle, &handle_stream)); + const int row_blocks = (m + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK; + const int col_blocks = (n + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK; + + bool random_start = !ctx->have_warm_start; + if (!random_start && d_row_mask) + { + elementwise_mask_kernel<<>>(eigenvector_d, d_row_mask, m); + /* Rows added to the mask after the previous estimate have zero entries in the + warm-start vector and could be omitted, causing the masked singular value to be + underestimated. Initialize these entries with small deterministic perturbations. + A masked warm vector that vanished (mask disjoint from the previous eigenvector's + support) carries no direction at all: start over from a random vector. */ + double warm_norm = 0.0; + CUBLAS_CHECK(cublasDnrm2_v2_64(blas_handle, m, eigenvector_d, 1, &warm_norm)); + if (isfinite(warm_norm) && warm_norm > 0.0) + { + double amplitude = 1e-3 * warm_norm / sqrt((double)m); + warm_start_fill_kernel<<>>( + eigenvector_d, d_row_mask, amplitude, m); + } + else + { + random_start = true; + } + } + if (random_start) + { + double *eigenvector_h = (double *)safe_malloc(m * sizeof(double)); + for (int i = 0; i < m; ++i) + { + eigenvector_h[i] = dist(gen); + } + CUDA_CHECK(cudaMemcpy(eigenvector_d, eigenvector_h, m * sizeof(double), cudaMemcpyHostToDevice)); + free(eigenvector_h); + if (d_row_mask) + { + elementwise_mask_kernel<<>>(eigenvector_d, d_row_mask, m); + } + } - for (int i = 0; i < max_iterations; ++i) + for (int i = 0; i < opts->max_iterations; ++i) { + result.iterations = i + 1; CUDA_CHECK(cudaMemcpy(next_eigenvector_d, eigenvector_d, m * sizeof(double), cudaMemcpyDeviceToDevice)); double eigenvector_norm; CUBLAS_CHECK(cublasDnrm2_v2_64(blas_handle, m, next_eigenvector_d, 1, &eigenvector_norm)); + if (!isfinite(eigenvector_norm) || eigenvector_norm <= 0.0) + { + result.status = SV_ESTIMATOR_DEGENERATE; + break; + } double inv_eigenvector_norm = 1.0 / eigenvector_norm; CUBLAS_CHECK(cublasDscal(blas_handle, m, &inv_eigenvector_norm, next_eigenvector_d, 1)); - cupdlpx_spmv_execute(sparse_handle, matAT, vecNextEigen, vecDual, dBufferAT, planAT); - cupdlpx_spmv_execute(sparse_handle, matA, vecDual, vecEigen, dBufferA, planA); + cupdlpx_spmv_execute(sparse_handle, ctx->matAT, ctx->vecNextEigen, ctx->vecDual, ctx->dBufferAT, ctx->planAT); + if (d_col_mask) + elementwise_mask_kernel<<>>( + ctx->dual_product_d, d_col_mask, n); + cupdlpx_spmv_execute(sparse_handle, ctx->matA, ctx->vecDual, ctx->vecEigen, ctx->dBufferA, ctx->planA); + if (d_row_mask) + elementwise_mask_kernel<<>>(eigenvector_d, d_row_mask, m); - CUBLAS_CHECK(cublasDdot(blas_handle, m, next_eigenvector_d, 1, eigenvector_d, 1, &sigma_max_sq)); + CUBLAS_CHECK(cublasDdot(blas_handle, m, next_eigenvector_d, 1, eigenvector_d, 1, &max_singular_value_squared)); + if (!isfinite(max_singular_value_squared) || max_singular_value_squared <= 0.0) + { + result.status = SV_ESTIMATOR_DEGENERATE; + break; + } + + if (opts->abort_singular_value_threshold > 0.0 && + sqrt(max_singular_value_squared) >= opts->abort_singular_value_threshold) + { + result.status = SV_ESTIMATOR_ABORTED; + break; + } - double neg_sigma_sq = -sigma_max_sq; - CUBLAS_CHECK(cublasDscal(blas_handle, m, &neg_sigma_sq, next_eigenvector_d, 1)); + double negative_max_singular_value_squared = -max_singular_value_squared; + CUBLAS_CHECK(cublasDscal(blas_handle, m, &negative_max_singular_value_squared, next_eigenvector_d, 1)); CUBLAS_CHECK(cublasDaxpy(blas_handle, m, &one, eigenvector_d, 1, next_eigenvector_d, 1)); double residual_norm; CUBLAS_CHECK(cublasDnrm2_v2_64(blas_handle, m, next_eigenvector_d, 1, &residual_norm)); + if (!isfinite(residual_norm)) + { + result.status = SV_ESTIMATOR_DEGENERATE; + break; + } - if (residual_norm < tolerance * fmin(1.0, sigma_max_sq)) + /* Use an absolute test above max_singular_value_squared = 1 and a relative test below it; + a purely absolute test is vacuously satisfied for operators with norm far below 1. */ + if (residual_norm < opts->tolerance * fmin(1.0, max_singular_value_squared)) + { + result.status = SV_ESTIMATOR_CONVERGED; break; + } } - CUDA_CHECK(cudaFree(dBufferAT)); - CUDA_CHECK(cudaFree(dBufferA)); - cupdlpx_spmv_release(descrAT, planAT); - cupdlpx_spmv_release(descrA, planA); - CUSPARSE_CHECK(cusparseDestroySpMat(matA)); - CUSPARSE_CHECK(cusparseDestroySpMat(matAT)); - CUSPARSE_CHECK(cusparseDestroyDnVec(vecEigen)); - CUSPARSE_CHECK(cusparseDestroyDnVec(vecNextEigen)); - CUSPARSE_CHECK(cusparseDestroyDnVec(vecDual)); - CUDA_CHECK(cudaFree(eigenvector_d)); - CUDA_CHECK(cudaFree(next_eigenvector_d)); - CUDA_CHECK(cudaFree(dual_product_d)); - - return sqrt(sigma_max_sq); + if (result.status != SV_ESTIMATOR_DEGENERATE) + { + result.max_singular_value = sqrt(max_singular_value_squared); + } + ctx->have_warm_start = result.status != SV_ESTIMATOR_DEGENERATE; + return result; } void compute_interaction_and_movement(pdhg_solver_state_t *state, double *interaction, double *movement) @@ -280,33 +405,36 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, const restart_parameters_t *restart_params, int termination_evaluation_frequency) { - bool do_restart = false; + const char *reason = NULL; if (solver_state->total_count == termination_evaluation_frequency) { - do_restart = true; + /* inner_count == total_count at the first check, so the long-inner-loop criterion holds */ + reason = "long inner loop"; } else if (solver_state->total_count > termination_evaluation_frequency) { if (solver_state->fixed_point_error <= restart_params->sufficient_reduction_for_restart * solver_state->initial_fixed_point_error) { - do_restart = true; + reason = "sufficient decay"; } - if (solver_state->fixed_point_error <= - restart_params->necessary_reduction_for_restart * solver_state->initial_fixed_point_error) + else if (solver_state->fixed_point_error <= + restart_params->necessary_reduction_for_restart * solver_state->initial_fixed_point_error && + solver_state->fixed_point_error > solver_state->last_trial_fixed_point_error) { - if (solver_state->fixed_point_error > solver_state->last_trial_fixed_point_error) - { - do_restart = true; - } + reason = "necessary decay + no local progress"; } - if (solver_state->inner_count >= restart_params->artificial_restart_threshold * solver_state->total_count) + else if (solver_state->inner_count >= restart_params->artificial_restart_threshold * solver_state->total_count) { - do_restart = true; + reason = "long inner loop"; } } solver_state->last_trial_fixed_point_error = solver_state->fixed_point_error; - return do_restart; + if (reason) + { + solver_state->last_restart_reason = reason; + } + return reason != NULL; } void set_default_parameters(pdhg_parameters_t *params) @@ -317,6 +445,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->pock_chambolle_alpha = 1.0; params->bound_objective_rescaling = true; params->verbose = true; + params->debug = false; params->termination_evaluation_frequency = 200; params->feasibility_polishing = false; params->reflection_coefficient = 1.0; @@ -342,6 +471,17 @@ void set_default_parameters(pdhg_parameters_t *params) params->presolve = true; params->matrix_zero_tol = 1e-9; params->infinite_bound = 1e20; + params->active_set_boost = true; + params->asb_activation_tol = 1e-4; + params->asb_window_iter = 10000; + params->asb_safety_factor = 0.9; + params->asb_max_reverts = 2; + params->asb_min_raise_ratio = 1.1; + params->asb_reestimate_change_ratio = 0.01; + params->asb_constraint_tol = 1e-8; + params->asb_variable_tol = 1e-8; + params->asb_divergence_ceiling_ratio = 0.7; + params->asb_divergence_margin = 0.05; } #define MATRIX_LARGE_VALUE 1e15 @@ -586,6 +726,26 @@ void restore_original_objective_sense(cupdlpx_result_t *result, objective_sense_ } \ } while (0) +/* Width of the iteration table for the given options; the banner uses it too. */ +static int iteration_table_width(const pdhg_parameters_t *params) +{ + const bool asb_columns = params->debug && params->active_set_boost; + return 88 + (params->debug ? 10 : 0) + (asb_columns ? 30 : 0); +} + +static void print_rule(int width) +{ + for (int i = 0; i < width; ++i) + putchar('-'); + putchar('\n'); +} + +static void print_centered(const char *text, int width) +{ + int pad = (width - (int)strlen(text)) / 2; + printf("%*s%s\n", pad > 0 ? pad : 0, "", text); +} + void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem) { pdhg_parameters_t default_params; @@ -594,17 +754,14 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro { return; } - printf("---------------------------------------------------------------------" - "------------------\n"); - printf(" cuPDLPx v%s " - " \n", - CUPDLPX_VERSION); - printf(" A GPU-Accelerated First-Order LP Solver " - " \n"); - printf(" (c) Haihao Lu, Massachusetts Institute of Technology, " - "2025 \n"); - printf("---------------------------------------------------------------------" - "------------------\n"); + const int width = iteration_table_width(params); + char version_line[64]; + snprintf(version_line, sizeof(version_line), "cuPDLPx v%s", CUPDLPX_VERSION); + print_rule(width); + print_centered(version_line, width); + print_centered("A GPU-Accelerated First-Order LP Solver", width); + print_centered("(c) Haihao Lu, Massachusetts Institute of Technology, 2025", width); + print_rule(width); printf("Problem: %d rows, %d columns, %d nonzeros\n", problem->num_constraints, @@ -640,8 +797,23 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro params->termination_criteria.eps_infeasible_relative, default_params.termination_criteria.eps_infeasible_relative); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); + PRINT_DIFF_BOOL("debug", params->debug, default_params.debug); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); PRINT_DIFF_DBL("infinite_bound", params->infinite_bound, default_params.infinite_bound); + PRINT_DIFF_BOOL("active_set_boost", params->active_set_boost, default_params.active_set_boost); + PRINT_DIFF_DBL("asb_activation_tol", params->asb_activation_tol, default_params.asb_activation_tol); + PRINT_DIFF_INT("asb_window_iter", params->asb_window_iter, default_params.asb_window_iter); + PRINT_DIFF_DBL("asb_safety_factor", params->asb_safety_factor, default_params.asb_safety_factor); + PRINT_DIFF_INT("asb_max_reverts", params->asb_max_reverts, default_params.asb_max_reverts); + PRINT_DIFF_DBL("asb_min_raise_ratio", params->asb_min_raise_ratio, default_params.asb_min_raise_ratio); + PRINT_DIFF_DBL( + "asb_reestimate_change_ratio", params->asb_reestimate_change_ratio, default_params.asb_reestimate_change_ratio); + PRINT_DIFF_DBL("asb_constraint_tol", params->asb_constraint_tol, default_params.asb_constraint_tol); + PRINT_DIFF_DBL("asb_variable_tol", params->asb_variable_tol, default_params.asb_variable_tol); + PRINT_DIFF_DBL("asb_divergence_ceiling_ratio", + params->asb_divergence_ceiling_ratio, + default_params.asb_divergence_ceiling_ratio); + PRINT_DIFF_DBL("asb_divergence_margin", params->asb_divergence_margin, default_params.asb_divergence_margin); } #undef PRINT_DIFF_INT @@ -652,8 +824,7 @@ void pdhg_final_log(const cupdlpx_result_t *result, const pdhg_parameters_t *par { if (params->verbose) { - printf("-------------------------------------------------------------------" - "--------------------\n"); + print_rule(iteration_table_width(params)); printf("Solution Summary\n"); printf(" Status : %s\n", termination_reason_to_string(result->termination_reason)); if (params->presolve) @@ -668,44 +839,76 @@ void pdhg_final_log(const cupdlpx_result_t *result, const pdhg_parameters_t *par printf(" Objective gap : %.3e\n", result->relative_objective_gap); printf(" Primal infeas : %.3e\n", result->relative_primal_residual); printf(" Dual infeas : %.3e\n", result->relative_dual_residual); + if (params->active_set_boost) + { + printf(" Active set boost : %d raises, %d reverts, %d PI iterations\n", + result->asb_raise_count, + result->asb_revert_count, + result->asb_pi_iterations); + } } +} - // if (stats != NULL && stats->n_rows_original > 0) { - // printf("\nPresolve Summary\n"); - // printf(" [Dimensions]\n"); - // printf(" Original : %d rows, %d cols, %d nnz\n", - // stats->n_rows_original, stats->n_cols_original, stats->nnz_original); - // printf(" Reduced : %d rows, %d cols, %d nnz\n", - // stats->n_rows_reduced, stats->n_cols_reduced, stats->nnz_reduced); - - // printf(" [Reduction Details (NNZ Removed)]\n"); - // printf(" Trivial : %d\n", stats->nnz_removed_trivial); - // printf(" Fast : %d\n", stats->nnz_removed_fast); - // printf(" Primal Propagation : %d\n", stats->nnz_removed_primal_propagation); - // printf(" Parallel Rows : %d\n", stats->nnz_removed_parallel_rows); - // printf(" Parallel Cols : %d\n", stats->nnz_removed_parallel_cols); - - // printf(" [Timing]\n"); - // printf(" Total Presolve : %.3g sec\n", stats->presolve_total_time); - // printf(" Init : %.3g sec\n", stats->ps_time_init); - // printf(" Fast : %.3g sec\n", stats->ps_time_fast); - // printf(" Medium : %.3g sec\n", stats->ps_time_medium); - // printf(" Primal Propagation : %.3g sec\n", stats->ps_time_primal_propagation); - // printf(" Parallel Rows : %.3g sec\n", stats->ps_time_parallel_rows); - // printf(" Parallel Cols : %.3g sec\n", stats->ps_time_parallel_cols); - // printf(" Postsolve : %.3g sec\n", stats->ps_time_post_solve); - // } +void display_iteration_header(const pdhg_parameters_t *params) +{ + if (!params->verbose) + { + return; + } + const bool asb_columns = params->debug && params->active_set_boost; + const int width = iteration_table_width(params); + printf("\n*: restart triggered\n"); + print_rule(width); + printf(" %s | %s | %s | %s", + " runtime ", + " objective ", + " absolute residuals ", + " relative residuals "); + if (params->debug) + { + printf(" | %s", " primal"); + } + if (asb_columns) + { + printf(" | %s", " active-set boost "); + } + printf(" \n"); + printf(" %s %s | %s %s | %s %s %s | %s %s %s", + " iter", + " time ", + " pr obj ", + " du obj ", + " pr res", + " du res", + " gap ", + " pr res", + " du res", + " gap "); + if (params->debug) + { + printf(" | %s", " weight"); + } + if (asb_columns) + { + printf(" | %s %s %s", " step ", "free var", "active con"); + } + printf(" \n"); + print_rule(width); } -void display_iteration_stats(const pdhg_solver_state_t *state, bool verbose) +void display_iteration_stats(pdhg_solver_state_t *state, const pdhg_parameters_t *params) { - if (!verbose) + if (!params->verbose) { return; } if (state->total_count % get_print_frequency(state->total_count) == 0) { - printf("%6d %.1e | %8.1e %8.1e | %.1e %.1e %.1e | %.1e %.1e %.1e \n", + /* a leading star marks that at least one restart happened since the previous row */ + char restart_marker = state->restart_count > state->logged_restart_count ? '*' : ' '; + state->logged_restart_count = state->restart_count; + printf("%c%6d %.1e | %8.1e %8.1e | %.1e %.1e %.1e | %.1e %.1e %.1e", + restart_marker, state->total_count, state->cumulative_time_sec, state->original_objective_sign * state->primal_objective_value, @@ -716,6 +919,22 @@ void display_iteration_stats(const pdhg_solver_state_t *state, bool verbose) state->relative_primal_residual, state->relative_dual_residual, state->relative_objective_gap); + if (params->debug) + { + printf(" | %.1e", state->primal_weight); + } + if (params->debug && params->active_set_boost) + { + long free_variables = state->asb_free_variables; + long binding_constraints = state->asb_binding_constraints; + if (state->total_count == 0) + { + free_variables = state->num_variables; + binding_constraints = state->num_constraints; + } + printf(" | %.1e %8ld %10ld", state->step_size, free_variables, binding_constraints); + } + printf(" \n"); } } @@ -1021,15 +1240,15 @@ void compute_infeasibility_information(pdhg_solver_state_t *state) THREADS_PER_BLOCK, 0, state->stream>>>(state->dual_product, - state->dual_slack, + state->infeasibility_dual_scratch, state->variable_lower_bound_finite_val, state->variable_upper_bound_finite_val, state->num_variables); double sum_primal_slack = get_vector_sum(state->blas_handle, state->num_constraints, state->ones_dual_d, state->primal_slack); - double sum_dual_slack = - get_vector_sum(state->blas_handle, state->num_variables, state->ones_primal_d, state->dual_slack); + double sum_dual_slack = get_vector_sum( + state->blas_handle, state->num_variables, state->ones_primal_d, state->infeasibility_dual_scratch); state->dual_ray_objective = (sum_primal_slack + sum_dual_slack) / (state->constraint_bound_rescaling * state->objective_vector_rescaling); @@ -1045,15 +1264,16 @@ void compute_infeasibility_information(pdhg_solver_state_t *state) state->variable_lower_bound, state->variable_upper_bound, state->num_variables, - state->dual_slack, + state->infeasibility_dual_scratch, state->variable_rescaling); state->max_primal_ray_infeasibility = get_vector_inf_norm(state->blas_handle, state->num_constraints, state->primal_slack); - double dual_slack_norm = get_vector_inf_norm(state->blas_handle, state->num_variables, state->dual_slack); - state->max_dual_ray_infeasibility = dual_slack_norm; + double dual_scratch_norm = + get_vector_inf_norm(state->blas_handle, state->num_variables, state->infeasibility_dual_scratch); + state->max_dual_ray_infeasibility = dual_scratch_norm; - double scaling_factor = fmax(dual_ray_inf_norm, dual_slack_norm); + double scaling_factor = fmax(dual_ray_inf_norm, dual_scratch_norm); if (scaling_factor > 0.0) { state->max_dual_ray_infeasibility /= scaling_factor; diff --git a/test/test_active_set_boost.py b/test/test_active_set_boost.py new file mode 100644 index 0000000..ffc9229 --- /dev/null +++ b/test/test_active_set_boost.py @@ -0,0 +1,224 @@ +# Copyright 2025 Haihao Lu +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +import pytest +import scipy.sparse as sp + +from cupdlpx import Model, PDLP +from cupdlpx._core import get_default_params, solve_once, validate_params + +# alias -> (backend key, default) +ASB_PARAMS = { + "ActiveSetBoost": ("active_set_boost", True), + "ASBActivationTol": ("asb_activation_tol", 1e-4), + "ASBWindowIter": ("asb_window_iter", 10000), + "ASBSafetyFactor": ("asb_safety_factor", 0.9), + "ASBMaxReverts": ("asb_max_reverts", 2), + "ASBMinRaiseRatio": ("asb_min_raise_ratio", 1.1), + "ASBReestimateChangeRatio": ("asb_reestimate_change_ratio", 0.01), + "ASBConstraintTol": ("asb_constraint_tol", 1e-8), + "ASBVariableTol": ("asb_variable_tol", 1e-8), + "ASBDivergenceCeilingRatio": ("asb_divergence_ceiling_ratio", 0.7), + "ASBDivergenceMargin": ("asb_divergence_margin", 0.05), +} +ASB_INFO_KEYS = ("ASBRaiseCount", "ASBRevertCount", "ASBPowerIterations") + + +def _model(base_lp_data, **params): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False, **params) + return model + + +def test_defaults_exposed_through_every_layer(base_lp_data): + defaults = get_default_params() + assert defaults["debug"] is False + model = _model(base_lp_data) + for alias, (key, default) in ASB_PARAMS.items(): + assert PDLP._PARAM_ALIAS[alias] == key + assert defaults[key] == default + assert model.getParam(alias) == default + assert model.getParam(key) == default + assert PDLP._PARAM_ALIAS["Debug"] == "debug" + assert model.getParam("Debug") is False + + +@pytest.mark.parametrize( + "alias, value", + [ + ("ASBWindowIter", 0), + ("ASBMaxReverts", -1), + ("ASBActivationTol", -1e-9), + ("ASBSafetyFactor", 0.0), + ("ASBMinRaiseRatio", 0.5), + ("ASBReestimateChangeRatio", -0.1), + ("ASBConstraintTol", -1e-9), + ("ASBVariableTol", -1e-9), + ("ASBDivergenceCeilingRatio", 0.0), + ("ASBDivergenceCeilingRatio", 1.5), + ("ASBDivergenceMargin", -0.1), + ("TermCheckFreq", 2), + ], +) +def test_out_of_range_values_fail_at_set_param(base_lp_data, alias, value): + """setParam runs the C validator, so a bad value is rejected before optimize() + with the backend's message and the stored parameters stay untouched.""" + model = _model(base_lp_data) + before = dict(model.Params.items()) + key = PDLP._PARAM_ALIAS[alias] + with pytest.raises(ValueError, match=key): + model.setParam(alias, value) + assert dict(model.Params.items()) == before + with pytest.raises(ValueError, match=key): + validate_params({key: value}) + + +def test_boundary_values_accepted(base_lp_data): + model = _model(base_lp_data) + model.setParam("ASBMinRaiseRatio", 1.0) + model.setParam("ASBDivergenceCeilingRatio", 1.0) + model.setParam("ASBMaxReverts", 0) + model.setParam("ASBActivationTol", 0.0) + model.setParam("TermCheckFreq", 3) + model.setParam("ActiveSetBoost", 0) + assert model.getParam("ActiveSetBoost") is False + model.optimize() + assert model.Status == PDLP.OPTIMAL + + +@pytest.mark.parametrize( + "params", + [ + {"asb_min_raise_ratio": 0.5}, + {"asb_divergence_ceiling_ratio": 1.5}, + {"asb_window_iter": 0}, + {"asb_safety_factor": float("nan")}, + ], +) +def test_core_rejects_out_of_range_values(base_lp_data, params): + c, A, l, u, lb, ub = base_lp_data + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params=params) + + +def _solve_info(base_lp_data, **params): + c, A, l, u, lb, ub = base_lp_data + base = {"verbose": False, "presolve": False} + base.update(params) + return solve_once(A, c, None, lb, ub, l, u, params=base) + + +def test_info_reports_controller_statistics(base_lp_data): + info = _solve_info(base_lp_data) + assert info["Status"] == "OPTIMAL" + for key in ASB_INFO_KEYS: + assert isinstance(info[key], int) + assert info[key] >= 0 + + +def test_boost_off_reports_zero_statistics(base_lp_data): + info = _solve_info(base_lp_data, active_set_boost=False) + assert info["Status"] == "OPTIMAL" + assert all(info[key] == 0 for key in ASB_INFO_KEYS) + + +def test_boost_on_and_off_reach_the_same_solution(base_lp_data, atol): + on = _solve_info(base_lp_data, active_set_boost=True) + off = _solve_info(base_lp_data, active_set_boost=False) + assert on["Status"] == off["Status"] == "OPTIMAL" + assert np.allclose(on["X"], off["X"], atol=atol) + assert abs(on["PrimalObj"] - off["PrimalObj"]) <= atol + + +def test_activated_controller_solves_to_tight_tolerance(base_lp_data, atol): + """Activate immediately and solve tightly so the controller runs its full path.""" + info = _solve_info( + base_lp_data, + eps_optimal_relative=1e-8, + eps_feasible_relative=1e-8, + asb_activation_tol=1.0, + asb_window_iter=1, + ) + assert info["Status"] == "OPTIMAL" + assert info["ASBRevertCount"] <= get_default_params()["asb_max_reverts"] + assert np.allclose(info["X"], [1.0, 2.0], atol=atol) + + +def test_debug_implies_verbose_and_solves(base_lp_data): + model = _model(base_lp_data, Debug=True) + assert model.getParam("Debug") is True + assert model.getParam("OutputFlag") is False + model.optimize() + assert model.Status == PDLP.OPTIMAL + + +def _hadamard(k): + H = np.array([[1.0]]) + while H.shape[0] < k: + H = np.block([[H, H], [H, -H]]) + return H + + +def _column_mask_lp(k=16): + """ + Rows H x + J y = b with H Hadamard (k x k) and J all-ones (k x k). + Minimize -sum(y), 0 <= x, y <= 1. + At the optimum x = 0.5 is interior (free) and y = 1 sits on its upper + bound (clamped). All rows are equalities, so every row is binding; the + dense block J lives only in the clamped columns. After rescaling the full + matrix has sigma ~ 0.71 while the masked matrix H alone has 1/sqrt(2k), + so a raise happens only if the column mask correctly drops y. + """ + H = _hadamard(k) + J = np.ones((k, k)) + x_star = np.full(k, 0.5) + y_star = np.ones(k) + b = H @ x_star + J @ y_star + A = sp.csr_matrix(np.hstack([H, J])) + c = np.concatenate([np.zeros(k), -np.ones(k)]) + lb = np.zeros(2 * k) + ub = np.ones(2 * k) + return c, A, b.copy(), b.copy(), lb, ub, np.concatenate([x_star, y_star]) + + +def test_raise_follows_the_column_mask(atol): + """The boost must raise the step on an instance whose gain comes from clamped columns. + + Regression guard: the infeasibility-certificate routine once reused + ``dual_slack`` as scratch, which wiped the clamping signal the column + mask reads; the masked matrix then equalled the full one and no raise + was possible. + """ + c, A, l, u, lb, ub, x_star = _column_mask_lp() + params = { + "verbose": False, + "presolve": False, + "eps_optimal_relative": 1e-8, + "eps_feasible_relative": 1e-8, + "termination_evaluation_frequency": 10, + "asb_activation_tol": 1.0, + "asb_window_iter": 1, + } + info = solve_once(A, c, None, lb, ub, l, u, params=params) + assert info["Status"] == "OPTIMAL" + assert np.allclose(info["X"], x_star, atol=atol) + assert info["ASBPowerIterations"] >= 1 + assert info["ASBRaiseCount"] >= 1 + + off = solve_once(A, c, None, lb, ub, l, u, params={**params, "active_set_boost": False}) + assert off["Status"] == "OPTIMAL" + assert np.allclose(off["X"], x_star, atol=atol) + assert off["ASBRaiseCount"] == 0 diff --git a/test/test_api_surface.py b/test/test_api_surface.py index b418ba5..29ec2ce 100644 --- a/test/test_api_surface.py +++ b/test/test_api_surface.py @@ -259,6 +259,12 @@ def test_direct_core_param_value_validation(base_lp_data): solve_once(A, c, None, lb, ub, l, u, params={"time_sec_limit": float("nan")}) with pytest.raises(ValueError): solve_once(A, c, None, lb, ub, l, u, params={"optimality_norm": "l1"}) + with pytest.raises(ValueError, match="Unknown parameter"): + solve_once(A, c, None, lb, ub, l, u, params={"eps_optimal_relatve": 1e-6}) + with pytest.raises(ValueError, match="must be a number"): + solve_once(A, c, None, lb, ub, l, u, params={"time_sec_limit": "60"}) + with pytest.raises(ValueError, match="int32"): + solve_once(A, c, None, lb, ub, l, u, params={"iteration_limit": 2**40}) def test_direct_core_model_data_validation(base_lp_data): diff --git a/test/test_interface.c b/test/test_interface.c index f62e773..a3935c0 100644 --- a/test/test_interface.c +++ b/test/test_interface.c @@ -227,7 +227,7 @@ int main() // Test 9: GPU solver path (presolve disabled) -- forces hipBLAS/hipSPARSE execution printf("\n=== Test 9: CSR Matrix (presolve disabled, GPU solver) ===\n"); { - lp_problem_t *prob9 = create_lp_problem(c, &A_csr, l, u, NULL, NULL, NULL); + lp_problem_t *prob9 = create_lp_problem(c, &A_csr, l, u, NULL, NULL, NULL, NULL); if (!prob9) { fprintf(stderr, "[test] create_lp_problem failed for Test 9.\n");