From 809c81092ff8913ac0704b8c2036152ea354df58 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 8 Sep 2026 11:03:56 +0200 Subject: [PATCH 1/5] De-duplicate channels where possible in `MPOHamiltonian` constructor --- src/operators/mpohamiltonian.jl | 154 ++++++++++++++++++++++++++----- test/operators/mpohamiltonian.jl | 82 ++++++++++++++++ 2 files changed, 214 insertions(+), 22 deletions(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 1ca4bacb6..974a5134b 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -425,6 +425,132 @@ function _find_channel(nonzero_keys; init = 2) return max(maximum(range) + 1, init) end +""" + _proportionality(O_new, O_old; tol) -> λ or nothing + +Return the scalar `λ` for which `O_new ≈ λ * O_old`, or `nothing` if the two operators are not +proportional. Both arguments are entries of a decomposed local MPO, i.e. either `MPOTensor`s or +scalar multiples of the identity. +""" +function _proportionality(O_new::Number, O_old::Number; tol = nothing) + iszero(O_old) && return nothing + return O_new / O_old +end +_proportionality(::Number, ::AbstractTensorMap; tol = nothing) = nothing +_proportionality(::AbstractTensorMap, ::Number; tol = nothing) = nothing +function _proportionality( + O_new::AbstractTensorMap, O_old::AbstractTensorMap; + tol = eps(real(promote_type(scalartype(O_new), scalartype(O_old))))^(3 / 4) + ) + space(O_new) == space(O_old) || return nothing + norm_new = norm(O_new) + norm²_old = inner(O_old, O_old) + (iszero(norm²_old) || iszero(norm_new)) && return nothing + # note that dividing by `inner(O_old, O_old)` instead of `norm(O_old)^2` avoids a + # roundtrip through `sqrt`, such that identical operators give `λ = 1` exactly + λ = inner(O_old, O_new) / norm²_old + norm(add(O_new, O_old, -λ)) ≤ tol * norm_new || return nothing + return λ +end + +""" + _instantiate_operators(lattice, local_operators) + +Instantiate all `local_operators` on `lattice`, decomposing every distinct operator into an MPO +only once. Operators that are proportional to one another share a single decomposition, with +the scalar factor absorbed into the final tensor. Terms that differ only in their prefactor +thus end up with literally the same operators, and can share their virtual channels in +[`_assign_channels!`](@ref) -- which the decomposition itself does not guarantee, as the SVDs +of two proportional operators need not be related by that same scalar. +""" +function _instantiate_operators(lattice, local_operators) + representatives = Tuple{Any, Any}[] + return map(collect(local_operators)) do term + return instantiate_operator(lattice, _decompose_once!(representatives, term)) + end +end + +# terms that are not supplied as an `inds => operator` pair -- e.g. a `LocalOperator` -- and +# operators that the caller already decomposed into an MPO are passed along untouched +_decompose_once!(representatives, term) = term +function _decompose_once!(representatives, term::Pair) + inds, O = term + return inds => _decompose_once!(representatives, O) +end +function _decompose_once!(representatives, O::AbstractTensorMap) + for (O_rep, mpo) in representatives + λ = _proportionality(O, O_rep) + isnothing(λ) && continue + isone(λ) && return mpo + Os = parent(mpo) + return FiniteMPO([i == lastindex(Os) ? λ * Os[i] : Os[i] for i in eachindex(Os)]) + end + mpo = FiniteMPO(O) + push!(representatives, (O, mpo)) + return mpo +end + +""" + _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) + +Distribute the instantiated local MPOs `local_mpos` over the virtual channels of the Jordan +block form, storing the resulting graph of operators in `nonzero_keys` and `nonzero_opps`. + +Terms are inserted one at a time, from left to right. Whenever the next operator of a term +coincides -- up to a scalar factor -- with an edge that is already present at that site and +starts from the same channel, that channel is reused instead of opening up a new one, and the +scalar is carried along to be absorbed into the final operator of the term. This is what keeps +the bond dimension of e.g. a long-range interaction linear instead of quadratic in its range: +all terms that start out the same way share a single channel until they part ways. + +Because a channel is only ever created for a unique combination of site, incoming channel and +operator, every channel is reached by exactly one sequence of operators. Sharing channels +between terms can therefore never generate paths that do not correspond to a requested term. +""" +function _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) + L = length(nonzero_keys) + # index the outgoing edges at every site by their incoming channel, such that a new + # operator only has to be compared against the few edges it could possibly share + outgoing = [Dict{Int, Vector{Int}}() for _ in 1:L] + + for (sites, local_mpo) in local_mpos + key_R = 1 + coeff = 1 + for (i, (site, O)) in enumerate(zip(sites, local_mpo)) + key_L = i == 1 ? 1 : key_R + if i == length(local_mpo) + # the final operator always drops back onto the last channel, and is where the + # accumulated scalar of all shared edges ends up + push!(nonzero_keys[site], (key_L, 0)) + push!(nonzero_opps[site], isone(coeff) ? O : coeff * O) + break + end + + edges = get!(Vector{Int}, outgoing[mod1(site, L)], key_L) + shared = nothing + for edge in edges + λ = _proportionality(O, nonzero_opps[site][edge]) + if !isnothing(λ) + shared = (last(nonzero_keys[site][edge]), λ) + break + end + end + + if isnothing(shared) + key_R = _find_channel(nonzero_keys[site]; init = key_L) + push!(nonzero_keys[site], (key_L, key_R)) + push!(nonzero_opps[site], O) + push!(edges, length(nonzero_opps[site])) + else + key_R, λ = shared + coeff = coeff * λ + end + end + end + + return nonzero_keys, nonzero_opps +end + function _find_first_mpotensor(Ws) for W in Ws for x in W @@ -450,19 +576,10 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_opera # partial sort by interaction range local_mpos = sort!( - map(Base.Fix1(instantiate_operator, lattice), collect(local_operators)); - by = x -> length(x[1]) + _instantiate_operators(lattice, local_operators); by = x -> length(x[1]) ) - for (sites, local_mpo) in local_mpos - local key_R # trick to define key_R before the first iteration - for (i, (site, O)) in enumerate(zip(sites, local_mpo)) - key_L = i == 1 ? 1 : key_R - key_R = i == length(local_mpo) ? 0 : _find_channel(nonzero_keys[site]; init = key_L) - push!(nonzero_keys[site], (key_L, key_R)) - push!(nonzero_opps[site], O) - end - end + _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) # construct the sparse MPO T = _find_tensortype(nonzero_opps) @@ -532,18 +649,11 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, local_ end # partial sort by interaction range - unsorted_mpos = map(Base.Fix1(instantiate_operator, lattice), [local_operators...]) - local_mpos = sort!(unsorted_mpos; by = x -> length(x[1])) + local_mpos = sort!( + _instantiate_operators(lattice, local_operators); by = x -> length(x[1]) + ) - for (sites, local_mpo) in local_mpos - local key_R # trick to define key_R before the first iteration - for (i, (site, O)) in enumerate(zip(sites, local_mpo)) - key_L = i == 1 ? 1 : key_R - key_R = i == length(local_mpo) ? 0 : _find_channel(nonzero_keys[site]; init = key_L) - push!(nonzero_keys[site], (key_L, key_R)) - push!(nonzero_opps[site], O) - end - end + _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) # construct the sparse MPO T = _find_tensortype(nonzero_opps) diff --git a/test/operators/mpohamiltonian.jl b/test/operators/mpohamiltonian.jl index 29e615e2d..08672f4c9 100644 --- a/test/operators/mpohamiltonian.jl +++ b/test/operators/mpohamiltonian.jl @@ -286,3 +286,85 @@ end @test expectation_value(mps2, H3) ≈ expectation_value(mps2, H4) end end + +@testset "MPOHamiltonian channel sharing" begin + # Terms that start out identically share a virtual channel, which is what keeps the bond + # dimension of a long-range interaction linear -- instead of quadratic -- in its range. + # Correctness is checked against the same Hamiltonian assembled by summing single-term + # Hamiltonians, which cannot share anything and is thus independent of how channels are + # assigned. + Z = TensorMap(Float64[1 0; 0 -1], ℂ^2, ℂ^2) + X = TensorMap(Float64[0 1; 1 0], ℂ^2, ℂ^2) + + infinite_ev(lattice, terms) = let L = length(lattice) + ψ = InfiniteMPS(collect(lattice), fill(ℂ^8, L)) + H = InfiniteMPOHamiltonian(lattice, terms...) + Href = sum(InfiniteMPOHamiltonian(lattice, t) for t in terms) + (H, sum(expectation_value(ψ, H)), sum(expectation_value(ψ, Href))) + end + + @testset "power-law couplings, unit cell L=$L" for L in (1, 2, 3) + R = 5 + lattice = fill(ℂ^2, L) + terms = Pair[(i, i + d) => (1 / d^3) * (Z ⊗ Z) for i in 1:L for d in 1:R] + append!(terms, [(i,) => 0.7 * X for i in 1:L]) + H, e, eref = infinite_ev(lattice, terms) + @test e ≈ eref + @test all(i -> dim(left_virtualspace(H, i)) ≤ R + 2, 1:L) + end + + @testset "proportional terms are merged" begin + lattice = fill(ℂ^2, 1) + H = InfiniteMPOHamiltonian(lattice, (1, 2) => Z ⊗ Z, (1, 2) => 0.5 * (Z ⊗ Z)) + Href = InfiniteMPOHamiltonian(lattice, (1, 2) => 1.5 * (Z ⊗ Z)) + ψ = InfiniteMPS([ℂ^2], [ℂ^8]) + @test sum(expectation_value(ψ, H)) ≈ sum(expectation_value(ψ, Href)) + @test dim(left_virtualspace(H, 1)) == dim(left_virtualspace(Href, 1)) + end + + @testset "negative prefactors" begin + lattice = fill(ℂ^2, 1) + terms = Pair[ + (1, 2) => Z ⊗ Z, (1, 3) => -2.0 * (Z ⊗ Z), (1, 4) => 0.25 * (Z ⊗ Z), + ] + _, e, eref = infinite_ev(lattice, terms) + @test e ≈ eref + end + + @testset "mixed operators and ranges" begin + lattice = fill(ℂ^2, 2) + terms = Pair[ + (1, 2) => Z ⊗ Z, (1, 3) => Z ⊗ X, (1, 4) => 3.0 * (Z ⊗ Z), + (2, 3, 4) => X ⊗ Z ⊗ X, (1,) => 0.7 * X, + ] + _, e, eref = infinite_ev(lattice, terms) + @test e ≈ eref + end + + @testset "finite chain" begin + L = 6 + lattice = fill(ℂ^2, L) + terms = Pair[ + (i, j) => (1 / (j - i)^3) * (Z ⊗ Z) for i in 1:L for j in (i + 1):L + ] + push!(terms, (1,) => 0.7 * X) + H = FiniteMPOHamiltonian(lattice, terms...) + Href = sum(FiniteMPOHamiltonian(lattice, t) for t in terms) + @test convert(TensorMap, H) ≈ convert(TensorMap, Href) + @test maximum(i -> dim(left_virtualspace(H, i)), 1:L) < + maximum(i -> dim(left_virtualspace(Href, i)), 1:L) + end + + @testset "symmetric tensors $(sectortype(pspace))" for (pspace, Dspace) in + zip(pspaces, vspaces) + O = rand(ComplexF64, pspace^2, pspace^2) + O += O' + terms = Pair[(1, 1 + d) => (1 / d^3) * O for d in 1:4] + lattice = fill(pspace, 1) + H = InfiniteMPOHamiltonian(lattice, terms...) + Href = sum(InfiniteMPOHamiltonian(lattice, t) for t in terms) + ψ = InfiniteMPS([pspace], [Dspace]) + @test sum(expectation_value(ψ, H)) ≈ sum(expectation_value(ψ, Href)) + @test dim(left_virtualspace(H, 1)) < dim(left_virtualspace(Href, 1)) + end +end From 6a7786dd9049e2cddd21fd4a7429ac5abbe18106 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 10 Sep 2026 09:51:24 +0200 Subject: [PATCH 2/5] Address review comments and update changelog --- docs/src/changelog.md | 4 ++ src/operators/mpohamiltonian.jl | 91 +++++++++++++++++++++++---------- 2 files changed, 69 insertions(+), 26 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index eeae4041f..7ba197845 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -41,6 +41,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Changed +- `FiniteMPOHamiltonian`/`InfiniteMPOHamiltonian` constructed from a set of local terms now share + virtual channels between terms that start out with the same operators, up to a scalar factor, + and add up terms that are linearly dependent. The resulting Hamiltonian is unchanged, but its + bond dimension is generally smaller ([#518](https://github.com/QuantumKitHub/MPSKit.jl/pull/518)) - Renormalization during time evolution is now controlled by an explicit `normalize` keyword on `timestep`/`time_evolve` (default `false`), decoupled from `imaginary_evolution`. By default the norm is preserved, so it retains useful information (the accumulated truncation error in real time, diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 974a5134b..a34b02ef5 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -426,30 +426,38 @@ function _find_channel(nonzero_keys; init = 2) end """ - _proportionality(O_new, O_old; tol) -> λ or nothing + _proportionality(O_new, O_old; tol, norm_squared_new, norm_squared_old) -> λ or nothing Return the scalar `λ` for which `O_new ≈ λ * O_old`, or `nothing` if the two operators are not proportional. Both arguments are entries of a decomposed local MPO, i.e. either `MPOTensor`s or scalar multiples of the identity. + +Operators that are the very same object are recognised without touching their entries, and +`norm_squared_new`/`norm_squared_old` allow a caller that compares one operator against many to +hoist the squared norms out of its loop. What remains is a single inner product per comparison. """ -function _proportionality(O_new::Number, O_old::Number; tol = nothing) +function _proportionality(O_new::Number, O_old::Number; kwargs...) iszero(O_old) && return nothing return O_new / O_old end -_proportionality(::Number, ::AbstractTensorMap; tol = nothing) = nothing -_proportionality(::AbstractTensorMap, ::Number; tol = nothing) = nothing +_proportionality(::Number, ::AbstractTensorMap; kwargs...) = nothing +_proportionality(::AbstractTensorMap, ::Number; kwargs...) = nothing function _proportionality( O_new::AbstractTensorMap, O_old::AbstractTensorMap; - tol = eps(real(promote_type(scalartype(O_new), scalartype(O_old))))^(3 / 4) + tol = eps(real(promote_type(scalartype(O_new), scalartype(O_old))))^(3 / 4), + norm_squared_new = real(inner(O_new, O_new)), + norm_squared_old = real(inner(O_old, O_old)) ) + # an operator is trivially proportional to itself, which is the common case as soon as a + # decomposition is shared between terms -- no arithmetic needed + O_new === O_old && return one(scalartype(O_new)) space(O_new) == space(O_old) || return nothing - norm_new = norm(O_new) - norm²_old = inner(O_old, O_old) - (iszero(norm²_old) || iszero(norm_new)) && return nothing + (iszero(norm_squared_old) || iszero(norm_squared_new)) && return nothing # note that dividing by `inner(O_old, O_old)` instead of `norm(O_old)^2` avoids a # roundtrip through `sqrt`, such that identical operators give `λ = 1` exactly - λ = inner(O_old, O_new) / norm²_old - norm(add(O_new, O_old, -λ)) ≤ tol * norm_new || return nothing + ip = inner(O_old, O_new) + λ = ip / norm_squared_old + norm(add(O_new, O_old, -λ)) ≤ tol * sqrt(norm_squared_new) || return nothing return λ end @@ -464,7 +472,7 @@ thus end up with literally the same operators, and can share their virtual chann of two proportional operators need not be related by that same scalar. """ function _instantiate_operators(lattice, local_operators) - representatives = Tuple{Any, Any}[] + representatives = Tuple{Any, Any, Any}[] return map(collect(local_operators)) do term return instantiate_operator(lattice, _decompose_once!(representatives, term)) end @@ -478,15 +486,25 @@ function _decompose_once!(representatives, term::Pair) return inds => _decompose_once!(representatives, O) end function _decompose_once!(representatives, O::AbstractTensorMap) - for (O_rep, mpo) in representatives - λ = _proportionality(O, O_rep) + # an operator that is handed in more than once -- the same object, as happens whenever a + # model reuses a single operator across its terms -- is recognised without any arithmetic + for (O_rep, _, mpo) in representatives + O === O_rep && return mpo + end + # otherwise every representative costs one inner product, with its own squared norm + # cached so that it is not recomputed for every term + norm_squared_new = real(inner(O, O)) + for (O_rep, norm_squared_rep, mpo) in representatives + λ = _proportionality( + O, O_rep; norm_squared_new, norm_squared_old = norm_squared_rep + ) isnothing(λ) && continue isone(λ) && return mpo Os = parent(mpo) return FiniteMPO([i == lastindex(Os) ? λ * Os[i] : Os[i] for i in eachindex(Os)]) end mpo = FiniteMPO(O) - push!(representatives, (O, mpo)) + push!(representatives, (O, norm_squared_new, mpo)) return mpo end @@ -503,44 +521,65 @@ scalar is carried along to be absorbed into the final operator of the term. This the bond dimension of e.g. a long-range interaction linear instead of quadratic in its range: all terms that start out the same way share a single channel until they part ways. +The same comparison is made for the operators that terminate a term, where a match means that +the two terms are linearly dependent: those are added together into a single edge, instead of +two edges that are only summed when the tensors are filled in. + Because a channel is only ever created for a unique combination of site, incoming channel and operator, every channel is reached by exactly one sequence of operators. Sharing channels between terms can therefore never generate paths that do not correspond to a requested term. """ function _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) L = length(nonzero_keys) - # index the outgoing edges at every site by their incoming channel, such that a new - # operator only has to be compared against the few edges it could possibly share + # index the edges at every site by their incoming channel, such that a new operator only + # has to be compared against the few edges it could possibly share. Edges terminating on + # `IdR` are kept separate, since those can only ever be merged with one another outgoing = [Dict{Int, Vector{Int}}() for _ in 1:L] + terminating = [Dict{Int, Vector{Int}}() for _ in 1:L] for (sites, local_mpo) in local_mpos key_R = 1 coeff = 1 for (i, (site, O)) in enumerate(zip(sites, local_mpo)) key_L = i == 1 ? 1 : key_R + keys_site, opps_site = nonzero_keys[site], nonzero_opps[site] + if i == length(local_mpo) - # the final operator always drops back onto the last channel, and is where the - # accumulated scalar of all shared edges ends up - push!(nonzero_keys[site], (key_L, 0)) - push!(nonzero_opps[site], isone(coeff) ? O : coeff * O) + # the final operator always drops back onto the last channel, and is where + # the accumulated scalar of all shared edges ends up + O_final = isone(coeff) ? O : coeff * O + edges = get!(Vector{Int}, terminating[mod1(site, L)], key_L) + merged = false + for edge in edges + λ = _proportionality(O_final, opps_site[edge]) + isnothing(λ) && continue + opps_site[edge] = (1 + λ) * opps_site[edge] + merged = true + break + end + if !merged + push!(keys_site, (key_L, 0)) + push!(opps_site, O_final) + push!(edges, length(opps_site)) + end break end edges = get!(Vector{Int}, outgoing[mod1(site, L)], key_L) shared = nothing for edge in edges - λ = _proportionality(O, nonzero_opps[site][edge]) + λ = _proportionality(O, opps_site[edge]) if !isnothing(λ) - shared = (last(nonzero_keys[site][edge]), λ) + shared = (last(keys_site[edge]), λ) break end end if isnothing(shared) - key_R = _find_channel(nonzero_keys[site]; init = key_L) - push!(nonzero_keys[site], (key_L, key_R)) - push!(nonzero_opps[site], O) - push!(edges, length(nonzero_opps[site])) + key_R = _find_channel(keys_site; init = key_L) + push!(keys_site, (key_L, key_R)) + push!(opps_site, O) + push!(edges, length(opps_site)) else key_R, λ = shared coeff = coeff * λ From 18bcd5eb5e5e40d7a6b418668df6ca58dd43b0d8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 10 Sep 2026 10:19:55 -0400 Subject: [PATCH 3/5] Track used channels incrementally in `_assign_channels!` `_find_channel` rebuilt the set of channels in use at a site from scratch on every call: `unique!(last.(nonzero_keys))` allocates a vector as long as that site's key list and hash-builds a set, which is then scanned linearly. As it is called once per channel created, the cost is quadratic in the bond dimension per site, and it dominated the channel assignment for every model measured -- by two orders of magnitude for a long-range interaction, where the scan is 66x to 250x more work than all of the operator comparisons combined. `ChannelPool` keeps the same set in a `BitSet` alongside the length of the contiguous prefix of used channels, so `claim_channel!` is amortised constant time and allocates nothing. It hands out the same index as before -- the smallest free channel not below the incoming one -- so the resulting Hamiltonian, bond dimension included, is unchanged. Full-constructor timings, versus the state before this commit: | model | before | after | | ------------------------------------ | ------ | ------ | | all pairs, one structure, L=100 | 66.8ms | 49.9ms | | all quartets, one structure, L=14 | 15.7ms | 8.6ms | | infinite power law, range 80 | 2.1ms | 1.7ms | | infinite, all distinct, range 64 | 200ms | 62.7ms | | all pairs, all distinct, L=36 | 202ms | 141ms | | all quartets, all distinct, L=13 | 775ms | 638ms | The channel assignment loses a full power of the term count: all-quartets all-distinct goes from `T^2.32` to `T^1.68`, the infinite all-distinct chain from `T^4.12` to `T^2.70`, and constructing an all-quartets Hamiltonian with a single operator structure is now linear in the number of terms. Co-Authored-By: Claude Opus 5 (1M context) --- src/operators/mpohamiltonian.jl | 49 +++++++++++++++++++++++++++------ 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index a34b02ef5..346cfe16e 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -415,14 +415,45 @@ function _find_tensortype(nonzero_operators::AbstractArray) end end -function _find_channel(nonzero_keys; init = 2) - init = max(init, 2) - range = unique!(last.(nonzero_keys)) - isempty(range) && return init - for i in init:max(maximum(range), 2) - i ∉ range && return i +""" + ChannelPool() + +The set of outgoing virtual channels in use at a single site of a Jordan block MPO. + +Channels are handed out by [`claim_channel!`](@ref), which returns the smallest index that is +still free, so that the channel indices stay dense and the resulting bond dimension is not +inflated by gaps. +""" +mutable struct ChannelPool + used::BitSet + # all of `2:contiguous` are in use, so `contiguous + 1` is the smallest free channel and + # answers every request that does not start above it. Tracking this incrementally is what + # keeps `claim_channel!` amortised constant time instead of linear in the bond dimension + contiguous::Int +end +ChannelPool() = ChannelPool(BitSet(), 1) + +""" + claim_channel!(pool::ChannelPool, init::Int) -> key + +Reserve and return the smallest channel index `≥ max(init, 2)` that is not yet in use. The +lower bound is what keeps the Jordan block form upper triangular: an edge leaving channel +`init` may not drop back to a lower channel. +""" +function claim_channel!(pool::ChannelPool, init::Int) + key = max(init, 2) + if key ≤ pool.contiguous + 1 + key = pool.contiguous + 1 + else + while key in pool.used + key += 1 + end + end + push!(pool.used, key) + while (pool.contiguous + 1) in pool.used + pool.contiguous += 1 end - return max(maximum(range) + 1, init) + return key end """ @@ -536,6 +567,8 @@ function _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) # `IdR` are kept separate, since those can only ever be merged with one another outgoing = [Dict{Int, Vector{Int}}() for _ in 1:L] terminating = [Dict{Int, Vector{Int}}() for _ in 1:L] + # the channels already handed out at every site + pools = [ChannelPool() for _ in 1:L] for (sites, local_mpo) in local_mpos key_R = 1 @@ -576,7 +609,7 @@ function _assign_channels!(nonzero_keys, nonzero_opps, local_mpos) end if isnothing(shared) - key_R = _find_channel(keys_site; init = key_L) + key_R = claim_channel!(pools[mod1(site, L)], key_L) push!(keys_site, (key_L, key_R)) push!(opps_site, O) push!(edges, length(opps_site)) From 2b05991c481f207c46c052eda1ff6de10e9f70cc Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 10 Sep 2026 10:20:13 -0400 Subject: [PATCH 4/5] Do not evaluate the norms in `_proportionality` before they are needed MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Julia evaluates the defaults of unsupplied keyword arguments in the method prologue, before the body runs, so spelling the squared norms out as defaults made them cost two inner products even on the two paths that return before reading them: the `===` shortcut and the mismatched-space check. Neither the compiler nor the caller can avoid that -- the work is emitted unconditionally. Verified on Julia 1.12.6: with the defaults in the signature, a call on the `===` path costs 2.0x to 2.2x a full inner product across local dimensions 2, 4 and 8 (5970us versus 0.036us at dimension 8), and a mismatched-space call costs one. Filling the norms in inside the body instead leaves both paths free, and every path that does reach the computation is unchanged -- `λ` agrees exactly. This is a latent rather than a measured win for the models benchmarked so far, whose operators all share one space and never hit the `===` shortcut, but it removes the only reason those two early returns were not already free. Co-Authored-By: Claude Opus 5 (1M context) --- src/operators/mpohamiltonian.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 346cfe16e..5bfb17723 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -476,16 +476,17 @@ _proportionality(::AbstractTensorMap, ::Number; kwargs...) = nothing function _proportionality( O_new::AbstractTensorMap, O_old::AbstractTensorMap; tol = eps(real(promote_type(scalartype(O_new), scalartype(O_old))))^(3 / 4), - norm_squared_new = real(inner(O_new, O_new)), - norm_squared_old = real(inner(O_old, O_old)) + norm_squared_new = nothing, norm_squared_old = nothing ) # an operator is trivially proportional to itself, which is the common case as soon as a # decomposition is shared between terms -- no arithmetic needed O_new === O_old && return one(scalartype(O_new)) space(O_new) == space(O_old) || return nothing + + # evaluate norms only after fast paths are checked + norm_squared_old = @something norm_squared_old real(inner(O_old, O_old)) + norm_squared_new = @something norm_squared_new real(inner(O_new, O_new)) (iszero(norm_squared_old) || iszero(norm_squared_new)) && return nothing - # note that dividing by `inner(O_old, O_old)` instead of `norm(O_old)^2` avoids a - # roundtrip through `sqrt`, such that identical operators give `λ = 1` exactly ip = inner(O_old, O_new) λ = ip / norm_squared_old norm(add(O_new, O_old, -λ)) ≤ tol * sqrt(norm_squared_new) || return nothing From c16ec235c36b5ed4eb9330a40f0e02b82f548dc3 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 10 Sep 2026 11:04:22 -0400 Subject: [PATCH 5/5] Only copy in `instantiate_operator` when the indices need canonicalising `_decompose_once!` hands the same decomposition to every term that shares an operator up to a scalar, but `instantiate_operator` then copied it -- and `Base.copy` of an MPO copies every tensor -- so no two terms ever held the same object and the `===` shortcut in `_proportionality` could not fire once. Measured: zero hits, in every model, including one that passes literally the same operator to all of its terms. The copy is there because `canonicalize_indices!` rewrites the MPO in place through `swap!` and `multiply_neighbours!`. Both of its loops are skipped for sorted, distinct indices though, which is the common case, so the copy is only needed for the others. This is safe because the assembled Hamiltonian never stores the caller's object: the tensors are filled in with `O[...] += o`, which stores the freshly computed sum. That matters, as `VectorInterface.scale!` mutates stored blocks in place and would scale a shared block once per slot it appeared in. Checked across seven models covering 56 to 448 filled slots: no block is reachable from more than one slot, and in-place scaling stays exact. With the shortcut reachable, the comparisons in `_assign_channels!` roughly halve -- 9800 to 4949 for all pairs of a chain of 100, 3640 to 1000 for all quartets of 14, and 4851 to none at all when one operator object is reused throughout. Construction gains 6% on those models and is unchanged on models whose operators are all distinct, where nothing can be shared in the first place; allocations drop by the copies that are no longer made (17.0 MiB to 7.9 MiB for all quartets of 14). The elision also drops what used to be an unconditional defensive copy of a caller-supplied `FiniteMPO`, so the accompanying test checks that neither path writes back into the operator it was given. Co-Authored-By: Claude Opus 5 (1M context) --- src/operators/mpohamiltonian.jl | 10 +++++++++- test/operators/mpohamiltonian.jl | 10 ++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 5bfb17723..4d26b165e 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -351,7 +351,6 @@ function instantiate_operator(state::AbstractMPS, O::Pair) end function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O)::Pair) inds = inds′ isa Int ? [inds′] : inds′ - mpo = O isa FiniteMPO ? copy(O) : FiniteMPO(O) # convert to linear index type indices = Vector{Int}(undef, length(inds)) @@ -359,6 +358,15 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O indices[i] = Base._to_linear_index(lattice, Tuple(inds[i])...) # this should mean all inds are valid... end + # _proportionality check wants `===` as fast path, so only copy if really needed, i.e. if canonicalize_indices does something. + mpo = if !(O isa FiniteMPO) + FiniteMPO(O) + elseif issorted(indices) && allunique(indices) + O + else + copy(O) + end + # sort indices and deduplicate indices, mpo = canonicalize_indices!(indices, mpo) operators = parent(mpo) diff --git a/test/operators/mpohamiltonian.jl b/test/operators/mpohamiltonian.jl index 08672f4c9..6083b74c1 100644 --- a/test/operators/mpohamiltonian.jl +++ b/test/operators/mpohamiltonian.jl @@ -206,6 +206,16 @@ end H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO(X ⊗ (X * Y) ⊗ Y)) H2 = FiniteMPOHamiltonian(chain, (1, 2, 2, 3) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y)) @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + # a caller-supplied MPO is only copied when the indices actually need canonicalising, + # so check that neither path writes back into it + for inds in ((1, 2, 3), (1, 2, 2, 3), (2, 3, 1, 1)) + mpo = FiniteMPO(X ⊗ X ⊗ Y ⊗ Y) + length(inds) == 3 && (mpo = FiniteMPO((X * X) ⊗ Y ⊗ Y)) + before = map(copy, parent(mpo)) + FiniteMPOHamiltonian(chain, inds => mpo) + @test all(splat(≈), zip(parent(mpo), before)) + end end @testset "InfiniteMPOHamiltonian $(sectortype(pspace))" for (pspace, Dspace) in zip(pspaces, vspaces)