Bayesian network structure learning using quantum annealing

,


Introduction
Bayesian networks are a widely used probabilistic graphical model in machine learning [1].A Bayesian network's structure encapsulates conditional independence within a set of random variables, and, equivalently, enables a concise, factored representation of their joint probability distribution.There are two broad classes of computational problems associated with Bayesian networks: inference problems, in which the goal is to calculate a probability distribution or the mode thereof given a Bayesian network and the state of some subset of the variables; and learning problems, in which the goal is to find the Bayesian network most likely to have produced a given set of data.Here, we focus on the latter, specifically the problem of Bayesian network structure learning.Bayesian network structure learning has been applied in fields as diverse as the short-term prediction of solar-flares [2] and the discovery of gene regulatory networks [3,4].The problem of learning the most likely structure to have produced a given data set, with reasonable formal assumptions to be enumerated later, is known to be NP-complete [5], so its solution in practice requires the use of heuristics.
Quantum annealing is one such heuristic.Though efficient quantum algorithms for certain problems are exponentially faster than their classical counter-part, it is believed that quantum computers cannot efficiently solve NP-complete problems [6].However, there exist quantum algorithms (for problems not known or believed to be NP-complete) that have a provable quadratic speedup over classical ones [7,8].There is therefore reason to believe quantummechanical effects such as tunneling could provide a polynomial speedup over classical computation for some sets of problems.The recent availability of quantum annealing devices from D-Wave Systems has sparked interest in the experimental determination of whether or not the current generation of the device provides such speedup [9,10].While there exists prior work related to "quantum Bayesian networks" [11] and the "quantum computerization" of classical Bayesian network methods [12], the results presented here are unrelated.
In this paper, we describe how to efficiently map a certain formulation of Bayesian Network Structure Learning (BNSL) to Quadratic Unconstrained Binary Optimization (QUBO).The QUBO formalism is useful because it is mathematically equivalent to that of a set Ising spins with arbitrary 2-body interactions, which can be mapped to the Ising spins with a limited 2-body interaction graph as implementable by physical quantum annealing devices.Similar mappings have been developed and implemented for lattice protein folding [13,14], planning and scheduling [15], fault diagnosis [16], graph isomorphism [17], training a binary classifier [18,19], and the computation of Ramsey numbers [20].
To map BNSL to QUBO, we first encode all digraphs using a set of Boolean variables, each of which indicates the presence or absence of an arc (i.e.directed edge), and define a pseudo-Boolean function on those variables that yields the score of the digraph encoded therein so long as it satisfies the necessary constraints.This function is not necessarily quadratic, and so we apply standard methods to quadratize (i.e.reduce the degree to two) using ancillary variables.We then introduce ancillary variables and add additional terms to the pseudo-Boolean function corresponding to constraints, each of which is zero when the corresponding constraint is satisfied and positive when it is not.The resulting QUBO instance is defined over O(n 2 ) Boolean variables when mapped from a BNSL instance with n Bayesian network variables.Interestingly, the structure of the QUBO is instance-independent for a fixed BNSL size.Because embedding the structure of QUBO into physical hardware is generally computationally difficult, this is an especially appealing feature of the mapping.
We also show sufficient lower bounds on penalty weights used to scale the terms in the Hamiltonian that penalize invalid states, like those containing a directed cycle or with parent sets larger than allowed.In a physical device, setting the penalty weights too high is counterproductive because there is a fixed maximum energy scale.The stronger the penalty weights, the more the logical energy spectrum is compressed, which is problematic for two reasons: first, the minimum gap, with which the running time of the algorithm scales inversely, is proportionally compressed, and, second, the inherently limited precision of a physical device's implementation of the interaction strengths prevents sufficient resolution of logical states close in energy as the spectrum is compressed.
The utility of the mapping from BNSL to QUBO introduced here is not limited to quantum annealing.Indeed, the methods used here were motivated by a previous mapping of the same problem to weighted MAX-SAT [21].Existing simulated annealing code is highly optimized [22] and may be applied to QUBO instances derived from our mapping.In that case, there is no need to quadratize, because simulated annealing does not have the limitation to 2-body interactions that physical devices do.With respect to penalty weights, while simulated annealing does not have the same gap and precision issues present in quantum annealing, there may still be reason to avoid setting the penalty weights too high.Be-cause the bits corresponding to arcs with different i.e. terminal vertices do not interact directly, many valid states are separated by invalid ones, and so penalty weights that are too strong may erect barriers that tend to produce basins of local optima.While simulated annealing directly on digraph structures is possible, mapping to QUBO and performing simulated annealing in that form has the advantage that it enables the exploitation of existing, highly optimized code, as well as providing an alternative topology of the solution space and energy landscape.
BNSL has a special property that makes it especially well-suited for the application of heuristics such as QA: Unlike in other problems where anything but the global minimum is undesirable or those in which an approximate solution is sufficient, in BNSL there is utility in having a set of high scoring DAGs.The scoring function encodes the posterior probability, and so sub-but near-optimal solution may be almost as probable as the global optimum.In practice, because quantum annealing is an inherently stochastic procedure, it is run many times for the same instance, producing a set of low-energy states.In cases where the BN structure is learned for the purpose of doing inference on it, a high-scoring subset of many quantum annealing runs can utilized by performing Bayesian model averaging, in which inference is done on the set of likely BNs and the results averaged proportionally.
In Section 2, we review the formalism of Bayesian networks and BNSL (2.1) and quantum annealing (2.2), elucidating the features that make the latter suitable for finding solutions of the former.In Section 3, we develop an efficient and instance-independent mapping from BNSL to QUBO.In Section 4, we provide sufficient lower bounds on the penalty weights in the aforementioned mapping.In Section 5, we discuss useful features of the mapping and conclude.In the Appendix, we prove the sufficiency of the lower bounds given; the methods used to do so may be useful in mappings for other problems.

