Reducing Binary Quadratic Forms for More Scalable Quantum Annealing
Georg Hahn Hristo N. Djidjev
quant-ph cs.DS




Reducing Binary Quadratic Forms for More Scalable Quantum Annealing
    Georg Hahn
Imperial College, London, UK
Hristo Djidjev (PI)
Los Alamos National Laboratory
    
==============================================================================================







Recent advances in the development of commercial quantum annealers such as the D-Wave 2X allow solving NP-hard optimization problems that can be expressed as quadratic unconstrained binary programs. However, the relatively small number of available qubits (around 1000 for the D-Wave 2X quantum annealer) poses a severe limitation to the range of problems that can be solved. This paper explores the suitability of preprocessing methods  for reducing the sizes of the input programs and thereby the number of qubits required for their solution on quantum computers. Such methods allow us to determine the value of certain variables that hold in either any optimal solution (called strong persistencies) or in at least one optimal solution (weak persistencies). We investigate preprocessing methods for two important NP-hard graph problems, the computation of a maximum clique and a maximum cut in a graph. We show that the identification of strong and weak persistencies for those two optimization problems is very instance-specific, but can lead to substantial reductions in the number of variables.







§ INTRODUCTION





 §.§ Quantum annealing and D-Wave

The recent availability of the first commercial quantum computers of D-Wave Systems Inc. <cit.> provides a novel tool to help finding,  by a process known as quantum annealing, global solutions to NP-hard problems that seem to be  difficult to obtain classically <cit.>.

As of 2017, the most recent version of the company's quantum annealer, D-Wave 2000Q, is capable of handling around 2000 qubits of quantum information. Each qubit consists of a superconducting loop on the D-Wave chip which simultaneously encodes a quantum superposition of -1 and +1 via two superimposed currents in both clockwise and counter-clockwise directions <cit.>. Upon completion of the annealing process, the system turns classical and each qubit takes a classical (binary) value that can be read and assigned to a variable. D-Wave is designed to minimize an  objective function consisting of a sum of linear and quadratic binary contributions, given as

    H=H(x_1,…,x_N)=∑_i ∈ V a_i x_i + ∑_(i,j) ∈ E a_ij x_i x_j,

where a_i, a_ij∈ℝ  and V = {1,…,N}, E=V × V for N ∈ qubits <cit.>. The functional type (<ref>) is called the Hamiltonian encoding the  energy of a quantum system. In (<ref>) qubits are encoded by variables x_i, E denotes the dependency structure between the qubits (the set of interactions), and a_i and a_ij encode physical characteristics of the qubits and the links (called couplers) between them, respectively. During annealing, the quantum system tries to reach a configuration of minimal energy of H. A quadratic form of type (<ref>) with x_i ∈{-1,1} is called an Ising model. By replacing x_i in (<ref>) with (x_i+1)/2, one gets a quadratic form of the same type as (<ref>), but with x_i ∈{0,1}. The resulting problem is called a quadratic unconstrained binary optimization (QUBO) problem. Any QUBO or Ising problem can be solved on D-Wave as long as there are enough qubits and couplers between those qubits available on the D-Wave architecture to represent the system (<ref>).



 §.§ Aim of this work

Although <cit.> show that certain instances designed to fit the D-Wave qubit interconnection network can be solved magnitudes faster with D-Wave than with its classical analogues, the limitation of only around 1000 qubits and about 3000 couplers available on D-Wave 2X (DW) at Los Alamos National Laboratory  poses a severe restriction. For most NP-hard optimization problems, due to the density of the Hamiltonian H, the largest problems able to fit DW are of size about only 45. 

To mitigate the current size limitations of the D-Wave architecture, this work tries to explore the usability of preprocessing methods allowing to determine the value of certain variables in a QUBO or Ising problem. Fixing all such variables then results in a new Hamiltonian of reduced size, consisting only of the unresolved variables. We investigate the usage of partitioning and roof duality methods for two important NP-hard graph problems, the Maximum Clique and the Maximum Cut problem, which have multiple applications including network analysis, bioinformatics, and computational chemistry. Given an undirected graph G = (V, E), a clique is a subset S of the vertices forming a complete subgraph, meaning that any two vertices of S are connected by an edge in G. The clique size is the number of vertices in S, and the Maximum Clique problem is to find a clique in G with a maximum number of vertices <cit.>. A cut S of a graph G is the set of (cut) edges between S and V ∖ S in G. The Maximum Cut problem is to find a maximum cut, which is a cut of maximum size.



 §.§ Previous work

