Low-rank representation of tensor network operators with long-range pairwise interactions

Lin Lin Department of Mathematics, University of California, Berkeley, and Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720. Email: linlin@math.berkeley.edu Yu Tong Department of Mathematics, University of California, Berkeley, CA 94720. Email: yu_tong@berkeley.edu

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.

keywords:

Matrix product operator, projected entangled-pair operator, upper-triangular low-rank matrix, hierarchical off-diagonal low-rank,

\mathcal{H}

-matrix, fast multipole method

AMS:

15A69, 41A99, 65Z05

1 Introduction

Tensor network states [37, 36], particularly presented by the matrix product states (MPS) [37, 28, 31] and projected entangled-pair states (PEPS) [34, 35, 25], are among the most promising classes of variational methods for approximating high-dimensional functions in quantum physics. Besides their successes for treating strongly correlated quantum systems [38, 23], the MPS (also known as the tensor train method (TT) [27, 26, 21]), have become useful in a wide range of applications [30, 32, 19]. A core component of tensor network algorithms is the efficient representation of linear operators. In MPS the operator representation is called the matrix product operator (MPO), and in PEPS the projected entangled-pair operator (PEPO). These forms of representation are naturally designed for short-range interactions, such as interactions of nearest-neighbor type on a regular lattice. For long-range interactions, such as those naturally appearing from electronic structure calculations due to the Coulomb interaction, a straightforward representation would lead to a large MPO/PEPO rank, which can be prohibitively expensive. Therefore efficient representation of tensor network operators with long-range interactions is crucial for the success of the methods.

In this paper, we focus on the tensor network operators with long-range pairwise interactions, such as the Coulomb interaction. A tensor network operator with pairwise interaction can be defined as

(1)

\hat{V}=\frac{1}{2}\sum_{1\leq i\neq j\leq N}\mathbf{V}(i,j)\hat{n}_{i}\hat{n}_{j}=\sum_{1\leq i<j\leq N}\mathbf{V}(i,j)\hat{n}_{i}\hat{n}_{j}.

Here $\mathbf{V}\in\mathbb{R}^{N\times N}$ is a symmetric matrix, called the coefficient matrix of $\hat{V}$ , and $\hat{n}_{i}$ is called a number operator. The precise definition of $\hat{n}_{i}$ , as well as the connection to the Coulomb interaction will be discussed in Section 2.

There have been a number of proposals to treat MPOs with long-range pairwise interactions. Using the exponential fitting techniques based on finite state machines (FSM) [12, 29, 14, 11, 13], the Coulomb interaction and other translation-invariant long-range interactions on one-dimensional lattice systems can be efficiently represented. For the Coulomb interaction, the MPO rank can be bounded by $\log(N/\epsilon)$ , where $N$ is the number of sites and $\epsilon$ is the target accuracy. For one-dimensional systems without the translation-invariance property, recently a method based on a sequence of singular value decomposition (SVD) compression operations [11, 33] has been developed. This method is observed to yield a compact MPO representation for the Coulomb interaction with basis functions arising from electronic structure calculations that are not translation-invariant. On the other hand, this is a greedy method, and there is no a priori upper bound for the MPO rank. For two-dimensional and higher dimensional lattice systems, PEPOs can provide more efficient representation of linear operators than MPOs. To our knowledge, the only efficient PEPO representation for long-range interaction is given by the recently developed correlation function valued PEPO (CF-PEPO) [24]. CF-PEPO uses a fitting method, which maps the original lattice system with translation-invariant long-range interactions to an effective system defined on a superlattice but with short-range interactions. There is no a priori upper bound for the PEPO rank. Neither is it clear whether the CF-PEPO method can be efficiently generalized to non-translation-invariant systems.

Contribution: The contribution of this paper is two-fold. First, we find that the success of both the exponential fitting method and the SVD compression method rests on the assumption that $\mathbf{V}$ satisfies the “upper-triangular low-rank” (UTLR) property, i.e. the upper-triangular part of $\mathbf{V}$ can be extended into a matrix that is approximately low-rank, while $\mathbf{V}$ can be a full-rank matrix. Therefore the problem can be viewed as a special type of matrix completion problem [10]. However, this matrix completion problem is highly ill-conditioned, and standard matrix completion algorithms suffer from numerical stability problems [9, 20, 22, 8, 2]. Based on the work of [8], we develop a modified incremental singular value decomposition method (ISVD) to simultaneously find the low-rank structure in a numerically stable fashion, and to construct the MPO representation. The algorithm yields an MPO representation equivalent to that in the SVD compression method [33]. The ISVD method is not restricted to translation-invariant systems, and numerical results indicate that the performance of ISVD is comparable to (in fact is marginally better than) that of the exponential fitting methods for translation-invariant systems.

Second, we propose a new way for representing long-range interactions based on the framework of hierarchical low-rank matrices. In particular, we focus on the hierarchical off-diagonal low-rank (HODLR) format [1], and the $\mathcal{H}$ -matrix format [15, 17]. The main advantage of the hierarchical low-rank format is that it allows us to construct efficient representation for both MPOs and PEPOs. For the Coulomb interaction, we can represent $\hat{V}$ using a linear combination of $\mathcal{O}(\log(N)\log(N/\epsilon))$ MPOs/PEPOs, and the bond dimension of each MPO/PEPO is bounded by a constant. From a computational perspective, such a format can be more preferable than a single MPO/PEPO operator with bond dimension bounded by $\mathcal{O}(\log(N)\log(N/\epsilon))$ . Furthermore, the hierarchical low-rank format can also be applied to systems without the translation-invariance property.

Notation:

Throughout the paper all vectors and operators are defined on finite-dimensional spaces. We use hatted letters such as $\hat{A},\hat{V}$ to denote high-dimensional linear operators. Vectors, matrices and tensor cores are denoted by bold font letters such as $\mathbf{A},\mathbf{X},\mathbf{V}$ . The tensor product of $\mathbf{A}$ and $\mathbf{B}$ is denoted by $\mathbf{A}\otimes\mathbf{B}$ . We use the matlab notation such as $\mathbf{A}(\mathcal{I},\mathcal{J})$ to denote sub-matrices, where $\mathcal{I},\mathcal{J}$ are index sets. Given a certain global ordering of all indices, the notation $\mathcal{I}\prec\mathcal{J}$ means that for any $i\in\mathcal{I}$ and $j\in\mathcal{J}$ , we have $i<j$ .

Organization: The rest of the paper is organized as follows. In Section 2 we briefly introduce the MPO and PEPO representation. In Section 3 we introduce and use the UTLR structure to construct MPO representation of long-range interaction for 1D systems. In Section 4 we introduce the method to construct MPO/PEPO representation using the hierarchical low-rank matrix format. Numerical results are presented in Section 6. We conclude our work in Section 7. The rules of finite state machine (FSM) for representing non-overlapping MPOs and PEPOs are given in the Appendix A and B, respectively.

2 Preliminaries

We first briefly review the construction of MPS/MPO, as well as of PEPS/PEPO. For more detailed information we refer readers to [31, 25]. A single vector $\mathbf{X}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ can be interpreted as an order- $d$ tensor, where the $d$ -tuple $(n_{1},\ldots,n_{d})$ is called the size of the tensor. We assume that the index follows a lexicographic ordering (a.k.a. row-major ordering). Each entry can be accessed using multi-indices $\mathbf{X}(i_{1},\ldots,i_{d}),1\leq i_{\alpha}\leq n_{\alpha},\alpha=1,\ldots,d$ . A linear operator $\hat{A}\in\mathbb{R}^{(m_{1}\times\cdots\times m_{d})\times(n_{1}\times\cdots\times n_{d})}$ is an order $2d$ -tensor, with each entry denoted by

\hat{A}(i_{1},\ldots,i_{d};j_{1},\ldots,j_{d}).

The application of $\hat{A}$ to $\mathbf{X}$ gives a tensor $\mathbf{Y}\in\mathbb{R}^{m_{1}\times\cdots\times m_{d}}$ as

\mathbf{Y}(i_{1},\ldots,i_{d})=\sum_{j_{1}=1}^{n_{1}}\cdots\sum_{j_{d}=1}^{n_{d}}\hat{A}(i_{1},\ldots,i_{d};j_{1},\ldots,j_{d})\mathbf{X}(j_{1},\ldots,j_{d}),

which can be written in short-hand notation as $\mathbf{Y}=\hat{A}\mathbf{X}$ .

In (1), we used $i,j$ as the site indices and $d=N$ . The sites can be organized into lattices of one, two and three dimensions. For Coulomb interaction in 1D $\mathbf{V}(i,j)=1/|i-j|$ for $i\neq j$ , and we set $\mathbf{V}(i,i)=0$ . In higher dimensions we use the notation $\mathbf{i},\mathbf{j}$ to represent a lattice site. For example, on an $N\times N$ 2D lattice we have $\mathbf{i}=(i_{x},i_{y})$ for $1\leq i_{x},i_{y}\leq N$ . For a pairwise interaction in the form of (1), the coefficient matrix has as its entries $\mathbf{V}(i_{x}+Ni_{y},j_{x}+Nj_{y})$ , using a row-major order, and this is written more compactly as $\mathbf{V}(\mathbf{i},\mathbf{j})$ hereafter. The Coulomb interaction takes the form $\mathbf{V}(\mathbf{i},\mathbf{j})=\frac{1}{\|\mathbf{i}-\mathbf{j}\|}$ , where $\|\cdot\|$ is the Euclidean distance. Using this notation, the pairwise interaction is called translation-invariant if $\mathbf{V}(\mathbf{i},\mathbf{j})$ is a function of $\mathbf{i}-\mathbf{j}$ (in the presence of periodic boundary conditions, we may simply redefine $\mathbf{i}-\mathbf{j}$ ). In Eq. (1), we may assume for simplicity that $n_{1}=\cdots=n_{d}=n$ . $\hat{n}_{i}$ acts only on a single-site $i$ and should be interpreted as $\mathbf{I}^{\otimes(i-1)}\otimes\mathbf{n}\otimes\mathbf{I}^{\otimes(n-i+1)}$ , where $\mathbf{I}$ is the identity matrix of size $n$ and $\mathbf{n}$ is an $n\times n$ matrix. $\hat{n}_{i}$ can be understood as a spin operator for quantum spin systems, or a number operator for quantum fermionic systems, though the precise form of $\mathbf{n}$ is not relevant for the purpose of this paper.

2.1 Matrix product operators

A vector $\mathbf{X}$ is represented by an MPS/TT format if each entry can be represented as as a matrix product

(2)

\mathbf{X}(i_{1},\ldots,i_{d})=\mathbf{U}_{1}(i_{1})\cdots\mathbf{U}_{d}(i_{d}),

where each $\mathbf{U}_{\alpha}(i_{\alpha})$ is a matrix of size $(r_{\alpha-1},r_{\alpha})$ . Since $\mathbf{X}(i_{1},\ldots,i_{d})$ is a scalar, by definition $r_{0}=r_{d}=1$ . Each $\mathbf{U}_{\alpha}$ , called a core tensor, can be viewed as a $3$ -tensor of size $(r_{\alpha-1},n_{\alpha},r_{\alpha})$ , and the matrix $\mathbf{U}_{\alpha}(i_{\alpha})$ is called the $i_{\alpha}$ -th slice of $\mathbf{U}_{\alpha}$ . Since the first and third indices of $\mathbf{U}_{\alpha}$ are to be contracted out via matrix-multiplication operations, they are called the internal indices, while the second index is called the external index. The $(d+1)$ -tuple $(r_{0},r_{1},\ldots,r_{d})$ is called the bond dimension (a.k.a. the MPO/TT rank). Sometimes the bond dimension also refers to a single scalar $\max_{1\leq\alpha\leq d}r_{\alpha}$ .

The MPO format of a linear operator $\hat{A}$ is

(3)

\hat{A}(i_{1},\ldots,i_{d};j_{1},\ldots,j_{d})=\mathbf{A}_{1}\left(i_{1},j_{1}\right)\mathbf{A}_{2}\left(i_{2},j_{2}\right)\cdots\mathbf{A}_{d}\left(i_{d},j_{d}\right).

Each $\mathbf{A}_{\alpha}$ is also called a core tensor, and is of size $(s_{\alpha-1},m_{\alpha},n_{\alpha},s_{\alpha})$ . The first and the fourth indices are called the internal indices, and the second and third indices the external indices. The matrix $\mathbf{A}_{\alpha}(i_{1},j_{1})\in\mathbb{R}^{m_{\alpha}\times n_{\alpha}}$ is called a slice. Again $s_{0}=s_{d}=1$ . The $(d+1)$ -tuple $(s_{0},\ldots,s_{d})$ , or sometimes simply $\max_{1\leq\alpha\leq d}s_{\alpha}$ , is called the bond dimension.

Consider the application of the linear operator $\hat{A}$ on $\mathbf{X}$ , denoted by $\mathbf{Y}=\hat{A}\mathbf{X}$ , then the vector $\mathbf{Y}$ has an MPS representation as

\mathbf{Y}(i_{1},\ldots,i_{d})=\mathbf{V}_{1}(i_{1})\cdots\mathbf{V}_{d}(i_{d}).

Using the MPO representation, we may readily find that

\mathbf{V}_{\alpha}(i_{\alpha})=\sum_{j_{\alpha}=1}^{n_{\alpha}}\mathbf{A}_{\alpha}(i_{\alpha},j_{\alpha})\otimes\mathbf{U}(j_{\alpha}).

The bond dimension of $\mathbf{Y}$ is $(s_{0}r_{0},s_{1}r_{1},\ldots,s_{d}r_{d})$ .

The simplest example of the MPO is an operator in the tensor product form, i.e.

