<figure><img src="image_1.png" /></figure>

HAL
open science

# Parameter identification through gradient flow on latent variables

Muriel Boulakia, Haibo Liu, Damiano Lombardi

► To cite this version:

Muriel Boulakia, Haibo Liu, Damiano Lombardi. Parameter identification through gradient flow on latent variables. 2023. hal-04364114v2

HAL Id: hal-04364114

https://inria.hal.science/hal-04364114v2

Preprint submitted on 8 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Parameter identification through gradient flow on latent variables

## Muriel Boulakia

Université Paris-Saclay, UVSQ, CNRS

Laboratoire de Mathématiques de Versailles, 45 avenue des Etats-Unis, 78035, Versailles, France

E-mail: muriel.boulakia@uvsq.fr

## Haibo Liu

Sorbonne Université, Inria, CNRS & NOTOCORD part of INSTEM

Laboratoire Jacques-Louis Lions (LJLL), 4 Pl. Jussieu, 75005, Paris, France

E-mail: haibojobapply@gmail.com

## Damiano Lombardi

Sorbonne Université, Inria, CNRS

Laboratoire Jacques-Louis Lions (LJLL), 4 Pl. Jussieu, 75005, Paris, France

E-mail: damiano.lombardi@inria.fr

**Abstract.** In this article, we consider a system of parametric ordinary differential equations (ODE) which involves unknown parameters and we seek to identify the values of the parameters associated with a given measurement. To do so, we place ourselves within the fairly usual framework that this single measurement is in fact taken from a population of data and we therefore want to take advantage of the statistical knowledge about the population to regularize the classical minimization problem associated with our identification problem. In the method that we propose and that we call the Latent Variable Gradient Flow method, the data set is represented by an autoencoder neural network which allows to associate with each element of the data set a latent variable. Then, introducing a non-linear mapping between the parameter space and the latent variable space allows to convexify the cost function and to demonstrate convergence properties. These properties are numerically illustrated by different tests on ODE systems.

<footnote>*Author to whom correspondence should be addressed. Electronic mail: damiano.lombardi@inria.fr</footnote>
# 1. Introduction and problem statement

## 1.1. Introduction

In the wide variety of inverse problems, our article focuses on the framework where a statistical learning method is combined with the knowledge of a mathematical model in order to identify parameters of the model. In experimental sciences, such as chemistry, biology or environmental engineering, mathematical models are used to describe the behavior of the dynamical systems. Those mathematical models consist of systems of ordinary differential equations (ODEs). For the most part, these ODEs contain parameters that are associated with physical, biological, or other properties of the system. To assess the relevance of the model, it is necessary to understand the role of these parameters and to set their values in order to generate signals that are close to reality. In addition to parameters, the system is completed by the data of the initial state whose values may be unknown or uncertain. When faced with an experimental observation, identifying the unknown parameters and the initial state to closely match the experimental results is a key step in the development and understanding of ODEs [47]. However, this identification problem is a complex task, especially when the ODEs system consists of a large number of equations and parameters.

Let us mention that machine learning methods have been used to enable or improve the resolution of inverse problems, in many different contexts. We give here a non exhaustive but significant list of references which combine machine learning and inverse problems. Machine learning has been in particular used in the design of regularizing terms in optimization formulations of inverse problems, as done for instance in [1], [42], [21]. In [28], machine learning is used to perform data-driven modeling (by random features map) and combined with an ensemble Kalman filter to perform sequential data assimilation. In [27], the propagator is learned sequentially and used in an ensemble method. These are also examples in which machine learning is used to learn the model. In [50], a deep neural network is trained to reproduce the parameter to an observable map of a model in electro-physiology. This can be exploited to efficiently solve inverse problems. In [37], neural ODEs are used in order to learn distributions and in the more recent [26], diffusion maps and Langevin dynamics are used in order to enabling sampling from an unknown (eventually conditional) distribution. These are ways in which machine learning can be used in order to improve the resolution of Bayesian inverse problems. On this same line, Generative Adversarial Networks (GANs) have been studied in [43] and conditional variational autoencoders (CVAE) in [23]. At last, among numerous applications, a very active field of applications where deep learning is applied to inverse problems corresponds to image restoration and reconstruction [38], [40], [42] and [21].

To define the context we are interested in more precisely, let us introduce some notations which will be defined in more detail later on. Let a system be described by a mathematical model, in which the output $u$ is affected by a certain number of parameters and by the initial conditions. Let us denote by $\theta$ the vector of these parameters and initial conditions. To simplify the formulation, in what follows, we will describe in general $\theta$ as a model parameters vector,
even if it also includes the initial conditions. Then, we introduce the parameter-to-output map $\varphi$ which associates with the vector $\theta$ the output $u: u = \varphi(\theta)$.

The general problem reads as follows: given some measurements $u^*$ of the output, we try to estimate the model parameters $\theta^*$ such that $u^* = \varphi(\theta^*)$.

Classically, the estimation of model parameters ([4, 7]) is cast as an optimization problem. In particular, a cost function is defined based on the data misfit and the parameter $\theta$ is estimated by minimizing the discrepancy between the model observations $\varphi(\theta)$ and the actual measurement $u^*$ of the system. This formulation often leads to a non-linear non-convex possibly high-dimensional optimization problem. This could be dealt with in two ways:

(i) Local (often gradient based) optimization methods. These methods include Quasi-Newton methods [19, 25], the Gauss-Newton method (incorporating explicitly derived and efficiently computed sensitivity equations [9]), as well as other gradient-based methods [6]. The local gradient-based optimization methods might encounter convergence issues when the initial guess is far from the true parameter values [35].

(ii) Global optimization approach. For example, commonly utilized stochastic search algorithms include evolutionary computation, adaptive stochastic methods, and clustering methods. Other global optimization methods are the genetic algorithms, as proposed for instance in [54] for the determination of rate constants in heterogeneous reaction systems, and collocation methods, as proposed for instance in [5] and [10]. In many cases, stochastic search algorithms are used to select good initial guesses for starting either a Quasi-Newton method or a gradient-based type minimization. Global optimization methods are often very expensive from a computational point of view and might become prohibitive when the number of parameters is large [31, 57].

Another way to try to solve the problem is purely based on data. Let us assume that, in certain situations, we have examples of pairs of parameters and observations, $\{\theta^{(i)}, u^{(i)}\}_{1\le i\le N_d}$. For example, the work in [24] used a supervised convolutional machine-learning method to learn the physical model parameters. The advantage of this kind of model is that no prior knowledge about the system is required. The main disadvantage consists in the fact that we often need a large number of training pairs parameters-observations.

Recently, many efforts were made in order to try to exploit both the knowledge coming from the mathematical model and from sets of data, such as physics-informed neural networks (PINNs) mentioned in [20, 29, 34, 55], deep operator network (DeepONet) find in [36] and physics-informed DeepONet (PI-DeepONet) find in [53].

PINNs operate by taking time as input and by training the network to generate corresponding observations of the dynamic system. This training process leverages both the available data and the underlying dynamics of the system. When addressing inverse problems, the PINN's loss function integrates two components: the discrepancies between predicted observations and given data, and the residual of the differential equations, which depends on the parameters. Then, the network optimizes the parameters of the network and estimates the parameters of the dynamical system. PINNs might suffer from the so-called spectral
bias [46], which refers to the fact that it tends to prioritize learning the low-frequency part of the solution, not rendering properly multiple frequency scales. To overcome this issue, Fourier Features-Neural Networks (FF-NN) have been proposed to solve dynamical systems and perform inverse problems [48]. Another limitation of PINN-based methods is that they often require a large number of iterations and evaluations to converge to an accurate solution. This computational intensity is a drawback, especially when dealing with dynamic systems consisting of a large number of equations and parameters. In the continuity of PINN-based methods, the method that we propose in this paper seeks to leverage both available data and underlining dynamical systems. However, as we will see, our method provides an alternative method of PINN for which the computational demanding part corresponding to the training on the data set is done in a preliminary part.

About the data set, we focus here on a specific situation: we consider that we do not have data composed of pairs of parameters-observations but that we only have observations

$$ \{u^{(i)}\}_{1 \le i \le N_d} $$

consisting in a collection of measurements taken on a population or a set of experiments and we assume that our (novel) observation $u^*$ comes from the same population. This situation, which is quite common in a number of realistic applications, is by far more difficult than the one in which pairs of parameters-observations are available (as we have no direct information about the parameters anymore).

In this paper, we will propose a new method of parameter estimation which basically relies on a regularization of a classical parameter estimation by means of the data set. A representation of the data set by an autoencoder will allow us to describe an element of the data set by latent variables. Then, the method performs a non-linear mapping between the parameter space, where the cost function associated with the parameter estimation is non-convex, and the latent variable space, where the cost function is convex. Since our method can be described as a gradient flow for the latent variables of the autoencoder, it will be called the Latent Variables Gradient Flow (LVGF) method. From a theoretical point of view, we show that the LVGF method has convergence properties. In the first result (stated in Proposition 2.1), we establish the local linear convergence. In a noticeable way, this result holds in a rather general framework since the main assumptions are that the composition of the encoder with the parameter-to-output map $\varphi$ is a submersion and that its gradient is locally Hölder continuous. The global version of this convergence result (which is stated in Proposition 2.3) adopts similar but more limited assumptions, as it is adapted to the global framework. The method has been numerically tested for identification problems applied to ODE systems of variable complexity (Van der Pol model, FitzHugh-Nagumo model and Hindmarsh-Rose model). In all the numerical tests, the LVGF method improves the results obtained with a gradient descent method applied to the classical misfit functional. Whereas the performance of this latter method is limited by the convergence towards local minimizers, the LVGF method is able to outperform it by reaching the global minimizer.

In addition to the presentation and study of the LVGF method, this paper presents,
as a secondary task, a novel method to estimate the intrinsic dimension of a data set $U = \{u^{(i)}\}_{1\le i\le N_d}$. Identifying this intrinsic dimension is an important point in the implementation of our method since it allows us to set the value of the latent dimension of the autoencoder trained on the data set. The proposed method will be based on the use of autoencoders and a criterion based on the approximation of the Lipschitz constant of the encoding and decoding maps. For an ODE system satisfying observability and identifiability properties, this intrinsic dimension is related to the dimension of the set of parameters which are varying in the population.

The structure of the work is the following: in Sections 1.2 and 1.3, we state our parameter identification problem and present the notions of identifiability and observability. Moving on to Section 2, we present our parameter estimation method which we call the LVGF method, and study its convergence properties. Next, in Section 3, we present a new criteria based on autoencoders to estimate the intrinsic dimension of a submanifold. At last, in Section 4, we assess the LVGF method's ability to solve the parameter estimation problem on the dynamical systems presented in the first section and compare its performance with the classical minimization method both for noiseless data and noisy data.

## 1.2. Presentation of the parameter identification problem

Let us define a continuous function $f: \mathbb{R}^N \times \mathbb{R}^q \to \mathbb{R}^N$. We consider an autonomous and parameterized system of ODEs of the form:

$$ \begin{cases} \dot{x}(t) = f(x(t), \zeta), \forall t \in (0, T), \\ x(0) = x_0 \end{cases} \qquad (1.1) $$

where $x: [0, T] \to \mathbb{R}^N$ is a real vector valued function, $\zeta \in \mathbb{R}^q$ represents a vector of constant parameters and $x_0 \in \mathbb{R}^N$ is the initial state. We assume that, for every $x_0 \in \mathbb{R}^N$ and $\zeta \in \mathbb{R}^q$, system (1.1) admits a unique solution $x$ in $C^1([0, T])^N$.

For this system, we consider that the values of the parameters and the initial conditions are unknown or uncertain and we are interested in the simultaneous identification of parameters and initial data. We denote by $\theta = (x_0, \zeta) \in \mathbb{R}^m$ with $m = N + q$ the vector that we want to identify and look for $\theta$ in a subset of $\mathbb{R}^m$ that we denote by $\Theta$.

The overall goal of our study is to identify the value of $\theta \in \Theta \subset \mathbb{R}^m$, by exploiting some discrete measurements that we denote by $u \in \mathbb{R}^n$ of the solution $x$ of system (1.1). We denote by $H$ the observation operator that models the measurement procedure:

$$ H: C^1([0, T])^N \to \mathbb{R}^n. \qquad (1.2) $$

For the sake of simplicity in the presentation, we assume that $H$ only gives one scalar information (for instance, one component of the vector valued function $x$) discretized in time.

In addition, we define the map which associates the measurement with the unknown parameters and initial data:

$$ \begin{cases} \varphi: \mathbb{R}^m \to \mathbb{R}^n \\ \theta \to H(x(\theta, \cdot)) \end{cases} \qquad (1.3) $$
where $x(\theta, \cdot)$ is the solution of system (1.1) associated with $\theta = (x_0, \zeta)$.

The classical deterministic identification problem can be formulated in the following way: we assume that a datum $u^* \in \mathbb{R}^n$ is given and we would like to identify $\theta^* \in \Theta \subset \mathbb{R}^m$ such that

