Document CTMRG, optimization and util things

src/PEPSKit.jl

@@ -31,7 +31,25 @@ include("algorithms/peps_opt.jl")
 
 include("utility/symmetrization.jl")
 
-# Default settings
+"""
+    module Defaults
+        const ctmrg_maxiter = 100
+        const ctmrg_miniter = 4
+        const ctmrg_tol = 1e-12
+        const fpgrad_maxiter = 100
+        const fpgrad_tol = 1e-6
+    end
+
+Module containing default values that represent typical algorithm parameters.
+
+- `ctmrg_maxiter = 100`: Maximal number of CTMRG iterations per run
+- `ctmrg_miniter = 4`: Minimal number of CTMRG carried out
+- `ctmrg_tol = 1e-12`: Tolerance checking singular value and norm convergence
+- `fpgrad_maxiter = 100`: Maximal number of iterations for computing the CTMRG
+fixed-point gradient
+- `fpgrad_tol = 1e-6`: Convergence tolerance for the fixed-point gradient iteration
+
+"""
 module Defaults
     const ctmrg_maxiter = 100
     const ctmrg_miniter = 4
@@ -49,7 +67,6 @@ export fixedpoint
 export InfinitePEPS, InfiniteTransferPEPS
 export InfinitePEPO, InfiniteTransferPEPO
 export initializeMPS, initializePEPS
-export PEPOOptimize, pepo_opt_environments
 export symmetrize, None, Depth, Full
 
 end # module

src/algorithms/ctmrg.jl

@@ -1,4 +1,18 @@
 # 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)
