Optimization of Quantum Systems Emulation via a Variant of the Bandwidth Minimization Problem
M. Yassine Naghmouchi Joseph Vovrosh Wesley da Silva Coelho Alexandre Dauphin











 

Optimization of Quantum Systems Emulation via a Variant of the Bandwidth Minimization Problem
    



1 M. Yassine Naghmouchi,
1 Joseph Vovrosh,
1 Wesley da Silva Coelho,
1 Alexandre Dauphin






1PASQAL SAS,7 rue Léonard de Vinci, 91300 Massy, France

    April 29, 2024
===============================================================================================================================================================

plain
plain


This paper introduces weighted-BMP, a variant of the Bandwidth Minimization Problem (BMP), with a significant application in optimizing quantum emulation. Weighted-BMP optimizes particles ordering to reduce the emulation costs, by designing a particle interaction matrix where strong interactions are placed as close as possible to the diagonal. We formulate the problem using a Mixed Integer Linear Program (MILP) and solve it to optimality with a state-of-the-art solver. To strengthen our MILP model, we introduce symmetry-breaking inequalities and establish a lower bound. Through extensive numerical analysis, we examine the impacts of these enhancements on the solver's performance. The introduced reinforcements result in an average CPU time reduction of 25.61%. Additionally, we conduct quantum emulations of realistic instances. Our numerical tests show that the weighted-BMP approach outperforms the Reverse Cuthill-McKee (RCM) algorithm—an efficient heuristic used for site ordering tasks in quantum emulation— achieving an average memory storage reduction of 24.48%. From an application standpoint, this study is the first to apply an exact optimization method, weighted-BMP, that considers interactions for site ordering in quantum emulation pre-processing, and shows its crucial role in cost reduction. From an algorithmic perspective, it contributes by introducing important reinforcements and lays the groundwork for future research on further enhancements, particularly on strengthening the weak linear relaxation of the MILP. 





Bandwidth minimization, tensor networks, quantum computing emulation.






§ INTRODUCTION


Given an undirected graph G = (V, E) with a set of vertices V and a set of edges E, the Bandwidth Minimization Problem (BMP) <cit.> is to find an ordering of the vertices that minimizes the maximum difference between the positions of adjacent vertices. Formally, let π: V →{1, 2, ..., |V|} be a bijection, where π(u) represents the position of vertex u in the ordering. The bandwidth, denoted as b, is defined as b = (u, v) ∈ Emax |π(u) - π(v)|.  The objective is to find a permutation π that minimizes b. In addition to its graph-based formulation, the BMP can be represented using matrices. Specifically, the objective is to reorder the rows and columns of a square matrix in a way that minimizes the largest distance between the positions of non-zero elements. This ordering minimizes the distance of the nonzero components from the main diagonal of the matrix. Achieving a minimal bandwidth in G is equivalent to minimizing the bandwidth in its adjacency matrix, where each vertex is assigned to the same row/column. 


In the field emulation of quantum many-body systems  <cit.>-involving numerous interacting sites- the BMP plays an essential role, particularly in improving the computational cost of handling Tensor Networks (TN) <cit.>. Tensors, multi-dimensional arrays that generalize scalars, vectors, and matrices, are fundamental in modeling quantum many-body systems and their inherent interactions. TNs provide a mathematical framework that use these tensors, connecting them to represent quantum states and operators. Matrix Product State (MPS), a specific case of TN, are used for representing quantum systems in a one-dimensional representation. MPS decompose quantum systems, initially represented by high-ranked tensors, into sequences of interconnected tensors, where each tensor corresponds to the quantum state of an individual site within the system. In quantum emulation, the conversion of a high-rank tensor into an MPS is a common pre-processing approach. This transformation is important as it computationally facilitates emulations that would otherwise be impractical, due to the exponential growth of the Hilbert space with the increasing number of sites.





