Implement full-infinite projector algorithm

README.md

@@ -44,7 +44,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t
 # configuring the parameters
 D = 2
 chi = 20
-ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi))
+ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),

docs/src/index.md

@@ -26,7 +26,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t
 # configuring the parameters
 D = 2
 chi = 20
-ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi))
+ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),

examples/boundary_mps.jl

@@ -32,7 +32,9 @@ 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), CTMRGEnv(peps, ComplexSpace(20)))
+ctm = leading_boundary(
+    peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))
+)
 N´ = abs(norm(peps, ctm))
 
 @show abs(N - N´) / N
@@ -53,7 +55,9 @@ 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), CTMRGEnv(peps2, ComplexSpace(20)))
+ctm2 = leading_boundary(
+    peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))
+)
 N2´ = abs(norm(peps2, ctm2))
 
 @show abs(N2 - N2´) / N2

examples/heisenberg.jl

@@ -11,7 +11,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1)
 # Parameters
 χbond = 2
 χenv = 20
-ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=2)
+ctm_alg = SimultaneousCTMRG(; tol=1e-10, verbosity=2)
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),

examples/hubbard_su.jl

@@ -1,5 +1,4 @@
 using Test
-using Printf
 using Random
 using PEPSKit
 using TensorKit
@@ -14,15 +13,16 @@ if symm == Trivial
 else
     error("Not implemented")
 end
-
 peps = InfiniteWeightPEPS(rand, Float64, Pspace, Vspace; unitcell=(N1, N2))
+
 # normalize vertex tensors
 for ind in CartesianIndices(peps.vertices)
     peps.vertices[ind] /= norm(peps.vertices[ind], Inf)
 end
+
 # Hubbard model Hamiltonian at half-filling
-t, U = 1.0, 6.0
-ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(N1, N2); t=t, U=U, mu=U / 2)
+t, U = 1, 6
+ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(N1, N2); t, U, mu=U / 2)
 
 # simple update
 dts = [1e-2, 1e-3, 4e-4, 1e-4]
@@ -31,32 +31,29 @@ maxiter = 10000
 for (n, (dt, tol)) in enumerate(zip(dts, tols))
     trscheme = truncerr(1e-10) & truncdim(Dcut)
     alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-    result = simpleupdate(peps, ham, alg; bipartite=false)
-    global peps = result[1]
+    peps, = simpleupdate(peps, ham, alg; bipartite=false)
 end
 
-# absort weight into site tensors
+# absorb weight into site tensors
 peps = InfinitePEPS(peps)
+
 # CTMRG
 χenv0, χenv = 6, 20
 Espace = Vect[fℤ₂](0 => χenv0 / 2, 1 => χenv0 / 2)
 envs = CTMRGEnv(randn, Float64, peps, Espace)
 for χ in [χenv0, χenv]
-    trscheme = truncerr(1e-9) & truncdim(χ)
-    ctm_alg = CTMRG(;
-        maxiter=40, tol=1e-10, verbosity=3, trscheme=trscheme, ctmrgscheme=:sequential
-    )
-    global envs = leading_boundary(envs, peps, ctm_alg)
+    ctm_alg = SequentialCTMRG(; maxiter=300, tol=1e-7)
+    envs = leading_boundary(envs, peps, ctm_alg)
 end
 
-"""
-Benchmark values of the ground state energy from
-Qin, M., Shi, H., & Zhang, S. (2016). Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method. Physical Review B, 94(8), 085103.
-"""
-# measure physical quantities
-E = costfun(peps, envs, ham) / (N1 * N2)
+# Benchmark values of the ground state energy from
+# Qin, M., Shi, H., & Zhang, S. (2016). Benchmark study of the two-dimensional Hubbard
+# model with auxiliary-field quantum Monte Carlo method. Physical Review B, 94(8), 085103.
 Es_exact = Dict(0 => -1.62, 2 => -0.176, 4 => 0.8603, 6 => -0.6567, 8 => -0.5243)
 E_exact = Es_exact[U] - U / 2
