Exact n-site measurement with opt_einsum#

yastn.tn.fpeps.EnvCTM.measure_nsite_exact() contracts the CTM environment around the rectangular window spanned by the requested sites exactly, row by row. yastn.tn.fpeps.EnvCTM.measure_nsite_exact_oe() computes the same number by handing the whole window, corners, edges, ket, bra and operators, to opt_einsum as one tensor network. The contraction path is optimized, individual bonds can be unrolled (sliced) to bound peak memory, and the norm and the numerator can be evaluated separately so that many operators on the same window share one norm contraction.

EnvCTM.measure_nsite_exact_oe(*operators, sites=None, unroll=None, checkpoint_loop=False, optimizer='default', devices=None, mp_workers_per_device=0, projectors=None, per_combo_path=False, combo_path_kwargs=None) → float#

Memory-efficient version of measure_nsite_exact() using opt_einsum contraction path optimization, optional block-sparse index unrolling, and checkpointing.

For DoublePepsTensor PEPS, ket, operator and bra of every site enter the network as separate tensors, with the fermionic crossings between them as ncon swap pairs; edge middle legs are unfused to match.

For single-layer PEPS, falls back to the fused double-layer approach.

Returns <psi| O0_s0 ... |psi> / <psi|psi>. See measure_nsite_norm_exact_oe() for the norm-only contraction and measure_nsite_numerator_exact_oe() for the unnormalized numerator – callers that share a single norm across multiple numerator evaluations should use those split functions to control the autograd graph lifetime explicitly.

Parameters:

  • operators (Sequence[yastn.Tensor] | yastn.tn.mps.MpsMpoOBC) – Local operators of <O0_s0 O1_s1 ...>, one plain two-leg operator per site. Alternatively a single MPO with one tensor per site, e.g. from yastn.tn.mps.generate_mpo(), which measures a sum of operator products in one contraction. Its chain is sites as listed, in any order and not necessarily adjacent; it applies no reordering sign of its own, the sign of each term being part of the MPO. See Operators as MPO tensors.

  • sites (Sequence[tuple[int, int]]) – A list of sites [s0, s1, …] matching corresponding operators.

  • unroll (dict or None) – Dict mapping bond labels to int (uniform slice size) or list[SlicedLeg]. See Network layout, bond labels and unrolling for the bond-label scheme.

  • checkpoint_loop (bool) – If True and unroll is not None, each unroll iteration is wrapped in torch.utils.checkpoint.checkpoint(), trading recomputation for lower peak memory.

  • optimizer (str or opt_einsum.paths.PathOptimizer) – Contraction-path optimizer passed to opt_einsum.contract_path(). "default" (also None, "dp", "dynamic-programming") uses opt_einsum.DynamicProgramming(minimize="write", search_outer=False, cost_cap=True), which minimizes the size of the intermediates. Any other value accepted by opt_einsum, e.g. "greedy" or "auto", is passed through unchanged.

  • devices (Sequence[str] or None) – Devices to spread the contraction over, e.g. ["cuda:0", "cuda:1"]. None (default) contracts on the device of the PEPS tensors. With unroll, the slice combinations are dispatched across the devices, which needs mp_workers_per_device >= 1; a single device with one worker contracts serially there and moves the result back. Without unroll, the network is moved to devices[0] and contracted there.

  • mp_workers_per_device (int) – Number of worker processes per device in the multiprocess pool that contracts the slice combinations. 0 (default) disables multiprocessing and requires devices to be None or the PEPS device.

  • per_combo_path (bool) – Only used with unroll. If True, search a separate contraction path for every slice combination, tuned to its slice dimensions and cached by shape, instead of reusing one path for all combinations. Default False.

  • combo_path_kwargs (dict or None) – Options of the per-combination path search when per_combo_path=True; keys may include optimizer, memory_limit, names and who. None (default) means {"optimizer": optimizer}.

EnvCTM.measure_nsite_norm_exact_oe(*, sites, unroll=None, checkpoint_loop=False, optimizer='default', devices=None, mp_workers_per_device=0, projectors=None, per_combo_path=False, combo_path_kwargs=None)#

Contract only the norm <psi|psi> over the bounding window of sites.

Same contraction backend and options as measure_nsite_exact_oe(), with operators omitted; operator-bond labels ('opb', k) in unroll are ignored. Use this when sharing a single norm value across multiple numerator evaluations (see measure_nsite_numerator_exact_oe()).

EnvCTM.measure_nsite_numerator_exact_oe(*operators, sites, unroll=None, checkpoint_loop=False, optimizer='default', devices=None, mp_workers_per_device=0, projectors=None, per_combo_path=False, combo_path_kwargs=None)#

Contract only the unnormalized numerator sign * <psi| O0_s0 ... |psi>.

Same contraction backend and options as measure_nsite_exact_oe(); the result is not divided by the norm. The caller is responsible for dividing by <psi|psi> (typically obtained via measure_nsite_norm_exact_oe()).

Network layout, bond labels and unrolling#

