Conditioning and degeneracy of the environment-differentiation linear systems in CTM-based projected entangled-pair state optimization

Variational optimization of infinite projected entangled-pair states (iPEPS) contracted with the corner transfer matrix renormalization group (CTMRG) requires differentiating the contraction environment. In the stable fixed-point differentiation scheme of Francuz, Schuch and Vanhecke, this reduces to Sylvester-type linear systems for the derivative of the truncation projector and its adjoints, which are uniquely solvable only when the retained and discarded blocks of the corner spectrum are strictly separated and are in practice poorly conditioned, so far handled only by heuristic Krylov solvers with a corner-inverse right preconditioner. The same collaboration's later implicit-differentiation framework removes the gauge and degeneracy artifacts of differentiating truncated eigendecompositions and improves cost and scaling, but retains the spectral-gap requirement and the heuristic preconditioning, and reports gradient-solver convergence fluctuations along optimization trajectories that are difficult to control uniformly. The open problem is a systematic numerical analysis of these environment-differentiation linear systems: conditioning bounds in terms of the corner spectrum and bond dimensions, provably robust solver and preconditioner designs, and a characterization of, and remedy for, the exactly degenerate corner-spectrum regime in which a multiplet straddles the truncation threshold, including whether the variational energy gradient remains well-defined there.

Background

Projected entangled-pair states (PEPS) generalize matrix product states (their one-dimensional counterpart) to two-dimensional lattices: a rank-five tensor is placed on every site, carrying one physical index and four virtual indices contracted with the neighboring tensors, and the dimension $D$ of the virtual bonds (the bond dimension) controls the amount of entanglement the ansatz can represent. In the infinite, translationally invariant version (iPEPS), the variational ground-state energy of a two-dimensional Hamiltonian is minimized over the tensor components with gradient-based optimizers such as L-BFGS (a memory-efficient quasi-Newton method), and this is currently the state-of-the-art approach to strongly correlated lattice systems in two dimensions.

The energy itself cannot be computed exactly, because contracting the infinite two-dimensional network of the PEPS norm is intractable. The standard approximation is the corner transfer matrix renormalization group (CTMRG) algorithm of Nishino and Okunishi: an iterative procedure that represents the half-infinite environment by corner matrices $C$ and edge (half-row transfer) tensors $T$, and repeatedly absorbs one further row and column of the lattice into an enlarged corner matrix $M$, which is then diagonalized (for an ansatz with the full square-lattice point-group symmetry $C_{4v}$) or factorized by a singular value decomposition (SVD, the general case). Only the dominant $\chi$ eigenvalues are kept, with $\chi$ the environment bond dimension setting the contraction accuracy, and the eigenvectors of the kept block define the projector $P$ used to renormalize the environment.

Because the environment depends implicitly on the PEPS tensor, the energy gradient requires differentiating this iterative procedure. Reverse-mode automatic differentiation (backpropagation) through all CTMRG iterations is memory-hungry, and the derivatives of a truncated eigendecomposition or SVD are intrinsically delicate: the classical formulas for eigenvector and singular-vector derivatives collected by Giles contain terms like $1/(\lambda_i-\lambda_j)$ that diverge when eigenvalues or singular values coincide. Francuz, Schuch and Vanhecke made this pipeline stable and correct: they differentiate the CTM fixed point directly (fixed-point differentiation, i.e. differentiating the condition that the CTMRG iterate be stationary rather than unrolling the iterations), include the discarded part of the spectrum that earlier backpropagation implementations silently dropped as the dominant source of gradient error, and obtain, for the derivative $dP$ of the truncation projector and for the adjoint variables of the backward pass, linear systems of Sylvester type (a matrix equation of the form $AY+YB=D$). Two limitations of this scheme are the subject of the present problem. First, the Sylvester equation for $dP$ has a unique solution only if the spectra of the retained block $C$ and the discarded block $C_\perp$ of the enlarged corner are strictly separated, so in practice one must avoid splitting multiplets (groups of exactly degenerate eigenvalues) of the corner spectrum across the truncation threshold when converging the CTM; what happens, and what should be done, when the spectra are exactly interleaved is not analyzed. Second, the adjoint Sylvester systems are, in the authors' words, 'usually poorly conditioned', and are solved by a Krylov-type iterative method (stabilized biconjugate gradients; a Krylov method builds the solution from successive applications of the operator to the right-hand side) with a heuristic right preconditioner $C^{-1}$ (a preconditioner is an easily invertible approximation of the operator that accelerates the iteration); no systematic convergence or stability analysis is given.

