E Bourreau
email: eric.bourreau@lirmm.fr

G Fleury

P Lacomme

Mixer Hamiltonian with QAOA for Max 𝒌-coloring: numerical evaluations

, numerous propositions in OR lies on 1-to-𝑛 mapping functions, including the split algorithm that transform one TSP solution into a VRP solution.

The objective is at first to give a compact and readable presentation of these Mixer Hamiltonians considering the functional analogies that exist between the OR community practices and the quantum field. Our experiments encompass numerical evaluations of circuit using the Qiskit library of IBM meeting the theoretical considerations.

Introduction

Quantum optimization is a new and attractive area with potential significant implication in operation research. Minimization problems can now be investigated using quantum metaheuristics with the promise of a strongly effective approach avoiding trapping into local minima that standard local search. Simulated Annealing based methods efficiency (commonly used in OR field) is performed by slowly reducing to 0 a parameter 𝑡 that permits to follow the potential barriers. The metaheuristic methods differ by the way used to avoid premature convergence to local minima and to conserve a strong capacity in search space investigation. The sequel simulated annealing is an example of a large class of metaheuristics that encompasses (but are not limited to): memetic algorithms, GRASP, VNS… From a quantum mechanics point of view, quantum fluctuations are similar to thermal fluctuation. Quantum mechanics differ from classical approaches considering that waves can tunnel through potential barriers (energy) [START_REF] Martoòák | Quantum annealing of the traveling-salesman problem[END_REF]. In the last few years several quantum metaheuristics have been introduced coming from the quantum physic community, defining a family of quantum approximate algorithms that encompasses for example the sequel Adiabatic based Algorithms that provides an approximate solution of the Schrödinger equation.

Lately (Farhi et al., 2014) have introduced a new class of algorithms based on alternation between two families of operators referred to as Hamiltonian and mixing Hamiltonian leading to Quantum Approximate Optimization Algorithms referred to as QAOA. Such algorithms are hybrid algorithms where the search space investigation is devoted to a classical computer to optimize a set of parameters and the distribution of probabilities evaluation is achieved by a quantum device. This algorithm that does no support local search consideration and that provides a full search space investigation, has been extended in the famed publication of [START_REF] Farhi | Quantum supremacy through the Quantum Approximate Optimization Algorithm[END_REF]. They defined new ansatz that support the exploration of the feasible subspace only, meaning that hard constraints are satisfied by definition. This approach is very similar to classical approaches of the OR community since the exploration of the feasible subspace is achieved by a careful definition of classical operation including for example permutation of qubit into the qubit-string used to model solutions;

QAOA based approach

Quantum Approximate Optimization Algorithms (Farhi et al, 2014) take advantage of alternations between the cost function investigation which is modeled by a Hamiltonian 𝐻 𝑃 from one side and a driver Hamiltonian operator 𝐻 𝐷 . The Quantum Alternating Operator Ansatz [START_REF] Hadfield | Quantum Algorithms for Scientific Computing and Approximate Optimization[END_REF] takes into consideration a general parameterized family of unitary operators to efficiently modelize the Hamiltonian. This algorithm supports a restricted number of states and create an efficient alternative to the Adiabatic Optimization

QAOA Principles

QAOA based approaches can be seen as metaheuristic-based methods taking advantages of timediscretization of adiabatic computing. It is a hybrid algorithm because it combines:  a Quantum parametrized circuit executed on a quantum machine,  a Classical meta-optimization loop executed on a classical machine. As introduced by [START_REF] Schrödinger | An Undulatory Theory of the Mechanics of Atoms and Molecules[END_REF] A specific resolution can be achieved considering an iterative approximation of the Schrödinger solution taking advantages of a very specific 𝐻 Hamiltonian (based on 𝑍 𝑖 , 𝑍 𝑖 𝑍 𝑗 operators). In Adiabatic Quantum Optimization, the system must be tuned in a ground state of one Hamiltonian 𝐻 𝐷 commonly referred to as the "driver Hamiltonian" where this Hamiltonian should not commute with the Hamiltonian that model the combinatorial problem. It requires to slowly tuned from 𝐻 𝐷 to 𝐻 𝑃 using an interpolation based on one parameter 𝑠(𝑡) assumed to smoothly decrease from 1 to 0: 𝐻(𝑡) = 𝑠(𝑡). 𝐻 𝐷 + [1 -𝑠(𝑡)]. 𝐻 𝑃 (Farih et al., 2000).

