Verification of Quantum Optimizers

# Verification of Quantum Optimizers

## Abstract

Methods for finding ground states of classical spin models are of great importance in optimization and are gaining additional relevance now for verifying the results quantum optimizers. We combine the state-of-the-art branch and bound method for solving such optimization problems via converging upper- and lower-bounds with ideas from polynomial optimization and semidefinite programming (SDP). The resulting chordal branch and bound (CBB) algorithm can exploit the locality and resulting sparsity in relevant Ising spin models in a systematic way. This yields certified solutions for many of the problems that are being used to benchmark quantum annealing devices more efficiently and for larger system sizes. We are able to verify the output of a D-Wave 2000Q device for the largest triangular lattice that can be embedded in the hardware and provide exact ground states for cases in which the quantum annealer returns a configuration with almost minimal energy but markedly different spin pattern. The method always yields provable polynomial time upper and lower bounds on the ground state energy. We benchmark our method against further planar and non-planar graphs and show that these bounds often converge after a small number of steps, even though the NP-hardness of general Ising models implies that exponentially many steps are required in the worst case. This new tool is a flexible and scalable solution for the verification and benchmarking of the next generation of quantum optimization devices.

## I Introduction

Classical Ising models are among the most paradigmatic and widely studied models in statistical physics. They are capable of describing an immense variety of interesting physics, ranging from ferromagnetic to frustrated and glassy phases. Moreover, they are important in fields as diverse as risk assessment in finance, logistics, machine learning Koller et al. (2007), and image de-noising Bhatt and Desai (1994). Many optimization and decision problems, such as partitioning, covering, and satisfiability can be mapped to such models Lucas (2014). The generality and exponentially growing configuration spaces of such models, however, precludes the existence of any efficient general purpose algorithm. It is hence no surprise that a wealth of approximate but more scalable classical techniques for the energy minimization in such models has been developed. Recently, novel approaches that leverage the power of near-term quantum devices such as quantum annealers, variational quantum eigensolvers Peruzzo et al. (2014), or networks of degenerate optical parametric oscillators QNNCloud (2017), are proposed for performing such tasks Denchev et al. (2016); Crosson and Harrow (2016).

The most scalable of the classical approaches based on simulated annealing Kirkpatrick et al. (1983) or variational ansatz classes, as well as the mentioned quantum mechanical methods, however, all have one thing in common: They only yield upper bounds on the ground state energy. A comparison of such upper bounds can provide a ranking of the methods’ performance for benchmark tasks. However, it cannot guarantee that the obtained upper bounds are actually close to the true ground state energy and, even if this is the case, there is no guarantee that the final spin configuration is anywhere near that of the true ground state, or just corresponds to a far away local minimum. It is thus a problem of great importance to develop schemes that provide lower bounds to the ground state energy, against which the results of such widely used upper-bound techniques can be compared.

The state-of-the art methods for doing this are the so-called branch and bound (BB) methods (see for instance Rendl et al. (2008) for a recent summary of these techniques). These methods generally yield both lower and upper bounds to ground state energies of classical spin models and they have been shown to yield exact solutions in reasonable time for many classes of models. In this work, we present a BB technique built upon a semidefinite programming (SDP) relaxation of the ground state energy problem and augmented with a trick known as the chordal extension Waki et al. (2006). In this way we obtain a more efficient classical algorithm that can find lower bounds to ground state energies as well as exact ground states of spin models of interest in physics. We will refer to our method as chordal branch and bound (CBB). As we will demonstrate, the CBB method can, for example, find exact ground states of disordered Ising models close to a phase transition with a runtime that scales substantially more favourable with the number of spins while at the same time using much less memory than standard SDP-based BB techniques. Another advantage of the CBB method, compared to, for instance variational methods, is that it can be made to directly yield physically interesting quantities, such as correlation functions or expectation values, instead of merely a prescription of an approximation to the ground state from which such quantities then need to computed separately.

The most timely application of our technique is the verification of quantum annealing devices such as the D-Wave machine. The complexity of the Ising ground state problem can change drastically when seemingly unimportant details of a given model are modified, which must be taken into account when benchmarking such devices (for details see Section III.2). In particular, some of the models that have been used in benchmarking Rønnow et al. (2014); Mandrá et al. (2017); Mandrá and Katzgraber (2018) of quantum annealing devices, while being difficult to solve with classical simulated annealing, can actually be provably solved exactly in polynomial time on a classical computer. We demonstrate that our CBB method is indeed capable of doing this, sometimes for the largest system size simulable on current quantum annealing devices. We also find that the quantum annealer is not always successful and can get stuck in local minima, but outperforms the CBB method in terms of absolute runtime.

How can one find the lowest energy configuration among the more than possible configurations of a typical Ising model with, say sites? The key is to compute a converging series of upper and lower bounds on the ground state energy and use these bounds to drastically reduce the search space that needs to be explored. Configurations that fall outside the band between the bounds no longer need to be explored. This is precisely the strategy of the BB approach (see from Figure 1). Obtaining good lower bounds is the non-trivial part and several proposals have been considered through the last yeas: One can generally divide them into techniques based on a linear-programming relaxation Liers et al. (2004) and ones based on semidefinite programming Rendl et al. (2008). The former approaches are known to perform better for very sparse graphs, such as 2D square lattices and planar graphs, but it is generally hard to adapt them to deal efficiently with more general problems. On the other hand, the SDP method introduced in Rendl et al. (2008) applies to general graphs and is always polynomial in the system size (per step), but does not outperform linear programming for very sparse graphs. The CBB keeps the polynomial dependence on the size, while being both more efficient than the BB approach from Rendl et al. (2008) for sparse instances and readily applicable to less sparse Ising models. It is therefore a good candidate for the verification of current and near future quantum optimization devices.

