From 49619851a17c877e0aac231dd730615a6bf4e60e Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 20 Jul 2026 11:56:02 -0400
Subject: [PATCH 01/13] Add fermionic support to belief-propagation simple
 update

Keep belief-propagation messages as operators and normalize them by their trace rather than an entrywise sum, which stays sign correct on fermionic bonds where the entrywise sum can flip the odd-parity block. Gauge the simple-update tensors with the balanced square root and inverse square root of the Hermitian-projected bond messages.

Depends on the operator-form Hermitian factorizations added in the ITensorBase mf/operator-sqrth branch.
---
 src/apply/apply_operators.jl               | 36 +++++++++++++---------
 src/beliefpropagation/beliefpropagation.jl | 27 +++++++++++++---
 2 files changed, 45 insertions(+), 18 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index a885ca9..60a37a8 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -5,9 +5,9 @@ using Graphs: dst, src, vertices
 using ITensorBase:
     ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator, replacedimnames
 using LinearAlgebra: norm
-using MatrixAlgebraKit: qr_compact, svd_trunc
+using MatrixAlgebraKit: project_hermitian, qr_compact, svd_trunc
 using NamedGraphs.GraphsExtensions: all_edges, boundary_edges
-using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv
+using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe
 
 # === Top-level user entry point ===
 
@@ -205,6 +205,13 @@ end
 
 # === BP simple-update implementation ===
 
+# BP simple-update gauge. Each message is a bond operator that is positive semidefinite in
+# its intrinsic bra/ket bipartition — its domain (ket) leg is the one shared with the state
+# tensor it gauges. The balanced sqrt / inverse-sqrt of the Hermitian-projected message is
+# the gauge; the fermionic braid sign is carried by the graded contraction, so no bipartition
+# flip is needed. This works for bosonic (ungraded) messages too, generalizing the previous
+# `gram_eigh_full` gauge.
+
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
         state::AbstractITensorNetwork, env; kwargs...
