The Deep Ritz method: A deep learning-based numerical algorithm for solving variational problems

Introduction

Deep learning has had great success in computer vision and other artificial intelligence tasks . Underlying this success is a new way to approximate functions, from an additive construction commonly used in approximation theory to a compositional construction used in deep neural networks. The compositional construction seems to be particularly powerful in high dimensions. This suggests that deep neural network based models can be of use in other contexts that involve constructing functions. This includes solving partial differential equations, molecular modeling, model reduction, etc. These aspects have been explored recently in .

In this paper, we continue this line of work and propose a new algorithm for solving variational problems. We call this new algorithm the Deep Ritz method since it is based on using the neural network representation of functions in the context of the Ritz method. The Deep Ritz method has a number of interesting and promising features, which we explore later in the paper.

The Deep Ritz Method

An explicit example of the kind of variational problems we are interested in is

and $H$ is the set of admissible functions (also called trial function, here represented by $u$), $f$ is a given function, representing external forcing to the system under consideration. Problems of this type are fairly common in physical sciences. The Deep Ritz method is based on the following set of ideas:

Deep neural network based approximation of the trial function.

A numerical quadrature rule for the functional.

An algorithm for solving the final optimization problem.

Our experience has suggested that the smoothness of the activation function $\phi$ plays a key role in the accuracy of the algorithm. To balance simplicity and accuracy, we have decided to use

The last term in (3), the residual connection, makes the network much easier to train since it helps to avoid the vanishing gradient problem . The structure of the two blocks, including two residual connections, is shown in Figure 1.

The full $n$-layer network can now be expressed as:

Here in the left-hand side and in what follows, we will use $\theta$ to denote the full parameter set $\{\theta,a,b\}$. Substituting this into the form of $I$, we obtain a function of $\theta$, which we should minimize.

For the functional that occurs in (2), denote:

then we are left with the optimization problem:

2 The stochastic gradient descent algorithm and the quadrature rule

To finish describing the algorithm, we need to furnish the remaining two components: the optimization algorithm and the discretization of the integral in $I$ in (2) or $L$ in (8). The latter is necessary since computing the integral in $I$ (or $L$) explicitly for functions of the form (6) is quite an impossible task.

In machine learning, the optimization problem that one encounters often takes the form:

where each term at the right-hand side corresponds to one data point. $n$, the number of data points, is typically very large. For this problem, the algorithm of choice is the stochastic gradient descent (SGD) method, which can be described as follows:

Here $\{\gamma^{k}\}$ are i.i.d random variables uniformly distributed over $\{1,2,\cdots,n\}$. This is the stochastic version of the gradient descent algorithm (GD). The key idea is that instead of computing the sum when evaluating the gradient of $L$, we simply randomly choose one term in the sum. Compared with GD, SGD requires only one function evaluation of $n$ function evaluations at each iteration. In practice, instead of picking one term, one chooses a ”mini-batch” of terms at each step.

At a first sight, our problem seems different from the ones that occur in machine learning since there are no data involved. The connection becomes clear once we view the integral in $I$ as a continuous sum, each point in $\Omega$ then becomes a data point. Therefore, at each step of the SGD iteration, one chooses a mini-batch of points to discretize the integral. These points are chosen randomly and the same quadrature weight is used at every point.

Note that if we use standard quadrature rules to discretize the integral, then we are bound to choose a fixed set of nodes. In this case, we run into the risk where the integrand is minimized on these fixed nodes but the functional itself is far from being minimized. It is nice that SGD fits naturally with the needed numerical integration in this context.

In summary, the SGD in this context is given by:

where for each $k$, $\{x_{j,k}\}$ is a set of points in $\Omega$ that are randomly sampled with uniform distribution. To accelerate the training of the neural network, we use the Adam optimizer version of the SGD .

Numerical Results

where $\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times\{0\}$. The solution to this problem suffers from the well-known ”corner singularity” caused by the nature of the domain . A simple asymptotic analysis shows that at the origin, the solution behaves as $u(x)=u(r,\theta)\sim r^{\frac{1}{2}}\sin\frac{\theta}{2}$ . Models of this type have been extensively used to help developing and testing adaptive finite element methods.

The network we used to solve this problem is a stack of four blocks (eight fully-connected layers) and an output layer with $m=10$. There are a total of $811$ parameters in the model. As far as we can tell, this network structure is not special in any way. It is simply the one that we used.

The boundary condition causes some problems. Here for simplicity, we use a penalty method and consider the modified functional

We choose $\beta=500$. The results from the Deep Ritz method is shown in see Figure 2(a). For comparison, we also plot the result of the finite difference method with $\Delta x_{1}=\Delta x_{2}=0.1$ ($1,681$ degrees of freedom), see Figure 2(b).

To analyze the error more quantitatively, we consider the following problem

where $\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times\{0\}$. This problem has an explicit solution $u^{*}(x)=r^{\frac{1}{2}}\sin\frac{\theta}{2}$ in polar coordinates. The error $e=\max|u^{*}(x)-u_{h}(x)|$, where $u^{*}$ and $u_{h}$ are the exact and approximate solutions respectively, is shown in Table 1 for both the Deep Ritz method and the finite difference method (on uniform grids). We can see that with fewer parameters, the Deep Ritz method gives more accurate solution than the finite difference method.

