[WIP] Split CTMRG

docs/src/assets/pepskit.bib

@@ -11,6 +11,31 @@ @article{francuz_stable_2025
   doi = {10.1103/PhysRevResearch.7.013237}
 }
 
+@article{naumann_variationally_2025,
+  title = {Variationally optimizing infinite projected entangled-pair states at large bond dimensions: A split corner transfer matrix renormalization group approach},
+  author = {Naumann, Jan and Weerda, Erik L. and Eisert, Jens and Rizzi, Matteo and Schmoll, Philipp},
+  journal = {Phys. Rev. B},
+  volume = {111},
+  pages = {235116},
+  year = {2025},
+  month = {Jun},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevB.111.235116},
+  url = {https://doi.org/10.1103/PhysRevB.111.235116}
+}
+
+@article{fishman_faster_2018,
+  title = {Faster methods for contracting infinite two-dimensional tensor networks},
+  author = {Fishman, M. T. and Vanderstraeten, L. and Zauner-Stauber, V. and Haegeman, J. and Verstraete, F.},
+  journal = {Phys. Rev. B},
+  volume = {98},
+  pages = {235148},
+  year = {2018},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevB.98.235148},
+  url = {https://doi.org/10.1103/PhysRevB.98.235148}
+}
+
 @article{evenbly_gauge_2018,
   title = {Gauge fixing, canonical forms, and optimal truncations in tensor networks with closed loops},
   author = {Evenbly, Glen},

docs/src/man/environments.md

@@ -1,4 +1,46 @@
 # Environments
 
 !!! note
-    This section is still under construction.
+    This section is still under construction.
+
+## Split CTMRG
+
+[`SplitCTMRGEnv`](@ref) stores separate ket and bra factors for each CTMRG edge, connected by an interlayer bond.
+[`SplitCTMRG`](@ref) absorbs and projects the two layers sequentially using the full projectors of [Naumann et al., Phys. Rev. B 111, 235116 (2025)](https://doi.org/10.1103/PhysRevB.111.235116).
+[`SplitProjector`](@ref) uses the Fishman construction for all four projector types, factorizing each environment half by SVD and keeping its square-root-weighted factor facing the cut before computing the projectors.
+This improves the conditioning of the projector calculation at the cost of two additional SVDs per projector pair, following [Fishman et al.](@cite fishman_faster_2018), Appendix A.
+The outer and interlayer spaces can be chosen independently:
+
+```julia
+env₀ = SplitCTMRGEnv(ψ, ComplexSpace(32); interlayer_space = ComplexSpace(24))
+alg = SplitCTMRG(; tol = 1.0e-8, maxiter = 100)
+env, info = leading_boundary(env₀, ψ, alg)
+energy = expectation_value(ψ, H, CTMRGEnv(env))
+```
+
+By default, both spaces remain fixed during iteration.
+Use `trunc` for the outer projectors and `interlayer_trunc` for the interlayer projectors to select other truncation strategies.
+When omitted, `interlayer_trunc` inherits `trunc`; the default `FixedSpaceTruncation()` preserves the existing spaces of each bond independently.
+`SplitCTMRGEnv(env::CTMRGEnv)` factorizes an existing environment exactly, while its `trunc` keyword can reduce the interlayer spaces.
+`CTMRGEnv(env::SplitCTMRGEnv)` contracts the interlayer bonds so the existing expectation-value routines can be reused.
+This version supports PEPS overlap networks, Fishman projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref).
+The half-environment factors use direct square roots of their singular values, masking weights whose original singular values fall below the relative cutoff `eps(real(scalartype(ψ)))^(3 / 4)`; the projector pseudoinverse uses the same cutoff.
+
+GMRES is the default solver for `FixedPointGradient` when optimization uses `SplitCTMRG`; select another solver with `gradient_alg = (; solver_alg = (; alg = :Arnoldi))` or supply a gradient-algorithm instance.
+
+Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref PEPSKit.ScramblingEnvGauge) on the alternating ket and bra boundary factors to align both outer and interlayer gauges.
+The ordinary and split environments use the same boundary-MPS gauge method, with each split edge contributing two MPS sites instead of one.
+Both sets of bond spaces are held fixed in the differentiated iteration, including when the forward run uses other truncation strategies.
+Use a complex environment for fixed-point gradients, as for conventional asymmetric CTMRG.
+The first version contracts split edges into `CTMRGEnv` for energy evaluation and local preconditioning, preserving the derivative through both factors:
+
+```julia
+env₀ = SplitCTMRGEnv(randn, ComplexF64, ψ, ComplexSpace(32); interlayer_space = ComplexSpace(24))
+opt_alg = PEPSOptimize(;
+    boundary_alg = SplitCTMRG(; tol = 1.0e-8),
+    gradient_alg = FixedPointGradient(; solver_alg = (; alg = :GMRES)),
+)
+ψ_opt, env_opt, energy, info = fixedpoint(H, ψ, env₀, opt_alg)
+```
+
+`ImplicitGradient` is not supported for split CTMRG.