+
+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`.
+
+"""
 @kwdef struct CTMRG
     trscheme::TruncationScheme = TensorKit.notrunc()
     tol::Float64 = Defaults.ctmrg_tol
@@ -8,8 +22,12 @@
     fixedspace::Bool = false
 end
 
-# Compute CTMRG environment for a given state
-# function MPSKit.leading_boundary(state, alg::CTMRG, envinit=CTMRGEnv(state))
+"""
+    MPSKit.leading_boundary(state, alg::CTMRG, envinit=CTMRGEnv(state))
+
+Contract `state` using CTMRG and return the CTM environment.
+Per default, a random initial environment is used.
+"""
 function MPSKit.leading_boundary(state, alg::CTMRG, envinit=CTMRGEnv(state))
     normold = 1.0
     CSold = map(x -> tsvd(x; alg=TensorKit.SVD())[2], envinit.corners)
@@ -68,7 +86,7 @@ function MPSKit.leading_boundary(state, alg::CTMRG, envinit=CTMRGEnv(state))
         ϵold = ϵ
     end
 
-    # do one final iteration that does not change the spaces
+    # 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
     )
@@ -79,7 +97,14 @@ function MPSKit.leading_boundary(state, alg::CTMRG, envinit=CTMRGEnv(state))
     return envfix
 end
 
-# Fix gauge of corner end edge tensors from last and second last CTMRG iteration
+"""
+    gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T}
+
+Fix the gauge of `envfinal` based on the previous environment `envprev`.
+This assumes that the `envfinal` is the result of one CTMRG iteration on `envprev`.
+Given that the CTMRG run is converged, the returned environment will be
+element-wise converged to `envprev`.
+"""
 function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T}
     # Check if spaces in envprev and envfinal are the same
     same_spaces = map(Iterators.product(axes(envfinal.edges)...)) do (dir, r, c)
@@ -238,7 +263,14 @@ end
 
 @non_differentiable check_elementwise_convergence(args...)
 
-# One CTMRG iteration x′ = f(A, x)
+"""
+    ctmrg_iter(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
+    
+Perform one iteration of CTMRG that maps the `state` and `env` to a new environment,
+and also return the truncation error.
+One CTMRG iteration consists of four `left_move` calls and 90 degree rotations,
+such that the environment is grown and renormalized in all four directions.
+"""
 function ctmrg_iter(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
     ϵ = 0.0
 
@@ -252,7 +284,12 @@ function ctmrg_iter(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
     return env, ϵ
 end
 
-# Grow environment, compute projectors and renormalize
+"""
+    left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
+
+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}
     corners::typeof(env.corners) = copy(env.corners)
     edges::typeof(env.edges) = copy(env.edges)
@@ -343,6 +380,7 @@ function northwest_corner(E4, C1, E1, peps_above, peps_below=peps_above)
         conj(peps_below[7; 6 -6 -3 2])
 end
 
+# Compute enlarged NE corner
 function northeast_corner(E1, C2, E2, peps_above, peps_below=peps_above)
     @tensor corner[-1 -2 -3; -4 -5 -6] :=
         E1[-1 1 2; 3] *
@@ -352,6 +390,7 @@ function northeast_corner(E1, C2, E2, peps_above, peps_below=peps_above)
         conj(peps_below[7; 2 6 -6 -3])
 end
 
+# Compute enlarged SE corner
 function southeast_corner(E2, C3, E3, peps_above, peps_below=peps_above)
     @tensor corner[-1 -2 -3; -4 -5 -6] :=
         E2[-1 1 2; 3] *
@@ -384,7 +423,11 @@ function grow_env_left(peps, Pl, Pr, C_sw, C_nw, T_s, T_w, T_n)
     return C_sw′, C_nw′, T_w′
 end
 
-# Compute norm of the entire CMTRG enviroment
+"""
+    LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv)
+
+Compute the norm of a PEPS contracted with a CTM environment.
+"""
 function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv)
     total = one(scalartype(peps))
 

src/algorithms/peps_opt.jl

@@ -1,14 +1,14 @@
 abstract type GradMode end
 
 """
-    NaiveAD <: GradMode
+    struct NaiveAD <: GradMode
 
 Gradient mode for CTMRG using AD.
 """
 struct NaiveAD <: GradMode end
 
 """
-    GeomSum <: GradMode
+    struct GeomSum <: GradMode
 
 Gradient mode for CTMRG using explicit evaluation of the geometric sum.
 """
@@ -19,7 +19,7 @@ Gradient mode for CTMRG using explicit evaluation of the geometric sum.
 end
 
 """
-    ManualIter <: GradMode
+    struct ManualIter <: GradMode
 
 Gradient mode for CTMRG using manual iteration to solve the linear problem.
 """
@@ -29,7 +29,19 @@ Gradient mode for CTMRG using manual iteration to solve the linear problem.
     verbosity::Int = 0
 end
 
-# Algorithm struct containing parameters for PEPS optimization
+"""
+    struct PEPSOptimize{G}(; boundary_alg = CTMRG(), optimizer::OptimKit.OptimizationAlgorithm = LBFGS()
+                           reuse_env::Bool = true, gradient_alg::G, verbosity::Int = 0)
+
+Algorithm struct that represent PEPS ground-state optimization using AD.
+Set the algorithm to contract the infinite PEPS in `boundary_alg`;
+currently only `CTMRG` is supported. The `optimizer` computes the gradient directions
+based on the CTMRG gradient and updates the PEPS parameters. In this optimization,
+the CTMRG runs can be started on the converged environments of the previous optimizer
+step by setting `reuse_env` to true. Otherwise a random environment is used at each
+step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm.
+Different levels of output verbosity can be activated using `verbosity` (0, 1 or 2).
+"""
 @kwdef struct PEPSOptimize{G}
     boundary_alg::CTMRG = CTMRG()  # Algorithm to find boundary environment
     optimizer::OptimKit.OptimizationAlgorithm = LBFGS(
@@ -40,7 +52,13 @@ end
     verbosity::Int = 0
 end
 
-# Find ground-state PEPS and energy
+"""
+    fixedpoint(ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize,
+               env₀::CTMRGEnv=CTMRGEnv(ψ₀; Venv=field(T)^20)) where {T}
+    
+Optimize `ψ₀` with respect to the Hamiltonian `H` according to the parameters supplied
+in `alg`. The initial environment `env₀` serves as an initial guess for the first CTMRG run.
+"""
 function fixedpoint(
     ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, env₀::CTMRGEnv=CTMRGEnv(ψ₀; Venv=field(T)^20)
 ) where {T}
@@ -74,6 +92,14 @@ Evaluating the gradient of the cost function for CTMRG:
 =#
 using Zygote: @showgrad
 
+@doc """
+    ctmrg_gradient((peps, envs), H, alg::PEPSOptimize)
+
+Compute the energy and its gradient for the Hamiltonian `H``.
+Here, `alg` determines the algorithm that is used for the CTMRG gradient computation.
+"""
+ctmrg_gradient
+
 function ctmrg_gradient((peps, envs), H, alg::PEPSOptimize{NaiveAD})
     E, g = withgradient(peps) do ψ
         envs′ = leading_boundary(ψ, alg.boundary_alg, envs)

src/operators/localoperator.jl

@@ -1,13 +1,38 @@
 abstract type AbstractInteraction end
 
+"""
+    struct OnSite <: AbstractInteraction
+
+Trivial interaction representing terms that act on one isolated site.
+"""
 struct OnSite <: AbstractInteraction end
+
+"""
+    struct NearestNeighbor <: AbstractInteraction
+
+Interaction representing nearest neighbor terms that act on two adjacent sites.
+"""
 struct NearestNeighbor <: AbstractInteraction end
 
+"""
+    struct NLocalOperator{I<:AbstractInteraction}
+    
+Operator in form of a `AbstractTensorMap` which is parametrized by an interaction type.
+Mostly, this is used to define Hamiltonian terms and observables.
+"""
 struct NLocalOperator{I<:AbstractInteraction}
     op::AbstractTensorMap
 end
 
-# Contraction of CTMRGEnv and PEPS tensors with open physical bonds
+@doc """
+    operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::AbstractInteraction)
+
+Contract a PEPS and a CTMRG environment to form an operator environment.
+The open bonds correspond to the indices of an operator with the specified
+`AbstractInteraction` type.
+"""
+operator_env
+
 function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::OnSite)
     return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c)
         @tensor ρ[-1; -2] :=
@@ -24,7 +49,6 @@ function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::OnSite)
     end
 end
 
-# Horizontally extended contraction of CTMRGEnv and PEPS tensors with open physical bonds
 function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::NearestNeighbor)
     return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c)
         cnext = _next(c, size(peps, 2))
@@ -46,6 +70,14 @@ function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::NearestNeighbor)
     end
 end
 
+@doc """
+    MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator)
+
+Evaluate the expectation value of any `NLocalOperator` on each unit-cell entry
+of `peps` and `env`.
+"""
+MPSKit.expectation_value
+
 # 1-site operator expectation values on unit cell
 function MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator{OnSite})
     result = similar(peps.A, eltype(O.op))
@@ -60,7 +92,6 @@ function MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator{OnS
     return result
 end
 
-# 2-site operator expectation values on unit cell
 function MPSKit.expectation_value(
     peps::InfinitePEPS, env, O::NLocalOperator{NearestNeighbor}
 )
@@ -76,7 +107,12 @@ function MPSKit.expectation_value(
     return result
 end
 
-# ⟨O⟩ from vertical and horizontal nearest-neighbor contributions
+"""
+    costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor})
+    
+Compute the expectation value of a nearest-neighbor operator.
+This is used to evaluate and differentiate the energy in ground-state PEPS optimizations.
+"""
 function costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor})
     oh = sum(expectation_value(peps, env, op))
     ov = sum(expectation_value(rotl90(peps), rotl90(env), op))

src/states/abstractpeps.jl

@@ -6,6 +6,16 @@ conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W.
 """
 const PEPSTensor{S} = AbstractTensorMap{S,1,4} where {S<:ElementarySpace}
 
