date+extrayear

# Deep Grey-Box Modeling With Adaptive Data-Driven Models
Toward Trustworthy Estimation of Theory-Driven Models

oya Takeishi<EMAIL_ADDRESS>
Geneva School of Business Administration
University of Applied Sciences and Arts Western Switzerland (HES-SO) exandros
Kalousis<EMAIL_ADDRESS>
Geneva School of Business Administration
University of Applied Sciences and Arts Western Switzerland (HES-SO)

###### Abstract

The combination of deep neural nets and theory-driven models, which we call
_deep grey-box modeling_ , can be inherently interpretable to some extent
thanks to the theory backbone. Deep grey-box models are usually learned with a
regularized risk minimization to prevent a theory-driven part from being
overwritten and ignored by a deep neural net. However, an estimation of the
theory-driven part obtained by uncritically optimizing a regularizer can
hardly be trustworthy when we are not sure what regularizer is suitable for
the given data, which may harm the interpretability. Toward a trustworthy
estimation of the theory-driven part, we should analyze regularizers’ behavior
to compare different candidates and to justify a specific choice. In this
paper, we present a framework that enables us to analyze a regularizer’s
behavior empirically with a slight change in the neural net’s architecture and
the training objective.

##  1 INTRODUCTION

Grey-box modeling in general refers to the combination of theory-driven
structures and data-driven components [17, see, e.g.,]. In this paper, we are
interested in combining theory-driven models such as mathematical models of
physical phenomena, and data-driven, machine-learning models. For example, for
regressing from $x$ to $y$, we are interested in models like

$y=f_{\mathrm{T}}(x;\theta_{\mathrm{T}})+f_{\mathrm{D}}(x;\theta_{D})+e,$

where $f_{\mathrm{T}}$ and $f_{\mathrm{D}}$ denote theory- and data-driven
models parameterized by $\theta_{\mathrm{T}}$ and $\theta_{\mathrm{D}}$,
respectively, and $e$ is noise. A more general model will appear in Section 2
and thereafter. We discuss how we should (or should not) cast the estimation
problem of a grey-box model’s parameters.

We argue over the interpretability of grey-box models. They can be
interpretable because a part of the model lies on a theory backbone and often
has a small number of parameters. They are valuable tools to glimpse insights
from data on which our theory is essentially incomplete. However, as we
discuss later, interpretability is not a free lunch, and we need to pay
attention to how we perform the estimation to secure the trustworthiness of
the interpretation.

We are interested in cases where the data-driven model, $f_{\mathrm{D}}$, is a
deep neural network. In such a case, care must be taken for the theory-based
model, $f_{\mathrm{T}}$, not to be overwritten and ignored by $f_{\mathrm{D}}$
due to the expressive power of the latter. Such models have been studied
recently [23, 18, 11, 21]. Deep grey-box models are typically learned by
solving optimization problems like

$\underset{\theta_{\mathrm{T}},\theta_{\mathrm{D}}}{\text{minimize}}~{}\mathcal{L}+\lambda\mathcal{R},$

where $\mathcal{L}$ is a prediction loss, and $\mathcal{R}$ is a regularizer
needed for preventing $f_{\mathrm{T}}$ from being ignored. Although this
optimization may result in a model with good prediction performance, we cannot
judge if such an estimator is worth being interpreted. There can be multiple
parameter values that achieve similar performance, and the optimization’s
solution tells us nothing about the property of $\mathcal{R}$, so it provides
only insufficient information to justify a choice of $\theta_{\mathrm{T}}$ to
be interpreted. For the same reason, we cannot compare different
$\mathcal{R}$s solely based on the optimization’s solution. Instead, we should
analyze the behavior of $\mathcal{R}$, at least empirically, toward obtaining
a trustworthy estimation and its interpretation.

Our idea is simple: We withhold the estimation of $\theta_{\mathrm{T}}$, at
least in the first trial of data analysis. To this end, we learn
$f_{\mathrm{D}}$ adaptively with realizations of $\theta_{\mathrm{T}}$,
essentially marginalizing out $\theta_{\mathrm{T}}$. Such a slight change in
the formulation allows us to analyze $\mathcal{R}$ empirically by examining
its landscape without any need for re-training. In this paper, we take up the
aforementioned argument about the estimation’s trustworthiness for discussion,
formulate the above idea, and conduct an empirical investigation of its
effectiveness.

A useful byproduct of the proposed formulation is that the optimization of
$\mathcal{L}$ and $\mathcal{R}$ can now be decoupled. It allows us to use
different optimizers for the two objectives, as well as to use unlabeled test
data to optimize $\mathcal{R}$. Moreover, we found that such a decoupled
optimization makes the optimization much less sensitive to the hyperparameter,
$\lambda$.

##  2 PRELIMINARY

###  2.1 Definition

We define _deep grey-box models_ as compositions of theory-driven models and
data-driven models, with the latter being deep neural networks. For the sake
of discussion, we suppose regression problems where $y\in\mathcal{Y}$ is to be
predicted from $x\in\mathcal{X}$, though the extension to other problems is
straightforward. We denote such a model in general by

$y=\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x),$ (1)

where $\mathcal{C}$ is a functional that takes the two types of functions and
an input variable as arguments. The two functions, $f_{\mathrm{T}}$ and
$f_{\mathrm{D}}$, are a theory-driven model and a deep neural network, with
unknown parameters $\theta_{\mathrm{T}}\in\Theta_{\mathrm{T}}$ and
$\theta_{\mathrm{D}}\in\Theta_{\mathrm{D}}$, respectively. We may write
$f_{\mathrm{T}}(x;\theta_{T})$ to manifest $f_{\mathrm{T}}$’s dependency on
$\theta_{\mathrm{T}}$ or write simply $f_{\mathrm{T}}(x)$ though it still
depends on $\theta_{\mathrm{T}}$ (and analogously for $f_{\mathrm{D}}$ and
$\theta_{\mathrm{D}}$). Note that not only $f_{\mathrm{D}}$ but also
$f_{\mathrm{T}}$ may have unknown parameters to be inferred. We usually expect
$\dim\theta_{\mathrm{T}}\ll\dim\theta_{\mathrm{D}}$. The functional,
$\mathcal{C}$, evaluates $f_{\mathrm{T}}$ and $f_{\mathrm{D}}$ with a given
$(\theta_{\mathrm{T}},\theta_{D},x)$ and then mixes up their outputs to give
the final output of the model.

We try to keep the generality of $\mathcal{C}$; it may include general
function compositions and their arbitrary transformations:

$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=\operatorname{SomeTransformation}[f_{\mathrm{D}}(f_{\mathrm{T}}(x),x)].$

Meanwhile, one of the most prevailing forms of $\mathcal{C}$ in the literature
is the additive grey-box ODEs like:

$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=\operatorname{ODESolve}[\dot{s}=f_{\mathrm{T}}(s)+f_{\mathrm{D}}(s)\mid
s_{0}=x],$

where $s$ is the state variable of the dynamics, $s_{0}$ is the initial
condition, and $\operatorname{ODESolve}$ denotes an operation that numerically
solves initial value problems. Such _grey-box_ (ordinary or partial)
differential equations have been studied by researchers such as [14, 23, 18,
11].

We are particularly interested in cases where $\mathcal{C}$ inherits the
expressive power of $f_{\mathrm{D}}$ as a function approximator. This is _not_
the case, for example, when $\mathcal{C}$ only contains compositions such as
$f_{\mathrm{T}}(f_{\mathrm{D}}(x),x)$, that is, $f_{\mathrm{T}}$ is “outside”
$f_{\mathrm{D}}$ [12, 2, 15, e.g.,]. Estimation of such models is less
challenging because $f_{\mathrm{T}}$ cannot be ignored by construction. In
contrast, we address more difficult cases where $f_{\mathrm{T}}$ is “inside”
$f_{\mathrm{D}}$, for which we should be careful so that $f_{\mathrm{T}}$ is
not overwritten and ignored by $f_{\mathrm{D}}$. We put the following
assumptions on the model:

######  Assumption 1.

$f_{\mathrm{D}}\colon\mathcal{X}\\!\to\\!\mathcal{Y}$ is a universal function
approximator; for any $\epsilon>0$ and a continuous function
$g\colon\mathcal{X}\to\mathcal{Y}$, there exists
$\theta_{\mathrm{D}}\in\Theta_{\mathrm{D}}$ satisfying $\sup_{x\in
S_{X}}\|f_{\mathrm{D}}(x;\theta_{\mathrm{D}})-g(x)\|<\epsilon$, where
$S_{X}\subset\mathcal{X}$ is some compact set.

######  Assumption 2.

$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};\cdot)\colon\mathcal{X}\to\mathcal{Y}$
is also a universal function approximator; that is, for any $\epsilon>0$,
$\theta_{\mathrm{T}}\in\Theta_{\mathrm{T}}$, and a continuous function
$g^{\prime}\colon\mathcal{X}\to\mathcal{Y}$, there exists
$\theta_{\mathrm{D}}\in\Theta_{\mathrm{D}}$ satisfying $\sup_{x\in
S_{X}}\|\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)-g^{\prime}(x)\|<\epsilon$.

######  Remark 1.

We assume the universal approximation property just for rigorously construct
the discussion. Even without the universal approximation property, as long as
$f_{\mathrm{D}}$ and $\mathcal{C}$ are much more expressive than
$f_{\mathrm{T}}$, discussions below would approximately hold in practice.

######  Assumption 3.

$f_{\mathrm{T}}$ and $f_{\mathrm{D}}$ are Lipschitz continuous with regard to
$(x,\theta_{\mathrm{T}})$ and $(x,\theta_{\mathrm{D}})$, respectively.

###  2.2 Why Grey-box?

Deep grey-box models are powerful function approximators with a certain level
of inherent interpretability owing to $f_{T}$, a human-understandable model
with a theory as a backbone. A typical use case would be to estimate a grey-
box model on data on which our theory is essentially incomplete and inspect
the estimated model to glimpse insights, e.g., when the incomplete theory is
correct or not, how the missing part approximated by $f_{\mathrm{D}}$ behaves,
and so on.

Deep grey-box models can also be advantageous in generalization capability and
robustness to extrapolation, as reported empirically so far [11, 23, 18, 21].
It is natural to expect such improvements because the presence of
$f_{\mathrm{T}}$ in the model would reduce the sample complexity of the
learning problem, and $f_{\mathrm{T}}$ is supposed to work well in the out-of-
data regime (in other words, it is a requirement for a model to be regarded as
theory-driven). However, rigorous analysis of generalization is challenging
for models involving deep neural nets. Anyway, we do not touch on such
performance aspects of deep grey-box models given the previous studies, so the
comparison to non-grey-box models is out of the paper’s scope.

###  2.3 Empirical Risk Minimization Cannot Select $\theta_{T}$

