# VS-PINN: A Fast and efficient training of physics-informed neural networks
using variable-scaling methods for solving PDEs with stiff behavior

Seungchan Ko and Sang Hyeon Park Department of Mathematics, Inha University,
Incheon, Republic of Korea. Email<EMAIL_ADDRESS>of Mathematics,
Inha University, Incheon, Republic of Korea. Email<EMAIL_ADDRESS>
( )

###### Abstract

Physics-informed neural networks (PINNs) have recently emerged as a promising
way to compute the solutions of partial differential equations (PDEs) using
deep neural networks. However, despite their significant success in various
fields, it remains unclear in many aspects how to effectively train PINNs if
the solutions of PDEs exhibit stiff behaviors or high frequencies. In this
paper, we propose a new method for training PINNs using variable-scaling
techniques. This method is simple and it can be applied to a wide range of
problems including PDEs with rapidly-varying solutions. Throughout various
numerical experiments, we will demonstrate the effectiveness of the proposed
method for these problems and confirm that it can significantly improve the
training efficiency and performance of PINNs. Furthermore, based on the
analysis of the neural tangent kernel (NTK), we will provide theoretical
evidence for this phenomenon and show that our methods can indeed improve the
performance of PINNs.

Keywords: Physics-informed neural networks, variable scaling, stiff behavior,
high frequency, spectral bias, neural tangent kernel

##  1 Introduction

Bridging the gap between classical scientific computing and machine learning,
the emerging field called scientific machine learning has introduced a
completely different framework to compute the solutions of partial
differential equations (PDEs). The forefront of this evolution lies within the
area of physics-informed neural networks (PINNs) [18, 10]. Based on the
universal approximation property of deep neural networks, the PINNs
approximate the solutions of PDEs by incorporating underlying physics into
deep neural networks. Due to their simplicity and flexibility in handling a
wide range of physical problems involving PDEs, PINNs have recently gained
great attention and have been applied to various fields in computational
science: bio-medical science [22, 11], fluids mechanics [20, 19, 9, 4, 21],
uncertainty quantification [30, 32, 29] and meta-material design [12, 3].
Moreover, since the PINNs utilize randomly-selected collocation points as
training samples in the spatio-temporal domain, the PINNs are available for
high-dimensional PDEs [24, 6], on the domains with complex geometries [13, 23,
26, 14].

However, despite the significant empirical success of PINNs, we only have
limited knowledge about the behavior of these constrained neural networks
during their training process, and the training of PINNs often fails. In
particular, since neural networks typically assume a smooth prior, it is often
challenging to train PINNs to learn a solution with a sharp transition which
poses significant obstacles for the model prediction. For example, some recent
studies have illustrated that training PINNs with fully-connected neural
networks usually suffers from so-called spectral bias or F-principle meaning
that it is difficult for PINNs to learn functions with high frequencies [17,
2, 25]. Furthermore, it has recently been observed that the PINNs have some
problems regarding convergence and accuracy when the target functions have
sharp spatial transitions or fast temporal evolution [27]. Some previous works
investigated the ways to handle these issues. For example, [28] theoretically
analyzed the training process of PINNs based on the recently developed neural
tangent kernel (NTK) [8], which explores the connection between deep neural
networks and the kernel regression, elucidating the training dynamics. From
the intuition obtained in the analysis, the authors proposed an NTK-guided
gradient-descent algorithm to efficiently train PINNs without spectral-bias
pathology and achieved remarkable improvements on various benchmark problems.
In order to handle stiff PDEs, the authors of [15] developed Self-Adaptive
PINNs, a novel paradigm for training PINNs using trainable weights based on
the attention mechanism. On the other hand, for the singular perturbation
problem which is another area in stiff PDEs, [5] developed a semi-analytic
PINN method, based on the so-called corrector functions obtained from the
boundary layer analysis.

In contrast to these methods, the present paper aims to provide a completely
different method for tackling the aforementioned problems. In particular, we
will propose a novel training method for PINNs using variable-scaling, which
we have named VS-PINNs (variable-scaling physics-informed neural networks).
Scaling is a common technique in the analysis of PDEs, but to the best of our
knowledge, the idea has never been applied to train PINNs. Throughout various
numerical experiments, we will show that this method can greatly improve the
performance of PINNs and it is very effective in solving stiff problems. As it
will be made clear soon, our method is simple and no significant extra
computational costs are incurred. Furthermore, it also offers computational
flexibility, resulting in significant advancements in accuracy and learning
efficiency for various computational problems.

The rest of the paper is organized as follows. In Section 2, we will collect
some preliminary materials including brief descriptions of the PINNs and the
NTK, which will be considered throughout the paper. In Section 3, we shall
describe the method that we mainly propose in this paper. And then in Section
4, we will perform comprehensive experiments to evaluate the performance of
the proposed method. Throughout the numerical simulation of a wide range of
PDEs with stiff behavior, high frequency and even for nonlinear problems, we
will illustrate the effectiveness of the method we propose. Once we
empirically confirm the better performance of the new method, in Section 5, we
will also rigorously investigate the proposed method by using the NTK. By
examining the eigenvalues of the NTK of VS-PINNs, we will provide some
theoretical evidence on why the proposed method can indeed improve the
performance of PINNs. Finally in Section 6, we will make a concluding remark,
addressing the limitation and further improvement of the method, as well as
the future research direction regarding the VS-PINNs.

##  2 Preliminary

###  2.1 Physics-Informed Neural Networks

Figure 1: Schematic illustration of the physics-informed neural network. The
left part visualizes a standard feed-forward neural network parameterized by
$\theta$, while the right part imposes the given physical laws to the neural
network.

In this section, we will briefly review the concept of physics-informed neural
networks (PINN) [18]. For the description of the method, let us consider the
following problem on a bounded open domain $\Omega\subset\mathbb{R}^{d}$:

$\displaystyle\mathcal{D}[u](\boldsymbol{x})$
$\displaystyle=f(\boldsymbol{x}),$ $\displaystyle\boldsymbol{x}\in\Omega,$
(2.1) $\displaystyle u(\boldsymbol{x})$ $\displaystyle=g(\boldsymbol{x}),$
$\displaystyle\boldsymbol{x}\in\partial\Omega,$

where $\mathcal{D}$ is a differential operator and
$u:\overline{\Omega}\rightarrow\mathbb{R}$ is the unknown function with the
variable $\boldsymbol{x}=(x_{1},x_{2},\cdots,x_{d})\in\Omega$. For the sake of
simplicity, for time-dependent problems, we regard the temporal variable $t$
as an additional coordinate of $\boldsymbol{x}$. Hence $\Omega$ may denote the
spatio-temporal domain, and the initial condition can be regarded as a special
type of boundary data on $\partial\Omega$. The main idea of the PINN is to
approximate the solution $u(\boldsymbol{x},t)$ with a deep neural network
$u(\boldsymbol{x},t;\boldsymbol{\theta})$, where $\boldsymbol{\theta}$ is a
collection of all parameters of the neural network. To train the model, we
minimize the so-called physics-informed loss function consisting of the
residual of the given PDE in the domain and the differences of the neural
network with the imposed conditions on the boundary. More precisely, the
neural network is trained by minimizing the total loss

$\mathcal{L}_{\rm total}=\lambda_{r}\mathcal{L}_{\rm
res}+\lambda_{d}\mathcal{L}_{\rm data},$ (2.2)

where each loss function is defined as follows:

$\displaystyle\mathcal{L}_{\rm res}$
$\displaystyle=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\left|\mathcal{D}[u](\boldsymbol{x}_{r}^{i};\boldsymbol{\theta})-f(\boldsymbol{x}_{r}^{i})\right|^{2},$
(2.3) $\displaystyle\mathcal{L}_{\rm data}$
$\displaystyle=\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left|u(\boldsymbol{x}_{b}^{j};\boldsymbol{\theta})-g(\boldsymbol{x}_{b}^{j})\right|^{2}.$
(2.4)

Here $N_{r}$, $N_{b}$ denote the batch sizes, $\lambda_{r}$, $\lambda_{b}$ are
the weights for the residual and boundary data respectively. Typically,
$\boldsymbol{x}_{r}^{i}$ and $\boldsymbol{x}_{b}^{j}$ are randomly sampled
from the uniform distributions $\mathcal{U}(\Omega)$ and
$\mathcal{U}(\partial\Omega)$ respectively, and all the derivatives involved
in the loss function can be computed accurately and quickly using reverse-mode
automatic differentiation [1, 16]. This ensures that information on
initial/boundary conditions propagates the interior domain while satisfying
the physics laws represented as PDEs, and is the main reason why we call it a
“physics-informed” neural network. If we have any sensor data for training, in
other words, the ground-truth data $(\boldsymbol{x}^{i}_{s},y^{i}_{s})$ with
$y^{i}_{s}=u(\boldsymbol{x}_{s}^{i})$ for all $i=1,2,\cdots,N_{s}$, we can
further add the so-called supervised loss $\mathcal{L}_{\rm
sup}=\sum^{N_{s}}_{i=1}|u(\boldsymbol{x}^{i}_{s};\boldsymbol{\theta})-y^{i}_{s}|^{2}$
to the total loss, which is known to improve the performance of the model.
However, in this paper, we will not consider this type of loss function, and
propose a novel method which can train the PINN efficiently without the use of
the training data. The schematic illustration of the PINN is depicted in
Figure 1.

