diff --git a/README.md b/README.md index 0cf6d0a..9ed7919 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,7 @@ # BPGauge - - - - +[![Build Status](https://github.com/CodingThrust/BPGauge.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/CodingThrust/BPGauge.jl/actions/workflows/CI.yml?query=branch%3Amain) + ## Introduction @@ -13,7 +11,7 @@ A simple example is to consider a cycle graph: ```julia using BPGauge, Graphs -g = chain(100) +g = path_graph(100) add_edge!(g, 1, 100) # generate a random state @@ -23,7 +21,7 @@ normalize_state!(tn) # generate a BP state, and do BP until convergence bp_state = BPState(tn) bp_path = BPPath(tn) -bp!(bp_state, bp_path, tn, err_bound = 1e-8) +bp!(bp_state, bp_path, tn, atol = 1e-8) # gauge transform the state gauge!(tn, bp_state) diff --git a/exm/Note.md b/exm/Note.md index 590d4dc..e62ad67 100644 --- a/exm/Note.md +++ b/exm/Note.md @@ -12,7 +12,7 @@ More complex geometry and property will reduce more the dimension of system by i $D=2^N \xrightarrow{\text{Constraint}} \alpha^N \xrightarrow{\text{Particle Number}}\alpha^{N-1}/k\xrightarrow{\text{translation}}\alpha^{N-1}/kN\xrightarrow{\text{spin flip/inversion}}\alpha^{N-2}/kN$. -where $\alpha \sim 1.618, 1.502$ for Fibonacci chain, hard core square lattice, $k \sim 5$. So for fully ED, $D \sim 10^6$. When doing dynamics, Krylov subspace $K_m​=\text{span}\{v,Hv,H^2v,⋯,H^{m−1}v\}$ is usually incorporated, which can scale up to $N=32$ even without anyon symmetry ($D \sim 10^{10}, m \sim 10^4$, as long as the locality of interaction ensuring sparse matrix). In such basis, the Its complexity decreases to $O(m^2D)$ +where $\alpha \sim 1.618, 1.502$ for Fibonacci path_graph, hard core square lattice, $k \sim 5$. So for fully ED, $D \sim 10^6$. When doing dynamics, Krylov subspace $K_m​=\text{span}\{v,Hv,H^2v,⋯,H^{m−1}v\}$ is usually incorporated, which can scale up to $N=32$ even without anyon symmetry ($D \sim 10^{10}, m \sim 10^4$, as long as the locality of interaction ensuring sparse matrix). In such basis, the Its complexity decreases to $O(m^2D)$ ## *Tensor Networks, TN* @@ -95,7 +95,7 @@ Supplemental Material](https://arxiv.org/pdf/1911.04882), which maybe helpful? Image [^Leshouches]: Les houch notes for Exact Diagonalizaton. https://indico.ictp.it/event/a14246/session/31/contribution/51/material/0/0.pdf -[^QMBSPRB]: Quantum scarred eigenstates in a Rydberg atom chain: Entanglement, breakdown of thermalization, and stability to perturbations. https://journals.aps.org/prb/pdf/10.1103/PhysRevB.98.155134 +[^QMBSPRB]: Quantum scarred eigenstates in a Rydberg atom path_graph: Entanglement, breakdown of thermalization, and stability to perturbations. https://journals.aps.org/prb/pdf/10.1103/PhysRevB.98.155134 [^LauchliKagome]: S = 1/2 kagome Heisenberg antiferromagnet revisited https://journals.aps.org/prb/pdf/10.1103/PhysRevB.100.155142 [^Dowling]: Monte Carlo techniques for real-time quantum dynamics https://arxiv.org/pdf/quant-ph/0507003 diff --git a/src/BPGauge.jl b/src/BPGauge.jl new file mode 100644 index 0000000..bd48c68 --- /dev/null +++ b/src/BPGauge.jl @@ -0,0 +1,30 @@ +module BPGauge + +using LinearAlgebra, Graphs, Random +using OMEinsum + +export TensorNetworkAnsatz +export BPState, BPPath, BPStep + +export zero_state, random_state, inner_product, normalize_state! +export bp!, bp_update! +export absorb! +export apply_gauge!, gauge! + +export square_lattice + +include("utils.jl") + +# define and construct the network +include("ansatz.jl") + +# bp on the tensor network ansatz +include("bp.jl") + +# operations about gauging +include("gauge.jl") + +# simple update +include("su.jl") + +end diff --git a/src/ansatz.jl b/src/ansatz.jl new file mode 100644 index 0000000..3a87f8f --- /dev/null +++ b/src/ansatz.jl @@ -0,0 +1,83 @@ +# the tensor network ansatz, |ket> +struct TensorNetworkAnsatz{TA <: AbstractArray, TB <: AbstractArray} + g::SimpleGraph{Int} # g is used to store the connections between the tensors + site_tensors::Vector{TA} # the tensors are stored as vector of arrays, corresponding to vertices of the graph. The virtual bounds are in order of its neighbors, the last dimension is the open dimension + gauge_tensors::Vector{TB} # corresponding to edges of the graph, listed in the order of edges(g), and enforce e.src < e.dst, real valued diagonal matrices + gauge_tensors_map::Dict{Tuple{Int, Int}, Int} # maps between the edge and the index of the gauge tensor + + function TensorNetworkAnsatz(g::SimpleGraph{Int}, d_virtual::Int, d_open::Int; type::Type = Float64) + TA = Array{Complex{type}} + TB = Array{type} + + site_tensors = Vector{TA}() + gauge_tensors = Vector{TB}() + + for v in vertices(g) + t = zeros(Complex{type}, [d_virtual for _ in 1:length(neighbors(g, v))]..., d_open) + push!(site_tensors, t) + end + + for e in edges(g) + t = (one(type) * I)(d_virtual) + push!(gauge_tensors, t) + end + + map = Dict{Tuple{Int, Int}, Int}() + for (i, e) in enumerate(edges(g)) + src, dst = minmax(e.src, e.dst) + map[(src, dst)] = i + end + + new{TA, TB}(g, site_tensors, gauge_tensors, map) + end + + function TensorNetworkAnsatz(g::SimpleGraph{Int}, site_tensors::Vector{TA}, gauge_tensors::Vector{TB}, map::Dict{Tuple{Int, Int}, Int}) where {TA, TB} + @assert length(site_tensors) == nv(g) + @assert length(gauge_tensors) == ne(g) + new{TA, TB}(g, site_tensors, gauge_tensors, map) + end +end + +Base.show(io::IO, ansatz::TensorNetworkAnsatz{TA, TB}) where {TA, TB} = print(io, "TensorNetworkAnsatz{$(TA), $(TB)} with $(nv(ansatz.g)) sites and $(ne(ansatz.g)) edges") + +Base.adjoint(ansatz::TensorNetworkAnsatz{TA, TB}) where {TA, TB} = TensorNetworkAnsatz(ansatz.g, TA[conj(t) for t in ansatz.site_tensors], ansatz.gauge_tensors, ansatz.gauge_tensors_map) + +function inner_product(bra::TensorNetworkAnsatz{TA, TB}, ket::TensorNetworkAnsatz{TA}; optimizer = GreedyMethod()) where {TA, TB} + @assert bra.g == ket.g + + raw_code = inner_product_eins(bra.g) + + xs = [bra.site_tensors..., bra.gauge_tensors..., ket.site_tensors..., ket.gauge_tensors...] + size_dict = OMEinsum.get_size_dict!(raw_code.ixs, xs, Dict{Int, Int}()) + + opt_code = optimize_code(raw_code, size_dict, optimizer) + # @info "contraction complexity: $(contraction_complexity(opt_code, size_dict))" + + res = opt_code(xs...) + return res[] +end + +# initialize the state for rydberg system, all sites are in the |0> state +function zero_state(g::SimpleGraph{Int}; d_virtual::Int = 1) + phi = TensorNetworkAnsatz(g, d_virtual, 2) + # set all to the |0> state + for tensor in phi.site_tensors + tensor[1] = ComplexF64(1.0) + end + return phi +end + +# generate non-normalized random state +function random_state(g::SimpleGraph{Int}; d_virtual::Int = 1) + phi = TensorNetworkAnsatz(g, d_virtual, 2) + for tensor in phi.site_tensors + tensor .= rand(ComplexF64, size(tensor)) + end + return phi +end + +function normalize_state!(phi::TensorNetworkAnsatz{TA, TB}) where {TA, TB} + n = inner_product(phi, adjoint(phi)) + phi.site_tensors .= phi.site_tensors ./ (sqrt(real(n))^(1 / nv(phi.g))) + return phi +end \ No newline at end of file diff --git a/src/bp.jl b/src/bp.jl new file mode 100644 index 0000000..06ac5be --- /dev/null +++ b/src/bp.jl @@ -0,0 +1,99 @@ +# bp on the TensorNetworkAnsatz + +# each message tensor have two legs, as follows +# --Ta -- |--1--Tb -- +# Mab +# --Ta*-- |--2--Tb*-- +struct BPState{TA} + messages::Dict{Tuple{Int, Int}, TA} + function BPState(ansatz::TensorNetworkAnsatz{TA, TB}) where {TA, TB} + messages = Dict{Tuple{Int, Int}, TA}() + TE = eltype(TA) + + # initialize the messages from neighbors to sites + for dst in vertices(ansatz.g) + for (i, src) in enumerate(neighbors(ansatz.g, dst)) + d_virtual = size(ansatz.site_tensors[dst])[i] + messages[(src, dst)] = ones(TE, d_virtual, d_virtual) + normalize_message!(messages[(src, dst)]) + end + end + + new{TA}(messages) + end +end + +# a single bp update step: \sum_{T_neighbors} T * T* * M_inputs... -> M_output +struct BPStep + eincode::AbstractEinsum # Eincode([T..., T*..., M_inputs...], [M_output]) + Tid::Int # the id of the tensor to interact with + inputs::Vector{Tuple{Int, Int}} # input messages + output::Tuple{Int, Int} # output message +end + +struct BPPath + bp_steps::Vector{BPStep} # steps of bp update, each step is a single bp update + function BPPath(ansatz::TensorNetworkAnsatz{TA, TB}; random_order::Bool = true, seed::Int = 1234) where {TA, TB} + g = ansatz.g + Random.seed!(seed) + bp_steps = BPStep[] + # the different step are taken in random order + outputs = Tuple{Int, Int}[] + for e in edges(g) + push!(outputs, (src(e), dst(e))) + push!(outputs, (dst(e), src(e))) + end + + # the steps are taken in random order + random_order && shuffle!(outputs) + + for (s, d) in outputs + inputs = Tuple{Int, Int}[] + Tid = s + for v in neighbors(g, s) + if v != d + push!(inputs, (v, s)) + end + end + + eincode = bp_eins(neighbors(g, s), inputs, (s, d)) + push!(bp_steps, BPStep(eincode, Tid, inputs, (s, d))) + end + + new(bp_steps) + end +end + +function bp!(state::BPState, path::BPPath, ansatz::TensorNetworkAnsatz; max_iter::Int = 1000, atol::Float64 = 1e-6, damping::Float64 = 0.2, verbose::Bool = false) + Ts = ansatz.site_tensors + Tcs = conj.(ansatz.site_tensors) + for iter in 1:max_iter + error = bp_update!(state, path, Ts, Tcs, damping) + verbose && @info "BP iteration $iter, error: $error" + error < atol && break + end + nothing +end + +function bp_update!(state::BPState, path::BPPath, Ts::Vector{TA}, Tcs::Vector{TB}, damping::Float64) where {TA, TB} + + error = 0.0 + for step in path.bp_steps + eincode = step.eincode + Tid = step.Tid + inputs = step.inputs + output = step.output + + T = Ts[Tid] + Tc = Tcs[Tid] + + M_inputs = [state.messages[input] for input in inputs] + + # update the message tensor + new_message = normalize_message!(eincode(T, Tc, M_inputs...)) + error = max(error, maximum(abs.(state.messages[output] - new_message))) + state.messages[output] .*= damping + state.messages[output] .+= (1 - damping) * new_message + end + return error +end \ No newline at end of file diff --git a/src/gauge.jl b/src/gauge.jl new file mode 100644 index 0000000..120965b --- /dev/null +++ b/src/gauge.jl @@ -0,0 +1,87 @@ +# absorb the gauge tensors (diagonal matrices) into the site tensors, for bp to work +# after the gauge is absorbed, set them to identity +function absorb!(ansatz::TensorNetworkAnsatz{TA, TB}) where {TA, TB} + for e in edges(ansatz.g) + s, d = minmax(e.src, e.dst) + Gamma = ansatz.gauge_tensors[ansatz.gauge_tensors_map[(s, d)]] + sqrt_Gamma = sqrt.(Gamma) + + eincode_s, eincode_d = absorb_eins(s, d, neighbors(ansatz.g, s), neighbors(ansatz.g, d)) + + s_xs = (ansatz.site_tensors[s], sqrt_Gamma) + sds = OMEinsum.get_size_dict(eincode_s.ixs, s_xs) + einsum!(eincode_s.ixs, eincode_s.iy, s_xs, ansatz.site_tensors[s], 1, 0, sds) + + d_xs = (ansatz.site_tensors[d], sqrt_Gamma) + sdd = OMEinsum.get_size_dict(eincode_d.ixs, d_xs) + einsum!(eincode_d.ixs, eincode_d.iy, d_xs, ansatz.site_tensors[d], 1, 0, sdd) + + # reset the gauge tensor to identity + Gamma .= I(size(Gamma, 1)) + end + nothing +end + +# assume that the existing gauge tensors are already absorbed into the site tensors +# the bp converge, using the bp messages to generate the new gauge tensors, which are real valued diagonal matrices +# according to Eqs.12~17 of https://scipost.org/SciPostPhys.15.6.222 +# site s and site d +# Ts -i- sqrt_Msd_inv (j, i) -j- sqrt_Msd (k, j) -k- sqrt_Mds (k, l) -l- sqrt_Mds_inv (l, m) -m- Td +# Ts -i- sqrt_Msd_inv (j, i) -j- U(j, n) -n- S(n, o) -o- Vt(o, l) -l- sqrt_Mds_inv (l, m) -m- Td +# Ts -i- As (i, n) -n- Gamma -o- Ad (m, o) -m- Td +# Ts -i- Gamma -m- Td +function gauge!(ansatz::TensorNetworkAnsatz{TA, TB}, state::BPState{TA}) where {TA, TB} + for e in edges(ansatz.g) + s, d = minmax(e.src, e.dst) + M_sd = state.messages[(s, d)] + M_ds = state.messages[(d, s)] + + apply_gauge!(ansatz, s, d, M_sd, M_ds) + end + nothing +end + +# given the message matrix M_sd and M_ds, generate the new gauge tensor Gamma, update T_s and T_d +function apply_gauge!(ansatz::TensorNetworkAnsatz{TA, TB}, s::Int, d::Int, M_sd::Matrix{T}, M_ds::Matrix{T}) where {TA, TB, T} + # square root the message matrix + sqrt_Msd, sqrt_Msd_inv = square_root(M_sd) + sqrt_Mds, sqrt_Mds_inv = square_root(M_ds) + + # the eigen decomposition is only correct when the message matrix is positive definite + # @show ein"ji, kj -> ik"(sqrt_Msd_inv, sqrt_Msd) ≈ I(size(M_sd, 1)) + # @show ein"kl, lm -> km"(sqrt_Mds, sqrt_Mds_inv) ≈ I(size(M_sd, 1)) + # @show ein"ji, kj, kl, lm -> im"(sqrt_Msd_inv, sqrt_Msd, sqrt_Mds, sqrt_Mds_inv) ≈ I(size(M_sd, 1)) + + # generate the new gauge tensor + M_mid = ein"kj, kl -> jl"(sqrt_Msd, sqrt_Mds) + + # @show ein"ji, jl, lm -> im"(sqrt_Msd_inv, M_mid, sqrt_Mds_inv) ≈ I(size(M_mid, 1)) + + res = svd(M_mid) + U = res.U # jn + S = diagm(res.S) # no + Vt = res.Vt # ol + + # @show ein"jn, no, ol -> jl"(U, S, Vt) ≈ M_mid + # @show maximum(abs.(ein"ji, jn, no, ol, lm -> im"(sqrt_Msd_inv, U, S, Vt, sqrt_Mds_inv) - I(size(M_mid, 1)))) + + As = ein"ji, jn -> in"(sqrt_Msd_inv, U) + Ad = ein"lm, ol -> mo"(sqrt_Mds_inv, Vt) + + # @show maximum(abs.(ein"in, no, mo -> im"(As, S, Ad) - I(size(As, 1)))) + + #update the gauge tensor + ansatz.gauge_tensors[ansatz.gauge_tensors_map[(s, d)]] .= S + + # absorb As and Ad into the site tensors + eincode_s, eincode_d = absorb_eins(s, d, neighbors(ansatz.g, s), neighbors(ansatz.g, d)) + s_xs = (ansatz.site_tensors[s], As) + sds = OMEinsum.get_size_dict(eincode_s.ixs, s_xs) + einsum!(eincode_s.ixs, eincode_s.iy, s_xs, ansatz.site_tensors[s], 1, 0, sds) + + d_xs = (ansatz.site_tensors[d], Ad) + sdd = OMEinsum.get_size_dict(eincode_d.ixs, d_xs) + einsum!(eincode_d.ixs, eincode_d.iy, d_xs, ansatz.site_tensors[d], 1, 0, sdd) + + nothing +end \ No newline at end of file diff --git a/src/su.jl b/src/su.jl new file mode 100644 index 0000000..7e66341 --- /dev/null +++ b/src/su.jl @@ -0,0 +1,4 @@ +# TODO +# 1. single site update +# 2. two site update +# 3. truncation \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..8ed3869 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,117 @@ +# generating graphs +function square_lattice(mx::Int, my::Int, p::Float64; seed::Int = 1234) + Random.seed!(seed) + n = Int(ceil(mx * my * p)) + g = grid((mx, my)) + selected_vertices = shuffle!(collect(vertices(g)))[1:n] + sub_g, _ = induced_subgraph(g, selected_vertices) + return sub_g +end + +# codes about constructing the einsum expression +function inner_product_eins(g::SimpleGraph{Int}) + # 1:nv(g) are the indices of open indices + count = nv(g) + 1 + ixs_bra, count = all_eins(g, count) + ixs_ket, count = all_eins(g, count) + + ixs = vcat(ixs_bra, ixs_ket) + + return EinCode(ixs, Int[]) +end + +function all_eins(g::SimpleGraph{Int}, count::Int) + indices_dict = Dict{Tuple{Int, Int}, Int}() + ixs = Vector{Vector{Int}}() + + # site tensors + for v in vertices(g) + ix = Int[] + for n in neighbors(g, v) + push!(ix, count) + indices_dict[(v, n)] = count + count += 1 + end + push!(ix, v) + push!(ixs, ix) + end + + # gauge tensors + for e in edges(g) + src, dst = minmax(e.src, e.dst) + push!(ixs, [indices_dict[(src, dst)], indices_dict[(dst, src)]]) + end + + return ixs, count +end + +function bp_eins(nebis::Vector{Int}, inputs::Vector{Tuple{Int, Int}}, output::Tuple{Int, Int}) + ixs = Vector{Vector{Int}}() + iy = Int[] + + # 1 is indices of the open indices + # 2:d + 1 are indices of the virtual indices of T + # d + 2:2d + 1 are indices of the virtual indices of T* + d = length(nebis) + T_ids = [2:d+1..., 1] + Tc_ids = [d+2:2d+1..., 1] + push!(ixs, T_ids) + push!(ixs, Tc_ids) + + output_leg = findfirst(x -> x == output[2], nebis) + iy = [output_leg + 1, output_leg + 1 + d] + + for (i, j) in inputs + input_leg = findfirst(x -> x == i, nebis) + push!(ixs, [input_leg + 1, input_leg + 1 + d]) + end + + raw_code = EinCode(ixs, iy) + return optimize_code(raw_code, uniformsize(raw_code, 2), GreedyMethod()) +end + + +# the eincode for absorbing the square root of the gauge tensor into the site +function absorb_eins(s::Int, d::Int, nebis_s::Vector{Int}, nebis_d::Vector{Int}) + id_open = 0 + + iy_s = [nebis_s..., id_open] + iy_d = [nebis_d..., id_open] + + ix_s = [nebis_s..., id_open] + ix_s[findfirst(==(d), ix_s)] = s + + ix_d = [nebis_d..., id_open] + ix_d[findfirst(==(s), ix_d)] = d + + return EinCode([ix_s, [s, d]], iy_s), EinCode([ix_d, [d, s]], iy_d) +end + +# a very simple strategy to uniform the message tensor +function normalize_message!(t::AbstractArray) + return normalize!(t, 1) +end + +# square_root of the message matrix M_12 via eigen decomposition +# note that the definition here is a little bit different from the one in the paper https://scipost.org/SciPostPhys.15.6.222 +# with eigen decomposition, M_12 = U * diagm(D) * adjoint(U) = SM_13 * SM_32 +# SM_13 = M_12^0.5 = U * sqrt(diagm(D)), SM_32 = M_12^(-0.5) = 1/sqrt(diagm(D)) * adjoint(U) +# it is assumed that the message matrix is positive definte, thus all elements of D should be positive +function square_root(M::Matrix{Complex{T}}) where T + @assert M ≈ adjoint(M) # the message matrix is Hermitian + M = LinearAlgebra.hermitian(M) + res = eigen(M) + D = res.values + U = res.vectors + + # it seems that eigen decomposition some time leads to small negative values, fix them to 10 times the machine precision + # see https://github.com/ITensor/ITensorNetworks.jl/blob/2d7db49c2f003c6588a1632fd6f0fd802c0ac9ca/src/gauging.jl#L70 + # this operation is extremely ill-conditioned, leading to huge drop of accuracy + D .= max.(D, 10 * eps(T)) + + # need to consider the pseudoinverse if D contains 0 + sqrt_D = sqrt.(D) + sqrt_M = U * diagm(sqrt_D) + sqrt_M_inv = diagm(inv.(sqrt_D)) * adjoint(U) + return sqrt_M, sqrt_M_inv +end \ No newline at end of file diff --git a/test/ansatz.jl b/test/ansatz.jl new file mode 100644 index 0000000..5e76ba9 --- /dev/null +++ b/test/ansatz.jl @@ -0,0 +1,14 @@ +using BPGauge +using Graphs, LinearAlgebra, OMEinsum +using Test + +@testset "states" begin + for g in [random_regular_graph(30, 3), BPGauge.square_lattice(10, 10, 0.8)] + phi = BPGauge.zero_state(g, d_virtual = 2) + @test inner_product(phi, adjoint(phi)) ≈ 1.0 + + phi = BPGauge.random_state(g, d_virtual = 2) + normalize_state!(phi) + @test inner_product(phi, adjoint(phi)) ≈ 1.0 + end +end \ No newline at end of file diff --git a/test/bp.jl b/test/bp.jl new file mode 100644 index 0000000..cc6f84f --- /dev/null +++ b/test/bp.jl @@ -0,0 +1,78 @@ +using BPGauge +using Graphs, LinearAlgebra, OMEinsum +using Test + +using Random +Random.seed!(1234) + +@testset "bp" begin + g = random_regular_graph(30, 3) + + tn = zero_state(g, d_virtual = 2) + bp_state = BPState(tn) + @test length(bp_state.messages) == ne(g) * 2 + + bp_path = BPPath(tn) + @test length(bp_path.bp_steps) == ne(g) * 2 + for step in bp_path.bp_steps + @test step.output[1] == step.Tid + @test step.output[2] ∈ neighbors(g, step.Tid) + for input in step.inputs + @test input[2] == step.Tid + @test input[1] ∈ neighbors(g, step.Tid) + end + end + + @test isnothing(bp!(bp_state, bp_path, tn)) +end + +@testset "1d path_graph" begin + g = path_graph(10) + for d in [2, 3, 10] + tn = random_state(g, d_virtual = 3) + normalize_state!(tn) + bp_state = BPState(tn) + bp_path = BPPath(tn) + + bp!(bp_state, bp_path, tn) + + T1 = tn.site_tensors[1] + T2 = tn.site_tensors[2] + T3 = tn.site_tensors[3] + + M12 = ein"li, ni -> ln"(T1, conj(T1)) + M23 = ein"li, ni, lmj, noj -> mo"(T1, conj(T1), T2, conj(T2)) + + r12 = M12[1] / bp_state.messages[(1, 2)][1] + @test isapprox(M12, r12 * bp_state.messages[(1, 2)], rtol = 1e-6) + + r23 = M23[1] / bp_state.messages[(2, 3)][1] + @test isapprox(M23, r23 * bp_state.messages[(2, 3)], rtol = 1e-6) + end +end + +@testset "tree" begin + g = SimpleGraph(6) + add_edge!(g, 1, 2) + add_edge!(g, 2, 3) + add_edge!(g, 2, 4) + add_edge!(g, 3, 5) + add_edge!(g, 3, 6) + + for d in [2, 3, 10] + tn = random_state(g, d_virtual = d) + normalize_state!(tn) + bp_state = BPState(tn) + bp_path = BPPath(tn) + + bp!(bp_state, bp_path, tn) + # check the message 2 -> 3 + T1 = tn.site_tensors[1] + T2 = tn.site_tensors[2] + T4 = tn.site_tensors[4] + + M23 = ein"li, oi, lnmj, oqpj, mk, pk -> nq"(T1, conj(T1), T2, conj(T2), T4, conj(T4)) + r23 = M23[1] / bp_state.messages[(2, 3)][1] + @test isapprox(M23, r23 * bp_state.messages[(2, 3)], rtol = 1e-6) + end +end \ No newline at end of file diff --git a/test/gauge.jl b/test/gauge.jl new file mode 100644 index 0000000..77ec76a --- /dev/null +++ b/test/gauge.jl @@ -0,0 +1,114 @@ +using BPGauge +using Graphs, LinearAlgebra, OMEinsum +using Test + +using Random +Random.seed!(1234) + +using BPGauge: square_root + +@testset "absorb" begin + for g in [random_regular_graph(30, 3), BPGauge.square_lattice(10, 10, 0.8)] + tn = random_state(g, d_virtual = 2) + for gt in tn.gauge_tensors + gt .= diagm(rand(size(gt, 1))) + end + normalize_state!(tn) + absorb!(tn) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 + + tn = zero_state(g, d_virtual = 2) + for gt in tn.gauge_tensors + gt .*= 2.0 + end + absorb!(tn) + for i in 1:nv(g) + @test tn.site_tensors[i][1] ≈ 2^(0.5 * degree(g, i)) + end + end +end + +@testset "svd" begin + A = rand(ComplexF64, 10, 10) + res = svd(A) + @test A ≈ ein"ij, jk, kl -> il"(res.U, diagm(res.S), res.Vt) +end + +@testset "gauge" begin + # gauge by message tensor + g = SimpleGraph(5) + add_edge!(g, 1, 3) + add_edge!(g, 2, 3) + add_edge!(g, 3, 4) + add_edge!(g, 4, 5) + tn = random_state(g, d_virtual = 10) + normalize_state!(tn) + origin_34 = ein"ijkm, kln -> ijmnl"(tn.site_tensors[3], tn.site_tensors[4]) + + t = rand(ComplexF64, 10, 10) + M_34 = t * adjoint(t) + t = rand(ComplexF64, 10, 10) + M_43 = t * adjoint(t) + + apply_gauge!(tn, 3, 4, M_34, M_43) + + gauged_34 = ein"ijom, pln, op -> ijmnl"(tn.site_tensors[3], tn.site_tensors[4], tn.gauge_tensors[tn.gauge_tensors_map[(3, 4)]]) + @test isapprox(gauged_34, origin_34, atol = 1e-8) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 +end + +@testset "gauge by bp" begin + g = random_regular_graph(10, 3) + tn = random_state(g, d_virtual = 3) + normalize_state!(tn) + bp_state = BPState(tn) + bp_path = BPPath(tn) + bp!(bp_state, bp_path, tn, verbose = true) + + gauge!(tn, bp_state) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 +end + +@testset "Vidal gauge mps" begin + g = path_graph(100) + tn = random_state(g, d_virtual = 8) + normalize_state!(tn) + + bp_state = BPState(tn) + bp_path = BPPath(tn) + bp!(bp_state, bp_path, tn, atol = 1e-14) + + gauge!(tn, bp_state) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 + + for i in 2:nv(g) - 2 + G = tn.gauge_tensors[tn.gauge_tensors_map[(i, i + 1)]] + T = tn.site_tensors[i + 1] + L = ein"ij, jk, jln, kmn -> lm"(G, G, T, conj(T)) + @test maximum(abs.(L ./ L[1, 1] - I(size(L, 1)))) < 1e-6 + end +end + +@testset "cycle graph gauge" begin + g = path_graph(100) + add_edge!(g, 1, 100) + + tn = random_state(g, d_virtual = 4) + normalize_state!(tn) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 + + bp_state = BPState(tn) + bp_path = BPPath(tn) + bp!(bp_state, bp_path, tn, atol = 1e-8) + + gauge!(tn, bp_state) + @test inner_product(tn, adjoint(tn)) ≈ 1.0 + + for i in 1:nv(g) - 2 + j = i == nv(g) ? 1 : i + 1 + G = tn.gauge_tensors[tn.gauge_tensors_map[(min(i, j), max(i, j))]] + T = tn.site_tensors[j] + L = ein"ij, jk, jln, kmn -> lm"(G, G, T, conj(T)) + @test maximum(abs.(L ./ L[1, 1] - I(size(L, 1)))) < 2 * 1e-6 + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..8483d53 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,18 @@ +using BPGauge +using Test + +@testset "utils" begin + include("utils.jl") +end + +@testset "ansatz" begin + include("ansatz.jl") +end + +@testset "bp" begin + include("bp.jl") +end + +@testset "gauge" begin + include("gauge.jl") +end diff --git a/test/utils.jl b/test/utils.jl new file mode 100644 index 0000000..e8b4f7a --- /dev/null +++ b/test/utils.jl @@ -0,0 +1,75 @@ +using BPGauge +using Graphs, LinearAlgebra, OMEinsum +using Test + +using Random +Random.seed!(1234) + +@testset "graphs" begin + g = BPGauge.square_lattice(10, 10, 0.8) + @test nv(g) == 80 +end + +@testset "eins" begin + g = random_regular_graph(30, 3) + ixs, count = BPGauge.all_eins(g, 31) + @test count == 31 + 2 * ne(g) + @test length(ixs) == nv(g) + ne(g) + + vec_ixs = ixs[1:nv(g)] + edge_ixs = ixs[nv(g) + 1:end] + + for (i, e) in enumerate(edges(g)) + src, dst = minmax(e.src, e.dst) + @test edge_ixs[i][1] == vec_ixs[src][findfirst(x -> x == dst, neighbors(g, src))] + @test edge_ixs[i][2] == vec_ixs[dst][findfirst(x -> x == src, neighbors(g, dst))] + end +end + +@testset "bp eins" begin + nebis = [2, 3, 4] + inputs = [(2, 1), (4, 1)] + output = (1, 3) + eincode_generated = BPGauge.bp_eins(nebis, inputs, output) + eincode_manul = ein"ijkl, mnol, im, ko -> jn" + T1 = rand(3, 4, 5, 2) + T2 = rand(3, 4, 5, 2) + M1 = rand(3, 3) + M3 = rand(5, 5) + M2_generated = eincode_generated(T1, T2, M1, M3) + M2_manul = eincode_manul(T1, T2, M1, M3) + @test M2_generated ≈ M2_manul +end + +@testset "absorb_eins" begin + g = SimpleGraph(5) + add_edge!(g, 1, 3) + add_edge!(g, 2, 3) + add_edge!(g, 3, 4) + add_edge!(g, 4, 5) + s,d = 3, 4 + eincode_s, eincode_d = BPGauge.absorb_eins(s, d, neighbors(g, s), neighbors(g, d)) + eincode_s_man = ein"ijkm, ko -> ijom" + eincode_d_man = ein"kln, ok -> oln" + + T3 = rand(3, 4, 5, 2) + T4 = rand(5, 6, 2) + G = diagm(rand(5)) + + @test eincode_s(T3, G) ≈ eincode_s_man(T3, G) + @test eincode_d(T4, G) ≈ eincode_d_man(T4, G) +end + +@testset "square_root" begin + # generate a random positive definite matrix + A = rand(ComplexF64, 10, 10) + A = A * adjoint(A) + sqrt_A, sqrt_A_inv = BPGauge.square_root(A) + + # sqrt + @test sqrt_A * adjoint(sqrt_A) ≈ A + + # inverse + @test sqrt_A * sqrt_A_inv ≈ I(10) + @test sqrt_A_inv * sqrt_A ≈ I(10) +end \ No newline at end of file