-@info @sprintf("Energy           = %.8f\n", E)
-@info @sprintf("Benchmark energy = %.8f\n", E_exact)
-@test isapprox(E, E_exact; atol=1e-2)
+
+# measure energy
+E = costfun(peps, envs, ham) / (N1 * N2)
+@info "Energy           = $E"
+@info "Benchmark energy = $E_exact"
+@test isapprox(E, E_exact; atol=5e-2)

src/PEPSKit.jl

@@ -3,7 +3,7 @@ module PEPSKit
 using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Base: @kwdef
 using Compat
-using Accessors
+using Accessors: @set
 using VectorInterface
 using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations
 using ChainRulesCore, Zygote
@@ -43,6 +43,9 @@ include("algorithms/contractions/ctmrg_contractions.jl")
 
 include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")
+include("algorithms/ctmrg/projectors.jl")
+include("algorithms/ctmrg/simultaneous.jl")
+include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 
 include("algorithms/time_evolution/gatetools.jl")
@@ -61,12 +64,16 @@ include("utility/symmetrization.jl")
         const ctmrg_tol = 1e-8
         const fpgrad_maxiter = 30
         const fpgrad_tol = 1e-6
-        const ctmrgscheme = :simultaneous
         const reuse_env = true
         const trscheme = FixedSpaceTruncation()
         const fwd_alg = TensorKit.SDD()
         const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
         const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
+        const projector_alg_type = HalfInfiniteProjector
+        const projector_alg = projector_alg_type(svd_alg, trscheme, 2)
+        const ctmrg_alg = SimultaneousCTMRG(
+            ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
+        )
         const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
         const gradient_linsolver = KrylovKit.BiCGStab(;
             maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
@@ -83,12 +90,14 @@ Module containing default values that represent typical algorithm parameters.
 - `ctmrg_tol`: Tolerance checking singular value and norm convergence
 - `fpgrad_maxiter`: Maximal number of iterations for computing the CTMRG fixed-point gradient
 - `fpgrad_tol`: Convergence tolerance for the fixed-point gradient iteration
-- `ctmrgscheme`: Scheme for growing, projecting and renormalizing CTMRG environments
 - `reuse_env`: If `true`, the current optimization step is initialized on the previous environment
 - `trscheme`: Truncation scheme for SVDs and other decompositions
 - `fwd_alg`: SVD algorithm that is used in the forward pass
 - `rrule_alg`: Reverse-rule for differentiating that SVD
 - `svd_alg`: Combination of `fwd_alg` and `rrule_alg`
+- `projector_alg_type`: Default type of projector algorithm
+- `projector_alg`: Algorithm to compute CTMRG projectors
+- `ctmrg_alg`: Algorithm for performing CTMRG runs
 - `optimizer`: Optimization algorithm for PEPS ground-state optimization
 - `gradient_linsolver`: Default linear solver for the `LinSolver` gradient algorithm
 - `iterscheme`: Scheme for differentiating one CTMRG iteration
@@ -97,19 +106,28 @@ Module containing default values that represent typical algorithm parameters.
 """
 module Defaults
     using TensorKit, KrylovKit, OptimKit, OhMyThreads
-    using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint
+    using PEPSKit:
+        LinSolver,
+        FixedSpaceTruncation,
+        SVDAdjoint,
+        HalfInfiniteProjector,
+        SimultaneousCTMRG
+    const ctmrg_tol = 1e-8
     const ctmrg_maxiter = 100
     const ctmrg_miniter = 4
-    const ctmrg_tol = 1e-8
     const fpgrad_maxiter = 30
     const fpgrad_tol = 1e-6
-    const ctmrgscheme = :simultaneous
     const sparse = false
     const reuse_env = true
     const trscheme = FixedSpaceTruncation()
     const fwd_alg = TensorKit.SDD()
     const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
     const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
+    const projector_alg_type = HalfInfiniteProjector
+    const projector_alg = projector_alg_type(svd_alg, trscheme, 2)
+    const ctmrg_alg = SimultaneousCTMRG(
+        ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
+    )
     const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
     const gradient_linsolver = KrylovKit.BiCGStab(;
         maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
@@ -165,9 +183,10 @@ end
 using .Defaults: set_scheduler!
 export set_scheduler!
 export SVDAdjoint, IterSVD, NonTruncSVDAdjoint
-export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length
+export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
+export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector
 export LocalOperator
-export expectation_value, costfun, product_peps
+export expectation_value, costfun, product_peps, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver
 export fixedpoint

src/algorithms/contractions/ctmrg_contractions.jl

@@ -210,13 +210,40 @@
 end
 
 """
-    halfinfinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3})
-    halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4,
-                             ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
-    halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x,
-                             ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
-    halfinfinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4,
-                             ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
+    contract_projectors(U, S, V, Q, Q_next)
+
+Compute projectors based on a SVD of `Q * Q_next`, where the inverse square root
+`isqS` of the singular values is computed.
+
+Left projector:
+```
+    -- |~~~~~~| -- |~~|
+       |Q_next|    |V'| -- isqS --
+    == |~~~~~~| == |~~|
+```
+
+Right projector:
+```
+               |~~| -- |~~~| --
+    -- isqS -- |U'|    | Q |
+               |~~| == |~~~| ==
+```
+"""
+function contract_projectors(U, S, V, Q, Q_next)
+    isqS = sdiag_pow(S, -0.5)
+    P_left = Q_next * V' * isqS  # use * to respect fermionic case
+    P_right = isqS * U' * Q
+    return P_left, P_right
+end
+
+"""
+    half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3})
+    half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4,
+                              ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
+    half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x,
+                              ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
+    half_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4,
+                              ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor}
 
 Contract two quadrants (enlarged corners) to form a half-infinite environment.
 
