Skip to content

Commit 6244f87

Browse files
Bambadejcarpent
Bambade
authored and
jcarpent
committed
add duality gap measure
1 parent d22f612 commit 6244f87

File tree

5 files changed

+152
-39
lines changed

5 files changed

+152
-39
lines changed

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

Lines changed: 36 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -985,6 +985,12 @@ qp_solve( //
985985
T primal_feasibility_in_lhs(0);
986986
T dual_feasibility_lhs(0);
987987

988+
T duality_gap(0);
989+
T rhs_duality_gap(0);
990+
isize max_dim = std::max(qpmodel.dim,qpmodel.n_eq);
991+
max_dim = std::max(max_dim,qpmodel.n_in);
992+
T sqrt_max_dim(std::sqrt(max_dim));
993+
988994
for (i64 iter = 0; iter < qpsettings.max_iter; ++iter) {
989995

990996
// compute primal residual
@@ -1002,13 +1008,19 @@ qp_solve( //
10021008

10031009
global_dual_residual(qpresults,
10041010
qpwork,
1011+
qpmodel,
10051012
ruiz,
10061013
dual_feasibility_lhs,
10071014
dual_feasibility_rhs_0,
10081015
dual_feasibility_rhs_1,
1009-
dual_feasibility_rhs_3);
1016+
dual_feasibility_rhs_3,
1017+
rhs_duality_gap,
1018+
duality_gap,
1019+
sqrt_max_dim
1020+
);
10101021
qpresults.info.pri_res = primal_feasibility_lhs;
10111022
qpresults.info.dua_res = dual_feasibility_lhs;
1023+
qpresults.info.duality_gap = duality_gap;
10121024

10131025
T new_bcl_mu_in(qpresults.info.mu_in);
10141026
T new_bcl_mu_eq(qpresults.info.mu_eq);
@@ -1059,8 +1071,9 @@ qp_solve( //
10591071
std::cout << std::scientific << std::setw(2) << std::setprecision(2)
10601072
<< "| primal residual=" << qpresults.info.pri_res
10611073
<< "| dual residual=" << qpresults.info.dua_res
1062-
<< " | mu_in=" << qpresults.info.mu_in
1063-
<< " | rho=" << qpresults.info.rho << std::endl;
1074+
<< "| duality gap=" << qpresults.info.duality_gap
1075+
<< "| mu_in=" << qpresults.info.mu_in
1076+
<< "| rho=" << qpresults.info.rho << std::endl;
10641077
ruiz.scale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
10651078
ruiz.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
10661079
ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
@@ -1079,8 +1092,10 @@ qp_solve( //
10791092
qpresults.info.rho = rho_new;
10801093
}
10811094
if (is_dual_feasible) {
1082-
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
1083-
break;
1095+
if (qpresults.info.duality_gap<=qpsettings.eps_abs* rhs_duality_gap ){
1096+
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
1097+
break;
1098+
}
10841099
}
10851100
}
10861101
qpresults.info.iter_ext += 1; // We start a new external loop update
@@ -1138,12 +1153,17 @@ qp_solve( //
11381153

11391154
global_dual_residual(qpresults,
11401155
qpwork,
1156+
qpmodel,
11411157
ruiz,
11421158
dual_feasibility_lhs_new,
11431159
dual_feasibility_rhs_0,
11441160
dual_feasibility_rhs_1,
1145-
dual_feasibility_rhs_3);
1161+
dual_feasibility_rhs_3,
1162+
rhs_duality_gap,
1163+
duality_gap,
1164+
sqrt_max_dim);
11461165
qpresults.info.dua_res = dual_feasibility_lhs_new;
1166+
qpresults.info.duality_gap = duality_gap;
11471167

11481168
is_dual_feasible =
11491169
dual_feasibility_lhs_new <=
@@ -1154,8 +1174,10 @@ qp_solve( //
11541174
std::max(dual_feasibility_rhs_1, qpwork.dual_feasibility_rhs_2)));
11551175

11561176
if (is_dual_feasible) {
1157-
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
1158-
break;
1177+
if (qpresults.info.duality_gap<=qpsettings.eps_abs * rhs_duality_gap ){
1178+
qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
1179+
break;
1180+
}
11591181
}
11601182
}
11611183
if (qpsettings.bcl_update) {
@@ -1190,12 +1212,17 @@ qp_solve( //
11901212

11911213
global_dual_residual(qpresults,
11921214
qpwork,
1215+
qpmodel,
11931216
ruiz,
11941217
dual_feasibility_lhs_new,
11951218
dual_feasibility_rhs_0,
11961219
dual_feasibility_rhs_1,
1197-
dual_feasibility_rhs_3);
1220+
dual_feasibility_rhs_3,
1221+
rhs_duality_gap,
1222+
duality_gap,
1223+
sqrt_max_dim);
11981224
qpresults.info.dua_res = dual_feasibility_lhs_new;
1225+
qpresults.info.duality_gap = duality_gap;
11991226

12001227
if (primal_feasibility_lhs_new >= primal_feasibility_lhs &&
12011228
dual_feasibility_lhs_new >= dual_feasibility_lhs &&

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

Lines changed: 35 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -345,11 +345,15 @@ template<typename T>
345345
void
346346
global_dual_residual(const Results<T>& qpresults,
347347
Workspace<T>& qpwork,
348+
const Model<T>& qpmodel,
348349
const preconditioner::RuizEquilibration<T>& ruiz,
349350
T& dual_feasibility_lhs,
350351
T& dual_feasibility_rhs_0,
351352
T& dual_feasibility_rhs_1,
352-
T& dual_feasibility_rhs_3)
353+
T& dual_feasibility_rhs_3,
354+
T& rhs_duality_gap,
355+
T& duality_gap,
356+
T sqrt_max_dim)
353357
{
354358
// dual_feasibility_lhs = norm(dual_residual_scaled)
355359
// dual_feasibility_rhs_0 = norm(unscaled(Hx))
@@ -362,8 +366,18 @@ global_dual_residual(const Results<T>& qpresults,
362366
qpwork.H_scaled.template selfadjointView<Eigen::Lower>() * qpresults.x;
363367
qpwork.dual_residual_scaled += qpwork.CTz;
364368
ruiz.unscale_dual_residual_in_place(
365-
VectorViewMut<T>{ from_eigen, qpwork.CTz });
369+
VectorViewMut<T>{ from_eigen, qpwork.CTz }); // contains unscaled Hx
366370
dual_feasibility_rhs_0 = infty_norm(qpwork.CTz);
371+
372+
ruiz.unscale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
373+
duality_gap = (qpmodel.g).dot(qpresults.x);
374+
rhs_duality_gap = std::abs(duality_gap);
375+
T xHx = (qpwork.CTz).dot(qpresults.x) ;
376+
duality_gap += xHx;// contains now 0.5*xHx+g.Tx
377+
rhs_duality_gap = std::max(rhs_duality_gap,std::abs(xHx));
378+
379+
ruiz.scale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
380+
367381
qpwork.CTz.noalias() = qpwork.A_scaled.transpose() * qpresults.y;
368382
qpwork.dual_residual_scaled += qpwork.CTz;
369383
ruiz.unscale_dual_residual_in_place(
@@ -383,6 +397,25 @@ global_dual_residual(const Results<T>& qpresults,
383397

384398
ruiz.scale_dual_residual_in_place(
385399
VectorViewMut<T>{ from_eigen, qpwork.dual_residual_scaled });
400+
401+
ruiz.unscale_dual_in_place_eq(
402+
VectorViewMut<T>{ from_eigen, qpresults.y });
403+
T by = (qpmodel.b).dot(qpresults.y);
404+
rhs_duality_gap = std::max(rhs_duality_gap,std::abs(by));
405+
duality_gap += by;
406+
ruiz.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
407+
ruiz.unscale_dual_in_place_in(
408+
VectorViewMut<T>{ from_eigen, qpresults.z });
409+
T zu= positive_part(qpresults.z).dot(qpmodel.u);
410+
rhs_duality_gap = std::max(rhs_duality_gap,std::abs(zu));
411+
duality_gap += zu;
412+
T zl = negative_part(qpresults.z).dot(qpmodel.l);
413+
rhs_duality_gap = std::max(rhs_duality_gap,std::abs(zl));
414+
duality_gap += zl;
415+
ruiz.scale_dual_in_place_in(
416+
VectorViewMut<T>{ from_eigen, qpresults.z });
417+
duality_gap /= sqrt_max_dim; // in order to get an a-dimensional duality gap
418+
rhs_duality_gap +=1.;
386419
}
387420

388421
} // namespace dense

‎include/proxsuite/proxqp/results.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ struct Info
4949
T objValue;
5050
T pri_res;
5151
T dua_res;
52-
52+
T duality_gap;
5353
//// sparse backend used by solver, either CholeskySparse or MatrixFree
5454
SparseBackend sparse_backend;
5555
};
@@ -106,6 +106,7 @@ struct Results
106106
info.objValue = 0.;
107107
info.pri_res = 0.;
108108
info.dua_res = 0.;
109+
info.duality_gap = 0.;
109110
info.status = QPSolverOutput::PROXQP_NOT_RUN;
110111
info.sparse_backend = SparseBackend::Automatic;
111112
}
@@ -132,6 +133,7 @@ struct Results
132133
info.rho_updates = 0;
133134
info.pri_res = 0.;
134135
info.dua_res = 0.;
136+
info.duality_gap = 0.;
135137
info.status = QPSolverOutput::PROXQP_MAX_ITER_REACHED;
136138
info.sparse_backend = SparseBackend::Automatic;
137139
}

