[WIP] Add support for symmetric boundary MPS contractions

src/Defaults.jl

@@ -17,6 +17,18 @@ Module containing default algorithm parameter values and arguments.
 * `ctmrg_tol_max=$(Defaults.ctmrg_tol_max)` : Maximal CTMRG tolerance used by `ctmrg_dynamic_tols`.
 * `ctmrg_tol_factor=$(Defaults.ctmrg_tol_factor)` : Tolerance scaling factor used by `ctmrg_dynamic_tols`.
 
+## Boundary contraction
+
+* `boundary_alg=:$(Defaults.boundary_alg)` : Default algorithm family used to contract a network, see [`PEPSKit.BoundaryAlgorithm`](@ref).
+    - `:SimultaneousCTMRG`, `:SequentialCTMRG`, `:C4vCTMRG` : CTMRG variants.
+    - `:SymmetricBoundaryMPS` : Boundary MPS contraction of a fully symmetric single-site network.
+
+## Boundary MPS
+
+* `boundarymps_mps_alg=:$(Defaults.boundarymps_mps_alg)` : Default MPS optimization algorithm driving a boundary MPS contraction.
+    - `:VUMPS` : Variational uniform MPS.
+    - `:VOMPS` : Variational optimization of the MPS through MPO-MPS overlap maximization.
+
 ## SVD forward & reverse
 
 * `trunc=:$(Defaults.trunc)` : Truncation scheme for SVDs and other decompositions.
@@ -150,6 +162,12 @@ const ctmrg_tol_min = 1.0e-12
 const ctmrg_tol_max = 1.0e-4
 const ctmrg_tol_factor = 1.0e-3
 
+# Boundary contraction
+const boundary_alg = :SimultaneousCTMRG # ∈ {:SimultaneousCTMRG, :SequentialCTMRG, :C4vCTMRG, :SymmetricBoundaryMPS}
+
+# Boundary MPS
+const boundarymps_mps_alg = :VUMPS # ∈ {:VUMPS, :VOMPS}
+
 # SVD forward & reverse
 const trunc = :FixedSpaceTruncation # ∈ {:FixedSpaceTruncation, :notrunc, :truncerror, :truncspace, :trunctol}
 const rrule_degeneracy_atol = 1.0e-13

src/PEPSKit.jl

@@ -84,6 +84,7 @@ include("operators/models.jl")
 
 include("environments/ctmrg_environments.jl")
 include("environments/vumps_environments.jl")
+include("environments/boundarymps_environments.jl")
 include("environments/suweight.jl")
 include("environments/bp_environments.jl")
 include("environments/product_state_environments.jl")
@@ -105,6 +106,7 @@ include("algorithms/contractions/absorb.jl")
 include("algorithms/contractions/absorb_weight.jl")
 include("algorithms/contractions/transfer.jl")
 include("algorithms/contractions/vumps_contractions.jl")
+include("algorithms/contractions/boundarymps_contractions.jl")
 include("algorithms/contractions/bp_messages.jl")
 include("algorithms/contractions/local_patch/expr_utils.jl")
 include("algorithms/contractions/local_patch/network_expr.jl")
@@ -120,6 +122,8 @@ include("algorithms/contractions/correlator/peps.jl")
 include("algorithms/contractions/correlator/pepo_purified.jl")
 include("algorithms/contractions/correlator/pepo_1layer.jl")
 
+include("algorithms/boundary_algorithm.jl")
+
 include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/projectors/projectors.jl")
@@ -161,6 +165,10 @@ include("algorithms/expectation_value/correlator_adapters.jl")
 include("algorithms/expectation_value/correlators.jl")
 include("algorithms/toolbox.jl")
 
+include("algorithms/boundarymps/symmetric_boundarymps.jl")
+include("algorithms/boundarymps/characteristic_equations.jl")
+include("algorithms/boundarymps/observables.jl")
+
 include("algorithms/optimization/implicit_differentiation.jl")
 include("algorithms/optimization/preconditioning.jl")
 include("algorithms/optimization/peps_optimization.jl")