There is a natural consequence of deep grey-box modeling; the theory-driven
model’s parameter, $\theta_{\mathrm{T}}$, cannot be chosen solely by
minimizing an empirical risk of prediction. For example, suppose we learn
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=f_{\mathrm{T}}(x)+f_{\mathrm{D}}(x)$
by minimizing the mean squared error,
$\mathcal{L}=\|y-(f_{\mathrm{T}}(x)+f_{\mathrm{D}}(x))\|_{2}^{2}$. The
empirical risk can be minimized to a similar extent for _any_
$\theta_{\mathrm{T}}\in\Theta_{\mathrm{T}}$ because $f_{\mathrm{D}}$, a deep
neural net, can, it alone, approximate any function on the training set (as
assumed in 1) and thus also the function $y-f_{\mathrm{T}}(x)$. We formally
state this fact as follows:

######  Proposition 1.

Let $S=\\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\\}$ be a training set. Let
$\mathcal{L}_{(x,y)}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})$ be a Lipschitz
continuous loss function between the prediction (i.e., the value of
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)$) and the target (i.e., $y$).
Let $\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})=\sum_{(x,y)\in
S}\mathcal{L}_{(x,y)}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})$ be the
empirical risk on the training set. Suppose that Assumptions 1–3 hold. Then,
for any $\epsilon^{\prime}>0$, $\theta_{\mathrm{D}}\in\Theta_{\mathrm{D}}$,
and $\theta_{\mathrm{T}},\theta^{\prime}_{\mathrm{T}}\in\Theta_{\mathrm{T}}$
where $\theta_{\mathrm{T}}\neq\theta^{\prime}_{\mathrm{T}}$, there exists
$\theta^{\prime}_{\mathrm{D}}\in\Theta_{\mathrm{D}}$ that satisfies

$\left|\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})-\mathcal{L}_{S}(\theta^{\prime}_{\mathrm{T}},\theta^{\prime}_{\mathrm{D}})\right|<\epsilon^{\prime}.$
(2)

###### Proof.

From the assumptions, for any $\epsilon>0$,
$\theta_{\mathrm{D}}\in\Theta_{D}$, and
$\theta_{\mathrm{T}}\neq\theta^{\prime}_{\mathrm{T}}\in\Theta_{\mathrm{T}}$,
there exists $\theta^{\prime}_{\mathrm{D}}\in\Theta_{D}$ that satisfies
$\sup_{x\in\\{x_{1},\dots,x_{n}\\}}\|\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)-\mathcal{C}(f^{\prime}_{\mathrm{T}},f^{\prime}_{\mathrm{D}};x)\|<\epsilon$,
where $f^{\prime}_{i}$ is parameterized by $\theta^{\prime}_{i}$ for
$i={\mathrm{T}},{\mathrm{D}}$. Since $\mathcal{L}$ is Lipschitz continuous,
$\sup_{x,y\in
S}|\mathcal{L}_{(x,y)}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})-\mathcal{L}_{(x,y)}(\theta^{\prime}_{\mathrm{T}},\theta^{\prime}_{\mathrm{D}})|<K\epsilon$
with where $K$ is $\mathcal{L}$’s Lipschitz constant. Therefore, with
$\epsilon^{\prime}\coloneqq|S|K\epsilon$,
$|\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})-\mathcal{L}_{S}(\theta^{\prime}_{\mathrm{T}},\theta^{\prime}_{\mathrm{D}})|<\epsilon^{\prime}$.
∎

###  2.4 Regularized Risk Minimization

Proposition 1 states that any $\theta_{\mathrm{T}}\in\Theta_{\mathrm{T}}$ can
be equally likely solely under the empirical risk. It necessitates
regularizing the problem; we should optimize
$\mathcal{L}_{S}+\lambda\mathcal{R}$ instead, where $\lambda\geq 0$ is a
regularization hyperparameter, and $\mathcal{R}$ is some regularizer that
should reflect our inductive biases on how we should combine the theory- and
data-driven models. Let us, for example, consider the linear combination case,
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=f_{\mathrm{T}}(x)+f_{\mathrm{D}}(x)$.
One of the common ways of thinking is that $f_{\mathrm{T}}$ should as
accurately explain the $x$–$y$ relation as possible, and $f_{\mathrm{D}}$
should have the least possible effect. This idea can be operationalized by
defining $\mathcal{R}=\|f_{\mathrm{D}}\|$, where the norm is the function’s
norm. Though such an $\mathcal{R}$ has been a popular choice, it is not the
only possibility. For example, when one wants the two models’ output to be
uncorrelated, one can use $\mathcal{R}=|\langle
f_{\mathrm{T}}(x),f_{\mathrm{D}}(x)\rangle|$.111Suggesting specific
$\mathcal{R}$ for each application or in general is out of the scope of this
paper; on contrary, our proposal in Section 3 is for cases where we cannot
specify $\mathcal{R}$ _a priori_. We assume that $\mathcal{R}$ depends only on
$x$. It is natural because the role of $\mathcal{R}$ is not to fit the $x$–$y$
relation. We will recall this assumption, if necessary, by writing
$\mathcal{R}_{S_{X}}$, where $S_{X}=\\{x_{1},\dots,x_{n}\\}$ is the extract of
$x$s from $S$. We do not suppose, at least explicitly, any more specifications
of $\mathcal{R}$ than this assumption.

The regularized estimation problem can be cast as follows:

##### Inductive Learning

The simplest formulation is

$\theta_{\mathrm{T}}^{*},\theta_{\mathrm{D}}^{*}=\arg\min_{\theta_{\mathrm{T}},\theta_{\mathrm{D}}}\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})+\lambda\mathcal{R}_{S_{X}}(\theta_{\mathrm{T}},\theta_{\mathrm{D}}).$
(3)

In this formulation, not only $\mathcal{L}$ but also $\mathcal{R}$ suffers a
generalization gap. Also, $\lambda$ needs to be tuned somehow.

##### Transductive Learning

Since we assume that $\mathcal{R}$ only depends on $x$, it is reasonable to
mention _transductive learning_ [8]. Let $S^{\prime}_{X}$ be some set of $x$
that may include $S_{X}$ as a subset. The idea is to minimize the unsupervised
part of the objective, not on the training data $S_{X}$ but rather on
$S^{\prime}_{X}$;

$\theta_{\mathrm{T}}^{*},\theta_{\mathrm{D}}^{*}=\arg\min_{\theta_{\mathrm{T}},\theta_{\mathrm{D}}}\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})+\lambda\mathcal{R}_{S^{\prime}_{X}}(\theta_{\mathrm{T}},\theta_{\mathrm{D}}).$
(4)

The generalization gap disappears for $\mathcal{R}$ when (a subset of)
$S^{\prime}_{X}$ is the test set, but still $\lambda$ needs to be tuned.

##  3 TOWARD TRUSTWORTHY ESTIMATION

###  3.1 Challenges of Deep Grey-box Model Estimation

Deep grey-box models have been studied mainly in terms of empirical
generalization and extrapolation capability [23, 18, 11, 21]. However, when
the model’s interpretation is concerned, the prediction performance does not
speak a lot; there can be multiple parameter values that perform similarly
(cf. _Rashomon_ sets), and we cannot judge which one we should interpret. The
solution of the optimization in Eq. 3 or (4) tells us nothing about the
analytical property of $\mathcal{R}$, so we can hardly understand the full
picture of how the optimization selects $\theta_{\mathrm{T}}$. Instead of
uncritically optimizing the regularizer, $\mathcal{R}$, we should know the
property of $\mathcal{R}$ in order to gain more information to explain the
choice of $\theta_{\mathrm{T}}$ to be interpreted. We would contrast the
situation with, for example, the estimators of linear regression models, which
have been extensively analyzed and thus are trustworthy in some sense. We do
not suggest analyzing our $\mathcal{R}$s analytically as it is too problem-
dependent, but analyzing them at least empirically would help us make
$\theta_{\mathrm{T}}$’s estimation more trustworthy.

The challenge due to not knowing $\mathcal{R}$’s property stands out more when
we do not know what $\mathcal{R}$ is suitable for the given data and need to
compare different candidate $\mathcal{R}$s, which is often the case as we do
not know the whole data-generating process. The point estimation via Eq. 3 or
(4) would not tell much about the goodness of $\mathcal{R}$, since different
$\mathcal{R}$s could achieve similar prediction performance. This viewpoint
also supports the need for analyzing $\mathcal{R}$ at least empirically for
gaining information to compare different $\mathcal{R}$s.

Another challenge, yet more technical, is the choice of the regularization
hyperparameter, $\lambda$. It can be tricky because, in Eq. 3 or (4), it
controls two things at the same time: “which $\theta_{\mathrm{T}}$ should be
selected” and “how much $f_{\mathrm{D}}$ should be regularized.” They are
different problems if interrelated, and thus decoupling them would be
beneficial.

###  3.2 Proposed Formulation

As we argued above, analyzing $\mathcal{R}$ empirically can be a useful first
step toward a trustworthy estimation of $\theta_{\mathrm{T}}$. More
specifically, we aim to explore the landscape of $\mathcal{R}$, that is, to
evaluate the values of $\mathcal{R}$ for different $\theta_{\mathrm{T}}$s. The
naïve way to do so is to re-run the optimization in Eq. 3 or (4) for many
times with different values of $\theta_{\mathrm{T}}$ fixed at each time, but
it is inefficient even for a low-dimensional $\theta_{\mathrm{T}}$, and
moreover, it does not allow to predict with a new value of
$\theta_{\mathrm{T}}$ that was not exactly tried during such many training
runs. Instead, we suggest an approach to exploring $\mathcal{R}$’s landscape
without the need for re-running the optimization.

At the core of our suggestion is “marginalizing out” $\theta_{\mathrm{T}}$ in
the training phase. To this end, we slightly modify $f_{\mathrm{D}}$ so that
it works adaptively with different values of $\theta_{\mathrm{T}}$ and then
minimize the objective, taking its expectation with regard to
$\theta_{\mathrm{T}}$. It should result in a model that can predict equally
well given whatever $\theta_{\mathrm{T}}$ in some feasible region. We first
estimate only $\theta_{\mathrm{D}}$ while leaving $\theta_{\mathrm{T}}$
undetermined; thus the idea is similar to probabilistic inference where some
variables are marginalized out. The remainder of this section explains each
step of the proposed formulation. Figure 1 depicts the overview of the
formulation.

Figure 1: Overview of the proposed formulation.

####  3.2.1 Step 1: Modification of $f_{\mathrm{D}}$

