diff --git a/Project.toml b/Project.toml index 08a513dc..8b4e9c75 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" -version = "0.17.0" +version = "0.18.0" authors = ["Matthew Fishman , Joseph Tindall and contributors"] [workspace] diff --git a/docs/Project.toml b/docs/Project.toml index 67f7ec0e..1920f12f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,7 +16,7 @@ path = ".." Documenter = "1.10" Graphs = "1" ITensorFormatter = "0.2.27" -ITensorNetworks = "0.17" +ITensorNetworks = "0.18" ITensors = "0.9" Literate = "2.20.1" NamedGraphs = "0.8.2" diff --git a/examples/Project.toml b/examples/Project.toml index ab5ad318..10dc8e59 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,4 +5,4 @@ ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" path = ".." [compat] -ITensorNetworks = "0.17" +ITensorNetworks = "0.18" diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index a5c44e79..dd7120f1 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -23,7 +23,6 @@ include("formnetworks/abstractformnetwork.jl") include("formnetworks/bilinearformnetwork.jl") include("formnetworks/linearformnetwork.jl") include("formnetworks/quadraticformnetwork.jl") -include("gauging.jl") include("utils.jl") include("update_observer.jl") diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 6ab598aa..a615c835 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -7,10 +7,10 @@ using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edge ne, neighbors, rem_edge!, src, vertices using ITensors: ITensors, @Algorithm_str, ITensor, addtags, combiner, commoninds, commontags, contract, dag, hascommoninds, inds, noprime, onehot, prime, replaceprime, - replacetags, setprime, settags, sim, swaptags, unioninds, uniqueinds + replacetags, setprime, settags, sim, swaptags, tags, unioninds, uniqueinds using LinearAlgebra: LinearAlgebra, factorize using MacroTools: @capture -using NDTensors: NDTensors, Algorithm, dim +using NDTensors: NDTensors, Algorithm, dim, scalartype using NamedGraphs.GraphsExtensions: directed_graph, incident_edges, rename_vertices, vertextype, ⊔ using NamedGraphs: NamedGraphs, NamedGraph, not_implemented, steiner_tree diff --git a/src/apply.jl b/src/apply.jl index c14f44ba..df5d254c 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -1,8 +1,9 @@ using .BaseExtensions: maybe_real using Graphs: has_edge +using ITensors.NDTensors: scalartype using ITensors: ITensors, ITensor, Index, Ops, apply, commonind, commoninds, contract, dag, denseblocks, factorize, factorize_svd, hasqns, isdiag, noncommoninds, noprime, prime, - replaceind, replaceinds, unioninds, uniqueinds + replaceind, replaceinds, tags, unioninds, uniqueinds using KrylovKit: linsolve using LinearAlgebra: eigen, norm, qr, svd using NamedGraphs: NamedEdge, has_edge @@ -304,80 +305,6 @@ function _contract_gate(o::AbstractEdge, ψv1, Λ, ψv2) return Qᵥ₁, Rᵥ₁, Qᵥ₂, Rᵥ₂, theta end -# In the future we will try to unify this into apply() above but currently leave it mostly as a separate function -# Apply() function for an ITN in the Vidal Gauge. Hence the bond tensors are required. -# Gate does not necessarily need to be passed. Can supply an edge to do an identity update instead. Uses Simple Update procedure assuming gate is two-site -function ITensors.apply( - o::Union{NamedEdge, ITensor}, ψ::VidalITensorNetwork; normalize = false, apply_kwargs... - ) - updated_ψ = copy(site_tensors(ψ)) - updated_bond_tensors = copy(bond_tensors(ψ)) - v⃗ = _gate_vertices(o, ψ) - if length(v⃗) == 2 - e = NamedEdge(v⃗[1] => v⃗[2]) - ψv1, ψv2 = ψ[src(e)], ψ[dst(e)] - e_ind = commonind(ψv1, ψv2) - - for vn in neighbors(ψ, src(e)) - if (vn != dst(e)) - ψv1 = noprime(ψv1 * bond_tensor(ψ, vn => src(e))) - end - end - - for vn in neighbors(ψ, dst(e)) - if (vn != src(e)) - ψv2 = noprime(ψv2 * bond_tensor(ψ, vn => dst(e))) - end - end - - Qᵥ₁, Rᵥ₁, Qᵥ₂, Rᵥ₂, theta = _contract_gate(o, ψv1, bond_tensor(ψ, e), ψv2) - - U, S, V = ITensors.svd( - theta, - uniqueinds(Rᵥ₁, Rᵥ₂); - lefttags = ITensorNetworks.edge_tag(e), - righttags = ITensorNetworks.edge_tag(e), - apply_kwargs... - ) - - ind_to_replace = commonind(V, S) - ind_to_replace_with = commonind(U, S) - S = replaceind(S, ind_to_replace => ind_to_replace_with') - V = replaceind(V, ind_to_replace => ind_to_replace_with) - - ψv1, updated_bond_tensors[e], ψv2 = U * Qᵥ₁, S, V * Qᵥ₂ - - for vn in neighbors(ψ, src(e)) - if (vn != dst(e)) - ψv1 = - noprime(ψv1 * ITensorsExtensions.inv_diag(bond_tensor(ψ, vn => src(e)))) - end - end - - for vn in neighbors(ψ, dst(e)) - if (vn != src(e)) - ψv2 = - noprime(ψv2 * ITensorsExtensions.inv_diag(bond_tensor(ψ, vn => dst(e)))) - end - end - - if normalize - ψv1 /= norm(ψv1) - ψv2 /= norm(ψv2) - updated_bond_tensors[e] /= norm(updated_bond_tensors[e]) - end - - setindex_preserve_graph!(updated_ψ, ψv1, src(e)) - setindex_preserve_graph!(updated_ψ, ψv2, dst(e)) - - return VidalITensorNetwork(updated_ψ, updated_bond_tensors) - - else - updated_ψ = apply(o, updated_ψ; normalize) - return VidalITensorNetwork(updated_ψ, updated_bond_tensors) - end -end - ### Full Update Routines ### # Calculate the overlap of the gate acting on the previous p and q versus the new p and q in the presence of environments. This is the cost function that optimise_p_q will minimise diff --git a/src/gauging.jl b/src/gauging.jl deleted file mode 100644 index a0819dfd..00000000 --- a/src/gauging.jl +++ /dev/null @@ -1,188 +0,0 @@ -using ITensors.NDTensors: dense, scalartype -using ITensors: tags -using NamedGraphs.PartitionedGraphs: quotientedge - -function default_bond_tensors(ψ::ITensorNetwork) - return DataGraph( - underlying_graph(ψ); vertex_data_eltype = Nothing, edge_data_eltype = ITensor - ) -end - -struct VidalITensorNetwork{V, BTS} <: AbstractITensorNetwork{V} - itensornetwork::ITensorNetwork{V} - bond_tensors::BTS -end - -site_tensors(ψ::VidalITensorNetwork) = ψ.itensornetwork -bond_tensors(ψ::VidalITensorNetwork) = ψ.bond_tensors -bond_tensor(ψ::VidalITensorNetwork, e) = bond_tensors(ψ)[e] - -function data_graph_type(TN::Type{<:VidalITensorNetwork}) - return data_graph_type(fieldtype(TN, :itensornetwork)) -end -data_graph(ψ::VidalITensorNetwork) = data_graph(site_tensors(ψ)) -function Base.copy(ψ::VidalITensorNetwork) - return VidalITensorNetwork(copy(site_tensors(ψ)), copy(bond_tensors(ψ))) -end - -default_norm_cache(ψ::ITensorNetwork) = BeliefPropagationCache(QuadraticFormNetwork(ψ)) - -function ITensorNetwork( - ψ_vidal::VidalITensorNetwork; (cache!) = nothing, update_gauge = false, update_kwargs... - ) - if update_gauge - ψ_vidal = update(ψ_vidal; update_kwargs...) - end - - ψ = copy(site_tensors(ψ_vidal)) - - for e in edges(ψ) - vsrc, vdst = src(e), dst(e) - root_S = ITensorsExtensions.sqrt_diag(bond_tensor(ψ_vidal, e)) - setindex_preserve_graph!(ψ, noprime(root_S * ψ[vsrc]), vsrc) - setindex_preserve_graph!(ψ, noprime(root_S * ψ[vdst]), vdst) - end - - if !isnothing(cache!) - bp_cache = default_norm_cache(ψ) - mts = messages(bp_cache) - - for e in edges(ψ) - vsrc, vdst = src(e), dst(e) - pe = quotientedge(bp_cache, (vsrc, "bra") => (vdst, "bra")) - set!(mts, pe, copy(ITensor[dense(bond_tensor(ψ_vidal, e))])) - set!(mts, reverse(pe), copy(ITensor[dense(bond_tensor(ψ_vidal, e))])) - end - - bp_cache = set_messages(bp_cache, mts) - cache![] = bp_cache - end - - return ψ -end - -# Use an ITensorNetwork ψ, its bond tensors and belief propagation cache to put ψ into the vidal gauge, return the bond tensors and updated_ψ. -function vidalitensornetwork_preserve_cache( - ψ::ITensorNetwork; - cache = default_norm_cache(ψ), - bond_tensors = default_bond_tensors, - message_cutoff = 10 * eps(real(scalartype(ψ))), - regularization = 10 * eps(real(scalartype(ψ))), - edges = NamedGraphs.edges(ψ), - svd_kwargs... - ) - ψ_vidal_site_tensors = copy(ψ) - ψ_vidal_bond_tensors = bond_tensors(ψ) - - for e in edges - vsrc, vdst = src(e), dst(e) - ψvsrc, ψvdst = ψ_vidal_site_tensors[vsrc], ψ_vidal_site_tensors[vdst] - - pe = quotientedge(cache, (vsrc, "bra") => (vdst, "bra")) - edge_ind = commoninds(ψvsrc, ψvdst) - edge_ind_sim = sim(edge_ind) - - X_D, X_U = - eigen(only(message(cache, pe)); ishermitian = true, cutoff = message_cutoff) - Y_D, Y_U = eigen( - only(message(cache, reverse(pe))); ishermitian = true, - cutoff = message_cutoff - ) - X_D, Y_D = ITensorsExtensions.map_diag(x -> x + regularization, X_D), - ITensorsExtensions.map_diag(x -> x + regularization, Y_D) - - rootX_D, rootY_D = - ITensorsExtensions.sqrt_diag(X_D), ITensorsExtensions.sqrt_diag(Y_D) - inv_rootX_D, inv_rootY_D = ITensorsExtensions.invsqrt_diag(X_D), - ITensorsExtensions.invsqrt_diag(Y_D) - rootX = X_U * rootX_D * prime(dag(X_U)) - rootY = Y_U * rootY_D * prime(dag(Y_U)) - inv_rootX = X_U * inv_rootX_D * prime(dag(X_U)) - inv_rootY = Y_U * inv_rootY_D * prime(dag(Y_U)) - - ψvsrc, ψvdst = noprime(ψvsrc * inv_rootX), noprime(ψvdst * inv_rootY) - - Ce = rootX - Ce = Ce * replaceinds(rootY, edge_ind, edge_ind_sim) - - U, S, V = svd(Ce, edge_ind; svd_kwargs...) - - new_edge_ind = Index[Index(dim(commoninds(S, U)), tags(first(edge_ind)))] - - ψvsrc = replaceinds(ψvsrc * U, commoninds(S, U), new_edge_ind) - ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_sim) - ψvdst = replaceinds(ψvdst * V, commoninds(V, S), new_edge_ind) - - setindex_preserve_graph!(ψ_vidal_site_tensors, ψvsrc, vsrc) - setindex_preserve_graph!(ψ_vidal_site_tensors, ψvdst, vdst) - - S = replaceinds( - S, - [commoninds(S, U)..., commoninds(S, V)...] => - [new_edge_ind..., prime(new_edge_ind)...] - ) - ψ_vidal_bond_tensors[e] = S - end - - return VidalITensorNetwork(ψ_vidal_site_tensors, ψ_vidal_bond_tensors) -end - -function VidalITensorNetwork( - ψ::ITensorNetwork; - (cache!) = nothing, - update_cache = isnothing(cache!), - cache_update_kwargs = (;), - kwargs... - ) - if isnothing(cache!) - cache! = Ref(default_norm_cache(ψ)) - end - - if update_cache - cache![] = update(cache![]; cache_update_kwargs...) - end - - return vidalitensornetwork_preserve_cache(ψ; cache = cache![], kwargs...) -end - -function update(ψ::VidalITensorNetwork; kwargs...) - return VidalITensorNetwork(ITensorNetwork(ψ; update_gauge = false); kwargs...) -end - -# Function to construct the 'isometry' of a state in the Vidal Gauge on the given edge -function vidal_gauge_isometry(ψ::VidalITensorNetwork, edge) - vsrc, vdst = src(edge), dst(edge) - ψ_vsrc = copy(ψ[vsrc]) - - for vn in setdiff(neighbors(ψ, vsrc), [vdst]) - ψ_vsrc = noprime(ψ_vsrc * bond_tensor(ψ, vn => vsrc)) - end - - ψ_vsrcdag = dag(ψ_vsrc) - ψ_vsrcdag = - replaceind(ψ_vsrcdag, commonind(ψ_vsrc, ψ[vdst]), commonind(ψ_vsrc, ψ[vdst])') - - return ψ_vsrcdag * ψ_vsrc -end - -function vidal_gauge_isometries(ψ::VidalITensorNetwork, edges::Vector) - return Dict([e => vidal_gauge_isometry(ψ, e) for e in edges]) -end - -function vidal_gauge_isometries(ψ::VidalITensorNetwork) - return vidal_gauge_isometries( - ψ, vcat(NamedGraphs.edges(ψ), reverse.(NamedGraphs.edges(ψ))) - ) -end - -# Function to measure the 'distance' of a state from the Vidal Gauge -function gauge_error(ψ::VidalITensorNetwork) - f = 0 - isometries = vidal_gauge_isometries(ψ) - for e in keys(isometries) - lhs = isometries[e] - f += message_diff(ITensor[lhs], ITensor[denseblocks(delta(inds(lhs)))]) - end - - return f / (length(keys(isometries))) -end diff --git a/test/Project.toml b/test/Project.toml index 786d3178..fe1222b5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -46,7 +46,7 @@ Glob = "1.3.1" Graphs = "1.12" GraphsFlows = "0.1.1" ITensorMPS = "0.3.6" -ITensorNetworks = "0.17" +ITensorNetworks = "0.18" ITensorPkgSkeleton = "0.3.42" ITensors = "0.7, 0.8, 0.9" KrylovKit = "0.8, 0.9, 0.10" diff --git a/test/test_apply.jl b/test/test_apply.jl index 86246f4f..70d447e1 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,8 +1,8 @@ @eval module $(gensym()) using Compat: Compat using Graphs: vertices -using ITensorNetworks: BeliefPropagationCache, ITensorNetwork, VidalITensorNetwork, apply, - environment, norm_sqr_network, random_tensornetwork, siteinds, update +using ITensorNetworks: BeliefPropagationCache, apply, environment, norm_sqr_network, + random_tensornetwork, siteinds, update using ITensors: ITensors, Algorithm, ITensor, inner, op using NamedGraphs.NamedGraphGenerators: named_grid using SplitApplyCombine: group @@ -23,7 +23,6 @@ using Test: @test, @testset bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter = 20) envsSBP = environment(bp_cache, [(v1, "bra"), (v1, "ket"), (v2, "bra"), (v2, "ket")]) - ψv = VidalITensorNetwork(ψ; cache_update_kwargs = (; maxiter = 20)) #This grouping will correspond to calculating the environments exactly (each column of the grid is a partition) bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1][1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter = 20) @@ -49,8 +48,6 @@ using Test: @test, @testset envisposdef = true, callback ) - ψOv = apply(o, ψv; maxdim = χ, normalize = true) - ψOVidal_symm = ITensorNetwork(ψOv) ψOGBP = apply( o, ψ; @@ -66,11 +63,6 @@ using Test: @test, @testset inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOSBP, ψOSBP; alg = inner_alg) ) - fVidal = - inner(ψOVidal_symm, ψOexact; alg = inner_alg) / sqrt( - inner(ψOexact, ψOexact; alg = inner_alg) * - inner(ψOVidal_symm, ψOVidal_symm; alg = inner_alg) - ) fGBP = inner(ψOGBP, ψOexact; alg = inner_alg) / sqrt( @@ -79,7 +71,6 @@ using Test: @test, @testset ) @test !iszero(truncerr) @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP)) - @test isapprox(real(fSBP * conj(fSBP)), real(fVidal * conj(fVidal)); atol = 1.0e-3) end end end diff --git a/test/test_gauging.jl b/test/test_gauging.jl deleted file mode 100644 index 2cf171b0..00000000 --- a/test/test_gauging.jl +++ /dev/null @@ -1,46 +0,0 @@ -@eval module $(gensym()) -using Compat: Compat -using ITensorNetworks: BeliefPropagationCache, ITensorNetwork, VidalITensorNetwork, - gauge_error, messages, random_tensornetwork, siteinds, update -using ITensors.NDTensors: Algorithm, vector -using ITensors: diag_itensor, inds, inner -using LinearAlgebra: diag -using NamedGraphs.NamedGraphGenerators: named_grid -using StableRNGs: StableRNG -using TensorOperations: TensorOperations -using Test: @test, @testset - -@testset "gauging" begin - n = 3 - dims = (n, n) - g = named_grid(dims) - s = siteinds("S=1/2", g) - χ = 6 - - rng = StableRNG(1234) - ψ = random_tensornetwork(rng, s; link_space = χ) - - # Move directly to vidal gauge - ψ_vidal = VidalITensorNetwork( - ψ; cache_update_kwargs = (; alg = "bp", maxiter = 30, verbose = true) - ) - @test gauge_error(ψ_vidal) < 1.0e-8 - - # Move to symmetric gauge - cache_ref = Ref{BeliefPropagationCache}() - ψ_symm = ITensorNetwork(ψ_vidal; (cache!) = cache_ref) - bp_cache = cache_ref[] - - # Test we just did a gauge transform and didn't change the overall network - @test inner(ψ_symm, ψ; alg = "exact") / - sqrt(inner(ψ_symm, ψ_symm; alg = "exact") * inner(ψ, ψ; alg = "exact")) ≈ 1.0 atol = - 1.0e-8 - - #Test all message tensors are approximately diagonal even when we keep running BP - bp_cache = update(bp_cache; maxiter = 10) - for m_e in values(messages(bp_cache)) - @test diag_itensor(vector(diag(only(m_e))), inds(only(m_e))) ≈ only(m_e) atol = - 1.0e-8 - end -end -end