REVIEW 3 major objections 6 minor 38 references
Low-rank representation of tensor network operators with long-range pairwise interactions
T0 review · 3 major / 6 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read Long-range pairwise interactions compress into O(log N · log(N/ε)) tensor-network operators of constant bond dimension, without any translation-invariance assumption; all efficient MPO constructions are traced to one upper-triangular…
desk verdict A genuinely useful unification of MPO compression ideas and a new hierarchical construction with a plausible O(log N log(N/epsilon)) bound, held back by a few unproved error-accumulation claims that should be fixed in revision. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The paper runs on two matrix structures plus a stabilization step. The first structure is the upper-triangular low-rank (UTLR) property (Definition 3.1): a symmetric matrix $A\in\mathbb{R}^{N\times N}$ is UTLR if for every $\epsilon>0$ there is a rank-$r$ matrix, $r \lesssim \log(N/\epsilon)$, that approximates $A$ to relative accuracy $\epsilon$ on every block of upper-triangular entries, even though $A$ itself is full-rank; this property converts MPO construction into the matrix-completion problem $\min_{\mathrm{rank}(\tilde V)\le M} \|P_\Omega(\tilde V)-P_\Omega(V)\|_F^2$ over the upper triangle. The second structure is hierarchical low-rank blocking: in the HODLR format an interval interacts only with its sibling (weak admissibility), while in the H-matrix format interacting boxes must be separated by a distance comparable to their size (strong admissibility); each admissible block is compressed by truncated SVD, and the rank bound $M = \mathcal{O}(\log(N/\epsilon))$ for the Coulomb kernel comes from the classical fast-multipole separation estimate for $1/r$. The assembly step that turns compressed blocks into operators is the finite-state-machine rule for sums of non-overlapping rank-one operators, which produces an MPO or PEPO of bond dimension 5; summing $L = \log N$ levels of blocks gives the total count $\mathcal{O}(LM)$. The stabilization step is the added QR factorization in the modified incremental SVD, which keeps the row-by-row matrix-completion updates well-conditioned even though the singular values of the completed matrix grow exponentially with $N$.
What would settle it
Run the hierarchical construction on a coefficient matrix engineered so that every off-diagonal block has the same low-rank structure and all truncation residuals point in the same direction; double the system size while holding the per-block threshold at $C\epsilon/N$ and measure the global 2-norm error. If the global error grows faster than linearly with the number of blocks, the $\mathcal{O}(\log(N)\log(N/\epsilon))$ accuracy claim fails. A weaker test for the ISVD route: push the paper's quasi-1D non-translation-invariant experiment to larger $N$ and check whether the achieved MPO rank stays bounded by $\log(N/\epsilon)$; the paper provides no proof that it does.
Extended reading notes
Core claim
On its own terms the paper makes one diagnostic and one constructive claim. The diagnostic claim is that every efficient existing method for building a matrix product operator for a long-range pairwise interaction — exponential fitting by finite-state machines and the sliced-basis SVD compression method alike — exploits the same hidden structure: the coefficient matrix $V\in\mathbb{R}^{N\times N}$ is full-rank, yet for any block of upper-triangular indices the block can be matched to relative accuracy $\epsilon$ by a rank-$r$ matrix with $r \lesssim \log(N/\epsilon)$. This 'upper-triangular low-rank' (UTLR) view recasts MPO construction as a matrix completion problem, and the paper's modified incremental SVD (an incremental SVD with an added QR factorization) solves it without ever forming entries that grow exponentially with system size; the MPO it produces is equivalent to the one from the sliced-basis compression method. The constructive claim is that hierarchical low-rank matrix formats give a representation valid for both MPOs and PEPOs: each off-diagonal block of $V$ is compressed by truncated SVD to rank $\mathcal{O}(\log(N/\epsilon))$, each compressed block contributes a few rank-one 'non-overlapping' operators, and the finite-state-machine assembly of those operators gives an MPO or PEPO of bond dimension 5. With $L = \log N$ levels, the Coulomb interaction is thereby written as a linear combination of $\mathcal{O}(\log N \cdot \log(N/\epsilon))$ such operators, each of constant bond dimension, with an a priori error bound inherited from fast-multipole separation estimates for the $1/r$ kernel. Both routes are stated and demonstrated without assuming the interaction is translation-invariant.
Load-bearing premise
The efficiency bound assumes that the tiny errors made when each matrix block is squeezed separately do not pile up as the number of blocks grows; the paper sets the per-block tolerance to $C\epsilon/N$ to make this plausible but gives no proof that the errors accumulate only linearly — and for the incremental-SVD route it openly states that no theoretical error bound exists.
Editorial extensions
If this is right
- One-dimensional systems with Coulomb or other long-range pairwise potentials obtain matrix product operators of rank $\mathcal{O}(\log(N/\epsilon))$ from the modified incremental SVD, with accuracy comparable to — and in the paper's tests slightly better than — exponential fitting, even when the interaction is not translation-invariant.
- Two-dimensional lattice systems obtain PEPO representations of the Coulomb interaction as linear combinations of $\mathcal{O}(\log N \cdot \log(N/\epsilon))$ operators of constant bond dimension with an a priori error bound, a guarantee the earlier correlation-function PEPO approach did not offer.
- Because the hierarchical representation is a sum of independent constant-rank operators, the operator can be applied to a state term by term, a structure the paper notes is naturally suited to parallel computation.
- Numerical experiments in the paper show the per-block rank of the H-matrix representation grows only mildly with system size (maximal rank 13 at $16\times16$ up to 15 at $128\times128$ for threshold $10^{-5}$), consistent with the logarithmic bound.
- The two routes distribute error differently: the modified ISVD concentrates error at the far upper-right corner of the coefficient matrix, while the HODLR and H-matrix formats spread it across blocks, as shown in the paper's error-distribution plots.
Reading between the lines
- Any interaction kernel with a separation property — screened Coulomb, dipole–dipole, Yukawa, or other Green's functions — should inherit the same $\mathcal{O}(\log N \cdot \log(N/\epsilon))$ representation, since the rank argument rests only on blockwise separability of the $1/r$ kernel; running the H-matrix construction on such kernels would test this directly.
- The upper-triangular low-rank view suggests a recipe for compact MPOs beyond pairwise interactions: any operator whose coefficient array is full-rank but admits a low-rank extension on its triangular support can be fed to the modified incremental SVD; three-body interactions are the obvious next test.
- The missing error-accumulation proof is empirically checkable: holding the per-block threshold at $C\epsilon/N$ while doubling $N$, a flat global error would validate linear accumulation, while growth with the number of levels would force a smaller threshold and a larger constant in the complexity bound.
- The paper notes the sum-of-constant-rank-operators form is parallel-friendly but does not quantify the payoff; if combined with parallel time-evolution algorithms, the practical cost of simulating Coulomb-interacting systems could drop from quadratic in $N$ to near-linear with a logarithmic-factor overhead.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes two strategies for constructing efficient MPO/PEPO representations of long-range pairwise interactions, with the Coulomb interaction as the main example. The first strategy formalizes the observation that many existing MPO constructions rely on an upper-triangular low-rank (UTLR) property of the coefficient matrix, casts the construction as an ill-conditioned matrix completion problem, and introduces a modified incremental SVD (ISVD) method that avoids explicit computation of the exponentially growing lower-triangular part. The second strategy uses hierarchical low-rank matrix formats: the HODLR format for 1D MPOs and the H-matrix format for 1D MPOs and 2D PEPOs. The authors claim that for the Coulomb interaction the resulting representation is a linear combination of O(log(N) log(N/epsilon)) MPOs/PEPOs, each with constant bond dimension, and that neither algorithm requires translation invariance. Numerical experiments in quasi-1D and 2D systems are presented in support of the scaling claims.
Significance. If the main quantitative claim is established, the hierarchical low-rank construction is a valuable contribution: it would give the first a priori bounded construction for non-translation-invariant long-range interactions in both MPO and PEPO formats, with an explicit logarithmic scaling in system size and accuracy. The paper is also useful for clearly identifying the UTLR property as the common mechanism behind earlier exponential-fitting and SVD-compression MPO methods. The numerical experiments are consistent with the claimed scaling, and the authors are honest in Section 7 about the absence of an error bound for the modified ISVD. However, the advertised accuracy guarantee for the hierarchical route currently rests on an unproved error-accumulation assertion, and the relation between the coefficient-matrix error and the operator error is not specified. The central result is therefore plausible but not yet fully supported.
major comments (3)
- [Section 4.4, after Eq. (4.4)] The global error control for the H-matrix construction is asserted rather than proved. The text states that choosing the blockwise SVD threshold epsilon' = C epsilon/N ensures a coefficient-matrix 2-norm error below epsilon, and this leads to M = O(log(N/epsilon)). But the coefficient matrix is assembled from O(N log N) independently truncated blocks plus symmetric counterparts, and the spectral norm of a sum of localized perturbations is not bounded by the maximum blockwise error. At a fixed level the interaction blocks are disjoint and could be analyzed as block-diagonal perturbations, but the accumulation over hierarchical levels and the role of the factor N are not explained. This step is load-bearing because the abstract's O(log(N) log(N/epsilon)) accuracy claim depends on it. Please provide an explicit error-accumulation proof, or clearly state the result as a numerical observation with a modified claim.
- [Section 7] The conclusion states that for the Coulomb interaction the MPO rank from the modified ISVD method can be bounded by log(N/epsilon), while the same section admits that no theoretical error bound exists for the method. The non-translation-invariant UTLR route is therefore supported only by numerical evidence. Please either supply an analysis of the incremental completion error or rephrase the claim in Section 7 and in the introduction as a numerical observation rather than a bound.
- [Section 4.4 and Abstract] The accuracy measure is ambiguous. The construction in Section 4.4 controls the 2-norm error of the coefficient matrix V, while the abstract refers to 'epsilon as the accuracy of the low-rank truncation' without specifying whether this is the coefficient-matrix error or the operator error. Since the operator is \hat V = \sum_{i<j} V(i,j) \hat n_i \hat n_j and the local operators \hat n_i are not orthonormal in general, the operator-norm error is not automatically equal to the coefficient-matrix 2-norm error; a dimension-dependent factor can enter. The scaling claim may still be correct, but the norm used for epsilon needs to be stated precisely throughout, and the conversion between the two errors needs to be justified if the abstract's accuracy language is retained.
minor comments (6)
- [Section 4.2, Definition 4.1] The inequality in Definition 4.1 says ||A(I,J) - \tilde A_{I,J}||_2 <= epsilon ||\tilde A(I,J)||_2; for consistency with Definition 3.1 and with the surrounding text, the right-hand side should be ||A(I,J)||_2.
- [Section 3, after Eq. (3.4)] The sentence 'The matrix completion problem (3.5) has several features' refers to the MPO tensor displayed in Eq. (3.5); the reference should be to the optimization problem in Eq. (3.4).
- [Section 4.4] There is a typo: 'Sill we divide the interactions' should be 'Still we divide the interactions'.
- [Figure 6.3(a)] The axis labels and legend in Fig. 6.3(a) are garbled in the current text ('number of MPOs 1 2 3 rel 2-norm error 1e 4'); the curves, axes, and threshold values should be labeled unambiguously.
- [Algorithms 1 and 2] The notation W_p is reused for matrices of different sizes in the algorithmic description after the 'abuse of notation' remark; for reproducibility it would help to define the dimensions and the submatrix convention explicitly in each algorithm.
- [Section 4.5] The counts 78, 12, and 8 interactions in the 2D H-matrix construction are stated without derivation; a short counting argument for the reduction from 120 pairwise box interactions to the listed numbers would make the constant prefactor transparent.
Circularity Check
No significant circularity: the Coulomb rank bounds come from external exponential-sum and multipole theorems, and the MPO/PEPO constructions are explicit translations of those low-rank blocks.
full rationale
The paper's central rank bounds are not circular. The M = O(log(N/epsilon)) estimates for the Coulomb interaction are imported from external results: the exponential-sum bound in [3,6] (Eq. 5.1) and the multipole expansion in [16, Theorem 3.2] (Eq. 5.2). The UTLR and hierarchical-low-rank definitions (Definitions 3.1 and 4.1) encode the existence of low-rank approximations, but the paper does not derive Coulomb's low-rankness from those definitions; it invokes the external theorems. The MPO/PEPO constructions in Sections 3 and 4 explicitly convert low-rank factors into FSM-based cores (Appendices A and B), and the count O(log(N) log(N/epsilon)) follows from the number of levels times the external block-rank bound. The statement that the modified ISVD yields an MPO equivalent to [33] is a comparative claim, not the source of any predicted complexity bound, and the paper does not rely on self-citations by the authors for its main conclusions. Two weaknesses are present but are correctness or completeness gaps, not circularity: Section 4.4 asserts global 2-norm error control by choosing epsilon' = C epsilon/N without proving error accumulation over the O(N log N) blocks, and Section 7 explicitly states that no theoretical error bound exists for the modified ISVD algorithm. These gaps mean the advertised accuracy bound is not fully proven, but they do not reduce a derived result to an input by construction.
Assumptions & free parameters
free parameters (1)
- local SVD truncation threshold epsilon' =
epsilon' = C epsilon / N, constant C unspecified
assumptions (5)
- standard math The inverse distance 1/r can be approximated to relative accuracy epsilon by a sum of O(log(N/epsilon)) exponentials.
- standard math Well-separated Coulomb blocks admit low-rank approximations of rank O(log(N/epsilon)), from multipole expansions.
- domain assumption The coefficient matrix V for the interactions considered satisfies the UTLR property (Definition 3.1).
- domain assumption Strong admissibility (separation comparable to box size) is required for the Coulomb rank bounds in 2D.
- domain assumption The operator has the pairwise form V = sum_{i<j} V(i,j) n_i n_j with a symmetric coefficient matrix V.
Cite this review
Pith. "Pith review of Low-rank representation of tensor network operators with long-range pairwise interactions." pith.science (2026). https://pith.science/paper/D4XNCYVJ
@misc{pith2026190902206,
author = {Pith},
title = {Pith review of: Low-rank representation of tensor network operators with long-range pairwise interactions},
year = {2026},
howpublished = {\url{https://pith.science/paper/D4XNCYVJ}},
note = {Machine review of arXiv:1909.02206}
}
abstract
Tensor network operators, such as the matrix product operator (MPO) and the projected entangled-pair operator (PEPO), can provide efficient representation of certain linear operators in high dimensional spaces. This paper focuses on the efficient representation of tensor network operators with long-range pairwise interactions such as the Coulomb interaction. For MPOs, we find that all existing efficient methods exploit a peculiar "upper-triangular low-rank" (UTLR) property, i.e. the upper-triangular part of the matrix can be well approximated by a low-rank matrix, while the matrix itself can be full-rank. This allows us to convert the problem of finding the efficient MPO representation into a matrix completion problem. We develop a modified incremental singular value decomposition method (ISVD) to solve this ill-conditioned matrix completion problem. This algorithm yields equivalent MPO representation to that developed in [Stoudenmire and White, Phys. Rev. Lett. 2017]. In order to efficiently treat more general tensor network operators, we develop another strategy for compressing tensor network operators based on hierarchical low-rank matrix formats, such as the hierarchical off-diagonal low-rank (HODLR) format, and the $\mathcal{H}$-matrix format. Though the pre-constant in the complexity is larger, the advantage of using the hierarchical low-rank matrix format is that it is applicable to both MPOs and PEPOs. For the Coulomb interaction, the operator can be represented by a linear combination of $\mathcal{O}(\log(N)\log(N/\epsilon))$ MPOs/PEPOs, each with a constant bond dimension, where $N$ is the system size and $\epsilon$ is the accuracy of the low-rank truncation. Neither the modified ISVD nor the hierarchical low-rank algorithm assumes that the long-range interaction takes a translation-invariant form.
Figures
Figures from the paper (8 more)
Reference graph
Works this paper leans on
-
[1]
S. Ambikasaran and E. Darve , AnO(n logn) fast direct solver for partial hierarchically semi-separable matrices, J. Sci. Comput., 57 (2013), pp. 477–501
work page 2013
-
[2]
L. Balzano and Stephen J Wright , On grouse and incremental svd , in 5th IEEE Interna- tional Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAM- SAP), IEEE, 2013, pp. 1–4
work page 2013
-
[3]
G. Beylkin and M. J. Mohlenkamp, Numerical operator calculus in higher dimensions, Proc. Natl. Acad. Sci. U.S.A., 99 (2002), pp. 10246–10251
work page 2002
-
[4]
G. Beylkin and L. Monz´on, On approximation of functions by exponential sums , Appl. Com- put. Harmon. A., 19 (2005), pp. 17–48. 30 rule num (l,u,d,r ) operator UL of R 1 (▽,∆,∆,∆) Ii 2 (∆,▽,∆,∆) Ii 3 (∆,∆,d,r ) Hi(1, 1,d,r ) U of R 4 (∆,∆,∆,∆) Ii 5 (l,∆,d,r ) Hi(l, 1,d,r ) UR of R 6 (∆,∆,∆,∆) Ii 7 (l,∆,d,∆) Hi(l, 1,d, 1) R of R 8 (∆,∆,∆,∆) Ii 9 (l,u,d, ∆...
work page 2005
-
[5]
J. Biamonte and V. Bergholm , Tensor networks in a nutshell , arXiv preprint arXiv:1708.00006, (2017)
arXiv 2017
-
[6]
D. Braess and W. Hackbusch , Approximation of 1/x by exponential sums in [1,∞), IMA J. Numer. Anal., 25 (2005), pp. 685–697
work page 2005
-
[7]
, On the efficient computation of high-dimensional integrals and the approximation by exponential sums , in Multiscale, nonlinear and adaptive approximation, Springer, 2009, pp. 39–74
work page 2009
-
[8]
Matthew Brand, Incremental singular value decomposition of uncertain data with missing values, in European Conference on Computer Vision, Springer, 2002, pp. 707–720
work page 2002
Show all 38 references
-
[9]
J.-F. Cai, E. J. Cand `es, and Z. Shen , A singular value thresholding algorithm for matrix completion, SIAM J. Optim., 20 (2010), pp. 1956–1982
2010
-
[10]
E. J. Cand`es and B. Recht , Exact matrix completion via convex optimization , Found. Com- put. Math., 9 (2009), p. 717
2009
-
[11]
G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White , Matrix product opera- tors, matrix product states, and ab initio density matrix renormalization group algorithms , J. Chem. Phys., 145 (2016), p. 014102
2016
-
[12]
G. M. Crosswhite and D. Bacon, Finite automata for caching in matrix product algorithms , Phys. Rev. A, 78 (2008), p. 012356
2008
-
[13]
M Crosswhite, A
G. M Crosswhite, A. C. Doherty, and G. Vidal , Applying matrix product operators to model systems with long-range interactions , Phys. Rev. B, 78 (2008), p. 035116
2008
-
[14]
Fr¨owis, V
F. Fr¨owis, V. Nebendahl, and W. D ¨ur, Tensor operators: Constructions and applications for long-range interaction systems , Phys. Rev. A, 81 (2010), p. 062337
2010
-
[15]
Grasedyck and W
L. Grasedyck and W. Hackbusch, Construction and arithmetics of h-matrices , Computing, 70 (2003), pp. 295–334
2003
-
[16]
Greengard and V
L. Greengard and V. Rokhlin , A new version of the fast multipole method for the laplace equation in three dimensions , Acta Numer., 6 (1997), pp. 229–269
1997
-
[17]
Hackbusch, A sparse matrix arithmetic based on H-matrices
W. Hackbusch, A sparse matrix arithmetic based on H-matrices. part i: Introduction to H- matrices, Computing, 62 (1999), pp. 89–108
1999
-
[18]
Hackbusch, B
W. Hackbusch, B. N. Khoromskij, and R. Kriemann, Hierarchical matrices based on a weak admissibility criterion, Computing, 73 (2004), pp. 207–243
2004
-
[19]
Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang , Unsupervised generative modeling 31 using matrix product states , Phys. Rev. X, 8 (2018), p. 031012
2018
-
[20]
Y. Hu, Y. Koren, and C. Volinsky , Collaborative filtering for implicit feedback datasets. , in ICDM, vol. 8, Citeseer, 2008, pp. 263–272
2008
-
[21]
Lubich, T
C. Lubich, T. Rohwedder, R. Schneider, and B. Vandereycken , Dynamical approxima- tion by hierarchical tucker and tensor-train tensors , SIAM J. Matrix Anal. A., 34 (2013), pp. 470–494
2013
-
[22]
Mnih and R
A. Mnih and R. R. Salakhutdinov , Probabilistic matrix factorization, in Advances in neural information processing systems, 2008, pp. 1257–1264
2008
-
[23]
M. R. Norman, Colloquium: Herbertsmithite and the search for the quantum spin liquid , Rev. Mod. Phys., 88 (2016), p. 041002
2016
-
[24]
M. J. O’Rourke, Z. Li, and G. K.-L. Chan, Efficient representation of long-range interactions in tensor network algorithms , Phys. Rev. B, 98 (2018), p. 205127
2018
-
[25]
Or ´us, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann
R. Or ´us, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys., 349 (2014), pp. 117–158
2014
-
[26]
I. V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317
2011
-
[27]
I. V. Oseledets and E. E. Tyrtyshnikov , Breaking the curse of dimensionality, or how to use svd in many dimensions , SIAM J. Sci. Comput., 31 (2009), pp. 3744–3759
2009
-
[28]
¨Ostlund and S
S. ¨Ostlund and S. Rommer , Thermodynamic limit of density matrix renormalization , Phys. Rev. Lett., 75 (1995), p. 3537
1995
-
[29]
Pirvu, V
B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, Matrix product operator representations, New J. Phys., 12 (2010), p. 025012
2010
-
[30]
Rakhuba and I
M. Rakhuba and I. Oseledets, Calculating vibrational spectra of molecules using tensor train decomposition, J. Chem. Phys., 145 (2016), p. 124101
2016
-
[31]
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., 326 (2011), pp. 96–192
2011
-
[32]
Stoudenmire and D
E. Stoudenmire and D. J. Schwab , Supervised learning with tensor networks , in Advances in Neural Information Processing Systems, 2016, pp. 4799–4807
2016
-
[33]
E. M. Stoudenmire and S. R. White , Sliced basis density matrix renormalization group for electronic structure, Phys. Rev. Lett., 119 (2017), p. 046401
2017
-
[34]
Verstraete and J
F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions , arXiv preprint cond-mat/0407066, (2004)
2004 arXiv
-
[35]
Verstraete, M
F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac , Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett., 96 (2006), p. 220601
2006
-
[36]
Vidal, Classical simulation of infinite-size quantum lattice systems in one spatial dimen- sion, Phys
G. Vidal, Classical simulation of infinite-size quantum lattice systems in one spatial dimen- sion, Phys. Rev. Lett., 98 (2007), p. 070201
2007
-
[37]
S. R. White , Density matrix formulation for quantum renormalization groups , Phys. Rev. Lett., 69 (1992), p. 2863
1992
-
[38]
S. Yan, D. A. Huse, and S. R. White, Spin-liquid ground state of the s= 1/2 kagome heisen- berg antiferromagnet, Science, 332 (2011), pp. 1173–1176
2011
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.