Houssam Alrachid
email: houssam.alrachid@u-pec.fr

Letif Mones
email: l.mones@warwick.ac.uk

Christoph Ortner
email: c.ortner@warwick.ac.uk

P-Langevin_Vh1 Some Remarks on Preconditioning Molecular Dynamics

Keywords: Preconditioned Overdamped Langevin algorithm, Asymptotic Variance, Central Limit Theorem, Model Hamiltonians, Lattice Model

We consider a Preconditioned Overdamped Langevin algorithm that does not alter the invariant distribution (up to controllable discretisation errors) and ask whether preconditioning improves the standard model in terms of reducing the asymptotic variance and of accelerating convergence to equilibrium. We present a detailed study of the dependence of the asymptotic variance on preconditioning in some elementary toy models related to molecular simulation. Our theoretical results are supported by numerical simulations.

Introduction

The problem of convergence to equilibrium for diffusion processes has attracted considerable attention in recent years. Due to the significant computational cost associated with MCMC type algorithms it is important to understand, and where possible accelerate, the convergence to equilibrium of systems in statistical physics, materials science, biochemistry, machine learning and many other areas.

In such applications it is typically necessary to compute expectations of the form

µ(f ) = R N f (x)µ(dx) (1.1)
of an observable f with respect to a target probability distribution µ(dx) on R N , where µ is of the form µ(dx) = Z -1 µ e -E N dx, November 7, 2017 11:9 P-Langevin_Vh1

2 HA,CO

where E N : R N → R is e.g. a potential energy and Z -1 µ the normalisation constant. Note that we have absorbed temperature into E N . The usual difficulty is that the integration (1.1) cannot be performed directly due to the high dimensionality of the problem and MCMC methods are instead employed.

Ill-conditioning of E N , which can be induced by a variety of mechanisms, but in particular high-dimensionality (large N ), is a common challenge to overcome in order to construct an efficient sampling scheme. An attractive approach is to precondition the MCMC algorithm. The algorithm is transformed by a well-chosen operator (the preconditioner) in a way that does not alter the invariant measure but (hopefully) accelerates convergence.

The purpose of this paper is to explore to what extent (or, whether at all) preconditioning of a Langevin-type algorithm helps to accelerate the computation of expectations. Our study is motivated by recent advances such as the Riemannian Manifold MALA [START_REF] Girolami | Riemann Manifold Langevin and Hamiltonian Monte Carlo methods[END_REF], Stochastic Newton Methods [START_REF] Martin | A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion[END_REF], non-reversible diffusions [START_REF] Duncan | Variance Reduction using Nonreversible Langevin Samplers[END_REF][START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF], optimal scaling for Langevin algorithms [START_REF] Beskos | Optimal scalings for local Metropolis-Hastings chains on nonproduct targets in high dimensions[END_REF][START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF] and affine-invariant sampler [START_REF] Hou | An affineinvariant sampler for exoplanet fitting and discovery in radial velocity data[END_REF].

All of these references present different variants of preconditioning or related modifications of MCMC methods and result in an improved performance. Reviewing all these different approaches would go beyond the scope of this introduction, however, we mention the Riemannian Manifold Monte Carlo method [START_REF] Girolami | Riemann Manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Xifara | Langevin diffusions and the Metropolis-adjusted Langevin algorithm[END_REF] as the main motivation for the present study. In this method, preconditioning is understood as performing MCMC on a Riemannian manifold defined via a "metric" P (x), the preconditioner. A range of preconditioners, ideally the hessian ∇ 2 E N (x) or a positive definite matrix P (x) closely approximating the Hessian, are tested. In a broad range of examples (including, e.g., logistic regression, a stochastic volatility model, and an ODE inference example) it has been shown in numerical tests that a well-chosen preconditioner improves both the mixing time as well as the convergence of the probability density functions to the target measure, with speed-ups ranging from moderate O(1) factors to orders of magnitude depending on the application.

Motivated by these promising results we applied analogous preconditioned sampling algorithms to some model molecular systems, but did not always (rarely in fact) observe the speed-ups we expected. For instance, in Figure 1 we show the convergence of reconstructed free energy profiles from metadynamics [START_REF] Laio | Escaping free-energy minima[END_REF] simulations comparing the unpreconditioned and preconditioned dynamics using a Hessianbased preconditioner with varying parameter sets. This type of preconditioner usually yields at least an order of magnitude speed-up in typical geometry optimisation tasks, but fails to accelerate the assembly of the free energy surface (see Appendix A.2 for the precise setup of the test). The question naturally arises whether this failure is due to a lack of fine-tuning or due to a more fundamental limitation.

Thus, to understand better these observations we will study some elementary analytical and numerical examples, which capture some essential characteristics of typical molecular systems, but where explicit results can still be obtained. The origin of the difficulty comes from the fact that it is highly dependent not only on µ but also on the observable f whether preconditioning can achieve a significant (or, any) speed-up. We will demonstrate that for some typical observables f , even moderate preconditioning can achieve significant speed-ups while for other, equally common, observables no speed-up should be at all expected even if E N is highly ill-conditioned. Moreover, we will show that the dimensionality (1D, 2D, 3D) of the molecular structure plays a crucial role. While our discussion is primarily motivated by applications in molecular simulations, it should straightforwardly adapt the arguments and findings to other application areas.

Langevin Algorithm

The most commonly employed algorithms in molecular simulation are based on discretising the Langevin equation, but for the sake of simplicity we will focus on the overdamped Langevin equation,

dX t = -∇E N (X t )dt + √ 2dW t , (1.2)
where W t is a N -dimensional standard Brownian motion. Under mild technical conditions on E N and X t it is known that the dynamics (X t ) t≥0 is ergodic with respect to the measure µ (see [START_REF] Roberts | Exponential convergence of langevin distributions and their discrete approximations[END_REF]). Discretising in time,