Due to the combinations of bounds from below and above, any BB technique ultimately converges to the lowest achievable energy and outputs a corresponding configuration that achieves it. It thereby gives an exact and certified solution to the ground state problem. Although the complexity of the problem implies that, at least in some cases, one has to run an exponential amount of steps to achieve convergence, in many scenarios convergence is reached after a polynomial number of steps and then the method is efficient (the complexity per step and number of required steps can be tuned by adjusting the level of the relaxation, as we will see below). Even in cases where convergence cannot be achieved with the available computational resources, the obtained lower and upper bounds still imply a interval in which the true ground state energy must provably lie. This fundamentally differentiates both BB and CBB from techniques such as simulated annealing or variational methods, which can reach larger system sizes, but provide not strong certificate for the solutions they propose.

## Ii Results

We first describe the results that we could obtain by combining the branch and bound approach with the chordal extension, before describing the CBB method in more detail in Section III. We start with a comparison of the standard SDP-based BB methods and CBB and then move on to the more interesting problem of verifying the results of two kinds of energy minimizations performed on an actual quantum annealing device. The numerical minimization with CBB and BB was run on a workstation with an Intel Xeon E5-1650v4 processor with six physical cores clocked at 3.60 GHz base frequency and 128 GByte RAM. Due to the polynomial scaling of the method, much larger system sizes can be reached with more powerful hardware. The sparse semidefinite relaxations were generated by Ncpol2sdpa Wittek (2015), and the semidefinite programs were solved by Mosek Mosek (2010). The code for the experiments is available under an open source license 1.

### ii.1 Ising model on a 2D square lattice

As a first test we compare the performance of the CBB method with a sparse Ising problem in the two cases of exploiting and not exploiting the chordal extension trick. In the absence of chordal extension, our BB method is comparable to the one introduced in Rendl et al. (2008) (for details, see Section III). As a benchmark of a sparse instance, we consider the standard 2D ferromagnetic Ising model in a statically disordered magnetic field (quenched disorder) that is picked independently from normal distributions of mean zero and variance for each site. As a function of the disorder strength , the model undergoes a phase transition from a ferromagnetic ground state (in which all states are aligned with each other) to a disordered phase (in which, for extremely large disorder, the spins are aligned with the local magnetic fields). For this model it is known that the ground state can in principle be found in polynomial time (for details see Section III.2).

Indeed, the non-chordal BB method is able to do that, but, especially in the interesting region close to the phase transition, fast growing memory requirements and runtime make the method impractical for systems that are larger than on the hardware we have at our disposal. The CBB method, in contrast, allows us to solve systems of over sites on the same hardware, due to both lower memory requirements and a very significantly reduced runtime, both in terms of absolute numbers and in terms of scaling (see Figure 2 for a comparison). While the runtime of non-chordal BB exhibits a roughly dependence on the number of spins , CBB scales roughly as .

### ii.2 Verifying the solution of a D-Wave quantum annealer

We now turn to the verification of a ground state energy search by quantum annealing. To show the flexibility of the CBB method and also to verify a quantum annealing solution for the largest system size simulable on a state-of-the art annealer, we move to a more interesting and slightly less spares lattice, namely the 2D triangular lattice. Spin models on the triangular lattices display a wealth of interesting physical phenomena, many driven by the possibility to have frustrated interactions. To remain in a regime that is comparable to the benchmarking we did before, we however concentrate on the interplay of ferromagnetic interactions with a disordered magnetic field.

We used a D-Wave 2000Q quantum annealer with 2040 functional qubits. The chip had 8 faulty qubits and the corresponding couplings were removed from a full-yield Chimera graph. We used the virtual full-yield Chimera graph abstraction to ensure consistent embeddings and improve the quality of the results. The coupling strengths were automatically scaled to the interval , and the logical qubits used a coupling strength of to hold a chain of physical spins together. The minor embedding was a heuristic method, yielding a chain length of . We also tried chains up to length , without significant change in the results, showing that the scaling in the couplings ensures that the chains do not break. For each data point, we sampled a thousand data points and chose the one with the lowest energy as the optimum. This takes constant time irrespective of the values set, in the range of milliseconds. The flux bias of the qubits was not offset.

Both the quantum annealer and the CBB simulation were done for the same disorder realizations (the disorder in the annealer is fully programmable) to obtain directly comparable results. We observe that for most disorder realizations, we can verify by CBB that the quantum annealer is able to find the exact ground state energy. This is true even for intermediate disorder strength, where the ground state spin pattern shows macroscopic islands of aligned spins whose precise shape and positions depends non-trivially on the disorder realization around . It does that typically in a very short time. However, there are also cases, in which, even after 1000 repetitions, the lowest energy found by the quantum annealer is still higher than the exact value computed with CBB. Here, the optimal spin configuration found by the quantum annealer typically differs markedly from the true ground state, which shows that it gets stuck in local minima and emphasizes the importance for exact methods such as CBB for comparison and benchmarking. We summarize results of simulations in this interesting regime of intermediate disorder strength in Figure 3.

### ii.3 Verifying solutions for a Chimera graph

Lastly, we consider the application of CBB to a denser graph. For this purpose we choose the Chimera architecture Bunyk et al. (2014), which is the natural graph on the D-wave 2000Q hardware. The corresponding graph is composed of fully connected bipartite unit cells, consisting of 8 spins â 4 horizontal and 4 vertical â with edges between each horizontal/vertical pair. These unit cells are arranged to form a 2D square lattice of size and a total number of spins. Because of the in-cell connectivity, such a graph is clearly non-planar and thus has the potential to encode NP-hard Ising models. We consider again a ferromagnetic model with a disordered magnetic field and we find that the CBB method is able to successfully identify the ground state energy for lattices sizes of up to a linear size of . This corresponds to architectures involving a total number of spins, which is a remarkable advancement compared to the corresponding limit of just spins for a standard SDP-based BB method on the same hardware.