In <cit.>,
optimization methods for quadratic pseudo-boolean functions, which are multi-linear polynomial functions of boolean variables, are considered. The authors analyze a variety of techniques based on roof duality, introduced in <cit.>, in order to determine partial assignments of values or relations between variables in such pseudo-boolean functions. A comprehensive overview of efficient preprocessing techniques for QUBO optimization can be found in <cit.>: Those techniques consists of methods based on first and second order derivatives, roof duality and probing techniques resulting in lower bounds on the minimum energy of a QUBO, the values in either any (strong persistencies) or at least one optimum (weak persistencies), binary relations between variables in some or every optimum, as well as decomposition techniques of the original problem into several smaller pairwise independent QUBO subproblems.

A more efficient implementation of the probing techniques of <cit.>, combined with the aforementioned roof duality, is presented in <cit.>, with a special focus on computer vision applications. The code for the methods investigated in <cit.> are made available for C++ in the Qpbo package of <cit.>. Bindings of this package for Python are available (package PyQpbo, which is part of the PyStruct package of <cit.>) and will be used in the experiments of this article.

The article is organized as follows. Section <ref> gives a detailed overview of the optimization methods used in this report. Moreover, details on the optimization problems for the two NP-hard graph problems we consider, Maximum Clique and Maximum Cut, are given. Section <ref> focuses on the Maximum Clique problem. We state its formulation as a QUBO, highlight the use of a decomposition method based on graph partitioning methods into smaller subproblems, and present experimental results for Qpbo. For this we employ several test graphs, both stand-alone and in connection with the decomposition approach. Section <ref> repeats the QUBO formulation and a Qpbo experimental study for the Maximum Cut problem. We conclude with a discussion in Section <ref>.

In the rest of the paper, we denote a graph as G=(V,E), where V = { 1,…,n } is a set of n vertices and E is a set of undirected edges.






§ BACKGROUND


We briefly overview the techniques used in <cit.> and <cit.> in order to identify strong and weak persistencies. Given a QUBO instance, it is first represented as a (quadratic) posiform 
    ϕ(x) = ∑_T ⊆ L a_T ∏_u ∈ T u,
 where a_T ≥ 0 for all T ≠∅ and L is a set of literals over the n variables {1,…,n} and their complements. It is shown in <cit.> that any posiform (<ref>) can be equivalently expressed as a quadratic form 
    Φ(x) = a_0 + ∑_u ∈ L a_u u + ∑_u,v ∈ L a_uv u v.
 Moreover, <cit.> shows that there is a one-to-one correspondence between posiforms and implication networks, which are capacitated networks G_Φ=(V,E) with n=|V| nodes and non-negative capacities γ_uv on the edges connecting any two nodes u,v ∈ V, u ≠ v, V={1,…,n}.

The implication network G_Φ has an important application in the computation of the roof dual, which is a lower bound on the value of a quadratic posiform <cit.>. The roof dual is computed with the help of a maximum flow through the implication network G_Φ.

More importantly for our work, the flow analysis of the implication network also identifies strong and weak persistencies for problem (<ref>), which are then transformed for persistencies for problem (<ref>).

As shown in <cit.>, the effectiveness of the roof duality method for finding persistencies can be improved by using probing, which is a method that assigns values of 0 and 1 to some variables and then compares any persistencies found for these alternative assignments. In some cases, such probing analysis can (a) determine additional persistencies for the original problem or (b) identify relationships between variables which must hold in a (local) optimum, meaning that some pairs of variables will have equal or inverted values. This then allows to substitute one of the variables in each pair with the other (or its inverse).






§ THE MAXIMUM CLIQUE PROBLEM





 §.§ QUBO formulations


