Quentin Clairon

J-B Brunel

State and Parameter Estimation of Partially Observed Linear Ordinary Dierential Equations with Deterministic Optimal Control

Keywords: Ordinary Dierential Equation, Optimal Control, Parameter Estimation, Smoothing, Riccati Equation, M-estimation. 1

Ordinary Dierential Equations are a simple but powerful framework for modeling complex systems. Parameter estimation from times series can be done by Nonlinear Least Squares (or other classical approaches), but this can give unsatisfactory results because the inverse problem can be ill-posed, even when the dierential equation is linear.

Following recent approaches that use approximate solutions of the ODE model, we propose a new method that converts parameter estimation into an optimal control problem: our objective is to determine a control and a parameter that are as close as possible to the data. We derive then a criterion that makes a balance between discrepancy with data and with the model, and we minimize it by using optimization in functions spaces: our approach is related to the so-called Deterministic Kalman Filtering, but dierent from the usual statistical Kalman ltering.

We show the root-n consistency and asymptotic normality of the estimators for the parameter and for the states. Experiments in a toy model and in a real case shows that our approach is generally more accurate and more reliable than Nonlinear Least Squares and Generalized Smoothing, even in misspecied cases.

Introduction

Ordinary Dierential Equations (ODE) are a widely used class of mathematical models in biology, physics, engineering, . . . Indeed, it is a relatively simple but powerful framework for expressing the main mechanisms and interactions of potentially complex systems. It is often a reference framework in population dynamics and epidemiology [START_REF] Ellner | Dynamic Models in Biology[END_REF], virology [START_REF] Nowak | Virus Dynamics: Mathematical Principles of Immunology and Virology[END_REF], or in genetics for describing gene regulation networks [START_REF] Costello | Wisdom of crowds for robust gene network inference[END_REF][START_REF] Wu | Sparse additive odes for dynamic gene regulatory network modeling[END_REF]. The model takes the form ẋ = f (t, x, θ), where f is a vector eld, x is the state, and θ is a parameter that can be partly known. The parameter θ is often of high interest, as it represents rates of changes, phenomenological constants needed for interpretability and analysis of the system. Typically, θ can be related to the sensitivity of a variable with respect to other variables.

Hence, the parameter estimation of ODEs from experimental data is a long-standing statistical subject that have been adressed with many dierent tools.

Estimation can be done with classical estimators such as Nonlinear Least Squares (NLS) and Maximum Likelihood Estimator (MLE) [START_REF] Li | Parameter estimation of ordinary dierential equations[END_REF][START_REF] Walter | Identication of parametric models[END_REF][START_REF] Pronzato | Optimal experimental design and some related control problems[END_REF] or

Bayesian approaches [START_REF] Huang | A bayesian approach for estimating antiviral ecacy in hiv dynamic models[END_REF][START_REF] Gelman | Physiological pharmacokinetic analysis using population modeling and informative prior distributions[END_REF][START_REF] Ghasemi | Bayesian parameter estimation for nonlinear modelling of biological pathways[END_REF][START_REF] Calderhead | Accelerating bayesian inference over nonlinear dierential equations with gaussian processes[END_REF]. Nevertheless, the statistical estimation of an ODE model by NLS leads to a dicult nonlinear estimation problem. Some diculties were pointed out by Ramsay et al. [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF] such as computational complexity, due to ODE integration and nonlinear optimization. These diculties are in fact reminiscent of intrinsic diculties in the parameter estimation problem, that makes it an ill-posed inverse problem, that needs some regularization [START_REF] Hein W Engl | Inverse problems in systems biology[END_REF][START_REF] Stuart | Inverse problems: A bayesian perspective[END_REF] .

Alternative statistical estimators have been developped to deal with this particular framework, such as Generalized Smoothing [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF][START_REF] Qi | Asymptotic eciency and nite-sample properties of the generalized proling estimation of parameters in ordinary dierential equations[END_REF][START_REF] Hooker | Parameter estimation in dierential equation models with constrained states[END_REF][START_REF] Campbell | Maximum prole likelihood estimation of dierential equation parameters through model based smoothing state estimates[END_REF] or Two-Step estimators [START_REF] Varah | A spline least squares method for numerical parameter estimation in differential equations[END_REF][START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF][START_REF] Liang | Parameter estimation for dierential equation models using a framework of measurement error in regression models[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF][START_REF] Brunel | Parameter estimation of ordinary dierential equations with orthogonality conditions[END_REF]. Two-step estimators use a nonparametric estimator

X and aim at minimizing quantities characterizing the dierential models, such as the weighted L 2 distance ´T 0 ˙ X(t) -f (t, X(t), θ) 2 w(t)dt. These estimators have a good computational eciency as they avoid repeated ODE integration.

In practice, the used criteria are also smoother and easier to optimize than the NLS criterion. Two-step estimators are consistent in general, but there is a trade-o with the statistical precision, and some care in the use of nonparametric estimate ˙ X has to be taken in order to keep a parametric rate [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF].

In the case of Generalized Smoothing [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF], the solution X * is approximated by a basis expansion that solves approximately the ODE model; hence, the parameter inference is performed by dealing with an imperfect model. Based on the Generalized Proling approach, Hooker proposed a criteria that estimates the lack-of-t through the estimation of a forcing function t → u(t) in the ODE ẋ -f (t, x, θ) = u(t), where θ is a previous estimate obtained by Generalized Proling.

In [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF], the authors have proposed a two-step estimator for linear models, that avoids the use of ˙ X and introduces a forcing function without the nite basis decomposition by using control theory. The principle is to transform the estimation problem into a control problem: we have to nd the best (or smallest) control u such that the ODE is close to the data. The limitations of the results provided in [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF] were the restriction to fully observed system with known initial condition. The objective of this paper is to provide a similar twostep estimate that permits the estimation of θ without knowing x 0 , that deals with the partially observed case and provides state estimates.

One interest of the approach used is to deal directly with the optimization in a function space without using of series expansion for function estimation.

Moreover, innite dimensional optimization tools give a powerful characterization of the solutions, useful in practice. This work can be seen as an extension of the previous one [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF], aiming to use control theory result for parameter inference.

We deal now with the partially observed case with unknown initial condition, that gives rise to a methodology close to the so-called Deterministic Kalman

Filter. Indeed, in that paper, we assume that the system is linear, with a linear observation function.

Our method provides a consistent parametric estimator when the model is correct. We show that it is root-n consistent and asymptotically normal. At the same time, we get a discrepancy measure between the model and the data under the form of an optimal control u analogous to the forcing function in [START_REF] Hooker | Forcing function diagnostics for nonlinear dynamics[END_REF],

and we show that we can estimate the nal and initial conditions and hence all the states if needed, in particular the hidden ones.

In the next section, we introduce the notations and we motivate our approach by discussing the Generalized Smoothing approach, and the link with Optimal Control Theory. In section 3, we investigate the existence and regularity of our new criterion; in particular, we derive necessary and sucient conditions for dening our approach in partially observed case. We show that the estimator is consistent under some regularity assumption about the model. Then in section 4, we show that we reach the root-n rate using regression splines for Y the nonparametric estimator of the observed signal. We derive then the consistency of the state estimator derived. Finally, we show the interest of our method on a toy model and in a real model used in chemical engineering, by a comparison with Nonlinear Least Squares and Generalized Smoothing.

Model and methodology

We introduce rst the statistical ODE model of interest, and the basic notations for dening our estimator. We relate this work to the Generalized Smoothing estimator and the Tracking estimator.

Model and Notations

We partially observe a true trajectory X * at random times 0 = t 1 < t

Y i = CX * (t i ) + i
where i is a random noise and C is the observation matrix of size d × d.

We assume that there is a true parameter θ * belonging to a subset Θ of R p , such that X * is the unique solution of the linear ODE

ẋ(t) = A θ (t)x(t) + r θ (t) (1)
with initial condition X * (0) = x * 0 ; where t → A θ (t) ∈ R d×d and t → r θ (t) ∈ R d . More generally, we denote X θ,x0 the solution of (1) for a given θ, and initial condition x 0 . We assume that x * 0 and θ * are unknown, and that they must be estimated from the data (y 1 , . . . , y n ). The parameter θ * is the main parameter of interest, whereas the initial condition is considered as a nuisance parameter, needed essentially for the computation of candidate trajectories X θ,x0 .

For linear equations, a central role is played by the solutions of the homogeneous ODE

ẋ(t) = A θ (t)x(t). (2)
Indeed, for each s in [0, T ], we denote t → Φ θ (t, s) the solution to the matrix ODE (2), with initial condition I d at time s (i.e Φ θ (s, s) = I d ). The function (t, s) → Φ θ (t, s) is a d × d matrix valued function, called the resolvant of the ODE. It permits to give an explicit dependence of the solutions of (1) in r θ and the initial condition x 0 , thanks to Duhamel's formula:

X θ,x0 (t) = Φ θ (t, 0)x 0 + ˆt 0 Φ θ (t, s)r θ (s)ds.
A consistent and classical method for the estimation of θ * is Nonlinear Least Squares (NLS), that minimizes

n i=1 Y i -CX θ,x0 (t i ) 2 2 .
A classical alternative is Generalized Smoothing (GS), that uses approximate solutions of the ODE (1). GS replaces the solutions X θ,x0 by splines that smooth data and solve approximately the ODE with a penalty based on the ODE model.

A basis expansion X(t, θ) = β(θ) T p(t) is computed for each θ, where β(θ) is obtained by minimizing in β the criterion

J n (β|θ, λ) = n i=1 Y i -Cβ T p(t) 2 2 + λ ˆT 0 β T ṗ(t) -A θ (t)β T p(t) + r θ (t) 2 2 dt (3)
This rst step is considered as proling along the nuisance parameter β, whereas the estimation of the parameter of interest is obtained by minimizing the sum of squared errors of the proxy X(t, θ):

θGS = arg min θ n i=1 Y i -C X(t i , θ) 2 (4)
In practice, the hyperparameter λ needs to be selected from the data with adaptive procedures, see [START_REF] Chernovena | Estimation of nonlinear dierential equation for glucose-insulin dynamics in type i diabetic patients using generalized smoothing[END_REF].

The essential dierence with NLS is the replacement of the exact solution X θ,x0 by the approximation X(•, θ) (that depends also on the data). This change induces a new source of error in the estimation of the true trajectory t → X * (t)

as the functions X(•, θ) are splines that do not solve exactly the ODE model

Y i -Cβ T p(t) 2 2
. For this reason, the ODE model ( 1) is not solved and it is useful to introduce the discrepancy term ûθ (t) = β T ṗ(t) -A θ (t)β T p(t) + r θ (t) that corresponds to a model error. In fact, the proxy X(•, θ) satises the perturbed ODE ẋ = A θ x + r θ + ûθ . This forcing function ûθ is an outcome of the optimization process and can be relatively hard to analyze or understand, but its analysis provides a good insight into the relevancy of the model [START_REF] Hooker | Forcing function diagnostics for nonlinear dynamics[END_REF][START_REF] Hooker | Goodness of t in nonlinear dynamics: Mis-specied rates or mis-specied states?[END_REF].

Based on these remarks, we introduce the perturbed linear ODE

ẋ(t) = A θ (t)x(t) + r θ (t) + u(t) (5)
where the function t → u(t) can be any function in L 2 . The solution of the

corresponding Initial Value Problem    ẋ(t) = A θ (t)x(t) + r θ (t) + u(t)
x(0) = x 0 is denoted X θ,x0,u . Instead of using the spline proxy X(•, θ) for approximating X * , we use the trajectories X θ,x0,u of the ODE (5) controlled by the additional functional parameter u.

