Computing Multi-Eigenpairs of High-Dimensional Eigenvalue Problems Using
  Tensor Neural Networks
In this paper, we propose a type of tensor-neural-network-based machine
learning method to compute multi-eigenpairs of high dimensional eigenvalue
problems without Monte-Carlo procedure. Solving multi-eigenvalues and their
corresponding eigenfunctions is one of the basic tasks in mathematical and
computational physics. With the help of tensor neural network and deep Ritz
method, the high dimensional integrations included in the loss functions of the
machine learning process can be computed with high accuracy. The high accuracy
of high dimensional integrations can improve the accuracy of the machine
learning method for computing multi-eigenpairs of high dimensional eigenvalue
problems. Here, we introduce the tensor neural network and design the machine
learning method for computing multi-eigenpairs of the high dimensional
eigenvalue problems. The proposed numerical method is validated with plenty of
numerical examples.
Yifan Wang, Hehi Xie
# Computing Multi-Eigenpairs of High-Dimensional Eigenvalue Problems Using Tensor Neural Networks+

###### Abstract

In this paper, we propose a type of tensor-neural-network-based machine learning method to compute multi-eigenpairs of high dimensional eigenvalue problems without Monte-Carlo procedure. Solving multi-eigenvalues and their corresponding eigenfunctions is one of the basic tasks in mathematical and computational physics. With the help of tensor neural network and deep Ritz method, the high dimensional integrations included in the loss functions of the machine learning process can be computed with high accuracy. The high accuracy of high dimensional integrations can improve the accuracy of the machine learning method for computing multi-eigenpairs of high dimensional eigenvalue problems. Here, we introduce the tensor neural network and design the machine learning method for computing multi-eigenpairs of the high dimensional eigenvalue problems. The proposed numerical method is validated with plenty of numerical examples.

**Keywords.** high-dimensional eigenvalue problem, tensor neural network, deep Ritz method, high-accuracy, machine learning, multi-eigenpairs.

**AMS subject classifications.** 65N30, 65N25, 65L15, 65B99

## 1 Introduction

In modern sciences and engineers, there exist many high-dimensional problems, which arise from quantum mechanics, statistical mechanics, and financial engineer. There have appeared more and more high dimensional eigenvalue problems along with the developments of science and engineer. Computing the eigenvalues and eigenfunctions of high dimensional operators becomes more and more important in modern scientific computing. The most famous and important high dimensional eigenvalue problem is the Schrodinger equation from quantum mechanics. The classical numerical
methods, such as finite difference, finite element, spectral, can only solve the low dimensional Schrodinger equations for some simple systems. Using these classical numerical methods to solve high dimensional eigenvalue problems will suffer from the so-called curse of dimensionality since the number of degrees of freedom and computational complexity grows exponentially as the dimension increases.

For the high dimensional problems, the Monte-Carlo based methods attract more and more attentions since it provides the possibility to solve high dimensional problems in the stochastic sense. For example, there are two well known variational Monte Carlo (VMC) [6] and diffusion Monte Carlo (DMC) [27] methods for high-dimensional eigenvalue problems in quantum mechanics. Recently, many numerical methods based on machine learning have been proposed to solve the high-dimensional PDEs ([3, 8, 9, 13, 14, 23, 24, 25, 26, 28, 30, 35, 36]). Among these machine learning methods, neural network-based methods attract more and more attention since it can be used to build approximations to the exact solutions of PDEs by machine learning methods. The essential reason is that neural networks can approximate any function given enough parameters. This type of method also provides a possible way to solve many useful high-dimensional PDEs from physics, chemistry, biology, engineering, and so on.

Due to its universal approximation property, the fully-connected neural network (FNN) is the most widely used architecture to build the functions for solving high-dimensional PDEs. There are several types of FNN-based methods such as the well-known deep Ritz [9], deep Galerkin method [30], PINN [28], and weak adversarial networks [35] for solving high-dimensional PDEs by designing different loss functions. Among these methods, the loss functions always include computing high-dimensional integration for the functions defined by FNN. For example, the loss functions of the deep Ritz method require computing the integrations on the high-dimensional domain for the functions which is constructed by FNN. Direct numerical integration for the high-dimensional functions also meets the "curse of dimensionality". Always, the Monte-Carlo method is adopted to do the high-dimensional integration with some types of sampling methods [9, 15]. Due to the low convergence rate of the Monte-Carlo method, the solutions obtained by the FNN-based numerical methods are difficult to obtain high accuracy and stable convergence process. In other words, the Monte-Carlo method decreases computational work in each forward propagation by decreasing the simulation efficiency and stability of the FNN-based numerical methods for solving high-dimensional PDEs.

In [32], based on the deep Ritz method, a type of tensor neural network (TNN) is proposed to build the trial functions for solving high-dimensional PDEs. The TNN is a function being designed by the tensor product operations on the neural networks or by correlated low-rank CAN-DECOMP/PARAFAC (CP) tensor decomposition [17, 22] approximations of FNNs. An important advantage is that we do not need to use Monte-Carlo method to do the integration for the functions which is constructed by TNN. The high dimensional integration of the TNN functions can be decomposed into one-dimensional integrations which can be computed by the highly accurate Gauss type quadrature scheme. The computational work for the integration of the functions by TNN is only a polynomial scale of the dimension, which means the TNN overcomes the "curse of dimensionality" in some sense for solving high-dimensional PDEs. Furthermore, the high accuracy
of high dimensional integration can improve the TNN-based machine learning methods for solving high-dimensional problems. The TNN-based machine learning methods have already been used to solve the smallest eigenvalue and its corresponding eigenfunction of high-dimensional eigenvalue problems [32] and Schrodinger equations in [33].

In this paper, we investigate the applications of TNN for computing multi-eigenpairs of high dimensional eigenvalue problems. The aim of this paper is to propose a universal machine learning solver for computing multi-eigenpairs of high dimensional eigenvalue problems based on the TNN and deep Ritz method without the a priori knowledge of the actual or guessed solutions. Furthermore, the proposed eigensolver can compute multi-eigenpairs of the high dimensional eigenvalue problem with obviously better accuracy than the Montor-Carlo based machine learning methods. In this paper, we will also find that the applications of TNN and the high accuracy of high dimensional integrations bring more choices to define the loss functions for the machine learning methods. For example, the loss function in this paper for computing multi-eigenpairs is more like the classical way for computing eigenvalue problems. The proposed eigensolver provides a new way to solve high dimensional eigenvalue problems from physics, chemistry, biology, engineering, and so on. We will provide plenty of numerical results for the examples from the physics and material sciences.

An outline of the paper goes as follows. In Section 2, we introduce the TNN architecture and its approximation property. In Section 3, the TNN-based machine learning method and corresponding numerical integration schemes are introduced for solving multi-eigenpairs of the high dimensional eigenvalue problems. Section 4 is devoted to providing numerical examples to validate the accuracy and efficiency of the proposed numerical methods. Finally, some concluding remarks are given in the last section.

## 2 Tensor neural network architecture

TNN structure, its approximation property and the computational complexity of related integration have been discussed and investigated in [32]. In this section, we introduce the architecture of TNN. In order to express clearly and facilitate the construction of the TNN method for solving multi-eigenpairs of high dimensional eigenvalue problems, here, we will also elaborate on some important definitions and properties.

TNN is a neural network of low-rank structure, which is built by the tensor product of of several one-dimensional input and multidimensional output subnetworks. Due to the low-rank structure of TNN, an efficient and accurate quadrature scheme can be designed for the TNN-related high dimensional integrations such as the inner product of two TNNs. In [32], we introduce TNN in detail and propose its numerical integration scheme with the polynomial scale computational complexity of the dimension. For each \(i=1,2,\cdots,d\), we use \(\Phi_{i}(x_{i};\theta_{i})=(\phi_{i,1}(x_{i};\theta_{i}),\phi_{i,2}(x_{i}; \theta_{i}),\cdots,\phi_{i,p}(x_{i};\theta_{i}))\) to denote a subnetwork that maps a set \(\Omega_{i}\subset\mathbb{R}\) to \(\mathbb{R}^{p}\), where \(\Omega_{i},i=1,\cdots,d\), can be a bounded interval \((a_{i},b_{i})\), the whole line \((-\infty,+\infty)\) or the half line \((a_{i},+\infty)\). The number of layers and neurons in each layer, the selections of activation functions and other hyperparameters can be different in different subnetworks. In this paper, in order to improve the numerical stability further,
the TNN is defined as follows:

\[\Psi(x;\Theta)=\sum_{j=1}^{p}c_{j}\widehat{\phi}_{1,j}(x_{1};\theta_{1})\widehat{ \phi}_{2,j}(x_{2};\theta_{2})\cdots\widehat{\phi}_{d,j}(x_{d};\theta_{d})=\sum_ {j=1}^{p}c_{j}\prod_{i=1}^{d}\widehat{\phi}_{i,j}(x_{i};\theta_{i}), \tag{2.1}\]

where \(c=\{c_{j}\}_{j=1}^{p}\) is a set of trainable parameters, \(\Theta=\{c,\theta_{1},\cdots,\theta_{d}\}\) denotes all parameters of the whole architecture. For \(i=1,\cdots,d,j=1,\cdots,p\), \(\widehat{\phi}_{i,j}(x_{i},\theta_{i})\) is a normalized functions as follows:

\[\widehat{\phi}_{i,j}(x_{i},\theta_{i})=\frac{\phi_{i,j}(x_{i},\theta_{i})}{ \|\phi_{i,j}(x_{i},\theta_{i})\|_{L^{2}(\Omega_{i})}}. \tag{2.2}\]

