Linus Bleistein
email: <li-nus.bleistein@inria.fr>.

Adeline Fermanian

Anne-Sophie Jannot

Agathe Guilloux

Learning the Dynamics of Sparsely Observed Interacting Systems

We address the problem of learning the dynamics of an unknown non-parametric system linking a target and a feature time series. The feature time series is measured on a sparse and irregular grid, while we have access to only a few points of the target time series. Once learned, we can use these dynamics to predict values of the target from the previous values of the feature time series. We frame this task as learning the solution map of a controlled differential equation (CDE). By leveraging the rich theory of signatures, we are able to cast this non-linear problem as a highdimensional linear regression. We provide an oracle bound on the prediction error which exhibits explicit dependencies on the individual-specific sampling schemes. Our theoretical results are illustrated by simulations which show that our method outperforms existing algorithms for recovering the full time series while being computationally cheap. We conclude by demonstrating its potential on real-world epidemiological data.

Introduction

Time series are ubiquitous in many areas such as finance, economics, robotics, agriculture, and healthcare. One is typically interested in modelling the evolution of a target quantity through time, which is known to be affected by a set of time-evolving features. For example, pollution levels in a city are driven by quantities such as temperature, pressure, traffic, or economic activity measured through time. Mathematically, one wishes to model the evolution of a quantity y t ∈ R p , p ≥ 1, as a function of some time evolving features x t ∈ R d , d ≥ 1, for t ∈ [0, 1]. In other words, the goal is to learn the dynamics that link the target to the features. Such an interaction is typically modelled via differential equations, which are a common choice of model in natural sciences [START_REF] Zwillinger | Handbook of differential equations[END_REF]. In this article, we assume that there exists a function G : R p × R d → R p such that

y t = y 0 + t 0 G(y s , x s )ds (1)
or equivalently

dy t = G(y t , x t )dt, y 0 ∈ R p .
The value y t depends on the trajectory of the features time series x s up to time t. Learning the dynamics of the system can be framed as learning the solution map of (1), i.e., a function Ψ which, given a time t, an initial point y 0 ∈ R p , and the history of the path up to time t, denoted by x [0,t] = (x s ) s∈[0,t] , outputs the value of y at time t.

If we know Ψ, we gain access to the values of y at any point in time provided we know the values of x up to this point ; this encompasses many tasks such as forecasting or interpolating between points of y. We specifically have in mind applications where we have an easy access to x but a limited one to y.

This problem is extremely common in healthcare. For example, in obstetrics, the lactic acidosis (LA) of the fetus, which is a proxy for fetal distress, is a quantity of high medical interest for predicting complications in the first hours after birth. This biomarker cannot be measured during pregnancy but only at birth because the measurement is highly invasive. Some vitals such as heart rate and fetal movement are however easy to measure during pregnancy. In this case, x is the non-invasive measurements made during pregnancy, while y is the invasive measurement of LA at birth. Predicting the value of y at any time t (both before and at birth) would allow for early diagnosis. Similarly, after surgery, patients are often monitored to detect hemorrhage. While some vitals such as heart rate of saturation are monitored in continuous time, haemoglobin-which is highly predictive of hemorrhage-is only measured by 1 arXiv:2301.11647v2 [stat.ML] 31 May 2023 blood samples taken a few times a day, which can significantly delay hemorrhage diagnostic.

Irregular data. In practice, the functions x and y are measured on discrete grids and take the form of time series. These often present a lot of heterogeneity, both within and across individuals.

(i) For every individual, the time between any two measurements can vary, and thus individuals may not be recorded on the same grid.

(ii) The number of total sampling points might vary between individuals.

(iii) Each measurement in time might be corrupted by measurement noise.

Mathematically, we consider n pairs of functions {(x 1 , y 1 ), . . . , (x n , y n )}. Each x i deterministically produces a specific y i through the Ordinary Differential Equation (ODE) (1). We call x i the feature path and y i the target path. Both x i and y i are only observed at a finite set of times specific to every individual. We denote by