We add two specifications of $f_{\mathrm{D}}$ to be used, without much loss of
generality. First, we suppose that $f_{\mathrm{D}}$ takes as arguments not
only $x$ but also $\theta_{\mathrm{T}}$ (and possibly the value of
$f_{\mathrm{T}}(x;\theta_{\mathrm{T}})$). Second, we define $f_{\mathrm{D}}$’s
architecture such that $f_{\mathrm{D}}$ can adapt to different values of
$\theta_{\mathrm{T}}$ given as an argument. We can realize them with general
techniques of conditional modeling; for example, by concatenating
$\theta_{\mathrm{T}}$ to the original input; or by using hypernetworks that
transform $\theta_{\mathrm{T}}$ to (a part of) $\theta_{\mathrm{D}}$. In the
experiments in Section 5, we simply concatenate $\theta_{\mathrm{T}}$ and
$f_{\mathrm{T}}(x;\theta_{\mathrm{T}})$ to $x$ and feed them to
$f_{\mathrm{D}}$ either with an extended input layer or with additional
channels. We may emphasize the dependency of $f_{\mathrm{D}}$ on
$\theta_{\mathrm{T}}$ (and thus on $f_{\mathrm{T}}$) by denoting it as
$f_{\mathrm{D}}(f_{\mathrm{T}})$ in what follows.

####  3.2.2 Step 2: Estimation of $\theta_{\mathrm{D}}$

We estimate $\theta_{\mathrm{D}}$ by optimizing the expectation of the
objective with regard to $p(\theta_{\mathrm{T}})$, some prior distribution of
$\theta_{\mathrm{T}}$.

$\theta_{\mathrm{D}}^{*}=\arg\min_{\theta_{\mathrm{D}}}\mathbb{E}_{p(\theta_{\mathrm{T}})}\\!\left[\mathcal{L}_{S}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})\\!+\\!\lambda\mathcal{R}_{S_{X}}(\theta_{\mathrm{T}},\theta_{\mathrm{D}})\right],$
(5)

where $\mathcal{L}_{S}$ and $\mathcal{R}_{S_{X}}$ are computed with regard to
the deep grey-box model
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}}(f_{\mathrm{T}});x)$ with the
adaptive data-driven model, $f_{\mathrm{D}}(f_{\mathrm{T}})$. We estimate the
objective in Eq. 5 with $m$ samples of $\theta_{\mathrm{T}}$ drawn from
$p(\theta_{\mathrm{T}})$. We empirically found that $m=1$ was sufficient in
our experiment.

The $\lambda\mathcal{R}_{S_{X}}$ term in Eq. 5 is present just for generality.
When one has multiple candidates of $\mathcal{R}$, which is one of our
motivating situations, it is unclear which $\mathcal{R}$ we should use in Eq.
5. We suggest two options: One is to set $\mathcal{R}$ as the sum of all the
candidates and use a small value of $\lambda$. This is reasonable because it
prevents the candidate $\mathcal{R}$s from getting overly large but still does
not strongly impose any one of them. Another option is to set $\lambda=0$,
with which we no longer need to set any $\mathcal{R}$. This is also a
reasonable choice because the minimization of $\mathcal{R}$, if it matters,
can happen later in the post-estimation phase as we will see.

We suppose that we have some $p(\theta_{\mathrm{T}})$ from which we can draw
$\theta_{\mathrm{T}}$s. It is typically available from domain knowledge
concerning the theory-based model. For example, it can be the uniform
distribution on a plausible range of the parameters. It is also technically
possible to let $p(\theta_{\mathrm{T}})$ have unknown parameters and estimate
them from data with the reparametrization trick.

####  3.2.3 Step 3: Post-Estimation Options

The optimization in Eq. 5 leaves $\theta_{\mathrm{T}}$ undetermined, while the
model works with any different values of $\theta_{\mathrm{T}}$ within a region
having a reasonable mass of $p(\theta_{\mathrm{T}})$. It allows us to take the
following options in the post-estimation phase.

##### Explore the Landscape of $\mathcal{R}$

We can compute the values of $\mathcal{R}$ for different
$\theta_{\mathrm{T}}$s, only with the cost of the forward evaluation of the
functions. If $\theta_{\mathrm{T}}$ is up to two-dimensional, we can directly
draw the landscape of $\mathcal{R}$. If $\dim\theta_{\mathrm{T}}>2$ but
remains moderate (say $\lessapprox 20$), we can watch the variations of
$\mathcal{R}$ by varying each pair of $\theta_{\mathrm{T}}$’s elements while
fixing the others at some reference values. The same discussion applies to the
analysis of the supervised loss, $\mathcal{L}$.

How should we utilize the landscape of $\mathcal{R}$? Although it is just up
to a user’s policy and belief in each application, as a general practice, we
suggest utilizing it for assessing the reliability of potential estimation
based on the $\mathcal{R}$. For example, if an $\mathcal{R}$ does not have
clear extrema along some axes of $\theta_{\mathrm{T}}$, it implies that those
elements of $\theta_{\mathrm{T}}$ are not quite identifiable under such an
$\mathcal{R}$. It can contribute to the estimation’s trustworthiness, e.g., by
withholding interpretation about some parameters. Another general usage is to
compare different candidates of $\mathcal{R}$. The landscapes can give an
intuition about the nature of each $\mathcal{R}$, which is useful for a user
to choose one (or more) out of the candidates; we will see concrete use cases
later in numerical experiments.

Models with high-dimensional $\theta_{\mathrm{T}}$ basically remain an open
challenge; analyzing and visualizing a high-dimensional parameter space are
very challenging in general. One of the options is to get an overview of the
landscape via dimensionality reduction of the parameter space by techniques
such as random projection and principal component analysis. With that being
said, even if $\theta_{\mathrm{T}}$ is high-dimensional, we can still benefit
from the proposed formulation for the point estimation or the posterior
inference discussed below.

##### Point Estimation of $\theta_{\mathrm{T}}$ (optional)

After analyzing $\mathcal{R}$s and choosing one to use, a user may want to
select $\theta_{\mathrm{T}}$ as

$\theta_{\mathrm{T}}^{*}=\arg\min_{\theta_{\mathrm{T}}}\mathcal{R}_{S^{\prime}_{X}}(\theta_{\mathrm{T}},\theta_{\mathrm{D}}^{*}).$
(6)

Depending on the application, $S^{\prime}_{X}$ may be a singleton of a test
sample, a set of test samples, the union of training and test sets, or a set
of grid points on $\mathcal{X}$. After choosing a specific
$\theta_{\mathrm{T}}$, one does not need to re-run the optimization in Eq. 5
because the estimated $f_{\mathrm{D}}$ works adaptively to
$\theta_{\mathrm{T}}$. With that being said, it is also an option to drop
$f_{\mathrm{D}}$’s dependency on $\theta_{\mathrm{T}}$ and run the
optimization in Eq. 3 or (4) only with regard to $\theta_{\mathrm{D}}$.

Although the end result of Eq. 6 (i.e., a point estimation) has the same form
with that of existing frameworks in Eq. 3 or (4), our framework can benefit
from the decoupled nature of the optimization, that is, $\theta_{\mathrm{D}}$
and $\theta_{\mathrm{T}}$ are optimized individually in Eq. 5 and Eq. 6,
respectively. In Eq. 5, $\mathcal{R}$ may be used solely for regularizing the
behavior of $f_{\mathrm{D}}$. In contrast, in Eq. 6, $\mathcal{R}$ is
minimized for selecting $\theta_{\mathrm{T}}$. Such a clear distinction of the
semantics of $\mathcal{R}$ in each scene can make the tuning of $\lambda$
easier. Moreover, the decoupled optimization allows us to use different
optimizers and datasets for estimating $\theta_{\mathrm{D}}$ and
$\theta_{\mathrm{T}}$.

When we should estimate different $\theta_{\mathrm{T}}$ for each query $x$, we
have two options: One is to solve Eq. 6 with $S^{\prime}_{X}$ being the
singleton, $S^{\prime}_{X}=\\{x\\}$, which is inefficient when there are a
large number of queries. Another option is to solve

$\theta_{h}^{*}=\arg\min_{\theta_{h}}\sum_{x}\mathcal{R}_{S^{\prime}_{X}}(\theta_{\mathrm{T}}=h(x),\theta_{\mathrm{D}}^{*}),$
(7)

where $h\colon\mathcal{X}\to\Theta_{\mathrm{T}}$ is a trainable model
parameterized with $\theta_{h}$ and is used for inferring
$\theta_{\mathrm{T}}$ given $x$ in an amortized manner (i.e., so-called an
encoder).

##### Distribution Estimation of $\theta_{\mathrm{T}}$ (optional)

We can also consider a distribution estimation of $\theta_{\mathrm{T}}$.
Recall that in deep grey-box models, $\mathcal{L}$ (i.e., the supervision from
labeled data) is no longer informative to decide the value of
$\theta_{\mathrm{T}}$ due to the expressive power of $f_{\mathrm{D}}$, and
only $\mathcal{R}$ dictates $\theta_{\mathrm{T}}$. Hence, it is reasonable to
define a distribution of $\theta_{\mathrm{T}}$ only with $\mathcal{R}$ as

$p(\theta_{\mathrm{T}}\mid\text{data})\propto
p(\theta_{\mathrm{T}})\exp\\{-\beta\mathcal{R}_{S^{\prime}_{X}}(\theta_{\mathrm{T}},\theta_{\mathrm{D}}^{*})\\},$
(8)

for some $\beta>0$. If we compute the full landscape of $\mathcal{R}$ for a
low-dimensional $\theta_{\mathrm{T}}$, it is equivalent to having this
distribution explicitly. Even though $p(\theta_{\mathrm{T}}\mid\text{data})$
cannot be normalized when $\theta_{\mathrm{T}}$ is high dimensional, we can
efficiently draw samples from the distribution using MCMC.

###  3.3 Discussion

The increased model complexity of $f_{\mathrm{D}}$ does not come with an
additional need for real training data because we can draw as many random
samples of $\theta_{\mathrm{T}}$ as the computational resources allow in the
estimation process. As a result, while the learning problem becomes (hopefully
slightly) more complicated than the original ones, it would not make the
problem significantly more challenging. We will empirically confirm it through
numerical experiments.

Though the full applicability of the proposed framework is limited to
moderate-dimensional $\theta_{\mathrm{T}}$, in practice, it would not
significantly limit applications because theory-based models often have (or
should have) a small number of parameters. Models with high-dimensional
$\theta_{\mathrm{T}}$ such as unknown fields are a challenging open problem.

As repeatedly argued, the main characteristic of the proposed formulation is
that it makes it easy to analyze $\mathcal{R}$s empirically. Here, care should
be taken not to reuse the same data both for such an analysis and the
estimation with $\mathcal{R}$ selected via the former analysis. We actually
commit “data reuse” in parts of the experiments in Section 5, which is
admittable only because our purpose is not to analyze the data but to compare
the different estimation schemes.

###  3.4 A Numerical Example

