Another look at estimating parameters in systems of ordinary - - PowerPoint PPT Presentation

Another look at estimating parameters in systems of

• rdinary differential equations via regularization

Ivan Vujaˇ ci´ c∗ Seyed Mahdi Mahmoudi∗∗, Ernst Wit∗∗

∗Department of Mathematics, Vrije Universiteit Amsterdam, The Netherlands ∗∗ Department of Statistics and Probability, University of Groningen, The Netherlands

Van Dantzig seminar, March 6, 2014

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 1 / 47

Introduction

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 2 / 47

Motivation

System of ordinary differential equations (ODEs) in the standard form x′(t) = f(x(t),t;θ), t ∈ [0,T], x(0) = ξ, (1) where x(t),ξ ∈ Rd and θ ∈ Rp. x(t;θ,ξ) denotes the solution of (1) for given ξ,θ. Many processes in science and engineering are modelled by (1).

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 3 / 47

Example: The FitzHugh-Nagumo neural spike potential equations

x′

1(t) = c{x1(t)−x1(t)3/3+x2(t)},

x′

2(t) = − 1 c{x1(t)−a+bx2(t)}.

x1 represents the voltage across an axon membrane. x2 summarizes outward currents. Example: ξ1 = −1, ξ2 = 1. a = 0.2, b = 0.2, c = 3.

5 10 15 20 −2 −1 1 2 time x1 5 10 15 20 −1.0 0.0 0.5 1.0 time x2

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 4 / 47

The problem

Noisy observations of x(t;θ0,ξ0) of some states of the system are available: yi(tj) = xi(tj;θ0,ξ0)+εi(tj), i = 1,...,d1;j = 1,...,n. where 0 ≤ t1 ≤ ··· ≤ tn ≤ T. For simplicity, we consider Gaussian errors.

Goal

Estimate θ0 from the data Y, where Y = (yi(ti))ij. This is inverse problem for the coefficients in a system of ODEs. If ξ0 is not known it is considered as parameter and estimated as well.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 5 / 47

FhNdata from R package ’CollocInfer’

5 10 15 20 −2 −1 1 2 time x1

• ● ●
• ●
• ● ●
• ● ●
• ● ●
• ●
• ●
• ●
• 5

10 15 20 −2 −1 1 2 time x2

• ● ●
• ●
• ●
• ● ● ●
• ●
• ● ●
• ● ●
• ● ●
• ●
• ● ● ● ●
• Ivan Vujaˇ

ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 6 / 47

Some existing approaches

1

Non-linear least squares (MLE)

2

Smooth and match estimators

3

Generalized profiling procedure

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 7 / 47

Non-linear least squares

1

Numerical solution x(t;θ,ξ) of the ODE system.

2

Criterion Mn(θ,ξ). Mn(θ,ξ) = −

d1

∑

i=1 n

∑

j=1

logp(yi(tj)| xi(tj;θ,ξ)), where p(yi(tj)| xi(tj;θ,ξ)) is the probability density function of the data. NLS estimator is √n-consistent and asymptotically efficient. Assumption: the maximum step size of the numerical solver goes to zero. Otherwise NLS is not consistent. [Xue et al., 2010]

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 8 / 47

Reference

Xue, H.,Miao, H. and Wu, Hulin (2010). Sieve estimation of constant and time-varying coefficients in nonlinear ordinary differential equation models by considering both numerical error and measurement error. Annals of statistics, 38:2351–2387.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 9 / 47

Smooth and match estimator

1

Smoother x(t)

2

Criterion Mn(θ) Mn(θ) =

T

x′(t)−f( x(t),θ)qw(t)dt. The √n-consistency was shown for: regression splines for 0 < q ≤ ∞. [Brunel et al., 2008] kernel estimator for q = 2. [Gugushvili and Klaassen, 2012] Asymptotic normality was shown for regression splines for q = 2. [Brunel et al., 2008]

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 10 / 47

References

Brunel, N. J. et al. (2008). Parameter estimation of ode’s via nonparametric estimators. Electronic Journal of Statistics, 2:1242–1267. Gugushvili, S. and Klaassen, C. A. J. (2012). √n-consistent parameter estimation for systems of ordinary differential equations: bypassing numerical integration via smoothing. Bernoulli, 18:1061–1098.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 11 / 47

Smooth and match estimator: integral criterion

1

Smoother x(t)

2

Criterion Mn(θ,ξ) Mn(θ,ξ) =

T

x(t)−ξ −

t