X m+1 = X m -δ∇E N (X m ) + √ 2δR m , (1.3)
where R m ∼ N (0, I N ×N ) and δ is a parameter quantifying the size of the discrete time increment (time-step), we obtain a Markov chain with invariant measure μ, where typically μ -µ ∼ O(t -1/2 ) + O(δ). Here, O(t -1/2 ) represents the statistical error due to the finite length of the simulation, while O(δ) represents the bias due to time discretisation. Adding metropolisation to (1.3) leads to the Metropolis Adjusted Langevin Algorithm (MALA). It is well known that the measure µ is invariant for the MALA [START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Roberts | Exponential convergence of langevin distributions and their discrete approximations[END_REF]. In molecular simulation it is common to assume that the bias in μ is negligible compared to modelling and statistical errors and therefore we will not consider metropolisation in the present paper. The resulting algorithm is called the Unadjusted Langevin Algorithm (ULA).

As a matter of fact, we will argue in Section 2.4 that the issues we are addressing are largely unrelated to the time-discretisation, hence we will focus on the continuous process (1.2) instead of the Markov chain (1.3).

Thus, given an observable f ∈ L 1 (µ), we are interested in quantifying the convergence

T (f ) := 1 T T 0 f (X t )dt → µ(f )
for µ-a.e X 0 .

(1.4)

Under additional assumptions on µ and f , this convergence result is accompanied by a central limit theorem which characterises the asymptotic distribution of the fluctuations, i.e.

√ t ( t (f ) -µ(f )) D -→ N (0, σ 2 f ), (1.5)
where σ 2 f is known as the asymptotic variance for the observable f (see [START_REF] Cattiaux | Central limit theorems for additive functionals of ergodic Markov diffusions processes[END_REF][START_REF] Kipnis | Central limit theorem for additive functionals of reversible Markov processes and applications to simple exclusions[END_REF]). In Section 4 we will explicitly compute σ 2 f for some simple energy functionals and observables, which mimic typical objects of interest in molecular simulations and demonstrate that for some typical observables, σ f may be strongly dependent on the conditioning of E N , in particular on system size N , while for others the dependence on N is negligible.

Preconditioned Langevin Algorithm

For a (fixed) preconditioner P ∈ R N ×N , symmetric positive definite, we consider the preconditioned Overdamped Langevin dynamics (P -Langevin),

dX P t = -P -1 ∇E N (X P t )dt + √ 2P -1/2 dW t . (1.6)
The standard overdamped Langevin dynamics is recovered by taking P = I. It is in principle possible to allow P = P (X) but for simplicity we will not consider this in November 7, 2017 11:9

P-Langevin_Vh1 pol 5
the present work. Note that, if P is fixed then the coordinate transform Z = P 1/2 X allows us to easily lift results from (1.2) to (1.6); see Section 2.1 for more details. The preconditioner does not affect the invariance property of the diffusion process, i.e., the target measure µ is still invariant for the P-Langevin process (1.6). However, it can affect the convergence of the process to the invariant measure. That is, for X P t ∼ µ P t = ϕ P t dx, we will characterise in Section 3 the rate of convergence of ϕ P t to ϕ ∞ (ϕ ∞ denotes the density of µ) and demonstrate how preconditioning improves this rate.

Moreover, we also obtain

P T (f ) := 1 T T 0 f (X P t )dt → µ(f ), for µ-a.e X 0 . (1.7)
and analogously to the standard Langevin dynamics, a central limit theorem characterizes the asymptotic distribution of the fluctuations,

√ t P t (f ) -µ(f ) D -→ N (0, σ 2 f,P ), (1.8)
where σ 2 f,P is the asymptotic variance of f under P -Langevin dynamics. The main aim of our paper is to present several simplified but still realistic examples at which we can observe whether or not preconditioning accelerates sampling in the sense that it achieves a reduction in the asymptotic variance, i.e., σ 2

f,P σ 2 f .

Outline

The rest of the paper is organised as follows. In Section 2 we present the model Hamiltonians, E N , to motivate some key assumptions that we make throughout our analysis. We also recall the coordinate transform that we use to reduce P -Langevin to standard Langevin, and we explain why the time-discretisation is a negligible component and can therefore be ignored for our analysis.

In Section 3 we show a long time convergence result to the invariant measure in the P -Langevin process. For some quadratic model Hamiltonians we can then precisely quantify the speed-up afforded through preconditioning.

Section 4 is devoted to describe how the central limit theorem (1.8) arises from the solution of the Poisson equation associated with the generator of the dynamics. This is then followed by a detailed analysis of σ 2 f and σ 2 f,P for some quadratic model Hamiltonians and observables.

In Section 5 we show a numerical application to illustrate the reduction in the asymptotic variance.

Preliminaries

Coordinate transformation

As mentioned above, a convenience afforded by our assumption that the preconditioner P is constant, is that a simple coordinate transformation can transform November 7, 2017 11:9 P-Langevin_Vh1 6 HA,CO the P-Langevin dynamics (1.6) into standard Langevin dynamics (1.2) by taking Z t = P 1/2 X t . We will now briefly review this transformation.

First, let us introduce the new coordinates and associated energy

z := P 1/2
x and E P N (z) := E N (x) = E N (P -1/2 z). The standard overdamped Langevin dynamics for z then reads

dZ t = -∇ z E N (P -1/2 Z t )dt + √ 2dW t = -P -1/2 ∇ x E N (P -1/2 Z t ) + √ 2dW t .
(2.1)

Upon multiplying the equation with P -1/2 we clearly recover (1.6); that is, (2.1) and (1.6) are indeed equivalent. In the new coordinate system the infinitesimal generator operator of the diffusion process Z t is given by:

L P = -∇ z E P N (z) • ∇ z + ∆ z . (2.2)
In terms of estimating observables, we obtain that

P T (f ) = 1 T T 0 f (X P t )dt = 1 T T 0 f (P -1/2 Z t )dt. (2.3)
The key observation in analysing how preconditioning changes the properties of E P N and hence of the Langevin dynamics is that ∇ 2 z E P N (z) = P -1/2 ∇ 2 x E N (x)P -1/2 . In the following, if H and P are symmetric positive definite, we will write H P := P -1/2 HP -1/2 . We also observe, for future reference, that the spectrum of H P satisfies inf σ(H P ) = inf

v T P v=1 v T Hv and sup σ(H P ) = sup v T P v=1
v T Hv.

Model Hamiltonians

Let x ∈ R N describe a system of N/d particles at positions (y α )

N/d α=1 ⊂ R d , then a simple model for potential energy is given by

E N (y) = N/d α=1 β =α φ(|y α -y β |),
where φ is e.g. a Lennard-Jones type potential, φ(r) = r -12 -2r -6 . Such systems exhibit complex meta-stable behaviour, which is an issue to be entirely separated from the ill-conditioning due to high dimension.

A much simpler situation is a mass-spring model, where u = (u α ) ∈ R N , with u α e.g. denoting out-of-plane displacement of a particle, and particle connectivity described by an equivalence relation α ∼ β. Then the energy can be written as where φ αβ could be taken as strictly convex to avoid meta-stability. We will use systems of this kind in our numerical experiments.

E N (u) = α∼β φ αβ (u α -u β ), (2.4
In order to admit explicit analytical calculations we simplify further E N by expanding it about an equilibrium (e.g., at u = 0) which then yields a quadratic energy

E N (u) = α∼β k αβ 2 |u α -u β | 2 =: 1 2 u T Hu. (2.5)
The spring constants k αβ could model how the interaction between different atomic species / environments differs. We assume throughout that there exist bounds k, k on the spring constants and a bound n on the number of neighbours, which are both independent of N . More precisely, we assume that

0 < k ≤ k αβ ≤ k and #{β : α ∼ β} ≤ n. (2.6)
Note that we have chosen a scaling of the energies where increasing system size N → ∞ does not yield a continuum limit but rather an infinite lattice system. The consequence is that if we have a bound in the spring constants 0 < k αβ ≤ k, and each atom α is connected with at most n neighbours, then it follows readily that

|E N (u)| nk |u| 2 or, equivalently H nk , (2.7)
where H denotes the 2 → 2 operator norm.

To be even more specific let us assume that E N is a d-dimensional lattice model, i.e., each atom index α corresponds to a coordinate α ∈ {0, . . . , M + 1} d and α ∼ β if and only if

| α -β | = 1, that is E N (u) := M α1=0 • • • M α d =0 d j=1 k α,α+ej u(α + e j ) -u(α) 2 .
(2.8)

By clamping the boundary sites at u = 0, we obtain N = M d free lattice sites. In this case, H is a (possibly inhomogeneous) discrete elliptic operator and employing (2.6) and using the min-max characterisation of eigenvalues (see, e.g., Sect XIII.1 in [START_REF] Reed | Methods of Modern Mathematical Physics. IV: Analysis of Operators[END_REF]) to compare H to the discrete Laplacian for which the spectrum can be computed explicitly [START_REF] Iserles | A First Course in the Numerical Analysis of Differential Equations[END_REF], we can readily show that there exist constants c 0 , c 1 such that

c 0 (j/N ) 2/d ≤ λ j ≤ c 1 (j/N ) 2/d , (2.9)
where σ(H) = {λ j | j = 1, . . . , N } is the ordered spectrum of H. This dimensiondependence of the eigenvalue distribution will be important later on.

A model preconditioner

Although for the simple model problems described in Section 2.2 it is straightforward to compute Hessians and use those as preconditioners, this would remove us from the practice of molecular simulations where Hessians are not normally computable. Instead, we will consider preconditioners that only roughly capture the structure of the energy functionals and their Hessians. Following [START_REF] Packwood | A universal preconditioner for simulating condensed phase materials[END_REF], we will use a preconditioner of the form

v T P v = c α∼β |v α -v β | 2 , (2.10)
where c is a free parameter to be fitted to the model. The idea is that it captures the connectivity information but not the fine details of the Hessian. For optimisation and saddle search, preconditioners of this kind have been shown to yield considerable speed-ups even for much more complex electronic structure type models [START_REF] Packwood | A universal preconditioner for simulating condensed phase materials[END_REF]. For example, considering E N given by (2.5) with hessian

H := ∇ 2 E N and choosing c := 1 2 (k + k), we observe that 2k k+k v T P v ≤ v T Hv ≤ 2 k k+k v T P v, (2.11)
which in particular implies that

κ H P ≤ k k .
Here, κ(A) := A A -1 denotes the condition number of the matrix A.

That is, provided the inhomogeneity k/k is not too severe, then the preconditioned energy landscape has only very moderate conditioning, independent of N , while typically κ(∇ 2 E) → ∞ as N → ∞, with a rate depending on the connectivity. For instance, for a lattice model (2.9) implies that κ(∇ 2 E N ) ∼ N 2/d .

Time Discretisation

To measure the cost/accuracy ratio of a practical sampling algorithm based on the Langevin process we need to account also for the time-discretisation,

X n+1 = X n -δP -1 ∇E N (X n ) + √ 2δP -1/2 R n , (2.12)
where

R n ∼ N (0, I N ×N ).
It is tempting to assume that one advantage afforded by preconditioning is to take larger time-steps. As we show in the following, this is a matter of scale, thus justifying our choice to focus purely on the time-continuous P-Langevin dynamics.

In the scaling that we have chosen in the model problems of Section 2.2, both the Hessians H and preconditioned Hessians H P are bounded operators (on 2 ), independent of system size N . Thus, the time-steps have similar restrictions for the preconditioned and unpreconditioned Langevin processes [START_REF] Hairer | Geometric numerical integration illustrated by the Stormer-Verlet method[END_REF]. More precisely, in the quadratic model problem (2.5) we may assume k = 1 without loss of generality. Then, using the model preconditioner (2.10) with c = k = 1, say, we have that sup σ(H) = sup σ(H P ). Proposition 1. Suppose that E N (x) = 1 2 x T Hx and P ∈ R N ×N with H and P both positive definite. Suppose further that δ 2 |P -1 H| < 1. Then the invariant measure for (2.12) is Gaussian with covariance matrix

C P δ = I -δ 2 P -1 H -1 H -1 .
Proof. See Section 6.1.

It follows that the error in the covariance operator is O(δ) with constant independence of N , and independent of the choice of preconditioner. Indeed, even with an "optimal" preconditioner P = H, one would not obtain an improvement in the bias: if there exist eigenmodes Hv = λv with λ 1, then

|C H δ v -Cv| = (1 -δ 2 ) -1 -1 |Cv|, whereas |C I δ v -Cv| = (1 -λδ 2 ) -1 -1 |Cv|, that is, |C H δ v -Cv| |C I δ v -Cv| ∼ λ -1
for δ, λ sufficiently small, where C = H -1 is the covariance operator of the unbiased measure. Therefore, for the remainder of the paper, we will not consider the effect of preconditioning on time-discretisation but only focus on the speed of convergence to equilibrium in the undiscretised (P-)Langevin dynamics. We only stress that this point of view is only valid as long as H and H P are comparable.

Exponential convergence to the invariant measure

In this section, we prove exponential convergence to the equilibrium. For the sake of simplicity we represent the probability density functions µ and µ P t by their respective densities ϕ ∞ and ϕ P t . Under the transformation z = P 1/2 x, we obtain transformed probability densities

ψ ∞ (z) := (det P ) -1/2 ϕ ∞ (P -1/2 z) and ψ P t (z) := (det P ) -1/2 ϕ P t (P -1/2 z).
Their evolution is described by the Fokker-Planck equation

∂ t ψ P t = L * P ψ P t := ∇ z • ∇ z E N (P -1/2 z)ψ P t + ∇ z ψ P t , (3.1)
where L * P is the dual operator of L P (for the L 2 (dz)-scalar product) defined in (2.2).

The proof of the following theorem is inspired from [START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF] or A.19 in [34].

Theorem 1. Suppose that E N ∈ C 2 (R N ), such that 1 2 |∇E N (x)| 2 -∆E N (x) → +∞ as |x| → +∞.
Then there exists λ P > 0 such that for all initial conditions

ψ 0 ∈ L 2 (ϕ -1 ∞ )
, and for all times t ≥ 0

ψ P t -ψ ∞ 2 L 2 (ψ -1 ∞ ) ≤ e -λ P t ψ 0 -ψ ∞ 2 L 2 (ψ -1 ∞ ) ,
or, equivalently,

ϕ P t -ϕ ∞ 2 L 2 (ϕ -1 ∞ ) ≤ e -λ P t ϕ 0 -ϕ ∞ 2 L 2 (ϕ -1 ∞ ) , ||.|| L 2 (ψ -1 ∞ ) denotes the norm in L 2 (R N , ψ -1 ∞ dz), namely ||f || 2 L 2 (ψ -1 ∞ ) = R N f 2 (z)ψ -1 ∞ dz.
The exponent λ P is the spectral gap of the Fokker-Planck operator L * P defined in (3.1) (i.e., the smallest non-zero eigenvalue of -L * P ).

Proof. See Section 6.2.

We now consider the quadratic case, E N (x) = 1 2 x T Hx with H symmetric and positive definite. In this case, we have the following characterisation of the spectrum of L * P and hence of λ P . The proof is based on [START_REF] Hitrik | Spectra and semigroup smoothing for non-elliptic quadratic operators[END_REF][START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Ottobre | Exponential return to equilibrium for hypoelliptic quadratic systems[END_REF].

Theorem 2. The spectrum of the operator

L * P is σ(L * P ) = σ(L P ) = -λ∈σ(H P ) k λ λ : k λ ∈ N . (3.2)
Proof. See Section 6.4.

An immediate consequence of Theorem 2 is that

λ P = inf σ(H P ) \ {0}, (3.3)
in other words, the smallest non-zero generalised eigenvalue of

Hv = λP v.
Returning to the quadratic model Hamiltonians and model preconditioners introduced in Sections 2.2 and 2.3, assume that min σ(H) ∼ N -s (which is consistent with (2.9)), while min σ(H P ) ∼ 1, then we obtain that λ I ∼ N -s with λ P ∼ 1. This result should give us confidence in the value of preconditioning.

Remark 1. The preconditioning ideas and results presented in Sections 2 and 3 are similar to the Brascamp-Lieb inequality [START_REF] Brascamp | On extensions of the Brunn-Minkowski and Prékopa-Leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation[END_REF]. In some sense, this inequality claims that a good preconditioner is the Hessian.

Precisely [20], if E N is strictly convex, for any function

f ∈ L 2 (e -E N ), f -f e -E N 2 e -E N ≤ ∇f (∇ 2 E N ) -1 • ∇f e -E N ,
where we assume here that the normalization e -E N = 1. This means that if one considers the Fokker-Planck equation

∂ t ψ = div[(∇ 2 E N ) -1 e -E N ∇(ψe E N )],
which is associated to the following overdamped Langevin dynamics [START_REF] Xifara | Langevin diffusions and the Metropolis-adjusted Langevin algorithm[END_REF]: then, if E N is strictly convex, we have

dX t = -(∇ 2 E N ) -1 ∇E N (X t )dt + div[(∇ 2 E N ) -1 ](X t )dt + 2(∇ 2 E N ) -1 (X t )dW t , November
1 2 d dt ψ e -E N -1 2 e -E N = -(∇ 2 E N ) -1 ∇ ψ e -E N • ∇ ψ e -E N e -E N ≤ ψ e -E N -1 2 e -E N .
And thus

ψ e -E N -1 2 e -E N ≤ ψ 0 e -E N -1 2 e -E N e -2t ,
whatever the potential E N is. It is in some sense "universal". For example, it does not depend on the temperature: if we multiply E N by a constant β (inverse of the temperature) it remains the same, whereas understanding the dependency of the spectral gap on the temperature is tricky in general.

Analysis of the asymptotic variance

In this section we present sufficient conditions under which the estimator

P T (f ) = 1 T T 0 f (X P t )dt = 1 T T 0 f (P -1/2 Z t )dt
satisfies a central limit theorem of the form (1.8) and we characterise the associated asymptotic variance.

Generalities

The fundamental requirements to prove the central limit theorem is establishing the well-posedness of the Poisson equation

-L P φ(z) = f (P -1/2 z) -µ(f ), µ(φ) = 0, (4.1)
for all bounded and continuous functions f : R N → R, where L P is defined by (2.2), and obtaining estimates on the growth of the unique solution φ. Recall that we shall assume that µ admits a smooth, strictly positive density denoted by ψ ∞ (z), such that R N ψ ∞ (z)dz = 1 and the SDE (2.1) has a unique strong solution. Referring to results in [START_REF] Glynn | A Liapounov bound for solutions of the Poisson equation[END_REF][START_REF] Meyn | A survey of Foster-Lyapunov techniques for general state space Markov processes[END_REF] we suppose that the process Z t admits a Lyapunov function (see the Definition 1 in Section 6.5), which is sufficient to ensure the geometric ergodicity of Z t (see [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Meyn | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF]). In terms of the potential energy E N and the preconditioner P , we require that there exists β ∈ (0, 1) such that

lim |z|→+∞ inf (1 -β)|∇E P N (z)| 2 + ∆E P N (z) > 0. (4.2)
It is straightforward to check that this condition holds whenever E N is strongly convex and in particular if it is of the form (2.5).

If condition (4.2) holds, then the process Z t will be geometrically ergodic. More specifically, the law of the process Z t started from a point z ∈ R N will converge exponentially fast in the total-variation norm to the equilibrium distribution µ (cf.

(1.7)).

Assuming (4.2), we also obtain the following well-posedness result for the Poisson equation (4.1).

Theorem 3. Suppose that (4.2) holds, then there exists c > 0, such that for any measurable observable f satisfying |f | 2 ≤ e -βE P N (z) , the Poisson equation (4.1) admits a unique strong solution satisfying the bound |φ(z)| 2 ≤ ce -βE P N (z) . In particular, φ(P 1/2 x) ∈ L 2 (µ).

Proof. See Section 6.5.

The technique of using a Poisson equation to obtain a central limit theorem for an additive functional of a Markov process is widely known (see e.g. [START_REF] Bhattacharya | A central limit theorem for diffusions with periodic coefficients[END_REF]). For linear and quadratic observables, we can in fact produce an analytic solution to this Poisson problem; see Section 4.2 below. Theorem 4. Under the conditions of Theorem 3, there exists a constant 0 < σ 2 f,P < ∞ such that the asymptotic distribution of the fluctuations of P t (f ) about µ(f ) are given by the central limit theorem

√ t P t (f ) -µ(f ) D -→ N (0, σ 2 f,P ), as t → +∞, (4.3)
where σ 2 f,P (the asymptotic variance for the observable f ) is given by

σ 2 f,P = 2 ∇ z φ(z) 2 µ P (dz) = 2 ∇ z φ(P 1/2 x) 2 µ(dx), (4.4)
where µ P (dz) = ψ ∞ (z)dz.

Proof. See Section 6.6.

Explicit solution for linear observables

In this section we exhibit explicit solutions φ of the Poisson equation (4.1) when E N is quadratic and f is linear, and compute the associated asymptotic variance σ 2 f,P . This simplest possible case is of course well-known but we summarise it nevertheless to prepare for more interesting cases.

Suppose, therefore, that

P = I, E N (x) = 1 2 x T Hx and f (x) = v • x,
where v ∈ R N . From symmetry it follows that µ(f ) = 0, hence the Poisson equation Seeking a solution of the form φ

(x) = d • x with d ∈ R N , we obtain d = H -1 v, i.e., φ(x) = x • H -1 v,
and in particular,

σ 2 f = σ 2 f,I = 2 |H -1 v| 2 µ(dx) = 2|H -1 v| 2 .
In particular, choosing v to be a normalised eigenmode of H with associated eigenvalue λ, we obtain σ 2 f = 2λ -2 . Focusing specifically on a d-dimensional lattice model, we know from (2.9) that min σ(H) ∼ N -2/d while max σ(H) ∼ 1. Hence, it follows that σ 2 f is moderate for the high-frequency eigenmodes, but large for the observables corresponding to low-frequency eigenmodes.

Next, we turn to the preconditioned dynamics. In this case we effectively replace H with H P = P -1/2 HP -1/2 and v with P -1/2 v and thus obtain

σ 2 f,P = 2 |H -1 P P -1/2 v| 2 µ(dx) = 2|P 1/2 H -1 v| 2 =: 2|H -1 v| 2 P . (4.6)
If we assume that c 0 = min σ(H P ), c 1 = max σ(H P ), then a simple rewrite yields

2c 0 v T H -1 v ≤ σ 2 f,P ≤ 2c 1 v T H -1 v. (4.7)
Comparing with σ 2 f = 2v T H -2 v and recalling our standing assumption (2.11) that c 0 ∼ 1, c 1 ∼ 1 (the spectrum is bounded above and below independently of N ), we conclude that preconditioning does not entirely remove ill-conditioning but it is (potentially) diminished.

More concretely, for a d-dimensional lattice system, we obtain that

N -2/d σ 2 f,P σ 2 f 1,
and both bounds are attained for specific observables. We conclude that preconditioned Langevin can be significantly more efficient than standard Langevin (lowfrequency observables) but that it will be comparable in efficiency for high-frequency observables.

Intuitively, low-frequency observables are "macroscopic" in nature and include e.g. energy, average bond-length etc., while high-frequency observables include in particular single bonds, bond angles and dihedral angles (in a large molecule) or a bond-length near a crack-tip. In the next sections, we consider three toy models mimicking "realistic" observables of these kinds, occurring in real-world simulations, to further substantiate our remarks.

Example 1: Energy per particle

We now consider E N (x) = 1 2 x T Hx and f (x) = N -1 E N (x). A straightforward computation yields

x • Bx µ = Tr(H -1 B) for B ∈ R N ×N , (4.8)
November 7, 2017 11:9 P-Langevin_Vh1

14 HA,CO which in particular implies that

µ(f ) = 1 N E N (x)e -E N (x) e -E N (x) = TrI 2N = 1 2 .
Thus, the Poisson equation becomes

Hx • ∇φ(x) -∆φ(x) = 1 2N x T Hx -1 2
We seek a solution of the form φ(x) = 1 2 x T Bx+l•x-TrB (to ensure that µ(φ) = 0), then this yields the equation

Hx • (Bx + l) -TrB = 1 2N x T Hx -1 2 .
This is satisfied for

B = 1 N I, l = 0, hence φ(x) = 1 2N |x| 2 -1 2 .
We can now compute the asymptotic variance as

σ 2 f = σ 2 f,I = 2 1 N x 2 µ(dx) = 2 N 2 TrH -1 .
Repeating the same argument in transformed coordinates z = P 1/2 x, we also obtain the asymptotic variance of the energy for the preconditioned Langevin dynamics:

σ 2 f,P = 2 N 2 TrH -1 P = 2 N 2 Tr H -1 P .
Let us now focus on a lattice model, where we have (2.9). Then we obtain that

1 2 σ 2 f ≈ N -2 N j=1 (j/N ) -2/d ≈ N -1 1 1/N s -2/d ds ≈        1, d = 1, N -1 log N, d = 2, N -1 , d = 3, while, clearly, 1 2 σ 2 f,P ≈ N -1 . In summary, σ 2 f σ 2 f,P ≈        N, d = 1, log N, d = 2, 1, d = 3;
that is, preconditioning only gives a significant speed-up in one-dimensional lattices but not in two-or three-dimensional lattices.

Example 2: Bond-length

In our second example we observe a single bond in the crystal or molecule. That is, we still use E N (x) = 1 2 x T Hx but the observable is now given by f (x) = x i -x j for some fixed bond i ∼ j. This is a linear observable, hence a special case of the discussion in Section 4.2. Hence, we obtain

σ 2 f = 2|H -1 l| 2 where l n =        1, n = i -1, n = j, 0, otherwise.
In order to estimate σ 2 f further we consider again the d-dimensional lattice model (2.8) and P given by (2.10). For d ≥ 2, since P is a homogeneous discrete elliptic operator, we know from [27] that

[P -1 ] ni -[P -1 ] nj ≤ C(1 + |n -i|) -d ,
where we note that i, j are now neighbouring lattice sites; i.e., [P -1 ] ni -[P -1 ] nj denotes a discrete gradient of the lattice Green's function.

Therefore, we obtain that

σ 2 f = 2|H -1 l| 2 2|P -1 l| 2 2 n∈Z d (1 + |n -i|) -2d < ∞; (4.9)
that is, σ 2 f has an upper bound that is independent of N . A lower bound follows simply from the fact that H ≤ 1 and hence v T H -1 v ≥ |v| 2 , which implies

σ 2 f ≥ |l| 2 = 1.
To obtain bounds on σ 2 f,P , we use (4.7) to estimate

2c 0 = 2c 0 |l| 2 ≤ σ 2 f,P ≤ 2c 1 l T H -1 ≤ 2c 1 |l||H -1 l|,
and we have already shown in (4.9) that this is bounded above, independently of N .

In summary, we obtain that

σ 2 f σ 2 f,P 1 for d ≥ 2,
that is, we expect no substantial (if any) speed-up for the bond-length observable from preconditioning for d ≥ 2. By contrast, for d = 1, the system P -1 l can be solved explicitly, and in this case one obtains

σ 2 f ∼ N as N → ∞ (specifically, |P -1 l| 2 = N/12), that is, σ 2 f σ 2 f,P ∼ N for d = 1.
Thus, we conclude that preconditioning helps to accelerate the computation of the bond-length observable only for one-dimensional structures. The intuitive explanation of this effect is that far-away regions of space have little influence on a single bond and hence only local equilibration matters. The difference in 1D is that elastic interaction is naturally more long-ranged than in dimension d > 1. Our final example is inspired by a technique called umbrella sampling [START_REF] Kästner | Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: Umbrella integration[END_REF]. Given a potential energy E N (x) and a reaction coordinate ξ(x), we wish to compute

A(ξ 0 ) := -log e -E N (x) δ(ξ(x) -ξ 0 ) dx.
Umbrella sampling achieves this by placing a restraint on the potential energy,

E N,K (x) := E N (x) + K 2 ξ(x) -ξ 0 2 ,
for some K > 0, ξ 0 ∈ R. Let µ K denote the corresponding equilibrium measure, then it can be shown [START_REF] Kästner | Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: Umbrella integration[END_REF] that, for K large,

∂A ∂ξ ξ= ξ0 ≈ -K( ξ0 -ξ 0 ) where ξ0 = ξ µ K .
Thus, ∂ ξ A and hence A can be reconstructed in this way. More sophisticated variations of the idea exist of course, but for the sake of simplicity of presentation we will focus on this particularly simple variant.

To construct an analytically accessible toy problem mimicking umbrella sampling we consider again a quadratic energy E N (x) = 1 2 x T Hx and a linear reaction coordinate ξ(x) := l • x where, for simplicity, we assume that |l| = 1 (for |l| = O(1), the argument is analogous). The restrained potential with penalty parameter K > 0 is then given by E N,K (x) = 1 2 x T Hx + K 2 (ξ(x) -ξ 0 ) 2 , for some ξ 0 ∈ R, while the observable from which we can reconstruct the mean force is simply

f (x) = K l • x -ξ 0 .
We are again in the context of Section 4.2 and therefore obtain

σ 2 f = 2K 2 H -1 K l 2 ,
where

H K = ∇ 2 E N,K = H + Kll T .
The Sherman-Morrison formula yields

H -1 K l = H -1 - KH -1 ll T H -1 1 + Kl T H -1 l l = H -1 l 1 + Kl T H -1 l , (4.10)
and hence,

σ 2 f = 2K 2 H -1 K l 2 = 2K 2 |H -1 l| 2 (1 + Kl T H -1 l) 2 .
Since our focus in the present example is the ill-conditioning induced by large K rather than ill-conditioning induced by H (e.g. through system size N ), let us assume that max σ(H) = 1 (as always) while min σ(H) ≥ c 0 , for some moderate pol 17 constant c 0 . This would, e.g., be the typical situation for a small molecule, or if we preconditioned H but without accounting for the umbrella. We then obtain

c 2 0 ≤ 2K 2 c 2 0 (1 + K) 2 ≤ σ 2 f ≤ 2K 2 (1 + Kc 0 ) 2 ≤ 2 c 2 0 ; that is, 1 2 σ 2 f ∼ 1 as K → ∞.
By contrast, suppose now that we choose a preconditioner

P K := P + Kll T ,
where P is a preconditioner for H satisfying c P 0 v T P v ≤ v T Hv ≤ v T P v. Then a straightforward calculation yields

c P 0 ≤ v T Hv + K(l • v) 2 v T P v + K(l • v) 2 = v T H K v v T P K v ≤ 1.
It follows from (4.7) that

σ 2 f,P = 2K 2 H -1 K l 2 P K ≈ 2K 2 |H -1/2 K l| 2 ,
where ≈ now indicates upper and lower bounds with constants that are independent of K. Using (4.10) we obtain

σ 2 f,P ≈ 2K 2 l T H -1 l 1 + Kl T H -1 l ∼ K as K → ∞.
We can therefore conclude that

σ 2 f,P σ 2 f ≈ (l T H -1 l)(1 + Kl T H -1 l) |H -1 l| 2 ∼ 1 + K as K → ∞;
that is, preconditioning the umbrella actually achieves a significant deterioration of the asymptotic variance and thus the P -Langevin dynamics actually becomes less efficient than the standard Langevin dynamics. However note the following crucial remark:

Remark 2 (Step-sizes revisited). The surprising result of the present section does not in fact fully fall within our starting assumptions. While H = 1, H K is in fact of order O(1 + K) which means that the time-step for the discretisation of the Langevin equation should be of order O(K -1 ), which exactly balances the lower mixing of the preconditioned dynamics and make the two schemes again comparable. Indeed, in our computational examples we will need to choose ∆t = O(K -1 ) to prevent instability. In practice, the restraint parameter K is chosen of the same order of magnitude of the stiffest bond in a molecule, while the reaction coordinate will normally be a function of the softest bonds and hence it would create no additional timestep restriction. In such a situation, it is indeed preferable to not precondition the restraint.

November 7, 2017 11:9 P-Langevin_Vh1 18 HA,CO However, we emphasise again that the interaction between preconditioning and time-stepping is an issue that we do not properly address in the present work and which will require further attention in the future.

Numerical Tests

We conclude our discussion by demonstrating the extension of our explicit computations to a mildly non-linear lattice model. As potential energy E M d : R {1,...,M } d → R, we choose

E N (u) := M α1=0 • • • M α d =0 d j=1 φ 1 √ d u(α + e j ) -u(α) ,
where u(α) := 0 if any α j ∈ {0, M + 1}, and with convex nearest-neighbour pair potential φ(r) = 1 8 r 2 + sin(r) 2 . Upon choosing an arbitrary linear labelling of indices α ∈ {1, . . . , M } d , this is a special case of (2.4). As preconditioner we choose (2.10), which we can also write as

P v, v = 1 4d M α1=0 • • • M α d =0 d j=1 |u(α + e j ) -u(α)| 2 .
The occurrence of d in the definitions of E N and P ensures that ∇E N ≈ P ≈ 1; moreover, since φ is strictly convex, we have that cond((∇E N (x)) P ) is bounded above independently of N ; cf. (2.11). For all simulations (with small modifications for the umbrella sampling example) we choose δ = 0.1, N steps = 10 5 , N runs = 400.

We then use a Cholesky factorisation to compute P , i.e. P = LL T , followed by

X P n+1 = X P n -δP -1 ∇E N (X P N ) + √ 2δL -T R n , for n = 1, . . . , N steps ,
where R n ∼ N (0, I) N . The estimate of the observable value is then given by f = N -1 steps Nsteps n=1 f (X P n ). We compute N runs trajectories in order to estimate the asymptotic variance from N runs independent samples of f . For the umbrella sampling test, we choose ξ(x) := l • x to be the bond-length observable again with ξ 0 = 0.33. The modification to the energy and observable is then as described in Section 4.5. With δ = 0.1, the discretised Langevin dynamics turns out to be unstable, hence we had to choose δ K = 0.1/K instead. To account for this (see also Remark 2), we study Kσ 2 f instead of σ 2 f in our tests. We perform the umbrella sampling test only for d = 2, N = 8 2 , since we focus here on the magnitude on the restraint parameter K rather than the system size. The results of the simulations are shown in table 1. The numbers closely match the analytical predictions of Sections 4.3, 4.4 and 4.5. To conclude this discussion we only remark that we did not fine-tune step-sizes, which means one could likely make small improvements to both the preconditioned and unpreconditioned processes. However, we believe that the trends across dimension, system size and restraint parameter are reliable.

Asymptotic Variance: Energy

d = 1 d = 2 d = 3 N 1 d σ 2 f σ 2 f,P N 1 d σ 2 f σ 2 f,P N 1 d σ 2 f σ 2 f,P 8
Asymptotic Variance: Bond-length d = 1 d = 2 d = 3 N 1 d σ 2 f σ 2 f,P N 1 d σ 2 f σ 2 f,P N 1 d σ 2 f σ 2 f,P 8
In particular, we stress that even though the preconditioned variants often have a smaller asymptotic variance, often (in particular for d = 3) this improvement is only by a moderate constant factor. Because only the trends are reliable indicators in these tests only a successive improvement with increasing N or K (e.g. as in the d = 1 tests) can be considered a success for the preconditioned algorithm.

Conclusion

In this paper we strived to develop an intuition what the effect of preconditioning has on molecular simulations. The results are very mixed: it is clear that preconditioning accelerates convergence of the probability density functions to equilibrium (see Theorem 1 as well as the discussion in Section 4.2), and this necessarily implies accelerated convergence for some observables. However, for many concrete observables of practical importance little (if any) benefit can be gained. This was a surprising outcome for us and indicates that alternative avenues need to be explored on how a priori information about the analytical structure of configuration space should be exploited in molecular simulation.

We emphasize again that our (partially negative) conclusion, contrary to much of the existing literature, is due to the fact that we test the convergence of specific observables. Moreover, we stress that we have only performed a limited set of tests on highly simplified toy models and a limited set of observables, while more realistic models may exhibit many features that we neglected. Proof. The covariance of the invariant measure associated to the dynamics (2.12) is given by the following identity:

C P δ = (I -δP -1 H)C P δ (I -δHP -1 ) + 2δP -1 C P δ = C P δ -δP -1 HC P δ -δC P δ HP -1 + δ 2 P -1 HC P δ HP -1 + 2δP -1 P -1 HC P δ + C P δ HP -1 = 2P -1 + δP -1 HC P δ HP -1 . Expanding C P δ , we have

C P δ ∼ C 0 + δC 1 + δ 2 C 2 + O(δ 3
), then one gets

P -1 HC 0 + C 0 HP -1 = 2P -1 ⇒ C 0 = H -1 P -1 HC 1 + C 1 HP -1 = P -1 HP -1 ⇒ C 1 = 1 2 P -1 P -1 HC 2 + C 2 HP -1 = 1 2 P -1 HP -1 HP -1 ⇒ C 2 = 1 4 P -1 HP -1 .
Proceeding by induction, one can therefore obtain:

∀k ∈ N, C k = 2 -k (P -1 H) k H -1 .
Therefore C P δ can be rewritten as:

C P δ = +∞ k=0 δ 2 P -1 H k H -1 .
November 7, 2017 11:9

P-Langevin_Vh1 pol 21
Indeed, the sum k≥0 C k converges since δ 2 |P -1 H| < 1. The covariance operator be can rewritten as

C P δ = ∞ k=0 δ 2 P -1 H k H -1 = I -δ 2 P -1 H -1 H -1
= (H -δ 2 HP -1 H) -1 , which concludes the proof.

Proof of Theorem 1

Proof. Under the assumptions on the potential E N , see A.19 in [34], the density ψ ∞ satisfies a Poincaré inequality: there exists λ P > 0 such that for all probability density functions φ, we have:

R N φ ψ ∞ -1 2 ψ ∞ dz ≤ 1 λ P R N ∇ φ ψ ∞ 2 ψ ∞ dz. (6.1)
The optimal parameter λ P in (6.1) is the opposite of the smallest (in absolute value) non-zero eigenvalue of the Fokker-Planck operator L * P , which is self-adjoint in L 2 (R N , ψ -1 ∞ dx). Thus the exponent λ P is the spectral gap of L * P . If ψ P t is a solution of (3.1), therefore for all initial condition ψ 0 ∈ L 2 (ψ -1 ∞ ), ∀t ≥ 0 :

d dt ψ P t -ψ ∞ 2 L 2 (ψ -1 ∞ ) = -2 R N ∇ ψ P t ψ ∞ 2 ψ ∞ dz. (6.2) Indeed, d dt ψ P t -ψ ∞ 2 L 2 (ψ -1 ∞ ) = d dt ψ P t -ψ ∞ 2 ψ -1 ∞ dz = 2 ∂ t ψ P t ψ P t -ψ ∞ ψ -1 ∞ dz = 2 ∇ • ∇E P N (z)ψ P t + ∇ψ P t ψ P t ψ ∞ -1 dz = -2 ∇E P N (z)ψ P t + ∇ψ P t • ∇ ψ P t ψ ∞ dz.
But we have which yields (6.2). Therefore, using (6.1),

∇E P N (z)ψ P t + ∇ψ P t = -∇ (ln (ψ ∞ )) ψ P t + ∇ψ P t = - ∇ψ ∞ ψ P t ψ ∞ + ∇ψ P t = ∇ ψ P t ψ ∞ ψ ∞ , November
d dt ψ P t -ψ ∞ 2 L 2 (ψ -1 ∞ ) ≤ -2λ P ψ P t ψ ∞ -1 2 ψ ∞ dz = -2λ P ψ P t -ψ ∞ 2 ψ -1 ∞ dz, then ψ P t -ψ ∞ 2 L 2 (ψ -1 ∞ ) ≤ e -λ P ψ 0 -ψ ∞ 2 L 2 (ψ - 1
∞ ) .

Preliminaries for spectral analysis

In the linear case, i.e E P N (z) = 1 2 z T H P z, where H and P are symmetric positive definite and H P := P -1/2 HP -1/2 , the analysis will be carried out in a suitable system of coordinates which simplifies the calculations and the proofs of the main theorems. For this reason, we will perform one conjugation and one additional change of variables .

From the partial differential equation point of view and in order to use standard techniques from the spectral analysis of partial differential equations, then it appears to be useful to work in L

2 (R N , dz; C) instead of L 2 (R N , ψ ∞ dz; C). The mapping φ → ψ -1/2 ∞ φ maps unitarily L 2 (R N , dz; C) into L 2 (R N ,
= ∇ z + 1 2 ∇ z H P .
Thus, the operator

L P = -∇ z E P N (z) • ∇ z + ∆ z is transformed into L P = e -1 2 H P L P e 1 2 H P = ∆ z - 1 4 |∇ z E P N (z)| 2 + 1 2 ∆ z E P N (z) = ∆ z - 1 4 z T H 2 P z + 1 2
Tr(H P ).

The kernel of L P is Ce -z T H P z 4 and the operator L P is unitarily equivalent to the operator L P .

In the goal of modifying the kernel of the operator L P into a centered Gaussian with identity covariance matrix, we perform a second change of variables. In the following, we introduce the new coordinates y = H . The last conjugation and change of variables are used to compute the spectrum of L * P needed to proof Theorem 2 (see Section 6.4). Let us now introduce some additional notations. Recall the space of rapidly decaying complex valued C ∞ functions

S(R N ) = f ∈ C ∞ (R N ), ∀α, β ∈ N N , ∃C α,β ∈ R + , sup x∈R N |x α ∂ β x f (x)| ≤ C α,β ,
and its dual is denoted S (R N ).

The Weyl-quantization q W (x, D x ) of a symbol q(x, ξ) ∈ S (R N ) is an operator defined by its Schwartz-kernel

q W (x, D x ) (x, y) = R N e i(x-y)•ξ q x + y 2 , ξ dξ (2π) N .
For instance, the Weyl symbol of of the operator

-LP + 1 2 Tr(H P ) = -∇ T y H P ∇ y + 1 4 y T H P y is q(y, ξ) = ξ T H P ξ + y T H P y 4 . (6.5)
Those tools are essential to proof Theorem 2 (see Section 6.4). For more details on Weyl-quantization, one can refer to [START_REF] Ottobre | Exponential return to equilibrium for hypoelliptic quadratic systems[END_REF].

Proof of Theorem 2

Proof. Referring to Theorem 1.2.2 in [START_REF] Hitrik | Spectra and semigroup smoothing for non-elliptic quadratic operators[END_REF], the spectrum of the operator q W (y, D y ) = -LP + 1 2 Tr(H P ) associated with the elliptic quadratic Weyl symbol q(y, ξ) defined by (6.5) is given by

σ(q W (y, D y )) =        λ∈σ(G) Imλ≥0 -iλ(r λ + 2k λ ), k λ ∈ N        .
where G is the so-called Hamilton map associated with q , and r λ is the algebraic multiplicity of λ ∈ σ(G) (the dimension of the characteristic space). The Hamilton map is the C-linear map G : C 2N → C 2N associated with the matrix

G = 0 H P -1 4 H P 0 ∈ C 2N ×2N .
The matrix G is similar to another matrix denoted G and defined by -iλr λ = Tr(H P ) 2 , which concludes the proof of the theorem. and U (z) ≥ 1, z ∈ R N , where 1 C is the indication function over a petite Borel subset C of R N (refer to [START_REF] Meyn | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF] for more details).

G = 1 √ 2 0 0 √ 2 G √ 2 0 0 1 √ 2 = 1
For the generator L P corresponding to (2.2), compact sets are always petite. In the following we prove that the Foster-Lyapunov criterion holds for (2.1). But first we need an assumption on the potential E P N . Assumption A There exists k > 0 such that µ P (z) is bounded from above for all |z| ≥ k and, for some 0 < β < 1, Proof. Recall the generator of (2.1)

L P = -∇ z E P N (z)∇ + ∆.
For U (z) = e βE P N (z) , one obtains: Finally, since ψ ∞ (z) is bounded, then U (z) is bounded away from zero uniformly. Then U (z) can be rescaled to satisfy the condition U (z) ≥ 1. Thus, we valid the Foster-Lyapunov criterion for (2.1).

L P U (z) = -∇U (
Using this lemma, the proof of the well-posedness result for the Poisson equation (4.1) come straightforward using Theorem 3.2 in [START_REF] Bhattacharya | A central limit theorem for diffusions with periodic coefficients[END_REF].

Proof of Theorem 4

Proof. We start the proof by decomposing µ P t (f ) -µ(f ) into a martingale and a remainder terms: Consider now the rescaling √ t( P t (f ) -µ(f )). Using the central limit theorem for the martingale term √ tM t (see [START_REF] Helland | Central limit theorems for martingales with discrete or continuous time[END_REF] , Theorem 5.3), one obtains the following convergence in distribution It remains to study the remainder term √ tR t . We consider the two cases: If Z 0 ∼ µ, then since φ ∈ L 2 (µ), we have that √ tR t D -→ 0 in L 2 (µ).

In the more general case, we must refer to a "propagation of chaos" argument (see for example [START_REF] Cattiaux | Central limit theorems for additive functionals of ergodic Markov diffusions processes[END_REF], Section 8 and [START_REF] Duncan | Variance Reduction using Nonreversible Langevin Samplers[END_REF]), to obtain the same result.

The result of this test is shown in Figure 1, where we plotted the RMS error of the reconstructed profiles based on 10 independent MTD simulations for each parameter set. We observe that in general, preconditioning does not improve the convergence.

As a matter of fact, we performed a variety of similar tests, e.g. with different parameters, or different observables, all of which led to similar conclusions. Indeed, depending on the choice of observable, preconditioning often has an even larger negative impact.

γFigure 1 .

 1 Figure 1. RMS errors of the reconstructed free energy surface profiles from metadynamics simulations without (black dashed line) and with preconditioning (coloured solid lines). Error bars represent one standard deviation based on 10 independent simulations. See Appendix A for the detailed setup.

  pol 9

( 4 .

 4 1) becomesHx • ∇φ(x) -∆φ(x) = v • x.

5 .

 5 Example 3: Umbrella sampling

1 .

 1 Proof of Proposition 1

6. 5 . 3 Definition 1 .

 531 Proof of Theorem (Foster-Lyapunov criterion) We say that the Foster-Lyapunov criterion holds for (2.1) if there exists a function U : R N → R and constants C > 0 and b ∈ R such that µ(U ) < ∞, L P U (z) ≤ -cU (z) + b1 C(6.6)

  lim |z|→+∞ inf (1 -β)|∇E P N (z)| 2 + ∆E P N (z) > 0.

pol 25 Lemma 1 .

 251 Under Assumption A, the Foster-Lyapunov criterion holds for (2.1) with:U (z) = e βE P N (z) , 0 < β < 1.(6.7)

  z) • ∇E P N (z) + ∆U (z) = -∇e βE P N (z) • ∇E P N (z) + ∇ • (∇e βE P N (z) ) = -β|∇E P N (z)| 2 e βE P N (z) + β∆E P N (z)e βE P N (z) + β 2 |∇E P N (z)| 2 e βE P N (z) = -β (1 -β)|∇E P N (z)| 2 + ∆E P N (z) U (z).(6.8)Therefore, by Assumption A, for ε > 0, ∃k > 0 such that ∀ |z| > k:(1 -β)|∇E P N (z)| 2 + ∆E P N (z) > ε,and so also,L P U (z) ≤ -βεU (z) + b1 C k ,where C k = {z ∈ R N ; |z| ≤ k} and b > 0.

1 t t 0 f 1 t t 0 -

 1010 -1/2 Z s )ds -µ(f ) = (P -1/2 Z s ) -µ(f ) ds = L P φ(Z s )ds = φ(Z 0 ) -φ(Z t ) s )dW s := R t + M t .

= 2

 2 |∇ z φ(z)| 2 µ P (dz) = 2 |∇ z φ(z)| 2 µ(dx).

Table 1 .

 1 Numerically estimated asymptotic variances of the energy observable (Section 4.3), the bond-length observable (Section 4.4) and the restraint observable occurring in umbrella sampling (Section 4.5). The nonlinear potential energy used in these tests is described in Section 5. All results match the analytical predictions of Sections 4.3, 4.4 and 4.5.

	2.4e+01 7.4e+00	4	3.2e+01 8.5e+00	4	1.9e+01 7.5e+00
	16 4.2e+01 7.2e+00	8	5.4e+01 9.5e+00	6	1.8e+01 7.6e+00
	32 7.6e+01 7.4e+00	16 7.2e+01 1.1e+01	9	1.9e+01 7.5e+00
	128 2.8e+02 7.0e+00	32 1.1e+02 1.1e+01	13 2.2e+01 7.7e+00
	Asymptotic Variance: Umbrella Sampling
			d = 2, N 1/2 = 8
		K	Kσ 2 f	σ 2 f,P
		10 1.6e+01 1.8e+01
		20 3.5e+01 4.0e+01
		40 6.7e+01 7.7e+01
		80 1.3e+02 1.8e+02

  ψ ∞ dz; C) with the associated transformation rules for the differential operators:

	e -1 2 H P ∇ z e	1 2 H P

  Now, the characteristic polynomial of G can be computed bydet(G -λI) = det(G -λI) = 2 -2N -2λI H P -H P -2λI -2N det(H P -i2λI)det(H P + i2λI).Since Re(σ(H P )) ≥ 0, one thus obtains that σ(G) ∩ {λ, Imλ ≥ 0} = i 2 σ(H P ).

	= 2 -2N	-2λI -H P -i2λI i(H P + i2λI) H P
	= 2 -2N i(-H P + i2λI) 0	H P i(H P + i2λI)
	= 2 In particular,
	-iλ2k λ =	k i 2 µ µ
	λ∈σ(G)	µ∈σ(H P )
	Imλ≥0
	and
	λ∈σ(G)
	Imλ≥0
			2	0 H P -H P 0	.