The contraction builds a tensor network over the Nx x Ny window enclosing the requested sites: the four CTM corners TL, TR, BL, BR, the edge tensors T[j], B[j], L[i], R[i], and one double-layer site * per window position. Every edge of that network carries a tuple label; the unroll dict refers to those labels. Coordinates i (row, 0 ... Nx-1) and j (column, 0 ... Ny-1) are window-local, not absolute lattice positions, and follow the Site(x, y) = (row, col) convention:

       j=-1            j=0             j=1          j=Ny-1     j=Ny
        :              :               :             :          :
i=-1   TL --h,-1,-1-- T[0] --h,-1,0-- T[1] -- ... -- h,-1,Ny-1 -- TR
        |              |               |             |          |
     v,0,-1         v,0,0           v,0,1        v,0,Ny-1     v,0,Ny
        |              |               |             |          |
i=0    L[0]-h,0,-1-----*---h,0,0-------*--- ... --h,0,Ny-1-----R[0]
        |              |               |             |          |
     v,1,-1         v,1,0           v,1,1        v,1,Ny-1     v,1,Ny
        |              |               |             |          |
i=1    L[1]-h,1,-1-----*---h,1,0-------*--- ... --h,1,Ny-1-----R[1]
        :              :               :             :          :
        |              |               |             |          |
     v,Nx,-1        v,Nx,0          v,Nx,1       v,Nx,Ny-1    v,Nx,Ny
        |              |               |             |          |
i=Nx   BL --h,Nx,-1-- B[0] --h,Nx,0-- B[1] -- ... -- h,Nx,Ny-1 -- BR

  • Horizontal bonds ('h', i, j) run left to right between columns j and j+1 at row i. Rows i = -1 and i = Nx are the boundary rows (chi bonds between edge tensors and corners); i = 0 ... Nx-1 are PEPS rows; j = -1 is the left-boundary column and j = Ny-1 the right-boundary column for the chi bonds attached to the side edges.

  • Vertical bonds ('v', i, j) run top to bottom in column j between rows i-1 and i: i = 0 connects the top row to the first PEPS row, i = Nx the last PEPS row to the bottom row, i = 1 ... Nx-1 are interior; j = -1 and j = Ny are the left and right boundary columns.

  • Ket / bra split. For DoublePepsTensor PEPS every PEPS-row bond, i.e. ('h', i, j) with 0 <= i < Nx and all ('v', i, j), carries two labels, (*, 'k') for the ket layer and (*, 'b') for the bra layer. An un-qualified label in unroll is expanded into both; a layer-qualified label like ('v', 1, 1, 'k') slices only the ket side. Boundary (chi) bonds are single-label.

  • Operator bonds ('opb', k) label the bond between MPO tensors k-1 and k (next section); they exist only in the numerator network and are ignored by the norm.

Examples for a 2 x 3 window (Nx = 2, Ny = 3):

# unroll the horizontal bond between columns 0 and 1 at the first PEPS
# row, one index at a time (an int is a uniform slice size):
unroll = {('h', 0, 0): 1}

# unroll the vertical bond in column 1 between rows 0 and 1:
unroll = {('v', 1, 1): 1}

# several bonds at once; a left-boundary (chi) bond:
unroll = {('h', 0, 0): 1, ('v', 1, 1): 1}
unroll = {('v', 0, -1): 1}

# one charge sector at a time:
from yastn.tensor.oe_blocksparse import make_sliced_legs
unroll = {('h', 0, 0): make_sliced_legs(leg)}

The ket, the operator and the bra of every site stay separate network tensors (_build_ketbra_separate()), with the fermionic crossings between them as swap pairs of the contraction.

Plain charged operators#

A product of plain operators \(O_{s_1} O_{s_2} \cdots O_{s_n}\), listed in any order, is measured as an MPO of bond dimension one if any of the operators is charged (a single \(c\) or \(c^\dagger\)): yastn.tn.mps.product_mpo() of the operators. The bond between two neighbouring operators of the chain carries the total charge of the operators after it, and the charge travels along the MPO bonds exactly as the bonds of an MPO passed by the caller do (next section). Operators listed on the same site are multiplied, with the sign of bringing them together. A product of operators of zero charge needs no bonds; they sit on their sites as they are.

Operators as MPO tensors#

Building the MPO#

A sum of operator products on one set of sites,

\[O = \sum_t c_t \, o^{(t)}_{s_0} \, o^{(t)}_{s_1} \cdots o^{(t)}_{s_{L-1}},\]

can be measured in a single contraction instead of one contraction per term. Build it with yastn.tn.mps.generate_mpo() and pass the MPO in place of the plain operators:

import yastn.tn.mps as mps

H = mps.generate_mpo(mps.product_mpo(I, N=len(sites)),
                     [mps.Hterm(coeff, list(range(len(sites))), ops) for coeff, ops in terms])
value = env.measure_nsite_exact_oe(H, sites=sites)

