The 7th International Symposium on Data Assimilation (ISDA2019) - PowerPoint PPT Presentation

The 7th International Symposium on Data Assimilation (ISDA2019) Efficient Implementations of Ensemble Based Methods In Sequential Data Assimilation: Accounting for Localization Elias D. Ni˜ no-Ruiz Applied Math and Computer Science Laboratory (AML-CS) Department of Computer Science Universidad del Norte BAQ 080001, Colombia E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Outline I Data Assimilation Components Ensemble Based Methods The Stochastic Ensemble Kalman Filter Localization Methods Precision Matrix Localization Efficient EnKF-MC Shrinkage Covariance Matrix Estimation Ensemble Kalman filter based on RBLW Efficient Implementation of the RBLW EnKF-RBLW EnKF-MC and EnKF-RBLW with the SPEEDY Model Accuracy of the EnKF-MC Local Estimation of B − 1 Accuracy of the EnKF-RBLW Parallel Implementations of Ensemble Based Methods Recent References References E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Components in DA [BS12] I � 10 8 � ◮ We want to estimate x ∗ ∈ ❘ n × 1 . n ∼ O . ◮ Imperfect numerical model: x next = M t current → t next ( x current ) , where x ∈ ❘ n × 1 . ◮ Noisy observations: y = H ( x ) + ǫ ∈ ❘ m × 1 , � 10 6 � where H : ❘ n → ❘ m and ǫ ∼ N ( 0 m , R ). m ∼ O . ◮ Prior estimate x b ∈ ❘ n × 1 with errors following N ( 0 , B ). E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Components in DA [BS12] II 0.5 0.6 20 20 0.7 40 0.8 40 0.9 60 60 1 1.1 80 80 1.2 100 100 1.3 1.4 120 120 1.5 20 40 60 80 100 120 2000 4000 6000 8000 10000 12000 14000 16000 20 40 60 80 100 120 (c) y = H · x ∗ + ǫ (a) x ∗ (b) H 20 40 60 80 100 120 20 40 60 80 100 120 (d) x b E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Components in DA [BS12] III ◮ By Bayes’ Theorem we know that: P ( x | y ) ∝ P ( x ) · L ( x | y ) where � � � � x − x b � − 1 2 � � P ( x ) ∝ exp 2 · � B − 1 � � − 1 2 · � y − H · x � 2 L ( x | y ) ∝ exp R − 1 and therefore, x a = arg max P ( x | y ) , x E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Components in DA [BS12] IV ◮ It can be easily shown that: � � x b + A · H T · R − 1 · d = A · B − 1 · x b + H T · R − 1 · y x a = � R + H · B · H T � − 1 x b + B · H T · · d = � � − 1 ∈ ❘ n × n , and B − 1 + H T · R − 1 · H where A = d = y − H · x b ∈ ❘ m × 1 . ◮ Posterior distribution: x ∼ N ( x a , A ) . ◮ How do we estimate x b and B ?. E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Ensemble Based Methods ◮ We can make use of an ensemble of model realizations: � x b [1] , x b [2] , . . . , x b [ N ] � X b = ∈ ❘ n × N ◮ Empirical moments of the ensemble: 1 x b ≈ x b N · X b · 1 N ∈ ❘ n × n , = 1 N − 1 · δ X · δ X T , B ≈ P b = and δ X = X b − x b · 1 T N ∈ ❘ n × N . E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

The Lorenz 96 Model - Toy Model I ◮ The Lorenz 96 model:   ( x 2 − x n − 1 ) · x n − x 1 + F for i = 1 ,  dx j dt = ( x i +1 − x i − 2 ) · x i − 1 − x i + F for 2 ≤ i ≤ n − 1 , (1)   ( x 1 − x n − 2 ) · x n − 1 − x n + F for i = n , where x i stands for the i -th model component, for 1 ≤ i ≤ n , usually n = 40. ◮ Each model component stands for a particle which fluctuates in the atmosphere. ◮ Exhibits chaotic behaviour when the external force F is set to 8. E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

The Lorenz 96 Model - Toy Model II 15 15 15 10 10 10 5 5 5 0 0 0 -5 -5 -5 -10 -10 -10 0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10 (e) x 5 (f) x 10 (g) x 20 15 15 15 10 10 10 5 5 5 0 0 0 -5 -5 -5 -10 -10 -10 0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10 (h) x 30 (i) x 35 (j) x 40 E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Estimation of B via N = 10 5 . 5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 (a) Structure (b) Surf Figure: Estimation of B via N = 10 5 . E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

The Stochastic Ensemble Kalman Filter [Eve03, Eve06] I ◮ Sequential Monte Carlo method for parameter and state estimation. ◮ Analysis ensemble (posterior ensemble): � � X b + P b · H T · R + H · P b · H X a · ∆Y = X b + P a · H T · R − 1 · ∆Y ∈ ❘ n × N , X a = � � � P b � − 1 P a · H T · R − 1 · Y s + X a · X b ∈ ❘ n × N , = � P b � − 1 � − 1 � H T · R − 1 · H + where P a = ∈ ❘ n × n , and the e -th column of ∆Y ∈ ❘ m × N and Y s ∈ ❘ n × N are: � x b [ e ] � d [ e ] = y + ǫ [ e ] − H ∈ ❘ m × 1 , and y s [ e ] = y + ǫ [ e ] , respectively, for 1 ≤ e ≤ N , and ǫ [ e ] ∼ N ( 0 m , R ). E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

