A clever neural network in solving inverse problems of Schrödinger
  equation
In this work, we solve inverse problems of nonlinear Schr\"{o}dinger
equations that can be formulated as a learning process of a special
convolutional neural network. Instead of attempting to approximate functions in
the inverse problems, we embed a library as a low dimensional manifold in the
network such that unknowns can be reduced to some scalars. The nonlinear
Schr\"{o}dinger equation (NLSE) is $i\frac{\partial \psi}{\partial
t}-\beta\frac{\partial^2 \psi}{\partial x^2}+\gamma|\psi|^2\psi+V(x)\psi=0,$
where the wave function $\psi(x,t)$ is the solution to the forward problem and
the potential $V(x)$ is the quantity of interest of the inverse problem. The
main contributions of this work come from two aspects. First, we construct a
special neural network directly from the Schr\"{o}dinger equation, which is
motivated by a splitting method. The physics behind the construction enhances
explainability of the neural network. The other part is using a library-search
algorithm to project the solution space of the inverse problem to a
lower-dimensional space. The way to seek the solution in a reduced
approximation space can be traced back to the compressed sensing theory. The
motivation of this part is to alleviate the training burden in estimating
functions. Instead, with a well-chosen library, one can greatly simplify the
training process. A brief analysis is given, which focuses on well-possedness
of some mentioned inverse problems and convergence of the neural network
approximation. To show the effectiveness of the proposed method, we explore in
some representative problems including simple equations and a couple equation.
The results can well verify the theory part. In the future, we can further
explore manifold learning to enhance the approximation effect of the
library-search algorithm.
Yiran Wang
# A clever neural network in solving inverse problems of Schrodinger equation

###### Abstract

In this work, we solve inverse problems of nonlinear Schrodinger equations that can be formulated as a learning process of a special convolutional neural network. Instead of attempting to approximate functions in the inverse problems, we embed a library as a low dimensional manifold in the network such that unknowns can be reduced to some scalars. The nonlinear Schrodinger equation (NLSE) is \(i\frac{\partial\psi}{\partial t}-\beta\frac{\partial^{2}\psi}{\partial x^{2}}+ \gamma|\psi|^{2}\psi+V(x)\psi=0\), where the wave function \(\psi(x,t)\) is the solution to the forward problem and the potential \(V(x)\) is the quantity of interest of the inverse problem. The main contributions of this work come from two aspects. First, we construct a special neural network directly from the Schrodinger equation, which is motivated by a splitting method. The physics behind the construction enhances explainability of the neural network. In particular, each convolution layer and activation function can be related to corresponding parts of the equation. Under this construction, it can be rigorously proved that the neural network has a convergence rate with respect to the length of input data and number of layers. The other part is using a library-search algorithm to project the solution space of the inverse problem to a lower-dimensional space. The way to seek the solution in a reduced approximation space can be traced back to the compressed sensing theory. The motivation of this part is to alleviate the training burden in estimating functions. Instead, with a well-chosen library, one can greatly simplify the training process. More specifically, in one of the experiments, we analysis the landscape of the loss function with respect to the training parameters to help us obtain the optimal solution to the inverse problem. A brief analysis is given, which focuses on well-possedness of some mentioned inverse problems and convergence of the neural network approximation. To show the effectiveness of the proposed method, we explore in some representative problems including simple equations and a couple equation. The results can well verify the theory part. In the future, we can further explore manifold learning to enhance the approximation effect of the library-search algorithm.

_Keywords--_ Physics-inspired neural network, Reduced order modeling, Inverse problems, Schrodinger equation, Compressed sensing

## 1 Introduction

Neural network approach has been widely explored in solving partial differential equations (PDE) since the pioneer work [30]. Most of them choose from some existing architectures and train the model to obtain satisfying results for the PDEs. The advantage is that it does not require too much information from the equation itself especially when it is a data-driving method. This provides great convenience in solving different types of equations hence it is user-friendly, which also brings some drawback on the other hand. Performing as a blackbox, this kind of neural network approach lacks sufficient explainability, which comes from three aspects. The first one is the existence of parameters that are optimal to some sense from theoretical aspect. This is somehow answered by the popular universal approximation theorem for neural
networks. However, it is not sufficient to choose only one single layer in many cases. For a network with a particular architecture, the existence of optimal parameters is unclear. Secondly, even though one can obtain the optimal solution from theory, it remains unclear whether the optimal solutions can be attained with the training process. In many cases, one can only find it with trial and error. The third uncertainty exists in the generalization. If one can obtain a set of "good" parameters to build an optimal neural network, whether the performance is satisfying when different dataset are used. If it performs well in some similar data as the training data, which function space it can approximate. Overall, there is not a general explanation about the above three problems, i.e., the existence of optimal parameters, how to train it, and the generalization power.

To enhance the explainability, one can equip the construction of neural network with information of equations instead of treating networks and PDEs separately. To this end, physics-informed neural newt (PINN) [46] serves as a representative example by utilizing a loss function derived from equation. Inspired by this work, we reconsider the relationship of neural network and PDEs and desire to figure out a novel structure of network that makes use of the equation itself. As for the interested equations, we focus on a class of nonlinear Schrodinger equation since they are significant in both scientific research and industry [24, 39, 28, 50]. Moreover, numerical treatments of Schrodinger equation have been investigated extensively [5, 15, 53, 36, 7], among which time splitting schemes are widely used [8, 9, 55, 14, 54]. Apart from forward process, inverse problems attract attention from many researchers due to its wide applications [38] and contribute some effective numerical techniques [35, 4, 48, 1, 52]. Based on the understanding of the time splitting method, we consider designing a novel network from the nonlinear Schrodinger equation. Hence, we call the derived network NLS-Net. In [59], we inverstigated an coefficient inverse problem in nonlinear Schrodinger equation without potential. In this work, we explore an inverse problem of finding potential from data with the NLS-Net.