In Section 3.2, we will discuss the architectures for each subnetwork \(\Phi_{i}(x_{i};\theta_{i})\) in detail according to different types of \(\Omega_{i}\). The TNN architecture (2.1) and the one defined in [32] are mathematically equivalent, but (2.1) has better numerical stability during the training process. Figure 1 shows the corresponding architecture of TNN. From Figure 1 and numerical tests, we can find the parameters for each rank of TNN are correlated by the FNN, which guarantee the stability of the TNN-based machine learning methods. This is also an important difference from the tensor finite element methods.

In order to show the reasonableness of TNN, we now introduce the approximation property from [32]. Since there exists the isomorphism relation between \(H^{m}(\Omega_{1}\times\cdots\times\Omega_{d})\) and the tensor product space \(H^{m}(\Omega_{1})\otimes\cdots\otimes H^{m}(\Omega_{d})\), the process of approximating the function \(f(x)\in H^{m}(\Omega_{1}\times\cdots\times\Omega_{d})\) by the TNN defined as (2.1) can be regarded as searching for a correlated CP decomposition structure to approximate \(f(x)\) in the space \(H^{m}(\Omega_{1})\otimes\cdots\otimes H^{m}(\Omega_{d})\) with the rank being not greater than \(p\). In [32] we introduce and prove the following approximation result to the functions in the space \(H^{m}(\Omega_{1}\times\cdots\times\Omega_{d})\) under the sense of \(H^{m}\)-norm.

**Theorem 1**.: _Assume that each \(\Omega_{i}\) is an interval in \(\mathbb{R}\) for \(i=1,\cdots,d\), \(\Omega=\Omega_{1}\times\cdots\times\Omega_{d}\), and the function \(f(x)\in H^{m}(\Omega)\). Then for any tolerance \(\varepsilon>0\), there exist a positive integer \(p\) and the corresponding TNN defined by (2.1) such that the following approximation property holds_

\[\|f(x)-\Psi(x;\theta)\|_{H^{m}(\Omega)}<\varepsilon. \tag{2.3}\]

The motivation for employing the TNN architectures is to provide high accuracy and high efficiency in calculating variational forms of high-dimensional problems in which high-dimensional integrations are included. The TNN itself can approximate functions in Sobolev space with respect to \(H^{m}\)-norm. Therefore we naturally put forward an approach to solve high-dimensional PDEs and the ground state eigenpair in [32, 33]. The major contribution of this paper is to propose a TNN-based machine learning method to compute the leading multi-eigenpairs of high-dimensional eigenvalue problems. We will find that the designing of the machine learning method in this paper for computing multi-eigenpairs also shows the advantages of TNN.

## 3 Machine learning method for computing multi-eigenpairs

This section is devoted to introducing the TNN-based machine learning method [11] to compute the multi-eigenpairs of high dimensional eigenvalue problems. We will introduce the way to build
the eigenfunction approximations by TNN functions and the quadrature schemes to compute the inner products of the TNN functions.

### Approximate eigen-subspace by TNNs

In this subsection, we present the TNN-based discretization of eigenvalue problems for solving multi-eigenpairs. Briefly speaking, analogous to the subspace projection method such as the finite element method, instead of using the finite element basis, our approach uses several TNNs as the basis to span a subspace of solution space and restrict the problem into this finite dimensional subspace. Through a machine learning process, an optimal approximation of the eigen-subspace represented by TNNs will be found. Then the final multi-eigenpair approximations are obtained by solving a finite-dimensional matrix eigenvalue problem, which is similar to the Rayleigh-Ritz step in the classical eigensolvers for the matrix eigenvalue problems [29].

Figure 1: Architecture of TNN. Black arrows mean linear transformation (or affine transformation). Each ending node of blue arrows is obtained by taking the scalar multiplication of all starting nodes of blue arrows that end in this ending node. The final output of TNN is derived from the summation of all starting nodes of red arrows.
For generality, we describe the eigenvalue and the TNN-based machine learning method in the abstract way. More specifically, assume Sobolev spaces \(\mathcal{V}\) and \(\mathcal{W}\) are two Hilbert spaces and satisfy \(\mathcal{V}\subset\mathcal{W}\). And let \(a(\cdot,\cdot)\) and \(b(\cdot,\cdot)\) denote two positive definite symmetric bilinear forms on \(\mathcal{V}\times\mathcal{V}\) and \(\mathcal{W}\times\mathcal{W}\), respectively. Furthermore, based on these bilinear forms, we can define the norms on the space \(\mathcal{V}\) and \(\mathcal{W}\) as follows

\[\|v\|_{a} = \sqrt{a(v,v)},\quad\forall v\in\mathcal{V}, \tag{3.1}\] \[\|w\|_{b} = \sqrt{b(w,w)},\quad\forall w\in\mathcal{W}. \tag{3.2}\]

Assume the norm \(\|\cdot\|_{a}\) is relatively compact with respect to the norm \(\|\cdot\|_{b}\)[7].

To describe our method and to build the loss function for computing leading \(k\) eigenpairs briefly, we focus on the following general eigenvalue problem: Find \((\lambda,u)\in\mathbb{R}\times\mathcal{V}\) such that \(b(u,u)=1\) and

\[a(u,v)=\lambda b(u,v),\quad\forall v\in\mathcal{V}. \tag{3.3}\]

It is well known that the eigenvalue problem (3.3) has an eigenvalue sequence \(\{\lambda_{j}\}\) (cf. [2]):

\[0<\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{k}\leq\cdots,\quad\lim_{k \rightarrow\infty}\lambda_{k}=\infty,\]

and associated eigenfunctions

\[u_{1},u_{2},\cdots,u_{k},\cdots,\]

where \(b(u_{i},u_{j})=\delta_{ij}\) (\(\delta_{ij}\) denotes the Kronecker function). In the sequence \(\{\lambda_{j}\}\), the \(\lambda_{j}\) are repeated according to their geometric multiplicity.

The eigenvalues satisfy the minimum-maximum principle [2]

\[\lambda_{k}=\min_{\begin{subarray}{c}\mathcal{V}_{k}\subset \mathcal{V}\\ \dim\mathcal{V}_{k}=k\end{subarray}}\max_{f\in\mathcal{V}_{k}}\mathcal{R}(f)= \max_{f\in\mathrm{span}\{u_{1},\cdots,u_{k}\}}\mathcal{R}(f), \tag{3.4}\]

where \(\mathcal{R}(f)=a(f,f)/b(f,f)\) denotes the Rayleigh quotient of the function \(f\). Let us define \(\mathcal{V}_{k}=\mathrm{span}\{v_{1},\cdots,v_{k}\}\) and \(v_{i}\in\mathcal{V}\) for \(i=1,\cdots,k\) are linearly independent of each other. Furthermore, we also define \(\mathcal{U}_{k}=\mathrm{span}\{u_{1},\cdots,u_{k}\}\) which denotes the eigensubspace with respect to the leading \(k\) eigenvalues. Use the terminology of the subspace approximation to the eigenvalue problems, we define the stiffness matrix \(\mathcal{A}(v_{1},\cdots,v_{k})\) and mass matrix \(\mathcal{B}(v_{1},\cdots,v_{k})\) as follows

\[\mathcal{A}(v_{1},\cdots,v_{k}) =\Big{(}\mathcal{A}_{ij}(v_{1},\cdots,v_{k})\Big{)}_{1\leq i,j \leq k}=\Big{(}a(v_{i},v_{j})\Big{)}_{1\leq i,j\leq k}\in\mathbb{R}^{k\times k}, \tag{3.5}\] \[\mathcal{B}(v_{1},\cdots,v_{k}) =\Big{(}\mathcal{B}_{ij}(v_{1},\cdots,v_{k})\Big{)}_{1\leq i,j \leq k}=\Big{(}b(v_{i},v_{j})\Big{)}_{1\leq i,j\leq k}\in\mathbb{R}^{k\times k}. \tag{3.6}\]

Then due to the minimum-maximum principle, by simple derivation, the summation of the leading \(k\) eigenvalues satisfies the following optimization problem

\[\sum_{i=1}^{k}\lambda_{i}=\min_{\begin{subarray}{c}\mathcal{V}_{k}=\mathrm{ span}\{v_{1},\cdots,v_{k}\}\\ v_{j}\in\mathcal{V},j=1,\cdots,k\end{subarray}}\mathrm{trace}\left(\mathcal{B }^{-1}(v_{1},\cdots,v_{k})\mathcal{A}(v_{1},\cdots,v_{k})\right), \tag{3.7}\]
and the eigensubspace with respect to the leading \(k\) eigenvalues consists with

\[\mathcal{U}_{k}=\text{span}\{u_{1},\cdots,u_{k}\}=\arg\min_{\mathcal{ V}_{k}=\text{span}\{v_{1}\cdots,v_{k}\}\subset\mathcal{V}}\text{trace}\Big{(} \mathcal{B}^{-1}\big{(}v_{1},\cdots,v_{k}\big{)}\mathcal{A}\big{(}v_{1},\cdots, v_{k}\big{)}\Big{)}. \tag{3.8}\]

We approximate \(\mathcal{U}_{k}\) by a space spanned by \(k\) TNNs \(\Psi_{1}(x;\Theta_{1})\), \(\cdots\), \(\Psi_{k}(x;\Theta_{k})\), where each \(\Psi_{\ell}(x;\Theta_{\ell})\) is defined by (2.1) with \(\Phi_{i,\ell}(x_{i};\theta_{i,\ell})=(\phi_{i,1,\ell}(x_{i};\theta_{i,\ell})\), \(\phi_{i,2,\ell}(x_{i};\theta_{i,\ell})\), \(\cdots\), \(\phi_{i,p,\ell}(x_{i};\theta_{i,\ell})\)) as the \(i\)-th subnetwork of the \(\ell\)-th TNN. Then, for \(\ell=1,\cdots,k\), the \(\ell\)-th TNN \(\Psi_{\ell}(x;\Theta_{\ell})\) is denoted as

