Choosing an algorithm¶

A circuit is applied to the state in blocks: gates are fused into an operator, and the operator is applied in one call. algorithm selects how that application is carried out.

sim = tw.TwSimulator(bonddim=128, algorithm="vmpoa")

All the algorithms compute the same product. They represent it exactly when the state fits under bonddim, and they differ only in how much of it survives when it does not. So the first thing to establish is whether the ceiling binds at all: if it does not, pick on speed alone.

The names say what each one is. The letter is the method and mpo is what it applies:

Setting	Stands for
`"dmpo"`	direct MPO application
`"cmpo"`	canonical MPO application
`"vmpoa"`	variational MPO application, variant a
`"vmpob"`	variational MPO application, variant b

A second setting, strict_bonddim, decides whether the bonds may grow past bonddim while the operator is being applied. Each algorithm has its own default for it, and the defaults keep memory bounded by bonddim. See Strict and relaxed bond dimension below.

The four¶

"dmpo", direct MPO application, applies the operator in a single pass over the chain, the zip-up. It is the cheapest per application. By default it cuts its intermediate to bonddim as it goes, like the reference MPS engine's zip-up, so its memory stays where bonddim puts it.

"cmpo", canonical MPO application, also applies the operator in a single pass, arranged so that the weight it discards is chosen against the state's own Schmidt spectrum. The name is the distinction: the truncation happens in canonical form, which is what makes those Schmidt values the right measure. It costs about what "dmpo" costs. When the ceiling binds it keeps more of the state, and it tends to leave the state narrower afterwards, which makes the gates that follow cheaper as well. It is also the accurate choice on GPU.

"vmpoa", variational MPO application variant a and the default, refines the result with single-site variational sweeps. It is the most accurate at a given bonddim on operators with long-range structure, and it costs roughly twice a single pass. It is the default because accuracy at a fixed bonddim is usually worth more than throughput: if a cheaper algorithm forces you to raise bonddim to recover the same fidelity, it was not the cheaper option.

"vmpob", variant b, refines with two-site sweeps. It converges best of the four and does the most work per sweep. Unlike the single-site variant it can widen a bond during refinement, which matters when the initial guess is too narrow to represent the answer.

Strict and relaxed bond dimension¶

bonddim caps the bond dimension of the state you get back. While an operator is being applied, though, the algorithm works with intermediates that are wider than the result. strict_bonddim decides whether those intermediates are held to bonddim too.

sim = tw.TwSimulator(bonddim=64, algorithm="dmpo", strict_bonddim=False)

Strict (strict_bonddim=True) keeps every intermediate bond at or below bonddim. The zip-up cuts its carry to bonddim as it sweeps, so the memory a run needs follows from bonddim alone and you can size a machine by it. The price is where those cuts happen: in a gauge where the singular values are not the state's Schmidt values. A cut there can drop weight the exact result did not need to lose. The fidelity the run reports includes that loss, so a strict run is less accurate, not silently wrong.

Relaxed (strict_bonddim=False) lets the zip-up's carry grow to max(4 * bonddim, D), where D is the bond dimension of the operator, and then cuts to bonddim in canonical form, where the singular values are the right ones to compare. The result is exact whenever it fits under bonddim, whatever entdim is. It costs more memory and 1.5 to 4 times the time.

Either way the state you get back has no bond above bonddim.

The defaults¶

Left at None, each algorithm takes its own default:

`algorithm`	Default	Strict	Relaxed
`"dmpo"`	strict	the zip-up cuts its carry to `bonddim`	the carry grows to `max(4 * bonddim, D)`
`"vmpoa"`, `"vmpob"`	strict	the sweeps start from the strict zip-up	the sweeps start from the relaxed zip-up
`"cmpo"`	relaxed	not available: `ValueError`	its sketch grows to `2 * bonddim`

"cmpo" has no strict form. Its random sketch is up to twice bonddim wide, so its memory is still tied to bonddim, just not equal to it. Asking for algorithm="cmpo", strict_bonddim=True is refused when the simulator is built. A strict run also needs bonddim >= entdim; see Strict needs bonddim at least entdim.

For the variational algorithms the setting chooses the starting point of the sweeps. A single-site sweep ("vmpoa") cannot change a bond dimension, so it keeps the bond profile it starts from: on 20-qubit random circuits at bonddim 16 to 128, relaxed "vmpoa" reached 5 to 15% lower infidelity than strict, at 1.5 to 4 times the time. Two-site sweeps ("vmpob") adapt the bonds on their own, and relaxed "vmpob" matched strict to within 1% at 1.4 to 3 times the time.

Strict needs `bonddim` at least `entdim`¶

A fused operator can reach bond dimension entdim (see entdim). Once the operator is wider than bonddim, a strict cut can do more than lose accuracy: it can wipe out the state altogether, driving even a state that is exactly a product state to zero norm.

