De-duplicate channels where possible in MPOHamiltonian constructor

[leburgel]

Description

In the current MPOHamiltonian construction, we naively add an additional channel across all relevant bonds for every term we encounter. While this has been mostly fine, for long-range interactions we very quickly end up with MPO Hamiltonians that are just much larger than they need to be. While we should probably add some support for after-the-fact compression (there's an attempt in #458, but I couldn't quite get that to work), the simplest improvement is to just not add a new channel for an interaction that shares a tensor product factor with something that's already there.

Since MPSKit doesn't use a symbolic form of the operators to construct the MPO, the checking if "things share a factor" required a bit more work, but here I attempted to implement the channel assignment to include the compression you would get from an exact left-to-right SVD compression of the corresponding symbolic finite state machine.

I think this should just always be a strict win, and for long-range interactions this can quite easily halve the resulting MPO bond dimension, so I think it's a logical first step before adding any dedicated compression algorithms.

Checklist

• Tests pass locally (julia --project=test test/runtests.jl, or the relevant subset)
• Documentation updated, if this PR changes public API (docstrings, docs/src/)
• Runic formatter is run
• Changelog entry added under [Unreleased] in docs/src/changelog.md, if this PR is user-facing (new feature, behavior change, bug fix, deprecation, or removal)

[codecov]

Codecov Report

❌ Patch coverage is 95.29412% with 4 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/operators/mpohamiltonian.jl	95.29%	4 Missing ⚠️

Files with missing lines	Coverage Δ
src/operators/mpohamiltonian.jl	92.64% <95.29%> (+0.85%)	⬆️

... and 3 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lkdvos]

Might be fair to try and insert a fast path for === tensors, such that if this ends up being a bit slow there is always a way out to have an efficient implementation that does not require a norm for each tensor? I think this entire approach scales more or less as (number of representatives) * (number of site terms), which looks a little scary to me, although this might not be all that relevant. Did you try this out on more problematic MPOs such as all-to-all quadratic or quartic?

[leburgel]

I added the === fast path that returns before any arithmetic, and reorganized so a comparison costs exactly one inner product: _decompose_once! does a pointer-only pass over the representatives first, then computes inner(O, O) once per term and reuses each representative's cached squared norm (norm_squared_new/norm_squared_old keywords). A model that reuses one operator object across its terms now evaluates no norms at all.

I checked the scaling for more and less favorable cases, and it doesn't seem too bad:

case	χ	build	alloc
quadratic, 1 structure, L=30	227 → 31	37.6 → 3.0 ms	74.2 → 2.9 MiB
quadratic, 4 structures, L=20	402 → 40	59.3 → 3.7 ms	116.7 → 4.6 MiB
quartic, 1 structure, L=12	467 → 167	79.1 → 4.3 ms	128.0 → 6.9 MiB
quadratic, all-distinct random, L=30	902 (unchanged)	37.1 → 50.2 ms	74.8 → 97.1 MiB
quartic, all-distinct random, L=12	4562 (unchanged)	104.4 → 165.3 ms	180.7 → 435.8 MiB

The last two are the worst case where literally all of the matching checks fail, and even then it doesn't seem very problematic.

[lkdvos]

See also https://github.com/lkdvos/OpSum.jl, which has to still be cleaned up a bit more but should then make all of this code somewhat obsolete :)

[leburgel]

See also https://github.com/lkdvos/OpSum.jl, which has to still be cleaned up a bit more but should then make all of this code somewhat obsolete :)

Happy to have this superseded at any point!

[lkdvos]

Scaling check for longer chains

Following up on the concern that this approach scales as (number of representatives) × (number of site terms): I measured the exponent, not just the prefactor, sweeping several sizes per model and fitting log-log over the largest three points (R² ≥ 0.98 throughout). Short answer: the exponent goes down in every regime a physical Hamiltonian produces, and the one regime where it goes up is bounded and identifiable.

Construction cost, exponent in the number of terms T

model	this PR	main	χ, PR vs main
all pairs, one operator structure, L ≤ 100	T^1.25	T^2.23	L^0.98 vs L^2.00
all quartets, one operator structure, L ≤ 14	T^0.95	T^1.55	L^3.59 vs L^4.55
infinite power law, range ≤ 80	T^0.88	T^3.42	R^0.95 vs R^1.97
all pairs, all-distinct operators, L ≤ 36	T^1.74	T^1.58	identical
all quartets, all-distinct operators, L ≤ 13	T^1.73	T^1.34	identical
infinite, all-distinct operators, range ≤ 64	T^1.94	T^3.15	identical

The load-bearing number is a count rather than a time: with a single operator structure, _decompose_once! performs exactly T-1 tensor comparisons — one per term. The number of representatives stays at 1, so the product in the cost model never materialises. With genuinely all-distinct operators it performs exactly T(T-1)/2 (verified: 2145 = C(66,2), 595 = C(35,2), 198135 = C(630,2)), which is the real quadratic term.

That is why a quantum-chemistry operator lands in the favourable column: the structure of V[pqrs] a⁺a⁺aa depends only on which of the four slots carry a creation operator — a handful of patterns — while the coefficient differs for every term. Proportionality catches exactly that, so the number of representatives is O(1) and the scan is linear in T. (Caveat: I benchmarked synthetic families with that structure, not an actual quantum-chemistry Hamiltonian.)

Absolute numbers at the largest size measured

