# Efficient And Scalable Branch-and-Bound Algorithm for Exact Qubit Allocation

Jean-Philippe Valois, Guillaume Helbecque, Nouredine Melab

► To cite this version:

Jean-Philippe Valois, Guillaume Helbecque, Nouredine Melab. Efficient And Scalable Branch-and-Bound Algorithm for Exact Qubit Allocation. Future Generation Computer Systems, 2026, 179, pp.108342. 10.1016/j.future.2025.108342. hal-05449040

HAL Id: hal-05449040

https://hal.science/hal-05449040v1

Submitted on 14 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Efficient and scalable branch-and-bound algorithm for exact qubit allocation

Jean-Philippe Valois ID*, Guillaume Helbecque ID, Nouredine Melab ID

Université de Lille, CNRS/CRISTAL UMR 9189, Centre Inria de l'Université de Lille, France

## ARTICLE INFO

**Keywords:**

Quantum computing

Qubit allocation

Noisy intermediate-scale quantum systems

Quadratic assignment problem

Parallel branch-and-bound

Ultra-scale computing

## ABSTRACT

Qubit allocation is a central step in adapting abstract quantum circuits to noisy intermediate-scale quantum devices, yet exact approaches for solving it face severe scalability limitations. In this work, we revisit the formulation of qubit allocation as a permutation-based quadratic assignment problem and develop a branch-and-bound algorithm for its exact resolution. We first establish a refined sequential implementation that achieves significantly faster runtimes than previous exact approaches on most problem instances, thereby setting a new state-of-the-art for this formulation. Building on this foundation, we extend the approach to a performance-aware parallel implementation that exploits both intra-node and inter-node parallelism on High-Performance Computing (HPC) infrastructures. Our experimental evaluation demonstrates near-linear strong scaling at the intra-node level and substantial scalability in distributed settings across nodes. Leveraging these capabilities, we provide reference optimal solutions for challenging benchmark circuits of up to 26 qubits—significantly larger than previously reported instances. These results show that large-scale parallelization can effectively extend the reach of exact methods for qubit allocation, thereby advancing the integration of combinatorial optimization and HPC techniques in quantum computing.

## 1. Introduction

Quantum computers have the potential to significantly accelerate the resolution of problems that are intractable for classical systems, by leveraging principles of quantum mechanics such as superposition and entanglement. Prominent examples include quantum search algorithms [1], quantum simulation of physical systems [2], quantum machine learning [3], and efficient integer factorization used in cryptography [4]. However, quantum computers remain in the realm of Noisy Intermediate-Scale Quantum (NISQ) devices. These systems typically feature tens to a few hundred qubits (quantum bits) and are characterized by limited reliability as prone to non-negligible error rates.

Recently, IBM has made quantum processing units with capacities ranging from 127 to 156 qubits publicly accessible through the IBM Q platform [5]. While these devices enable users to execute quantum programs, the execution process, much like in classical computing, requires an initial compilation stage. To be precise, in the quantum context this process is a transpilation, and involves two major phases. First, the target algorithm is transformed into a logical circuit composed of universal quantum gates. Then, this logical circuit is further translated into a physical circuit that conforms to the hardware constraints. This latter phase, known as *qubit mapping*, plays a critical role in the overall error rate of

the resulting quantum circuit [6]. Qubit mapping itself consists of two consecutive sub-tasks: *qubit allocation*, which assigns logical qubits to physical ones, and *dynamic routing*, which inserts additional operations to satisfy connectivity constraints during execution. In this work, our main focus is given to the qubit allocation problem.

Existing exact approaches to qubit allocation all face a fundamental limitation: poor scalability. In this work, we revisit an approach that formulates qubit allocation as a permutation problem, enabling its exact resolution through a more efficient Branch-and-Bound (B&B) algorithm. Our main objective is to overcome the scalability barrier by enhancing this method with large-scale parallelism. The resulting parallel B&B implementation enables exact reference solutions for problem instances involving quantum circuits of substantially larger sizes than previously studied, all while keeping computational requirements within practical limits. Overall, this study lies at the intersection of quantum computing, combinatorial optimization, and High-Performance Computing (HPC).

The primary contributions of this work are summarized as follows.

* We revisit an existing approach that formulates qubit allocation as a combinatorial optimization problem, and develop a refined sequential B&B algorithm for its exact resolution. This implementation achieves substantially faster runtimes than previous exact

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: jean-philippe.valois@inria.fr (J.-P. Valois), guillaume.helbecque@univ-lille.fr (G. Helbecque), nouredine.melab@univ-lille.fr (N. Melab).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of standard quantum gates: Pauli-X, Hadamard and CNOT.</figcaption></figure>

approaches on most problem instances, establishing new state-of-the-art performance for this formulation, while also serving as a baseline for evaluating the parallel variant.

* We extend the previous design to a performance-aware parallel B&B algorithm, leveraging the Chapel-based P3D-DFS framework [7], with an implementation specifically tailored to the considered problem formulation. The performance and scalability of this design are evaluated on challenging qubit allocation instances, demonstrating its ability to effectively exploit the computational power of modern HPC clusters.

* We further showcase the capabilities of our parallel B&B implementation through large-scale experiments, providing exact reference solutions for the qubit allocation of quantum circuits with up to 26 qubits.

The remainder of this paper is organized as follows. Section 2 introduces the qubit allocation problem, including a background on quantum computing and a review of related works. Section 3 presents the design of our B&B algorithm applied to qubit allocation, and describes its parallelization at scale. These algorithms are then evaluated and compared through extensive experiments in Section 4. Finally, Section 5 concludes the study and outlines promising directions for future work.

## 2. The qubit allocation problem

### 2.1. Background on quantum computing

This section briefly introduces the main concepts of quantum computing for a better understanding of the qubit allocation problem.

**Qubit.** In contrast with classical bits, quantum computers use *qubits* (quantum bits), which also have two basis states, denoted $|0\rangle$ and $|1\rangle$, but with a crucial difference: a qubit can exist in a linear combination of these basis states. Mathematically, the state of a qubit is written as $|\Psi\rangle = \alpha|0\rangle + \beta|1\rangle$, where $\alpha, \beta \in \mathbb{C}$. The complex amplitudes $\alpha$ and $\beta$ define the probability of measuring the qubit in a particular basis state, with $|\alpha|^2$ corresponding to $|0\rangle$ and $|\beta|^2$ to $|1\rangle$. In other words, unlike classical bits, a qubit can be in both the 'zero' and the 'one' states simultaneously. This phenomenon is known as *superposition*.

**Quantum Gates.** Quantum computers use *quantum logic gates* to perform unitary operations on qubits. These gates operate on either single qubits or multiple qubits. Commonly used single-qubit gates include the Pauli-X gate (also known as the quantum NOT gate) and the Hadamard gate (*H*) (cf. Fig. 1). Mathematically, single-qubit gates are represented as $2 \times 2$ unitary matrices that act on the qubit's state vector $(\alpha, \beta)$. For instance, applying an *X* gate to a qubit in the state $|\Psi\rangle = \alpha|0\rangle + \beta|1\rangle$ yields the new state $|\Psi'\rangle = \beta|0\rangle + \alpha|1\rangle$. Particular attention will be given to two-qubit gates, the most prominent being the Controlled-NOT (CNOT). Acting on a control and a target qubit, the CNOT applies an *X* gate to the target if the control is in state $|1\rangle$, and leaves it unchanged otherwise (cf. Fig. 1). Together with all single-qubit gates, the CNOT forms a universal gate set [8], meaning any quantum computation can be decomposed into such operations. For this reason, the CNOT is regarded as the canonical multi-qubit gate in both theory and practice.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of a quantum circuit composed of Hadamard and CNOT gates, followed by measurements.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Connectivity graph of the IBM Boeblingen [5] quantum device featuring $N = 20$ physical qubits. Nodes represent qubits and edges indicate available two-qubit interactions.</figcaption></figure>

**Quantum Circuit.** A quantum circuit is typically represented as a diagram in which each horizontal line is associated to a qubit and defines the time sequence in which quantum gates are executed from left to right. Fig. 2 illustrates a simple 3-qubit circuit composed of two Hadamard gates and two CNOT gates. At the right end of each qubit line, measurement operations are shown. These represent a special class of operations, usually performed at the end of a computation, to extract classical information from the quantum state. In essence, quantum circuits manipulate the probability amplitudes of multi-qubit systems through sequences of quantum gates, steering the system toward a configuration from which the desired result can be obtained upon measurement.

