Gilles Vilmart

Postprocessed integrators for the high order integration of ergodic SDEs

Keywords: stochastic differential equations, effective order, postprocessor, modified differential equations, invariant measure, ergodicity AMS subject classification (2010): 65C30, 60H35, 37M25

The concept of effective order is a popular methodology in the deterministic literature for the construction of efficient and accurate integrators for differential equations over long times. The idea is to enhance the accuracy of a numerical method by using an appropriate change of variables called the processor. We show that this technique can be extended to the stochastic context for the construction of new high order integrators for the sampling of the invariant measure of ergodic systems. The approach is illustrated with modifications of the stochastic θ-method applied to Brownian dynamics, where postprocessors achieving order two are introduced. Numerical experiments, including stiff ergodic systems, illustrate the efficiency and versatility of the approach.

Introduction

A popular technique in the deterministic literature for the construction of efficient and high order integrators for differential equations is the notion of effective order of Butcher [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF]. Inspired by this idea, we introduce in this paper a new methodology for the construction of high order integrators for the invariant measure of ergodic systems of stochastic differential equations (SDEs). Given a deterministic problem

dy(t) dt = f (y(t)), (1)
with a smooth vector field f : R d → R d , the idea is to search for a computationally efficient integrator K h for (1), named the kernel, and perturbations of the identity map χ h and χ -1 h , inverse of each other, and called respectively the postprocessor and preprocessors, such that the integrator y n+1 = Φ h (y n ) defined by

Φ h = χ h • K h • χ -1 h
has a higher order of accuracy than K h . Considering a constant stepsize h yields

Φ N h = χ h • K N h • χ -1 h , (2)
which means that N compositions of Φ h can be efficiently implemented by applying N steps of the kernel K h and only once the processor maps χ -1 h , χ h , respectively at the beginning and the end of the integration. This approach is efficient if the computational cost of the kernel K h is cheaper than the cost of Φ h . It was first introduced in [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF] as a mean to construct 5th order explicit Runge-Kutta methods with only 5 internal stages, corresponding to the number of function evaluations per time step. A renewed interest for this idea appeared in the context of deterministic geometric numerical integration [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF][START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF][START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF], where the use of a constant stepsize is required in general, for the construction of high order integrators with a favourable long time behaviour (see e.g. [13, Sect. V. 3.3] and references therein).

In this paper, we consider a system of (Itô) stochastic differential equations

dX(t) = f (X(t))dt + g(X(t))dW (t), X(0) = X 0 , (3)
where f : R d → R d and g : R d → R d×m are smooth fields assumed globally Lipchitz for simplicity, and W (t) = (W 1 (t), . . . , W m (t)) T is a standard m-dimensional Wiener process. Under appropriate smoothness and growth assumptions on the fields f, g, one can show that (9) has ergodic dynamics, i.e. for an arbitrary initial condition X 0 = x, assumed deterministic for simplicity, the time average of the trajectory satisfies with probability 1, lim

T →∞ 1 T T 0 φ(X(t)) = R d φ(y)dµ(y).
for all smooth test functions φ with derivatives at all orders with polynomial growth. In addition, in many situations, one can show the following exponential ergodicity property for all t ≥ 0, E(φ(X(t))) -

R d φ(y)dµ(y) ≤ K(x)e -λt (4)
where λ > 0 is independent of x, t and K(x) depends on φ and the initial condition X 0 = x (assume deterministic for simplicity), but is independent of t. Consider a one-step numerical integrator for [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF],

X n+1 = Ψ(X n , h, ξ n ). (5)
We say that numerical method X n for the approximation of (9) at time t n = nh is ergodic if it has a unique invariant probability law µ h with finite moments of any order and with probability 1, lim

N →∞ 1 N + 1 N n=0 φ(X n ) = R d φ(y)dµ h (y), (6)
for all initial condition X 0 = x and all test functions φ. We say that the scheme has order r ≥ 1 with respect to the invariant measure if |e(φ, h)| ≤ Ch r with e(φ, h) := lim

N →∞ 1 N + 1 N n=0 φ(X n ) - R d φ(y)dµ(y), (7)
where C is independent of h small enough and X 0 . Analogously to (4), in many situations, one can show for all t n = nh and all initial condition X 0 = x, the estimate

E(φ(X n )) - R d φ(y)dµ(y) ≤ K(x)e -λtn + Ch r (8)
where C, K(x) are independent of n, h.

