Concurrent sweeps and error control

The search sweeps variables in a fixed order tied to the contraction order, so rotating and reflecting the lattice starts the sweep from a different corner and materially stabilizes the result. The standard protocol therefore solves the same instance once per element of all_lattice_transformations and keeps the best outcome.

Those solves are independent. Since the contraction caches became contractor-owned there is nothing to stop them running at once, and sweep_transformations does exactly that.

potts_h = potts_hamiltonian(
    ising_graph(instance);
    spectrum = full_spectrum,
    cluster_assignment_rule = super_square_lattice((m, n, t)),
)
params = MpsParameters{Float64}(; bond_dim = 16, num_sweeps = 1)

sweep = sweep_transformations(
    transform -> MpsContractor(
        SVDTruncate,
        PEPSNetwork{KingSingleNode{GaugesEnergy},Dense,Float64}(m, n, potts_h, transform),
        params;
        onGPU = true, beta = 2.0, graduate_truncation = true,
    ),
    SearchParameters(; max_states = 2^8, cutoff_prob = 1e-4);
    merge_strategy = ctr -> merge_branches(ctr; merge_prob = :none),
)

sol = best_solution(sweep)

build_contractor is called once per transformation inside the task that will run it, so each solve gets its own network, projector workspace, and contraction cache. Julia must be started with more than one thread (julia -t auto) for the solves to overlap; with a single thread the sweep runs serially and says so.

How much did concurrency buy on the tested systems? Up to ~3× on CPU, up to ~1.7× on H100

The answer depends entirely on the device, so both halves are given below. Same protocol throughout: all eight transformations, interleaved A/B rounds with a full GC.gc (plus CUDA.reclaim() on the device) before every timed section, reported as the median of the per-round paired ratios against the serial loop.

On CPU it scales

Xeon Platinum 8462Y+, CUDA_VISIBLE_DEVICES="", 5 rounds:

case	serial	c=1	c=2	c=4	c=8
chimera 3×4×3, D=16, `SVDTruncate`	0.04 s	0.88×	1.41×	2.14×	3.09×
chimera 128power, D=32, `Zipper`	16.0 s	0.89×	1.11×	1.31×	1.64×

Monotonic in concurrency, so :auto uses min(n, Threads.nthreads()) here. c=1 is below 1.0 because the driver's fixed overhead (calibration, bookkeeping) is not amortized when there is nothing to overlap — it costs most on the 70 ms case.

On a GPU it depends on the card

On the tested RTX 5080 fanning out did not pay (7 rounds):

case	serial (median)	sweep, c=2	sweep, c=4
chimera 3×4×3, D=16, `SVDTruncate`	3.97 s	0.92×	0.80×
chimera 128power, D=32, `Zipper`	14.30 s	0.88×	0.89×

the solves overlap but per-solve time degrades at the same rate, so no level beats the serial loop. The device is not the constraint: utilization stays near 10%, and what saturates is serialization inside the CUDA API and allocator, which this solver provokes because its kernels are small and numerous. The tested H100 did have headroom to overlap them (7 rounds):

case	serial (median)	sweep, c=2	sweep, c=4
chimera 3×4×3, D=16, `SVDTruncate`	1.38 s	1.69×	1.44×
chimera 128power, D=32, `Zipper`	18.18 s	1.22×	1.43×

The gain begins at c=2; c=1 is a slight net loss because the driver's fixed overhead is not amortized when nothing overlaps.

Because the common case is the smaller card, concurrency = :auto stays at 1 on any GPU as a conservative default. Set concurrency = 2–4 explicitly on a large device, for a CPU-only run, for several devices, or for an instance you have measured yourself.

Measuring this is harder than it looks

A naive protocol — warm up, time the serial arm, then time the concurrent one — reported 1.68× and 1.06× for these same two cases. That was entirely an artifact: timing the serial arm straight after a concurrent warm-up leaves the CUDA pool in a state that inflated the baseline by up to 2.7×. Interleave the arms and reclaim before every timed section, and report paired per-round ratios rather than a min over separate batches.

Below ~2000 spins, this solver is faster on the CPU