In [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF], the same perturbed model is introduced but the cost function is simpler as the observation matrix C is the identity, and the initial condition is xed. In that framework, an M-estimator for θ is proposed, based on the optimization of the criterion

S( Y ; x 0 , θ, λ) = inf u∈L 2 { Y -X θ,x0,u 2
L 2 + u 2 L 2 }. (6)
The proper denition of S and the derivation of its properties were obtained by using some classical results of Optimal Control Theory. Essentially, the computation of S corresponds to the classical "tracking problem" that can be solved by the Linear-Quadratic theory (LQ theory). LQ theory solves the minimization problem in L 2 of the cost function

C(u) = X θ,x0,u (t) 2 L 2 + u(t) 2 L 2 + X θ,x0,u (T ) QX θ,x0,u (T ) (7)
The criteria S used for parameter estimation is associated to the value function dened in Optimal Control as S(t, x) = inf{C(u)|X θ,x,u (t) = x}. The value function plays a critical role in the analysis of optimal control problems, typically for the computation of an optimal policy. Under regularity assumptions, the value function S is the solution of the Hamilton-Jacobi-Bellman Equation, which is a rst order Partial Dierential Equation [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF]. Quite remarkably, for a linear ODE with a quadratic cost such as [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF], the value function is a quadratic form in the state x, i.e S(t, x) = -x E(t)x, where E(t) is the solution of a matrix ODE (the Riccati equation), which makes its computation very tractable in practice.

LQ theory can be adapted for tracking of an output signal Y = CX * + with a perturbed linear ODE, see chapter 7 in [START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF]. When we do not know the initial condition, some adaptations are required. Indeed, as the initial condition can have a strong inuence on the optimal control and the optimal cost; it seems much harder to solve the control problem when the initial condition is not known: the current state x(t) is unknown and all the admissible trajectories must be considered. Nevertheless, this problem is solved by the Deterministic Kalman Filter (DKF) by using the fact that the value function S is a quadratic form on the state.

We show in the next section that the Deterministic Kalman Filtering (DKF) is well adapted for developing parameter estimation, as it enables to prole on x 0 , considered as a nuisance parameter. In a two-step approach, it is critical as we need to control the inuence of the nonparametric estimate of Y on the convergence rate. As we use Y (0) as a proxy for Cx * 0 , we need to show that the rate of the two-step estimator is not polluted by the use of nonparametric estimates of the boundary conditions, and that we keep a parametric rate for θ * and x * 0 . This property was carefully checked in [5, 6, 25]; in that paper, as we do not use implicitly or explicitly the derivative of the nonparametric estimate, the mechanics of the proof are dierent.

In the next section, we give some details on LQ theory and on the criterion S.

The classical costs in optimal control consist of an integral term plus a penalty term on the nal state, such as X θ,x0,u (T ) QX θ,x0,u (T ). A preliminary timereversing transformation is used for introducing properly the initial state in the cost C, rather than the nal state. In a second step, we derive the criterion S, and we give a tractable expression for estimation. Finally, we discuss the importance of identiability and observability in the denition on our criterion.

The Deterministic Kalman Filter and the proled cost

Following the Tracking estimator, we look for a candidate X θ,x0,u that minimizes at the same time the discrepancy with the data and the size of the perturbations u L 2 . We consider nearly the same cost as in [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF] C Ŷ ;

x 0 , u, θ, λ = ˆT 0 Ŷ (t) -CX θ,x0,u (t) 2 2 dt + λ ˆT 0 u(t) 2 2 dt (8)
for given λ > 0. We can also add a positive quadratic form x T 0 Qx 0 , where Q is a positive symmetric matrix Q. This additional term permits to introduce easily some prior knowledge on x 0 such that we have a cost dened as

C Ŷ ; u, x 0 , θ, λ = x T 0 Qx 0 + C Ŷ ; x 0 , u, θ, λ .
Moreover, the matrix Q avoids some technical problems in the denition of our criterion S.

For each θ in Θ, we denote

S Ŷ ; θ, λ = inf x0 inf u∈L 2 C Ŷ ; x 0 , u, θ, λ (9)
obtained by proling on the function u and then in the initial condition x 0 .

The function S Ŷ ; x 0 , θ, λ = inf u∈L 2 C Ŷ ; x 0 , u, θ, λ is the criterion used in the case of xed and known initial conditions. Our approach is rather "natural" as we simply prole the regularized criterion x T 0 Qx 0 + S Ŷ ; x 0 , θ, λ . The denition of S mimics the minimization of J n (β|θ, λ) except that GS uses a discretized solution, dened on a B-splines basis. Nevertheless, our estimator possesses two other essential dierences with Generalized Smoothing. As it was already mentioned in [START_REF] Nicolas | A tracking approach to parameter estimation in linear ordinary dierential equations[END_REF], we dene our estimator as the global minimum of the proled cost:

θ K = arg min θ∈Θ S Ŷ ; θ, λ (10)
whereas the GS estimator minimizes a dierent criterion

n i=1 Y i -C X(t i , θ) 2 .
This means that in our methodology, we try to nd a parameter θ that maintain a reasonable trade-o between model and data, whereas the Generalized Smoothing Estimator θGS is dedicated to t the data with the proxy X(•, θ),

without considering the size of model error represented by ūθ . Another impor- tant dierence is in the way we deal with the unobserved part of the system.

For simplicity, let us consider that we observe only the rst k < p components of X, such that the state vector can be written X = X obs , X unobs . For Generalized Smoothing, both functions X obs and X unobs are decomposed in a B-spline basis, and the corresponding coecients β obs and β unobs are obtained by mini-mizing J n β obs , β unobs |θ, λ . Because β unobs does not have to make a trade-o between the data and the ODE model, the estimated missing part Xunobs (•, θ)

is the exact solution to the ODE [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF]. At the contrary, even in the case of partial observations, the perturbed solution X θ,x0,u is used for estimating the missing states and a perturbation exists for each component. Consequently, the estimated hidden states are not solution of the initial ODE. We think that this an advantage for state and parameter estimation with respect to Generalized Smoothing (and NLS) because it avoids to rely too strongly on a uncertain model during estimation. This uncertainty can be caused by errors in parameter estimation, or it can be due model misspecication, such as the presence of a forcing function u * . In our experiments, we show that imposing model uncertainty for the unobserved variables is benecial for error prediction.

Before going deeper into the interpretation and analysis of our estimator, we need to show that the criterion S Ŷ ; θ, λ is properly dened and that we can obtain a tractable expression for computations and for the theoretical analysis of [START_REF] Campbell | Maximum prole likelihood estimation of dierential equation parameters through model based smoothing state estimates[END_REF]. We use the Deterministic Kalman Filter (DKF) to obtain a closed-form expression for the minimal cost w.r.t the control u and x 0 (9).

The initial aim of the DKF is to propose an estimation of the nal state X * (T ) by making a balance between the information brought by the noisy signal Y and the ODE model (see [START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF] for an introduction). We recall the two steps necessary for the lter construction, more details are given in appendix:

1. For a given initial condition x 0 , we nd the minimum cost thanks the fundamental theorem in LQ Theory (presented in A.1), 2. We minimize the quadratic form w.r.t the nal condition.

We give now the main theorem of that section about the existence of the criterion dened in equation [START_REF] Calderhead | Accelerating bayesian inference over nonlinear dierential equations with gaussian processes[END_REF].

Theorem and Denition of S (ζ; θ, λ). Let t → ζ(t) be a function belonging to L ∞ ([0, T ] , R d ) and X θ,x0,u be the solution to the controlled ODE [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF].

For any θ in Θ, λ > 0, Q > 0, there exists a unique optimal control ūθ,λ and initial condition x 0 that minimizes the cost function

C (ζ; u, x 0 , θ, λ) = x T 0 Qx 0 + ˆT 0 ζ(t) -CX θ,x0,u (t) 2 2 + λ u(t) 2 2 dt (11)
The optimal control ūθ,λ is

ūθ,λ (t) = 1 λ E θ (t)X θ, x0,ū θ,λ (t) + h θ (t, ζ) (12)
where E θ and h θ are solutions of the Initial Value Problem

   Ėθ (t) = C T C -A T θ E θ -E θ A θ -1 λ E 2 θ , ḣθ (t, ζ) = -α θ (t)h θ (t, ζ) -β θ (t, ζ) (13)
with

(E θ (0), h θ (0, ζ)) = (Q, 0).
The functions α θ and β θ are dened by

α θ (t) = A θ (t) T + E θ (t) λ β θ (t, ζ) = C T ζ + E θ r θ For all t ∈ [0, T ],
the matrix E θ (t) is symmetric, and the ODE dening the matrix-valued function t → E θ (t) is called the Matrix Riccati Dierential Equation of the ODE [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF].

Finally, the Proled Cost S has the closed form:

S (ζ; θ, λ) = ´T 0 ζ(t) 2 -2r θ (t) T h θ (t, ζ) -1 λ h θ (t, ζ) 2 dt -h θ (T, ζ) T E θ (T ) -1 h θ (T, ζ). ( 14
)
and the nal state is estimated by

X θ, x0,ū θ,λ (T ) = -E θ (T ) -1 h θ (T, ζ). (15)
Remark 2.1. The functions t → (E(t), h(t)) are classically called the adjoint model. They depend also on θ, λ and ζ because of their denition via equation [START_REF] Ellner | Dynamic Models in Biology[END_REF]. Nevertheless, we do not write it systematically for notational brevity. As mentioned in the theorem, it is possible to compute X θ, x0,u θ,λ in a closed-loop form as we can solve in a preliminary stage the adjoint model ( 13) that gives the function E and h for all t ∈ [0, T ]. Thanks to equation ( 12), the closed-form expression of the optimal control ūθ,λ can be plugged into (5). We can compute X θ, x0,u θ,λ by solving the following Final Value Problem:

   ẋ(t) = A θ (t) + E(t) λ x(t) + r θ (t) + h(t, Y ) λ x(T ) = -E(T ) -1 h(T, Y ). ( 16
)
The estimate of x 0 = X θ, x0,u θ,λ (0) of the initial condition is simply the initial value of the Backward ODE [START_REF] Ghasemi | Bayesian parameter estimation for nonlinear modelling of biological pathways[END_REF]. Then by using X θ, x0,u θ,λ , we can compute eectively the control ūθ thanks to (12).

The existence of the criterion S and the fundamental expression ( 14) heavily relies on the nonsingularity of the nal value of the Riccati solution E θ (T ). In particular, the nal state is estimated by

X θ, x0,ū θ,λ (T ) = -E θ (T ) -1 h θ (T, ζ),
and it is then critical to identify the assumptions that could prevent E θ to be singular. Our "Theorem and Denition" is legitimate (and proved in the appendix), because the assumption Q > 0 ensures that E θ (t) is nonsingular for all t in [0, T ]. Moreover, the matrix Q can be thought as a kind of prior for helping the state inference. In our basic denition of the cost (2.2), we put a prior on the norm of the initial condition and our regularization penalizes "huge" solutions. Nevertheless, we can have a more rened prior and use a preliminary guess µ ∈ R d . The modication of the criterion is straightforward by setting

C µ (ζ; u, x 0 , θ, λ) = (x 0 -µ) T Q(x 0 -µ) + ζ(t) -CX θ,x0,u (t) 2 L 2 + λ u(t) 2 L 2 .
By re-parameterizing the initial conditions with y 0 = x 0 -µ and exploiting the relation X θ,x0-µ,u (t) = X θ,x0-µ,u (t) -Φ θ (t, 0)µ (consequence of the linearity of the ODE) , we get that

inf x0 inf u C µ (ζ; x 0 , u, θ, λ) = S (ζ -CΦ θ (•, 0)µ; θ, λ) .
At the opposite, it might be inappropriate in some circumstances to impose such kind of information for the initial condition. This can be the case if the number of observations tends to innity and Ŷ becomes quite close to the truth.

