C4v QR-CTMRG

src/Defaults.jl

@@ -107,9 +107,10 @@ const eigh_rrule_alg = :trunc # ∈ {:trunc, :full}
 const eigh_rrule_verbosity = 0
 
 # QR forward & reverse
-# const qr_fwd_alg = :something # TODO
-# const qr_rrule_alg = :something
-# const qr_rrule_verbosity = :something
+const qr_fwd_alg = :qr
+const qr_fwd_positive = true
+const qr_rrule_alg = :qr
+const qr_rrule_verbosity = 0
 
 # Projectors
 const projector_alg = :halfinfinite # ∈ {:halfinfinite, :fullinfinite}

src/PEPSKit.jl

@@ -42,6 +42,7 @@ 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")
@@ -70,6 +71,7 @@ include("environments/bp_environments.jl")
 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/projector.jl")
 include("algorithms/contractions/ctmrg/halfinf_env.jl")
 include("algorithms/contractions/ctmrg/fullinf_env.jl")
@@ -125,7 +127,7 @@ export SVDAdjoint, FullSVDReverseRule, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
-export EighAdjoint, IterEigh, C4vCTMRG, C4vEighProjector, C4vQRProjector
+export EighAdjoint, IterEigh, QRAdjoint, C4vCTMRG, C4vEighProjector, C4vQRProjector
 export initialize_random_c4v_env, initialize_singlet_c4v_env
 export LocalOperator, physicalspace
 export product_peps

src/algorithms/contractions/ctmrg/enlarge_corner_column.jl

@@ -0,0 +1,24 @@
+# Column-enlarged corner contractions
+# ----------------------------
+
+# Northwest corner
+# ----------------
+"""
+$(SIGNATURES)
+
+Contract the half-enlarged northwest corner of the CTMRG environment.
+```
+    C₁-←-E₁-←-
+    ↓    |
+```
+"""
+@generated function column_enlarge_northwest_corner(
+        C_northwest::CTMRGCornerTensor, edge::CTMRGEdgeTensor{T, S, N}
+    ) where {T, S, N}
+    CE_e = tensorexpr(:CE, -(1:N), -(N + 1))
+    C_e = tensorexpr(:C_northwest, -1, 1)
+    E_e = tensorexpr(:edge, (1, -(2:N)...), -(N + 1))
+    return macroexpand(@__MODULE__, :(return @tensor $CE_e := $C_e * $E_e))
+end
+
+# TODO: Other column-enlarged corners when QR-CTMRG for arbitrary unit cell is possible

src/algorithms/ctmrg/c4v.jl

@@ -73,27 +73,88 @@ function C4vEighProjector(; kwargs...)
 end
 PROJECTOR_SYMBOLS[:c4v_eigh] = C4vEighProjector
 
-# struct C4vQRProjector{S, T} <: ProjectorAlgorithm
-#     decomposition_alg::S
-#     verbosity::Int
-# end
-# function C4vQRProjector(; kwargs...)
-#     return ProjectorAlgorithm(; alg = :c4v_qr, kwargs...)
-# end
-# PROJECTOR_SYMBOLS[:c4v_qr] = C4vQRProjector
+"""
+$(TYPEDEF)
+
+Projector algorithm implementing the `qr` decomposition of a column-enlarged corner.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+## Constructors
+
+    C4vQRProjector(; kwargs...)
+
+Construct the C₄ᵥ `qr`-based projector algorithm
+based on the following keyword arguments:
+
+* `decomposition_alg=QRAdjoint()` : `left_orth` algorithm including the reverse rule. See [`QRAdjoint`](@ref).
+"""
+struct C4vQRProjector{S} <: ProjectorAlgorithm
+    # TODO: support all `left_orth` algorithms
+    decomposition_alg::S
+end
+function C4vQRProjector(; kwargs...)
+    return ProjectorAlgorithm(; alg = :c4v_qr, kwargs...)
+end
+PROJECTOR_SYMBOLS[:c4v_qr] = C4vQRProjector
+
+# no truncation
+_set_truncation(alg::C4vQRProjector, ::TruncationStrategy) = alg
+
+function check_input(
+        network::InfiniteSquareNetwork, env::CTMRGEnv, alg::C4vCTMRG; atol = 1.0e-10
+    )
+    # check unit cell size
+    length(network) == 1 || throw(ArgumentError("C4v CTMRG is only compatible with single-site unit cells."))
+    O = network[1, 1]
+    # check for fermionic braiding statistics
+    BraidingStyle(sectortype(spacetype(network))) != Bosonic() &&
+        throw(ArgumentError("C4v CTMRG is currently only implemented for networks consisting of tensors with bosonic braiding."))
+    # check sufficient condition for equality and duality of spaces that we assume
+    west_virtualspace(O) == _elementwise_dual(north_virtualspace(O)) ||
+        throw(ArgumentError("C4v CTMRG requires south and west virtual space to be the dual of north and east virtual space."))
+    # check rotation invariance of the local tensors, with the exact spaceflips we assume
+    _isapprox_localsandwich(O, flip_virtualspace(_rotl90_localsandwich(O), [EAST, WEST]); atol = atol) ||
+        throw(ArgumentError("C4v CTMRG requires the local tensors to be invariant under 90° rotation."))
+    # check the hermitian reflection invariance of the local tensors, with the exact spaceflips we assume
+    _isapprox_localsandwich(
+        O,
+        flip_physicalspace(flip_virtualspace(herm_depth(O), [EAST, WEST]));
+        atol = atol
+    ) ||
+        throw(ArgumentError("C4v CTMRG requires the local tensors to be invariant under hermitian reflection."))
+    # TODO: check compatibility of network and environment spaces in general?
+    return nothing
+end
 
 #
 ## C4v-symmetric CTMRG iteration (called through `leading_boundary`)
 #
 