‎include/proxsuite/proxqp/sparse/solver.hpp‎

Lines changed: 41 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -724,6 +724,11 @@ qp_solve(Results<T>& results,
724724
break;
725725
}
726726
}
727+
T rhs_duality_gap(0);
728+
isize max_dim = std::max(data.dim,data.n_eq);
729+
max_dim = std::max(max_dim,data.n_in);
730+
T sqrt_max_dim(std::sqrt(max_dim));
731+
727732
for (isize iter = 0; iter < settings.max_iter; ++iter) {
728733

729734
results.info.iter_ext += 1;
@@ -777,7 +782,8 @@ qp_solve(Results<T>& results,
777782
VEG_BIND( // ?
778783
auto,
779784
(primal_feasibility_lhs, dual_feasibility_lhs),
780-
detail::unscaled_primal_dual_residual(primal_residual_eq_scaled,
785+
detail::unscaled_primal_dual_residual(results,
786+
primal_residual_eq_scaled,
781787
primal_residual_in_scaled_lo,
782788
primal_residual_in_scaled_up,
783789
dual_residual_scaled,
@@ -786,12 +792,14 @@ qp_solve(Results<T>& results,
786792
dual_feasibility_rhs_0,
787793
dual_feasibility_rhs_1,
788794
dual_feasibility_rhs_3,
795+
rhs_duality_gap,
796+
sqrt_max_dim,
789797
precond,
790798
data,
791799
qp_scaled.as_const(),
792-
detail::vec(x_e),
793-
detail::vec(y_e),
794-
detail::vec(z_e),
800+
detail::vec_mut(x_e),
801+
detail::vec_mut(y_e),
802+
detail::vec_mut(z_e),
795803
stack));
796804
/*put in debug mode
797805
if (settings.verbose) {
@@ -824,8 +832,9 @@ qp_solve(Results<T>& results,
824832
std::cout << std::scientific << std::setw(2) << std::setprecision(2)
825833
<< "| primal residual=" << primal_feasibility_lhs
826834
<< "| dual residual=" << dual_feasibility_lhs
827-
<< " | mu_in=" << results.info.mu_in
828-
<< " | rho=" << results.info.rho << std::endl;
835+
<< "| dualality gap=" << results.info.duality_gap
836+
<< "| mu_in=" << results.info.mu_in
837+
<< "| rho=" << results.info.rho << std::endl;
829838
results.info.pri_res = primal_feasibility_lhs;
830839
results.info.dua_res = dual_feasibility_lhs;
831840
precond.scale_primal_in_place(VectorViewMut<T>{ from_eigen, x_e });
@@ -834,10 +843,12 @@ qp_solve(Results<T>& results,
834843
}
835844
if (is_primal_feasible(primal_feasibility_lhs) &&
836845
is_dual_feasible(dual_feasibility_lhs)) {
837-
results.info.pri_res = primal_feasibility_lhs;
838-
results.info.dua_res = dual_feasibility_lhs;
839-
results.info.status = QPSolverOutput::PROXQP_SOLVED;
840-
break;
846+
if (results.info.duality_gap<=settings.eps_abs*rhs_duality_gap){
847+
results.info.pri_res = primal_feasibility_lhs;
848+
results.info.dua_res = dual_feasibility_lhs;
849+
results.info.status = QPSolverOutput::PROXQP_SOLVED;
850+
break;
851+
}
841852
}
842853

843854
LDLT_TEMP_VEC_UNINIT(T, x_prev_e, n, stack);
@@ -1202,7 +1213,8 @@ qp_solve(Results<T>& results,
12021213
VEG_BIND(
12031214
auto,
12041215
(primal_feasibility_lhs_new, dual_feasibility_lhs_new),
1205-
detail::unscaled_primal_dual_residual(primal_residual_eq_scaled,
1216+
detail::unscaled_primal_dual_residual(results,
1217+
primal_residual_eq_scaled,
12061218
primal_residual_in_scaled_lo,
12071219
primal_residual_in_scaled_up,
12081220
dual_residual_scaled,
@@ -1211,20 +1223,24 @@ qp_solve(Results<T>& results,
12111223
dual_feasibility_rhs_0,
12121224
dual_feasibility_rhs_1,
12131225
dual_feasibility_rhs_3,
1226+
rhs_duality_gap,
1227+
sqrt_max_dim,
12141228
precond,
12151229
data,
12161230
qp_scaled.as_const(),
1217-
detail::vec(x_e),
1218-
detail::vec(y_e),
1219-
detail::vec(z_e),
1231+
detail::vec_mut(x_e),
1232+
detail::vec_mut(y_e),
1233+
detail::vec_mut(z_e),
12201234
stack));
12211235

12221236
if (is_primal_feasible(primal_feasibility_lhs_new) &&
12231237
is_dual_feasible(dual_feasibility_lhs_new)) {
1224-
results.info.pri_res = primal_feasibility_lhs_new;
1225-
results.info.dua_res = dual_feasibility_lhs_new;
1226-
results.info.status = QPSolverOutput::PROXQP_SOLVED;
1227-
break;
1238+
if (results.info.duality_gap<=settings.eps_abs*rhs_duality_gap){
1239+
results.info.pri_res = primal_feasibility_lhs_new;
1240+
results.info.dua_res = dual_feasibility_lhs_new;
1241+
results.info.status = QPSolverOutput::PROXQP_SOLVED;
1242+
break;
1243+
}
12281244
}
12291245

12301246
// vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
@@ -1259,7 +1275,8 @@ qp_solve(Results<T>& results,
12591275
VEG_BIND(
12601276
auto,
12611277
(_, dual_feasibility_lhs_new_2),
1262-
detail::unscaled_primal_dual_residual(primal_residual_eq_scaled,
1278+
detail::unscaled_primal_dual_residual(results,
1279+
primal_residual_eq_scaled,
12631280
primal_residual_in_scaled_lo,
12641281
primal_residual_in_scaled_up,
12651282
dual_residual_scaled,
@@ -1268,12 +1285,14 @@ qp_solve(Results<T>& results,
12681285
dual_feasibility_rhs_0,
12691286
dual_feasibility_rhs_1,
12701287
dual_feasibility_rhs_3,
1288+
rhs_duality_gap,
1289+
sqrt_max_dim,
12711290
precond,
12721291
data,
12731292
qp_scaled.as_const(),
1274-
detail::vec(x_e),
1275-
detail::vec(y_e),
1276-
detail::vec(z_e),
1293+
detail::vec_mut(x_e),
1294+
detail::vec_mut(y_e),
1295+
detail::vec_mut(z_e),
12771296
stack));
12781297
proxsuite::linalg::veg::unused(_);
12791298

0 commit comments

Comments
 (0)