A standard approach for the construction of integrators with high order for the invariant measure is to design methods with high weak order of accuracy. This is known for large classes of ergodic SDEs in particular for a degenerate noise and locally Lipschitz fields in [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. This can be combined with extrapolation techniques from Talay [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]. It is known, however, that high order can still be achieved using low weak order methods, as shown in particular for certain classes of splitting methods for nonlinear Langevin dynamics [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF][START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF][START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for second order stochastic dynamics[END_REF]. We also mention the exact schemes based on Markov Chain Monte-Carlo methods, see e.g. the survey [START_REF] Sanz-Serna | Chain Monte Carlo and Numerical Differential Equation[END_REF]. An approach proposed in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF][START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for second order stochastic dynamics[END_REF] permits to design methods with a high order for the invariant measure despite a low deterministic order of accuracy. It relies on recent developments in the stochastic context [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped langevin processes[END_REF] of the theory of backward error analysis and modified differential equations, a fundamental tool in the field of deterministic geometric numerical integration [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF][START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF][START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. In the same spirit, the proposed approach of postprocessed integrators, based on Butcher's effective order [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF], relies on many ideas from [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF].

We shall illustrate our approach at the example of Brownian dynamics generated by the overdamped Langevin equation

dX(t) = -∇V (X(t))dt + σdW (t), (9)
where

V : R d → R is a smooth potential, σ > 0 is a constant, and W = (W 1 , . . . , W d ) T is a standard d-dimensional
Wiener process. This model describes the motion of a particle in a potential subject to thermal noise [START_REF] Risken | The Fokker-Planck equation[END_REF][START_REF] Gardiner | Handbook of stochastic methods[END_REF]. Under appropriate assumptions on the potential V , the system is ergodic and has a unique invariant measure µ which is characterized by the Gibbs density function (sometimes called the Boltzmann density function)

ρ(x) = Ze -2 σ 2 V (x) , ( 10
)
where Z is a renormalization constant such that R d ρ(y)dy = 1. We show that our approach is closely related to an efficient non-Markovian method of order two proposed in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] for [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], with convergence properties first analyzed in [START_REF] Leimkuhler | On the long-time integration of stochastic gradient systems[END_REF]. It was obtained considering the large friction limit in a particular splitting method for Langevin dynamics. We show in this paper that this scheme can be interpreted as Markovian up to a change of variable, which is precisely the postprocessor map. This paper is organized as follows. In Section 2, we introduce our new methodology to extend Butcher's effective order technique to the stochastic context and we provide examples of postprocessed integrators. In Section 3, we recall standard assumptions that guaranty the ergodicity of the system and we state the assumptions on the numerical methods needed for our analysis. In Section 4, we state our main results and derive conditions on a postprocessor to achieve high order for the invariant measure of a given integrator with emphasis on the case of Brownian dynamics. Finally, the Section 5 is devoted to the numerical experiments, and shows the efficiency and versatility of the approach for various stiff and nonstiff Brownian dynamics models.

Postprocessed integrators

In this section, we describe the idea of postprocessed integrators for an ergodic system of SDEs [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF]. Given a smooth test function φ, we recall that there are two standard strategies when computing ergodic averages R d φ(y)dµ(y) using a one step integrator [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] to compute numerical realizations of (3).

1. Given an initial condition X 0 , one can consider a single numerical trajectory X 0 , X 1 , X 2 , . . . over a long time interval of size T = N h, and then compute the approximation

1 N + 1 N k=0 φ(X k ) ≃ R d φ(y)dµ(y),
taking advantage of the estimate (7).

2. Alternatively, one can compute a large number M of independent numerical trajectories {X i 0 , X i 1 , X i 2 , . . . , X i N }, i = 1 . . . , M on a relatively shorter time interval T = N h and use the approximation

1 M M i=1 φ(X i N ) ≃ R d φ(y)dµ(y),
based on the exponential estimate [START_REF] Burrage | Low rank Runge-Kutta methods, symplecticity and stochastic Hamiltonian problems with additive noise[END_REF]. The advantage of this second approach, as shown in [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF] in the context of ergodic SDEs, is that it permits to treat easily the case of non-globally Lipchitz vector fields using standard integrators (in particular explicit ones) by using the concept of rejecting exploding trajectories, a rigorous methodology introduced in [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally lipschitz coefficients[END_REF]. This procedure is often preceded by an equilibration period to make the law of X i 0 close to the invariant measure µ, see e.g. [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF] for details.

The idea of postprocessed integrators is to consider an integrator [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] with order r for the invariant measure in [START_REF] Burrage | Numerical methods for strong solutions of stochastic differential equations: an overview[END_REF] and [START_REF] Burrage | Low rank Runge-Kutta methods, symplecticity and stochastic Hamiltonian problems with additive noise[END_REF], and to search for a postprocessor X n = G n (X n ), which can be a deterministic or a stochastic perturbation of X n , such that X n has a higher order of accuracy r > r for the invariant measure. Precisely, for all h assumed small enough, and all t n = nh,

lim N →∞ 1 N + 1 N k=0 φ(X k ) - R d φ(y)dµ(y) ≤ Ch r (11)
E(φ(X n )) - R d φ(y)dµ(y) ≤ K(x)e -λtn + Ch r , (12)
where C, λ are independent of h, n and the initial condition X 0 = x, and K(x) is independent of h, n.

If ergodic averages are computed using the estimate [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF], then the postprocessor X k = G k (X k ) is applied at every time step k, but it has a negligible overcost in some situations where no additional function evaluations are needed as shown below. In addition X k can always be computed in parallel with X k+1 , independently. If alternatively, ergodic averages are approximated based on the exponential estimate [START_REF] Gardiner | Handbook of stochastic methods[END_REF], then the postprocessor X N = G N (X N ) is applied only once at the end time t N = N h of each numerical trajectory, resulting in a negligible overhead.

Remark 2.1. Compared to the deterministic setting (2), G n plays the role of the postprocessor map χ h . In contrast, notice however that no preprocessor χ -1 h is involved in the stochastic setting. This is not surprising because the limit E(φ(X(t))) → R d φ(y)µ(y)dy for t → ∞ is independent of the initial condition X 0 = x. In other words, the ergodic systems "forgets" the initial condition. The advantage compared to the deterministic case is that the postprocessor map G n needs not to be reversible, it can be in particular the flow of a stiff deterministic problem, or even a stochastic problem.

To illustrate our new approach, we focus on the problem of Brownian dynamics [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], and we consider the standard stochastic θ-method

X n+1 = X n + h(1 -θ)f (X n ) + hθf (X n+1 ) + σ √ hξ n , (13)
where f = -∇V and ξ n ∼ N (0, I) are independent Gaussian random vectors with dimension d. For θ = 0, it reduces to the Euler-Maruyama method, while for θ = 0, the method is implicit and requires in general the resolution of a non-linear system at each timestep. For θ ≥ 1/2, it is shown in [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] that the θ-method is mean-square A-stable, a property which makes it well-suited for stiff mean-square stable stochastic problems. This scheme has weak order one of accuracy, except in the special case θ = 1/2 for which the weak order of the method is two. Is was shown in [START_REF] Schurz | Preservation of probabilistic laws through Euler methods for Ornstein-Uhlenbeck process[END_REF] that the θ-method preserves exactly the invariant measure in the linear case if and only if θ = 1/2 (i.e. for V (x) a quadratic function in [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], which corresponds to an Orstein-Uhlenbeck process). This result was extended in [START_REF] Burrage | Low rank Runge-Kutta methods, symplecticity and stochastic Hamiltonian problems with additive noise[END_REF] for linear second order problems with additive noise, where it was shown that the θ = 1/2-method is in fact the only measure exact scheme in a certain class of stochastic Runge-Kutta methods involving only one Wiener increment per step.

For multi-dimensional nonlinear problems of the form (9), we introduce the following postprocessed methods of order 2 for the invariant measure based on modifications of the θ-method [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. Precisely, we prove in Section 4 that for each method below, X n has order r = 2 for the invariant measure in [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF], [START_REF] Gardiner | Handbook of stochastic methods[END_REF]. For θ = 0, we obtain the following modification of the Euler-Maruyama method, which is equivalent to a method first proposed in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] (see [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF] below)

X n+1 = X n + hf X n + 1 2 σ √ hξ n + σ √ hξ n , X n = X n + 1 2 σ √ hξ n . (14)
For θ = 1, we derive the new scheme

X n+1 = X n + hf (X n+1 + aσ √ hξ n ) + σ √ hξ n , a = - 1 2 + √ 2 2 , (15)
together with two possible postprocessors

X n = X n + cσ √ hJ -1 n ξ n , c = 2 √ 2 -1 2 (16)
and alternatively,

X n = X n + hbf (X n ) + cσ √ hξ n , b = √ 2 2 , c = 4 √ 2 -1 2 . ( 17
)
The matrix J -1 n in ( 16) is defined as the inverse of the d × d Jacobian matrix

J n = I -hf ′ (X n + aσ √ hξ n-1 ).
Remark 2.2. The matrix J n is involved in standard Newton iterations needed to solve the nonlinear system [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF]. For an efficient implementation, it is known that one should never compute the inverse J -1 n itself but rather solve the appropriate linear systems. Using a precomputed LU decomposition of this matrix, the evaluation of J -1 n ξ n in (16) therefore has a negligible cost. This term is introduced only to improve the numerical stability and does not influence the order 2 of accuracy of the postprocessed method (15),( 16), as discussed in Section 4 (Theorem 4.6). The use of Jacobian matrices such as J n for stabilization is classical for stiff deterministic problems [14, Sect. IV.8], and reveals efficient also in the context of mean-square stable stiff SDEs in [START_REF] Abdulle | Mean-square A-stable diagonally driftimplicit integrators with high order for stiff Ito systems of stochastic differential equations with noncommutative noise[END_REF]. 14), a calculation shows that the scheme ( 14) is equivalent to the method

Substituting X n = X n -1 2 σ √ hξ n in (
X n+1 = X n + hf (X n ) + 1 2 σ √ h(ξ n + ξ n+1 ). ( 18
)
The scheme (18) of order two, which has the same cost as the standard Euler-Maruyama method of order one, was first introduced in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] as a non-Markovian method. A rigorous analysis of its order two of convergence for the numerical invariant measure of ( 9) was first proposed in [START_REF] Leimkuhler | On the long-time integration of stochastic gradient systems[END_REF]. Notice that the point of view of postprocessed integrators provides here an alternative proof of the order two of accuracy of the scheme [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF]. Analogously, notice that a non-Markovian formulation of the postprocessed method ( 15),( 16) can also be written,

X n+1 = X n + hf (X n+1 + σ √ h(aξ n -cξ n+1 )) + σ √ h((1 -c)ξ n + cξ n+1 ),
where a, c are defined in ( 15),( 16), and we omit for simplicity the stabilization term J n .

Background

We recall in this section standard results and assumptions that will be useful for the analysis of postprocessed integrators in Section 4. We assume that f, g in ( 3) are smooth fields and we focus for simplicity on the case of globally Lipchitz vector fields.

We denote C ∞ P (R d , R) the set of C ∞ functions whose derivatives up to any order have a polynomial growth of the form

|φ(x)| ≤ C(1 + |x| s ) (19)
for some constants s and C independent of x. A classical fundamental tool for the analysis of

(3) is the backward Kolmogorov equation [27, Chap. 2]. Given φ ∈ C ∞ P (R d , R) and setting u(t, x) = E (φ(X(t))|X 0 = x), it yields that u(x, t) is the solution of the deterministic parabolic PDE ∂u ∂t = Lu, u(x, 0) = φ(x), x ∈ R d , t > 0, ( 20
)
where the generator L is defined by

Lφ := f • ∇φ + 1 2 gg T : ∇ 2 φ, (21)
where ∇φ and ∇ 2 φ are the gradient vector and the Hessian matrix of φ, and A : B = trace(A T B) is a scalar product on matrices.

Assumption 3.1. We assume the following.

1. f, g are of class C ∞ , with bounded derivatives of any order, and g is bounded;

2. the generator L in (21) is a uniformly elliptic operator, i.e. there exists α > 0 such that for all x, ξ ∈ R d , x T g(ξ)g(ξ) T x ≥ αx T x;

3. there exist C, β > 0 such that for all

x ∈ R d , x T f (x) ≤ -βx T x + C.
Assumption 3.1 implies automatically that (3) is ergodic with a unique invariant measure (see e.g. [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF]) and satisfies the exponential ergodicity property [START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for second order stochastic dynamics[END_REF]. Notice that hypo-elliptic generators L could also be considered, to include the case of Langevin dynamics. We obtain for all

φ ∈ C ∞ P (R d , R) the identity lim t→∞ u(x, t) = φ(x) + ∞ 0 Lu(t, x)dt = R d φ(y)ρ(y)dy. ( 22
)
Considering L 2 -adjoint of the generator L, given by

L * φ = -div(f φ) + 1 2 gg T : ∇ 2 φ, (23)
we recall that the density ρ of the invariant measure µ of ( 3) is the unique solution of

L * ρ = 0. ( 24
)
We make the following assumptions on the one step integrator [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF]. For

φ ∈ C ∞ P (R d , R), consider the quantity U (x, h) = E(φ(X 1 )|X 0 = x)). ( 25
)
We assume that the scheme has local weak order p i.e. it satisfies for all initial condition

X 0 = x, |E(φ(X 1 )) -E(φ(X(h)))| ≤ C(x)h p+1 , (26)
for all h sufficiently small, where C(x) is independent of h but depends on x, φ.

Assumption 3.2. We assume that (25) has a weak Taylor series expansion of the form,

U (x, h) = φ(x) + hA 0 φ(x) + h 2 A 1 φ(x) + . . . , (27)
for all φ ∈ C ∞ P (R d , R), where A i , i = 0, 1, 2, . . . are linear differential operators with coefficients depending smoothly on the drift and diffusion functions f, g, and their derivatives (and depending on the choice of the integrator). In addition, we assume that A 0 coincides with the generator L given in [START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF], which means that the method has (at least) local order one in the weak sense,

A 0 (f, g) = L. ( 28
)
For instance, for the θ-method (13) applied to (9), a calculation shows

A 1 φ = 1 2 φ ′′ (f, f ) + σ 2 2 d i=1 φ ′′′ (e i , e i , f ) + σ 4 8 d i,j=1 φ (4) (e i , e i , e j , e j ) + θφ ′ f ′ f + θ σ 2 2 φ ′ d i=1 f ′′ (e i , e i ) + θσ 2 d i=1 φ ′′ (f ′ e i , e i ), (29)
where e 1 , . . . , e d is the canonical basis of R d , and φ ′ (x), φ ′′ (x), φ ′′′ (x), . . . are the derivatives of φ, which are linear, bilinear, trilinear, . . . , forms, respectively (and analogously for f ). We refer to [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF] for examples of similar calculations. We shall also use the following assumption [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2.2]) which permits to guaranty that the numerical moments E(|X n | k ) up to any order k are bounded along time. This also permits to infer automatically "a global weak order p" from the local weak error (26) of a given integrator, see [27, Chap. 2.2].

Assumption 3.3. We assume that the numerical integrator (3) satisfies for all

x ∈ R d that |E(X 1 -X 0 |X 0 = x)| ≤ C(1 + |x|)h, |X 1 -X 0 | ≤ M (1 + |X 0 |) √ h, ( 30
)
where C is independent of h small enough and M is a random variable that has bounded moments of all orders independent of h and X 0 .

Theorem 3.4 below is standard and combines results derived by Talay and Milstein. Precisely, it has been proved in [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] for specific methods (the Euler-Maruyama and the Milstein methods), while the general procedure to infer the global weak order from the local weak order is due to Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] and can be found in [27, Chap. 2.2, 2.3]. The proof of Theorem 3.4 is thus omitted. Theorem 3.4. [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF][START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] Let X N be an ergodic numerical solution of (3) on [0, T ] and assume Assumption 3.1, Assumption 3.2, Assumption 3.3. Assume further the local weak order p estimate [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] where C(x) has a polynomial growth of the form [START_REF] Kopec | Weak backward error analysis for overdamped langevin processes[END_REF]. Then, the invariant measure error in [START_REF] Burrage | Numerical methods for strong solutions of stochastic differential equations: an overview[END_REF]

satisfies for all φ ∈ C ∞ P (R d , R) and h → 0, e(φ, h) = λ p h p + O(h p+1 ) (31)
for any deterministic initial condition, with λ p defined as

λ p = +∞ 0 R d A p - 1 (p + 1)! L p+1 u(y, t)ρ(y)dydt (32)
where u(x, t) is the solution of (20). In addition, there exists constants λ, C, K(x) > 0 such that for all n ≥ 0,

E(φ(X n )) - R d φ(y)dµ(y) -λ p h p ≤ K(x)e -λtn + Ch p+1 (33)
where t n = nh, the constants C, K(x), λ are independent of n, h, and h is assumed small enough (C, K(x) depend on φ and K(x) depends on x).

In fact, the assumption that the integrator has weak order p in Theorem 3.4 can be relaxed, as explain in the following remark. Remark 3.5. Using L * ρ = 0 and considering the L 2 -adjoint 1 A * p of the linear differential operator A p , we deduce that λ p in (32) satisfies

λ p = +∞ 0 R d u(y, t)A * p ρ(y)dydt.
1 For instance, for the Euler-Maruyama method, setting θ = 0 in (29) and using the notation ∆ for the Laplace operator, we obtain

A1φ = 1 2 φ ′′ (f, f ) + σ 2 2 (∆φ) ′ f + σ 4 8 ∆ 2 φ and A * 1 φ = 1 2 div φ(f ′ f + div(f )f ) + f 2 ∇φ -σ 2 f ∆φ + σ 4 8 ∆ 2 φ.
We see that if A * p ρ = 0, then the scheme has order p + 1 for the invariant measure. More generally, it is shown in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] based on backward error analysis results in [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF] on the torus and generalizations [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped langevin processes[END_REF] in R d that the (local) weak order p assumption (26) in Theorem 3.4 (equivalent to A j+1 = L j /j! for all j < p in Assumption 3.2) can be replaced by the more general condition

A * j ρ = 0, for j = 1, • • • p -1. ( 34
)
This permits to derive in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] new high order integrators based on modified differential equations while it serves in [START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for second order stochastic dynamics[END_REF] as a tool for the analysis of the order of convergence of Lie-Trotter splitting methods applied to Langevin dynamics.

Analysis of postprocessed integrators

This section is devoted to the analysis of postprocessed integrators. To illustrate our approach, we first study the simplest linear case of the Orstein-Uhlenbeck process in dimension one. We then focus on the analysis for nonlinear problems, with special emphasis on the case of Brownian dynamics.

Accuracy and stability analysis in the linear case

We consider for simplicity the Orstein-Uhlenbeck process in dimensions

d = m = 1, dX = -γXdt + σdW (t), (35)
where γ, σ > 0 are fixed constants. It corresponds to Brownian dynamics (9) with V (x) = γ 2 x 2 . For all initial condition X 0 = x, the exact solution X(t) is a Gaussian random variable with lim t→∞ E(X(t) 2 ) = σ 2 2γ . Indeed, this system is ergodic and its invariant measure is Gaussian, with density

ρ(x) = 1 σ √ 2π e -γ 2σ 2 x 2 .
Since the invariant measure is Gaussian, the error is obtained simply by comparing the exact and numerical variances. Applying a one step numerical method to [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF], in particular a Runge-Kutta method, yields in general a recursion of the form

X n+1 = A(z)X n + B(z) √ hσξ n , (36)
where z = -γh and A(z) = 1 + z + O(z 2 ) and B(z) = 1 + O(z) are analytic functions depending on the choice of the method and ξ n ∼ N (0, 1) are independent Gaussian random variables. Then, a calculation using the recursion

E(X 2 n+1 ) = E(X 2 n )A(z) 2 + B(z) 2 σ 2 h yields lim n→∞ E(X 2 n ) = σ 2 2γ |R(z)|, where R(z) = 2zB(z) 2 A(z) 2 -1 , z = -γh, (37)
for all h small enough such that |A(z)| < 1 (stable stepsize), and

lim n→∞ E(X 2 n ) = ∞ if |A(z)| > 1 and E(X 2 0 ) = 0 (unstable stepsize).
Stability for linear problems We recall here useful standard stability results. The scalar linear SDE [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF] can be considered as a test equation which gives insight on the behavior for nonlinear problem with additive noise. Notice that other test equations are used in the literature in the case of multiplicative noise, see [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF][START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF][START_REF] Saito | Mean-square stability of numerical schemes for stochastic differential systems[END_REF][START_REF] Burrage | Numerical methods for strong solutions of stochastic differential equations: an overview[END_REF][START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF][START_REF] Rathinasamy | Mean-square stability of second-order Runge-Kutta methods for multi-dimensional linear stochastic differential systems[END_REF][START_REF] Buckwar | Towards a systematic linear stability analysis of numerical methods for systems of stochastic differential equations[END_REF]. We see from the above discussion that E(X 2 n ) remains bounded as n → ∞ (mean-square stability) if |A(z)| < 1. For the θ method, it can be checked [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF][START_REF] Higham | A-stability and stochastic mean-square stability[END_REF] that this is satisfied for all stepsize h > 0 if and only if θ ≥ 1/2, while for the Euler-Maruyama method (θ = 0), we have the stepsize restriction hγ < 2 which can be a severe restriction for large values of γ. This stability property for all complex γ with positive real part, unconditionally on the stepsize h, is called A-stability in the framework of deterministic Runge-Kutta methods [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]. This makes the θ-method with θ ≥ 1/2 well suited for stiff stochastic problems, in contrast to the Euler-Maruyama method. We have the additional feature lim γ→+∞ E(X 2 1 ) = 0 for all X 0 = x if and only if θ = 1, a property called L-stability in the setting of deterministic Runge-Kutta-methods (σ = 0), and that is also important in the context of mean-square stable stiff SDEs [START_REF] Abdulle | Mean-square A-stable diagonally driftimplicit integrators with high order for stiff Ito systems of stochastic differential equations with noncommutative noise[END_REF]. The L-stability property is desirable for severely stiff problems with very large eigenvalues on the negative real axis, because it permits to damp the higher modes. We emphasis that the L-stability property lim γ→+∞ E(X 2 1 ) = 0 is also satisfied for the postprocessed methods ( 15),( 16) and ( 15), [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF].

Accuracy analysis We see from [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] that the method (without postprocessor) has order p for the Gaussian invariant measure if and only if

R(z) = 1 + O(z p ) for z → 0. Applying a corrector X n = G n (X n ) of the form G n (x) = C(z)x + D(z) √ hσξ n , z = -γh,
where C(z) = 1 + O(z) does not affect the mean-square stability of X n , we obtain

lim n→∞ E(X 2 n ) = σ 2 2γ C(z) 2 |R(z)| + hσ 2 D(z) 2 , z = -γh.
Comparing with the exact solution for which lim t→∞ E(X(t) 2 ) = σ 2 2γ , we deduce that the method has order p for the invariant measure of the Orstein-Uhlenbeck process if and only if the above limit is σ 2 2γ + O(h p ), equivalently

R(z)C(z) 2 -2zD(z) 2 = 1 + O(z p ), z → 0. ( 38
)
It turns out that the above identity holds exactly for the postprocessed scheme [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF] for which

A(z) = 1 + z, B(z) = 1 + z/2, C(z) = 1, D(z) = 1/2
. This yields the following remark which states that the explicit postprocessed scheme (14) (equivalently ( 18) from [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF]) has no bias in sampling the invariant measure for all stepsize h, a feature in the linear case shared with the implicit θ = 1/2-method [START_REF] Schurz | Preservation of probabilistic laws through Euler methods for Ornstein-Uhlenbeck process[END_REF].

Remark 4.1. The postprocessed scheme (14) applied to Brownian dynamics (9) samples exactly the invariant measure in the case of a linear drift (i.e. if V (x) is a quadratic function on R d ).

In [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], in the context of linear second order problems, new s-stage explicit Runge-Kutta methods are introduced with high order in the velocity correlation matrix with accuracy O(h s ) for s = 3, 4, 5, respectively. In the same spirit, we make the following remark which states using [START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] that there always exist postprocessors achieving arbitrarily high order for sampling the invariant measure for all consistent integrators applied to the Orstein-Uhlenbeck process. Remark 4.2. Given a consistent method of the form (36) for [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF], for arbitrary integer s, there exists a deterministic postprocessor of the form X n = C(-γh)X n where C(z) = R(z) -1/2 + O(z s ) is a polynomial of degree strictly less than s (corresponding in particular to an (s -1)-stage explicit Runge-Kutta method) such that the corresponding postprocessed scheme has order s for the invariant measure. For instance, for the Euler-Maruyama method (see [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] with θ = 0), we obtain C(z) = 1 + z/2 = 1 + z/4z 2 /32 + . . ..

Postprocessed integrators for nonlinear ergodic SDEs

We may now state the main result of this paper for the analysis of postprocessed integrators. We shall use the commutator notation [A, B] = AB -BA on linear differential operators. Theorem 4.3. Assume the hypotheses of Theorem 3.4 and consider a numerical method {X n } n≥0 with weak order p. Let G n denote independent and identically distributed random maps in R d , independent of {X j } j≤n and satisfying (30) with X 1 replaced by G 1 (X 0 ) and a weak Taylor expansion of the form

E(φ(G n (x))) = φ(x) + h p A p φ(x) + O(h p+1 ), ( 39
)
for all φ ∈ C ∞ P (R d , R)
, where the constant in O has a polynomial growth with respect to x of the form [START_REF] Kopec | Weak backward error analysis for overdamped langevin processes[END_REF]. Assuming further

(A p + [L, A p ]) * ρ = 0, ( 40
)
then X n = G n (X n ) yields an approximation of order p + 1 for the invariant measure, precisely

lim N →∞ 1 N + 1 N k=0 φ(X k ) - R d φ(y)dµ(y) ≤ Ch p+1 (41) E(φ(X n )) - R d φ(y)dµ(y) ≤ K(x)e -λtn + Ch p+1 , (42)
for all φ ∈ C ∞ P (R d , R) with t n = nh, where C, K(x) are independent of n and h assumed small enough but depend of φ (and K(x) depends on the initial condition X 0 = x).

Proof. Using ( 39) yields

E R d φ(G n (y))ρ(y)dy = R d φ(y)ρ(y)dy + h p R d φ(y)A * p ρ(y)dy + O(h p+1 ),
and we obtain from [START_REF] Leimkuhler | On the long-time integration of stochastic gradient systems[END_REF],

R d φ(y)A * p ρ(y)dy = - ∞ 0 R d Lu(t, y)A * p ρ(y)dydt = ∞ 0 R d [L, A p ]u(t, y)ρ(y)dydt = - R d A p φ(y)ρ(y)dy = -λ p
where λ p is defined in [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF], and we used assumption (40) and L * ρ = 0 (see ( 24)). We deduce

E R d φ(G n (y))ρ(y)dy = R d φ(y)ρ(y)dy -h p λ p + O(h p+1 ), (43)
We now apply Theorem 3.4 where φ is replaced by φ • G n = φ + O(h). We deduce using [START_REF] Saito | Mean-square stability of numerical schemes for stochastic differential systems[END_REF] and properties of conditional expectations that

X n = G n (X n ) satisfies E(φ(X n )) -E R d φ(G n (y))ρ(y)dy -h p λ p ≤ K(x)e -λtn + Ch p+1 .
This combined with (43) concludes the proof of (42) and analogously (41).

Remark 4.4. Similarly to Remark 3.5, observe in Theorem 4.3 that the (local) weak order p assumption (26) on the scheme X n can be replaced by the more general condition (34).

Remark 4.5. Notice that in the assumptions of Theorem 4.3, the postprocessor map G n needs not neither to involve the same random variables as X n → X n+1 nor to be independent from these random variables. In particular, the postprocessors X n in ( 14),( 16),( 17) could be defined using other independent Gaussian random variables with law N (0, I) instead of ξ n . Our numerical investigations indicate that there is no benefit to generate additional random variables.

Application to Brownian dynamics

To illustrate the proposed approach, we focus on the Brownian dynamics model ( 9) where f = -∇V is the gradient of a smooth potential. We assume2 Assumption 3.1, which guaranties the ergodicity of the system and the exponential estimate (4). We now consider the following modification with negligible overhead of the θ-method,

X n+1 = X n + h(1 -θ)f (X n + aσ √ hξ n ) + hθf (X n+1 + aσ √ hξ n ) + σ √ hξ n , ( 44
)
where a is a fixed scalar parameter, together with the postprocessor

X n = X n + bh((1 -θ)f (X n ) + θf (X n )) + cσ √ hξ n , (45)
where ξ n ∼ N (0, I) are independent Gaussian random vectors with dimension d, and b, c are fixed parameters. Notice that for all a, the modified θ-method (44) with θ = 1/2 has the same order one of accuracy for the invariant measure as the standard version [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. For a = b = c = 0, we recover the standard θ-method (13) and X n = X n . The following Theorem provides order conditions on a, b, c for X n to achieve second order of accuracy for the invariant measure.

Theorem 4.6. Assume Assumption 3.1 and consider the scheme (44) and the postprocessor (45) applied to Brownian dynamics [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF]. If

a 2 + 2θa + b -c 2 = 0, a + b -c 2 - 1 4 + θ 2 = 0, (46)
then, assuming the ergodicity of {X n }, we have lim

N →∞ 1 N + 1 N k=0 φ(X k ) - R d φ(y)dµ(y) ≤ Ch 2 E(φ(X n )) - R d φ(y)dµ(y) ≤ K(x)e -λtn + Ch 2 ,
for all φ ∈ C ∞ (R d , R) with t n = nh, where C, K(x) are independent of n and h assumed small enough. In particular, assuming ergodicity, the postprocessed schemes ( 14), ( 15)-( 16), and (15)-( 17) satisfy the above estimates.

The proof of Theorem 4.6 is postponed to the Appendix. We explain here why it yields the order two of accuracy for the invariant measure of the postprocessed schemes ( 14), ( 15)-( 16), and ( 15)- [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF].

If θ = 0 or θ ≥ 1/2, and for a given c, the order conditions (46) can be shown to have real solutions given by

a = 1 2 -θ ± 2θ(2θ -1) 2 , b = c 2 - 1 4 + θ 2 ∓ 2θ(2θ -1) 2 .
We now focus on special values of θ to construct various methods of order two for the invariant measure of (9).

• For θ = 0 (Euler-Maruyama method), we obtain a = 1 4 , b = c 2 -1 4 . Imposing respectively b = 0 and c = 0 yields the following postprocessed methods of order two: the method ( 14) from [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF],

X n+1 = X n + hf (X n + 1 2 σ √ hξ n ) + σ √ hξ n , X n = X n + 1 2 σ √ hξ n .
and the variant

X n+1 = X n + hf (X n + 1 2 σ √ hξ n ) + σ √ hξ n , X n = X n + h 4 f (X n ).
• For θ = 1 2 , we obtain a = 0 and b = c 2 . This is not surprising because the standard θ-method (13) with θ = 1 2 is already exact for sampling the invariant measure as noted in [START_REF] Schurz | Preservation of probabilistic laws through Euler methods for Ornstein-Uhlenbeck process[END_REF], and it has weak order two when applied to SDEs with additive noise such as [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF]. Postprocessors for order two are therefore useless in this case.

• For θ = 1, the assumption b = 0 yields the method (15),( 16) (without the matrix

J -1 n ) of the form X n+1 = X n + hf (X n+1 + aσ √ hξ n ) + σ √ hξ n , X n = X n + cσ √ hξ n .
Since the matrix J -1 n in ( 16) is a O(h) perturbation of the identity, its does not change the order p = 1 weak expansion of the postprocessor in assumption (39) of Theorem 4.3. The method (15),( 16) therefore has order two of accuracy for the invariant measure.

Imposing alternatively the condition of order 3 for linear problems (i.e. [START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] with p = 3) yields the method (15), [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] of the form

X n+1 = X n + hf (X n+1 + aσ √ hξ n ) + σ √ hξ n , X n = X n + hbf (X n ) + cσ √ hξ n ,
which has order two of accuracy for the invariant measure, and order three for linear problems (i.e. for a quadratic potential V ). Figure 1: Linear problem with potential (47). Implicit Euler method ((13) with θ = 1, order 1), θ = 1/2-method (order 2), and postprocessed versions (( 15),( 16) and ( 15),( 17)) of order 2. The variant of the method ( 15),( 16) without the stabilization term J -1 n is also included. Error for the second moment R (x 2 1 + x 2 2 )ρ(x)dx versus time stepsize h obtained using 10 trajectories on a long time interval of length T = 10 5 . The vertical bars indicate the standard deviation intervals due to the Monte-Carlo error.

Numerical experiments

In this section, we illustrate numerically the convergence properties of the new implicit schemes introduced in Section 2 for the high order sampling of the invariant measure of Brownian dynamics [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF]. Although our analysis applies only to globally Lipschitz vector fields, we shall consider numerical examples with non-globally Lipschitz fields, since such a regularity often arises in applications. We highlight once again that the use of the proposed schemes can still be made rigorous in this context, using the concept of rejecting exploding trajectories proposed in [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally lipschitz coefficients[END_REF][START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF].

We focus on stiff problems for which implicit methods are needed to avoid the severe step size restrictions of standard explicit integrators. We consider the standard implicit Euler method ((13) with θ = 1), the standard θ = 1/2 method (13), and the new postprocessed methods [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF],( 16) and ( 15), [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF]. We have implemented all methods in Fortran, using a Newton method with tolerance 10 -13 to solve the nonlinear implicitness of the methods. The results for the explicit postprocessed scheme [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF] from [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] are not included because a severe stepsize restriction h ≪ 1 makes explicit methods unusable for such stiff problems (analogously to the deterministic case).

Linear test problem

We consider first the linear problem in 2D, with a quadratic potential defined by

V (x) = x 2 1 2 + x 2 2 2ε + x 1 x 2 4 . ( 47
)
For the arbitrary initial condition x 1 = x 2 = -3, we plot in Figure 1 the error for the second moment R (x 2 1 + x 2 2 )ρ(x)dx for many time stepsizes h, using the average over 10 long trajectories of length 10 5 . We compare in Figure 1 the nonstiff case (ε = 1) and the stiff case (ε = 10 -6 ). [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF] for the nonlinear problems with potential (48) (nonstiff case, left pictures) and (49) (stiff case, right pictures). Four level curves are represented in solid lines, the three extrema are represented with crosses. In the bottom pictures, a numerical trajectory of the method ( 15) is included, with length T = 10 3 and stepsize h = 0.05.

• In the nonstiff case (left picture of Figure 1), we observe the expected order 2 for the postprocessed method ( 15),( 16) (we also include the variant without matrix J -1 n in [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF], which yields nearly identical results), while the postprocessed method ( 15), [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] has order 3 for linear problems and is significantly more accurate by a factor of about 10. A line of slope 1 can be observed for the standard implicit Euler method. The θ = 1/2 method is exact for linear Brownian dynamics and the observed bias with magnitude 10 3 is due to the Monte-Carlo error (known [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF] to decay as O(T -1/2 )).

• In the stiff case (right picture of Figure 1), the implicit Euler method and the postprocessed methods ( 15),( 16) and ( 15),( 17) do not loose their efficiency. In contrast, the θ = 1/2 method becomes poorly accurate in the very stiff regime. This phenomena is well known already for severely stiff deterministic problems (the θ = 1/2 method is not L-stable), and this is highlighted in the stochastic context in [START_REF] Li | Effectiveness of implicit methods for stiff stochastic differential equations[END_REF]. A loss of accuracy is also observed for the variant of the method (15), [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF] where the term J -1 n providing L-stability is omitted in [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF]. Figure 3: Nonlinear nonstiff problem (48) and stiff problem(49). Implicit Euler method ((13) with θ = 1, order 1) and postprocessed versions ( 15),( 16) and ( 15),( 17) of order 2. The θ = 1/2-method (order 2) is also included in the nonstiff case (left picture). Error in d R (x 2 + x 2 1 ) 2 ρ(x)dx versus time stepsize h obtained using 10 trajectories on a long time interval of length T = 10 5 . The vertical bars indicate the standard deviation intervals due to the Monte-Carlo error.

Nonlinear test problems

We next consider two nonlinear problems. On the one hand, we consider a nonstiff problem in dimension d = 2 with potential

V (x) = (1 -x 2 1 ) 2 + x 4 2 -x + x 1 cos(x 2 ) + (x 2 + x 2 1 ) 2 , (48)
and on the other hand, we consider the stiff problem in dimension d = 3 with potential

V (x) = (1 -x 2 1 ) 2 + x 4 2 -x + x 3 cos(x 2 ) + 100(x 2 + x 2 1 ) 2 + 10 6 2 (x 1 -x 3 ) 2 . (49)
In Figure 2, we plot the associated Gibbs density measure [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF] together with five level curves for the components3 x 1 , x 2 for the nonstiff problem (48) (left pictures) and the stiff problem (49) (right pictures). It can be seen that the distribution is much narrower in the stiff case and concentrated close to the curve {x 2 + x 2 1 = 0}, due to the stiff term 100(x 2 + x 2 1 ) 2 in (49). In the bottom pictures of Figure 2, we also include one numerical trajectory of length T = 10 3 with stepsize h = 0.05 for the L-stable method [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF]. It fits well the invariant distribution for both the nonstiff and stiff cases.

For both problems, we take the initial conditions x i = i (i = 1, 2, 3) and consider ten long trajectories with length T = 10 4 . Analogously to the linear case, we compare in Figure 3 the results in the nonstiff and stiff cases. We plot for many constant stepsizes h the errors for the evaluation of the quantity R d (x 2 + x 2 1 ) 2 ρ(x)dx (a reference value for this integral is computed using an accurate quadrature method).

• In the nonstiff case (left picture of Figure 3), we observe the expected convergence lines of slope 1 and 2 for the implicit Euler method and the variant postprocessed method 13) with θ = 1, order 1) and postprocessed versions (new method ( 15),( 17) and variant ( 15),( 16)) of order 2. The θ = 1/2-method (order 2) is also included. Error for R d (|x| -1) 2 ρ(x)dx versus time stepsize h obtained using 10 trajectories on a long time interval of length T = 10 4 . ( 15),( 16), while the θ = 1/2 method and the postprocessed method ( 15),(17) (order 2 for nonlinear problems, order 3 for linear problems) reveal much more accurate.

• In the stiff case (right picture of Figure 3), the Newton iterations of the θ = 1/2method fail to converge for the considered stepsizes, a convergence curve is thus not included. The implicit Euler method still exhibits a convergence curve 1. The best accuracy is observed for the postprocessed method ( 15),( 17) of order two.

Finally, as a more challenging problem, we consider the following potential in dimension d = 10,

V (x) = 25(1 -|x|) 4 , |x| = d i=1 x 2 i . (50)
This is a modification of a standard nonlinear model for the motion of a spring in d dimensions considered as an illustrative expample in [START_REF] Sanz-Serna | Chain Monte Carlo and Numerical Differential Equation[END_REF]Sect. 7]. We plot in Figure 4 (left picture) the Gibbs density measure [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF] in dimension d = 2, and it can be seen that the distribution is concentrated close to the sphere {|x| = 1} centred at the origin with radius 1.

