Issue with `SimultaneousCTMRG` convergence starting from random initial environment

[leburgel]

I ran into an issue when running leading_boundary on a $2 \times 2$ InfinitePartitionFunction using SimultaneousCTMRG (originally pointed out by @dartsushi).

Usually, when running leading_boundary at a fixed environment bond dimension $\chi$ using either trscheme = FixedSpaceTruncation() or trscheme = truncdim(χ), I would initialize a CTMRG environment with a total bond dimension $\chi$ first and then plug this into the leading_boundary call. However, it turns out that the default randomly initialized environments can give very bad convergence if you directly intialize them with a nontrivial virtual space. Initializing with a trivial virtual space and growing it solves this, but in this case you can't use FixedSpaceTruncation.

For example, for the 2D classical Ising model with a $2 x 2$ unit cell at criticality:

using TensorKit
using PEPSKit
using Random

using MPSKitModels: classical_ising

Random.seed!(1234) # same problem for many seeds...

# toggle symmetry, but same issue for both
S = Z2Irrep
# S = Trivial

χ = 20
tol = 1.0e-4
maxiter = 1000
verbosity = 2

# initialize
T = classical_ising(S)
O = T[1]
n = InfinitePartitionFunction([O O; O O])
Venv = S == Z2Irrep ? Z2Space(0 => χ / 2, 1 => χ / 2) : ℂ^χ
P = space(O, 2)
trscheme = FixedSpaceTruncation()

# simultaneous CTMRG: very sensitive to initial guess
simultaneous_alg = (;
    alg = :simultaneous,
    tol,
    verbosity,
    trscheme,
    maxiter,
)
env0 = CTMRGEnv(n, Venv)
env_sim, = leading_boundary(env0, n; simultaneous_alg...)

gives

