A QUBO Formulation of Minimum Multicut Problem Instances in Trees for D-Wave Quantum Annealers

Quantum annealing algorithms were introduced to solve combinatorial optimization problems by taking advantage of quantum fluctuations to escape local minima in complex energy landscapes typical of NP − hard problems. In this work, we propose using quantum annealing for the theory of cuts, a field of paramount importance in theoretical computer science. We have proposed a method to formulate the Minimum Multicut Problem into the QUBO representation, and the technical difficulties faced when embedding and submitting a problem to the quantum annealer processor. We show two constructions of the quadratic unconstrained binary optimization functions for the Minimum Multicut Problem and we review several tradeoffs between the two mappings and provide numerical scaling analysis results from several classical approaches. Furthermore, we discuss some of the expected challenges and tradeoffs in the implementation of our mapping in the current generation of D-Wave machines.

Quantum Annealing algorithms (QA) constitute a paradigm of quantum computation focused on solving combinatorial optimization problems [1][2][3][4][5][6][7][8][9] . QA are based on quantum effects in order to escape local minima of a cost function by the effect of tunneling through barriers separating local minima. An example of physical realizations of quantum annealing are the D-Wave quantum computers which have proved to be advantageous to solve instances of some NP − hard problems 10 .
To solve a problem using the D-Wave architecture, we must express it as a quadratic unconstrained Boolean optimization problem (QUBO) or an equivalent Ising function defined on logical variables. Then, we embed the logical problem in the physical architecture of the quantum annealer by mapping logical variables and qubits. The final step consists of performing an annealing process and obtaining the results. There are some limitations of the D-Wave hardware such as the maximum number of physical qubits available (which is 2000 qubits in its most recent computer, the D-Wave 2000Q system); this limitation imposes a restriction on the size of the logical problem that can be embedded into the hardware. Another limitation is the connectivity of the current architecture which is in the form of a Chimera graph as, in order to represent a logical variable, it is necessary to use a chain of physical qubits. Finally, numerical precision is another key feature in the mathematical definition and computation of Ising model parameters due to the analog nature of the device and the presence of additive noise.
In this article, we present a QUBO transformation of a family of graphs, to be later used on a QA approach, in order to solve the Minimum Multicut Problem (MCC problem) which is of paramount importance in theoretical computer science and other disciplines such as computer vision. The MMC problem is defined as follows: given a graph G = (V, E) with associated positive weights to the edges and a list of vertex pairs (s i , t i ), 1 ≤ i ≤ k, a minimum weight set of edges separating each pair of vertices in the list must be found.
The MMC problem is an NP − hard problem. In this paper, we study a QUBO formulation of the MCC problem on a family of connected trees, being this variant also NP − hard. We show two QUBO constructions for the studied problem. The first one is based on an intuitive linear programming formulation used to obtain a function of degree equal to the largest unique path between each pair (s i , t i ), 1 ≤ i ≤ k in a tree. The second one is based on the basic idea that in order to separate each path between terminals s i and t i , it is necessary to remove at most k edges. This basic observation was used to define a QUBO function of degree at most k − 1. Also, a general penalty approach for the QUBO function based on Boolean circuits is shown.
−  The Adiabatic theorem asserts that for sufficiently large t a , for some solution ψ(t) to the Schrödinger equation with Hamiltonian H(t). Consequently, state ψ(t a ) will be very close to the ground state of Hamiltonian H p with high probability. A sufficient condition for the algorithm running time needed to satisfy the Adiabatic theorem is The degree of f, denoted as deg(f), is the cardinality of the largest subset S ⊆ {1, …, N} for which c(S) ≠ 0. Similarly, the size of f denoted as size(f) is the total number of variable occurrences in it, i.e., size(f) = Σ S:c(S)≠0 |S|.
Of our particular interest are the quadratic pseudo-Boolean functions f : {0, 1} ue N  → with deg(f ue ) ≤ 2 expressed by polynomials of the form: . Also, the quadratic pseudo-Boolean map coincides with the Ising model via variable substitution It is often the case that a problem can be formulated in terms of pseudo Boolean expressions of degree greater than two, which can subsequently be reduced to a QUBO function. However, this reduction commonly implies adding new variables to the former problem [12][13][14] (an introduction to NP−hard problems and quantum annealing-based algorithms can be found in 15 ).
In the following sections, we address the QUBO formulation of the MMC problem which is of utmost theoretical and practical importance in the computer science community.
Quantum formulation of the MMc problem. An undirected weighted graph has the form G For k = 1, the MMC problem reduces to the min-cut/max-flow problem that can be solved in polynomial time 16 . The problem is also tractable when k = 2, by using two applications of the min-cut/max-flow algorithm 17 . It becomes NP − hard when k ≥ 3 for general graphs, but can be solved in polynomial time for planar graphs for any fixed k 18 . It has been proved that the MMC problem is NP − hard and MAXSNP − hard on trees 18,19 . The MAXSNP − hardness of the MCC problem implies that no polynomial time approximation scheme exists unless P = NP.
The MMC problem can be applied in many areas such as telecommunication, routing, VLSI design and circuit partitioning 20,21 . In the following, we study the QUBO formulation of the MMC problem for the case when G is a tree.
A direct mapping of MMC into QUBO. Let G = (V, E, w) be a weighted graph and (s i , t i ), 1 ≤ i ≤ k a list of vertex pairs. For each edge e ∈ E, we introduce a Boolean variable x e such that x e = 0(1), if e is (not) considered for a MMC in G. If graph G is restricted to be a tree, then there exists a unique path in G for every pair of vertices. We define the unique path from s i to t i in G as p i where its length l i is equal to the number of edges that it crosses. The diameter of a tree T = (V, E), denoted as diam(T), is equal to the maximum path length between every pair of vertices in T.
Let = ∈ In (6), H w expresses the weight of the selected edges to be removed and H pen serves to add a penalty value when the considered edges do not correspond to a MMC. Based on this construction, the MMC problem on trees is equivalent to minimizing H over all possible assignments to the Boolean variables x e .
The penalty term satisfies 0 ≤ H pen ≤ kλ P , where λ P is a positive constant. In particular, when all paths are disconnected, H pen = 0 which means that at least one edge was removed in every path p i for i = 1, …, k. We must ensure that λ P is big enough so that the weight of an invalid multicut is greater than the weight of any valid multicut. The value of λ P can be upper bounded by λ = ∑ ∈ w e ( ) P e P . The degree of H pen is equal to the maximum length of the paths p i for i = 1, …, k. In other words, In order to optimize (6) via QA, we need to write it in QUBO form (possibly by using section 1.1). For instance, assume that each path p i does not have a trivial length, i.e. l i > 2; using the reduction method in 14 , the penalty term H pen can be transformed into a quadratic function adding a total of i k l 1 This transformation is as follows: • If l i is even then www.nature.com/scientificreports www.nature.com/scientificreports/ • If l i is odd then we find that is a vector of k 1 new variables, and k l 1 The proposed QUBO function obtained from (6) using the Ishikawa reduction method presented in 14 has as main advantage that its coefficients are small. This is a highly desirable property when programming the D-Wave quantum annealer because of its limited hardware precision to specify the values of Ising model parameters h i and J ij .
A disadvantage of function H given in (6) is that its degree depends on the length of the paths. In subsection 2.2 we present another method to build a low degree H pen , based on the number of edges shared between paths.
Construction of H pen based on crossing paths. Although several techniques exist for degree reduction of an arbitrary pseudo-Boolean function into a quadratic one, it is preferable to construct an initial expression of degree as low as possible. In this method, the key idea behind buliding a low degree function H pen is to notice that a multicut consists of a subset of edges of cardinality of at most k. This condition imposes a restriction on the number of edges that could be removed in order to disconnect each path. Based on this observation, our goal will be to construct a new penalty function H pen ′ such that it will penalize all subsets of edges of cardinality larger than k. For any i, j ∈ {1, …, k}, let ζ ij be a characteristic function defined as We say that two paths p i and is equal to the number of paths p i that intersect with path p j . Thus, the maximum number of paths than can intersect with a path p j is k − 1. It is assumed that the number of intersections cannot be greater than the length of a given path.
For any j ∈ {1, …, k} and η ∈ {1, …, k − 1}, let C ηj be given by where C ηj = 0 if the number of edges to be removed from path p j equals η, and C ηj > 0 otherwise. The new penalty term is ′ has as objective to penalize sets of edges of cardinality greater than k rather than penalizing sets of edges that do not correspond to a multicut as the term H pen does. H 0 pen ′ = if and only if every path p j is disconnected, and ′ > H 0 pen otherwise. The degree of H pen ′ is at most twice the maximum number of intersections in all paths, i.e.
pen j which does not depend on the length of paths p i in G.
The penalty function H pen ′ can be reduced into a QUBO function using the method in 12 , producing a quadratic expression, the size of which is polynomially bounded in size( ′ H pen ) and the number of new variables is O(n 2logd ) where d = deg( ′ H pen ) and n = |P|. A disadvantage of this method is that the resulting quadratic function has many large coefficients and also introduces many positive quadratic terms. These two effects make the minimization of the resulting function a hard problem 22,23 .

