[WIP] Environment initialization for CTMRG

src/PEPSKit.jl

@@ -2,7 +2,7 @@ module PEPSKit
 
 using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Compat
-using Accessors: @set, @reset
+using Accessors: @set, @reset, @insert
 using VectorInterface
 import VectorInterface as VI
 
@@ -64,6 +64,7 @@ include("algorithms/ctmrg/projectors.jl")
 include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
+include("algorithms/ctmrg/initialization.jl")
 
 include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
@@ -87,6 +88,8 @@ using .Defaults: set_scheduler!
 export set_scheduler!
 export SVDAdjoint, FullSVDReverseRule, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
+export initialize_environment,
+    RandomInitialization, ProductStateInitialization, ApplicationInitialization
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export LocalOperator, physicalspace

src/algorithms/ctmrg/ctmrg.jl

@@ -88,6 +88,7 @@ containing the following fields:
 
 * `truncation_error` : Last (maximal) SVD truncation error of the CTMRG projectors.
 * `condition_number` : Last (maximal) condition number of the enlarged CTMRG environment.
+* `convergence_error` : Convergence error of the CTMRG algorithm at termination.
 
 In case the `alg` is a `SimultaneousCTMRG`, the last SVD will also be returned:
 
@@ -120,6 +121,7 @@ function leading_boundary(
         for iter in 1:(alg.maxiter)
             env, info = ctmrg_iteration(network, env, alg)  # Grow and renormalize in all 4 directions
             η, CS, TS = calc_convergence(env, CS, TS)
+            info = @insert info.convergence_error = η
 
             if η ≤ alg.tol && iter ≥ alg.miniter
                 ctmrg_logfinish!(log, iter, η, network, env)

src/algorithms/ctmrg/initialization.jl

@@ -0,0 +1,59 @@
+abstract type InitializationStyle end
+struct ProductStateInitialization <: InitializationStyle end
+struct RandomInitialization{F} <: InitializationStyle
+    f::F
+    RandomInitialization(f::F = randn) where {F} = new{F}(f)
+end
+struct ApplicationInitialization <: InitializationStyle end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        alg::RandomInitialization,
+        virtual_spaces... = oneunit(spacetype(n)),
+    )
+    return CTMRGEnv(alg.f, elt, n, virtual_spaces...)
+end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::ProductStateInitialization,
+        virtual_spaces... = oneunit(spacetype(n)),
+    )
+    i = one(sectortype(n))
+    env = CTMRGEnv(ones, elt, n, virtual_spaces...)
+    for (dir, r, c) in Iterators.product(axes(env)...)
+        @assert i in blocksectors(env.corners[dir, r, c])
+        for (c, b) in blocks(env.corners[dir, r, c])
+            b .= 0
+            c == i && (b[1, 1] = 1)
+        end
+    end
+    return env
+end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::ApplicationInitialization,
+        trscheme::TruncationScheme;
+        boundary_alg = (;
+            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1,
+        )
+    )
+    boundary_alg = (; boundary_alg..., trscheme) # merge trscheme with optional alg definition
+    env = initialize_environment(elt, n, ProductStateInitialization())
+    env, = leading_boundary(env, n; boundary_alg...)
+    return env
+end
+
+function initialize_environment(n::InfiniteSquareNetwork, args...; kwargs...)
+    return initialize_environment(scalartype(n), n, args...; kwargs...)
+end
+function initialize_environment(A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
+    return initialize_environment(scalartype(A), A, args...; kwargs...)
+end
+function initialize_environment(elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
+    return initialize_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
+end

src/networks/infinitesquarenetwork.jl

@@ -55,6 +55,7 @@ end
 
 ## Spaces
 
+TensorKit.spacetype(::Type{T}) where {T <: InfiniteSquareNetwork} = spacetype(eltype(T))
 virtualspace(n::InfiniteSquareNetwork, r::Int, c::Int, dir) = virtualspace(n[r, c], dir)
 
 ## Vector interface

test/ctmrg/initialization.jl

@@ -0,0 +1,51 @@
+using Test
+using TensorKit
+using PEPSKit
+using Random
+
+using MPSKitModels: classical_ising
+
+sd = 12345
+
+# toggle symmetry, but same issue for both
+symmetries = [Z2Irrep, Trivial]
+
+χ = 20
+tol = 1.0e-4
+maxiter = 1000
+verbosity = 2
+trscheme = FixedSpaceTruncation()
+boundary_alg = (;
+    alg = :simultaneous,
+    tol,
+    verbosity,
+    trscheme,
+    maxiter,
+)
+
+@testset "CTMRG environment initialization for critical ising with $S symmetry (#255)" for S in symmetries
+    # initialize
+    T = classical_ising(S)
+    O = T[1]
+    n = InfinitePartitionFunction([O O; O O])
+    Venv = S == Z2Irrep ? Z2Space(0 => χ / 2, 1 => χ / 2) : ℂ^χ
+    P = space(O, 2)
+
+    # random, doesn't converge
+    Random.seed!(sd)
+    env0_rand = initialize_environment(n, RandomInitialization(), Venv)
+    env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
+    @test_broken info.convergence_error ≤ tol
+
+    # embedded product state, converges
+    Random.seed!(sd)
+    env0_prod = initialize_environment(n, ProductStateInitialization(), Venv)
+    env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
+    # grown product state, converges
+    Random.seed!(sd)
+    env0_appl = initialize_environment(InfiniteSquareNetwork(n), ApplicationInitialization(), truncdim(χ))
+    env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+end

test/runtests.jl

@@ -35,6 +35,9 @@ end
         @time @safetestset "correlation length" begin
             include("ctmrg/correlation_length.jl")
         end
+        @time @safetestset "initialization" begin
+            include("ctmrg/initialization.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "GRADIENTS"
         @time @safetestset "CTMRG gradients" begin