1

I am trying to solve a simple quadratic program using CVXOPT and am troubled by the fact that I can guess a feasible solution better than the optimum provided by the solver. The optimisation is of the form:

enter image description here

I will provide the definitions of P,q,G,h,A, and b at the end. When I import and run:

from cvxopt import matrix, spmatrix, solvers
# Code that creates matrices goes here
sol = solvers.qp(P, q, G, h, A, b)

the result is:

   pcost       dcost       gap    pres   dres
 0:  0.0000e+00 -5.5000e+00  6e+00  6e-17  4e+00
 1:  0.0000e+00 -5.5000e-02  6e-02  1e-16  4e-02
 2:  0.0000e+00 -5.5000e-04  6e-04  3e-16  4e-04
 3:  0.0000e+00 -5.5000e-06  6e-06  1e-16  4e-06
 4:  0.0000e+00 -5.5000e-08  6e-08  1e-16  4e-08
Optimal solution found.

Objective = 0.0

However I can define a different solution guessed_solution that is feasible and further minimises the objective:

guessed_solution = matrix([0.5,0.5,0.0,0.0,0.0,0.0,0.5,0.5,0.0,0.0,1.0])

# Check Ax = b; want to see zeroes
print(A * guessed_solution - b)
>>>
[ 0.00e+00]
[ 0.00e+00]
[ 2.78e-17]

# Check Gx <= h; want to see non-positive entries
print(G * guessed_solution - h)
>>>
[-5.00e-01]
[-5.00e-01]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[-5.00e-01]
[-5.00e-01]
[-1.00e+00]
[-1.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[-1.00e+00]

# Check objective
print(guessed_solution.T * P * guessed_solution + q.T * guessed_solution)
>>>[-6.67e-01]

This results in an objective that is -2/3, clearly less than 0. I presume that the 2.78e-17 error in the Ax=b test is not relevant.

Any help on resolving this would be appreciated! And below is the definition of relevant matrices in code (biggest matrix is 11 by 11).

P = matrix([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/3.0, 0.0, 2.0/3.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/3.0, 2.0/3.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0],[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0/3.0, 0.0, -2.0/3.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0/3.0, -2.0/3.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]).T
q = matrix([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

A = matrix([[0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],[0.0, 1.0, 1.0/3.0, 2.0/3.0, 0.0, -1.0, -1.0/3.0, -2.0/3.0, 0.0, 0.0, 0.0]]).T
b = matrix([1.0, 1.0, 0.0])

G = spmatrix([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0], [0,1,2,3,4,5,6,7,8,9,10,11,12,13], [0,1,2,3,4,5,6,7,8,9,10,8,9,10])
h = matrix([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0])
N. Mao
  • 499
  • 1
  • 10
  • 19

1 Answers1

1

Your quadratic form is not valid in regards to the assumptions.

It needs to be PSD (and symmetric).

Making it symmetric:

P = (P + P.T) / 2

will lead to cvxopt showing an error, which is due to P being indefinite:

import numpy as np

np_matrix = np.array(P)
print(np.linalg.eigvalsh(np_matrix))

#[-8.16496581e-01 -7.45355992e-01 -5.77350269e-01 -2.40008780e-16 -6.33511351e-17 -4.59089160e-17 -3.94415555e-22  5.54077304e-17  5.77350269e-01  7.45355992e-01  8.16496581e-01]

You got a solver designed for convex-optimization problems (if and only if P is PSD) feeded by some non-convex optimization problem. This won't work (in general).

sascha
  • 32,238
  • 6
  • 68
  • 110