J. M. Sanz‐Serna

Teresa Klatzer, Paul Dobson, Yoann Altmann et al.

This paper presents a new accelerated proximal Markov chain Monte Carlo methodology to perform Bayesian inference in imaging inverse problems with an underlying convex geometry. The proposed strategy takes the form of a stochastic relaxed proximal-point iteration that admits two complementary interpretations. For models that are smooth or regularised by Moreau-Yosida smoothing, the algorithm is equivalent to an implicit midpoint discretisation of an overdamped Langevin diffusion targeting the posterior distribution of interest. This discretisation is asymptotically unbiased for Gaussian targets and shown to converge in an accelerated manner for any target that is $κ$-strongly log-concave (i.e., requiring in the order of $\sqrtκ$ iterations to converge, similarly to accelerated optimisation schemes), comparing favorably to [M. Pereyra, L. Vargas Mieles, K.C. Zygalakis, SIAM J. Imaging Sciences, 13,2 (2020), pp. 905-935] which is only provably accelerated for Gaussian targets and has bias. For models that are not smooth, the algorithm is equivalent to a Leimkuhler-Matthews discretisation of a Langevin diffusion targeting a Moreau-Yosida approximation of the posterior distribution of interest, and hence achieves a significantly lower bias than conventional unadjusted Langevin strategies based on the Euler-Maruyama discretisation. For targets that are $κ$-strongly log-concave, the provided non-asymptotic convergence analysis also identifies the optimal time step which maximizes the convergence speed. The proposed methodology is demonstrated through a range of experiments related to image deconvolution with Gaussian and Poisson noise, with assumption-driven and data-driven convex priors. Source codes for the numerical experiments of this paper are available from https://github.com/MI2G/accelerated-langevin-imla.

Ander Murua, J. M. Sanz-Serna

We study word series and extended word series, classes of formal series for the analysis of some dynamical systems and their discretizations. These series are similar to but more compact than B-series. They may be composed among themselves by means of a simple rule. While word series have appeared before in the literature, extended word series are introduced in this paper. We exemplify the use of extended word series by studying the reduction to normal form and averaging of some perturbed integrable problems. We also provide a detailed analysis of the behaviour of splitting numerical methods for those problems.

Ander Murua, J. M. Sanz-Serna

We study a model problem describing vibrational resonance by means of a high-order averaging technique based on so-called word series. With the tech- nique applied here, the tasks of constructing the averaged system and the associ- ated change of variables are divided into two parts. It is first necessary to build recursively a set of so-called word basis functions and, after that, all the required manipulations involve only scalar coefficients that are computed by means of sim- ple recursions. As distinct from the situation with other approaches, with word- series, high-order averaged systems may be derived without having to compute the associated change of variables. In the system considered here, the construction of high-order averaged systems makes it possible to obtain very precise approxima- tions to the true dynamics.

A. Murua, J. M. Sanz-Serna

In the first part of the present work we consider periodically or quasiperiodically forced systems of the form $(d/dt)x = εf(x,t ω)$, where $ε\ll 1$, $ω\in\mathbb{R}^d$ is a nonresonant vector of frequencies and $f(x,θ)$ is $2π$-periodic in each of the $d$ components of $θ$ (i.e.\ $θ\in\mathbb{T}^d$). We describe in detail a technique for explicitly finding a change of variables $x = u(X,θ;ε)$ and an (autonomous) averaged system $(d/dt) X = εF(X;ε)$ so that, formally, the solutions of the given system may be expressed in terms of the solutions of the averaged system by means of the relation $x(t) = u(X(t),tω;ε)$. Here $u$ and $F$ are found as series whose terms consist of vector-valued maps weighted by suitable scalar coefficients. The maps are easily written down by combining the Fourier coefficients of $f$ and the coefficients are found with the help of simple recursions. Furthermore these coefficients are {\em universal} in the sense that they do not depend on the particular $f$ under consideration. In the second part of the contribution, we study problems of the form $(d/dt) x = g(x)+f(x)$, where one knows how to integrate the "unperturbed" problem $(d/dt)x = g(x)$ and $f$ is a perturbation satisfying appropriate hypotheses. It is shown how to explicitly rewrite the system in the "normal form" $(d/dt) x = \bar g(x)+\bar f(x)$, where $\bar g$ and $\bar f$ are {\em commuting} vector fields and the flow of $(d/dt) x = \bar g(x)$ is conjugate to that of the unperturbed $(d/dt)x = g(x)$. In Hamiltonian problems the normal form directly leads to the explicit construction of formal invariants of motion. Again, $\bar g$, $\bar f$ and the invariants are written as series consisting of known vector-valued maps and universal scalar coefficients that may be found recursively.