Recent studies <cit.> highlight the benefits of efficiently mapping high-ranked tensors to MPS for quantum emulations, in terms of memory storage and computational time. An efficient strategy for mapping usually involves using the interaction matrix of the quantum many-body system, which details the strength and nature of interactions between particles in a quantum many-body system. The goal is to reorder its rows and columns in a way that minimizes the bandwidth, and then use this ordering for the MPS emulations. This ordering process is equivalent to address the BMP. Nonetheless, these studies mainly utilize the BMP in its classical definition, neglecting the magnitude of the interactions which might impact the efficiency of the ordering. Moreover, they often depend on heuristic methods such as the Reverse Cuthill-McKee (RCM) algorithm <cit.>, which, while reducing bandwidth, yield sub-optimal orderings that may not induce the best computational efficiency in emulations. Given the significance of incorporating site interactions in the ordering task, and the benefit of an optimal ordering through exact methods rather than heuristics, this paper introduces an exact resolution method for a BMP variant, named weighted-BMP, that includes interactions (weights). This variant, aims at finding an ordering of the rows and columns of the interaction matrix, such that elements corresponding to stronger interactions are as close as possible to the main diagonal.


This paper presents two main contributions: an algorithmic/optimization contribution and an application contribution. From an algorithmic perspective, we model the weighted-BMP using a Mixed Integer Linear Program (MILP), which we solve to optimality. To address the computational complexity of this problem, we incorporate reinforcements such as symmetry-breaking inequalities and establish a theoretical lower bound. Furthermore, we conduct extensive numerical analysis in order to evaluate the effect of these reinforcements on the solver's performance. On the application side, we perform quantum emulations for significant physical systems, in particular, amorphous solids <cit.>, in the context of a neutral atom processor <cit.>. We investigate the enhancement in quantum emulation efficiency through the optimal weighted-BMP ordering. Throughout this numerical study, we benchmark our weighted-BMP solutions against those obtained using the RCM method. From an application standpoint, the primary goal of this work is to underscore the role of weighted-BMP in quantum emulation. On the other hand, from an algorithmic viewpoint, our objective is to propose effective preliminary reinforcements, via symmetry-breaking inequalities and a lower bound.  

This paper is organized as follows: Section <ref> reviews related literature. In Section <ref>, we provide background on tensor networks. Section <ref> introduces the problem and its MILP formulation. Section <ref> discusses formulation reinforcements. Our numerical results are presented in Section <ref>.  Section <ref> concludes the paper.








§ RELATED WORK



The BMP is known to be NP-hard, even for trees <cit.>. In response, researchers utilize heuristics to address near-optimal solutions that balance computational efficiency with solution quality, as explored in studies<cit.>. Applications of genetic algorithms and ant-based systems in this context can be found in the work <cit.>, and a Tabu search method is employed in <cit.>. A particularly significant heuristic, which has seen extensive use in quantum emulation studies, is the RCM algorithm. First introduced by Cuthill and McKee <cit.>, the RCM algorithm starts with an arbitrary vertex and uses a breadth-first search to assign levels to vertices. It then sorts vertices based on these levels and further organizes them according to the number of neighbors in earlier levels. This process, repeated across all vertices, clusters them in a manner that minimizes the maximum distance between adjacent vertices in the graph, thereby reducing bandwidth, although it does not always guarantee the optimal solution. 

Another research direction in addressing the BMP focuses on developing theoretical findings, such as lower bounds, and exact solutions. This research direction is primarily motivated by the challenge posed by BMP's weak lower bounds and linear relaxations, which significantly complicates the task for exact computational methods. Initiated by Chvátal's density lower bound <cit.>, a variety of lower bounds have been proposed for estimating the minimum bandwidth. Among the notable contributions, two methods stand out for solving this problem. The first, introduced by Del Corso <cit.>, is particularly effective for small to medium-sized instances. The second, developed by Caprara and Salazar <cit.>, extends the first by implementing tighter lower bounds, thus making it capable of handling larger instances. Furthermore, a branch-and-bound algorithm for the BMP is proposed in <cit.>, and several formulations of the problem are studied in <cit.>. 