In Figure 4 (right picture), we plot the convergence curves for R d (|x| -1) 2 ρ(x)dx. Again, excellent performances of the postprocessed method ( 15),( 17) can be observed. For the considered stepsizes, it is about 20 times more accurate than the standard θ = 1/2 method, also of order two for the invariant measure. As a reference solution, we took the result of the θ = 1/2 method with stepsize h = 2 -4.5 ≃ 4.4 • 10 -2 .

Summary

In this paper, we extended to the stochastic context the methodology of effective order of Butcher. The approach of stochastic postprocessed integrators was illustrated in the case of Brownian dynamics, for which new schemes of order two were constructed based on the θ-method, in particular implicit schemes with favourable stability properties for the sampling of the invariant measure of ergodic systems. We rediscovered the second order non-Markovian scheme from [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] and analyzed in [START_REF] Leimkuhler | On the long-time integration of stochastic gradient systems[END_REF] which can be interpreted as a postprocessed integrator. This new insight permitted us to construct new high order integrators based on the θ-method, and to provide a rigorous error analysis based on the standard theory of Markovian schemes. Our approach is in principle not restricted to Brownian dynamics and could apply to larger classes of SDEs (Theorem 4.3). This will be investigated in future works.

Acknowledgements. This work was partially supported by the Fonds National Suisse, project No. 200020 144313/1.

