Skip to content

DMRG

dmrg_ sweeps a state to the ground state of an operator, in place, with two-site updates. Everything on this page runs on the core install: no scipy, no quimb, no jax.

>>> from tenet.models import spin_half
>>> from tenet.network import MPO, MPS, dmrg_
>>> from tenet.symmetry import U1Sector
>>> site, n = spin_half(), 8
>>> terms = []
>>> for i in range(n - 1):
...     terms.append((1.0, [(site.ops["Sz"], i), (site.ops["Sz"], i + 1)]))
...     terms.append((0.5, [(site.ops["S+"], i), (site.ops["S-"], i + 1)]))
...     terms.append((0.5, [(site.ops["S-"], i), (site.ops["S+"], i + 1)]))
>>> h = MPO.from_terms(n, terms)
>>> psi = MPS.product(site.phys, [U1Sector(1 if i % 2 else -1) for i in range(n)])

Three terms per bond, a Néel product state as the seed, and no bond space written down anywhere: the MPO's grading comes out of the operators' own charges, and the state's comes out of the seed's.

>>> out = dmrg_(psi, h, chi=32)
>>> round(out.energy, 9)
-3.374932599

One call sweeps to convergence and mutates psi in place. Everything below is a lens on that call.

What is being minimized

DMRG minimizes the Rayleigh quotient over matrix product states of a bounded bond dimension:

\[ E[\psi] \;=\; \frac{\langle\psi\vert H\vert\psi\rangle}{\langle\psi\vert\psi\rangle}, \qquad \Psi_{s_1\cdots s_N} \;=\; \sum_{\{\alpha\}} A^{s_1}_{\alpha_0\alpha_1} A^{s_2}_{\alpha_1\alpha_2}\cdots A^{s_N}_{\alpha_{N-1}\alpha_N}, \]

with \(\alpha_0\) and \(\alpha_N\) one-dimensional — that is the D=1 boundary leg every MPS carries. Site \(n\)'s tensor is psi[n], whose legs are (left bond OUT, physical OUT, right bond IN), so \(\alpha_{n-1}\), \(s_n\), \(\alpha_n\) are axes 0, 1, 2 in that order. Charge flows left to right, \(\alpha_{n-1} \otimes s_n \to \alpha_n\), which is why a non-unit sector on bond 0 is the targeted total charge.

The minimization is not global. A sweep fixes every tensor but two, which makes \(E\) a Rayleigh quotient in the remaining variables alone: with the rest of the chain in canonical form the norm denominator is the identity, and the stationarity condition is an ordinary Hermitian eigenproblem for the two-site block \(\Theta_{\alpha_{n-1} s_n s_{n+1} \alpha_{n+1}}\),

\[ H^{\mathrm{eff}}_{(n,n+1)}\,\Theta \;=\; E\,\Theta, \]

where \(H^{\mathrm{eff}}\) is the left environment, the two MPO tensors \(W_n\), \(W_{n+1}\) and the right environment, contracted around the open block. Env.heff2 is that operator as a matrix-vector product — it is never formed as a matrix — and lanczos finds its lowest eigenpair. The updated \(\Theta\) is then split back into two site tensors by a truncated SVD, and that split is where the bond dimension bound enters: Truncation is what it means.

Two consequences worth keeping in view. The energy is variational, so a converged out.energy is an upper bound on the true ground-state energy. And \(H^{\mathrm{eff}}\) is built from the same graded tensors as everything else, so it is block diagonal in the symmetry sector: the eigensolver cannot leave the sector the boundary legs fixed, and dmrg_ needs no constraint machinery to keep it there.

The driver

dmrg_ right-canonicalizes psi first, so a freshly seeded random MPS or a product state is the expected input. It mutates psi and returns a DMRG_out whose psi is that same object:

field what it holds
psi the converged state — the object you passed in
sweeps how many sweeps ran
energy, denergy the last sweep's energy and energy change
max_dSchmidt the last sweep's worst-cut Schmidt change
max_discarded_weight the last sweep's maximum per-bond discarded weight
history one (energy, denergy, dSchmidt, discarded) tuple per sweep
schedule the realized schedule, one Sweep per sweep run

zip(out.schedule, out.history) is exact, so out.schedule alone answers whether a run reached its final chi or converged earlier.

>>> len(out.history) == out.sweeps == len(out.schedule)
True

Targeting a sector

The state's boundary legs fix the symmetry sector, and the site tensors' invariance keeps it there. MPS.product takes one physical sector per site and derives the bonds backwards from those charges, so the total lands on bond 0 where it is printable:

>>> neel = MPS.product(site.phys, [U1Sector(1), U1Sector(-1)] * 4)
>>> neel[0].legs[0].space.sectors          # S^z_tot = 0
((U1Sector(charge=0), 1),)
>>> up = MPS.product(site.phys, [U1Sector(1)] * 4)
>>> up[0].legs[0].space.sectors            # S^z_tot = 2
((U1Sector(charge=4), 1),)

