Skip to content

Commit 42578ad

Browse files
jcarpentjcarpent
authored
Merge pull request #241 from Bambade/add_infeasibility_solving
Add infeasibility solving feature for dense and sparse backends
2 parents 97cb913 + 1fe0e3e commit 42578ad

26 files changed

+718
-257
lines changed

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@ We are ready to integrate **ProxSuite** within other optimization ecosystems.
4242
- Python and Julia bindings for easy code prototyping without sacrificing performance.
4343

4444
**Proxsuite** has a dedicated feature for solving batch of QPs.
45+
**Proxsuite** has a dedicated feature for solving closest feasible QPs (in $$\ell_2$$ sense) if they appear to be primal infeasible.
4546
**Proxsuite** is extensible.
4647
**Proxsuite** is reliable and extensively tested, showing the best performances on the hardest problems of the literature.
4748
**Proxsuite** is supported and tested on Windows, Mac OS X, Unix and Linux.

‎bindings/python/src/expose-results.hpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,8 @@ exposeResults(pybind11::module_ m)
2525
.value("PROXQP_SOLVED", QPSolverOutput::PROXQP_SOLVED)
2626
.value("PROXQP_MAX_ITER_REACHED", QPSolverOutput::PROXQP_MAX_ITER_REACHED)
2727
.value("PROXQP_PRIMAL_INFEASIBLE", QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE)
28+
.value("PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE",
29+
QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE)
2830
.value("PROXQP_DUAL_INFEASIBLE", QPSolverOutput::PROXQP_DUAL_INFEASIBLE)
2931
.value("PROXQP_NOT_RUN", QPSolverOutput::PROXQP_NOT_RUN)
3032
.export_values();
@@ -66,6 +68,14 @@ exposeResults(pybind11::module_ m)
6668
Results<T>,
6769
z,
6870
"The dual solution associated to the inequality constraints.")
71+
.PROXSUITE_PYTHON_EIGEN_READWRITE(
72+
Results<T>,
73+
se,
74+
"Optimal shift to the closest feasible problem wrt equality constraints.")
75+
.PROXSUITE_PYTHON_EIGEN_READWRITE(Results<T>,
76+
si,
77+
"Pptimal shift to the closest feasible "
78+
"problem wrt inequality constraints.")
6979
.def_readwrite("info", &Results<T>::info)
7080
.def(pybind11::self == pybind11::self)
7181
.def(pybind11::self != pybind11::self)

‎bindings/python/src/expose-settings.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,10 @@ exposeSettings(pybind11::module_ m)
8383
.def_readwrite("bcl_update", &Settings<T>::bcl_update)
8484
.def_readwrite("merit_function_type", &Settings<T>::merit_function_type)
8585
.def_readwrite("alpha_gpdal", &Settings<T>::alpha_gpdal)
86+
.def_readwrite("primal_infeasibility_solving",
87+
&Settings<T>::primal_infeasibility_solving)
88+
.def_readwrite("frequence_infeasibility_check",
89+
&Settings<T>::frequence_infeasibility_check)
8690
.def(pybind11::self == pybind11::self)
8791
.def(pybind11::self != pybind11::self)
8892
.def(pybind11::pickle(

‎bindings/python/src/expose-solve.hpp‎

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@ solveDenseQp(pybind11::module_ m)
4242
proxsuite::proxqp::InitialGuessStatus,
4343
bool,
4444
optional<T>,
45-
optional<T>>(&dense::solve<T>),
45+
optional<T>,
46+
bool>(&dense::solve<T>),
4647
"Function for solving a QP problem using PROXQP sparse backend directly "
4748
"without defining a QP object. It is possible to set up some of the solver "
4849
"parameters (warm start, initial guess option, proximal step sizes, "
@@ -98,7 +99,11 @@ solveDenseQp(pybind11::module_ m)
9899
pybind11::arg_v("eps_duality_gap_rel",
99100
nullopt,
100101
"relative accuracy threshold used for the duality-gap "
101-
"stopping criterion."));
102+
"stopping criterion."),
103+
pybind11::arg_v("primal_infeasibility_solving",
104+
false,
105+
"solves the closest feasible problem in L2 sense "
106+
"if the QP problem appears to be infeasible."));
102107