###  2.2 Neural Tangent Kernel

In the seminar work [8], the authors proposed a novel approach, the so-called
neural tangent kernel (NTK), that opened up new avenues in the analysis of
training neural networks. The authors showed that under a certain parameter
initialization, as the width of the layer increases infinitely, the kernel
becomes constant during training so that training neural networks is
equivalent to performing a deterministic kernel regression. In other words,
under appropriate conditions, analysis of training using gradient descent can
be regarded as a kernel regression analysis. In this section, we will briefly
introduce the theory of NTK, which will be used in the later section to
theoretically investigate the method proposed in this paper.

Let us consider the scenario where we need to minimize the square loss
function on a dataset $(\boldsymbol{x}_{i},y_{i})^{m}_{i=1}$:
$\frac{1}{2}\sum^{m}_{i=1}(y_{i}-u(\boldsymbol{x}_{i};\boldsymbol{\theta}))^{2}$,
where $u$ is parametrized by $\boldsymbol{\theta}$. A continuous version of
gradient descent (gradient flow) can be written as

$\frac{{\rm d}\boldsymbol{\theta}(t)}{{\rm
d}t}=-\nabla_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta}(t))=-\sum_{i=1}^{m}(y_{i}-u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t)))\nabla_{\boldsymbol{\theta}}u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t)).$
(2.5)

Note further that

$\frac{{\rm d}u(\boldsymbol{x};\boldsymbol{\theta}(t))}{{\rm
d}t}=\bigg{[}{\frac{{\rm d}u(\boldsymbol{x};\boldsymbol{\theta}(t))}{{\rm
d}\boldsymbol{\theta}}}\bigg{]}\cdot\bigg{[}\frac{{\rm
d}\boldsymbol{\theta}(t)}{{\rm
d}t}\bigg{]}=\sum_{i=1}^{m}(y_{i}-u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t)))[\nabla_{\boldsymbol{\theta}}u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t))]^{T}[\nabla_{\boldsymbol{\theta}}u(\boldsymbol{x};\boldsymbol{\theta}(t))].$
(2.6)

This motivates us to define a component of NTK in the following way:

$\widehat{\boldsymbol{\Theta}}(\boldsymbol{x},\boldsymbol{x}^{\prime})(t)=\left\langle\frac{\partial
u(\boldsymbol{x};\boldsymbol{\theta}(t))}{\partial\boldsymbol{\theta}},\frac{\partial
u(\boldsymbol{x}^{\prime};\boldsymbol{\theta}(t))}{\partial\boldsymbol{\theta}}\right\rangle.$

Then over the dataset $(\boldsymbol{x}_{i},y_{i})^{m}_{i=1}$, (2.6) can be
rewritten as

$\frac{\rm d}{{\rm
d}t}{\left[\begin{matrix}u(\boldsymbol{x}_{1};\boldsymbol{\theta}(t))\\\
u(\boldsymbol{x}_{2};\boldsymbol{\theta}(t))\\\ \vdots\\\
u(\boldsymbol{x}_{m};\boldsymbol{\theta}(t))\end{matrix}\right]}=-\widehat{\boldsymbol{\Theta}}(t)\left[\begin{matrix}y_{1}-u(\boldsymbol{x}_{1};\boldsymbol{\theta}(t))\\\
y_{2}-u(\boldsymbol{x}_{2};\boldsymbol{\theta}(t))\\\ \vdots\\\
y_{3}-u(\boldsymbol{x}_{m};\boldsymbol{\theta}(t))\end{matrix}\right],$ (2.7)

where $\widehat{\boldsymbol{\Theta}}(t)\in\mathbb{R}^{m\times m}$ and the
$(i,j)$-th component of $\widehat{\boldsymbol{\Theta}}(t)$ is

$\widehat{\boldsymbol{\Theta}}_{ij}(t)=\widehat{\boldsymbol{\Theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=[\nabla_{\boldsymbol{\theta}}u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t))]^{T}[\nabla_{\boldsymbol{\theta}}u(\boldsymbol{x}_{j};\boldsymbol{\theta}(t))].$

If we write $\boldsymbol{U}(t)=\\{U_{i}(t)\\}^{m}_{i=1}$ with
$U_{i}(t)=y_{i}-u(\boldsymbol{x}_{i};\boldsymbol{\theta}(t))$, and suppose the
NTK remains constant and positive definite during training (i.e.
$\widehat{\boldsymbol{\Theta}}(t)=\widehat{\boldsymbol{\Theta}}(0)$ and all
the eigenvalues of $\widehat{\boldsymbol{\Theta}}(0)$ are positive). Then from
(2.7), we have

$\frac{\rm d}{{\rm
d}t}\left(\frac{1}{2}\|\boldsymbol{U}(t)\|^{2}_{2}\right)=-\boldsymbol{U}^{T}(t)\widehat{\boldsymbol{\Theta}}(0)\boldsymbol{U}(t)\leq-\frac{\lambda_{0}}{2}\|\boldsymbol{U}(t)\|^{2}_{2},$

where $\lambda_{0}$ is the smallest eigenvalue of
$\widehat{\boldsymbol{\Theta}}(0)$, which leads us to the convergence result

$\|\boldsymbol{U}(t)\|^{2}_{2}\leq
e^{-\lambda_{0}t}\|\boldsymbol{U}(0)\|^{2}_{2}.$ (2.8)

Next, we shall discuss the assumption that the NTK remains constant during the
training process, and we will briefly review the relevant theory. To do this,
let us consider an $L$-layer fully-connected neural network
$f(\boldsymbol{x})=g_{L}(\boldsymbol{x})$ defined recursively as

$g_{\ell}(\boldsymbol{x})=\frac{1}{\sqrt{n_{\ell-1}}}W_{\ell}\boldsymbol{x}_{\ell-1}(\boldsymbol{x}),\quad\boldsymbol{x}_{\ell-1}(\boldsymbol{x})=\sigma(g_{\ell-1}(\boldsymbol{x})),\quad\boldsymbol{x}_{0}(x)=\boldsymbol{x}$

for $\ell=1,2,\cdots,L$ where $n_{\ell}$ and
$W^{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}$ are the number of neurons
and the weight matrix in $L$-the layer respectively. We suppose that all
weights are initialized with standard Gaussian distributions independently.
Let us further assume the circumstance that we aim to minimize a generic
differentiable loss function $\mathcal{L}$:

$\frac{{\rm d}\boldsymbol{\theta}(t)}{{\rm
d}t}=-\nabla_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta}(t))=-\sum_{i=1}^{m}\frac{\partial\mathcal{L}(y_{i},z)}{\partial
z}\bigg{|}_{z=f(\boldsymbol{x}_{i};\boldsymbol{\theta}(t))}\nabla_{\boldsymbol{\theta}}f(\boldsymbol{x}_{i};\boldsymbol{\theta}(t)),$

where $\boldsymbol{\theta}$ is a collection of all weights $\\{W_{\ell}\\}$.
One of the main results of the paper [8] is encapsulated in the following
theorem.

######  Theorem 2.1.

Assume that the activation function $\sigma$ is of class $C^{2}$ and
Lipschitz, and the loss function $\mathcal{L}$ is Lipschitz. Then the neural
tangent kernel $\widehat{\boldsymbol{\Theta}}(t)$ converges in probability to
a deterministic kernel independent of $t>0$ as $n_{1}$, $n_{2},\cdots
n_{L-1}\rightarrow\infty.$

The above theorem demonstrates that training a neural network under Gaussian
parameterization is equivalent to conducting a certain kernel method as widths
go to infinity.

###  2.3 Neural Tangent Kernel for PINNs

In [28], the authors investigated the NTK theory for PINNs. They have shown
that when the parameters are initialized by the Gaussian distribution and the
number of neurons in each layer goes to infinity, the NTK for PINNs is
determined by the initial state and it remains constant during the training.
To illustrate the relevant theory in more detail, let us consider the problem
(2.1), and the corresponding physics-informed loss function (2.2) with
(2.3)-(2.4). Then it was proved in [28] that the evolution of $u$ and
$\mathcal{D}[u]$ can be described by the gradient flow:

$\left[\begin{matrix}\frac{{\rm{d}}u(\boldsymbol{x}_{b};\boldsymbol{\theta}(t))}{{\rm{d}}t}\\\
\frac{{\rm{d}}\mathcal{D}[u](\boldsymbol{x}_{r};\boldsymbol{\theta}(t))}{{\rm{d}}t}\end{matrix}\right]=-\left[\begin{matrix}\boldsymbol{K}_{uu}(t)&\boldsymbol{K}_{ur}(t)\\\
\boldsymbol{K}_{ru}(t)&\boldsymbol{K}_{rr}(t)\end{matrix}\right]\left[\begin{matrix}u(\boldsymbol{x}_{b};\boldsymbol{\theta}(t))-g(\boldsymbol{x}_{b})\\\
\mathcal{D}[u](\boldsymbol{x}_{r};\boldsymbol{\theta}(t))-f(\boldsymbol{x}_{r})\end{matrix}\right],$