Adiabatic quantum optimization and all quantum optimization based annealing approaches are limited to small-scale instances due to the slow evolution that is required from 𝐻 𝐷 to 𝐻 𝑃 . In 2014 (Farih et al., 2014) have introduced a new trend of approaches referred to as quantum heuristics leading to a more compact decomposition scheme denoted Quantum Approximate Optimization Algorithm (QAOA). QAOA received a considerable amount of attention including but not limited to (Farih and Harrow, 2019), [START_REF] Rahmani | Optimizing variational quantum algorithms using Pontryagin's minimum principle[END_REF], [START_REF] Jiang | Near-optimal quantum circuit for Grover's unstructured search using a transverse field[END_REF], [START_REF] Wecker | Training a quantum optimizer[END_REF] and [START_REF] Wang | Quantum approximate optimization algorithm for MaxCut: A fermionic view[END_REF].

Lately [START_REF] Farhi | Quantum supremacy through the Quantum Approximate Optimization Algorithm[END_REF] prove that the output distribution obtained at the end of one QAOA algorithm cannot be efficiently approximated using one classical algorithm because of the undue complexity of the required algorithms. Over the past few years QAOA has attracted attention to many researchers, because it seems to be one candidate for demonstration of quantum supremacy as stressed by [START_REF] Farhi | Quantum supremacy through the Quantum Approximate Optimization Algorithm[END_REF]. A quantum supremacy demonstration consists in defining one circuit with results requiring an exponential number of operations into a classical algorithm.

Modelling function to Hamiltonian

QAOA seeks to solve a hard optimization problem i. Note that Hadfield in 2021 (Hadfield, 2021) gives a concise description of rules for composing Hamiltonians representing clauses that model clauses or functions including Boolean formulas. A Hamiltonian is implemented into a quantum circuit by deriving 𝑈 𝐻 (𝑡) = 𝑒 -𝑖.𝐻.𝑡 with 𝑡 ∈ [0; 2𝜋] and using both CNOT and Z-rotations. 𝑡 refers to the weight in the iterative search process of QAOA.

Search space investigation

The 𝛽 ⃗ and 𝛾 ⃗ weights parametrized a quantum state |𝜑(𝛽 ⃗ , 𝛾 ⃗ )⟩ that defines a solution 𝑦 with probability We encompass that QAOA efficiency strongly relies on some key-points (figure 2):