Although the D-Wave 200Q quantum annealer is currently implementing a Chimera graph with functional physical qubits, they are seldom actually used as logical spins. Most recent studies encode each cell as a single logical spin, in order to suppress errors due to the finite size and qubit quality of the system Mandrá and Katzgraber (2018). This results in effectively solving Ising models on a 2D square lattice which, being planar, is actually proven to be polynomially solvable (see also Section III.2). The numerical test was performed on the actual Chimera graph. This opens up the way to benchmarking future annealing devices, once their physical qubit quality has improved to a point that makes the individual spins useful, in the much more interesting regime of non-planar graphs.

## Iii Methods

### iii.1 Setting and notation

We consider classical spin systems whose configurations are vectors of spin variables to each of which a Hamiltonian assigns an energy . We will mostly be interested in Hamiltonians of Ising type, that can be written in the form

 H(→σ)=−∑1≤i

with couplings and local fields , but the method we develop is more general and can also be applied to Hamiltonians that are higher order polynomials of the and couple three or more spins in a single term.

Many spin models of interest in physics and beyond are characterized by a locality structure, i.e., not all the spins interact with each other so that some are zero. Such locality of interactions implies a sparsity of the resulting Hamiltonian and hence optimization problem. This locality and sparsity can be most efficiently captured in the language of graph theory: We define the interaction (hyper-)graph as the graph whose vertex set is the set of indices of the spin variables and whose (hyper-)edge set is the set of all pairs (or larger subsets) of indices for which the Hamiltonian contains an interaction, i.e., in the above example . More important for us is the dependency graph , which is a related concept but always a graph and not a hyper-graph, namely the one that contains an edge for any pair of indices of spin variables that appear together in the same Hamiltonian term. For the example Hamiltonian (1) and are identical, but in general one obtains by replacing every hyper-edge in by a clique, that is a fully connected subgraph. Typical examples of interaction and dependency graphs are regular cubic grids, such as for example in the 1D and 2D Ising model.

Among all the configurations of such a system there are those that achieve the minimal possible energy, also known as the ground state energy and defined as

 Eg\coloneqqmin→σ∈{−1,1}NH(→σ). (2)

If only one configuration achieves energy we call this configuration the ground state and say that it is unique, otherwise we call the collection of all configurations with energy the ground state space. For our purposes, solving the ground state problem for a given Hamiltonian means finding and outputting a configuration that achieves it. Obviously, the ground state problem is an optimization problem that can, in principle, be solved by brute force, by simply computing for all possible configurations . This however quickly becomes infeasible as the number of configurations grows like , restricting this approach to systems of size roughly .

Therefore, several methods have been proposed to solve the ground state problem more efficiently. Among them, we can distinguish three main approaches: the ones that lead to upper bounds to , the ones that seek for lower bounds, and the ones that yield exact solutions. To the fist family belong all techniques based on Markov chain Monte Carlo, such as simulated annealing, both in its classical Kirkpatrick et al. (1983) and quantum form Denchev et al. (2016); Crosson and Harrow (2016). Lower bounds can be obtained by considering convex relaxations of the ground state problem Lasserre (2001a, b); Pironio et al. (2010). A trick from graph theory known as the chordal extension allows exploiting sparsity in the problem to be relaxed Waki et al. (2006); Lasserre (2006); Erdogdu et al. (2017). The state-of-the-art algorithms to obtain exact solutions to the ground state problem rely on a combination of lower and upper bound techniques, following the so-called branch and bound (BB) method Rendl et al. (2008). In this work we introduce a SDP-based branch and bound method augmented with the chordal extension, which allows us to systematically exploit the locality present inherent in many relevant Ising type optimization problems and efficiently obtain exact solutions to the ground state problem in certain cases.

Before describing the techniques we combine in more detail, we first review what is known about the complexity of the ground state problem and point out the consequences this has for the benchmarking of quantum annealing devices. We then describe the relaxation of the ground state problem that we will use to obtain lower bounds on the ground state energy and how sparsity can be exploited by means of the chordal extension. We then show how this lower bound technique can be combined with a methods that yields upper bounds in a branch and bound algorithm.

### iii.2 Complexity of finding Ising ground states

There is a wealth of results on the worst case complexity of finding the ground state of various Ising models Grötschel et al. (1987). Thereby “worst case complexity” is the complexity of the hardest instances within a class of families of problems of increasing size. How hard it is to solve the ground states problem of such a family varies with the interaction graph and can crucially depend on seemingly unimportant details. We consider Hamiltonians that are polynomials (with fixed finite degree and finite precision coefficients) in the spin variables (such as those given in (1)) and interaction (hyper-)graph . The size of a problem is the number of vertices . We say that an Ising model of the form (1) has no fields if for all , we say it has an external field if all for all and some , and we say it has on-site fields if all can be chosen independently.

Finding the ground state of Hamiltonians of the form (1) for arbitrary graphs is in general NP-hard Barahona (1982), even without any fields. This is still true for and a finite 3D cubic grid graph and even for a cubic two-layer 2D grid Barahona (1982). In contrast, for planar graphs and without local fields, the ground state can be found efficiently even without the restriction  Grötschel et al. (1987). With the restriction this even holds for toroidal graphs (grids on a torus, i.e., systems with periodic boundary conditions) Grötschel et al. (1987). Similarly, if , then even some systems with local fields can be solved in polynomial time Grötschel et al. (1987). On the other hand, for general planar graphs with interactions and uniform external field , finding the ground state is again NP-hard Barahona (1982). Ref. Grötschel et al. (1987) contains a list of further concrete models whose complexity is either known to be in P or proven to be NP-hard.

