Tensor Network Simulations

• Tensor network simulations represent high-dimensional quantum states by contracting lower-rank tensors, enabling scalable approximations of many-body systems.
• They use heuristic and optimized contraction paths to reduce computational complexity by exploiting entanglement locality and algebraic factorization.
• Modern implementations leverage GPUs, FPGAs, and distributed systems to simulate large-scale quantum circuits while mitigating exponential memory demands.

Tensor network simulations provide a computational paradigm for representing and contracting high-rank tensors arising in quantum many-body physics, quantum circuit simulation, statistical mechanics, and related high-dimensional numerical domains. The tensor network (TN) formalism enables scalable approximation of quantum states, operators, and partition functions, exploiting entanglement locality and algebraic factorization to avoid the exponential memory and computational costs of explicit Hilbert-space representations. Modern implementations leverage architectural features of CPUs, GPUs, and FPGAs, as well as parallel and distributed systems, to evaluate TN contractions for systems exceeding the limits of brute-force state-vector methods.

1. Mathematical Formulation and Classes of Tensor Networks

Tensor networks encode a high-rank tensor (e.g., a quantum state $$|\psi\rangle = \sum_{i_1,\ldots,i_N} c_{i_1,\ldots,i_N}|i_1\ldots i_N\rangle$$) as the contraction of a network of lower-rank tensors ("nodes") connected by "bond" indices. The design of the TN graph is dictated by geometry, physical locality, or problem structure.

Key classes:

• Matrix Product States (MPS): 1D chain of rank-3 tensors, bond dimension $\chi$, representing quantum states or operators efficiently when entanglement is area-law bounded. Storage scales as $O(N\chi^2)$, operations (TEBD, DMRG) as $O(N\chi^3)$ (García et al., 2024).
• Tree Tensor Networks (TTN): Hierarchically organized, loop-free networks matching tree-like entanglement. Suitable for problems where the entanglement structure is not linear (Li et al., 2024).
• Projected Entangled Pair States (PEPS): 2D generalization for lattice systems. Each site is a rank-5 tensor connected to up to four neighbors. Exact contraction is $\#P$-hard; practical contraction uses boundary MPS or other approximations (Pang et al., 2020, García et al., 2024).
• Multiscale Entanglement Renormalization Ansatz (MERA): Layered network combining unitaries and isometries to efficiently represent critical states with power-law entanglement.

Quantum circuits can be mapped to TNs by expressing each gate as a small tensor and each wire as an index, with contractions corresponding to summing over the intermediary indices (Levental, 2021, Brennan et al., 2021).

2. Contraction Algorithms and Path Optimization

The defining computational task in TN simulations is evaluating the full contraction, typically a sum over shared indices yielding a scalar (for partition functions), a vector (for quantum states), or specific amplitudes (for quantum circuits).

Contraction path selection:

• The order of pairwise contractions ("contraction path") profoundly affects computational and memory complexity; finding the optimal path is NP-hard (Pan et al., 2023, Levental, 2021).
• Contraction is formalized via Einstein summation ($\texttt{einsum}$) notation, mapping to sequences of GEMMs (general matrix multiplies) for efficiency on modern linear algebra hardware (Pan et al., 2023).
• Heuristic methods (simulated annealing, graph partitioning, treewidth minimization, Bayesian or greedy optimization) are used in practice to find low-cost paths. The path cost models incorporate estimates of total FLOPs, peak memory, and largest intermediate tensor (Pan et al., 2023, Levental, 2021, Brennan et al., 2021).
• For 2D and higher-dimensional networks, contraction paths may use boundary-MPS techniques, implicit randomized SVD truncations, or coarse-graining approaches (e.g., TRG, HOTRG) (Pang et al., 2020, Shimizu et al., 2024).

Approximate contraction:

• Intractable contractions are rendered tractable by truncating bond dimensions (e.g., via SVD thresholding), which introduces controlled errors but keeps the simulation scalable as long as the underlying entanglement is not volume-law (McCaskey et al., 2018, Pang et al., 2020).

3. Hardware-Accelerated and Distributed Simulation

Achieving scalable tensor network simulation requires hardware-tailored algorithms:

• GPUs: Mapping critical contraction steps to GEMM operations enables effective use of GPU tensor cores. Techniques include permuting and reshaping input tensors to maximize GEMM size, mixed-precision scheduling, and careful monitoring of numerical errors (Pan et al., 2023, Schieffer et al., 27 Jan 2025). Mixed-precision (e.g., 3×TF32) and extended-precision approaches achieve high throughput while maintaining fidelity (Pan et al., 2023).
• FPGAs: Matrix-multiplies are implemented by custom systolic arrays or naive mat-mul engines; output-stationary designs reduce required bandwidth and improve performance; block RAM constraints necessitate careful network factorization and slice management (Levental, 2021).
• Distributed HPC: Task-based parallelization schemes (layer decomposition, dynamic work-stealing) and distributed tensor frameworks (e.g., TAMM, Cyclops, QuantEx) enable scaling to exascale platforms. Slicing and load balancing strategies are essential for efficiency (Hoyt et al., 7 Jan 2026, Brennan et al., 2021, Pang et al., 2020).
• CPU/GPU Hybrid: Taskflow-based libraries (e.g., Jet) manage work reuse, efficient scheduling, and memory optimization by constructing DAGs for dependency-aware contraction, yielding 2–4× speedups on single-node systems (Vincent et al., 2021).

Scalability: Strong and weak scaling is routinely demonstrated up to hundreds of nodes or GPUs for model systems (e.g., 60–90 qubit circuits, 2D PEPS on $8\times8$ grids) with practical memory usage (Hoyt et al., 7 Jan 2026, Schieffer et al., 27 Jan 2025, Pang et al., 2020, Brennan et al., 2021).

4. Leading Algorithms and Software Implementations

Ground state and time evolution: TEBD, DMRG, TDVP, and VUMPS are widely used in TN simulations for static and dynamical properties (Hauschild et al., 2018, Hauschild et al., 2024, Lacroix et al., 2024).

Observables and error estimation: Measurement observables and multi-time correlations are computed by operator insertion and contraction. Error metrics include truncation errors, L2 norm error, and fidelity-based estimations (e.g., LXEB) (Pan et al., 2023, Lacroix et al., 2024).

Software: TeNPy, ITensor, ExaTENSOR, Koala, TAMM, QuantEx, Jet, and NWQ-Sim represent a significant shared codebase, supporting CPU, GPU, and distributed platforms with extensible APIs for network topologies (MPS, TTN, PEPS), symmetry handling, and observable measurement (Hauschild et al., 2024, Hauschild et al., 2018, Vincent et al., 2021, Brennan et al., 2021, Hoyt et al., 7 Jan 2026, Pang et al., 2020). Libraries often support U(1) and abelian symmetries for block-sparse acceleration (Hauschild et al., 2018, Hauschild et al., 2024).

Randomized algorithms: Randomized SVD (R-SVD) provides order-of-magnitude speedups over deterministic SVD in DMRG, TEBD, and PEPS contraction without loss of accuracy, especially in high local dimension or critical systems (Kohn et al., 2017, Pang et al., 2020).

5. Exemplary Applications

Quantum circuit simulation

• Large-scale random quantum circuits (Sycamore, QAOA, QFT) are simulated using TN contraction, with simulations of 53–90 qubit circuits tractable using mixed-precision, hardware-optimized, and parallelized approaches (Pan et al., 2023, Schieffer et al., 27 Jan 2025, Hoyt et al., 7 Jan 2026, Brennan et al., 2021).
• Surface code error correction under arbitrary local noise is modeled by expressing code states and CPTP maps as PEPS; both exact contraction and boundary-MPS approximations can efficiently deliver threshold and subthreshold logical error probabilities (Darmawan et al., 2016).
• Efficient simulation of weakly-measured circuits is achieved by MPS+Markov chain sampling, revealing measurement-induced phase transitions in entanglement scaling (Pereira et al., 8 Oct 2025).

Open quantum system dynamics

