GPU-Accelerated Tensor Contractions

• GPU-accelerated tensor contractions are computational methods that transform multilinear summations into GPU-efficient kernels like GEMM and Tensor Core instructions.
• They employ advanced techniques such as contraction-order optimization, memory-layout control, and communication-aware scheduling to tackle bottlenecks from transposes, decompositions, and bandwidth limitations.
• Applications in finite element methods, tensor networks, and sparse tensor operations demonstrate significant performance gains, with speedups reaching up to two orders of magnitude.

Searching arXiv for papers on GPU-accelerated tensor contractions and related tensor-network/GPU implementations. GPU-accelerated tensor contractions are computational methods that map multilinear summations, reductions, and tensor-network contractions onto graphics processing units in order to raise throughput, reduce data-movement overheads, and extend feasible problem sizes. Across finite element methods, tensor renormalization, projected entangled-pair states, sparse tensor decomposition, quantum-circuit tensor networks, and distributed contraction engines, the central technical theme is the reformulation of contractions into GPU-efficient kernels such as GEMM, batched GEMM, GEMV, or Tensor Core instructions, together with contraction-order optimization, memory-layout control, and communication-aware scheduling (Cui, 2024, Zhang et al., 1 May 2025, Pan et al., 1 Jun 2026). The literature emphasizes that performance is often limited not only by floating-point throughput but also by transposes, decompositions, atomics, memory bandwidth, and inter-GPU communication, so successful implementations typically combine algebraic reformulation with architecture-specific execution strategies.

1. Tensor contractions as a GPU workload

Tensor contractions generalize matrix multiplication by summing over one or more shared indices. In the single-index case, one representative form is

$$C_{\mathcal{C}} \leftarrow \alpha \sum_k A_{\mathcal{A}} B_{\mathcal{B}} + \beta C_{\mathcal{C}},$$

with free-index set $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$ (Shi et al., 2016). In sparse tensor decomposition, sparse MTTKRP is written as

$$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$

and elementwise as a sum over nonzeros multiplied by rows of dense factor matrices (Wijeratne et al., 2024). In dense tensor-vector contraction, the $k$-mode TVC is defined elementwise by contraction over index $j$, and after matricization reduces either to one large matrix-vector multiplication or to independent vector-matrix multiplies (Martinez-Ferrer et al., 6 Jan 2025).

Quantum-circuit tensor networks provide another canonical formulation. A single-qubit gate contraction is