Another situation is when the initial conditions of the unobserved part are largely unknown. Hence, we extend our estimator to the case Q = 0, that corresponds also to our framework for studying the asymptotics of θK . In order to derive relevant and tractable conditions for ensuring the existence of S, we need to ensure that only one trajectory, with a unique initial condition (or nal condition), is the global minimum of C(ζ; u, x 0 , θ, λ). The nonsingularity of E θ (t) is in fact related to the concept of observability in control theory. In the next proposition, we will pave the way to the assumptions on C and the vector eld A θ that can guarantee the general existence of our method. Proposition 2.2. For a given parameter θ ∈ Θ and observation matrix C, the properties 1 and 2 are equivalent:

1. The system outputs Y (t) = CΦ(t, 0)x 0 satisfy ˆT 0 CΦ θ (t, 0)x 0,1 -CΦ θ (t, 0)x 0,2 2 dt = 0 =⇒ x 0,1 = x 0,2 (17)
2. The (nal) observability matrix O θ (T ) is nonsingular

O θ (T ) = ˆT 0 Φ θ (t, 0) C CΦ θ (t, 0)dt (18)
If one of the properties is satised, then E θ (T ) is nonsingular and S is dened for Q ≥ 0.

An important feature of that proposition is that the criterion does not depend on r θ . Moreover, if C is full rank, the matrix E θ (T ) is always nonsingular for all θ in Θ. The criterion 1 means that for a given θ, any solution X θ,x1,0 and X θ,x2,0 of (1) can be distinguished by their partial observation

Y i θ (t) := CX θ,xi,0 (t), i = 1, 2.
The matrix C "gives" enough information about the system so that the observed part is sucient to uniquely characterize the whole system's state.

The next section is dedicated to the derivation of the regularity properties of S. Thanks to the dierent possible expressions for the criterion S, we can show the smoothness in ζ and θ, and compute directly the needed derivatives.

3 Consistency of the Deterministic Kalman Filter Estimator

Properties of the criterion S( Y ; θ, λ)

We have a tractable expression of the cost function S( Y ; θ, λ) for a given θ, but we still need to derive the properties of θ → S( Y ; θ, λ) and θ → S(Y * ; θ, λ) on Θ, and shows some convergence properties. First of all, we need to ensure the

existence of S( Y ; θ, λ); this is the case if the non-parametric estimator Y belongs to L ∞ ([0, T ] , R d ) (more explanations are given in appendix A). We show that for all Y in L ∞ ([0, T ] , R d ), the function θ → S(Y ; θ, λ) is well dened and C 1
on Θ, under some regularity and identiability assumptions, detailed below:

C1: Θ is a compact subset of R p and θ * is in the interior Θ, C2a: Q = 0 and for all θ in Θ, O θ (T ) is nonsingular, C2b: The model is identiable at (θ * , x * 0 ) i.e ∀ (θ, x 0 ) ∈ Θ × X ; CX θ,x0 = CX θ * ,x * 0 =⇒ (θ, x 0 ) = (θ * , x * 0 ), C3: ∀ (t, θ) ∈ [0 , T ] × Θ, (t, θ) → A θ (t) and (t, θ) → r θ (t) are continuous, C4: ∀ (t, θ) ∈ [0, T ] × Θ, (t, θ) -→ ∂A θ ∂θ and (t, θ) -→ ∂r θ
∂θ are continuous.

Condition 2 is about identiability condition: condition 2a is needed for the existence of the criterion S, and is related to the identiability of the initial condition. But C2a is not sucient, and we need Condition 2b for structural identiability, based on the joint identiability at (θ * , x * 0 ). We require that the observed output CX θ * ,x * 0 can be generated on by the couple (θ * , x * 0 ). The identiability problem of systems can be dicult (more than observability). For linear system, several approaches can be used, such as Laplace Transform [START_REF] Bellman | On structural identiability[END_REF], or

Power Expansions [START_REF] Pohjanpalo | System identiability based on the power series expansion of the solution[END_REF], see [START_REF] Miao | On identiability of nonlinear ode models and applications in viral dynamics[END_REF] for a review. So far, most of existing methods are poorly used because they rely on (heavy) formal computations, which limit their interest to low dimensional system. Nonetheless, progress in automatic formal computation has promoted new methods based on dierential algebra and the Ritt's algorithm, that improves identiability checking, [START_REF] Hubert | Essential components of an algebraic dierential equation[END_REF][START_REF] Enqvist | Dierence algebra and system identication[END_REF][START_REF] Noiret | Some eective approaches to check the identiability of uncontrolled nonlinear systems[END_REF].

According to the context, the norm

2 will denote the Euclidean norm in R d , X 2 = d i=1 X 2 i or the Frobenius matrix norm A 2 = i,j |a i,j | 2 . We use the functional norm in L 2 [0, T ] , R d dened by: f L 2 = ´T 0 f (t) 2 2 dt.
Continuity and dierentiability have to be understood according to these norms.

Proposition 3.1. Under conditions 1, 2a and 3 we have: Proposition 3.2. Under conditions 1, 2, 3

A = sup θ∈Θ A θ L 2 < +∞ r = sup θ∈Θ r θ L 2 < +∞ X = sup θ∈Θ X θ L 2 < +∞ Ē = sup θ∈Θ E θ L 2 < +∞ and ∀Y ∈ L ∞ ([0, T ] , R d ), hY = sup θ∈Θ h θ (., Y ) L 2 < +∞ Hence, for all Y in L ∞ ([0, T ] , R d ), the map θ -→ S(Y ; θ, λ) is well dened on Θ (i.e sup θ∈Θ S(Y ; θ, λ) < +∞) We have shown that for all Y in L ∞ [0 T ] , R d , the maps θ → S(Y ; θ, λ)
∀Y ∈ L ∞ ([0, T ] , R d ), θ -→ S(Y ; θ, λ)
is continuous on Θ. Under conditions 1, 2a, 3, 4 it is C 1 on Θ.

In proposition 3.1 we have shown that our criteria θ -→ S(Y ; θ, λ) is well dened ( i.e 0 ≤ S(Y ; θ, λ) < +∞) and here we have demonstrated (using regularity assumptions on the model) that our nite and asymptotic criteria are continuous or even C 1 on Θ. Theses regularity properties justify the use of classical optimization method to retrieve the minimum of S( Y ; ., λ).

Consistency

We show the consistency of the parameter estimator θK when the model is well- specied. As already mentioned, we have dened an M -estimator, and we can proove the consistency (see [START_REF] Van Der | Asymptotic Statistics. Cambridge Series in Statistical and Probabilities Mathematics[END_REF]), by showing

Asymptotics of θ K

The aim of this section is to derive the rate of convergence and asymptotic law of θ K . For this reason, we need more precise assumptions on Ŷ . The way we proceed is based on the plug-in properties of nonparametric estimates, when the functionals of interest are relatively smooth. In the case of series expansion, these properties are well understood [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF][START_REF] Bickel | Nonparametric estimators which can be plugged-in[END_REF]. We focus here on regression splines, as they are well-used in practice and relatively simple to study, although more rened nonparametric estimators can be used in the same context, such as Penalized Splines. We assume that Ŷ has a B-Spline expansion

Y (t) = K k=1 β kK p kK (t) = β T K p K (t)
where β K is computed by linear least-squares, and the dimension K increases with n. We introduce then additional regularity conditions on the ODE model, and on the distribution of observations:

C5: ∀ (t, θ) ∈ [0, T ] × Θ, (t, θ) -→ ∂ 2 A θ ∂θ T ∂θ and (t, θ) -→ ∂r 2
θ ∂θ T ∂θ are continuous, C6:

∂ 2 S(Y * ;θ * ,λ) ∂θ T ∂Y is nonsingular, C7: The observations (t i , Y i ) are i.i.d with V ar(Y i | t i ) = σI d with σ < +∞, C8:
The observation times t i are uniformly distributed on [0 , T ],

C9: It exists s ≥ 1 such that t -→ A θ * (t), t -→ r θ * (t) are C s-1 [0 , T ])R d and √ nK -s -→ 0, K 4 n -→ 0 C10: The meshsize max i |τ i+1,K+1 -τ i,K | -→ 0 whenK -→ ∞
The proofs of the rate and asymptotic normality are somewhat technical, and they are relegated in the Supplementary Materials. We obtain a parametric convergence rate, and the asymptotic normality, by using two facts:

1. θ K -θ * behaves like the dierence Γ( Y ) -Γ(Y * ),
X θ, x0,ū (T ) = -E θ (T ) -1 h θ (T, Ŷ ).
Ẋ(t) = A θ X(t) + r θ (t) + L θ,λ (t) CX(t) -Ŷ (t)
. The matrix L θ,λ is the continuous counterpart of the classical Kalman Gain Matrix, derived from the Filtering Riccati Dierential Equation, see page 313 in [START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF]. In that recursive form, the Deterministic Kalman Filter is somehow similar to the Kalman-Bucy Filter, which is the continuous version of the usual Kalman Filter. Nevertheless, there is a huge dierence in the assumptions because the Kalman-Bucy Filter assumes that X(t) is a Stochastic Dierential Equation, driven by a Brownian Motion W (t). This means that the deterministic perturbation u(t) is replaced by a random pertubation σdW (t). The state estimate is then dierent from the one we consider as it can be shown that the lter is the solution of a stochastic dierential equation driven by the stochastic process CX(t) -Ŷ (t) , see for instance [START_REF] Bensoussan | Stochastic Control of Partially Observable Systems[END_REF]. The state estimate t → X θ, x0,ū (t) is solution of the pertubed ODE, with the control ū computed from all the data Ŷ (s), s ∈ [0, T ] : hence, our state estimation is based on Kalman Smoothing and not on Filtering, as we have a backward integration step. In the rest of that section, we show that the estimator X θ, x0,ū (t) is also a consistent estimator of the state X * (t). In order to do that, we show rst that X(T ) = -E θ (T ) -1 h θ (T, Y ) is a consistent estimator of the nal state.

Final state estimation

In a way, the consistency of the nal state estimator is a rather obvious conclusion. The Deterministic Kalman Filter is initially designed for getting the best possible estimate of the nal state, starting from any initial condition x 0 . It is then normal that we have a good estimator of X * (T ) when Ŷ is close to Y * and θK is close to θ * .

Proposition 5.1. We assume that conditions C1-C4 are satised and that Y is a consistent estimator of Y * . Then, the nal state estimator

X(T ) = -E θ (T ) -1 h θ (T, Y ) converges in probability to X * (T ).
Proof. We show rst that the true nal state value is reached for Y = Y * and

θ = θ * i.e X * (T ) = -E θ * (T ) -1 h θ * (T, Y * ). We recall that S (Y ; θ, λ) = inf x0∈R d inf u∈L 2 C(Y ; x 0 , u, θ, λ) ,
and that S (Y * ; θ * , λ) = 0. The identiability condition 2b implies that the reconstructed state is the exact one. In our case, the minimum is reached when the optimal control u is equal to 0, i.e.

ūθ * ,λ (t) = 1 λ (E θ * (t)X * (t) + h θ * (t, Y * )) which implies that X * (T ) = -E θ * (T ) -1 h θ * (T, Y * ) (E θ * (T ) is nonsingular). We can decompose the dierence X(T ) -X * (T ): X(T ) -X * (T ) = E θ * (T ) -1 h θ * (T, Y * ) -E θ (T ) -1 h θ (T, Y ) = E θ * (T ) -1 h θ * (T, Y * ) -h θ (T, Y ) + E θ * (T ) -1 -E θ (T ) -1 h θ (T, Y )
The convergence will come from the consistency of h θ (T, Y ) and E θ (T ) -1 :