where $\boldsymbol{K}_{uu}(t)\in\mathbb{R}^{N_{b}\times N_{b}}$,
$\boldsymbol{K}_{rr}\in\mathbb{R}^{N_{r}\times N_{r}}$ and
$\boldsymbol{K}_{ru}(t)=[\boldsymbol{K}_{ur}(t)]^{T}\in\mathbb{R}^{N_{r}\times
N_{b}}$ whose $(i,j)$-th entries are

$\displaystyle(\boldsymbol{K}_{uu})_{ij}(t)=\left<\frac{{\rm{d}}u(\boldsymbol{x}_{b}^{i};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}},\frac{{\rm{d}}u(\boldsymbol{x}_{b}^{j};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}}\right>,$
(2.9)
$\displaystyle(\boldsymbol{K}_{ru})_{ij}(t)=\left<\frac{{\rm{d}}\mathcal{D}[u](\boldsymbol{x}_{r}^{i};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}},\frac{{\rm{d}}u(\boldsymbol{x}_{b}^{j};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}}\right>,$
$\displaystyle(\boldsymbol{K}_{rr})_{ij}(t)=\left<\frac{{\rm{d}}\mathcal{D}[u](\boldsymbol{x}_{r}^{i};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}},\frac{{\rm{d}}\mathcal{D}[u](\boldsymbol{x}_{r}^{j};\boldsymbol{\theta}(t))}{{\rm{d}}\boldsymbol{\theta}}\right>.$

The matrix

$\boldsymbol{K}(t)=\left[\begin{matrix}\boldsymbol{K}_{uu}(t)&\boldsymbol{K}_{ur}(t)\\\
\boldsymbol{K}_{ru}(t)&\boldsymbol{K}_{rr}(t)\end{matrix}\right]$

is called the NTK of a PINN. Motivated by the fact that an infinitely-wide
neural network becomes a Gaussian process and its NTK is constant during the
training [8], [28] showed that, if the number of neurons of each layer goes to
infinity, the NTK of a PINN $\boldsymbol{K}(t)$ indeed converges in
probability to a deterministic kernel $\boldsymbol{K}^{*}$ and
$\sup_{t\in[0,T]}\|\boldsymbol{K}(t)-\boldsymbol{K}(0)\|_{2}\rightarrow 0$ for
any $T>0$, under some appropriate assumptions.

Since the kernel $\boldsymbol{K}(t)$ is positive semi-definite, by the
spectral decomposition, we can see that the eigenvalues of the kernel
$\boldsymbol{K}(t)$ provide some information on the convergence rate for
training as presented in (2.8) for given data points (see [28] for details).
However, as the number of training data increases, it is not easy to
explicitly compute all the eigenvalues individually. Therefore, based on the
fact that the trace is equal to the sum of eigenvalues, motivated by [28], we
shall used the following averaged quantity to measure the approximate rate of
convergence for the training of PINNs

$\frac{{\rm Tr}({\boldsymbol{K}_{uu}}(t))}{N_{b}}+\frac{{\rm
Tr}({\boldsymbol{K}_{rr}}(t))}{N_{r}}.$ (2.10)

In the later section, we will explicitly compute the above quantity when the
variables are scaled to illustrate the faster convergence of the proposed
method in this paper.

##  3 Variable-Scaling Physics-Informed Neural Networks

In this section, we will describe the method we mainly propose in this paper.
The key idea of our method is to use a variable-scaling technique in the
training process. More precisely, we exploit the change of variable
$\boldsymbol{x}\mapsto\overline{\boldsymbol{x}}/N$ to define a loss function,
in order to relieve the potential stiff behavior of solutions. After the
training is completed, we scale back to the original variable for the solution
prediction.

To illustrate the methodology in more detail, let us consider the following
PDE of the form:

$\displaystyle F(\boldsymbol{x};u,Du,D^{2}u,\cdots,D^{k}u)$ $\displaystyle=0$
$\displaystyle\boldsymbol{x}\in\Omega,$ (3.1) $\displaystyle
u(\boldsymbol{x})$ $\displaystyle=g(\boldsymbol{x})$
$\displaystyle\boldsymbol{x}\in\partial\Omega,$

where $F(\boldsymbol{x})=\mathcal{D}[u](\boldsymbol{x})-f(\boldsymbol{x})$ in
(2.1) and $D^{j}$ denotes the mixed partial derivatives of degree
$j\in\mathbb{N}$. Now the corresponding residual loss function for the PINN is
then given by

$\mathcal{L}_{\rm
original}=\frac{\lambda_{r}}{N_{r}}\sum^{N_{r}}_{i=1}\bigg{|}F(\boldsymbol{x}^{i}_{r};u(\boldsymbol{x}^{i}_{r}),Du(\boldsymbol{x}^{i}_{r}),D^{2}u(\boldsymbol{x}^{i}_{r}),\cdots,D^{k}u(\boldsymbol{x}^{i}_{r}))\bigg{|}^{2}+\frac{\lambda_{b}}{N_{b}}\sum^{N_{b}}_{j=1}\left|u(\boldsymbol{x}^{j}_{b})-g(\boldsymbol{x}^{j}_{b})\right|^{2},$
(3.2)

where the $\boldsymbol{x}^{i}_{r}$ and $\boldsymbol{x}^{j}_{b}$ are i.i.d.
random samples chosen from the uniform distribution
$\boldsymbol{x}^{i}_{r}\sim\mathcal{U}(\Omega)$ and
$\boldsymbol{x}^{j}_{b}\sim\mathcal{U}(\partial\Omega)$ respectively. Next,
let us introduce a new variable $\overline{\boldsymbol{x}}=N\boldsymbol{x}\in
N\Omega:=\\{N\boldsymbol{x}:\boldsymbol{x}\in\Omega\\}$ with the notation
$u(\overline{\boldsymbol{x}}/N)=v(\overline{\boldsymbol{x}})$. Intuitively, as
we can see from Figure 2, the change of variable
$\boldsymbol{x}=\overline{\boldsymbol{x}}/N$ can be interpreted as an
$N$-times “zoom-in” for the solution profile, and thus it has an effect of
reducing any sharp transitions of solution dynamics on the designated domain.

Figure 2: An illustration of the effect of the variable scaling when the
scaling factor $N=50$ is applied to the function $\exp(-10000x^{2})$. The
scaling has the effect of “zoom-in” and the resulting function exhibits
relatively less stiff behavior on the designated domain.

By the chain rule, we know that

$D_{\boldsymbol{x}}u(\boldsymbol{x})=D_{\overline{\boldsymbol{x}}}u\left(\frac{\overline{\boldsymbol{x}}}{N}\right)\frac{{\rm
d}\overline{\boldsymbol{x}}}{{\rm
d}\boldsymbol{x}}=ND_{\overline{\boldsymbol{x}}}v(\overline{\boldsymbol{x}}),$

and this can be easily extended to
$D^{k}_{\overline{\boldsymbol{x}}}u(\boldsymbol{x})=N^{k}D_{\overline{\boldsymbol{x}}}v(\overline{\boldsymbol{x}})$.
Therefore, the given PDE (3.1) can be rewritten as

$\displaystyle
F(\overline{\boldsymbol{x}}/N;v,NDv,N^{2}D^{2}v,\cdots,N^{k}D^{k}u)$
$\displaystyle=0$ $\displaystyle\overline{\boldsymbol{x}}\in N\Omega,$ (3.3)
$\displaystyle v(\overline{\boldsymbol{x}})$
$\displaystyle=g(\overline{\boldsymbol{x}}/N)$ $\displaystyle\boldsymbol{x}\in
N\partial\Omega.$

The associate loss function for the representation (3.3) is

$\mathcal{L}_{\rm
scaled}=\frac{\lambda_{r}}{N_{r}}\sum^{N_{r}}_{i=1}\bigg{|}F(\overline{\boldsymbol{x}}^{i}_{r}/N;v(\overline{\boldsymbol{x}}^{i}_{r}),NDv(\overline{\boldsymbol{x}}^{i}_{r}),N^{2}D^{2}v(\overline{\boldsymbol{x}}^{i}_{r}),\cdots,N^{k}D^{k}v(\overline{\boldsymbol{x}}^{i}_{r}))\bigg{|}^{2}+\frac{\lambda_{b}}{N_{b}}\sum^{N_{b}}_{j=1}\bigg{|}v(\overline{\boldsymbol{x}}^{j}_{b})-g\bigg{(}\frac{\overline{\boldsymbol{x}}^{j}_{b}}{N}\bigg{)}\bigg{|}^{2},$
(3.4)

with $\overline{\boldsymbol{x}}^{i}_{r}\sim\mathcal{U}(N\Omega)$ and
$\overline{\boldsymbol{x}}^{j}_{b}\sim\mathcal{U}(N\partial\Omega)$. In this
scaled loss function, the weight parameters $\lambda_{r}$ and $\lambda_{b}$
should be carefully chosen according to the result of the scaling. For
example, they can be chosen to cancel a common factor or extreme coefficients
(see Section 4 for details). The key novel part of our method is to minimize
the loss function (3.4) instead of (3.2). Once the training is finished and we
obtain the solution prediction $v(\overline{\boldsymbol{x}})$, we pull the
variable back to the original scale by the formula

$u(\boldsymbol{x})=v(N\boldsymbol{x}),$ (3.5)

which is the original solution we aimed to predict at the initial stage. Since
the domain $\Omega$ has become larger to $N\Omega$ in the method, one may
wonder whether we need more random samples for the training. However, as we
will illustrate in the numerical experiments in the next section, this is not
the case, and we can significantly improve the training efficiency even if we
train the model with the same number of random samples from the scaled domain.