We are interested in finding a maximum clique of a graph G=(V,E), where V={1,…,n}. We assign a decision variable x_i to each vertex i, i ∈{1,…,n}, indicating whether or not vertex i is part of the maximum clique. One way to formulate the Maximum Clique problem as a QUBO is to use the equivalence between Maximum Clique and the maximum independent set problem. The latter asks for a maximum independent set of vertices of G, that is a set such that no two vertices from it are connected by an edge in G. Specifically, an independent set of H = (V, E) defines a clique in graph G=(V,E), where E is the complement of set E. This leads to a constrained formulation of Maximum Clique given by

    x_v ∈{0,1}maximize ∑_v ∈ V x_v     subject to ∑_(u,v) ∈E x_u x_v=0.

The constraint formulation can equivalently be written in QUBO form as

    H = -A∑_v ∈ V x_v + B∑_(u,v) ∈E x_u x_v,

where the penalty weights can be chosen as A=1, B=2 <cit.>. One disadvantage of (<ref>) lies in the fact that H contains O(n^2) quadratic terms for sparse graphs (of O(n) edges), limiting the size of problems that can fit DW.

An alternative QUBO formulation of <cit.> assumes the clique size K ≥ 1 is known. Its Hamiltonian

    H_K = A ( K - ∑_v ∈ V x_v )^2 + B(K2  - ∑_(u,v) ∈ E x_u x_v )

is designed to attain the value zero only for an assignment defining a clique of size K, which is achieved if one chooses weights A=K+1 and B=1  <cit.>.



 §.§ Problem decomposition





In <cit.> we developed several algorithms for decomposing an input graph too large for DW into subgraphs small enough that the (Maximum Clique) QUBO corresponding to each generated subproblem fits DW. After solving each subproblem, the solutions are combined again into a solution to Maximum Clique for the original graph.

We describe here an improved version of the most universal of those algorithms, in the sense that it can be applied to all input graphs regardless of their structure. It works as follows. Given a graph G=(V,E), pick any vertex v, e.g, with the smallest number of neighbors as in <cit.>, and define two new graphs: a graph G_1 induced by all neighbors of v (but without v itself) and the graph G_2 induced by V∖{v}, that is G without v and its incident edges (Figure <ref>). Then find recursively the maximum cliques C_1  of G_1 and C_2 of G_2. If |C_1|>|C_2|, return as a maximum clique for G the subgraph induced by C_1∪{v}, otherwise return C_2. The recursion stops when G is small enough to fit DW.

It is easy to see the recursion is finite, as both G_1 and G_2 contain no more than |V|-1 vertices each, since neither of them contains v, and that it correctly finds a clique of maximum size for G. Unfortunately, in the worst case it may generate an exponential number of subgraphs (exponential in the size of G). We will study in the next subsection how persistency analysis can help to reduce the number of subgraphs generated.



 §.§ Results with Qpbo

In our experiments, we used the C++ package Qpbo of <cit.> for finding persistencies, developed for computer vision applications. We access Qpbo through a Python interface that we adapted by  modifying the PyQpbo package <cit.>. This section presents results on the effectiveness of the method for Maximum Clique. We investigate the application of Qpbo to QUBOs (<ref>) and (<ref>) on random as well as special graphs (c-fat and Hamming graphs), both without and in connection with the graph splitting algorithm of Section <ref>. Moreover, we compare the two QUBO formulations (<ref>) and (<ref>) presented in Section <ref>.



  §.§.§ Qpbo on c-fat and Hamming graphs





We apply Qpbo to the QUBO formulation (<ref>). The test graphs we use are from the 1993 DIMACS Challenge on maximum cliques, coloring and satisfiabilty <cit.> and have also been used in <cit.>. Both graph families depend on two parameters: the number of vertices n and an additional internal parameter, precisely the partition parameter c for c-fat graphs and the Hamming distance d for Hamming graphs. We use the generation algorithms of <cit.> for both graph families.

Table <ref> shows the results (where both parameters c and d are summarized as a generic parameter q). We see that, in six out of eight cases, Qpbo can find persistencies for 100% of the variables, while in the remaining two cases it cannot find any persistency. Also we observe that probing is the most effective algorithm, while finding strong persistencies is the least effective.




In order to observe and compare a wider range of percentages of solved variables with Qpbo, we insert (or delete) edges from c-fat or Hamming graphs with a certain probability p which we will vary.