X(T ) -X * (T ) ≤ √ d E θ * (T ) -1 2 h θ * (T, Y * ) -h θ (T, Y ) 2 + √ d h θ (T, Y ) 2 E θ * (T ) -1 -E θ (T ) -1 2
The two right-hand side terms can be controlled easily by the Ŷ -Y * and θ -θ * , as it is shown in Lemma B.2 and B.3. We end up with the following inequalities:

h θ (t, Y ) -h θ (t, Y ) 2 ≤ K 6 e L1 E λ λ Y -Y L 2 + K 7 + K 8 E λ K 4 + K5 λ E λ e L1 E λ λ e 2L1 E λ λ θ -θ + K 9 e 2L1 E λ λ + K 10 e L1 E λ λ E λ θ -θ E -1 θ (T ) -E -1 θ (T ) 2 ≤ K12 λ + K 11 e K13+ K 14 λ θ -θ E -1 θ (T ) ≤ K15 λ h θ (T, Y ) 2 ≤ √ T d 2 C 2 e √ d A+ E λ λ T Y L 2 + √ dE λ r θ Under our conditions, we have θ, Y -→ (θ * , Y * ), which implies that X(T )
converges also in probability.

By plug-in principle, we can also derive the asymptotic normality and the rate of X(T ) as described in the next proposition.

Proposition 5.2. Under conditions C1-C10, the nal state estimator X(T ) is asymptotically normal and

X(T ) -X * (T ) = O P (n -1/2 )
Proof. We have the following decomposition:

X(T ) -X * (T ) = E θ * (T ) -1 h θ * (T, Y * ) -E θ (T ) -1 h θ (T, Y ) = E θ * (T ) -1 h θ * (T, Y * ) -h θ * (T, Y ) + E θ * (T ) -1 h θ * (T, Y ) -h θ (T, Y ) + E θ * (T ) -1 -E θ (T ) -1 h θ (T, Y )
According to Theorem 7 in [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF] Y is a consistent estimator of Y * hence using proposition 5.1 and continuous mapping theorem we have:

X(T ) -X * (T ) = E θ * (T ) -1 h θ * (T, Y * ) -h θ * (T, Y ) + o p (1)
Using the linear representation for h θ * we obtain:

E θ * (T ) -1 h θ * (T, Y * ) -h θ * (T, Y ) = -E θ * (T ) -1 ˆT 0 R θ * (T, s)C T Y (s) -Y * (s) ds
We dene

H(t, θ).Y = E θ (T ) -1 R θ (T, t)C T Y (t)
the linear form such that

-E θ * (T ) -1 h θ * (T, Y * ) -h θ * (T, Y ) = ˆT 0 H(s, θ * ).Y * -H(s, θ * ). Y ds
As for the normality of θK , we can use theorem 9 in [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF] in order the obtain the asymptotic normality of ´T 0 H(s, θ * ).Y * -H(s, θ * ). Y ds with √ n-rate.

Estimation of the states on [0, T ] and inuence of λ

We can estimate the trajectory X * (t) with the smoothed trajectory t → X θ,x0,ū (t)

or with the exact model t → X θ,x0,0 , without the perturbation ū. We need then to have a better understanding of the quality of these two estimates, and in particular of the relevancy of x0 , dened as the unknown initial condition of the Final Value Problem [START_REF] Ghasemi | Bayesian parameter estimation for nonlinear modelling of biological pathways[END_REF]. We have proled the initial condition in the denition of S, in order to separate the estimation of θ * from the estimation of the initial condition. Nevertheless, the estimation of the states is a by-product of the parameter estimation, and the remaining point in our analysis is to ensure that x0 is really a good estimator for x * 0 . This is the case, and we will show more generally that X θ,x0,ū is a good estimator of the trajectory X * . Quite remarkably, the consistency of X θ,x0,ū is the rst result that relies on a assumption on the hyperparameter λ. This is due to the fact that ū is a perturbation computed for tracking Ŷ , while taking into account the model uncertainty estimated by θ instead of θ * . The convergence of Ŷ to Y * and the identiability conditions 2a and 2b (plus regularity conditions) are sucient to ensure the convergence of θ to θ * , without particular assumptions on λ. This is possible because the true model X θ,x0,0 is included into the perturbed model X θ,x0,u .

If λ is not big enough, the size of the perturbation ū 2 L 2 is not highly constrained in the cost function S, and we can have overtting: the estimator CX θ,x0,ū can be quite close to Ŷ with a big ū that makes X θ,x0,ū far from of X * . This problem can be even more important, if we have errors on θ, because ū will have to compensate the errors in the parameter estimation. In that case, we cannot guarantee to have a consistent estimate for x 0 , if we don't have λ -→ ∞. Indeed, the trajectory X θ,x0,ū is the solution to the pertubed initial value problem

   ẋ(t) = A θ (t) + 1 λ E θ,λ (t) x(t) + r θ (t) + 1 λ h(t, Y ) x(T ) = -E θ,λ (T ) -1 h(T, Y ) (19)
Because of the convergence of θ, Y -→ (θ * , Y * ), we can ensure the convergence to the right trajectory if we control λ.

Proposition 5.3. Under conditions C1-C10 and if

λ n -→ ∞, then X θ,x0,ū (t) -→ X * (t) for all t ∈ [0, T ]. Moreover, x 0 -x * 0 = O P (n -1/2 ).
Proof. We rst need to show that for all θ ∈ Θ, and ζ, the functions E θ,λ and h θ,λ (•, ζ) are bounded (as they converge) when λ -→ ∞. As E θ,λ is solution of the matrix equation Ė = C C -A E -EA -1 λ E 2 that depends smoothly in λ -1 ; hence E θ,λ -→ E θ,∞ dened as the solution of the linear matrix ODE Ė =

C C -A E -EA (with E(0) = 0). Moreover, h θ,λ (•, ζ) is solution of the linear ODE ḣ = -α θ,λ h -β θ,λ (•, ζ) with α θ,λ -→ A θ and β θ,λ (•, ζ) -→ β θ,∞ = C ζ + E θ,∞ r θ . As the dependency in λ -1 is smooth, the solution h θ,λ (•, ζ) converges to h θ,∞ , solution of ḣ = -α θ,∞ h -β θ,∞ (•, ζ). Additionally, the dependency in (λ, ζ) is smooth on R + × L 2 and h θ,λ (•, ζ) converges to h θ,∞ (•, Y * ) as (λ, ζ) tends to (∞, Y * ). This means that if (θ, Y ) converges to (θ * , Y * ) as λ -→ ∞, then X θ,x0,λ converges to the solution of the nal value problem    ẋ(t) = A θ * (t)x(t) + r θ * (t) x(T ) = X θ * ,x * 0 ,0 (T ) (20)
as λ -1 E θ,λ and λ -1 h θ,λ (t, Y ) tends to zero when λ -→ ∞, and if x(T ) -→ X * (T ). Because of the uniqueness of the solutions to Initial or Final Value Problem, we have X θ,x0,λ -→ X * for all t ∈ [0, T ].

In proposition 5.1, under conditions C1-C8, we have shown that ( θK n , Ŷn ) converges in probability to (θ * , Y * ) for all λ on Θ × L 2 . By the continuous mapping theorem applied to ( θK n , Ŷn , λ n ), with λ n -→ ∞, we have X θ,x0,ū λn (t) -→ X * (t) for all t ∈ [0, T ] in probability. In particular, we obtain the convergence of x0 to x * 0 .

The asymptotic normality and root-n rate of x0 comes from the asymptotic normality and rates of θ and X(T ). If ψ θ,λ (t, 0) is the resolvant of the ODE ẋ(t) = A θ (t) + 1 λ E θ,λ (t) x(t), we have a closed form for the smoother

X θ,x0,ū λn (t) = ψ θ,λ (t, 0)x 0 + ψ θ,λ ˆt 0 ψ θ,λ (s, 0) -1 r θ (s) + 1 λ h θ,λn (s, Ŷ ) ds
When we evaluate at t = T , we obtain the following formula for the initial

state x0 = ψ θ,λ (T, 0) -1 X(T ) -´T 0 ψ θ,λ (s, 0) -1 r θ (s) + 1 λn h θ,λn (s, Ŷ ) ds.
Hence, x0 is a smooth transformation of θ, X(T ) , and we can conclude by the parametric delta-method.

Choice of λ and cross-validation

Our theoretical analysis shows that when n tends to innity, we have a family of good estimates ( θK λn , x 0,λn ), with λ n -→ ∞. The remaining question is to dene an appropriate selection procedure for λ, that could be used in practice with a nite number of observations (y 1 , . . . , y n ). A straightforward way of selecting λ is to use a cross-validation selection procedure. Indeed, our criterion S Ŷ ; θ, λ is based on a balance between data delity and model delity, and a rough analysis shows that when λ -→ 0, we can select any u in order to interpolate Y and θ has almost no inuence on S Ŷ ; θ, λ . Whereas when λ -→ ∞, the optimal perturbation ū -→ 0, and we get a NLS-like criterion where the observations Y i 's are replaced by the proxy Y .

A good hyperparameter λ n should give a good estimate of the states X * (t) (and of the output Y * (t)), even if we are only interested in parameter estimation.

Anyway, if we want to use the minimization of prediction error for selecting λ and θK λ , we need to have a good estimate of the initial condition x 0 as it is necessary for computing the predictions. We propose then to select λ by minimizing the Sum of Squared Errors

SSE(λ) = n i=1 Y i -CX θ λ , x 0,λ ,0 (t i ) 2 2 . ( 21
)

22

Moreover, this criterion gives a way to reduce the inuence of the nonparametric estimate Ŷ , as we use the original noisy data. This is the selection procedure that we implemented in the experiments part.

Experiments

We use two test beds for evaluating the practical eciency of the deterministic Kalman lter estimator θK ; we compare it with the NLS estimator θNLS and the estimator obtained by Generalized Smoothing θGS . The two models are linear in the states, and they can be linear or nonlinear w.r.t parameters. We use several sample size and several variance error for comparing robustness and eciency. Runge-Kutta algorithm (ode45 in Matlab), and a centered Gaussian noise (with variance σ) is added, in order to obtain the Y i 's. We compare the accuracy of the 3 parameters θK , θ GS and θ N LS , but we are also interested in their mean prediction error dened as

Experimental design

E P X = E (Y1,...,Yn) E θ * ,σ Y * -X (22)
where Y * is a new observation generated with the parameters (θ * , x * 0 , σ), and X is an estimator of the trajectory, based on one of the three estimates θK , θ GS and θ N LS . For the three estimators, the initial condition is estimated consistently: NLS: x 0 is obtained simultaneously with the parameter estimation (as an additional parameter), Kalman: x 0 and λ n are selected as described in section 5.3, Generalized Smoothing: x 0 is the initial value of the estimated curve corresponding to the estimated parameter θ GS , with smoothing parameter λ n selected adaptively as described in [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF].

We insist on the fact that parameter estimation and prediction are two dierent statistical tasks, that are evaluated by dierent criteria. Parameter estimation is required when the parameter has an interest by itself or when the model has an explicative purpose, whereas the prediction error is dedicated to estimation of the state X, in the most ecient way. Our primary interest is parameter estimation but we also discuss prediction for the three methods; as we have seen in section 5, parameter estimation and state estimation are tightly related in particular for the selection of λ. We will consider two possible estimators for the state: a parametric estimator X θ, x0 and a smoothed (or corrected) estimator X θ, x0,ū . For the NLS estimator, the parametric and corrected state estimator are the same, whereas X θGS , x0 and X(•, θGS ) are dierent, as X θ, x0 diers from X θ, x0,ū .