Appendix

Proof of Theorem 4.6. For the scheme (44), a calculation analogous to [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally lipschitz coefficients[END_REF] shows that [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF] holds with Using integration by parts, the following identities where proved in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF],

A 1 φ = 1 2 φ ′′ (f, f ) + σ 2 2
φ ′′ (f, f ) = -φ ′ (f ′ f + div(f )f + 2 σ 2 f 2 f ) , σ 2
i φ ′′′ (f, e i , e i ) = φ ′ (σ 2 i f ′′ (e i , e i ) + 4f ′ f + 2 div(f )f + 4 σ 2 f 2 f ) , σ 2 ij φ (4) (e i , e i , e j , e j ) =i 2φ ′′′ (f, e i , e i ) , σ 2 i φ ′′ (f ′ e i , e i ) = -φ ′ (σ 2 i f ′′ (e i , e i ) + 2f ′ f ) ,

where we use the notation u = R d u(x)ρ(x)dx and the sums are for i, j = 1, . . . , d. For instance, to prove the first identity in (51), we integrate by parts and use ∂ j ρ = 2 σ 2 f j ρ, where ∂ j denotes the derivative ∂/(∂x j ),

φ ′′ (f, f ) = R d d i,j=1 ∂ i ∂ j φf i f j ρdy = - R d d i,j=1 ∂ i φ∂ j (f i f j ρ)dy = - d i,j=1 ∂ i φ(f i ∂ j f j + ∂ j f i f j + f i f j 2 σ 2 f j ) = -φ ′ (f ′ f + div(f )f + 2 σ 2 f 2 f ) .
The other identities in (51) can be proved analogously. We deduce

