Problems with edge case (Maros Meszaros)

[migarstka]

I came across the following problem when I tried to debug a Maros Meszaros benchmark. The following problem is HS51 from the problem set. It has the following form:

1/2 x' P x + q' x + r
s.t. l <= Ax <= u

We can solve the problem in COSMO with a Box-constraint:

using COSMO


P =[2.0 -2.0  0.0  0.0  0.0;
  -2.0   4.0  2.0  0.0  0.0;
   0.0   2.0  2.0  0.0  0.0;
   0.0   0.0  0.0  2.0  0.0;
   0.0   0.0  0.0  0.0  2.0]

q = [0.; -4; -4; -2; -2];

A = [ 1.0  3.0  0.0  0.0   0.0;
      0.0  0.0  1.0  1.0  -2.0;
      0.0  1.0  0.0  0.0  -1.0;
      1.0  0.0  0.0  0.0   0.0;
      0.0  1.0  0.0  0.0   0.0;
      0.0  0.0  1.0  0.0   0.0;
      0.0  0.0  0.0  1.0   0.0;
      0.0  0.0  0.0  0.0   1.0]

b = zeros(8)

l = [ 4.0;
  0.0;
  0.0;
 -1.0e20;
 -1.0e20;
 -1.0e20;
 -1.0e20;
 -1.0e20]

u = [4.0;
 0.0;
 0.0;
 1.0e20;
 1.0e20;
 1.0e20;
 1.0e20;
 1.0e20;]


 # solve with Box constraint l <= Ax <= u
model = COSMO.Model()
constraint = COSMO.Constraint(A, b, COSMO.Box(l, u))
assemble!(model, P, q, constraint)
result = COSMO.optimize!(model);

This yields the correct result obj_val = -6.0.

Now, lets assume we want to solve the same problem using one-sided constraints:
l <= Ax <= u <==> -Ax + u >= 0 and Ax - l >= 0

model = COSMO.Model()
constraint1 = COSMO.Constraint(-A, u, COSMO.Nonnegatives)
constraint2 = COSMO.Constraint(A, -l, COSMO.Nonnegatives)
assemble!(model, P, q, [constraint1, constraint2])
result = COSMO.optimize!(model);

The solver output:

------------------------------------------------------------------
             COSMO - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2018
------------------------------------------------------------------

Problem:  x ∈ R^{5},
          constraints: A ∈ R^{16x5} (24 nnz), b ∈ R^{16},
          matrix size to factor: 21x21 (441 elem, 73 nnz)
Sets:     Nonnegatives{Float64} of dim: 16
Settings: ϵ_abs = 1.0e-04, ϵ_rel = 1.0e-04,
          ϵ_prim_inf = 1.0e-06, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1.0e-6, α = 1.6,
          max_iter = 2500,
          scaling iter = 10 (on),
          check termination every 1 iter,
          check infeasibility every 40 iter