The two test beds are partially observed models with one missing state variable. We compare the ability of the dierent methods to accurately reconstruct the hidden state. Thus, we compute for each estimator the L 2 -distance between the true missing state and the obtained reconstruction after parameter estimation:

∆ Xunobs = E X unobs θ * ,x * 0 -Xunobs L 2 (23)
The nonparametric estimate Ŷ is a regression spline, with a B-spline basis dened on a uniform knot sequence ξ k , k = 1, . . . , K. For each run and each state variables, the number of knots is selected by minimizing the GCV criterion, [START_REF] Ruppert | Semiparametric regression. Cambridge series on statistical and probabilistic mathematics[END_REF]. For optimizing the criterion S, we use the Matlab function 'fminunc' that implements a trust region algorithm for which gradient expression is required.

The computation of the gradient of S w.r.t the parameter θ is computationally involved and is based on the sensitivity equations of the ODE model. The computational details are left in appendix C.

Toy Examples: Partially Observed ODE in 3 D

We consider the autonomous ODE

     ẋ1 = -(k 1 + k 2 )x 1 ẋ2 = k 1 x 1 ẋ3 = k 2 x 1 (24)
where we observe only the variables x 2 and x 3 . Using the notation introduced in this paper, we have

A θ =    -(k 1 + k 2 ) 0 0 k 1 0 0 k 2 0 0    and C = 0 1 0 0 0 1 .
With that model, we show that the conditions introduced in the statistical analysis are workable on some simple models, in particular the conditions for identiability C2a and C2b that needs to be checked. In the case of autonomous system (i.e when A θ and r θ do not depend on time), a simple sucient and necessary criteria is the so-called Kalman criterion: Proposition 6.1. In the case of an autonomous model, the matrix E θ (T ) is nonsingular if and only if the matrix

K A,C =       C CA θ . . . CA d-1 θ       ( 25
)
has a rank equals to d.

The matrix K A,C is usually called the Kalman matrix. In order to dene properly our criterion S, we need to check that condition 2b is also satised (joint identiability of θ and x 0 ). For this model, the analysis is relatively easy and we can use the characterisation proposed by [START_REF] Pohjanpalo | System identiability based on the power series expansion of the solution[END_REF] based on the power series expansion. As, the Kalman matrix ( 25) is

   C CA θ CA 2 θ    =            0 1 0 0 0 1 k 1 0 0 k 2 0 0 -k 1 (k 1 + k 2 ) 0 0 -k 2 (k 1 + k 2 ) 0 0           
, the Kalman condition is fulllled (i.e the matrix rank is 3) if k 1 = 0 or k 2 = 0.

Hence, C2a holds for all relevant cases ( k 1 = 0 or k 2 = 0 correspond to the case where x 2 and x 3 variations are disconnected from x 1 which makes the model useless for explanation or prediction purposes).

For condition C2b, we use the result shown by Pohjanpalo et

h (j) (x 0 , θ, x) = a j (x 0 ) j = 0, 1, . . . ( 26
)
has a unique solution θ. Pohjanpalo et al. showed that for linear autonomous system, this condition is sucient and necessary. In our case, the equation ( 26) can be written as CA j θ X 0 = a j j = 0, 1, . . .

Since the initial condition X 0 = (X 0,1 , X 0,2 , X 0,3 ) is unknown, we have to consider the extended parameter θ = (k 1 k 2 , X T 0 ). The equations for j = 0 CX 0 = a j j = 0, 1, . . . allow us to identify X 0,2 = a 0,2 , X 0,3 = a 0,3 . For j = 1, we have

k 1 X 0,1 = a 1,2 k 2 X 0,1 = a 1,3
and the solutions are X 0,1 = a1,3 k2 and k 1 = a1,2 a1,3 k 2 . Finally, we have a unique solution for k 2 , if we consider the additional equation ( 26) for j = 2. In that case, the system

-k 1 (k 1 + k 2 )X 0,1 = a 2,2 -k 2 (k 1 + k 2 )X 0,1 = a 2,3 has a unique solution k 2 = - a2,3 a1,2+a1,3 .
Well-specied model (Toy Model 1) We test two sample sizes n = 200 and n = 100 (observations times are uniformely sampled between t = 0 and t = 100) and two noise levels σ = 3 and σ = 6. For the computation of the regression splines Ŷ , we select manually the knots location instead of using the GCV driven selection (to avoid overtting). We have placed four equispaced knots respectively at time t = 0, 33, 66 and 100. The true parameter is θ * = (k * 1 , k * 2 ) = (0.0593, 0.0296) and the initial condition x * 0 equals (0, 0, 100). For the Kalman estimator, we select λ n by cross-validation among the values λ v = 10 k k∈[ 5 16] .

(n, σ) Misspecied model (Toy Model 2) In our simulation, we give also some insight in the case of misspecied models. Indeed, our perturbed ODE framework permits to consider naturally the problem of model misspecication, when the

MSE (10 -6 ) ARE (10 -2 ) E P X θ, x0 ∆ X θ, x0 (200, 3) θ N
true model is ẋ(t) = A θ (t)x(t)+r θ (t)+v(t), with v ∈ L 2 ([0, T ] , R d ) an unknown function.
We do not provide any theoretical analysis for this kind of model misspecication. The Kalman estimator gives more accurate estimation than the NLS estimator in that case, as we consider pertubations of the initial model.

Moreover, the optimal control u obtained along the parameter estimation can be used as a correction term to add to the initial model to counter-balance misspecication. This implies potentially a better prediction power. The true model is nearly the same model as above

Ẋ = A θ X + v(t) (27)
with θ * = (k * 1 , k * 2 ) = (0.0593, 0.0296) and x * 0 = (0, 0, 100), but we add a pertubation v : [0, T ] -→ R 3 with entries equal to 0.4 × sin( t 5 ). Nevertheless for parameter estimation, we still use the unperturbed model Ẋ = A θ X.

In the case of the Kalman estimator θK , the optimal control u can be used for correcting the model and for dening a new model

Ẋ = A θ X + u. ( 28
)
We are then interested in evaluating the prediction error of X θ, x0,u , dened as E P X θ, x0,u . We also estimate the error between the true rst state value and the obtained reconstruction with the corrected model. As shown in the introduction, Generalized Smoothing can also evaluate a correction term for θ GS , dened as u(t) = ˙ X(t, θ GS ) -A θ X(t, θ GS ) (where X(t, θ GS ) is the spline corresponding to the estimated parameter θ GS with adaptive λ). In the case of NLS, we cannot compute a correction ū, as the estimated are exactly solution of the ODE for θ N LS . In the case of Generalized Smoothing, we have ∆ X * 1 ; X θ, x0 ≈ ∆ X * 1 ; X θ, x0,ū because the hidden parts are (almost) exactly trajectories of the ODE with parameter θGS . The estimates that change is the Kalman-based one.

(n, σ) The GS parameter estimator is outperformed by the Kalman and the NLS estimator. Our approach improves the estimation accuracy for θ (lower MSE and ARE in every cases)on the NLS estimator. This dierence is bigger than in the well specied case (Toy Model 1), as we are more robust to the presence of a perturbation than the NLS. The Kalman estimator gives also better prediction error in every cases and the correction u slightly improves the prediction errors.

MSE (10 -5 ) ARE (10 -2 ) E P X θ, x0 E P X θ, x0,ū ∆ X unobs θ, x0 ∆ X unobs θ, x0,ū (200, 3) θ N
Nevertheless, the NLS estimator provides the smallest ∆ X * 1 ; X θ, x0 among all estimation methods in every cases but the rst one. Nonetheless using u minimizes in most of case the error for X 1 estimation for our approach and allows us to obtain slightly better result than the NLS estimator.

The correction term ū is related (correlated) to the perturbation t → v(t) as we can in gure 1, where we plot the mean of each component of u, when (n, σ) = (200, 3). Even though the scale is not the same (we need to rescale by 10 -5 for easing comparisons), the correction ū exhibits some important features of the true one, such as oscillations with a period close to the period of v. The analysis of ū is beyond the scope of that paper, but the presence of strong patterns in u can be used to detect misspecication, in the same way that the analysis of residuals permits to detect lack of t in regression models.

Real case example: Methanation reaction

We consider an ODE model introduced in [START_REF] Happel | Multiple isotope tracing of methanation over nickel catalyst[END_REF] for describing the dynamics of carbon monoxide and hydrogen methanation over a supported nickel catalyst by transient isotopic tracer in a gradientless circulating reactor. This Methanation reaction model is a linear autonomous equation in R 4 , with a forcing term. A important dierence w.r.t the previous is the nonlinearity in parameters as we have βC C0 /W +C COl , 0, 0, 0 . The state X is dened as X = X CO , X H 2 O , X CO 2 , X Os , and represents the quantity of the chemical species involved in the reaction. A constant inlet CO ow rate with constant and known fraction of isotope 18 O is introduced within the reactor; the fraction of 18 O present in oxygen atoms for each component is measured at dierent timeframe using a mass spectrometer.

A θ =       - V +V +F C0 0 /W βC C0 /W +C COl 0 0 0 V +V βC H 2 0 /W -V +V +v5 βC H 2 0 /W 0 v5 βC H 2 0 /W V βC CO 2 /W 0 -V +v6 βC CO 2 /W v6 βC CO 2 /W 0 v5 C Os v6 C Os -v5+v6 C Os      
In the model, X j (t) represents the measured fraction of 18 O present in oxygen atoms of the chemical species j at time t. The total amount of oxygen X Os cannot be measured. Some of the parameter are already known:

• F CO i /F CO 0 : inlet/outlet ow rates of CO (0.59/0.45)

• z CO i : the constant fraction of 18 O present in oxygen atoms of the CO inlet ow rate (0.132)

• V /V : rates of production (0.124/0.01)

• C j : concentrations of gas phases in the reaction system (j = CO, H 2 0, CO 2 )

• W total weight of catalyst within system (0.744)

• β volume of dead space (206.1) Our aim is to estimate the parameter θ = C COl , C Os , v 5 , v 6 .

For simulating the datasets, we use two sample sizes n = 100 and n = 50 (observations are uniformely sampled the time interval [0, 40]), with 2 noise levels σ = 0.002 and σ = 0.004. The true parameter value is the estimate provided in [START_REF] Happel | Multiple isotope tracing of methanation over nickel catalyst[END_REF], i.e θ * = (0.1, 11.1, 0.35, 0.008) and with initial condition equals to x * 0 = (0, 0, 0, 0). For the computation of the Kalman estimator, we select λ among 1, 5, 20, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000. Finally, the nonparametric estimate Ŷ is a regression splines, with knots selected manually (instead of GCV selection, because of overtting): we use three equispaced knots at times t = 0, 20, 40.

(n, σ) The dierence can be partly explained by the nonlinearity in parameters that makes their estimation more dicult. We can have estimates that are far from the true parameter value, but in the case of the Kalman estimator, the important errors for the parameter are balanced by a more important correction term ū that ameliorates signicantly state estimation and prediction.

MSE ARE E P X θ, x0 E P X θ, x0,ū ∆ X unobs θ, x0

Discussions

We have considered the statistical problem of parameter and state estimation of a linear Ordinary Dierential Equations as an Optimal Control problem. By doing this, we follow the lines drawn in [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF] or in the two-step approaches, that consist in dening a statistical criterion more adapted to ordinary dierential equation than the likelihood. A new theory was needed in order to assess the statistical eciency of this new estimator, that heavily relies on the Linear-Quadratic Theory. Indeed, the linear structure of the model gives a closed-form for the criterion S which permits to establish the needed regularity properties for statistical analysis. An important question is to determine the conditions under which we can apply the same methodology for nonlinear ODEs. It is probably more involved but the characterization used here is directly generalized by the Pontryagin's Maximum Principle, that gives also a tractable way to solve the optimal control problem.