\[\Psi_{\ell}(x;\Theta_{\ell})=\sum_{j=1}^{p_{\ell}}c_{j,\ell}\prod _{i=1}^{d}\widehat{\phi}_{i,j,\ell}(x_{i};\theta_{i,\ell}), \tag{3.9}\]

where each \(c_{\ell}=\{c_{j,\ell}\}_{j=1}^{p}\) is a set of trainable parameters, each \(\Theta_{\ell}=\{c_{\ell},\theta_{1,\ell},\cdots,\theta_{d,\ell}\}\) denotes all parameters of the \(\ell\)-th TNN. Similar to (2.2), each \(\widehat{\phi}_{i,j,\ell}\) is normalized function. We can select the appropriate activation function such that all \(\Psi_{\ell}(x;\Theta_{\ell}),\ell=1,\cdots,k\) belong to the space \(\mathcal{V}\). Let us define \(p=\max\{p_{1},\cdots,p_{k}\}\). These \(k\) TNNs are trained using the following loss function

\[\text{Loss}\Big{(}\Psi_{1}(x;\Theta_{1}),\cdots,\Psi_{k}(x;\Theta_ {k})\Big{)}\] \[=\text{trace}\Big{(}\mathcal{B}^{-1}\big{(}\Psi_{1}(x;\Theta_{1} ),\cdots,\Psi_{k}(x;\Theta_{k})\big{)}\mathcal{A}\big{(}\Psi_{1}(x;\Theta_{1}),\cdots,\Psi_{k}(x;\Theta_{k})\big{)}\Big{)}. \tag{3.10}\]

Since we do the inner-products of TNN in (3.10), the way to build the loss function is deep Ritz type. In order to assemble the matrices \(\mathcal{A}\) and \(\mathcal{B}\) in (3.10), we need to compute the high dimensional integrations included in the bilinear forms \(a(\cdot,\cdot)\) and \(b(\cdot,\cdot)\) in (3.3). The detailed method for assembling the matrices \(\mathcal{A}\) and \(\mathcal{B}\) in (3.10) will be introduced in Section 3.2. The loss function (3.10) is automatically differentiable thanks to packages that support backpropagation such as TensorFlow and PyTorch. In this paper, the gradient descent (GD) method is adopted to update all trainable parameters

\[\Theta_{\ell}-\eta\nabla_{\Theta_{\ell}}\text{Loss}\to\Theta_{ \ell},\quad\ell=1,\cdots,k, \tag{3.11}\]

where \(\eta\) is the learning rate and adjusted by the ADAM optimizer [20].

After sufficient training steps, we obtain a sequence of parameters \(\{\Theta_{1}^{*},\cdots,\Theta_{k}^{*}\}\) such that the loss function (3.10) arrives its minimum value under the required tolerance. By solving the following finite-dimensional matrix eigenvalue problem

\[\mathcal{A}\big{(}\Psi_{1}(x;\Theta_{1}^{*}),\cdots,\Psi_{k}(x; \Theta_{k}^{*})\big{)}\mathbf{y}=\lambda\mathcal{B}\big{(}\Psi_{1}(x;\Theta_{1 }^{*}),\cdots,\Psi_{k}(x;\Theta_{k}^{*})\big{)}\mathbf{y}, \tag{3.12}\]

we obtain a sequence of eigenvalues

\[0<\widehat{\lambda}_{1}\leq\widehat{\lambda}_{2}\leq\cdots\leq \widehat{\lambda}_{k} \tag{3.13}\]

and corresponding eigenvectors

\[\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{k}, \tag{3.14}\]
where \(\mathbf{y}_{j}=[y_{j,1},\cdots,y_{j,k}]^{\top}\) for \(j=1,\cdots,k\). Then \(\widehat{\lambda}_{1}\), \(\cdots\), \(\widehat{\lambda}_{k}\) can be chosen as the approximations to the first \(k\) eigenvalues \(\lambda_{1}\), \(\cdots\), \(\lambda_{k}\) and the corresponding approximations \(\widehat{u}_{1}\), \(\cdots\), \(\widehat{u}_{k}\) to the first \(k\) eigenfunctions of problem (3.3) can be obtained by the following linear combination procedure

\[\widehat{u}_{j}(x) = \sum_{\ell=1}^{k}y_{j,\ell}\Psi_{\ell}(x;\Theta_{\ell}^{*}). \tag{3.15}\]

**Remark 1**.: _In [36], the loss function is defined by the summation of \(k\) Rayleigh quotients and a penalty term which constrains the \(k\) neural networks to be mutually orthogonal and normalized. Different form [36], in this paper, the loss function (3.10) is defined without the penalty term and the orthogonalization condition of \(k\) TNNs is not imposed directly. The reason is that the use of TNN architectures can provide a high-accuracy and high-efficiency quadrature scheme, as we will see in Section 3.2, for assembling the matrices in the loss function (3.10) and eigenvalue problem (3.12). The penalty term makes optimization process much more difficult and thus affects the final accuracy. The orthogonalization condition is implicitly guaranteed by solving the optimization problem in the machine learning method with the loss function (3.10). The main innovation of the present paper is for the first time to apply the TNN architecture to compute multi-eigenpairs of high dimensional eigenvalue problems with high accuracy. In Section 4.1, we will give an example from [36] to show the advantages of proposed approach._

### Quadrature scheme for inner product

In this subsection, we provide a detailed description to calculate inner products in the loss function (3.10). For simplicity, we are concerned with the following model problem as an example: Find \((\lambda,u)\in\mathbb{R}\times H_{0}^{1}(\Omega)\) such that

\[\begin{cases}-\Delta u+V(x)u=\lambda u,&\text{in}\ \ \Omega,\\ \hskip 14.226378ptu=0,&\text{on}\ \partial\Omega,\end{cases} \tag{3.16}\]

where \(\Omega=\Omega_{1}\times\cdots\times\Omega_{d}\), each \(\Omega_{i},i=1,\cdots,d\), can be a bounded interval \((a_{i},b_{i})\), the whole line \((-\infty,+\infty)\) or the half line \((a_{i},+\infty)\), \(V(x)\in L^{2}(\Omega)\) is a potential function. We assume that the potential \(V(x)\) has separated representation in the tensor product space \(L^{2}(\Omega_{1})\otimes\cdots\otimes L^{2}(\Omega_{d})\) as follows

\[V(x)=\sum_{j=1}^{q}\prod_{i=1}^{d}V_{i,j}(x_{i}), \tag{3.17}\]

where \(V_{i,j}(x_{i})\in L^{2}(\Omega_{i})\). Then the equivalent variational form of the eigenvalue problem (3.16) can be defined as follows: Find \((\lambda,u)\in\mathbb{R}\times H_{0}^{1}(\Omega)\) such that

\[a(u,v)=\lambda b(u,v),\ \ \ \ \forall v\in H_{0}^{1}(\Omega), \tag{3.18}\]

where the bilinear forms \(a(\cdot,\cdot)\) and \(b(\cdot,\cdot)\) are defined as follows

\[a(u,v)=\int_{\Omega}(\nabla u\cdot\nabla v+Vuv)d\Omega,\ \ \ \ \ \ b(u,v)=\int_{\Omega}uvd\Omega. \tag{3.19}\]
Due to the definition of \(k\) TNNs (3.9) and (3.17), entries of the matrix \(\mathcal{A}\) in the loss function (3.10) have following expansions

\[\mathcal{A}_{mn} = a\big{(}\Psi_{m}(x;\Theta_{m}),\Psi_{n}(x;\Theta_{n})\big{)} \tag{3.20}\] \[= \sum_{s=1}^{d}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{ \ell,n}\prod_{i\neq s}^{d}\int_{\Omega_{i}}\phi_{i,j,m}(x_{i};\theta_{i,m}) \phi_{i,\ell,n}(x_{i};\theta_{i,n})dx_{i}\] \[\qquad\cdot\int_{\Omega_{s}}\frac{\partial\phi_{s,j,m}}{\partial x _{s}}(x_{s};\theta_{s,m})\frac{\partial\phi_{s,\ell,n}}{\partial x_{s}}(x_{s };\theta_{s,n})dx_{s}\] \[\quad+\sum_{s=1}^{q}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j, m}c_{\ell,n}\prod_{i=1}^{d}\int_{\Omega_{i}}V_{i,s}(x_{i})\phi_{i,j,m}(x_{i}; \theta_{i,m})\phi_{i,\ell,m}(x_{i};\theta_{i,n})dx_{i},\]

and entries of the matrix \(\mathcal{B}\) can be expanded as follows

\[\mathcal{B}_{mn} = b\big{(}\Psi_{m}(x;\Theta_{m}),\Psi_{n}(x;\Theta_{n})\big{)} \tag{3.21}\] \[= \sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{\ell,n}\prod_{i=1 }^{d}\int_{\Omega_{i}}\phi_{i,j,m}(x_{i};\theta_{i,m})\phi_{i,\ell,n}(x_{i}; \theta_{i,n})dx_{i}.\]

Since only one-dimensional integrations are involved in (3.20) and (3.21), there is no need to use the Monte Carlo procedure to do these high dimensional integrations. In order to guarantee the high accuracy of the high dimensional integrations included in the loss functions, we should use the high order one dimensional quadrature scheme, such as Gauss-type rules [10], for the integrations in (3.20) and (3.21).

#### 3.2.1 Legendre-Gauss quadrature scheme for bounded domain

Without loss of generality, we decompose each \(\Omega_{i}\) into \(M_{i}\) equal subintervals with length \(h_{i}=|\Omega_{i}|/M_{i}\) and choose \(N_{i}\) Legendre-Gauss points in each subinterval. Denote the total quadrature points and corresponding weights as follows

