From 5f606d7d17506a663471876a43b08fd75ccc4225 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 12:21:25 -0400 Subject: [PATCH 1/6] Working on it --- .../LoopExpansion/free_energy.jl | 84 +++++++++ .../LoopExpansion/graph_clusters.jl | 177 ++++++++++++++++++ .../{ => LoopExpansion}/loopcorrection.jl | 4 +- src/expect.jl | 53 ++++++ 4 files changed, 315 insertions(+), 3 deletions(-) create mode 100644 src/MessagePassing/LoopExpansion/free_energy.jl create mode 100644 src/MessagePassing/LoopExpansion/graph_clusters.jl rename src/MessagePassing/{ => LoopExpansion}/loopcorrection.jl (94%) diff --git a/src/MessagePassing/LoopExpansion/free_energy.jl b/src/MessagePassing/LoopExpansion/free_energy.jl new file mode 100644 index 00000000..4d541e5b --- /dev/null +++ b/src/MessagePassing/LoopExpansion/free_energy.jl @@ -0,0 +1,84 @@ +# Free-energy (generating-function) estimate of a single-site observable, exposed through +# `expect(...; alg = "loopcorrections")`. +# +# ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_{ε=0} ≈ [F(ε) − F(−ε)] / (2ε), +# F(t) = ln⟨ψ|e^{t Ô}|ψ⟩ = ln‖e^{t Ô/2}|ψ⟩‖² (Hermitian Ô), +# +# Each F is the loop-corrected free energy of a *genuine norm network*: the gate +# e^{±ε Ô/2} is absorbed (un-normalized) into the ket at the observable site and BP is +# re-solved, so the loop expansion prunes ALL leaves — there is no protected operator +# vertex, hence no "anomalous" leaf-containing cluster, which is what tends to make this +# estimator converge smoothly. BP is re-solved on each shifted network so the loop +# corrections also pick up the linear response of the messages. + +# F(α-shifted) = loop-corrected free energy of e^{α Ô}|ψ⟩, built by absorbing the +# (un-normalized) one-site gate into the ket and re-solving BP for the perturbed network. +# `F = ln Z_BP + Σ_C w_C` is the additive linked-cluster free energy (`loopcorrected_free_energy`), +# the genuine extensive log-partition-function whose smooth ε-dependence makes this estimator +# converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which +# only agrees with it to O(w) and resums the same clusters multiplicatively instead. +function _gated_loop_free_energy( + ψ::TensorNetworkState, op_string::String, v, α, max_configuration_size::Integer; + cache_update_kwargs, + ) + s = only(siteinds(ψ)[v]) + G = ITensors.exp(α * ITensors.op(op_string, s); ishermitian = true) + bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) + # `normalize_tensors = false`: the un-normalized gated tensor is exactly what makes the + # squared norm equal the partition function ⟨ψ|e^{2α Ô}|ψ⟩ we want to differentiate. + bpc, _ = apply_gate!(G, bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) + bpc = update(bpc; cache_update_kwargs...) + return loopcorrected_free_energy(bpc, max_configuration_size) +end + +""" + expect(ψ, obs; alg = "loopcorrections", max_configuration_size, ε = 1e-4, cache_update_kwargs...) + +Free-energy / generating-function estimate of a **single-site Hermitian** observable +`obs = (op, vertex[, coeff])`, + + ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_0 ≈ [F(ε) − F(−ε)] / (2ε), + +with `F(t) = ln‖e^{t Ô/2}|ψ⟩‖²` the loop-corrected free energy of that norm network +(`loopcorrected_free_energy`, the additive linked-cluster form `ln Z_BP + Σ_C w_C`; +`max_configuration_size` counts EDGES). Because each `F` is the free energy of a genuine +norm network, every loop cluster is leaf-free and there is no protected operator vertex, so +this typically converges smoothly (little oscillation). BP is re-solved for each shifted +network so the loop corrections include the linear response of the messages. + +`ε` is the central finite-difference step (default `1e-4`). `Ô` must be Hermitian (the +generating-function identity `⟨ψ|e^{εÔ}|ψ⟩ = ‖e^{εÔ/2}|ψ⟩‖²` relies on `Ô† = Ô`). Accepts a +`TensorNetworkState` or a `BeliefPropagationCache`, and a single observable or a vector of +them; multi-site observables are not supported. +""" +function expect( + ::Algorithm"loopcorrections", + ψ::TensorNetworkState, + obs::Tuple; + max_configuration_size::Integer, ε::Real = 1e-4, + cache_update_kwargs = default_bp_update_kwargs(ψ), + ) + op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ)) + iszero(coeff) && return zero(coeff) + length(obs_vs) == 1 || + error("expect(...; alg = \"loopcorrections\") supports single-site observables only; got $(length(obs_vs)) sites.") + v = only(obs_vs) + op_string = only(op_strings) + + Fp = _gated_loop_free_energy(ψ, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) + Fm = _gated_loop_free_energy(ψ, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) + return coeff * (Fp - Fm) / (2ε) +end + +function expect(alg::Algorithm"loopcorrections", cache::BeliefPropagationCache, obs::Tuple; kwargs...) + return expect(alg, network(cache), obs; kwargs...) +end + +function expect( + alg::Algorithm"loopcorrections", + ψ::Union{TensorNetworkState, BeliefPropagationCache}, + observables::Vector{<:Tuple}; + kwargs..., + ) + return map(obs -> expect(alg, ψ, obs; kwargs...), observables) +end diff --git a/src/MessagePassing/LoopExpansion/graph_clusters.jl b/src/MessagePassing/LoopExpansion/graph_clusters.jl new file mode 100644 index 00000000..036e07d7 --- /dev/null +++ b/src/MessagePassing/LoopExpansion/graph_clusters.jl @@ -0,0 +1,177 @@ +# Edge-induced no-leaf subgraph enumeration for the loop / linked-cluster expansion. +# +# This file is PURE GRAPH CODE: it never touches a tensor and operates only on a +# NamedGraph. A "generalized loop" is a connected edge set in which EVERY vertex has +# induced degree >= 2 (no leaves). These are exactly the polymers of the hard-core loop +# gas whose log is the free-energy correction `F − ln Z_BP`, and the configurations whose +# antiprojector weights sum to `Z / Z_BP − 1` (see `weight`/`weights` and +# `loopcorrected_free_energy` in `MessagePassing/loopcorrection.jl`). +# +# NamedGraphs' `edgeinduced_subgraphs_no_leaves` builds candidates as unions of simple +# cycles, so it is connected-only AND silently drops any no-leaf subgraph that needs a +# non-cycle (bridge) edge to be no-leaf — e.g. a "dumbbell" of two triangles joined by a +# single bridge edge (the bridge lies on no cycle, so it is never added; the two bare +# triangles fail the connectivity filter). Those diagrams have all degrees >= 2 and a +# nonzero antiprojector weight, so omitting them under-counts the loop series. The +# enumerator here grows connected edge sets directly, so it is connected *by construction* +# and catches the bridge diagrams. +# +# Enumeration uses ESU (Wernicke, "A faster algorithm for detecting network motifs", +# WABI 2006) on the line graph: connected edge sets of `g` are connected vertex sets of the +# line graph, and ESU visits each connected vertex set EXACTLY ONCE via the "smallest-index +# vertex is the root, only grow with strictly-larger exclusive neighbours" rule. That makes +# the de-duplication structural — no `seen` hash set, and no subgraph is ever generated more +# than once — and the per-node bookkeeping (induced degrees, leaf count) is maintained +# incrementally on integer-indexed arrays instead of rebuilt from a dictionary each step. + +# Apply edge `i` to the running induced-degree state, updating the leaf counter +# (`nleaf` = number of vertices whose induced degree is exactly 1). +@inline function _add_edge!(i, esrc, edst, vdeg, nleaf) + @inbounds for v in (esrc[i], edst[i]) + d = vdeg[v] + vdeg[v] = d + 1 + d == 0 ? (nleaf[] += 1) : d == 1 ? (nleaf[] -= 1) : nothing + end + return nothing +end + +# Inverse of `_add_edge!`: remove edge `i` from the running state. +@inline function _del_edge!(i, esrc, edst, vdeg, nleaf) + @inbounds for v in (esrc[i], edst[i]) + d = vdeg[v] + vdeg[v] = d - 1 + d == 1 ? (nleaf[] -= 1) : d == 2 ? (nleaf[] += 1) : nothing + end + return nothing +end + +# ESU recursion over the line graph. `Vsub` is the current edge-index set (a stack), `ext` +# its extension list (candidate edges, all with index > `root`), and `marked[u]` records +# whether edge `u` is already in `Vsub` or has entered the extension somewhere on the path +# from the root (so it is never re-queued — this is what enforces exact-once visiting). +function _esu_extend!( + out, Vsub, ext, root, eadj, esrc, edst, vdeg, nleaf, marked, max_edges, + ) + # `nleaf == 0` with a non-empty edge set <=> every touched vertex has degree >= 2. + (nleaf[] == 0 && !isempty(Vsub)) && push!(out, copy(Vsub)) + length(Vsub) >= max_edges && return + # Leaf-prune: each added edge clears at most 2 degree-1 vertices, so a set with more + # leaves than 2*(remaining budget) can never become leaf-free within budget. + nleaf[] > 2 * (max_edges - length(Vsub)) && return + added = Int[] + while !isempty(ext) + w = pop!(ext) + newext = copy(ext) # siblings still to be processed (ESU: V_ext \ {w}) + empty!(added) + @inbounds for u in eadj[w] # plus the exclusive neighbours of w, index > root + if u > root && !marked[u] + marked[u] = true + push!(newext, u) + push!(added, u) + end + end + _add_edge!(w, esrc, edst, vdeg, nleaf) + push!(Vsub, w) + _esu_extend!(out, Vsub, newext, root, eadj, esrc, edst, vdeg, nleaf, marked, max_edges) + pop!(Vsub) + _del_edge!(w, esrc, edst, vdeg, nleaf) + @inbounds for u in added # unmark only what THIS w introduced; w stays + marked[u] = false # marked for its siblings (parent unmarks it) + end + end + return nothing +end + +""" + connected_edgeinduced_subgraphs_no_leaves(g, max_edges; anchor = nothing) -> Vector + +All connected edge-induced subgraphs of `g` with at most `max_edges` edges in which every +vertex has induced degree `>= 2` (no leaves) — the "generalized loops" of the loop/linked- +cluster series. Unlike `edgeinduced_subgraphs_no_leaves` (which unions simple cycles and so +misses bridge-joined no-leaf subgraphs such as two triangles linked by a single edge), this +grows connected edge sets directly and is therefore both connected-by-construction and +bridge-complete. + +Enumeration is ESU on the line graph: every connected edge set is visited exactly once (no +deduplication hashing), and induced degrees / the leaf count are tracked incrementally. With +`anchor` set, only subgraphs containing that vertex are produced: the anchor-incident edges +are given the smallest indices and used as the ESU roots, so a subgraph is generated from +its minimum edge iff that edge touches the anchor — i.e. iff the subgraph touches the anchor. +With `anchor = nothing` the whole graph is enumerated. +""" +function connected_edgeinduced_subgraphs_no_leaves( + g::AbstractGraph, max_edges::Integer; anchor = nothing, + ) + E0 = collect(edges(g)) + m = length(E0) + (m == 0 || max_edges < 3) && return typeof(g)[] # smallest no-leaf subgraph is a 3-cycle + + # Integer-index the vertices for O(1) degree bookkeeping. + V = collect(vertices(g)) + vidx = Dict{eltype(V), Int}() + for (k, v) in enumerate(V) + vidx[v] = k + end + + # Order the edge list so that, when anchored, the anchor-incident edges occupy the + # smallest indices `1:nseed`; ESU rooted at `1:nseed` then yields exactly the subgraphs + # whose minimum edge touches the anchor == the subgraphs that touch the anchor. + if anchor === nothing + E = E0 + nseed = m + else + ai = vidx[anchor] + anchored_first = sort(eachindex(E0); by = i -> (vidx[src(E0[i])] != ai && vidx[dst(E0[i])] != ai)) + E = E0[anchored_first] + nseed = count(e -> vidx[src(e)] == ai || vidx[dst(e)] == ai, E) + end + + esrc = Vector{Int}(undef, m) + edst = Vector{Int}(undef, m) + inc = [Int[] for _ in 1:length(V)] # incident edge indices per vertex + for (i, e) in enumerate(E) + a = vidx[src(e)] + b = vidx[dst(e)] + esrc[i] = a + edst[i] = b + push!(inc[a], i) + push!(inc[b], i) + end + # Line-graph adjacency: two edges are adjacent iff they share a vertex. + eadj = [Int[] for _ in 1:m] + for ev in inc, i in ev, j in ev + i != j && push!(eadj[i], j) + end + for i in 1:m + unique!(eadj[i]) + end + + vdeg = zeros(Int, length(V)) + nleaf = Ref(0) + marked = falses(m) + Vsub = Int[] + out = Vector{Int}[] + added = Int[] + for r in 1:nseed + _add_edge!(r, esrc, edst, vdeg, nleaf) + push!(Vsub, r) + marked[r] = true + ext = Int[] + empty!(added) + @inbounds for u in eadj[r] + if u > r && !marked[u] + marked[u] = true + push!(ext, u) + push!(added, u) + end + end + _esu_extend!(out, Vsub, ext, r, eadj, esrc, edst, vdeg, nleaf, marked, max_edges) + pop!(Vsub) + _del_edge!(r, esrc, edst, vdeg, nleaf) + marked[r] = false + @inbounds for u in added + marked[u] = false + end + end + return [edge_subgraph(g, E[S]) for S in out] +end diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/LoopExpansion/loopcorrection.jl similarity index 94% rename from src/MessagePassing/loopcorrection.jl rename to src/MessagePassing/LoopExpansion/loopcorrection.jl index e3b82188..793d892f 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/LoopExpansion/loopcorrection.jl @@ -6,10 +6,8 @@ function loopcorrected_partitionfunction( ) zbp = partitionfunction(bp_cache) bp_cache = rescale(bp_cache) - #TODO: Fix edgeinduced_subgraphs_no_leaves for PartitionedGraphView type - #Count the cycles using NamedGraphs egs = - edgeinduced_subgraphs_no_leaves(graph(bp_cache), max_configuration_size) + connected_edgeinduced_subgraphs_no_leaves(graph(bp_cache), max_configuration_size) isempty(egs) && return zbp ws = weights(bp_cache, egs) return zbp * (1 + sum(ws)) diff --git a/src/expect.jl b/src/expect.jl index c3cb0da0..6c77d22b 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -155,6 +155,59 @@ function expect( return expect(alg, ψ_bmps, observable; bmps_messages_up_to_date = true, kwargs...) end +""" + expect(ψ, obs; alg = "loopcorrections", max_configuration_size, ε = 1e-4, cache_update_kwargs...) + +Free-energy / generating-function estimate of a **single-site Hermitian** observable +`obs = (op, vertex[, coeff])`, + + ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_0 ≈ [F(ε) − F(−ε)] / (2ε), + +with `F(t) = ln‖e^{t Ô/2}|ψ⟩‖²` the loop-corrected free energy of that norm network +(`loopcorrected_free_energy`, the additive linked-cluster form `ln Z_BP + Σ_C w_C`; +`max_configuration_size` counts EDGES). Because each `F` is the free energy of a genuine +norm network, every loop cluster is leaf-free and there is no protected operator vertex, so +this typically converges smoothly (little oscillation). BP is re-solved for each shifted +network so the loop corrections include the linear response of the messages. + +`ε` is the central finite-difference step (default `1e-4`). `Ô` must be Hermitian (the +generating-function identity `⟨ψ|e^{εÔ}|ψ⟩ = ‖e^{εÔ/2}|ψ⟩‖²` relies on `Ô† = Ô`). Accepts a +`TensorNetworkState` or a `BeliefPropagationCache`, and a single observable or a vector of +them; multi-site observables are not supported. +""" +function expect( + ::Algorithm"loopcorrections", + ψ::TensorNetworkState, + obs::Tuple; + max_configuration_size::Integer, ε::Real = 1e-4, + cache_update_kwargs = default_bp_update_kwargs(ψ), + ) + op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ)) + iszero(coeff) && return zero(coeff) + length(obs_vs) == 1 || + error("expect(...; alg = \"loopcorrections\") supports single-site observables only; got $(length(obs_vs)) sites.") + v = only(obs_vs) + op_string = only(op_strings) + + Fp = _gated_loop_free_energy(ψ, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) + Fm = _gated_loop_free_energy(ψ, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) + return coeff * (Fp - Fm) / (2ε) +end + +function expect(alg::Algorithm"loopcorrections", cache::BeliefPropagationCache, obs::Tuple; kwargs...) + return expect(alg, network(cache), obs; kwargs...) +end + +function expect( + alg::Algorithm"loopcorrections", + ψ::Union{TensorNetworkState, BeliefPropagationCache}, + observables::Vector{<:Tuple}; + kwargs..., + ) + return map(obs -> expect(alg, ψ, obs; kwargs...), observables) +end + + #Process an observable into more readable form function collectobservable(obs::Tuple, g::NamedGraph) From 257f675a82a42181f8850f41acdd63bc113db3f8 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 12:24:31 -0400 Subject: [PATCH 2/6] Stuff --- .../LoopExpansion/free_energy.jl | 84 ------------------- .../LoopExpansion/loopcorrection.jl | 63 ++++++++++++++ src/TensorNetworkQuantumSimulator.jl | 3 +- 3 files changed, 65 insertions(+), 85 deletions(-) delete mode 100644 src/MessagePassing/LoopExpansion/free_energy.jl diff --git a/src/MessagePassing/LoopExpansion/free_energy.jl b/src/MessagePassing/LoopExpansion/free_energy.jl deleted file mode 100644 index 4d541e5b..00000000 --- a/src/MessagePassing/LoopExpansion/free_energy.jl +++ /dev/null @@ -1,84 +0,0 @@ -# Free-energy (generating-function) estimate of a single-site observable, exposed through -# `expect(...; alg = "loopcorrections")`. -# -# ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_{ε=0} ≈ [F(ε) − F(−ε)] / (2ε), -# F(t) = ln⟨ψ|e^{t Ô}|ψ⟩ = ln‖e^{t Ô/2}|ψ⟩‖² (Hermitian Ô), -# -# Each F is the loop-corrected free energy of a *genuine norm network*: the gate -# e^{±ε Ô/2} is absorbed (un-normalized) into the ket at the observable site and BP is -# re-solved, so the loop expansion prunes ALL leaves — there is no protected operator -# vertex, hence no "anomalous" leaf-containing cluster, which is what tends to make this -# estimator converge smoothly. BP is re-solved on each shifted network so the loop -# corrections also pick up the linear response of the messages. - -# F(α-shifted) = loop-corrected free energy of e^{α Ô}|ψ⟩, built by absorbing the -# (un-normalized) one-site gate into the ket and re-solving BP for the perturbed network. -# `F = ln Z_BP + Σ_C w_C` is the additive linked-cluster free energy (`loopcorrected_free_energy`), -# the genuine extensive log-partition-function whose smooth ε-dependence makes this estimator -# converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which -# only agrees with it to O(w) and resums the same clusters multiplicatively instead. -function _gated_loop_free_energy( - ψ::TensorNetworkState, op_string::String, v, α, max_configuration_size::Integer; - cache_update_kwargs, - ) - s = only(siteinds(ψ)[v]) - G = ITensors.exp(α * ITensors.op(op_string, s); ishermitian = true) - bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) - # `normalize_tensors = false`: the un-normalized gated tensor is exactly what makes the - # squared norm equal the partition function ⟨ψ|e^{2α Ô}|ψ⟩ we want to differentiate. - bpc, _ = apply_gate!(G, bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) - bpc = update(bpc; cache_update_kwargs...) - return loopcorrected_free_energy(bpc, max_configuration_size) -end - -""" - expect(ψ, obs; alg = "loopcorrections", max_configuration_size, ε = 1e-4, cache_update_kwargs...) - -Free-energy / generating-function estimate of a **single-site Hermitian** observable -`obs = (op, vertex[, coeff])`, - - ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_0 ≈ [F(ε) − F(−ε)] / (2ε), - -with `F(t) = ln‖e^{t Ô/2}|ψ⟩‖²` the loop-corrected free energy of that norm network -(`loopcorrected_free_energy`, the additive linked-cluster form `ln Z_BP + Σ_C w_C`; -`max_configuration_size` counts EDGES). Because each `F` is the free energy of a genuine -norm network, every loop cluster is leaf-free and there is no protected operator vertex, so -this typically converges smoothly (little oscillation). BP is re-solved for each shifted -network so the loop corrections include the linear response of the messages. - -`ε` is the central finite-difference step (default `1e-4`). `Ô` must be Hermitian (the -generating-function identity `⟨ψ|e^{εÔ}|ψ⟩ = ‖e^{εÔ/2}|ψ⟩‖²` relies on `Ô† = Ô`). Accepts a -`TensorNetworkState` or a `BeliefPropagationCache`, and a single observable or a vector of -them; multi-site observables are not supported. -""" -function expect( - ::Algorithm"loopcorrections", - ψ::TensorNetworkState, - obs::Tuple; - max_configuration_size::Integer, ε::Real = 1e-4, - cache_update_kwargs = default_bp_update_kwargs(ψ), - ) - op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ)) - iszero(coeff) && return zero(coeff) - length(obs_vs) == 1 || - error("expect(...; alg = \"loopcorrections\") supports single-site observables only; got $(length(obs_vs)) sites.") - v = only(obs_vs) - op_string = only(op_strings) - - Fp = _gated_loop_free_energy(ψ, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) - Fm = _gated_loop_free_energy(ψ, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) - return coeff * (Fp - Fm) / (2ε) -end - -function expect(alg::Algorithm"loopcorrections", cache::BeliefPropagationCache, obs::Tuple; kwargs...) - return expect(alg, network(cache), obs; kwargs...) -end - -function expect( - alg::Algorithm"loopcorrections", - ψ::Union{TensorNetworkState, BeliefPropagationCache}, - observables::Vector{<:Tuple}; - kwargs..., - ) - return map(obs -> expect(alg, ψ, obs; kwargs...), observables) -end diff --git a/src/MessagePassing/LoopExpansion/loopcorrection.jl b/src/MessagePassing/LoopExpansion/loopcorrection.jl index 793d892f..bc14e5c9 100644 --- a/src/MessagePassing/LoopExpansion/loopcorrection.jl +++ b/src/MessagePassing/LoopExpansion/loopcorrection.jl @@ -13,6 +13,69 @@ function loopcorrected_partitionfunction( return zbp * (1 + sum(ws)) end +# Linked-cluster (free-energy) form of the loop correction. +# +# F = ln Z = ln Z_BP + ln(Z / Z_BP), Z / Z_BP = Σ_{polymer configs} ∏ w_C , +# +# where the polymers are the *connected* no-leaf edge subgraphs (generalized loops) and the +# hard-core exclusion is vertex-sharing. The linked-cluster theorem makes F extensive: its +# expansion is a sum over CONNECTED clusters only — the vertex-disjoint products that the +# partition-function series `loopcorrected_partitionfunction` would have to enumerate beyond +# total size `2·girth − 1` are resummed by the exponential and never appear here. To leading +# cumulant order this is `ln Z_BP + Σ_C w_C` over the connected no-leaf clusters (overlapping- +# cluster corrections, `−½ Σ_{C∼C'} w_C w_{C'} + …`, are higher order). Both this and +# `loopcorrected_partitionfunction` enumerate the same generalized loops via +# `connected_edgeinduced_subgraphs_no_leaves` (so both capture the bridge "dumbbell" +# diagrams); they differ only in how the cluster weights are resummed. +# +# Note `loopcorrected_free_energy` and `log(loopcorrected_partitionfunction)` differ at +# O(w²): they agree exactly only when no loops fit (Σw = 0), where both reduce to `ln Z_BP`. +function loopcorrected_free_energy( + bp_cache::BeliefPropagationCache, + max_configuration_size::Integer, + ) + zbp = partitionfunction(bp_cache) + F = log(complex(zbp)) + bp_cache = rescale(bp_cache) + g = graph(bp_cache) + egs = connected_edgeinduced_subgraphs_no_leaves(g, max_configuration_size) + isempty(egs) && return F + return F + sum(weights(bp_cache, egs)) +end + +# Free-energy (generating-function) estimate of a single-site observable, exposed through +# `expect(...; alg = "loopcorrections")`. +# +# ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_{ε=0} ≈ [F(ε) − F(−ε)] / (2ε), +# F(t) = ln⟨ψ|e^{t Ô}|ψ⟩ = ln‖e^{t Ô/2}|ψ⟩‖² (Hermitian Ô), +# +# Each F is the loop-corrected free energy of a *genuine norm network*: the gate +# e^{±ε Ô/2} is absorbed (un-normalized) into the ket at the observable site and BP is +# re-solved, so the loop expansion prunes ALL leaves — there is no protected operator +# vertex, hence no "anomalous" leaf-containing cluster, which is what tends to make this +# estimator converge smoothly. BP is re-solved on each shifted network so the loop +# corrections also pick up the linear response of the messages. + +# F(α-shifted) = loop-corrected free energy of e^{α Ô}|ψ⟩, built by absorbing the +# (un-normalized) one-site gate into the ket and re-solving BP for the perturbed network. +# `F = ln Z_BP + Σ_C w_C` is the additive linked-cluster free energy (`loopcorrected_free_energy`), +# the genuine extensive log-partition-function whose smooth ε-dependence makes this estimator +# converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which +# only agrees with it to O(w) and resums the same clusters multiplicatively instead. +function _gated_loop_free_energy( + ψ::TensorNetworkState, op_string::String, v, α, max_configuration_size::Integer; + cache_update_kwargs, + ) + s = only(siteinds(ψ)[v]) + G = ITensors.exp(α * ITensors.op(op_string, s); ishermitian = true) + bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) + # `normalize_tensors = false`: the un-normalized gated tensor is exactly what makes the + # squared norm equal the partition function ⟨ψ|e^{2α Ô}|ψ⟩ we want to differentiate. + bpc, _ = apply_gate!(G, bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) + bpc = update(bpc; cache_update_kwargs...) + return loopcorrected_free_energy(bpc, max_configuration_size) +end + #Transform the indices in the given subgraph of the tensornetwork so that antiprojectors can be inserted without duplicate indices appearing function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) bpc = copy(bpc) diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 6d0c752b..1fc05706 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -16,7 +16,8 @@ include("Forms/quadraticform.jl") include("MessagePassing/abstractbeliefpropagationcache.jl") include("MessagePassing/beliefpropagationcache.jl") include("MessagePassing/boundarympscache.jl") -include("MessagePassing/loopcorrection.jl") +include("MessagePassing/LoopExpansion/graph_clusters.jl") +include("MessagePassing/LoopExpansion/loopcorrection.jl") include("graph_ops.jl") include("utils.jl") From 0bfabf8f139b87bb570bce63b958bcb22068c312 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 12:50:55 -0400 Subject: [PATCH 3/6] Stuff --- examples/loopcorrections.jl | 15 ++++++++++++++- src/Apply/apply_gates.jl | 5 +++++ .../LoopExpansion/loopcorrection.jl | 10 +++++----- src/expect.jl | 18 ++++++++++-------- src/imports.jl | 1 + src/utils.jl | 2 +- 6 files changed, 36 insertions(+), 15 deletions(-) diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index 8b599a94..7c137032 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -19,8 +19,11 @@ function main() (named_grid((nx, ny)), "square", 4), ] for (g, g_str, smallest_loop_size) in gs + println("\n") + println("-----------------------") + obs = ("Z", first(center(g))) println("Testing for $g_str lattice with $(nv(g)) vertices") - ψ = random_tensornetworkstate(ComplexF32, g, "S=1/2"; bond_dimension = χ) + ψ = random_tensornetworkstate(Float64, g, "S=1/2"; bond_dimension = χ) ψ = normalize(ψ; alg = "bp") @@ -31,6 +34,16 @@ function main() println("Bp Value for norm is $norm_bp") println("1st Order Loop Corrected Value for norm is $norm_loopcorrected") println("Exact Value for norm is $norm_exact") + + sz_bp = expect(ψ, obs; alg = "bp") + sz_loopcorrected = expect(ψ, obs; alg = "loopcorrections",max_configuration_size = 2 * (smallest_loop_size) - 1) + sz_exact = expect(ψ, obs; alg = "exact") + + println("\n") + + println("Bp Value for sz is $sz_bp") + println("1st Order Loop Corrected Value for sz is $sz_loopcorrected") + println("Exact Value for sz is $sz_exact") end return end diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index e9d7beeb..a7f0270f 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -127,4 +127,9 @@ function apply_gate!( return ψ_bpc, err end +function apply_gate(gate::ITensor, ψ_bpc::BeliefPropagationCache; kwargs...) + ψ_bpc = copy(ψ_bpc) + return apply_gate!(gate, ψ_bpc; kwargs...) +end + const apply_circuit = apply_gates diff --git a/src/MessagePassing/LoopExpansion/loopcorrection.jl b/src/MessagePassing/LoopExpansion/loopcorrection.jl index bc14e5c9..7de88c06 100644 --- a/src/MessagePassing/LoopExpansion/loopcorrection.jl +++ b/src/MessagePassing/LoopExpansion/loopcorrection.jl @@ -63,17 +63,17 @@ end # converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which # only agrees with it to O(w) and resums the same clusters multiplicatively instead. function _gated_loop_free_energy( - ψ::TensorNetworkState, op_string::String, v, α, max_configuration_size::Integer; + ψ_bpc::BeliefPropagationCache, op_string::String, v, α, max_configuration_size::Integer; cache_update_kwargs, ) + ψ = network(ψ_bpc) s = only(siteinds(ψ)[v]) G = ITensors.exp(α * ITensors.op(op_string, s); ishermitian = true) - bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) # `normalize_tensors = false`: the un-normalized gated tensor is exactly what makes the # squared norm equal the partition function ⟨ψ|e^{2α Ô}|ψ⟩ we want to differentiate. - bpc, _ = apply_gate!(G, bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) - bpc = update(bpc; cache_update_kwargs...) - return loopcorrected_free_energy(bpc, max_configuration_size) + ψ_bpc, _ = apply_gate(G, ψ_bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) + ψ_bpc = update(ψ_bpc; cache_update_kwargs...) + return loopcorrected_free_energy(ψ_bpc, max_configuration_size) end #Transform the indices in the given subgraph of the tensornetwork so that antiprojectors can be inserted without duplicate indices appearing diff --git a/src/expect.jl b/src/expect.jl index 6c77d22b..a65fba9b 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -176,26 +176,28 @@ generating-function identity `⟨ψ|e^{εÔ}|ψ⟩ = ‖e^{εÔ/2}|ψ⟩‖²` r them; multi-site observables are not supported. """ function expect( - ::Algorithm"loopcorrections", - ψ::TensorNetworkState, + alg::Algorithm"loopcorrections", + ψ_bpc::BeliefPropagationCache, obs::Tuple; max_configuration_size::Integer, ε::Real = 1e-4, - cache_update_kwargs = default_bp_update_kwargs(ψ), + cache_update_kwargs = default_bp_update_kwargs(ψ_bpc), ) - op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ)) + op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ_bpc)) iszero(coeff) && return zero(coeff) length(obs_vs) == 1 || error("expect(...; alg = \"loopcorrections\") supports single-site observables only; got $(length(obs_vs)) sites.") v = only(obs_vs) op_string = only(op_strings) - Fp = _gated_loop_free_energy(ψ, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) - Fm = _gated_loop_free_energy(ψ, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) + Fp = _gated_loop_free_energy(ψ_bpc, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) + Fm = _gated_loop_free_energy(ψ_bpc, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) return coeff * (Fp - Fm) / (2ε) end -function expect(alg::Algorithm"loopcorrections", cache::BeliefPropagationCache, obs::Tuple; kwargs...) - return expect(alg, network(cache), obs; kwargs...) +function expect(alg::Algorithm"loopcorrections", + ψ::TensorNetworkState, obs::Tuple; cache_update_kwargs= default_bp_update_kwargs(ψ), kwargs...) + ψ_bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) + return expect(alg, ψ_bpc, obs; kwargs...) end function expect( diff --git a/src/imports.jl b/src/imports.jl index 910796db..044c1ead 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -31,6 +31,7 @@ using NamedGraphs.GraphsExtensions: a_star, add_edge!, edgetype, + edge_subgraph, leaf_vertices, post_order_dfs_edges, decorate_graph_edges, diff --git a/src/utils.jl b/src/utils.jl index 58844c5e..0327c1a0 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -56,7 +56,7 @@ function algorithm_check(tns::Union{AbstractBeliefPropagationCache, TensorNetwor return error("Expected BeliefPropagationCache or TensorNetworkState for 'loop correction' algorithm, got $(typeof(tns))") end - if f ∈ ["normalize", "expect", "sample", "truncate", "rdm"] + if f ∈ ["normalize", "sample", "truncate", "rdm"] return error("Loop correction-based contraction not supported for this functionality yet") end elseif alg == "boundarymps" From 32fbf2dd447186061796b93b4a164d8a7450f939 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 12:58:34 -0400 Subject: [PATCH 4/6] Bump --- Project.toml | 2 +- examples/loopcorrections.jl | 2 +- src/MessagePassing/LoopExpansion/loopcorrection.jl | 2 +- src/expect.jl | 8 +++----- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index d8915558..8b6b3c90 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" -version = "0.3.9" +version = "0.3.10" authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbitrary topology." diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index 7c137032..adadc536 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -23,7 +23,7 @@ function main() println("-----------------------") obs = ("Z", first(center(g))) println("Testing for $g_str lattice with $(nv(g)) vertices") - ψ = random_tensornetworkstate(Float64, g, "S=1/2"; bond_dimension = χ) + ψ = random_tensornetworkstate(ComplexF64, g, "S=1/2"; bond_dimension = χ) ψ = normalize(ψ; alg = "bp") diff --git a/src/MessagePassing/LoopExpansion/loopcorrection.jl b/src/MessagePassing/LoopExpansion/loopcorrection.jl index 7de88c06..3969b86b 100644 --- a/src/MessagePassing/LoopExpansion/loopcorrection.jl +++ b/src/MessagePassing/LoopExpansion/loopcorrection.jl @@ -62,7 +62,7 @@ end # the genuine extensive log-partition-function whose smooth ε-dependence makes this estimator # converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which # only agrees with it to O(w) and resums the same clusters multiplicatively instead. -function _gated_loop_free_energy( +function gated_lc_free_energy( ψ_bpc::BeliefPropagationCache, op_string::String, v, α, max_configuration_size::Integer; cache_update_kwargs, ) diff --git a/src/expect.jl b/src/expect.jl index a65fba9b..b39893c2 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -165,9 +165,7 @@ Free-energy / generating-function estimate of a **single-site Hermitian** observ with `F(t) = ln‖e^{t Ô/2}|ψ⟩‖²` the loop-corrected free energy of that norm network (`loopcorrected_free_energy`, the additive linked-cluster form `ln Z_BP + Σ_C w_C`; -`max_configuration_size` counts EDGES). Because each `F` is the free energy of a genuine -norm network, every loop cluster is leaf-free and there is no protected operator vertex, so -this typically converges smoothly (little oscillation). BP is re-solved for each shifted +`max_configuration_size` counts EDGES). BP is re-solved for each shifted network so the loop corrections include the linear response of the messages. `ε` is the central finite-difference step (default `1e-4`). `Ô` must be Hermitian (the @@ -189,8 +187,8 @@ function expect( v = only(obs_vs) op_string = only(op_strings) - Fp = _gated_loop_free_energy(ψ_bpc, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) - Fm = _gated_loop_free_energy(ψ_bpc, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) + Fp = gated_lc_free_energy(ψ_bpc, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) + Fm = gated_lc_free_energy(ψ_bpc, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) return coeff * (Fp - Fm) / (2ε) end From e3fd28655bfd100c58102f92179b7023f3eb04e1 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 14:08:40 -0400 Subject: [PATCH 5/6] Stuff --- examples/loopcorrections.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index adadc536..945a2d41 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -10,7 +10,7 @@ using Random Random.seed!(1634) function main() - nx, ny = 4, 4 + nx, ny = 4,4 χ = 3 ITensors.disable_warn_order() gs = [ @@ -19,6 +19,7 @@ function main() (named_grid((nx, ny)), "square", 4), ] for (g, g_str, smallest_loop_size) in gs + max_configuration_size = 2*smallest_loop_size -1 println("\n") println("-----------------------") obs = ("Z", first(center(g))) @@ -28,7 +29,7 @@ function main() ψ = normalize(ψ; alg = "bp") norm_bp = norm(ψ; alg = "bp") - norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size = 2 * (smallest_loop_size) - 1) + norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size) norm_exact = norm(ψ; alg = "exact") println("Bp Value for norm is $norm_bp") @@ -36,7 +37,7 @@ function main() println("Exact Value for norm is $norm_exact") sz_bp = expect(ψ, obs; alg = "bp") - sz_loopcorrected = expect(ψ, obs; alg = "loopcorrections",max_configuration_size = 2 * (smallest_loop_size) - 1) + sz_loopcorrected = expect(ψ, obs; alg = "loopcorrections",max_configuration_size) sz_exact = expect(ψ, obs; alg = "exact") println("\n") From a1eec94918ae684ec2fd28315f018dcec6edf367 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Jun 2026 15:30:37 -0400 Subject: [PATCH 6/6] Docs --- docs/src/api.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/src/api.md b/docs/src/api.md index 1f3ba7bc..ee8e9ca6 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -83,6 +83,12 @@ update update_iteration! ``` +## Loop Cluster Expansion + +```@docs +connected_edgeinduced_subgraphs_no_leaves +``` + ## Utilities ```@docs