Pith. sign in

REVIEW 3 major objections 3 minor 1 cited by

PEPS energy gradients become one linear solve instead of nested backpropagation.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · deepseek-v4-flash

2026-08-02 00:21 UTC pith:UOYZKEKW

load-bearing objection Honest, useful framework with a real flaw: the boundary-MPS characteristic equations look scale-singular, so the advertised linear solve is not well-posed as written. the 3 major comments →

arxiv 2607.15030 v1 pith:UOYZKEKW submitted 2026-07-16 quant-ph cond-mat.str-elphysics.comp-ph

Implicit differentiation of tensor network algorithms

Lander Burgelman , Anna Francuz , Paul Brehmer , Lukas Devos , Jutho Haegeman , Frank Verstraete , Bram Vanhecke This is my paper
classification quant-ph cond-mat.str-elphysics.comp-ph
keywords implicit differentiationtensor networksPEPSgradient optimizationCTMRGboundary MPSautomatic differentiationcharacteristic equation
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper aims to make gradient-based optimization of projected entangled-pair states (PEPS) cheaper and more stable. It argues that the energy gradient does not need to be computed by differentiating the iterative contraction algorithm (CTMRG or boundary MPS) subroutine by subroutine. Instead, the converged contraction environment can be characterized as the root of a purely algebraic equation F(y,p)=0, and the gradient follows from solving a single linear system. This un-nests the gradient problem, removes instabilities caused by differentiating eigensolvers and SVDs, and makes implementation easier because only tensor contractions need to be differentiated.

Core claim

On the paper's own terms, the central claim is that the derivative de/dp = xbar ∂p x + pbar can be evaluated without ever differentiating the contraction algorithm's inner subroutines. One writes the converged environment as the unique root y*(p) of a characteristic equation F(y,p)=0 chosen so that it is purely algebraic and so that its variables parametrize only the physical variations of the environment — isometry variations along the fixed-point subspace are removed by construction. Differentiating F at the root gives the adjoint equation Fbar ∂yF = ybar, and the desired vector–Jacobian product is −Fbar ∂pF. The paper works out explicit characteristic equations for the three most common c

What carries the argument

The key object is the characteristic equation F(y,p)=0 — a set of tensor-contraction equations that implicitly define the converged contraction environment, augmented with the defining relations of the truncated eigenspace or SVD and parametrized by variables y that live only in the null spaces of the fixed-point isometries. This parametrization removes the gauge-induced zero modes and divergent terms: variations along the kept subspace (skew-Hermitian ω terms) are dropped, while the corner/edge tensors are treated as generic complex Hermitian matrices in the reverse pass. The adjoint linear system Fbar ∂yF = ybar is then solved once with a Krylov solver; evaluating ∂yF and ∂pF requires only

Load-bearing premise

The load-bearing premise is that the algebraic equations F(y,p)=0 have a unique, isolated root at the contraction fixed point and that the Jacobian ∂yF is invertible there for every parameter p encountered during optimization — if a residual gauge or normalization freedom leaves a zero mode, the single linear solve fails and the gradient is undefined.

What would settle it

Take a PEPS parameter p where the CTMRG truncation step has exactly degenerate leading eigenvalues or singular values and check whether the GMRES solution of Eq. (17) converges and matches a finite-difference gradient. A second check: scan the optimization trajectory and monitor the smallest singular value of ∂yF; if it touches zero away from phase-transition-like gaps, the core well-posedness assumption is violated.

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

  • Gradient evaluation for PEPS becomes a single un-nested linear solve, asymptotically subleading to the contraction itself in the tested settings.
  • No subroutine of the contraction algorithm (eigensolvers, SVDs, QR) needs to be made differentiable; only tensor contractions are differentiated.
  • Numerical instabilities from degenerate eigenvalues or singular values in the pullback are removed by the gauge-aware parametrization, since the divergent in-subspace rotations are not part of the variables.
  • The same characteristic-equation approach applies on top of any existing CTMRG or boundary MPS implementation, including QR-based and pulling-through variants, by mapping the forward environment to the root y*.
  • For poorly converged environments, gradient accuracy matches black-box reverse-mode AD when the linear solver tolerance is adjusted accordingly.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same 'characteristic equation plus one linear solve' recipe likely transfers to other tensor network algorithms beyond PEPS, such as time evolution or excitation ansätze, whenever the fixed point can be characterized algebraically.
  • The alternative boundary-MPS decomposition into three smaller linear problems suggests a general trade-off: one large linear solve can be exchanged for smaller parallel ones, which may matter at very large environment bond dimensions.
  • Because the conditioning of the linear system is set by the spectral gap between kept and truncated eigenvalues, near a phase transition where that gap closes the gradient may become ill-conditioned — but this shows up visibly as solver non-convergence rather than hidden in AD pullbacks.
  • A direct test of the zero-mode concern would be to intentionally construct a degenerate truncation subspace and check whether GMRES still converges for Eq. (17).

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

