mps-benchmarks

A standalone, reproducible benchmark comparing ground-state DMRG across four tensor-network libraries:

library	language	role
MPSKit	Julia	two-site DMRG2
ITensorMPS	Julia	two-site dmrg
TenPy	Python	TwoSiteDMRGEngine
YASTN	Python	two-site dmrg_

Every library solves the identical problem — same model, system size, symmetry sector, initial product state, bond-dimension schedule, and sweep count — and reports the full per-sweep trajectory (wall time, energy, bond dimension, and, where the library exposes it, local-eigensolver matvec count). The design follows the convention shared by the community's dedicated benchmark repos (tenpy_benchmarks, ITensorBenchmarks.jl, yastn_benchmarks): a per-library subdirectory with a native runner, a shared model/config layer, committed raw results, and a portable plotting step.

Results

Figures from a full run on one Flatiron gen node (12 BLAS threads per arm, nsweeps = 16): MPSKit main (0.13.13) · ITensorMPS 0.4.1 · TenPy 1.1.0 · YASTN 1.6.2. The per-symmetry χ grid runs the dense cases to 768 and the block-sparse symmetric cases to 2048.

For each case: the first figure is the default arms — energy error vs cumulative wall time (time-to-accuracy); the second is the fixed ncv=3 arms — single-sweep wall time vs χ (raw per-sweep cost, matched eigensolver work). BLAS caveat: MPSKit/ITensor route through MKL, TenPy/YASTN's numpy through OpenBLAS — the thread count is matched (12) but the backend differs, so absolute Julia-vs-Python wall times carry a backend difference (see Fairness protocol). Energy is a correctness gate, not the headline: at the largest χ the well-conditioned cases agree across all six arms to ≲1e-8.

Spin-1 Heisenberg chain, L=100 (canonical_s1)

The community-canonical DMRG benchmark. U(1) (Sᶻ), χ→2048:

Dense (no symmetry), χ→768:

Spin-1/2 Heisenberg chain, L=50 (heisenberg_s12)

Gapless/critical, so the bond dimension keeps mattering up the whole schedule. U(1) (Sᶻ), χ→2048:

Dense (no symmetry), χ→768:

Half-filled Hubbard chain, L=50 (hubbard)

Mott-gapped — the eigensolver acid test. U(1)×U(1) (N and Sᶻ), χ→2048:

Fermion parity only (trivial), χ→768:

Headline finding. On the Mott-gapped Hubbard case the fixed ncv=3 arms (ITensor, YASTN, and the *_fixed arms) plateau above the ground state even at χ=2048, while the adaptive defaults (MPSKit-default, TenPy-default) converge to it — exactly the eigensolver-settings effect this benchmark is built to expose. The gapped/critical Heisenberg cases converge identically across all six arms, so the difference is specific to the hard, near-degenerate local problems of the Mott gap.

Comparison arms

Each library runs at its out-of-the-box DMRG default, plus (where the library allows it) a fixed ITensor-like Krylov arm (krylovdim = 3, maxiter = 1) so the eigensolver-settings axis is exposed:

arm	library	eigensolver
mpskit_default	MPSKit	default (adaptive)
mpskit_fixed	MPSKit	fixed Krylov dim 3, 1 iteration, tol 1e-14
tenpy_default	TenPy	default (ED for small blocks, adaptive Lanczos above)
tenpy_fixed	TenPy	fixed Lanczos N_min = N_max = 3
itensor	ITensorMPS	default — already krylovdim 3 (= the fixed setting)
yastn	YASTN	default — already ncv = 3 (= the fixed setting)

