Cluster expansions

.gitignore

@@ -1,9 +1,12 @@
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
+*.~
 /docs/Manifest*.toml
 /docs/build/
 
 Manifest.toml
 .vscode/
-.DS_Store
+.DS_Store
+.*
+*.ipynb

Project.toml

@@ -10,8 +10,10 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5"
+GraphIO = "aa1b3936-2fda-51b9-ab35-c553d3a640a2"
 GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
@@ -32,8 +34,10 @@ Adapt = "4.3.0"
 Combinatorics = "1.0.3"
 Dictionaries = "0.4"
 EinExprs = "0.6.4"
+GraphIO = "0.7.1"
 GraphRecipes = "0.5.13"
 Graphs = "1.8.0"
+HyperDualNumbers = "4.0.10"
 ITensorMPS = "0.3.17"
 ITensors = "0.9"
 KrylovKit = "0.10.2"

cluster/expect-corrected.jl

@@ -0,0 +1,128 @@
+using TensorNetworkQuantumSimulator
+using NamedGraphs
+using NamedGraphs: AbstractNamedGraph
+using Graphs
+
+const G = Graphs
+const NG = NamedGraphs
+const TN = TensorNetworkQuantumSimulator
+using HyperDualNumbers
+using Adapt: adapt
+using Dictionaries
+
+function prep_insertions(obs)
+    if isnothing(obs)
+        return (coeffs = identity, op_strings = v->"I")
+    end
+    op_strings, verts, _ = TN.collectobservable(obs)
+    @assert length(verts) <= 2
+
+    function hyper_coeff(v)
+        if v==verts[1]
+	    return Hyper(0,1,0,0)
+	elseif length(verts)==2 && v==verts[2]
+	    return Hyper(0,0,1,0)
+	else
+	    return 1
+	end
+    end
+
+    function insertion_operator(v)
+        if v==verts[1]
+	    return op_strings[1]
+	elseif length(verts)==2 && v==verts[2]
+	    return op_strings[2]
+	else
+	    return "I"
+	end
+    end
+    return (coeffs = hyper_coeff, op_strings = insertion_operator)
+end
+
+"""
+Cluster expansion. See clustercorrections.jl
+"""
+function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector, interaction_graph; obs = nothing)
+
+    kwargs = prep_insertions(obs)
+
+    logZbp = TN.free_energy(bpc; kwargs...)
+    isempty(egs) && return [0], [[logZbp]], [[1]]
+
+    circuit_lengths = sort(unique([c.weight for c=clusters]))
+
+    # Rescale the messages, but deal with the vertices separately
+    TN.rescale_messages!(bpc)
+    vns = Dictionary(TN.vertex_scalar(bpc, v; use_epsilon = true, kwargs...) for v=vertices(network(bpc).tensornetwork.graph))
+
+    # calculate weight of each generalized loop first
+    wts = TN.weights(bpc, egs; rescales = vns, kwargs...)
+
+    logZs = Array{Array}(undef, length(circuit_lengths) + 1)
+    logZs[1] = [logZbp]
+
+    coeffs = Array{Array}(undef, length(circuit_lengths) + 1)
+    coeffs[1] = [1]
+
+    # now calculate contribution to logZ from each cluster
+    for (cl_i, cl)=enumerate(circuit_lengths)
+        clusters_cl = filter(c->c.weight==cl, clusters)
+	logZs[cl_i + 1] = [prod([prod(fill(wts[l],c.multiplicities[l])) for l=c.loop_ids]) for c=clusters_cl]
+	coeffs[cl_i + 1] = [TN.ursell_function(c, interaction_graph) for c=clusters_cl]
+    end
+
+    return vcat([0],circuit_lengths), logZs, coeffs
+end	    	    
+
+"""
+Cluster cumulant expansion. See cumulant-clustercorrections.jl
+"""
+function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict; obs = nothing, rescale::Bool = false)
+
+    kwargs = prep_insertions(obs)
+
+    use_g = findall(gg->counting_nums[gg] != 0, regions)
+    egs = [induced_subgraph(network(bpc).tensornetwork.graph, gg)[1] for gg=regions[use_g]]
+
+    isempty(egs) && return logZbp, [], []
+
+    # Rescale the messages, but deal with the vertices separately
+    if rescale
+        TN.rescale_messages!(bpc)
+    	vns = Dictionary(TN.vertex_scalar(bpc, v; use_epsilon = true, kwargs...) for v=vertices(network(bpc).tensornetwork.graph))
+    else
+        vns = Dictionary(1 for v=vertices(network(bpc).tensornetwork.graph))
+    end
+
+    # calculate weight of each cluster first
+    wts = TN.weights(bpc, egs; rescales = vns, project_out = false, kwargs...)
+
+    return log.(wts), [counting_nums[gg] for gg=regions[use_g]]
+end
+
+"""
+onepoint or twopoint connected correlation function, using cluster cumulant expansion
+"""
+function cc_correlation(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict, obs)
+    logZs, cnums = cc_weights(bpc, regions, counting_nums; obs = obs)
+    op_strings, verts, _ = TN.collectobservable(obs)
+    if length(verts)==1
+        return sum(logZs .* cnums).epsilon1
+    else
+        return sum(logZs .* cnums).epsilon12
+    end
+end
+
+"""
+onepoint or twopoint connected correlation function, using cluster expansion
+"""
+function cluster_correlation(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector, interaction_graph, obs)
+    cluster_wts, logZs, ursells = cluster_weights(bpc, clusters, egs, interaction_graph; obs = obs)
+    op_strings, verts, _ = TN.collectobservable(obs)
+    cumul_dat = cumsum([sum([logZs[i][j] * ursells[i][j] for j=1:length(logZs[i])]) for i=1:length(logZs)])
+    if length(verts)==1
+        return cluster_wts, [d.epsilon1 for d=cumul_dat]
+    else
+        return cluster_wts, [d.epsilon12 for d=cumul_dat]
+    end
+end