Compared with scalars, it is far more complicated to seek functions. To this end, we introduce a library-search approach to reduce the problem to a coefficient inverse problem. This method is motivated from an algorithm in sparse identification of nonlinear dynamics (SINDy-PI) [27]. The priciple idea is to choose some functions from a reference set, which can be understood as a library. We further use a linear combination of the chosen functions to represent the unknown potential to the inverse problem. In other words, we project the original solution to a fixed approximation space, hence we only need to determine the corresponding coefficients. However, from the perspective of expressive power, we need to incoporate sufficient modes in the library, which performs similarly as utilizing an either deep or wide neural network. Put differently, it is obvious that approximation power comes from sufficient degree of freedom. On the other hand, increasing the number of unknowns means a growth in computation cost, which is apparently a side effect. For example, everyone in the area of scientific computing is familiar with the word "the curse of dimensionality" [29], which may bring a blow up in computation cost. Hence, model reduction methods [2, 12, 44, 11, 43, 45] are introduced to seek a "cheaper" model based on the original model to attain the quantity of interest. In particular, there are considerable number of works in reduced order methods of solving PDEs numerically [42, 13, 60, 3, 40, 57, 47, 23]. For inverse problems in PDEs, many researchers devote good ideas of applying model reuction [19, 16, 21, 32, 20, 22]. It is understandable that one needs to sacrifice accuracy in order to reduce the computation cost since there is no free lunch. However if one has prior information about the quantity of interest, it is possible to use some well-designed models that can attain a relatively high accuracy with reduced workload. In image or signal processing, compressed sensing [37, 17, 56, 18, 34, 51, 49] is utilized for enforcing sparsity in some incorporated variables based on the knowledge that the solution lies in a low-dimensional subspace. Therefore, it is reasonable to apply this technique in the area of numerical PDE as a other reduced order approach, which becomes another motivation of this work where we desire to approximate the unknown potential in the nonlinear Schrodinger equation. The prior information is that the potential function should have a relatively simple expression, which can be a linear combination of a few different basis functions. Hence, it is reasonable to assume the target coefficient vectors with respect to bases exploits a sparse form, which explains the motivation of applying compressed sensing in this context. Overall, contributed by above motivations, we propose a library-search approach combined with compressed sensing to train a special network called NLS-Net that is derived from NLSE. The proposed method is potential to be developed a new framework of solving inverse problems. To evaluate the performance of the proposed method, we use three reprentative examples, two of which are single equations with unknown potentials while the rest one is coupled system of two NLSE. In all the examples, we use NLS-Net to obtain the forward solutions. From [6], one can derive an upper bound for the approximation errors of forward solutions, which can in turn guarantee the convergence of
NLS-Net with respect to number of layers and data amount. In our numerical results, we show that the performance of the NLS-Net is consistent to the theory. For the inverse part, we use the library-search combined with compressed sensing to seek unknown potentials. In view of the specificity of the third case, we apply a landscape analysis to help us find the optimal solution more efficiently. It is satisfying that the concerned loss as well as errors present an exponential decay in all the three examples. Moreover, we compare the numerical solutions of the inverse part, which are approximately identical to the groundtruth. Both of them shows the effect and accuracy of the proposed method.

The paper is organized in six parts. We present the concerned problem in the first part and review the conventional split-step methods in the next part. The proposed method is introduced in the third section and it is compsed of four subsections. We first introduce the mechanism of designing a convolutional neural network from NLSE. Following this, we provide some background in compressed sensing. The library-search algorithm is then proposed, which is combined with compressed sensing to form the final approach in the last part. We alsp present some relevant theroy in the fifth section. Numerical simulations are demonstated in the following part. We conclude this work and offer some remark and potential interested problems in the last place.

## 2 Problem description

From an initial value problem of Schrodinger equation as follows:

\[i\frac{\partial\psi}{\partial t}+\beta\frac{\partial^{2}\psi}{ \partial x^{2}}+\gamma|\psi|^{2}\psi+V(x)\psi =0, \tag{1}\] \[\psi(0,x) =h(x),\]

where \(\psi(x,t):\mathbb{R}\times[0,T]\rightarrow\mathbb{C}\) is the wave function and \(V(x)\) is the potential to be estimated. Moreover, \(\beta\),\(\gamma\) are constants and \(h(x):\mathbb{R}\rightarrow\mathbb{C}\) is a given as the initial condition. The observed dataset composed of continuous functions is defined to be \(\mathcal{D}_{\text{cont}}:=\{\psi_{0}(x),\psi_{T}(x)\}\) where \(\psi_{0}(x)\) and \(\psi_{T}(x)\) is the observed wave at the initial time \(0\) and final time \(T\), respectively. Based on the observed data, we estimate the potential \(V(x)\). In particular, we define the \(S(t;V)\) to be Schrodinger operator in (1) as \(\psi(\cdot,t)=S(t;V)\psi(0,\cdot)\). To identify the potential \(V\) from \(\mathcal{D}_{\text{cont}}\), we formulate the following optimization problem:

\[\min_{V}\|S(t;V)\psi_{0}-\psi_{T}\|_{2}^{2}. \tag{2}\]

Numerically, we can not use the exact \(S(t;V)\) operator that requires directly solving (1) which is infeasible. Instead, we use operator splitting [25, 31, 26] as an approximation method.

## 3 Conventional split-step methods

In this part, we review some theoretical results on the well-posedness and stability of the inverse problem. Besides, we state the conventional split-step method. First, we rewrite (1) into the form

\[\frac{\partial}{\partial t}\psi(x,t)=(\mathcal{L}+\mathcal{N}) \psi(x,t),\quad x\in\mathbb{R}. \tag{3}\]

Here \(\mathcal{L}:=-i\beta\frac{\partial}{\partial x^{2}}\) is the linear operator and \(\mathcal{N}:=i\gamma|\psi(x,t)|^{2}+iV(x,t)\) is the nonlinear operator. Instead of directly solving (3), one can split it into two subproblems corresponding to two operators \(\mathcal{L}\) and \(\mathcal{N}\), respectively. In particular, one can solve a linear problem

\[\frac{\partial\psi}{\partial t}=\mathcal{L}\psi_{L}(x,t), \tag{4}\]

and a nonlinear problem

\[\frac{\partial\psi}{\partial t}=\mathcal{N}\psi_{N}(x,t). \tag{5}\]

Both of the above two problems can be solved sequentially to obtain a final solution. We denote \(\delta t\) to be the time step size. In terms of performing the splitting, there are different orders of the method. The simplest version is the first-order method, in which is each advancement of \(\Delta\) is composed of two steps.
First, one solves either the linear (4) or the nonlinear problem (5) to generate a intermediate solution. Then one apply the other operator to finish this time block. More precisely, if one choose \(\mathcal{L}\) prior to \(\mathcal{N}\) for each \(\Delta t\), then a direct result after forwarding one time block is