\[\Big{\{}x_{i}^{(n_{i})}\Big{\}}_{n_{i}=1}^{M_{i}N_{i}},\quad\Big{\{}w_{i}^{(n _{i})}\Big{\}}_{n_{i}=1}^{M_{i}N_{i}},\quad i=1,2,\cdots,d, \tag{3.22}\]

and define \(N=\max\{N_{1},\cdots,N_{d}\}\), \(\underline{N}=\min\{N_{1},\cdots,N_{d}\}\), \(M=\max\{M_{1},\cdots,M_{d}\}\) and \(\underline{M}=\min\{M_{1},\cdots,M_{d}\}\). Then using the quadrature scheme (3.22), entries of the matrix \(\mathcal{A}\) defined by (3.20) have the following numerical format

\[\mathcal{A}_{mn}=a\big{(}\Psi_{m}(x;\Theta_{m}),\Psi_{n}(x; \Theta_{n})\big{)}\] \[\approx\sum_{s=1}^{d}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j, m}c_{\ell,n}\prod_{i\neq s}^{d}\sum_{n_{i}=1}^{M_{i}N_{i}}w_{i}^{(n_{i})} \phi_{i,j,m}(x_{i}^{(n_{i})};\theta_{i,m})\phi_{i,\ell,n}(x_{i}^{(n_{i})}; \theta_{i,n})\] \[\cdot\sum_{n_{s}=1}^{M_{s}N_{s}}w_{s}^{(n_{s})}\frac{\partial \phi_{s,j,m}}{\partial x_{s}}(x_{s}^{(n_{s})};\theta_{s,m})\frac{\partial\phi _{s,\ell,n}}{\partial x_{s}}(x_{s}^{(n_{s})};\theta_{s,n})\] \[+\sum_{s=1}^{q}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{ \ell,n}\prod_{i=1}^{d}\sum_{n_{i}=1}^{M_{i}N_{i}}w_{i}^{(n_{i})}V_{i,s}(x_{i}^ {(n_{i})})\phi_{i,j,m}(x_{i}^{(n_{i})};\theta_{i,m})\phi_{i,\ell,m}(x_{i}^{(n_ {i})};\theta_{i,n}), \tag{3.23}\]
and entries of the matrix \(\mathcal{B}\) defined by (3.21) can be computed as follows

\[\mathcal{B}_{mn} = b\big{(}\Psi_{m}(x;\Theta_{m}),\Psi_{n}(x;\Theta_{n})\big{)} \tag{3.24}\] \[\approx \sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{\ell,n}\prod_{i=1 }^{d}\sum_{n_{i}=1}^{M_{i}N_{i}}w_{i}^{(n_{i})}\phi_{i,j,m}(x_{i}^{(n_{i})}; \theta_{i,m})\phi_{i,\ell,n}(x_{i}^{(n_{i})};\theta_{i,n}).\]

Then the complete TNN-based algorithm for computing the first \(k\) eigenpairs can be summarized in Algorithm 1.

```
1. Initialization step: Build \(k\) initial TNNs \(\Psi_{1}^{(0)}(x;\Theta_{1}^{(0)}),\cdots,\Psi_{k}^{(0)}(x;\Theta_{k}^{(0)})\) as in (3.9), maximum training steps \(L\), learning rate \(\eta\), quadrature points and weights (3.22).
2. Assemble matrices \(\mathcal{A}\big{(}\Psi_{1}^{(\ell)}(x;\Theta_{1}^{(\ell)}),\cdots,\Psi_{k}^{( \ell)}(x;\Theta_{k}^{(\ell)})\big{)}\) and \(\mathcal{B}\big{(}\Psi_{1}^{(\ell)}(x;\Theta_{1}^{(\ell)}),\cdots,\Psi_{k}^{( \ell)}(x;\Theta_{k}^{(\ell)})\big{)}\) according to quadrature schemes (3.23) and (3.24), respectively.
3. Compute the value of the loss function (3.10), where the matrix \(\mathcal{B}^{-1}\mathcal{A}=:\mathcal{C}\) is obtained by solving the matrix equation \(\mathcal{BC}=\mathcal{A}\).
4. Compute the gradient of the loss function with respect to parameters \(\{\Theta_{1},\cdots,\Theta_{k}\}\) by automatic differentiation and update parameters \(\{\Theta_{1},\cdots,\Theta_{k}\}\) by (3.11).
5. Set \(\ell=\ell+1\) and go to Step 2 for the next step until \(\ell=L\).
6. Post-processing step: Solve the matrix eigenvalue problrm (3.12) and compute the first \(k\) eigenpair approximations by (3.15).
```

**Algorithm 1**TNN-based method for the first \(k\) eigenpairs

It is worth mentioning that the quadrature schemes (3.23) and (3.24) which use 1-dimensional \(M_{i}N_{i}\) quadrature points in each dimension are actually equivalent to implementing the following \(d\)-dimensional full tensor quadrature scheme

\[\mathcal{A}_{mn} \approx \sum_{n\in\mathcal{N}}w^{(n)}\nabla\Psi_{m}(x^{(n)};\Theta_{m}) \cdot\nabla\Psi_{n}(x^{(n)};\Theta_{n})\] \[+\sum_{n\in\mathcal{N}}w^{(n)}V(x^{(n)})\Psi_{m}(x^{(n)};\Theta_ {m})\Psi_{n}(x^{(n)};\Theta_{n}),\] \[\mathcal{B}_{mn} \approx \sum_{n\in\mathcal{N}}w^{(n)}\Psi_{m}(x^{(n)};\Theta_{m})\Psi_{n }(x^{(n)};\Theta_{n}), \tag{3.25}\]

where \(d\)-dimensional quadrature points and weights are defined as follows

\[\left\{x^{(n)}\right\}_{n\in\mathcal{N}} = \left\{\left\{x_{1}^{(n_{1})}\right\}_{n_{1}=1}^{M_{1}N_{1}} \times\left\{x_{2}^{(n_{2})}\right\}_{n_{2}=1}^{M_{2}N_{2}}\times\cdots\times \left\{x_{d}^{(n_{d})}\right\}_{n_{d}=1}^{M_{d}N_{d}}\right\},\] \[\left\{w^{(n)}\right\}_{n\in\mathcal{N}} = \left\{\left\{w_{1}^{(n)}\right\}_{n_{1}=1}^{M_{1}N_{1}}\times \left\{w_{2}^{(n_{2})}\right\}_{n_{2}=1}^{M_{2}N_{2}}\times\cdots,\times \left\{w_{d}^{(n_{d})}\right\}_{n_{d}=1}^{M_{d}N_{d}}\right\}.\]
The quadrature points (3.26) for general \(d\)-dimensional integrand has accuracy \(\mathcal{O}(h^{2\underline{N}}/(2\underline{N})!)\) but \(\mathcal{N}\) has \(\mathcal{O}(M^{d}N^{d})\) elements. This becomes too large for \(d\gg 1\) and makes the approximation of integration intractable. Thanks to TNN having the low-rank tensor type structure, using the full tensor quadrature scheme for the inner product of two TNNs has the splitting schemes (3.23) and (3.24). In quadrature schemes (3.23) and (3.24), assembling matrices \(\mathcal{A}\) and \(\mathcal{B}\) costs only \(\mathcal{O}(d^{2}p^{2}MN+dp^{2}qMN)\) and \(\mathcal{O}(dp^{2}MN)\) operations, respectively. Furthermore, due to the equivalence of using full tensor quadrature points (3.26), quadrature schemes (3.23) and (3.24) also have accuracy \(\mathcal{O}(h^{2\underline{N}}/(2\underline{N})!)\). This means the TNN-based method proposed in this paper overcomes the "curse of dimensionality" in the sense of numerical integration.

For generality, in the next two sections, we will also introduce the integration of TNN on unbounded domains. The corresponding complexity analysis will be omitted since they are similar to that in this section.

#### 3.2.2 Hermite-Gauss quadrature scheme for the whole line

For case \(\Omega_{i}=(-\infty,+\infty)\), we use Hermite-Gauss quadrature to assemble matrix \(\mathcal{A}\) and \(\mathcal{B}\). Hermite-Gauss quadrature scheme satisfies following property.

**Lemma 1**.: _[_31_, Theorem 7.3]_ _Let \(\{x^{(k)}\}_{k=0}^{N}\) be the zeros of the \((N+1)\)-th order Hermite polynomial \(H_{N+1}(x)\), and let \(\{w^{(k)}\}_{k=1}^{N}\) be given by_

\[w^{(k)}=\frac{\sqrt{\pi}2^{N}N!}{(N+1)H_{N}^{2}(x^{(k)})},\quad 0 \leq k\leq N. \tag{3.27}\]

_Then we have_

\[\int_{-\infty}^{+\infty}p(x)e^{-x^{2}}dx=\sum_{k=0}^{N}p(x^{(k)})w ^{(k)},\quad\forall p\in P_{2N+1}, \tag{3.28}\]

_where \(P_{2N+1}\) denotes the set of polynomial functions of degree less than \(2N+1\)._

For using Hermite-Gauss quadrature, the \(i\)-th subnetwork of \(\ell\)-th TNN is defined as follows

\[\Phi_{i,\ell}(x_{i};\theta_{i,\ell}) = \big{(}\phi_{i,1,\ell}(x_{i};\theta_{i,\ell}),\phi_{i,2,\ell}(x_ {i};\theta_{i,\ell}),\cdots,\phi_{i,p,\ell}(x_{i};\theta_{i,\ell})\big{)}\] \[= e^{-\frac{\beta_{i}^{2}x_{i}^{2}}{2}}\big{(}\varphi_{i,1,\ell}( \beta_{i}x_{i};\theta_{i,\ell}),\varphi_{i,2,\ell}(\beta_{i}x_{i};\theta_{i, \ell}),\cdots,\varphi_{i,p,\ell}(\beta_{i}x_{i};\theta_{i,\ell})\big{)},\]