**NISQ Hardware and Circuit Mapping.** This work considers *Noisy Intermediate Scale Quantum* hardware, which characterizes the current generation of quantum devices. These systems are inherently error-prone: gate operations and measurements are subject to noise, introducing a significant probability of incorrect outcomes. They also remain limited in scale, usually featuring on the order of tens to hundreds of qubits, with only sparse connectivity. The physical architecture of an NISQ device is commonly described by a connectivity graph, where vertices represent physical qubits and edges denote possible two-qubit interactions (cf. Fig. 3). Connectivity is particularly critical for implementing two-qubit gates such as the CNOT. While theoretical circuit models assume arbitrary qubit interactions (cf. Fig. 2), real devices only permit them along physical connections. Consequently, most circuits cannot be executed directly and must be transformed, typically through the insertion of SWAP gates to enable required interactions. Addressing this challenge defines the qubit mapping problem: adapting an abstract quantum circuit to the physical and topological constraints of the device.

### 2.2. Qubit allocation

The qubit mapping problem consists of two consecutive sub-tasks: qubit allocation and dynamic circuit routing. Qubit allocation involves assigning each logical qubit in the input quantum circuit to a physical qubit on the target quantum device. Based on this initial assignment, the subsequent task, dynamic routing, modifies the mapping by inserting SWAP gates to ensure that each two-qubit gate in the circuit is executable given the hardware's connectivity constraints. A SWAP gate
<figure><img src="image_6.png" /><figcaption>Fig. 4. Decomposition of a SWAP gate into an equivalent sequence of three CNOT gates.</figcaption></figure>

exchanges the states of two adjacent qubits and is typically decomposed into a sequence of three CNOT gates as shown in Fig. 4.

In this work, we focus primarily on the qubit allocation task, which can be framed as an optimization problem. The objective is to determine an initial mapping of logical to physical qubits that minimizes the number of SWAP gates required during the routing phase, thereby improving the fidelity and efficiency of the resulting circuit. This is particularly important because SWAP operations introduce additional error sources and significantly increase the overall circuit depth, which corresponds to the length of the longest path of sequential gate operations in the circuit. A deeper circuit is more susceptible to qubit decoherence, the process by which quantum states lose their coherence (or reliability) over time due to environmental interactions. As a result, minimizing the number of SWAP gates is essential for maintaining computational accuracy and is thus a central objective in quantum compilation.

The example in Fig. 5 illustrates the principle of initial qubit mapping, also known as qubit allocation, and highlights its impact on the reliability and efficiency of the resulting quantum circuit. The goal is to assign the logical qubits {$q_0, q_1, q_2$} of the 3-qubit circuit shown in Fig. 5a to the physical qubits {$P_0, P_1, P_2$} of a linear quantum device, whose coupling graph is depicted in Fig. 5b. The assignment must be bijective, mapping each logical qubit to a unique physical qubit. For instance, the mapping {$q_0 \to P_2, q_1 \to P_0, q_2 \to P_1$} is one valid configuration, which can be compactly represented as the permutation {$2, 0, 1$}. We compare hardware-compliant versions of the original circuit obtained from two different initial allocations. In Fig. 5c, the allocation {$0, 1, 2$} results in a situation where the CNOT gate between $q_0$ and $q_2$ cannot be directly executed due to the device's connectivity constraints. As a result, a SWAP gate must be inserted, increasing the total circuit depth to six, given that a SWAP consists of three consecutive CNOT gates (cf. Fig. 4). In contrast, the allocation {$1, 0, 2$}, shown in Fig. 5d, avoids this issue, allowing all two-qubit gates to be executed natively without any additional SWAP operations. Consequently, the final circuit retains a depth of only three. This simple example clearly demonstrates how an effective qubit allocation can significantly reduce both the gate count and depth of an executable circuit, thereby improving the overall reliability of quantum computations.

## 2.3. Related works

The qubit mapping problem has been the focus of numerous studies, many of which rely on heuristic-based approximate methods to cope with its computational complexity. Li et al. [9] introduced the SABRE heuristic search method, designed for NISQ devices with arbitrary qubit connectivity. Their successful approach was later integrated into IBM's Qiskit framework [10], widely used for designing and executing quantum programs on both simulated and real hardware. More recently, Rouzé et al. [11] proposed a parallel genetic algorithm for qubit mapping, which produced higher-quality solutions than SABRE on a variety of instances, while exploiting parallelism to accelerate the transilation phase.

In contrast, exact methods for qubit mapping have also been the subject of numerous studies. While approximate approaches remain indispensable for scaling to large problem sizes, exact methods serve a complementary purpose. Beyond providing reference solutions to assess heuristic quality, exact methods also offer insights into the structural properties of the mapping problem, establish rigorous performance bounds, and inspire the design of more effective heuristics. Among these

exact approaches, Zulehner et al. [12] employed an A\* search strategy to develop an efficient solver for the qubit allocation problem, though their method was restricted to linear nearest-neighbor architectures. Subsequently, Zhu et al. [13] reformulated qubit allocation as a Quadratic Assignment Problem (QAP) and solved it using a B&B algorithm. Building on this direction, Zhang et al. [14] extended the scope by addressing the full qubit mapping pipeline within a more advanced guided search framework.

A common limitation across existing exact methods is their lack of scalability, primarily due to the exponential growth of the search space and the inherent complexity of the problem. In particular, studies [13] and [14] were restricted to circuits involving no more than 16 qubits. Motivated by this challenge, our work revisits the qubit allocation problem through a B&B approach, with the aim of substantially improving scalability via ultra-scale parallelism. To the best of our knowledge, no large-scale parallel optimization algorithms currently exist for this problem. Bridging this gap would make it possible to compute exact reference solutions for quantum circuits beyond the 16-qubit threshold. To this end, we adopt the problem formulation of Zhu et al. [13], framing qubit allocation as a combinatorial optimization task, specifically a permutation-based instance of the QAP. We likewise address the problem using a B&B algorithm, but with a distinct algorithmic design.

## 2.4. Mathematical problem formulation

Before formulating the qubit allocation problem, we begin by introducing key definitions that capture the essential properties of the quantum circuit and the quantum device on which it is to be executed.

**Definition 1 (Interaction frequency matrix).** Given a quantum circuit involving $n$ logical qubits, the interaction frequency matrix is the $n \times n$ symmetric matrix $\mathbb{F} = (f_{ij})$, where each entry $f_{ij}$ denotes the total number of two-qubit gates (i.e., interactions) between logical qubits $q_i$ and $q_j$ throughout the circuit.

**Definition 2 (Coupling distance matrix).** Given a NISQ device with $N$ physical qubits, the coupling distance matrix is the $N \times N$ matrix $\mathbb{D} = (d_{ij})$, where $d_{ij}$ represents the minimum number of SWAP operations required to enable a two-qubit interaction between physical qubits $P_i$ and $P_j$. As the coupling graph of an NISQ device is typically undirected, $\mathbb{D}$ is symmetric.

The matrices $\mathbb{F}$ and $\mathbb{D}$ encapsulate all the necessary information derived from the quantum circuit and the target hardware, respectively. These matrices serve as inputs to the qubit allocation problem, enabling a formulation as a QAP. In practice, to ensure dimensional compatibility with the $N \times N$ device matrix, we pad $\mathbb{F}$ with zeros whenever $N > n$, thereby introducing $N - n$ dummy logical qubits that do not participate in any two-qubit interactions. Fig. 6 provides an illustrative example of a 3-qubit circuit and a 4-qubit device, along with their corresponding interaction frequency and coupling distance matrices. The coupling distance matrix $\mathbb{D}$ can be computed using a standard Breadth-First Search (BFS) traversal of the unweighted and undirected coupling graph. For any $(i, j) \in [0, N-1]^2$, the value $d_{ij}$ equals the shortest path length between $P_i$ and $P_j$ (measured in number of edges) minus one, reflecting the number of SWAP gates required to connect the two qubits.

With these definitions established, we are now equipped to formally express the qubit allocation problem as a QAP.