So a run that resolves to strict with bonddim < entdim is refused with a ValueError when the simulator is built, before any circuit runs. The check uses entdim rather than the width each operator actually reaches, so whether a run is accepted does not depend on the circuit. To run it anyway:

Raise bonddim to entdim or above. With the defaults, bonddim=256 and entdim=16, the operator never outgrows the state.
Lower entdim to bonddim or below. Fewer gates fuse into each operator, so more applications run, each on a narrower operator.
Pass strict_bonddim=False, or use "cmpo", which is relaxed by default. Relaxed accepts any entdim at the cost in memory and time above.

import tensorweaver as tw

tw.TwSimulator(bonddim=4, entdim=16, algorithm="dmpo")   # ValueError
tw.TwSimulator(bonddim=4, entdim=16, algorithm="dmpo",
               strict_bonddim=False)                     # exact when it fits
tw.TwSimulator(bonddim=16, entdim=16, algorithm="dmpo")  # strict, accepted

The low-level MPS.apply_mpo makes the same check against the bond dimension of the operator it is given, and raises BondDimensionError.

How much strict costs in accuracy¶

From entdim up, strict still cuts in the wrong gauge, and the cost depends on the algorithm:

"dmpo" suffers most. With bonddim only a little above entdim, the strict zip-up loses weight as a matter of routine. Deferring the compression sweep (post_compress=False) makes it worse.
The variational algorithms repair most of it. Their sweeps start from the strict zip-up and then optimise the result, so most of what the seed got wrong is recovered.

The fidelity a strict run reports includes what it lost. It can overestimate a fidelity that is already low, so read a low reported fidelity as "this run is not accurate" rather than as a precise figure. When it is too low, raise bonddim or pass strict_bonddim=False.

A high reported fidelity is still not proof of convergence. Check it as described in Deciding it by measurement.

Choosing between them¶

Choose strict (the default for "dmpo", "vmpoa", "vmpob") when memory has to be predictable from bonddim: a shared machine, a large bonddim close to what fits, or a batch of runs sized in advance. Watch the reported fidelity, which shows what it cost.
Choose relaxed when memory has headroom and accuracy at a given bonddim matters more than time: a small bonddim, an entdim above bonddim (where strict is refused), or a circuit whose state you expect to stay narrow and want exactly.
If you want relaxed behaviour at single-pass cost, "cmpo" already runs that way.

Picking one¶

Situation	Reach for
The state stays well below `bonddim`	`"cmpo"`, or `"dmpo"` with `strict_bonddim=False`. Both are exact here, and both are single passes.
A long run at a `bonddim` that binds	`"cmpo"`. Single-pass cost, and the narrower state it leaves compounds over the run.
Highest accuracy at a fixed `bonddim`	`"vmpoa"`, or `"vmpob"` where the operator has structure a single-site sweep converges poorly on.
Reproducing a MIMIQ or MPSSim baseline	`"vmpoa"`, which starts from the same zip-up as MPSSim's default.
A little more accuracy than `"vmpoa"`, at a higher cost	`"vmpoa"` with `strict_bonddim=False`.
`bonddim` below `entdim`	`"cmpo"`, or any algorithm with `strict_bonddim=False`. Strict is refused here.
`usegpu=True`	`"cmpo"`. See GPU engine.
You are not sure	Leave the default. Then sweep.

Deciding it by measurement¶

The reliable way to choose is to run the same circuit both ways and compare the observable you actually care about, not the reported fidelity.

import tensorweaver as tw

for algorithm in ("dmpo", "cmpo", "vmpoa"):
    sim = tw.TwSimulator(bonddim=64, algorithm=algorithm, seed=1)
    state, fidelity = sim.evolve(sim.zerostate(n), sim.compile(circuit))
    print(f"{algorithm:>6}  chi {state.mps.max_used_bond_dim:>4}  "
          f"fidelity {fidelity.value:.6f}")

Two columns are worth reading, not one.

The fidelity says how much weight each algorithm kept. A larger number means less was discarded.

The final bond dimension says how wide the state it produced is. This is the column that is easy to overlook and often the more useful of the two: an algorithm that finishes narrower is not only more compact, it makes every subsequent gate cheaper, and the effect compounds across a long circuit. Two algorithms can report similar fidelities and leave states of very different width.

If the two algorithms give the same answer, take the cheaper. If they disagree, the ceiling is binding and the disagreement is telling you the run is not converged. Raise bonddim until they agree, and then decide on cost.

On the GPU¶

A state moved onto the device runs "dmpo" and "cmpo" there. The variational algorithms run on the host, and the simulator refuses them rather than transferring the state back and forth for every application, which would cost more than the refinement saves. Between the two that do run, prefer "cmpo".

What the choice does not affect¶

algorithm and strict_bonddim do not change the ceiling, the cutoffs, or the layout. It changes which weight is discarded once something has to be, and nothing else. If a run is inaccurate because bonddim is far too small for the circuit, no algorithm will rescue it; see Cutoffs and bond dimension.