+function ctmrg_iteration(network, env::CTMRGEnv, ::C4vCTMRG{P}) where {P}
+    throw(ArgumentError("Unknown C4v projector algorithm $P"))
+end
 function ctmrg_iteration(
         network,
         env::CTMRGEnv,
-        alg::C4vCTMRG,
+        alg::C4vCTMRG{<:C4vEighProjector},
     )
     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)
+    corner′, projector, info = c4v_projector!(enlarged_corner, alg.projector_alg)
+    edge′ = c4v_renormalize_edge(network, env, projector)
+    return CTMRGEnv(corner′, edge′), info
+end
+function ctmrg_iteration(
+        network,
+        env::CTMRGEnv,
+        alg::C4vCTMRG{<:C4vQRProjector},
+    )
+    enlarged_corner = c4v_enlarge(env, alg.projector_alg)
+    projector, info = c4v_projector!(enlarged_corner, alg.projector_alg)
+    edge′ = c4v_renormalize_edge(network, env, projector)
+    corner′ = c4v_qr_renormalize_corner(edge′, projector, info.R)
     return CTMRGEnv(corner′, edge′), info
 end
 
@@ -104,19 +165,26 @@ Compute the normalized and Hermitian-symmetrized C₄ᵥ enlarged corner.
 """
 function c4v_enlarge(network, env, ::C4vEighProjector)
     enlarged_corner = TensorMap(EnlargedCorner(network, env, (NORTHWEST, 1, 1)))
-    return 0.5 * (enlarged_corner + enlarged_corner') / norm(enlarged_corner)
+    # TODO: replace by `project_hermitian`
+    enlarged_corner = 0.5 * (enlarged_corner + enlarged_corner')
+    return enlarged_corner / norm(enlarged_corner)
 end
-# function c4v_enlarge(enlarged_corner, alg::C4vQRProjector)
-#     # TODO
-# end
+"""
+    c4v_enlarge(env, ::C4vQRProjector)
 
+Compute the normalized column-enlarged northeast corner for C₄ᵥ QR-CTMRG.
 """
-    c4v_projector(enlarged_corner, alg::C4vEighProjector)
+function c4v_enlarge(env, ::C4vQRProjector)
+    return TensorMap(ColumnEnlargedCorner(env, (NORTHWEST, 1, 1)))
+end
+
+"""
+    c4v_projector!(enlarged_corner, alg::C4vEighProjector)
 
 Compute the C₄ᵥ projector from `eigh` decomposing the Hermitian `enlarged_corner`.
 Also return the normalized eigenvalues as the new corner tensor.
 """
-function c4v_projector(enlarged_corner, alg::C4vEighProjector)
+function c4v_projector!(enlarged_corner, alg::C4vEighProjector)
     trunc = truncation_strategy(alg, enlarged_corner)
     D, V, info = eigh_trunc!(enlarged_corner, decomposition_algorithm(alg); trunc)
 
@@ -130,21 +198,91 @@ function c4v_projector(enlarged_corner, alg::C4vEighProjector)
 
     return D / norm(D), V, (; D, V, info...)
 end
-# function c4v_projector(enlarged_corner, alg::C4vQRProjector)
-#     # TODO
-# end
+"""
+    c4v_projector!(enlarged_corner, alg::C4vQRProjector)
+
+Compute the C₄ᵥ projector by decomposing the column-enlarged corner with `left_orth`.
+```
+                   R--←--
+                   ↓
+    C-←-E-←-  =  [~Q~]    
+    ↓   |        ↓   |
+```
+"""
+function c4v_projector!(enlarged_corner, alg::C4vQRProjector)
+    Q, R = left_orth!(enlarged_corner, decomposition_algorithm(alg))
+    return Q, (; Q, R)
+end
 
 """
-    c4v_renormalize(network, env, projector)
+    c4v_renormalize_edge(network, env, projector)
 
 Renormalize the single edge tensor.