Figure <ref> (left) shows that when adding random edges (with edge probability p) to the (500,2) c-fat graph, no strong and weak persistencies are found for both the original graph as well as for almost all dense modified graphs (since the x-axis starts at p=0.97). Persistencies (as a percentage of solved variables out of the n=500 variables in the QUBO) are only observed for p ≥ 0.97 and increase to 100% over a very short range. As expected, weak persistencies are found (slightly) earlier than strong ones. The behavior for the Hamming graph with n=8 and Hamming distance d=4 is similar.



Similarly to Figure <ref>, Figure <ref> shows strong and weak persistencies as a function of p for the Hamming graph with n=8 and d=4. However, the behavior of strong and weak persistencies is different in this case. Figure <ref> shows that Qpbo is able to assign a value to 100% of all variables in the QUBO of the original Hamming graph via weak persistencies (see Table <ref>). Nevertheless, strong persistencies are not found. We therefore add edges with probability p to observe a progression of strong persistencies (Fig. <ref>, left) and remove edges with probability p for weak persistencies (Fig. <ref>, right).

Figure <ref> shows that strong persistencies increase gradually  over the entire interval p ∈ [0,1], whereas the found weak persistencies disappear almost instantly as soon as edges are removed from the graph.



  §.§.§ Comparison of alternative QUBO formulations for Maximum Clique

The QUBO formulation (<ref>) was used for all previous investigations. This is due to the fact that, in contrast to (<ref>), the former does not need the (usually unknown) clique size K as an input parameter. Nevertheless, this section investigates whether the choice of formulation has an effect on the persistencies found and whether changing the QUBO formulation can result in more or fewer persistencies. If this is the case, then it may be worth investing time to search for formulations that are more suitable for the persistency algorithm.

In order to find the clique size that is needed for (<ref>), we first run an alternative algorithm to find the maximum clique size K for each tested graph. We then compute the QUBO coefficients for (<ref>) with parameter K and identify strong and weak persistencies using Qpbo.

Table <ref> shows results for a Hamming graph with parameters n=8 and d=2 and a c-fat graph with parameters n=200 and c=1. As can be seen from the table, formulation (<ref>) results in much larger QUBO problems (which are likely to also have a more complex structure, thus being more difficult to solve for Qpbo) and thus in fewer persistencies compared to (<ref>), so clearly the latter formulation is preferable, at least for the type of graphs tested.





  §.§.§ Combining Qpbo with the decomposition algorithm


We apply the graph splitting algorithm described in Section <ref> to solve the Maximum Clique problem on graphs of various sizes. However, before splitting a graph, we look for persistencies (in particular, for strong and weak persistencies since applying probing in each time step took prohibitively long) and, if found, reduce the size of the graph by removing the corresponding variables. We stop splitting when all generated (sub-)graphs are of size small enough to fit DW, which in our case is 45 vertices.

We test the resulting algorithm on random graphs with n=500 vertices and an expected number of edges in the interval [10000,40000].



Figure <ref> visualizes the proportional reduction in the number of generated subproblems. For this, we plot the ratio (n_qpbo-n_no-qpbo)/n_no-qpbo as a function of p, where n_qpbo and n_no-qpbo denote the number of solver calls when employing graph splitting with and without Qpbo, respectively. The figure shows that, especially for sparse graphs, substantial reductions with respect to the number of generated subproblems can be achieved.






§ THE MAXIMUM CUT PROBLEM


This section takes a closer look at the Maximum Cut problem introduced in Section <ref>. We again state a QUBO formulation for Maximum Cut and present simulation results on g and U graphs, which have already been defined and previously been used for benchmarking Maximum Cut problems in <cit.> and <cit.>.



 §.§ QUBO and Ising formulation

We use the QUBO formulation for Maximum Cut of <cit.> who show that the value f of any cut in a graph G=(V,E) is given as

    f(x) = ∑_(u,v)∈ E( x_u (1-x_v) + (1-x_u) x_v ),

where x_u ∈{0,1} for all u ∈ V is the indicator signalling which side of the cut vertex u belongs to. Maximizing (<ref>) over x ∈{0,1}^|V| leads to a maximum cut.