\hat{A}=\mathbf{a}_{1}\otimes\cdots\otimes\mathbf{a}_{n},

where $\mathbf{a}_{\alpha}\in\mathbb{R}^{n_{\alpha}\times n_{\alpha}}$ is a matrix, i.e. $\mathbf{a}_{\alpha}(i_{\alpha},j_{\alpha})$ is a scalar. In the component form

\hat{A}(i_{1},\ldots,i_{d};j_{1},\ldots,j_{d})=\mathbf{a}_{1}\left(i_{1},j_{1}\right)\mathbf{a}_{2}\left(i_{2},j_{2}\right)\cdots\mathbf{a}_{d}\left(i_{d},j_{d}\right).

Clearly the MPO bond dimension is $1$ .

The second example is an operator with nearest-neighbor interaction, e.g.

\hat{A}=\sum_{\alpha=1}^{d-1}(\mathbf{I}_{1}\otimes\cdots\otimes\mathbf{I}_{\alpha-1})\otimes\mathbf{a}_{\alpha}\otimes\mathbf{a}_{\alpha+1}\otimes(\mathbf{I}_{\alpha+2}\otimes\cdots\otimes\mathbf{I}_{d}).

where $\mathbf{I}_{\alpha}$ is the identity matrix of size $n_{\alpha}$ . In the component form

(4)

\begin{split}&\hat{A}(i_{1},\ldots,i_{d};j_{1},\ldots,j_{d})\\ =&\sum_{\alpha=1}^{d-1}(\delta_{i_{1},j_{1}}\cdots\delta_{i_{\alpha-1},j_{\alpha-1}})\mathbf{a}_{\alpha}\left(i_{\alpha},j_{\alpha}\right)\mathbf{a}_{\alpha+1}\left(i_{\alpha+1},j_{\alpha+1}\right)(\delta_{i_{\alpha+2},j_{\alpha+2}}\cdots\delta_{i_{d},j_{d}}).\end{split}

$\hat{A}$ can be viewed as the linear combination of $(d-1)$ MPOs in the tensor product form, and if so the resulting MPO rank would be $(d-1)$ . However, note that we may define

\mathbf{A}_{1}(i_{1},j_{1})=\begin{pmatrix}0&\mathbf{a}_{1}(i_{1},j_{1})&\delta_{i_{1},j_{1}}\end{pmatrix},\quad\mathbf{A}_{d}(i_{d},j_{d})=\begin{pmatrix}\delta_{i_{d},j_{d}}\\ \mathbf{a}_{d}(i_{d},j_{d})\\ 0\end{pmatrix},

and

\mathbf{A}_{\alpha}(i_{\alpha},j_{\alpha})=\begin{pmatrix}\delta_{i_{\alpha},j_{\alpha}}&0&0\\ \mathbf{a}_{\alpha}(i_{\alpha},j_{\alpha})&0&0\\ 0&\mathbf{a}_{\alpha}(i_{\alpha},j_{\alpha})&\delta_{i_{\alpha},j_{\alpha}}\end{pmatrix},\quad 2\leq\alpha\leq d-1.

Then it can be readily verified that the MPO of the form $\eqref{eqn:MPO}$ that agrees with the component form $\eqref{eqn:MPONNcomponent}$ . Hence the MPO rank is independent of $d$ and is always $3$ .

When the context is clear, we may also identify $\hat{a}_{\alpha}$ with $(\mathbf{I}_{1}\otimes\cdots\otimes\mathbf{I}_{\alpha-1})\otimes\mathbf{a}_{\alpha}\otimes(\mathbf{I}_{\alpha+1}\otimes\cdots\otimes\mathbf{I}_{d})$ . Then

\hat{A}=\sum_{\alpha=1}^{d-1}\hat{a}_{\alpha}\hat{a}_{\alpha+1}.

In the MPO form, we may also omit the $(i,j)$ indices and write

\mathbf{A}_{1}=\begin{pmatrix}0&\mathbf{a}_{1}&\mathbf{I}_{1}\end{pmatrix},\quad\mathbf{A}_{d}=\begin{pmatrix}\mathbf{I}_{d}\\ \mathbf{a}_{d}\\ 0\end{pmatrix},\quad\mathbf{A}_{\alpha}=\begin{pmatrix}\mathbf{I}_{\alpha}&0&0\\ \mathbf{a}_{\alpha}&0&0\\ 0&\mathbf{a}_{\alpha}&\mathbf{I}_{\alpha}\end{pmatrix},\quad 2\leq\alpha\leq d-1.

The third example is an operator with a special form of long-range interaction:

\hat{A}=\sum_{1\leq\alpha<\beta\leq d}e^{\lambda(\alpha-\beta)}\hat{a}_{\alpha}\hat{a}_{\beta}.

We assume $\lambda>0$ , and omit the component form for simplicity. The length of the interaction is characterized by $1/\lambda$ . A straightforward term-by-term representation would lead to an MPO of rank $\mathcal{O}(d^{2})$ . Nonetheless, using the fact that

e^{-\lambda n}=\underbrace{e^{-\lambda}e^{-\lambda}\cdots e^{-\lambda}}_{n~\text{times}},\quad n\in\mathbb{N}_{+},

$\hat{A}$ can also be written as an MPO with rank- $3$ as

\mathbf{A}_{1}=\begin{pmatrix}0&\mathbf{a}_{1}&\mathbf{I}_{1}\end{pmatrix},\quad\mathbf{A}_{d}=\begin{pmatrix}\mathbf{I}_{d}\\ e^{-\lambda}\mathbf{a}_{d}\\ 0\end{pmatrix},\quad\mathbf{A}_{\alpha}=\begin{pmatrix}\mathbf{I}_{\alpha}&0&0\\ e^{-\lambda}\mathbf{a}_{\alpha}&e^{-\lambda}\mathbf{I}_{\alpha}&0\\ 0&\mathbf{a}_{\alpha}&\mathbf{I}_{\alpha}\end{pmatrix},\quad 2\leq\alpha\leq d-1.

In physics literature, this is a special case of the finite state machine (FSM) [13]. For more complex settings, the FSM rule often specifies the row/column indices of each core tensor as input/output signals. We will use the language of input/output signals only in the Appendices A and B.

In order to proceed with the discussion of PEPS/PEPO, it is no longer productive to keep using the component form. Instead, the graphical representation of tensors is preferred. Fig. 1 shows an example of the graphical representation of tensors and tensor contraction operations. A tensor is visualized as a vertex and each index (internal and external) is represented as an edge. When two vertices are linked by an edge, the corresponding index is contracted out. A more detailed introduction can be found in [5]. Using this representation, MPS and MPO can be represented using graphs in Fig. 2 (a) and (b) respectively.

2.2 Projected entangled-pair operators

The MPS and MPO impose an intrinsically one-dimensional structure on the tensor following the index $1,\ldots,d$ . This is a very efficient representation for certain problems defined on a one-dimensional lattice. For problems defined on a two-dimensional lattice and beyond, the PEPS and PEPO often provide a more efficient representation of vectors and linear operators, respectively. The PEPS/PEPO can be most conveniently represented using the graphical representation (see Fig. 3 (a) and (b)).

In the two-dimensional case, each core tensor in the PEPS/PEPO format has up to 5/6 indices, respectively. Similar to what we have done for MPS/MPO, the (up to 4) contracted indices are called internal indices, while remaining ones are called external indices.

3 Long-range interaction MPO in 1D systems

In this section we focus on a 1D system consisting of $N$ sites. The simplest example of constructing a long-range interaction MPO in 1D system is the exponential fitting method [29, 12, 14]. We briefly describe this method used for the Coulomb interaction. First we approximate the inverse-distance by the sum of exponentials

(5)

\frac{1}{r}\approx\sum_{k=1}^{M}a_{k}e^{-\lambda_{k}r},\quad 1\leq r\leq N.

Then the tensor corresponding to the Coulomb interaction can be approximated by

\hat{V}\approx\sum_{k=1}^{M}a_{k}\sum_{i<j}e^{-\lambda_{k}(j-i)}\hat{n}_{i}\hat{n}_{j},

where $a_{k},\lambda_{k}\geq 0$ . This can be done analytically through the quadrature of an integral representation [7], or numerically through a least squares procedure. Following the discussion of the finite state machine in Section 2.1, the operator on the right-hand side can be represented by an MPO of bond dimension $M+2$ , the MPO tensor of which at site $i$ has the form

(6)

\left(\begin{array}[]{ccc}\mathbf{I}&0&0\\ \mathbf{b}\otimes\mathbf{n}_{i}&e^{-\mathbf{\Lambda}}\otimes\mathbf{I}&0\\ 0&\mathbf{a}^{T}\otimes\mathbf{n}_{i}&\mathbf{I}\end{array}\right),

where $\mathbf{a}=(a_{1},\cdots,a_{M})^{T}$ , $\mathbf{b}=(e^{-\lambda_{1}},\cdots,e^{-\lambda_{M}})^{T}$ , and $\mathbf{\Lambda}=\mathop{\operator@font diag}\nolimits(\lambda_{1},\cdots,\lambda_{M})$ . This is true for $2\leq i\leq N-1$ . For the tensors at the beginning and end of the MPO we only need to use the last row and first column of the above matrix respectively. Then in order to achieve a target accuracy $\epsilon$ for a system of size $N$ , we only need $M=\mathcal{O}(\log(N/\epsilon))$ [3, 6].

For the Coulomb interaction, there are two main methods for performing the exponential fitting procedure. One way is to represent $1/r$ as an integral, e.g.

\frac{1}{r}=\int^{\infty}_{0}e^{-\lambda r}\mathrm{d}\lambda,

and then approximate the right hand side with a quadrature scheme with $M$ points, such as the Gauss quadrature. This procedure results in pointwise error of $\epsilon$ in a fixed size finite interval $[a,b]$ where $0<a<b<\infty$ with $M$ upper bounded by $\mathcal{O}(\log(1/\epsilon))$ terms [4, 7, 6]. Another way is to solve a nonlinear least squares problem to find the fitting parameters. However, this optimization problem can have many stationary points, making finding the global minimum difficult. We can only start from different reasonable initial guesses to increase the chance of obtaining the global minimum.

The exponential fitting method is efficient when $\mathbf{V}$ is strictly decaying and translation-invariant. When either of these two conditions is not satisfied, this method needs some modification. There are methods to deal with the former by introducing complex exponents [29]. The generalization to non-translation-invariant systems can be more difficult.

In this work we propose a new perspective based on matrix completion, which naturally generalizes the exponential fitting method to situations where $\mathbf{V}$ is neither strictly decaying nor translation-invariant. We define a rank $M$ matrix $\tilde{\mathbf{V}}$ by

(7)

\tilde{\mathbf{V}}(i,j)=\sum_{k=1}^{M}a_{k}e^{-\lambda_{k}(j-i)}=\sum_{k=1}^{M}a_{k}e^{-\lambda_{k}j}e^{\lambda_{k}i}.

The success of the exponential fitting method rests on the fact that $\tilde{\mathbf{V}}$ can agree very well with $\mathbf{V}$ in the upper-triangular part, but can differ with $\mathbf{V}$ arbitrarily on the diagonal part and the lower-triangular part. Indeed, the lower-triangular part of $\tilde{\mathbf{V}}$ grows exponentially with respect to $(i-j)$ , while that of $\mathbf{V}$ decays with respect to $(i-j)$ due to the symmetry of $\mathbf{V}$ . Despite the difference, the tensor $\sum_{i<j}\tilde{\mathbf{V}}(i,j)\hat{n}_{i}\hat{n}_{j}$ very well approximates $\sum_{i<j}\mathbf{V}(i,j)\hat{n}_{i}\hat{n}_{j}$ .

More generally, we shall demonstrate that it is possible to construct an efficient MPO representation of $\hat{V}$ , if the matrix $\mathbf{V}$ satisfies the following upper triangular low rank (UTLR) property.

Definition 1 (Upper-triangular low-rank matrix).

A symmetric matrix $\mathbf{A}\in\mathbb{R}^{N\times N}$ is an upper-triangular low-rank (UTLR) matrix, if for any $\epsilon>0$ , there exists $r\lesssim\log(N/\epsilon)$ , and a rank- $r$ matrix $\tilde{\mathbf{A}}$ , such that for any index sets $\mathcal{I}\prec\mathcal{J}$ ,

\lVert\mathbf{A}(\mathcal{I},\mathcal{J})-\tilde{\mathbf{A}}(\mathcal{I},\mathcal{J})\rVert_{2}\leq\epsilon\lVert\mathbf{A}(\mathcal{I},\mathcal{J})\rVert_{2}.

For a UTLR matrix $\mathbf{V}$ , we may find its rank- $M$ approximation by solving the following optimization problem

(8)		$\displaystyle\min_{\tilde{\mathbf{V}}}$	$\displaystyle\\|P_{\Omega}(\tilde{\mathbf{V}})-P_{\Omega}(\mathbf{V})\\|_{F}^{2}$
(8)		$\displaystyle\mathrm{s.t.}$	$\displaystyle\mathrm{rank}(\tilde{\mathbf{V}})\leq M,$

where $\Omega$ denotes the set of indices corresponding to the upper-triangular part of the matrix, and $P_{\Omega}$ is an restriction operation that returns only the elements in the set $\Omega$ . Eq. (8) is a matrix completion problem [10].