**Definition 3 (Qubit allocation as a QAP).** Given a quantum circuit involving $n$ logical qubits, characterized by an interaction frequency matrix $\mathbb{F} \in \mathcal{M}_n(\mathbb{N})$, and a quantum device with $N \ge n$ physical qubits, described by a coupling distance matrix $\mathbb{D} \in \mathcal{M}_N(\mathbb{N})$, the objective is to find a partial permutation $\pi \in \mathcal{S}_n^N$ that minimizes the route cost function $RC(\pi)$:

$$ \min_{\pi \in \mathcal{S}_n^N} RC(\pi). \qquad (1) $$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Illustration of the qubit allocation problem and its impact on circuit transformation: (a) An example quantum circuit. (b) Linear coupling graph of the target device. (c) Allocation {0, 1, 2} introduces an additional SWAP, resulting in a depth of 6. (d) Allocation {1, 0, 2} yields a more efficient circuit with depth 3.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Inputs of the qubit allocation problem and their mathematical representation: (a) quantum circuit with two-qubit gates, (b) linear coupling graph of the target device, (c) interaction matrix $\mathbb{F}$ capturing gate frequencies, and (d) distance matrix $\mathbb{D}$ encoding device connectivity constraints.</figcaption></figure>

Where $S_n^N$ denotes the set of all injective mappings (partial permutations) from $n$ logical qubits to $N$ physical qubits, and the cost function $RC$ is defined as:

$$RC(\pi) = 2 \sum_{i=0}^{n-1} \sum_{j=i+1}^{n-1} f_{ij} d_{\pi(i)\pi(j)}. \qquad (2)$$

The QAP (1) is not only NP-hard to solve [15], but is also NP-hard to approximate within any constant factor [16], making it among the most challenging combinatorial optimization problems. The search space $S_n^N$ has cardinality $N!/(N-n)!$, corresponding to the number of injective assignments of $n$ logical qubits into $N$ available physical qubits.

Because the matrices $\mathbb{F}$ and $\mathbb{D}$ are symmetric with zero diagonal entries, the double summation in Eq. (2) is taken only over the upper triangular part of the matrices (i.e., $i < j$), thereby avoiding redundancy. Each term $f_{ij}d_{\pi(i)\pi(j)}$ represents the cost of executing the interactions between logical qubits $q_i$ and $q_j$, given that they are mapped to physical qubits $P_{\pi(i)}$ and $P_{\pi(j)}$, respectively. The distance $d_{\pi(i)\pi(j)}$ reflects the number of SWAP operations required to enable the corresponding two-qubit gate.

For any candidate solution $\pi \in S_n^N$, the route cost $RC(\pi)$ can be interpreted as an estimate of the total number of SWAP operations required to make the quantum circuit executable on the target hardware, starting from the initial qubit mapping defined by $\pi$. This estimate is based on a naive dynamic routing strategy, which operates under a simple rule: whenever a two-qubit gate involves non-adjacent physical qubits, SWAP gates are inserted to bring the qubits next to each other, the gate is executed, and then the same number of SWAP gates restore the original positions of the qubits. This round-trip movement justifies the factor of two in the definition of the route cost function in Eq. (2).

It is important to recognize that, for a given circuit and initial mapping, multiple routing strategies may exist, some more efficient than others. This exposes a key limitation of the QAP formulation, which relies on a naive routing strategy to estimate mapping cost. In practice, more sophisticated routing methods can reduce the number of required SWAPs, yielding more efficient compiled circuits than those predicted by the QAP model alone. Combining the QAP-based allocation approach with techniques that handle the dynamic routing phase—such as heuristic routing algorithms—could extend this work toward a complete qubit mapping pipeline, although such an investigation lies beyond the scope of this study.

## 3. Branch-and-bound algorithm design

Branch-and-bound is an exact combinatorial optimization method that systematically explores the space of feasible solutions by recursively dividing the original problem into smaller subproblems. During exploration, it discards (or prunes) branches of the search tree that cannot contain the optimal solution, based on bounding evaluations. More precisely, the algorithm relies on four key operators<fnref n="1" />:

* **Branching.** Decompose the current problem into smaller subproblems, forming nodes of a search tree, commonly referred to as the state-space tree.

* **Selection.** Determine the next node (or subproblem) to explore based on a predefined strategy (e.g., depth-first, best-first).

* **Bounding.** Compute a lower bound for the objective function in the current subproblem. This bound represents the best possible value that any complete solution descending from this node could achieve.

* **Pruning.** Eliminate subproblems whose lower bounds are no better than the best-known solution, as they cannot lead to an improved result.

The bounding step is one of the most critical operators in limiting the size of the search space. A well-designed bounding function can drastically reduce the number of explored nodes by efficiently identifying and pruning subproblems that cannot lead to an optimal solution. However, designing such a function is often a trade-off between tightness (how close the bound is to the true minimum) and computational cost. Tighter bounds typically enable more pruning, but may be more expensive to compute. In the following, we provide further detail on the bounding operator in the context of this work.

### 3.1. Bounding operator for qubit allocation

For our B&B algorithm applied to the qubit allocation problem, we need to be able to produce lower bound estimations of the route cost function $RC$ (2). To this end, we adopt the Hightower-Hahn Bound (HHB) introduced in [17], known for its high precision in estimating lower bounds for QAP instances. The HHB relies on a cost representation that consolidates all pairwise interaction and distance costs into a four-dimensional structure called the cost tensor.

<footnote n="1">¹ Without loss of generality the minimization case is considered in this work.</footnote>
$$
\mathbb{C} = \left[ \begin{bmatrix} C_{00} & \overline{C_{01}} & \overline{C_{02}} \\ \overline{C_{10}} & C_{11} & \overline{C_{12}} \\ \overline{C_{20}} & \overline{C_{21}} & C_{22} \end{bmatrix} \right] = \left[ \begin{array}{cccccccccc} c_{0000} & * & * & * & * & c_{0101} & * & * & * & c_{0202} \\ * & c_{0011} & c_{0012} & * & * & c_{0110} & c_{0112} & c_{0201} & c_{0211} & * \\ * & * & * & c_{0021} & - & c_{0022} & c_{0122} & c_{0220} & c_{0221} & * \\ * & * & c_{1001} & c_{1002} & - & c_{1100} & c_{1102} & c_{1200} & c_{1201} & * \\ c_{1010} & * & * & * & c_{1111} & * & * & * & c_{1212} & * \\ * & * & - & c_{1021} & c_{1022} & c_{1120} & c_{1122} & c_{1220} & c_{1221} & * \\ * & * & * & c_{2001} & c_{2002} & c_{2100} & c_{2102} & c_{2200} & c_{2201} & * \\ * & c_{2011} & c_{2012} & c_{2110} & c_{2112} & c_{2210} & c_{2211} & c_{2211} & c_{2211} & * \\ c_{2020} & * & * & * & c_{2121} & * & * & * & * & c_{2222} \end{array} \right]
$$

<caption>Fig. 7. Cost matrix representation (N = 3).</caption>

$$
\mathbb{L} = \begin{bmatrix} c_{0000} & c_{0101} & c_{0202} \\ c_{1010} & c_{1111} & c_{1212} \\ c_{2020} & c_{2121} & c_{2222} \end{bmatrix}
$$

<caption>Fig. 8. Leader matrix representation (N = 3).</caption>

**Definition 4 (Cost tensor).** Given an interaction frequency matrix $\mathbb{F} \in \mathcal{M}_N(\mathbb{N})$ and a coupling distance matrix $\mathbb{D} \in \mathcal{M}_N(\mathbb{N})$, the cost tensor $\mathbb{C} = (c_{ijk\ell})$ is defined for all $(i, j, k, \ell) \in \{0, \dots, N-1\}^4$ as:

$$ c_{ijk\ell} = f_{ik} d_{j\ell}. \qquad (3) $$

Each entry $c_{ijk\ell}$ quantifies the routing cost associated with assigning logical qubits $q_i$ and $q_k$ to physical qubits $P_j$ and $P_\ell$, respectively. This formulation captures the quadratic nature of the cost function.

The tensor $\mathbb{C}$ satisfies two key properties:

$$ \begin{gathered} \forall p \in \{0, \dots, N-1\}, i \neq k, j \neq \ell, \\ c_{ipkp} \text{ and } c_{pj\ell} \text{ are undefined elements,} \end{gathered} \qquad (4) $$

$$ \forall (i, j, k, \ell) \in \{0, N-1\}^4, c_{ijk\ell} = c_{k\ell ij}. \qquad (5) $$