HA,CO

Acknowledgments

We thank Tony Lelièvre and Jonathan Weare for helpful discussions. The work was supported by ERC Starting Grant 335120. November 7, 2017 11:9 P-Langevin_Vh1 pol 27

Appendix A. Setup of the Metadynamics example

In the example shown in Figure 1, we investigated how preconditioning the second order Langevin dynamics affects the sampling of different observables when applied to adaptive potential of mean force (PMF) techniques such metadynamics (MTD) [START_REF] Laio | Escaping free-energy minima[END_REF]. Unlike the previously discussed umbrella sampling method, adaptive PMF techniques build their bias and mean force or potential of mean force estimates on-the-fly during the dynamics. In a general form, the corresponding biased and preconditioned equations of motion with preconditioner P = P (q) can be written as

where M ∈ R N ×N is the diagonal mass matrix, β is the inverse temperature, ξ is the collective variable (reaction coordinate) and F b (ξ, t) is the biasing force. In the case of metadynamics,

is some history dependent biasing potential composed by Gaussians regularly deposited in the collective variable space:

where δ t is the height of the corresponding Gaussian and w is a diagonal matrix including the widths of the collective variable components.

The molecular system we chose for this test was the 2-(formylamino) propionaldehyde in gas phase using the Amber99SB [START_REF] Hornak | Comparison of multiple Amber force fields and development of improved protein backbone parameters[END_REF] force field.