where \(\beta_{i}\) is trainable parameter and \(\varphi_{i,\ell}=\big{(}\varphi_{i,1,\ell}(\beta_{i}x_{i};\theta_{i,\ell}), \varphi_{i,2,\ell}(\beta_{i}x_{i};\theta_{i,\ell}),\cdots,\varphi_{i,p,\ell}( \beta_{i}x_{i};\theta_{i,\ell})\big{)}\) is a fully-connected neural network which maps \(\mathbb{R}\) to \(\mathbb{R}^{p}\). Then the \(\ell\)-th TNN can be written as

\[\Psi_{\ell}(x;\Theta_{\ell})=\sum_{j=1}^{p}c_{j,\ell}\prod_{i=1}^ {d}e^{-\frac{\beta_{i}^{2}x_{i}^{2}}{2}}\widehat{\varphi}_{i,j,\ell}(\beta_{i }x_{i};\theta_{i,\ell}), \tag{3.29}\]

where \(\widehat{\varphi}_{i,j,\ell}\) satisfies normalization property \(\|e^{-\frac{\beta_{i}^{2}x_{i}^{2}}{2}}\widehat{\varphi}_{i,j,\ell}(\beta_{i} x_{i};\theta_{i,\ell})\|_{L^{2}(\Omega_{i})}=1\) and is defined as follows

\[\widehat{\varphi}_{i,j,\ell}(\beta_{i}x_{i};\theta_{i,\ell})= \frac{\varphi_{i,j,\ell}(\beta_{i}x_{i};\theta_{i,\ell})}{\Big{\|}e^{-\frac{ \beta_{i}^{2}x_{i}^{2}}{2}}\varphi_{i,j,\ell}(\beta_{i}x_{i};\theta_{i,\ell}) \Big{\|}_{L^{2}(\Omega_{i})}}. \tag{3.30}\]
For \(i=1,\cdots,d\), let \(\{z_{i}^{(k_{i})}\}_{k_{i}=1}^{N_{i}}\) and \(\{w_{i}^{(k_{i})}\}_{k_{i}=1}^{N_{i}}\) be the Hermite-Gauss quadrature points and weights, respectively. Under coordinate transformation \(z_{i}=\beta_{i}x_{i}\), using the quadrature scheme (3.28), entries of the matrix \(\mathcal{A}\) defined by (3.20) have the following numerical format

\[\mathcal{A}_{mn}=\sum_{i=1}^{d}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p _{n}}c_{j,m}c_{\ell,n}\prod_{i\neq s}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_ {i}}\widehat{\varphi}_{i,j,m}(z_{i}^{(k_{i})};\theta_{i,m})\widehat{\varphi}_{ i,\ell,n}(z_{i}^{(k_{i})};\theta_{i,n})w_{i}^{(k_{i})}\] \[\cdot\frac{1}{\beta_{s}}\sum_{k_{s}}^{N_{s}}\big{(}\widehat{ \varphi}^{\prime}_{s,j,m}(z_{s}^{(k_{s})};\theta_{s,m})-\beta_{s}\widehat{ \varphi}_{s,j,m}(z_{s}^{(k_{s})};\theta_{s,m})\big{)}\] \[\cdot\big{(}\widehat{\varphi}^{\prime}_{s,j,m}(z_{s}^{(k_{s})}; \theta_{s,m})-\beta_{s}\widehat{\varphi}_{s,j,m}(z_{s}^{(k_{s})};\theta_{s,m}) \big{)}w_{s}^{(k_{s})}\] \[+\sum_{s=1}^{q}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{ \ell,n}\prod_{i=1}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_{i}}V_{i,s}(\frac{ z_{i}^{(k_{i})}}{\beta_{i}})\widehat{\varphi}_{i,j,m}(z_{i}^{(k_{i})};\theta_{i,m}) \widehat{\varphi}_{i,\ell,n}(z_{i}^{(k_{i})};\theta_{i,n})w_{i}^{(k_{i})}. \tag{3.31}\]

And entries of the matrix \(\mathcal{B}\) defined by (3.21) have the following numerical format

\[\mathcal{B}_{mn}=\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{\ell,n}\prod _{i=1}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_{i}}\widehat{\varphi}_{i,j,m}(z _{i}^{(k_{i})};\theta_{i,m})\widehat{\varphi}_{i,\ell,n}(z_{i}^{(k_{i})}; \theta_{i,n})w_{i}^{(k_{i})}. \tag{3.32}\]

#### 3.2.3 Laguerre-Gauss quadrature scheme for the half line

For case \(\Omega_{i}=(0,+\infty)\), we use Laguerre-Gauss quadrature to assemble matrices \(\mathcal{A}\) and \(\mathcal{B}\). Laguerre-Gauss quadrature satisfies the following theorem

**Lemma 2**.: _[_31_, Theorem 7.1]_ _Let \(\{x^{(k)}\}_{k=0}^{N}\) be the zeros of the \((N+1)\)-th order Laguerre polynomial \(\mathcal{L}_{N+1}(x)\), and let \(\{w^{(k)}\}_{k=0}^{N}\) be given by_

\[w^{(k)}=\frac{\Gamma(N+1)}{(N+1)(N+1)!}\frac{x^{(k)}}{\big{[} \mathcal{L}_{N}(x^{(k)})\big{]}^{2}},\quad 0\leq 0\leq N. \tag{3.33}\]

_Then we have_

\[\int_{-\infty}^{+\infty}p(x)e^{-x}dx=\sum_{k=0}^{N}p(x^{(k)})w^{(k)},\quad \forall p\in P_{2N+1}, \tag{3.34}\]

_where \(P_{2N+1}\) denotes the set of polynomial functions of degree less than \(2N+1\)._

For using Laguerre-Gauss quadrature, the \(i\)-th subnetwork of the \(\ell\)-th TNN is defined as follows

\[\Phi_{i,\ell}(x_{i};\theta_{i,\ell}) = \big{(}\phi_{i,1,\ell}(x_{i};\theta_{i,\ell}),\phi_{i,2,\ell}(x_{i };\theta_{i,\ell}),\cdots,\phi_{i,p,\ell}(x_{i};\theta_{i,\ell})\big{)}\] \[= e^{-\frac{\beta_{i}x_{i}}{2}}\big{(}\varphi_{i,1,\ell}(\beta_{i} x_{i};\theta_{i,\ell}),\varphi_{i,2,\ell}(\beta_{i}x_{i};\theta_{i,\ell}),\cdots, \varphi_{i,p,\ell}(\beta_{i}x_{i};\theta_{i,\ell})\big{)},\]

where \(\beta_{i}\) is trainable parameter and \(\varphi_{i,\ell}=\big{(}\varphi_{i,1,\ell}(\beta_{i}x_{i};\theta_{i,\ell}), \varphi_{i,2,\ell}(\beta_{i}x_{i};\theta_{i,\ell}),\cdots,\varphi_{i,p,\ell}( \beta_{i}x_{i};\theta_{i,\ell})\big{)}\) is a fully-connected neural network which maps \(\mathbb{R}\) to \(\mathbb{R}^{p}\). Then the \(\ell\)-th TNN can be written as

\[\Psi_{\ell}(x;\Theta_{\ell})=\sum_{j=1}^{p}c_{j,\ell}\prod_{i=1}^{ d}e^{-\frac{\beta_{i}x_{i}}{2}}\widehat{\varphi}_{i,j,\ell}(\beta_{i}x_{i}; \theta_{i,\ell}), \tag{3.35}\]
where \(\widehat{\varphi}_{i,j,\ell}\) satisfies normalization property \(\|e^{-\frac{\beta_{i}x_{i}}{2}}\widehat{\varphi}_{i,j,\ell}(\beta_{i}x_{i};\theta _{i,\ell})\|_{L^{2}(\Omega_{i})}=1\) and is defined as follows

\[\widehat{\varphi}_{i,j,\ell}(\beta_{i}x_{i};\theta_{i,\ell})=\frac{\varphi_{i, j,\ell}(\beta_{i}x_{i};\theta_{i,\ell})}{\left\|e^{-\frac{\beta_{i}x_{i}}{2}} \varphi_{i,j,\ell}(\beta_{i}x_{i};\theta_{i,\ell})\right\|_{L^{2}(\Omega_{i})}}. \tag{3.36}\]

For \(i=1,\cdots,d\), let \(\{z_{i}^{(k_{i})}\}_{k_{i}=1}^{N_{i}}\) and \(\{w_{i}^{(k_{i})}\}_{k_{i}=1}^{N_{i}}\) be Laguerre-Gauss quadrature points and weights, respectively. Under coordinate transformation \(z_{i}=\beta_{i}x_{i}\), using the quadrature scheme (3.28), entries of the matrix \(\mathcal{A}\) defined by (3.20) have the following numerical format