We show a scenario of how the framework could be used. We generate data from
$y=\sin(x)+\cos(x)+e$, where $x\in\mathcal{X}=[-\pi,\pi]$, and
$e\sim\mathcal{N}(0,0.1^{2})$. We use a deep grey-box model
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=f_{\mathrm{T}}(x)+f_{\mathrm{D}}(f_{\mathrm{T}},x)$
with $f_{\mathrm{T}}(x)=a\sin(x+c)$ and
$f_{\mathrm{D}}(f_{\mathrm{T}},x)=\operatorname{MLP}(x,\theta_{\mathrm{T}},f_{\mathrm{T}}(x))$,
where $\theta_{\mathrm{T}}=[a,c]$ is the $f_{\mathrm{T}}$’s parameters.
$\operatorname{MLP}(\cdot)$ is a feed-forward neural net. We feed
$\operatorname{MLP}$ with
$[x;\theta_{\mathrm{T}};f_{\mathrm{T}}(x)]\in\mathbb{R}^{4}$, and it returns
values in $\mathbb{R}$. Note that our $f_{\mathrm{D}}$ is aware of
$\theta_{\mathrm{T}}$ and can work adaptively to different values of
$\theta_{\mathrm{T}}$, as advised in Step 1 (Section 3.2.1).

The next step, Step 2 (Section 3.2.2), is to optimize the expected objective
in Eq. 5. We run it with the squared error,
$\mathcal{L}=\sum(y-\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x))^{2}$, and
$\lambda=0$ (i.e., we do not specify $\mathcal{R}$). We let
$p(\theta_{\mathrm{T}})$ be the uniform distribution over
$[0,2]\times[-\pi,\pi]$. As a result, we obtain $\theta_{\mathrm{D}}^{*}$.

As Step 3 (Section 3.2.3), we overview the landscapes of some candidates of
$\mathcal{R}$. In Fig. 2, we visualize the values of the following
$\mathcal{R}$s over a grid of $\theta_{\mathrm{T}}=[a,c]$:

  * •

$\mathcal{R}_{\text{normD}}\\!=\\!\sum_{x}f_{\mathrm{D}}(x)^{2}$ (i.e.,
$f_{\mathrm{D}}$ should work minimally);

  * •

$\mathcal{R}_{\text{corr}}=\left|\sum_{x}f_{\mathrm{T}}(x)\cdot
f_{\mathrm{D}}(x)\right|$ (i.e., $f_{\mathrm{T}}$ and $f_{\mathrm{D}}$ should
work uncorrelatedly); and

  * •

$\mathcal{R}_{\text{normdif}}=\left|\sum_{x}f_{\mathrm{T}}(x)^{2}-\sum_{x}f_{\mathrm{D}}(x)^{2}\right|$
(i.e., $f_{\mathrm{T}}$ and $f_{\mathrm{D}}$ should work to the same extent).

The leftmost heatmap shows $\mathcal{R}_{\text{normD}}$; we see that it
reaches the minimum around $[a,c]=[\sqrt{2},\pi/4]$, which is natural because
$\sin(x)+\cos(x)=\sqrt{2}\sin(x+\pi/4)$. It serves as a “truth” for those who
believe that $f_{\mathrm{D}}$ should work minimally. However, the case is not
over if one has prior knowledge dictating that $f_{\mathrm{T}}$ and
$f_{\mathrm{D}}$ should be uncorrelated. Such a user would like to use
$\mathcal{R}_{\text{corr}}$, which is visualized in the second heatmap from
the left. We see that multiple $\theta_{T}$s are practically equivalent under
this regularizer. If we want to further narrow down the possible choices of
$\theta_{\mathrm{T}}$, we need to consider more criteria. Let us see
$\mathcal{R}_{\text{normdif}}$, which is shown in the center heatmap. We find
that it is not enough, either, to select a (few) $\theta_{\mathrm{T}}$
value(s). Summing $\mathcal{R}_{\text{corr}}$ and
$\mathcal{R}_{\text{normdif}}$ results in the landscape shown in the fourth
heatmap from the left, where we see two minima. In practice, this would be the
most reasonable regularizer that could be designed and explained through these
empirical analyses of $\mathcal{R}$s.

$\mathcal{R}_{\text{normD}}$ |  $\mathcal{R}_{\text{corr}}$ |  $\mathcal{R}_{\text{normdif}}$ |  $\mathcal{R}_{\text{corr}}+\mathcal{R}_{\text{normdif}}$ |  $\mathcal{R}_{\text{corr}}\\!+\\!\mathcal{R}_{\text{normdif}}\\!+\\!c^{2}$
---|---|---|---|---
|  |  |  |
Figure 2: Landscapes of some regularizers and their combinations, for
$f_{\mathrm{D}}$ trained by Eq. 5 with $\lambda=0$ on the toy dataset in
Section 3.4. The horizontal and vertical axes correspond to $a$ and $c$ of
$f_{\mathrm{T}}$, respectively.

We provide one more visualization in the rightmost panel of Fig. 2, just for
completeness. It shows the values of
$\mathcal{R}_{\text{corr}}+\mathcal{R}_{\text{normdif}}+c^{2}$, which has a
single extremum. However, the $c^{2}$ term can hardly be specified in practice
without the knowledge of the data-generating process.

A message of the above story is that we should analyze the landscape of
$\mathcal{R}$, instead of uncritically optimizing it. Suppose that we follow
the prior knowledge dictating that $f_{\mathrm{T}}$ and $f_{\mathrm{D}}$
should be uncorrelated. Applying a standard learning scheme (i.e., optimize
$\mathcal{L}+\lambda\mathcal{R}_{\text{corr}}$ without the adaptivity of
$f_{\mathrm{D}}$) would result in $\theta_{\mathrm{T}}$ being practically a
random choice from one of the many local minima of $\mathcal{R}_{\text{corr}}$
that we see in Fig. 2 (2nd from the left). This is hardly meaningful since
these parameters are supposed to have an interpretation within the context of
the relevant domain knowledge.

Figure 3: Main loss $\mathcal{L}$ and regularizer $\mathcal{R}$ on test data
for the models learned with different values of regularization hyperparameter,
$\lambda$, on the toy dataset in Section 3.4. Error bars show the standard
errors by 20 random trials.

Let us take one of the further options of Step 3, the point estimation.
Suppose using
$\mathcal{R}=\mathcal{R}_{\text{corr}}+\mathcal{R}_{\text{normdif}}+c^{2}$. We
compare our framework, Eqs. 5 and 6, with inductive learning, Eq. 3, and
transductive learning, Eq. 4. We vary the value of $\lambda$ from $0.001$ to
$0.1$. Figure 3 reports the squared error, $\mathcal{L}$, and the regularizer,
$\mathcal{R}$, computed on a test set for each configuration. Our framework
achieves small values of $\mathcal{R}$ practically for any value of $\lambda$;
this is not the case for inductive and transductive learning. The slightly
large $\mathcal{L}$ for the proposed framework with $\lambda=0.1$ is probably
because the training of the adaptive $f_{\mathrm{D}}$ was inhibited by the
large value of $\lambda$. This result supports the use of a small value of
$\lambda$ in the optimization of Eq. 5. We also confirmed that our framework
worked well with $\lambda=0$.

##  4 RELATED WORK

##### Deep Grey-box Modeling

Combination of deep neural nets and theory-driven models has been studied in
various contexts [23, 18, 21, 13, 3, 4, 6, 22, 9, 19, 20, 10, e.g., recent
studies include]; we cannot enumerate more from a large number of studies in
older times and in application domains, so readers are recommended consulting
the references in the aforementioned papers.

Among most related studies, [23] formulate the learning of deep grey-box
models as a constrained optimization problem where they minimize $\mathcal{R}$
under some constraint on $\mathcal{L}$ (e.g., $\mathcal{L}=0$). [18] study a
family of deep grey-box models in the context of variational autoencoders. In
both of these studies, the theory- and data-driven models are learned once and
together. [21] proposes a method to deal with distribution shift utilizing
deep grey-box models. Their method is applicable to our framework, too.

##### Other Relevant Contexts

When combining the models, a theory-driven model may constitute the final
layer (in contrast, we are interested in the opposite, i.e., the data-driven
model at the final layer). [2] developed an epidemiological model whose
parameters are predicted by neural nets. Physics-informed neural nets [12] and
their variants belong also to this category because their loss, the residual
of a differential equation, can be regarded as the final layer. [15] propose
an optimization method for neural nets when the loss is computed based on some
physics-based models.

##### Learning Schemes

Our method can minimize a part of the objective at prediction time. This is
also the case with model-agnostic meta-learning (MAML) [7]. Differently from
our setting, MAML is meant for settings where one has many tasks and wants to
adapt to a new task at prediction time. Meta-tailoring [1] can be thought of
as a variant of MAML with each sample being one task. Nonetheless, these
methods do not necessarily allow us to efficiently explore objectives’
landscape.

[16] proposed a method to learn deep neural nets with supervision from both
data and rules. They suggest “marginalizing out” a parameter during training
so that we can choose it freely in the inference time. Despite such
similarities, there are several differences from our framework, e.g., in the
target of the marginalization and in how theory is incorporated, which prevent
direct comparison.

##  5 EXPERIMENTS

We demonstrate the capability of the proposed framework to explore the
landscape of regularizers. We also compare the proposed framework with
existing learning schemes when we take the option of point estimation. The
source codes of the experiments are available at
https://github.com/n-takeishi/deepgreybox.

###  5.1 Datasets, Models, and Optimization

In all cases, $\mathcal{L}$ is the mean squared error.

##### Controlled Pendulum

We use time-series data of a frictionless compound pendulum controlled by an
unknown regulator. The dataset comprises pairs
$(x=s_{t},y=(s_{t+1},\dots,s_{t+10}))$, where $s_{t}=[s_{t,1},s_{t,2}]$ is the
state (angle $s_{t,1}$ and angular velocity $s_{t,2}$) of the pendulum. Hence,
the task is to predict the next 10 steps future given a state.

We use the following model:
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=\operatorname{ODESolve}[$
$\dot{s}_{t}=f_{\mathrm{T}}(s_{t})+f_{\mathrm{D}}(s_{t},\theta_{T},f_{\mathrm{T}}(s_{t}))\mid
s_{0}=x]$, where
$f_{\mathrm{T}}(s_{t})=[s_{t,2},\frac{3\theta_{\mathrm{T}}}{2}\sin(s_{t,1})]$,
and $f_{\mathrm{D}}(s_{t},\theta_{\mathrm{T}},f_{\mathrm{T}}(s_{t}))$ is a
network with fully-connected layers. $f_{\mathrm{D}}$ should mimic the
behavior of the unknown controller. It is not obvious what $\mathcal{R}$ we
should use because we do not know the controller’s nature. In Eq. 5, we set
$\mathcal{R}=\mathcal{R}_{\text{normD}}+\mathcal{R}_{\text{corr}}$, where
$\mathcal{R}_{\text{normD}}=\sum_{x}\|f_{\mathrm{D}}(x)\|_{2}^{2}$ and
$\mathcal{R}_{\text{corr}}=|\sum_{x}f_{\mathrm{T}}(x)\cdot f_{\mathrm{D}}(x)|$
are candidate $\mathcal{R}$s, with a small coefficient $\lambda=0.001$.

