Experimental analysis of quantum annealers and hybrid solvers using benchmark optimization problems

This paper studies the Hamiltonian Cycle Problem (HCP) and the Traveling Salesman Problem (TSP) on D-Wave's quantum systems. Initially, motivated by the fact that most libraries present their benchmark instances in terms of adjacency matrices, we develop a novel matrix formulation for the HCP and TSP Hamiltonians, which enables the seamless and automatic integration of benchmark instances in quantum platforms. our extensive experimental tests have led us to some interesting conclusions. D-Wave's {\tt Advantage\_system4.1} is more efficient than {\tt Advantage\_system1.1} both in terms of qubit utilization and quality of solutions. Finally, we experimentally establish that D-Wave's Hybrid solvers always provide a valid solution to a problem, without violating the QUBO constraints, even for arbitrarily big problems, of the order of $120$ nodes. When solving TSP instances, the solutions produced by the quantum annealer are often invalid, in the sense that they violate the topology of the graph. To address this use we advocate the use of \emph{min-max normalization} for the coefficients of the TSP Hamiltonian. Finally, we present a thorough mathematical analysis on the precise number of constraints required to express the HCP and TSP Hamiltonians. This analysis, explains quantitatively why, almost always, running incomplete graph instances requires more qubits than complete instances. It turns out that incomplete graph require more quadratic constraints than complete graphs, a fact that has been corroborated by a series of experiments.


Introduction
Quantum computers are the strongest contenders against conventional, silicon-based systems and gain traction by the day.Proposed by Richard Feynman in 1982, the idea of using quantum "parallelism" [1] to solve computational problems has been proven true in quite a few cases, like Shor's [2] and Grover's [3] algorithms, which are more efficient than their conventional counterparts.Further information about quantum computing and information is available in [4].
The adiabatic theorem is a staple in quantum mechanics [5,6].According to the theorem, when a quantum system with a starting ground state of H init is subject to a gradually changing Hamiltonian, with the same starting ground state and a final ground state of H f in , the system itself will adapt to the changes, ending at H f in .In the first decade of the 2000s, Farhi et al. introduced adiabatic quantum computation [7,8], basing it on the aforementioned adiabatic theorem.Adiabatic quantum computation has been shown to be equivalent to standard quantum computation [9].
In the early 2000s, Kadowaki and Nishimori introduced quantum annealing, a metaheuristic based on the principles of quantum adiabatic computation [10].Since then, it has gained prominence for the tackling of combinatorial optimization problems.Quantum annealing, similar to conventional simulated annealing [11], uses a multivariable function to create an energy landscape of all the possible states in a quantum superposition, with its ground state representing the problem's solution.Then, the system goes through the quantum annealing process repeatedly until it encounters an optimal solution with a high probability.The difference between quantum and simulated annealing is the multi-qubit tunneling used by quantum annealing [10], which also improves performance due to the high degree of parallelism.Quantum annealers try to calculate the optimal solutions by analyzing every possible input simultaneously, an attribute that might be vital when handling NP-complete problems.However, as of now, these methods are still in development.Thus, it is recommended for them to be used as an automatic heuristic-finding program than as a precise solving program [12].
When solving large scale optimization problems, one of the most effective classical strategies is to search among the nearest neighbors and search for paths between the initial and final configurations to improve upon an initial guess y t .This search takes place locally among neighboring configurations similar to y t .By finding the best solution within the local neighborhood of y t , the next candidate solution y t+1 appears.Then a new local search starts with y t+1 as the starting point in the neighborhood.Unfortunately, greedy local improvement, in case of hard problems, may be deceptively leading the solution into a local minima whose energy may be much higher than the globally minimum value.
Optimization problems are one of the areas where using quantum computers is considered to be advantageous [13,14,15,16].To solve an optimization problem with a quantum computer, we first need to create a Hamiltonian, the ground state of which represents an optimal solution for the problem.In most cases using quantum annealing, a system starts in an equal superposition for all its states with said Hamiltonian applied.Over time the system evolves according to the Schrödinger equation, and the system's state changes depending on the strength of the local transverse field as it changes over time.To get the problem's solution, we slowly turn off the transverse field and, the system settles in its ground state.If we have chosen the correct Hamiltonian, then the ground state of our system will also be its optimal solution.D-Wave offers quantum computing systems using quantum annealing for solving optimization and probabilistic sampling problems.One of the optimization problems that the D-Wave computers have demonstrated the ability to solve is the quadratic unconstrained binary optimization problem.The quantum processing units (QPUs) inside the D-Wave machines handle the quantum annealing process.The lowest energy states of the superconducting loops are the quantum bits (which we shall refer to as qubits) and are the analog of conventional bits for the QPU [17,18].The quantum annealing process takes a system with all of its qubits in superposition and guides it towards a new state in which every qubit collapses into a value of either 0 or 1.The resulting system will be a system in a classical state, which will also be the optimal solution to the problem.In 2020 D-Wave announced their new generation of quantum computers, the Advantage System.The Advantage System uses the Pegasus topology, improving from the previous generation's Chimera topology.The Advantage QPU contains more than 5000 qubits, with each qubit having 15 couplers to other qubits, totaling more than 35000 couplers.Qubits in the Advantage QPU are mapped to a P16 Pegasus graph, meaning they are logically mapped into a 15 × 15 × 2 matrix of unit cells on a diagonal grid.A Pegasus unit cell contains twenty-four qubits, with each qubit coupled to one similarly aligned qubit in the cell and two similarly aligned qubits in adjacent cells [19].The percentage of total working qubits is called the working graph and is a subgraph of the total connected cells that physically exist in the QPU.
One of the types of problems that D-Wave computers can solve is the quadratic unconstrained binary optimization (QUBO for short) problem.A QUBO model is a pattern matching technique with many applications, ranging from machine learning to solving optimization problems [20].Its basis is minimizing a quadratic polynomial function over binary variables [21,22,23,24,25,18,26,21] and thus is an NPhard problem [27].The Schrödinger equation is a linear partial differential equation.It expresses the wave function or state function of a quantum-mechanical system.By knowing the necessary variables in a system, we can model most physical systems.The Ising model is a model that is equivalent to the QUBO model.Proposed in the 1920s by Ernst Ising and Wilhelm Lenz, the Ising model is a well-researched model in ferromagnetism [28,29].It assumes that qubits represent the model's variables, and their interactions stand for the costs associated with each pair of qubits.The model's design makes it possible to formulate it into an undirected graph with qubits as vertices and couplers as edges among them.In 2017, D-Wave introduced qbsolv, an open-source software that, according to D-Wave, can handle problems of arbitrary size.Qbsolv is a hybrid system meaning it breaks down the problem's graph into smaller subgraphs using its CPU and solves each partition using its QPU.It then repeats the process and uses a Tabu-based search to examine if a better solution is available [30].