An important feature of our approach is that we can cope with model misspecication, and the estimation process gives a way to evaluate the lack-of-t thanks to the analysis of the control ū. Thanks to that, we are able to estimate properly the parameters, but also to do prediction and state estimation. Our experiments show that we can have better performance than the classical NLS and Generalized Smoothing and that it is benecial to account for possible perturbation. A good choice for the trade-o hyperparameter λ is then necessary, and our selection methodology is satisfying in practice but needs more insight to explain its inuence for the selection of good predictors, in particular for hidden states. The penalty term u L 2 is an energy related to the degrees of freedom of the predictor X θ,x0,u , but it is not related to the usual criterion of model complexity for smoothing.

In our analysis, we assume that the observability O θ (T ) is nonsingular, which avoids the use of the quadratic form x 0 Qx 0 in the criterion. If this regularization term is used, the mechanics of the proof would be the same, but with Q = Q n that should tend to 0 , as n tends to innity. Nevertheless, it can have consequences on the asymptotics of the estimators, as it corresponds to cases where the loss of information is too big and needs additional information. Quite interestingly, our criterion about identiability remains tractable, and can be relatively easy to check in practice.

It is the main point of the following theorem for a given θ and Z(T ) it ensures the existence, the uniqueness of this control ū and gives a closed form for both ū and C( Ŷ ; Z(T ), ū, θ, λ).

Theorem A.1. Let A ∈ L 2 ([0, T ] , R d×d ) and B ∈ L 2 ([0, T ] , R d×d ) We consider z u the solution of the following ODE:

żu (t) = A(t)z u (t) + B(t)u(t), z(t 0 ) = z 0
and the cost:

C(t 0 , u, U ) = z u (T ) T Qz u (T ) + ˆT t0 z u (t) T W (t)z u (t) + u(t) T U (t)u(t)dt with Q positive,W ∈ L ∞ ([0, T ] , R d×d )
positive matrix for all t ∈ [0, T ] and U (t) denite positive matrix for all t ∈ [0, T ] respecting the coercivity condition:

∃α > 0 s.t ∀u ∈ L 2 ([0, T ] , R d ) : ˆT 0 u(t) T U (t)u(t)dt ≥ α ˆT 0 u(t) 2 2 dt
For a given t 0 we want to minimize the cost

C(t 0 , u, U ) on L 2 ([0, T ] , R d ).
We know it exists an unique control ū, called optimal control, associated to the trajectory z ū, called optimal trajectory, minimizing this cost. Moreover ū is under the closed-feedback loop form u(t) = U -1 (t)E(t)B(t)z u (t) where E is the matricial solution of the ODE:

Ė(t) = W (t) -A(t) t E(t) -E(t)A(t) -E(t)B(t)U (t) -1 B(t) T E(t) E(T ) = -Q
this ODE its called Ricatti equation associated to LQ problem composed of the cost C(t 0 , u, U ) and the ODE żu (t) = A(t)z u (t)+B(t)u(t), z(t 0 ) = z 0 . Moreover E(t) is symetric and the minimal cost is equal to:

C(t 0 , u, U ) = -z T 0 E(t 0 )z 0 . By identifying in the last theorem A with A θ (t) r θ (t) 0 0 , Q with Q 1 , W
with W 1 and U with λI d we obtain the corresponding minimal cost reached for the optimal cost u for a given initial condition x 0

C( Ŷ ; Z θ,x0,u (T ), u, θ, λ) = -X θ,x0,u (0) T E θ (0) X θ,x0,u (0)=2 X θ,x0,u (0) T h θ (0)= α θ (0). (32)
with the associated ODE:

˙ E θ (t) = C T C -A θ E θ -E θ A θ -1 λ E θ 2 , E θ (T ) = -Q ˙ h θ (t) = -C T Y -A θ T h θ -E θ r θ -1 λ E θ h θ , h θ (T ) = 0 ˙ α θ = Y T Y -2 rθ T hθ -1 λ h T θ h θ , α θ (T ) = 0
To be able to apply Theorem A.

1 we need W 1 to belong to L ∞ ([0, T ] , R d×d ), that is why we require Y ∈ L ∞ ([0, T ] , R d ).
A.

-E θ (0) -1 h θ (0) (33)
we denote x 0 the unique initial condition such that X θ, x0,u (T ) = -E θ (0) -1 h θ (0).

In that case the minimal cost is equal to:

C( Ŷ ; Z θ, x0,u (T ), u, θ, λ) = h θ (0) T E θ (0) -1 h θ (0)= α θ (0)
and for a given parameter θ we have:

S Ŷ ; θ, λ = min x0∈R d min u∈L 2 C( Ŷ ; x 0 , u, θ, λ) = h θ (0) T E θ (0) -1 h θ (0)= α θ (0) = h θ (0) T E θ (0) -1 h θ (0) + ´T 0 Y (t) T Y (t) -2 rθ (t) T h θ (t) -1 λ h θ (t) T h θ (t) dt A.3 Minimal cost expression By posing E θ (t) = -E θ (T -t), h θ (t) = -h θ (T -t)
S (Y ; θ, λ) = -h θ (T, Y ) T E θ (T ) -1 h θ (T, Y ) + ´T 0 Y (t) T Y (t) -2r θ (t) T h θ (t, Y ) -1 λ h θ (t, Y ) T h θ (t, Y ) dt (34)
the associated ODE:

Ėθ (t) = C T C -A T θ E θ -E θ A θ -1 λ E 2 θ , ḣθ (t, Y ) = -α θ (t)h θ (t, Y ) -β θ (t, Y ) (E θ (0), h θ (0, Y )) = (Q, 0) (35)
and the functions α and β dened by:

α θ (t) = A θ (t) T + E θ (t) λ β θ (t, Y ) = C T Y + E θ r θ
Hence we have obtained the expression for the optimal control, the minimal cost and the nal state value presented in Theorem 2.2.

B State Estimation: Controls of the variations of the adjoint variables

Lemma B.1. We have

E θ (t) -E θ (t) 2 ≤ K 1 E λ e L1 E λ λ θ -θ by denoting E λ = sup t,θ∈[0, T ]×Θ E θ (t) 2 and E λ ≤ K 2 e L 1 λ
Proof. Thanks to condition 3 ∀θ ∈ Θ t -→ A θ (t) is continuous on [0 , T ] and ∀θ ∈ Θ t -→ E θ (t) is dened on [0 , T ] and obviously continuous on the same interval as an ODE solution.

∀(θ, θ ) ∈ Θ 2 we have:

Ėθ (t) -Ėθ (t) = A θ (t) T E θ (t) -A θ (t) T E θ (t) + E θ (t)A θ (t) -E θ (t)A θ (t) + 1 λ E 2 θ (t) -E 2 θ (t)
and by integrating between 0 and t, taking the norm gives us:

E θ (t) -E θ (t) 2 ≤ ´t 0 A θ (s) T E θ (s) -A θ (s) T E θ (s) 2 ds + ´t 0 E θ (s)A θ (s) -E θ (s)A θ (s) 2 ds + 1 λ ´t 0 E 2 θ (s) -E 2 θ (s) 2 ds
and:

´t 0 E 2 θ (s) -E 2 θ (s) 2 ds ≤ √ d ´t 0 E θ (s) 2 E θ (s) -E θ (s) 2 ds + √ d ´t 0 E θ (s) 2 E θ (s) -E θ (s) 2 ds ≤ 2 √ dE λ ´t 0 E θ (s) -E θ (s) 2 ds by denoting E λ = sup t,θ∈[0, T ]×Θ E θ (t) 2 .
Now we bound the remaining term: F (Q θ , θ, t) = G(Q θ , θ, t)