|⟨𝑦||𝜑(𝛽 ⃗ , 𝛾 ⃗) ⟩|
 providing a good ratio between the estimation of quality of 𝐶 𝑝 (𝛽 ⃗⃗⃗⃗⃗⃗ , 𝛾 ⃗⃗⃗⃗⃗⃗ ) versus the number of shots is required and must be tuned carefully.  The last computed distribution |𝜓(𝛽 * ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ must be collected on a subset of solutions strongly smaller than the total number of solutions to avoid a costly inefficient enumerations i.e. the algorithm has converged to the optimal solutions and quasi solutions.

In the following Figure 2 we illustrate two separated search spaces computed in a hybrid way on two separated machines (one classical, one quantum computer).

First, we have to compute analytically a Hamiltonian 𝐻 with 2 angles (𝛽 ⃗⃗⃗⃗⃗⃗ , 𝛾 ⃗⃗⃗⃗⃗⃗ ). Second, we should be capable to convert this as a circuit. Finally we have to estimate by sampling the expected value 𝐶 𝑝 (𝛽 ⃗⃗⃗⃗⃗⃗ , 𝛾 ⃗⃗⃗⃗⃗⃗ ). This give information to investigate new value parameter to compute a new Hamiltonian.

The success of QAOA relies on the algorithm capacity to find a good set of 2. 𝑝 parameters and by consequence an efficient search strategy is required in this variational phase.

Mixer Hamiltonian with QAOA

A wide majority of Operational Research problems can be defined to both minimization of one objective function and satisfaction of problem dependent constraints. Resolution approaches are based on exact methods, heuristics or metaheuristics depending on the problem and the computation time available to provide a solution. For a practical point of view, heuristics and metaheuristics are commonly used in resolution of large instances and they address the following key-points during the search space investigation:

 generation of good initial solutions by powerful constructive heuristics dedicated to the problem;  application of a local search to solutions to find local minima and by consequence favoring convergence;

 diversification mechanisms to avoid premature convergence and search space trap;  Indirect representation of solutions to relax constraints and alternate between several search spaces.

In the next section, we will focus on representations of solutions.

Mapping function in O.R. field

Quantum annealing optimization principles (QA) consists into a smooth decreases of quantum fluctuations to traverse barriers in the energy landscapes investigating the global minima of the function.

Quantum annealers are limited in resolution of optimization problems by the Hamiltonian circuit size and specifically by constraints included in the objective function using penalty terms. Driven Hamiltonian based approaches are designed to avoid this penalty term and the all-to-all connectivity and to defined new Hamiltonian that commutes with the constraints.

Mixer Hamiltonians are tuned to favor investigation of operational researches where the feasible subspace is strongly lower that the full space (set of 𝑛-bit strings used for 𝐻) using a mapping from initial states to final solutions. Investigating only feasible subspace has received attention to the Operation Research (O.R.) community for decades and such approaches have been successful in several research areas leading to very efficient metaheuristics.

In Job-Shop Scheduling a common indirect representation is based on the Bierwith vector (Bierwith, 1996) that can be transform (in 𝑂(𝑛)) into one oriented disjunctive graph that model one solution of the job-shop. We can also cite scheduling problems based on disjunctive graph, the local search approaches take advantages of the longest path (see the block definition of Grabowski for example) [START_REF] Grabowski | A block approach for single machine scheduling with release date"s and due dates[END_REF]. In flow based scheduling problems (RCPSP for example) the local search is based on the cut [START_REF] Artigues | Insertion techniques for static and dynamic resource-constrained project scheduling[END_REF], and in routing problems based on the geometry of solutions. The diversification mechanisms are metaheuristic dependent, for example they are based on the temperature for the simulated annealing (SA) or mutation in memetic algorithms (MA) [START_REF] Moscato | Memetic Algorithms: A Short Introduction[END_REF], and permits to accept transitions from one solution to worst solutions with a non-zero probability to cross the barriers in the objective function landscape.

In VRP a well-known representation is the giant trip that can be transformed into one VRP solution using the Split Algorithm [START_REF] Lacomme | Competitive genetic algorithms for the Capacitated Arc Routing Problem and its extensions[END_REF]. The Split algorithm defines a mapping from the set of giant trip of a TSP (Traveling Salesman Problem) to solutions of the VRP and the metaheuristic based approaches manipulate efficiently the set of giant trips only. [START_REF] Cheng | A tutorial survey of job-shop scheduling problems using genetic algorithms -I representation[END_REF][START_REF] Cheng | A tutorial survey of job-shop scheduling problems using genetic algorithms -I representation[END_REF] are the very first who made a full analysis of non-string coding approach and decoding mechanisms in the global context of constraint optimization.

Their remarks encompass both scheduling and routing by emphasizing the principle of a mapping perspective from the problem constraint in order to have indirect representation that generate either an initial phase and evolutionary process production of feasible solution only.

The split algorithms (see [START_REF] Prins | Order-first split-second methods for vehicle routing problems: A review[END_REF]) are one classical approach in routing by transforming a giant trip of the TSP into (for example) a solution of the VRP.

In this context most order-first/split-second methods generate in a first phase an indirect solution representation, referred to as giant tour or task ordering. Feasible routes are deduced from this indirect representation in a second phase. This approach offers the following advantages: (i) any solution of the routing problem has an indirect string representation (ISR); (ii) using a splitting procedure, each ISR can be decoded into a solution to the original problem, and this splitting can be made optimally and (iii) there exists at least one "optimal" ISR, i.e., one that gives an optimal solution to the original problem after splitting. It has been proved by [START_REF] Prins | A simple and effective evolutionary algorithm for the. Vehicle Routing Problem[END_REF]) that the local search operator must be defined both in modification of a VRP solution including basic moves (2-OPT, insert….) and also moves applied to the giant trip. The global performance of the local search is based on alternation between the two spaces.

The Mixer Hamiltonians are on the same trend of research considering that a careful definition of the mapping function, should permit to transform one qubit-strings into one solution to favor exploration of the "feasible" search space only. Contrary to the OR classical field, the explicit local search dynamic is embedded by the Mixer Hamiltonien.

Mixer Hamiltonian: Mapping function to feasible solutions

A wide and important class of optimization problems seeks to minimize an objective function 𝑓: 𝐹 → ℝ to be optimized under a set of constraints. The major drawback of the classical Hamiltonian formulation concerns the weighted quadratic terms in the objective function meaning that the search space investigation is not limited to "feasible" solutions but also encompasses 𝑥 𝑖 assignments where the constraints do not hold. The second drawback concerns theweighted parameters 𝐴 𝑗 that must be large enough to favor convergence to the feasible solution set, but must permit to avoid to be trapped into local minima.