A 1 φ = a 2 + 2θa σ 2 2 φ ′ d i=1
f ′′ (e i , e i ) + -

2 Figure 2 :

 22 Figure2: Plots of the Gibbs density[START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF] for the nonlinear problems with potential (48) (nonstiff case, left pictures) and (49) (stiff case, right pictures). Four level curves are represented in solid lines, the three extrema are represented with crosses. In the bottom pictures, a numerical trajectory of the method (15) is included, with length T = 10 3 and stepsize h = 0.05.

Figure 4 :

 4 Figure 4: Nonlinear problem in dimension d = 10 with potential (50). Left picture: plot of the Gibbs density (10) in the case d = 2. Right picture: comparison of the implicit Euler method ((13) with θ = 1, order 1) and postprocessed versions (new method (15),(17) and variant (15),(16)) of order 2. The θ = 1/2-method (order 2) is also included. Error for

φ 8 d i,j=1 φ ( 4 )f

 84 ′′′ (e i , e i , f ) + σ 4 (e i , e i , e j , e j ) + θφ′ f ′ f + (1θ)a 2 + θ(a + 1) ′′ (e i , e i ) + a + θ σ 2 d i=1 φ ′′ (f ′ e i , e i ).

φσ 2 2 φ ′ d i=1 fφ

 i=1 ′′ (f ′ e i , e i )To improve the accuracy of the scheme (44), we now consider the postprocessor (45) where b and c ≥ 0 are fixed parameters. UsingA 1 φ = bφ ′ f + c 2 2 σ 2 ∆φ, we have [L, A 1 ]φ = (bc 2 ) 2 σ 2 φ ′ d i=1 f ′′ (e i , e i ) + (bc 2 )σ 2 d i=1 φ ′′ (f ′ e i , e i )Summing up then yields the value ofA 1 φ + [L, A 1 ]φ , a 2 + 2θa + bc 2 ′′ (e i , e i ) + -′′ (f ′ e i , e i )and this quantity is zero if (46) holds. We conclude the proof of Theorem 4.6 using Theorem 4.3.

Université de Genève, Section de mathématiques,

2-4 rue du Lièvre, CP 64, 1211 Genève 4. Gilles.Vilmart@unige.ch

Notice that the item 2. of Assumption

3.1 is automatically satisfied since g(ξ) T g(ξ) : ∇ 2 = σ 2 ∆ for the Brownian dynamics model[START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], where ∆ is the Laplacian operator.

Since (49) is in 3D, the right pictures of Figure2correspond to the Gibbs density measure for the marginal law of the components (x1, x2).


Julien Mairal

Piotr Koniusz

Zaid Harchaoui

Cordelia Schmid

