Mathias Rousset

Mathias Rousset, Giovanni Samaey

We discuss variance reduced simulations for an individual-based model of chemotaxis of bacteria with internal dynamics. The variance reduction is achieved via a coupling of this model with a simpler process in which the internal dynamics has been replaced by a direct gradient sensing of the chemoattractants concentrations. In the companion paper \cite{limits}, we have rigorously shown, using a pathwise probabilistic technique, that both processes converge towards the same advection-diffusion process in the diffusive asymptotics. In this work, a direct coupling is achieved between paths of individual bacteria simulated by both models, by using the same sets of random numbers in both simulations. This coupling is used to construct a hybrid scheme with reduced variance. We first compute a deterministic solution of the kinetic density description of the direct gradient sensing model; the deviations due to the presence of internal dynamics are then evaluated via the coupled individual-based simulations. We show that the resulting variance reduction is \emph{asymptotic}, in the sense that, in the diffusive asymptotics, the difference between the two processes has a variance which vanishes according to the small parameter.

Mathias Rousset, Giovanni Samaey

We discuss velocity-jump models for chemotaxis of bacteria with an internal state that allows the velocity jump rate to depend on the memory of the chemoattractant concentration along their path of motion. Using probabilistic techniques, we provide a pathwise result that shows that the considered process converges to an advection-diffusion process in the (long-time) diffusion limit. We also (re-)prove using the same approach that the same limiting equation arises for a related, simpler process with direct sensing of the chemoattractant gradient. Additionally, we propose a time discretization technique that retains these diffusion limits exactly, i.e., without error that depends on the time discretization. In the companion paper \cite{variance}, these results are used to construct a coupling technique that allows numerical simulation of the process with internal state with asymptotic variance reduction, in the sense that the variance vanishes in the diffusion limit.

Mathias Rousset

This paper considers Schrödinger operators, and presents a probabilistic interpretation of the variation (or shape derivative) of the Dirichlet groundstate energy when the associated domain is perturbed. This interpretation relies on the distribution on the boundary of a stopped random process with Feynman-Kac weights. Practical computations require in addition the explicit approximation of the normal derivative of the groundstate on the boundary. We then propose to use this formulation in the case of the so-called fixed node approximation of fermion groundstates, defined by the bottom eigenelements of the Schrödinger operator of a fermionic system with Dirichlet conditions on the nodes (the set of zeros) of an initially guessed skew-symmetric function. We show that the shape derivative of the fixed node energy vanishes if and only if either (i) the distribution on the nodes of the stopped random process is symmetric; or (ii) the nodes are exactly the zeros of a skew-symmetric eigenfunction of the operator. We propose an approximation of the shape derivative of the fixed node energy that can be computed with a Monte-Carlo algorithm, which can be referred to as Nodal Monte-Carlo (NMC). The latter approximation of the shape derivative also vanishes if and only if either (i) or (ii) holds.

Petr Plechac, Mathias Rousset

We propose and analyze an implicit mass-matrix penalization (IMMP) technique which enables efficient and exact sampling of the (Boltzmann/Gibbs) canonical distribution associated to Hamiltonian systems with fast degrees of freedom (fDOFs). The penalty parameters enable arbitrary tuning of the timescale for the selected fDOFs, and the method is interpreted as an interpolation between the exact Hamiltonian dynamics and the dynamics with infinitely slow fDOFs (equivalent to geometrically corrected rigid constraints). This property translates in the associated numerical methods into a tunable trade-off between stability and dynamical modification. The penalization is based on an extended Hamiltonian with artificial constraints associated with each fDOF. By construction, the resulting dynamics is statistically exact with respect to the canonical distribution in position variables. The algorithms can be easily implemented with standard geometric integrators with algebraic constraints given by the expected fDOFs, and has no additional complexity in terms of enforcing the constraint and force evaluations. The method is demonstrated on a high dimensional system with non-convex interactions. Prescribing the macroscopic dynamical timescale, it is shown that the IMMP method increases the time-step stability region with a gain that grows linearly with the size of the system. The latter property, as well as consistency of the macroscopic dynamics of the IMMP method is proved rigorously for linear interactions. Finally, when a large stiffness parameter is introduced, the IMMP method can be tuned to be asymptotically stable, converging towards the heuristically expected Markovian effective dynamics on the slow manifold.

Patrick Héas, Frédéric Cérou, Mathias Rousset

Atmospheric motion vectors (AMVs) extracted from satellite imagery are the only wind observations with good global coverage. They are important features for feeding numerical weather prediction (NWP) models. Several Bayesian models have been proposed to estimate AMVs. Although critical for correct assimilation into NWP models, very few methods provide a thorough characterization of the estimation errors. The difficulty of estimating errors stems from the specificity of the posterior distribution, which is both very high dimensional, and highly ill-conditioned due to a singular likelihood. Motivated by this difficult inverse problem, this work studies the evaluation of the (expected) estimation errors using gradient-based Markov Chain Monte Carlo (MCMC) algorithms. The main contribution is to propose a general strategy, called here chilling, which amounts to sampling a local approximation of the posterior distribution in the neighborhood of a point estimate. From a theoretical point of view, we show that under regularity assumptions, the family of chilled posterior distributions converges in distribution as temperature decreases to an optimal Gaussian approximation at a point estimate given by the Maximum A Posteriori, also known as the Laplace approximation. Chilled sampling therefore provides access to this approximation generally out of reach in such high-dimensional nonlinear contexts. From an empirical perspective, we evaluate the proposed approach based on some quantitative Bayesian criteria. Our numerical simulations are performed on synthetic and real meteorological data. They reveal that not only the proposed chilling exhibits a significant gain in terms of accuracy of the point estimates and of their associated expected errors, but also a substantial acceleration in the convergence speed of the MCMC algorithms.

Petr Plechac, Mathias Rousset

An implicit mass-matrix penalization (IMMP) of Hamiltonian dynamics is proposed, and associated dynamical integrators, as well as sampling Monte-Carlo schemes, are analyzed for systems with multiple time scales. The penalization is based on an extended Hamiltonian with artificial constraints associated with some selected DOFs. The penalty parameters enable arbitrary tuning of timescales for the selected DOFs. The IMMP dynamics is shown to be an interpolation between the exact Hamiltonian dynamics and the dynamics with rigid constraints. This property translates in the associated numerical integrator into a tunable trade-off between stability and dynamical modification. Moreover, a penalty that vanishes with the time-step yields order two convergent schemes for the exact dynamics. Moreover, by construction, the resulting dynamics preserves the canonical equilibrium distribution in position variables, up to a computable geometric correcting potential, leading to Metropolis-like unbiased sampling algorithms. The algorithms can be implemented with a simple modification of standard geometric integrators with algebraic constraints imposed on the selected DOFs, and has no additional complexity in terms of enforcing the constraints and force evaluations. The properties of the IMMP method are demonstrated numerically on the $N$-alkane model, showing that the time-step stability region of integrators and the sampling efficiency can be increased with a gain that grows with the size of the system. This feature is mathematically analyzed for a harmonic atomic chain model. When a large stiffness parameter is introduced, the IMMP method is shown to be asymptotically stable and to converge towards the heuristically expected Markovian effective dynamics on the slow manifold.