# Adaptive importance sampling for Deep Ritz

Xiaoliang Wan
&Tao Zhou
&Yuancheng Zhou
Department of Mathematics and Center for Computation and Technology, Louisiana
State University, Baton Rouge 70803, USA. Email<EMAIL_ADDRESS>LSEC,
Institute of Computational Mathematics and Scientific/Engineering Computing,
AMSS, Chinese Academy of Sciences, Beijing, China. Email<EMAIL_ADDRESS>LSEC, Institute of Computational Mathematics and Scientific/Engineering
Computing, AMSS, Chinese Academy of Sciences, Beijing, China. Email:
<EMAIL_ADDRESS>
###### Abstract

We introduce an adaptive sampling method for the Deep Ritz method aimed at
solving partial differential equations(PDEs). Two deep neural networks are
used. One network is employed to approximate the solution of PDEs, while the
other one is a deep generative model used to generate new collocation points
to refine the training set. The adaptive sampling procedure consists of two
main steps. The first step is solving the PDEs using the Deep Ritz method by
minimizing an associated variational loss discretized by the collocation
points in the training set. The second step involves generating a new training
set, which is then used in subsequent computations to further improve the
accuracy of the current approximate solution. We treat the integrand in the
variational loss as an unnormalized probability density function(PDF) and
approximate it using a deep generative model called bounded KRnet. The new
samples and their associated PDF values are obtained from the bounded KRnet.
With these new samples and their associated PDF values, the variational loss
can be approximated more accurately by importance sampling. Compared to the
original Deep Ritz method, the proposed adaptive method improves accuracy,
especially for problems characterized by low regularity and high
dimensionality. We demonstrate the effectiveness of our new method through a
series of numerical experiments.

_Keywords_ Deep Ritz method $\cdot$ generative model $\cdot$ importance
sampling

##  1 Introduction

Recent years have witnessed fast development in the field of machine learning.
There has been a growing interest in utilizing machine learning, particularly
deep learning, to address scientific computing problems. Solving partial
differential equations (PDEs) is one of the most important tasks in scientific
computing. Consequently, many studies on deep learning have been recently
conducted to solve problems associated with PDEs. Classical numerical methods
such as finite difference method and finite element method face the challenge
known as the “curse of dimensionality” when dealing with high-dimensional
PDEs. In contrast to classical numerical methods, deep learning methods
utilize Monte Carlo approximation instead of computational mesh for
discretization, which makes them more suitable for solving high-dimensional
PDEs. Available deep learning methods can be categorized based on various
formulations of PDEs, such as PINNs [8] and DGM [10] for the strong form of
PDEs, WANs [15] for the weak form of PDEs and Deep Ritz [14] for the
variational form of PDEs. In this work, our primary focus will be on the Deep
Ritz method.

As deep learning methods are meshless approaches, a critical concern is how to
select training points to effectively approximate the integral for the loss
functional. The accuracy of the discrete loss functional is pivotal to
achieving an accurate approximation of the PDE solution. In the original
methods of PINNs, DGM, WANs, and Deep Ritz, the loss functionals are
approximated using the Monte Carlo method with training points randomly
generated from a uniform distribution. However, the uniform sampling strategy
is often inefficient if the PDE solution is of low regularity and the
situation deteriorates as dimensionality increases. In order to alleviate this
problem, some adaptive sampling strategies have been developed to improve the
Monte Carlo approximation of the loss functional. However, it’s important to
note that most adaptive sampling methods are primarily designed for PINNs/DGM,
where a change of measure for the integration of the residual does not affect
the well-posedness of the formulation. These adaptive sampling strategies
include RAR [5], RAD [13], RAR-G [13], DAS-G [12], FI-PINNs [3]. All of these
adaptive sampling strategies provide a certain way to update the current
training set by adding some important new training points, which results in a
weighted integration of the residual, i.e., the Lebesgue measure is replaced
by a new one. In the Deep Ritz method, the residual term is substituted with a
variational term whose minimum value is in general not zero. A straightforward
replacement of the Lebesgue measure will result in an incorrect variational
term such that minimizing the modified loss functional leads to an erroneous
solution. As a consequence, the adaptive sampling strategies RAR, RAD, RAR-G,
DAS-G, and FI-PINNs cannot be directly applied to the Deep Ritz method. In the
Deep Ritz method, it is essential to modify the measure without changing the
integral, in other words, importance sampling is a necessity. Importance
sampling was also considered in DAS-R [12] for PINNs.

In this work, we introduce a deep adaptive sampling method for the Deep Ritz
method based on importance sampling. Assuming no prior knowledge is available,
we intend to use the approximation of the optimal measure for the change of
measure required by importance sampling. To approximate the optimal measure
induced by the absolute value of the integrand, we need a generic density
model. In this work we employ a recently developed deep generative model
called bounded KRnet [16] to take advantage of its modeling capability for
high-dimensional distributions. In our adaptive sampling methods, we utilize
two deep neural networks: one for approximating the PDE solution and the other
one for refining the training set. The approximate PDE solution is attained by
minimizing the discretized loss functional, which comprises a variational term
and a boundary term. The bounded KRnet defines an invertible transport map
from the target distribution on a bounded rectangular domain to a prior
(uniform) distribution on a hypercube $[-1,1]^{d}$, which is achieved by
minimizing the Kullback-Leibler (KL) divergence between the probability
density function (PDF) induced by bounded KRnet and the PDF given by the
optimal measure. We have proposed three adaptive schemes. Our basic adaptive
scheme is to completely replace the previous training set with new training
points generated from the learned bounded KRnet. We also devise two
alternative PDF models, which correspond to two different schemes to update
the previous training set, to enhance the performance of the basic adaptive
scheme. In all provided adaptive schemes, the updated training set is
subsequently employed to further enhance the accuracy of the Deep Ritz method.
In summary, the approximate PDE solution and the bounded KRnet will be updated
alternately until a certain criterion is satisfied.

The rest of the paper is organized as follows. In the next section, we briefly
introduce the Deep Ritz method for the approximation of PDEs. The adaptive
sampling method for Deep Ritz is presented in detail in section 3. Then, in
section 4, we demonstrate the performance of our method with numerical
experiments. We conclude our work in section 5.

##  2 Deep Ritz method for PDEs

In the Deep Ritz method, our objective is to minimize the following
variational problem

$\min_{u\in V}I(u)$

where $I(u)$, as a characteristic functional of an elliptic problem defined on
a compact domain $\Omega\subset\mathbb{R}^{d}$ with $d\in\mathbb{N}_{+}$,
takes the form

$I(u)=\int_{\Omega}G(u(x))\,dx,$

with $G(u)$ being a function of $u(x)$. To deal with the boundary conditions
$B(x,u(x))=0$, $x\in\partial\Omega$, we include a penalty term into the
objective functional and consider the following optimization problem

$\min_{u\in V}J(u)$ (1)

with

$J(u)=\int_{\Omega}G(u(x))\,dx+\beta\|B(x,u(x))\|_{\partial\Omega,2}^{2},$ (2)

where $\beta>0$ is a penalty parameter and $\|\cdot\|_{\partial\Omega,2}$
indicates the $L_{2}$ norm on $\partial\Omega$.