Bayesian Network Structure Learning
A Bayesian network (BN) is a probabilistic graphical model for a set of random variables that encodes their joint probability distribution in a more compact way and with fewer parameters than would be required otherwise by taking into account conditional independences among the variables.It consists of both a directed acyclic graph (DAG) whose vertices correspond to the random variables and an associated set of conditional probabilities for each vertex.Here and throughout the literature, the same notation is used for both a random variable and its corresponding vertex, and the referent will be clear from context.
Formally, a BN B for n random variables X = (X i ) n i=1 is a pair (B S , B P ), where B S is a DAG representing the structure of the network and B P is the set of conditional probabilities {p(X i |Π i (B S ))|1 ≤ i ≤ n} that give the probability distribution for the state of a variable X i conditioned on the joint state of its parent set Π i (B S ) (those variables for which there are arcs in the structure B S from the corresponding vertices to that corresponding to X i ; we will write simply Π i where the structure is clear from context).Let r i denote the number of states of the variable X i and q i = j∈Πi r j denote the number of joint states of the parent set Π i of X i (in B S ).Lowercase variables indicate realizations of the corresponding random variable; x ik indicates the k-th state of variable X i and π ij indicates the j-th joint state of the parent set Π i .The set of conditional probabilities B P consists of n probability distributions (θ ij ) qi j=1 n i=1 , where θ ij = (θ ijk ) ri k=1 is the conditional probability distribution for the states (x ik ) ri k=1 of the variable X i given the joint state π ij of its parents Π i (i.e.p(x ik |π ij ) = θ ijk ).
Given a database D = {x i |1 ≤ i ≤ N } consisting of N cases, where each x i denotes the state of all variables X, the goal is to find the structure that maximizes the posterior distribution p(B S |D) out of all possible structures.By Bayes's Theorem, The marginal probability of the database p(D) is the same for all structures, so assuming that each structure is equally likely, this simplifies to In Section 3.5, we describe how to account for certain types of non-uniform prior distributions over the graph structures.With certain further reasonable assumptions, namely multinomial sampling, parameter independence and modularity, and Dirichlet priors, the latter conditional probability is where N ijk is the number of cases in D such that variable X i is in its k-th state and its parent set Π i is in its j-th state, N ij = ri k=1 N ijk , α ijk is the hyperparameter for θ ijk in the Dirichlet distribution from which θ ij is assumed to be drawn, and α ij = ri k=1 α ijk [23].
Given a database D, our goal is equivalent to that of finding the structure with the largest likelihood, i.e. the structure that yields the largest probability of the given database conditioned on that structure.We do this by encoding all structures into a set of bits and defining a quadratic pseudo-Boolean function on those bits and additional ancillary bits whose minimizing bitstring encodes the structure with the largest posterior probability.