##  4 Numerical Experiments

In this section, we will perform a series of numerical experiments to evaluate
the performance of the proposed training method against the standard PINNs. As
mentioned earlier, the choice of parameters $\lambda_{r}$ and $\lambda_{b}$ in
(3.4) as well as the scaling factor $N$ is important for the overall
performance and they should be carefully selected. As we will see in the
numerial experiments in this section, as $N$ increases, the parameters related
to training become more sensitive, which tends to make the training somewhat
unstable with more fluctuation in learning curves. Therefore, to leverage the
advantages of the VS-PINN while ensuring adequate training stability, it is
necessary to choose an appropriate $N$. We will also discuss this issue in
more details in the later section of this paper, and here we shall demonstrate
some basic strategies for choosing these parameters in the experiments.

Throughout the section, the model trained with the original PINN will be
compared with the model trained by the proposed variable-scaling technique,
and all the other settings including the model architecture, initialization,
and the optimizer will be the same. And most importantly, even though the
domain is $N$-times extended, the number of collocation points for the
training is the same for both the original method and the variable-scaling
method.

###  4.1 Wave Equation

As a starting point, we shall conduct some experiments solving the one-
dimensional wave equation, where the initial data with a high-frequency part
is imposed. Namely, we will solve the following problem:

$\displaystyle u_{tt}-u_{xx}=0$ $\displaystyle(x,t)\in(0,1)\times(0,1),$ (4.1)
$\displaystyle u(0,t)=u(1,t)=0$ $\displaystyle t\in[0,1],$ $\displaystyle
u(x,0)=0$ $\displaystyle x\in[0,1],$ $\displaystyle u_{t}(x,0)=2\pi\sin(2\pi
x)+10\pi\sin(10\pi x)$ $\displaystyle x\in[0,1].$

The exact solution for the above problem can be represented as

$u(x,t)=\sin(2\pi x)\sin(2\pi t)+\sin(10\pi x)\sin(10\pi t),$

where the $\sin(10\pi x)\sin(10\pi t)$ part exhibit a high frequency.

For the baseline model, we use a neural network with 4 hidden layers with 128
neurons for each, activated by the hyperbolic tangent function. By scaling the
variables to expand the domain, we can relatively mitigate the high-frequency
regions, which could potentially improve the training process. To demonstrate
this, we conducted experiments with three different scaling factors
$N=1,4,10$, which leads us to compute the solutions in three different
computational domains $[0,1]^{2},[0,4]^{2}$ and $[0,10]^{2}$. The scaled
problem of (4.1) can be written as

$\displaystyle N^{2}u_{tt}-N^{2}u_{xx}=0$
$\displaystyle(x,t)\in(0,N)\times(0,N),$ (4.2) $\displaystyle u(0,t)=u(N,t)=0$
$\displaystyle t\in[0,N],$ $\displaystyle u(x,0)=0$ $\displaystyle x\in[0,N],$
$\displaystyle u_{t}(x,0)=\frac{2\pi}{N}\sin\left(\frac{2\pi
x}{N}\right)+\frac{10\pi}{N}\sin\left(\frac{10\pi x}{N}\right)$ $\displaystyle
x\in[0,N],$

and the corresponding scaled loss is

$\mathcal{L}_{\rm{scaled}}=2\mathcal{L}_{\rm
data}+\frac{1}{N^{4}}\mathcal{L}_{\rm res}.$ (4.3)

As mentioned earlier, note that it is reasonable to choose
$\lambda_{r}=\frac{1}{N^{4}}$ to cancel the $N^{4}$ factor in the residual
loss function. For the training, we randomly sample 6400 collocation points
for $\mathcal{L}_{\rm res}$ and 512 points for $\mathcal{L}_{\rm data}$. The
learning-curve comparison for the methods with different scales $N=1$
(standard PINN) and $N=4,10$ (VS-PINN) is presented in Figure 3. For each
epoch, we compute the relative $L^{2}$ error between the true solution and the
PINN to evaluate the training efficiency. As the scale increases, we can
observe in the figure that the error decreases more quickly and rapidly,
meaning that the VS-PINN can learn the target solution with high frequency
more efficiently.

Figure 3: Learning curves for the VS-PINNs with the scales $N=1,4,10$.

We train the model in a zoom-in state, followed by an evaluation in a zoom-out
state with the process: $u(x_{\rm test},t_{\rm test})=u_{\rm model}(Nx_{\rm
test},Nt_{\rm test})$. After the 60000 epochs of training, a comparison
between the true solution and the predicted solution is illustrated in Figure
4. As we can see from Figure 4, the solution prediction with variable scaling
is much more accurate for the high-frequency solution compared to the solution
computed by the standard PINN.

Figure 4: One-dimensional wave equation: (a) The exact solution versus the
predicted solution by training the standard PINN. The relative $L^{2}$ error
is 6.31e-1. (b) The exact solution versus the predicted solution by training
the VS-PINN with the scale $N=10$. The relative $L^{2}$ error is 1.20e-2.

###  4.2 Allen–Cahn Equation

We shall also apply the variable-scaling technique to the PINNs for other
types of equations. In particular, it is noteworthy that the change of
variable can be used to offset the effect of the extreme coefficient that
appeared in the given equation. This idea will be examined through the
numerical experiments for the Allen–Cahn equation in this subsection and the
singular perturbation problem in the next subsection in order.

The Allen–Cahn equation arises in mathematical physics, which describes the
process of phase separation in multi-component alloy systems. From a
mathematical perspective, it is the reaction-diffusion-type equation with
nonlinear terms, typically assumed to have a small diffusion parameter. We
shall perform the experiment in the same setting adopted in [18], which is
represented by the following formulation:

$\displaystyle u_{t}-0.0001u_{xx}+5u^{3}-5u=0$ $\displaystyle
x\in(-1,1),\;\;\;t\in(0,1),$ (4.4) $\displaystyle u(x,0)=x^{2}\cos(\pi x)$
$\displaystyle x\in[-1,1],$ $\displaystyle u(-1,t)=u(1,t)$ $\displaystyle
t\in[0,1],$ $\displaystyle u_{x}(-1,t)=u_{x}(1,t)$ $\displaystyle t\in[0,1].$

For the model architecture, we exploit the neural network consisting of 4
hidden layers with 64 neurons for each, and the corresponding parameters are
initialized with random initialization using a fixed seed. For training, we
first proceed with supervised learning and train the neural network with pre-
computed training data as performed in [18]. In addition, in order to further
highlight the effectiveness of our method, we will also train the VS-PINN in
an unsupervised manner and confirm that we can obtain satisfactory results
even for unsupervised learning.

For the supervised learning part, data points that are randomly chosen from
the exact solution at $t=0.1$ as conducted in [18]. For the other collocation
points (both for supervised and unsupervised learning) we select $128$ random
samples for the periodic boundary conditions, $200$ points for the initial
condition and $10000$ data points for the interior. For the scaling, we
consider the change of variable only for the spatial variable
$x\mapsto\overline{x}/N$. Through a derivative, we can cancel the diffusion
coefficient which might cause a stiff behavior. More precisely, to define the
scaled loss function, we shall consider the following scaled problem:

$\displaystyle u_{t}-N^{2}\cdot 0.0001u_{xx}+5u^{3}-5u=0$ $\displaystyle
x\in(-N,N),\;\;\;t\in(0,1),$ (4.5) $\displaystyle u(x,0)=(x/N)^{2}\cos(\pi
x/N)$ $\displaystyle x\in[-N,N],$ $\displaystyle u(-N,t)=u(N,t)$
$\displaystyle t\in[0,1],$ $\displaystyle u_{x}(-N,t)=u_{x}(N,t)$
$\displaystyle t\in[0,1].$

As we can see from the scaled problem (4.5), it would be reasonable to choose
some large scaling factor $N$ so that we can cancel the diffusion coefficient
$0.0001$, which mitigates the potential extreme behavior of solution caused by
the biased coefficient. In the simulation, we shall examine the cases when $N$
is $1$ (the standard PINN), $10$ and $100$ (the VS-PINN). Finally, we use the
training loss function with adjusted parameters $\lambda_{\rm data}=2$ and
$\lambda_{\rm res}=0.3$:

$\mathcal{L}_{\rm scaled}=2\mathcal{L}_{\rm data}+0.3\mathcal{L}_{\rm res}.$

As before, the model is trained in a zoom-in state and the performance is
evaluated in a zoom-out state which is the original scale: $u(x_{\rm
test},t_{\rm test})=u_{\rm model}(Nx_{\rm test},t_{\rm test})$.

(a) Supervised learning

(b) Unsupervised learning

Figure 5: Traning dynamics for the Allen–Cahn equation using VS-PINN with the
scales $N=1$ (standard PINN), $10$, $100$. (a) is for the model trained with
the pre-computed training data at $t=0.1$ and (b) is for the model trained in
an unsupervised manner.