The Deep Ritz method introduces a neural network $u(\cdot,\theta)$ to
approximate the solution $u$ and uses the Monte Carlo method to discretize the
integrals in $J(u)$. Then solving the original variational problem becomes
minimizing a loss function over the neural network parameters. Let
$S_{\Omega}=\left\\{x_{i}\right\\}_{i=1}^{N_{v}}$ and
$S_{\partial\Omega}=\left\\{x_{j}\right\\}_{i=1}^{N_{b}}$ be two sets of
uniformly distributed collocation points on $\Omega$ and $\partial\Omega$
respectively. The discretized optimization problem is

$\min_{\theta}J_{N}(u(x,\theta))$

with

$J_{N}(u(x,\theta))=\frac{1}{N_{v}}\sum_{i=1}^{N_{v}}G(u(x_{i},\theta))+\frac{\beta}{N_{b}}\sum_{j=1}^{N_{b}}B(x_{j},u(x_{j},\theta))^{2},$

which is then solved by stochastic gradient-based methods. Compared to
traditional methods, the Deep Ritz method is reasonably simple and less
sensitive to the dimensionality of the problem [14].

Recently, many adaptive sampling methods have been established to improve the
performance of PINNs, most of which are residual-based methods since the
objective of PINNs is to minimize the residual on selected collocation points.
The Deep Ritz method is not involved with the residual, implying that the
residual-based adaptive sampling methods designed for PINNs cannot be directly
applied to the Deep Ritz method. In this work, we develop an adaptive sampling
method tailored for the Deep Ritz method, with a particular focus on reducing
the statistical error in the approximate solution.

##  3 Adaptive sampling for Deep Ritz method

The accuracy of the Deep Ritz method depends on various factors involved in
the computation. We are interested in the error induced by discretization,
which contains two parts: statistical error and approximation error. The
approximation error comes from the difference between the approximation space
and the true space. The statistical error arises from discretizing the loss
functional through Monte Carlo approximation. The statistical error depends on
the random samples used by the Monte Carlo method while the approximation
error relies on the modeling capability of the neural network for the
approximation of the PDE solution [12]. We aim to reduce the statistical error
in the Deep Ritz method by seeking a more accurate Monte Carlo approximation
of the loss functional, which will be achieved by variance induction through
adaptive sampling.

###  3.1 Importance sampling for Deep Ritz method

For an integrand $G(u)$ with low regularity, the integral can be dominated by
the values of $G(u)$ in a small region within $\Omega$. In such cases,
uniformly increasing the number of training points across the entire domain
$\Omega$ may be not effective in reducing the discretization error of the loss
functional. To illustrate this point, we provide the following example.
Suppose $G(u(x))=1/\sqrt{2\pi}\exp(-x^{2}/2\sigma^{2})$ on $\Omega=[-1,1]$ and
$\sigma\ll 1<1/3$. As the parameter $\sigma$ decreases, the shape of $G(u(x))$
becomes sharper around zero. Using $3\sigma$ principle, the integral can be
computed as

$\int_{\Omega}G(u(x))\,dx=\int_{\Omega}\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^{2}}{2\sigma^{2}})\,dx\approx\sigma.$

The Monte Carlo approximation of the integral is

$\hat{I}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\sqrt{2\pi}}\exp(-\frac{x_{i}^{2}}{2\sigma^{2}})$

where $\left\\{x_{i}\right\\}$, $i=1,\cdots,N$ is a set of uniformly samples
in $\Omega$. The relative error of the approximation is

$\frac{(\operatorname{Var}{I})^{\frac{1}{2}}}{\sigma}\approx\frac{1}{\sigma}\sqrt{\frac{\sqrt{2}\Gamma(\frac{5}{4})\sigma}{2\pi
N}}=C(\sigma N)^{-1/2}$

where $\Gamma(\cdot)$ is the Gamma function and $C$ is a constant. Then, in
order to obtain a relative error of $O(1)$, we need around $O(1/\sigma)$
number of samples. This means the uniform samples are not effective for this
example. In order to improve the accuracy, more effective samples generated by
variance reduction techniques are required.

One of the most commonly used variance reduction techniques is importance
sampling. The idea of importance sampling is to generate a training set from a
well-chosen distribution and get an unbiased estimator,

$\displaystyle I(u)=\int_{\Omega}G(u(x,\theta))\,dx$
$\displaystyle=\int_{\Omega}\frac{G(u(x,\theta))}{p(x)}p(x)\,dx$ (3)
$\displaystyle=\operatorname{\mathbb{E}}_{p}\left[\frac{G(u(X,\theta))}{p(X)}\right]\approx\frac{1}{N_{v}}\sum_{i=1}^{N_{v}}\frac{G(u(x_{i},\theta))}{p(x_{i})}$

where $\left\\{x_{i}\right\\}$, $i=1,\cdots,N_{v}$ is a set of samples from
the probability density function $p(x)$. If the variance of
$G(u(X,\theta)/p(X))$ with respect to density function of $p$

$\displaystyle\sigma_{p}=\operatorname{Var}_{p}\left[\frac{G(u(X,\theta))}{p(X)}\right]$
$\displaystyle=\operatorname{\mathbb{E}}_{p}\left[\frac{G^{2}(u(X,\theta))}{p^{2}(X)}\right]-\left(\operatorname{\mathbb{E}}_{p}\left[\frac{G(u(X,\theta))}{p(X)}\right]\right)^{2}$
(4)
$\displaystyle=\int_{\Omega}\frac{G^{2}(u(x,\theta))}{p(x)}\,dx-\left(\int_{\Omega}G(u(x,\theta))\,dx\right)^{2}$

is smaller than that of $G(u(x,\theta))$ with respect to the uniform
distribution on $\Omega$, the performance of the Monte Carlo method will be
improved.

From importance sampling (3), if we assume $G(u(x,\theta))\geq 0$ for all
$x\in\Omega$, one can find the best choice of the density function $p(x)$ is

$p^{\star}(x)=\frac{G(u(x,\theta))}{\mu}$

where $\mu=\int_{\Omega}G(u(x,\theta))\,\mathrm{dx}$. We can check this fact
by plugging the $p^{\star}$ into equation (4). The variance is zero which
means the error is also zero. In practice, $\mu$ is what we need to compute
and therefore we do not know $p^{\star}$, however, results suggest variance
reduction can be achieved if we can find a $p$ close to $p^{\star}$ [12].

In the Deep Ritz method, one may encounter the situation when the integrand
$G(u(x,\theta))$ is not always non-negative for $x\in\Omega$. For this case,
the best choice for the selection of $p(x)$ is

$p^{\star}(x)=\frac{|G(u(x,\theta))|}{\mu},$ (5)

where $\mu=\int_{\Omega}|G(u(x,\theta))|\,dx$. More specifically, we have for
any density function $p(x)$

$\displaystyle\operatorname{\mathbb{E}}_{p^{\star}}\left[\frac{G^{2}(u(X,\theta))}{{p^{\star}}^{2}(X)}\right]$
$\displaystyle=\int_{\Omega}\frac{G^{2}(u(x,\theta))}{p^{\star}(x)}\,dx$ (6)
$\displaystyle=\left(\int_{\Omega}|G(u(x,\theta))|\,dx\right)^{2}$
$\displaystyle=\left(\int_{\Omega}\frac{|G(u(x,\theta))|}{p(x)}p(x)\,dx\right)^{2}$
$\displaystyle\leq\int_{\Omega}\frac{|G(u(x,\theta))|^{2}}{p^{2}(x)}p(x)\,dx\cdot\int_{\Omega}p(x)\,dx$
$\displaystyle=\operatorname{\mathbb{E}}_{p}\left[\frac{G^{2}(u(X,\theta))}{p^{2}(X)}\right],$