The molecule has a single slow degree of freedom, a dihedral angle, that was selected as the one-dimensional collective variable (we note that this variable is highly associated to one of the two dihedrals of alanine dipeptide, a test system widely used in the computational chemistry field).

We performed 200 ps long molecular dynamics simulations at 300 K using the BBK integrator scheme [START_REF] Brünger | Stochastic boundary conditions for molecular dynamics simulations of ST2 water[END_REF] with 0.5 fs timestep. In the case of unpreconditioned dynamics P = γM was used with γ = 5.0 ps -1 , while for the preconditioned dynamics we applied P = γM + τ H, where H is a Hessian-based preconditioner, whose positiveness is guaranteed by rebuilding the matrix using the spectral decomposition of the Hessian with the absolute values of the eigenvalues [START_REF] Mones | Regularity and locality of point defects in multilattices, to appear in AMRX[END_REF]. We used γM as diagonal a stabiliser and varied γ and τ parameters.

A deposition frequency of 1/50 fs -1 and a starting height of δ = 0.004 eV were used in a well-tempered variant of MTD [START_REF] Barducci | Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method[END_REF] with T w = 10000 K.

Free energy profiles were reconstructed simply as the negative of the actual history dependent biasing potential. The reference for computing the RMS error of the profiles was obtained from a 200 ns long unpreconditioned constrained dynamics simulation [START_REF] Mones | Exploration, Sampling, And Reconstruction of Free Energy Surfaces with Gaussian Process Regression[END_REF].