Setup Time: 0.46ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-8.9378e-01	1.1155e+20	4.3746e+00	1.0000e-01
2	-3.6748e-01	6.6932e+19	5.1392e+00	1.0000e-01
3	-1.5205e+00	4.0159e+19	5.0813e+00	1.0000e-01
4	-1.1850e+00	2.4095e+19	5.1194e+00	1.0000e-01
5	-1.8224e+00	1.4457e+19	5.0379e+00	1.0000e-01
6	-1.7447e+00	8.6743e+18	5.0732e+00	1.0000e-01
7	-2.1014e+00	5.2046e+18	5.0594e+00	1.0000e-01
8	-2.1504e+00	3.1228e+18	5.0766e+00	1.0000e-01
9	-2.3635e+00	1.8737e+18	5.0702e+00	1.0000e-01
10	-2.4491e+00	1.1242e+18	5.0760e+00	1.0000e-01
11	-2.5902e+00	6.7452e+17	5.0741e+00	1.0000e-01
12	-2.6777e+00	4.0471e+17	5.0770e+00	1.0000e-01
13	-2.7802e+00	2.4283e+17	5.0769e+00	1.0000e-01
14	-2.8583e+00	1.4570e+17	5.0785e+00	1.0000e-01
15	8.7145e-01	8.7417e+16	4.3610e+00	1.0000e-01
16	-3.3878e-01	5.2450e+16	8.3304e+00	1.0000e-01
17	-5.2683e+00	3.1470e+16	3.6624e+00	1.0000e-01
18	-5.0070e+00	1.8882e+16	1.6080e+00	1.0000e-01
19	-5.4227e+00	1.1329e+16	3.2440e+00	1.0000e-01
20	-5.7013e+00	6.7976e+15	1.3145e+00	1.0000e-01
21	-5.5454e+00	4.0785e+15	1.6094e+00	1.0000e-01
22	-5.7007e+00	2.4471e+15	1.4807e+00	1.0000e-01
23	-5.5601e+00	1.4683e+15	1.7154e+00	1.0000e-01
24	-5.6145e+00	8.8097e+14	1.5812e+00	1.0000e-01
25	-5.5645e+00	5.2858e+14	1.6702e+00	1.0000e-01
26	-5.5581e+00	3.1715e+14	1.5066e+00	1.0000e-01
27	-5.5191e+00	1.9029e+14	1.5663e+00	1.0000e-01
28	-5.5048e+00	1.1417e+14	1.5274e+00	1.0000e-01
29	-5.4773e+00	6.8504e+13	1.5490e+00	1.0000e-01
30	-5.4648e+00	4.1102e+13	1.5495e+00	1.0000e-01
31	-5.4466e+00	2.4661e+13	1.5535e+00	1.0000e-01
32	-5.4342e+00	1.4797e+13	1.5506e+00	1.0000e-01
33	-5.4222e+00	8.8781e+12	1.5484e+00	1.0000e-01
34	-5.4121e+00	5.3269e+12	1.5501e+00	1.0000e-01
35	-5.4033e+00	3.1961e+12	1.5493e+00	1.0000e-01
36	-5.3957e+00	1.9177e+12	1.5496e+00	1.0000e-01
37	-5.3890e+00	1.1506e+12	1.5492e+00	1.0000e-01
38	-5.3834e+00	6.9036e+11	1.5490e+00	1.0000e-01
39	-5.3783e+00	4.1422e+11	1.5489e+00	1.0000e-01
40	-5.3741e+00	2.4853e+11	1.5489e+00	1.0000e-01
41	-5.7532e+00	1.4912e+11	6.4953e-01	7.3505e-06
42	-5.8355e+00	8.9471e+10	3.8156e-01	7.3505e-06
43	-5.9076e+00	5.3683e+10	2.2797e-01	7.3505e-06
44	-5.9087e+00	3.2210e+10	1.3640e-01	7.3505e-06
45	-5.9238e+00	1.9326e+10	8.1474e-02	7.3505e-06
46	-5.9205e+00	1.1595e+10	4.8892e-02	7.3505e-06
47	-5.9246e+00	6.9573e+09	2.9152e-02	7.3505e-06
48	-5.9229e+00	4.1744e+09	1.7588e-02	7.3505e-06
49	-5.9242e+00	2.5046e+09	1.0414e-02	7.3505e-06
50	-5.9236e+00	1.5027e+09	6.3671e-03	7.3505e-06
51	-5.9240e+00	9.0164e+08	3.6916e-03	7.3505e-06
52	-5.9238e+00	5.4101e+08	2.3388e-03	7.3505e-06
53	-5.9240e+00	3.2458e+08	1.2771e-03	7.3505e-06
54	-5.9239e+00	1.9475e+08	8.9129e-04	7.3505e-06
55	-5.9240e+00	1.1685e+08	4.0917e-04	7.3505e-06

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 55
Optimal objective: -5.924
Runtime: 0.026s (26.27ms)

which yields the wrong result. It seems like the large values in l and u which internally get set as b mess up (or at least dilute) the residual computations.
We can solve the same problem by removing the unnecessary constraint rows 4-8:

l_short = [ 4.0; 0.0; 0.0]

u_short = [4.0; 0.0; 0.0]

A_short = A[1:3, :]

model = COSMO.Model()
constraint1 = COSMO.Constraint(-A_short, u_short, COSMO.Nonnegatives)
constraint2 = COSMO.Constraint(A_short, -l_short, COSMO.Nonnegatives)
assemble!(model, P, q, [constraint1, constraint2], settings = COSMO.Settings(verbose = true, check_termination = 1))
result = COSMO.optimize!(model);

unsurprisingly, this again yields the correct result of obj_val = -6.0