where the Cauchy-Schwarz inequality is applied. We add
$-(\int_{\Omega}G(u(x,\theta))\,\mathrm{dx})^{2}$ to both side of inequality
(6) and obtain

$\sigma_{p^{\star}}=\operatorname{Var}_{p^{\star}}\left[\frac{G(u(X,\theta))}{p^{\star}(X)}\right]\leq\operatorname{Var}_{p}\left[\frac{G(u(X,\theta))}{p(X)}\right]=\sigma_{p}.$

Also, we know from the proof

$\sigma_{p^{\star}}=\left(\int_{\Omega}|G(u(x,\theta))|\,dx\right)^{2}-\left(\int_{\Omega}G(u(x,\theta))\,dx\right)^{2}>0,$

if the integrand $G(u(x,\theta))$ contains both positive part and negative
part. Then the error of the Monte Carlo approximation satisfies

$\operatorname{Var}_{p}^{\frac{1}{2}}\left[\frac{1}{N_{v}}\sum_{i=1}^{N_{v}}\frac{G(u(X_{i},\theta))}{p(X_{i})}\right]=\frac{\sigma_{p}}{\sqrt{N_{v}}}\geq\frac{\sigma_{p^{\star}}}{\sqrt{N_{v}}}>0.$

Unlike the case when the integrand $G(u(x,\theta))$ is non-negative, the
variance will not reduce to zero but to a positive number. If we want to
further reduce the variance after $p$ reaches to $p^{\star}$, we have to
increase the number of samples.

###  3.2 Bounded KRnet: PDF approximation and sample generation

From the previous section, we know that $p^{\star}=|G(u(x,\theta))|/\mu$ for a
fixed parameter $\theta$ is theoretically the optimal density function for
implementing importance sampling in the Deep Ritz method. In order to
calculate the estimator in $(\ref{importance sampling})$, we may choose $p$
that is sufficiently close to the $p^{\star}$. However, PDF approximation is a
challenging task, especially in high-dimensional spaces. Since the shape of
the graph $|G(u(x,\theta))|/\mu$ can be complex in high-dimensional cases, it
is hard to use the exponential family of distributions and mixture models to
approximate this function. Thus, deep generative models which have stronger
representation ability are taken into consideration. In particular, we use a
recently developed model, called bounded KRnet [16], to approximate
$p^{\star}$ and generate samples as well. Bounded KRnet is a kind of
normalizing flow [7]. It is an invertible transport map between a prior
distribution and the target distribution. The bounded KRnet is defined on a
compact domain while the original KRnet is defined on the whole space. This
difference makes bounded KRnet a more suitable tool than the original KRnet to
approximate $|G(u(x,\theta))|/\mu$ since in most cases, the computation domain
of a PDE is compact.

Suppose $\Omega\subset\mathbb{R}^{d}$ is a compact set in $\mathbb{R}^{d}$.
Let $X\in\Omega$ be a random vector associated with the target PDF $p_{X}(x)$.
Let $Z\in\Omega$ be another random vector with a known PDF $p_{Z}(z)$, e.g.,
PDF of the uniform distribution. Our goal is to find an invertible transport
map $f:\Omega\rightarrow\Omega$ from $X$ to $Z$. Once we obtain the map $f$,
we can calculate the PDF of $X=f^{-1}(Z)$ using the change of variables

$p_{X}(x)=p_{Z}(z)|\det\nabla_{x}f|.$ (7)

This is an explicit density model for $X$. Using the inverse map $f^{-1}$, we
can also easily obtain exact samples of $X$ by $X=f^{-1}(Z)$. The idea of
KRnet originates from Knothe-Rosenblatt rearrangement [9, 1]. It integrates
the K-R rearrangement into the definition of the invertible transport map $f$.
In practice, the invertible map $f$ is composed of a sequence of simple
invertible maps and each map in the sequence is a shallow neural network. In
bounded KRnet, these simple invertible map is constructed by a non-linear map
from $[-1,1]^{d}$ to $[-1,1]^{d}$ while in original KRnet, the map is from
$\mathbb{R}^{d}$ to $\mathbb{R}^{d}$.

Suppose $f_{\text{bKRnet}}(\cdot,\theta_{f})$ is the invertible transport map
induced by the bounded KRnet where $\theta_{f}$ includes the model parameters.
According to equation (7), we can explicitly write down the PDF with respect
to $f_{\text{bKRnet}}$

$p_{\text{bKRnet}}(x,\theta_{f})=p_{Z}(f_{\text{bKRnet}}(x,\theta_{f}))|\det\nabla_{x}f_{\text{bKRnet}}|$

The sample of $p_{\text{bKRnet}}(x,\theta_{f})$ can be generated through
$X=f_{\text{bKRnet}}^{-1}(Z)$ by sampling $Z$. In bounded KRnet, the random
vector $Z$ is typically assumed to be a uniform random variable on the cube
$[-1,1]^{d}$. Our aim is to find out the optimal parameters
$\theta_{f}^{\star}$ such that the PDF
$p_{\text{bKRnet}}(x,\theta_{f}^{\star})$ is close enough to a target PDF $p$.
In order to find out $\theta_{f}^{\star}$, we minimize the KL divergence from
$p$ to $p_{\text{BKRnet}}(x,\theta_{f}^{\star})$, i.e.,

$\theta_{f}^{\star}=\mathop{\arg\min}_{\theta_{f}}D_{\text{KL}}(p\,||\,p_{\text{bKRnet}},(x,\theta_{f})).$
(8)

where

$\displaystyle D_{\text{KL}}(p\,||\,p_{\text{bKRnet}}(x,\theta_{f}))$ (9)
$\displaystyle=$
$\displaystyle\int_{\Omega}p(x)\log(p(x))\,{dx}-\int_{\Omega}p(x)\log(p_{\text{bKRnet}}(x,\theta_{f})))\,{dx}.$

It can be seen the first term in equation (9) is irrelevant to the density
model and therefore can be ignored. Hence the minimization problem (8) is
equivalent to minimizing the second term in equation (9) which is the cross
entropy $H(p,p_{\text{bKRnet}})$ from $p$ to
$p_{\text{bKRnet}}(x,\theta_{f}^{\star})$. That is

$\displaystyle\theta_{f}^{\star}=$
$\displaystyle\mathop{\arg\min}_{\theta_{f}}H(p,p_{\text{bKRnet}}).$ (10)

For the Deep Ritz method, we set $p=p^{\star}$ in equation (5). We also apply
importance sampling to the computation of the cross entropy in problem (10)
using a prior PDF model $\tilde{p}$, i.e.

$H(p^{\star},p_{\text{bKRnet}})\approx-\frac{1}{N_{v}}\sum_{i=1}^{N_{v}}\frac{|G(u(x_{i},\theta))|\log(p_{\text{bKRnet}}(x_{i},\theta_{f}))}{\mu\tilde{p}(x_{i})},$
(11)