Figure 5 displays the relative $L^{2}$ error for the standard PINN and the VS-
PINN as a function of the training epoch. From the training curves, we can
confirm that we can significantly improve the training efficiency and obtain
higher accuracy both for supervised and unsupervised scenarios. In the case of
unsupervised learning, as presented in Figure 5 (b), it was observed that as
the domain expanded, the point at which the error began to decrease was
delayed when $N=100$. However, as the training proceeds further, we can
finally obtain a faster convergence with higher accuracy. For both cases, the
solution profiles trained in an unsupervised manner are presented in Figure 6.
In particular, the snapshots at $t=0.5$ and $t=0.9$ are illustrated in the
bottom line for both cases. The error was observed to be 2.05e-2 for the
standard PINN ($N=1$) and 3.36e-3 for the VS-PINN ($N=100$). As we can see
from these pictures, the solution prediction via the VS-PINN better captures
the stiff behavior of the solution.

Figure 6: Top left: Plot of the prediction $\widehat{u}(x,t)$ by the standard
PINN. Top right: Plot of the prediction $\widehat{u}(x,t)$ via the VS-PINN.
Bottom left: Snapshots of the prediction $\widehat{u}(x,t)$ using the original
PINN (blue solid curve) vs. the exact solution $u(x,t)$ (red dashed curve) at
$t=0.5$ and $t=0.9$. Bottom right: Snapshots of the prediction
$\widehat{u}(x,t)$ via the VS-PINN (blue solid curve) vs. the exact solution
$u(x,t)$ (red dashed curve) at $t=0.5$ and $t=0.9$.

###  4.3 Boundary-Layer Problem

The next example focuses on the singularly perturbed boundary-layer problem,
which is a well-known challenge in scientific computing. If the diffusion
coefficient is small, the solution exhibits a steep increase in thin regions
near boundaries which may lead to failures in overall prediction. Typically,
in such a case, neural-network-based solution prediction often has
difficulties in finding optimal solutions. This is also mainly due to the
spectral bias phenomenon mentioned in the introduction of this paper [17]. We
aimed to alleviate such sharp changes through a variable scaling.
Specifically, as mentioned earlier, if we make a change of variable, a small
diffusion coefficient can be canceled by the Jacobian factor. To illustrate
this idea, we shall consider the following one-dimensional singular
perturbation problem: for a small parameter $0<\varepsilon\ll 1$,

$\displaystyle\varepsilon u_{xx}+u_{x}+1=0\;\;\;x\in[0,1],$ (4.6)
$\displaystyle u(0)=u(1)=0.$

The exact solution can be computed analytically, and is presented as follows:

$u(x)=\frac{1}{1-e^{-1/\varepsilon}}-x-\frac{1}{1-e^{-1/\varepsilon}}e^{-x/\varepsilon}.$

When $\varepsilon$ is sufficiently small, we can ignore the term
$e^{-1/\varepsilon}$ and the solution is reduced to

$u(x)=1-x-e^{-x/\varepsilon}.$

This solution has a slope of $1/\varepsilon-1$ at $x=0$ and it attains its
maximum at $x=-\varepsilon\ln(\varepsilon)$, indicating there is an explosive
change near the boundary. We set $\varepsilon=$1e-6 and perform the change of
variable with the scale $N=1000$ so that we can remove the effect of
$\varepsilon>0$ by multiplying the Jacobian factor. Then, the original
equation is transformed as follows:

$\displaystyle u_{xx}+1000u_{x}+1=0\;\;\;x\in[0,1000],$ (4.7) $\displaystyle
u(0)=u(1000)=0,$

and we set the loss as

$\mathcal{L}_{\rm scaled}=20\mathcal{L}_{\rm data}+\mathcal{L}_{\rm res}.$

Figure 7: Left: Exact solution of (4.6) compared with the solution prediction
by the standard PINN and the VS-PINN. Right: Learning curves for both the
original PINN and the VS-PINN.

We shall use the model with 8 hidden layers, with 20 neurons per layer and the
model is initialized with a fixed seed and is trained using 1000 interior
collocation data. When the model is trained with a variable scaling, the
highest accuracy achieved was 1.17e-2, while in the original case, there is
almost no progress in training and the learning curve remains mostly
unchanged. Therefore, we can conclude that the mitigation of stiff behavior
through the variable scaling has a significant impact on training and can make
a much more accurate solution prediction for the boundary-layer problem.

###  4.4 The Navier–Stokes Equations

Another advantage of using the VS-PINN is that the neural network can learn
the nonlinear structure more efficiently by alleviating the potentially
complicated dynamics of solutions. As a last experiment, we will discuss the
experimental results of the challenging Navier–Stokes equations. The numerical
computation of solutions to fluid equations has a long and rich history, and
it has seen significant advancements across various fields. Furthermore, in
recent years, there has been a wide range of studies to predict fluid flow
based on machine-learning methods, and this topic has become of independent
interest. For example [21] proposed a mixed-variable PINN to compute the
solution of the incompressible Navier–Stokes equations, [4] proposed the
enhanced PINN to solve the three-dimensional Navier–Stokes equations and [9]
performed PINN computations by considering both the velocity-pressure
formulation and the vorticity-velocity to simulate the flows from laminar to
turbulent formulation.

As an extension of these prior research papers, in this paper, we shall
conduct the experiments for the steady and transient laminar flows passing
over a circular cylinder in the two-dimensional domain without any training
data as investigated in [21]. More precisely, we shall consider the following
incompressible Newtonian fluid flow problem, described by the following
Navier–Stokes equations:

$\displaystyle u_{x}+v_{y}$ $\displaystyle=0$ $\displaystyle(x,y)\in\Omega,$
(4.8) $\displaystyle uu_{x}+vu_{y}$
$\displaystyle=-\frac{1}{\rho}p_{x}+\frac{\mu}{\rho}(u_{xx}+u_{yy})$
$\displaystyle(x,y)\in\Omega,$ $\displaystyle uv_{x}+vv_{y}$
$\displaystyle=-\frac{1}{\rho}p_{y}+\frac{\mu}{\rho}(v_{xx}+v_{yy})$
$\displaystyle(x,y)\in\Omega,$

where $(u,v)\in\mathbb{R}^{2}$ is the velocity vector, $p\in\mathbb{R}$ is the
pressure and the viscosity and the density of the fluid are set to be
$\mu=0.02$ and $\rho=1$ respectively. The computational domain $\Omega$ is
depicted in Figure 8 which will be expanded by a factor of $N$ in the training
for the VS-PINN.

Figure 8: An illustration of the computational domain where we compute the
solution of the incompressible Navier–Stokes equations.

A velocity profile is applied on the inlet while the condition of zero
pressure is imposed on the outlet, as shown in Figure 8. No-slip boundary
conditions are assumed on the wall and cylinder boundaries. Namely, we enforce
the following conditions to the equations (4.8):

$\displaystyle u(x,y)=v(x,y)=0$ $\displaystyle\text{on wall and cylinder
boundaries},$ (4.9) $\displaystyle p(1.1,y)=0$ $\displaystyle y\in[0,0.41],$
$\displaystyle u(0,y)=4(0.41-y)\frac{y}{0.41^{2}}$ $\displaystyle
y\in[0,0.41],$ $\displaystyle v(0,y)=0$ $\displaystyle y\in[0,0.41].$

As before, the corresponding scaled boundary-value problem can be formulated
as

$\displaystyle Nu_{x}+Nv_{y}$ $\displaystyle=0$ $\displaystyle(x,y)\in
N\Omega,$ (4.10) $\displaystyle uNu_{x}+vNu_{y}$
$\displaystyle=-\frac{1}{\rho}Np_{x}+\frac{\mu}{\rho}(N^{2}u_{xx}+N^{2}u_{yy})$
$\displaystyle(x,y)\in N\Omega,$ $\displaystyle uNv_{x}+vNv_{y}$
$\displaystyle=-\frac{1}{\rho}Np_{y}+\frac{\mu}{\rho}(N^{2}v_{xx}+N^{2}v_{yy})$
$\displaystyle(x,y)\in N\Omega,$ $\displaystyle u(x,y)$
$\displaystyle=v(x,y)=0$ on wall and cylinder boundaries, $\displaystyle
p(N\cdot 1.1,y)$ $\displaystyle=0$ $\displaystyle y\in[0,N\times 0.41],$
$\displaystyle u(0,y)$ $\displaystyle=4(0.41-y/N)\frac{y/N}{0.41^{2}}$
$\displaystyle y\in[0,N\times 0.41],$ $\displaystyle v(0,y)$ $\displaystyle=0$
$\displaystyle y\in[0,N\times 0.41],$

and we define the associated scaled loss function by

$\mathcal{L}=\frac{1}{N^{2}}\mathcal{L}_{\rm res}+2\mathcal{L}_{\rm data}.$

Figure 9: The standard PINN vs. the VS-PINN trained by the Adam+L-BFGS
optimizer.

We train the models both for the standard PINN and the VS-PINN with two
different optimizing strategies to highlight the efficiency of the proposed
training method. For the comparison, we first exploit the strategy used in
[21], that is, Adam optimizer with first $10000$ epochs and followed by the
extra $10000$ epochs using the L-BFGS. For the Adam optimizer, we randomly
sampled $1200$ collocation points from the boundary and $6000+600$ (near the
cylinder) interior data points. For the L-BFGS optimizer, 1200 random samples
were chosen from the boundary and $8000+1000$ (near the cylinder) data points
were used for the interior. The Adam optimizer randomly samples the points
from the interior at each epoch, while L-BFGS samples are fixed throughout
training. The test errors are 1.83e-2 for the velocity and 2.37e-2 for the
pressure, and the training curves are presented in Figure 9.

