Sophie Carneiro Esteves

Antoine Vacavant
email: antoine.vacavant@uca.fr

Odyssée Merveille
email: odyssee.merveille@creatis.insa-lyon.fr

Learning a reconnecting regularization term for blood vessel variational segmentation

Keywords: blood vessel segmentation, reconnection, deep learning, variational approach

The segmentation of blood vessels in medical images is a challenging task as they are thin, connected and tortuous. The detection of a connected vascular network is of the utmost importance in clinical applications (e.g. blood flow simulations, vascular network modeling and analysis). Deep learning approaches have been developed to tackle this issue, but they require a large annotated dataset for each new application of interest, which is very challenging to build for vascular networks. In this work, rather than learning the segmentation task, we propose to learn a reconnecting regularization term that learns geometric properties of vascular networks independent of the image modality. Therefore, this term generalizes better than deep learning segmentation models, and can be easily plugged into variational segmentation frameworks to detect vascular networks in different datasets without requiring annotations. We apply this approach on retinal images by training our reconnecting term on the STARE dataset and applying it on the DRIVE dataset. We show that our approach better preserves the connectivity of vascular networks than classic regularization terms in the literature. Finally, we illustrate the generalization power of our reconnecting term by applying it to other types of data.

I. INTRODUCTION

The detection of blood vessels in medical images is the key task of most vascular network analysis tools that help clinicians diagnose and treat vascular diseases. Despite being studied for the past twenty years [START_REF] Moccia | Blood vessel segmentation algorithms-review of methods, datasets and evaluation metrics[END_REF], the segmentation of blood vessels is still an open problem as vascular networks are thin, connected, tortuous, sparse and low-contrasted structures.

In this work, we address the important problem of blood vessel connectivity loss upon segmentation of vascular networks. Preserving connectivity is key, as many subsequent processing tasks rely on a connected network such as blood flow simulation or vascular network modeling. Previous work have focused on this issue. Tracking strategies (see [START_REF] Moccia | Blood vessel segmentation algorithms-review of methods, datasets and evaluation metrics[END_REF] Section 7) intrinsically preserve the connectivity of blood vessels by recruiting neighbors pixels recursively, but usually require to be initialized with manual information. Other methods estimate the local direction of blood vessels in order to guide the segmentation in the vessel direction [START_REF] Merveille | nD variational restoration of curvilinear structures with prior-based directional regularization[END_REF][START_REF] Wu | Deep vessel tracking: A generalized probabilistic approach via deep learning[END_REF].

Recently, deep learning strategies have been proposed to segment vascular networks. Some of them aim at enforcing the network connectivity by using topological loss [START_REF] Haft-Javaherian | A topological encoding convolutional neural network for segmentation of 3D multiphoton images of brain vasculature using persistent homology[END_REF] or by computing both the segmentation and centerline of the vessels with the same network [START_REF] Tetteh | Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-d angiographic volumes[END_REF]. These strategies yield interesting results, but require a large annotated dataset to provide good results. By contrast, variational segmentation strategies do not require annotated data, but are less flexible to represent complex information such as connectivity. To combine the strengths of both methods, decoupled approaches [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF] were proposed. They consist in plugging an external operator to act as a regularization term in a variational framework. Early works used denoising algorithms such as BM3D or Non Local Means [START_REF] Heide | Flexisp: A flexible camera image processing framework[END_REF] to replace classic denoising regularization terms such as the Total Variation (TV). More recent strategies learn the external regularization term [START_REF] Kofler | Neural networks-based regularization for large-scale medical image reconstruction[END_REF][START_REF] Meinhardt | Learning proximal operators: Using denoising networks for regularizing inverse imaging problems[END_REF] to take advantage of the power representation of deep learning.

In this work, we propose to learn a regularization term designed to preserve the connectivity of vascular networks, and plug it into a variational iterative segmentation framework which does not require annotated data. We show that this regularization term better preserves the connectivity of blood vessels than classic regularization terms such as TV, and more complex ones including directional information.

II. MATERIALS AND METHODS

A. Learning a reconnecting term

We propose to learn a reconnecting network by training a model to reconnect blood vessels from images with artificial disconnections. This network takes as input a binary image of disconnected blood vessels and outputs an image with reconnected vessels. We chose to work on binary images for two reasons. First, creating realistic disconnections is easier on binary images. Second, our goal is to learn a reconnecting term independent of the dataset to apply it on new data for which ground-truths are not available. Learning on binary images is then a better strategy, as the model cannot learn any intensitybased feature of the training dataset.

We first created an algorithm to simulate realistic disconnections. It creates random disconnections in blood vessels with a thickness-dependent probability (the thinner the vessel, the higher the disconnection probability). The disconnections have a random shape to emulate disconnection due to image Fig. 1. Pipeline of our proposed method. A reconnecting network is first trained on a separate dataset based on a residual U-Net architecture. Then, a variational segmentation approach is used with the Chan [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] data fidelity term and either the total variation regularization term only or both the total variation and our proposed reconnecting term.

noise. Finally, some artefacts have been added to make the model distinguish parts of vessels that should be reconnected and noise that should not be connected to the vascular network. An example of disconnected image generated by this algorithm is presented in Figure 1.

To learn the reconnecting term, the residual U-Net [START_REF] Kerfoot | Left-ventricle quantification using residual U-Net[END_REF] from Monai1 , a Pytorch-based framework has been used. U-Net shows great performance in medical applications and the residual units increase its robustness to image variations. 20 segmentation ground-truths from the STARE [START_REF] Hoover | Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response[END_REF] dataset were used to create pairs of connected/disconnected vascular images to learn our reconnecting regularization term. We performed data augmentation by rotating each image 4 times and applying on each one different disconnections with our algorithm. This resulted in a dataset of 80 images. 2 images were kept for validation and the rest was used for training. We chose to input patches of size 96 × 96 instead of full images to the model to avoid learning the global architecture of the vascular network instead of the reconnecting feature. Additional random rotations and flips were applied on each patch and at each epoch to increase variability and prevent overfitting. A maximum of 1000 epochs, a batch size of 64, an Adam optimizer with a learning rate of about 10 -3 and a Dice loss were chosen.

B. Variational approach

Segmentation may be formulated as a variational problem consisting in the optimization of a two-term energy:

û = argmin u∈[0,1] N E data (u, f ) + λE reg (u), (1)
where f ∈ R N is the initial image, E data is the data fidelity, E reg the regularisation term and λ ∈ R a regularization coefficient acting as a trade-off between both terms. The form of E data depends on the choice of application. For segmentation the Chan et al. term [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] is classically used:

E data (u, f ) = u, c f F , c f → (c 1 -f ) 2 -(c 2 -f ) 2 , (2)
where c 1 and c 2 are constant values corresponding to the mean background and foreground intensity of the initial image, ., . F is the Frobenius product.

E reg constrains the problem with regularization and/or prior knowledge. A well-known term is the Total Variation (TV), which reduces the noise in the image, and is defined as follows:

T V (u) = λ||∇u|| 2,1 , (3)
with ||.|| 2,1 the l 1 norm of l 2 norm. To solve Eq (1) with a forward-backward primal-dual algorithm (FBPD) [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF], we can reformulate the model as follows:

û = argmin u h(u, f ) + g(Lu) + k(u), (4)
h(u, f ) = u, c f F , g(.) = λ||.|| 2,1 , L = ∇ and k(u) = ι u∈[0,1] N is the indicator function of the set [0, 1] N .
The FBPD algorithm is then :

u i+1 = prox τ k (u i -τ (∇h(x i ) + L T v i )) v i+1 = prox σg (v i + σL(2u i+1 -u i )), (5)
where prox g is the proximity operator of g, g is the conjugate of g, τ ∈ R + and σ ∈ R + are step-size parameters.

C. Decoupled approach

Our reconnecting term takes binary input images and returns a more connected vascular network. We thus propose to plug this reconnecting term after a few iterations of the algorithm presented in Equation 5, when the result u is close to a binary image. We thus replace, in Eq 5 prox σk (.) by Φ(.) defined as follows:

Φ(u) = prox σk (u) if i < α G reco (u) otherwise (6)
where G reco is the learnt reconnecting term, α ∈ N + * is the iteration number from which the reconnecting term is applied. Algorithm 1 present our proposed algorithm.

III. EXPERIMENTS

We applied our approach on the test set of the DRIVE dataset [START_REF]Ridge-based vessel segmentation in color images of the retina[END_REF], another retinal image dataset. We compared the results with the classic Chan et al. model [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF], and the more recent directional total variation model [START_REF] Merveille | nD variational restoration of curvilinear structures with prior-based directional regularization[END_REF]. We optimized the parameters of both compared methods on the DRIVE training set and selected the parameters yielding the best mean results. We set the parameters of our approach experimentally to α = 350 and λ = 0.005. Algorithm 1: Decoupled segmentation with learnt reconnecting term