These hardness results are typically obtained by reducing the ground state problem to the so-called Max-Cut Problem, which is known to the NP-hard. The polynomial time algorithms to solve the other families of systems, in turn, work by finding perfect matchings Barahona (1982) or rely on so-called max-flow/min-cut methods Grötschel et al. (1987).

### iii.3 Ingredients for the chordal branch and bound method

In this section we describe in detail the methods we build upon and the resulting chordal branch and bound algorithm. We begin by introducing additional notation. A state of an spin system is a probability distribution over the set of configurations . For every function we can then define its expectation value in the state as

 ⟨f⟩P=∑→σ∈{−1,1}Nf(→σ)P(→σ). (3)

We call any state that is supported only on the ground state space of a model a ground state and such are manifestly those that achieve the minimal possible expectation value for the Hamiltonian, i.e, .

#### A hierarchy of SDP relaxations

We first discuss how relaxations of the ground state problem can be used as a means of obtaining lower bounds on the ground state energy. Specifically we will be using a method pioneered by Lassere Lasserre (2001a, b), that yields a hierarchy of tighter and tighter lower bounds to the ground state energy. To make this work self-contained, we revise the general method by using the notation here, instead of introducing it in the more typical framework of polynomial optimization.

Let us consider a vector of monomials of the spin variables. For such vector and state we define its moment matrix as the matrix of expectation values . We can now generalize this concept and think of any real symmetric matrix as an assignment to expectation values of monomials , that need not necessarily be achievable by any physical state . In particular, it is not hard to see that for any state , the moment matrix is positive semidefinite, i.e., , and that, depending on what the elements of are, it further obeys certain linear constrains. The constraints reflect the two basic properties of classical spin variables of taking dichotomic values and commuting with each other. It is not hard to see that these properties imply conditions on the expectation values such as, for example, .

If the vector contains at least a suitable subset of polynomials of spin variables, one can further compute the expectation value of the Hamiltonian from the entries of the moment matrix in the sense that there is a matrix (depending on the and in case the Hamiltonian is of the form in (1)), such that for any physical state it holds that . If this is the case, one can relax the ground state problem by computing the energy in this way and optimizing over all matrices that are positive semidefinite and fulfil the above mentioned linear constraints, rather than over those that can actually arise from a physical state .

A systematic way of constructing a hierarchy of such relaxations is as follows: Let be the vector of all monomials of spin variables of degree up to . Then (for any sufficiently large so that , implicitly defined via , exists), the minimum achievable in the minimization problem

 E(ν)g=minΓ(ν)tr(hΓ(ν)) (4) s.t.Γ(ν) ⪰0, (5) ∀m∈{1,…,k}:tr(FmΓ(ν)) =0, (6)

where the represent the linear constrains mentioned before, is a lower bound on the true ground state energy, i.e., . It is further obvious that the are ordered in the sense that for any . Remarkably, if all the relevant are taken into account, the bounds actually converge to the true ground state energy for any fixed finite system size and Hamiltonian , in the sense that  Lasserre (2001b); Pironio et al. (2010).

We illustrate this with an example. If we are working at level for a system of spins, the corresponding matrix take the following form:

 Γ(2)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1⟨σ1⟩⟨σ2⟩⟨σ3⟩⟨σ1σ2⟩⟨σ1σ3⟩⟨σ2σ3⟩⟨σ1⟩1⟨σ1σ2⟩⟨σ1σ3⟩⟨σ2⟩⟨σ3⟩⟨σ1σ2σ3⟩⟨σ2⟩⟨σ1σ2⟩1⟨σ2σ3⟩⟨σ1⟩⟨σ1σ2σ3⟩⟨σ3⟩⟨σ3⟩⟨σ1σ3⟩⟨σ2σ3⟩1⟨σ1σ2σ3⟩⟨σ1⟩⟨σ2⟩⟨σ1σ2⟩⟨σ2⟩⟨σ2⟩⟨σ1σ2σ3⟩1⟨σ2σ3⟩⟨σ1σ3⟩⟨σ1σ3⟩⟨σ3⟩⟨σ1σ2σ3⟩⟨σ1⟩⟨σ2σ3⟩1⟨σ1σ2⟩⟨σ2σ3⟩⟨σ1σ2σ3⟩⟨σ3⟩⟨σ2⟩⟨σ1σ3⟩⟨σ1σ2⟩1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (7)

The expectation value of any Hamiltonian of the form given in (1) can be expressed as a function of the entries of as . Similarly, one can see how the linear constraints on the entries reflect the properties of the spin variables. Dichotomity directly imposes the conditions on the diagonal variables and, combined with commutations, allows to identify some of the entries, such as for example .

#### Exploiting sparsity via the chordal extension

Recall that, depending on the kind of system considered, the optimisation problem defined by the Hamiltonian (1) can be sparse. As is shown in Refs. Waki et al. (2006); Lasserre (2006), we can exploit this sparsity to derive a more scalable relaxation than (4). Intuitively, the idea behind the modification is the following: for any pair of non-interacting spins , the corresponding two-body expectation value is not needed for the computation of the energy. Thus, a moment matrix with all two-body correlations is including some potentially unnecessary information. Hence, finding the minimal amount of such information that is sufficient to perform the optimisation of the energy, helps define a more efficient relaxation.