0 f(x(t),θ)ds2dt.

For f(x(t),θ) = g(x(t))θ, g : Rd → Rd×p √n-consistency was shown for: local polynomials [Dattner and Klaassen(2013)]. certain step function estimator in [Vujacic et al.(2014)].

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 12 / 47

References

Dattner, I., Klaassen, C.A.: Estimation in systems of ordinary differential equations linear in the parameters. arXiv preprint arXiv:1305.4126, (2013) Vujaˇ ci´ c, I., Dattner, I., Gonz´ alez, J., Wit, E. : Time-course window estimator for ordinary differential equations linear in the parameters. Statistics and Computing, (2014) (To appear in Statistics and Computing. Published

• nline. )

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 13 / 47

Generalized profiling procedure

1

Model based smoother x(t;θ,ξ), where x = argminx∈XmJ(x).

2

Criterion Mn(θ,ξ) Inner criterion J(x) = −

d1

∑

i=1 n

∑

j=1

logp(yi(tj)|xi(tj;θ,ξ))+λ

d

∑

i=1

wi

T

0 {x′ i(t)−fi(x(t),t,θ)}2dt,

Outer criterion Mn(θ,ξ) = −

d1

∑

i=1 n

∑

j=1

logp(yi(tj)| xi(tj;θ,ξ)). The estimator is consistent and asymptotically efficient. [Ramsay et al.(2007)] The only frequentist approach that can handle partially observed systems.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 14 / 47

Reference

Ramsay, J.O., Hooker, G., Campbell, D., Cao, J.: Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5): 741–796, (2007)

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 15 / 47

Summary

The framework: Stochastic or deterministic approximation x of the solution. Criterion function Mn.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 16 / 47

This talk

For simplicity let ξ0 be known. Otherwise, define augmented vector θ ∗ = (θ,ξ). The framework: 1.

• x(θ) = argminx∈XmTα,γ(x|θ),

2.

• θn = argminθ∈ΘMn(θ|

x(θ),Y). We consider log-likelihood criterion Mn.

Aim

Define Tα,γ such that: It yields asymptotically efficient estimator. It can handle partially observed systems.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 17 / 47

Structure of the rest of the presentation

1

Background on regularization theory.

2

Applying the regularization theory to ODE problem.

3

Asymptotic results.

4

Conceptual comparison with the generalized profiling procedure. Only theory in this talk; no simulation studies.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 18 / 47

• 1. Background on regularization

theory.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 19 / 47

References

Vasin, V. V. and Ageev, A. L. (1995). Ill-posed problems with a priori information, volume 3. Walter de Gruyter. Engl, H. W., Hanke, M., and Neubauer, A. (1996). Regularization of inverse problems, volume 375. Springer. P¨

• schl, C. (2008).

Tikhonov regularization with general residual term. University Innsbruck.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 20 / 47

Well-posedness in the sense of Hadamard

Let F : X → Y where X,Y are linear normed spaces and consider the equation F(x) = y, (2) x ∈ X, y ∈ Y. The problem (2) is well-posed in the sense of Hadamard on (X,Y) if:

1

The solution of (2) exists.

2

It is unique.

3

It is continuous with respect to y. The problem (2) is ill-posed on (X,Y) if it is not well-posed.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 21 / 47

Objective functional

Equation F(x) = y, (3) can be solved on a set S ⊂ X by minimizing objective functional J(x) = F(x)−y2,

• n S.

Quasisolution of equation (3) on S ⊂ X is any minimizer of J on S. It is also called pseudo solution or least squares solution. Remark: This idea dates back to the beginning of the 19th century (Gauss, Legendre).

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 22 / 47

Stabilizing functional and Tikhonov regularization

Ω - stabilizing functional Ω incorporates a priori information on the smoothness of the solution x. Ω is usually given by a norm or a semi-norm on X. Tikhonov regularization involves minimization of the Tikhonov functional Tα(x) = J(x)+αΩ(x−x0), where x0 is trial solution α ≥ 0 is regularization parameter

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 23 / 47

Similarity functional and generalized Tikhonov regularization

Similarity functional S incorporates a priori information on values of x. S measures the closeness of the solution to this a priori information. Generalized Tikhonov regularization involves minimization of Tα,γ(x) = J(x)+αΩ(x−x0)+γS(x), where γ ≥ 0 is the penalty parameter. We will call Tα,γ generalized Tikhonov functional. We will call any minimizer of Tα,γ generalized Tikhonov regularizer.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 24 / 47