where $\left\\{x_{i}\right\\}$, $i=1,\cdots,N_{v}$ are samples generated from
the prior PDF $\tilde{p}$. The prior PDF $\tilde{p}$ will be specified in our
adaptive sampling procedure presented in the next section.

###  3.3 Adaptive sampling method

We now present our adaptive sampling algorithms. We first present a basic
algorithm, where all collocation points in the training set will be replaced
at each adaptivity iteration in terms of a new density model. Then we propose
two less aggressive but more robust strategies, where the training set is
partially updated at each adaptivity iteration in terms of a mixture density
model.

####  3.3.1 Basic adaptive algorithm(replace all training points)

Let $S_{\Omega}^{k}=\left\\{x_{\Omega,i}^{k},i=1,\cdots,N_{v}\right\\}$ and
$S_{\partial\Omega}^{k}=\left\\{x_{\partial\Omega,j}^{k},j=1,\cdots,N_{b}\right\\}$
be two sets of collocation points at the $k$th adaptivity iteration for the
Deep Ritz method, which respectively deal with the variational term and the
penalty term for the boundary conditions. Let
$S_{B}^{k}=\left\\{x_{B,l}^{k},l=1,\cdots,N_{B}\right\\}$ be the training
points for bounded KRnet at the $k$th adaptivity iteration. In the $(k+1)$th
adaptivity iteration, we first solve the following optimization problem from
the Deep Ritz method

$\theta^{k+1}=\mathop{\arg\min}_{\theta}\frac{1}{N_{v}}\sum_{i=1}^{N_{v}}\frac{G(u(x_{\Omega,i}^{k},\theta))}{p_{\text{bKRnet}}(x_{\Omega,i}^{k},\theta_{f}^{k})}+\frac{\beta}{N_{b}}\sum_{j=1}^{N_{b}}B^{2}(x_{\partial\Omega,j}^{k},u(x_{\partial\Omega,j}^{k},\theta))$
(12)

where $p_{\text{bKRnet}}(\,\cdot\,,\theta_{f}^{k})$ is the PDF given by
bounded KRnet when $k>0$, and chosen as a unfirom distribution when $k=0$.
Correspondingly, $S_{\Omega}^{k}$ is sampled from a uniform distribution when
$k=0$ while it is sampled from $p_{\text{bKRnet}}(\,\cdot\,,\theta_{f}^{k})$
when $k>0$. Once we have computed $\theta^{k+1}$, we can update the bounded
KRnet by

$\theta_{f}^{k+1}=\mathop{\arg\min}_{\theta_{f}}-\frac{1}{N_{B}}\sum_{l=0}^{N_{B}}\frac{|G(u(x_{B,l}^{k},\theta^{k+1}))|\log(p_{\text{bKRnet}}(x_{B,l}^{k},\theta_{f}))}{\mu^{k+1}p_{\text{bKRnet}}(x_{B,l}^{k},\theta_{f}^{k})}$
(13)

where $S_{B}^{k}$ is sampled from a uniform distribution when $k=0$ while it
is sampled from $p_{\text{bKRnet}}(\,\cdot\,,\theta_{f}^{k})$ when $k>0$. This
way, the training set $S_{B}^{k}$ can be the same as the training set
$S_{\Omega}^{k}$ for the Deep Ritz method. $\mu^{k+1}$ in equation (13) can be
omitted since it is a constant. After $\theta_{f}^{k+1}$ is computed, we
obtain the bounded KRnet at the $(k+1)$th adaptivity iteration and then use it
to sample new training sets $S_{\Omega}^{k+1}$ and $S_{B}^{k+1}$ for the next
iteration. The whole algorithm is shown in Algorithm 1.

Algorithm 1 Adaptive sampling method for Deep Ritz

0: Initial $u(x,\theta^{0})$, $p_{\text{bKRnet}}(x,\theta_{f}^{0})$, maximum
epoch number $N_{e}$, batch size $m$, initial Deep Ritz training set
$S_{\Omega}^{0}$ and $S_{\partial\Omega}^{0}$, initial bounded KRnet training
set $S_{B}^{0}$

1: for $k=0$ to $N_{\text{adaptive}}-1$ do

2: // Solve PDE by Deep Ritz method

3: for $i=0$ to $N_{e}$ do

4: for $j=0$ to $J$ do

5: Sample $m$ samples from $S_{\Omega}^{k}$.

6: Sample $m$ samples from $S_{\partial\Omega}^{k}$.

7: Update $u(x,\theta)$ by descending the stochastic gradient of the objective
in $(\ref{A1})$.

8: end for

9: end for

10: // Train bounded KRnet

11: for $i=0$ to $N_{e}$ do

12: for $j=0$ to $J$ do

13: Sample $m$ samples from $S_{B}^{k}$

14: Update $p_{\text{bKRnet}}(x,\theta_{f})$ by descending the stochastic
gradient of the objective in $(\ref{A2})$.

15: end for

16: end for

17: // Generate new training set

18: Generate $S_{\Omega}^{k+1}$ and $S_{B}^{k+1}$ using
$p_{\text{bKRnet}}(x,\theta_{f}^{k+1})$.

19: end for

####  3.3.2 Alternative PDF models and sample generation

In practice, using the numerical approximation of the optimal PDF for the
change of measure may cause a problem that the ratio
$\frac{G(u(x,\theta))}{p(x)}$ in equation (3) is too large for some samples.
This is possible since the numerical approximation of $|G(u(x))|/\mu$
introduces uncertainties where we do not have any control over the ratio. We
then consider two alternative PDF models to enhance both the robustness and
effectiveness:

  1. 1.

We build a mixture model by combining the PDF induced by bounded KRnet and the
PDF of the uniform distribution on the domain.

  2. 2.

We build a mixture model by combining the PDF induced by bounded KRnet in the
current adaptivity iteration and the mixture model in the previous adaptivity
iteration.

We look at the first alternative PDF model, which is a mixture model defined
as

$p_{\text{mixture}}(x,\theta_{f}^{k+1})=\epsilon\,p_{\text{bKRnet}}(x,\theta_{f}^{k+1})+(1-\epsilon)\,p_{\text{uniform}}(x),\quad\forall
x\in S_{\Omega}^{k+1},$ (14)

where $0<\epsilon<1$. In contrast to the basic algorithm, a portion of the
training points are from the uniform distribution instead of bounded KRnet.
Noting that

$p_{\text{mixture}}(x,\theta_{f}^{k+1})\geq\frac{1-\epsilon}{|\Omega|},$

we then avoid the problem that the ratio
$\frac{G(u(x,\theta_{f}^{k+1}))}{p_{\text{mixture}}(x,\theta_{f}^{k+1})}$
becomes too large. The same approach can also be applied to the training set
$S_{B}^{k}$. The algorithm is presented in Algorithm 2.

Algorithm 2 Adaptive sampling with alternative PDF model $1$

0: Initial $u(x,\theta^{0})$, $p_{\text{bKRnet}}(x,\theta_{f}^{0})$, maximum
epoch number $N_{e}$, batch size $m$, initial Deep Ritz training set
$S_{\Omega}^{0}$ and $S_{\partial\Omega}^{0}$, initial bounded KRnet training
set $S_{B}^{0}$, mixture parameter $\epsilon$