A. Alamo, J. M. Sanz-Serna

We present an analysis based on word combinatorics of splitting integrators for Ito or Stratonovich systems of stochastic differential equations. In particular we present a technique to write down systematically the expansion of the local error; this makes it possible to easily formulate the conditions that guarantee that a given integrator achieves a prescribed strong or weak order. This approach bypasses the need to use the Baker-Campbell-Hausdorff (BCH) formula and shows the existence of an order barrier of two for the attainable weak order. The paper also provides a succinct introduction to the combinatorics of words.

J. M. Sanz-Serna, Beibei Zhu

We propose and analyze a heterogenous multiscale method for the efficient integration of constant-delay differential equations subject to fast periodic forcing. The stroboscopic averaging method (SAM) suggested here may provide approximations with $\(\mathcal{O}(H^2+1/Ω^2)\)$ errors with a computational effort that grows like $\(H^{-1}\)$ (the inverse of the stepsize), uniformly in the forcing frequency Omega.

M. P. Calvo, J. M. Sanz-Serna, Beibei Zhu

We introduce and analyze a family of heterogeneous multiscale methods for the numerical integration of highly oscillatory systems of delay differential equations with constant delays. The methodology suggested provides algorithms of arbitrarily high accuracy.

J. M. Sanz-Serna

This article is an abridged and commented translation into Spanish of the 1815 memoir where Gauss introduced the quadrature rules now associated with his name. Gauss' work does not resemble at all the stardard text-book treatment of Gaussian quadrature. The original memoir is an example of mathematical virtuosity, based on a superb use of series, where the problem is reformulated as a problem in functional approximation that is solved by means of continued fractions.

J. M. Sanz-Serna, Konstantinos C. Zygalakis

We present a framework that allows for the non-asymptotic study of the $2$-Wasserstein distance between the invariant distribution of an ergodic stochastic differential equation and the distribution of its numerical approximation in the strongly log-concave case. This allows us to study in a unified way a number of different integrators proposed in the literature for the overdamped and underdamped Langevin dynamics. In addition, we analyse a novel splitting method for the underdamped Langevin dynamics which only requires one gradient evaluation per time step. Under an additional smoothness assumption on a $d$--dimensional strongly log-concave distribution with condition number $κ$, the algorithm is shown to produce with an $\mathcal{O}\big(κ^{5/4} d^{1/4}ε^{-1/2} \big)$ complexity samples from a distribution that, in Wasserstein distance, is at most $ε>0$ away from the target distribution.

J. M. Sanz-Serna, Konstantinos C. Zygalakis

In this manuscript, we study the properties of a family of second-order differential equations with damping, its discretizations and their connections with accelerated optimization algorithms for $m$-strongly convex and $L$-smooth functions. In particular, using the Linear Matrix Inequality LMI framework developed by \emph{Fazlyab et. al. $(2018)$}, we derive analytically a (discrete) Lyapunov function for a two-parameter family of Nesterov optimization methods, which allows for the complete characterization of their convergence rate. In the appropriate limit, this family of methods may be seen as a discretization of a family of second-order ordinary differential equations for which we construct(continuous) Lyapunov functions by means of the LMI framework. The continuous Lyapunov functions may alternatively, be obtained by studying the limiting behaviour of their discrete counterparts. Finally, we show that the majority of typical discretizations of the family of ODEs, such as the Heavy ball method, do not possess Lyapunov functions with properties similar to those of the Lyapunov function constructed here for the Nesterov method.