L − 2 Error Norms in Time, N = 10 5 2 1.5 1 1 0 0.5 -1 0 -0.5 -2 0 5 10 15 0 5 10 15 (a) p = 50% (b) p = 100% Figure: L − 2 error norms in time, N = 10 5 . But too many samples!!! In practice, model realizations are constrained by the hundreds... E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

L − 2 error norms in time, N = 10 1.5 1.5 1.48 1.46 1.45 1.44 1.4 1.42 1.4 1.38 1.35 0 5 10 15 0 5 10 15 (a) p = 50% (b) p = 100% Figure: L − 2 error norms in time, N = 10. What is going on here? ... E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Estimation of B via N = 10 5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 (a) Structure (b) Surf Figure: Estimation of B via N = 10. What can we do? Localization methods... E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Localization Methods ◮ Avoid the impact of spurious correlations. ◮ Increase the rank of P b . ◮ Three different flavors: 1. Covariance Matrix Localization. (Precision Localization) [NRSD15, NRSD17, NR17, NRSD18]. 2. Spatial Domain Localization [OHS + 04]. 3. Observation Localization [AND07, AND09]. E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Precision Matrix Localization I ◮ Component-wise products are prohibitive in high-dimensional spaces. ◮ When two model components are conditional independent, their corresponding entry in the precision covariance matrix is zero . (a) r = 0 (b) r = 1 (c) r = 3 E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Precision Matrix Localization II ◮ Modified Cholesky Decomposition [BL + 08]: B − 1 = T T · D − 1 · T � where the non-zero elements from T ∈ ❘ n × n are given by fitting models of the form: � x [ i ] = x [ q ] · {− T } i , q + ǫ [ i ] ∈ ❘ N × 1 , for 1 ≤ i ≤ n , q ∈ P ( i , r ) � ǫ [ i ] � and { D } i , i = var . (a) N (6 , 1) (b) P (6 , 1) E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Precision Matrix Localization III ◮ An estimate: 0 0 5 5 10 10 15 15 20 20 25 25 30 30 35 35 40 40 0 10 20 30 40 0 10 20 30 40 nz = 160 nz = 40 (a) P b (b) T (c) D 0 5 10 15 20 25 30 35 40 0 10 20 30 40 nz = 298 B − 1 Str (d) � (e) � B − 1 (f) � B Results: E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Precision Matrix Localization IV 2 2 2 1 1 1 0 0 0 -1 -1 -1 -2 -2 -2 0 5 10 15 0 5 10 15 0 5 10 15 (a) N = 30, r = 1, (b) N = 30, r = 3, (c) N = 30, r = 5, p = 100% p = 100% p = 100% 1.4 1.5 1.4 1.2 1 1.2 1 0.5 0.8 1 0 0.6 0.8 -0.5 0.4 0.2 -1 0.6 0 5 10 15 0 5 10 15 0 5 10 15 (d) N = 30, r = 1, (e) N = 30, r = 3, (f) N = 30, r = 5, p = 50% p = 50% p = 50% E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Efficient EnKF-MC I Consider � � − 1 X b + B − 1 + H T · R − 1 / 2 · R − 1 / 2 · H · H T · R − 1 · ∆Y � X a = � B − 1 + Z · Z T � − 1 · H T · R − 1 · ∆Y X b + � =   − 1 � z [ j ] � T m � X b + B − 1 + z [ j ] · · H T · R − 1 · ∆Y , �  = j =1 z [ j ] ∈ ❘ n × 1 is the j -th column of Z = H T · R − 1 / 2 ∈ ❘ n × m . E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Efficient EnKF-MC II � T (0) � T � T (0) � · D (0) · = T T · D · T = � B − 1 , A (0) = � z [1] � T � T (1) � T � T (1) � A (0) + z [1] · · D (1) · A (1) = = , � z [2] � T � T (2) � T � T (2) � A (1) + z [2] · · D (2) · A (2) = = , . . . � z [ m ] � T A ( m − 1) + z [ m ] · A ( m ) = � T ( m ) � T � T ( m ) � · D ( m ) · T T · � T = A − 1 , = � D · � = E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS

Efficient EnKF-MC III at any intermediate step j , for 1 ≤ j ≤ m , we have, � T ( j − 1) � T � T ( j − 1) � � z [ j ] � T · D ( j − 1) · + z [ j ] · A ( j ) = � p ( j ) � T � � T ( j − 1) � T � � T ( j − 1) � D ( j − 1) + p ( j ) · · · = , � T ( j − 1) � T · p ( j ) = z [ j ] ∈ ❘ n × 1 . By computing the Cholesky where decomposition of, � p ( j ) � T � T ( j − 1) � T � T ( j − 1) � D ( j − 1) + p ( j ) · · D ( j ) · � � = , therefore, � T ( j − 1) · T ( j − 1) � T � T ( j − 1) · T ( j − 1) � · D ( j ) · � � A ( j ) = � T ( j ) � T � T ( j ) � · D ( j ) · = , E. Ni˜ no-Ruiz, ISDA2019 - RIKEN R-CCS