1. No category

Workshop on Data Assimilation
Oxford-Man Institute, September 24–28, 2012
Large sample central limit theorem
for the ensemble Kalman filter (EnKF)
Fran¸cois Le Gland
INRIA Rennes + IRMAR
Val´erie Monbet
Universit´e de Rennes 1 + IRMAR
supported by ANR project PREVASSEMBLE
coordinated by Olivier Talagrand
outline
• EnKF as particle system with mean–field interaction
• identification of the limit
• large sample asymptotics of EnKF
• toy example
• central limit theorem for EnKF
1
(EnKF as particle system with mean–field interactions : 1)
non–linear state–space model of a special form
Xk = fk (Xk−1 ) + Wk
Y k = Hk X k + V k
with
with
Wk ∼ N(0, Qk )
Vk ∼ N(0, Rk )
with non–necessarily Gaussian initial condition X0 ∼ η0
observation noise covariance matrix Rk assumed invertible
in such a model, conditional probability distribution (aka Bayesian filter)
of hidden state Xk given past observations Y0:k = (Y0 , · · · , Yk )
is not Gaussian
ensemble Kalman filter provides a Bayes–free Monte Carlo approach to
numerically evaluate the Bayesian filter, as an alternative to particle filters
2
(EnKF as particle system with mean–field interactions : 2)
I initialization : initial ensemble (X01,f , · · · , X0N,f ) is simulated as i.i.d. random
vectors with probability distribution η0 , i.e. with same statistics as initial condition
X0
1,a
N,a
I prediction (forecast) step : given analysis ensemble (Xk−1
, · · · , Xk−1
) each
member is propagated independently according to (set of decoupled equations)
i,a
) + Wki
Xki,f = fk (Xk−1
with
Wki ∼ N(0, Qk )
notice that i.i.d. random vectors (Wk1 , · · · , WkN ) are simulated here, with same
statistics as additive Gaussian noise Wk in original state equation : in particular
1,a
N,a
(Wk1 , · · · , WkN ) are independent of forecast ensemble (Xk−1
, · · · , Xk−1
)
define empirical covariance matrix
PkN,f
N
1 ∑ i,f
i,f
N,f ∗
=
(Xk − mN,f
)
(X
−
m
k
k
k )
N i=1
of forecast ensemble (Xk1,f , · · · , XkN,f )
with
mN,f
k
N
1 ∑ i,f
=
Xk
N i=1
3
(EnKF as particle system with mean–field interactions : 3)
I correction (analysis) step : given forecast ensemble (Xk1,f , · · · , XkN,f ) each
member is updated independently according to (set of equations with mean–field
interactions)
Xki,a = Xki,f + Kk (PkN,f ) (Yk − Hk Xki,f − Vki )
with
Vki ∼ N(0, Rk )
in terms of Kalman gain mapping defined by
P 7−→ Kk (P ) = P Hk∗ (Hk P Hk∗ + Rk )−1
for any m × m covariance matrix P ,
and in terms of empirical covariance matrix PkN,f of forecast ensemble
notice that i.i.d. random vectors (Vk1 , · · · , VkN ) are simulated here, with same
statistics as additive Gaussian noise Vk in original observation equation : in
particular (Vk1 , · · · , VkN ) are independent of forecast ensemble (Xk1,f , · · · , XkN,f )
4
(EnKF as particle system with mean–field interactions : 4)
mean–field interaction : given that
Xki,a = Xki,f + Kk (PkN,f ) (Yk − Hk Xki,f − Vki )
each analysis member depends on whole forecast ensemble (Xk1,f , · · · , XkN,f )
through empirical probability distribution
µN,f
k
N
1 ∑
δ i,f
=
N i=1 Xk
actually only through empirical covariance matrix PkN,f of forecast ensemble
−→ results in dependent analysis ensemble (Xk1,a , · · · , XkN,a )
5
(EnKF as particle system with mean–field interactions : 5)
question : does ensemble empirical probability distribution converge to Bayesian
filter, defined as
µ−
k (dx) = P[Xk ∈ dx | Y0:k−1 ]
and
µk (dx) = P[Xk ∈ dx | Y0:k ]
i.e. does
µN,f
k
N
1 ∑
δ i,f −→ µ−
=
k
N i=1 Xk
and
µN,a
k
N
1 ∑
δ i,a −→ µk
=
N i=1 Xk
hold in some sense, as N ↑ ∞ ?
answer in general is negative
however, in the linear Gaussian case (i.e. if fk (x) depends linearly on x), then
answer is positive and in particular ensemble empirical mean vector converges to
Kalman filter, i.e.
b N,f
X
k
as N ↑ ∞
N
1 ∑ i,f
b−
=
Xk −→ X
k
N i=1
and
b N,a
X
k
N
1 ∑ i,a
bk
=
Xk −→ X
N i=1
6
(EnKF as particle system with mean–field interactions : 6)
decoupling approach : to study asymptotic behaviour of empirical probability
distributions
µN,f
k
N
1 ∑
=
δ i,f
N i=1 Xk
µN,a
k
and
N
1 ∑
=
δ i,a
N i=1 Xk
of forecast ensemble and analysis ensemble, respectively, approximating i.i.d.
random vectors are introduced
¯ i,f = X i,f , i.e. initial set of i.i.d. random vectors coincides exactly with
initially X
0
0
initial forecast ensemble
these vectors are propagated independently according to (set of fully decoupled
equations)
¯ i,f = fk (X
¯ i,a ) + Wki
X
k
k−1
with
Wki ∼ N(0, Qk )
and
¯ i,a = X
¯ i,f + Kk (P¯ f ) (Yk − Hk X
¯ i,f − Vki )
X
k
k
k
k
with
Vki ∼ N(0, Rk )
¯ i,f
where P¯kf denotes (true) covariance matrix of i.i.d. random vectors X
k
7
(EnKF as particle system with mean–field interactions : 7)
heuristics : these i.i.d. random vectors are close (contiguous) to members in
ensemble Kalman filter, since they
• start from same initial values exactly
• use same i.i.d. random vectors (Wk1 , · · · , WkN ) and (Vk1 , · · · , VkN ) exactly,
already simulated and used in ensemble Kalman filter
esentially a theoretical (not practical) concept
• large sample asymptotics is simple to analyze, because of independance
• true covariance matrix P¯kf is unknown, hence these i.i.d. random vectors are
not computable in practice
in contrast, members in ensemble Kalman filter are computable but dependent,
because they all contribute to / use empirical covariance matrix PkN,f which
results in mean–field interaction
outline
• EnKF as particle system with mean–field interaction
• identification of the limit
• large sample asymptotics of EnKF
• toy example
• central limit theorem for EnKF
8
(identification of the limit : 1)
intuition : limiting probability distributions µ
¯fk and µ
¯ak are probability distributions
¯ i,f and X
¯ i,a respectively, and are completely
of i.i.d. random vectors X
k
k
characterized by integrals of arbitrary test functions
¯ i,f = X i,f and X i,f ∼ η0 , hence µ
I initialization : recall that X
¯f = η0
0
0
0
0
I forecast (expression of µ
¯fk in terms of µ
¯ak−1 ) : recall that
¯ i,f = fk (X
¯ i,a ) + W i
X
k
k
k−1
with
Wki ∼ N(0, Qk )
a
¯ i,a has probability distribution µ
and since X
¯
k−1 (by definition), then
k−1
∫
¯ i,f )] = E[ϕ(fk (X
¯ i,a ) + Wki ))]
ϕ(x′ ) µ
¯fk (dx′ ) = E[ϕ(X
k
k−1
Rm
∫
=
∫
ϕ(fk (x) + w) pW
¯ak−1 (dx)
k (dw) µ
Rm Rm
|
{z
}
Tk ϕ(x)
where pW
k (dw) is Gaussian probability distribution with zero mean vector and
covariance matrix Qk , i.e. probability distribution of random vector Wki
9
(identification of the limit : 2)
I analysis (expression of µ
¯ak in terms of µ
¯fk ) : recall that
¯ i,a = X
¯ i,f + Kk (P¯ f ) (Yk − Hk X
¯ i,f − Vki )
X
k
k
k
k
with
Vki ∼ N(0, Rk )
sufficient conditions on drift function fk can be given, under which µ
¯fk has finite
second order moments, in which case covariance matrix P¯kf is finite
f
¯ i,f has probability distribution µ
¯
and since X
k (by definition), then
k
∫
¯ i,a )] = E[ϕ(X
¯ i,f + Kk (P¯ f ) (Yk − Hk X
¯ i,f − V i ))]
ϕ(x′ ) µ
¯ak (dx′ ) = E[ϕ(X
k
k
k
k
k
Rm
∫
=
∫
ϕ(x + Kk (P¯kf ) (Yk − Hk x − v)) qkV (v) dv µ
¯fk (dx)
Rm Rd
|
{z
}
TkKF (¯
µfk ) ϕ(x)
where qkV (v) is Gaussian density with zero mean vector and invertible covariance
matrix Rk , i.e. probability density of random vector Vki
10
(identification of the limit : 3)
on the other hand, Bayesian filter, defined as
µ−
k (dx) = P[Xk ∈ dx | Y0:k−1 ]
satisfies recurrent relation
∫
∫
′
ϕ(x′ ) µ−
(dx
)=
k
and
µk (dx) = P[Xk ∈ dx | Y0:k ]
∫
ϕ(fk (x) + w) pW
k (dw) µk−1 (dx)
Rm Rm
|
{z
}
Tk ϕ(x)
Rm
and (Bayes rule)
∫
Rm
∫
ϕ(x′ ) µk (dx′ ) =
Rm
∫
′
ϕ(x′ ) qkV (Yk − Hk x′ ) µ−
(dx
)
k
Rm
′
qkV (Yk − Hk x′ ) µ−
k (dx )
with initial condition µ−
0 = η0
clearly, limiting probability distributions of forecast / analysis ensemble do not
a
¯fk ̸= µ−
and
µ
¯
coincide with Bayesian predictor / filter, i.e. µ
k ̸= µk , except in the
k
linear Gaussian case
outline
• EnKF as particle system with mean–field interaction
• identification of the limit
• large sample asymptotics of EnKF
• toy example
• central limit theorem for EnKF
11
(large sample asymptotics of EnKF : 1)
indeed, intuition is correct : ensemble empirical probability distributions
µN,•
k
N
1 ∑
=
δ i,•
N i=1 Xk
do converge (in some sense) as N ↑ ∞ to the probability distribution µ
¯•k of i.i.d.
¯ i,• (hence, not to the Bayesian filter)
random vectors X
k
reference
F. Le Gland, V. Monbet and Vu–Duc Tran Large sample asymptotics for
the ensemble Kalman filter , chapter 22 in The Oxford Handbook of
Nonlinear Filtering, 2011
12
(large sample asymptotics of EnKF : 2)
Theorem (law of large numbers) under mild assumptions on drift function fk
and on test function ϕ
∫
N
∑
1
ϕ(Xki,• ) −→
ϕ(x) µ
¯•k (dx)
N i=1
Rm
in probability as N ↑ ∞
Theorem (Lp –convergence and rate of convergence) under mild assumptions on
drift function fk and on test function ϕ, and provided initial condition X0 has
finite moments of any order p
sup
N ≥1
for any order p
√
∫
N
∑
1
N ( E|
ϕ(Xki,• ) −
ϕ(x) µ
¯•k (dx) |p )1/p < ∞
N i=1
Rm
13
(large sample asymptotics of EnKF : 3)
to summarize : ensemble Kalman filter
• gain matrix depends on empirical covariance matrix
• ensemble empirical probability distribution converges to the wrong limit
(different from Bayesian filter), except for linear Gaussian model
√
• rate of convergence 1/ N
vs. (any brand of) particle filter
• weighted empirical probability distribution of particle system converges to the
correct limit (Bayesian filter)
√
• rate of convergence 1/ N , central limit theorem
question : is there any advantage to use ensemble Kalman filter ?
idea : prove central limit theorem (and compare asymptotic error variances)
outline
• EnKF as particle system with mean–field interaction
• identification of the limit
• large sample asymptotics of EnKF
• toy example
• central limit theorem for EnKF
14
(toy example : 1)
toy example
linear Gaussian system
the target distribution, i.e. the Bayesian filter, is known explicitly as a
Gaussian distribution, with mean and covariance provided by the Kalman
filter
hidden state
√
Xk = a Xk−1 + 1 − a2 Wk
with
Wk ∼ N(0, σ 2 )
initial condition X0 ∼ N(0, σ 2 ) so that stationarity holds
observations
Yk = X k + Vk
Vk ∼ N(0, s2 )
with
numerical values
a
σ
s
0.5
1
1
15
(toy example : 2)
true state and Kalman filter
2.5
true state
Kalman filter
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−2.5
0
10
20
30
40
50
60
70
80
90
100
16
(toy example : 3)
true state, Kalman filter and EnKF empirical mean : 10000 elements
2.5
true state
Kalman filter
EnKF
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−2.5
0
10
20
30
40
50
60
70
80
90
b N,a with N = 10000 members
EnKF empirical mean vector X
k
100
17
(toy example : 4)
EnKF histogram + Kalman filter density : time 100
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−3
−2
−1
0
EnKF histogram with N = 10000 members
1
2
3
4
18
(toy example : 5)
conclusion : not only does the EnKF empirical mean vector
b N,a
X
k
N
1 ∑ i,a
=
X
N i=1 k
bk , but more generally the EnKF empirical
converge to the Kalman filter X
probability distribution
N
∑
1
µN,a
=
δ i,a
k
N i=1 Xk
converges to the Gaussian distribution with moments given by the Kalman filter
next different question : how fast does the empirical mean vector converge to the
Kalman filter, e.g. is the normalized difference
√
N
∑
1
i,a
b N,a − X
bk ) = √
b
N (X
(X
k
k − Xk )
N i=1
asymptotically normally distributed and how to compute the asymptotic variance ?
19
(toy example : 6)
toy example (continued)
numerical simulations : for EnKF / bootstrap particle filter / particle filter with
optimal importance distribution
• M Monte Carlo runs
• each Monte Carlo run evaluates one ensemble / particle average, based on N
members / particles and compares this average with the (known) limit
• histogram of the M normalized differences is shown
same toy example : stationary linear Gaussian system
20
(toy example : 7)
histogram of normalized error of EnKF empirical mean : 10000 Monte Carlo runs
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−3
−2
−1
0
1
√
2
3
b N,a − X
bk ) for N = 1000
histogram of EnKF normalized differences N (X
k
members and M = 10000 Monte Carlo runs
21
(toy example : 8)
histogram of normalized error of PF empirical mean : 10000 Monte Carlo runs
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−3
−2
−1
0
1
2
histogram of (bootstrap) particle filter normalized differences
N = 1000 particles and M = 10000 Monte Carlo runs
3
√
¯ N − mk ) for
N (X
k
22
(toy example : 9)
histogram of normalized error of PF empirical mean : 10000 Monte Carlo runs
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−3
−2
−1
0
1
2
histogram of (optimal) particle filter normalized differences
N = 1000 particles and M = 10000 Monte Carlo runs
3
√
¯ N − mk ) for
N (X
k
23
(toy example : 10)
conclusion (as drawn from these simulations) : in terms of speed of convergence
(a smaller asymptotic variance means a faster convergence)
PF with optimal importance distribution ≫ bootstrap PF ≫ EnKF
outline
• EnKF as particle system with mean–field interaction
• identification of the limit
• large sample asymptotics of EnKF
• toy example
• central limit theorem for EnKF
24
(central limit theorem for EnKF : 1)
Theorem (central limit theorem) under mild assumptions on drift function fk
and on test function ϕ
1
√
N
N
∑
∫
[ϕ(Xki,• ) −
i=1
Rm
ϕ(x) µ
¯•k (dx) ] =⇒ N(0, v¯k• (ϕ))
in distribution as N ↑ ∞, with (more or less explicit) expression for asymptotic
variance v¯k• (ϕ)
beyond the qualitative statement
• recurrence relations for the asymptotic variance ?
• practical computations ?
25
(central limit theorem for EnKF : 2)
because of the recursive nature of the ensemble Kalman filter, it seems natural to
prove the CLT by induction, hence to rely on the following strategy, already used
in K¨
unsch (Annals of Statistics, 2005)
Lemma
if
′
• conditionally w.r.t. FN , the r.v. ZN
converges in distribution to a Gaussian
r.v. with zero mean and variance V ′ , in the sense that for any fixed u
′
E[exp{j u ZN
} | FN ] −→ exp{− 12 u2 V ′ }
in probability, and in L1 by the Lebesgue dominated convergence thorem,
′′
• the r.v. ZN
is measurable w.r.t. FN , and converges in distribution to a
Gaussian r.v. with zero mean and variance V ′′ , i.e. for any fixed u
′′
E[exp{j u ZN
}] −→ exp{− 12 u2 V ′′ }
′
′′
then the r.v. ZN = ZN
+ ZN
converges in distribution to a Gaussian r.v. with
zero mean and variance V = V ′ + V ′′ , as N ↑ ∞
26
(central limit theorem for EnKF : 3)
I initialization : recall that X0i,f ∼ η0 and µ
¯f0 = η0 , hence
1
√
N
N
∑
i=1
∫
[ϕ(X0i,f ) −
Rm
ϕ(x) µ
¯f0 (dx) ] =⇒ N(0, v¯0f (ϕ))
in distribution as N ↑ ∞, with asymptotic variance
∫
∫
v¯0f (ϕ) = var(ϕ, η0 ) =
|ϕ(x)|2 η0 (dx) − |
Rm
Rm
ϕ(x) η0 (dx) |2
27
(central limit theorem for EnKF : 4)
I forecast step : recall that µ
¯fk = µ
¯ak−1 Tk , where
∫
Tk ϕ(x) =
ϕ(fk (x) + w) pW
k (dw)
Rm
Proposition asymptotic variance of forecast approximation
a
v¯kf (ϕ) = v¯k−1
(Tk ϕ) + σk2,f (ϕ)
in terms of
• asymptotic variance of analysis approximation at previous step,
evaluated for a transformed test function
• asymptotic Monte Carlo variance
∫
∫
Tk |ϕ|2 (x) µ
¯ak−1 (dx) −
σk2,f (ϕ) =
Rm
Rm
|Tk ϕ(x)|2 µ
¯ak−1 (dx)
28
(central limit theorem for EnKF : 5)
hint :
ZN
∫
N
1 ∑
[ϕ(Xki,f ) −
ϕ(x′ ) µ
¯fk (dx′ ) ]
= √
N i=1
Rm
1
= √
N
N
∑
∫
i,a
[ϕ(fk (Xk−1
) + Wki ) −
i=1
Rm
ϕ(x′ ) µ
¯fk (dx′ ) ]
N
1 ∑
i,a
i,a
[ ϕ(fk (Xk−1
) + Wki ) − Tk ϕ(Xk−1
)]
= √
N i=1
1
+√
N
′
′′
= ZN
+ ZN
N
∑
i=1
∫
i,a
[ Tk ϕ(Xk−1
)−
Rm
Tk ϕ(x) µ
¯ak−1 (dx) ]
29
(central limit theorem for EnKF : 6)
conditionnally w.r.t. the analysis ensemble at previous stage, the random variables
i,a
i,a
eik = ϕ(fk (Xk−1
) + Wki ) − Tk ϕ(Xk−1
)
are independent (not identically distributed), with mean zero and variance
∫
∫
i,a
i,a
2
sik =
|ϕ(fk (Xk−1
) + w)|2 pW
ϕ(fk (Xk−1
) + w) pW
k (dw) − |
k (dw)|
Rm
Rm
i,a
i,a
= Tk |ϕ|2 (Xk−1
) − |Tk ϕ(Xk−1
)|2
respectively, for i = 1 · · · N , and it follows from our law of large numbers that
N
N
1 ∑ i
1 ∑
i,a
i,a
sk =
[ Tk |ϕ|2 (Xk−1
) − |Tk ϕ(Xk−1
)|2 ]
N i=1
N i=1
∫
−→ σk2,f (ϕ) =
Rm
∫
Tk |ϕ|2 (x) µ
¯ak−1 (dx) −
in probability and in L1 as N ↑ ∞
Rm
|Tk ϕ(x)|2 µ
¯ak−1 (dx)
30
(central limit theorem for EnKF : 7)
therefore
N
∑
u
N,a
N,a
′
} | Fk−1
] = E[exp{j √
eik } | Fk−1
] −→ exp{− 12 u2 σk2,f (ϕ)}
E[exp{j u ZN
N i=1
in L1 as N ↑ ∞
moreover, it follows from the induction assumption that
1
′′
ZN
=√
N
N
∑
i=1
∫
i,a
[ Tk ϕ(Xk−1
)−
Rm
a
Tk ϕ(x) µ
¯ak−1 (dx) ] =⇒ N(0, v¯k−1
(Tk ϕ))
hence
′′
a
E[exp{j u ZN
}] −→ exp{− 12 u2 v¯k−1
(Tk ϕ)}
as N ↑ ∞
. . . this concludes the proof 2
31
(central limit theorem for EnKF : 8)
I analysis step : recall that µ
¯ak = TkKF (¯
µfk ) µ
¯fk , where
∫
TkKF (¯
µfk ) ϕ(x) =
ϕ(x + Kk (P¯kf ) (Yk − Hk x − v)) qkV (v) dv
Rd
introduce new transform
f
f ∗
f
KF f
KF f
,
ϕ)
(x
−
m
¯
)
M
(¯
µ
)
ϕ(x)
+
(x
−
m
¯
)
ϕ(x)
=
T
(¯
µ
QKF
(¯
µ
k
k
k
k)
k
k
k
k
where the matrix MkKF (¯
µfk , ϕ) is defined as
f
f
¯
MkKF (¯
µfk , ϕ) = Hk∗ (Hk P¯kf Hk∗ + Rk )−1 LKF
(¯
µ
,
ϕ)
(I
−
K
(
P
k
k
k
k ) Hk )
f
and where the matrix LKF
(¯
µ
k
k , ϕ) is defined as
∫ ∫
f
f
′
¯
(Y
−
H
x
−
v)
ϕ
(x
+
K
(
P
LKF
(¯
µ
,
ϕ)
=
k
k
k
k
k ) (Yk − Hk x − v))
k
Rm
Rd
qkV (v) dv µ
¯fk (dx)
. . . motivation will come later
32
(central limit theorem for EnKF : 9)
Proposition asymptotic variance of analysis approximation
f
2,a
v¯ka (ϕ) = v¯kf (QKF
(¯
µ
)
ϕ)
+
σ
k
k
k (ϕ)
in terms of
• asymptotic variance of analysis approximation at previous step,
evaluated for a transformed test function
• asymptotic Monte Carlo variance
∫
∫
σk2,a (ϕ) =
TkKF (¯
µfk ) |ϕ|2 (x) µ
¯fk (dx) −
Rm
Rm
|TkKF (¯
µfk ) ϕ(x)|2 µ
¯fk (dx)
33
(central limit theorem for EnKF : 10)
hint :
ZN
∫
N
1 ∑
ϕ(x′ ) µ
¯ak (dx′ ) ]
= √
[ϕ(Xki,a ) −
N i=1
Rm
∫
N
1 ∑
= √
[ϕ(Xki,f + Kk (PkN,f ) (Yk − Hk Xki,f − Vki )) −
ϕ(x′ ) µ
¯ak (dx′ ) ]
N i=1
Rm
N
1 ∑
i,f
)
ϕ(X
[ ϕ(Xki,f + Kk (PkN,f ) (Yk − Hk Xki,f − Vki )) − TkKF (µN,f
= √
k )]
k
N i=1
N
1 ∑ KF N,f
+√
[ Tk (µk ) ϕ(Xki,f ) − TkKF (¯
µfk ) ϕ(Xki,f ) ]
N i=1
1
+√
N
N
∑
i=1
′
′′
′′′
= ZN
+ ZN
+ ZN
∫
µfk ) ϕ(Xki,f ) −
[ TkKF (¯
Rm
TkKF (¯
µfk ) ϕ(x) µ
¯fk (dx) ]
34
(central limit theorem for EnKF : 11)
conditionnally w.r.t. the forecast ensemble at current stage, the random variables
eik = ϕ(Xki,f + Kk (PkN ) (Yk − Hk Xki,f − Vki ))
are independent (not identically distributed), with mean zero and variance
∫
sik =
|ϕ(Xki,f + Kk (PkN ) (Yk − Hk Xki,f − v)) |2 qkV (v) dv
Rm
∫
−|
Rm
ϕ(Xki,f + Kk (PkN ) (Yk − Hk Xki,f − v)) qkV (v) dv |2
i,f
i,f 2
2
KF N,f
= TkKF (µN,f
)
|ϕ|
(X
)
−
|T
(µ
)
ϕ(X
k
k
k
k
k )|
respectively, for i = 1 · · · N , and it follows from our law of large numbers that
N
N
1 ∑ i
1 ∑ KF N,f
i,f 2
sk =
[ Tk (µk ) |ϕ|2 (Xki,f ) − |TkKF (µN,f
)
ϕ(X
k
k )| ]
N i=1
N i=1
∫
∫
−→ σk2,a (ϕ) =
TkKF (¯
µfk ) |ϕ|2 (x) µ
¯fk (dx) −
|TkKF (¯
µfk ) ϕ(x)|2 µ
¯fk (dx)
Rm
in probability and in L1 as N ↑ ∞
Rm
35
(central limit theorem for EnKF : 12)
therefore
N
∑
u
′
} | FkN,f ] = E[exp{j √
eik } | FkN,f ] −→ exp{− 12 u2 σk2,a (ϕ)}
E[exp{j u ZN
N i=1
in L1 as N ↑ ∞
′′′
refrain from jumping directly to the last term ZN
...
′′
. . . indeed, more is still to come from the second term ZN
36
(central limit theorem for EnKF : 13)
contribution of nonlinear transformation : Taylor expansion
′′
ZN
N
1 ∑ KF N,f
= √
[ Tk (µk ) ϕ(Xki,f ) − TkKF (¯
µfk ) ϕ(Xki,f ) ]
N i=1
N ∫
1 ∑
[ ϕ(x + Kk (PkN,f ) (Yk − Hk Xki,f − v))
= √
N i=1 Rd
− ϕ(x + Kk (P¯kf ) (Yk − Hk Xki,f − v)) ] qkV (v) dv
N ∫
1 ∑
= √
ϕ′ (Xki,f + Kk (P¯kf ) (Yk − Hk Xki,f − v))
N i=1 Rd
(I − Kk (P¯kf ) Hk ) (PkN,f − P¯kf ) Hk∗ (Hk P¯kf Hk∗ + Rk )−1
(Yk − Hk Xki,f − v) qkV (v) dv + εN
k
in terms of the derivative of the Kalman gain mapping P 7−→ Kk (P )
37
(central limit theorem for EnKF : 14)
further approximation : law of large numbers
N ∫
∑
1
′′
ZN
=
ϕ′ (Xki,f + Kk (P¯kf ) (Yk − Hk Xki,f − v))
N i=1 Rd
(I −
Kk (P¯kf )
Hk )
√
N (PkN,f − P¯kf ) Hk∗ (Hk P¯kf Hk∗ + Rk )−1
(Yk − Hk Xki,f − v) qkV (v) dv + εN
k
N ∫
∑
√
1
i,f
urow (Xk , v) N (PkN,f − P¯kf ) ucol (Xki,f , v) qkV (v) dv + εN
=
k
N i=1 Rd
N
1 ∑
= √
[ (Xki,f − m
¯ fk )∗ MkKF (¯
µfk , ϕ) (Xki,f − m
¯ fk )
N i=1
∫
(x − m
¯ fk )∗ MkKF (¯
µfk , ϕ) (x − m
¯ fk ) µ
¯fk (dx) ] + εN
−
k
Rm
38
(central limit theorem for EnKF : 15)
hint for the last identity : relies on the following argument
N ∫
√
1 ∑
i,f
urow (Xk , v) N (PkN,f − P¯kf ) ucol (Xki,f , v) qkV (v) dv
N i=1 Rd
= trace[
√
N ∫
∑
1
ucol (Xki,f , v) urow (Xki,f , v) qkV (v) dv ]
− P¯kf )
N i=1 Rd
N (PkN,f
√
≈ trace[ N (PkN,f − P¯kf ) Mk ]
∫
with
Mk =
∫
Rm
Rd
ucol (x, v) urow (x, v) qkV (v) dv µ
¯fk (dx)
and on
trace[
√
N (PkN,f
N
∑
1
− P¯kf ) Mk ] ≈ √
[ (Xki,f − m
¯ fk )∗ Mk (Xki,f − m
¯ fk )
N i=1
∫
−
(x − m
¯ fk )∗ Mk (x − m
¯ fk ) µ
¯fk (dx) ]
Rm
39
(central limit theorem for EnKF : 16)
upon identifying
urow (x, v) = ϕ′ (x + Kk (P¯kf ) (Yk − Hk x − v)) (I − Kk (P¯kf ) Hk )
and
ucol (x, v) = Hk∗ (Hk P¯kf Hk∗ + Rk )−1 (Yk − Hk x − v)
this yields
∫
MkKF (¯
µfk , ϕ) =
∫
Rm
Rd
ucol (x, v) urow (x, v) qkV (v) dv µ
¯fk (dx)
f
f
¯
= Hk∗ (Hk P¯kf Hk∗ + Rk )−1 LKF
(¯
µ
,
ϕ)
(I
−
K
(
P
k
k
k
k ) Hk )
where
LKF
µfk , ϕ) =
k (¯
∫
∫
Rm
Rd
(Yk − Hk x − v) ϕ′ (x + Kk (P¯kf ) (Yk − Hk x − v))
qkV (v) dv µ
¯fk (dx)
40
(central limit theorem for EnKF : 17)
therefore
′′
′′′
+ ZN
ZN
∫
N
1 ∑ KF f
f
f
N
= √
[ Qk (¯
µk ) ϕ(Xki,f ) −
QKF
(¯
µ
)
ϕ(x)
µ
¯
(dx)
]
+
ε
k
k
k
k
N i=1
Rm
where new transform (now better motivated) is defined by
f
f ∗
f
KF f
KF f
QKF
(¯
µ
)
ϕ(x)
=
T
(¯
µ
)
ϕ(x)
+
(x
−
m
¯
)
M
(¯
µ
,
ϕ)
(x
−
m
¯
k
k
k
k
k
k
k
k)
and where εN
k converges to zero in probability as N ↑ ∞
41
(central limit theorem for EnKF : 18)
it follows from the induction assumption that
′′
′′′
ZN
+ ZN
∫
N
1 ∑ KF f
f
f
N
QKF
=√
[ Qk (¯
µk ) ϕ(Xki,f ) −
(¯
µ
)
ϕ(x)
µ
¯
(dx)
]
+
ε
k
k
k
k
N i=1
Rm
f
a
=⇒ N(0, v¯k−1
(QKF
(¯
µ
k
k ) ϕ))
hence
f
′′
′′′
a
E[exp{j u (ZN
+ ZN
)}] −→ exp{− 12 u2 v¯k−1
(QKF
(¯
µ
k
k ) ϕ)}
as N ↑ ∞
recall that
N
∑
u
′
E[exp{j u ZN
} | FkN,f ] = E[exp{j √
eik } | FkN,f ] −→ exp{− 12 u2 σk2,a (ϕ)}
N i=1
in L1 as N ↑ ∞
. . . this concludes the proof 2
42
(central limit theorem for EnKF : 19)
I practical computations : iterating the recurrence relations
f
2,a
v¯ka (ϕ) = v¯kf (QKF
(¯
µ
)
ϕ)
+
σ
k
k
k (ϕ)
and
a
v¯kf (ϕ) = v¯k−1
(Tk ϕ) + σk2,f (ϕ)
yields
f
2,a
v¯ka (ϕ) = v¯kf (QKF
(¯
µ
)
ϕ)
+
σ
k
k
k (ϕ)
2,f
f
2,a
a
KF f
= v¯k−1
(Tk QKF
ϕ)
+
σ
(Q
(¯
µ
)
(¯
µ
)
ϕ)
+
σ
k
k
k
k
k
k (ϕ)
|
{z
}
|
{z
}
σk2 (ϕ)
RkKF
with initialization
f
2,a
2,a
KF
(¯
µ
)
ϕ)
+
σ
(ϕ)
=
var(Q
(η
)
ϕ,
η
)
+
σ
v¯0a (ϕ) = v¯0f (QKF
0
0
0
0
0
0
0 (ϕ)
43
(central limit theorem for EnKF : 20)
a
v¯ka (ϕ) = v¯k−1
(RkKF ϕ) + σk2 (ϕ)
a
a
KF
2
v¯k−1
(RkKF ϕ) = v¯k−2
(Rk−1
RkKF ϕ) + σk−1
(RkKF ϕ)
..
.
KF
a
KF
v¯la (Rl+1
· · · RkKF ϕ) = v¯l−1
(RlKF · · · RkKF ϕ) + σl2 (Rl+1
· · · RkKF ϕ)
..
.
v¯1a (R2KF · · · RkKF ϕ) = v¯0a (R1KF · · · RkKF ϕ) + σ12 (R2KF · · · RkKF ϕ)
hence
v¯ka (ϕ) = v¯0a (R1KF · · · RkKF ϕ) +
k
∑
KF
· · · RkKF ϕ)
σl2 (Rl+1
l=1
in terms of backward–propagated functions
KF
KF
Rl+1:k
ϕ = Rl+1
· · · RkKF ϕ
44
(central limit theorem for EnKF : 21)
further simplifications occur in the special case of
• linear (and quadratic) test functions ϕ
• linear drift function fk
indeed
• forward–propagated distributions µ
¯fk and µ
¯ak are Gaussian distributions with
moments given by the Kalman filter
KF
• backward–propagated functions Rl+1
· · · RkKF ϕ remain quadratic at all steps
and explicit calculations can be obtained
forthcoming event
International Conference on Ensemble Methods in Geophysical Sciences
Toulouse, November 11-16, 2012
http://www.meteo.fr/cic/meetings/2012/ensemble.conference/
organized within the ANR project PREVASSEMBLE