Compare the serial columns above: 0.04 s versus 1.38 s, and 16.0 s versus 18.18 s. Solving identical configurations on each device across three instance sizes and four bond dimensions puts the crossover well up the size range (CPU/GPU wall-clock ratio; below 1 means the host is faster):

spins	D=8	D=16	D=32	D=64
36 (dense)	0.02	0.02	0.02	0.02
128 (dense)	0.14	0.21	0.36	0.46
2048 (dense)	0.50	0.64	0.90	—
2048 (sparse)	—	0.78	1.45	—

The device wins only in the largest sparse case (2048 spins, bond 32). On the tested Xeon/H100 system, the CPU is preferable for exploratory work in the smaller or dense configurations. MpsContractor defaults to onGPU = true, which suits the D-Wave-scale sparse regime the package targets; measure the intended system before assuming. Energies agree in every cell.

So: the sweep was a real speed-up on both tested devices — up to ~3× on CPU and 1.2–1.7× on the H100; on the tested RTX 5080 the throughput gain was marginal, which is why :auto stays at a conservative 1 on any GPU. Its other values hold regardless of device — one call instead of a hand-written loop, a device-memory budget that keeps concurrency from exhausting the card, deterministic per-transformation seeding (a Zipper sweep was not reproducible under concurrency before), failure isolation, and the cross-transformation agreement diagnostics below.

Why the governor counts bytes

The constraint on fanning out eight solves is not CPU cores, it is device memory. A single bond-32 solve on a Pegasus-scale instance can hold several GiB, so a fixed eight-way fan-out reliably exhausts a consumer card, while a small bond-16 Chimera instance would happily run all eight. The number of solves that fit is a property of the instance, not of the machine — so concurrency is gated by a byte budget:

One transformation is solved alone while its peak device usage is measured. This is not overhead: its solution counts toward the sweep, and it absorbs the compilation that eight concurrent solves would otherwise duplicate under contention.
That measured peak, times reservation_headroom, becomes the per-solve reservation. As many reservations as fit in vram_fraction of free device memory are admitted at once; the rest queue.
Each admitted solve runs with its reservation installed as DEVICE_MEMORY_BUDGET, so kernel_batch_size sizes kernel intermediates against that slice.

Step 3 is what makes the reservation real rather than advisory. Without it every concurrently running solve would measure the same free pool and size its batches as though it owned all of it.

A reservation larger than the whole budget is admitted alone rather than deadlocking, so a budget that turns out to be too small degrades to serial execution instead of hanging. Pass reservation = <bytes> to skip calibration when you already know the figure, or reservation = :none to disable the governor.

r = sweep.report
r.calibrated_peak    # bytes measured during the solo run
r.reservation        # bytes reserved per concurrent solve
r.max_concurrency    # how many were admitted at once
r.waits              # how many solves blocked -> the sweep was VRAM-bound

Reproducibility

The Zipper strategy draws a random sketch for its randomized range finder, so without explicit seeding a concurrent sweep would return results that depend on task scheduling. Each transformation is seeded from hash((seed, index)); both Julia's and CUDA.jl's default RNGs are task-local, so every transformation gets an independent, reproducible stream whatever order the tasks run in. Pass seed = nothing to opt out.

A bare low_energy_spectrum does no seeding — that is the caller's job. The sketch is drawn from the global RNG, so seed it yourself if you need results that are bit-identical across sessions:

using Random
Random.seed!(1234)
sol, _ = low_energy_spectrum(ctr, search_params)

In practice the range finder is not fragile: on a 2500-spin instance, six different seeds and five BLAS thread counts all returned the same energy to the last bit. Seed anyway when a number is going to be published, since the guarantee costs nothing.

Error control

A heuristic contraction is only as trustworthy as the weight it keeps. Truncating factorizations report what they discard into a task-scoped TruncationLog; sweep_transformations installs one per transformation and reports the result.

t = sweep.report.per_transform[1].truncation
t.discarded_sum   # Σᵢ εᵢ — leading-order accumulated fidelity loss
t.discarded_max   # the worst single truncation
t.saturated       # how many truncations the bond dimension forced
t.dims_kept, t.dims_offered