Property (4) enforces the one-to-one mapping constraint: each logical qubit must be assigned to a unique physical qubit, and vice versa. The associated costs that would violate this constraint are disregarded.

Property (5) follows from the symmetry of $\mathbb{F}$ and $\mathbb{D}$, ensuring that the cost contributions are symmetric with respect to the interacting qubit pairs.

The tensor $\mathbb{C}$ can be conceptually reorganized into a large cost matrix of size $N^2 \times N^2$, consisting of $N^2$ sub-matrices $C_{ij}$ of dimension $N \times N$. Each sub-matrix $C_{ij}$ corresponds to a particular interaction between logical qubit $q_i$ assigned to physical qubit $P_j$ and all other possible assignments of the remaining qubits. We will henceforth refer to $\mathbb{C}$ as the cost matrix rather than the cost tensor.

In each sub-matrix $C_{ij}$, the element $c_{ijij}$, referred to as the leader element, serves as the diagonal anchor. Collecting all leader elements across the $N^2$ sub-matrices gives rise to a new matrix: $\mathbb{L} = (\ell_{ij})_{0 \le i, j \le N} = (c_{ijij})_{0 \le i, j \le N}$ which we refer to as the leader matrix. This matrix plays a central role in the computation of the HHB lower bound. Fig. 7 illustrates, for the case $N=3$, how $\mathbb{C}$ is structured as a block matrix, where undefined elements are marked as asterisks. In practice, these entries are set to $+\infty$, representing unmanageable costs. Fig. 8 shows the corresponding leader matrix.

We now describe the procedure for computing the HHB, which unfolds in three main steps. Each step incrementally transforms the original cost matrix $\mathbb{C}$, extracting numerical information that enables the derivation of tight lower bound estimates.

* **Distributing the leaders.** For each sub-matrix $C_{ij}$, distribute the leader cost $c_{ijij}$ by incrementing each allowed element in $C_{ij}$ by $c_{ijij}/(N-1)$, excluding the leader itself.

* **Halve complementary pairs.** In order to maintain symmetry in the cost matrix $\mathbb{C}$, all complementary pairs are equaled. This means that for all $(i, j, k, \ell) \in [0, N-1]^4$, $c_{ijk\ell}$ and $c_{k\ell ij}$ are set to $(c_{ijk\ell} + c_{k\ell ij})/2$. If the sum is odd, the rounding ambiguity is resolved deterministically: the entry with an even total index $i + j + k + \ell$ receives the extra unit, preserving numerical balance across the matrix.

* **Apply Hungarian Algorithm to the sub-matrices.** The Hungarian Algorithm (HA), introduced by Kuhn [18] and later refined by Munkres [19], makes up the core of the HHB lower bound procedure. Each sub-matrix $C_{ij}$ and the leader matrix $\mathbb{L}$ can be treated as a Linear Assignment Problem (LAP) instance of size $N$, to which the HA is applied. The algorithm outputs a cost corresponding to the optimal solution of the LAP by successively performing row and column reductions and identifying a minimum-cost perfect matching, modifying the input $N \times N$ matrix in the process. In our case, these matrices may contain undefined entries, which are implicitly handled since they were previously set to $+\infty$. Solving these LAP instances yields incremental lower-bound contributions that collectively tighten the bound for the full QAP.

The combination of the three previously described steps, leader distribution, symmetry enforcement, and linear assignment resolution, defines the complete lower-bound operator for the route cost function in the QAP formulation (1). This operator is formalized in the pseudo-code presented in Algorithm 1. The procedure is iterative, and the three main steps are repeated for a maximum of $it\_max$ iterations, a configurable parameter of the approach. A larger $it\_max$ can improve the refinement of the lower-bound estimate in cases where convergence is slow (i.e., when $incre \neq 0$ for many iterations), though at the cost of increased computational effort.

<table><caption>Algorithm 1: Lower-bound procedure.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>C</td><td>cost matrix</td></tr><tr><td>L</td><td>leader matrix</td></tr><tr><td>it_max</td><td>iteration threshold</td></tr><tr><td>lb</td><td>0, it &lt; 0</td></tr><tr><td colspan="2">while it &lt; it_max do</td></tr><tr><td></td><td>it &lt; it + 1</td></tr><tr><td></td><td>DistributeLeader(C)</td></tr><tr><td></td><td>HalveComplementary(C)</td></tr><tr><td></td><td>for each sub-matrix C<sub>ij</sub> in C do</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;cost &lt; Hungarian(C<sub>ij</sub>)</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;c<sub>ijij</sub> &lt; c<sub>ijij</sub> + cost</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;ℓ<sub>ij</sub> &lt; ℓ<sub>ij</sub> + cost</td></tr><tr><td></td><td>end</td></tr><tr><td></td><td>incre &lt; Hungarian(L)</td></tr><tr><td></td><td>if incre = 0 then</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;| break</td></tr><tr><td></td><td>end</td></tr><tr><td></td><td>lb &lt; lb + incre</td></tr><tr><td></td><td>end</td></tr><tr><td></td><td>return lb</td></tr></tbody></table>
Convergence can be detected early: if the total cost increment *incre* becomes zero at the end of an iteration, the procedure is considered to have converged and terminates before reaching the iteration limit. In practice, it is also terminated early whenever the estimated lower bound exceeds the current best known cost during the B&B search, thereby avoiding unnecessary computations.

Applying the HA to all sub-matrices constitutes the most computationally intensive step of the bounding procedure. Since the cost matrix $C$ is composed of $N^2$ sub-matrices and the HA has a time complexity of $\mathcal{O}(N^3)$, the overall time complexity of Algorithm 1 is $\mathcal{O}(N^5)$, assuming the maximum number of iterations $it\_max$ is bounded.

Initially, Algorithm 1 is applied to the complete cost matrix $C$ to compute a lower bound for the entire problem instance. Subsequently, during the B&B search, this lower-bound operator is applied to each node of the permutation tree. For these intermediate evaluations, the cost matrix $C$ is modified (or reduced) according to the current partial assignment.

## 3.2. Branch-and-bound for qubit allocation

We now present the design of the B&B algorithm for solving Eq. (1). The algorithm adopts a pool-based strategy, in which the nodes of the search tree are maintained in a data structure referred to as the *pool*. This choice is particularly well suited to parallelization, providing a natural basis for extending the sequential algorithm into a parallel one. The pool orchestrates the traversal of the state-space tree, where each node represents a partial permutation: a partial assignment of logical qubits to physical qubits. A *leaf node* corresponds to a complete allocation of all $n$ logical qubits, evaluated using the routing cost function *RC* defined in Eq. (2). Throughout the search, the algorithm maintains the best solution found so far, with a cost denoted *min_cost*.

During exploration, the algorithm incrementally builds candidate allocations by expanding each internal node into its associated *child nodes*, effectively assigning one additional logical qubit to a physical qubit at each step. A reduction operator is then applied to shrink the dimensionality of the cost matrix $C$, producing a reduced subproblem consistent with the new assignment constraints.

A crucial design choice is the order in which logical qubits are allocated. Following Zhu et al.'s [13] approach, we prioritize qubits according to their interaction activity, defined as the total number of two-qubit gates in which each qubit participates. This quantity is computed from the interaction matrix $F$ by summing the weights of all two-qubit interactions involving the given qubit. Logical qubits with the highest activity are placed earlier in the allocation sequence, so that their assignments more quickly tighten the lower bounds. This heuristic increases the likelihood of early pruning and significantly reduces the number of subproblems explored.

To control the search effort, the algorithm employs two pruning mechanisms:

* **Bounding-based pruning.** For every internal child, the bounding operator from Algorithm 1 estimates the minimal achievable cost of extending the partial solution. If this bound exceeds the current *min_cost*, the child is discarded.

* **Decomposition-based pruning.** During expansion, certain assignments are eliminated immediately by incrementing the lower bound using leader elements in the cost matrix, preventing unpromising branches from being explored.

The surviving child nodes are inserted into the pool for further exploration. The pool is implemented as a stack operating in Last-In-First-Out (LIFO) order, corresponding to a Depth-First Search (DFS) of the state-space tree. This design improves memory efficiency and often accelerates the discovery of complete candidate solutions.