The general recipe is a D=1 boundary leg carrying U1Sector(q), targeting \(S^z_{\mathrm{tot}} = q/2\); the boundary charge and the tensors' invariance hold it. MPS.product is Abelian-only: a single sector is not a non-Abelian multiplet, so under SU(2) the route is MPS.random with a charged boundary leg.

Schedules

dmrg_'s flat chi/cutoff keywords are one spelling of a one-entry schedule. The other is schedule=, a list of Sweep entries — one record per sweep, whose last entry repeats until convergence or max_sweeps:

from tenet.network import Sweep, dmrg_

dmrg_(psi, h, chi=64)                              # flat
dmrg_(psi, h, schedule=[Sweep(chi=64)])            # the same run, exactly
dmrg_(psi, h, schedule=[Sweep(32, noise=1e-4)] * 4
              + [Sweep(64, noise=1e-5)] * 4
              + [Sweep(64)], max_sweeps=20)        # ramp, cool down, converge

A Sweep carries chi, cutoff, noise and noise_type. The loop tolerances — energy_tol, schmidt_tol, max_sweeps, ncv — are properties of the loop, not of a sweep, and stay flat keywords on dmrg_.

Passing schedule together with chi or cutoff raises, and so does an empty schedule.

A callback= is invoked once per sweep with that sweep's DMRG_out, after history is appended, so it sees the sweep that just finished. Its return value is ignored: there is no early-stop protocol.

Convergence

The loop stops when both criteria are met in one sweep: the energy change is below energy_tol (default 1e-12) and the worst-cut Schmidt change is below schmidt_tol (default 1e-8). The Schmidt criterion is the sensitive one, and it is what catches a run whose energy has plateaued on a wrong bond structure.

Convergence is never declared on a sweep that is still noisy or still inside the schedule: the exit requires the schedule's last entry and noise == 0.0. An energy that stopped moving under noise at a ramp's intermediate chi has converged to the wrong thing.

Converged, or only plateaued?

Both criteria are change tests, and a change test is satisfied by a run stuck on a wrong bond structure: nothing moved because nothing could. The check that is not a change test is the energy variance:

>>> round(abs(h.variance(out.psi)), 9)
0.0

\(\langle\psi\vert H^2 \vert\psi\rangle/\langle\psi\vert\psi\rangle - E^2\) is zero for an exact eigenstate. Run it at two bond dimensions: a state converging on an eigenstate has a variance falling towards zero as chi grows; a state plateaued on the wrong structure has one that does not.

Extrapolating in the discarded weight

The truncation error is linear in the discarded weight near convergence, so a fit of energy against discarded weight extrapolates to the \(\chi \to \infty\) energy. Fit a reverse schedule run on the already-converged state, not the forward run's history — the forward energies are not converged at their bond dimensions:

out = dmrg_(psi, h, chi=64)                                  # converge first
rev = dmrg_(psi, h, schedule=[Sweep(16)] * 2 + [Sweep(12)] * 2 + [Sweep(8)] * 2,
            energy_tol=0.0, max_sweeps=6)                    # nothing exits early
dws      = [record[3] for record in rev.history]
energies = [record[0] for record in rev.history]
e_extrapolated = np.polyfit(dws, energies, 1)[1]

The conventions worth keeping: zero noise throughout the reverse run; energy_tol=0.0 so no sweep exits early; an even number of sweeps per chi, since odd and even half-sweeps report different discarded weights; the first reverse chi slightly below the last forward one; and the error bar quoted as one fifth of the extrapolation distance. The recipe needs two-site sweeps, which is what dmrg_ runs.

Noise

noise on a Sweep mixes a perturbation in at each split, at relative strength noise. noise_type says which perturbation, and therefore which split runs:

(noise, noise_type) the split
noise == 0.0, any noise_type svd_truncated of the two-site tensor
> 0, "wavefunction" svd_truncated of a perturbed two-site tensor
> 0, "perturbative" eigh of a perturbed density matrix

Nothing else decides it — no bond width, no chi, no runtime probe. noise=0.0 draws no random number and builds no density matrix.

Wavefunction noise (the default) adds a random symmetric tensor over the two-site tensor's own legs after the eigensolver and before the split, then renormalizes. It fills every structurally allowed coupled sector of the two-site map, including the ones the eigensolver left numerically empty and which the truncation therefore dropped from the bond. That is the local minimum a symmetric DMRG falls into: a sector that is zero stays zero otherwise. It cannot reach outside \(\mathrm{bond}_l \otimes \mathrm{phys}\).