Figure 10: The standard PINN vs. the VS-PINN trained by the Adam optimizer
only with a smaller number of collocation points.

Note that even though we utilized the smallest model (5 hidden layers with 40
neurons for each) among the ones used in [21], the achieved accuracy is
comparable to the accuracy obtained by using the larger model (7 hidden layers
with 40 neurons each) in [21], which trained on $50000$ data points. Hence, we
can confirm that the VS-PINN is effective in computing the approximate
solutions for the Navier–Stokes equations.

On the other hand, the second strategy only utilizes the Adam optimizer with
$1200$ collocation points on the boundary and $6000+600$ (near the cylinder)
data points in the interior and proceeds $40000$ iteration. The final test
errors are 2.18e-2 for the velocity and 5.75e-2 for the pressure. It is
noteworthy that even though we only use much smaller number of training
samples compared with the previous option using the L-BFGS as done in [21], we
can still obtain the result with high accuracy, which was not observed in
other previous works.

The exact solutions and the predicted velocity and pressure computed by both
the standard PINN and the VS-PINN are presented in Figure 11. From these
contour plots, we can observe that the VS-PINN makes a much more accurate
solution compared with the solution predicted by the original PINN. Moreover,
since our method is universal, we can combine the proposed VS-PINN method with
the previous ones to obtain better performance. For example, we can use
various formulations proposed in the previous works [9, 4, 21] and apply the
variable-scaling technique to these additional formulations.

Figure 11: Contour plots for the solution and the approximate solutions. Left
column: Exact solutions for the problem (4.8)-(4.9). Middle column: Solution
predictions using the standard PINN ($N=1$). Right column: Solution
predictions using the VS-PINN ($N=10$).

##  5 Neural Tangent Kernel Analysis

In this section, we shall investigate the effect of the variable scaling on
training from an NTK perspective. Let us first describe the relevant problem
that we shall consider for the analysis. As we discussed in Section 2.3, the
averaged convergence rate can be expressed in terms of the quantities ${\rm
Tr}({\boldsymbol{K}_{uu}(0)})$ and ${\rm Tr}({\boldsymbol{K}_{rr}(0)})$. In
order to convey the idea clearly, we will focus on a simple computational
setting and explicitly compute the values ${\rm Tr}({\boldsymbol{K}_{uu}(0)})$
and ${\rm Tr}({\boldsymbol{K}_{rr}(0)})$. More precisely, we will consider the
one-dimensional Poisson problem with a homogeneous Dirichlet boundary
condition

$\displaystyle-u_{xx}$ $\displaystyle=f(x)$ $\displaystyle x\in[0,1],$ (5.1)
$\displaystyle u(0)$ $\displaystyle=u(1)=0,$

whose scaled problem can be written as

$\displaystyle-u_{xx}$ $\displaystyle=f(x/N)/{N^{2}}\;\;$ $\displaystyle
x\in[0,N],$ (5.2) $\displaystyle u(0)$ $\displaystyle=u(N)=0.$

We then define the corresponding scaled loss as

$\mathcal{L}_{\rm
scaled}=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}|u_{xx}(x_{r}^{i};\boldsymbol{\theta})+f(x_{r}^{i}/N)/{N^{2}}|^{2}+\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}|u(x_{b}^{j};\boldsymbol{\theta})|^{2}.$

Regarding the model architecture, for the sake of simplicity, we consider a
neural network $u(x;\boldsymbol{\theta})$ with one hidden layer, activated by
the function $\sigma=\max\\{0,x^{3}\\}$, which was considered in some previous
literature (see, e.g., [31]):

$u(x;\boldsymbol{\theta})=\frac{1}{\sqrt{d_{1}}}\sum_{k=1}^{d_{1}}W_{2}^{k}\sigma(W_{1}^{k}x+b_{1}^{k})+b_{2},$
(5.3)

where $d_{1}$ denotes the number of neurons in the hidden layer, and the
parameters $\\{W^{k}_{1}\\}_{k=1}^{d_{1}}$, $\\{W^{k}_{2}\\}_{k=1}^{d_{1}}$,
$\\{b^{k}_{1}\\}_{k=1}^{d_{1}}$ and $b_{2}$ are i.i.d. randomly chosen from
the Gaussian distribution $\mathcal{N}(0,1)$. Before proceeding further, we
introduce the following lemma which will be used to compute the averaged
eigenvalues (see, e.g., [7]).

######  Lemma 5.1.

If $X\sim\mathcal{N}(0,\delta^{2})$, then there holds

$\mathbb{E}\left[X^{2n}\right]=\frac{(2n)!}{2^{n}n!}\delta^{2n}\quad{\rm{and}}\quad\mathbb{E}\left[X^{2n-1}\right]=0\qquad\forall
n\in\mathbb{N}.$

###  5.1 Averaged Convergence Rate

In this section, we will explicitly compute the quantities ${\rm
Tr}(\boldsymbol{K}_{uu}(0))/{N_{b}}$ and ${\rm
Tr}(\boldsymbol{K}_{rr}(0))/{N_{r}}$ for the scaled problem (5.2) on the given
collocation points $\\{x_{b}^{j}\\}_{j=1}^{N_{b}}$ and
$\\{x_{r}^{i}\\}_{i=1}^{N_{r}}$ with
$x_{b}^{j}\sim\,{\rm{i.i.d.}}\,\,\mathcal{U}(\partial[0,N])$ and
$x_{r}^{i}\sim\,{\rm{i.i.d.}}\,\,\mathcal{U}((0,N))$. We first compute the
convergence rate regarding the boundary collocation points. From the
definition (2.9), we see that

$\frac{{\rm
Tr}(\boldsymbol{K}_{uu}(0))}{N_{b}}=\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left\langle\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right\rangle,$

where each summand can be decomposed as

$\displaystyle\left\langle\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right\rangle$
$\displaystyle=\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{1}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{1}}\right>+\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{2}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{2}}\right>$
$\displaystyle\hskip
11.38109pt+\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{1}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{1}}\right>+\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{2}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{2}}\right>.$

Since the parameters $W_{1}^{k}$, $W_{2}^{k}$ and $b_{1}^{k}$ are i.i.d.
random samples from $\mathcal{N}(0,1)$ for each $k=1,\cdots d_{1}$, we may
write them as random variables without the index $k$

$X_{1}:=W^{k}_{1}\sim\,\mathcal{N}(0,1),\,\,X_{2}:=W^{k}_{2}\sim\,\mathcal{N}(0,1),\,\,Z:=b^{k}_{1}\sim\,\mathcal{N}(0,1).$

Furthermore, since the random variables $X_{1}$ and $Z$ are independent, it is
easy to verify that

$Y:=X_{1}x^{j}_{b}+Z\sim\mathcal{N}(0,|x_{b}^{j}|^{2}+1).$

Then if the number of neurons in the hidden layer $d_{1}$ goes to infinity, we
have from the law of large numbers and Lemma 5.1 that

$\displaystyle\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{1}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{1}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left|W_{2}^{k}\sigma^{\prime}(W_{1}^{k}x_{b}^{j}+b_{1}^{k})x_{b}^{j}\right|^{2}\xrightarrow{\mathcal{P}}|x_{b}^{j}|^{2}\mathbb{E}\left[\sigma^{\prime}(Y)^{2}\right],$
$\displaystyle\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{2}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{dW_{2}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left|\sigma(W_{1}^{k}x_{b}^{j}+b_{1}^{k})\right|^{2}\xrightarrow{\mathcal{P}}\mathbb{E}\left[\sigma(Y)^{2}\right],$
$\displaystyle\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{1}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{1}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left|W_{2}^{k}\sigma^{\prime}(W_{1}^{k}x_{b}^{j}+b_{1}^{k})\right|^{2}\xrightarrow{\mathcal{P}}\mathbb{E}\left[\sigma^{\prime}(Y)^{2}\right],$
$\displaystyle\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{2}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{db_{2}}\right>$
$\displaystyle=1.$

By recalling the definition of the activation function
$\sigma=\max\\{0,x^{3}\\}$ and using Lemma 5.1, there holds

$\displaystyle\mathbb{E}\left[\sigma^{\prime}(Y)^{2}\right]$
$\displaystyle=\frac{9}{2\sqrt{2\pi(|x^{j}_{b}|^{2}+1)}}\int^{\infty}_{-\infty}y^{4}\exp\left(-\frac{y^{2}}{2(|x^{j}_{b}|^{2}+1)}\right)\,\mathrm{d}y=\frac{9}{2}\mathbb{E}\left[Y^{4}\right]=\frac{27}{2}\left(|x^{j}_{b}|^{2}+1\right)^{2}.$

Next, by Lemma 5.1 again, the second term can be computed as

$\displaystyle\mathbb{E}\left[\sigma(Y)^{2}\right]$
$\displaystyle=\frac{1}{2\sqrt{2\pi(|x^{j}_{b}|^{2}+1)}}\int^{\infty}_{-\infty}y^{6}\exp\left(-\frac{y^{2}}{2(|x^{j}_{b}|^{2}+1)}\right)\,\mathrm{d}y=\frac{1}{2}\mathbb{E}\left[Y^{6}\right]=\frac{15}{2}\left(|x^{j}_{b}|^{2}+1\right)^{3}.$

Altogether, we obtain as $d_{1}\rightarrow\infty$ that

$\displaystyle\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left<\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right>$
$\displaystyle\xrightarrow{\mathcal{P}}\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left(21|x^{j}_{b}|^{6}+63|x^{j}_{b}|^{4}+63|x^{j}_{b}|^{2}+21\right).$
(5.4)

Therefore, for sufficiently large $N_{b}$, by the law of large numbers and the
fact $x^{j}_{b}\sim\mathcal{U}(\partial[0,N])$, we can conclude that

$\frac{{\rm
Tr}(\boldsymbol{K}_{uu}(0))}{N_{b}}=\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left\langle\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du(x_{b}^{j};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right\rangle\approx\mathcal{O}(N^{6}).$
(5.5)

As a next step, we shall compute the value ${\rm Tr}(\boldsymbol{K}_{\rm
rr}(0))$ on the data set $\\{x_{r}^{i}\\}_{i=1}^{N_{r}}$ with
$x^{i}_{r}\sim\,{\rm{i.i.d.}}\,\,\mathcal{U}((0,N))$. By recalling the
definition (2.9) and (5.3), we see that

$\frac{{\rm
Tr}(\boldsymbol{K}_{rr}(0))}{N_{r}}=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right>$
(5.6)

where we can verify that

$u_{xx}(x_{r}^{i};\boldsymbol{\theta})=\frac{1}{\sqrt{d_{1}}}\sum_{k=1}^{d_{1}}W_{2}^{k}\sigma^{\prime\prime}(W_{1}^{k}x_{r}^{i}+b_{1}^{k})|W_{1}^{k}|^{2}.$
(5.7)

In (5.6), for each $i=1,\cdots,N_{r}$, the summand can be written as

$\displaystyle\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right>$
$\displaystyle=\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{1}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{1}}\right>+\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{2}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{2}}\right>$
$\displaystyle\hskip
11.38109pt+\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{1}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{1}}\right>+\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{2}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{2}}\right>.$

