A Deep Learning Method for Computing Eigenvalues of the Fractional Schrödinger Operator^∗

Fractional partial differential equations have proven effective in modeling anomalous diffusion phenomena^[1], turbulent flows^[2], porous media flows^[3] and many others. In the field of quantum mechanics, Laskin introduced the fractional Schrödinger equation^{[4, 5]} by using Levy flights to replace the classical Brownian motion. This equation has been used to reveal some novel phenomena^{[6, 7]}, which could not be explained by the standard Schrödinger equation.

Because of the non-locality and the singularity, most of the numerical methods for solving fractional partial differential equations focus on one-dimensional problems. For two-dimensional problems, the most popular method is an adaptive finite element method proposed by Ainsworth and Glusa^{[8, 9]}, while [10] also presents some modified finite difference methods. In addition, the authors in [11, 12] studied the spectral methods to solve the fractional partial differential equations, but applying them in irregular domains presents significant challenges. Another innovative method is the walk-on-spheres (WOS) method^{[13, 14]} and its extension^[15]. Those methods use the Feynman-Kac formula to calculate the solution at any given point by simulating a large number of paths of the $2s$-stable Lévy process. We refer to [16, 17, 18] for a comprehensive review of the numerical methods for fractional partial differential equations.

In this paper, we study the eigenvalue problem of the fractional Schrödinger operator in a bounded domain.

Here, $V(x)$ is a real-valued potential energy function. The precise definition of the operator $(-\Delta)^{s}$ with $s\in(0,1)$ will be presented in the subsequent section. When $V(x)=0$, the problem simplifies to the eigenvalue problem of the fractional Laplace operator:

To the best of our knowledge, the exact eigenvalue for any of the aforementioned problems is still unknown, and no analytical formulas even exist for calculating them under any specific circumstance. However, several numerical methods have been developed to tackle this challenge.

In [19, 20], the author presents a spectral method capable of computing very tight bounds for the eigenvalues of the fractional Laplacian within a unit ball in any dimension. Additionally, [21] introduces a spectral method that approximates the eigenvalues of the one-dimensional fractional Schrödinger operator. This method employs the same basis functions as [11] and can be extended to handle problems in hypercubes. Furthermore, a finite element method is displayed in [22] for approximating the eigenvalues within arbitrary bounded domains in one and two dimensions.

The use of neural networks for computing eigenvalues and eigenvectors dates back to the 1990s^{[23, 24]}. While the recent efforts focus on solving the eigenvalue problem of many-body quantum systems^{[25, 26, 27]}. These methods use neural networks to represent the underlying functions and incorporate techniques such as the variational Monte Carlo method and stochastic optimizers to approach the ground state of the quantum system. [28, 29, 30, 31, 32, 33] also exploit neural networks to solve eigenvalue problems. However, there is little literature on the numerical analysis of solving eigenvalue problems with neural networks except^[34].

In this work, we propose a deep learning method to solve the eigenvalue problem associated with the fractional Schrödinger operator. We convert the original eigenvalue problem into a sequence of minimization problems with constraints. By solving these minimization problems sequentially, we can determine the eigenmodes in ascending order of their eigenvalues. In addition, we introduce a novel loss function for solving these minimization problems. This loss function incorporates penalty terms to efficiently handle the orthogonal constraints. Moreover, we design a new neural network architecture that incorporates various feature functions. These feature functions are derived from prior knowledge of the underlying problem, particularly focusing on aspects such as singularity and boundary conditions.

To evaluate the accuracy of our method, we firstly solve a range of examples that could be addressed by the spectral methods also. The dimension of these problems varies from $1$ to $9$ and our network employs approximately 3,500 parameters. The relative error of our method is less than $0.1\%$ for the first $5$ eigenvalues. We further calculate up to $30$ eigenvalues while maintaining an error of less than $1\%$. Subsequently, we test our method in general domains where the spectral methods could not be used and compare the results with those obtained by the finite element method. The results demonstrate that our method produces more accurate results than the finite element method over the finest mesh.

