Optimized QUBO formulation methods for quantum computing
Abstract
Quantum computers have strict requirements for the problems that they can efficiently solve. One of the principal limiting factor for the performances of noisy intermediate-scale quantum (NISQ) devices is the number of qubits required by the running algorithm. Several combinatorial optimization problems can be solved with NISQ devices once that a corresponding quadratic unconstrained binary optimization (QUBO) form is derived. Several techniques have been proposed to achieve such reformulations and, depending on the method chosen, the number of binary variables required, and therefore of qubits, can vary considerably. The aim of this work is to drastically reduce the variables needed for these QUBO reformulations in order to unlock the possibility to efficiently obtain optimal solutions for a class of optimization problems with NISQ devices. This goal is achieved by introducing novel tools that allow an efficient use of slack variables, even for problems with non-linear constraints, without the need to approximate the starting problem. We divide our new techniques in two independent parts, called the iterative quadratic polynomial and the master-satellite methods. Hence, we show how to apply our techniques in case of an NP-hard optimization problem inspired by a real-world financial scenario called Max-Profit Balance Settlement (MPBS). We follow by submitting several instances of this problem to two D-wave quantum annealers, comparing the performances of our novel approach with the standard methods used in these scenarios. Moreover, this study allows to appreciate several performance differences between the D-wave Advantage and next-generation Advantage2 quantum annealers. We show that the adoption of our techniques in this context allows to obtain QUBO formulations with significantly fewer slack variables, i.e., around 90% less, and D-wave annealers provide considerably higher correct solution rates, which moreover do not decrease with the input size as fast as when adopting standard techniques.
I Introduction
In recent years, many efforts have been devoted to understand to what extent modern noisy intermediate-scale quantum (NISQ) devices can help to solve optimization problems of any sort (see Refs. [1, 2, 3] for recent reviews). Particular attention has been paid to their potential to solve NP-hard and NP-complete problems, which most of the times can be formulated as quadratic unconstrained binary optimizations (QUBO). For instance, all famous NP Karp’s problems [4] can be casted in this form [5].
The possibility to solve QUBO problems with quantum devices has been studied on several supports, especially on quantum annealers [2, 3, 6, 1, 7, 8] but also on logical-gate quantum computers [9, 10, 11, 12] and Rydberg atom arrays [13, 14, 15]. Although these technologies are not expected to solve NP problems efficiently [16, 17], namely providing exponential speed-ups compared to classical strategies, possible polynomial in-time advantages attract great interest. Indeed, given the importance of the real-world applications involved with the solutions of these hard combinatorial problems [2], any achievable improvement is well-received. Exemplary scenarios come from finance [18, 19, 20, 21, 22, 23], logistics [24, 25, 26] and drug discovery [27].
Consider the situation where we aim to solve a generic constrained optimization problem with a QUBO solver. The first step is to reformulate it in a quadratic unconstrained form, namely as a QUBO. We call logical variables those defining the initial problem. The constraints attached to the original problem can be enforced in a QUBO form by employing additional slack variables. The larger is the total number of variables, logical and slack, in a QUBO problem, the harder is in general to solve it.
Whereas the well-established procedures to translate optimization problems into QUBOs [28] can be efficient in several scenarios, their implementation can be extremely inefficient at times. To be more precise, we say that a method to obtain these formulations is inefficient for a problem whenever it requires too many slack variables, namely their number is comparable with that of logical variables. In these cases, the practical usefulness of QUBO solvers can be highly limited. For instance, if we have a maximum size for the problems that our solver can receive, any reduction of slack variables allows to increase the total number of logical variables associated to the initial problem, namely larger problems can be tackled. Moreover, the more slack variables we implement, the higher is the connectivity that we require in our QUBO and therefore in the QUBO solver that we aim to use. This consequence can be a very limiting factor for the performances of modern NISQ devices.
The main motivation of this work is to unlock the possibility to solve certain classes of optimization problems with those quantum technologies specialized to solve QUBO problems. For this purpose, we introduce the iterative quadratic polynomial and master-satellite methods. Their goal is to provide QUBO forms requiring a minimal employment of slack variables. The former consists in a new paradigm to translate problem constraints into corresponding quadratic penalty forms, while the latter allows the simultaneous enforcement of different constraints sharing the same restricted set of binary variables. Noteworthy, our techniques can also be used to enforce non-linear equality and inequality constraints defined through integer or non-integer parameters, without the need to approximate the constraints or to employ slack variables solely to obtain corresponding linearizations. Hence, the same computational effort is required both for linear and non-linear constraints, where no approximate enforcing of the constraints is required. Moreover, the QUBO variables employed by our methods for each constraint are not necessarily completely connected. These features allow improving the output quality of NISQ devices used as QUBO solvers as fewer and less-connected qubits are in general required.
We apply these techniques in the context of an NP-hard problem coming from finance, namely the Max-Profit Balance Settlement (MPBS) problem [29]. We provide a drastic reduction of the slack variables necessary to enforce the corresponding inequality constraints. We generate several instances of this problem reflecting some main features of realistic datasets [30] and we show that our methods provide the corresponding QUBO forms by employing around 90% less slack variables, if compared with the QUBO forms obtainable with a standard procedure [28].
We solve several instances of MPBS with two different quantum annealers manufactured by D-Wave Systems, Inc. [31], namely the D-wave Advantage_system4.1 and Advantage2_prototype1.1. The annealer Advantage2_prototype1.1, being a prototype of the next-generation D-wave annealers, has a reduced number of qubits but is characterized by higher connectivity and lower noise. Hence, not only do we compare the outputs obtained when the QUBO problems are generated with different approaches, but we also provide a comparison between the performances of these annealers. We show that, while by using the standard approach the solution quality drops quite fast with the input size, our methods unlock the possibility to consider much larger instances. The number of successes that we obtain when the QUBO forms are generated with our methods are always significantly higher: the multiplicative factor between the successes obtainable with the two methods ranges from (smallest instances with Advantage_system4.1) to (largest instances with Advantage2_prototype1.1).
We point out the existence of other techniques attempting to reduce the slack variables employed for QUBO reformulations. For instance, the alternating direction method of multipliers [32] and the augmented Lagrangian method [33, 34] have recently been adopted in this context [35, 36, 37], where preliminary QUBO problems have to be solved in order to get the final QUBO formulation of the original optimization problem. Another approach consists in enforcing constraints in an approximate fashion via a finely-tuned rescaling of the parameters defining the target linear inequalities [38, 39, 40]. Differently from these strategies, neither we require to solve preliminary QUBOs nor we approximate the starting optimization problem. Rather, we change the paradigm under which constraints are enforced.
The content of this work is structured as follows. In Section II, we introduce the fundamental tools needed throughout this work, as the definition of graph-based problems, QUBO forms and the standard technique to achieve these quadratic forms. The first part of Section III is devoted to describe in an intuitive way the main ideas behind of our novel methods. Instead, the iterative quadratic polynomial and the master-satellite methods are respectively described in details in Sections III.1 and III.2. The latter method requires a tuning of the QUBO penalty multipliers that is radically different from those found in standard QUBO formulations. Different strategies for this tuning are presented in Section III.3. The last heuristic that we introduce is a generalization of the master-satellite approach, suitable for those cases where three or more constraints are defined over the same variables, which is proposed in Section III.4. We identify a class of optimization problems, namely the locally-constrained problems, in Section IV, which are particularly suitable for the application of our techniques. The definition of the NP-hard Max-Profit Balance Settlement (MPBS) problem is given in Section V, where we also analyse its QUBO formulation and the corresponding needs in terms of slack variables when using standard techniques. Instead, we show how to apply our novel techniques to MPBS step-by-step in Section VI. The corresponding improved performance obtainable in the context of quantum annealing is shown in Section VII, where the novel and the standard strategies are put in comparison.
II Main tools
QUBO problems are defined as optimizations of quadratic forms defined over a set of binary variables , without the possibility to attach any constraint [28]. Hence, we have the following equivalent notations:
(1) |
where is an upper-triangular real-valued matrix and the corresponding optimization can either be a minimization or a maximization. Although these two cases are computationally equivalent, we consider the maximization case from now on.
Quantum annealers are the QUBO solvers that we adopt later in this work to solve these optimization problems, which in general are NP-hard [5]. The explanation of how quantum annealers work in order to obtain optimal solutions of QUBO problems can be found, e.g., in Refs. [1, 2, 3]. We point out that we do not consider other heuristics from classical computation to solve these problems. An example of such techniques is called simulated annealing [41, 42, 43], the classical counterpart of quantum annealing, which searches into the space of possible solutions, while slowly decreasing the probability to accept worse alternatives of the current solution. A comparison between quantum and simulated annealing can be found in Ref. [44]. A different heuristic designed for the solution of QUBO problems is called tabu search, which is put in contrast with simulated annealing in Ref. [45].
We proceed by explaining how several optimizations can be translated into a QUBO problem and where does the necessity to introduce new methods for this translation come from. Consider an optimization problem defined through the maximization of an objective function, which from now on we call , where some constraints on the possible solutions must be satisfied. Usually, from the very definition of the objective function, we obtain a natural set of logical variables associated to the units that describe the problem itself. For instance, if the problem requires to find a particular sub-graph of an input graph, the binary variables may be associated to its arcs and/or nodes, or, if we aim to subdivide a particular quantity among different parties, the binary variables may be used to describe the corresponding fractions. Whenever can be written as a linear or quadratic form of the variables , its QUBO form is obtained straightforwardly. Nonetheless, in order to translate the constraints attached to our optimization problem in a QUBO form, we usually need additional binary variables , called slack variables, which are employed solely to enforce constraints [28].
This work is devoted to expose novel techniques to implement constraints in QUBOs. We point out that, even if some constraints could be particularly demanding in terms of slack variables, others may be not. Hence, whenever we use the expression “target constraints”, we refer to those on which we decide to apply our novel tools.
We study the generic situation of problems defined on directed, weighted and node-attributed multigraphs, namely:
Definition 1 (R-multigraph).
An R-multigraph consists of the quadruple , where is a set of nodes, is a multiset of ordered pairs of nodes called arcs, is a function assigning positive real numbers, called weights, to arcs and is a function assigning real-valued -dimensional vectors, called attributes, to nodes.
Hence, we allow multiple arcs between two nodes , where each arc has an origin , a target node and a positive weight . Moreover, each node may have attributes given by a real-valued -dimensional vector . Naturally, what follows can be applied to simpler graphs, e.g., undirected and/or unweighted.
The majority of the optimization problems considered in the literature can be casted as follows: find the sub-graph of a given R-multigraph that maximizes the quantity under some given constraints, where is a generic sub-graph of . In order to write these problems in a more compact form, we introduce a binary variable for each arc . Hence, we obtain a bijection , where and . Each sub-graph is identified by a single combination of , where if the corresponding arc is included in and otherwise. Given this bijection, each combination of defines a different and , where if there exists at least one arc in such that is either the origin or the target node. For instance, .
We start by considering the class of optimization problems that have objective functions and constraints that are, respectively, (at most) quadratic and linear in the binary variables . In case of equality constraints and inequality constraints , an optimization problem of this kind can be written as:
(2) |
Moreover, we assume the parameters defining the linear inequality constraints to be integer (positive or negative), so that can only assume integer values. The linearity assumption of all the constraints and the requirement that the inequality ones are integer valued are necessary in order to introduce the standard method to obtain QUBO reformulations, as we show in Section II.1. Later, in Section III, we introduce novel methods that do not require the starting optimization problem to satisfy none of these assumptions: non-linear and non-integer constraints can be enforced without additional effort with respect to the linear integer-valued constraint cases.
II.1 Standard QUBO formulation
The goal of a QUBO reformulation of the optimization problem (2) is to provide a quadratic unconstrained binary optimization problem
(3) |
that provides the same optimal configurations provided by Eq. (2), where are additional slack variables that are employed to enforce inequality constraints in a quadratic form. For the sake of clarity, we dropped the notation that makes explicit the dependence of , and from and . As we see later, linear equality constraints do not need any slack variables to be enforced in this form. In order for this reformulation to have the meaning that we required, we have to find quadratic forms that provide penalties, namely are smaller or equal than -1, whenever the corresponding equality or inequality constraint are violated by , otherwise, if satisfied, they have to be zero valued for at least one value of . For instance, if and only if violates the first equality constraint and if and only if satisfies the same constraint. Here, is a multiplier that has to be chosen large enough in order to make sure that all the combinations of violating one or more constraints cannot be the optimal configurations selected by Eq. (3).
The standard approach to enforce constraints in a quadratic form [28] reads as follows:
(4) |
where are linear functions of the binary slack variables . The quadratic term that enforces the equality constraint in the QUBO form is simply given by . The reason of this choice is straightforward: any that violates provides a negative contribution. Therefore, as this optimization is formulated as a maximization, those leading to a violation of an equality constraint are discouraged. Hence, if is set large enough, we are sure that the solution of Eq. (4) satisfies all the equality constraints.
For what concerns the inequality constraints, we build some that help to translate these constraints into equality ones [28]. Before explaining how this goal is achieved, we notice that, if the inequality constraints are satisfied for at least one combination of , namely if there exists a solution of Eq. (2), then . Hence, if assumes all the values in the interval , checking whether satisfies the inequality constraint is equivalent to check whether satisfies the equality constraint for at least one value of . Hence, the enforcement of this new equality constraint in the QUBO form can be made similarly to by considering the quadratic term . Finally, in order for to have the properties described above, we use the definition:
(5) |
where and are such that assumes all the integer values in the interval . Hence, is the number of slack variables needed by this method to enforce in a QUBO form, namely Eq. (4). Notice that in Eq. (5) we implied a relabelling of the slack variable indices such that those associated to are the first , where . It is possible to apply similar techniques also in the case of non-integer-valued inequality constraints, but in general their enforcement requires more slack variables.
In summary, each term of the first and second sum of Eq. (4) is respectively associated to an equality or inequality constraint which provides a negative contribution, or penalty, if and only if violates the corresponding condition. More precisely, if violates a constraint, the contribution of the corresponding constraint term is strictly negative for all . On the other side, the maximizations over the same terms are zero-valued if and only if satisfies the corresponding conditions. Hence, a large-enough value of guarantees that the optimal solution does not violate any constraint and therefore Eqs. (2) and (4) have the same solution. Whereas a straightforward choice of this parameter is , we provide other techniques later in this work for this tuning.
We underline that, in order to achieve this QUBO formulation, we have to make a disjoint use of slack variables, namely each slack variable appearing in , does not appear in any other for . Moreover, the total number of slack variables used to enforce solely depends on . As a result, the number of slack variables deployed is:
(6) |
III Iterative quadratic polynomial and master-satellite methods
The main goal of this work is to introduce a novel set of tools to obtain QUBO formulations of constrained binary optimization problems with the least possible amount of slack variables. This can help to reduce the number of qubits required to solve these problems on NISQ devices, potentially unlocking the possibility to tackle larger and more complex optimizations. Our strategy consists of two partially independent methods: the iterative quadratic polynomial (IQP) method and the master-satellite (MS) method. While here we provide an intuitive explanation of these techniques, a more detailed description can be found in the next sections.
The IQP method provides a new way to enforce constraints in a QUBO form. This method is particularly well-suited for constraints defined on few binary variables. Moreover, it can be easily generalized to enforce non-linear equality and inequality constraints into a QUBO form defined with integer or non-integer parameters. We summarize the main idea behind this method in Figure 1, while its technical explanation is presented in Section III.1.
The second main contribution of this work is a novel interplay between the penalty quadratic forms of a QUBO problem that enforce constraints, called the MS method. While the standard approach enforces constraints independently from each other, our method exploits a synergy between different constraints sharing the same binary variables. In case of two or more such constraints, we can divide them into master and satellite. Master constraints are enforced for each possible binary variable combination, while satellite constraints are enforced only for those combinations that already satisfy the master constraints. As a result, the penalty quadratic forms associated with master constraints provide penalties/no-penalties whenever the original constraints are violated/satisfied, namely they perfectly enforce their corresponding constraints. Instead, the penalty quadratic forms associated with satellite constraints are enforced only for those binary variables combinations that already satisfy the master constraints. A detailed description of the MS method is presented in Section III.2.
A consequence of this approach is that the quadratic forms associated to satellite constraints may provide accidental incentives, namely a positive contribution, for those binary variables combinations violating the master constraints. Hence, in order to avoid the promotion of those combinations, we amplify the influence of master constraints penalty terms until we have a net penalty even for those combinations receiving accidental incentives. This is done by a careful calibration of the QUBO penalty multipliers, namely some quantities analogous to the parameter from Eq. (4). This tuning is explained in detail in Section III.3. We show that, through the adoption of this strategy, we are able to reduce the employment of slack variables. We summarize the idea behind the MS method in Figure 2.
It is evident that, given a particular optimization problem, the task of deciding which constraints have to be enforced as master and which as satellite is not a trivial one. Remember that, this method requires the presence of at least one master constraint and the more master constraints we have, the more efficient the implementation of satellite constraints is going to be. We can either decide to consider as master those constraints that require the least amount of slack variables, with or without the IQP method, or maybe those that are the most regular among the graph, e.g., do not depend on and but solely on the arcs and nodes in the multigraph. This situation is examined through an exemplary problem in Section V. We point out that, the separation of master and satellite constraints described here is called bipartite, namely when constraint can either be master or satellite. We introduce a concatenated version of the MS method in Section III.4, which results to be even more efficient in terms of slack variables employed.
The implementation of the MS method is particularly straightforward when satellite constraints are enforced with the IQP method. Indeed, this method gives us the possibility to decide whether to associate a penalty or not to precise binary variables combinations. Moreover, when satellite constraints are enforced with the IQP method, even less slack variables are necessary, if compared to the enforcement of the same constraints as master. Indeed, the slack variables needed by the IQP method increase with the number of binary variables combinations for which we require a penalty or not. Hence, we anticipate that the benefits of the conjunction between the IQP and the MS (IQPMS) methods are maximized when dealing with constraints defined on few common binary variables. Nonetheless, in Appendix A we show how to apply the MS method when constraints are enforced with the standard method explained in Section II.1.
In summary, while the IQP method is more efficient for constraints defined over few variables, the MS approach can be applied for sets of constraints defined on the same subsets of binary variables. Depending on the features of the combinatorial optimization problem studied, either one or both methods can be applied. For the sake of clarity, we introduce the class of locally-constrained optimizations problems in Section IV, for which a simultaneous adoption of the IQPMS methods can be efficiently implemented, namely all constraints can be enforced with these new methods. Moreover, the IQPMS methods perform the best when the graphs defining these problems are lowly connected.
Finally, we point out that, while both the IQP and the MS methods offer several advantages, it is important to note their potential limitations. The IQP method can become computationally complex for constraints defined over several variables, and the effectiveness of the MS method may depend on the particular selection of master and satellite constraints. Nonetheless, compared to existing approaches in the literature, such as the alternating direction method of multipliers [32], the augmented Lagrangian method [33, 34] or approximate QUBO enforcing, which have recently been adopted in Refs. [35, 36, 37, 38, 39, 40], our methods do not require solving preliminary QUBOs or approximating the problem constraints.
III.1 IQP method
In this section we introduce a novel technique to enforce constraints in the QUBO form. We start by focusing on the construction of a quadratic polynomial that enforces an inequality constraint of the target optimization problem, see Eq. (2). As we explain later, this technique can also be used for more general constraints, but for the sake of simplicity we begin with the linear inequality case. We consider the case where the constraint is defined over binary variables, namely over a subset of all the binary variables defining the starting optimization problem. Hence, differently from the previous section (see Figure 1), we relabel the binary variables associated to the target inequality constraint by using .
Consider a generic quadratic polynomial (QP) in the following variables , where are the binary variables associated to the studied constraint and are slack variables:
(7) |
The number of free parameters in are:
(8) |
The slack variables used to define are disjoint from those considered to define QPs associated to other constraints 111If and are enforced by and , respectively, they would make use of and independent slack variables..
QP with no-slack variables: We start by considering the scenario without slack variables, namely . In order to obtain a quadratic form that enforces the constraint , we have to solve the following linear system in the parameters :
(9) |
If such a solution exists, then we managed to translate the constraint into a QUBO form without using slack variables. We say that each value of identifies a different sub-constraint, which can be either assigned to a no-penalty or to a penalty . Hence, in order for to represent , it has to satisfy sub-constraints, one for each value of . It may sound that we can use no slack variables only when is not smaller than the number of sub-constraint . Nonetheless, as we see later, it is not uncommon to find quadratic polynomials that enforce a number of sub-constraints larger than the number of its free parameters. This is because inequality sub-constraints, namely , are less restrictive than equality ones. Finally, notice that, if the constraint is satisfied for all , then we do not even need to enforce this constraint via a penalty polynomial: it is not violated for any combination of and therefore no penalty has to be applied. Indeed, in this case we would simply obtain the null polynomial, namely for all and .
QP with slack variables: The implementation of slack variables into this scenario cannot be made by simply considering the linear system (9), where is replaced by . Indeed, even if has several more degrees of freedom coming from the presence of slack variables, these are not exploited in the system as there is no reference of them in the linear system that we should solve. Consider for instance a value of associated to a no-penalty sub-constraint. In this case, we would require for all the values of . Hence, even if we have more degrees of freedom compared to the no-slack variable method, in this way we would solely have more conditions in the linear system. Therefore, it is easy to convince ourselves that this approach provides no improvement with respect to the scenario of a QP with no slack variables.
In order to exploit the potential of the slack variables inside a generic QP in this context, we start by noting the following. If we replace with in Eq. (9), the penalty (inequality) sub-constraints are less restrictive than the no-penalty (equality) ones. Indeed, the former type allows to assume any value smaller or equal to -1 and for different values of and the same the polynomial can assume different penalties. Instead, for no-penalty sub-constraints we require to be exactly 0, for all and . Nonetheless, if is associated to a no-penalty, we do not need to be true for all the values of , it is enough to require it for at least one , while at the same time no incentives are provided for the other values of . Indeed, is one component of the QUBO formulation of an optimization problem, which undergoes a maximization over (which contains all the variables in ) and . Hence, our minimal requirements on the parameters , and defining are:
(10) |
or, equivalently,
(11) |
The enforcement of the sub-constraints proposed in Eqs. (10) and (11) can be casted as a feasibility problem. In case of a single slack variable , namely for , it assumes the form:
(12) |
where the symbols and represent, respectively, the AND and OR logical operations. Hence, the solution of this feasibility problem provides the parameters and that defines the quadratic polynomial that enforces our constraint.
IQP method: Consider the following iterative procedure to enforce constraints in a QUBO form, where in each step we add a slack variable to the QP used, until this enforcement is possible. Hence, given a target constraint, we execute the following step-by-step algorithm:
-
•
STEP 0: Apply the method with no-slack variables given in Eq. (9). If it provides a solution, halt.
- •
- •
-
•
…
- •
where the AND and OR logical operations are performed over .
In the following, we use the expression “IQP method” to refer to this algorithm. Notice that, if this procedure halts at a certain step, then the same amount of slack variables are enough for the IQP method to encode the target constraint as a QUBO penalty term. Remember that we start from STEP 0, which corresponds to the no-slack variables scenario. Hence, this method is designed to enforce our constraint via a QP (7) employing the least number of slack variables. In the following, the QP obtained through this method could be simply called IQP.
We point out that, as we increase the number of slack variables throughout these steps, gets more and more free-parameters at disposal (see Eq. (8)), while at the same time we apply the same amount of “strong” equality conditions and more “soft” inequality conditions (see Eq. (10)) in the corresponding linear system. An indicative estimate of the slack variables necessary for this kind of method is given by requiring that the free parameters in (see Eq. (8)) is larger than the number of sub-constraints, which in the general case are , namely the conditions in the linear system (9) considered with no slack variables. Nonetheless, as we show it later with an example, a lower number of slack variables is often sufficient. Hence, it should already be clear that this method is particularly useful for small . Indeed, while for the standard method the number of slack variables do not change with , in general the number of slack variables needed by the IQP method increases with .
Finally, we underline two more qualities of this technique. First, this method can be extended to any linear/non-linear equality/inequality constraint defined with integer/non-integer parameters. Indeed, consider the linear systems (10) and (13) defining the parameters of the quadratic polynomial enforcing the constraint. The proposed procedure makes no reference to the form of and, more importantly, its linearity is never exploited. The only operation that we have to perform is to verify whether is violated or not for each . It follows that, in case we need to enforce any other type of constraint , the same procedure described in this section can be adopted, where we simply replace “ satisfies/violates ” with “ satisfies/violates ”. A second merit of the IQP method is that, differently from the standard approach, it does not require all the variables to be coupled to each other: some of the parameters can be null. We encounter this scenario for several constraints defining the optimization problem studied later in this work, namely the MPBS problem (see Section VI). This scenario can improve the performance of NISQ devices used as QUBO solvers as sparser formulations require fewer and less-connected qubits.
III.2 MS methods
Consider the scenario where two (or more) constraints are defined over the same subset of variables , where the first has been enforced in a quadratic form either with the IQP method or any other method, e.g., as in Section II.1. We call this constraint the master constraint, while the second, which has yet to be enforced, as satellite. We follow by presenting a variation of the IQP method suited to enforce satellite constraints. This technique, as we show below, employs even less slack variables, if compared with the method presented in Section III.1, which instead is better suited for master constraints.
The main difference between the standard IQP method and its variation for satellite constraints is that we do not impose sub-constraints for all the possible combinations of the binary variables , but only for the strictly smaller subset of combinations satisfying the master constraints. The final discussion of the previous section should already convince the reader that having less than sub-constraints allows to implement satellite constraints with even fewer slack variables.
We start by considering two constraints sharing the same variables , where we call “ True” the constraint considered as master, and “ True” the satellite one, where both can either be linear/non-linear equality/inequality constraints defined with integer/non-integer parameters. Requiring , obtained form Eq. (7), to enforce the satellite constraint with slack variables corresponds to verify whether there exists a solution of the following linear system in the parameters , and :
(14) |
or, equivalently,
(15) |
The step-by-step IQP method described in the previous section can be generalized to the satellite scenario, where we simply restrict the corresponding conditions for to the sub-constraints that satisfy the master constraints, as shown in Eqs. (14) and (15), where is the number of combinations of that violate the master constraint. Hence, all the corresponding steps can be formulated similarly. For instance, the feasibility problem formulation (13) reads:
(16) |
where the AND and OR logical operations are performed over . Remember that this procedure can be straightforwardly applied to constraints that are non-linear or not defined through integer parameters, as described in Sec. III.1.
Notice that, in case all the combinations of satisfying the master constraint satisfy the satellite constraint too, then there is no need to enforce the satellite constraint with a quadratic polynomial at all. Indeed, in this case a solution of Eq. (16) would be given (without employing any slack variable) by for all , , and .
In order to make the notation clearer, we considered the case of one master and one satellite constraint. In case of multiple master constraints, we have to require sub-constraints only for those satisfying all the master constraints. Hence, in such a scenario we would impose even less sub-constraints, namely , where now is the number of combinations that violate at least one master constraint. Naturally, if we have more than one satellite constraint, the solution of this linear system has to be performed independently for each one of these constraints.
The natural consequence of not imposing any sub-constraint for those violating the master constraints is that we may obtain an accidental incentive for some , where violates the master constraint. In this situation such a combination would be promoted, and therefore we need a strategy to avoid that the optimal solution of the corresponding QUBO problem violates one or more constraints. This goal is achieved by amplifying the weight of master constraints as shown in Section III.3.
We say that this MS method adopts a bipartite approach as, even in case of more than two constraints involved, each constraint is either considered as a master or a satellite constraint. Instead, in case of three or more constraints, we can either adopt the strategy introduced in this section or its concatenated version, presented in Section III.4.
III.3 Tuning of the relative multipliers
Consider the QUBO form (3) of a generic optimization problem (2), where a single constraint has been enforced with the IQP method presented in Section III.1. Since the corresponding quadratic form provides penalties/no-penalties if and only if violates/satisfies the corresponding constraint, no adjustments have to be applied to the form given in Eq. (3). The same applies if multiple constraints are enforced with the standard IQP method, namely when the MS method is not exploited.
In case two or more constraints are enforced with the MS method, some adjustments to Eq. (3) have to be applied. In practice, we need to amplify the role of the master constraint quadratic form as follows:
(17) |
where the master and satellite constraints considered here are two constraints picked from the inequality and equality constraints of Eq. (2) that we decided to enforce with the MS approach, hence sharing the same subset of binary variables . Hence, replaces two of the quadratic polynomials on the r.h.s. of Eq. (3). The multiplier can be defined as follows:
(18) |
where and, if the maximization on the r.h.s. is greater than zero, we say that at least one of the satellite constraints provide an accidental incentive. In case of no accidental incentives, .
Thanks to this transformation, we finally achieve a QUBO form of these two constraints that faithfully represents the corresponding constraints. Indeed, the r.h.s. of Eq. (17) is smaller or equal than -1 if violates either the master or the satellite constraint and is zero-valued for at least one value of if both constraints are satisfied:
(19) |
Consider the case where we enforce multiple constraints as master and satellite (either with the IQP method or not) defined over the same shared set of variables . Both master and satellite constraints can be chosen among the equality and inequality constraints. We apply the following transformation of the QUBO form:
(20) |
where, even if not made explicit, the slack variables employed for the master and satellite quadratic polynomials are independent. Therefore, this transformation involves quadratic polynomials picked from Eq. (3). Similarly to Eq. (18), we use:
(21) |
For what concerns the multiplier that appears in Eq. (3), we can apply the ordinary strategies that can be found in the literature. Several exemplary problems for which this parameter is calculated are given in Ref. [5], where at times a different parameter is associated to each quadratic penalty term . In what follows, we define a class of optimization problems, called locally-constrained, for which we propose a different tuning of these parameters.
III.4 Master-satellite concatenation
We briefly describe a method to concatenate the MS approach in order to reduce even more the number of considered sub-constraints, and therefore of slack variables. This approach can be applied solely for those cases where three or more constraints share the same subset of variables . For instance, consider the four constraints , , and (we omit the “” notation) defined over the same binary variables.
The approach of the previous sections requires to identify two groups of constraints, where all the constraints in the second group are satellite of the constraints in the first group. For instance, we could choose the master constraints to be and , while their satellite constraints are respectively and . Each constraint in the satellite group is actively enforced solely for those combinations of that satisfy the constraints from the master group. On the other side, all the master constraints have to be enforced for all the combinations of .
We concatenate the MS method among all the constraints defined over . In particular, we start by choosing two constraints, e.g., and , and enforce as the master constraint of . Hence, would be enforced for all the combinations of , while only for those satisfying (as the master-satellite approach described in the previous sections in the case of two constraints). Nonetheless, differently from before, we continue this concatenation of constraints by considering as the master constraint of . Hence, would be enforced only for those satisfying . Finally, would be enforced as a satellite constraint of all the other constraints, namely with being its master constraint. Therefore, has to be enforced for those already satisfying all the other constraints. This hierarchy of constraint enforcement allows to consider fewer sub-constraints if compared with the non-concatenated case, and therefore, as show before, less slack variables.
This approach requires the introduction of multiple relative multipliers that have the role to counterbalance the accidental incentives coming from the enforcement of , and . If we consider the previously described example of four constraints, we would need three different multipliers similar to from Eq. (18) via the following replacement of the quadratic polynomials enforcing in Eq. (3):
(22) |
where is the IQP obtained while considering the constraint in this concatenated scheme, namely as a satellite of all the constraints with and master of those with . Remember that each can either be an equality or an inequality constraint. The tuning of the parameters is given by the following generalization of Eq. (18):
(23) |
where .
IV Locally-constrained optimization problems
A class of optimization problems that is particularly suited for the application of our techniques is given by the locally-constrained ones, namely those having constraints that are defined node-by-node and depend solely on the variables associated to the arcs incident to each node (see Definition 1). The optimization problems of this type assume the form:
(24) |
where the equality and inequality constraints and are defined solely over those binary variables associated to the arcs incident to the corresponding node and is the subset of nodes identified by the arcs selected by . For instance, if , where corresponds to an arc connecting and , would consists of the nodes and . Hence, for each node in we impose equality constraints and inequality constraints 222Extending our approach to node-dependent and is straightforward but, for the sake of simplicity, we consider each node having the same number of equality and inequality constraints.. Notice that, while in Eq. (2) and referred to the total number of equality and inequality constraints in the whole graph, respectively, in Eq. (24) they refer to the number of equality/inequality constraints per node 333If Eq. (2) is locally-constrained, it can be written in the form of Eq. (24) via the rescaling ..
Locally-constrained optimization problems are particularly appropriate for the application of the IQPMS methods. These problems present a natural subdivision of the constraints defined over the same set of binary variables . Indeed, each node has constraints defined solely over the variables associated to the arcs incident to . Therefore, the application of the IQPMS methods is extremely direct.
Consider a locally-constrained problem with at least two constraints per node and divide them into master constraints and satellite constraints , where . All the linear equality constraints, since they do not need the employment of slack variables, are promoted to master constraints, while the remaining ones can either be satellite or master. We can adopt the IQPMS methods introduced in Sections III.1 and III.2 to enforce the master and satellite constraints as quadratic penalty terms and , which are defined over the variables associated to the arcs incident to the node . Therefore, the QUBO formulation of the optimization problems in this class, thanks to the IQPMS methods, assume the form:
(25) |
where the quantities and have two different roles. The relative multipliers ensure that, for each combination of the binary variables associated to the node , if violates one of the constraints, the net contribution of the term inside the parenthesis is smaller than or equal to -1. Indeed, we need to amplify the effect of master constraints accordingly to the accidental incentives that provide. Instead, if satisfies the constraints, the term inside the parenthesis has to be zero for at least one value of and non-positive for the remaining combinations of . These quantities have the same role of the parameter defined in Eqs. (18) and (20). Instead, the parameters replace the role of in Eq. (3). Hence, they have to amplify enough the penalties in order to suppress those that violate one or more constraints.
IV.1 Tuning of the relative multipliers and the MS concatenation
We first discuss the tuning of the relative multipliers and later we focus on . A straightforward application of the results presented in Section III.3 allows us to consider the following definition
(26) |
where . Notice that all the relative multipliers of the same node have the same value.
In what follows, we describe three proposals for the tuning of , which we call local, neighbour and global. Our first approach is to choose a value of such that the minimum penalty applied when at least one constraint in is violated is (in modulo) greater than the maximum contribution provided by the arcs incident in to . In order to do this, we define to be the restriction of to the node . More precisely, if we write , in we include only those terms where either , or both belong to , namely corresponds to arcs in . Hence, we define:
(27) |
where .
The next two approaches that we define naturally provide larger values of . On the one side it helps preventing the “non-local” violation of constraints that may occur in presence of a local tuning. In Appendix C we describe this phenomenon in the context of the Max-Profit Balanced Settlement optimization problem: an NP-hard problem that we analyse in Section V. Nonetheless, the presence of penalty terms that are too large compared to those found in can suppress the output quality of several QUBO solvers, especially quantum ones. Moreover, depending on the structure of , the evaluation of these two extra methods may be more demanding (this is not the case for diagonal ). Therefore, in this work we privilege the use of the local tuning, and, in case of necessity, we pick larger values of .
The second possible tuning that we discuss extends the previous approach to those nodes neighbouring , namely those sharing an arc with . Similarly to , we now define to be the restriction of to the binary variables associated to arcs incident to and its neighbouring nodes. Hence, we define:
(28) |
where .
The last discussed tuning is labelled as global because it is not node-dependent and it makes sure that the penalty applied is always large enough whenever any constraint is violated. Indeed, it consists in setting:
(29) |
where . Whereas this is a sure choice to make our QUBO formulation exactly encoding the starting optimization problem, namely they provide the same solution, it has some downturns. Depending on the structure of , the optimization may have a computational cost that is similar to the whole optimization problem, and therefore evaluating may become excessively demanding. More importantly, this choice can easily produce a QUBO problem where the penalty terms are (in modulo) extremely larger than the objective function terms. This scenario can deeply affect the solution quality of QUBO solvers and, moreover, we expect that in many common situations this brute-force approach is in fact unnecessary.
To conclude this section, we give an exemplary reformulation of the results regarding the concatenated IQPMS methods to the case of locally-constrained problems. Imagine to have four constraints per node, hence defined over the same binary variables , and concatenate them as in Section III.4. We can reformulate the QUBO form (25) of this locally-constrained optimization problem as:
(30) |
where is the IQP obtained while considering the constraint in this concatenated scheme, namely as a satellite of all the constraints with and master of those with , for . The tuning of the parameters is given by the following generalization of Eq. (23):
(31) |
Regarding the multipliers appearing in Eq. (30), we can adopt the same techniques discussed above in this section.
V Case of study: MPBS
We present the locally-constrained optimization problem that we focus on in this work, which is inspired by a real-world financial scenario. Consider a network of users having a set of pending receivables that still have to be paid for. Each receivable indicates the details of a precise pending payment and therefore it is defined by its value, the corresponding creditor and debtor users and some additional temporal features, such as the dates the receivable has been submitted and when the corresponding payment falls due. These scenarios are often handled by receivable funders who have the role of anticipating the resolution of receivables: users are requested to pay a fee in order to have instant access to a lump sum of capital, which significantly eases the cash-flow issues associated with receivables. Nonetheless, a different strategy for receivable resolution has been proposed in Ref. [29], where a do-ut-des mechanism is exploited. The main goal in the design of this network-based receivable funding is to select the largest-amount subset of receivables for which autonomous payment among users is possible, without involving the funder. This optimization has to be performed under two main constraints. First, the account balance of each user has to be upper- and lower-bounded by some pre-fixed values before and after the resolution of the chosen receivables. We call this constraint CAP/FLOOR, where the CAP (FLOOR) value corresponds to the upper (lower) bound of the account balance. Secondly, each user accepts to realize at least one anticipated payment for which they are debtor if and only if they also receive at least one payment corresponding to a receivable for which they are creditor. Hence, this constraint, called IN/OUT, prevents users from only receiving or only executing the payment of receivables.
To reformulate this optimization problem in a more mathematical framework, we say that each day, we have a different set of receivables that constitutes the pending transactions, namely those not resolved the day before together with those added in the meantime between the two days. Therefore, each receivable is defined by its amount , the corresponding debtor and creditor , together with several temporal features which we do not exploit in this analysis [29]. Indeed, we focus on the single-day resolution problem defined by the R-multigraph induced by , formulated as:
Definition 2 (Receivables R-multigraph).
The R-multigraph induced by the set of receivables consists of the quadruple , where is a set of nodes, is a multiset of ordered pairs of nodes called arcs, and assigns positive numerical values to arcs and assigns four attributes to each node. The nodes represent the users in , the arcs represent the receivables in , where represents the corresponding debtor and the creditor, and is the value of the receivable . Each node is assigned attributes , respectively the corresponding receivable balance, actual balance, cap and floor values.
Hence, given a set of receivables , we can identify a corresponding R-multigraph. The aim of the optimization problem that we are about to define is to maximize the overall value of the receivables paid under the condition that the constraints CAP/FLOOR and IN/OUT described above are satisfied by all users. More precisely:
Definition 3 (Max-Profit Balanced Settlement).
Given the R-multigraph induced by the receivables , find
(32) | |||||
s.t. | (34) | ||||
where is a generic subset of all the arcs/transactions and . Condition (34) confines the balance of the users that take part to the solution multigraph to be inside a finite region, so that floor and cap conditions are not violated after the realization of the selected transactions:
(35) |
Condition (34) imposes that each user taking part to the solution multigraph has to be debtor and creditor for at least one transaction:
(36) |
We start by underlying the locally-constrained nature of this optimization problem: each constraint is defined node-by-node and only depends on the arcs incident to that same node. Secondly, the Max-Profit Balanced Settlement (MPBS) problem is NP-hard as it can be reduced from the SubsetSum problem [29], which is a special case of the Knapsack problem. This problem presents the same structure of other well-known minimum-cost circulation problems, which have been extensively studied for decades (see, e.g., [49]). In the following, we show how the IQPMS methods can be used to reduce the amount of slack variables needed for a QUBO formulation of MPBS with respect to the standard procedure of Section II.1. A simple example of this problem is depicted in Figure 3.
It is obvious that, if a node has either no incoming or outgoing arc in the input R-multigraph, the optimal solution of the corresponding MPBS problem cannot include any arc incident to , otherwise IN/OUT() would be violated. Moreover, the exclusion of all the arcs incident to may cause another node to not have any incoming or outgoing arc. In the following, in order to make our exposition clearer, we assume that all the nodes of the considered input R-multigraphs have at least one incoming and one outgoing arc.
V.1 QUBO formulation of MPBS
In order to cast MPBS into a QUBO form, we proceed as follows. First, we associate to each arc a binary variable . Hence, we obtain a string of binary variables , where is the total number of arcs in . Given the one-to-one relation we can rewrite any expression in Definition 3 with the binary variables . For instance, if we use the definition , the maximization of the objective function of MPBS can be formulated as:
(37) |
where we forced the previous notation by defining: . From now on, we solely consider the quadratic polynomial form of QUBO elements, e.g., , instead of the corresponding matricial form, e.g., .
Our goal is to achieve a QUBO form of MPBS defined over the logical variables and slack binary variables , hence considering the optimization of a quadratic form defined over variables that corresponds to the instance of MPBS identified by :
(38) |
Given the one-to-one relation between subsets of arcs and binary variables combinations , we define:
(39) |
where is the solution of MPBS. We say that the quadratic polynomial provides a QUBO formulation of MPBS if
(40) |
Hence, in order to achieve this formulation, we must translate the constraints of MPBS, namely CAP/FLOOR and IN/OUT, through quadratic penalty terms and , respectively, in order to obtain:
(41) |
MPBS is locally-constrained. Indeed, constraints are clearly defined node-by-node. As we show later, we can achieve the following QUBO structure similar to Eq. (25):
(42) |
where IN/OUT() is considered the master constraint of the node , while CAP/FLOOR() as the corresponding satellite constraint.
We define to be the subset of arcs incident to the node , where is the subset of incoming arcs where represents the creditor and is the subset of outgoing arcs where represents the debtor. Hence, . We use the notation and . Moreover, we use and to indicate specific subsets of the binary variables . For instance, “” stands for all the values of such that corresponds to a arc. Instead, corresponds to the amount of slack variable needed to enforce a constraint associated to the node , where the specific constraint will be clear from the context. We remember that each constraint is enforced with independent slack variables.
We conclude this section by expressing the CAP/FLOOR and IN/OUT conditions with the binary variables associated to the arcs/pending transactions . First, it is straightforward to check that Eq. (35) corresponds to:
(43) |
where and . Moreover, the IN/OUT condition for the node given in Eq. (36) is equivalent to the following pair of inequality constraints:
(44) |
V.2 Standard QUBO formulation of IN/OUT
We start by showing how the standard approach discussed in Section II.1 applies to cast the IN/OUT constraints into a QUBO form. Each IN/OUT constraint in Eq. (44) can be considered as a different inequality constraint , with , of a locally-constrained problem (see Eq. (24)).
It is easy to check that, for each , the IN/OUT inequality constraints have the same minimum . A straight application of the standard procedure described in Section II.1 would involve the introduction of a quantity analogous to (see Eq. (4)), where the amount of slack variables is chosen such that can assume all the integer values in the interval . In general, the larger is this interval, the more slack variables are necessary (see Eq. (6)). For instance, we could enforce into a QUBO form through:
(45) |
where is given by Eq. (5) and assumes all the integer values in . Nonetheless, this path would employ more slack variables than necessary.
Indeed, the combinations of for which violate and, similarly, the combinations of for which violate . Hence, requiring to be able to assume the value , so that Eq. (45) can be zero-valued when , is not necessary, as there would be anyway a violation of the other inequality constraint . Instead, the same is not true for those such that or are equal to the second-lowest possible value, namely .
As a consequence, we consider the QUBO quadratic penalty that corresponds to as Eq. (45), where is given by Eq. (5) and assumes all the integer values in . Hence, the implementation of the first IN/OUT inequality constraint requires slack variables (see Eq. (6)).
In summary, in order to implement a quadratic penalty term that enforces IN/OUT() in a QUBO formulation of MPBS, we can use the standard technique described in Section II.1 and construct , where is given in Eq. (45) and can be obtained similarly. We underline that, while these two addendum share the same logical binary variables, the slack variables are disjoint, namely those involved in do not enter in the definition of and vice-versa. The slack variable cost of this implementation for each node is:
Hence, the number of slack variables needed to implement the IN/OUT constraint in the whole R-multigraph is:
(46) |
In Table 1 we report the amount of slack variables required by the standard method to enforce IN/OUT() in a QUBO.
We notice that the amount of slack variables needed for the QUBO formulation of this constraint strictly depends on the topology of the multigraph considered. Indeed, the amount of slack variables given by Eq. (46) is a function of the number of incoming and outgoing arcs of each node.
V.3 Standard QUBO formulation of CAP/FLOOR
We proceed by applying the standard method discussed in Section II.1 to convert CAP/FLOOR() into a QUBO form. This constraint can be considered as one of the linear inequality constraints in Eq. (2), where now has to assume all the integer values in . Hence, to obtain the corresponding QUBO form, we use Eq. (5) to define the matrix that appears in Eq. (41) and we obtain:
(47) |
Hence, the implementation of this quadratic penalty term requires slack variables for each node (see Eq. (6)) and therefore a total of
(48) |
slack variables to implement the CAP/FLOOR constraint in the whole R-multigraph.
Hence, when adopting the standard method presented in Section II.1, differently from the IN/OUT case, the slack variable cost for this constraint does not depend on the topology of the graph, i.e., the amount of incoming and outgoing arcs per node, but it only depends on the values of and . Notice that, in case of multigraphs taken from real-world financial scenarios, the order of magnitude of can easily be around . Therefore, the amount of slack variables needed per node, even in case of few incoming and outgoing transactions, can be larger than . Later, when considering real instances of MPBS, we set in order to need 4 slack variables per node for the enforcement of this constraint.
V.4 MPBS on quantum annealers
We can plug the newly formulated quadratic penalties for CAP/FLOOR and IN/OUT given respectively in Eqs. (47) and (45), into Eq. (42) and obtain the standard QUBO formulation of MPBS. As a result, the number of binary variables over which we optimize is given by variables , which are associated to the arcs of the input R-multigraph, slack variables to enforce the IN/OUT (see Eq. (46)) and for CAP/FLOOR (see Eq. (48)).
One possibility to solve a QUBO problem with quantum technologies is via a quantum annealer, e.g., those commercially produced by D-wave [31]. In the best-case scenario, namely when each QUBO variable is not coupled with too many variables and the topology of multigraph of arcs is particularly fortunate, D-wave annealers employ around one qubit per QUBO variable. Hence, by using the standard technique described above, the qubits employed to include CAP/FLOOR and IN/OUT are in general much more than those representing the actual transactions. As a result, any attempt to use classical or quantum QUBO solvers in order to solve MPBS for larger and larger R-multigraphs is highly limited by the number of slack variables required by this procedure.
While in the near future quantum processing units with more qubits and higher connectivities may be introduced and therefore better computing performances may be achieved, any reduction of the slack variables needed to implement MPBS as a QUBO makes the use of quantum annealers more suitable to solve instances of this problem. In case of more advantageous strategies, such reductions could give the opportunity to approach real-world problems with a larger input size. An approach that could be adopted together with the slack variable reduction in order to efficiently use quantum annealers is given by splitting large QUBO problems in smaller sub-QUBOs, e.g., by considering the algorithm qbsolv [50]. This method allows to provide sub-optimal solutions to QUBO problems that would not fit entirely on a quantum chip.
In the following sections we first analyse the features of the realistic R-multigraphs studied in Ref. [30] and later we show why they are particularly suited for the application of our newly introduced techniques, namely the IQPMS methods. Real-world instances of MPBS are indeed characterized by a low number of arcs per node and therefore, as discussed in Section III.1, the adoption of IQPMS is particularly effective. Step-by-step, we show how a drastic slack variables reduction is possible. Thanks to our approach, all the nodes with 3 or less arcs require no slack variables in order for IN/OUT() and CAP/FLOOR() to be imposed in a QUBO form. Moreover, for 4-arc nodes, we need only 1 slack variable to implement IN/OUT(). Instead, for nodes with 5 arcs, we need either 1 or 2 slack variables for IN/OUT() and from 0 to 2 slack variables for CAP/FLOOR(). We summarize these slack variable costs in Table 2. Hence, compared to the standard method presented in Section II.1, our methods implement significantly less slack variables and this amount depends only on the multigraph topology: no dependency from the problem parameters occurs. We remind that the standard method enforces CAP/FLOOR() in a QUBO form with a slack variable cost that directly depends on (see Eq. (48)). Finally, we provide various approaches to tune the corresponding penalty multipliers.
V.5 Realistic datasets
In [29, 30] the authors show how to obtain sub-optimal solutions of MPBS for real scenarios obtained by the pan-European bank company UniCredit. The authors consider time-dependent R-multigraph, which vary day-by-day depending on the transactions resolved each day and the newly submitted ones. The average number of arcs per node obtained by averaging over a time span of two years is (see Table 2 from [30]):
(49) |
where these averages vary depending on the particular problem parameters chosen, e.g., cap(), fl() and the lifetime that each receivable is allowed to stay in the R-multigraph. Instead, the single days with the largest ratios are such that:
(50) |
The authors of [30] also analyse the cases where the CAP constraint is dropped (while FLOOR is still enforced), namely for cap(). Nonetheless, the corresponding values of Eqs. (49) and (50) do not change significantly.
It is evident that a technique that allows a reduction of the number of slack variables needed to enforce MPBS constraints for nodes with a number of arcs inside the ranges described here would be of fundamental importance in order to use QUBO solvers. As anticipated, we indeed realize such a protocol which works nicely for nodes with 5 or less arcs. All the nodes with 6 or more arcs should be treated separately and would in general need more slack variables, but the techniques introduced in this work, namely the IQPMS methods, can still be exploited.
VI IQPMS methods for MPBS
We show how the techniques introduced in this work, namely the IQPMS methods, apply to obtain a QUBO representation of MPBS. We apply the IQPMS methods, where IN/OUT() is enforced as master constraint and CAP/FLOOR() as satellite. We summarize the results of this section in Table 2, where we provide the number of slack variables needed by our techniques to enforce the constraints of MPBS in the QUBO form. As we show, the quadratic polynomials found by the IQPMS methods enforcing IN/OUT() and CAP/FLOOR() can be non-completely connected, namely some of the parameters defining the corresponding polynomials (7) can be null. Hence, if the quadratic forms enforcing IN/OUT() and CAP/FLOOR() do not couple the same two or more variables for one or more node , we get a QUBO where the logical variables are not completely connected and therefore solving such a problem is in general easier than in the completely-connected case. Similarly, the same applies to the slack variables, which are not completely connected to each other and to the logical variables. Remember that the standard method in general employs several more slack variables and it completely connects all the logical and slack variables to each other.
VI.1 Reduction of slack variables for IN/OUT
An interesting feature of this constraint is that it is homogenous among the graph, namely, modulo a relabelling of the binary variables involved in , all the nodes having the same values of and have the same . Moreover, the same is true if we swap the values of and : the IN/OUT penalty term for a node with and can be simply obtained from the IN/OUT penalty given to another node having and . Indeed, this penalty requires each node to have either none or at least one incoming and one outgoing arc, while the values associated to each arc and the attributes and do not influence this constraint. This regularity is the main reason why we pick IN/OUT as master: once that we enforce this constraint for a given pair of and , we automatically obtain the IN/OUT() penalty term for all the nodes with either the same and or swapped values of and .
We proceed by exposing how the IQP method for master constraints applies to IN/OUT. As anticipated, this method is particularly effective when is not too large and therefore we fix our attention to those nodes with arcs.
Nodes with 2 arcs: Consider a node with . Without loss of generality, we associate the binary variable to the incoming arc and to the outgoing one. The IQP method provides the following solution with no slack variables
(51) |
where are the combinations that satisfy IN/OUT, while violate IN/OUT.
Nodes with 3 arcs: Consider a node with . We either have and or and . Hence, we simply say that we have a 1vs2 scenario and associate to the type of arc that appears only once. More precisely, if , we associate to the incoming arc and to the outgoing ones. Otherwise, we associate to the outgoing arc and to the incoming ones when .
Interestingly, even if the free parameters in the IQP with are fewer than the independent conditions that have to be applied on , respectively , the IQP method provides infinite many solutions with no slack variables. Remember that the standard method required 2 slack variables to enforce the same constraint (see Section V.2). Hence, we replace some inequality sub-constraints with equality sub-constraints in Eq. (9) in order to have a larger penalty for the string that violates IN/OUT and is likely to provide a large incentive to the objective function of the QUBO formulation of MPBS, namely for . Hence, the following linear system of 8 equations in 7 variables
(52) |
provides the following unique solution:
(53) |
Nodes with 4 arcs: Consider a node with . In the following, the 1vs3 case refers to the and scenario and the and scenario. Instead, the 2vs2 case refers to the and scenario.
The IQP method needs 1 slack variable in order to enforce IN/OUT as a QUBO penalty term for this type of nodes, which we call . Remember that the standard method required 4 slack variables to enforce the same constraint (see Section V.2). Before providing the explicit solutions obtained for , we explain the ordering chosen for the binary variables. Concerning the 1vs3 case, we associate to the type of arc that appears only once in , while are associated to the other type of arcs 444For instance, if and , we associate to the only outgoing arc of the node.. Instead, for the 2vs2 case, are associated to the incoming arcs and to the outgoing ones.
Finally, the corresponding quadratic forms are:
(1vs3) | (54) | ||||
(2vs2) | (55) |
Concerning the sparsity of these penalty terms, we notice that the zero-valued parameters in the solution IQPs (7) are in the 1vs3 scenario and in the 2vs2 scenario.
Nodes with 5 arcs: Consider a node with arcs. In the following, the 1vs4 case refers to the and scenario and the and scenario. Similarly, the 2vs3 case refers to the and scenario and the and scenario.
The IQP method applied to the 1vs4 case allows a solution with 1 slack variable, which we call , while for the 2vs3 case the IQP method requires 2 slack variables, which we call and . Remember that the standard method required 4 and 6 slack variables, respectively, to enforce the same constraint (see Section V.2). Before providing the explicit solutions obtained for , we explain the ordering chosen for the binary variables. Concerning the 1vs4 case, we associate to the type of arc that appears only once in , while are associated to the other type of arcs. Instead, for the 2vs3 case, are associated to the type of arcs that appears twice in , while are associated to the other type of arcs.
Finally, the corresponding quadratic forms are:
Concerning the sparsity of these penalty terms, we notice that the zero-valued parameters in the solution IQPs (7) are in the 1vs4 scenario and in the 2vs3 scenario.
VI.2 Reduction of slack variables for CAP/FLOOR
We follow by enforcing CAP/FLOOR in the QUBO form as a satellite constraint of IN/OUT, which has already been treated as a master constraint. In this section, we construct the quadratic penalty terms with the IQP method by requiring that it provides penalties if satisfies IN/OUT() and violates CAP/FLOOR(). Moreover, this quadratic form has to give no-penalties for at least one combination of whenever satisfies IN/OUT() and CAP/FLOOR(). In other terms, we require to penalize the combinations violating CAP/FLOOR() solely when they already satisfy IN/OUT() (see Section III.2). It follows that such a could provide accidental incentives when IN/OUT() is violated and therefore we have to proceed as described in Section III.3. Hence, we look for a IQP (7) that satisfies the linear system (14), which now assumes the form:
(58) |
Before proceeding, we briefly analyse the slack variables cost expected for IQP method in this scenario. The total number of free parameters inside a IQP grows quadratically with (see Eq. (8)), where is the number of slack variables employed. In case of master constraints, the total number of combinations of for which we have to enforce sub-constraints are . The advantage of considering CAP/FLOOR() a satellite constraint of IN/OUT() is that we can ignore those sub-constraints violating IN/OUT(), which are 555The number of combinations of such that there are no outgoing (incoming) arcs activated are . The number of combinations of such that there are no incoming arcs activated are . Since the null combination is allowed by IN/OUT(), we obtain that the number of combinations violating IN/OUT() are .. Hence, the total number of sub-constraints that have to be imposed for CAP/FLOOR() are:
(59) |
As discussed in Section III.1, an indicative estimate for the minimum number of slack variables necessary to IQP method to enforce a given constraint is given by that such that the number of free parameters (see Eq. (8)) is larger than the number of sub-constraints (see Eq. (59)), namely for:
From the values of and reported in Table 3, we expect that one or more slack variables are necessary to enforce CAP/FLOOR() in a quadratic form with the IQP method only for . As we show later, this is indeed true. Moreover, by generating several random graphs, a major part of the 5-arc nodes studied required either zero or one slack variable and very few needed two slack variables.
In case of we need no slack variables to enforce CAP/FLOOR() in the QUBO form with our methods and moreover we have a closed form for the corresponding polynomials . Moreover, since the linear systems (58) leave free many parameters of the IQP with no slack variables, we can insert back some sub-constraints violating IN/OUT() while still obtaining a solution in a closed form. In particular, for we can add 2 conditions back, and therefore we do not need to assume IN/OUT() to be satisfied, for we can add 3 conditions back and for we can add 3 conditions back in the 1vs3 scenario and 1 in the 2vs2 scenario. Instead, for what concerns the case, we could not find a closed form of . Therefore, in the numerical simulations presented below, the corresponding IQPs have been calculated node by node. In addition, while looking for these polynomials, we incentivised solutions with fewer slack variables and with sparser representative matrices.
In the following, we provide the quadratic forms that are the solutions of the linear system (58) for . These penalties solely depend on the particular combinations that satisfy IN/OUT() and violate CAP/FLOOR(). Hence, given a node , we simply have to check which are these combinations and make a simple substitution inside a parametrized polynomial, without the need to solve any linear system. We remember that for this constraint the standard method required slack variables for each node. Consider that, if we consider realistic scenarios as those analysed in [30], can easily be greater than .
Nodes with 2 arcs: In this case CAP/FLOOR() and IN/OUT() can be enforced at the same time with a single polynomial. It is enough to find an IQP in two variables satisfying the following linear system:
where if satisfies CAP/FLOOR() and otherwise. Remember that always satisfies both CAP/FLOOR() and IN/OUT() while and always violate IN/OUT(). If satisfies CAP/FLOOR(), we obtain . Instead, if violates CAP/FLOOR, we have . Indeed, we can express the solution IQP as follows:
(60) |
Nodes with 3 arcs: We consider the ordering introduced in Section VI.1 for which is associated to the type of arc that appears only once in , while correspond to the other type of arcs. Similarly to the case, we express the linear system (58) as a function of the CAP/FLOOR() penalties/no-penalties that have to be assigned to those combinations satisfying IN/OUT(). In particular, we consider the system:
(61) |
where if satisfies CAP/FLOOR() and otherwise. Notice that, as anticipated, we added three extra sub-constraints in order to completely determine all the free parameters of . Indeed, we required penalties for , which violate IN/OUT(). The solution of such system is:
(62) | |||||
Nodes with 4 arcs: We consider the ordering introduced in Section VI.1 for which, in the 1vs3 scenario is associated to the type of arc that appears only once in , while correspond to the other type of arcs. Instead, for the 2vs2 scenario we consider associated to the incoming arcs, while are associated to the outgoing ones. Similarly to the case, we express the linear system (58) as a function of the CAP/FLOOR() penalties/no-penalties that have to be assigned to those combinations satisfying IN/OUT(). In particular, for the 1vs3 case we consider the system:
(63) |
where if satisfies CAP/FLOOR() and otherwise. Notice that we added three extra penalty sub-constraints for , which violate IN/OUT(). The solution of such system is given in Appendix B.
Instead, for the 2vs2 case we consider the system:
(64) |
where is defined analogously. Notice that we added the extra penalty sub-constraint for , which violates IN/OUT(). The solution of such system is given in Appendix B.
Nodes with 5 arcs: We randomly generated multigraphs and, by using the IQP method for satellite constraints, we never needed more than 1 slack variable to enforce CAP/FLOOR() in the 1vs4 scenario and 2 slack variables for the 2vs3 scenario. Moreover, we often encountered instances where even less slack variables were needed, namely 0 slack variables in the 1vs4 scenario and 0 or 1 slack variables in the 2vs3 scenario. Differently from the cases, we could not achieve a closed form for parametrized by some quantities analogous to the considered above. Hence, node by node, we had to solve the linear system (58).
VI.3 Tuning of MPBS penalty multipliers
We showed how to implement CAP/FLOOR() and IN/OUT() conditions with a reduced number of slack variables, if compared to the standard method, for nodes with arcs. In this section we show how to set the multipliers and in order to be sure that the combinations that solve the maximization in the QUBO encoding the MPBS problem satisfy CAP/FLOOR() and IN/OUT(). In other terms, we want to faithfully represent MPBS, and therefore to satisfy Eq. (41). In the following, we use the corresponding single-node contribution:
(65) |
First, we show how to set and then we proceed with . From now on, we assume that our QUBO formulation of MPBS is made with the IQPMS methods.
The tuning of these multipliers has to take into account the possibility of particularly unfortunate scenarios that may occur. Hence, and have to be set large enough so that also these scenarios get the right penalty factors that allow to obtain the correct optimal result (see Eq. (39)). We explained a possible approach to this problem in Section III.3, where we did not assume any particular structure for the target locally-constrained problem. By studying a particular problem, such as MPBS, it may occur that its particular features allow to set smaller multipliers. While in the following we limit to apply the formulas obtained in Section III.3, in Appendix C we prove that the worst-case requirements that could be find within instances of MPBS require these same multiplier bounds.
First of all, we can set the IN/OUT() relative multiplier by using Eq. (26). The role of the satellite master and satellite constraints are now played by IN/OUT and CAP/FLOOR, respectively. Therefore, we obtain:
(66) |
where . For what concerns , the local, neighbour and global tunings described in Section III.3 leads, respectively, to:
(67) | |||||
(68) | |||||
(69) |
where and neigh() is the set of nodes sharing at least one arc with . We defined the total amount of the transactions associated to as . Instead, is proportional to the sum of the transactions in and its neighbouring nodes. Finally, is node-independent and is proportional to the total amount of transactions in the graph.
Concerning the possibility to consider a global tuning of , this choice surely imposes CAP/FLOOR() and IN/OUT() for all the nodes of the R-multigraph and therefore, if we modify accordingly to these new multipliers, we get . Indeed, as soon as one node violates either CAP/FLOOR() or IN/OUT(), a penalty greater (in modulo) than the value of all the transactions is applied and therefore we are ensured that CAP/FLOOR() or IN/OUT() are satisfied for all the graph nodes.
Nonetheless, this approach generates very large penalties which may compromise the quality of the solution due to the evident unbalance between the weight of the incentives and the penalties for each . Indeed, large unbalances between the terms of a QUBO may be problematic when using quantum annealers, e.g., D-wave. For these reasons, we believe that Eqs. (67) and (68) are the best pick to start studying the MPBS problem within this QUBO framework.
VII MPBS simulations on D-wave quantum annealers
We tested how the standard approach to translate constraints in a quadratic form compares to the IQPMS methods when they are adopted to solve instances of MPBS with D-wave quantum annealers, which are mainly designed to solve QUBO problems. We generated several random R-multigraphs with different numbers of arcs and nodes. Each identified a different MPBS problem, which we solved exactly by exhaustion, namely by checking all the possible combinations of in Eq. (2). Hence, we used the standard and the IQPMS methods to realize a QUBO formulation of each MPBS problem and solved them with two different D-wave quantum annealers, Advantage_system4.1 (Adv1) and the prototype Advantage2_prototype1.1 (Adv2), respectively with 5627 and 563 operating qubits. The corresponding connectivities, namely the maximum number of qubits interacting with each qubit, are equal to 15 for Adv1 (Pegasus topology) and 20 for Adv2 (Zephyr topology) [31]. Finally, we compared the rates of success for which these annealers were able to find the corresponding optimal solutions together with the analysis of other performance figure of merits. Since all the instances have been solved on both quantum annealers, our results allow to compare their performances. For another comparative study of D-wave annealers, see, e.g., Ref. [53].
We follow by describing the features of the R-multigraphs generated. We considered 5 different arc-settings, each one with a different number of arcs, which were . We made this choice because, how we see later, within this range we can appreciate how the standard method goes from being satisfactory to almost useless. For each setting, we generated 4 instances of R-multigraphs with different numbers of nodes, for a total of 20 instances. Hence, fixed a setting, or number of arcs, the different instances had different arcs/nodes ratios, ranging from 1.5 to 2.2. Hence, the average number of arcs incident to each node was . In order to generate an R-multigraph with a given number of arcs and nodes, we first randomly distributed the minimal amount of arcs that allowed each node to have at least one incoming and one outgoing arc. The remaining arcs were then distributed randomly until the desired number of arcs was achieved, where we prevented this procedure to generate nodes with . We assigned integer values to arcs which were randomly sampled from the interval . Finally, we chose and for all the generated instances.
After the generation of these MPBS problems, we exactly solved them by exhaustion. Hence, for each instance we collected the corresponding optimal combination and the optimal value (see Section V.1). Moreover, in order to estimate how restrictive the constraints of this problem are and what is the probability to randomly sample a combination of that simply satisfies our constraints, we explicitly count how many combinations satisfy IN/OUT() and CAP/FLOOR() for all . While the possible combinations grow exponentially with , the amount of that do not lead to a constraint violation did not grow appreciably. After this analysis, we generated the QUBO representations of all the problems both with the standard and the IQPMS methods. Finally, the tuning of the penalty multipliers was executed as described in the previous sections, where we adopted the local approach and we chose . The main features of the instances generated together with the corresponding amount of combinations satisfying the constraints and the number of variables needed to generate the described QUBO representations can be found in Table 4. Moreover, in Figure 4(left) we show how the average number of QUBO variables needed changes when we increase the number of arcs in the instances considered. We show that the standard method needs approximately 4.88 variables per arc added to the graph, while the IQPMS methods need only 1.30 variables.
The number of variables needed by the IQPMS methods are much fewer, in general less than a third. Moreover, the ratio between the variables needed with the standard and the IQPMS methods increases with . One could object that these results are highly influenced by our arbitrary choice of . This is in true, for better or for worse. Indeed, this choice caused the standard method to need 4 slack variables per node, as described in Section V.3. Nonetheless, as discussed above, in realistic scenarios could easily assume values in the order of and therefore the standard method would need more than 15 slack variables per node only to implement CAP/FLOOR. Notice that our methods do not depend on the particular values of and , as shown in Table 2, and therefore its performances would not change in realistic scenarios. Hence, our pick tried to find a compromise between the possibility to obtain QUBO formulations of MPBS with the standard method that were solvable by D-wave, namely not too large, and the idea of being somehow faithful with the possible scenarios encountered in realistic scenarios, namely not too small. Nonetheless, even if we would pick smaller interval , and therefore requiring even less slack variables for the standard QUBO implementation of CAP/FLOOR, we remind that the IQPMS methods always require zero slack variables for this constraint when we have arcs. Moreover, the IQPMS methods need anyway less slack variables to implement the master constraint IN/OUT (compare Tables 1 and 2).
We followed by submitting these QUBO problems to the D-wave quantum annealers Advantage_system4.1 and Advantage2_prototype1.1, where the latter is the latest quantum annealer prototype having enhanced performance. In the following, we simply call these two annealers Adv1 and Adv2, respectively. For each one of the 20 instances of MPBS described above, we performed 4k annealing runs with Adv1 and 4k annealing runs with Adv2. These runs were performed both for the QUBO forms obtained with the standard and the IQPMS methods. Moreover, the 4k runs dedicated to each instance, annealer and method were subdivided in 16 blocks of 250 runs. The motivation of this choice is the following. Every time we ask the D-wave online services to perform a given number of runs, they generate a particular embedding that maps our QUBO problem into an initialization of their quantum chip and all the runs are then performed with that particular embedding. Hence, our motivation to split the runs in 16 different blocks was to average out the possible advantages/disadvantages coming from particularly fortunate/unfortunate embeddings generated by D-wave. As a result, we performed a total 320k runs on D-wave quantum annealers. The annealing times have not been tuned and therefore all the runs had the same default duration of .
We show the average number of qubits needed for these embeddings depending on the method and the annealer adopted, for different number of arcs in Figure 4 (right). We notice that, when considering the standard method, D-wave embeddings need approximately 10.5 qubits per arc introduced into the graph. Instead, when considering our method, we only have approximately 2 qubits per arc introduced into the graph. Several details regarding the embeddings of these problems can be found in Table 5. In particular, depending on the number of arcs in our instance, annealer used and method adopted, we report the average number of qubits employed, the average number of qubits employed per QUBO variable of the instance and the corresponding average maximum chain lengths. The maximum chain length of a given embedding is the maximum number of qubits needed by D-wave to embed a single QUBO variable into the quantum chip. Indeed, since D-wave annealers have limited connectivity, which are 15 for Adv1 and 20 for Adv2, it often happens that more than one qubit is needed to embed a single QUBO variable. Naturally, since the IQPMS methods require far fewer slack variables per node, if compared with the standard method, each binary variable is connected with fewer other variables and therefore our methods allow a much more efficient consumption of qubits. This is one of the main goals achieved by this work, namely providing an encouraging scaling for the qubits deployed by the embedding of locally-constrained QUBO problems into D-wave quantum annealers.
We recorded the number of successes obtained by Adv1 and Adv2 in the individuation of the MPBS optimal solutions. Moreover, we performed the same analysis by counting how many times the D-wave annealers were able to provide outputs that satisfied the MPBS constraints and gave a value of the objective function that was at least the of the optimal value , where is the optimal solution, which are referred to as the “95% solutions”. We summarize these results in Table 6, where we also include the ratio between the number of successes obtained with the IQPMS and the standard methods. These gains, both for the optimal and the 95% solutions, result to be monotonically increasing with the number of arcs in the MPBS instances. Finally, in Figure 5 we plot the number of successes obtained in the described scenarios.
VIII Discussion
We provided a novel set of techniques to translate optimization problems into a QUBO form. In particular, we showed their potential to drastically reduce the number of slack variables and to increase the sparsity of the generated QUBO problems. The first of the two techniques, the IQP method, is particularly useful whenever the constraints to be enforced are defined over few logical variables and/or the values of the parameters inside inequality constraints are large. This method does not necessary require the locally-constrained graph problem structure. Nonetheless, in this case we can easily adopt the MS method, which provides the simultaneous enforcement of multiple constraints defined over the same set of binary variables. Moreover, we underline that our method has the same efficiency whether the constraints considered are equality/inequality linear/non-linear constraints defined with integer/non-integer parameters. Finally, we showed how to tune the corresponding constraint multipliers and we gave a generalization of the MS method for those cases where at least three constraints share the same set of variables.
We followed by considering an NP-hard problem from finance, namely MPBS, and we showed how to apply our methods to this real-world optimization problem. Noteworthy, we obtained closed forms for the implementations of its inequality constraints INOUT() and CAPFLOOR() for all the nodes with arcs, with the only exception of CAPFLOOR() when , which we evaluated case by case. Hence, all the constraints for which we achieved a closed form will not need any additional calculation in order to find the corresponding QUBO implementation: the only required operation is the identifications of the node binary variables combinations that violate/satisfy the constraints.
We generated several instances of this problem for different numbers of nodes and edges involved. First, we quantified how better our methods work in terms of slack variables employed and later we showed how these improvements provide a much higher output quality when quantum annealers are the chosen QUBO solver. Indeed, we considered two D-wave quantum annealers, namely Advantage_system4.1 and Advantage2_prototype1.1, comparing our methods with a standard procedure that provides QUBO forms. In addition, this study provides a performance comparative study between the considered annealers.
We believe that our approach is particularly useful in the context of quantum computation as all the pre-processing applied is devoted to make the resulting QUBO as light as possible and we do not require additional rounds of computation on the QUBO solver. Indeed, quantum annealers in general provide higher-quality outputs when problems are defined in a more compact form, namely when fewer and less-connected qubits are necessary. Hence, focusing on the refinement of QUBO formulations is clearly one key ingredient to obtain better near-term quantum computation results. Classical computing strategies that aim to obtain solutions for these kind of combinatorial problems within shorter times, usually undertake completely different paths [29], avoiding QUBO formulations.
There are several research directions that could be explored starting from our results. First of all, the adoption of IQPMS to solve other combinatorial problems. Possible candidates are the MultipleKnapsack, the Multiple-ConstrainedKnapsack and the MultipleSubsetSum problems. Moreover, given the high performances obtainable in case of locally-constrained optimization problems, its employment could easily fit for real-world combinatorial problems defined over sparse networks (see, e.g., [40]).
Another interesting direction would be given by comparing our methods with those from Ref. [35]. Both output quality from quantum annealing and running times could be compared in order to understand under which regimes one could outperform the other. Indeed, while our methods may employ more slack variables, a fair comparison would allow the IQPMS methods to consider much longer quantum computation times. Moreover, it is not clear whether the number of rounds necessary to the iterative strategies are more sensitive to an increase of the input size or the number of constraints. If the latter would be the case, given realistic instances of MPBS are highly constrained, namely given by sparse graphs with node-by-node constraints, we believe that our methods should provide higher success rates. A similar comparison with the techniques presented in Ref. [40] would be fruitful as well. In summary, we expect that, if we compare the already available strategies with our novel methods in the context of solving locally-constrained lowly-connected optimization problems, we should obtain optimal solutions in shorter quantum computation times.
Future applications of the IQPMS methods could be employed with more sophisticated techniques, e.g., the tuning of the annealing times, the adoption of D-wave classical-quantum hybrid QUBO solvers or particular pre/post-processing procedures, e.g., the variable reduction proposed in Ref. [54], where an iterative algorithm fixes several binary variables to values that have a high probability of being optimal. Another path is given by the use of different NISQ devices, such as in the context of the Quantum Approximate Optimization Algorithm [12] or the Quantum Alternating Operator Ansatz [55], which are designed for logical-gate quantum computers. We remind that the aim of this work is to compare in the fairest way our methods with a well-established standard procedure. Hence, the setup considered has not the intent to prove how far our algorithms can be pushed to achieve optimal solutions via quantum annealing for larger and larger instances, but it rather focuses on giving a clear motivation of when and why we should adopt our novelties. Therefore, we are curious to discover up to which input size our methods, in conjunction with other techniques, allow to provide optimal solutions of MPBS instances defined via realistic datasets.
Acknowledgements
We thank Francesco Gullo and Ilaria Bordino for useful discussions. Dario De Santis acknowledges support from the research project “Dynamics and Information Research Institute - Quantum Information, Quantum Technologies” within the agreement between UniCredit Bank and Scuola Normale Superiore di Pisa (CI14_UNICREDIT_MARMI). Salvatore Tirone acknowledges support from projects PRIN 2017 Taming complexity via Quantum Strategies: a Hybrid Integrated Photonic approach (QUSHIP) Id. 2017SRN-BRK and PRO3 Quantum Pathfinder. Vittorio Giovannetti acknowledges financial support by MUR (Ministero dell’Istruzione, dell’Università e della Ricerca) through the following projects: PNRR MUR project PE0000023-NQSTI.
References
- Hauke et al. [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Perspectives of quantum annealing: methods and implementations, Reports on Progress in Physics 83, 054401 (2020).
- Yarkoni et al. [2022] S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, Quantum annealing for industry applications: introduction and review, Reports on Progress in Physics 85, 104001 (2022).
- Abbas et al. [2023] A. Abbas, A. Ambainis, B. Augustino, A. Bärtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, N. Franco, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, V. Markov, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. O’Meara, E. P. Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. C. B. Symons, S. Tornow, V. Valls, S. Woerner, M. L. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, Quantum optimization: Potential, challenges, and the path forward (2023), arXiv:2312.02279 [quant-ph] .
- Karp [1972] R. M. Karp, Reducibility among combinatorial problems, in Complexity of Computer Computations: Proceedings of a symposium on the Complexity of Computer Computations (Springer US, Boston, MA, 1972) pp. 85–103.
- Lucas [2014] A. Lucas, Ising formulations of many np problems, Frontiers in Physics 2 (2014).
- Rajak et al. [2023] A. Rajak, S. Suzuki, A. Dutta, and B. K. Chakrabarti, Quantum annealing: an overview, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 381 (2023).
- Das and Chakrabarti [2008] A. Das and B. K. Chakrabarti, Colloquium: Quantum annealing and analog quantum computation, Rev. Mod. Phys. 80, 1061 (2008).
- Farhi et al. [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem, Science 292, 472 (2001).
- Farhi et al. [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014), arXiv:1411.4028 [quant-ph] .
- Khumalo et al. [2022] M. T. Khumalo, H. A. Chieza, K. Prag, and M. Woolway, An investigation of ibm quantum computing device performance on combinatorial optimisation problems (2022), arXiv:2107.03638 [quant-ph] .
- Chieza et al. [2020] H. A. Chieza, M. T. Khumalo, K. Prag, and M. Woolway, On the computational performance of ibm quantum devices applied to combinatorial optimisation problems, in 2020 7th International Conference on Soft Computing & Machine Intelligence (ISCMI) (2020) pp. 260–264.
- Zhou et al. [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices, Phys. Rev. X 10, 021067 (2020).
- Ebadi et al. [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J.-G. Liu, R. Samajdar, X.-Z. Luo, B. Nash, X. Gao, B. Barak, E. Farhi, S. Sachdev, N. Gemelke, L. Zhou, S. Choi, H. Pichler, S.-T. Wang, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum optimization of maximum independent set using rydberg atom arrays, Science 376, 1209 (2022), https://www.science.org/doi/pdf/10.1126/science.abo6587 .
- Nguyen et al. [2023] M.-T. Nguyen, J.-G. Liu, J. Wurtz, M. D. Lukin, S.-T. Wang, and H. Pichler, Quantum optimization with arbitrary connectivity using rydberg atom arrays, PRX Quantum 4, 010316 (2023).
- Cain et al. [2023] M. Cain, S. Chattopadhyay, J.-G. Liu, R. Samajdar, H. Pichler, and M. D. Lukin, Quantum speedup for combinatorial optimization with flat energy landscapes (2023), arXiv:2306.13123 [quant-ph] .
- Altshuler et al. [2010] B. Altshuler, H. Krovi, and J. Roland, Anderson localization makes adiabatic quantum optimization fail, Proceedings of the National Academy of Sciences 107, 12446 (2010), https://www.pnas.org/doi/pdf/10.1073/pnas.1002116107 .
- Dickson and Amin [2011] N. G. Dickson and M. H. S. Amin, Does adiabatic quantum optimization fail for np-complete problems?, Phys. Rev. Lett. 106, 050502 (2011).
- Orús et al. [2019] R. Orús, S. Mugel, and E. Lizaso, Quantum computing for finance: Overview and prospects, Reviews in Physics 4, 100028 (2019).
- Mugel et al. [2022] S. Mugel, C. Kuchkovsky, E. Sánchez, S. Fernández-Lorenzo, J. Luis-Hita, E. Lizaso, and R. Orús, Dynamic portfolio optimization with real datasets using quantum processors and quantum-inspired tensor networks, Phys. Rev. Res. 4, 013006 (2022).
- Venturelli and Kondratyev [2019] D. Venturelli and A. Kondratyev, Reverse quantum annealing approach to portfolio optimization problems, Quantum Machine Intelligence 1, 17 (2019).
- Herman et al. [2022] D. Herman, C. Googin, X. Liu, A. Galda, I. Safro, Y. Sun, M. Pistoia, and Y. Alexeev, A survey of quantum computing for finance (2022), arXiv:2201.02773 [quant-ph] .
- Buonaiuto et al. [2023] G. Buonaiuto, F. Gargiulo, G. De Pietro, M. Esposito, and M. Pota, Best practices for portfolio optimization by quantum computing, experimented on real quantum devices, Scientific Reports 13, 19434 (2023).
- Egger et al. [2020] D. J. Egger, C. Gambella, J. Marecek, S. McFaddin, M. Mevissen, R. Raymond, A. Simonetto, S. Woerner, and E. Yndurain, Quantum computing for finance: State-of-the-art and future prospects, IEEE Transactions on Quantum Engineering 1, 1–24 (2020).
- Weinberg et al. [2023] S. J. Weinberg, F. Sanches, T. Ide, K. Kamiya, and R. Correll, Supply chain logistics with quantum and classical annealing algorithms, Scientific Reports 13, 4770 (2023).
- Sales and Araos [2023] J. F. A. Sales and R. A. P. Araos, Adiabatic quantum computing for logistic transport optimization (2023), arXiv:2301.07691 [quant-ph] .
- Hernández et al. [2020] F. Hernández, K. Díaz, M. Forets, and R. Sotelo, Application of quantum optimization techniques (qubo method) to cargo logistics on ships and airplanes, in 2020 IEEE Congreso Bienal de Argentina (ARGENCON) (2020) pp. 1–1.
- Zinner et al. [2021] M. Zinner, F. Dahlhausen, P. Boehme, J. Ehlers, L. Bieske, and L. Fehring, Quantum computing’s potential for drug discovery: Early stage industry dynamics, Drug Discovery Today 26, 1680 (2021).
- Glover et al. [2022] F. Glover, G. Kochenberger, R. Hennig, and Y. Du, Quantum bridge analytics i: a tutorial on formulating and using qubo models, Annals of Operations Research 314, 141 (2022).
- Bordino and Gullo [2018] I. Bordino and F. Gullo, Network-based receivable financing, in Proceedings of the 27th ACM International Conference on Information and Knowledge Management, CIKM ’18 (ACM, 2018).
- Bordino et al. [2020] I. Bordino, F. Gullo, and G. Legnaro, The power of connection: Leveraging network analysis to advance receivable financing (2020), arXiv:2006.13738 [cs.DS] .
- [31] D. S. Inc., D-wave technical report: Technical Description of the D-Wave Quantum Processing Unit, Tech. Rep.
- Boyd et al. [2011] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine Learning 3, 1 (2011).
- Hestenes [1969] M. R. Hestenes, Multiplier and gradient methods, Journal of Optimization Theory and Applications 4, 303 (1969).
- Powell [1969] M. J. Powell, A method for nonlinear constraints in minimization problems, Optimization , 283 (1969).
- Djidjev [2023] H. N. Djidjev, Quantum annealing with inequality constraints: the set cover problem (2023), arXiv:2302.11185 [quant-ph] .
- Yonaga et al. [2020] K. Yonaga, M. J. Miyama, and M. Ohzeki, Solving inequality-constrained binary optimization problems on quantum annealer (2020), arXiv:2012.06119 [quant-ph] .
- Mücke and Gerlach [2023] S. Mücke and T. Gerlach, Efficient Light Source Placement using Quantum Computing, MDPI Proc. L 3630, 478 (2023), arXiv:2312.01156 [quant-ph] .
- Verma and Lewis [2021] A. Verma and M. Lewis, Variable reduction for quadratic unconstrained binary optimization (2021), arXiv:2105.07032 [math.OC] .
- van der Linde et al. [2023] S. G. van der Linde, W. van der Schoot, and F. Phillipson, Efficient quantum solution for the constrained tactical capacity problem for distributed electricity generation, in Innovations for Community Services, edited by U. R. Krieger, G. Eichler, C. Erfurth, and G. Fahrnberger (Springer Nature Switzerland, Cham, 2023) pp. 203–221.
- Colucci et al. [2023] G. Colucci, S. v. d. Linde, and F. Phillipson, Power network optimization: A quantum approach, IEEE Access 11, 98926 (2023).
- Kirkpatrick et al. [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671 (1983), https://www.science.org/doi/pdf/10.1126/science.220.4598.671 .
- van Laarhoven and Aarts [1987] P. J. M. van Laarhoven and E. H. L. Aarts, Simulated annealing, in Simulated Annealing: Theory and Applications (Springer Netherlands, Dordrecht, 1987) pp. 7–15.
- Wilhelm and Ward [1987] M. R. Wilhelm and T. L. Ward, Solving quadratic assignment problems by ‘simulated annealing’, IIE Transactions 19, 107 (1987), https://doi.org/10.1080/07408178708975376 .
- Wang et al. [2016] C. Wang, H. Chen, and E. Jonckheere, Quantum versus simulated annealing in wireless interference network optimization, Scientific Reports 6, 25797 (2016).
- Paul [2010] G. Paul, Comparative performance of tabu search and simulated annealing heuristics for the quadratic assignment problem, Operations Research Letters 38, 577 (2010).
- Note [1] If and are enforced by and , respectively, they would make use of and independent slack variables.
- Note [2] Extending our approach to node-dependent and is straightforward but, for the sake of simplicity, we consider each node having the same number of equality and inequality constraints.
- Note [3] If Eq. (2) is locally-constrained, it can be written in the form of Eq. (24) via the rescaling .
- Goldberg and Tarjan [1990] A. V. Goldberg and R. E. Tarjan, Finding minimum-cost circulations by successive approximation, Mathematics of Operations Research 15, 430 (1990).
- Inc. [2017] D. S. Inc., D-wave technical report: Partitioning Optimization Problems for Hybrid Classical/Quantum Execution, Tech. Rep. (2017).
- Note [4] For instance, if and , we associate to the only outgoing arc of the node.
- Note [5] The number of combinations of such that there are no outgoing (incoming) arcs activated are . The number of combinations of such that there are no incoming arcs activated are . Since the null combination is allowed by IN/OUT(), we obtain that the number of combinations violating IN/OUT() are .
- Hahn et al. [2023] G. Hahn, E. Pelofske, and H. N. Djidjev, Posiform planting: generating qubo instances for benchmarking, Frontiers in Computer Science 5, 10.3389/fcomp.2023.1275948 (2023).
- Karimi et al. [2017] H. Karimi, G. Rosenberg, and H. G. Katzgraber, Effective optimization using sample persistence: A case study on quantum annealers and various monte carlo optimization methods, Phys. Rev. E 96, 043312 (2017).
- Hadfield et al. [2019] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12 (2019).
- Note [6] A tighter version of could be formulated, but we do not see any particular advantage in performing the extra computation required.
Appendix A The MS method without the IQP method
We showed how to apply the MS method when constraints are enforced with the IQP method in Section III.2. Despite their combination, namely the IQPMS method, provides various advantages, the MS method can be adopted independently from the IQP method. Here, we show how to employ this method when constraints are enforced in a QUBO form with the standard techniques described in Section II.1. This approach could be useful when the constraints considered are defined over a number of variables that is too large to be handled by the IQP method.
Consider an optimization problem of the generic form (2). Suppose that, among all the equality and inequality constraints and , two of them are defined over the same three binary variables. For instance, consider the case where one equality and one inequality constraints are of the form:
where it is required and . Our goal is to achieve quadratic forms of , and that provides penalties when these constraints are violated. In the context of the MS approach, we are going to construct a penalty function corresponding to the master constraints, which provides a penalty iff it is violated, otherwise it is zero-valued. Instead, the penalty function corresponding to the satellite constraint provides penalties when satisfy the master and violates the satellite constraint and is zero-valued when both constraints are satisfied. Hence, we do not impose any requirement to the satellite penalty function when the master constraint is violated.
We pick as master and as satellite constraint. Hence, in order to have a simultaneous reformulation of these constraints within the MS formalism, we start by enforcing with the standard method (see Section II.1), hence with the quadratic penalty function .
The quadratic penalty function relative to has not to be enforced by considering all the possible values of , but solely those combinations already satisfying , namely . Notice that satisfy , while leads to a violation. Hence, in order to construct the quantity needed to enforce our inequality/satellite constraint (see Eq. (5)), we solely have to worry to build a quadratic form such that for and .
If would be enforced outside the MS formalism, namely as in Section II.1, we would require to assume all the integer values in the interval , where the minimization is performed over all . Hence, in the standard case we would obtain : two slack variables would be required. Now, since has to be enforced as a satellite of , we have to ignore those strings violating the master constraint and therefore the same minimization is now performed solely over , which leads to:
Hence, by using the MS method, can be enforced with a quadratic form without employing any slack variable via , which is equal to iff satisfies the master and the satellite constraints, namely for and is smaller or equal to otherwise. Hence, in this case, there are no accidental incentives and the relative multiplier of the master constraint can be fixed to (see Section III.3). Notice that penalises also those combinations of satisfying the satellite and violating the master constraint. Hence, the sum of the two penalty functions:
is equal to zero iff satisfies both constraints and is smaller or equal to -1 iff one or both constraints are violated, where no slack variables have been employed.
Appendix B IQPMS solutions for CAP/FLOOR() when
Appendix C Worst-case scenarios for MPBS penalty multipliers
Local tuning: Consider the worst-case scenario where the optimal solution (see Eq. (39)) is unique and for the node provides the minimum contribution to the term , while there is a different combination that violates CAP/FLOOR() and provides the maximum contribution to . Hence, we would have that the binary variables of associated to the transactions involving are and . We choose such that is larger for than for :
(72) |
where we used because satisfies the node constraints, and the fact that because satisfies CAP/FLOOR() and IN/OUT(). The exclamation mark above the inequality symbol underlines that, if satisfied, can be considered large enough in this local setting.
From Eq. (72), we derive the lower-bound , where is the total amount of the transactions involved with the node . Notice that this is the same result obtained in Eq. (67). Hence, we fix the notation
(73) |
An estimate for can be obtained similarly. The worst-case scenario is obtained when the restriction of on is , namely providing a null contribution to , while is the combination providing the largest contribution to among those violating IN/OUT() and provides the largest accidental incentive for , namely . We choose such that the term is larger for than for :
The first inequality is justified because satisfies IN/OUT() and therefore , , and for combinations violating IN/OUT(). Moreover, we used and the fact that satisfies the node constraints. A consequence of this first inequality is that we are imposing a looser condition than previously declared 666A tighter version of could be formulated, but we do not see any particular advantage in performing the extra computation required.. Hence, we have the lower-bound Since we assumed , for the general case we can consider:
(74) |
Non-local tuning: We set some lower-bounds (73) and (74) for, respectively, and by imposing that the sub-QUBO that corresponds to the node assumes the largest value when satisfies CAP/FLOOR() and IN/OUT(). We called it local because we are not considering the possible consequences that a violation of CAP/FLOOR() and/or IN/OUT() may have on for some node . Indeed, the possible incentives that these situation may generate could require to set higher values of and .
Imagine the following unfortunate scenario: , in order to respect CAP/FLOOR() and IN/OUT(), has many s for the variables corresponding to the transactions in . As a consequence, the transactions of the nodes neighbouring are not activated too (otherwise CAP/FLOOR() or IN/OUT() may be violated), resulting in many other for the arcs belonging to the nodes sharing at least one arc with . Hence, it may happen that a violation of CAP/FLOOR() and/or IN/OUT() in allows the activation of many other transactions for the nodes neighbouring , where instead there are no violations of CAP/FLOOR() and IN/OUT(). Therefore, there may exists a string that receives a penalty for the violation of CAP/FLOOR() and/or IN/OUT() in , but this penalty is overcompensated by the incentives activated by the neighbouring nodes (without penalties). We define the nodes adjacent to as follows:
(75) |
Therefore, to prevent that violates CAP/FLOOR() and/or IN/OUT() as described, we consider:
(76) |
where and is the the sum of all the transactions involving and the neighbouring nodes. Notice that the multipliers involved in Eq. (76) are obtained by replacing with from the local formulation, namely we have:
(77) |
Global tuning: In order to get a global tuning, we consider the following choice of :
(78) |
and is the sum of all the transaction defined on .