For this case, we should be more careful since $X_{1}$ and $X_{1}x^{i}_{r}+Z$
appear simultaneously in (5.7) and they may not be independent. As
$d_{1}\rightarrow\infty$, by Lemma 5.1 and the law of large numbers, we obtain

$\displaystyle\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{1}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{1}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left||W_{1}^{k}|^{2}W_{2}^{k}\sigma^{(3)}(W_{1}^{k}x_{r}^{i}+b_{1}^{k})x_{r}^{i}+2W_{1}^{k}W_{2}^{k}\sigma^{\prime\prime}(W_{1}^{k}x_{r}^{i}+b_{1}^{k})\right|^{2}$
$\displaystyle\xrightarrow{\mathcal{P}}|x^{i}_{r}|^{2}\mathbb{E}\left[X_{1}^{4}\sigma^{(3)}(X_{1}x^{i}_{r}+Z)^{2}\right]+4x_{r}^{i}\mathbb{E}\left[X_{1}^{3}\sigma^{(3)}(X_{1}x^{i}_{r}+Z)\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)\right]$
$\displaystyle\hskip
17.07164pt+4\mathbb{E}\left[X_{1}^{2}\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)^{2}\right],$
$\displaystyle\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{2}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{dW_{2}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left|\sigma^{\prime\prime}(W_{1}^{k}x_{r}^{i}+b_{1}^{k})|W_{1}^{k}|^{2}\right|^{2}\xrightarrow{\mathcal{P}}\mathbb{E}\left[X_{1}^{4}\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)^{2}\right],$
$\displaystyle\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{1}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{1}}\right>$
$\displaystyle=\frac{1}{d_{1}}\sum_{k=1}^{d_{1}}\left|W_{2}^{k}\sigma^{(3)}(W_{1}^{k}x_{r}^{i}+b_{1}^{k})|W_{1}^{k}|^{2}\right|^{2}\xrightarrow{\mathcal{P}}\mathbb{E}\left[X_{1}^{4}\sigma^{(3)}(X_{1}x^{i}_{r}+Z)^{2}\right],$
$\displaystyle\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{2}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{db_{2}}\right>$
$\displaystyle=0.$

Regarding the derivative of $\sigma$ in the sense of distribution to define
$\sigma^{(3)}$, we see by Lemma 5.1 and the fact
$|\sigma^{(3)}(X_{1}x^{i}_{r}+Z)|\leq 6$ that

$\mathbb{E}\left[X_{1}^{4}\sigma^{(3)}(X_{1}x^{i}_{r}+Z)^{2}\right]\leq\frac{36}{{2\pi}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x^{4}\exp\left(-\frac{x^{2}+z^{2}}{2}\right)\,\mathrm{d}x\,\mathrm{d}z=36\mathbb{E}\left[X_{1}^{4}\right]=108.$

Next, by Lemma 5.1 and Fubini’s theorem,

$\displaystyle\mathbb{E}\left[X_{1}^{3}\sigma^{(3)}(X_{1}x^{i}_{r}+Z)\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)\right]$
$\displaystyle\leq\frac{6}{2\pi}\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}x^{3}\sigma^{\prime\prime}(xx^{i}_{r}+z)\exp\left(-\frac{x^{2}+z^{2}}{2}\right)\,\mathrm{d}z\,\mathrm{d}x$
$\displaystyle\leq\frac{36}{2\pi}\int_{-\infty}^{\infty}\int_{-x^{i}_{r}x}^{\infty}(x_{r}^{i}x^{4}+x^{3}z)\exp\left(-\frac{x^{2}+z^{2}}{2}\right)\,\mathrm{d}z\,\mathrm{d}x$
$\displaystyle\leq\frac{36}{2\pi}\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}x^{i}_{r}x^{4}\exp\left(-\frac{x^{2}+z^{2}}{2}\right)\,\mathrm{d}z\,\mathrm{d}x$
$\displaystyle\hskip
14.22636pt+\frac{36}{2\pi}\int^{\infty}_{-\infty}x^{3}\exp\left(-\frac{x^{2}}{2}\right)\int_{-x^{i}_{r}x}^{\infty}z\exp\left(-\frac{z^{2}}{2}\right)\,\mathrm{d}z\,\mathrm{d}x$
$\displaystyle=36x_{r}^{i}\mathbb{E}\left[X_{1}^{4}\right]+\frac{36}{2\pi}\int^{\infty}_{-\infty}x^{3}\exp\left(-\frac{(|x^{i}_{r}|^{2}+1)x^{2}}{2}\right)=108x^{i}_{r}.$

Moreover, by Lemma 5.1 again and the fact
$|\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)|^{2}\leq 36|X_{1}x^{i}_{r}+Z|^{2}$,
we have

$\displaystyle\mathbb{E}\left[X_{1}^{2}\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)^{2}\right]$
$\displaystyle\leq\frac{36}{2\pi}\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}x^{2}(x^{2}|x_{r}^{i}|^{2}+2xzx_{r}^{i}+z^{2})\exp\left(-\frac{x^{2}+z^{2}}{2}\right)\,\mathrm{d}z\,\mathrm{d}x$
$\displaystyle=36|x_{r}^{i}|^{2}\mathbb{E}\left[X_{1}^{4}\right]+36\mathbb{E}\left[X_{1}^{2}\right]\mathbb{E}\left[Z^{2}\right]=36\left(3|x^{i}_{r}|^{2}+1\right),$

and similarly we obtain

$\displaystyle\mathbb{E}\left[X_{1}^{4}\sigma^{\prime\prime}(X_{1}x^{i}_{r}+Z)^{2}\right]\leq
36|x_{r}^{i}|^{2}\mathbb{E}\left[X_{1}^{6}\right]+36\mathbb{E}\left[X_{1}^{4}\right]\mathbb{E}\left[Z^{2}\right]=36\left(15|x^{i}_{r}|^{2}+3\right).$

From above, by collecting all the terms, we have as $d_{1}\rightarrow\infty$
that

$\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\left<\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right>\lesssim\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\left(|x^{i}_{r}|^{2}+1\right).$

Finally, for sufficiently large $N_{r}$, from the law of large numbers and the
fact $x^{i}_{r}\sim\mathcal{U}((0,N))$, we may conclude that

$\frac{{\rm
Tr}(\boldsymbol{K}_{rr}(0))}{N_{r}}=\frac{1}{N_{r}}\sum_{j=1}^{N_{r}}\left\langle\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}},\frac{du_{xx}(x_{r}^{i};\boldsymbol{\theta})}{d\boldsymbol{\theta}}\right\rangle\lesssim\mathcal{O}(N^{2}).$
(5.8)

Figure 12: Log-log plots of $\frac{{\rm Tr}(\boldsymbol{K}_{uu})}{N_{b}}$ and
$\frac{{\rm Tr}(\boldsymbol{K}_{rr})}{N_{r}}$ against the input samples
$x^{j}_{b}$ and $x^{i}_{r}$. We can confirm that the experimental results
align well with the theoretical results. Figure 13: The averaged convergence
rates $\frac{{\rm Tr}(\boldsymbol{K}_{uu})}{N_{b}}$ and $\frac{{\rm
Tr}(\boldsymbol{K}_{rr})}{N_{r}}$ during training.