The Mixer Hamiltonians extend the initial proposal of [START_REF] Hen | Quantum Annealing for Constrained Optimization[END_REF] in the adiabatic context where they identify the 𝐻 𝑗𝑘 = 𝑋 𝑗 . 𝑋 𝑘 + 𝑌 𝑗 . 𝑌 𝑘 permitting to restrict the state evolution to feasible subspace for some problems. Note that 𝐻 𝑗𝑘 and 𝐻 𝑖𝑘 do not commute in general and that trotterization is required using the Masuo's proposition [START_REF] Suzuki | Generalized Trotter's Formula and Systematic Approximants of Exponential Operators and Inners Derivations with Applications to Many-Body Problems[END_REF]. A more compact and readable demonstration has been introduced lately by Barthel et Zhang in 2019 (Barthel and Zhang, 2019) and it concludes by the following proposition.

Proposition. 𝑙𝑖𝑚 𝑛→∞ (𝑒 𝑖.𝐴.𝑡/𝑛 . 𝑒 𝑖.𝐵.𝑡/𝑛 ) 𝑛 = 𝑒 𝑖.(𝐴+𝐵).𝑡

And as a consequence 𝑒 𝑖.(𝐴+𝐵).𝑡 ≈ 𝑒 𝑖.(𝐴).𝑡 . 𝑒 𝑖.(𝐵).𝑡 where 𝐴 and 𝐵 are matrices.

Defining one Mixer Hamiltonian consists in including restriction in search space investigation in the mixing operator, is supposed to reduce the efficiency required in terms of gates and avoiding inclusion of extra weighted terms in the objective.

XY-Hamiltonian analysis

The 𝑋𝑌-Hamiltonian

𝐻(𝛾) = ∑[(1 + 𝛾). 𝑋 𝑖 . 𝑋 𝑖+1 + (1 -𝛾). 𝑌 𝑖 . 𝑌 𝑖+1 ] 𝑖=𝑁 𝑖=1
can be expressed as a quadratic form in creation and annihilation operators and can be diagonalized to full describe the complete set of states, excitation energies.

It is very similar to the generalized Heisenberg model described by the Hamiltonian:

𝐻(𝛾) = ∑[(1 + 𝛾). 𝑋 𝑖 . 𝑋 𝑖+1 + (1 -𝛾). 𝑌 𝑖 . 𝑌 𝑖+1 + (1 -𝛾). 𝑍 𝑖 . 𝑍 𝑖+1 ] 𝑖=𝑁 𝑖=1
Note that when 𝛾 → 1, this Hamiltonian tends to the Ising model where the 𝑥-componants are fully ordered.

The transverse terms (i.e. 𝑌 𝑖 . 𝑌 𝑖+1 and 𝑍 𝑖 . 𝑍 𝑖+1 ) consist in favoring the ordered of both 𝑦-components and 𝑧components. A complete description of this operator has been introduced in antiferromagnetic chain by [START_REF] Lieb | Two Soluble Models of an Antiferromagnetic Chain[END_REF].

Qudit operators and basic principles

Let us denote |𝜑⟩ = |𝜑 1 𝜑 2 𝜑 3 ⟩ a qubit composed of 3 qubits |𝜑 1 ⟩ = ( 𝑎 𝑏 ), |𝜑 2 ⟩ = ( 𝑐 𝑑 ), |𝜑 3 ⟩ = ( 𝑒 𝑓 ). The operator 𝑋 1 𝑋 2 + 𝑌 1 𝑌 2 is the application of operator 𝑋 1 + 𝑌 1 to qubit 1 and 𝑋 2 + 𝑌 2 to qubit 2.

𝑋 + 𝑌 can be rewritten:

1 √2 (𝑋 + 𝑌) = 1 √2 [( 0 1 1 0 ) + ( 0 -𝑖 +𝑖 0 )] = ( 0 1 √2 (1 -𝑖) 1 √2 (1 + 𝑖) 0 )
For

( 1 √2 (1 -𝑖). 𝑏 1 √2 (1 + 𝑖). 𝑎 𝑐 𝑑 1 √2 (1 -𝑖). 𝑓 1 √2 (1 + 𝑖). 𝑒 )
If we successively apply 𝑋 𝑖 𝑋 𝑖+1 + 𝑌 𝑖 𝑌 𝑖+1 to a basis qubit, this will allow to superpose the current state with a new one where the 1 value is successively switched to the next position in the qubit.

To implement this operator, we will manipulate: 𝑒 -𝑖.𝑡.𝛾.𝑋 1 .𝑋 2 = cos(𝑡. 𝛾) . 𝐼𝑑 -𝑖. sin(𝑡. 𝛾) . 𝑋 1 . 𝑋 2 𝑒 -𝑖.𝑡.𝛾.𝑋 1 .𝑋 2 = cos(𝑡. 𝛾) . ( 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ) -𝑖. sin(𝑡. 𝛾) . ( 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 ) 𝑒 -𝑖.𝑡.