D i = t i 1 , .
. . , t i ki , i = 1, . . . , n, the sampling grid of x i and by Di the sampling grid of y i . We stress that both the number of sampling times k i and the sampling times t i 1 , . . . , t i ki themselves are individual specific, as described in (i) and (ii). Moreover, the observations are corrupted by additive noise, such that we observe X i t = x i t + ξ i t for all t ∈ D i , and similarly Y i t = y i t + ε i t for every t ∈ Di , where the ξ i t and ε i t are sub-gaussian i.i.d. random vectors. Each input may therefore be written as a matrix X i = (X i t ) t∈D i ∈ R ki×d which we call the feature time series. Similarly, the quantity of interest is a matrix Y i = (Y i t ) t∈ Di ∈ R mi×p (where m i is the length of Di ) and is called the target time series. The grid Di is assumed to be a subset of D i : in our setup y i is hard to sample and therefore measured at only a few points (and sometimes only one) while x i is easy to access and measured at high frequency. Our goal is to approximate the dynamics linking x and y from the irregular, heterogeneous, and fuzzy data X i and Y i . Such heterogeneity is difficult to handle by classical machine learning algorithms such as Long short-term memory networks (LSTM, [START_REF] Hochreiter | Long short-term memory[END_REF] which assume that the data is regularly sampled. Some more recent approaches (Rubanova et al., 2019;De Brouwer et al., 2019;[START_REF] Kidger | Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu[END_REF][START_REF] Herrera | Neural jump ordinary differential equations: Consistent continuous-time prediction and filtering[END_REF] have adapted these models by introducing continuously evolving hidden states to account for the irregular spacing between observation times.

We build upon the approach of Neural Controlled Differential Equations (Neural CDE, [START_REF] Kidger | Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu[END_REF]Morrill et al., 2021b), which have proven to be very successful for time series classification and online prediction tasks (Morrill et al., 2021a). The key idea of Neural CDE is that under some fairly mild assumptions, any general ordinary differential equation of the form (1) can be rewritten as

y t = y 0 + t 0 F(y s )dx s , (2)
where F : R p → R p×d is a matrix-valued vector field, such that the right-hand-side of ( 2) is a matrix-vector product (see, e.g., [START_REF] Fermanian | Embedding and learning with signatures[END_REF], Proposition 2, for a proof). The function x is often called the driver of the CDE. In a Neural CDE setting, the driver x is a continuous interpolation of the feature times series, y corresponds to a continuously-evolving state, and F is chosen to be a neural network. This network is then trained such that the values of (y t ) can be used as features for classification or regression tasks. While Neural CDE have been shown to outperform other architectures with limited memory usage, their training time is considerable and no statistical guarantees exist.

Model. We model the interactions between the target and the feature paths through a CDE of the form (2). This modelling choice encapsulates a broad variety of settings, since the vector field F can be any (regular enough) function. A priori, the solution map Ψ of this CDE is a complex function of time and the history of x up to t; however, by linearizing the model, we are able to approximate Ψ by a simple scalar product between a deterministic transformation of the history of x, called the signature of x at order N ≥ 1 and denoted by S N (x [0,t] ), and a time independent matrix θ * N . Informally, we have

Ψ x [0,t] , t ≈ S N x [0,t] ⊤ θ * N .
Two striking features of this linearized model are (i) that θ * N can be learned on any time horizon [0, t] since it is independent of time, and (ii) that once it has been learned, the model can be called at any time t.

Contributions. Our contributions are threefold. First, we frame the task of learning the interactions between two time series as learning the flow of a CDE, which can be linearized in the signature space. While the connection between CDEs and signatures is well-known, this is the first time CDEs are used as a statistical model. We then leverage this linearization to derive statistical guarantees on the prediction error with an explicit dependence on both sampling irregularities and the noise affecting measurements.

To our knowledge, this is the first bound of this type for signature-based models, allowing for better understanding of the dependencies between prediction performance and sampling roughness. Finally, the resulting algorithm, called SigLasso, is shown to be computationally cheap and competitive compared to existing baselines on a wide range of simulated data and a real-world example of hospitalization growth rate prediction during the Covid pandemic.

Related works. Signatures originated as a prominent tool in stochastic analysis [START_REF] Chen | Integration of paths-a faithful representation of paths by non-commutative formal power series[END_REF][START_REF] Lyons | Differential Equations driven by Rough Paths[END_REF][START_REF] Friz | Multidimensional Stochastic Processes as Rough Paths: Theory and Applications[END_REF] and have proven to be a powerful feature extraction method in machine learning in various domains such as healthcare (Morrill et al., 2020b;[START_REF] Wang | Learning to detect bipolar disorder and borderline personality disorder with language and speech in non-clinical interviews[END_REF], human action recognition [START_REF] Yang | Developing the path signature methodology and its application to landmark-based human action recognition[END_REF], or financial modelling [START_REF] Lyons | A feature set for streams and an application to high-frequency financial tick data[END_REF][START_REF] Buehler | A data-driven market simulator for small data environments[END_REF]. Their appealing properties include a capacity to handle irregular data, to capture dependence between coordinates, and their links with the theory of CDE. We refer to [START_REF] Lyons | Signature methods in machine learning[END_REF] for a recent survey on their use cases. However, the statistical properties of signatures based algorithms have received little attention so far, with a few notable exceptions [START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF][START_REF] Lemercier | Distribution regression for sequential data[END_REF][START_REF] Fermanian | Functional linear regression with truncated signatures[END_REF].

On the other hand, the interplay between dynamical systems and machine learning has received considerable attention in the recent years. A first line of work has focused on approximating the solution of ODE and Partial Differential Equations (PDE) with neural networks [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF][START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF][START_REF] Zubov | NeuralPDE: Automating physicsinformed neural networks (PINNs) with error approximations[END_REF] and directly learning dynamical systems [START_REF] Long | PDE-net: Learning PDEs from data[END_REF][START_REF] Fattahi | Learning sparse dynamical systems from a single sample trajectory[END_REF]. Recent approaches have been interested in combining deep learning algorithms with physical knowledge (Greydanus et al., 2019;[START_REF] Brunton | Machine learning for fluid mechanics[END_REF][START_REF] Willard | Integrating physics-based modeling with machine learning: A survey[END_REF]. Finally, dynamical systems, seen as continuous versions of neural network architectures, have also been a great source of inspiration for analysing and designing machine learning algorithms in the recent years (Chen et al., 2018;[START_REF] Fermanian | Embedding and learning with signatures[END_REF][START_REF] Marion | Scaling resnets in the large-depth regime[END_REF]. We refer to [START_REF] Kidger | On neural differential equations[END_REF] for an extensive review.

We stress that our problem is different in nature from most problems encountered in the time series literature, since we seek to model the relationship between two sparsely observed systems with heterogeneous sampling. We do not model the evolution of the feature time series, and take it as an input, contrarily to methods such as Gaussian Process. Most models either focus on the case where one time series is observed and forecasted, or on regular sampling, or on univariate time series. Our problem bears close resemblance to frameworks encountered in sequence-tosequence learning (Sutskever et al., 2014;[START_REF] Gehring | Convolutional sequence to sequence learning[END_REF] and functional regression [START_REF] Ramsay | Some tools for functional data analysis[END_REF][START_REF] Marx | Generalized linear regression on sampled signals and curves: a P-spline approach[END_REF].

Overview. Section 2 introduces the CDE model for interacting systems, the mathematical context and the learning procedure. Our main theoretical result is presented in Section 3. We conclude by an empirical study on synthetic and real-world data in Section 4. All proofs are postponed to the appendix and the code to reproduce the experiments is available at https://github.com/LinusBleistein/SigLasso.

Model and Assumptions

A summary table of all notations introduced in the main body of this article is provided to the reader in Appendix A.

A CDE-Based Model on the Dynamics

We start by describing our assumptions on the feature and target paths, which are linked by Equation (2). To correctly define the integral of Equation ( 2), we need to impose some conditions on the x i and on F. Note that we consider that the x i are defined on [0, 1] but our results extend easily to any compact time interval [a, b].

Assumption 1. All paths (x i ) 1≤i≤n are continuous and there exists 0 < L < 1 such that, for all i = 1, . . . , n,

x i 1-var,[0,1] = sup D k x i t k+1 -x i t k ≤ L,
where ∥•∥ is the Euclidean norm and the supremum is taken over all finite dissections

D = {0 = t 1 < • • • < t k = 1}.
We write C 1-var L ([0, 1], R d ) for the set of continuous paths of total variation bounded by L. Outside the statistical context, when referring to general paths, we will drop the superscripts i and simply write x and y to alleviate notations.

We assume that the target path y is the solution of the ODE (1). This modelization choice means that the evolution of y is governed by a dynamical system whose dynamics itself are allowed to vary with the current value of the feature path. Observe that this model can be seen as a generalized form of a non-autonomous system [START_REF] Lyons | Differential Equations driven by Rough Paths[END_REF], which we recover by taking x t = t. Since Equation (1) can be rewritten as a CDE, the starting point of our work is to assume that the true dynamics of the data follow such a CDE, as stated in the following assumption.

Assumption 2. There exists a smooth vector field F : R p → R p×d such that, for all i = 1, . . . , n, y i is the solution of the CDE (2) driven by x i with initial condition y i 0 = y 0 ∈ R p homogeneous amongst individuals.

By "smooth" we mean that each coordinate of F is infinitely differentiable, that is, is C ∞ . The vector field F and the initial condition y 0 are common to all individuals, which can be seen as homogeneity assumptions on our sample. On the other hand, since every individual i has her own feature path x i , the target paths y i are individual specific. In other words, there exists a solution map Ψ that depends only on y 0 and F and is such that, for any t ∈ [0, 1], Ψ(x i [0,t] , t) = y i t . The vector field F encapsulates the common physical dynamics governing the evolution of y i , which are affected by the changes in x i . Note that there is no parametric model on F (although some strong smoothness requirements will be needed) contrarily to functional or traditional time series models [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Morris | Functional regression[END_REF].

Linearizing the CDE with Signatures

Before defining the Taylor expansion of the CDE (2), which will allow us to linearize the solution map Ψ, we need to introduce the notion of signature, which have emerged as a powerful tool to model time series [START_REF] Levin | Learning from the past, predicting the statistics for the future, learning an evolving system[END_REF]Kidger et al., 2019).

From now on, for any feature path x with values in R d , we denote by x (j) its jth coordinate, for j = 1, . . . , d. Definition 2.1.

Let x ∈ C 1-var L ([0, 1], R d ).
Take a word of length k from the alphabet {1, . . . , d}, that is, an element I = (i 1 , . . . , i k ) ∈ {1, . . . , d} k . For all t ∈ [0, 1], the signature coefficient associated to this word is the scalar

S I x [0,t] = • • • 0<u1<•••<u k <t dx (i1) u1 • • • dx (i k ) u k .
We introduce a series of notation on signature coefficients, grouping them by the length k of the words. For any k ∈ N and t ∈ [0, 1], the signature of order k of x [0,t] is

X k,[0,t] = S I (x [0,t] ) I∈{1,...,d} k ∈ R d k ,
where the words (i 1 , . . . , i k ) are in lexicographic order. We then denote the full signature by

S(x [0,t] ) = 1, X 1 [0,t] , . . . , X k [0,t] , . . .
and, for any N ≥ 1, the signature truncated at order N by

S N (x [0,t] ) = 1, X 1 [0,t] , . . . , X N [0,t] ⊤ .
The size of the signature truncated at order N grows exponentially with N , since, for d ≥ 2, it is equal to

s d (N ) = 1 + d + d 2 + • • • + d N = d N +1 -1 d -1 .
When computing signatures, it is common practice to add time as a coordinate to the path [START_REF] Fermanian | Embedding and learning with signatures[END_REF], that is, consider the path (t, x t ) ⊤ . From now on, we assume that the first dimension of x always corresponds to time, so that d ≥ 2.

Signatures encode geometric properties of paths and have numerous appealing properties as a feature set for time series. We refer to [START_REF] Chevyrev | A primer on the signature method in machine learning[END_REF]; Fermanian (2021); [START_REF] Lyons | Signature methods in machine learning[END_REF] for more detailed introductions to signatures. In order to provide supplementary intuition, we first give three insights on signatures.

A geometric insight. Signatures are a geometric alternative to representations based on frequency such as the Fourier transform. Indeed, consider a function f ∈ C ∞ ([0, 1], [0, 1]) and the two dimensional path (t, f (t)).

For simplicity assume that f (0) = 0. Then, the first order signature coefficients are equal to the last positions t and f (t). The second order signature coefficients are equal to

t 0 f (s)ds, 1 2 t 2 , 1 2 f (t) 2
, and f (t)t -t 0 f (s)ds, and capture how the area under the curve evolves with time.

A computational insight. Consider the linear path

x t = at + b for t ∈ [0, 1], where a = (a 1 , . . . , a d ) ⊤ ∈ R d . For any word (i 1 , . . . , i k ) of size k the associated signature coefficient is a i1 . . . a i k k! t k .
In the linear case, signatures are therefore simply polynomials in t with path-specific coefficients. This result generalizes nicely to piecewise linear paths (via a result known as Chen's Lemma) and allows for computational efficiency when computing the signature. We refer to Kidger & Lyons (2020) for further computational details.

A functional insight. Recall that we are interested in approximating functions f (x t , t) which depend both on time and on the values of the feature path (x t ). When computing signatures, we always consider the time-augmented path (t, x t ). The coefficients will thus be divided in two parts: a first set of coefficients related to the time dimension, and a set of coefficients related to the path dimensions. The time-specific coefficients are simply

t, t 2 2! , t 3 3! , . . . , t N N ! (3)
and thus form a polynomial basis. Roughly speaking, these coefficients approximate the time dependant part of the function f . The path-specific coefficients, which can be thought of as polynomials of a path, approximate the part of f depending on (x t ). This highlights that the Taylor development of a CDE, which is the cornerstone of our work, is very similar in nature to the approximation of a function by its classical Taylor development (see Appendix B.4).

With these insights in hand, we are now ready to properly define the Taylor expansion of a CDE.

Definition 2.2. Let N ≥ 1. The Taylor expansion of order N of the solution y of Equation ( 2) is defined by

y N,t = y 0 + N k=1 I∈{1,...,d} k S I x [0,t] × Φ I F (y 0 ) , (4)
where Φ I F (•) ∈ R p is the differential product of the vector field F along I, whose definition is postponed to Appendix B.3.

The differential products Φ I F (y 0 ) are essentially a combination of multiplication and summation of derivatives of the different components of F, evaluated at y 0 .

The Taylor expansion crucially allows to write the solution map as a product between a time-varying term (the signature of the feature path) and a constant-over-time term (the differential product). This is similar to a regular Taylor expansion; we discuss this analogy in Appendix B.4. Note that Equation (4) may also be written as a matrix-vector product

y N,t ⊤ = S N (x [0,t] ) ⊤ θ * N ∈ R p , (5)
where θ * N ∈ R s d (N )×p is a matrix collecting all differential products Φ I F (y 0 ) up to order N and the offset y 0 . Since θ * N depends neither on x nor on t, the Taylor expansion of y at order N is simply a linear function of the truncated signature.

For this expansion to become exact in the sens of pointwise convergence, we need quite strong regularity assumption on F. This is the price to pay for the non-parametric nature of our model. Assumption 3. The vector field F has fast decaying derivatives. In other words, defining

Λ k (F) = sup I∈{1,...,d} k Φ I F (y 0 ) , we assume that ∞ k=0 d k Λ k (F)/k! < ∞.
Let y be the solution of the CDE (2), and let (y N,t ) t∈[0,1] be its N -th Taylor expansion. Under Assumptions 1, 2, and 3, we then have, for any t ∈ [0, 1],

y t -y N,t -→ N →+∞ 0. (6)
We refer the reader to [START_REF] Friz | Euler estimates for rough differential equations[END_REF] and Fermanian et al. (2021, Proposition 4) for more details on this result.

While we impose relatively strong regularity assumptions on the vector field F, the assumptions on x are mild (see Assumption 1). Therefore, our assumptions, while enforcing a fairly high amount of regularity on the dynamical structure, still accommodate most of the real world data. Also note that [START_REF] Fermanian | Embedding and learning with signatures[END_REF] give conditions under which Assumption 3 is satisfied when F is a layer of a neural network with smooth activation functions.

The Learning Problem

We go back to the statistical learning problem. Recall that in practice, we do not observe the continuous paths {(x 1 , y 1 ), . . . , (x n , y n )} but their discretized and fuzzy counterparts {(X 1 , Y 1 ), . . . , (X n , Y n )} measured on a set of individual grids D = {D 1 , . . . , D n } and D = { D1 , . . . , Dn }. We allow for high data heterogeneity: two individuals can be sampled at very different frequencies and at different time-points, and therefore have different observation grids.

The meshsize of a sampling grid D, denoted by |D| is defined as the largest gap between two successive sampling times, that is,

|D| = max ti∈D |t i+1 -t i |.
Its cardinality, denoted by #D, is the number of sampling points in D. We make the following assumption on the sampling procedure.

Assumption 4. There exists η ∈ [0, 1] such that for all i = 1, . . . , n, one has

0 ∈ D i , #D i ≥ 2, t i ki ≥ η and Di ⊂ D i .
Let us briefly comment on this assumption. We require that the measurements on all individuals start at 0. We also do not allow for arbitrarily short time series: for every individual i, the last observation time t i ki must be greater than η. Finally, we require the sampling grid of Y i to be coarser than the sampling grid of X i . We also let |D| = max i |D i | be the biggest gap between two successive sampling times within the whole set of individual grids and #D = n i=1 #D i the total number of sampling points of the feature time series. We recall that the random vectors ξ i t , t ∈ D i , are the noises affecting the measurements of x i . The random vectors ε i t affect the measurements of y i . We end with assumptions on the law of these measurement noises.

Assumption 5 (Noise on the feature time series). The noises (ξ i t ) i∈{1,...,n},t∈D i are i.i.d. v ξ -subgaussian random vectors.

Assumption 6 (Noise on the target time series). The noises (ε i t ) i∈{1,...,n},t∈ Di are i.i.d. v ε -subgaussian random vectors and independent from (ξ i t ) i∈{1,...,n},t∈D i .

The goal of the learning procedure is to learn the solution map Ψ, in order to infer the value of y at any time t, given observations of x up to time t. We have seen previously that this problem boils down to estimating the matrix θ * N via Equation (5). The truncation order N ≥ 1 is an hyperparameter and will be selected by cross-validation. Single target measurement. For the sake of simplicity, we first present the case where Y i is measured only at the end of the observation period, as in the example of LA measurements in obstetrics. In this case, for all i = 1, . . . , n, we have

m i = 1 such that M = n, D i = (t i 1 , . . . , t i ki )
, and Di = (t i ki ).

Then Y i ∈ R p , and we denote by

Y = Y 1 , . . . , Y n ⊤ ∈ R n×p
the matrix containing all target measurements. Since we do not have access to the feature paths x i but only to the discrete measurements X i , we compute the signature of its linear interpolation normalized by its total variation, sampled up to final time t i ki , and denote by S D N ∈ R n×s d (N ) the matrix of stacked signatures. Note that signatures of piecewise linear functions are fast to compute with packages such as signatory (Kidger & Lyons, 2020) or iisignature [START_REF] Reizenstein | Algorithm 1004: The iisignature library: Efficient calculation of iteratedintegral signatures and log signatures[END_REF]. The complexity to compute the signature truncated at order N of the ith feature time series is of the order O(#D i d N ). Finally, we approximate θ * N ∈ R s d (N )×p by solving the optimisation problem

min θ∈R s d (N )×p 1 2n Y -S D N θ 2 F + Ω(θ), (7)
where Ω : R s d (N )×p → R + is a regularization term and ∥•∥ F is the Frobenius norm. We have reduced the complex problem of learning the solution map Ψ to a simple penalized linear regression in the signature space. This linear model on the signature is close to the one studied by [START_REF] Levin | Learning from the past, predicting the statistics for the future, learning an evolving system[END_REF]; [START_REF] Lyons | A feature set for streams and an application to high-frequency financial tick data[END_REF][START_REF] Fermanian | Functional linear regression with truncated signatures[END_REF].

Multiple target measurements. We also cover the case when # Di > 1, that is, the target is measured at multiple times for every individual, as in the example of hemorrhage detection. In this case, we have Y i ∈ R mi×p and we stack the different measurement matrices

Y i to obtain a matrix Y of size M × p, where M = m 1 + • • • + m n .
For any i = 1, . . . , n and every t ∈ Di , we predict Y i t using the signature of the linear interpolation of the normalized (X i 0 , . . . , X i t ). In this manner, we will be able to predict Y i t at every point where X i t is sampled. The exact workflow of our model is described in Figure 1.

Theoretical Guarantees

Mathematical Setup

We consider a general multiple target measurements setting. To simplify the exposure of our results, we consider a univariate target path, i.e., p = 1. In this case, the true parameter θ * N is a vector of size s d (N ) and not a matrix. The general case p ≥ 1, which our algorithm handles as running p Lasso regressions in parallel, is considered in Appendix C, and all theoretical results are proved in this general case. In addition, to lighten the presentation of the oracle inequality, we also focus in this section on the case of ω-Lipschitz feature paths, that is, for every i = 1, . . . , n and for all s, t ∈ [0, 1], x i tx i s ≤ ω|t -s|. We stress that our results are valid for continuous paths of bounded variations. We let y ∈ R M be the matrix collecting all unobserved values of the target paths at measurement times such that y = E Y , where the expectation is taken over the noises ε i t , and define θ N,M as θ N,M ∈ arg min

θ∈R s d (N ) 1 2M Y -S D N θ 2 2 + Ω(θ). (8)
For δ ∈ (0, 1), we define the set

A ξ (δ) = max ξ i t ≤ v ξ √ d + v ξ 1 c log #D δ =:C δ (9)
where the maximum is taken on all i = 1, . . . , n and t ∈ D i , and c is a universal constant. This set is of probability greater than 1δ under Assumption 5 (see Appendix B.5).

Similarly, for k ≥ 0 and δ ∈ (0, 1), let

C k ( δ) = v ε log(2N d k / δ)
and define

A ε ( δ) = N k=0 ε ⊤ S D •,[k] ∞ ≤ M 1 2 C k ( δ) k! , (10)
where

S D •,[k] is the sub-matrix of size M × d k of S D
N associated to the signature coefficients of order k, and ε ∈ R M is a vector of i.i.d. noise terms satisfying Assumption 6 (see Appendix C.1). Under Assumptions 5 and 6, A ε ( δ) is of probability at least 1 -δ, and

A ξ (δ) ∩ A ε ( δ) is of proba- bility at least (1 -δ)(1 -δ). Let Ω(θ) = N k=0 C k ( δ) k! √ M θ [k] 1 , (11)
where θ [k] is the subvector of size d k that collects all elements of θ associated to words of size k (see Appendix C.1). This penalization can be implemented by rescaling the feature matrix S D N and solving a standard ℓ 1 -penalized regression problem (see Appendix D.1). Our result extends to more general penalties by adapting existing techniques from [START_REF] Chesneau | Some theoretical results on the grouped variables lasso[END_REF] for the group-lasso, or [START_REF] Lederer | Oracle inequalities for high-dimensional prediction[END_REF] for the hierarchical lasso.

Main Results

The error made when learning θ * N by θ N,M comes from three different sources. (i) Truncating the signature used in the regression at depth N ≥ 1 results in a truncation bias. (ii) Discretization of the feature path and the noise affecting each measurement point induce a discretization error. In particular, there is a trade-off between sampling frequency and variance of the noise. (iii) The measurement error on y i and the finite-sample setting induce a classical estimation error.

The following lemmas bound each of those errors. We first bound the variance of the estimator with arguments borrowed from [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF].

Lemma 3.1. Under Assumptions 1 and 2, on the set A ε ( δ) ∩ A ξ (δ), the prediction error

1 2M y -S D N θ N,M 2 2
is bounded above by

1 2M y -S D N θ * N 2 2 + 2C N ( δ) √ M N k=0 d k Λ k (F) k! .
See Appendix C.2 for a proof. This inequality decomposes the error into a bias and a variance term. We denote the signature matrix of the unobserved paths N ) .

x i by S N ∈ R M ×s d (
As for the bias term, notice that one can write

1 2M y -S D N θ * N 2 2 ≤ 1 M ∥y -S N θ * N ∥ 2 2 Truncation bias + 1 M S N θ * N -S D N θ * N 2 2
Discretization error

.

Bounding each of these terms corresponds to, respectively, Lemmas 3.2 and 3.3. We stress that the truncation bias is of a different nature than the discretization error since it depends on a choice of hyperparameter while the latter is inherent to the data at hand.

Lemma 3.2. Under Assumptions 1 and 2, for any

N ≥ 1, 1 M ∥y -S N θ * N ∥ 2 2 ≤ d N +1 Λ N +1 (F) (N + 1)! 2 .
Under Assumption 3, the right-hand-side decays exponentially fast with N . This lemma is an immediate consequence of [START_REF] Fermanian | Embedding and learning with signatures[END_REF] (see Appendix C.3). We now turn to the error induced by the discretization of the feature path.

Lemma 3.3. Under Assumptions 1, 4, and 5, on the set

A ξ (δ), one has 1 M (S N -S D N )θ * N 2 2 ≤ C D,N (δ) N k=0 d k Λ k (F) 2 k! 2 , where C D,N (δ) is equal to 4e 2 L 2 N ! 2 × ω|D| + C δ + 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ 2 .
This lemma relies on a fine analysis of the distance between two signature layers. The dependence of C D,N (δ) on sampling mechanisms and noise is of particular interest. First, The second part of C D,N (δ) is a consequence of the noises ξ i t affecting the measurement points of the feature time series and does not vanish with n. The emergence of such a bias is a well studied phenomenon in errors-in-variable models, and cannot be corrected without precise knowledge of the noise's variance (Loh & Wainwright, 2011).

The last term corresponds to a bias coming from the interplay of the normalisation of the signatures by the total variation of the path and the noise, which is a standard practice (Morrill et al., 2020a). We found that this normalization performs best empirically. Note that this bias term is equal to 0 if the path has total variation exactly equal to 1 and is observed without noise.

From Lemmas 3.1 to 3.3 and the definitions of A ξ and A ε finally we get the following oracle inequality. The proof is given in Appendix C.6. Theorem 3.4 (An oracle inequality for learning with signatures). Under Assumptions 1, 2, 4, 5, and 6, with probability at least (1δ)(1 -δ), the prediction error is bounded above by

1 2M y -S D N θ N,M 2 2 ≤ d N +1 Λ N +1 (F) (N + 1)! 2 (12) + C D,N (δ) N k=0 d k Λ k (F) 2 k! 2 (13) + 2C N ( δ) √ M N k=0 d k Λ k (F) k! . ( 14
)
All terms depend on the regularity of the vector field F via the constants Λ k (F): the bigger these constants, the faster the vector field F may vary, making the CDE harder to predict. The convergence speed in 1/ √ M is classical. Also note that our bound is non-asymptotic and is valid for any M ≥ 1.

The dependence of the bound on N is highly non-trivial and requires an in-depth analysis of the regularity of F in order to bound Λ k (F), which is out of the scope of this paper. The asymptotic behaviour of this oracle inequality is discussed in Appendix C.7.

Experiments

We study the performance of SigLasso obtained by solving the optimization problem (8), where Ω(θ) is defined by Equation ( 11). All details are given in Appendix D.

Simulations

We consider several settings of data generation. First, in the well-specified setting, the data is generated from a model with regular feature paths x (piecewise polynomials) and target paths y wich are solutions to the CDE

dy t = tanh(Ay t )dx t ,
where A is a randomly drawn matrix. In the ill-specified setting, the target y is equal to

y t = log ∥ 10 h=1 x t-h ∥ for any t ∈ [0, 1].
In the third setting, called OU setting, the feature paths are realizations of Brownian motions and the target paths are Ornstein-Uhlenbeck processes [START_REF] Borodin | Handbook of Brownian Motion-Facts and Formulae[END_REF] driven by the feature paths. The last setting corresponds to the tumor growth model from [START_REF] Simeoni | Predictive pharmacokineticpharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents[END_REF]. The feature path represents the concentration of a treatment drug, generated as the squared value of the smooth paths used in the well-specified setting, and the target path y, the weight of the tumor, is governed by a system of differential equations given in Appendix D.6.

We compare SigLasso to a GRU and Neural CDE [START_REF] Kidger | Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu[END_REF]. We measure the performance of the models with two metrics on a test set: the mean squared error for predicting the last observation point of the target paths and the L 2 error for predicting the full path on a fine grid.

The results are shown in Table 1 andFigure 2. In Figure 2 we consider the well-specified setting and vary the number of sampling points of the target paths between 1 and 20. In Table 1 it is fixed to 5 but the simulation settings change. Sampling of both the target and the feature time series is highly irregular. SigLasso outperforms Neural CDE and GRU models in generalizing from a few learning points of the target to its full trajectory in all settings. We conduct supplementary experiments with RNN and LSTM, which are also outperformed by SigLasso (see Table 2 in Appendix D.9). An additional byproduct of SigLasso's simple form is its training speed: its is approximately 10 times faster than GRU and 100 times faster than Neural CDE, including cross-validation to select N and regularization strength (see Appendix D.9).

Forecasting the Growth Rate of Hospitalizations in

France During the Covid-19 Pandemic

Forecasting hospitalizations in real time during the Covid-19 pandemic is a notably difficult task. In this experiment, we train our model to learn the dynamics linking population data related to mobility, vaccination, and weather, and the hospitalization growth rate (HGR) in each of the 9 metropolitan regions of France based on the data of [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF]. The feature time series is regularly sampled and specific to each region and 12-dimensional. We consider prediction horizons h = 1, . . . , 14, meaning that we predict the hospital saturation at time t using the history of the feature time series up to th. Using our notations, we have for each region d = 12, p = 1, n = 1. Concretely, this means that we train one Siglasso model per region and per horizon. The target is sampled every day during the training period left of the dotted line in Figure 3, and the models are fitted to those values: on this time span, the models learns to interpolate the target time series. The model is then asked to predict the HGR on the days right of the dotted line seing only the feature time series. This means that it performs a prediction task on this time span.

Both GRU and SigLasso learn smooth and precise dynamics, which generalize well above the learning horizon and yield similar prediction performance (see Appendix D.8). Figure 3 shows an example of reconstruction and prediction of HGR obtained with SigLasso, GRU and NCDE.

Conclusion

We have introduced a novel CDE-based model for interacting systems. Drawing on the theory of signatures, we derive an oracle bound that depends explicitly on the roughness of the data sampling. We illustrate the high performance of our approach on synthetic and real-world data.

The obtained theoretical guarantees rely on strong regularity assumptions on the vector field F. The exact approximation properties of this class of vector field are a very interesting direction for future work. Considering other penalties that take into account the underlying structure of θ * N would also be an interesting extension of our work.

Supplementary Material

A. Summary of used notations

The following table provides an exhaustive list of notations used in the main body of the paper. Notations are grouped by subsections. Matrix of measurements of the target time series for individual i -

∥x∥ 1-var,[0,t]
Total variation of the path x on the interval [0, t] Section 2, page 3, Assumption 1

C 1-var([0,1],R d ) L
Set of continuous paths of total variation bounded by L -F Unknown smooth generative vector field Section 2, page 3, Assumption 2

S I (x [0,t] ) Signature coefficient of the path x on [0, t] associated to the word I Section 2, page 4, Definition 2.1 X k,[0,t] Signature of order k - S(x [0,t] )
Full signature at order N -SN Siglasso estimator Section 3, page 6, Equation (7)

A ξ (δ)
Set bounding the magnitude of noises (ξ i t ) Section 3, page 6, Equation ( 9)

C k ( δ) Constant Section 3, page 6 Aε(δ)
Set bounding the magnitude of the noises (ε i t ) Section 3, page 7, Equation ( 10) Ω(θ)

Penalty evalutated at θ Section 3, page 7, Equation (11)

B. Mathematical details B.1. The Riemann-Stieltjes integral

We fitst recall two key properties on the Riemann-Stieljes integral. For a general presentation of the Riemann-Stieltjes integral, we refer to [START_REF] Friz | Multidimensional Stochastic Processes as Rough Paths: Theory and Applications[END_REF].

Proposition B.1. Let x ∈ C 1-var L ([0, 1], R d ) and y : [0, 1] → R d be a continuous path. Then t 0 y s dx s ≤ ∥y∥ ∞,[0,t] ∥x∥ 1-var,[0,t]
We refer the reader to Friz & Victoir (2010, Proposition 2.2) for a proof.

B.2. The truncated tensor algebra

This section introduces notations and definitions on the space in which signatures are defined, namely, the tensor algebra.

While for the exposition of our main results, the truncated signature of a path N ) , it is often useful to place ourselves in the tensor algebra to obtain finer bounds or technical results.

x ∈ C 1-var L ([0, 1], R d ) at depth N ≥ 1 can be assimilated to an element of R s d (
Let x ∈ C 1-var L ([0, 1], R d
) be a path of bounded variation. For a word I = (i 1 , . . . , i k ) ∈ {1, . . . , d} k of size k, the signature coefficient S I (x [0,1] ) can be seen as an element of the k-th tensor product of R d with itself, denoted by R d ⊗k .

For instance, the coefficients of order k = 1 can be written as a vector and the coefficients of order k = 2 as a matrix, and so on, i.e.,

X 1 [0,1] =     1 0 dx (1) s
. . .

1 0 dx (d) s     and X 2 [0,1] =     1 0 dx (1) s dx (1) s
. . .

1 0 dx (1) s dx (d) s . . . . . . 1 0 dx (d) s dx (1) s . . . 1 0 dx (d) s dx (d) s     .
We now define a norm on (R d ) ⊗k . Let a ∈ (R d ) ⊗k and (e 1 , . . . , e d ) be the canonical basis of R d . Then (e i1 ⊗ • • • ⊗ e i k ) (i1,...,i k )∈{1,...,d} k is a basis of (R d ) ⊗k . We can thus write a as a = (a I ) I∈{1,...,d} k . For every k ≥ 0, the vector space (R d ) ⊗k is naturally endowed with the norm

∥a∥ 2 (R d ) ⊗k = I∈{1,...,d} k a I 2 .
Remark that this norm satisfies for any x ∈ (R d ) ⊗k and y

∈ (R d ) ⊗m , ∥x ⊗ y∥ (R d ) ⊗(k+m) = ∥x∥ (R d ) ⊗k ∥y∥ (R d ) ⊗m . ( 15
)
We refer to [START_REF] Fermanian | Embedding and learning with signatures[END_REF] for further deails. The signature truncated at depth N ≥ 1 collects elements from R, R d ⊗2 , . . . , R d ⊗N . It can thus be seen as an element of the truncated tensor algebra

T N (R d ) = R ⊕ R d ⊗2 ⊕ • • • ⊕ R d ⊗N .
Let a = (a 0 , . . . , a N ) ∈ T N (R d ), where every a k ∈ (R d ) ⊗k . We define the norm

∥a∥ T N (R d ) = N k=0 ∥a k ∥ 2 (R d ) ⊗k 1/2 .
To clarify, if we consider the truncated signature of x at depth N ≥ 1, which is an element of T N (R d ), then

S N (x [0,t] ) T N (R d ) = N k=0 X k [0,t] 2 (R d ) ⊗k 1/2 = N k=0 I∈{1,...,d} k S I (x [0,t] ) 2 1/2 .
Note that this norm is exactly equivalent to the Euclidian norm of R s d (N ) , which is the space we consider in the exposition of our main results for the sake of simplicity.

We are now ready to define the tensor product on the truncated tensor algebra. For two elements a = (a 0 , . . . For any k = 0, . . . , N , we let

π k : T N (R d ) → (R d ) ⊗k be the canonical projection of T N (R d ) onto (R d ) ⊗k . More precisely, for every a = (a 0 , . . . , a N ) ∈ T N (R d ), π k (a) = a k
We also define the canonical projection

Π k : T N (R d ) → T k (R d ) defined by Π k (a) = (a 0 , . . . , a k ) .

B.3. The differential product

We first define the differential product, in order to give a precise statement of Assumption 3. Definition B.3. Let F, G : R p → R p be two smooth vector fields, i.e., each of their components is C ∞ . Denote by J(•) the Jacobian matrix. The differential product F ⋆ G : R p → R p is the smooth vector field defined for any h ∈ R p

(F ⋆ G)(h) = e j=1 ∂G ∂h j (h)F j (h) = J(G)(h)F (h).
The differential product is not associative. We therefore use the convention to evaluate it from right to left, that is,

F 1 ⋆ F 2 ⋆ F 3 = F 1 ⋆ F 2 ⋆ F 3 .
Let F : R p → R p×d be a smooth vector field. We write F 1 , . . . , F d the columns of F. Every F i , for i = 1, . . . , d, can thus be seen as a map from R p to R p . Recall that y 0 ∈ R p is the initial condition of the CDE defined in Assumption 2. Let I = (i 1 , . . . , i k ) ∈ {1, . . . , d} k . We now define

Φ I F (y 0 ) = F i1 ⋆ • • • ⋆ F i k (y 0 ) ∈ R p .
We refer to [START_REF] Fermanian | Embedding and learning with signatures[END_REF] for greater details on the differential product. We now define for all k ≥ 1

Λ k (F) = sup 1≤i1,...,i k ≤d Φ I F (y 0 ) ∈ R, (16)
and use the convention Λ 0 (F) = ∥y 0 ∥. Remark that Assumption 3 implies that

d N +1 Λ N +1 (F) (N + 1)! -→ N →+∞ 0.
As an immediate consequence, the truncation bias

d N +1 Λ N +1 (F) (N + 1)!

B.4. Analogy with the Taylor extension

The definition of the Taylor expansion of a CDE exposed in the previous subsection is technical. However, it can simply be thought of as a generalization of the the classical Taylor expansion of a C ∞ function f : R → R. Recall that in this case, the Taylor expansion at 0 evaluated at t ∈ R of order N ∈ N writes as a power series

f (t) ≈ f (0) + f ′ (0) 1! t + • • • + f (N ) (0) N ! t N . ( 17
)
Every element of this power series is a product of two terms: the derivatives of f encode some information about the regularity of f at the initial point 0 and do not depend on t, while the polynomial terms t k allow this linearized form to evolve with time t. Remark that these polynomial terms do not depend on f .

Similarly, Equation ( 4) is also a sum of products of two terms. On the one hand, the evolving nature of the system, instead of being handled by the polynomial terms t k , are now captured by the signature coefficients S I (x [0,t] ). As the polynomial terms, they do not depend on F. On the other hand, the information about the initial value of the system at time t = 0 and the dynamics of F are summarized by the differential product Φ I F (y 0 ), which play the same role as the successive derivatives in Equation ( 17). To capture the multivariate nature of the paths, the Taylor expansion is summed over multi-indexes, or words, I = (i 1 , . . . , i k ) ∈ {1, . . . , d} k of size k for k ∈ N.

B.5. Properties of subgaussian random vectors

We start with the definition of a subgaussian random variable, see [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF] for more details.

Definition B.4. A real-valued random variable X is said to be σ 2 -subgaussian if for all t > 0

P(X > t) ≤ exp(-t 2 /σ 2 ), or, equivalently, if for all t ∈ R E(e tX ) ≤ exp(-ct 2 σ 2 ),
where c is an universal constant. A random vector Z is subgaussian if, for any vector c of norm 1, ⟨Z, c⟩ is subgaussian.

The norm of a sequence of d subgaussian random variables concentrates around √ d, as stated by the following lemma.

Lemma B.5. Let X 1 , . . . , X n be a sequence of i.i.d. σ 2 -subgaussian random variables. Let X = (X 1 , . . . , X d ) ∈ R d . There exists a universal constant c such that for all t > 0

P(∥X∥ 2 ≥ t + σ √ d) ≤ exp(-ct 2 /σ 2 ).
Proof. We refer to Vershynin (2018, Theorem 3.1.1) for a proof.

We can use this lemma to bound the maximum of n sequences of d subgaussian random variables with high probability.

Lemma B.6. Let X 1 , . . . , X n be a sequence of i.i.d. σ 2 -subgaussian random variables, such that for all i = 1, . . . , n, X i = (X i1 , . . . , X id ). Then there exists a universal constant c such that for all δ ∈ (0, 1)

P max i=1,...,n ∥X i ∥ ≤ σ √ d + σ 1 c log(n/δ) ≥ 1 -δ.
Proof. Using Lemma B.5 and a union bound, we have

P max i=1,...,n ∥X i ∥ ≥ σ √ d + σ 1 C log(n/δ) = P n i=1 ∥X i ∥ 2 ≥ σ √ d + σ 1 C log(n/δ) ≤ n i=1 P ∥X i ∥ 2 ≥ σ √ d + σ 1 c log(n/δ) ≤ δ,
which yields the desired inequality.

Notice that the universal constant is identical between both lemmas. As a consequence of this last lemma, under Assumption 5, the set

A ξ (δ) = max i=1,...,n,t∈D i ξ i t ≤ v ξ √ d + v ξ 1 c log(#D/δ) (18)
where

#D = n i=1
#D i is of probability at least 1δ.

We also need the following lemma.

Lemma B.7. Let X 1 , . . . , X n be a sequence of i.i.d. σ 2 -subgaussian random variables. Let Z 1 , . . . , Z n be random variables such that for all i = 1, . . . , n,

|Z i | ≤ α almost surely. Then n i=1 X i Z i is nσ 2 α 2 -subgaussian.
Proof. We use the characterization of subgaussian random variables by their characteristic function. For all t > 0,

E e t n i=1 XiZi = E n i=1 E e tXiZi | Z 1 , . . . , Z n ≤ E n i=1 E e tXiα ≤ E e ct 2 nα 2 σ 2 .
This finally yields that

E e t XiZi ≤ E e ct 2 nα 2 σ 2 ,
which concludes the proof.

C. Proofs

C.1. Preliminary notations

Let (E, ∥•∥ E ) be a normed vector space and x : [0, 1] → E. The supremum norm of x is defined for all t ∈ [0, 1] as

∥x∥ ∞,[0,t] = sup s∈[0,t] ∥x s ∥ E .
When referring to the total variation ∥x∥ 1-var,[0,1] of a path x : [0, 1] → R d over the whole domain, depending on the mathematical context, we will sometimes drop the time subscript and simply write ∥x∥ 1-var .

When referring to a matrix A = (A ij ) ∈ R n×p , we define classicaly the infinite and Frobenius norms by We now introduce some notations to take advantage of the structure of θ * N . The true parameter of the Taylor expansion of the model CDE, defined in Equation ( 5), can be written in block notation as

∥A∥ ∞ = max
θ * N =                                       θ * [0],1 • • • θ * [0],p θ * [1],1 • • • θ * [1],p θ * [2],1 • • • θ * [2],p . . . θ * [N ],1 • • • θ * [N ],p                                       ∈ R s d (N )×p , where θ * [k],ℓ ∈ R d k ×1 , k = 0, . . . , N, ℓ = 1, . . . , p. (19)
Every column of θ * N corresponds to a dimension of the target, while blocks of lines correspond to signatures layers. Thus for every k = 0, . . . , N and ℓ = 1, . . . , p, θ *

[k],ℓ is a column vector of size d k .

Similarly, for a general θ ∈ R s d (N )×p and the SigLasso estimator θ N,M , we will refere to the blocks forming these matrices as respectively θ [k],ℓ and θ [k],ℓ , for k = 0, . . . , N and ℓ = 1, . . . p.

Likewise, the signature feature matrix S D N ∈ R M ×s d (N ) can be written in block notation as

S D N = 1 S D •,[1] S D •,[2] • • • S D •,[N ] =    1 S D 1,[1] S D 1,[2] S D 1,[N ] . . . . . . . . . • • • . . . 1 S D n,[1] S D n,[2] S D n,[N ]    ,
where for any k = 1, . . . , N ,

S D •,[k] ∈ R M ×d k and, for every individual i = 1, . . . , n, S D i,[k] ∈ R mi×d k
(recall that m i is the number of measurements of the target path y i ). More precisely, given her target sampling grid Di = ( ti 1 , . . . , ti mi ), the individual-specific signature block of depth k is equal to

S D i,[k] =     1 S (1) (X i [0, ti 1 ] ) • • • S (d,...,d) (X i [0, ti 1 ] ) . . . . . . . . . 1 S (1) (X i [0, ti m i ] ) • • • S (d,...,d) (X i [0, ti m i ] )     ,
where the path t → X i t is a linear interpolation of the observed time series X i . The same notations will be used for the true signature feature matrix S N . We use the bracket notation [•] both in θ * N and S D N to emphasise that both the columns of the feature matrix and the lines of learned parameter correspond to words of the alphabet {1, . . . , d}.

The unobserved matrix of true values of the target writes as

y =    y 1 . . . y n    =    y 1 1 • • • y 1 p . . . . . . y n 1 • • • y n p    =                y 1 1, t1 1 • • • y 1 p, t1 1 . . . . . . y 1 1, t1 m 1 • • • y 1 p, t1 m 1 . . . . . . y n 1, tn 1 • • • y n p, tn 1 . . . . . . y n 1, tn mn • • • y n p, tn mn                ∈ R M ×p (20)
and the measurement matrix Y ∈ R M ×p can be written in a similar fashion as

Y =    Y 1 . . . Y n    =    Y 1 1 • • • Y 1 p . . . . . . Y n 1 • • • Y n p    =                Y 1 1, t1 1 • • • Y 1 p, t1 1 . . . . . . Y 1 1, t1 m 1 • • • Y 1 p, t1 m 1 . . . . . . Y n 1, tn 1 • • • Y n p, tn 1 . . . . . . Y n 1, tn mn • • • Y n p, tn mn                =                y 1 1, t1 1 + ε 1 1, t1 1 • • • y 1 p, t1 1 + ε 1 p, t1 1 . . . . . . y 1 1, t1 m 1 + ε 1 1, t1 m 1 • • • y 1 p, t1 m 1 + ε 1 p, t1 m 1 . . . . . . y n 1, tn 1 + ε n 1, tn 1 • • • y n p, tn 1 + ε n p, tn 1 . . . . . . y n 1, tn mn + ε n 1, tn mn • • • y n p, tn mn + ε n p, tn mn                (21) C.2. Proof of Lemma 3.1
Using the definition of Λ k (F) (see Equation ( 16)), we get the following proposition which allows to obtain an explicit dependence of the oracle bound on the regularity of F.

Proposition C.1. Let θ * N be defined as in Equation (5). Then

∥θ * N ∥ 2 F ≤ N k=0 d k Λ k (F) 2 ,
and, for all k = 0, . . . , N and ℓ = 1, . . . , p, θ *

[k],ℓ 1 ≤ d k Λ k (F).
Proof. By definition,

∥θ * N ∥ 2 F = N k=0 1≤i1,...,i k ≤d F i1 ⋆ • • • ⋆ F i k (y 0 ) 2 2 .
Since for all (i 1 , . . . , i k ) ∈ {1, . . . , d} k ,

F i1 ⋆ • • • ⋆ F i k (y 0 ) 2 2 ≤ Λ k (F) 2 , we get N k=0 1≤i1,...,i k ≤d F i1 ⋆ • • • ⋆ F i k (y 0 ) 2 2 ≤ N k=0 1≤i1,...,i k ≤d Λ k (F) 2 ≤ N k=0 d k Λ k (F) 2 .
We now turn to the second inequality. For k = 0, the inequality holds by definition. For k = 1, . . . , N and ℓ = 1, . . . , p, by definition of the ℓ 1 norm,

θ * [k],• 1 = 1≤i1,...,i k ≤d Φ I F (y 0 ) 1 , This yields θ * [k],• 1 ≤ d k Λ k (F)
and thus θ *

[k],ℓ 1

≤ d k Λ k (F) for ℓ = 1, . . . , p.
The following lemma is needed to leverage classical proof techniques to bound the prediction error of the Lasso estimator.

Lemma C.2. Let x ∈ C 1-var L ([0, 1], R d ).
Then conditionally on A ξ (δ), for a given signature layer k ≥ 1, the maximum among all signature coefficients and individuals is bounded from above, that is

S D •,[k] ∞ ≤ 1 k! .
Proof. It is well known (see, e.g., [START_REF] Fermanian | Functional linear regression with truncated signatures[END_REF]

, Proposition 3) that if X k is the signature of a path x ∈ C 1-var L ([0, 1], R d ), then X k (R d ) ⊗k ≤ ∥x∥ k 1-var k! .
As a consequence, for every word I of size k, one gets

S I (x)) ≤ ∥x∥ k 1-var k! .
The matrix S D N is constructed by taking signatures of linear interpolations of the X i s normalized by their total variation. It therefore contains only signatures of paths of total variation bounded by 1. Taking the maximum on I ∈ {1, . . . , d} k and individuals i = 1, . . . , n, we get

S D •,[k] ∞ ≤ 1 k! .
This final inequality being stated, we can now go back to the proof of Lemma 3.1. We prove it in full generality for p ≥ 1.

In this proof, we make extensive use of the notations introduced in Subsection C.1 and refer the reader to it if a notation is unclear.

Proof. In all the proof, we place ourselves on the set A ξ (δ) defined by Equation ( 9), which ensures that the matrix S D N , seen as a random quantity, is well defined. Recall that we have two sources of randomness: the feature noises ξ i t on the X i s and the target noises ε i t on the Y i s. The feature noises appear only in S D N and make it a random quantity. For S D N to be well-defined, we then need the total variation of the linear interpolation of the feature time series X i to be finite. This holds on the set A ξ (δ) since all noises are then bounded.

Recall that we have defined θ N,M as θ N,M ∈ arg min

θ∈R s d (N )×p 1 2M Y -S D N θ 2 F + Ω(θ).
Note that

1 2M Y -S D N θ 2 F + Ω(θ) = p ℓ=1 1 2M Y ℓ -S D N θ [•],ℓ
where Y ℓ ∈ R M is the ℓ-th column of the target measurement matrix defined in Equation ( 21). The quantity θ

[•],ℓ ∈ R s d (N )
is the ℓ-th column of the parameter matrix defined in Equation ( 19).

By definition, for any θ ∈ R s d (N ) , we have

Y ℓ -S D N θ [•],ℓ 2 2 ≤ Y ℓ -S D N θ [•],ℓ 2 2 + Ω(θ [•],ℓ ) -Ω( θ [•],ℓ ). Moreover, letting ε ℓ = (ε 1 ℓ, t1 1 , . . . , ε n ℓ, tn mn
) ⊤ ∈ R M be a vector of i.i.d. noises (see Equation ( 21)), we have Y ℓ = y ℓ + ε ℓ .

The Pythagorean theorem then yields for any θ ∈ R s d (N ) ,

Y ℓ -S D N θ 2 2 = y ℓ -S D N θ 2 2 + ∥ε ℓ ∥ 2 + 2⟨ε ℓ , y ℓ -S D N θ⟩.
Applying this equation to θ [•],ℓ and θ [•],ℓ , we obtain

1 2M y ℓ -S D N θ [•],ℓ 2 2 ≤ 1 2M y ℓ -S D N θ [•],ℓ 2 2 + 1 M ⟨ε ℓ , S D N ( θ [•],ℓ -θ [•],ℓ )⟩ + Ω(θ [•],ℓ ) -Ω( θ [•],ℓ ). ( 22
)
We now work at each layer of the signature matrix S D N . Towards that end, we rewrite

S D N θ [•],ℓ -θ [•],ℓ = N k=0 S D •,[k] θ [k],ℓ -θ [k],ℓ ,
and bound

ε ℓ , S D N ( θ [•],ℓ -θ [•],ℓ ) = N k=0 ε ℓ , S D •,[k] ( θ [k],ℓ) -θ [k],ℓ ) ≤ N k=0 ∥ε ⊤ ℓ S D •,[k] ∥ ∞ ∥ θ [k],ℓ -θ [k],ℓ ∥ 1 by ℓ 1 -ℓ ∞ norms duality. We fix k and study the term ∥ε ⊤ ℓ S D •,[k] ∥ ∞ . Lemma C.
2 ensures that each of the words of the signature layer of depth k is bounded by 1/k!. As a consequence, by Lemma B.7, under Assumption 6, every element of the vector ε

⊤ ℓ S D •,[k] is v ε M/k! 2 -subgaussian.
It follows that, for any real number µ > 0,

P ∥ε ⊤ ℓ S D •,[k] ∥ ∞ > µ ≤ 2d k exp - (k!) 2 µ 2 v ε M .
We furthermore place ourselves on A ε ( δ) defined by

A ε ( δ) = p ℓ=1 N k=0 ∥ε ⊤ ℓ S D •,[k] ∥ ∞ ≤ 1 k! v ε M log(2pN d k / δ) .
We have just seen that, under Assumption 6 (and still conditionally on A ξ (δ)), one has P(A ε ( δ)) ≥ 1 -δ. Putting together all terms in Equation ( 22) and plugging the definition of Ω given in Equation ( 11), we obtain that, on the set

A ε ( δ) ∩ A ξ (δ), for all θ ∈ R s d (N )×p , 1 2M y -S D N θ N,M 2 F ≤ 1 2M y -S D N θ 2 F + p ℓ=1 N k=0 1 M ∥ε ⊤ ℓ S D •,[k] ∥ ∞ ∥ θ [k],ℓ -θ [k],ℓ ∥ 1 + C k ( δ) k! √ M ∥θ [k],ℓ ∥ 1 -∥ θ [k],ℓ ∥ 1 ≤ 1 2M y -S D N θ 2 F + p ℓ=1 N k=0 1 k! √ M v ε log(2pN d k / δ) ∥ θ [k],ℓ -θ [k],ℓ ∥ 1 + ∥θ [k],ℓ ∥ 1 -∥ θ [k],ℓ ∥ 1 .
Choosing θ = θ * N , by the triangular inequality,

∥ θ [k],ℓ -θ * [k],ℓ ∥ 1 + ∥θ * [k],ℓ ∥ 1 -∥ θ [k],ℓ ∥ 1 ≤ 2∥θ * [k],ℓ ∥ 1 ,
which finally gives us

1 2M y -S D N θ N,M 2 F ≤ 1 2M y -S D N θ * N 2 F + 2 √ M v ε log(2pN d N / δ) p ℓ=1 N k=0 ∥θ * [k],ℓ ∥ 1 k! ≤ 1 2M y -S D N θ * N 2 F + 2p √ M v ε log(2pN d N / δ) N k=0 d k Λ k (F) k! = 1 2M y -S D N θ * N 2 F + 2pC N ( δ) √ M N k=0 d k Λ k (F) k! ,
where the second inequality comes from Proposition C.1. To conclude the proof, we just need to compute the probability of the set A ξ (δ) ∩ A ε ( δ). It is an immediate consequence of Lemma B.6 that P(A ξ (δ)) ≥ 1δ, and we have seen that

P(A ε ( δ)|A ξ (δ)) ≥ 1 -δ, which yields that P(A ξ (δ) ∩ A ε ( δ)) ≥ (1 -δ)(1 -δ). C.3. Proof of Lemma 3.2
This proof relies on bouding the remainder of the Taylor expansion of the CDE.

Proof. For every i = 1, . . . , n and a given point t i ∈ Di , one has, using the upper bound of the approximation error of a CDE by its Taylor expansion provided by Fermanian et al. (2021, Proposition 4)

y i ti -S N (x i [0,ti] )θ * N ≤ d N +1 Λ N +1 (F) (N + 1)! .
This immediately gives

1 M ∥y -S N θ * N ∥ 2 F = 1 M n i=1 ti∈ Di y i ti -S N (x i [0,ti] )θ * N 2 ≤ 1 M M i=1 d N +1 Λ N +1 (F) (N + 1)! 2 = d N +1 Λ N +1 (F) (N + 1)! 2 ,
which concludes the proof.

C.4. A layer-wise bound on the signature

We now prove that signature layers are locally Lipschitz mappings. We start with the following proposition.

Proposition C.3. Let x ∈ C 1-var L ([0, 1], R d ).
Then for all t ∈ [0, 1], the path t → X k [0,t] has 1-variation bounded by

X k 1-var,[0,t] ≤ L k k! .
Proof. By definition of the total variation,

X k 1-var,[0,t] = sup D m i=1 X k [0,ti+1] -X k [0,ti] (R d ) ⊗k = sup D m i=1 X k [ti,ti+1] (R d ) ⊗k , since X k [0,t] = t 0 dx u1 ⊗ • • • ⊗ dx u k ,
and where the supremum is taken over finite dissections D = {0 = t 1 , . . . , t m = 1} of [0, 1]. Notice that the signature layer of depth k is here written as an element of (R d ) ⊗k , which is more convenient for this proof. Then

sup D m i=1 X k [ti,ti+1] (R d ) ⊗k ≤ sup D m i=1 ∥x∥ k 1-var,[ti,ti+1] k! ≤ 1 k! sup D m i=1 ∥x∥ 1-var,[ti,ti+1] k = 1 k! sup D ∥x∥ k 1-var,[0,1] ≤ L k k! ,
where the second inequality follows from the multinomial theorem and the last equality comes from the fact that for all s < u < t, ∥x∥ 1-var,[s,u] + ∥x∥ 1-var,[u,t] = ∥x∥ 1-var, [s,t] . This ends our proof.

Note that this inequality implies that if z is chosen as the linear interpolation of a discretization of x on a grid D, and if the grid gets finer, all signature layers converge at speed ∥x -z∥ ∞,[0,t] but the multiplicative constant increases with depth (if L ≥ 1). Figure 4 illustrates this phenomenon. Figure 4: Difference between the signature of a continuous path x and the signature of its discretized and noisy counterpart X, without noise on the discretization points (left), with noise of variance v ξ = 0.08 2 (middle) and with noise of variance v ξ = 0.5 2 . For every number of sampling points, we average the distance between the two signature over 50 randomly chosen discretizations of the interval [0, 1]. The discretized path is generated as in the well-specified setting (see Appendix D.4).

C.5. Proof of Lemma 3.3

First, recall that for a generic path x : [0, 1] → R d , a modulus of continuity is a continuous function ω x : R ≥0 → R ≥0 vanishing at 0 such that for all s, t ∈

[0, 1] ∥x t -x s ∥ ≤ ω x (|t -s|).
Also recall that by Heine's theorem, we can define such a modulus of continuity for every continuous mapping [0, 1] to R d .

We start by giving a general lemma that bounds the difference between the signature layers of a path and its discretized version. Its proof is based on the results of the previous section.

Lemma C.5. Let x ∈ C 1-var L ([0, 1], R d ) and ω x : R ≥0 → R ≥0 its modulus of continuity. Let x D : [0, 1] → R d be the path obtained by linear interpolation of the discretization of x on a grid D corrupted by additive noise ξ.

Let X k [0,t] and X k,D [0,t]
be their respective k-th layers of signature on [0, t]. Then for all k ≥ 2

X k -X k,D ∞,[0,1] ≤ 2L k-1 k-1 j=1 1 j! max 0≤s≤|D| ω x (s) + max t∈D ∥ξ t ∥ ,
and for k = 1

X 1 -X 1,D ∞,[0,1] ≤ 2 max 0≤s≤|D| ω x (s) + max t∈D ∥ξ t ∥ . Proof. Theorem C.4 yields for k ≥ 2 X k -X k,D ∞,[0,1] ≤ 2L k-1 k-1 j=1 1 j! x -x D ∞,[0,t] Now, remark that x -x D ∞,[0,1] ≤ ∥x -x∥ ∞,[0,1] + max t∈D ∥ξ t ∥ (23)
from the triangular inequality, where x is the piecewise linear path obtained by linear interpolation of x 0 , x t1 , . . . , x tj . Now, since the paths x and x coincide on 0, t 1 , . . . , t j , we have

∥x -x∥ ∞,[0,1] = max i=0,...,j-1 ∥x -x∥ ∞,[ti,ti+1] ≤ max i=0,...,j-1 ω x |t i+1 -t i | = max 0≤s≤|D| ω x (s).
This gives us

X k -X k,D ∞,[0,1] ≤ 2L k-1 k-1 j=1 1 j! max 0≤s≤|D| ω x (s) + max t∈D ∥ξ t ∥ .
For the case k = 1, we immediately get

X 1 -X 1,D ∞,[0,1] ≤ 2 max 0≤s≤|D| ω x (s) + max t∈D ∥ξ t ∥
using the same technique as above.

This result is illustrated in Figure 4. One can notice that as predicted by our theoretical bounds, the convergence of signature of high order happens at the same rate than the convergence of signatures of lower order. However, the multiplicative constant controlling the tightness of the bound increases with N , leading to a slower convergence when N increases. Strong noise hinders the convergence of the signature of the discretized path since in this case, the noise's variance is independent of the number of sampling points : adding more sampling points means adding more noise. There are therefore two tradeoffs when learning with signatures. A first trade-off is between sampling frequency and order: with paths sampled at low resolution, one should prefer lower order signatures, which trade model complexity against precise features. A second trade-off is between sampling and noise: if the feature time series are very noisy, the precision of the features increases up to a certain point, past which noise prevails.

With this result in hand, we can now prove Lemma 3.3.

Proof. We restrict ourselves to the ω-Lipschitz case.

In our setup, after linearly interpolation the time series to obtain x D , we normalize it by its total variation x D 1-var,[0,1] , which is a standard practice when learning with signatures (Morrill et al., 2020a). This means that we compute the signature of the path

1 ∥x D ∥ 1-var,[0,tj ] x D . (24)
Theorem C.4 gets us for k ≥ 2

X k -X k,D ∞,[0,1] ≤ 2L k-1 k-1 j=1 1 j! x - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,1] (25)
≤ 2L k-1 k-1 j=1 1 j! x -x D ∞,[0,1] + 2L k-1 k-1 j=1 1 j! x D - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,1] . (26)
The first term can be bounded by using the fact that in our setting, ω x (s) = ωs, and we thus get

2L k-1 k-1 j=1 1 j! x -x D ∞,[0,1] ≤ 2L k-1 k-1 j=0 1 j! ω|D| + max t∈D ∥ξ t ∥ ≤ 2L k-1 e ω|D| + max t∈D ∥ξ t ∥ (27)
The second term can be bounded by

2L k-1 k-1 j=1 1 j! x D - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,t] ≤ 2L k-1 k-1 j=1 1 j! 1 - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,t] (28)
≤ 2L k-1 k-1 j=1 1 j! 1 - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,t] . (29)
In order to bound

1 - 1 ∥x D ∥ 1-var,[0,1] = x D 1-var,[0,1] -1 ∥x D ∥ 1-var,[0,1] ,
we need both an upper and a lower bound on

x D 1-var,[0,1] . Remark that x D 1-var,[0,tj ] = tu,tu-1∈D x tu + ξ tu -x tu-1 -ξ tu-1 . (30)
Recall that we assume the path (x t ) to be time-augmented, and that the measurement times are not noisy. This means that tu,tu-1∈D

x tu ξ tu -x tu-1 -ξ tu-1 ≥ tu,tu-1∈D |t u -t u-1 | ≥ t 2 -t 1 = t 2 (31)
since t 1 = 0 and Assumption 4 guarantees that there are at least two sampling points in every grid. This gives us that

1 ∥x D ∥ 1-var,[0,1] ≤ 1 t 2 ≤ 1 η , (32)
since we have required that the last sampling time is at least η in Assumption 4. Turning to the upper bound, we get that

1 -x D 1-var,[0,1] ≤ 1 -L + tu,tu-1∈D ∥ξ u -ξ u-1 ∥
by definition of the total variation of a piecewise linear path. Finally,

tu,tu-1∈D ∥ξ u -ξ u-1 ∥ ≤ 2#D max t∈D ∥ξ t ∥ . (33)
Putting everything together gives us

1 - 1 ∥x D ∥ 1-var,[0,1] ≤ 1 -L + 2#D max t∈D ∥ξ t ∥ η . (34)
In the end, we get that

2L k-1 k-1 j=1 1 j! 1 - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,1] ≤ 2L k-1 e 1 -L + 2#D max t∈D ∥ξ t ∥ η x D ∞,[0,1] . (35)
Now, remark that since the path x D is piecewise linear,

x D ∞,[0,1] = x D -x 0 + x 0 ∞,[0,1] ≤ max t∈D ∥x t + ξ t -x 0 ∥ + ∥x 0 ∥ (36) ≤ max t∈D ∥x t -x 0 ∥ + max t∈D ∥ξ t ∥ + ∥x 0 ∥ (37)
≤ ∥x 0 ∥ + L + max t∈D ∥ξ t ∥ (38)
where the inequality max

t∈D ∥x t -x 0 ∥ ≤ L
follows from the definition of the total variation. This means that

2L k-1 k-1 j=1 1 j! 1 - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,1] ≤ 2L k-1 e 1 -L + 2#D max t∈D ∥ξ t ∥ η ∥x 0 ∥ + L + max t∈D ∥ξ t ∥ . (39)
We have written these inequalities for a generic random variable. Let us now consider individual observations of our dataset.

On the set A ξ (δ), one has

max i=1,...,n,t∈D i ξ i t ≤ v ξ √ d + v ξ c -1 log(δ -1 #D). (40)
To simplify notations, let us write

C δ := v ξ √ d + v ξ c -1 log(δ -1 #D). (41)
We get

2L k-1 k-1 j=1 1 j! x i -x i,D ∞,[0,1] ≤ 2L k-1 e ω|D| + C δ , (42)
where we recall that D is the collection of individual grids, and |D| is the biggest sampling gap among individuals. Similarly,

2L k-1 k-1 j=1 1 j! 1 - 1 ∥x D ∥ 1-var,[0,1] x D ∞,[0,1] ≤ 2L k-1 e 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ (43)
Now moving to the feature matrices, we have

1 M (S N -S D N )θ * N 2 F ≤ 1 M n i=1 N k=0 (S i,[k] -S D i,[k] )θ * [k],• 2 F ≤ 1 M n i=1 t∈ Di N k=0 d k Λ k (F) 2 2eL k-1 ω|D| + C δ + 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ 2 ≤ 4e 2 ω|D| + C δ + 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ 2 L 2 N k=0 d k Λ k (F) 2 k! 2 × k! 2 ≤ 4e 2 N ! 2 ω|D| + C δ + 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ 2 L 2 N k=0 d k Λ k (F) 2 k! 2 . Writing C D,N (δ) = 4e 2 L 2 N ! 2 ω|D| + C δ + 1 -L + 2#DC δ η ∥x 0 ∥ + L + C δ 2 ,
one finally gets with probability 1δ that

1 M (S N -S D N )θ * N 2 F ≤ C D,N (δ) N k=0 d k Λ k (F) 2 k! 2 .
Dimension p of the target path. Our oracle bound only depends on p through the statistical error term. This term is proportional to p √ log p, which is expected in multitask regression.

Dimension d of the feature path. Our oracle bound exhibits multiple dependencies in d. First, the truncation bias grows polynomially with d. Similarly, the discretization bias also depends polynomially on d. Finally, the statistical error term is proportional to log d times a polynomial term.

D. Algorithms, experiments, and supplementary results

D.1. Implementation details

Recall that the SigLasso estimator θ N,M is defined as

θ N,M ∈ arg min θ∈R s d (N )×p 1 2M Y -S D N θ 2 F + Ω(θ),
where

Ω(θ) = N k=0 M 1 2 C k ( δ) k! θ [k],• 1 , and
C k ( δ) = v ε log(2pN d k / δ) for δ ∈ [0, 1].
Our goal is first to rewrite the penalty Ω(θ) as

Ω(θ) = C N k=0 λ k θ [k],
• 1 , such that training will only require to scale each layer of θ and to crossvalidate the multiplicative constant C. Since for k ≥ 1,

C k ( δ) = √ k × v ε ( δ/k + log(pN )/k + log d) ≤ √ k × v ε ( δ + log(pN ) + log d),
we let λ k = √ k k! . We now show that the minimization problem with layer-specific penalty can be written as a standard regression problem with ℓ 1 penalization by rescaling the feature matrix, that is, multiplying S D N by a well-chosen diagonal matrix. Consider the ℓ 1 -penalized problem

min θ∈R s d (N )×p 1 2M Y -S D N θ 2 F + C N k=0 λ k θ [k],• 1 ,
where C > 0 controls the strength of the penalization.

Making the change of variable

θ = diag(1,
min θ∈R s d (N )×p 1 2M Y -S D N W θ 2 F + C N k=0 θ[k],• 1 .
We can thus obtain the SigLasso estimator by (i) multiplying the feature matrix S D N by W and solving the associated ℓ 1 -penalized problem (ii) multiplying the obtained solution by W .

The Learn-And-Reconstruct algorithm is the generic algorithm used in our work. It is applicable for a wide variety of tasks such as missing values inference, trajectory reconstruction, forecasting and many more. It is described in Algorithm 1. We two metrics used to asses the importance of the different dimensions of the feature path.

Given a truncation depth N and a dimension i ∈ {1, . . . , d}, we define its pure feature importance (PFI) as the sum of the norm of the coefficients (or vectors in the case the target is multivariate) of θ N,M that are associated to signatures taken on the words I 1 = (i), I 2 = (i, i), and so forth until I N = (i, . . . , i). Mathematically,

P F I(i) = 1 N θ I1 2 + θ I2 2 + • • • + θ I N 2 .
Since signatures also capture interactions between dimensions of the feature path, we also define the cross feature importance (CFI) as the sum of norms of the coefficients (or vectors) of θ N,M that are associated to signatures coefficients of words of length ≤ N in which the letter i appears. Mathematically,

CF I(i) = 1 s d (N ) -s d-1 (N ) I s.t. i ∈ I θ I 2 .
For a given truncation depth N , note that there are

s d (N ) -s d-1 (N ) = N k=0 d k - N k=0
(d -1) k terms in the last sum, which justifies our choice of normalization.

D.3. Details on model implementation and evaluation

SigLasso. The SigLasso model is implemented using the CVLasso class in scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. This implementation optimises the objective function using coordinate descent and features automatic cross-validation of the penalty strength. We use iisignature [START_REF] Reizenstein | Algorithm 1004: The iisignature library: Efficient calculation of iteratedintegral signatures and log signatures[END_REF] to compute the signature of the feature time series. Every time series is standardized prior to this through division by its own total variation, as suggested by Morrill et al. (2020a). The depth of the signature is a hyperparameter chosen between 2 and 9 or 6 depending on the experiment. An intercept is added.

GRU. The GRU is of width 128 and systematically trained with 100 epoches using a learning rate of 0.001.

Neural CDE. We use the implementation of Neural CDE provided by torchcde [START_REF] Kidger | Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu[END_REF]. We use the original vector field described in the documentation of this package, with the small tweak that we use a smoother nonlinearity (tanh instead of ReLU). We observed that using the rk4 solver instead of dopri5 significantly accelerates the training time of the Neural CDE without affecting the model's performances. The learning rate is hand-tuned to either 0.001 or 0.0001 depending on the experiment. We train the model for 100 epochs and asses its convergence by using a standard stopping criteria.

Metrics. The MSE is computed in a classical fashion. To compute the integrate MSE, we compute the L 2 distance between the piecewise constant interpolations of the true y t and the predicted y t .

D.4. Details on the well specified model

Generation of the training data. We generate a two-dimensional feature path by interpolating for every dimension 15 points in [0, 1], each of them being draw randomly for a normal distribution N (0, 1). The interpolation is done with Hermite cubic splines with backward differences using the package torchcde [START_REF] Kidger | Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu[END_REF]. Time is added as a supplementary channel, which is a standard practice when learning with signatures and Neural CDEs. These paths are then downsampled by randomly drawing sampling points for the target and the feature time series specific to every individual. The target path is the solution of a CDE of the form

dy t = σ Ay t dx t
where σ : R p → R p is the hyperbolic tangent, A ∈ R d×p is a matrix drawn randomly from N (0, I d×p ) and x t is the feature path constructed as above. The solution of this CDE is computed using torchcde.

Generation of the test data. We generate the test data in the same way than the training data. However, this data is not downsampled as we wish to assess the generalization capacities of our model-i.e., is our model capable of approximating the dynamics and extrapolating to continuous feature paths.

D.5. Details on Ornstein-Uhlenbeck experiment

We take (x t ) t∈[0,1] to be a 1-dimensional Brownian motion with variance σ 2 = 0.1, and generate (y t ) as a 1-dimensional Ornstein-Uhlenbeck process driven by (x t ), that is, for all t ∈ [0, 1]

dy t = θ(µ -y t )dt + dx t .
Simulation of (y t ) is done using a standard Euler-Maruyama simulation scheme. We let θ = 3 and µ = 1. The training data is then downsampled as in the well specified experiment.

D.6. Details on the tumor growth experiment

We consider the following tumor growth model taken from [START_REF] Simeoni | Predictive pharmacokineticpharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents[END_REF]. Let x ∈ C 1-var ([0, 1], R). The weight y ∈ C([0, 1], R + ) under the concentration of a treatment drug x is governed by the differential system

du 1 t = λ 0 u 1 t 1 + ( λ 0 λ 1 y t ) ψ -1/ψ -k 2 x t u 1 t dt du 2 t = k 2 x t u 1 t -k 1 u 2 t dt du 3 t = k 1 (u 2 t -u 3 t ) dt du 4 t = k 1 (u 3 t -u 4 t ) dt y t = u 1 t + u 2 t + u 3 t + u 4 t
with initial condition (u 1 0 , u 2 0 , u 3 0 , u 4 0 , y 0 ) = (2, 0, 0, 0, 2) and parameters (k 1 , k 2 , λ 0 , λ 1 , ψ) = (10, 0.5, 0.9, 0.7, 20). The concentration (x t ) is chosen to be the squared value of the paths used for the well-specified experiment. Notice that this system is non-linear w.r.t. x. Indeed, writing dy t = G(y t , x t )dt, one has, for α ∈ R, G(y t , αx t ) ̸ = αG(y t , x t ). The training data is then downsampled as in the well specified experiment.

D.7. Supplementary results

Table 2: Performance of SigLasso, GRU, Neural CDE, RNN and LSTM in different simulation settings, averaged over 10 iterations. In every setting, n = 50, # Di = 5 for all i = 1, . . . , n (and therefore M = 250). We illustrate the performance of our method and competitors on French Covid data from 2021-03-31 to 2021-07-07 available on Gitlab. Hospital data was obtained from the SI-VIC database, the national inpatient surveillance system.

Target path. Following [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF], we chose the predict the growth rate of incident hospitalisations in each of the 9 metropolitan regions of France. The exponential growth rate was computed from raw data using a 2 days rolling window and then smoothed using local polynomial regression as in [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF]. Mathematically, our target time series is the R-valued growth rate, and we fit a different model for every of the 12 regions. It is displayed for all 12 regions in Figure 6.

Feature path. As in [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF], we consider a set of 12 time-dependant predictors of different types summarized in Table 3 and plotted in Figure 5. Both SIDEP ("Système d'Information de Dépistage Populationnel") and VAC-SI datasets are publicly available. The mobility data was obtained from Google. The mobility-related predictors describe travel trends for different kind of public spaces such as such as shops and leisure spaces, food stores and pharmacies, parks, public transport stations, workplaces and residential areas. The meteorological data was obtained from Météo France.

Models. SigLasso, NCDE and GRU algorithms were trained on the period from 2021-03-31 to 2021-06-23 and tested on the period from 2021-06-24 to 2021-07-07. We included a history of 10 days at each point and performed prediction for different horizons ranging from 1 to 14. In others words, at horizon h, features values from day th -10 to day th to were used to compute the prediction at time t. All feature time series are normalized to have total variation equal to 1. 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 Figure 11: RMSE accross all regions on the training period (left) and the testing period (right) for the ensemble methode [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF], NCDE, GRU, and SigLasso. See Figure 12 for a zoom-in on GRU and SigLasso performances. Moreover, we show in Figure 13 the results of the L 2 reconstruction error in the well-specified setting, when we vary the number of sampling points of the feature paths between 10 and 10 3 . We see that SigLasso always outperforms GRU and Neural CDE but that the difference of performance is more important when there are only a few sampling points. In this regime SigLasso is moreover more stable.

0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 Growth rate ARA 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 Days - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 GES 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 NAQ 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 Growth rate BFC 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 Days - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 HDF 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 OCC 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 Growth rate BRE 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 Days - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 IDF 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 PDL 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 Growth rate CVL 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 Days - 0 .1 0 - 0 .0 5 0 .0 0 0 .0 5 0 .1 0 NOR 2 0 2 1 -0 4 -0 7 2 0 2 1 -0 4 -1 9 2 0 2 1 -0 5 -0 1 2 0 2 1 -0 5 -1 3 2 0 2 1 -0 5 -2 5 2 0 2 1 -0 6 -0 6 2 0 2 1 -0 6 -1 8 2 0 2 1 -0 6 -3 0 - 0 .1 0 - 0 .0 5

  Proceedings of the 40 th International Conference on Machine Learning, Honolulu, Hawaii, USA. PMLR 202, 2023. Copyright 2023 by the author(s).

Figure 1 :

 1 Figure 1: The workflow of our model. Starting from the left, the first panel describes our modelling hypothesis: the target and feature time series are linked through an unknown CDE. The second panel shows the observed data. The third panel shows how every observation of the target is mapped to the signature of the corresponding path, how to construct the dataset S D N , Y from this data and finally how to learn the SigLasso estimator.

Figure 2 :

 2 Figure 2: L 2 reconstruction error of SigLasso, GRU and Neural CDE in the well-specified setting, for varying number of target samples.

Figure 3 :

 3 Figure3: Interpolation (left of dotted line) and prediction (right of dotted line) of HGR in region Île de France for Si-gLasso, GRU and NCDE. The lighter the blue, the smaller the horizon h. Ground truth is in red. NCDE overfits and is unable to predict the HGR during the test period.

Fη

  (x [0,t] ) Truncated signature at order N s d (N ) Size of the signature truncated at order N of a d dimensional path -ȳN,t Taylor expansion of the solution of a CDE of order N evaluated in t Section 2, page 4, Definition 2.2 Φ I Differential product of the vector field F along I Appendix B.3, page 15, Definition B.3 θ ⋆ N Matrix collecting all differential products up to order N Section 2, page 5, Equation (5) Λ k (F) Norm on the differential product of F Section 2, page 5, Assumption 3 Minimal sampling time Section 2, page 5, Assumption 4 D Set of individual specific sampling grids of features Section 2, page 5 D Set of individual specific sampling grids of targets -|D| Meshisze of a sampling grid D -#D Number of sampling points in a sampling grid -Y Matrix collecting all target measurements of the sample -

Proposition B. 2 (

 2 Integration by parts). Let x, y ∈ C 1-var L ([0, 1], R d ). Then t s y u dx u + t sx u dy u = y t x ty s x s SeeFriz & Victoir (2010, Proposition 2.4) for a proof.

  , a N ) and b = (b 0 , . . . , b N ) both in T N (R d ), we define a ⊗ b = (c 0 , . . . , c j , . . . , c N ), where c j = j k=0 a k ⊗ b j-k .

  |A ij | and ∥A∥ F = i=1,...,n j=1,...,p |A ij | 2 .

Figure 6 :

 6 Figure 6: Hospitalization growth rate through time during the full period for the 12 different regions of France.

Figure 8 :

 8 Figure 8: Interpolation (left of dotted line) and prediction (right of dotted line) of hospitalization growth rate for all 12 french regions using NCDE.

Figure 9 :

 9 Figure 9: Interpolation (left of dotted line) and prediction (right of dotted line) of hospitalization growth rate for all 12 french regions using GRU.

PACFigure 10 :

 10 Figure10: Interpolation (left of dotted line) and prediction (right of dotted line) of hospitalization growth rate for all 12 french regions using ensemble methods[START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF].

Figure 12 :

 12 Figure12: RMSE accross all regions on the training period (left) and the testing period (right) for GRU and SigLasso.

Figure 13 :

 13 Figure 13: L 2 reconstruction error of SigLasso, GRU and Neural CDE in the well-specified setting, for varying number of feature samples.

Table 1 :

 1 Performance of SigLasso, GRU and Neural CDE in different simulation settings, averaged over 10 iterations. In every setting, n = 50, # Di = 5 for all i = 1, . . . , n (and therefore M = 250). term ω|D| refers to the longest time between two observations amongst individuals. Not sampling an individual during a long period of time causes a loss in information, which is bounded by the Lipschitz control ω of the feature path.

			L 2 error			MSE on last point
	Setting	SigLasso	GRU	Neural CDE	SigLasso	GRU	Neural CDE
	Well-specified	0.13 ± 0.07	1.05 ± 0.42	0.61 ± 0.38	0.73 ± 0.56	3.32 ± 1.60	1.46 ± 1.20
	Ill-specified	0.15 ± 0.02	0.24 ± 0.11	0.29 ± 0.15	0.09 ± 0.05	0.19 ± 0.09	0.22 ± 0.15
	OU	0.01 ± 0.02	0.05 ± 0.06	0.17 ± 0.12	0.018 ± 0.025	0.014 ± 0.020	0.013 ± 0.016
	Tumor growth	0.16 ± 0.02	0.66 ± 0.09	5.29 ± 1.38	0.35 ± 0.12	2.00 ± 0.38	8.76 ± 9.26

the

  λ 1 , . . . , λ 1 , . . . , 1/λ 1 , . . . , 1/λ k , . . . , 1/λ k ) θ, and denoting by W this last weight matrix, we get the equivalent minimization problem

		, λ 2 , . . . , λ 2	, . . . , λ k , . . . , λ k	)θ,
	d repetitions	d 2 repetitions	d N repetitions
	which is equivalent to
	θ = diag(1, 1/λ 1

  Algorithm 1 Learn-and-Reconstruct Algorithm. The algorithm infers for every individual in the test set a reconstructed time series Ŷ i t . 1. Learn the dynamics Input: train dataset of normalized paths (X 1 , Y 1 ), . . . , (X n , Y n ) sampled on (D 1 , D1 ), . . . , (D n , Dn ). test dataset X1 , . . . , Xn sampled on D1 , . . . , Dn .

	Construct the feature matrix S D N and the target vector Y
	for i = 1 to n do for t in Di do
	S D N ← Append S N (X i [0,t] ) t Y ← Append Y i end for
	end for Compute θN,M by solving (7) with Y, S D N using coordinate descent.
	2. Reconstruct trajectories
	Input: for i = 1 to n do for t in Di do Ŷ i t = θN,M S N ( Xi [0,t] )
	end for
	end for
	D.2. Assessing feature importance

Table 3 :

 3 The set of time-dependant predictors used to predict the hospital admission growth rate

			14										0.40
								0.10
													0.35
			12
													0.30
								0.05
	ARA BFC BRE	abs_humidity	8 10				P_r	0.00				couv_complet	0.20 0.25
	CVL												0.15
	GES
	HDF							0.05
	IDF												0.10
	NOR		6
	NAQ
	OCC												0.05
	PDL							0.10
	PAC
			4
			2021-03-14	2021-04-13	2021-05-13 Date	2021-06-12	2021-07-12	2021-03-14	2021-04-13	2021-05-13 Date	2021-06-12	2021-07-12	2021-03-14	2021-04-13	2021-05-13 Date	2021-06-12	2021-07-12
	Figure 5: The 12 different feature time series used to forecast the hospitalization growth rate. Every different color
	corresponds to a given region of France.

introduced in Lemma

3.2, vanishes as N grows.

+ Ω(θ [•],ℓ ),

Acknowledgement. We thank anonymous ICML reviewers for their remarks, which helped us improve this paper. LB thanks Gérard Biau and Claire Boyer for supervising a previous internship which sparked his interest in signatures. LB and AF thank the Sorbonne Center for Artificial Intelligence (SCAI) and its team.

We now state a bound on the difference between the k-th layer of the signatures of two different paths. Theorem C.4. Let x, z ∈ C 1-var L ([0, 1], R d ). Then for all k ≥ 2, the difference in supremum norm between the paths t → X k [0,t] and t → Z k [0,t] is bounded by

Proof. Our proof works by induction. Let x, z ∈ C 1-var L ([0, 1], R d ), and for t ∈ [0, 1] denote by X k [0,t] (resp. Z k [0,t] ) the k-th layer of the signature of x (resp. z). For k = 1 and t ∈ [0, 1], remark that

Consider now k ≥ 2. We have

and thus

We now bound each of these terms separately. First,

Moving to the first integral, integration by parts yields

We stress that Proposition (B.2) applies since the integral over the tensor product is taken coordinate-wise. Since X k-1 [0,0] = 0, we are left with

Using Lemma C.3 and submultiplicativity of the tensor norms, this can thus be bounded by

Finally, we are left with

which can be recursively bounded by

C.6. Proof of the main Theorem

We finally combine all Lemmas to obtain the desired oracle bound.

Proof. First, we have from Lemma 3.1 that on A ε ( δ),

The first term of the right-hand side of this inequality is bounded by

By Lemma 3.2 and Lemma 3.3, this can in turn be bounded on

Combining all the pieces, this finally gives us, on

C.7. Asymptotics

We briefly discuss the asymptotic behaviour of the upper bound of the oracle inequality.

Truncation depth N . A natural question is whether the bias of our estimator vanishes as N → ∞. If we have perfect sampling, i.e. the limit case where D = 0 and v ξ = 0, our bound on the prediction error becomes on A ε ( δ)

The first term of this bound vanishes as an immediate consequence of Assumption 3, while the second term is a statistical error term that behaves like

. In order to obtain an asymptotic convergence, we thus need that N log(dN ) = o(M ).

In the more realistic setting where |D| > 0, the discretization bias behaves like L N -1 N !|D|. It is thus sufficient to assume that |D| = o(1/N !). If v ξ > 0, our estimator is durably biased due to the measurement noise, and this bias increases with N → ∞. This is due to a "propagation of chaos" phenomenon: the difference between the unobserved feature path and the interpolated feature time series is amplified by taking the successive iterated integrals that define the signature. This advocates for using simple, low-order signature models in the presence of noise, as the gain in precision obtained when taking higher N and reducing the truncation bias will at some point be lost because of the amplified noise. Architectural details. The GRU has width 128 and is trained for 100 epochs with a learning rate of 0.0001. The NCDE is trained for 30 epochs with a learning rate of 0.001. It has 2 hidden layers of width 128, an intermediate Tanh(•) nonlinearity and a final linear readout. This architecture is identical to the one proposed in [START_REF] Kidger | On neural differential equations[END_REF]. Penalty strenght of the SigLasso is crossvalided using the internal implementation LassoCV of scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].

All details, in particular the features used for each individual prediction, can be found in [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF].

We refer to the supplementary information file of [START_REF] Paireau | An ensemble model based on early predictors to forecast covid-19 health care demand in france[END_REF] and our code for more details.

Results. Figure 11 displays the RMSE (on all regions) of NCDE, SigLasso, GRU, and the Ensemble method for all prediction horizons h = 1, . . . , 14. Figures 7, 9 and 8 display the obtained interpolation for SigLasso, GRU and NCDE at different horizons (corresponding to different line colors in winter matplotlib palette). The lighter the blue, the smaller the time horizon: the lightest curve corresponds to a time horizon equal to h = 1. Truth is in red.

D.9. Additional results

We give in Table 4 some additional results on the experiments described above.