$$ \varphi(\theta^*) = u^*. \qquad (1.4) $$

Most of the time, this problem is numerically studied by reformulating it as an optimization problem. In a classical way, the functional to minimize can be for instance of the form:

$$ J(\theta) = \|\varphi(\theta) - u^*\|_n^2 \qquad (1.5) $$

where $\| \cdot \|_n$ corresponds to the Euclidean norm in $\mathbb{R}^n$.

Such a problem is known to be cumbersome in a number of situations, due to the fact that it is usually non-convex and therefore classical gradient descent methods may fail to converge (we will give a practical example of those situations at the end of the next section).

In the present work, we consider that, in addition to the map $\varphi$ and the datum $u^*$, we have access to a set of $N_d \in \mathbb{N}^*$ data, denoted $U = \{u^{(i)}\}_{1 \le i \le N_d} \in \mathbb{R}^n$ which $u^*$ is a part. We assume that these data come from a population (of individuals, or experiments), characterized by a certain (unknown) distribution of parameters $\theta$. This setting is justified by the fact that it is indeed common to have access not to a single measurement but to a large number of data. Therefore, our objective is to study how we can take advantage of this whole set of measurements $U$ in the identification of parameters associated with a single measurement $u^*$ and if adding information coming from the knowledge of the set $U$ can allow to construct a regularisation for the parameter estimation problem. This idea will be the starting point of the LVGF method presented in 2.2.

### 1.3. Identifiability and observability of the inverse problem

Before testing numerical methods for the resolution of the identification problem (1.4), it is essential to ensure the uniqueness of a solution to (1.4) or, in other words, the injectivity of $\varphi$ in $\Theta$. Not having this type of theoretical property can in fact explain numerical convergence issues independently of the numerical method chosen.

Moreover, contrary to Tykhonov methods that can overcome a lack of uniqueness in the identification problem through the add of a prior, in our case the only information that we add to regularize our minimization problem comes from the knowledge of a data set. Thus, the underdetermined nature of the problem would remain despite the regularization process.

Since $\theta$ is composed of both model parameters and initial data, this identification property is related to two distinct notions in inverse problems: first, the notion of identifiability which refers to the fact that model parameters are uniquely determined by the output and, second, the notion of observability which refers to the fact that the initial conditions are uniquely determined by the output.
Let us notice that the identifiability and observability of the problem introduced above can be reduced to the observability of the following autonomous augmented system (as presented for instance in [17] and [39]):

$$ \begin{cases} \dot{X}(t) = F(X(t)), \forall t \in (0, T), \\ X(0) = (x_0, \zeta) = \theta \end{cases} \qquad (1.6) $$

where $\theta$ is taken in $\Theta$, $X = (x, \zeta) : [0, T] \to \mathbb{R}^m$ is a regular function and $F$ is given by

$$ \begin{aligned} F : \mathbb{R}^m & \to \mathbb{R}^m \\ (x, \zeta) & \to (f(x, \zeta), 0). \end{aligned} \qquad (1.7) $$

Contrary to the observation operator given by (1.2) and adapted to the numerical framework where the solution is discretized in time, we consider in the theoretical framework that we have access to measurements at each time in $[0, T]$ and we denote by

$$ h : C^1([0, T])^N \to C^1([0, T]) \qquad (1.8) $$

the observation map which associates the time-dependent measurement with the solution $X$ of the ODE system (1.6).

Numerous works address the question of identifiability and observability for autonomous ODE systems and there is a wide variety of methods for doing so. We can quote [2], [8], [16], [51] and [52] among many references. Following these last two references, we are interested by the structural observability property for system (1.6) (or by the structural identifiability and observability property for system (1.1)) which means that for almost all initial conditions $X_1(0)$ and $X_2(0)$ in $\Theta$

$$ h(X_1) = h(X_2) \text{ in } [0, T] \Rightarrow X_1(0) = X_2(0). $$

**Remark 1.1.** Since the notions of identifiability and observability are in general studied for autonomous ODE systems, we have chosen to present our work in this framework. We could also have considered non-autonomous ODE systems and, in this more general framework, the methods and results presented in this work still apply.

## 2. Latent Variables Gradient Flow (LVGF) for parameter estimation

### 2.1. Principle of the LVGF method

In what follows, we represent the given data set $U$ as a sampling of an embedded submanifold $U \subset \mathbb{R}^n$. We consider $u^* \in U$ and $\theta^* \in \Theta$ which are related by $u^* = \varphi(\theta^*)$ and we will explain the main principles of the proposed LVGF method to numerically identify the value of $\theta^*$.

The first ingredient of the LVGF method relies on the construction of an autoencoder approximating the given data set $U$. It consists in an encoder map $\Psi : \mathbb{R}^n \to \mathbb{R}^p$ which compresses an input data $u \in U$ into a latent variable $\alpha = \Psi(u) \in \mathbb{R}^p$ with $p < n$ and a
decoder map $\Psi^\dagger : \mathbb{R}^p \to \mathbb{R}^n$ which recreates an approximation of the input data from the latent variable. For all $u \in \mathcal{U}$, $\Psi^\dagger \circ \Psi(u)$ is a good approximation of $u$. The set-up of the autoencoder will be specified in Section 3.1.

The rationale behind the introduction of the autoencoder is the following: the latent variable $\alpha \in \mathbb{R}^p$ is, at the same time, observable and related to the parameter $\theta \in \Theta$. Indeed, thanks to the encoder $\Psi$, we can easily compute the latent variable $\alpha^* = \Psi(u^*)$ associated with $u^*$. This implies that we have the following relationship between $\theta^*$ and $\alpha^*$:

$$g(\theta^*) = \alpha^* \qquad (2.1)$$

where $g$ is given by $g = \Psi \circ \varphi : \mathbb{R}^m \to \mathbb{R}^p$.

This equality naturally leads to a rewriting of the classical identification problem given by (1.4). Therefore, a first optimization problem based on the encoder map could be the following : minimize with respect to $\theta \in \Theta$ the functional

$$\|g(\theta) - \alpha^*\|_p^2.$$

However, due to the lack of convexity, classical numerical methods associated with this functional would have the same convergence issues as (1.5). Therefore, the second ingredient of our method is to introduce the following convex functional in the latent variable space

$$\|\alpha^* - \alpha\|_p^2 \qquad (2.2)$$

and to design a numerical method based on the gradient descent method associated with this functional and the use of the functional $g$ to move from the latent variable space to the parameter space. Therefore, since this method can be viewed as the discrete counterpart of a gradient flow dynamics in the latent variable $\alpha$, this has motivated the method name: Latent Variable Gradient Flow (LVGF).

## 2.2. Description of the LVGF method

In what follows, we assume that $g$ is a $C^1$ submersion in $\Theta$ that is to say:

$$g \in C^1(\Theta) \text{ and } \mathrm{rank} \nabla g(\theta) = p, \forall \theta \in \Theta. \qquad (2.3)$$

Under these hypotheses, we can consider the generalized inverse (also called the Moore-Penrose pseudoinverse [41]) of $\nabla g(\theta)$ defined for a matrix $M$, by $M^\dagger := M^T (MM^T)^{-1}$.

The LVGF algorithm is based on the following iterative process: let a step $0 < s < 1$ and an initial value $\theta_0 \in \Theta$ be given. We denote by $\alpha_0$ the initial value associated with $\theta_0$, that is $\alpha_0 = g(\theta_0)$. Then, we define the sequence $(\theta_k)_{k \in \mathbb{N}}$ iteratively by: for $k \in \mathbb{N}$