3 major / 3 minor

Summary. The paper proposes an implicit differentiation framework for computing energy gradients in variational PEPS optimization. Instead of backpropagating through the iterative CTMRG or boundary-MPS contraction routines, the authors introduce an algebraic characteristic equation F(y,p)=0 whose root characterizes the converged contraction environment, and evaluate the gradient by solving a single linear system \bar F ∂_yF = \bar y. Explicit characteristic equations are given for C4v-symmetric CTMRG, boundary MPS, and general asymmetric CTMRG, and benchmarks compare the approach with fixed-point differentiation. The central claim is that this un-nests the gradient computation and removes gauge-induced instabilities from the differentiation of eigensolvers and SVDs.

Significance. If the framework is correct, it is a valuable methodological contribution: it replaces nested differentiation of iterative subroutines with one linear solve, simplifies the AD implementation, and ships with open-source code and benchmarks. The C4v and general CTMRG derivations are carefully constructed, and the numerical speedups for those cases are plausible. However, the boundary-MPS formulation, which is one of the three main examples, has a load-bearing singularity issue that is not addressed in the manuscript. The claim that the characteristic equations are well-posed for implicit differentiation is therefore not established for that section.

major comments (3)
  1. [§IV B, Eqs. (50)–(57)] The characteristic system is scale-invariant. For any α, the transformation G_tL → α G_tL, G_tR → α^{-1}G_tR leaves λ_t in Eq. (50) unchanged and leaves each of Eqs. (51)–(55) satisfied at the root; the same holds for the bottom pair (G_bL, G_bR) and for the mixed pair (G_L, G_R) with λ in Eq. (58). Thus the 12-variable system has at least three independent null directions, ∂_yF is singular at the root, and Eq. (17) is not well-posed as written. The normalization statements around Eq. (59) are forward fixed-point conventions, not constraints included in F. The authors must either add explicit gauge-fixing/normalization equations to F, or prove that \bar y lies in the range of ∂_yF and that any solution of Eq. (17) yields the same physical gradient. As it stands, the central boundary-MPS claim is not defined.
  2. [Appendix B and §II C, Eq. (17)] Appendix B analyzes only a scalar eigenvalue problem and explicitly shows that without a normalization/phase-fixing condition the tangent is not unique. The manuscript asserts that the analogous zero-eigenvalue elimination works for the full 7-, 12-, and 20-variable systems, but no corresponding argument is given for the boundary-MPS scale modes. This is not merely a missing proof: the equations as displayed are singular, so the application of the implicit function theorem and the linear solve in Eq. (17) requires additional hypotheses that are not stated.
  3. [§VI, Fig. 1 and Figs. 4–5] The gradient accuracy benchmark in Fig. 1 compares implicit differentiation with finite differences only for C4v CTMRG and asymmetric CTMRG. The boundary-MPS section is validated only by timing comparisons in Figs. 4–5. Given the singular Jacobian identified above, a numerical correctness test for the boundary-MPS gradient is necessary to support the claim that the approach computes the correct gradient in that case.