Related work
The Hamiltonian Cycle Problem (HCP for short) asks a simple yet complex in answering question: "Given a graph of n nodes, is there a path that passes through each node only once and returns to the starting node?"It is an NP-Complete problem, first proposed in the 1850s, and is a fascinating problem for computer scientists.Flinders Hamiltonian Cycle Project (FHCP) provides many benchmarking graph datasets for both the HCP and TSP [31].In this paper we have used their dataset of small graphs, generated by GENREG [32].Their small number of nodes made them suitable for solving by the QPU without needing a hybrid solver.The Traveling Salesman Problem (TSP from now on) is an NP-hard combinatorial optimization problem that builds upon the HCP [33].Given a weighted graph, instead of finding the existence of a Hamiltonian Cycle, the TSP asks what's the cycle with the minimum total cost.
The complexity of the TSP makes it a great topic of interest for many researchers, and leads them to pursue other avenues.One such alternative approach advocates the use of metaheuristics, i.e., high-level heuristics designed to select a lower-level heuristic that can produce a fairly good solution with limited computing capacity [34].The term "metaheuristics" was suggested by Glover.Of course metaheuristic procedures, in contrast to exact methods, do not guarantee a global optimal solution [35].Papalitsas et al. [36] designed a metaheuristic based on VNS for the TSP with emphasis on Time Windows.This quantum-inspired procedure was also applied successfully to the solution of real-life problems that can be modeled as TSP instances [15].More recently, [37] applied a quantum-inspired metaheuristic for tackling the practical problem of garbage collection with time windows that produced particularly promising experimental results, as further comparative analysis demonstrated in [38].A thorough statistical and computational analysis on asymmetric, symmetric, and national TSP benchmarks from the well known TSPLIB benchmark library was conducted in [39].
Another state-of-the-art approach for combinatorial optimization problem like TSP is to employ unconventional computing [40,41].It turns out that HCP and TSP are suitable candidates for formulation into the QUBO model so as to be solved by a quantum annealer.An excellent reference for many Ising formulations for NP-complete and NP-hard optimization problems is the paper by Lucas [29].In this paper, Lucas also presented the QUBO formulations for the Hamiltonian Cycle Problem and the Traveling Salesman Problem, giving particular emphasis to the minimization of the usage of qubits.Another formulation of the Traveling Salesman Problem, focusing on handling the symmetric version of the problem, was presented by Martonák et al. [42].The basis for the QUBO formulation used in this paper was proposed by Papalitsas et al. in [43].A recent work [44] conducted an experimental analysis of the performance of two hybrid systems provided by D-Wave, the Kerberos solver and the LeapHybridSampler solver, using the TSP as a benchmark.The conclusion of the tests in [44] is that the Kerberos solver is superior to the LHS, as Kerberos consistently yields solutions closer to the optimal route for every TSP instance.However, the LHS still produced routes that obeyed the imposed constraints for the TSP, even if the quality of its solutions wasn't as high as Kerberos's.
Contribution.Motivated by the fact that most libraries present their benchmark instances in terms of adjacency matrices, and in order to facilitate their execution by the quantum annealer and by tools like qubovert [45], we set out convert the HCP and TSP Hamiltonians in matrix form.This formulation, which we believe has not appeared previously in the literature, enables the seamless and automatic integration of the benchmark instances in quantum testing platforms.
We also present a thorough mathematical analysis on the precise number of constraints required to express the HCP and TSP Hamiltonians.This analysis, explains quantitatively why, almost always, running incomplete graph instances requires more qubits than complete instances.It turns out that incomplete graph require more quadratic constraints than complete graphs.This theoretical finding has been corroborated by a series of experiments outlined in subsection 5.3.Its importance is not only theoretical, but practical as well, since in the current technological stage, qubits still remain a precious commodity.
When solving TSP instances, the solutions produced by the quantum annealer are often invalid, in the sense that they violate the topology of the graph.This is to be expected because if the weight of an edge is higher than the constraint penalty, the annealer, in order to reach a lower energy state, ignores the constraint.To address this use we advocate the use of min-max normalization of the coefficients of the TSP Hamiltonian.This well-known and easy to implement technique was experimentally tested and found to be considerably helpful.
Furthermore, our extensive experimental tests have led us to some interesting conclusions.We found out that D-Wave's Advantage system4.1 is more efficient than Advantage system1.1 both in terms of qubit utilization and quality of solutions.Using our proposed matrix formulation, it is possible to run the Burma14 instance of the TSPLIB library on a D-Wave's Advantage system4.1, although without obtaining a "correct" solution.Finally, we have established experimentally that D-Wave's Hybrid solvers always provide a valid solution to a problem and never break the constraints of a QUBO model, even for arbitrarily big problems, of the order of 120 nodes.