@@ -236,27 +263,26 @@
      |       ||          ||           |
 ```
 
-Alternatively, contract environment with a vector `x` acting on it
+Alternatively, contract the environment with a vector `x` acting on it
 
 ```
     C_1 --  E_2      --  E_3      -- C_2
      |       ||          ||           | 
     E_1 == ket_bra_1 == ket_bra_2 == E_4
      |       ||          ||           |
                          [~~~~~~x~~~~~~]
-                         ||           |
 ```
 
 or contract the adjoint environment with `x`, e.g. as needed for iterative solvers.
 """
-function halfinfinite_environment(
+function half_infinite_environment(
     quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}
 ) where {S}
     return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] :=
         quadrant1[χ_in D_inabove D_inbelow; χ D1 D2] *
         quadrant2[χ D1 D2; χ_out D_outabove D_outbelow]
 end
-function halfinfinite_environment(
+function half_infinite_environment(
     C_1, C_2, E_1, E_2, E_3, E_4, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2
 ) where {P<:PEPSTensor}
     return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] :=
@@ -271,7 +297,7 @@
         C_2[χ4; χ5] *
         E_4[χ5 D7 D8; χ_out]
 end
-function halfinfinite_environment(
+function half_infinite_environment(
     C_1,
     C_2,
     E_1,
@@ -297,7 +323,7 @@
         E_4[χ5 D7 D8; χ6] *
         x[χ6 D11 D12]
 end
-function halfinfinite_environment(
+function half_infinite_environment(
     x::AbstractTensor{S,3},
     C_1,
     C_2,
@@ -310,7 +336,7 @@
     bra_1::P=ket_1,
     bra_2::P=ket_2,
 ) where {S,P<:PEPSTensor}
-    return @autoopt @tensor env_x[χ_in D_inabove D_inbelow] :=
+    return @autoopt @tensor x_env[χ_in D_inabove D_inbelow] :=
         x[χ1 D1 D2] *
         conj(E_1[χ1 D3 D4; χ2]) *
         conj(C_1[χ2; χ3]) *