\[\psi(x,t+\Delta t)\approx\exp(\mathcal{N}\Delta t)\exp(\mathcal{L} \Delta t)\psi(x,t).\]

Similarly, the second-order method also be decoupled into alternatively solving linear (4) and (5). The only difference is we split each \(\Delta\) into three steps. The key idea is that the solution of one step serves as the initial condition for the next step. The most popular second-order splitting is the standard Strang splitting with the form

\[\psi(x,t+\Delta t)\approx\exp(\mathcal{L}\Delta t/2)\exp(\mathcal{ N}\Delta t)\exp(\mathcal{L}\Delta t/2)\psi(x,t).\]

The advantage of the splitting method is that both (4) and (5) can be much easier than (3). In particular, one can explicitly solve (5) and for (4), many choices can be chosen, for example, finite difference schemes or fourier method. The final combination of splitting and spatial discretization ways form the different numerical methods like split-step finite difference (SSFD) and split-step fourier method (SSFM).

## 4 A physics-inspired neural network with a library-search approach

Compared the last section, we now introduce an new approach with a special-designed neural network (NN) [59] and a library-search way to solve some interested inverse problems. First of all, we first give an introduction of the special NN. Then, we introduce the library-search approach which serves as a low-rank representation of the solution to the inverse problem. Finally, we summarize the overall method in the last subsection.

### Learning process of a convolutional neural network

Motivated by SSFM, we can formulate a convolutional neural network to solve the inverse problem. In particular, for linear and nonlinear operator in (4) and (5), we use a convolutional operator and a nonlinear operator, whose can be explicitly derived from the equation. As the solution of the inverse problem is incorporated into the NN layers, the backpropagation essentially serves as a process of the seeking the ground truth of the inverse problem. We define the linear and nonlinear part as follows.

For linear problem, we apply

\[\psi_{L}(x,t+\Delta t)=\kappa\int K(s;\beta\Delta t)\psi_{L}(x-s,t), \tag{6}\]

where \(K(s;\eta)\) is a convolution kernel defined by

\[K(s;\eta):=\sqrt{\frac{i}{\pi\eta}}\exp\left(-\frac{is^{2}}{4 \eta}\right),\forall s\in\mathbb{R}.\]

Furthermore, we can introduce a convolution operator \(\mathcal{C}(\cdot;\beta\Delta t)\) such that

\[\psi_{L}(\cdot,t+\Delta t)=\mathcal{C}(\cdot;\beta\Delta t)\psi( \cdot,t):=K(s;\beta\Delta t)*\psi_{L}(\cdot,t),, \tag{7}\]

where \(*\) refers to convolution with respect to \(x\). For nonlinear problem, we use

\[\psi_{N}(x,t+\Delta t)=\mathcal{G}(\psi_{N};\gamma\Delta t,V\Delta t), \tag{8}\]

where the nonlinear operator \(\mathcal{G}\) is defined as

\[\mathcal{G}(u;g,h):=u\exp i(g|u|^{2}+h).\]

Finally we formulate the forward propagation as

\[\psi(\cdot,t+N\Delta t)\approx\mathcal{S}(\cdot;\beta,\gamma,V, \Delta t)\psi(\cdot,t):=\left(\mathcal{C}\left(\cdot;\frac{\beta\Delta t}{2} \right)\circ\mathcal{G}(\cdot;\gamma\Delta t,V\Delta t)\circ\mathcal{C}\left( \cdot;\frac{\beta\Delta t}{2}\right)\right)^{\diamond N}\psi(\cdot,t) \tag{9}\]
Corresponding to the continuous version of linear and nonlinear operators in (6) and (8), we further define the discrete version of \(\mathcal{L}\) and \(\mathcal{N}\). Before this, we need to consider a corresponding equtaion with periodic boundary condition as follows.

\[i\frac{\partial\psi}{\partial t}+\beta\frac{\partial^{2}\psi}{ \partial x^{2}}+\gamma|\psi|^{2}\psi+V(x)\psi =0,\quad x\in\Omega,t\in[0,T], \tag{10}\] \[\psi(x,0) =h(x),\quad x\in\Omega,\] \[\psi(a,t) =\psi(b,t),\quad t\in[0,T]\] \[\frac{\partial\psi}{\partial x}(a,t) =\frac{\partial\psi}{\partial x}(b,t),\quad t\in[0,T],\]

where \(\Omega\) is specified as \([a,b]\). We choose homogeneous grid \(\{a_{j}\}_{j=1}^{M}\) with \(a_{j}=a+\frac{j(b-a)}{M}\) for \(j=1,\ldots,M\). We define \(\Delta x=\frac{b-a}{M}\). Since we enforce the periodic boundary condition, i.e., \(\psi(t,a)=\psi(t,b)\), we choose \(b\) in the grid and \(a\) is not in the grid. Then, we define a discrete convolution kernel as

\[\tilde{K}(q;\eta)[j]:=\sqrt{\frac{i}{\pi\eta}}\exp\left(-\frac{i(jq)^{2}}{4 \eta}\right),\quad\forall j\in[M],\]

where \(M\) is the number of collocation points in the spatial domain and we use a combinatorial notation \([M]=1,\ldots,M\). We define \(t_{n}=n\Delta t\). With the spatial discretization into the grid \(\{a_{j}\}_{j=1}^{M}\), we define \(\tilde{\psi}_{M,L}^{m}\) to be the approximate solution to (6) at time \(t_{m}\), from which we update the solution at \(t_{m+1}\) with discrete convolution as follows,

\[\psi_{L}(x_{n},t_{m+1})\approx\tilde{\psi}_{M,L}^{m+1}[n]=\tilde{K}(q;\beta \Delta t)*\tilde{\psi}_{M,L}^{m}:=\sum_{j=1}^{M}\tilde{K}(q;\beta\Delta t)[j] \tilde{\psi}_{M,L}^{m}[n-j],\]