Finite-dimensional approximation

Numerical minimization - on some finite-dimensional subspace Xm ⊂ X. Minimal assumptions:

1

X1 ⊂ X2 ⊂ ...

2

∪∞

m=1Xm is dense in X.

Remarks: In statistics literature Xms are called sieves. Finite-dimensional approximation is a form of regularization. It is called self regularization or regularization by projection.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 25 / 47

Summary

Generalized Tikhonov functional Tα,γ(x) = J(x)+αΩ(x−x0)+γS(x).

1

Objective functional J.

2

Stabilizing functional Ω.

3

Similarity functional S.

4

Finite-dimensional approximation.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 26 / 47

• 2. Applying the regularization

theory to ODE problem.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 27 / 47

Is the problem x′(t) = f(x(t),t;θ), t ∈ [0,T], x(0) = ξ, ill-posed?

NO.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 28 / 47

Is the problem x′(t) = f(x(t),t;θ), t ∈ [0,T], ill-posed?

YES.

Even if the initial conditions are known, non-uniqueness can still be introduced through finite dimensional approximation.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 29 / 47

Finite-dimensional approximation

The construction is for fixed θ. We suppress dependence on θ for notational simplicity. Solution of the system belongs to (C1[0,T])d. Xm ⊂ C1[0,T] linear subspace of dimension m with basis {h1,...,hm}. Each component of x is approximated by an element of Xm. xi(t) =

m

∑

k=1

βikhk(t) = β ⊤

i h(t),

where βi = (βi1,...,βim)⊤ h(t) = (h1(t),...,hm(t))⊤

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 30 / 47

J - objective functional

Consider x′(t) = f(x(t),t;θ), t ∈ [0,T], for fixed θ. Define F(x(·)) = x′(·)−f(x(·),·,θ), ODE system is equivalent to the equation F(x) = 0d. The corresponding objective functional is J(x) = x′ −f(x,·,θ)2

2,w.

where w = (w1,...,wd), wi > 0 for i = 1,...,d, x2,w =

• ∑d

i=1 wi

T

0 x2 i (t)dt.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 31 / 47

Ω - stabilizing functional

Here we list two options common in the literature. Norm in (L2[0,T])d Ω(x) = x2

2,w = d

∑

i=1

wi

T

0 x2 i (t)dt.

Norm in Sobolev space (H2[0,T])d Ω(x) =

d

∑

i=1

vi

T

0 {x′′ i (t)}2dt.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 32 / 47

S - similarity functional

The observations Y represent: the data for the problem of the estimation of θ0. a priori information for the problem of finding the solution x(t;θ0,ξ0). We have: The true distribution of the data g. Postulated, a priori distribution of the solution p(·|x(·;θ,ξ)). ”Distance” between g and p(·|x(·;θ,ξ)) should be small. Taking KL divergence yields: S(x) = KL(g(·);p(·|x)) ≈ −

d1

∑

i=1 n

∑

j=1

logp(yi(tj)|xi(tj)).

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 33 / 47

Tα,γ - generalized Tikhonov functional

For fixed θ the generalized Tikhonov functional is Tα,γ(x(β)) = J(x(β))+αΩ(x(β)−x0)+γS(x(β)), (4) where the functionals J, Ω and S are defined in previous slides. The regularized solution is found by optimizing (4) over Xd

m.

This can be achieved by optimizing (4) with respect to β over Rdm:

• β = argminβ∈RdmTα,γ(x(β)),

and applying basis expansion xi(t) = ∑m

k=1

βikhk(t) = β ⊤

i h(t).

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 34 / 47

Artificial example:smooth and match estimators fit into the proposed framework

Tα,γ(x) = J(x)+αΩ(x−x0)+γS(x). Take trial solution x0 to be some smoother of the data.

• x = argminx∈Xd

mT∞,0(x) = x0.

Mn(θ) =

T

x′(t)−f( x(t),θ)qw(t)dt, Remark: Similarly, taking trial solution x0 to be numerical solution yields NLS.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 35 / 47

• 3. Asymptotics

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 36 / 47

The estimator

1.

• x(θ) = argminx∈Xd

mTα,γ(x|θ),

2.

• θn = argminθ∈ΘMn(θ|

x(θ),Y). We consider log-likelihood criterion Mn and Ω(x) =

d

∑

i=1

vi

T

0 {x′′ i (t)}2dt.

Result for Ω(x) = x2

2,w

carries over without any modification.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 37 / 47

Reference