1: for $k=0$ to $N_{\text{adaptive}}-1$ do

2: // Solve PDE by Deep Ritz method

3: for $i=0$ to $N_{e}$ do

4: for $j=0$ to $J$ do

5: Sample $m$ samples from $S_{\Omega}^{k}$.

6: Sample $m$ samples from $S_{\partial\Omega}^{k}$.

7: Update $u(x,\theta)$ by descending the stochastic gradient of the objective
in $(\ref{A1})$.

8: end for

9: end for

10: // Train bounded KRnet

11: for $i=0$ to $N_{e}$ do

12: for $j=0$ to $J$ do

13: Sample $m$ samples from $S_{B}^{k}$

14: Update $p_{\text{bKRnet}}(x,\theta_{f})$ by descending the stochastic
gradient of the objective in $(\ref{A2})$.

15: end for

16: end for

17: // Generate new training set

18: Generate uniform training sets $S_{\Omega,\text{uniform}}^{k+1}$ and
$S_{B,\text{uniform}}^{k+1}$ using the uniform distribution on $\Omega$.

19: Generate training sets $S_{\Omega,\text{bKRnet}}^{k+1}$ and
$S_{B,\text{bKRnet}}^{k+1}$ using $p_{\text{bKRnet}}(x,\theta_{f}^{k+1})$.

20: Set $S_{\Omega}^{k+1}=S_{\Omega,\text{uniform}}^{k+1}\cup
S_{\Omega,\text{bKRnet}}^{k+1}$ and
$S_{B}^{k+1}=S_{B,\text{uniform}}^{k+1}\cup S_{B,\text{bKRnet}}^{k+1}$.

21: Calculate the PDF values of training sets $S_{\Omega}^{k+1}$ and
$S_{B}^{k+1}$ using $(\ref{alternative_model_1})$.

22: end for

We now present the second alternative PDF model, which is defined as

$\displaystyle p_{\text{mixture}}(x,\left\\{\theta_{f}^{t},t\leq
k+1\right\\})$ (15) $\displaystyle=$
$\displaystyle\epsilon\,p_{\text{bKRnet}}(x,\theta_{f}^{k+1})+(1-\epsilon)\,p_{\text{mixture}}(x,\left\\{\theta_{f}^{t},t\leq
k\right\\})$ $\displaystyle=$
$\displaystyle\epsilon\,p_{\text{bKRnet}}(x,\theta_{f}^{k+1})$
$\displaystyle+(1-\epsilon)\left(\epsilon\,p_{\text{bKRnet}}(x,\theta_{f}^{k})+(1-\epsilon)\,p_{\text{mixture}}(x,\left\\{\theta_{f}^{t},t\leq
k-1\right\\})\right)$ $\displaystyle=$ $\displaystyle\cdots$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{k+1}\epsilon(1-\epsilon)^{k+1-t}\,p_{\text{bKRnet}}(x,\theta_{f}^{t})+(1-\epsilon)^{k+1}/|\Omega|$

where we form a mixture model using the PDF model
$p_{\text{mixture}}(x,\left\\{\theta_{f}^{t},t\leq k\right\\})$ at adaptivity
iteration $k$ and the PDF model $p_{\text{bKRnet}}(x,\theta_{f}^{k+1})$ at
adaptivity iteration $k+1$. The initial data set $S_{\Omega}^{0}$ is sampled
from a uniform distribution. As shown in equation (15), the training points
from the last mixture model and their PDF values depend on all the previous
bounded KRnets. The value of the proposed mixture model is at least
$(1-\epsilon)^{k+1}/|\Omega|$. The same manipulation can also be implemented
to the training set $S_{B}^{k+1}$.

The second alternative PDF model is an extension of the first one. In the
second alternative PDF model, all previous bounded KRnets are saved and reused
every iteration. Although method $2$ needs more storage and computation cost,
it is more flexible than method $1$ and the extra cost can be neglected
compared to the total cost. The algorithm is presented in Algorithm 3.

Algorithm 3 Adaptive sampling with alternative PDF model $2$

0: Initial $u(x,\theta^{0})$, $p_{\text{bKRnet}}(x,\theta_{f}^{0})$, maximum
epoch number $N_{e}$, batch size $m$, initial Deep Ritz training set
$S_{\Omega}^{0}$ and $S_{\partial\Omega}^{0}$, initial bounded KRnet training
set $S_{B}^{0}$, mixture parameter $\epsilon$

1: for $k=0$ to $N_{\text{adaptive}}-1$ do

2: // Solve PDE by Deep Ritz method

3: for $i=0$ to $N_{e}$ do

4: for $j=0$ to $J$ do

5: Sample $m$ samples from $S_{\Omega}^{k}$.

6: Sample $m$ samples from $S_{\partial\Omega}^{k}$.

7: Update $u(x,\theta)$ by descending the stochastic gradient of the objective
in $(\ref{A1})$.

8: end for

9: end for

10: // Train bounded KRnet

11: for $i=0$ to $N_{e}$ do

12: for $j=0$ to $J$ do

13: Sample $m$ samples from $S_{B}^{k}$

14: Update $p_{\text{bKRnet}}(x,\theta_{f})$ by descending the stochastic
gradient of the objective in $(\ref{A2})$.

15: end for

16: end for

17: // Save bounded KRnet

18: Restore current bound KRnet parameters $\theta_{f}^{k+1}$.

19: // Generate new training set

20: Initialize new training sets with uniform training sets
$S_{\Omega}^{k+1}=S_{\Omega,\text{uniform}}^{k+1}$ and
$S_{B}^{k+1}=S_{B,\text{uniform}}^{k+1}$.

21: for $t=0$ to $k+1$ do

22: Load bounded KRnet with parameters $\theta_{f}^{t+1}$.

23: Generate training sets $S_{\Omega,\text{bKRnet}}^{t+1}$,
$S_{B,\text{bKRnet}}^{t+1}$ using $p_{\text{bKRnet}}(x,\theta_{f}^{t+1})$.

24: Set $S_{\Omega}^{k+1}=S_{\Omega}^{k+1}\cup S_{\Omega,\text{bKRnet}}^{t+1}$
and $S_{B}^{k+1}=S_{B}^{k+1}\cup S_{B,\text{bKRnet}}^{t+1}$ .

25: end for

26: Calculate the PDF values of training sets $S_{\Omega}^{k+1}$ and
$S_{B}^{k+1}$ using $(\ref{method_2})$.

27: end for

##  4 Numerical experiments

In this section, we present some numerical experiments to illustrate the
effectiveness of our methods. We test four examples; among these examples, two
have a solution with low regularity, another has a singularity in the
solution, and the last one is a high-dimensional problem. For the Deep Ritz
method, we use a deep network structure similar to ResNet [4]. In [2], the
authors compared the performance of several different activation functions.
Their results show that $\sin^{3}(x)$ is the best among all tested activation
functions. Hence, we use $\sin^{3}(x)$ as the activation function in all
experiments. For bounded KRnet, we use the same architecture as described in
[16]. For both Deep Ritz method and bounded KRnet, we train the neural
networks by the ADAM method. In order to test the validity of our method, we
use the following $L_{2}$ error

$\text{Error}=\frac{\sqrt{\sum_{i=1}^{N}|u(x_{i},\theta)-u(x_{i})|^{2}}}{\sqrt{\sum_{i=1}^{N}|u(x_{i})|^{2}}}$

