feat: active-set stepsize boosting

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 <output_directory>. The filenames are derived from the input file's basename. For an input `INSTANCE.mps.gz`, the output will be:

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);

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

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

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

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;

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,14 +141,18 @@ 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);
 
     int get_print_frequency(int iter);
 
     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);

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:
 

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",
 }

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
+from ._cupdlpx_core import solve_once, get_default_params, validate_params, read_mps