The formulation (<ref>) can be simplified when switching the input range of x_u to an Ising model. In (<ref>), each term x_u (1-x_v) + (1-x_u) x_v is zero if and only if u and v belong to the same side of the cut, and one otherwise. In an Ising formulation, that is when x_u ∈{-1,+1} for all u ∈ V, the same term can be expressed as (1-x_u x_v)/2, leading to the Hamiltonian

    H_0(x) = ∑_(u,v)∈ E1/2 (1-x_u x_v),

which has to be maximized (just like (<ref>)). Equivalently, we can minimize

    H(x)    = - ∑_(u,v)∈ E (1-x_u x_v)= -|E| + ∑_(u,v)∈ E x_u x_v 
       ∼∑_(u,v)∈ E x_u x_v,

since |E| is constant. In the following we will solely use the simpler Ising formulation (<ref>).



 §.§ Ising on Lg and Lg graphs




We test Qpbo on the Ising formulation (<ref>) for the Maximum Cut problem applied to g and U graphs. Table <ref> shows the results. As seen in the table, some  weak persistencies are found, but no strong ones. Unlike the results for Maximum Clique, where the found weak persistencies are either 0% or 100%, here we observe more varied results. Again, probing is the most effective algorithm. Our results for Maximum Clique are roughly similar to the ones presented in <cit.> (who, however, use a different and possibly more elaborate algorithm than Qpbo, which they did not make available).

Additionally, we again studied how the performance of Qpbo varies as the graph becomes denser by adding or removing random edges. These results are similar to the ones already described for Maximum Clique (Figures <ref> and <ref>) and are thus not presented here.






§ CONCLUSIONS


This paper studies the use of preprocessing methods to reduce the size of QUBO and Ising integer programs. Such models have recently attracted increased attention since they allow constant time solutions on newly available quantum annealers. However, the small number of available qubits on such devices severely limits the range of solvable problems. This motivates the use of preprocessing methods to reduce the number of required qubits, thus extending the application range of such devices.

We use the package Qpbo <cit.> to compute strong and weak persistencies in integer programs and test the suitability of preprocessing methods on two important NP-hard graph problems, the Maximum Clique and the Maximum Cut problems.

Our main findings can be summarized as follows:

  
  * The observed reduction in variables seems to be very problem-specific, both in terms of the QUBO/ Ising formulation considered as well as the employed test graphs. Test graphs on which Qpbo achieved 100% reductions for Maximum Clique did not lead to any reduction for Maximum Cut or vice versa. We observe 0% and 100% reductions considerably more often than reductions between those extremes.
  
  * When applied to a decomposition algorithm which divides up a maximum clique computation on graphs of arbitrary size into many smaller subproblems, Qpbo can be used to reduce the number of generated subproblems. This reduction seems especially significant for sparser graphs.
  
  * Using two different QUBO formulations for Maximum Clique, we show that the performance of Qpbo is not only dependent on the type of problem solved and the input graph, but also on the specific QUBO formulation used. This, therefore, justifies research into finding formulations for different optimization problems that are more suitable for identifying persistencies. This is consistent with the much more general observation in combinatorial optimization practice that the choice of formulation for a given optimization problem can have a huge effect on the time it will take to solve the resulting problem using popular solvers.


The effectiveness of Qpbo varied widely between seemingly similar problem instances. An especially interesting and practically important open problem is to be able to characterize input graphs and problem formulations that are more suitable for identification of persistencies.

Another interesting open problem is how to best speed up the persistency finding algorithms. For instance, in the case where persistency finding (via Qpbo) is combined with the decomposition algorithm (for Maximum Clique), we currently run Qpbo on every subgraph generated during the graph splitting process. A better way is to try to use information collected during the persistency analysis for the larger graph in order to infer persistencies for the smaller ones, rather than computing them from scratch, thereby saving computation time.







[Balas and Yu, 1986]Balas1986
Balas, E. and Yu, C. (1986).
Finding a maximum clique in an arbitrary graph.
SIAM J Comp, 15:1054–1068.

[Boros and Hammer, 1991]BorosHammer1991
Boros, E. and Hammer, P. (1991).
The max-cut problem and quadratic 0-1 optimization; Polyhedral
  aspects, relaxations and bounds.
