Nonconvex QP

This is an example of how to find the global optimum to a non-convex quadratic programming (QP). The problem set-up is done using tomSym, and a mixed-integer solver, such as CPLEX or GUROBI does the numeric work.

We want to minimize:

  0.5*x'*Q*x'

Subject to:

  Aeq*x == beq
  0 <= x <= 1

If Q were positive semidefinite, this would be a Quadratic Programming. However, Q is indefinite, making the problem much harder. We could easily find a local optimum using a nonlinear solver, but finding a global optimum is hard.

To do this, we set up the first-order Karush-Kuhn-Tucker (KKT) conditions, and obtain a mixed integer programming. This will yield a solution that is proven to be within tolerances of the global optimum.

% Input parameters.
% (Q is a randomly generated 20x20 symmetric matrix.)
Q = setSymmetric([0.99813;-0.9328;-0.41602;0.75425;0.36975;-0.16011;...
    0.20505;-0.22754;0.21575;-0.63597;-0.29342;0.27946;-0.91458;...
    -0.33758;0.72371;0.91155;-0.95604;-0.58809;0.11226;-0.60359;...
    1.3152;-0.45212;-0.0081002;0.73651;-0.89605;1.0736;1.0053;...
    -0.089598;-0.86427;-0.83047;-0.35381;-0.10368;-0.20189;-0.5286;...
    0.74354;1.142;0.11788;0.57156;0.16935;-1.854;0.80398;1.0658;...
    -1.0921;0.4812;0.17866;0.26736;0.68732;-0.04921;0.51998;...
    -0.4118;0.20122;0.62257;-0.48735;-1.2293;-0.27855;0.020442;...
    -1.1031;-0.85323;0.084205;1.3643;-0.095637;0.12529;0.071885;...
    0.74562;-0.36408;-1.0812;0.32137;0.79026;0.4365;1.536;0.73981;...
    1.5732;1.0358;0.86736;0.70902;-1.5217;-0.93981;-0.94568;0.77429;...
    -0.77916;1.1569;0.26974;-0.00055197;-0.13087;0.25253;-0.83893;...
    0.090814;0.025564;0.3809;-0.079573;0.041911;0.12882;1.0671;...
    -0.42695;-1.1162;-0.00038358;0.37305;-0.047864;-0.83988;...
    -0.074761;0.11005;0.33679;-0.91941;0.13811;1.2319;-0.077739;...
    0.40736;-0.70684;0.24293;-0.64002;1.631;-0.36573;0.34084;...
    -0.046399;-0.85519;0.065744;-0.99468;-0.47995;0.91457;-0.30563;...
    1.2142;-0.9352;-1.1335;0.32376;-0.68342;0.92591;-0.13294;-0.77732;...
    0.41786;-0.38176;1.4631;-0.452;-0.54978;0.91733;-0.17072;0.79189;...
    0.22303;0.64006;-0.22713;-0.78101;0.17289;1.4549;-1.1471;-2.4195;...
    -0.0197;0.11204;0.41687;0.079953;-0.79564;0.55421;0.084912;...
    -0.051431;-0.25787;0.29086;0.14115;-1.1786;-0.62126;-0.97309;...
    0.33165;-0.59094;-0.1878;-0.20598;-0.48744;0.16878;-1.0183;...
    -0.12057;-0.18849;-0.16661;0.20649;-0.53152;0.81186;0.86475;...
    0.43589;-0.68085;0.59335;0.49119;-0.75623;0.21621;-1.1849;...
    0.29998;-0.5439;0.3483;-0.8631;0.30066;0.092835;0.97066;0.085208;...
    -1.0018;-0.36108;-1.0288;0.28215;0.23805;-0.016555;-1.1906;...
    -0.23784;-0.19766;1.0435;0.39938;-0.19682;0.33933;-0.67583;...
    -0.34775;-0.1395;0.17759;0.082428;-0.43612;1.6515;-1.2706;...
    -0.99358;-1.4966]);
Aeq = [0 1 1 1 1 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0; ...
    0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0];
beq = [5;2];

N = size(Q,1);

% Lower and upper bounds on x, expressed as A*x <= b
% We could also put arbitrary linear inequalities here.
A = [eye(N); -eye(N)];
b = [ones(N,1); zeros(N,1)];

M = size(A,1);

% Upper bound on s and y
% If these values are too small, then no solution will be found.
% Too large numbers may make the solver less efficient, or lead to
% incorrect solutions due to rounding.
bnds = 1; 
bndy = 100;

% Symbolic variables
x = tom('x',N,1);
s = tom('s',M,1); % Slack variables
y = tom('y',M,1); % Lagrange multipliers for inequality constraints
r = tom('r',M,1,'integer'); % 0 for inactive constraints, 1 for active.
yeq = tom('yeq',size(Aeq,1)); % Multipliers for equality constraints

kktconditions = {
    A'*y+Aeq'*yeq+Q*x==0
    Aeq*x==beq
    A*x+s == b
    0 <= y
    0 <= s
    y <= bndy*r
    s <= bnds*(1-r)
    };
    
% This linear objective will be equal to 0.5*x'*C*x at optimum
obj = -0.5*(y'*b+yeq'*beq);

%% Solve the problem
options = struct;
options.Prob.MIP.cpxControl.EPGAP = 1e-4; % Tolerance relative to global optimum.
solution = ezsolve(obj,kktconditions,[],options);