Once the low-rank solution $\tilde{\mathbf{V}}=\mathbf{L}\mathbf{R}^{T}$ is obtained, where $\mathbf{L}$ and $\mathbf{R}$ are both $N\times M$ matrices, we let $\mathbf{l}_{i}$ , $\mathbf{r}_{i}$ denote the transpose of the $i$ -th row of $\mathbf{L}$ and $\mathbf{R}$ respectively. Then $\tilde{\mathbf{V}}(i,j)=\mathbf{l}_{i}^{T}\mathbf{r}_{j}$ . Therefore we can construct an MPO with the following tensor at site $i$ ( $2\leq i\leq N-1$ ):

(9)

\left(\begin{array}[]{ccc}\mathbf{I}&0&0\\ \mathbf{r}_{i}\otimes\mathbf{n}_{i}&\mathbf{I}\otimes\mathbf{I}&0\\ 0&\mathbf{l}_{i}^{T}\otimes\mathbf{n}_{i}&\mathbf{I}\end{array}\right).

Again the core tensor for $i=1$ and $i=N$ are given by the last row and the first column of Eq. (9). In order to recover the MPO construction in the exponential fitting method (6), we may define

\mathbf{l}_{i}=e^{-(N-i)\mathbf{\Lambda}}\mathbf{a},\quad\mathbf{r}_{i}=e^{(N-i+1)\mathbf{\Lambda}}\mathbf{b}.

The extra term $e^{N\mathbf{\Lambda}}$ is introduced to keep the notation consistent with a more general case discussed later.

3.1 An online matrix completion algorithm

The matrix completion problem (9) has several features:

1.

The sparsity pattern $\Omega$ is always fixed to be all the upper-triangular entries.
2.

We are actually not concerned about restoring the missing data. The missing data only serves to support a low-rank structure, and they do not necessarily need to be explicitly computed.
3.

The matrix completion problem is ill-conditioned, i.e. the ratio between the largest singular value and the smallest nonzero singular value scales exponentially with respect to the system size. Therefore explicit computation of the missing entries can easily lead to an overflow error.

The last feature requires some explanation. Consider the matrix restored by the exponential fitting procedure as in (7). The lower left corner of the matrix is $\tilde{\mathbf{V}}(N,1)=\sum_{k}a_{k}e^{\lambda_{k}(N-1)}$ . Therefore this entry, and nearby entries, grow exponentially with respect to the system size. This also results in an exponential growth of the largest singular value.

For the matrix completion problem, both convex optimization-based approaches [9] and alternating least squares (ALS)-based approaches [20, 22] require some regularization term. In our case, the regularization term will become dominant even if the regularization parameter is chosen to be very small. Online matrix completion algorithms, which applies SVD compression operations incrementally [8, 2], also suffer from the exponential growth of matrix entries in the low-rank approximation. More importantly, the largest singular value, which grows exponentially with system size, also makes the matrix completion problem difficult to solve.

We introduce a modified incremental SVD (ISVD) algorithm with missing data, which avoids the computation of the SVD of the whole matrix and therefore improves the numerical stability. This method is modified from the ISVD algorithm with missing data introduced in [2, 8]. The main modification is that we add a QR factorization step which makes the method stable even for extremely ill-conditioned matrices. We also find that the resulting MPO representation is equivalent to that in [33]. For clarity, we first introduce how to apply the original ISVD with missing data on the present matrix completion problem in this section and Section 3.2. In Section 3.3 we introduce the modified algorithm.

The algorithm proceeds in a row-by-row fashion. Define $\Omega_{p}:=\{p+1,\ldots,N\}$ the set of column indices corresponding to the sparsity pattern $\Omega$ in the $p$ -th row, and correspondingly $\Omega^{c}_{p}:=\{1,\ldots,p\}$ . Suppose we have already completed the first $p$ rows, and the completed matrix is $\tilde{\mathbf{V}}_{p}$ . This matrix is rank- $M$ and has an SVD $\tilde{\mathbf{V}}_{p}=\mathbf{U}_{p}\mathbf{S}_{p}\mathbf{W}_{p}^{T}$ . Now we want to complete the first $p+1$ rows denoted by the matrix $\mathbf{\bar{V}}_{p+1}$ . Note that the first $p$ rows of $\mathbf{\bar{V}}_{p+1}$ are just $\tilde{\mathbf{V}}_{p}$ , and on the $(p+1)$ -th row we have $\mathbf{\bar{V}}(p{+}1,\Omega_{p+1})=\mathbf{V}(p{+}1,\Omega_{p+1})$ . Our goal is to choose the unknown entries $\mathbf{\bar{V}}(p{+}1,\Omega_{p+1}^{c})$ , so that $\mathbf{\bar{V}}_{p+1}$ can also be approximated by a rank- $M$ matrix. Now assume $\mathbf{\bar{V}}(p{+}1,:)=\mathbf{x}^{T}\mathbf{W}_{p}^{T}+\mathbf{r}^{T}$ , where $\mathbf{r}^{T}\mathbf{W}_{p}=0$ . If $\mathbf{r}=0$ , then $\mathbf{\bar{V}}$ is already a rank- $M$ matrix. Otherwise

(10)

\begin{split}\mathbf{\bar{V}}_{p+1}&=\left(\begin{array}[]{c}\tilde{\mathbf{V}}_{p}\\ \mathbf{\bar{V}}(p{+}1,:)\end{array}\right)=\left(\begin{array}[]{c}\mathbf{U}_{p}\mathbf{S}_{p}\mathbf{W}_{p}^{T}\\ \mathbf{x}^{T}\mathbf{W}_{p}^{T}+\mathbf{r}^{T}\end{array}\right)\\ &=\left(\begin{array}[]{cc}\mathbf{U}_{p}\mathbf{S}_{p}&0\\ \mathbf{x}^{T}&\|\mathbf{r}\|\end{array}\right)\left(\begin{array}[]{c}\mathbf{W}_{p}^{T}\\ \mathbf{r}^{T}/\|\mathbf{r}\|\end{array}\right)\\ &=\left(\begin{array}[]{cc}\mathbf{U}_{p}&0\\ 0&1\end{array}\right)\left(\begin{array}[]{cc}\mathbf{S}_{p}&0\\ \mathbf{x}^{T}&\|\mathbf{r}\|\end{array}\right)\left(\begin{array}[]{c}\mathbf{W}_{p}^{T}\\ \mathbf{r}^{T}/\|\mathbf{r}\|\end{array}\right).\end{split}

Note that on the third line of Eq. (10), both the first and the third matrices are orthogonal matrices. Hence the smallest singular value of $\mathbf{\bar{V}}_{p+1}$ is upper bounded by $\|\mathbf{r}\|$ . In order to approximate $\mathbf{\bar{V}}_{p+1}$ by a rank- $M$ matrix, we want $\|\mathbf{r}\|$ to be as small as possible. This corresponds to the following optimization problem:

\min_{\mathbf{\bar{V}}(p{+}1,\Omega_{p+1}^{c})}\min_{\mathbf{x}}\|\mathbf{\bar{V}}(p{+}1,:)-\mathbf{x}^{T}\mathbf{W}_{p}^{T}\|^{2}.

It is easy to see that the solution is

	$\displaystyle\mathbf{x}=\mathbf{W}_{p}(:,\Omega_{p+1})^{+}\mathbf{\bar{V}}(p{+}1,\Omega_{p+1})^{T}$
	$\displaystyle\mathbf{\bar{V}}(p{+}1,\Omega_{p+1}^{c})=\mathbf{x}^{T}\mathbf{W}_{p}(:,\Omega_{p+1}^{c}),$

where $\mathbf{A}^{+}$ denotes the Moore–Penrose inverse of $\mathbf{A}$ . Then we obtain the truncated SVD

(11)

\left(\begin{array}[]{cc}\mathbf{S}_{p}&0\\ \mathbf{x}^{T}&\|\mathbf{r}\|\end{array}\right)\approx\mathbf{X}_{p+1}\mathbf{S}_{p+1}\mathbf{Y}_{p+1}^{T}

where $\mathbf{S}_{p+1}$ is an $M\times M$ matrix, $\mathbf{X}_{p+1},\mathbf{Y}_{p+1}$ are column orthogonal matrices of size $(M+1)\times M$ . Thus we may approximate $\mathbf{\bar{V}}_{p+1}$ by

\displaystyle\tilde{\mathbf{V}}_{p+1}

\displaystyle=\left[\left(\begin{array}[]{cc}\mathbf{U}_{p}&0\\ 0&1\end{array}\right)\mathbf{X}_{p+1}\right]\mathbf{S}_{p+1}\left[\mathbf{Y}_{p+1}^{T}\left(\begin{array}[]{c}\mathbf{W}_{p}^{T}\\ \mathbf{r}^{T}/\|\mathbf{r}\|\end{array}\right)\right].

Thus letting

\mathbf{U}_{p+1}=\left(\begin{array}[]{cc}\mathbf{U}_{p}&0\\ 0&1\end{array}\right)\mathbf{X}_{p+1},\quad\mathbf{W}_{p+1}=\left(\begin{array}[]{cc}\mathbf{W}_{p}&\mathbf{r}/\|\mathbf{r}\|\end{array}\right)\mathbf{Y}_{p+1},

we can proceed to the next row. The $\mathbf{U}_{N-1}$ , $\mathbf{S}_{N-1}$ and $\mathbf{W}_{N-1}$ we get in the end gives the SVD of the completed matrix (noting that $\Omega_{N}=\emptyset$ ).

3.2 Constructing the MPO

We now apply the above procedure to the coefficient matrix $\mathbf{V}$ . Let $\mathbf{l}_{i}=\mathbf{U}_{N-1}(i,:)^{T}$ and $\mathbf{r}_{i}=\mathbf{S}_{N-1}\mathbf{W}_{N-1}^{T}(:,i)$ , and use (9), we get an MPO that represents the interaction $\hat{V}$ . However, this procedure is not numerically stable. Let us consider the Coulomb interaction approximated with exponential fitting method. Then the ratio $\|\mathbf{l}_{N-1}\|/\|\mathbf{l}_{1}\|$ grows exponentially with respect to the system size, and so is $\|\mathbf{r}_{2}\|/\|\mathbf{r}_{N}\|$ .

To solve this problem, we will make use of the interior part of the MPO tensor. Write

\mathbf{X}_{p}=\left(\begin{array}[]{c}\tilde{\mathbf{X}}_{p}\\ \mathbf{a}_{p}^{T}\end{array}\right),

then we have

	$\displaystyle\mathbf{U}_{p}$	$\displaystyle=\left(\begin{array}[]{cc}\mathbf{U}_{p-1}&0\\ 0&1\end{array}\right)\left(\begin{array}[]{c}\tilde{\mathbf{X}}_{p}\\ \mathbf{a}_{p}^{T}\end{array}\right)=\left(\begin{array}[]{c}\mathbf{U}_{p-1}\tilde{\mathbf{X}}_{p}\\ \mathbf{a}_{p}^{T}\end{array}\right)$
		$\displaystyle=\left(\begin{array}[]{c}\mathbf{U}_{p-2}\tilde{\mathbf{X}}_{p-1}\tilde{\mathbf{X}}_{p}\\ \mathbf{a}_{p-1}^{T}\tilde{\mathbf{X}}_{p}\\ \mathbf{a}_{p}^{T}\end{array}\right)=\cdots$

Thus

\mathbf{l}_{i}^{T}=\mathbf{U}_{N-1}(i,:)=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i+1}\tilde{\mathbf{X}}_{i+2}\cdots\tilde{\mathbf{X}}_{N-1}

where $\mathbf{a}_{1}=\mathbf{U}_{1}=(1)$ . Therefore we have

\tilde{\mathbf{V}}(i,j)=\mathbf{l}_{i}^{T}\mathbf{r}_{j}=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i+1}\tilde{\mathbf{X}}_{i+2}\cdots\tilde{\mathbf{X}}_{N-1}\mathbf{r}_{j}.

Now define

(12)

\mathbf{b}_{j}=\tilde{\mathbf{X}}_{j}\tilde{\mathbf{X}}_{j+1}\cdots\tilde{\mathbf{X}}_{N-1}\mathbf{r}_{j},

we then have

(13)

\tilde{\mathbf{V}}(i,j)=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i+1}\tilde{\mathbf{X}}_{i+2}\cdots\tilde{\mathbf{X}}_{j-1}\mathbf{b}_{j}.

Therefore we can construct an MPO whose tensors are of the form

(14)

\left(\begin{array}[]{ccc}\mathbf{I}&0&0\\ \mathbf{b}_{i}\otimes\mathbf{n}_{i}&\tilde{\mathbf{X}}_{i}\otimes\mathbf{I}&0\\ 0&\mathbf{a}_{i}^{T}\otimes\mathbf{n}_{i}&\mathbf{I}\end{array}\right)

to represent $\sum_{i<j}\tilde{\mathbf{V}}(i,j)\hat{n}_{i}\hat{n}_{j}$ . This is analogous to (6) in the exponential fitting. The algorithm is summarized in Algorithm 1. Due to the fact that $\Omega_{p+1}\subset\Omega_{p}$ , as a practical algorithm we do not need to keep track of $\mathbf{\bar{V}}_{p}$ or $\tilde{\mathbf{V}}_{p}$ , but only $\mathbf{W}_{p}(:,\Omega_{p})$ .

Algorithm 1 Constructing the MPO representation via matrix completion

\mathbf{V}\in\mathbb{R}^{N\times N},M

\mathbf{S}_{1}\leftarrow(\|\mathbf{V}(1,2{:}N)\|)

\mathbf{W}_{1}\leftarrow(\mathbf{V}(1,2{:}N)^{T}/\|\mathbf{V}(1,2{:}N)\|)\in\mathbb{R}^{(N-1)\times 1}

\mathbf{a}_{1}\leftarrow(1)

\mathbf{b}_{2}\leftarrow(\mathbf{V}(1,2))

5: for

p=1{:}N-2