The same collaboration subsequently reformulated the entire gradient computation (Burgelman et al., 2026). Instead of differentiating the CTMRG subroutines, they characterize the contraction environment by algebraic characteristic equations and differentiate only those (implicit differentiation). This removes the need to differentiate the truncation step at all: rotations of the basis inside the kept subspace are a gauge freedom of the environment, all physical observables being invariant under them, so the divergent in-space terms can be dropped, which eliminates the instabilities associated with degeneracies within the kept block. It also turns the nested sequence of Sylvester solves of fixed-point differentiation into a single linear problem with better asymptotic scaling, and it has been implemented in Julia (PEPSKit.jl together with ImplicitDifferentiationPEPS.jl) and generalized to non-symmetric CTMRG, non-trivial unit cells, and boundary matrix product state contraction (an alternative contraction scheme in which the environment is a matrix product state optimized variationally). Crucially for the present problem, this framework does not remove the two limitations above: the adjoint for genuine changes of the truncated subspace is, in their words, 'determined by the spectral gap between the smallest kept eigenvalue and the largest truncated eigenvalue' and is 'stable as long as we do not truncate through nearly degenerate eigenvalues'; the corner-inverse right preconditioner is retained heuristically; and their benchmarks report gradient-solver convergence that 'may fluctuate across the optimization trajectory, making it difficult to control uniformly across' the bond dimensions $D$ and $\chi$.

Adjacent work addresses conditioning and gauge questions at other levels of the pipeline: Zhang, Yang, Corboz, Haegeman and Tang (2026) precondition the outer optimization loop with the leading term of the metric tensor, attacking the ill-conditioned optimization landscape rather than the backward-pass linear systems; Tang, Vanderstraeten and Haegeman (2025) show that gradient-based optimizers exploit the gauge freedom of the PEPS tensors themselves, which degrades contraction accuracy and produces artificially low variational energies, and introduce a gauge-fixed optimization strategy. Neither touches the linear systems that differentiate the contraction environment.

Problem Statement

Make the gradient of a CTM-contracted iPEPS variational energy well-defined and stably computable at controlled cost for all spectral regimes of the corner transfer matrix. The objects under study are the linear systems that determine the derivative of the contraction environment with respect to the PEPS parameters, namely the Sylvester equation for the truncation-projector derivative $dP$ and the vectorized adjoint Sylvester systems of Francuz, Schuch and Vanhecke, the Jacobian of the characteristic equations of Burgelman et al., or an equivalent formulation; all of them share the same stability limits. The problem comprises three jointly supporting parts.

(1) Conditioning theory. Establish bounds on the condition numbers and on the iteration counts of Krylov-type solvers for these systems in terms of accessible quantities: the corner spectrum and its gaps, the PEPS bond dimension $D$, the environment bond dimension $\chi$, and the position along the optimization trajectory. The bounds should make explicit how they degrade as the gap between the smallest retained and the largest discarded corner eigenvalue closes.

(2) Robust solvers. Design solvers and preconditioners for these systems that improve on the heuristic corner-inverse right preconditioner, with guaranteed or uniformly demonstrated convergence across parameter regimes and along full optimization trajectories, at a per-gradient cost comparable to one forward CTMRG sweep.

(3) The exactly degenerate regime. Characterize the case in which the corner spectrum is exactly degenerate across the truncation threshold, i.e. a multiplet straddling the kept/discarded boundary or, more generally, interleaved retained and discarded spectra: determine whether the variational energy gradient with respect to the PEPS tensor is still well-defined (for example as a limit, or in a gauge- or subspace-invariant sense), and provide a stable numerical procedure that computes it, thereby removing the current operational constraint that one must avoid splitting multiplets when converging the CTM.

Answers should demonstrate the claimed behavior on recognized iPEPS benchmarks, for example Heisenberg and Fermi-Hubbard optimization runs, including instances in which exact degeneracies straddle the truncation threshold. The common objective of all three parts is that the environment gradient be well-defined and stably computable in every spectral regime of the corner matrix, so an answer that covers only the strictly gapped regime does not resolve the problem.

References