To accelerate convergence, the algorithm is initialized with a greedy solution obtained using the heuristic of [13]. This procedure sequentially allocates each logical qubit by selecting the physical qubit that

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Overall workflow of the parallel B&B algorithm in a distributed setting across $K$ compute nodes, each leveraging $p$ CPU cores. Each computing thread $t_i$ maintains a private pool of tree nodes and concurrently performs branching and pruning operations. Intra-node and inter-node communication mechanisms used for dynamic load balancing and information sharing are also illustrated.</figcaption></figure>

yields the smallest increase in the routing cost function *RC*. The cost of the resulting mapping defines the initial value of *min_cost*, enabling early pruning of suboptimal nodes during the search.

The algorithm terminates once all nodes have been either explored or discarded. The final output is the optimal qubit allocation(s) with the lowest routing cost *min_cost*, representing the most favorable solution(s) for the given problem instance.

## 3.3. Parallel design

The design of parallel B&B algorithms is widely studied in the literature [20,21]. One of the most successful strategies, known for its generality and high degree of parallelism, is called the *parallel tree exploration* model. It consists of exploring several disjoint sub-trees concurrently and asynchronously in parallel on different computing threads. The design challenges raised by this model are the irregularity between B&B workers, which usually requires efficient load balancing strategies, global termination detection, and global knowledge sharing (e.g., the cost of the best solution found so far).

In this work, the B&B algorithm is parallelized using the P3D-DFS framework [7], implemented in the Chapel parallel programming language.<fnref n="2" /> The framework provides several levels of abstraction for implementing parallel B&B algorithms on CPU-based systems. The first is the *Node type*, representing subproblems. It typically stores the data needed to compute the bound for a subproblem—such as the (reduced) $C$ and $\ell$ matrices—and also contains the partial solution associated with that subproblem. The second level of abstraction is dedicated to the definition of the problem itself. This includes the implementation of B&B operators, such as the bounding procedure defined in Section 3.1.

Once the *Node type* and the *Problem class* are implemented, the framework enables the use of different versions of the B&B scheme: a sequential version, a multi-core parallel version, and a distributed version. Parallelism is managed through a massively parallel and decentralized data structure, where each computing thread maintains its own local work pool. Load balancing is achieved through a Work-Stealing (WS) mechanism, while global consistency is ensured by an advanced synchronization scheme. This combination enhances the scalability of the approach, making it suitable for architectures ranging from single multi-core machines to distributed clusters.

Fig. 9 illustrates the overall workflow of the presented parallel B&B algorithm.

<footnote n="2">² The parallel B&B design is general and can be adapted to other environments such as MPI, OpenMP, or C++ threads, ensuring broad applicability beyond Chapel.</footnote>
# 4. Experimental results

## 4.1. Experimental protocol and setup

In the following, two B&B implementations are evaluated:

* SBB. A sequential implementation that primarily serves as a baseline for evaluating the parallel variants;

* PBB. A parallel implementation featuring two levels of parallelism: intra-node and inter-node.

The experimental protocol is organized as follows. First, we verify the accuracy and efficiency of our approach by comparing results with those reported in the literature on previously solved instances. We then evaluate the scalability of the parallel implementation across both levels of parallelism. Finally, we test the limits of PBB by targeting the exact solution of large-scale instances. All reported execution times are averaged over five independent runs.

The experiments are conducted on an HPC cluster, using up to 64 compute nodes, where each node is equipped with two 64-core AMD Epyc ROME 7H12 @ 2.6 GHz CPUs and 256 GB of RAM. The nodes are interconnected via a fast InfiniBand HDR100 network configured in a Fat-Tree topology and run Red Hat Linux 8.10. B&B implementations are based on Chapel 2.5.0 in a fine-tuned configuration environment, using GCC 13.2.0 as the back-end compiler.

## 4.2. Comparison against state-of-the-art

We begin our experimental evaluation by reproducing the results obtained in the study by Zhu et al. [13], which represents the current state-of-the-art for the considered problem formulation. To this end, we use our sequential B&B baseline implementation SBB for comparison. While both approaches are conceptually similar, the implementation in Zhu et al. [13], referred to as qa_opt, is written in C++ and is based on a recursive B&B formulation and employs a modified HHB-based bounding operator designed to reduce the number of levels explored in the state-space tree. In contrast, as detailed in Section 3, our approach relies on a pool-based strategy and uses the regular HHB as the lower-bounding procedure.

The evaluation consists of solving the qubit allocation problem (1) for benchmark circuits sourced from RevLib [22], targeting the 16-qubit IBM Melbourne NISQ architecture, whose connectivity graph is shown in Fig. 10. Table 1 reports the results for each benchmark, including the number of logical qubits *n* and the optimal routing cost obtained. We also compare the number of visited nodes in the permutation tree and the normalized execution time of our B&B algorithm with the corresponding results reported in [13] using qa_opt. To ensure a fair comparison, execution times are normalized by CPU frequency, corresponding to 2.5 × time for qa_opt (2.5 GHz Intel i7-4710 in [13]) and 2.6 × time for SBB (2.6 GHz CPU, cf. Section 4.1).

Across all benchmarks, our execution times are consistently much lower, with most instances being solved in less than five seconds (considering the raw execution times). The exception is the Quantum Fourier Transform (QFT) circuit (listed as qft_10 in Table 1), which highlights a limitation of our implementation. Due to its highly regular structure, where every logical qubit interacts equally with all others, any permutation of physical qubits yields the same routing cost. This greatly reduces pruning effectiveness, since many subproblems appear equally promising and must be explored. Zhu et al. [13] address this challenge through a symmetry pruning mechanism, which likely accounts for the pronounced gap in execution times observed for this case. In fact, this symmetry test is computationally heavy, as it must be performed at every node of the B&B tree, and its benefits are essentially restricted to the QFT instance. For all other benchmarks—20 out of the 21 considered—our approach, unburdened by this additional overhead, achieves substantially lower runtimes. In a sense, we have alleviated the B&B algo-

<figure><img src="image_10.png" /><figcaption>Fig. 10. Connectivity graph of the IBM Melbourne [5] quantum device featuring *N* = 16 physical qubits.</figcaption></figure>

<caption>Table 1</caption>

Comparison of sequential branch-and-bound results (SBB) with Zhu et al.'s qa_opt approach on benchmark circuits up to 16 qubits targeting the IBM Melbourne architecture (cf. Fig. 10). Optimal costs, number of visited nodes, and normalized runtimes are reported. Best runtimes are shown in bold.

<table><thead><tr><th rowspan="2">circuit</th><th rowspan="2">n</th><th rowspan="2">optimal cost</th><th colspan="2">visited nodes</th><th colspan="2">normalized time</th></tr><tr><th>qa_opt</th><th>SBB</th><th>qa_opt</th><th>SBB</th></tr></thead><tbody><tr><td>qft_10</td><td>10</td><td>240</td><td>6280</td><td>49615236</td><td>22.0</td><td>28469.0</td></tr><tr><td>sqn_258</td><td>10</td><td>6140</td><td>16074</td><td>235</td><td>1642.8</td><td>3.1</td></tr><tr><td>sym9_148</td><td>10</td><td>13904</td><td>10324</td><td>208</td><td>1357.3</td><td>1.8</td></tr><tr><td>sym9_193</td><td>11</td><td>23936</td><td>6817</td><td>267</td><td>2644.0</td><td>4.7</td></tr><tr><td>wim_266</td><td>11</td><td>480</td><td>12493</td><td>644</td><td>658.0</td><td>4.9</td></tr><tr><td>z4_268</td><td>11</td><td>2140</td><td>36743</td><td>770</td><td>2623.8</td><td>15.9</td></tr><tr><td>cycle10_2_110</td><td>12</td><td>4688</td><td>2004</td><td>356</td><td>264.3</td><td>7.5</td></tr><tr><td>rd84_253</td><td>12</td><td>9864</td><td>1115</td><td>213</td><td>288.0</td><td>3.1</td></tr><tr><td>sym9_146</td><td>12</td><td>196</td><td>3103</td><td>1153</td><td>99.3</td><td>12.7</td></tr><tr><td>dist_223</td><td>13</td><td>31682</td><td>1250</td><td>481</td><td>732.5</td><td>9.9</td></tr><tr><td>radd_250</td><td>13</td><td>2484</td><td>2982</td><td>475</td><td>338.8</td><td>7.3</td></tr><tr><td>root_255</td><td>13</td><td>13514</td><td>1381</td><td>457</td><td>611.5</td><td>9.6</td></tr><tr><td>clip_206</td><td>14</td><td>33412</td><td>4453</td><td>1400</td><td>1946.3</td><td>30.4</td></tr><tr><td>cm42a_207</td><td>14</td><td>836</td><td>1323</td><td>2460</td><td>150.5</td><td>3.8</td></tr><tr><td>cm85a_209</td><td>14</td><td>10678</td><td>901</td><td>709</td><td>341.5</td><td>14.6</td></tr><tr><td>co14_215</td><td>15</td><td>10824</td><td>33</td><td>323</td><td>37.0</td><td>2.3</td></tr><tr><td>misex1_241</td><td>15</td><td>3188</td><td>335</td><td>517</td><td>110.8</td><td>4.2</td></tr><tr><td>sqrt_7</td><td>15</td><td>3072</td><td>725</td><td>437</td><td>516.8</td><td>9.1</td></tr><tr><td>inc_237</td><td>16</td><td>8318</td><td>100</td><td>379</td><td>41.3</td><td>3.1</td></tr><tr><td>ising_model_16</td><td>16</td><td>0</td><td>227</td><td>1380</td><td>2.8</td><td>0.8</td></tr><tr><td>mlp4_245</td><td>16</td><td>21408</td><td>1435</td><td>2429</td><td>827.8</td><td>53.8</td></tr></tbody></table>