\[\mathcal{A}_{mn}=\sum_{i=1}^{d}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p _{n}}c_{j,m}c_{\ell,n}\prod_{i\neq s}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_ {i}}\widehat{\varphi}_{i,j,m}(z_{i}^{(k_{i})};\theta_{i,m})\widehat{\varphi}_ {i,\ell,n}(z_{i}^{(k_{i})};\theta_{i,n})w_{i}^{(k_{i})}\] \[\cdot\frac{1}{\beta_{s}}\sum_{k_{s}}^{N_{s}}\big{(}\widehat{ \varphi}_{s,j,m}^{\prime}(z_{s}^{(k_{s})};\theta_{s,m})-\frac{1}{2}\widehat{ \varphi}_{s,j,m}(z_{s}^{(k_{s})};\theta_{s,m})\big{)}\] \[\cdot\big{(}\widehat{\varphi}_{s,j,m}^{\prime}(z_{s}^{(k_{s})}; \theta_{s,m})-\frac{1}{2}\widehat{\varphi}_{s,j,m}(z_{s}^{(k_{s})};\theta_{s,m })\big{)}w_{s}^{(k_{s})}\] \[+\sum_{s=1}^{q}\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{ \ell,n}\prod_{i=1}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_{i}}V_{i,s}(\frac{z _{i}^{(k_{i})}}{\beta_{i}})\widehat{\varphi}_{i,j,m}(z_{i}^{(k_{i})};\theta_{i,m})\widehat{\varphi}_{i,\ell,n}(z_{i}^{(k_{i})};\theta_{i,n})w_{i}^{(k_{i})}. \tag{3.37}\]

And entries of the matrix \(\mathcal{B}\) defined by (3.21) have the following numerical format

\[\mathcal{B}_{mn}=\sum_{j=1}^{p_{m}}\sum_{\ell=1}^{p_{n}}c_{j,m}c_{\ell,n}\prod _{i=1}^{d}\frac{1}{\beta_{i}}\sum_{k_{i}=1}^{N_{i}}\widehat{\varphi}_{i,j,m}( z_{i}^{(k_{i})};\theta_{i,m})\widehat{\varphi}_{i,\ell,n}(z_{i}^{(k_{i})}; \theta_{i,n})w_{i}^{(k_{i})}. \tag{3.38}\]

## 4 Numerical examples

In this section, we provide several examples to investigate the performance of the TNN-based eigensolver proposed in this paper. To show the convergence behavior and accuracy of our method, we define the relative errors for the approximated eigenvalues \(\widehat{\lambda}_{\ell}\) and eigenfunctions \(\mathbf{y}_{\ell}\) as follows

\[\mathrm{err}_{\lambda,\ell}:=\frac{|\widehat{\lambda}_{\ell}-\lambda_{\ell}|} {|\lambda_{\ell}|},\ \mathrm{err}_{L^{2},\ell}:=\frac{\|u_{\ell}-\mathcal{Q}_{\ell}u_{\ell}\|_{L^{2}( \Omega)}}{\|u_{\ell}\|_{L^{2}(\Omega)}},\ \mathrm{err}_{H^{1},\ell}:=\frac{|u_{\ell}- \mathcal{P}_{\ell}u_{\ell}|_{H^{1}(\Omega)}}{|u_{\ell}|_{H^{1}(\Omega)}},\ \ell=1, \cdots,k, \tag{4.1}\]

where \(\lambda_{\ell}\) and \(u_{\ell}\) are reference eigenvalues and eigenfunctions. In the first two examples, the reference eigenpairs are obtained by high order finite element methods on the meshes with a sufficiently small mesh size. As for the harmonic oscillator problems, the exact eigenpairs are chosen as references. In (4.1), \(\mathcal{P}_{\ell}:H_{0}^{1}(\Omega)\to M_{\ell}\) and \(\mathcal{Q}:H_{0}^{1}(\Omega)\to M_{\ell}\) denote \(L^{2}(\Omega)\) projection operator and \(H^{1}(\Omega)\) projection operator, respectively, to the approximate eigenspace corresponding to the eigenvalue \(\lambda_{\ell}\). These two projection operators are defined as follows

\[\left\langle\mathcal{P}_{\ell}u,v\right\rangle_{L^{2}} = \left\langle u,v\right\rangle_{L^{2}}:=\int_{\Omega}uvdx,\quad \forall v\in M_{\ell}\ \ \mathrm{for}\ u\in H_{0}^{1}(\Omega), \tag{4.2}\] \[\left\langle\mathcal{Q}_{\ell}u,v\right\rangle_{H^{1}} = \left\langle u,v\right\rangle_{H^{1}}:=\int_{\Omega}\nabla u\cdot \nabla vdx,\quad\forall v\in M_{\ell}\ \ \mathrm{for}\ u\in H_{0}^{1}(\Omega). \tag{4.3}\]
In implementation, we use the quadrature scheme similar to that in (3.23) and (3.24), (3.31) and (3.32), (3.37) and (3.38), to compute \(\operatorname{err}_{L^{2},\ell}\) and \(\operatorname{err}_{H^{1},\ell}\) with the same tensor product quadrature points and weights as computing the loss functions if the reference solution \(u_{\ell}(x)\) has a low-rank representation, otherwise we only report \(\operatorname{err}_{\lambda,\ell}\). With the help of Theorem 2 in [32], the high efficiency and accuracy for computing \(\operatorname{err}_{L^{2},\ell}\) and \(\operatorname{err}_{H^{1},\ell}\) can be guaranteed.

### Infinitesimal generators of metastable diffusion processes

In this subsection, we study the eigenvalue problem from [36] associated with the operator

\[\mathcal{L}_{d}=\nabla V_{d}\cdot\nabla-\Delta. \tag{4.4}\]

The potential functions \(V_{d}:\mathbb{R}^{d}\to\mathbb{R}\) for \(d=2,50,100\) are defined as follows

\[V_{d}(x)=V(\theta)+2(r-1)^{2}+5e^{-5r^{2}}+5\sum_{i=3}^{d}x_{i}^{2},\quad \forall x=(x_{1},x_{2},\cdots,x_{d})\in\mathbb{R}^{d}, \tag{4.5}\]

where \((r,\theta)\in[0,+\infty)\times[-\pi,\pi)\) denotes the polar coordinates which relate to the first two dimensional Euclid space \((x_{1},x_{2})\in\mathbb{R}^{2}\) by

\[x_{1}=r\cos\theta,\quad x_{2}=r\sin\theta, \tag{4.6}\]

and \(V:[-\pi,\pi)\to\mathbb{R}\) is a double-well potential function which is defined as follows

\[V(\theta)=\begin{cases}\left[1-\left(\frac{3\theta}{\pi}+1\right)^{2}\right]^ {2},&\theta\in\left[-\pi,-\frac{\pi}{3}\right),\\ \frac{1}{5}\left(3-2\cos(3\theta)\right),&\theta\in\left[-\frac{\pi}{3},\frac{ \pi}{3}\right),\\ \left[1-\left(\frac{3\theta}{\pi}-1\right)^{2}\right]^{2},&\theta\in\left[ \frac{\pi}{3},\pi\right).\end{cases} \tag{4.7}\]

The reference eigenvalue for \(d=2\) is obtained by using the third order conforming finite element method with the mesh size \(h=\frac{3}{1024}\sqrt{2}\) on the domain \([-3,3]^{2}\subset\mathbb{R}^{2}\). Here we use the open parallel finite element package OpenPFEM [34] to do the discretization and then using Krylovschur method from SLEPc [16] to solve the corresponding algebraic eigenvalue problem 1. In this way, we obtain the first three eigenvalues

Footnote 1: We express our thanks to Yangfei Liao for this computation

\[\lambda_{1}=0.21881493133369,\quad\lambda_{2}=0.76371970025476,\quad\lambda_{ 3}=2.79019347384363,\]

as the reference values for our numerical investigation. The corresponding eigenfunctions are shown in the first column of Figure 2.

In implementation, we use 3 TNNs to learn the lowest 3 eigenvalues. For \(d=2,50,100\) cases, each TNN has depth 3 and width 20 and the rank is chosen to be \(p=10\). The Adam optimizer is
employed with a learning rate \(0.003\). We use the Adam optimizer in the first \(100000\) steps and then the LBFGS in the subsequent \(10000\) steps. The final eigenvalue approximations are represented in Table 1 and the corresponding eigenfunction approximations are shown in Figure 2, where we can find the proposed method has obviously better accuracy than that in [36].
### Harmonic oscillator problems

In this subsection, we consider the Schrodinger equation associated with the \(d\)-dimensional harmonic oscillator and the corresponding Hamiltonian operator reads as

\[H=-\frac{1}{2}\sum_{i=1}^{d}\nabla_{i}^{2}+\frac{1}{2}x^{T}Ax, \tag{4.8}\]

where \(x=(x_{1},x_{2},\cdots,x_{d})^{T}\), \(A=(a_{ij})_{d\times d}\in\mathbb{R}^{d\times d}\) is a symmetric positive definite matrix.

The exact wavefunctions and the corresponding energy (i.e. eigenvalue) of different states can be obtained with the similar way in [4]. According to the property of the symmetric positive definite matrix, there exists an orthogonal matrix \(Q=(q_{ij})_{d\times d}\in\mathbb{R}^{d\times d}\) such that

\[Q^{T}AQ=\text{diag}\big{\{}\mu_{1},\mu_{2},\cdots,\mu_{d}\big{\}}, \tag{4.9}\]

where \(\mu_{j}\) and \(\mathbf{q}_{j}=(q_{1j},q_{2j},\cdots,q_{d})^{T}\), \(j=1,2,\cdots,d\) are eigenvalues and corresponding normalized eigenvectors of the matrix \(A\). Then under the rotation transformation \(y=Q^{T}x\), the Hamiltonian operator (4.8) in the coordinate system \((y_{1},y_{2},\cdots,y_{d})\) can be written as the following decoupled harmonic form

\[H=-\frac{1}{2}\sum_{i=1}^{d}\nabla_{i}^{2}+\frac{1}{2}\sum_{i=1 }^{d}\mu_{i}y_{i}^{2}. \tag{4.10}\]

Then the exact wavefunction of state \((n_{1},n_{2},\cdots,n_{d})\) can be immediately obtained as follows

\[\Psi_{n_{1},n_{2},\cdots,n_{d}}(y_{1},y_{2},\cdots,y_{d})=\prod_{i =1}^{d}\mathcal{H}_{n_{i}}(\mu_{i}^{1/4}y_{i})e^{-\mu_{i}^{1/2}y_{i}^{2}/2}, \tag{4.11}\]