+"""
+    PEPSTensor(f, ::Type{T}, Pspace::S, Nspace::S, Espace::S=Nspace,
+               Sspace::S=Nspace', Wspace::S=Espace') where {T,S<:ElementarySpace}
+    PEPSTensor(f, ::Type{T}, Pspace::Int, Nspace::Int, Espace::Int=Nspace,
+               Sspace::Int=Nspace, Wspace::Int=Espace) where {T}
+                
+Construct a PEPS tensor based on the physical, north, east, west and south spaces.
+Alternatively, only the space dimensions can be provided and ℂ is assumed as the field.
+The tensor elements are generated based on `f` and the element type is specified in `T`.
+"""
 function PEPSTensor(
     f,
     ::Type{T},
@@ -17,7 +27,6 @@ function PEPSTensor(
 ) where {T,S<:ElementarySpace}
     return TensorMap(f, T, Pspace ← Nspace ⊗ Espace ⊗ Sspace ⊗ Wspace)
 end
-
 function PEPSTensor(
     f,
     ::Type{T},
@@ -41,14 +50,14 @@ const PEPOTensor{S} = AbstractTensorMap{S,2,4} where {S<:ElementarySpace}
 """
     abstract type AbstractPEPS end
 
-Abstract supertype for a 2D projected entangled pairs state.
+Abstract supertype for a 2D projected entangled-pair state.
 """
 abstract type AbstractPEPS end
 
 """
     abstract type AbstractPEPO end
 
-Abstract supertype for a 2D projected entangled pairs operator.
+Abstract supertype for a 2D projected entangled-pair operator.
 """
 abstract type AbstractPEPO end
 

src/utility/rotations.jl

@@ -8,5 +8,9 @@ const NORTHEAST = 2
 const SOUTHEAST = 3
 const SOUTHWEST = 4
 
-# Rotate tensor to any direction by successive application of Base.rotl90
+"""
+    rotate_north(t, dir)
+
+Rotate north direction of `t` to `dir` by successive applications of `rotl90`.
+"""
 rotate_north(t, dir) = mod1(dir, 4) == NORTH ? t : rotate_north(rotl90(t), dir - 1)

src/utility/util.jl

@@ -49,7 +49,13 @@ function is_degenerate_spectrum(
     return false
 end
 
-# Create empty projectors for given state without recomputing transpose
+"""
+    projector_type(T::DataType, size)
+
+Create two arrays of specified `size` that contain undefined tensors representing
+left and right acting projectors, respectively. The projector types are inferred
+from the TensorMap type `T` which avoids having to recompute tranpose tensors.
+"""
 function projector_type(T::DataType, size)
     Pleft = Array{T,length(size)}(undef, size)
     Prtype = tensormaptype(spacetype(T), numin(T), numout(T), storagetype(T))