## 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 N2 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 N2 + 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```