Quantum Annealing
Quantum annealing is a method for finding the minimum value of a given objective function.It is the quantum analogue of classical simulated annealing, where the computation is driven by quantum, rather than thermal, fluctuations [24].A similar procedure is called adiabatic quantum computation, in which the adiabatic interpolation of a Hamiltonian whose ground state is easily prepared to one whose ground state encodes the solution to the desired optimization problem guarantees that final state is indeed the ground state of the latter [25].The formalism for both is similar, and the methods described here are useful for both.Specifically, the time-dependent Hamiltonian is for 0 ≤ t ≤ T , where H 0 is the initial Hamiltonian, H 1 is the final Hamiltonian, A(t) is a real monotonic function such that A(0) = 1 and A(T ) = 0, and B(t) is a real monotonic function such that B(0) = 0 and B(T ) = 1.The adiabatic theorem states that if the system starts in the ground state of H 0 and H(t) varies slowly enough, then the system will be in the ground state of H 1 at time T .Using this procedure to solve an optimization problem entails the construction of H 1 such that its ground state encodes the optimal solution.In practice, arbitrary Hamiltonians are difficult to implement, but this is ameliorated by results showing the ability to effectively implement arbitrary Hamiltonians using physically-realizable connectivity through various gadgetry with reasonable overhead [26,27].
The main contribution of this paper is a construction of H 1 such that its ground state encodes the solution for a given instance of BNSL.Specifically, we construct an instance of QUBO whose solution is the score-maximizing DAG; there is a simple transformation between a classically defined QUBO instance and a diagonal quantum 2-local Hamiltonian consisting of only Pauli Z and ZZ terms [28].
When the desired Hamiltonian is diagonal and 2local an embedding technique called graph-minor embedding can be used [29,30].A graph G is a minor of another graph H if there exists a mapping from vertices of G to disjoint, individually connected subgraphs of H such that for every edge e in G there is an edge in H whose adjacent vertices are mapped to by the adjacent vertices of the edge e.The desired Hamiltonian and hardware are considered as graphs, called the logical and physical respectively, where qubits correspond to vertices and edges correspond to a 2-body interaction, desired or available.Graph-minor embedding consists of two parts: finding a mapping of the logical vertices to sets of physical as described, and setting the parameters of the physical Hamiltonian such that the logical fields are distributed among the appropriate physical qubits and there is a strong ferromagnetic coupling between physical qubits mapped to my the same logical qubit so that they act as one.Determining the graph-minor mapping, or even if the logical graph is a minor of the physical one, is itself NP-hard, and so in practice heuristics are used [31].

Mapping BNSL to QUBO
We use n(n − 1) bits d = (d ij ) 1≤i<j≤n i =j to encode each of the possible arcs in a directed graph, where d ij = 1 indicates the presence of the arc from vertex X i to vertex X j and d ij = 0 indicates its absence.In this way, the matrix whose entries are {d ij } is the adjacency matrix of a directed graph (where d ii = 0).Let G(d) be that directed graph encoded in some d.The mapping consists of the construction of a function of these "arc bits" that is equal to the logarithm of the score of the structure they encode, as well as a function that penalizes states that encode graphs with directed cycles.Additionally, due to resource constraints, we add a function that penalizes structures in which any node has more than m parents and allow that the scoring function only works on states that encode structures in which each vertex has at most m parents.

