Algorithms for high-dimensional non-linear filtering and smoothing - PowerPoint PPT Presentation

Algorithms for high-dimensional non-linear filtering and smoothing problems Jana de Wiljes, Sahani Pathiraja, Sebastian Reich Kobe 2019

Setting Model: x k +1 = f ( x k , . . . , x k − n ) + ǫ k , ǫ k ∼ N (0 , Q ) (1) where x k +1 ∈ R N x and Q ∈ R N x × N x . Observations: (2) y k +1 = h ( x k +1 ) + ν k , ν k ∼ N (0 , R ) where y k +1 ∈ R N y and R ∈ R N y × N y . 1

Filtering and Smoothing Problem Goal: Approximate • p ( x 0: T | y 1: T ) (joint smoothing density) • p ( x k | y 1: T ) (marginal smoothing density) • p ( x k | y 1: k ) (filtering density) Here focus on: Fixed-lag smoothing: p ( x k − L : k | y 1: k ) (3) Remark: • For L = 0 Fixed-lag smoothing = filtering • For L > 0 Fixed-lag smoothing can be interpreted as filtering with augmented state space 2

State of the art Linear Model : Kalman Filter/Smoother • Gaussian posterior • deterministic formulas Nonlinear Model : approximate with empirical measure, i.e., M � 1 M δ ( x − x ( i ) p ( x 0: T | y 1: T ) ≈ 0: T | 1: T ) . i =1 → This approach leads to a variety of smoothers • Particle smoother • Ensemble Kalman smoother 3

Particle collapse Degeneracy y k of Particle Smoother x ( i ) k − L : k w ( i ) k − L : k > 0 4

Alternative approach Move particles y k towards observation x ( i ) k − L : k D k − L : k | k 5

Linear Ensemble Transform Smoothers

Linear Ensemble Transform Smoothers Update via linear transformation: X k − L : k | k = X k − L : k | k − 1 D k − L : k | k = X k − L : k | k − L D k − L : k | k − L +1 · · · · · D k − L : k | k where X k − L : k | k = [ x (1) k − L : k | k , ..., x ( M ) k − L : k | k ] ∈ R N X L × M and x ( i ) k − L : k | s ∼ p ( x k − L : k | y 1: s ) (4) with s ∈ { k − L , . . . , k } 6

Motivation: LETS Benefits: • localisation [Reich, 2013][Poterjoy, 2015] • hybrid [Frei and K¨ unsch, 2013][Chustagulprom et al., 2016] • moment matching [T¨ odter and Ahrens, 2015, Bishop et al., 2001, Xiong et al., 2006, Lei and Bickel, 2011, Acevedo et al., 2017] Motivation for Hybrid formulations • Ensemble Kalman Filters (EnKFs) + Robust and moderately affordable – Gaussian assumptions • Traditional Particle Filters (PFs) + No Gaussianity assumptions – Liable to the ”Curse of Dimensionality” + consistent in the ensemble limit 7

Example: Ensemble Kalman smoother Analysis step: X k | k = X k | k − 1 D k | k with 1 M − 1( x ( i ) { D k | k } ij = δ ij − k | k − 1 − m k | k − 1 ) · H ⊤ ( H P k | k − 1 H ⊤ + R ) − 1 ( H x ( j ) k | k − 1 + ξ ( j ) − y k ) Pros and Cons : + Robust and moderately affordable + Works well in practice – Gaussian assumptions – Mathematical foundation for nonlinear models largely missing 8