Convolutional Kernel Networks

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We have recently seen a revival of attention given to convolutional neural networks (CNNs) [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] due to their high performance for large-scale visual recognition tasks [START_REF] Donahue | DeCAF: A deep convolutional activation feature for generic visual recognition[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Wan | Regularization of neural networks using dropconnect[END_REF]. The architecture of CNNs is relatively simple and consists of successive layers organized in a hierarchical fashion; each layer involves convolutions with learned filters followed by a pointwise non-linearity and a downsampling operation called "feature pooling". The resulting image representation has been empirically observed to be invariant to image perturbations and to encode complex visual patterns [START_REF] Zeiler | Visualizing and understanding convolutional networks[END_REF], which are useful properties for visual recognition. Training CNNs remains however difficult since high-capacity networks may involve billions of parameters to learn, which requires both high computational power, e.g., GPUs, and appropriate regularization techniques [START_REF] Goodfellow | Maxout networks[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Wan | Regularization of neural networks using dropconnect[END_REF].

The exact nature of invariance that CNNs exhibit is also not precisely understood. Only recently, the invariance of related architectures has been characterized; this is the case for the wavelet scattering transform [START_REF] Bruna | Invariant scattering convolution networks[END_REF] or the hierarchical models of [START_REF] Bouvrie | On invariance in hierarchical models[END_REF]. Our work revisits convolutional neural networks, but we adopt a significantly different approach than the traditional one. Indeed, we use kernels [START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF], which are natural tools to model invariance [START_REF] Decoste | Training invariant support vector machines[END_REF]. Inspired by the hierarchical kernel descriptors of [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF], we propose a reproducing kernel that produces multi-layer image representations.

Our main contribution is an approximation scheme called convolutional kernel network (CKN) to make the kernel approach computationally feasible. Our approach is a new type of unsupervised convolutional neural network that is trained to approximate the kernel map. Interestingly, our network uses non-linear functions that resemble rectified linear units [START_REF] Bengio | Learning deep architectures for AI[END_REF][START_REF] Wan | Regularization of neural networks using dropconnect[END_REF], even though they were not handcrafted and naturally emerge from an approximation scheme of the Gaussian kernel map.

By bridging a gap between kernel methods and neural networks, we believe that we are opening a fruitful research direction for the future. Our network is learned without supervision since the label information is only used subsequently in a support vector machine (SVM). Yet, we achieve competitive results on several datasets such as MNIST [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], CIFAR-10 [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF] and STL-10 [START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF] with simple architectures, few parameters to learn, and no data augmentation. Open-source code for learning our convolutional kernel networks is available on the first author's webpage.

Related Work

There have been several attempts to build kernel-based methods that mimic deep neural networks; we only review here the ones that are most related to our approach.

Arc-cosine kernels. Kernels for building deep large-margin classifiers have been introduced in [START_REF] Cho | Large-margin classification in infinite neural networks[END_REF]. The multilayer arc-cosine kernel is built by successive kernel compositions, and each layer relies on an integral representation. Similarly, our kernels rely on an integral representation, and enjoy a multilayer construction. However, in contrast to arc-cosine kernels: (i) we build our sequence of kernels by convolutions, using local information over spatial neighborhoods (as opposed to compositions, using global information); (ii) we propose a new training procedure for learning a compact representation of the kernel in a data-dependent manner.

Multilayer derived kernels. Kernels with invariance properties for visual recognition have been proposed in [START_REF] Bouvrie | On invariance in hierarchical models[END_REF]. Such kernels are built with a parameterized "neural response" function, which consists in computing the maximal response of a base kernel over a local neighborhood. Multiple layers are then built by iteratively renormalizing the response kernels and pooling using neural response functions. Learning is performed by plugging the obtained kernel in an SVM. In contrast to [START_REF] Bouvrie | On invariance in hierarchical models[END_REF], we propagate information up, from lower to upper layers, by using sequences of convolutions. Furthermore, we propose a simple and effective data-dependent way to learn a compact representation of our kernels and show that we obtain near state-of-the-art performance on several benchmarks. Hierarchical kernel descriptors. The kernels proposed in [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF][START_REF] Bo | Kernel descriptors for visual recognition[END_REF] produce multilayer image representations for visual recognition tasks. We discuss in details these kernels in the next section: our paper generalizes them and establishes a strong link with convolutional neural networks.

Convolutional Multilayer Kernels

The convolutional multilayer kernel is a generalization of the hierarchical kernel descriptors introduced in computer vision [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF][START_REF] Bo | Kernel descriptors for visual recognition[END_REF]. The kernel produces a sequence of image representations that are built on top of each other in a multilayer fashion. Each layer can be interpreted as a non-linear transformation of the previous one with additional spatial invariance. We call these layers image feature maps 1 , and formally define them as follows: Definition 1. An image feature map ϕ is a function ϕ : Ω → H, where Ω is a (usually discrete) subset of [0, 1] d representing normalized "coordinates" in the image and H is a Hilbert space.

For all practical examples in this paper, Ω is a two-dimensional grid and corresponds to different locations in a two-dimensional image. In other words, Ω is a set of pixel coordinates. Given z in Ω, the point ϕ(z) represents some characteristics of the image at location z, or in a neighborhood of z. For instance, a color image of size m × n with three channels, red, green, and blue, may be represented by an initial feature map ϕ 0 : Ω 0 → H 0 , where Ω 0 is an m × n regular grid, H 0 is the Euclidean space R 3 , and ϕ 0 provides the color pixel values. With the multilayer scheme, non-trivial feature maps will be obtained subsequently, which will encode more complex image characteristics. With this terminology in hand, we now introduce the convolutional kernel, first, for a single layer. Definition 2 (Convolutional Kernel with Single Layer). Let us consider two images represented by two image feature maps, respectively ϕ and ϕ ′ : Ω → H, where Ω is a set of pixel locations, and H is a Hilbert space. The one-layer convolutional kernel between ϕ and ϕ ′ is defined as

K(ϕ, ϕ ′ ) := z∈Ω z ′ ∈Ω ϕ(z) H ϕ ′ (z ′ ) H e -1 2β 2 z-z ′ 2 2 e -1 2σ 2 φ(z)-φ′ (z ′ ) 2 H , (1)
where β and σ are smoothing parameters of Gaussian kernels, and φ(z

) := (1/ ϕ(z) H ) ϕ(z) if ϕ(z) = 0 and φ(z) = 0 otherwise. Similarly, φ′ (z ′ ) is a normalized version of ϕ ′ (z ′ ). 2
It is easy to show that the kernel K is positive definite (see Appendix A). It consists of a sum of pairwise comparisons between the image features ϕ(z) and ϕ ′ (z ′ ) computed at all spatial locations z and z ′ in Ω. To be significant in the sum, a comparison needs the corresponding z and z ′ to be close in Ω, and the normalized features φ(z) and φ′ (z ′ ) to be close in the feature space H. The parameters β and σ respectively control these two definitions of "closeness". Indeed, when β is large, the kernel K is invariant to the positions z and z ′ but when β is small, only features placed at the same location z = z ′ are compared to each other. Therefore, the role of β is to control how much the kernel is locally shift-invariant. Next, we will show how to go beyond one single layer, but before that, we present concrete examples of simple input feature maps ϕ 0 : Ω 0 → H 0 .

Gradient map. Assume that H 0 = R2 and that ϕ 0 (z) provides the two-dimensional gradient of the image at pixel z, which is often computed with first-order differences along each dimension. Then, the quantity ϕ 0 (z) H0 is the gradient intensity, and φ0 (z) is its orientation, which can be characterized by a particular angle-that is, there exists θ in [0; 2π] such that φ0 (z) = [cos(θ), sin(θ)]. The resulting kernel K is exactly the kernel descriptor introduced in [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF][START_REF] Bo | Kernel descriptors for visual recognition[END_REF] for natural image patches.

Patch map. In that setting, ϕ 0 associates to a location z an image patch of size m × m centered at z. Then, the space H 0 is simply R m×m , and φ0 (z) is a contrast-normalized version of the patch, which is a useful transformation for visual recognition according to classical findings in computer vision [START_REF] Jarrett | What is the best multi-stage architecture for object recognition?[END_REF]. When the image is encoded with three color channels, patches are of size m × m × 3.

We now define the multilayer convolutional kernel, generalizing some ideas of [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF]. (i) choose a patch shape P k defined as a bounded symmetric subset of [-1, 1] d , and a set of coordinates Ω k such that for all location z k in Ω k , the patch

Definition 3 (Multilayer Convolutional Kernel). Let us consider a set Ω k-1 ⊆ [0, 1] d and
{z k } + P k is a subset of Ω k-1 ; 3 In other words, each coordinate z k in Ω k corresponds to a valid patch in Ω k-1 centered at z k .
(ii) define the convolutional kernel K k on the "patch" feature maps P k → H k-1 , by replacing in [START_REF] Bengio | Learning deep architectures for AI[END_REF]: Ω by P k , H by H k-1 , and σ, β by appropriate smoothing parameters σ k , β k . We denote by H k the Hilbert space for which the positive definite kernel K k is reproducing.

An image represented by a feature map

ϕ k-1 : Ω k-1 → H k-1 at layer k-1 is now encoded in the k-th layer as ϕ k : Ω k → H k , where for all z k in Ω k , ϕ k (z k ) is the representation in H k of the patch feature map z → ϕ k-1 (z k + z) for z in P k .
Concretely, the kernel K k between two patches of ϕ k-1 and ϕ ′ k-1 at respective locations z k and z ′ k is

z∈P k z ′ ∈P k ϕ k-1 (z k + z) ϕ ′ k-1 (z ′ k + z ′ ) e -1 2β 2 k z-z ′ 2 2 e -1 2σ 2 k φk-1 (z k +z)-φ′ k-1 (z ′ k +z ′ ) 2 , (2)
where . is the Hilbertian norm of H k-1 . In Figure 1(a), we illustrate the interactions between the sets of coordinates Ω k , patches P k , and feature spaces H k across layers. For two-dimensional grids, a typical patch shape is a square, for example

P := {-1/n, 0, 1/n} × {-1/n, 0, 1/n} for a 3 × 3 patch in an image of size n × n.
Information encoded in the k-th layer differs from the (k-1)-th one in two aspects: first, each point ϕ k (z k ) in layer k contains information about several points from the (k-1)-th layer and can possibly represent larger patterns; second, the new feature map is more locally shift-invariant than the previous one due to the term involving the parameter β k in (2).

The multilayer convolutional kernel slightly differs from the hierarchical kernel descriptors of [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF] but exploits similar ideas. Bo et al. [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF] define indeed several ad hoc kernels for representing local information in images, such as gradient, color, or shape. These kernels are close to the one defined in (1) but with a few variations. Some of them do not use normalized features φ(z), and these kernels use different weighting strategies for the summands of (1) that are specialized to the image modality, e.g., color, or gradient, whereas we use the same weight ϕ(z) H ϕ ′ (z ′ ) H for all kernels. The generic formulation (1) that we propose may be useful per se, but our main contribution comes in the next section, where we use the kernel as a new tool for learning convolutional neural networks.

Ω 0 ϕ 0 (z 0 ) ∈ H 0 {z 1 } + P 1 ϕ 1 (z 1 ) ∈ H 1 Ω 1 {z 2 } + P 2 Ω 2 ϕ 2 (z 2 ) ∈ H 2
(a) Hierarchy of image feature maps.

Ω ′ k-1 ξ k-1 (z) ψ k-1 (z k-1 ) (patch extraction) {z k-1 }+P ′ k-1 convolution + non-linearity p k ζ k (z k-1 ) Ω k-1 Gaussian filtering + downsampling = pooling Ω ′ k ξ k (z)

Training Invariant Convolutional Kernel Networks

Generic schemes have been proposed for approximating a non-linear kernel with a linear one, such as the Nyström method and its variants [START_REF] Bo | Efficient match kernel between sets of features for visual recognition[END_REF][START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF], or random sampling techniques in the Fourier domain for shift-invariant kernels [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. In the context of convolutional multilayer kernels, such an approximation is critical because computing the full kernel matrix on a database of images is computationally infeasible, even for a moderate number of images (≈ 10 000) and moderate number of layers. For this reason, Bo et al. [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF] use the Nyström method for their hierarchical kernel descriptors.

In this section, we show that when the coordinate sets Ω k are two-dimensional regular grids, a natural approximation for the multilayer convolutional kernel consists of a sequence of spatial convolutions with learned filters, pointwise non-linearities, and pooling operations, as illustrated in Figure 1(b). More precisely, our scheme approximates the kernel map of K defined in (1) at layer k by finite-dimensional spatial maps ξ k : Ω ′ k → R p k , where Ω ′ k is a set of coordinates related to Ω k , and p k is a positive integer controlling the quality of the approximation. Consider indeed two images represented at layer k by image feature maps ϕ k and ϕ ′ k , respectively. Then,

(A) the corresponding maps ξ k and ξ ′ k are learned such that K(ϕ k-1 , ϕ ′ k-1 ) ≈ ξ k , ξ ′ k ,
(z k ) in H k admit finite-dimensional approximations ψ k (z k ) in R |P ′ k |p k ; as illustrated in Figure 1(b), ψ k (z k ) is a patch from ξ k centered at location z k with shape P ′ k ; (C) an activation map ζ k : Ω k-1 → R p k is computed from ξ k-1
by convolution with p k filters followed by a non-linearity. The subsequent map ξ k is obtained from ζ k by a pooling operation.

We call this approximation scheme a convolutional kernel network (CKN). In comparison to CNNs, our approach enjoys similar benefits such as efficient prediction at test time, and involves the same set of hyper-parameters: number of layers, numbers of filters p k at layer k, shape P ′ k of the filters, sizes of the feature maps. The other parameters β k , σ k can be automatically chosen, as discussed later. Training a CKN can be argued to be as simple as training a CNN in an unsupervised manner [START_REF] Ranzato | Unsupervised learning of invariant feature hierarchies with applications to object recognition[END_REF] since we will show that the main difference is in the cost function that is optimized.

Fast Approximation of the Gaussian Kernel

A key component of our formulation is the Gaussian kernel. We start by approximating it by a linear operation with learned filters followed by a pointwise non-linearity. Our starting point is the next lemma, which can be obtained after a simple calculation.

Lemma 1 (Linear expansion of the Gaussian Kernel). For all x and x ′ in R m , and σ > 0,

e -1 2σ 2 x-x ′ 2 2 = 2 πσ 2 m 2 w∈R m e -1 σ 2 x-w 2 2 e -1 σ 2 x ′ -w 2 2 dw. (3)
The lemma gives us a mapping of any x in R m to the function

w → √ Ce -(1/σ 2 ) x-w 2 2 in L 2 (R m )
, where the kernel is linear, and C is the constant in front of the integral. To obtain a finite-dimensional representation, we need to approximate the integral with a weighted finite sum, which is a classical problem arising in statistics (see [START_REF] Wahba | Spline models for observational data[END_REF] and chapter 8 of [START_REF] Bottou | Large-Scale Kernel Machines (Neural Information Processing[END_REF]). Then, we consider two different cases.

Small dimension, m ≤ 2. When the data lives in a compact set of R m , the integral in (3) can be approximated by uniform sampling over a large enough set. We choose such a strategy for two types of kernels from Eq. ( 1): (i) the spatial kernels e

-1 2β 2 z-z ′ 2 2 ; (ii) the terms e -( 1 2σ 2 ) φ(z)-φ′ (z ′ ) 2 H
when ϕ is the "gradient map" presented in Section 2. In the latter case, H = R 2 and φ(z) is the gradient orientation. We typically sample a few orientations as explained in Section 4.

Higher dimensions. To prevent the curse of dimensionality, we learn to approximate the kernel on training data, which is intrinsically low-dimensional. We optimize importance weights η

= [η l ] p l=1 in R p + and sampling points W = [w l ] p l=1 in R m×p on n training pairs (x i , y i ) i=1,...,n in R m × R m : min η∈R p + ,W∈R m×p 1 n n i=1 e -1 2σ 2 xi-yi 2 2 - p l=1 η l e -1 σ 2 xi-w l 2 2 e -1 σ 2 yi-w l 2 2 2 . (4)
Interestingly, we may already draw some links with neural networks. When applied to unit-norm vectors x i and y i , problem (4) produces sampling points w l whose norm is close to one. After learning, a new unit-norm point x in R m is mapped to the vector [ √ η l e -(1/σ 2 ) x-w l 2 2 ] p l=1 in R p , which may be written as [f (w ⊤ l x)] p l=1 , assuming that the norm of w l is always one, where f is the function u → e (2/σ 2 )(u-1) for u = w ⊤ l x in [-1, 1]. Therefore, the finite-dimensional representation of x only involves a linear operation followed by a non-linearity, as in typical neural networks. In Figure 2, we show that the shape of f resembles the "rectified linear unit" function [START_REF] Wan | Regularization of neural networks using dropconnect[END_REF].

u f (u) f (u) = e (2/σ 2 f (u) = max(u, 0) 0 1 -1
Figure 2: In dotted red, we plot the "rectified linear unit" function u → max(u, 0). In blue, we plot non-linear functions of our network for typical values of σ that we use in our experiments.

Approximating the Multilayer Convolutional Kernel

We have now all the tools in hand to build our convolutional kernel network. We start by making assumptions on the input data, and then present the learning scheme and its approximation principles.

The zeroth layer. We assume that the input data is a finite-dimensional map ξ 0 : Ω ′ 0 → R p0 , and that ϕ 0 : Ω 0 → H 0 "extracts" patches from ξ 0 . Formally, there exists a patch shape P ′ 0 such that

Ω ′ 0 = Ω 0 + P ′ 0 , H 0 = R p0|P ′ 0 |
, and for all z 0 in Ω 0 , ϕ 0 (z 0 ) is a patch of ξ 0 centered at z 0 . Then, property (B) described at the beginning of Section 3 is satisfied for k = 0 by choosing ψ 0 = ϕ 0 . The examples of input feature maps given earlier satisfy this finite-dimensional assumption: for the gradient map, ξ 0 is the gradient of the image along each direction, with p 0 = 2, P ′ 0 = {0} is a 1×1 patch, Ω 0 = Ω ′ 0 , and ϕ 0 = ξ 0 ; for the patch map, ξ 0 is the input image, say with p 0 = 3 for RGB data. The convolutional kernel network. The zeroth layer being characterized, we present in Algorithms 1 and 2 the subsequent layers and how to learn their parameters in a feedforward manner. It is interesting to note that the input parameters of the algorithm are exactly the same as a CNN-that is, number of layers and filters, sizes of the patches and feature maps (obtained here via the subsampling factor). Ultimately, CNNs and CKNs only differ in the cost function that is optimized for learning the filters and in the choice of non-linearities. As we show next, there exists a link between the parameters of a CKN and those of a convolutional multilayer kernel.

Algorithm 1 Convolutional kernel network -learning the parameters of the

k-th layer. input ξ 1 k-1 , ξ 2 k-1 , . . . : Ω ′ k-1 → R p k-1
(sequence of (k-1)-th maps obtained from training images); P ′ k-1 (patch shape); p k (number of filters); n (number of training pairs); 1: extract at random n pairs (x i , y i ) of patches with shape P ′ k-1 from the maps ξ 1 k-1 , ξ 2 k-1 , . . .; 2: if not provided by the user, set σ k to the 0.1 quantile of the data ( x iy i 2 ) n i=1 ; 3: unsupervised learning: optimize (4) to obtain the filters

W k in R |P ′ k-1 |p k-1 ×p k and η k in R p k ; output W k , η k ,
k : Ω k-1 → R p k as ζ k : z → ψ k-1 (z) 2 √ η kl e -1 σ 2 k ψk-1 (z)-w kl 2 2 p k l=1 , (5)
where ψ k-1 (z) is a vector representing a patch from ξ k-1 centered at z with shape P ′ k-1 , and the vector ψk-1 (z) is an ℓ 2 -normalized version of ψ k-1 (z). This operation can be interpreted as a spatial convolution of the map ξ k-1 with the filters w kl followed by pointwise non-linearities; 2: set β k to be γ k times the spacing between two pixels in Ω k-1 ; 3: feature pooling: Ω ′ k is obtained by subsampling Ω k-1 by a factor γ k and we define a new map ξ k : Ω ′ k → R p k obtained from ζ k by linear pooling with Gaussian weights:

ξ k : z → 2/π u∈Ω k-1 e -1 β 2 k u-z 2 2 ζ k (u). ( 6
)
output ξ k : Ω ′ k → R p k (new map);
Approximation principles. We proceed recursively to show that the kernel approximation property (A) is satisfied; we assume that (B) holds at layer k-1, and then, we show that (A) and (B) also hold at layer k. This is sufficient for our purpose since we have previously assumed (B) for the zeroth layer. Given two images feature maps ϕ k-1 and ϕ ′ k-1 , we start by approximating K(ϕ k-1 , ϕ ′ k-1 ) by replacing ϕ k-1 (z) and ϕ ′ k-1 (z ′ ) by their finite-dimensional approximations provided by (B):

K(ϕ k-1 , ϕ ′ k-1 ) ≈ z,z ′ ∈Ω k-1 ψ k-1 (z) 2 ψ ′ k-1 (z ′ ) 2 e -1 2β 2 k z-z ′ 2 2 e -1 2σ 2 k ψk-1 (z)-ψ′ k-1 (z ′ ) 2 2 . (7)
Then, we use the finite-dimensional approximation of the Gaussian kernel involving σ k and

K(ϕ k-1 , ϕ ′ k-1 ) ≈ z,z ′ ∈Ω k-1 ζ k (z) ⊤ ζ ′ k (z ′ )e -1 2β 2 k z-z ′ 2 2 , (8)
where ζ k is defined in ( 5) and ζ ′ k is defined similarly by replacing ψ by ψ′ . Finally, we approximate the remaining Gaussian kernel by uniform sampling on Ω ′ k , following Section 3.1. After exchanging sums and grouping appropriate terms together, we obtain the new approximation

K(ϕ k-1 , ϕ ′ k-1 ) ≈ 2 π u∈Ω ′ k z∈Ω k-1 e -1 β 2 k z-u 2 2 ζ k (z) ⊤ z ′ ∈Ω k-1 e -1 β 2 k z ′ -u 2 2 ζ ′ k (z ′ ) , (9)
where the constant 2/π comes from the multiplication of the constant 2/(πβ 2 k ) from (3) and the weight β 2 k of uniform sampling orresponding to the square of the distance between two pixels of Ω ′ k . 4As a result, the right-hand side is exactly ξ k , ξ ′ k , where ξ k is defined in [START_REF] Bottou | Large-Scale Kernel Machines (Neural Information Processing[END_REF], giving us property (A). It remains to show that property (B) also holds, specifically that the quantity (2) can be approximated by the Euclidean inner-product ψ k (z k ), ψ ′ k (z ′ k ) with the patches ψ k (z k ) and ψ ′ k (z ′ k ) of shape P ′ k ; we assume for that purpose that P ′ k is a subsampled version of the patch shape P k by a factor γ k .

We remark that the kernel (2) is the same as (1) applied to layer k-1 by replacing Ω k-1 by {z k }+P k . By doing the same substitution in (9), we immediately obtain an approximation of (2). Then, all Gaussian terms are negligible for all u and z that are far from each other-say when u-z 2 ≥ 2β k . Thus, we may replace the sums u∈Ω

′ k z,z ′ ∈{z k }+P k by u∈{z k }+P ′ k z,z ′ ∈Ω k-1
, which has the same set of "non-negligible" terms. This yields exactly the approximation ψ k (z k ), ψ ′ k (z ′ k ) . Optimization. Regarding problem (4), stochastic gradient descent (SGD) may be used since a potentially infinite amount of training data is available. However, we have preferred to use L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] on 300 000 pairs of randomly selected training data points, and initialize W with the K-means algorithm. L-BFGS-B is a parameter-free state-of-the-art batch method, which is not as fast as SGD but much easier to use. We always run the L-BFGS-B algorithm for 4 000 iterations, which seems to ensure convergence to a stationary point. Our goal is to demonstrate the preliminary performance of a new type of convolutional network, and we leave as future work any speed improvement.

Experiments

We now present experiments that were performed using Matlab and an L-BFGS-B solver [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] interfaced by Stephen Becker. Each image is represented by the last map ξ k of the CKN, which is used in a linear SVM implemented in the software package LibLinear [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF]. These representations are centered, rescaled to have unit ℓ 2 -norm on average, and the regularization parameter of the SVM is always selected on a validation set or by 5-fold cross-validation in the range 2 i , i = -15 . . . , 15.

The patches P ′ k are typically small; we tried the sizes m × m with m = 3, 4, 5 for the first layer, and m = 2, 3 for the upper ones. The number of filters p k in our experiments is in the set {50, 100, 200, 400, 800}. The downsampling factor γ k is always chosen to be 2 between two consecutive layers, whereas the last layer is downsampled to produce final maps ξ k of a small size-say, 5×5 or 4×4. For the gradient map ϕ 0 , we approximate the Gaussian kernel e (1/σ 2 1 ) ϕ0(z)-ϕ ′ 0 (z ′ ) H 0 by uniformly sampling p 1 = 12 orientations, setting σ 1 = 2π/p 1 . Finally, we also use a small offset ε to prevent numerical instabilities in the normalization steps ψ(z) = ψ(z)/ max( ψ(z) 2 , ε).

Discovering the Structure of Natural Image Patches

Unsupervised learning was first used for discovering the underlying structure of natural image patches by Olshausen and Field [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. Without making any a priori assumption about the data except a parsimony principle, the method is able to produce small prototypes that resemble Gabor wavelets-that is, spatially localized oriented basis functions. The results were found impressive by the scientific community and their work received substantial attention. It is also known that such results can also be achieved with CNNs [START_REF] Ranzato | Unsupervised learning of invariant feature hierarchies with applications to object recognition[END_REF]. We show in this section that this is also the case for convolutional kernel networks, even though they are not explicitly trained to reconstruct data.

Following [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF], we randomly select a database of 300 000 whitened natural image patches of size 12 × 12 and learn p = 256 filters W using the formulation (4). We initialize W with Gaussian random noise without performing the K-means step, in order to ensure that the output we obtain is not an artifact of the initialization. In Figure 3, we display the filters associated to the top-128 largest weights η l . Among the 256 filters, 197 exhibit interpretable Gabor-like structures and the rest was less interpretable. To the best of our knowledge, this is the first time that the explicit kernel map of the Gaussian kernel for whitened natural image patches is shown to be related to Gabor wavelets.

Digit Classification on MNIST

The MNIST dataset [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] consists of 60 000 images of handwritten digits for training and 10 000 for testing. We use two types of initial maps in our networks: the "patch map", denoted by CNK-PM and the "gradient map", denoted by CNK-GM. We follow the evaluation methodology of [START_REF] Ranzato | Unsupervised learning of invariant feature hierarchies with applications to object recognition[END_REF] Figure 3: Filters obtained by the first layer of the convolutional kernel network on natural images. gives better results than one layer. More details about the network architectures are provided in the supplementary material. In general, our method achieves a state-of-the-art accuracy for this task since lower error rates have only been reported by using data augmentation [START_REF] Ciresan | Multi-column deep neural networks for image classification[END_REF].

Visual Recognition on CIFAR-10 and STL-10

We now move to the more challenging datasets CIFAR-10 [20] and STL-10 [START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF]. We select the best architectures on a validation set of 10 000 examples from the training set for CIFAR-10, and by 5-fold cross-validation on STL-10. We report in Table 2 results for CKN-GM, defined in the previous section, without exploiting color information, and CKN-PM when working on raw RGB patches whose mean color is subtracted. The best selected models have always two layers, with 800 filters for the top layer. Since CKN-PM and CKN-GM exploit a different information, we also report a combination of such two models, CKN-CO, by concatenating normalized image representations together. The standard deviations for STL-10 was always below 0.7%. Our approach appears to be competitive with the state of the art, especially on STL-10 where only one method does better than ours, despite the fact that our models only use 2 layers and require learning few parameters. Note that better results than those reported in Table 2 have been obtained in the literature by using either data augmentation (around 90% on CIFAR-10 for [START_REF] Goodfellow | Maxout networks[END_REF][START_REF] Wan | Regularization of neural networks using dropconnect[END_REF]), or external data (around 70% on STL-10 for [START_REF] Swersky | Multi-task Bayesian optimization[END_REF]). We are planning to investigate similar data manipulations in the future.

Conclusion

In this paper, we have proposed a new methodology for combining kernels and convolutional neural networks. We show that mixing the ideas of these two concepts is fruitful, since we achieve near state-of-the-art performance on several datasets such as MNIST, CIFAR-10, and STL10, with simple architectures and no data augmentation. Some challenges regarding our work are left open for the future. The first one is the use of supervision to better approximate the kernel for the prediction task.

The second consists in leveraging the kernel interpretation of our convolutional neural networks to better understand the theoretical properties of the feature spaces that these networks produce.

A Positive Definiteness of K

To show that the kernel K defined in (1) is positive definite (p.d.), we simply use elementary rules from the kernel literature described in Sections 2.3.2 and 3.4.1 of [START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF]. A linear combination of p.d. kernels with nonnegative weights is also p.d. (see Proposition 3.22 of [START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF]), and thus it is sufficient to show that for all z, z ′ in Ω, the following kernel on Ω → H is p.d.:

(ϕ, ϕ ′ ) → ϕ(z) H ϕ ′ (z ′ ) H e -1 2σ 2 φ(z)-φ′ (z ′ ) 2 H .
Specifically, it is also sufficient to show that the following kernel on H is p.d.:

(φ, φ ′ ) → φ H φ ′ H e -1 2σ 2 φ φ H - φ ′ φ ′ H 2 H .
with the convention φ/ φ H = 0 if φ = 0. This is a pointwise product of two kernels and is p.d. when each of the two kernels is p.d. The first one is obviously p.d.: (φ, φ ′ ) → φ H φ ′ H . The second one is a composition of the Gaussian kernel-which is p.d.-, with feature maps φ/ φ H of a normalized linear kernel in H. This composition is p.d. according to Proposition 3.22, item (v) of [START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF] since the normalization does not remove the positive-definiteness property.

B List of Architectures Reported in the Experiments

We present in details the architectures used in the paper in Table 3. Table 3: List of architectures reported in the paper. N is the number of layers; p 1 and p 2 represent the number of filters are each layer; m 1 and m 2 represent the size of the patches P ′ 1 and P ′ 2 that are of size m 1 × m 1 and m 2 × m 2 on their respective feature maps ζ 1 and ζ 2 ; γ 1 is the subsampling factor between layer 1 and layer 2; S is the size of the output feature map, and the last column indicates the number of parameters that the network has to learn.

  (b) Zoom between layer k-1 and k of the CKN.

Figure 1 :

 1 Figure 1: Left: concrete representation of the successive layers for the multilayer convolutional kernel. Right: one layer of the convolutional neural network that approximates the kernel.

Algorithm 2

 2 and σ k (smoothing parameter); Convolutional kernel network -computing the k-th map form the (k-1)-th one.input ξ k-1 : Ω ′ k-1 → R p k-1 (input map); P ′ k-1 (patch shape); γ k ≥ 1 (subsampling factor); p k (number of filters); σ k (smoothing parameter); W k = [w kl ] p k l=1 and η k = [η kl ] p k l=1 (layer parameters); 1: convolution and non-linearity: define the activation map ζ

  where ., . is the Euclidean inner-product acting as if ξ k and ξ ′

k were vectors in R |Ω ′ k |p k ; (B) the set Ω ′ k is linked to Ω k by the relation Ω ′ k = Ω k + P ′ k where P ′

k is a patch shape, and the quantities ϕ k

Table 1 :

 1 Test error in % for various approaches on the MNIST dataset without data augmentation. The numbers in parentheses represent the size p 1 and p 2 of the feature maps at each layer.for comparison when varying the training set size. We select the regularization parameter of the SVM by 5-fold cross validation when the training size is smaller than 20 000, or otherwise, we keep 10 0000 examples from the training set for validation. We report in Table1the results obtained for four simple architectures. CKN-GM1 is the simplest one: its second layer uses 3 × 3 patches and only p 2 = 50 filters, resulting in a network with 5 400 parameters. Yet, it achieves an outstanding performance of 0.58% error on the full dataset. The best performing, CKN-GM2, is similar to CKN-GM1 but uses p 2 = 400 filters. When working with raw patches, two layers (CKN-PM2)

	Tr. CNN Scat-1 Scat-2 CKN-GM1 CKN-GM2 CKN-PM1 CKN-PM2 [32] [18] [19] size [25] [8] [8] (12/50) (12/400) (200) (50/200)
	300 7.18	4.7	5.6	4.39	4.24	5.98	4.15	NA
	1K 3.21	2.3	2.6	2.60	2.05	3.23	2.76	NA
	2K 2.53	1.3	1.8	1.85	1.51	1.97	2.28	NA
	5K 1.52	1.03	1.4	1.41	1.21	1.41	1.56	NA
	10K 0.85	0.88	1	1.17	0.88	1.18	1.10	NA
	20K 0.76	0.79	0.58	0.89	0.60	0.83	0.77	NA
	40K 0.65	0.74	0.53	0.68	0.51	0.64	0.58	NA
	60K 0.53	0.70	0.4	0.58	0.39	0.63	0.53	0.47 0.45 0.53

Table 2 :

 2 Classification accuracy in % on CIFAR-10 and STL-10 without data augmentation.

In the kernel literature, "feature map" denotes the mapping between data points and their representation in a reproducing kernel Hilbert space (RKHS)[START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF]. Here, feature maps refer to spatial maps representing local image characteristics at everly location, as usual in the neural network literature[START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF].

When Ω is not discrete, the notation in (1) should be replaced by the Lebesgue integral in the paper.

For two sets A and B, the Minkowski sum A + B is defined as {a + b : a ∈ A, b ∈ B}.

The choice of β k in Algorithm 2 is driven by signal processing principles. The feature pooling step can indeed be interpreted as a downsampling operation that reduces the resolution of the map from Ω k-1 to Ω k by using a Gaussian anti-aliasing filter, whose role is to reduce frequencies above the Nyquist limit.

Acknowledgments

This work was partially supported by grants from ANR (project MACARON ANR-14-CE23-0003-01), MSR-Inria joint centre, European Research Council (project ALLEGRO), CNRS-Mastodons program (project GARGANTUA), and the LabEx PERSYVAL-Lab (ANR-11-LABX-0025).