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 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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