$$\left\{ \begin{array}{l} \Delta\theta_k = sM_k^\dagger(\alpha^* - \alpha_k) \text{ where } M_k = \nabla g(\theta_k) \\ \theta_{k+1} = \theta_k + \Delta\theta_k \\ \alpha_{k+1} = g(\theta_{k+1}). \end{array} \right. \qquad (2.4)$$
By definition, $\Delta\theta_k$ is thus the unique vector of smallest norm solution in $\mathbb{R}^m$ of the system

$$M_k \Delta\theta_k = s(\alpha^* - \alpha_k).$$

Therefore, in the right hand side, we retrieve the classical increment of the sequence defined by the gradient descent method applied to the functional (2.2).

Let us notice that the hypothesis (2.3) on the rank of the gradient of $g$ can be easily computationally verified. Moreover, since this hypothesis may fail if the hyperparameter $p$ is taken too large with respect to the intrinsic dimension of the manifold, it is important to accurately identify this intrinsic dimension. By using some recent results in approximation theory, we will present a new method to identify it in Section 3.2.

### 2.3. Convergence properties of the LVGF method

In this subsection, we give two convergence results for the sequence $(\theta_k)$ defined by (2.4). The first result is stated in Proposition 2.1 and gives a local convergence result. Compared to classical convergence results for the gradient descent method established under convexity assumptions, our result holds within a general framework since the only requirement in addition to (2.3) is that the gradient of $g$ is locally Hölder continuous and that $\theta^*$ belongs to the interior of $\Theta$. The second result is stated in Proposition 2.3 and gives a global convergence result in the case when $\Theta = \mathbb{R}^m$. To establish this strong result, we assume that the gradient of $g$ is globally Hölder continuous and that $g$ satisfies (2.12). Both results state that the sequence $(\theta_k)$ converges to some limit denoted by $\tilde{\theta}$. In addition to these local and global results, Corollaries 2.2 and 2.4 establish that, under additional injectivity hypotheses on $\Psi$ and $\varphi$, the limit $\tilde{\theta}$ coincides with the parameter $\theta^*$.

Let us first state a local convergence result on the sequence $(\theta_k)$.

**Proposition 2.1.** Let $\theta^*$ satisfy (2.1) and belong to $\text{Int}(\Theta)$. We assume that $g$ satisfies (2.3) and $\nabla g$ is locally Hölder continuous. Then, there exists $\eta > 0$ such that, for $s > 0$ small enough, for all $\theta_0 \in \bar{B}(\theta^*, \eta)$ the sequence $(\theta_k)_{k \in \mathbb{N}}$ given by (2.4) converges to some $\tilde{\theta} \in \Theta$ which satisfies:

$$\alpha^* = g(\tilde{\theta}). \qquad (2.5)$$

Moreover, there exists a constant $C > 0$ such that:

$$\|\tilde{\theta} - \theta_k\|_m \le \frac{C}{s}(1 - \frac{s}{2})^k.$$

*Proof.* Let $M > 0$ be given such that the closed ball $\bar{B}(\theta^*, M)$ of center $\theta^*$ and of radius $M$ is included in $\Theta$. By hypothesis, $\nabla g$ is Hölder continuous in $\bar{B}(\theta^*, M)$, thus there exist $\beta > 0$ and $\tilde{C} > 0$ such that, for all $\theta_1, \theta_2 \in \bar{B}(\theta^*, M)$

$$\|\nabla g(\theta_1) - \nabla g(\theta_2)\|_p \le \tilde{C} \|\theta_1 - \theta_2\|_m^\beta$$

We also denote by $C_L$ the Lipschitz constant of $g$ in $\bar{B}(\theta^*, M)$. Since $\bar{B}(\theta^*, M)$ is compact, using hypothesis (2.3), there exists $C_M > 0$ such that, for all $\theta \in \bar{B}(\theta^*, M)$

$$\|\nabla g(\theta)^{\dagger}\| \le C_M. \qquad (2.6)$$
In what follows, we take $0 < s < 1$ such that

$$ \tilde{C}C_M^{1+\beta}C_L^{\beta}M^{\beta}s^{\beta} \le \frac{1}{2} \qquad (2.7) $$

and next we take $\eta > 0$ such that

$$ \eta \left( 1 + \frac{2C_M C_L}{s} \right) \le M. \qquad (2.8) $$

Let us take $\theta_0 \in \bar{B}(\theta^*, \eta)$. We will show by induction that the sequences $(\alpha_k)$ and $(\theta_k)$ defined by (2.4) satisfy

$$ \left\{ \begin{array}{l} \|\alpha^* - \alpha_k\|_p \le (1 - \frac{s}{2})^k \|\alpha^* - \alpha_0\|_p \\ \\ \|\theta^* - \theta_k\|_m \le \|\theta^* - \theta_0\|_m + C_M \|\alpha^* - \alpha_0\|_p \sum_{j=0}^{k-1} (1 - \frac{s}{2})^j. \end{array} \right. \qquad (2.9) $$

This property obviously holds for $k = 0$. Now we assume that it holds for $k$ and want to prove that it holds for $k + 1$. By hypothesis, we have

$$ \|\theta^* - \theta_k\|_m \le \eta + \frac{2C_M}{s} \|\alpha^* - \alpha_0\|_p \le \eta + \frac{2C_M C_L}{s} \|\theta^* - \theta_0\|_m \le \eta \left(1 + \frac{2C_M C_L}{s}\right) \le M $$

according to (2.8). Thus $\theta_k$ belongs to $\bar{B}(\theta^*, M)$ and satisfies (2.6). This implies that

$$ \|\Delta\theta_k\|_m \le sC_M \|\alpha^* - \alpha_k\|_p \le C_M \|\alpha^* - \alpha_k\|_p \le C_M \left(1 - \frac{s}{2}\right)^k \|\alpha^* - \alpha_0\|_p. $$

Therefore, we deduce that

$$ \|\theta^* - \theta_{k+1}\| \le \|\theta^* - \theta_k\| + \|\Delta\theta_k\| \le \|\theta^* - \theta_0\|_m + C_M \|\alpha^* - \alpha_0\|_p \sum_{j=0}^{k} \left(1 - \frac{s}{2}\right)^j. \qquad (2.10) $$

In particular, we deduce that $\theta_{k+1}$ belongs to $\bar{B}(\theta^*, M)$. This allows us to write that

$$ \alpha_{k+1} = g(\theta_k + \Delta\theta_k) = g(\theta_k) + M_k \Delta\theta_k + \xi_k \qquad (2.11) $$

where $\xi_k \in \mathbb{R}^p$ satisfies

$$ \|\xi_k\|_p \le \tilde{C} \|\Delta\theta_k\|_m^{1+\beta} \le \tilde{C} C_M^{1+\beta} s^{1+\beta} \|\alpha^* - \alpha_k\|_p^{1+\beta} \le \frac{s}{2} \|\alpha^* - \alpha_k\|_p $$

where we have used (2.7). We deduce that

$$ \alpha^* - \alpha_{k+1} = \alpha^* - \alpha_k - M_k \Delta\theta_k - \xi_k = (1 - s)(\alpha^* - \alpha_k) - \xi_k $$

and therefore

$$ \|\alpha^* - \alpha_{k+1}\|_p \le (1 - \frac{s}{2}) \|\alpha^* - \alpha_k\|_p \le (1 - \frac{s}{2})^{k+1} \|\alpha^* - \alpha_0\|_p. $$
Combining this estimate with (2.10), we deduce that property (2.9) holds for $k+1$ and this concludes the proof by induction.

From the first estimate of (2.9), we deduce that the sequence $(\theta_k)_{k \in \mathbb{N}}$ satisfies

$$ \|\theta_{k+1} - \theta_k\|_m = \|\Delta\theta_k\|_m \le C_M \|\alpha^* - \alpha_k\|_p \le C_M \left(1 - \frac{s}{2}\right)^k \|\alpha^* - \alpha_0\|_p. $$

Thus $(\theta_k)_{k \in \mathbb{N}}$ is a Cauchy sequence and it converges to some limit that we denote by $\tilde{\theta} \in \Theta$. Moreover, we have

$$ \|\tilde{\theta} - \theta_k\|_m \le \frac{2C_M}{s} \left(1 - \frac{s}{2}\right)^k \|\alpha^* - \alpha_0\|_p. $$

At last, passing to the limit in the expression $\alpha_k = g(\theta_k)$, we deduce the formula (2.5). □

A consequence of this proposition is that, under additional injectivity hypotheses on $\Psi$ and $\varphi$, the sequence $(\theta_k)$ converges to $\theta^*$:

**Corollary 2.2.** Under the same hypotheses as Proposition 2.1, we assume in addition that $\Psi: \mathcal{U} \to \mathbb{R}^p$ and $\varphi: \Theta \to \mathbb{R}^n$ are one-to-one and that $\tilde{\theta} \in \varphi^{-1}(\mathcal{U})$. Then, the sequence $(\theta_k)_{k \in \mathbb{N}}$ obtained in Proposition 2.1 converges to $\theta^* \in \Theta$.

This corollary directly comes from the fact that, by definition of $\alpha^*$ and according to equation (2.5), the limit $\tilde{\theta}$ of the sequence $(\theta_k)_{k \in \mathbb{N}}$ satisfies:

$$ \Psi \circ \varphi(\tilde{\theta}) = \Psi \circ \varphi(\theta^*). $$

Thus, under the hypotheses of the corollary, $\tilde{\theta} = \theta^*$.

Let us add a few comments on the additional assumptions of Corollary 2.2 (this comment also applies to Corollary 2.4 below). First, we remark that the hypothesis that $\Psi: \mathcal{U} \to \mathbb{R}^p$ is one-to-one is naturally related to the encoding properties of this function for which a given latent variable can correspond to only one vector in $\mathcal{U}$. The hypothesis of injectivity of $\varphi: \mathbb{R}^m \to \mathbb{R}^n$ is also a natural assumption related to the identifiability and observability discussed in Section 1.3.

For the global convergence result stated below, we consider the case where the parameter space $\Theta$ is the full space $\mathbb{R}^m$. Moreover, we assume, in addition to (2.3), that there exists a constant $\bar{C} > 0$ such that, for all $\theta \in \mathbb{R}^m$

$$ \|\nabla g(\theta)^{\dagger}\| \le \bar{C}. \qquad (2.12) $$

Let us mention that this bound is related to the property (2.3). Indeed, if we were to place ourselves in a compact, property (2.12) would be a direct consequence of (2.3). However, to ensure that the sequence stays in a compact set, it is necessary to add a projection step (we refer to Remark 2.5) and the convergence study still needs to be adapted to this modified algorithm.
**Proposition 2.3.** We assume that $\Theta = \mathbb{R}^m$. Moreover, we assume that $g$ satisfies the hypotheses (2.3) and (2.12) and that $\nabla g$ is Hölder continuous in $\mathbb{R}^m$. Then, for $s > 0$ small enough, the sequence $(\theta_k)_{k \in \mathbb{N}}$ defined by (2.4) converges to some $\tilde{\theta} \in \mathbb{R}^m$ which satisfies

$$ \alpha^* = g(\tilde{\theta}). $$

Moreover, there exists a constant $C > 0$ depending on $\alpha^* - \alpha_0$ and $\bar{C}$ such that

$$ \|\tilde{\theta} - \theta_k\|_m \le \frac{C}{s} \left(1 - \frac{s}{2}\right)^k $$

*Proof.* The arguments of the proof are similar to the ones in the proof of Proposition 2.1, except that it is no longer necessary to ensure that the sequence $(\theta_k)$ stays in a given ball, so we need only part of the proof in what follows.

Let $k \in \mathbb{N}$ be given. According to (2.4), we have

$$ \alpha^* - \alpha_{k+1} = \alpha^* - g(\theta_k + \Delta\theta_k). $$

Since $g \in C^1(\mathbb{R}^m)$ and since the gradient of $g$ is Hölder continuous with a Hölder exponent that we denote $\beta > 0$, we can write

$$ g(\theta_k + \Delta\theta_k) = g(\theta_k) + M_k \Delta\theta_k + \xi_k = \alpha_k + s(\alpha^* - \alpha_k) + \xi_k $$

where $\xi_k \in \mathbb{R}^p$ satisfies

$$ \|\xi_k\|_p \le \tilde{C} \|\Delta\theta_k\|_m^{1+\beta} \le \tilde{C} \bar{C}^{1+\beta} s^{1+\beta} \|\alpha^* - \alpha_k\|_p^{1+\beta} $$

where we have used (2.12) and denoted by $\tilde{C}$ the Hölder constant of $\nabla g$. This implies that

$$ \|\alpha^* - \alpha_{k+1}\|_p = \|(1-s)(\alpha^* - \alpha_k) - \xi_k\|_p \\ \le (1-s)\|\alpha^* - \alpha_k\|_p + \tilde{C}\bar{C}^{1+\beta}s^{1+\beta}\|\alpha^* - \alpha_k\|_p^{1+\beta} $$

Assume now that the step $0 < s < 1$ is chosen such that

$$ \tilde{C}\bar{C}^{1+\beta}s^{\beta}\|\alpha^* - \alpha_0\|_p^{\beta} \le \frac{1}{2}. $$

Then, an argument by induction allows us to deduce from the previous inequality that, for all $k \in \mathbb{N}$

$$ \|\alpha^* - \alpha_k\|_p \le \|\alpha^* - \alpha_0\|_p. $$

Therefore, we deduce that, for all $k \in \mathbb{N}$

$$ \|\alpha^* - \alpha_{k+1}\|_p \le \left(1 - \frac{s}{2}\right) \|\alpha^* - \alpha_k\|_p. $$

So the sequence $(\alpha_k)_{k \in \mathbb{N}}$ linearly converges to $\alpha^*$ and

$$ \|\alpha^* - \alpha_k\|_p \le \left(1 - \frac{s}{2}\right)^k \|\alpha^* - \alpha_0\|_p. $$
At last, following the end of the proof of Proposition 2.1, we obtain that $(\theta_k)_{k \in \mathbb{N}}$ converges to some $\tilde{\theta}$ and that

$$ \|\tilde{\theta} - \theta_k\|_m \le \frac{2\bar{C}}{s} \left(1 - \frac{s}{2}\right)^k \|\alpha^* - \alpha_0\|_p $$

□

In a similar way as Corollary 2.2, injectivity hypotheses allow to identify the limit of the sequence $(\theta_k)_{k \in \mathbb{N}}$ as the parameter $\theta^* \in \mathbb{R}^m$.

**Corollary 2.4.** Under the same hypotheses as Proposition 2.3, we assume in addition that $\Psi : \mathcal{U} \to \mathbb{R}^p$ and $\varphi : \mathbb{R}^m \to \mathbb{R}^n$ are one-to-one and that $\tilde{\theta} \in \varphi^{-1}(\mathcal{U})$. Then, for $s > 0$ small enough, the sequence $(\theta_k)_{k \in \mathbb{N}}$ defined by (2.4) converges to $\theta^* \in \mathbb{R}^m$.

Since the algorithm given by equation (2.4) does not allow to ensure that the property $\tilde{\theta} \in \varphi^{-1}(\mathcal{U})$ is satisfied, in practice, we correct it by adding a projection procedure. Let us denote by $\mathcal{P} : \mathbb{R}^m \to \mathbb{R}^m$ a projection operator on $\varphi^{-1}(\mathcal{U})$. For a given $\theta_{-1} \in \mathbb{R}^m$, we compute $\theta_0 = \mathcal{P}(\theta_{-1}) \in \varphi^{-1}(\mathcal{U})$ and replace system (2.4) by

$$ \left\{ \begin{array}{l} \Delta\theta_k = sM_k^\dagger(\alpha^* - \alpha_k) \quad \text{where } M_k = \nabla g(\theta_k) \\ \theta_{k+1} = \mathcal{P}(\theta_k + \Delta\theta_k) \\ \alpha_{k+1} = g(\theta_{k+1}). \end{array} \right. \qquad (2.13) $$

For this new algorithm, further investigations would be needed to determine the convergence properties of this algorithm.

**Remark 2.5.** In most realistic situations, $\Theta$ is not the whole space $\mathbb{R}^m$ but a compact set in $\mathbb{R}^m$. The method can be written in a parametric domain (say, a box) containing this domain, and such that $\theta^*$ is an interior point of the domain. In the case in which the iterates are all interior points, the result presented holds true without any modification. If some iterates fall out of the domain, they can be reprojected into it.

**Remark 2.6.** The LVGF method comes with a way to verify, a posteriori, whether or not the model solution $\varphi(\theta)$ belongs to the population. This can be done by exploiting the decoder $\psi^\dagger$. Let us consider: $\|\psi^\dagger \circ \psi \circ \varphi(\theta) - \varphi(\theta)\|_n^2$. If this quantity is significantly larger than the error obtained on the validation set when training the autoencoder, it implies that the obtained solution is far from the ones belonging to the population.

## 2.4. Numerical implementation

For the numerical implementation of the method described by (2.13), it is necessary to describe how the projection on $\varphi^{-1}(\mathcal{U})$ is defined and computed. By noticing that the elements $\theta$ which belong to $\varphi^{-1}(\mathcal{U})$ satisfy

$$ \Psi^\dagger \circ \Psi \circ \varphi(\theta) - \varphi(\theta) = 0, $$

we introduce the functional

$$ E(\theta) = \|\Psi^\dagger \circ \Psi \circ \varphi(\theta) - \varphi(\theta)\|_n^2 \qquad (2.14) $$
and define the projection $\mathcal{P}(\theta)$ as the solution of the minimisation of $E$ starting from the initial data $\theta$. The minimization of the function $E$ is achieved by using the gradient descent method.

In practice, to initialize the sequence $(\theta_k)_{k \in \mathbb{N}}$ with a $\theta_0$ such that $\theta_0 \in \varphi^{-1}(\mathcal{U})$, the gradient descent method may fail if we start with a value far from the submanifold $\varphi^{-1}(\mathcal{U})$. So, we consider different initial values for $\theta_0$, run for each initial value the gradient descent method and keep a $\theta_0$ for which the optimization procedure gives a value of $E$ close to zero. In addition, we do not run the projection step at each iteration but only if $\theta_k$ begins to move away from $\varphi^{-1}(\mathcal{U})$ which is tested by comparing the value $E(\theta_k)$ to a threshold value. The details of the LVGF method are shown in Algorithm 1.

<caption>**Algorithm 1** The LVGF algorithm</caption>

1: **Input:** $N_{max}$ = maximum number of iterations; $\sigma$ = tolerance for $||\alpha^* - \alpha||_p$; $S = [s_1, s_2, ..., s_n]_{1 \le i \le n}$
2: start with a random initial guess $\theta_0$
3: **if** $E(\theta_0) > \epsilon$ **then**
4: run the projection $\theta = \mathcal{P}(\theta_0)$
5: **else**
6: $\theta = \theta_0$
7: **end if**
8: compute $\alpha = \Psi(\varphi(\theta))$
9: **while** $j \le N_{max}$ and $||\alpha^* - \alpha||_p \ge \sigma$ **do**
10: compute $M = [\nabla\Psi(\varphi(\theta))][\nabla\varphi(\theta)]$
11: compute $\Delta\theta = [M]^\dagger(\alpha^* - \alpha)$
12: select $1 \le i \le n$ the minimum point of $\{||\alpha^* - \Psi(\varphi(\theta + s_i\Delta\theta))||_p, 1 \le i \le n\}$
13: update $\theta = \theta + s_i\Delta\theta$
14: **if** $E(\theta) > \epsilon$ **then**
15: run the projection $\theta = \mathcal{P}(\theta)$
16: **end if**
17: compute $\alpha = \Psi(\varphi(\theta))$
18: **end while**
19: **return** the value of the last iteration of $\theta$

## 2.5. A first numerical illustration of the LVGF method

In this section, we give a first numerical illustration of the LVGF method by applying it to the identification of two scalar parameters in a simple ODE system. The fact of being in dimension 2 makes it easy to represent the successive iterations of the LVGF method and we present a comparison with the gradient descent method applied to the minimization of the classical functional (1.5).
The ODE system that we have considered is the Van der Pol model given by:

$$ 
\begin{cases}
\dot{x} = v, \\
\dot{v} = \mu(1 - x^2)v - x, \\
(x, v)(0) = (x_0, v_0).
\end{cases}
\qquad (2.15) $$

In this equation, $\mu > 0$ is a parameter which reflects the degree of nonlinearity of this system, $x_0$ is the initial position and $v_0$ is the initial velocity. We assume that we measure the state variable $x$ on the whole time interval and we are interested in the identification of $\theta = (v_0, \mu)$ while $x_0$ is assumed to be known and its value is given by $x_0 = 0.5$. We consider a measurement $u^*$ defined by $u^* = \varphi(\theta^*)$ with $\theta^* = (0.5, 0.5)$ and our objective is to identify this value $\theta^*$. In addition to $u^*$, we consider that we have access to a data set $U$ (which contains $N_d = 1200$ simulations) which has been generated in a preliminary step by varying the values of $\mu$, randomly drawn from a uniform distribution in $[0.05, 2.0]$ and setting $v_0 = 0.5$. So in this case, the intrinsic dimension of the submanifold $U$ is equal to 1. The latent dimension of the autoencoder trained on this data set is fixed at $p = 1$ (this dimension can also be rediscovered following the method presented in Section 3.2). We used AE1 (mentioned in Table 1) to approximate the data set $U$. There are 1000 samples used in the training set and 200 in the validation set. Moreover, the parameters given as input of Algorithm 1 have been fixed at $N_{max} = 2000$, $\sigma = 10^{-4}$ and $S = [0.0005, 0.0025, 0.0125, 0.0625]$.

In addition to the LVGF method, we have implemented the gradient descent method applied to the minimization of the classical functional (1.5) in order to compare both methods. In Figure 1, we plot the contour of the classical cost function given by equation (1.5). We observe that the function is non-convex and presents local minima in addition to the global minimum at $\theta^* = (0.5, 0.5)$ represented by a red sign "+". The different iterations for the classical gradient descent method applied to equation (1.5) and for the LVGF method have also been depicted starting with the same initialization $\theta_0 = (0.8, 1.9)$.

It is interesting to observe that the gradient descent method converges as expected to a local minimizer whereas the LVGF method is able to go against the direction of the steepest slope in order to change valley and reach the global minimizer. Thus, the information coming from the data set and added through the use of the latent variable made it possible to correct the convexity defect of the classical functions and to identify the right value of the unknown parameters. Further numerical tests to better assess the performance of the LVGF method will be presented in Section 4.

### 3. Lipschitz-stable autoencoders and intrinsic dimension

As presented in Section 2.2, our LVGF method relies on a description of the data set $U = \{u^{(i)}\}_{1\le i\le N_d}$ by an autoencoder. These neural networks are dimensionality reduction algorithms that learn two functions, called encoder and decoder. The encoder function associates with an input vector in $\mathbb{R}^n$ a vector called a code or a latent variable of smaller dimension $p$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Contour plot of the classical function $J(\theta) = \|u^* - \varphi(\theta)\|_n^2$ with $\theta = (v_0, \mu)$: comparison between the iterations of the classical gradient descent method (blue dots) and the LVGF steps (red dots) with the same initial guess (blue cross). The value of the true parameter, which is the minimum of $J$ is $\theta^* = (0.5, 0.5)$, marked with a red cross.</figcaption></figure>

It is important to note that the dimension $p$ is a hyperparameter that has to be set before training the autoencoder, although, in general, the intrinsic dimension of a data set is unknown. That is why the first question that naturally arises in the representation of a data set by an autoencoder is to understand how to settle the dimension of the latent variable. In this section, we propose a criterion based on the notion of stable manifold width to estimate the intrinsic dimension of the data set. This criterion is a consequence of the result given by Corollary 3.2 which highlights the fact that, if the dimension of the latent variable is smaller than the intrinsic dimension, the product of the Lipschitz constants of the encoder and decoder functions will blow up. In Section 3.4, this criteria is applied through numerical tests.

### 3.1. Presentation of autoencoders

In this section, we describe the autoencoder neural network and give some details on the implementation of the autoencoder that we have considered in this work.

As described in [18] and [44], an autoencoder is a feed-forward neural network with a systemic structure that is composed of multiple hidden layers which include an encoding and a decoding part. The input and output layers have the same size as shown in Figure 2. The
<figure><img src="image_4.png" /><figcaption>Figure 2. autoencoder systemic structure.</figcaption></figure>

encoder compresses the input data ($u \in \mathbb{R}^n$) to some latent representation ($\alpha \in \mathbb{R}^p$):

$$ \Psi: \begin{cases} U \to \mathbb{R}^p \\ u \mapsto \alpha := \Psi(u) \end{cases} $$

while the decoder decodes the compressed latent representation $\alpha$ and reconstructs the input in $\tilde{u} \in \mathbb{R}^n$:

$$ \Psi^{\dagger}: \begin{cases} \mathbb{R}^p \to \mathbb{R}^n \\ \alpha \mapsto \tilde{u} := \Psi^{\dagger}(\alpha). \end{cases} $$

Among the key aspects that have to be set in order to define an autoencoder, we need to specify the loss function to be minimized. In our paper, we have considered the following loss function:

$$ \frac{1}{N_d} \sum_{i=1}^{N_d} \|u^{(i)} - \Psi^{\dagger}(\Psi(u^{(i)}))\|_n^2 \quad (3.1) $$

where $U = \{u^{(i)}\}_{1 \le i \le N_d}$ corresponds to the data set. In addition, we chose the exponential linear unit (ELU) [14] as the activation function:

$$ ELU(x) = \begin{cases} x & \text{if } x \ge 0 \\ \rho(\exp(x) - 1) & \text{otherwise} \end{cases} \quad (3.2) $$

where $\rho > 0$ is a constant number (we have taken $\rho = 1.0$ in our numerical tests). This activation function is of class $C^1$, a property which will allow us to estimate the Lipschitz constants associated with the encoder and decoder functions. The weights and biases of the hidden layers were optimized by using an Adam optimizer [56].

In the numerical tests described in Sections 3.4 and 4, different architectures will be considered and, for each of them, we will specify the values of the hyperparameters (numbers of hidden layers, number of neurons per layer).
## 3.2. A criterion based on stable manifold widths to determine the latent variable dimension

In this section, we consider $\mathcal{U} \subset \mathbb{R}^n$ a compact embedded submanifold of topological dimension $\bar{p}$ and we propose a criterion based on autoencoders to identify the dimension $\bar{p}$. We refer to [12] and [22] for an overview of the intrinsic dimension estimation methods, the methods aim to project the original data $\mathcal{U} \subset \mathbb{R}^n$ to a lower M-dimensional submanifold of $\mathbb{R}^n$ that $M < n$ in a way that we will not lose the information of the original data. When the minimum M-dimension is necessary to represent the observed properties in $\mathcal{U}$, $M$ can be called intrinsic dimension.

Having in mind autoencoder applications, we are interested in nonlinear methods of approximation of $\mathcal{U}$ depending on $p$ parameters and built on two Lipschitz mappings that correspond to the encoder map and the decoder map. Following [15], we then introduce the quantity: for $p \in \mathbb{N}$ and for given constants $\gamma > 0$ and $\gamma^\dagger > 0$,

$$ \delta_{p,\gamma,\gamma^\dagger}(\mathcal{U}) = \inf_{\Psi, \Psi^\dagger} \sup_{u \in \mathcal{U}} \|u - \Psi^\dagger \circ \Psi(u)\|_n $$

where the infimum is taken over the $\gamma$ Lipschitz functions $\Psi : \mathcal{U} \to \mathbb{R}^p$ and the $\gamma^\dagger$ Lipschitz functions $\Psi^\dagger : \mathbb{R}^p \to \mathbb{R}^n$.

This definition is introduced in [15] under the name of stable manifold width with the slight nuance that the definition of [15] also involves the infimum over all the norms in $\mathbb{R}^p$. As explained in [15], the concept of stable manifold width compared to the concept of manifold width (where the functions $\Psi$ and $\Psi^\dagger$ are only assumed to be continuous) is motivated by the fact that the Lipschitz regularity allows to explicitly control perturbations coming from noise or numerical approximation.

In the following proposition, we are interested by the case where $p < \bar{p}$. This corresponds to the case where the dimension of the latent representation space is smaller than the topological dimension of the submanifold that the autoencoder has to approximate. Under this assumption, the proposition gives a bound from below for the stable manifold width.

**Proposition 3.1.** Let $\mathcal{U}$ be a compact embedded submanifold of $\mathbb{R}^n$ of topological dimension $\bar{p}$. We assume that $p < \bar{p}$. Then, there exists a constant $C > 0$ which only depends on $\mathcal{U}$ such that

$$ \delta_{p,\gamma,\gamma^\dagger}(\mathcal{U}) \geq \left( \frac{C}{\gamma^\dagger \gamma} \right)^{p/(\bar{p}-p)} R \quad (3.3) $$

where $R > 0$ is such that $\mathcal{U} \subset B_n(0, R)$.

For $m \in \mathbb{N}^*$, we have denoted by $B_m(y, r)$ with $y \in \mathbb{R}^m$ and $r > 0$ the closed ball of center $y$ and radius $r$.

*Proof.* The proof of this proposition will rely on a result given in [45]. This paper introduces and studies the notion of Lipschitz width, a concept close to stable manifold width. Let $\epsilon$ be such that $\epsilon > \delta_{p,\gamma,\gamma^\dagger}(\mathcal{U})$. Then, there exist a $\gamma$ Lipschitz function $\Psi : \mathcal{U} \to \mathbb{R}^p$ and a $\gamma^\dagger$ Lipschitz function $\Psi^\dagger : \mathbb{R}^p \to \mathbb{R}^n$ such that

$$ \sup_{u \in \mathcal{U}} \|u - \Psi^\dagger \circ \Psi(u)\|_n \leq \epsilon. $$
Since $\Psi$ is $\gamma$ Lipschitz and $\mathcal{U} \subset B_n(0, R)$, we have $\Psi(\mathcal{U}) \subset B_p(\Psi(0), \gamma R)$. Thus, we get that

$$ \inf_{\Psi_p^\dagger} \sup_{u \in \mathcal{U}} \inf_{\alpha \in B_p(\Psi(0), \gamma R)} \|u - \Psi_p^\dagger(\alpha)\|_n \le \epsilon $$

where the infimum is taken over all the maps $\Psi_p^\dagger : (B_p(\Psi(0), \gamma R), \| \cdot \|_p) \to \mathbb{R}^n$ such that

$$ \sup_{\alpha_1, \alpha_2 \in B_p(\Psi(0), \gamma R)} \frac{\|\Psi_p^\dagger(\alpha_1) - \Psi_p^\dagger(\alpha_2)\|_n}{\|\alpha_1 - \alpha_2\|_p} \le \gamma^\dagger $$

Let us now introduce the scaled norm $\| \cdot \|_{p,s}$ in $\mathbb{R}^p$ defined by $\|\alpha\|_{p,s} = \frac{1}{\gamma R} \|\alpha\|_p$. Then we have

$$ \inf_{\Psi_p^\dagger} \sup_{u \in \mathcal{U}} \inf_{\alpha \in B_{p,s}(\Psi(0),1)} \|u - \Psi_p^\dagger(\alpha)\|_n \le \epsilon $$

where the infimum is taken over all the maps $\Psi_p^\dagger : (B_{p,s}(\Psi(0),1), \| \cdot \|_{p,s}) \to \mathbb{R}^n$ such that

$$ \sup_{\alpha_1, \alpha_2 \in B_{p,s}(\Psi(0),1)} \frac{\|\Psi_p^\dagger(\alpha_1) - \Psi_p^\dagger(\alpha_2)\|_n}{\|\alpha_1 - \alpha_2\|_{p,s}} \le \gamma^\dagger \gamma R. $$

This property coincides with the fact that the fixed Lipschitz width of $\mathcal{U}$ associated with the norm $\| \cdot \|_{p,s}$ for the Lipschitz constant $\gamma^\dagger \gamma R$ is smaller than $\epsilon$. Now, let us introduce the Lipschitz width of $\mathcal{U}$ for the Lipschitz constant $\gamma^\dagger \gamma R$ which is defined as the infimum of the fixed Lipschitz widths over all the norms in $\mathbb{R}^p$ and denote it by $d_{p,\gamma^\dagger \gamma R}(\mathcal{U})$. We thus have the following property:

$$ d_{p,\gamma^\dagger \gamma R}(\mathcal{U}) < \epsilon. $$

Thus according to Proposition 3.5 in [45], this implies that we have an upper bound on the Lipschitz constant which is given by:

$$ \gamma^\dagger \gamma R \ge \frac{1}{3} \epsilon N_{2\epsilon}^{1/p}(\mathcal{U}) $$

where $N_\epsilon(\mathcal{U})$ is the $\epsilon$-covering number of $\mathcal{U}$. Since $\mathcal{U} \subset B_n(0, R)$ is a compact embedded submanifold of dimension $\bar{p}$, we have that

$$ N_\epsilon(\mathcal{U}) \ge C \frac{R^{\bar{p}}}{\epsilon^{\bar{p}}}, $$

where $C$ depends on the Lipschitz constants of the local maps of the finite atlas describing $\mathcal{U}$. Thus, we get that

$$ \gamma^\dagger \gamma \ge C R^{\bar{p}/p-1} \epsilon^{1-\bar{p}/p} $$

which implies that

$$ \epsilon \ge \left( \frac{C}{\gamma^\dagger \gamma} \right)^{p/(\bar{p}-p)} R $$

for every $\epsilon > \delta_{p,\gamma,\gamma^\dagger}(\mathcal{U})$. This property allows us to conclude the proof. $\square$
In the following corollary, we give another formulation of the previous proposition which will be more useful for our numerical tests.

**Corollary 3.2.** Under the same notations and hypotheses as in Proposition 3.1, we assume that there exist a $\gamma$ Lipschitz function $\Psi: \mathcal{U} \to \mathbb{R}^p$ and a $\gamma^\dagger$ Lipschitz function $\Psi^\dagger: \mathbb{R}^p \to \mathbb{R}^n$ such that

$$ \sup_{u \in \mathcal{U}} \|u - \Psi^\dagger \circ \Psi(u)\|_n \le \epsilon. $$

Then, there exists a constant $C > 0$ which only depends on $\mathcal{U}$ such that

$$ \gamma^\dagger \gamma \ge C R^{\bar{p}/p-1} \epsilon^{1-\bar{p}/p}. \quad (3.4) $$

This result can be interpreted in the following way: if the dimension $p$ is smaller than the topological dimension of $\mathcal{U}$, then getting an accurate approximation of the elements of $\mathcal{U}$ with the latent dimension $p$ may be achieved only if the Lipschitz constants of the approximation mappings are sufficiently large. It is related to the fact that, as long as the latent dimension is smaller than the intrinsic dimension, the encoder and decoder functions have to compensate this with strong variations, in the same vein as the space-filling curves. Our criterion thus relies on the representation of the variation of the product of the Lipschitz constants of the encoder and decoder functions with respect to the latent dimension $p$ at a given accuracy $\epsilon$. Its relevance will be illustrated in Section 3.4.

### 3.3. Estimation of the Lipschitz constants

The criterion that we propose to estimate the intrinsic dimension relies on the inequality (3.4). Therefore, to be able to implement it numerically, it is necessary to propose a numerical method to approximate the Lipschitz constant $\gamma^\dagger$ of the encoder $\Psi$ and the Lipschitz constant $\gamma$ of the decoder $\Psi^\dagger$. In this section, we explain how this is numerically achieved.

Since the activation function used to build these functions is the ELU, $\Psi$ and $\Psi^\dagger$ are regular functions, in particular, they are $C^1$ functions. Let us introduce the gradient of the function $\Psi: \mathcal{U} \to \mathbb{R}^p$ and denote it by $\nabla \Psi$:

$$ \nabla \Psi: \begin{cases} \mathcal{U} \to \mathbb{R}^{p \times n} \\ u \mapsto \nabla \Psi(u) \end{cases} $$

Similarly, we denote by $\nabla \Psi^\dagger$ the gradient of $\Psi^\dagger: \mathbb{R}^p \to \mathbb{R}^n$:

$$ \nabla \Psi^\dagger: \begin{cases} \mathbb{R}^p \to \mathbb{R}^{n \times p} \\ \alpha \mapsto \nabla \Psi^\dagger(\alpha) \end{cases} $$

The estimation of the product $\gamma^\dagger \gamma$ of the Lipschitz constants relies on the following lemma whose proof uses classical arguments.
**Lemma 3.3.** For a function $f: \mathcal{U} \subset \mathbb{R}^n \to \mathbb{R}^n$ of class $C^1$, we have the following estimate: for all $u^{(1)}, u^{(2)} \in \mathcal{U}$,

$$ \frac{\|f(u^{(1)}) - f(u^{(2)})\|_{l^{q,n}}}{\|u^{(1)} - u^{(2)}\|_{l^{q,n}}} \leq \left[ \sum_{k=1}^{n} \left( \sum_{j=1}^{n} \sup_{u \in \mathcal{U}} \left| \frac{\partial f_k}{\partial u_j}(u) \right|^z \right)^{\frac{q}{z}} \right]^{\frac{1}{q}} $$

where $z, q \in ]1, +\infty[$ are such that $\frac{1}{z} + \frac{1}{q} = 1$.

In particular, if we set $z = q = 2$, we get that, for all $u^{(1)}, u^{(2)} \in \mathcal{U}$,

$$ \frac{\|f(u^{(1)}) - f(u^{(2)})\|_n}{\|u^{(1)} - u^{(2)}\|_n} \leq \left( \sum_{j,k=1}^{n} \sup_{u \in \mathcal{U}} \left| \frac{\partial f_k}{\partial u_j}(u) \right|^2 \right)^{\frac{1}{2}} $$

Thus, if we apply this estimate to $f = \Psi^\dagger \circ \Psi$, we can have an approximation of a lower bound of $\gamma^\dagger \gamma$ by computing

$$ \left( \sum_{j,k=1}^{n} \sum_{j=1}^{n} \max_{1 \leq m \leq N} \left| \sum_{l=1}^{p} \frac{\partial \Psi_k^\dagger}{\partial \alpha_l}(\Psi(u^{(m)})) \frac{\partial \Psi_l}{\partial u_j}(u^{(m)}) \right|^2 \right)^{\frac{1}{2}} $$

where we have selected $N$ samples $(u^{(m)})_{1 \leq m \leq N}$ to have a discrete approximation of the supremum value of the gradient.

## 3.4. Numerical tests of intrinsic dimension estimations

In this subsection, we will present some numerical tests on the estimation of the intrinsic dimension of the data set by using an autoencoder, as described in Section 3.2. In particular, we will compare the criterion that we propose which is based on the estimation of the Lipschitz constants with a more classical criterion (we refer to [33] for a presentation of other criteria) based solely on the evolution of the accuracy with respect to the dimension of the latent dimension.

For the tests which we will present hereafter, we consider the architectures of three autoencoders, named AE1, AE2, and AE3, and described in Table 1, in which only half of the architecture is written, since we consider symmetric autoencoders (the encoder and the decoder have the same architectures). The number of layers corresponds to the number of hidden layers.

Each autoencoder is trained for different choices of $p$, the dimension of the latent layer. Then, the product of the Lipschitz constants $\gamma^\dagger \gamma$ is evaluated following the method presented in Section 3.3 and the error $\epsilon$ (which corresponds to the value of the loss function given by (3.1)) is computed in a validation set.
<table><caption>**Table 1.** Architectures of the autoencoder models</caption><thead><tr><td></td><td>Number of layers</td><td>Number of hidden units</td></tr></thead><tbody><tr><td>AE1</td><td>7</td><td>150, 80, 20, p</td></tr><tr><td>AE2</td><td>11</td><td>250, 150, 80, 20, 10, p</td></tr><tr><td>AE3</td><td>15</td><td>250, 180, 120, 80, 50, 30, 15, p</td></tr></tbody></table>

### 3.4.1. Tests on Van der Pol model
To build a data set from Van der Pol model (2.15), among the three parameters of the model, we fixed the initial velocity $x_0 = 0.5$ while the parameters $v_0$ and $\mu$ were randomly drawn from a uniform distribution in [0.05, 1.5] and [0.05, 2.0] respectively. Then, for each parameter value, we computed an approximated solution of the ODE system in [0, T], with $T = 30$ thanks to the Crank-Nicolson scheme with the time step given by $\Delta t = 0.075$. A total of $N_d = 1200$ solutions were generated, of which 1000 were used in the training set and 200 in the validation set.

The training of the autoencoder was performed by optimizing the loss function defined in equation (3.1) by using the Adam optimiser and by taking 1000 epochs with a batch size equal to 40. Each training for a given $p$ was repeated 10 times and the results presented hereafter correspond to the average of these 10 tests.

In Figure 3(a), we have represented the variations of the error with respect to $p$ for different architectures. As explained in [33], these curves allow us to identify the intrinsic dimension which corresponds to the value of $p$ from which the curve begins to stagnate. In Figure 3(b), we have represented the variations of the products of the Lipschitz constants to apply the criterion that we proposed. For both methods, we can observe a stagnation for $p = 2$. Therefore, for this simple data set, both criteria are able to identify the right intrinsic dimension $\bar{p} = 2$.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Intrinsic dimension estimation for Van der Pol model and $\bar{p} = 2$</figcaption></figure>
### 3.4.2. Tests on the FitzHugh-Nagumo model
In this section, we present the results obtained on the intrinsic dimension estimation for the solutions of the FitzHugh-Nagumo model. This model which describes the dynamics of a spiking neuron and is given by

$$ 
\left\{ 
\begin{aligned}
\dot{v} &= v - \frac{v^3}{3} - w + I_{ext}, \\
\dot{w} &= \frac{1}{\tau}(v + a - bw), \\
(v, w)(0) &= (v_0, w_0)
\end{aligned}
\right. 
\qquad (3.5) 
$$

where $v$ corresponds to the membrane potential, $w$ to the recovery variable and $I_{ext}$ to a stimulus current. We assume that we only observe the state variable $v$ and that the value of $I_{ext}$ is known.

The numerical approximation of FitzHugh-Nagumo model is carried out by using the Crank-Nicolson scheme on a time interval $(0, T)$ with $T = 200$ and with the time step given by $\Delta t = 0.5$.

We consider here a scenario in which we test the dimension estimation when $\bar{p} = 5$. To generate the data set, we have kept $I_{ext}$ fixed at the value $I_{ext} = 0.325$ whereas $b, \tau, a, v_0$ and $w_0$ are drawn from a uniform distribution respectively in the intervals $[0.05, 0.5]$, $[12.0, 13.0]$, $[0.05, 0.5]$, $[0.05, 1.0]$, and $[0.05, 1.0]$.

If we consider the curve of the accuracy with respect to the latent dimension $p$, reported in Figure 4(a), we observe a stagnation from $p = 3$, which would lead to a wrong conclusion. Let us mention that the difficulty to identify the latent dimension thanks to this simple and natural criterion has already been highlighted in several papers (we refer for instance to [33]). In particular, in relatively high dimensions, it is often observed that the curve gradually decreases and starts to stagnate before reaching the intrinsic dimension. On the other hand, if we consider the curve of the product of the Lipschitz constants with respect to the latent dimension, reported in Figure 4(b), we can observe that a stagnation only occurs from $p = 5$. Thus the criterion that we propose is able to correctly identify the value of the intrinsic dimension.

## 4. Numerical results for the parameter identification with the LVGF method

In this section, we will assess the performances of the LVGF method for the identification of parameters in different ODE models. For each numerical test, we will compare the results given by the LVGF method with what we will call the classical method: it consists in minimizing the classical misfit functional $J$ given by (1.5) thanks to a gradient descent algorithm with a fixed step size.

Before describing each test, let us start with some practical details for the implementation of the LVGF method and the classical method. For the training process in the preliminary step of the LVGF method, we have trained the autoencoders in 1000 epochs with 40 batch sizes. In Algorithm 1 which describes the LVGF method, the maximum number of iterations has been set to $N_{max} = 2000$, the tolerance for $\|\alpha^* - \alpha_k\|_p$ has been set to $\sigma = 10^{-4}$ and the steps were given by $S = [0.0005, 0.0025, 0.0125, 0.0625]$.
<figure><img src="image_6.png" /><figcaption>Figure 4. Intrinsic dimension estimation for FitzHugh-Nagumo model and $\bar{p} = 5$</figcaption></figure>

To quantify the errors in the parameter estimation, we evaluated the total error by computing the root relative sum of squares (RRSSE):

$$RRSSE(\theta^*, \theta) = \frac{\|\theta^* - \theta\|_m}{\|\theta^*\|_m}$$

and, to evaluate the errors in each of the parameters, we also computed the following relative error (RE) for each parameter:

$$RE(\theta_i^*, \theta_i) = \frac{|\theta_i^* - \theta_i|}{|\theta_i^*|}, \quad 1 \le i \le m.$$

At last, the error between the observed signal $u^*$ and the signal associated with the estimated parameters $u = \varphi(\theta)$ is defined as the root relative mean squared error (RRMSE):

$$RRMSE(u^*, u) = \frac{\|u^* - u\|_n}{\|u^*\|_n}.$$

### 4.1. Parameter identification for Van der Pol model

In this section, we consider Van der Pol model and are interested by the identification of the parameter vector $\theta = (x_0, v_0, \mu)$.

For the LVGF method, we have taken the same data set of 1200 samples and the same autoencoder setup as in Section 3.4.1. In particular, in the data set, $\theta_1$ has a valued fixed at $\theta_1 = 0.5$, whereas $\theta_2$ and $\theta_3$ are, respectively, drawn from uniform distributions in [0.05, 1.5] and in [0.05, 2.0]. The architecture of the autoencoder corresponds to AE1 described in Table 1 and its latent dimension has been identified in Section 3.4.1 to be equal to $p = 2$ thanks to the criterion presented in Section 3.2.

Before showing numerical tests on more complex models for which we will present an in-depth and statistical comparison between the classical method and the LVGF method, we
illustrate the results with a case which is representative of what is observed in general for the Van der Pol model. In this example, we consider the case where the parameter to recover is given by $\theta^* = (0.5, 0.5, 0.5)$ and where we start from the initial guess $\theta = (0.6615, 1.23, 1.6994)$. The results obtained in this specific case are reported in Table 2. We can see that the values obtained with the LVGF method are significantly closer to the right values than the classical method. To illustrate the convergence issues of the classical method, we have also represented the parameter values at each iteration (to better visualize them, we have chosen to represent only $\theta_2$ and $\theta_3$) in Figure 5 (a) and the values taken by the misfit function $J$ given by (1.5) in Figure 5 (b). We note that the gradient method experiences convergence issues and that the iterations stagnate to a local minimum which is far away from $\theta^*$.

<table><caption>**Table 2.** Comparison between the classical method and the LVGF method in a specific case for Van der Pol model: parameter to identify $\theta^* = (0.5, 0.5, 0.5)$ and initial guess $\theta = (0.6615, 1.23, 1.6994)$</caption><thead><tr><th></th><th>estimated θ</th><th>RRSSE of the estimated θ</th><th>RRMSE of the estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(0.44, -0.16, 2.32)</td><td>3.8666</td><td>1.81</td></tr><tr><td>LVGF method</td><td>(0.52, 0.52, 0.502)</td><td>0.0541</td><td>0.011</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Illustration of the test presented in Table 2: representations of the iterations with the classical method</figcaption></figure>

On the other hand, with the LVGF method, the iterations converge to the correct values of the parameters, as illustrated in Figure 6. In addition, Figure 7 also depicts the iterations of the latent variable $\alpha = (\alpha_1, \alpha_2)$ and we observe that the descent direction for this variable allows us to go in a relatively direct way towards the value $\alpha^*$. Therefore, in the latent variable space, we observe that the evolution of the iterations $(\alpha_k)_{k \in \mathbb{N}}$ of the LVGF method can be assimilated to the evolution of iterations of the classical gradient descent method for the functional $\|\alpha - \alpha^*\|_p^2$. This observation perfectly conforms with the principles of the LVGF method presented in Section 2.1.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Illustration of the test presented in Table 2: representations of the iterations with the LVGF method</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Illustration of the test presented in Table 2: values of $\alpha$ at each iteration</figcaption></figure>

To conclude this first series of numerical results for the Van der Pol model, let us compare the computational costs between the LVGF method and the classical gradient descent method. A first comment is that, compared to the classical method, a pre-processing step is necessary for the LVGF method, in order to train the autoencoder on the data set. This point aside, the computational costs of the LVGF method are comparable to the classical gradient descent method. More precisely, an iteration of LVGF lasts between two and three times longer with respect to an iteration of the classical gradient descent method. Since the stopping criteria for both methods are different (for LVGF, the method stops when $||\alpha - \alpha^*||$ is smaller than a threshold value whereas, for the gradient descent method, the method stops when the norm of the gradient of the misfit functional is smaller than a threshold value), it is difficult to compare the number of iterations. However, we observe that LVGF needs about two times less iterations to reach convergence. For instance, for a threshold value of $10^{-4}$, we need about 1000 iterations compared with 2000 iterations for the classical method. Moreover, a key property of the LVGF method is that, unlike the classical method, it does not encounter convergence
issues.

## 4.2. Parameter identification for FitzHugh-Nagumo model

Let us now present the numerical results obtained for the identification of parameters in FitzHugh-Nagumo model (3.5). In these tests, the values of $I_{ext}$ and of the initial condition of $v$ are fixed at $I_{ext} = 0.325$ and $v_0 = 1$ and we want to identify the four remaining parameters $\theta = (w_0, a, b, \tau)$.

For the LVGF method, we have considered a data set of 1200 samples corresponding to a fixed value of $\tau$ given by $\tau = 12.5$ whereas the values of $w_0, a$ and $b$ are respectively drawn from a uniform distribution in $[0.05, 1.0], [0.05, 0.8], \text{ and } [0.05, 0.8]$. With regard to the setup and training process of the autoencoder, we made the same choices as in the previous section for Van der Pol model, except that the dimension of the latent layer is given by $p = 3$.

To test the performances of the LVGF method and compare them with the classical method (which is presented in the beginning of Section 4), we made statistical tests by considering 10 random values of $\theta^*$ and 20 initial guesses for each value of $\theta^*$. We took the parameters in the same range as the data set for $w_0, a, b$ and in the range of $[12, 13]$ for $\tau$.

The results are presented in three tables: Table 3 gives the averaged results for noiseless data whereas Table 4 details the results obtained for two values of $\theta^*$ (corresponding to the best case and the worst case) still averaged on the different initial guesses. In addition, Table 5 corresponds to the averaged results with different levels of noise in the data.

As for the Van der Pol model, we observe that the relative error made on the parameters is smaller with the LVGF method than with the classical method. More precisely, we see that the error is reduced by a factor of 2 or more and that the error on the signal is reduced by a factor of about 7 for the noiseless case. If we look at the relative error parameter by parameter, we also observe that the error on the reconstruction of the parameter $b$ (third parameter in $\theta$) is relatively large with both methods still with a clear improvement with the LVGF method (around 54% with the classical method and 27% for the LVGF method). Since this lack of accuracy has a rather low impact on the reconstruction of the signal, this difficulty to identify the parameter $b$ compared to $w_0, a$, and $\tau$ is related to the differences of sensitivities of the signal with respect to the parameters: its sensitivity with respect to $b$ is smaller than its sensitivity with respect to the other parameters and so a relatively rough identification still allows to accurately reconstruct the signal. For the noisy data, we can see that we lose precision in the parameters estimation and signals reconstruction but the LVGF method still reduces by a factor of about 2 the errors both on the signal and on the parameter compared with the classical method.

At last, the results for a specific example corresponding to $\theta^* = (0.5, 0.5, 0.5, 12.5)$ with the initial guess $\theta = (0.4815, 0.223, 0.689, 12.52)$ are detailed in Table 6. In that case, the reconstruction of the signal associated with the parameters is presented in Figure 8. We can see that the signal generated by the parameters obtained with the LVGF method perfectly fits the measured signal.
<table><caption>**Table 3.** Comparison between the classical method and the LVGF method: averaged results on 10 tests (FitzHugh-Nagumo model)</caption><thead><tr><th>Method</th><th>RE of each estimated θ</th><th>RRSSE of the estimated θ</th><th>RRMSE of the estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(0.174, 0.322, 0.544, 0.036)</td><td>0.703</td><td>0.3662</td></tr><tr><td>LVGF method</td><td>(0.074, 0.08, 0.27, 0.014)</td><td>0.3026</td><td>0.0529</td></tr></tbody></table>

<table><caption>**Table 4.** Comparison between the classical method and the LVGF method for FitzHugh-Nagumo model: detailed results for two specific cases (the best case and the worst case) and averaged results over 10 different parameter values. At each line, the results are averaged over the initial guesses.</caption><thead><tr><th></th><th>Methods</th><th>RE of each estimated θ</th><th>RRSSE of the estimated θ</th><th>RRMSE of the estimated signal</th></tr></thead><tbody><tr><td rowspan="2">Best case</td><td>Classical</td><td>(0.06, 0.14, 0.25, 0.02)</td><td>0.304</td><td>0.152</td></tr><tr><td>LVGF</td><td>(0.052, 0.035, 0.16, 0.014)</td><td>0.168</td><td>0.0142</td></tr><tr><td rowspan="2">Worst case</td><td>Classical</td><td>(0.33, 0.48, 1.74, 0.04)</td><td>1.867</td><td>0.579</td></tr><tr><td>LVGF</td><td>(0.1, 0.11, 0.47, 0.019)</td><td>0.506</td><td>0.076</td></tr></tbody></table>

<table><caption>**Table 5.** Comparison between the classical method and the LVGF method: averaged results on 10 tests for FitzHugh-Nagumo model with 0.33%, 1%, and 3% noise</caption><thead><tr><th>Noise level</th><th>Methods</th><th>RE of each estimated θ</th><th>RRSSE of the estimated θ</th><th>RRMSE of the estimated signal</th></tr></thead><tbody><tr><td rowspan="2">0.33%</td><td>classical method</td><td>(0.28, 0.27, 0.56, 0.028)</td><td>0.7604</td><td>0.3711</td></tr><tr><td>LVGF method</td><td>(0.099, 0.13, 0.29, 0.016)</td><td>0.3768</td><td>0.1884</td></tr><tr><td rowspan="2">1%</td><td>classical method</td><td>(0.21, 0.27, 0.6, 0.025)</td><td>0.735</td><td>0.341</td></tr><tr><td>LVGF method</td><td>(0.14, 0.11, 0.33, 0.017)</td><td>0.4091</td><td>0.185</td></tr><tr><td rowspan="2">3%</td><td>classical method</td><td>(0.21, 0.25, 0.57, 0.025)</td><td>0.6988</td><td>0.3337</td></tr><tr><td>LVGF method</td><td>(0.094, 0.096, 0.34, 0.018)</td><td>0.3805</td><td>0.1393</td></tr></tbody></table>

<table><caption>**Table 6.** Comparison between the classical method and the LVGF method in a specific case for FitzHugh-Nagumo model: parameter to identify $\theta^* = (0.5, 0.5, 0.5, 12.5)$ and initial guess $\theta = (0.4815, 0.223, 0.689, 12.52)$</caption><thead><tr><th></th><th>Estimated θ</th><th>RRSSE of estimated θ</th><th>RRMSE of estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(0.59, 0.05, 0.59, 12.99)</td><td>0.9344</td><td>0.6695</td></tr><tr><td>LVGF method</td><td>(0.503, 0.502, 0.508, 12.48)</td><td>0.0203</td><td>0.00624</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>(a) Signal generated with the parameters obtained (b) Signal generated with the parameters obtained with the classical method with the LVGF method</figcaption></figure>

<caption>**Figure 8.** Illustration of the tests presented in Table 6: comparison of the signals with the classical method and with the LVGF method (FitzHugh-Nagumo model)</caption>

In terms of comparison of computational costs between the LVGF method and the classical gradient descent method, we have observations similar to those with the Van der Pol model. More precisely, an iteration of the LVGF method is about three times longer than an iteration of the classical method.

### 4.3. Parameter identification for the Hindmarsh-Rose model

Let us now present the results that we have obtained for a more complex ODE system. Hindmarsh-Rose model is an ODE system that simulates the spiking-bursting behavior of the membrane potential of the neuron [32]. The system of ordinary differential equations is the following:

$$ 
\left\{ 
\begin{aligned}
\dot{x} &= y - ax^3 + bx^2 - z + I \\
\dot{y} &= c - dx^2 - y \\
\dot{z} &= r(s(x - x_r) - z) \\
(x, y, z)(0) &= (x_0, y_0, z_0)
\end{aligned}
\right. 
\qquad (4.1) $$

For this model, we assume that the parameters $I, r, a, y_0$, and $z_0$ are fixed and given by $I = 1.6, r = 0.005, a = 1, y_0 = 1$, and $z_0 = 1$ and our objective is to identify the values of the 6 parameters $x_0, b, c, d, s$, and $x_r$.

We assume that we have a data set of 1200 samples that has been generated by considering that the parameter $x_r$ is fixed at $x_r = -1.6$ and by taking values of $x_0, b, c, d$, and $s$, respectively drawn from a uniform distribution in [0.5, 1.5], [2.56, 2.76], [0.9, 1.1], [4.5, 5.0], and [3.5, 4.0]. The procedure for training of the autoencoder and the LVGF method is similar to the one implemented for Van der Pol and Fitzhugh-Nagumo models, except that the latent dimension is now $p = 5$.

Similarly as in the previous subsection, we have performed statistical tests by randomly generating 5 different $\theta^*$ and considering 10 initial guesses for each $\theta^*$. The results are summarized in Tables 7 and 8. For noiseless data, the results highlight that the LVGF method
<table><caption>**Table 7.** Comparison between the classical method and the LVGF method: averaged results on 5 tests (Hindmarsh-Rose model)</caption><thead><tr><th>Method</th><th>RE of each<br>estimated θ</th><th>RRSSE of the<br>estimated θ</th><th>RRMSE of the<br>estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(0.252, 0.039, 0.067, 0.036, 0.044, 0.122)</td><td>0.319</td><td>0.895</td></tr><tr><td>LVGF method</td><td>(0.054, 0.0033, 0.053, 0.024, 0.037, 0.037)</td><td>0.108</td><td>0.059</td></tr></tbody></table>

<table><caption>**Table 8.** Comparison between the classical method and the LVGF method: averaged results on 5 tests for Hindmarsh-Rose model with 0.33%, 1%, and 3% noise</caption><thead><tr><th>Noise level</th><th>Methods</th><th>RE of each<br>estimated θ</th><th>RRSSE of the<br>estimated θ</th><th>RRMSE of the<br>estimated signal</th></tr></thead><tbody><tr><td rowspan="2">0.33%</td><td>classical method</td><td>(0.196, 0.035, 0.058, 0.039, 0.04, 0.12)</td><td>0.272</td><td>0.37</td></tr><tr><td>LVGF method</td><td>(0.06, 0.0055, 0.035, 0.018, 0.037, 0.036)</td><td>0.1104</td><td>0.0892</td></tr><tr><td rowspan="2">1%</td><td>classical method</td><td>(0.192, 0.035, 0.06, 0.038, 0.041, 0.115)</td><td>0.2638</td><td>0.8238</td></tr><tr><td>LVGF method</td><td>(0.078, 0.01, 0.044, 0.017, 0.036, 0.052)</td><td>0.1367</td><td>0.109</td></tr><tr><td rowspan="2">3%</td><td>classical method</td><td>(0.192, 0.034, 0.061, 0.038, 0.04, 0.114)</td><td>0.267</td><td>0.922</td></tr><tr><td>LVGF method</td><td>(0.124, 0.009, 0.046, 0.018, 0.038, 0.037)</td><td>0.1667</td><td>0.141</td></tr></tbody></table>

can provide a parameter estimation which has 3 times less error compared to the classical method. In addition, the reconstructed signal is 10 times closer to the true signal than with the classical descent method. When we have noisy data, the LVGF method is still able to reduce 50% of the error in the parameter estimation and gives 5 times less error in the signal reconstruction.

For a specific case study, the results are presented in Table 9 where we have considered the case $\theta^* = (1.0, 2.6, 1.0, 4.7, 3.8, -1.6)$ with an initial guess given by $\theta = (0.9, 2.8, 0.9, 4.8, 3.7, -1.8)$. For this specific case, we see that the LVGF method provides an estimation of the parameter with an error nearly 5 times lower than the classical method and that the reconstructed signals given by Figure 9 are also much closer to the true one.

<table><caption>**Table 9.** Comparison between the classical method and the LVGF method in a specific case for Hindmarsh-Rose model: parameter to identify $\theta^* = (1.0, 2.6, 1.0, 4.7, 3.8, -1.6)$ and initial guess $\theta = (0.9, 2.8, 0.9, 4.8, 3.7, -1.8)$</caption><thead><tr><th></th><th>Estimated θ</th><th>RRSSE of<br>estimated θ</th><th>RRMSE of<br>estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(1.45, 2.59, 0.99, 4.58, 3.53, -1.34)</td><td>0.1201</td><td>1.149</td></tr><tr><td>LVGF method</td><td>(0.95, 2.6, 1.01, 4.74, 3.65, -1.66)</td><td>0.024</td><td>0.0117</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>(a) Signal generated with the parameters obtained (b) Signal generated with the parameters obtained with the classical method with the LVGF method</figcaption></figure>

<caption>**Figure 9.** Illustration of the test presented in Table 7: comparison of the signals with the classical method and with the LVGF method (Hindmarsh-Rose model)</caption>

### 4.4. Parameter identification for an influenza model

To complete this series of numerical tests, we present the results obtained for an influenza model ([3, 11]) given by:

$$ 
\left\{ 
\begin{array}{l}
S' = -\beta S[I + \delta A] \\
L' = \beta S[I + \delta A] - \kappa L \\
I' = p\kappa L - \alpha I \\
A' = (1 - p)\kappa L - \eta A \\
(S, L, I, A)(0) = (S_0, L_0, I_0, A_0)
\end{array}
\right.
\qquad (4.2) $$

In this model, $S$ represents the number of individuals who are susceptible to contract the disease, $L$ represents the number of individuals who are latent, i.e. infected but not yet symptomatic, $I$ represents the number of infected and symptomatic individuals, and $A$ represents the number of asymptomatic but infected individuals. We assume that the initial conditions $S_0, A_0$, and $L_0$ are fixed and given by $S_0 = 1000, L_0 = 20$ and $A_0 = 5$.

All the parameters appearing in this model correspond to positive constants and, in addition, $p$ satisfies $0 < p < 1$. We have considered that the value of $\beta$ is known and have set its value at $\beta = 4 \times 10^{-4}$. Then, we try to estimate the 6 following scalar parameters $\theta = [\delta, \kappa, p, \alpha, \eta, I_0]$.

To generate a data set of 1200 samples, we have taken into account the following relationship between $p, \alpha, \eta$ and $\delta$ which is given in [3]:

$$ 0.4 \left( \frac{p}{\alpha} + \frac{\delta(1-p)}{\eta} \right) = 1.37. $$

In addition, we have considered that $\kappa, \alpha$, and $\eta$ have fixed values given by: $\kappa = 3, \alpha = 5$, and $\eta = 6$. Thus, it follows that $\delta$ can be expressed with respect to $p$: $\delta = \frac{20.55-1.2p}{1-p}$. Finally, we have taken the values of $p$ and of the integer $I_0$, respectively, drawn from a uniform distribution in
[0.0, 0.5] and [5, 15]. The latent dimension associated with this data set is 2 and the procedure for training the autoencoder and the LVGF method is similar to that presented for Van der Pol and Fitzhugh-Nagumo models.

As for the previous models, we performed statistical tests by taking 5 different $\theta^*$ and by randomly generating 10 initial guesses for each $\theta^*$. The results obtained are summarized in Table 10 for noiseless data and in Table 11 for noisy data. In the noiseless case, we observe that the LVGF method gives smaller errors in the estimation of $\theta^*$ and errors three times smaller on the estimated signal compared to the classical method. When considering noisy data, we also observe that the results given by the LVGF method are better than those with the classical method, both in terms of the estimation of $\theta^*$ and of the estimated signal.

<table><caption>**Table 10.** Comparison between the classical method and the LVGF method: averaged results on 5 tests (Influenza model)</caption><thead><tr><th>Method</th><th>RE of each<br>estimated θ</th><th>RRSSE of the<br>estimated θ</th><th>RRMSE of the<br>estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(0.1655, 0.0867, 0.171, 0.1228, 0.1195, 0.088)</td><td>0.3477</td><td>0.0265</td></tr><tr><td>LVGF method</td><td>(0.1344, 0.067, 0.1169, 0.1054, 0.117, 0.02868)</td><td>0.27</td><td>0.0089</td></tr></tbody></table>

<table><caption>**Table 11.** Comparison between the classical method and the LVGF method: averaged results on 5 tests for Influenza model with 0.33%, 1%, and 3% noise</caption><thead><tr><th>Noise level</th><th>Methods</th><th>RE of each<br>estimated θ</th><th>RRSSE of the<br>estimated θ</th><th>RRMSE of the<br>estimated signal</th></tr></thead><tbody><tr><td rowspan="2">0.33%</td><td>classical method</td><td>(0.1665, 0.0889, 0.1655, 0.122, 0.119, 0.143)</td><td>0.372</td><td>0.038</td></tr><tr><td>LVGF method</td><td>(0.143, 0.072, 0.129, 0.125, 0.1244, 0.035)</td><td>0.2762</td><td>0.0147</td></tr><tr><td rowspan="2">1%</td><td>classical method</td><td>(0.168, 0.0896, 0.166, 0.125, 0.126, 0.156)</td><td>0.381</td><td>0.042</td></tr><tr><td>LVGF method</td><td>(0.143, 0.0668, 0.151, 0.1301, 0.1138, 0.0405)</td><td>0.3052</td><td>0.0187</td></tr><tr><td rowspan="2">3%</td><td>classical method</td><td>(0.169, 0.09, 0.152, 0.121, 0.1188, 0.169)</td><td>0.39</td><td>0.052</td></tr><tr><td>LVGF method</td><td>(0.138, 0.064, 0.148, 0.15, 0.123, 0.057)</td><td>0.312</td><td>0.035</td></tr></tbody></table>

To complete this study, we also present in Table 12 the results obtained for a specific case. In this case, we have considered $\theta^* = (2.51, 3, 0.092, 5, 6, 11)$ with an initial guess given by $\theta = (25.74, 3.1, 0.21, 5.52, 5.06, 10)$. In this specific case, we see that the LVGF method can estimate $\theta^*$ with an error nearly 6 times smaller than the classical method. In Figure 10, we observe that the signal associated with the parameter estimated by the LVGF method is also much closer to the true one.

## 5. Conclusion and discussion

In this study, we have presented an innovative method called Latent Variables Gradient Flow (LVGF) which leverages the available data and the underlying dynamics of the system in order to identify parameters of the ODE model. This approach entails a dual-phase process. In the
<table><caption>**Table 12.** Comparison between the classical method and LVGF method in a specific case for influenza model: parameter to identify $\theta^* = (22.51, 3, 0.092, 5, 6, 11)$ and initial guess $\theta = (25.74, 3.1, 0.21, 5.52, 5.06, 10)$</caption><thead><tr><th></th><th>Estimated θ</th><th>RRSSE of estimated θ</th><th>RRMSE of estimated signal</th></tr></thead><tbody><tr><td>classical method</td><td>(25.096, 2.72, 0.19, 5.92, 6.98, 11.15)</td><td>1.074</td><td>0.057</td></tr><tr><td>LVGF method</td><td>(25.14, 2.85, 0.093, 5.001, 6.72, 11.01)</td><td>0.1743</td><td>0.00185</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>(a) Signal ($I$) generated with the parameters obtained with the classical method
(b) Signal ($I$) generated with the parameters obtained with the LVGF method</figcaption></figure>

<caption>**Figure 10.** Illustration of the test presented in Table 12: comparison of the signals with the classical method and with the LVGF method (Influenza model)</caption>

first phase, an autoencoder is trained on the data set in order to represent it in a compressed way. The description of the data set by a latent variable is then exploited in the second phase which corresponds to an optimization method that can be described as a gradient flow for the latent variables.

For this new method, we presented numerical tests and a comparison with the classical method corresponding to a gradient descent method applied to the minimisation of the classical data misfit functional (1.5). A noticeable property of the LVGF method is that it is able to avoid convergence towards local minimum points. In this way, the LVGF method generally gives much better results than the classical method for the parameters identification problem.

In addition, we have shown that the sequence defined by the LVGF method satisfies convergence properties. These results still need to be extended to the algorithm actually implemented where there is an additional projection step. Nevertheless, our study, even if it is incomplete, allows us to understand how the use of the non-linear mapping between the parameter space and the latent variable space allows us to convexify our initial minimization problem.

In a complementary manner to this parameter identification method, we also proposed a new criterion to identify the intrinsic dimension of a sub-manifold thanks to autoencoders. This criterion is justified by a theoretical result based on the notion of stable manifold width.
It describes the behavior of the product of the Lipschitz constants of the encoder and decoder when the latent variable is too small. We have compared our criterion to the classical criterion based on the evolution of the error with respect to the latent dimension and observed that it leads to better results. Indeed, in the case of the classical criterion, the error gradually decreases when the latent dimension increases, making it difficult to pick a single value for the dimension. On the contrary, our criterion clearly highlights two phases in the evolution of the product of the Lipschitz constants and the intrinsic dimension corresponds to the elbow of the curve.

In addition to the theoretical study which must be in-depth, this work opens a certain number of perspectives. First of all, we intend to continue exploring its use in complex situations corresponding to larger dynamical systems that involve a large number of parameters. Moreover, in the same spirit of the methods presented in [49], we could also try to test the LVGF method on spatially heterogeneous signals modeled by Partial Differential Equations instead of Ordinary Differential Equations.

At last, the autoencoder used to approximate the available population data is based on neural networks. This is not a necessary choice, and other manifold learning methods could be used ([13]). For instance, Principal Component Analysis (PCA), kernel-PCA, Locally Linear Embedding, Isomap, Laplacian Eigenmaps, and Semidefinite Embedding are all methods that could be used in order to build an encoder-decoder pair.

## Acknowledgements

We would like to thank the support from INSPIRE European Training Network which receives funding from the EU Horizon 2020 Research and Innovation programme, under the Marie Skłodowska-Curie GA 858070 [30]. We also appreciated all the help from Sylvain Bernasconi and Christophe Bleunven at NOTOCORD, an Instem company.

## Data availability statement

No new data were created or analyzed in this study.
References

[1] J. Adler and O. Öktem. "Solving ill-posed inverse problems using iterative deep neural networks". In: *Inverse Problems* 33.12 (2017), pp. 124007, 24. ISSN: 0266-5611,1361-6420. DOI: 10.1088/1361-6420/aa9581. URL: https://doi.org/10.1088/1361-6420/aa9581.

[2] F. Anstett-Collin, L. Denis-Vidal, and G. Millérioux. "A priori identifiability: An overview on definitions and approaches". In: *Annual Reviews in Control* 50 (2020), pp. 139-149. ISSN: 1367-5788. DOI: https://doi.org/10.1016/j.arcontrol.2020.10.006. URL: https://www.sciencedirect.com/science/article/pii/S1367578820300729.

[3] J. Arino, F. Brauer, P. van den Driessche, J. Watmough, and J. Wu. "Simple models for containment of a pandemic". In: *Journal of the Royal Society Interface* 3.8 (2006), pp. 453-457.

[4] R. C. Aster, B. Borchers, and C. H. Thurber. *Parameter estimation and inverse problems*. Elsevier, 2018.

[5] N. Baden and J. Villadsen. "A family of collocation based methods for parameter estimation in differential equations". In: *The Chemical Engineering Journal* 23.1 (1982), pp. 1-13.

[6] Y. Bard. "Comparison of gradient methods for the solution of nonlinear parameter estimation problems". In: *SIAM Journal on Numerical Analysis* 7.1 (1970), pp. 157-186.

[7] J. V. Beck and K. J. Arnold. *Parameter estimation in engineering and science*. James Beck, 1977.

[8] R. Bellman and K. J. Åström. "On structural identifiability". In: *Mathematical biosciences* 7.3-4 (1970), pp. 329-339.

[9] A. Bjorck. *Numerical Methods for Least Squares Problems*. Vol. 51. SIAM, 1996.

[10] B. V. D. Bosch and L. Hellinckx. "A new method for the estimation of parameters in differential equations". In: *AIChE Journal* 20.2 (1974), pp. 250-255.

[11] F. Brauer, C. Castillo-Chavez, Z. Feng, et al. *Mathematical models in epidemiology*. Vol. 32. Springer, 2019.

[12] F. Camastra and A. Staiano. "Intrinsic dimension estimation: Advances and open problems". In: *Information Sciences* 328 (2016), pp. 26-41.

[13] L. Cayton. *Algorithms for manifold learning*. eScholarship, University of California, 2008.

[14] D.-A. Clevert, T. Unterthiner, and S. Hochreiter. "Fast and accurate deep network learning by exponential linear units (elus)". In: *arXiv preprint arXiv:1511.07289* (2015).

[15] A. Cohen, R. DeVore, G. Petrova, and P. Wojtaszczyk. "Optimal stable nonlinear approximation". In: *Foundations of Computational Mathematics* 22.3 (2022), pp. 607-648.
[16] D. Csercsik, G. Szederkényi, and K. M. Hangos. "Identifiability of a Hodgkin-Huxley type ion channel under voltage step measurement conditions". In: *IFAC Proceedings Volumes* 43.5 (2010). 9th IFAC Symposium on Dynamics and Control of Process Systems, pp. 332-337. ISSN: 1474-6670. DOI: https://doi.org/10.3182/20100705-3-BE-2011.00055. URL: https://www.sciencedirect.com/science/article/pii/S1474667016300568.

[17] N. Cunniffe, F. Hamelin, A. Iggidr, A. Rapaport, and G. Sallet. "Observability, Identifiability and Epidemiology-A primer". In: (2023).

[18] H. A. Dau, V. Ciesielski, and A. Song. "Anomaly detection using replicator neural networks trained on examples of one class". In: *Simulated Evolution and Learning: 10th International Conference, SEAL 2014, Dunedin, New Zealand, December 15-18, 2014. Proceedings 10*. Springer. 2014, pp. 311-322.

[19] J. E. Dennis Jr and J. J. Moré. "Quasi-Newton methods, motivation and theory". In: *SIAM review* 19.1 (1977), pp. 46-89.

[20] T. T. Dufera, Y. C. Seboka, C. Fresneda Portillo, et al. "Parameter Estimation for Dynamical Systems Using a Deep Neural Network". In: *Applied Computational Intelligence and Soft Computing* 2022 (2022).

[21] M. Duff, N. D. Campbell, and M. J. Ehrhardt. "Regularising inverse problems with generative machine learning models". In: *Journal of Mathematical Imaging and Vision* 66.1 (2024), pp. 37-56.

[22] K. Fukunaga. *Introduction to statistical pattern recognition*. Elsevier, 2013.

[23] H. Gabbard, C. Messenger, I. S. Heng, F. Tonolini, and R. Murray-Smith. "Bayesian parameter estimation using conditional variational autoencoders for gravitational-wave astronomy". In: *Nature Physics* 18.1 (2022), pp. 112-117.

[24] L. Gabrielli, S. Tomassetti, S. Squartini, and C. Zinato. "Introducing deep machine learning for parameter estimation in physical modelling". In: *Proceedings of the 20th international conference on digital audio effects*. 2017.

[25] P. E. Gill and W. Murray. "Quasi-Newton methods for unconstrained optimization". In: *IMA Journal of Applied Mathematics* 9.1 (1972), pp. 91-108.

[26] G. Gottwald, F. Li, Y. Marzouk, and S. Reich. "Stable generative modeling using diffusion maps". In: *arXiv preprint arXiv:2401.04372* (2024).

[27] G. A. Gottwald and S. Reich. "Combining machine learning and data assimilation to forecast dynamical systems from noisy partial observations". In: *Chaos: An Interdisciplinary Journal of Nonlinear Science* 31.10 (2021).

[28] G. A. Gottwald and S. Reich. "Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation". In: *Physica D: Nonlinear Phenomena* 423 (2021), p. 132911.
[29] V. Grimm, A. Heinlein, A. Klawonn, M. Lanser, and J. Weber. "Estimating the time-dependent contact rate of SIR and SEIR models in mathematical epidemiology using physics-informed neural networks". In: *Electronic Transactions on Numerical Analysis* 56 (2022), pp. 1-27.

[30] P.-J. Guns. ""INSPIRE"": A European training network in safety pharmacology creating opportunities for 15 PhD students". In: *Journal of Pharmacological and Toxicological Methods* 105 (2020), p. 106838.

[31] B. Hartke. "Global optimization". In: *Wiley Interdisciplinary Reviews: Computational Molecular Science* 1.6 (2011), pp. 879-887.

[32] J. L. Hindmarsh and R. Rose. "A model of neuronal bursting using three coupled first order differential equations". In: *Proceedings of the Royal society of London. Series B. Biological sciences* 221.1222 (1984), pp. 87-102.

[33] T. Kärkkäinen and J. Hänninen. "Additive autoencoder for dimension estimation". In: *Neurocomputing* 551 (2023), p. 126520.

[34] X. Kong, K. Yamashita, B. Foggo, and N. Yu. "Dynamic parameter estimation with physics-based neural ordinary differential equations". In: *2022 IEEE Power & Energy Society General Meeting (PESGM)*. IEEE. 2022, pp. 1-5.

[35] H. Liang and H. Wu. "Parameter Estimation for Differential Equation Models Using a Framework of Measurement Error in Regression Models". In: *Journal of the American Statistical Association* 103.484 (2008), pp. 1570-1583.

[36] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis. "Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators". In: *Nature machine intelligence* 3.3 (2021), pp. 218-229.

[37] Y. Marzouk, Z. Ren, S. Wang, and J. Zech. "Distribution learning via neural differential equations: a nonparametric statistical perspective". In: *arXiv preprint arXiv:2309.01043* (2023).

[38] M. T. McCann, K. H. Jin, and M. Unser. "Convolutional Neural Networks for Inverse Problems in Imaging: A Review". In: *IEEE Signal Processing Magazine* 34.6 (Nov. 2017), pp. 85-95. ISSN: 1053-5888. DOI: 10.1109/msp.2017.2739299. URL: http://dx.doi.org/10.1109/MSP.2017.2739299.

[39] H. Miao, X. Xia, A. S. Perelson, and H. Wu. "On identifiability of nonlinear ODE models and applications in viral dynamics". In: *SIAM review* 53.1 (2011), pp. 3-39.

[40] A. Mohammad-Djafari. "Regularization, Bayesian inference, and machine learning methods for inverse problems". In: *Entropy* 23.12 (2021), p. 1673.

[41] E. H. Moore. *On the reciprocal of the general algebraic matrix*. Vol. 26. 9. Bulletin of the American Mathematical Society, 1920, pp. 394-395.

[42] S. Mukherjee, S. Dittmer, Z. Shumaylov, S. Lunz, O. Öktem, and C.-B. Schönlieb. "Learned convex regularizers for inverse problems". In: *arXiv preprint arXiv:2008.02839* (2020).
[43] D. V. Patel, D. Ray, and A. A. Oberai. "Solution of physics-based Bayesian inverse problems with deep generative priors". In: *Computer Methods in Applied Mechanics and Engineering* 400 (2022), p. 115428.

[44] J. Pereira and M. Silveira. "Unsupervised representation learning and anomaly detection in ECG sequences". In: *International Journal of Data Mining and Bioinformatics* 22.4 (2019), pp. 389-407.

[45] G. Petrova and P. Wojtaszczyk. "Lipschitz widths". In: *Constructive Approximation* 57.2 (2023), pp. 759-805.

[46] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville. "On the spectral bias of neural networks". In: *International Conference on Machine Learning*. PMLR. 2019, pp. 5301-5310.

[47] J. R. Raol, G. Girija, and J. Singh. *Modelling and parameter estimation of dynamic systems*. Vol. 65. Iet, 2004.

[48] O. Sallam and M. Fürth. "On the use of Fourier Features-Physics Informed Neural Networks (FF-PINN) for forward and inverse fluid mechanics problems". In: *Proceedings of the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime Environment* 237 (4 2023), pp. 846-866.

[49] D. N. Tanyu, J. Ning, T. Freudenberg, N. Heilenkötter, A. Rademacher, U. Iben, and P. Maass. "Deep learning methods for partial differential equations and related parameter identification problems". In: *Inverse Problems* 39.10 (2023), p. 103001.

[50] R. Tenderini, S. Pagani, A. Quarteroni, and S. Deparis. "PDE-aware deep learning for inverse problems in cardiac electrophysiology". In: *SIAM Journal on Scientific Computing* 44.3 (2022), B605-B639.

[51] N. Verdière and C. Jauberthie. "Parameter estimation procedure based on input-output integro-differential polynomials. Application to the Hindmarsh-Rose model". In: *2020 European Control Conference (ECC)*. 2020, pp. 220-225. DOI: 10.23919/ECC51009.2020.9143670.

[52] A. F. Villaverde and G. Massonis. "On testing structural identifiability by a simple scaling method: Relying on scaling symmetries can be misleading". In: *PLoS computational biology* 17.10 (2021), e1009032.

[53] S. Wang, H. Wang, and P. Perdikaris. "Learning the solution operator of parametric partial differential equations with physics-informed DeepONets". In: *Science advances* 7.40 (2021), eabi8605.

[54] D. Wolf and R. Moros. "Estimating rate constants of heterogeneous catalytic reactions without supposition of rate determining surface steps—an application of a genetic algorithm". In: *Chemical Engineering Science* 52.7 (1997), pp. 1189-1199.
[55] A. Yazdani, L. Lu, M. Raissi, and G. E. Karniadakis. "Systems biology informed deep learning for inferring parameters and hidden dynamics". In: *PLOS Computational Biology* 16.11 (Nov. 2020), pp. 1–19. URL: https://doi.org/10.1371/journal.pcbi.1007575.

[56] A. Zhang, Z. C. Lipton, M. Li, and A. J. Smola. "Dive into deep learning". In: *arXiv preprint arXiv:2106.11342* (2021).

[57] A. Žilinskas and A. Zhigljavsky. "Stochastic global optimization: a review on the occasion of 25 years of Informatica". In: *Informatica* 27.2 (2016), pp. 229–256.