Skip to content

Commit 3d2001b

Browse files
authored
Merge pull request #169 from stephane-caron/feature/set_dual_eps
Add duality-gap thresholds
2 parents 0d3629d + ef80075 commit 3d2001b

File tree

13 files changed

+85
-32
lines changed

13 files changed

+85
-32
lines changed

bindings/python/src/expose-settings.hpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,8 @@ exposeSettings(pybind11::module_ m)
7171
&Settings<T>::compute_preconditioner)
7272
.def_readwrite("update_preconditioner", &Settings<T>::update_preconditioner)
7373
.def_readwrite("check_duality_gap", &Settings<T>::check_duality_gap)
74+
.def_readwrite("eps_duality_gap_abs", &Settings<T>::eps_duality_gap_abs)
75+
.def_readwrite("eps_duality_gap_rel", &Settings<T>::eps_duality_gap_rel)
7476
.def_readwrite("verbose", &Settings<T>::verbose)
7577
.def_readwrite("bcl_update", &Settings<T>::bcl_update)
7678
.def(pybind11::self == pybind11::self)

bindings/python/src/expose-solve.hpp

Lines changed: 21 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,9 @@ solveDenseQp(pybind11::module_ m)
4040
bool,
4141
optional<isize>,
4242
proxsuite::proxqp::InitialGuessStatus,
43-
bool>(&dense::solve<T>),
43+
bool,
44+
optional<T>,
45+
optional<T>>(&dense::solve<T>),
4446
"Function for solving a QP problem using PROXQP sparse backend directly "
4547
"without defining a QP object. It is possible to set up some of the solver "
4648
"parameters (warm start, initial guess option, proximal step sizes, "
@@ -88,7 +90,15 @@ solveDenseQp(pybind11::module_ m)
8890
"check_duality_gap",
8991
false,
9092
"if set to true, include the duality gap in absolute and relative "
91-
"stopping criteria."));
93+
"stopping criteria."),
94+
pybind11::arg_v("eps_duality_gap_abs",
95+
nullopt,
96+
"absolute accuracy threshold used for the duality-gap "
97+
"stopping criterion."),
98+
pybind11::arg_v("eps_duality_gap_rel",
99+
nullopt,
100+
"relative accuracy threshold used for the duality-gap "
101+
"stopping criterion."));
92102
}
93103

94104
} // namespace python
@@ -151,7 +161,15 @@ solveSparseQp(pybind11::module_ m)
151161
pybind11::arg_v("check_duality_gap",
152162
false,
153163
"if set to true, include the duality gap in absolute and "
154-
"relative stopping criteria."));
164+
"relative stopping criteria."),
165+
pybind11::arg_v("eps_duality_gap_abs",
166+
nullopt,
167+
"absolute accuracy threshold used for the duality-gap "
168+
"stopping criterion."),
169+
pybind11::arg_v("eps_duality_gap_rel",
170+
nullopt,
171+
"relative accuracy threshold used for the duality-gap "
172+
"stopping criterion."));
155173
}
156174

157175
} // namespace python

doc/2-PROXQP_API/2-ProxQP_api.md

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -67,11 +67,11 @@ It is important to note that this stopping criterion on primal and dual residual
6767

6868
$$\begin{equation}\label{eq:approx_dg_sol}
6969
\begin{aligned}
70-
r_g := | x^T H x + g^T x + b^T y + u^T [z]_+ + l^T [z]_- | \leq \epsilon_{\text{abs}} + \epsilon_{\text{rel}} \max(\|x^T H x\|, \|g^T x\|, \|b^T y\|, \|u^T [z]_+\|, \|l^T [z]_-\|), \\
70+
r_g := | x^T H x + g^T x + b^T y + u^T [z]_+ + l^T [z]_- | \leq \epsilon^{\text{gap}}_{\text{abs}} + \epsilon^{\text{gap}}_{\text{rel}} \max(\|x^T H x\|, \|g^T x\|, \|b^T y\|, \|u^T [z]_+\|, \|l^T [z]_-\|), \\
7171
\end{aligned}
7272
\end{equation}$$
7373

74-
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). Note finally that meeting all these criteria can be difficult for ill-conditioned problems.
74+
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

7676
\section OverviewAPIstructure ProxQP unified API for dense and sparse backends
7777