minor comments (3)
  1. [§III B, Eq. (28)] The use of C^{*-1} as a preconditioner should mention that C* is assumed invertible; if the corner tensor can have zero or small eigenvalues, the conditioning of Eq. (17) needs discussion.
  2. [§IV A, text after Eq. (42)] The notation λ_t ≡ λ_t^L = λ_t^R = λ_t^AC is introduced informally; it would help to define the normalization convention exactly once, since normalization plays a role in the boundary-MPS zero modes.
  3. [Algorithm 2, line 6] The notation "˘y←(˘x ∂_y x,0)" is slightly confusing because x is a function of y; rewriting as the adjoint of the parametrization map y→x would improve clarity.

Circularity Check

0 steps flagged

No significant circularity: the gradient is obtained by standard implicit differentiation of self-constructed characteristic equations, with independent benchmark validation; the same-group citation is not load-bearing.

full rationale

The central derivation is a textbook implicit-function-theorem argument: given F(y,p)=0 characterizing the contraction fixed point, differentiating yields ∂yF ∂p y = −∂p F and the gradient is obtained by solving ˘F ∂yF = ˘y. The paper constructs F explicitly from the algebraic fixed-point and optimality conditions of CTMRG and boundary MPS; the λ scalars are defined as inner products (e.g. Eq. (50)) and are therefore differentiable functions of y, not fitted parameters. The root y* is explicitly the fixed-point environment; no data are fitted and no external prediction is derived from measurements. The benchmarks compare the implicit gradient to black-box reverse-mode AD and finite differences, providing independent confirmation. Ref. [29] is a same-group citation; it is cited for the prior stable-derivative construction that the C4v equations reproduce and as a source of the known eigenvalue-gap instabilities, but the F equations here are independently formulated and the stability claim is tested numerically. The boundary-MPS scaling zero-mode issue raised in the skeptic note concerns whether ∂yF is invertible — a mathematical well-posedness condition — not whether the derivation is circular; if valid it would make the gradient ambiguous, not tautological. Overall, no equation reduces to its own input by construction.

Axiom & Free-Parameter Ledger

0 free parameters · 4 axioms · 0 invented entities

No free parameters are fitted to data; the λ values in the characteristic equations are inner-product normalizers that eliminate zero eigenvalues, not tunable inputs. The central claim rests on standard linear algebra/IFT and on domain-specific assumptions about gauge symmetry and fixed-point isolation, listed above. The auxiliary variables (u, l, r, modified corners/edges) are internal parameterizations, not new physical entities.

axioms (4)
  • domain assumption Implicit function theorem applies: ∂yF is invertible at the root y* for the constructed 7-, 12-, and 20-variable systems.
    Invoked in Section II C to convert F(y,p)=0 into the VJP ˘y∂py = -˘F ∂pF; no proof of invertibility is given for the coupled systems, only for the scalar eigenvalue case in Appendix B.
  • domain assumption The energy density e(x,p) is invariant under in-space gauge rotations of the isometries (Q rotations), so variations along U* and A* can be dropped.
    Sections III C and IV C use this to justify the null-space parameterization (Eqs. 25, 47-48, 71-72); it holds for CTMRG and VUMPS environments but is not automatic for all possible contraction schemes.
  • domain assumption The contraction algorithms converge to an isolated fixed point with a gap between kept and truncated subspace, so the truncated environment is differentiable along the kept/null decomposition.
    Needed for the root y* to exist and for the gradient to be finite; the text notes instability appears when truncating through nearly degenerate eigenvalues (Section III C).
  • domain assumption For boundary MPS with reflection symmetry, the bottom eigenvector is the Hermitian conjugate of the top; without it, the 12-variable system is needed.
    Section IV A/B; if the transfer matrix is non-Hermitian and the bottom MPS is not captured by the top equations, the reduced 5-equation set would be incomplete.

pith-pipeline@v1.3.0-alltime-deepseek · 34862 in / 12904 out tokens · 133488 ms · 2026-08-02T00:21:25.598530+00:00 · methodology

0 comments
read the original abstract