The method works as follows: take the dependency graph of the problem and check if it is chordal. A graph is said to be chordal if all its cycles of four or more vertices have a chord, i.e. an edge that is not part of the cycle but connects two vertices of the cycle. If is not chordal, construct a so-called chordal extension of by suitably adding edges until the graph is chordal (see Figure 4 for an example). The chordal extension of a graph is not unique, but a chordal extension can always be found, simply because any fully connected graph is chordal. The challenge is to find a chordal extension that is still relatively sparse. A method that works well in this respect for all the cases studied here is to compute an approximate minimum degree ordering of the graph nodes, followed by Cholesky factorization Vandenberghe and Andersen (2015).

Once a specific chordal extension is constructed, it will contain a number of maximal cliques . A clique, that is a fully connected subgraph, is maximal if it cannot be extended by including any other adjacent vertex. Since the graph represents a sparse Hamiltonian, and is obtained from by simply adding some edges, the function (1) can be decomposed into a sum of terms that each contain only variables contained in a given maximal clique .

One can then modify the optimization problem in (4) as follows: replace the big matrix by a direct sum of smaller matrices , one for each clique, constructed from the spin variables belonging to the clique . Some spin variables appear in more than one clique, which can be captured with additional linear constraints that involve variables of the blocks . Writing this explicitly, the chordal version of the SDP relaxation then reads

 minΓ(ν)l ∑ntr(hnΓ(ν)n) (8) s.t. Γ(ν)l⪰0,∀l=1,…nC, (9) tr(Fm,lΓ(ν)l)=0m=1,…,kl,l=1,…nC, (10) ∑ntr((Gl,n)Γ(ν)n)=0l=1,…,k, (11)

where the are the intra-block constraints coming from the properties of the spin variables, while the correspond to the constraints identifying expectation values belonging to different blocks. Interestingly this relaxation still converges to the exact result, as was shown in Lasserre (2006). Depending on the sparsity of the graph (and its chordal extension ), substituting the original optimization relaxation (4) – (6) by (8) – (11) leads to a substantial simplification and improved scaling in runtime and memory. In practical applications, the latter are typically dominated by the the largest block, i.e., the largest maximal clique in .

Let us illustrate how this chordal extended relaxation works in practice, by going back to the three spins example introduced in the previous subsection. Imagine one wants to solve the 1D Ising model with Hamiltonian . The corresponding dependency graph is already chordal and is composed of two cliques and . Then, for a relaxation at level the matrix (7) can be substituted by the two blocks:

 Γ(2)C1 =⎛⎜ ⎜ ⎜ ⎜⎝1⟨σ1⟩⟨σ2⟩⟨σ1σ2⟩⟨σ1⟩1⟨σ1σ2⟩⟨σ2⟩⟨σ2⟩⟨σ1σ2⟩1⟨σ1⟩⟨σ1σ2⟩⟨σ2⟩⟨σ1⟩1⎞⎟ ⎟ ⎟ ⎟⎠ (12) Γ(2)C2 =⎛⎜ ⎜ ⎜ ⎜⎝1⟨σ2⟩⟨σ3⟩⟨σ2σ3⟩⟨σ2⟩1⟨σ2σ3⟩⟨σ3⟩⟨σ3⟩⟨σ2σ3⟩1⟨σ2⟩⟨σ2σ3⟩⟨σ3⟩⟨σ2⟩1⎞⎟ ⎟ ⎟ ⎟⎠ (13)

The constraints derive from that the variable belongs to both cliques, hence many expectation values appear in both blocks. Some of the unnecessary expectation values in (7), such as and no longer appear in the two smaller blocks and . Such a simplification is particularly useful because it reduces the number of variables involved in the SDP. Although the constraints do not allow to split the problem into independent ones, one can still see that the scaling of the computational effort is dominated by the size of the largest block alone.

#### Branch and bound techniques to solve the ground state energy problem

The state-of-the-art techniques to compute exact solutions to the ground state problem rely on the BB method. This is a general iteration strategy that has been applied in several different ways (see, for instance, Ref. Rendl et al. (2008) for a review). However, all different applications build upon three main ingredients:

1. Lower bound: a method to find a lower bound to the ground state energy. This generally implies the use of some relaxation of the ground state energy problem.

2. Upper bound: some heuristics that gets a spin configurations that constitutes an approximation from above to the ground state energy. A simple example would be classical annealing.

3. Branching procedure: a branching consist in dividing the original problem into two sub-problems that corresponds to the opposite cases of a dichotomic choice. In the ground state search, it can be done by choosing a vertex and considering the two sub-sets of spin configurations that have fixed. Solving the ground state problem in both of the sub-cases can be cast as another ground state problem for a modified graph where the vertex has been removed and the couplings have been modified accordingly. Then, clearly, the solution of the original ground state energy is just the minimum between the solutions of the two sub-cases.

The trick is now to use the upper and lower bounds to reduce the number of branches to explore. The BB procedure does that as follows: (i) start with the original graph and compute a lower and upper bound to the ground state energy; (ii) if the bounds differ, choose a branching and compute the lower bounds for the two cases; (iii) pick the branch with the lowest an compute the corresponding . Discard all branches in which the lower bound is higher that the newly found ; (iv) if the lower and upper bounds still differ, go back to point (ii) and perform another branching; (v) keep repeating until the upper and lower bounds converge.