We also implement our method for a fractional Schrödinger operator with an inverse square potential function in three dimensions. By computing the first $30$ eigenvalues, we observe that the order of the eigenvalues exchanges as the fractional order varies in this example. Additionally, we provide our estimations of eigenvalues for problems with different potential functions that have never been tested and exhibit eigenfunctions for various problems. All these examples collectively demonstrate the accuracy and efficiency of our method.

As an application, we apply the method to the fractional version of the isospectral problem. Based on our numerical result, we conjecture that even if the spectra of two domains are identical for the Laplacian, they would not be the same for the fractional Laplacian.

The rest of the paper is as follows. In Section 2, we begin with a discussion about the fractional operator and present the newly devised loss function. In Section 3, we display a deep learning scheme for solving the eigenvalue problem with a novel neural network architecture. We show the results of numerous numerical experiments in Section 4 and compare them with other existing methods. We apply our method to the fractional version of the isospectral problem in Section 5 and draw the conclusions in Section 6.

The fractional Laplace operator $(-\Delta)^{s}$ has many different definitions^{[17, 35]} and these definitions are not equivalent in all circumstances. In this paper, we adopt the Ritz definition as follows.

where

and $\|y-x\|$ represents the distance between $x$ and $y$, i.e., $\|y-x\|=\|y-x\|_{\ell_{2}}$. The integral in (3) should be interpreted in the principal value sense. This definition is also known as the integral definition in the literature and it is equivalent to the definition via the Fourier transform^[36]

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ represent the Fourier transform and the inverse Fourier transform, respectively.

To discuss the variational form of the eigenvalue problem (1), we first define the standard fractional order Sobolev space as

With the equivalent definition of the fractional Laplacian, it can be expressed as

where the norm is

with the seminorm is

For the eigenvalue problem with homogeneous Dirichlet boundary condition, we find a solution in

A bilinear form associated with the fractional Schrödinger operator is derived from the nonlocal Green’s first identity [16, (1.22)]

Using this bilinear form and the variational principle, the $k$th smallest eigenvalue is given by

where $E$ is a $k$-dimensional subspace of $H^{s}_{c}({\it\Omega})$. [37, 38, 39] display some analytical results for the eigenvalue of the fractional Schrödinger operator.

Numerically solving the min-max problem to derive the eigenvalues is challenging and requires significant effort. Therefore, we propose a new formulation that provides a more convenient approach to solve the problem by solely minimizing a loss function. Given the first $k$ eigenmodes $(\lambda_{1},u_{1}),\cdots,(\lambda_{k},u_{k})$, the subsequent eigenvalue can be characterized as

where

We use a neural network to approximate the next eigenfunction and construct a loss function by incorporating penalty terms to handle the orthogonal constraints. The loss function for computing the $k$th eigenvalue is

where $\beta$ is a penalty parameter and its value should be greater than $\lambda_{k}-\lambda_{1}$. When $k=1$, the loss function reduces to

and the penalty term becomes unnecessary. We remark that the loss function (15) originates from the variational principle and does not rely on any specific properties of the fractional differential operators. Thus, this loss function can be employed for the standard Schrödinger operator and other ordinary differential operators that admit Dirichlet forms.

We use the deep learning scheme introduced in the next section to minimize the loss. After minimizing the loss, the neural network provides an approximation of the eigenfunction $\mathaccent 866{u}_{k}$, while the eigenvalue $\lambda_{k}$ is approximated as

We iterate this process to obtain the eigenmodes one by one. Since the exact eigenfunction is unknown, the approximating function $\mathaccent 866{u}_{k}$ will be used as a substitute in (15). Consequently, any numerical errors from previous calculations can impact the accuracy of eigenmodes with higher eigenvalues. As the training progresses, the numerical error becomes larger and larger. However, subsequent experiments demonstrate that our method can calculate the first dozens of eigenmodes with great precision.

In this section, we present our deep learning scheme for solving the eigenvalue problem (1). We introduce a new architecture and describe the Monte Carlo sampling method used to calculate the loss.

When applying deep learning to solve PDE problems, the fully connected neural network (FCNN) and the residual network (ResNet) are commonly employed to approximate functions. However, these architectures have some weaknesses. They could not ensure the boundary conditions, and penalty terms are often added to the loss function to constrain the functions in general. Moreover, these architectures tend to underperform when the solutions contain singular terms. According to [40], solutions of fractional partial differential equations exhibit an $s$-order singularity near the boundary, i.e.,

