Skip to content

Commit 81be8c0

Browse files
pbrehmerYue-Zhengyuanlkdvos
authored
Add C4v CTMRG (#321)
* Add `C4vCTMRG` contraction algorithm * To that end, implement a `EighAdjoint` struct that wraps around MatrixAlgebraKit's `eigh_trunc!` * Add all required rrules to make the C4v CTMRG algorithm differentiable * Add a dedicated C4v gauge fixing routine * Add a couple of tests for the contraction, gauge fixing and optimization --------- Co-authored-by: Yue Zhengyuan <yuezy1997@icloud.com> Co-authored-by: Lukas Devos <ldevos98@gmail.com>
1 parent cf49bda commit 81be8c0

26 files changed

Lines changed: 1196 additions & 147 deletions

src/Defaults.jl

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,9 +101,20 @@ const svd_rrule_alg = :full # ∈ {:full, :gmres, :bicgstab, :arnoldi}
101101
const svd_rrule_broadening = 1.0e-13
102102
const krylovdim_factor = 1.4
103103

104+
# eigh forward & reverse
105+
const eigh_fwd_alg = :qriteration # ∈ {:qriteration, :bisection, :divideandconquer, :multiple, :lanczos, :blocklanczos}
106+
const eigh_rrule_alg = :trunc # ∈ {:trunc, :full}
107+
const eigh_rrule_verbosity = 0
108+
109+
# QR forward & reverse
110+
# const qr_fwd_alg = :something # TODO
111+
# const qr_rrule_alg = :something
112+
# const qr_rrule_verbosity = :something
113+
104114
# Projectors
105115
const projector_alg = :halfinfinite # ∈ {:halfinfinite, :fullinfinite}
106116
const projector_verbosity = 0
117+
const projector_alg_c4v = :c4v_eigh # ∈ {:c4v_eigh, :c4v_qr (TODO)}
107118

108119
# Fixed-point gradient
109120
const gradient_tol = 1.0e-6

src/PEPSKit.jl

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,24 @@ import VectorInterface as VI
99

1010
using MatrixAlgebraKit
1111
using MatrixAlgebraKit: TruncationStrategy, LAPACK_DivideAndConquer, LAPACK_QRIteration
12+
using MatrixAlgebraKit: NoTruncation, LAPACK_EighAlgorithm, truncate
13+
using MatrixAlgebraKit: eigh_pullback!, eigh_trunc_pullback!, findtruncated, diagview
14+
1215
using TensorKit
16+
using TensorKit: AdjointTensorMap, SectorDict
17+
using TensorKit: throw_invalid_innerproduct, similarstoragetype
18+
using TensorKit.Factorizations: TruncationSpace, _notrunc_ind
19+
20+
using KrylovKit
21+
using KrylovKit: Lanczos, BlockLanczos
1322

14-
using KrylovKit, OptimKit, TensorOperations
23+
using TensorOperations, OptimKit
1524
using ChainRulesCore, Zygote
1625
using LoggingExtras
1726

1827
using MPSKit
1928
using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, ProductTransferMatrix
29+
using MPSKit: InfiniteEnvironments
2030
import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
2131
import MPSKit: infinite_temperature_density_matrix
2232

@@ -29,6 +39,7 @@ include("Defaults.jl") # Include first to allow for docstring interpolation wit
2939

3040
include("utility/util.jl")
3141
include("utility/diffable_threads.jl")
42+
include("utility/eigh.jl")
3243
include("utility/svd.jl")
3344
include("utility/rotations.jl")
3445
include("utility/hook_pullback.jl")
@@ -42,6 +53,8 @@ include("networks/infinitesquarenetwork.jl")
4253
include("states/infinitepeps.jl")
4354
include("states/infinitepartitionfunction.jl")
4455

56+
include("utility/symmetrization.jl")
57+
4558
include("operators/infinitepepo.jl")
4659
include("operators/transfermatrix.jl")
4760
include("operators/localoperator.jl")
@@ -81,6 +94,7 @@ include("algorithms/ctmrg/projectors.jl")
8194
include("algorithms/ctmrg/simultaneous.jl")
8295
include("algorithms/ctmrg/sequential.jl")
8396
include("algorithms/ctmrg/gaugefix.jl")
97+
include("algorithms/ctmrg/c4v.jl")
8498

8599
include("algorithms/truncation/truncationschemes.jl")
86100
include("algorithms/truncation/fullenv_truncation.jl")
@@ -99,8 +113,6 @@ include("algorithms/transfermatrix.jl")
99113
include("algorithms/toolbox.jl")
100114
include("algorithms/correlators.jl")
101115

102-
include("utility/symmetrization.jl")
103-
104116
include("algorithms/optimization/fixed_point_differentiation.jl")
105117
include("algorithms/optimization/peps_optimization.jl")
106118

@@ -112,6 +124,8 @@ export SVDAdjoint, FullSVDReverseRule, IterSVD
112124
export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
113125
export FixedSpaceTruncation, SiteDependentTruncation
114126
export HalfInfiniteProjector, FullInfiniteProjector
127+
export EighAdjoint, IterEigh, C4vCTMRG, C4vEighProjector, C4vQRProjector
128+
export initialize_random_c4v_env, initialize_singlet_c4v_env
115129
export LocalOperator, physicalspace
116130
export product_peps
117131
export reduced_densitymatrix, expectation_value, network_value, cost_function

src/algorithms/ctmrg/c4v.jl

Lines changed: 238 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,238 @@
1+
"""
2+
$(TYPEDEF)
3+
4+
CTMRG algorithm assuming a C₄ᵥ-symmetric PEPS, i.e. invariance under 90° spatial rotation and
5+
Hermitian reflection. This requires a single-site unit cell. The projector is obtained from
6+
`eigh` decomposing the Hermitian enlarged corner.
7+
8+
## Fields
9+
10+
$(TYPEDFIELDS)
11+
12+
## Constructors
13+
14+
C4vCTMRG(; kwargs...)
15+
16+
Construct a C₄ᵥ CTMRG algorithm struct based on keyword arguments.
17+
For a full description, see [`leading_boundary`](@ref). The supported keywords are:
18+
19+
* `tol::Real=$(Defaults.ctmrg_tol)`
20+
* `maxiter::Int=$(Defaults.ctmrg_maxiter)`
21+
* `miniter::Int=$(Defaults.ctmrg_miniter)`
22+
* `verbosity::Int=$(Defaults.ctmrg_verbosity)`
23+
* `trunc::Union{TruncationStrategy,NamedTuple}=(; alg::Symbol=:$(Defaults.trunc))`
24+
* `decomposition_alg::Union{<:EighAdjoint,NamedTuple}`
25+
* `projector_alg::Symbol=:$(Defaults.projector_alg_c4v)`
26+
"""
27+
struct C4vCTMRG{P <: ProjectorAlgorithm} <: CTMRGAlgorithm
28+
tol::Float64
29+
maxiter::Int
30+
miniter::Int
31+
verbosity::Int
32+
projector_alg::P
33+
end
34+
function C4vCTMRG(; kwargs...)
35+
return CTMRGAlgorithm(; alg = :c4v, kwargs...)
36+
end
37+
CTMRG_SYMBOLS[:c4v] = C4vCTMRG
38+
39+
"""
40+
$(TYPEDEF)
41+
42+
Projector algorithm implementing the `eigh` decomposition of a Hermitian enlarged corner.
43+
44+
## Fields
45+
46+
$(TYPEDFIELDS)
47+
48+
## Constructors
49+
50+
C4vEighProjector(; kwargs...)
51+
52+
Construct the C₄ᵥ `eigh`-based projector algorithm based on the following keyword arguments:
53+
54+
* `decomposition_alg::Union{<:EighAdjoint,NamedTuple}=EighAdjoint()` : `eigh` algorithm including the reverse rule. See [`EighAdjoint`](@ref).
55+
* `trunc::Union{TruncationStrategy,NamedTuple}=(; alg::Symbol=:$(Defaults.trunc))` : Truncation strategy for the projector computation, which controls the resulting virtual spaces. Here, `alg` can be one of the following:
56+
- `:fixedspace` : Keep virtual spaces fixed during projection
57+
- `:notrunc` : No singular values are truncated and the performed SVDs are exact
58+
- `:truncerror` : Additionally supply error threshold `η`; truncate to the maximal virtual dimension of `η`
59+
- `:truncrank` : Additionally supply truncation dimension `η`; truncate such that the 2-norm of the truncated values is smaller than `η`
60+
- `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space
61+
- `:trunctol` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η`
62+
* `verbosity::Int=$(Defaults.projector_verbosity)` : Projector output verbosity which can be:
63+
0. Suppress output information
64+
1. Print singular value degeneracy warnings
65+
"""
66+
struct C4vEighProjector{S <: EighAdjoint, T} <: ProjectorAlgorithm
67+
decomposition_alg::S
68+
trunc::T
69+
verbosity::Int
70+
end
71+
function C4vEighProjector(; kwargs...)
72+
return ProjectorAlgorithm(; alg = :c4v_eigh, kwargs...)
73+
end
74+
PROJECTOR_SYMBOLS[:c4v_eigh] = C4vEighProjector
75+
76+
# struct C4vQRProjector{S, T} <: ProjectorAlgorithm
77+
# decomposition_alg::S
78+
# verbosity::Int
79+
# end
80+
# function C4vQRProjector(; kwargs...)
81+
# return ProjectorAlgorithm(; alg = :c4v_qr, kwargs...)
82+
# end
83+
# PROJECTOR_SYMBOLS[:c4v_qr] = C4vQRProjector
84+
85+
#
86+
## C4v-symmetric CTMRG iteration (called through `leading_boundary`)
87+
#
88+
89+
function ctmrg_iteration(
90+
network,
91+
env::CTMRGEnv,
92+
alg::C4vCTMRG,
93+
)
94+
enlarged_corner = c4v_enlarge(network, env, alg.projector_alg)
95+
corner′, projector, info = c4v_projector(enlarged_corner, alg.projector_alg)
96+
edge′ = c4v_renormalize(network, env, projector)
97+
return CTMRGEnv(corner′, edge′), info
98+
end
99+
100+
"""
101+
c4v_enlarge(network, env, ::C4vEighProjector)
102+
103+
Compute the normalized and Hermitian-symmetrized C₄ᵥ enlarged corner.
104+
"""
105+
function c4v_enlarge(network, env, ::C4vEighProjector)
106+
enlarged_corner = TensorMap(EnlargedCorner(network, env, (NORTHWEST, 1, 1)))
107+
return 0.5 * (enlarged_corner + enlarged_corner') / norm(enlarged_corner)
108+
end
109+
# function c4v_enlarge(enlarged_corner, alg::C4vQRProjector)
110+
# # TODO
111+
# end
112+
113+
"""
114+
c4v_projector(enlarged_corner, alg::C4vEighProjector)
115+
116+
Compute the C₄ᵥ projector from `eigh` decomposing the Hermitian `enlarged_corner`.
117+
Also return the normalized eigenvalues as the new corner tensor.
118+
"""
119+
function c4v_projector(enlarged_corner, alg::C4vEighProjector)
120+
trunc = truncation_strategy(alg, enlarged_corner)
121+
D, V, info = eigh_trunc!(enlarged_corner, decomposition_algorithm(alg); trunc)
122+
123+
# Check for degenerate eigenvalues
124+
Zygote.isderiving() && ignore_derivatives() do
125+
if alg.verbosity > 0 && is_degenerate_spectrum(D)
126+
vals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(D))
127+
@warn("degenerate eigenvalues detected: ", vals)
128+
end
129+
end
130+
131+
return D / norm(D), V, (; D, V, info...)
132+
end
133+
# function c4v_projector(enlarged_corner, alg::C4vQRProjector)
134+
# # TODO
135+
# end
136+
137+
"""
138+
c4v_renormalize(network, env, projector)
139+
140+
Renormalize the single edge tensor.
141+
"""
142+
function c4v_renormalize(network, env, projector)
143+
new_edge = renormalize_north_edge(env.edges[1], projector, projector', network[1, 1])
144+
new_edge = _project_hermitian(new_edge) # additional Hermitian projection step for numerical stability
145+
return new_edge / norm(new_edge)
146+
end
147+
148+
# TODO: this should eventually be the constructor for a new C4vCTMRGEnv type
149+
function CTMRGEnv(
150+
corner::AbstractTensorMap{T, S, 1, 1}, edge::AbstractTensorMap{T′, S, N, 1}
151+
) where {T, T′, S, N}
152+
corners = fill(corner, 4, 1, 1)
153+
edge_SW = physical_flip(edge)
154+
edges = reshape([edge, edge, edge_SW, edge_SW], (4, 1, 1))
155+
return CTMRGEnv(corners, edges)
156+
end
157+
158+
#
159+
## utility
160+
#
161+
162+
# Adjoint of an edge tensor, but permutes the physical spaces back into the codomain.
163+
# Intuitively, this conjugates a tensor and then reinterprets its 'direction' as an edge tensor.
164+
function _dag(A::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
165+
return permute(A', ((1, (3:(N + 1))...), (2,)))
166+
end
167+
168+
function physical_flip(A::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
169+
return flip(A, 2:N)
170+
end
171+
172+
# call it `_project_hermitian` to avoid type piracy with MAK's exported project_hermitian
173+
function _project_hermitian(E::AbstractTensorMap{T, S, N, 1}) where {T, S, N}
174+
= (E + physical_flip(_dag(E))) / 2
175+
return
176+
end
177+
function _project_hermitian(C::AbstractTensorMap{T, S, 1, 1}) where {T, S}
178+
= (C + C') / 2
179+
return
180+
end
181+
182+
# should perform this check at the beginning of `leading_boundary` really...
183+
function check_symmetry(state, ::RotateReflect; atol = 1.0e-10)
184+
@assert length(state) == 1 "check_symmetry only works for single site unit cells"
185+
@assert norm(state[1] - _fit_spaces(rotl90(state[1]), state[1])) /
186+
norm(state[1]) < atol "not rotation invariant"
187+
@assert norm(state[1] - _fit_spaces(herm_depth(state[1]), state[1])) /
188+
norm(state[1]) < atol "not hermitian-reflection invariant"
189+
return nothing
190+
end
191+
192+
#
193+
## environment initialization
194+
#
195+
196+
"""
197+
initialize_random_c4v_env([f=randn, T=scalartype(state)], state, Venv::ElementarySpace)
198+
199+
Initialize a C₄ᵥ-symmetric `CTMRGEnv` on virtual spaces `Venv` with random entries created
200+
by `f` and scalartype `T`.
201+
"""
202+
function initialize_random_c4v_env(state, Venv::ElementarySpace)
203+
return initialize_random_c4v_env(randn, scalartype(state), state, Venv)
204+
end
205+
function initialize_random_c4v_env(f, T, state::InfinitePEPS, Venv::ElementarySpace)
206+
Vpeps = north_virtualspace(state, 1, 1)'
207+
return initialize_random_c4v_env(f, T, Vpeps Vpeps', Venv)
208+
end
209+
function initialize_random_c4v_env(f, T, state::InfinitePartitionFunction, Venv::ElementarySpace)
210+
Vpf = north_virtualspace(state, 1, 1)'
211+
return initialize_random_c4v_env(f, T, Vpf, Venv)
212+
end
213+
function initialize_random_c4v_env(f, T, Vstate::VectorSpace, Venv::ElementarySpace)
214+
corner₀ = DiagonalTensorMap(randn(real(T), Venv Venv))
215+
edge₀ = f(T, Venv Vstate Venv)
216+
edge₀ = _project_hermitian(edge₀)
217+
return CTMRGEnv(corner₀, edge₀)
218+
end
219+
220+
"""
221+
initialize_singlet_c4v_env([T=scalartype(state)], state::InfinitePEPS, Venv::ElementarySpace)
222+
223+
Initialize a C₄ᵥ-symmetric `CTMRGEnv` with a singlet corner of dimension `dim(Venv)` and an
224+
identity edge from `id(T, Venv ⊗ Vpeps)`.
225+
"""
226+
function initialize_singlet_c4v_env(state::InfinitePEPS, Venv::ElementarySpace)
227+
return initialize_singlet_c4v_env(scalartype(state), state, Venv)
228+
end
229+
function initialize_singlet_c4v_env(T, state::InfinitePEPS, Venv::ElementarySpace)
230+
Vpeps = north_virtualspace(state, 1, 1)'
231+
return initialize_singlet_c4v_env(T, Vpeps, Venv)
232+
end
233+
function initialize_singlet_c4v_env(T, Vpeps::ElementarySpace, Venv::ElementarySpace)
234+
corner₀ = DiagonalTensorMap(zeros(real(T), Venv Venv))
235+
corner₀.data[1] = one(real(T))
236+
edge₀ = permute(id(T, Venv Vpeps), ((1, 2, 4), (3,)))
237+
return CTMRGEnv(corner₀, edge₀)
238+
end

src/algorithms/ctmrg/ctmrg.jl

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -19,17 +19,19 @@ function CTMRGAlgorithm(;
1919
maxiter = Defaults.ctmrg_maxiter, miniter = Defaults.ctmrg_miniter,
2020
verbosity = Defaults.ctmrg_verbosity,
2121
trunc = (; alg = Defaults.trunc),
22-
svd_alg = (;),
22+
decomposition_alg = (;),
2323
projector_alg = Defaults.projector_alg, # only allows for Symbol/NamedTuple to expose projector kwargs
2424
)
2525
# replace symbol with projector alg type
2626
haskey(CTMRG_SYMBOLS, alg) || throw(ArgumentError("unknown CTMRG algorithm: $alg"))
2727
alg_type = CTMRG_SYMBOLS[alg]
2828

2929
# parse CTMRG projector algorithm
30-
30+
if alg == :c4v && projector_alg == Defaults.projector_alg
31+
projector_alg = Defaults.projector_alg_c4v
32+
end
3133
projector_algorithm = ProjectorAlgorithm(;
32-
alg = projector_alg, svd_alg, trunc, verbosity
34+
alg = projector_alg, decomposition_alg, trunc, verbosity
3335
)
3436

3537
return alg_type(tol, maxiter, miniter, verbosity, projector_algorithm)
@@ -64,8 +66,9 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st
6466
3. Iteration info
6567
4. Debug info
6668
* `alg::Symbol=:$(Defaults.ctmrg_alg)` : Variant of the CTMRG algorithm. See also [`CTMRGAlgorithm`](@ref).
67-
- `:simultaneous`: Simultaneous expansion and renormalization of all sides.
68-
- `:sequential`: Sequential application of left moves and rotations.
69+
- `:simultaneous` : Simultaneous expansion and renormalization of all sides.
70+
- `:sequential` : Sequential application of left moves and rotations.
71+
- `:c4v` : CTMRG assuming C₄ᵥ-symmetric PEPS and environment.
6972
7073
### Projector algorithm
7174
@@ -76,10 +79,11 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st
7679
- `:truncrank` : Additionally supply truncation dimension `η`; truncate such that the 2-norm of the truncated values is smaller than `η`
7780
- `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space
7881
- `:trunctol` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η`
79-
* `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.
82+
* `decomposition_alg` : Tensor decomposition algorithm for computing projectors. See e.g. [`SVDAdjoint`](@ref).
8083
* `projector_alg::Symbol=:$(Defaults.projector_alg)` : Variant of the projector algorithm. See also [`ProjectorAlgorithm`](@ref).
8184
- `:halfinfinite` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments.
8285
- `:fullinfinite` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments.
86+
- `:c4v_eigh` : Projection via `eigh` of the Hermitian enlarged corner.
8387
8488
## Return values
8589
@@ -101,6 +105,8 @@ set of vectors and values will be returned as well:
101105
* `U_full` : Last unit cell of all left singular vectors.
102106
* `S_full` : Last unit cell of all singular values.
103107
* `V_full` : Last unit cell of all right singular vectors.
108+
109+
For `C4vCTMRG` instead the last eigendecomposition `V` and `D` (and `V_full`, `D_full`) will be returned.
104110
"""
105111
function leading_boundary(env₀::CTMRGEnv, network::InfiniteSquareNetwork; kwargs...)
106112
alg = select_algorithm(leading_boundary, env₀; kwargs...)

0 commit comments

Comments
 (0)