[ Info: CTMRG init:     obj = -2.222242140629e+01 +1.949978774954e+01im err = 1.0000e+00
┌ Warning: CTMRG cancel 1000:   obj = +3.739239477627e+01 -5.522835410994e-01im err = 2.5066866523e-01  time = 11.77 sec
└ @ PEPSKit ~/ctmbench/lib/PEPSKit.jl/src/algorithms/ctmrg/ctmrg.jl:152

Initializing with a trivial virtual space and growing until truncdim(χ) to grow the environment first seems to avoid this,

env0 = CTMRGEnv(n, oneunit(P))
env_sim, = leading_boundary(env0, n; trscheme = truncdim(χ), simultaneous_alg...)

[ Info: CTMRG init:     obj = +3.397056274848e+01       err = 1.0000e+00
[ Info: CTMRG conv 534: obj = +4.121413131732e+01       err = 9.9855540314e-05  time = 6.24 sec

but this can give a space that is slightly different from the one we want,

space(env_sim.corners[1])

Rep[ℤ₂](0=>11, 1=>9) ← Rep[ℤ₂](0=>11, 1=>9)

For some reason, SequentialCTMRG seems to suffer much less from this problem:

sequential_alg = (;
    alg = :sequential,
    tol,
    verbosity,
    trscheme = FixedSpaceTruncation(),
    maxiter,
)
env0 = CTMRGEnv(n, Venv)
env_seq, = leading_boundary(env0, n; sequential_alg...)

[ Info: CTMRG conv 318: obj = +4.121402200669e+01 -7.312051507224e-05im err = 9.9978661272e-05  time = 4.70 sec
[ Info: CTMRG init:     obj = +3.397056274848e+01       err = 1.0000e+00

To circumvent the problem for SimultaneousCTMRG, sometimes initializing with ones instead of randn stabilizes things, but seems a bit too black magic to actually change the default.

A better solution that seems to really solve things is to initialize an environment with the proper virtualspace Venv by growing it from a trivial virtual space using a few dummy CTMRG iteration with trscheme = truncspace(Venv), and then starting from there,

dummy_alg = (;
    alg = :simultaneous,
    trscheme = truncspace(Venv),
    maxiter = 5,
)
env0 = CTMRGEnv(n, oneunit(P))
env1, = leading_boundary(env0, n; dummy_alg...)
env_sim_better, = leading_boundary(env1, n; simultaneous_alg...)

[ Info: CTMRG init:     obj = +3.397056274848e+01       err = 1.0000e+00
┌ Warning: CTMRG cancel 5:      obj = +4.095118941961e+01       err = 9.9825766951e-02  time = 0.02 sec
└ @ PEPSKit ~/ctmbench/lib/PEPSKit.jl/src/algorithms/ctmrg/ctmrg.jl:152
[ Info: CTMRG init:     obj = +4.095118941961e+01       err = 1.0000e+00
[ Info: CTMRG conv 190: obj = +4.121361334018e+01       err = 9.8857674692e-05  time = 0.82 sec

So long story short, I was wondering if it's sensible to just not use random initializations, but always initialize an environment of with a requested virtual space by growing it from a smaller one?

[dartsushi]

To circumvent the problem for SimultaneousCTMRG, sometimes initializing with ones instead of randn stabilizes things, but seems a bit too black magic to actually change the default.

CTMRG is practically computing 2N x 2N open-boundary partition funciton. The initial condition of the boundary tensors then represents the bounadry condition thereof. If we start from bond dimension 1 with ones this corresponds to a fixed boundary condition with a spin in a specific direction. Starting from radn with bond dimension chi, on the other hand, means that your boundary spins are highly entangled in a random way. This seems to be very problematic to me.

[pbrehmer]

Thanks for documenting this! I feel that I encountered similar issues with contracting PEPS and randomly initialized environments before but I've never tried to pin them down. Also not sure where the difference between :simultaneous and :sequential comes from and whether that is a systematic thing.

I'm wondering whether growing the initial environment from a smaller one generically is the better approach. It seems that using randn initialized environment with fixed virtual spaces does work reasonably well in many cases. It might still be that in these cases growing the initial environment improves convergence.

In any case I think we should have the option to (automatically) initialize from a smaller environment, so I'm linking the roadmap where we actually discussed that it would be useful to have more environment initialization techniques.

[Yue-Zhengyuan]

Using random initializations is fine if you first converge for an environment with small χ, and then grow it by relaxing the truncation scheme (using dynamic ones like truncerr, truncbelow together with truncdim(χ)). For sequential CTMRG it is also likely to fail to converge if you directly set χ to a large value that you want to obtain in the end. The example on fermi-Hubbard simple update demonstrates my usual workflow (here for simplicity I only used truncdim).

PEPSKit.jl/examples/hubbard_su/main.jl

Lines 75 to 93 in 1ac5857

	md"""
	## Computing the ground-state energy
	In order to compute the energy expectation value with evolved PEPS, we need to converge a
	CTMRG environment on it. We first converge an environment with a small enviroment dimension
	and then use that to initialize another run with bigger environment dimension. We'll use
	`trscheme=truncdim(χ)` for that such that the dimension is increased during the second CTMRG
	run:
	"""
	χenv₀, χenv = 6, 16
	env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2)
	normalize!.(peps.A, Inf)
	env = CTMRGEnv(rand, Float64, peps, env_space)
	for χ in [χenv₀, χenv]
	global env, = leading_boundary(
	env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trscheme = truncdim(χ)
	)
	end

[Yue-Zhengyuan]

BTW, besides of being completely unable to converge, do you sometimes observe that the obj and err jumps back and forth between two (or a bit more) values? like the following when I try to measure a Gutzwiller-projected Gaussian fermionic PEPS (with SequentialCTMRG):

[Info] CTMRG  34:	obj = +2.142606851162e+00	err = 2.5476124790e-02	time = 2.08 sec
[Info] CTMRG  35:	obj = +2.142606851663e+00	err = 2.9776378201e-02	time = 2.31 sec
[Info] CTMRG  36:	obj = +2.142606851684e+00	err = 1.5221951224e-02	time = 1.95 sec
[Info] CTMRG  37:	obj = +2.142606851162e+00	err = 2.5476125480e-02	time = 1.86 sec
[Info] CTMRG  38:	obj = +2.142606851663e+00	err = 2.9776376717e-02	time = 1.92 sec
[Info] CTMRG  39:	obj = +2.142606851684e+00	err = 1.5221951225e-02	time = 1.79 sec
[Info] CTMRG  40:	obj = +2.142606851162e+00	err = 2.5476127930e-02	time = 2.16 sec
[Info] CTMRG  41:	obj = +2.142606851663e+00	err = 2.9776377441e-02	time = 2.00 sec
[Info] CTMRG  42:	obj = +2.142606851684e+00	err = 1.5221951020e-02	time = 2.00 sec
[Info] CTMRG  43:	obj = +2.142606851162e+00	err = 2.5476123501e-02	time = 2.21 sec
[Info] CTMRG  44:	obj = +2.142606851663e+00	err = 2.9776375638e-02	time = 2.09 sec
[Info] CTMRG  45:	obj = +2.142606851684e+00	err = 1.5221951187e-02	time = 1.96 sec
[Info] CTMRG  46:	obj = +2.142606851162e+00	err = 2.5476126710e-02	time = 1.76 sec
[Info] CTMRG  47:	obj = +2.142606851663e+00	err = 2.9776376392e-02	time = 1.78 sec
[Info] CTMRG  48:	obj = +2.142606851684e+00	err = 1.5221950872e-02	time = 1.83 sec
[Info] CTMRG  49:	obj = +2.142606851162e+00	err = 2.5476123580e-02	time = 2.08 sec
[Warn] CTMRG cancel 50:	obj = +2.142606851663e+00	err = 2.9776376416e-02	time = 1.85 min

Here the period is 3.

[pbrehmer]

I have definitely seen this before (I think with both :sequential and :simultaneous) and I've always wondered how these cycles work and what the underlying reason is that CTMRG can get stuck in such a loop.

[lkdvos]

CTMRG is practically computing 2N x 2N open-boundary partition funciton. The initial condition of the boundary tensors then represents the bounadry condition thereof. If we start from bond dimension 1 with ones this corresponds to a fixed boundary condition with a spin in a specific direction. Starting from radn with bond dimension chi, on the other hand, means that your boundary spins are highly entangled in a random way. This seems to be very problematic to me.

While I somewhat agree with this intuitive picture, keep in mind that if you start from an MPS with random entries, the spins aren't highly entangled at all, in fact the entanglement basically becomes zero with increasing system size. This is different from starting from a random boundary (dense) state, where you indeed would expect a volume law of entanglement.

---

Something that might also be relevant, especially when dealing with symmetries is that for iMPS simulations there actually is a very similar kind of pattern, especially when the initialization happens with the wrong distribution of sectors. It seems that in these cases, the local updates are never large enough to redistribute the weight in the different sectors, and therefore the entire algorithm becomes stuck in a local basin. In that sense I would also expect that just starting from a product state in the right phase and gradually growing it helps since the distribution of weights between the different sectors is at least based on (small patches of) the actual PEPS.
It might just be more important that the state is low bond dimension, such that it doesn't prioritize one of the sectors, rather than whether or not the entanglement is high.

---

I have definitely also noticed these cycles. It might be useful to have a look at the corner entanglement spectra in these cases and see how they change. Sometimes it is really clear that this comes from one value being truncated or not, or multiplets that are being cut, or ...

[leburgel]

Okay, setting aside the strange cycles and how to diagnose them, I think the conclusion is that some more robust initialization routines would be useful to many people. The roadmap mostly discusses the use of environments during a full optimization, specifically between iterations. However, I think a good place to start is to just provide some initialization tools in the first place, after which we can think how to integrate this into the optimization.

Specifically, we could provide functions like initialize_environment and initialize_peps (or initialize_fixedpoint which combines both?) that take some specification of what we want to end up with (say virtual spaces) and how to get there (say a contraction algorithm, or at least parts that specify one to the degree that we need). This would solve the original problem brought up here, and if properly demonstrated could save people a lot of headaches figuring out why their contractions aren't converging.

[lkdvos]

I like the idea of designing an initialization interface. I would also like this in MPSKit, do you think there is a way we can share similar call signatures/functions?

[leburgel]

I'm sure we can make the signatures match, should be very similar. I can make a PR here which wraps the fix outlined in the issue here into an initialization interface. Then we have an idea how this would look like, and we can discuss what would have to change to make it usable elsewhere and where to best implement it.