Primal infeasibility detection in EMOSQP

Many thanks for the OSQP developers for sharing this beautiful piece of code !

I have ~1’500 instances of QP’s with identical structure, where I compare the solvers
quadprog (Matlab, interior point)
OSQP (class based)
EMOSQP (code generated)

7 problem instances are known to be primal infeasible.

quadprog and OSQP class based detect properly these 7 primal infeasibility condition.
But EMOSQP (code generated), with identical setup runs into MAX_ITER instead of detecting primal infeasibility !?

Code snippet :

prob.setup(P, q, AA, l, u, ‘eps_abs’, 1e-3, ‘eps_rel’, 1e-3, ‘eps_prim_inf’, 1e-4, ‘eps_dual_inf’, 1e-4, ‘scaling’, 1, ‘rho’, 5e-4, ‘adaptive_rho’, 1, ‘adaptive_rho_interval’, 25, ‘polish’, 0, ‘warm_start’, 1, ‘verbose’, 1, ‘sigma’, 1e-6, ‘delta’, 1e4, ‘max_iter’, 4000);

prob.codegen(‘OSQP’, ‘parameters’, ‘matrices’);

Many thanks in advance for advice,

Raoul

It could be that the non-codegen version is updating rho at different points, since the default behaviour is to update after a set amount of the factorization time. I think the codegen version you will get everything pre-factored, so the rho-update happens instead at fixed intervals (25 iterations as you have it).

You might see this in the output log, since I think we report the current rho value as we go.

Thanks for this remark,
indeed the problem might be related to rho updating at different points.

For the class based (non-codegen) there happen 6 factorizations :


       OSQP v0.6.2  -  Operator Splitting QP Solver
          (c) Bartolomeo Stellato,  Goran Banjac
    University of Oxford  -  Stanford University 2021

problem: variables n = 310, constraints m = 3122
nnz(P) + nnz(A) = 12660
settings: linear system solver = qdldl,
eps_abs = 1.0e-03, eps_rel = 1.0e-03,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 5.00e-04 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: off, time_limit: off

iter objective pri res dua res rho time
1 -8.1989e+05 1.79e+04 9.41e+03 5.00e-04 6.15e-03s
200 -3.0524e+06 9.21e+03 2.56e+00 1.55e+00 2.28e-02s
400 -3.0523e+06 9.21e+03 1.25e+01 5.47e+01 3.59e-02s
425 1.0000e+30 9.21e+03 1.69e+01 5.47e+01 3.74e-02s

status: primal infeasible
number of iterations: 425
run time: 3.75e-02s
optimal rho estimate: 1.31e+02

K>> results.info

ans =

struct with fields:

         iter: 425
       status: 'primal infeasible'
   status_val: -3
status_polish: 0
      obj_val: Inf
      pri_res: 9.2089e+03
      dua_res: 16.8749
   setup_time: 0.0060
   solve_time: 0.0316
  update_time: 6.1885e-04
  polish_time: 0
     run_time: 0.0375
  rho_updates: 6
 rho_estimate: 131.3437

It is clear that the codegen version does many more rho updates.

Is it possible that the rho update scheme becomes unstable ?

I will try to increase adaptive_rho_interval …

No change when increasing adaptive_rho_interval from 25 to 100.
Primal infeasibility is still not detected properly in EMOSQP.

Strange since the solver parameters are supposed to be identical to the class based version of OSQP …

Are you able to post the progress log that OSQP produces an output in each case, i.e. also for the codegen version?

It is almost certainly a settings issue, since the the code that is produced by codegen doesn’t actually touch the solver code, it just pre-allocates memory and populates the initial KKT matrix factorization.

I agree.
Progress log of EMOSQP is not possible for the moment, although verbose was 1 when the codegen method was called.

I will implement “current_settings” and “update_settings” in the EMOSQP wrapper and keep you updated.

Many thanks,

Raoul

I checked the update points and values of work->settings->rho in the codegen version : they are identical to the class based version of osqp.
6 rho updates at the same points to the same values.

The class based osqp stops after 425 iterations with “primal infeasibility”, which is correct.

EMOSQP continues and runs into max_iter (4’000) while further increasing rho from rho=5e1 @iter 425 to rho=6e5 @iter 4’000.

Unclear for me why the behaviours are different. I don’t see any difference in the settings.

Any clue ?

Thanks in advance and best regards,

Raoul

@paul.goulart : of course the osqp C code is identical.
But there are quite a lot of #if EMBEDDED in the code which could potentially cause different behaviours, e.g. this one in auxil.c :

#ifndef EMBEDDED

// Normalize infeasibility certificates if embedded is off
// NB: It requires a division
if ((work->info->status_val == OSQP_PRIMAL_INFEASIBLE) ||
    ((work->info->status_val == OSQP_PRIMAL_INFEASIBLE_INACCURATE))) {
  norm_vec = vec_norm_inf(work->delta_y, work->data->m);
  vec_mult_scalar(work->delta_y, 1. / norm_vec, work->data->m);
}

if ((work->info->status_val == OSQP_DUAL_INFEASIBLE) ||
    ((work->info->status_val == OSQP_DUAL_INFEASIBLE_INACCURATE))) {
  norm_vec = vec_norm_inf(work->delta_x, work->data->n);
  vec_mult_scalar(work->delta_x, 1. / norm_vec, work->data->n);
}

#endif /* ifndef EMBEDDED */