Tensor Network Methods

July 8, 2026 · View on GitHub

Two modules for tensor-network-based simulation of the Kuramoto-XY system:

MPS/DMRG/TEBD (phase/mps_evolution.py) — ground state and time evolution via quimb
Contraction optimiser (phase/contraction_optimiser.py) — optimal einsum paths via cotengra

These extend the package beyond exact diagonalisation limits to $n = 32$ –64 qubits (system-size-dependent on bond dimension and entanglement).

Caveat: Both modules require optional dependencies (quimb, cotengra). Install via pip install scpn-quantum-control[quimb] or pip install scpn-quantum-control[cotengra].

The QWC-4.2 N=30-40 baseline design is tracked in TN/MPS Baseline Design and regenerated with scpn-bench s2-tn-mps-baseline-design. It is a CPU-first plan, not a measured N=40 tensor-network benchmark.

Part 1: MPS Evolution

scpn_quantum_control.phase.mps_evolution

Theory

Matrix Product States

An MPS represents an $n$ -qubit state as a chain of tensors:

$|\psi\rangle = \sum_{s_1 \ldots s_n} A^{s_1} A^{s_2} \cdots A^{s_n} |s_1 \ldots s_n\rangle$

where each $A^{s_i}$ is a matrix of dimension $\chi \times \chi$ ( $\chi$ = bond dimension). Memory scales as $O(n \chi^2 d)$ instead of $O(d^n)$ for exact states.

MPS is efficient for states with bounded entanglement (area-law states). Ground states of 1D gapped Hamiltonians satisfy this. Highly entangled states (e.g., volume-law after long time evolution) require large $\chi$ .

DMRG

Density Matrix Renormalisation Group finds the ground state by variationally optimising the MPS tensors site-by-site. Converges to the true ground state for 1D systems with moderate entanglement.

TEBD

Time-Evolving Block Decimation propagates $|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle$ by Trotterising the evolution into nearest-neighbour gates and applying them to the MPS, truncating the bond dimension after each gate.

Limitations

Nearest-neighbour only: quimb's SpinHam1D supports NN couplings. Direct DMRG/TEBD calls reject longer-range terms from $K_{nm}$ unless allow_long_range_truncation=True is passed. Truncated runs report coupling_scope="nearest_neighbour_truncated" and omitted_coupling_l1 so they cannot be mistaken for full-K simulations.
Bond dimension: Higher $\chi$ → more accurate but slower. Typical values: $\chi = 32$ –128 for ground states, $\chi = 64$ –256 for dynamics.
Entanglement growth: TEBD accuracy degrades at long times as entanglement grows beyond what the bond dimension can capture.

API Reference

from scpn_quantum_control.phase.mps_evolution import (
    is_quimb_available,
    dmrg_ground_state,
    tebd_evolution,
)

`dmrg_ground_state`

result = dmrg_ground_state(
    K: np.ndarray,          # (n, n) coupling matrix
    omega: np.ndarray,      # (n,) natural frequencies
    bond_dim: int = 64,     # maximum bond dimension
    cutoff: float = 1e-10,  # SVD truncation cutoff
    max_sweeps: int = 20,   # DMRG sweep limit
) -> dict

Returns:

{
    "energy": float,         # ground state energy
    "mps": quimb.MatrixProductState,  # ground state MPS
    "converged": bool,       # True if energy converged
    "bond_dims": list[int],  # bond dimensions along the chain
}

`tebd_evolution`

result = tebd_evolution(
    K: np.ndarray,
    omega: np.ndarray,
    t_max: float = 1.0,     # total evolution time
    dt: float = 0.05,       # Trotter step size
    bond_dim: int = 64,     # maximum bond dimension
    cutoff: float = 1e-10,  # SVD truncation
    order: int = 2,         # Trotter order (2 or 4)
) -> dict

Returns:

{
    "times": np.ndarray,          # time points
    "R": np.ndarray,              # order parameter R(t)
    "bond_dims_final": list[int], # bond dimensions at final time
    "mps_final": quimb.MatrixProductState,  # final MPS
}

`is_quimb_available`

is_quimb_available() -> bool

Returns True if quimb is installed.

Tutorial

import numpy as np
from scpn_quantum_control.phase.mps_evolution import (
    dmrg_ground_state, tebd_evolution, is_quimb_available
)

if not is_quimb_available():
    raise ImportError("Install quimb: pip install quimb")

n = 16
K = 0.45 * np.exp(-0.3 * np.abs(np.subtract.outer(range(n), range(n))))
np.fill_diagonal(K, 0.0)
omega = np.linspace(0.8, 1.2, n)