Ann Oper Res, 33:151–180.

[Boros and Hammer, 2001]BorosHammer2001
Boros, E. and Hammer, P. (2001).
Pseudo-Boolean Optimization.
Rutcor Report, pages 1–83.

[Boros et al., 2006]Boros2006
Boros, E., Hammer, P., and Tavares, G. (2006).
Preprocessing of Unconstrained Quadratic Binary Optimization.
Rutcor Research Report, RRR 10-2006:1–58.

[Bunyk et al., 2014]Bunyk2014
Bunyk, P., Hoskinson, E., Johnson, M., Tolkacheva, E., Altomare, F., Berkley,
  A., Harris, R., Hilton, J., Lanting, T., Przybysz, A., and Whittaker, J.
  (2014).
Architectural considerations in the design of a superconducting
  quantum annealing processor.
IEEE Trans on Appl Superconductivity, 24(4):1–10.

[D-Wave, 2016]dwave2016
D-Wave (2016).
Introduction to the D-Wave quantum hardware.

[Djidjev et al., 2017]Djidjev2017
Djidjev, H., Chapuis, G., Hahn, G., and Rizk, G. (2017).
Finding maximum cliques on a quantum annealer.
ACM Computing Frontiers, 1(1):1–8.

[Hammer et al., 1984]Hammer1984
Hammer, P., Hansen, P., and Simeone, B. (1984).
Roof duality, complementation and persistency in quadratic 0-1
  optimization.
Mathematical Programming, 28(2):121–155.

[Hasselberg et al., 1993]Hasselberg1993
Hasselberg, J., Pardalos, P., and Vairaktarakis, G. (1993).
Test Case Generators and Computational Results for the Maximum
  Clique Problem.
Journal of Global Optimization, 3:463–482.

[Johnson and Trick, 1996]johnson-trick-96
Johnson, D. and Trick, M., editors (1996).
Clique, Coloring, and Satisfiability: Second DIMACS
  Implementation Challenge, DIMACS, volume 26.
American Mathematical Society.

[Johnson et al., 2011]Johnson2011
Johnson, M., Amin, M., Gildert, S., Lanting, T. Hamze, F., Dickson, N., Harris,
  R., Berkley, A., Johansson, J., Bunyk, P., Chapple, E., Enderud, C., Hilton,
  J., Karimi, K., Ladizinsky, E., Ladizinsky, N., Oh, T., Perminov, I., Rich,
  C., Thom, M., Tolkacheva, E., Truncik, C., Uchaikin, S., Wang, J., B., W.,
  and Rose, G. (2011).
Quantum annealing with manufactured spins.
Nature, 473:194–198.

[Kim et al., 2001]Kim2001
Kim, S.-H., Kim, Y.-H., and Moon, B.-R. (2001).
A Hybrid Genetic Algorithm for the MAX CUT Problem.
Proceeding GECCO'01 Proceedings of the 3rd Annual Conference on
  Genetic and Evolutionary Computation, pages 416–423.

[King et al., 2015]King2015
King, J., Yarkoni, S., Nevisi, M. M., Hilton, J. P., and McGeoch, C. C. (2015).
Benchmarking a quantum annealing processor with the time-to-target
  metric.
arXiv:1508.05087, pages 1–29.

[Kolmogorov, 2014]Kolmogorov2014
Kolmogorov, V. (2014).
QPBO v1.4 software package.
C++ package manual.

[Lucas, 2014]Lucas2014
Lucas, A. (2014).
Ising formulations of many np problems.
Frontiers in Physics, 2(5):1–27.

[Müller and Behnke, 2014]JMLR:v15:mueller14a
Müller, A. C. and Behnke, S. (2014).
Pystruct - Learning Structured Prediction in Python.
Journal of Machine Learning Research, 15:2055–2060.

[Rother et al., 2007]Rother2007
Rother, C., Kolmogorov, V., Lempitsky, V., and Szummer, M. (2007).
Optimizing Binary MRFs via Extended Roof Duality.
Proceedings CVPR, IEEE Computer Society Conference on Computer
  Vision and Pattern Recognition, pages 1–15.


