How to hint so that polish has higher success chance without reducing eps?

[colinfang]

I have many small problems of various sizes (number of variables & constraints both between (5 ~ 200)
The quadratic term in the objective is typically quite small.

I want to find a good balance betwen number of iterations and quality of solutions. (I have to use OSQP throughout instead of switching to other optimizers)

The following is one example of the data, for simplification I removed the quadratic term so it is basically a LP

P = scipy.sparse.csc_matrix((9, 9))

A = scipy.sparse.csc_matrix(np.array([
       [ 1.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  1.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  1.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  1.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  1.      ,  0.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  1.      ,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  1.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  1.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  1.      ],
       [-0.34254 ,  0.545208,  0.213992,  0.133316,  0.557486,  0.153901,  0.216159,  0.155356, -0.34134 ],
       [ 0.11394 , -0.255139,  0.446662,  0.371519, -0.260884,  0.386901,  0.449502,  0.389477,  0.113541]
]))

q = np.array([-0.01735301,  0.01893665,  0.02478311,  0.00815631, -0.00315658,  0.02802239,  0.02184102,  0.02329007,  -0.01379083])
lb = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
ub = np.array([ 50, 50, 50, 47.996, 9.8307, 50, 119, 50, 50, 0, 0])

m = osqp.OSQP()
m.setup(P=P, q=q, A=A, l=lb, u=ub, polish=True, eps_abs=1e-4)

results = m.solve()

results.x

# array([4.99868160e+01, 5.00154694e+01, 4.65649348e-02, 1.04621379e+01, 9.84657281e+00, 3.83493845e-02, 4.69118973e-02,
#       3.86318732e-02, 4.99867442e+01])

-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 9, constraints m = 11
          nnz(P) + nnz(A) = 27
settings: linear system solver = qdldl,
          eps_abs = 1.0e-04, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (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: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -3.2322e-01   4.63e+00   2.21e-02   1.00e-01   8.87e-05s
 100  -5.5125e-01   4.69e-02   4.88e-05   1.07e-02   1.58e-04s

status:               solved
solution polish:      unsuccessful
number of iterations: 100
optimal objective:    -0.5513
run time:             2.15e-04s
optimal rho estimate: 7.88e-03

As you can see under this settings, the solution has low accuracy. In fact the solutions sometimes live outside of the bounds.

I know that if I reduce eps, I have a higher chance for a successful polishing

m = osqp.OSQP()
m.setup(P=P, q=q, A=A, l=lb, u=ub, polish=True, eps_abs=1e-4, eps_rel=1e-4)
results = m.solve()

# array([5.00000000e+01, 5.00000000e+01, 1.78439505e-23, 1.06547467e+01, 9.83070000e+00, 1.05272596e-22, 8.61113317e-23,
#        1.32792822e-22, 4.99042768e+01])

# iter   objective    pri res    dua res    rho        time
#    1  -3.2322e-01   4.63e+00   2.21e-02   1.00e-01   1.68e-04s
#  200  -5.5366e-01   2.11e-02   2.76e-07   1.07e-02   2.71e-04s
#  250  -5.5318e-01   7.25e-04   1.24e-05   7.03e-02   3.24e-04s
# plsh  -5.5317e-01   7.99e-16   1.06e-16   --------   3.74e-04s
# 
# status:               solved
# solution polish:      successful
# number of iterations: 250
# optimal objective:    -0.5532
# run time:             3.74e-04s
# optimal rho estimate: 1.27e-02

My problems all have a structure that all variables bounded between 0 to ub. (that's the identity block in the constraints) The number of other constraints are small, perhaps just 1 or 2. and they are equalities to 0. I wonder if there is any smart way to exploit this structure?

I observe that in many cases, when the polishing is unsucessful, the solutions are quite close to optimal. E.g. if a variable ends with 1e-2, then the exact solution is very likely to be 0.
If a variable ends with 49.99 and upper bound is 50, then the exact solution is very likely to be 50.

Therefore, I would like to somehow manually guess a better solution based on the unsucessful polished one, and warm start in osqp so that it has higher chance of polishing.

lb0 = np.array([50, 50, 0,   0,       0,      0, 0, 0,  0,   0, 0])
ub0 = np.array([50, 50, 0,  47.996 ,   9.8307,  0., 0, 0, 50 , 0 0])
x0 = results.x.clip(min=lb0[:-2], max=ub0[:-2])

m.warm_start(x=x0)
# Not sure if I should also update the bounds to hint it harder
m.update(l=lb0, u=ub0)

# iter   objective    pri res    dua res    rho        time
#    1  -5.5285e-01   1.16e-01   2.88e-05   1.07e-02   5.11e-05s
#   25  -5.5422e-01   4.17e-02   8.11e-10   1.07e-02   1.06e-04s
# 
# status:               solved
# solution polish:      unsuccessful
# number of iterations: 25
# optimal objective:    -0.5542
# run time:             1.58e-04s
# optimal rho estimate: 1.82e+00

Unfortunately, it doesn't quite work.

Any suggestions how I can achieve high quality solutions in a small number of iterations here?

(I also tried hint with y, but never really understood how to guess a sensible value for dual)

[colinfang]

Another idea is to change basis of my problem.

In the above case my constraints excluding bounds is basically C^T x = 0

[-0.34254 ,  0.545208,  0.213992,  0.133316,  0.557486,  0.153901,  0.216159,  0.155356, -0.34134 ],
[ 0.11394 , -0.255139,  0.446662,  0.371519, -0.260884,  0.386901,  0.449502,  0.389477,  0.113541]

So I can find the basis of my feasible set of x (Suppose I can always do it efficiently and cached among many similar problems)

# 9 x 7
basis = np.array([
       [ 0.43284408,  0.32325158,  0.30895512,  0.34903597,  0.43625998,  0.35171408, -0.21571652],
       [ 0.24103442,  0.22610833, -0.53221465,  0.22682918,  0.24210448,  0.22809261,  0.28824535],
       [ 0.76326669, -0.18799347,  0.00210348, -0.19878278, -0.23839942, -0.20019214,  0.01877965],
       [-0.18772955,  0.84996649,  0.01625166, -0.15830703, -0.18903354, -0.15941993,  0.00643587],
       [-0.00192829,  0.01305232,  0.77723529,  0.00865042, -0.00220509,  0.00856684,  0.12933445],
       [-0.19859693, -0.15838106,  0.0120342 ,  0.83276659, -0.1999825 , -0.16841241,  0.00979947],
       [-0.23840418, -0.1893031 ,  0.00185509, -0.20017366,  0.75991764, -0.20159305,  0.01906483],
       [-0.2000076 , -0.15949657,  0.01197464, -0.16841461, -0.20140318,  0.83039797,  0.00995318],
       [ 0.0211175 ,  0.00831333,  0.12899384,  0.01177728,  0.02141875,  0.01194481,  0.92341807]
])

Hence I can re-parameter the problem in terms of the coefficients of the basis.

P = scipy.sparse.csc_matrix((7, 7))
A = scipy.sparse.csc_matrix(basis)
# original_q.dot(basis)
q = np.array([-0.0012773 , -0.01149745, -0.01883068,  0.0068729 , -0.0044216 ,  0.00198029, -0.0025005 ])
# Same as original
lb = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0.])
ub = np.array([ 50, 50, 50, 47.996, 9.8307, 50, 119, 50, 50])

m = osqp.OSQP()
m.setup(P=P, q=q, A=A, l=lb, u=ub, polish=True, eps_abs=1e-4, eps_rel=1e-4)

results = m.solver()

-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 7, constraints m = 9
          nnz(P) + nnz(A) = 63
settings: linear system solver = qdldl,
          eps_abs = 1.0e-04, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (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: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -2.8360e-01   5.04e+00   1.23e-02   1.00e-01   8.67e-05s
  75  -5.5176e-01   3.88e-02   5.31e-05   1.42e-02   2.23e-04s

status:               solved
solution polish:      unsuccessful
number of iterations: 75
optimal objective:    -0.5518
run time:             2.89e-04s
optimal rho estimate: 5.79e-03

This problem has 2 less variables and 2 less constraints, but the constraints is no longer sparse.

What do you think of this? Would it typically speed up the convergence or not?