The chain of the MPO is sites as listed: position k of an yastn.tn.mps.Hterm is the operator acting on sites[k], and I is the local identity, which fills the positions a term does not list. The sites may be listed in any order, and neighbouring positions of the chain need not be neighbouring sites of the lattice.

yastn.tn.mps.generate_mpo() returns the operator in the Fock basis: the entries of its tensors are the matrix elements of \(O\), the string of the chain among them. On the chain it is then an ordinary MPO, bonds running straight from one tensor to the next and crossing nothing:

    i0             i1             i2             i3
     |              |              |              |
   +-+--+         +-+--+         +-+--+         +-+--+
   | M0 |----1----| M1 |----2----| M2 |----3----| M3 |
   +-+--+         +-+--+         +-+--+         +-+--+
     |              |              |              |
    o0             o1             o2             o3

matrix element  <o0 o1 o2 o3| O |i0 i1 i2 i3>
1, 2, 3 = network labels ('opb', 1), ('opb', 2), ('opb', 3)

In the window each bond joins its two MPO tensors directly. Its swap gates are the ones the same operator would pick up if it were applied to the ket as a gate, the way yastn.tn.fpeps.Peps.apply_gate_() applies an MPO in time evolution: along a path of nearest-neighbour sites, each MPO tensor is contracted with the ket of its site and each bond becomes part of the ket’s virtual leg towards the next site of the path, with an identity on every site the path only passes. The measurement keeps the MPO tensors separate from the ket instead, and gives each bond the swap gates with the legs it would cross in that application.

What the measurement requires#

  • Either plain operators, one per listed site, or one MPO, passed as a single yastn.tn.mps.MpsMpoOBC; the two kinds are never mixed in one call.

  • Plain operators may list their sites in any order, and a site more than once. The product is taken as written, \(O_{s_1} O_{s_2} \cdots O_{s_n}\) with the rightmost operator acting first, whatever the lattice’s fermionic order.

  • The chain of an MPO is sites as listed: H[k] acts on sites[k]. Any order works, neighbouring positions of the chain need not be neighbouring sites, and each site appears once.

  • The measurement applies no sign of its own to an MPO: generate_mpo puts the sign of every term into the MPO, from the order in which that term lists its operators, with the same convention as plain operators.

A term may list its sites in any order, independently of the chain and of the other terms: each operator takes the position of its site in the chain, and generate_mpo sorts the term by those positions and multiplies it by the sign of that commutation:

H = mps.generate_mpo(mps.product_mpo(I, N=len(sites)),
                     [mps.Hterm(coeff, [sites.index(s) for s in term_sites], term_ops)
                      for coeff, term_sites, term_ops in terms])
value = env.measure_nsite_exact_oe(H, sites=sites)

Operator bonds and unrolling#

The bond between MPO tensors k-1 and k is labelled ('opb', k). It is unrolled like any other bond, e.g. unroll={('opb', 1): 2}, or one charge sector at a time with yastn.tensor.oe_blocksparse.make_sliced_legs(). The dimension-one bonds at the two ends of the chain are dropped and have no label. A product of plain charged operators has such bonds too, of dimension one, along its sites in the order listed. The norm network contains no operator bonds, so these entries are dropped when the norm is contracted.

Example#

Hopping in both directions plus density-density interaction on a diagonal pair of sites of a fermionic PEPS with CTM environment env. The sites shared by all terms are listed in sites, in any order. Every term names the site of each of its operators; to pass the terms to yastn.tn.mps.generate_mpo(), those sites are turned into positions in sites with sites.index:

import yastn
import yastn.tn.mps as mps

ops = yastn.operators.SpinlessFermions(sym='U1')
c, cp, n, I = ops.c(), ops.cp(), ops.n(), ops.I()
sites = [(1, 0), (0, 1)]
terms = [(-1.0, [(0, 1), (1, 0)], [cp, c]),   # c+_(0,1) c_(1,0)
         (-1.0, [(1, 0), (0, 1)], [cp, c]),   # c+_(1,0) c_(0,1)
         (0.5, [(0, 1), (1, 0)], [n, n])]
H = mps.generate_mpo(mps.product_mpo(I, N=len(sites)),
                     [mps.Hterm(co, [sites.index(s) for s in ss], oo) for co, ss, oo in terms])

norm = env.measure_nsite_norm_exact_oe(sites=sites)
num = env.measure_nsite_numerator_exact_oe(H, sites=sites, unroll={('opb', 1): 1})
value = num / norm

# the same as the plain measurements, each term with its sites as it lists them
ref = sum(co * env.measure_nsite_numerator_exact_oe(*oo, sites=ss) for co, ss, oo in terms) / norm

Listing sites in another order, with the same terms, gives the same value:

sites_r = [(0, 1), (1, 0)]
H_r = mps.generate_mpo(mps.product_mpo(I, N=len(sites_r)),
                       [mps.Hterm(co, [sites_r.index(s) for s in ss], oo) for co, ss, oo in terms])
value_r = env.measure_nsite_numerator_exact_oe(H_r, sites=sites_r) / norm   # == value