Merge branch 'master' into lb/ctmrgenv_fix

.github/workflows/CI.yml

@@ -21,7 +21,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8' # LTS version
+          - '1.9'
           - '1' # automatically expands to the latest stable 1.x release of Julia
         os:
           - ubuntu-latest

Project.toml

@@ -23,18 +23,18 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 Accessors = "0.1"
 ChainRulesCore = "1.0"
 Compat = "3.46, 4.2"
-KrylovKit = "0.6, 0.7, 0.8"
+KrylovKit = "0.8"
 LinearAlgebra = "1"
 MPSKit = "0.11"
 OptimKit = "0.3"
 Printf = "1"
 Random = "1"
 Statistics = "1"
-TensorKit = "0.12.4"
+TensorKit = "0.12.5"
 TensorOperations = "4"
 VectorInterface = "0.4"
 Zygote = "0.6"
-julia = "1.8"
+julia = "1.9"
 
 [extras]
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"

README.md

@@ -36,19 +36,12 @@ For example, in order to obtain the groundstate of the 2D Heisenberg model, we c
 using TensorKit, PEPSKit, KrylovKit, OptimKit
 
 # constructing the Hamiltonian:
-Jx, Jy, Jz = (-1, 1, -1) # sublattice rotation to obtain single-site unit cell
-physical_space = ComplexSpace(2)
-T = ComplexF64
-σx = TensorMap(T[0 1; 1 0], physical_space, physical_space)
-σy = TensorMap(T[0 im; -im 0], physical_space, physical_space)
-σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space)
-H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz)
-Heisenberg_hamiltonian = NLocalOperator{NearestNeighbor}(H / 4)
+H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) # sublattice rotation to obtain single-site unit cell
 
 # configuring the parameters
 D = 2
 chi = 20
-ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1)
+ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
@@ -60,7 +53,7 @@ opt_alg = PEPSOptimize(;
 # ground state search
 state = InfinitePEPS(2, D)
 ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
-result = fixedpoint(state, Heisenberg_hamiltonian, opt_alg, ctm)
+result = fixedpoint(state, H, opt_alg, ctm)
 
 @show result.E # -0.6625...
 ```

docs/src/index.md

@@ -21,19 +21,12 @@ For example, in order to obtain the groundstate of the 2D Heisenberg model, we c
 using TensorKit, PEPSKit, KrylovKit, OptimKit
 
 # constructing the Hamiltonian:
-Jx, Jy, Jz = (-1, 1, -1) # sublattice rotation to obtain single-site unit cell
-physical_space = ComplexSpace(2)
-T = ComplexF64
-σx = TensorMap(T[0 1; 1 0], physical_space, physical_space)
-σy = TensorMap(T[0 im; -im 0], physical_space, physical_space)
-σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space)
-H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz)
-Heisenberg_hamiltonian = NLocalOperator{NearestNeighbor}(H / 4)
+H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) # sublattice rotation to obtain single-site unit cell
 
 # configuring the parameters
 D = 2
 chi = 20
-ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1)
+ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
@@ -44,8 +37,13 @@ opt_alg = PEPSOptimize(;
 
 # ground state search
 state = InfinitePEPS(2, D)
+<<<<<<< HEAD
 ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
 result = fixedpoint(state, Heisenberg_hamiltonian, opt_alg, ctm)
+=======
+ctm = leading_boundary(CTMRGEnv(state; Venv=ComplexSpace(chi)), state, ctm_alg)
+result = fixedpoint(state, H, opt_alg, ctm)
+>>>>>>> master
 
 @show result.E # -0.6625...
 ```

examples/boundary_mps.jl

@@ -32,9 +32,7 @@ mps, envs, ϵ = leading_boundary(mps, T, VUMPS())
 N = abs(prod(expectation_value(mps, T)))
 
 # This can be compared to the result obtained using the CTMRG algorithm
-ctm = leading_boundary(
-    peps, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(peps, ComplexSpace(20))
-)
+ctm = leading_boundary(peps, CTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20)))
 N´ = abs(norm(peps, ctm))
 
 @show abs(N - N´) / N
@@ -55,9 +53,7 @@ mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2))
 mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS())
 N2 = abs(prod(expectation_value(mps2, T2)))
 
-ctm2 = leading_boundary(
-    peps2, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(peps2, ComplexSpace(20))
-)
+ctm2 = leading_boundary(peps2, CTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20)))
 N2´ = abs(norm(peps2, ctm2))
 
 @show abs(N2 - N2´) / N2

