Ergodicity in data assimilation methods David Kelly Andy Majda Xin Tong Courant Institute New York University New York NY www.dtbkelly.com April 7, 2016 Statistics and Applied Math seminar, U Chicago . David Kelly (CIMS) Data assimilation April 7, 2016 1 / 34
Part I: What is data assimilation? David Kelly (CIMS) Data assimilation April 7, 2016 2 / 34
What is data assimilation? Suppose u satisfies du dt = F ( u ) with some unknown initial condition u 0 . We are most interested in geophysical models, so think high dimensional, nonlinear, possibly stochastic. Suppose we make partial, noisy observations at times t = h , 2 h , . . . , nh , . . . y n = Hu n + ξ n where H is a linear operator (think low rank projection), u n = u ( nh ), and ξ n ∼ N (0 , Γ) iid. The aim of data assimilation is to say something about the conditional distribution of u n given the observations { y 1 , . . . , y n } David Kelly (CIMS) Data assimilation April 7, 2016 3 / 34
Illustration (Initialization) y Ψ Figure: The blue circle represents our guess of u 0 . Due to the uncertainty in u 0 , this is a probability measure. x obs David Kelly (CIMS) Data assimilation April 7, 2016 4 / 34
Illustration (Forecast step) y Ψ Figure: Apply the time h flow map Ψ. This produces a new probability measure which is our forecasted estimate of u 1 . This is x called the forecast step. obs David Kelly (CIMS) Data assimilation April 7, 2016 4 / 34
Illustration (Make an observation) y Ψ Figure: We make an observation y 1 = Hu 1 + ξ 1 . In the picture, we only observe the x variable. x obs David Kelly (CIMS) Data assimilation April 7, 2016 4 / 34
Illustration (Analysis step) y Figure: Using Bayes formula we compute the conditional distribution of u 1 | y 1 . This new Ψ measure (called the posterior) is the new estimate of u 1 . The uncertainty of the estimate is reduced by x incorporating the obs observation. The forecast distribution steers the update from the observation. David Kelly (CIMS) Data assimilation April 7, 2016 4 / 34
Bayes’ formula filtering update Let Y n = { y 0 , y 1 , . . . , y n } . We want to compute the conditional density P ( u n +1 | Y n +1 ), using P ( u n | Y n ) and y n +1 . By Bayes’ formula, we have P ( u n +1 | Y n +1 ) = P ( u n +1 | Y n , y n +1 ) ∝ P ( y n +1 | u n +1 ) P ( u n +1 | Y n ) But we need to compute the integral � P ( u n +1 | Y n ) = P ( u n +1 | Y n , u n ) P ( u n | Y n ) du n . David Kelly (CIMS) Data assimilation April 7, 2016 5 / 34
The ‘optimal filtering’ framework is impractical for high dimensional models, as the integrals are impossible to compute and densities impossible to store. In numerical weather prediction , we have O (10 9 ) variables for ocean-atmosphere models (discretized PDEs). David Kelly (CIMS) Data assimilation April 7, 2016 6 / 34
The Kalman Filter For a linear model u n +1 = Mu n + η n +1 , the Bayesian integral is Gaussian and can be computed explicitly. The conditional density is characterized by its mean and covariance m n +1 = (1 − K n +1 H ) � m n + K n +1 y n +1 C n +1 = ( I − K n +1 H ) � C n +1 , where m n +1 , � • ( � C n +1 ) is the forecast mean and covariance. C n +1 H T ) − 1 is the Kalman gain . • K n +1 = � C n +1 H T (Γ + H � The procedure of updating ( m n , C n ) �→ ( m n +1 , C n +1 ) is known as the Kalman filter . David Kelly (CIMS) Data assimilation April 7, 2016 7 / 34
Extended Kalman filter Suppose we have a nonlinear model: u n +1 = Φ( u n ) + Σ 1 / 2 η n +1 where Φ is a nonlinear map, η n Gaussian. The Extended Kalman filter is given by the same update formulas m n +1 = (1 − K n +1 H ) � m n +1 + K n +1 y n +1 C n +1 = ( I − K n +1 H ) � C n +1 , C n +1 = D Φ( m n ) C n D Φ( m n ) T + Σ. m n +1 = Φ( m n ) and � where � Thus we approximate the forecast distribution with a Gaussian. Still too expensive for O (10 9 ) variables.... David Kelly (CIMS) Data assimilation April 7, 2016 8 / 34
Ensemble Kalman filter (Evensen 94) y Ψ Figure: Start with K ensemble members drawn from some distribution. Empirical representation of u 0 . The ensemble members x are denoted u ( k ) 0 . obs Only KN numbers are stored. Better than Kalman if K < N . David Kelly (CIMS) Data assimilation April 7, 2016 9 / 34
Ensemble Kalman filter (Forecast step) y Ψ Figure: Apply the dynamics Ψ to each ensemble member. x obs David Kelly (CIMS) Data assimilation April 7, 2016 9 / 34
Ensemble Kalman filter (Make obs) y Ψ Figure: Make an observation. x obs David Kelly (CIMS) Data assimilation April 7, 2016 9 / 34
Ensemble Kalman filter (Analysis step) y Ψ Figure: Approximate the forecast distribution with a Gaussian. Fit the Gaussian using the empirical statistics of x the ensemble. obs David Kelly (CIMS) Data assimilation April 7, 2016 9 / 34
How to implement the Gaussian approximation The naive method is to simply write: P ( y 1 | u 1 ) P ( u 1 ) ∝ exp( − 1 2 | Γ − 1 / 2 ( y 1 − Hu 1 ) | 2 ) exp( − 1 − 1 / 2 ( u 1 − � 2 | � m 1 ) | 2 ) C with the empirical statistics K � m 1 = 1 Ψ ( k ) ( u ( k ) � 0 ) K k =1 � � � � T � K 1 Ψ ( k ) ( u ( k ) Ψ ( k ) ( u ( k ) � C 1 = 0 ) − � m 1 0 ) − � m 1 . K − 1 k =1 In the linear model case Ψ( u n ) = Mu n + η n , this produces an unbiased estimate of the posterior mean, but a biased estimate of the covariance. David Kelly (CIMS) Data assimilation April 7, 2016 10 / 34
How to implement the Gaussian approximation A better approach is to sample using Randomized Maximum Likelihood (RML) method: Draw the sample u ( k ) by minimizing the functional 1 1 − Hu ) | 2 + 1 − 1 / 2 2 | Γ − 1 / 2 ( y ( k ) ( u − Ψ ( k ) ( u ( k ) 2 | � 0 )) | 2 C 1 1 where y ( k ) = y 1 + ξ ( k ) is a perturbed observation. 1 1 In the linear case Ψ( u n ) = Mu n + η n , this produces iid Gaussian samples with mean and covariance satisfying the Kalman update equations, with � C in place of the true forecast covariance. We end up with u ( k ) = (1 − K 1 H )Ψ ( k ) ( u ( k ) 0 ) + K 1 Hy ( k ) 1 1 David Kelly (CIMS) Data assimilation April 7, 2016 11 / 34
Ensemble Kalman filter (Perturb obs) y Ψ Figure: Turn the observation into K artificial observations by perturbing by the same source of observational x noise. obs y ( k ) = y 1 + ξ ( k ) 1 1 David Kelly (CIMS) Data assimilation April 7, 2016 12 / 34
Ensemble Kalman filter (Analysis step) y Ψ Figure: Update each member using the Kalman update formula. The Kalman gain K 1 is computed using the x ensemble covariance. obs u ( k ) = (1 − K 1 H )Ψ( u ( k ) 0 ) + K 1 Hy ( k ) K 1 = � C 1 H T (Γ + H � C 1 H T ) − 1 1 1 � K 1 (Ψ( u ( k ) m n +1 )(Ψ( u ( k ) � m n +1 ) T C 1 = 0 ) − � 0 ) − � K − 1 David Kelly (CIMS) k =1 Data assimilation April 7, 2016 12 / 34
Part II: Ergodicity for DA methods: Why are we interested in ergodicity? And what kind? David Kelly (CIMS) Data assimilation April 7, 2016 13 / 34
Two types of ergodicity 1 - Signal-Filter Ergodicity : Ergodicity for the homogeneous Markov chain ( u n , u (1) n , . . . , u ( K ) ). n n ( · ; u (1) 0 , . . . , u ( K ) Conditional ergodicity : Let P Y n 2 - ) be the law of 0 ( u (1) n , . . . , u ( K ) ), given the observations Y n and initialization n ( u (1) 0 , . . . , u ( K ) ). We call the filter conditionally ergodic when any two 0 such measures with different initialization, but same observations Y n , converge as n → ∞ . David Kelly (CIMS) Data assimilation April 7, 2016 14 / 34
Two types of ergodicity Signal-filter ergodicity tells us that: • The filter will not blow-up (catastrophic filter divergence). • The long-time statistics of the filter are stable with respect to initialization. • The filter inherits ‘physical properties’ from the underlying model. Conditional ergodicity suggests that the method is accurate , since all measures are synchronizing, they should be shadowing the true trajectory. David Kelly (CIMS) Data assimilation April 7, 2016 15 / 34
Animation I The model du = −∇ V ( u ) dt + σ dW , where V ( x , y ) = (1 − x 2 − y 2 ) 2 . Observation u x ( t ) + ξ ( t ). For EnKF, the process is signal-filter ergodic, but it is not conditionally ergodic. David Kelly (CIMS) Data assimilation April 7, 2016 16 / 34
Animation II The model du = −∇ V ( u ) dt + σ dW where V ( x , y ) = (1 − x 2 − y 2 ) 2 . Observation y = u 1 + ξ . This filter (a type of particle filter) is both signal-filter ergodic and conditionally ergodic. David Kelly (CIMS) Data assimilation April 7, 2016 17 / 34
Today we focus on signal-filter ergodicity . Conditional ergodicity is much more difficult (but work in progress!). David Kelly (CIMS) Data assimilation April 7, 2016 18 / 34
Recommend
More recommend