Organization
This paper is structured as follows.Section 1 gives an introduction to the subject along with the most relevant references.Section 2 introduces and explains the standard QUBO formulation of the HCP and TSP problems.Section 3, through a detailed mathematical analysis, presents the matrix QUBO formulation of the HCP and TSP problems.Section 4 demonstrates the graph weight normalization procedure.Section 5 contains the experimental results obtained from our tests on D-Waves's QPU.Section 6 presents the results from the tests contacted on D-Waves's Hybrid Solver, and section 7 summarizes and discusses the conclusions drawn from the experiments.

The standard QUBO formulation
We recall that given a graph G = (V, E), where V is the set of vertices and E is the set of edges, both HCP and TSP problems involve finding a tour where n is the number of vertices, and p i is the vertex in the i th position of the tour.Typically, QUBO expressions are formed from binary decision variables, i.e., variables that can only take the values 0 or 1.In this context, these binary variables will be designated by x v,p and they will have the usual meaning: The HCP and TSP QUBO forms, as given by [29], are built-up by the following simpler Hamiltonians. (2.4) ) These 4 Hamiltonians, H v,p , H p,v , H E c and H W , encode fundamental properties that capture the essence of the two problems.
1. H v,p asserts that every vertex must appears in exactly one position of the tour.
2. H p,v states that each position of the tour is occupied by precisely one vertex.
3. H E c requires that the tour is comprised of edges that "really" exist.If the tour, mistakenly, contains a "phantom" edge, i.e., an edge belonging to E c , then this will incur an energy penalty.
4. H W computes the cost or weight of the tour.A tour minimizing this Hamiltonian is an optimal tour.
The complete QUBO forms of the HCP and TSP problems are given below.
In the above formulas c 1 and c 2 are positive constants, different, in general.We shall have more to say later about the significance of these constants.The first 3 Hamiltonians, H v,p , H p,v , H E c , express integrity constraints; the slightest violations renders the solution invalid.The fourth Hamiltonian, H W , is the minimization objective, that is necessary in order to converge to the tour associated with the minimal cost.Example 2.1.To provide a clearer picture, we will give a small scale example of formulating a TSP into QUBO, using the graph G 1 shown in Figure 1.The problem's graph has n = 4 nodes, with its adjacency matrix being Figure 2. Using equations (2.3)-(2.8)we can construct the Hamiltonian of the QUBO model for solving the TSP in this graph, as shown below.To emphasize the fact the Hamiltonians computed in this example refer to the particular graph G 1 , we use the notation T SP when referring to them.
In this example, the general expression 1 This expression can be further simplified by taking into account the identity and the fact that the square of a binary decision variable is just the variable to give Therefore, In a similar way, we may construct H G1 p,v .H E c depends on the topology of the graph.In our example, where the graph is complete, H G1 E c would only involve self-loops, that is To avoid any potential confusion, let us point out that in the tour, which is a cycle, position n+1 coincides with position 1, a fact which was used in the formulation of (2.15).As a result of the completeness of the graph, all the above 16 quadratic constraints are already present in H G1 v,p , as can be verified by recalling (2.13).Of course, their presence in H G1 E c increases by 1 their corresponding coefficient, amplifying their relative significance.It is useful to briefly consider what happens when the graph is not complete, i.e., at least one edge that is not a self-loop is missing.If, for instance, G = (V, E) is an undirected graph and (1, 3) ∈ E, which implies that also ( H W captures the topology of the graph and the cost, also referred to as weight, associated with each edge.In this particular example, it will produce the constrains listed below, where we have again identified position n + 1 with position 1.The expansion of H G1 W shown in (2.18) creates 48 new quadratic constrains.It is instructive to see what changes if the graph is not complete, assuming, as before, that the missing edges are not self-loops.If the are m edges present, then each such edge will involve n (here n = 4 quadratic constraints).
In view of (2.7), the Hamiltonian for solving the HCP on Using the previous example as a starting point, it is straightforward to derive the number of constraints that for a general graph with n nodes.At this point, it is expedient to make the following remark.In the QUBO formulation, the binary variables corresponding to an edge (or arc) (u, v) are distinct from the binary variables corresponding to the edge (or arc) (v, u).Therefore, they give rise to distinct constraints in the Hamiltonians H HCP and H T SP .Hence, from this perspective, irrespective of whether a complete graph G = (V, E) with n nodes is directed or undirected, it should be considered as having n(n − 1) edges.In other words, even in a undirected graph, we should view the edges (u, v) and (v, u) as different.• If the graph is not complete, then H E c creates for each missing edge (u, v) = E that is not a selfloop, i.e., u = v, n additional constraints.If k such edges are missing, there will be nk additional quadratic constraints in total.This fact is experimentally validated from the results presented in subsection 5.3.
• For a complete graph with n nodes H W introduces n 2 (n − 1) quadratic constraints.One could, equivalently, assert that in a complete graph H W requires n|E| quadratic constraints, where |E| stands for the number of edges in the graph.This last formula can be generalized further to the case where the graph is not complete and has m edges.In such a case H W will add nm quadratic constraints.
• For a complete graph, H HCP will require n 2 binary variables and n 3 constraints in total (recall that the n 2 linear constraints are the same).If the graph is not complete, it will require even more constraints, specifically n 3 + nk, where k is the number of missing edges.
• For a complete graph, H T SP will require n 2 binary variables, 2n 2 (n − 1) = 2n|E| quadratic constraints, and n 2 (2n−1) constraints in total.If the graph is not complete, it will require n 3 +n(m+k), where m is the number of existing edges and k is the number of missing edges.
• In the typical and practically important case where we have a simple graph, that is a graph with no self-loops and no parallel edges, then m + k = n(n − 1).This leads to the simplification of the previous formulas, which now become 2n 2 (n − 1) for the number of quadratic constraints and n 2 (2n − 1) for the total number of constraints.Let us point out the interesting fact that in this case the number of constraints required for the H T SP Hamiltonian is equal to the number of constraints required when the graph is complete.
Tables 1, 2, and 3 summarize the above observations regarding the number of constraints involved in the HCP and TSP Hamiltonians.
3 The matrix QUBO formulation In this section we shall present the matrix form of the Hamiltonians (2.3)-(2.8)so as to enable their seamless execution of HCP and TSP instances by the quantum annealer and by tools like qubovert [45].
For this purpose, we define for each v ∈ V the column vector X v containing the binary decision variables x v,1 , x v,2 , . . ., x v,n corresponding to the possible positions of vertex v in the tour.Then, by combining the n vectors X v , 1 ≤ v ≤ n, we construct the column vector X, which of course has n 2 rows.
The matrices expressing (2.7) and (2.8) are based on the matrix form of (2.3)-(2.6).To define the latter in the easiest way, we need a few auxiliary n × n matrices.As usual, the n × n identity matrix is denoted by I n , its negation by −I n , and the n × n zero matrix whose entries are all 0 by O n .
Two other auxiliary n × n matrices are J n and K n defined below.
Let us assume that W = [w u,v ] is the n × n adjacency matrix characterizing the topology of the given graph.In this matrix if w u,v = 0 then the edge (or arc) (u, v) does not exist in the graph, whereas if w u,v > 0, then w u,v expresses the cost or weight associated with (u, v).To each entry w u,v of W , we associate two n × n matrices W (u,v) and W (u,v) defined as follows: We proceed now to give the matrix form of the Hamiltonians (2.3)-(2.6).Specifically we establish the correspondence between Hamiltonians and n 2 × n 2 block matrices shown below, where the symbol → is used to convey association and not equality.
In view of the relations (3.10)-(3.13), the sum of the 3 Hamiltonians H v,p + H p,v + H E c can be associated with the linear expression 2n + X (M v,p + M p,v + M E c ) X, where X is the column vector of the binary variables given by (3.2).The constant term 2n is irrelevant in the whole quantum annealing process, and can be omitted.By defining for succinctness the matrix M HCP , the Hamiltonian H HCP can be expressed in matrix form as shown below.
Analogously, H T SP can be cast in a compact matrix form with the help of the auxiliary matrix M T SP .
To demonstrate how the previous machinery can be utilized, we return to the graph G 1 of Example 2.1 and express the HCP and TSP Hamiltonians as matrices.For clarity of presentation and simplicity we will take the 2 positive constants c 1 and c 2 equal to 1.To emphasize the fact the X and the matrices computed in this example refer to the particular graph G 1 , we use the notation HCP and M G1 T SP when referring to them.In this example, the column vector X G1 , containing the n 2 = 16 binary variables, is listed below.Note that in order to save space we actually show X G1 .Matrices M G1 v,p , M G1 p,v , M G1 E c and M G1 W , adopted to the specific graph G 1 , now become as follows.
Before we proceed, let us remark that a straightforward computation of the expressions X M G1 v,p X and X M G1 p,v X will produce exactly the same constraints (minus the constant terms) as those given in equations (2.13) and (2.14), respectively.M G1 v,p and M G1 p,v are rather generic, in the sense that they do not convey any specific information about the graph at hand.The topology of the graph is captured by M G1 E c and M G1 W . Since G 1 is complete, matrix M G1 E c assumes a particularly simple form.As expected, by explicitly computing X M G1 E c X we obtain exactly the same 16 constraints encountered in (2.15).
In this example, the quantitative information contained in the adjacency matrix (2.9), which gives the topology of the graph and the weight of each edge, is incorporated in the M G1 W matrix.It is easy to ascertain that X M G1 W X creates exactly the same 24 constraints as those given by equation (2.18).
Finally, by adding the matrices M G1 HCP and M G1 W as per equation (3.17), and taking c 2 = 1 for simplicity, we derive the matrix M G1 T SP and the matrix form of the Hamiltonian H G1 T SP .

Normalizing graph weights
During the extensive experimental tests of various instances of the TSP, a certain pattern consistently emerged.The aforementioned pattern demonstrates emphatically the importance of the numerical values of the weights used in the TSP Hamiltonian.In almost every occasion where the weights of the graph were above a certain threshold, e.g., 10, the annealer preferred to break one or more of the imposed validity constraints rather than to strictly adhere to the topology of the graph.In hindsight, one could say that this is indeed the rational and expected behavior.If the weight of an edge is higher than the constraint penalty, the annealer, in order to reach a lower energy state, ignores the constraint, thus converging to solutions that violate the topology of the graph.On the bright side, the predictability of this behavior, gives us the opportunity to remedy this problem.The solution that worked best in our experimental tests was the min-max normalization [46] applied to the weights of the matrix M W before adding the constraints.After the normalization process, the weight values are still proportional but lower than the constraint penalties.Specifically, given the n 2 × n 2 matrix M W = [m i,j ], we denote by m min and m max the minimal and maximal elements of M W , where without loss of generality we assume that m min = m max .All elements of M W are modified according to the formula (4.3) and the resulting normalized matrix is designated by N W .Therefore, by defining appropriately the normalized matrix N T SP , encompassing all the normalized constraints, we derive the normalized matrix form of the Hamiltonian H T SP , as shown below.From now on we shall always employ this normalized matrix, without further mention.
Figure 3: The min-max normalization formula.
Figure 4: The min-max normalization algorithm.
Example 4.1.We explain the normalization process by examining how it can be applied to the matrix M G1 W , as given by (3.24).Attempting to solve the TSP using H G1 T SP given by equation (3.26), which is based on the M G1 T SP matrix of equation (3.25), will typically yield a solution that violates some integrity constraints.We can counteract this by applying the normalization Algorithm 1.The final normalized QUBO model is shown below.
After the normalization, the constraints imposed on the QUBO model by the weights of the TSP graph retain their respective importance while assuming a range from 0 (for the minimum weight) to 1 (for the maximum weight).Thus, we can impose the rest of our constraints without conflicts, as the higher values in the previous QUBO model no longer overshadow them, forcing the annealer to ignore them.Figure 5 demonstrates the importance of the normalization process as a tool that amplifies the probability that the solutions produced by the annealer will not violate any of the constraints imposed by the QUBO model.To test its effectiveness We ran a total of 20 experiments on D-Wave's LeapHybridSampler.The normalization procedure was applied to exactly half of the instances, e.g., 10, and the other instances were executed without normalization.As can be seen in Figure 5, LeapHybridSampler failed to yield a valid solution satisfying the graph topology in all 10 instances that weren't normalized, whereas LeapHybridSampler succeeded in producing a valid solution in all 10 normalized instances.

HCP using QPU
As explained in section 1, the Hamiltonian Cycle Problem (HCP) is the basis for the Traveling Salesman Problem (TSP), and its solution is an integral part of the solution of the TSP.We solved the HCP using both QPU systems available by D-Wave, namely Advantage system1.1 and Advantage system4.1.The Flinders Hamiltonian Cycle Project (FHCP) provides many benchmarking graph datasets for both the HCP and the TSP [31].We used their dataset of small graphs, namely the graph instances H 6, H 8, H 10, H 12, H 14, generated by GENREG [32].Their small number of nodes made them fit to run on D-Wave's Quantum Annealer's available qubits, without requiring the use of a hybrid solver.The average metrics that were obtained during the experiments are summarized in Table 4.The "Nodes" column contains the number of nodes of the graphs, the "Solver" column shows the selected QPU that was used to solve the problem, and the "Runs" column gives the total number of times the experiment was executed.The "QAT" column shows the QPU Access Time, that is the average time it took the Quantum Processing Unit to solve the problem over all runs.The "Qubits" column contains the average number of qubits used in the embedding for each run.Finally, the "Success Rate" column depicts the percentage of runs that provided a solution that didn't violate any of the imposed constraints.As can be seen, violations of the imposed constraints start to appear when the number of nodes n ≥ 8.This is in line with similar findings in a recent work [50].We also point out the clear improvements of Advantage system4.1 compared to Advantage system1.1, both with respect to the success rate and the qubit usage.While Advantage system4.1 seems to be relatively slower, it more than makes up for it by offering better overall performance.
Figure 6 gives the average number of qubits used in the embedding as a function of the number of nodes per graph.This visualization shows the graphs for the Advantage system1.1 solver, the Advantage system4.1, and the expected qubits per instance.Both the Advantage system1.1 and the Advantage system4.1 demonstrate a rapidly increasing growth in the number of qubits as n increases.One can also see that near and after the intersection points between the solvers' lines and the expected qubits' line is the point where the QPU starts breaking the imposed constraints.It is also noteworthy that the intersection point of Advantage system4.1 is to the right of Advantage system1.1, something that reflects the first solver's superior success rate.In both cases, after the intersection point, every run on the solver returns a result that breaks at least one of the imposed constraints.Figure 7 provides a comparison of the average QPU Access Time between the two solvers.It is evident that Advantage system4.1 needs more time to solve the problems, even if the difference is minuscule for real-world metrics.However, this relative loss of speed is compensated by the decrease of qubits in the embedding, meaning that Advantage system4.1 can handle bigger instances in a more satisfactory way.

TSP using QPU
Moving from the HCP to the TSP is in principle quite straightforward, as their main difference is the introduction of weights on the graph.In our experiments, we relied on the TSPLIB.The TSPLIB is a library of TSP graph instances, commonly used for benchmarking algorithms that aim to solve the TSP [51].The TSPLIB instance with the least number of nodes is the burma14 graph, with just n = 14 nodes.Table 5 below contains the average metrics that were gathered during the experiments.The "Nodes" column contains the number of nodes of the graphs, the "Solver" column shows the selected QPU that was used to solve the problem, and the "Runs" column gives the total number of times the experiment was executed.The "QAT" column shows the QPU Access Time, that is the average time it took the Quantum Processing Unit to solve the problem over all runs.The "Qubits" column contains the average number of qubits used in the embedding for each run.The "Cost" column shows the average cost of the solution path from each run.Finally, the "Success Rate" column depicts the percentage of runs that provided a solution that didn't violate any of the imposed constraints.The results in Table 5 are quite similar to those obtained in the HCP case.Every run produced a solution that violated at least one of the integrity constraints.The qubit usage on the Advantage system4.1 is visually depicted in Figure 8.The number of qubits used is slightly higher than the optimal number given in Table 1 of section 2. This, coupled with the improvement of the results produced by Advantage system4.1 shown in Table 4, is a clear indicator of the degree of technological progress achieved.

The effect of graph connectivity on qubit usage
When comparing the results from the experiments, a consistent pattern emerged.Specifically, solving the HCP on graphs that were not complete, i.e., there were "missing" edges, always required more qubits than solving the HCP on complete graphs with the same number of nodes.We believe that this fact is the experimental validation of the theoretical results derived in Example 2.1 and quantified in Tables 1, 2, and 3, which stipulate that when solving the HCP, incomplete graphs require more integrity constraints than complete graphs.In particular, the HCP on a graph with n nodes and k missing edges requires nk additional quadratic constraints compared to a complete graph with n nodes.To verify this experimentally, we created graphs of varying connectivity based on the FHCP 14 nodes graph and measured the number of qubits required to solve the HCP with respect to the number of existing edges.Figure 9 shows the number of qubits used in the embedding of the HCP Hamiltonian, as a function of the number of edges of the graph.It is evident that the more connected the graph is the fewer qubits its embedding needs.This is in full compliance with the theoretical analysis outlined in Table 2 of section 2.

Experimenting on D-Wave's Leap's Hybrid Solvers
To determine the viability of solving the HCP and TSP problems on a hybrid solver, we used D-Wave's LeapHybridSampler [52] because it can sample arbitrarily big problems [53].

HCP using QPU-Hybrid
Running the HCP problems using D-Wave's LeapHybridSolver resulted in the metics found in Table 6.The "Nodes" column contains the number of nodes of the graphs, the "Solver" column shows the selected QPU that was used to solve the problem, and the "Runs" column gives the total number of times the experiment was executed.The "QAT" column shows the QPU Access Time, that is the average time it took the Quantum Processing Unit to solve the problem over all runs.The "Run Time" column contains is the average total time (in miliseconds) it took the solver to solve the problem instance.Finally, the "Success Rate" column depicts the percentage of runs that provided a solution that didn't violate any of the imposed constraints.In stark contrast with the results obtained from running the problems directly on the QPU, as shown in Table 4, on the hybrid solver all the runs returned a valid solution that did not violate any of the imposed constraints.No matter the size of the problem, the hybrid solver always returns an acceptable solution.Another noteworthy discrepancy is the variance of the values in the QPU Access Time column.It does not seem to follow any distribution or linear growth but varies wildly.The way D-Wave's hybrid solver work, it is possible that in some runs, the QPU Access Time will be 0, as the computer might solve the problem without using the actual QPU [54].

TSP using QPU-Hybrid
Running the TSP problems using D-Wave's LeapHybridSolver resulted in the metrics found in Table 7.
The "Nodes" column contains the number of nodes of the graphs, the "Solver" column shows the selected QPU that was used to solve the problem, and the "Runs" column gives the total number of times the experiment was executed.The "QAT" column shows the QPU Access Time, that is the average time it took the Quantum Processing Unit to solve the problem over all runs.The "Run Time" column contains is the average total time (in miliseconds) it took the solver to solve the problem instance.The "Cost" column shows the average cost of the solution path from each run.Finally, the "Success Rate" column depicts the percentage of runs that provided a solution that didn't violate any of the imposed constraints.Table 7 shows that while the QPU Access Time is an increasing function of n, the results are again different from table 4 since the hybrid computer did not access the QPU every time.We particularly emphasize that the LeapHybridSolver successfully solved the gr120 problem from the TSPLIB, an instance with 120 nodes.The fact that the resulting paths were far from optimal, can be attributed to the way LeapHybridSolver breaks the problems down to smaller instances, possibly losing important "shortcut" paths.Figure 10 presents the average cost of the solution produced by the LeapHybridSampler, compared to the optimal solution for each problem.One can see that as the number of nodes of the graph increases, the solution produced by LeapHybridSampler deviates from the optimal even more.Additionally, the average cost of the solutions reached by the annealer exhibits significant discrepancies compared to the optimal solutions.This is in line with the findings from a very recent study [50] that also highlighted the difficulty the D-Wave solvers face in finding optimal solutions for TSP instances.

Discussion and conclusions
In this work, we have considered various implementations, computational analyses, and comparisons on different versions of the D-Wave quantum annealer for the HCP and TSP problems.We extended our previous work by providing a complete matrix formulation for the HCP and TSP problems, and by introducing a normalization technique on the weight matrices.This technique helped avoid the constraint generated by the TSP graph's weights overshadowing the other imposed constraints, causing D-Wave's annealers to misbehave.We also ran some experiments on D-Wave's annealers and came to the following conclusions.
• D-Wave's Advantage system4.1 is more efficient than Advantage system1.1 in its use of qubits for solving a problem and provides more consistently correct solutions.
• It is possible to run the Burma14 instance of the TSPLIB library on a quantum annealer using the aforementioned methods, although it can't provide a "correct" solution because of the annealer's limitations.
• The more connected a graph is, the fewer qubits are needed for it to be solved by the quantum annealers, as fewer constraints need to be imposed.
• Hybrid solvers always provide a correct solution to a problem and never break the constraints of a QUBO model, even for arbitrarily big problems.
The above conclusions are particularly significant as they indicate that quantum annealers can solve TSP and HCP problems.With the current advancements in quantum computing, the prospect of using a quantum annealer for solving real-world TSP problems looks very promising, especially if quantum computers continue to scale at the grade we see today.
Our computational results indicate that our method is very promising.We intend to continue this line of research and extend its usage to a variety of problem classes, with emphasis on other variations of TSP problems, such as the Traveling Salesperson Problem.Other possible improvements we intend to evaluate experimentally include the introduction of suitable "phantom edges" on graphs that are not fully connected, in order to circumvent the addition of more constraints, as well as experimenting with more normalization methods for the TSP graph weights.

Figure 5 :
Figure 5: LeapHybridSampler returned a valid solution in all 10 normalized instances, but failed to do so in all 10 instances that weren't normalized.

Figure 6 :
Figure 6: The number of qubits used in the embedding as a function of the number of nodes.

Figure 7 :
Figure 7: The average QPU Access Time as a function of the number of nodes.

Figure 8 :
Figure 8: The average number of qubits used for Burma14 on the QPU.

Figure 9 :
Figure 9: The average number of qubits used as a function of the number of edges in a graph with 14 nodes.

Figure 10 :
Figure 10: Comparison between the average costs of the solution produced by LeapHybridSampler and the optimal cost.
Each of the expansions of H G1 v,p and H G1 p,v , as given in (2.13) and (2.14), creates a constant term, namely 4, and 40 constraints in total: 16 linear constraints and 24 quadratic constraints.The linear constraints are the same, but the quadratic constraints are different.
3, 1) ∈ E, then H G1 E c , apart from the constraints shown in (2.15), would additionally, contain the following constraints.
If, on the other hand, G = (V, E) is a directed graph and (1, 3) ∈ E, which does not necessarily imply that (3, 1) ∈ E, then, H G1 E c , apart from the constraints shown in (2.15), would additionally, contain the following constraints.
where c 1 is a positive constant.H G1 HCP is constructed by combining (2.13) and (2.14), it requires 16 binary variables, and it involves 64 constraints.Analogously, the Hamiltonian for solving the TSP on G 1 is H G1 T SP = H G1 HCP + c 2 H G1 W , where c 2 is also a positive constant.Since H G1 T SP is constructed by combining (2.13), (2.14) and (2.18), it contains 112 constraints.