Perturbative noise builds the density matrix and splits with eigh instead. Squaring the two-site tensor into rho resolves a singular value \(\sigma\) through \(\sigma^2\), so the split's accuracy floor is the square root of machine epsilon, which is why a noiseless sweep — including the cooling tail of a ramp — takes the SVD split.

Noise is not variational: a noisy sweep's energy can sit above its clean twin. Taper it — 1e-4 early, 1e-5 in the middle, 0.0 at the end.

seed= makes the draw at bond n reproducible as seed + n, distinctly per sweep.

Excited states

orthogonal_to= takes already-converged states and holds psi orthogonal to them, turning the run into an excited-state search:

ground = dmrg_(psi1, h, chi=64)
first  = dmrg_(psi2, h, chi=64, orthogonal_to=[ground.psi])

The given states are not modified and need no particular gauge. The machinery is one two-state Env per given state over MPO.identity, swept alongside the main environment; what each contributes at a bond is a projection vector handed to the eigensolver. The reported energy is the projected operator's own Ritz value, so it is the excited energy directly, with no shift to subtract.

Sector targeting composes with it: a charged D=1 boundary leg fixes the sector, orthogonality walks up the spectrum inside it, and a converged state whose boundary legs put it in a different sector is dropped from the projection — the symmetry already made it orthogonal.

Restarting

Save the state, load it back, and re-enter the schedule at a slice:

out = dmrg_(psi, h, schedule=schedule[:2])   # ... interrupted after two sweeps
psi = MPS.load("checkpoint")
out = dmrg_(psi, h, schedule=schedule[2:])   # matches the uninterrupted run

The slice is the position, and a sweep is a full round trip, so no direction or sweep-index argument is needed. Saving and loading covers the files.

Measuring the converged state

A measurement is an expectation on out.psi, and it is a different act from reading out.energy. out.energy is the effective eigenvalue the two-site eigensolver returned at the orthogonality centre of the final sweep, that is \(E = \langle \psi \vert H \vert \psi \rangle\); a measurement contracts an operator you choose against the state afterwards.

The bond energy is the clearest case. The Hamiltonian is a sum of two-site terms, one per bond of the open chain,

\[ H = \sum_{n=1}^{N-1} h_{n,n+1}, \qquad h_{n,n+1} = \mathbf{S}_n \cdot \mathbf{S}_{n+1}, \]

and the bond energy is the expectation of one such term, \(e_n = \langle \psi \vert h_{n,n+1} \vert \psi \rangle\). By linearity \(\sum_n e_n = E\), which makes the sum a cross-check that the two readings are on one scale:

>>> from tenet.network import expectation_2site, local_op
>>> ss = local_op(site.matrices["S.S"], phys=site.phys)
>>> bonds = [expectation_2site(out.psi, ss, i) for i in range(n - 1)]
>>> round(sum(bonds), 9) == round(out.energy, 9)
True
>>> [round(e, 4) for e in bonds]
[-0.6612, -0.2843, -0.5874, -0.3092, -0.5874, -0.2843, -0.6612]

The canonical form of the state enters only as the evaluation method, not as part of the definition: with the orthogonality centre on bond \(n\) every environment outside the bond is an identity, so \(e_n\) is a contraction two sites wide rather than a pass over the whole chain. The profile is not flat — open boundaries induce an alternating, dimerized pattern, strongest at the two edges and decaying inward towards the uniform bulk value, which in the thermodynamic limit is \(1/4 - \ln 2 \approx -0.4431\) per bond.

The rest of the measurement set reads the same way:

>>> import numpy as np
>>> import tenet
>>> from tenet import IN, OUT, Leg
>>> from tenet.network import expectation_profile, expectation_1site, overlap
>>> sz = tenet.SymmetricTensor.from_dense(
...     np.diag([-0.5, 0.5]), (Leg(site.phys, OUT), Leg(site.phys, IN))
... )
>>> profile = expectation_profile(out.psi, sz)
>>> max(abs(v) for v in profile) < 1e-9
True
call what it returns
expectation_1site \(\langle\psi\vert o_n \vert\psi\rangle/\langle\psi\vert\psi\rangle\) at one site
expectation_2site the same for a rank-4 operator on (n, n+1)
expectation_profile \(\langle\psi\vert o_n \vert\psi\rangle/\langle\psi\vert\psi\rangle\) at every site, in one pass
overlap \(\langle \phi \vert\psi\rangle\), undivided
measure_mpo \(\langle\phi\vert H \vert\psi\rangle\), undivided
correlation_function {(i, j): value} for the pairs you ask for

