Skip to content

Commit 4c5f078

Browse files
pbrehmerYue-Zhengyuanlkdvosleburgel
authored
Add C4v CTMRG and QR-CTMRG example (#336)
* Allow for real-valued C4v optimization with `:fixed` mode and add real-valued C4v Heisenberg test * Add dummy truncation error and condition number return to `c4v_projector` and add QR-CTMRG optimization tests * Add rendered C4v eaxmple * Fix typos * Apply suggestions Co-authored-by: Yue Zhengyuan <yuezy1997@icloud.com> * Rerender C4v example * Write Hamiltonian without `rmul!` Co-authored-by: Lukas Devos <ldevos98@gmail.com> * Update rendered C4v example --------- Co-authored-by: Yue Zhengyuan <yuezy1997@icloud.com> Co-authored-by: Lukas Devos <ldevos98@gmail.com> Co-authored-by: leburgel <lander.burgelman@gmail.com>
1 parent 98fb536 commit 4c5f078

8 files changed

Lines changed: 717 additions & 3 deletions

File tree

‎docs/make.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ mathengine = MathJax3(
4242

4343
# examples pages
4444
examples_optimization = joinpath.(
45-
["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard"], Ref("index.md")
45+
["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard", "c4v_ctmrg"], Ref("index.md")
4646
)
4747
examples_time_evolution = joinpath.(
4848
["heisenberg_su", "hubbard_su", "j1j2_su"], Ref("index.md")

‎docs/src/assets/pepskit.bib‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -147,3 +147,13 @@ @article{liu_gapless_2022
147147
url = {https://www.sciencedirect.com/science/article/pii/S2095927322001001},
148148
author = {Wen-Yuan Liu and Shou-Shu Gong and Yu-Bin Li and Didier Poilblanc and Wei-Qiang Chen and Zheng-Cheng Gu},
149149
}
150+
151+
@misc{zhang_accelerating_2025,
152+
title={Accelerating two-dimensional tensor network contractions using QR-decompositions},
153+
author={Yining Zhang and Qi Yang and Philippe Corboz},
154+
year={2025},
155+
eprint={2505.00494},
156+
archivePrefix={arXiv},
157+
primaryClass={cond-mat.str-el},
158+
url={https://arxiv.org/abs/2505.00494},
159+
}

‎docs/src/examples/c4v_ctmrg/index.md‎

Lines changed: 248 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,248 @@
1+
```@meta
2+
EditURL = "../../../../examples/c4v_ctmrg/main.jl"
3+
```
4+
5+
[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/PEPSKit.jl/gh-pages?filepath=dev/examples/c4v_ctmrg/main.ipynb)
6+
[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](https://nbviewer.jupyter.org/github/QuantumKitHub/PEPSKit.jl/blob/gh-pages/dev/examples/c4v_ctmrg/main.ipynb)
7+
[![](https://img.shields.io/badge/download-project-orange)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumKitHub/PEPSKit.jl/examples/tree/gh-pages/dev/examples/c4v_ctmrg)
8+
9+
10+
# C₄ᵥ CTMRG and QR-CTMRG
11+
12+
In this example we demonstrate specialized CTMRG variants that exploit the $C_{4v}$ point group
13+
symmetry of the physical model and PEPS at hand. This allows us to significantly reduce the
14+
computational cost of contraction and optimization. To that end, we will consider the Heisenberg
15+
Hamiltonian on the square lattice
16+
17+
```math
18+
H = \sum_{\langle i,j \rangle} \left ( J_x S^{x}_i S^{x}_j + J_y S^{y}_i S^{y}_j + J_z S^{z}_i S^{z}_j \right )
19+
```
20+
21+
We want to treat the model in the antiferromagnetic regime where the ground state exhibits
22+
bipartite sublattice structure. To be able to represent this ground state in a $C_{4v}$-invariant
23+
manner on a single-site unit cell, we perform a unitary sublattice rotation by setting
24+
$(J_x, J_y, J_z)=(-1, 1, -1)$.
25+
26+
Let's get started by seeding the RNG and doing the imports:
27+
28+
````julia
29+
using Random
30+
using TensorKit, PEPSKit
31+
Random.seed!(123456789);
32+
````
33+
34+
## Defining a specialized Hamiltonian for C₄ᵥ-symmetric PEPS
35+
36+
Since the model under consideration and its ground state are invariant under 90° rotation and
37+
Hermitian reflection, evaluating the expectation values of the horizontal and vertical energy
38+
contributions are exactly equivalent. This allows us to effectively halve the computational cost
39+
by evaluating only half of the terms and multiplying by 2. In practice, we implement this
40+
using a specialized [`LocalOperator`](@ref) that contains only the relevant terms:
41+
42+
````julia
43+
using MPSKitModels: S_xx, S_yy, S_zz
44+
45+
# Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term
46+
function heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)
47+
return heisenberg_XYZ_c4v(ComplexF64, Trivial, lattice; kwargs...)
48+
end
49+
function heisenberg_XYZ_c4v(
50+
T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;
51+
Jx = -1.0, Jy = 1.0, Jz = -1.0, spin = 1 // 2,
52+
)
53+
@assert size(lattice) == (1, 1) "only trivial unit cells supported by C4v-symmetric Hamiltonians"
54+
term =
55+
S_xx(T, S; spin = spin) * Jx +
56+
S_yy(T, S; spin = spin) * Jy +
57+
S_zz(T, S; spin = spin) * Jz
58+
spaces = fill(domain(term)[1], (1, 1))
59+
return LocalOperator( # horizontal and vertical contributions are identical
60+
spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term
61+
)
62+
end;
63+
````
64+
65+
## Initializing C₄ᵥ-invariant PEPSs and environments
66+
67+
In order to use $C_{4v}$-symmetric algorithms, it is of course crucial to use initial guesses
68+
with $C_{4v}$ symmetry. First, we create a real-valued random PEPS that we explicitly
69+
symmetrize using [`symmetrize!`](@ref) and the $C_{4v}$ symmetry [`RotateReflect`](@ref):
70+
71+
````julia
72+
symm = RotateReflect()
73+
D = 2
74+
T = Float64
75+
peps_random = InfinitePEPS(randn, T, ComplexSpace(2), ComplexSpace(D))
76+
peps₀ = symmetrize!(peps_random, symm);
77+
````
78+
79+
Initializing an $C_{4v}$-invariant environment is a bit more subtle and there is no one-size-fits-all
80+
solution. As a good starting point one can use the initialization function [`initialize_random_c4v_env`](@ref)
81+
(or also [`initialize_singlet_c4v_env`](@ref)) where we construct a diagonal corner with random
82+
real entries and a random Hermitian edge tensor.
83+
84+
````julia
85+
χ = 16
86+
env_random_c4v = initialize_random_c4v_env(peps₀, ComplexSpace(χ));
87+
````
88+
89+
Then contracting the PEPS using $C_{4v}$ CTMRG is as easy as just calling [`leading_boundary`](@ref)
90+
but passing the initial PEPS and environment as well as the `alg = :c4v` keyword argument:
91+
92+
````julia
93+
env₀, = leading_boundary(env_random_c4v, peps₀; alg = :c4v, tol = 1.0e-10);
94+
````
95+
96+
````
97+
[ Info: CTMRG init: obj = -1.430301957018e-02 err = 1.0000e+00
98+
[ Info: CTMRG conv 36: obj = +8.685181513863e+00 err = 6.8459865700e-11 time = 1.30 sec
99+
100+
````
101+
102+
## C₄ᵥ-symmetric optimization
103+
104+
We now take `peps₀` and `env₀` as a starting point for a gradient-based energy
105+
minimization where we contract using $C_{4v}$ CTMRG such that the energy gradient will also
106+
exhibit $C_{4v}$ symmetry. For that, we call `fixedpoint` and specify `alg = :c4v`
107+
as the boundary contraction algorithm:
108+
109+
````julia
110+
H = real(heisenberg_XYZ_c4v(InfiniteSquare())) # make Hamiltonian real-valued
111+
peps, env, E, = fixedpoint(
112+
H, peps₀, env₀; optimizer_alg = (; tol = 1.0e-4), boundary_alg = (; alg = :c4v),
113+
);
114+
````
115+
116+
````
117+
[ Info: LBFGS: initializing with f = -5.047653728981e-01, ‖∇f‖ = 1.9060e-01
118+
[ Info: LBFGS: iter 1, Δt 1.41 s: f = -5.056459159937e-01, ‖∇f‖ = 1.3798e-01, α = 1.00e+00, m = 0, nfg = 1
119+
[ Info: LBFGS: iter 2, Δt 737.1 ms: f = -6.375541333037e-01, ‖∇f‖ = 1.7202e-01, α = 2.79e+01, m = 1, nfg = 5
120+
[ Info: LBFGS: iter 3, Δt 26.1 ms: f = -6.486427009452e-01, ‖∇f‖ = 1.3183e-01, α = 1.00e+00, m = 2, nfg = 1
121+
[ Info: LBFGS: iter 4, Δt 27.6 ms: f = -6.520903819553e-01, ‖∇f‖ = 1.2693e-01, α = 1.00e+00, m = 3, nfg = 1
122+
[ Info: LBFGS: iter 5, Δt 19.3 ms: f = -6.543775422131e-01, ‖∇f‖ = 8.4368e-02, α = 1.00e+00, m = 4, nfg = 1
123+
[ Info: LBFGS: iter 6, Δt 23.7 ms: f = -6.574414623345e-01, ‖∇f‖ = 9.2421e-02, α = 1.00e+00, m = 5, nfg = 1
124+
[ Info: LBFGS: iter 7, Δt 24.8 ms: f = -6.589599949721e-01, ‖∇f‖ = 4.1336e-02, α = 1.00e+00, m = 6, nfg = 1
125+
[ Info: LBFGS: iter 8, Δt 19.2 ms: f = -6.593158985369e-01, ‖∇f‖ = 1.6527e-02, α = 1.00e+00, m = 7, nfg = 1
126+
[ Info: LBFGS: iter 9, Δt 24.1 ms: f = -6.594942583549e-01, ‖∇f‖ = 1.3210e-02, α = 1.00e+00, m = 8, nfg = 1
127+
[ Info: LBFGS: iter 10, Δt 22.2 ms: f = -6.598272997895e-01, ‖∇f‖ = 1.2343e-02, α = 1.00e+00, m = 9, nfg = 1
128+
[ Info: LBFGS: iter 11, Δt 18.5 ms: f = -6.600089784768e-01, ‖∇f‖ = 8.5851e-03, α = 1.00e+00, m = 10, nfg = 1
129+
[ Info: LBFGS: iter 12, Δt 21.3 ms: f = -6.601648097143e-01, ‖∇f‖ = 3.1456e-03, α = 1.00e+00, m = 11, nfg = 1
130+
[ Info: LBFGS: iter 13, Δt 20.4 ms: f = -6.601883355292e-01, ‖∇f‖ = 2.2842e-03, α = 1.00e+00, m = 12, nfg = 1
131+
[ Info: LBFGS: iter 14, Δt 17.2 ms: f = -6.602036974772e-01, ‖∇f‖ = 2.8361e-03, α = 1.00e+00, m = 13, nfg = 1
132+
[ Info: LBFGS: iter 15, Δt 20.3 ms: f = -6.602112757039e-01, ‖∇f‖ = 2.0029e-03, α = 1.00e+00, m = 14, nfg = 1
133+
[ Info: LBFGS: iter 16, Δt 16.5 ms: f = -6.602199349095e-01, ‖∇f‖ = 1.2421e-03, α = 1.00e+00, m = 15, nfg = 1
134+
[ Info: LBFGS: iter 17, Δt 19.7 ms: f = -6.602252742030e-01, ‖∇f‖ = 7.3767e-04, α = 1.00e+00, m = 16, nfg = 1
135+
[ Info: LBFGS: iter 18, Δt 15.8 ms: f = -6.602292481916e-01, ‖∇f‖ = 6.4543e-04, α = 1.00e+00, m = 17, nfg = 1
136+
[ Info: LBFGS: iter 19, Δt 19.8 ms: f = -6.602308444814e-01, ‖∇f‖ = 3.7333e-04, α = 1.00e+00, m = 18, nfg = 1
137+
[ Info: LBFGS: iter 20, Δt 15.6 ms: f = -6.602310765298e-01, ‖∇f‖ = 2.5277e-04, α = 1.00e+00, m = 19, nfg = 1
138+
[ Info: LBFGS: converged after 21 iterations and time 1.42 m: f = -6.602310926956e-01, ‖∇f‖ = 2.8135e-05
139+
140+
````
141+
142+
We note that this energy is slightly higher than the one obtained from an
143+
[optimization using asymmetric CTMRG](@ref examples_heisenberg) with equivalent settings.
144+
Indeed, this is what one would expect since the $C_{4v}$ symmetry restricts the PEPS ansatz
145+
leading to fewer free parameters, i.e. an ansatz with reduced expressivity.
146+
Comparing against Juraj Hasik's data from $J_1\text{-}J_2$
147+
[PEPS simulations](https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.0/state_1s_A1_j20.0_D2_chi_opt48.dat),
148+
we find very good agreement:
149+
150+
````julia
151+
E_ref = -0.6602310934799577 # Juraj's energy at D=2, χ=16 with C4v symmetry
152+
@show (E - E_ref) / E_ref;
153+
````
154+
155+
````
156+
(E - E_ref) / E_ref = -1.1879467582794218e-9
157+
158+
````
159+
160+
As a consistency check, we can compute the vertical and horizontal correlation lengths,
161+
and should find that they are equal (up to the sparse eigensolver tolerance):
162+
163+
````julia
164+
ξ_h, ξ_v, = correlation_length(peps, env)
165+
@show ξ_h ξ_v;
166+
````
167+
168+
````
169+
ξ_h = [0.6625965820483917]
170+
ξ_v = [0.6625965820483924]
171+
172+
````
173+
174+
## QR-CTMRG
175+
176+
The conventional $C_{4v}$ CTMRG algorithm works by performing an eigendecomposition of
177+
the enlarged corner $A = V D V^\dagger$, taking the diagonal eigenvalue tensor as the new
178+
corner $C' = D$ and then renormalizing the edge using the isometry $V$. There exist modifications
179+
to this standard algorithm, most notably *QR-CTMRG* as presented in a recent publication by
180+
[Zhang, Yang and Corboz](@cite zhang_accelerating_2025). There the idea is to replace the
181+
eigendecomposition by a QR decomposition of a lower-rank approximation of the enlarged corner.
182+
While less accurate in some ways, the QR-CTMRG approach substantially accelerates contraction
183+
and optimization times, and also has vastly improved GPU performance. Notably, it is found
184+
that QR-CTMRG converges to the same fixed point as regular $C_{4v}$ CTMRG.
185+
186+
In PEPSKit terms, using QR-CTMRG just amounts to switching out the projector algorithm that is
187+
used by the [`C4vCTMRG`](@ref) algorithm to `projector_alg = :c4v_qr` (as opposed to `:c4v_eigh`).
188+
QR-CTMRG tends to need significantly more iterations to converge while still being much faster,
189+
hence we need to increase `maxiter`:
190+
191+
````julia
192+
env_qr₀, = leading_boundary(
193+
env_random_c4v, peps; alg = :c4v, projector_alg = :c4v_qr, maxiter = 500,
194+
);
195+
````
196+
197+
````
198+
[ Info: CTMRG init: obj = +5.600046917739e-03 err = 1.0000e+00
199+
┌ Warning: CTMRG cancel 500: obj = +5.924386039247e-01 err = 3.1337720153e-05 time = 0.23 sec
200+
└ @ PEPSKit ~/repos/PEPSKit.jl/src/algorithms/ctmrg/ctmrg.jl:153
201+
202+
````
203+
204+
To optimize using QR-CTMRG we proceed analogously by specifiying `projector_alg = :c4v_qr` and
205+
increasing the `maxiter` when setting the boundary algorithm parameters. We make sure to supply
206+
the `env_qr₀` initial environment because it does not use `DiagonalTensorMap`s as its corner
207+
type (only regular `eigh`-based $C_{4v}$ CTMRG produces diagonal corners):
208+
209+
````julia
210+
peps_qr, env_qr, E_qr, = fixedpoint(
211+
H, peps₀, env_qr₀;
212+
optimizer_alg = (; tol = 1.0e-4),
213+
boundary_alg = (; alg = :c4v, projector_alg = :c4v_qr, maxiter = 500),
214+
gradient_alg = (; alg = :linsolver)
215+
);
216+
@show (E_qr - E_ref) / E_ref;
217+
````
218+
219+
````
220+
[ Info: LBFGS: initializing with f = -5.047653728981e-01, ‖∇f‖ = 1.9060e-01
221+
[ Info: LBFGS: iter 1, Δt 689.5 ms: f = -5.056459386582e-01, ‖∇f‖ = 1.3798e-01, α = 1.00e+00, m = 0, nfg = 1
222+
[ Info: LBFGS: iter 2, Δt 505.7 ms: f = -6.375600534372e-01, ‖∇f‖ = 1.7220e-01, α = 2.79e+01, m = 1, nfg = 5
223+
[ Info: LBFGS: iter 3, Δt 25.8 ms: f = -6.486459057961e-01, ‖∇f‖ = 1.3187e-01, α = 1.00e+00, m = 2, nfg = 1
224+
[ Info: LBFGS: iter 4, Δt 34.8 ms: f = -6.520930551047e-01, ‖∇f‖ = 1.2680e-01, α = 1.00e+00, m = 3, nfg = 1
225+
[ Info: LBFGS: iter 5, Δt 21.3 ms: f = -6.543790524877e-01, ‖∇f‖ = 8.4444e-02, α = 1.00e+00, m = 4, nfg = 1
226+
[ Info: LBFGS: iter 6, Δt 22.8 ms: f = -6.576333383072e-01, ‖∇f‖ = 8.5748e-02, α = 1.00e+00, m = 5, nfg = 1
227+
[ Info: LBFGS: iter 7, Δt 34.4 ms: f = -6.589645527955e-01, ‖∇f‖ = 4.1392e-02, α = 1.00e+00, m = 6, nfg = 1
228+
[ Info: LBFGS: iter 8, Δt 26.4 ms: f = -6.593253867913e-01, ‖∇f‖ = 1.6340e-02, α = 1.00e+00, m = 7, nfg = 1
229+
[ Info: LBFGS: iter 9, Δt 175.3 ms: f = -6.595005575055e-01, ‖∇f‖ = 1.3085e-02, α = 1.00e+00, m = 8, nfg = 1
230+
[ Info: LBFGS: iter 10, Δt 20.1 ms: f = -6.598308924091e-01, ‖∇f‖ = 1.2318e-02, α = 1.00e+00, m = 9, nfg = 1
231+
[ Info: LBFGS: iter 11, Δt 16.3 ms: f = -6.600131668884e-01, ‖∇f‖ = 8.7100e-03, α = 1.00e+00, m = 10, nfg = 1
232+
[ Info: LBFGS: iter 12, Δt 18.9 ms: f = -6.601658264061e-01, ‖∇f‖ = 3.0216e-03, α = 1.00e+00, m = 11, nfg = 1
233+
[ Info: LBFGS: iter 13, Δt 14.0 ms: f = -6.601880521970e-01, ‖∇f‖ = 2.4696e-03, α = 1.00e+00, m = 12, nfg = 1
234+
[ Info: LBFGS: iter 14, Δt 19.6 ms: f = -6.602022384904e-01, ‖∇f‖ = 2.3369e-03, α = 1.00e+00, m = 13, nfg = 1
235+
[ Info: LBFGS: iter 15, Δt 14.6 ms: f = -6.602097424282e-01, ‖∇f‖ = 1.9888e-03, α = 1.00e+00, m = 14, nfg = 1
236+
[ Info: LBFGS: iter 16, Δt 19.8 ms: f = -6.602207442903e-01, ‖∇f‖ = 1.3522e-03, α = 1.00e+00, m = 15, nfg = 1
237+
[ Info: LBFGS: iter 17, Δt 161.0 ms: f = -6.602279508700e-01, ‖∇f‖ = 6.7499e-04, α = 1.00e+00, m = 16, nfg = 1
238+
[ Info: LBFGS: iter 18, Δt 99.1 ms: f = -6.602306993252e-01, ‖∇f‖ = 3.3521e-04, α = 1.00e+00, m = 17, nfg = 1
239+
[ Info: LBFGS: iter 19, Δt 185.0 ms: f = -6.602310628576e-01, ‖∇f‖ = 1.8734e-04, α = 1.00e+00, m = 18, nfg = 1
240+
[ Info: LBFGS: converged after 20 iterations and time 15.13 s: f = -6.602310908749e-01, ‖∇f‖ = 9.5413e-05
241+
(E_qr - E_ref) / E_ref = -3.945689570922226e-9
242+
243+
````
244+
245+
---
246+
247+
*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
248+

0 commit comments

Comments
 (0)