saturated == 0 means the bond dimension was never the binding constraint — the truncations dropped only numerically negligible singular values, so raising bond_dim will not help. A large discarded_sum means the opposite.

Two caveats on discarded_sum. It sums over every truncation in the solve, so once it approaches 1 the "accumulated fidelity loss" reading breaks down and the number only tells you the contraction is untrustworthy — a bond-4 solve of a 2500-spin instance reaches ≈ 1.3 over 2162 truncations. And it describes the contraction, not the answer: see the β warning below.

Two sweep-level numbers need no oracle to interpret:

sweep.report.energy_spread   # best-to-worst energy across transformations
sweep.report.consensus       # how many transformations reached the best energy

Eight distinct transformed contraction orders that agree provide a useful consistency check; a spread that is a sizeable fraction of the energy scale means they do not agree and flags transformation-order-sensitive outcomes, which may reflect contraction or search error. Agreement does not by itself establish overall solver convergence or solution quality.

To record truncation error for a single solve, install a log yourself:

using Base.ScopedValues: with

log = TruncationLog()
with(TRUNCATION_LOG => log) do
    low_energy_spectrum(ctr, search_params; show_progress = false)
end
truncation_stats(log)

Recording costs two device reductions per truncating factorization, so it is opt-in; with no log installed it costs nothing.

Choosing β

β is the solver's most consequential parameter. It sets how sharply the represented Boltzmann distribution concentrates on low-energy states: too small and the conditional probabilities the search branches on say little about the ground state. The optimal value depends on the instance, which the original documentation conceded without offering a way to find it.

beta_ladder walks an increasing schedule, and each step reuses the previous step's boundary MPS as the starting point for variational compression rather than rebuilding W * ψ exactly and truncating it:

ladder = beta_ladder(ctr, [2.0, 3.0, 4.0, 6.0], search_params)
sol = selected_solution(ladder)
[(s.beta, s.energy, s.truncation.discarded_sum) for s ∈ ladder.steps]

With a finite max_discarded, the ladder minimizes energy over successful rungs whose truncation.discarded_sum does not exceed the threshold. If none qualifies, it returns the minimum-energy successful rung as an untrusted fallback. The flag is available as steps[selected_index].trusted. With the default Inf, selection is simply by minimum energy. Each rung reports its discarded weight, so a scan yields evidence about both the answer and the contraction behind it.

Discarded weight does not select β

It is tempting to read max_discarded as a quality criterion. It is not, and on at least one family it points the wrong way. Ten 2500-spin square-lattice instances, bond 8, 500 states:

β	2	3	4	6	8
Σε (median)	2.3e-4	2.2e-3	3.4e-3	7.3e-4	1.4e-4
energy error (median)	7.8e-4	2.1e-4	6.2e-5	0	0

Σε is not monotone in β. It rises while the distribution is still sharpening — more structure for the boundary MPS to carry — then falls once the distribution concentrates enough to sit close to a product state, which truncates easily. Solution quality keeps improving throughout, so here the best βs carry among the lowest discarded weight. Stopping at the first rung above the threshold could miss those later rungs.

Σε answers "how much of the distribution did the contraction throw away?". That is worth knowing, and a large value is a real warning. It is not a proxy for whether the search then found a good state.

stop_when_untrusted follows from the same caveat: use it to bound cost when a rung blows up, not to conclude that higher β cannot help — a later rung may well come back under the guard.

The contractor is mutated in place at every step, so pass one the call may own.

Warm starting and the error guard do not mix

Use max_discarded with warm_start = false. The two measure different things:

a cold build forms W * ψ exactly, then truncates it — the weight it drops is a genuine discarded weight, and svd_fact records it;
a warm start optimizes the previous β's MPS within a fixed bond dimension and never performs a truncating factorization, so it reports ~zero discarded weight however accurate it actually is. Its error is a variational optimization gap, which the truncation log does not measure.

Measured on the 2048power instance over β = 1.5, 2.25, 3.0, the cold ladder reports Σε = [3.2e-3, 1.7e-4, 1.9e-5] while the warm one reports [3.2e-3, 0, 0] for identical energies — the zeros are an artifact of where the truncation happens, not better accuracy. Setting both warns for this reason.