The current leading approach to the variational optimization of projected entangled-pair states (PEPS) is based on automatic differentiation, which allows for a convenient evaluation of the energy gradient with respect to the local variational degrees of freedom. However, evaluating the energy gradient not only remains a major computational bottleneck of the optimization procedure, but also suffers from frequent numerical instabilities. In this work, we adopt recent advances in implicit differentiation techniques to address these challenges in PEPS optimization. By reformulating the core step of the gradient computation in terms of a single characteristic equation for the contraction environment, we reduce the cost of the gradient computation and improve its scaling with the problem size. By choosing a suitable parametrization of this characteristic equation based on the intrinsic symmetries of the contraction environment, we can directly remove instabilities from the global gradient computation that would otherwise arise from the derivatives of subroutines of the contraction algorithm. Finally, we demonstrate how this approach drastically simplifies the practical implementation of stable gradient-based PEPS optimization.

Figures

Figures reproduced from arXiv: 2607.15030 by Anna Francuz, Bram Vanhecke, Frank Verstraete, Jutho Haegeman, Lander Burgelman, Lukas Devos, Paul Brehmer.

Figure 1
Figure 1. Figure 1: FIG. 1. Relative distance ∆ of the directional derivative of the [PITH_FULL_IMAGE:figures/full_fig_p017_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: FIG. 2. Speedup [PITH_FULL_IMAGE:figures/full_fig_p018_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: FIG. 3. Gradient linear problem convergence times [PITH_FULL_IMAGE:figures/full_fig_p018_3.png] view at source ↗
Figure 6
Figure 6. Figure 6: FIG. 6. Ratio of the fixed-point and implicit gradient linear [PITH_FULL_IMAGE:figures/full_fig_p019_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: FIG. 7. Fixed-point and implicit gradient linear-solving [PITH_FULL_IMAGE:figures/full_fig_p019_7.png] view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

  1. On the origin of finite entanglement scaling

    quant-ph 2026-07 conditional novelty 7.0

    Finite-bond-dimension MPS approximations of two decoupled critical Ising chains induce entangling inter-chain perturbations beyond the most relevant CFT operator.

Reference graph

Works this paper leans on

62 extracted references · 7 linked inside Pith · cited by 1 Pith paper

  1. [1]

    A real-valued cost functioneof the complex parametersp is given by a mape:C n →R

    Automatic differentiation in a nutshell Consider a set of variational parametersp, which we can interpret as ann-component complex vectorp∈C n. A real-valued cost functioneof the complex parametersp is given by a mape:C n →R. We can take derivatives of such a cost function with respect to the complex variables pusing the formalism of Wirtinger derivatives...

  2. [2]

    the order in which the functions were applied,

  3. [3]

    at which value of the parameters the functions were evaluated,

  4. [4]

    22 Points 1

    the VJP action ˘p= (∂pg)† ˘gof every function used in the composition. 22 Points 1. and 2. can be easily handled automatically by any off-the-shelf AD engine. The VJP action in point

  5. [5]

    Having introduced the notions of JVP, VJP and reverse-mode accumulation, we can take a step back and phrase things a bit more informally

    needs to be manually specified for any function in the sequence that does not have a standard implementation. Having introduced the notions of JVP, VJP and reverse-mode accumulation, we can take a step back and phrase things a bit more informally. In particular, we will drop all complex conjugates, and generalize our no- tation from simple vectors in Eucl...

  6. [6]

    symmetric

    V ariational PEPS optimization using AD Returning to the problem of variational PEPS opti- mization of Section II A, we can now discuss how to eval- uate the gradient of the energy density Eq. (6) using reverse-mode accumulation. The most straightforward way to evaluate the VJP ˘x ∂pxin Eq. (6) is to examine how the contraction environment is computed in ...

  7. [7]

    Theeigh-basedC 4v CTMRG algorithm as de- scribed in Section III A [14]

  8. [8]

    The QR-basedC 4v CTMRG algorithm as intro- duced in Ref. 38

  9. [9]

    Crucially, all of these forward schemes converge to the same fixed point, meaning they can use the same char- acteristic equation in the backwards pass

    A Krylov-based contraction approach, known as thepulling-throughorfixed-point CTMRGalgo- rithm [17, 51], which is equivalent to the VUMPS algorithm in the presence of aC 4v symmetry. Crucially, all of these forward schemes converge to the same fixed point, meaning they can use the same char- acteristic equation in the backwards pass. To conclude this Appe...

  10. [10]

    Or´ us, Tensor networks for complex quantum systems, Nat

    R. Or´ us, Tensor networks for complex quantum systems, Nat. Rev. Phys.1, 538–550 (2019)

    2019

  11. [11]

    M. C. Ba˜ nuls, Tensor Network Algorithms: A Route Map, Annu. Rev. Condens. Matter Phys.14, 173–191 (2023)

    2023

  12. [12]

    J. I. Cirac, D. P´ erez-Garc ´ ıa, N. Schuch, and F. Ver- straete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys.93, 045003 (2021)

    2021

  13. [13]

    S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett.69, 2863 (1992)

    1992

  14. [14]

    S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B48, 10345 (1993)

    1993

  15. [15]

    Verstraete, D

    F. Verstraete, D. Porras, and J. I. Cirac, Density Matrix Renormalization Group and Periodic Boundary Condi- tions: A Quantum Information Perspective, Phys. Rev. Lett.93, 227205 (2004)

    2004

  16. [16]

    Schollw¨ ock, The density-matrix renormalization group in the age of matrix product states, Ann

    U. Schollw¨ ock, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. (N.Y.) 326, 96 (2011)

    2011

  17. [17]

    ¨Ostlund and S

    S. ¨Ostlund and S. Rommer, Thermodynamic Limit of Density Matrix Renormalization, Phys. Rev. Lett.75, 3537 (1995)

    1995

  18. [18]

    Verstraete and J

    F. Verstraete and J. I. Cirac, Renormalization algorithms for Quantum-Many Body Systems in two and higher di- mensions (2004), arXiv:cond-mat/0407066

    Pith/arXiv arXiv 2004

  19. [20]

    H. C. Jiang, Z. Y. Weng, and T. Xiang, Accurate De- termination of Tensor Network State of Quantum Lat- tice Models in Two Dimensions, Phys. Rev. Lett.101, 090603 (2008)

    2008

  20. [21]

    Corboz, Variational optimization with infinite pro- jected entangled-pair states, Phys

    P. Corboz, Variational optimization with infinite pro- jected entangled-pair states, Phys. Rev. B94, 035133 28 (2016)

    2016

  21. [22]

    Vanderstraeten, J

    L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Ver- straete, Gradient methods for variational optimization of projected entangled-pair states, Phys. Rev. B94, 155123 (2016)

    2016

  22. [23]

    Nishino and K

    T. Nishino and K. Okunishi, Corner Transfer Matrix Renormalization Group Method, J. Phys. Soc. Jpn.65, 891 (1996)

    1996

  23. [24]

    Or´ us and G

    R. Or´ us and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Cor- ner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009)

    2009

  24. [25]

    Jordan, R

    J. Jordan, R. Or´ us, G. Vidal, F. Verstraete, and J. I. Cirac, Classical Simulation of Infinite-Size Quantum Lat- tice Systems in Two Spatial Dimensions, Phys. Rev. Lett. 101, 250602 (2008)

    2008

  25. [26]

    Haegeman and F

    J. Haegeman and F. Verstraete, Diagonalizing Transfer Matrices and Matrix Product Operators: A Medley of Exact and Computational Methods, Annu. Rev. Con- dens. Matter Phys.8, 355 (2017)

    2017

  26. [27]

    Zauner-Stauber, L

    V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Variational optimization algorithms for uniform matrix product states, Phys. Rev. B97, 045145 (2018)

    2018

  27. [28]

    Liao, J.-G

    H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Differen- tiable Programming Tensor Networks, Phys. Rev. X9, 031041 (2019)

    2019

  28. [29]

    Hasik, D

    J. Hasik, D. Poilblanc, and F. Becca, Investigation of the N´ eel phase of the frustrated Heisenberg antiferromag- net by differentiable symmetric tensor networks, SciPost Phys.10, 012 (2021)

    2021

  29. [30]

    Ponsioen, F

    B. Ponsioen, F. F. Assaad, and P. Corboz, Automatic differentiation applied to excitations with projected en- tangled pair states, SciPost Phys.12, 006 (2022)

    2022

  30. [31]

    S. Niu, J. Hasik, J.-Y. Chen, and D. Poilblanc, Chiral spin liquids on the kagome lattice with projected entan- gled simplex states, Phys. Rev. B106, 245119 (2022)

    2022

  31. [32]

    Zhang, S

    X.-Y. Zhang, S. Liang, H.-J. Liao, W. Li, and L. Wang, Differentiable programming tensor networks for Kitaev magnets, Phys. Rev. B108, 085103 (2023)

    2023

  32. [33]

    Y. He, R. Rausch, M. Peschke, C. Karrasch, P. Corboz, N. Bultinck, and S. A. Parameswaran, Itinerant mag- netism in the triangular-lattice Hubbard model at half doping: Application to twisted transition metal dichalco- genides, Phys. Rev. B113, L041107 (2026)

    2026

  33. [34]

    J. Hu, R. Chi, Y. Guo, B. Normand, H.-J. Liao, and T. Xiang, Dynamical Spectral Function of the Kagome Quantum Spin Liquid (2025), arXiv:2512.18831

    arXiv 2025

  34. [35]

    H. Chen, T. Neupert, and J. Hasik, Simulating fermionic fractional Chern insulators with infinite pro- jected entangled-pair states (2026), arXiv:2512.20697

    Pith/arXiv arXiv 2026

  35. [36]

    Christianson, Reverse accumulation and attractive fixed points, Optim

    B. Christianson, Reverse accumulation and attractive fixed points, Optim. Methods Softw.3, 311 (1994)

    1994

  36. [37]

    Naumann, E

    J. Naumann, E. L. Weerda, M. Rizzi, J. Eisert, and P. Schmoll, An introduction to infinite projected entangled-pair state methods for variational ground state simulations using automatic differentiation, SciPost Phys. Lect. Notes86(2024)

    2024

  37. [38]

    Francuz, N

    A. Francuz, N. Schuch, and B. Vanhecke, Stable and effi- cient differentiation of tensor network algorithms, Phys. Rev. Res.7, 013237 (2025)

    2025

  38. [39]

    A. V. Fiacco, Sensitivity analysis for nonlinear program- ming using penalty methods, Math. Program.10, 287 (1976)

    1976

  39. [40]

    Blondel, Q

    M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-L´ opez, F. Pedregosa, and J.-P. Vert, Ef- ficient and Modular Implicit Differentiation (2022), arXiv:2105.15183

    Pith/arXiv arXiv 2022

  40. [41]

    Nocedal and S

    J. Nocedal and S. J. Wright,Numerical Optimization, Springer Series in Operations Research and Financial En- gineering (Springer New York, 2006)

    2006

  41. [42]

    Here we have chosen to denote adjoints as ˘xas opposed to the more conventional bar-notation ¯x, to avoid confusion with our notation for complex conjugation of tensors

  42. [43]

    M. B. Giles, Collected Matrix Derivative Results for For- ward and Reverse Mode Algorithmic Differentiation, in Advances in Automatic Differentiation(Springer Berlin Heidelberg, 2008) pp. 35–44

    2008

  43. [44]

    Vanderstraeten, L

    L. Vanderstraeten, L. Burgelman, B. Ponsioen, M. Van Damme, B. Vanhecke, P. Corboz, J. Haegeman, and F. Verstraete, Variational methods for contracting projected entangled-pair states, Phys. Rev. B105, 195140 (2022)

    2022

  44. [45]

    Nietner, B

    A. Nietner, B. Vanhecke, F. Verstraete, J. Eisert, and L. Vanderstraeten, Efficient variational contraction of two-dimensional tensor networks with a non-trivial unit cell, Quantum4, 328 (2020)

    2020

  45. [46]

    Hauru, M

    M. Hauru, M. Van Damme, and J. Haegeman, Rieman- nian optimization of isometric tensor networks, SciPost Phys.10, 040 (2021)

    2021

  46. [47]

    Zhang, Q

    Y. Zhang, Q. Yang, and P. Corboz, Accelerating two- dimensional tensor network contractions using QR- decompositions (2025), arXiv:2505.00494

    Pith/arXiv arXiv 2025

  47. [48]

    Hodgkiss, L

    L. Hodgkiss, L. Lootens, A. Ueda, B. Vanhecke, and F. Verstraete, On the origin of finite entanglement scal- ing, in preparation (2026)

    2026

  48. [49]

    Vanhecke, M

    B. Vanhecke, M. Van Damme, J. Haegeman, L. Vander- straeten, and F. Verstraete, Tangent-space methods for truncating uniform MPS, SciPost Phys. Core4(2021)

    2021

  49. [50]

    Vanderstraeten, J

    L. Vanderstraeten, J. Haegeman, and F. Verstraete, Tangent-space methods for uniform matrix product states, SciPost Phys. Lect. Notes7(2019)

    2019

  50. [51]

    Corboz, J

    P. Corboz, J. Jordan, and G. Vidal, Simulation of fermionic lattice models in two dimensions with projected entangled-pair states: Next-nearest neighbor hamiltoni- ans, Phys. Rev. B82, 245119 (2010)

    2010

  51. [52]

    Brehmer, L

    P. Brehmer, L. Burgelman, Z.-Y. Yue, G. Fedorovich, J. Haegeman, and L. Devos, PEPSKit.jl: A Julia pack- age for projected entangled-pair state simulations (2026), arXiv:2605.19960

    Pith/arXiv arXiv 2026

  52. [53]

    Devos, M

    L. Devos, M. Van Damme, and J. Haegeman, MPSKit.jl, Zenodo (2026)

    2026

  53. [54]

    Devos and J

    L. Devos and J. Haegeman, TensorKit.jl: A Julia pack- age for large-scale tensor computations, with a hint of category theory (2025), arXiv:2508.10076

    Pith/arXiv arXiv 2025

  54. [55]

    Burgelman and P

    L. Burgelman and P. Brehmer, ImplicitDifferen- tiationPEPS.jl,https://github.com/leburgel/ ImplicitDifferentiationPEPS.jl(2026)

    2026

  55. [56]

    Some implementations evaluate the Sylvester equation in the pullback using the LAPACK-based Bartels–Stewart algorithm [54]. However, significantly improved perfor- mance can be achieved already at moderate dimensions by solving the Sylvester equation using appropriately pre- conditioned Krylov-based iterative methods [29] such as GMRES

  56. [57]

    Vanderstraeten, J

    L. Vanderstraeten, J. Haegeman, and F. Verstraete, Sim- ulating excitation spectra with projected entangled-pair states, Phys. Rev. B99, 165121 (2019). 29

    2019

  57. [58]

    R. N. C. Pfeifer, G. Evenbly, S. Singh, and G. Vi- dal, NCON: A tensor network contractor for MATLAB (2015), arXiv:1402.0939

    Pith/arXiv arXiv 2015

  58. [59]

    This approach is entirely equivalent, and can even be implemented by only tracking the ad- joint variables as defined above

    As an alternative to the approach outlined here, one could incorporate complex numbers into AD by treating all variables and their complex conjugates as separate degrees of freedom. This approach is entirely equivalent, and can even be implemented by only tracking the ad- joint variables as defined above

  59. [60]

    M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for con- tracting infinite two-dimensional tensor networks, Phys. Rev. B98, 235148 (2018)

    2018

  60. [61]

    Xie, J.-G

    H. Xie, J.-G. Liu, and L. Wang, Automatic differentiation of dominant eigensolver and its applications in quantum physics, Phys. Rev. B101, 245139 (2020)

    2020

  61. [62]

    Corboz, S

    P. Corboz, S. R. White, G. Vidal, and M. Troyer, Stripes in the two-dimensionalt-Jmodel with infinite projected entangled-pair states, Phys. Rev. B84, 041108 (2011)

    2011

  62. [63]

    R. H. Bartels and G. W. Stewart, Algorithm 432 [C2]: Solution of the matrix equation AX + XB = C [F4], Commun. ACM15, 820–826 (1972)

    1972