´t 0 A θ (s) T E θ (s) -A θ (s) T E θ (s) 2 ds ≤ √ d ´t 0 A θ (s) 2 E θ (s) -E θ (s) 2 ds + √ d ´t 0 E θ (s) 2 A θ (s) -A θ (
h θ (t, Y ) -h θ (t, Y ) 2 ≤ √ d C 2 K 3 e L1 E λ λ Y -Y L 2 + √ d C 2 Y L 2 + d √ T rE λ K 4 + K5 λ E λ e L1 E λ λ e 2L1 E λ λ θ -θ + drK 3 e L1 E λ λ K 1 E λ e L1 E λ λ θ -θ + d √ T r θ K 3 e L1 E λ λ E λ θ -θ 2 ≤ K 6 e L1 E
H(Q θ , θ)
with G and H dened by: G(Q θ , θ, t) := -A θ T + E θ λ h θ -C T Y + E θ r θ H (j-1)d+i (Q θ , θ) := δ i,j -(A T θ,i E j + A T θ,j E θ,i + 1 λ E T θ,i E θ,j )

and A θ,i beeing the i -th column of A θ .

For the next subsections we will drop dependence in θ for A θ , r θ , E θ , h θ ∂ E(T ) -1 ∂θ j = -E(T ) -1 ∂E(T ) ∂θ j E(T ) -

( 1 )

 1 . The ODE constraint is relaxed into an inequality constraint dened on the interval [0, T ]. The model constraint is never set to 0 because of the tradeo with the data-tting term n i=1

  is well dened and so are θ → S( Y ; θ, λ) and θ → S(Y * ; θ, λ) as long as the non-parametric estimator Y are well-dened on [0, T ].

For a given

  sample size n and noise level σ, we estimate the Mean Square Error and the mean Absolute Relative Error (ARE) E θ * |θ * -θ| |θ * | by Monte Carlo, based on N M C = 100 runs. For each run, we simulate an ODE solution with a

Fig. 1 :

 1 Fig. 1: Toy Model 2, (n, σ) = (200, 3): Mean correction ū (red curve); rescaled true perturbation 10 -5 × v (blue curve)

  upper bound in the previous inequality we obtain the desired result.

C. 2 ∂r

 2 Gradient computation by sensitivity equationStraightforward computation gives us :∇ θ S( Y ; θ, λ) = -2 ´T 0

  [START_REF] Bellman | On structural identiability[END_REF]

  where Γ is a linear functional,2. if Γ is smooth enough, Γ( Y -Y * ) is asymptotically normal in the case of If conditions C1-C10 are satised, then θ K -θ * = O P (n -1/2 )and θ is asymptotically normal.

	regression splines.
	Conditions C5 and C6 ensures the suciency of second order optimality condi-
	tions for the criteria S. Conditions C7 to C10 are sucient for the consistency
	of Y , as well as for the consistency and the asymptotic normality of the plug-in
	estimators of linear functionals.
	Theorem 4.1. 5 State Estimation

Once the unknown model has been estimated with θK , we focus on the problem of state estimation. From the denition (2.2), the criterion S is built with an estimation of the state based on the solution of the pertubed ODE X θ, x0,ū . The estimate of the initial condition x 0 is derived from a Final Value Problem with the nal state

  al. If the model is ẋ = f (t, x, θ) and the observation function is h(t, θ, x), condition C2b is satised if the nonlinear system

  The results are presented in table3, that gathers the statistics about the parameter estimation accuracy, and the prediction of the complete state, and in particular the estimation of the hidden variable X Os . The Kalman estimator gives more accurate parameter estimates than Nonlinear Least Squares or Generalized Smoothing. The dramatic dierence for the MSE comes from the estimation of C Os that is of greater magnitude than the other parameters, thus ARE seems more relevant for comparisons. However, the MSE enlighten the difculty for NLS and GS estimator to correctly estimate C Os ; moreover, a great number of outliers for C Os estimates have been removed for the NLS estimation before computing ARE and MSE. Additionaly, state estimation improves dramatically, as the prediction error E P and missing state reconstruction ∆ of the Kalman estimator outperforms the two others. This improvement is even more signicative when the correction ū is used.

								∆ X unobs θ, x0,ū
		θ N LS 17.28	1.09	19.43	19.43	6.15	6.15
	(100, 0.002)	θ K	3.60	1.06	8.22	1.11	2.16	0.70
		θ GS	21.54	1.14	19.45	19.45	6.41	6.41
		θ N LS 57.23	2.19	31.49	31.49	12.62	12.62
	(50, 0.002)	θ K	21.38	2.05	9.26	3.49	2.58	1.38
		θ GS	58.05	2.11	44.40	44.39	12.40	12.40
		θ N LS 50.98	1.55	41.60	41.60	12.00	12.00
	(100, 0.004)	θ K	26.76	1.44	7.54	2.11	2.50	1.35
		θ GS	55.61	1.59	33.66	33.66	15.96	15.96
		θ N LS 80.03	2.25	43.13	43.13	14.75	14.75
	(50, 0.004)	θ K	35.87	2.16	28.69	3.06	7.79	1.57
		θ GS	94.30	2.29	44.59	44.59	17.86	17.86
			Tab. 3: Methanation Model

  2 Optimal x 0 selection For a given x 0 we have obtained the minimal cost expression w.r.t control. How can we choose x 0 in order to minimize this minimal cost?We recall that X θ,x0,u (0) = X θ,x0,u (T )

so C( Ŷ ; Z θ,x0,u (T ), u, θ, λ) dened by (

32

) is a quadratic form w.r.t the nal condition ( α θ (0) do not depend on X θ,x0,u (0)). Since it makes no dierence to minimize C( Ŷ ; Z θ,x0,u (T ), u, θ, λ) w.r.t the nal condition instead of x 0 because of unicity of ODE solution we look for the nal condition minimizing

[START_REF] Qi | Asymptotic eciency and nite-sample properties of the generalized proling estimation of parameters in ordinary dierential equations[END_REF]

. Hence if E θ (0) is invertible the minimum is reached for

  we dene our estimator as:

	with the functional criteria:
	θ = arg min θ∈Θ	S Ŷ ; θ, λ

  s) 2 ds ≤ √ dA ´t 0 E θ (s) -E θ (s) 2 + √ dE λ ´t 0 A θ (s) -A θ (s) 2 dUsing these bounds in the main inequality drive us to the following inequal-√ dT E λ A θ -A θ L 2 e 2which gives us the proper results using θ -→ A θ continuity.The bound for E λ is obtained by a direct application of Gronwall's lemma:E θ (t) 2 ≤ |C| 2 + ´t 0 A θ (s) T E θ (s) + A θ (s) T E θ (s) + 1 λ E θ (s) T E θ (s) 2 ds ≤ |C| 2 + 2 √ d ´t 0 A θ (s) 2 + E λ Lemma B.2.Assuming condition C3 and C4 we know it exists constants K i such that:h θ (t, Y ) -h θ (t, Y ) 2 ≤ K 6 e L1 E λ λ Y -Y L 2 + K 7 + K 8 E λ K 4 + K5 λ E λ e L1 E λ Proof. We know that h θ (t, Y ) = -´t 0 R θ (t, s)C T Y (s)ds-´t 0 R θ (t, s)E θ (s)r θ (s)dsWe straightforwardly bound R θ (t, s) 2 by application of Gronwall's lemma:Using successively norm inequalities and Gronwall's lemma we obtain:R θ (t, s) -R θ (t, s) 2 ≤ √ d α θ (t) -α θ (t) 2 R θ (t, s) 2 + √ d R θ (t, s) -R θ (t, s) 2 α θ (t) 2 E θ (t) -E θ (t) 2 R θ (t, s) 2 + √ d R θ (t, s) -R θ (t, s) 2 A + E λ

	and	h θ (t, Y ) 2 ≤	λ E λ θ -θ e 2L1 E λ λ T Y L 2 + λ + K 10 e L1 E λ K 9 e 2L1 E λ λ T d 2 C 2 e + √ √ d A+ E λ λ √ ≤ √ d A θ θ -θ 2 + 1 λ λ ≤ √ d A θ + 1 λ K 1 E λ e L1 E λ λ λ θ -θ K 3 e L1 E λ + √ d A + E λ λ R θ (t, s) -R θ (t, s) 2 ≤ √ de √ dAT K 3 A θ + 1 λ e 2L1 E λ λ λ K 1 E λ e L1 E λ dE λ r θ := K 4 + K5 λ E λ e L1 E λ λ e 2L1 E λ λ θ -θ	θ -θ θ -θ
						√	dE λ	ˆT 0	A √ d(	E λ λ +A)dt
	and we obtain thanks to Cauchy-Schwarz inequality:
									√	d(	E λ λ +A)T
		R θ (t, s) 2 ≤ ≤		√ √	d + de √	√ d A+ E d A + E λ ´t s R θ (u, s) 2 du λ T := K 3 e L1 E λ λ
		≤ |C|	2 e 2	√	d A+	E λ λ	T	λ	E θ (s) 2 ds
	hence					E λ ≤ |C| 2 e 2 √	dAT e 2 √	d	E λ λ

ity:

E θ (t) -E θ (t) 2 ≤ 2 √ d( E λ λ + A) ´t 0 E θ (s) -E θ (s) 2 ds + 2 √ dE λ ´t 0 A θ (s) -A θ (s) 2 ds then Gronwall's lemma gives us E θ (t) -E θ (t) 2 ≤ 2 θ (s) -A θ (s) 2 dt.e ´t 0 2 E θ (t) -E θ (t) 2 ≤ 2 hence ∀ (Y, Y ) ∈ L ∞ ([0, T ] , R d ) we have: h θ (t, Y ) -h θ (t, Y ) 2 ≤ √ d C 2 ´t 0 R θ (t, s) 2 Y (s) -Y (s) 2 ds + √ d C 2 ´t 0 R θ (t, s) -R θ (t, s) 2 Y (s) 2 ds + dr ´t 0 R θ (t, s) -R θ (t, s) 2 E θ (s) 2 ds + dr ´t 0 R θ (t, s) 2 E θ (s) -E θ (s) 2 ds + dr θ θ -θ 2 ´t 0 R θ (t, s) 2 E θ (s) 2 ds

Cauchy Schwarz inequality gives us:

h θ (t, Y ) -h θ (t, Y ) 2 ≤ √ d C 2 R θ (., s) L 2 Y -Y L 2 + √ d C 2 Y L 2 + d √ T rE λ R θ (., s) -R θ (., s) L 2 + dr R θ (., s) L 2 E θ -E θ L 2 + d √ T r θ θ -θ 2 R θ (., s) L 2 E λ

(36)