In Ref. Rendl et al. (2008), a versatile BB has been proposed, which was shown to be applicable to any ground state energy problem for any dependency graph . Interestingly, it exploits a lower bound method that is essentially equivalent to the first level of the relaxation (4), with the addition of some hand crafted linear constraints. The hierarchy discussed here, allows to systematically construct an infinite family of increasingly precise relaxations that yield better and better bounds at an increasing computational cost. This increases the computational effort per branching step, but can drastically reduce the number of steps and the overall memory requirements. In fact, we have checked numerically that all the constraints introduced in Ref. Rendl et al. (2008) are automatically implied by increasing the level of the hierarchy to . As we will see below, also an upper bound (and the corresponding spin configuration) can be obtained directly from the moment matrix resulting from the SDP, so that from just one method, namely the chordal SDP, one can obtain both the upper and lower bound, with a tunable tightness and in a systematic fashion.

### iii.4 The chordal branch and bound method

We now describe in detail the three ingredients of our chordal branch and bound (CBB) strategy: the lower bound, the upper bound, and the branching rule.

#### Lower bound

The cheapest method to get a bound on the ground state energy from below would be use the relaxation in (8) at its lowest level, namely . However, in practical applications this leads to a lower bound that can be more than away from the corresponding upper bound. As already mentioned in Ref. Rendl et al. (2008), having such a big initial gap slows down the convergence of the branch and bound, making it very difficult to reach a point in the branching where the lower bound is high enough to start excluding the first branches. This problem is overcome by tightening the relaxation, which can be done in several ways.

The option considered in Ref. Rendl et al. (2008) was to tighten the relaxation at level by adding hand crafted linear inequalities, so-called triangle inequalities, between entries of the matrix corresponding to the two-body correlations of triples of spin variables. Since the amount of all these possible constraints scales as , usually only part of them is introduced. In our CBB we can exploit the structure of the problem to obtain more systematic efficient improvements.

The chordal extension reduces the amount of meaningful constraints that can be added. Indeed, the resulting block structure implies that the only two-body expectation values that appear in the moment matrices correspond to spins belonging to the same block. Hence, all the triangle inequalities that can actually be imposed have to involve triples that appear in the same clique.

The numerical effort for one step in the CBB method is mostly determined by the largest block in the moment matrix. Therefore, we choose to take a hybrid approach, introducing an intermediate level with for all the blocks involving less than variables, while keeping all the bigger blocks at level . This devises a test that corresponds at least to the case of level plus the addition of all triangle inequalities between variables in the smaller blocks. Taking such a hybrid level yields a significant improvement in the initial lower-upper bound gap already for smaller values of .

Moreover, we also allow for additional triangle constraints between two-body correlations belonging to bigger blocks. In particular, we add them in an iterative way, as shown in Ref. Rendl et al. (2008), until the improvement on the lower bound is smaller than some numerical precision. In most cases we tested, there was actually no need to introduce these additional constraints.

#### Upper bound

For the upper bound one needs a good guess for a spin configurations that is close to the ground state energy. Here we develop an improvement over the method proposed in Ref. Rendl et al. (2008). A nice feature of the method is that it extracts an upper bound directly from the moment matrix that is obtained by solving the SDP to get the lower bound. Intuitively, this can be seen as a way to obtain the spin configurations “closest” to the optimal (but typically unphysical) solution achieved by the relaxation.

Our method has a very simple interpretation: assign the deterministic value to the spin corresponding to the sign of the expectation value as obtained by the optimization performed to obtain the lower bound. Let us explain it in more detail: take to be the moment matrices resulting from the lower bound optimization for clique . This matrix always contains the sub-matrix , involving the variables occurring already at level , namely that generated by the set of variables . is a matrix, where stands form the number of spins contained in the clique (compare also the examples in Eqs. (12) and (13)).

Because they are positive semi-definite, these matrices can be identified with a collection of vectors in a dimensional space that reproduce their entries by the scalar product . A way to obtain such vectors is to perform a Cholesky decomposition of the moment matrices and reduce the columns of the resulting up to the point where can still be recovered up to the desired numerical precision.

Inspired by a similar argument used in Ref. Navascués et al. (2008), we interpret the first vector of each block as a representation of a state of the spin systems, while the remaining vectors represent the spin variables in the clique . Indeed, one can see that the first row of the moment matrices recovers the spin expectation values . Now we can define a spin configuration for the clique from the sign of these scalar products, i.e., we set the variables of the th spin in the th subset to be .

This approach is different from what had been done in earlier approaches. In Ref. Rendl et al. (2008), the moment matrix is missing the first row and column vector, and thus exactly the entries we need. In this case, a configuration had to be extracted via the of inner products between and randomly chosen vectors . Constructing spin configurations via random vectors in this way is incompatible with the chordal extension, since the cliques of a graph overlap on some vertices, it will happen that the same spin index appears in more than one list . Therefore, given that the aim it to extract a complete and consistent spin configuration , one must avoid having conflicting assignments for any of the spin variables. If one takes the randomized construction from Ref. Rendl et al. (2008), such inconsistencies will appear, even if the random vectors are taken to be the same for each block. On the contrary, since we essentially set each deterministic guess to be the sign of a one-body expectation values, our method ensures consistency.

To conclude, once a valid spin configuration has been extracted, we simply set the upper bound to be its corresponding energy . Surprisingly, we noticed that by following the above procedure, the exact ground states is usually recovered very soon in the branching (cfr. Figure 1 as an example). It then just takes time to find a matching lower bound to verify that this is indeed the lowest achievable energy. This makes us believe that our procedure is very efficient in finding the optimal deterministic configuration, even in cases where constrains on the computational resources might prevent one from achieving full convergence of the upper and lower bound.

#### Branching rule

Here we follow the same method outlined in Rendl et al. (2008), but with a different choice of branching procedure. Indeed, the authors in Rendl et al. (2008) choose the dichotomic choice to be the relative alignment of a pair of connected spins. That is, given a choice of indices , the two branches correspond to the two cases . However, as mentioned before, we prefer to branch on the value of the single spin, by choosing between the two values . The reason for this is that, in the latter case the number of possible branching steps depends only on the number of spins in the system. On the contrary, the former method involves an amount of branching choices that depends on the number of edges in the dependency graph , which can be much higher, often as high as , or even more.