Beyond the quantum emulation application discussed in the introduction, the BMP has significant relevance in a wide range of other applications. Particularly in the context of sparse matrices, reducing bandwidth enhances cache utilization and minimizes memory access patterns. This leads to more efficient matrix operations, such as matrix multiplication. The BMP is also known role in solving linear systems of equations like Ax = b <cit.>, where pre-processing the matrix A to minimize its bandwidth substantially reduces computational efforts. In particular, Gaussian elimination on matrices with bandwidth b can be executed in O(nb^2), offering a considerable speed advantage over the standard O(n^3) algorithm when b ≪ n <cit.>. The BMP also finds practical applications in diverse fields such as electromagnetic industry, chemical kinetics, numerical geophysics <cit.>, information retrieval in hypertext systems <cit.>, and extends to graph embedding in the context of physical implementation <cit.>.





§ TENSOR NETWORKS FOR THE EMULATION OF QUANTUM MANY-BODY SYSTEM



A tensor is a multi-dimensional array that generalizes the concepts of scalars, vectors, and matrices. The rank of a tensor, indicating the number of indices or dimensions it possesses, determines its type: rank-0 for scalars, rank-1 for vectors, and rank-2 for matrices. Higher-rank tensors expand upon these basic structures. Tensors can be represented graphically through Penrose diagrams, here, tensors are represented as nodes, with the incident edges are indices, as seen in Figure <ref>. The rank is reflected by the number of edges connected to each node. Index contraction, a key operation in tensor calculus, involves summing over shared indices between tensors, analogous to the dot product in vector algebra. For instance, the matrix product C_αγ = ∑_β=1^D A_αβ B_βγ contracts the index β, merging dimensions of matrices A and B. Graphically, a contraction of two tensors is represented as an edge connecting nodes. 

More generally, Tensor Networks (TNs) represent complex tensor operations through graphical representations with multiple nodes corresponding to tensors, and multiple edges, each representing contracted indices. This visualization simplifies complex computations and clarifies multidimensional relationships. Figure <ref> illustrates TNs of various ranks. TNs include MPS, Projected Entangled Pair States
(PEPS), Tree Tensor Networks (TTN), and Multiscale Entanglement Renormalization Ansatz (MERA) which can be
defined in any dimension.





In quantum computing, a qubit is the fundamental unit of quantum information, analogous to a bit in classical computing. A qubit can exist in a superposition of two basis states, often denoted as |0⟩ and |1⟩. A quantum state of a qubit is a complex linear combination of these basis states, described by a vector in a two-dimensional Hilbert space. The composite state |Ψ⟩ of a many-body quantum system represents the collective state of all sites in the system. The wavefunction is a mathematical function that describes the quantum state, containing all possible states and their corresponding probabilities.

Emulating a quantum many-body system <cit.> often involves computing the systems final state after applying unitary operators <cit.>. This process models the systems evolution under different conditions. A significant challenge in emulating these systems is the exponential growth of the state space with the number of particles. Specifically, the wavefunction of an N-qubit system is represented by a tensor of rank 2^N. The computational complexity of handling such a tensor grows exponentially with N, making exact emulations intractable for large systems.


MPS offers a solution for the above challenge by transforming the high-rank tensor into a 1-dimensional chain TN, as seen in the Penrose diagram shown in  Figure <ref>. The bond dimension of an MPS, denoting the number of singular values retained in the tensor decomposition, is a critical factor in determining the accuracy and computational cost of the representation. Higher bond dimensions generally allow for more accurate representations but increase memory usage and computational time. The order of contractions in MPS, crucial for emulation's computational efficiency, refers to the sequence in which tensor contractions are performed. This sequence significantly impacts the computational resources required and can affect the bond dimension. In quantum computing, this often involves mapping qubits to the nodes of an MPS, dictating how tensors associated with different qubits interact.  
An optimal contraction order is key for efficient computation, particularly in emulating entangled states in quantum systems, and can help manage the bond dimension to balance accuracy and computational cost.











