L-infinity Norm Termination Criteria

.gitignore

@@ -41,6 +41,7 @@
 *.su
 *.idb
 *.pdb
+Testing/*
 
 # Kernel Module Compile Results
 *.mod*
@@ -60,3 +61,5 @@ test/*
 /.vscode
 /.venv
 /_b
+*.whl
+*.pyc

README.md

@@ -84,6 +84,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the
 | `-v`, `--verbose` | `flag` | Enable verbose logging. | `false` |
 | `--time_limit` | `double` | Time limit in seconds. | `3600.0` |
 | `--iter_limit` | `int` | Iteration limit. | `2147483647` |
+| `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` |
 | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` |
 | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` |
 | `--eps_infeas_detect` | `double` | Infeasibility detection tolerance. | `1e-10` |

include/cupdlpx_types.h

@@ -37,6 +37,12 @@ extern "C"
 		TERMINATION_REASON_FEAS_POLISH_SUCCESS
 	} termination_reason_t;
 
+	typedef enum
+    {
+        NORM_TYPE_L2 = 0,
+        NORM_TYPE_L_INF = 1
+    } norm_type_t;
+
 	typedef struct
 	{
 		int num_variables;
@@ -93,6 +99,7 @@ extern "C"
 		restart_parameters_t restart_params;
 		double reflection_coefficient;
 		bool feasibility_polishing;
+		norm_type_t optimality_norm;
 		bool presolve;
 	} pdhg_parameters_t;
 

internal/internal_types.h

@@ -103,6 +103,7 @@ typedef struct
 	double initial_fixed_point_error;
 	double last_trial_fixed_point_error;
 	int inner_count;
+	norm_type_t optimality_norm;
 
 	cusparseHandle_t sparse_handle;
 	cublasHandle_t blas_handle;

python/README.md

@@ -141,6 +141,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio
 | `IterationLimit` | `iteration_limit` | int | `2147483647` | Maximum number of iterations. |
 | `OutputFlag`, `LogToConsole` | `verbose` | bool | `False` | Enable (`True`) or disable (`False`) console logging output. |
 | `TermCheckFreq` | `termination_evaluation_frequency` | int | `200` | Frequency (in iterations) at which termination conditions are evaluated. |
+| `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. |
 | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. |
 | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. |
 | `InfeasibleTol` | `eps_infeasible` | float | `1e-10` | Threshold for declaring infeasibility. |

python/cupdlpx/PDLP.py

@@ -53,6 +53,8 @@
     # feasibility polishing
     "FeasibilityPolishing": "feasibility_polishing",
     "FeasibilityPolishingTol": "eps_feas_polish_relative",
+    # termination criteria
+    "OptimalityNorm": "optimality_norm",
     # singular value estimation (power method)
     "SVMaxIter": "sv_max_iter",
     "SVTol": "sv_tol",

python_bindings/_core_bindings.cpp

@@ -137,6 +137,21 @@ static const int32_t *get_index_ptr_i32(py::object obj, const char *name,
     throw std::invalid_argument(std::string(name) + " must be int32 or int64.");
 }
 
+// helper function to convert norm string to enum
+static norm_type_t parse_norm_string(const std::string& s)
+{
+    std::string lower = s;
+    std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower);
+
+    if (lower == "l2") {
+        return NORM_TYPE_L2;
+    } else if (lower == "linf") {
+        return NORM_TYPE_L_INF;
+    } else {
+        throw std::invalid_argument("Unknown norm type: " + s + ". Use 'l2' or 'linf'.");
+    }
+}
+
 // ensure 1D array or None with expected length
 static void ensure_len_or_null(py::object obj, const char *name, int expect_len)
 {
@@ -247,6 +262,8 @@ static py::dict get_default_params_py()
     d["feasibility_polishing"] = p.feasibility_polishing;
     d["eps_feas_polish_relative"] = p.termination_criteria.eps_feas_polish_relative;
 
+    // Termination criteria norm
+    d["optimality_norm"] = (p.optimality_norm == NORM_TYPE_L_INF) ? "linf" : "l2";
     // power method for singular value estimation
     d["sv_max_iter"] = p.sv_max_iter;
     d["sv_tol"] = p.sv_tol;
@@ -270,6 +287,19 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
     { if (d.contains(k)) tgt = py::cast<int>(d[k]); };
     auto getb = [&](const char *k, bool &tgt)
     { if (d.contains(k)) tgt = py::cast<bool>(d[k]); };
+    auto get_norm = [&](const char *k, norm_type_t &tgt)
+    { 
+        if (d.contains(k)) {
+            py::object val = d[k];
+            if (py::isinstance<py::str>(val)) {
+                std::string sval = py::cast<std::string>(val);
+                tgt = parse_norm_string(sval);
+            }
+            else {
+                throw std::invalid_argument("optimality_norm must be a string ('l2'/'linf')");
+            }
+        }
+    };
 
     // verbosity
     getb("verbose", p->verbose);
@@ -303,6 +333,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
     getb("feasibility_polishing", p->feasibility_polishing);
     getf("eps_feas_polish_relative", p->termination_criteria.eps_feas_polish_relative);
 
+    // Termination criteria norm
+    get_norm("optimality_norm", p->optimality_norm);
     // power method for singular value estimation
     geti("sv_max_iter", p->sv_max_iter);
     getf("sv_tol", p->sv_tol);

src/cli.c

@@ -201,6 +201,8 @@ void print_usage(const char *prog_name)
                     "Enable feasibility use feasibility polishing (default: false).\n");
     fprintf(stderr, "      --eps_feas_polish <tolerance>   Relative feasibility "
                     "polish tolerance (default: 1e-6).\n");
+    fprintf(stderr, "      --opt_norm <norm_type>          "
+                    "Norm for optimality criteria: l2 or linf (default: l2).\n");
     fprintf(stderr, "      --no_presolve                   "
                     "Disable presolve (default: enabled).\n");
 }
@@ -227,7 +229,8 @@ int main(int argc, char *argv[])
         {"sv_max_iter", required_argument, 0, 1011},
         {"sv_tol", required_argument, 0, 1012},
         {"eval_freq", required_argument, 0, 1013},
-        {"no_presolve", no_argument, 0, 1014},
+        {"opt_norm", required_argument, 0, 1014},
+        {"no_presolve", no_argument, 0, 1015},
         {0, 0, 0, 0}};
 
     int opt;
@@ -283,7 +286,20 @@ int main(int argc, char *argv[])
         case 1013: // --eval_freq
             params.termination_evaluation_frequency = atoi(optarg);
             break;
-        case 1014: // --no_presolve
+        case 1014: // --opt_norm
+            {
+                const char *norm_str = optarg;
+                if (strcmp(norm_str, "l2") == 0) {
+                    params.optimality_norm = NORM_TYPE_L2;
+                } else if (strcmp(norm_str, "linf") == 0) {
+                    params.optimality_norm = NORM_TYPE_L_INF;
+                } else {
+                    fprintf(stderr, "Error: opt_norm must be 'l2' or 'linf'\n");
+                    return 1;
+                }
+            }
+            break;
+        case 1015: // --no_presolve
             params.presolve = false;
             break;
         case '?': // Unknown option

src/solver.cu

@@ -116,6 +116,8 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params,
     rescale_info_t *rescale_info = rescale_problem(params, working_problem);
     pdhg_solver_state_t *state = initialize_solver_state(params, working_problem, rescale_info);
 
+    state->optimality_norm = params ? params->optimality_norm : NORM_TYPE_L2;
+
     rescale_info_free(rescale_info);
     initialize_step_size_and_primal_weight(state, params);
     clock_t start_time = clock();
@@ -408,32 +410,59 @@ initialize_solver_state(const pdhg_parameters_t *params,
     free(temp_host);
 
     double sum_of_squares = 0.0;
+    double max_val = 0.0;
+    double val = 0.0;
 
     for (int i = 0; i < n_vars; ++i)
     {
-        sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i];
+        if (state->optimality_norm == NORM_TYPE_L_INF) {
+            val = fabs(working_problem->objective_vector[i]);
+            if (val > max_val) max_val = val;
+        } else {
+            sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i];
+        }
+    }
+
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->objective_vector_norm = max_val;
+    } else {
+        state->objective_vector_norm = sqrt(sum_of_squares);
     }
-    state->objective_vector_norm = sqrt(sum_of_squares);
 
     sum_of_squares = 0.0;
+    max_val = 0.0;
+    val = 0.0;
 
     for (int i = 0; i < n_cons; ++i)
     {
         double lower = working_problem->constraint_lower_bound[i];
         double upper = working_problem->constraint_upper_bound[i];
 
-        if (isfinite(lower) && (lower != upper))
-        {
-            sum_of_squares += lower * lower;
+        if (state->optimality_norm == NORM_TYPE_L_INF) {
+            if (isfinite(lower) && (lower != upper)) {
+                val = fabs(lower);
+                if (val > max_val) max_val = val;
+            }
+            if (isfinite(upper)) {
+                val = fabs(upper);
+                if (val > max_val) max_val = val;
+            }
+        } else {
+            if (isfinite(lower) && (lower != upper)) {
+                sum_of_squares += lower * lower;
+            }
+            if (isfinite(upper)) {
+                sum_of_squares += upper * upper;
+            }
         }
+    }
 
-        if (isfinite(upper))
-        {
-            sum_of_squares += upper * upper;
-        }
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->constraint_bound_norm = max_val;
+    } else {
+        state->constraint_bound_norm = sqrt(sum_of_squares);
     }
 
-    state->constraint_bound_norm = sqrt(sum_of_squares);
     state->num_blocks_primal =
         (state->num_variables + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
     state->num_blocks_dual =

src/utils.cu

@@ -224,8 +224,8 @@ bool optimality_criteria_met(const pdhg_solver_state_t *state,
                              double rel_opt_tol, double rel_feas_tol)
 {
     return state->relative_dual_residual < rel_feas_tol &&
-           state->relative_primal_residual < rel_feas_tol &&
-           state->relative_objective_gap < rel_opt_tol;
+            state->relative_primal_residual < rel_feas_tol &&
+            state->relative_objective_gap < rel_opt_tol;
 }
 
 bool primal_infeasibility_criteria_met(const pdhg_solver_state_t *state,
@@ -349,6 +349,7 @@ void set_default_parameters(pdhg_parameters_t *params)
     params->restart_params.k_d = 0.0;
     params->restart_params.i_smooth = 0.3;
 
+    params->optimality_norm = NORM_TYPE_L2;
     params->presolve = true;
 }
 
@@ -707,13 +708,26 @@ void compute_residual(pdhg_solver_state_t *state)
         state->constraint_upper_bound_finite_val, state->num_constraints,
         state->num_variables);
 
-    CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints,
-                                   state->primal_residual, 1,
-                                   &state->absolute_primal_residual));
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, 
+                                          state->num_constraints, state->primal_residual);
+    } else {
+        CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, 
+                                       state->primal_residual, 1, 
+                                       &state->absolute_primal_residual));
+    }
+
     state->absolute_primal_residual /= state->constraint_bound_rescaling;
-    CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables,
-                                   state->dual_residual, 1,
-                                   &state->absolute_dual_residual));
+
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, 
+                                        state->num_variables, state->dual_residual);
+    } else {
+        CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables,
+                                       state->dual_residual, 1,
+                                       &state->absolute_dual_residual));
+    }
+
     state->absolute_dual_residual /= state->objective_vector_rescaling;
 
     CUBLAS_CHECK(cublasDdot(
@@ -736,12 +750,15 @@ void compute_residual(pdhg_solver_state_t *state)
                                        state->objective_vector_rescaling) +
                                   state->objective_constant;
 
-    state->relative_primal_residual =
+    state->relative_primal_residual = 
         state->absolute_primal_residual / (1.0 + state->constraint_bound_norm);
+
     state->relative_dual_residual =
         state->absolute_dual_residual / (1.0 + state->objective_vector_norm);
+
     state->objective_gap =
         fabs(state->primal_objective_value - state->dual_objective_value);
+
     state->relative_objective_gap =
         state->objective_gap / (1.0 + fabs(state->primal_objective_value) +
                                 fabs(state->dual_objective_value));
@@ -1230,8 +1247,17 @@ void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_
         state->constraint_upper_bound, state->constraint_rescaling,
         state->num_constraints);
 
-    CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->primal_residual, 1, &state->absolute_primal_residual));
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, 
+                                          state->num_constraints, state->primal_residual);
+    } else {
+        CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, 
+                                       state->primal_residual, 1, 
+                                       &state->absolute_primal_residual));
+    }
+
     state->absolute_primal_residual /= state->constraint_bound_rescaling;
+
     state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm);
 
     CUBLAS_CHECK(cublasDdot(state->blas_handle, state->num_variables, ori_state->objective_vector, 1, state->pdhg_primal_solution, 1, &state->primal_objective_value));
@@ -1256,8 +1282,16 @@ void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_so
         ori_state->constraint_upper_bound_finite_val,
         state->num_variables, state->num_constraints);
 
-    CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->dual_residual, 1, &state->absolute_dual_residual));
+    if (state->optimality_norm == NORM_TYPE_L_INF) {
+        state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, 
+                                        state->num_variables, state->dual_residual);
+    } else {
+        CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, 
+                                       state->dual_residual, 1, &state->absolute_dual_residual));
+    }
+
     state->absolute_dual_residual /= state->objective_vector_rescaling;
+
     state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm);
 
     double base_dual_objective;