and applying this bound in 36 gives the following inequality:

  + K 7 + K 8 E λ K 4 + K5 λ E λ e L1 E λLemma B.3. Assuming condition C3 and C4 we know it exists constants K i such that:E -1 θ (T ) -E -1 θ (T ) 2 ≤ K12 λ + K 11 e K13+ K 14Proof. We have already shown∀θ ∈ Θ t -→ E θ (t) is dened on [0 , T ] and obviously continuous on the same interval as an ODE solution.When E -1 θ (t) is dened we know it follows the ODE:

				λ λ	Y -Y L 2
										λ	e 2L1 E λ λ	θ -θ
		+	K 9 e 2L1 E λ λ + K 10 e L1 E λ λ	E λ θ -θ
	By a similar computation we obtain:
	h θ (t, Y ) 2 ≤	√	T d 2 C 2 e	√	d A+	E λ λ	T Y L 2 +	√	dE λ r θ
								λ		θ -θ and E -1 θ (T ) 2 ≤
	K15
	λ

  and we know that R θ (0) = (0, Q r ) so ∂Q θ (0)∂θ = 0, hence we can obtain ∂Q θ (t) ∂Q θ (t) ∂θ ) = ∂F ∂Q (Q θ (t), θ, t) ∂Q θ (t)In order to compute sensitivity equation we need to compute ∂F ∂Q (Q θ , θ, t)and ∂F ∂θ (Q θ , θ, t), for ∂F ∂Q (Q θ , θ, t) and ∂F ∂θ (Q θ , θ, t) we obtain:We also need to compute H(Q θ , θ) partial derivative w.r.t E r and θ, we -th position and E t i is in j -th position.

			∂H(Q θ , θ) ∂θ	(j-1)d+i	= -E t i	∂A j ∂θ	-E t j	1 ∂A i ∂θ
	thus we need to compute because:	∂Q θ (t), ∂θ	solution of the sensitivity equation:
	d dt j E i + A t ( ∂Q θ (t) ∂θ ∂θ A t • ∂ i E j = E t ) = ∂F ∂Q i ∂Aj (Q θ (t), θ, t) ∂θ + E t j ∂Ai ∂Q θ (t) ∂θ ∂θ where matrix	+ ∂Ai ∂F ∂θ ∂θ = ∂Ai ∂θ1 • • • ∂Ai ∂θp (Q θ (t), θ, t)	a d × p
										∂θ
	by solving the Cauchy problem:
			d dt ( ∂θ	+ ∂F ∂θ (Q θ (t), θ, t)
			∂Q θ (0) ∂θ	= 0
	with:
			∂Gi ∂E r (k-1)d+h	(Q θ , θ, t) = -δ i,h r(t) + h λ k
			∂G ∂θ (Q θ , θ, t) = -h T ∂Ai(t) ∂θ	1≤i≤d	-E ∂r(t) ∂θ
	have:
	∂H(E r , θ) ∂E r	(j-1)d+i	= -0 A t j	0 A t i	0 -	1 λ	0 E t j	0 E t i	0
	because:
	• ∂ ∂E r A t j E i + A t i E j = 0 A t j	0 A t i	0 where A t j is in i-th position
	and A t i is in j -th position.
	• 1 λ	∂ ∂E E t j E i = 0 1 λ E t j		0 0 0 + 0 0 1 λ E t i	0 0 where E t j
	is in i and:

∂F ∂Q (Q θ , θ, t) = -A(t) T + E λ ∂Gi ∂E r j (Q θ , θ, t) 0 d 2 ,d ∂H(Q θ ,θ) ∂E r ∂F ∂θ (Q θ , θ, t) = ∂G ∂θ (Q θ , θ, t) ∂H ∂θ (E r , θ)

Appendix: State and Parameter Inference for Partially Observed ODE A Derivation of deterministic Kalman lter estimator using Linear-Quadratic Theory

In this section we describe more precisely how the deterministic Kalman Filter is constructed (see [START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF] for an introduction), it involves two steps:

1. For a given initial condition x 0 we determine the minimum cost expression thanks to theorem A.1 (subsection A.1).

2. inimal cost is a quadratic form w.r.t nal condition and hence it exists an unique nal condition (and hence a unique initial condition by ODE solution uniqueness) minimizing this minimal cost (subsection A.2).

A.1 x 0 xed, minimal cost expression

To derive a closed form for the minimal cost for a given x 0 . For that we dene the reverse time functions:

And by denoting

we can rewrite our cost under the form:

The issue here is to minimize [START_REF] Pronzato | Optimal experimental design and some related control problems[END_REF] in a non-nite dimensional space but thanks to results coming from Optimal control and Riccati theory we know that for a given θ and a given Z(T ) it exists a unique control ū such that C( Ŷ ; Z(T ), ū, θ, λ) = min u∈L 2 C( Ŷ ; Z(T ), u, θ, λ)


Mauro Dragoni
email: dragoni@fbk.euandrea.tettamanzi|celia.pereira@unice.fr

Andrea G B Tettamanzi

Célia Da Costa Pereira

A Fuzzy System for Concept-Level Sentiment Analysis

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction And Related Work

Sentiment Analysis is a kind of text categorization task that aims to classify documents according to their opinion (polarity) on a given subject [START_REF] Pang | Thumbs up? sentiment classification using machine learning techniques[END_REF]. This task has created a considerable interest due to its wide applications. However, in the classic Sentiment Analysis the polarity of each term of the document is computed independently of the domain which the document belongs to. Recently, the idea of adapting term polarities to different domains emerged [START_REF] Blitzer | Biographies, bollywood, boom-boxes and blenders: Domain adaptation for sentiment classification[END_REF]. The rationale behind the idea of such investigation is simple. Let's consider the following example concerning the adjective "small":

1. The sideboard is small and it is not able to contain a lot of stuff. 2. The small dimensions of this decoder allow to move it easily.

In the first text, we considered the Furnishings domain and, within it, the polarity of the adjective "small" is, for sure, "negative" because it highlights an issue of the described item. On the other hand, in the second text, where we considered the Electronics domain, the polarity of such adjective can be considered "positive".

In the literature, different approaches related to the Multi-Domain Sentiment Analysis have been proposed. Briefly, two main categories may be identified: (i) the transfer of learned classifiers across different domains [3] [4], and (ii) the use of propagation of labels through graph structures [5] [6]. Independently of the kind of approach, works using concepts rather than terms for representing different sentiments have been proposed.

Unlike the approaches already discussed in the literature, we address the multidomain sentiment analysis problem by applying the fuzzy logic theory for modeling membership functions representing the relationships between concepts and domains.

Moreover, the proposed system exploits the use of semantic background knowledge for propagating information represented by the learned fuzzy membership functions to each element of the network. To the best of our knowledge, the proposed approach is innovative with respect to the state of the art of Multi-Domain Sentiment Analysis.

The paper is structured as follows. Section 2 introduces the background knowledge and tools used during the development of the system that is described in detail in Section 3. Finally, Section 4 provides a description of how the tasks of the challenge have been addressed and concludes the paper.

Preliminaries

The system is implemented on top of background knowledge used for representing the linguistic connections between "concepts" described in several resources. Below, it is possible to find the list of such resources and the links where further information about them may be found.

WordNet3 [START_REF] Fellbaum | WordNet: An Electronic Lexical Database[END_REF] is one of the most important resources available to researchers in the field of text analysis, computational linguistics, and many related areas. In the implemented system, WordNet has been used as a starting point for the construction of the semantic graph used by the system (see Section 3) However, due to some coverage limitations occurring in WordNet, it has been extended by linking further terms coming from the Roget's Thesaurus [START_REF] Kipfer | Roget's 21st century thesaurus[END_REF].

SenticNet4 [START_REF] Cambria | SenticNet: A publicly available semantic resource for opinion mining[END_REF] is a publicly available resource for opinion mining, which exploits both Artificial Intelligence and Semantic Web techniques to infer the polarity associated with common-sense concepts and represent it in a semantic-aware format. In particular, SenticNet uses dimensionality reduction to calculate the affective valence of a set of Open Mind concepts and represent it in a machine-accessible and machine-processable format.

All resources have been connected by exploiting links contained in ConceptNet5 [START_REF] Liu | ConceptNet-a practical commonsense reasoning tool-kit[END_REF] in order to build a single graph for representing the entire background knowledge exploitable by the system.

System

The main aim of the implemented system is the learning of fuzzy membership functions representing the degree of membership of a concept to a domain in terms of both sentiment polarity and aboutness. The two pillars on which the system has been built are: (i) the use of fuzzy logic for modeling the polarity of a concept with respect to a domain as well as its aboutness, and (ii) the creation of a two-levels graph where the top level represents the semantic relationships between concepts and the bottom level contains the links between all concept membership functions and the domains.

Figure 1 shows the conceptualization of the two-levels graph. Relationships between the concepts of the Level 1 (the Semantic Level) are described by the background knowledge exploited by the system as described in Section 2. The type of relationships are the same generally used in linguistic resource: for example, concepts C 1 and C 3 may be connected through an Is-A relationship rather than the Antonym one. Instead, each connection of the Level 2 (the Sentiment Level) describes the membership of each concept in the different domains considered.

The system has been trained by using the Blitzer dataset6 in two steps: first, the fuzzy membership functions have been initially estimated by analyzing only explicit information present within the dataset (Section 3.1); then, (ii) explicit information has been propagated through the Sentiment Level graph by exploiting the connections defined in the Semantic Level.

Preliminary Learning Phase

The Preliminary Learning (PL) phase aims at estimating the initial polarity of each concept with respect to a domain. The estimation of this value is done by analyzing only explicit information provided by the training set. This phase allows to define the preliminary fuzzy membership functions between the concepts defined in the Semantic Level of the graph and the domains that are defined in the Sentiment Level. Such a value is computed by Equation 1:

polarity * i (C) = k i C T i C ∈ [-1, 1] ∀i = 1, . . . , n, (1)
where C is the concept taken into account, index i refers to domain D i which the concept belongs to, n is the number of domains available in the training set, k i C is the arithmetic sum of the polarities observed for concept C in the training set restricted to domain D i , and T i C is the number of instances of the training set, restricted to domain D i , in which concept C occurs. The shape of the fuzzy membership function generated during this phase is a triangle with the top vertex in the coordinates (x, 1), where x = polarity * i (C), and with the two bottom vertices in the coordinates (-1, 0) and (1, 0), respectively. The rationale is that while we have one point (x) in which we have full confidence, our uncertainty covers the entire space because we do not have any information concerning the remaining polarity values.

Information Propagation Phase

The Information Propagation (IP) phase aims at exploiting the explicit information learned in the PL phase in order to both (i) refine the fuzzy membership function of the known concepts, as well as (ii) to model such functions for concepts that are not specified in the training set, but that are semantically related to the specified ones. Figure 1 presents how the two-levels graph evolves before and after the execution of the IP phase. After the PL phase only four membership functions are modeled: C 1 and C 2 for the domain D 1 , and C 1 and C 5 for the domain D 2 (Figure 1a). However, as we may observe, in the Semantic Level there are concepts that are semantically related to the ones that were explicitly defined in the training set, namely C 3 and C 4 ; furthermore, there are also concepts for which a fuzzy membership function has not been modeled for some domains (i.e. C 2 for the domain D 2 and C 5 for the domain D 1 ).

Such fuzzy membership functions may be inferred by propagating the information modeled in the PL phase. Similarly, existing fuzzy membership functions are refined by the influence of the other ones. Let's consider the polarity between the concept C 3 and the domain D 2 . The fuzzy membership function representing this polarity is strongly influenced by the ones representing the polarities of concepts C 1 and C 5 with respect to the domain D 2 .

The propagation of the learned information through the graph is done iteratively where, in each iteration, the estimated polarity value of the concept x learned during the PL phase is updated based on the learned values of the adjoining concepts. At each iteration, the updated values are saved in order to exploit them for the re-shaping of the fuzzy membership function associating the concept x to the domain i.

The resulting shapes of the inferred fuzzy membership functions will be trapezoids where the extension of the upper base is proportional to the difference between the value learned during the PL phase (V pl ) and the value obtained at the end of the IP phase (V ip ), while the support is proportional to both the number of iterations needed by the concept x to converge to the V ip and the variance with respect to the average of the values computed after each iteration of the IP phase.

Polarity Aggregation And Decision Phases

The fuzzy polarities of different concepts, resulting from the IP phase, are finally aggregated by a fuzzy averaging operator obtained by applying the extension principle (for the technical details see [START_REF] Zadeh | The concept of a linguistic variable and its application to approximate reasoning -i[END_REF]) in order to compute fuzzy polarities for complex entities, like texts, which consist of a number of concepts and thus derive, so to speak, their polarity from them. When a crisp polarity value is needed, it may be computed from a fuzzy polarity by applying one of the defuzzification methods proposed in the literature [START_REF] Zadeh | The concept of a linguistic variable and its application to approximate reasoning -i[END_REF].

Let µ C : [-1, 1] → [0, 1] be the fuzzy interval (i.e., a convex fuzzy set) representing the fuzzy polarity of concept C resulting from the IP phase. Let T be a text (or any other entity that may be regarded as a combination of concepts) related to concepts C 1 , . . . , C n . The fuzzy polarity of T , µ T : [-1, 1] → [0, 1], may be defined as the average of the fuzzy polarities of concepts C 1 , . . . , C n , by applying the extension principle, as follows, for all x ∈ [-1, 1]:

µ T (x) = sup x= 1 n n i=1 xi min i=1,...,n µ Ci (x i ). (2)
The result of the polarity aggregation phase is a fuzzy polarity, whose membership function reflects the uncertainty of the available estimate obtained by the system. In this sense, µ T may be regarded as a possibility distribution of the actual polarity of T . Given x ∈ [-1, 1], the membership degree µ T (x) represent the degree to which it is possible that the polarity of T is x. Here, we are making the assumption that polarity is gradual, i.e., that a text may be more or less negative or positive.

At some point, if a decision must be made based on the polarity of T , some criterion has to be adopted, which takes the uncertainty of the estimate into account. As a matter of fact, a criterion can be defined only with reference to a given application scenario. For instance, if we can afford any desired number of texts and what we want is to pick a few of them whose polarity is certain, we can look for T such that either d T < 0 or a T > 0, i.e., the support of µ T lies entirely on the left or on the right of zero, because in those cases it is certain that polarity is negative (in the former case) or positive (in the latter). In other scenarios, where what we want is to classify each and every text as either negative or positive as accurately as possible, we will have to be less picky and rely on a defuzzification method to transform µ T into a crisp polarity value.

Challenge Tasks and Conclusion

In this paper, we have presented a fuzzy concept-based sentiment analysis system able to model fuzzy membership functions representing the polarities and the aboutness of concepts with respect to a particular domain. The system has been implemented in the context of the ESWC 2014 Concept-Level Sentiment Analysis Challenge. The Tasks proposed by the challenge have been addressed as follows.

Elementary Task: the polarity of each text is computed by aggregating the fuzzy membership functions associated with the extracted concepts. The aggregation operation is performed by applying the extension principle as described in Section 3.3.

Advanced Task #1 and #2: both aspects and concepts (simple and complex) are extracted by exploiting the built knowledge base (as explained in Section 2) and, concerning the Advanced Task #1, its polarity is computed by applying the approach used in the Elementary Task.

Advanced Task #3: similarly to the Elementary Task, the classification of each text is done by analyzing the associations between concepts and domains (independently from the polarity); therefore, the domain of each text is extracted by applying the extension principle of fuzzy sets. Finally, the system has been preliminarily tested on the full version of the Blitzer dataset as shown in Table 1 7 . The system has been compared with three different baselines representing the most well-known machine learning techniques available today demonstrating the feasibility of the proposed approach for addressing the multi-domain sentiment analysis problem.

SVN [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] Naive-Bayes [START_REF] Mccallum | Mallet: A machine learning for language toolkit[END_REF] Max-Entropy [START_REF] Mccallum | Mallet: A machine learning for language toolkit[END_REF]

Fig. 1 :

 1 Fig. 1: The two-layer graph initialized during the Preliminary Learning Phase (a) and its evolution after the execution of the Information Propagation Phase (b).

Table 1 :

 1 Results obtained on the full version of the Blitzer dataset.

	MDFSA Precision MDFSA Recall

https://wordnet.princeton.edu/

http://sentic.net/

http://conceptnet5.media.mit.edu/

http://www.cs.jhu.edu/ mdredze/datasets/sentiment/

Detailed results and tool demo are available at http://dkmtools.fbk.eu/moki/ demo/mdfsa/mdfsa_demo.html