•
Each of the Hamiltonians H v,p and H p,v , as given in (2.13) and (2.14), require a constant term, namely n, n 2 linear constraints, and n 2 (n−1).It important to clarify that the linear constraints are the same in both H v,p and H p,v , but the quadratic constraints are different.•If the graph is complete, then H E c does not add any new constraints.To be precise, H E c involves n 2 quadratic constraints, which are also present in the H v,p Hamiltonian.Hence, H E c does not add new constraints, it only enhances the relative weight of existing constraints.

Table 1 :
This table summarizes the number of constraints that are involved in each of the Hamiltonians (2.3)-(2.8),when applied to a complete graph with n nodes.

Table 2 :
This table summarizes the number of constraints that are involved in each of the Hamiltonians (2.3)-(2.8),when applied to an arbitrary graph with n nodes, m edges and k missing edges (none of which is a self-loop).

Table 3 :
This table shows the number of constraints that are involved in the H T SP Hamiltonian (2.8), when m + k = n(n − 1).

Table 4 :
This Table contains the average metrics that were obtained from solving HCP instances using D-Wave's QPU.

Table 5 :
This Table contains the average metrics that were obtained from solving TSP instances using D-Wave's QPU.

Table 6 :
This Table contains the average metrics that were obtained from solving HCP instances using D-Wave's LeapHybridSampler.

Table 7 :
This Table contains the average metrics that were gathered from experiments on TSP instances using D-Wave's LeapHybridSampler.