A Sampling Filter for Non-Gaussian Data Assimilation

Computer Science T ec hnical Rep ort CSTR-4/2014 Decem b er 9, 2014 Ahmed Attia and Adrian Sandu “ A Hybrid Monte Carlo Sampling Filter for Non-Gaussian Data Assimilation ” Computational Science Lab oratory Computer Science Departmen t Virginia P olytec hnic Institute and State Univ ersit y Blac ksburg, V A 24060 Phone: (540)-231-2193 F ax: (540)-231-6075 Email: sandu@cs.vt.edu W eb: http://csl.cs.vt.edu Innovative Computational Solutions A Hybrid Mon te Carlo Sampling Filter for Non-Gaussian Data Assimilation Ahmed A ttia a , Adrian Sandu a a Computational Scienc e L ab or atory Dep artment of Computer Scienc e Vir ginia Polyte chnic Institute and State University 2201 Know le dgeworks II, 2202 Kr aft Drive, Blacksbur g, V A 24060, USA Phone: 540-231-2193, F ax: 540-231-9218 E-mail: sandu@cs.vt.e du Abstract Data assimilation combines information from mo dels, measurements, and priors to estimate the state of a dynamical system suc h as the atmosphere. The Ensemble Kalman ﬁlter (EnKF) is a family of ensemble-based data assimilation approaches that has gained wide p opularit y due its simple form ulation, ease of implemen tation, and go od practical results. Most EnKF algorithms assume that the underlying probabilit y distributions are Gaussian. Although this assumption is well accepted, it is to o restrictiv e when applied to large non- linear models, nonlinear observ ation operators, and large levels of uncertain ty . Sev eral ap- proac hes ha ve b een proposed in order to av oid the Gaussianit y assumption. One of the most successful strategies is the maxim um likelihoo d ensem ble ﬁlter (MLEF) which computes a maxim um a p osteriori estimate of the state assuming the p osterior distribution is Gaussian. MLEF is designed to w ork with nonlinear and ev en non-diﬀeren tiable observ ation op erators, and shows go od practical p erformance. Ho w ever, there are limits to the degree of nonlin- earit y that MLEF can handle. This pap er prop oses a new ensemble-based data assimilation metho d, named the “ sampling ﬁlter ”, which obtains the analysis by sampling directly from the p osterior distribution. The sampling strategy is based on a Hybrid Monte Carlo (HMC) approac h that can handle non-Gaussian probabilit y distributions. Numerical exp erimen ts are carried out using the Lorenz-96 mo del and observ ation op erators with diﬀerent lev els of non-linearity and diﬀeren tiabilit y . The prop osed ﬁlter is also tested with shallo w water mo del on a sphere with linear observ ation op erator. The results show that the sampling ﬁlter can p erform w ell ev en in highly nonlinear situations w ere EnKF and MLEF ﬁlters div erge. Keywor ds: Data assimilation, v ariational metho ds, ensem ble ﬁlters, Mark ov chain, hybrid Mon te-Carlo Pr eprint submitte d to QJRMS De c emb er 9, 2014 Con ten ts 1 In tro duction 3 2 Data Assimilation 4 2.1 Problem form ulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 The ensem ble Kalman ﬁlter . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 The maxim um lik eliho o d ensem ble ﬁlter . . . . . . . . . . . . . . . . . . . . 6 3 Hybrid Mark ov Chain Mon te Carlo 7 3.1 Hamiltonian dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 HMCMC sampling algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4 The Sampling Filter for Data Assimilation 10 5 Numerical Results 13 5.1 The Lorenz-96 mo del . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 5.2 Observ ations and observ ation op erators . . . . . . . . . . . . . . . . . . . . . 13 5.3 Exp erimen tal setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.4 Linear observ ation op erator exp erimen ts . . . . . . . . . . . . . . . . . . . . 15 5.5 Quadratic observ ation op erator exp erimen ts . . . . . . . . . . . . . . . . . . 18 5.6 Cubic observ ation op erator exp erimen ts . . . . . . . . . . . . . . . . . . . . . 20 5.7 Absolute v alue observ ation op erator exp erimen ts . . . . . . . . . . . . . . . . 23 5.8 Quadratic observ ation op erator with threshold exp erimen ts . . . . . . . . . . 25 5.9 Exp onen tial observ ation op erator (with factor r = 0 . 2) exp erimen ts . . . . . 28 5.10 MLEF performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.11 T uning the num b er of MC steps b et ween successive state selections . . . . . 37 5.12 A highly nonlinear observ ation op erator . . . . . . . . . . . . . . . . . . . . 39 5.13 Shallow w ater mo del on a sphere . . . . . . . . . . . . . . . . . . . . . . . . 42 5.14 Results for shallow water mo del with linear observ ations . . . . . . . . . . . 44 6 Conclusions and F uture W ork 45 A Symplectic n umerical integrators 46 A.1 P osition V erlet in tegrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.2 Tw o-stage in tegrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.3 Three-stage in tegrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.4 F our-stage integrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 A.5 General in tegrator deﬁned on Hilb ert space . . . . . . . . . . . . . . . . . . . 48 2 1. In tro duction Data assimilation is the pro cess of com bining information from models, measurements, and priors - all with asso ciated uncertain ties - in order to obtain the b est estimate of the state of a physical system. Tw o families of metho ds, v ariational and ensemble based ﬁlters, ha v e pro ved very successful in real applications. V ariational metho ds, rooted in control theory , require costly developmen ts of tangent linear and adjoin t mo dels [ 20 ]. Ensem ble- based sequential data assimilation sc hemes are ro oted in statistical estimation theory . The ensem ble Kalman Filter was introduced by Evensen [ 10 ] and has undergone considerable dev elopmen ts since then. EnKF form ulations fall in one of tw o classes, namely sto chastic or deterministic form ulations [ 35 ]. In the stochastic approach, eac h ensem ble mem b er is up dated using a p erturb ed version of the observ ation vector [ 6 , 16 ]. In the deterministic form ulation (which leads to square ro ot ensemble ﬁlters [ 1 , 4 , 30 , 35 , 36 ] no observ ation noise is added, but transformations of the co v ariance matrix are applied such as to recov er the correct analysis statistics. All v ariants of the EnKF work w ell in case of linear observ ations [ 12 ], how ever in real applications the observ ation op erators are in general nonlinear. EnKF can accommo date nonlinear observ ation op erators using linearization, in the spirit of the extended Kalman ﬁlter [ 37 ]. An alternative approach to handle the non-linearit y of observ ation op erators is to use the diﬀerence b etw een nonlinear op erators ev aluated at t w o states instead of the linearized v ersion; this approach can result in mathematical inconsistencies [ 37 ]. A diﬀeren t approac h to deal with nonlinear observ ations is to p ose a nonlinear estimation problem in a subspace spanned b y the ensem ble mem b ers, and to compute the maxim um a p osteriori estimate in that subspace. This leads to the maximum lik eliho o d ensemble ﬁlter (MLEF) prop osed by Zupanski [ 37 ]. MLEF minimizes a cost function that dep ends on nonlinear observ ation operators. MLEF do esn’t require the observ ation op erator to b e diﬀeren tiable and uses a diﬀerence approximation of the Jacobian of the observ ation operator. Ho w ever, this approach may diverge if the observ ation op erator is highly nonlinear. In addition it is inheren tly assumed that the p osterior distribution is Gaussian; the MLEF maximum a p osteriori probabilit y estimate ma y face diﬃculties in case of multimodal distributions. The current adv ances in sampling algorithms make it feasible to directly sample from the p osterior probabilit y distribution of the system state. A promising step to w ards eﬃcient sequen tial Mon te Carlo sampling from the p osterior densit y is the implicitly particle ﬁlter [ 7 ]. This algorithm directs the sampling tow ards the regions of high densit y areas in the p osterior. This helps to control the n umber of particles in case of of v ery high dimensional state spaces. The implicit sampling ﬁlter, how ever, is exp ensive: it requires an optimization step for each particle and eac h ensem ble member is generated by solving a set of algebraic equations. This w ork seeks to develop an ensemble-based data assimilation ﬁltering tec hnique that can accommo date non-Gaussian p osterior distributions and can b e eﬃcien tly applied in op- erational situations. Our approach is based on directly sampling the p osterior probabilit y densit y using a Marko v Chain Monte Carlo (MCMC) strategy that generates a Marko v c hain whose inv arian t (stationary) distribution is the target probability densit y . Sp eciﬁcally , we emplo y the h ybrid Marko v Chain Mon te Carlo (HMCMC) algorithm, a v ariant of MCMC 3 sampling that incorp orates an auxiliary v ariable and takes adv antage of the prop erties of Hamiltonian system dynamics [ 9 ].This sampling sc heme turns out to be v ery useful in case of complex high dimensional distribution. The new fully nonlinear sampling ﬁlter can ac- commo date nonlinear observ ation op erators and it do es not require the target probabilit y distribution to b e Gaussian. The pap er is organized as follows. An ov erview of data assimilation problem and widely- used solution strategies is given in Section 2 . Sampling MCMC and HMC algorithms are summarized in Section 3 . The prop osed sampling ﬁlter is presen ted in Section 4 . Numerical exp erimen ts, and a comparison of the sampling ﬁlter against traditional EnKF and MLEF metho ds, are giv en in Section 5 . Conclusions are drawn in Section 6 . 2. Data Assimilation This section pro vides a brief o verview of the data assimilation (D A) problem and of sev eral solution strategies, and highlights the motiv ation b ehind the presen t research. 2.1. Pr oblem formulation Data assimilation com bines information from prior (bac kground) kno wledge, a n umerical mo del, and observ ations, all with asso ciated errors, to obtain a statistically b est estimate of the state x true of a ph ysical system. The background represen ts the b est estimate of the true state prior to any measuremen t b eing a v ailable. The bac kground errors (uncertain ties) are generally assumed to ha ve a Gaussian distribution x b − x true ∈ N (0 , B ), where B is the background error co v ariance matrix. The Gaussian assumption is widely used and we will follow it as well. The n umerical mo del propagates the initial mo del state (initial condition) x 0 ∈ R n v ar at time t 0 to future states x k ∈ R n v ar at times t k : x k = M t 0 → t k ( x 0 ) , t 0 ≤ t k ≤ t F , (1) where t 0 and t F are the b eginning and the end points of the simulation time interv al. The mo del solution op erator M represen ts, for example, a discrete appro ximation of the partial diﬀeren tial equations that go vern the evolution of the dynamical system (e.g., the atmo- sphere). The state space is typically large, e.g., n v ar ∼ 10 6 − 10 9 v ariables for atmospheric sim ulations. Small p erturbations δ x of the state of the system evolv e according to the tangen t linear mo del: δ x k = M t 0 → t k ( x 0 ) · δ x 0 , t 0 ≤ t k ≤ t F , (2) where M = M 0 is the linearized mo del solution op erator. Observ ations of the true state are a v ailable at discrete time instants t k , t 0 ≤ t k ≤ t F , y k = y ( t k ) = H k ( x k ) + ε k , k = 0 , 1 , . . . , n obs − 1 . The observ ation op erator H k maps the state space to the observ ation space at time t k . The observ ations are corrupted by measuremen t and representativ eness errors [ 8 ], which are also 4 assumed to hav e a normal distribution, ε k ∈ N (0 , R k ), where R k is the observ ation error co v ariance matrix at time t k . Data assimilation combines the background estimate, the measurements , and the mo del to obtain an improv ed estimate x a , called the “ analysis ” (or p osterior ), of the true state x true . Two approaches for solving the data assimilation problem hav e gained widespread p opularit y , v ariational and ensem ble-based metho ds. The sampling ﬁlter prop osed in this pap er b elongs to the latter family . W e will compare the new metho dology with tw o existing algorithms in this family , the ensemble Kalman ﬁlter and the maximum lik eliho o d ensemble ﬁlter, whic h are review ed next. 2.2. The ensemble Kalman ﬁlter Kalman ﬁlters (KF) [ 18 , 19 ] are sequen tial data ass imilation metho dologies, where mea- suremen ts are incorp orated at the time momen t when they b ecome av ailable. Sequen tial data assimilation algorithms proceed in tw o steps, namely , for e c ast and analysis . In the forecast step, the state of the system is propagated forw ard b y the model equations ( 1 ) to the next time p oin t where observ ations are av ailable, pro ducing a forecast of the state of the system, and a forecast error co v ariance matrix is presen ted to quan tify the uncertain ty of the forecast. The ensemble Kalman ﬁlter (EnKF) [ 6 , 10 , 11 , 16 ] tak es a Monte-Carlo approac h to represen ting the uncertain t y . An ensem ble of n ens states ( x a k − 1 ( e ), e = 1 , . . . , n ens ) is used to sample the analysis probabilit y distribution at time t k − 1 . Eac h mem b er of the ensem ble is propagated to t k using the nonlinear mo del ( 1 ) to obtain the ”forecast” ensem ble x f k ( e ) = M t k − 1 → t k ( x a k − 1 ( e )) + η k ( e ) , e = 1 , . . . , n ens . (3a) T o sim ulate the fact that the mo del is an imp erfect represen tation of realit y mo del errors are added. They are t ypically considered Gaussian random v ariables, η k ∈ N (0 , Q k ). The ensem ble mean and cov ariance approximate the background estimate and the background error co v ariance of the state at the next time p oin t t k : x f k = 1 n ens n ens X e =1 x f k ( e ) , (3b) X f k = [ x f k (1) − x f k , . . . , x f k ( n ens ) − x f k ] , (3c) B k =  1 n ens − 1  X f k  X f k  T   ◦ ρ. (3d) T o reduce sampling error due to the small ensemble size, lo calization [ 15 , 17 , 36 ] is p erformed b y taking the p oin t-wise pro duct of the ensem ble cov ariance and a decorrelation matrix ρ . Eac h mem b er of the forecast (ensemble of forecast states { x f k ( e ) } e =1 ,..., n ens ) is analyzed separately using the Kalman ﬁlter form ulas [ 6 , 10 ] x a k ( e ) = x f k ( e ) + K k  [ y k + ζ k ( e )] − H k ( x f k ( e ))  , (4a) K k = B k H T k  H k B k H T k + R k  − 1 . (4b) 5 The sto c hastic (“p erturb ed observ ations” ) version [ 6 ] of the ensem ble Kalman ﬁlter adds a diﬀeren t realization of the observ ation noise ζ k ∈ N (0 , R k ) to eac h individual assimilation. The Kalman gain matrix K k mak es use of the linearized observ ation op erator H k = H 0 k ( x f k ). The same Kalman gain is used for all ensem ble members. Square root versions (deterministic formulations) of EnKF [ 35 ] a void adding random noise to observ ations, and th us a void additional sampling errors. They also av oid the explicit construction of the full co v ariance matrices and work b y up dating only a matrix of state deviations from the mean. A detailed discussion of EnKF and v ariants can b e found in [ 12 ]. The main shortcomings of the ensem ble Kalman ﬁlter are the Gaussianity assumption on which the Kalman up dates are based. The ﬁlter is optimal only when the observ ation op erators are linear, and b oth the forecast and the observ ation errors are Gaussian. 2.3. The maximum likeliho o d ensemble ﬁlter The maxim um likelihoo d ensem ble ﬁlter (MLEF) [ 37 ] seeks to alleviate the limitations of the Gaussian assumptions b y computing the maxim um lik eliho o d estimate of the state in the ensem ble space. Sp eciﬁcally , it maximizes the p osterior probability densit y , or equiv alently , minimizes the follo wing nonlinear ob jective function ov er the ensemble subspace [ 22 , 37 ]: x opt k = arg min x J ( x ) , (5a) J ( x ) = 1 2 ( x − x b k ) T B − 1 k ( x − x b k ) (5b) + 1 2 ( y k − H k ( x )) T R − 1 k ( y k − H k ( x )) , and then updates the analysis error co v ariance matrix based on the fact that it is approxi- mately equal to the in verse of the Hessian matrix at the minimum [ 13 ]. The MLEF algorithms op erates sequentially by applying a forecast step and an analysis step. Let x opt k − 1 b e the optimal solution at the previous time p oin t t k − 1 , and let A 1 / 2 k − 1 = [ a k − 1 (1) , a k − 1 (2) , , . . . , a k − 1 ( n ens )] , (6) b e the matrix of scaled perturbations corresp onding to the analysis ensem ble at t k − 1 , suc h that the analysis co v ariance matrix is A k − 1 = A 1 / 2 k − 1 A T / 2 k − 1 . The forecast step pro vides the bac kground state x b k and a square ro ot of the background co v ariance matrix at the curren t time p oint t k as follo ws: x b k = M t k − 1 → t k ( x opt k − 1 ) , (7a) b k ( e ) = M t k − 1 → t k ( x opt k − 1 + a k − 1 ( e )) (7b) − M t k − 1 → t k ( x opt k − 1 ); e = 1 , . . . , n ens , B 1 / 2 k = [ b k (1) , b k (2) , , . . . , b k ( n ens )] . (7c) T o sp eed up the optimization problem ( 5a ) Hessian preconditioning is carried out through the c hange of v ariables x k ( ξ ) = x b k + B 1 / 2 k  I + C (0)  − T 2 ξ , (8) 6 where ξ is a vector of control v ariables in the ensemble space and C ( ζ ) = Z ( ζ ) T Z ( ζ ) , (9a) Z ( ζ ) = [ z ( ζ , 1) , z ( ζ , 2) , . . . , z ( ζ , n ens )] , (9b) z ( ζ , e ) = R − 1 2 k  y k − H k  x k ( ζ )  (9c) − R − 1 2 k  y k − H k  x k ( ζ ) + b k ( e )  . The matrix C (0) in ( 8 ) is obtained b y using x k (0) ≡ x b k in form ula ( 9d ). After replacing ( 8 ) in ( 5b ) the optimal solution is found b y solving the following mini- mization problem in the ensem ble subspace: ξ opt = arg min J ( ξ ) , (10) J ( ξ ) = 1 2 ξ T  I + C (0)  − 1 ξ (11) + 1 2 ( y k − H k ( x k ( ξ ))) T R − 1 k ( y k − H k ( x k ( ξ ))) . The gradien t reads: ∇ ξ J ( ξ ) =  I + C (0)  − 1 ξ (12) −  I + C (0)  − 1 / 2 Z ( ξ ) T R − 1 2 k ( y k − H k ( x k ( ξ ))) . The optimal solution in the mo del subspace is giv en b y: x opt k = x b k + B 1 / 2 k  I + C (0)  − T 2 ξ opt . (13) The matrix of scaled p erturbations represen ting the analysis is up dated as: A 1 / 2 k = B 1 / 2 k  I + C ( ξ opt )  − T 2 . (14a) An imp ortant adv antage of the algorithm is that the observ ation operator is not lin- earized. Consequently MLEF can w ork eﬃciently with non-linear observ ation op erators (without the requirement of diﬀeren tiability and without using ﬁnite-diﬀerence appro xima- tions of the Jacobian of the observ ation op erators) [ 37 ].) The cost function ( 5b ) to minimize implicitly assumes that the p osterior distribution is Gaussian. The metho d is unlik ely to giv e go o d results when the p osterior distributions are multimodal. 3. Hybrid Mark ov Chain Mon te Carlo Mark o v Chain Monte Carlo (MCMC) algorithms [ 26 ], in tro duced b y Metrop olis et. al [ 25 ], can sample from distributions with complex probabilit y densities π ( x ). They gener- ate a Marko v chain { x ( i ) } i ≥ 0 for whic h π ( x ) is the inv ariant (stationary) distribution, given that π ( x ) is known up to a multiplicativ e constant [ 26 ]. MCMC metho ds w ork b y generat- ing a random walk using a prop osal PDF and an “acceptance/rejection” criterion to decide 7 whether prop osed samples should b e accepted as part of the Mark o v c hain or should just b e rejected. These algorithms are generally p o werful, but may take a long time to explore the whole state space or ev en to con verge [ 34 ]. This section starts with a review of the Hybrid MCMC sampling (HMCMC) then presen ts the sampling ﬁlter algorithm for data assimilation. Hybrid Monte Carlo (HMC) metho ds, also kno wn an Hamiltonian Monte Carlo, orig- inated in the physics literature [ 9 ]. They attempt to handle the dra wbac ks of MCMC algorithms by incorp orating an auxiliary v ariable such as to reduce the correlation b et w een successiv e samples, to explore the entire space in very few steps, and to ensure high proba- bilit y of acceptance for prop osed samples in high dimensions [ 31 ]. 3.1. Hamiltonian dynamics Hamiltonian dynamical systems op erate in a phase space of p oin ts ( p , x ) ∈ R 2 n v ar , where the individual v ariables are the p osition x ∈ R n v ar and the momen tum p ∈ R n v ar . The total energy of the system is describ ed b y the Hamiltonian function H ( p , x ). The dynamics of the system in time is describ ed b y the follo wing ordinary diﬀerential equations: d x dt = ∇ p H , d p dt = −∇ x H . (15) The time ev olution of the system ( 15 ) in state space is describ ed b y the ﬂow [ 27 , 32 ] Φ T : R 2 n v ar → R 2 n v ar , Φ T  p (0) , x (0)  =  p ( T ) , x ( T )  , (16) whic h maps the initial state of the system ( p (0) , x (0)) to ( p ( T ) , x ( T )), the state of the system at time T . In practical computations the analytic ﬂo w Φ T is replaced by a numerical solution using a time reversible and symplectic n umerical integration method [ 32 , 31 ]. In this pap er we use ﬁv e diﬀerent high order symplectic integrators based on Strang’s splitting formula [ 32 ]: V erlet (St¨ ormer, Leapfrog) algorithm ( 42 ) [ 32 , 31 ], higher order integrators namely , t wo-stage ( 43 ), three-stage ( 44 ), and four-stage ( 45 ) p osition splitting in tegrators from [ 5 ], and the Hilb ert space in tegrator ( 46 ) from [ 3 ]. The methods are summarized in A . T o approximate Φ T the integrator at hand takes m steps of size h = T /m . With a sligh t abuse of notation w e will also denote by Φ T the ﬂo w of the numerical solution. 3.2. HMCMC sampling algorithm In order to draw samples { x ( e ) } e ≥ 0 from a giv en probabilit y distribution π ( x ) HMC mak es the following analogy with a Hamiltonian mec hanical system ( 15 ). The state x is view ed as a “p osition v ariable”, and an auxiliary “momen tum v ariable” p is included. The Hamiltonian function of the system is: H ( p , x ) = 1 2 p T M − 1 p − log( π ( x )) = 1 2 p T M − 1 p + J ( x ) . (17) The negative logarithm of the target probabilit y density J ( x ) = − log( π ( x )) is view ed as the potential energy of the system. The kinetic energy of the system is giv en b y the auxiliary 8 momen tum v ariable p . The constan t p ositiv e deﬁnite symmetric “mass matrix” M is y et to b e deﬁned [ 32 ]. Based on the Hamiltonian equations ( 15 ) the dynamics of the system is giv en b y d x dt = M − 1 p , d p dt = −∇ x J ( x ) . (18) The canonical probability distribution of the state of the system ( p , x ) in the phase space R 2 n v ar is, up to a constan t, equal to exp ( − H ( p , x )) = exp  − 1 2 p T M − 1 p − J ( x )  (19) = exp  − 1 2 p T M − 1 p  · π ( x ) . The pro duct form of this join t probabilit y distribution shows that the t w o v ariables p , x are indep enden t [ 31 ]. The distribution of the momen tum v ariable is Gaussian, p ∼ N (0 , M ), while the distribution of the p osition v ariable is the target probability density , x ∼ π [ 31 ]. The HMC sampling algorithm builds a Mark ov c hain starting from an initial state x 0 = x (0). Algorithm 1 summarizes the transition from the curren t Mark o v c hain state x k to a new state x k +1 [ 31 ]. Practical issues are related to the choice of the numerical integrator, the time step, and the c hoice of the function J ( x ) that represents the PDF we wish to sample from. The construction of the mass matrix M do es not impact the ﬁnal distribution, but do es aﬀect the computational p erformance of the algorithm [ 14 ]. The mass matrix M is symmetric and p ositive deﬁnite and is a parameter that is tuned b y the user. It can be for example, a constant multiple of the identit y [ 27 ], or a diagonal matrix whose entries are the bac kground error v ariances [ 3 , 21 ]. W e found that the latter approac h is more eﬃcien t for the curren t application and used it in all numerical exp erimen ts rep orted here. 9 Algorithm 1 HMCMC Sampling [ 31 ]. 1: Dra w a random v ector p k ∼ N (0 , M ). 2: Use a symplectic numerical integrator (from A ) to adv ance the curren t state ( p k , x k ) b y a time incremen t T to obtain a pr op osal state ( p ∗ , x ∗ ): ( p ∗ , x ∗ ) = Φ T  ( p k , x k )  . (20) 3: Ev aluate the loss of energy based on the Hamiltonian function. F or the standard V erlet ( 42 ), t wo-stage ( 43 ), three-stage ( 44 ), and four-stage ( 45 ) in tegrators [ 5 , 31 ] the loss of energy is computed as: ∆ H = H ( p ∗ , x ∗ ) − H ( p k , x k ) . (21) F or the Hilb ert space integrator ( 46 ) [ 3 ] the loss of energy is computed as: ∆ H = φ ( x ∗ ) − φ ( x k ) (22) + h 2 8  | M − 1 2 ( −∇ φ ( x k )) | 2 − | M − 1 2 ( −∇ φ ( x ∗ )) | 2  + h m − 1 X i =1  p T k ( −∇ φ ( x k ))  + h 2  p T k ( −∇ φ ( x k )) + ( p ∗ ) T ( −∇ φ ( x ∗ ))  , where φ ( x ) = − log ( π ( x )) and h = T /m is the integration time step [ 31 ]. 4: Calculate the probabilit y: a ( k ) = 1 ∧ e − ∆ H . (23) 5: Discard b oth p ∗ , p k . 6: (Acceptance/Rejection) Draw a uniform random v ariable u ( k ) ∼ U (0 , 1): i- If a ( k ) > u ( k ) accept the prop osal as the next sample: x k +1 := x ∗ ; ii- If a ( k ) ≤ u ( k ) reject the prop osal and con tin ue with the current state: x k +1 := x k . 7: Rep eat steps 1 to 6 un til suﬃcien tly many distinct samples are drawn. 4. The Sampling Filter for Data Assimilation The goal of this ﬁlter is to replace the analysis step in the traditional EnKF with a resampling pro cedure that draws represen tative ensemble mem b ers from the p osterior dis- tribution π ( x ) = P a ( x ). Even if the p osterior may in general b e non-Gaussian we assume, as most of the curren t ensemble-based data assimilation algorithms, that the posterior has 10 the form: π ( x ) = P a ( x ) ∝ exp  −J ( x )  , (24) J ( x ) = 1 2  x − x b  T B − 1  x − x b  (25) + 1 2  y − H ( x )  T R − 1  y − H ( x )  . where x b is the background state (forecast), y is the observ ation v ector, and H is the observ ation op erator that is generally non-linear. F or sampling at time t k the corresp onding J ( x ) is: J ( x ) = − log  P a ( x )  (26) = 1 2  x − x b k  T B − 1 k  x − x b k  (27) + 1 2  y k − H k ( x )  T R − 1 k  y k − H k ( x )  , and its gradien t has the form ∇ x J ( x ) = B − 1 k ( x − x b k ) + H T k R − 1 k  y k − H k ( x )  , (28) where H k = H 0 k ( x ) is the linearized observ ation op erator. Algorithm ( 1 ) is used to generate n ens ensem ble mem b ers dra wn from the p osterior distribution { x a k ( e ) ∼ P a ( x ) } e =1 , 2 ,..., n ens . The mean of this ensemble is an estimate of the analysis state, and the ensem ble cov ariance estimates the analysis error co v ariance matrix. Note that the prop osed sampling ﬁlter is not restricted to a sp eciﬁc form of the p osterior PDF, and the Gaussian assumption ( 26 ) can in principle b e remov ed. The remaining issue is to represen t non-Gaussian probability density functions and their logarithm. In the next section w e describ e the prop osed sampling ﬁlter as an alternativ e to the EnKF. diagonal The sampling ﬁlter is describ ed in Algorithm 2 . Like most of the ensemble-based se- quen tial data assimilation algorithms the sampling ﬁlter consists of t wo stages, namely , the for e c ast step and the analysis step. Start with an ensemble { x a k − 1 ( e ) } e =1 ,..., n ens describing the analysis PDF at time t k − 1 . In the forecast step eac h ensemble mem b er is propagated by the full mo del to the next time t k − 1 where observ ations are av ailable, resulting in the forecast ensemble. In the analysis step the HMCMC algorithm is simply used to sample from the posterior PDF of the state, pro viding the new analysis ensemble { x a k ( e ) } e =1 ,..., n ens . 11 Algorithm 2 Sampling Filter 1: F orecast step: giv en an analysis ensemble { x a k − 1 ( e ) } e =1 , 2 ,..., n ens at time t k − 1 ; generate the forecast ensem ble b y via the mo del M : x b k ( e ) = M t k − 1 → t k  x a k − 1 ( e )  , e = 1 , 2 , . . . , n ens . (29) 2: Analysis step: giv en the observ ation vector y k at time p oin t t k , follo w the steps: i- Set the initial state x 0 of the Mark ov Chain to b e to the best estimate av ailable, e.g., the mean of the forecast ensemble. One can use the EnKF analysis if the cost is acceptable, and this choice is exp ected to result in a faster conv ergence of the c hain to the stationary distribution. ii- Calculate the ensemble-based forecast error co v ariance matrix B k (and p ossibly balance it by a ﬁxed (or frequently up dated) cov ariance matrix B 0 ), and apply lo calization as in equation ( 3d ). It is imp ortan t to emphasize that building the full background error co v ariance matrix is not necessary for the curren t algorithm to w ork. iii- Cho ose a p ositiv e deﬁnite diagonal mass matrix M . One c hoice that fav ors the p erformance of the sampling algorithm is the diagonal of the matrix B − 1 k [ 27 ] whic h scales the comp onents of the state vector v ary . Ideally , M should b e set to the diagonal of the in verse p osterior co v ariance matrix. iv- Apply Algorithm 1 with initial state x 0 and generate n ens ensem ble mem b ers. In practice one starts accepting samples after a w arm-up phase (of, say , 30 steps), to guaran tee that selected mem b ers explore the entire state space. v- Use the generated samples { x a k ( e ) } e =1 , 2 ,..., n ens as an analysis ensem ble and calculate the b est estimate of the state (e.g. the mean), and the analysis error cov ariance matrix. 3: Increase time k := k + 1 and rep eat steps 1 and 2. As stated in step ii of Algorithm 2 , the explicit represen tation of the matrix B k is not necessary - one only needs to apply its inv erse to a vector in ( 26 ), ( 28 ). T ypically B k is formed as a linear a combination betw een a ﬁxed matrix B 0 and the ensem ble cov ariance. The calculation requires to ev aluate the pro ducts u = B − 1 k ( x − x b k ) (30) = γ B 0 + 1 − γ n ens − 1 n ens X e =1 ∆ x ( e ) (∆ x ( e )) T ! − 1 ( x − x b k ) , where ∆ x ( e ) is the deviation of the ensem ble mem b er x ( e ) from the mean of the ensemble. The linear system γ B 0 + 1 − γ n ens − 1 n ens X e =1 ∆ x ( e ) (∆ x ( e )) T ! · u = x − x b k , (31) 12 can b e solv ed without ha ving to build the full matrix B k as discussed in [ 29 ]. In our n umerical exp eriments w e build ﬂo w-dep enden t bac kground error co v ariance ma- trices B k at eac h time step. W e set M to be equal to the diagonal of B k in case of Lorenz-96 mo del following ([ 3 , 21 ]. T aking M equal to the diagonal of B − 1 k lead to similar results for the Lorenz-96 model. F or the shallo w-water mo del on the sphere we set M to b e equal to the diagonal of B − 1 k . 5. Numerical Results 5.1. The L or enz-96 mo del Numerical tests are primarily p erformed using the 40-v ariables Lorenz-96 mo del [ 23 ] whic h is describ ed b y the equations: dx i dt = x i − 1 ( x i +1 − x i − 2 ) − x i + F , (32) where x = ( x 1 , x 2 , . . . , x 40 ) T ∈ R 40 is the state v ector. The indices w ork in a circular fashion, e.g., x 0 ≡ x 40 . The forcing parameter is set to F = 8 in our exp erimen t. These settings mak e the system c haotic [ 24 ]. The initial condition is obtained b y integrating a v ector of equidistan t comp onen ts ranging from − 2 to 2 for 10 time units before the b eginning of the exp erimen t time interv al. The sim ulation time in terv al is [0 , 10] units with observ ations a v ailable at time p oin ts { t k = 0 . 1 × k } k =1 , 2 ,..., 101 . T o study the b eha vior of the sampling algorithm with small ensemble, the n umber of ensemble mem b ers is c hosen to b e 30. All observ ations are synthetic, created by applying the observ ation op erator to the reference tra jectory (b y applying the corresp onding observ ation op erator) and adding Gaussian noise with a standard deviation equal to 5% of the a v erage magnitude of the corresp onding ob- serv ation along the reference tra jectory . The bac kground error is Gaussian with a diagonal co v ariance matrix B 0 ; the standard deviation of each comp onen t is 8% of the av erage mag- nitude of the initial condition of the system. 5.2. Observations and observation op er ators W e choose six diﬀerent observ ation op erators of diﬀerent complexities and v arying lev els of non-linearit y to test the performance of the sampling ﬁlter. All the six op erator w ere used with the Lorenz-96 mo del. Both quadratic and cubic observ ation operators used here were emplo y ed by Zupanski [ 37 , 38 ] in the simple case of one dimensional state space. Syn thetic observ ations are obtained by applying them to a reference tra jectory and adding Gaussian random noise with a standard deviation of 5% of the magnitude of the reference observ ation v alues. Line ar observation op er ator. The ﬁrst observ ation op erator is a linear op erator that selects a sp eciﬁc subset of the comp onen ts of the state v ector. This op erator mak es J ( x ) diﬀeren- tiable. In our experiments w e observe eac h third component of the state, starting with the ﬁrst comp onen t H ( x ) = Hx = ( x 1 , x 4 , x 7 , . . . , x 37 , x 40 ) T ∈ R 14 . (33) 13 Quadr atic observation op er ator. This is a non-linear but diﬀeren tiable observ ation op erator that squares selected comp onen ts ( 33 ) of the state. In our exp erimen ts we use: H ( x ) = ( x 2 1 , x 2 4 , x 2 7 , . . . , x 2 37 , x 2 40 ) T ∈ R 14 . (34) Cubic observation op er ator. This is another non-linear but diﬀeren tiable observ ation oper- ator that squares selected comp onen ts ( 33 ) of the state. In our exp erimen ts we use: H ( x ) = ( x 3 1 , x 3 4 , x 3 7 , . . . , x 3 37 , x 3 40 ) T ∈ R 14 . (35) Magnitude observation op er ator. This non-diﬀeren tiable observ ation op erator returns the absolute v alues of selected comp onen ts ( 33 ). The observ ation vector reads: H ( x ) = ( | x 1 | , | x 4 | , | x 7 | , . . . , | x 37 | , | x 40 | ) T ∈ R 14 . (36) Quadr atic observation op er ator with a thr eshold. This observ ation op erator is similar to the simple v ersion used b y Zupanski et al in [ 38 ]. The observ ation vector is: H ( x ) = ( x 0 1 , x 0 4 , x 0 7 , . . . , x 0 37 , x 0 40 ) T ∈ R 14 , (37) where x 0 i =  x 2 i : x i ≥ 0 . 5 − x 2 i : x i < 0 . 5 , This op erator is non-linear and discon tin uous. Exp onential observation op er ator. This is a highly nonlinear, diﬀerentiable observ ation op- erator: H ( x ) = ( e r · x 1 , e r · x 4 , e r · x 7 , . . . , e r · x 37 , e r · x 40 ) T ∈ R 14 , (38) where r ∈ R is a scaling factor that controls the degree of nonlinearity . 5.3. Exp erimental setting W e p erform t w o sampling ﬁlter data assimilation exp erimen ts with eac h observ ation op erator describ ed in Section 5.2 . Both share the same mo del parameters but use diﬀerent step sizes of the symplectic in tegration during the HMC sampling. This is found to hav e a great impact on the p erformance of the sampling ﬁlter. In the ﬁrst exp erimen t a time T = 0 . 1 with h = 0 . 01 and m = 10 is used for all inte- grators tested. This choice guarantees that the standard p osition V erlet integrator yields satisfactory results with the linear observ ation operator, but the performance on nonlinear observ ation op erators remains to b e c hec ked. The second exp eriment analyzes the p erfor- mance of the sampling ﬁlter when all time integrators take roughly the same computational cost. The parameters, m , h , are tuned b y trial and error such as to make the V erlet in- tegrator successful, if p ossible, with the nonlinear observ ation op erators. The other time in tegrators use the same total time step T as the V erlet in tegrator; the v alues of m and h are c hosen for eac h metho d suc h that the n umber of gradien t calculations don e b y all in tegrators 14 is the same. In general, how ever, the time-stepping parameters of eac h symplectic in tegrator should b e set individually to get the b est p erformance of the sampling ﬁlter. Eac h n umerical exp erimen t p erforms 100 realizations of the sampling ﬁlter. Eac h realiza- tion uses the same settings but the sequence of random num b er generated b y the sampling ﬁlter, for b oth the p otential v ariable and the acceptance/rejection rule, was diﬀerent. The ro ot mean squared error (RMSE) metric is used to compare the analyses against the reference solution at observ ation time p oin ts: RMSE = v u u t 1 n v ar n v ar X i =1 ( x i − x true i ) 2 , where x true is the reference state of the system. The RMSE is calculated at all assimilation time p oin ts along the tra jectory o ver the time span of the exp erimen t. T o guaran tee that the Marko v c hain reac hes the stationary distribution b efore starting the sampling process a set of 200 steps are perform as burn-in stage. W e noticed that the c hain alwa ys conv erges in a small num b er (10 − 20) burn-in steps. Stationarity tests will b e giv en sp ecial atten tion in our future work. After the burn-in stage an ensem ble mem b er is selected after eac h 30 generated states; this choice decreases correlation b et w een generated ensemble members since the c hain is not memoryless. The n umber of ensem bles that are not retained will be referred to as the n um b er of in ter-c hain steps. I n our exp erimen ts the acceptance probability is high (usually o v er 0 . 9) with this sampling strategy . The n umber of in ter-chain steps is a parameter that can b e tuned b y the user to con trol the p erformance of the sampling ﬁlter. Stabilit y requiremen ts imp ose tigh t upp er bounds on the step size h of the V erlet in te- grator. The step size should b e decreased with the increasing dimension of the system in order to main tain O (1) acceptance probability [ 2 ]. On the other hand large steps of the symplectic integrator are needed in order to explore the space eﬃciently . There is no precise rule av ailable to select the optimal step size v alues [ 31 ] and consequen tly h should b e tuned for each problem. The higher-order in tegrators ( 43 ), ( 44 ), ( 45 ) a re expected to be more stable than V erlet for larger time steps [ 5 , 27 ]. T o guarantee ergo dicity of the Mark o v chain, which is a prop ert y required for the chain to conv erge to its inv ariant distribution, we follo w [ 5 , 27 ] and change the step length at the b eginning of each Marko v step (once at the b eginning of the Hamiltonian tra jectory) to h = (1 + r ) h ref where h ref is a reference step size and r ∼ U ( − 0 . 2 , 0 . 2) is a uniformly distributed random v ariable. Randomizing the step size of the symplectic integrator, in addition to other beneﬁts, ensures that the results obtained are not en trusted with sp eciﬁc c hoice of the step size [ 27 ]. 5.4. Line ar observation op er ator exp eriments Figure 1 sho ws the analysis results of diﬀerent ﬁlters when the system uses linear obser- v ation op erators ( 33 ). The accuracy of the analyses pro vided b y diﬀeren t ﬁlters is plotted at diﬀerent time moments. Results are rep orted for 100 instances of the sampling ﬁlter. 15 The red line represents the median RMSE v alues across all instances and the cen tral blue b o x represen ts the v ariance. The tw o vertical lines (whisk ers) extend up to 1 . 5 times the heigh t of the cen tral b o x. The v alues exceeding the length of the whisk ers are considered outliers (extremes) and are plotted as red crosses. All symplectic in tegrators sho w outliers (the red crosses) with the exception of the Hilb ert space integrator. The RMSE errors of the sampling in tegrator are larger than those of EnKF, how ever they remain small o verall. The analysis follo ws closely the reference tra jectory as seen in Figure 22 . (a) P osition V erlet integrator ( 42 ) (b) Two-stage integrator ( 43 ) (c) Three-stage in tegrator ( 44 ) (d) F our-stage integrator ( 45 ) (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 1: Data assimilation results with the linear observ ation operator ( 33 ). The symplectic integrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 inter-c hain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as box plots. The red line represen ts the median RMSE v alues across all instances and the cen tral blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. Figure 2 shows results with the time step parameters of symplectic integrators tuned 16 to provide equalized w ork. The n um b er of steps m for V erlet is increased compared to the tests in Figure 1 for t wo reasons: to test the capabilities of p osition V erlet with diﬀeren t step sizes, and to allo w for more accuracy using V erlet in tegrators. The sampling ﬁlters p erform w ell and are as accurate as EnKF with any choice of integrator, except for the Hilb ert space one whic h yields larger RMSE errors. (a) Position V erlet integrator ( 42 ); h = 0 . 01 , m = 24 (b) Two-stage in tegrator ( 43 ); h = 0 . 02 , m = 12 (c) Three-stage in tegrator ( 44 ); h = 0 . 03 , m = 8 (d) F our-stage integrator ( 45 ); h = 0 . 04 , m = 6 (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 01 , m = 24 Figure 2: Data assimilation results with the linear observ ation operator ( 33 ). The symplectic integrator used is indicated under each panel. The time step for all in tegrators is T = 0 . 24 (units), and h and m are c hosen suc h as to equalize the computational eﬀort. The n umber of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as b ox plots. The red line represents the median RMSE v alues across all instances and the central blue b ox represents the v ariance. The t wo v ertical lines (whisk ers) extend up to 1 . 5 times the height of the cen tral b o x. The v alues exceeding the length of the whisk ers are considered outliers (extremes) and are plotted as red crosses. 17 5.5. Quadr atic observation op er ator exp eriments Figure 3 sho ws the results with quadratic observ ation op erator ( 34 ). All symplectic in tegrators use the parameters h = 0 . 01 and m = 10. The sampling ﬁlter with V erlet in tegrator fails to con v erge and pro duce represen tative samples from the analysis PDF. When high-order integrators are used the sampling ﬁlter gives satisfactory analysis RMSE, comparable to that obtained by MLEF, except for occasional failures represented in the plots as outliers (red crosses). Section 5.11 will discuss strategies to handle p ossible failures and av oid these outliers. The ﬁlter with the Hilb ert space integrator has a larger RMSE error than EnKF, ho w ever it do es not suﬀer from outliers as muc h as the other in tegrators. 18 (a) P osition V erlet integrator ( 42 ) (b) Two-stage integrator ( 43 ) (c) Three-stage in tegrator ( 44 ) (d) F our-stage integrator ( 45 ) (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 3: Data assimilation results with the quadratic observ ation op erator ( 34 ). The symplectic integrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 inter-c hain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as box plots. The red line represen ts the median RMSE v alues across all instances and the cen tral blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. Figure 4 shows results with the time parameters tuned suc h as to obtain the b est possible results with the V erlet scheme; the step sizes of other in tegrators are c hosen suc h that their w ork p er step is equal to V erlet’s work p er step. The V erlet integrator results in high uncertain t y in the RMSE whic h mak es the divergence of the ﬁlter v ery likely . High-order in tegrators contin ue to give go o d results but the num b er outliers seem to increase. The in tegrator deﬁned on Hilb ert space fails completely and yields large RMSE in many cases. W e conclude that the step sizes should b e tuned indep enden tly for eac h integrator. The exp erimen ts indicate that each of the in tegrators, except p erhaps V erlet, can b e tuned to 19 giv e v ery satisfactory ﬁltering results. (a) Position V erlet integrator ( 42 ); h = 0 . 01 , m = 40 (b) Two-stage in tegrator ( 43 ); h = 0 . 02 , m = 20 (c) Three-stage integrator ( 44 ); h = 0 . 03 , m = 13 (d) F our-stage in tegrator ( 45 ); h = 0 . 04 , m = 10 (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 01 , m = 40 Figure 4: Data assimilation results with the quadratic observ ation op erator ( 34 ). The symplectic integrator used is indicated under eac h panel. The time step for all in tegrators is T = mh , with h, m , as indicated under each panel. The num b er of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as box plots. The red line represents the median RMSE v alues across all instances and the central blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. 5.6. Cubic observation op er ator exp eriments Figure 5 shows the results with cubic observ ation op erator ( 35 ). Both EnKF and MLEF fail to con ver ge due to high non-linearity of the observ ation op erator. The MLEF failure w as unexpected and ma y be due to its sensitivit y to the uncertaint y lev els of either the 20 bac kground, or the observ ations or b oth. The level of nonlinearit y of the observ ation op er- ator has a ma jor impact on the success of the MLEF ﬁlter as sho wn in 5.10 . The sampling ﬁlter with p osition V erlet integrator fails to con v erge. The results are better for high-stage in tegrators, and the four-stage integrator provides satisfactory results that are similar to those obtained with linear observ ation op erators. (a) P osition V erlet integrator ( 42 ) 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast EnKF MLEF (b) Tw o-stage integrator ( 43 ) (c) Three-stage in tegrator ( 44 ) (d) F our-stage integrator ( 45 ) (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 5: Data assimilation results with the cubic observ ation operator ( 35 ). The symplectic in tegrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 inter-c hain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as box plots. The red line represen ts the median RMSE v alues across all instances and the cen tral blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. Figure 6 shows results with tuned parameters suc h that the work is equal for all in te- grators. Position V erlet requires more w ork and ﬁner step sizes to pro vide conv ergence of the sampling ﬁlter, and ev en in this case there are man y outliers that sho w divergence, as 21 seen in Figure 6(a) . The high-order integrators give go o d results, but reducing the step size increases their computational costs. As sho wn in Figure 6(e) , the Hilb ert in tegrator leads to large RMSE with this setting of step size. Again, it is advisable to tune the step size of this in tegrator indep enden tly . (a) Position V erlet integrator ( 42 ); h = 0 . 001 , m = 40 (b) Two-stage integrator ( 43 ); h = 0 . 002 , m = 20 (c) Three-stage in tegrator ( 44 ); h = 0 . 003 , m = 13 (d) F our-stage integrator ( 45 ); h = 0 . 004 , m = 10 (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 001 , m = 40 Figure 6: Data assimilation results with the cubic observ ation operator ( 35 ). The symplectic in tegrator used is indicated under eac h panel. The time step for all in tegrators is T = mh , with h, m , as indicated under each panel. The num b er of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as box plots. The red line represents the median RMSE v alues across all instances and the central blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. 22 5.7. Absolute value observation op er ator exp eriments Figure 7 sho ws the results with absolute observ ation op erator ( 36 ). The Jacobian of this observ ation op erator is taken as the sign of the measured comp onen ts of the state vector. Similar to the case of linear observ ation operator, MLEF conv erges in the beginning, since the observ ation op erator is w eekly non-linear, but div erges later in the exp eriment, mostly due to large observ ation errors and low observ ation frequency . The sampling ﬁlter using Hilb ert in tegrator shows impro vemen t ov er the forecast, but its analysis is less accurate than MLEF or EnKF analyses (when they conv erge). V erlet and the high-order integrators b eha ve almost iden tically . The distribution of outliers is similar to that for quadratic observ ation op erator. W e will discuss how to deal with the o ccasional ﬁlter divergence in Section 5.11 . 23 (a) P osition V erlet integrator ( 42 ) 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 Time RMSE Forecast EnKF MLEF (b) Tw o-stage integrator ( 43 ) (c) Three-stage in tegrator ( 44 ) (d) F our-stage integrator ( 45 ) (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 7: Data assimilation results with the magnitude observ ation operator ( 36 ). The symplectic integrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 inter-c hain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as box plots. The red line represen ts the median RMSE v alues across all instances and the cen tral blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. Figure 8 sho ws results with a larger step size and with equalized integrator work. The results with V erlet are similar to those rep orted in Figure 7 , ho wev er the results obtained using high-order in tegrators are worse than before. The use of Hilb ert space in tegrator results in large RMSE, ho wev er these errors are stable (do not increase) with time. 24 (a) Position V erlet integrator ( 42 ); h = 0 . 01 , m = 24 (b) Two-stage in tegrator ( 43 ); h = 0 . 02 , m = 12 (c) Three-stage in tegrator ( 44 ); h = 0 . 03 , m = 8 (d) F our-stage integrator ( 45 ); h = 0 . 04 , m = 6 (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 01 , m = 24 Figure 8: Data assimilation results with the magnitude observ ation operator ( 36 ). The symplectic integrator used is indicated under eac h panel. The time step for all in tegrators is T = mh , with h, m , as indicated under each panel. The num b er of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as box plots. The red line represents the median RMSE v alues across all instances and the central blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. 5.8. Quadr atic observation op er ator with thr eshold exp eriments Figure 9 shows the results with quadratic observ ation operator ( 37 ) with threshold a = 0 . 5. Even if MLEF was successfully tested with one dimensional mo dels with this version of observ ation op erator [ 38 ], it do es not p erform w ell with the Lorenz mo del. The sampling ﬁlter using V erlet in tegrator fails due to the high non-linearity of the observ ation op erator and/or the uncertain t y lev els. The high order integ rators show goo d results and the analysis 25 RMSE has the level obtained in case of linear observ ation op erator. W e can conclude that the ensemble pro duced b y the ﬁlter is represen tativ e to the p osterior PDF as b oth the mean and the cov ariance are incorporated in the analysis steps.The likelihoo d of outliers is small and decreases using higher-order integrators. The Hilb ert integrator gives reasonable results. (a) P osition V erlet integrator ( 42 ) (b) Two-stage integrator ( 43 ) (c) Three-stage in tegrator 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 Time RMSE Forecast MLEF EnKF (d) F our-stage integrator ( 45 ) (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 9: Data assimilation results with the quadratic observ ation op erator ( 37 ) with a threshold a = 0 . 5. The symplectic integrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 in ter-chain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as b o x plots. The red line represen ts the median RMSE v alues across all instances and the central blue b o x represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the cen tral b o x. The v alues exceeding the length of the whisk ers are considered outliers (extremes) and are plotted as red crosses. Figure 10 sho ws results with obtained with step sizes and equalized work.The ﬁlter with four stage integrator is sup erior with such level of non-linearity as it suﬀers the least from outliers and giv es a small RMSE. The V erlet in tegrator also gives satisfactory results. The 26 Hilb ert integrator results in large analysis RMSE but it performs robustly ev en if the analysis is v ery far from the true solution and large step sizes are selected. (a) Position V erlet integrator ( 42 ); h = 0 . 01 , m = 24 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 Time RMSE Forecast EnKF MLEF (b) Two-stage in tegrator ( 43 ); h = 0 . 02 , m = 12 (c) Three-stage in tegrator ( 44 ); h = 0 . 03 , m = 8 (d) F our-stage integrator ( 45 ); h = 0 . 04 , m = 6 (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 01 , m = 24 Figure 10: Data assimilation results with the quadratic observ ation op erator ( 37 ) with a threshold a = 0 . 5. The symplectic in tegrator used is indicated under each panel. The time step for all in tegrators is T = mh , with h, m , as indicated under each panel. The num b er of in ter-chain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as b o x plots. The red line represents the median RMSE v alues across all instances and the central blue box represen ts the v ariance. The tw o vertical lines (whisk ers) extend up to 1 . 5 times the height of the central b ox. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. The Jacobian of this observ ation op erator is appro ximated using ﬁnite diﬀerences. Al- ternativ es will b e considered in the future. 27 5.9. Exp onential observation op er ator (with factor r = 0 . 2 ) exp eriments Figure 11 sho ws the results with the exponential observ ation op erator ( 38 ) with factor r = 0 . 2. This observ ation operator is diﬀeren tiable, ho wev er small perturbations in the state migh t result in relativ ely large changes in the measured v ales. Under strongly non- linear conditions the sampling ﬁlter performs b etter than either MLEF and EnKF. The p erformance of the sampling ﬁlter in this exp erimen t is similar to its p erformance in case of linear observ ation op erators. 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE EnKF Forecast MLEF (a) P osition V erlet integrator ( 42 ) 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (b) Tw o-stage integrator ( 43 ) 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE EnKF Forecast MLEF (c) Three-stage in tegrator ( 44 ) 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (d) F our-stage integrator ( 45 ) 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (e) In tegrator deﬁned on Hilb ert space ( 46 ) Figure 11: Data assimilation results with the exp onential observ ation op erator ( 38 ) with a factor r = 0 . 2. The symplectic integrator used is indicated under each panel. The time step for all integrators is T = 0 . 1 with h = 0 . 01, m = 10, and 30 in ter-chain steps. The RMSE for 100 instances of the sampling ﬁlter results are sho wn as b o x plots. The red line represen ts the median RMSE v alues across all instances and the central blue b o x represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the cen tral b o x. The v alues exceeding the length of the whisk ers are considered outliers (extremes) and are plotted as red crosses. 28 As shown in Figure 12(a) , further tuning of the step size for V erlet in tegrator do es not result in notable impro v ements o ver the results in Figure 11(a) . Figures 12(b) , 12(c) , and 12(d) sho w that the tw o-stage, three-stage, and four-stage in tegrators b eha ve similarly , and giv e slightly b etter results than those rep orted in Figures 11(b) , 11(c) , and 11(d) , resp ectiv ely . The inﬁnite dimensional in tegrator p erformance do es not c hange with the c hange in step size, as can b e seen in Figures 11(e) , and 12(e) . 29 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (a) Position V erlet integrator ( 42 ); h = 0 . 01 , m = 16 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (b) Tw o-stage integrator ( 43 ); h = 0 . 02 , m = 8 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (c) Three-stage in tegrator ( 44 ); h = 0 . 03 , m = 6 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (d) F our-stage integrator ( 45 ); h = 0 . 04 , m = 4 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 Time RMSE Forecast EnKF MLEF (e) Integrator deﬁned on Hilbert space ( 46 ); h = 0 . 01 , m = 16 Figure 12: Data assimilation results with the exp onen tial observ ation op erator ( 38 ). The symplectic in tegra- tor used is indicated under each panel. The time step for all integrators is T = mh , with h, m , as indicated under each panel. The num b er of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as box plots. The red line represents the median RMSE v alues across all instances and the central blue b ox represen ts the v ariance. The tw o vertical lines (whiskers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. 5.10. MLEF p erformanc e The standard v ersion of MLEF seems to b e sensitiv e to the lev el of uncertainties in observ ations and bac kground state, and to the degree of nonlinearit y of the observ ation op erator. Figures 15 , 17 , and 18 show the results of MLEF applied to the tests with cubic, quadratic with a threshold, and exponential (with a factor r = 0 . 2) observ ation op erators, 30 resp ectiv ely . In these tests all v ariables of the model are observ ed (unlik e observing only each third comp onen t of the state v ector as in the previous tests). Also, sev eral uncertain ty lev els are considered. The results indicate that the MLEF p erformance degrades considerably when the observ ations are sparser (when only each second or third v ariables are observed). Also, the p erformance degrades for higher uncertain ty lev els and for higher degrees of nonlinearit y of the observ ation op erators. 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 13: MLEF data assimilation results with the linear observ ation op erator ( 33 ). MLEF is applied with v arying observ ation frequencies, and with diﬀeren t noise levels for b oth observ ations and the bac kground state. The frequency of observ ations is indicated under each panel. The background error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 31 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 14: MLEF data assimilation results with the quadratic observ ation op erator ( 34 ). MLEF is applied with v arying observ ation frequencies, and with diﬀerent noise lev els for both observ ations and the bac kground state. The frequency of observ ations is indicated under each panel. The background error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 32 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 15: MLEF data assimilation results with the cubic observ ation op erator ( 35 ). MLEF is applied with v arying observ ation frequencies, and with diﬀeren t noise levels for b oth observ ations and the bac kground state. The frequency of observ ations is indicated under each panel. The background error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 33 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 16: MLEF data assimilation results with the magnitude observ ation op erator ( 36 ). MLEF is applied with v arying observ ation frequencies, and with diﬀerent noise lev els for both observ ations and the bac kground state. The frequency of observ ations is indicated under each panel. The background error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 34 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 17: MLEF data assimilation results with the quadratic observ ation operator ( 37 ) with a threshold a = 0 . 5. MLEF is applied with v arying observ ation frequencies, and with diﬀerent noise lev els for both observ ations and the bac kground state. The frequency of observ ations is indicated under eac h panel. The bac kground error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 35 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (a) All comp onen ts are observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (b) Eac h second comp onent is observed 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 14 Time RMSE Forecast MLEF: σ x0 =4%, σ obs =1% MLEF: σ x0 =6%, σ obs =3% MLEF: σ x0 =8%, σ obs =5% (c) Eac h third comp onent is observed Figure 18: MLEF data assimilation results with the exp onen tial observ ation operator ( 38 ). MLEF is applied with v arying observ ation frequencies, and with diﬀerent noise lev els for both observ ations and the bac kground state. The frequency of observ ations is indicated under each panel. The background error standard deviation is σ x 0 , and the observ ation error standard deviation is σ obs . 36 5.11. T uning the numb er of MC steps b etwe en suc c essive state sele ctions This section discusses the preven tion of outliers (ﬁlter div ergence) that can happ en for nonlinear observ ation op erators (e.g., in case of quadratic observ ation op erator in our ex- p erimen ts). This is done by tuning in tegration parameters. In addition to selecting the mass matrix M and the n umber of burn-in steps, there t wo more parameters to b e tuned. They are the step size of the symplectic integrator as discussed b efore, and the num b er of steps skipp ed b etw een selected states at stationarity (referred to as inter-c hain steps). T o study the eﬀect of tuning the last tw o parameters the quadratic observ ation op erator is re-tested with the high-order integrators and with the optimal step sizes suggested by Blanes [ 5 ]. Figure 19 shows the a verage and the standard deviation of analysis RMSE for 25 realizations of the sampling ﬁlter. V arious settings of the num b er of inter-c hain steps are used to study its eﬀect on the p erformance of the prop osed ﬁlter. T uning the step size of the symplectic in tegration results in a notable reduction in the av erage RMSE compared to the results obtained with the empirical settings and presented in Section 5.5 . Outliers are still presen t as inferred from Figures 19(b) , 19(d) , and 19(f ) . T uning the n umber of inter-c hain steps can in principle greatly enhance b oth the p er- formance of the ﬁlter and the reliabilit y of the results. Setting the num b er of inter-c hain steps to 30 is not optimal for the quadratic observ ation op erator, and b etter results can b e obtained with 40 steps, as seen in the results reported in Figure 20 . These results indicate that a careful tuning of both the step size and the n um b er of inter-steps in the c hain ma y o v ercome the problem of outliers and lead to the desired p erformance of the ﬁlter. 37 (a) RMSE mean: t wo-stage in tegrator; h = 2 / n v ar , m = n v ar (b) RMSE standard deviation: tw o-stage inte- grator; h = 2 / n v ar , m = n v ar (c) RMSE mean: three-stage integrator; h = 3 / n v ar , m = n v ar / 2 (d) RMSE standard deviation: three-stage inte- grator; h = 3 / n v ar , m = n v ar / 2 (e) RMSE mean: four-stage integrator; h = 4 / n v ar , m = n v ar / 2 (f ) RMSE standard deviation: t wo-stage inte- grator; h = 4 / n v ar , m = n v ar / 2 Figure 19: Data assimilation results with the quadratic observ ation op erator ( 34 ). The sampling ﬁlter is applied with sev eral settings of the n umber of in ter-chain steps. The av erage and standard deviation of RMS errors are b oth ev aluated, at each observ ation time p oin t, ov er the 25 realization of the sampling ﬁlter. The step size h and the num b er of steps m for the symplectic integrator are indicated under each panel, where n v ar = 40 is the dimension of the state v ector. 38 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 Time RMSE Forecast EnKF MLEF Figure 20: Data assimilation results with the quadratic observ ation operator ( 34 ). The symplectic integrator used is the three-stage in tegrator ( 44 ). The time step for the in tegrator is T = 3 / 2 with h = 3 / n v ar , m = n v ar / 2, and 40 inter-c hain steps. The RMSE for 100 instances of the sampling ﬁlter results are shown as b o x plots. The red line represents the median RMSE v alues across all instances and the central blue box represen ts the v ariance. The t wo vertical lines (whiskers) extend up to 1 . 5 times the height of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. In addition to con trolling the time step settings of the in tegrator, and tuning the num b er of steps of the c hain, we can use the Hilb ert integrator (with tuned step size) to p eriodically v alidate the ensembles obtained using other in tegrators, since the Hilbert integrator suﬀers less from outliers. A simple solution is to run the assimilation pro cess sev eral times and exclude outlier states b y creating a combined ensemble. Care must b e exercised, how ever, to not c hange the probabilit y densit y . These alternatives will be insp ected in depth in future w ork in the con text of more complex mo dels. 5.12. A highly nonline ar observation op er ator W e ha ve also tested the sampling ﬁlter capabilities in a very challenging setting: the exp onen tial observ ation op erator ( 38 ) is considered with a factor of r = 0 . 5. This factor leads to a large range of observ ation v alues (from e − 3 . 7 to e 6 . 2 ). In addition, small perturbations in the state v ariables cause v ery large changes in the corresp onding measurements. Both traditional metho ds EnKF and MLEF div erge when applied to this test, and consequently their results are not rep orted here. This test problem is challenging for the sampling ﬁlter as well and the symplectic in te- gration step sizes need to be tuned to ac hiev e con vergence. F or example, the num b er of steps tak en b y V erlet integrator has to b e increased to m = 60 while keeping the step-size ﬁxed to h = 0 . 01, to result in go od p erformance. The length of the tra jectory of the Hamiltonian system has to b e increased as w ell. F or the three-stage integrator a shorter tra jectory of the Hamiltonian system w orks well if the step size is suﬃcien tly reduced, e.g., h = 0 . 001, and m = 30. Empirical tuning of the V erlet, tw o-stage, and four-stage integrators pro ved to b e c hallenging with this observ ation op erator. Ho wev er, the three-stage integrator pro duced v ery satisfactory results with larger step sizes, as shown in Figure 21(a) . 39 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 Time RMSE Forecast (a) Three-stage integrator ( 44 ); h = 0 . 01 , m = 60 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 Time RMSE Forecast (b) F our-stage integrator ( 45 ); h = 0 . 001 , m = 30 Figure 21: Data assimilation results with the exp onential observ ation op erator ( 38 ) with a factor r = 0 . 2. The symplectic integrator used is indicated under each panel. The step size h and the n umber of steps m are indicated under eac h panel. The n um b er of inter-c hain steps is 30. The RMSE for 100 instances of the sampling ﬁlter results are shown as b o x plots. The red line represen ts the median RMSE v alues across all instances and the central blue b o x represents the v ariance. The tw o v ertical lines (whisk ers) extend up to 1 . 5 times the heigh t of the central b o x. The v alues exceeding the length of the whiskers are considered outliers (extremes) and are plotted as red crosses. The Hilb ert space in tegrator p erforms robustly and yields analyses that are as accurate as the ones for the simpler observ ation op erators; see Figure 22 . While the RMSE v alue ac hiev ed by the ﬁlter using the Hilb ert space in tegrator is relatively large, one can argue that this lev el is acceptable when dealing with large systems, nonlinear op erators where all other ﬁlters fail. The results in Figure 22 sho w that the analysis (of selected comp onen ts) follo ws the truth reasonably closely . 40 (a) x 1 (b) x 4 (c) x 7 (d) x 11 Figure 22: Integrator deﬁned on Hilb ert space is used; h = 0 . 001 , m = 30. The comp onents x 1 , x 4 , x 7 , x 11 of the state vector x are plotted. The num b er of inter-c hain steps is 30. Figure 23 plots the analyses obtained with the three-stage integrator sampling ﬁlter. A large num b er of steps is required to ac hiev e go o d results due to the large magnitude of observ ations. The Hilbert in tegrator op erates at a m uch low er cost, and can be used to p erio dically c heck the results obtained with the three-stage integrator to safeguard against outliers. 41 (a) x 1 (b) x 4 (c) x 7 (d) x 11 Figure 23: Three-stage integrator is used; h = 0 . 01 , m = 60. The comp onen ts x 1 , x 4 , x 7 , x 11 of the state v ector x are plotted. The n umber of inter-c hain steps is 30. The statistics of the results with Lorenz-96 mo del are summarized in T ables 1 through 4 . The results for 100 instances of the sampling ﬁlter, EnKF, and MLEF, o v er the time in terv al [8 , 10], are summarized in T able 1 and 2 . The results obtained with the exp onential observ ation op erator ( 38 ) with r = 0 . 5 are shown in T able 2 . In T able 1 the columns named “ Fixe d step ” present statistics obtained from exp erimen ts with the ﬁxed step size settings h = 0 . 01 and m = 10. The columns named “ Diﬀer ent step ” report statistics from the exp erimen ts where the w ork was equalized among the symplectic integrators. T ables 3 and 4 are shorter v ersions of T ables 1 , and 4 resp ectiv ely; only the results of the sampling ﬁlter with ﬁxed time step of the symplectic integrators are included and only the av erages and the standard deviations o v er the time interv al [8 , 10] are summarized. 5.13. Shal low water mo del on a spher e As a ﬁrst step tow ards large mo dels w e test the prop osed sampling ﬁlter on the shal- lo w w ater mo del on a sphere, using linear observ ation op erator where all comp onen ts are observ ed. The shallow w ater equations provide a simpliﬁed model of the atmosphere whic h de- scrib es the essential wa ve propagation mec hanisms found in general circulation mo dels 42 (GCMs) [ 33 ]. The shallow water equations in spherical co ordinates are giv en as ∂ u ∂ t + 1 a cos θ  u ∂ u ∂ λ + v cos θ ∂ u ∂ θ  −  f + u tan θ a  v + g a cos θ ∂ h ∂ λ = 0 , (39a) ∂ v ∂ t + 1 a cos θ  u ∂ v ∂ λ + v cos θ ∂ v ∂ θ  +  f + u tan θ a  u + g a ∂ h ∂ θ = 0 , (39b) ∂ h ∂ θ + 1 a cos θ  ∂ ( hu ) ∂ λ + ∂ ( hv cos θ ) ∂ θ  = 0 . (39c) The Coriolis parameter is giv en b y f = 2Ω sin θ , where Ω is the angular sp eed of the rotation of the Earth, and θ is latitudinal direction. The longitudinal direction is λ . The heigh t of the homogeneous atmosphere is represen ted b y h , the zonal and meridional wind comp onen ts are giv en b y u and v resp ectiv ely . The radius of the earth is a , and the gra vitational constan t is given by g . The space discretization follows the unstaggered T urkel-Zw as scheme [ 28 ]. The discretization has nl on = 72 no des in longitudinal direction and nl at = 36 no des in the latitudinal direction. The semi-discretization in space results in the follo wing discrete mo del: x k +1 = M t k → t k +1 ( x k , θ ) , k = 0 , . . . , N , (40) x 0 = x 0 ( θ ) . (41) The state space vector x in ( 40 ) com bines the zonal wind, the meridional wind, and the heigh t v ariables into the vector x ∈ R n v ar with n v ar = 3 × nlat × nlon. The time integration is conducted using an adaptive time-stepping algorithm. A reference initial condition is used to generate a reference tra jectory . Syn thetic observ ations are created from the reference tra jectory b y adding Gaussian noise with zero mean and ﬁxed standard deviation for each of the three comp onents. The lev el of observ ation noise for heigh t comp onen t is set to 1 . 5% of the av erage magnitude of the reference heigh t component in the reference initial condition. The lev el of observ ation noise for wind components is set to 10% of the a v erage magnitude of the reference wind comp onen t in the initial condition. The initial background state is created b y p erturbing the reference initial condition b y a Gaussian noise drawn from a mo delled bac kground error co v ariance matrix B 0 . The standard deviation of the background errors for the height comp onen t is 2% of the a verage magnitude of the reference heigh t comp onen t in the reference initial condition. The standard deviation of the bac kground errors for the wind comp onen ts is 15% of the a verage magnitude of the reference wind comp onent in the reference initial condition. The mo deled version of the background error cov ariance, B 0 , that accoun ts for correlations b et ween state v ariables is created as follows: • Start with a diagonal bac kground error cov ariance matrix with uncertain t y lev els as men tioned previously . • Apply the ensem ble Kalman ﬁlter for 48 hours. Synthetic initial ensemble is created b y adding zero-mean Gaussian noise to the reference initial condition with co v ariances set to the initial (diagonal) bac kground error cov ariance matrix. 43 • Decorrelate the ensemble-based co v ariances using a decorrelation matrix ρ with decor- relation distance L = 1000 k m . • Calculate B 0 b y a veraging the cov ariances ov er the last 6 hours. This metho d of creating a synthetic initial background error cov ariance matrix is totally empirical, but we found that the resulting bac kground error cov ariance matrix p erforms w ell for several algorithms including 4D V AR. Enhancing the quality of this bac kground error co v ariance matrix can b e done b y making use of the ensembles generated b y the sampling ﬁlter. In future work, we will in v estigate the p ossibility of estimating the background error co v ariances using the prop osed sampling ﬁlter. 5.14. R esults for shal low water mo del with line ar observations The assimilation time in terv al is 6 hours and there are hourly observ ations a v ailable. The n um b er of burn-in steps in the Marko v c hain is set to 50. W e use the t wo-stage symplectic in tegrator ( 43 ) with step size h = 0 . 01 and num b er of steps m = 10. The num b er of in ter-c hain steps is 10. The EnKF and the sampling ﬁlter results are sho wn in Figure 24 . As shown in Figures 24(b) , 24(d) , and 24(f ) the analysis is noisy and further tuning of the sampling ﬁlter pa- rameters is needed in order to outp erform the EnKF analysis. Parameter tuning for the sampling ﬁlter with this model will b e studied in the future. Moreov er ensemble based forecast co v ariances need to used for all analyses to improv e results. 44 5.2 5.4 5.6 5.8 6 6.2 6.4 x 10 4 (a) EnKF analysis: h 5.2 5.4 5.6 5.8 6 6.2 6.4 x 10 4 (b) Sampling ﬁlter analysis: h −20 −15 −10 −5 0 5 10 15 20 (c) EnKF analysis: u −25 −20 −15 −10 −5 0 5 10 15 20 (d) Sampling ﬁlter analysis: u −20 −15 −10 −5 0 5 10 15 20 (e) EnKF analysis: v −25 −20 −15 −10 −5 0 5 10 15 20 25 (f ) Sampling ﬁlter analysis: v Figure 24: Data assimilation results for SWE on the sphere with linear observ ations where all comp onen ts are observed. The plotted comp onen t of the state vector is indicated under each panel. The analysis shown is obtained after sequential assimilation of hourly observ ations for seven hours. Only state at the seven th hour is plotted. The sampling ﬁlter analysis is an o verage of analysis states obtained from 50 instances of the sampling ﬁlter. The symplectic integrator used is the tw o-stage integrator ( 43 ). The length of the Hamiltonian tra jectory is T = mh , with h = 0 . 01 , m = 10. The num b er of in ter-chain steps is 10. 6. Conclusions and F uture W ork This pap er prop oses a sampling ﬁlter for data assimilation where the analysis sc heme is replaced by sampling directly from the p osterior distribution. A Hybrid MCMC technique is employ ed to generate a representativ e analysis ensemble at eac h time. The sampling ﬁlter a v oids the need to dev elop tangen t linear or adjoint mo dels of the mo del solution op era- tor. The sampling ﬁlter can work with highly nonlinear observ ation op erators and pro vides analysis ensem bles that describ e non-Gaussian p osterior probability densities. The mean 45 of the generated posterior ensem ble pro vides the analysis (a minimum v ariance estimate of the state). The ensem ble cov ariance oﬀers an estimate of the analysis error co v ariance matrix and can b e used to quan tify the uncertaint y asso ciated with the analysis state. The implemen tation does not require the construction of full cov ariance matrices, whic h mak es the metho d attractiv e for large scale data assimilation problems with op erational mo dels and complex observ ation op erators. Numerical exp eriments are carried out with the Lorenz-96 mo del with sev eral observ ation op erators with diﬀeren t lev els of non-linearit y and smo othness. The sampling ﬁlter comp etes with EnKF for linear observ ations. F or nonlinear observ ations the results are very promising, and the new ﬁlter outp erforms b oth EnKF and MLEF. In addition, sampling ﬁlter contin ues to pro duce satisfactory results in cases where EnKF and MLEF fail. Large scale ensemble ﬁltering data assimilation problems are typically run on large par- allel machines. One imp ortan t challenge is the failure of subsets of no des, whic h terminates some of the ensemble mem b er runs, and leads to fewer ensem ble members b eing av ailable at the next time. Ov er sev eral cycles the num b er of ensem ble mem b ers can decrease consider- ably . The sampling strategy prop osed herein can b e used to replace dead ensem ble members in any parallel implemen tation of the EnKF. In addition, the sampling ﬁlter can b e used in com bination with classical ﬁlters b y building analysis ensembles that ha ve mem b ers given b y EnKF analysis (these members retain the history of the system) mixed with sampling mem b ers (whic h are consistent with the p osterior probability density , but add new directions to explore and can therefore a void ﬁlter divergence). The computational p erformance of the sampling ﬁlter dep ends on tuning its parameters, esp ecially the symplectic integration time step and the n um b er of steps tak en in the Mark ov c hain b et ween successive accepted ensemble mem b ers. F uture work will fo cus on reﬁning the strategies for parameter tuning in the con text of large op erational mo dels at high- resolution. W e also plan to p erform a side-b y-side comparison betw een the prop osed ﬁlter and the implicit sampling ﬁlter. Ac kno wledgmen ts This work was supp orted in part by aw ards NSF CCF–1218454 and AF OSR F A9550– 12–1–0293–DEF, and b y the Computational Science Lab oratory at Virginia T ech. A. Symplectic numerical in tegrators Here w e present the ﬁv e numerical integrators employ ed in this w ork. W e start with the standard p osition V erlet in tegrator in A.1 . The results of the standard V erlet are very sensitiv e to the choice of the time step. Three higher order integrators namely , t wo-stage ( A.2 ), three-stage ( A.3 ), and four-stage ( A.4 ) p osition splitting in tegrators, are tak en from [ 5 ]. These higher-order in tegrators lead to ﬁlters that are more stable and eﬃcient than V erlet. The last integrator tested ( A.5 ) is from [ 3 ] and is designed to w ork eﬃciently in inﬁnite dimensional state spaces, and to a v oid problems resulting from subtracting inﬁnitely 46 large num b ers related to the total energy of the Hamiltonian system for inﬁnite dimensional state spaces. All the in tegrators are applied to a Hamiltonian system of the form ( 18 ) [ 31 , 32 ]. A.1. Position V erlet inte gr ator One step of the p osition V erlet algorithm adv ances the solution of the Hamiltonian equations ( 18 ) from time t k to time t k +1 = t k + h as follo ws [ 31 ]: x k +1 / 2 = x k + h 2 M − 1 p k , (42a) p k +1 = p k − h ∇ x J ( x k +1 / 2 ) , (42b) x k +1 = x k +1 / 2 + h 2 M − 1 p k +1 . (42c) The optimal time step h is h ∝ (1 / n v ar ) 1 / 4 [ 2 ]. The exp eriments sho w that the step size should be small (close to zero) to mak e this in tegrator stable. It may still fail for high dimensionalit y and whenever complications are presen t in the target distributions. The w eakness of this simple integrator is illustrated in our n umerical exp erimen ts with highly nonlinear observ ation op erators. A.2. Two-stage inte gr ator One step of the tw o-stage algorithm adv ances the solution of the Hamiltonian equations ( 18 ) from time t k to time t k +1 = t k + h as follo ws [ 5 ]: x 1 = x k + ( a 1 h ) M − 1 p k , (43a) p 1 = p k − ( b 1 h ) ∇ x J ( x 1 ) , (43b) x 2 = x 1 + ( a 2 h ) M − 1 p 1 , (43c) p k +1 = p 1 − ( b 1 h ) ∇ x J ( x 2 ) , (43d) x k +1 = x 2 + ( a 2 h ) M − 1 p k +1 , (43e) where a 1 = 0 . 21132, a 2 = 1 − 2 a 1 , and b 1 = 0 . 5. The stabilit y of this time integrator is ac hiev ed for time step that lies in the interv al (0 , 2 . 6321480259) (units), that is, h should b e c hosen suc h that 0 < h < 2 . 6321480259 [ 5 ]. A.3. Thr e e-stage inte gr ator One step of the three-stage algorithm adv ances the solution of the Hamiltonian equations ( 18 ) from time t k to time t k +1 = t k + h b y the set of equations [ 5 ]: x 1 = x k + ( a 1 h ) M − 1 p k , (44a) p 1 = p k − ( b 1 h ) ∇ x J ( x 1 ) , (44b) x 2 = x 1 + ( a 2 h ) M − 1 p 1 , (44c) p 2 = p 1 − ( b 2 h ) ∇ x J ( x 2 ) , (44d) x 3 = x 2 + ( a 2 h ) M − 1 p 2 , (44e) p k +1 = p 2 − ( b 1 h ) ∇ x J ( x 3 ) , (44f ) x k +1 = x 3 + ( a 1 h ) M − 1 p k +1 , (44g) 47 where: a 1 = 0 . 11888010966548, a 2 = 0 . 5 − a 1 , b 1 = 0 . 29619504261126, and b 2 = 1 − 2 b 1 . The stability in terv al of the time step asso ciated with this time integrator is of length ≈ 4 . 67, that is, h should b e c hosen suc h that 0 < h < 4 . 67 [ 5 ]. A.4. F our-stage inte gr ator One step of the four-stage algorithm adv ances the solution of the Hamiltonian equations ( 18 ) from time t k to time t k +1 = t k + h as follo ws [ 5 ]: x 1 = x k + ( a 1 h ) M − 1 p k , (45a) p 1 = p k − ( b 1 h ) ∇ x J ( x 1 ) , (45b) x 2 = x 1 + ( a 2 h ) M − 1 p 1 , (45c) p 2 = p 1 − ( b 2 h ) ∇ x J ( x 2 ) , (45d) x 3 = x 2 + ( a 3 h ) M − 1 p 2 , (45e) p 3 = p 2 − ( b 2 h ) ∇ x J ( x 3 ) , (45f ) x 4 = x 3 + ( a 2 h ) M − 1 p 3 , (45g ) p k +1 = p 3 − ( b 1 h ) ∇ x J ( x 4 ) , (45h) x k +1 = x 4 + ( a 1 h ) M − 1 p k +1 , (45i) where: a 1 = 0 . 071353913450279725904, a 2 = 0 . 268458791161230105820, a 3 = 1 − 2 a 1 − 2 a 2 , b 1 = 0 . 1916678, and b 2 = 0 . 5 − b 1 . This in tegrator has a stabilit y in terv al of length ≈ 5 . 35, that is, h should b e c hosen suc h that 0 < h < 5 . 35 [ 5 ]. The time here has unsp eciﬁed units. Generally sp eaking, the high order in tegrators ( 43 , 44 , 45 ), pro vide more fa vorable and wider stability ranges for the time step. F or more on the stabilit y interv als of the time step settings of these high-order in tegrators, see [ 5 ]. A.5. Gener al inte gr ator deﬁne d on Hilb ert sp ac e One step of the Hilb ert integrator adv ances the solution of the Hamiltonian equations ( 18 ) from time t k to time t k +1 = t k + h as follo ws [ 3 ]: p 1 = p k − h 2 M − 1 ∇ x J ( x k ) , (46a) x k +1 = cos ( h ) x k + sin ( h ) p 1 , (46b) p 2 = − sin ( h ) x k + cos ( h ) p 1 , (46c) p k +1 = p 2 − h 2 M − 1 ∇ x J ( x k +1 ) . (46d) As with the standard position V erlet in tegrator the selection criterion of step size is not precisely deﬁned, ho w ev er, it is designed to w ork with ﬁnite (non-zero) steps in inﬁnite dimensional settings. Numerical results presen ted in Section 5 sho w that with careful tuning this in tegrator pro vides satisfactory results. 48 References [1] Anderson JL. 2001. A n ensemble adjustment Kalman ﬁlter for data assimilation . Monthly W eather Review, 129(12):2884–2903. [2] Besk os A, Pillai N, Rob erts G, Sanz-Serna JM, Stuart A. 2013. Optimal tuning of the hybrid Monte Carlo algorithm . Bernoulli, 19(5A):1501–1534. [3] Besk os A, Pinski FJ, Sanz-Serna JM, Stuart A. 2011 Hybrid Monte Carlo on Hilb ert sp ac es . Sto c hastic Pro cesses and their Applications, 121(10). [4] Bishop CH, Etherton BJ, and Ma jumdar SJ. 2001. A daptive sampling with the ensemble tr ansform Kalman ﬁlter. Part I: The or etic al asp e cts . Monthly W eather Review, 129:420–436. [5] Blanes S, Casas F, Sanz-Serna JM. 2014 Numeric al inte gr ators for the hybrid Monte Carlo metho d . arXiv Preprint [6] Burgers G, V an Leeu wen PJ, Ev ensen G. 1998. Analysis scheme in the ensemble Kalman ﬁlter . Monthly W eather Review, 126:1719–1724. [7] Chorin A, Morzfeld M, T u X. 2010. Implicit p article ﬁlters for data assimilation . Comm unications in Applied Mathematics and Computational Science, 5(2):221–240. [8] Cohn SE. 1997. An intr o duction to estimation the ory . Journal of the Meteorological So ciet y of Japan, 75:257–288. [9] Duane S, Kennedy AD, Pendleton BJ, and Row eth D. 1987. Hybrid Monte Carlo . Physics Letters B, 195(2):216–222. [10] Ev ensen G. 1994. Se quential data assimilation with a nonline ar quasi-ge ostr ophic mo del using Monte Carlo metho ds to for e c ast err or statistics . Journal of Geoph ysical Research, 99(C5):10143–10162. [11] Ev ensen G. 2003. The ensemble Kalman ﬁlter: The or etic al formulation and pr actic al implementation . Ocean Dynamics, 53. [12] Ev ensen G. 2007. Data assimilation: The ensemble Kalman ﬁlter . Springer. [13] Fisher M, Courtier P . 1995. Estimating the c ovarianc e matric es of analysis and for e c ast err or in variational data assimilation . Europ ean Cen ter for Medium-Range W eather F orecasts. [14] Girolami M, Calderhead B. 2011. R iemann manifold L angevin and Hamiltonian Monte Carlo metho ds . Journal of the Roy al Statistical So ciet y: Series B (Statistical Metho dology), 73(2):123–214. [15] Hamill TM, Whitaker JS, Snyder C. 2001. Distanc e-dep endent ﬁltering of b ackgr ound err or c ovarianc e estimates in an ensemble Kalman ﬁlter . Monthly W eather Review, 129:2776–2790. [16] Houtek amer PL, Mitchell HL. 1998. Data assimilation using an ensemble Kalman ﬁlter te chnique . Mon thly W eather Review, 126:796–811. [17] Houtek amer PL, Mitchell HL. 2001. A se quential ensemble Kalman ﬁlter for atmospheric data assimi- lation . Monthly W eather Review, 129:123–137. [18] Kalman RE. 1960. A new appr o ach to line ar ﬁltering and pr e diction pr oblems . T ransaction of the ASME- Journal of Basic Engineering, 82:35–45. [19] Kalman RE, Bucy RS. 1961. New r esults in line ar ﬁltering and pr e diction the ory . Journal of Basic Engineering, 83(1):95–108. [20] Kalna y E. 2002. Atmospheric mo deling, data assimilation and pr e dictability . Cambridge Univ ersity Press. [21] Liu JS. 2008. Monte Carlo str ate gies in scientiﬁc c omputing . Springer. [22] Lorenc AC. 1986. A nalysis metho ds for numeric al we ather pr e diction . Quarterly Journal of the Ro yal Meteorological So ciet y , 112(474). [23] Lorenz EN. 1996. Pr e dictability: A pr oblem p artly solve d . Pro c. Seminar on Predictability , volume 1. [24] Lorenz EN, Eman uel KA. 1998. Optimal sites for supplementary we ather observations: Simulation with a smal l mo del . Journal of the Atmospheric Sciences, 55(3):399–414. [25] Metrop olis N, Rosenbluth A W, Rosenbluth MN, T eller AH, T eller E. 1953. Equation of state c alculations by fast c omputing machines . The Journal of Chemical Physics, 21(6):1087–1092. [26] Neal RM. 1993. Pr ob abilistic infer enc e using Markov chain Monte Carlo metho ds . Departmen t of Computer Science, Universit y of T oronto T oronto, Ontario, Canada. [27] Neal RM. 2011. MCMC using Hamiltonian dynamics . Handb o ok of Marko v Chain Monte Carlo. 49 [28] Neta B, Giraldo FX, Nav on IM. 1997. Analysis of the Turkel-Zwas scheme for the two-dimensional shal low water e quations in spheric al c o or dinates . Journal of Computational Physics, Elsevier, 102–112. [29] Nino Ruiz ED, Sandu A, Anderson JL. 2014. An eﬃcient implementation of the ensemble Kalman ﬁlter b ase d on an iter ative Sherman–Morrison formula . Statistics and Computing, 1–17. [30] Ott E, Hun t BR, Szun yogh I, Zimin A V, Kostelic h EJ, Corazza M, Kalna y E, P atil DJ, Y orke JA. 2004. A lo c al ensemble kalman ﬁlter for atmospheric data assimilation . T ellus A, 56(5):415–428. [31] Sanz-Serna JM. 2014. Markov chain Monte Carlo and numeric al diﬀer ential e quations . Current Chal- lenges in Stability Issues for Numerical Diﬀerential Equations, 39–88. Springer. [32] Sanz-Serna JM , Calv o MP . 1994. Numeric al Hamiltonian pr oblems , volume 7. Chapman & Hall London. [33] St-Cyr A, Jablonowski C, Dennis JM, T ufo HM, Thomas SJ. 2007. A c omp arison of two shal low water mo dels with nonc onforming adaptive grids . Monthly W eather Review, 136:1898–1922. [34] Tierney L. 1994. Markov chains for exploring p osterior distributions . The Annals of Statistics, pages 1701–1728. [35] Tipp ett MK, Anderson JL, Bishop CH, Hamill TM, Whitaker JS. 2003. Ensemble squar e r o ot ﬁlters . Mon thly W eather Review, 131:1485–1490. [36] Whitak er JS, Hamill TM. 2002. Ensemble data assimilation without p erturb e d observations . Mon thly W eather Review, 130:1913–1924. [37] Zupanski M. 2005. Maximum likeliho o d ensemble ﬁlter: The or etic al asp e cts . Mon thly W eather Review, 133(6). [38] Zupanski M, Na von IM, Zupanski D. 2008. The maximum likeliho o d ensemble ﬁlter as a non- diﬀer entiable minimization algorithm . Quarterly Journal of the Roy al Meteorological So ciety , 134(633). 50 T able 1: RMS error statistics of exp erimen ts for assimilation time p oints 8 ≤ t ≤ 10 (after ﬁlter stabilizes) Observ ation Operator Statistics Integrator used with sampling ﬁlter T raditional Filters V erlet Two-stage Three-stage F our-stage Integrator on Hilb ert space Fixe d step Diﬀer ent step Fixe d step Diﬀerent step Fixe d step Diﬀerent step Fixe d step Diﬀer ent step Fixe d step Diﬀer ent step EnKF MLEF Linear observ ation operator Min 0.20266 0.192911 0.22185 0.19728 0.218158 0.193404 0.209657 0.194968 1.176385 1.225406 0.153021 1.167385 Max 0.334077 0.275264 0.371612 0.271278 0.345005 0.27841 0.350564 0.262573 1.740247 1.749126 0.477702 5.341185 Mean 0.266498 0.2247 0.269914 0.225638 0.263563 0.230312 0.263878 0.227828 1.413573 1.411903 0.305293 3.318066 Std 0.024486 0.016978 0.026238 0.013978 0.02482 0.017086 0.028621 0.016688 0.097418 0.094974 0.091245 1.413271 Mean+2 ∗ std 0.31547 0.258656 0.32239 0.253594 0.313203 0.264484 0.32112 0.261204 1.608409 1.601851 0.487783 6.144607 Mean − 2 ∗ std 0.217526 0.190744 0.217438 0.197682 0.213923 0.19614 0.206636 0.194452 1.218737 1.221955 0.122803 0.491525 Quadratic observ ation operator Min 0.293207 0.223915 0.229499 0.227357 0.257983 0.225245 0.26389 0.216563 1.803022 1.804532 3.340558 4.669170 Max 5.437528 5.710717 3.525457 5.148642 3.286706 5.412255 3.489366 5.042837 3.449921 3.414033 4.555472 6.005837 Mean 4.492176 2.802898 0.65887 1.161938 0.577684 1.132138 0.645417 0.80452 2.487458 2.328798 3.927612 5.118004 Std 0.836105 2.326364 0.828443 1.711295 0.662738 1.642162 0.8385 1.295648 0.432504 0.401300 0.323017 0.377321 Mean+2 ∗ std 6.164386 7.455626 2.315756 4.584528 1.90316 4.416462 2.322417 3.395816 3.352466 3.131399 4.573646 5.872646 Mean − 2 ∗ std 2.819966 -1.84983 -0.998016 -2.260652 -0.747792 -2.152186 -1.031583 -1.786776 1.62245 1.526198 3.281578 4.363361 Cubic observ ation operator Min 3.005206 0.453015 0.295747 0.475533 0.310439 0.526365 0.293922 0.426434 1.673592 1.249802 5.747101 4.720723 Max 5.195681 5.083948 3.904205 2.310545 2.606434 2.510949 1.32128 2.505575 2.659555 2.595746 11.620567 6.672957 Mean 4.159783 1.514528 1.300494 1.006574 0.606089 1.235655 0.454872 1.071686 2.213129 2.078567 8.768142 5.691641 Std 0.382443 0.970227 1.153967 0.482539 0.482256 0.519058 0.166534 0.481047 0.220747 0.248998 1.61792 0.665145 Mean+2 ∗ std 4.924669 3.454982 3.608428 1.971652 1.570601 2.273771 0.78794 2.03378 2.654623 2.576563 12.003982 7.021931 Mean − 2 ∗ std 3.394897 -0.425926 -1.00744 0.041496 -0.358423 0.197539 0.121804 0.109592 1.771635 1.580570 5.532302 4.361352 Absolute v alue observ ation operator Min 0.223854 0.215914 0.230685 0.215965 0.221819 0.207496 0.213172 0.21736 1.594316 1.350040 0.156432 2.814780 Max 3.770117 4.278896 4.221607 3.440444 3.783911 3.414783 3.255274 4.605238 3.186142 3.809884 0.489043 5.822207 Mean 0.390026 0.693854 0.488504 0.591975 0.514439 0.569576 0.401922 0.761355 2.240063 2.184913 0.235906 4.566194 Std 0.472502 0.938676 0.70202 0.74757 0.719887 0.721405 0.52001 0.993143 0.410658 0.493300 0.101655 0.794186 Mean+2 ∗ std 1.33503 2.571206 1.892544 2.087115 1.954213 2.012386 1.441942 2.747641 3.061379 3.171513 0.439216 6.154566 Mean − 2 ∗ std -0.554978 -1.183498 -0.915536 -0.903165 -0.925335 -0.873234 -0.638098 -1.224931 1.418747 1.198312 0.032596 2.977822 Quadratic observ ation operator with threshold Min 0.256795 0.208163 0.25226 0.201484 0.253754 0.207151 0.229083 0.203665 1.348305 1.326491 2.018643 4.691330 Max 4.585849 4.775978 0.512091 3.665279 0.402913 2.117683 0.415055 0.283818 1.91278 3.690276 3.296156 7.047109 Mean 3.406461 1.920009 0.303039 0.448207 0.295141 0.290576 0.303422 0.240194 1.579326 1.853404 2.801247 5.918242 Std 1.067004 1.927817 0.044306 0.686636 0.03201 0.235245 0.037552 0.016292 0.119657 0.456134 0.311357 0.603508 Mean+2 ∗ std 5.540469 5.775643 0.391651 1.821479 0.359161 0.761066 0.378526 0.272778 1.81864 2.765672 3.423961 7.125258 Mean − 2 ∗ std 1.272453 -1.935625 0.214427 -0.925065 0.231121 -0.179914 0.228318 0.20761 1.340012 0.941136 2.178533 4.711226 Exponential observ ation operator with r = 0 . 2 Min 0.318723 0.315953 0.313345 0.296624 0.309291 0.277073 0.321489 0.291672 1.42193 1.408169 2.028516 4.460268 Max 2.563976 3.075637 0.688663 0.43783 0.643646 0.612475 0.674257 0.439936 1.97735 2.081439 3.911317 7.736204 Mean 0.433829 0.453889 0.408104 0.348357 0.405423 0.349503 0.408271 0.348964 1.610456 1.661009 3.153979 5.713367 Std 0.267084 0.320384 0.05718 0.028835 0.055325 0.041517 0.05946 0.02787 0.10526 0.139997 0.551218 0.934198 Mean+2 ∗ std 0.967997 1.094657 0.522464 0.406027 0.516073 0.432537 0.527191 0.404704 1.820976 1.941003 4.256415 7.581764 Mean − 2 ∗ std -0.100339 -0.186879 0.293744 0.290687 0.294773 0.266469 0.289351 0.293224 1.399936 1.381015 2.051543 3.844971 51 T able 2: RMS error statistics of exp eriments for assimilation time p oints 8 ≤ t ≤ 10 (after ﬁlter stabilizes). The exp onen tial observ ation op erator with factor r = 0 . 5 is used. Statistics In tegrator used with sampling ﬁlter Three-stage ; h = 0 . 01 , m = 60 Hilb ert ; h = 0 . 001 , m = 30 Min 0.304178 1.234498 Max 2.671971 2.350684 Mean 0.439776 1.699096 Std 0.274643 0.250088 Mean+2 ∗ std 0.989062 2.199272 Mean − 2 ∗ std -0.109510 1.198920 T able 3: RMS error statistics of exp eriments for assimilation time p oints 8 ≤ t ≤ 10 (after ﬁlter stabilizes). Fixed step size h = 0 . 01 , m = 10. Observ ation Op erator Statistics In tegrator used with sampling ﬁlter T raditional Filters V erlet Tw o-stage Three-stage F our-stage Hilbert EnKF MLEF Linear Mean 0.266498 0.269914 0.263563 0.263878 1.413573 0.305293 3.318066 Std 0.024486 0.026238 0.02482 0.028621 0.097418 0.091245 1.413271 Quadratic Mean 4.492176 0.65887 0.577684 0.645417 2.487458 3.927612 5.118004 Std 0.836105 0.828443 0.662738 0.8385 0.432504 0.323017 0.377321 Cubic Mean 4.159783 1.300494 0.606089 0.454872 2.213129 8.768142 5.691641 Std 0.382443 1.153967 0.482256 0.166534 0.220747 1.61792 0.665145 Absolute v alue Mean 0.390026 0.488504 0.514439 0.401922 2.240063 0.235906 4.566194 Std 0.472502 0.70202 0.719887 0.52001 0.410658 0.101655 0.794186 Quadratic with threshold Mean 3.406461 0.303039 0.295141 0.303422 1.579326 2.801247 5.918242 Std 1.067004 0.044306 0.03201 0.037552 0.119657 0.311357 0.603508 Exponential with r = 0 . 2 Mean 0.433829 0.408104 0.405423 0.408271 1.610456 3.153979 5.713367 Std 0.267084 0.05718 0.055325 0.05946 0.10526 0.551218 0.934198 T able 4: RMS error statistics of exp eriments for assimilation time p oints 8 ≤ t ≤ 10 (after ﬁlter stabilizes). The exp onen tial observ ation op erator with factor r = 0 . 5 is used. Statistics In tegrator used with sampling ﬁlter Three-stage ; h = 0 . 01 , m = 60 Hilb ert ; h = 0 . 001 , m = 30 Mean 0.439776 1.699096 Std 0.274643 0.250088 52

A Sampling Filter for Non-Gaussian Data Assimilation

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment