11institutetext: Centre of Mathematics (CMAT), University of Minho, Braga,
4710-057, Portugal
22institutetext: Department of Mechanical Engineering (Section of Mathematics)
- FEUP, University of Porto, Porto, 4200-465, Portugal
22email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# A Self-Adaptive Penalty Method for Integrating Prior Knowledge Constraints
into Neural ODEs

C. Coelho 11 0009-0009-4502-937X M. Fernanda P. Costa 11 0000-0001-6235-286X
L.L. Ferrás 1122 0000-0001-5477-3226

###### Abstract

The continuous dynamics of natural systems has been effectively modelled using
Neural Ordinary Differential Equations (Neural ODEs). However, for accurate
and meaningful predictions, it is crucial that the models follow the
underlying rules or laws that govern these systems. In this work, we propose a
self-adaptive penalty algorithm for Neural ODEs to enable modelling of
constrained natural systems. The proposed self-adaptive penalty function can
dynamically adjust the penalty parameters. The explicit introduction of prior
knowledge helps to increase the interpretability of Neural ODE -based models.
We validate the proposed approach by modelling three natural systems with
prior knowledge constraints: population growth, chemical reaction evolution,
and damped harmonic oscillator motion. The numerical experiments and a
comparison with other penalty Neural ODE approaches and _vanilla_ Neural ODE,
demonstrate the effectiveness of the proposed self-adaptive penalty algorithm
for Neural ODEs in modelling constrained natural systems. Moreover, the self-
adaptive penalty approach provides more accurate and robust models with
reliable and meaningful predictions.

###### Keywords:

Neural ODEs Constrained Optimisation Natural Systems.

##  1 Introduction

Accurately predicting the behaviour of natural systems requires mathematical
equations that describe the relationships between their variables over time.
Due to its continuous-time dynamics, these systems are often formulated by
Ordinary Differential Equations (ODEs). ODEs provide a robust framework for
understanding and analysing the continuous changes that occur in natural
systems.

However, traditional methods of modelling natural systems using ODEs can be
challenging due to the complexity of the interactions between variables. In
recent years, Neural Networks (NNs) have been used to model the dynamics of
natural systems based on experimental data. However, the resulting models are
typically discrete and may not be suitable for modelling continuous-time
dynamics. Specifically, the discrete nature of NN models limits their ability
to capture the complex interactions between variables in natural systems,
which are often continuous [3].

Neural ODEs [3] are a NN architecture that introduces the concept of adjusting
the solution curve of a continuous-time function (ODE)
$\boldsymbol{f_{\theta}}$ to the data. The function $\boldsymbol{f_{\theta}}$
is defined by a NN with parameters $\boldsymbol{\theta}$, and the result of
training a Neural ODE is an ODE. To make predictions, an Initial Value Problem
(IVP) is solved [3]. Some natural systems modelled with Neural ODEs are
reported in the literature [11, 14]. In [11] the authors use a Neural ODE to
optimise the kinetic parameters of chemical reactions, demonstrating the
effectiveness of this approach in modelling complex combustion phenomena. In
[14] the authors propose using a Neural ODE to model glucose levels based on
sparse training data, which exhibits high precision in predicting blood
glucose trends, demonstrating its potential for practical applications in
patient care.

Natural systems often have governing laws that are mathematically expressed as
constraints. However, when we create mathematical models for these systems, we
may encounter complex and nonlinear dependencies among variables, which can
make the modelling process difficult and lead to inaccurate predictions.

Neural networks (NNs) are commonly referred to as ”black-box” models because
they are complex and opaque. It can be challenging to understand how the model
uses the training data to make predictions, which raises concerns about their
interpretability and reliability. This lack of transparency can lead to
distrust in the scientific community, as there is no guarantee that NN-based
models satisfy to the underlying governing laws of these natural systems.

Incorporating prior knowledge into NNs can reduce the need for large amounts
of training data, improve generalisation, predictive performance and reduce
the risk of overfitting.

In the literature, several strategies have been proposed to incorporate prior
knowledge constraints into NNs [13]. A popular and well-established constraint
optimisation approach is the use of penalty methods. Penalty methods have
already been used in NNs and Neural ODEs [8, 12]. For instance in [8] the
state of health of a lithium-ion battery is predicted by introducing
constraints related to the charging profile of the battery. In [12] the
authors use a penalty method to introduce prior knowledge to guarantee safety
and performance when modelling industrial systems.

Penalty methods are very straightforward to implement and easy to incorporate
into standard NN training algorithms. However, a major challenge is the choice
of the penalty parameter $\mu$. The penalty parameter plays a crucial role in
controlling the trade-off between the original loss function and the penalty
term, influencing the overall training process and the model’s final
performance. Thus, selecting an appropriate value for $\mu$ is not a trivial
task, as it requires a careful balance. To address the issue of selecting an
appropriate $\mu$, adaptive penalty functions have been proposed. These
functions are designed to dynamically adjust penalty parameters during the
optimisation process [1, 7].

This work presents a new approach to modelling constrained natural systems
using Neural ODEs. We introduce a new self-adaptive penalty function and
algorithm that adjusts the penalty parameters $\mu$ based on constraints
violation (at each iteration). The proposed method applies stronger penalties
to heavily violated constraints and weaker penalties to slightly violated
ones, balancing the objective function minimisation with constraint
satisfaction. By adapting $\mu$ over time, our approach is more robust and
adaptable to possible changes in input data distribution. Our contributions
are significant steps towards more effective and flexible modelling of
constrained natural systems using Neural ODEs.

This paper is organised as follows. Section 2 provides a brief introduction to
Neural ODEs and penalty methods. Section 3 presents the new self-adaptive
penalty function and Section 4 outlines the self-adaptive penalty algorithm
for Neural ODEs. Section 5 presents the results obtained for some numerical
experiments. The conclusions and future work are presented in Section 6.

##  2 Background

Consider a time series describing a natural system
$\boldsymbol{X}=(\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{N})$,
with $\boldsymbol{x}_{n}\in\mathbb{R}^{d}$ at time step $t_{n}$
($n=1,\dots,N$). Let
$\boldsymbol{Y}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\dots,\boldsymbol{y}_{N})$
be the corresponding ground-truth output time series, with
$\boldsymbol{y}_{n}\in\mathbb{R}^{d^{*}}$ and
$\boldsymbol{\hat{Y}}=(\boldsymbol{\hat{y}}_{1},\boldsymbol{\hat{y}}_{2},\dots,\boldsymbol{\hat{y}}_{N})$
is the prediction, with $\boldsymbol{\hat{y}}_{n}\in\mathbb{R}^{d^{*}}$, at
time step $t_{n}$.

###  2.1 Neural ODEs

Chen et al. proposed Neural ODEs, a NN architecture that adjusts a continuous-
time function (ODE) to the dynamics of the training data [3].

A Neural ODE is composed of two components, a NN that builds the ODE dynamics,
$\boldsymbol{f}_{\boldsymbol{\theta}}$, and a numerical ODE solver. During
training, the parameters $\boldsymbol{\theta}$ of the NN are adjusted by
comparing the ground-truth values with the predictions made by using the ODE
solver to solve an IVP over a time interval, $(t_{0},t_{N})$:

$\\{\boldsymbol{\hat{y}}_{n}\\}_{n=1\dots
N}=ODESolve(\boldsymbol{f_{\theta}},\boldsymbol{x}_{0},(t_{0},t_{N})),$

where $(\boldsymbol{x}_{0},t_{0})$ is the initial condition of the IVP,
$x(t_{0})=x_{0}$, $\boldsymbol{\hat{y}}_{n}$ is the solution at an arbitrary
time step $t_{n}$ with $n=1,\dots,N$. Therefore, the solver produces
predictions that are adjusted by the NN to be as close as possible to the
ground-truth data $\boldsymbol{Y}$.

###  2.2 Approaches for Handling Constrained Systems

In general, a constrained natural system can be formulated as the following
constrained optimisation problem:

—l—[0] θ ∈R^n_θl(θ) c_t^i(θ)= 0, i ∈ε, t=t_1, …, t_N c_t^j(θ) ≤0, j ∈I ,
t=t_1, …, t_N,

where $l:\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}$ is the objective
function,
$\boldsymbol{c}_{t}^{i},\boldsymbol{c}_{t}^{j}:\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}$
are the equality and inequality constraint functions, respectively, with
$\varepsilon$ the equality and $\mathcal{I}$ inequality index sets of
constraints. The set of points that satisfy all the constraints defines the
feasible set
$\mathcal{S}=\\{\boldsymbol{\theta}\in\mathbb{R}^{n_{\theta}}:\boldsymbol{c}_{t}^{i}(\boldsymbol{\theta})=0,\,i\in\varepsilon;\,\,\boldsymbol{c}_{t}^{j}(\boldsymbol{\theta})\leq
0,\,j\in\mathcal{I},t=t_{1},\dots,t_{N}\\}$.

The constrained optimisation problem (2.2), can be rewritten and solved as an
unconstrained problem using penalty methods. Penalty methods combine the
objective function $l$ and the constraints into a penalty function [9]. For
instance, using a quadratic penalty function,

—l—[0] θ∈R^n_θl(θ) + μ2 P(θ)

with
$P(\boldsymbol{\theta})=\sum_{i\in\varepsilon}\dfrac{1}{N}\sum_{n=1}^{N}|c_{t_{n}}^{i}(\boldsymbol{\theta})|^{2}+\sum_{j\in\mathcal{I}}\dfrac{1}{N}\sum_{n=1}^{N}\left([c_{t_{n}}^{j}(\boldsymbol{\theta})]^{+}\right)^{2}$,
where $\mu>0$ is the penalty parameter and $[z]^{+}$ denotes $\max(z,0)$.

Note that, to generalise penalty methods for time series data, we compute
penalty terms as the average of the quadratic constraints violation for all
time steps $t_{n}$ with $n=1,\dots,N$.

Generally, the penalty function is minimised for a sequence of increasingly
larger values of $\mu$ until a solution is found [10]. Selecting an
appropriate initial value of $\mu$ can be challenging. Large values of $\mu$
enforce constraints more strictly, but can result in slower convergence.
Conversely, smaller values of $\mu$ penalize constraint violations more
lightly, potentially leading to infeasible solutions [1]. Thus, finding the
right balance for $\mu$ is crucial to achieve both constraint satisfaction and
efficient training. However, this process often involves a tedious trial and
error approach due to the trade-off between constraint adherence and
convergence speed.

When applied to NNs, penalty methods incorporate constraints by adding a
penalty term into the loss function, resulting in a penalty loss function. The
penalty function is defined by a measure of the distance between the ground-
truth and the predicted values, $l$, and a measure of the violation of the
constraints, $\mu P$. In practice, it is common to use a fixed value of $\mu$
during all the NN training process [2, 8]. However, as mentioned previously,
selecting an appropriate value of $\mu$ is challenging, and using a fixed
value that remains unchanged during the optimisation process can lead to some
problems. One of the primary issues is that a fixed $\mu$ may not be optimal
for all stages of training. In the early stages, when the NN is far from the
optimal solution, using a large $\mu$ could lead to over-penalisation of
constraints, impeding the model’s ability to explore the solution space
effectively. On the other hand, in later stages when the NN is close to
convergence, a fixed $\mu$ might not provide enough constraint enforcement,
leading to sub-optimal solutions.

##  3 Self-Adaptive Penalty Function

The aim of this study is to develop a self-adaptive penalty function,
$\phi(\boldsymbol{\theta})$, which dynamically adjusts the penalty parameters
$\mu$ taking into account the constraints violation during the training
process of the NN. The proposed self-adaptive penalty function adapts the
penalty parameters to the degree of constraint violation observed at the
current iteration.

Problem (2.2) is rewritten as an unconstrained problem,

—l—[0] θ∈R^n_θϕ(θ) ,

where the self-adaptive penalty function $\phi$ relates the objective function
$l$ and the constraints $\boldsymbol{c}_{t}^{i},\boldsymbol{c}_{t}^{j}$.

To define $\phi$, the objective function $l$ and the constraint violation
values at each point $\theta$, are normalised using a bounded function $\psi$,
to ensure that all values have the same order of magnitude. In this work, we
propose the following normalisation function $\psi(x)=1-\dfrac{1}{1+x}$ where
$0\leq\psi(x)\leq 1$, for all $x\in\mathbb{R}_{0}^{+}$.

The objective function value $l$ at each point $\boldsymbol{\theta}$ is
normalised using $\psi$ resulting in the new objective function
$F_{\boldsymbol{\theta}}$:

$F_{\boldsymbol{\theta}}=\psi(l(\boldsymbol{\theta})).$ (1)

For each equality constraint $c^{i}_{t}$, the constraint violation vector
$\boldsymbol{v}^{i}\in\mathbb{R}^{\varepsilon}$, at each point
$\boldsymbol{\theta}$, is given by,
$\boldsymbol{v}^{i}=|c_{t}^{i}(\boldsymbol{\hat{y}}_{t})|^{2}_{t=t_{1}\dots
t_{N}}.$ Then, the normalisation step $\psi$ is applied and the total
violation $P_{i}$ is given by the average of the violations $v_{t}^{i}$, with
$t=t_{1},\dots,t_{N}$:

$P_{i}=\dfrac{1}{N}\sum_{n=1}^{N}\psi(v^{i}_{t_{n}}).$ (2)

Likewise, for each inequality constraint $c^{j}_{t}$ the constraint violation
vector $\boldsymbol{v}^{j}\in\mathbb{R}^{\mathcal{I}}$, at each point
$\boldsymbol{\theta}$, is defined by,
$\boldsymbol{v}^{j}=([c_{t}^{j}(\boldsymbol{\hat{y}}_{t})]^{+})^{2}_{t=t_{1},\dots,t_{N}}.$
Then, the normalisation step $\psi$ is applied and the total violation $P_{j}$
is given by the average of the violations $v^{j}_{t}$, with
$t=t_{1},\dots,t_{N}$:

$P_{j}=\dfrac{1}{N}\sum_{n=1}^{N}\psi(v^{j}_{t_{n}}).$ (3)

For each equality and inequality constraints, penalty parameters $\mu_{i}$ and
$\mu_{j}$, respectively, are defined by the proportion of predictions that
violate the constraint at all time steps $t=t_{1},\dots,t_{N}$, at the current
iteration:

$\mu_{i}=\dfrac{\\#\\{t:v_{t}^{i}\neq
0\\}}{N},\,\,\,\mu_{j}=\dfrac{\\#\\{t:v_{t}^{j}\neq 0\\}}{N},$ (4)

where $\\#\\{z\\}$ denotes the cardinality of set $z$.

One of the advantages of computing self-adaptive penalty parameters is that it
does not require to be provided by the user. Instead, they are computed using
the information, gathered at each iteration, on how many predictions violate a
constraint. A constraint that is violated in a higher number of time steps
(than any other constraint) will have a larger penalty parameter.
Additionally, the normalisation step $\psi$ prevents numerical instability and
improves the accuracy of the training process.

Finally, the self-adaptive penalty function $\phi$ is dynamically defined at
each iteration as follows. The $\phi$ function is given by the
$F_{\boldsymbol{\theta}}$ function, if the current point $\boldsymbol{\theta}$
is a feasible point. Otherwise, $\phi$ is defined by $F_{\boldsymbol{\theta}}$
plus the penalty terms $P_{i}$ and $P_{j}$ multiplied by the self-adaptive
penalty parameters $\mu_{i}$ and $\mu_{j}$. Specifically, a penalty term per
constraint violation, given by the normalised sum of the constraint violation
vector $\boldsymbol{v}^{i},\boldsymbol{v}^{j}$ multiplied by its corresponding
penalty parameter $\mu_{i},\mu_{j}$ for equality and inequality constraints,
respectively.

$\phi(\boldsymbol{\theta})=\begin{cases}F_{\boldsymbol{\theta}},&\text{if
}\boldsymbol{\theta}\in\mathcal{S},\cr
F_{\boldsymbol{\theta}}+\dfrac{1}{\\#\\{\varepsilon\\}}\sum_{i\in\varepsilon}\mu_{i}P_{i}+\dfrac{1}{\\#\\{\mathcal{I}\\}}\sum_{j\in\mathcal{I}}\mu_{j}P_{j},&\text{if
}\boldsymbol{\theta}\notin\mathcal{S},\end{cases}$ (5)

with $F_{\theta}$, $P_{i}$, $P_{j}$, $\mu_{i}$ and $\mu_{i}$ given by (1)-(4),
respectively.

We now demonstrate the equivalence of problems (2.2) and (3), thereby
establishing that they possess identical global minimisers.

###### Theorem 3.1

Let $\boldsymbol{\theta^{*}}$ be a global solution of (2.2) and let
$\boldsymbol{\omega}\in\mathcal{S}$ be such that $l(\boldsymbol{\omega})\geq
l(\boldsymbol{\theta}^{*})$. Then, $\boldsymbol{\theta^{*}}$ is a global
solution to (3), where $\phi$ is the penalty function defined in (5).

###### Proof

Let $\boldsymbol{\theta^{*}}\in\mathcal{S}$ be a global solution of (2.2). By
definition, we have $l(\boldsymbol{\theta^{*}})\leq l(\boldsymbol{\theta})$
for all $\boldsymbol{\theta}\in\mathcal{S}$. Therefore, for all
$\boldsymbol{\theta}\in\mathcal{S}$, we have:

$\phi(\boldsymbol{\theta^{*}})=1-\dfrac{1}{1+l(\boldsymbol{\theta^{*}})}\leq
1-\dfrac{1}{1+l(\boldsymbol{\theta})}=\phi(\boldsymbol{\theta}).$

Now, consider the case when
$\boldsymbol{\theta}\in\mathbb{R}^{n_{\theta}}{\setminus}\mathcal{S}$.
Assuming that $l(\boldsymbol{\theta})\leq l(\boldsymbol{\omega})$, we have,

$\leavevmode\resizebox{433.62pt}{}{$\phi(\boldsymbol{\theta}^{*})=\psi(l(\boldsymbol{\theta^{*}}))\leq\psi(l(\boldsymbol{\omega}))\leq\psi(l(\boldsymbol{\theta}))\leq\psi(l(\boldsymbol{\theta}))+\dfrac{1}{\\#\\{\varepsilon\\}}\sum_{i\in\varepsilon}\mu_{i}P_{i}+\dfrac{1}{\\#\\{\mathcal{I}\\}}\sum_{j\in\mathcal{I}}\mu_{j}P_{j}=\phi(\boldsymbol{\theta}),$
}$ (6)

which implies $\phi(\boldsymbol{\theta^{*}})\leq\phi(\boldsymbol{\theta})$ for
all $\boldsymbol{\theta}\in\mathbb{R}^{n_{\theta}}$, i.e.,
$\boldsymbol{\theta^{*}}$ is a global solution of (3).

###### Lemma 1

If $\boldsymbol{\theta^{*}}$ is a global solution of (3), where $\phi$ is the
penalty function defined in (5), then $\boldsymbol{\theta^{*}}$ is a feasible
point of (2.2).

###### Proof

By contradiction, we assume that
$\boldsymbol{\theta^{*}}\in\mathbb{R}^{n_{\theta}}{\setminus}\mathcal{S}$.
When $l(\boldsymbol{\theta^{*}})\leq l(\boldsymbol{\omega})$ and
$\boldsymbol{\omega}\in\mathcal{S}$ we get, from (5),

$\leavevmode\resizebox{433.62pt}{}{$\phi(\boldsymbol{\theta^{*}})=\psi(l(\boldsymbol{\theta^{*}}))+\dfrac{1}{\\#\\{\varepsilon\\}}\sum_{i\in\varepsilon}\mu_{i}P_{i}+\dfrac{1}{\\#\\{\mathcal{I}\\}}\sum_{j\in\mathcal{I}}\mu_{j}P_{j}>\psi(l(\boldsymbol{\theta^{*}}))>\psi(l(\boldsymbol{\omega}))=\phi(\boldsymbol{\omega}),$
}$ (7)

which contradict the definition of a global solution of (3). Therefore,
$\boldsymbol{\theta^{*}}\in\mathcal{S}$.

We are now able to establish the reciprocal of Theorem 3.1.

###### Theorem 3.2

Let $\boldsymbol{\theta^{*}}\in\mathcal{R}^{n_{\theta}}$ be a global solution
of (3), where $\phi$ is the penalty function defined by (5). Then,
$\boldsymbol{\theta^{*}}$ is a global solution of (2.2).

###### Proof

By Lemma 1
$\boldsymbol{\theta^{*}}\in\mathcal{S}\subset\mathbb{R}^{n_{\theta}}$. We have
$\psi(l(\boldsymbol{\theta^{*}}))=\phi(\boldsymbol{\theta^{*}})\leq\phi(\boldsymbol{\theta})$
for all $\boldsymbol{\theta}\in\mathbb{R}^{n_{\theta}}$, and in particular,
for all $\boldsymbol{\theta}\in\mathcal{S}$, we have
$\psi(l(\boldsymbol{\theta^{*}}))\leq\psi(l(\boldsymbol{\theta}))$, which
implies $l(\boldsymbol{\theta^{*}})\leq l(\boldsymbol{\theta})$. Therefore,
$\boldsymbol{\theta^{*}}$ is a global solution of (2.2).

##  4 Self-adaptive Penalty Algorithm for Neural ODEs

In this section, we introduce an algorithm that implements the self-adaptive
penalty method for Neural ODEs to model natural systems with prior knowledge
constraints. The algorithm is straightforward and requires minimal
modification to the traditional Neural ODE training algorithm presented in
[3].

Throughout the process of optimising the neural network, we keep track of the
point $\boldsymbol{\theta}$ that produces the smallest value of $\phi$
(denoted as $\phi_{\text{best}}$) and store it as
$\boldsymbol{\theta}_{\text{best}}$. This serves as a guide during the
optimization process to prevent the solution from deteriorating. If a new
point does not improve upon $\boldsymbol{\theta}_{\text{best}}$, we discard it
and continue the optimisation process using
$\boldsymbol{\theta}_{\text{best}}$.

The stopping criteria is the maximum number of iterations, $k_{\max}$. When
$k_{\max}$ is reached, the final parameters of the NN that build the ODE
dynamics are $\boldsymbol{\theta}_{\text{best}}$.

The proposed self-adaptive penalty algorithm for Neural ODEs is presented in
Algorithm 1.

Algorithm 1 : The self-adaptive penalty algorithm for Neural ODEs.

Input: Initial condition $(\boldsymbol{y}_{0},t_{0})$, start time $t_{0}$, end
time $t_{N}$, maximum number of iterations $k_{max}$;

$tol\leftarrow 1e-4$;

$\psi(x)=1-\dfrac{1}{1+x}$;

$\boldsymbol{f_{\theta}}=DynamicsNN()$;

Initialise parameters $\boldsymbol{\theta}$;

$\phi_{\text{best}}\leftarrow+\infty$;

for $k=1:k_{max}$ do

{$\boldsymbol{\hat{y}}_{n}\\}_{n=1,\dots,N}=ODESolve(\boldsymbol{f_{\theta}},\boldsymbol{y}_{0},(t_{0},t_{N}))$;

Evaluate normalised objective function $F_{\boldsymbol{\theta}}$;

$P_{\theta}\leftarrow 0$;

for $i\in\varepsilon$ do

Compute penalty term $P_{i}$ using (2);

Compute self-adaptive penalty parameter $\mu_{i}$ using (4);

$P_{\boldsymbol{\theta}}\leftarrow P_{\boldsymbol{\theta}}+P_{i}$;

end for

for $j\in\mathcal{I}$ do

Compute penalty term $P_{j}$ using (3);

Compute self-adaptive penalty parameter $\mu_{j}$ using (4);

$P_{\boldsymbol{\theta}}\leftarrow P_{\boldsymbol{\theta}}+P_{j}$;

end for

if $P_{\boldsymbol{\theta}}\leq tol$ then

$\phi\leftarrow F_{\boldsymbol{\theta}}$;

else

$\phi\leftarrow
F_{\boldsymbol{\theta}}+\dfrac{1}{\\#\\{\varepsilon\\}}\sum_{i\in\varepsilon}\mu_{i}P_{i}+\dfrac{1}{\\#\\{\mathcal{I}\\}}\sum_{j\in\mathcal{I}}\mu_{j}P_{j}$;

end if

if $\phi_{\theta}<\phi_{\text{best}}$ then

$\boldsymbol{\theta}_{\text{best}}\leftarrow\boldsymbol{\theta}$;

$\phi_{\text{best}}\leftarrow\phi_{\boldsymbol{\theta}}$;

end if

$\nabla\phi_{\text{best}}\leftarrow
Optimiser.BackpropCall(\phi_{\text{best}})$;

$\boldsymbol{\theta}\leftarrow
Optimiser.Step(\nabla\phi_{\text{best}},\boldsymbol{\theta}_{\text{best}})$;

end for

return $\boldsymbol{\theta}_{\text{best}}$;

##  5 Numerical Experiments

To assess the effectiveness of using the self-adaptive penalty algorithm for
Neural ODEs, three synthetic datasets describing natural systems with prior
knowledge constraints were used. Namely, the World Population Growth (WPG)
[6], Chemical Reaction (CR) [5], and Damped Harmonic Oscillator (DHO) [4].
Note that the selected datasets were chosen because they contained valuable
information regarding the constraints of the system, which is typically not
readily available. Three distinct experiments were carried out to test and
validate the models trained using the proposed approach: (a) reconstructing
the data, where the same dataset was used for both training and testing, (b)
extrapolating to a longer time horizon, where the model was tested on a larger
time interval, (c) completion, where the test set had a different sampling
frequency in the same time interval.

As a form of comparison, we consider a _vanilla_ Neural ODE (a traditional
Neural ODE without incorporating constraints) and a Neural ODE with a
quadratic penalty function, using three different $\mu$ values
($\mu=1,10,100$) to demonstrate the influence of $\mu$ on the model’s
performance, which serve as our baselines. The assessment of the models was
based on the Mean Squared Error (MSE) and mean total constraints violation
$P$.

Tables 1-3 show the average values of MSE and $P$, and their respective
standard deviation (std), after three independent runs, for self-adaptive
penalty algorithm for Neural ODEs and two baselines, _vanilla_ Neural ODE and
a Neural ODE with quadratic penalty function. The best MSE and $P$ values at
each task are written in bold. Figures 1-3 show the predicted and real curves
for the best run at each task.

In the experiments, the objective function $l$ is defined by,
$l(\boldsymbol{\theta})=\sum_{n=1}^{N}\left(\boldsymbol{\hat{y}}_{n}-\boldsymbol{y}_{n}\right)^{2}.$

###  5.1 World Population Growth

The WPG dataset is publicly available on _Kaggle_ [6]. It is a time series
composed of two features (time step $t$, population $p(t)$) which describes a
system with an inequality constraint defined by the carrying capacity at each
time step $t$. Ideally, all data points of the adjusted model must satisfy
this constraint, $\\{p(t)\leq 12\\}_{t=t_{1},\dots,t_{N}}.$

The trained NNs have $4$ hidden layers: linear with $50$ neurons; hyperbolic
tangent (tanh); linear with $50$ neurons; Exponential Linear Unit (ELU). The
input and output layers have $1$ neuron. The Adam optimiser was used with a
learning rate of $1e-5$. Training was done for $10000$ iterations
($k_{\max}$). The reconstruction task uses the same set of $200$ data points
in the time interval $[0,300]$ for both training and testing. In the
extrapolation task, the training set comprises $200$ data points in the time
interval $[0,300]$, while the testing set consists of $200$ data points in the
interval $[0,400]$. The completion task involves $200$ training points and
$300$ testing points in the interval $[0,300]$.

From Table 1, as expected, the vanilla Neural ODEs show higher $P$ and MSE
values than the quadratic penalty ($\mu=1$) and self-adaptive versions. Hence,
incorporating constraints into Neural ODEs not only ensures they are
satisfied, but also contribute to better data fitting.

In general, in Neural ODEs with a quadratic penalty function, a larger $\mu$
results in models with smaller constraints violation value $P$ but with a
higher MSE. Inversely, a smaller value of $\mu$ produces models with lower MSE
values but a higher $P$ value.

Neural ODEs with self-adaptive penalty algorithm exhibit the best performance
in all experiments, with the smallest MSE and $P$ values. Thus, effectively
balancing the minimisation of the objective function (fit to the training
data) and the constraints violation.

Figure 1 displays the predicted and real curves of the testing set for each
task’s best run. The top row shows the results for the _vanilla_ Neural ODE,
which was unable to capture the inequality constraint from the training data.
This is evident from the exponential increase in population values over time.
In contrast, the middle row shows the results for Neural ODEs with a quadratic
penalty function, which do not exhibit this behaviour. The bottom row,
displays the results for the Neural ODE with a self-adaptive penalty
algorithm, which closely fits the real curve and demonstrates its superiority.

Table 1: Performance on WPG dataset of _vanilla_ Neural ODE, Neural ODE with a
quadratic penalty function ($\mu=1,10,100$), and self-adaptive penalty
algorithm for Neural ODE.

_Vanilla_ Neural ODE Quadratic Penalty Function Self-adaptive Penalty
Algorithm $\mu=1$ $\mu=10$ $\mu=100$ Experiment $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$
$\text{MSE}\pm\text{std}$ $P\pm\text{std}$ $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$ Reconstruction
1.9e-1 $\pm$ 1.5e-1 7.8e-2 $\pm$ 5.1e-2 4.7e-2 $\pm$ 5e-2 9.1e-3 $\pm$ 1.5e-2
4.6e-2 $\pm$ 4e-2 1.6e-4 $\pm$ 5e-05 4.1e-2 $\pm$ 1.9e-2 2e-7 $\pm$ 1e-7
1.3e-3 $\pm$ 4.3e-4 6.8e-4 $\pm$ 2.3e-4 Extrapolation 1.3e-1 $\pm$ 7.2e-2
1.5e-1 $\pm$ 5.3e-2 2.6e-2 $\pm$ 1.9e-2 3.2e-2 $\pm$ 1.5e-2 1.4e-1 $\pm$
1.3e-1 4.5e-2 $\pm$ 4.e-2 8.8e-2 $\pm$ 4.3e-2 2.5e-2 $\pm$ 7.3e-3 1.4e-3 $\pm$
1.8e-3 1e-2 $\pm$ 3.9e-4 Completion 4.3e-2 $\pm$ 2.9e-3 2.8e-2 $\pm$ 1.2e-2
2.6e-2 $\pm$ 1.8e-2 2.8e-3 $\pm$ 2.5e-3 3.1e-2 $\pm$ 1e-2 6e-5 $\pm$ 7e-5
2.5e-1 $\pm$ 1e-1 2e-5 $\pm$ 3e-5 8.9e-4 $\pm$ 5e-5 4.4e-4 $\pm$ 5e-5

Figure 1: Plot of the real (dashed line) and predicted (solid line) curves,
for the WPG, on different tasks trained with with _vanilla_ Neural ODE (top),
Neural ODE with a quadratic penalty function $\mu=1$ (middle), self-adaptive
penalty algorithm for Neural ODE (bottom): (a)(d)(g) Reconstruction; (b)(e)(h)
Extrapolation; (c)(f)(i) Completion.

###  5.2 Chemical Reaction

The CR dataset is publicly available on _Kaggle_ , which describes the
evolution of the species in a chemical reaction [5]. It is a time-series
composed of five features (time step $t$, and the masses of species A,
$m_{A}$, B, $m_{B}$, C, $m_{C}$, and D, $m_{D}$). At each time step $t$
($t=t_{1},\dots,t_{N}$), this system has an equality constraint defined by the
conservation of mass
$\\{m_{A}(t)+m_{B}(t)+m_{C}(t)+m_{D}(t)=m_{\text{total}}\\}_{t=t_{1},\dots,t_{N}},$
with $m_{\text{total}}$ the total constant mass of the system. Ideally, all
data points of the adjusted model must satisfy this constraint,

The NNs were trained for $10000$ iterations ($k_{\max}$) and have $6$ hidden
layers: linear with $50$ neurons; tanh; linear with $64$ neurons; ELU; linear
with $50$ neurons; tanh. The input and output layers have $4$ neurons. The
Adam optimiser was used with a learning rate of $1e-5$. For the reconstruction
task, both the training and testing sets consist of $100$ data points in the
time interval $[0,100]$. The extrapolation task involves training with $100$
data points in the time interval $[0,100]$ and testing with $100$ data points
in the interval $[0,200]$. The completion task uses $100$ training points and
$200$ testing points in the interval $[0,100]$.

Table 2 indicates that the performance of the _vanilla_ Neural ODE is
comparable to that of Neural ODE with a quadratic penalty function and
$\mu=1$. Increasing $\mu$ to $10$ and $100$ slightly raises MSE values but
also improves constraint violation, except for the extrapolation task with
$\mu=100$.

The Neural ODE with self-adaptive penalty algorithm was found to be the best
model for the CR dataset, with very low values of both mean squared error
(MSE) and constraint violation across all three tasks. These results suggest
that the equality constraint is effectively incorporated into the model while
also providing the best fit to the data. This is supported by the comparison
of plots in Figure 2, which demonstrate that the Neural ODE with self-adaptive
penalty algorithm 2-(h) performs well even for the more challenging
extrapolation task, accurately modelling the dynamics of the data. In
contrast, the other two models 2-(b) and 2-(e) deviate from the true dynamics
for a time horizon beyond that seen during training ($100s$).

Figure 2 shows the difficulty of _vanilla_ Neural ODE (top row) to fit the
real curves, achieving a better result for the completion task. The plots of
the predicted curves of Neural ODE with a self-adaptive penalty algorithm
(bottom row) show the best modelling performance for the three tasks.

Table 2: Performance on CR dataset of _vanilla_ Neural ODE, Neural ODE with a
quadratic penalty function ($\mu=1,10,100$), and self-adaptive penalty
algorithm for Neural ODE.

Vanilla Neural ODE Quadratic Penalty Function Self-adaptive Penalty Algorithm
$\mu=1$ $\mu=10$ $\mu=100$ Experiment $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$
$\text{MSE}\pm\text{std}$ $P\pm\text{std}$ $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$ Reconstruction
2.1e-3 $\pm$ 1.4e-3 7.8e-3 $\pm$ 6.7e-3 1.2e-2 $\pm$ 2.1e-2 5.3e-3 $\pm$
1.2e-2 2.9e-2 $\pm$ 3.5e-2 3.1e-4 $\pm$ 4e-4 5.7e-2 $\pm$ 3.8e-2 8e-5 $\pm$
1.1e-4 2e-4 $\pm$ 2e-4 2.2e-6 $\pm$ 1.4e-6 Extrapolation 4.1e-2 $\pm$ 5.7e-2
1.6e-2 $\pm$ 1.6e-2 2.3e-1 $\pm$ 2.7e-1 2e-1 $\pm$ 2.6e-1 2.1e-3 $\pm$ 1.8e-3
1.8e-3 $\pm$ 1.6e-3 7.9e-1 $\pm$ 1.1 6.7 $\pm$ 9.3 8e-5 $\pm$ 2e-5 2e-5 $\pm$
2e-5 Completion 9e-4 $\pm$ 5.9e-4 2.6e-3 $\pm$ 3.1e-3 1.7e-2 $\pm$ 2.3e-2
1.e-3 $\pm$ 1.6e-3 2.5e-2 $\pm$ 3.5e-2 3e-4 $\pm$ 4.2e-4 1.2e-1 $\pm$ 9.7e-2
2.9e-4 $\pm$ 3.8e-4 1e-3 $\pm$ 1.4e-3 1.4e-4 $\pm$ 2e-4

Figure 2: Plot of the real (dashed line) and predicted (solid line) curves,
for the CR, on different tasks trained with _vanilla_ Neural ODE (top), Neural
ODE with a quadratic penalty function $\mu=1$ (middle) and self-adaptive
penalty algorithm for Neural ODE (bottom): (a)(d)(g) Reconstruction; (b)(e)(h)
Extrapolation; (c)(f)(i) Completion.

###  5.3 Damped Harmonic Oscillator

The DHO dataset, available on _Kaggle_ , describes the motion of a mass
attached to a spring under the influence of a dissipative force described by a
second-order ODE [4]. The dataset consists of a time series with four
features: time step $t$, displacement $x$, velocity $v$, and acceleration $a$.
At each time step $t$, the system has one inequality constraint and one
equality constraint. The inequality constraint requires the system’s energy to
decrease due to damping $\\{E(t)>E(t+1)\\}_{t=t_{1},\dots,t_{N}},$ where
$E(t)=\dfrac{1}{2}m\ v(t)^{2}+\dfrac{1}{2}kx(t)^{2}$, and $k$ is the spring
constant. The equality constraint requires the conservation of the rate at
which energy is dissipated due to damping $\\{\Delta
P(t)=0\\}_{t=t_{1},t_{N}},$ where $P(t)=-cv(t)x(t)$, and $c$ is the damping
coefficient. To obtain meaningful results, all data points of the adjusted
model should satisfy these constraints.

The DHO dataset is unique among the datasets used in this study because it
involves a system with two constraints. As shown in Figure 3, both the
_vanilla_ Neural ODE and the Neural ODE with a quadratic penalty function
struggled to model this system, with predicted curves significantly deviating
from the actual curves in all three runs. In contrast, the Neural ODE with the
self-adaptive penalty algorithm accurately adjusted to the data dynamics,
particularly in the extrapolation task (Figures 3-(b)(e)(h)). The self-
adaptive penalty algorithm effectively incorporates the constraints into the
model, as demonstrated by the model’s ability to capture the physical dynamics
of the system even when predicting for longer time horizons not seen during
training. This behaviour was not observed with either the _vanilla_ Neural ODE
or the Neural ODE with a quadratic penalty function.

The trained NNs have $4$ hidden layers: linear with $50$ neurons; hyperbolic
tangent (tanh); linear with $50$ neurons; Exponential Linear Unit (ELU). The
input and output layers have $2$ neurons. The Adam optimiser was used with a
learning rate of $1e-5$ and the optimisation process was done until $10000$
iterations ($k_{\max}$) were done. In the reconstruction task, the training
and testing sets contain the same $400$ data points in the time interval
$[0,50]$. The extrapolation task involves training with $400$ data points in
the time interval $[0,50]$ and testing with $400$ data points in the interval
$[0,400]$. For the completion task, there are $400$ training points and $600$
testing points in the interval $[0,50]$.

The numerical results presented in Table 3 support the findings shown in
Figure 3. The use of a Neural ODE with a quadratic penalty function performed
the worst, failing completely at the extrapolation task. However, when
compared to a _vanilla_ Neural ODE, constraint violation values were lower for
the reconstruction and completion tasks.

In contrast, the Neural ODE with the self-adaptive penalty algorithm performed
significantly better, exhibiting a notable difference in magnitude for both
MSE and constraint violation $P$ when compared to the baselines. Moreover, the
standard deviation values indicate that this method produced the most stable
models across all three runs.

Table 3: Performance on DHO dataset of _vanilla_ Neural ODE, Neural ODE with a
quadratic penalty function ($\mu=1,10,100$), and self-adaptive penalty
algorithm for Neural ODE.

Vanilla Neural ODE Quadratic Penalty Function Self-adaptive Penalty Algorithm
$\mu=1$ $\mu=10$ $\mu=100$ Experiment $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$
$\text{MSE}\pm\text{std}$ $P\pm\text{std}$ $\text{MSE}\pm\text{std}$
$P\pm\text{std}$ $\text{MSE}\pm\text{std}$ $P\pm\text{std}$ Reconstruction
1.2e-1 $\pm$ 2.9e-2 1.2e-3 $\pm$ 1.5e-3 2.7e-1 $\pm$ 2.3e-1 5e-6 $\pm$ 7.1e-6
2.2e-1 $\pm$ 8.3e-2 1.3e-4 $\pm$ 1.8e-4 2.3e-1 $\pm$ 1.1e-1 1e-5 $\pm$ 1e-5
8e-5 $\pm$ 1e-4 6.4e-3 $\pm$ 3e-5 Extrapolation 8.5 $\pm$ 11.9 13.9 $\pm$ 19.7
220.5 $\pm$ 250.2 78.7 $\pm$ 85.1 331.4 $\pm$ 442.9 79.5 $\pm$ 111.2 24.6
$\pm$ 34.6 3.1 $\pm$ 4.4 1e-1 $\pm$ 2e-5 8e-4 $\pm$ 1e-8 Completion 1e-1 $\pm$
9.3e-3 1.4e-4 $\pm$ 6e-5 1.2e-1 $\pm$ 1.6e-2 1e-8 $\pm$ 1e-8 15e-1 $\pm$
7.2e-3 2e-8 $\pm$ 4e-8 4.6e-1 $\pm$ 1.9e-1 1e-7 $\pm$ 3e-8 5e-5 $\pm$ 3e-5
6.3e-3 $\pm$ 1e-5

Figure 3: Plot of the real (dashed line) and predicted (solid line) curves,
for the DHO, on different tasks trained with _vanilla_ Neural ODE (top),
Neural ODE with a quadratic penalty function $\mu=1$ (middle), self-adaptive
penalty algorithm for Neural ODE (bottom): (a)(d)(g) Reconstruction; (b)(e)(h)
Extrapolation; (c)(f)(i) Completion.

##  6 Conclusion

This work proposes a self-adaptive penalty algorithm for Neural ODEs to enable
modelling of constrained natural systems. The proposed self-adaptive penalty
function can dynamically adjust the penalty parameters

Our approach is an improvement over traditional penalty methods which require
an appropriate initialisation and tuning of penalty parameters $\mu$ during
the optimization process. This selection is challenging and time-consuming. In
general, in the context of NNs, a fixed penalty parameter is used, and
consequently an optimal solution is not guaranteed to be found.

The proposed self-adaptive penalty function dynamically adjusts penalty
parameters taking into account the degree of constraints violation, resulting
in more efficient and accurate predictions, especially for complex systems
with constraints. The self-adaptive penalty function employs a normalisation
step to ensure all values (objective and constraints violations) have the same
order of magnitude, improving training stability, critical for complex systems
prone to numerical instability.

To evaluate the self-adaptive penalty algorithm for Neural ODEs, we used three
constrained natural systems and tested the performance on three tasks:
reconstruction, extrapolation, and completion. We compared the results
obtained by our algorithm with two baselines, a _vanilla_ Neural ODE and three
Neural ODEs with quadratic penalty functions with the penalty parameter values
$\mu=1,10,100$, respectively. We measured performance in terms of mean squared
error (MSE) and total constraints violation ($P$). Our algorithm achieved
remarkably higher performance, exhibiting lower MSE and reduced total
constraint violation values, indicating that the models not only fit data but
also satisfy the governing laws of constrained systems. Additionally, our
algorithm showed the ability to extrapolate and generalise well.

The proposed self-adaptive penalty function and self-adaptive penalty
algorithm for Neural ODEs can be applied to any NN architecture and represents
a promising approach for incorporating prior knowledge constraints into NNs
while enhancing interpretability. Future work will focus on modelling more
complex scenarios with more complex constraints using Neural ODEs with our
self-adaptive penalty method, as well as conducting theoretical studies on the
convergence of our approach.

####  6.0.1 Acknowledgements

The authors acknowledge the funding by Fundação para a Ciência e Tecnologia
(FCT) (Portuguese Foundation for Science and Technology) through CMAT projects
UIDB/00013/2020 and UIDP/00013/2020. C. Coelho would like to thank FCT for the
funding through the scholarship with reference 2021.05201.BD.

## References

  * [1] Ali, M.M., Zhu, W.X.: A penalty function-based differential evolution algorithm for constrained global optimization. Computational Optimization and Applications 54(3), 707–739 (Apr 2013). https://doi.org/10.1007/s10589-012-9498-3
  * [2] Cai, S., Mao, Z., Wang, Z., Yin, M., Karniadakis, G.E.: Physics-informed neural networks (PINNs) for fluid mechanics: A review. Acta Mechanica Sinica 37(12), 1727–1738 (Dec 2021). https://doi.org/10.1007/s10409-021-01148-1
  * [3] Chen, R.T., Rubanova, Y., Bettencourt, J., Duvenaud, D.K.: Neural ordinary differential equations. Advances in neural information processing systems 31 (2018)
  * [4] Coelho, C., Costa, M.F.P., Ferrás, L.: Damped harmonic oscillator dataset (2023). https://doi.org/10.34740/KAGGLE/DS/3147798, https://www.kaggle.com/ds/3147798
  * [5] Coelho, C., Costa, M.F.P., Ferrás, L.: Synthetic chemical reaction (2023). https://doi.org/10.34740/KAGGLE/DS/3010478, https://www.kaggle.com/ds/3010478
  * [6] Coelho, C., Costa, M.F.P., Ferrás, L.: World population growth (2023). https://doi.org/10.34740/KAGGLE/DS/3010437, https://www.kaggle.com/ds/3010437
  * [7] Costa, M.F.P., Francisco, R.B., Rocha, A.M.A.C., Fernandes, E.M.G.P.: Theoretical and Practical Convergence of a Self-Adaptive Penalty Algorithm for Constrained Global Optimization. Journal of Optimization Theory and Applications 174(3), 875–893 (Sep 2017). https://doi.org/10.1007/s10957-016-1042-7
  * [8] Dai, H., Zhao, G., Lin, M., Wu, J., Zheng, G.: A Novel Estimation Method for the State of Health of Lithium-Ion Battery Using Prior Knowledge-Based Neural Network and Markov Chain. IEEE Transactions on Industrial Electronics 66(10), 7706–7716 (Oct 2019). https://doi.org/10.1109/TIE.2018.2880703
  * [9] Fletcher, R., Leyffer, S., Toint, P.L., et al.: A brief history of filter methods. Preprint ANL/MCS-P1372-0906, Argonne National Laboratory, Mathematics and Computer Science Division 36 (2006)
  * [10] Nocedal, J., Wright, S.J.: Numerical Optimization. Springer Series in Operation Research and Financial Engineering, Springer, New York, NY, second edition edn. (2006)
  * [11] Su, X., Ji, W., An, J., Ren, Z., Deng, S., Law, C.K.: Kinetics Parameter Optimization via Neural Ordinary Differential Equations (arXiv:2209.01862) (Sep 2022)
  * [12] Tuor, A., Drgona, J., Vrabie, D.: Constrained Neural Ordinary Differential Equations with Stability Guarantees (arXiv:2004.10883) (Apr 2020)
  * [13] von Rueden, L., Mayer, S., Beckh, K., Georgiev, B., Giesselbach, S., Heese, R., Kirsch, B., Pfrommer, J., Pick, A., Ramamurthy, R., Walczak, M., Garcke, J., Bauckhage, C., Schuecker, J.: Informed Machine Learning – A Taxonomy and Survey of Integrating Prior Knowledge into Learning Systems. IEEE Transactions on Knowledge and Data Engineering 35(1), 614–633 (Jan 2023). https://doi.org/10.1109/TKDE.2021.3079836
  * [14] Xing, Y., Ye, H., Zhang, X., Cao, W., Zheng, S., Bian, J., Guo, Y.: A continuous glucose monitoring measurements forecasting approach via sporadic blood glucose monitoring. In: 2022 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). pp. 860–863 (Dec 2022). https://doi.org/10.1109/BIBM55620.2022.9995522


EPJ Web of Conferences Lattice2017  11institutetext: HISKP (Theory) and Bethe
Center for Theoretical Physics, University of Bonn, Germany 22institutetext:
Institut fÃr Kernphysik, Johannes Gutenberg-UniversitÃt Mainz, Germany
33institutetext: ITP, Albert Einstein Center for Fundamental Physics,
University of Bern, Switzerland

# The $\eta^{\prime}$ meson at the physical point with $N_{f}=2$ Wilson
twisted mass fermions

Christopher Helmes 11 Bastian Knippschild 11 Bartosz Kostrzewa 11 Liuming Liu
11 Christian Jost 11 Konstantin Ottnad 1122 Carsten Urbach Speaker, 11
<EMAIL_ADDRESS>Urs Wenger 33 Markus Werner
for the ETM collaboration 11

###### Abstract

We present results for the $\eta^{\prime}$ meson and the topological
susceptibility in $N_{f}=2$ flavour lattice QCD. The results are obtained
using Wilson twisted mass fermions at maximal twist with pion masses ranging
from $340\ \mathrm{MeV}$ down to the physical point. A comparison to
literature values is performed giving a handle on discretisation effects.

##  1 Introduction

Due to the persisting $3-5\ \sigma$ deviation in the anomalous magnetic moment
of the muon $a_{\mu}$ between theory and experiment there is considerable
interest in the decays $\eta\to\gamma^{\star}\gamma^{\star}$ and
$\eta^{\prime}\to\gamma^{\star}\gamma^{\star}$. A better knowledge of the
corresponding transition form factors could help to reduce the uncertainty in
the hadronic light-by-light contribution to $a_{\mu}$, see for instance Ref.
Jegerlehner:2015stw . Moreover, $\eta$ and $\eta^{\prime}$ mesons are
interesting from a theoretical point of view because the large mass of the
$\eta^{\prime}$ meson is explained by the anomalously broken $U_{A}(1)$ axial
symmetry in QCD. The $\eta,\eta^{\prime}$ mixing pattern and the
aforementioned transition form factors can be computed non-perturbatively
using lattice techniques. There has been considerable progress in studying
$\eta$ and $\eta^{\prime}$ mesons from lattice QCD. In Ref. Michael:2013gka
the corresponding mixing has been studied for three values of the lattice
spacing and a large, but still unphysical range of pion mass values in
$N_{f}=2+1+1$ flavour QCD. After extrapolation to the physical pion mass value
excellent agreement to experiment was found. Further lattice results for
$\eta,\eta^{\prime}$ can be found in Refs. Christ:2010dd ; Dudek:2011tt ;
Gregory:2011sg ; Dudek:2013yja , however, all at a single lattice spacing
value and unphysically large pion mass values. In Refs. Fukaya:2015ara ;
Aoki:2017paw  the continuum limit has been studied.

In this proceeding we attempt to study the $\eta^{\prime}$ meson directly at
the physical point, in a first step in $N_{f}=2$ flavour QCD. In $N_{f}=2$
flavour QCD there exists a pion triplet and one flavour singlet, which is
related to the aforementioned anomaly. We will denote it as $\eta_{2}$ meson
to distinguish it from the $\eta^{\prime}$ meson in full QCD, which is only
approximately a flavour eigenstate. The $\eta_{2}$ and the $\eta^{\prime}$
meson have in common that both receive significant fermionic disconnected
contributions to their correlation functions. Their masses are expected to
differ only by $200\ \mathrm{MeV}$, since the strange quark is expected to
contribute significantly less than the light quarks McNeile:2000hf . In
particular, both are expected to have a similar dependence on the light quark
mass. The most recent lattice QCD study of the $\eta_{2}$ meson can be found
in Ref. Jansen:2008wv .

Hence, studying the $\eta_{2}$ meson at the physical point will reveal on the
one hand important qualitative information on the realisation of the anomaly
in QCD. On the other hand it represents a feasibility study for a later
investigation of $\eta$ and $\eta^{\prime}$ in $N_{f}=2+1+1$ QCD at the
physical point finkenrath:2017lat . The results obtained here are also
important prerequisites for an exploratory study of
$\eta_{2}\to\gamma^{\star}\gamma^{\star}$.

Since the $\eta_{2}$ meson is tightly connected to topology, we augment this
study by measuring the topological susceptibility as well. For this purpose we
apply gradient flow techniques Luscher:2010iy  and compare to results
available in the literature.

##  2 Lattice Action and Operators

Table 1: The gauge ensembles used in this study. The labelling of the ensembles follows the notations in Ref. Abdel-Rehim:2015pwa . In addition to the relevant input parameters we give the lattice volume $(L/a)^{3}\times T/a$ and the number of evaluated configurations $N_{\mathrm{conf}}$. ensemble | $\beta$ | $c_{\mathrm{sw}}$ | $a\mu_{\ell}$ | $(L/a)^{3}\times T/a$ | $N_{\mathrm{conf}}^{\eta_{2}}$ | $N_{\mathrm{conf}}^{\chi_{t}}$
---|---|---|---|---|---|---
$cA2.09.48$ | 2.10 | 1.57551 | 0.009 | $48^{3}\times 96$ | $615$ | $1277$
$cA2.30.48$ | 2.10 | 1.57551 | 0.030 | $48^{3}\times 96$ | $352$ | $702$
$cA2.30.24$ | 2.10 | 1.57551 | 0.030 | $24^{3}\times 48$ | $-$ | $703$
$cA2.60.32$ | 2.10 | 1.57551 | 0.060 | $32^{3}\times 64$ | $337$ | $2590$
$cA2.60.24$ | 2.10 | 1.57551 | 0.060 | $24^{3}\times 48$ | $-$ | $1611$

The results presented in this paper are based on the gauge configurations
generated by the ETMC with Wilson clover twisted mass quark action at maximal
twist Frezzotti:2000nk . We employ the Iwasaki gauge action Iwasaki:1985we .
The measurements are performed on up to five $N_{f}=2$ ensembles with pion
mass at its physical value, 240 MeV and 340 MeV, respectively. The lattice
spacing is $a=0.0931(2)\ \mathrm{fm}$ for all three ensembles. In Table 1 we
list the five ensembles with the relevant input parameters, the lattice volume
and the number of configurations. More details about the ensembles can be
found in Ref. Abdel-Rehim:2015pwa .

As a smearing scheme we use the stochastic Laplacian Heavyside (sLapH) method
Peardon:2009gh ; Morningstar:2011ka  for our computation. The details of the
sLapH parameter choices for a set of $N_{f}=2+1+1$ Wilson twisted mass
ensembles are given in Ref. Helmes:2015gla . The parameters for the ensembles
used in this work are the same as those for $N_{f}=2+1+1$ ensembles with the
corresponding lattice volume.

In $N_{f}=2$ flavour QCD there is the neutral pion, corresponding to the
neutral one of the three pions in the triplet, and the $\eta_{2}$, the flavour
singlet pseudo-scalar meson related to the axial $U(1)$ anomaly. The
(hermitian) interpolating operator projected to zero momentum reads

$\mathcal{O}^{s}(t)\ =\
\frac{1}{\sqrt{2}}\sum_{\mathbf{x}}\bar{\psi}(\mathbf{x},t)i\gamma_{5}\,\mathds{1}_{f}\,\psi(\mathbf{x},t)\,.$
(1)

Here, $\mathds{1}_{f}$ is the unit matrix acting in flavour space. From the
operator one builds the correlation functions

$C_{\eta_{2}}(t-t^{\prime})\ =\
\langle\mathcal{O}^{s}(t)\,(\mathcal{O}^{s}(t^{\prime}))^{\dagger}\rangle\,,$
(2)

which allows one to determine the mass $M_{\eta_{2}}$ from its decay in
Euclidean time.

The correlation function in Eq. 2 has both a fermionic connected and a
fermionic disconnected contribution. The connected contribution can be
estimated with high statistical accuracy. It decays asymptotically
exponentially with the mass of the (connected neutral) pion. The disconnected
part is statistically noisy. The signal for the $\eta_{2}$ meson arises from
the difference of connected and disconnected contribution,

$C_{\eta_{2}}(t)\ =\ C^{\mathrm{conn}}(t)\ -\ C^{\mathrm{disc}}\,,$ (3)

where the disconnected contribution has to cancel the pion contribution coming
from the connected part. Since the pion is significantly lighter than the
$\eta_{2}$, the statistical analysis is rather delicate.

The topological susceptibility is a measure for the fluctuations of the
topological charge. It is defined as

$\chi_{t}\ =\ \frac{1}{V}\int\mathrm{d}x\ \int\mathrm{d}y\ \langle
q(x)q(y)\rangle$ (4)

with the volume $V$ and the topological charge density $q(x)$ for which we use
a clover-type discretisation of the field strength tensor, following the
conventions in Ref. Bruno:2014ova . For the topological susceptibility
$\chi_{t}$ we employ the gradient flow as introduced for lattice QCD in Ref.
Luscher:2010iy . This setup has the advantage of yielding a renormalised
topological susceptibility at any finite flow time $t$. Since the renormalised
susceptibility is scale invariant, it becomes independent of $t$ at
sufficiently large flow times. This is indeed what we observe in our
calculation and we choose $t=3t_{0}$ with $\sqrt{t_{0}}=0.1535(12)$ fm from
Ref. Abdel-Rehim:2015pwa  obtained at the physical point.

##  3 Analysis

The data is analysed using a blocked bootstrap procedure with $R=1500$
bootstrap samples. Depending on the ensemble, we have chosen the block size
such that roughly $150$ blocked data points where left. For the analysis of
the topological susceptibility we use $R=1000$ bootstrap samples and block
sizes yielding around $100$ data blocks depending on the ensemble under
consideration. The resulting error is compared to the naive one appropriately
rescaled with the integrated autocorrelation time of $Q^{2}$, and the larger
of the two is chosen as the final error.

###  3.1 Excited State Subtraction

In particular for the $\eta_{2}$ meson, the fermionic disconnected
contributions are very noisy. As a consequence, the signal is lost relatively
early in Euclidean time. For this reason we have in the past applied a method
to subtract excited states Jansen:2008wv ; Michael:2013gka ; Liu:2016cba ,
originally proposed in Ref. Neff:2001zr . Because it works very well, we will
apply it here again for the $\eta_{2}$ meson. It consists of subtracting
excited states from the connected contribution only. This is feasible, because
the connected part – representing a pion correlation function – has a signal
for all Euclidean time values. Therefore, we can fit it at large enough
Euclidean times such that excited states have decayed sufficiently. Next, we
replace the connected correlation function at small times by the fitted
function. Thereafter, the such subtracted connected contribution is summed
according to Eq. 3 to the full $\eta_{2}$ correlation functions.

The underlying assumption is that disconnected contributions are large for the
ground state, i.e. the $\eta_{2}$, but not for excited states. If this
assumption is correct, the effective mass

$M_{\mathrm{eff}}\ =\ -\log\frac{C_{\eta_{2}}(t)}{C_{\eta_{2}}(t+1)}$ (5)

should show a plateau from very early Euclidean times on. We have found in
Refs. Jansen:2008wv ; Michael:2013gka  that this approach works very well for
the $\eta_{2}$ meson in $N_{f}=2$ flavour QCD as well as for $\eta$ and
$\eta^{\prime}$ mesons in $N_{f}=2+1+1$ flavour QCD.

Phenomenologically, the connected contribution gives the pion as the ground
state for large Euclidean times. And, due to symmetries, it cannot couple to
the $\eta_{2}$. The pion contribution needs to be cancelled by the
disconnected contribution, which, therefore, reads

$C_{\eta_{2}}^{\mathrm{disc}}\ \approx\ A\cdot e^{-M_{\pi^{0}}}-B\cdot
e^{-M_{\eta_{2}}t}$

if no higher excited states contribute.

###  3.2 Shifted Correlation Functions

The expected time dependence of the correlation functions considered here
reads as follows

$C(t)=|\langle 0|\mathcal{O}^{s}|0\rangle|^{2}+\sum_{n}\frac{|\langle
0|\mathcal{O}^{s}|n\rangle|^{2}}{2E_{n}}\left(e^{-E_{n}t}+e^{-E_{n}(T-t)}\right)\,,$
(6)

where $n$ labels the states with the corresponding quantum numbers. The first,
time independent term on the right hand side corresponds to the vacuum
expectation value (vev). Using the symmetries of our action one can show that
for $\mathcal{O}^{s}$ the vev must be zero. Potential deviations from this
expectation due to finite statistics and volume can be dealt with by building
the shifted correlation function

$\tilde{C}(t)\ =\ C(t)-C(t+1)\,.$ (7)

The difference cancels any constant contributions in the correlation function,
while also changing the time dependence to be anti-symmetric in time

$\tilde{C}(t)\ \propto\ \left(e^{-E_{n}t}-e^{-E_{n}(T-t)}\right)\,.$ (8)

In addition, there can be finite volume effects due to incorrectly sampled
topological charge sectors. The corresponding formulae for fixed topological
sectors have been worked out in Ref. Aoki:2007ka . They lead to a finite
volume effect in $C_{\eta_{2}}$ constant in Euclidean time of the form
Bali:2014pva

$\propto\ \frac{a^{5}}{T}\left(\chi_{t}-\frac{Q^{2}}{V}\right)$ (9)

proportional to the difference of the topological susceptibility, $\chi_{t}$,
and the expectation value of the squared topological charge $Q^{2}/V$. If
present, such a term will cause the $\eta_{2}$ correlation function to stay
finite at large Euclidean times. Depending on the sign of the coefficient in
front of the finite volume effect, the correlation function may even turn
negative at relatively small Euclidean times. Clearly, a finite volume effect
of this type can be subtracted again using the shifting procedure.

The shifting procedure has the additional advantage of removing statistical
correlations between neighbouring time slices $t/a$ and $t/a+1$. Therefore,
shifting can also be advantageous for the analysis procedure, because the
inverse variance-covariance matrix is often difficult to estimate reliably
Michael:1994sz . Of course, when the correlation between adjacent time slices
can be taken into account correctly, there is no difference in using the
unshifted or the shifted correlation function and both give exactly identical
results (unless there is a constant contribution).

###  3.3 Scale Setting

In order to compare the $\eta_{2}$ masses to older publications available in
the literature, we use the Sommer scale $r_{0}/a$ for the $M_{\eta_{2}}$. The
value $r_{0}/a=5.317(48)$ from Ref. Abdel-Rehim:2015pwa  for the $cA2.09.48$
ensemble at the physical point is used for all ensembles considered here. The
corresponding physical value reads $r_{0}=0.4907(5)\ \mathrm{fm}$.

Since we use the gradient flow to determine the topological susceptibility, it
is natural to use a scale from the gradient flow for $\chi_{t}$. In order to
compare our results with the ones from Ref. Bruno:2014ova  we use the scale
$t_{1}(m_{\pi})$. Its physical value at the physical point reads $t_{1}=0.061\
\mathrm{fm}^{2}$ Bruno:2014ova . Note, however, that $\chi_{t}$ is evaluated
at flow time $t=3t_{0}$ as discussed above.

##  4 Results

Figure 1: Left: $\eta_{2}$ correlation function $C_{\eta_{2}}$ and its shifted
counterpart $\tilde{C}_{\eta_{2}}$ for ensemble cA2.30.48 as a function of
$t/a$. Right: connected, disconnected and full shifted correlation function
$\tilde{C}_{\eta_{2}}$ as a function of $t/a$ with logarithmic $y$-axis.

In the left panel of Figure 1 we show the $\eta_{2}$ correlation function
$C_{\eta_{2}}$ and the corresponding shifted one $\tilde{C}_{\eta_{2}}$ as
functions of the Euclidean time for ensemble cA2.30.48. For $C_{\eta_{2}}$
errors are larger than the typical fluctuations for $t/a>10$, pointing towards
large correlation between separate time slices. Still, the correlation
function is not incompatible with zero at large times. In the shifted
correlation function $\tilde{C}_{\eta_{2}}$ the correlation is largely removed
and the error bars are of the typical fluctuation size. Also
$\tilde{C}_{\eta_{2}}$ is compatible with zero for large times $t/a>15$.

In the right panel of Figure 1 we show the different contributions to the
$\eta_{2}$ shifted correlation function for ensemble cA2.30.48 in a
logarithmic plot. While the signal for the connected part (red circles) is
excellent up to $t/a=48$, the disconnected contribution (blue squares) starts
to become noisy around $t/a=15$. We observe that connected and disconnected
contributions are very similar in magnitude from around $t/a=10$ on. The full
$\eta_{2}$ correlator (green diamonds) is given by the difference of connected
and disconnected parts. Consequently, the signal deteriorates around $t/a=10$.

From the left panel of Figure 1 it is also clear that the ground state energy
and amplitude of the connected contribution (in twisted mass the connected
neutral pion) can be extracted reliably. The connected contribution is then
replaced by only this ground state and the full disconnected part is
subtracted from it. As a result we obtain plateaus from $t/a=3$ on.

Table 2: Results for the $\eta_{2}$ meson in lattice units and for the topological susceptibility in units of $t_{1}^{2}$ for the ensembles investigated here. Ensemble | $aM_{\eta_{2}}$ | $10^{3}t_{1}^{2}\ \chi_{t}$
---|---|---
$cA2.09.48$ | $0.358(9)$ | $0.53(4)$
$cA2.30.48$ | $0.372(8)$ | $0.63(6)$
$cA2.30.24$ | — | $0.52(6)$
$cA2.60.32$ | $0.374(7)$ | $0.91(5)$
$cA2.60.24$ | — | $0.91(8)$

The such extracted energy levels in lattice units are compiled in Table 2. The
$\eta_{2}$ masses are shown in units of the Sommer parameter $r_{0}$ in Figure
2 as a function of the squared pion mass $(r_{0}M_{\pi^{\pm}})^{2}$. The red
circles represent the results presented here. In addition we show results for
Wilson twisted mass fermions with $N_{f}=2$ dynamical quark flavours without
clover term taken from Ref. Jansen:2008wv . The blue squares correspond to a
lattice spacing of about $0.09\ \mathrm{fm}$, while the green diamonds
correspond to about $0.08\ \mathrm{fm}$. These results from Ref. Jansen:2008wv
have larger statistical errors and the smallest pion mass corresponds to about
$300\ \mathrm{MeV}$. However, the agreement with the results presented here is
good. In particular, the almost constant extrapolation to the physical point
is confirmed. An estimate of lattice artefacts is difficult, but they seem to
be not larger than our statistical uncertainties.

Figure 2: Compilation of $N_{f}=2$ Wilson twisted mass results with and
without clover term for the $\eta_{2}$ meson. The twisted mass (tm) results
have been taken from Ref. Jansen:2008wv .

Figure 3: Topological susceptibility $\chi_{t}$ as a function of the squared
pion mass, both in appropriate units of $t_{1}$. The solid line with shaded
error band indicates a fit to the data according to Eq. 10.

In Table 2 we have also compiled our results for the topological
susceptibility in units of $t_{1}$. They are displayed in Figure 3 as a
function of $t_{1}M_{\pi^{\pm}}^{2}$. In leading order Wilson ChPT one expects
the following dependence of $\chi_{t}$ on the lattice spacing and the pion
mass Bruno:2014ova

$t_{1}^{2}\,\chi_{t}\ =\ \frac{1}{8}\,t_{1}^{2}\,f_{\pi}^{2}\,M_{\pi^{0}}^{2}\
+\ a^{2}\,\frac{c_{2}}{t_{1}}\,.$ (10)

Apart from the ensemble with too small volume cA2.30.24, our data is nicely
compatible with this expectation: the solid line in Figure 3 represents a fit
of the function

$g(M_{\pi^{\pm}}^{2})\ =\ c_{1}t_{1}M_{\pi^{\pm}}^{2}\ +\
a^{2}\,\frac{c_{2}}{t_{1}}$

to the data with fit parameters $c_{1}$ and $c_{2}$. Note that we use the
charged pion mass, because charged and neutral pion masses are degenerate
within errors Abdel-Rehim:2015pwa . The best fit parameter for $c_{1}$ is
compatible with $t_{1}f_{\pi}^{2}/8$. Note that ensemble cA2.30.24 has a very
small volume explaining the outlier in Figure 3.

The fitted value for $c_{2}$ can be compared to the results of Ref.
Bruno:2014ova . We obtain $c_{2}=3.3(3)\cdot 10^{-3}$, while Bruno et al. have
$c_{2}=5.1(7)\cdot 10^{-3}$.

##  5 Summary

In this proceeding contribution we have presented results for the $\eta_{2}$
meson related to the axial anomaly and the topological susceptibility in two-
flavour QCD. The results have been obtained using $N_{f}=2$ lattice QCD
ensembles generated by ETMC with the Wilson twisted clover discretisation
Abdel-Rehim:2015pwa . Pion mass values reach from its physical value up to
$340\ \mathrm{MeV}$ at a single lattice spacing value of $a=0.0931(2)\
\mathrm{fm}$. For the $\eta_{2}$ we could confirm the almost constant
extrapolation in $M_{\pi^{\pm}}^{2}$ towards the physical point. Errors are
significantly reduced compared to previous calculations. Lattice artefacts
seem to be not larger than our statistical uncertainty.

The topological susceptibility has been computed using the gradient flow. As
expected, $\chi_{t}$ is proportional to $M_{\pi^{\pm}}^{2}$ up to a lattice
artefact independent of $M_{\pi^{\pm}}$. Even if we are not able to finally
confirm this with a single lattice spacing at hand, this constant term should
be of $\mathcal{O}(a^{2})$. The size of this artefact appears to be
significantly smaller than what is observed with Wilson clover fermions in
Ref. Bruno:2014ova .

We thank the members of ETMC for the most enjoyable collaboration. This
project was funded by the DFG as a project in the Sino-German CRC110. The
computer time for this project was made available to us by the John von
Neumann-Institute for Computing (NIC) on the Juqueen system in Jülich and by
the ITP Bern on their HPC clusters. This work was granted access to the HPC
resources of CINES and IDRIS under the allocation 52271 made by GENCI.

## References

  * (1) F. Jegerlehner, EPJ Web Conf. 118, 01016 (2016), 1511.04473
  * (2) C. Michael, K. Ottnad, C. Urbach (ETM), Phys.Rev.Lett. 111, 181602 (2013), 1310.1207
  * (3) N. Christ, C. Dawson, T. Izubuchi, C. Jung, Q. Liu et al., Phys.Rev.Lett. 105, 241601 (2010), 1002.2999
  * (4) J.J. Dudek, R.G. Edwards, B. Joo, M.J. Peardon, D.G. Richards et al., Phys.Rev. D83, 111502 (2011), 1102.4299
  * (5) E.B. Gregory, A.C. Irving, C.M. Richards, C. McNeile (UKQCD), Phys.Rev. D86, 014504 (2012), 1112.4384
  * (6) J.J. Dudek, R.G. Edwards, P. Guo, C.E. Thomas (Hadron Spectrum), Phys. Rev. D88, 094505 (2013), 1309.2608
  * (7) H. Fukaya, S. Aoki, G. Cossu, S. Hashimoto, T. Kaneko, J. Noaki (JLQCD), Phys. Rev. D92, 111501 (2015), 1509.00944
  * (8) S. Aoki, G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko (JLQCD) (2017), 1705.10906
  * (9) C. McNeile, C. Michael (UKQCD), Phys.Lett. B491, 123 (2000), hep-lat/0006020
  * (10) K. Jansen, C. Michael, C. Urbach (ETM), Eur.Phys.J. C58, 261 (2008), 0804.3871
  * (11) J. Finkenrath et al. (ETM), this conference (2017)
  * (12) M. LÃscher, JHEP 08, 071 (2010), [Erratum: JHEP03,092(2014)], 1006.4518
  * (13) A. Abdel-Rehim et al. (ETM) (2015), 1507.05068
  * (14) R. Frezzotti, P.A. Grassi, S. Sint, P. Weisz (ALPHA), JHEP 08, 058 (2001), hep-lat/0101001
  * (15) Y. Iwasaki, Nucl. Phys. B258, 141 (1985)
  * (16) M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R.G. Edwards, B. Joo, H.W. Lin, D.G. Richards, K.J. Juge (Hadron Spectrum), Phys. Rev. D80, 054506 (2009), 0905.2160
  * (17) C. Morningstar, J. Bulava, J. Foley, K.J. Juge, D. Lenkner, M. Peardon, C.H. Wong, Phys. Rev. D83, 114505 (2011), 1104.3870
  * (18) C. Helmes, C. Jost, B. Knippschild, C. Liu, J. Liu, L. Liu, C. Urbach, M. Ueding, Z. Wang, M. Werner (ETM), JHEP 09, 109 (2015), 1506.00408
  * (19) M. Bruno, S. Schaefer, R. Sommer (ALPHA), JHEP 08, 150 (2014), 1406.5363
  * (20) L. Liu et al., Phys. Rev. D96, 054516 (2017), 1612.02061
  * (21) H. Neff, N. Eicker, T. Lippert, J.W. Negele, K. Schilling, Phys.Rev. D64, 114509 (2001), hep-lat/0106016
  * (22) S. Aoki, H. Fukaya, S. Hashimoto, T. Onogi, Phys. Rev. D76, 054508 (2007), 0707.0396
  * (23) G.S. Bali, S. Collins, S. DÃrr, I. Kanamori, Phys. Rev. D91, 014503 (2015), 1406.5449
  * (24) C. Michael, A. McKerrell, Phys. Rev. D51, 3745 (1995), hep-lat/9412087