examples/heisenberg.jl

@@ -11,9 +11,9 @@ H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1)
 # Parameters
 χbond = 2
 χenv = 20
-ctmalg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1)
-alg = PEPSOptimize(;
-    boundary_alg=ctmalg,
+ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(χenv))
+opt_alg = PEPSOptimize(;
+    boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
     gradient_alg=GMRES(; tol=1e-6, maxiter=100),
     reuse_env=true,
@@ -28,5 +28,5 @@ alg = PEPSOptimize(;
 # Of course there is a noticable bias for small χbond and χenv.
 ψ₀ = InfinitePEPS(2, χbond)
 env₀ = leading_boundary(CTMRGEnv(ψ₀, ℂ^χenv), ψ₀, ctmalg)
-result = fixedpoint(ψ₀, H, alg, env₀)
+result = fixedpoint(ψ₀, H, opt_alg, env₀)
 @show result.E

src/PEPSKit.jl

@@ -9,7 +9,7 @@ using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations
 using ChainRulesCore, Zygote
 
 include("utility/util.jl")
-include("utility/eigsolve.jl")
+include("utility/svd.jl")
 include("utility/rotations.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
@@ -59,12 +59,14 @@ module Defaults
     const fpgrad_tol = 1e-6
 end
 
-export CTMRG, CTMRGEnv
+export SVDAdjoint, IterSVD, NonTruncSVDAdjoint
+export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv
 export LocalOperator
 export expectation_value, costfun
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolve
 export fixedpoint
+
 export InfinitePEPS, InfiniteTransferPEPS
 export InfinitePEPO, InfiniteTransferPEPO
 export initializeMPS, initializePEPS

src/algorithms/ctmrg.jl

@@ -1,24 +1,60 @@
+"""
+    FixedSpaceTruncation <: TensorKit.TruncationScheme
+
+CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD
+is performed fixed. Since different environment directions and unit cell entries might
+have different spaces, this truncation style is different from `TruncationSpace`.
+"""
+struct FixedSpaceTruncation <: TensorKit.TruncationScheme end
+
+# TODO: add option for different projector styles (half-infinite, full-infinite, etc.)
+"""
+    struct ProjectorAlg{S}(; svd_alg = TensorKit.SVD(), trscheme = TensorKit.notrunc(),
+                           fixedspace = false, verbosity = 0)
+
+Algorithm struct collecting all projector related parameters. The truncation scheme has to be
+a `TensorKit.TruncationScheme`, and some SVD algorithms might have further restrictions on what
+kind of truncation scheme can be used. If `fixedspace` is true, the truncation scheme is set to
+`truncspace(V)` where `V` is the environment bond space, adjusted to the corresponding
+environment direction/unit cell entry.
+"""
+@kwdef struct ProjectorAlg{S<:SVDAdjoint,T}
+    svd_alg::S = SVDAdjoint()
+    trscheme::T = FixedSpaceTruncation()
+    verbosity::Int = 0
+end
+
 # TODO: add abstract Algorithm type?
 """
-    struct CTMRG(; trscheme = TensorKit.notrunc(), tol = Defaults.ctmrg_tol,
-                 maxiter = Defaults.ctmrg_maxiter, miniter = Defaults.ctmrg_miniter,
-                 verbosity = 0, fixedspace = false)
+    CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter,
+          miniter=Defaults.ctmrg_miniter, verbosity=0,
+          svd_alg=TensorKit.SVD(), trscheme=FixedSpaceTruncation())
 
 Algorithm struct that represents the CTMRG algorithm for contracting infinite PEPS.
-The projector bond dimensions are set via `trscheme` which controls the truncation
-properties inside of `TensorKit.tsvd`. Each CTMRG run is converged up to `tol`
-where the singular value convergence of the corners as well as the norm is checked.
-The maximal and minimal number of CTMRG iterations is set with `maxiter` and `miniter`.
-Different levels of output information are printed depending on `verbosity` (0, 1 or 2).
-Regardless of the truncation scheme, the space can be kept fixed with `fixedspace`.
+Each CTMRG run is converged up to `tol` where the singular value convergence of the
+corners as well as the norm is checked. The maximal and minimal number of CTMRG iterations
+is set with `maxiter` and `miniter`. Different levels of output information are printed
+depending on `verbosity` (0, 1 or 2). The projectors are computed from `svd_alg` SVDs
+where the truncation scheme is set via `trscheme`.
 """
-@kwdef struct CTMRG
-    trscheme::TruncationScheme = TensorKit.notrunc()
-    tol::Float64 = Defaults.ctmrg_tol
-    maxiter::Int = Defaults.ctmrg_maxiter
-    miniter::Int = Defaults.ctmrg_miniter
-    verbosity::Int = 0
-    fixedspace::Bool = false
+struct CTMRG
+    tol::Float64
+    maxiter::Int
+    miniter::Int
+    verbosity::Int
+    projector_alg::ProjectorAlg
+end
+function CTMRG(;
+    tol=Defaults.ctmrg_tol,
+    maxiter=Defaults.ctmrg_maxiter,
+    miniter=Defaults.ctmrg_miniter,
+    verbosity=0,
+    svd_alg=SVDAdjoint(),
+    trscheme=FixedSpaceTruncation(),
+)
+    return CTMRG(
+        tol, maxiter, miniter, verbosity, ProjectorAlg(; svd_alg, trscheme, verbosity)
+    )
 end
 
 """
@@ -39,6 +75,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
 
     for i in 1:(alg.maxiter)
         env, ϵ = ctmrg_iter(state, env, alg)  # Grow and renormalize in all 4 directions
+
         conv_condition, normold, CSold, TSold, ϵ = ignore_derivatives() do
             # Compute convergence criteria and take max (TODO: How should we handle logging all of this?)
             Δϵ = abs((ϵold - ϵ) / ϵold)
@@ -86,9 +123,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
     end
 
     # Do one final iteration that does not change the spaces
-    alg_fixed = CTMRG(;
-        alg.trscheme, alg.tol, alg.maxiter, alg.miniter, alg.verbosity, fixedspace=true
-    )
+    alg_fixed = @set alg.projector_alg.trscheme = FixedSpaceTruncation()
     env′, = ctmrg_iter(state, env, alg_fixed)
     envfix = gauge_fix(env, env′)
     check_elementwise_convergence(env, envfix; atol=alg.tol^(1 / 2)) ||
@@ -288,7 +323,7 @@ function ctmrg_iter(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
     ϵ = 0.0
 
     for _ in 1:4
-        env, _, _, ϵ₀ = left_move(state, env, alg)
+        env, _, _, ϵ₀ = left_move(state, env, alg.projector_alg)
         state = rotate_north(state, EAST)
         env = rotate_north(env, EAST)
         ϵ = max(ϵ, ϵ₀)
@@ -303,14 +338,15 @@ end
 Grow, project and renormalize the environment `env` in west direction.
 Return the updated environment as well as the projectors and truncation error.
 """
-function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
+function left_move(state, env::CTMRGEnv{C,T}, alg::ProjectorAlg) where {C,T}
     corners::typeof(env.corners) = copy(env.corners)
     edges::typeof(env.edges) = copy(env.edges)
     ϵ = 0.0
     Pleft, Pright = Zygote.Buffer.(projector_type(T, size(state)))  # Use Zygote.Buffer instead of @diffset to avoid ZeroTangent errors in _setindex
 
     for col in 1:size(state, 2)
         cprev = _prev(col, size(state, 2))
+        cnext = _next(col, size(state, 2))
 
         # Compute projectors
         for row in 1:size(state, 1)
@@ -332,15 +368,15 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
             )
 
             # SVD half-infinite environment
-            trscheme = if alg.fixedspace == true
-                truncspace(space(env.edges[WEST, row, col], 1))
+            trscheme = if alg.trscheme isa FixedSpaceTruncation
+                truncspace(space(env.edges[WEST, row, cnext], 1))
             else
                 alg.trscheme
             end
             @autoopt @tensor QQ[χ_EB D_EBabove D_EBbelow; χ_ET D_ETabove D_ETbelow] :=
                 Q_sw[χ_EB D_EBabove D_EBbelow; χ D1 D2] *
                 Q_nw[χ D1 D2; χ_ET D_ETabove D_ETbelow]
-            U, S, V, ϵ_local = tsvd!(QQ; trunc=trscheme, alg=TensorKit.SVD())
+            U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, alg.svd_alg; trunc=trscheme)
             ϵ = max(ϵ, ϵ_local / norm(S))
             # TODO: check if we can just normalize enlarged corners s.t. trunc behaves a bit better