arXiv:2606.01852v1 [cs.DC] 1 Jun 2026
Parallelizing Large-Scale Tensor Network Contraction on Multiple GPUs Feng Pan∗
Hanfeng Gu∗
Paul Springer
Xipeng Li
Singapore University of Technology and Design
NVIDIA Corporation
NVIDIA Corporation
NVIDIA Corporation
Abstract—Exact tensor network contraction underpins quantum circuit simulation, quantum error correction, combinatorial optimization, and many-body dynamics. The dominant parallelization strategy, slicing, scales exponentially and incurs redundant computation. We present a multi-GPU framework that instead distributes intermediate tensors across devices with explicit communication, converting a fixed contraction path into a communication-efficient schedule via GEMM-oriented mode reordering and communication-aware mode distribution planning. Within a single DGX H100 node (8 GPUs, NVLink), distribution delivers 7–173× extra speedup beyond embarrassingly parallel slicing, capturing nearly all of the available compute reduction (87–101%) because NVLink’s high bandwidth keeps communication small relative to compute. Scaling the same four workloads to 1024 H100 GPUs over InfiniBand, the extra speedup beyond slicing ranges from 42× to 67,869×, demonstrating that communication-aware distributed contraction far surpasses slicing-based scaling limits for frontier tensor networks. Index Terms—tensor networks, multi-GPU, distributed contraction, performance modeling, cuTENSORMp
I. I NTRODUCTION Classical simulation of high-dimensional systems remains a central challenge in computational science. Tensor networks (TNs) address this challenge by factorizing a highorder tensor into a structured network of smaller tensors that can be stored and manipulated much more efficiently [1], [2]. Their importance is no longer limited to condensedmatter physics. They now form a broadly useful computational language for exact and approximate algorithms in quantum circuit simulation, quantum error correction, combinatorial optimization, quantum many-body dynamics, and tensorized machine learning [3]–[8]. Across these domains, the same core operation appears repeatedly: TN contraction, i.e., summing over shared indices to produce amplitudes, coset probabilities, counting results, partition functions, or other observables of interest. The need for scalable TN contraction is increasingly urgent. Random circuit sampling continues to test the classical boundary of quantum advantage, logical-qubit experiments demand stronger exact decoding baselines, exact counting problems on irregular graphs remain computationally difficult, and quantum hardware is beginning to probe many-body dynamical regimes that are hard to verify classically [9]–[12]. For large ∗ Both authors contributed equally to this research. Corresponding authors: F. Pan (feng [email protected]) and H. Gu ([email protected]).
instances, however, contraction rapidly becomes the dominant bottleneck. Its cost depends acutely on the contraction path: different orders of pairwise contractions can lead to orders-ofmagnitude differences in floating-point work, memory traffic, and peak intermediate size. Finding a globally optimal path is NP-hard in general [13], so practical software relies on highquality heuristics together with memory-management techniques such as slicing and path reconfiguration [3]. Even then, frontier-scale problems often remain limited not by arithmetic throughput but by memory capacity. Once an intermediate tensor no longer fits in the high-bandwidth memory of a single GPU, the contraction simply becomes infeasible. This paper focuses on that systems-level barrier. Our goal is not only to accelerate tensor contractions on multiple GPUs, but to make previously infeasible contractions executable by combining communication-aware distribution with tensoraware data layout. The central idea is that multi-GPU execution should be decided using the structure of the contraction path itself: which modes are kept, which are reduced, where large intermediate tensors appear, and where redistribution can avoid expensive communication later in the path. This paper makes three main contributions. First, we present a systematic framework for converting a fixed TN contraction path into a multi-GPU execution schedule, explicitly addressing data partitioning, communication, and execution ordering. Second, we introduce an offline scheduling strategy that combines GEMM-oriented mode reordering with communicationaware mode distribution planning; a dynamic-programming search guided by a hardware-aware cost model selects partitioned modes and redistribution points for a target GPU count. Third, we evaluate the framework experimentally on quantum circuit simulation and three many-body dynamics benchmarks, demonstrating generality across application families by scaling all four workloads from a single NVLink node to 1024 GPUs over InfiniBand. II. BACKGROUND A. Tensor Networks TNs represent a large tensor as a collection of lower-order tensors connected by shared indices [1], [2]. This factorized view is especially powerful when the target problem has locality, limited entanglement, or other structure that allows lowrank compression. In practice, the same formalism supports both approximate methods, such as matrix product state time
evolution, and exact contraction algorithms used in quantum circuit simulation and counting problems. A TN can be viewed as a graph whose nodes are tensors (multidimensional arrays) and whose edges denote shared indices. Closed indices connect two tensors and are summed over during contraction; open indices survive in the final output. We focus on exact contraction: summing over all closed indices P while preserving the open ones. For two tensors, Cij = k Aik Bkj ; larger networks apply the same rule repeatedly. B. Contraction Path and Complexity In most implementations, a TN is evaluated as a sequence of pairwise contractions. Although the final result is independent of the contraction order, the computational cost is not. This ordered sequence is called the contraction path and can be represented as a binary contraction tree whose leaves are the original tensors and whose internal nodes are intermediate tensors. The quality of a contraction path is usually measured by three related metrics: 1) Time complexity (Ct ): the total number of floatingpoint operations (FLOPs) required by the full contraction. 2) Space complexity (Cs ): the size of the largest intermediate tensor that appears along the path. 3) Memory complexity (Cm ): a proxy for total data movement, obtained by summing the sizes of the tensors read and written at each contraction step. Optimizing a contraction path for general TNs is NPhard [13]. Exhaustive search is therefore limited to very small instances, and practical software relies on heuristics, hypergraph partitioning, local tree rewrites, or stochastic search to obtain high-quality paths [3]. To define these costs, consider one pairwise contraction step a written as AM K × BKN → CM N , where M and N denote the retained modes and K denotes the reduced modes. Let Y Y Y m= dx , n= dx , k= dx . x∈M
x∈N
x∈K
The contraction-operation count of this step is then (a)
Ct
= mnk,
the corresponding peak space requirement is Cs(a) = max(mk, kn, mn), and the memory complexity is (a) Cm = mk + kn + mn.
For a full contraction path, the total time and memory complexities are additive over steps, X (a) Ct = Ct (1) a
Cm =
X
(a) Cm
(2)
a
whereas the space complexity is determined by the largest intermediate tensor encountered anywhere along the path: Cs = max Cs(a) . a
(3)
On modern accelerators, Cs often determines feasibility, while Ct and Cm determine performance once the contraction fits in memory. For clarity, the formulas above omit the constant factor associated with complex arithmetic. In the experimental FLOP and throughput numbers, tensors use complex64 arithmetic and one complex multiply-add is counted as eight real FP32 FLOPs (four multiplies and four additions), matching the operation counter used by the path logs. C. Contraction Techniques Even a good contraction path may still require an intermediate tensor larger than device memory. In that case, a standard remedy is slicing: one or more indices are fixed to specific values, the resulting smaller contractions are executed independently, and the partial results are summed at the end. Slicing reduces peak memory at the cost of repeating work across many subproblems. State-of-the-art optimizers such as cotengra [3] interleave slicing with path refinement (dynamic slicing), substantially reducing overhead. However, slicing does not remove the underlying source of hardness: contraction complexity grows exponentially with the treewidth of the network graph [13], and the largest intermediate tensor is also exponential in treewidth. For such large-treewidth problems, multi-GPU execution becomes more than a speed optimization. It is often the only way to overcome the single-device memory wall while preserving exactness. The challenge then shifts from path optimization alone to the coordinated design of data layout, communication, and redistribution across devices, which is the focus of this paper. D. Multi-GPU Hardware Platform We evaluate on the NVIDIA DGX H100 [14] (Table I): eight H100 GPUs per node, each with 80 GB HBM3, interconnected by NVLink (900 GB/s per GPU) within a node and InfiniBand (400 Gb/s per GPU) across nodes. This two-tier hierarchy— fast intra-node NVLink vs. slower inter-node InfiniBand— motivates our emphasis on limiting communication volume and avoiding fine-grained transfers. E. cuTENSORMp: Distributed Tensor Contraction Engine Translating a multi-GPU distribution plan into efficient device-level execution requires coordinating local tensor arithmetic with inter-device data movement. cuTENSORMp [15] is a multi-process extension of the cuTENSOR library that provides exactly this capability: given a distributed tensor contraction specified in Einstein summation notation, it automatically orchestrates local computation and inter-process communication across an arbitrary number of GPUs.
TABLE I K EY HARDWARE PARAMETERS OF THE DGX H100 PLATFORM . A LL BANDWIDTH FIGURES ARE BIDIRECTIONAL PER GPU.
DGX H100 GPUs per NVLink domain GPU memory (HBM) HBM bandwidth Peak FP32 FLOP/s per GPU Aggregate domain memory NVLink generation NVLink bandwidth / GPU Inter-node interconnect
8 80 GB 3.35 TB/s 67 TFLOP/s 640 GB 4th 900 GB/s 400 Gb/s IB
1) Programming Model: cuTENSORMp exposes a highlevel interface in which the caller specifies: 1) An einsum equation string defining the contraction, e.g., abcde,bdfg → acefg. 2) The global extents of each tensor mode. 3) A per-mode ranksPerMode array for each operand and the output, indicating how many ranks share each mode. A mode with ranksPerMode = 1 is replicated; a mode with ranksPerMode = dx is evenly partitioned across dx ranks along that dimension. From these descriptors the library constructs distributed tensor layouts, determines which data each rank owns, and builds an execution plan—all without requiring the caller to manage point-to-point transfers or collective schedules explicitly. 2) Pipelined Computation–Communication Overlap: For our purposes, the central systems feature of cuTENSORMp is its ability to overlap local computation with inter-device communication. Internally, the library decomposes a distributed contraction into a pipeline of three concurrent activities: • Local compute. Each GPU executes cuTENSOR tensor contraction kernels on its local shard of the operands. • Communication. NCCL [16] collectives and point-topoint transfers redistribute, broadcast, or reduce tensor shards as dictated by the distribution layout. • Tiling and double-buffering. Large contractions are decomposed into tiles; while one tile is being computed on the Streaming Multiprocessors, the results of the previous tile are simultaneously being transferred over NVLink or InfiniBand, and the operands for the next tile are being prefetched. This pipelining is important for the workloads in this paper, where a single distributed contraction step may move tens of gigabytes of intermediate data. 3) NCCL Transport and Network Awareness: cuTENSORMp delegates all inter-process communication to NCCL, which automatically selects the best available transport: NVLink and NVSwitch for intra-node transfers and InfiniBand RDMA for cross-node links. 4) Role in the Framework: The availability of cuTENSORMp creates a clean separation of concerns in our system. The offline planner (Section IV) decides which modes to distribute, when to redistribute, and when to gather—a combinato-
rial scheduling problem. cuTENSORMp then determines how each scheduled contraction is realized on the target platform, including tiling, workspace management, and communication orchestration. This separation lets the planner reason at the level of communication volume and transfer granularity, while the library maps those choices onto the available NVLink, NVSwitch, and InfiniBand substrate. F. Prior Work on Large-Scale TN Simulation The rapid growth of quantum processors has driven a parallel effort in large-scale classical simulation using TNs. We organize the most relevant prior work into three themes: contraction path optimization and slicing, supercomputer-scale quantum circuit simulation, and multi-GPU simulation systems. 1) Contraction Path Optimization and Slicing: At the algorithmic level, the performance of a TN simulator is largely determined by the contraction path. Markov and Shi [13] connected simulation complexity to treewidth, providing the graph-theoretic basis for path search and width reduction. More recently, path finders combine tree search, local rewrites, and dynamic slicing to find contraction schemes that trade peak memory for additional work. Huang et al. [17] showed that slicing can also expose substantial coarse-grained parallelism by decomposing a contraction into many independent subproblems, enabling more than 105 × acceleration on representative quantum-circuit benchmarks. However, these techniques optimize the contraction itself; they do not by themselves define how an individual contraction step should be partitioned and executed efficiently across many GPUs. 2) Supercomputer-Scale Quantum Circuit Simulation: Random circuit sampling (RCS) has been the highest-profile target for distributed TN simulation. Liu et al. [18] demonstrated real-time Sunway-based simulation of the Sycamore task using three-level parallelism across 42 million cores, reporting one million correlated samples in 304 seconds. Pan and Zhang [19] introduced the big-batch TN method for computing exact amplitudes or probabilities of large sets of correlated bitstrings, while Pan et al. [20] used 512 GPUs to generate one million uncorrelated samples for Sycamore circuits. Zhao et al. [21] later used 1432 NVIDIA A100 GPUs to generate three million uncorrelated samples with an XEB comparable to Sycamore in 86.4 seconds. Their implementation maps large stem tensors across 8 intra-node GPUs along the leading three dimensions and assigns sliced subtasks independently to nodes. Fu et al. [22] extended this line of work to 2304 GPUs at SC24, achieving a time-to-solution of 14.22 seconds through a system-level three-tier parallelization scheme with hybrid communication and low-precision quantization. Chen et al. [23] demonstrated the SW-TNC framework on Sunway, using data-reuse and step-fusion strategies to accelerate Zuchongzhi-60-24 simulation by over 10× on 1024 nodes. Viewed together, these RCS simulators rely primarily on slicing to expose coarse-grained parallelism and typically pair it with bespoke mappings of large intermediates onto fixed device hierarchies. Our work instead targets the execution
layer: we automatically choose which modes to partition, when to redistribute, and when to gather across an arbitrary number of GPUs, and then realize those decisions through cuTENSORMp mode-distributed contractions with pipelined computation–communication overlap. 3) Multi-GPU Simulation Systems and Libraries: On the systems side, Atlas [24] performs Schrödinger-style (statevector) quantum circuit simulation on multi-GPU clusters using hierarchical partitioning and ILP-based staging to reduce communication. In contrast to TN simulators, state-vector methods still scale with the full 2n state size, even when extended with host-memory or DRAM offload, whereas TN contraction can exploit circuit structure to reach larger qubit counts. Brown et al. [25] benchmarked multi-GPU quantum simulation across several GPU generations and showed that interconnect advances have contributed more than 16× improvements in time-to-solution. At the library level, cuTensorNet provides high-performance TN contraction, including distributed execution built mainly around slice-based parallelism [26], whereas our framework builds on the complementary cuTENSORMp library to support mode-distributed contractions with explicit computation–communication overlap. III. M OTIVATING A PPLICATIONS The need for distributed exact TN contraction extends well beyond quantum circuit simulation. Fig. 1 shows the computeonly complexity reduction (the FLOP reduction, ignoring communication) from distributing intermediates across GPUs for six workloads; in all cases, contraction cost drops rapidly as more devices contribute memory. Section V provides end-toend measured results confirming that these reductions translate into real wall-clock gains. A. Quantum Circuit Simulations Random circuit sampling (RCS) is one of the most demanding benchmarks for exact TN contraction because it intentionally generates highly entangled networks with large effective treewidth. As circuit depth increases, the computation becomes dominated by a small number of very large intermediates, making memory capacity and communication at least as important as arithmetic throughput. This benchmark has become a focal point in the race between quantum hardware and classical simulation algorithms [9], [27], [28]. On the classical side, TN methods have progressed from approximate contraction on arbitrary networks to big-batch amplitude evaluation, direct sampling, and large-scale GPU and supercomputer implementations [17]–[22], [29], [30]. This is the primary experimental focus of our paper, because Sycamore- and Zuchongzhi-style circuits sit squarely in the regime where a good path alone is insufficient and efficient multi-GPU execution becomes decisive. Fig. 1(a) shows the complexity reduction for Zuchongzhi n60m24 as the number of GPUs grows from 1 to 1024.
B. Quantum Error Correction Quantum Error Correction (QEC) is essential for realizing fault-tolerant quantum computing. Exact maximum-likelihood decoding can be written as a TN contraction over all error configurations consistent with an observed syndrome [6], [31], [32]. For surface codes under circuit-level noise, the resulting networks become effectively three-dimensional and quickly exceed single-GPU memory as code distance grows. This matters increasingly in practice because logical-qubit experiments and learned decoders both require strong exact baselines [10], [33]–[35]. Fig. 1(b) shows that exact distance7 rotated-surface-code decoding benefits strongly from distributed TN contraction, with contraction cost dropping rapidly as more GPUs participate. C. Combinatorial Optimization Problems Many counting and optimization problems can be reformulated as exact TN contractions by encoding local constraints into small tensors. These workloads are especially useful because they often generate irregular contraction trees and nonuniform communication patterns rather than the regular structures seen in lattice physics. Tropical and semiring formulations, including recent work on independent-set enumeration on King’s graphs, show that exact TN methods can push well beyond naive search [7], [11]. Fig. 1(c) shows that this irregular workload also benefits substantially from distributed contraction, reinforcing that the need for multi-GPU execution is not limited to circuit-like networks. D. Simulation of Quantum Many-Body Dynamics Real-time quantum dynamics can also be mapped to deep spacetime TNs through Suzuki-Trotter decomposition and related time-evolution methods [36], [37]. In two-dimensional lattices or long-time simulations, bond dimensions and intermediate tensors grow rapidly, so access to distributed highbandwidth memory becomes just as important as local compute throughput. This is especially timely because IBM, DWave, and Google experiments are already probing many-body or gauge-theory regimes that are increasingly difficult to verify classically [8], [12], [38]–[41]. Fig. 1(d–f) show the same rapid decrease in contraction cost across rectangular, hexagonal, and triangular lattices once the dominant intermediates can be distributed across more GPUs. All results in this section are estimated from fixed contraction paths produced by a state-of-the-art path finder. In other words, the challenge is no longer only to find a good contraction path in the abstract, but to parallelize that path across multiple GPUs with hardware-efficient mode layouts, distribution choices, and communication schedules. The next sections address this missing execution layer and show how our techniques convert such fixed paths into high-performance multi-GPU execution plans. IV. F ROM C ONTRACTION PATHS TO M ULTI -GPU E XECUTION We assume a fixed binary contraction path produced by an existing optimizer and focus on converting that path into
35 30 25 20 1
log10(Time Complexity)
(b) QEC (d7r7)
4
16
24 23 22 21 1
4
16
25
26
40
24
24
30
23
22
20
22
64 256 1K
(d) Rectangular 49×20
25
28
50
64 256 1K
Number of GPUs
1 40 35 30 25 20 15
4
16
64 256 1K
(e) Hexagonal 8×8
22 20 20 10 4
16
64 256 1K
Number of GPUs
20 15 8
32
64 128 512
(f) Triangular 49×24
40 30
30 25
4
24
1
(c) King's subgraph
Number of Slices
(a) QCS (n60m24)
25.5 25.0 24.5 24.0 23.5 23.0
Number of Slices
23
35
Number of Slices
log10(Time Complexity)
Time Complexity
30
22
25
21
20
20
15 1
4
16
64 256 1K
Number of GPUs
Fig. 1. Theoretical complexity reduction from distributing intermediate tensors across GPUs for six workloads: (a) quantum circuit simulation (Zuchongzhi n60m24), (b) exact QEC decoding (distance-7 rotated surface code), (c) independent-set enumeration on King’s subgraphs, and many-body dynamics on (d) rectangular, (e) hexagonal, and (f) triangular lattices. Red curves show log10 (total FLOPs); blue curves show the number of sliced indices. In all cases, contraction cost drops rapidly—often close to exponentially—as more GPUs contribute shared memory. Many-body panels show mean ± std over 15 independent path-finding runs.
an efficient multi-GPU execution. Our framework performs two offline transformations—GEMM-Oriented Mode Reordering and Communication-Aware Mode Distribution Planning— and then replays the resulting annotated schedule with a cuTENSORMp-based executor. Fig. 2 illustrates the end-toend workflow. The input is a contraction-path descriptor produced by an upstream path finder, containing the initial tensors, the pairwise contraction sequence, and any slicing metadata. The planner operates entirely offline: mode reordering first rewrites mode orders so that every pairwise contraction maps to a GEMM-like memory layout, then the distribution planner decides which tensor modes to partition across GPUs and when to redistribute or gather data. The executor replays the annotated schedule at runtime: one MPI rank per GPU loads its input, applies host-side slicing, and delegates each contraction step—including any inter-device communication— to cuTENSORMp [15]. A. GEMM-Oriented Mode Reordering Given a fixed contraction path of S pairwise einsum equations, mode reordering permutes the modes of each tensor so that every pairwise contraction admits a regular, matrixmultiplication-like layout—without changing the contraction tree or the mathematical result. Existing libraries such as cuTENSOR [42] and TBLIS [43] handle arbitrary mode orders
at runtime for individual contractions, inserting transpose or permutation kernels as needed. Mode reordering takes the complementary approach: it eliminates transposes ahead of time by choosing globally consistent mode orders across the entire contraction tree. Each pairwise contraction A × B → C partitions its modes into retained modes (surviving in C) and reduced modes (summed over). Grouping the retained modes as the GEMM m/n dimensions and the reduced modes as k, the canonical operand layout is [retained∥reduced], which maps directly to a matrix multiplication. The reordering pass seeks to realize this canonical structure at every step simultaneously. 1) Backward pass and emergent lifetime ordering: The key observation is that each mode in a contraction path has a well-defined remaining lifetime: the number of contraction steps until it is summed over. Reduced modes at step a have lifetime 0 (they are contracted at a); retained modes survive to some later step and therefore have lifetime ≥ 1; open modes that persist to the final output have the longest lifetime. A single backward pass over the contraction tree—from the last step to the first—propagates the [retained∥reduced] constraint from each consumer to its producers: 1) The output mode order of each step has already been set when the backward pass processed its downstream consumer. (The root’s output order is fixed by the problem specification, e.g., the open indices of the TN.)
trary local mode orders (gray) together with the downstreamimposed order I5 = gahe, which seeds the backward pass. Panel B shows the actual rewrite rule applied at each visited contraction:
Offline planner
Fixed Contraction Path einsums, tensor shapes, and slicing metadata
input = [shared in consumer order | reduced here].
1. GEMM-Oriented Reordering lifetime-ordered [retained | reduced]
2.DPDistribution Planning selects distributed modes and redistribution points
Annotated Schedule per-step layouts, distribution states, and gather / redistribute actions
Runtime
cuTENSORMp Executor host slicing
distributed einsums
overlap
Fig. 2. End-to-end workflow. The offline planner takes a fixed contraction path and applies two transformations: mode reordering produces GEMMready tensor layouts via lifetime ordering, and distribution planning selects distributed modes and redistribution points via DP search. The resulting annotated schedule is replayed at runtime by the cuTENSORMp executor with compute–communication overlap.
2) For each input operand of step a, extract the retained modes in the order they appear in the output (preserving the ordering already established) and append the reduced modes. This yields the operand layout [retained in output order∥reduced]. 3) The operand is the output of some earlier producer step p. Apply the same permutation to the output of p so that the data layout stays consistent. Each producer is modified at most once (the first consumer to request it wins). An important structural property emerges: after the backward pass, every tensor’s modes are sorted by remaining lifetime—longest-lived modes leftmost, shortest-lived rightmost. This is not an explicit sorting step; it falls out naturally from the recursive propagation of the [retained∥reduced] constraint. At each step, reduced modes (lifetime 0) are placed rightmost, and within the retained group, modes destined for later contraction are pre-sorted for their downstream consumers. The result is a globally consistent layout in which: 1) Every operand is in [retained∥reduced] form, mapping directly to a GEMM call. 2) Modes that survive many future steps sit in leading (leftmost) positions, while modes about to be contracted sit in trailing positions. 3) Each downstream consumer receives its inputs already in the layout it expects, eliminating runtime transposes. Fig. 3 illustrates this on a two-step subtree of a larger contraction tree. Dashed edges indicate that I1 –I3 are outputs of earlier contractions and I5 feeds into a later one; none are initial input tensors or the final result. Panel A shows arbi-
Starting from the fixed output order of I5 , the pass first rewrites the step-2 inputs as I4 = ae|bf and I3 = gh|bf, and then propagates the same consumer-imposed permutation backward to the producer outputs I1 = ab|cd and I2 = ef|cd. Panel C shows the resulting operand layouts, where retained modes precede reduced modes, and the dashed split line marks the GEMM boundary for each contraction input. Note the interleaved output of step 1: I4 = aebf rather than abef, because mode e (lifetime 2, from I2 ) outlives mode b (lifetime 1, from I1 ). Thus an output order may interleave modes from the two operands; the invariant is that each downstream consumer receives its inputs in [retained | reduced] form without any transpose. The backward pass is deterministic: the contraction tree uniquely determines every mode’s lifetime, so there is exactly one lifetime ordering and therefore one reordered scheme. All transformations act exclusively on the ordering of modes within each tensor; no mode labels are created, deleted, or reassigned, and the mathematical result is unchanged. B. Communication-Aware Mode Distribution Planning Given the reordered contraction path and a target GPU count P , the distribution planner decides (i) which contraction steps require multi-GPU execution, (ii) which tensor modes to partition across devices at each such step, and (iii) when to redistribute or gather data between steps. This problem is distinct from slicing (Fig. 4). Slicing eliminates a mode from the contraction entirely and creates independent sub-problems that share no data; slicing nslice indices yields 2nslice independent tasks that can be parallelized trivially but with redundant computation across tasks. In contrast, distribution partitions one or more modes of an intermediate tensor across devices so that each device holds a shard of the same logical tensor, reducing the total FLOP count at the cost of inter-device communication whose volume must be controlled. 1) Large Tensor Identification and Distribution: We first identify the large intermediate tensors that drive the need for distribution. A contraction step is considered large if any of its operands or its output exceeds a configurable memory threshold s (8 GiB in our experiments): L ARGE(a) ≡ max |A(a) |, |B (a) |, |C (a) | ≥ s. For each large tensor, we collect its use-chain: the sequence of subsequent contractions in which the tensor or its descendants participate. Any operand of a large contraction step may be distributed; cuTENSORMp supports distributing both operands and the output simultaneously.
A. Before: arbitrary mode orders I5
hgae
I4
g a h e downstream fixes output
sum c, d
cadb I1
f hbg
dfce I2
I3
B. Backward pass: rebuild each input from its consumer input = [ retained |
Rule
(consumer order)
backward (step 2 step 1)
Step 2: I4 × I3 I5
reduced ]
consumer I5 (fixed)
ghbf I3 consumer I4
efcd
I1
I2
C. After: every operand is [ retained | reduced ] I5
gahe sum b, f
I4 a e b f sum c, d
abcd I1
ef cd I2
Retained, long-lived Retained, short-lived
ghb f I3
Reduced at step Before reordering
Fig. 3. GEMM-oriented mode reordering on a two-step subtree. Dashed edges show connections to the rest of the contraction tree. Panel A: arbitrary local orders plus the downstream request I5 = gahe. Panel B: the backward pass visits step 2 and then step 1; at each visit, each input is rebuilt as “shared modes in consumer order, then modes reduced here.” Panel C: the resulting tree has [retained | reduced] operand layouts (dashed split line), with modes shaded by remaining lifetime and reduction status. The output I4 interleaves modes from I1 and I2 by lifetime, so step 2 can use I4 as [retained | reduced] without a transpose.
When a large tensor is first distributed across P devices, we select the minimum prefix of its leading modes whose product of dimensions is at least P : ( ) Y ′ D = arg min D ⊆ prefix(α) dx ≥ P . (4) ′ |D |
slice 1
slice 2
slice 3
contr. 4
contr. 4
contr. 4
contr. 4
contr. 3
contr. 3
contr. 3
contr. 3
contr. 2
contr. 2
contr. 2
contr. 2
contr. 1
contr. 1
contr. 1
contr. 1
contr. 0
contr. 0
contr. 0
contr. 0
GPU 0
GPU 1
GPU 2
GPU 3
slice i all GPUs cooperate contr. 4 (local) contr. 3 (local)
all local
one tensor, shared across GPUs
contr. 2 (redistribution required) contr. 1 (local) contr. 0 (local)
GPU 0
GPU 1
GPU 2
GPU 3
reduced FLOP count (fewer sliced indices)
Fig. 4. Slicing vs. distribution. (a) Slicing fixes indices to create independent sub-tasks with no communication but redundant computation. (b) Distribution partitions a single tensor across GPUs, reducing total FLOPs at the cost of communication that our planner minimizes.
reduce c, d
aebf
abcd
e outlives b (lifetime 2 vs 1), so the two operands interleave.
slice 0
reduce b, f
gahe
I4
I4 = ae bf, not abe f:
Full contraction (2m slices, m < n)
no communication redundant computation across slices
(operand-A order)
aebf Step 1: I1 × I2 I4
(b) Distribution (ours)
fix indices
sum b, f
be f a
(a) Slicing Full contraction (2n slices)
x∈D ′
The modes in D are called the distributed modes, and each device receives a contiguous shard of the tensor along these outermost dimensions. This prefix-based distribution interacts directly with the GEMM-oriented mode reordering of Section IV-A. Because all tensors are stored in row-major order, the leading (leftmost) modes correspond to the outermost memory dimensions. The backward pass sorts modes by remaining lifetime—longest-
lived leftmost—so the leading prefix modes selected for distribution are precisely the modes that survive the most contraction steps. This has two important consequences. First, because the distributed modes are outermost in row-major layout, each device’s shard is a single contiguous block in memory (no stride gaps), which maximizes the granularity of inter-device transfers and avoids the overhead of many small, non-contiguous messages. Second, the distributed modes are stable across consecutive contractions: because they are the longest-lived, they tend not to be contracted away in the near term, reducing the frequency of forced redistributions. 2) Distributed Mode Tracking: Once a tensor becomes distributed, its distributed modes must be tracked through subsequent contractions despite mode relabeling. Each mode label in the contraction path is tagged with the equation at which it will be contracted (or marked as an open index if it survives to the final output). At each step involving a distributed tensor, the planner identifies which of the current operand’s modes correspond to the previously distributed modes by matching these tags, and propagates the distributed status accordingly. Each step along a use-chain falls into one of four states: • Activate: the tensor is first distributed; prefix modes are selected per Eq. (4). • Keep: the output inherits the distributed modes from its operands; no communication is needed. • Redistribute: fresh distributed modes are selected for the output, triggering an all-to-all data shuffle. A redistribution is structurally forced when a currently distributed mode is about to be contracted at the next step; it may also be elective, chosen by the planner to maintain large contiguous blocks per device. • Gather: the tensor is small enough to fit on a single device; distributed modes are cleared and the data is replicated. 3) Redistribution Point Selection via Dynamic Programming: The key planning question is where along each usechain to place elective redistributions. For a use-chain of L equations, each presents a binary choice: keep or redistribute. The raw search space is 2L , but we collapse it with dynamic
Da−1
with a gather cost added at the final equation. Backtracing yields the optimal redistribution set R⋆ . In practice, the state space is small. When all distributed mode extents equal 2, the number of distributed modes is ⌈log2 P ⌉, and at each equation the DP tracks only the distinct
256 GB plateau (no elective redistribution) R
300
Keep (no communication) Redistribute (DP-chosen) Gather (exit)
valley
250 200 R
150
R
100 R
50 entry
R
R
exit (gather)
R
44 454 1
7 40
6 37
9
5
Equation index along use-chain
34
0
32
3
30
28
24
5
0 5
where P is the GPU count, Bnet the interconnect bandwidth, λ the per-message latency, nblk the number of contiguous blocks per device, and sblk the base block size. The block-granularity term naturally penalizes redistributions that produce many small messages (latency-bound) while treating large-block transfers as essentially free beyond their bandwidth cost. A natural consequence of this cost model is that the DP prefers to redistribute at points where the tensor is small: redistributing a 16 GB tensor costs ∼16× less than a 256 GB tensor at the same step. The DP thus automatically discovers that redistribution should concentrate at size valleys along the use-chain. b) DP formulation: Let Da denote the set of distributed modes at the output of equation a. The DP state at equation a is simply Da . At each equation, the planner propagates Da−1 to the current operand modes via mode-label tracking and evaluates two transitions: keep (inherit the propagated modes, (a) cost tkeep ) or redistribute (select fresh leading prefix modes, (a) cost tredist ). Because mode reordering has already sorted modes by remaining lifetime (Section IV-A), the leading prefix modes are the longest-lived, making them the most stable choice for distribution and minimizing the frequency of forced redistributions. When a distributed mode is nonetheless about to be contracted, redistribution is forced and only the redistribute transition is considered. The DP selects the minimumcost path through the use-chain: dp[a][Da ] = min dp[a−1][Da−1 ] + t(a) ,
21
block-granularity term
9
bandwidth term
16
where Drw is the total bytes read and written per device, Bdev is the device memory bandwidth, and Fdev is the peak FLOP rate. The redistribution time captures both bulk transfer and perblock overhead: sblk |C (a) | (P − 1) + nblk · max λ, , (7) tcomm = P · Bnet Bnet {z } | {z } |
2
The GEMM time is estimated from per-device tensor sizes after distributing the selected modes: Drw FLOPs , , (6) tgemm = max Bdev Fdev
13
(5)
98
(a) t(a) = t(a) gemm + ⊮[redist] · tcomm .
sets of ⌈log2 P ⌉ modes that can occupy the leading prefix positions. The total number of state evaluations across a typical use-chain is a few hundred, making the DP essentially instantaneous. c) Why deferring redistribution can hurt: A subtle failure mode motivates the DP over simpler heuristics. If redistribution is deferred too long, the distributed modes may drift into unfavorable stride positions as subsequent contractions reshape the tensor. When a forced redistribution finally becomes unavoidable—because a distributed mode is about to be contracted—it may occur at a point where the tensor is large (hundreds of gigabytes) and the stride pattern produces millions of tiny memory blocks, making the transfer latencybound and extremely expensive. The DP avoids this by redistributing proactively at points where tensors are small and block granularity is favorable, rather than waiting until a forced redistribution occurs at a costly peak. Fig. 5 illustrates the DP result on the Zuchongzhi n60m24 benchmark with 8 GPUs. The distributed tensor’s output size varies from 4 GB to 256 GB across a 50-equation use-chain. During the initial 256 GB plateau (equations 132–205), no elective redistribution occurs—only one forced redistribution at equation 169, where the current distributed modes are about to be contracted. The DP then places proactive redistributions in the valley (equations 223–300, where tensors are 16– 32 GB) and at local minima within the second peak (equations 349 and 376, both at 128 GB). The total redistributed volume is 596 GB—only 4.6% of the overall data movement— demonstrating that the DP concentrates communication at the cheapest points along the use-chain.
Output tensor size (GB)
programming by observing that the future cost depends only on which modes are currently distributed, not on the full history of past decisions. a) Cost model: Each transition from equation a to a+1 incurs a cost composed of a local-compute term and, if redistributing, a communication term:
Fig. 5. DP redistribution-point selection for the Zuchongzhi n60m24 benchmark on 8 H100 GPUs. The shaded area shows the distributed tensor’s output size along the use-chain. Blue dots are equations where distributed modes are kept (no communication); red diamonds mark DP-chosen redistributions. The DP avoids redistributing on the 256 GB plateau, instead concentrating transfers in the valley and at local minima of the second peak.
V. E XPERIMENTS All experiments use DGX H100 nodes (Section II-D). All tensors use complex64 (single-precision complex) storage with FP32 arithmetic for contraction kernels. For each GPU count, the path finder is run with a fixed time budget (approximately one hour) to produce a contraction path optimized for that configuration’s memory capacity. Because path optimization is NP-hard, the quality of the best path found within this
budget varies across GPU counts. Once a path is selected, it is held fixed and our framework applies mode reordering (Section IV-A) and communication-aware distribution planning (Section IV-B) on top of it. Thus, the scaling results below measure the end-to-end, capacity-aware workflow: increasing the GPU count changes both the executable distributed schedule and the amount of slicing required by the selected path. They should not be interpreted as fixed-path strong scaling; the non-monotonic trends below expose the remaining pathquality variation. A. Metrics For a configuration using P GPUs, let bP denote the number of sliced bonds in the selected path and let tP denote the measured wall-clock runtime for one sliced subproblem executed by the distributed contraction engine. When the full contraction contains 2bP independent slices, we report the projected full-contraction time TPproj = tP 2bP .
(8)
The projected full-contraction speedup relative to the one-GPU configuration is SP =
T1proj TPproj
=
t1 2b1 . tP 2bP
(9)
We do not execute every slice when 2bP is very large; instead, we measure per-slice execution time and combine it with the path’s sliced-bond count using the standard slicing cost model above. This projection does not include external slice scheduling overhead or the final accumulation across slices. The baseline is ideal embarrassingly parallel slicing, in which the contraction is decomposed into independent subtasks distributed evenly across GPUs with no inter-device communication. In this regime, the speedup equals the GPU count: P GPUs yield a P × speedup. The extra speedup is SP , (10) P quantifying the benefit of distributed contraction beyond ideal slicing. The complexity reduction is EP =
RP =
Ct,1 , Ct,P
(11)
where Ct,P is the full-contraction FLOP count, including all slices, for the path selected at P GPUs. B. Scaling Across Applications We evaluate four workloads: Zuchongzhi n60m24 quantum circuit simulation and three many-body dynamics benchmarks on hexagonal, rectangular, and triangular lattices. The 1–8 GPU points run within a single DGX H100 NVLink domain; larger GPU counts span multiple DGX H100 nodes connected by InfiniBand. Table II reports the 8-GPU single-node point, and Fig. 6 shows the full sweep from 1 to 1024 GPUs. At the 8-GPU NVLink point, all four workloads exhibit significant super-linear speedup beyond the 8× slicing baseline, with extra speedups ranging from 7.0× (triangular)
TABLE II S INGLE - NODE RESULTS ON 8 H100 GPU S (NVL INK ). S PEEDUP IS PROJECTED TO THE FULL CONTRACTION USING E Q . (9). E XTRA SPEEDUP IS PROJECTED SPEEDUP DIVIDED BY THE 8× EMBARRASSINGLY PARALLEL BASELINE . C OMPLEXITY REDUCTION IS THE COMPUTE - ONLY FLOP REDUCTION FROM DISTRIBUTING INTERMEDIATES ( COMMUNICATION - FREE , HENCE AN OPTIMISTIC ESTIMATE RATHER THAN AN ACHIEVABLE BOUND ).
Workload
Full Extra Complexity TFLOP/s speedup speedup reduction per GPU
Circuit (n60m24) 148× 18.5× Hexagonal 8×8 1,383× 172.9× Rectangular 49×20 75× 9.4× Triangular 49×24 56× 7.0×
18.5× 197.8× 9.3× 7.4×
28.1 31.6 32.9 31.8
to 172.9× (hexagonal). Because NVLink’s high bandwidth (900 GB/s per GPU) keeps the communication term small relative to compute, distribution converts essentially all of the available compute reduction into wall-clock gains: the achieved extra speedup reaches 87–101% of the compute-only FLOP reduction, and all workloads sustain 28–33 TFLOP/s per GPU. To our knowledge, these are the first distributed multi-GPU TN contraction benchmarks for real-time quantum many-body dynamics on two-dimensional lattices. C. Multi-Node Scaling to 1024 GPUs We next scale all four workloads to 1024 GPUs across multiple DGX H100 nodes connected by InfiniBand, where inter-node bandwidth is substantially lower than intra-node NVLink. Table III summarizes the 1024-GPU results. TABLE III M ULTI - NODE RESULTS ON 1024 H100 GPU S OVER I NFINI BAND . RUNTIME IS MEASURED PER SLICED SUBPROBLEM ; PROJECTED FULL SPEEDUP COMBINES SLICING PARALLELISM AND DISTRIBUTED CONTRACTION . E XTRA SPEEDUP ISOLATES THE GAIN BEYOND THE 1024× EMBARRASSINGLY PARALLEL BASELINE AND IS OUR PRIMARY METRIC ; COMPLEXITY REDUCTION IS THE COMPUTE - ONLY FLOP REDUCTION ( COMMUNICATION - FREE ), NOT AN ACHIEVABLE BOUND .
Workload Circuit (n60m24) Hexagonal 8×8 Rectangular 49×20 Triangular 49×24
Per-slice Sliced Projected full Extra Complexity runtime (s) bonds speedup speedup reduction 20.19 113.27 34.70 12.19
20 6 14 14
42,759× 41.8× 69.5M× 67,869× 221,212× 216.0× 135,728× 132.6×
418× 1.49M× 3,154× 986×
At 1024 GPUs, all four workloads remain well above the linear slicing baseline. We take the extra speedup—the gain beyond the 1024× embarrassingly parallel baseline—as the primary, technique-attributable metric: it ranges from 41.8× on the quantum-circuit benchmark to 67,869× on the hexagonal many-body benchmark, the latter driven by a reduction in sliced bonds from 37 to 6. The corresponding projected fullcontraction speedups in Table III are far larger because they also fold in the path’s slicing parallelism; since that factor is governed by the path finder and the available memory rather than by our execution layer, we report it for completeness but treat the extra speedup as the headline result.
End-to-end speedup
Quantum circuit (Zuchongzhi n60m24)
Many-body dynamics (hexagonal 8x8)
42.8K×
104
105
102
103
101
Slicing Ours
100 1
8
64
GPUs
1K
69.5M×
107
103
Many-body dynamics (rectangular 49x20)
101 1
8
64
GPUs
1K
105
Many-body dynamics (triangular 49x24)
221.2K×
105
104
104
103
103
102
102
101
101
100
100 1
8
64
GPUs
1K
135.7K×
1
8
64
GPUs
1K
Fig. 6. Projected full-contraction speedup vs. embarrassingly parallel slicing (dashed) on 1–1024 H100 GPUs, computed from measured per-slice runtime and sliced-bond count. The dotted vertical line marks one DGX H100 node (8 GPUs); larger configurations use InfiniBand across nodes.
1) Non-monotonic scaling: Neither the extra speedup nor the complexity reduction increases monotonically with GPU count. For example, on Zuchongzhi n60m24, the extra speedup reaches 18.49× at 8 GPUs, drops to 5.44× at 16 GPUs, and later recovers to 41.76× at 1024 GPUs; similarly, the complexity reduction peaks at 522.54× at 512 GPUs but falls to 418.03× at 1024 GPUs. Two factors contribute to this nonmonotonicity. First, contraction-path quality varies across GPU counts: because path optimization is NP-hard, we allocate the same fixed time budget (approximately one hour) to the path finder for each configuration, and the best path found within that budget can differ substantially—this directly affects the complexity reduction, which depends solely on the path. Second, 8 GPUs reside within a single NVLink domain on DGX H100, where redistribution is fast and cheap; at 16 GPUs, the execution spans two nodes connected by InfiniBand, introducing higher-latency cross-node communication that further reduces the achieved extra speedup. Notably, even when scaling beyond the NVLink domain onto InfiniBand interconnect, the framework continues to deliver substantial extra speedups. At 1024 GPUs, the extra speedup ranges from 41.8× on the quantum-circuit workload to 67,869× on the hexagonal many-body workload, demonstrating that communication-aware distributed contraction remains highly effective even over lower-bandwidth inter-node networks. 2) Communication, not arithmetic, bounds the achieved speedup: The complexity reduction column in Table III reports the compute-only FLOP reduction from distributing intermediate tensors. It deliberately excludes communication and is therefore not an achievable performance bound: distribution is fundamentally a compute–communication trade, and on InfiniBand the communication term is substantial. The meaningful comparison is thus between interconnect regimes, not against this communication-free estimate. Within a single NVLink node, where inter-GPU bandwidth is high, the extra speedup tracks the compute-only reduction almost exactly (87–101%).
Across nodes over InfiniBand, redistribution, cross-node collectives, and broadcast costs raise the communication term, yet the extra speedup remains large in absolute terms (41.8– 67,869×). In short, the benefit of distribution is limited by communication rather than by arithmetic, and it grows with interconnect bandwidth. VI. C ONCLUSION We presented a multi-GPU framework for exact TN contraction that moves beyond embarrassingly parallel slicing. The two core techniques—GEMM-oriented mode reordering and communication-aware mode distribution planning— are fully general: they depend only on the contraction tree, mode lifetimes, and tensor shapes, not on any applicationspecific property. A compute-only analysis across six diverse workloads—quantum circuit simulation, quantum error correction, combinatorial optimization on King’s subgraphs, and many-body dynamics on rectangular, hexagonal, and triangular lattices (Fig. 1)—shows that the FLOP reduction from distributed memory is consistent across application families. Per-slice execution measurements on four of these workloads sustain 28–33 TFLOP/s per GPU (42–49% of H100 FP32 peak). Within a single NVLink node, distribution yields 7.0– 173× extra speedup beyond slicing, capturing 87–101% of the compute-only reduction because NVLink’s high intra-node bandwidth leaves communication well below the compute time. Across 1024 H100 GPUs over InfiniBand—where internode bandwidth is an order of magnitude lower than NVLink—the extra speedup beyond embarrassingly parallel slicing remains large: 41.8× on Zuchongzhi n60m24, 67,869× on hexagonal many-body dynamics, 216× on rectangular many-body dynamics, and 132.6× on triangular many-body dynamics. These results establish a clear principle: distribution trades arithmetic for communication, so its benefit is limited by interconnect bandwidth rather than by the available FLOP
reduction. On NVLink, where bandwidth is high, this communication cost is small and the extra speedup tracks the compute-only reduction; on InfiniBand it becomes the dominant cost—yet in both regimes the extra speedup stays far above linear slicing. This makes the framework particularly well-suited to emerging platforms such as the NVIDIA GB200 NVL72, which extends NVLink-class bandwidth to 72 GPUs (13.8 TB aggregate memory, 1.8 TB/s per GPU), and future multi-node NVLink (MNNVL) systems that further widen the high-bandwidth domain. As interconnect technology continues to advance, a larger share of the compute reduction will translate into wall-clock gains, making communication-aware distributed contraction an increasingly powerful alternative to slicing for frontier TNs. D ISCLOSURE AI writing assistants (Claude Code and Codex) were used for drafting and refining text, generating figure-plotting scripts, and restructuring sections. All technical content, experimental data, algorithm design, and scientific claims originate from the authors. All AI-generated output was reviewed and validated by the authors. R EFERENCES [1] U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics, vol. 326, no. 1, pp. 96–192, 2011. [Online]. Available: https://doi.org/10.1016/j.aop.2010.09.012 [2] R. Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Annals of Physics, vol. 349, pp. 117–158, 2014. [Online]. Available: https: //doi.org/10.1016/j.aop.2014.06.013 [3] J. Gray and S. Kourtis, “Hyper-optimized tensor network contraction,” Quantum, vol. 5, p. 410, 2021. [Online]. Available: https://doi.org/10. 22331/q-2021-03-15-410 [4] E. M. Stoudenmire and D. J. Schwab, “Supervised learning with tensor networks,” in Advances in Neural Information Processing Systems 29, 2016, pp. 4799–4807. [5] A. Novikov, D. Podoprikhin, A. Osokin, and D. P. Vetrov, “Tensorizing neural networks,” in Advances in Neural Information Processing Systems 28, 2015, pp. 442–450. [6] A. J. Ferris and D. Poulin, “Tensor networks and quantum error correction,” Physical review letters, vol. 113, no. 3, p. 030501, 2014. [Online]. Available: https://doi.org/10.1103/PhysRevLett.113.030501 [7] J.-G. Liu, L. Wang, and P. Zhang, “Tropical tensor network for ground states of spin glasses,” Physical Review Letters, vol. 126, no. 9, p. 090506, 2021. [Online]. Available: https://doi.org/10.1103/PhysRevLett. 126.090506 [8] B. Fauseweh, “The state of quantum computing applications in challenging many-body quantum dynamics,” Nature Communications, vol. 15, p. 2123, 2024. [Online]. Available: https://doi.org/10.1038/ s41467-024-46402-9 [9] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell et al., “Quantum supremacy using a programmable superconducting processor,” Nature, vol. 574, no. 7779, pp. 505–510, 2019. [Online]. Available: https://doi.org/10.1038/s41586-019-1666-5 [10] Google Quantum AI, “Suppressing quantum errors by scaling a surface code logical qubit,” Nature, vol. 614, no. 7949, pp. 676–681, 2023. [Online]. Available: https://doi.org/10.1038/s41586-022-05434-1 [11] K. Liang, “Independent set enumeration in king’s graphs by tensor network contractions,” arXiv preprint arXiv:2505.12776, 2025. [Online]. Available: https://arxiv.org/abs/2505.12776 [12] A. D. King, A. Nocera, M. Rams, J. Dziarmaga, R. Wiersema et al., “Beyond-classical computation in quantum simulation,” Science, pp. 199–204, 2025. [Online]. Available: https://doi.org/10.1126/science. ado6285
[13] I. L. Markov and Y. Shi, “Simulating quantum computation by contracting tensor networks,” SIAM Journal on Computing, vol. 38, no. 3, pp. 963–981, 2008. [14] NVIDIA Corporation, “NVIDIA H100 Tensor Core GPU Architecture,” https://resources.nvidia.com/en-us-tensor-core, 2022. [15] ——, “cuTENSORMp: Multi-process tensor contraction library,” https: //docs.nvidia.com/cuda/cutensor/latest/user guide cutensorMp.html, 2024. [16] ——, “NCCL: NVIDIA collective communications library,” https:// developer.nvidia.com/nccl, 2024. [17] C. Huang, F. Zhang, M. Newman, X. Ni, D. Ding, J. Cai, X. Gao, T. Wang, F. Wu, G. Zhang, H.-S. Ku, Z. Tian, J. Wu, H. Xu, H. Yu, B. Yuan, M. Szegedy, Y. Shi, H.-H. Zhao, C. Deng, and J. Chen, “Efficient parallelization of tensor network contraction for simulating quantum computation,” Nature Computational Science, vol. 1, pp. 578–587, 2021. [Online]. Available: https://doi.org/10.1038/s43588-021-00119-7 [18] Y. Liu, X. Liu, F. Li, Y. Yang, J. Song, P. Zhao, Z. Wang, D. Peng, H. Fu, D. Chen, W. Wu, H. Huang, and C. Guo, “Closing the “quantum supremacy” gap: Achieving real-time simulation of a random quantum circuit using a new Sunway supercomputer,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. ACM, 2021, pp. 1–12. [Online]. Available: https://doi.org/10.1145/3458817.3487399 [19] F. Pan and P. Zhang, “Simulation of quantum circuits using the big-batch tensor network method,” Physical Review Letters, vol. 128, no. 3, p. 030501, 2022. [Online]. Available: https: //doi.org/10.1103/PhysRevLett.128.030501 [20] F. Pan, K. Chen, and P. Zhang, “Solving the sampling problem of the Sycamore quantum circuits,” Physical Review Letters, vol. 129, no. 9, p. 090502, 2022. [Online]. Available: https://doi.org/10.1103/ PhysRevLett.129.090502 [21] X.-H. Zhao, H.-S. Zhong, F. Pan et al., “Leapfrogging Sycamore: Harnessing 1432 GPUs for 7x faster quantum random circuit sampling,” National Science Review, vol. 12, no. 3, p. nwae317, 2025. [Online]. Available: https://doi.org/10.1093/nsr/nwae317 [22] R. Fu, Z. Su, H.-S. Zhong, X.-H. Zhao, J. Zhang, F. Pan et al., “Surpassing Sycamore: Achieving energetic superiority through systemlevel circuit simulation,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024. [Online]. Available: https://doi.org/10.1109/SC41406.2024. 00085 [23] Y. Chen, Z. Sun, C. Qiu, Z. Li, Y. Liu, L. Gan, X. Duan, and G. Yang, “SW-TNC: Reaching the most complex random quantum circuit via tensor network contraction,” arXiv preprint arXiv:2504.09186, 2025. [Online]. Available: https://arxiv.org/abs/2504.09186 [24] M. Xu, S. Cao, X. Miao, U. A. Acar, and Z. Jia, “Atlas: Hierarchical partitioning for quantum circuit simulation on GPUs,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC24), 2024. [Online]. Available: https://doi.org/10.1109/SC41406.2024.00072 [25] W. M. Brown, A. Ramesh, T. Lubinski, T. Nguyen, and D. E. B. Neira, “Multi-GPU quantum circuit simulation and the impact of network performance,” arXiv preprint arXiv:2511.14664, 2025. [Online]. Available: https://arxiv.org/abs/2511.14664 [26] A. Cicero, M. A. Maleki, M. W. Azhar, A. F. Kockum, and P. Trancoso, “Simulation of quantum computers: Review and acceleration opportunities,” ACM Transactions on Quantum Computing, vol. 7, no. 1, p. 3, 2025. [Online]. Available: https://doi.org/10.1145/ 3701725 [27] Y. Wu, W.-S. Bao, S. Cao, F. Chen, Y. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan et al., “Strong quantum computational advantage using a superconducting quantum processor,” Physical Review Letters, vol. 127, no. 18, p. 180501, 2021. [Online]. Available: https://doi.org/10.1103/PhysRevLett.127.180501 [28] A. Zlokapa, F. Fuchs, L. Schaeffer, A. M. Dalzell, E. Lau, E. T. Holland et al., “The computational boundaries of quantum advantage,” npj Quantum Information, vol. 9, p. 36, 2023. [Online]. Available: https://doi.org/10.1038/s41534-023-00744-7 [29] F. Pan, P. Zhou, S. Li, and P. Zhang, “Contracting arbitrary tensor networks: General approximate algorithm and applications in graphical models and quantum circuit simulations,” Physical Review Letters, vol. 125, no. 6, p. 060503, 2020. [Online]. Available: https://doi.org/10.1103/PhysRevLett.125.060503
[30] F. Pan, H. Gu, L. Kuang, B. Liu, and P. Zhang, “Efficient quantum circuit simulation by tensor network methods on modern GPUs,” ACM Transactions on Quantum Computing, vol. 5, no. 4, 2024. [Online]. Available: https://doi.org/10.1145/3696465 [31] S. Bravyi, M. Suchara, and A. Vargo, “Efficient algorithms for maximum likelihood decoding of quantum error-correcting codes,” Physical Review A, vol. 90, no. 3, p. 032326, 2014. [Online]. Available: https://doi.org/10.1103/PhysRevA.90.032326 [32] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, “Surface codes: Towards practical large-scale quantum computation,” Physical Review A, vol. 86, no. 3, p. 032324, 2012. [Online]. Available: https://doi.org/10.1103/PhysRevA.86.032324 [33] C. Piveteau, C. T. Chubb, and J. M. Renes, “Tensor-network decoding beyond 2d,” PRX Quantum, vol. 5, no. 4, p. 040303, 2024. [Online]. Available: https://doi.org/10.1103/PRXQuantum.5.040303 [34] J. Bausch, M. S. Kesselring, A. Elben, V. Swaroop, B. Yao, A. Molle et al., “Learning high-accuracy error decoding for quantum processors,” Nature, vol. 635, pp. 834–840, 2024. [Online]. Available: https://doi.org/10.1038/s41586-024-08148-8 [35] Google Quantum AI et al., “Quantum error correction below the surface code threshold,” Nature, vol. 638, pp. 920–926, 2025. [Online]. Available: https://doi.org/10.1038/s41586-024-08449-y [36] M. Suzuki, “Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems,” Communications in Mathematical Physics, vol. 51, no. 2, pp. 183–190, 1976. [Online]. Available: https: //doi.org/10.1007/BF01609348 [37] G. Vidal, “Efficient simulation of one-dimensional quantum many-body systems,” Physical Review Letters, vol. 93, no. 4, p. 040502, 2004. [Online]. Available: https://doi.org/10.1103/PhysRevLett.93.040502 [38] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, “Evidence for the utility of quantum computing before fault tolerance,” Nature, vol. 618, pp. 500–505, 2023. [Online]. Available: https://doi.org/10.1038/s41586-023-06096-3 [39] O. Shtanko, D. S. Wang, H. Zhang, N. Harle, A. Seif, R. Movassagh, and Z. Minev, “Uncovering local integrability in quantum manybody dynamics,” Nature Communications, 2025. [Online]. Available: https://doi.org/10.1038/s41467-025-57623-x [40] A. D. King et al., “Quantum critical dynamics in a 5,000-qubit programmable spin glass,” Nature, vol. 617, pp. 61–66, 2023. [41] J. Mildenberger, Z. Jiang, W. Mruczkiewicz, J. C. Halimeh, and P. Hauke, “Confinement in a Z2 lattice gauge theory on a quantum computer,” Nature Physics, 2025. [Online]. Available: https://doi.org/10.1038/s41567-024-02723-6 [42] NVIDIA Corporation, “cuTENSOR: A high-performance CUDA library for tensor primitives,” https://developer.nvidia.com/cutensor, 2024. [43] D. A. Matthews, “High-performance tensor contraction without transposition,” SIAM Journal on Scientific Computing, vol. 40, no. 1, pp. C1–C24, 2018. [Online]. Available: https://doi.org/10.1137/ 16M108968X