Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization Chuan-Chi Wang
Yan-Jie Wang
Chia-Heng Tu
Shih-Hao Hung
National Taiwan University Taipei, Taiwan [email protected]
National Taiwan University Taipei, Taiwan [email protected]
National Cheng Kung University Tainan, Taiwan [email protected]
National Taiwan University Taipei, Taiwan [email protected]
arXiv:2604.12256v1 [quant-ph] 14 Apr 2026
Abstract
states simultaneously. This capability allows for the processing of numerous possibilities all at once. A quantum circuit is a sequence of quantum gates, each operating on qubits, and represents a quantum algorithm designed for a specific computation. By mapping target problems to quantum circuit programs, quantum computers hold the potential to revolutionize various fields, such as cryptography, drug discovery, and combinatorial optimization problems. Quantum circuit simulation mimics the behaviors of quantum systems with classical computers. It is a crucial way to facilitate the development of quantum algorithms, given the prohibitive expense and noise issue (Noisy Intermediate Scale Quantum; NISQ [6]) that restricts the use of physical quantum computers. Quantum circuit simulation can be categorized by its internal representations for quantum states and simulation algorithms. In particular, full-state simulation and amplitude sampling methods are popular choices for building quantum circuit simulators. In this work, we prioritize full-state simulation for its accurate final probability distribution, suitability for debugging intermediate results, and efficiency in handling complex problems with deep quantum circuits, in contrast to amplitude sampling methods. The major challenge of full-state quantum circuit simulation is the exponential growth of memory space and computations. In particular, simulating an 𝑁 -qubit quantum circuit requires the memory space to keep 2N quantum states. Mathematical operations are required to compute the effects of the quantum gates in the circuit, and they need 2N state updates for each simulated gate. A widely adopted strategy to address this challenge involves the use of clustered computing environments [12, 19, 21, 22], which distribute both quantum states and computations across multiple nodes. Research efforts have also investigated leveraging storage capacity to extend the number of simulatable qubits [26]. To utilize the computational power of multicore CPUs and GPUs, several studies have incorporated multithreading techniques [10, 11, 14, 22, 31, 32, 34, 35]. Furthermore, various techniques for optimizing quantum circuits have been developed to transform input circuits into more efficient configurations, thereby enhancing overall simulation performance. Gate fusion, guided by cost functions, consolidates multiple gates into a single equivalent operation, significantly reducing computational overhead [15]. Qubit reordering involves remapping qubit indices to minimize memory movement, a particularly effective strategy in multi-device simulations [10, 12, 21]. To improve data locality in both single-node and multi-node environments, quantum gates of an input circuit can be partitioned into gate blocks, accompanied by enhanced data access schemes, providing an effective
Quantum circuit simulation is crucial for the development of quantum algorithms, particularly given the high cost and noise limitations of physical quantum hardware. While full-state quantum circuit simulation is commonly employed for prototyping and debugging, it poses challenges because of the exponential increase in simulation time for large quantum systems. In this work, we propose an extensible framework designed to enhance simulation performance by optimizing both data locality and computational efficiency, thereby addressing these challenges. This framework is seamlessly integrated with an optimizer that restructures quantum circuits and a simulator that adjusts execution strategies for various quantum operations. For the newly developed components, merge booster and diagonal detector, the underlying algorithms are inspired by the principles of quantum entanglement and gate fusion, as well as by the limitations identified in existing third-party simulation libraries. The experiments were conducted on eight DGX-H100 workstations, each equipped with eight NVIDIA H100 GPUs, employing both gate-level and circuit-level benchmarks. The results indicate a speedup of up to 160 times for circuit-level benchmarks and an acceleration of up to 34 times for diagonal-heavy gate-level benchmarks compared to existing simulators. The proposed methodologies are anticipated to deliver more robust and faster quantum circuit simulations, thereby fostering the advancement of novel quantum algorithms.
Keywords quantum computing, quantum circuit simulation, quantum circuit optimization, parallel programming, high-performance computation
1
Introduction
Quantum computing is an emerging technology exploiting quantum mechanics phenomena to address complex problems that pose challenges to classical computers. Different from classical computing, which represents a ‘0’ or a ‘1’ state with a bit, a quantum bit (also known as qubit) in quantum computing can represent both Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from [email protected]. Conference acronym ’XX, Woodstock, NY © 2018 Copyright held by the owner/author(s). Publication rights licensed to ACM. ACM ISBN 978-1-4503-XXXX-X/2018/06 https://doi.org/XXXXXXX.XXXXXXX 1
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
• A series of benchmarks is employed to evaluate the efficiency of the proposed methodologies. The simulator [37] is extended with our new methodologies and executed on a HPC cluster comprising a total of 64 NVIDIA H100 GPUs. The enhanced simulator consistently delivers superior performance, particularly in large-scale simulations.
solution [37, 38, 40, 41]. Detailed knowledge on these optimizations and limitations of prior approaches is described in Sections 2 and 3. To extend previous approaches with novel optimization techniques, we revisit quantum circuit partitioning and gate fusion methods, resulting in a systematic framework designed to enable efficient simulation. The proposed methodologies are anchored in two primary components: swarm optimization and adaptive simulation. The swarm optimization component integrates most of the leading optimization strategies for quantum circuits, while further incorporating innovative boosters applicable to a wide range of scenarios and an overhead-free gate fusion technique tailored for diagonal gates. The adaptive simulation component executes the optimized circuits generated by the swarm optimization, employing the most suitable simulation strategy at the most appropriate execution point. Consistent with the principle of extensibility, both components are designed to support additional customized optimization approaches. To validate our approach, the compact and lightweight quantum circuit simulator [37] is extended with the proposed methodologies and an MPI-based communication across the entire cluster. The experimental evaluation was conducted on a high-performance computing (HPC) cluster equipped with eight DGX-H100 workstations, demonstrating substantial performance improvements over leading simulators, including QuEST [22], AerSimulator [20], cuQuantum [10], and HyQuas [40]. The experiments encompass both gate-level and circuit-level benchmarks, with a detailed exploration of specific optimizations. Notably, our improved simulator achieves a 34× speedup over the mainstream simulators in executing the 38-qubit quantum Fourier transform program. The contributions of this work are summarized as follows. • A systematic and extensible simulation framework, incorporating circuit optimization and simulation, is proposed to expedite full-state quantum circuit simulations.
2
Background and Related work
This section introduces the necessary background information and related works on quantum circuit simulations, with particular emphasis on full-state simulation. The quantum bit and state representations are described in Section 2.1. Mathematical formulations of quantum gate operations are provided in Section 2.2. A qubit reordering technique for improving data locality is presented in Section 2.3. Subsequently, the data access schemes employed in quantum circuit simulation are detailed in Section 2.4. Ultimately, existing efforts are reviewed in Section 2.5.
2.1
Quantum Bit and State
Quantum bits (or qubits) are basic units of maintaining information for quantum computing. The state of a qubit can be characterized by a column vector [𝛼 𝛽]𝑇 = 𝛼 |0⟩ + 𝛽 |1⟩, where 𝛼 and 𝛽 are complexvalued numbers and equate to 1 in 𝐿2 norm. This is also referred to as a quantum state vector in this work. The basis states |0⟩ and |1⟩, corresponding to the column vectors [1 0]𝑇 and [0 1]𝑇 in Dirac notation, form the foundation for qubit state representation. In the full-state quantum circuit simulation, each state vector is represented as two 64-bit floating-point numbers, known as a complex amplitude. Extending this concept to an 𝑁 -qubit system, the quantum state can be formulated as following notation. This representation encapsulates the superposition of all possible basis states in the 𝑁 -qubit system.
• An adaptive simulator not only supports the most efficient implementation of quantum gates but also selects the most appropriate computation strategy based on the optimized circuit. In particular, a proprietary implementation of fused quantum gates is developed to overcome the practical limitations of existing libraries, as presented in Section 4.3.
|𝜙⟩ =
∑︁
𝑎𝑖 |𝑖⟩
(1)
𝑖 ∈ [0,2𝑁 )
2.2
Quantum Gate Operation
A gate operation in quantum computing performs a specific mathematical transformation on one or more qubits within a quantum circuit. A sequence of quantum gate operations can be represented as a series of matrix operations. Considering a single-qubit gate operation 𝑈 applied on the 𝑗th qubit of the state vector, the transformation can be denoted as 𝑈 𝑗 = 𝐼 ⊗𝑁 − 𝑗 −1 ⊗𝑈 ⊗ 𝐼 ⊗ 𝑗 , where 𝑈 𝑗 is a 2𝑁 ×2𝑁 matrix multiplication on the state vector. The matrix operation on a quantum state |𝜓 ⟩ can be formulated as follows, where 𝛼 and 𝛼 ′ denote the amplitude before and after the operation.
• A swarm optimizer incorporates the most effective techniques for quantum circuits. Given that this procedure can be categorized as an NP-hard problem [38], the complexity of all employed techniques is strictly constrained to polynomialtime to prevent transferring execution time from the simulation to the optimization stage, as described in Section 4.2. • A novel algorithm merge booster is developed to take advantage of the inter-gate block entanglement-free property. It improves computational efficiency while minimizing memory overhead, as will be detailed in Section 4.4. • A novel algorithm diagonal detector is developed to consider the commuting gates property while performing gate fusion. It significantly improves performance when simulating quantum circuits with a high proportion of diagonal gates, as described in Section 4.5.
" ′ # 𝛼𝑏 𝑏 ...0 𝑗 ...𝑏 0 𝛼 𝑁 −1 𝑁 −2 = 𝑈 𝑏 𝑁 −1𝑏 𝑁 −2 ...0 𝑗 ...𝑏 0 ′ 𝛼𝑏 𝑏 ...1 𝑗 ...𝑏 0 𝛼𝑏 𝑁 −1𝑏 𝑁 −2 ...1 𝑗 ...𝑏 0
(2)
𝑁 −1 𝑁 −2
A two-qubit 4 × 4 unitary gate 𝑉 performing the matrix multiplication on the 𝑗-th and 𝑘-th qubits can be written as follows, where 𝑗 is strictly greater than 𝑘. 2
Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization
𝛼𝑏′ 𝑏 ...0 ...0 ...𝑏 𝛼𝑏 𝑁 −1𝑏 𝑁 −2 ...0 𝑗 ...0𝑘 ...𝑏 0 𝑁 −1 𝑁 −2 𝑗 𝑘 0 𝛼 ′ 𝛼𝑏 𝑏 ...0 ...1 ...𝑏 𝑏 𝑁 −1𝑏 𝑁 −2 ...0 𝑗 ...1𝑘 ...𝑏 0 ′ = 𝑉 𝑁 −1 𝑁 −2 𝑗 𝑘 0 𝛼𝑏 𝑁 −1𝑏 𝑁 −2 ...1 𝑗 ...0𝑘 ...𝑏 0 𝛼𝑏 𝑁 −1𝑏 𝑁 −2 ...1 𝑗 ...0𝑘 ...𝑏 0 ′ 𝛼𝑏 𝑏 ...1 ...1 ...𝑏 𝛼 𝑁 −1 𝑁 −2 𝑗 𝑘 0 𝑏 𝑏 ...1 𝑗 ...1 ...𝑏 0 𝑁 −1 𝑁 −2
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Quantum Circuit File (e.g., OpenQASM)
(e.g., Qiskit, QuEST)
Simulation Results (a00 = 0.707 a01 = 0.000 ...)
Figure 1: Workflow of quantum circuit simulation.
𝑘
technical terminology and simulator configurations in this work are provided in Section 4.1. Gate-by-Gate Simulation Scheme. In this naive scheme, quantum gates specified in the input quantum circuit are evaluated one by one, and the entire state vector space is traversed and updated after each quantum gate evaluation. The pseudocode of the gateby-gate simulation scheme is given in Algorithm 1. This scheme is intuitive to implement and is widely adopted by numerous commercial simulators [20, 22]. However, it would incur a significant amount of cache misses when the size of the quantum states exceeds the size of the storage space (e.g., the last-level cache of the processor). In this case, accessing the quantum states dominates performance, and it would greatly improve performance if the states to be accessed could be cached to avoid accessing lower-level memory frequently.
Qubit Reordering and Qubit Permutation
In an N -qubit simulation, all 2N states are constructed by combining the individual states of each qubit through a tensor product. When the qubit order is reordered, it solely impacts the arrangement of qubits within a state vector. Equation 4 expresses the operation of a swap gate for the 𝑖-th qubit that exchanges the quantum states of the 𝑖-th and 𝑗-th qubits. Based on the results, the state of the 𝑖-th qubit in state vector 𝑎 is identical to the 𝑗-th qubit in statevector 𝑎 ′ . " # " # ′ 𝑎 ∗...∗1 𝑎 ∗...∗0𝑖 ∗...∗1 𝑗 ∗...∗ 𝑖 ∗...∗0 𝑗 ∗...∗ → ↦ ′ 𝑎 ∗...∗0 𝑎 ∗...∗1𝑖 ∗...∗0 𝑗 ∗...∗ 𝑖 ∗...∗1 𝑗 ∗...∗
(4) Algorithm 1 Gate-by-gate simulation scheme. 1: procedure gateByGateSim(stateVec, circuit)
This technique can directly influence the permutation within a quantum circuit. For instance, consider a 5-qubit system with the arrangement (𝑞 4𝑞 3𝑞 2𝑞 1𝑞 0 ). To swap the positions of 𝑞 3 and 𝑞 2 , it is essential to exchange the pairs of amplitudes corresponding to |∗03 12 ∗ ∗⟩ and |∗13 02 ∗ ∗⟩ in the mathematical representation. In quantum hardware, this is typically achieved by applying a SWAP gate between 𝑞 3 and 𝑞 2 to meet connectivity requirements and minimize the accumulation of gate errors [24, 30, 39]. The intentional integration of supplementary swaps within the circuit highlights the requirement to modify the indexing of subsequent gates. In state vector-based simulation, while this method does demand additional operations, the qubits for the following gates have already been rearranged to occupy the least significant bits for the classical computer. This arrangement signifies that the state represented by that qubit has been positioned closer to the computational unit, resulting in substantial improvements [10, 12].
2.4
Quantum Circuit Simulator
(3)
In a multithreaded simulation, each thread is typically designated to execute a matrix operation. As the number of target and control qubits for a given gate increases, the number of threads that can be utilized decreases, and each thread manages an exponentially larger workload. To alleviate this challenge for certain specialized gates, finer-grained optimization can be implemented [22].
2.3
Architecturedependent Circuit Optimizer
for gate in circuit do stateVec ← operate(gate, stateVec) 4: end for 5: end procedure 2: 3:
Block-by-Block Simulation Scheme. The core concept of blockby-block simulation can be brought forward in earlier research [18]. It is developed to simulate an input quantum circuit in a blockby-block style, thereby alleviating the memory access overhead while updating the quantum states. Its simulation methodology is specified in Algorithm 2, where the regions subject to optimization and evaluation are collectively referred to as gate blocks. These gate blocks are identified and optimized prior to the simulation. During the simulation, the quantum gates are evaluated in a block-by-block style to achieve data locality. Upon the evaluation of a gate block, the state vectors of these quantum gates can be accommodated in the higher-level memory of classical processors and processed consecutively, and are designated as chunks. The technique that transforms an input quantum circuit to facilitate the block-by-block scheme is referred to as cache blocking in this work. In contrast to the prior gate-by-gate simulation, the caching blocking approach requires an optimizer to accurately determine appropriate gate blocks for consecutive chunk operations in subsequent simulations. The caching blocking approach can be divided into two classes: the on-demand strategy [12, 21] and the strategy aimed at maximizing the size of the gate block [37, 38, 40, 41]. While the former class is straightforward and executes rapidly during the optimization step, the latter allows for more efficient overall simulation, which constitutes the primary focus of this work.
Simulation Scheme
Figure 1 illustrates a typical workflow adopted by modern quantum circuit simulation. The input consists of a structured file in a specific format (e.g., Quil [33] and OpenQASM [9]) to represent a raw quantum circuit. Subsequently, a quantum circuit optimizer, similar to a quantum compiler, is proficient in performing various quantum circuit optimizations, such as combining sequential quantum gates to reduce circuit depth or overall gate count. The optimizer tends to implement optimizations tailored for the characteristics of a quantum circuit simulator. The simulator mimics the behavior of quantum gates specified in the input quantum circuit on a classical computer and generates all amplitudes that comprise the final simulation result. The following paragraphs present three major schemes for full-state quantum circuit simulation. The details of the 3
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
2.5
Algorithm 2 Block-by-block simulation scheme. 1: procedure blockByBlockSim(stateVec, gateBlocks)
Related Work
Extensive research has been conducted on quantum circuit simulations. Various internal data representations, such as state vectors, tensor network states, and density matrices, have been employed by quantum circuit simulators, each exhibiting distinct advantages and drawbacks. Among these, the state vector method is particularly valued for its noise-free characteristics, attracting significant attention and optimization efforts from both quantum and classical computing researchers. However, the full potential performance of state vector-based simulation still needs to be realized, consolidated, and integrated. Hence, the fundamental acceleration techniques and scalability enhancements for both the quantum circuit and simulation, as presented below, have been systematically organized.
for block in gateBlocks do for chunk in stateVec do 4: for gate in block do 5: chunk ← operate(gate, chunk) 6: end for 7: end for 8: end for 9: end procedure 2: 3:
To realize this process, three feasible methods have been proposed to ensure that the required qubits, which are operated by the gates within a gate block, are kept in chunks. The first method involves the use of multiple fused-swap gates [12, 19]. The second method transforms the gate operations to bypass quantum matrix computations, thereby enabling simultaneous swapping for N pairs of qubits [37]. The third method leverages direct indexing to efficiently store data in the cache [40], which is recognized as the most efficient approach, especially for single-machine-node scenarios.
Parallel Simulation. Numerical efforts have been made for parallelizing state vector-based quantum circuit simulations that exploit data-level parallelism. Remarkable examples, such as qHiPSTER [32], qsim [35], and QuEST [22], take advantage of CPU multithreading support to boost simulation performance. There are works, such as Cirq [11] and qHiPSTER [32], leveraging the singleinstruction-multiple-data (SIMD) support on multicore processors for simulation acceleration. More recently, GPU-based approaches have been used to improve simulation efficiency, with key examples including Aer-Simulator [20], cuQuantum [10], HyQuas [40], and Queen [37]. The performance of cuQuantum and HyQuas is impressive in both cross-node and single-node simulations, respectively. Queen demonstrates both aspects but has yet to explore the potential benefits of third-party libraries.
Hybrid Simulation Scheme. A hybrid simulation scheme is proposed to further accelerate quantum circuit simulations by switching between the two aforementioned schemes, as detailed in Algorithm 3. This approach is particularly advantageous when specialized hardware is available for computational acceleration, such as Tensor Cores. In general, the scheme employs a block-by-block simulation strategy, referred to as SharedMem in the original work [40]. However, when a gate acts on a sufficiently large number of qubits, the evaluation of its operations becomes analogous to a general matrix multiplication (GeMM). In such cases, the hybrid scheme switches to a gate-by-gate simulation mode, termed BatchMV, and leverages high-performance matrix multiplication libraries, such as cuBLAS [1] and cuQuantum [10], to efficiently execute these operations. It is noteworthy that an additional matrix transpose, e.g., via cuTT [17], is required to reshape the data layout into a format compatible with cuBLAS1 [40].
Scalable Simulation. Typically, state vector-based quantum circuit simulations rely on primary memory to maintain full quantum states. To extend the capacity for simulating a larger number of quantum bits, one possible approach is to incorporate memories from different machine nodes. This allows a larger number of quantum state vectors to be distributed to these nodes. The Message Passing Interface (MPI) is commonly used to facilitate the necessary communication involving data transfers between nodes [16, 21, 22]. For efficient GPU-based simulation, NVIDIA Collective Communications Library (NCCL) [27] utilizes NVLink to achieve superior interconnect throughput across multiple GPUs [37, 38, 40]. An alternative, cost-effective solution has been suggested, which involves leveraging secondary storage devices, such as SSDs, to support larger-scale quantum circuit simulations. In comparison to NVLink methods, intra-machine disk I/O operations can also enhance the simulation efficiency [26].
Algorithm 3 Hybrid simulation scheme. 1: procedure hybridSim(stateVec, gateBlocks)
for block in gateBlocks do if isFusedGate then ⊲ Gate-by-gate Scheme (BatchMV) 4: transVec ← transpose(stateVec) 5: stateVec ← operateGeMM(fusedGate, transVec) 6: else ⊲ Block-by-block Scheme (ShareMem) 7: blockByBlockSim(stateVec, gateBlocks) 8: end if 9: end for 10: end procedure 2: 3:
Quantum Circuit Optimization. Additionally, quantum circuit optimizations can be employed to further accelerate simulations. These optimizations transform input circuits into more streamlined versions before the simulation begins, thereby reducing simulation time through techniques such as gate fusion and qubit reordering. Gate fusion is a well-established technique that combines consecutive quantum gates into a single generic quantum gate to minimize memory access overhead. Despite its advantages, fused general gates typically require a cost function to maintain parallel efficiency [15, 40]. Qubit reordering modifies the data access pattern of qubits used in quantum algorithms to improve data
1 According to HyQuas [40], the state can then remain in this configuration without
requiring reordering back to its original layout. 4
Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY State Vector
locality [10, 12, 21, 37]. The comprehensive implementation strategies for both intra-node and inter-node qubit exchanges have been thoroughly explored in Queen [37]. While potentially lowering readability for developers, encoding gate information into a bitmap can confer a substantial advantage in optimizer efficiency [40, 41]. These diverse optimization scopes have established a robust foundation, demonstrating the effectiveness of full-state simulation when combined with other supportive methodologies.
Sub-State Vector
Single-Qubit Gate
Controlled Gate
Qubit 0 1 2 3
3
Motivation
4
This section addresses the limitations of existing works, concentrating on three main aspects. First, the challenges related to portability and efficiency of the existing library in the hybrid scheme are presented. Second, the constraints of current gate block search algorithms in the block-by-block scheme are examined. Third, the shortcomings of existing gate fusion algorithms in the gate-by-gate scheme are addressed.
5 6 7
Figure 2: Memory allocation and quantum gate arrangements for the standard state vector-based simulation. The blue and yellow regions indicate the quantum operations applied to the full-state vector and the sub-state vector, respectively.
Limitation of the Adoption of Existing Library. State vectorbased simulation inherently provides a structured quantum representation, enabling the integration of diverse optimization techniques suited to various scenarios. A wide spectrum of advanced approaches has been proposed across numerous significant studies, making it impractical for any single framework to encompass them all. Conversely, the indiscriminate integration of every available technique can result in excessive complexity, ultimately hindering maintainability. Specifically, the use of third-party libraries [1, 17] to accelerate the computations for fused gates forces a reversion to the gate-by-gate scheme, as shown in Algorithm 3. These methods not only fail to preserve the cached property but often necessitate the eviction of cached data from shared memory to accommodate input loading. Furthermore, relying on such libraries significantly amplifies memory consumption and data transfer overhead, potentially leading to unmanageable out-of-memory issues and excessive data transfers during simulations. In contrast, implementing the All-in-one simulation [37] can help mitigate these challenges by reducing unnecessary resource overhead and enhancing the capacity for flexible optimization mechanisms.
are not entangled with the qubits outside of this partitioned circuit. This scenario is referred to as an inter-gate block entanglementfree situation. In the aforementioned example, since there is no entanglement between the qubits within the highlighted region and those external to it, the required number of state updates can be diminished to 19 × 24 . Capitalizing on this characteristic can result in a substantial reduction in memory access overhead; conversely, overlooking it may restrict potential performance enhancements. The implementation of this optimization entails additional costs, such as the need for extra memory buffers and supplementary tensor product operations. That is, when a parallel thread processes a partitioned quantum circuit, an auxiliary buffer becomes necessary for storing intermediate results, alongside a tensor product operation for reshaping the computed data to align with the layout of the original state vector. If this optimization is not meticulously managed, it may result in degraded computational efficiency and increased memory consumption. To address these trade-offs comprehensively, this work proposes the merge booster to identify suitable gate blocks within an input quantum circuit, striking a balance between memory utilization and computational overhead.
Limitation of Gate Block Search Algorithms. A widely recognized strategy for accelerating large-scale quantum circuit simulations involves utilizing a block-by-block scheme through gate block search algorithms, which effectively decompose a raw circuit into multiple gate blocks. These search algorithms are typically designed to maintain polynomial-time complexity [12, 37, 40], ensuring that simulations can be conducted within a practical timeframe for online applications. Regardless of how the circuit is divided into various blocks, full access to the state vector remains essential. Figure 2 illustrates an 8-qubit quantum circuit consisting of 40 gates, including 26 single-qubit gates and 14 controlled gates. Simulating this circuit requires a total of 40×28 state vector updates. However, although the highlighted (yellow) gate block operates on only 4 qubits, it still incurs 21 × 28 state vector updates. This access pattern represents a common limitation of contemporary full-state vector simulators [10, 11, 20, 22, 38, 40, 41]. One of our key findings is that we can further reduce the number of state vector updates when the quantum gates in a partitioned quantum circuit
Limitation of Gate Fusion Algorithms. In addition to modifying the access pattern, gate fusion is a commonly adopted optimization technique in state vector-based quantum circuit simulations for reducing computational overhead. It works by combining multiple consecutive quantum gates into a single equivalent matrix operation (e.g., arbitrary unitary gates), thereby minimizing the number of matrix multiplications required to simulate the original consecutive quantum gates. In particular, diagonal fusion is a specialized form of gate fusion, targeting consecutive diagonal gates, such as the Pauli-Z (Z) gate, the controlled-phase (CP) gate, and the two-qubit ZZ rotation (RZZ) gate. This optimization is particularly useful for quantum circuits containing a significant number of diagonal gates, such as the standard Quantum Approximate Optimization Algorithm (QAOA) [5, 23]. 5
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
4.2
Existing gate fusion algorithms overlook the fusion opportunities arising from the property of commuting gates and miss potential performance improvements enabled by these gates. Two quantum gates are said to commute if swapping their order of application does not affect the final quantum state. Leveraging this property, such gates can be efficiently detected and combined whenever feasible. As illustrated in Figure 4(a), a naive fusion algorithm stops the fusion process upon encountering the RX7 gate while scanning for potential fusion opportunities from top-left to bottom-right of a given quantum circuit. However, when the commuting gate property is taken into account, the RZZ8 and CP9 gates can be repositioned and fused with earlier gates before the RX7 gate, as shown in Figure 4(b). In this work, we devise a diagonal detector algorithm to take advantage of commuting gates to enlarge the quantum gates to be fused and accelerate simulation speed.
4
Methodology
In this section, we present the fundamental components of our proposed framework. Section 4.1 provides an overview of the essential parameter configurations required for quantum circuit simulations. In Section 4.2, we outline an integrated collection of efficient circuit optimizations. Furthermore, we discuss the design of the simulation module, which is adaptable to various hardware scenarios, in Section 4.3. Section 4.4 explains the merge booster, designed to streamline computations by eliminating unnecessary state updates. Lastly, Section 4.5 introduces the diagonal detector, employed for implementing gate fusion specifically on diagonal gates.
4.1
Swarm Optimization Module
Algorithm 4 illustrates the overall workflow of the swarm optimization module, which is strictly constrained to operate using polynomial-time algorithms and typically accounts for less than one-thousandth of the total execution time. This module is designed to enhance quantum circuit performance by systematically applying various optimization techniques across multiple hierarchical levels, including rank-level, machine-level, and computation-level operations. The gate block search algorithm (GBSA) builds upon the design of the fundamental simulator. This algorithm utilizes a unified function to separate circuits into several non-overlapping subcircuits, effectively tackling similar partitioning optimization problems. To simplify the overall workflow, the proposed GBSA is repeatedly invoked to jointly enable gate fusion and cache blocking techniques. Additionally, to improve performance across various dimensions, this algorithm can be replaced or extended as noted in previous works [12, 28, 38, 40]. To support the new optimization techniques, the merge booster produces entanglement-free blocks with efficiency, avoiding unnecessary quantum state updates. The diagonal detector consolidates consecutive diagonal gates into a single gate without specific ordering, simultaneously reducing both the number and depth of gates. Algorithm 4 Swarm optimization. Require: N , R, C, F, B, circuit 1: devGBs ← GBSA(N, N -R, circuit, 0) ⊲ Rank-level 2: for subCircuit in devGBs.size() do ⊲ Machine-level 3: GBs ← Booster(N, C, B, F, subCircuit) 4: GBs ← DiagonalDetector(N, D, subCircuit) 5: GBs ← GBSA(N, C, subCircuit, 0) 6: for GB in GBs do ⊲ Computation-level 7: fusedGBs ← GBSA(N, F, GB, 1) 8: end for 9: GBs ← fusedGBs 10: end for
Parameter Configurations
The quantum circuit simulation requires several configurations and system parameters. In the simulation framework, the data structure env encapsulates the environment, including critical parameters such as N, B, and C, as well as other relevant metrics. The parameter N signifies the total number of qubits within the simulation. The boosting optimization process is impacted by B, where 2𝐵 determines the minimum division size. The parameter C denotes the size of the sub-states extracted from the 2𝑁 amplitudes, which must be carefully selected to stay within the cache capacity for optimal performance. Additionally, the parameter R indicates the total number of ranks within the simulation system. For fusion strategies, the parameters D and F specify the maximum number of target qubits processed through diagonal and cost-based gate fusion, respectively. The simulation also relies on stateVec, which holds the quantum state composed of 2𝑁 amplitudes. The quantum operations applied during the simulation are represented by gate, while G signifies the total number of quantum gates within the circuit. The parameter targs specifies the target qubits associated with each gate. The quantum circuit itself, consisting of a sequence of gates, is labeled as circuit. To enhance efficiency, an optimized arrangement of operations, referred to as gateBlocks (GBs), is derived by considering the original circuit, the applied optimizations, and the system environment.
To summarize, the swarm-inspired optimization approach systematically applies circuit transformations, enhancing simulation efficiency and preserving flexibility for future extensions.
4.3
Adaptive Simulation Module
The adaptive simulation scheme judiciously selects the most suitable simulation strategy for efficient quantum circuit simulations, as demonstrated in Algorithm 5. This module takes as input the quantum state vector stateVec and an optimized sequence of gateBlocks. Each block within gateBlocks is processed based on one of three primary strategies: tensor product computation, diagonal fused gate simulation, or general gate simulation. When the tensor product condition is satisfied, the simulator computes the relevant sub-vectors, subVec1 and subVec2, and updates stateVec accordingly. In the case of fused diagonal gates, the simulator applies the fused diagonal operation directly to stateVec, thereby circumventing the overhead associated with generalpurpose matrix operations. For conditions requiring the rearrangement of the state vector, the task is delegated to the qubit reordering 6
Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY State Vector
function. In instances appropriate for an all-in-one simulation strategy, stateVec undergoes processing on a chunk-wise basis. This approach optimizes data locality and minimizes unnecessary computations, regardless of whether the gates are fused. For all other promising scenarios, this module provides opportunities for further optimizations through an extensible function region, thereby facilitating future enhancements.
TP
Tensor Product
Single-Qubit Gate
Controlled Gate
Qubit 0 1 TP 2 3 TP
Algorithm 5 Adaptive simulation scheme.
4
1: procedure adaptiveSim(stateVec, gateBlocks)
5
for block in gateBlocks do 3: if isTensorProduct then ⊲ Tensor Product 4: stateVec ← tensorProduct(subVec1, subVec2) 5: else if isFusedDiag then ⊲ Diagonal Simulation 6: stateVec ← operate(fusedDiagGate, stateVec) 7: else if isQubitReorder then ⊲ Qubit Reordering 8: stateVec ← reorder(qubits, stateVec) 9: else if isBlock then ⊲ All-in-one Simulation 10: for chunk in stateVec do 11: if isFusedGate then 12: chunk ← operate(fusedGate, chunk) 13: else 14: chunk ← operate(gate, chunk) 15: end if 16: end for 17: else ⊲ Applying further extensible optimization 18: doExtensibleOptimization(stateVec) 19: end if 20: end for 21: end procedure 2:
4.4
TP 6 7
Figure 3: Memory allocation and quantum gate arrangements for the proposed state vector-based simulation. Blue region represents the quantum operation applied to the state vector, and the grey block indicates the tensor product operation.
qubits increases, with cache blocking and gate fusion seamlessly integrated into the process. Overall, the merge booster achieves the time complexity of O (N𝑙𝑜𝑔(N) + G + C), where N represents the number of qubits, G denotes the total number of gates, and C is the overhead introduced when fusion and blocking techniques are enabled. The space complexity is O(N ), driven by a queue and a one-dimensional vector. Algorithm Design. The core idea can be considered as an extension of the classical merge sort algorithm, primarily divided into two components: division and merging, as depicted in Algorithm 7. The divide component in Algorithm 6 involves reducing the maximum memory usage by powers of two, which can ensure minimal additional memory consumption. For a 32-qubit simulation, only two 216 sub-state vectors are required to handle the most significant memory block. Besides, based on the microbenchmarking results for quantum gates, the reduction in qubits no longer yields exponential gains once the state vector can be cached. Therefore, the divide component provides a mechanism for controlling the minimum size of gate blocks.
Merge Booster
In Section 3, the idea of reducing the number of quantum state updates for the quantum gates without inter-gate block entanglement incurs additional memory and computation overheads. The merge booster algorithm is devised to balance the simulation efficiency and the incurred overhead. A larger gate block can improve simulation efficiency by reducing the number of extra computations required to merge the sub-state vector, but at the cost of allocating a larger memory buffer to store temporary quantum states. Conversely, a smaller gate block demands less memory but incurs a higher computational cost for merging sub-state vectors back into the original state vector. This algorithm aims to improve simulation time, memory consumption, and computational overhead.
Algorithm 6 Divide component of Merge Booster. 1: procedure Divider(N, divSize, que)
if N ≤ divSize then que.push(N ) 4: return 5: end if 6: Divider(N >> 1, divSize, que) 7: N ← (N & 1) ? N + 1 : N 8: Divider(N >> 1, divSize, que) 9: end procedure 2: 3:
An Illustrative Example. Given the raw quantum circuit in Figure 2, Figure 3 provides a clear illustration of the circuit optimized by the merge booster. With 40 gates in the 8-qubit circuit, a naive approach requires 40×28 accesses to update the 28 state vector for each gate operation. In contrast, the merge booster can reduce the access count to 1,336. The reduced access count is calculated based on quantum gates in 𝑛-qubit gate blocks: (6+6+7+7)×22 +(7+4)×24 +3×28 for state vector, with an additional 2 × 24 + 28 for tensor product overheads. Based on this empirical analysis, it can result in a 7-fold improvement in the number of state updates. Furthermore, the benefits become even more pronounced as the number of simulated
⊲ Halve the queue size
The merge component primarily relies on the collaboration between a queue and an array through their member functions for assigning specific indices to gate blocks. Moreover, it can be divided 7
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
Algorithm 7 Pseudocode of the Merge Booster.
Algorithm 8 Pseudocode of Diagonal Detector.
1: procedure MergeBooster(circuit, gateBlocks)
1: procedure diagonalDetector(env, gateBlocks)
que ← queue() 3: Divider(N, (N + B - 1) / B, que) 4: while que.size() > 1 do ⊲ Start merge component 5: accIdx ← 1 6: mergedVec ← [] 7: mergedSum ← 0 8: for i in que.size() do 9: gBlock ← [] ⊲ 1. Prepare data 10: boostSet ← set() 11: for j in range(accIdx, accIdx + que.front()) do 12: boostSet.insert(j) 13: end for ⊲ 2. Generate the gate blocks 14: gBlock ← genGateBlock(boostSet, circuit) 15: gBlock ← doFusion(gBlock) 16: gateBlocks ← genGateBlock(gBlock) ⊲ 3. Merge required data 17: mergedSum ← mergedSum + que.front() 18: if i & 1 then 19: mergedVec.push(mergedSum) 20: mergedSum ← 0 21: end if 22: accIdx ← accIdx + que.front() 23: que.pop() 24: end for ⊲ 4. Update merged states and gate block 25: for i in mergedVec.size() do 26: que.push(mergedVec[i]) 27: qubits ← getAllElements(mergedVec, i) 28: gateBlocks.insert(TensorProduct(qubits)) 29: end for 30: end while 31: end procedure
resList ← [] ⊲ A list of processed gates it ← gateBlocks.begin() 4: while it != gateBlocks.end() do 5: if !isDiagonal(it) then ⊲ Skip other types 6: resList.push_back(*it) 7: ++it 8: continue 9: end if 10: diagList ← [*it] ⊲ A list of diagonal gates 11: uList ← [] ⊲ A list of unprocessed gates 12: stopTable ← zeros(env.N ) 13: depSet ← set(it->targs) 14: diagSize ← 1 15: fIt ← next(it) 16: while fIt != gateBlocks.end() do 17: intSet ← intersect(fIt->targs, depSet) 18: if !isDiagonal(fIt) then ⊲ Not a diagonal gate 19: if !intSet.empty() then 20: stopTable ← update(stopTable, fIt) 21: uList.push_back(*fIt) 22: else 23: resList.push_back(*fIt) 24: end if 25: else ⊲ A diagonal gate is detected 26: depSet.insert(fIt->targs) 27: if checkStop(stopTable, fIt) then 28: break 29: end if 30: diagList.push_back(*fIt) 31: ++diagSize 32: end if 33: ++fIt 34: end while 35: if diagSize > 1 then ⊲ Do diagonal fusion 36: fusedGate ← doDiagonalFusion(diagList) 37: resList.push_back(fusedGate) 38: end if 39: resList.extend(uList) 40: it = fIt 41: end while 42: gateBlocks ← resList 43: end procedure
2:
2: 3:
into four distinct steps and is fully compatible with both gate fusion and cache blocking, unlike the hybrid strategy [40]. (1) Preparation Step: Extract the required qubits from the queue and place the corresponding data into the set (Lines 9-13). This step ensures that all necessary qubits are available for subsequent processing. (2) Generation Step: Generate the gate block within the circuit by identifying gates that involve qubits from the set. When the number of required qubits exceeds the available cache size or gate fusion optimization is feasible, the fine-grained gate blocks are generated to optimize both memory usage and computational efficiency (Lines 14-16). (3) Merge Step: Merge the data pairwise, compute the subsequent round of required qubits, and remove the corresponding elements from the queue (Lines 17-24). (4) Update Step: Update the queue with the newly required qubits, and insert the additional tensor product operations to the gateblock (Lines 25-29).
4.5
Diagonal Detector
Diagonal gates, such as Z, CP, and RZZ gates, possess unique properties that can be exploited to optimize quantum circuit simulations. Unlike general quantum gates, diagonal gates only alter the phase of quantum states without affecting probability amplitudes, which makes them ideal candidates for specialized fusion techniques. To efficiently identify and fuse the largest gate blocks composed of diagonal gates, we introduce the diagonal detector algorithm, which restructures quantum circuits using a linked-list-based approach that captures the dependencies of target and control qubits. 8
Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
The matrix operation of an N -qubit diagonal gate on the N qubit state vector is illustrated below. In a multithreaded execution environment, thread𝑖 is tasked only with performing complex multiplication between 𝜆𝑖 and 𝛼𝑖 of the following equation. 𝜆1 0 |𝜓 ′ ⟩ = 𝐷 N |𝜓 ⟩ = . .. 0
0 𝜆2 .. . 0
0 0 .. . 0
··· ··· .. . ···
stopping conditions are met, diagonal gates continue to be recorded in the diagList (Lines 30-31). This functionality enables algorithm can merge RZZ2 , RZZ5 , CP6 , RZZ8 , and CP9 into D44 successfully. In terms of complexity analysis, despite the presence of two while loops in the diagonal detector, the efficient organization of multiple lists ensures that each gate is accessed only once. As a result, both the time and space complexity are O(G), where G denotes the total number of gates in the circuit.
0 𝛼 1 0 𝛼 2 .. .. . . 𝜆2𝑁 𝛼 2𝑁
5
An Illustrative Example. The diagonal detector maximizes fusion opportunities without violating circuit constraints by analyzing gate dependencies and applying set intersections. As shown in Figure 4, our method significantly improves upon conventional approaches, allowing gates like RZZ2 , RZZ5 , CP6 , RZZ8 , and CP9 to be merged efficiently into a single operation D44 . On the contrary, a conventional approach merges RZZ5 and CP6 gates into D35 , and the combination of RZZ8 and CP9 gates leads to the synthesis of D37 . The conventional approach is straightforward to implement, but it hinders the potential performance improvements. Note that the conventional approach stops merging the gates when encountering a barrier, RX7 , during the visiting of the gates for fusion. Conversely, the diagonal detector further checks the property of commuting gates beyond the barrier to maximize the size of the gate block. In practice, excessive tensor products of diagonal matrices may still lead to increased optimizer runtime, especially when more qubits are involved. Therefore, a practical upper bound is introduced as a heuristic to mitigate this limitation.
Evaluation
The experimental setup for evaluating the proposed algorithms and our quantum circuit simulator is described in Section 5.1. Two types of benchmark sets are considered in our simulation environment: gate-level benchmarks and circuit-level benchmarks. The experimental results for the entire HPC system are presented in Section 5.2 and Section 5.3, respectively. Our performance results are compared with a range of state-of-the-art works, specifically those utilizing circuit optimizers that maintain moderate transpilation costs. To investigate single-GPU performance, Section 5.4 provides an analysis of boosting techniques for common circuits, while Section 5.5 reports ablation studies across various simulators.
5.1
Experimental Setup
The quantum circuit simulation is performed on an HPC cluster comprising eight NVIDIA H100 workstations [3], each equipped with eight NVIDIA H100 GPUs, resulting in a total GPU memory capacity exceeding 5 TB. Within each workstation, the GPUs are interconnected via NVSwitch, providing an intra-server bandwidth of 900 GB/s. The workstations are connected externally through InfiniBand, with an inter-server bandwidth of 200 GB/s. The system is powered by two Intel Xeon Platinum 8480 processors and runs on Red Hat Linux 8.5.0 with kernel version 4.18.0. Given this configuration, simulations of up to 38 qubits can be performed utilizing double-precision floating-point operations. In our experiments, simulation performance is measured in seconds, utilizing double-precision floating-point numbers for data representation across all 64 GPUs. To effectively demonstrate the impact of the various optimizations, the gate and circuit benchmarks are chosen to represent the complex scenarios. The gatelevel benchmark employs RZZ gates with full connectivity. The circuit-level benchmark focuses on programs with and without the potential for diagonal fusion optimizations, specifically the 5-level fully connected Quantum Approximate Optimization Algorithm (QAOA) [13] and the Quantum Fourier Transform (QFT) [7]. Due to the stable performance, the IBM Aer simulator with the cuQuantum [10] backend is used for single-node simulations, while the cusvaer simulator is employed for multi-node simulations. These simulators without gate fusion serve as the baseline and are referred to as Naive. HyQuas [40], which achieves the most efficiency under the hybrid scheme in GPU, is treated as the representative method and denoted as Hybrid 2 .
Algorithm Design. As shown in Algorithm 8, the diagonal detector restructures quantum circuits using a linked-list-based approach. Different types of gates are categorized into separate lists: resList stores the final optimized circuit, diagList collects diagonal gates for subsequent fusion, and uList tracks unprocessed gates. Additionally, a stopTable serves as an array and ensures compliance with state vector simulation rules. Initially, the algorithm processes gates in gateBlock sequentially. Gates that are not diagonal are immediately added to the resList (Lines 5–8) and bypass further processing. Diagonal gates, in contrast, are subjected to a set of conditional checks to identify and link mergeable gates for potential fusion. Accordingly, the detection focus shifts to the subsequent gates (Lines 16-34). If the subsequent gate is not a diagonal gate, instead of terminating the merge immediately, it is preferable to use the set intersection for further evaluation. In Lines 18-24, if the required qubits of the gate overlap with those of any previously encountered gates, this information is recorded in the stopTable, and the gates are then linked to the uList. This consideration can allow the algorithm to continue detecting subsequent gates. Specifically, as shown in Figure 4(b), the RZZ2 gate can disregard H3 and RY4 , continuing with further detection. The RX7 also adheres to the same bypassing rule. In the case of encountering a diagonal gate, the evaluation against the stopTable becomes a crucial step. Should the gate violate the established simulation rules, the process is terminated (Lines 27-29). Subsequently, the gates in the diagList that can be merged are processed (Lines 35-38). A case in point is the detection of RZZ11 following RZZ2 , which illustrates this condition. When no
2 The RZZ gates in the benchmarks are decomposed into equivalent CNOT and RZ gates
to account for their absence in HyQuas. Although this decomposition increases the total number of gates, it allows leveraging the advantages of Tensor Cores optimization, thereby enabling performance comparison under this configuration. 9
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Qubit
OPENQASM 2.0; include "qelib1.inc"; qreg q[4]; h q[0]; h q[1]; rzz q[0], q[1]; h q[2]; ry(0.2) q[3]; rzz(0.5) q[2], q[3]; cp(0.4) q[1], q[2]; rx(0.5) q[2]; rzz(0.2) q[0], q[3]; cp(0.2) q[0], q[1]; h q[3]; rzz(0.1) q[1], q[3];
0
H0
1
H1
2
H3
3
RY4
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
RZZ2
Qubit
RZZ8 CP9 RZZ11
CP6 RZZ5
RX7 H10
0
H0
1
H1
2
H3
RZZ2
D37 RZZ9
D35 RX6
3 RY4
H8
(a) Diagonal Gate Fusion via Naive Approach
Qubit 0
H0
1
H1
2
H3
3
RY4
RZZ2
Qubit
RZZ8 CP9 RZZ11
CP6 RZZ5
RX7 H10
OPENQASM 2.0; include "qelib1.inc"; qreg q[4]; h q[0]; h q[1]; h q[2]; ry(0.2) q[3]; D4(...) q[0], q[1], q[2], q[3]; rx(0.5) q[2]; h q[3]; rzz(0.1) q[1], q[3];
D44
0
H0
1
H1
2
H2
RX5
3
RY3
H6
OPENQASM 2.0; include "qelib1.inc"; qreg q[4]; h q[0]; h q[1]; rzz q[0], q[1]; h q[2]; ry(0.2) q[3]; D3(...) q[1], q[2], q[3]; rx(0.5) q[2]; D3(...) q[0], q[1], q[3]; h q[3]; rzz(0.1) q[1], q[3];
RZZ7
(b) Diagonal Gate Fusion via Diagonal Detector Optimization
Figure 4: Different ways to enable diagonal gate fusion: (a) a naive approach, (b) the proposed diagonal detector optimization. The left side represents the circuit before optimization, while the right side shows the optimized circuit. The D represents a diagonal gate, and the subscript in the gate notation indicates the execution order.
150 100 50 0
Naive Naivef Hybrid Oursall
1000 Elasped Time (s)
Elasped Time (s)
1200
Naive Naivef Hybrid Oursall
200
800 600 400 200
31
32
33
34
Qubit
35
36
37
0
38
(a) QFT.
31
32
33
34
Qubit
35
36
37
38
(b) 5-level fully-connected QAOA.
Figure 5: Elapsed time of the circuit-level benchmarks ranges from 31 to 38 qubits.
5.2
Gate-Level Benchmark
above. In contrast, our work demonstrates substantial performance gains, with Ours𝑎𝑙𝑙 completing the 38-qubit simulation in under one second and achieving a speedup of over 164× relative to the baseline. These results indicate that the proposed optimizations can significantly improve computational efficiency while maintaining scalability for large-scale quantum circuit simulations.
To explicitly indicate which optimizations are enabled, as illustrated in Table 1, the subscript f denotes the activation of gate fusion, b enables the merge booster, and all refers to the activation of all optimizations, including caching blocking, gate fusion, and merge booster. The performance speedup is indicated in parentheses following the simulation time. Overall, the Naive results exhibit exponential growth in runtime as the number of qubits increases. Applying gate fusion on the Naive simulator can yield up to 2× improvement. Hybrid provides notable speedups of up to 6.5× for larger qubit counts but may encounter out-of-memory (OOM) issues on the GPU for 37 qubits and
5.3
Circuit-Level Benchmark
Figure 5 presents the simulation times for QFT and QAOA, respectively, with qubits ranging from 31 to 38. Due to recurring out-ofmemory issues of the hybrid method, the 37-qubit and 38-qubit simulations are not specifically addressed and discussed here. 10
Large-Scale Quantum Circuit Simulation on HPC Cluster via Cache Blocking, Boosting, and Gate Fusion Optimization
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Table 1: Elapsed time of RZZ gates ranging from 30 to 38 qubits (unit: seconds), with speedup shown in parentheses.
Table 2: Comparison of per-gate performance and speedup on 31-qubit circuits (unit: milliseconds). Circuit BV [4] HS [36] QAOA [13] QFT [7] QV [8] SC [25] VC [29]
Qubit Naive Naive 𝑓 Hybrid Ours𝑏 Ours𝑎𝑙𝑙 30 2.68 2.29 (1.2x) 1.55 (1.7x) 0.02 (115x) 0.02 (157x) 31 3.59 3.16 (1.1x) 2.04 (1.8x) 0.03 (110x) 0.02 (191x) 32 3.98 3.25 (1.2x) 1.77 (2.2x) 0.06 (67x) 0.03 (130x) 33 6.08 4.26 (1.4x) 2.23 (2.7x) 0.10 (61x) 0.04 (157x) 34 10.06 6.26 (1.6x) 2.92 (3.4x) 0.20 (50x) 0.07 (140x) 35 18.69 10.37 (1.8x) 3.65 (5.1x) 0.39 (48x) 0.11 (164x) 36 37.83 18.53 (2.0x) 5.79 (6.5x) 0.80 (47x) 0.22 (169x) 37 75.45 37.37 (2.0x) OOM 1.69 (45x) 0.45 (166x) 38 151.52 74.24 (2.0x) OOM 3.50 (43x) 0.93 (164x)
Analysis of Boosting Optimization
Upon completing large-scale, multinode simulations, it is essential to understand the performance impact of various optimizations on individual hardware units. For a 31-qubit simulation on a single device, Table 2 reports the per-gate execution time of different simulation approaches across benchmark circuits drawn from various application domains [2]. The baseline Naive exhibits a consistent per-gate cost of approximately 33 to 43 milliseconds. Depending on the quantum circuit, enabling gate fusion typically yields a 3- to 5-fold performance improvement, but it also increases the variability of the average per-gate execution time. The cache blocking mechanism consistently reduces the per-gate cost to below 5 microseconds, resulting in a performance improvement of up to 10 times across all benchmarks. This indicates that modifying the simulation scheme should be prioritized over general gate fusion. With boosting optimization enabled (Ours𝑏 ), further improvements are observed, yielding speedups ranging from 1.4× to 3.5×. The Hidden Shift (HS) and Supremacy Circuit (SC) circuits achieve up to 2.9× and 3.5× speedups, respectively, as the entire circuit only requires full entanglement at later stages. These results indicate that integrating the booster also provides substantial performance benefits compared to gate fusion and cache blocking.
5.5
Naive 𝑓 12.27 7.19 12.60 11.74 4.33 6.93 8.31
Ours 4.76 3.78 4.26 3.44 3.39 3.18 3.47
Ours𝑏 3.51 1.31 2.35 1.26 2.21 0.91 1.28
Speedup𝑏 1.4x 2.9x 1.8x 2.7x 1.5x 3.5x 2.7x
different simulators. QuEST is slightly faster due to its specialized implementation for the typical gates. When gate fusion is utilized, the Aer Simulator, using cuTrust and cuQuantum as the backends, achieves a speedup of at least threefold in single-machine. A specialized simulator [23] achieves a speedup of more than 9.7 times by utilizing preprocessed and optimized code for diagonal gates. HyQuas employs its hybrid approach and leverages the cuBLAS library, resulting in a 5.3x optimization. It is worth noting that for programs exhibiting great potential for diagonal fusion optimization, such as QAOA, a simulator specifically optimized for diagonal fusion optimization is faster than a highly optimized one. Under typical conditions, our simulator outperforms all costbased fusion approaches. When the merge booster optimization is enabled, Ours𝑏 delivers an additional speedup. However, when only diagonal fusion is activated, Ours 𝑓 has a slightly lower performance than that of the specialized simulator. This discrepancy can be attributed to the restricted fusion size to fit the processor cache, which is necessary to maintain the efficiency of the tensor product during the optimization step for the circuit. Lastly, when all optimizations are enabled, our approach remains fully compatible and demonstrates a notable performance improvement of 29.6 times.
In the case of the QFT circuit, Naive 𝑓 achieves a 2.4× performance improvement over Naive when gate fusion is enabled. In the 36-qubit simulation, Hybrid attains speedups of 8.5× and 3.5× relative to Naive and Naive 𝑓 , respectively. Our framework can outperform the others, achieving a remarkable 34× speedup over Naive and operating 5.8× faster than Hybrid. For the QAOA simulation, Hybrid can consistently outperform Naive and Naive 𝑓 across all scenarios. With a robust implementation ensuring efficient crossnode communication, our simulator delivers superior performance across all tested scenarios, achieving speedups of up to 34 times.
5.4
Naive 43.21 42.95 41.99 33.94 41.25 41.30 41.39
Table 3: Elapsed time of 30-qubit fully-connected QAOA with different optimizations (unit: seconds). Simulator cuQuantum AerSimulator QuEST cuQuantum 𝑓 AerSimulator 𝑓 Lin et al 𝑓 [23] HyQuas Ours Ours𝑏 Ours 𝑓 Ours𝑎𝑙𝑙
Ablation Experiment
To systematically assess the impact of various optimizations on simulation performance in a fundamental setup, we conducted ablation experiments. The evaluation was conducted on the most complex 30-qubit QAOA circuit using a single NVIDIA H100 GPU, with a particular emphasis on our optimization techniques. In the simulation scheme without the caching block technique, there is generally minimal variation in execution times between the
6
Cache Fusion Boosting Time Speedup No No No 29.41 No No No 28.04 1.0x No No No 26.68 1.1x No Yes No 9.09 3.2x No Yes No 8.86 3.3x No Yes No 3.03 9.7x Hybrid No 5.57 5.3x Yes No No 5.02 5.8x Yes No Yes 2.76 10.6x Yes Yes No 3.07 9.6x Yes Yes Yes 0.99 29.6x
Conclusion
This work aims to advance a full-state quantum circuit simulation framework, inspired by our empirical findings of two key circuit properties—inter-gate block entanglement-free behavior and commuting gates—as well as a practical limitation in using existing third-party libraries to handle fused gate operations. 11
Conference acronym ’XX, June 03–05, 2018, Woodstock, NY
Chuan-Chi Wang, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung
To address these findings, the methodologies (merge booster and diagonal detector) are proposed and integrated into a quantum circuit optimizer, a common preprocessing step for accelerating quantum circuit simulation prior to execution. The developed simulator outperforms state-of-the-art systems by exploiting data locality and computational efficiency through the proposed merge booster and diagonal detector. On a single device, performance can exceed the baseline by more than ten times across all benchmarks. By extending the MPI implementation, our proposed approaches still achieve a speedup of more than one order of magnitude. Owing to the strong extensibility of our simulator, future work will explore incorporating the proposed algorithms into existing simulators to further broaden the impact of the methodologies.
[20] Ali Javadi-Abhari, Matthew Treinish, Kevin Krsulich, Christopher J. Wood, Jake Lishman, Julien Gacon, Simon Martiel, Paul D. Nation, Lev S. Bishop, Andrew W. Cross, Blake R. Johnson, and Jay M. Gambetta. 2024. Quantum computing with Qiskit. arXiv:2405.08810 [quant-ph] doi:10.48550/arXiv.2405.08810 [21] Chenyang Jiao, Weihua Zhang, and Li Shen. 2023. Communication Optimizations for State-vector Quantum Simulator on CPU+GPU Clusters. In Proceedings of the 52nd International Conference on Parallel Processing (, Salt Lake City, UT, USA,) (ICPP ’23). Association for Computing Machinery, New York, NY, USA, 203–212. https://doi.org/10.1145/3605573.3605631 [22] Tyson Jones, Anna Brown, Ian Bush, and Simon Benjamin. 2019. QuEST and High Performance Simulation of Quantum Computers. Scientific Reports 9 (07 2019). doi:10.1038/s41598-019-47174-9 [23] Yu-Cheng Lin, Chuan-Chi Wang, Chia-Heng Tu, and Shih-Hao Hung. 2024. Towards Optimizations of Quantum Circuit Simulation for Solving Max-Cut Problems with QAOA. In Proceedings of the 39th ACM/SIGAPP Symposium on Applied Computing (SAC ’24). ACM, 1487–1494. doi:10.1145/3605098.3635897 [24] Ji Liu, Peiyi Li, and Huiyang Zhou. 2022. Not All SWAPs Have the Same Cost: A Case for Optimization-Aware Qubit Routing. arXiv:2205.10596 [quant-ph] https://arxiv.org/abs/2205.10596 [25] Igor L. Markov, Aneeqa Fatima, Sergei V. Isakov, and Sergio Boixo. 2018. Quantum Supremacy Is Both Closer and Farther than It Appears. arXiv:1807.10749 [26] Hsu Nai-Wei, Chuan-Chi Wang, Chia-Hsin Hsu, Chia-Heng Tu, and Hung Shih-Hao. 2024. Toward cost-effective quantum circuit simulation with performance tuning techniques. Connection Science 36, 1 (2024), 2349541. arXiv:https://doi.org/10.1080/09540091.2024.2349541 doi:10.1080/09540091.2024. 2349541 [27] NVIDIA Corporation. 2024. NVIDIA NCCL. https://developer.nvidia.com/nccl [28] Daeyoung Park, Heehoon Kim, Jinpyo Kim, Taehyun Kim, and Jaejin Lee. 2022. SnuQS: scaling quantum circuit simulation using storage devices. In Proceedings of the 36th ACM International Conference on Supercomputing. 1–13. [29] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien. 2014. A variational eigenvalue solver on a photonic quantum processor. Nature Communications 5, 1 (July 2014). doi:10.1038/ncomms5213 [30] Vicente Pina-Canelles, Adrian Auer, and Inés de Vega. 2025. Improving and benchmarking NISQ qubit routers. arXiv:2502.03908 [quant-ph] https://arxiv. org/abs/2502.03908 [31] Qiskit contributors. 2023. Qiskit: An Open-source Framework for Quantum Computing. doi:10.5281/zenodo.2573505 [32] Mikhail Smelyanskiy, Nicolas P. D. Sawaya, and Alán Aspuru-Guzik. 2016. qHiPSTER: The Quantum High Performance Software Testing Environment. arXiv:1601.07195 [quant-ph] [33] Robert S. Smith, Michael J. Curtis, and William J. Zeng. 2017. A Practical Quantum Instruction Set Architecture. arXiv:1608.03355 [quant-ph] https://arxiv.org/abs/ 1608.03355 [34] Yasunari Suzuki, Yoshiaki Kawase, Yuya Masumura, Yuria Hiraga, Masahiro Nakadai, Jiabao Chen, Ken M. Nakanishi, Kosuke Mitarai, Ryosuke Imai, Shiro Tamiya, Takahiro Yamamoto, Tennin Yan, Toru Kawakubo, Yuya O. Nakagawa, Yohei Ibe, Youyuan Zhang, Hirotsugu Yamashita, Hikaru Yoshimura, Akihiro Hayashi, and Keisuke Fujii. 2021. Qulacs: a fast and versatile quantum circuit simulator for research purpose. Quantum 5 (Oct. 2021), 559. doi:10.22331/q-202110-06-559 [35] Quantum AI team and collaborators. 2020. qsim. doi:10.5281/zenodo.4023103 [36] Wim van Dam, Sean Hallgren, and Lawrence Ip. 2002. Quantum Algorithms for some Hidden Shift Problems. arXiv:quant-ph/0211140 [37] Chuan-Chi Wang, Yu-Cheng Lin, Yan-Jie Wang, Chia-Heng Tu, and Shih-Hao Hung. 2024. Queen: A quick, scalable, and comprehensive quantum circuit simulation for supercomputing. arXiv:2406.14084 [quant-ph] https://arxiv.org/ abs/2406.14084 [38] Mingkuan Xu, Shiyi Cao, Xupeng Miao, Umut A. Acar, and Zhihao Jia. 2024. Atlas: Hierarchical Partitioning for Quantum Circuit Simulation on GPUs (Extended Version). arXiv:2408.09055 [cs.DC] https://arxiv.org/abs/2408.09055 [39] Ge Yan, Wenjie Wu, Yuheng Chen, Kaisen Pan, Xudong Lu, Zixiang Zhou, Yuhan Wang, Ruocheng Wang, and Junchi Yan. 2025. Quantum Circuit Synthesis and Compilation Optimization: Overview and Prospects. arXiv:2407.00736 [quant-ph] https://arxiv.org/abs/2407.00736 [40] Chen Zhang, Zeyu Song, Haojie Wang, Kaiyuan Rong, and Jidong Zhai. 2021. HyQuas: hybrid partitioner based quantum circuit simulation system on GPU. In Proceedings of the 35th ACM International Conference on Supercomputing (Virtual Event, USA) (ICS ’21). Association for Computing Machinery, New York, NY, USA, 443–454. doi:10.1145/3447818.3460357 [41] Chen Zhang, Haojie Wang, Zixuan Ma, Lei Xie, Zeyu Song, and Jidong Zhai. 2022. UniQ: A Unified Programming Model for Efficient Quantum Circuit Simulation. In 2022 SC22: International Conference for High Performance Computing, Networking, Storage and Analysis (SC). IEEE Computer Society, 692–707.
References [1] 2024. cuBLAS: Basic Linear Algebra on NVIDIA GPUs. https://developer.nvidia. com/cublas [2] 2024. Qibojit Benchmarks: Benchmarking quantum simulation. https://github. com/qiboteam/qibojit-benchmarks [3] 2026. NVIDIA DGX H100. https://www.nvidia.com/zh-tw/data-center/dgxh100/ Quantum Com[4] Ethan Bernstein and Umesh Vazirani. 1997. plexity Theory. SIAM J. Comput. 26, 5 (1997), 1411–1473. arXiv:https://doi.org/10.1137/S0097539796300921 doi:10.1137/S0097539796300921 [5] Shin-Wei Chiu, Chuo-Min Yang, Shan-Jung Hou, Po-Hsuan Huang, Chuan-Chi Wang, Chia-Heng Tu, and Shih-Hao Hung. 2025. FOR-QAOA: Fully Optimized Resource-Efficient QAOA Circuit Simulation for Solving the Max-Cut Problems. In Practice and Experience in Advanced Research Computing 2025: The Power of Collaboration (PEARC ’25). Association for Computing Machinery, New York, NY, USA, Article 2, 11 pages. doi:10.1145/3708035.3736006 [6] Jerry Chow, Oliver Dial, and Jay Gambetta. 2021. IBM Quantum breaks the 100-qubit processor barrier. https://research.ibm.com/blog/127-qubit-quantumprocessor-eagle. [7] D. Coppersmith. 2002. An approximate Fourier transform useful in quantum factoring. arXiv:quant-ph/0201067 [quant-ph] [8] Andrew W. Cross, Lev S. Bishop, Sarah Sheldon, Paul D. Nation, and Jay M. Gambetta. 2019. Validating quantum computers using randomized model circuits. Physical Review A 100, 3 (Sept. 2019). doi:10.1103/physreva.100.032328 [9] Andrew W. Cross, Lev S. Bishop, John A. Smolin, and Jay M. Gambetta. 2017. Open Quantum Assembly Language. doi:10.48550/ARXIV.1707.03429 [10] The cuQuantum development team. 2023. cuQuantum. doi:10.5281/zenodo. 7806810 [11] Cirq development team. 2022. Cirq is a Python library for writing, manipulating, and optimizing quantum circuits and running them against quantum computers and simulators. https://github.com/quantumlib/Cirq [12] Jun Doi and Hiroshi Horii. 2020. Cache Blocking Technique to Large Scale Quantum Computing Simulation on Supercomputers. In 2020 IEEE International Conference on Quantum Computing and Engineering (QCE). IEEE. doi:10.1109/ qce49297.2020.00035 [13] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann. 2014. A Quantum Approximate Optimization Algorithm. arXiv:1411.4028 [quant-ph] [14] Vlad Gheorghiu. 2018. Quantum++: A modern C++ quantum computing library. PLOS ONE 13, 12 (dec 2018), e0208073. [15] Hiroshi Horii and Jun Doi. 2021. Optimization of Quantum Computing Simulation with Gate Fusion. https://ipsj.ixsq.nii.ac.jp/record/210570/files/IPSJ-QS21002023. pdf [16] Chia-Hsin Hsu, Chuan-Chi Wang, Nai-Wei Hsu, Chia-Heng Tu, and Shih-Hao Hung. 2023. Towards Scalable Quantum Circuit Simulation via RDMA. In Proceedings of the 2023 International Conference on Research in Adaptive and Convergent Systems (Gdansk, Poland) (RACS ’23). Association for Computing Machinery, New York, NY, USA, Article 3, 8 pages. doi:10.1145/3599957.3606215 [17] Antti-Pekka Hynninen and Dmitry I. Lyakh. 2017. cuTT: A High-Performance Tensor Transpose Library for CUDA Compatible GPUs. arXiv:1705.01598 [cs.MS] https://arxiv.org/abs/1705.01598 [18] Thomas Häner and Damian S. Steiger. 2017. 0.5 petabyte simulation of a 45qubit quantum circuit. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. ACM. doi:10.1145/ 3126908.3126947 [19] Satoshi Imamura, Masafumi Yamazaki, Takumi Honda, Akihiko Kasagi, Akihiro Tabuchi, Hiroshi Nakao, Naoto Fukumoto, and Kohta Nakashima. 2022. mpiQulacs: A Distributed Quantum Computer Simulator for A64FX-based Cluster Systems. arXiv:2203.16044 [cs.DC] 12