\mathbf{x}\leftarrow\mathrm{argmin}_{\mathbf{x}}\|\mathbf{W}_{p}(2{:}N{-}p,:)\mathbf{x}-\mathbf{V}(p{+}1,p{+}2{:}N)^{T}\|

\mathbf{r}\leftarrow\mathbf{V}(p{+}1,p{+}2{:}N)^{T}-\mathbf{W}_{p}(2{:}N{-}p,:)\mathbf{x}

8: if

\|r\|>\mathrm{tol}

then

\mathbf{X}_{p+1},\mathbf{S}_{p+1},\mathbf{Y}_{p+1}\leftarrow\mathrm{SVD\ of\ }\left(\begin{array}[]{cc}\mathbf{S}_{p}&0\\ \mathbf{x}^{T}&\|\mathbf{r}\|\end{array}\right)

{keep

M

singular values}

10:

\mathbf{W}_{p+1}\leftarrow\left(\begin{array}[]{cc}\mathbf{W}_{p}(2{:}N{-}p,:)&\mathbf{r}/\|\mathbf{r}\|\end{array}\right)\mathbf{Y}_{p+1}

11: else

12:

\mathbf{X}_{p+1},\mathbf{S}_{p+1},\mathbf{Y}_{p+1}\leftarrow\mathrm{SVD\ of\ }\left(\begin{array}[]{c}\mathbf{S}_{p}\\ \mathbf{x}^{T}\end{array}\right)

{keep

M

singular values}

13:

\mathbf{W}_{p+1}\leftarrow\mathbf{W}_{p}(2{:}N{-}p,:)\mathbf{Y}_{p+1}

14: end if

15:

\mathbf{b}_{p+2}\leftarrow\mathbf{S}_{p+1}\mathbf{W}_{p+1}^{T}(:,1)

16: if

p\neq N-1

then

17:

\mathbf{a}_{p+1}\leftarrow\mathbf{X}_{M+1}(p{+}1,:)^{T}

\tilde{\mathbf{X}}_{p+1}\leftarrow\mathbf{X}_{p+1}(1{:}M,:)

18: end if

19: end for

19:

\mathbf{a}_{i}

\mathbf{b}_{j}

\tilde{\mathbf{X}}_{k}

for

1\leq i\leq N-1

2\leq j\leq N

2\leq k\leq N-1

, and MPO according to (14).

Note that Eq. (12) is not the only possible choice for $\mathbf{b}_{j}$ . If the matrix $\mathbf{V}$ is exactly of rank $M$ , then

\mathbf{b}_{j}=\tilde{\mathbf{X}}_{j}\mathbf{S}_{j}\mathbf{W}_{j}^{T}(:,j)=\mathbf{S}_{j-1}\mathbf{W}_{j-1}^{T}(:,j).

This relation becomes approximately correct for the case when $\mathbf{V}$ is numerically low-rank. Therefore we may set

(15)

\mathbf{b}_{j}=\mathbf{S}_{j-1}\mathbf{W}_{j-1}^{T}(:,j).

We argue that Eq. (15) is in fact a better choice. Just like we can obtain $\mathbf{a}_{i}$ without computing $\mathbf{l}_{i}$ that grows exponentially with the system size, we want to do the same with $\mathbf{b}_{j}$ so that it does not rely on computing $\mathbf{r}_{j}$ either. And this is achieved by our new choice for the value of $\mathbf{b}_{j}$ . Also, the new $\mathbf{b}_{j}$ can be expected to give a more accurate approximation than the original one, because for an element of $\mathbf{V}$ on the $p$ -th column, the error of the approximation using this $\mathbf{b}_{j}$ only comes from the first $p$ SVDs, and does not depend on any future steps.

Furthermore, using the fact that $\Omega_{p+1}\subset\Omega_{p}$ , we do not need to explicitly fill out the missing data either, thereby avoiding computing and storing the lower triangular part of the matrix $\mathbf{\tilde{V}}_{p}$ . We can avoid computing the missing part $\mathbf{\bar{V}}(p{+}1,\Omega_{p+1}^{c})$ for each row as well. The pseudocode for the resulting algorithm is presented in Algorithm 1. Note that there is some abuse of notation involved in the algorithm, as we use $\mathbf{W}_{p}$ to denote the last $N-p$ rows of the $\mathbf{W}_{p}$ matrix we described previously.

3.3 A modified practical algorithm

Algorithm 2 A modified ISVD method for robust construction of the MPO representation

\mathbf{V}\in\mathbb{R}^{N\times N},M

\mathbf{S}_{1}\leftarrow(\|\mathbf{V}(1,2{:}N)\|)

\mathbf{W}_{1}\leftarrow(\mathbf{V}(1,2{:}N)^{T}/\|\mathbf{V}(1,2{:}N)\|)\in\mathbb{R}^{(N-1)\times 1}

\mathbf{a}_{1}\leftarrow(1)

\mathbf{b}_{2}\leftarrow(\mathbf{V}(1,2))

5: for

p=1{:}N-2

\mathbf{Q},\mathbf{R}\leftarrow

QR factorization of

\mathbf{W}_{p}(2{:}N{-}p,:)

\mathbf{x}\leftarrow\mathrm{argmin}_{\mathbf{x}}\|\mathbf{W}_{p}(2{:}N{-}p,:)\mathbf{x}-\mathbf{V}(p{+}1,p{+}2{:}N)^{T}\|

\mathbf{r}\leftarrow\mathbf{V}(p{+}1,p{+}2{:}N)^{T}-\mathbf{W}_{p}(2{:}N{-}p,:)\mathbf{x}

9: if

\|r\|>\mathrm{tol}

then

10:

\mathbf{X}_{p+1},\mathbf{S}_{p+1},\mathbf{Y}_{p+1}\leftarrow\mathrm{SVD\ of\ }\left(\begin{array}[]{cc}\mathbf{S}_{p}\mathbf{R}^{T}&0\\ \mathbf{x}^{T}\mathbf{R}^{T}&\|\mathbf{r}\|\end{array}\right)

{keep

M

singular values}

11:

\mathbf{W}_{p+1}\leftarrow\left(\begin{array}[]{cc}\mathbf{Q}&\mathbf{r}/\|\mathbf{r}\|\end{array}\right)\mathbf{Y}_{p+1}

12: else

13:

\mathbf{X}_{p+1},\mathbf{S}_{p+1},\mathbf{Y}_{p+1}\leftarrow\mathrm{SVD\ of\ }\left(\begin{array}[]{c}\mathbf{S}_{p}\mathbf{R}^{T}\\ \mathbf{x}^{T}\mathbf{R}^{T}\end{array}\right)

{keep

M

singular values}

14:

\mathbf{W}_{p+1}\leftarrow\mathbf{Q}\mathbf{Y}_{p+1}

15: end if

16:

\mathbf{b}_{p+2}\leftarrow\mathbf{S}_{p+1}\mathbf{W}_{p+1}^{T}(:,1)

17: if

p\neq N-1

then

18:

\mathbf{a}_{p+1}\leftarrow\mathbf{X}_{M+1}(p{+}1,:)^{T}

\tilde{\mathbf{X}}_{p+1}\leftarrow\mathbf{X}_{p+1}(1{:}M,:)

19: end if

20: end for

20:

\mathbf{a}_{i}

\mathbf{b}_{j}

\tilde{\mathbf{X}}_{k}

for

1\leq i\leq N-1

2\leq j\leq N

2\leq k\leq N-1

, and MPO according to (14).

Despite that Algorithm 1 computes $\mathbf{a},\mathbf{b},\tilde{\mathbf{X}}$ directly, in practice it can still become unstable when $N$ becomes large. This is because the singular values of $\tilde{\mathbf{V}}$ grows exponentially with respect to $N$ . Then in Algorithm 1, the SVD step (11) also increasingly ill-conditioned as $p$ increases.

We now further modify the ISVD algorithm as follows. Instead of trying to approximate the SVD of $\mathbf{V}(1{:}p,:)$ for each $p$ as is done in the previous sections, we compute the approximate SVD for each off-diagonal block $\mathbf{V}(1{:}p,\Omega_{p})=\mathbf{V}(1{:}p,p{+}1{:}N)$ , which avoids the exponential growth of singular values. This can be implemented by introducing only an extra QR factorization step.

Denoting the approximate SVD obtained for $\mathbf{V}(1{:}p,p{+}1{:}N)$ as $\mathbf{U}_{p}\mathbf{S}_{p}\mathbf{W}_{p}^{T}$ (note this $\mathbf{W}_{p}$ is different from the one defined previously), we perform a QR factorization on $\mathbf{W}_{p}(2{:}N{-}p,:)=\mathbf{QR}$ . Then we have

\left(\begin{array}[]{c}\mathbf{U}_{p}\mathbf{S}_{p}[\mathbf{W}_{p}(2{:}N{-}p,:)]^{T}\\ \mathbf{V}(p{+}1,p{+}2{:}N)\end{array}\right)=\left(\begin{array}[]{cc}\mathbf{U}_{p}&0\\ 0&1\end{array}\right)\left(\begin{array}[]{cc}\mathbf{S}_{p}\mathbf{R}^{T}&0\\ \mathbf{x}^{T}\mathbf{R}^{T}&\|\mathbf{r}\|\end{array}\right)\left(\begin{array}[]{c}\mathbf{Q}^{T}\\ \mathbf{r}^{T}/\|\mathbf{r}\|\end{array}\right)

Then we perform SVD on the matrix in the middle of the right-hand side, and use this to obtain $\mathbf{U}_{p+1}$ , $\mathbf{S}_{p+1}$ , $\mathbf{W}_{p+1}$ . The pseudocode for the algorithm is described in Algorithm 2. The extra QR factorization step is in Line 6 in Alg. 2 We remark that the MPO representation obtained from Algorithm 2 is equivalent to that obtained in [33].

Below we demonstrate that one can reconstruct the UTLR low-rank factors from the $\mathbf{a}_{i}$ , $\tilde{\mathbf{X}}_{i}$ , $\mathbf{b}_{i}$ obtained from Algorithm 2. First, note that not all $\tilde{\mathbf{X}}_{i}$ ’s are square matrices, though all $\tilde{\mathbf{X}}_{i}$ ’s are full-rank matrices. Each $\tilde{\mathbf{X}}_{i}$ for $i=2,3,\cdots,M$ is an $(i-1)\times i$ matrix, and each $\tilde{\mathbf{X}}_{i}$ for $i=N-M+1,N-M+2,\cdots,N-1$ is an $(N-i+1)\times(N-i)$ matrix. Each $\tilde{\mathbf{X}}_{i}$ ’s in the middle is an $M\times M$ matrix. Therefore the first $M-1$ $\tilde{\mathbf{X}}_{i}$ ’s have right inverses, the last $M-1$ $\tilde{\mathbf{X}}_{i}$ ’s have left inverses, and every $\tilde{\mathbf{X}}_{i}$ in the middle has an inverse. We denote these left or right inverses as $\tilde{\mathbf{X}}_{i}^{+}$

Now we define

	$\displaystyle\mathbf{l}_{i}^{T}$	$\displaystyle=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i+1}\tilde{\mathbf{X}}_{i+2}\cdots\tilde{\mathbf{X}}_{N-M},\quad i=1,2,\cdots,N-M$
	$\displaystyle\mathbf{l}_{i}^{T}$	$\displaystyle=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i}^{+}\tilde{\mathbf{X}}_{i-1}^{+}\cdots\tilde{\mathbf{X}}_{N-M+1}^{+}\quad i=N-M+1,N-M+2,\cdots,N-1$
	$\displaystyle\mathbf{r}_{j}$	$\displaystyle=\tilde{\mathbf{X}}^{+}_{N-M}\tilde{\mathbf{X}}^{+}_{N-M-1}\cdots\tilde{\mathbf{X}}^{+}_{j}\mathbf{b}_{j}\quad j=2,3,\cdots,N-M$
	$\displaystyle\mathbf{r}_{j}$	$\displaystyle=\tilde{\mathbf{X}}_{N-M+1}\tilde{\mathbf{X}}_{N-M+2}\cdots\tilde{\mathbf{X}}_{j-1}\mathbf{b}_{j}\quad j=N-M+1,N-M+2,\cdots,N$

and it can be checked that $\mathbf{l}_{i}^{T}\mathbf{r}_{j}=\mathbf{a}_{i}^{T}\tilde{\mathbf{X}}_{i+1}\cdots\tilde{\mathbf{X}}_{j-1}\mathbf{b}_{j}$ for any $1\leq i<j\leq N$ . Therefore we have a upper-triangular low-rank approximation for the matrix $\mathbf{V}$ .

4 Representing long-range interaction via the hierarchical low-rank format

In this section we present an alternative method to construct tensor network operators with long-range pairwise interactions using the hierarchical low-rank format. The advantage of this method is that it can be naturally generalized to higher dimensional tensors represented by PEPOs. Compared to the CF-PEPO approach [24] which relies on a fitting procedure and therefore is difficult to establish the a priori error bound, we demonstrate that our algorithm yields a sum of $\mathcal{O}(\log(N)\log(N/\epsilon))$ PEPOs whose bond dimension is bounded by a constant for the Coulomb interaction.

We study the case when $\mathbf{V}$ can be represented as a hierarchical off-diagonal low-rank (HODLR) matrix and a $\mathcal{H}$ -matrix respectively. Furthermore, the hierarchical low-rank representation is not restricted to translation-invariant operators.

4.1 Rank-one MPOs and PEPOs

We first introduce some terminologies with respect to the MPOs and PEPOs we are going to use in this section. There are mainly two types of operators we are going to use extensively, which can be efficiently represented by MPOs and PEPOs. Consider a rank-one operator in 1D