$$T'_{…\,i'_k\,…} = \sum_{i_k =0}^1 G_{i'_k,i_k}\;T_{…\,i_k\,…},$$

while a two-qubit gate contraction is

$$T'_{…\,i'_a,i'_b\,…} = \sum_{i_a,i_b=0}^1 G_{i'_a,i'_b\,;\,i_a,i_b}\;T_{…\,i_a,i_b\,…}$$

(Xiao et al., 15 Apr 2026). In tensor renormalization, higher-order tensor renormalization group (HOTRG) coarse-graining is expressed by explicit multi-tensor contractions such as

$$M_{ijkl}=\sum_{a,b,c,d,e,f} T_{\,jabe}\,T_{\,iecd}\,T_{\,labf}\,T_{\,kfcd},$$

followed by SVD-based truncation and further four-way contractions (Jha et al., 2023).

These formulations share a structural property: the arithmetic can often be expressed as matrix multiplication after reshaping and index grouping. This suggests why GPU implementations repeatedly target GEMM-like kernels, since such kernels provide high utilization when tensor modes can be ordered so that contraction and free indices map cleanly to matrix dimensions (Shi et al., 2016, Pan et al., 1 Jun 2026).

2. Kernelization, data layout, and Tensor Core execution

A recurrent strategy in GPU tensor contraction is to reduce contractions to GEMM or batched GEMM while avoiding explicit copy and transpose operations. The STRIDEDBATCHEDGEMM primitive was introduced precisely to perform a batch of small GEMMs in one call, with a constant stride between consecutive matrices in the batch, thereby avoiding explicit out-of-place copies or transposes of tensor data (Shi et al., 2016). In that formulation, each tensor contraction is mapped mechanically to GEMM arguments $(M,N,K)$ plus batch and stride parameters, and the GPU implementation organizes work by a three-dimensional grid whose third axis is the batch dimension (Shi et al., 2016).

High-level frameworks can exploit the same principle. In HOTRG with PyTorch and CUDA, no custom CUDA kernels were written; contractions were dispatched through the opt_einsum “EinsumPlanner” backend to PyTorch on CUDA, while SVD, reshape, and permute used torch, and the contraction-path optimizer minimized FLOPs before calling cuBLAS under the hood (Jha et al., 2023). Memory movement was handled by PyTorch, and in practice the entire coarse-graining loop lived on the GPU (Jha et al., 2023).

At the other end of the optimization spectrum, finite element tensor-product operators were implemented directly with NVIDIA Tensor Cores using two programming approaches: the C++ Warp Matrix Functions in nvcuda::wmma and the inline Parallel Thread Execution instructions mma.sync.aligned (Cui, 2024). That work placed significant focus on versatile inline PTX instructions combined with a conflict-free shared memory access pattern, identified as a key to unlocking superior performance (Cui, 2024). On the NVIDIA A100 GPU Tensor Cores, this yielded a 2.3-fold increase in double precision performance over traditional CUDA Cores, reaching 8 TFLOPS/s, corresponding to 45% of the theoretical maximum (Cui, 2024).

For tensor-network simulation in CUDA-Q, cuTensorNet provided automatic tensor-network graph construction, path finding, fusion of small-tensor contractions into high-utilization GEMM calls, and low-level kernel launches, while cuBLAS supplied high-throughput matrix-multiply primitives when cuTensorNet fused a contraction into batched GEMM (Xiao et al., 15 Apr 2026). All tensor legs were stored in GPU DRAM with a compact, stride-ordered layout, column-major for compatibility with cuBLAS, and cuTensorNet allocated a reusable pool of intermediate buffers sized to the maximum peak of any contraction in the chosen path (Xiao et al., 15 Apr 2026). This design eliminates per-step cudaMalloc and cudaFree, thereby avoiding fragmentation and overhead (Xiao et al., 15 Apr 2026).

These implementations collectively show that GPU acceleration depends as much on index ordering, shared-memory behavior, and avoidance of transposes as on raw compute capability. A plausible implication is that the most effective “tensor contraction” on GPU is often one that has first been transformed into a dense linear algebra kernel with favorable memory access patterns.

3. Bottlenecks beyond FLOPs: decompositions, atomics, and bandwidth

Although contractions are often the nominal target, several studies identify other operations as the true bottlenecks. In iPEPS contraction with standard corner transfer matrix renormalization group (CTMRG), the computational bottleneck typically lies in the singular value or eigenvalue decompositions involved in the renormalization step; on modern GPUs, tensor-matrix multiplies reach hundreds of Tflop/s, whereas SVD/EVD kernels remain comparatively slow and can dominate runtime, up to 99% of wall time in standard CTMRG (Zhang et al., 1 May 2025). The proposed remedy was to combine CTMRG with QR decompositions, replacing the enlarged-corner EVD/SVD with a reduced corner and a GPU-friendly QR decomposition (Zhang et al., 1 May 2025).

In sparse MTTKRP, the central challenges were different: eliminating global atomic operations across GPU thread blocks, avoiding communication of intermediate values between GPU thread blocks and GPU global memory, and ensuring balanced workload distribution across GPU thread blocks (Wijeratne et al., 2024). The proposed algorithm partitioned nonzero elements by output-mode index so that one thread block handled each partition; because only that block updated the corresponding rows of the output factor matrix, global atomic operations became unnecessary (Wijeratne et al., 2024). Intermediate partial products were accumulated in per-thread registers and/or shared memory, so no inter-block exchange of partial sums was required (Wijeratne et al., 2024).

Dense tensor-vector contraction exposes still another limiting regime. TVC was described as the most memory-bound operation of its class, with arithmetic intensity of approximately 1–2 FLOP/byte, making memory bandwidth rather than peak FLOPs the governing resource (Martinez-Ferrer et al., 6 Jan 2025). On three architectures featuring multi-core and accelerators, dTVC and dHOPM3 remained relatively close to peak system memory bandwidth, specifically 50%–80% depending on the architecture, and on par with STREAM benchmark figures (Martinez-Ferrer et al., 6 Jan 2025). On GPUs using cuBLAS, performance varied between 20% and 90% of peak depending on batch sizes and contraction mode, with occupancy limitations for fine one-dimensional splits (Martinez-Ferrer et al., 6 Jan 2025).

The finite element Tensor Core results illustrate a different precision-performance trade-off. In half precision, numerical experiments demonstrated a fourfold enhancement in solving the Poisson equation using the flexible GMRES method, preconditioned by a multigrid method in 3D, while maintaining the same discretization error as observed in double precision computations (Cui, 2024). This directly ties numerical fidelity to contraction kernel design rather than treating precision selection as an isolated concern.

4. Contraction-order optimization and structure-aware reformulation

Contraction order is a primary determinant of performance and memory use. In GPU tensor-network simulation for quantum molecular generation, cuTensorNet used a heuristic path optimizer inspired by Markov & Shi (2008) to choose the sequence of pairwise contractions, estimate the peak memory of each intermediate, and minimize total FLOPs under the available GPU memory cap (Xiao et al., 15 Apr 2026). Its built-in path finder used a Greedy + Branch-and-Bound heuristic: it first scored all possible pairwise contractions by estimated symbolic FLOP count and memory usage, greedily picked low-cost edges, and backtracked when a peak-memory estimate exceeded the device limit (Xiao et al., 15 Apr 2026). No user-tuned contraction order was required; the default cuTensorNet settings sufficed to keep intermediate tensors small, with typical rank at most 12 for SQMG (Xiao et al., 15 Apr 2026).

For multi-GPU exact tensor-network contraction, the contraction path was assumed fixed, but the problem became the transformation of that path into a communication-efficient distributed schedule (Pan et al., 1 Jun 2026). The workflow consisted of GEMM-oriented mode reordering and communication-aware mode distribution planning, executed by cuTENSORMp with automatic communication-compute overlap (Pan et al., 1 Jun 2026). The backward-pass mode-reordering algorithm permuted each tensor’s modes so that every pairwise contraction became a contiguous matrix-matrix multiplication in memory, eliminating runtime reshapes and transposes and storing tensors in row-major order as “[retained | reduced]”, with longest-lived modes first (Pan et al., 1 Jun 2026). A dynamic-programming planner then selected whether to activate, keep, redistribute, or gather distributed modes along a use-chain of large intermediates, based on a cost model combining per-step GEMM time and communication time (Pan et al., 1 Jun 2026).

Structure-aware reformulation also appears in iPEPS contraction. In QR-CTMRG, instead of absorbing a full extra row and column and building an enlarged corner of size $(\chi D^2)\times(\chi D^2)$, the method constructs a reduced corner $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$0 of size $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$1 and factors it as

$$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$2

where $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$3 is automatically an isometry because $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$4 already has rank at most $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$5 (Zhang et al., 1 May 2025). The updated corner and edge tensors are then formed by

$$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$6

This reduces the decomposition cost to $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$7 with a small prefactor and excellent GPU efficiency (Zhang et al., 1 May 2025).

A related reformulation underlies STRIDEDBATCHEDGEMM. There, contractions are reduced either to a single sb_gemm or to nested sb_gemm calls, with outer loops over batch indices and direct mapping from tensor modes to matrix dimensions and stride parameters (Shi et al., 2016). In all of these cases, the decisive step is not merely “running on a GPU” but reorganizing algebraic structure so the GPU sees contiguous, reusable, and large-enough matrix operations.

5. Distributed and multi-GPU contraction

Single-GPU acceleration does not remove memory-capacity limits, particularly for exact tensor-network contractions and high-order dense tensors. Distributed execution therefore introduces communication as a first-class constraint. In the multi-GPU framework for large-scale tensor-network contraction, slicing was identified as the dominant parallelization strategy but one that scales exponentially and incurs redundant computation (Pan et al., 1 Jun 2026). The alternative was to distribute intermediate tensors across devices with explicit communication, thereby reducing total FLOPs by a factor of $$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$8 ideally, while introducing communication whose cost depends on network bandwidth and latency (Pan et al., 1 Jun 2026).

The framework formalized each pairwise contraction step as

$$\mathcal{C} = (\mathcal{A} \cup \mathcal{B}) \setminus (\mathcal{A} \cap \mathcal{B})$$9

then estimated compute time as

$$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$0

and redistribution communication time by a bandwidth term plus a latency term involving tensor size, number of GPUs, contiguous memory blocks, block size, and per-message latency (Pan et al., 1 Jun 2026). The end-to-end projected time for one sliced subtask was written as

$$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$1

and the full contraction time as

$$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$2

with speedup and extra speedup beyond ideal slicing defined accordingly (Pan et al., 1 Jun 2026).

Within a single DGX H100 node with 8 GPUs connected by NVLink, this distribution approach delivered 7–173× extra speedup beyond embarrassingly parallel slicing and captured 87–101% of the available compute reduction because NVLink’s high bandwidth kept communication small relative to compute (Pan et al., 1 Jun 2026). On 1024 H100 GPUs over InfiniBand, the extra speedup beyond slicing ranged from 42× to 67,869× across four workloads (Pan et al., 1 Jun 2026). Sustained single-node performance reached 28.1, 31.6, 32.9, and 31.8 TFLOP/s per GPU on the Zuchongzhi n60m24 circuit, hexagonal 8×8 lattice, rectangular 49×20 lattice, and triangular 49×24 lattice respectively (Pan et al., 1 Jun 2026).

Distributed tensor-vector contraction exposes a complementary pattern. When the contraction mode does not coincide with the tensor splitting dimension, each rank performs a local contraction and the global result is formed by disjoint union; when the contraction mode equals the splitting dimension, each rank computes a local contraction and the result requires an MPI_Allreduce (Martinez-Ferrer et al., 6 Jan 2025). The associated dHOPM3 algorithm used three buffers to reduce streamed memory and required one MPI_Allreduce or MPI_Allgather on the converged vector after each mode (Martinez-Ferrer et al., 6 Jan 2025). Communication among GPUs used NCCL to exchange partial vectors, while TVC itself mapped to cuBLAS strided-batched GEMV (Martinez-Ferrer et al., 6 Jan 2025).

These results indicate that multi-GPU tensor contraction is fundamentally a scheduling problem in which contraction order, tensor layout, distribution choices, and network topology jointly determine efficiency. This suggests that future gains may depend less on isolated kernel speedups than on integrated planners that reason simultaneously about compute reduction and communication overhead.

6. Representative application domains and empirical results

The breadth of application domains illustrates that GPU-accelerated tensor contractions are not a single technique but a family of domain-specific implementations.

Domain	Method	Representative result
Finite element methods	Tensor-product operators on NVIDIA A100 Tensor Cores	2.3-fold increase in double precision performance, 8 TFLOPS/s, 45% of theoretical maximum (Cui, 2024)
Quantum molecular generation	GPU tensor-network simulation with cuTensorNet	Exact simulation extended to $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$3 heavy atoms (Xiao et al., 15 Apr 2026)
iPEPS / CTMRG	QR-based renormalization on H100	Up to two orders of magnitude speedup compared to standard CTMRG (Zhang et al., 1 May 2025)
HOTRG / TRG	PyTorch + CUDA + opt_einsum	At $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$4, single GPU about 8× faster than CPU (Jha et al., 2023)
Sparse tensor decomposition	spMTTKRP on GPU	Geometric mean speedup of 1.5×, 2.0×, and 21.7× versus three baselines (Wijeratne et al., 2024)
Multi-GPU tensor networks	Communication-aware distributed contraction	42× to 67,869× extra speedup beyond slicing on 1024 H100 GPUs (Pan et al., 1 Jun 2026)

In finite element methods, tensor-product operations were accelerated using Tensor Cores for operators with tensor products, with half-precision Poisson solves showing a fourfold enhancement under FGMRES with multigrid preconditioning while maintaining the same discretization error as in double precision (Cui, 2024). In quantum molecular generation, the state-vector simulator on GPU achieved speeds of up to $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$5 over the state-vector CPU baseline at $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$6 heavy atoms, while the GPU tensor-network simulator achieved up to $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$7 over the same baseline and remained viable to $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$8, where state-vector methods became memory-limited (Xiao et al., 15 Apr 2026).

In iPEPS, the QR-based scheme achieved up to 100× speedup in pure contraction for the Heisenberg model at $$\widetilde U^{(d)}=\mathcal{X}_{(d)}\Bigl(U^{(N-1)}\odot\cdots\odot U^{(d+1)}\odot U^{(d-1)}\odot\cdots\odot U^{(0)}\Bigr),$$9, $k$0, reducing step time on H100 from about $k$1 for standard CTMRG to about $k$2 for QR-CTMRG, and delivered 30–50× speedups per optimization step in end-to-end automatic-differentiation-based energy optimization up to $k$3 (Zhang et al., 1 May 2025). For the $k$4-$k$5 model at $k$6, one hour on an H100 GPU yielded variational energy $k$7 at $k$8, $k$9, and $j$0 (Zhang et al., 1 May 2025).

In HOTRG, runtime on a two-dimensional generalized XY model scaled as approximately $j$1 on CPU and empirically as approximately $j$2 on one RTX 2080 Ti GPU, with about 12× speedup using two RTX 2080 Ti cards and an additional approximately 1.3× with four Titan RTX GPUs (Jha et al., 2023). For the Ising model at criticality, an A100 GPU outperformed an RTX 2080 and a 4× CPU baseline across bond dimensions from 84 to 109, with runtimes explicitly reported in seconds (Jha et al., 2023).

In sparse MTTKRP, the proposed approach was described as the only GPU implementation that can support tensors with modes greater than 4, because state-of-the-art works had implementation constraints for large numbers of modes (Wijeratne et al., 2024). The workload-balancing procedure sorted vertices in each mode by degree and assigned them cyclically to $j$3 partitions, guaranteeing that no partition overloaded by more than $j$4 of the perfectly balanced ideal (Wijeratne et al., 2024). Dynamic remapping overhead ranged from 5% to 35% of total time (Wijeratne et al., 2024).

7. Design principles, limitations, and common misconceptions

Several design principles recur across the literature. One is that contraction-path quality and data layout are often as important as kernel choice. The HOTRG study states that choosing a good contraction order is as crucial on the GPU as on the CPU (Jha et al., 2023), while the multi-GPU tensor-network framework begins from a fixed path and still extracts major gains through mode reordering and distribution planning (Pan et al., 1 Jun 2026). Another is that specialized libraries provide substantial leverage: cuTensorNet and cuBLAS in tensor-network simulation (Xiao et al., 15 Apr 2026), cuTENSORMp in distributed contraction (Pan et al., 1 Jun 2026), and opt_einsum with PyTorch in HOTRG (Jha et al., 2023).

A common misconception is that GPU tensor contraction performance is determined chiefly by nominal tensor-core throughput. The surveyed work shows instead that performance may be dominated by decomposition kernels in CTMRG (Zhang et al., 1 May 2025), by memory bandwidth in TVC (Martinez-Ferrer et al., 6 Jan 2025), by atomics and irregularity in sparse MTTKRP (Wijeratne et al., 2024), or by communication in multi-GPU execution over InfiniBand (Pan et al., 1 Jun 2026). Another misconception is that slicing is the natural or sufficient scaling strategy for multi-GPU tensor networks. The distributed framework explicitly characterizes slicing as incurring redundant computation and shows that communication-aware distribution can far surpass slicing-based scaling limits (Pan et al., 1 Jun 2026).

The limitations are likewise domain-specific. STRIDEDBATCHEDGEMM does not directly cover all single-index contractions; eight out of thirty-six single-index contractions between a second-order and third-order tensor were identified as exceptional cases requiring an extended operation mode (Shi et al., 2016). For very large matrices, standard flattened GEMM can outperform strided batched GEMM, motivating heuristics that switch between the two (Shi et al., 2016). In HOTRG, device RAM remains the dominant bottleneck, and very large bond dimensions can still trigger out-of-memory conditions despite splitting heuristics (Jha et al., 2023). In distributed TVC and dHOPM3 on GPU, performance drops below 10% of peak at 128 ranks because of small-vector collectives and occupancy issues (Martinez-Ferrer et al., 6 Jan 2025). In sparse MTTKRP, the published implementation was tuned for one GPU architecture, and block dimensions may require retuning on other devices (Wijeratne et al., 2024).

Taken together, these studies define GPU-accelerated tensor contractions as a mature but heterogeneous research area in which algebraic restructuring, memory-system awareness, precision management, and topology-aware parallel scheduling are inseparable. The consistent outcome is that tensor contractions become substantially more tractable on modern GPUs when reformulated to minimize transposes, atomics, decompositions, and communication while exposing regular GEMM-like structure or Tensor Core-compatible microkernels (Cui, 2024, Shi et al., 2016, Pan et al., 1 Jun 2026).