# Ground state
gs = dmrg_ground_state(K, omega, bond_dim=64, max_sweeps=20)
print(f"DMRG energy: {gs['energy']:.6f}")
print(f"Converged: {gs['converged']}")
print(f"Bond dimensions: {gs['bond_dims']}")

# Time evolution
dyn = tebd_evolution(K, omega, t_max=2.0, dt=0.05, bond_dim=64, order=2)
print(f"R(0) = {dyn['R'][0]:.3f}")
print(f"R(T) = {dyn['R'][-1]:.3f}")

Bond Dimension Selection

System	Recommended $\chi$	Memory per MPS
$n \leq 16$ ground state	32–64	< 1 MB
$n = 32$ ground state	64–128	~10 MB
$n = 64$ ground state	128–256	~100 MB
Short-time dynamics ( $t < 1$ )	64–128	~10 MB
Long-time dynamics ( $t > 5$ )	256+	~1 GB

Part 2: Contraction Optimiser

scpn_quantum_control.phase.contraction_optimiser

Theory

Tensor contraction order affects computational cost dramatically. For a chain of matrix multiplications $A \cdot B \cdot C$ , the parenthesisation matters: $(AB)C$ vs $A(BC)$ can differ by orders of magnitude in flops.

cotengra (Gray & Kourtis, 2021) finds near-optimal contraction paths for arbitrary tensor networks using hyper-optimisation over graph partitioning algorithms.

This module provides a drop-in np.einsum replacement that uses cotengra when available, with numpy fallback.

API Reference

from scpn_quantum_control.phase.contraction_optimiser import (
    is_cotengra_available,
    optimal_contraction_path,
    contract,
    benchmark_contraction,
)

`contract`

result = contract(
    subscripts: str,        # einsum subscripts, e.g. "ij,jk->ik"
    *operands: np.ndarray,  # input tensors
    optimiser: str = "auto",  # "auto", "greedy", "optimal"
) -> np.ndarray

Drop-in replacement for np.einsum. Uses cotengra path optimisation when available.

`optimal_contraction_path`

path, info = optimal_contraction_path(
    subscripts: str,
    *operands: np.ndarray,
    optimiser: str = "auto",
) -> tuple[list[tuple[int, ...]], dict]

Returns the contraction path and metadata (flops, memory estimates) without performing the contraction.

`benchmark_contraction`

result = benchmark_contraction(
    subscripts: str,
    *operands: np.ndarray,
    n_repeats: int = 10,
) -> dict

Returns:

{
    "naive_ms": float,       # wall time with default path (ms)
    "optimised_ms": float,   # wall time with cotengra path (ms)
    "speedup": float,        # naive / optimised
}

Example

import numpy as np
from scpn_quantum_control.phase.contraction_optimiser import (
    contract, benchmark_contraction, is_cotengra_available
)

print(f"cotengra available: {is_cotengra_available()}")

# Simple matrix chain
A = np.random.randn(100, 200)
B = np.random.randn(200, 300)
C = np.random.randn(300, 50)

# Drop-in replacement for np.einsum
result = contract("ij,jk,kl->il", A, B, C)
print(f"Result shape: {result.shape}")  # (100, 50)

# Benchmark
bench = benchmark_contraction("ij,jk,kl->il", A, B, C, n_repeats=20)
print(f"Naive: {bench['naive_ms']:.1f} ms")
print(f"Optimised: {bench['optimised_ms']:.1f} ms")
print(f"Speedup: {bench['speedup']:.1f}×")

Comparison

Feature	This module (quimb)	quimb standalone	ITensor	TeNPy
DMRG	Yes	Yes	Yes	Yes
TEBD	Yes	Yes	Yes	Yes
2D PEPS	No	Yes	Yes	No
Hamiltonian	Kuramoto-XY (built-in)	Any	Any	Any
Language	Python	Python	Julia/C++	Python
cotengra	Yes (via contraction_optimiser)	Yes	No	No

Our modules wrap quimb with Kuramoto-XY-specific Hamiltonian construction and order parameter extraction. They are not a quimb replacement.

References

White, S. R. "Density matrix formulation for quantum renormalization groups." PRL 69, 2863 (1992). (DMRG)
Vidal, G. "Efficient simulation of one-dimensional quantum many-body systems." PRL 93, 040502 (2004). (TEBD)
Gray, J. & Kourtis, S. "Hyper-optimized tensor network contraction." Quantum 5, 410 (2021). (cotengra)
Gray, J. "quimb: A Python library for quantum information and many-body calculations." JOSS 3, 819 (2018).