Verify Optimality by Solving QUBO as MILP

Note

Installation Required: This functionality requires MATLAB Support Package for Quantum Computing.

A Quadratic Unconstrained Binary Optimization (QUBO) problem can be difficult to solve exactly. The solve function might return solutions that are not globally optimal. However, for small enough QUBO problems, you can convert the QUBO problem to a mixed-integer linear programming (MILP) problem. To solve an MILP problem, use intlinprog (Optimization Toolbox), which requires an Optimization Toolbox™ license. If intlinprog returns a solution with gap 0, the solution is guaranteed to be optimal. This approach is practical for up to 100 or 200 variables. For larger problems, the N² problem size slows the solution too much.

Theory: Convert Quadratic Problem to Linear Integer Program

Suppose that your binary variable vector is a column vector x with N elements. Let xij be an N-by-N binary matrix with entries xij(i,j) = x(i)*x(j). You can write the objective function

x'*Q*x + c'*x + d = sum(Q.*xij,"all") + c'*x + d.

This formulation converts the problem from quadratic in x to linear in xij and x, which is a simplification. However, the number of variables increases from N to N² + N, which is a complication.

The following three linear inequalities tie the matrix xij to the vector x and ensure that, at the solution, xij(i,j) = x(i)*x(j). For all i and j,

xij(i,j) >= x(i) + x(j) - 1

xij(i,j) <= x(i)

xij(i,j) <= x(j).

Represent these inequality constraints in a structured way. Define xij as an N-by-N matrix of optimization variables, and define x as a column vector of optimization variables with N entries. Use the following code to create the constraints for an optimization problem named prob.