1. A. Francuz, N. Schuch, B. Vanhecke, Stable and efficient differentiation of tensor network algorithms, Phys. Rev. Research 7, 013237 (2025), DOI: 10.1103/PhysRevResearch.7.013237, arXiv:2311.11894, https://arxiv.org/abs/2311.11894
2. L. Burgelman, A. Francuz, P. Brehmer, L. Devos, J. Haegeman, F. Verstraete, B. Vanhecke, Implicit differentiation of tensor network algorithms, arXiv:2607.15030 (2026), https://arxiv.org/abs/2607.15030
3. P. Brehmer, L. Burgelman, Z.-Y. Yue, G. Fedorovich, J. Haegeman, L. Devos, PEPSKit.jl: A Julia package for projected entangled-pair state simulations, arXiv:2605.19960 (2026), https://arxiv.org/abs/2605.19960
4. X.-Y. Zhang, Q. Yang, P. Corboz, J. Haegeman, W. Tang, Accelerating two-dimensional tensor network optimization by preconditioning, Phys. Rev. B 113, 125111 (2026), DOI: 10.1103/h396-yc28, arXiv:2511.09546, https://arxiv.org/abs/2511.09546
5. W. Tang, L. Vanderstraeten, J. Haegeman, Gauging the variational optimization of projected entangled-pair states, arXiv:2508.10822 (2025), https://arxiv.org/abs/2508.10822
6. Y. Sugimoto, Forward-mode automatic differentiation for the tensor renormalization group and its relation to the impurity method, Phys. Rev. D 113, 094502 (2026), DOI: 10.1103/9wb1-rrcn, arXiv:2602.08987, https://arxiv.org/abs/2602.08987
7. Y. Zhang, P. Corboz, Fast two-dimensional tensor-network contraction via subspace iteration, arXiv:2607.15158 (2026), https://arxiv.org/abs/2607.15158
8. Y. Zhang, Q. Yang, P. Corboz, Accelerating two-dimensional tensor network contractions using QR decompositions, Phys. Rev. B 113, L201106 (2026), DOI: 10.1103/ptls-kr9z, arXiv:2505.00494, https://arxiv.org/abs/2505.00494
9. A. D. S. Richards, E. S. Sørensen, Ace-TN: GPU-Accelerated Corner-Transfer-Matrix Renormalization of Infinite Projected Entangled-Pair States, arXiv:2503.13900 (2025), https://arxiv.org/abs/2503.13900
10. J. Naumann, E. L. Weerda, J. Eisert, M. Rizzi, P. Schmoll, Variationally optimizing infinite projected entangled-pair states at large bond dimensions: A split corner transfer matrix renormalization group approach, Phys. Rev. B 111, 235116 (2025), DOI: 10.1103/PhysRevB.111.235116, arXiv:2502.10298, https://arxiv.org/abs/2502.10298
11. J. Naumann, J. Eisert, P. Schmoll, Variational optimization of projected entangled-pair states on the triangular lattice, Phys. Rev. B 113, 045117 (2026), DOI: 10.1103/g5gm-tzf8, arXiv:2510.04907, https://arxiv.org/abs/2510.04907
12. I. N. M. Le, R. Wiersema, C. B. Mendl, Hessian-vector products for tensor networks via recursive tangent-state propagation, arXiv:2604.20384 (2026), https://arxiv.org/abs/2604.20384
13. J. Naumann, E. L. Weerda, M. Rizzi, J. Eisert, P. Schmoll, An introduction to infinite projected entangled-pair state methods for variational ground state simulations using automatic differentiation, SciPost Phys. Lect. Notes 86 (2024), DOI: 10.21468/SciPostPhysLectNotes.86, arXiv:2308.12358, https://arxiv.org/abs/2308.12358
14. H. Chen, T. Neupert, J. Hasik, Simulating Fermionic Fractional Chern Insulators with Infinite Projected Entangled-Pair States, arXiv:2512.20697 (2025), https://arxiv.org/abs/2512.20697
15. T. Nishino, K. Okunishi, Corner Transfer Matrix Renormalization Group Method, J. Phys. Soc. Jpn. 65, 891 (1996), DOI: 10.1143/JPSJ.65.891, https://doi.org/10.1143/JPSJ.65.891
16. M. B. Giles, Collected Matrix Derivative Results for Forward and Reverse Mode Algorithmic Differentiation, in Advances in Automatic Differentiation, Lecture Notes in Computational Science and Engineering (Springer, Berlin Heidelberg, 2008), pp. 35-44, DOI: 10.1007/978-3-540-68942-3_4, https://doi.org/10.1007/978-3-540-68942-3_4