rithm, setting state-of-the-art performance for the vast majority of the benchmark suite.

These runtime comparisons should, however, be interpreted with care, as the two implementations differ not only in programming environments and hardware platforms but also in their algorithmic strategies. In both cases, execution time is tightly linked to the number of explored nodes, since each node requires an expensive lower-bound evaluation (cf. Algorithm 1), carried out here with a convergence threshold of *it_max* = 10. Ultimately, pruning effectiveness remains the decisive factor in performance.

## 4.3. Parallel branch-and-bound performance evaluation

In this section, we evaluate the performance of our parallel B&B implementation. To do so, we target the 20-qubit IBM Boeblingen NISQ architecture, whose connectivity graph is shown in Fig. 3, and consider random circuits from the MQT Bench library [23] as benchmarks. Specifically, we consider circuits with logical qubit counts ranging from *n* = 11 to 20, denoted rd_11 to rd_20. These circuits yield particularly challenging QAP instances due to their dense and relatively uniform interaction patterns between logical qubits. We first provide a calibration of the bounding parameter *it_max*, followed by a scalability analysis of PBB considering both intra-node and inter-node parallel levels.

### 4.3.1. Bounding parameter calibration

Calibration is performed with the sequential B&B implementation SBB, evaluated on each instance with *it_max* values ranging from 2 to 6. The corresponding execution times are reported in Table 2, spanning from several minutes to multiple days of computation, and thus reflecting both the breadth of the benchmark and the challenging nature of
<caption>**Table 2**</caption>

Impact of the bounding iteration threshold *it_max* on SBB execution times (in seconds). For each circuit, the best performance is highlighted in bold.

<table><thead><tr><th rowspan="2">circuit</th><th colspan="5">it_max</th></tr><tr><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>rd_11</td><td>1092.7</td><td><strong>1074.2</strong></td><td>1091.3</td><td>1110.9</td><td>1111.4</td></tr><tr><td>rd_12</td><td>1752.8</td><td><strong>1635.9</strong></td><td>1654.9</td><td>1699.2</td><td>1720.8</td></tr><tr><td>rd_13</td><td>1971.7</td><td><strong>1881.5</strong></td><td>1920.0</td><td>1988.8</td><td>2024.8</td></tr><tr><td>rd_14</td><td>4602.2</td><td>4247.0</td><td><strong>4238.9</strong></td><td>4360.1</td><td>4401.1</td></tr><tr><td>rd_15</td><td>14859.8</td><td>12907.8</td><td><strong>12522.8</strong></td><td>12830.7</td><td>13042.6</td></tr><tr><td>rd_16</td><td>7161.5</td><td><strong>6511.5</strong></td><td>6543.6</td><td>6798.5</td><td>6952.5</td></tr><tr><td>rd_17</td><td>36258.8</td><td>31823.0</td><td><strong>31521.9</strong></td><td>32477.6</td><td>33332.0</td></tr><tr><td>rd_18</td><td>33548.9</td><td><strong>30738.4</strong></td><td>30936.4</td><td>32006.0</td><td>33016.8</td></tr><tr><td>rd_19</td><td>223996.0</td><td><strong>217967.5</strong></td><td>225757.0</td><td>239942.0</td><td>249269.0</td></tr><tr><td>rd_20</td><td>138682.4</td><td><strong>133001.4</strong></td><td>136525.2</td><td>141992.8</td><td>148010.4</td></tr></tbody></table>

the instances. On average, execution time is minimized at *it_max* = 3. Larger values offer no further benefit, as the reduction in visited nodes is insufficient to offset the increased cost of the bounding operator, which typically requires more iterations when *it_max* grows. For example, setting *it_max* = 6 on the rd_20 instance increases execution time by nearly 10% compared to the optimal setting, adding roughly five hours of computation. Consequently, we retain *it_max* = 3 for the remainder of this section. However, this value might not be optimal for all classes of circuits, as the best calibration can depend on the structural characteristics of each benchmark instance. Incorporating adaptive mechanisms for dynamically tuning *it_max*—for example, based on instance-specific features—represents a promising direction for future work.

### 4.3.2. *Scalability and performance evaluation*

We now evaluate the performance of our parallel B&B algorithm, PBB, beginning with its strong-scaling behavior at both the intra-node and inter-node levels.

To ensure fairness, each experiment is performed with the initial upper bound *min_cost* set to the known optimal routing cost. For instance, in the case of rd_11, *min_cost* was initialized to 232, rather than to the heuristic estimate used in previous experiments. This choice enforces exploration of the critical B&B tree, i.e., the minimal search tree required to prove optimality. Under this configuration, the number of visited nodes is identical for both SBB and PBB, regardless of the parallel execution strategy, thus enabling a consistent comparison.

Fig. 11a presents the intra-node scalability results, showing the relative speed-ups obtained when executing PBB on a single compute node in a multi-core setting. Three benchmark instances are considered: the qubit allocation of the random circuits rd_12, rd_14, and rd_16, all targeting the IBM Boeblingen architecture (cf. Fig. 3). These circuits were selected as representative moderately challenging instances of increasing size. The intra-node speed-up values are reported for an increasing number of CPU cores, relative to the sequential baseline SBB. Fig. 11a highlights the strong intra-node scalability of our parallel B&B implementation PBB. The observed speed-ups remain consistently close to the ideal values, even as the number of CPU cores increases. In all cases, when utilizing the full capacity of a node with 128 cores, PBB achieves more than 87% of the linear speed-up, underscoring the efficiency of the intra-node parallelization strategy.

Fig. 11b presents the inter-node scalability results, showing the relative speed-ups obtained when executing PBB in a distributed setting across multiple compute nodes, each leveraging 128 CPU cores. Three benchmark instances are considered: the qubit allocation of the random circuits rd_18, rd_20, and rd_19, all targeting the IBM Boeblingen architecture (cf. Fig. 3). These circuits were selected as representative challenging instances of increasing size. The inter-node speed-up values are reported for an increasing number of nodes, relative to the baseline performance of PBB on a single node with 128 CPU cores. In this

<caption>**Table 3**</caption>

Profiling results at the inter-node level using 64 compute nodes, summarizing the distribution of total execution time across key components for each benchmark circuit.

<table><caption>**Table 4**</caption><thead><tr><th rowspan="2">circuit</th><th colspan="4">% of total execution time</th></tr><tr><th>initialization</th><th>B&B operators</th><th>termination</th><th>WS</th></tr></thead><tbody><tr><td>rd_18</td><td>36%</td><td>48%</td><td>1%</td><td>9%</td></tr><tr><td>rd_20</td><td>12%</td><td>65%</td><td>&lt;1%</td><td>17%</td></tr><tr><td>rd_19</td><td>8%</td><td>79%</td><td>&lt;1%</td><td>11%</td></tr></tbody></table>

Results of SBB and PBB on random circuits rd_11 to rd_20 targeting the IBM Boeblingen architecture (cf. Fig. 3). For each instance, optimal cost, number of visited nodes, runtime (in seconds), and achieved speed-up are reported.