103108
m.def(
104109
"solve",
@@ -126,7 +131,8 @@ solveDenseQp(pybind11::module_ m)
126131
proxsuite::proxqp::InitialGuessStatus,
127132
bool,
128133
optional<T>,
129-
optional<T>>(&dense::solve<T>),
134+
optional<T>,
135+
bool>(&dense::solve<T>),
130136
"Function for solving a QP problem using PROXQP sparse backend directly "
131137
"without defining a QP object. It is possible to set up some of the solver "
132138
"parameters (warm start, initial guess option, proximal step sizes, "
@@ -184,7 +190,11 @@ solveDenseQp(pybind11::module_ m)
184190
pybind11::arg_v("eps_duality_gap_rel",
185191
nullopt,
186192
"relative accuracy threshold used for the duality-gap "
187-
"stopping criterion."));
193+
"stopping criterion."),
194+
pybind11::arg_v("primal_infeasibility_solving",
195+
false,
196+
"solves the closest feasible problem in L2 sense "
197+
"if the QP problem appears to be infeasible."));
188198
}
189199

190200
} // namespace python
@@ -255,7 +265,11 @@ solveSparseQp(pybind11::module_ m)
255265
pybind11::arg_v("eps_duality_gap_rel",
256266
nullopt,
257267
"relative accuracy threshold used for the duality-gap "
258-
"stopping criterion."));
268+
"stopping criterion."),
269+
pybind11::arg_v("primal_infeasibility_solving",
270+
false,
271+
"solves the closest feasible problem in L2 sense "
272+
"if the QP problem appears to be infeasible."));
259273
}
260274

261275
} // namespace python

‎doc/2-PROXQP_API/2-ProxQP_api.md‎

Lines changed: 47 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,21 @@ $$\begin{equation}\label{eq:approx_dg_sol}
7373