@@ -177,6 +185,7 @@ export corner, edge, setcorner!, setedge!
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
+export SymmetricBoundaryMPS, SymmetricBoundaryMPSEnv
 export initialize_random_c4v_env, initialize_singlet_c4v_env
 export LocalOperator, physicalspace
 export product_peps

src/algorithms/boundary_algorithm.jl

@@ -0,0 +1,45 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for all algorithms that contract an infinite square network by computing
+a boundary fixed point, such as CTMRG and boundary MPS algorithms.
+
+This is the type of the `boundary_alg` field of [`PEPSOptimize`](@ref), and hence the
+supertype every algorithm which can be used to contract a network during a variational
+optimization must belong to.
+"""
+abstract type BoundaryAlgorithm end
+
+const BOUNDARY_ALGORITHM_SYMBOLS = IdDict{Symbol, Type{<:BoundaryAlgorithm}}()
+
+"""
+    _tol(alg::BoundaryAlgorithm)
+
+Effective convergence tolerance of a boundary contraction algorithm. Defaults to the `tol`
+field; algorithms which keep their tolerance elsewhere (such as boundary MPS algorithms,
+where it lives on the wrapped MPS optimization algorithm) must overload this.
+"""
+_tol(alg::BoundaryAlgorithm) = alg.tol
+
+"""
+    BoundaryAlgorithm(; alg=:$(Defaults.boundary_alg), kwargs...)
+
+Keyword argument parser returning the appropriate [`BoundaryAlgorithm`](@ref) struct, where
+`alg` selects the boundary contraction *family*:
+
+* `:SimultaneousCTMRG`, `:SequentialCTMRG`, `:C4vCTMRG` : dispatch to [`CTMRGAlgorithm`](@ref)
+* `:SymmetricBoundaryMPS` : dispatch to [`SymmetricBoundaryMPS`](@ref)
+
+All remaining keyword arguments are forwarded to the corresponding parser. Note that the
+underlying MPS optimization algorithm of a boundary MPS contraction is *not* selected here,
+but through that parser's own `mps_alg` keyword, since several boundary MPS families can be
+driven by the same MPS algorithm.
+"""
+function BoundaryAlgorithm(; alg = Defaults.boundary_alg, kwargs...)
+    # CTMRG variants go through their own parser, which handles the projector and
+    # decomposition keyword arguments
+    haskey(CTMRG_SYMBOLS, alg) && return CTMRGAlgorithm(; alg, kwargs...)
+    haskey(BOUNDARY_ALGORITHM_SYMBOLS, alg) ||
+        throw(ArgumentError("unknown boundary algorithm: $alg"))
+    return BOUNDARY_ALGORITHM_SYMBOLS[alg](; kwargs...)
+end

src/algorithms/boundarymps/characteristic_equations.jl

@@ -0,0 +1,68 @@
+#
+# Characteristic equation used in implicit differentiation of symmetric boundary MPS
+# contractions
+#
+
+"""
+    generate_boundary_mps_characteristic_equation(
+        env::SymmetricBoundaryMPSEnv, (VLfp, VRfp)::Tuple{<:EdgeTensor, <:RightProjector}
+    )
+
+Takes the fixed-point values of a converged symmetric boundary MPS contraction, given by the
+environment `env`, along with the left null space `VLfp` of its left-gauged MPS tensor and
+the right null space `VRfp` of its right-gauged MPS tensor, and generates a function
+``F(s, l, r, c, gl, gr)`` which characterizes the convergence of the boundary MPS contraction
+in terms of the characteristic equation ``F(s, l, r, c, gl, gr) = 0``.
+
+Here, ``s`` corresponds to a state variable (e.g. an `InfinitePEPS` that is being optimized),
+``c``, ``gl`` and ``gr`` directly represent the MPS bond tensor and the left and right
+environments, while ``l`` and ``r`` parametrize differentiable left- and
+right-gauged MPS tensors as ``AL = AL_{fp} + V_{L,fp} * l`` and
+``AR = AR_{fp} + (r * V_{R,fp})``.
+
+``F`` returns a tuple of five tensors, corresponding to an equation for each of ``l``, ``r``,
+``c``, ``gl`` and ``gr``. The first three are obtained by projecting the effective site
+operator applied to the center-gauged MPS tensor onto the respective tangent directions, and
+the last two express that the environments are fixed points of the MPS-network-MPS transfer
+matrix. All equations are normalized by the eigenvalue ``λ`` of the effective site operator.
+
+See also [`generate_symmetric_characteristic_equation`](@ref).
+"""
+function generate_boundary_mps_characteristic_equation(
+        env::SymmetricBoundaryMPSEnv, (VLfp, VRfp)::Tuple{TE, TP}
+    ) where {TE <: EdgeTensor, TP <: RightProjector}
+    ALfp, ARfp, Cfp, = _unpack(env)
+
+    # constant preconditioner
+    # NOTE: this relies on the bond tensor having been diagonalized, which
+    # `leading_boundary` takes care of
+    iCfp = sdiag_pow(real(DiagonalTensorMap(Cfp)), -1)
+
+    function boundary_mps_characteristic_equation(state, l, r, c, gl, gr)
+        network = InfiniteSquareNetwork(state)
+        O = network[1, 1]
+
+        # prepare appropriately parametrized MPS tensors
+        AL = ALfp + VLfp * l
+        AR = ARfp + repartition_left(r * VRfp)
+        ARR = repartition_right(AR) # 'right isometry' form
+
+        # construct the center tensor in a symmetric way
+        AC = (absorb_right(AL, c) + absorb_left(AR, c)) / 2
+
+        # main partial contractions to reuse
+        AC´L = ∂AC(AC, gl, O, gr) # 'left isometry' form
+        AC´R = repartition_right(AC´L) # 'right isometry' form
+        λ = dot(AC, AC´L)
+
+        F1 = VLfp' * AC´L * iCfp / λ - l
+        F2 = iCfp * AC´R * VRfp' / λ - r
+        F3 = (AL' * AC´L + AC´R * ARR') / (2 * λ) - c
+        F4 = MPSKit.transfer_left(gl, O, AL, AL) / λ - gl
+        F5 = MPSKit.transfer_right(gr, O, AR, AR) / λ - gr
+
+        return F1, F2, F3, F4, F5
+    end
+
+    return boundary_mps_characteristic_equation
+end

src/algorithms/boundarymps/observables.jl

@@ -0,0 +1,34 @@
+#
+# Observables evaluated using a symmetric boundary MPS environment
+#
+
+"""
+    network_value(network::InfiniteSquareNetwork, env::SymmetricBoundaryMPSEnv)
+
+Return the value (per unit cell) of a contractible network contracted using a symmetric
+boundary MPS environment.
+"""
+function network_value(network::InfiniteSquareNetwork, env::SymmetricBoundaryMPSEnv)
+    size(network) == (1, 1) ||
+        throw(ArgumentError("symmetric boundary MPS environments require a single-site unit cell"))
+    AC = get_AC(env)
+    AC´ = PEPS_AC_Hamiltonian(env.GL, network[1, 1], env.GR) * AC
+    return dot(AC, AC´)
+end
+function network_value(state, env::SymmetricBoundaryMPSEnv)
+    return network_value(InfiniteSquareNetwork(state), env)
+end
+
+function LinearAlgebra.norm(peps::InfinitePEPS, env::SymmetricBoundaryMPSEnv)
+    return network_value(InfiniteSquareNetwork(peps), env)
+end
+
+## Partition function tensor insertions
+
+function contract_local_tensor(
+        ::Union{CartesianIndex{2}, Tuple{Int, Int}}, O::PartitionFunctionTensor,
+        env::SymmetricBoundaryMPSEnv,
+    )
+    # the index is irrelevant here: symmetric boundary MPS environments are single-site
+    return _contract_site(get_AC(env), env.GL, env.GR, O)
+end