Working on optimal orders

examples/belief_propagation/bpexample.jl

@@ -4,6 +4,7 @@ using Metis
 using ITensorNetworks
 using Random
 using SplitApplyCombine
+using NamedGraphs
 
 using ITensorNetworks:
   belief_propagation,
@@ -35,6 +36,7 @@ function main()
   )
 
   mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+
   numerator_network = approx_network_region(
     ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])])
   )

src/beliefpropagation/beliefpropagation.jl

@@ -81,8 +81,8 @@ function belief_propagation_iteration(
   mts::DataGraph;
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
   compute_norm=false,
-  update_sequence::String="parallel",
-  edges=Graphs.edges(mts),
+  update_sequence::String="sequential",
+  edges = edge_update_order(undirected_graph(underlying_graph(mts))),
 )
   new_mts = copy(mts)
   if update_sequence != "parallel" && update_sequence != "sequential"
@@ -97,7 +97,6 @@ function belief_propagation_iteration(
       incoming_mts[e_in] for
       e_in in setdiff(boundary_edges(incoming_mts, [src(e)]; dir=:in), [reverse(e)])
     ]
-
     new_mts[src(e) => dst(e)] = update_message_tensor(
       tn, incoming_mts[src(e)], environment_tensornetworks; contract_kwargs
     )
@@ -120,17 +119,16 @@ function belief_propagation(
   niters=20,
   target_precision::Union{Float64,Nothing}=nothing,
   update_sequence::String="parallel",
+  edges = edge_update_order(undirected_graph(underlying_graph(mts)))
 )
   compute_norm = target_precision == nothing ? false : true
   for i in 1:niters
     mts, c = belief_propagation_iteration(
-      tn, mts; contract_kwargs, compute_norm, update_sequence
+      tn, mts; contract_kwargs, compute_norm, update_sequence, edges
     )
     if compute_norm && c <= target_precision
       println(
-        "Belief Propagation finished. Reached a canonicalness of " *
-        string(c) *
-        " after $i iterations. ",
+        "BP converged to desired precision after $i iterations.",
       )
       break
     end

src/utils.jl

@@ -25,3 +25,15 @@ function line_to_tree(line::Vector)
   end
   return [line_to_tree(line[1:(end - 1)]), line[end]]
 end
+
+#Find an optimal ordering of the edges in an undirected graph
+function edge_update_order(g::AbstractNamedGraph)
+  es = []
+  for v in vertices(g)
+      new_es = reverse(reverse.(edges(bfs_tree(g, v))))
+      push!(es, setdiff(new_es, es)...)
+  end
+
+  @assert Set(es) == Set(vcat(edges(g), reverse.(edges(g))))
+  return es
+end

test/test_belief_propagation.jl

@@ -22,7 +22,7 @@ ITensors.disable_warn_order()
 
 @testset "belief_propagation" begin
 
-  #FIRST TEST SINGLE SITE ON AN MPS, SHOULD BE EXACT
+  #First test on an MPS, should be exact
   g_dims = (1, 6)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -44,7 +44,7 @@ ITensors.disable_warn_order()
     ψψ,
     mts;
     contract_kwargs=(; alg="exact"),
-    target_precision=1e-8,
+    niters = 1,
     update_sequence="sequential",
   )
 
@@ -56,8 +56,40 @@ ITensors.disable_warn_order()
 
   @test abs.(bp_sz - exact_sz) <= 1e-14
 
-  #NOW TEST TWO_SITE_EXPEC TAKING ON THE PARTITION FUNCTION OF THE RECTANGULAR ISING. SHOULD BE REASONABLE AND 
-  #INDEPENDENT OF INIT CONDITIONS, FOR SMALL BETA
+  #Now test on a tree, should also be exact
+  g = named_comb_tree((6,6))
+  s = siteinds("S=1/2", g)
+  χ = 2
+  Random.seed!(1564)
+  ψ = randomITensorNetwork(s; link_space=χ)
+
+  ψψ = ψ ⊗ prime(dag(ψ); sites=[])
+
+  v = (1, 3)
+
+  Oψ = copy(ψ)
+  Oψ[v] = apply(op("Sz", s[v]), ψ[v])
+  exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ)
+
+  Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
+  mts = message_tensors(Z)
+  mts = belief_propagation(
+    ψψ,
+    mts;
+    contract_kwargs=(; alg="exact"),
+    niters = 1,
+    update_sequence="sequential",
+  )
+
+  numerator_network = approx_network_region(
+    ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])])
+  )
+  denominator_network = approx_network_region(ψψ, mts, [(v, 1)])
+  bp_sz = contract(numerator_network)[] / contract(denominator_network)[]
+
+  @test abs.(bp_sz - exact_sz) <= 1e-14
+
+  #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close
   g_dims = (3, 4)
   g = named_grid(g_dims)
   s = IndsNetwork(g; link_space=2)
@@ -84,7 +116,7 @@ ITensors.disable_warn_order()
 
   @test abs.(bp_szsz - actual_szsz) <= 0.05
 
-  #TEST FORMING OF A TWO SITE RDM. JUST CHECK THAT IS HAS THE CORRECT SIZE, TRACE AND IS PSD
+  #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD
   g_dims = (4, 4)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -118,7 +150,7 @@ ITensors.disable_warn_order()
   @test all(>=(0), real(eigs)) && all(==(0), imag(eigs))
   @test size(rdm) == (2^length(vs), 2^length(vs))
 
-  #TEST MORE ADVANCED BP WITH ITENSORNETWORK MESSAGE TENSORS. IN THIS CASE IT SHOULD BE LIKE BOUNDARY MPS
+  #Test more advanced block BP with MPS message tensors on a grid 
   g_dims = (5, 4)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -144,6 +176,8 @@ ITensors.disable_warn_order()
     contract_kwargs=(;
       alg="density_matrix", output_structure=path_graph_structure, cutoff=1e-16, maxdim
     ),
+    niters = 1,
+    update_sequence="sequential",
   )
 
   numerator_network = approx_network_region(ψψ, mts, [v]; verts_tn=ITensorNetwork(ψOψ[v]))
@@ -155,4 +189,4 @@ ITensors.disable_warn_order()
     contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16)
 
   @test abs.(bp_sz - exact_sz) <= 1e-5
-end
+end