+```
+        |~~~|-←-E-←-|~~~|
+    -←--| P'|   |   | P |--←-
+        |~~~|---A---|~~~|
+                |
+```
 """
-function c4v_renormalize(network, env, projector)
+# TODO: possible missing twists for fermions
+function c4v_renormalize_edge(network, env, projector)
     new_edge = renormalize_north_edge(env.edges[1], projector, projector', network[1, 1])
-    new_edge = _project_hermitian(new_edge) # additional Hermitian projection step for numerical stability
+    # additional Hermitian projection step for numerical stability
+    new_edge = _project_hermitian(new_edge)
     return new_edge / norm(new_edge)
 end
 
+"""
+    c4v_qr_renormalize_corner(new_edge, projector, R)
+
+Renormalize the single corner tensor
+```
+    C-←-E-←-|~~~|
+    |   |   | P |-←-
+    E---A---|~~~|
+    |   |
+    [~P']
+      ↓
+```
+Using the already calculated QR decomposition
+```
+                   R--←--
+                   ↓
+    C-←-E-←-  =  [~P~]    
+    ↓   |        ↓   |
+```
+we rewrite the renormalized corner as
+```
+    R-←-|~~~|
+    ↓   | P |-←-
+    E′--|~~~|
+    ↓
+```
+which reuses the renormalized edge `E′` (`new_edge`).
+(Credit: https://github.com/qiyang-ustc/QRCTM/blob/dd160116c3d7b02076691ceaf0a9833511ae532d/heisenberg.py#L80)
+"""
+# TODO: possible missing twists for fermions
+function c4v_qr_renormalize_corner(new_edge::CTMRGEdgeTensor, projector, R)
+    # contract edge and R
+    edge′ = physical_flip(new_edge)
+    ER = edge′ * twistdual(R, 1)
+    # contract (edge, R) with projector
+    new_corner = contract_edges(ER, projector)
+    new_corner = _project_hermitian(new_corner)
+    return new_corner / norm(new_corner)
+end
+
+# auxilary function: contract two CTMRG edge tensors
+@generated function contract_edges(
+        EL::CTMRGEdgeTensor{T, S, N}, ER::CTMRGEdgeTensor{T, S, N}
+    ) where {T, S, N}
+    C´_e = tensorexpr(:C´, -1, -2)
+    EL_e = tensorexpr(:EL, (-1, (2:N)...), 1)
+    ER_e = tensorexpr(:ER, 1:N, -2)
+    return macroexpand(@__MODULE__, :(return @tensor $C´_e := $EL_e * $ER_e))
+end
+
 # TODO: this should eventually be the constructor for a new C4vCTMRGEnv type
 function CTMRGEnv(
         corner::AbstractTensorMap{T, S, 1, 1}, edge::AbstractTensorMap{T′, S, N, 1}
@@ -159,6 +297,8 @@ end
 ## utility
 #
 
+# TODO: re-examine these for fermions
+
 # 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}
@@ -179,15 +319,6 @@ function _project_hermitian(C::AbstractTensorMap{T, S, 1, 1}) where {T, S}
     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

src/algorithms/ctmrg/ctmrg.jl

@@ -44,6 +44,14 @@ Perform a single CTMRG iteration in which all directions are being grown and ren
 """
 function ctmrg_iteration(network, env, alg::CTMRGAlgorithm) end
 
+"""
+    check_input(network, env, alg::CTMRGAlgorithm)
+
+Check compatibility of a given network and environment with a specified CTMRG algorithm.
+"""
+function check_input(network, env, alg::CTMRGAlgorithm) end
+@non_differentiable check_input(args...)
+
 """
     leading_boundary(env₀, network; kwargs...) -> env, info
     # expert version:
@@ -115,6 +123,7 @@ end
 function leading_boundary(
         env₀::CTMRGEnv, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm
     )
+    check_input(network, env₀, alg)
     log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG"))
     return LoggingExtras.withlevel(; alg.verbosity) do
         env = deepcopy(env₀)

src/algorithms/ctmrg/gaugefix.jl

@@ -11,7 +11,7 @@ end
 function gauge_fix(alg::ProjectorAlgorithm, signs, info)
     decomposition_alg_fixed = gauge_fix(decomposition_algorithm(alg), signs, info)
     alg_fixed = @set alg.decomposition_alg = decomposition_alg_fixed # every ProjectorAlgorithm needs an `decomposition_alg` field?
-    alg_fixed = @set alg_fixed.trunc = notrunc()
+    alg_fixed = _set_truncation(alg_fixed, notrunc()) # potentially set no truncation
     return alg_fixed
 end
 

src/algorithms/ctmrg/projectors.jl

@@ -33,6 +33,11 @@ function ProjectorAlgorithm(;
         _alg_or_nt(QRAdjoint, decomposition_alg)
     end # TODO: how do we solve this in a proper way?
 
+    # qr-ctmrg does not need truncation or degeneracy checks
+    if alg in [:c4v_qr]
+        return alg_type(decomposition_algorithm)
+    end
+
     # parse truncation scheme
     truncation_strategy = if trunc isa TruncationStrategy
         trunc
@@ -84,6 +89,17 @@ function truncation_strategy(alg::ProjectorAlgorithm, edge)
     end
 end
 
+"""
+    _set_truncation(alg::ProjectorAlgorithm, trunc::TruncationStrategy)
+
+Update the truncation strategy of a given projector algorithm, keeping all other settings
+the same.
+"""
+function _set_truncation(alg::ProjectorAlgorithm, trunc::TruncationStrategy)
+    alg´ = @set alg.trunc = trunc
+    return alg´
+end
+
 """
 $(TYPEDEF)