Finally there is the question of which spin to choose for the next branching. The way we do this here, is based on the corresponding expectation value recovered from the moment matrices obtained by solving the SDP. The intuition is the following: spins with an expectation value close to zero are “difficult” choices, because flipping the value of such a spin is likely to lead to a slight change in the energy of the system; conversely, expectations values very close to are identified as “easy” choices, because flip such a spin is likely to lead to a significant change in the energy of the system. We set the branching rule to “easy-first”, that is, at the end of each optimization round, the next branching is performed on the most deterministic spin in the .

### iii.5 Possible improvements

There is some extent of freedom in the choices we outlined in the previous subsections. Given the huge difference in complexity that can be exhibited by various instances of the Ising model, we expect the optimal choice to be model dependent. Here we briefly discuss which modifications we imagine to be most useful for practical applications.

Let us start by recalling that, in order to accelerate the convergence process and to keep memory requirements low, it is crucial to reduce the initial lower-upper bound gap as much as possible and as early as possible. One way to do that is to modify the hybrid hierarchy level introduced above. In our applications, it was always enough to set the threshold to at most . However, such value can be significantly increased without affecting too much the scalability of the CBB. Indeed, the main bottleneck of our method is the memory required to solve the largest SDP. This depends mainly on the block leading to the largest matrix . Therefore, as long as increasing the level of the smaller blocks does not lead to bigger matrices that the one for the largest clique, the SDP will still have the same memory requirements – although the solving time will clearly increase.

Other branching rules can be also be adopted. For instance, one can replace the “easy-first” approach with a “difficult-first”. In this case, one picks the next branching from the least deterministic spin in the We expect the choice of the most effective branching rule to depend on the system under consideration.

Lastly, there are instances in which CBB does not outperform other methods. This is true for specific cases of very sparse problems, where linear-programming relaxations were shown to work very well Liers et al. (2004), or some hand-crafted models for which exact polynomial algorithms are known Mandrá et al. (2017). Nevertheless, it would still be interesting to see if one could combine the construction based on the chordal extension with those methods and provide some further advantage.

## Iv Conclusion

In this work we propose the chordal branch and bound method, which combines ideas form the state-of-the-art branch and bound method for finding exact ground states of classical spin model with such from convex optimization and allows. It allows to exploit the locality of physical interactions to solve Ising models more efficiently and up to larger system sizes with a systematic approach. The resulting chordal branch and bound (CBB) method yields a hierarchy of polynomial time upper and lower bounds on the ground state energy. If these bounds converge, it yields a certified and exact ground state energy, together with a ground state configuration. The favourable scaling in memory and runtime of our method enables the verification of quantum annealing devices, such as the D-Wave machine Lanting et al. (2014), degenerate optical parametric oscillators QNNCloud (2017), and variational circuits Farhi et al. (2014).

## V Acknowledgement

F.B acknowledges useful discussions with J. Tura and M. Navascués and support from the Spanish Ministry MINECO (a Severo Ochoa PhD fellowship).

C. G. acknowledges support by the European Union’s Marie Skłodowska-Curie Individual Fellowships (IF-EF) programme under GA: 700140 as well as financial support from the Spanish Ministry MINECO (National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No. SEV-2015-0522), FundaciÃ³ Cellex, Generalitat de Catalunya (Grants No. SGR 874 and No. 875, CERCA Programme, AGAUR Grant No. 2017 SGR 1341 and CERCA/Program), ERC (CoG QITBOX and AdG OSYRIS), EU FETPRO QUIC, EU STREP program EQuaM (FP7/2007-2017, Grant No. 323714), and the National Science Centre, Poland-Symfonia Grant No. 2016/20/W/ST4/00314

Parts of this work were completed while P.W. was employed in ICFO. He acknowledges financial support from the ERC (Consolidator Grant QITBOX), Spanish Ministry of Economy and Competitiveness (Severo Ochoa Programme for Centres of Excellence in R&D SEV-2015-0522 and QIBEQI FIS2016-80773-P), Generalitat de Catalunya (CERCA Programme and SGR 875), and Fundació Privada Cellex, computational resources granted by the High Performance Computing Center North (SNIC 2015/1-162 and SNIC 2016/1-320), and a hardware donation by Nvidia. This research was supported by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation.

A. A. acknowledges support from the ERC CoG QITBOX, the AXA Chair in Quantum Information Science, the Spanish MINECO (QIBEQI FIS2016-80773-P and Severo Ochoa SEV-2015-0522), the Generalitat de Catalunya (SGR1381 and CERCA Program), and the Fundació Privada Cellex.

### Footnotes

1. https://gitlab.com/FBaccari/Chordal_BB

### References

