Correct adjoint for truncated SVD

Project.toml

@@ -22,16 +22,16 @@ 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.10"
+MPSKit = "0.11"
 OptimKit = "0.3"
 Printf = "1"
 Statistics = "1"
-TensorKit = "0.12"
+TensorKit = "0.12.5"
 VectorInterface = "0.4"
 Zygote = "0.6"
-julia = "1.6"
+julia = "1.8"
 
 [extras]
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"

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; Venv=ComplexSpace(20))
-)
+ctm = leading_boundary(peps, CTMRG(; verbosity=1), CTMRGEnv(peps; Venv=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; Venv=ComplexSpace(20))
-)
+ctm2 = leading_boundary(peps2, CTMRG(; verbosity=1), CTMRGEnv(peps2; Venv=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,
@@ -27,6 +27,6 @@ alg = PEPSOptimize(;
 # E/N = −0.6694421, which is a QMC estimate from https://arxiv.org/abs/1101.3281.
 # Of course there is a noticable bias for small χbond and χenv.
 ψ₀ = InfinitePEPS(2, χbond)
-env₀ = leading_boundary(CTMRGEnv(ψ₀; Venv=ℂ^χenv), ψ₀, ctmalg)
-result = fixedpoint(ψ₀, H, alg, env₀)
+env₀ = leading_boundary(CTMRGEnv(ψ₀; Venv=ℂ^χenv), ψ₀, ctm_alg)
+result = fixedpoint(ψ₀, H, opt_alg, env₀)
 @show result.E

src/PEPSKit.jl

@@ -9,7 +9,7 @@ using TensorKit, KrylovKit, MPSKit, OptimKit
 using ChainRulesCore, Zygote
 
 include("utility/util.jl")
-include("utility/eigsolve.jl")
+include("utility/svd.jl")
 include("utility/rotations.jl")
 include("utility/hook_pullback.jl")
 
@@ -58,12 +58,14 @@ module Defaults
     const fpgrad_tol = 1e-6
 end
 
-export CTMRG, CTMRGEnv
+export FullSVD, IterSVD, OldSVD
+export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv
 export NLocalOperator, AnisotropicNNOperator, OnSite, NearestNeighbor
 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,T}
+    svd_alg::S = TensorKit.SVD()
+    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=TensorKit.SVD(),
+    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)) ||
@@ -146,8 +181,8 @@ function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T}
 
         # Decompose and multiply
         Up, _, Vp = tsvd!(ρ_prev)
-        Uf, _, Vf = tsvd!(ρ_final)
         Qprev = Up * Vp
+        Uf, _, Vf = tsvd!(ρ_final)
         Qfinal = Uf * Vf
         σ = Qprev * Qfinal'
 
@@ -283,7 +318,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(ϵ, ϵ₀)
@@ -298,7 +333,7 @@ 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
@@ -328,13 +363,13 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
             )
 
             # SVD half-infinite environment
-            trscheme = if alg.fixedspace == true
+            trscheme = if alg.trscheme isa FixedSpaceTruncation
                 truncspace(space(env.edges[WEST, row, cnext], 1))
             else
                 alg.trscheme
             end
             @tensor QQ[-1 -2 -3; -4 -5 -6] := Q_sw[-1 -2 -3; 1 2 3] * Q_nw[1 2 3; -4 -5 -6]
-            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