• Non-Markovian environments are handled using chain mapping (e.g., TEDOPA, T-TEDOPA), TTN representation for system+bath, and TDVP for real-time or thermalized evolution, with full measurement and observables support (Lacroix et al., 2024, Fux et al., 2022).
• Tree topology networks with optimal operators (minimum vertex cover via bipartite graphs) are used for general sum-of-product Hamiltonians in ML-MCTDH simulations (Li et al., 2024).

Machine learning and optimization

• MPS and TTN architectures serve as classical models for supervised learning and generative modeling, attaining MNIST-scale accuracies and supporting quantum-inspired stochastic models (García et al., 2024, Pereira et al., 8 Oct 2025).
• Quantum annealing for classical neural network optimization is efficiently simulated as TN time evolution, supporting problems up to $N\sim100$ with low bond dimension (Lami et al., 2022).

6. Computational Complexity and Performance Analysis

The computational complexity of TN simulations is determined by the TN class, contraction path, and bond dimension growth:

• MPS: $O(N\chi^3)$ cost, $O(N\chi^2)$ memory. For volume-law entanglement, $\chi \sim 2^{N/2}$, rendering the simulation exponential (Schieffer et al., 27 Jan 2025, McCaskey et al., 2018).
• PEPS (2D): Approximate contraction via boundary-MPS, randomized SVD, or directional CTM handles $8\times8$–$10\times10$ grids, with cost scaling $O(n^2 m^3 D^4)$ prior to optimization, and $O(n^2 m^2 D^4 + n^2 m^3 D^2)$ under randomized refactorization (Pang et al., 2020).
• TTN/TTNS: For $k$-ary trees, cost per step is $O(S^{k+1}O^{k-1}d)$ for state and operator contraction, scaling linearly with modes for fixed bond (Li et al., 2024).

Accuracy and error control: Truncations (e.g., via SVD or power thresholding), mixed-precision operations, and systematic error analysis (L2, fidelity) serve to control simulation accuracy (Pan et al., 2023, Macmahon et al., 2018).

Performance figures:

• GPU-sustained TFLOPs: on A100, up to 21 TFLOPS with mixed-precision TN simulation, reducing verification time by 4–12× compared to state-of-the-art CPU or older GPU results (Pan et al., 2023).
• Parallelized brickwork circuits (depth=80, $\chi=250$): strong scaling up to 4 nodes for 80-qubit MPS contraction; weak scaling shown up to 256 nodes for 2D PEPS (Hoyt et al., 7 Jan 2026, Pang et al., 2020).

7. Limitations and Prospects

Intrinsic limitations:

• The dependence of simulation cost on entanglement entropy: volume-law growth leads to exponential cost regardless of network topology (McCaskey et al., 2018, Schieffer et al., 27 Jan 2025, Pang et al., 2020).
• Contraction-order optimization remains an NP-hard bottleneck; heuristic methods are necessary for large and complex graphs (Pan et al., 2023, Levental, 2021).
• PEPS and higher-dimensional networks present significant scaling barriers; while randomized and distributed methods alleviate cost, approximations introduce trade-offs in accuracy (Pang et al., 2020).
• Benchmarking and cross-software reproducibility are complicated by varied interfaces, choices of TN topology, and contraction schemes (García et al., 2024).

Emerging directions:

• Automated contraction path optimization integrated with hardware-aware scheduling (Pan et al., 2023, Hoyt et al., 7 Jan 2026).
• Batched, fused, and mixed-precision kernels, as well as parallel-SVD backends for $\chi \gg 10^3$ (Schieffer et al., 27 Jan 2025, Hoyt et al., 7 Jan 2026).
• Expanding TN classes (PEPS, TTN, MERA, fermionic/anyonic TNs) and hybridization with classical/quantum resources.
• Application-driven innovations: quantum error correction, quantum optimization, high-dimensional machine learning, and nontrivial manifold (e.g., non-orientable surface) simulations (Darmawan et al., 2016, Pereira et al., 8 Oct 2025, Shimizu et al., 2024).

Tensor network simulations have proven to be a universal tool for the compression and manipulation of large quantum states, providing practical—and in many cases asymptotically optimal—scalability for a wide spectrum of problems, limited essentially by the underlying entanglement and tensor contraction complexity (García et al., 2024, Pan et al., 2023, Brennan et al., 2021).