7474
where $[z]_+$ and $[z]_-$ stand for the projection of z onto the positive and negative orthant. ProxQP provides the ``check_duality_gap`` option to include this duality gap in the stopping criterion. Note that it is disabled by default, as other solvers don't check in general this criterion. Enable this option if you want a stronger guarantee that your solution is optimal. ProxQP will then check the same termination condition as SCS (for more details see, e.g., SCS's [optimality conditions checks](https://www.cvxgrp.org/scs/algorithm/index.html#optimality-conditions) as well as [section 7.2](https://doi.org/10.1137/20M1366307) in the corresponding paper). The absolute and relative thresholds $\epsilon^{\text{gap}}_{\text{abs}}, \epsilon^{\text{gap}}_{\text{rel}}$ for the duality gap can differ from those $\epsilon_{\text{abs}}, \epsilon_{\text{rel}}$ for residuals because, contrary to residuals which result from an infinite norm, the duality gap scales with the square root of the problem dimension (thus it is numerically harder to achieve a given duality gap for larger problems). A recommended choice is $\epsilon^{\text{gap}}_{\text{abs}} = \epsilon_{\text{abs}} \sqrt{\max(n, n_{\text{eq}}, n_{\text{ineq}})}$. Note finally that meeting all residual and duality-gap criteria can be difficult for ill-conditioned problems.
7575

76+
Finally, note that ProxQP has a specific feature for handling primal infeasibility. More precisely, if the problem appears to be primal infeasible, it will solve the closest primal feasible problem in $$\ell_2$$ sense, and (x,y,z) will satisfy
77+
78+
$$\begin{equation}\label{eq:approx_closest_qp_sol_rel}
79+
\begin{aligned}
80+
&\left\{
81+
\begin{array}{ll}
82+
\|Hx+g+A^Ty+C^Tz\|_{\infty} \leq \epsilon_{\text{abs}} + \epsilon_{\text{rel}}\max(\|Hx\|_{\infty},\|A^Ty\|_{\infty},\|C^Tz\|_{\infty},\|g\|_{\infty}), \\\
83+
\|A^\top(Ax-b)+C^\top[Cx-u]_+\|_{\infty} \leq \|A^\top 1 + C^\top 1\|_{\infty}*\epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(\|Ax\|_{\infty},\|b\|_{\infty}), \\\
84+
\end{array}
85+
\right.
86+
\end{aligned}
87+
\end{equation}$$
88+
89+
You can find more details on these subjects (and how to activate this feature with ProxQP) in the subsections describing the Settings and Results classes. You can also find more technical references with [this work](https://hal.science/hal-01057577).
90+
7691
\section OverviewAPIstructure ProxQP unified API for dense and sparse backends
7792

7893
ProxQP algorithm is implemented in two versions specialized for dense and sparse matrices. One simple and unified API has been designed for loading the dense and sparse backends. Concretely, it contains three methods:
@@ -395,7 +410,8 @@ In this table you have on the three columns from left to right: the name of the
395410
| safe_guard | 1.E4 | Safeguard parameter ensuring global convergence of the scheme. More precisely, if the total number of iteration is superior to safe_guard, the BCL scheme accept always the multipliers (hence the scheme is a pure proximal point algorithm).
396411
| preconditioner_max_iter | 10 | Maximal number of authorized iterations for the preconditioner.
397412
| preconditioner_accuracy | 1.E-3 | Accuracy level of the preconditioner.
398-
| HessianType | Dense | Defines the type of problem solved (Dense, Zero or Diagonal). In case Zero or Diagonal option are used, the solver optimize perform internally linear algebra operations involving the Hessian H.
413+
| HessianType | Dense | Defines the type of problem solved (Dense, Zero or Diagonal). In case Zero or Diagonal option are used, the solver optimize perform internally linear algebra operations involving the Hessian H.
414+
| primal_infeasibility_solving | False | If set to true, it solves the closest primal feasible problem if primal infeasibility is detected.
399415

400416
\subsection OverviewInitialGuess The different initial guesses
401417

@@ -464,6 +480,8 @@ The result subclass is composed of:
464480
* x: a primal solution,
465481
* y: a Lagrange optimal multiplier for equality constraints,
466482
* z: a Lagrange optimal multiplier for inequality constraints,
483+
* se: the optimal shift in $$\ell_2$$ with respect to equality constraints,
484+
* si: the optimal shift in $$\ell_2$$ with respect to inequality constraints,
467485
* info: a subclass which containts some information about the solver's execution.
468486

469487
If the solver has solved the problem, the triplet (x,y,z) satisfies:
@@ -481,6 +499,33 @@ $$\begin{equation}\label{eq:approx_qp_sol_rel}
481499
\end{equation}$$
482500
accordingly with the parameters eps_abs and eps_rel chosen by the user.
483501

502+
If the problem is primal infeasible and you have enable the solver to solve the closest feasible problem, then (x,y,z) will satisfies
503+
$$\begin{equation}\label{eq:approx_closest_qp_sol_rel}
504+
\begin{aligned}
505+
&\left\{
506+
\begin{array}{ll}
507+
\|Hx+g+A^Ty+C^Tz\|_{\infty} \leq \epsilon_{\text{abs}} + \epsilon_{\text{rel}}\max(\|Hx\|_{\infty},\|A^Ty\|_{\infty},\|C^Tz\|_{\infty},\|g\|_{\infty}), \\\
508+
\|A^\top(Ax-b)+C^\top[Cx-u]_+\|_{\infty} \leq \|A^\top 1 + C^\top 1\|_{\infty}*\epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(\|Ax\|_{\infty},\|b\|_{\infty}), \\\
509+
\end{array}
510+
\right.
511+
\end{aligned}
512+
\end{equation}$$
513+
514+
(se, si) stands in this context for the optimal shifts in $$\ell_2$$ sense which enables recovering a primal feasible problem. More precisely, they are derived such that
515+
516+
$$\begin{equation}\label{eq:QP_primal_feasible}\tag{QP_feas}
517+
\begin{aligned}
518+
\min_{x\in\mathbb{R}^{d}} & \quad \frac{1}{2}x^{T}Hx+g^{T}x \\\
519+
\text{s.t.}&\left\{
520+
\begin{array}{ll}
521+
Ax = b+se, \\\
522+
Cx \leq u+si, \\\
523+
\end{array}
524+
\right.
525+
\end{aligned}
526+
\end{equation}\\\
527+
defines a primal feasible problem.
528+
484529
Note that if you use the dense backend and its specific feature for handling box inequality constraints, then the first $$n_in$$ elements of z correspond to multipliers associated to the linear inequality formed with $$C$$ matrix, whereas the last $$d$$ elements correspond to multipliers associated to the box inequality constraints (see for example solve_dense_qp.cpp or solve_dense_qp.py).
485530
486531
\subsection OverviewInfoClass The info subclass
@@ -529,6 +574,7 @@ The solver has five status:
529574
* PROXQP_SOLVED: the problem is solved.
530575
* PROXQP_MAX_ITER_REACHED: the maximum number of iterations has been reached.
531576
* PROXQP_PRIMAL_INFEASIBLE: the problem is primal infeasible.
577+
* PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE: the closest feasible problem in L2 sense is solved.
532578
* PROXQP_DUAL_INFEASIBLE: the problem is dual infeasible.
533579
* PROXQP_NOT_RUN: the solver has not been run yet.
534580

‎include/proxsuite/proxqp/dense/helpers.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,7 @@ setup_equilibration(Workspace<T>& qpwork,
157157
qpsettings.preconditioner_accuracy,
158158
hessian_type,
159159
box_constraints,
160+
qpsettings.primal_infeasibility_solving,
160161
stack);
161162
qpwork.correction_guess_rhs_g = infty_norm(qpwork.g_scaled);
162163
}

‎include/proxsuite/proxqp/dense/linesearch.hpp‎

Lines changed: 9 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,7 @@ gpdal_derivative_results(const Model<T>& qpmodel,
8080
qpwork.primal_residual_in_scaled_up_plus_alphaCdx =
8181
qpwork.primal_residual_in_scaled_up + qpwork.Cdx * alpha;
8282
qpwork.primal_residual_in_scaled_low_plus_alphaCdx =
83-
qpwork.primal_residual_in_scaled_low + qpwork.Cdx * alpha;
83+
qpresults.si + qpwork.Cdx * alpha;
8484

8585
T a(qpwork.dw_aug.head(qpmodel.dim).dot(qpwork.Hdx) +
8686
qpresults.info.mu_eq_inv * (qpwork.Adx).squaredNorm() +
@@ -103,14 +103,13 @@ gpdal_derivative_results(const Model<T>& qpmodel,
103103
(qpwork.err.head(qpmodel.dim)).dot(qpwork.dw_aug.head(qpmodel.dim)) +
104104
qpresults.info.mu_eq_inv *
105105
(qpwork.Adx)
106-
.dot(qpwork.primal_residual_eq_scaled +
106+
.dot(qpresults.se +
107107
qpresults.y * qpresults.info.mu_eq)); // contains now: b =
108108
// dx.dot(H.dot(x) +
109109
// rho*(x-xe) + g) +
110110
// mu_eq_inv * Adx.dot(res_eq)
111111

112-
qpwork.rhs.segment(qpmodel.dim, qpmodel.n_eq) =
113-
qpwork.primal_residual_eq_scaled;
112+
qpwork.rhs.segment(qpmodel.dim, qpmodel.n_eq) = qpresults.se;
114113
b += qpresults.info.mu_eq_inv *
115114
qpwork.err.segment(qpmodel.dim, qpmodel.n_eq)
116115
.dot(qpwork.rhs.segment(
@@ -141,7 +140,7 @@ gpdal_derivative_results(const Model<T>& qpmodel,
141140
.select(qpwork.primal_residual_in_scaled_up,
142141
Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_constraints)) +
143142
(qpwork.primal_residual_in_scaled_low_plus_alphaCdx.array() < T(0.))
144-
.select(qpwork.primal_residual_in_scaled_low,
143+
.select(qpresults.si,
145144
Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_constraints));
146145

147146
b +=
@@ -209,7 +208,7 @@ primal_dual_derivative_results(const Model<T>& qpmodel,
209208
qpwork.primal_residual_in_scaled_up_plus_alphaCdx =
210209
qpwork.primal_residual_in_scaled_up + qpwork.Cdx * alpha;
211210
qpwork.primal_residual_in_scaled_low_plus_alphaCdx =
212-
qpwork.primal_residual_in_scaled_low + qpwork.Cdx * alpha;
211+
qpresults.si + qpwork.Cdx * alpha;
213212

214213
T a(qpwork.dw_aug.head(qpmodel.dim).dot(qpwork.Hdx) +
215214
qpresults.info.mu_eq_inv * (qpwork.Adx).squaredNorm() +
@@ -232,14 +231,13 @@ primal_dual_derivative_results(const Model<T>& qpmodel,
232231
(qpwork.err.head(qpmodel.dim)).dot(qpwork.dw_aug.head(qpmodel.dim)) +
233232
qpresults.info.mu_eq_inv *
234233
(qpwork.Adx)
235-
.dot(qpwork.primal_residual_eq_scaled +
234+
.dot(qpresults.se +
236235
qpresults.y * qpresults.info.mu_eq)); // contains now: b =
237236
// dx.dot(H.dot(x) +
238237
// rho*(x-xe) + g) +
239238
// mu_eq_inv * Adx.dot(res_eq)
240239

241-
qpwork.rhs.segment(qpmodel.dim, qpmodel.n_eq) =
242-
qpwork.primal_residual_eq_scaled;
240+
qpwork.rhs.segment(qpmodel.dim, qpmodel.n_eq) = qpresults.se;
243241
b += qpresults.info.nu * qpresults.info.mu_eq_inv *
244242
qpwork.err.segment(qpmodel.dim, qpmodel.n_eq)
245243
.dot(qpwork.rhs.segment(
@@ -268,7 +266,7 @@ primal_dual_derivative_results(const Model<T>& qpmodel,
268266
.select(qpwork.primal_residual_in_scaled_up,
269267
Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_constraints)) +
270268
(qpwork.primal_residual_in_scaled_low_plus_alphaCdx.array() < T(0.))
271-
.select(qpwork.primal_residual_in_scaled_low,
269+
.select(qpresults.si,
272270
Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_constraints));
273271

274272
b += qpresults.info.mu_in_inv *
@@ -377,8 +375,7 @@ primal_dual_ls(const Model<T>& qpmodel,
377375
if (alpha_ > machine_eps) {
378376
qpwork.alphas.push(alpha_);
379377
}
380-
alpha_ = -qpwork.primal_residual_in_scaled_low(i) /
381-
(qpwork.Cdx(i) + machine_eps);
378+
alpha_ = -qpresults.si(i) / (qpwork.Cdx(i) + machine_eps);
382379
if (alpha_ > machine_eps) {
383380
qpwork.alphas.push(alpha_);
384381
}
@@ -526,7 +523,6 @@ primal_dual_ls(const Model<T>& qpmodel,
526523
* [last_alpha_neg,first_alpha_pos] and can be computed exactly as phi'
527524
* is an affine function in alpha
528525
*/
529-
530526
qpwork.alpha = std::abs(alpha_last_neg -
531527
last_neg_grad * (alpha_first_pos - alpha_last_neg) /
532528
(first_pos_grad - last_neg_grad));

‎include/proxsuite/proxqp/dense/preconditioner/ruiz.hpp‎

Lines changed: 34 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ ruiz_scale_qp_in_place( //
3737
Symmetry sym,
3838
HessianType HessianType,
3939
const bool box_constraints,
40+
bool preconditioning_for_infeasible_problems,
4041
proxsuite::linalg::veg::dynstack::DynStackMut stack) -> T
4142
{
4243
T c(1);
@@ -165,25 +166,29 @@ ruiz_scale_qp_in_place( //
165166
}
166167
break;
167168
}
168-
for (isize k = 0; k < n_eq; ++k) {
169-
T aux = sqrt(infty_norm(A.row(k)));
170-
if (aux == T(0)) {
171-
delta(n + k) = T(1);
172-
} else {
173-
delta(n + k) = T(1) / (aux + machine_eps);
169+
if (preconditioning_for_infeasible_problems) {
170+
delta.tail(n_eq + n_constraints).setOnes();
171+
} else {
172+
for (isize k = 0; k < n_eq; ++k) {
173+
T aux = sqrt(infty_norm(A.row(k)));
174+
if (aux == T(0)) {
175+
delta(n + k) = T(1);
176+
} else {
177+
delta(n + k) = T(1) / (aux + machine_eps);
178+
}
174179
}
175-
}
176-
for (isize k = 0; k < n_in; ++k) {
177-
T aux = sqrt(infty_norm(C.row(k)));
178-
if (aux == T(0)) {
179-
delta(k + n + n_eq) = T(1);
180-
} else {
181-
delta(k + n + n_eq) = T(1) / (aux + machine_eps);
180+
for (isize k = 0; k < n_in; ++k) {
181+
T aux = sqrt(infty_norm(C.row(k)));
182+
if (aux == T(0)) {
183+
delta(k + n + n_eq) = T(1);
184+
} else {
185+
delta(k + n + n_eq) = T(1) / (aux + machine_eps);
186+
}
182187
}
183-
}
184-
if (box_constraints) {
185-
for (isize k = 0; k < n; ++k) {
186-
delta(k + n + n_eq + n_in) = T(1) / sqrt(i_scaled[k] + machine_eps);
188+
if (box_constraints) {
189+
for (isize k = 0; k < n; ++k) {
190+
delta(k + n + n_eq + n_in) = T(1) / sqrt(i_scaled[k] + machine_eps);
191+
}
187192
}
188193
}
189194
}
@@ -394,19 +399,22 @@ struct RuizEquilibration
394399
const T epsilon,
395400
const HessianType& HessianType,
396401
const bool box_constraints,
402+
const bool preconditioning_for_infeasible_problems,
397403
proxsuite::linalg::veg::dynstack::DynStackMut stack)
398404
{
399405
if (execute_preconditioner) {
400406
delta.setOnes();
401-
c = detail::ruiz_scale_qp_in_place({ proxqp::from_eigen, delta },
402-
logger_ptr,
403-
qp,
404-
epsilon,
405-
max_iter,
406-
sym,
407-
HessianType,
408-
box_constraints,
409-
stack);
407+
c =
408+
detail::ruiz_scale_qp_in_place({ proxqp::from_eigen, delta },
409+
logger_ptr,
410+
qp,
411+
epsilon,
412+
max_iter,
413+
sym,
414+
HessianType,
415+
box_constraints,
416+
preconditioning_for_infeasible_problems,
417+
stack);
410418
} else {
411419

412420
auto H = qp.H.to_eigen();

0 commit comments

Comments
 (0)