and the corresponding exact energy is

\[E_{n_{1},n_{2},\cdots,n_{d}}=\sum_{i=1}^{d}\Big{(}\frac{1}{2}+n_ {i}\Big{)}\mu_{i}^{1/2}, \tag{4.12}\]

where \(y_{i}=\sum_{j=1}^{d}q_{ij}x_{j}\) and \(\mathcal{H}_{n}\) is physicists' Hermite polynomials [1].

In all three examples of this subsection, we adopt the same parameters \(a_{ij}\) as in [25] and calculate the lowest 16 energy states as [25] done. The aim here is to demonstrate the performance of the proposed approach through these comparisons.

#### 4.2.1 Two-dimensional harmonic oscillator

We first examine a simple case of the two-dimensional harmonic oscillator problem

\[-\frac{1}{2}\Big{(}\nabla_{1}^{2}+\nabla_{2}^{2}\Big{)}\Psi+ \frac{1}{2}(x_{1}^{2}+x_{2}^{2})\Psi=E\Psi. \tag{4.13}\]

Since the problem is essentially decoupled, the exact eigenfunction has low-rank representation in coordinate \((x_{1},x_{2})\) as follows

\[\Psi_{n_{1},n_{2}}=\mathcal{H}_{n_{1}}(x_{1})e^{-x_{1}^{2}/2} \mathcal{H}_{n_{2}}(x_{2})e^{-x_{2}^{2}/2}. \tag{4.14}\]
The corresponding exact energy is

\[E_{n_{1},n_{2}}=\Big{(}\frac{1}{2}+n_{1}\Big{)}+\Big{(}\frac{1}{2}+n_{2}\Big{)}. \tag{4.15}\]

We use 16 TNNs to learn the lowest 16 energy states and each subnetwork of a single TNN with depth 3 and width 50, and \(p=20\). The Adam optimizer is employed with a learning rate 0.001 in the first 500000 epochs and then the L-BFGS in the subsequent 10000 steps to produce the final result. The Hermite-Gauss quadrature scheme with 99 points are adopted in each dimension.

The corresponding numerical results are shown in Table 2, where we can find the proposed TNN-based machine learning method has obvious better accuracy than that in [25], where the accuracy is 1.0e-2 by the Monte-Carlo based machine learning methods.

#### 4.2.2 Two-dimensional coupled harmonic oscillator

Since TNN-based methods carry out operations separately in each dimension, it is not unexpected that our method has an impressive performance in the completely decoupled case. To show the generalities of TNN-based machine learning method, the next two examples are to compute the multi-states of operators with coupled oscillators.

First, we consider the following two-dimensional eigenvalue problem of the operator with coupled

\begin{table}
\begin{tabular}{c c c c c c c} \hline \(n\) & \((n_{1},n_{2})\) & Exact \(E_{n}\) & Approx \(E_{n}\) & err\({}_{E}\) & err\({}_{L^{2}}\) & err\({}_{H^{1}}\) \\ \hline
0 & (0,0) & 1.0 & 1.00000000000441 & 4.414e-13 & 2.935e-07 & 7.314e-07 \\
1 & (0,1) & 2.0 & 2.0000000013887 & 6.944e-11 & 5.889e-06 & 1.021e-05 \\
2 & (1,0) & 2.0 & 2.00000000369235 & 1.846e-10 & 9.371e-06 & 1.550e-05 \\
3 & (0,2) & 3.0 & 3.000000000851601 & 2.839e-10 & 1.524e-05 & 2.180e-05 \\
4 & (1,1) & 3.0 & 3.000000001492529 & 4.975e-10 & 1.970e-05 & 2.958e-05 \\
5 & (2,0) & 3.0 & 3.000000005731970 & 1.911e-09 & 3.964e-05 & 5.852e-05 \\
6 & (0,3) & 4.0 & 4.00000000287978 & 7.200e-11 & 8.690e-06 & 1.312e-05 \\
7 & (1,2) & 4.0 & 4.000000000727938 & 1.820e-10 & 1.374e-05 & 1.838e-05 \\
8 & (2,1) & 4.0 & 4.000000001748148 & 4.370e-10 & 2.272e-05 & 3.100e-05 \\
9 & (3,0) & 4.0 & 4.000000005090556 & 1.273e-09 & 3.645e-05 & 5.110e-05 \\
10 & (1,3) & 5.0 & 5.000000000117699 & 2.354e-11 & 5.199e-06 & 2.009e-05 \\
11 & (2,2) & 5.0 & 5.000000000746078 & 1.492e-10 & 1.821e-05 & 3.091e-05 \\
12 & (3,1) & 5.0 & 5.000000001093248 & 2.186e-10 & 1.973e-05 & 3.467e-05 \\
13 & (0,4) & 5.0 & 5.000000001562438 & 3.125e-10 & 2.651e-05 & 2.520e-05 \\
14 & (4,0) & 5.0 & 5.000000004861336 & 9.723e-10 & 4.059e-05 & 4.161e-05 \\
15 & (3,2) & 6.0 & 6.000000043151862 & 7.192e-09 & 3.095e-05 & 3.760e-05 \\ \hline \end{tabular}
\end{table}
Table 2: Errors of two-dimensional harmonic oscillator problem for the 16 lowest energy states.
harmonic oscillator

\[-\frac{1}{2}\Big{(}\nabla_{1}^{2}+\nabla_{2}^{2}\Big{)}\Psi+\frac{1}{2}(a_{11}x_{ 1}^{2}+2a_{12}x_{1}x_{2}+a_{22}x_{2}^{2})\Psi=E\Psi. \tag{4.16}\]

The coefficients \(a_{11},a_{12},a_{22}\) are chosen as in [25]. Since in [25] the coefficients are rounded to four decimal places, we use rounded values as exact coefficients, that is, \(a_{11}=0.8851\), \(a_{12}=-0.1382\), \(a_{22}=1.1933\). The exact eigenfunction has the following representation in coordinate \((y_{1},y_{2})\)

\[\Psi_{n_{1},n_{2}}(y_{1},y_{2})=\mathcal{H}_{n_{1}}(\mu_{1}^{1/4}y_{1})e^{-\mu _{1}^{1/2}y_{1}^{2}/2}\cdot\mathcal{H}_{n_{2}}(\mu_{2}^{1/4}y_{2})e^{-\mu_{2}^ {1/2}y_{2}^{2}/2}, \tag{4.17}\]

where \(y_{1}=-0.9339352418x_{1}-0.3574422527x_{2}\) and \(y_{2}=0.3574422527x_{1}-0.9339352418\), \(\mu_{1}=0.8322071257\), \(\mu_{2}=1.2461928742\), the two quantum numbers \(n_{1},n_{2}\) take the values \(0,1,2,\cdots\). The exact energy for the state \((n_{1},n_{2})\) is

\[E_{n_{1},n_{2}}=\Big{(}\frac{1}{2}+n_{1}\Big{)}\mu_{1}^{1/2}+\Big{(}\frac{1}{ 2}+n_{2}\Big{)}\mu_{2}^{1/2}. \tag{4.18}\]

We use 16 TNNs to learn the lowest 16 energy states. In each TNN, the rank is chosen to be \(p=20\), the subnetwork is built with depth 3 and width 50. The Adam optimizer is employed with a learning rate 0.001. We use the Adam optimizer in the first 500000 steps and then the L-BFGS in the subsequent 10000 steps. In each direction, 99 points Hermite-Gauss quadrature scheme are used to do the integration.

The corresponding numerical results are collected in Table 3, where we can find the proposed numerical method can also obtain the higher accuracy for the Schrodinger equation with coupled harmonic oscillator.

The corresponding approximate wavefunctions obtained by TNN-based machine learning method and the exact wavefunctions are shown in Figure 3, which implies that the approximate wavefunctions have very good accuracy even for the coupled harmonic oscillator.

#### 4.2.3 Five-dimensional coupled harmonic oscillator

Then, we investigate the performance of the proposed method for the case of five-dimensional coupled harmonic oscillator. Here, the Hamiltonian operator is defined as (4.8) with the matrix \(A\) being replaced with the following matrix

\[A=\begin{bmatrix}1.05886042&0.01365034&0.09163945&0.11975290&0.05625013\\ 0.01365034&1.09613742&0.10887930&0.07448974&0.07407652\\ 0.09163945&0.10887930&1.00935913&0.05588543&0.08968956\\ 0.11975290&0.07448974&0.05588543&1.17627129&0.06049045\\ 0.05625013&0.07407652&0.08968956&0.06049045&0.94969417\end{bmatrix}. \tag{4.19}\]

Then the exact energy is

\[E_{n_{1},n_{2},\cdots,n_{5}}=\sum_{i=1}^{5}\Big{(}\frac{1}{2}+n_{i}\Big{)}\mu_ {i}^{1/2}, \tag{4.20}\]
where \(\mu_{1}=0.88021303\), \(\mu_{2}=0.90973982\), \(\mu_{3}=1.02312382\), \(\mu_{4}=1.10243017\), \(\mu_{5}=1.37481559\), each of the five quantum numbers \(n_{i}\), \(i=1,2,\cdots,5\) takes the values \(0,1,2,\cdots\).

We use 16 TNNs to learn the lowest 16 energy states. A larger TNN structure than the 2-dimensional example is used, the rank is chosen to be \(p=50\), the subnetwork is built with depth 3 and width 100. The Adam optimizer is employed with a learning rate 0.001 and epochs of 500000. Then the final result is given by the subsequent 10000 steps LBFGS. The same 99 points Hermite-Gauss quadrature scheme as that in the last two examples is sufficient to this example.