<table><thead><tr><th rowspan="2">circuit</th><th rowspan="2">n</th><th rowspan="2">optimal cost</th><th rowspan="2">visited nodes</th><th colspan="2">time (s)</th><th rowspan="2">speed-up</th></tr><tr><th>SBB</th><th>PBB</th></tr></thead><tbody><tr><td>rd_11</td><td>11</td><td>232</td><td>103 008</td><td>1033.3</td><td>5.1</td><td>x202</td></tr><tr><td>rd_12</td><td>12</td><td>322</td><td>122 002</td><td>1595.5</td><td>5.0</td><td>x319</td></tr><tr><td>rd_13</td><td>13</td><td>398</td><td>68 870</td><td>1135.4</td><td>5.2</td><td>x218</td></tr><tr><td>rd_14</td><td>14</td><td>428</td><td>341 368</td><td>4109.6</td><td>5.6</td><td>x733</td></tr><tr><td>rd_15</td><td>15</td><td>682</td><td>702 928</td><td>9398.3</td><td>6.6</td><td>x1423</td></tr><tr><td>rd_16</td><td>16</td><td>822</td><td>424 674</td><td>5618.9</td><td>5.8</td><td>x968</td></tr><tr><td>rd_17</td><td>17</td><td>826</td><td>2 530 810</td><td>30 917.5</td><td>10.3</td><td>x3001</td></tr><tr><td>rd_18</td><td>18</td><td>1140</td><td>2 507 128</td><td>29 569.1</td><td>9.9</td><td>x2986</td></tr><tr><td>rd_19</td><td>19</td><td>1400</td><td>24 465 752</td><td>200 456.0</td><td>35.2</td><td>x5694</td></tr><tr><td>rd_20</td><td>20</td><td>1582</td><td>10 693 810</td><td>98 856.0</td><td>20.4</td><td>x4845</td></tr></tbody></table>

configuration, more noticeable limitations arise compared to the intra-node setting. Nonetheless, scalability remains robust, particularly for the most demanding instances, reaching up to 74% of the ideal speed-up for rd_19 at maximum scale—64 compute nodes, totaling 64 × 128 = 8192 CPU cores—highlighting the efficiency of the distributed parallelization strategy.

The profiling data summarized in Table 3 provide further insight into the inter-node scalability trends observed in Fig. 11b. For smaller circuits such as rd_18, a significant fraction of total runtime is devoted to initialization (36%) and work-stealing (9%), indicating that communication and ramp-up overheads dominate at large node counts. As the problem size increases (rd_19 and rd_20), the relative cost of these phases decreases, while the proportion of time spent in B&B operators grows markedly (up to 79%). This shift reflects higher computational intensity per node, which mitigates the effects of communication latency and workload imbalance. Overall, the profiling confirms that inter-node scalability is constrained mainly by initialization and synchronization overheads for smaller benchmarks, whereas larger circuits are increasingly computation-bound, leading to improved parallel efficiency.

Having exhibited the scaling capabilities of PBB, we now turn to its overall performance on the qubit allocation of the random benchmark circuits rd_11 to rd_20. Table 4 reports, for each circuit, the optimal cost, the number of visited nodes in the B&B tree, and the execution times of both SBB and PBB. The PBB runs were executed in a distributed setting on 64 compute nodes, each leveraging 128 CPU cores. The table also lists the speed-ups obtained by PBB over its sequential counterpart SBB. Together, these results highlight the exponential growth of computational effort in the sequential approach and the substantial performance gains enabled by our parallel implementation. They further demonstrate PBB’s ability to effectively exploit the resources of modern HPC clusters, particularly for larger problem instances.

### 4.4. *Large-scale experiments*

We conclude our experimental evaluation by testing the limits of our parallel B&B implementation PBB, aiming to obtain exact qubit
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Strong-scaling evaluation of PBB for various benchmark instances on the IBM Boeblingen architecture (cf. Fig. 3), considering both levels of parallelism. (a) Intra-node strong-scaling performance of PBB. Speed-ups are measured relative to the sequential baseline SBB. (b) Inter-node strong-scaling performance of PBB. Speed-ups are measured relative to multi-core PBB using 128 CPU cores.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Connectivity graph of the IBM Cairo [5] quantum device featuring $N = 27$ physical qubits.</figcaption></figure>

allocation solutions for larger problem instances. Specifically, we target the 27-qubit IBM Cairo architecture, whose connectivity graph is shown in Fig. 12. As benchmarks, we use real quantum circuits from the RevLib [22] library, ranging from 17 to 26 logical qubits, decomposed with the Real2QASM framework [10,22,24,25] into the IBM-style basis gate set (u3, CX). For each benchmark circuit, Table 5 reports the results obtained with PBB executed in a distributed setting on 64 compute nodes, each leveraging 128 CPU cores. The table provides the optimal routing costs, the number of visited B&B tree nodes, and the corresponding execution times.

The results confirm the ability of our approach to deliver exact qubit allocation solutions for circuits of up to 26 qubits, with execution times ranging from a few tens of seconds to several hours depending on the benchmark. The size of the explored B&B trees varies widely, from $1.5 \times 10^5$ to $1.1 \times 10^9$ nodes, illustrating the exponential growth in problem complexity. The results also suggest that the difficulty of solving a given instance is influenced not only by its size (i.e., the number of qubits) but also by intrinsic circuit characteristics, making benchmark landscaping an interesting direction for future analysis.

For the largest instance, rd73_312, the parallel execution took about 6795 seconds. Assuming a parallel speed-up comparable to the best observed on smaller benchmarks (up to $\times 5694$, cf. Table 4), this implies a sequential runtime lower bound of $5694 \times 6795$ seconds, i.e., more than one year of computation. This estimate underscores the prohibitive cost of solving such instances sequentially.

Overall, these findings demonstrate that large-scale parallelization can successfully counterbalance the exponential growth of the

<caption>**Table 5**</caption>

Results of PBB on RevLib benchmark circuits targeting the 27-qubit IBM Cairo architecture (cf. Fig. 12). For each instance, the optimal cost, number of visited nodes, and runtime (in seconds) are reported.

<table><thead><tr><td>circuit</td><td>n</td><td>optimal cost</td><td>visited nodes</td><td>time (s)</td></tr></thead><tbody><tr><td>5xp1_194</td><td>17</td><td>3034</td><td>448 319</td><td>36.9</td></tr><tr><td>ryy6_256</td><td>17</td><td>14 942</td><td>2 005 834</td><td>123.1</td></tr><tr><td>alu3_200</td><td>18</td><td>8306</td><td>1 112 126</td><td>48.4</td></tr><tr><td>sqrt6_259</td><td>18</td><td>2222</td><td>239 359</td><td>29.1</td></tr><tr><td>add6_196</td><td>19</td><td>15 408</td><td>306 030</td><td>20.0</td></tr><tr><td>cmb_214</td><td>20</td><td>4428</td><td>3 060 927</td><td>46.9</td></tr><tr><td>ex1010_230</td><td>20</td><td>224 476</td><td>156 164</td><td>21.3</td></tr><tr><td>decod_217</td><td>21</td><td>3644</td><td>448 807</td><td>14.7</td></tr><tr><td>dk17_224</td><td>21</td><td>3600</td><td>710 611</td><td>26.1</td></tr><tr><td>apla_203</td><td>22</td><td>8004</td><td>492 185</td><td>31.7</td></tr><tr><td>cm105a_210</td><td>22</td><td>1780</td><td>4 965 734</td><td>52.3</td></tr><tr><td>plus63mod4096</td><td>23</td><td>52</td><td>32 171 779</td><td>88.7</td></tr><tr><td>arb8_323</td><td>24</td><td>1880</td><td>152 951 233</td><td>1363.7</td></tr><tr><td>cu_219</td><td>25</td><td>2652</td><td>10 131 120</td><td>335.1</td></tr><tr><td>rd73_312</td><td>25</td><td>654</td><td>1 110 112 856</td><td>6795.0</td></tr><tr><td>in0_235</td><td>26</td><td>44 912</td><td>1 798 354</td><td>60.8</td></tr></tbody></table>

search space, enabling exact qubit allocation for circuits of up to 26 logical qubits—instances that would otherwise be computationally intractable.
# 5. Conclusions and perspectives

