Convergence failed in LP when using primalstart

[Ben]

Hi,

I would like to thank the developers for working on cvxopt, a project I use
frequently.

Let me describe a problem I'm trying to solve. I have
c. 1800 LPs, with the same constraints, but different objective functions. Each
has a small number of variables. The objective functions are quite similar.
That is, if we denote the objective vectors by c_1, ..., c_{1800}, then
||c_{i+1} - c_i|| is small for each i.

Since each problem is similar to the one that proceeds it, I thought it should
substantially reduce the total computation time to reuse the primal variables from problem i as a starting point for
problem i + 1. It seems, though, that specifying the primalstart argument can
cause cvxopt to fail to find an optimal solution for even very small problems.

1. Why is this the case?
   
2. Is there a small change I could make to the primalstart values that would ensure convergence [while retaining their feasibility]?
3. I gather that the conelp solver uses an interior point method. Would it be better to use the simplex method for this problem?

A minimal example is appended.

Many thanks,
Ben

def test_problem0(use_primalstart):
    G = np.array([[  15.9893514126975997,   15.9893514126975997,   -1.                ],
                         [  15.9893514126975997,   94.5471809862181942,   -1.                ],
                         [  94.5471809862181942,   15.9893514126975997,   -1.                ],
                         [  94.5471809862181942,   94.5471809862181942,   -1.                ],
                         [ 118.4264850360312948,  118.4264850360312948,   -1.                ],
                         [  11.7649639986445642,   11.7649639986445642,   -1.                ],
                         [  11.7649639986445642,  118.4264850360312948,   -1.                ],
                         [ 118.4264850360312948,   11.7649639986445642,   -1.                ]])
    f = np.array([ 86.2044757760202032,  85.3971784683824069,  85.1577896088097646,
                          17.4115147061936213,   2.6339183736003102,  97.0445533548508195,
                          97.0445533548508195,  97.0445533548508195])
    c0 = np.array([-16.1753152677003733,  -17.1247066863988913, 1.])
    sol0 = lp(matrix(c0), matrix(G), matrix(f))
    c1 = np.array([-18.3559242032701846, -18.9669653736389954,   1.                ])
    ps = {"x": sol0["x"], "s": sol0["s"]} if use_primalstart else None
    sol1 = lp(matrix(c1), matrix(G), matrix(f), primalstart=ps)
    return sol1


In [22]: cds.test_problem0(False)
     pcost       dcost       gap    pres   dres   k/t
 0: -1.0809e+02 -7.6562e+02  3e+02  3e-01  3e+01  1e+00
 1: -8.6188e+01 -1.1545e+02  1e+01  1e-02  1e+00  4e-02
 2: -8.5698e+01 -8.7438e+01  7e-01  8e-04  9e-02  4e-03
 3: -8.6030e+01 -8.6377e+01  1e-01  2e-04  2e-02  5e-03
 4: -8.6028e+01 -8.6063e+01  1e-02  2e-05  2e-03  1e-03
 5: -8.6032e+01 -8.6032e+01  2e-04  3e-07  3e-05  2e-05
 6: -8.6032e+01 -8.6032e+01  2e-06  3e-09  3e-07  2e-07
 7: -8.6032e+01 -8.6032e+01  2e-08  3e-11  3e-09  2e-09
Optimal solution found.
     pcost       dcost       gap    pres   dres   k/t
 0: -1.0637e+02 -7.5648e+02  3e+02  3e-01  3e+01  1e+00
 1: -8.4358e+01 -1.1187e+02  1e+01  1e-02  1e+00  3e-02
 2: -8.3893e+01 -8.5159e+01  5e-01  6e-04  6e-02  4e-03
 3: -8.4121e+01 -8.4375e+01  9e-02  1e-04  1e-02  4e-03
 4: -8.4123e+01 -8.4148e+01  9e-03  1e-05  1e-03  7e-04
 5: -8.4126e+01 -8.4126e+01  9e-05  1e-07  1e-05  8e-06
 6: -8.4126e+01 -8.4126e+01  9e-07  1e-09  1e-07  8e-08
 7: -8.4126e+01 -8.4126e+01  9e-09  1e-11  1e-09  8e-10