examples/clustercorrections.jl

@@ -0,0 +1,97 @@
+using TensorNetworkQuantumSimulator
+const TN = TensorNetworkQuantumSimulator
+
+using ITensors
+
+using NamedGraphs
+using Graphs
+const NG = NamedGraphs
+const G = Graphs
+using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph
+
+using LinearAlgebra: norm
+
+using EinExprs: Greedy
+
+using Random
+Random.seed!(1634)
+
+include("../cluster/expect-corrected.jl")
+
+function main(nx,ny)
+    χ = 3
+    ITensors.disable_warn_order()
+    gs = [
+        (named_grid((nx, 1)), "line", 0,-1),
+        (named_hexagonal_lattice_graph(nx, ny), "hexagonal", 6,11),
+        (named_grid((nx, ny)), "square", 4,11),
+    ]
+
+    states = []
+    for (g, g_str, smallest_loop_size, wmax) in gs
+        println("*****************************************")
+        println("Testing for $g_str lattice with $(NG.nv(g)) vertices")
+	wmax = min(wmax, NG.nv(g))
+        ψ = TN.random_tensornetworkstate(ComplexF32, g, "S=1/2"; bond_dimension = χ)
+
+        ψ = normalize(ψ; alg = "bp")
+	ψIψ = BeliefPropagationCache(ψ)
+	ψIψ = update(ψIψ)
+
+	# BP expectation value
+	v = first(center(g))
+	expect_bp = real(expect(ψIψ, ("Z", [v])))
+        expect_exact_v = real(expect(ψ, ("Z", [v]); alg = "exact"))
+	clusters, egs, ig = TN.enumerate_clusters(g, wmax; must_contain=[v], min_deg = 1, min_v = smallest_loop_size)
+	cluster_wts, expects = cluster_correlation(ψIψ,clusters, egs, ig, ("Z", [v]))
+
+
+	regs = Dict()
+	cnums = Dict()
+
+	cc_wts = [1; smallest_loop_size:wmax;]
+	for w=cc_wts
+            regs[w],_,cnums[w]=TN.build_region_family_correlation(g,v,v,w)
+        end
+
+	expects_cc = Dict()
+	for w=cc_wts
+    	    expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], ("Z", [v])))
+        end
+
+        println("Bp expectation value for Z on site $(v) is $expect_bp")
+	println("Cluster expansion expectation values: $(cluster_wts), $(real.(expects))")
+	println("Cluster cumulant expansion: $(cc_wts), $([expects_cc[w] for w=cc_wts])")
+        println("Exact expectation value is $expect_exact_v")
+
+	println("***********************************")
+	u = neighbors(g, v)[1]
+	obs = (["Z","Z"], [u,v])
+	expect_exact_u = real(expect(ψ, ("Z", [u]); alg = "exact"))
+	expect_exact = real(expect(ψ,obs; alg = "exact")) - expect_exact_u * expect_exact_v
+	println("Calculating connected correlation function between $(v) and $(u)")
+
+	clusters, egs, ig = TN.enumerate_clusters(g, max(1,min(wmax,2*smallest_loop_size)); must_contain=[u,v], min_deg = 1, min_v = 2)
+
+	cluster_wts, expects = cluster_correlation(ψIψ,clusters, egs, ig, obs)
+
+	regs = Dict()
+	cnums = Dict()
+
+	cc_wts = [2;3:wmax;]
+	for w=cc_wts
+            regs[w],_,cnums[w]=TN.build_region_family_correlation(g,u,v,w)
+        end
+
+	expects_cc = Dict()
+	for w=cc_wts
+    	    expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], obs))
+        end
+	println("Cluster expansion expectation values: $(cluster_wts), $(real.(expects))")
+	println("Cluster cumulant expansion: $(cc_wts), $([expects_cc[w] for w=cc_wts])")
+
+        println("Exact expectation value is $expect_exact")
+	push!(states, ψ)
+    end
+    return states
+end

src/Apply/full_update.jl

@@ -107,6 +107,7 @@ function optimise_p_q(
         nfullupdatesweeps = 10,
         print_fidelity_loss = false,
         envisposdef = true,
+	verbose = false,
         apply_kwargs...,
     )
     p_cur, q_cur = factorize(
@@ -138,16 +139,17 @@ function optimise_p_q(
         b_vec = b(p, q, o, envs, q_cur)
         M_p_partial = partial(M_p, envs, q_cur, qs_ind)
 
-        p_cur, info = linsolve(
+        p_cur, info1 = linsolve(
             M_p_partial, b_vec, p_cur; isposdef = envisposdef, ishermitian = false
         )
 
         b_tilde_vec = b(p, q, o, envs, p_cur)
         M_p_tilde_partial = partial(M_p, envs, p_cur, ps_ind)
 
-        q_cur, info = linsolve(
+        q_cur, info2 = linsolve(
             M_p_tilde_partial, b_tilde_vec, q_cur; isposdef = envisposdef, ishermitian = false
         )
+	verbose && println("Linsolve info, iteration $(i): $(info1), $(info2)")
     end
 
     fend = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0