A. Murua, J. M. Sanz-Serna

In a series of papers the present authors and their coworkers have developed a family of algebraic techniques to solve a number of problems in the theory of discrete or continuous dynamical systems and to analyze numerical integrators. Given a specific problem, those techniques construct an abstract, {\em universal} version of it which is solved algebraically; then, the results are tranferred to the original problem with the help of a suitable morphism. In earlier contributions, the abstract problem is formulated either in the dual of the shuffle Hopf algebra or in the dual of the Connes-Kreimer Hopf algebra. In the present contribution we extend these techniques to more general Hopf algebras, which in some cases lead to more efficient computations.

A. Alamo, J. M. Sanz-Serna

We present a technique, based on so-called word series, to write down in a systematic way expansions of the strong and weak local errors of splitting algorithms for the integration of Stratonovich stochastic differential equations. Those expansions immediately lead to the corresponding order conditions. Word series are similar to, but simpler than, the B-series used to analyze Runge-Kutta and other one-step integrators. The suggested approach makes it unnecessary to use the Baker-Campbell-Hausdorff formula. As an application, we compare two splitting algorithms recently considered by Leimkuhler and Matthews to integrate the Langevin equations. The word series method bears out clearly reasons for the advantages of one algorithm over the other.

J. M. Sanz-Serna

It is well known that symplectic Runge-Kutta and Partitioned Runge-Kutta methods exactly preserve {\em quadratic} first integrals (invariants of motion) of the system being integrated. While this property is often seen as a mere curiosity (it does not hold for arbitrary first integrals), it plays an important role in the computation of numerical sensitivities, optimal control theory and Lagrangian mechanics, as described in this paper, which, together with some new material, presents in a unified way a number of results now scattered or implicit in the literature. Some widely used procedures, such as the direct method in optimal control theory and the computation of sensitivities via reverse accumulation imply "hidden" integrations with symplectic Partitioned Runge-Kutta schemes.

Jesus Maria Sanz-Serna, Ander Murua

This paper provides a brief history of B-series and the associated Butcher group and presents the new theory of word series and extended word series. B-series (Hairer and Wanner 1976) are formal series of functions parameterized by rooted trees. They greatly simplify the study of Runge-Kutta schemes and other numerical integrators. We examine the problems that led to the introduction of B-series and survey a number of more recent developments, including applications outside numerical mathematics. Word series (series of functions parameterized by words from an alphabet) provide in some cases a very convenient alternative to B-series. Associated with word series is a group G of coefficients with a composition rule simpler than the corresponding rule in the Butcher group. From a more mathematical point of view, integrators, like Runge-Kutta schemes, that are affine equivariant are represented by elements of the Butcher group, integrators that are equivariant with respect to arbitrary changes of variables are represented by elements of the word group G.

Alexandros Beskos, Natesh S. Pillai, Gareth O. Roberts et al.

We investigate the properties of the Hybrid Monte-Carlo algorithm (HMC) in high dimensions. HMC develops a Markov chain reversible w.r.t. a given target distribution $Π$ by using separable Hamiltonian dynamics with potential $-\logΠ$. The additional momentum variables are chosen at random from the Boltzmann distribution and the continuous-time Hamiltonian dynamics are then discretised using the leapfrog scheme. The induced bias is removed via a Metropolis-Hastings accept/reject rule. In the simplified scenario of independent, identically distributed components, we prove that, to obtain an $\mathcal{O}(1)$ acceptance probability as the dimension $d$ of the state space tends to $\infty$, the leapfrog step-size $h$ should be scaled as $h= l \times d^{-1/4}$. Therefore, in high dimensions, HMC requires $\mathcal{O}(d^{1/4})$ steps to traverse the state space. We also identify analytically the asymptotically optimal acceptance probability, which turns out to be 0.651 (to three decimal places). This is the choice which optimally balances the cost of generating a proposal, which {\em decreases} as $l$ increases, against the cost related to the average number of proposals required to obtain acceptance, which {\em increases} as $l$ increases.