where $v(x)$ is a smooth function over $\overline{{\it\Omega}}$. Furthermore, [21] also conjectures that the eigenfunctions of the fractional Schrödinger operator exhibit an $s$-order singularity near the boundary. To overcome these drawbacks, we design a new architecture that incorporates prior knowledge about the singular term of the functions and the boundary conditions.

Our architecture is based on FCNN, and a similar one is shown in [41]. It consists of several hidden layers and an output layer, while we mainly change the formulation of the output layer. Each hidden layer is a fully connected layer with the same width. We denote the number of hidden layers as $l$ and the width of each layer as $m$. The activation function we used is $\sigma=\tanh$ and we denote the vectorized function of $\sigma(x)$ as $\phi(x)$, i.e.,

The input layer contains a linear transformation that maps input values from $\mathbb{R}^{d}$ to $\mathbb{R}^{m}$, followed by an activation introducing non-linearity. The output of the first layer is given by

where $W_{1}\in\mathbb{R}^{m\times d}$ and $b_{1}\in\mathbb{R}^{m}$. The other hidden layers also contain similar transformations, mapping values from $\mathbb{R}^{m}$ to $\mathbb{R}^{m}$, and the output of the $i$th layer is represented as

where $W_{i}\in\mathbb{R}^{m\times m}$ and $b_{i}\in\mathbb{R}^{m}$. The output layer combines the outputs of the $l$th layer and a series of feature functions. The final output of the entire network is given by

where $W_{l+1}\in\mathbb{R}^{m}$, $W_{l+1}^{(j)}$ and $r_{l}^{(j)}$ is the $j$th component of $W_{l+1}$ and $r_{l}$, respectively. The feature functions are the key part of our method and play a critical role in ensuring the boundary conditions and capturing the singularity of the eigenfunctions. The selection of these functions is based on prior knowledge and significantly influences the accuracy and efficiency of our method. We enforce that the feature functions belong to $H^{s}_{c}({\it\Omega})$ to ensure the final output resides in the appropriate function space. The choice of the feature functions for different examples will be provided in the next section.

The set of all parameters is defined as

During each epoch, we calculate the loss function $L_{k}[u(x;\theta)]$ and employ stochastic optimization methods, such as the adaptive moment estimation (ADAM) optimizer, to update the parameters. After multiple epochs, the loss function will diminish to a value small enough and we will derive the approximated eigenmode. To balance accuracy and efficiency, we gradually reduce the learning rate and increase the number of sampling points as the training progresses.

Next, we present the details of the sampling technique used to estimate the loss function. The calculation of $(u,u_{j})$ and $\|u\|^{2}_{L^{2}({\it\Omega})}$ is straightforward. We uniformly sample $N$ points $x_{1},x_{2},\dots,x_{N}$ in ${\it\Omega}$ and calculate these values unbiasedly by

and

However, calculating the quadratic form $a(u,u)$ requires a more sophisticated approach since it is a double integral over $\mathbb{R}^{d}$ and the denominator becomes extremely small when the two variables are in close proximity. To address these difficulties, we pick a convex domain $D$ that contains ${\it\Omega}$ and separate the integral into two parts.

For the first part $A_{1}$, the integral domain is bounded, but its value approaches infinity as $y$ approaches $x$. We alleviate this by applying a coordinate transformation to change the integral formulation

Here, $S^{d-1}$ is the unit $(d-1)$-dimensional sphere and $w^{+}$ denotes the distance from the point $x$ to the boundary $\partial D$ along the direction $\xi$. The shape of $D$ should facilitate this calculation and its size should be as small as possible to enhance sampling efficiency. It is crucial that $D$ is a convex domain, which guarantees that any ray starting from $x$ intersects the boundary $\partial D$ only once. This makes the last equal sign in (26) and simplifies the calculation.

For the numerical implementation, we first uniformly sample $x_{1},x_{2},\dots,x_{N}$ in $D$, and uniformly sample $\xi_{1},\xi_{2},\dots,\xi_{N}$ in $S^{d-1}$. Then, we calculate $w^{+}_{i}$ which depends on $x_{i}$ and $\xi_{i}$, and sample $w_{i}$ with the probability