The interaction matrix U represents the strength and nature of interactions between different sites in a quantum many-body system. Inspired by a Rydberg atom array <cit.>, in this work the components of U are defined by u_ij = C_6/d_i,j^6, where C_6 is a device-dependent constant, and d_i,j denotes the physical distance between sites i, j. It has been demonstrated that minimizing the bandwidth of the interaction matrix – solving a BMP – can lead to a reduction in the bond dimension of MPS <cit.>. This is a significant insight because reducing the bond dimension can greatly decrease the computational complexity and resources needed for simulating quantum many-body systems. Typically, researchers have used heuristics like RCM to minimize the bandwidth of the interaction matrix. However, these heuristics do not, generally, take into account interaction strengths, and do not guarantee optimality. Thus, they might not always yield the most efficient reduction in bond dimension. 

Our study aims to advance beyond existing heuristics by exploring the exact resolution of the weighted-BMP problem. This variant of the BMP focuses on an ordering for the rows/columns of the interaction matrix that places the largest interactions closer to the main diagonal, thereby taking into account the interactions between sites. While this method may necessitate increased computational effort, it has the promise of improving the efficiency of tensor-to-MPS mapping through further reduction of the bond dimension. In the next section, we study the weighted-BMP problem. The reader is encouraged to consult <cit.> for comprehensive insights into TNs, and <cit.> for a foundational background in quantum computing. 



§ PROBLEM STATEMENT AND FORMULATION


In this section, we introduce the weighted-BMP problem, and present its MILP formulation.



 §.§ Problem statement



Let G = (V, E) be an undirected complete graph representing the quantum many-body system. Here, V denotes the set of qubits of size n, and E represents the interactions between them. Each edge uv ∈ E is associated with a distance d_uv.The interaction matrix U is then defined such that u_uv = 1/d_uv^6 for all uv ∈ E, intentionally excluding the constant C_6 for model simplification. Let π: V →{1, 2, …, n} be a bijective function that assigns each vertex u to a unique ordered position π(u). The bandwidth b of the system is expressed as b = max_(u, v) ∈ E u_uv |π(u) - π(v)|. Weighted-BMP is to find the bijection π that minimizes the bandwidth b. Such optimal ordering ensure that the vertices with the largest interactions are positioned as close as possible to each other. In terms of matrix representation, this ordering positionate the largest interactions as near to the matrix's diagonal as possible.



 §.§ MILP formulation

We define the following decision variables:


    b: a positive continuous variable representing the bandwidth,
        x_v^i = 
          1    if the position of  v  is  i, 
    
          0    otherwise,  ∀ v ∈ V,  ∀ i ∈{1, …, n}.

The weighted-BMP is equivalent to the following MILP


min      b        

s.t.      ∑_v ∈V x_v^i = 1            ∀i = 1, …, n,           

    ∑_i=1^n x_v^i = 1            ∀v ∈V,           

    ∑_i=1^n i(x_u^i - x_v^i) ≤b ·d_uv^6            ∀(u, v) ∈E: u ≠v,           

    b ≥0,        

    x_v^i ∈{0, 1}            ∀v ∈V, i = 1, …, n.


Constraints (<ref>) and (<ref>) define an ordering of V by ensuring that a unique vertex is assigned to each position, and exactly one position is assigned to each vertex, respectively. Constraint (<ref>) calculates the bandwidth for each distinct pair of vertices (u, v) ∈ V^2. Here, the term ∑_i=1^n i(x_u^i - x_v^i)  is designed to compute the difference in positions between vertices
u and v, indicating how far apart they are positioned within the ordering. The objective is to minimize the maximum bandwidth. This formulation is solved to optimality. 

While the MILP formulation of weighted-BMP is compact, involving n^2 + 1 variables and n(n+2) constraints, it encounters limitations due to its weak linear relaxation and the presence of inherent symmetries that affect the efficiency of the resolution algorithm within the solver. To address these challenges, the subsequent section will detail the introduction of symmetry-breaking inequalities, and the establishment of a theoretical lower bound, aimed at fortifying the formulation. These enhancements are designed to improve the algorithmic performance of the solver by mitigating the issues associated with the initial model weaknesses.