@@ -312,6 +312,8 @@ In this table you have on the three columns from left to right: the name of the
312312
| eps_abs | 1.E-5 | Asbolute stopping criterion of the solver.
313313
| eps_rel | 0 | Relative stopping criterion of the solver.
314314
| check_duality_gap | False | If set to true, include the duality gap in absolute and relative stopping criteria.
315+
| eps_duality_gap_abs | 1.E-4 | Asbolute duality-gap stopping criterion of the solver.
316+
| eps_duality_gap_rel | 0 | Relative duality-gap stopping criterion of the solver.
315317
| VERBOSE | False | If set to true, the solver prints information at each loop.
316318
| default_rho | 1.E-6 | Default rho parameter of result class (i.e., for each initial guess, except WARM_START_WITH_PREVIOUS_RESULT, after a new solve or update, the solver initializes rho to this value).
317319
| default_mu_eq | 1.E-3 | Default mu_eq parameter of result class (i.e., for each initial guess, except WARM_START_WITH_PREVIOUS_RESULT, after a new solve or update, the solver initializes mu_eq to this value).

doc/3-ProxQP_solve.md

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -67,11 +67,11 @@ It is important to note that this stopping criterion on primal and dual residual
6767

6868
$$\begin{equation}\label{eq:approx_dg_sol}
6969
\begin{aligned}
70-
r_g := | x^T H x + g^T x + b^T y + u^T [z]_+ + l^T [z]_- | \leq \epsilon_{\text{abs}} + \epsilon_{\text{rel}} \max(\|x^T H x\|, \|g^T x\|, \|b^T y\|, \|u^T [z]_+\|, \|l^T [z]_-\|), \\
70+
r_g := | x^T H x + g^T x + b^T y + u^T [z]_+ + l^T [z]_- | \leq \epsilon^{\text{gap}}_{\text{abs}} + \epsilon^{\text{gap}}_{\text{rel}} \max(\|x^T H x\|, \|g^T x\|, \|b^T y\|, \|u^T [z]_+\|, \|l^T [z]_-\|), \\
7171
\end{aligned}
7272
\end{equation}$$
7373

74-
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). Note finally that meeting all these criteria can be difficult for ill-conditioned problems.
74+
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

7676
\section OverviewAsingleSolveFunction A single solve function for dense and sparse backends
7777