Internet-Draft Google Intended status: Standards Track

This document presents the motivation and requirements for the use of Network Level Packet Erasure Coding to improve the performance of the QUIC protocol that is proposed a new transport protocol. The document does not specify a specific code but lists the salient features that a code should have in order to deal with know loss patterns on QUIC paths.

Terminology

The key words "MUST", "MUST NOT", "REQUIRED", "SHALL", "SHALL NOT", "SHOULD", "SHOULD NOT", "RECOMMENDED", "MAY", and "OPTIONAL" in this document are to be interpreted as described in [RFC2119].

In addition, while most of the the terminology in this document is conform to the taxonomy presented in [[NC-Taxonomy]] for clarity and comparison with existing QUIC documents we continue to use the word packet to indication the entity that will be encoded vs. symbol in the taxonomy document. NOTE: while using drafts in references is not compliant with IETF/ IRTF rules they will be replaced by RFCs as they become available.

Introduction

The QUIC (Quick UDP-based Internet Connection)protocol is currently being proposed as new transport protocol than multiplexes connections over UDP. The major elements have been defined and are being implemented by the QUIC IETF working group [QUIC-WG] including wire format, connection establishment, stream multiplexing, stream and connection-level flow control, and data encryption [numerous draft references]. This document addresses an outstanding element of the QUIC protocol, namely how to account and correct for packet losses at that allow to vary with the size of the burst. That will be better suited to recover the losses with statistically approximately the same overhead as the average packet loss without limitations on the loss pattern.