Data: α ∈ N + * , u 0 ∈ R N 2 , v 0 ∈ R 2N 2 , (τ, σ) ∈ ]0, +∞[ 2 , λ n ∈]0, +∞[ for n ≥ 1 do p n = (u n -τ (∇h(u n ) + L T v n ) if n < α then p n = prox ι [0,1] (p n ) else p n = G reco (p n ) q n = prox σg (v n + σL(2p n -u n )) (u n+1 , v n+1 ) = (u n , v n ) + λ n ((p n , q n ) -(u n , v n ))
We performed a quantitative evaluation of the results by computing the mean sensitivity (TPR), specificity (TNR) and accuracy (Acc) over the 20 images of the test set (see Eq. 7). All metrics were computed in the field of view provided in the dataset for a fair comparison with the literature results.

with TP (resp. TN) the number of true positives (resp. negatives), and FP (resp. FN) the number of false positives (resp. negatives).

The quantitative results are shown in Table III. We also included the quantitative results from a recent deep learning approach [START_REF] Arias | A new deep learning method for blood vessel segmentation in retinal images based on convolutional kernels and modified U-Net model[END_REF] as a comparison. Illustrations of the results are presented in Fig. 2 and Fig. 3.

Qualitatively, we observe that our reconnecting term increases the connectivity of the whole vascular network. It also yields a better connectivity than the directional total variation. Moreover, blood vessels seem smoother than the ones obtained with the other two methods, which is for instance an appealing feature for modeling and blood flow simulations. Our results also exhibit less noise, as evidenced by the higher TNR. However, the mean TPR of our approach is lower than for directional TV. We observed that directional TV tends to detect fragmented vessel parts at the extremities of the network. Such fragments are irrelevant for applications requiring connectivity. By contrast, our approach tends to only preserve connected vessels: we found that our approach detected 230% less connected components than the directional TV approach. This translates into a lower number of true positives, but a better overall quantitative performance and qualitative results for connectivity-requiring applications.

We added the quantitative results from a recent deep learning approach [START_REF] Arias | A new deep learning method for blood vessel segmentation in retinal images based on convolutional kernels and modified U-Net model[END_REF] as a baseline for comparison. These approaches are trained with images from the same dataset and thus yield better results. However, these methods do not generalize well on new data (e.g. blood vessels from another organ or another modality). They need to be retrained on each new dataset with their corresponding annotations. Conversely, our approach only requires annotations to train the reconnecting term. Once it is trained, it can be applied in the variational segmentation framework on new data without requiring annotations. To show the power of generalization of our method, we applied the same reconnecting term used on the DRIVE dataset to a different type of image containing blood vessel-like structures. Figure 4 shows the segmentation results on this image. We observe that even though the initial image is different from the STARE dataset, our reconnecting term generalizes well and is able to reconnect the structures of interest whereas a trained U-Net has more difficulties to generalize to this new image type. In this case it seems that the network is not able to generalize on images with a larger range of vessel diameters.

IV. CONCLUSION

We developed a reconnecting term for blood vessel segmentation. This term learns to reconnect blood vessels by focusing on their geometry instead of learning grey-level features in the image. By plugging this term into a variational segmentation framework as a decoupled approach, we are able to better preserve the connectivity of vascular networks than standard regularization terms. Our approach can be applied to unseen data and performs well, unlike deep learning segmentation approaches that have to be trained on each new dataset. Future work will extend our reconnecting term in 3D for the segmentation of various 3D vascular networks, where annotations are far more complicated to obtain. We will also modify our reconnecting network to learn a maximally monotone operator to ensure the convergence of our algorithm [START_REF] Pesquet | Learning Maximally Monotone Operators for Image Recovery[END_REF].

Fig. 2 . 3 (Fig. 3 .

 233 Fig. 2. Comparison of segmentation results on one DRIVE image. Zooms in the red boxes are presented in Figure 3

  TP + TN + FN + FP ,

Fig. 4 .

 4 Fig. 4. Illustration of the generalization power of our reconnecting term.

TABLE I QUANTITATIVE

 I SEGMENTATION RESULTS ON THE DRIVE DATABASE.

		TPR	TNR	Acc
	Chan model [10]	0.6615 0.9833 0.9423
	Directional TV model [2]	0.6822	0.9809 0.9427
	Our approach	0.6631	0.9839	0.9429
	Deep learning approach (Arias et al. [15])	0.8597 0.9690	0.9563

This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, and by the IDEXLYON project within the program "Investissements d'Avenir" (ANR-11-IDEX-0007 and ANR-16-IDEX-0005) operated by the French National Research Agency (ANR).

https://monai.io/start.html


Christophe Zhang

Enrique Zuazua

A quantitative analysis of Koopman operator methods for system identification and predictions

Keywords: Koopman operator, system identification, finite elements spaces, data-driven approximation, Extended Dynamic Mode Decomposition

We give convergence and cost estimates for a data-driven system identification method: given an unknown dynamical system, the aim is to recover its vector field and its flow from trajectory data. It is based on the so-called Koopman operator, which uses the well-known link between differential equations and linear transport equations. Data-driven methods recover specific finite-dimensional approximations of the Koopman operator, which can be understood as a transport operator. We focus on such approximations given by classical finite-elements spaces, which allow us to give estimates on the approximation of the Koopman operator as well as the solutions of the associated linear transport equation. These approximations are thus relevant objects to solve the system identification problem.

We then analyze the convergence of a variant of the generator Extended Dynamic Mode Decomposition (gEDMD) algorithm, one of the main algorithms developed to compute approximations of the Koopman operator from data. We find however that, when combining this algorithm with classical finite elements spaces, the results are not satisfactory numerically, as the convergence of the data-driven approximation is too slow for the method to benefit from the accuracy of finite elements spaces. In particular, for problems in dimension 1 it is less efficient than direct interpolation methods to recover the vector field. We provide some numerical examples to illustrate this last point.

Introduction

Unknown dynamics and system identification

Many problems in applied mathematics rely on the existence of a mathematical model for the system under consideration. Given the potential complexity of systems that arise in all scientific fields, from econometry to molecular dynamics, finding relevant and numerically tractable models for applications such as forecasting and control is a central challenge.

The traditional approach to mathematical modelling consists in deriving governing equations from fundamental laws of physics such as Newton's second law or the principle of least action, after a series of assumptions and ideal simplifications. Models derived in this way have been thoroughly studied, giving rise to rich theories and a profound understanding of some of the phenomena they describe. But these remain ideal models, which sometimes fall short of real-life observations and applications. Indeed, as real-life systems become more complex, the relevance of such models for applications comes into question, as there is a risk that usual, convenient simplifications could lead to essential features of the dynamics being overlooked.

On the other hand, the development of data-driven techniques has brought a new perspective to this challenge (see for example [START_REF] Lennart | System identification: theory for the user[END_REF][START_REF] Bamieh | Identification of linear parameter varying models[END_REF][START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF]): how can one complement traditional modelling approaches with data-driven considerations? Can relevant features of the dynamics of the system be recovered from data? More generally, how can one learn a model from data, and propose a system of governing equations that allow for reliable predictions, and effective control?

Mathematically speaking, we can cast this as follows: suppose that, for a given system, its behaviour can be relevantly modelled with an autonomous, nonlinear differential equation

ẋ = f (x), x ∈ R d . (1)
How can one learn an approximation of the vector field f from observational data, that is, samples of trajectories? Given initial conditions, can one produce predictions on the evolution of the system in a data driven way? This is a form of inverse problem, which differs from more classical inverse problems (see for example [START_REF] Alberto P Calderón | On an inverse boundary value problem[END_REF][START_REF] Isakov | Inverse problems for partial differential equations[END_REF][START_REF] Yamamoto | Carleman estimates and an inverse heat source problem for the thermoelasticity system[END_REF]) in that the latter rely on the a priori choice of a class of models (heat equation, polynomial differential equations, balance laws...) so that one has to recover parameters of the given model (diffusivity, coefficients...).

Regression, interpolation and operator theoretic approach

The task of recovering a function f : R d → R d from a data sample of values {x i , y i = f (x i )} N i=1 in order to make predictions lies at the heart of statistical regression analysis (linear regression being perhaps the simplest example), with examples dating back as early as Legendre's works on planetary orbits using a least squares method ( [START_REF] Marie | Nouvelles méthodes pour la détermination des orbites des comètes[END_REF]). More recently, a variety of supervised learning and deep learning techniques have been developed to tackle this problem, such as kernel-based methods, neural networks, or deep neural networks. We refer to [START_REF] Binev | Universal algorithms for learning theory part i: piecewise constant functions[END_REF][START_REF] Binev | Universal algorithms for learning theory. part ii: Piecewise polynomial functions[END_REF][START_REF] Goodfellow | Deep Learning[END_REF] for a comprehensive overview and further references.

The core idea of these regression techniques is to find a function f belonging to a certain class of functions S (linear functions, polynomials, wavelets...), that minimizes the loss

N i=1 f (x i ) -y i 2 , (2)
in such a way that f will also be a good approximation of f on new data. The first approach to solving (2) would be to find, if it exists, a function in S that interpolates the data, that is,

f ∈ S, N i=1 f (x i ) -y i 2 = 0 i.e. f (x i ) = y i = f (x i ), i = 1, • • • , N. (3)
A typical example of functions used for interpolation are polynomials, as there exists a unique polynomial of degree at most N -1 that satisfies [START_REF] Arbabi | Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator[END_REF]. We refer to [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]Chapter 3] for more details on these interpolation results. However, although f fits the given data perfectly, it has no reason a priori to be a good approximation of f on new data. Nonetheless under some assumptions on the x i and on f , there are some well-known estimates for interpolation by polynomials and continuous piecewise polynomial functions (we refer again to [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]Chapter 3]).

In general, regression techniques are not designed to interpolate the data exactly, as in the above, but rather to reach a compromise between a small loss function in [START_REF] Arbabi | A Data-Driven Koopman Model Predictive Control Framework for Nonlinear Partial Differential Equations[END_REF] and a good overall approximation of the target function f . The goal is to ensure a good generalization of f to new data, without resorting to the aforementioned general interpolation estimates. There are then two important aspects to consider: first, the choice of the family of functions S and its representation, which can benefit from prior knowledge on f . Second, the method used to minimize the loss, which can involve for example regularization to avoid overfitting f to the data to the detriment of generalization.

In particular, if one has the prior knowledge that f is linear, then the above problem amounts to a multivariate linear regression, and the least squares problem (2) has a well-known solution. This, among other things, has motivated an alternative approach system identification, the so-called operator theoretic approach. Instead of treating system identification directly as a general nonlinear problem, one can recast this problem as a linear one, using the well-known connection between differential equations and linear transport equations. This is appealing both on a theoretical level, as it gives an elegant framework to understand some intricate properties of nonlinear dynamical systems, and on a practical level, due to the many tools available both for forecasting and control of linear systems. As we will see, with this approach system identification then becomes a linear regression, but in a potentially much larger space.

Koopman operators

Definitions

Consider a vector field f : R d → R d . One can view it as the vector field of a nonlinear ordinary differential equation:

ẋ = f (x), x ∈ R d , (4)
but also as the velocity field of the linear transport equation:

∂ϕ ∂t (t, x) = f (x) • ∇ϕ. (5)
The equivalence between these two points of view has been long known. Existence and uniqueness results for nonlinear differential equations with Sobolev vector fields were obtained by studying the corresponding transport equation ( [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF][START_REF] Hauray | Deux remarques sur les flots généralisés d'équations différentielles ordinaires[END_REF]). Conversely, the method of characteristics solves (5) using ( 4): if f is regular enough, it generates a flow Φ t f which solves (4), and the solutions of ( 5) are constant along the characteristics (t, Φ t f (x)), t ≥ 0, x ∈ R d . The solution of ( 5) with initial condition

ϕ 0 ∈ L 2 (R d ) is then given by ϕ(t, x) := ϕ 0 Φ t f (x) . (6)
This corresponds to the following semigroup:

ϕ 0 ∈ L 2 (R d ) → ϕ 0 • Φ t f , (7)
which is a simple composition operator with the following infinitesimal generator

K := f (x) • ∇, (8)
which is simply the transport operator with velocity field f , on the domain

D(K) = {ϕ ∈ L 2 (R d ), f • ∇ϕ ∈ L 2 (R d )}. (9)
In what follows, we will denote by e tK the semigroup [START_REF] Binev | Universal algorithms for learning theory. part ii: Piecewise polynomial functions[END_REF], generated by the Koopman operator K.

The semigroup and its generator were studied and used by B.O.Koopman in [START_REF] Koopman | Hamiltonian Systems and Transformation in Hilbert Space[END_REF] to study properties of fluid dynamics. They have since then been named after him: what is known as the Koopman operator of f is the operator K given by [START_REF] Binev | Universal algorithms for learning theory part i: piecewise constant functions[END_REF], and the Koopman semigroup of f is the semigroup e tK given by [START_REF] Binev | Universal algorithms for learning theory. part ii: Piecewise polynomial functions[END_REF]. Practically speaking, for a function ϕ ∈ L 2 (R d ) and a trajectory x(t) of the differential equation ( 4), the Koopman semigroup describes the evolution of t → ϕ(x(t)).

Remark 2.1. In applications involving discrete-time dynamical systems

x + = T (x), x ∈ R d , (10)
the corresponding discrete-time Koopman operator is simply the composition operator

U : g → g • T. (11)
Equivalently, when dealing with continuous-time systems, the Koopman operator is sometimes defined as the bounded operator e τ K with the notations of (7), for some time increment τ > 0.

Historically, B.O.Koopman used the linear composition operator [START_REF] Binev | Universal algorithms for learning theory. part ii: Piecewise polynomial functions[END_REF] to study measure-preserving systems. It was notably instrumental in J. Von Neumann's proof of the mean ergodic theorem [START_REF] Neumann | Proof of the quasi-ergodic hypothesis[END_REF], and a first connection of chaotic behavior to the existence of continuous parts in the spectrum of the Koopman semigroup was noted in [START_REF] Bo Koopman | Dynamical systems of continuous spectra[END_REF]. More recently, spectral properties of the Koopman semigroup have been shown to be related to a variety of aspects of the underlying dynamical system such as ergodicity, stability, existence of different time scales, mixing and non-mixing properties (see for example [START_REF] Mezić | Spectral Properties of Dynamical Systems, Model Reduction and Decompositions[END_REF][START_REF] Mezić | Analysis of fluid flows via spectral properties of the koopman operator[END_REF][START_REF] Lan | Linearization in the large of nonlinear systems and koopman operator spectrum[END_REF][START_REF] Mauroy | Global Stability Analysis Using the Eigenfunctions of the Koopman Operator[END_REF][START_REF] Igor | Spectrum of the Koopman Operator, Spectral Expansions in Functional Spaces, and State-Space Geometry[END_REF][START_REF] Kvalheim | Existence and uniqueness of global koopman eigenfunctions for stable fixed points and periodic orbits[END_REF]). An important feature in these works that the spectral analysis of the Koopman operator is done in a rather geometrical spirit in combination with linearization results and conjugation techniques, focusing on local behaviour near equilibria or attractors.

Koopman operator and data-driven modelling

In addition to this renewed theoretical exploration, data-driven applications of Koopman operator theory have gained quite some interest in the last decade, in particular for system identification, modelling, predictions and control problems [START_REF] Steven | Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control[END_REF][START_REF] Korda | Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control[END_REF][START_REF] Klus | Data-driven approximation of the koopman generator: Model reduction, system identification, and control[END_REF]. In these works, among others, the idea is that one can identify the underlying nonlinear dynamical system, (or at least some features of this system), by identifying its Koopman operator, or its Koopman semigroup, or the discrete-time Koopman operator, depending on the nature of the system.

There are two potential benefits to this approach:

1. The linearity of the Koopman operators, which offers the possibility to describe nonlinear dynamics in a linear way.

2. As a consequence, the possibility to make predictions by computing solutions to the linear equations given by these operators, instead of integrating a nonlinear differential equation.

Numerically speaking, the focus of these methods lies on finite-dimensional reductions of the Koopman operator, as it is not possible to recover an infinite-dimensional operator from a finite amount of data.

The question then becomes: which finite-dimensional reductions are accurate and relevant for system identification and predictions? What is their dimension? How does one recover them from data? At what cost?

• Exact reductions. In some cases, there exist finite-dimensional subspaces V of functions which are invariant under the action of the Koopman operator. Then, the linear transport equation ( 5) reduces exactly to a linear differential equation on these subspaces:

φ = Kϕ, ϕ ∈ V.
The main example of such a situation can be found in [START_REF] Steven | Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control[END_REF], where a system of 2 polynomial differential equations are lifted to a system of 3 linear differential equations by considering a basis of multinomials (note that this is a particular case where the Carleman linearization [START_REF] Carleman | Application de la théorie des équations intégrales linéaires aux systèmes d'équations différentielles non linéaires[END_REF] of a polynomial system is exact and does not need to be truncated).

Although the existence of such subspaces holds great interest for applications, and there have been efforts to recover them from data [START_REF] Takeishi | Learning koopman invariant subspaces for dynamic mode decomposition[END_REF], generally speaking, this seems to be a very rare situation.

Let us point out also that even though in this example, the dimension of the invariant subspace is small, in general it can be much larger the dimension d of the underlying problem.

Another way to find invariant subspaces is to find eigenfunctions of the Koopman operator. Then, any subspace spanned by a finite number of eigenfunctions ϕ 1 , • • • , ϕ N with eigenvalues λ 1 , • • • , λ N will be invariant under the action of the Koopman operator, and the linear transport equation ( 5) reduces exactly to the following diagonal system of linear differential equations:

φi = λ i ϕ i , i ∈ 1, N .
However, even if there exist eigenfunctions to the Koopman operator, there is no reason a priori for them to span an approximation of the Koopman operator that is relevant for system identification.

• Spectral properties of the Koopman operator. On the other hand, the complete spectral analysis of the Koopman operator is a delicate topic, and depends heavily on the choice of a function space. The transport operator

K = f (x) • ∇, D(K) = H 1 (R d )
does not have a natural basis of eigenfunctions, sometimes no eigenfunctions at all. This is related to the fact that in general transport operators do not have a compact resolvent, and are not generally self-adjoint unless the velocity field f has zero divergence.

Let us nevertheless point out an interesting recent development for the discrete time Koopman operator introduced in Remark 2.1. In the case where Ω is the unit circle in R2 , [START_REF] Slipantschuk | Dynamic mode decomposition for analytic maps[END_REF] shows that an adequate choice of a function space (namely, Hilbert-Hardy spaces) leads to interesting spectral results. Indeed, by restricting the function space to a class of analytical functions, the adjoint of the Koopman operator (the Perron-Frobenius operator) becomes compact. By duality, this means that in this case, by extending the discrete time Koopman operator to a larger space containing L 2 , it becomes compact. When it is self-adjoint, which is equivalent to the existence of an invariant measure, it then has a basis of eigenfunctions. This result however concerns a certain class of discrete-time dynamical systems, and analogous results in the continuous time setting, regarding the Koopman generator [START_REF] Binev | Universal algorithms for learning theory part i: piecewise constant functions[END_REF], have yet to be explored.

In terms of data-driven methods, spectral analysis of the Koopman operator from data has received a lot of attention (see [START_REF] Kaiser | Data-driven discovery of koopman eigenfunctions for control[END_REF][START_REF] Korda | Optimal construction of Koopman eigenfunctions for prediction and control[END_REF][START_REF] Govindarajan | On the Approximation of Koopman Spectra for Measure Preserving Transformations[END_REF][START_REF] Korda | Data-driven spectral analysis of the Koopman operator[END_REF]), in particular the data-driven recovery of eigenfunctions when they exist. As is often the case in numerical analysis, only the smallest eigenvalues can be accurately approximated. This is of interest to isolate components of the observed system that evolve slowly in time (or are invariant, if the eigenvalue is 0). However, these give only a very partial picture of how the system actually evolves.

As we have pointed out earlier, in general, transport operators (Koopman operators) do not have nice spectral properties. For system identification and predictions, there seems to be better hope in finding accurate and appropriate finite-dimensional approximations of the Koopman operator independently of such properties.

• Galerkin projections and Extended Dynamic Mode Decomposition. Classical Galerkin methods used in numerical analysis, which approximate solutions of partial differential equations using approximation spaces, seem a natural solution to provide reliable approximations of the Koopman operator. They have indeed proven to be an efficient way to obtain reduced-order models of partial differential equations without having to use spectral decompositions.

To recover a data-driven approximation of such a Galerkin projection, one of the main algorithms is the so-called Extended Dynamic Mode Decomposition (EDMD) (introduced in [START_REF] Matthew | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF]). Its predecessing algorithm, Dynamic Mode Decomposition (DMD), computes the "best fit" system of linear differential equations for the unknown system (4) in the state space R d . EDMD does this in the subspace of functions given by the Galerkin projection. As this corresponds to recovering the best linear approximation of a linear infinite-dimensional operator, this has a better chance of yielding convincing results than DMD.

EDMD was first introduced in the discrete-time framework, to recover approximations of the discrete-time Koopman operator, that is, the classical composition operator. It has been actively explored as a way to perform data-driven system identification and predictions [START_REF] Korda | On Convergence of Extended Dynamic Mode Decomposition to the Koopman Operator[END_REF][START_REF] Lusch | Deep learning for universal linear embeddings of nonlinear dynamics[END_REF][START_REF] Klus | Kernel-Based Approximation of the Koopman Generator and Schrödinger Operator[END_REF]. It can also provide approximations of so-called coarse-grain models (see [START_REF] Kmiecik | Coarse-Grained Protein Models and Their Applications[END_REF]) of highly complex systems such as molecular dynamics ( [START_REF] Klus | Data-driven approximation of the koopman generator: Model reduction, system identification, and control[END_REF]). It has been used in combination with model predictive control to achieve data-driven control of systems with unknown dynamics [START_REF] Korda | Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control[END_REF], and combined with deep learning techniques [START_REF] Arbabi | A Data-Driven Koopman Model Predictive Control Framework for Nonlinear Partial Differential Equations[END_REF][START_REF] Ping | Deep Koopman model predictive control for enhancing transient stability in power grids[END_REF][START_REF] Korda | Power grid transient stabilization using koopman model predictive control[END_REF][START_REF] Otto | Linearly-recurrent autoencoder networks for learning dynamics[END_REF][START_REF] Yeung | Learning deep neural network representations for koopman operators of nonlinear dynamical systems[END_REF].

In continuous time, there are two main approaches using EDMD to recover an approximation of the Koopman operator.

In this article, we implement a variant of gEDMD, using classical finite elements spaces to define finite-dimensional approximations of the Koopman operator. Our central aim is to quantify the quality and cost of system identification and predictions by this method. More precisely, we will focus on the case of an unknown vector field f of a certain regularity, on a bounded domain Ω ∈ R d , as it would be unreasonable to expect in general a good approximation of f on all of R d from a finite number of observations. Moreover, to simplify matters and alleviate notations, we will focus on the case where the flow

Φ t f satisfies Φ t f (Ω) ⊂ Ω. (12)
In Section 3, we give estimates on the finite-dimensional approximations of the Koopman operator, thanks to properties of finite elements spaces. We show how these approximations can be used to solve the system identification problem. The estimates depend on the size of the mesh and on the degree of the finite elements functions we consider, which are related to the dimension of the finite elements space. Accurate Galerkin projections require finite elements spaces of considerable dimension, which increases exponentially with the dimension d of the system under consideration.

In Section 4 we devise a variant of gEDMD to compute a data-driven approximation of these Galerkin projections. The convergence of this approximation depends on the number of data samples, and is essentially based on Monte-Carlo approximation of integrals.

In Section 5 we then analyze the total approximation error of gEDMD with finite elements, and estimate its computational cost, in order to quantify whether this method can be efficient for system identification and predictions in practice.

This allows us to give some clear insights on the challenges and limitations of this method. The overall conclusion is that, even when using accurate approximation spaces such as finite elements, the approximation error decreases relatively slowly, because on one hand the convergence of the data-driven approximation is quite slow (in m -1 2 , m being the number of samples), but also because this convergence is further slowed if one chooses a finite elements space of greater dimension (i.e. refines the mesh or considers higher-order finite elements).

Finally in Section 6 we give some numerical illustrations in dimension d = 1 that illustrate the clear advantage of direct interpolation methods, given the same data sample, due to the slowness of the convergence of gEDMD.

Notations

We denote the usual Lebesgue measure by µ.

For the rest of this article, when there is no ambiguity we will note

L 2 (R d , µ) = L 2 to alleviate notations. Notice that if f ∈ L ∞ (R d , R d )
, by the definition [START_REF] Steven | Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control[END_REF] we have

H 1 (R d ) ⊂ D(K). ( 13
)
We denote by

• the usual L 2 norm, • H s (R d ) the usual Sobolev norms for s ≥ 0, and | • | H k (R d ) the Sobolev seminorms for k ∈ N: |ϕ| 2 H k (Ω) = α∈N d |α|=k D α ϕ 2 , ϕ ∈ H k (Ω). (14)
In Euclidean spaces, • denotes the usual Euclidean norm.

For matrices ||| • ||| denotes the operator norm induced by the Euclidean norm, • F denotes the Frobenius norm, • ∞ denotes the max norm

A ∞ = max ij |A ij | .
The sign † denotes the Moore-Penrose pseudo-inverse, and will denote the usual transposition.

3 Finite-dimensional approximations of the Koopman operator

Restricting the Koopman operator to a subdomain Ω

As we have mentioned in the introduction, we will focus on studying the Koopman operator of a vector field f on a bounded domain Ω ⊂ R d . Indeed, from a finite data sample one cannot expect to obtain a good approximation of f on all of R d .

When dealing with a vector field that is indeed defined on all of R d , this means we are studying the Koopman operator K Ω of the restriction f |Ω :

K Ω = f |Ω • ∇. (15)
Again, K Ω is an unbounded operator, with domain

D(K Ω ) = {ϕ ∈ L 2 (Ω), f • ∇ϕ ∈ L 2 (Ω)}. (16)
However, with this domain, in general K Ω does not generate a semigroup as in [START_REF] Bamieh | Identification of linear parameter varying models[END_REF]. Indeed in general the flow Φ t f does not preserve Ω. Noting

∂Ω out = {x ∈ ∂Ω, f (x) • ν(x) > 0}, (17)
where ν(•) denotes the outward normal unit vector, it is necessary to specify boundary conditions on ∂Ω out in order for the corresponding linear transport equation to be well posed. Typically, the partial differential equation

   ∂ϕ ∂t = f (x) • ∇ϕ, x ∈ Ω, ϕ(t, x) = 0 x ∈ ∂Ω out , (18)
is well posed. Accordingly, the unbounded operator

KΩ := f • ∇, D( KΩ ) = {ϕ ∈ H 1 (Ω), ϕ = 0 on ∂Ω out }, (19)
generates the C 0 -semigroup

e t KΩ ϕ 0 (x) = ϕ 0 (Φ t f (x)) if Φ t f (x) ∈ Ω, 0 otherwise, ∀ϕ 0 ∈ L 2 (Ω), ∀x ∈ Ω. (20)

System identification and predictions with the Koopman operator

We now show how the Koopman operator can be used to recover the vector field f , and in some cases, its flow Φ t f . Let us note the coordinate functions

g i (x) = x i , x ∈ Ω, i = 1, • • • , d. (21)
Then, as they are polynomial, and Ω is a bounded domain,

g i ∈ D(K Ω ), i = 1, • • • , d. (22)
Moreover,

K Ω g i = f • ∇g i = f i , i = 1, • • • , d, (23)
where the f i denote the coordinates of the vector field f in (4). It is thus possible, with the restricted Koopman operator K Ω , to recover the vector field f . On the other hand, in general

g i / ∈ D( KΩ )
due to the boundary conditions in [START_REF] Goodfellow | Deep Learning[END_REF]. However, we can still apply the semigroup defined in ( 20):

e t KΩ g i (x) = (Φ t f (x)) i if Φ t f (x) ∈ Ω, 0 otherwise, ∀x ∈ Ω, (24)
where the (Φ t f (•)) i denote the coordinates of the flow (4). Thus, it is possible, with the restricted Koopman operator KΩ , to recover the flow Φ t f , and thus make predictions, as long as the flow does not exit Ω.

To make matters simpler, in order to bring forth our conclusions with more clarity, we will for the rest of this article focus on the case where µ ∂Ω ∂Ω out = 0, i.e., the case where Ω is forward invariant under Φ t f . Then, we will set, to alleviate notations,

K := K Ω = f • ∇, D(K) := D(K Ω ) = D( KΩ ).
With this domain, K generates the semigroup:

e tK ϕ 0 (x) = ϕ 0 Φ t f (x) , ∀ϕ 0 ∈ L 2 (Ω), ∀x ∈ Ω. (25)
i.e. the linear transport equation

   ∂ϕ ∂t (t, x) = f (x) • ∇ϕ(t, x), t ≥ 0, x ∈ Ω, ϕ(0, •) = ϕ 0 (•) ∈ L 2 (Ω), (26)
is well-posed.

In particular, (23) remains unchanged, and ( 24) becomes

e tK g i (x) = Φ t f (x) i , ∀x ∈ Ω. ( 27
)

Galerkin projection of the Koopman operator

We now consider a bounded domain Ω ⊂ R d , and we focus on the case where f is such that

Φ t f (Ω) ⊂ Ω.
On this domain, the Koopman operator of f

K = f • ∇ with domain D(K) = H 1 (Ω)
generates the semigroup

e tK ϕ 0 = ϕ 0 • Φ t f , ∀ϕ 0 ∈ L 2 .
The idea of Koopman operator methods is to observe trajectories of the dynamical system (4), but instead of trying to identify the vector field f directly, to apply linear regression methods to identify the linear operator K.

In practice, identifying K numerically is impossible as it is infinite dimensional. Instead, the goal is to identify a Galerkin projection of K, which corresponds to a reduced-order model of the linear transport equation [START_REF] Bamieh | Identification of linear parameter varying models[END_REF]. The projection is determined by a family of linearly independent functions

Ψ := {ψ 1 , • • • , ψ N } ⊂ L 2 , which span a finite-dimensional subspace V N ⊂ L 2 which should satisfy V N ⊂ D(K). (28)
Then, the Galerkin projection of K is given by

K N := Π N KΠ N , (29)
where Π N denotes the L 2 -orthogonal projection on V N . In particular,

K N ϕ = Π N Kϕ, ϕ ∈ V N . (30)
In this article we will focus on a specific choice of V N , given by continuous finite elements. Classical approximation results for these spaces then ensure the Galerkin projection is a relevant approximation of the actual Koopman operator.

Galerkin projection on finite elements spaces

In order to alleviate notations, and simplify matters pertaining to the choice and regularity of meshes, we fix for the remainder of this article

Ω = Q := [0, L] d (31)
We then define the family of rectangulations of Q:

T M = d i=1 p i L M , (p i + 1) L M , 0 ≤ p i ≤ M -1, M ∈ N * . ( 32
)
In what follows, we will note

T h := T M where h = L M .
These rectangulations form a regular family of rectangulations of Q.

For k ≥ 1, we define Q k the space of polynomials of degree at most k in each variable, and the following continuous finite elements spaces:

V k h = v ∈ C 0 (Q), v |K ∈ Q k , ∀K ∈ T h . (33)
We consider the nodes

xα = h k α, α ∈ 0, k L h d ,
and note

N k h = (1 + kM ) d = 1 + kL h d ( 34
)
in order to re-index the nodes by xj , j ∈ 1, N k h . Then the classical shape functions ψ j which satisfy

ψ j (x k ) = δ jk ∀j, k ∈ 1, N k h , form a basis of V k h , so that dim V k h = N k h , (35)
i.e.,

h = kL N k h 1 d -1 . ( 36
)
One can see from [START_REF] Kunsch | Breaking the curse for uniform approximation in Hilbert spaces via Monte Carlo methods[END_REF] that both refining the mesh (h → 0) and taking finite elements of higher degree (k → ∞) will increase the dimension of the approximation subspace at the same rate.

We now recall classical projection inequalities for finite elements spaces. Note r k h the classical finite elements interpolation operator

r k h (ϕ) = N k h j=1 ϕ(x j )ψ j , ϕ ∈ C 0 (Ω). (37)
The classical interpolation inequalities for triangular finite elements (see for example Quarteroni and Valli [53,Section 3.4]) can be adapted to the rectangulation T h : Proposition 3.1. Let T h be the rectangulation defined by [START_REF] Korda | Optimal construction of Koopman eigenfunctions for prediction and control[END_REF]. Let s > d/2, ν ∈ {0, 1}, and define l = min(k, s -1). Then, for ϕ ∈ H s (Ω), ϕ ∈ C 0 (Ω) so that r k h (ϕ) is well defined, and there exists C > 0 independent of h such that

r k h (ϕ) -ϕ H ν (Q) ≤ Ch l+1-ν |ϕ| H l+1 (Q) , ∀ϕ ∈ H s (Q). ( 38
)
Denoting Π k h the L 2 -orthogonal projection on V k h , one gets, by property of orthogonal projections in Hilbert spaces: Corollary 3.1. With the notations above,

Π k h (ϕ) -ϕ H ν (Q) ≤ Ch l+1-ν |ϕ| H l+1 (Q) , ∀ϕ ∈ H s (Q). (39)
Now, we have

V k h ⊂ D(K). (40)
Following [START_REF] Bo Koopman | Dynamical systems of continuous spectra[END_REF], the Galerkin projection of K is then given by

K k h := Π k h KΠ k h .

System identification and predictions

We know show how the approximation power of finite elements can be combined with ( 23) and ( 27) to obtain approximations of the vector field f and its flow Φ t f . First, as the g i are linear, we clearly have

g i ∈ V k h , ∀i ∈ 1, d . (41)
We then have, from [START_REF] Isakov | Inverse problems for partial differential equations[END_REF],

K k h g i = Π k h Kg i = Π k h f i , (42)
which can be understood as a finite elements approximation of [START_REF] Isakov | Inverse problems for partial differential equations[END_REF].

In the same spirit, we obtain a finite elements approximation of ( 26) by the Galerkin method:

   d dt ϕ k,h i = K k h ϕ k,h i , ϕ k,h i ∈ V k h , ϕ k,h i (0, •) = g i (•), i ∈ 1, d . (43)
Remark 3.1. In accordance with the assumptions we have made on the domain Ω and the behavior of the flow, there is no need here to specify boundary conditions for the Galerkin method, as the actual linear transport equation ( 26) is well-posed without having to specify any boundary conditions.

Noting

ϕ k,h i (t, •) = N k h j=1 (ϕ k,h i (t)) j ψ j , g i = N k h j=1 (g i ) j ψ j ,
We have the following system of linear differential equations

   M k h d dt ϕ k,h i = R k h ϕ k,h i , ϕ k,h i (0) = g i , (44)
where

M k h = ( ψ i , ψ j ) i,j , R k h = ψ i , K k h ψ j i,j
are the mass and stiffness matrices of the basis (ψ j ), for the operator K k h : noting K k h its matrix in the basis (ψ j ), we have

K k h = M k h -1 R k h . (45)
We see from ( 42) and ( 43) that solving the system identification problem with a Galerkin projection corresponds to projecting the vector field on the finite elements space V k h , and approximating the solutions of the linear transport equation [START_REF] Klus | Data-driven approximation of the koopman generator: Model reduction, system identification, and control[END_REF] as in classical numerical analysis. With this in mind, the choice of classical approximation spaces used in numerical analysis for the approximation of partial differential equations seems adequate.

Accordingly, under some conditions on the order k of the finite elements, and the regularity of the vector field f , we have the following estimates both on the approximation of f and its flow Φ t f :

Proposition 3.2. Let k + 1 > d/2. There exist constants C 1 , C 2 (f ) > 0 such that, for f ∈ C k+1 (Ω, R d ),
the following error estimates hold:

K k h g i -f i ≤ C 1 h k+1 |f i | k+1 . ( 46
)
Φ t f i -ϕ k,h i (t, •) ≤ C 2 (f )h k . (47)
Proof. From (42), we have

K k h g i -f i = Π k h f i -f i . (48)
To estimate the right hand-term, we apply Corollary 3.1 with s = k + 1 to the f i ∈ C k+1 (Ω) ⊂ H k+1 (Ω) (as Ω is bounded), which directly yields [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF]. Now consider the solutions of the following linear transport equations:

   ∂ϕ i ∂t (t, x) = f (x) • ∇ϕ i (t, x), t ≥ 0, x ∈ Ω, ϕ i (0, •) = g i (•), i ∈ 1, d , (49)
Recalling [START_REF] Kmiecik | Coarse-Grained Protein Models and Their Applications[END_REF], these solutions are given by

ϕ i (t, x) = g i Φ t f (x) = Φ t f (x) i . (50)
As f ∈ C k+1 (Ω, R d ), by the Cauchy-Lipschitz theorem, for all i ∈ 1, d ,

(t, x) → Φ t f (x) i (51) is in C k+1 ([0, T ] × Ω). Hence, ϕ i ∈ L 2 0, T ; H k+1 (Ω) ∩ H 1 0, T ; H k (Ω)
. Adapting [53, Chapter 14, Section 3] to our setting, which is simplified by the absence of boundary conditions and the fact that ϕ k,h i (0, •) ∈ V k h , we can then estimate the convergence of solutions of (43) to the ϕ i : there exists a constant C 2 (f ) > 0 such that

ϕ i (t, •) -ϕ k,h i (t, •) ≤ C 2 (f )h k . ( 52
)
Together with [START_REF] Penrose | A generalized inverse for matrices[END_REF], this directly yields [START_REF] Neumann | Proof of the quasi-ergodic hypothesis[END_REF].

Thus, Galerkin projections of the Koopman operator on finite elements spaces are relevant objects to approximately solve the system identification problem, both to recover the vector field f and to approximate its flow Φ t f . However the accuracy obtained in the estimates ( 46) and ( 47) comes at a cost, as we will see in the next section.

Curse of dimensionality

In Proposition 3.2 we have given the error estimate [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF] for system identification using the Galerkin projection K k h of K on finite elements spaces. In terms of the dimension of the finite elements space, we have the following: Proposition 3.3. The dimension N k h of the finite elements spaces needed to ensure the quality of approximation ε > 0 for system identification grows exponentially with the dimension d of the underlying space:

N k h ≥ C (L, f )k d ε -d k+1 , (53)
for some constant C (L, f ) > 0.

Proof. If one requires the quality of approximation in [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF]:

C 1 h k+1 |f | k+1 ≤ ε, (54)
then h ≤ C (f )ε 1 k+1 . ( 55
)
From ( 34) we then get

N = 1 + kL h d ≤ C (L, f )k d ε d k+1 , (56)
which proves [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF].

This is a form of curse of dimensionality, a phenomenon first pointed out in [START_REF] Bellman | Dynamic Programming[END_REF] in the context of dynamic programming. It is present in a wide range of situations such as function approximation ( [START_REF] Devore | Nonlinear approximation[END_REF][START_REF] Kunsch | Breaking the curse for uniform approximation in Hilbert spaces via Monte Carlo methods[END_REF]), integral computation ( [START_REF] Novak | The curse of dimension and a universal method for numerical integration[END_REF]), approximation of parametric PDEs ( [START_REF] Bungartz | Sparse grids[END_REF][START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF]), where the complexity of a problem increases exponentially with its underlying dimension. On the other hand, for a fixed dimension d, the required dimension N k h is polynomial in the error ε. This curse of dimensionality is a common feature in approximation theory. It stems from the use of a mesh to define the approximation spaces (here, continuous functions that are polynomial on each cell of the mesh). This allows a very accurate approximation, but the complexity of the mesh increases exponentially with the dimension, hence the curse of dimensionality.

Recalling that predictions are made by integrating the system of N k h linear differential equations [START_REF] Mezić | Analysis of fluid flows via spectral properties of the koopman operator[END_REF], this means that in high dimension d, making predictions by using Galerkin projections of the Koopman operator becomes prohibitively expensive.

In the next sections we will also see that this curse of dimensionality has an impact on the data-driven approximation of the Galerkin projections.

A variant of generator Extended Dynamic Mode Decomposition

Having given estimates on how the Galerkin projection K k h converges to the Koopman operator K under suitable regularity assumptions on the vector field f , it remains to compute a data-driven approximation of K k h . To achieve this, we present a modified version of the generator Extended Dynamic Mode Decomposition (gEDMD) algorithm.

The algorithm

Recalling the identity (45), we see that there two distinct components in the Galerkin projection K k h : the mass matrix M k h , which depends on the space V k h and the basis (ψ j ) (thus, essentially, on the domain Ω), and the stiffness matrix R k h which depends on the vector field f , the space V k h , and the basis (ψ j ). In general, the mass matrix can be approximated by numerical integration methods such as quadrature formulae, or Monte-Carlo integration. In the case of finite elements, in particular linear finite elements, it can even be given explicitly. In any case, this does not require trajectory data, as it does not depend on f , and can thus be done separately.

On the other hand, the stiffness matrix R k h depends on f . The algorithm we now lay out mainly consists in approximating this matrix from observations of trajectories of the system (4). To our knowledge, this approach, consisting in separating the approximation of the mass and stiffness matrices, is new. As we will see, it has the advantage of allowing the gEDMD method to be extended to localized bases such as finite elements, or wavelets.

Practically speaking, we consider a collection of points {x l ∈ Ω, l ∈ 1, m }, and for each l ∈ 1, m , the solution x l (t) of (4) with initial condition x l .

To account for the potentially random nature of measurements for real-life systems, we make the common assumption that the x l are drawn independently and uniformly with respect to the Lebesgue measure on Ω. From an adequate sampling of the corresponding trajectories x l (t), l ∈ 1, m , the initial velocities ẋk (0) = f (x k ) can be approximated by total variation regularized differentiation, which is also used in other methods such as SINDy (see [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Chartrand | Numerical differentiation of noisy, nonsmooth, multidimensional data[END_REF]). We then proceed as follows:

Step 1 From the data set {x l , y l := f (x l ), l ∈ 1, m }, define the following matrices:

G m k,h = ψ i (x j ) i∈ 1,N j∈ 1,m , A m k,h = y j • ∇ψ i (x j ) i∈ 1,N j∈ 1,m , (57)
and compute the following:

R m k,h := 1 m G m k,h A m k,h . (58)
As we will establish later on, R m k,h is an approximation of the stiffness matrix R k h of (44).

Step 2 Solve the following least-squares problem:

min K∈R N ×N M k h K -R m k,h 2
F . (59)
The solution is then given by

K m k,h := M k h -1 R m k,h . (60)

Convergence

We will now prove that K m k,h converges to the matrix K k h given by [START_REF] Igor | Spectrum of the Koopman Operator, Spectral Expansions in Functional Spaces, and State-Space Geometry[END_REF]. Now, as we have considered random points (x l ) m l=1 , this means that the coefficients of

K k h -K m k,h
are random variables. The notion of convergence of random variables we consider here is then the notion of convergence in distribution:

Definition 4.1 (Convergence in distribution of a sequence of random variables). A sequence (X n ) n∈N of real random variables is said to converge in distribution to a real random variable X:

X n D ----→ n→∞ X
if their cumulative distribution functions converge pointwise to the cumulative distribution function of X on its continuity set.

We can now state the following convergence result (see also [START_REF] Williams | Extending Data-Driven Koopman Analysis to Actuated Systems[END_REF][START_REF] Korda | On Convergence of Extended Dynamic Mode Decomposition to the Koopman Operator[END_REF][START_REF] Koltai | On the numerical approximation of the perronfrobenius and koopman operator[END_REF]), which focuses on the convergence of the coefficients of the approximate stiffness matrix R m k,h : Proposition 4.1 (Convergence and complexity of gEDMD). The computational complexity of gEDMD has the upper bound O(mN 2 + N 3 ).

For f ∈ L ∞ (Ω, R d ), the following convergence in distribution holds:

√ m R m k,h ij -R k h ij D ----→ m→∞ N 0, σ k,h ij , (61)
where

σ k,h ij 2 = Ω ψ 2 i (x) (f (x) • ∇ψ j (x)) 2 dx - Ω ψ i (x)f (x) • ∇ψ j (x)dx 2 , ( 62
)
and N (0, σ k,h ij ) is the normal distribution with mean 0 and standard deviation σ k,h ij . Proof. The upper bound for the complexity is given by the complexity of the linear regression that is performed to compute K m k,h . It is well-known (see [START_REF] Tan | Introduction to Data Mining[END_REF]Appendix D]) that this complexity is bounded by O(mN 2 + N 3 ) due to the matrix product operation, and the pseudo-inversion that follows.

As for the convergence in distribution, notice that the coefficients of R m k,h are given by

R m k,h ij = 1 m m l=1 ψ i (x l )f (x l ) • ∇ψ j (x l ),
which is none other than the classical Monte-Carlo approximation of the integral

ψ i , f • ∇ψ j = Ω ψ i (x)f (x) • ∇ψ j (x)dx = R k h ij ,
as the sampling points x l are drawn independently and uniformly with respect to the Lebesgue measure. Convergence [START_REF] Matthew | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF] then follows directly from a classical result of Monte-Carlo integration (see [START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF]Theorem 2.1]:

√ m 1 m m l=1 ψ i (x l )f (x l ) • ∇ψ j (x l ) - Ω ψ i (x)f (x) • ∇ψ j (x)dx D ----→ m→∞ N (0, σ k,h ij ). ( 63
)
Remark 4.1. In the usual implementation of gEDMD, the mass matrix is not computed offline, but also approximated with the Monte-Carlo approximation

M m k,h := 1 m G m k,h G m k,h . The matrix K m k,h is then given by K m k,h = M m k,h † R m k,h , which converges to K k
h provided the pseudo-inverses converge. By a classical result due to R.Penrose [START_REF] Penrose | A generalized inverse for matrices[END_REF], this is the case if, for m large enough, the M m k,h are invertible with probability 1. This corresponds to the following condition:

ϕ(x 1 ) = • • • = ϕ(x m ) = 0 =⇒ ϕ = 0, ∀ϕ ∈ V k h . ( 64
)
When the (x j ) j are drawn uniformly, this is equivalent to

∀ϕ ∈ V N \ {0}, µ{x ∈ Ω, ϕ(x) = 0} = 0, (65)
which is not satisfied by the shape functions (ψ j ) j , nor by other classical localized approximation functions such as wavelets.

In our case, to ensure (64), it would be necessary to draw sample points in every element of the rectangulation T h , which would be extremely impractical.

Convergence in matrix norms

Proposition 4.1 essentially means that, for a given confidence rate α ∈ [0, 1), there exists a constant C k,h ij (α) given by the normal distribution N (0, σ k,h ij ), such that for large enough m the following estimate holds with probability at least α:

R m k,h ij -R k h ij ≤ C k,h ij (α)m -1 2 . ( 66
)
This probabilistic estimate tells us that the coefficients converge in m -1 2 , with a constant which depends on the standard deviation σ k,h ij , that is, recalling its expression [START_REF] Xiao | A deep learning framework based on koopman operator for data-driven modeling of vehicle dynamics[END_REF], on the shape functions (ψ j ) and f . We then note

C k h (f, α) := max ij (C k,h ij (α)). ( 67
)
Now as we work with L 2 functions, we need to relate this matrix norm to the operator norm on L 2 (Ω).

We define the following norm on R N k h :

c 2 L 2 = N k h j=1 c j ψ j 2 = c M k h c = (M k h ) 1 2 c 2 . ( 68
)
The matrix M k h is symmetric definite positive, we note ρ + k,h and ρ - k,h its largest and smallest eigenvalue respectively. We then have

ρ - k,h c ≤ c L 2 ≤ ρ + k,h c , ∀c ∈ R N k h . (69)
With this norm we then define the induced operator norm on

R N k h ×N k h : |||A||| L 2 = sup c L 2 =1 Ac L 2 . ( 70
)
Applying (69), we then have

ρ - k,h ρ + k,h |||A||| L 2 ≤ |||A||| ≤ ρ + k,h ρ - k,h |||A||| L 2 , ∀A ∈ R N k h ×N k h . ( 71
)
We can now state the following proposition:

Proposition 4.2 (Convergence of gEDMD in matrix norm). Let α ∈ [0, 1). For f ∈ L ∞ (Ω, R d ), there exists a constant C k h (f, α
) such that, for m large enough and with probability at least α,

|||K k h -K m k,h ||| L 2 ≤ N k h ρ - k,h ρ + k,h ρ - k,h C k h (f, α)m -1 2 . ( 72
)
Proof. From the probabilistic estimate (66) we can then deduce the following estimate on the max matrix norm: for m large enough, with probability at least α,

R m k,h -R k h ∞ ≤ C k h (f, α)m -1 2 (73)
Now, we recall the result of [START_REF] Tonge | Equivalence constants for matrix norms: a problem of goldberg[END_REF] on sharp equivalence constants for matrix norms:

1 N k h |||A||| ≤ A ∞ ≤ |||A|||, ∀A ∈ R N k h ×N k h . ( 74
)
Putting ( 74), ( 71) and ( 73) together, we get, for m large enough and with probability at least α:

|||R k h -R m k,h ||| L 2 ≤ N k h ρ + k,h ρ - k,h C k h (f, α)m -1 2 (75)
Finally, by submultiplicativity of operator norms, and the expression of K m k,h , we have

|||K k h -K m k,h ||| L 2 = |||(M k h ) -1 (R k h -R m k,h )||| L 2 ≤ |||(M k h ) -1 ||| L 2 |||R k h -R m k,h ||| L 2 . ( 76
)
Now, as

|||(M k h ) -1 ||| L 2 = |||(M k h ) -1 ||| = 1 ρ - k,h , (77)
we finally deduce (72) from ( 75).

This shows that the convergence of K k h to K m k,h , i.e. the convergence of gEDMD, is asymptotically in m -1 2 , which is quite slow but does not depend on the underlying dimension d of the problem (4). It is however further slowed down by numerous factors: the standard deviation σ k,h ij , which represents the oscillatory behavior of the shape functions, the conditioning number and smallest eigenvalue of the mass matrix M k h , and the dimension of the finite elements space V k h .

5 Quantitative analysis

Using gEDMD for system identification and predictions: convergence analysis

We now use the matrix K m k,h computed by gEDMD to give an approximation of f and its flow Φ t f . Recalling [START_REF] Isakov | Inverse problems for partial differential equations[END_REF] and [START_REF] Mauroy | Global Stability Analysis Using the Eigenfunctions of the Koopman Operator[END_REF], an approximation of the f i is given by

K m k,h g i , i ∈ 1, d . (78)
Recalling ( 26), ( 27), ( 43), ( 45) and ( 60), an approximation of the Φ t f i is given by the solutions to

   d dt ϕ k,h i,m = K m k,h ϕ k,h i,m , ϕ k,h i,m ∈ V k h , ϕ k,h i,m (0, •) = g i (•), i ∈ 1, d . (79)
In the basis (ψ j ), noting

g i = N k h j=1 γ i j ψ j , the components of the K m k,h g i are given by K m k,h γ i ∈ R N k h , i ∈ 1, d (80)
and the components of the ϕ k,h i,m (t, •) are given by the solutions to

   d dt ϕ k,h i,m = K m k,h ϕ k,h i,m , ϕ k,h i,m ∈ R N k h , ϕ k,h i,m (0, •) = γ i , i ∈ 1, d . (81)
We now give quantitative estimates for system identification and predictions, using gEDMD with finite elements spaces: Proposition 5.1. Let f ∈ C k+1 (Ω). There exists constants C 1 , C 2 (f ), C 3 > 0 such that, for a given probability rate α ∈ [0, 1], one has for m large enough and with probability at least α:

K m k,h g i -f i ≤ N k h ρ - k,h ρ + k,h ρ - k,h C k h (f, α)m -1 2 g i + C 1 h k+1 |f i | k+1 , ϕ k,h i,m (t, •) -Φ t f i ≤ C 3 te t|||K k h ||| L 2 N k h ρ - k,h ρ + k,h ρ - k,h C k h (f, α)m -1 2 g i + C 2 (f )h k , ∀i ∈ 1, d . (82)
Proof. Recalling ( 23) and ( 42), we have the following:

for i ∈ 1, d , K m k,h g i -f i ≤ K m k,h g i -K k h g i + K k h g i -f i ≤ |||K m k,h -K k h ||| L 2 g i + K k h g i -f i . (83)
Applying ( 72) and ( 46) yields the first inequality of (82). Recalling ( 27) and ( 43), we have the following:

for i ∈ 1, d , ϕ k,h i,m (t, •) -Φ t f i ≤ ϕ k,h i,m (t, •) -ϕ k,h i (t, •) + ϕ k,h i (t, •) -ϕ i (t, •) ≤ |||e tK m k,h -e tK k h ||| L 2 g i + Φ t f i -ϕ k,h i (t, •) . (84)
By property of the matrix exponential,

|||e tK m k,h -e tK k h ||| L 2 ≤ t|||K m k,h -K k h ||| L 2 e t|||K k h ||| L 2 e t|||K m k,h -K k h ||| L 2 (85)
Combining ( 85) with ( 84) and (72) then yields the second inequality of (82).

With these inequalities it is now clear that the "slowest" term is the data-driven approximation error term in m -1 2 . It depends both on the number of samples m and the dimension of the Galerkin projection space V k h . More precisely, for given k and h, we see from the first inequality of (82) that for the data-driven approximation error to reach the same accuracy as the finite elements approximation error, m has to satisfy

m ≥ N k h ρ - k,h 2 ρ + k,h ρ - k,h C k h (f, α) 2 h -2k-2 . (86)
For a fixed k, choosing more accurate finite elements spaces by decreasing the mesh size h has a double impact on the required number of samples m: on the h -2k-2 term, on one hand, and on the quantity

(N k h /ρ - k,h ) 2 (ρ + k,h /ρ - k,h )C k h (f, α) 2
on the other hand. Moreover, due to the presence of N k h we see that this number of samples must also increase exponentially with d. Thus, indirectly, the convergence of gEDMD is also subject to the curse of dimensionality.

Overall, the slowness of convergence of classical Monte-Carlo approximations usually employed in gEDMD conceals the high accuracy one can obtain from classical approximation spaces, which in a way defeats the purpose of choosing good approximation spaces.

In the next section we will focus on the case d = 1 and fix k = 1, in order to give precise estimates of the terms in (86). This will in turn allow us to give a precise estimate of the computational cost of gEDMD and give a clear comparison to direct interpolation methods.

Comparison with direct interpolation in 1-D for linear finite elements

In this subsection, d = 1 and k = 1. That is, we focus on gEDMD for systems on the interval [0, L], using Galerkin projections on linear finite elements.

System identification using the Koopman operator appears as a quite indirect method, as it consists in observing trajectories of the differential equation, working with the associated transport equation, then coming back to the differential equation.

In the gEDMD algorithm, the values f (x l ) at the sample points are approximated. Instead of using them to compute an approximation of the Koopman operator, one could also use them to interpolate f directly, without resorting to Koopman operator methods. In that case, predictions are then made by integrating the nonlinear ODE

ẋ = f (x), x ∈ R d ,
where f is the interpolant of f , instead of a Galerkin method of the linear transport equation such as (81). We will focus here on linear and cubic spline interpolation.

In theory, the linear interpolant (f, x i ) and the natural cubic spline interpolant σ(f, x i ) (with derivative 0 at the endpoints) satisfy the following inequalities

(f, x i ) -f ≤ C 1 δ 2 |f | 2 , f ∈ H 2 (0, L), σ(f, x i ) -f ≤ C 2 δ 4 |f | 4 , f ∈ H 4 (0, L), (87)
where

δ = max{x n+1 -x n , x 1 + L, L -x m , n ∈ 1, m } (88)
The inequality for the linear interpolant follows from classical inequalities on C 2 functions, for the cubic spline interpolation error bounds, we refer to [START_REF] Hall | Optimal error bounds for cubic spline interpolation[END_REF].

The m samples are drawn in the whole interval, and for large m we have (see [START_REF] Pyke | Spacings. (With discussion.)[END_REF])

δ ∝ log(m) m so that (3) becomes (f, x i ) -f ≤ C 2 log(m) m 2 |f | 2 , f ∈ H 2 (0, L), σ(f, x i ) -f ≤ C 2 log(m) m 4 |f | 4 , f ∈ H 4 (0, L). (89)
Comparing with the m -1 2 rate of convergence of gEDMD, we get

Hence, for a given number of samples m, gEDMD is less accurate than direct interpolation methods.

To give a more thorough comparison, let ε > 0. We now compare the number m of samples required in each method in order to to achieve an accuracy of at most ε.

For linear interpolation, the required number of samples satisfies:

log(m) m ∼ ε 1 2 =⇒ ε -1 2 ≤ m ≤ Cε -1 2 + , ∀ > 0.
For cubic spline interpolation, the required number of samples satisfies:

log(m) m ∼ ε 1 4 =⇒ ε -1 4 ≤ m ≤ Cε -1 4 + , ∀ > 0.
For modified gEDMD, first consider the Galerkin projection. From [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF], h must satisfy:

h ≤ Cε 1 2 . (91)
Recalling (86):

m ≥ C 1 h (f, α) 1 + 1 h ρ - 1,h 2 c 1 h h -4
and using the estimates (110), (111), and (120) from Appendix A, one has:

m ≥ C C(α) f L ∞ h -2 1 + 1 h h 2 ε -2 ≥ C C(α) 2 f 2 L ∞ h -8 ε -2 ≥ C C(α) 2 f 2 L ∞ ε -6 (92)
Now, it is well-known that interpolation methods have linear complexity (see [START_REF] Toraichi | Computational complexity of spline interpolation[END_REF]). On the other hand, we have seen in Proposition 4.1 that the computational complexity of gEDMD is in O(mN 2 +N 3 ). Regarding N , (91) implies

N ≥ C ε -1 2 .
Together with (92), this yields the following complexity in terms of ε:

O ε -7 + ε -3 2 = O(ε -7 )
. We summarize the above in the following table:

Parameters Cost Koopman h operator k ε -7 method m Linear m ε -1 2 + interpolation Splines m ε -1 4 + interpolation
This clearly shows that interpolation methods are more cost-effective than gEDMD with linear finite elements.

Numerical illustration

We show the results of some implementations of gEDMD, focusing on the approximation of the vector field f . Given our remarks on the curse of dimensionality in Section 3.6, we focus on the case d = 1. We have given precise estimates in the previous section that show how the convergence of the data-driven approximation is affected by the choice of the Galerkin projection. Given [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF] and (86), this problem worsens in higher dimensions.

We take f given by combinations of trigonometric functions, polynomials, and exponentials on the interval [0, 1]. We implement gEDMD with linear finite elements to recover an approximation of f . We compare the results with linear interpolation and spline interpolation on the same data sample.

Linear finite elements

Data sample we form the data set:

{x n , f (x n ), n ∈ 1, m } (93)
To simplify matters we directly take the values of the nonlinearity f we have chosen, at the sample points x n . In practice, f is unknown, so one obtains approximate values of f at the sample points by observing the trajectories of the system with initial conditions x n , and approximating the initial velocities by finite differences (see Section 4 and [15]).

gEDMD. Compute the following matrices from the data:

G m 1,h := (ψ i (x j )) i∈ 1,N h j∈ 1,m A m 1,h := ( ẋj • ∇ψ i (x j )) i∈ 1,N h j∈ 1,m , (94)
and define

R m 1,h = 1 m G m 1,h A m 1,h
.

Perform the linear least squares regression:

K m 1,h := argmin A∈R N ×N M 1 h A -R m 1,h 2 = M m 1,h -1 R m 1,h . (95)
Identifying the nonlinearity in the differential equation. Now consider the identity function g(x) = x: g ∈ P 1 as it is linear, and

g(x) = M j=0 jhψ j (x) = c ψ(x), ∀x ∈ [0, L]. (96)
To recover f , we compute

c := c K m 1,h (97)
to get the following approximation of f on the interval [0, L]:

f m 1,h (x) = c ψ(x) = M j=0 cj ψ j (x), ∀x ∈ [0, L]. (98)

Numerical simulations

Figure 1: Identification of different functions using gEDMD with linear finite elements.

Figure 1 shows approximations of f (given by the thicker black curve) computed with gEDMD for a fixed dimension of the finite elements space, and different numbers of samples m. One can see that for a small number of samples, there are sizeable discrepancies. For a larger number of samples the curves obtained by gEDMD seem to match the curve rather well.

In Figure 2, we plot the relative approximation error

f m N -f f
as a function of the number of samples m, for several values of N : as expected, increasing the number of samples decreases the error. This decrease appears to be slower for greater values of N (see the red dots in the figure below), which is consistent with (86). Finally, in Figure 3 we compare the relative approximation error of gEDMD with that of direct interpolation methods. We implement gEDMD for increasing values of N , with m = 100N samples. We compare the results with linear interpolation and spline interpolation with m random samples. There is a clear hierarchy between the three methods, gEDMD being the least accurate by far.

Conclusion

In this article, we have presented the implementation of a system identification method, based on finite-dimensional, data-driven approximations K m N of the so-called Koopman operator, which is the transport operator associated to the unknown vector field f :

f (x) • ∇.
From trajectory data on the system ẋ = f (x), measured at randomly drawn sample points (x n ) n∈ 1,m one can compute the sample values

f (x n ), x n ∈ Ω, n ∈ 1, m .
From this data set, the method provides an approximation of the vector field f and its flow Φ t f , using the coordinate functions g i (x) = x i :

f i ≈ K m N g i , Φ t f i ≈ e tK m N g i .
The potential advantage of this method lies in its linearity of the Koopman operator: as can be seen above, the nonlinear flow Φ t f is approximated thanks to the linear flow e tK m N , which is however defined on a subspace of functions rather than the state space Ω ⊂ R d .

We have focused on the case of Galerkin projections of the Koopman operator on finite elements spaces. They are generic approximation spaces, widely used in numerical analysis to approximate PDEs. In this case, the approximating properties of these spaces guarantee that the corresponding Galerkin projection of the Koopman operator is a reliable approximation, without having to rely on hypothetical spectral properties. Indeed the Koopman operator (which is essentially a transport operator) does not have nice spectral properties in general.

For 1-dimensional problems, we have compared this method with direct interpolation of the vector field f from the sample values f

(x n ), x n ∈ Ω, n ∈ 1, m .
Interpolation yields an approximation f of f . Predictions can then be made by integrating the nonlinear differential equation ẋ = f (x).

We have found that this Koopman operator based system identification method is generally not advantageous in practice, due to two drawbacks.

Dimension of the Galerkin projections With Galerkin projections on classical approximation

spaces, high accuracy requires high-dimensional spaces of functions. This drawback, due to the well-known curse of dimensionality, is a recurring feature of approximation spaces. On these spaces, the approximation of the linear transport equation associated to the Koopman operator then comes at a high computational cost.

2. Data-driven approximation of the Galerkin projections Galerkin projection on finite elements spaces provide accurate, albeit high-dimensional, approximations of the Koopman operator for system identification and predictions. However, for the system identification problem we have formulated, these projections must be recovered from a finite amount of data. To that effect, we have brought a novel modification to the usual gEDMD algorithm. As we are working with classical finite elements spaces, the mass matrix can be computed offline. Then, recovering the Galerkin projection of the Koopman operator reduces to approximating the stiffness matrix from data. To approximate the integrals in the stiffness matrix, given the nature of the data set, the natural path is to use Monte-Carlo integration (as is the case for regular EDMD in the literature).

In addition to the inherent slowness of Monte-Carlo integration, we furthermore observe that the rate of convergence of the data-driven approximation has multiple dependencies in the parameters h (mesh size) and k (degree) of the finite elements space. As a consequence, the more accurate Galerkin projection we take, the harder it becomes to approximate from data by the gEDMD method.

We have quantified this precisely in 1-D, which allows for a clear comparison between recovery of the vector field f by the gEDMD method and direct interpolation (linear or cubic splines) of f .

We have based our analysis on a specific choice of approximation space, but our estimates illustrate what is probably a general phenomenon: the Galerkin projection of the Koopman operator on generic approximation spaces is prohibitively slow to approximate from data.

Let us conclude with some prospects that could be explored to overcome these drawbacks:

Galerkin projections of low dimension Considering that classical interpolation methods (or SINDy with sparsity constraints) provide a very good approximation of the vector field, the goal for the gEDMD method is to provide a numerically interesting approximation of the flow Φ t f in order to make predictions. This will be the case if adequate lower-dimensional Galerkin projections can be found, which still allow for reliable estimates as in Proposition 3.2 with finite elements spaces. Finding such spaces would also reduce the impact of the choice of the Galerkin projection on the Monte-Carlo convergence, which in addition to being inherently slow, is further slowed down when taking higher-dimensional finite elements spaces.

It seems very likely that there are many situations where this is not possible (chaotic systems are a frequently cited example [START_REF] Arbabi | Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator[END_REF][START_REF] Arbabi | Generative stochastic modeling of strongly nonlinear flows with non-gaussian statistics[END_REF]). Between the universal solution of classical approximation spaces, and the ideal situation of subspaces of functions that are invariant under the Koopman operator (see [START_REF] Steven | Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control[END_REF]), the complete picture is still missing: what characterizes the existence of these relevant low-dimensional Galerkin projections? How "rare" is it? Is there a data-driven way to determine if they exist and find them?

Use of deep neural networks

There have been many recent attempts to find such low-dimensional Galerkin projections, using deep neural networks [START_REF] Ping | Deep Koopman model predictive control for enhancing transient stability in power grids[END_REF][START_REF] Yeung | Learning deep neural network representations for koopman operators of nonlinear dynamical systems[END_REF][START_REF] Xiao | A deep learning framework based on koopman operator for data-driven modeling of vehicle dynamics[END_REF][START_REF] Korda | Power grid transient stabilization using koopman model predictive control[END_REF][START_REF] Takeishi | Learning koopman invariant subspaces for dynamic mode decomposition[END_REF]. The output of such networks is a family of functions ψ 1 , • • • , ψ N , where N is fixed beforehand. These functions are taken from a large class of functions (essentially given by the architecture of the network) and the network is trained to minimize the prediction error made when performing EDMD with the subspace V N = Span{ψ 1 , • • • , ψ N }, while additionally penalizing situations where V N does not contain the g i , for the purpose of system identification (see Section 5.1). Although numerical implementations yield promising results, there is still much work to be done to interpret the outputs of these deep neural networks. Moreover, as is mentioned in [START_REF] Takeishi | Learning koopman invariant subspaces for dynamic mode decomposition[END_REF], there is no guarantee that the output is not a local minimizer. We note

N := dim V 1 h = M + 1 = 1 h + 1. (102)
We denote by ψ j , j ∈ 0, M the node functions, which satisfy

ψ j (kh) = δ jk , j, k 0, M , (103)
they form a basis of V 1 h . Explicitly, for a given node jh, ψ j is defined by

j ∈ 1, M -1 , ψ j (x) =                  0 on [-L, (j -1)h] ∪ [(j + 1)h, L]
x -(j -1)h h on [(j -1)h, jh] (j + 1)h -x h on [jh, (j + 1)h] (104)

Ψ 0 (x) =        h -x h on [0, h] 0 on [h, L] (105)
Ψ M (x) =        0 on [-L, L -h] x -L + h h on [L -h, L] (106)
Classically, the mass matrix is given by:

M 1 h = h                1/3 1/6 0 • • • • • • • • • 0 1/6 2/
• • • • • • • • • 0 1/6 1/3                (107)
and its eigenvalues are given by (after a slight modification of [1, Exercise 7.4.1]):

λ j = h 3 2 + cos jhπ 1 + 2h , j ∈ 1, N . (108)
Consequently,

ρ + 1,h = h 3 2 + cos πh 1 + 2h , (109)
Let us now give estimates on the variances of the coefficients of the rigidity matrix, which is involved in the convergence of the method in Section 4.2. First, an obvious remark is that

σ k,h ij = 0, |i -j| > 1. (112)
On the other hand, assuming that f ∈ C 2 (Ω): for i ∈ 0, M -1 ,

Ω ψ i (x)f (x)ψ i+1 (x) 2 dx - f (x i ) 2 h 2 xi+1 xi ψ i (x) 2 dx = 1 h 2 xi+1 xi ψ i (x) 2 (f (x) 2 -f (x i ) 2 )dx ≤ f ∞ C 1 h , (113)
for some constant C 1 > 0. Similarly, Ω (ψ i+1 (x)f (x)ψ i (x))

2 dx-f

(x i ) 2 h 2 xi+1 xi ψ i+1 (x) 2 dx = 1 h 2 xi+1 xi ψ i+1 (x) 2 (f (x) 2 -f (x i ) 2 )dx ≤ f ∞ C 2 h , (114) and
Ω ψ i (x)f (x)ψ i+1 (x)dx 2 - f (x i ) 2 h 2 xi+1 xi ψ i (x)dx 2 ≤ f ∞ C 3 h , ( 115
) Ω ψ i+1 (x)f (x)ψ i (x)dx 2 - f (x i ) 2 h 2 xi+1 xi ψ i+1 (x)dx 2 ≤ f ∞ C 4 h , (116)
for some constants C 2 , C 3 , C 4 > 0.

Putting (113) and (115) together, we get, for i ∈ 0, M -1 ,

σ k,h i,i+1 2 = 1 3h 2 f (x i ) 2 + O 1 h , (117)
and similarly with (114) and (116),

σ k,h i+1,i 2 = 1 3h 2 f (x i ) 2 + O 1 h , (118)
and, for i ∈ 0, M ,

σ k,h ii 2 = 2 3h 2 f (x i ) 2 + O 1 h . (119)
Recalling (67), this yields

C 1 h (f, α) = 2 3h 2 f 2 ∞ + O 1 h . (120)

Figure 2 :

 2 Figure 2: Relative error of gEDMD with linear elements for different values of N . Dependence on the number m of samples.

Figure 3 :

 3 Figure 3: Comparison with direct interpolation methods. m = 100(N -1).23

Figure 4 :

 4 Figure 4: The shape function ψ j

The "Koopman lifting method", proposed in[START_REF] Mauroy | Koopman-based lifting techniques for nonlinear systems identification[END_REF], where the authors apply EDMD to a continuous time system to obtain an approximation of the Koopman semigroup U (τ ) for a sampling time τ > 0. Then, they compute its matrix logarithm to produce an approximation of the infinitesimal generator K.

2. Generator Extended Dynamic Mode Decomposition (gEDMD), proposed in[START_REF] Klus | Data-driven approximation of the koopman generator: Model reduction, system identification, and control[END_REF], where EDMD is adapted to directly produce an approximation of K by processing data in different way.

Acknowledgements We would like to thank Pr.Lars Grüne for fruitful discussions on this topic.

We are greatly indebted to Jesus Oroya Villalta, for his help on the numerical experiments and his astute comments.

A Linear 1-D finite elements

We fix a mesh size

we define the mesh


N Loiseau-Witon

R Kéchichian

S Valette

A Bartoli
email: adrien.bartoli@gmail.com

Learning 3D Medical Image Patch Descriptors with the Triplet Loss

Keywords: Medical imaging, Triplet loss, Convolution Neural Network

Purpose

Computational anatomy focuses on the analysis of the human anatomical variability. Typical applications are the discovery of differences across healthy and sick subjects and the classification of anomalies. A fundamental tool in computational anatomy, which forms the central focus of this paper, is the computation of point correspondences across volumes (3D images) such as Computed Tomography (CT) volumes, for multiple subjects. More specifically, we consider automatically detected keypoints and their local descriptors, computed from the image or volume patch surrounding each keypoint. Theses descriptors are essential because they must be discriminant and repeatable [START_REF] Bay | Surf: Speeded up robust features[END_REF][START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. Learned descriptors based on Convolutional Neural Networks (CNN) have recently shown great success for 2D images [START_REF] Balntas | Learning local feature descriptors with triplets and shallow convolutional neural networks[END_REF]. However, while classical 2D image descriptors were extended to volumes [START_REF] Agier | Hubless 3d medical image bundle registration[END_REF], recent learning-based approaches have been limited to 2D detection and description. The extension to 3D descriptors was only proposed in [START_REF] Blendowski | 3d-cnns for deep binary descriptor learning in medical volume data[END_REF], in the context of image retrieval. We propose a methodology to construct these learned volume keypoint descriptors. The main difficulty is to define a sound training approach, combining a training dataset and a loss function. In short, we propose to generate semi-synthetic data by transforming real volumes and to use a triplet loss inspired by 2D descriptor learning. Our experimental results show that our learned descriptor outperforms the hand-crafted descriptor 3D-SURF [START_REF] Agier | Hubless 3d medical image bundle registration[END_REF], a 3D extension of SURF, with similar runtime.

Methods

Our first goal is to create a reference dataset defining keypoint correspondences between multiple volumes. In 2D, these correspondences can be established using Structure-from-Motion [START_REF] Altwaijry | Learning to detect and match keypoints with deep architectures[END_REF]. In 3D medical images, keypoints could be defined as anatomical landmarks placed by medical experts. However, no such large annotated dataset is publicly available. We thus propose to create a semi-synthetic reference dataset by transforming real volumes.

Constructing a semi-synthetic dataset

We use two subsets from the Visceral dataset [START_REF] Langs | Visceral: towards large data in medical imaging-challenges and directions[END_REF]. The first subset, named Gold, contains 20 CT volumes, each annotated with about 40 landmarks. The second subset, Silver, contains 60 CT volumes without landmarks.

In order to generate new volumes, we estimate the probability density of possible inter-volume transformations, and sample from this density to warp CT volumes and define keypoint correspondences. We first compute inter-volume local affine transformations in the least-squares sense. For each landmark in each volume of the Gold subset, we use the landmark, its three closest landmarks and the four corresponding landmarks in another volume to estimate a local transformation. When the landmark and its three neighbours are almost collinear (e.g. vertebral landmarks), the least-squares problem is ill-conditioned and we therefore discard the transformation. Thus we obtain at most Lk(k -1)/2 affine transformation matrices of size 4 × 4, where L = 40 is the number of landmarks and k = 20 is the number of volumes. The Pearson test shows that the elements of these matrices are independent, allowing us to sample each element from its distribution independently. We apply Kernel Density Estimation (KDE) to these matrices to estimate the density of inter-volume transformations. Student's ttest shows that a Gaussian Kernel is a good fit for the KDE. We estimate the kernel bandwidth via Scott's rules. To generate our semi-synthetic dataset, we sample transformations from this density and apply them to volumes in the Silver subset. More specifically, for a sampled transformation t and a silver volume V i , we obtain the volume V t i . We detect the keypoints in V i and V t i using 3D-SURF and obtain two keypoint sets P i and P t i . Note that 3D-SURF is both a detector and a descriptor. We then apply the inverse transform t -1 to the keypoints from V t i . Finally we use a k-d tree to construct the set of corresponding keypoints between the volumes as the set of pairs: (p ∈ P i , q ∈ t -1 (P t i )), whose p to q distance is lower than 8 mm. This threshold was chosen to obtain a large number of correct correspondences.

Training the descriptor with the triplet loss

We learn a descriptor CNN mapping a 3D patch of 10 3 voxels surrounding a keypoint to a descriptor vector. Recent work in 2D has shown that learning descriptors using triplets yields better results than using pairs [START_REF] Hoffer | Deep metric learning using triplet network[END_REF]. Triplet learning requires forming triplets of patches {a, p, n} where a is an anchor, p a positive representing a different patch of the same class as a, and n a negative representing a patch of a different class. In our case a and p are two patches around corresponding keypoints from different volumes and n is a patch around a different keypoint. The aim is to optimize the CNN parameters in order to bring a and p close together in descriptor space and to push n away from a. Thus the triplet loss is defined by L(a, p, n) = max( f (a)-f (p) 2 -f (a)-f (n) 2 +α, 0), where f (•) is the CNN and α the margin parameter.

Our CNN is defined by two 3D convolution layers, one maxpooling layer between the convolutions, and a fully-connected layer which gives the final descriptor. The network architecture is illustrated by Figure 1. Passing a patch of size 10 3 through this CNN gives us a descriptor vector of the desired size. We use a size of 48 for direct comparability with 3D-SURF.

Results

Data split. We divide the Silver dataset into training and validation subsets.

The training data consist of 55 subjects with 10 transformed volumes each, following our procedure of semi-synthetic data generation. The validation data consist of 5 subjects with 10 transformed volumes each. For testing, we use the Gold subset with 20 subjects and the associated anatomical landmarks. Training. Optimization is performed via Stochastic Gradient Descent, with a batch size of 1000 patches, a learning rate of 0.1, a momentum of 0.9, a weight decay of 10 -6 and a loss margin of 0.2. We also use online triplet mining to find the best triplets for learning. Our CPU-based implementation uses the PyTorch library. The training of a single epoch with 10 6 triplets takes about 30 minutes and approximately 10 GB of memory on a Linux 64-bit platform running on an Intel Xeon 2.6 GHz CPU. Our model is light enough to be trained on the CPU; using the GPU did not significantly reduce training time. Evaluation. We evaluate the descriptor using two different metrics. The first metric is the false positive rate at point 0.95 of true positive recall (FPR95) [START_REF] Brown | Discriminative learning of local image descriptors[END_REF]. We compute FPR95 based on descriptor distances of randomly selected 10 5 keypoint pairs with 50% corresponding and 50% non corresponding pairs. A low FPR95 indicates good results. The second metric is the mean landmark distance calculated on ground-truth landmarks in Gold volumes after registering them to a common space using the keypoint-based FROG registration algorithm [START_REF] Agier | Hubless keypoint-based 3d deformable groupwise registration[END_REF]. For comparison, we replace the 3D-SURF descriptor used in this algorithm with our learned descriptor. Low mean distance between landmarks indicates good results. Table 1 shows that the proposed descriptor yields better results in terms of both FPR95 and mean landmark distance compared to the 3D-SURF descriptor.

Fig. 1 .

 1 Fig. 1. Architecture of our CNN and the triplet-loss. The input patches are 10 3 voxels. The CNN has two convolutions with tanh activation, one maxpooling and a fullyconnected layer. The triplet {a, p, n} is passed through the network, and descriptor vectors are fed to the triplet loss.

Table 1 .

 1 Performance comparison of the 3D-SURF and our learned descriptors.4 Conclusions and future workOur results, although preliminary, show that a learned 3D descriptor, trained on semi-synthetic data, can outperform a carefully hand-crafted one. We intend to further explore these promising results by extending our training dataset and conducting more experiments. Future research will address training a 3D keypoint detector.

	Type of descriptor	FPR95	Mean landmark distance
	3D-SURF	0.077	8.74
	Learned	0.022	8.54

Acknowledgements

This work was funded by the TOPACS ANR-19-CE45-0015 project of the French National Research Agency (ANR).