Multi-Objective Optimization and Network Routing with Near-Term Quantum Computers

Multi-objective optimization is a ubiquitous problem that arises naturally in many scientific and industrial areas. Network routing optimization with multi-objective performance demands falls into this problem class, and finding good quality solutions at large scales is generally challenging. In this work, we develop a scheme with which near-term quantum computers can be applied to solve multi-objective combinatorial optimization problems. We study the application of this scheme to the network routing problem in detail, by first mapping it to the multi-objective shortest path problem. Focusing on an implementation based on the quantum approximate optimization algorithm (QAOA) -- the go-to approach for tackling optimization problems on near-term quantum computers -- we examine the Pareto plot that results from the scheme, and qualitatively analyze its ability to produce Pareto-optimal solutions. We further provide theoretical and numerical scaling analyses of the resource requirements and performance of QAOA, and identify key challenges associated with this approach. Finally, through Amazon Braket we execute small-scale implementations of our scheme on the IonQ Harmony 11-qubit quantum computer.

Multi-objective optimization problems (MOOPs) arise naturally in many scientific and industrial areas, where the interplay between multiple conflicting objectives gives rise to a set of optimal solutions, rather than a unique one.They are especially prevalent in engineering contexts, where complex systems involving multiple objectives are encountered [1][2][3].In particular, we focus on the case where the state space and hence the set of optimal solutions is discrete, where they are referred to as multi-objective combinatorial optimization problems (MCOPs).An example of such a MCOP is the multi-objective network routing problem, which asks for paths between specified source and destination nodes in a graph that are optimal with respect to multiple objectives.This optimization problem arises in contexts such as the design of Wireless Ad-Hoc/Sensor Networks [4][5][6][7] and next-generation communication networks [8][9][10], where large-scale networks with multiple requirements have to be simultaneously satisfied to meet performance demands.
However, the solution of such problems is generally difficult, as finding the global optima -also known as Pareto-optimal solutions -of general MCOPs is NP-hard [11,12].Physics-and biology-inspired meta-heuristic classical algorithms have been particularly successful in this area in the past two decades [13][14][15] due to their scalable computational costs, and are the go-to methods for industrial use-cases.Even so, due to the importance and complexity of MCOPs, there is a need to develop theoretical and algorithmic tools to solve large-scale problems in more memory and time efficient ways.
Besides developments in classical computers, quantum computers are currently experiencing an explosive growth, both in theory and realization.Noisy intermediate-scale quantum (NISQ) computers with qubit counts in the small hundreds are now available, and there is a growing body of work investigating their ability to solve challenging optimization problems in areas ranging from chemistry [16] to finance [17], with the hope that they can outperform existing classical algorithms in the near future.In particular, there has been a strong focus on variational quantum algorithms (VQAs) such as the Quantum Approximate Optimization Algorithm (QAOA) [18], which can be used to tackle combinatorial optimization problems with existing quantum computers.A natural question is whether similar methods can be applied to to efficiently obtain high quality solutions to MCOPs such as the network routing problem, to satisfy the performance demands of next-generation wireless networks.
In this context, a significant challenge is the fact that, for a given use-case, investigating the performance of VQAs such as QAOA typically requires an empirical approach.However, given the relatively small size of currently available quantum computers and the complexity of simulating quantum computations classically, the scope for experimentation is currently rather limited.Nonetheless, in this work we adopt a pragmatic approach, using the standard form of QAOA [18] as a solution method for MCOPs.This permits us to: (a) use small-scale problem instances to draw insights on how the QAOA solution relates to Pareto-optimal solutions of the MCOP, and; (b) analyse the scalability of the QAOA approach against known generic limitations of VQAs.
The key results of our work are summarized as follows: 1. We develop a general framework with which nearterm quantum computers can be used to solve MCOPs, by producing multiple Pareto-optimal solutions in both a priori (with the preferences of the decision maker taken into account prior to optimization) and a posteriori (independent of the decision maker's preferences) manners.This is achieved by casting the multiple objectives and constraints of a MCOP in Quadratic Unconstrained Binary Optimization (QUBO) form, and scalarizing it to obtain a cost function which can be variationally optimized with a VQA.For small problem instances, visualization of the output quantum state on Pareto plots provides insights which we analyze and explain in a qualitative manner.
2. We provide a formulation of the network routing problem that is amenable to implementation with QAOA on near-term quantum computers.Using results from graph theory, we determine its resource requirements, and provide a scaling analysis.In particular, we show that this encoding scheme possesses resource requirements that scale mildly with the connectivity of the underlying graph, which is in principle compatible with resources available on current NISQ hardware.
3. Numerical simulations of small problem instances using standard QAOA (as described in Ref. [18]) show that this framework can produce high-quality solutions efficiently.Concretely, by increasing the circuit depth, we observe a correspondingly proportional increase in the probability of successfully obtaining Pareto-optimal solutions.However, as we explain in Section IV A, in terms of problem size the efficacy of the QUBO-based QAOA approach is limited by the large fraction of infeasible solutions in the underlying search space, resulting from the presence of optimization constraints.
4. We run a number of demonstrative test cases on the 11-qubit IonQ Harmony quantum computer, accessed through Amazon Braket.The results obtained are in clear agreement with those of numerical simulations.
This article is structured as follows.We begin by fully describing our approach to obtain optimal solutions of a MCOP and discuss interpretations of resulting Pareto plots in Section.II.We then describe the network routing problem, a concrete example of a generally difficult MCOP which we study throughout the article, along with relevant objectives that arise from reasonable assumptions in Section.III.This is followed by theoretical and numerical analyses on the scaling of resources and performance of the procedure in Section.IV, and experiments on quantum computers in Section.V.
We refer the reader to the appendices for a survey of relevant work in the usage of near-term quantum algorithms in solving MCOPs and the shortest path problem (Appendix.A), a review of MOOPs (Appendix.B), the mathematical formulation of the multiobjective shortest-path problem (Appendix.C), graph theoretic ideas used in this work (Appendix.D), the formulation of the network routing problem along with complete specifications of physical parameters used to model its objectives (Appendix.E), specifications of problem instances considered in our numerical calculations and experiments (Appendix.F and H), additional numerical results for the chosen parameter initialization scheme (Appendix.G), and additional discussions on the Pareto plot (Appendix.I).
Numerical simulations and the set-up of quantum computations in this work were performed through Open-QAOA [19], an open-source Python package tailored for QAOA and its variants.

II. SOLVING MCOPS WITH QUANTUM COMPUTERS
In this section, we describe our approach to solve MCOPs using near-term quantum computers with VQAs such as the QAOA, starting with brief introductions to the QAOA, MCOPs, and combinatorial optimization/QUBO problems.For more in-depth discussions on the above topics, we refer the reader to [18], [20] and [2,3] respectively, along with the Appendices.
A key component of our procedure is the QAOA.Belonging to the class of VQAs, it was introduced to allow NISQ computers to provide approximate solutions to combinatorial optimization problems such as graph partitioning, coloring, and constraint satisfaction problems [18].A QUBO problem can be defined by an upper triangular matrix with real values.It is specified within QAOA through a corresponding cost Hamiltonian H c expressed as a sum of l local terms H c = l i=1 h i , which can be obtained by converting the binary variables of the QUBO problem to Ising variables [20].This allows the quantum state |ψ( ⃗ β, ⃗ γ)⟩ to be prepared by a p-layer state-preparation ansatz U QAOA ( ⃗ β, ⃗ γ): where the real vectors ⃗ β = (β 1 , ..., β p ) and ⃗ γ = (γ 1 , ..., γ p ) are variational parameters and |+ • • • +⟩ is the initial quantum state.By sampling from a quantum computer, the expectation values ⟨ψ( ⃗ β, ⃗ γ)| H c |ψ( ⃗ β, ⃗ γ)⟩ can be computed and taken as the cost function for the optimization problem: In practice, the optimization problem is solved in a variational manner with the usage of a classical optimization algorithm, yielding a quantum state |ψ( ⃗ β * , ⃗ γ * )⟩ that approximates the ground state of H c .This in turn allows the solution of the underlying QUBO problem to be extracted.Increasing the number of layers p leads to an increase in the expressibility of the ansatz, and hence the potential quality of the output solution, at the expense of longer computation times [21].We have chosen to investigate the original form of QAOA described in Ref. [18] in our work, motivated by its relative simplicity and amenability to implementation on currently available quantum computers, with a few tens of qubits.
Next, we provide a general description of MCOPs.Given L objective functions that map states from a state space to real numbers, multi-objective optimization asks for states that are optimal with respect to all of the L objectives.Since the objectives generally produce competing effects with one another, states that are simultaneously optimal in all objectives generally do not exist.In this case, we ask instead for a set of Paretooptimal/efficient states, which are optimal in the sense that no other states which improve on at least one individual objective without deteriorating in others can be found.On top of the objectives, a number of additional constraints may also be present, imposing further complications.
More concretely, we require that the L individual objectives and K constraints of the problem can be cast into a QUBO problem.This yields an encoding of the states x = (x 1 , ..., x n ) as vectors of n binary decision variables x 1 , ..., x n ∈ {0, 1} in a state space S = {0, 1} n , L quadratic objective functions C i : S → R of the form: and K quadratic penalty functions P j : S → R that enforce the K constraints: where the Q i 's and P j 's are upper triangular matrices with real entries.To compute these functions with quantum computers, we write the objective and penalty functions in terms of Ising variables s = (s 1 , ..., , to eventually convert them into equivalent Hamiltonians, and denote the corresponding objective and penalty functions in terms of s as E C i (s) and E P j (s) respectively.The MCOP is then: subject to the minimization of E P 1 (s), ..., E P K (s).Numerous classical methods exist to solve the above optimization problem.Of particular interest to us is linear scalarization, which seeks for a solution of Eq. ( 5) by first solving a simpler problem obtained by aggregating the objective and penalty functions of the full MCOP in a linear manner: where the scalarization weights w i ≥ 0, and the penalty weight w P > 0 controls the degree at which infeasible solutions are penalized.A choice of w P that is sufficiently large ensures that the solution of Eq. ( 6) is necessarily a Pareto-optimal solution of Eq. ( 5) that satisfies the P constraints imposed (i.e. a feasible solution) [3,22,23].In practice, smaller values can be used, and we set w P = 1 throughout this work in an empirical manner, and set w i ∈ [0, 1] without loss of generality (since only the relative weight w P /w i matters).
For further details and references on MCOPs and scalarization, we refer the reader to Appendix.B.
A. Obtaining the Pareto Front with the QAOA We now describe our procedure for solving MCOPs with near-term quantum computers, which is summarized in Fig. (1).
Following Eq. ( 6), the QUBO objective vector ⃗ E C (s) = (E C 1 (s), ..., E C L (s)) is first scalarized by aggregating them as a convex sum with weights ⃗ w = (w 1 , ..., w L ) with w i ∈ [0, 1] and L i=1 w i = 1.This is added to the penalty terms K i=1 E P i (s) with each penalty function weighted equally by 1, resulting in a scalarized QUBO cost function: Minimization of this cost function then yields a feasible solution that corresponds to a Pareto-optimal solution s * of the MCOP Eq. ( 5).This optimization can be performed on a quantum computer by constructing equivalent objective Hamiltonians H C 1 , ..., H C L and penalty Hamiltonians H P 1 , ..., H P K via transformation to Ising variables, and aggregating them to result in an anologous scalarized cost Hamiltonian: whose ground state encodes the same solution s * , where ⃗ H C = (H C 1 , ..., H C L ).We perform the ground state search with a VQA, chosen to be QAOA in our work due to its simplicity and amenability by current quantum computers.That is, we variationally optimize the QAOA circuit of Eq. (1) to produce an optimized output quantum state |ψ( ⃗ β * , ⃗ γ * )⟩ that approximates the ground state of H scalar , which is then sampled from k times to yield a set of candidate solutions B = {b 1 , ..., b k }, with k always chosen to be at most polynomial in the number of qubits n.The set B then constitutes a reduced search space within which Pareto-optimality can be efficiently checked with any classical method (such as brute force).The output of this procedure is therefore expected to be a set of feasible solutions that are Pareto-optimal within B. As the number of QAOA layers is increased to ∞, the QAOA ansatz supports an increasingly better approximation of the ground state of Eq. 8. Provided that it can be found in the variational procedure, sampling this state would then return solutions that are also Pareto-optimal within S, with high probability.
The above scalarization procedure involves the linear aggregation of objectives.The convex weights w i can be interpreted as a priori preferences that the decision maker can select, which biases the resulting solution according to ⃗ w.While other scalarization schemes that involve higher order terms (such as quadratic scalarization) and inequality constraints (such as Chebyshev scalarization) are available, we focus on the linear scheme, as it preserves the QUBO form of the individual objectives and constraints.The resulting scalarized cost function Eq. ( 7) then contains only linear and quadratic terms.We discuss possible limitations and extensions of this approach in Appendix.B.
So far, this procedure belongs to the class of a priori methods, where the decision makers' preferences are specified before the optimization to bias the output of the optimization procedure (c.f.Appendix.B for details).To recover the entire Pareto front in an unbiased way, this procedure can be converted into a completely a posteriori method by repeating the procedure with numerous choices of ⃗ w in a problem-agnostic manner, with e.g.uniform random sampling of the weights or a discretization over all possible weights.The full procedure to recover the Pareto front with QAOA is illustrated in Fig.

B. Visualization of Solutions in the Pareto Plot
The state space of an L−objective MCOP can be visualized in a Pareto plot [2,3], which is an L−dimensional The graph (with solid lines denoting its edges and circles denoting its nodes) is converted to the Hamiltonian Eq. ( 8) that incorporates the objective and constraints of the problem.The Hamiltonian itself can be represented by a graph, G H (bottom right graph, with dotted lines with dotted lines denoting edges and circles denoting qubits).The illustration describes the conversion for the multi-objective routing problem, which is elaborated in Section.III.(b) Quantum-classical optimization with QAOA to solve the optimization problem Eq. ( 2) for the scalarized cost Hamiltonian Eq. ( 8).The p-layer QAOA circuit is run at different angles according to the feedback of a classical optimizer, navigating the cost landscape until it converges to a minima.(c) Aggregation over results of one or more scalarization choices, determination of Pareto front with classical solver, and extraction of final result(s).
plot with 2 |S| points.Each point corresponds to a state/bitstring s, with the value of the L objectives as its coordinates.Pareto-optimal solutions then constitute points located at the boundary of the region populated by the 2 |S| points [2].
In our case, we must also account for the penalty terms E P 1 (s), ..., E P K (s), whose contribution should be reflected in the Pareto plot in a consistent way.This can be achieved by defining a point ⃗ r = (r 1 , ..., r L ) in the Ldimensional Pareto plot as follows: This definition has the property that the cost function arising from linear scalarization with weights ⃗ w can be interpreted as the projection of points in the Pareto plot onto the vector ⃗ w, since: which is exactly the scalarized cost function Eq. ( 7) that is minimized in our procedure.We now give a brief qualitative discussion of some notable features of the Pareto plot obtained from the application of QAOA to a MCOP, chosen to be a small network routing problem (to be defined in Section.III) defined on 13 qubits for illustration.While the precise details of a Pareto plot will naturally differ depending on the specific problem instance, the observations we make in the following were general to all the cases we have considered in this work.2)), the sampled states can be observed to be squeezed vertically, towards the bottom of the Pareto plot.
• Points in the Pareto front corresponding to solutions of other scalarization choices can be sampled : If there are other Pareto-optimal states close to the solution of the scalarized problem in the Pareto plot, in the sense that they are low-energy states of the scalarized problem, the set of candidate states B i sampled from |ψ( ⃗ β * , ⃗ γ * )⟩ may contain these nearby Pareto-optimal solutions with relatively high probability.In principle, this allows multiple Pareto-optimal states to be obtained from one iteration of the procedure.For example, in Fig. (2) we observe that the Pareto-optimal solution corresponding to the rightmost red point is sampled with relatively high probability in all three scalarization choices, even though it is not a solution of any of these scalarized Hamiltonians.We observe that this is often the case in problem instances we consider.
• Points in the Pareto front located in locally concave regions can be sampled : Again, because states near the target Pareto optimal point can be sampled with high probability from |ψ( ⃗ β * , ⃗ γ * )⟩, states located at locally concave regions of the Pareto front can in principle be sampled as long as they are near the target (in the Pareto/objective space), even with linear scalarization.This allows Pareto-optimal points located at concave regions to be identified, even though they can never correspond to solutions of scalarized cost functions.Using the rightmost red point in Fig. (2) as an example again, we observe that it is sampled with relatively high probability, despite being located at a locally concave region, which cannot be obtained as the solution of any scalarized Hamiltonian.
The latter two observations indicate a degree of robustness to the choice of scalarization and the concavity of the Pareto region, due to QAOA being a sampling-based algorithm, which requires sampling from the optimized quantum state in order to extract information it.This introduces a trade-off between the ability of QAOA to approximate the ground state (by increasing p), and for |ψ( ⃗ β * , ⃗ γ * )⟩ to remain as a superposition of energy eigenstates so that low-energy states can be sampled (by keeping p small), potentially imposing a less stringent scaling of p with problem size.
Finally, we remark that most points in the Pareto plots correspond to infeasible solutions, which violate at least one of the problem constraints encoded in the penalty terms of Eq. ( 8).As we explain in Appendix.I, in general the fraction of feasible to infeasible solutions decreases rapidly with the network size.This fact represents a challenge for QUBO approaches where the encoding of constraints as penalty terms is unable to restrict the search through solutions to remain solely within the feasible subspace.We discuss this point further in Section.IV A in the specific context of QAOA.