in our experiments where the set of test points is either generated using a
tensor grid or drawn from a uniform distribution (for high dimensional
problem). In the original Deep Ritz method, the training points are reselected
from a uniform distribution at each training epoch. Hence, the concept total
number of training points is not explicitly defined for the Deep Ritz method
since it depends on the training process. One typically uses the batch size
(the number of training points used in one training epoch) to describe the
scale of the training set in the Deep Ritz method. Similar to the original
Deep Ritz method, in our adaptive approach, new training points can be
efficiently sampled from a bounded KRnet at each training epoch for the basic
algorithm and the first alternative PDF model. For the second alternative PDF
model, the dynamic sampling is not efficient since all previous models are
involved. In this work, at the $k$th adaptive iteration, we generate a finite
training set $S_{\Omega}^{k}$($S_{B}^{k}$) of size $N_{v}$($N_{B}$), based on
which mini-batches are formed for each training epoch. Across all numerical
experiments, we compare the performance of the following four sampling
strategies: (1) The conventional Deep Ritz, where the training set will be
completely updated with new uniform samples at each adaptivity iteration; (2)
The basic adaptive algorithm based on bounded KRnet, i.e., Algorithm 1; (3)
The adaptive algorithm based on the first alternative PDF model, i.e.,
Algorithm 2; and (4) The adaptive algorithm based on the second alternative
PDF model, i.e., Algorithm 3. For a fair comparison, we use the same uniformly
distributed training set, denoted as $S_{\Omega}^{0}$ for Deep Ritz at
iteration $0$ across all four sampling strategies. Similarly, we employ the
same uniformly distributed training set, denoted as $S_{B}^{0}$, to train
bounded KRnet at iteration $0$ for all three bKRnet-based methods.

###  4.1 Two dimensional peak problem

We consider the following 2D elliptic problem with a single peak:

$\displaystyle-\Delta u(x_{1},x_{2})$ $\displaystyle=f(x_{1},x_{2}),$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\Omega$ $\displaystyle
u(x_{1},x_{2})$ $\displaystyle=g(x_{1},x_{2}),$ $\displaystyle(x_{1},x_{2})$
$\displaystyle\in\partial\Omega$

where the computational domain is $\Omega=[-1,1]^{2}$. The exact solution is
given by

$u(x_{1},x_{2})=\exp(-1000((x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}))$

which has a peak at $(0.5,0.5)$ and decreases rapidly away from $(0.5,0.5)$.
This is a classical benchmark problem commonly used to assess the performance
of adaptive techniques in both deep learning methods [12, 3] and finite
element methods [6].

(a) Calculated solution, True solution, Absolute error Figure 1: Numerical
results of the alternative model $1$ with $\epsilon=0.5$: Numerical solution,
Exact solution, Absolute error at iteration $7$.

We choose a ResNet $u(x,\theta)$ with a stack of $4$ blocks(eight fully-
connected layers) and $32$ neurons in each layers to approximate the solution.
For bounded KRnet, we use $6$ CDF layers and two fully connected layers with
$24$ neurons for each CDF layers. The training process for the Deep Ritz
method consists of $10000$ epochs for the initial iteration, followed by an
additional $1000$ epochs for each subsequent iteration. The number of epochs
for training the bounded KRnet is $3000$. The balance parameter $\beta$ in the
Deep Ritz method is configured to be $100$. The learning rate for ADAM
operator is set to $0.0001$ and the batch size is set to $m=5000$. We test our
problem on a $101\times 101$ tensor grid.

We first apply the alternative PDF model $1$ with $\epsilon=0.5$ to the Deep
Ritz method. The numerical results are presented in Figure 1. We show the
evolution of the training sets with respect to adaptivity iterations
$k=1,3,5,7$ in Figure 2. It can be observed that the peak density in each
training set consistently occurs around the coordinates $[0.5,0.5]$. This is
because we use bKRnet to approximate the normalized absolute
integrand(computed from the Deep Ritz solution), and its shape remains
relatively stable during the training process. Note that this is different
from the distributions of training points generated by the residual based
DAS-R for PINNs. The distribution of training points generated by DAS-R
becomes less concentrated as the number of adaptive iteration increases.
Another difference is that computing the residual during bKRnet training in
DAS-R requires second order derivatives while only first order derivative are
needed in our approach.

(a) Training points at iteration $1$

(b) Training points at iteration $3$

(c) Training points at iteration $5$

(d) Training points at iteration $7$

Figure 2: Training points of the alternative model $1$ with $\epsilon=0.5$

(a) Alternative model $1$ with $\epsilon=0.5$

(b) Alternative model $2$ with $\epsilon=0.3$

(c) Uniform

(d) Pure bounded KRnet

Figure 3: Comparison of numerical solutions using different sampling schemes

We compare the performances of different adaptive sampling strategies in
Figure 3. The size of the training set is 200000. It is evident that all
strategies based on importance sampling can capture information in the region
of low regularity more effectively than the original Deep Ritz method.

In Figure 4, we present the $L_{2}$ errors for four sampling strategies. As
all adaptive sampling methods require replacing all training points at each
iteration, the relative $L_{2}$ errors do not exhibit a monotonic decrease.
However, a decreasing trend can still be observed, allowing us to compare the
performances of different methods. Among all the methods, the alternative
model $1$ provides the most accurate solution. It is also worth noting that
the alternative method $2$ provides a reasonably accurate solution at
iteration $3$ which means we can stop our adaptive method earlier to save the
computational cost.

(a) Relative $L_{2}$ error Figure 4: Relative $L_{2}$ errors

###  4.2 Two dimensional problem with two peaks

We consider the following two dimensional problem with two peaks

$\displaystyle-\nabla^{2}u(x_{1},x_{2})+\nabla\cdot[u(x_{1},x_{2})\nabla(x_{1}^{2}+x_{2}^{2})]$
$\displaystyle=s(x_{1},x_{2}),$ $\displaystyle(x_{1},x_{2})$
$\displaystyle\in\Omega$ $\displaystyle u(x_{1},x_{2})$
$\displaystyle=g(x_{1},x_{2}),$ $\displaystyle(x_{1},x_{2})$
$\displaystyle\in\partial\Omega.$

where the computation domain is $\Omega=[-1,1]^{2}$. We choose the exact
solution of the problem to be

$\displaystyle u(x_{1},x_{2})=$
$\displaystyle\exp(-1000[(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}])$
$\displaystyle+\exp(-1000[(x_{1}+0.5)^{2}+(x_{2}+0.5)^{2}])$

which has two peaks at the points $(0.5,0.5)$ and $(-0.5,-0.5)$. The Dirichlet
boundary condition $g(x_{1},x_{2})$ and the right-hand side term are also
given by the exact solution.

(a) Calculated solution, True solution, Absolute error Figure 5: Numerical
results of the alternative model $2$ with $\epsilon=0.1$: Numerical solution,
Exact solution, Absolute error at iteration $7$

We apply the same ResNet and bounded KRnet as in the one peak problem. We also
utilize the same optimizer and learning rate for training these deep neural
networks. However, in this experiment, we use fewer training epochs and
points.

(a) Training points at iteration $1$

(b) Training points at iteration $3$