Architecture

In order to define the potential NC solution, a detailed architecture is necessary. Hence, the main QUIC/NC architecture topics to be addressed includes the following (each will become a subsection in the future).

o Type of connection: while unicast and/or multicast/broadcast communications are possible over QUIC it is assumed that an initial implementation will be limited to unicast. o Addressing single source flow or multiple flows: at the the coding level, if there is a multiplexing on top of the coding level, already managed by QUIC machinery, it may be totally transparent. We can assume that individual packets and connections can be individually identified. o Use of feedback: since the code will need to deal with correlated losses can it benefit from feedback to manage the window as opposed to a fully unidirectional source-destination mechanism. This will allow not to lose any packet part of a current generation before the loss burst ends (we need a reference on window growth and maximum size) o Minimization of latency:Latency is key for the solution design.

This includes reducing extra delay due to encoding/decoding at the ingress, egress and intermediate nodes (middleboxes) especially on delay sensitive paths. At the same time long bandwidth-delay product networks coding should reduce the overall end-to-end delay experienced by an application significantly by minimizing the effect of packet losses and retransmission on TCP congestion control and throughput. o Throughput aspects: it is expected that the QUIC flows will include high throughput flows, very low throughput flows and mixed sizes flows. o Interactions with other functionalities: Interactons with congestion control and encryption will also be key. Directions will be taken from the beginning it should be possible to turn off the coding to prevent catastrophic congestion collapse for example.

Use-cases

Note: this will be detailed in the next version of the document

Requirements

The initial requirements for the QUIC NC are presented below. This list will help to chose the best solution amongst existing codes.

Requirements:

o Simplicity/low complexity: both encoding and decoding operations should be simple and ass little complexity to the QUIC operations; the use of systematic coding will be encouraged. o Low overhead: the NC overhead to compensate for all losses should be as close as possible to the average loss on the path as to not create additional congestion condition. o In network coding: there should be ways to create additional coded symbols inside the network either directly or via partial or full decoding. o Multipath: there should be ways to take advantage of multipath communications for example to send packets and coded symbols on different paths to reduce delay and overhead on some delay or loss sensitive paths. o Licensing/IPR: the solution should be license/patent free.

Next Steps

Besides adding the sections missing in the document based on future discussion it is proposed to define a strawman architecture based on existing codes and using the standard APIs being developed in the RG.

  9. IANA Considerations XX RFC ED -PLEASE REMOVE THIS SECTION XXX This memo includes no request to IANA. 10. Security Considerations Security: While NC will not impact security in itself it will be important to verify how NC interacts with current encryption used in QUIC and presented in [[QUIC-TLS]]. Swett, et al. Expires May 3, 2018 [Page 5]

  completed in a future version. For the needs of the current document we need to know that the QUIC packet format contains a unique packet identifier (ID) and a connection ID.

	Internet-Draft Internet-Draft	QUIC and NC QUIC and NC	October 2017 October 2017
	the network layer that will have a very negative impact on transport
	delay, throughput and reliability. This document presents the
	salient features and requirements for a network coding (NC) protocol
	to provide the QUIC packet loss recovery it requires. NC provides a
	structured, algebraic mechanism to recover lost packets based on a
	vast heritage of Forward Error Correction (FEC) and has shown better
	performance of packet recovery than XORs or repetition codes to deal
	with the losses in the Internet. The Network Coding taxonomy
	document [[NC-Taxonomy]] contains an overview of top NC concepts.
	Note: we need a small NC draft that explains how it works.
	3. QUIC Background
	This section will o Queuing losses across multiple flows
	o Intermittent timeouts
	o Connection losses
	o Residual physical or MAC layer losses
	o Misrouting
	o other?
	The main feature of all the patterns associated with the loss events
	above is the fact that losses appear in clusters (burst or correlated
	losses). Hence they are not the 'random losses' that can be
	recovered by non structured mechanisms like XOR or repetitions codes
	even with high overhead or simple block codes with fixed window
	sizes. Hence because of the correlated losses, the first requirement
	for a good code for QUIC is one that allows variable window sizes
	Swett, et al. Swett, et al.	Expires May 3, 2018 Expires May 3, 2018	[Page 2] [Page 3]