prob.Constraints.f = xij >= repmat(x,1,N) + repmat(x',N,1) - 1;
prob.Constraints.g = xij <= repmat(x,1,N);
prob.Constraints.h = xij <= repmat(x',N,1);

Example: Solve QUBO Using MILP

Take a matrix Q of size 100-by-100 generated by the makeq Helper Function given at the end of this example.

Q = makeq(100,1);

Create a QUBO problem from Q and solve the problem using the default tabu search algorithm.

qprob = qubo(Q);
rng default % For reproducibility
result = solve(qprob)

result = 

  quboResult with properties:

                BestX: [100×1 double]
    BestFunctionValue: -9610
      AlgorithmResult: [1×1 tabuSearchResult]

The returned objective function value is –9610. Check whether this answer is reliable by resolving the problem.

result = solve(qprob)

result = 

  quboResult with properties:

                BestX: [100×1 double]
    BestFunctionValue: -9610
      AlgorithmResult: [1×1 tabuSearchResult]

The returned objective function value does not change. However, the tabu search algorithm does not necessarily give consistent results.

To obtain a reliable result, formulate the problem as an MILP using optimization variables.

N = size(Q,1);
x = optimvar("x",N,Type="integer",LowerBound=0,UpperBound=1);
% Need N^2 variables for MILP
xij = optimvar("xij",N,N,Type="integer",LowerBound=0,UpperBound=1);
prob = optimproblem;
% The next three constraint arrays are the connections between xij and x
prob.Constraints.f = xij >= repmat(x,1,N) + repmat(x',N,1) - 1;
prob.Constraints.g = xij <= repmat(x,1,N);
prob.Constraints.h = xij <= repmat(x',N,1);
% Formulate the objective in terms of xij
% If you have a linear term c and a constant term d, the objective is
% prob.Objective = sum(Q.*xij,"all") + dot(c,x) + d;
prob.Objective = sum(Q.*xij,"all");
% Solve calls intlinprog
[solxij,fv] = solve(prob);

Solving problem using intlinprog.
Running HiGHS 1.5.3: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
30000 rows, 10100 cols, 69900 nonzeros
23742 rows, 8014 cols, 55398 nonzeros
18978 rows, 6426 cols, 44282 nonzeros
18978 rows, 5701 cols, 44282 nonzeros
Objective function is integral with scale 1

Solving MIP model with:
   18978 rows
   5701 cols (5701 binary, 0 integer, 0 implied int., 0 continuous)
   44282 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   -22706          inf                  inf        0      0      0         0     0.6s
 R       0       0         0   0.00%   -11353          1368             929.90%        0      0      0      5043     0.8s
 C       0       0         0   0.00%   -11353          -36                Large        8      5      0      5049     1.3s
         0       0         0   0.00%   -10745          -36                Large      295     98      0      7419     6.3s
 L       0       0         0   0.00%   -10745          -878            1123.80%      340    104      0      7425    13.8s

0.5% inactive integer columns, restarting
Model after restart has 16713 rows, 5669 cols (5669 bin., 0 int., 0 impl., 0 cont.), and 38997 nonzeros

         0       0         0   0.00%   -10745          -878            1123.80%       86      0      0     13993    13.9s
         0       0         0   0.00%   -10739          -878            1123.12%       86     80      9     15524    14.4s
 L       0       0         0   0.00%   -10657          -1134            839.77%      192    103      9     16083    16.9s

Symmetry detection completed in 0.4s
Found 619 full orbitope(s) acting on 1238 columns

 L       0       0         0   0.00%   -10657          -1216            776.40%      192    102      9     17007    18.5s
 T       0       0         0   0.00%   -10657          -8078             31.93%      192    102      9     21224    19.4s
 T      21       0         5   0.02%   -10657          -8158             30.63%      197    102    154     26697    19.6s
 T      25       0         6   0.03%   -10657          -8182             30.25%      198    102    185     26925    19.7s
 T      30       0         9   0.03%   -10657          -8268             28.89%      201    102    256     27082    19.7s
 T      38       0        14   0.07%   -10657          -8362             27.45%      206    102    342     28022    20.0s
 T      64       0        26   0.21%   -10657          -8370             27.32%      218    102    590     31876    20.8s
 T      69       1        29   0.23%   -10657          -8408             26.75%      221    102    652     32418    21.0s
 T      82       1        36   0.33%   -10657          -8472             25.79%      228    102    728     34629    21.7s
 T      84       1        37   0.35%   -10657          -8594             24.01%      229    102    749     35443    21.8s
 T      86       1        38   0.37%   -10657          -8616             23.69%      230    102    750     36173    21.9s
 L     102       9        44   0.44%   -10657          -8766             21.57%      239    105    970     38098    23.4s

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

 T     129      10        55   1.61%   -10657          -8878             20.04%      250    105   1148     49252    25.8s
 T     141       9        59   1.75%   -10657          -8926             19.39%      254    105   1180     52528    26.1s
 T     164      11        71   2.81%   -10657          -9014             18.23%      266    105   1374     59822    27.6s
 L     202      12        88   4.79%   -10508          -9224             13.92%      313    109   1698     71134    31.3s
 T     212      10        91   6.54%   -10508          -9282             13.21%      316    109   1792     77167    32.2s
 T     237       9       104   9.96%   -10508          -9298             13.01%      329    109   1936     82816    33.5s
 T     239       9       105  10.06%   -10508          -9328             12.65%      330    109   1942     83049    33.6s
 T     286       8       129  15.43%   -10508          -9404             11.74%      354    109   2335     96674    36.7s
 T     289       8       130  15.48%   -10508          -9408             11.69%      355    109   2339     97832    36.7s
 T     298       8       135  16.80%   -10508          -9426             11.48%      360    109   2408    100588    37.4s
 T     299       8       136  17.19%   -10508          -9610              9.34%      361    109   2416    100607    37.7s
       332       7       154  61.72%   -10333          -9610              7.52%      356    136   2826    117870    42.8s
       340       4       159  85.16%   -10333          -9610              7.52%     1794    148   2858    125406    47.9s

Solving report
  Status            Optimal
  Primal bound      -9610
  Dual bound        -9610
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    -9610 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            50.64 (total)
                    0.65 (presolve)
                    0.00 (postsolve)
  Nodes             351
  LP iterations     133983 (total)
                    0 (strong br.)
                    6663 (separation)
                    14353 (heuristics)

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06.
The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

The returned gap is 0, meaning the solution is guaranteed to be optimal. The returned objective function value is –9610, which is the same as the tabu search result. You can see this value at the end of the iterative display in the Solving report, in the Solution status section.

Check that the returned variable solxij.x gives the same objective function value in the quadratic expression x'*Q*x.

myx = round(solxij.x); % Round any fractional entries
% Not needed here because the reported int. viol. = 0
disp(myx'*Q*myx)

       -9610

Helper Function

This code creates the makeq helper function.

function Q = makeq(N,seed)
% N must be a positive integer, seed is a random stream seed
% Q is an N-by-N sparse symmetric integer matrix, values -100 through 100
% Q has about N^2/10 nonzeros
% Q is modeled after Beasley 1998

% Copyright 2022 The MathWorks, Inc.
if nargin == 1
    seed = 0;
end
strm = RandStream("mt19937ar",Seed=seed);
Q = randi(strm,201,N);
Q = Q - 101; % Random integers uniform from -100 through 100
JJ = rand(strm,N) >= 0.1;
Q(JJ) = 0; % 10% density of nonzeros
for i = 1:N
    for j = 1:i
        Q(i,j) = 0; % Make upper triangular, 0 on diagonal
    end
end
Q = sparse(Q + Q'); % Symmetrize
end