On this instance, warm starting reduced the warmed-rung times by 24% to 25% and the complete ladder time by about 16%. Use a cold ladder to audit the contraction.

Examples

Three runnable scripts in examples/, in increasing order of scale:

beta_ladder.jl — 18 spins, runs in seconds. Error control and the β ladder, annotated inline; the place to start.
concurrent_sweep.jl — 128 spins. The transformation sweep with the device-memory governor.
square_50x50.jl — 2500 spins, the size used by the article's figures. All three features, including the agreement diagnostics on a hard instance.

API

SpinGlassPEPS.SpinGlassEngine.sweep_transformations — Function

sweep_transformations(
    build_contractor,
    sparams::SearchParameters;
    transformations,
    merge_strategy,
    symmetry,
    concurrency,
    reservation,
    reservation_headroom,
    vram_fraction,
    seed,
    blas_threads,
    diagnostics,
    show_progress
) -> SweepSolution

Solve one instance once per lattice transformation, concurrently, under a device memory budget, and return every solution together with diagnostics.

Replaces the serial for transform ∈ all_lattice_transformations loop that the published examples spell out by hand. The transformations are independent, so the only thing that limits the fan-out is device memory; this function measures what one solve costs and admits as many as fit. See the file header for why the governor is denominated in bytes.

Arguments

build_contractor: transformation -> MpsContractor. Called once per transformation, inside the task that will run it, so each solve gets its own network, projector workspace, and contraction cache.
sparams::SearchParameters: forwarded to low_energy_spectrum.

Keyword arguments

transformations = all_lattice_transformations: which transformations to run.
merge_strategy = _ -> no_merge: ctr -> strategy, because the published merge strategies close over the contractor (merge_branches(ctr; ...)).
symmetry::Symbol = :noZ2: forwarded to low_energy_spectrum.
concurrency = :auto: cap on simultaneously running solves. :auto is 1 on a GPU — a conservative default: on the tested RTX 5080 fanning these solves out over one device did not beat the serial loop, while the tested H100 did benefit. Measure before setting concurrency = 2–4 explicitly on another GPU. On CPU it is min(length(transformations), Threads.nthreads()). Raise it explicitly for a CPU-only run, several devices, or an instance you have measured. The byte budget may admit fewer.
reservation = :calibrate: bytes to reserve per solve. :calibrate measures it from a solo run; pass an integer to skip calibration (useful when you already know the figure and want all transformations to start at once), or :none to disable the governor entirely.
reservation_headroom = 1.3: multiplier on the calibrated peak, covering instance-to-instance variation between transformations.
vram_fraction = 0.85: fraction of free device memory the governor may hand out, leaving room for fragmentation and the driver context.
seed = 1234: base seed; transformation i is seeded with hash((seed, i)). nothing leaves RNG state alone, which makes a Zipper sweep non-reproducible.
blas_threads = :auto: BLAS threads per solve during the parallel phase. :auto divides the current setting by the admission limit — several solves each calling multi-threaded LAPACK (qr_fact falls back to the CPU below its shape threshold) otherwise oversubscribe the machine badly.
diagnostics = true: record truncation error. Costs two device reductions per truncating factorization.
show_progress = false: per-solve progress bars. Off by default because concurrent bars interleave into noise.

Returns

A SweepSolution. Use best_solution for the winner and .report for the diagnostics.

Example

potts_h = potts_hamiltonian(ising_graph(instance); spectrum = full_spectrum,
                            cluster_assignment_rule = super_square_lattice((m, n, t)))
params  = MpsParameters{Float64}(; bond_dim = 16, num_sweeps = 1)

sweep = sweep_transformations(
    t -> MpsContractor(
        SVDTruncate,
        PEPSNetwork{KingSingleNode{GaugesEnergy},Dense,Float64}(m, n, potts_h, t),
        params; onGPU = true, beta = 2.0, graduate_truncation = true,
    ),
    SearchParameters(; max_states = 2^8, cutoff_prob = 1e-4);
    merge_strategy = ctr -> merge_branches(ctr; merge_prob = :none),
)

sol = best_solution(sweep)
sweep.report.energy_spread   # do the transformations agree?