III. APPLICATION TO THE NETWORK ROUTING PROBLEM
With a general description of our algorithm for MCOPs established, we now consider the concrete example of the multi-objective network routing problem.
We consider a generic multi-hop wireless network with relay stations distributed throughout a geographical area.This can be modelled as an undirected weighted graph G = (V (G), E(G)) (which we also refer to as the network's graph) with nodes V (G) that correspond to relay stations and edges E(G) defined by possible transmission paths between stations (Directed graphs can be included with a more general formulation [22]).Two nodes are specified to be the source and destination nodes (denoted s and d respectively).The routing problem then asks for data transmission routes between s and d that are Pareto optimal with respect to multiple objectives, which can be written as functions of the node and edge weights of the graph.This is an instance of the multiobjective shortest path problem, which is a combinatorial optimization problem known to be NP-hard [24].Its solution is relevant across a wide range of network design and optimization tasks [4][5][6][8][9][10], especially beyond current wireless network protocols based on individual objectives such as the Optimized Link State Routing protocol [25] which only considers hop count, and is thus unable to maximize network resource utilization.
The form of the objectives depends on specific performance demand requirements.We consider the following four objectives in our work: 1. Path loss: Transmission between two stations incurs an energy cost that is dependent on their physical distance.
2. Node Delay: Signal processing at each station incurs a time delay, which sums up to an overall delay.
3. Data rate: The total data output rate of a transmission path is determined by the minimum data rate along its path.
4. Bit error: Bit errors occur with finite probability during transmission between two stations, which depends on the path loss and the channel's noise profile.
To solve this problem with our approach, a QUBO encoding of the objectives and constraints of the multiobjective shortest-path problem is needed.We provide a complete description of this encoding in terms of binary variables together with the form of the quadratic cost function in Appendix.C (following [22]).This is followed by QUBO formulations of the above four objectives in Appendix.E. Along with reasonable assumptions on network parameters based on software-defined radio use cases [26], the associated objective and penalty Hamiltonians H C i and H P i can then be constructed explicitly, thereby fully specifying the inputs for our procedure in Fig. (1).
Importantly, we observe that for a network with graph G, the graph of its corresponding scalarized Hamiltonian G H -that is, the graph with connectivity defined by the quadratic terms of Eq. ( 8) -can be interpreted as the network graph's middle graph [27], denoted as M (G).Using known properties of M (G), we will exploit this correspondence in the next sections to bound the resources needed to execute the algorithm.The Hamiltonian's construction is illustrated in Fig. (1).(a).The network's graph G (top left graph, with solid lines denoting its connectivity, red circles denoting source and destination nodes, and blue circles denoting remaining nodes) is converted to a Hamiltonian H, with connectivity defined by the graph G H , which takes the form of M (G) (bottom right graph, with dotted lines denoting the connectivity of G H or M (G), and circles denoting nodes).We elaborate on points regarding M (G), including its definition and properties, in Appendix.D.

IV. THEORETICAL AND NUMERICAL SCALING ANALYSES
We now turn to analyze in detail the scalability and performance of the QAOA approach to the network routing problem.
In Section IV A, we begin with a short qualitative discussion of the role of the QAOA mixer and initial state in determining the efficacy of the algorithm.Subsequently, in Section IV B, we discuss and provide a pragmatic analysis of the hardware resources required to implement a problem instance of given size, and compare these requirements against known limitations of VQAs.This allows us to distinguish between problem instances that are infeasible for our approach, from those that are potentially amenable.Finally, in Section IV C we present the results of numerical simulations that explore the scaling of several success metrics as a function of problem size and the QAOA depth p.We remark that due to limitations on the system sizes that can be simulated on a classical computer, strong conclusions on the asymptotic performance of the algorithm cannot be drawn at this stage.
All analyses performed here are based on a single scalarization (i.e.we apply the algorithm in an a priori manner by pre-specifying the scalarization weights ⃗ w).As described in Section.II A, if the intent is to recover the entire Pareto front instead, the total runtime of the procedure depends on the number of scalarizations, which introduces an additional multiplicative overhead in the time required to obtain the Pareto front.As is also the case with classical approaches, this additional complexity is highly problem-dependent.

A. QAOA mixer and initial state
We begin with a qualitative discussion of the role of the QAOA mixer and circuit initial state in determining the efficacy of the approach.In this work, we have used the standard QAOA initial state and mixer pair [18], where the circuit is initially prepared in an equal superposition of all solutions, i.e. |+⟩ ⊗n , and the mixer Hamiltonian H m = − n i X i drives bit flips across the register.Constraints are enforced through energy penalties in the cost Hamiltonian Eq. ( 8), with the search ideally converging towards low-energy feasible solutions.However, the solution space is increasingly dominated by infeasible configurations as the network size grows (see Appendix.I), motivating the need for alternative ways of searching the solution space.
In the context of the network routing problem, we leave the question of designing improved initial state and mixer Hamiltonian pairings for future work (see also our remarks in Section.VI).However, we note that the goal of such strategies is to reduce the size of the solution space to be searched, either by entirely avoiding infeasible solutions, or by avoiding some subset of them.In the former case, where the search takes place through feasible solutions only, the penalty terms in Eq. ( 8) can be eliminated.In the latter case (which may arise if a strategy to search only feasible solutions cannot be found, or carries impractical resource requirements) the penalty terms would still be necessary to enforce constraints indirectly (i.e. through the objective function) With these considerations in mind, in the following subsections we give a detailed analysis of the resource requirements of Eq. (8).We also remark that in the standard QAOA approach (Eq.( 1)), the initial state |+⟩ ⊗n can be easily prepared by applying a Hadamard gate to each qubit.Since the mixer Hamiltonian e −iβi n j=1 Xj also involves only single-qubit operations, its implementation is also easy.

B. Cost Hamiltonian: Theoretical Resource Estimations
To determine the resource requirements of our QAOA approach, we examine the runtime and number of qubits necessary for its implementation on quantum hardware.This is determined by considering the resources needed to implement each layer of the QAOA circuit, and the number of repeated executions required to estimate the cost function up to an error ϵ.
We remark that for a given error tolerance ϵ, a complete performance analysis would involve exposing the dependence of the number of required QAOA layers p on problem size, which is challenging from both an analytic and numerical perspective.This is further related to the issue of trainability and existence of barren plateaus of VQAs [28][29][30], which depend strongly on problem class, presence of noise, and ansatz choice, and is a subject of ongoing study.Specifically, if the physical circuit depth for an application scales super-linearly in n, the optimization procedure will suffer from a noise-induced barren plateau [28], implying expensive gradient computations that scale exponentially in n, a phenomenon which is conceptually similar to the issue of vanishing gradients that previously plagued the training of neural networks in classical machine learning [31].
We begin by outlining a few assumptions and simplifications that we will make.Firstly, the quantum circuit of a p-layer QAOA consists of p alternating mixer and cost unitaries (Eq.( 1)).Since elementary gates in the cost unitary e −iγiHc fully commute with one another, the 1qubit RZ gates can be scheduled to be executed first, in parallel.Together with the fact that the mixer unitary consists only of 1-qubit RX gates which can also be executed simultaneously, and that 1-qubit gates execute significantly faster than 2-qubit gates, we only need to consider contributions from 2-qubit gates.Secondly, we neglect compilation overheads arising from qubit routing, i.e. the need to include additional SWAP gates to carry out 2-qubit gates that are not natively executable on a quantum processor.While this is a valid assumption for quantum processors with a fully-connected topology (such as currently available ion trap devices, which we utilise in Section.V), the mismatch between the quantum processor and the cost Hamiltonian's topologies will generally require routing.The routed output depends on the degree of mismatch and the routing strategy employed, among many other factors which are beyond the scope of our discussion.Nonetheless, we remark that the inclusion of n − 1 additional layers of SWAP gates is a naive upper bound [32,33], incurring at most O(n) additional layers of quantum gates.Thirdly, to streamline our argument we ignore quadratic terms arising from the source and destinations constraints, which only incurs a negligible, constant number of terms (detailed in Appendix.C).
Finally, we remark that in our formulation of the problem, increasing the number of objectives solely results in additional linear terms in the cost Hamiltonian Eq. ( 8), since only the penalties contain quadratic terms.An increase in the number of objectives therefore does not directly incur additional time and qubit costs.

Number of qubits
Since the graph of the cost Hamiltonian G H can be mapped to the middle graph M (G), the number of qubits involved is equal to the number of nodes of M (G), re- , where ∆ G is the maximum degree of the network's graph G (see Appendix.D, Eq. (D5)).This implies that problem instances with maximum degree independent of problem size are expected to be more efficient, with n = O(|V (G)|).

Cost computation time
The computation time of the p-layer QAOA cost function is determined by two factors: (1) the number of repetitions n rep required to estimate the cost function expectation value, given some specified error tolerance, and; (2) the depth of each executed circuit, which we denote D p for p layers.These two factors combine to give an execution time of O(D p n rep ), or O(pD 1 n rep ) in terms of the depth of a single QAOA layer D 1 .
D 1 depends on the number and connectivity of the 2-qubit terms in the cost Hamiltonian.As they are Ising gates that are diagonal in the computational basis, they fully commute with one another, and can be scheduled to maximize parallelization.The determination of such a schedule amounts to the solution of an edge-coloring problem on G H , where edges (i.e.gates) of the same color are scheduled to be executed simultaneously.The minimal number of colors required is commonly termed the edge chromatic number χ ′ (G H ), so we have D 1 = χ ′ (G H ) layers in total.To estimate χ ′ (G H ), note that the structure of G H is completely determined by Eqs.(C3)-(C5).Intuitively, we expect G H to inherit the local structure of G, since the squares in Eqs.(C3)-(C5) only lead to interactions between nodes (of G H ) representing nodes and edges incident to them (i.e. the x i x ij terms), and between edges that are incident on the same node (i.e. the x ij x ik terms).It turns out that this is precisely the relationship between a graph G and its middle graph M (G) -in other words, G H = M (G), and hence χ ′ (G H ) = χ ′ (M (G)).This realization allows us to leverage on known properties of M (G) to show that: where we delegate the proof and a review of relevant ideas to Appendix.D. We conclude that the depth D 1 of each QAOA layer scales linearly with the maximum degree of the network's graph G, i.e.D 1 = O(∆ G ). On the other hand, it is known that n rep = O(L/ϵ 2 ) measurements are needed to compute expectation values up to an error tolerance of ϵ by operator averaging [34], where L is equal to the number of Pauli terms in the operator (which is the cost Hamiltonian in our case).Of the L terms, terms are linear (since the number of node and edge weights is equal to the total number of nodes and edges of G), while G /ϵ 2 ).Summarizing, the computation time of a QAOA cost function up to an error ϵ scales as O(p|V (G)|∆ 3  G /ϵ 2 ), which can be parallelized up to a factor of n rep .
We remark that this implies efficiency for problem classes with a connectivity or maximum degree that is independent of the problem size (such as k-regular/latticetype graphs), with only a linear dependence of the computation time in problem size.This is likely the sce-nario for large-scale applications, since the cost of supporting a dense network (such as fully-connected mesh networks) over large geographical regions will otherwise be high, rendering it practically infeasible.In that case, the computation time (and the number of qubits needed) becomes at most quadratic in problem size.

Performance Metrics
We now introduce several measures that quantify the degree to which our scheme has successfully solved the task at hand, which is to obtain Pareto-optimal points of a multi-objective problem.
A standard measure of the quality of solution output by QAOA is the approximation ratio, defined as the ratio between the energy of the output state and the ground state: An approximation ratio of 1 thus implies that QAOA has found the exact ground state of H.
In our multi-objective context, however, we ask for the set of Pareto-optimal solutions, which generally cannot be encoded as ground states of a single Hamiltonian.As a more pragmatic measure in a multi-objective context, we supplement analyses of r approx by defining the success probability, which is the probability that a Pareto-optimal state is sampled from |ψ( ⃗ β * , ⃗ γ * )⟩, A success probability of 1 implies that sampling from |ψ( ⃗ β * , ⃗ γ * )⟩ will always yield a Pareto-optimal state.Note that it does not contain information about whether it contains points from the entire Pareto front.

Scaling with System Size and Number of Layers
In this section, we present results from classical simulations of QAOA to study the scalability of the approach with respect to different problem sizes, graph geometries, and the number of QAOA layers p.These simulations assume an absence of both statistical and hardware noise.Taking into account limitations on the classical simulation and optimization of general QAOA circuits for large n and p, and the need to average over multiple problem instances, we consider problems of sizes up to 16 qubits and p = 10 layers.
Fig. (3) displays the scaling of the approximation ratio and success probability as a function of p, for different problem sizes chosen from different graph geome-tries, namely triangular lattices, square lattices, and cycle graphs.For all examples, we consider two linear objectives, one involving only edge variables (i.e. of the type Eq.(C2)) and another only involving node variables (i.e. of the type Eq.(C1)), with node/edge weights uniformly randomized between -1 and 1, averaged over 50 instances.The scalarization weights are also chosen to be equal for both objectives, i.e ⃗ w = (0.5, 0.5).Finally, initial parameters for QAOA are chosen according to a linear ramp initialization scheme, which is a heuristic choice that linearly ramps up γ's and ramps down β's, based on the analogy between QAOA and quantum annealing [21].Additional numerical results for the performance of this initialization scheme can be found in Appendix.(G).
As first observations, we see that both metrics increase monotonically as a function of p.This is in accordance with the fact that as p → ∞ with QAOA, we recover the limit of infinitesimally slow adiabatic quantum annealing, and recover the ground state/optimal solution of scalarized cost Hamiltonians with unit probability.For any fixed value of p, both metrics also achieve higher values for smaller problems, again mirroring the increase in annealing time needed for larger problems.Finally, we observe that problem classes with higher average connectivities tend to reach saturation at a slower rate (number of layers to reach saturation, in descending order : cycle graphs, triangular lattice, square lattice).These conclusions hold for generic problem instances of the same class, due to the average over numerous instances.
Focusing on the success probability (lower row of Fig. (3)), we find that it increases monotonically with p, with indications of saturation at the maximum value of 1 for the limited problem sizes we have considered.These results suggest that, at least for the problem instances considered here, the procedure indeed outputs Pareto-optimal solutions with increasing frequency as we increase p.This justifies the pragmatic approach of increasing p to improve the quality of the procedure.
Nonetheless, we note that the aforementioned difficulty in classically simulating and optimizing large QAOA circuits prevents further scaling conclusions to be drawn.We leave a more complete investigation of these issues, involving the large-scale benchmarking of our scheme on actual quantum computers, as further work.As a final remark, a positive indication is provided in the closely related context of quantum annealing, by the polynomially vanishing energy gap between the ground and first excited state of the Hamiltonian Eq. ( 8) [23], which is strongly correlated to a better performance of QAOA [21].

V. COMPUTATIONS ON QUANTUM COMPUTERS
In this section, we describe an implementation of our scheme for a small-scale network routing problem on the 11-qubit IonQ Harmony quantum computer, accessed through the Amazon Braket cloud quantum computing service with the OpenQAOA Python SDK [19].These experiments serve to empirically verify that the scheme can yield Pareto-optimal solutions to the network routing problem, in a manner consistent with our description in Section.(II).We consider a 4-node fully-connected network and a 6-node square-lattice network, requiring 8 and 11 qubits respectively to encode in QUBO form.For both problem instances, we consider 4 objectives relevant to the network routing problem: data rate, path loss, node delay, and bit-error-rate.Further details on the problem instances can be found in Appendix.(F).
We run the algorithm for the scalarization choice ⃗ w = (1/4, 1/4, 1/4, 1/4), with p = 1 and 2000 shots per circuit execution, and the standard QAOA circuit involving standard single-qubit RX mixers and uniform computational basis state initialisation.Starting at initially suboptimal parameters, COBYLA [35] is chosen as the classical optimization algorithm.
The results of the optimization are shown in Fig. ( 4), which displays the trajectory of the parameters β and γ (black line), plotted on top of the cost function landscape obtained by classical simulation.The inset at each plot shows the evolution of the cost function value during optimization.Starting from suboptimal initial parameters, we observe convergence to local minima for both problems.
To obtain Pareto-optimal solutions of the network routing problem, we sample from the optimized quantum states.Fig. (5) displays the Pareto plot for two pairs of objectives obtained from the IonQ device, with the top k = 125 most probable states marked as coloured crosses (we limit k for visual clarity).Consistent with our discussions in Section.(II), we observe the clustering of high-probability states close to the Pareto front in these plots.The observed success probabilities are also computed to be ≈ 0.008 for the 11-qubit problem and ≈ 0.04 for the 8-qubit problem, which are higher compared to random sampling (3/2 11 ≈ 0.001 for the 11-qubit problem and 4/2 8 ≈ 0.016 for the 8-qubit problem).

VI. CONCLUSION AND OUTLOOK
This work proposes and studies a procedure to obtain Pareto-optimal solutions of MCOPs in a manner amenable to near-term quantum computers via VQAs.By scalarizing the multi-objective cost function in QUBO form and variationally optimizing it with a quantum computer, the procedure allows the Pareto front of MCOPs to be recovered and visualized in an intuitive manner.Focusing on the practically relevant network routing problem, which is an instance of the multi-objective shortestpath problem, efficient QUBO formulations for practically relevant objectives were detailed.Analytical scaling analyses on the resources required for this procedure point towards implementation feasibility for real-world  problem instances, for networks with bounded maximum degrees.This is supplemented by numerical simulations in small scales, which show that high-quality results can be obtained by systematically increasing the depth of the QAOA ansatz.Finally, we tested our approach on an ion trap quantum computer, obtaining consistent results that verify its applicability, constituting one of the first studies of multi-objective optimization on actual quantum computers.
Our suggested framework can be readily generalized to incorporate other MCOPs and VQAs, as long as QUBO formulations of the individual objectives are available.Firstly, to further tailor the scheme we have developed for MCOPs, more general scalarization choices beyond the simple linear scalarization approach could be explored.Examples are the Chebyshev and epsilonconstraint scalarization schemes, which are more powerful and versatile, but incur higher computational costs in the form of additional inequality constraints [1,3].
We can also ask whether it is possible to more directly leverage the fact that QAOA works with quantum superpositions of candidate solutions to the optimization problem.Rather than optimizing based on the energy directly, it may be fruitful to consider modifying the cost function to one that favours states which contain Pareto-optimal configurations with high probability, which shares conceptual similarities with recently developed methods such as the Conditional Value-at-Risk technique [36].Such a modification can be performed efficiently on the sample, since the Hamiltonian corresponding to QUBO problems only contains terms that fully commute with one another.
Another improvement is in the direction of tailored state initializations and mixer designs.Due to the large fraction of infeasible solutions in the state space, the default initial state of the QAOA procedure -an equal superposition of all computational basis states -is highly suboptimal.Ideally one would like to be able to (a) prepare an initial state that is an equal superposition over all feasible solutions, and (b) implement a mixer Hamiltonian that does not induce transitions to infeasible solutions, so that the search is fully restricted to the feasible subspace.In the context of the network routing  problem we leave this as a question for future work, and note that [37] describes QAOA mixer Hamiltonians for network flow optimization.We also remark that the two aforementioned requirements (a) and (b) may carry a significant resource footprint, depending on the complexity of preparing the initial state and the polynomial degree of the mixer Hamiltonian.Even on classical computers, MCOPs are generally NP-hard optimization problems that are solved with metaheuristic optimization algorithms such as particle swarm and genetic algorithms (c.f.Appendix.B).Can our scheme, or VQAs in general, outperform them in some respect at large scales that are challenging for classical approaches?Due to the difficulty to theoretically analyze the performance of QAOA at large scales, and the intractability of classical simulations of quantum computations, general claims about performance likely must be made empirically, on real quantum computers.While we were able to implement and obtain solutions for the procedure for small-scale problems, applications in practical settings will require larger and more reliable quantum computers than those available today.Our work serves as a first step towards answering this question, and motivates the need for large-scale experiments on quantum computers.
Appendix A: Survey of relevant work This section briefly discusses relevant work in the usage of near-term quantum algorithms in solving MCOPs and the shortest path problem.
Applications of quantum annealing to the shortest path problem were first studied in [22], which explored different QUBO formulations of the shortest path problem.The authors of that work provide numerical and empirical evidence on its performance on a quantum annealer.With the same formalism, [23] studied a related problem in the context of chemistry, by again solving a relatively large problem on a quantum annealer, providing numerical evidence on its efficiency with increasing system size through a polynomially (rather than exponentially) vanishing ground state energy gap.
On the other hand, the application of near-term quantum algorithms to MCOPs is a relatively unexplored area.[38] provides theoretical arguments verifying that adiabatic quantum algorithms can be applied to find Pareto-optimal points corresponding to the solution of scalarized objective functions.Most recently, [39] applied QAOA to the multi-objective routing problem in the context of 6G communication networks by sequentially solving single-objective problems with a lexicographic ordering method.
In relation to the above references, our work explores the solution of the multi-objective routing problem with QAOA through scalarization, and provides in-depth analysis on the approach's scalability and performance.The approach can be generalised to VQAs beyond QAOA, MCOPs beyond the routing problem, and techniques beyond linear scalarization.To the best of our knowledge, experiments on actual quantum hardware to solve multi-objective optimization problems (described in Section.V) are new.

Appendix B: Multi-objective Optimization
This section provides a brief review of relevant notions in MOOPs and their solutions.
Given m objective functions f i : S → R, i = 1, ..., m that map objects from a state space S to real numbers, multi-objective optimization asks for states that are optimal with respect to all of the objectives {f i }.Depending on their forms, the objectives generally produce competing effects with one another, so states that are simultaneously optimal in all objectives generally do not exist.In this case, we ask instead for a set of Paretooptimal/efficient states, which are optimal in the sense that no other states which improve on at least one individual objective without deteriorating in others can be found.The set of all Pareto-optimal states is also conventionally referred to as the Pareto front, and they are always located at the boundary of the region occupied by the set of states [2].When the state space S is finite -as is the case for problems we consider -this is called a multi-objective combinatorial optimization problem (MCOP).Depending on the form of the objectives, MCOP's are generally difficult, belonging to the class of NP-hard problems [11,12].
Classically, the most common approach to solve multiobjective optimization problems is with scalarization techniques.Instead of optimizing the set of objectives {f i } directly, scalarization techniques aggregate them into a single objective function, thereby converting a multi-objective problem into a single-objective one.This amounts to projecting an m dimensional vector onto a line, and finding the optimal state within it.Depending on the forms of {f i } and the resources required to compute the objectives, different scalarization techniques are employed -this is a well-studied subject in the literature [1,3].For our purposes, we consider linear scalarization, where the final aggregated objective is a convex sum of the m objectives: where the weights w i are real numbers.This corresponds to the projection of the m dimensional vector onto a straight line with gradient ⃗ w.A property of this procedure is that a solution of the linearised problem is also a point on the Pareto front, guaranteeing the Paretooptimality of the linearised problem's solution.Therefore, by optimizing over different choices of the weights, we can in principle build up a set of different points in the Pareto front.
While this is one of the simplest scalarization techniques, it preserves the quadratic form of the objectives and constraints in our MCOP, and also allows the scalarization weights w i to be interpreted as a priori preferences that the decision maker can select.In the case where the objectives are convex, it can be shown that every state in the Pareto front corresponds to a solution of a linearised problem.However, this simple linear scalarization may not be sufficient to capture the entire Pareto front if the objectives are concave, where some points in the Pareto front may never be captured as solutions of any linearised problem.This can be overcome with more advanced scalarization techniques such as Chebychev scalarization [1,3], at the expense of introducing inequality constraints.In the context of our work, if a non-linear scalarization scheme is used, the resulting cost Hamiltonian may contain higher-order terms, resulting in cost unitary circuit with greater depth due to the need to implement higher-order gates.Alternatively, one can consider a more general variational ansatz, and optimize over a non-linear cost Hamiltonian.
Finally, to optimize the scalarized objective classically, suitable heuristic optimization algorithms are usually employed.Algorithms that work very well in practice are population (e.g.ant colony methods) and evolutionary or genetic-based algorithms, such as NSGA-II [1,3].
More broadly, approaches to solve MOOPs can be broadly categorized into three classes: • A priori methods take the preferences (e.g. the weights/relative importance of each objective) of the engineer/decision maker into account prior to the optimization, and adapts the optimization process based on this preference.Scalarization falls into this class of methods, where the weights of the objectives are specified in an a priori manner.
• A posteriori methods aim to solve for either a representative subset or all possible Pareto-optimal solutions, only taking preferences of the decision maker into account after the optimization process.
• Interactive methods are adaptive and iterative methods that require the continuous interaction of the decision maker at each step of the optimization process.
As discussed in the main text, our scheme belongs to the class of a priori or a posteriori methods, depending on whether the choice of scalarization weights is made explicitly.
Appendix C: QUBO Formulation of the multi-objective shortest path problem Following the problem statement in Section III, this appendix details the QUBO formulation of the multiobjective shortest path problem.We begin by providing an encoding of the state space of the problem in terms of binary variables, before detailing the form of the quadratic cost function to be optimized.More detailed analyses in the context of annealing, including an extension to directed graphs, can be found in [22].
An encoding of the problem can be achieved with |V | + |E| − 2 binary variables representing each state x in the state space [22] (We will use x i 's and x ij 's interchangeably with x to denote states when the context is clear from now on.).The first |V | − 2 variables x 1 , x 2 , ..., x |V |−2 correspond to the nodes in the graph excluding the source and destination nodes, while the remaining |E| variables x ij , where (i, j) ∈ E, correspond to the edges.With this encoding, we consider objectives that take the following forms: • Node cost: Associates each node in the graph with a cost, depending on the node's weight V i .
• Edge cost: Associates each edge in the graph with a cost, depending on the edge's weight E ij .
With only one objective (for instance, the minimization of the distance between two points of a graph, which is an edge cost), the problem reduces to a single-objective shortest path problem which can be solved efficiently in polynomial time by Dijkstra's algorithm [40], which takes O(|E| + |V | log |V |) steps.The presence of multiple objectives constitute a MCOP.Eqs.(C1) and (C2) allow us to compute the cost vector associated with a state x.However, of the 2 |V |+|E|−2 possible states, not all of them represent actual paths.A valid path is one that (1) starts from a specified source node s, (2) ends at a specified destination node d, and (3) has no broken links or branches along the path from s to d. States that satisfy these criteria are called feasible, and infeasible otherwise.These 3 constraints can be enforced as quadratic penalty terms added to our cost function, so that infeasible solutions have higher total energies than feasible ones: • Source constraint: Penalizes paths that do not have exactly one edge connected to the source node: where s is the index of the source node, and the sum is over all edges that are connected to node s.This term has a minimum value of −1, which occurs for states where the source node is used (x s = 1), and there is only one way of leaving the source node (the bracketed term is equal to zero).
• Destination constraint: Penalizes paths that do not have exactly one edge connected to the destination node: where d is the index of the destination node, and the sum is over all edges that are connected to node d.This term has a minimum value of −1, similar to the source constraint.
• Path constraint: Penalizes paths that do not have exactly 2 edges connected to intermediate nodes: E path = i∈V E i , such that for each intermediate node i: where the sum is over all edges that are connected to node i.This has a minimum value of 0, which occurs for states where for all intermediate nodes used (for which x i = 1), the node degree is equal to 2.
A solution of the multi-objective shortest path problem is then a path x that achieves the minimum of the penalties Eqs.(C3),(C4), and (C5) and is Pareto-optimal with respect to the node and edge costs.
As mentioned in the main text, the connectivity of the graph associated with the sum of the path constraints is precisely that of the problem's middle graph, which we review and exploit in Appendices D 2 and D 3 respectively.Furthermore, the constraints Eqs.(C3) and (C4) only incur a small number of linear and quadratic terms (dependent on the local connectivity of the source and destination nodes) which we safely ignore for conciseness.
Appendix D: Details on Graph Theory, the Middle Graph, and the Proof of Eq. ( 11) Here, we list down notions in graph theory used in the resource estimation part of the main text, and provide a proof for Eq.(11).
As mentioned in Section.III, the correspondence G H = M (G) implies that the solution of a shortest path problem defined on a network G can be mapped to the ground state of a Hamiltonian/Ising model H with graph G H = M (G).Ignoring the quadratic terms resulting from the source and destination nodes, which only contributes a negligible constant overhead, Eq. ( 11) implies that the implementation of the cost unitary is efficient, depending only linearly on the maximum degree of the network G.

Notions in graph theory
For an undirected graph G = (V (G), E(G)), the degree of a vertex u ∈ V (G), denoted deg(u), is the number of edges that are incident to the vertex.The maximum degree, denoted ∆ G , is the maximum of its vertices' degrees ∆ G = max u∈V (G) deg(u).The degree sum formula relates them to the number of edges of G: The line graph of G = (V (G), E(G)) is the graph L(G) = (E(G), E ′ ) with the edges of G as vertices, such that they are adjacent if and only if their corresponding edges are incident on the same vertex.The maximum degrees of G and L(G) are related by: ∆ L(G) ≤ 2∆ G − 2. (D2) The middle graph of G = (V (G), E(G)) [27] is the graph M (G) = (V (G)∪E(G), E ′ ), with vertices u, v that are adjacent if either: (a) u is a vertex in G and v is an edge in G incident to u, or (b) u and v are edges in G that both incident on the same vertex.
The endline graph of G, denoted G + , is defined as the graph obtained from G by adding to each of its nodes an end-vertex (i.e. an edge that is connected to a single node).The attachment of an end-vertex to each node immediately implies that: An edge coloring of a graph G is an assignment of colors to its edges, such that no two edges sharing a vertex can have the same color.The edge chromatic number/chromatic index of G, denoted χ ′ (G), is the minimum number of such colors needed.Vizing's theorem relates the edge chromatic number of G with the maximum degree for any graph: (D4)

Properties of the middle graph
We introduce some properties on M (G) that will be used: This immediately implies that their maximum degree are also equal: The following chain of (in)equalities are true: where applying Vizing's theorem yields the first inequality, M (G) ∼ = L(G) + yields the second equality, Eq. (D3) yields the third inequality, and Eq.(D2) yields the final equality.
this metric in QUBO form as a linear sum of edge costs, we use a negative exponential function to weight the data rates Γ ij : x ij • e −βΓij , (E6) so that sub-optimal links with low data rates are penalized heavily.Here, β is a graph-dependent coefficient that should be chosen such that for any choice of transmission path x, the dominant contribution to E Γ corresponds to the minimum data rate along the path.Figure 6: Scaling of the approximation ratio and success probability with the number of layers for different initializations, for one specific problem.To produce the random initializations curves, we perform 100p random initializations and optimize with QAOA.The lighter region displays standard deviation around the mean.The average result from 100p random initializations is labelled as 'random mean', the best result from the 100p random initializations is labelled as 'random max' (to refer to the maximum approximation ratio achieved), and the (one) result from the linear ramp initialization is labelled as 'linear ramp'.Table.VI provides the complete specification of the triangular lattice graph and the multi-objective cost function for the problem instance used in Fig. (2).That is, each of the 4 objective functions takes the form: where h i corresponds to node weights and J ij corresponds to edge weights with values as assigned in the table.The problem can then be specified via Eq.(H1).
tion, we can minimally 'perturb' it by flipping the value of individual variables to obtain infeasible configurations.However, we could build 'even worse' infeasible configurations by introducing such perturbations in multiple locations simultaneously, not only around a single vertex or edge.The number of combinations of infeasible config-urations we can generate in this way grows exponentially in the number of variables (i.e. in the number of vertices and edges).We therefore conclude that the fraction of feasible to infeasible configurations shrinks exponentially in the problem size.
(1) -M different scalarization weights are chosen corresponding to M QAOA runs, resulting in M sets of candidate solutions B 1 , .., B M which are aggregated as M i=1 B i .A classical method then checks for Pareto-optimality within this set, yielding a set of feasible solutions that approximates the Pareto front.

Figure 1 :
Figure 1: Summary of workflow to approximate the Pareto front of a MCOP.(a) Visualization of the conversion from the problem's underlying graph to the Hamiltonian corresponding to the QUBO formulation of the MCOP.The graph (with solid lines denoting its edges and circles denoting its nodes) is converted to the Hamiltonian Eq. (8) that incorporates the objective and constraints of the problem.The Hamiltonian itself can be represented by a graph, G H (bottom right graph, with dotted lines with dotted lines denoting edges and circles denoting qubits).The illustration describes the conversion for the multi-objective routing problem, which is elaborated in Section.III.(b) Quantum-classical optimization with QAOA to solve the optimization problem Eq. (2) for the scalarized cost Hamiltonian Eq. (8).The p-layer QAOA circuit is run at different angles according to the feedback of a classical optimizer, navigating the cost landscape until it converges to a minima.(c) Aggregation over results of one or more scalarization choices, determination of Pareto front with classical solver, and extraction of final result(s).

Fig. 2 .•
(a)  shows the Pareto plot for our 13-qubit network routing problem (with parameters specified in Table.H).The problem involves 4 objective functions, and the figures show a projection onto the plane of two of those objectives.Each point in the plot corresponds to a possible state (a bitstring), with feasible states marked as blue and Pareto-optimal solutions marked as red.The top-k most probable states from |ψ( ⃗ β * , ⃗ γ * )⟩ are marked with crosses and color coded according to their probability (colored crosses).In Fig. 2.(b),(c),(d) we illustrate the solution returned by QAOA with different scalarization weights ( ⃗ w = (0, 1, 0, 0) for (a), ⃗ w = (1/4, 1/4, 1/4, 1/4) for (b), and ⃗ w = (1, 0, 0, 0) for (c)).A few important observations from these plots are the following: The states with highest sampled probability in the output state |ψ( ⃗ β * , ⃗ γ * )⟩ are in the lowenergy sector of the Pareto plot : This is a direct consequence of the quantumclassical optimization loop that results in an output quantum state |ψ( ⃗ β * , ⃗ γ * )⟩ that has minimal scalarized energy when optimized.Sampling from |ψ( ⃗ β * , ⃗ γ * )⟩ then leads to a set of candidate states B i that are biased to contain high-quality/lowenergy states of the scalarized problem, which we can observe from the clustering of the colored crosses in the bottom left region of the Pareto plot.• States sampled from |ψ( ⃗ β * , ⃗ γ * )⟩ are squeezed in the direction of the scalarization line : Scalarization amounts to projecting points on the Pareto plot to the line passing through the origin with gradient equal to the scalarization weight vector (black dotted line).Different choices of scalarization therefore result in different regions in the

Figure 2 :
Figure 2: Pareto plot of a 4-objective problem, projected on two of the metrics.Each point in the plot corresponds to a possible state/solution (for a total of 2 n qubit points), with feasible solutions marked as blue and Pareto-optimal solutions marked as red.We also select the top k = 100 most probable states from |ψ( ⃗ β * , ⃗ γ * )⟩ and color code them according to their probability (colored crosses).We plot the results for 3 different scalarization choices, with the gradient of the scalarization weight vector visualized (black dotted lines passing through origin).The scalarization vectors are ⃗ w = (0, 1, 0, 0) for (a), ⃗ w = (1/4, 1/4, 1/4, 1/4) for (b), and ⃗ w = (1, 0, 0, 0) for (c).The underlying network is a triangular lattice graph with 6 nodes and 9 edges, requiring 13 qubits to encode.

Figure 3 :
Figure 3: Scaling of the approximation ratio and success probability with the number of layers for triangular lattices, square lattices, and cycle graphs (left to right columns).Different colors indicate different problem sizes, with standard deviations displayed in shaded regions.For cycle graphs (rightmost column), the legend indicates the number of nodes besides the soure and destination nodes.The data was produced after optimization with linear ramp initialization, two metrics (node based and edge based respectively), with node/edge weights uniformly randomized between -1 and 1, averaged over 50 instances.scalarization weights are chosen to be equal for both objectives, i.e ⃗ w = (0.5, 0.5).

Figure 4 :
Figure 4: Trajectory of the QAOA parameters during optimization (black line) for both the 8-qubit fully-connected network (left plot) and 11-qubit square lattice networks (right plot), overlaid on the full cost function landscape obtained by classical simulation of the QAOA circuits.Initial parameters are displayed as black crosses.Inset shows cost vs iteration plots of the optimization.

Figure 5 :
Figure 5: Pareto plots obtained from sampling the optimized quantum state of QAOA with the IonQ device for both the 8-qubit fully-connected network (left two plots) and 11-qubit square lattice networks (right two plots).Each point in the plot corresponds to a possible state, with feasible solutions marked as blue and Pareto-optimal solutions marked as red.The top k = 125 most probable states sampled are plotted as crosses and color coded according to its probability.The gradient of the scalarisation weight vector is also visualised (black dotted lines).

Appendix H:
Specification of problem instance of Fig. (2) Table VI: Weights of cost terms of the triangular lattice problem instance of Fig. (2) for all 4 objectives, labelled according to their corresponding edges or nodes.6 node weights and 9 edge weights are present.

Table V :
It can be estimated by either considering the link in G with the smallest possible data rate, or by determining the value of c(β) = x∈Paths e −βR x min − j̸ =min e −βR x j by sampling over different paths.R x min is the minimum data rate of path x, and R x j are the data rates on all other edges in x. edge distance data rate path loss Rij bit error (i, j) dij [m] Γij [kbit/s] Lij [dB] pij [dB] Edge properties and parameters of the square lattice network.