Optimal solution found.
Out[22]:
{'dual infeasibility': 1.1555468652274997e-09,
 'dual objective': -84.12579258329689,
 'dual slack': 1.6938597579539248e-11,
 'gap': 8.690176188238481e-09,
 'iterations': 7,
 'primal infeasibility': 1.2456718081131636e-11,
 'primal objective': -84.12579255851365,
 'primal slack': 9.337570810781302e-10,
 'relative gap': 1.0329978385872602e-10,
 'residual as dual infeasibility certificate': None,
 'residual as primal infeasibility certificate': None,
 's': <8x1 matrix, tc='d'>,
 'status': 'optimal',
 'x': <3x1 matrix, tc='d'>,
 'y': <0x1 matrix, tc='d'>,
 'z': <8x1 matrix, tc='d'>}

In [23]: cds.test_problem0(True)
     pcost       dcost       gap    pres   dres   k/t
 0: -1.0809e+02 -7.6562e+02  3e+02  3e-01  3e+01  1e+00
 1: -8.6188e+01 -1.1545e+02  1e+01  1e-02  1e+00  4e-02
 2: -8.5698e+01 -8.7438e+01  7e-01  8e-04  9e-02  4e-03
 3: -8.6030e+01 -8.6377e+01  1e-01  2e-04  2e-02  5e-03
 4: -8.6028e+01 -8.6063e+01  1e-02  2e-05  2e-03  1e-03
 5: -8.6032e+01 -8.6032e+01  2e-04  3e-07  3e-05  2e-05
 6: -8.6032e+01 -8.6032e+01  2e-06  3e-09  3e-07  2e-07
 7: -8.6032e+01 -8.6032e+01  2e-08  3e-11  3e-09  2e-09
Optimal solution found.
     pcost       dcost       gap    pres   dres   k/t
 0: -8.4126e+01 -7.5648e+02  2e+02  3e-11  3e+01  1e+00
 1: -8.4126e+01 -9.0849e+01  2e+00  3e-13  3e-01  1e-02
 2: -8.4126e+01 -8.4193e+01  2e-02  3e-15  3e-03  1e-04
 3: -8.4126e+01 -8.4126e+01  2e-04  1e-16  3e-05  1e-06
 4: -8.4126e+01 -8.4126e+01  2e-06  1e-16  4e-07  1e-08
 5: -8.4126e+01 -8.4286e+01  2e-08  2e-16  2e-03  1e-10
 6: -8.4126e+01 -8.4213e+01  2e-10  3e-16  3e-04  1e-12
 7: -8.4126e+01 -8.4127e+01  2e-12  2e-16  2e-05  1e-14
 8: -8.4126e+01 -8.3329e+01  4e-14  5e-16  1e-02  3e-15
 9: -8.4126e+01 -8.3899e+01  9e-15  2e-16  2e-03  7e-16
10: -8.4126e+01 -8.8961e+01  5e-15  1e-16  3e-02  7e-16
11: -8.4126e+01 -8.8907e+01  4e-15  2e-16  3e-02  7e-16
12: -8.4126e+01 -8.8910e+01  4e-15  2e-16  3e-02  7e-16
13: -8.4126e+01 -8.8910e+01  4e-15  2e-16  3e-02  7e-16
14: -8.4126e+01 -8.8910e+01  4e-15  3e-16  3e-02  7e-16
Terminated (singular KKT matrix).
Out[23]:
{'dual infeasibility': 0.03353162360643785,
 'dual objective': -88.91047464421752,
 'dual slack': 6.3520619921953765e-18,
 'gap': 4.4212678256667316e-15,
 'iterations': 14,
 'primal infeasibility': 2.8100540907245236e-16,
 'primal objective': -84.12579256099686,
 'primal slack': 7.090794118668737e-16,
 'relative gap': 5.255543741190926e-17,
 'residual as dual infeasibility certificate': 0.011886960818525817,
 'residual as primal infeasibility certificate': None,
 's': <8x1 matrix, tc='d'>,
 'status': 'unknown',
 'x': <3x1 matrix, tc='d'>,
 'y': <0x1 matrix, tc='d'>,
 'z': <8x1 matrix, tc='d'>}

[Martin]

Hi Ben,

1) Warm-starting an interior-point method can be difficult. Using a point close to or on the boundary of the feasible set as a starting point can lead to numerical problems.

2) See this paper by Skajaa et al.: Warmstarting the Homogeneous and Self-Dual Interior Point Method for Linear and Conic Quadratic Problems

3) Yes, CVXOPT uses a path-following interior-point method. It may be easier to get started with warm-starting if you use a simplex method. You can use the GLPK simplex solver in CVXOPT by setting the solver option (if you installed CVXOPT's GLPK module): http://cvxopt.org/userguide/coneprog.html?highlight=simplex#optional-solvers

Kind regards, 

Martin