Add C4v CTMRG

Project.toml

@@ -1,7 +1,7 @@
 name = "PEPSKit"
 uuid = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
-authors = ["Paul Brehmer", "Lander Burgelman", "Lukas Devos <ldevos98@gmail.com>"]
 version = "0.7.0"
+authors = ["Paul Brehmer", "Lander Burgelman", "Lukas Devos <ldevos98@gmail.com>"]
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

src/Defaults.jl

@@ -101,9 +101,20 @@ const svd_rrule_alg = :full # ∈ {:full, :gmres, :bicgstab, :arnoldi}
 const svd_rrule_broadening = 1.0e-13
 const krylovdim_factor = 1.4
 
+# eigh forward & reverse
+const eigh_fwd_alg = :qriteration
+const eigh_rrule_alg = :trunc
+const eigh_rrule_verbosity = 0
+
+# QR forward & reverse
+# const qr_fwd_alg = :something # TODO
+# const qr_rrule_alg = :something
+# const qr_rrule_verbosity = :something
+
 # Projectors
 const projector_alg = :halfinfinite # ∈ {:halfinfinite, :fullinfinite}
 const projector_verbosity = 0
+const projector_alg_c4v = :c4v_eigh
 
 # Fixed-point gradient
 const gradient_tol = 1.0e-6

src/PEPSKit.jl

@@ -29,7 +29,9 @@ include("Defaults.jl")  # Include first to allow for docstring interpolation wit
 
 include("utility/util.jl")
 include("utility/diffable_threads.jl")
+include("utility/eigh.jl")
 include("utility/svd.jl")
