Implement inner(x,A,y) for TTN via BP.

src/abstractitensornetwork.jl

@@ -721,6 +721,48 @@ function contract_inner(
   return contract(tn; sequence)[]
 end
 
+function contract_with_BP(
+  itn::AbstractITensorNetwork; outputlevel=1, partitioning=group(v -> v, vertices(itn))
+)
+  return contract_with_BP(
+    ComplexF64, itn::AbstractITensorNetwork; outputlevel, partitioning
+  )
+end
+
+function contract_with_BP(
+  T::Type,
+  itn::AbstractITensorNetwork;
+  outputlevel=1,
+  partitioning=group(v -> v, vertices(itn)),
+)
+  @assert isempty(externalinds(itn))
+  bp_cache = BeliefPropagationCache(copy(itn), partitioning)
+
+  bp_cache = update(bp_cache)
+
+  pg = partitioned_itensornetwork(bp_cache)
+
+  if !is_tree(partitioned_graph(pg)) && outputlevel > 0
+    println("Partitioned graph is not a tree, result will be approximate!!")
+  end
+
+  log_numerator, log_denominator = 0, 0
+  for pv in partitionvertices(pg)
+    incoming_mts = incoming_messages(bp_cache, [pv])
+    local_state = ITensor[itn[v] for v in vertices(pg, pv)]
+    log_numerator += log(complex(ITensors.contract(vcat(incoming_mts, local_state))[]))
+  end
+  for pe in partitionedges(pg)
+    log_denominator += log(
+      complex(
+        ITensors.contract(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))))[]
+      ),
+    )
+  end
+  res = exp(log_numerator - log_denominator)
+  return T(res)
+end
+
 # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl,
 # or `norm_sqr` to match `LinearAlgebra.norm_sqr`
 norm_sqr(ψ::AbstractITensorNetwork; sequence) = contract_inner(ψ, ψ; sequence)

src/treetensornetworks/abstracttreetensornetwork.jl

@@ -365,20 +365,14 @@ end
 # Inner products
 #
 
-# TODO: implement using multi-graph disjoint union
-function inner(
-  y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; root_vertex=default_root_vertex(x, A, y)
-)
-  traversal_order = reverse(post_order_dfs_vertices(x, root_vertex))
-  check_hascommoninds(siteinds, A, x)
-  check_hascommoninds(siteinds, A, y)
+# TODO: Remove dispatch on AbstractTTN (unless we want to trigger warning for trees if BP is nonexact)  
+function inner(y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; kwargs...)
   ydag = sim(dag(y); sites=[])
-  x = sim(x; sites=[])
-  O = ydag[root_vertex] * A[root_vertex] * x[root_vertex]
-  for v in traversal_order[2:end]
-    O = O * ydag[v] * A[v] * x[v]
-  end
-  return O[]
+  blf = BilinearFormNetwork(A, ydag, x)
+  blf_scalar_bp = contract_with_BP(
+    tensornetwork(blf); outputlevel=1, partitioning=group(v -> first(v), vertices(blf))
+  )
+  return blf_scalar_bp
 end
 
 # TODO: implement using multi-graph disjoint

test/test_forms.jl

@@ -1,4 +1,5 @@
 using ITensors
+using ITensors: contract
 using Graphs
 using NamedGraphs
 using ITensorNetworks
@@ -15,7 +16,7 @@ using ITensorNetworks:
 using Test
 using Random
 
-@testset "FormNetworkss" begin
+@testset "FormNetworks" begin
   g = named_grid((1, 4))
   s_ket = siteinds("S=1/2", g)
   s_bra = prime(s_ket; links=[])
@@ -49,3 +50,4 @@ using Random
   @test underlying_graph(ket_network(qf)) == underlying_graph(ψket)
   @test underlying_graph(operator_network(qf)) == underlying_graph(A)
 end
+nothing

test/test_treetensornetworks/test_treetensornetwork.jl

@@ -0,0 +1,53 @@
+using ITensors
+using ITensors: contract
+using ITensorNetworks
+using ITensorNetworks: contract_with_BP, group
+using Test
+
+@testset "TTN constructor defaulting to link_space=1" begin
+  tooth_lengths = fill(5, 6)
+  c = named_comb_tree(tooth_lengths)
+  s = siteinds("S=1/2", c)
+  d = Dict()
+  for (i, v) in enumerate(vertices(s))
+    d[v] = isodd(i) ? "Up" : "Dn"
+  end
+  states = v -> d[v]
+  #test a few signatures
+  state = TTN(s, states)
+  lds = edge_data(linkdims(state))
+  @test all([isone(lds[k]) for k in keys(lds)])
+  state = TTN(s)
+  lds = edge_data(linkdims(state))
+  @test all([isone(lds[k]) for k in keys(lds)])
+end
+
+@testset "Inner products for TTN via BP" begin
+  Random.seed!(1234)
+  Lx, Ly = 3, 3
+  χ = 2
+  g = named_comb_tree((Lx, Ly))
+  s = siteinds("S=1/2", g)
+  y = TTN(randomITensorNetwork(ComplexF64, s; link_space=χ))
+  x = TTN(randomITensorNetwork(ComplexF64, s; link_space=χ))
+
+  A = TTN(ITensorNetworks.heisenberg(s), s)
+  #First lets do it with the flattened version of the network
+  xy = inner_network(x, y; combine_linkinds=true)
+  xy_scalar = contract(xy)[]
+  xy_scalar_bp = contract_with_BP(xy)
+
+  @test_broken xy_scalar ≈ xy_scalar_bp
+
+  #Now lets keep it unflattened and do Block BP to keep the partitioned graph as a tree
+  xy = inner_network(x, y; combine_linkinds=false)
+  xy_scalar = contract(xy)[]
+  xy_scalar_bp = contract_with_BP(xy; partitioning=group(v -> first(v), vertices(xy)))
+
+  @test xy_scalar ≈ xy_scalar_bp
+  # test contraction of three layers for expectation values
+  # for TTN inner with this signature passes via contract_with_BP
+  @test inner(prime(x), A, y) ≈
+    inner(x, apply(A, y; nsweeps=10, maxdim=16, cutoff=1e-10, init=y)) rtol = 1e-6
+end
+nothing