Therefore, from (5.5) and (5.8), we have shown that the averaged convergence
rate (2.10) becomes larger when we apply the variable scaling. Figure 12
displays the averaged convergence rates $\frac{{\rm
Tr}(\boldsymbol{K}_{uu})}{N_{b}}$ and $\frac{{\rm
Tr}(\boldsymbol{K}_{rr})}{N_{r}}$ as functions of the input samples
$x^{j}_{b}$ and $x^{i}_{r}$ respectively. From the figures, we can confirm
that the experimental results go along with the theoretical findings.
Furthermore, from Figure 13, we can observe that the convergence rates of the
VS-PINN are maintained higher as N increases during the training process,
which supports the faster convergence of the VS-PINN observed in the
experiments. Note, however, that this result only deals with the theoretical
convergence rates and holds for the ideal case. In practice, large scaling
factor $N$ does not always guarantee faster convergence, and the choice of
suitable $N$ is important. We will discuss this issue in the following
section.

###  5.2 Proper choice of the scaling factor $N$

In the previous sections, we have seen that applying the variable scaling to
PINNs significantly improves the performance, especially for PDE problems with
stiff behavior and high frequency. Furthermore, throughout the NTK analysis,
we have provided some theoretical evidence of why the VS-PINN can achieve
better performance as the scaling factor $N$ increases. However, simply
increasing $N$ without careful consideration may deteriorate the performance
of the VS-PINN in practice. We shall briefly discuss this issue in this
section through a simple experiment.

Let us consider the one-dimensional Poisson problem

$\displaystyle-u_{xx}$ $\displaystyle=\pi^{2}\sin(\pi x)\;\;\;\;x\in[0,1],$
$\displaystyle u(0)$ $\displaystyle=u(1)=0,$

and compute the approximate solution using the VS-PINN with various scaling
factors. We shall use a two-layer neural network with $40000$ neurons in the
hidden layer, where the weight and bias parameters are initialized according
to the normal distribution $\mathcal{N}(0,0.1^{2})$. The training proceeds
with two boundary data points and $50$ uniformly-distributed interior data
points.

As we can see from Figure 14, when $N=2$, the performance of the VS-PINN is
better than the standard PINN. However, if we set $N=4$, we can observe that
the training becomes unstable and accuracy decreases after a certain point.
Furthermore, if we choose a much larger scaling factor $N=1000$, the
performance deteriorates even further and the VS-PINN is no longer promising
compared to the original PINN. In fact, when the scaling factor is excessively
large, the training dynamic becomes very sensitive and we need to carefully
adjust the parameters related to the training. For example, since the domain
is getting extended as $N$ increases, we need a larger number of collocation
points to guarantee satisfactory performance. This observation aligns with the
situations encountered when applying the law of large numbers in the
theoretical analysis in the previous section. In other words, if the scaling
factor is chosen appropriately, we can effectively train the model only with a
reasonable number of sample points leveraging the advantages of the VS-PINN,
while this is no longer the case for the extreme choice of $N$. Moreover, if
$N$ is too large, the corresponding scaled loss function might become
excessively large. In this case, the change of parameters in a single step of
gradient descent becomes larger, which hinders us from finding the global
minima and makes the training unstable. Therefore, as $N$ increases, we need
to put in a lot of effort to find the proper settings for training and need to
carefully adjust the factors related to training such as the weight parameters
in the loss, variance of the initialization distribution and the learning
rate. In conclusion, to effectively utilize VS-PINN without excessive effort
to find a proper circumstance for training, it is important to choose a
reasonable scaling factor $N$, which can indeed significantly improve the
performance of the PINNs as we saw in the experiments in Section 4.

Figure 14: Learning curves of the VS-PINN with various scaling factors.

##  6 Conclusion

In this paper, we proposed a novel training method for PINNs based on the
variable-scaling technique. This method is simple and it can be applied to a
wide range of PDE problems with stiff behavior or high frequency. We applied
this new strategy to various problems including the wave equation, the Allen-
Cahn equation, the singular perturbation problem and the Navier–Stokes
equations. Throughout these numerical experiments, we confirmed that the VS-
PINN provides a significant performance improvement both in accuracy and
training efficiency for stiff problems and even for nonlinear problem.
Moreover, we performed the NTK analysis to provide theoretical evidence of why
the VS-PINN works well. More precisely, we explicitly computed the eigenvalues
of the NTK of VS-PINNs which concerns the rate of convergence for training,
and showed that the eigenvalues corresponding to the VS-PINN are indeed larger
when the variable is scaled.

An interesting future research direction is to theoretically analyze the
convergence properties of the proposed method from the viewpoint of
approximation and generalization. Furthermore, the systematic guidelines for
choosing the scaling factor $N$ and the parameters related to training should
also be investigated. Since increasing $N$ causes instability in training,
establishing a general strategy that can address this issue while further
enhancing the overall performance of the VS-PINN is an important research
question. Since the choice of the parameters depends on the structure of the
equation under consideration, we will provide a more detailed answer to this
question with a general principle and theoretical intuition by examining a
wider variety of PDE problems using the VS-PINN, which will be addressed in
the forthcoming paper.

## References

  * [1] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind.  Automatic differentiation in machine learning: a survey.  J. Mach. Learn. Res., 18:153:1–153:43, 2017.
  * [2] Y. Cao, Z. Fang, Y. Wu, D.-X. Zhou, and Q. Gu.  Towards understanding the spectral bias of deep learning.  2020\.
  * [3] Y. Chen, L. Lu, G. E. Karniadakis, and L. D. Negro.  Physics-informed neural networks for inverse problems in nano-optics and metamaterials.  Opt. Express, 28(8):11618–11633, Apr 2020.
  * [4] A. Farkane, M. Ghogho, M. Oudani, and M. Boutayeb.  Epinn-nse: Enhanced physics-informed neural networks for solving navier-stokes equations.  2023\.
  * [5] Y. H. Gung-Min Gie and C.-Y. Jung.  Semi-analytic pinn methods for singularly perturbed boundary value problems.  Applicable Analysis, 0(0):1–18, 2024.
  * [6] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [7] R. Hogg, J. McKean, and A. Craig.  Introduction to Mathematical Statistics.  Pearson, 2019.
  * [8] A. Jacot, F. Gabriel, and C. Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
  * [9] X. Jin, S. Cai, H. Li, and G. E. Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics, 426:109951, 2021.
  * [10] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [11] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, and P. Perdikaris.  Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks.  Computer Methods in Applied Mechanics and Engineering, 358:112623, 2020.
  * [12] D. Liu and Y. Wang.  Multi-Fidelity Physics-Constrained Neural Network and Its Application in Materials Modeling.  Journal of Mechanical Design, 141(12):121403, 09 2019.
  * [13] L. Lu, M. Dao, P. Kumar, U. Ramamurty, G. E. Karniadakis, and S. Suresh.  Extraction of mechanical properties of materials through deep learning from instrumented indentation.  Proceedings of the National Academy of Sciences, 117(13):7052–7062, 2020.
  * [14] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  DeepXDE: a deep learning library for solving differential equations.  SIAM Rev., 63(1):208–228, 2021.
  * [15] L. D. McClenny and U. M. Braga-Neto.  Self-adaptive physics-informed neural networks.  Journal of Computational Physics, 474:111722, 2023.
  * [16] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  2017\.
  * [17] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville.  On the spectral bias of neural networks.  In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5301–5310. PMLR, 09–15 Jun 2019\.
  * [18] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [19] M. Raissi, Z. Wang, M. S. Triantafyllou, and G. E. Karniadakis.  Deep learning of vortex-induced vibrations.  Journal of Fluid Mechanics, 861:119–137, 2019.
  * [20] M. Raissi, A. Yazdani, and G. E. Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  Science, 367(6481):1026–1030, 2020.
  * [21] C. Rao, H. Sun, and Y. Liu.  Physics-informed deep learning for incompressible laminar flows.  Theoretical and Applied Mechanics Letters, 10(3):207–212, 2020\.
  * [22] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, and E. Kuhl.  Physics-informed neural networks for cardiac activation mapping.  Frontiers in Physics, 8, 2020.
  * [23] K. Shukla, P. C. D. Leoni, J. L. Blackshire, D. Sparkman, and G. E. Karniadakis.  Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks.  CoRR, abs/2005.03596, 2020.
  * [24] J. Sirignano and K. Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [25] M. Tancik, P. P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. T. Barron, and R. Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  2020\.
  * [26] S. Wang, S. Sankaran, and P. Perdikaris.  Respecting causality for training physics-informed neural networks.  Computer Methods in Applied Mechanics and Engineering, 421:116813, 2024.
  * [27] S. Wang, Y. Teng, and P. Perdikaris.  Understanding and mitigating gradient flow pathologies in physics-informed neural networks.  SIAM Journal on Scientific Computing, 43(5):A3055–A3081, 2021.
  * [28] S. Wang, X. Yu, and P. Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  Journal of Computational Physics, 449:110768, 2022.
  * [29] L. Yang, X. Meng, and G. E. Karniadakis.  B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data.  Journal of Computational Physics, 425:109913, 2021.
  * [30] Y. Yang and P. Perdikaris.  Adversarial uncertainty quantification in physics-informed neural networks.  Journal of Computational Physics, 394:136–152, 2019.
  * [31] B. Yu et al.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [32] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, and P. Perdikaris.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  Journal of Computational Physics, 394:56–81, 2019.