Qi, X. and Zhao, H. (2010). Asymptotic efficiency and finite-sample properties of the generalized profiling estimation

• f parameters in ordinary differential equations.

The Annals of Statistics, 38(1):435–481.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 38 / 47

Union of sieves is dense in (C1[0,T])d

An(θ,ξ) = xo(θ,ξ,·)−w∞ ∨

• dxo

dt (θ,ξ,·)− dw dt

• ∞

∨

• d2xo

dt2 (θ,ξ,·)− d2w dt2

• ∞

Bn(θ,ξ) = xu(θ,ξ,·)−v∞ ∨

• dxu

dt (θ,ξ,·)− dv dt

• ∞

∨

• d2xu

dt2 (θ,ξ,·)− d2v dt2

• ∞

.

Lemma

Under Assumption 2 of [Qi and Zhao, 2010], there exist a sequence of finite-dimensional subspaces Xn of C1[0,T] such that for any compact subset Θ0 of Θ and any compact subset Ξ0

• f Ξ, it holds

lim

n→∞rn = 0,

where rn = max

• sup

(θ,ξ)∈Θ0×Ξ0

inf

w∈Xn,w(0)=ξ o

An(θ,ξ), sup

(θ,ξ)∈Θ0×Ξ0

inf

v∈Xn,v(0)=ξ u

Bn(θ,ξ)

• .

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 39 / 47

Consistency and asymptotic efficiency

Theorem (Consistency)

Let Assumptions 1-5 from [Qi and Zhao, 2010] hold. If as n → ∞

1

rn → 0

2

αn → 0

3

γn → 0 then θn −θ0 = oP(1). Tα,γ(x) = J(x)+αΩ(x−x0)+γS(x).

Theorem (Asymptotic efficiency)

Let Assumptions 1-6 from [Qi and Zhao, 2010] hold. If rn = o(n−1), αn = o(n−2) and γn = o(n−2) as n → ∞ then θn is asymptotically normal with the same asymptotic covariance matrix as that of the maximum likelihood estimation.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 40 / 47

• 4. Conceptual comparison with

the generalized profiling procedure.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 41 / 47

Generalized profiling fits into the proposed framework

Inner criterion of the generalized profiling procedure J(x) = −

d1

∑

i=1 n

∑

j=1

logp(yi(tj)|xi(tj;θ))+λ

d

∑

i=1

wi

T

0 {x′ i(t)−fi(x(t),t,θ)}2dt

can be written as J(x) = λ 1 λ S(x)+J(x)

• = λT0,1/λ(x).

Thus, model based smoother x is

• x = argminx∈Xd

mT0,1/λ(x). Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 42 / 47

Smoothing VS Generalized Tikhonov regularization

”For solutions to the dynamic systems, however, the roles of goodness of fit and ’roughness penalty’ seems more likely reversed, with fidelity to the ODE the major concern and the ’error distribution’ of the data an afterthought (Chong Gu - in the discussion section of [Ramsay et al.(2007)]). In the generalized profiling: Fidelity to the ODE term is the penalty. λ must approach ∞: leads to ill conditioning in the optimization. In the regularization formulation Fidelity to the ODE term is the main term— objective functional. γ must approach 0: no ill conditioning in the optimization.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 43 / 47

Generalized Tikhonov regularizer and its special cases

Parameters Tα,γ(x)

• x = argminx∈Xd

mTα,γ(x)

α > 0, γ > 0 J(x)+αΩ(x−x0)+γS(x)

• Gen. Tikhonov’s regularizer

α = 0, γ = 0 J(x) Ivanov’s quasi solution α > 0, γ = 0 J(x)+αΩ(x−x0) Tikhonov’s regularizer α = 0, γ > 0 J(x)+γS(x) model based smoother α = ∞, γ = 0 J(x0)/δ(x−x0) trial solution x0

Table: The last row should be interpreted as Tα,0(x) → J(x0)/δ(x−x0) as α → +∞, where δ is the Dirac’s delta function.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 44 / 47

Conclusion

Regularization provides a coherent and principled framework for defining an approximation of the solution of ODE. ODE system is solved in the least square sense.

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 45 / 47

Acknowledgments

Bartek Knapik Department of mathematics, Vrije Universiteit Amsterdam, The Netherlands Itai Dattner Department of statistics, University of Haifa, Israel

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 46 / 47

Questions, comments,...

Ivan Vujaˇ ci´ c (VU) Generalized Tikhonov regularization for ODEs Van Dantzig seminar, March 6, 2014 47 / 47