Fig. 3. Example of proper coloring

The graph coloring is the corner stone of many OR problems including the timetabling, scheduling and it has received a considerable amount of attention for its utility.

The classical formulation relies on binary variables with 𝑥 𝑖𝑘 = 1 if the color 𝑘 is assigned to node 𝑖. With such a representation, the solution of figure 5 can be modelized by the following bit-string 1000 -0100 -1000 -0001 -0010 Each group of 4 digits is a Dirac representation of the color number k, each 5 group of digits represent a vertex. Here 𝑐 1. = 1000, 𝑐 2. = 0100, 𝑐 3. = 1000, 𝑐 4. = 0001 and 𝑐 5. = 0010.

The objective of a max K-coloring dealing with proper and improper colorings is to maximize the number of adjacent nodes having different color between 1 to K The driver operator

𝑀𝑎𝑥 ∑ ∑(𝑥

𝐻 𝐷 = -∑ ∑ 𝑋 𝑖𝑘 𝐾 𝑘=1 𝑛 𝑖=1
is applied to |𝜓 0 ⟩ = 𝐻. |0⟩ ⨂𝑛 to define the initial fundamental ground state.

The qubits-string encoding a solution is 𝑛 × 𝐾 qubit long and the Hamiltonian mapping function encodes a total search space of 2 𝑛×𝐾 colorings of 𝐺 that encompasses both proper and improper colorings. For an operational research point of view, the mapping function is not an indirect coding function and it is strongly different than the Split procedure of the VRP or than the efficient Bierwith vector encoding for the Job-Shop.

The previous remarks push us to consider that the Driven Hamiltonian must have afforded restrictions of exploration to coloring with only one color per node, so that the promise to only switch over coloring that satisfies the constraint C1, enabling us to control which qubit-strings must be iteratively investigated.

A Mixer Hamiltonian definition

The mixing operator

𝐻 𝐷 = -∑ ∑(𝑋 𝑖𝑗 𝑋 𝑖𝑗+1 + 𝑌 𝑖𝑗 𝑌 𝑖𝑗+1 ) 𝐾 𝑗=1 𝑛 𝑖=1
permits to have Driven Hamiltonian that is a 1-to-1 mapping function i.e. an operator that maps one qubitstring on the desired subspace where the C1 constraint holds (figure 4).

Qudit-strings proper coloring

Coloring with one color per node

Assignements of color to node

Driven Hamiltonian

Fig. 4. Example of proper coloring

The next sections are dedicated to numerical experiments achieved on the simulator provided by IBM using the Qiskit library. The objective is to evaluate and to illustrate the efficiency for a Mixer Hamiltonian on few instances of graph coloring.

Experimental validation

With the following widget examples, the goal is to illustrate the bonus achieve by a Mixer Hamiltonian but it is not a performance analysis of quantum v.s. classical optimization since only tiny examples are tractable on quantum simulator. Of course, any classical metaheuristic will overfit the QAOA results. In a near future, if scaling in the number of (denoised) qubits and long enough coherent qubits time occurs, this evaluation will give insight on the way to model efficiency combinatorial optimization problems.

Numerical experiments with QAOA with 5 nodes / 6 edges and 4 colors

First experiment is achieved on the graph of figure 5 with the following set of parameters:  𝑝 = 3 (iterations of Hamiltonian and mixer Hamiltonian);  𝐴 = 100 is weight assigned to the quadratic term in the objective function (it can be lower by this value improve the readability of the numeric solutions) ;  𝛽 0 ⃗⃗⃗⃗⃗ = (0.0, … 0.0) and 𝛾 0 ⃗⃗⃗⃗ = (0.0, … 0.0) ;  Cobyla is the iterative numerical method used to minimize the estimated expectation value with 500 iterations;

 100 samplings are used to estimate 𝐶 𝑝 (𝛽 ⃗⃗⃗⃗⃗⃗ , 𝛾 ⃗⃗⃗⃗⃗⃗ ) at step 𝑖;

 100 samplings of |𝜓(𝛽 * ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ are achieved at the end of QAOA to compute |𝑠 * ⟩.

Results are display in table 1. A very low number of samplings at the end of the QAOA algorithm collects subset of solutions strongly smaller than the total number of solutions and avoid a costly inefficient total enumeration. The very small value 100 used permits to have the optimal solution and 100 samplings of |𝜓(𝛽 * ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ samplings is sufficient in this very specific run. Note that with very low number of samplings significant differences could appear in the estimation of the distribution and that the optimal solution could not appear. The overall trend of the probability distribution can be estimated with larger number of samplings and conclusively push into considering that the emphasis of probabilities has a global trend to the improve probability of solutions where the constraint C1 holds. Note that it pushes us into considering that QAOA has aggregated probability on solutions closed to the optimal ones. The experiment of figure 7 has been achieved with same parameters and 1000 samplings of |𝜓(𝛽 * ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ To identify the combinatoric problems more closely and the consequence with the objective to propose the adequate solutions, we may consider a graph 𝐺 with 5 nodes and 4 colors for which the total number of qubit-strings is 2 20 = 1 048 576. For one node, there is only 4 assignments of exactly one color, over the 2 4 = 16 possible assignments. The probability to have one color only assigned to each node is about A careful analysis of the table 2 shows that about 0.6% of the distribution is aggregated to proper or improper solution without penalty cost (> 100) that is about 60 times bigger than the previous value 0.098% which is the value of a uniform state distribution.

The fact that the optimal solutions (solutions with cost 0) have a lowest probability than solutions with cost 1 or 2, is basically linked to the limited number of optimal solutions.

The representation of figure 6 shows how the probabilities of each cost vary between the uniform distribution (in red) used at the beginning of QAOA and the distribution found at the end of QAOA (blue).

Higher probabilities are assigned to lower costs at the end of QAOA and lower probabilities are assigned to higher cost proving that QAOA has be successful in controlling the distribution over iterations to high quality states.

The modifications of the distribution shapes of figure 5 can be analyzed considering the between the previous curves (after/before) display in figure 6. This ratio is higher on small cost (that are related to high quality solutions) that reassures that the theoretical based affirmation of QAOA promises meet our experiments.

Fig. 6. Overall deformation of probabilities

Numerical experiments with 8 nodes graph

We consider a graph with 8 nodes and 11 arcs introduced in figure 7. The qubit-string space grows to 4 294 967 296 with 4 colors since length are 32 qubits long (maximum size available on the IBM Qiskit simulator). Note that the qubit-strings search space is much larger than the number of graph coloring solutions which does not exceed 4 8 = 65 536 for the 8 nodes/4 colors instances (0.0015% of the qubit string search space).

To put it in more concrete terms, there exists a probability of 0.0015% to have a qubit-string in the form 1000 0100 …0001 0010 where only one color is assigned to each node. With 8 colors the qubit-string space is about 167 777 216 strings and the quantum circuit requires only 24 qubits. The number of graph coloring solutions is 3 8 = 6 561 for the 8 nodes/3 colors i.e. 0.00391% of the qubit-string search space. We can note in table 3 that the very low value of 𝑝 and the small number of samplings affect the QAOA capacity in aggregating the distribution on the optimal and quasi-optimal solutions when the problem scale increases.

Resolution of the graph coloring problem with the 8-nodes graph and 4 colors has a modelization based on 32 qubits which is the larger circuit that can be simulated on the IBM Qiskit simulator. The huge memory consumption in the exponentiation of the operators leads to a large time consuming simulation. One hour has been experienced for the 4 colors instance, on an Intel(R) Xeon(R) at 3.40GHz with 128GO of memory.

Numerical experiments of the "Driven Hamiltonian": 5 nodes graph

The previous remarks push us into considering that the Driven Hamiltonian should afforded restrictions of exploration to coloring with one color only per node, so that the promise to switch over coloring only that satisfy the constraint C1, enabling us to control which qubit-strings have to be iteratively investigated.

QAOA based on the Driven Hamiltonian has been used with the same set of parameters of section 4.2 and the results of table 4 show that 100% of the distribution is now aggregated on the solution space with one color per node meeting the theoretical considerations. The comparison of QAOA with a Mixer Hamiltonian and a classical Hamiltonian proves the competitive advantages of the Mixer Hamiltonian v.s. the classical one and establishes a numerical validation of the theoretical considerations on the mapping functions. For example, the split algorithm uses such a trick by transforming one giant trip into a VRP solution allowing next definition of specific moves in the neighbors. Note that the Bierwith's vector (Bierwith, 1996) for the shop-shop is the corner store of the best ever published methods because it permits to investigate only the set of acyclic disjunctive graphs that model solution only avoiding exploration of all disjunctive graph. The best ever published method on the Job-Shop [START_REF] Nowicki | A Fast Taboo Search Algorithm for the Job Shop Problem[END_REF] takes advantages of this kind of exploration.

The differences between both distributions is strongly significant considering that the probability of a solution of cost 0 is 20 times higher with the Driven Hamiltonian (table 5).

The differences between the two distributions that are apparent in table 4 illustrate that the Mixer Hamiltonian has a substantial advantage by providing a strongly better aggregation on the optimal and quasi-optimal solutions as highlighted on figure 8. The Driven Hamiltonian offers competitive advantages since it possible to use lower value of 𝑝 as stressed in table 5. The global trend of probabilities remains very similar to the trend obtained with 𝑝 = 3 and there is no significant differences between the two probabilities distributions (figure 9) pushing into considering that more compact circuit (low depth circuit) could be used with Driven Hamiltonian. The previous remark must be put in the context of optimization of strongly large scale instances where quantum algorithms make sense. The Driven Hamiltonian enormously simplifies the optimal solution investigation, reducing it to a simple and well-defined feasible search-subspace investigation where a strong deformation of probabilities has been applied to a very short number of strongly high quality solutions.

Concluding remarks

This work addresses the question of how to implement Mixer Hamiltonians to restrict investigation on a specific subspace for the well-known graph coloring problem that is one of the corner-stone of the operation researcher community. We provided general considerations on Mixer Hamiltonians that are mapping functions defining an indirect representation which is a quantum equivalent of the well-known indirect representation scheme widely used in OR field. The numerical experiments bring into considering the large applicability and efficiency such approaches. To conclude we should emphasis that theoretical considerations and the capability in numerical experiment should push the researchers of the OR community to investigate this new field but keeping in mind that the current quantum computers are currently limited in the number of qubits.

2

  Fig. 1. QAOA principles

Fig. 2 .

 2 Fig. 2. Links between the weighted Hamiltonian and the expected cost.

  𝑓(𝑥) = ∑ 𝑐 𝑖 . 𝑥 𝑖 𝑛 𝑖=1 ∀𝑗 = 1. . 𝑚 𝑔 𝑗 (𝑥) ≥ 𝑏 𝑗 Despite the fact that it is possible to create a QAOA circuit for an objective function of order higher that quadratic, and for a large majority of classical O.R. quadratic expressions can be used: Hamiltonian construction, an initial state |𝑠⟩ = | + ⟩ ⨂𝑛 gives an initial Hamiltonian decomposition composed only of 𝐼 and 𝑋 Pauli gates applied to 𝐻 to |0⟩ ⨂𝑛 leading to one uniform distribution over the 2 𝑛 values of 𝑓.

  The problem can be rewritten into a minimization of a quadratic formulation by considering: 𝑆 𝑘 = ∑ 𝑥 𝑖𝑘 . 𝑛 (𝑖,𝑗)∈𝐸𝑥 𝑗𝑘 as the number of adjacent nodes with the same color  the weighted penalized term, 𝐴. ∑ (1 the eigenvalue 1 and 0, that model |1⟩. Reformulation of the function of our coloring problem with 𝑍 𝑗 variables leads to:

Fig. 5 .

 5 Fig. 5. Distribution of Cost Function values

Fig

  Fig. 7. A graph with 8 nodes, 11 edges to color with 4 colors

%Fig. 8 .

 8 Fig. 8. Comparison of the two probabilities distribution with 𝑝 = 3 for both Hamiltonians

Fig. 9 .

 9 Fig. 9. Comparison of the two probabilities distribution with 𝑝 = 2 for the Driven Hamiltonian

  e. minimizing or maximizing one objective function 𝑓(𝑥) that is assumed to act on 𝑛 -𝑏𝑖𝑡𝑠 strings 𝑥. QAOA is based on 𝑝 consecutive iterations of one Hamiltonian 𝐻 𝑃 cumulated with a driver Hamiltonian 𝐻 𝐷 , where this weighted sum of Hamiltonian terms varies in time. The Hamiltonian 𝐻 maps the function 𝑓 with 2 𝑛 eigenvalues that model the 2 𝑛 values of 𝑓. The optimal solution (i.e. the extremal value of 𝑓) is an eigenvalue of 𝐻 and 𝐻 is satisfying: 𝐻. |x⟩ = 𝑓(𝑥)|x⟩ Because a Hamiltonian is a Hermitian operator it has a spectral decomposition: 𝐻 = ∑ 𝑒 𝑖 𝑖 |𝑒 𝑖 ⟩⟨𝑒 𝑖 | where |𝑒 𝑖 ⟩ is the 𝑖 𝑡ℎ basis vector.

	The Hamiltonian is defined with Pauli operator-basis and takes advantages of the Pauli 𝑍 leading to
	expression
	𝐻 = 𝛼 0 . 𝐼 + 𝛼 1 . 𝑍 0 . 𝑍 1 + 𝛼 2 . 𝑍 0 . 𝑍 2 …
	where 𝛼 𝑖 are real numbers.

  If we apply the operator 𝑋 2 𝑋 3 + 𝑌 2 𝑌 3 , |𝜑⟩ is now

			1 √2	(1 -𝑖). 𝑏
			1 √2	(1 + 𝑖). 𝑎
			1 √2	(1 -𝑖). 𝑑
			1 √2	(1 + 𝑖). 𝑐 𝑒
			(	𝑓	)
		𝑎
		𝑏
	our qudit |𝜑⟩ =	𝑐 𝑑	, the operator 𝑋 1 𝑋 2 + 𝑌 1 𝑌 2 transforms |𝜑⟩ into a new state
		𝑒
	(	𝑓)

  If there exists (𝑖, 𝑗) ∈ 𝐸(𝐺), with 𝜑(𝑖) = 𝜑(𝑗), then 𝜑 is an improper coloring of 𝐺. Figure3gives a proper coloring solution for a 5 nodes graph.

				Color 2			Color 1
		Color 1	1	2				3
		Color 4		Color 3
			4	5
			cos(𝑡. 𝛾)	0		0	-𝑖. sin(𝑡. 𝛾)
	𝛾.𝑋 1 .𝑋 2 = (	0 0	cos(𝑡. 𝛾) -𝑖. sin(𝑡. 𝛾)	-𝑖. sin(𝑡. 𝛾) cos(𝑡. 𝛾)	0 0	)
		-𝑖. sin(𝑡. 𝛾)	0		0	cos(𝑡. 𝛾)
	This is easily converted to the following tensor product, using quantum gates CNOT (CXij), Hadamard (H)
	and Rotation on 𝑋 or on 𝑍 axis:
		𝐻 ⊗ 𝐻. 𝐶𝑋 1,2 . 𝑅 𝑍 2 (𝑡. 2. 𝛾). 𝐶𝑋 1,2 . 𝐻 ⊗ 𝐻 for 𝑋 1 𝑋 2 and
	𝑅 𝑋 1 ( 𝜋 2	) ⊗ 𝑅 𝑋 2 ( 𝜋 2	) . 𝐶𝑋 1,2 . 𝑅 𝑍 2 (𝑡. 2. 𝛾). 𝐶𝑋 1,2 . 𝑅 𝑋 1 (-	𝜋 2	) ⊗ 𝑅 𝑋 2 (-	𝜋 2	) for 𝑌 1 𝑌 2 .
	4. Graph coloring and max k-coloring

4.1 Definition

Coloring problem want to find the minimum value of the chromatic number 𝜒(𝐺) where there exists a 𝜒(𝐺) coloring for a graph 𝐺. More formally, for a graph 𝐺(𝑉, 𝐸), a proper coloring 𝜑: 𝑉(𝐺) ⟶ 𝐶 exists with 𝜑(𝑖) ≠ 𝜑(𝑗) if (𝑖, 𝑗) ∈ 𝐸(𝐺).

Table 1

 1 One possible run with 𝐴 = 100 and 100 samplings of |𝜓(𝛽

	Solution cost Probability
	0	1
	101	1
	102	3
	103	3
	201	7
	202	4
	203	4
	204	1
	300	1
	301	1
	…	…

* ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩

Table 4

 4 Sampling of sampling of |𝜓(𝛽Driven Hamiltonian performances meet the classical well know results in Operation Research community where significant improvements in results have been obtained by both dedicated local search that restrict neighborhood to feasible solutions and by metaheuristic based approaches where the search space investigation is driven by specific information on the solution characteristics.

		Probabilities with the
		classical mixer operator	Mixer Hamiltonian
	Solution cost	𝒑 = 𝟑	𝒑 = 𝟑
	0	0.2	4.8
	1	0.3	26.9
	2	0.1	34.6
	3	0.0	22.6
	4	0.0	6.8
	5	0.0	2.8
	6	0.0	1.4
	100	0.3	0.0
	101	1.6	0.0
	102	1.3	0.0
	103	0.9	0.0
	104	0.2	0.0
	106	0.1	0.0
	…	…	…

* ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ with the two Hamiltonians

Table 5

 5 Sampling of sampling of |𝜓(𝛽

		Probabilities with the	Probabilities with the
		classical mixer operator	Driven Hamiltonian
	Solution cost	𝒑 = 𝟑	𝒑 = 𝟐
	0	4.8	3.2
	1	26.9	26.7
	2	34.6	37.1
	3	22.6	21.8
	4	6.8	5.4
	5	2.8	3.0
	6	1.4	2.6

* ⃗⃗⃗⃗⃗ , 𝛾 * ⃗⃗⃗⃗⃗ )⟩ with the two Hamiltonians