model	main	PR at 6a7786d	PR now	allocations, main → now
all pairs, 1 structure, L=100, T=4950	7419 ms, χ=2502	66.8 ms	49.2 ms, χ=101	7.4 GiB → 36 MiB
all quartets, 1 structure, L=14, T=1001	304 ms, χ=933	15.7 ms	7.8 ms, χ=288	272 → 7.9 MiB
infinite power law, range 80	484 ms, χ=3242	2.1 ms	1.7 ms, χ=82	523 → 1.2 MiB
infinite, all distinct, range 64	200 ms, χ=8322	200.3 ms	28.7 ms, χ=8322	211 → 27.8 MiB
all pairs, all distinct, L=36	75.7 ms	201.6 ms	139.4 ms	96.7 → 149 MiB
all quartets, all distinct, L=13	236 ms	775.4 ms	743.7 ms	200 → 794 MiB

Where the time was actually going, and the three follow-up commits

Profiling put the cost somewhere unexpected: not in the new operator comparisons at all, but in the pre-existing _find_channel, which rebuilt the set of channels in use at a site from scratch on every call (unique!(last.(nonzero_keys)) allocates and hash-builds) and is called once per channel created — quadratic in the bond dimension per site. In the long-range models that scan was 66× to 253× more work than every operator comparison combined (328k scan steps versus 4.9k comparisons at L=100).

• 18bcd5eb tracks the used channels incrementally in a BitSet plus the length of the contiguous used prefix, so claiming a channel is amortised constant time and allocation-free, handing out the same index as before. The channel assignment loses a full power of T: all-quartets all-distinct T^2.32 → T^1.68, the infinite all-distinct chain T^4.13 → T^2.70.
• 2b05991c stops evaluating the norms in _proportionality before the fast paths. Julia evaluates unsupplied keyword defaults in the method prologue, so writing them out as defaults cost two inner products even on the === and mismatched-space returns — 2.0× to 2.2× a full inner product, or 5970 µs versus 0.036 µs at local dimension 8.
• 7657a1fd copies in instantiate_operator only when the indices actually need canonicalising. canonicalize_indices! rewrites the MPO in place, but both its loops are no-ops for sorted, distinct indices, and the unconditional copy meant the === shortcut never fired even once — including for a model that hands literally the same operator to all of its terms. With it reachable, comparisons roughly halve (9800 → 4949 at L=100; 3640 → 1000 for all quartets of 14) or vanish entirely (4851 → 0 when one object is reused).

Safety on that last one: the assembled Hamiltonian never stores the caller's object, because the tensors are filled in with O[...] += o, which stores the freshly computed sum. This matters because VectorInterface.scale! mutates stored blocks in place and would scale a shared block once per slot it appeared in. Checked across seven models spanning 56 to 448 filled slots: no block is reachable from more than one slot, and in-place scaling stays exact. The accompanying test also pins down that a caller-supplied FiniteMPO is left untouched on both paths.

Bond dimensions are identical before and after all three commits, in every model. operators/mpohamiltonian.jl passes 232/232.

Two things I tried that made it slower, recorded so they are not re-attempted

Caching each edge's squared norm so a comparison costs one inner product instead of three. Slower than the channel-pool alone in every model (1.1× to 2.0×): the operators in an MPO decomposition are small enough that an inner is cheaper than the bookkeeping needed to avoid it.

Rescaling the stored operators to unit norm, so norm_squared_old is a known constant with no bookkeeping at all. Also slower — 0.31× to 1.0× — and worst for large local operators, which is the opposite of what one would expect: for d=4 four-index terms it went from 15.7 ms to 95.3 ms. Rescaling in place rather than allocating narrows the gap but does not close it.

The reason is structural rather than a benchmark artifact. Normalising pays a cost per stored edge, one pass over the whole tensor, to make each comparison cheaper. But the favourable case — the one where the de-duplication works — is precisely the case with many edges and a single comparison per edge that succeeds immediately. So the price scales with χ·L while the saving scales with the comparison count, and in the models that matter the latter is the smaller number.

The one remaining superlinear term

The gap to main for all-distinct operators is entirely the _decompose_once! scan, untouched by the above: exactly T(T-1)/2 comparisons at roughly 2.75 µs each for a four-site operator, so T = 10⁴ is about 2.3 minutes, T = 3·10⁴ about 21 minutes, and T = 10⁵ about 3.8 hours. Nothing physical lands there — it needs every term's operator to be non-proportional to every other — but if a guard seems worthwhile, bucketing representatives by space(O) is free and would also help mixed-structure models, and a cap on the scan length would bound the worst case to main's behaviour rather than a quadratic.

_{Method: Julia 1.12.6, two scratch environments with byte-identical resolved versions of all 83 dependencies, MPSKit taken from this branch versus main at 4b57994. Physical space ℂ^2 except where d=4 is noted. Timings are best-of-n on a shared machine, so treat sub-millisecond entries as noise-dominated; the exponents and the comparison counts are the robust part.}

🤖 Generated with Claude Code

[lkdvos]

I did some more benchmarks myself, mostly because I was curious about the scaling, rather than the absolute timings for a given L. The results should be commented here somewhere, but in any case it seems like this is overall a huge improvement, and I can't really think of a case where someone wants a huge number of completely unrelated terms, and in that case they can still be manually deduplicated externally, so for me this looks great. I slightly improved some of the overhead for the exact case where there are a lot of representatives with the channelpool, giving a minor speedup, and did some work to make sure the === fast path actually is hit, since there is a copy in the instantiate_operator before this was otherwise never hit.

Thanks for both investigating and solving this issue!

[lkdvos]

Documentation error is unrelated, mostly due to v1.13 changes