\sum_{i<j}\hat{a}_{i}\hat{b}_{j},

where $\hat{a}_{i}$ and $\hat{b}_{j}$ are local operators defined on sites $i$ and $j$ , and as a result they commute with each other. This operator can be represented by an MPO with bond dimension 3, since we can treat it as a special case of the interaction in (2.1) in which $\lambda=0$ .

The notion of rank-one operators can be easily generalized to a 2D system, for operators in the form

\sum_{\mathbf{i}\prec\mathbf{j}}\hat{a}_{\mathbf{i}}\hat{b}_{\mathbf{j}},

where $\prec$ is some order assigned to the system. This can be represented by a snake-shaped MPO [24] with bond dimension 3. The coefficient matrix is a UTLR matrix of rank 1. We add some bonds with bond dimension 1 to make it a PEPO. The bond dimension of the PEPO is still 3. A graphical illustration of the PEPO is provided in Fig. 5.

Sometimes the operator acts trivially everywhere except for a given region (described by an index set $\mathcal{I}$ ) in the system, we may order the sites in such a way that all sites in $\mathcal{I}$ are labeled before the rest of the sites. Then the operator is $\hat{O}=\hat{O}_{\mathcal{I}}\otimes\hat{I}_{\mathcal{I}^{c}}$ . In such cases we call the set $\mathcal{I}$ the support of operator $\hat{O}$ . Now suppose $\hat{O}_{1},\hat{O}_{2},\cdots,\hat{O}_{m}$ are $m$ rank-one operators in a 1D system. Integer intervals $\mathcal{I}_{1},\mathcal{I}_{2},\cdots,\mathcal{I}_{m}$ are the supports of each operator respectively. If $\mathcal{I}_{i}\cap\mathcal{I}_{j}=\emptyset$ , then the operators $\hat{O}_{i},\hat{O}_{j}$ are non-overlapping. The sum of these operators can be expressed as an MPO with bond dimension 5. We give the MPO rules in Appendix A.

This notion of two operators being non-overlapping can be easily extended to 2D systems. Instead of intervals the supports are replaced by boxes. The sum of these operators can be expressed as a PEPO with bond dimension 5. We will also give the PEPO rules in Appendix B. As mentioned before, in 2D we need to define an order $\prec$ for the rank-one operators. However when we consider a sum of non-overlapping rank-one operators, the order only needs to be defined locally within each support, i.e. we do not need to impose a global ordering valid for all rank-one operators.

4.2 Hierarchical low-rank matrix

In this section we introduce the hierarchical low-rank matrix, and its use to construct MPO and PEPO representation of tensor network operators with long-range interactions. Consider the 1D system first. We denote by $\mathcal{I}_{0}=\{0,1,\cdots,N-1\}$ the set of all sites in the system. Note that we use zero-based indexing for sites in the system from now on. We divide the domain hierarchically. At the first level, we divide the system into two equal parts $\mathcal{I}_{1;0}=\{0,1,\cdots,N/2-1\}$ and $\mathcal{I}_{1;1}=\{N/2,N/2+1,\cdots,N-1\}$ . At the second level there will be four such intervals, and third level eight intervals, and so on. In general we define $\mathcal{I}_{\ell;\alpha}=\{\alpha 2^{N-\ell},\alpha 2^{N-\ell}+1,\cdots,(\alpha+1)2^{N-\ell}-1\}$ , where $0\leq\alpha\leq 2^{\ell}-1$ . There are in total $L=\log_{2}(N)$ levels. Each interval has a neighbor list $\mathrm{NL}(\mathcal{I})$ containing all the other intervals on the same level that are neighbors of $\mathcal{I}$ .

The intervals defined in this way give rise to a tree structure. When an interval $\mathcal{I}_{\ell;\alpha}\subset\mathcal{I}_{\ell-1;\beta}$ , then we say $\mathcal{I}_{\ell;\alpha}$ is a child of $\mathcal{I}_{\ell-1;\beta}$ , or equivalently $\mathcal{I}_{\ell-1;\beta}$ is a parent of $\mathcal{I}_{\ell;\alpha}$ . When two intervals are children of the same parent we call them siblings.

A key component of the definition of the hierarchical low-rank matrix is the interaction list, which we denote by $\mathrm{IL}(\mathcal{I})$ for each interval $\mathcal{I}$ . If $\mathbf{A}$ is a hierarchical low-rank matrix, and $\mathcal{J}\in\mathrm{IL}(\mathcal{I})$ , then $\mathbf{A}(\mathcal{I},\mathcal{J})$ can be approximated by a low-rank matrix. The two hierarchical low-rank matrix formats considered in this paper, the HODLR and the $\mathcal{H}$ -matrix formats, differ only in the choice of interaction lists. Given a set of interaction lists, we are ready to define the hierarchical low-rank matrix:

Definition 2 (Hierarchical low-rank matrix).

A symmetric matrix $\mathbf{A}\in\mathbb{R}^{N\times N}$ is a hierarchical low-rank matrix, if for any index set $\mathcal{I}$ and any index set $\mathcal{J}\in\mathrm{IL}(\mathcal{I})$ , there exists a rank $r\lesssim\log(N/\epsilon)$ matrix $\tilde{\mathbf{A}}_{\mathcal{I},\mathcal{J}}$ such that

\lVert\mathbf{A}(\mathcal{I},\mathcal{J})-\tilde{\mathbf{A}}_{\mathcal{I},\mathcal{J}}\rVert_{2}\leq\epsilon\lVert\tilde{\mathbf{A}}(\mathcal{I},\mathcal{J})\rVert_{2}.

From Definition 2, we immediately recognize that a UTLR matrix must be a hierarchical low-rank matrix. To see this, we simply pick $\tilde{\mathbf{A}}_{\mathcal{I},\mathcal{J}}$ to be the $(\mathcal{I},\mathcal{J})$ -th matrix block of $\tilde{\mathbf{A}}$ given by the UTLR format. The reverse statement is not true, since the matrix blocks $\tilde{\mathbf{A}}_{\mathcal{I},\mathcal{J}}$ in a hierarchical low-rank matrix may not be related at all to each other.

The concept above can be naturally extends to two-dimensional and higher dimensional systems. Instead of intervals, the interaction lists are defined using boxes. We will first use hierarchical low-rank structure in 1D to build MPOs, and then discuss its generalization to 2D for constructing PEPOs.

4.3 HODLR format

A matrix $\mathbf{A}$ of the HODLR format can be constructed as follows. For an interval $\mathcal{I}$ on level $\ell>1$ , we define the interaction list $\mathrm{IL}(\mathcal{I})$ to be the set of all $\mathcal{J}$ such that $\mathcal{J}$ is a sibling of $\mathcal{I}$ . For level 1 the two intervals have each other as the only member in their respective interaction list. Each matrix block $\mathbf{A}(\mathcal{I},\mathcal{J})$ for $\mathcal{J}\in\mathrm{IL}(\mathcal{I})$ then has a low-rank structure and therefore admits a low-rank approximation. Using this interaction list we can divide the whole matrix into non-overlapping blocks each with a low-rank structure. The division is shown in Fig. 6 (a).

Now we construct the MPOs from this format. At level $\ell$ , define

(16)

\hat{V}^{(\ell)}_{\alpha\beta}=\sum_{i\in\mathcal{I}_{\ell;\alpha}}\sum_{j\in\mathcal{I}_{\ell;\beta}}\mathbf{V}(i,j)\hat{n}_{i}\hat{n}_{j},

for $\mathcal{I}_{\ell;\beta}\in\mathrm{IL}(\mathcal{I}_{\ell;\alpha})$ . The sum of all the interactions at level $\ell$ can be written as

\hat{V}^{(\ell)}=\sum_{\alpha=0}^{2^{\ell-1}}\hat{V}_{2\alpha,2\alpha+1}^{(\ell)}.

Note only interactions allowed in the interaction list is covered on this level, which is an important fact that enables us to construct MPOs with a small bond dimension. Then we can decompose the tensor as

\hat{V}=\sum_{\ell=1}^{L}\hat{V}^{(\ell)}.

This sum covers all the interactions in $\hat{V}$ .

Because of the hierarchical low-rank structure, we can yield some approximate low-rank decomposition $\mathbf{V}(\mathcal{I}_{\ell;2\alpha},\mathcal{I}_{\ell;2\alpha+1})\approx\mathbf{A}\mathbf{B}^{T}$ . Then

\sum_{i\in\mathcal{I}_{\ell;2\alpha}}\sum_{j\in\mathcal{I}_{\ell;2\alpha+1}}\mathbf{V}(i,j)\hat{n}_{i}\hat{n}_{j}\approx\sum_{r=1}^{M}\sum_{i\in\mathcal{I}_{\ell;2\alpha}}\sum_{j\in\mathcal{I}_{\ell;2\alpha+1}}(\mathbf{A}(i,r)\hat{n}_{i})(\mathbf{B}(j,r)\hat{n}_{j}).

Note that on the right-hand side each $r$ indexes a rank-one operator. The support of each rank-one operator is $\mathcal{I}_{\ell;2\alpha}\cup\mathcal{I}_{\ell;2\alpha+1}=\mathcal{I}_{\ell-1;\alpha}$ .

Therefore, for each $\alpha=0,1,\cdots,2^{\ell-1}-1$ , we have $M$ rank-one operators. For each $\alpha$ , we pick one from these $M$ operators correspond to a given index $r$ , and sum them up to get,

\hat{V}^{(\ell)}_{r}=\sum_{\alpha}\sum_{i\in\mathcal{I}_{\ell;2\alpha}}\sum_{j\in\mathcal{I}_{\ell;2\alpha+1}}(\mathbf{A}(i,r)\hat{n}_{i})(\mathbf{B}(j,r)\hat{n}_{j}),

and then we have

\hat{V}^{(\ell)}=\sum_{r=1}^{M}\hat{V}^{(\ell)}_{r}.

Hence all the interactions at level $\ell$ can be collected into $M$ operators $\hat{V}^{(\ell)}_{r}$ , each one is a sum of non-overlapping rank-one operators. From this procedure, we have collected all interactions at level $\ell$ into a sum of $M$ MPOs each with bond dimension 5. $M$ can be chosen to be $\mathcal{O}(\log(N/\epsilon))$ for the coefficient matrix to have a 2-norm error of $\epsilon$ . The system contains $L$ levels so in total we need $ML=\mathcal{O}(\log(N)\log(N/\epsilon))$ MPOs, each with a bounded rank.

4.4 $\mathcal{H}$ -matrix format

In the HODLR matrix, a matrix block $\mathbf{A}(\mathcal{I},\mathcal{J})$ is low-rank as long as $\mathcal{I},\mathcal{J}$ are non-overlapping and are at the same level. This is also called the weak-admissibility condition [18]. The $\mathcal{H}$ -matrix is a more general class of hierarchical low-rank matrices, which requires the strong-admissibility condition, i.e. $\mathcal{I},\mathcal{J}$ are not only non-overlapping, but their distance should be comparable to their sizes, i.e. $\mathrm{dist}(\mathcal{I},\mathcal{J})\gtrsim\max\{|\mathcal{I}|,|\mathcal{J}|\}$ .

In the $\mathcal{H}$ -matrix format, we define the interaction list $\mathrm{IL}(\mathcal{I})$ of an interval $\mathcal{I}$ in the following way: an interval $\mathcal{J}\in\mathrm{\mathcal{I}}$ if and only if the parent of $\mathcal{J}$ and the parent of $\mathcal{I}$ are neighbors, and $\mathcal{J}$ is not a neighbor of $\mathcal{I}$ . Therefore in 1D it is guaranteed that $\mathrm{dist}(\mathcal{I},\mathcal{J})\geq\max\{|\mathcal{I}|,|\mathcal{J}|\}$ if $\mathcal{J}\in\mathrm{IL}(\mathcal{I})$ . A hierarchical low-rank matrix with interaction list defined in this way is called an $\mathcal{H}$ -matrix. A graphical illustration is shown in Fig. 6 (b).

We first consider a 1D system. The system is divided into intervals hierarchically the same way as in Section 4.3. Sill we divide the interactions into different levels

(17)

\hat{V}=\sum_{\ell=2}^{L}\hat{V}^{(\ell)}+\hat{D},

where $\hat{D}$ contains only nearest-neighbor interactions. Note that the lowest level is $2$ in this format, as the interaction list at level $1$ is empty. Based on the interaction list of the $\mathcal{H}$ -matrix, $\hat{V}^{(\ell)}$ is defined as

\hat{V}^{(\ell)}=\sum_{\alpha=0}^{2^{\ell-1}-2}\left(\hat{V}^{(\ell)}_{2\alpha,2\alpha+2}+\hat{V}^{(\ell)}_{2\alpha,2\alpha+3}+\hat{V}^{(\ell)}_{2\alpha+1,2\alpha+3}\right),

where $\hat{V}^{(\ell)}_{\alpha\beta}$ is defined in Eq. (16). One can see how this corresponds to the $\mathcal{H}$ -matrix interaction list from Fig. 6 (b).

Now we want to efficiently represent the operator $\hat{V}^{(\ell)}$ as linear combination of MPOs. Similar to the previous section we first write

\hat{V}^{(\ell)}_{\alpha\beta}\approx\sum_{r=1}^{M}\hat{V}^{(\ell)}_{\alpha\beta,r},\quad\hat{V}^{(\ell)}_{\alpha\beta,r}=\sum_{i\in\mathcal{I}_{\ell;\alpha}}\sum_{j\in\mathcal{I}_{\ell;\beta}}(\mathbf{A}(i,r)\hat{n}_{i})(\mathbf{B}(j,r)\hat{n}_{j}),