(c) Training points at iteration $5$

(d) Training points at iteration $7$

Figure 6: Training points of the alternative model $2$ with $\epsilon=0.1$

We aim to find out if our adaptive sampling methods can compute a reasonably
accurate solution when the uniform sampling strategy fails. In this
experiment, the training process for the Deep Ritz method consists of only
$3000$ epochs for the initial iteration, followed by an additional $100$
epochs for each subsequent iteration. The number of epochs for training the
bounded KRnet is reduced to $1000$. We have reduced the size of the training
set to $20000$ and the batch size to $2000$.

We present the numerical solution and the evolution of training points of the
adaptive sampling method based on alternative PDF model $2$ in Figure 5 and
Figure 6 respectively. It is seen that our adaptive methods can effectively
capture the information around the two peaks and therefore result in a
solution with reasonable accuracy. We compare the performances of four
different adaptive sampling strategies in Figure 7. Due to the complexity of
the two peak problem and the reduced number of training epochs and points, the
uniform sampling strategy produces an incorrect numerical solution. In
contrast, all three adaptive methods yield reasonably accurate solutions.
Among these methods, the two alternative models outperform the pure bounded
KRnet.

(a) Alternative model $1$ with $\epsilon=0.5$

(b) Alternative model $2$ with $\epsilon=0.1$

(c) Uniform

(d) Pure bounded KRnet

Figure 7: Comparison of the absolute errors of numerical solution using
different sampling schemes

###  4.3 Two-dimensional problem with singularity

We consider the following 2D elliptic problem

$\displaystyle-$ $\displaystyle\Delta u(x_{1},x_{2})=0,$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\Omega$ $\displaystyle
u(x_{1},x_{2})=u(r,\theta)=r^{\frac{1}{2}}\sin{\frac{\theta}{2}},$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\partial\Omega$

where the domain
$\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times\left\\{0\right\\}$. The
boundary is shown in the Figure 8. The solution to this problem suffers from
the ’corner singularity’ caused by the nature of the domain [11]. The explicit
solution $u^{*}(x)=r^{\frac{1}{2}}\sin\frac{\theta}{2}$ in polar coordinates.

Figure 8: Training points on the boundary of an elliptic problem with
singularity.

Considering the problem in the coordinates of $x_{1}$ and $x_{2}$, we obtain
the following PDE and its associated boundary conditions:

$\displaystyle-$ $\displaystyle\Delta u(x_{1},x_{2})=0,$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\Omega,$ $\displaystyle
u(x_{1},x_{2})=(x_{1}^{2}+x_{2}^{2})^{\frac{1}{4}}\sin\left({\frac{1}{2}\arctan\frac{x_{2}}{x_{1}}}\right),$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\partial\Omega\cap Q_{1},$
$\displaystyle
u(x_{1},x_{2})=(x_{1}^{2}+x_{2}^{2})^{\frac{1}{4}}\sin\left(\pi+{\frac{1}{2}\arctan\frac{x_{2}}{x_{1}}}\right),$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\partial\Omega\cap Q_{2},$
$\displaystyle
u(x_{1},x_{2})=(x_{1}^{2}+x_{2}^{2})^{\frac{1}{4}}\sin\left(\pi+{\frac{1}{2}\arctan\frac{x_{2}}{x_{1}}}\right),$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\partial\Omega\cap Q_{3},$
$\displaystyle
u(x_{1},x_{2})=(x_{1}^{2}+x_{2}^{2})^{\frac{1}{4}}\sin\left(2\pi+{\frac{1}{2}\arctan\frac{x_{2}}{x_{1}}}\right),$
$\displaystyle(x_{1},x_{2})$ $\displaystyle\in\partial\Omega\cap Q_{4}$

where $x_{1}\neq 0$ and $Q_{1},\cdots,Q_{4}$ represent the four quadrants in
2D space. We then use a ResNet $u(x,\theta)$ to approximate the solution
$u(x)$. In the experiment, we choose a ResNet $u(x,\theta)$ with a stack of
$4$ blocks(eight fully-connected layers) and $10$ neurons in each layers. For
bounded KRnet, we use $6$ CDF layers and two fully connected layers with $24$
neurons for each CDF layers. We implement $9$ adaptivity iterations for this
problem The number of epochs for training bounded KRnet is set to $3000$ and
the number of epochs for training $u(x,\theta)$ is chosen as $20000$ for the
first iteration and followed by an additional $1000$ epochs for each
subsequent iteration. The balance parameter $\beta$ in the Deep Ritz method is
configured to be $1000$. The learning rate for ADAM operator is set to $0.001$
and the batch size is set to $m=2000$. In order to avoid testing the
singularity at the origin, we test our problem on a $102\times 102$ tensor
grid.

(a) Iteration $0$

(b) Iteration $2$

(c) Iteration $4$

(d) Iteration $6$

(e) Iteration $8$

Figure 9: The evolution of the normalized integrand of the loss functional.

(a) Iteration $0$

(b) Iteration $2$

(c) Iteration $4$

(d) Iteration $6$

(e) Iteration $8$

Figure 10: The evolution of the training points generated from the bounded
KRnet.

(a) Iteration $0$

(b) Iteration $2$

(c) Iteration $4$

(d) Iteration $6$

(e) Iteration $8$

Figure 11: The evolution of the absolute errors.

We present a set of numerical results using the pure KRnet sampling approach
with a large training set containing $200000$ points. In Figure 9, we display
the evolution of the normalized integrand during the training process. It can
be observed that the peak (where the singularity occurs) becomes sharper as
the number of iterations increases. Since the bounded KRnet learns a
distribution from the normalized integrand, the samples generated from the
bounded KRnet become increasingly concentrated around the central peak as
shown in Figure 10. In contrast to the problem with a low regularity solution,
our adaptive method recognizes that sampling from the region around the
singularity becomes increasingly important as the number of iterations grows.
Figure 11 illustrates the evolution of the absolute error during the training
process. The absolute error around the singularity decreases as the iterations
progress.

(a) Uniform

(b) Alternative model $1$ with $\epsilon=0.5$

(c) Alternative model $2$ with $\epsilon=0.1$

(d) Pure bounded KRnet

Figure 12: Comparison of three adaptive sampling methods. From left to right:
numerical solution, exact solution, absolute error

Figure 13: Comparison of $L_{2}$ errors

We compare the absolute errors and $L_{2}$ errors of the four adaptive
sampling methods in Figure 12 and Figure 13. Unlike the problem with a low
regularity solution, the uniform samples based method works well for this
problem. However, its performance is not as good as the other three methods
based on deep generative models. The absolute error of the uniform method is
almost twice as large as the error of the alternative model $1$ based method.
From iteration $3$, the errors of alternative model $1$ based method and pure
bounded KRnet based method drop faster than that of the uniform samples based
method.

###  4.4 High dimensional Poisson problem

We consider the following $10$D elliptic equation

$\displaystyle-\Delta u(x)$ $\displaystyle=f(x),$ $\displaystyle
x\in\Omega=[-1,1]^{10}$

with an exact solution

$u(x)=e^{-10\|x\|_{2}^{2}}.$