@@ -106,6 +106,8 @@ Different options are available for the solve function. In the table below you h
106106
| eps_abs | 1.E-5 | Asbolute stopping criterion of the solver.
107107
| eps_rel | 0 | Relative stopping criterion of the solver.
108108
| check_duality_gap | False | If set to true, include the duality gap in absolute and relative stopping criteria.
109+
| eps_duality_gap_abs | 1.E-4 | Asbolute duality-gap stopping criterion of the solver.
110+
| eps_duality_gap_rel | 0 | Relative duality-gap stopping criterion of the solver.
109111
| mu_eq | 1.E-3 | Proximal step size wrt equality constraints multiplier.
110112
| mu_in | 1.E-1 | Proximal step size wrt inequality constraints multiplier.
111113
| rho | 1.E-6 | Proximal step size wrt primal variable.
@@ -157,4 +159,4 @@ Note that if some elements of your QP model are not defined (for example a QP wi
157159
\include initializing_with_none_without_api.py
158160
</td>
159161
</tr>
160-
</table>
162+
</table>

include/proxsuite/proxqp/dense/solver.hpp

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1075,8 +1075,9 @@ qp_solve( //
10751075
}
10761076
if (is_primal_feasible && is_dual_feasible) {
10771077
if (qpsettings.check_duality_gap) {
1078-
if (std::abs(qpresults.info.duality_gap) <=
1079-
qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
1078+
if (std::fabs(qpresults.info.duality_gap) <=
1079+
qpsettings.eps_duality_gap_abs +
1080+
qpsettings.eps_duality_gap_rel * rhs_duality_gap) {
10801081
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
10811082
break;
10821083
}
@@ -1161,8 +1162,9 @@ qp_solve( //
11611162

11621163
if (is_dual_feasible) {
11631164
if (qpsettings.check_duality_gap) {
1164-
if (std::abs(qpresults.info.duality_gap) <=
1165-
qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
1165+
if (std::fabs(qpresults.info.duality_gap) <=
1166+
qpsettings.eps_duality_gap_abs +
1167+
qpsettings.eps_duality_gap_rel * rhs_duality_gap) {
11661168
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
11671169
break;
11681170
}

include/proxsuite/proxqp/dense/utils.hpp

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -347,9 +347,6 @@ global_dual_residual(Results<T>& qpresults,
347347
// dual_feasibility_rhs_3 = norm(unscaled(CTz))
348348
//
349349
// dual_residual_scaled = scaled(Hx + g + ATy + CTz)
350-
const isize max_dim =
351-
std::max(qpmodel.dim, std::max(qpmodel.n_eq, qpmodel.n_in));
352-
const T sqrt_max_dim(std::sqrt(max_dim)); // for normalizing scalar products
353350

354351
qpwork.dual_residual_scaled = qpwork.g_scaled;
355352
qpwork.CTz.noalias() =
@@ -361,7 +358,7 @@ global_dual_residual(Results<T>& qpresults,
361358

362359
ruiz.unscale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
363360
duality_gap = (qpmodel.g).dot(qpresults.x);
364-
rhs_duality_gap = std::abs(duality_gap);
361+
rhs_duality_gap = std::fabs(duality_gap);
365362
const T xHx = (qpwork.CTz).dot(qpresults.x);
366363
duality_gap += xHx; // contains now xHx+g.Tx
367364
rhs_duality_gap = std::max(rhs_duality_gap, std::abs(xHx));
@@ -409,9 +406,6 @@ global_dual_residual(Results<T>& qpresults,
409406
duality_gap += zl;
410407

411408
ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
412-
413-
duality_gap /= sqrt_max_dim; // in order to get an a-dimensional duality gap
414-
rhs_duality_gap /= sqrt_max_dim;
415409
}
416410

417411
} // namespace dense

include/proxsuite/proxqp/dense/wrapper.hpp

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -393,6 +393,10 @@ struct QP
393393
* initial iterate values.
394394
* @param check_duality_gap If set to true, include the duality gap in absolute
395395
* and relative stopping criteria.
396+
* @param eps_duality_gap_abs absolute accuracy threshold for the duality-gap
397+
* criterion.
398+
* @param eps_duality_gap_rel relative accuracy threshold for the duality-gap
399+
* criterion.
396400
*/
397401
template<typename T>
398402
proxqp::Results<T>
@@ -418,7 +422,9 @@ solve(
418422
optional<isize> max_iter = nullopt,
419423
proxsuite::proxqp::InitialGuessStatus initial_guess =
420424
proxsuite::proxqp::InitialGuessStatus::EQUALITY_CONSTRAINED_INITIAL_GUESS,
421-
bool check_duality_gap = false)
425+
bool check_duality_gap = false,
426+
optional<T> eps_duality_gap_abs = nullopt,
427+
optional<T> eps_duality_gap_rel = nullopt)
422428
{
423429
isize n(0);
424430
isize n_eq(0);
@@ -449,6 +455,12 @@ solve(
449455
if (max_iter != nullopt) {
450456
Qp.settings.max_iter = max_iter.value();
451457
}
458+
if (eps_duality_gap_abs != nullopt) {
459+
Qp.settings.eps_duality_gap_abs = eps_duality_gap_abs.value();
460+
}
461+
if (eps_duality_gap_rel != nullopt) {
462+
Qp.settings.eps_duality_gap_rel = eps_duality_gap_rel.value();
463+
}
452464
Qp.settings.compute_timings = compute_timings;
453465
Qp.init(H, g, A, b, C, l, u, compute_preconditioner, rho, mu_eq, mu_in);
454466
Qp.solve(x, y, z);

include/proxsuite/proxqp/settings.hpp

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,10 @@ struct Settings
8484
bool update_preconditioner;
8585
bool compute_preconditioner;
8686
bool compute_timings;
87+
8788
bool check_duality_gap;
89+
T eps_duality_gap_abs;
90+
T eps_duality_gap_rel;
8891

8992
isize preconditioner_max_iter;
9093
T preconditioner_accuracy;
@@ -139,6 +142,8 @@ struct Settings
139142
* time).
140143
* @param check_duality_gap If set to true, duality gap will be calculated and
141144
* included in the stopping criterion.
145+
* @param eps_duality_gap_abs absolute duality-gap stopping criterion.
146+
* @param eps_duality_gap_rel relative duality-gap stopping criterion.
142147
* @param preconditioner_max_iter maximal number of authorized iterations for
143148
* the preconditioner.
144149
* @param preconditioner_accuracy accuracy level of the preconditioner.
@@ -187,6 +192,8 @@ struct Settings
187192
bool compute_preconditioner = true,
188193
bool compute_timings = false,
189194
bool check_duality_gap = false,
195+
T eps_duality_gap_abs = 1.e-4,
196+
T eps_duality_gap_rel = 0,
190197
isize preconditioner_max_iter = 10,
191198
T preconditioner_accuracy = 1.e-3,
192199
T eps_primal_inf = 1.E-4,
@@ -223,6 +230,8 @@ struct Settings
223230
, compute_preconditioner(compute_preconditioner)
224231
, compute_timings(compute_timings)
225232
, check_duality_gap(check_duality_gap)
233+
, eps_duality_gap_abs(eps_duality_gap_abs)
234+
, eps_duality_gap_rel(eps_duality_gap_rel)
226235
, preconditioner_max_iter(preconditioner_max_iter)
227236
, preconditioner_accuracy(preconditioner_accuracy)
228237
, eps_primal_inf(eps_primal_inf)
@@ -269,6 +278,8 @@ operator==(const Settings<T>& settings1, const Settings<T>& settings2)
269278
settings1.compute_preconditioner == settings2.compute_preconditioner &&
270279
settings1.compute_timings == settings2.compute_timings &&
271280
settings1.check_duality_gap == settings2.check_duality_gap &&
281+
settings1.eps_duality_gap_abs == settings2.eps_duality_gap_abs &&
282+
settings1.eps_duality_gap_rel == settings2.eps_duality_gap_rel &&
272283
settings1.preconditioner_max_iter == settings2.preconditioner_max_iter &&
273284
settings1.preconditioner_accuracy == settings2.preconditioner_accuracy &&
274285
settings1.eps_primal_inf == settings2.eps_primal_inf &&

include/proxsuite/proxqp/sparse/solver.hpp

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -843,8 +843,9 @@ qp_solve(Results<T>& results,
843843
if (is_primal_feasible(primal_feasibility_lhs) &&
844844
is_dual_feasible(dual_feasibility_lhs)) {
845845
if (settings.check_duality_gap) {
846-
if (std::abs(results.info.duality_gap) <=
847-
settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
846+
if (std::fabs(results.info.duality_gap) <=
847+
settings.eps_duality_gap_abs +
848+
settings.eps_duality_gap_rel * rhs_duality_gap) {
848849
results.info.pri_res = primal_feasibility_lhs;
849850
results.info.dua_res = dual_feasibility_lhs;
850851
results.info.status = QPSolverOutput::PROXQP_SOLVED;
@@ -1244,7 +1245,8 @@ qp_solve(Results<T>& results,
12441245
is_dual_feasible(dual_feasibility_lhs_new)) {
12451246
if (settings.check_duality_gap) {
12461247
if (std::fabs(results.info.duality_gap) <=
1247-
settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
1248+
settings.eps_duality_gap_abs +
1249+
settings.eps_duality_gap_rel * rhs_duality_gap) {
12481250
results.info.pri_res = primal_feasibility_lhs_new;
12491251
results.info.dua_res = dual_feasibility_lhs_new;
12501252
results.info.status = QPSolverOutput::PROXQP_SOLVED;

include/proxsuite/proxqp/sparse/utils.hpp

Lines changed: 1 addition & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -637,9 +637,6 @@ unscaled_primal_dual_residual(
637637
{
638638
isize n = x_e.rows();
639639

640-
const isize max_dim = std::max(data.dim, std::max(data.n_eq, data.n_in));
641-
const T sqrt_max_dim(std::sqrt(max_dim)); // for normalizing scalar products
642-
643640
LDLT_TEMP_VEC_UNINIT(T, tmp, n, stack);
644641
dual_residual_scaled = qp_scaled.g.to_eigen();
645642
{
@@ -652,7 +649,7 @@ unscaled_primal_dual_residual(
652649
dual_feasibility_rhs_0 = infty_norm(tmp);
653650
precond.unscale_primal_in_place({ proxqp::from_eigen, x_e });
654651
results.info.duality_gap = x_e.dot(data.g); // contains gTx
655-
rhs_duality_gap = std::abs(results.info.duality_gap);
652+
rhs_duality_gap = std::fabs(results.info.duality_gap);
656653

657654
const T xHx = (tmp).dot(x_e);
658655
results.info.duality_gap += xHx;
@@ -681,10 +678,6 @@ unscaled_primal_dual_residual(
681678
rhs_duality_gap = std::max(rhs_duality_gap, std::abs(zu));
682679

683680
precond.scale_dual_in_place_in({ proxsuite::proxqp::from_eigen, z_e });
684-
685-
results.info.duality_gap /=
686-
sqrt_max_dim; // in order to get an a-dimensional duality gap
687-
rhs_duality_gap /= sqrt_max_dim;
688681
}
689682

690683
{

0 commit comments

Comments
 (0)