where $\hat{V}^{(\ell)}_{\alpha\beta,r}$ is a rank-one operator defined in (4.1). The low-rank decomposition is done by performing truncated SVD and cutting off singular values smaller than $\epsilon^{\prime}\|\mathbf{V}(\mathcal{I}_{\ell;\alpha},\mathcal{I}_{\ell;\beta})\|_{2}$ for some $\epsilon^{\prime}>0$ . Then

\hat{V}^{(\ell)}=\sum_{r=1}^{M}\sum_{\alpha=0}^{2^{\ell-1}-2}\left(\hat{V}^{(\ell)}_{2\alpha,2\alpha+2,r}+\hat{V}^{(\ell)}_{2\alpha,2\alpha+3,r}+\hat{V}^{(\ell)}_{2\alpha+1,2\alpha+3,r}\right).

However, unlike in the HODLR case, these operators can be overlapping for different $\alpha,\beta$ . For example $\hat{V}^{(\ell)}_{2\alpha,2\alpha+2,r}$ and $\hat{V}^{(\ell)}_{2\alpha,2\alpha+3,r}$ are overlapping. Also $\hat{V}^{(\ell)}_{2\alpha,2\alpha+2,r}$ and $\hat{V}^{(\ell)}_{2\alpha+2,2\alpha+4,r}$ are overlapping as well. Therefore we need to further separate the terms as

(18)		$\displaystyle\hat{V}^{(\ell)}$	$\displaystyle=\sum_{r=1}^{M}\left(\sum_{\alpha=0}^{2^{\ell-2}-1}\hat{V}^{(\ell)}_{4\alpha,4\alpha+2,r}+\sum_{\alpha=0}^{2^{\ell-2}-1}\hat{V}^{(\ell)}_{4\alpha,4\alpha+3,r}+\sum_{\alpha=0}^{2^{\ell-2}-1}\hat{V}^{(\ell)}_{4\alpha+1,4\alpha+3,r}\right)$
(19)			$\displaystyle+\sum_{r=1}^{M}\left(\sum_{\alpha=0}^{2^{\ell-2}-2}\hat{V}^{(\ell)}_{4\alpha+2,4\alpha+4,r}+\sum_{\alpha=0}^{2^{\ell-2}-2}\hat{V}^{(\ell)}_{4\alpha+2,4\alpha+5,r}+\sum_{\alpha=0}^{2^{\ell-2}-2}\hat{V}^{(\ell)}_{4\alpha+3,4\alpha+5,r}\right)$

For each of the six terms on the right-hand side we construct $M$ MPOs. Each term represents a sum of non-overlapping rank-one operators, and therefore has bond dimension 5. Operators on the first line we call type-0 and on the second line we call type-1. They can be visualized on the coefficient matrix, as shown in Fig. 6.

From the discussion we construct $6M$ MPOs each with bond dimension 5. Since the system contains $L$ levels in total we have $6LM$ MPOs. To ensure the 2-norm error of the coefficient matrix to be smaller than $\epsilon$ , we need to choose $\epsilon^{\prime}=C\epsilon/N$ for some constant $C$ . Then we have $M=\mathcal{O}(\log(1/\epsilon^{\prime}))=\mathcal{O}(\log(N/\epsilon))$ . The total number of MPOs is $\mathcal{O}(\log(N)\log(N/\epsilon))$ .

In practical calculation, we may further reduce the constant prefactor from $6$ to $4$ without increasing $M$ . This is because the strong admissibility condition allows one to concatenate certain blocks together only with a mild increase of the numerical rank [18]. For example, matrices $\mathbf{V}(\mathcal{I}_{\ell;4\alpha},\mathcal{I}_{\ell;4\alpha+2})$ and $\mathbf{V}(\mathcal{I}_{\ell;4\alpha},\mathcal{I}_{\ell;4\alpha+3})$ can be concatenated into a single matrix $\mathbf{V}(\mathcal{I}_{\ell;4\alpha},\mathcal{I}_{\ell;4\alpha+2}\cup\mathcal{I}_{\ell;4\alpha+3})$ . This allows us to combine some of the operators and thereby reduce the prefactor. The reduction of the prefactor becomes more significant in 2D and higher-dimensional cases.

4.5 Higher dimensional cases

When the coefficient matrix $\mathbf{V}$ satisfies the hierarchical low-rank format, the tensor $\hat{V}$ can be efficiently represented by PEPOs when the lattice dimension is larger than $1$ . We consider the 2D case below, and the 3D case follows similarly. We consider an $N\times N:=2^{L}\times 2^{L}$ lattice, with each site indexed by $\mathbf{i}=(\alpha_{x},\alpha_{y})\in\{1,2,\cdots,N\}^{2}$ .

In the 2D case, the intervals are replaced by boxes. We define a box on level $\ell$ , indexed by $\bm{\alpha}=(\alpha_{x},\alpha_{y})$ , by

	$\displaystyle\mathcal{I}_{\ell;\bm{\alpha}}$	$\displaystyle=\{\alpha_{x}2^{L-\ell},\alpha_{x}2^{L-\ell}+1,\cdots,(\alpha_{x}+1)2^{L-\ell}-1\}$
		$\displaystyle\times\{\alpha_{y}2^{L-\ell},\alpha_{y}2^{L-\ell}+1,\cdots,(\alpha_{y}+1)2^{L-\ell}-1\},$

and interaction between two level- $\ell$ boxes $\bm{\alpha}$ , $\bm{\beta}$ is denoted as $\hat{V}^{(\ell)}_{\bm{\alpha}\bm{\beta}}$ , similar to (16).

Similar to $\mathcal{H}$ -matrices in 1D, the interaction list $\mathrm{IL}(\mathcal{I})$ is the set of all boxes $\mathcal{J}$ such that

1.

The parent of $\mathcal{J}$ is a neighbor of the parent of $\mathcal{I}$ ,
2.

$\mathcal{I}$ and $\mathcal{J}$ are not neighbors.

Fig. 7 (a) shows an example of the allowed interactions at level $4$ . All the interactions at level $\ell$ are collected into an operator $\hat{V}^{(\ell)}$ , and we have the same decomposition as in (17).

In the previous section we have divided the interactions at level $\ell$ into type-0 and type-1, and further into six terms in (18) and (19) respectively. Here we divide all interactions into four types: type- $(p,q)$ for $p,q\in\{0,1\}$ . And we write

\hat{V}^{(\ell)}=\sum_{pq}\sum_{\alpha_{x}=0}^{2^{\ell-2}-1-p}\sum_{\alpha_{y}=0}^{2^{\ell-2}-1-q}\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})},

where $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}$ covers most, but not all interactions at level $\ell$ (to avoid double counting) between boxes $\mathcal{I}_{\ell;\bm{\gamma}}$ for $\bm{\gamma}=(4\alpha_{x}+2p+\beta_{x},4\alpha_{y}+2q+\beta_{y})$ , $\beta_{x},\beta_{y}\in\{0,1,2,3\}$ . The boxes can be divided into a $4\times 4$ array. Some examples are given in Fig. 7 (b). In general, $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}$ with the same $p,q$ but different $\alpha_{x},\alpha_{y}$ are non-overlapping, and therefore can be summed up without increasing the bond dimension beyond 5. A graphical description of the rules is given in Fig. 7 (c).

Using the simplest scheme each $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}$ contains at most 78 interactions between pairs of boxes at level $\ell$ . The number 78 is obtained as follows: there are $16\times 15/2=120$ pairwise interactions between boxes in $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}$ , as can be seen in Fig. 7 (b). 42 of these pairwise interactions are between adjacent boxes (sharing an edge or vertex), and therefore in total there are 78 permissible ones. Note that there will be double counting if we include all 78 interactions in each $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}$ , so in fact this maximum is only attained at the lower-right corner according to the assignment rule in Fig. 7 (c).

However, using the fact that for $\mathcal{I}_{\ell;\bm{\beta}_{1}},\mathcal{I}_{\ell;\bm{\beta}_{2}}\in\mathrm{IL}(\mathcal{I}_{\ell;\bm{\alpha}})$ , we may directly compress the interactions between $\mathcal{I}_{\ell;\bm{\alpha}}$ and $\mathcal{I}_{\ell;\bm{\beta}_{1}}\cup\mathcal{I}_{\ell;\bm{\beta}_{2}}$ . This is the same as what we have done for the 1D case in the previous section. Using this method, all pairwise interactions involving one box in the $4\times 4$ box array can be merged into one. Any interaction included has to involve one of the 12 boxes that are not among the 4 boxes at the lower-right corner. The result is that we can reduce the number of interactions to 12 for $p=q=0$ and 8 for all other cases.

Each of the 12 or 8 interactions, which we index by $k$ , is then approximated by the linear combination of rank-one operators. Therefore

\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y})}=\sum_{k=1}^{K_{pq}}\sum_{r=1}^{M}\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y}),kr},

where $K_{pq}=12$ for $p=q=0$ and 8 otherwise, and each $\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y}),kr}$ is a rank-one operator. We simply add zero operators when there is not enough $k$ . It is clear that this is not the most efficient method, but it gives us the desired scaling in the end. The rank-one operators are obtained using a low-rank decomposition similar to (4.4), which is in turn obtained from a truncated SVD.

Then we can define

\hat{V}^{(l,p,q)}_{kr}=\sum_{\alpha_{x}=0}^{2^{\ell-2}-1-p}\sum_{\alpha_{y}=0}^{2^{\ell-2}-1-q}\hat{V}^{(l,p,q)}_{(\alpha_{x},\alpha_{y}),kr}.

Everything we have done above is to make this a sum of non-overlapping rank-one operator introduced in Section 4.1. It can be represented by a PEPO with bond dimension 5.

Therefore all the interactions at level $\ell$ can be written as

\hat{V}^{(\ell)}=\sum_{p=0}^{1}\sum_{q=0}^{1}\sum_{k=1}^{K_{pq}}\sum_{r=1}^{M}\hat{V}^{(l,p,q)}_{kr}.

This is a linear combination of $36M$ PEPOs with bond dimension 5.

The total number of PEPOs needed for the system will be $36ML$ . To ensure the coefficient matrix 2-norm error to be below $\epsilon$ we need $M=\mathcal{O}(\log(N/\epsilon))$ and $L=\log(N)$ . The number of PEPOs needed will scale like $\mathcal{O}(\log(N/\epsilon)\log(N))$ . We remark that this is in fact a rather loose upper bound. Our numerical experiments indicate that for small-to-medium sized systems (e.g. $N=128$ ), the actual number of PEPOs needed is usually much smaller.

5 Error analysis for the Coulomb matrix

The modified ISVD algorithm and the hierarchical low-rank representation can be applied to general pairwise interactions with or without the translation-invariance property. To illustrate the efficiency of the compressed representation, in this section we briefly list some known results for the Coulomb interaction.

For the Coulomb interaction, the error bound for the exponential fitting method using the quadrature scheme is given in [3, 6]. In particular, $0<\delta<r<D$ and any $\epsilon>0$ , there exist $M=\mathcal{O}\{\log[D/(\delta\varepsilon)]\}$ numbers $\alpha_{l},\lambda_{l}>0$ such that

(20)

\left|\frac{1}{r}-\sum_{l=1}^{M}\alpha_{l}\exp\left(-\lambda_{l}r\right)\right|<\epsilon.

A more general approximation can be obtained from classical results in the fast multipole method (e.g. [16, Theorem 3.2]). Consider the two-dimensional case for example, we may obtain the following approximate low-rank decomposition of the off-diagonal part of the Coulomb kernel: for two regions $\mathcal{A}\subset\mathbb{R}^{2}$ and $\mathcal{B}\subset\mathbb{R}^{2}$ , such that $\mathcal{A}\subset B(\mathbf{r}_{0},a)$ and $\mathcal{B}\subset B(\mathbf{r}_{0},r)^{c}$ , $a<r$ , then for any $p\geq 0$ , there exists $a_{k}\in C^{\infty}(\mathcal{A})$ and $b_{k}\in C^{\infty}(\mathcal{B})$ , $1\leq k\leq 2p+1=:M$ such that for any $\mathbf{r}\in\mathcal{A}$ , $\mathbf{r}^{\prime}\in\mathcal{B}$ ,

(21)

\left|\frac{1}{\lVert\mathbf{r}-\mathbf{r^{\prime}}\rVert}-\sum_{k=1}^{M}a_{k}(\mathbf{r})b_{k}(\mathbf{r}^{\prime})\right|\leq\frac{1}{r-a}\left(\frac{a}{r}\right)^{p+1}.

This immediately implies that in order to control the pointwise error up to $\epsilon$ , the rank $M$ can be bounded by $\log(1/\epsilon)$ . In order to control the 2-norm of each low-rank block, $M$ can be bounded by $\log(N/\epsilon)$ . Note that the separation between the sets $\mathcal{A},\mathcal{B}$ is crucial for the error estimate, which corresponds to the strong admissibility condition. For 2D and higher dimensional Coulomb interaction, the separation is indeed crucial to bound the PEPO rank. Hence we only consider the $\mathcal{H}$ -matrix format for 2D systems in the numerical examples below.

6 Numerical examples

In this section we present numerical results for 1D and 2D systems using the methods we introduced above. As before, we will only focus on long-range pairwise interactions of the form (1).

6.1 Quasi 1D system

We compare the numerical performance of the exponential fitting method and the modified ISVD described in Algorithm 2 on Coulomb interaction of the form

\hat{V}=\sum_{1\leq i<j\leq N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}\hat{n}_{i}\hat{n}_{j},