Three things to know:

  • expectation_profile is the one to reach for over a list comprehension. A per-site loop costs two full-chain transfer passes per site; the profile moves the orthogonality centre once along the chain and reads the operator off it. Same numbers, \(O(N)\) instead of \(O(N^2)\).
  • overlap and measure_mpo do not divide. A fidelity is overlap(phi, psi) / (phi.norm() * psi.norm()), and you write the division. The divided readings are the ones whose names say expectation.
  • correlation_function takes the rank-3 charged operators MPO.from_terms takes, which is the form a fermionic c has. It returns {(i, j): value} for i < j, every pair by default, and the Jordan-Wigner string between i and j is the fZ2 braiding the term builder inserts. It costs one MPO build and one pass per pair, so pass pairs= for the row or the distance you want.

Env(psi, h, bra=phi) is the object all of this stands on, and measure_mpo(phi, h, psi) is its one-line spelling. Env.heff2 refuses on a two-state environment — the prepared matvec reads the IdL/IdR channels as gauge identities, true of a canonical chain against itself and false of a mixed transfer — which is why measurement and the sweep are different entry points into one cache. Env.measure returns the unnormalized \(\langle\psi\vert H \vert\psi\rangle\).

The entanglement profile

>>> entropy = out.psi.entanglement_entropy()          # {bond: S}, keyed by the bond's left site
>>> sorted(entropy) == list(range(n - 1))
True
>>> renyi2 = out.psi.entanglement_entropy(alpha=2)    # the Renyi family, same keys
>>> values = out.psi.schmidt_values()                 # the spectrum the entropy comes from
>>> sectors = out.psi.schmidt_sectors()[3]            # bond 3's spectrum, by symmetry sector
  • The unit is nats. \(S = (c/6) \log x\) on an open chain wants the natural log; divide by \(\log 2\) for bits.
  • The key is the bond's left site, \(0 \dots N-2\) — the same key the sweep's schmidt dict uses. The two trivial boundary cuts are zero and are not returned.
  • schmidt_sectors is the read a graded bond is for. On an SU(2) bond a single \(j\) multiplet stands for \(2j + 1\) dense Schmidt values, and the entropy accounts for that, which is why a two-site singlet reports \(\log 2\) under SU(2) and under U(1) alike.

Each of the three readers canonizes a copy of the state and runs its own SVD sweep, so they never re-gauge the state you hand them. Keep the result rather than calling twice on a large state.

DMRG_out carries no spectrum field: the spectrum is a property of the state, and out.psi answers for it exactly and in any gauge.

Compressing

>>> discarded = out.psi.compress_(chi=8)
>>> out.psi[4].legs[0].space.dim <= 8
True

compress_ truncates in place and returns the total discarded weight, \(\sqrt{\sum_{\mathrm{bonds}} \mathrm{dw}}\), where a sweep reports the per-bond maximum. One answers "how much of my state did I throw away", the other "which bond is the convergence diagnostic".

compile=

The two-site matvec is a pure function of a fixed contraction structure, so it can be handed to jax.jit. dmrg_ takes the callable; this layer names no accelerator, so you supply it and the jax extra:

import jax                       # pip install "symtenet[jax]"
import tenet

tenet.enable_jax()               # registers SymmetricTensor as a JAX pytree

dmrg_(psi, h, chi=64)                    # the plain run, NumPy
dmrg_(psi, h, chi=64, compile=jax.jit)   # the same run, matvec compiled

It changes the run's performance regime, not its accuracy. Read the shape of it before reaching for it: the sweep around the matvec is not traceable, because the truncating SVD re-decides each bond space every sweep, so on the JAX backend a compiled run is slower end to end than the plain NumPy one; and Env rebuilds its prepared operator at every bond visit and so invokes compile again each time, one XLA trace per bond per sweep.

The sweep's fixed choices

Every sweep is a two-site update, which is what the extrapolation recipe above consumes. The split is svd_truncated, or eigh of a density matrix under perturbative noise, per the table above. The contractions inside the sweep run in hand-written pairwise orders, which is what keeps them right on a U(1) bond with unevenly filled sectors, where a cost model reading physical leg sizes misprices the network.

The matvec the eigensolver calls is Env.heff2, and it has two paths. Which one runs is decided when the operator is built, not at run time:

the MPO was built heff2 runs
at the default the site-tensor contraction — the two W tensors between the two environments
symbolic=True the prepared matvec — complementary operators assembled per bond, the sum dispatched term family by term family

The first is the cheapest thing that can happen to a lattice model whose MPO bond is five or eight wide. The second is what a Hamiltonian with a bond in the thousands needs — in practice, quantum chemistry, where it is the route that fits in memory at all. Nothing probes a threshold and nothing switches mid-run: the representation the operator is in when you hand it to dmrg_ is the choice. MPO.materialize() moves an operator from the second path to the first, and Building a Hamiltonian is where the choice is made.

heff2 refuses outright on a two-state environment — see the measurement section above — because the prepared path reads the IdL/IdR channels as gauge identities, which is true of a canonical chain against itself and false of a mixed transfer.

Where next