Score Hamiltonian
For numerical efficiency, it is the logarithm of the likelihood for a given structure that is actually computed in practice.The likelihood given in Equation 3 decomposes into a product of likelihoods for each vari-able, which we exploit here.Let i.e. the negation of the "local" score function, and The negation is included because while we wish to maximize the likelihood, in QUBO the objective function is minimized.We wish to define a quadratic pseudo-Boolean function and define Any pseudo-Boolean such as score has a unique multinomial form and s i (Π i (G(d))) depends only on arcs whose head is X i (i.e.those encoded in d i ), so we write without loss of generality From this it is clear that w i (∅) = s i (∅).If X i has a single parent X j , then the above simplifies to which yields w i ({j}) = s i ({X j }) − s i (∅) for arbitrary j.Similarly, if X i has two parents X j and X k , then Extrapolating this pattern, we find that Note that the general form given in Equation 9includes terms of order (n − 1).Ultimately, we require a quadratic function and reducing high-order terms to quadratic requires many extra variables.Therefore, we limit the number of parents that each variable has to m via H max , described below, and allow that the score Hamiltonian actually gives the score only for structures with maximum in-degree m: which is equal to

Max Hamiltonian
Now we define a function max whose value is zero if variable X i has at most m parents and positive otherwise.This is done via a slack variable y i for each node.Define i.e. d i is the in-degree of x i , i.e. µ is the number of bits needed to represent an integer in [0, m], i.e. y i ∈ Z is encoded using the µ bits y i = (y il ) µ l=1 ∈ B µ , and where δ max > 0 is the weight of the penalty.For convenience, we also write max takes its minimal value of zero when However, when d i > m, because y i ≥ 0, we cannot set y i in that way.By taking the derivative with respect to y i , (18) we see that max takes its minimum value over the domain of y i when y i = 0, and that value is Noting that max is added.