where $\mathbf{r}_{i}=(x_{i},y_{i})$ . We always choose $x_{i}=i$ , $|y_{i}|<1$ so that this mimics a quasi-1D system. Define the error matrix as $\mathbf{E}(i,j)=\tilde{\mathbf{V}}(i,j)-\mathbf{V}(i,j)$ , and we measure the error using the 2-norm of the error matrix.

In the first example we let all $y_{i}=0$ . The resulting system is translation invariant, and we can apply the exponential fitting method as well as modified ISVD. The numerical results are shown in Fig. 8.

Fig. 8 shows that the accuracy of the two methods is comparable, and the modified ISVD outperforms exponential fitting by a small margin in terms of 2-norm error. However, the distribution of the error of the two methods is very different. The error of the modified ISVD is concentrated at the upper-right corner of the matrix, which is the interaction between the two ends of the system. This can be expected from Eq. (13), where the error accumulation becomes more significant as more matrices $\tilde{\mathbf{X}}_{k}$ are multiplied together. On the other hand, the error of exponential fitting is concentrated some distance away from the diagonal.

In the second example we set $y_{i}$ to be random numbers drawn uniformly from $[-0.2,0.2]$ . system is then no longer translation-invariant, and the exponential fitting method is not applicable. However, the modified ISVD still performs well. In Fig. 9 we can see that the convergence of error with respect to $M$ and error distribution are both very similar to the translation-invariant case.

We also test the performance of HODLR and $\mathcal{H}$ -matrix based MPO construction on this model with $y_{i}=0$ . The number of MPOs needed with fixed singular value threshold $\epsilon^{\prime}=10^{-3}$ and the relative 2-norm error is plotted in Fig. 10 (a). The $\mathcal{H}$ -matrix format requires a large number of MPOs than the HODLR format. The error distribution is plotted for system size $N=128$ , in Fig. 10 (b) (c).

system size	number of PEPOs	maximal rank
$16\times 16$	328	13
$32\times 32$	756	14
$64\times 64$	1221	15
$128\times 128$	1696	15

Table 1: Performance of the

\mathcal{H}

-matrix format for the 2D system. The truncation threshold

\epsilon^{\prime}=10^{-5}

6.2 2D system

We then consider a 2D system with sites on a uniform $N\times N$ grid. For the 2D system we only show the numerical results for the $\mathcal{H}$ -matrix format, since the Coulomb interaction in 2D only satisfies the strong admissibility condition. In other words, the growth of the rank $M$ becomes unfavorable using the HODLR format. The error distribution is plotted in Fig. 10 (d). Table 1 shows the number of PEPOs, as well as the maximal PEPO rank for different system sizes in order to achieve the same truncation criterion $\epsilon^{\prime}$ (we truncate singular values for block $\mathbf{V}(\mathcal{I},\mathcal{J})$ ) at threshold $\epsilon^{\prime}\|\mathbf{V}(\mathcal{I},\mathcal{J})\|_{2}$ ). We observe that the maximal rank $M$ grows very mildly with respect to the system size. The number of PEPOs increases more noticeably when the system size is small, and the growth becomes slower as the side length $N$ increases from $64$ to $128$ . The number of PEPOs should eventually become only logarithmic with respect to $N$ .

7 Conclusion

We introduced two methods for constructing efficient MPO and PEPO representation of tensors with pairwise long-range interactions, without assuming that the system is translation-invariant. Using the upper-triangular low-rank (UTLR) structure, the construction of the MPO representation can be solved as a matrix completion problem. We developed a modified ISVD algorithm to solve the matrix completion problem with in a numerically stable fashion. Although we do not yet have a theoretical error bound for the accuracy of the modified ISVD algorithm, numerical results indicate that the error of the modified ISVD method is comparable to that of the exponential fitting method, which indicates that for Coulomb interaction, the MPO rank from the modified ISVD method can be bounded by $\log(N/\epsilon))$ . The hierarchical low-rank format can be used to construct both MPOs and PEPOs. For the Coulomb interaction, the tensor can be approximated by the linear combination of $\mathcal{O}(\log(N)\log(N/\epsilon))$ MPOs/PEPOs, and the rank of each MPO/PEPO is bounded by a constant. The hierarchical low-rank format leads to an MPO representation with a higher MPO rank than that from the modified ISVD method. However, the linear combination representation may naturally facilitate parallel computation, and this is particularly the case for the PEPO representation where a relatively large number of PEPOs can be needed. On the other hand, the MPOs/PEPOs from our hierarchical low-rank format have many zero entries. It would be of practical interest to develop alternative algorithms to reduce the number of zeros and to obtain a more compact representation with reduced preconstants.

Acknowledgments

This work was partially supported by the Department of Energy under Grant No. DE-SC0017867, No. DE-AC02-05CH11231 (L. L.), by the Air Force Office of Scientific Research under award number FA9550-18-1-0095 (L.L. and Y. T.). We thank Garnet Chan, Matthew O’Rourke and Miles Stoudenmire for helpful discussions.

Appendix A FSM rules for sum of non-overlapping MPOs

For a 1D system we divide into intervals $\{\mathcal{I}_{\alpha}\}$ such as done in our construction of the hierarchical low-rank matrix. First we clarify some notations. For an operator defined in an interval we use the notation $\hat{O}$ . For example, the identity operator defined in an interval $\{\mathcal{I}_{\alpha}\}$ is $\hat{I}_{\alpha}$ . Tensors and operators defined on a single site are still denoted as $\mathbf{T}$ , $\mathbf{O}$ , etc. We want to represent an operator of the form

\hat{H}=\sum_{\alpha}\left(\bigotimes_{\beta<\alpha}\hat{I}_{\beta}\right)\otimes\hat{H}_{\alpha}\otimes\left(\bigotimes_{\beta>\alpha}\hat{I}_{\beta}\right),

using an MPO, where we already have an MPO representation of each operator $\hat{H}_{\alpha}$ defined on interval $\mathcal{I}_{\alpha}$ .

Using the finite state machine (FSM) description of MPOs [12, 14], we assume there are the following general finite state machine rules for each $\hat{H}_{\alpha}$ shown in Table 2. It is equivalent to the equation

\hat{H}_{\alpha}=\sum_{r_{1},r_{2},\cdots,r_{\|\mathcal{I}_{\alpha}\|-1}}\bigotimes_{i=1}^{|\mathcal{I}_{\alpha}|}\mathbf{H}_{i}(r_{i-1},r_{i}),

where $r_{0}=r_{\|\mathcal{I}_{\alpha}\|}=1$ , and $\mathbf{H}(l,r)$ is an operator for each $l,r$ . The MPO rules for $\hat{H}$ are shown in Table 3.

position	rule num	input state	operator	output state
ALL	1	$l$	$\mathbf{H}_{i}(l,r)$	$r$

Table 2: MPO rules for the

\hat{H}_{\alpha}

on interval

\mathcal{I}_{\alpha}

	rule num	input state	operator	output state
HEAD	1	$\bigtriangleup$	$\mathbf{I}_{i}$	$\bigtriangleup$
	2	$\bigtriangledown$	$\mathbf{I}_{i}$	$\bigtriangledown$
	3	$\bigtriangleup$	$\mathbf{H}_{i}(1,r)$	$r$
BODY	4	$\bigtriangleup$	$\mathbf{I}_{i}$	$\bigtriangleup$
	5	$\bigtriangledown$	$\mathbf{I}_{i}$	$\bigtriangledown$
	6	$l$	$\mathbf{H}_{i}(l,r)$	$r$
TAIL	7	$\bigtriangleup$	$\mathbf{I}_{i}$	$\bigtriangleup$
	8	$\bigtriangledown$	$\mathbf{I}_{i}$	$\bigtriangledown$
	9	$l$	$\mathbf{H}_{i}(l,1)$	$\bigtriangledown$

Table 3: MPO rules for the sum of operators

\left(\bigotimes_{\beta<\alpha}\hat{I}_{\beta}\right)\otimes\hat{H}_{\alpha}\otimes\left(\bigotimes_{\beta>\alpha}\hat{I}_{\beta}\right)

For the FSM rules of the operators $\hat{H}_{\alpha}$ , we assume the FSM is in state 1 at the very beginning and will be in state 1 at the end. This means on the first site in interval $\mathcal{I}_{\alpha}$ an operator $\mathbf{H}_{i}(1,r)$ is activated and then the FSM goes to state $r$ , and on the last site an operator $\mathbf{H}_{i}(l,1)$ is activated, and the FSM goes into state 1. The first site in the interval is called the head and the last site is called the tail. Note that the FSM rules on the head and tail sites are different from on other sites inside the interval.

Now we construct the FSM rules for $\hat{H}$ . We add two states $\bigtriangleup$ and $\bigtriangledown$ into the FSM. $\bigtriangleup$ is the default state meaning a non-trivial operator, i.e. $\hat{H}_{\alpha}$ , has not been activated, and $\bigtriangledown$ means such an operator has been activated, and only identity operators should be picked afterwards. The FSM is in state $\bigtriangleup$ both at the very beginning and at the end.

At the first site, i.e. the head, of each interval $\mathcal{I}_{\alpha}$ , the input signal from the left is either $\bigtriangleup$ or $\bigtriangledown$ . If the signal is $\bigtriangleup$ , then we choose either to not activate $\hat{H}_{\alpha}$ and pass along the signal $\bigtriangleup$ , as in rule 1 of Table 3, or to activate $\hat{H}_{\alpha}$ , as in rule 3 of the same table. If the signal is $\bigtriangledown$ , then $\hat{H}_{\alpha}$ cannot be activated and a signal $\bigtriangledown$ is passed down. This is done in rule 2.

In the middle, or body, of each interval $\mathcal{I}_{\alpha}$ , if $\hat{H}_{\alpha}$ has been activated at the first site, then we follow the FSM rules for the individual $\hat{H}_{\alpha}$ (rule 6). If we receive either $\bigtriangleup$ or $\bigtriangledown$ then we pass the signal down (rules 4 and 5).

At the last site, or the tail, of each interval, if $\hat{H}_{\alpha}$ has been activated, then the output signal is $\bigtriangledown$ (rule 9). Otherwise the output signal is whatever this interval received from the beginning (rules 7 and 8).

An important fact is that we have only added two signals into the original set of signals. This means if the original MPOs $\{\hat{H}_{\alpha}\}$ have a maximum bond dimension $D$ , then the new MPO for the sum has a bond dimension $D+2$ . In the case of the sum of non-overlapping rank-one operators as discussed in 4.1, the bond dimension of the sum is 5.

Appendix B FSM rules for sum of non-overlapping PEPOs

The FSM rules for PEPOs are more complex. To describe our FSM rules for PEPOs, we first divide the sites inside a given box into different parts, as shown in the right part of Fig. 11. Different FSM rules should be applied to different parts. Sites in a box are labeled based on their position in the box, as shown in the left part of Fig. 11. For example, the box at the upper-right corner of a box is labeled UR.

We assume for each individual PEPO $\hat{H}_{\bm{\alpha}}$ , we have the FSM rules specified in Table 4. Note that not all PEPO tensors have 4 internal indices and 2 physical indices, as can be seen from Fig. 3 (b). For example, the tensors on the edges have 3 internal indices. In this case, we simply add bonds with bond dimension 1 to make them have 4 internal indices, similar to what we have done in the snake-shaped PEPO construction in 4.1. These extra indices all take the value 1.

Now we design a set of FSM rules to represent

\hat{H}=\sum_{\bm{\alpha}}\left(\bigotimes_{\bm{\beta}\prec\bm{\alpha}}\hat{I}_{\bm{\beta}}\right)\otimes\hat{H}_{\bm{\alpha}}\otimes\left(\bigotimes_{\mathbf{\bm{\beta}}\succ\bm{\alpha}}\hat{I}_{\bm{\beta}}\right),

where $\prec$ gives the lexicographical order defined in Section 1. One key difference between 2D and 1D is that in 2D we no longer have a fixed direction for the signal of FSM to be passed along. Therefore we require that signals between boxes (there are signals passing within boxes according to the original FSM rules in Table 4, which we leave unchanged) only pass rightward or downward, making it easier to design FSM rules.

position	rule num	$(l,u,d,r)$	operator
ALL	1	$(l,u,d,r)$	$\mathbf{H}_{i}(l,u,d,r)$

Table 4: PEPO rules for

\hat{H}_{\bm{\alpha}}

in the box

\mathcal{I}_{\bm{\alpha}}

We first design the FSM rules for all parts of the system except for the upper-right corner, the right edge, and the lower-right corner (the colored parts in Fig. 11 (b)). The rules are shown in Table 5. Just like in the 1D case, we introduce two signals $\bigtriangleup$ and $\bigtriangledown$ . $\bigtriangleup$ is the default signal. We change the value of the extra indices we added at the edges and corners from 1 to $\bigtriangleup$ .

At the UL site of each box $\mathcal{I}_{\bm{\alpha}}$ , there are two input signals from left and up. When both are $\bigtriangleup$ , then it means $\hat{H}_{\bm{\alpha}}$ can be activated in this box (rule 3), but we can also choose not to activate it and pass along the signal (rule 1). If the site receives a signal $\bigtriangledown$ from the left, and $\bigtriangleup$ from above, then it means an $\hat{H}_{\bm{\alpha}^{\prime}}$ has already been activated for some $\bm{\alpha}^{\prime}$ to the left, and the only thing that can be done now is to pass down the signal $\bigtriangledown$ to the right until the right edge (rule 2), where other FSM rules will process it.

	rule num	$(l,u,d,r)$	operator
