EnKF-based particle filters Jana de Wiljes, Sebastian Reich, Wilhelm - PowerPoint PPT Presentation

EnKF-based particle filters Jana de Wiljes, Sebastian Reich, Wilhelm Stannat, Walter Acevedo June 20, 2017

Filtering Problem Signal 1 √ Signal d X t = f ( X t ) d t + 2 C d W t Obs 0.5 Observations 0 dY t = h ( X t ) dt + R 1 / 2 dV t . −0.5 time Goal : determine π ( x | Y 0 : t )

State of the art Linear Model : Kalman- Bucy Filter x M x M t d t + b d t − P M t H T R − 1 ( H ¯ x M d ¯ = A ¯ t d t − d Y t ) t Non-linear Model : approximate with empirical measure, i.e., M � π ( x | Y 0 : t ) ≈ 1 δ ( x − X i t ) . M i = 1 Ansatz : define modified evolution equation for particles X i t

Ensemble Kalman Filter (EnKF) P ( x | Y 0: t ) P ( x ) time t t H T R − 1 � � t − 1 d X i t = f ( X i t ) d t + C d W i 2 P M HX i x M t d t + H ¯ t d t − 2 d Y t M M � � = 1 1 x M X i P M ( X i x M t )( X i x M t ) T ¯ = t − ¯ t − ¯ t t t M M − 1 i = 1 i = 1 Works remarkably well in practice : meteorology, oil reservoir exploration But : theoretical understanding is largely missing