Details will be obtained from [[QUIC-Connect]] and [[QUIC-Trans]]

4. Motivation

The QUIC protocol from its early implementations, wanted to address packet losses in the Internet as they can greatly impact protocol performance and impact the performance congestion control mechanisms [[QUIC-Loss]]. For example TCP goodput goes below 20% with 3% loss [are there any other references besides the famous Mathis curve?]. In this section we review the motivations behind the use of a network coding approach to reduce the impact of packet losses on QUIC. It is important to note that we limit the sources of the losses to IP layer and above and losses in lower layers are addressed by other standardization organization.

It is known (is it?) that the main sources of losses in the Internet include (but are not limited to):

  [[QUIC-Loss]] and [[QUIC-TLS]] and other relevant documents o Code changes and future proofing: any protocol designed within QUIC should be able to maybe use more than one code to change codes easily by without major impact this is to address different network conditions or improve performance if a new code was to be developed. It is assumed that the code remain the same for the full QUIC session lifetime but that within a session at least for

	Internet-Draft	QUIC and NC	October 2017
	Swett, et al.	Expires May 3, 2018	[Page 4]

Internet-Drafts. The list of current Internet-Drafts is at https://datatracker.ietf.org/drafts/current/.

Internet-Drafts are draft documents valid for a maximum of six months and may be updated, replaced, or obsoleted


Nawel Medjkoune

Frédéric Armetta

Mathieu Lefort

Stefan Duffner

Autonomous object recognition in videos using Siamese Neural Networks

For a robot to be deployed in unconstrained real world environments, it needs to be autonomous. In this preliminary work, we focus on the capacity of an autonomous robot to discover and recognize objects in its visual field. Current existing solutions mainly employ complex deep neural architectures that need to be pre-trained using large datasets in order to be effective. We propose a new model for autonomous and unsupervised object learning in videos that does not require supervised pre-training and uses relatively simple visual filtering. The main idea relies on the saliency-based detection and learning of objects considered similar (thanks to a spatiotemporal continuity). For this purpose the learning of objects is based on a Siamese Neural Network (SNN). We demonstrate the capacity of the SNN to learn a good feature representation despite the deliberately simple and noisy process used to extract candidate objects.

I. INTRODUCTION

Nowadays, service and social robots are still limited in the quality and quantity of tasks they can accomplish. These limitations are mainly due to the poorly structured or unconstrained open environments they are operating in. In order to be efficient in these real-world settings, they would need to become more autonomous. Ideally, a robot progressively acquires informations from the outside world during its lifetime, similarly to a baby that is learning. This gain in terms of autonomy raises, amongst others, the problems of unsupervised, on-line and life-long learning in non-stationary environments [START_REF] Mazac | On bootstrapping sensori-motor patterns for a constructivist learning system in continuous environments[END_REF]. In this work, we focus on the capacity of an autonomous robot to discover and recognize objects localized in its visual field. Two main issues arise in this case: how to detect and localize new objects in this open environment and how to learn meaningful visual representations incrementally in an unsupervised way.

Multiple robust solutions exist for detecting, tracking and recognizing objects or persons in images and video streams, e.g. [START_REF] Leibe | Coupled object detection and tracking from static cameras and moving vehicles[END_REF], [START_REF] Ess | A mobile vision system for robust multi-person tracking[END_REF], [START_REF] Jafari | Real-time RGB-D based people detection and tracking for mobile robots and head-worn cameras[END_REF], [START_REF] Leal-Taixé | Learning by tracking: Siamese CNN for robust target association[END_REF], [START_REF] Borji | iLab-20M: A large-scale controlled object dataset to investigate deep learning[END_REF]. However, these methods are designed and pre-trained for specific object categories, like faces, persons or cars, and in a dedicated environment, and they are not able to accommodate new object categories during operation. To provide the visual perception algorithms with more autonomy, other works proposed unsupervised learning methods, based on neural networks, that are able to learn discriminative internal representations from unlabelled data and to adapt to a given context using transfer learning. For example, Wang et al. [START_REF] Wang | Unsupervised learning of visual representations using videos[END_REF] use a large unlabelled dataset of videos and specific low-level object patch extraction and tracking techniques to learn visual representations and then exploit this model to learn specific objects with labelled data. Another approach proposed by Liang et al. [START_REF] Liang | Towards computational baby learning: A weaklysupervised approach for object detection[END_REF] consists in continuously updating the learned visual representation using recognized objects from videos. However, their method incorporates a considerable amount of prior knowledge by using a Convolutional Neural Network (CNN), pre-trained on a large labelled dataset.

These existing methods are similar to our proposed approach in that they use unsupervised or weakly supervised learning techniques for visual representation learning exploiting the temporal coherence in videos. Although showing a high accuracy of object detection and recognition, they have some inherent limitations for an autonomous agent, e.g. they require a supervised pre-training on a given large image dataset in order to be effective. Moreover, the pre-training as well as their operation are computationally expensive due to the used complex Deep Neural Networks architectures. In our work, we address these issues by proposing a weakly supervised approach for object recognition using a convolutional SNN not requiring any labelled data or supervised pre-training. Our method uses no prior knowledge about the environment and only employs relatively simple visual filtering techniques inspired by human cognition (like saliency) to extract candidate objects.

II. MODEL

We propose a new model for autonomous and unsupervised object learning and recognition in videos using neural networks, illustrated in Fig. 1.

In our approach, we consider an autonomous agent equipped with a visual sensor (the camera) and an internal structure (the cognitive system). The agent is exposed to a visual flow coming from the environment and processes it frame by frame using low-level filters.

The main idea relies on the learning of objects considered similar (thanks to the saliency-based detection [START_REF] Itti | A model of saliencybased visual attention for rapid scene analysis[END_REF] and spatiotemporal continuity [START_REF] Hawkins | On intelligence: How a new understanding of the brain will lead to the creation of truly intelligent machines[END_REF]) as detailed below. For this purpose, we form pairs of images representing similar and dissimilar objects and learn this similarity metric using a convolutional SNN [START_REF] Bromley | Signature verification using a" siamese" time delay neural network[END_REF], [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF]. In a bootstrapping process, this learned similarity space can then be used to recognize previously seen objects. However, here, we are particularly interested in the first placement used for data. Indeed, without any knowledge of the objects and due to the relatively simple temporal filtering to detect candidate objects, the first extraction is coarse and highly sensitive to environment noise.

First, proto-object proposals are extracted. As we want the system to be functional without prior knowledge, we suggest a simple method for object proposals inspired by the human visual system. In a common visual scene, objects differ from their surroundings in color, intensity, texture, orientation, depth and other features. The human visual system is able detect regions that differ from their surroundings using bottomup saliency [START_REF] Treisman | A feature-integration theory of attention[END_REF], [START_REF] Nothdurft | Salience of feature contrast[END_REF]. Based on this observation, we propose to compute saliency as a cue to quickly form candidates of possible objects that we want to recognize, as illustrated in Fig. 2. A saliency map is computed from the original image and a thresholding is applied resulting in a binary map, where potential object candidates are highlighted with white blobs. Bounding boxes are then drawn around these blobs, representing object candidates.

Given that the input images are not labeled, we use the principle of spatio-temporal coherence of the objects present in the video sequences in order to provide a weakly supervised signal. The idea is that detections that are close in space and time are likely to correspond to the same object. Therefore, we associate detections in consecutive frames using spatial locality (see Fig. 3) and create "tubelets" of patches each one representing one object. Then, positive and negative pairs of patches (i.e. supposedly coming from the same and from different objects) are created based on these tubelets by associating respectively two patches from the same tubelet, and two patches from different tubelets. The learning is performed off-line by extracting all pairs from input videos.

The extracted pairs form the input of a SNN [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF], which transforms the data into a representation that brings closer similar data points and pulls apart dissimilar data points in the feature space by minimizing a contrastive loss function:

L(W ) = M i=1 L(W, Y i , X i 1 , X i 2 ) = (Y i ) 1 2 D 2 (X i 1 , X i 2 ) 2 + (1 -Y i ) 1 2 (max(0, m -D 2 (X i 1 , X i 2 ))) 2 , (1)
where (X i 1 , X i 2 ) are pair elements and Y i their label (positive or negative), M is the minibatch size, W the shared parameters, m an arbitrary margin and D 2 (X i 1 , X i 2 ) is the Euclidean distance between feature vectors X i 1 and X i 2 . In our experiments, we used a very light architecture for the SNN consisting of a 3 blocks of 10, 20 and 30 channel 3 × 3 convolution layer with input size 32 × 32 and ReLU activation function followed by a 2 × 2 max-pooling, and a fully-connected layer of size 20.

III. RESULTS

We evaluate the proposed approach on the RGB-D Scenes Dataset Version 2 [START_REF] Lai | A large-scale hierarchical multi-view rgb-d object dataset[END_REF]. The evaluation is done on both manually labeled classes of objects and automatically extracted training pairs which we augment randomly using data generation. The transformations made on the augmented pairs are: rotation, horizontal and vertical translation, shearing, zoom and horizontal flip. We train the SNN using several thousands of pairs and report the results after convergence.

First, we use the representations learned by the SNN in order to get the feature representation of images in four manually labeled classes. A TSNE dimensionality reduction is then applied to project the transformed data points, illustrated in Fig. 4. In this example, training is done on data pairs extracted from video 12 (a representative example) and the projection on manually labeled data from video 12 and 14 respectively (both videos present the same objects in different backgrounds). The plots show that there is a certain structure when projecting manually labeled data images, with four separated clusters with a small confusion when generalizing. Secondly, we compare the classification accuracy of the manually labeled dataset using two different feature representations: one learned from a supervised CNN trained on the same dataset, and the other learned using our model. The two models have similar architectures and are trained with the same parameters. Table I shows the results of three experiments where we chose 2, 4 or 6 different objects as follow: first we take two very distinct classes (cup and hat) and then add gradually different classes of objects that can look alike (e.g. cups with different colors, cup and bowl with the same color). The SNN gives comparable results to the supervised CNN with few classes, and falls to 74.15% accuracy in the presence of similar objects by color or shape.