Acyclicity
Lastly, we must ensure that the structure encoded in {d ij } has no directed cycles.We do so by introducing additional Boolean variables r = (r ij ) 1≤i<j≤n that will encode a binary relation on the set of variables.Every directed acyclic graph admits at least one topological order of the vertices, and a graph with a directed cycle admits none.A topological order "≤" of the vertices {X i } of a digraph is a total order thereon such that for every edge (i, j) in the digraph X i ≤ X j .Such an order is not unique in general.Let r ij = 1 represent x i ≤ x j and r ij = 0 represent x i ≥ x j .
To ensure acyclicity, we define a function H trans (r) such that H trans (r) is zero if the relation encoded in {r ij } is transitive and is positive otherwise, as well as a function H consist such that H consist (d is zero if the order encoded in {r ij } is consistent with the directed edge structure encoded by {d ij } and positive otherwise.First, we ensure that {r ij } is transitive.Because if a tournament has any cycle, it has a cycle of length three, it is sufficient to penalize directed 3-cycles.Define where and δ (ijk) trans is the penalty weight added if r encodes either 3-cycle containing {x i , x j , x k }.Note that the superscripted indices on the penalty weight variable are unordered so that δ trans for all permutations (i , j , k ) of (i, j, k).
Second, we must penalize any state that represents an order and a directed graph that are inconsistent with each other, i.e. in which r ij = 1 and (x j , x i ) ∈ E(G(d)) or r ij = 0 and ((x i , x j ) ∈ E(G(d)).Equivalently, we want to ensure that neither and which has the desired features.Again the superscripted indices on the penalty weight variable are unordered, so that δ which takes on its minimal value of zero if G(d) is a DAG and is strictly positive otherwise.

Total Hamiltonian
Putting together the parts of the Hamiltonian defined above, define In the next section, we show lower bounds on the penalty weights therein that ensure that the ground state of the total Hamiltonian H encodes the highestscoring DAG with a maximum parent set size of m.The sets of variables described above have the following sizes: A quadratic pseudo-Boolean function can be identified with a graph whose vertices correspond to its arguments and whose edges correspond to nonzero quadratic terms.The graph identified with the Hamiltonian described above for m = 2 has several features that indicate it may be difficult to embed in sparsely connected physical devices.First, for each variable X i there is a clique consisting of the variables {d ji } ∪ {y il }, whose order is (n − 1) + µ.Second, the set of variables {r ij } are almost fully connected.

Utilizing Prior Information
The mapping so far described assumes a uniform prior distribution over all possible DAGs of the appropriate size and with the given maximum number of parents.However, there are situations in which it may be desirable to fix the presence or absence of an arc in the search space.This could be because of domain knowledge or because hardware limitations prevent the implementation of the mapping for all arcs, in which case resort can be made to iterative search procedures such as the bootstrap method [32].To realize the reduction in qubits needed by accounting for such a reduced search space, suppose that we wish to only consider network structures that include the arc (i, j), where i < j without loss of generality.We then set d ij = 1, d ji = 0, and r ij = 1.Similarly, if (i, j) is to be excluded, we set d ij = d ji = 0 and keep r ij as a free variable.This can be done for any number of arcs.The Hamiltonian remains unchanged once these substitutions are made, and the lower bounds on the penalty weight remain sufficient, with the exception of the terms used in quadratization in the case m > 2, in which case the quadratization should be done after substitution to utilize the reduction in degree of some terms.

Penalty Weights
In the expression above, there are several sets of free parameters called penalty weights: ≤ n, i = j}, and {δ ijk trans |1 ≤ i < j < k}.They are associated with penalty terms, i.e. parts of the Hamiltonian whose value is zero on states satisfying the corresponding constraint and is positive on states violating it.The purpose of their inclusion is to ensure that the energy-minimizing state of the total Hamiltonian satisfies the requisite constraints by increasing the energy of those that do not.More strongly, the penalty weights must be set such  Each vertex corresponds to a bit in the original QUBO and an edge between two vertices indicates a non-zero quadratic term containing the corresponding bits.The central cluster is the order bits used to enforce acyclicity; it is highly connected but not complete.Each "spike" corresponds to a variable X i in the Bayesian network.
The outer two vertices are the corresponding slack bits {y il } and the remaining inner vertices are the arc bits {d ji } representing those arcs for which the corresponding Bayesian network variable is the head.Each spike is a clique, due to H max (independent of which the arc bits for a given BN variable are fully connected due to H score ).Each arc bit is connected to a single order bit and each order bit is connected to two arc bits, due to H consist .Bottom: Schematic of the Hamiltonian.The row of disks represents all of the bits in the original QUBO problem, colored consistently with the logical graph above.They are grouped into three sets: the arc bits representing the presence of the possible arcs, the order bits representing a total ordering by which we enforce acyclicity, and the slack bits used to limit the size of the parent sets.An arrow from a group of bits to a part of the Hamiltonian indicates that that part of the Hamiltonian is a function of that set of bits.
that the ground state of the total Hamiltonian is the lowest energy state of H score that satisfies the constraints.Here we provide sufficient lower bounds on the penalty weight necessary to ensure that this purpose is met.No claim is made to their necessity, and tighter lower bounds may exist.It is important to note that these bounds are mathematical, i.e. they ensure their purpose is met as stated above.In pure adiabatic quantum computation, in which the quantum system is in its ground state for the duration of the algorithm, this is sufficient (though the computation time necessary for the conditions of the adiabatic theorem to hold may be longer than otherwise if a penalized state has lower energy than the first excited unpenalized state).In practical quantum annealing, however, a combination of physical effects may cause the optimal value (in the sense of minimizing the energy of the lowest-energy state found, which may or may not be the global ground state) of the penalty weights to be less than these bounds.This remains the case even for bounds shown to be tight with respect to their mathematical properties.The bound presented for each of the three sets of penalty weights is based on the notion that only the addition of an arc (i.e.changing some d ij from 0 to 1) can lead to the violation of two of the constrains we are concerned with: the maximum number of parents and the consistency of the arc bits and the order bits.Therefore, we can use a basis for how strongly the associated penalty needs to be the greatest difference in the energy of H score adding each arc can contribute.The penalty for the third constraint, the absence of directed 3-cycles among the order bits, will then be a function of the penalty for the consistency constraint.
Formally, we wish to set the penalties such that for any d violating at least one of the constraints, we have min where This is achieved by showing that for any such d violating at least one constraint, there is another d that satisfies all the constraints such that min Because d satisfies all the constraints, min y,r which implies the inequality in (28).In this section, we state the bounds and provide brief justification, but relegate the proofs to Appendix B.

Auxiliary Quantity
In this section, we briefly define an auxiliary quantity, , (32) that will allow us to define the maximum penalty weights associated with the bounds described previously.For details of the calculation of this quantity, see Appendix A. In general it is possible that ∆ ji < 0 as defined above.We thus define the quantity In the proof of the bounds, the two following facts will be useful. Claim

"Maximum" Penalty Weights
Here we show a lower bound for {δ This idea can be applied iteratively to show that if, for some d and i, |d i | > m, there is some d with lesser energy such that d j = d j for j = i and |d i | ≤ m.This idea in turn can be applied iteratively to show that if for some

"Reduction" Penalty Weights
The degree of the "score" Hamiltonian H score is natively m-local as constructed.If m = 2, as it often will be in practice, the total Hamiltonian is natively quadratic.If m > 2, additional ancilla bits are needed to reduce the locality.The general method for doing this is to replace the conjunction of a pair bits with an ancilla bit and to add a penalty term with sufficiently strong weighting that penalizes states in which the ancillary bit is not equal to the conjunction to which it should be.For m = 3, this can be done using n (n − 2) 2 /4 ancilla bits, but no more, where each H (i) score containing n − 1 arc bits is quadratized independently; furthermore, heuristic methods have been developed that reduce needed weight of the penalty terms [33].For m = 4, at most n n−1 2 ancilla bits are needed.More generally, O(n 2 log d ) is ancilla bits are needed [34].Because the proof of the bounds on the other penalty weights are secular as to the degree of H score so long as {∆ ij } is computed appropriately, the quadratization of H score , including the addition of penalty terms and the needed weights therefor, can be done using the standard methods described in the literature independent of the other penalties described here.

"Cycle" Penalty Weights
First, we that if the consistency penalty is set high enough, for any d encoding a graph with a 2-cycle, there is some d encoding one whose minimal value of H over y, r is strictly less than that of d.

Claim 5 (Removal of 2-cycles.). If δ (ij)
consist > max{∆ ij , ∆ ji } for all 1 ≤ i < j ≤ n , then for all d such that G(d) contains a 2-cycle, there is some d ≤ d such that G(d ) does not contain a 2-cycle and min y,r H(d, y, r) > min y,r H(d , y, r).
Second, we show that for any d that encodes a digraph without a 2-cycle, the minimal value of H consist over all r is zero.

Claim 6 (Sufficiency of "Consistency" Penalty Weights
trans for 1 ≤ i < j ≤ n, then for all d such that G(d) contains no 2-cycle, H consist (d, r * ) = 0, where r * = arg min r H cycle (d, r).
Third, we show that for any d that encodes a digraph not containing a 2-cycle but that is not a DAG, there is some d that does encode a DAG and whose minimal value of H over all y, r is strictly less than that of d.

Claim 7 (Sufficiency of "Transitivity" Penalty Weights
trans for 1 ≤ i < j ≤ n and δ does not contain a 2-cycle but does contain a directed cycle there is some d such that G(d ) is a DAG and min y,r H(d, y, r) > min y,r H(d , y, r).
Lastly, we show that for all d that encode a digraph that is not a DAG, there is some d that does encode a DAG and whose minimal value of H over all y, r is strictly less than that of d.Claim 8 (Sufficiency of "Cycle" Penalty Weights).If δ trans for all 1 ≤ i < j ≤ n and δ and min y,r H(d , y, r) < min y,r H(d, y, r).

Overall Sufficiency
Finally, we show that the digraph encoded in the ground state of the total Hamiltonian H is a DAG and has a maximum parent set size of at most m, and that it is the solution to the corresponding BNSL instance.
Claim 9 (Overall Sufficiency).If δ for all 1 ≤ i < j ≤ n and δ The strict inequalities used in the specification of the lower bounds ensures that the global ground state is a score-maximizing DAG with maximum parent set size m, but replacing them with weak inequalities is sufficient to ensure that the ground state energy is the greatest score over all DAGs with maximum parent set size m.However, the latter is of little interest in the present situation because it is the DAG itself that is of interest, not its score per se.

Conclusion
We have introduced a mapping from the native formulation of BNSL to QUBO that enables the solution of the former using novel methods.
The mapping is unique amongst known mappings of optimization problems to QUBO in that the logical structure is instance-independent for a given problem size.This enables the expenditure of considerably more computational resources on the problem of embedding the logical structure into a physical device because such an embedding need only be done once and reused for new instances.The problem addressed, BNSL, is special among optimization problems in that approximate solutions thereto often have value rivaling that of the exact solution.This property, along with the general intractability of exact solution, implies the great value of efficient heuristics such as SA or QA implemented using this mapping.
At present, only problems of up to seven BN variables can be embedded in existing quantum annealing hardware (i.e. the D-Wave Two chip installed at NASA Ames Research Center), whereas classical methods are able to deal with many of tens of BN variables.Nevertheless, the quantum state of the art is quickly advancing, and it is conceivable that quantum annealing could be competitively applied to BNSL in the near future.Given the already advanced state of classical simulated annealing code, it is similarly conceivable that its application to the QUBO form described here could be competitive with other classical methods for solving BNSL.
where the first term is independent of d ji and thus cancels in the argument of the minimization on the right-hand side of Equation 32.Thus ∆ ji simplifies to For m = 1, ∆ ji is trivially −w i ({j}), the constant value of the expression to be extremized in Equation 36 regardless of the values of {d ki |k = i, j}.For m = 2, ∆ ji can still be calculated exactly: However, for m ≥ 3, calculation of the extremum in Equation 36 is an intractable optimization problem in its own right and therefore we must resort to a reasonable bound.Because ∆ ji will be used in finding a lower bound on the necessary penalty weights, caution ditates that we use, if needed, a greater value than necessary.A reasonable upper bound on the true value is: B Proofs of Penalty Weight Lower Bounds Proof.In the statement of the claim, we implicitly assume that δ consist > 0 for all 1 ≤ i < j ≤ n.Let r * = arg min r H cycle (d, r).Because d ji ≥ d ji for all i, j such thati = j and 1 ≤ i, j ≤ n, and because max > max j =i ∆ ji for all 1 ≤ i ≤ n, then for all d such that, for some i * , and min y,r H(d, y, r) > min y,r H(d , y, r).
Proof.We prove the existence of such a d by construction.Let , where j * = arg min j∈{j|d ji * =1} ∆ ji * .
First, we note that by design min yi H which rearranges to In context, By Claim 1 and the fact that d ≤ d, Proof.We prove the sufficiency of the given bound by iterative application of Claim 3. Let d (0,0) ≡ d.
Consider an arbitrary d (l) .If G(d (l) ) does not contain a directed 2-cycle, then l = l * and so we set d = d (l * ) to complete the proof.Otherwise, choose some 2-cycle in G(d (l) ) arbitrarily, i.e. some pair {i, j} such that (i, j), (j, i) ∈ E(G(d)), or, equivalently, that d ij = d ji = 1.Without loss of generality, assume i < j.Let r * ≡ arg min r H cycle (d (l) , r) and i.e. the arc in G(d (l) inconsistent with G(r * ).
Proof.We prove the claim via its contrapositive: for all d, r, if H consist (d, r) > 0, there is some r such that H cycle (d, r) > H cycle (d, r ), so r = arg min r H cycle (d, r).
Consider an arbitrary d and some r such that H consist (d, r) > 0. The positivity of H consist (d, r) indicates that there is at least one inconsistency between d and r, i.e. there is some (i * , j * ) such that For convenience, we prove the claim for the case in which i * < j * ; the proof provided can be easily modified for the case in which i * > j * .Let r be the same as r exept in the bit corresponding to this inconsistency: . Then Furthermore, Together, the above imply Claim 7 (Sufficiency of "Transitivity" Penalty Weights).If δ trans for 1 ≤ i < j ≤ n and δ (ijk) trans = δ trans > max 1≤i ,j ≤n i =j ∆ i j for 1 ≤ i < j < k ≤ n, then for all d such that G(d) does not contain a 2-cycle but does contain a directed cycle there is some d such that G(d ) is a DAG and min y,r H(d, y, r) > min y,r H(d , y, r).
Proof.Consider an arbitrary d (l) such that G(d (l) ) does not contain a 2-cycle but does contain a directed cycle.Let r (l) ≡ arg min r H cycle (d (l) , r).By Claim 6, H consist (d (l) , r (l) ) = and so H cycle (d (l) , r (l) ) = H trans (d (l) , r (l) ).
If δ (ijk) trans = δ trans for 1 ≤ i < j < k ≤ n, i.e. the trasitivity penalty weight is uniform for all directed triangles, then H trans (r) is equal to the product of δ trans and the number of directed triangles in the tournament G(r) for all r.
In any tournament with a positive number of directed triangles, there is always some arc whose switch of direction lowers the number of directed triangles.Let (i * , j * ) be such such an arc for r (l) .Define Define some r such that r ij = r ij , (i, j) = (i * , j * ), 1 − r ij , (i, j) = (i * , j * ), , which would have the properties that H consist (d (l) , r ) = 0 and, by construction, H trans (r ) < H trans (r (l) , so that H cycle (d (l) , r ) < H cycle (d (l) , r (l) ) = min r H cycle (d (l) , r).Because G(d (l) ) does not contain a 2-cycle, d Because of this, it must be that H trans (r (l) ) ≥ H trans (r (l+1) ), whose negation contradicts the definition of r (l+1) .Therefore, ) Furthermore, while H max and H cycle are natively 2local, H score is m-local.For each variable x i there are n−1 l possible parent sets of size l and the same number of corresponding l-local terms in H score .If m = 3, the full set of n−1 3 high-local terms j∈J d ji ||J| = 3 corresponding to parent sets of the variable x i can be reduced using (n−2) 2 4 ancilla variables.In total, n (n−2) 2 4 ancilla variables are needed to reduce H score to 2-local.

Figure 1 :
Figure1: Top: Logical Graph for n = 7 BN variables with a maximum of m = 2 parents.Each vertex corresponds to a bit in the original QUBO and an edge between two vertices indicates a non-zero quadratic term containing the corresponding bits.The central cluster is the order bits used to enforce acyclicity; it is highly connected but not complete.Each "spike" corresponds to a variable X i in the Bayesian network.The outer two vertices are the corresponding slack bits {y il } and the remaining inner vertices are the arc bits {d ji } representing those arcs for which the corresponding Bayesian network variable is the head.Each spike is a clique, due to H max (independent of which the arc bits for a given BN variable are fully connected due to H score ).Each arc bit is connected to a single order bit and each order bit is connected to two arc bits, due to H consist .Bottom: Schematic of the Hamiltonian.The row of disks represents all of the bits in the original QUBO problem, colored consistently with the logical graph above.They are grouped into three sets: the arc bits representing the presence of the possible arcs, the order bits representing a total ordering by which we enforce acyclicity, and the slack bits used to limit the size of the parent sets.An arrow from a group of bits to a part of the Hamiltonian indicates that that part of the Hamiltonian is a function of that set of bits.
} that guarantees that if d is such that max 1≤i≤n |d i | > m there exists a d with lesser total energy such that max 1≤i≤n |d i | ≤ m.To do so, we show that if, for some d and i, |d i | > m, there is a d such that |d i | = |d i | − 1 and min y,r H(d, y, r) ≥ min y,r H(d , y, r).
and min y,r H(d, y, r) > min y,r H(d , y, r).Claim 4 (Sufficiency of "Maximum" Penalty Weight).If δ (i) max > max j =i ∆ ji for all i, then for all d such that max i |d i | > m, there is a d ≤ d such that max i |d i | ≤ m and min y,r H(d, y, r) > min y,r H(d , y, r).
a DAG, and max i |d * i | ≤ m, where d * = arg min d {min y,r H(d, y, r)}.

=Claim 4 (
H score (d) + min y H max (d, y) + min r H cycle (d, r) > H score (d ) + min y H max (d , y) + min r H cycle (d , r) Sufficiency of "Maximum" Penalty Weight).If δ (i) max > max j =i ∆ ji for all i, then for all d such that max i |d i | > m, there is a d ≤ d such that max i |d i | ≤ m and min y,r H(d, y, r) > min y,r H(d , y, r).
, and min y,r H(d , y, r) < min y,r H(d, y, r).Setting d ≡ d (n,0) completes the proof.Claim 5 (Removal of 2-cycles.).If δ (ij) consist > max{∆ ij , ∆ ji } for all 1 ≤ i < j ≤ n ,then for all d such that G(d) contains a 2-cycle, there is some d ≤ d such that G(d ) does not contain a 2-cycle and min y,r H(d, y, r) > min y,r H(d , y, r).