[WIP] Add fermionic support to belief-propagation simple update

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]
@@ -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"

src/apply/apply_operators.jl

@@ -1,13 +1,14 @@
 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, replacedimnames
+using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
+    outputnames, 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 +206,25 @@ 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.
+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...
@@ -233,7 +253,7 @@ function apply_gate_bp_nsite!(
     ψv = ITB.apply(op, state[v])
     if normalize
         gauges = [
-            conj(gram_eigh_full(env[e]))
+            first(message_gauge(env[e]))
                 for e in boundary_edges(state, vs; dir = :in)
         ]
         ψv /= norm(prod([[ψv]; gauges]))
@@ -249,12 +269,10 @@ 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))
+    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])
@@ -267,17 +285,16 @@ 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
 
     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,))
-    env[v2 => v1] = operator(
-        conj(replacedimnames(S, name_v1 => name_v2, name_v2 => name_v1)),
-        (name_v2,), (name_v1,)
-    )
+    # 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

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, inputnames, operator, outputnames, 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 = state.(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, outputnames(old_message), inputnames(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

src/beliefpropagation/messagecache.jl

@@ -209,13 +209,14 @@ 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 = linknames(braview, edge)
 
-            # Message axis is conj to the tensor it points to.
-            message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
+            # 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
         end