Recent accuracy results for EnKF Interacting particle representation of the model error : d X i t = f ( X i t ) d t + CC ⊤ ( P M t ) − 1 ( X i x M t − ¯ t ) d t t H T R − 1 � � − 1 2 P M HX i x M t d t + H ¯ t d t − 2 d Y t Setting: d Y t = X t d t + R 1 / 2 d W t with R = ε I Results: ([dWRS16]) ◮ Control of spectrum of covariance matrix P M over t ∈ [ 0 , ∞ ) . t t � 2 in expectation ◮ Control of estimation error e t = � X ref x M − ¯ t E [ E t ] = O ( ǫ 1 / 2 ) (1) and pathwise over fixed interval t ∈ [ 0 , T ] . E t ≥ c 1 ǫ q ] ≤ O ( ǫ 1 / 2 − η − q ) P [ sup (2) t ≤ T

Numerical confirmation for L63 time − averaged mean squared error 1 10 error theoretical order 0 10 − 1 10 − 2 10 − 4 − 3 − 2 − 1 10 10 10 10 measurement noise variance �

EnKF-based particle filters Ansatz: modified evolution equation for the particles e.g., of the form � d X i t = f ( X i t ) d t + C d W i X j t ds ji t − t j Aims: ◮ achieve first/second order accuracy ◮ to go beyond Gaussian assumption ◮ consistency for M → ∞ ◮ hybrid formulation to combine different interacting particle filters ([CRR16])

Linear ensemble transform filters (LETF) Given : M samples x f i ∼ π ( x ) (prior) Analysis step: � x a x f j = i d ij (3) i with transformation matrix D = { d ij } subject to M � d ij = 1 (4) i = 1 Examples: ENKF, ESRF, NETF, ETPF (see [RC15])

Ensemble transform Particle Filter (ETPF) Given : ◮ M samples x f i ∼ π ( x ) (prior ensemble) ◮ normalized importance weights w i ∝ π ( y | x f i ) (likelihood) Desired: M samples x a i ∼ π ( x | y ) Ansatz: replace resampling step with linear transformation by interpreting it as discrete Markov chain given by transition matrix D ∈ R M × M s.t. d ij ≥ 0 and � � 1 d ij = 1 , d ij = w i . M i j Benefits : localization, hybrid

Ensemble transform Particle Filter (ETPF) Then the posterior ensemble members are distributed according to the columns of the transformation   d 1 j   d 2 j �   ˜ j = E [˜ X a and x a X a x a   j ∼ j ] = i d ij (5) .  .  .   i d Mj Example. Monomial resampling   w 1 w 1 · · · w 1   w 2 w 2 · · · w 2     D Mono := w ⊗ 1 = . . . . ...   . . . . . .   w M w M · · · w M

Ensemble transform Particle Filter (ETPF) Then the posterior ensemble members are distributed according to the columns of the transformation   d 1 j   d 2 j �   ˜ j = E [˜ X a and x a X a x a   j ∼ j ] = i d ij (5) .  .  .   i d Mj Example. Monomial resampling   w 1 w 1 · · · w 1   w 2 w 2 · · · w 2     D Mono := w ⊗ 1 = . . . . ...   . . . . . .   w M w M · · · w M Here : X f and X a are independent. Idea : increase correlation between X f and X a

Ensemble transform Particle Filter (ETPF) Solve optimization problem M � t ij || x f i − z f j || 2 D ETPF = arg min M i , j = 1 to find transformation matrix D ETPF that increases correlation ([Rei13]). Remarks : ◮ consistent for M → ∞ ◮ first-order accurate for finite M , i.e, M M � � x a = 1 x a w i x f i = (6) i M i = 1 i = 1 ◮ But : not second-order accurate

Second-order accuracy The analysis covariance matrix M � P a = 1 � ( x a i − x a )( x a i − x a ) T M i = 1 is equal to the covariance matrix defined by the importance weights, i.e. M � P a = w i ( t k )( x f i − x a )( x f i − x a ) T . i = 1

First-order accurate LETFs Notation : X f = ( x f M ) , X a = ( x a 1 , . . . , x f 1 , . . . , x a M ) ∈ R N x × M and analogously w = ( w 1 , . . . , w M ) T ∈ R M × 1 and W = diag ( w ) LETF is first-order accurate if 1 M X a 1 = X f w . This holds if D satisfies 1 M D 1 = w .

First-order accurate LETFs Class of first-order accurate LETFs D 1 = { D ∈ R M × M | D T 1 = 1 , D 1 = M w } Examples : ◮ D EnKF / ∈ D 1 ◮ D ESRF / ∈ D 1 ◮ D ETPF ∈ D 1 ◮ D Mono ∈ D 1

Second-order accurate LETF Analysis covariance matrix: P a = 1 � M X f ( D − w1 T )( D − w1 T ) T ( X f ) T (7) Importance sampling covariance matrix: P a = X f ( W − ww T )( X f ) T . (8) Thus the following equation has to hold for second-order accuracy: ( D − w1 T )( D − w1 T ) T = W − ww T (9) Class of second-order accurate LETFs D 2 = { D ∈ D 1 | ( D − w1 T )( D − w1 T ) T = W − ww T } . (10)

Second-order corrected LETF Given : D ∈ D 1 Goal : correct transformation � D ∈ D 2 Ansatz : � D = D + ∆ with ∆ ∈ R M × M such that ∆ 1 = 0 , ∆ T 1 = 0 ([dWAR17]). Need to solve algebraic Riccati equation: M ( W − ww T ) − ( D − w1 T )( D − w1 T ) T =( D − w1 T )∆ T + ∆( D − w1 T ) T + ∆∆ T .

Numerical example I Gaussian prior, non-Gaussian likelihood: Figure: Prior and posterior distribution for the single Bayesian inference step

Numerical example I Figure: Absolute errors in the first four moments of the posterior distribution.

Numerical example II Lorenz-63 model, first component observed infrequently ( ∆ t = 0 . 12) and with large measurement noise ( R = 8): Figure: RMSEs for various second-order accurate LETFs compared to the ETPF, the ESRF, and the SIR PF as a function of the sample size, M .

Numerical example II Hybrid filter: D := D ESRF D ETPF . Figure: RMSEs for hybrid ESRF ( α = 0 ) and 2nd-order corrected LETF/ETPF ( α = 1 ) as a function of the sample size, M .

Data Assimilation Center Potsdam ◮ PhD and Postdoc positions in 11 projects ◮ Mathematical foundation and applications in geophysics, biologie and cognitive sciences ◮ Support for external projects and long term visitors

Data Assimilation Center Potsdam ◮ PhD and Postdoc positions in 11 projects ◮ Mathematical foundation and applications in geophysics, biologie and cognitive sciences ◮ Support for external projects and long term visitors Thank you!

References I Nawinda Chustagulprom, Sebastian Reich, and Maria Reinhardt. A hybrid ensemble transform filter for nonlinear and spatially extended dynamical systems. SIAM/ASA J Uncertainty Quantification , 4:592–608, 2016. J. de Wiljes, W. Acevedo, and S. Reich. A second-order accurate ensemble transform particle filter. accepted in SIAM Journal on Scientific Computing , (ArXiv:1608.08179), 2017. J. de Wiljes, S. Reich, and W. Stannat. Long-time stability and accuracy of the ensemble Kalman-Bucy filter for fully observed processes and small measurement noise. Technical Report https://arxiv.org/abs/1612.06065, 2016. S. Reich and C. Cotter. Probabilistic Forecasting and Bayesian Data Assimilation . Cambridge University Press, 2015.

References II Sebastian Reich. A non-parametric ensemble transform method for Bayesian inference. SIAM J Sci Comput , 35:A2013–A2014, 2013.