Feat: infeasibility detection

include/cupdlpx_types.h

@@ -87,6 +87,7 @@ extern "C"
         double eps_optimal_relative;
         double eps_feasible_relative;
         double eps_feas_polish_relative;
+        double eps_infeasible_relative;
         double time_sec_limit;
         int iteration_limit;
     } termination_criteria_t;

python/README.md

@@ -163,6 +163,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio
 | `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_relative` | float | `1e-14` | Relative tolerance on the ray certificate for declaring primal/dual infeasibility. |
 | `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. |
 | `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations of L∞ Ruiz scaling. Improves numerical conditioning. |
 | `UsePCAlpha` | `has_pock_chambolle_alpha` | bool | `True` | Whether to use the Pock–Chambolle α step size adjustment. |

python/cupdlpx/PDLP.py

@@ -41,6 +41,7 @@
     # tolerances
     "OptimalityTol": "eps_optimal_relative",
     "FeasibilityTol": "eps_feasible_relative",
+    "InfeasibleTol": "eps_infeasible_relative",
     # scaling / step size
     "GeoMeanIters": "geometric_mean_iterations",
     "RuizIters": "l_inf_ruiz_iterations",

python/cupdlpx/model.py

@@ -53,6 +53,7 @@
     {
         "eps_optimal_relative",
         "eps_feasible_relative",
+        "eps_infeasible_relative",
         "time_sec_limit",
         "pock_chambolle_alpha",
         "artificial_restart_threshold",
@@ -70,6 +71,7 @@
     {
         "eps_optimal_relative",
         "eps_feasible_relative",
+        "eps_infeasible_relative",
         "eps_feas_polish_relative",
         "sv_tol",
         "infinite_bound",

python_bindings/_core_bindings.cpp

@@ -273,6 +273,7 @@ static py::dict get_default_params_py()
     // tolerances
     d["eps_optimal_relative"] = p.termination_criteria.eps_optimal_relative;
     d["eps_feasible_relative"] = p.termination_criteria.eps_feasible_relative;
+    d["eps_infeasible_relative"] = p.termination_criteria.eps_infeasible_relative;
 
     // limits
     d["time_sec_limit"] = p.termination_criteria.time_sec_limit;
@@ -384,6 +385,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
     // tolerances
     getf("eps_optimal_relative", p->termination_criteria.eps_optimal_relative);
     getf("eps_feasible_relative", p->termination_criteria.eps_feasible_relative);
+    getf("eps_infeasible_relative", p->termination_criteria.eps_infeasible_relative);
 
     // limits
     getf("time_sec_limit", p->termination_criteria.time_sec_limit);
@@ -435,6 +437,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
         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)

src/solver.cu

@@ -203,6 +203,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t *
         compute_fixed_point_error(state);
 
         compute_residual(state, params->optimality_norm);
+        compute_infeasibility_information(state);
         state->inner_count += params->termination_evaluation_frequency;
         state->total_count += params->termination_evaluation_frequency;
 

src/utils.cu

@@ -254,6 +254,16 @@ void check_termination_criteria(pdhg_solver_state_t *solver_state, const termina
         solver_state->termination_reason = TERMINATION_REASON_OPTIMAL;
         return;
     }
+    if (primal_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative))
+    {
+        solver_state->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE;
+        return;
+    }
+    if (dual_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative))
+    {
+        solver_state->termination_reason = TERMINATION_REASON_DUAL_INFEASIBLE;
+        return;
+    }
     if (solver_state->total_count >= criteria->iteration_limit)
     {
         solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT;
@@ -318,6 +328,7 @@ void set_default_parameters(pdhg_parameters_t *params)
     params->termination_criteria.time_sec_limit = 3600.0;
     params->termination_criteria.iteration_limit = INT32_MAX;
     params->termination_criteria.eps_feas_polish_relative = 1e-6;
+    params->termination_criteria.eps_infeasible_relative = 1e-14;
 
     params->restart_params.artificial_restart_threshold = 0.36;
     params->restart_params.sufficient_reduction_for_restart = 0.2;
@@ -625,6 +636,9 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro
     PRINT_DIFF_DBL("eps_feas_polish_relative",
                    params->termination_criteria.eps_feas_polish_relative,
                    default_params.termination_criteria.eps_feas_polish_relative);
+    PRINT_DIFF_DBL("eps_infeasible_relative",
+                   params->termination_criteria.eps_infeasible_relative,
+                   default_params.termination_criteria.eps_infeasible_relative);
     PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve);
     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);

test/test_infeasible_unbounded.py

@@ -35,13 +35,15 @@ def test_infeasible_lp(base_lp_data, atol):
     model.setConstraintUpperBound(u)
     # turn off output
     model.setParams(OutputFlag=False, Presolve=False)
+    # set infeasibility tolerance
+    model.setParams(InfeasibleTol=1e-6)
     # optimize
-    #model.optimize()
+    model.optimize()
     # check status
-    #assert hasattr(model, "Status"), "Model.Status not exposed."
-    #assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}"
+    assert hasattr(model, "Status"), "Model.Status not exposed."
+    assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}"
     # check dual ray
-    #assert model.DualRayObj > atol, f"DualRayObj should be positive for dual infeasible, got {model.DualRayObj}"
+    assert model.DualRayObj > atol, f"DualRayObj should be positive for primal infeasible, got {model.DualRayObj}"
 
 
 def test_unbounded_lp(base_lp_data, atol):
@@ -64,10 +66,12 @@ def test_unbounded_lp(base_lp_data, atol):
     model.setVariableLowerBound(lb)
     # turn off output
     model.setParams(OutputFlag=False, Presolve=False)
+    # set infeasibility tolerance
+    model.setParams(InfeasibleTol=1e-6)
     # optimize
-    #model.optimize()
+    model.optimize()
     # check status
-    #assert hasattr(model, "Status"), "Model.Status not exposed."
-    #assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}"
+    assert hasattr(model, "Status"), "Model.Status not exposed."
+    assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}"
     # check primal ray
-    #assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}"
+    assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}"