§ FORMULATION STRENGTHENING


From constraints (<ref>), it is easy to see that


    max_u,v ∈ V^2: u ≠ v1/d_uv^6≤ b.


Inequality (<ref>) gives a theoretical lower bound that we add to the MILP formulation, in order to strengthen its linear relaxation, thereby improving the convergence of the algorithm. In addition to that, we add symmetry-breaking inequalities (<ref>) to reduce the number of nodes in the resolution tree. These inequalities are based on the following result.



Let u be an arbitrary vertex in V, then there exists an optimal solution such that 

    ∑_i=1^n i x_u^i ≤⌈n/2⌉.




Denote by x(π) the incidence vector of any weighted-BMP solution π, where π is a bijective function defining an ordering: π: V →{1, 2, …, |V|}, that assigns to each u ∈ V a position π(u) in the ordering.

We prove the proposition by contradiction. Assume that for any optimal solution π, it holds that ∑_i=1^n ix_u^i(π) ≥⌈n/2⌉ + 1. Now, define the solution π' such that π'(u) = n + 1 - π(u) for all u ∈ V, which creates an inverse ordering of π.

Evaluating the bandwidths of solutions π and π', denoted b(π) and b(π') respectively, we observe that due to the symmetry of the inverse ordering, for every pair of vertices u, v ∈ V, the difference in positions |π'(u) - π'(v)| remains unchanged. Thus,
b(π') = max_u,v ∈ V^21/d_uv^6 |π'(u) - π'(v)| = b(π).

Therefore, π' also constitutes an optimal solution. However, as ∑_i=1^n ix_u^i(π) ≥⌈n/2⌉ + 1, π' verifies 


    ∑_i=1^n i x_u^i(π') = n + 1 - ∑_i=1^n x_u^i(π) ≤ n-⌈n/2⌉≤⌈n/2⌉,


contradicting our initial assumption and concluding the proof.






§ PERFORMANCE EVALUATION


In this section, we conduct numerical tests to evaluate the efficiency of the algorithmic reinforcements introduced previously, and assess the effect of weighted-BMP optimal solutions on optimizing the computational costs of quantum emulation, using a set of realistic instances.



  §.§.§ Instances description and implementation features


We generate multiple graphs that represent amorphous solids using <cit.>. Amorphous solids <cit.>, characterized by their lack of a crystalline structure, do not exhibit long-range order but maintain short-range order through specific bond lengths and angles akin to their crystalline counterparts. This feature leads to a well defined coordination number, C, denoting the number of nearest neighbors around each site (node degree). Our generated instances simulate amorphous solids with a coordination number C=4, spanning site counts of 10, 15, 20, 25, 30, 35, and 40. For each site count category, 10 distinct graphs were randomly generated to model the structure of amorphous solids. 

Our performance evaluation is based on three principal metrics: CPU time, the RCM gap to optimality, and the bond dimension. The RCM gap, expressed as a percentage, is defined by RCM   gap(%) = Obj(RCM) - Opt/Opt× 100,
where Opt is the best MILP objective value, and Obj(RCM) is the RCM objective value. The bond dimension, quantifies the MPS complexity, influencing the computational resources required for quantum emulation. 



The optimization process was executed in the Python programming language. The MILP formulation is solved to optimality using CPLEX 22.1 <cit.>, while fixing 10 hours as a time limit. On the other hand, the emulation process was executed in the Julia programming language using ITensors package <cit.>. The entire experiments were conducted on a system with an AMD EPYC 7282 16-Core architecture, in 64-bit modes, with a frequency of 2800 MHz.

To provide a preliminary visual understanding before delving into detailed numerical results, Figure <ref> displays how a weighted-BMP solution affects the interaction matrix in a 30-site instance. The figure highlights the repositioning of highest interactions around the matrix main diagonal, which results in a reduced bandwidth.







  §.§.§ Effect of formulation reinforcements


Figure <ref> illustrates the impact of formulation reinforcements on average CPU time for solving the MILP. A clear trend is observed: as the number of sites increases, so does the CPU time for both scenarios. The scenario with reinforcements consistently outperforms the one without, as evidenced by the lower average CPU times at all site numbers. The most significant difference occurs at 35 sites, where the time is nearly halved: 10782.55 seconds with reinforcements against 20442.85 seconds without. On average, incorporating reinforcements improves CPU time by 25.61%.  Yet, at 40 sites, the limitations of reinforcements are revealed, as they fail to confirm optimality within the time limit. This is essentially caused by the weakness of the linear relaxation. Future work will focus on introducing valid inequalities to address this.




The Figure <ref> presents the average RCM gap percentages, for various numbers of sites. The median values of these gaps, around 18% for 10 sites, 43% for 30 sites, and approximately 35% for 40 sites, indicates a consistent fail of the RCM heuristic from achieving optimality. Although some variability is observed in the data, the results identifies an overall average RCM gap to the best integer solution of 33.74%, which underscores the general trend of optimal solutions that clearly outperform the RCM heuristic. 





  §.§.§ Impact on quantum emulation



Figure <ref> provides a comparative analysis on the average bond dimension as a function of the number of sites, between RCM and Weighted-BMP in MPS simulations. For technical details on our MPS simulations, see Appendix <ref>. The orange curve, representing RCM, shows a consistent rise in the average bond dimension with more sites, which indicates the growing complexity and resource demands. Notably, the shaded area around the RCM curve expands with the number of sites, indicating greater variability and unpredictable performance in larger systems. In contrast, the Weighted-BMP method, represented by the green curve, not only achieves a lower average bond dimension across all site counts—with a notable difference at 35 sites where Weighted-BMP bond dimension is approximately 261.60 versus RCM's 417.10—but also displays reduced variability, as seen in the shaded area. This shows a more stable and predictable scaling performance. The overall improvement in average bond dimension offered by Weighted-BMP over RCM, quantified at 24.48%, underscores its clear advantage in managing the complexities of quantum emulations more efficiently.









§ CONCLUSION AND PERSPECTIVES


In this paper we have studied weighted-BMP, and its application for optimizing the tensors-to-MPS mapping in quantum emulation. We have modeled the problem using a MILP formulation, and have incorporated formulation reinforcements, such as symmetry-breaking inequalities and a theoretical lower bound. Experimental results have shown the benefits of introducing these reinforcements from an algorithmic standpoint. From an application perspective, we have shown the significant role of weighted-BMP in reducing quantum emulation costs for complex quantum systems. An interesting direction in the future is to generalize the weighted-BMP for mapping tensors to more complex structures like PEPS, TTN, and MERA. On the other hand, from an algorithmic viewpoint, this work has laid the foundation for future research, particularly in strengthening the MILP models linear relaxation, using valid inequalities. 





§ ACKNOWLEDGMENT

AD and JV acknowledge funding from the European Union under Grant Agreement 101080142 and the project EQUALITY. We also would like to acknowledge Kemal Bidzhiev for fruitful discussion. 


IEEEtran




§ APPENDIX



 §.§ Details of MPS simulations


We perform MPS simulations on the antiferromagnetic transverse field Ising model, applicable to Rydberg atom arrays as outlined in <cit.>. Specifically, we analyze N interacting spin S=1/2 sites within a two-dimensional amorphous solid, governed by the Hamiltonian:

    H=∑_i<j1/|𝐫_i-𝐫_j|^6σ^z_iσ^z_j+∑_i σ^x_i,

where the operators σ^α_i denote the Pauli matrices at the i-th site, positioned at 𝐫_i. We focus on simulating the non-equilibrium dynamics starting from an initial ferromagnetic state |000⋯0⟩, observing the system up to a time tJ=1 with time increments of dt = 0.01/J. We maintain a truncation error of 10e^-12 throughout the simulations to ensure energy conservation within an accuracy of 10e^-5.