These experiments show two results: when performing weakly supervised training of an SNN, with a relatively simple architecture and no prior pre-training, good visual representation can be learned despite the noise in the unlabeled input data, with a fair generalization to objects present in different backgrounds. Moreover, these results also show that this weakly supervised approach, only based on the feedback from spatio-temporal coherence, can be used to bootstrap the learning process.

This work is a preliminary step in order to design an autonomous system for object recognition using neural networks. We showed with these first experiments the capacity of a simple neural network in learning visual feature representations of objects in simple visual scenes. In the next step, we aim to validate these results in more complex realworld scenes before directing it towards on-line incremental learning with a feedback loop (Fig. 1, arrow 4) to guide the localization of objects using the trained neural network.

IV. ACKNOWLEDGEMENT

Fig. 1 .

 1 Fig.1. Overview of the proposed approach. The agent captures the visual field frame by frame and processes it to discover proto-objects (arrow 1). Using spatio-temporal coherence, pairs of similar and dissimilar objects are created and given to the neural network (arrow 2) that learns an internal visual representation of the objects (arrow 3).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Steps of object discovery using saliency attention model. a: original image, b: saliency map, c: binary map after thresholding the saliency map, d: bouding boxes around proto-objects

Fig. 4 .

 4 Fig. 4. TSNE dimensionality reduction of feature representations learned by the Siamese Neural Network. Each color represents one object: Cup, Bowl, hat and soda can. a: manually labaled dataset from the video 12. b: manually labeled dataset from video 14 (same objects, different background).

We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research.


Lazhar Dhaouadi
email: lazhardhaouadi@yahoo.fr

Hardy's theorem for the q-Bessel Fourier transform

In this paper we give a q-analogue of the Hardy's theorem for the q-Bessel Fourier transform. The celebrated theorem asserts that if a function f and its Fourier transform

1 The q-Foureir Bessel transform Throughout this paper we consider 0 < q < 1 and we adopt the standard conventional notations of [START_REF] Gasper | Basic hypergeometric series, Encycopedia of mathematics and its applications 35[END_REF]. We put R + q = {q n , n ∈ Z}, and for complex a (a; q) 0 = 1, (a; q) n = n-1 i=0

(1 -aq i ), n = 1...∞.

Jackson's q-integral (see [START_REF] Jackson | On a q-Definite Integrals[END_REF]) in the interval [0, ∞[ is defined by

∞ 0 f (x)d q x = (1 -q) ∞ n=-∞ q n f (q n ).
We introduce the following functional spaces L q,1,ν of even functions f defined on R + q such that

f q,1,v = ∞ 0 |f (x)|x 2v+1 d q x < ∞.
The normalized Hahn-Exton q-Bessel function of order ν > -1 (see [START_REF] Swarttouw | The Hahn-Exton q-Bessel functions[END_REF])is defined by

j ν (z, q) = ∞ n=0 (-1) n q n(n-1) 2
(q, q) n (q ν+1 , q) n z n .

It's an entire analytic function in z.

Lemma 1 For every p ∈ N, there exist σ p > 0 for which

|z| 2p |j ν (z, q 2 )| < σ p e |z| , ∀z ∈ C.
Proof. In fact

|z| 2p |j ν (z, q 2 )| ≤ 1 (q 2 , q 2 ) ∞ (q 2ν+2 , q 2 ) ∞ ∞ n=0 q n(n-1) |z| 2n+2p ≤ q p(p+1) (q 2 , q 2 ) ∞ (q 2ν+2 , q 2 ) ∞ ∞ n=p q n(n-2p-1) |z| 2n .

Now using the Stirling's formula

n! ∼ √ 2πn n n e n ,
we see that there exist an entire n 0 ≥ p such that

q n(n-2p-1) < 1 (2n)! , ∀n ≥ n 0 , which implies ∞ n=n 0 q n(n-2p-1) |z| 2n < ∞ n=n 0 1 (2n)! |z| 2n < e |z| .
Finally there exist σ p > 0 such that

|z| 2p |j ν (z, q 2 )| e |z| < σ p , ∀z ∈ C.
This finish the proof.

The q-Bessel Fourier transform F q,ν introduced in [START_REF] Dhaouadi | Inequalities in q-Fourier Analysis[END_REF][START_REF] Koornwinder | On q-Analogues of the Hankel and Fourier Transform[END_REF] as follow

F q,ν f (x) = c q,ν ∞ 0 f (t)j ν (xt, q 2 )t 2ν+1 d q t.
where

c q,ν = 1 1 -q (q 2ν+2 , q 2 ) ∞ (q 2 , q 2 ) ∞ .
The following theorem was proved in [START_REF] Dhaouadi | Inequalities in q-Fourier Analysis[END_REF] Theorem 1 Given f ∈ L q,1,ν then we have

F 2 q,ν (f )(x) = f (x), ∀x ∈ R + q .
Proof. See [1] p 3.

Hardy's theorem

The following Lemma from complex analysis is crucial for the proof of our main theorem.

Lemma 2 Let h be an entire function on C such that

|h(z)| ≤ Ce a|z| 2 , z ∈ C, |h(x)| ≤ Ce -ax 2 , x ∈ R,
for some positive constants a and C. Then h(z) = Const.e -ax 2 .

Proof. See [START_REF] Sitaram | An analogue of Hardy's theorem for very rapidly decreasing functions on semi-simple Lie groups[END_REF] p 4. Now we are in a position to state and prove the q-analogue of the Hardy's theorem Theorem 2 Suppose f ∈ L q,1,ν satisfying the following estimates

|f (x)| ≤ Ce -1 2 x 2 , ∀x ∈ R + q , |F q,ν f (x)| ≤ Ce -1 2 x 2 , ∀x ∈ R, where C is a positive constant. Then there exist A ∈ R such that f (z) = A.c q,ν F q,ν e -1 2 x 2 (z), ∀z ∈ C.
Proof. We claim that F q,ν f is an analytic function and there exist C ′ > 0 such that

|F q,ν f (z)| ≤ C ′ e 1 2 |z| 2 , ∀z ∈ C. We have |F q,ν f (z)| ≤ c q,ν ∞ 0 |f (x)||j ν (zx, q 2 )|x 2ν+1 d q x.
From the Lemma 1, if |z| > 1 then there exist σ 1 > 0 such that

x 2ν+1 |j ν (zx, q 2 )| = 1 |z| 2ν+1 (|z| x) 2ν+1 |j ν (zx, q 2 )| < σ 1 1 + |z| 2 x 2 e x|z| , ∀x ∈ R + q .
Then we obtain

|F q,ν f (z)| ≤ Cσ 1 c q,ν ∞ 0 e -1 2 (x-|z|) 2 1 + |z| 2 x 2 d q x e 1 2 |z| 2 < Cσ 1 c q,ν ∞ 0 1 1 + x 2 d q x e 1 2 |z| 2 .
Now, if |z| ≤ 1 then there exist σ 2 > 0 such that

x 2ν+1 |j ν (zx, q 2 )| ≤ σ 2 e x , ∀x ∈ R + q . Therefore |F q,ν f (z)| ≤ Cσ 2 c q,ν ∞ 0 e -1 2 x 2 +x d q x ≤ Cσ 2 c q,ν ∞ 0 e -1 2 x 2 +x d q x e 1 2 |z| 2 .
Which leads to the estimate (1). Using Lemma 2, we obtain

F q,ν f (z) = const.e -1 2 z 2 , ∀z ∈ C,
and by theorem 1, we conclude that

f (z) = const.F q,ν e -1 2 t 2 (z), ∀z ∈ C.
This finish the proof.

Corollary 1 Suppose f ∈ L q,1,ν satisfying the following estimates

|f (x)| ≤ Ce -px 2 , ∀x ∈ R + q , |F q,ν f (x)| ≤ Ce -σx 2 , ∀x ∈ R,
where C, p, σ are a positive constant and pσ = 1 4 . We suppose that there exist a ∈ R + q such that a 2 p = 1 2 . Then there exist A ∈ R such that

f (z) = A.c q,ν F q,ν e -σt 2 (z), ∀z ∈ C.
Proof. Let a ∈ R + q , and put

f a (x) = f (ax), then F q,ν f a (x) = 1 a 2ν+2 F q,ν f (x/a).
In the end, applying Theorem 2 to the function f a .

Remark 1 Using the q-Central limit Theorem (see [START_REF] Dhaouadi | Inequalities in q-Fourier Analysis[END_REF]) we give a probability interpretation of the function c q,ν F q,ν e -σt 2 . In fact if (ξ n ) n≥0 be a sequence of positive probability measures of R + q , satisfying lim n→∞ nσ n = (q 2 , q 2 ) 1 (q 2v+2 , q 2 ) 1 q 2 σ, where

σ n = ∞ 0 t 2 t 2v+1 d q ξ n (t),
and

lim n→∞ n σ n = 0, where σ n = ∞ 0 t 4 1 + t 2 t 2v+1 d q ξ n (t),
then the nth q-convolution product ξ * n n converge strongly toward a measure ξ defined by d q ξ(x) = c q,ν F q,ν e -σt 2 (x)d q x.

Corollary 2 Suppose f ∈ L q,1,ν satisfying the following estimates

|f (x)| ≤ Ce -px 2 , ∀x ∈ R + q , |F q,ν f (x)| ≤ Ce -σx 2 , ∀x ∈ R,
where C, p, σ are a positive constant and pσ > 1 4 .We suppose that there exist a ∈ R + q such that a 2 p = 1 2 . Then f ≡ 0.

Proof. In fact there exist σ ′ < σ such that pσ ′ = 1 4 . Then the function f satisfying the estimates of Corollary 1, if we replacing σ by σ ′ . Which implies

F q,ν f (x) = const.e -σ ′ x 2 , ∀x ∈ R.
On the other hand, f satisfying the estimates of Corollary 2, then |const.e -σ ′ x 2 | ≤ Ce -σx 2 , ∀x ∈ R. This implies F q,ν f ≡ 0, and by Theorem 1 we conclude that f ≡ 0 Remark 2 Hardy's theorem asserts the impossibility of a function and its q-Fourier Bessel transform to be simultaneously "very rapidly decreasing". Hardy's theorem can also be viewed as a sort of "Qualitative uncertainty principles". One such example can be the fact that a function and its q-Bessel Fourier transform cannot both have compact support.