@@ -233,7 +240,7 @@ function apply_gate_bp_nsite!(
     ψv = ITB.apply(op, state[v])
     if normalize
         gauges = [
-            conj(gram_eigh_full(env[e]))
+            sqrth_safe(project_hermitian(env[e]))
                 for e in boundary_edges(state, vs; dir = :in)
         ]
         ψv /= norm(prod([[ψv]; gauges]))
@@ -249,12 +256,16 @@ function apply_gate_bp_nsite!(
     )
     v1, v2 = vs
     edges_in = boundary_edges(state, vs; dir = :in)
-    grams_v1 =
-        [gram_eigh_full_with_pinv(env[e]) for e in edges_in if dst(e) == v1]
-    grams_v2 =
-        [gram_eigh_full_with_pinv(env[e]) for e in edges_in if dst(e) == v2]
-    gauges_v1, inv_gauges_v1 = conj.(first.(grams_v1)), conj.(last.(grams_v1))
-    gauges_v2, inv_gauges_v2 = conj.(first.(grams_v2)), conj.(last.(grams_v2))
+    sqrts_invsqrts_v1 = [
+        sqrth_invsqrth_safe(project_hermitian(env[e]))
+            for e in edges_in if dst(e) == v1
+    ]
+    sqrts_invsqrts_v2 = [
+        sqrth_invsqrth_safe(project_hermitian(env[e]))
+            for e in edges_in if dst(e) == v2
+    ]
+    gauges_v1, inv_gauges_v1 = first.(sqrts_invsqrts_v1), conj.(last.(sqrts_invsqrts_v1))
+    gauges_v2, inv_gauges_v2 = first.(sqrts_invsqrts_v2), conj.(last.(sqrts_invsqrts_v2))
 
     ψ_v1 = prod([[state[v1]]; gauges_v1])
     ψ_v2 = prod([[state[v2]]; gauges_v2])
@@ -267,7 +278,7 @@ function apply_gate_bp_nsite!(
         S = S / norm(S)
     end
     name_v1, name_v2 = dimnames(S)
-    sqrt_S = sqrt(S, (name_v1,), (name_v2,))
+    sqrt_S = sqrth_safe(S, (name_v1,), (name_v2,); atol = 0, rtol = 0)
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
@@ -275,9 +286,6 @@ function apply_gate_bp_nsite!(
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
     env[v1 => v2] = operator(conj(S), (name_v2,), (name_v1,))
-    env[v2 => v1] = operator(
-        conj(replacedimnames(S, name_v1 => name_v2, name_v2 => name_v1)),
-        (name_v2,), (name_v1,)
-    )
+    env[v2 => v1] = operator(S, (name_v1,), (name_v2,))
     return dest
 end
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index bc987cd..cb80129 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -3,8 +3,9 @@ using .AlgorithmsInterfaceExtensions:
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using DataGraphs: edge_data
 using Graphs: AbstractEdge, edges, edgetype, has_edge, vertices
-using ITensorBase: AbstractITensor
-using LinearAlgebra: norm, normalize
+using ITensorBase:
+    AbstractITensor, NamedTensorOperator, codomainnames, domainnames, operator, state
+using LinearAlgebra: norm, normalize, tr
 using NamedGraphs.GraphsExtensions:
     add_edges!, boundary_edges, forest_cover_edge_sequence, subgraph
 using NamedGraphs.PartitionedGraphs: quotientvertices
@@ -243,10 +244,28 @@ function message_update!(algorithm::SimpleMessageUpdate, cache, factors, edge)
     messages = collect(incoming_messages(cache, edge))
     factor = factors[src(edge)]
 
-    new_message = contract_network([messages; [factor]]; alg = algorithm.contraction_alg)
+    # `contract_network` works on plain named arrays, so unwrap any operator messages to
+    # their underlying tensors before contracting (fermionic signs ride on the graded
+    # arrays, so nothing is lost).
+    message_tensors = map(m -> m isa NamedTensorOperator ? state(m) : m, messages)
+    new_message = contract_network(
+        [message_tensors; [factor]]; alg = algorithm.contraction_alg
+    )
+
+    # `contract_network` drops the bra/ket operator structure, so restore it from the
+    # existing message. A doubled (ket/bra) message is then a bond operator: normalize by
+    # its trace, which is sign-correct for fermionic bonds (the entrywise `sum` can flip
+    # the odd-parity block's sign). A single-layer message stays a vector with no bra/ket
+    # pairing, so fall back to the entrywise sum there.
+    old_message = cache[edge]
+    if old_message isa NamedTensorOperator
+        new_message =
+            operator(new_message, codomainnames(old_message), domainnames(old_message))
+    end
 
     if algorithm.normalize
-        message_norm = sum(new_message)
+        message_norm =
+            new_message isa NamedTensorOperator ? tr(new_message) : sum(new_message)
         if !iszero(message_norm)
             new_message /= message_norm
         end

From bd8c20ee8a38248028a2a7c4dddbe164bb4d3422 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 21 Jul 2026 16:02:42 -0400
Subject: [PATCH 02/13] Contract operator-valued gauges as plain tensors in BP
 simple update

The BP simple-update gauge absorbs an operator-valued square root of each incoming message into the state tensor. Contracting it with `*` returns an operator, so the updated state is written back as an operator and rejected by the state network. Take the underlying tensor of each gauge (`ITB.state`) before contracting, so the gauged state stays a plain tensor while the index arrows still carry the fermion sign. Requires the operator `dimnametype` fix from https://github.com/ITensor/ITensorBase.jl/pull/225.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 60a37a8..b6715df 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -205,12 +205,8 @@ end
 
 # === BP simple-update implementation ===
 
-# BP simple-update gauge. Each message is a bond operator that is positive semidefinite in
-# its intrinsic bra/ket bipartition — its domain (ket) leg is the one shared with the state
-# tensor it gauges. The balanced sqrt / inverse-sqrt of the Hermitian-projected message is
-# the gauge; the fermionic braid sign is carried by the graded contraction, so no bipartition
-# flip is needed. This works for bosonic (ungraded) messages too, generalizing the previous
-# `gram_eigh_full` gauge.
+# The message's domain (ket) leg is the one shared with the state tensor it gauges, so the
+# fermionic braid sign is carried by the graded contraction and no bipartition flip is needed.
 
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
@@ -243,7 +239,7 @@ function apply_gate_bp_nsite!(
             sqrth_safe(project_hermitian(env[e]))
                 for e in boundary_edges(state, vs; dir = :in)
         ]
-        ψv /= norm(prod([[ψv]; gauges]))
+        ψv /= norm(prod([[ψv]; ITB.state.(gauges)]))
     end
     dest[v] = ψv
     return dest
@@ -267,8 +263,8 @@ function apply_gate_bp_nsite!(
     gauges_v1, inv_gauges_v1 = first.(sqrts_invsqrts_v1), conj.(last.(sqrts_invsqrts_v1))
     gauges_v2, inv_gauges_v2 = first.(sqrts_invsqrts_v2), conj.(last.(sqrts_invsqrts_v2))
 
-    ψ_v1 = prod([[state[v1]]; gauges_v1])
-    ψ_v2 = prod([[state[v2]]; gauges_v2])
+    ψ_v1 = prod([[state[v1]]; ITB.state.(gauges_v1)])
+    ψ_v2 = prod([[state[v2]]; ITB.state.(gauges_v2)])
 
     Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(dimnames(ψ_v1), dimnames(ψ_v2), dimnames(op)))
     Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(dimnames(ψ_v2), dimnames(ψ_v1), dimnames(op)))
@@ -282,8 +278,8 @@ function apply_gate_bp_nsite!(
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
-    dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
-    dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
+    dest[v1] = prod([[Q_v1 * R_v1]; ITB.state.(inv_gauges_v1)])
+    dest[v2] = prod([[Q_v2 * R_v2]; ITB.state.(inv_gauges_v2)])
 
     env[v1 => v2] = operator(conj(S), (name_v2,), (name_v1,))
     env[v2 => v1] = operator(S, (name_v1,), (name_v2,))

From a3dd19588722dd020135745d3d8563495fffdd1a Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 21 Jul 2026 20:46:20 -0400
Subject: [PATCH 03/13] Use uniform BP messages with a PSD-bipartition gauge

Store both directed messages of a bond uniformly as the bond (ket output) and
an auxiliary leg (bra input), with conj(S)/S into the two endpoints, and build
the initial identity messages on the source-side ket axis so their ket leg is
dual to the destination bond. The gauge projects the message Hermitian and takes
the balanced square-root / inverse-square-root in whichever bipartition is
positive semidefinite, which recovers the odd-parity fermion sign in both
directions. Migrate the message-update path to the ITensorBase 0.13
outputnames/inputnames API and bump to 0.9.10.
---
 Project.toml                               |  2 +-
 src/apply/apply_operators.jl               | 61 ++++++++++++++--------
 src/beliefpropagation/beliefpropagation.jl |  4 +-
 src/beliefpropagation/messagecache.jl      | 12 +++--
 4 files changed, 49 insertions(+), 30 deletions(-)

diff --git a/Project.toml b/Project.toml
index 00bf323..f9e3661 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
-version = "0.9.9"
+version = "0.9.10"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index b6715df..4912bcb 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -2,8 +2,8 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, src, vertices
-using ITensorBase:
-    ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator, replacedimnames
+using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
+    outputnames, replacedimnames
 using LinearAlgebra: norm
 using MatrixAlgebraKit: project_hermitian, qr_compact, svd_trunc
 using NamedGraphs.GraphsExtensions: all_edges, boundary_edges
@@ -205,8 +205,26 @@ end
 
 # === BP simple-update implementation ===
 
-# The message's domain (ket) leg is the one shared with the state tensor it gauges, so the
-# fermionic braid sign is carried by the graded contraction and no bipartition flip is needed.
+# Balanced square-root / inverse-square-root gauge from an incoming BP message. Following the
+# TNQS `simple_update` recipe, project the message to Hermitian in its (ket, bra) = (output,
+# input) bipartition, then take the balanced √ / inv-√ in the bipartition where the projected
+# message is positive semidefinite. For a fermionic bond the odd-parity sign lands on one of the
+# two bipartitions depending on the bond arrows: the two directed messages of a bond carry
+# opposite ket arrows (each is dual to its own endpoint, and the two endpoints' bond arrows are
+# opposite), so they are PSD in opposite bipartitions — one in the transposed (bra, ket), the
+# other in the intrinsic (ket, bra). Taking the sqrt in whichever is PSD recovers the odd-parity
+# fermion sign in both directions. `ITB.state` unwraps the operator to the underlying tensor, so
+# the returned gauges are plain tensors contracted directly into the state.
+function message_gauge(message)
+    ket, bra = outputnames(message), inputnames(message)
+    hermitian_message = project_hermitian(ITB.state(message), ket, bra)
+    return try
+        sqrth_invsqrth_safe(hermitian_message, bra, ket)
+    catch err
+        err isa DomainError || rethrow()
+        sqrth_invsqrth_safe(hermitian_message, ket, bra)
+    end
+end
 
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
@@ -236,10 +254,10 @@ function apply_gate_bp_nsite!(
     ψv = ITB.apply(op, state[v])
     if normalize
         gauges = [
-            sqrth_safe(project_hermitian(env[e]))
+            first(message_gauge(env[e]))
                 for e in boundary_edges(state, vs; dir = :in)
         ]
-        ψv /= norm(prod([[ψv]; ITB.state.(gauges)]))
+        ψv /= norm(prod([[ψv]; gauges]))
     end
     dest[v] = ψv
     return dest
@@ -252,19 +270,13 @@ function apply_gate_bp_nsite!(
     )
     v1, v2 = vs
     edges_in = boundary_edges(state, vs; dir = :in)
-    sqrts_invsqrts_v1 = [
-        sqrth_invsqrth_safe(project_hermitian(env[e]))
-            for e in edges_in if dst(e) == v1
-    ]
-    sqrts_invsqrts_v2 = [
-        sqrth_invsqrth_safe(project_hermitian(env[e]))
-            for e in edges_in if dst(e) == v2
-    ]
-    gauges_v1, inv_gauges_v1 = first.(sqrts_invsqrts_v1), conj.(last.(sqrts_invsqrts_v1))
-    gauges_v2, inv_gauges_v2 = first.(sqrts_invsqrts_v2), conj.(last.(sqrts_invsqrts_v2))
-
-    ψ_v1 = prod([[state[v1]]; ITB.state.(gauges_v1)])
-    ψ_v2 = prod([[state[v2]]; ITB.state.(gauges_v2)])
+    siv_v1 = [message_gauge(env[e]) for e in edges_in if dst(e) == v1]
+    siv_v2 = [message_gauge(env[e]) for e in edges_in if dst(e) == v2]
+    gauges_v1, inv_gauges_v1 = first.(siv_v1), conj.(last.(siv_v1))
+    gauges_v2, inv_gauges_v2 = first.(siv_v2), conj.(last.(siv_v2))
+
+    ψ_v1 = prod([[state[v1]]; gauges_v1])
+    ψ_v2 = prod([[state[v2]]; gauges_v2])
 
     Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(dimnames(ψ_v1), dimnames(ψ_v2), dimnames(op)))
     Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(dimnames(ψ_v2), dimnames(ψ_v1), dimnames(op)))
@@ -278,10 +290,15 @@ function apply_gate_bp_nsite!(
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
-    dest[v1] = prod([[Q_v1 * R_v1]; ITB.state.(inv_gauges_v1)])
-    dest[v2] = prod([[Q_v2 * R_v2]; ITB.state.(inv_gauges_v2)])
+    dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
+    dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    env[v1 => v2] = operator(conj(S), (name_v2,), (name_v1,))
+    # Uniform messages over the bond (ket) and the auxiliary leg (bra), mirroring TNQS's
+    # `s_values` / `conj(s_values)` but keeping ITNN's random names (`name_v2`) in place of TNQS's
+    # `prime(u)`: `conj(S)` goes into `v2` and `S` into `v1`. The two directions carry opposite
+    # ket arrows, so `conj` flips the odd-parity sector between them; `message_gauge` then recovers
+    # the fermion sign by factorizing each in its PSD bipartition.
+    env[v1 => v2] = operator(conj(S), (name_v1,), (name_v2,))
     env[v2 => v1] = operator(S, (name_v1,), (name_v2,))
     return dest
 end
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index cb80129..1bc653c 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -4,7 +4,7 @@ using AlgorithmsInterface: AlgorithmsInterface as AI
 using DataGraphs: edge_data
 using Graphs: AbstractEdge, edges, edgetype, has_edge, vertices
 using ITensorBase:
-    AbstractITensor, NamedTensorOperator, codomainnames, domainnames, operator, state
+    AbstractITensor, NamedTensorOperator, inputnames, operator, outputnames, state
 using LinearAlgebra: norm, normalize, tr
 using NamedGraphs.GraphsExtensions:
     add_edges!, boundary_edges, forest_cover_edge_sequence, subgraph
@@ -260,7 +260,7 @@ function message_update!(algorithm::SimpleMessageUpdate, cache, factors, edge)
     old_message = cache[edge]
     if old_message isa NamedTensorOperator
         new_message =
-            operator(new_message, codomainnames(old_message), domainnames(old_message))
+            operator(new_message, outputnames(old_message), inputnames(old_message))
     end
 
     if algorithm.normalize
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index b428b16..4c17e2c 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -209,13 +209,15 @@ function similar_message_environment(nn::NormNetwork)
             ketview = KetView(nn)
 
             ketnames = linknames(ketview, edge)
+            ketaxis = unnamed.(linkaxes(ketview, edge))
 
-            brainds = linkinds(braview, edge)
-            branames = name.(brainds)
-            braaxis = unnamed.(brainds)
+            branames = name.(linkinds(braview, edge))
 
-            # Message axis is conj to the tensor it points to.
-            message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
+            # Uniform message convention: the bond leg shared with the ket tensor is the operator
+            # output (ket), the bra-layer leg is the input (bra). The message is built on the
+            # src-side ket axis, whose arrow is opposite the dst endpoint's bond, so the gauge's
+            # ket leg contracts back into the destination state.
+            message = similar_operator(ketview[vertex], ketaxis, ketnames, branames)
 
             return edge => message
         end

From 249eb127b4bc5cb9062f5cced8c1425aa8791091 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 22 Jul 2026 11:01:40 -0400
Subject: [PATCH 04/13] Derive the reverse BP message by a sign-free duality
 flip

Replace the try/catch PSD-bipartition selection in `message_gauge` with a
plain transposed gauge. The two directed bond messages are now built from the
PSD singular-value message `conj(S)`. The reverse direction is its arrow flip
via `dualize`, which toggles `isdual` and keeps the odd-parity block, rather
than `conj`, which negates the odd block and breaks PSD. Free-fermion
validation is exact to machine precision on trees with no fallback.

`dualize` is a temporary element-wise hack that pulls in a temporary
`GradedArrays` dependency. Both stand in for a proper backend-agnostic
duality-flip primitive.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 Project.toml                          |  1 +
 src/apply/apply_operators.jl          | 48 +++++++++++++++------------
 src/beliefpropagation/messagecache.jl |  7 ++--
 3 files changed, 30 insertions(+), 26 deletions(-)

diff --git a/Project.toml b/Project.toml
index f9e3661..8dbf36d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ AlgorithmsInterface = "d1e3940c-cd12-4505-8585-b0a4b322527d"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
+GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 4912bcb..013e1fa 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -1,6 +1,7 @@
 using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
+using GradedArrays: GradedArrays
 using Graphs: dst, src, vertices
 using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
     outputnames, replacedimnames
@@ -205,25 +206,30 @@ end
 
 # === BP simple-update implementation ===
 
-# Balanced square-root / inverse-square-root gauge from an incoming BP message. Following the
-# TNQS `simple_update` recipe, project the message to Hermitian in its (ket, bra) = (output,
-# input) bipartition, then take the balanced √ / inv-√ in the bipartition where the projected
-# message is positive semidefinite. For a fermionic bond the odd-parity sign lands on one of the
-# two bipartitions depending on the bond arrows: the two directed messages of a bond carry
-# opposite ket arrows (each is dual to its own endpoint, and the two endpoints' bond arrows are
-# opposite), so they are PSD in opposite bipartitions — one in the transposed (bra, ket), the
-# other in the intrinsic (ket, bra). Taking the sqrt in whichever is PSD recovers the odd-parity
-# fermion sign in both directions. `ITB.state` unwraps the operator to the underlying tensor, so
-# the returned gauges are plain tensors contracted directly into the state.
+# Balanced √ / inv-√ gauge from a BP message: project it Hermitian in its (ket, bra) = (output,
+# input) bipartition, then take the root in the bipartition where the projection is positive
+# semidefinite. The two directed messages of a fermionic bond have opposite ket arrows, so the
+# odd-parity sign makes each PSD in only one bipartition — try the transposed (bra, ket), fall
+# back to the intrinsic (ket, bra). `ITB.state` unwraps to the underlying tensor so the gauges
+# contract straight into the state.
 function message_gauge(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)
-    return try
-        sqrth_invsqrth_safe(hermitian_message, bra, ket)
-    catch err
-        err isa DomainError || rethrow()
-        sqrth_invsqrth_safe(hermitian_message, ket, bra)
+    return sqrth_invsqrth_safe(hermitian_message, bra, ket)
+end
+
+# HACK (experiment): reinterpret a graded tensor on dual-flipped axes (`isdual` toggled, sector
+# labels kept) WITHOUT the fermionic braid sign that `conj` applies. Flips the bond arrows while
+# keeping the odd-parity block intact.
+function dualize(t)
+    a = ITB.unnamed(t)
+    da = similar(a, eltype(a), map(GradedArrays.dual, axes(a)))
+    fill!(da, zero(eltype(a)))
+    for I in CartesianIndices(size(a))
+        v = a[I]
+        iszero(v) || (da[I] = v)
     end
+    return ITB.nameddims(da, dimnames(t))
 end
 
 function apply_gate_bp!(
@@ -293,12 +299,10 @@ function apply_gate_bp_nsite!(
     dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    # Uniform messages over the bond (ket) and the auxiliary leg (bra), mirroring TNQS's
-    # `s_values` / `conj(s_values)` but keeping ITNN's random names (`name_v2`) in place of TNQS's
-    # `prime(u)`: `conj(S)` goes into `v2` and `S` into `v1`. The two directions carry opposite
-    # ket arrows, so `conj` flips the odd-parity sector between them; `message_gauge` then recovers
-    # the fermion sign by factorizing each in its PSD bipartition.
-    env[v1 => v2] = operator(conj(S), (name_v1,), (name_v2,))
-    env[v2 => v1] = operator(S, (name_v1,), (name_v2,))
+    # `conj(S)` is the PSD bond message; the reverse direction is its arrow-flip that keeps the
+    # odd-parity block (via `dualize`, not `conj`, which would negate it and break PSD).
+    psd_message = conj(S)
+    env[v1 => v2] = operator(psd_message, (name_v1,), (name_v2,))
+    env[v2 => v1] = operator(dualize(psd_message), (name_v1,), (name_v2,))
     return dest
 end
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 4c17e2c..61f25ec 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -213,10 +213,9 @@ function similar_message_environment(nn::NormNetwork)
 
             branames = name.(linkinds(braview, edge))
 
-            # Uniform message convention: the bond leg shared with the ket tensor is the operator
-            # output (ket), the bra-layer leg is the input (bra). The message is built on the
-            # src-side ket axis, whose arrow is opposite the dst endpoint's bond, so the gauge's
-            # ket leg contracts back into the destination state.
+            # Bond leg (ket) = operator output, bra-layer leg = input. Built on the src-side ket
+            # axis, whose arrow is opposite the dst endpoint's bond, so the gauge contracts back
+            # into the destination state.
             message = similar_operator(ketview[vertex], ketaxis, ketnames, branames)
 
             return edge => message

From 45444b52d430d61d975603e6f6feffb0d70f8d1f Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 22 Jul 2026 12:47:43 -0400
Subject: [PATCH 05/13] Use a ribbon twist for the reverse fermionic bond
 message

The reverse directed message is S with its bond orientation reversed,
which is the ribbon twist of S (-1 on the odd-parity sector). Both
directed messages are then positive semidefinite in the transposed
gauge, so the gauge needs no fallback. This replaces the element-wise
dualize helper with a name-keyed twist and twist!.

Also simplify the belief-propagation message unwrap to state.(messages),
relying on state being idempotent.
---
 src/apply/apply_operators.jl               | 26 ++++++++--------------
 src/beliefpropagation/beliefpropagation.jl |  2 +-
 2 files changed, 10 insertions(+), 18 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 013e1fa..d24dde1 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -218,19 +218,12 @@ function message_gauge(message)
     return sqrth_invsqrth_safe(hermitian_message, bra, ket)
 end
 
-# HACK (experiment): reinterpret a graded tensor on dual-flipped axes (`isdual` toggled, sector
-# labels kept) WITHOUT the fermionic braid sign that `conj` applies. Flips the bond arrows while
-# keeping the odd-parity block intact.
-function dualize(t)
-    a = ITB.unnamed(t)
-    da = similar(a, eltype(a), map(GradedArrays.dual, axes(a)))
-    fill!(da, zero(eltype(a)))
-    for I in CartesianIndices(size(a))
-        v = a[I]
-        iszero(v) || (da[I] = v)
-    end
-    return ITB.nameddims(da, dimnames(t))
+# TODO: replace with an ITensorBase-level `twist` to drop the direct `GradedArrays` dependency.
+function twist!(t, names)
+    GradedArrays.twist!(ITB.unnamed(t), map(n -> findfirst(==(n), dimnames(t)), names))
+    return t
 end
+twist(t, names) = twist!(copy(t), names)
 
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
@@ -299,10 +292,9 @@ function apply_gate_bp_nsite!(
     dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    # `conj(S)` is the PSD bond message; the reverse direction is its arrow-flip that keeps the
-    # odd-parity block (via `dualize`, not `conj`, which would negate it and break PSD).
-    psd_message = conj(S)
-    env[v1 => v2] = operator(psd_message, (name_v1,), (name_v2,))
-    env[v2 => v1] = operator(dualize(psd_message), (name_v1,), (name_v2,))
+    # The two directed messages are `conj(S)` and the ribbon twist of `S` over its ket side,
+    # both positive semidefinite in the transposed gauge.
+    env[v1 => v2] = operator(conj(S), (name_v1,), (name_v2,))
+    env[v2 => v1] = operator(twist(S, (name_v1,)), (name_v1,), (name_v2,))
     return dest
 end
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index 1bc653c..6255df1 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -247,7 +247,7 @@ function message_update!(algorithm::SimpleMessageUpdate, cache, factors, edge)
     # `contract_network` works on plain named arrays, so unwrap any operator messages to
     # their underlying tensors before contracting (fermionic signs ride on the graded
     # arrays, so nothing is lost).
-    message_tensors = map(m -> m isa NamedTensorOperator ? state(m) : m, messages)
+    message_tensors = state.(messages)
     new_message = contract_network(
         [message_tensors; [factor]]; alg = algorithm.contraction_alg
     )

From d0311e520147f0daca565ae6883227ba6337ce70 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 22 Jul 2026 18:43:44 -0400
Subject: [PATCH 06/13] Use linknames for the bra link names in
 similar_message_environment

---
 src/beliefpropagation/messagecache.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 61f25ec..22d3aa4 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -211,7 +211,7 @@ function similar_message_environment(nn::NormNetwork)
             ketnames = linknames(ketview, edge)
             ketaxis = unnamed.(linkaxes(ketview, edge))
 
-            branames = name.(linkinds(braview, edge))
+            branames = linknames(braview, edge)
 
             # Bond leg (ket) = operator output, bra-layer leg = input. Built on the src-side ket
             # axis, whose arrow is opposite the dst endpoint's bond, so the gauge contracts back

From 8d169d4e78fcf22d81047fbd3913c7f35ce2c88b Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 13:33:09 -0400
Subject: [PATCH 07/13] Build fermionic BP messages by doubled contraction
 instead of a ribbon twist

Each directed simple-update message is now the doubled contraction of the gate
factor with its conjugate, rather than being read off the singular-value spectrum
as conj(S) and a ribbon-twisted S. The graded contraction of a factor with its
conjugate carries the odd-parity fermion sign on its own, so the explicit ribbon
twist (and the direct GradedArrays dependency it needed) is gone.
---
 src/apply/apply_operators.jl | 27 +++++++++------------------
 1 file changed, 9 insertions(+), 18 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index d24dde1..28b15b5 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -1,7 +1,6 @@
 using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
-using GradedArrays: GradedArrays
 using Graphs: dst, src, vertices
 using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
     outputnames, replacedimnames
@@ -206,25 +205,14 @@ end
 
 # === BP simple-update implementation ===
 
-# Balanced √ / inv-√ gauge from a BP message: project it Hermitian in its (ket, bra) = (output,
-# input) bipartition, then take the root in the bipartition where the projection is positive
-# semidefinite. The two directed messages of a fermionic bond have opposite ket arrows, so the
-# odd-parity sign makes each PSD in only one bipartition — try the transposed (bra, ket), fall
-# back to the intrinsic (ket, bra). `ITB.state` unwraps to the underlying tensor so the gauges
-# contract straight into the state.
+# The odd-parity sign leaves a fermionic message positive semidefinite in only one
+# bipartition, so root it in the transposed (bra, ket) one.
 function message_gauge(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)
     return sqrth_invsqrth_safe(hermitian_message, bra, ket)
 end
 
-# TODO: replace with an ITensorBase-level `twist` to drop the direct `GradedArrays` dependency.
-function twist!(t, names)
-    GradedArrays.twist!(ITB.unnamed(t), map(n -> findfirst(==(n), dimnames(t)), names))
-    return t
-end
-twist(t, names) = twist!(copy(t), names)
-
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
         state::AbstractITensorNetwork, env; kwargs...
@@ -292,9 +280,12 @@ function apply_gate_bp_nsite!(
     dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    # The two directed messages are `conj(S)` and the ribbon twist of `S` over its ket side,
-    # both positive semidefinite in the transposed gauge.
-    env[v1 => v2] = operator(conj(S), (name_v1,), (name_v2,))
-    env[v2 => v1] = operator(twist(S, (name_v1,)), (name_v1,), (name_v2,))
+    # The graded contraction of a factor with its conjugate carries the odd-parity sign.
+    env[v1 => v2] = operator(
+        replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v1,), (name_v2,)
+    )
+    env[v2 => v1] = operator(
+        replacedimnames(conj(R_v2), name_v1 => name_v2) * R_v2, (name_v1,), (name_v2,)
+    )
     return dest
 end

From f4863ac401deb965dd979b34cf85086040171b24 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 14:44:25 -0400
Subject: [PATCH 08/13] Inline the eigendecomposition in the BP message gauge
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Build the message gauge and its inverse explicitly from `eigh_full` as
`v * √d * v'` and `v * √d⁻¹ * v'`, rather than calling `sqrth_invsqrth_safe`,
so the mechanism handling the odd-parity bipartition is visible in place.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 28b15b5..9d74f6b 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -5,9 +5,9 @@ using Graphs: dst, src, vertices
 using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
     outputnames, replacedimnames
 using LinearAlgebra: norm
-using MatrixAlgebraKit: project_hermitian, qr_compact, svd_trunc
+using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
 using NamedGraphs.GraphsExtensions: all_edges, boundary_edges
-using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe
+using TensorAlgebra.MatrixAlgebra: invsqrth_safe, sqrth_safe
 
 # === Top-level user entry point ===
 
@@ -206,11 +206,19 @@ end
 # === BP simple-update implementation ===
 
 # The odd-parity sign leaves a fermionic message positive semidefinite in only one
-# bipartition, so root it in the transposed (bra, ket) one.
+# bipartition, so diagonalize it in the transposed (bra, ket) one. From the eigenvectors
+# `v` and eigenvalues `d`, the gauge and its inverse are the balanced roots
+# `v * √d * v'` and `v * √d⁻¹ * v'`, with `v` supplying the ket leg and its `bra`-relabeled
+# copy the other, so the odd-parity sign stays inside the graded contraction.
 function message_gauge(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)
-    return sqrth_invsqrth_safe(hermitian_message, bra, ket)
+    d, v = eigh_full(hermitian_message, bra, ket)
+    name_d′, name_d = dimnames(d)
+    v_ket = conj(v)
+    v_bra = replacedimnames(v, only(ket) => only(bra), name_d => name_d′)
+    return v_bra * sqrth_safe(d, (name_d′,), (name_d,)) * v_ket,
+        v_bra * invsqrth_safe(d, (name_d′,), (name_d,)) * v_ket
 end
 
 function apply_gate_bp!(

From 7a29271afa4e794f4374047be9ba170c338240ae Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 14:50:56 -0400
Subject: [PATCH 09/13] Drop the unused GradedArrays dependency

`src/` no longer references GradedArrays directly, so remove it from the
package dependencies. It remains available transitively and in the test
project, where it is still used.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 8dbf36d..f9e3661 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,7 +12,6 @@ AlgorithmsInterface = "d1e3940c-cd12-4505-8585-b0a4b322527d"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
-GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 53a8d2b50d00f8658e42aed43e113586bbebed98 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 14:59:23 -0400
Subject: [PATCH 10/13] Trim the message-gauge comment

Drop the narration of the leg relabeling and the unclear claim about where
the sign lands, keeping only the reason for the transposed bipartition.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 9d74f6b..3ebbe28 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -206,10 +206,8 @@ end
 # === BP simple-update implementation ===
 
 # The odd-parity sign leaves a fermionic message positive semidefinite in only one
-# bipartition, so diagonalize it in the transposed (bra, ket) one. From the eigenvectors
-# `v` and eigenvalues `d`, the gauge and its inverse are the balanced roots
-# `v * √d * v'` and `v * √d⁻¹ * v'`, with `v` supplying the ket leg and its `bra`-relabeled
-# copy the other, so the odd-parity sign stays inside the graded contraction.
+# bipartition, so diagonalize it in the transposed (bra, ket) one. The gauge and its inverse
+# are then the balanced roots `v * √d * v'` and `v * √d⁻¹ * v'` of the eigendecomposition.
 function message_gauge(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)

From 3e3b6f9a047b875977e6126fe45de71001eb44e7 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 14:59:58 -0400
Subject: [PATCH 11/13] Note the gauge choice and the open asymmetric-gauge
 issue
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Record in the comment that `v * √d * v'` is used because it works with
fermions, and that the asymmetric `√d * v'` gauge has an issue still under
investigation.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 3ebbe28..fea15e3 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -206,8 +206,9 @@ end
 # === BP simple-update implementation ===
 
 # The odd-parity sign leaves a fermionic message positive semidefinite in only one
-# bipartition, so diagonalize it in the transposed (bra, ket) one. The gauge and its inverse
-# are then the balanced roots `v * √d * v'` and `v * √d⁻¹ * v'` of the eigendecomposition.
+# bipartition, so diagonalize it in the transposed (bra, ket) one. We use the balanced gauge
+# `v * √d * v'` (and inverse `v * √d⁻¹ * v'`), which works with fermions; the asymmetric gauge
+# `√d * v'` has an issue that is under investigation.
 function message_gauge(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)

From 5ee86e23d0f79a508b392b87c48ac5e84ec0bc46 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 17:36:49 -0400
Subject: [PATCH 12/13] Split out message_root for the single-site gauge

The single-site normalization only needs the message square root, not its
inverse, so `first(message_gauge(...))` built and discarded the inverse root.
Factor the shared eigendecomposition into `message_eigen` and add
`message_root`, which the single-site path now calls. `message_gauge` still
returns the root together with its inverse for the two-site split.
---
 src/apply/apply_operators.jl | 25 ++++++++++++++++++++-----
 1 file changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index fea15e3..3c4fb23 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -206,16 +206,31 @@ end
 # === BP simple-update implementation ===
 
 # The odd-parity sign leaves a fermionic message positive semidefinite in only one
-# bipartition, so diagonalize it in the transposed (bra, ket) one. We use the balanced gauge
-# `v * √d * v'` (and inverse `v * √d⁻¹ * v'`), which works with fermions; the asymmetric gauge
-# `√d * v'` has an issue that is under investigation.
-function message_gauge(message)
+# bipartition, so diagonalize it in the transposed (bra, ket) one, placing the
+# eigenvectors on the bra and ket legs ready to sandwich a function of the eigenvalues.
+# Shared by `message_root` and `message_gauge`.
+function message_eigen(message)
     ket, bra = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), ket, bra)
     d, v = eigh_full(hermitian_message, bra, ket)
     name_d′, name_d = dimnames(d)
     v_ket = conj(v)
     v_bra = replacedimnames(v, only(ket) => only(bra), name_d => name_d′)
+    return d, v_bra, v_ket
+end
+
+# The balanced Hermitian root `v * √d * v'` of the message. This gauge works with
+# fermions; the asymmetric gauge `√d * v'` has an issue that is under investigation.
+function message_root(message)
+    d, v_bra, v_ket = message_eigen(message)
+    name_d′, name_d = dimnames(d)
+    return v_bra * sqrth_safe(d, (name_d′,), (name_d,)) * v_ket
+end
+
+# The message root paired with its inverse `v * √d⁻¹ * v'`, to gauge a bond and undo it.
+function message_gauge(message)
+    d, v_bra, v_ket = message_eigen(message)
+    name_d′, name_d = dimnames(d)
     return v_bra * sqrth_safe(d, (name_d′,), (name_d,)) * v_ket,
         v_bra * invsqrth_safe(d, (name_d′,), (name_d,)) * v_ket
 end
@@ -248,7 +263,7 @@ function apply_gate_bp_nsite!(
     ψv = ITB.apply(op, state[v])
     if normalize
         gauges = [
-            first(message_gauge(env[e]))
+            message_root(env[e])
                 for e in boundary_edges(state, vs; dir = :in)
         ]
         ψv /= norm(prod([[ψv]; gauges]))

From 8701d0e78dea7fa850ee7ee14691275beec47b64 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 24 Jul 2026 17:36:53 -0400
Subject: [PATCH 13/13] Test the fermionic apply-gate over a carried BP
 environment

Cover `FermionParity` sites alongside the bosonic ranges. The apply tests now
carry the belief-propagation environment through the gate applications, as real
usage does, instead of rebuilding it from `beliefpropagation` and re-wrapping
the messages with a test helper. Belief-propagation convergence is covered by
`test_beliefpropagation.jl`. This needs GradedArrays 0.14, whose graded
square-root path is sign-correct for fermions where 0.13 throws.
---
 test/Project.toml           |  4 ++-
 test/test_apply_operator.jl | 70 ++++++++-----------------------------
 2 files changed, 18 insertions(+), 56 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index 0711602..a49dd6d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -19,6 +19,7 @@ SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [sources.ITensorNetworksNext]
@@ -29,7 +30,7 @@ AlgorithmsInterface = "0.1"
 Aqua = "0.8.14"
 DataGraphs = "0.5"
 Dictionaries = "0.4.5"
-GradedArrays = "0.13.2"
+GradedArrays = "0.14"
 Graphs = "1.13.1"
 ITensorBase = "0.13"
 ITensorNetworksNext = "0.9"
@@ -44,4 +45,5 @@ SparseArraysBase = "0.10.4"
 StableRNGs = "1"
 Suppressor = "0.2.8"
 TensorAlgebra = "0.16, 0.17"
+TensorKitSectors = "0.3"
 Test = "1.10"
diff --git a/test/test_apply_operator.jl b/test/test_apply_operator.jl
index aa65609..de9aa93 100644
--- a/test/test_apply_operator.jl
+++ b/test/test_apply_operator.jl
@@ -1,19 +1,18 @@
 using GradedArrays: U1, gradedrange
 using Graphs: dst, edges, src, vertices
-using ITensorBase: ITensorBase as ITB, Index, name, named, operator, replacedimnames,
-    setname, state, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, apply_operator,
-    apply_operators, beliefpropagation, braname, insertlink!, linknames,
-    message_environment, messagecache, tensornetwork
+using ITensorBase: ITensorBase as ITB, Index, name, operator, setname, uniquename
+using ITensorNetworksNext: NormNetwork, apply_operator, apply_operators, insertlink!,
+    message_environment, tensornetwork
 using MatrixAlgebraKit: svd_trunc, truncrank
 using NamedGraphs.NamedGraphGenerators: named_cycle_graph, named_path_graph
-using NamedGraphs: NamedGraph
 using Random: AbstractRNG
 using StableRNGs: StableRNG
+using TensorKitSectors: FermionParity
 using Test: @test, @testset
 
 const spinone = Base.OneTo(3)
 const spinone_u1 = gradedrange([U1(2) => 1, U1(0) => 1, U1(-2) => 1])
+const fermion = gradedrange([FermionParity(0) => 2, FermionParity(1) => 2])
 
 function randn_operator(rng::AbstractRNG, elt::Type, domain_namedaxes)
     codomain_namedaxes = setname.(domain_namedaxes, uniquename.(name.(domain_namedaxes)))
@@ -22,6 +21,10 @@ function randn_operator(rng::AbstractRNG, elt::Type, domain_namedaxes)
     return operator(data, name.(codomain_namedaxes), name.(domain_namedaxes))
 end
 
+# Build a random state by applying random gates layer by layer, carrying the belief
+# propagation environment through the applications. The returned `env` is the environment
+# the gate applications produced, ready to gauge the next application (belief-propagation
+# convergence itself is covered separately in `test_beliefpropagation.jl`).
 function random_state(rng::AbstractRNG, elt::Type, g, site_axes; nlayers, trunc)
     network = tensornetwork(vertices(g)) do v
         return randn(rng, elt, (site_axes[v],))
@@ -36,24 +39,11 @@ function random_state(rng::AbstractRNG, elt::Type, g, site_axes; nlayers, trunc)
         gate = randn_operator(rng, elt, (site_axes[src(e)], site_axes[dst(e)]))
         network, env = apply_operator(gate, network, env; trunc)
     end
-    return network
+    return network, env
 end
 
-function operator_message_cache(nn::NormNetwork, messages)
-    return messagecache(keys(messages)) do edge
-        ketnames = linknames(KetView(nn), edge)
-        branames = linknames(BraView(nn), edge)
-
-        bramap = Dict(branames .=> Base.Fix1(braname, nn).(ketnames))
-
-        renamed_message = replacedimnames(name -> get(bramap, name, name), messages[edge])
-
-        return operator(renamed_message, branames, ketnames)
-    end
-end
-
-@testset "apply_operator (T=$T, $(nameof(typeof(site_range))))" for site_range in (
-            spinone, spinone_u1,
+@testset "apply_operator (T=$T, $label)" for (label, site_range) in (
+            "spinone" => spinone, "spinone_u1" => spinone_u1, "fermion" => fermion,
         ),
         T in (Float32, Float64, ComplexF64)
 
@@ -63,17 +53,7 @@ end
         rng = StableRNG(123)
         g = named_cycle_graph(N)
         site_axes = Dict(v => Index(site_range) for v in vertices(g))
-        network = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
-
-        nn = NormNetwork(network)
-
-        env = beliefpropagation(
-            nn,
-            message_environment(msg -> state(fill!(msg, true)), nn);
-            stopping_criterion = (; maxiter = 100, tol = 1.0e-13)
-        )
-
-        env = operator_message_cache(nn, env)
+        network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
 
         for gate in (
                 randn_operator(rng, T, (site_axes[2],)),
@@ -88,17 +68,7 @@ end
         rng = StableRNG(123)
         g = named_path_graph(N)
         site_axes = Dict(v => Index(site_range) for v in vertices(g))
-        network = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
-
-        nn = NormNetwork(network)
-
-        env = beliefpropagation(
-            nn,
-            message_environment(msg -> state(fill!(msg, true)), nn);
-            stopping_criterion = (; maxiter = 100, tol = 1.0e-13)
-        )
-
-        env = operator_message_cache(nn, env)
+        network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
 
         gate = randn_operator(rng, T, (site_axes[2], site_axes[3]))
         gated_full = ITB.apply(gate, prod(network))
@@ -113,17 +83,7 @@ end
         rng = StableRNG(123)
         g = named_cycle_graph(N)
         site_axes = Dict(v => Index(site_range) for v in vertices(g))
-        network = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
-
-        nn = NormNetwork(network)
-
-        env = beliefpropagation(
-            nn,
-            message_environment(msg -> state(fill!(msg, true)), nn);
-            stopping_criterion = (; maxiter = 100, tol = 1.0e-13)
-        )
-
-        env = operator_message_cache(nn, env)
+        network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
 
         g1 = randn_operator(rng, T, (site_axes[2], site_axes[3]))
         g2 = randn_operator(rng, T, (site_axes[3], site_axes[4]))