src/PEPSKit.jl

@@ -63,6 +63,7 @@ include("utility/twistdual.jl")
 include("utility/eigh.jl")
 include("utility/svd.jl")
 include("utility/qr.jl")
+include("utility/lq.jl")
 include("utility/rotations.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
@@ -86,6 +87,7 @@ include("operators/lattices/squarelattice.jl")
 include("operators/models.jl")
 
 include("environments/ctmrg_environments.jl")
+include("environments/split_ctmrg_environments.jl")
 include("environments/vumps_environments.jl")
 include("environments/suweight.jl")
 include("environments/bp_environments.jl")
@@ -95,11 +97,14 @@ include("algorithms/contractions/ctmrg/types.jl")
 include("algorithms/contractions/ctmrg/expr.jl")
 include("algorithms/contractions/ctmrg/enlarge_corner.jl")
 include("algorithms/contractions/ctmrg/enlarge_corner_column.jl")
+include("algorithms/contractions/ctmrg/enlarge_corner_split.jl")
 include("algorithms/contractions/ctmrg/projector.jl")
 include("algorithms/contractions/ctmrg/halfinf_env.jl")
 include("algorithms/contractions/ctmrg/fullinf_env.jl")
+include("algorithms/contractions/ctmrg/split_env.jl")
 include("algorithms/contractions/ctmrg/renormalize_corner.jl")
 include("algorithms/contractions/ctmrg/renormalize_edge.jl")
+include("algorithms/contractions/ctmrg/renormalize_split.jl")
 include("algorithms/contractions/ctmrg/network_value.jl")
 include("algorithms/contractions/ctmrg/gaugefix.jl")
 include("algorithms/contractions/ctmrg/characteristic_equations.jl")
@@ -128,11 +133,14 @@ include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/projectors/projectors.jl")
 include("algorithms/ctmrg/projectors/halfinfinite.jl")
 include("algorithms/ctmrg/projectors/fullinfinite.jl")
+include("algorithms/ctmrg/projectors/split.jl")
 include("algorithms/ctmrg/projectors/c4v_eigh.jl")
 include("algorithms/ctmrg/projectors/c4v_qr.jl")
 include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
+include("algorithms/ctmrg/split.jl")
 include("algorithms/ctmrg/gaugefix.jl")
+include("algorithms/ctmrg/gaugefix_split.jl")
 include("algorithms/ctmrg/c4v.jl")
 include("algorithms/ctmrg/initialization.jl")
 
@@ -172,13 +180,14 @@ include("algorithms/select_algorithm.jl")
 
 using .Defaults: set_scheduler!
 export set_scheduler!
-export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint
+export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint, LQAdjoint
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
+export SplitCTMRGEnv, SplitCTMRG
 export initialize_ctmrg_environment,
     RandomInitialization, ProductStateInitialization, ApplicationInitialization, IdentityInitialization
 export corner, edge, setcorner!, setedge!
 export FixedSpaceTruncation, SiteDependentTruncation
-export HalfInfiniteProjector, FullInfiniteProjector
+export HalfInfiniteProjector, FullInfiniteProjector, SplitProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
 export initialize_random_c4v_env, initialize_singlet_c4v_env
 export LocalOperator, physicalspace

src/algorithms/contractions/ctmrg/enlarge_corner_split.jl

@@ -0,0 +1,30 @@
+# Corner enlargements
+# -------------------
+
+"""Build the northwest enlarged corner for the environment of the teal projectors."""
+function enlarge_northwest_corner((Wk, Wb)::Tuple, C, N, A::PEPSSandwich, green)
+    return @autoopt @tensor Q[χout Dsb; χne Dek Deb] :=
+        green[χout; χsw Dsk] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C[χnw; χnnw] * N[χnnw Dnk Dnb; χne] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+end
+
+"""Build the northeast enlarged corner for the environment of the teal projectors."""
+function enlarge_northeast_corner(N, C, (Ek, Eb)::Tuple, A::PEPSSandwich, green)
+    return @autoopt @tensor Q[χnw Dwk Dwb; χout Dsb] :=
+        N[χnw Dnk Dnb; χnne] * C[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * green[χse Dsk; χout] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+end
+
+"""Build the southwest enlarged corner for the environment of the teal projectors."""
+function enlarge_southwest_corner(S, C, (Wk, Wb)::Tuple, A::PEPSSandwich, green)
+    return @autoopt @tensor Q[χse Dek Deb; χout Dnb] :=
+        S[χse Dsk Dsb; χssw] * C[χssw; χsw] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * green[χnw Dnk; χout] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+end
+
+"""Build the southeast enlarged corner for the environment of the teal projectors."""
+function enlarge_southeast_corner((Ek, Eb)::Tuple, C, S, A::PEPSSandwich, green)
+    return @autoopt @tensor Q[χout Dnb; χsw Dwk Dwb] :=
+        green[χout; χne Dnk] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * C[χse; χsse] * S[χsse Dsk Dsb; χsw] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+end

src/algorithms/contractions/ctmrg/gaugefix.jl

@@ -28,7 +28,7 @@ $(SIGNATURES)
 
 Apply `fix_gauge_corner` to the northwest corner with appropriate row and column indices.
 """
-function fix_gauge_northwest_corner((row, col), env::CTMRGEnv, signs)
+function fix_gauge_northwest_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs)
     return fix_gauge_corner(
         corner(env, NORTHWEST, row, col),
         signs[WEST, row, col],
@@ -41,7 +41,7 @@ $(SIGNATURES)
 
 Apply `fix_gauge_corner` to the northeast corner with appropriate row and column indices.
 """
-function fix_gauge_northeast_corner((row, col), env::CTMRGEnv, signs)
+function fix_gauge_northeast_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs)
     return fix_gauge_corner(
         corner(env, NORTHEAST, row, col),
         signs[NORTH, row, col],
@@ -54,7 +54,7 @@ $(SIGNATURES)
 
 Apply `fix_gauge_corner` to the southeast corner with appropriate row and column indices.
 """
-function fix_gauge_southeast_corner((row, col), env::CTMRGEnv, signs)
+function fix_gauge_southeast_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs)
     return fix_gauge_corner(
         corner(env, SOUTHEAST, row, col),
         signs[EAST, row, col],
@@ -67,7 +67,7 @@ $(SIGNATURES)
 
 Apply `fix_gauge_corner` to the southwest corner with appropriate row and column indices.
 """
-function fix_gauge_southwest_corner((row, col), env::CTMRGEnv, signs)
+function fix_gauge_southwest_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs)
     return fix_gauge_corner(
         corner(env, SOUTHWEST, row, col),
         signs[SOUTH, row, col],

src/algorithms/contractions/ctmrg/renormalize_split.jl

@@ -0,0 +1,30 @@
+# Renormalization
+# ---------------
+
+"""Apply the green and teal projectors sequentially when absorbing a north boundary edge."""
+function renormalize_northwest_corner(C::CTMRGCornerTensor, N::CTMRG_PEPS_EdgeTensor, green, teal)
+    @tensor tmp[χg Dnb; χe] := green[χg; χw Dnk] * C[χw; χn] * N[χn Dnk Dnb; χe]
+    return @tensor χout[χw; χe] := teal[χw; χg Dnb] * tmp[χg Dnb; χe]
+end
+
+"""Apply the green and teal projectors sequentially when absorbing a south boundary edge."""
+function renormalize_southwest_corner(C::CTMRGCornerTensor, S::CTMRG_PEPS_EdgeTensor, green, teal)
+    @tensor tmp[χe Dnb; χg] := S[χe Dnk Dnb; χs] * C[χs; χw] * green[χw Dnk; χg]
+    return @tensor χout[χe; χw] := tmp[χe Dnb; χg] * teal[χg Dnb; χw]
+end
+
+"""Absorb the ket layer and apply the split projectors to a west-edge factor."""
+function renormalize_west_edge(K::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor)
+    return @autoopt @tensor E[χinew Dek; χnnew] :=
+        K[χiold Dwk; χnold] * A[d; Dnk Dek Dsk Dwk] * green[χnold Dnk; χng] *
+        yellow[χiy; χiold Dsk] * teal[χng Dnb; χnnew] * red[χinew; χiy d Dnb]
+end
+
+"""Absorb the bra layer and apply the split projectors to the other west-edge factor."""
+function renormalize_west_edge(B::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor, ::Val{:bra})
+    @autoopt @tensor BY[χsg Dwb; χiy] :=
+        green[χsg; χsold Dsk] * B[χsold Dwb; χiold] * yellow[χiold Dsk; χiy]
+    @tensor BYT[χsnew Dwb Dsb; χiy] := teal[χsnew; χsg Dsb] * BY[χsg Dwb; χiy]
+    @tensor BYTA[χsnew Dnb Deb; χiy d] := BYT[χsnew Dwb Dsb; χiy] * conj(A[d; Dnb Deb Dsb Dwb])
+    return @tensor E[χsnew Deb; χinew] := BYTA[χsnew Dnb Deb; χiy d] * red[χiy d Dnb; χinew]
+end

src/algorithms/contractions/ctmrg/split_env.jl

@@ -0,0 +1,80 @@
+# Projector environments (green, teal, yellow, red)
+# -------------------------------------------------
+
+"""Collect the tensors surrounding a site in a split-CTMRG left move."""
+function _split_local_environment((r, c), network, env::SplitCTMRGEnv)
+    return (;
+        C1 = corner(env, NORTHWEST, r - 1, c - 1),
+        C2 = corner(env, NORTHEAST, r - 1, c + 1),
+        C3 = corner(env, SOUTHEAST, r + 1, c + 1),
+        C4 = corner(env, SOUTHWEST, r + 1, c - 1),
+        N = edge(env, NORTH, r - 1, c), S = edge(env, SOUTH, r + 1, c),
+        Wk = edge(env, :ket, WEST, r, c - 1), Wb = edge(env, :bra, WEST, r, c - 1),
+        Ek = edge(env, :ket, EAST, r, c + 1), Eb = edge(env, :bra, EAST, r, c + 1),
+        A = network[r, c],
+    )
+end
+
+"""Build the two halves of the environment for the green projectors (Fig. 3)."""
+function _split_green_environment((r, c), network, env::SplitCTMRGEnv)
+    (; C1, C2, C3, C4, N, S, Wk, Wb, Ek, Eb, A) = _split_local_environment((r, c), network, env)
+    top = @autoopt @tensor Q[χsw Dsk; χse Dsb] :=
+        Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χnne] *
+        C2[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+    bottom = @autoopt @tensor Q[χne Dnb; χnw Dnk] :=
+        Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C4[χssw; χsw] * S[χsse Dsk Dsb; χssw] *
+        C3[χse; χsse] * Ek[χne Dek; χie] * Eb[χie Deb; χse] *
+        ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+    return top, bottom
+end
+
+"""Build the two halves of the environment for the teal projectors (Fig. 7)."""
+function _split_teal_environment((r, c), network, env::SplitCTMRGEnv, green_west, green_east)
+    r = mod1(r, size(env, 2))
+    l = _split_local_environment((r, c), network, env)
+    h = _split_local_environment((r, c + 1), network, env)
+    nextrow = _next(r, size(env, 2))
+    NW = enlarge_northwest_corner((l.Wk, l.Wb), l.C1, l.N, l.A, green_west[2][nextrow])
+    NE = enlarge_northeast_corner(h.N, h.C2, (h.Ek, h.Eb), h.A, green_east[1][nextrow])
+    SW = enlarge_southwest_corner(l.S, l.C4, (l.Wk, l.Wb), l.A, green_west[1][r])
+    SE = enlarge_southeast_corner((h.Ek, h.Eb), h.C3, h.S, h.A, green_east[2][r])
+    return NW * NE, SE * SW
+end
+
+"""Build the two halves of the environment for the yellow projectors (Fig. 8)."""
+function _split_yellow_environment((r, c), network, env::SplitCTMRGEnv)
+    (; C1, C2, C3, C4, N, S, Wk, Wb, Ek, Eb, A) = _split_local_environment((r, c), network, env)
+    top = @autoopt @tensor Q[χiw Dsk; χie Dnb d] :=
+        Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χnne] * C2[χnne; χne] *
+        Ek[χne Dek; χie] * ket(A)[d; Dnk Dek Dsk Dwk]
+    bottom = @autoopt @tensor Q[χie Dnb d; χiw Dsk] :=
+        Wb[χsw Dwb; χiw] * C4[χssw; χsw] * S[χsse Dsk Dsb; χssw] * C3[χse; χsse] *
+        Eb[χie Deb; χse] * conj(bra(A)[d; Dnb Deb Dsb Dwb])
+    return top, bottom
+end
+
+"""Build the quarter-environments for the red projectors before multiplying their factors (Fig. 9)."""
+function _split_red_quarters((r, c), network, env::SplitCTMRGEnv, yellow_west, yellow_east)
+    l = _split_local_environment((r, c), network, env)
+    h = _split_local_environment((r, c + 1), network, env)
+    Ywl, Ywr = yellow_west[1][r], yellow_west[2][r]
+    Yel, Yer = yellow_east[1][r], yellow_east[2][r]
+    C1, N, Wk, A = l.C1, l.N, l.Wk, l.A
+    @autoopt @tensor NW[χout d Dnb; χne Dek] :=
+        Ywr[χout; χiw Dsk] * Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χne] *
+        ket(A)[d; Dnk Dek Dsk Dwk]
+    C2, N, Ek, A = h.C2, h.N, h.Ek, h.A
+    @autoopt @tensor NE[χnw Dwk; χout d Dnb] :=
+        N[χnw Dnk Dnb; χnne] * C2[χnne; χne] * Ek[χne Dek; χie] * Yel[χie Dsk; χout] *
+        ket(A)[d; Dnk Dek Dsk Dwk]
+    C4, S, Wb, A = l.C4, l.S, l.Wb, l.A
+    @autoopt @tensor SW[χse Deb; χout d Dnb] :=
+        S[χse Dsk Dsb; χssw] * C4[χssw; χsw] * Wb[χsw Dwb; χiw] * Ywl[χiw Dsk; χout] *
+        conj(bra(A)[d; Dnb Deb Dsb Dwb])
+    C3, S, Eb, A = h.C3, h.S, h.Eb, h.A
+    @autoopt @tensor SE[χout d Dnb; χsw Dwb] :=
+        Yer[χout; χie Dsk] * Eb[χie Deb; χse] * C3[χse; χsse] * S[χsse Dsk Dsb; χsw] *
+        conj(bra(A)[d; Dnb Deb Dsb Dwb])
+    return (; NW, NE, SE, SW)
+end