A Boolean circuit construction for H pen
Function φ cj can be constructed if we consider a Boolean circuit g: {0, 1} N → {0, 1} such that g(x) = 1 when φ cj (x) > 0, and g(x) = 0 otherwise. This Boolean circuit can always be expressed as a pseudo Boolean function φ g = c j φ ′ such that they both have the same values at every point. Let us formalize this result. A disjunctive form (DF) is an expression of the form for all k, l ∈ {1, …, m} with k ≠ l. An orthogonal DF ϕ represents a Boolean circuit g : {0, 1} N → {0, 1} if the truth value points of g coincide with the truth value points of ϕ. Furthermore, every Boolean circuit g : {0, 1} N → {0, 1} represented by an orthogonal DF has an associated multilinear polynomial introduced in 24 and presented in (18), The Boolean circuit g: {0, 1} N → {0, 1} can be constructed using the standard Karnaugh map procedure 25 . This technique receives as input a truth table on N Boolean variables and returns a Boolean circuit expressed as a sum of products (SOP) of the involved variables which can be implemented using AND, OR and NOT gates. After obtaining the Boolean formulae in terms of the elementary Boolean gates, the NOT gate is interpreted as For instance, the corresponding pseudo-Boolean function obtained through Karnaugh maps as SOP can be expressed as As it can be seen in (19), the degree of function φ ′ c j is equal to the length of path p j . This penalty function has the same degree as the expression given in (8), however, the former has the property that penalizes multicuts of cardinality greater than k.
Experiments and Discussion. We now present our algorithms and corresponding outcomes, as well as the limitations of our proposal. First, let us review the road map (steps) for solving an optimization problem using a quantum annealing approach: (i) Select an optimization combinatorial problem (in our case, the MMC problem).
(ii) Construct a pseudo-Boolean function on binary variables for the selected optimization problem, so that those assignments that minimize the expression also correspond to solutions of the given problem. Since the D-Wave quantum processor has a specific architecture, it only supports pairwise interaction between qubits. Hence, the constructed pseudo-Boolean function must be of degree two, i.e., a QUBO expression. It is not always possible to directly obtain a QUBO expression for a given optimization problem, so in practice we have a high degree pseudo-Boolean function which is to be transformed into a QUBO expression at a later stage. (iii) Once we have computed a QUBO expression for the selected problem, the next step is to embed the logical problem into the fixed architecture of the quantum processor. This architecture is represented as a Chimera graph which has a limited interconnectivity (see Fig. 1 for an example of a Chimera graph). At this step, a key feature to remember is physical resources vs logical variables as, in order to embed a logical problem with arbitrary interconnectivity, we will frequently need a larger number of physical qubits than the number of logical variables; moreover, if the number of logical variables scales up rapidly, we may run out of physical qubits. The estima tion of the minim um number of physical qubits required to embed a logical problem is an active area of research 26,27 . The Chimera graph architecture implements an Ising model, hence the embedding process maps logical variables into Ising variables. Another important consideration is that the D-Wave processor has limited precision to represent Ising coefficients; consequently, the c oefficients of the Q UBO expression should not be too large so that they can be correctly mapped into the quantum pro cessor and results can be discriminated from noise level. (iv) Finally, the annealing proc ess is initiated to find the minimum energy configuration in which the solutio n of the problem is codified. www.nature.com/scientificreports www.nature.com/scientificreports/ For steps (i) and (ii), we consider the QUBO formulations of the MMC pr oblem present ed in section 2: the formulatio n with direct reduction, the one base d on crossing of paths and the one based on logical circu its. For the last two cases, the reduction methods cited in Section 1.1 are use d to obtain the QUBO functions.
In the following section, we consider the mapping problem presented above in step (iii), namely mapping logical variables in a QUBO function to physical qubits in the actual hardware architecture.
Minor embedding. The D-Wave processor can be represented by an architecture known as a Chimera graph which consists of an M × N-lattice of blo cks each one having 2L physical qubits for a total of 2MNL qubits (see Fig. 1). Each block in the Chimera graph is a L-bipartite graph and each physical qubit is connected with at most six other qubits. To solve a problem using a D-Wave processor, it is necessary to represent an Ising/QUBO problem as a subgraph of the Chimera. However, it is seldom possible to find a one-to-one mapp ing of logical variables with physical qubits.
The method to find an equivalent subgraph into the Chimera to a given logical problem expressed as an Ising/ QUBO function is called minor embeddi ng and it is stated as follows:   www.nature.com/scientificreports www.nature.com/scientificreports/ Finding a minor embedding which uses the minimum number of qubits is a hard problem in general 26 . However, there are heuristic algorithms that obtain approximated minor embeddings in polynomial time. One of these techniques is the one proposed in 28 which maps a logical variable to a chain of qubits. The main disadvantage of these techniques is that they create long chains of qubits to allow the connectivity of the logical problem into the Chimera.
In this paper, we use the D-Wave solver application programming interface (SAPI) which provides heuristic algorithms for minor embedding. Heuristic algorithms are constrained by a peculiar property: their time complexity (i.e., their worst and average runtime) is usually lower than exhaustive search algorithms at the cost of producing approximative solutions most of the time and optimal results only rarely. In our case, SAPI heuristic algorithms produce different embeddings that are not always optimal with respect to the number of physical qubits. For each or our experiments, we have run the heuristic embedding algorithms twenty times and have  www.nature.com/scientificreports www.nature.com/scientificreports/ chosen only those embeddings that have shorter qubit chains. This process is carried out on a conventional digital computer, before the adiabatic process.
Instance problem generation. For this paper, we have generated random instances of the MMC problem, namely random trees, according to the following procedure: (1) generate a connected random graph 16 , (2) compute their corresponding spanning tree (ST) and (3) assign integer positive weights to the edges of the obtained ST. Here, we consider the generation of connected random graphs using the Erdös-Rényi (ER) and the Watts-Strogatz (WS) models. In the ER model, a random graph G(n, p) with n vertices is constructed by including each possible edge with a probability p independently from every other edge. In the WS model, a random graph G(n, r, β) with n vertices is constructed by creating an initial ring of vertices with each vertex connected to its r nearest neighbors. Thus, replace every possible edge {u, v} with a new edge {u, v′} with probability β, duplicated edges are forbidden, but original edges may end up being reinstated. Figure 2 shows some examples of trees generated using the ER and WS models. We chose the ER and WS models because they allow us to generate trees with small and large diamete rs, respectively. For instance, Fig. 2 (left) shows a tree with a diameter of 4 using the ER model and Fig. 2 (right) shows a tree with a diameter of 10 using the WS model, in both cases the number of vertices in the random graph are the same.
Let us consider an example of how the QUBO formulations given in Section 2 for the MMC problem are constructed. Figure 3(a) shows a random ER graph for 50 vertices and probability p = 0.3, and Fig. 3(b) shows its corresponding random spanning tree. Assume an instance of the MMC problem with unitary weights for k = 3 and set of pairs {(s 1 , t 1 ), (s 2 , t 2 ), (s 3 , t 3 )}; from Fig. 3(b), it can be seen that paths p j , j = 1, 2, 3 contain edges {1, 3, 4}, {4, 5, 6} and {2, 6, 7}, respectively. Notice that the length of paths p j is l j = 3 for j = 1, 2, 3, and paths p 1 and p 2 share the edge 4, and paths p 2 and p 3 share the edge 6.
The direct mapping given in (6) for the instance of the MMC problem in Fig. 3(b) is The logical graph obtained from the QUBO reduction of (20) can be seen in Fig. 3(c) where the blue color vertices are the original problem variables and the brown color vertices are the introduced new variables. Notice that the induced subgraphs from vertices {1, 3, 4} and {4, 5, 6} share vertex 4, and the induced subgraphs from vertices {4, 5, 6} and {2, 6, 7} share vertex 6. Also, notice that vertices {1, 3, 4} together with the new variable 8 form a complete graph; the same is also true for vertices {4, 5, 6} and {2, 6, 7} together for new variables 9 and 10, respectively. From the above, it can be concluded that if the paths do not intersect each other, then their corresponding logical graphs will be disconnected.
On the other hand, the mapping based on crossing paths given in (14) for the instance of the MMC problem in Fig. 3(b) is www.nature.com/scientificreports www.nature.com/scientificreports/ 91 22  22  22  43  22  43  22  14  14  14  14  14  14  14  14 14 which is already in QUBO form since we penalize cuts of cardinality at most three. In other words, each path p j , j = 1, 2, 3 will be disconnected by removing at most one edge. Figure 3(d) shows the logical graph of expression (21) where no new variables are used. Notice that only the induced subgraphs of vertices {1, 3, 4}, {4, 5, 6} and {2, 6, 7} are shown. Finally, the Karnaugh mapping given in (19) for the instance of the MMC problem in Fig. 3(b) is whose logical graph coincides with the QUBO reduction of (20).
Problem scaling. Figure 4 shows a comparison of the number of quadratic terms generated by the proposed QUBO formulations presented in sections 2.1, 2.2 and 2.3, that will be named as the direct, intersection and Karnaugh methods, respectively. We generated random trees using the ER and WS models with probabilities p = 0.3 and β = 0.12, for n = 50 and number of pairs k = 3, 4, 5. We generate random WS trees with paths between the pairs (s i , t i ) of length less than 8. For each random tree model used, the corresponding QUBO function was constructed using the direct, intersection or Karnaugh methods. The energy function given in (6) was used for the three methods with different penalty terms. The direct method corresponds to the penalty term as presented in (8), the intersection method has as penalty term the expression shown in (14), and the Karnaugh method has as penalty term the expression given in (19). We use the reduction method of Ishikawa 14 explained in section 2.1 as well as Freedman's method 13 , to obtain the corresponding QUBO function. Figure 4(a) presents is the average number of quadratic terms over 100 ER random trees, generated using the direct, intersection and Karnaugh methods for k = 3, 4, 5. Figure 4(b) shows their corresponding average number of variables after the degree reduction of the constructed QUBO functions shown in Fig. 4(a). The vertical lines in Fig. 4(b,d) indicate a standard deviation. As it can be seen in the Fig. 4(a,b), the intersection method uses the minimum number of variables to represent the MMC problem and the Karnaugh method requires more variables to represent the same problem. It is important to mention that the average diameter of the generated ER random trees was 4. We choose n = 50 since the number of quadratic terms is an increasing function with respect to n. Similarly, www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 4(c,d) show the average number of quadratic terms using WS random trees for the direct, intersection and Karnaugh methods. It is remarkable that Karnaugh's method requires a huge number of quadratic terms and variables in comparison with the direct and intersection methods. On the other hand, the direct method uses the minimum number of quadratic terms. In this case, the average diameter of the generated WS random trees was 10. Figure 5 presents an even more descriptive comparison of the number of quadratic terms and number of variables obtained after degree reduction, for the direct, intersection and Karnaugh methods. The area in each box represents the interquartile range, the statistical median as a horizontal line in the box, and its vertical lines the lower and upper whiskers. The data shown in Fig. 5(a-d) correspond to the same data shown in Fig. 4. As can be noted from Fig. 5, the lower and upper whiskers for the Karnaugh method are distant from their median, which explains the standard deviation in Fig. 4(b,d).
In view of these scaling results, it can be seen that Karnaugh's construction given in (19) has terms of degree equal to the length of the paths between pairs (s i , t i ). Thus, when applying a reduction method on (19), the number of quadratic terms increases, as opposed to the method given in (14) whose degree depends on the number of crossing paths. As a consequence, the dimension of the search space of its corresponding optimization problems also increases. www.nature.com/scientificreports www.nature.com/scientificreports/ Classical solution performance. In order to test the performance of obtaining solutions using the proposed formulations, the qbsolv tool was used to solve large QUBO problems by partitioning into subproblems targeted for execution on a D-Wave system 29 . The qbsolv tool returns approximated solutions of large QUBO problems. The experiments were performed on a Desktop computer MacBook Air with a Intel Core i5 processor at 1.3 GHz and 4 GB of RAM. We generated 80 random trees for n = 50 and k = 3 using the ER and WS models, and for each random tree, their corresponding QUBO functions were constructed using the direct, intersection and Karnaugh methods.
Before mapping a QUBO expression to Ising, constant terms are omitted without changing the original problem because constant terms cannot be represented in the graph topology. For the D-Wave SAPI software, an auto-scaling mode is used which utilizes the largest possible parameter ranges. In the case of the qbsolv tool we do not have precision limitations on the parameters, so each QUBO instance is directly submitted. Figure 6(a,b) show the energy solutions obtained using qbsolv for ER and WS random trees, respectively. To compare the approximated energies solutions, they were sorted in an increasing order over the same set of instances for the direct, intersection and Karnaugh methods. Figure 6(a,b) show that the direct and intersection methods have similar energy solutions for the ER and WS random trees. On the other hand, the Karnaugh method has larger energy solutions. The latter can be explained since the large number of used variables in the Karnaugh method, as can be seen in Fig. 4. We were unable to calculate the optimal energy solutions of the random instances since the large number of variables. The cases for k = 4, 5 were omitted since they have similar behavior.
The D-Wave SAPI also provides classical algorithms to solve QUBO/Ising problems by using simulated quantum annealing (SQA) 3 on a Chimera graph of dimension 4 × 4 blocks. Figure 6(c,d) show the energy solutions obtained using D-Wave SAPI for ER and WS random trees, respectively. In this case, we use the same 80 random trees for n = 50 and k = 3 as in Fig. 6(a,b). For each instance of the MMC problem, 1000 readouts was requested to the SQA algorithm. Figure 6(c) shows the minimum energies among the requested 1000 readout per instance, for ER random trees. As can be seen in Fig. 6(c), the obtained energies have a similar behavior as in Fig. 6(a). Figure 6(a) also shows the energies for the case of WS random trees using the SQA algorithm. In this latter case, it was not possible to obtain the energies for the Karnaugh method since the embedding algorithm fails to embed large QUBO functions in a Chimera graph of a 4 × 4 dimension. Figure 7 shows how the sorted energies in Fig. 6 are correlated with respect to the number of quadratic terms in the set of instances. Figure 7(a,b) present the distribution of quadratic terms for the direct, intersection and Karnaugh methods using random ER and WS trees in Fig. 6. Also, Fig. 7(c-f) show the number of quadratic terms of the direct, intersection and Karnaugh methods, where the instances were arranged based on the same energy sorting shown in Fig. 6(a-d), respectively. For instance, the higher energies in Fig. 6(a) are correlated with the higher number of quadratic terms generated by the Karnaugh method. The same correlation is observed for the classical solvers qbsolv and D-Wave simulator.
We have proposed QUBO formulations of the MMC problem restricted to the family of trees. Our simulation results show that the direct and intersection methods have a similar scaling, in the number of quadratic terms, and performance. On the other hand, the Karnaugh method has the worst results, mainly by the required large number of variables and quadratic terms in its construction.

Discussion
We have proposed a method to formulate the Minimum Multicut Problem into the QUBO representation, and the technical difficulties faced when embedding and submitting a problem to the quantum annealer processor. We have considered a special NP − hard case of the Minimum Multicut problem based on random trees representations. This special case allows us to formulate a QUBO expression that can be embedded into the Chimera graph using a moderate number of qubits after the degree reduction of the high degree expression. Moreover, we have proposed the Karnaugh method to analytically construct instances of the MMC problem and the preliminary results of our algorithms are promising.