Table 4 shows the corresponding numerical results, where we can find the proposed numerical method can obtain obviously better accuracy than that in [25].

Figure 4 shows the corresponding approximate wavefunctions obtained by TNN-based machine learning method. From Figure 4, the proposed numerical method here also has good accuracy for the five-dimensional eigenvalue problem of the coupled harmonic oscillator.

### Energy states of hydrogen atom

In this section, we study energy states of hydrogen atom. The wave function \(\Psi(x,y,z)\) of the hydrogen atom satisfies the following Schrodinger equation

\[-\frac{1}{2}\Delta\Psi-\frac{\Psi}{|\mathbf{r}|}=E\Psi, \tag{4.21}\]

\begin{table}
\begin{tabular}{c c c c c c c} \hline \(n\) & \((n_{1},n_{2})\) & Exact \(E_{n}\) & Approx \(E_{n}\) & err\({}_{E}\) & err\({}_{L^{2}}\) & err\({}_{H^{1}}\) \\ \hline
0 & (0,0) & 1.014291981649766 & 1.014291988589516 & 6.842e-09 & 2.801e-05 & 9.650e-05 \\
1 & (1,0) & 1.926545852963290 & 1.926545854461407 & 7.776e-10 & 1.264e-05 & 3.167e-05 \\
2 & (0,1) & 2.130622073635773 & 2.130622076362180 & 1.280e-09 & 1.766e-05 & 4.081e-05 \\
3 & (2,0) & 2.838799724276814 & 2.838799728095222 & 1.345e-09 & 2.099e-05 & 4.402e-05 \\
4 & (1,1) & 3.042875944949297 & 3.042875947694890 & 9.023e-10 & 1.633e-05 & 3.673e-05 \\
5 & (0,2) & 3.246952165621781 & 3.246952166999784 & 4.244e-10 & 1.171e-05 & 2.485e-05 \\
6 & (3,0) & 3.751053595590338 & 3.751053597870394 & 6.078e-10 & 1.500e-05 & 3.116e-05 \\
7 & (2,1) & 3.955129816262821 & 3.955129818022521 & 4.449e-10 & 1.289e-05 & 2.784e-05 \\
8 & (1,2) & 4.159206036935306 & 4.159206038606588 & 4.018e-10 & 1.281e-05 & 2.388e-05 \\
9 & (0,3) & 4.363282257607788 & 4.363282258514639 & 2.078e-10 & 9.008e-06 & 2.002e-05 \\
10 & (4,0) & 4.663307466903863 & 4.663307470243584 & 7.162e-10 & 1.856e-05 & 3.314e-05 \\
11 & (3,1) & 4.867383687576346 & 4.867383691191087 & 7.426e-10 & 1.987e-05 & 3.338e-05 \\
12 & (2,2) & 5.071459908248830 & 5.071459911990555 & 7.378e-10 & 1.992e-05 & 3.220e-05 \\
13 & (1,3) & 5.275536128921312 & 5.275536131659159 & 5.190e-10 & 1.741e-05 & 2.794e-05 \\
14 & (0,4) & 5.479612349593796 & 5.479612351630867 & 3.718e-10 & 1.421e-05 & 2.538e-05 \\
15 & (5,0) & 5.575561338217387 & 5.575561344662695 & 1.156e-09 & 2.769e-05 & 3.627e-05 \\ \hline \end{tabular}
\end{table}
Table 3: Errors of two-dimensional coupled harmonic oscillator problem for the 16 lowest energy states.
where \(|{\bf r}|=(x^{2}+y^{2}+z^{2})^{1/2}\). The exact energy of the hydrogen atom are \(E_{n}=-\frac{1}{2n^{2}}\) and there are \(n^{2}\) states consist with energy \(E_{n}\).

In order to compute the singular integrals of the Column potential terms \(1/|{\bf r}|\), we adopt spherical coordinates \((r,\theta,\varphi)\) with density \(r^{2}\sin\theta\). Then the wave function \(\Psi({\bf r})\) should be written as \(\Psi(r,\theta,\varphi)\). The Laplace \(\Delta\) has following expression

\[\Delta\Psi = \frac{\partial^{2}\Psi}{\partial r^{2}}+\frac{2}{r}\frac{\partial \Psi}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}\Psi}{\partial\theta^{2}}+ \frac{\cos\theta}{r^{2}\sin\theta}\frac{\partial\Psi}{\partial\theta}+\frac{ 1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}\Psi}{\partial\varphi^{2}} \tag{4.22}\] \[= \frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{ \partial\Psi}{\partial r}\right)+\frac{1}{r^{2}\sin\theta}\frac{\partial}{ \partial\theta}\left(\sin\theta\frac{\partial\Psi}{\partial\theta}\right)+ \frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}\Psi}{\partial\varphi^{2}}.\]

Figure 3: The contour plots of the first 16 eigenfunctions for two-dimensional coupled harmonic oscillator example in coordinate \((x_{1},x_{2})\). The two dashed lines are \(y_{1}=0\) and \(y_{2}=0\), respectively.
The 99 points Laguerre-Gauss quadrature is used in the direction \(r\) and 16 points Legendre-Gauss quadrature with subintervals length \(\frac{\pi}{64}\) in directions \(\theta,\varphi\). The TNN structure is defined as follows

\[\Psi(r,\theta,\varphi)=\sum_{j=1}^{p}c_{j}\phi_{r,j}(\beta r)e^{-\frac{\beta r}{ 2}}\cdot\phi_{\theta,j}(\theta)\cdot\big{(}\phi_{\varphi,j}(\varphi)\sin( \varphi/2)+\gamma_{j}\big{)}, \tag{4.23}\]

where \(\phi_{r}=(\phi_{r,1},\cdots,\phi_{r,p})\), \(\phi_{\theta}=(\phi_{\theta,1},\cdots,\phi_{\theta,p})\) and \(\phi_{\varphi}=(\phi_{\varphi,1},\cdots,\phi_{\varphi,p})\) are three FNNs with depth 3 and width 50, and \(p=20\). The activation function is selected as \(\sin(x)\). The trainable parameter \(\gamma_{j}\) is introduced to satisfy periodic boundary conditions \(\Psi(r,\theta,0)=\Psi(r,\theta,2\pi)\).

In implement, we use 15 TNNs to learn the lowest 15 energy states, each TNN is defined as (4.23). The Adam optimizer is employed with a learning rate 0.0003 and epochs of 100000 and then the L-BFGS in the subsequent 10000 steps to produce the final result. Table 5 shows the final energy approximations and corresponding errors.

## 5 Conclusions

In this paper, based on the deep Ritz method, we design a type of TNN-based machine learning method to compute the leading multi-eigenpairs of high dimensional eigenvalue problems. The most important advantage of TNN is that the high dimensional integrations of TNN functions can

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(n\) & \((n_{1},n_{2},n_{3},n_{4},n_{5})\) & Exact \(E_{n}\) & Approx \(E_{n}\) & err\({}_{E}\) \\ \hline
0 & (0,0,0,0,0) & 2.562993697776131 & 2.562993699775476 & 7.801e-10 \\
1 & (1,0,0,0,0) & 3.501190387362160 & 3.501190399601748 & 3.496e-09 \\
2 & (0,1,0,0,0) & 3.516796517763949 & 3.516796531293733 & 3.847e-09 \\
3 & (0,0,1,0) & 3.574489532179441 & 3.574489543933587 & 3.288e-09 \\
4 & (0,0,0,1,0) & 3.612960443213281 & 3.612960454781468 & 3.202e-09 \\
5 & (0,0,0,0,1) & 3.735519003914085 & 3.735519020687214 & 4.490e-09 \\
6 & (2,0,0,0,0) & 4.439387076948189 & 4.439387114214263 & 8.394e-09 \\
7 & (1,1,0,0,0) & 4.454993207349979 & 4.454993243128091 & 8.031e-09 \\
8 & (0,2,0,0,0) & 4.470599337751768 & 4.470599382176197 & 9.937e-09 \\
9 & (1,0,1,0,0) & 4.512686221765470 & 4.512686265259870 & 9.638e-09 \\
10 & (0,1,1,0,0) & 4.528292352167259 & 4.528292394259157 & 9.295e-09 \\
11 & (1,0,0,1,0) & 4.551157132799310 & 4.551157174501202 & 9.163e-09 \\
12 & (0,1,0,1,0) & 4.566763263201100 & 4.566763304302355 & 9.000e-09 \\
13 & (0,0,2,0,0) & 4.585985366582751 & 4.585985402475187 & 7.827e-09 \\
14 & (0,0,1,1,0) & 4.624456277616591 & 4.624456313384630 & 7.735e-09 \\
15 & (0,0,0,2,0) & 4.662927188650432 & 4.662927231227110 & 9.131e-09 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Errors of five-dimensional coupled harmonic oscillator problem for the 16 lowest energy states.
be calculated with high accuracy and efficiency. Based on the high accuracy and efficiency of the high dimensional integration, we can build the corresponding machine learning method for solving high dimensional problems with the high accuracy. The presented numerical examples show that the proposed machine learning method in this paper can obtain obviously better accuracy than the Monte-Carlo-based machine learning methods.

In our numerical implementation, we also find that the accuracy and stability of the machine learning process should be paid more attention. These are necessary to get the final high accuracy for solving high dimensional problems by using machine learning methods.

Figure 4: The contour plots of the first 16 eigenfunctions for five-dimensional coupled harmonic oscillator example.
Actually, the proposed TNN and the corresponding machine learning method can be extended to other high dimensional problems such as Schrodinger equations, Boltzmann equations, Fokker-Planck equations, stochastic equations, multiscale problems and so on. This means TNN-based machine learning method can bring more practical applications in physics, chemistry, biology, material science, engineering and so on. These will be our future work.