In this work, we addressed the scalability challenge of exact qubit allocation by revisiting its formulation as a permutation-based QAP and solving it with a B&B algorithm. We first established a refined sequential implementation that achieves substantially faster runtimes than a previous exact approach on most problem instances, thereby setting a new state-of-the-art for this formulation. Building upon this foundation, we extended the approach to a performance-aware parallel version capable of effectively exploiting both intra-node and inter-node parallelism on modern HPC infrastructures. Our experimental evaluation demonstrated excellent strong-scaling behavior at the intra-node level, with speed-ups approaching ideal values, and substantial scalability in distributed settings, reaching up to 74% of the ideal performance across nodes. Large-scale experiments on a 27-qubit architecture further confirmed the ability of our method to compute exact solutions for benchmark circuits of up to 26 qubits, establishing new reference points for the community and showing that large-scale parallelization can effectively extend the reach of exact approaches beyond previously tractable problem sizes.

Looking ahead, several directions for future work emerge. First, scalability could be further enhanced by extending our implementation to a GPU-accelerated B&B algorithm capable of exploiting modern heterogeneous HPC clusters [26]. In this context, we could greatly benefit from GPU acceleration, as our results show that for large problem instances, the B&B operators increasingly dominate the total runtime compared to communication phases. Another promising direction is the exploration of multi-objective formulations for qubit allocation. In particular, incorporating error-aware criteria would allow the framework to account not only for connectivity constraints but also for hardware reliability—since, in practice, physical qubits on a device can exhibit varying error rates. Finally, combining our exact allocation approach with heuristic techniques designed for the dynamic routing phase [9,11] could further extend this work toward a complete and practically applicable qubit mapping framework.

## CRediT authorship contribution statement

**Jean-Philippe Valois:** Writing – review & editing, Writing – original draft, Visualization, Validation, Investigation, Formal analysis, Conceptualization; **Guillaume Helbecque:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Investigation, Conceptualization; **Nouredine Melab:** Writing – review & editing, Supervision, Project administration, Funding acquisition.

## Data availability

The code used for this article is publicly available at https://github.com/Guillaume-Helbecque/P3D-DFS.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Jean-Philippe Valois reports financial support was provided by the French National Research Agency. Guillaume Helbecque reports financial support was provided by the Luxembourg National Research Fund. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported by the Agence Nationale de la Recherche (ANR) [grant number ANR-22-CE46-0011] and the Luxembourg National Research Fund [grant number INTER/ANR/22/17133848] under the UltraBO project, and by the ANR [grant number ANR-21-ESRE-0051] under the EquipEx+ MesoNET project. The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg.

## References

[1] L.K. Grover, A fast quantum mechanical algorithm for database search, in: Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing, 1996, p. 212-219. https://doi.org/10.1145/237814.237866

[2] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P.J. Love, A. Aspuru-Guzik, J.L. O'Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5 (1) (2014). https://doi.org/10.1038/ncomms5213

[3] M. Schuld, I. Sinayskiy, F. Petruccione, The quest for a quantum neural network, Quantum Inf. Process. 13 (11) (2014) 2567-2586. https://doi.org/10.1007/s11128-014-0809-8

[4] P.W. Shor, Algorithms for quantum computation: discrete logarithms and factoring, in: Proceedings 35th Annual Symposium on Foundations of Computer Science, 1994, pp. 124-134. https://doi.org/10.1109/SFCS.1994.365700

[5] M. Bozzo-Rey, R. Loredo, Introduction to the IBM Q experience and quantum computing, in: Proceedings of the 28th Annual International Conference on Computer Science and Software Engineering, 2018, p. 410-412.

[6] S.S. Tannu, M.K. Qureshi, Not all qubits are created equal: a case for variability-aware policies for NISQ-era quantum computers, in: Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems, Association for Computing Machinery, 2019, p. 987-999. https://doi.org/10.1145/3297858.3304007

[7] G. Helbecque, J. Gmys, N. Melab, T. Carneiro, P. Bouvry, Parallel distributed productivity-aware tree search using chapel, Concurr. Comput. Pract. Exper. 35 (27) (2023) e7874. https://doi.org/10.1002/cpe.7874

[8] A. Barenco, C.H. Bennett, R. Cleve, D.P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J.A. Smolin, H. Weinfurter, Elementary gates for quantum computation, Phys. Rev. A 52 (5) (1995) 3457-3467. https://doi.org/10.1103/physreva.52.3457

[9] G. Li, Y. Ding, Y. Xie, Tackling the qubit mapping problem for NISQ-Era quantum devices, in: Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems, 2019, p. 1001-1014. https://doi.org/10.1145/3297858.3304023

[10] IBM, Qiskit, Open Source Quantum Information Science Kit, 2018, https://qiskit.org/.

[11] J. Rouzé, N. Melab, D. Tuytens, A parallel genetic algorithm for qubit mapping on noisy intermediate-scale quantum machines, in: Optimization and Learning, Optimization and Learning, 2025, p. 305-320. https://doi.org/10.1007/978-3-031-77941-1_24

[12] A. Zulehner, S. Gasser, R. Wille, Exact global reordering for nearest neighbor quantum circuits using A*, in: Reversible Computation, 2017, pp. 185-201. https://doi.org/10.1007/978-3-319-59936-6_15

[13] P. Zhu, C. Xueyun, Z. Guan, An exact qubit allocation approach for NISQ architectures, Quant. Inf. Process. 19 (2020). https://doi.org/10.1007/s11128-020-02901-4

[14] C. Zhang, A.B. Hayes, L. Qiu, Y. Jin, Y. Chen, E.Z. Zhang, Time-optimal qubit mapping, in: 26th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, 2021, pp. 360-374. https://doi.org/10.1145/3445814.3446706

[15] M.R. Garey, D.S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness, W. H. Freeman & Co., USA, USA, 1990.

[16] S. Sahni, T. Gonzalez, P-Complete Approximation problems, J. ACM 23 (3) (1976) 555-565. https://doi.org/10.1145/321958.321975

[17] P. Hahn, T. Grant, Lower bounds for the quadratic assignment problem based upon a dual formulation, Oper. Res. 46 (6) (1998) 912-922. https://doi.org/10.1287/opre.46.6.912

[18] H.W. Kuhn, The hungarian method for the assignment problem, Naval Research Logistics Quarterly 2 (1-2) (1955) 83-97. https://doi.org/10.1002/nav.3800020109

[19] J. Munkres, Algorithms for the assignment and transportation problems, J. Soc. Industrial. Appl. Math. 5 (1) (1957) 32-38. https://doi.org/10.1137/0105003

[20] B. Gendron, T.G. Grainic, Parallel branch-and-bound algorithms: survey and synthesis, Oper. Res. 42 (6) (1994) 1042-1066. https://doi.org/10.1287/opre.42.6.1042

[21] T.G. Grainic, B. Le Cun, C. Roucairol, Parallel Branch-and-Bound Algorithms, John Wiley & Sons, Ltd, 2006, pp. 1-28. https://doi.org/10.1002/9780470053928.ch1

[22] R. Wille, D. Große, L. Teuber, G.W. Dueck, R. Drechsler, RevLib: an online resource for reversible functions and reversible circuits, in: Proceedings of the 38th International Symposium on Multiple Valued Logic, 2008, pp. 220-225. https://doi.org/10.1109/ISMVL.2008.43

[23] N. Quetschlich, L. Burgholzer, R. Wille, MQT Bench: benchmarking software and design automation tools for quantum computing, Quantum 7 (1062) (2023) 1-14. https://doi.org/10.22331/q-2023-07-20-1062

[24] K.-Y. Chang, C.-Y. Lee, Mapping nearest neighbor compliant quantum circuits onto a 2-D hexagonal architecture, IEEE Trans. Comput. Aided Des. Integr. Circuits Syst. 41 (10) (2022) 3373-3386. https://doi.org/10.1109/TCAD.2021.3127868

[25] A.W. Cross, L.S. Bishop, J.A. Smolin, J.M. Gambetta, Open quantum assembly language, arXiv e-prints, arXiv:1707.03429, (2017), https://doi.org/10.48550/arXiv.1707.03429.

[26] G. Helbecque, PGAS-based Parallel Branch-and-Bound for Ultra-Scale GPU-powered Supercomputers, Ph.D. thesis, Université de Lille ; Université du Luxembourg, 2025.