Being a naturally nonlinear variational method, the Deep Ritz method is also naturally adaptive. We believe that this contributes to the better accuracy of the Deep Ritz method.

2 Poisson equation in high dimension

Experiences in computer vision and other artificial intelligence tasks suggest that deep learning-based methods are particularly powerful in high dimensions. This has been confirmed by the results of the Deep BSDE method . In this subsection, we investigate the performance of the Deep Ritz method in relatively high dimension.

The solution of this problem is simply $u(x)=\sum_{k=1}^{5}x_{2k-1}x_{2k}$, and we will use the exact solution to compute the error of our model later.

For the network structure, we stack six fully-connected layers with three skip connections and a final linear layer, and there are a total of 671 parameters. For numerical integration, at each step of the SGD iteration, we sample 1,000 points in $\Omega$ and 100 points at each hyperplane that composes $\partial\Omega$. We set $\beta=10^{3}$. After 50,000 iterations, the relative $L_{2}$ error was reduced to about $0.4\%$. The training process is shown in Figure 3(a).

Also shown in Figure 3(b) is the training process for the problem:

with $d=100$ with a similar network structure (stack 3 blocks of size m=100). The solution of this problem is $u(x)=\sum_{k}x_{k}^{2}$. After 50000 iterations, the relative error is reduced to about $2.2\%$.

3 An example with the Neumann boundary condition

The exact solution is $u(x)=\sum_{k}\cos(\pi x_{k})$

without any penalty function for the boundary.

With a similar network structure the relative $L_{2}$ error reaches $1.3\%$ for $d=5$ and $1.9\%$ for $d=10$. The training process is shown in Figure 4.

4 Transfer learning

An important component of the training process is the initialization. Here we investigate the benefit of transferring weights in the network when the forcing function $f$ is changed.

where $\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times\{0\}$. Here we used a mixture of rectangular and polar coordinates. The exact solution is

The network consists of a stack of 3 blocks with m=10, that is, six fully-connected layers and three residual connections and a final linear transformation layer to obtain $u$. We show how the error and the weights in the layers change during the training period in Figure 5.

We also transfer the weights from the problem:

where $\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times\{0\}$.

The error and the weights during the training period are also shown in Figure 5. We see that transferring weights speeds up the training process considerably during the initial stage of the training. This suggests that transferring weights is a particularly effective procedure if the accuracy requirement is not very strigent.

5 Eigenvalue problems

Problems of this kind occur often in quantum mechanics where $v$ is the potential function.

There is a well-known variational principle for the smallest eigenvalue:

The functional we minimize here is called the Rayleigh quotient.

To avoid getting the trivial optimizer $u=0$, instead of using the functional

One might suggest that with the last penalty term, the denominator in the Rayleigh quotient is no longer necessary. It turns out that we found in practice that this term still helps in two ways: (1) In the presence of this denominator, there is no need to choose a large value of $\gamma$. For the harmonic oscillator in $d=5$, we choose $\beta=2000$, $\gamma$ to be 100 and this seems to be large enough. (2) This term helps to speed up the training process.

To solve this problem, we build a deep neural network much like the Densenet . There are skip connections between every pairwise layers, which help gradients flow through the whole network. The network structure is shown in Figure 6.

We use an activation function $\phi(x)=max(0,x)^{2}$. If we use the same activation function as before, we found that the gradients can become quite large and we may face the gradient explosion problem.

The remaining components of the algorithm are very much the same as before.

The problem is then equivalent to solving:

The smallest eigenvalue is $\lambda_{0}=d\pi^{2}$.

The results of the Deep Ritz method in different dimensions are shown in Table 2.

The results in different dimensions are shown in Table 3.

The results deteriorate substantially as the dimension is increased. We believe that there is still a lot of room for improving the results. We will leave this to future work.

Discussion

We proposed a variational method based on representing the trial functions by deep neural networks. Our limited experience with this method suggests that it has the following advantages:

It is less sensitive to the dimensionality of the problem and has the potential to work in rather high dimensions.

The method is reasonably simple and fits well with the stochastic gradient descent framework commonly used in deep learning.

We also see a number of disadvantages that need to be addressed in future work:

The variational problem that we obtain at the end is not convex even when the initial problem is. The issue of local minima and saddle points is non-trivial.

At the present time, there is no consistent conclusion about the convergence rate.

The treatment of the essential boundary condition is not as simple as for the traditional methods.

In addition, there are still interesting issues regarding the choice of the network structure, the activation function and the minimization algorithm. The present paper is far from being the last word on the subject.

Acknowledgement: We are grateful to Professor Ruo Li and Dr. Zhanxing Zhu for very helpful discussions. The work of E and Yu is supported in part by the National Key Basic Research Program of China 2015CB856000, Major Program of NNSFC under grant 91130005, DOE grant DE-SC0009248, and ONR grant N00014-13-1-0338.

References