@article{HastermannReinhardtKleinetal.2021,
author = {Hastermann, Gottfried and Reinhardt, Maria and Klein, Rupert and Reich, Sebastian},
title = {Balanced data assimilation for highly oscillatory mechanical systems},
series = {Communications in applied mathematics and computational science : CAMCoS},
volume = {16},
journal = {Communications in applied mathematics and computational science : CAMCoS},
number = {1},
publisher = {Mathematical Sciences Publishers},
address = {Berkeley},
issn = {1559-3940},
doi = {10.2140/camcos.2021.16.119},
pages = {119 -- 154},
year = {2021},
abstract = {Data assimilation algorithms are used to estimate the states of a dynamical system using partial and noisy observations. The ensemble Kalman filter has become a popular data assimilation scheme due to its simplicity and robustness for a wide range of application areas. Nevertheless, this filter also has limitations due to its inherent assumptions of Gaussianity and linearity, which can manifest themselves in the form of dynamically inconsistent state estimates. This issue is investigated here for balanced, slowly evolving solutions to highly oscillatory Hamiltonian systems which are prototypical for applications in numerical weather prediction. It is demonstrated that the standard ensemble Kalman filter can lead to state estimates that do not satisfy the pertinent balance relations and ultimately lead to filter divergence. Two remedies are proposed, one in terms of blended asymptotically consistent time-stepping schemes, and one in terms of minimization-based postprocessing methods. The effects of these modifications to the standard ensemble Kalman filter are discussed and demonstrated numerically for balanced motions of two prototypical Hamiltonian reference systems.},
language = {en}
}
@article{LeungLeutbecherReichetal.2020,
author = {Leung, Tsz Yan and Leutbecher, Martin and Reich, Sebastian and Shepherd, Theodore G.},
title = {Impact of the mesoscale range on error growth and the limits to atmospheric predictability},
series = {Journal of the atmospheric sciences},
volume = {77},
journal = {Journal of the atmospheric sciences},
number = {11},
publisher = {American Meteorological Soc.},
address = {Boston},
issn = {0022-4928},
doi = {10.1175/JAS-D-19-0346.1},
pages = {3769 -- 3779},
year = {2020},
abstract = {Global numerical weather prediction (NWP) models have begun to resolve the mesoscale k(-5/3) range of the energy spectrum, which is known to impose an inherently finite range of deterministic predictability per se as errors develop more rapidly on these scales than on the larger scales. However, the dynamics of these errors under the influence of the synoptic-scale k(-3) range is little studied. Within a perfect-model context, the present work examines the error growth behavior under such a hybrid spectrum in Lorenz's original model of 1969, and in a series of identical-twin perturbation experiments using an idealized two-dimensional barotropic turbulence model at a range of resolutions. With the typical resolution of today's global NWP ensembles, error growth remains largely uniform across scales. The theoretically expected fast error growth characteristic of a k(-5/3) spectrum is seen to be largely suppressed in the first decade of the mesoscale range by the synoptic-scale k(-3) range. However, it emerges once models become fully able to resolve features on something like a 20-km scale, which corresponds to a grid resolution on the order of a few kilometers.},
language = {en}
}
@article{AkhmatskayaBouRabeeReich2009,
author = {Akhmatskaya, Elena and Bou-Rabee, Nawaf and Reich, Sebastian},
title = {A comparison of generalized hybrid Monte Carlo methods with and without momentum flip},
issn = {0021-9991},
doi = {10.1016/j.jcp.2008.12.014},
year = {2009},
abstract = {The generalized hybrid Monte Carlo (GHMC) method combines Metropolis corrected constant energy simulations with a partial random refreshment step in the particle momenta. The standard detailed balance condition requires that momenta are negated upon rejection of a molecular dynamics proposal step. The implication is a trajectory reversal upon rejection, which is undesirable when interpreting GHMC as thermostated molecular dynamics. We show that a modified detailed balance condition can be used to implement GHMC without momentum flips. The same modification can be applied to the generalized shadow hybrid Monte Carlo (GSHMC) method. Numerical results indicate that GHMC/GSHMC implementations with momentum flip display a favorable behavior in terms of sampling efficiency, i.e., the traditional GHMC/GSHMC implementations with momentum flip got the advantage of a higher acceptance rate and faster decorrelation of Monte Carlo samples. The difference is more pronounced for GHMC. We also numerically investigate the behavior of the GHMC method as a Langevin-type thermostat. We find that the GHMC method without momentum flip interferes less with the underlying stochastic molecular dynamics in terms of autocorrelation functions and it to be preferred over the GHMC method with momentum flip. The same finding applies to GSHMC.},
language = {en}
}
@article{AkhmatskayaBouRabeeReich2009,
author = {Akhmatskaya, Elena and Bou-Rabee, Nawaf and Reich, Sebastian},
title = {Erratum to "A comparison of generalized hybrid Monte Carlo methods with and without momentum flip" [J. Comput. Phys. 228 (2009), S. 2256 - 2265]},
issn = {0021-9991},
doi = {10.1016/j.jcp.2009.06.039},
year = {2009},
abstract = {The generalized hybrid Monte Carlo (GHMC) method combines Metropolis corrected constant energy simulations with a partial random refreshment step in the particle momenta. The standard detailed balance condition requires that momenta are negated upon rejection of a molecular dynamics proposal step. The implication is a trajectory reversal upon rejection, which is undesirable when interpreting GHMC as thermostated molecular dynamics. We show that a modified detailed balance condition can be used to implement GHMC without momentum flips. The same modification can be applied to the generalized shadow hybrid Monte Carlo (GSHMC) method. Numerical results indicate that GHMC/GSHMC implementations with momentum flip display a favorable behavior in terms of sampling efficiency, i.e., the traditional GHMC/GSHMC implementations with momentum flip got the advantage of a higher acceptance rate and faster decorrelation of Monte Carlo samples. The difference is more pronounced for GHMC. We also numerically investigate the behavior of the GHMC method as a Langevin-type thermostat. We find that the GHMC method without momentum flip interferes less with the underlying stochastic molecular dynamics in terms of autocorrelation functions and it to be preferred over the GHMC method with momentum flip. The same finding applies to GSHMC.},
language = {en}
}
@article{LeimkuhlerReich2009,
author = {Leimkuhler, Benedict and Reich, Sebastian},
title = {A metropolis adjusted Nos{\´e}-Hoover thermostat},
issn = {0764-583X},
doi = {10.1051/M2an/2009023},
year = {2009},
abstract = {We present a Monte Carlo technique for sampling from the canonical distribution in molecular dynamics. The method is built upon the Nose-Hoover constant temperature formulation and the generalized hybrid Monte Carlo method. In contrast to standard hybrid Monte Carlo methods only the thermostat degree of freedom is stochastically resampled during a Monte Carlo step.},
language = {en}
}
@article{Reich2013,
author = {Reich, Sebastian},
title = {A nonparametric ensemble transform method for bayesian inference},
series = {SIAM journal on scientific computing},
volume = {35},
journal = {SIAM journal on scientific computing},
number = {4},
publisher = {Society for Industrial and Applied Mathematics},
address = {Philadelphia},
issn = {1064-8275},
doi = {10.1137/130907367},
pages = {A2013 -- A2024},
year = {2013},
abstract = {Many applications, such as intermittent data assimilation, lead to a recursive application of Bayesian inference within a Monte Carlo context. Popular data assimilation algorithms include sequential Monte Carlo methods and ensemble Kalman filters (EnKFs). These methods differ in the way Bayesian inference is implemented. Sequential Monte Carlo methods rely on importance sampling combined with a resampling step, while EnKFs utilize a linear transformation of Monte Carlo samples based on the classic Kalman filter. While EnKFs have proven to be quite robust even for small ensemble sizes, they are not consistent since their derivation relies on a linear regression ansatz. In this paper, we propose another transform method, which does not rely on any a priori assumptions on the underlying prior and posterior distributions. The new method is based on solving an optimal transportation problem for discrete random variables.},
language = {en}
}
@article{BergemannGottwaldReich2009,
author = {Bergemann, Kay and Gottwald, Georg and Reich, Sebastian},
title = {Ensemble propagation and continuous matrix factorization algorithms},
issn = {0035-9009},
doi = {10.1002/qj.457},
year = {2009},
abstract = {We consider the problem of propagating an ensemble of solutions and its characterization in terms of its mean and covariance matrix. We propose differential equations that lead to a continuous matrix factorization of the ensemble into a generalized singular value decomposition (SVD). The continuous factorization is applied to ensemble propagation under periodic rescaling (ensemble breeding) and under periodic Kalman analysis steps (ensemble Kalman filter). We also use the continuous matrix factorization to perform a re-orthogonalization of the ensemble after each time-step and apply the resulting modified ensemble propagation algorithm to the ensemble Kalman filter. Results from the Lorenz-96 model indicate that the re-orthogonalization of the ensembles leads to improved filter performance.},
language = {en}
}
@article{CotterHamPainetal.2009,
author = {Cotter, Colin J. and Ham, David A. and Pain, Christopher C. and Reich, Sebastian},
title = {LBB stability of a mixed Galerkin finite element pair for fluid flow simulations},
issn = {0021-9991},
doi = {10.1016/j.jcp.2008.09.014},
year = {2009},
abstract = {We introduce a new mixed finite element for solving the 2- and 3-dimensional wave equations and equations of incompressible flow. The element, which we refer to as P1(D)-P2, uses discontinuous piecewise linear functions for velocity and continuous piecewise quadratic functions for pressure. The aim of introducing the mixed formulation is to produce a new flexible element choice for triangular and tetrahedral meshes which satisfies the LBB stability condition and hence has no spurious zero-energy modes. The advantage of this particular element choice is that the mass matrix for velocity is block diagonal so it can be trivially inverted; it also allows the order of the pressure to be increased to quadratic whilst maintaining LBB stability which has benefits in geophysical applications with Coriolis forces. We give a normal mode analysis of the semi-discrete wave equation in one dimension which shows that the element pair is stable, and demonstrate that the element is stable with numerical integrations of the wave equation in two dimensions, an analysis of the resultant discrete Laplace operator in two and three dimensions on various meshes which shows that the element pair does not have any spurious modes. We provide convergence tests for the element pair which confirm that the element is stable since the convergence rate of the numerical solution is quadratic.},
language = {en}
}
@article{BergemannReich2010,
author = {Bergemann, Kay and Reich, Sebastian},
title = {A localization technique for ensemble Kalman filters},
issn = {0035-9009},
doi = {10.1002/Qj.591},
year = {2010},
abstract = {Ensemble Kalman filter techniques are widely used to assimilate observations into dynamical models. The phase- space dimension is typically much larger than the number of ensemble members, which leads to inaccurate results in the computed covariance matrices. These inaccuracies can lead, among other things, to spurious long-range correlations, which can be eliminated by Schur-product-based localization techniques. In this article, we propose a new technique for implementing such localization techniques within the class of ensemble transform/square-root Kalman filters. Our approach relies on a continuous embedding of the Kalman filter update for the ensemble members, i.e. we state an ordinary differential equation (ODE) with solutions that, over a unit time interval, are equivalent to the Kalman filter update. The ODE formulation forms a gradient system with the observations as a cost functional. Besides localization, the new ODE ensemble formulation should also find useful application in the context of nonlinear observation operators and observations that arrive continuously in time.},
language = {en}
}
@article{BergemannReich2010,
author = {Bergemann, Kay and Reich, Sebastian},
title = {A mollified ensemble Kalman filter},
issn = {0035-9009},
doi = {10.1002/Qj.672},
year = {2010},
abstract = {It is well recognized that discontinuous analysis increments of sequential data assimilation systems, such as ensemble Kalman filters, might lead to spurious high-frequency adjustment processes in the model dynamics. Various methods have been devised to spread out the analysis increments continuously over a fixed time interval centred about the analysis time. Among these techniques are nudging and incremental analysis updates (IAU). Here we propose another alternative, which may be viewed as a hybrid of nudging and IAU and which arises naturally from a recently proposed continuous formulation of the ensemble Kalman analysis step. A new slow-fast extension of the popular Lorenz-96 model is introduced to demonstrate the properties of the proposed mollified ensemble Kalman filter.},
language = {en}
}
@article{ShinZoellerHolschneideretal.2011,
author = {Shin, Seoleun and Z{\"o}ller, Gert and Holschneider, Matthias and Reich, Sebastian},
title = {A multigrid solver for modeling complex interseismic stress fields},
series = {Computers \& geosciences : an international journal devoted to the publication of papers on all aspects of geocomputation and to the distribution of computer programs and test data sets ; an official journal of the International Association for Mathematical Geology},
volume = {37},
journal = {Computers \& geosciences : an international journal devoted to the publication of papers on all aspects of geocomputation and to the distribution of computer programs and test data sets ; an official journal of the International Association for Mathematical Geology},
number = {8},
publisher = {Elsevier},
address = {Oxford},
issn = {0098-3004},
doi = {10.1016/j.cageo.2010.11.011},
pages = {1075 -- 1082},
year = {2011},
abstract = {We develop a multigrid, multiple time stepping scheme to reduce computational efforts for calculating complex stress interactions in a strike-slip 2D planar fault for the simulation of seismicity. The key elements of the multilevel solver are separation of length scale, grid-coarsening, and hierarchy. In this study the complex stress interactions are split into two parts: the first with a small contribution is computed on a coarse level, and the rest for strong interactions is on a fine level. This partition leads to a significant reduction of the number of computations. The reduction of complexity is even enhanced by combining the multigrid with multiple time stepping. Computational efficiency is enhanced by a factor of 10 while retaining a reasonable accuracy, compared to the original full matrix-vortex multiplication. The accuracy of solution and computational efficiency depend on a given cut-off radius that splits multiplications into the two parts. The multigrid scheme is constructed in such a way that it conserves stress in the entire half-space.},
language = {en}
}
@article{GreweBrinkopJoeckeletal.2014,
author = {Grewe, Volker and Brinkop, Sabine and Joeckel, Patrick and Shin, Seoleun and Reich, Sebastian and Yserentant, Harry},
title = {On the theory of mass conserving transformations for Lagrangian methods in 3D atmosphere-chemistry models},
series = {Meteorologische Zeitschrift},
volume = {23},
journal = {Meteorologische Zeitschrift},
number = {4},
publisher = {Schweizerbart},
address = {Stuttgart},
issn = {0941-2948},
doi = {10.1127/0941-2948/2014/0552},
pages = {441 -- 447},
year = {2014},
language = {en}
}
@article{ShinReichFrank2012,
author = {Shin, Seoleun and Reich, Sebastian and Frank, Jason},
title = {Hydrostatic Hamiltonian particle-mesh (HPM) methods for atmospheric modelling},
series = {Quarterly journal of the Royal Meteorological Society},
volume = {138},
journal = {Quarterly journal of the Royal Meteorological Society},
number = {666},
publisher = {Wiley-Blackwell},
address = {Hoboken},
issn = {0035-9009},
doi = {10.1002/qj.982},
pages = {1388 -- 1399},
year = {2012},
abstract = {We develop a hydrostatic Hamiltonian particle-mesh (HPM) method for efficient long-term numerical integration of the atmosphere. In the HPM method, the hydrostatic approximation is interpreted as a holonomic constraint for the vertical position of particles. This can be viewed as defining a set of vertically buoyant horizontal meshes, with the altitude of each mesh point determined so as to satisfy the hydrostatic balance condition and with particles modelling horizontal advection between the moving meshes. We implement the method in a vertical-slice model and evaluate its performance for the simulation of idealized linear and nonlinear orographic flow in both dry and moist environments. The HPM method is able to capture the basic features of the gravity wave to a degree of accuracy comparable with that reported in the literature. The numerical solution in the moist experiment indicates that the influence of moisture on wave characteristics is represented reasonably well and the reduction of momentum flux is in good agreement with theoretical analysis.},
language = {en}
}
@article{Reich2011,
author = {Reich, Sebastian},
title = {A dynamical systems framework for intermittent data assimilation},
series = {BIT : numerical mathematics ; the leading applied mathematics journal for all computational mathematicians},
volume = {51},
journal = {BIT : numerical mathematics ; the leading applied mathematics journal for all computational mathematicians},
number = {1},
publisher = {Springer},
address = {Dordrecht},
issn = {0006-3835},
doi = {10.1007/s10543-010-0302-4},
pages = {235 -- 249},
year = {2011},
abstract = {We consider the problem of discrete time filtering (intermittent data assimilation) for differential equation models and discuss methods for its numerical approximation. The focus is on methods based on ensemble/particle techniques and on the ensemble Kalman filter technique in particular. We summarize as well as extend recent work on continuous ensemble Kalman filter formulations, which provide a concise dynamical systems formulation of the combined dynamics-assimilation problem. Possible extensions to fully nonlinear ensemble/particle based filters are also outlined using the framework of optimal transportation theory.},
language = {en}
}
@article{SkeelReich2011,
author = {Skeel, R. D. and Reich, Sebastian},
title = {Corrected potential energy functions for constrained molecular dynamics},
series = {European physical journal special topics},
volume = {200},
journal = {European physical journal special topics},
number = {1},
publisher = {Springer},
address = {Heidelberg},
issn = {1951-6355},
doi = {10.1140/epjst/e2011-01518-8},
pages = {55 -- 72},
year = {2011},
abstract = {Atomic oscillations present in classical molecular dynamics restrict the step size that can be used. Multiple time stepping schemes offer only modest improvements, and implicit integrators are costly and inaccurate. The best approach may be to actually remove the highest frequency oscillations by constraining bond lengths and bond angles, thus permitting perhaps a 4-fold increase in the step size. However, omitting degrees of freedom produces errors in statistical averages, and rigid angles do not bend for strong excluded volume forces. These difficulties can be addressed by an enhanced treatment of holonomic constrained dynamics using ideas from papers of Fixman (1974) and Reich (1995, 1999). In particular, the 1995 paper proposes the use of "flexible" constraints, and the 1999 paper uses a modified potential energy function with rigid constraints to emulate flexible constraints. Presented here is a more direct and rigorous derivation of the latter approach, together with justification for the use of constraints in molecular modeling. With rigor comes limitations, so practical compromises are proposed: simplifications of the equations and their judicious application when assumptions are violated. Included are suggestions for new approaches.},
language = {en}
}
@article{AmezcuaIdeKalnayetal.2014,
author = {Amezcua, Javier and Ide, Kayo and Kalnay, Eugenia and Reich, Sebastian},
title = {Ensemble transform Kalman-Bucy filters},
series = {Quarterly journal of the Royal Meteorological Society},
volume = {140},
journal = {Quarterly journal of the Royal Meteorological Society},
number = {680},
publisher = {Wiley-Blackwell},
address = {Hoboken},
issn = {0035-9009},
doi = {10.1002/qj.2186},
pages = {995 -- 1004},
year = {2014},
abstract = {Two recent works have adapted the Kalman-Bucy filter into an ensemble setting. In the first formulation, the ensemble of perturbations is updated by the solution of an ordinary differential equation (ODE) in pseudo-time, while the mean is updated as in the standard Kalman filter. In the second formulation, the full ensemble is updated in the analysis step as the solution of single set of ODEs in pseudo-time. Neither requires matrix inversions except for the frequently diagonal observation error covariance. We analyse the behaviour of the ODEs involved in these formulations. We demonstrate that they stiffen for large magnitudes of the ratio of background error to observational error variance, and that using the integration scheme proposed in both formulations can lead to failure. A numerical integration scheme that is both stable and is not computationally expensive is proposed. We develop transform-based alternatives for these Bucy-type approaches so that the integrations are computed in ensemble space where the variables are weights (of dimension equal to the ensemble size) rather than model variables. Finally, the performance of our ensemble transform Kalman-Bucy implementations is evaluated using three models: the 3-variable Lorenz 1963 model, the 40-variable Lorenz 1996 model, and a medium complexity atmospheric general circulation model known as SPEEDY. The results from all three models are encouraging and warrant further exploration of these assimilation techniques.},
language = {en}
}
@article{BergemannReich2012,
author = {Bergemann, Kay and Reich, Sebastian},
title = {An ensemble Kalman-Bucy filter for continuous data assimilation},
series = {Meteorologische Zeitschrift},
volume = {21},
journal = {Meteorologische Zeitschrift},
number = {3},
publisher = {Schweizerbart},
address = {Stuttgart},
issn = {0941-2948},
doi = {10.1127/0941-2948/2012/0307},
pages = {213 -- 219},
year = {2012},
abstract = {The ensemble Kalman filter has emerged as a promising filter algorithm for nonlinear differential equations subject to intermittent observations. In this paper, we extend the well-known Kalman-Bucy filter for linear differential equations subject to continous observations to the ensemble setting and nonlinear differential equations. The proposed filter is called the ensemble Kalman-Bucy filter and its performance is demonstrated for a simple mechanical model (Langevin dynamics) subject to incremental observations of its velocity.},
language = {en}
}
@article{Reich2012,
author = {Reich, Sebastian},
title = {A Gaussian-mixture ensemble transform filter},
series = {Quarterly journal of the Royal Meteorological Society},
volume = {138},
journal = {Quarterly journal of the Royal Meteorological Society},
number = {662},
publisher = {Wiley-Blackwell},
address = {Malden},
issn = {0035-9009},
doi = {10.1002/qj.898},
pages = {222 -- 233},
year = {2012},
abstract = {We generalize the popular ensemble Kalman filter to an ensemble transform filter, in which the prior distribution can take the form of a Gaussian mixture or a Gaussian kernel density estimator. The design of the filter is based on a continuous formulation of the Bayesian filter analysis step. We call the new filter algorithm the ensemble Gaussian-mixture filter (EGMF). The EGMF is implemented for three simple test problems (Brownian dynamics in one dimension, Langevin dynamics in two dimensions and the three-dimensional Lorenz-63 model). It is demonstrated that the EGMF is capable of tracking systems with non-Gaussian uni- and multimodal ensemble distributions.},
language = {en}
}
@article{AcevedoReichCubasch2016,
author = {Acevedo, Walter and Reich, Sebastian and Cubasch, Ulrich},
title = {Towards the assimilation of tree-ring-width records using ensemble Kalman filtering techniques},
series = {Climate dynamics : observational, theoretical and computational research on the climate system},
volume = {46},
journal = {Climate dynamics : observational, theoretical and computational research on the climate system},
publisher = {Springer},
address = {New York},
issn = {0930-7575},
doi = {10.1007/s00382-015-2683-1},
pages = {1909 -- 1920},
year = {2016},
abstract = {This paper investigates the applicability of the Vaganov-Shashkin-Lite (VSL) forward model for tree-ring-width chronologies as observation operator within a proxy data assimilation (DA) setting. Based on the principle of limiting factors, VSL combines temperature and moisture time series in a nonlinear fashion to obtain simulated TRW chronologies. When used as observation operator, this modelling approach implies three compounding, challenging features: (1) time averaging, (2) "switching recording" of 2 variables and (3) bounded response windows leading to "thresholded response". We generate pseudo-TRW observations from a chaotic 2-scale dynamical system, used as a cartoon of the atmosphere-land system, and attempt to assimilate them via ensemble Kalman filtering techniques. Results within our simplified setting reveal that VSL's nonlinearities may lead to considerable loss of assimilation skill, as compared to the utilization of a time-averaged (TA) linear observation operator. In order to understand this undesired effect, we embed VSL's formulation into the framework of fuzzy logic (FL) theory, which thereby exposes multiple representations of the principle of limiting factors. DA experiments employing three alternative growth rate functions disclose a strong link between the lack of smoothness of the growth rate function and the loss of optimality in the estimate of the TA state. Accordingly, VSL's performance as observation operator can be enhanced by resorting to smoother FL representations of the principle of limiting factors. This finding fosters new interpretations of tree-ring-growth limitation processes.},
language = {en}
}
@article{GregoryCotterReich2016,
author = {Gregory, A. and Cotter, C. J. and Reich, Sebastian},
title = {MULTILEVEL ENSEMBLE TRANSFORM PARTICLE FILTERING},
series = {SIAM journal on scientific computing},
volume = {38},
journal = {SIAM journal on scientific computing},
publisher = {Society for Industrial and Applied Mathematics},
address = {Philadelphia},
issn = {1064-8275},
doi = {10.1137/15M1038232},
pages = {A1317 -- A1338},
year = {2016},
abstract = {This paper extends the multilevel Monte Carlo variance reduction technique to nonlinear filtering. In particular, multilevel Monte Carlo is applied to a certain variant of the particle filter, the ensemble transform particle filter (EPTF). A key aspect is the use of optimal transport methods to re-establish correlation between coarse and fine ensembles after resampling; this controls the variance of the estimator. Numerical examples present a proof of concept of the effectiveness of the proposed method, demonstrating significant computational cost reductions (relative to the single-level ETPF counterpart) in the propagation of ensembles.},
language = {en}
}
@article{EscribanoAkhmatskayaReichetal.2015,
author = {Escribano, Bruno and Akhmatskaya, Elena and Reich, Sebastian and Azpiroz, Jon M.},
title = {Multiple-time-stepping generalized hybrid Monte Carlo methods},
series = {Journal of computational physics},
volume = {280},
journal = {Journal of computational physics},
publisher = {Elsevier},
address = {San Diego},
issn = {0021-9991},
doi = {10.1016/j.jcp.2014.08.052},
pages = {1 -- 20},
year = {2015},
abstract = {Performance of the generalized shadow hybrid Monte Carlo (GSHMC) method [1], which proved to be superior in sampling efficiency over its predecessors [2-4], molecular dynamics and hybrid Monte Carlo, can be further improved by combining it with multi-time-stepping (MTS) and mollification of slow forces. We demonstrate that the comparatively simple modifications of the method not only lead to better performance of GSHMC itself but also allow for beating the best performed methods, which use the similar force splitting schemes. In addition we show that the same ideas can be successfully applied to the conventional generalized hybrid Monte Carlo method (GHMC). The resulting methods, MTS-GHMC and MTS-GSHMC, provide accurate reproduction of thermodynamic and dynamical properties, exact temperature control during simulation and computational robustness and efficiency. MTS-GHMC uses a generalized momentum update to achieve weak stochastic stabilization to the molecular dynamics (MD) integrator. MTS-GSHMC adds the use of a shadow (modified) Hamiltonian to filter the MD trajectories in the HMC scheme. We introduce a new shadow Hamiltonian formulation adapted to force-splitting methods. The use of such Hamiltonians improves the acceptance rate of trajectories and has a strong impact on the sampling efficiency of the method. Both methods were implemented in the open-source MD package ProtoMol and were tested on a water and a protein systems. Results were compared to those obtained using a Langevin Molly (LM) method [5] on the same systems. The test results demonstrate the superiority of the new methods over LM in terms of stability, accuracy and sampling efficiency. This suggests that putting the MTS approach in the framework of hybrid Monte Carlo and using the natural stochasticity offered by the generalized hybrid Monte Carlo lead to improving stability of MTS and allow for achieving larger step sizes in the simulation of complex systems.},
language = {en}
}
@article{AizingerKornGiorgettaetal.2015,
author = {Aizinger, Vadym and Korn, Peter and Giorgetta, Marco and Reich, Sebastian},
title = {Large-scale turbulence modelling via alpha-regularisation for atmospheric simulations},
series = {Journal of turbulence},
volume = {16},
journal = {Journal of turbulence},
number = {4},
publisher = {Routledge, Taylor \& Francis Group},
address = {Abingdon},
issn = {1468-5248},
doi = {10.1080/14685248.2014.991443},
pages = {367 -- 391},
year = {2015},
abstract = {We study the possibility of obtaining a computational turbulence model by means of non-dissipative regularisation of the compressible atmospheric equations for climate-type applications. We use an -regularisation (Lagrangian averaging) of the atmospheric equations. For the hydrostatic and compressible atmospheric equations discretised using a finite volume method on unstructured grids, deterministic and non-deterministic numerical experiments are conducted to compare the individual solutions and the statistics of the regularised equations to those of the original model. The impact of the regularisation parameter is investigated. Our results confirm the principal compatibility of -regularisation with atmospheric dynamics and encourage further investigations within atmospheric model including complex physical parametrisations.},
language = {en}
}
@article{AcevedoFallahReichetal.2017,
author = {Acevedo, Walter and Fallah, Bijan and Reich, Sebastian and Cubasch, Ulrich},
title = {Assimilation of pseudo-tree-ring-width observations into an atmospheric general circulation model},
series = {Climate of the past : an interactive open access journal of the European Geosciences Union},
volume = {13},
journal = {Climate of the past : an interactive open access journal of the European Geosciences Union},
publisher = {Copernicus},
address = {G{\"o}ttingen},
issn = {1814-9324},
doi = {10.5194/cp-13-545-2017},
pages = {545 -- 557},
year = {2017},
abstract = {Paleoclimate data assimilation (DA) is a promising technique to systematically combine the information from climate model simulations and proxy records. Here, we investigate the assimilation of tree-ring-width (TRW) chronologies into an atmospheric global climate model using ensemble Kalman filter (EnKF) techniques and a process-based tree-growth forward model as an observation operator. Our results, within a perfect-model experiment setting, indicate that the "online DA" approach did not outperform the "off-line" one, despite its considerable additional implementation complexity. On the other hand, it was observed that the nonlinear response of tree growth to surface temperature and soil moisture does deteriorate the operation of the time-averaged EnKF methodology. Moreover, for the first time we show that this skill loss appears significantly sensitive to the structure of the growth rate function, used to represent the principle of limiting factors (PLF) within the forward model. In general, our experiments showed that the error reduction achieved by assimilating pseudo-TRW chronologies is modulated by the magnitude of the yearly internal variability in themodel. This result might help the dendrochronology community to optimize their sampling efforts.},
language = {en}
}
@article{SeeligRabeMalemShinitskietal.2020,
author = {Seelig, Stefan A. and Rabe, Maximilian Michael and Malem-Shinitski, Noa and Risse, Sarah and Reich, Sebastian and Engbert, Ralf},
title = {Bayesian parameter estimation for the SWIFT model of eye-movement control during reading},
series = {Journal of mathematical psychology},
volume = {95},
journal = {Journal of mathematical psychology},
publisher = {Elsevier},
address = {San Diego},
issn = {0022-2496},
doi = {10.1016/j.jmp.2019.102313},
pages = {32},
year = {2020},
abstract = {Process-oriented theories of cognition must be evaluated against time-ordered observations. Here we present a representative example for data assimilation of the SWIFT model, a dynamical model of the control of fixation positions and fixation durations during natural reading of single sentences. First, we develop and test an approximate likelihood function of the model, which is a combination of a spatial, pseudo-marginal likelihood and a temporal likelihood obtained by probability density approximation Second, we implement a Bayesian approach to parameter inference using an adaptive Markov chain Monte Carlo procedure. Our results indicate that model parameters can be estimated reliably for individual subjects. We conclude that approximative Bayesian inference represents a considerable step forward for computational models of eye-movement control, where modeling of individual data on the basis of process-based dynamic models has not been possible so far.},
language = {en}
}
@article{MalemShinitskiOpperReichetal.2020,
author = {Malem-Shinitski, Noa and Opper, Manfred and Reich, Sebastian and Schwetlick, Lisa and Seelig, Stefan A. and Engbert, Ralf},
title = {A mathematical model of local and global attention in natural scene viewing},
series = {PLoS Computational Biology : a new community journal},
volume = {16},
journal = {PLoS Computational Biology : a new community journal},
number = {12},
publisher = {PLoS},
address = {San Fransisco},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1007880},
pages = {21},
year = {2020},
abstract = {Author summary
Switching between local and global attention is a general strategy in human information processing. We investigate whether this strategy is a viable approach to model sequences of fixations generated by a human observer in a free viewing task with natural scenes. Variants of the basic model are used to predict the experimental data based on Bayesian inference. Results indicate a high predictive power for both aggregated data and individual differences across observers. The combination of a novel model with state-of-the-art Bayesian methods lends support to our two-state model using local and global internal attention states for controlling eye movements.
Understanding the decision process underlying gaze control is an important question in cognitive neuroscience with applications in diverse fields ranging from psychology to computer vision. The decision for choosing an upcoming saccade target can be framed as a selection process between two states: Should the observer further inspect the information near the current gaze position (local attention) or continue with exploration of other patches of the given scene (global attention)? Here we propose and investigate a mathematical model motivated by switching between these two attentional states during scene viewing. The model is derived from a minimal set of assumptions that generates realistic eye movement behavior. We implemented a Bayesian approach for model parameter inference based on the model's likelihood function. In order to simplify the inference, we applied data augmentation methods that allowed the use of conjugate priors and the construction of an efficient Gibbs sampler. This approach turned out to be numerically efficient and permitted fitting interindividual differences in saccade statistics. Thus, the main contribution of our modeling approach is two-fold; first, we propose a new model for saccade generation in scene viewing. Second, we demonstrate the use of novel methods from Bayesian inference in the field of scan path modeling.},
language = {en}
}
@article{deWiljesPathirajaReich2020,
author = {de Wiljes, Jana and Pathiraja, Sahani Darschika and Reich, Sebastian},
title = {Ensemble transform algorithms for nonlinear smoothing problems},
series = {SIAM journal on scientific computing},
volume = {42},
journal = {SIAM journal on scientific computing},
number = {1},
publisher = {Society for Industrial and Applied Mathematics},
address = {Philadelphia},
issn = {1064-8275},
doi = {10.1137/19M1239544},
pages = {A87 -- A114},
year = {2020},
abstract = {Several numerical tools designed to overcome the challenges of smoothing in a non-linear and non-Gaussian setting are investigated for a class of particle smoothers. The considered family of smoothers is induced by the class of linear ensemble transform filters which contains classical filters such as the stochastic ensemble Kalman filter, the ensemble square root filter, and the recently introduced nonlinear ensemble transform filter. Further the ensemble transform particle smoother is introduced and particularly highlighted as it is consistent in the particle limit and does not require assumptions with respect to the family of the posterior distribution. The linear update pattern of the considered class of linear ensemble transform smoothers allows one to implement important supplementary techniques such as adaptive spread corrections, hybrid formulations, and localization in order to facilitate their application to complex estimation problems. These additional features are derived and numerically investigated for a sequence of increasingly challenging test problems.},
language = {en}
}
@article{PathirajaReichStannat2021,
author = {Pathiraja, Sahani Darschika and Reich, Sebastian and Stannat, Wilhelm},
title = {McKean-Vlasov SDEs in nonlinear filtering},
series = {SIAM journal on control and optimization : a publication of the Society for Industrial and Applied Mathematics},
volume = {59},
journal = {SIAM journal on control and optimization : a publication of the Society for Industrial and Applied Mathematics},
number = {6},
publisher = {Society for Industrial and Applied Mathematics},
address = {Philadelphia},
issn = {0363-0129},
doi = {10.1137/20M1355197},
pages = {4188 -- 4215},
year = {2021},
abstract = {Various particle filters have been proposed over the last couple of decades with the common feature that the update step is governed by a type of control law. This feature makes them an attractive alternative to traditional sequential Monte Carlo which scales poorly with the state dimension due to weight degeneracy. This article proposes a unifying framework that allows us to systematically derive the McKean-Vlasov representations of these filters for the discrete time and continuous time observation case, taking inspiration from the smooth approximation of the data considered in [D. Crisan and J. Xiong, Stochastics, 82 (2010), pp. 53-68; J. M. Clark and D. Crisan, Probab. Theory Related Fields, 133 (2005), pp. 43-56]. We consider three filters that have been proposed in the literature and use this framework to derive Ito representations of their limiting forms as the approximation parameter delta -> 0. All filters require the solution of a Poisson equation defined on R-d, for which existence and uniqueness of solutions can be a nontrivial issue. We additionally establish conditions on the signal-observation system that ensures well-posedness of the weighted Poisson equation arising in one of the filters.},
language = {en}
}
@article{LeungLeutbecherReichetal.2021,
author = {Leung, Tsz Yan and Leutbecher, Martin and Reich, Sebastian and Shepherd, Theodore G.},
title = {Forecast verification},
series = {Quarterly journal of the Royal Meteorological Society},
volume = {147},
journal = {Quarterly journal of the Royal Meteorological Society},
number = {739},
publisher = {Wiley},
address = {Hoboken},
issn = {0035-9009},
doi = {10.1002/qj.4120},
pages = {3124 -- 3134},
year = {2021},
abstract = {The philosophy of forecast verification is rather different between deterministic and probabilistic verification metrics: generally speaking, deterministic metrics measure differences, whereas probabilistic metrics assess reliability and sharpness of predictive distributions. This article considers the root-mean-square error (RMSE), which can be seen as a deterministic metric, and the probabilistic metric Continuous Ranked Probability Score (CRPS), and demonstrates that under certain conditions, the CRPS can be mathematically expressed in terms of the RMSE when these metrics are aggregated. One of the required conditions is the normality of distributions. The other condition is that, while the forecast ensemble need not be calibrated, any bias or over/underdispersion cannot depend on the forecast distribution itself. Under these conditions, the CRPS is a fraction of the RMSE, and this fraction depends only on the heteroscedasticity of the ensemble spread and the measures of calibration. The derived CRPS-RMSE relationship for the case of perfect ensemble reliability is tested on simulations of idealised two-dimensional barotropic turbulence. Results suggest that the relationship holds approximately despite the normality condition not being met.},
language = {en}
}
@article{WormellReich2021,
author = {Wormell, Caroline L. and Reich, Sebastian},
title = {Spectral convergence of diffusion maps},
series = {SIAM journal on numerical analysis / Society for Industrial and Applied Mathematics},
volume = {59},
journal = {SIAM journal on numerical analysis / Society for Industrial and Applied Mathematics},
number = {3},
publisher = {Society for Industrial and Applied Mathematics},
address = {Philadelphia},
issn = {0036-1429},
doi = {10.1137/20M1344093},
pages = {1687 -- 1734},
year = {2021},
abstract = {Diffusion maps is a manifold learning algorithm widely used for dimensionality reduction. Using a sample from a distribution, it approximates the eigenvalues and eigenfunctions of associated Laplace-Beltrami operators. Theoretical bounds on the approximation error are, however, generally much weaker than the rates that are seen in practice. This paper uses new approaches to improve the error bounds in the model case where the distribution is supported on a hypertorus. For the data sampling (variance) component of the error we make spatially localized compact embedding estimates on certain Hardy spaces; we study the deterministic (bias) component as a perturbation of the Laplace-Beltrami operator's associated PDE and apply relevant spectral stability results. Using these approaches, we match long-standing pointwise error bounds for both the spectral data and the norm convergence of the operator discretization. We also introduce an alternative normalization for diffusion maps based on Sinkhorn weights. This normalization approximates a Langevin diffusion on the sample and yields a symmetric operator approximation. We prove that it has better convergence compared with the standard normalization on flat domains, and we present a highly efficient rigorous algorithm to compute the Sinkhorn weights.},
language = {en}
}
@article{EngbertRabeKliegletal.2021,
author = {Engbert, Ralf and Rabe, Maximilian Michael and Kliegl, Reinhold and Reich, Sebastian},
title = {Sequential data assimilation of the stochastic SEIR epidemic model for regional COVID-19 dynamics},
series = {Bulletin of mathematical biology : official journal of the Society for Mathematical Biology},
volume = {83},
journal = {Bulletin of mathematical biology : official journal of the Society for Mathematical Biology},
number = {1},
publisher = {Springer},
address = {New York},
issn = {0092-8240},
doi = {10.1007/s11538-020-00834-8},
pages = {16},
year = {2021},
abstract = {Newly emerging pandemics like COVID-19 call for predictive models to implement precisely tuned responses to limit their deep impact on society. Standard epidemic models provide a theoretically well-founded dynamical description of disease incidence. For COVID-19 with infectiousness peaking before and at symptom onset, the SEIR model explains the hidden build-up of exposed individuals which creates challenges for containment strategies. However, spatial heterogeneity raises questions about the adequacy of modeling epidemic outbreaks on the level of a whole country. Here, we show that by applying sequential data assimilation to the stochastic SEIR epidemic model, we can capture the dynamic behavior of outbreaks on a regional level. Regional modeling, with relatively low numbers of infected and demographic noise, accounts for both spatial heterogeneity and stochasticity. Based on adapted models, short-term predictions can be achieved. Thus, with the help of these sequential data assimilation methods, more realistic epidemic models are within reach.},
language = {en}
}
@article{HuangHuangReichetal.2022,
author = {Huang, Daniel Zhengyu and Huang, Jiaoyang and Reich, Sebastian and Stuart, Andrew M.},
title = {Efficient derivative-free Bayesian inference for large-scale inverse problems},
series = {Inverse problems : an international journal of inverse problems, inverse methods and computerised inversion of data},
volume = {38},
journal = {Inverse problems : an international journal of inverse problems, inverse methods and computerised inversion of data},
number = {12},
publisher = {IOP Publ. Ltd.},
address = {Bristol},
issn = {0266-5611},
doi = {10.1088/1361-6420/ac99fa},
pages = {40},
year = {2022},
abstract = {We consider Bayesian inference for large-scale inverse problems, where computational challenges arise from the need for repeated evaluations of an expensive forward model. This renders most Markov chain Monte Carlo approaches infeasible, since they typically require O(10(4)) model runs, or more. Moreover, the forward model is often given as a black box or is impractical to differentiate. Therefore derivative-free algorithms are highly desirable. We propose a framework, which is built on Kalman methodology, to efficiently perform Bayesian inference in such inverse problems. The basic method is based on an approximation of the filtering distribution of a novel mean-field dynamical system, into which the inverse problem is embedded as an observation operator. Theoretical properties are established for linear inverse problems, demonstrating that the desired Bayesian posterior is given by the steady state of the law of the filtering distribution of the mean-field dynamical system, and proving exponential convergence to it. This suggests that, for nonlinear problems which are close to Gaussian, sequentially computing this law provides the basis for efficient iterative methods to approximate the Bayesian posterior. Ensemble methods are applied to obtain interacting particle system approximations of the filtering distribution of the mean-field model; and practical strategies to further reduce the computational and memory cost of the methodology are presented, including low-rank approximation and a bi-fidelity approach. The effectiveness of the framework is demonstrated in several numerical experiments, including proof-of-concept linear/nonlinear examples and two large-scale applications: learning of permeability parameters in subsurface flow; and learning subgrid-scale parameters in a global climate model. Moreover, the stochastic ensemble Kalman filter and various ensemble square-root Kalman filters are all employed and are compared numerically. The results demonstrate that the proposed method, based on exponential convergence to the filtering distribution of a mean-field dynamical system, is competitive with pre-existing Kalman-based methods for inverse problems.},
language = {en}
}
@article{MolkenthinDonnerReichetal.2022,
author = {Molkenthin, Christian and Donner, Christian and Reich, Sebastian and Z{\"o}ller, Gert and Hainzl, Sebastian and Holschneider, Matthias and Opper, Manfred},
title = {GP-ETAS: semiparametric Bayesian inference for the spatio-temporal epidemic type aftershock sequence model},
series = {Statistics and Computing},
volume = {32},
journal = {Statistics and Computing},
number = {2},
publisher = {Springer},
address = {Dordrecht},
issn = {0960-3174},
doi = {10.1007/s11222-022-10085-3},
pages = {25},
year = {2022},
abstract = {The spatio-temporal epidemic type aftershock sequence (ETAS) model is widely used to describe the self-exciting nature of earthquake occurrences. While traditional inference methods provide only point estimates of the model parameters, we aim at a fully Bayesian treatment of model inference, allowing naturally to incorporate prior knowledge and uncertainty quantification of the resulting estimates. Therefore, we introduce a highly flexible, non-parametric representation for the spatially varying ETAS background intensity through a Gaussian process (GP) prior. Combined with classical triggering functions this results in a new model formulation, namely the GP-ETAS model. We enable tractable and efficient Gibbs sampling by deriving an augmented form of the GP-ETAS inference problem. This novel sampling approach allows us to assess the posterior model variables conditioned on observed earthquake catalogues, i.e., the spatial background intensity and the parameters of the triggering function. Empirical results on two synthetic data sets indicate that GP-ETAS outperforms standard models and thus demonstrate the predictive power for observed earthquake catalogues including uncertainty quantification for the estimated parameters. Finally, a case study for the l'Aquila region, Italy, with the devastating event on 6 April 2009, is presented.},
language = {en}
}