where the Dirichlet boundary condition on $\partial\Omega$ is given by the
exact solution. In a high-dimensional hypercube, a large portion of uniformly
distributed samples tend to cluster near the surface of the hypercube.
However, from the exact solution of this 10D Poisson problem, we know the
important computation region is mainly in the neighborhood of the origin.
Hence, if we employ uniformly distributed samples to train the neural
networks, the majority of these samples will contribute little to error
reduction.

(a) Calculated solution, True solution, Absolute error Figure 14: Numerical
result: Calculated solution, True solution, Absolute error on first two
dimension at iteration $7$

We use a ResNet $u(x,\theta)$ with a stack of $4$ blocks(eight fully-connected
layers) and $64$ neurons in each layer. For bounded KRnet, we use $8$ CDF
layers and two fully connected layers with $32$ neurons for each CDF layers.
We conduct a total of $7$ adaptive iterations for this problem. The number of
epochs for training bounded KRnet is set to $2000$, while for training
$u(x,\theta)$, we initially use $10000$ epochs for the first iteration,
followed by an additional $1000$ epochs for each subsequent iteration. The
balance parameter $\beta$ in the Deep Ritz method is selected to be $1000$.
The learning rate for ADAM operator is set to $0.001$ and the batch size is
set to $m=20000$. We apply the alternative PDF model $1$ with $\epsilon=0.5$.
We conduct our tests on a test set with $1000000$ points sampled from a
uniform distribution on $\Omega$.

We present the numerical solution and the absolute error in Figure 14. We
illustrate the evolution of training points for iterations $1$,$3$, $5$, and
$7$, focusing on the components $x_{1}$ and $x_{2}$ for visualization
purposes. We also examined the other components and did not observe
significantly different results. Based on the numerical results and the
distribution of training points, our adaptive method effectively captures the
crucial computation region around the origin and therefore reduces the error
peak.

(a) Training points on first two dimension at iteration $1$

(b) Training points on first two dimension at iteration $3$

(c) Training points on first two dimension at iteration $5$

(d) Training points on first two dimension at iteration $7$

Figure 15: The evolution of training points on first two dimensions

We compare the $L_{2}$ errors of four different adaptive sampling strategies
in Figure 16. It can be seen that the $L_{2}$ error of the uniform sampling
strategy does not decrease as the iteration progresses. The $L_{2}$ error of
the other three adaptive methods exhibits decreasing trends. Among them, the
alternative method $2$ with $\epsilon=0.1$ performs best. It achieves the
smallest relative $L_{2}$ error $5.2451e-02$ at iteration $4$. Both
alternative model based methods can be terminated earlier at iteration $4$ to
reduce computational cost.

Figure 16: Relative $L_{2}$ errors

##  5 Conclusion

In this work, we propose an adaptive sampling strategy for the Deep Ritz
method. We use a deep neural network to approximate the solution of the PDE
and a deep generative model called bounded KRnet to approximate the integrand
in the variational loss of the Deep Ritz method. Next, we employ the bounded
KRnet to generate new samples and calculate the associated probability density
function (PDF) values for these samples. The variational loss can then be more
accurately calculated using importance sampling and therefore a more precise
solution can be obtained from the Deep Ritz method. We present four numerical
experiments to illustrate the performance of our adaptive sampling strategy.
These numerical experiments include two problems with low regularity, a
problem with singularity and a high dimensional problem. In all four
experiments, our adaptive sampling strategy consistently achieves a more
accurate solution compared to the uniform sampling strategy used in the
original Deep Ritz method. We also develop two alternative PDF models based on
bounded KRnet. In our experiments, these alternative PDF models can outperform
the PDF directly from a bounded KRnet. There are several possible ways to
further improve our adaptive sampling strategy based on bounded KRnet. First,
in this work, our adaptive sampling strategies require to replace entire
training set at each adaptive iteration. The disadvantage of this approach is
that the relative $L_{2}$ errors do not monotonically decrease. For PINNs,
gradually adding training points at each iteration is an effective way to
mitigate $L_{2}$ error oscillations caused by the randomness associated with
replacing entire training set [12]. Therefore, one potential approach to
improve current adaptive strategies is to develop a suitable adding points
strategy for Deep Ritz method. Second, other variance reduction techniques
such as control variates can be taken into consideration to further improve
the approximation of the variational loss. Research on these issues will be
reported in future studies.

##  6 Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * [1] Guillaume Carlier, Alfred Galichon and Filippo Santambrogio  “From Knothe’s transport to Brenier’s map and a continuation method for optimal transport”  In _SIAM Journal on Mathematical Analysis_ 41.6 SIAM, 2010, pp. 2554–2576
  * [2] Jingrun Chen, Rui Du and Keke Wu  “A comparison study of deep Galerkin method and deep Ritz method for elliptic problems with different boundary conditions”  In _arXiv preprint arXiv:2005.04554_ , 2020
  * [3] Zhiwei Gao, Liang Yan and Tao Zhou  “Failure-informed adaptive sampling for PINNs”  In _SIAM Journal on Scientific Computing_ 45.4 SIAM, 2023, pp. A1971–A1994
  * [4] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun  “Deep residual learning for image recognition”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778
  * [5] Lu Lu, Xuhui Meng, Zhiping Mao and George Em Karniadakis  “DeepXDE: A deep learning library for solving differential equations”  In _SIAM review_ 63.1 SIAM, 2021, pp. 208–228
  * [6] Pedro Morin, Ricardo H Nochetto and Kunibert G Siebert  “Convergence of adaptive finite element methods”  In _SIAM review_ 44.4 SIAM, 2002, pp. 631–658
  * [7] George Papamakarios et al.  “Normalizing flows for probabilistic modeling and inference”  In _The Journal of Machine Learning Research_ 22.1 JMLRORG, 2021, pp. 2617–2680
  * [8] Maziar Raissi, Paris Perdikaris and George E Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational physics_ 378 Elsevier, 2019, pp. 686–707
  * [9] Filippo Santambrogio  “Optimal transport for applied mathematicians”  In _Birkäuser, NY_ 55.58-63 Springer, 2015, pp. 94
  * [10] Justin Sirignano and Konstantinos Spiliopoulos  “DGM: A deep learning algorithm for solving partial differential equations”  In _Journal of computational physics_ 375 Elsevier, 2018, pp. 1339–1364
  * [11] Gilbert Strang, George J Fix and DS Griffin  “An analysis of the finite-element method”, 1974
  * [12] Kejun Tang, Xiaoliang Wan and Chao Yang  “DAS-PINNs: A deep adaptive sampling method for solving high-dimensional partial differential equations”  In _Journal of Computational Physics_ 476 Elsevier, 2023, pp. 111868
  * [13] Chenxi Wu et al.  “A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks”  In _Computer Methods in Applied Mechanics and Engineering_ 403 Elsevier, 2023, pp. 115671
  * [14] Bing Yu  “The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems”  In _Communications in Mathematics and Statistics_ 6.1 Springer, 2018, pp. 1–12
  * [15] Yaohua Zang, Gang Bao, Xiaojing Ye and Haomin Zhou  “Weak adversarial networks for high-dimensional partial differential equations”  In _Journal of Computational Physics_ 411 Elsevier, 2020, pp. 109409
  * [16] Li Zeng, Xiaoliang Wan and Tao Zhou  “Bounded KRnet and its applications to density estimation and approximation”  In _arXiv preprint arXiv:2305.09063_ , 2023