# Extreme learning machine collocation for the numerical solution of elliptic
PDEs with sharp gradients

Francesco Calabrò<EMAIL_ADDRESS>Gianluca Fabiani
<EMAIL_ADDRESS>Constantinos Siettos<EMAIL_ADDRESS>Dipartimento di Matematica e Applicazioni “Renato Caccioppoli”, Università
degli Studi di Napoli Federico II  Scuola Superiore Meridionale di Napoli,
Università degli Studi di Napoli Federico II

###### Abstract

We introduce a new numerical method based on machine learning to approximate
the solution of elliptic partial differential equations with collocation using
a set of sigmoidal functions. We show that a feedforward neural network with a
single hidden layer with sigmoidal functions and fixed, random, internal
weights and biases can be used to compute accurately a collocation solution.
The choice to fix internal weights and bias leads to the so-called Extreme
Learning Machine network. We discuss how to determine the range for both
internal weights and biases in order to obtain a good underlining
approximating space, and we explore the required number of collocation points.
We demonstrate the efficiency of the proposed method with several one-
dimensional diffusion-advection-reaction problems that exhibit steep
behaviors, such as boundary layers. The boundary conditions are imposed
directly as collocation equations. We point out that there is no need of
training the network, as the proposed numerical approach results to a linear
problem that can be easily solved using least-squares. Numerical results show
that the proposed method achieves a good accuracy. Finally, we compare the
proposed method with finite differences and point out the significant
improvements in terms of computational cost, thus avoiding the time-consuming
training phase.

###### keywords:

Partial differential equations, Collocation methods, Artificial Neural
Networks, Extreme Learning Machine, Boundary layer, Sigmoidal functions.

##  1 Introduction

Over the last few years various algorithms have been addressed for the
numerical solution of Partial Differential Equations (PDEs) with the aid of
Artificial Neural Networks (ANN) (see e.g. [6, 17, 38, 41]). The use of ANN
for the numerical solutions of systems of differential equations goes back to
[27, 28], and the interest on new applications is in nowadays at the forefront
of research interest [2, 3, 8, 13, 15, 18, 29, 31, 32, 35, 40]. The main focus
of such techniques is on the numerical solution of some very difficult
problems: for example high dimensional systems with difficult geometries and
non-linear. In all cases, ANN seem to outperform conventional techniques being
extremely flexible and computationally fast. The celebrated theorem of the
universal approximation of ANNs gives some justification to such behavior,
although some results are still unavailable: i.e. it is proved that the
approximation error vanishes asymptotically, but there are few available
results on explicit error bounds for the approximation error or about
polynomial reproduction (see e.g. [1, 9, 14, 36]). Among all types of ANNs,
recently great attention is given to the so-called Extreme Learning Machines
(ELM) [22, 24]. However, the efficiency and applicability of ELM for the
approximation of differential problems is still unexplored. Indeed, to the
best of our knowledge, the only study on the subject is that of [13] where the
authors however report a failure of ELM to deal with sharp gradient problems.

Here, we show how ELM can be exploited to deal with boundary problems, in
particular 1D second order elliptic linear equations that exhibit steep
behaviors, such as boundary layers. It is known that such problems if solved
by classical methods such as finite differences may lead to several numerical
instabilities such as spurious oscillations (see e.g. [11, 34]) even if the
solution is regular. Thus, we focus on linear advection-reaction problems with
constant coefficients where the solution is known analytically and it is
regular, but the advection-dominated or the reaction-dominated cases may
result to boundary layers. For such problems, several approaches (for example
artificial diffusion, upwind schemes or mesh adaptivity [33, 37]) have been
addressed to deal with the emerged steep gradients. Compared to the above
techniques, our approach is free from problem-dependent modifications being
able to detect the steep gradients that arise. The results reveal how ELM can
serve as a robust numerical approach to obtain accurate solutions without the
need of stabilization techniques. In particular, we propose an under-
determined collocation method applied to an ELM with one hidden layer which
they don’t need to be trained. Our findings reveal the ability of ELM to
approximate the solutions efficiently, thus extending and giving a new insight
on the use of such types of ANN. What we propose here, by taking both internal
weights and biases a-priori fixed is that we do not need to train the network
in order to obtain a sufficiently accurate solution, at least if the number of
neurons fits the scale of the layer. Obviously, such solutions can be taken as
initial choices for a training process if the problem is non-linear and/or
demands a higher approximation accuracy. The paper is organized as follows. In
Section 2, we introduce our approximation space on the base of ANN and ELM,
and in Section 3 we present the proposed ELM collocation method. The
approximation efficiency and convergence property of the scheme is analysed in
Section 4. In Section 5, we present the numerical results obtained the
proposed approach for the solution of several 1D boundary-layer elliptic
problems and compare them with both the exact-analytical and FD numerical
solution.

##  2 The proposed ANN architecture: the ELM network

In this section, we briefly introduce the ANN that we propose and the main
results that give the properties of such universal approximators. We consider
an ANN with a single hidden layer with $n$ neurons. We choose the activation
function to be the standard logistic sigmoid function defined by

$\sigma(x)=\frac{1}{1+\text{exp}(-x)}\ .$

The output function is the linear activation function, so that the final
approximation is achieved via a linear combination of sigmoid functions. In
the network one has to fix the biases $\beta$, the connections between the
input layer and the hidden layer ${\alpha}$, represented by a vector of
internal weights, and the connections of the hidden layer with the output
layer ${w}$, denoted by the external weights. Therefore, this type of ANN
describes a map $G:\mathbb{R}\rightarrow\mathbb{R}$ that can be written as:

$G(x;\bm{\alpha},\bm{\beta},\bm{w})=\sum_{i=1}^{n}w_{i}\sigma(\alpha_{i}x+\beta_{i}).$
(2.1)

As it is well-known, this kind of ANN is a universal approximator for any
$L^{1}$ function as stated, for example in [10, 19, 20].

######  Theorem 2.1.

For any function $f\in L^{1}([a,b])$ and for all $\varepsilon>0$ there exist a
choice of $\bm{\alpha},\bm{\beta},\bm{w}$ such that

$\|G-f\|_{L^{1}}<\varepsilon.$

where $\|\cdot\|_{L^{1}}$ is the usual $L^{1}$ norm, i.e.
$||\psi||_{L^{1}}=\int_{a}^{b}|\psi(x)|dx$.

Similar results apply for continuous functions and for the derivatives of
differentiable functions [4, 19, 36].

These results open the way to the use of ANN for the numerical solution of
differential equations. The most frequent approach when using ANN for this
purpose is the minimization of a cost function in order to compute the optimal
parameters $\bm{\alpha},\bm{\beta},\bm{w}$. This optimization procedure is
referred in ANN as training. The optimization method can be also exploited for
the resolution of differential problems because for many of these, the exact
jacobian can be computed easily, see e. g. [27]. Nevertheless, two main
questions raise when training an ANN via optimization: the initial choice of
the parameters and the overall computational cost. When dealing with time-
dependent problems, the initial choice is related to the previous time step
(see e.g. [38]); in other problems where the main difficulties are related to
the geometry of the boundary, the initial choice is efficiently chosen as a
boundary lift (see e.g. [6, 27]). In the problem we consider here, there is no
such a good indication by the problem itself that may help to choose “good”
initial guesses. Here, we present a procedure that computes a solution in the
framework of ELM networks that do not need training unless a big accuracy is
needed; in the later case it can be used for providing an initial guess for
training purposes. Toward to the above aim, we first show how one can set the
internal weights and biases to get good approximating functions.

Consider the function $G$ in (2.1) that depends non-linearly on the internal
weights and biases. Let us now take the function:

$\sigma_{i}(x)=\sigma(\alpha_{i}x+\beta_{i})=\frac{1}{1+\text{exp}(-\alpha_{i}x-\beta_{i})}\
.$ (2.2)

The first and second derivatives of the above with respect to the independent
variable $x$ read:

$\begin{split}\sigma_{i}^{\prime}(x)=&\alpha_{i}\frac{\text{exp}(z)}{(1+\text{exp}(z))^{2}}\\\
\sigma_{i}^{\prime\prime}(x)=&\alpha_{i}^{2}\frac{\text{exp}(z)(1-\text{exp}(z))}{(1+\text{exp}(z))^{3}},\end{split}$
(2.3)

where $z=\alpha_{i}x+\beta_{i}$. Note, that if one takes two functions
$\sigma_{i}$ where at least one of the parameters $\alpha,\beta$ is different,
then these are linearly independent (see [26]). Finally, the ratio between
$\alpha_{i}$ and $\beta_{i}$ gives the location of the inflection point, while
the internal weights $\alpha_{i}$ govern the variation of the amplitude of the
S–shape:

  * 1.

$\sigma_{i}$ has an inflection point at $x=-\dfrac{\beta_{i}}{\alpha_{i}}$,
that we call the center $C_{i}$ of the sigmoid function;

  * 2.

$\sigma_{i}$ is monotone, $\lim_{x\to-\infty}=0;\lim_{x\to+\infty}=1$ if
$\alpha_{i}$ is positive, the other way if negative. Moreover the range where
the values are between 0.05 and 0.95 is
$\left[C_{i}-\frac{2.945}{\alpha_{i}},C_{i}+\frac{2.945}{\alpha_{i}}\right]$.
We denote this interval as the S–shape amplitude $I_{S}$ of the function.

Figure 1: The functions $\sigma_{i}$ of (2.2) with varying parameters. On the
top left panel, we have set $\alpha_{i}=100$, on the top right,
$\alpha_{i}=1$. On the bottom left functions with a fixed center $C_{i}=1/2$
are depicted. On the bottom right, a set of 10 functions obtained by varying
the coefficients as stated in Section 3.

If $\alpha_{i},\beta_{i}$ are fixed once randomly, we keep their values
constant, without updating them by tuning. This is the case of a ANN network
referred as Extreme Learning Machine (ELM) network [23, 24, 25]. The only free
parameters need to be learned are the connections (weights) between the hidden
layer and the output layer. The analysis reviewed in [22] gives the conditions
that are required for ensuring the interpolation property, thus the non-
singularity of the collocation matrices; we report this result in Section 4.

Outside $I_{S}$, the function $\sigma_{i}$ is almost constant; then, we focus
on the values of the parameters $\alpha_{i},\beta_{i}$ such that the
intersection of the S–shape amplitude $I_{S}$ and the domain of interest is
non-null, thus leading to functions $\sigma_{i}$ with $C_{i}$ inside the
interval. The choice of $\alpha_{i}$ is related to the behavior of the
derivative of the function. When taking all internal weights $\alpha_{i}\equiv
K$, one obtains Heaviside-like functions if $K$ is big (as used in [16]), or
almost-linear functions if $K$ is small (as used in [14, 30]). In our case,
the use of both kind of functions can help to approximate both steep gradients
and global behaviors. In figure 1, we plot different cases of such functions:
on the top panels $\alpha_{i}$ is a constant, on the bottom left panel $C_{i}$
is set equal to the mid-point of the interval and finally on the bottom right
panel, we plot functions with random values of $\alpha_{i},\beta_{i}$. Note
that some functions can have similar shapes: interpolation properties are not
changed but the problem can become ill-conditioned.