To reduce numerical instability caused by dividing a very small amount, we set a minimum value $w_{c}$ for $w$. In our numerical experiments, we take $w_{c}=10^{-4}$ and

Through these steps, we derive a practicable and efficient sampling method for calculating $A_{1}$,

The calculation of $A_{2}$ is comparatively simpler, as $u(y)=0$ when $y\in{\it\Omega}^{c}$. We simplify the equation by calculating part of the integral in advance and this prevents sampling in an infinite domain.

By employing the same method to uniformly sample $x_{i}$ and $\xi_{i}$ in $D$ and $S^{d-1}$, we can approximate $A_{2}$ unbiasedly by

In this section, we present the numerical results for various problems. Firstly, we introduce the training configuration. We employ the networks introduced in the last section. Each of them contains $3$ hidden layers with a width of $40$, unless otherwise stated. The total number of parameters is approximately $3,500$. For each eigenvalue, we conduct $120,000$ epochs. Initially, the learning rate is set to $5$e$\text{-}03$, and $1,000$ points are sampled. After every $20,000$ epochs, we reduce the learning rate to one-fourth of its current value, and double the number of sampling points. To identify the $k$th eigenvalue $(k\geq 2)$, we set the penalty parameter $\beta$ to $4$ times the maximum eigenvalue we have found. This selection is sufficiently large to discover the subsequent correct eigenvalue. Using a too large penalty parameter would cause the loss function to become sharp, leading to inefficient training. All experiments are conducted on a single NVIDIA A$100$ GPU.

We remark that there are several potential enhancements for the aforementioned setup, such as using larger and deeper networks, training more epochs and so on. However, implementing these enhancements would require much more effort and resources. The current configuration is chosen by balancing efficiency and accuracy. Generally, it takes approximately $5$ minutes to discover a new eigenvalue in most cases.

In this section, we explore the fractional order isospectral problem. In 1966, Kac posed the famous isospectral problem^[44], “Can one hear the shape of a drum”, which asks whether the Laplace operator with Dirichlet boundary conditions on two different domains can have the same spectrum? In 1992, Gordan, Webb, and Wolpert gave a negative answer to this question with a counterexample^{[45, 46]}, proving that it is possible for two domains to have the same spectrum. Figure 6 represents their counterexample. Since then, researchers have discovered numerous pairs of domains with identical spectra, and the eigenvalues of many of them are calculated by some numerical works^{[47, 48, 49]}.

Now, we wonder whether two different domains that have the same spectrum for the Laplace operator will also have the same spectrum for the fractional Laplace operator. We solve the eigenvalue problem of the fractional Laplace operator in these two domains to draw a conjecture to this question. The relative differences between two eigenvalues are calculated by

Here, $\lambda_{A,k}^{(s)}$ and $\lambda_{B,k}^{(s)}$ are the $k$th eigenvalues for the fractional Laplacian with order $s$ in the domains of drum A and drum B, respectively. The previous conclusion stated that

But based on the following experiments, we speculate that when $0<s<1$,

Since these two domains are not convex, we select a convex domain $D$ for sampling and the shape of the domain $D$ is plotted in Figure 6 also. We construct two types of feature functions similar to the previous example. The first-type feature functions are used to capture the singularity near the boundary while the second-type captures the singularity at the corners. The network has a width of $60$ and consists of $40$ first-type feature functions and $20$ second-type feature functions. All other settings remain the same as the previous examples.

Table 10 reports the first two eigenvalues in these two domains. It is evident that the first eigenvalue in the domain of drum A is smaller than that in drum B, whereas the second eigenvalue in drum A is greater than that in drum B.

We further compute more eigenvalues and calculate the relative difference between them. The results for different $k$ and $s$ are shown in Figure 7. The value $R_{k}^{(s)}$ for these $k$ and $s$ we calculated are significantly different from $0$ and the maximum relative difference reaches $1.8\%$. These discrepancies cannot be explained solely by the sampling error of the Monte Carlo method. Therefore, we conjecture that even if the spectra of two domains are identical when $s=1$, they would not be the same for $0<s<1$.