##### Reaction–Diffusion System

We generated data from the two-dimensional two-component reaction–diffusion
system of the FitzHugh–Nagumo type: $\partial u/\partial t=0.0015\Delta
u+u-u^{3}-v+0.005$ and $\partial v/\partial t=0.005\Delta v+u-v$, where
$u,v\in\mathbb{R}^{32\times 32}$ are the concentration of two substances. The
dataset comprises pairs $(x=s_{0},y=(s_{1},s_{2},\dots,s_{15}))$ where
$s_{t}=[u_{t},v_{t}]$. We note that similar configurations of the same system
have been used in previous related studies [23, 21].

We use a grey-box neural ODE as in the previous example, with
$f_{\mathrm{T}}(s_{t})=[a\Delta u_{t},b\Delta v_{t}]$ whose parameter is
$\theta_{\mathrm{T}}=[a,b]$ and
$f_{\mathrm{D}}(s_{t},\theta_{\mathrm{T}})=\operatorname{ConvNet}_{\theta_{\mathrm{T}}}(s_{t})$,
a convolutional net whose filter is partially parameterized by
$\theta_{\mathrm{T}}$. In the previous studies [23, 21], it has been reported
that $\mathcal{R}_{\text{normD}}=\sum_{x}\|f_{\mathrm{D}}(x)\|_{2}^{2}$ works
to some extent, yet without detailed analysis of its behavior. Hence, we use
$\mathcal{R}=\mathcal{R}_{\text{normD}}+\mathcal{R}_{\text{corr}}$ with a
small value of $\lambda$ to examine the both regularizers as candidates.

##### Predator–Prey System

We use real data of a planktonic predator–prey system [5]. From the original
data, we extracted the measurements of the population densities of the prey
(unicellular algae) and the predator (rotifer). We split the original time-
series into subsequences of length $11$ [days] and created a dataset
comprising $x=y=(s_{t},\dots,s_{t+10})$ where
$s_{t}=[s_{t,1},s_{t,2}]\in\mathbb{R}^{2}$ is the state (prey density
$s_{t,1}$ and predator density $s_{t,2}$) at day $t$. Hence, the task is
autoencoding.

The model comprises a decoder and two encoders. The decoder is again a grey-
box neural ODE:
$\mathcal{C}(f_{\mathrm{T}},f_{\mathrm{D}};x)=\operatorname{ODESolve}\left[\dot{s}_{t}=f_{\mathrm{T}}(s_{t})+f_{\mathrm{D}}(z,\\!s_{t},\\!\theta_{T},\\!f_{\mathrm{T}}(s_{t}))\mid
s_{0}=x_{0}\right]$, where $z$ is a latent variable. $f_{\mathrm{T}}$ is the
Lotka–Volterra equations: $f_{\mathrm{T}}(s_{t})=[\alpha s_{t,1}-\beta
s_{t,1}s_{t,2},-\gamma s_{t,2}+\delta s_{t,1}s_{t,2}]$, where
$\theta_{\mathrm{T}}=[\alpha,\beta,\gamma,\delta]$ ($\alpha$: prey’s growth
rate, $\beta$: prey’s decay rate, $\gamma$: predator’s decay rate, and
$\delta$: predator’s growth rate). One of the encoders is a neural net with
fully-connected layers that takes $x$ as input and outputs $z$. We run Eq. 5
with $\mathcal{R}=\sum_{x}\|f_{\mathrm{D}}(x)\|_{2}^{2}$ and $\lambda=0.001$
to explore the $\mathcal{R}$’s landscape. Then, we run Eq. 7 to get another
encoder that infers $\theta_{\mathrm{T}}$ for each $x$.

###  5.2 Explore the Landscape of $\mathcal{R}$

Figure 4: (_Left_ & _center_) Landscapes of the regularization terms,
$\mathcal{R}_{\text{normD}}$ and $\mathcal{R}_{\text{corr}}$ computed on the
test set of the controlled pendulum dataset. (_Right_) Test NRMSEs.

Figure 4 shows, for the controlled pendulum dataset, the values of
$\mathcal{R}_{\text{normD}}$ and $\mathcal{R}_{\text{corr}}$ for different
$\theta_{\mathrm{T}}$s. We can see that the two regularizers have slightly
different peaks. While it is not easy to decide which one is better in some
sense, it is notable that our method can provide such empirical clues on the
differences between the regularizers. The rightmost plot of Fig. 4 shows the
normalized root mean squared errors (NRMSEs) of the prediction on the test
set. The model can predict well regardless of the value of
$\theta_{\mathrm{T}}$.

$\mathcal{R}_{\text{normD}}$ |  $\mathcal{R}_{\text{corr}}$ |  $\mathcal{R}_{\text{normD}}+\mathcal{R}_{\text{corr}}$ |  $\mathcal{R}_{\text{normD}}\cdot\mathcal{R}_{\text{corr}}$ | RMSE
---|---|---|---|---
|  |  |  |
Figure 5: (_Left four_) Landscapes of $\mathcal{R}_{\text{normD}}$,
$\mathcal{R}_{\text{corr}}$, and their combinations, for the
reaction–diffusion system dataset. The horizontal and vertical axes correspond
to $a$ and $b$ of $f_{\mathrm{T}}$, respectively. (_Right_) Test RMSEs.

Figure 5 visualizes the $\mathcal{R}$ candidates and the test RMSE for the
reaction–diffusion system dataset, analogously to what was explained in the
previous paragraph. $\mathcal{R}_{\text{normD}}$ has a peak at some value,
while the extrema of $\mathcal{R}_{\text{corr}}$ are much less clear. We found
that the product of these two terms,
$\mathcal{R}_{\text{normD}}\cdot\mathcal{R}_{\text{corr}}$ , also have clear,
if not unique, peaks. Interestingly, one of the local minima of
$\mathcal{R}_{\text{normD}}\cdot\mathcal{R}_{\text{corr}}$ nicely points the
data-generating value of $\theta_{\mathrm{T}}=[a,b]$, that is,
$[a,b]=[0.0015,0.005]$. On the other hand, $\mathcal{R}_{\text{normD}}$, which
has been used in the previous studies, does not point to this value as the
minimum. This suggests that, again, we should empirically analyze
$\mathcal{R}$ to know how different $\mathcal{R}$ may result in different
$\theta_{\mathrm{T}}$ estimations, instead of optimizing it uncritically. We
are _not_ suggesting
$\mathcal{R}_{\text{normD}}\cdot\mathcal{R}_{\text{corr}}$ is a “better”
regularizer; we never know the truth of $\theta_{\mathrm{T}}$ in practice and
thus cannot evaluate anything from such a viewpoint.

prey |
---|---
predator |
$\alpha$ |
$\beta$ |
$\gamma$ |
$\delta$ |
$\|f_{\mathrm{D}}(x)\|_{2}^{2}$ |
NRMSE [%] |

$\beta$ |  |  — |  —
---|---|---|---
$\gamma$ |  |  |  —
$\delta$ |  |  |
|  $\alpha$ |  $\beta$ |  $\gamma$

Figure 6: (_Left_) A test sequence of the predator–prey system dataset,
inferred $\alpha,\beta,\gamma,\delta$, $\|f_{\mathrm{D}}(x)\|_{2}^{2}$, and
prediction NRMSE. (_Right_) Partial landscape of
$\mathcal{R}=\sum_{x}\|f_{\mathrm{D}}(x)\|_{2}^{2}$ at $t=65$.

Figure 6 (left) shows an outcome of the estimated model on a test sequence of
the predator–prey system dataset. We took a sliding window of length $11$ from
the long test sequence, applied the model to the subsequences from the window,
and calculated the averages of the outputs for the overlapping steps. The
value of $\|f_{\mathrm{D}}(x)\|_{2}^{2}$ at each time is informative to assess
how well the $f_{\mathrm{T}}$ could explain the data. We can observe that the
latter half of the sequence is relatively well explained by $f_{\mathrm{T}}$.
We can also analyze the landscape of the regularizer,
$\mathcal{R}=\sum_{x}\|f_{\mathrm{D}}(x)\|_{2}^{2}$. Figure 6 (right) shows
the slices of the values of the $\mathcal{R}$. Each slice is computed by
fixing the values of the remaining elements of $\theta_{\mathrm{T}}$ at the
encoder’s outputs. An observation is that while $\mathcal{R}$ seems to have
some extrema along the directions of $\beta$, $\gamma$, $\delta$, it looks
quite flat along the direction of $\alpha$. It implies the possible difficulty
of inferring $\alpha$, prey’s growth rate without interaction, from the data,
and thus we should not immediately interpret at least the results of
$\alpha$’s inference.

###  5.3 Point Estimation of $\theta_{\mathrm{T}}$

Figure 7: Loss $\mathcal{L}$ and regularizer $\mathcal{R}$ on test set of the
controlled pendulum dataset. Error bars show the standard errors by 20 random
trials.

Figure 8: Loss $\mathcal{L}$ and regularizer $\mathcal{R}$ on test set of the
reaction–diffusion system dataset. Error bars show the standard errors by 20
random trials.

We performed one of the further options in the proposed framework. For the
controlled pendulum dataset, we set the regularizer as
$\mathcal{R}=\mathcal{R}_{\text{corr}}$ and performed the point estimation by
Eq. 3 (inductive), Eq. 4 (transductive), and Eq. 6 (proposed). Figure 7 shows
the values of $\mathcal{L}$ and $\mathcal{R}$ on the test set. $\mathcal{L}$
is not well minimized with a large $\lambda$ for all the methods, which is
natural. A notable difference is that by the proposed method, $\mathcal{R}$ is
well minimized even with small $\lambda$, though it is not surprising because
the optimization is decoupled. In contrast, the inductive and transductive
learning schemes result in larger $\mathcal{R}$ values with small $\lambda$s,
due to which striking a good value of $\lambda$ may be difficult for these
learning schemes. Figure 8 reports basically the same thing but for the
reaction–diffusion system dataset.

##  6 CONCLUSION

Deep grey-box models are combinations of deep neural networks and theory-
driven models. We argued that, toward trustworthy estimation of the theory-
driven part of the model, we should empirically analyze the regularizers we
use. We suggested a slight modification of the network architecture and the
optimization objective used in deep grey-box models such that we can
empirically analyze regularizers after training. The proposed formulation is
also useful as it decouples the training of the data-driven and theory-driven
models. The main limitation is the applicability to theory-driven models with
high-dimensional parameters. They pose issues in two aspects: computation of
the expectation with regard to the theory parameters and visualization of a
regularizer’s landscape.

#### Acknowledgements

This work was supported by the Innosuisse project _Industrial artificial
intelligence for intelligent machines and manufacturing digitalization_
(39453.1 IP-ICT) and the Swiss National Science Foundation Sinergia project
_Modeling pathological gait resulting from motor impairments_ (CRSII5_177179).

## References

  * [1] Ferran Alet, Maria Bauza, Kenji Kawaguchi, Nurullah Giray Kuru, Tomas Lozano-Perez and Leslie Pack Kaelbling  “Tailoring: Encoding Inductive Biases by Optimizing Unsupervised Objectives at Prediction Time”  In _Advances in Neural Information Processing Systems 34_ , 2021, pp. 29206–29217
  * [2] Sercan Ö. Arık, Chun-Liang Li, Jinsung Yoon, Rajarishi Sinha, Arkady Epshteyn, Long T. Le, Vikas Menon, Shashank Singh, Leyou Zhang, Nate Yoder, Martin Nikoltchev, Yash Sonthalia, Hootan Nakhost, Elli Kanal and Tomas Pfister  “Interpretable Sequence Learning for COVID-19 Forecasting”  In _Advances in Neural Information Processing Systems 33_ , 2020, pp. 18807–18818
  * [3] Abigail R. Azari, Jeffrey W. Lockhart, Michael W. Liemohn and Xianzhe Jia  “Incorporating Physical Knowledge into Machine Learning for Planetary Space Physics”  In _Frontiers in Astronomy and Space Sciences_ 7, 2020, pp. 36
  * [4] Timur Bikmukhametov and Johannes Jäschke  “Combining Machine Learning and Process Engineering Physics towards Enhanced Accuracy and Explainability of Data-Driven Models”  In _Computers & Chemical Engineering_ 138, 2020, pp. 106834
  * [5] Bernd Blasius, Lars Rudolf, Guntram Weithoff, Ursula Gaedke and Gregor F. Fussmann  “Long-Term Cyclic Persistence in an Experimental Predator–Prey System”  In _Nature_ 577, 2020, pp. 226–230
  * [6] Weinan E, Jiequn Han and Linfeng Zhang  “Integrating Machine Learning with Physics-Based Modeling”, 2020  arXiv:2006.02619
  * [7] Chelsea Finn, Pieter Abbeel and Sergey Levine  “Model-Agnostic Meta-Learning for Fast Adaptation of Deep Networks”  In _Proceedings of the 34th International Conference on Machine Learning_ , 2017, pp. 1126–1135
  * [8] Alex Gammerman, Volodya Vovk and Vladimir Vapnik  “Learning by Transduction”  In _Proceedings of the 14th Conference on Uncertainty in Artificial Intelligence_ , 1998, pp. 148–155
  * [9] George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang and Liu Yang  “Physics-Informed Machine Learning”  In _Nature Reviews Physics_ 3, 2021, pp. 422–440
  * [10] Johan Kon, Dennis Bruijnen, Jeroen Wijdeven, Marcel Heertjes and Tom Oomen  “Unifying Model-Based and Neural Network Feedforward: Physics-guided Neural Networks with Linear Autoregressive Dynamics”, 2022  arXiv:2209.12489
  * [11] Zhaozhi Qian, William R. Zame, Lucas M. Fleuren, Paul Elbers and Mihaela Schaar  “Integrating Expert ODEs into Neural ODEs: Pharmacology and Disease Progression”  In _Advances in Neural Information Processing Systems 34_ , 2021, pp. 11364–11383
  * [12] Maziar Raissi, Paris Perdikaris and George E. Karniadakis  “Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations”  In _Journal of Computational Physics_ 378, 2019, pp. 686–707
  * [13] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais and Prabhat  “Deep Learning and Process Understanding for Data-Driven Earth System Science”  In _Nature_ 566.7743, 2019, pp. 195–204
  * [14] Riku Sasaki, Naoya Takeishi, Takehisa Yairi and Koichi Hori  “Neural Gray-Box Identification of Nonlinear Partial Differential Equations”, Lecture Notes in Computer Science 11671, 2019, pp. 309–321
  * [15] Patrick Schnell, Philipp Holl and Nils Thuerey  “Half-Inverse Gradients for Physical Deep Learning”  In _Proceedings of the 10th International Conference on Learning Representations_ , 2022
  * [16] Sungyong Seo, Sercan O. Arik, Jinsung Yoon, Xiang Zhang, Kihyuk Sohn and Tomas Pfister  “Controlling Neural Networks with Rule Representations”  In _Advances in Neural Information Processing Systems 34_ , 2021, pp. 11196–11207
  * [17] B. Sohlberg and E.W. Jacobsen  “Grey Box Modelling – Branches and Experiences”  In _IFAC Proceedings Volumes_ 41.2, 2008, pp. 11415–11420
  * [18] Naoya Takeishi and Alexandros Kalousis  “Physics-Integrated Variational Autoencoders for Robust and Interpretable Generative Modeling”  In _Advances in Neural Information Processing Systems 34_ , 2021, pp. 14809–14821
  * [19] Laura Rueden, Sebastian Mayer, Katharina Beckh, Bogdan Georgiev, Sven Giesselbach, Raoul Heese, Birgit Kirsch, Julius Pfrommer, Annika Pick, Rajkumar Ramamurthy, Michal Walczak, Jochen Garcke, Christian Bauckhage and Jannis Schuecker  “Informed Machine Learning – A Taxonomy and Survey of Integrating Knowledge into Learning Systems”  In _IEEE Transactions on Knowledge and Data Engineering_ , 2021
  * [20] Rui Wang  “Physics-Guided Deep Learning for Dynamical Systems: A Survey”, 2021  arXiv:2107.01272
  * [21] Antoine Wehenkel, Jens Behrmann, Hsiang Hsu, Guillermo Sapiro, Gilles Louppe and Jörn-Henrik Jacobsen  “Robust Hybrid Learning with Expert Augmentation”, 2022  arXiv:2202.03881
  * [22] Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach and Vipin Kumar  “Integrating Physics-Based Modeling with Machine Learning: A Survey”, 2020  arXiv:2003.04919
  * [23] Yuan Yin, Vincent Le Guen, Jérémie Dona, Emmanuel Bézenac, Ibrahim Ayed, Nicolas Thome and Patrick Gallinari  “Augmenting Physical Models with Deep Networks for Complex Dynamics Forecasting”  In _Proceedings of the 9th International Conference on Learning Representations_ , 2021

##  Appendix A ADDITIONAL INFORMATION ON EXPERIMENTS

###  A.1 Common

##### Optimization

In all the experiments including the one in Section 3.4, the objective of Eq.
5 and its gradients were estimated using only one sample of
$\theta_{\mathrm{T}}$ per an instance of $x$, that is, $m=1$. We also tried
larger $m$ up to $m=100$ per a single $x$, but the final improvement of
performance was marginal. We used the AdamW optimizer for the optimization of
inductive learning, Eq. 3, transductive learning, Eq. 4, and the training-time
optimization of the proposed method, Eq. 5. The parameters of the optimizer
were set to the default values of the library unless stated otherwise.

##### Others

  * •

All implementations were done with PyTorch 1.11.0.

  * •

The NRMSEs were computed by dividing RMSEs by $y_{\text{max}}-y_{\text{min}}$.

###  A.2 Toy Dataset (in Section 3.4)

##### Data

In the example in Section 3.4, we generated the data by sampling $x$ from the
uniform distribution on $\mathcal{X}=[-\pi,\pi]$. We created training,
validation, and test datasets, each of which was with $40$ samples of $(x,y)$.

##### Model

In the model, $f_{\mathrm{D}}$ is a neural net with fully-connected layers
(with two hidden layers of size $16$) and the leaky ReLU activation function.

##### Optimization

We optimized the model parameters ($\theta_{\mathrm{T}}=[a,c]$ and
$\theta_{\mathrm{D}}$ in the inductive and transductive learning schemes, Eqs.
3 and 4; only $\theta_{\mathrm{D}}$ in the proposed method, Eq. 5) with
learning rate varied from $0.01$ to $0.0001$ exponentially. We set the mini-
batch size to $10$ and ran the optimization for $2000$ epochs. In the
inductive and transductive learning schemes, the value of
$\theta_{\mathrm{T}}=[a,c]$ was constrained to be $0\leq a\leq 2$ and
$-\pi\leq c\leq\pi$. In the training-time optimization of the proposed method,
Eq. 5, we sampled $\theta_{\mathrm{T}}=[a,c]$ from the uniform distribution on
$0\leq a\leq 2$ and $-\pi\leq c\leq\pi$. The prediction-time optimization of
the proposed method, Eq. 6, was performed with the Adam optimizer (without
weight decay) with learning rate varied from $0.01$ to $0.0001$. We ran it for
$2000$ epochs with the full batch.

##### Test prediction errors

Figure A.1: Test RMSEs in the result reported in Fig. 2.

In Fig. 2 of the main text, we had no room to show the value of $\mathcal{L}$
evaluated on the test set for the same range of $\theta_{\mathrm{T}}$ as that
in the plots of Fig. 2. We report it here in Fig. A.1. Although some extreme
values of $\theta_{\mathrm{T}}=[a,c]$ tend to have slightly large
$\mathcal{L}$, we can see that the model predicts similarly well with
different $\theta_{\mathrm{T}}$s.

###  A.3 Controlled Pendulum

##### Data

We used the dataset available222https://github.com/Stable-Baselines-
Team/stable-
baselines/blob/f877c85b6e45084da1d8ccf73e7e730dc2001c3f/stable_baselines/gail/dataset/expert_pendulum.npz
on the GitHub repository of the Stable Baselines
package333https://github.com/Stable-Baselines-Team/stable-baselines/. It is a
collection of sequences measuring the state of a frictionless compound
pendulum, with a rigid uniform rod of length $1$. The gravitational
acceleration is $g=10$ in the data-generation environment. In each sequence,
the pendulum starts at a random position and is controlled to stay vertically
upright by some controller. The detail of the controller does not matter, so
we just leave it unknown. The original dataset comprises sequences of length
$200$. We disposed of the last $100$ steps of each sequence because most of
the time, the pendulum is already almost stabilized at the goal state by that
time. We ran a sliding window of size $11$ on the remaining sequences and
created pairs $(x,y)$ with $x$ being the first snapshot in the window and $y$
being the last $10$ snapshots in the window. We split the original dataset
into training, validation, and test sets of size $3600$, $2700$, and $2700$,
respectively.

##### Model

In the model, the neural net,
$\operatorname{MLP}(t,s_{t},\theta_{\mathrm{T}},f_{\mathrm{T}}(s_{t}))$, is a
feed-forward network with fully-connected layers with three hidden layers of
size $128$. It takes the concatenation of $s_{t}$, $\theta_{\mathrm{T}}$, and
$f_{\mathrm{T}}(s_{t})$ as the input (without using $t$ directly). We used the
torchdiffeq library444https://github.com/rtqichen/torchdiffeq/ for the
numerical integration of the ODE with the 4th-order Runge–Kutta method.

##### Optimization

In training, we varied the learning rate from $0.001$ to $0.00001$
exponentially, set the mini-batch size to $50$, and ran the optimization for
$500$ epochs. In the inductive and transductive learning schemes, the value of
$\theta_{\mathrm{T}}$ was constrained to be in $[8,12]$. In the training-time
optimization of the proposed method, Eq. 5, we sampled $\theta_{\mathrm{T}}$
from the uniform distribution on $[8,12]$. The prediction-time optimization of
the proposed method, Eq. 6, was performed simply by the grid search because
$\theta_{\mathrm{T}}$ is one-dimensional.

###  A.4 Reaction–Diffusion System

##### Data

We generated the data from the two-dimensional two-component
reaction–diffusion system of the FitzHugh–Nagumo type:

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=0.0015\Delta
u+u-u^{3}-v+0.005,$ $\displaystyle\frac{\partial v}{\partial t}$
$\displaystyle=0.005\Delta v+u-v.$

The two variables of the system, $u$ and $v$, are defined over the two-
dimensional space spanning $[-1,1]\times[-1,1]$. We discretized the space with
the $32\times 32$ even grid. In generating the data and in the trained model,
the Laplacian operator, $\Delta$, was computed using the five-point stencil.
We generated $1000$ sequences of length $16$ from $t=0$ to $t=1.5$ using the
4th-order Runge–Kutta method with step size $0.001$. We prepared training,
validation, and test sets of size $400$, $300$, and $300$, respectively.

##### Model

In the model, the neural net,
$\operatorname{ConvNet}_{\theta_{\mathrm{T}}}(s_{t})$, returns the sum of the
outputs of two subnetworks. One of them is a convolutional network with two
hidden layers having $16$ channels and with the leaky ReLU activation
function. Another subnetwork is a two-layer convolutional network whose
filters’ weights are the discrete Laplacian operator (with the five-point
stencil) multiplied by scalars for each filter. Those scalars are the output
of a feed-forward network that takes $\theta_{\mathrm{T}}$ as input. It has
two hidden layers of size $128$ with the leaky ReLU activation function
applied to the intermediate layers and the hyperbolic tangent applied to the
final layer. When the output of this network is below
$-\theta_{\mathrm{T}}/2$, it is clipped to be $-\theta_{\mathrm{T}}/2$ so that
the second subnetwork does not completely cancel the output of
$f_{\mathrm{T}}$.

##### Optimization

In training, we fixed the learning rate at $0.001$, set the mini-batch size to
$20$, and ran the optimization for $1000$ epochs. In the inductive and
transductive learning schemes, the value of $\theta_{\mathrm{T}}=[a,b]$ was
constrained to be in $[0.001,0.002]\times[0.001,0.01]$. In the training-time
optimization of the proposed method, Eq. 5, we sampled $\theta_{\mathrm{T}}$
from the uniform distribution on $[0.001,0.002]\times[0.001,0.01]$. The
prediction-time optimization of the proposed method, Eq. 6, was performed with
the Adam optimizer (without weight decay) with learning rate $0.001$. We ran
it for $100$ epochs with the full batch.

###  A.5 Predator–Prey System

##### Data

We used the data available
online555https://doi.org/10.6084/m9.figshare.10045976.v1, which contain
measurements of the population densities of the prey (unicellular algae) and
the predator, (rotifer). In order to adjust the scale of the data values, we
roughly normalized the data by multiplying $0.5$ and $0.02$ by the population
densities of the two species, respectively. We excluded the part of data
containing a value larger than $7$ (after the normalization) from the training
and validation sets; the test set might include such parts, but we did not
check it.

##### Model

The model has three neural networks in total: $f_{\mathrm{D}}$, the encoder of
$z$, and the encoder of $\theta_{\mathrm{T}}$. These networks have the same
architecture: three hidden layers of size $128$ and the leaky ReLU activation
function.

##### Optimization

In training, we fixed the learning rate at $0.001$, set the mini-batch size to
$100$, and ran the optimization for $500$ epochs. In the inductive and
transductive learning schemes, the value of
$\theta_{\mathrm{T}}=[\alpha,\beta,\gamma,\delta]$ was constrained to be in
$[0,1.5]\times[0,3]\times[0,1.5]\times[0,3]$. In the training-time
optimization of the proposed method, Eq. 5, we sampled $\theta_{\mathrm{T}}$
from the uniform distribution on $[0,1.5]\times[0,3]\times[0,1.5]\times[0,3]$.
The prediction-time optimization of the proposed method, Eq. 6, was performed
with the Adam optimizer (without weight decay) with learning rate $0.001$. We
ran it for $50$ epochs with the mini-batch of size $100$.


* Schneider et al. (2020) Schneider, F.R.N.; Ohlmann, S.T.; Podsiadlowski, P.; Röpke, F.K.; Balbus, S.A.; Pakmor, R.  Long-term evolution of a magnetic massive merger product.  Mon. Not. R. Astron. Soc. 2020, 495, 2796–2812.  https://doi.org/10.1093/mnras/staa1326.
  * Caiazzo et al. (2021) Caiazzo, I.; Burdge, K.B.; Fuller, J.; Heyl, J.; Kulkarni, S.R.; Prince, T.A.; Richer, H.B.; Schwab, J.; Andreoni, I.; Bellm, E.C.; et al.  A highly magnetized and rapidly rotating white dwarf as small as the Moon.  Nature 2021, 595, 39–42.  https://doi.org/10.1038/s41586-021-03615-y.
  * Keszthelyi et al. (2021) Keszthelyi, Z.; Meynet, G.; Martins, F.; de Koter, A.; David-Uraz, A.  The effects of surface fossil magnetic fields on massive star evolution—III. The case of $\tau$ Sco.  Mon. Not. R. Astron. Soc. 2021, 504, 2474–2492.  https://doi.org/10.1093/mnras/stab893.
  * Shultz et al. (2019) Shultz, M.E.; Johnston, C.; Labadie-Bartz, J.; Petit, V.; David-Uraz, A.; Kochukhov, O.; Wade, G.A.; Pepper, J.; Stassun, K.G.; Rodriguez, J.E.; et al.  MOBSTER—III. HD 62658: A magnetic Bp star in an eclipsing binary with a non-magnetic ‘identical twin’.  Mon. Not. R. Astron. Soc. 2019, 490, 4154–4165.  https://doi.org/10.1093/mnras/stz2846.
  * Shultz et al. (2015) Shultz, M. et al. [BinaMIcS Collaboration]  Detection of magnetic fields in both B-type components of the $\epsilon$ Lupi system: A new constraint on the origin of fossil fields?  Mon. Not. R. Astron. Soc. 2015, 454, L1–L5.  https://doi.org/10.1093/mnrasl/slv096.
  * Moe and Di Stefano (2017) Moe, M.; Di Stefano, R.  Mind Your Ps and Qs: The Interrelation between Period (P) and Mass-ratio (Q) Distributions of Binary Stars.  Astrophys. J. Suppl. Ser. 2017, 230, 15.  https://doi.org/10.3847/1538-4365/aa6fb6.
  * de Mink et al. (2013) de Mink, S.E.; Langer, N.; Izzard, R.G.; Sana, H.; de Koter, A.  The Rotation Rates of Massive Stars: The Role of Binary Interaction through Tides, Mass Transfer, and Mergers.  Astrophys. J. Lett. 2013, 764, 166.  https://doi.org/10.1088/0004-637X/764/2/166.
  * Mathys (1988) Mathys, G.  A search for magnetic fields in blue stragglers of M67.  Messenger 1988, 53, 39–42.
  * Freytag et al. (2002) Freytag, B.; Steffen, M.; Dorch, B.  Spots on the surface of Betelgeuse—Results from new 3D stellar convection models.  Astron. Nachrichten 2002, 323, 213–219.  https://doi.org/10.1002/1521-3994.
  * Tessore et al. (2017) Tessore, B.; Lèbre, A.; Morin, J.; Mathias, P.; Josselin, E.; Aurière, M.  Measuring surface magnetic fields of red supergiant stars.  Astron. Astrophys. 2017, 603, A129.  https://doi.org/10.1051/0004-6361/201730473.
  * Keszthelyi (2022) Keszthelyi, Z.  Reproduction package for the paper “The effects of surface fossil magnetic fields on massive star evolution: IV. Grids of models at Solar, LMC, and SMC metallicities”.  Zenodo Software Package, 2022; id. 7069766.  https://doi.org/10.5281/zenodo.7069766.


# OPTIMIZATION OF THE ANALYSIS OF DATA OF VITAL EVENTS USING THREADS

Coasaca Callacondo Rubi melania
<EMAIL_ADDRESS>
Chata Iscarra Grylia Yaneth
<EMAIL_ADDRESS>
Fred Torres-Cruz
<EMAIL_ADDRESS>
(28 November 2022)

## Abstra

Objective: Optimize the time of data analysis of Vital Events (births, deaths
and marriages) using Threads.

Methodology: A code was created in python without threads and another with
threads, after He performed 5 tests with a single attribute and with
1,2,3,4,5,6,7,8 and 9 attributes this was done in both codes with the same
amount of data, to know if the python code with threads is optimal when
parsing Vital Facts data.

Results: The python code with Threads turned out to be the most optimal since
it optimized the compilation time of the 5 tests with 1 attribute by
16attributes was obtained as a result that the more attributes you group, the
more effective the use of threads.

Conclusion: Python code with threads is more optimal than code without threads
Therefore, it is concluded that the implementation of threads is recommended
in the analysis of data in similar works.

## Introduction:

Peru generates millions of data records daily, one of those records are the
Vital Facts (births, deaths and marriages), this record is increasing day by
day due to the fact that record events are generated daily. As there is such a
large amount of data, the analysis of said data takes time and as a solution
to this problem, what is intended in this study is to optimize the analysis of
500,000 data from the Vital Events registry, implementing a code in python
where the data analysis will be optimized through parallel computing,
implementing threads. This optimization will be demonstrated by the read time
and compile count of 9 attributes.

## Methodology:

Parallel computing: Computing parallel divides a big problem into several
parts small, separate, sometimes similar can be executed simultaneously by
multiple processors that communicate through shared memory, the results are
finally they execute as if the work had been done in a single general
algorithm.

Threads: Threads are multiple threads running simultaneously in parallel and
performing different tasks in a single program, [7] Memory-sharing threads
will be created to speed up the process of reading and counting the 500000
data with 9 attributes.

Data: The data corresponds to the National Open Data Platform, Governance
Vital Facts option, which can be viewed at the following link:bit.ly/3udD9Wh
This information is provided by the National Institute of Statistics and
Information ethics - INEI, source-RENIEC.

Data description: We will work with 500,000 vital event data that are in CSV
format, containing 9 attributes. The last time this data was updated was in
September 2022. The data contains registrations of births, deaths, and
marriages registered in civil status registry offices nationwide, broken down
by year and month of registration, sex, vital event, type of registration,
department, province and district where the registration was made.

Data Dictionary:

Implementation Logic: A code will be made in python language where parallel
computation will be applied implementing threads [8]

\- Load the data of Vital Facts that are in CSV format, the code is in python
language but it is executed in google colab.
\- The code starts reading and counting the data until it finishes all the
traversal of the requested attribute.
\- Performs a reading and counting of 500,000 data, with 1, 2 or up to 9
attributes that are indicated to the code.
\- The reading and counting of the 500000 data is done in parallel, shortening
the execution time through the implementation of threads.

Implementation: The following libraries, functions and coherence were used to
make the python codes with threads and without threads. -the pandas library is
used for the export and processing of the 500000 data.

-The time library is in charge of controlling the execution time of the data.

-The threading library allows code to execute multiple operations simultaneously.

-The Numpy library is used to create matrices faster and store a large amount of information.[2]

-The data is imported, we continue with the extraction of the data attributes; and select the attributes to work with.

-The timer starts once the attributes have been assigned to it.

-The function count data is in charge of counting all the data that belongs to the requested attribute.

-It begins with the first attribute that was assigned and its respective reading and count.

-Successively, it is passed from attribute to attribute with its respective count according to the requested attribute.

-Prints the amount of data and the time it took to count and read it.

Pseudocode: Pseudocode is one of the most widely used methods for designing
algorithms, and it is independent of the programming language that is going to
be used.[11]

Pseudocode: [3]

Libraries and threadless code:

libraries and code with threads:

The python code is shown where it performs a reading and a count of the data
according to the requested attributes, then only the part of lines of code
where the threads are implemented.

## Results:

In both codes the 500000 were compiled data of vital events first with 1
attribute, 2 attributes,…, until completing the 9 attributes, likewise 5 tests
were carried out with 5 independent attributes (year, month of registration,
sex, event and registration 1) .

Time difference with 9 attributes:

Graphic:

Interpretation: In the Bar Graph it is visualized that the difference that
exists between attribute 1 and attribute 9 is wide, it is also observed that
the bars are in ascending order, this leads us to say that the more attributes
are compiled, the more effective the use of threads.

mean test:

Interpretation:

According to the table, it can be seen that in the data analyzed without
optimization we have an average of 42.55872547 seconds while with optimization
we have an average of 16.56789342 seconds, which indicates that it is better
to use optimization through threads. Therefore, we say that we can obtain
response times and the percentage of improvement according to the number of
attributes analyzed, resulting in a percentage of between 60% to 65%
improvement in each test, as can be seen in the following table.

Interpretation: Applying the means test on the data obtained with an Alpha
equal to 0.05, we reject the null hypothesis, which means that if there is a
significant difference in the response when comparing both in the year and
registration of a total of 500000 data.

Time difference of 5 tests:

Graphic:

Interpretation: The bar graph shows that the difference between test 1 and
test 5 is almost similar, this leads us to say that by individually compiling
each attribute, a 16% improvement can be seen when using threads.

average: To do this we will use the following formula:

Interpretation: The average of the differences of the 5 tests with independent
attributes between the python code with threads and the python code without
threads is 16%.

## Conclution:

A noticeable difference is seen between the difference of the 9 attributes as
well as in the 5 tests carried out on individual attributes using threads.
Therefore, it is concluded that using threads is effective for similar cases
where data analysis is required to be more optimal.

## Observations on the data:

Graphic attribute Sex:

Interpretation: The pie chart shows us that the holder of the document was
44.16% male, 43.36% female, and 12.48% not applicable. in the registry of
vital events with 500,000 data. From which it can be deduced that it was men
and women who were the most holders of the minutes while it did not apply, it
was very low.

Graph attribute Done:

Interpretation: The pie chart shows us that there were 61.54% births, 25.98%
deaths and 12.48% marriages. in the registry of vital events with 500,000
data. From which it can be deduced that there were more births than deaths and
marriages and it is also visualized that the percentage of marriages is half
of the deaths.

## Discussion:

It was proposed to reduce the processing time of large amounts of data in a
case of three measures of semantic similarity in the field of biology using
parallel processing, as a result a reduction in time was achieved that
gradually increases with the increase in data.[1]

Visual information data processing is important in homes when making a
comparison of which language adjusts to the processing conditions, as a result
it was obtained that the code is more effective in Python processing images in
parallel using threads since it receives the same image but processing it on
different threads.[5]

In general, data processing also usually includes queries, algorithms for
query processing of context-free routes, this can include an endless number of
symbols which reflects execution time. It was evaluated through tests in Go
and Python, a biology database was used and as a result, performance gains
were demonstrated when using threads and not only working with a single
process.[12]

Projecting billions and millions of longitude and latitude location data onto
hundreds of thousands of highly irregular population census block polygons is
computationally challenging and can be implemented in any scripting language
such as Python. the tests were divided into simple and fast as a result it is
appreciated that the fast approach exploits the optimizations through threads,
achieving quick integration of location and demographic data.[9]

If it is proposed to implement and evaluate two error propagation metrics
first error propagation through structured data second location corruption
fraction both software and hardware, both metrics were evaluated by running
single threads and multiple threads, which demonstrated that using
multithreading is more effective detecting both metrics.[14]

It is observed that the efficiency of the journey optimizes and speeds up the
processing time of large amounts of data in the broader field of analysis
through open source projects and industry projects, using parallel processing
and as a result a potential reduction was achieved. of time that is
beneficially increased by the increase in data mining.[6]

It has been carried out in order to process the large data on the Markov chain
Monte Carlo Algorithm, it was possible to observe the speed of time was
considerable and satisfactory at the moment of parallelization processing
through the implementation of threads that optimize and is the best handling
of reads at compile time.[10]

an explosion of software tools have been created to process these data files
in order to facilitate this, a library was produced from the original SAMtools
implementation, with a focus on performance and robustness performance has
improved considerably, with a BAM read and write loop that runs 5 times faster
and a conversion from BAM to SAM 13 times faster and at the time of processing
the large data it was possible to observe the speed of time was satisfactory
at the time of parallelization processing through the threads implementation
that optimizes and is the best when compiling.[4]

Especially in data development parsing algorithms While true when processing
data and compiling code with threads they provide powerful tools to speed up
the development process, on the other hand, by specializing the interpreter
for a given script. Loops in the specialized code can then be parallelized to
further improve performance time.[13]

A code was implemented in order to optimize the time in parallel with threads
when processing a large amount of data that have made possible important
advances in computer science thanks to the computing power and the consequent
impossibility of analyzing and understanding the results with the current
techniques. The Neurolytics implementation is primarily focused on simulation
neuroscience to help scientists structure and understand the results of their
experiments in a fast and efficient way. [15]

## References

  * [1] A.M. Almasoud, H.S. Al-Khalifa, and A.S. Al-Salman.  Handling big data scalability in biological domain using parallel and distributed processing: A case of three biological semantic similarity measures.  BioMed Research International, 2019, 2019.  cited By 2.
  * [2] Alexander David Andrango Ayo.  Implementación de un sistema de detección de monedas por visión artificial: Sistema de detección de patrones simples de monedas.  B.S. thesis, Quito: EPN, 2022., 2022.
  * [3] Matías S Ávalos.  1\. resolución de problemas.
  * [4] James K Bonfield, John Marshall, Petr Danecek, Heng Li, Valeriu Ohan, Andrew Whitwham, Thomas Keane, and Robert M Davies.  Htslib: C library for reading/writing high-throughput sequencing data.  Gigascience, 10(2):giab007, 2021.
  * [5] A. Bustamante, L.M. Belmonte, R. Morales, A. Pereira, and A. Fernández-Caballero.  Video processing from a virtual unmanned aerial vehicle: Comparing two approaches to using opencv in unity.  Applied Sciences (Switzerland), 12(12), 2022.  cited By 0.
  * [6] Fabian Heseding, Willy Scheibel, and Jürgen Döllner.  Tooling for time-and space-efficient git repository mining.  arXiv preprint arXiv:2205.01351, 2022.
  * [7] Mario Ibáñez Bolado et al.  Paralelización de técnicas de toma de decisiones en empresas de acuicultura.  2022\.
  * [8] Oriel Keba Manianga et al.  Aplicación de cómputo en paralelo basado en threads en el cálculo de métricas usadas en el análisis armónico.  2021\.
  * [9] J. Kepner, A. Kipf, D. Engwirda, N. Vembar, M. Jones, L. Milechin, V. Gadepally, C. Hill, T. Kraska, W. Arcand, D. Bestor, W. Bergeron, C. Byun, M. Hubbell, M. Houle, A. Kirby, A. Klein, J. Mullen, A. Prout, A. Reuther, A. Rosa, S. Samsi, C. Yee, and P. Michaleas.  Fast mapping onto census blocks.  2020\.  cited By 2.
  * [10] Botao Li, Synge Todo, Anthony C Maggs, and Werner Krauth.  Multithreaded event-chain monte carlo with local times.  Computer Physics Communications, 261:107702, 2021.
  * [11] Antonio López García and Jaime Urquiza-Fuentes.  Experiencias de uso del pseudocódigo y java en la enseñanza de programación en ciclos formativos y bachillerato.  2022\.
  * [12] C.M. Medeiros, M.A. Musicante, and U.S. Costa.  Querying graph databases using context-free grammars.  Journal of Computer Languages, 68, 2022.  cited By 0.
  * [13] Brandon Neth and Michelle Mills Strout.  Automatic parallelization of irregular x86-64 loops.  In 2019 IEEE/ACM International Symposium on Code Generation and Optimization (CGO), pages 266–266. IEEE, 2019.
  * [14] Z. Ozturk, H.R. Topcuoglu, and M.T. Kandemir.  Studying error propagation on application data structure and hardware.  Journal of Supercomputing, 78(17):18691–18724, 2022.  cited By 0.
  * [15] Judit Planas, Fabien Delalondre, and Felix Schürmann.  Accelerating data analysis in simulation neuroscience with big data technologies.  In International Conference on Computational Science, pages 363–377. Springer, 2018.