1. Daphne Koller, Nir Friedman, Lise Getoor,  and Ben Taskar, “Graphical models in a nutshell,” in Introduction to Statistical Relational Learning, edited by L. Getoor and B. Taskar (MIT Press, 2007).
2. Mahabaleswara R. Bhatt and Uday B. Desai, “Robust image restoration algorithm using markov random field model,” CVGIP: Graphical Models and Image Processing 56, 61–74 (1994).
3. Andrew Lucas, “Ising formulations of many NP problems,” Front. Phys. 2, 1–15 (2014)arXiv:1302.5843 .
4. Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik,  and Jeremy L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nat. Commun. 5, 4213 (2014)arXiv:1304.3061 .
5. Vasil S. Denchev, Sergio Boixo, Sergei V. Isakov, Nan Ding, Ryan Babbush, Vadim Smelyanskiy, John Martinis,  and Hartmut Neven, “What is the computational value of finite range tunneling?” Physical Review X 6, 031015 (2016)arXiv:1512.02206 .
6. Elizabeth Crosson and Aram W. Harrow, “Simulated Quantum Annealing Can Be Exponentially Faster Than Classical Simulated Annealing,” in 2016 IEEE 57th Annu. Symp. Found. Comput. Sci., Vol. 2016-Decem (IEEE, 2016) pp. 714–723, arXiv:1601.03030 .
7. Scott Kirkpatrick, C. Daniel Gelatt,  and Mario P. Vecchi, “Optimization by Simulated Annealing,” Science 220, 671–680 (1983).
8. Franz Rendl, Giovanni Rinaldi,  and Angelika Wiegele, “Solving max-cut to optimality by intersecting semidefinite and polyhedral relaxations,” Mathematical Programming 121, 307–335 (2008).
9. Hayato Waki, Sunyoung Kim, Masakazu Kojima,  and Masakazu Muramatsu, “Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity,” SIAM Journal on Optimization 17, 218–242 (2006).
10. Troels F. Rønnow, Zhihui Wang, Joshua Job, Sergio Boixo, Sergei V. Isakov, David Wecker, John M. Martinis, Daniel A. Lidar,  and Matthias Troyer, “Defining and detecting quantum speedup,” Science 345, 420–424 (2014)arXiv:1401.2910 .
11. Salvatore Mandrá, Helmut G. Katzgraber,  and Creighton Thomas, “The pitfalls of planar spin-glass benchmarks: raising the bar for quantum annealers (again),” Quantum Science and Technology 2, 038501 (2017).
12. Salvatore Mandrá and Helmut G Katzgraber, “A deceptive step towards quantum speedup detection,” Quantum Science and Technology  (2018).
13. Frauke Liers, Michael Jünger, Gerhard Reinelt,  and Giovanni Rinaldi, “Computing exact ground states of hard ising spin glass problems by branch-and-cut,” New optimization algorithms in physics 50, 6 (2004).
14. Peter Wittek, “Algorithm 950: Ncpol2sdpa—sparse semidefinite programming relaxations for polynomial optimization problems of noncommuting variables,” ACM Transactions on Mathematical Software 41, 21 (2015)arXiv:1308.6029 .
15. Mosek, “The MOSEK optimization software,”  (2010).
16. Paul I Bunyk, Emile M Hoskinson, Mark W Johnson, Elena Tolkacheva, Fabio Altomare, Andrew J Berkley, Richard Harris, Jeremy P Hilton, Trevor Lanting, Anthony J Przybysz, et al., “Architectural considerations in the design of a superconducting quantum annealing processor,” IEEE Transactions on Applied Superconductivity 24, 1–10 (2014).
17. Jean B Lasserre, “An explicit exact SDP relaxation for nonlinear 0-1 programs,” in International Conference on Integer Programming and Combinatorial Optimization (Springer, 2001) p. 293â303.
18. Jean B. Lasserre, “Global optimization with polynomials and the problem of moments,” SIAM Journal on Optimization 11, 796–817 (2001b).
19. Stefano Pironio, Miguel Navascués,  and Antonio Acín, ‘‘Convergent relaxations of polynomial optimization problems with noncommuting variables,” SIAM Journal on Optimization 20, 2157–2180 (2010).
20. Jean B. Lasserre, “Convergent SDP-relaxations in polynomial optimization with sparsity,” SIAM Journal on Optimization 17, 822–843 (2006).
21. Murat A. Erdogdu, Yash Deshpande,  and Andrea Montanari, “Inference in graphical models via semidefinite programming hierarchies,” Advances in Neural Information Processing Systems 30 (2017), arXiv:1709.06525 .
22. Martin Grötschel, Michael Jünger,  and Gerhard Reinelt, “Calculating exact ground states of spin glasses: A polyhedral approach,” in Heidelberg Colloquium on Glassy Dynamics, Vol. 275 (1987) pp. 325–353.
23. Francisco Barahona, “On the computational complexity of Ising spin glass models,” Journal of Physics A: Mathematical and General 15, 3241–3253 (1982).
24. Lieven Vandenberghe and Martin S. Andersen, “Chordal graphs and semidefinite optimization,” Foundations and Trends\textregistered in Optimization 1, 241–433 (2015).
25. Miguel Navascués, Stefano Pironio,  and Antonio Acín, “A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations,” New Journal of Physics 10, 073013 (2008)arXiv:0803.4290 .
26. Trevor Lanting, Anthony J. Przybysz, A. Yu. Smirnov, Frederico M. Spedalieri, Mohammad H. Amin, Andrew J. Berkley, Richard Harris, Fabio Altomare, Sergio Boixo, Paul Bunyk, Neil Dickson, C. Enderud, Jeremy P. Hilton, Emile Hoskinson, Mark W. Johnson, Eric Ladizinsky, N. Ladizinsky, R. Neufeld, T. Oh, Ilya Perminov, Christopher Rich, Murray C. Thom, Elena Tolkacheva, Sergey Uchaikin, A. B. Wilson,  and Geordie Rose, “Entanglement in a quantum annealing processor,” Physical Review X 4, 021041 (2014)arXiv:1401.3500 .
27. QNNCloud, “Quantum neural network: Optical neural networks operating at the quantum limit,”  (2017).
28. Edward Farhi, Jeffrey Goldstone,  and Sam Gutmann, “A quantum approximate optimization algorithm,”  (2014), arXiv:1411.4028 .
264113