Ensemble Square Root Smoother Transform coefficients: k − L : k | k − 1 w ( i ) ( D ESRS k − L : k | k ) ij := ( S k − L : k | k ) ij + ˆ M , (5) with Square root matrix � � − 1 / 2 1 HA k − L : k | k − 1 ) T R − 1 ˜ M − 1(˜ S k − L : k | k = I + HA k − L : k | k − 1 , where � � 1 x (1) m k − L : k | k − 1 , . . . , x ( M ) √ k − L : k | k − 1 − � k − L : k | k − 1 − � A k − L : k | k − 1 = m k − L : k | k − 1 M − 1 and k − L : k | k = 1 1 w ( i ) i S 2 k − L : k | k (˜ HA k − L : k | k − 1 ) T R − 1 (˜ M − 1 e T ˆ M − H � m k − L : k | k − 1 − y k ) 9

Nonlinear Ensemble Transform Smoother [T¨ odter and Ahrens, 2015] NETF: √ D NETS M [ W k | k − w k | k ( w k | k ) T ] 1 / 2 Ω k − L : k | k = w k | k 1 + (6) where k | k ) T ∈ R M × 1 w = ( w (1) k | k , . . . , w ( M ) and W k | k = diag ( w k | k ) 10

Ensemble Transform Particle Smoother [de Wiljes et al., 2019] Given: • M samples x ( i ) k − L : k | k − 1 ∼ p ( x k − L : k | y 1: k − 1 ) (prior ensemble) • normalized importance weights w ( j ) = p ( y k | x ( j ) k − L : k | k − 1 ) (likelihood) k Ansatz: replace resampling step with linear transformation by interpreting it as discrete Markov chain given by transition matrix [Reich, 2013] D ETPS k − L : k | k ∈ R M × M subject to: = ( D ETPS ( A ) d ETPS k − L : k | k ) ij ≥ 0 ∀ i , j ij ( B ) D ETPS k − L − k | k 1 = M w k (7) ( C ) ( D ETPS k − L : k | k ) T 1 = 1 11

Ensemble Transform Particle Smoother [de Wiljes et al., 2019] Then the posterior ensemble members are distributed according to the columns of the transformation   d ETPS 1 j   d ETPS �   X ( j ) ∼ 2 j and x ( j ) x ( i ) ˜   k − L : k | k − 1 = E [ ˜ X ( j ) ] = k − L : k | k − 1 d ETPS .   ij .  .  i d ETPS Mj Example. Monomial resampling   w (1) w (1) w (1) · · ·   w (2) w (2) w (2) · · ·   D Mono   k − L : k | k := w ⊗ 1 = . . . . ...   . . .  . . .  w ( M ) w ( M ) w ( M ) · · · Here: X 1: k | k − 1 and X 1: k | k are independent. Idea: increase correlation between X 1: k | k − 1 and X 1: k | k 12

Ensemble Transform Particle Smoother [de Wiljes et al., 2019] Idea: increase correlation between X 1: k | k − 1 and X 1: k | k Underlying ansatz: approximate a transfer map between X 1: k | k − 1 and X 1: k | k Note: choose transfer map that induces a coupling µ ∗ Z that minimizes � µ ∗ = arg E ||X 1: k | k − 1 − X 1: k | k || 2 inf µ ∈ Π( p X 1: k | k − 1 , p X 1: k | k ) 13

Ensemble Transform Particle Smoother [de Wiljes et al., 2019] Discretization: Solve optimization problem M � k − L : k | k ) ij || x ( i ) k − L : k | k − 1 − x ( j ) D ETPS ( D ETPS k − L : k | k − 1 || 2 k − L : k | k = arg min i , j =1 to find transformation matrix D ETPS k − L : k | k that increases correlation ([Reich, 2013]). Remarks: • consistent for M → ∞ • maximizes correlation between prior and posterior ensembles • is deterministic, preserving the regularity of the state fields 14

OT Solvers Exact solution • Computational complexity O ( M 3 ln( M )) • Efficient Earth Mover Distance algorithms available, e.g. FastEMD [Pele and Werman, 2009]. Entropic Regularized Approximation [Cuturi, 2013] M � k − L : k | k ) ij || x ( i ) k − L : k | k − 1 − x ( j ) ( D ETPS k − L : k | k − 1 || 2 J ( D ) = i , j =1 + 1 λ ( D ETPS k : k + L | k + L ) ij ln( D ETPS k : k + L | k + L ) ij • Sinkhorn’s fixed point iteration can be used. • Computational complexity O ( M 2 · C ( λ )) 1D Approximation(Each variable independently updated) • OT problem reduces to reordering. • Computational complexity O ( M ln( M )) 15 • No particle distance needed

Combating challenges of nonlinear smoothing

Adaptive spread correction and rotation

Moment Matching Idea: design filter to match specific moments for a finite ensemble [Xiong et al., 2006][Lei and Bickel, 2011][T¨ odter and Ahrens, 2015] [Acevedo et al., 2017] Linear case: Kalman Formula m k − L : k | k = m k − L : k | k − 1 − K ( H m k − L : k | k − 1 − y k ) (8) P k − L : k | k = P k − L : k | k − 1 − KH P k − L : k | k − 1 (9) Nonlinear case: Particle filter, e.g., first and second moment M � w ( i ) k x ( i ) m k − L : k | k = k − L : k | k − 1 i =1 � M w ( i ) k ( x ( i ) k − L : k | k − m k − L : k | k )( x ( i ) k − L : k | k − m k − L : k | k ) T P k − L : k | k = i =1 16

First-order accurate LETSs LETS is first-order accurate if � M w ( i ) k x ( i ) m k − L : k | k = (10) k − L : k | k − 1 i =1 is equal to � M m k − L : k | k = 1 x ( i ) � (11) k − L : k | k M i =1 17

First-order accurate LETSs Class of first-order accurate LETSs D 1 = { D ∈ R M × M | D T 1 = 1 , D1 = M w } Examples : • D EnKS / ∈ D 1 • D ESRS / ∈ D 1 • D ETPS ∈ D 1 • D Mono ∈ D 1 • D NETS ∈ D 1 18

First-order accurate LETSs Linear transformations D λ replacing resampling step : M � k − L : k | k ) ij || x ( i ) k − L : k | k − 1 − x ( j ) D λ ( D λ k − L : k | k − 1 || 2 k − L : k | k = arg min (12) i , j =1 + 1 λ ( D λ k − L : k | k ) ij ln( D λ k − L : k | k ) ij (13) for given λ > 0 subject to � � 1 ( D λ ( D λ ( D λ k − L : k | k ) ij = w ( i ) k − L : k | k ) ij ≥ 0 , k − L : k | k ) ij = 1 , k | k M i j Remark. • λ → 0: D 0 = w ⊗ 1 (monomial resampling). • λ → ∞ : D ∞ solves the optimal coupling/transport problem. • Effective iterative solvers are available [Cuturi, 2013]. 19

Moment Matching The analysis covariance matrix � M P k − L : k | k = 1 � ( x ( i ) m k − L : k | k )( x ( i ) m k − L : k | k ) T . k − L : k | k − � k − L : k | k − � M i =1 is equal to the covariance matrix defined by the importance weights, i.e. � M w ( i ) k | k ( x ( i ) k − L : k | k − m k − L : k | k )( x ( i ) k − L : k | k − m k − L : k | k ) T P k − L : k | k = i =1 → the following equation has to be satisfied ( D k − L : k | k − w k | k 1 T )( D k − L : k | k − w k | k 1 T ) T = W k | k − w k | k w T k | k . 20

Second-order corrected LETSs Given: D ∈ D 1 Goal: correct transformation � D ∈ D 2 Ansatz: � D = D + ∆ with ∆ ∈ R M × M such that ∆1 = 0 , ∆ T 1 = 0 ([Acevedo et al., 2017]). 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 21