##  3 The proposed numerical algorithm

In this section, we propose a simple algorithm to get a numerical
approximation of a linear differential problem with the aid of ELM networks
with $n$ neurons in the hidden layer. The main idea is to fix the range of
values of weights and biases, so that we construct a “good” basis of $n$
sigmoid functions $\sigma_{i}$, and use a collocation method to fix the
weights of the linear output activation function.

For our illustrations, we have chosen a stationary one-dimensional PDE with
constant coefficients.

$\left\\{\begin{array}[]{l}-\mu u^{\prime\prime}(x)+\gamma
u^{\prime}(x)+\lambda u(x)=f(x),\qquad x\in I:=[0,1]\subseteq\mathbb{R}\\\
\nu_{0}u^{\prime}(0)+\rho_{0}u(0)=g_{0}\,,\
\nu_{1}u^{\prime}(1)+\rho_{1}u(1)=g_{1},\end{array}\right.$ (3.1)

where $\mu$, $\gamma$ and $\lambda$ are respectively the diffusion, advection
and reaction terms. The domain is fixed to $[0,1]$ for convenience. The
boundary conditions are in the general Robin form; Dirichlet and Neumann
boundary conditions are derived by setting $\nu_{i}=0$ or $\rho_{i}=0$.

Input:  coefficients $\mu,\gamma,\lambda;\nu,\rho$ ; number of neurons $n$,
number of collocation points $M$, function evaluations
$[\bm{f},\bm{g}]\in\mathbb{R}^{M}$

1\. set $\bm{\alpha}\in\mathbb{R}^{n}$ randomly according to (3.3) ;

2\. set $\bm{\beta}\in\mathbb{R}^{n}$ such that $C_{i}$ are uniformly
distributed, i.e. $\beta_{i}=-\alpha_{i}*C_{i}$;

3\. compute matrix $\mathbb{S}$ according to (3.5)-(3.6) ;

4\. compute $\mathbb{B}$ according to (3.7)-(3.8) ;

5\. solve the (eventually under/over -determined) linear system
$[\mathbb{S};\mathbb{B}]\bm{w}=[\bm{f};\bm{g}]$ ;

Output: weights $\bm{w}\in\mathbb{R}^{n}$ that fix the ELM approximate
solution $\tilde{u}(x)$ in (3.2);

Algorithm 1 Pseudo-code that describe the ELM collocation method

We consider an approximate solution $\tilde{u}(x,\bm{w})$ given by the ELM:

$\tilde{u}(x;\bm{w})=\sum_{i=1}^{n}w_{i}\sigma_{i}(x).$ (3.2)

We choose $\bm{\alpha}$ to be random uniformly distributed in a range that
depends on the number of neurons $n$; if we fix the domain of the differential
problem to have unitary length, we choose:

$\bm{\alpha}:={\text{rand}\left(\left[-\frac{n-10}{10}-4,\frac{n-10}{10}+4\right]\right)}\
.$ (3.3)

In this way, if $n=10$ we have a maximum weight e.g. $4$ and as $n$ grows by
10 then the maximum weight grows by $1$. The biases $\bm{\beta}$ are fixed so
that the centers $C_{i}=-\beta_{i}/\alpha_{i}$ of the sigmoid functions are
located on random uniformly distributed points inside the interval.

Then, the ELM network is collocated in (3.1) on $M-2$ equidistant collocation
points $x_{j}\in(0,1)$ getting:

$-\mu\tilde{u}^{\prime\prime}(x_{j})+\gamma\tilde{u}^{\prime}(x_{j})+\lambda\tilde{u}(x_{j})=f(x_{j})\qquad\text{for
}j=2,\dots,M-1.$ (3.4)

Analytical derivatives of the sigmoid function $\sigma_{i}$ can be computed by
(2.3). Then:

$-\mu\sum_{i=1}^{n}w_{i}\sigma_{i}^{\prime\prime}(x_{j})+\gamma\sum_{i=1}^{n}w_{i}\sigma_{i}^{\prime}(x_{j})+\lambda\sum_{i=1}^{n}w_{i}\sigma_{i}(x_{j})=f(x_{j})\
\text{for }j=2,\dots,M-1.$

Finally, we observe that we can write separately the external weights $\bm{w}$
and write this in a matrix form as:

$\mathbb{S}\bm{w}:=(-\mu\mathbb{S}^{(2)}+\gamma\mathbb{S}^{(1)}+\lambda\mathbb{S}^{(0)})\bm{w}=\bm{f},$
(3.5)

where $\mathbb{S}_{2}$, $\mathbb{S}_{1}$ and $\mathbb{S}_{0}$ are matrices of
elements:

$\begin{split}\mathbb{S}^{(2)}=\bigg{(}\sigma_{i}^{\prime\prime}(x_{j})\bigg{)}_{i,j}\
,\ \mathbb{S}^{(1)}=\bigg{(}\sigma_{i}^{\prime}(x_{j})\bigg{)}_{i,j}\
\text{and}\ \mathbb{S}^{(0)}=\bigg{(}\sigma_{i}(x_{j})\bigg{)}_{i,j}\\\
i=1,\dots,n\,,\ j=2,\dots,M-1\end{split}$ (3.6)

For the boundary conditions defined in (3.1), we augment the system by two
equations collocating on the boundary points $0$ and $1$

$\nu\sum_{i=1}^{n}w_{i}\sigma_{i}^{\prime}(x_{k})+\rho\sum_{i=1}^{n}w_{i}\sigma_{i}(x_{k})=g(x_{k})\qquad
k=0,M\,;\ x_{1}=0,x_{M}=1\ .$

As before, we rewrite the above in a matrix form:

$\mathbb{B}\bm{w}:=(\nu\mathbb{B}^{(1)}+\rho\mathbb{B}^{(0)})\bm{w}=\bm{g}\,,$
(3.7)

where $\mathbb{B}^{(1)}$ and $\mathbb{B}^{(0)}$ are matrices of elements:

$\begin{split}\mathbb{B}^{(1)}=\bigg{(}\sigma_{i}^{\prime}(x_{k})\bigg{)}_{i,{k}}\
\text{and}\ \mathbb{B}^{(0)}=\bigg{(}\sigma_{i}(x_{k})\bigg{)}_{i,{k}}\\\
i=1,\dots,n\,,\ k=0,M\ .\end{split}$ (3.8)

To this end, we have a linear system of $M$ equations and $n$ unknowns. The
overall system can be constructed to be over-determined ($M>n$), under-
determined ($M<n$) or squared ($M=n$). In the over- and under-determined
cases, we consider the solution to be the one obtained in the least-square
sense (in the under-determined problem the minimum-norm regularized solution).
In the next section, we show that the under-determined solutions reach the
same accuracy when compared with squared or over-determined cases. Moreover,
we observe that the squared case can lead to ill-conditioned linear problems,
so we propose under-determined collocation.

##  4 Analysis of the collocation approximation

In this section, building on previous works [5, 21, 25], we prove the
consistency of the proposed ELM collocation method. First, based on the
Theorem 2.1 presented in [25], we state the following Theorem that fits to our
proposed framework.

######  Theorem 4.2.

Let $(x_{i},y_{i})\,,\ i=1,\dots,M$ be a set of points such that
$x_{i}<x_{i+1}$, and take the ELM network with $n<M$ neurons
$\tilde{u}(x;\bm{w})$ in (3.2) such that the internal weights $\bm{\alpha}$
and the biases $\bm{\beta}$ are randomly generated independently from the data
according to any continuous probability distribution. Then,
$\forall\varepsilon>0$, there exists a choice of $\bm{w}$ such that
$\|(\tilde{u}(x_{i};\bm{w})-y_{i})_{i}\|<\varepsilon$ . Here $\|\cdot\|$
denotes the $L^{2}$ Euclidean norm of vectors, the analogous of the Frobenius
norm.

Moreover, if $n=M$ then $\bm{w}$ can be found such that
$\|(\tilde{u}(x_{i};\bm{w})-y_{i})_{i}\|=0$.

The above theorem states that if the number of hidden neurons are at least
equal to the number of data points, then the interpolation error is zero.
Starting from the above theorem, we consider the properties of ELM network for
the approximation of functions. To do this, we consider the construction of a
sequence of nested ELM families: $\mathbb{E}^{(n)}:=\\{\tilde{u}(x;\bm{w})\,,\
\bm{w}\in\mathbb{R}^{n}\\}\subset\mathbb{E}^{(n+1)}:=\\{\tilde{u}(x;\bm{w})\,,\
\bm{w}\in\mathbb{R}^{n+1}\\}$, where we add one neuron at time.

Based on the above theorem, we can now state the following theorem (see e.g.
Theorem 2 of [22]).

######  Theorem 4.3.

Let $\phi(x)$ be a continuous function. Then there exist a sequence of ELM
network functions $\tilde{u}^{(n)}\in\mathbb{E}^{(n)}$ such that:

$\|\phi-\tilde{u}^{(n)}\|_{L^{2}}\to 0\,,$

where by $\|\cdot\|_{2}$ we denote the $L^{2}$ norm.

To prove the above, we consider a family of nested sets of ordered distinct
points $\\{x^{(M)}_{i}\,,\ i=1,\dots,M\\}\subset\\{x^{(M+1)}_{i}\,,\
i=1,\dots,M+1\\}$ and call diameter of the set of points $H>0$ the maximum
$H=max_{i=1,\dots,M-1}(x^{(M)}_{i+1}-x^{(M)}_{i})$. Then take in the space
$\mathbb{E}^{(n)}$ and call $\tilde{u}^{(n)}$ the ELM network that
interpolates $\phi$ on a set of points of diameter $H$: from Theorem 4.2 this
construction exists. Finally, we can say that
$\|(\phi(x_{i})-\tilde{u}^{(n)}_{H}(x_{i}))\|\to_{H\to 0}0$, using the vector
norm: theorem 4.3 extends the convergence to a general continuous function
$\phi$ in the $L^{2}$ sense, thus states that the finite-dimensional discrete
space of ELM networks is consistent in the $L^{2}$ setting.

Theorem 4.2 and 4.3 also imply that the collocation matrices
$\mathbb{S}^{(0)}$ are of full rank.

We should note that for the problems that we consider here, in some cases we
obtain ill-conditioned matrix problems. In these cases the solution of the
under-determined collocation (number of neurons greater then number of points)
with least squares is more suitable.

The collocation solution can now be sought as a solution of a variational
formulation tested against delta-type functions or, equivalently, a
discretized formulation of the problem with $C^{2}$ test functions. We follow
the classic approach, see e.g. [5, 42] and write the weak formulation for the
equation (3.1) that seeks a $u(x)$ in the trial space (denoted by $U$) such
that:

$\int_{0}^{1}\left[-\mu u^{\prime\prime}(x)+\gamma u^{\prime}(x)+\lambda
u(x)\right]\,v(x)\ dx=\int_{0}^{1}f(x)v(x)\ dx$

for all test function $v\in V$. For simplicity, we suppose that the boundary
conditions are imposed on the space $U$ and are satisfied in all the subspaces
introduced.

The general Galerkin method is a weak formulation as above where we look for a
function $\tilde{u}_{n}\in\mathbb{U}_{n}$ for which the equations hold true
$\forall v_{m}\in\mathbb{V}_{m}$ and both $\mathbb{U}_{n}$ and
$\mathbb{V}_{m}$ are finite-dimensional. In our case, we fix
$\mathbb{U}_{n}=\mathbb{E}^{(n)}$, i.e. the ELM space with $n$ neurons, so
that the degrees of freedom of the finite dimensional space are the weights
$\bm{w}\in\mathbb{R}^{n}$.

Now, consider the test space to
$\mathbb{V}_{m}=span\\{\delta_{i,\varepsilon}(x)\,,\ i=1,\dots,m\\}$ defined
starting from a set of $m=n$ distinct points $\\{x_{i}\\}$. We also consider
function $f$ to be approximated by its interpolant in the space of ELM, i.e.
$\tilde{f}\approx f\,,\ \tilde{f}\in\mathbb{E}^{(n)}\,,\
\tilde{f}(x_{i})=f(x_{i})$. Then the Galerkin problem reads:

$\begin{split}\int_{0}^{1}\left[-\mu\tilde{u}_{\varepsilon}^{\prime\prime}(x)+\gamma\tilde{u}_{\varepsilon}^{\prime}(x)+\lambda\tilde{u}_{\varepsilon}(x)\right]\,\delta_{x_{i},\varepsilon}(x)\
dx=\int_{0}^{1}\tilde{f}(x)\delta_{x_{i},\varepsilon}(x)\ dx\\\ \quad\forall
i=1,\dots,M\end{split}$ (4.1)

where functions $\delta$ are defined by

$\displaystyle\delta_{x_{i},\varepsilon}(x)=\phi(|x-x_{i}|/\varepsilon)\,,$
$\displaystyle\phi(\zeta)=\left\\{\begin{array}[]{ll}\dfrac{exp\\{1/(\zeta^{2}-1)\\}}{\int
exp\\{1/(\zeta^{2}-1)\\}\,d\zeta}&\zeta\in[0,1)\\\
0&\zeta\in[1,+\infty)\end{array}\right..$

This Galerkin problem is related to the above collocation problem by the
following theorem.

######  Theorem 4.4.

Let $\tilde{u}_{n,\varepsilon}(x)$ be a solution to the Galerkin problem
(4.1). Then $\lim_{\varepsilon\to
0}\tilde{u}_{n,\varepsilon}(x)=\tilde{u}_{n}(x)$ where $\tilde{u}_{n}(x)$ is
the solution in the ELM space of the collocation problem (3.4) on the points
$\\{x_{i}\\}$:

$-\mu\tilde{u}_{n}^{\prime\prime}(x_{i})+\gamma\tilde{u}_{n}^{\prime}(x_{i})+\lambda\tilde{u}_{n}(x_{i})=\tilde{f}(x_{i})\qquad\text{for
}i=1,\dots,M.$ (4.2)

This result follows from convergence of the integrals to evaluations of
function on sites, due to the convergence of $\delta$ functions to $\delta$
distributions; a proof of this can be found in [5].

Based on the Theorem 2.2 presented in [42], we are now ready to prove the
following convergence theorem for our proposed scheme.

######  Theorem 4.5.

Let $\\{x^{(n)}_{i}\,,\ i=1,\dots,n\\}_{n\in\mathbb{N}}$ be an increasing
family of collocation points such that the diameter $H$ tends to zero.
Moreover take $\mathbb{E}^{(n)}$ a family of ELM. where we assume that matrix
$\mathbb{S}$ of equation (3.5) constructed on the collocation points is full
rank and uniformly invertible.
Then, if $\tilde{u}_{n}\in\mathbb{E}^{(n)}$ is ELM network solution to the
collocation problem (4.2) and $u(x)$ is the solution to equation (3.1), we
have that:

$\lim_{n\to\infty}\|\tilde{u}_{n}-u\|_{L^{2}}=0$

###### Proof.

Denote by $L$ the differential operator in equation (3.1), such that $Lu=f$.
Call $Q_{n}:\phi\in\mathbb{R}^{n}$ the operator that gives evaluation of
function collocated on points $\\{x_{i}\\}$, i.e.
$Q_{n}(\phi)=(\phi(x_{i}))_{i}$. Then call $P_{n}:\phi\in\mathbb{R}^{n}$ the
operator that gives the ELM network that interpolates functions at collocation
nodes, i.e. $P_{n}(\phi)=\tilde{\phi}:=\sum w_{j}\sigma_{j}(x)$ such that
$\tilde{\phi}(x_{i})=\phi(x_{i})$. Then the solution $\tilde{u}_{n}$ to the
collocation problem can be seen as the solution to
$Q_{n}LP_{n}\tilde{u}_{n}=Q_{n}f$, and matrix $\mathbb{S}$ of equation (3.5)
is the discrete representation of the operator $Q_{n}LP_{n}$. With this we
have:

$\displaystyle Lu=f\Rightarrow Q_{n}Lu=Q_{n}f\Rightarrow$ $\displaystyle
Q_{n}Lu-Q_{n}LP_{n}u+Q_{n}LP_{n}u=Q_{n}f\Rightarrow$ $\displaystyle
Q_{n}LP_{n}u+Q_{n}L(u-P_{n}u)=Q_{n}f\Rightarrow$ $\displaystyle
u-(Q_{n}LP_{n})^{-1}Q_{n}f=(Q_{n}LP_{n})^{-1}Q_{n}L(P_{n}u-u)$

At last notice that from $Q_{n}LP_{n}\tilde{u}_{n}=Q_{n}f$ we can write
$\tilde{u}_{n}=(Q_{n}LP_{n})^{-1}Q_{n}f$. Finally:

$\|u-\tilde{u}_{n}\|\leq\|(Q_{n}LP_{n})^{-1}\|\,\|Q_{n}\|\,\|L\|\,\|u-P_{n}u\|\
.$

In this, we have that the first term is bounded from the hypothesis on the
operator $\mathbb{S}^{-1}$; the operator $Q_{n}$ is trivially bounded in our
setting; the linearity of operator $L$ imply that also the third term is
bounded. The thesis then follows from the convergence properties stated in
theorems 4.2-4.3. ∎

We remark that the collocation equations can be derived also via quadrature
rules applied to equations (4.1). In this case there is no need of Theorem
4.4, and the relation with collocation is obtained by the requirement that the
quadrature rule is supported on the collocation points. In this case, the
functions $\delta$ can be substituted with any family of locally supported
functions whose support includes only one collocation point.

##  5 Numerical results

In this section, we apply the proposed method for the solution of benchmark
problems exhibiting steep gradients as suggested in [33]. More particularly,
in Section 5.1, we first consider regular problems, in Section 5.2 boundary
layer problems, and in Section 5.3 internal layer problems. For our
computations, we have used Matlab2020b. In order to get reproducible results,
we have used rng(5) to create the random values for the vectors $\bm{\alpha}$,
$\bm{\beta}$ before calling the random value generator of Matlab (rand). The
solution of the least squares problems is achieved via the backslash command.

For each problem, the number $n$ of neurons ranged from 10 to 1280, doubling
the number of neurons at each execution. For each fixed $n$, we compute
solutions with various numbers $M$ of collocation points. The computational
time for all the above problems was of the order of 0.5 sec on a single intel
core i7-10750H with 16GB RAM 2.60Ghz. In particular the maximum computational
time was 0.6 sec for $n=1280$. To evaluate the approximation error, we used
two different metrics: the absolute difference between computed and exact
solution, denoted by $E_{u}$, and the residual of the equation, denoted by
$R_{f}$:

$\begin{split}E_{u}(x)=&|u(x)-\tilde{u}(x)|\\\
R_{f}(x)=&|f(x)+\mu\tilde{u}^{\prime\prime}(x)-\gamma\tilde{u}^{\prime}(x)-\lambda\tilde{u}(x)|\
.\end{split}$

Note, that all errors are absolute, not normalized. Then, we computed the
$L^{2}$ norm of $E_{u}$ and $R_{f}$; when computing the $L^{2}$-norm, we used
5000 equispaced points and the trapezoidal rule for integration. As reference
numerical method, we considered the finite difference (FD) scheme obtained
with a maximal order on 7 nodes computed on $M$ equispaced points. The FD
function considered here is the piecewise linear polynomial one interpolating
the computed values. We emphasize that in all computations the number of
points considered for the FD solution is equal to $M$, i.e. the number of
collocation equations considered for the ELM implementation. In two cases, we
also considered the solution only at collocation points, thus considering
$L^{\infty}$ errors between vectors in $\mathbb{R}^{M}$.

###  5.1 Regular problems

In this section, we analyze sinusoidal-bump problems as well as a high-order
polynomial problem to show that the proposed approach can result to a high
approximation accuracy.

The first example is a simple 1D boundary value problem with homogeneous
Dirichlet boundary conditions:

$\left\\{\begin{array}[]{l}u^{\prime\prime}+(4k^{2}\pi^{2}-1)u=4k\pi
e^{x}\cos(2k\pi x),\qquad 0<x<1\\\ u(0)=0\,,\ u(1)=0\qquad\text{ with
}k\in\mathbb{N}\end{array}\right.$ (5.1)

The above has the exact analytical solution:

$u=\text{exp}(x)\sin(2k\pi x)$

The coefficient $k$ represents the number of oscillations in the domain. We
consider the simple case with $k=1$ and a more challenging one with $k=5$ (see
Figure 2. The computed errors with respect to the exact solution are reported
in Figure 2 and in Table 1. In both cases, we note that the proposed ELM
network results to a good approximation with a modest number of neurons.

Figure 2: Applying the proposed ELM network for the numerical solution of the
problem given in (5.1): approximation with $n$ neurons, the number of
collocation points $M$ fixed to be $n/2$ (on the left panel, $k=1$, on the
right panel, $k=5$.)

Figure 3: Error and residual convergence for the ELM numerical solution of the
problem (5.1). On the top panels, $k=1$: on the left top panel $M=n/2$, on the
right top panel $n=40$. On the bottom panels, $k=5$: on the left bottom panel,
$M=n/2$, on the right bottom panel, $n=160$.

Figure 4: Converge plots of the $L^{\infty}$ norm at the collocation points
for the solution of the problem (5.1): on the left panel, $k=1$, on the right
panel, $k=5$.

k=1 |  |  |
---|---|---|---
neurons | nodes | Error_L2 | Residual
10 | 3 | 1.2825e+00 | 1.7566e+01
| 4 | 1.1336e+00 | 1.4846e+01
| 5 | 1.6084e+00 | 1.5970e+01
| 6 | 6.4758e+00 | 4.7969e+01
| 8 | 1.9045e-01 | 6.7843e+00
| 10 | 1.1092e-02 | 1.5801e-01
20 | 6 | 4.1912e-01 | 2.1314e+01
| 8 | 4.1948e-01 | 7.0691e+00
| 10 | 7.5525e-03 | 4.5484e-01
| 13 | 2.1932e-03 | 4.7934e-02
| 16 | 9.4571e-05 | 4.7279e-03
| 20 | 3.4250e-07 | 2.3290e-05
40 | 13 | 4.5810e-02 | 1.0488e+00
| 16 | 8.2667e-04 | 7.3745e-02
| 20 | 3.6740e-06 | 6.1176e-04
| 26 | 3.6832e-07 | 2.9772e-05
| 33 | 8.2306e-08 | 1.0047e-05
| 40 | 1.9954e-10 | 3.9535e-08
80 | 26 | 1.4968e-04 | 1.1856e-02
| 32 | 1.4893e-05 | 1.9670e-03
| 40 | 2.3953e-07 | 4.7799e-05
| 53 | 1.9385e-10 | 5.7391e-08
| 66 | 3.0444e-10 | 6.2147e-08
| 80 | 2.3971e-11 | 5.9213e-09

k=5 |  |  |
---|---|---|---
neurons | nodes | Error_L2 | Residual
40 | 13 | 1.0291e+00 | 5.4833e+02
| 16 | 9.0266e-01 | 2.8111e+02
| 20 | 8.9477e-01 | 3.7287e+02
| 26 | 1.0524e+00 | 2.7019e+01
| 33 | 1.1745e+00 | 1.1760e+00
| 40 | 3.0679e-01 | 8.3835e+00
80 | 26 | 1.0283e+00 | 2.7422e+01
| 32 | 2.3444e-02 | 1.4636e+00
| 40 | 6.3788e-03 | 2.4238e-01
| 53 | 6.1162e-04 | 3.7449e-02
| 66 | 2.7654e-06 | 1.3766e-03
| 80 | 1.0742e-06 | 1.9191e-04
160 | 53 | 3.2723e-03 | 1.3777e-01
| 64 | 3.0106e-05 | 1.7641e-03
| 80 | 6.4715e-06 | 5.2040e-04
| 106 | 7.7111e-07 | 1.1040e-04
| 133 | 1.1303e-06 | 2.0768e-05
| 160 | 2.8612e-07 | 1.0990e-05
320 | 106 | 5.1808e-06 | 9.5524e-04
| 128 | 6.5518e-06 | 1.1001e-03
| 160 | 5.4403e-07 | 1.7879e-04
| 213 | 1.3397e-07 | 1.0850e-05
| 266 | 9.3901e-08 | 4.1553e-06
| 320 | 2.6074e-09 | 8.7342e-08

Table 1: Error and residuals for the ELM numerical solution of problem (5.1):
on the left table, $k=1$, on the right table $k=5$.

The second boundary value problem that we consider here is a problem
containing a high-order polynomial:

$\left\\{\begin{array}[]{l}u^{\prime\prime}=2^{2p}px(1-x)^{p-2}x^{p-2}(-1+2x-2x^{2}+p(1-4x+4x^{2})),\
0<x<1\\\ u(0)=0\,,\ u(1)=0\end{array}\right.$ (5.2)

The above BV problem has the exact solution:

$u(x)=2^{2p}x^{p}(1-x)^{p}$ (5.3)

Here, we consider the case $p=10$, so that to get a non-exact solution with
FD.

Figure 5: Numerical tests for the solution of problem (5.2) with $p=10$. On
the top panel are depicted the exact (5.3) and the ELM numerical solutions
with $n=40,80$; the number of collocation points $M$ was set to $n/2$. On the
bottom panels are shown the error and residual convergence: on the left panel,
we fix $M=n/2$ and vary $n$ and on the right panel we fix $n=80$ and vary $M$

.

The computed errors arising from the proposed ELM network with respect to the
exact analytical solution are reported in Figure 5.

###  5.2 Boundary Layer problems

In this section, we consider benchmark boundary layer problems, namely an
advection-dominated and a reaction-dominated problem that lead to sharp
gradients near the boundary (see also [37]). We show how the proposed scheme
can deal properly with the steep gradients that arise in these cases.

When $\lambda=0$ and $f\equiv 0$ in (3.1), we have a diffusion-advection
problem. If we impose Dirichlet boundary conditions with $g_{0}=0\,,\
g_{1}=1$, the solution reads:

$u(x)=\dfrac{\text{exp}\left(\dfrac{\gamma}{\mu}x\right)-1}{\text{exp}\left(\dfrac{\gamma}{\mu}\right)-1}\
.$ (5.4)

If $|\frac{\gamma}{\mu}|<<1$, the solution approaches the line connecting the
boundary conditions while if $|\frac{\gamma}{\mu}|>>1$ the solution is near to
zero in almost the whole domain, except in a neighborhood of the right
boundary where the function has a steep gradient. In this last case, we say
that the solution has a boundary layer of width order
$\mathcal{O}\left(\dfrac{\mu}{\gamma}\right)$. From a numerical point of view,
for the particular problem, in which advection dominates diffusion, we need to
catch the behavior on a small scale $\frac{\mu}{\gamma}$, but it is well known
that using finite difference or finite element methods, the approximate
solution can oscillate while the exact solution is monotone (see [37]).

In the advection-dominated problem, the presence of the boundary layer is
usually related to the so-called global Péclet number:

$\mathbb{P}e_{g}=\frac{|\gamma|\cdot|I|}{2\mu}$ (5.5)

where $|I|$ is the domain length.

Figure 6: Numerical tests for the solution of the problem (3.1) with
$\mu=1,\gamma=100,\lambda=0$. In this case, the Péclet number is
$\mathbb{P}e_{g}=50$, the boundary layer is of the order of $10^{-2}$. On the
top panel are shown both the exact (5.4) and the ELM numerical solutions with
$n=40,80$; the number of collocation points is $M=n/2$. On the bottom panels
are shown the error and residual convergence plots; on the left panel we fix
$M=n/2$ and vary $n$ and on the right panel we fix $n=320$ and vary $M$.

We considered the case with $\gamma=100$ and $\mu=1$, with a Péclet number
$\mathbb{P}e_{g}=50$. Figure 5 depicts the exact analytical solution and some
of the ELM numerical solutions along with the corresponding approximation
errors. When $\gamma=0$, $\lambda>0$ and $f\equiv 0$ in the problem (3.1), we
have a diffusion-reaction problem. As before, we impose Dirichlet boundary
conditions with $g_{0}=0\,,\ g_{1}=1$. The problem has analogous difficulties
as the previous one. The exact analytical solution is:

$u(x)=\frac{\text{sinh}(\theta x)}{\text{sinh}(\theta)},\qquad\text{ where
}\theta=\sqrt{\lambda/\mu}.$ (5.6)

Also in this case the solution can give steep gradients if $\lambda/\mu>>1$.
We call this behavior as boundary layer at the right boundary of width order
$\mathcal{O}(\sqrt{\mu/\lambda})$. In this case, the usual definition of
global Péclet number is:

$\mathbb{P}e_{g}=\frac{|\lambda|\cdot|I|^{2}}{6\mu}\ .$ (5.7)

We have set $\lambda=300$ and $\mu=1$ so that the Péclet number is
$\mathbb{P}e_{g}=50$ as before. The computed errors with respect to the exact
solution are reported in Figure 7.

Figure 7: Numerical tests for the solution of the problem (3.1) with
$\mu=1,\gamma=0,\lambda=300$. In this case, the Péclet number is
$\mathbb{P}e_{g}=50$, the boundary layer is of the order of $10^{-1}$. On the
top panel are shown the exact analytical (5.6) and the ELM numerical solutions
for $n=40,80$; the number of collocation points was set to $M=n/2$. On the
bottom panels are depicted the error and the residual convergence: on the left
panel we fix $M=n/2$ and vary $n$ and on the right panel we fix $n=80$ and
vary $M$.

Figure 8: $L^{2}$ errors of computed ELM solutions with fixed $n$ and $M=n/2$.
On abscissae, we report the Péclet number as computed with equation (5.5) on
the left and equation (5.7) on the right. In particular, on the left panel we
consider $\gamma=[0,1,2,5,10,20,50,100,200,500]$ and on the right panel
$\lambda=[0,3,6,15,30,60,150,300,600,1500]$.

Finally, we consider both the advection- and the reaction- dominated problems
taking different $\gamma$ and $\lambda$ so that we get appropriate values of
the Péclet number. In Figure 8, we plot the $L^{2}$ errors of the computed ELM
solutions versus the Péclet numbers for two different choices of number of
neurons.

###  5.3 Internal Layer Problems

In this section, we consider high transient problems that lead to an internal
layer. For the solutions of these problems an adaptive mesh is usually
proposed. For a description of the numerical problems that can arise, we refer
to [33]. In this section, we show that the ELM network can provide good
approximations. First, we consider the $atan$ problem:

$\left\\{\begin{array}[]{l}-u^{\prime\prime}=\frac{(2\alpha^{3}(x-x_{0})}{(1+\alpha^{2}(x-x_{0})^{2})^{2}},\qquad
0<x<1\\\ u(0)=\theta_{0}\,,\ u(1)=\theta_{1}\end{array}\right.$ (5.8)

$\theta_{0},\theta_{1}$ are fixed in order to have an exact solution that
reads:

$u(x)=\text{atan}(\alpha(x-x_{0}))$ (5.9)

In our tests, we fix $\alpha=60$ and $x_{0}=4/9$ that leads to a non-symmetric
internal layer as in [12]. Figure 9 depicts the exact-analytical solution and
some of the numerical ones. The approximation errors with respect to the
exact-analytical solution are also given.

Figure 9:  Numerical tests for the solution of problem (5.8) with $\alpha=60$
and $x_{0}=4/9$. On the top panel it is depicted the exact-analytical (5.9)
and the ELM numerical solutions for $n=640,1280$; the number of collocation
points $M$ was fixed to $n/2$. On the bottom panel, we depict the error and
residual convergence: on the left panel, we fix $M=n/2$ and vary $n$ and on
the right panel we fix $n=1280$ and vary $M$.

The second test with internal peak is a rescaled sinusoidal problem suggested
in [39]:

$\left\\{\begin{array}[]{l}-u^{\prime\prime}=\big{(}\frac{-4x^{2}}{\varepsilon^{2}}+\frac{2}{\varepsilon}\big{)}e^{-x^{2}/\varepsilon},\qquad
0<x<1\\\ u(0)=\theta_{0}\,,\ u(1)=\theta_{1}\end{array}\right.$ (5.10)

$\theta_{0},\theta_{1}$ are fixed in order to have as exact solution:

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

We consider the case $\varepsilon=10^{-3}$. Figure 10 depicts the exact-
analytical solution and some of the numerical ones. The approximation errors
with respect to the exact-analytical solution are also given.

Figure 10: Numerical results for the solution of the problem (5.10) with
$\varepsilon=10^{-3}$. On the top panel are depicted the exact-analytical
(5.11) and the computed ELM solutions with number of neurons $n=640,1280$ and
number of collocation points $M=n/2$. On the bottom panel, are depicted the
error and residual convergence: on the left panel, we fix $M=n/2$ and vary $n$
and on the right panel, we fix $n=1280$ and vary $M$.

To this end, we consider a very difficult numerical problem that has as exact-
analytical solution a comb-like profile. For this problem, we compute the
approximation errors using several number of neurons and the behavior only on
collocation points. The equation we consider is the following:

$\left\\{\begin{array}[]{l}-u^{\prime\prime}=-2(\varepsilon+x)\cos\bigg{(}\frac{1}{\varepsilon+x}\bigg{)}-\frac{\sin\big{(}\frac{1}{\varepsilon+x}\big{)}}{(\varepsilon+x)^{4}},\qquad
0<x<1\\\ u(0)=\theta_{0}\,,\ u(1)=\theta_{1}\end{array}\right.$ (5.12)

We have set $\theta_{0},\theta_{1}$ in order to have the following exact-
analytical solution:

$u(x)=\sin\Biggl{(}\frac{1}{\varepsilon+x}\Biggr{)}\ .$ (5.13)

Then, we considered the case $\varepsilon=1/10\pi$ that gives a solutions with
10 zeros and 5 oscillations in the domain. Figure 11 depicts the exact-
analytical solution as well as the ELM numerical solution for $n=1280,2560$
and the resulting approximation errors with respect to the exact-analytical
solution considering only the approximation errors at the collocation points
using the $L^{\infty}$ norm.

Figure 11: Numerical results for the solution of the problem (5.12) with
$\varepsilon=1/10\pi$. On the top panels it is depicted the exact-analytical
(5.13) and the numerical solutions with the number of neurons $n=1280,2560$
and number of collocation points fixed at $M=n/2$. On the bottom panels are
shown the error and residual convergence: on the left panel, it is shown the
$L^{2}$ norms with fixed $M=n/2$ and on the right panel it is shown the
$L^{\infty}$ on the collocation points.

##  6 Conclusions and future work

In the present paper, we address collocation scheme based on ELM networks for
the numerical approximation of 1D elliptic problems with steep gradients. The
weights between the input layer and the hidden layer as well as the biases are
fixed appropriately in advance and kept fixed. In this way, the only
coefficients that describe the approximation space are the weights of the
linear combination used for the element of the linear activation of the ELM.
First, we have explored the underlying space, so as to construct it in a
random but proper way. Then, the numerical solution to the boundary value
problems was reduced to the computation of the $n$ weights of the ELM network
via $M$ collocation equations. We have proved the convergence of the scheme in
the square case, and proposed the use of an underdetermined collocation scheme
to avoid ill-conditioning.

Numerical tests were performed for the evaluation of the approximation
accuracy of the proposed method with the aid of benchmark boundary value
problems exhibiting sharp gradients. We also compared the proposed approach
with a high order finite difference (FD) scheme. We show that for both
internal and boundary layer problems the proposed method outperforms FD, when
the number of neurons is taken to be big enough to catch the steep gradient of
the layer. We emphasize also that in our case, the solutions are $C^{\infty}$
functions.

The proposed method can be further developed by exploring some alternatives,
such as the use of optimal activation functions, as proposed in [14], or the
selection of the solution of the under-determined system by using the
properties of null rules, as explored in [7]. Finally, the application of the
proposed ELM collocation method can be extended in various ways. Among others,
we are currently considering the use for:

  * 1.

Non linear differential problems;

  * 2.

Multidimensional problems;

  * 3.

Time-dependent problems;

  * 4.

Inverse problems with overfitting.

## Acknowledgements

The authors would like to thank Prof. Ferdinando Auricchio, Prof. Bert Jüttler
for fruitful discussions on the subject of the paper. Francesco Calabrò are
partially supported by INdAM, through GNCS research projects. This support is
gratefully acknowledged.

## References

  * [1] J. Almira, P. Lopez-de Teruel, D. Romero-Lopez, and F. Voigtlaender.  Negative results for approximation using single layer and multilayer feedforward neural networks.  Journal of Mathematical Analysis and Applications, page 124584, 2020\.
  * [2] C. Anitescu, E. Atroshchenko, N. Alajlan, and T. Rabczuk.  Artificial neural network methods for the solution of second order boundary value problems.  Computers, Materials and Continua, 59(1):345–359, 2019.
  * [3] H. Arbabi, J. E. Bunder, G. Samaey, A. J. Roberts, and I. G. Kevrekidis.  Linking machine learning with multiscale numerics: Data-driven discovery of homogenized equations.  JOM, pages 1–14, 2020.
  * [4] J.-G. Attali and G. Pagès.  Approximations of functions by a multilayer perceptron: a new approach.  Neural networks, 10(6):1069–1081, 1997.
  * [5] F. Auricchio, L. B. Da Veiga, T. J. Hughes, A. Reali, and G. Sangalli.  Isogeometric collocation for elastostatics and explicit dynamics.  Computer methods in applied mechanics and engineering, 249:2–14, 2012.
  * [6] J. Berg and K. Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  Neurocomputing, 317:28–41, 2018.
  * [7] F. Calabrò, D. Bravo, C. Carissimo, E. Di Fazio, A. Di Pasquale, A. Eldray, C. Fabrizi, J. Gerges, S. Palazzo, and J. Wassef.  Null rules for the detection of lower regularity of functions.  Journal of Computational and Applied Mathematics, 361:547–553, 2019\.
  * [8] Q. Chan-Wai-Nam, J. Mikael, and X. Warin.  Machine learning for semi linear pdes.  Journal of Scientific Computing, 79(3):1667–1712, 2019.
  * [9] D. Costarelli and R. Spigler.  Approximation results for neural network operators activated by sigmoidal functions.  Neural Networks, 44:101–106, 2013.
  * [10] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314, 1989\.
  * [11] C. De Falco and E. O’Riordan.  Interior layers in a reaction-diffusion equation with a discontinuous diffusion coefficient.  Int. J. Numer. Anal. Model, 7(3):444–461, 2010.
  * [12] L. Demkowicz, W. Rachowicz, and P. Devloo.  A fully automatic hp-adaptivity.  Journal of Scientific Computing, 17(1-4):117–142, 2002.
  * [13] V. Dwivedi and B. Srinivasan.  Physics informed extreme learning machine (PIELM) - A rapid method for the numerical solution of partial differential equations.  Neurocomputing, 391:96 – 118, 2020.
  * [14] N. J. Guliyev and V. E. Ismailov.  On the approximation by single hidden layer feedforward neural networks with fixed weights.  Neural Networks, 98:296–304, 2018.
  * [15] H. Guo, X. Zhuang, X. Meng, and T. Rabczuk.  Analysis of three dimensional potential problems in non-homogeneous media with deep learning based collocation method.  arXiv preprint arXiv:2010.12060, 2020.
  * [16] N. Hahm and B. I. Hong.  An approximation by neural networks with a fixed weight.  Computers & Mathematics with Applications, 47(12):1897–1903, 2004\.
  * [17] J. Han, A. Jentzen, and E. Weinan.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [18] J. Han, M. Nica, and A. R. Stinchcombe.  A derivative-free method for solving elliptic partial differential equations with deep neural networks.  Journal of Computational Physics, 419:109672, 2020.
  * [19] K. Hornik, M. Stinchcombe, and H. White.  Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks.  Neural networks, 3(5):551–560, 1990.
  * [20] K. Hornik, M. Stinchcombe, H. White, et al.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [21] H.-Y. Hu and Z.-C. Li.  Collocation methods for Poisson’s equation.  Computer Methods in Applied Mechanics and Engineering, 195(33-36):4139–4160, 2006.
  * [22] G. Huang, G.-B. Huang, S. Song, and K. You.  Trends in extreme learning machines: A review.  Neural Networks, 61:32 – 48, 2015.
  * [23] G.-B. Huang, D. H. Wang, and Y. Lan.  Extreme learning machines: a survey.  International journal of machine learning and cybernetics, 2(2):107–122, 2011.
  * [24] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew.  Extreme learning machine: a new learning scheme of feedforward neural networks.  In 2004 IEEE international joint conference on neural networks (IEEE Cat. No. 04CH37541), volume 2, pages 985–990. IEEE, 2004.
  * [25] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew.  Extreme learning machine: theory and applications.  Neurocomputing, 70(1-3):489–501, 2006.
  * [26] Y. Ito.  Nonlinearity creates linear independence.  Advances in Computational Mathematics, 5(1):189–203, 1996.
  * [27] I. E. Lagaris, A. Likas, and D. I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9(5):987–1000, 1998.
  * [28] I. E. Lagaris, A. C. Likas, and D. G. Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks, 11(5):1041–1049, 2000.
  * [29] S. Lee, M. Kooshkbaghi, K. Spiliotis, C. I. Siettos, and I. G. Kevrekidis.  Coarse-scale pdes from fine-scale observations via machine learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(1):013141, 2020.
  * [30] S. Lin, X. Guo, F. Cao, and Z. Xu.  Approximation by neural networks with scattered data.  Applied Mathematics and Computation, 224:29–35, 2013.
  * [31] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  Deepxde: A deep learning library for solving differential equations.  arXiv preprint arXiv:1907.04502, 2019.
  * [32] C. Michoski, M. Milosavljević, T. Oliver, and D. Hatch.  Solving differential equations using deep neural networks.  Neurocomputing, 2020.
  * [33] W. F. Mitchell.  A collection of 2D elliptic problems for testing adaptive grid refinement algorithms.  Applied Mathematics and Computation, 220:350 – 364, 2013.
  * [34] E. Oñate, J. Miquel, and P. Nadukandi.  An accurate FIC-FEM formulation for the 1D advection–diffusion–reaction equation.  Computer Methods in Applied Mechanics and Engineering, 298:373–406, 2016.
  * [35] G. Pang, L. Yang, and G. E. Karniadakis.  Neural-net-induced gaussian process regression for function approximation and pde solution.  Journal of Computational Physics, 384:270–288, 2019.
  * [36] A. Pinkus.  Approximation theory of the MLP model in neural networks.  Acta numerica, 8(1):143–195, 1999.
  * [37] A. Quarteroni.  Diffusion-transport-reaction equations.  In Numerical Models for Differential Problems, pages 315–365. Springer, 2017.
  * [38] K. Rudd and S. Ferrari.  A constrained integration (CINT) approach to solving partial differential equations using artificial neural networks.  Neurocomputing, 155:277–285, 2015.
  * [39] A. Schmidt and K. G. Siebert.  A posteriori estimators for the h–p version of the finite element method in 1d.  Applied numerical mathematics, 35(1):43–66, 2000.
  * [40] J. Sirignano, J. F. MacArt, and J. B. Freund.  DPM: A deep learning PDE augmentation method with application to large-eddy simulation.  Journal of Computational Physics, 423:109811, 2020.
  * [41] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [42] G. Strang and G. J. Fix.  An analysis of the finite element method.  Prentice-hall, 1973.


# Free Energy Calculations using Smooth Basin Classification

Sander Vandenhaute  Center for Molecular Modeling, Ghent University (BE)  Tom
Braeckevelt  Center for Molecular Modeling, Ghent University (BE)  Pieter
Dobbelaere  Center for Molecular Modeling, Ghent University (BE)  Massimo
Bocus  Center for Molecular Modeling, Ghent University (BE)  Veronique Van
Speybroeck<EMAIL_ADDRESS>Center for Molecular Modeling,
Ghent University (BE)

###### Abstract

The efficiency of atomic simulations of materials and molecules can rapidly
deteriorate when large free energy barriers exist between local minima. We
propose smooth basin classification, a universal method to define reaction
coordinates based on the internal feature representation of a graph neural
network. We achieve high data efficiency by exploiting their built-in symmetry
and adopting a transfer learning strategy. We benchmark our approach on
challenging chemical and physical transformations, and show that it matches
and even outperforms reaction coordinates defined based on human intuition.

Dynamic simulations of materials and molecules can provide both qualitative
and quantitative information on chemical and physical transformations. While
explicit integration of Newton’s equation of motion requires timesteps in the
order of a femtosecond, the time required for many transformations to complete
is often beyond the micro- or millisecond range because of large energy
barriers between initial and final state. As such, simulation of a single
transformation can easily require billions of time steps [1, 2].

Over the past decades, a variety of advanced sampling techniques have been
proposed in order to improve the sampling efficiency of dynamic simulations. A
particularly intuitive approach is to bias the dynamics along certain
directions in phase space, effectively lifting the system over obstructive
free energy barriers. The bias energy is defined along a certain predefined
reaction coordinate or collective variable (CV). The CV is a function of the
atomic coordinates and should naturally discriminate between initial and final
state without violating the physical symmetries in the system. Specifically,
the CV should remain invariant with respect to global translations and
rotations, as well as permutations of atoms which are chemically identical. In
the vast majority of practical applications, such CVs are manually defined
based on chemical and physical intuition, and they usually contain simple
invariants such as specific unit cell parameters or coordination numbers of
one chemical species with respect to another. However, manual definition can
be highly challenging and time-consuming, because the metastable
intermediate(s) and/or transition state(s) are not always known in advance, or
because the required invariances are difficult to satisfy [3, 4, 5].

Considering the numerous achievements of machine learning in atomistic
modeling, it is sensible to anticipate that suitable reaction coordinates can
also be learned from (unbiased) simulation data, and ideally without requiring
human intuition about the transition mechanism(s). A wide variety of such
methods have been proposed over the past decade, and we briefly mention a few
notable examples while referring the reader to specialized reviews for a more
comprehensive discussion [6, 7]. Supervised methods learn a continuous
discriminant function which separates free energy minima based on user-defined
structural descriptors, for example using (deep) linear discriminant analysis
[8, 9, 10]. Unsupervised methods do not directly learn a discriminant function
but instead identify so-called slow modes, either within a given set of
descriptors or entirely end-to-end; see e.g. time-lagged independent component
analysis or time-lagged autoencoders [11, 12, 13]. While these methods
demonstrate the potential of using machine learning for defining reaction
coordinates, their widespread adoption has remained rather limited. We
hypothesize that this is due to inherent limitations regarding the required
volume of simulation data, the difficulty of constructing an appropriate set
of input features, or the difficulty in interpreting and validating the
obtained CVs.

Figure 1:  (a) Illustration of the partitioning of the phase space into basins
of attraction. Each point $\bm{x}$ can be uniquely associated to a basin via
its associated potential energy minimum. (b) Base GNN which is pretrained to
reference energies and forces. The interaction layers effectively generate a
learned feature representation $\bm{f}$ for the atomic environments. (c) The
GNN is augmented with a second readout module; a phase classifier which is
trained separately on atomic geometries which are labeled with the
corresponding basin as defined in (a). (d) Schematic illustration of the
feature space embedding of atomic environments based on their phase, together
with the isolines of $\Delta q$. (e) Active learning procedure for the phase
readout layers; the log probabilities $\bm{q}$ are used to bias the dynamics
from one basin to another, in order to sample geometries along the transition
path(s). The gathered data is labeled with the correct basin of attraction (as
determined by minimizing the energy) after which the readout is retrained. See
Section S2 in the Supplemental Material for more information.

In this work, we propose smooth basin classification (SBC); a universal method
to construct CVs and compute free energy differences for chemical and physical
transformations. SBC builds upon the successful development of graph neural
networks (GNNs) as effective interatomic potentials by using their learned
feature space as ansatz for constructing physically meaningful CVs [14, 15,
16]. We apply SBC for complex rare events in both materials science and
chemistry, and demonstrate that it matches or outperforms CVs based on human
intuition. Relative free energy predictions based on learned CVs are within 1
kJ/mol of reference values as computed based on established CVs for each of
the case studies. In addition, we observe that SBC requires orders of
magnitude less training data as compared to existing approaches [6, 7, 8, 9,
10, 11, 12, 13]. This is a direct consequence of the fact that the constructed
CVs are embedded in the learned feature space of the GNN, thereby inheriting
all of their established benefits regarding data efficiency and universality.

Consider a system with two basins of attraction, named A and B. Each point
$\bm{r}$ in the system’s coordinate space can be uniquely labeled with A or B
depending on the trajectory of a geometry optimization; $\bm{r}$ belongs to
basin A whenever it optimizes to any of the local minima within A, and vice-
versa for B (Figure 1, a). In SBC, we aim to construct a mapping between each
configuration $\bm{r}$ and its associated basin of attraction. Instead of
defining this mapping directly in terms of the atomic coordinates – which are
not invariant with respect to the symmetry operations in the system – we embed
it into the learned feature space $\bm{f}$ of a GNN which has been pretrained
to reproduce the potential energy surface of the system (Figure 1, b). In
general, the interaction layers in the GNN parse the chemical environment of
atom $i$ into a feature vector $\bm{f}_{i}$ in this space (which typically
contains anywhere between 8 and 64 dimensions). Importantly, this feature
vector is rigorously invariant with respect to all symmetry operations in the
system, and is specifically trained to capture the physical interactions in
the system. Both the atomic energy $E_{i}$ and the total energy $E$ are
computed from the atomic feature vectors $\bm{f}_{i}$ using a simple readout
function:

$\displaystyle\bm{r}\quad\xrightarrow[]{\text{GNN}}\quad\\{\bm{f}_{i}\\}_{i=1}^{n}\quad\xrightarrow[]{\text{energy
readout}}\quad E=\sum_{i=1}^{n}E_{i}$ (1)

Because the atomic features $\bm{f}_{i}$ are maximally informative for
determining the interaction energy of a given atomic environment, we may
assume that they are also highly informative in predicting the basin of
attraction, i.e. the final destination of an energy minimization. In full
analogy with equation 1, we can therefore construct a prediction for the basin
of attraction using the same GNN but with a different readout function (Figure
1, c):

$\displaystyle\bm{r}\quad\xrightarrow[]{\text{GNN}}\quad\\{\bm{f}_{i}\\}_{i=1}^{n}\quad\xrightarrow[]{\text{phase
readout}}\quad\bm{q}=\sum_{i=1}^{n}\bm{q}_{i}$ (2)

in which we introduce the per-atom log probabilities $\bm{q}_{i}$ as an
invariant tuple which measures how strongly the model associates a given
atomic environment $\bm{f}_{i}$ to basin A or B – as commonly done in a
classification setting [17]:

$\displaystyle\bm{q}_{i}=\Big{(}q_{i}^{A},q_{i}^{B}\Big{)}$ (3)

For example, if the model believes that a given atomic environment $i$ is more
likely to appear in phase $A$ as opposed to phase $B$, it will predict log
probabilities $\bm{q}_{i}$ such that $q_{i}^{A}>q_{i}^{B}$. After summing the
per-atom log probabilities $\bm{q}_{i}$ to obtain the total log probabilities
$\bm{q}=\left(q^{A},q^{B}\right)$, the predicted basin of attraction
$\mathcal{C}^{*}$ of the total configuration is then determined by the phase
which received the largest log probability in $\bm{q}$:

$\displaystyle\mathcal{C}^{*}=\operatorname*{argmax}_{\mathcal{C}}\bm{q}$ (4)

A more detailed overview of the model architecture can be found in the
Supplemental Material (Section S1).

Figure 2: (a) The nonreactive phase transition in MIL-53(Al). The scatter plot
shows the distribution of both the unit cell volume $V$ as well as the log
probability difference $\Delta q=q_{\text{op}}-q_{\text{cp}}$. Samples were
obtained from enhanced sampling simulations with umbrellas centered on $\Delta
q$. Each dot represents a single atomic geometry, and its color represents
whether it minimizes to the $op$ or $cp$ phases. Black dots indicate potential
energy minima. (b) Free energy profiles as obtained using umbrella sampling
and multistate Bennett acceptance ratio (MBAR) analysis [18, 19, 20, 21],
using the unit cell volume (bottom) or the log probability difference (top) as
biasing coordinate. While free energy profiles are intrinsically dependent on
the mathematical form of the CV, the total relative free energy difference is
intrinsic to the system and can be obtained by integrating over the minima.
Both approaches correctly yield the same value (within 1 kJ/mol). See the
Supplemental Material (Section S5) for additional computational details.

While the prediction of the basin of attraction is essentially a discrete
classification task, the actual log probabilities $\bm{q}_{i}$ are smooth
functions of the atomic feature vectors $\bm{f}_{i}$ and hence the atomic
coordinates – which means that they can be used to introduce biasing forces
and/or virial stress. Ideally, we would like to design this bias such that it
can drive transitions from and to any of the basins. To this end, consider the
difference between two log probabilities:

$\displaystyle\Delta q$ $\displaystyle=q^{B}-q^{A}$ (5)

In essence, $\Delta q$ measures how much more likely the model thinks a given
configuration belongs to basin B rather than basin A. It is therefore strictly
positive in basin B and strictly negative in basin A provided that the phase
readout is properly trained (Figure 1, d). Clearly, the quantity $\Delta q$
correctly discriminates between two free energy minima while also being
continuously differentiable with respect to the atomic coordinates (using
standard backpropagation). It is therefore, by definition, a valid collective
variable for the transition between basins A and B. Note that the invariance
of $\Delta q$ with respect to global translations and rotations as well as
permutation of chemically identical atoms is guaranteed by construction of the
GNN. The built-in symmetry of SBC towards translations, rotations, and
permutations incorporates important physical behavior into the model, and
should be contrasted with both common practice CVs and most other ML-based
approaches [9, 13, 12, 22]. which often use functions and/or input features
that depend explicitly on the atomic indices (thereby violating permutation
invariance) or specific components of the unit cell vectors (thereby violating
rotational as well as periodic invariance).

The phase readout layers are analogous to the existing energy readouts found
in conventional GNNs [14, 15, 16] (Figure 1, c), and essentially consist of a
multilayer perceptron (MLP) containing at most two nonlinearities and about
$\sim$1000 weights in total (compared to $\sim$100,000 weights in the GNN).
The classifier is trained while keeping the GNN and energy readouts fixed,
similar to conventional transfer learning approaches. Its training data
consists of atomic geometries which are labeled with their corresponding basin
of attraction (e.g. A or B). To ensure that the classifier remains accurate
for out-of-equilibrium geometries, we employ a combination of enhanced
sampling (in $\Delta q$) and active learning to explore transition regions
between phases in an iterative fashion (Figure 1, e). More information
regarding model training and the active learning approach can be found in the
Supplemental Material (Sections S1 and S2).

Figure 3: (a) A proton transfer reaction in a zeolite. The scatter plot shows
the distribution of both the coordination number difference as well as the log
probability difference $\Delta q=q_{\text{cat}}-q_{\text{iso}}$. Samples were
obtained from enhanced sampling simulations with umbrellas centered on $\Delta
q$. Black dots indicate potential energy minima; red dots indicate isobutene
configurations; blue dots indicate cation configurations. Remarkably, the
learned CV even outperforms the coordination number. (b) Free energy profiles
as obtained using umbrella sampling and MBAR analysis, using the coordination
number (bottom) or $\Delta q$ (top) as biasing coordinate. Free energy
differences were computed by integrating $F$ over the respective minima. See
the Supplemental Material (Section S5) for additional computational details.

Figure 2 (a) shows a nonreactive phase transition in a solid-state system
named MIL-53(Al), a flexible metal-organic framework based on a winerack
topology [23]. The system exhibits two basins of attraction which are relevant
for tuning the material towards gas storage and/or sensing applications; a
closed pore ($cp$) and open pore ($op$) phase. In these types of transitions,
i.e. those for which the unit cell parameters vary significantly between the
different basins, researchers often employ some combination of unit cell
parameters as collective variable – in this case, the unit cell volume [24].
Figure 2 (a) demonstrates the performance of the unit cell volume as formal
discriminator between the two basins, and compares it with a difference of
learned log probabilities as proposed in equation 5. Both criteria perform
excellent in terms of classification accuracy, and in addition exhibit a very
strong correlation. For both collective variables, the misclassified
geometries are all localized near the transition region (such they will not
cause substantial differences in the relative stability). Note that the
network learns to distinguish low- and high-volume configurations based on
atomic environments only, i.e. without explicitly relying on unit cell
parameters as input.

Using the log probability difference $\Delta q=q_{\text{op}}-q_{\text{cp}}$ as
CV (see equation 5), we can employ conventional enhanced sampling techniques
and evaluate the free energy profile $F(\Delta q)$. Here, we choose to employ
umbrella sampling with harmonic restraints on the CV:

$\displaystyle U_{\text{bias}}(\bm{r},\bm{h})=\frac{K}{2}\Big{[}\Delta
q\left(\bm{r},\bm{h}\right)-q_{0}\Big{]}^{2}$ (6)

in which we emphasize that the bias energy contains an implicit dependence on
both the atomic coordinates $\bm{r}$ and the unit cell parameters $\bm{h}$
through the learned atomic environment features $\bm{f}_{i}$ [25]. We use the
multistate Bennett acceptance ratio (MBAR) equations to combine samples from
multiple umbrellas in a statistically optimal manner and obtain free energy
profiles as shown in Figure 2 (b, top) [18, 19, 20, 21]. The reference free
energy profile using the unit cell volume as bias coordinate is added for
completeness in Figure 2 (b, bottom). The total relative free energy
difference between the $op$ and $cp$ basin, as obtained by integrating each of
the profiles in their respective minima, agrees to within 1 kJ/mol between
both CVs. Note that while relative stabilities are intrinsic to the system,
the actual free energy profiles depend on the particular choice of CV, so the
learned and reference free energy profiles need not have the same shape
(although there exist ways to transform one profile into another [26]).

For global phase transformations in solids such as the one in Figure 2, a
size-extensive definition of $\bm{q}$ as a sum over per-atom log probabilities
$\bm{q}_{i}$ is intuitively sensible because it gives the model the ability to
assign nonzero log probability contributions to all atoms simultaneously.
However, in chemical transformations, the rare event is most often a local
process whose free energy barrier does not necessarily depend on the size of
the surrounding environment (e.g. the catalyst). To demonstrate the
applicability of our approach to those cases, we investigate a proton transfer
reaction from a Brønsted acid site in a zeolite to an alkene guest (Figure 3,
a). The reference collective variable in this case is a coordination number
difference which is defined explicitly in terms of the indices of the oxygen
atoms of the active site as well as the indices of the carbon atoms of the
molecule and all of the hydrogens – see the Supplemental Material (Section S5)
for more details. While such definitions are extremely common, they violate
permutation symmetry and make each collective variable highly system-specific.
Figure 3 (a) shows the classification performance of the coordination number
versus a learned CV based on equation 5. Again, classification accuracies are
excellent for both the proposed CV based on chemical intuition as well as the
learned log probability difference. Misclassified geometries are localized at
or near the transition region, which means that they will not contribute to
the relative free energy of both phases. Indeed, as indicated in Figure 3 (b),
the free energy difference between both basins as computed from either of the
two profiles agrees to within 1 kJ/mol.

Figure 4: Free energy landscape of alanine dipeptide at the PBE-D3(BJ) level
of theory, as obtained from umbrella sampling in the conventional backbone
dihedral angles $\phi$ and $\psi)$. The scattered dots indicate the training
data distribution as obtained during active learning of the SBC model; the
color-coding indicates the log probability difference (CV value) which was
assigned to them after final training of the log probabilities. See the
Supplemental Material (Section S5) for additional computational details.

A final case study concerns the conformational free energy landscape of
alanine dipeptide, given that it is the most commonly studied system in the
context of free energy calculations and the automated development of
collective variables [12, 27, 13, 28, 29, 30, 9, 31, 32]. At the PBE-D3(BJ)
level of theory and room temperature conditions, the molecule is observed to
exhibit two main stable conformations; the c7eq and c7ax conformations. There
exist multiple transition paths between both conformations, and each comprises
multiple potential energy minima. The free energy of the system is most often
visualized using a two-dimensional free energy landscape of the molecule based
on two conventional backbone dihedral angles $\phi$ and $\psi$. Figure 4
presents the free energy landscape as computed using traditional $(\phi,\psi)$
umbrella sampling, and superimposes a scatter plot of the samples gathered
during the active learning loop of SBC. Note that the active learning loop is
bootstrapped from equilibrium samples only, i.e. the model does not a priori
know what the transition regions look like. Nevertheless, iterative
application of harmonic bias potentials centered on the learned log
probabilities (see Figure 1) yields a complete sampling of the phace space. In
the Supplemental Material (Section S5), we further demonstrate that the
learned log probabilities can be used to perform 1D umbrella sampling in order
to estimate the relative stability between c7eq and c7ax, which is much more
efficient than having to perform 2D umbrella sampling across the entire
($\phi,\psi$) domain.

By embedding SBC into the learned feature space of a GNN, we exploit the
intrinsic overlap that exists between a quantitative understanding of atomic
interactions and free energy minima that arise as a result of those
interactions. To demonstrate the efficiency of this deliberate design choice,
Table 1 enumerates and compares the most prominent CV learning methods with
SBC for alanine dipeptide. We evaluate a number of key characteristics,
including the required amount of input data to train a given CV, its
transferability towards symmetrically equivalent geometries, whether or not it
requires manual feature selection, and whether or not an initial transition
path is required. In essence, Table 1 shows that SBC outperforms every single
method in all of these aspects. It requires the fewest structures and smallest
total simulation time; it is the only method that is rigorously invariant with
respect to all symmetry operations; and it does not require any kind of manual
input feature selection or initial transition path. While not highlighted in
Table 1, it should also be noted that SBC is competitive in terms of its
computational cost during inference, i.e. during production (enhanced
sampling) simulations. Evaluation of the log probability difference $\Delta q$
as a learned CV during biased simulations is essentially trivial since the
node features $\bm{f}_{i}$ have already been computed and used for the
(unbiased) energy and forces.

| training cost |  input descriptors | symmetries |  initialized without prior transition path(s)
---|---|---|---|---
| snapshots | time [ns] |  | trans | rot | perm |
time-lagged AE [12] |  | 750 ns | learned111based on alignment with respect to a reference atomic geometry | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
FEBILAE [27] |  | 1500 ns | learned1 | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
VAMPnets [13] | 250,000 | 250 ns | learned1 | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | ✓
LINES [28] |  | 100 ns | manual | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
LED [29] |  | 100 ns | manual | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
Ref. [30] |  $>$300,000 |  | learned1 | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | ✓
Deep-LDA [9] | 10,000 |  | manual | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
Deep-TDA [22] |  | 40 ns | learned222based on a list of interatomic distances | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | ✓
Ref. [31] | 15,000 | 15 ns | manual | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | ✓
Ref. [32] | 8,000 | 512 ns | manual | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ |
SBC | 3,694 | 11 ns | learned | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | $\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$ | ✓
Table 1: Performance of various CV learning methods on alanine dipeptide. Key
characteristics are the required amount of input data (lower is better), the
built-in symmetry as inductive bias, and whether or not an initial transition
path is required for the model to learn the transition. Missing values
indicate that the corresponding quantity was not directly reported in either
the main manuscript or any of its appendices. For each model, we indicate
whether it rigorously satisfies the translational, rotational, and
permutational symmetry of the underlying potential energy surface. Yellow dots
indicate that the required invariance was achieved by means of a manual
structural alignment procedure – see the corresponding references for more
information.

We conclude with a discussion on the current limitations of SBC as well as
directions for future research. Clearly, the success of our method is
critically dependent on the availability of ML potentials which can predict
the atomic interactions in the system with high fidelity. With the
introduction of equivariant feature representations, specialized interaction
layers, and efficient online learning workflows, such potentials can now be
routinely developed without requiring massive amounts of QM reference data,
even for challenging systems [16, 33, 34, 35]. In addition, the availability
of foundation models – which are applicable to essentially the entire periodic
table – can alleviate some of the burden that comes with training and
validating the base GNN [36, 37, 38]. In this work, we considered both self-
developed and foundation models as base GNN in which the SBC model is inserted
– the alanine dipeptide and MIL-53(Al) case studies were performed with custom
ML potentials trained using psiflow (see the Supplemental Material, Section
S4), while a MACE foundation model was used for the proton transfer reaction.
As such, we can conclude that SBC is fundamentally compatible with both.
However, more challenging scenarios might require multi-task training of the
entire GNN (including its interaction layers), in which each sample in the
training set is labeled with QM reference data and additionally a phase or
basin label – this will incentivize the model to design the per-atom feature
embedding to be informative for both potential energy prediction as well as
basin classification. Finally, we mention that the proposed method can be
extended towards the more general case of $n$ free energy minima simply by
increasing the number of classes in the classification layer. This will be
necessary when considering more extreme scenarios such as large chemical
reaction networks or protein conformational landscapes.

###### Acknowledgements.

The authors acknowledge funding from the Research Board of Ghent University
(V.V.S) and iBOF-21-085 PERSIST (V.V.S. and T.B.). S.V. and P.D. wish to thank
the Research Foundation – Flanders (FWO) for doctoral fellowships (grant nos.
11H6821N and 11O2123N respectively). M.B. acknowledges financial support from
the Fund for Scientific Research Flanders and the Excellence of Science (EOS)
Project BioFact (EOS ID 30902231). We acknowledge the EuroHPC Joint
Undertaking for awarding this project access to the EuroHPC supercomputer
LUMI, hosted by CSC (Finland) and the LUMI consortium through a EuroHPC
Regular Access call. Part of the simulations were performed on the Luxembourg
national supercomputer MeluXina. The authors gratefully acknowledge the
LuxProvide teams for their expert support.

## References

  * Frenkel and Smit [2023] D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_ (Elsevier Science, 2023).
  * Tuckerman [2010] M. Tuckerman,  _Statistical Mechanics: Theory and Molecular Simulation_ , Oxford Graduate Texts (OUP Oxford, 2010).
  * Grifoni  _et al._ [2019] E. Grifoni, G.-M. Piccini, and M. Parrinello, Microscopic description of acid–base equilibrium, Proc. Natl. Acad. Sci. U.S.A. 116, 4054 (2019).
  * Ahalawat and Mondal [2018] N. Ahalawat and J. Mondal, Assessment and optimization of collective variables for protein conformational landscape: GB1 $\beta$-hairpin as a case study, J. Chem. Phys. 149, 094101 (2018).
  * Niu  _et al._ [2018] H. Niu, P. M. Piaggi, M. Invernizzi, and M. Parrinello, Molecular dynamics simulations of liquid silica crystallization, Proc. Natl. Acad. Sci. U.S.A. 115, 5348 (2018).
  * Chen [2021] M. Chen, Collective variable-based enhanced sampling and machine learning, Eur. Phys. J. B 94, 211 (2021).
  * Bhakat [2022] S. Bhakat, Collective variable discovery in the age of machine learning: reality, hype and everything in between, RSC Adv. 12, 25010 (2022).
  * Elishav  _et al._ [2023] O. Elishav, R. Podgaetsky, O. Meikler, and B. Hirshberg, Collective variables for conformational polymorphism in molecular crystals, J. Phys. Chem. Lett. 14, 971 (2023).
  * Bonati  _et al._ [2020] L. Bonati, V. Rizzi, and M. Parrinello, Data-driven collective variables for enhanced sampling, J. Phys. Chem. Lett. 11, 2998 (2020).
  * Mendels  _et al._ [2018] D. Mendels, G. Piccini, and M. Parrinello, Collective variables from local fluctuations, J. Phys. Chem. Lett. 9, 2776 (2018).
  * Pérez-Hernández  _et al._ [2013] G. Pérez-Hernández, F. Paul, T. Giorgino, G. De Fabritiis, and F. Noé, Identification of slow molecular order parameters for markov model construction, J. Chem. Phys. 139, 015102 (2013).
  * Wehmeyer and Noé [2018] C. Wehmeyer and F. Noé, Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics, J. Chem. Phys. 148, 241703 (2018).
  * Mardt  _et al._ [2018] A. Mardt, L. Pasquali, H. Wu, and F. Noé, Vampnets for deep learning of molecular kinetics, Nat. Commun. 9, 5 (2018).
  * Thomas  _et al._ [2018] N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, Tensor field networks: Rotation- and translation-equivariant neural networks for 3d point clouds (2018), arXiv:1802.08219 [cs.LG] .
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun. 13, 2453 (2022).
  * Kovács  _et al._ [2023] D. P. Kovács, I. Batatia, E. S. Arany, and G. Csányi, Evaluation of the MACE force field architecture: From medicinal chemistry to materials science, J. Chem. Phys. 159, 044118 (2023).
  * Bishop [2006] C. M. Bishop,  _Pattern Recognition and Machine Learning (Information Science and Statistics)_ (Springer-Verlag, Berlin, Heidelberg, 2006).
  * Shirts and Chodera [2008] M. R. Shirts and J. D. Chodera, Statistically optimal analysis of samples from multiple equilibrium states, J. Chem. Phys. 129, 124105 (2008).
  * Shirts [2017] M. R. Shirts, Reweighting from the mixture distribution as a better way to describe the multistate bennett acceptance ratio (2017), arXiv:1704.00891 [cond-mat.stat-mech] .
  * Shirts and Ferguson [2020] M. R. Shirts and A. L. Ferguson, Statistically optimal continuous free energy surfaces from biased simulations and multistate reweighting, J. Chem. Theory Comput. 16, 4107 (2020).
  * Shirts  _et al._ [2023] M. Shirts, K. Beauchamp, L. Naden, J. Chodera, jaimergp, S. Martiniani, C. Stern, M. Henry, J. Fass, R. Gowers, R. T. McGibbon, B. Dice, C. Jones, D. L. Dotson, F. York, I. Pulido, M. Thompson, and T. Burgin, choderalab/pymbar: 4.0.2 optional jax build (2023).
  * Trizio and Parrinello [2021] E. Trizio and M. Parrinello, From enhanced sampling to reaction profiles, The Journal of Physical Chemistry Letters 12, 8621 (2021).
  * Serre  _et al._ [2002] C. Serre, F. Millange, C. Thouvenot, M. Noguès, G. Marsolier, D. Louër, and G. Férey, Very large breathing effect in the first nanoporous chromium(iii)-based solids: MIL-53 or CrIII(OH)·{O2C-C6H4-CO2}·{HO2C-C6H4-CO2H}x·H2Oy, Journal of the American Chemical Society 124, 13519 (2002).
  * Demuynck  _et al._ [2018] R. Demuynck, J. Wieme, S. M. J. Rogge, K. D. Dedecker, L. Vanduyfhuys, M. Waroquier, and V. Van Speybroeck, Protocol for identifying accurate collective variables in enhanced molecular dynamics simulations for the description of structural transformations in flexible metal–organic frameworks, Journal of Chemical Theory and Computation 14, 5511 (2018).
  * Tan  _et al._ [2024] A. R. Tan, J. C. B. Dietschreit, and R. Gomez-Bombarelli, Enhanced sampling of robust molecular datasets with uncertainty-based collective variables (2024), arXiv:2402.03753 [cs.LG] .
  * Bailleul  _et al._ [2020] S. Bailleul, K. Dedecker, P. Cnudde, L. Vanduyfhuys, M. Waroquier, and V. Van Speybroeck, Ab initio enhanced sampling kinetic study on mto ethene methylation reaction, Journal of Catalysis 388, 38 (2020).
  * Belkacemi  _et al._ [2022] Z. Belkacemi, P. Gkeka, T. Lelièvre, and G. Stoltz, Chasing collective variables using autoencoders and biased trajectories, J. Chem. Theory Comput. 18, 59 (2022).
  * Odstrcil  _et al._ [2022] R. E. Odstrcil, P. Dutta, and J. Liu, Lines: Log-probability estimation via invertible neural networks for enhanced sampling, J. Chem. Theory Comput. 18, 6297 (2022).
  * Vlachas  _et al._ [2022] P. R. Vlachas, J. Zavadlav, M. Praprotnik, and P. Koumoutsakos, Accelerated simulations of molecular systems through learning of effective dynamics, J. Chem. Theory Comput. 18, 538 (2022).
  * Sun  _et al._ [2022] L. Sun, J. Vandermause, S. Batzner, Y. Xie, D. Clark, W. Chen, and B. Kozinsky, Multitask machine learning of collective variables for enhanced sampling of rare events, J. Chem. Theory Comput. 18, 2341 (2022).
  * Bonati  _et al._ [2021] L. Bonati, G. Piccini, and M. Parrinello, Deep learning the slow modes for rare events sampling, Proc. Natl. Acad. Sci. U.S.A. 118, e2113533118 (2021).
  * Müllender  _et al._ [2023] L. Müllender, A. Rizzi, M. Parrinello, P. Carloni, and D. Mandelli, Effective data-driven collective variables for free energy calculations from metadynamics of paths (2023), arXiv:2311.05571 [physics.comp-ph] .
  * Vandenhaute  _et al._ [2023] S. Vandenhaute, M. Cools-Ceuppens, S. DeKeyser, T. Verstraelen, and V. Van Speybroeck, Machine learning potentials for metal-organic frameworks using an incremental learning approach, npj Comput. Mater. 9, 19 (2023).
  * Liao and Smidt [2023] Y.-L. Liao and T. Smidt, Equiformer: Equivariant graph attention transformer for 3d atomistic graphs (2023), arXiv:2206.11990 [cs.LG] .
  * Pelaez  _et al._ [2024] R. P. Pelaez, G. Simeon, R. Galvelis, A. Mirarchi, P. Eastman, S. Doerr, P. Thölke, T. E. Markland, and G. D. Fabritiis, Torchmd-net 2.0: Fast neural network potentials for molecular simulations (2024), arXiv:2402.17660 [cs.LG] .
  * Batatia  _et al._ [2023] I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, W. J. Baldwin, N. Bernstein, A. Bhowmik, S. M. Blau, V. Cărare, J. P. Darby, S. De, F. D. Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, E. Fako, A. C. Ferrari, A. Genreith-Schriever, J. George, R. E. A. Goodall, C. P. Grey, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. Holm, J. Jaafar, S. Hofmann, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O’Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. L. Schaaf, C. Schran, E. Sivonxay, T. K. Stenczel, V. Svahn, C. Sutton, C. van der Oord, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, F. Zills, and G. Csányi, A foundation model for atomistic materials chemistry (2023), arXiv:2401.00096 [physics.chem-ph] .
  * Deng  _et al._ [2023] B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel, and G. Ceder, Chgnet as a pretrained universal neural network potential for charge-informed atomistic modelling, Nature Machine Intelligence 5, 1031 (2023).
  * Chen and Ong [2022] C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nature Computational Science 2, 718 (2022).

See pages 1 of SI.pdf See pages 2 of SI.pdf See pages 3 of SI.pdf See pages 4
of SI.pdf See pages 5 of SI.pdf