ITensorMPS and YASTN have a single arm each: their DMRG default local eigensolver is the fixed ncv-3 Krylov solve the others pin their *_fixed arm to, so each serves as both the library default and an external fixed-Krylov reference. (This was verified against YASTN's dmrg_ docstring, whose default opts_eigs is {'hermitian': True, 'ncv': 3, 'which': 'SR'}.)

Cases

A case is a physical model; each is run in a dense (no-symmetry) and a symmetric variant.

case	model	L	symmetries
canonical_s1	spin-1 Heisenberg XXX (Haldane-gapped)	100	dense, U(1) Sz
heisenberg_s12	spin-1/2 Heisenberg XXX (gapless)	50	dense, U(1) Sz
hubbard	half-filled Hubbard, t=1, U=8, μ=U/2 (Mott)	50	Z2 fermion parity, U(1)×U(1) N⊗Sz

canonical_s1 is the community-canonical DMRG benchmark — spin-1 Heisenberg, L=100, two-site DMRG, sweeping χ, dense and U(1) — used verbatim by the ITensor↔TeNPy and YASTN comparisons, so its numbers are directly comparable to published results. Cases are declared in cases/*.yml and read identically by all four runners.

Fairness protocol

• Same problem: all libraries start from the same bond-dimension-1 product state, pinned to the ground-state charge sector for the symmetric runs (Néel, total Sz = 0; Hubbard: one electron per site ⇒ Nf = L, Sz = 0), so the symmetric and dense runs converge to the same energy. The MPSKit / ITensor / TenPy / YASTN model conventions are matched exactly (Hubbard at the particle-hole-symmetric μ = U/2).
• Same schedule: a matched χ grid and a fixed sweep count (nsweeps) for every library, so all four do equal work and their trajectories have the same number of points. The grid is set per-symmetry in each cases/*.yml (chis: {trivial: [...], u1: [...]}, or a flat list applied to all): the block-sparse U(1) runs extend to much higher χ (up to 2048) than the dense/trivial runs (up to 768), since symmetry makes large χ far cheaper. TenPy is pinned with min_sweeps = max_sweeps, ITensor and YASTN run a fixed sweep loop, and MPSKit is driven one sweep at a time (its DMRG2 has a truncation-aware early stop that a single fixed-sweep call would trigger once converged; stepping it sweep-by-sweep forces the full count).
• Clean timing: the per-sweep hook stamps the wall-time boundary before measuring energy and subtracts the measurement cost, so neither the per-sweep times nor the total ever charge for it.
• Energy is a correctness gate, not the headline: at large χ all arms should agree on the energy. Wall time is the primary metric; matvec count (numops), where available, is the secondary metric.
• Reproducibility: each library env pins package versions; every results JSON records the resolved versions, hostname, thread count, and BLAS backend in its meta block.

BLAS-backend caveat

MPSKit and ITensorMPS route dense linear algebra through MKL; TenPy's and YASTN's numpy link OpenBLAS (no MKL-numpy is available without conda). All arms are pinned to the same thread count (BENCH_BLAS_THREADS), so the thread budget is equal, but the BLAS backend differs for the Python libraries. Absolute Julia-vs-Python wall times therefore carry a backend difference; the within-library default-vs-fixed and the numops comparisons are unaffected.

Layout

mps-benchmarks/
├── cases/                     # library-agnostic case definitions (YAML)
├── libraries/
│   ├── mpskit/  {Project.toml, common.jl, run.jl}   # isolated Julia env
│   ├── itensor/ {Project.toml, common.jl, run.jl}   # isolated Julia env
│   ├── tenpy/   {requirements.txt, common.py, run.py}
│   └── yastn/   {requirements.txt, common.py, run.py}
├── results/                   # committed raw JSON (one per library × case × symmetry)
├── plots/plot.py              # matplotlib aggregator → plots/figures/*.png (numpy + matplotlib only)
├── cluster/                   # optional Flatiron templates (run_case.sh, all_cases.disbatch, run_all.sbatch)
├── scripts/setup_envs.sh      # instantiate all four environments
└── bench.py                   # portable orchestrator (no Slurm needed)

Running

1. Set up the environments (once)

bash scripts/setup_envs.sh

This instantiates the two Julia projects and creates the two Python venvs from their requirements.txt. Julia resolves under 1.12.6 (via juliaup); override the binary with JULIA_BIN. YASTN is not on PyPI and is installed from its git repository (see libraries/yastn/requirements.txt).

2. Run the benchmark (portable, no Slurm)

# one case, one symmetry, all four libraries (sequential, fair timing)
python bench.py --case canonical_s1 --symmetry u1

# everything
python bench.py --case all --symmetry all

# a subset of libraries
python bench.py --case hubbard --libraries mpskit,tenpy

bench.py runs under a bare system Python (it shells out to each library's own environment). Environment overrides forwarded to every runner: BENCH_CHIS, BENCH_L, BENCH_NSWEEPS, BENCH_SEED, BENCH_REPS, BENCH_BLAS_THREADS.

3. Plot

python plots/plot.py canonical_s1 u1     # or:  python plots/plot.py all

Requires only numpy + matplotlib (present in the TenPy venv: libraries/tenpy/venv/bin/python plots/plot.py …).

On the cluster (Flatiron)

Submit from the repo root (the sbatch files are provided; you run sbatch):

sbatch cluster/run_all.sbatch                    # all units via disBatch in one allocation
# or one unit:
bash cluster/run_case.sh canonical_s1 u1         # inside an allocation

Each (case, symmetry) unit runs the four library runners in parallel with a matched per-runner thread budget (default 12 ⇒ 4 × 12 = 48 cores/unit). Runs are resumable: every runner writes results atomically after each χ point and skips completed points on restart.

Outputs

Per (case, symmetry) the plotter writes two figures into plots/figures/, split by eigensolver-arm family (where the interesting comparisons live):

• …__sweeptime_fixed.png — single-sweep wall time vs χ for the fixed ITensor-like arms (mpskit_fixed, itensor, tenpy_fixed, yastn). All do the same ncv=3 Krylov work, so this isolates raw per-sweep cost (implementation + BLAS backend), not eigensolver strategy.
• …__convergence_default.png — energy error vs cumulative wall time, per sweep, faceted by χ, for the default arms (mpskit_default, itensor, tenpy_default, yastn) — each library's out-of-the-box time-to-accuracy.

ITensor and YASTN appear in both families (their single arm is both the library default and the fixed ncv=3 reference). The plotter also prints a text summary (walltime, E − E_ref, t/sweep, Σ numops per arm/χ) covering all arms.

Adding a case

Drop a new cases/<name>.yml (copy an existing one), extend build_case in each libraries/*/common.{jl,py} if it introduces a new model, and regenerate cluster/all_cases.disbatch (command in that file's header). The repo name is deliberately algorithm-neutral so other algorithms (TDVP, VUMPS, excitations) can be added as sibling case/runner families later.