UL	1	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	2	$(\bigtriangledown,\bigtriangleup,\bigtriangleup,\bigtriangledown)$	$\mathbf{I}_{i}$
	3	$(\bigtriangleup,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(1,1,d,r)$
U	4	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	5	$(\bigtriangledown,\bigtriangleup,\bigtriangleup,\bigtriangledown)$	$\mathbf{I}_{i}$
	6	$(l,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(l,1,d,r)$
UR	7	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	8	$(\bigtriangledown,\bigtriangleup,\bigtriangleup,\bigtriangledown)$	$\mathbf{I}_{i}$
	9	$(l,\bigtriangleup,d,\bigtriangledown)$	$\mathbf{H}_{i}(l,1,d,1)$
R	10	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
R	11	$(l,u,d,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,d,1)$
DR	12	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
DR	13	$(l,u,\bigtriangleup,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,1,1)$
D	14	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
D	15	$(l,u,\bigtriangleup,r)$	$\mathbf{H}_{i}(l,u,1,r)$
DL	16	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
DL	17	$(\bigtriangleup,u,\bigtriangleup,r)$	$\mathbf{H}_{i}(1,u,1,r)$
L	18	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
L	19	$(\bigtriangleup,u,d,r)$	$\mathbf{H}_{i}(1,u,d,r)$
BODY	20	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
BODY	21	$(l,u,d,r)$	$\mathbf{H}_{i}(l,u,d,r)$

Table 5: PEPO rules for the sum of operators

\left(\bigotimes_{\bm{\beta}\prec\bm{\alpha}}\hat{I}_{\bm{\beta}}\right)\otimes\hat{H}_{\bm{\alpha}}\otimes\left(\bigotimes_{\mathbf{\bm{\beta}}\succ\bm{\alpha}}\hat{I}_{\bm{\beta}}\right)

in the whole system except for the upper-right corner, the right edge, and the lower-right corner.

Other rules in Table 5 are either identical to the FSM rules for $\hat{H}_{\bm{\alpha}}$ , or to pass along signals $\bigtriangleup$ and $\bigtriangledown$ . Note that the signal $\bigtriangledown$ is only passed along the upper edge (UL, U, and UR in Fig. 11 (a)) as specified in rules 2, 5, 8. If $\hat{H}_{\bm{\alpha}}$ has been activated then a signal $\bigtriangledown$ will be passed to the box to the right from the UR site (rule 9).

	rule num	$(l,u,d,r)$	operator
UL	1	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	2	$(\bigtriangledown,\bigtriangleup,\bigtriangledown,\bigtriangleup)$	$\mathbf{I}_{i}$
	3	$(\bigtriangleup,\bigtriangledown,\bigtriangledown,\bigtriangleup)$	$\mathbf{I}_{i}$
	4	$(\bigtriangleup,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(1,1,d,r)$
U	5	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
U	6	$(l,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(l,1,d,r)$
UR	7	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
UR	8	$(l,\bigtriangleup,d,\bigtriangleup)$	$\mathbf{H}_{i}(l,1,d,1)$
R	9	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
R	10	$(l,u,d,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,d,1)$
DR	11	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
DR	12	$(l,u,\bigtriangleup,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,1,1)$
D	13	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
D	14	$(l,u,\bigtriangleup,r)$	$\mathbf{H}_{i}(l,u,1,r)$
DL	15	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	16	$(\bigtriangleup,\bigtriangledown,\bigtriangledown,\bigtriangleup)$	$\mathbf{I}_{i}$
	17	$(\bigtriangleup,u,\bigtriangledown,r)$	$\mathbf{H}_{i}(1,u,1,r)$
L	18	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	19	$(\bigtriangleup,\bigtriangledown,\bigtriangledown,\bigtriangleup)$	$\mathbf{I}_{i}$
	20	$(\bigtriangleup,u,d,r)$	$\mathbf{H}_{i}(1,u,d,r)$
BODY	21	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
BODY	22	$(l,u,d,r)$	$\mathbf{H}_{i}(l,u,d,r)$

Table 6: PEPO rules for the upper-right corner and the right edge. Main differences are in the rules for UL, UR, L, DL, U.

We then describe the FSM rules for the upper-right corner and the right edge of the system as shown in Table 6. Together with the FSM rules for the lower-right corner in Table 7, these FSM rules are to ensure that either only one signal $\bigtriangledown$ is passed down from the left, or exactly one $\hat{H}_{\bm{\alpha}}$ is activated in this part of the system.

Here we explain the rules in Table 6. For the UL site of each box $\mathcal{I}_{\bm{\alpha}}$ , there are four possible incoming signal combinations it can receive. The first possibility is that it receives two $\bigtriangleup$ signals. In this case it means no $\hat{H}_{{\bm{\alpha}}^{\prime}}$ has been activated to the left or from above, we can choose whether to activate $\hat{H}_{\bm{\alpha}}$ or not (rules 1 and 4). If only one of the signals is $\bigtriangledown$ then it means an $\hat{H}_{{\bm{\alpha}}^{\prime}}$ has been activated elsewhere, and it is only allowed to pass the signal $\bigtriangledown$ downward (instead of rightward as in Table 5). This is specified in rule 2. Other rules exist to pass along the signals or are the same as those for the individual $\hat{H}_{\bm{\alpha}}$ . At the DL site, if $\hat{H}_{\bm{\alpha}}$ has been activated, then a signal $\bigtriangledown$ is passed to the next box below (rule 16). The signal $\bigtriangledown$ is only passed along the left boundary, as opposed to along the upper boundary in Table 5.

	rule num	$(l,u,d,r)$	operator
UL of R	1	$(\bigtriangledown,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	2	$(\bigtriangleup,\bigtriangledown,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
	3	$(\bigtriangleup,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(1,1,d,r)$
U of R	4	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
U of R	5	$(l,\bigtriangleup,d,r)$	$\mathbf{H}_{i}(l,1,d,r)$
UR of R	6	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
UR of R	7	$(l,\bigtriangleup,d,\bigtriangleup)$	$\mathbf{H}_{i}(l,1,d,1)$
R of R	8	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
R of R	9	$(l,u,d,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,d,1)$
DR of R	10	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
DR of R	11	$(l,u,\bigtriangleup,\bigtriangleup)$	$\mathbf{H}_{i}(l,u,1,1)$
D of R	12	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
D of R	13	$(l,u,\bigtriangleup,r)$	$\mathbf{H}_{i}(l,u,1,r)$
DL of R	14	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
DL of R	15	$(\bigtriangleup,u,\bigtriangleup,r)$	$\mathbf{H}_{i}(1,u,1,r)$
L of R	16	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
L of R	17	$(\bigtriangleup,u,d,r)$	$\mathbf{H}_{i}(1,u,d,r)$
BODY of R	18	$(\bigtriangleup,\bigtriangleup,\bigtriangleup,\bigtriangleup)$	$\mathbf{I}_{i}$
BODY of R	19	$(l,u,d,r)$	$\mathbf{H}_{i}(l,u,d,r)$

Table 7: PEPO rules for the lower-right corner of the system.

Now all that is left is to formulate the FSM rules for the lower-right corner of the system. The rules are shown in Table 7. If the UL site receives exactly one signal $\bigtriangleup$ and exactly one signal $\bigtriangledown$ form left and above respectively, then this means everything is good and we only need to pass the signal $\bigtriangleup$ to the boundaries without activating the operator (rules 1 and 2). If it receives two signals $\bigtriangleup$ then it means no $\hat{H}_{{\bm{\alpha}}^{\prime}}$ has been activated in any other box, and therefore $\hat{H}_{\bm{\alpha}}$ has to be activated in the current box (rule 3). No signal $\bigtriangledown$ will be passed down beyond this point.

The rules in Tables 5, 6, 7 ensure that exactly one operator $\hat{H}_{\bm{\alpha}}$ is going to be activated in the system. Otherwise some box on the right boundary (upper-right corner, right edge, and lower-right corner) will have two signals $\bigtriangledown$ incoming at the UL site, which is not allowed in the FSM rules and therefore will return a zero operator, making the whole tensor product zero.

Similar to the 1D case, the bond dimension of each PEPO $\hat{H}$ we constructed in this way is the bond dimension of the original PEPO $\hat{H}_{\bm{\alpha}}$ plus 2. In the case of rank-one operators the sum PEPO has a bond dimension of 5.

References

[1] S. Ambikasaran and E. Darve, An $\mathcal{O}(n\log n)$ fast direct solver for partial hierarchically semi-separable matrices, J. Sci. Comput., 57 (2013), pp. 477–501.
[2] L. Balzano and Stephen J Wright, On grouse and incremental svd, in 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), IEEE, 2013, pp. 1–4.
[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.
[4] G. Beylkin and L. Monzón, On approximation of functions by exponential sums, Appl. Comput. Harmon. A., 19 (2005), pp. 17–48.
[5] J. Biamonte and V. Bergholm, Tensor networks in a nutshell, arXiv preprint arXiv:1708.00006, (2017).
[6] D. Braess and W. Hackbusch, Approximation of 1/x by exponential sums in $[1,\infty)$ , IMA J. Numer. Anal., 25 (2005), pp. 685–697.
[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.
[8] Matthew Brand, Incremental singular value decomposition of uncertain data with missing values, in European Conference on Computer Vision, Springer, 2002, pp. 707–720.
[9] J.-F. Cai, E. J. Candès, and Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM J. Optim., 20 (2010), pp. 1956–1982.
[10] E. J. Candès and B. Recht, Exact matrix completion via convex optimization, Found. Comput. Math., 9 (2009), p. 717.
[11] G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White, Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms, J. Chem. Phys., 145 (2016), p. 014102.
[12] G. M. Crosswhite and D. Bacon, Finite automata for caching in matrix product algorithms, Phys. Rev. A, 78 (2008), p. 012356.
[13] 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.
[14] F. Fröwis, V. Nebendahl, and W. Dür, Tensor operators: Constructions and applications for long-range interaction systems, Phys. Rev. A, 81 (2010), p. 062337.
[15] L. Grasedyck and W. Hackbusch, Construction and arithmetics of h-matrices, Computing, 70 (2003), pp. 295–334.
[16] 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.
[17] W. Hackbusch, A sparse matrix arithmetic based on $\mathcal{H}$ -matrices. part i: Introduction to $\mathcal{H}$ -matrices, Computing, 62 (1999), pp. 89–108.
[18] W. Hackbusch, B. N. Khoromskij, and R. Kriemann, Hierarchical matrices based on a weak admissibility criterion, Computing, 73 (2004), pp. 207–243.
[19] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, Unsupervised generative modeling using matrix product states, Phys. Rev. X, 8 (2018), p. 031012.
[20] Y. Hu, Y. Koren, and C. Volinsky, Collaborative filtering for implicit feedback datasets., in ICDM, vol. 8, Citeseer, 2008, pp. 263–272.
[21] C. Lubich, T. Rohwedder, R. Schneider, and B. Vandereycken, Dynamical approximation by hierarchical tucker and tensor-train tensors, SIAM J. Matrix Anal. A., 34 (2013), pp. 470–494.
[22] A. Mnih and R. R. Salakhutdinov, Probabilistic matrix factorization, in Advances in neural information processing systems, 2008, pp. 1257–1264.
[23] M. R. Norman, Colloquium: Herbertsmithite and the search for the quantum spin liquid, Rev. Mod. Phys., 88 (2016), p. 041002.
[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.
[25] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys., 349 (2014), pp. 117–158.
[26] I. V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317.
[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.
[28] S. Östlund and S. Rommer, Thermodynamic limit of density matrix renormalization, Phys. Rev. Lett., 75 (1995), p. 3537.
[29] B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, Matrix product operator representations, New J. Phys., 12 (2010), p. 025012.
[30] M. Rakhuba and I. Oseledets, Calculating vibrational spectra of molecules using tensor train decomposition, J. Chem. Phys., 145 (2016), p. 124101.
[31] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys., 326 (2011), pp. 96–192.
[32] E. Stoudenmire and D. J. Schwab, Supervised learning with tensor networks, in Advances in Neural Information Processing Systems, 2016, pp. 4799–4807.
[33] E. M. Stoudenmire and S. R. White, Sliced basis density matrix renormalization group for electronic structure, Phys. Rev. Lett., 119 (2017), p. 046401.
[34] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions, arXiv preprint cond-mat/0407066, (2004).
[35] 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.
[36] G. Vidal, Classical simulation of infinite-size quantum lattice systems in one spatial dimension, Phys. Rev. Lett., 98 (2007), p. 070201.
[37] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett., 69 (1992), p. 2863.
[38] S. Yan, D. A. Huse, and S. R. White, Spin-liquid ground state of the s= 1/2 kagome heisenberg antiferromagnet, Science, 332 (2011), pp. 1173–1176.

Low-rank representation of tensor network operators with long-range pairwise interactions

Abstract

keywords:

AMS:

1 Introduction

2 Preliminaries

2.1 Matrix product operators

2.2 Projected entangled-pair operators

3 Long-range interaction MPO in 1D systems

Definition 1 (Upper-triangular low-rank matrix).

3.1 An online matrix completion algorithm

3.2 Constructing the MPO

3.3 A modified practical algorithm

4 Representing long-range interaction via the hierarchical low-rank format

4.1 Rank-one MPOs and PEPOs

4.2 Hierarchical low-rank matrix

Definition 2 (Hierarchical low-rank matrix).

4.3 HODLR format

4.4 ℋ\mathcal{H}-matrix format

4.5 Higher dimensional cases

5 Error analysis for the Coulomb matrix

6 Numerical examples

6.1 Quasi 1D system

6.2 2D system

7 Conclusion

Acknowledgments

Appendix A FSM rules for sum of non-overlapping MPOs

Appendix B FSM rules for sum of non-overlapping PEPOs

References

4.4 $\mathcal{H}$ -matrix format