where \(\tilde{K}(q;\eta)\) are replaced by zero when \(n-j\notin[M]\). Again, we define an associated discrete convolution operator \(\tilde{\mathcal{C}}(\cdot;\eta)u:=\tilde{K}(q;\eta)*u\). As for the nonlinear operator, we use \(\tilde{\mathcal{G}}\) for the discrete version but the definition of \(\tilde{\mathcal{G}}\) is the same as \(\mathcal{G}\) but we can restrict the domain to set of discrete vector. Corresponding to (9), we can summarize the discrete version by combining the \(\mathcal{C}\) and \(\mathcal{G}\) as

\[\tilde{\psi}_{M}^{N}=\tilde{\mathcal{S}}(\cdot;\beta,\gamma,V,\Delta t)\tilde {\psi}_{M}^{0}:=\left(\tilde{\mathcal{C}}\left(\cdot;\frac{\beta\Delta t}{2} \right)\circ\tilde{\mathcal{G}}(\cdot;\gamma\Delta t,V\Delta t)\circ\tilde{ \mathcal{C}}\left(\cdot;\frac{\beta\Delta t}{2}\right)\right)^{\odot N}\tilde{ \psi}_{M}^{0}, \tag{11}\]

where we define \(\tilde{\psi}_{M}^{N}\) as the discrete solution to (10) with the operator splitting method. This can be expressed as a \(N-\) layer complex-valued CNN with activation function \(\tilde{\mathcal{G}}\) as below:

\[\tilde{\mathcal{W}}(\Theta)=\widetilde{C}_{h}\circ\mathcal{G}\circ\tilde{ \mathcal{C}}\circ\ldots\circ\tilde{\mathcal{C}}\circ\mathcal{G}\circ\widetilde {\mathcal{C}}_{h}. \tag{12}\]

Here \(\Theta\) is the set of weights contained in the operators \(\tilde{\mathcal{C}}\) and \(\tilde{\mathcal{G}}\) and \(\widetilde{C}_{h}\) refers to \(\frac{\Delta t}{2}\) of forwarding in the linear step.

### Compressed sensing background

In [17], the problem is to reconstruct unknown image or signal \(x\in\mathbb{R}^{m}\) from \(n\) measurements. Based on the prior information that \(x\) has a sparse representation in some orthonormal basis, the number of measurements \(n\) can be dramatically smaller than the size \(m\). For example, some natural classes of images with \(m\) pixels need only \(n=O(m^{1/4}\log^{\frac{5}{4}}(m))\) samples for relatively accurate recovery. Here, we briefly state the problem and focus on the information that can be used in the proposed method.

Suppose \(X\) is a class of objects of interest and one needs to recover \(X\) from some measurements. First they define an information operator \(I_{n}:X\rightarrow\mathbb{R}^{n}\) as follows:

\[I_{n}(x)=(\langle\xi_{1},x\rangle,\ldots,\langle\xi_{n},x\rangle)\]

where \(\xi_{i}\) are sampling kernels. Furthermore, \(A_{n}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) is an algorithm that reconstruct \(x\) from \(I_{n}(x)\). Overall, the problem is formulated as minimizing the \(l_{2}\) error of the reconstruction

\[E_{n}(X)=\inf_{A_{n},I_{n}}\sup_{x\in X}\|x-A_{n}(I_{n}(x))\|_{2}.\]
Here, \(A_{n}\), \(I_{n}\) can be any possible methods of linear sampling and reconstruction. The goal is to obtain a faithful recovery with limited measurements. In particular, it can be split into two questions: how to choose \(I_{n}\) and once it is fixed, how to determine an effective way to reconstruct \(x\), i.e., finding \(A_{n}\). First of all, the following assumption is used.

**Assumption 1**.: _For each \(x\in X\), there exists a set of orthogonal bases \((\psi_{i}\in\mathbb{R}^{m},i=1,\ldots,m)\). Define \(\theta_{i}:=\langle x,\psi_{i}\rangle\) to be the corresponding transform coefficients as \(\theta_{i}:=\langle x,\psi_{i}\rangle\). Then \(x\) is sparse in the sense that for some \(0<p<2\) that_

\[\|\theta\|_{p}\equiv\left(\sum_{i}|\theta_{i}|^{p}\right)^{\frac{1}{p}}\leq R\]

_for some \(R>0\)._

The authors are only interested in near-optimal performance, i.e., optimal up to some constants. In [17], they provided the condition of near-optimal information operator and the reconstruction procedure based on linear programming. Briefly speaking, \(\Phi\) is an \(n\times m\) matrix generated by randomly sampling the columns and \(J\subset\{1,\ldots,m\}\) is a set of index. Let \(\Phi_{j}\) denote a submatrix of \(\Phi\), which is a result of choosing indicated columns of \(\Phi\). This matrix contains the information of sampling. We pay special attention to one condition of \(\Phi_{j}\) for the desired approximation power: for strict positive parameter \(\rho\), the minimal singular value of \(\Phi_{J}\) exceeds \(\eta_{1}>0\) uniformly in \(|J|<\frac{\rho n}{\log(m)}\). Here, the columns of \(\Phi_{j}\) are used as bases functions. To attain an effective approximation, it is significant to obtain a particular quantitative degree of linear independence among bases. Based on the information matrix \(I_{n}\), one can solve the following convex optimization problem:

\[\min\|\Psi^{T}\|_{1}\text{ subject to }y_{n}=\Phi\Psi^{T}(x),\]

where \(\Psi\in\mathbb{R}^{m\times n}\) is an orthorgonal matrix. In other words, a desired form is to find \(x\) with smallest \(l_{1}\) norm of transform coefficients and consistent with the measurements \(y_{n}\).

### A library-search algorithm

Motivated by the [27], we can formulate a library-search algorithm in solving the inverse problem. Suppose we have a set of functions denoted by \(\mathcal{P}=\{g_{1},\ldots,g_{N_{L}}\}\), where \(N_{L}\) is the cardinality of the function set \(\Phi\). The key idea of this library-search method is that we solve \(V\) in a space spanned by \(\mathcal{S}\). Mathematically, we can formulate the process as follows. Given a set of collocation points \(\{x_{1},\ldots,x_{N_{x}}\}\), we define vectors \(g_{i,d}\in\mathbb{R}^{N_{x}}\) by \(g_{i,d}[j]:=g_{i}(x_{j})\), \(\forall j\in[N_{x}]\), \(\forall i\in[N_{L}]\). The discretized target solution \(\psi_{T}\) can be similarly defined as \(\psi_{T,d}[j]=\psi_{T}[x_{j}]\), \(\forall j\in[N_{x}]\). Now we have obtained a set of bases \(\{g_{i,d}\}\), which can generate a subspace \(H\) in \(\mathbb{R}^{N_{x}}\). We further assemble a matrix \(\Phi\in\mathbb{R}^{N_{x}\times N_{L}}\) define by \(\Phi e_{j}=g_{i,d}\), where \(e_{j}\) are the \(j\)-th
standard basis for \(\mathbb{R}^{N_{L}}\). Given a coordinate vector \(\mathbf{c}\in\mathbb{R}^{N_{L}}\), one can generate a vector \(y\in H\) defined by \(y=\Phi\mathbf{c}\). Now we replace (2) by the following problem:

\[\min_{\mathbf{c}}\|\tilde{\mathcal{S}}(\cdot;\beta,\gamma,\Phi\mathbf{c},\Delta t )\tilde{\psi}_{L}^{0}-\psi_{T,d}\|_{2}^{2}. \tag{13}\]

### Our approach

From first subsection, we explore some details of the construction of the convolutional neural network derived from the Schrodinger equation, which is therefore denoted by NLS-Net. To represent the potential, we introduce the library-search algorithm, which indeed serves as an effective way of parameterizing the NLS-Net. Inspired by the compressed sensing, we can seek an approximation with relatively high sparsity. With the notation in the compressed sensing, we desire to find a way to find a corresponding transform coefficients with relatively small \(l_{1}\) norm.

For the training process of NLS-Net, we use a summation of \(l^{2}\) error and \(l^{1}\) norm of \(\mathbf{c}\) for the loss function defined below:

\[J_{\mathbf{c}}=\frac{1}{2}\frac{\|\tilde{\mathcal{S}}(\cdot;\beta,\gamma, \Phi\mathbf{c},\Delta t)\tilde{\psi}_{L}^{0}-\psi_{T,d}\|_{2}^{2}}{\|\psi_{T, d}\|_{2}^{2}}+\lambda\|c\|_{1}. \tag{14}\]

To minimize \(J_{\mathbf{c}}\), we use proximal gradient descent [41]. We summarize the method of combining the NLS-Net with the library search algorithm in the following Table.

## 5 Theoretical background

In this section, we present the theoretical background of the inverse problem. We briefly mention two parts as follows:

1. The well-posedness of the inverse problem.
2. The convergence analysis of the NLS-Net

In particular, in the first part, we review some important results about uniqueness and stability of a special form of inverse problem. We remark that the existing result is limited to linear case however, it offers insights of the nonlinear case which deserves exploration. For the second part, we present the convergence rate of the constructed neural network with respect to the number of collocation points and layers.

In [10], Baudouin et al. (2007) proved uniqueness and stability in an inverse problem for the Schrodinger equation. They consider the following equation:

\[\left\{\begin{array}{l}iy^{\prime}(x,t)+\Delta y(x,t)+q(x)y(x,t)=0,x\in \Omega,t\in(0,T)\\ y(x,t)=h(x,t),x\in\partial\Omega,t\in(0,T)\\ y(x,0)=y_{0}(x),x\in\Omega.\end{array}\right. \tag{15}\]

For this equation, they proved the following theorem.

**Theorem 2**.: _Let \(\mathcal{U}\) be a bounded subset of \(L^{\infty}(\Omega)\), \(q\in L^{\infty}(\Omega)\) and \(y\) be a solution of (15). Assuming_

\[\exists x_{0}\in\mathbb{R}^{N}\backslash\bar{\Omega}\text{ such that }\Gamma_{0}\supset\{x\in\partial\Omega;(x-x_{0})\cdot\nu(x)\geq 0\}\] \[y(q)\in W^{1,2}\left(0,T,L^{\infty}(\Omega)\right)\]
\(y_{0}\) is real valued and \(|y_{0}|\geq r_{0}>0\), as in \(\bar{\Omega}\). There exists a constant \(C=C(\Omega,T,\Gamma_{0},\|q\|_{L^{\infty}(\Omega)},y_{0},h,\mathcal{U})>0\) such that if_

\[\frac{\partial y(q)}{\partial\nu}-\frac{\partial y(p)}{\partial\nu}\in H^{1}(0, T;L^{2}(\Gamma_{0}))\]

_then \(\forall p\in\mathcal{U}\),_

\[\|q-p\|_{L^{2}(\Omega)}\leq C\left\|\frac{\partial y(q)}{\partial\nu}-\frac{ \partial y(p)}{\partial\nu}\right\|_{H^{1}(0,T;L^{2}(\Gamma_{0}))}.\]

This theorem characterizes the fact that the solutions to the inverse problem change continuously with the normal derivative of the forward solutions.

In the following part, we review some important results of the convergence rate of the split-step method, which is the core of the NLS-Net. In [33], Lubich (2008) proved the convergence rate of the splitting methods for Schrodinger-Poisson and cubic nonlinear Schrodinger equations. Here we only consider the latter one as follows:

\[i\frac{\partial\psi}{\partial t}=-\nabla\psi+V\psi,\quad x\in\mathbb{R}^{3}, \quad t>0, \tag{16}\]

where \(V=V[\psi]=\pm|\psi|^{2}.\) Assuming the solution \(\psi(t)\) to the (16) is in \(H^{4}\) and let \(m=\max_{0\leq t\leq T}\|\psi(t)\|_{H^{4}}\), the following theorem is proved.

**Theorem 3**.: _The numerical solution \(\psi_{n}\) given by the strang-split scheme with step size \(\tau\) for equation (16) has a first-order error bound in \(H^{2}\) and a second-order error bound in \(L^{2}\), for \(t_{n}=n\tau\leq T\),_

\[\|\psi_{n}-\psi(t_{n})\|_{H^{2}}\leq C(m,T)\tau,\] \[\|\psi_{n}-\psi(t_{n})\|_{L^{2}}\leq C(m,T)\tau^{2}.\]

This theorem provides the convergence rate of strang-split scheme with respect to the time step \(\tau\) under \(H^{2}\) and \(L^{2}\) norm. However, it is a continuous case in the spatial domain and they didn't consider a discretized case. Hence, we proceed to present the result given in [6] as follows.

**Theorem 4**.: _Let \(\psi^{n}\in X_{M}\) be the numerical approximation obtained by the TSSP [6]. Under assumptions (4.84) and (4.85), there exist constants \(0<\tau_{0},h_{0}\leq 1\), such that if \(0<h\leq h_{0},0<\tau\leq\tau_{0}\) and \(m\geq 5\), we have_

\[\|\psi\left(x,t_{n}\right)-I_{M}\left(\psi^{n}\right)(x)\|_{L^{2}(U)}\lesssim h ^{m}+\tau^{2},\quad\|\psi^{n}\|_{\infty}\leq M_{1}+1,\]

\[\|\nabla\left(\psi\left(x,t_{n}\right)-I_{M}\left(\psi^{n}\right)(x)\right)\| _{L^{2}(U)}\lesssim h^{m-1}+\tau^{2},\quad 0\leq n\leq\frac{T}{\tau},\]

_where the interpolation operator \(I_{M}\) is given below and \(M_{1}=\max_{t\in[0,T]}\|\psi(\cdot,t)\|_{L^{\infty}}\)._

The mentioned notations are defined as follows.

1. \(X_{M}=\{v=(v_{j})_{j\in\mathcal{T}_{M}^{0}}|v_{0}=v_{M}=0\}\subset C^{M+1}\) with \(\mathcal{T}_{M}^{0}=\{j|j=0,1,\ldots,M\}\).
2. TSSP refers to time-splitting sine pseudospectral method introduced in [6]. The difference of this method from split step fourier method is that the spatial discretization is sine spectral method while the latter one is fourier method.
3. \(I_{M}\) is an interpolation operator defined as \[(I_{M}\phi)(x)=\sum\limits_{l=1}^{M-1}\tilde{\phi}_{l}\sin(\mu_{l}(x-a)),\] where \(\tilde{\phi}_{l}=\frac{2}{M}\sum\limits_{j=1}^{M-1}\phi_{j}\sin(\frac{d\pi}{M })\) with \(l\in\mathcal{T}_{M}\) and \(\phi_{j}=\phi(x_{j})\).

Motivated from this theorem, we can have the following upper bound for NLS-Net.

**Theorem 5**.: _We define that \(\tilde{\psi}_{M}^{n}\) is the solution at time \(t_{n}=n\Delta t\) generated by (11) with \(M\) as the input data length and \(N\) is the number of layers. We recall that \(\psi(x,t_{n})\) is the analytical solution to (1) at time \(t_{n}\). Suppose that the initial value \(\psi(\cdot,0)\in H^{s}(\mathbb{R}\to\mathbb{C})\) for \(s\geq 5\). Then there exist parameters \(\Theta\) of the neural network \(\tilde{\mathcal{W}}(\Theta)\), which is defined by (12) such that_

\[\|\tilde{\mathcal{W}}(\Theta)\tilde{\psi}_{M}^{n}-\psi_{T,d}\|\lesssim(\Delta x )^{s}+(\Delta t)^{2}, \tag{17}\]

_where \(\Delta x=\frac{b-a}{M}\) and \(\Delta t=\frac{T}{N}\)._
We remark that this theorem specifies that the approximation effect of the neural network converges respect to the spatial and time step sizes.

## 6 Numerical experiments

In this section, we demonstrate the performance of the proposed method with three examples. In the first example, we use the NLS-Net to approximate the one-soliton solution for the generalized nonlinear Schrodinger equation (NLSE), where the potential is exponential. For the second one, we further consider another form of NLSE with a periodic potential, which is called Gross-Pitaevskii (GP) equation. In both examples, only a single equation is involved. For the third one, we explore a coupled system with two equations, where the potential is a function of solutions to both equations. In all the examples, we transform the problem of function estimating to seeking some coefficients, which are further used to recover the potentials. In particular, in the first examples with single equations, we use the library-search algorithm to construct a low-dimensional approxiamtion of the true potential. More specifically, we use a library as follows

\[\Phi:=\{x,\sin(x),\cos(x),x^{2},\sin^{2}(x),\cos^{2}(x),e^{-x},e^{-2x},e^{-2x^{ 2}}\}.\]

In addition, we make use of compressed sensing to obtain a relatively sparse coefficient vector, with which only very few functions are involved in the final numerical solution. For the third one, owing to the special setting of the problem, one can easily relate the unknown function to some scalars hence the problem can well simplified. To present the performance of the proposed method, we define \(e_{\phi}\) and \(e_{V}\) as below:

\[e_{\psi}=\frac{1}{2}\frac{\|\psi^{\rm numerical}-\psi^{\rm exact}\|_{2}^{2}}{ \|\psi^{\rm exact}\|_{2}^{2}},\quad e_{V}=\frac{\|V^{\rm numerical}-V^{\rm exact }\|_{2}}{\|V^{\rm exact}\|_{2}}, \tag{18}\]

where \(\psi^{\rm numerical}\) and \(V^{\rm numerical}\) are the numerical results where \(\psi^{\rm exact}\) and \(V^{\rm exact}\) are corresponding references. In particular, \(e_{\psi}\) is from the \(l^{2}\) term in loss function \(J_{\rm e}\), which we show more interest in compared with the whole loss \(J_{\rm e}\).

**Example 1**.: One-soliton solution for NLSE.

In this example, we consider (10) with \(\beta=-1\), \(\gamma=1\), and \(\Omega=[-10,10]\). Moreover, the true potential has the exponential form as follows,

\[V(x)=4x^{2}-\exp(-2x^{2}).\]

Then the exact solution \(\psi(x,t)\) is has the following form:

\[\psi(x,t)=\exp(-x^{2}+2it).\]

We first show the convergence of \(e_{\psi}\) as the number of collocation points (\(M\)) and layers (\(N\)) in Figure 2. As for the left one, we fix the number of layers \(N\) to be \(200\) and increase the number of collocation points \(M\) in data. One can see that the by increasing \(M\) to \(20\), the error first shaply decrease to about \(10^{-9}\) and become stable. On the other hand, with \(M\) fixed as \(200\), one can observe a steady decay in \(e_{\psi}\) as using more layers, equivalently, smaller \(\Delta t\). Hence the performance of the NLS-Net is consistent with the Theorem 5. Since in our case, we deal with the inverse problem by iteratively solving the forward process, the effect of NLS-Net in approximating exact solution \(\psi\) serve as a basis in finding the true potential.

We further present the training performance in Figure 3. The left subfigure describes the dynamics of \(e_{\psi}\) in the training process while the right one plots \(e_{V}\). Here we set the regularization coefficient \(\lambda\) to be zero, which means \(J_{\rm e}=e_{\psi}\). Hence, we simply use gradient descent instead of proximal gradient descent. In particular, the initial learning rate is \(1.5\) and for each \(2000\) epochs, it decays by multiplying \(0.99\). Besides, we set a tolerence \(\tau\) to be \(10^{-10}\) as a creterion to further decrease the learning rate. Once \(e_{\psi}\) is below this tolerance, we multiple learning rate by \(0.4\) for each remaining step. For the network, it is composed of four layers and with input \(\tilde{\Psi}_{M}^{0}\) containing \(4000\) uniform collocation points in \([-10,10]\), i.e, \(M=4000\). Both the dynamics of \(e_{\psi}\) and \(e_{V}\) are plotted with logarithmic scale in y axis, hence one can observe a linear decay. We remark that a change in slope of lines near the end of training is resulted by exponential decreasing the learning rate since \(e_{\psi}\) already drops below a given small tolerance.

Finally we compare the numerical (Left) and exact (Right) potential in Figure 4. The domain \([-10,10]\) is discretized by \(4000\) uniform collocation points and we use \(6000\) epochs to training the potential. After that, one can observe that the apprximation of potential is almost impossible to be distinguished from the
Figure 3: Dynamics of loss function \(e_{\psi}\) (Left) and potential error \(e_{V}\) (Right) during training in **Example 1**. Here, the regularization coefficient \(\lambda\) is zero, i.e., \(J_{\textbf{e}}\) is the same as \(e_{\psi}\). Gradient descent is used with initial learning rate 1.5, which is multiplied by 0.99 every 2000 epochs. When \(e_{\psi}\) decreases to be lower than \(10^{-10}\), learning rate is multipied by 0.4 each step.

Figure 2: Convergence of \(e_{\psi}\) as the increase of length of input data (Left) and number of layers (Right). For the left one, we fix the number of layers to be 200 and for the right one, the number of neurons in each hidden layer is chosen to be 20.
goundtruth.

**Example 2**.: Gross-Pitaevskii (GP) equation.

In this example, we consider the GP equation, which is special case of (10) with periodic potential. In particular, we set \(\beta=\frac{1}{2},\gamma=-1\), the spatial domain \(\Omega=\)and \(V(x)=-\cos^{2}(x)\). For this equation, the exact solution [58] is

\[\psi(x,t)=\sin(x)\exp(-\frac{3ti}{2}).\]

First of all, we show the dynamics of \(e_{\psi}\) and \(e_{V}\) in the training process. To obtain a relatively good approximation, we choose \(\lambda=20\) in (14). As a comparison, we also present the case without regularization term in Figure 5. In this case, we simply use SGD. Even with a very small learning rate \(5*10^{-6}\), the loss \(e_{\psi}\) still oscillates, which may imply a relatively high level of nonconvexity in \(e_{\psi}\). Corresponding, the error \(e_{V}\) also presents a similar trend during training. Both of them shows that it is relatively difficult to attain the optimal solution with the only term \(e_{\psi}\) in \(J_{\text{e}}\). On the other hand, with a nonzero regulatization term, we can obtain a near exponential decay in both \(e_{\psi}\) and \(e_{V}\), which is shown in Figure 6. To obtain this result, we use proximal gradient descent. In particular, we use SGD to tackle the \(l^{2}\) term and a proximal operator is utilized to deal with the \(l^{1}\) norm, i.e. the second term in \(J_{\text{e}}\). For SGD, we begin with a learning rate \(5*10^{-3}\) and it decays by rate \(0.1\) every \(3000\) epochs. Similarly to the implementation in the last example, we shrink the learning by multiplying it with \(0.95\) each step after \(e_{\psi}\) has satisfied the condition that \(e_{\psi}<10^{-10}\). We remark that even though there are still regular oscillations in the dynamics because of SGD, both \(e_{\psi}\) and \(e_{V}\) display a steady decrease, which shows the effect of regularization term in \(J_{\text{e}}\).

In Figure 7, we compare the approximation of potential which is plotted on the left with the groundtruth on the other side. One can easily see that the numerical potential is approximately identical to the reference, which shows that the training process is effective.

**Example 3**.: A coupled system of NLSE.

In this example, we explore a coupled nonlinear Schrodinger equations as follows. For \(x\in[a,b]\), \(t\in[0,T]\),

\[\mathrm{i}\frac{\partial\psi_{j}}{\partial t}=-\frac{1}{2}\frac {\partial^{2}\psi_{j}}{\partial x^{2}}+V_{j}(x)\psi_{j}+i\alpha_{j}\frac{ \partial\psi_{j}}{\partial x}+f_{j}\left(\left|\psi_{1}\right|^{2},\left|\psi _{2}\right|^{2}\right)\psi_{j},\quad j=1,2,\] \[\psi_{j}(x,0)=\psi_{j}^{0}(x),\quad j=1,2\] \[\psi_{j}(a,t)=\psi_{j}(b,t)=0,\quad j=1,2,\]

where \(f_{1}(x,y)=x+\zeta_{1}y\), \(f_{2}(x,y)=\zeta_{2}x+y.\) Here, \(\zeta_{1}\) and \(\zeta_{2}\) are scalars. We can know that it admits one

Figure 4: A comparison of numerical (Left)and exact (Right) potential in Example 1 after training \(6000\) epochs. The spatial domain \([-10,10]\) is descritized by \(4000\) uniform collocation points and the time domain \([0,T]\) is partitioned by \(\Delta t=T/4\).
Figure 5: Dynamics of loss function \(e_{\psi}\) (Left) and potential error \(e_{V}\) (Right) during training in Example 2. Here, the regularization coefficient \(\lambda\) is zero, i.e., \(J_{\mathbf{c}}\) is the same as \(e_{\psi}\). Gradient descent is used with initial learning rate \(5*10^{-6}\), which is multiplied by \(0.9\) every \(3000\) epochs.

Figure 6: Dynamics of loss function \(e_{\psi}\) (Left) and potential error \(e_{V}\) (Right) during training in **Example 2**. Here, the regularization coefficient \(\lambda\) is \(20\). Proximal gradient descent is used with initial learning rate \(5*10^{-3}\), which is multiplied by \(0.9\) every \(3000\) epochs. When \(e_{\psi}\) decreases to be lower than \(10^{-10}\), learning rate is multiplied by \(0.95\) each step.
exact solution [58] as follows when \(\zeta_{1}=\zeta_{2}\), \(V_{1}(x)=V_{2}(x)=0\), and \(-\alpha_{1}=\alpha_{2}=0.5\). The overall loss is defined as \(J_{\mathbf{c}}=e_{\psi_{1}}+e_{\psi_{2}}\), where \(\mathbf{c}=[\zeta_{1},\zeta_{2}]\).

\[\begin{split}\psi_{1}^{\text{exact}}\ (x,t)&=\sqrt{\frac{2 \alpha}{1+\zeta}}\operatorname{sech}(\sqrt{2\alpha}(x-vt))\mathrm{e}^{ \mathrm{i}\left[\left(v-\delta\right)x-\left(\frac{2^{2}-\delta^{2}}{2}-\alpha \right)t\right]}\\ \psi_{2}^{\text{exact}}\ (x,t)&=\sqrt{\frac{2\alpha}{1+ \zeta}}\operatorname{sech}(\sqrt{2\alpha}(x-vt))\mathrm{e}^{\mathrm{i}\left[ \left(v+\delta\right)x-\left(\frac{2^{2}-\delta^{2}}{2}-\alpha\right)t\right] }\end{split} \tag{19}\]

In our numerical simulation, we use the following setting:

\[a=-20,\quad b=80,\quad v=0.,\quad\delta=0.5,\quad\alpha=1.0,\quad\zeta_{1}= \zeta_{2}=\frac{2}{3}.\]

The difference of this example is that the potential is a function of forward solutions \(\psi_{1}\) and \(\psi_{2}\), while in the first two examples, we have the prior information that the unknown potential is only a function of \(x\). However, as a compensation, the form of potential is specified and we only need to determine \(\zeta_{1}\) and \(\zeta_{2}\) rather than find a unknown function. We remark that even though the reference solutions defined in (19) correponds to the special case \(\zeta_{1}=\zeta_{2}\), it is not reasonable to assume this information hence we need to determine both \(\zeta_{1}\) and \(\zeta_{2}\). However, since there are still only two unknown scalars, we can easily plot the landscape of loss function \(J_{\mathbf{c}}\) as a guide to seek the optimal solution. In Figure 6, we show the contour of loss with respect to \(\zeta_{1}\) and \(\zeta_{2}\). First of all, it presents a reflection symmetry along the line \(\zeta_{1}=\zeta_{2}\), which can offer the insight that the optimal \(\zeta_{1}\) and \(\zeta_{2}\) are almost equal. Besides, the loss function shows the property of local convexity. Hence once the \([\zeta_{1},\zeta_{2}]\) arrives the basin of attraction, where we borrow the terminology in ecology, it is not demanding to attain the optimal solution. Both two observations can help us find the minimizer in a more efficient way.

Following the loss landscape, we further show the training performance. As before, we arrange the loss dynamics on the left and errors of inverse problems on the other side. Since in this situation, we only need to determine the unknowns scalars, we use the absolute error defined as

\[e_{\zeta_{j}}=\frac{|\zeta_{j}^{\text{numerical}}-\zeta_{j}^{\text{exact}}|}{| \zeta_{j}^{\text{exact}}|},\quad j=1,2.\]

Owing to the landscape analysis shown above, we can simply use gradient descent. Here we fix the learning rate to be \(100\). The reason for such a big learning rate may come from a small magnititude of the gradient with respect to the parameters \(\zeta_{1}\) and \(\zeta_{2}\). However, in spite of the big learning rate, we still obtain an exponeital decay of loss as well as the concerned errors. One can observe a convergence before about \(1000\) epochs.

Figure 7: A comparison of numerical (Left)and exact (Right) potential in Example 2 at final time \(T=1\). The spatial domain \([-10,10]\) is descritized by \(4000\) uniform collocation points and the time domain \([0,T]\) is partitioned by \(\Delta t=T\).
Figure 8: Dynamics of loss function \(e_{\psi}\) (Left), errors of \(\zeta_{1}\) and \(\zeta_{2}\) (Right) during training in Example 3. Gradient descent is used with fixed learning rate 100.
Conclusion remarks

In this work, we solve some inverse problems of nonlinear Schrodinger equations as a learning process of a special type of convolutional neural network (CNN). The first novelty of this work comes from the special architecture of the neural network, which is designed based on the equation and the operator splitting method. In particular, the convolution layers serve as a linear operator while the activation functions constitude the nonlinear part. A composition of these two forms the time-marching process. Owing to this design, layers and activation functions have physical interpretations. An intermediate convenience is that one can easily choose the type of activation functions from the information of equations. Besides, the existence of an optimal network with this special structure can be guaranteed since the neural network can be directly derived from the equation. We remark that the existence is not a direct result of the famous universal approximation theorem for neural networks since the proposed network has a special architecture while the theorem does not specify the structure of the optimal network. Moreover, based on the existing theoretical result of operator splitting, one can easily guarantee the convergence of the neural network solution with respect to the depth and width. Overall, the particular structure derived from the equation enhances the explainability of network. Instead of using neural network as a black box, the inherited physics sheds light on the approximation effect of the neural network.

For the inverse problem, we are interested in approximating the potential function in the Schrodinger equations. Different from the coefficient inverse problem, function approximation is more complicated. For example, one needs a neural network to approximate a function while for coefficients, some neurons are sufficient. To this end, we use a library-search method to transform the original problem to a coefficient inverse problem, which is the second contribution. In particular, we construct a library that contains a set of basis functions then we can use a linear combination of the bases to express the unknown potential. In other words, we project the solution of inverse problems to an approximation space so we only need to determine the coefficients with respect to the basis functions. In this way, the complexity of neural network can be effectively reduced. However, to approximate a function space, we need to incorporate sufficiently many bases to obtain a relatively good accuracy. On the other hand, the form of potential is not so complicated, which may be a linear combination of a small number of different functions. Hence, the coefficient vector correponding to basis functions should have some level of sparsity, which motivates the need of using compressed sensing.
