From cd93a8b49edbdf41c74d7e0aecb1fed8fef6c409 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 20 Jul 2026 15:06:55 -0400 Subject: [PATCH 1/3] Support fermionic sites in simple update and belief propagation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The simple-update gauge reads each environment's ket and bra legs from the state tensor and takes the Hermitian square root in the transposed bipartition, which carries the fermionic braid sign on the odd-parity sector. Doubled ket/bra messages normalize by their trace rather than an entrywise sum, which is sign-correct on fermionic bonds. Operator construction routes through project_aux, so a parity-odd operator gets an auxiliary charge leg instead of throwing, and project_pair builds operator pairs that share a contractible aux leg for fermion strings like c†c. --- src/Apply/simple_update.jl | 15 ++-- .../abstractbeliefpropagationcache.jl | 11 ++- src/Ops.jl | 9 ++- src/itensors.jl | 68 ++++++++++++++++--- 4 files changed, 84 insertions(+), 19 deletions(-) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 05f15cd..6bae780 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -37,12 +37,15 @@ function simple_update( # The environments are hermitian only up to numerical noise, so project before # the square roots (which require hermitian input). - sqrt_invsqrt = env -> sqrth_invsqrth_safe( - project_hermitian(env, (inds(env)[1],), (inds(env)[2],)), - (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0 - ) - sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1) - sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2) + sqrt_invsqrt = function (env, ψᵥ) + ket, bra = commonind(env, ψᵥ), uniqueind(env, ψᵥ) + return sqrth_invsqrth_safe( + project_hermitian(env, (ket,), (bra,)), + (bra,), (ket,); atol = sqrt_cutoff, rtol = 0 + ) + end + sqrt_inv_sqrt_envs_v1 = map(env -> sqrt_invsqrt(env, ψ⃗[1]), envs_v1) + sqrt_inv_sqrt_envs_v2 = map(env -> sqrt_invsqrt(env, ψ⃗[2]), envs_v2) sqrt_envs_v1, inv_sqrt_envs_v1 = first.(sqrt_inv_sqrt_envs_v1), last.(sqrt_inv_sqrt_envs_v1) sqrt_envs_v2, inv_sqrt_envs_v2 = first.(sqrt_inv_sqrt_envs_v2), last.(sqrt_inv_sqrt_envs_v2) diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index dfb8f04..4202cec 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -180,7 +180,16 @@ function updated_message( updated_message = contract_network(contract_list; sequence) if alg.kwargs.normalize - message_norm = sum(updated_message) + codomain, domain = operator_inds(updated_message) + # A doubled (ket/bra) network message is a bond operator: normalize by its + # trace, which is sign-correct for fermionic bonds (the entrywise `sum` can + # flip the message's sign). A single-layer network message is a vector with + # no bra/ket pairing and no trace, so fall back to the entrywise sum there. + message_norm = if all(in(inds(updated_message)), codomain) + tr(updated_message, codomain, domain) + else + sum(updated_message) + end if !iszero(message_norm) updated_message = updated_message / message_norm end diff --git a/src/Ops.jl b/src/Ops.jl index 1d2a67c..8f6e724 100644 --- a/src/Ops.jl +++ b/src/Ops.jl @@ -59,14 +59,17 @@ end # Embed a `d^n × d^n` operator matrix (computational basis, first site most # significant) into an `ITensor` with codomain `prime.(sites)` (outputs) and -# domain `sites` (inputs). `project` is checked, so on graded sites an operator -# that is not symmetric under the site index's grading throws an `InexactError`. +# domain `sites` (inputs). Routed through `project_aux` (like `state`): a graded +# operator that is odd under the site grading (e.g. `X`/`Y`, bare `c`/`c†`) has no +# parity-block-diagonal `(prime.(sites), sites)` map, so it gets a trailing auxiliary +# charge leg instead of throwing an `InexactError`. Even operators and dense (ungraded) +# sites project cleanly, so they are unaffected (no aux leg). function _op_matrix_to_itensor(M::AbstractMatrix, sites::Tuple) ds = length.(sites) n = length(sites) A = reshape(Matrix{ComplexF64}(M), (reverse(ds)..., reverse(ds)...)) A = permutedims(A, (reverse(1:n)..., reverse((n + 1):(2n))...)) - return project(A, ITensorBase.prime.(sites), sites) + return project_aux(A, ITensorBase.prime.(sites), sites) end # Top-level `op(name, sites...; kwargs...)`. The identity is dimension-general (used diff --git a/src/itensors.jl b/src/itensors.jl index c31d5de..44fcedd 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -6,16 +6,66 @@ using Adapt: Adapt using ITensorBase: ITensorBase, ITensor, Index, inds, noprime, plev, prime, space, unnamed using TensorAlgebra: TensorAlgebra, project, tryproject -function project_aux(v::AbstractVector{<:Number}, i::Index) - length(v) == length(i) || - error( - "state vector has dimension $(length(v)) but the site index has dimension $(length(i))" +# Project `a` onto the symmetry-restricted space given by `codomain`/`domain`, mirroring +# `TensorAlgebra.project`/`tryproject` (three-argument operator form; `a` is indexed positionally +# as `(codomain..., domain...)`). When the plain projection is charge-forbidden — a parity-odd +# state or operator — the residual charge is absorbed into a trailing auxiliary index instead of +# throwing, so the returned ITensor carries that extra aux leg last. +function project_aux(a::AbstractArray, codomain, domain) + return @something tryproject(a, codomain, domain) begin + raw = project(reshape(a, (size(a)..., 1)), space.(codomain), space.(domain)) + aux = Index(TensorAlgebra.axes(raw, ndims(a) + 1)) + ITensor(raw, (codomain..., domain..., aux)) + end +end +project_aux(a::AbstractArray, codomain) = project_aux(a, codomain, ()) +project_aux(v::AbstractVector{<:Number}, i::Index) = project_aux(v, (i,)) + +# Build a pair of tensors that share one contractible auxiliary leg, so `t1 * t2` contracts the +# pair (e.g. the fermion string of `c†ᵢcⱼ` from a `c†` and a `c`) with no `flip`. `codomain`/`domain` +# give each operand's own axes; the shared aux is appended to each `domain`. Unlike `project_aux`, we +# always want the aux, so we go straight to `project` (no `tryproject` fallback). Two layers mirror +# `project`'s own split: a name-agnostic core, and a named wrapper over it. +# +# Name-agnostic core: mint the aux by projecting `a1` with a trailing dummy axis — `project` absorbs +# `a1`'s residual charge into it — then attach that same aux axis to `a2`. Returns the two raw arrays, +# each carrying the aux as its trailing axis. Uses only `project`/`reshape`/`axes` (no names), so it +# could live beside `TensorAlgebra.project`. +function project_pair(a1::AbstractArray, codomain1, domain1, a2::AbstractArray, codomain2, domain2) + p1 = project(reshape(a1, (size(a1)..., 1)), codomain1, domain1) + aux = TensorAlgebra.axes(p1, ndims(a1) + 1) + p2 = project(reshape(a2, (size(a2)..., 1)), codomain2, (domain2..., aux)) + return p1, p2 +end + +# Named worker: strip to bare axes, call the core, and reattach names — the one extra step vs +# `project`'s wrapper is minting a single fresh aux `Index` shared by both outputs, so `t1 * t2` +# contracts them by name. +function project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2) + p1, p2 = project_pair( + a1, space.(codomain1), space.(domain1), + a2, space.(codomain2), space.(domain2), + ) + aux = Index(TensorAlgebra.axes(p1, ndims(a1) + 1)) + t1 = ITensor(p1, (codomain1..., domain1..., aux)) + t2 = ITensor(p2, (codomain2..., domain2..., aux)) + return t1, t2 +end + +# Named layer (mirrors ITensorBase's named `project`): two entries read operand 1's flavor from +# whichever of its sides is non-empty (an empty `codomain1` is the all-domain/state mirror of the +# empty-domain form), so a named operand never falls through to the bare-axis core. +function project_pair( + a1::AbstractArray, codomain1::Tuple{Index, Vararg{Index}}, domain1::Tuple{Vararg{Index}}, + a2::AbstractArray, codomain2::Tuple{Vararg{Index}}, domain2::Tuple{Vararg{Index}} + ) + return project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2) +end +function project_pair( + a1::AbstractArray, codomain1::Tuple{}, domain1::Tuple{Index, Vararg{Index}}, + a2::AbstractArray, codomain2::Tuple{Vararg{Index}}, domain2::Tuple{Vararg{Index}} ) - ψ = tryproject(v, (i,)) - isnothing(ψ) || return ψ - raw = project(reshape(v, (length(v), 1)), (space(i),), ()) - aux = Index(TensorAlgebra.axes(raw, 2)) - return ITensor(raw, (i, aux)) + return project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2) end function onehot(eltype::Type, (i, p)::Pair{<:Index}) From 4b42755ece834e44207b54cc1ed2769b3f9ff049 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 20 Jul 2026 15:29:22 -0400 Subject: [PATCH 2/3] Clean up the fermion changes after review Rename local index variables (ket_ind/bra_ind, codomain_inds/domain_inds, projected_a), factor an is_operator predicate out of the belief-propagation message normalization instead of inlining the primed-index check, and tighten comments. --- src/Apply/simple_update.jl | 8 ++++-- .../abstractbeliefpropagationcache.jl | 12 ++++---- src/itensors.jl | 28 ++++++++++--------- 3 files changed, 25 insertions(+), 23 deletions(-) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 6bae780..b9ce458 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -38,10 +38,12 @@ function simple_update( # The environments are hermitian only up to numerical noise, so project before # the square roots (which require hermitian input). sqrt_invsqrt = function (env, ψᵥ) - ket, bra = commonind(env, ψᵥ), uniqueind(env, ψᵥ) + ket_ind, bra_ind = commonind(env, ψᵥ), uniqueind(env, ψᵥ) + # Factorize from bra to ket (codomain = bra), the bipartition in which the + # fermionic message is PSD. return sqrth_invsqrth_safe( - project_hermitian(env, (ket,), (bra,)), - (bra,), (ket,); atol = sqrt_cutoff, rtol = 0 + project_hermitian(env, (ket_ind,), (bra_ind,)), + (bra_ind,), (ket_ind,); atol = sqrt_cutoff, rtol = 0 ) end sqrt_inv_sqrt_envs_v1 = map(env -> sqrt_invsqrt(env, ψ⃗[1]), envs_v1) diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 4202cec..e16fcfc 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -180,13 +180,11 @@ function updated_message( updated_message = contract_network(contract_list; sequence) if alg.kwargs.normalize - codomain, domain = operator_inds(updated_message) - # A doubled (ket/bra) network message is a bond operator: normalize by its - # trace, which is sign-correct for fermionic bonds (the entrywise `sum` can - # flip the message's sign). A single-layer network message is a vector with - # no bra/ket pairing and no trace, so fall back to the entrywise sum there. - message_norm = if all(in(inds(updated_message)), codomain) - tr(updated_message, codomain, domain) + # A doubled (ket/bra) network message is a bond operator: normalize by its trace, + # which is sign-correct for fermionic bonds (the entrywise `sum` can flip the sign). + # A single-layer network message is a vector with no bra/ket pairing, so use `sum`. + message_norm = if is_operator(updated_message) + tr(updated_message, operator_inds(updated_message)...) else sum(updated_message) end diff --git a/src/itensors.jl b/src/itensors.jl index 44fcedd..1b9c8da 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -13,24 +13,18 @@ using TensorAlgebra: TensorAlgebra, project, tryproject # throwing, so the returned ITensor carries that extra aux leg last. function project_aux(a::AbstractArray, codomain, domain) return @something tryproject(a, codomain, domain) begin - raw = project(reshape(a, (size(a)..., 1)), space.(codomain), space.(domain)) - aux = Index(TensorAlgebra.axes(raw, ndims(a) + 1)) - ITensor(raw, (codomain..., domain..., aux)) + projected_a = project(reshape(a, (size(a)..., 1)), space.(codomain), space.(domain)) + aux = Index(TensorAlgebra.axes(projected_a, ndims(a) + 1)) + ITensor(projected_a, (codomain..., domain..., aux)) end end project_aux(a::AbstractArray, codomain) = project_aux(a, codomain, ()) project_aux(v::AbstractVector{<:Number}, i::Index) = project_aux(v, (i,)) -# Build a pair of tensors that share one contractible auxiliary leg, so `t1 * t2` contracts the -# pair (e.g. the fermion string of `c†ᵢcⱼ` from a `c†` and a `c`) with no `flip`. `codomain`/`domain` -# give each operand's own axes; the shared aux is appended to each `domain`. Unlike `project_aux`, we -# always want the aux, so we go straight to `project` (no `tryproject` fallback). Two layers mirror -# `project`'s own split: a name-agnostic core, and a named wrapper over it. -# -# Name-agnostic core: mint the aux by projecting `a1` with a trailing dummy axis — `project` absorbs -# `a1`'s residual charge into it — then attach that same aux axis to `a2`. Returns the two raw arrays, -# each carrying the aux as its trailing axis. Uses only `project`/`reshape`/`axes` (no names), so it -# could live beside `TensorAlgebra.project`. +# Build two tensors that share one contractible auxiliary leg, so `t1 * t2` gives the fermion +# string of e.g. `c†ᵢcⱼ` (a `c†` and a `c`) with no `flip`. The aux is minted by projecting `a1` +# with a trailing dummy axis, which absorbs `a1`'s residual charge, then attaching it to `a2`. +# This bare-axis core takes spaces (no names), so it could live beside `TensorAlgebra.project`. function project_pair(a1::AbstractArray, codomain1, domain1, a2::AbstractArray, codomain2, domain2) p1 = project(reshape(a1, (size(a1)..., 1)), codomain1, domain1) aux = TensorAlgebra.axes(p1, ndims(a1) + 1) @@ -115,6 +109,14 @@ delta(is::Tuple) = delta(Float64, is) delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) +# Whether `a` is an operator tensor: its indices are exactly the plev-0 indices paired with +# their primes. A plain state/vector has the plev-0 indices but not their primed partners. +function is_operator(a::ITensor) + domain = filter(i -> plev(i) == 0, inds(a)) + codomain = prime.(domain) + return issetequal(inds(a), [codomain; domain]) +end + # The codomain/domain bipartition of an operator tensor: each plev-0 index paired with its # prime. Viewing the operator as this square map is what `tr` factors through. function operator_inds(a::ITensor) From bc2ee80bc2dba44db4324f0e6935fcdf321954db Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 21 Jul 2026 14:33:41 -0400 Subject: [PATCH 3/3] Raise the TensorAlgebra compat floor to 0.17.6 The fermionic simple-update and belief-propagation path relies on the corrected fermionic operator contraction from GradedArrays 0.14, which requires TensorAlgebra 0.17.6. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9e7d1ca..8cc9f16 100644 --- a/Project.toml +++ b/Project.toml @@ -41,7 +41,7 @@ SimpleGraphConverter = "0.1.0" SplitApplyCombine = "1.2.3" Statistics = "1.11.1" StatsBase = "0.34.4" -TensorAlgebra = "0.17.5" +TensorAlgebra = "0.17.6" TypeParameterAccessors = "0.3.10, 0.4" VectorInterface = "0.5" julia = "1.10"