+# include("utility/qr.jl")
 include("utility/rotations.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
@@ -42,6 +44,8 @@ include("networks/infinitesquarenetwork.jl")
 include("states/infinitepeps.jl")
 include("states/infinitepartitionfunction.jl")
 
+include("utility/symmetrization.jl")
+
 include("operators/infinitepepo.jl")
 include("operators/transfermatrix.jl")
 include("operators/localoperator.jl")
@@ -71,6 +75,7 @@ include("algorithms/ctmrg/projectors.jl")
 include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
+include("algorithms/ctmrg/c4v.jl")
 
 include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
@@ -89,8 +94,6 @@ include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
 include("algorithms/correlators.jl")
 
-include("utility/symmetrization.jl")
-
 include("algorithms/optimization/fixed_point_differentiation.jl")
 include("algorithms/optimization/peps_optimization.jl")
 
@@ -102,6 +105,7 @@ export SVDAdjoint, FullSVDReverseRule, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
+export EighAdjoint, IterEigh, C4vCTMRG, C4vEighProjector, C4vQRProjector
 export LocalOperator, physicalspace
 export product_peps
 export reduced_densitymatrix, expectation_value, network_value, cost_function

src/algorithms/ctmrg/c4v.jl

@@ -0,0 +1,144 @@
+struct C4vCTMRG{P <: ProjectorAlgorithm} <: CTMRGAlgorithm
+    tol::Float64
+    maxiter::Int
+    miniter::Int
+    verbosity::Int
+    projector_alg::P
+end
+function C4vCTMRG(; kwargs...)
+    return CTMRGAlgorithm(; alg = :c4v, kwargs...)
+end
+CTMRG_SYMBOLS[:c4v] = C4vCTMRG
+
+struct C4vEighProjector{S <: EighAdjoint, T} <: ProjectorAlgorithm
+    alg::S
+    trunc::T
+    verbosity::Int
+end
+function C4vEighProjector(; kwargs...)
+    return ProjectorAlgorithm(; alg = :c4v_eigh, kwargs...)
+end
+PROJECTOR_SYMBOLS[:c4v_eigh] = C4vEighProjector
+
+# struct C4vQRProjector{S, T} <: ProjectorAlgorithm
+#     alg::S
+#     verbosity::Int
+# end
+# function C4vQRProjector(; kwargs...)
+#     return ProjectorAlgorithm(; alg = :c4v_qr, kwargs...)
+# end
+# PROJECTOR_SYMBOLS[:c4v_qr] = C4vQRProjector
+
+#
+## C4v-symmetric CTMRG iteration (called through `leading_boundary`)
+#
+
+function ctmrg_iteration(
+        network,
+        env::CTMRGEnv,
+        alg::C4vCTMRG,
+    )
+    enlarged_corner = c4v_enlarge(network, env, alg.projector_alg)
+    corner′, projector, info = c4v_projector(enlarged_corner, alg.projector_alg)
+    edge′ = c4v_renormalize(network, env, projector)
+    return CTMRGEnv(corner′, edge′), info
+end
+
+function c4v_enlarge(network, env, ::C4vEighProjector)
+    return TensorMap(EnlargedCorner(network, env, (NORTHWEST, 1, 1)))
+end
+# function c4v_enlarge(enlarged_corner, alg::C4vQRProjector)
+#     # TODO
+# end
+
+function c4v_projector(enlarged_corner, alg::C4vEighProjector)
+    hermitian_corner = 0.5 * (enlarged_corner + enlarged_corner') / norm(enlarged_corner)
+    trunc = truncation_strategy(alg, enlarged_corner)
+    D, U, info = eigh_trunc!(hermitian_corner, decomposition_algorithm(alg); trunc)
+
+    # Check for degenerate eigenvalues
+    Zygote.isderiving() && ignore_derivatives() do
+        if alg.verbosity > 0 && is_degenerate_spectrum(D)
+            vals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(D))
+            @warn("degenerate eigenvalues detected: ", vals)
+        end
+    end
+
+    return D / norm(D), U, (; D, U, info...)
+end
+# function c4v_projector(enlarged_corner, alg::C4vQRProjector)
+#     # TODO
+# end
+
+function c4v_renormalize(network, env, projector)
+    new_edge = renormalize_north_edge(env.edges[1], projector, projector', network[1, 1])
+    return new_edge / norm(new_edge)
+end
+
+function CTMRGEnv(
+        corner::AbstractTensorMap{T, S, 1, 1}, edge::AbstractTensorMap{T′, S, N, 1}
+    ) where {T, T′, S, N}
+    corners = fill(corner, 4, 1, 1)
+    edge_SW = physical_flip(edge)
+    edges = reshape([edge, edge, edge_SW, edge_SW], (4, 1, 1))
+    return CTMRGEnv(corners, edges)
+end
+
+#
+## utility
+#
+
+# Adjoint of an edge tensor, but permutes the physical spaces back into the codomain.
+# Intuitively, this conjugates a tensor and then reinterprets its 'direction' as an edge tensor.
+function _dag(A::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
+    return permute(A', ((1, (3:(N + 1))...), (2,)))
+end
+function physical_flip(A::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
+    return flip(A, 2:N)
+end
+function project_hermitian(E::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
+    E´ = (E + physical_flip(_dag(E))) / 2
+    return E´
+end
+function project_hermitian(C::AbstractTensorMap{T, S, 1, 1}) where {T, S}
+    C´ = (C + C') / 2
+    return C´
+end
+
+# should perform this check at the beginning of `leading_boundary` really...
+function check_symmetry(state, ::RotateReflect; atol = 1.0e-10)
+    @assert length(state) == 1 "check_symmetry only works for single site unit cells"
+    @assert norm(state[1] - _fit_spaces(rotl90(state[1]), state[1])) /
+        norm(state[1]) < atol "not rotation invariant"
+    @assert norm(state[1] - _fit_spaces(herm_depth(state[1]), state[1])) /
+        norm(state[1]) < atol "not hermitian-reflection invariant"
+    return nothing
+end
+
+#
+## environment initialization
+#
+
+# TODO: rewrite this using `initialize_environment` and C4v-specific initialization algorithms
+# environment with dummy corner singlet(V) ← singlet(V) and identity edge V ← V, initialized at dim(Venv)
+# function initialize_singlet_c4v_env(Vpeps::ElementarySpace, Venv::ElementarySpace, T = ComplexF64)
+#     corner₀ = DiagonalTensorMap(zeros(real(T), Venv ← Venv))
+#     corner₀.data[1] = one(real(T))
+#     edge₀ = permute(id(T, Venv ⊗ Vpeps), ((1, 2, 4), (3,)))
+#     return CTMRGEnv(corner₀, edge₀)
+# end
+
+function initialize_random_c4v_env(Vstate::VectorSpace, Venv::ElementarySpace, T = ComplexF64)
+    corner₀ = DiagonalTensorMap(randn(real(T), Venv ← Venv))
+    edge₀ = randn(T, Venv ⊗ Vstate ← Venv)
+    edge₀ = project_hermitian(edge₀)
+    return CTMRGEnv(corner₀, edge₀)
+end
+function initialize_random_c4v_env(state::InfinitePEPS, Venv::ElementarySpace, T = scalartype(state))
+    Vpeps = domain(state[1])[1]
+    return initialize_random_c4v_env(Vpeps ⊗ Vpeps', Venv, T)
+end
+function initialize_random_c4v_env(state::InfinitePartitionFunction, Venv::ElementarySpace, T = scalartype(state))
+    Vpf = domain(state[1])[1]
+    return initialize_random_c4v_env(Vpf, Venv, T)
+end

src/algorithms/ctmrg/ctmrg.jl

@@ -19,17 +19,19 @@ function CTMRGAlgorithm(;
         maxiter = Defaults.ctmrg_maxiter, miniter = Defaults.ctmrg_miniter,
         verbosity = Defaults.ctmrg_verbosity,
         trunc = (; alg = Defaults.trunc),
-        svd_alg = (;),
+        decomposition_alg = (;),
         projector_alg = Defaults.projector_alg, # only allows for Symbol/NamedTuple to expose projector kwargs
     )
     # replace symbol with projector alg type
     haskey(CTMRG_SYMBOLS, alg) || throw(ArgumentError("unknown CTMRG algorithm: $alg"))
     alg_type = CTMRG_SYMBOLS[alg]
 
     # parse CTMRG projector algorithm
-
+    if alg == :c4v && projector_alg == Defaults.projector_alg
+        projector_alg = Defaults.projector_alg_c4v
+    end
     projector_algorithm = ProjectorAlgorithm(;
-        alg = projector_alg, svd_alg, trunc, verbosity
+        alg = projector_alg, decomposition_alg, trunc, verbosity
     )
 
     return alg_type(tol, maxiter, miniter, verbosity, projector_algorithm)
@@ -76,7 +78,7 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st
     - `:truncrank` : Additionally supply truncation dimension `η`; truncate such that the 2-norm of the truncated values is smaller than `η`
     - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space 
     - `:trunctol` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η`
-* `svd_alg::Union{<:SVDAdjoint,NamedTuple}` : SVD algorithm for computing projectors. See also [`SVDAdjoint`](@ref). By default, a reverse-rule tolerance of `tol=1e1tol` where the `krylovdim` is adapted to the `env₀` environment dimension.
+* `decomposition_alg::Union{<:DecompositionAdjoint,NamedTuple}` : Tensor decomposition algorithm for computing projectors. See e.g. [`SVDAdjoint`](@ref). 
 * `projector_alg::Symbol=:$(Defaults.projector_alg)` : Variant of the projector algorithm. See also [`ProjectorAlgorithm`](@ref).
     - `:halfinfinite` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments.
     - `:fullinfinite` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments.

src/algorithms/ctmrg/gaugefix.jl

@@ -1,3 +1,39 @@
+"""
+    gauge_fix(alg::CTMRGAlgorithm, signs, info)
+    gauge_fix(alg::ProjectorAlgorithm, signs, info)
+
+Fix the free gauges of the tensor decompositions associated with `alg`.
+"""
+function gauge_fix(alg::CTMRGAlgorithm, signs, info)
+    alg_fixed = @set alg.projector_alg = gauge_fix(alg.projector_alg, signs, info)
+    return alg_fixed
+end
+function gauge_fix(alg::ProjectorAlgorithm, signs, info)
+    decomposition_alg_fixed = gauge_fix(alg.alg, signs, info) # every ProjectorAlgorithm needs an `alg` field
+    alg_fixed = @set alg.alg = decomposition_alg_fixed
+    alg_fixed = @set alg_fixed.trunc = notrunc()
+    return alg_fixed
+end
+
+"""
+$(TYPEDEF)
+
+CTMRG environment gauge fixing algorithm implementing the "general" technique from
+https://arxiv.org/abs/2311.11894. This works by constructing a transfer matrix consisting
+of an edge tensor and a random MPS, thus scrambling potential degeneracies, and then
+performing a QR decomposition to extract the gauge signs. This is adapted accordingly for
+asymmetric CTMRG algorithms using multi-site unit cell transfer matrices.
+"""
+struct ScramblingEnvGauge end
+
+"""
+$(TYPEDEF)
+
+C4v-symmetric equivalent of the [ScramblingEnvGauge`](@ref) environment gauge fixing
+algorithm.
+"""
+struct ScramblingEnvGaugeC4v end
+
 """
 $(SIGNATURES)
 
@@ -6,15 +42,14 @@ This assumes that the `envfinal` is the result of one CTMRG iteration on `envpre
 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}
+function gauge_fix(envfinal::CTMRGEnv{C, T}, ::ScramblingEnvGauge, envprev::CTMRGEnv{C, T}) where {C, T}
     # Check if spaces in envprev and envfinal are the same
     same_spaces = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c)
         space(envfinal.edges[dir, r, c]) == space(envprev.edges[dir, r, c]) &&
             space(envfinal.corners[dir, r, c]) == space(envprev.corners[dir, r, c])
     end
     @assert all(same_spaces) "Spaces of envprev and envfinal are not the same"
 
-    # Try the "general" algorithm from https://arxiv.org/abs/2311.11894
     signs = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c)
         # Gather edge tensors and pretend they're InfiniteMPSs
         if dir == NORTH
@@ -54,6 +89,42 @@ function gauge_fix(envprev::CTMRGEnv{C, T}, envfinal::CTMRGEnv{C, T}) where {C,
     return fix_global_phases(envprev, CTMRGEnv(cornersfix, edgesfix)), signs
 end
 
+# C4v specialized gauge fixing routine with Hermitian transfer matrix
+function gauge_fix(envfinal::CTMRGEnv{C, T}, ::ScramblingEnvGaugeC4v, envprev::CTMRGEnv{C, T}) where {C, T}
+    # Check if spaces in envprev and envfinal are the same
+    same_spaces = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c)
+        space(envfinal.edges[dir, r, c]) == space(envprev.edges[dir, r, c]) &&
+            space(envfinal.corners[dir, r, c]) == space(envprev.corners[dir, r, c])
+    end
+    @assert all(same_spaces) "Spaces of envprev and envfinal are not the same"
+
+    # "general" algorithm from https://arxiv.org/abs/2311.11894
+    Tprev = envprev.edges[1, 1, 1]
+    Tfinal = envfinal.edges[1, 1, 1]
+
+    # Random Hermitian MPS of same bond dimension
+    # (make Hermitian such that T-M transfer matrix has real eigenvalues)
+    M = project_hermitian(randn(scalartype(Tfinal), space(Tfinal)))
+
+    # Find right fixed points of mixed transfer matrices
+    ρinit = randn(
+        scalartype(T), MPSKit._lastspace(Tfinal)' ← MPSKit._lastspace(M)'
+    )
+    ρprev = c4v_transfermatrix_fixedpoint(Tprev, M, ρinit)
+    ρfinal = c4v_transfermatrix_fixedpoint(Tfinal, M, ρinit)
+
+    # Decompose and multiply
+    Qprev, = left_orth!(ρprev)
+    Qfinal, = left_orth!(ρfinal)
+
+    σ = Qprev * Qfinal'
+
+    @tensor cornerfix[χ_in; χ_out] := σ[χ_in; χ1] * envfinal.corners[1][χ1; χ2] * conj(σ[χ_out; χ2])
+    @tensor edgefix[χ_in D_in_above D_in_below; χ_out] :=
+        σ[χ_in; χ1] * envfinal.edges[1][χ1 D_in_above D_in_below; χ2] * conj(σ[χ_out; χ2])
+    return CTMRGEnv(cornerfix, edgefix), fill(σ, (4, 1, 1))
+end
+
 # this is a bit of a hack to get the fixed point of the mixed transfer matrix
 # because MPSKit is not compatible with AD
 # NOTE: the action of the transfer operator here is NOT the same as that of
@@ -82,6 +153,13 @@ function transfermatrix_fixedpoint(tops, bottoms, ρinit)
     end
     return first(vecs)
 end
+function c4v_transfermatrix_fixedpoint(top, bottom, ρinit)
+    _, vecs, info = eigsolve(ρinit, 1, :LM, Lanczos()) do ρ
+        PEPSKit.mps_transfer_right(ρ, top, bottom)
+    end
+    info.converged > 0 || @warn "eigsolve did not converge"
+    return first(vecs)
+end
 
 # Explicit fixing of relative phases (doing this compactly in a loop is annoying)
 function fix_relative_phases(envfinal::CTMRGEnv, signs)
@@ -168,7 +246,9 @@ end
 Check if the element-wise difference of the corner and edge tensors of the final and fixed
 CTMRG environments are below `atol` and return the maximal difference.
 """
-function calc_elementwise_convergence(envfinal::CTMRGEnv, envfix::CTMRGEnv; atol::Real = 1.0e-6)
+function calc_elementwise_convergence(
+        envfinal::CTMRGEnv, envfix::CTMRGEnv; atol::Real = 1.0e-6
+    )
     ΔC = envfinal.corners .- envfix.corners
     ΔCmax = norm(ΔC, Inf)
     ΔCmean = norm(ΔC)
@@ -194,5 +274,18 @@ function calc_elementwise_convergence(envfinal::CTMRGEnv, envfix::CTMRGEnv; atol
 
     return max(ΔCmax, ΔTmax)
 end
+function calc_elementwise_convergence(
+        envfinal::CTMRGEnv{C}, envfix::CTMRGEnv{C′}; kwargs...
+    ) where {C <: DiagonalTensorMap, C′} # case where one of the environments might have diagonal corners
+    return calc_elementwise_convergence( # make corners non-diagonal TensorMaps such that you can compute difference
+        CTMRGEnv(convert.(TensorMap, envfinal.corners), envfinal.edges; kwargs...),
+        CTMRGEnv(envfix.corners, envfix.edges; kwargs...)
+    )
+end
+function calc_elementwise_convergence(
+        envfinal::CTMRGEnv{C}, envfix::CTMRGEnv{C′}; kwargs...
+    ) where {C, C′ <: DiagonalTensorMap}
+    return calc_elementwise_convergence(envfix, envfinal)
+end
 
 @non_differentiable calc_elementwise_convergence(args...)