Virginie Ehrlacher

Tony Lelièvre

Pierre Monmarché

Adaptive force biasing algorithms: new convergence results and tensor approximations of the bias

Keywords: Monte Carlo methods, tensor, free energy, importance sampling, molecular dynamics. MSC class (2010): 65C05, 65N12

We analyze and propose variants of the Adaptive Biasing Force method. First, we prove the convergence of a version of the algorithm where the biasing force is estimated using a weighted occupation measure, with an explicit asymptotic variance. Second, we propose a new flavour of the algorithm adapted to high dimensional reaction coordinates, for which the standard approaches suffer from the curse of dimensionality. More precisely, the free energy is approximated by a sum of tensor products of one-dimensional functions. The consistency of the tensor approximation is established. Numerical experiments on 5dimensional reaction coordinates demonstrate that the method is indeed able to capture correlations between them.

Introduction

Consider x ∈ T D a vector representing the positions of particles with periodic boundary conditions (T = R/Z), and a potential energy V ∈ C ∞ T D . Our first motivation is the computation, in multimodal cases, of expectations of the form 1 T D e -βV (x) dx T D ϕ(x)e -βV (x) dx =:

T D ϕdµ V,β (1)
where ϕ : T D → R is called an observable and dµ V,β := e -βV (x) dx is the Gibbs law with potential V and inverse temperature β > 0. Our second motivation is the estimation of the free energy associated to this Gibbs law and a given high-dimensional reaction coordinate, as will be explained below. Our main contributions are, first, the long-time convergence of an adaptive Importance Sampling method for estimating (1) (Theorem 2) with an explicit asymptotic variance (Theorem 3) and, second, the definition of a variant of this algorithm based on a tensor approximation (Algorithms 2 and 3) to tackle high-dimensional reaction coordinates and the proof of its consistency (Theorem 5).

The large dimension D is so significant that, in practice, the quantities [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF] have to be computed with Markov Chain Monte Carlo (MCMC) algorithms, which consist in approximating the average of ϕ with respect to µ V,β along dynamics that are ergodic with respect to µ V,β . A typical sampler is the overdamped Langevin dynamics dX t = -∇V (X t )dt + 2β -1 dB t where (B t ) t 0 is a Brownian motion over T D . It is ergodic with invariant measure µ V,β , so that

1 t t 0 ϕ(X s )ds -→ t→∞ T D ϕdµ V,β a.s.
for all measurable bounded ϕ, see e.g. [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF] and references therein. Nevertheless, the convergence of the process (or, in practice, of any alternative Markov process with invariant measure µ V,β ) toward its equilibrium in the long-time limit may be very slow. This is due to the so-called metastability phenomenon, according to which the process remains for long times in some region of the space, with very rare transitions from one of these metastable regions to another. This is related to the multi-modality of the Gibbs measure and the fact MCMC algorithms typically perform local moves, so that leaving a mode of the target measure µ V,β is a rare event. We refer to [START_REF] Lelièvre | Two mathematical tools to analyze metastable stochastic processes[END_REF] for more details on this topic. For this reason, several adaptive methods have been developed in order to force the process to leave the metastable traps faster.

Among those, we focus on the adaptive biasing force (ABF) algorithm, which may be seen as a particular Importance Sampling method. The general idea is to run a biased process d Xt = -∇V ( Xt )dt + ∇V bias,t ( Xt )dt + 2β -1 dB t

where the biasing potential V bias,t is adaptively constructed from the past trajectory ( Xs ) s∈[0,t] in such a way that it is expected to converge to some V bias,∞ . Expectations with respect to µ V,β are then recovered through a reweighting step, assuming that ergodicity still holds:

1 t t 0 ϕ( Xs )e -βV bias,s ( Xs) ds 1 t t 0 e -βV bias,s ( Xs) ds

-→ t→∞ T D ϕe -βV bias,∞ dµ V -V bias,∞ ,β T D e -βV bias,∞ dµ V -V bias,∞ ,β = T D ϕdµ V,β . (3)
Classically, in such an Importance Sampling scheme, the aim is to design a target bias V bias,∞ such that two conditions are met: 1) sampling the biased equilibrium µ V -V bias,∞ ,β is simpler than the initial problem (i.e. the corresponding overdamped Langevin process is less metastable) and 2) the biased equilibrium is not too far from the initial target so that the exponential weights in (3) do not cause the asymptotical variance of the estimator to skyrocket.

In the ABF algorithm, this issue is addressed with the use of so-called reaction coordinates (or collective variables) and the associated free energy as a bias. Reaction coordinates consist of a small number d D of macroscopic coordinates of the whole microscopic system x ∈ T D . These coordinates are defined through a map ξ : T D → M where M is a manifold of dimension d. In molecular dynamics, for example, x ∈ T D is a vector which gathers the positions of all the different atoms of the system of interest, and ξ(x) typically represents some distances between particular pairs of atoms, or angles formed by some triplets of atoms. These reaction coordinates should be chosen to capture the main causes of the metastability of the system. More precisely, ξ(X t ) should converge to equilibrium as slowly as X t , while the conditional laws L(X | ξ(X) = z) for fixed z ∈ M when X ∼ µ V,β should be easier to sample (see [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF] or Section 2.4 for more detailed considerations). In other words, ξ(x) should be a low-dimensional representation of x that captures the slow variables of the system.

To these reaction coordinates ξ is associated the corresponding free energy A : M → R, given by

A(z) = - 1 β ln {x∈T D , ξ(x)=z}
e -βV (x) δ ξ(x)-z (dx) ,

where δ ξ(x)-z is the so-called delta measure, which can be defined from the Lebesgue measure on the submanifold {x ∈ T D , ξ(x) = z} through the co-area formula, see for example [START_REF] Lelièvre | Free energy computations: A mathematical perspective[END_REF]Section 3.2.1]. This definition ensures that, if X is a random variable with law µ V,β on T D , then ξ(X) is a random variable with law µ A,β on M. The heuristic of the ABF algorithm is the following. Suppose that M is compact. If we were to sample from the process

dY t = -∇ (V -A • ξ) (Y t )dt + 2β -1 dB t , (4)
the equilibrium would be µ V -A•ξ,β , whose image through ξ, by definition of A, is the uniform measure on M. This means that there would be no more metastability along ξ, since all the regions of M would be equally visited by ξ(Y t ). Unfortunately, it is not possible to use directly this free-energy biased dynamics in practice, since it would require the knowledge of A and thus the computation of expectations in large dimension. The idea of the ABF method is to learn A on the fly, i.e. to run a process ( Xt ) t 0 solving (2) with a biasing potential V bias,t constructed from ( Xs ) s∈[0,t] and designed to target A • ξ in the longtime limit.

As mentioned at the beginning of this introduction, in fact, independently from the question of obtaining MCMC estimations for [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF], estimating the free energy A associated to some reaction coordinates ξ can be an objective per se. Indeed, in molecular dynamics (for either chemistry, biophysics or material science), it is a key quantity to analyze phenomena such as conformational stability, solvation properties or ligand binding affinity, see e.g. [START_REF] Guo | Efficient and flexible algorithm for free energy calculations using the λ-dynamics approach[END_REF][START_REF] Lelièvre | Free energy computations: A mathematical perspective[END_REF][START_REF] Jorgensen | Monte carlo simulation of differences in free energies of hydration[END_REF][START_REF] Bash | Free energy calculations by computer simulation[END_REF] and references within. In these cases, the reaction coordinates work as a small dimensional representation of the macroscopic state of the system. This is also the case in coarse-graining problems, where the goal is to replace numerically intensive high-dimensional molecular simulations of a microscopic system X t by cheaper simulations of an approximate effective dynamics for the low-dimensional ξ(X t ). In that case, the effective force is usually given by the free energy, which has to be estimated on an all-atom high dimensional simulation, see e.g. [START_REF] Lu | The multiscale coarse-graining method. vii. free energy decomposition of coarse-grained effective potentials[END_REF][START_REF] Monticelli | The martini coarse-grained force field: Extension to proteins[END_REF] and references within.

In practice, the choice of good reaction coordinates, both in term of enhanced sampling and of low-dimensional representation of the system, is a difficult problem. Up to recently, their definition has been based on the knowledge and intuition of specialists. The question of the automatic learning of suitable reaction coordinates is currently a vivid research area, see for instance [START_REF] Bittracher | Data-driven computation of molecular reaction coordinates[END_REF][START_REF] Brandt | Machine learning of biomolecular reaction coordinates[END_REF] and the recent review [START_REF] Gkeka | Machine learning force fields and coarsegrained variables in molecular dynamics: application to materials and biological systems[END_REF]. Moreover, some techniques like the orthogonal space random walk [START_REF] Lu | Predictive sampling of rare conformational events in aqueous solution: Designing a generalized orthogonal space tempering method[END_REF] provide a general way to construct new reaction coordinates from previous ones. Due to these recent progresses, one would like to consider a relatively large d. In ABF, V bias,t is a function of the d reaction coordinates. From a numerical point of view, since V bias,t is adaptively learned on the fly, its values have to be kept in memory, which requires a grid whose size typically scales exponentially with d. This limits the application of ABF to small dimensional reaction coordinates (d 4). The aim of the present work is to lift this limitation by approximating V bias,t using a sum of tensor products of one-dimensional functions, which reduces the size of the memory to O(dm) where m is the number of tensor terms. Remark that this can in turn help for the definition of good reaction coordinates, by considering as candidates a relatively large number of reaction coordinates and then conduct a statistical study to select or combine some of them. A basic idea would be to conduct a sensitivity analysis of the free energy, computing for instance for each reaction coordinate ξ i the best approximation in the least square sense of the (estimated) free energy by a function only of the other reaction coordinates (ξ j ) j =i (which is easily done for a function given as a sum of tensor products, see [START_REF] Konakli | Global sensitivity analysis using low-rank tensor approximations[END_REF]) and then discarding the reaction coordinate whose disparition gives the lowest error. Nevertheless, this question exceeds the scope of the present work, in which ξ is supposed to be given.

Note that the question of increasing the number of reaction coordinates in adaptive biasing algorithms has also been considered in the Bias-Exchange algorithm introduced in [START_REF] Piana | A bias-exchange approach to protein folding[END_REF], where several replicas of the system are run in parallel, each associated with a one-dimensional reaction coordinate. The replicas exchange their bias according to some Metropolis-Hastings probability, so that each replica eventually feels the bias in all the different directions of the reaction coordinates. Nevertheless, in this case where one-dimensional reaction coordinates are treated independently one from the others, correlations between reaction coordinates cannot be measured (the same goes for the generalized ABF introduced in [START_REF] Zhao | The extended generalized adaptive biasing force algorithm for multidimensional free-energy calculations[END_REF]), contrary to the algorithm introduced in the present work.

Besides, let us mention that numerical methods involving both tensor approximation and Monte Carlo methods for molecular dynamics are also introduced in [START_REF] Klus | Towards tensor-based methods for the numerical approximation of the Perron-Frobenius and Koopman operator[END_REF][START_REF] Nüske | Variational tensor approach for approximating the rare-event kinetics of macromolecular systems[END_REF] for other purposes.

The rest of the paper is organized as follows. The algorithms are introduced in Section 2, where our main results are stated and discussed. Section 3 is devoted to the proofs of Theorems 2 and 3 which are related to the long-time convergence of the ABF algorithm. In Section 4, we prove Theorem 5, which concerns the tensor approximation scheme. Section 5 provides a detailed discussion on practical considerations and possible variations of the algorithm. Finally, some numerical experiments with the tensor ABF algorithm are reported in Section 6.

Algorithms and results

In this section we provide a presentation of the ABF algorithm which is considered in this work in a simple framework, and refer to Section 5 for generalizations. The presentation is divided into two parts. In Section 2.1, we present the reference ABF algorithm that we consider, without the tensor-product approximation, and state a long-time convergence result and a CLT. In Section 2.2, we introduce the tensor-product approximation of the bias and state its consistency. These two ingredients are then combined to yield the Tensor-ABF algorithm in Section 2.3. The two algorithms and associated convergence proofs of the reference ABF algorithm and of the tensor-product approximation are presented separately since they have their own interest.

Free energy and the ABF algorithm

Let us first present the ABF algorithm in a simple framework (see [START_REF] Hénin | Overcoming free energy barriers using unconstrained molecular dynamics simulations[END_REF][START_REF] Darve | Calculating free energies using average force[END_REF][START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF] for more general settings). From now on, we write µ = µ V,β , seen both as a probability law and as the density of the latter with respect to the Lebesgue measure.

Let us assume that M = T d and that, for all x = (q, z)

∈ T D = T p × T d , ξ(x) = ξ(q, z) = z where p = D -d.
At first sight, this may seem a very restrictive choice of reaction coordinates. But, using extended variables (see [START_REF] Fu | Extended adaptive biasing force algorithm. an onthe-fly implementation for accurate free-energy calculations[END_REF]), this can actually be applied in very general contexts. We refer the reader to Section 5 for more details on this point.

The associated free energy for z ∈ T d is then

A(z) = - 1 β ln T p
e -βV (q,z) dq .

Following the previous discussion, our aim is then to define for all time t 0 a function A t on T d and to sample the process

dQ t = -∇ q V (Q t , Z t )dt + 2β -1 dB 1 t dZ t = -∇ z V (Q t , Z t )dt + ∇ z A t (Z t )dt + 2β -1 dB 2 t , (5)
where B 1 and B 2 are independent Brownian motions respectively of dimension p and d, in such a way that A t gets close to A in large time.

Note that the free energy A satisfies

∇ z A(z) = T p ∇ z V (q, z)e -βV (q,z) dq T p e -βV (q,z) dq = E µ [∇ z V (Q, Z) | Z = z] .
The following alternative equivalent characterization of A will be useful in the sequel. Denoting by H 1 (T) the set of functions of L 2 (T d ) with a weak gradient in L 2 (T d ), define

H := f ∈ H 1 (T d ) : T d f (z)dz = 0 , (6)
and let us denote by P(T p × T d ) the set of probability measures on T p × T d = T D . For all ν ∈ P(T p × T d ) and f ∈ H, let us define

E ν (f ) := T p ×T d |∇ z V (q, z) -∇ z f (z)| 2 dν(q, z).
As detailed in [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF], up to an additive constant (like the potential V , the free energy is in fact always defined up to an additive constant), A is the unique minimizer in H of the functional E µ , i.e.

A -

T d A(z)dz = argmin f ∈H E µ (f ). ( 7
)
At time t 0, a trajectory (Q s , Z s ) s∈[0,t] of ( 5) is available. Let ν t be the probability measure on T p × T d = T D defined as follows: for all ϕ ∈ C(T p × T d ),

T p ×T d ϕdν t = t 0 e -βAs(Zs) ds -1 t 0 ϕ (Q s , Z s ) e -βAs(Zs) ds . (8)
We call ν t the unbiased occupation distribution of the process. By the ergodic limit (3), ν t is expected to converge weakly to µ as t goes to infinity almost surely (at least if A t does not change too fast with t). However, note that ν t is a singular probability measure, so that the minimization problem

inf f ∈H E νt (f )
is ill-posed. To circumvent this difficulty, one may consider two different alternatives to regularize the problem which we detail hereafter. Consider a smooth symmetric positive density kernel

K ∈ C ∞ (T d × T d , R + ) with T d K(y, z)dz = 1 and K(y, z) = K(z, y) ∀y, z ∈ T d . (9)
In practice, K(y, •) should be close to a Dirac mass at y (see Theorem 2 below). For instance, a possible choice for K would be the so-called von-Mises kernel for a given small parameter ε > 0, i.e.

K(y, z)

∝ d i=1 exp - 1 ε 2 /2 sin 2 z i -y i 2 . (10)
Now, consider also a regularization parameter λ 0. For all ν ∈ P(T p × T d ) and all f ∈ H, we define

J ν (f ) := T p ×T d ×T d |∇ y V (q, y) -∇ z f (z)| 2 K(y, z)dzdν(q, y) + λ T d |∇ z f (z)| 2 dz , (11)
Note that, as K(y, •) converges weakly toward the Dirac mass at y and λ goes to 0, for all f ∈ H, J ν (f ) converges towards E ν (f ). The interest of introducing J ν is that, thanks to the regularization, the minimization problem is now well-posed:

Proposition 1. Assume that either K > 0 on T d ×T d or λ > 0.
Then, for all ν ∈ P(T p ×T d ), J ν admits a unique minimizer in H.

This is a direct consequence of the strict convexity of J , see Section 4. In summary, in the whole article, we work under the following conditions.

Assumption 1. V ∈ C ∞ (T D ), D 3, β > 0, λ 0 and K ∈ C ∞ (T d × T d , R + ) satisfies (9). Moreover, either K > 0 or λ > 0.
We now have all the elements to define the reference ABF algorithm in this work, see Algorithm 1 below.

Algorithm 1 ABF algorithm

1: Input: 2: Initial condition (q 0 , z 0 ) ∈ T p × T d 3: Brownian motion (B 1 t , B 2 t ) t 0 on T p × T d 4: Regularization parameters K, λ 5: Update period T up > 0, number of updates N up ∈ N * , total simulation time T tot = T up N up 6: Output:

7: Estimated free energy A Ttot ∈ H 8: Trajectory (Q t , Z t ) t∈[0,Ttot] ∈ C [0, T tot ], T p × T d 9: Begin: 10: Set (Q 0 , Z 0 ) = (q 0 , z 0 ). 11: Set A 0 (z) = 0 for all z ∈ T d .
12: Set t k = kT up for all k ∈ 0, N up . 13: for k ∈ 1, N up do 14:

Set A t = A t k-1 for all t ∈ [t k-1 , t k ).
Set (Q t , Z t ) t∈[t k-1 ,t k ]
to be the solution of (5) with initial condition (Q t k-1 , Z t k-1 ) at time t k-1 .

16:

Set A t k to be the minimizer in H of J νt k given by ( 8) and [START_REF] Brandt | Machine learning of biomolecular reaction coordinates[END_REF]. 17: end for

18: Return A Ttot and (Q t , Z t ) t∈[0,Ttot] .
Remark that, contrary to the cases studied in other theoretical works like [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF][START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF][START_REF] Benaïm | Convergence of adaptive biasing potential methods for diffusions[END_REF], in Algorithm 1, the bias A t is piecewise constant in time, with updates at the times t k , k ∈ 1, N up . This is due to the fact that, as will be detailed in Section 2.2, in addition to the classical case, we are also interested in a case where a bias update is numerically demanding, and thus cannot be performed at each time step.

We prove in Section 3 the long-time convergence of Algorithm 1:

Theorem 2. Under Assumption 1, let (Q t , Z t , A t ) t 0 be given by Algorithm 1 (with a fixed T up > 0 and N up = +∞ so that T tot = +∞ and the process is defined for all positive times). Then, as t → +∞, almost surely, ν t given by (8) weakly converges toward µ and

∇A t -∇A * ∞ -→ t→∞ 0,
where A * is the unique minimizer in H of J µ . Moreover, A * satisfies

T d |∇A(z) -∇A * (z)| 2 T p ×T d K(y, z)µ(q, y)dqdy dz 4 ∇ 2 A ∞ sup y∈T d T d |y -z| 2 K(y, z)dz + 2λ T d |∇A(z)| 2 dz . ( 12
)
Note that [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF] implies that, as λ and sup y∈T d T d |y -z| 2 K(y, z)dz go to zero, A * converges in H to A -A (which corresponds to λ = 0 and K(y, z) = δ y (z)).

The almost sure weak convergence of ν t toward µ implies, of course, the almost sure convergence of the importance sampling estimator T D ϕdν t toward the target T D ϕdµ for all continuous observable ϕ.

The long-time convergence of a similar ABF algorithm has been established in [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF] but in a case where, instead of its occupation measure, the process interacts with its law at time t. Rather than a self-interacting process (i.e. a single trajectory with memory), this corresponds to a system of N interacting particles (with no memory), and more precisely to the mean-field limit as N goes to infinity of this system. The techniques to study such a non-linear process is completely different from our non-Markovian case. Moreover, a result similar to Theorem 2 has been established in [START_REF] Benaïm | Convergence of adaptive biasing potential methods for diffusions[END_REF] for a closely related self-interacting process, the adaptive biasing potential algorithm. In addition, in the recent work [START_REF] Benaïm | Analysis of an Adaptive Biasing Force method based on self-interacting dynamics[END_REF], a similar result is established for the ABF algorithm but when the occupation measure is not unbiased (see the discussion in Section 5.3).

The previous qualitative result states that the algorithm is consistent, but gives no information on its efficiency. We now state that the asymptotic variance of the estimators obtained from the ABF algorithm is the same as in the case of a process with constant biasing potential equal to A * . More precisely, consider

X * = (Q * , Z * ) the solution of dQ * t = -∇ q V (Q * t , Z * t )dt + 2β -1 dB 1 t dZ * t = -∇ z V (Q * t , Z * t )dt + ∇ z A * (Z * t )dt + 2β -1 dB 2 t
, where A * is given by Theorem 2, and let

ν * t = t 0 e -βA * (Zs) ds -1 t 0 δ X * s e -βA * (Zs) ds .
Theorem 3. Under the settings of Theorem 2, there exist C > 0 such that for all t 0 and all ϕ ∈ C T D ,

E |ν t (ϕ) -µ(ϕ)| 2 C t ϕ 2 ∞ .
Moreover, tE (|ν t (ϕ) -µ(ϕ)| 2 ) converges as t → +∞ to a limit σ 2 ∞ (ϕ) ∈ R + , which is also the limit of tE (|ν * t (ϕ) -µ(ϕ)| 2 ). As detailed in the proof of Theorem 3 (in Section 3), the asymptotic variance is given as follows:

σ 2 ∞ (ϕ) := 2 β T D e βA * •ξ |∇ψ| 2 dµ T D e βA * •ξ dµ
where ψ solves the Poisson equation

1 β ∆ -∇ (V -A * • ξ) ∇ ψ = e -βA * •ξ ϕ - T D ϕdµ .
The consequences of Theorem 3 in term of efficiency of the algorithm are discussed in Section 2.4.

Tensor approximation

This section focuses on the minimization step of Algorithm 1. Assumption 1 is enforced. Fix ν ∈ P(T p × T d ). For all f ∈ H, the cost function J ν (f ) defined by ( 11) is equal to

J ν (f ) = C ν + (1 + λ) T d |F ν (z) -∇f (z)| 2 θ ν (z)dz ( 13
)
with some constant C ν independent from f and where, for all z ∈ T d ,

θ ν (z) := 1 λ + 1 λ + T p ×T d K(y, z)dν(q, y) (14)
F ν (z) := 1 (λ + 1)θ ν (z) T p ×T d ∇ y V (q, y)K(y, z)dν(q, y) .

Note that, under Assumption 1, θ ν is the density of a probability measure, bounded from below by (λ + min K)/(1 + λ) > 0. Moreover, since K is smooth and bounded, so are θ ν and F ν . Note that neither the additive constant C ν nor the multiplication by 1 + λ affect the problem of minimizing J ν . As a consequence, the unique minimizer f * of J ν on H (see Proposition 1) is equivalently the unique minimizer of

H f → J ν (f ) := T d |F ν (z) -∇f (z)| 2 θ ν (z)dz. ( 16
)
The gradient of J ν at f * is the Helmholtz projection in L 2 (θ ν ) of F ν . The Euler-Lagrange equation associated to the minimization problem of J ν over H is

∇ • (θ ν (∇f * -F ν )) = 0 , (17)
where ∇• denotes the divergence operator. When d is small (d = 2 in [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF]), as t increases, the functions θ νt and θ νt F νt are updated and kept in memory on a discrete grid of dimension M d for some M ∈ N * , and the Euler equation is solved with standard PDE techniques. However this is not sustainable if one wants to consider a larger number of reaction coordinates. For this reason, we now present a method to approximate f * by a sum of tensor products, namely by a function f m ∈ H which reads as follows

∀z := (z 1 , • • • , z d ) ∈ T d , f m (z) = m k=1 d j=1 r k,j (z j )
for some m ∈ N * and some functions r k,j : T → R for 1 ≤ j ≤ d and 1 ≤ k ≤ m. See [START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF] for a general overview on tensor methods. Let g be a simple tensor product function, i.e. a function such that for all z = (z

1 , • • • , z d ) ∈ T d , g(z) = d j=1 r j (z j ) for some r 1 , • • • , r d ∈ H 1 (T)
. Such a simple tensor product function will be denoted herefater by g = d j=1 r j . If g belongs to H, its (Lebesgue) integral vanishes, which is equivalent to the fact there exists i ∈ 1, d such that the (Lebesgue) integral of r i vanishes. This motivates the introduction of the following subspaces of H: for i ∈ 1, d , define

Σ i := g ∈ H, g = d j=1
r j with r j ∈ H 1 (T) for all j ∈ 1, d and

T r i (z i )dz i = 0 .
Proposition 4. Under Assumption 1, for all ν ∈ P(T p × T d ), i ∈ 1, d and f ∈ H, there always exists at least one minimizer in Σ i to the optimization problem

min g∈Σ i J ν (f + g). (18)
This is proven in Section 4. From Proposition 4, the greedy algorithm described in Algorithm 2 below is well-defined.

Algorithm 2 Greedy(ν, f 0 , m) 1: Input: 2: Probability measure ν ∈ P(T p × T d ) 3: Initial guess f 0 ∈ H 4: number of tensor terms m ∈ N * 5: Output: 6: f m ∈ H.
7: Begin: 8: n = 0 9: while n < m do 10:

for i ∈ 1, d do 11:
Find

g n := d j=1 r n,j a minimizer of g → J ν (f n + g) over g ∈ Σ i (i.e. with r n,j ∈ H 1 (T) for all 1 ≤ j ≤ d and T r n,i = 0) 12:
Set f n+1 = f n + g n .

13:

Increment n ← n + 1.

14:

end for 15: end while

16: Return f m .
In Section 4 is established the following: Theorem 5. Under Assumption 1, let f * be the minimizer of J ν in H and f m = Greedy(ν, f 0 , m) as given by Algorithm 2 for some ν

∈ P(T p × T d ), f 0 ∈ H and m ∈ N. Then f m -f * H 1 -→ m→+∞ 0 .
The interest of Algorithm 2 is that at each iteration, one only has to compute d onedimensional functions, which makes it possible to implement even if d is relatively large (say 4 < d < 10). Notice that the price to pay when going from the original problem of minimizing J ν over H to the problem [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF] is that the Euler-Lagrange equations associated to the initial problem are linear (since J ν is a quadratic functional) whereas the Euler-Lagrange equations associated to [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF] are nonlinear. This is due to the fact that the quadratic functional is minimized over a non-linear space in [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF].

In practice, a minimizer of J ν (f + g) over g = d j=1 r j ∈ Σ i is approximated through the Alternating Least Square method [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF], which is a fixed point procedure on the Euler-Lagrange equation ( 17): the r j 's are optimized one after the other, the others being fixed, repeatedly. This amounts to solving a system of one-dimensional elliptic PDEs of the form

∂ z j a j ∂ z j r j (z j ) -b j (z j )r j (z j ) = c j (z j ) (19)
with

a j (z j ) = T d-1 l =j r l (z l ) 2 θ ν (z)dz =j b j (z j ) = h =j T d-1 ∂ z h l =j r l (z l ) 2 θ ν (z)dz =j c j (z j ) = T d-1 l =j r l (z l ) ∂ z j (F ν,j θ ν ) (z)dz =j - h =j T d-1 ∂ z h l =j r l (z l ) F ν,h (z)θ ν (z)dz =j ,
where dz =j means that all variables except the j th are integrated and F ν,j denotes the j th component of 10)), for ν = ν t given by ( 8),

F ν . If for all y = (y 1 , • • • , y d ), z = (z 1 , • • • , z d ) ∈ T d , K(y, z) = Π d i=1 K i (y i , z i ) for some functions K i : T × T → R for all 1 ≤ i ≤ d (like the kernel (
(1 + λ)a j (z j ) = λ l =j r l 2 L 2 (T) + T D l =j T r 2 l (z l )K l (y l , z l )dz l K j (y j , z j )dν t (q, y) ,
which can be computed without computing F νt (z) and θ νt (z) for all z ∈ T d (which would be impossible in practice). The same holds for b j and c j .

The tensor ABF algorithm

As already explained above, the main objective of this work is to introduce a new algorithm to adapt the standard ABF approach to multi-dimensional reaction coordinates. Combining Algorithms 1 and 2, the Tensor ABF (TABF) algorithm is described in Algorithm 3 below. Note that, for the sake of clarity, it has been kept relatively simple. In particular, we haven't addressed here the question of time and space discretization. Moreover, the proofs of convergence of Algorithm 1 and Algorithm 2 also have their own interest. The convergence of Algorithm 1 is based on the so-called ordinary differential equation method [START_REF] Benaïm | Self-interacting diffusions[END_REF], and requires specific contractivity bounds. The convergence of Algorithm 2 is an adaptation of the proof of convergence of greedy algorithms [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF], the main difficulty being to deal with the zero average constraint in H.

Discussion on the results and efficiency

In view of our theoretical results, we now discuss separately the efficiency of the three algorithms: first, the usual ABF algorithm (Algorithm 1) associated with Theorems 2 and 3;

Algorithm 3 TABF algorithm 1: Input: 2: Initial condition (q 0 , z 0 ) ∈ T p × T d 3: Brownian motion (B 1 t , B 2 t ) t 0 on T p × T d 4: Regularization parameters K, λ 5: Update period T up > 0, number of updates N up ∈ N * , total simulation time T tot = T up N up 6: Number of tensor terms by update m ∈ N * 7: Output: 8: Estimated free energy A Ttot ∈ H 9: Trajectory (Q t , Z t ) t∈[0,Ttot] ∈ C [0, T tot ], T p × T d 10: Begin: 11: Set (Q 0 , Z 0 ) = (q 0 , z 0 ). 12: Set A 0 (z) = 0 for all z ∈ T d . 13: Set t k = kT up for all k ∈ 0, N up . 14: for k ∈ 1, N up do 15: Set A t = A t k-1 for all t ∈ [t k-1 , t k ). 16: Set (Q t , Z t ) t∈[t k-1 ,t k ] to be the solution of (5) with value (Q t k-1 , Z t k-1 ) at time t k-1 . 17: Set f m = Greedy(ν t k , A t k-1 , m)
given by Algorithm 2 where ν t k is given by ( 8).

18:

Set A t k = f m . 19: Return A Ttot and (Q t , Z t ) t∈[0,Ttot] .
second, the greedy tensor approximation of the free energy (Algorithm 2), associated with Theorem 5; finally, the TABF algorithm (Algorithm 3).

Enhanced sampling with the ABF algorithm

As Algorithm 1 is meant to tackle metastability issues, a natural frame to discuss its efficiency is the low temperature regime β → +∞. For Markov processes, obtaining an equivalent in this regime of the convergence rate of the law of the process toward its equilibrium is a classical topic, but the case of non-Markovian self-interacting dynamics or similar stochastic algorithms is known to be much more difficult, and there are much less results. Theorem 3 states that, in term of asymptotic variance, the efficiency of the adaptive scheme is approximately (as A * is close to A) the same as the efficiency of the importance sampling scheme based on (4). As discussed earlier, the question is thus whether (4) would be efficient. One way to quantify the interest of using ( 4) is to discuss the spectral gap of the associated infinitesimal generator: it is indeed known that the larger the spectral gap, the smaller the asymptotic variance, and the faster the convergence to equilibrium. In fact, at low temperature, the spectral gap of the overdamped Langevin process is well known to scale as exp(-βc * ) where c * is the socalled critical depth of the potential, see [START_REF] Holley | Asymptotics of the spectral gap with applications to the theory of simulated annealing[END_REF]. On the other hand, applying the results of [START_REF] Lelièvre | A general two-scale criteria for logarithmic sobolev inequalities[END_REF], we see that the spectral gap of (4) can be obtained, on the one hand, from the Poincaré inequality satisfied by the marginal law of the reaction coordinates and, on the other hand, by the conditional laws for fixed values of the reaction coordinates. The marginal law being uniform on the torus for all β, the scaling in β of the spectral gap of ( 4) is given by the scaling of the Poincaré inequality of the conditional laws, which means only the "orthogonal" metastability intervenes. The spectral gap of ( 4) then scales at most as exp(-β sup z∈T d c * (z)) where c * (z) is the critical depth of q → V (q, z). This gives a precise criterion (although difficult to use in practice) for selecting reaction coordinates: a reaction coordinate is good if sup z∈T d c * (z) < c * . In the toy problem studied in Section 6.1, for instance, sup z∈T d c * (z) = 0 (there is no orthogonal metastability).

A numerical comparison of a classical overdamped Langevin sampler and of a TABF algorithm in this case at low temperature is given in Figure 4. We can see that, in the same physical time, the adaptive process successfully visits the whole space, while the classical sampler remains trapped in its initial well. Of course this is not a fair comparison of the practical algorithms since the numerical cost of the adaptive algorithm is higher, but it illustrates the difference of the sampling rates of the continuous-time processes. In fact, for a reaction coordinate of dimension less than 3, a standard ABF method, for which the free energy is recorded on a grid rather than approximated in a tensor form, has a numerical cost in term of computational time which is very similar to a basic unbiased simulation of an overdamped Langevin process (there is also a small additional memory). This means that, in that case, a comparison at a given computational time between the two algorithms yields a result similar to a comparison at a given physical time (as in our theoretical result of Theorem 3 and in the numerical results displayed in Figure 4), i.e. the ABF method is efficient in practice at a given computational cost. From a numerical point of view, this is a well-known fact, indeed this algorithm has proven to be useful and efficient for nearly 20 years in a large number of empirical studies, see e.g. [START_REF] Darve | Calculating free energies using average force[END_REF][START_REF] Hénin | Overcoming free energy barriers using unconstrained molecular dynamics simulations[END_REF][START_REF] Fu | Extended adaptive biasing force algorithm. an onthe-fly implementation for accurate free-energy calculations[END_REF][START_REF] Comer | The adaptive biasing force method: Everything you always wanted to know but were afraid to ask[END_REF] and references within. For this reason, in the present work, only the TABF variant has been implemented and used for the numerical experiments of Section 6. However, our contribution, i.e. Theorems 2 and 3, gives the first theoretical proof of the consistency and efficiency of the self-interacting ABF method.

Tensor approximation of the free energy

First, it should be noticed that the convergence of the greedy tensor method, Algorithm 2, is proven only in the case where the problem of minimizing J ν (f + g) over single tensor terms g = d j=1 r j is exactly solved, which is in fact not the case with the Alternating Least Square method (see [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF][START_REF] Uschmajew | Local convergence of the alternating least squares algorithm for canonical tensor approximation[END_REF] for convergence results for this algorihm). Moreover, by contrast with Theorem 3 for the long-time convergence, the convergence of the tensor approximation stated in Theorem 5 is only qualitative.

That being said, let us focus on the efficiency of Algorithm 2. Whereas the precise convergence rate is not known from a theoretical point of view, it is numerically observed in numerous applications related to computational mechanics [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations: a primer[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF], rheology [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF], biology [START_REF] Barone | Efficient estimation of cardiac conductivities: A proper generalized decomposition approach[END_REF] or uncertainty quantification [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF], that this method, also called Proper Generalised Decomposition enables to give very good numerical approximations of the functions to be computed with a low number of tensor terms or iterations. In our context, we provide in Section 6 a numerical experiment where interesting non-trivial results are obtained when approximating a 5-dimensional free energy with m = 140 tensor terms, each one-dimensional function being piecewise linear on a grid with M = 30 points. Hence, the total memory cost is 5mM , orders of magnitude smaller than M 5 . We refer to the end of Section 6 (in particular Figure 10) for a more detailed numerical discussion on the efficiency of the tensor approximation.

TABF and free energy estimation

Recall that there are two distinct motivations for the TABF algorithm: first, enhancing the sampling and, second, estimating the free energy associated to a high-dimensional reaction coordinate.

Concerning the first objective, what has been said in Section 2.4.1 for the usual ABF algorithm is still true for the TABF variant as far as we consider the continuous-time processes in physical time (again, see Figure 4, which is obtained with the TABF algorithm). However, in contrast with the usual case where the free energy is recorded on a grid, since the greedy tensor approximation is numerically demanding, the numerical cost of the TABF algorithm is increased with respect to a non-biased standard MCMC simulation. In the context of the numerical experiments conducted in Section 6, which have a purely illustrative purpose, the TABF algorithm does not yield any gain with respect to a non-biased simulation in term only of sampling at a given computational cost, for two reasons:

As discussed in Section 5, to be competitive with unbiased algorithms, a more efficient implementation of the TABF algorithm would be necessary.

In the models used in the numerical experiments (e.g. 100 particles with pairwise interactions), the numerical cost for computing the forces ∇V is low, which implies that the construction of the tensor approximation becomes the bottleneck. This is not representative of the real-life test cases where evaluating the force is much more expensive.

As a consequence, the current implementation of the algorithm which has been used in Section 6 does not allow for a real comparison in term of sampling efficiency of the TABF algorithm with a standard MCMC method. This would require the interfacing of a high-performance TABF module within a general molecular dynamics software, in order to test the algorithm on applications in computational chemistry where the computation of the forces at each time step is the overwhelming limiting factor, so that the extra cost of the tensor algorithm does not impair the overall complexity. This exceeds the scope of the present work.

Concerning the second objective, namely the estimation of the free energy, the numerical experiments of Section 6 show that a simple implementation of Algorithm 3 is able to estimate the free energy of 5-dimensional highly-correlated reaction coordinates (see in particular Figure 7). Let us emphasize that, to our knowledge, there are no other implementable provably convergent algorithm to perform such a task due to the curse of dimensionality.

Proof of the long-time convergence

In the whole Section 3 we consider the ABF process (Q t , Z t , A t ) t 0 obtaind through Algorithm 1 (with N up = +∞), and Assumption 1 holds. Lemma 6. For all r ∈ N * and all multi-index α ∈ N r , there exists a constant C α > 0 such that, for all t 0, ∂ α A t ∞ C α .

Proof. Since R + t → A t is piecewise constant, we may assume that t = t k = kT up for some k ∈ N without loss of generality. Using the notation of Section 2.2, A t is then the minimizer over H of J νt defined in [START_REF] Comer | The adaptive biasing force method: Everything you always wanted to know but were afraid to ask[END_REF]. Recall that for all f ∈ H,

J νt (f ) = T d |F νt (z) -∇f (z)| 2 θ νt (z)dz.
Remark that θ νt is bounded from below uniformly in t and z by (λ + min K)/(1 + λ) > 0, and similarly all the derivatives in z of θ νt and of F νt are bounded in L ∞ (T d ) by constants which depend on K and V but not on t. The Euler-Lagrange equation associated to the minimization of J νt reads

∇ • (θ νt ∇A t ) = ∇ • (θ νt F νt ) . ( 20
)
By elliptic regularity (cf. [START_REF] Ambrosio | Lectures on elliptic partial differential equations[END_REF]), A t is thus C ∞ and, differentiating [START_REF] Fort | Convergence and efficiency of adaptive importance sampling techniques with partial biasing[END_REF], multiplying it by derivatives of ∇A t and integrating, we classically get by induction that

T d |∂ α ∇A t | 2 θ νt C α
where α ∈ N r is any multi-index for any r ∈ N * , for some constant C α > 0 which does not depend on t. Conclusion follows from Sobolev embeddings.

Theorem 2 will be a direct corollary of:

Proposition 7. Almost surely, ν t weak -→ t→∞ µ.
The proof of Proposition 7 is postponed to the end of this section. Let us prove that indeed, given the latter, Theorem 2 holds: Proof of Theorem 2. By the arguments of the previous proof, for all t ≥ 0, the function T d z → θ νt (z) is bounded and Lipschitz with constants which are uniform in t. Hence, for any ε > 0, we can find N ε ∈ N * and a finite set of points z 1 , • • • , z Nε ∈ T d such that for all z ∈ T d , there exists i z ∈ 1, N ε such that, for all t > 0, |θ νt (z iz ) -θ νt (z)| ε. The same holds for θ µ . On the other hand, according to Proposition 7, almost surely,

sup i∈ 1,Nε |θ νt (z i ) -θ µ (z i )| -→ t→∞ 0,
so that θ νt -θ µ ∞ goes to zero as t → ∞. Similar arguments enable us to obtain the same results for all the derivatives of θ νt and for F νt and all its derivatives. Note that A * is the minimizer of

J µ (f ) = T d |F µ (z) -∇f (z)| 2 θ µ (z)dz.
Let t = t n for some n ∈ N * . The associated Euler-Lagrange equations associated with the two minimization problems on A t and A * lead to

∇ • (θ µ ∇ (A t -A * )) = ∇ • (θ νt F νt -θ µ F µ -(θ νt -θ µ ) ∇A t ) . (21)
Multiplying this equality by A t -A * , integrating and using the uniform control on ∇A t established in Lemma 6 (and the lower bound on θ µ ), we get that

T d |∇ (A t (z) -A * (z)) | 2 dz -→ t→∞ 0.
More generally, differentiating [START_REF] Fu | Extended adaptive biasing force algorithm. an onthe-fly implementation for accurate free-energy calculations[END_REF], multiplying it by derivatives of A t -A * , integrating and using the uniform controls of the derivatives of A t , we obtain by induction that

T d |∇∂ α (A t -A * ) (z)| 2 dz -→ t→∞ 0 for all multi-index α ∈ N d .
The first statement of Theorem 2 then follows from Sobolev embeddings. Finally, inequality [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF] stems from the fact that J µ (A * ) J µ (A). More precisely, using that

T d ∇A(y)µ(q, y)dq = T d ∇ y V (q, y)µ(q, y)dq , we get that for all f ∈ H, J µ (f ) = J µ (f ) + |∇ y V | 2 dµ -|∇A| 2 dµ where J µ (f ) = T p ×T d ×T d |∇A(y) -∇f (z)| 2 K(y, z)dzdµ(q, y) + λ T d |∇f (z)| 2 dz.
In other words, J µ and J µ only differ by an additive constant, so that A * is the minimizer of J µ over H. Then

T p ×T d ×T d |∇A(z) -∇A * (z)| 2 K(y, z)dzdµ(q, y) 2 J µ (A * ) + 2 T p ×T d ×T d |∇A(y) -∇A(z)| 2 K(y, z)dzdµ(q, y) 2 J µ (A) + 2 T p ×T d ×T d |∇A(y) -∇A(z)| 2 K(y, z)dzdµ(q, y) 2λ T d |∇A(z)| 2 dz + 4 T p ×T d ×T d |∇A(y) -∇A(z)| 2 K(y, z)dzdµ(q, y) 2λ T d |∇A(z)| 2 dz + 4 ∇ 2 A 2 ∞ sup y∈T d T d |y -z| 2 K(y, z)dz
The rest of the section is dedicated to the proof of Proposition 7 and Theorem 3. We start with a presentation of the the so-called ordinary differential equation (ODE) method of [START_REF] Benaïm | Self-interacting diffusions[END_REF], which introduces some general ideas of the proof of Proposition 7 (although, as we will see, we are in a very simple case so that we won't really use the fully general method).

Time change and the ODE method

Following an idea of [START_REF] Benaïm | Convergence of adaptive biasing potential methods for diffusions[END_REF], we introduce the (random) time change:

τ (t) := t 0
e -βAs(Zs) ds, so that

ν t = 1 τ (t) t 0 δ Qs,Zs τ (s)ds = 1 τ (t) τ (t) 0 δ Q τ -1 (s) ,Z τ -1 (s) ds .
In other words, considering the time-changed process X t := Q τ -1 (t) , Z τ -1 (t) and its occupation measure

νt = 1 t t 0 δ Xs ds , (22)
then ν t = ντ(t) . Since, at a fixed time t 0, A t is smooth and with Lebesgue integral zero, there always exists z ∈ T d such that A t (z) = 0, so that

A t ∞ √ d/2 ∇A t ∞ (23)
where we used that √ d/2 is the diameter of T d . Together with Lemma 6, this implies that in particular, τ (t) goes to infinity with t.

Denoting S t (x) := A τ -1 (t) (z) for all x = (q, z) ∈ T p × T d , the inhomogeneous Markov process X solves the SDE

dX t = -e βSt(Xt) ∇ (V -S t ) X t + 2β -1 e βSt(Xt) dB t , (24)
where (B t ) t 0 is a standard Brownian motion on T D , obtained from (B t ) t 0 through rescaling. We denote by (L t ) t 0 its infinitesimal generator, defined by: for all

ϕ ∈ C 2 (T p × T d ) = C 2 (T D ) and all x ∈ T p × T d , L t ϕ(x) = lim h→0 E ϕ(X t+h ) | X t = x -ϕ(x) h
whenever the limit exists. Here,

L t ϕ(x) = -∇ (V -S t ) (x) • ∇ϕ(x) + 1 β ∆ϕ(x) e βSt(x) .
We denote by (P

(t)
s ) s 0 the Markov semi-group generated by L t for a fixed t. Formally, P (t) s = e sLt . For all t 0 the unique invariant measure of (P

(t) s ) s 0 is µ (see [6, Proposition 3.1]
), which is a natural consequence of the fact we consider a process interacting with its unbiased occupation measure. From Lemma 6 and the bound [START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF], we consider C 0 > 0 such that S t ∈ B C 0 for all t 0 where

B C 0 := S ∈ C ∞ T D , T D S(x)dx = 0, S C 2 (T D ) C 0 .
The principle of the ODE method is the following: for large values of the time t, the evolution of νt is slow (because of the t -1 factor in ( 22

In other words, the evolution of νt approximately follows the deterministic flow

∂ t m t = 1 t (µ -m t ) ,
which converges to µ, so that νt (hence ν t ) should also converge to µ.

In general cases of self-interacting processes, as those studied in [START_REF] Benaïm | Self-interacting diffusions[END_REF], the asymptotic deterministic flow may be more complicated (see in particular [START_REF] Benaïm | Analysis of an Adaptive Biasing Force method based on self-interacting dynamics[END_REF] for the ABF algorithm with the non-reweighted occupation measure). Here, as mentioned above, we are in a very simple case since for all t 0, µ is the invariant measure of the time-homogenous Markov process with generator L t , so that the flow is simply a relaxation toward this equilibrium. For this reason, in order to make rigorous the previous heuristic, instead of applying the technical arguments of [START_REF] Benaïm | Self-interacting diffusions[END_REF], we may use a shortcut that yields a simpler proof and more explicit estimates (allowing in particular to tackle the question of the asymptotic variance, which may be much more intricate in other cases), similarly to [7] for the ABP algorithm.

We will need some quantitative estimates. Indeed, note that, for the approximation (25) to hold, the speed of convergence of P (t) toward µ should be uniform in t, and the time evolution of A t should be controlled in some sense. As we will see below, these are direct consequences of the estimates of Lemma 6.

Preliminary estimates

For a fixed S ∈ C ∞ (T D ), consider L S defined for ϕ ∈ C ∞ (T D ) by

L S ϕ(x) = -∇ (V -S) (x) • ∇ϕ(x) + 1 β ∆ϕ(x) e βS(x) ,
which is the infinitesimal generator of the SDE dX S t = -e βS(X S t ) ∇ (V -S) (X S t )dt + 2β -1 e βS(X S t ) dB t .

Denote P S t t 0 the associated (homogeneous) semi-group and Γ S the associated carré-duchamps operator, defined for ϕ, ψ ∈ C ∞ (T D ) and all x ∈ T D by

Γ S (ϕ, ψ)(x) := 1 2 (L S (ϕψ) -ϕL S ψ -ψL S ϕ) (x) = β -1 e βS(x) ∇ϕ(x) • ∇ψ(x) ,
and Γ S (ϕ) := Γ S (ϕ, ϕ). By classical elliptic regularity arguments, if ϕ ∈ C ∞ (T D ) then P S t ϕ ∈ C ∞ (T D ), in particular C ∞ (T D ) is a core for L S , see [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF]Section 1.13]. More precisely each derivative of P S t ϕ is uniformly bounded over all finite time interval, which ensures the validity of the computations in the proofs of the next lemmas. Integrating twice by parts, it can be easily seen that for all ϕ, ψ ∈ C ∞ (T D ),

T D ϕ(x)L S ψ(x)µ(dx) = T D ψ(x)L S ϕ(x)µ(dx) ,
in other words L S is a self-adjoint operator on L 2 (µ).

Lemma 8. Let us assume that D ≥ 3. Then, there exists C 1 > 0 such that for all S ∈ B C 0 , (µ, L S ) satisfies a Poincaré inequality and a Sobolev inequality both with constant C 1 , in the sense that for all ϕ ∈ C ∞ (T D ),

ϕ 2 L 2 (µ) C 1 T D Γ S (ϕ)dµ ϕ 2 L p (µ) C 1 ϕ 2 L 2 (µ) + T D Γ S (ϕ)dµ ,
where p = 2D D-2 . Proof. For S = 0, the first inequality is the classical Poincaré inequality, which holds here since the density of µ with respect to the Lebesgue measure is bounded above and below away from zero, see [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF]Proposition 5.1.6]. As a consequence, there exists c > 0 such that for all S ∈ B C 0 and ϕ ∈ C ∞ (T D ),

ϕ 2 L 2 (µ) c T D |∇ϕ| 2 dµ ce βC 0 T D Γ S (ϕ)dµ .
Similarly, from the Sobolev inequality satisfied by the Lebesgue measure on T D [3, Section 6],

ϕ 2 L p (µ) µ 2/p ∞ ϕ 2 L p (T D ) C µ 2/p ∞ ϕ 2 L 2 (T D ) + ∇ϕ 2 L 2 (T D ) C µ 2/p ∞ µ -1 2 ∞ ϕ 2 L 2 (µ) + ∇ϕ 2 L 2 (µ) Ce βC 0 µ 2/p ∞ µ -1 2 ∞ ϕ 2 L 2 (µ) + T D Γ S (ϕ)dµ .
These inequalities, in turn, yield the following estimates:

Lemma 9. There exist C 2 > 0 such that, for all S ∈ B C 0 , t 0 and ϕ ∈ C ∞ (T D ),

P S t Πϕ L 2 (µ) e -t/C 2 Πϕ L 2 (µ) P S t ϕ ∞ C 2 min(1, t d/2 ) ϕ L 2 (µ) ∇P S t ϕ ∞ C 2 min(1, √ t) ϕ ∞ ,
with Πϕ := ϕ -T D ϕdµ.

Proof. The first estimate is a usual consequence of the Poincaré inequality, see [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF]Proposition 5.1.3]. The second one, namely the ultracontractivity of the semi-group, is a consequence of the Sobolev inequality (see [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF]Theorem 6.3.1]). The last one can be established thanks to the Bakry-Emery calculus (see [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF]Section 1.16] for an introduction), by showing that L S satisfies a curvature estimate, as we now detail. We would like to compare |∇P t ϕ| 2 and P t (ϕ 2 ). A seminal idea of the Bakry-Emery calculus is that quantities of the form Θ(P t ϕ) and P t Θ(ϕ), where Θ is some operator can be linked through the interpolation P t-s Θ(P s ϕ), s ∈ [0, t], so that Θ(P t ϕ) -P t Θ(ϕ) = t 0 ∂ s (P t-s Θ(P s ϕ)) ds. When differentiating with respect to s, we obtain quantities of the form -2P t-s Γ Θ (P s ϕ) for some operator Γ Θ , which is of a form similar to the interpolation (Θ being replaced by Γ Θ ).

More precisely, when Θ(ϕ) = ϕ 2 , then Γ Θ is the usual carré-du-champ operator, and when Θ(ϕ) = |∇ϕ| 2 we end up with

Γ ∇,S (ϕ) = 1 2 L S |∇ϕ| 2 -∇ϕ • ∇L S ϕ , for ϕ ∈ C ∞ (T D ). Writing [ϕ, ψ] = ϕψ -ψϕ, we compute Γ ∇,S (ϕ) = D i=1 (Γ S (∂ x i ϕ) + ∂ x i ϕ[∂ x i , L S ]ϕ) D i=1 β -1 e -β S ∞ |∇∂ x i ϕ| 2 -β -1 e β S ∞ |∇∂ x i (V -S)||∇ϕ||∂ x i ϕ| -β|∂ x i S|e β S ∞ |∂ x i ϕ||∇ (V -S) • ∇ϕ + 1 β ∆ϕ| -c|∇ϕ| 2
for some c > 0 which is uniform over S ∈ B C 0 . Now, following [42, Lemma 4], we want to consider the interpolation between α(t)|∇P t ϕ| 2 + (P t ϕ) 2 and α(0)P t |∇ϕ| 2 + P t (ϕ 2 ) for some α with α(0) = 0 < α(t). For fixed ϕ ∈ C ∞ (T D ), x ∈ T D and t 0, we set for all s ∈ [0, t] Ψ(s) = α(s)P S t-s |∇P S s ϕ| 2 (x) + e βC 0 P S t-s P S s ϕ 2 (x) with α(s) = (1 -exp(-2ct))/c, so that

∂ s Ψ(s) = P S t-s -2α(s)Γ ∇,S + α (s)|∇ • | 2 -2e βC 0 Γ S P S s ϕ (x) (2α(s)c + α (s) -2) P S t-s |∇P S s ϕ| 2 (x) = 0 .
In particular,

α(t)|∇P S t ϕ| 2 (x) Ψ(t) Ψ(0) = e βC 0 P S t ϕ 2 (x) e βC 0 ϕ 2 ∞
which yields the desired estimate.

Lemma 10. There exists C 3 > 0 such that for all S ∈ B C 0 , the operator R S defined for all ϕ ∈ C ∞ (T D ) by

R S ϕ = - ∞ 0 P S t Πϕdt satisfies L S R S = R S L S = Π and, for all ϕ ∈ C ∞ (T D ), R S ϕ ∞ + ∇R S ϕ ∞ + ∆R S ϕ ∞ C 3 ϕ ∞ . (26)
Proof. We follow the proof of [9, Section 5.2 and Lemma 5.1]. First, from Lemma 9 (and using the fact that P S t ϕ ∞ ϕ ∞ for all t 0),

∞ 0 P S t Πϕ ∞ dt 1 0 Πϕ ∞ dt + ∞ 1 P S t Πϕ ∞ dt 2 ϕ ∞ + C 2 ∞ 1 P S t-1 Πϕ L 2 (µ) dt 2 ϕ ∞ + C 2 ∞ 1 e -(t-1)/C 2 Πϕ L 2 (µ) dt 2(1 + C 2 2 ) ϕ ∞ ,
and similarly, using the fact that

∇P S t Πϕ ∞ = ∇P 1 P S t-1 Πϕ ∞ ≤ C 2 P S t-1 Πφ ∞ for t 1, ∞ 0 ∇P S t Πϕ ∞ dt 1 0 C 2 √ t Πϕ ∞ dt + C 2 ∞ 1 P S t-1 Πϕ ∞ dt C 2 4 + 2(1 + C 2 2 ) ϕ ∞ .
In particular R S ϕ and ∇R S ϕ are well defined in L ∞ (T D ) for ϕ ∈ C ∞ (T D ). Moreover, using the fact that, from Lemma 9,

P S t Πϕ ∞ C 2 e -(t-1)/C 2 Πϕ L 2 (µ) → 0 as t → +∞, L S R S ϕ = - ∞ 0 L S P S t Πϕdt = - ∞ 0 ∂ t P S t Πϕ dt = Πϕ .
The case of R S L S is similar: since µ is invariant for L S , L S Π = L S = ΠL S , and thus P S t ΠL S = P S t L S Π = ∂ t (P S t Πϕ) for all t 0. As a consequence, |∆R S ϕ| e βC 0 |e βS ∆R S ϕ| e βC 0 e βS ∇(V -S)

• ∇R S ϕ ∞ + Πϕ ∞ C 3 ϕ ∞
for some C 3 > 0 uniform over S ∈ B C 0 , which yields the desired result.

Lemma 11. There exist C 4 > 0 such that for all S 1 , S 2 ∈ B C 0 and ϕ ∈ C ∞ (T D ),

R S 1 ϕ -R S 2 ϕ ∞ + ∇R S 1 ϕ -∇R S 2 ϕ ∞ C 4 ∇S 1 -∇S 2 ∞ ϕ ∞ . Proof. From R S L S = Π, (R S 1 -R S 2 ) L S 1 + R S 2 (L S 1 -L S 2 ) = 0.
Multiplying this equality by R S 1 on the right, and using that R S Π = R S , we get for all

ϕ ∈ C ∞ (T D ), (R S 1 -R S 2 ) ϕ = R S 2 (L S 2 -L S 1 ) R S 1 ϕ .
Thus, from [START_REF] Hénin | Overcoming free energy barriers using unconstrained molecular dynamics simulations[END_REF],

(R S 1 -R S 2 ) ϕ ∞ + ∇R S 1 ϕ -∇R S 2 ϕ ∞ C 3 (L S 2 -L S 1 ) R S 1 ϕ ∞ C 3 e βS 2 ∇(S 1 -S 2 ) • ∇R S 1 ϕ ∞ + C 3 (1 -e β(S 2 -S 1 ) )L S 1 R S 1 ϕ ∞ C 3 e βC 0 ∇(S 1 -S 2 ) ∞ ∇R S 1 ϕ ∞ + C 3 e 2βC 0 S 2 -S 1 ∞ L S 1 R S 1 ϕ ∞ .
Conclusion follows from [START_REF] Hénin | Overcoming free energy barriers using unconstrained molecular dynamics simulations[END_REF]. Indeed, notice that S 1 -S 2 is a continuous function with integral zero, so that there exist x ∈ T D such that (S 1 -S 2 )(x) = 0, and then for all

y ∈ T D |(S 1 -S 2 )(y)| |x -y| ∇(S 1 -S 2 ) ∞ so that S 1 -S 2 ∞ √ D ∇S 1 -∇S 2 ∞ .
Lemma 12. There exists C 5 > 0 such that for all k 1 and ϕ ∈ C ∞ (T D )

R At k ϕ -R At k-1 ϕ ∞ C 5 k ϕ ∞ .
Proof. From Lemma 6, A t ∈ B C 0 for all t 0, so that Lemma 11

ν t k+1 = (1 -p)ν t k + pm.
In particular, for some c, c > 0, for all z ∈ T d and k ∈ N,

|θ k+1 (z) -θ k (z)| = p 1 + λ (q,y)∈T p ×T d K(z, y)(dm(q, y) -dν t k (q, y)) cp c k ,
where we used that A s ∈ B C 0 for all s ∈ [0, t k+1 ]. The same argument also works for the derivatives of θ νt , for F νt and its derivatives, so that for any multi-index α ∈ N d , there exists a constant C α such that for all k 1,

∂ α F k+1 -∂ α F k ∞ + ∂ α θ k+1 -∂ α θ k ∞ C α k .
Now, from the Euler equations satisfied by A t k and A t k+1 , we get

∇ • θ k ∇ A t k -A t k+1 = ∇ • ∇A t k+1 (θ k+1 -θ k ) -∇ • (θ k+1 F k+1 -θ k F k ) . (27)
Multiplying this equation by A t k -A t k+1 , integrating and using Lemma 6 and the lower bound on θ k , we get

T d |∇ A t k -A t k+1 (z)| 2 dz c k 2
for some c > 0. Next, differentiating [START_REF] Holley | Asymptotics of the spectral gap with applications to the theory of simulated annealing[END_REF], multiplying it by derivatives of A t k -A t k+1 , integrating and using by induction the previous estimates, we obtain in fact that

T d |∇∂ α A t k -A t k+1 (z)| 2 dz c α k 2
for some c α > 0 for all α ∈ N d , and Sobolev embeddings then yield the conclusion.

Proof of the main results

In this section we denote ν(ϕ) = T D ϕdν the expectation of an observable ϕ with respect to a probability measure ν.

Proof of Proposition 7. In the following, we use the same notation C for various constants. For all t 0 and ϕ ∈ C ∞ T D ,

ν t (ϕ) -µ(ϕ) = 1 τ (t) t 0 e -βAs(Zs) Πϕ(X s )ds = 1 τ (t)
t 0 e -βAs(Zs) L As R As ϕ(X s )ds .

To alleviate notations, write φs = R As ϕ. For k ∈ N and t ∈ [t k , t k+1 ), from ( 5) by Itô's formula,

φt (X t ) -φt k (X t k ) = t t k e -βAs(Zs) L As φs (X s )ds + 2 β t t k ∇ φs (X s )dB s , so that τ (t) (ν t (ϕ) -µ(ϕ)) = t 0 e -βAs(Zs) L As φs (X s )ds = φt (X t ) -φ0 (X 0 ) (28)
+ 0<t k t φt k-1 (X t k ) -φt k (X t k ) - 2 β t 0 ∇ φs (X s )dB s .
Recall that, from Lemma 6, there exists C > 0 such that almost surely τ (t) t/C for all t 0. Together with Lemma 12, we get that there exists C > 0 such that, almost surely, for all t > 0 and ϕ ∈ C ∞ T D ,

1 τ (t) | φt (X t ) -φ0 (X 0 )| + 0<t k <t | φt k (X t k ) -φt k-1 (X t k )| C ln(1 + t) t ϕ ∞ . (29)
Moreover, applying Itô's isometry,

E t 0 ∇ φs (X s )dB s 2 = E t 0 |∇ φs (X s )| 2 ds tC 2 3 ϕ 2 ∞ ,
where we used [START_REF] Hénin | Overcoming free energy barriers using unconstrained molecular dynamics simulations[END_REF]. As a consequence, there exists C > 0 such that

E |ν t (ϕ) -µ(ϕ)| 2 C t ϕ 2 ∞ ( 30
)
for all t > 0 and ϕ ∈ C ∞ T D . As in the proof of [7, Lemma 5.1], this implies the almost sure weak convergence of ν t to µ as follows. Indeed, for all r > 0, the Borel-Cantelli Lemma yields the almost sure convergence of ν exp(nr) (ϕ) toward µ(ϕ) as n → +∞, n ∈ N, so that

P ∀k ∈ N * , ν exp(n/k) (ϕ) -→ n→+∞ µ(ϕ) = 1 . (31)
Moreover, using the almost sure bounds t/C τ (t) Ct and exp(-βA s ) ∞ C for some C > 0, we get that there exists C, C > 0 such that for all t s > 0,

|ν t (ϕ) -ν s (ϕ)| C 1 τ (t) - 1 τ (s) s + |t -s| τ (t) ϕ ∞ C |τ (t) -τ (s)|s τ (t)τ (s) + |t -s| τ (t) ϕ ∞ C |t -s| t ϕ ∞ .
As a consequence, for all ϕ ∈ C ∞ T D , t → ν exp(t) (ϕ) is almost surely C ϕ ∞ -Lipschitz, and in particular

P ∀k ∈ N * , ∀t 0, |ν exp(t) (ϕ) -ν exp( tk /k) (ϕ)| C ϕ ∞ k = 1 . ( 32
)
The almost sure convergence of ν t (ϕ) to µ(ϕ) for a given ϕ then follows from

ν t (ϕ) -→ t→+∞ µ(ϕ) = k∈N * lim sup t→∞ |ν exp(t) -µ(ϕ)| C ϕ ∞ k ⊃ k∈N * sup t 0 |ν exp(t) (ϕ) -ν exp( tk /k) (ϕ)| C ϕ ∞ k ∩ ν exp( tk /k) -→ t→+∞ µ(ϕ) ,
the last event having probability 1 from ( 31) and [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF]. Considering a sequence (ϕ k ) k∈N of C ∞ T D functions that is dense in C T D , we get that

P ν t (ϕ k ) -→ t→+∞ µ(ϕ k ) ∀k ∈ N = 1 ,
so that almost surely ν t converges weakly to µ as t → +∞. This concludes the proof of Proposition 7, hence of Theorem 2.

Proof of Theorem 3. The first claim of the theorem has already been established in the proof of Proposition 7, see [START_REF] Klus | Towards tensor-based methods for the numerical approximation of the Perron-Frobenius and Koopman operator[END_REF]. Fix ϕ ∈ C ∞ T D . We have seen in the proof of Proposition 7 (see [START_REF] Jorgensen | Monte carlo simulation of differences in free energies of hydration[END_REF] and ( 29)) that, denoting φs = R As ϕ,

ν t (ϕ) -µ(ϕ) = ε t - 1 τ (t) 2 β t 0 ∇ φs (X s )dB s
for some ε t such that almost surely |ε t | C ln(1 + t)/t for all t > 0 for some C > 0. The martingale part having zero expectation, the bias is bounded as

|E(ν t (ϕ) -µ(ϕ))| 2 = |E(ε t )| 2 C 2 ln 2 (1 + t) t 2 = o t→+∞ 1 t .
In other words, the asymptotic mean-square error is only due to the asymptotic variance, which is itself only due to the martingale part of ν t (ϕ) -µ(ϕ).

Remark that, as a corollary of Theorem 2, A t -A * ∞ → 0 almost surely. Together with the uniform bounds of Lemma 6 and the weak convergence of ν t to µ, we get that where ψ = R A * ϕ solves L A * ψ = Πϕ, which reads

t τ (t) = t - √ t τ (t) + o t→+∞ (1) = 1 τ (t)
tE |ν t (ϕ) -µ(ϕ)| 2 = 2t β E 1 τ 2 (t) t 0 ∇ φs (X s )dB s 2 + o t→+∞ (1) = 2κ 2 βt E t 0 ∇ φs (X s )dB s 2 + o t→+∞ (1) = 2κ 2 βt t 0 E |∇ φs (X s )| 2 ds + o t→+∞ (1) . From Lemma 11, ∇ φs -∇R A * ϕ ∞ C 4 ∇A s -∇A * ∞ ϕ ∞ ,
tE |ν t (ϕ) -µ(ϕ)| 2 = 2κ 2 βt t 0 E |∇R A * ϕ(X s )| 2 ds + o t→+∞ (1)
1 β ∆ -∇ (V -A * • ξ) ∇ ψ = e -βA * •ξ Πϕ .
It remains to see that we get the same formula for the asymptotic variance of √ t (ν * t (ϕ) -µ(ϕ)). The computations are similar, so we only sketch the main points. Denoting τ * (t) = t 0 exp(-βA * (Z * s ))ds, as in Proposition 7,

τ * (t) (ν * t (ϕ) -µ(ϕ)) = t 0 e -βA * (Z * s ) Πϕ(X * s )ds = t 0 e -βA * (Z * s ) L A * R A * ϕ(X * s )ds = R A * ϕ(X * t ) -R A * ϕ(X * 0 ) - 2 β t 0 ∇R A * ϕ(X s )dB s .
Again, from the almost sure bound τ * (t) t/C for some C > 0, Itô's isometry and the bounds on R A * ϕ given by ( 26), we get that there exists C > 0 such that for all t > 0 and ϕ ∈ C(T D ),

E |ν * t (ϕ) -µ (ϕ)| 2 C t ϕ 2 ∞ .
The proof that this implies the almost sure weak convergence of ν * t to µ is similar to the end of the proof of Proposition 7. As a corollary,

t τ * (t) = ν * t e βA * •ξ -→ t→+∞ κ .
Finally,

tE |ν t (ϕ) -µ(ϕ)| 2 = 2t β E 1 τ 2 * (t) t 0 ∇R A * ϕ(X * s )dB s 2 + o t→+∞ (1) = 2κ 2 βt E t 0 ∇R A * ϕ(X * s )dB s 2 + o t→+∞ (1) = 2κ 2 βt t 0 E |∇R A * ϕ(X * s )| 2 ds + o t→+∞ (1) = 2κ 2 βt E τ * (t)ν * t |∇R A * ϕ| 2 e βA * •ξ ds + o t→+∞ (1) -→ t→+∞ 2κ β µ |∇R A * ϕ| 2 e βA * •ξ .

Consistency of the tensor approximation

This section is devoted to the proof of Propositions 1 and 4 and Theorem 5. In all this section, Assumption 1 holds and we write θ = θ ν , F = F ν , J = J ν for some fixed ν ∈ P(T p × T d ). Recall that F, θ ∈ C ∞ (T d ), that θ is a positive probability density on T d , and that the minimizers of J in H are exactly the minimizers of J in H, where for all f ∈ H,

J (f ) = T d |F (z) -∇f (z)| 2 θ(z)dz ,
the link between J and J begin given by ( 13). Since θ is bounded from above and below by positive constants, the weighted spaces L 2 (T d ; θ) and H 1 (T d ; θ) are equal to the flat spaces L 2 (T d ; dz) and H 1 (T d ; dz). We endow H (whose definition is given in [START_REF] Benaïm | Convergence of adaptive biasing potential methods for diffusions[END_REF], with the norm

f = (1 + λ) T d |∇f (z)| 2 θ(z)dz,
which is indeed a norm, equivalent to the usual H 1 norm from the Poincaré-Wirtinger inequality: there exists C > 0 such that for all f ∈ H,

T d f 2 (z)θ(z)dz θ ∞ T d f 2 (z)dz C θ ∞ T d |∇f (z)| 2 dz C θ ∞ θ -1 ∞ T d |∇f (z)| 2 θ(z)dz.
The scalar product associated with • is denoted by • . The choice of such a norm is motivated by the fact that, denoting by J the differential of J , then for all f, g ∈ H,

J (f ) = J (0) + J (0) • f + f 2 J (f ) • g = J (0) • g + 2 f, g .
Proposition 1 is then a direct consequence of the strict convexity of J . The unique minimizer f * of J over H being a minimizer of J , it satisfies J (f * ) = 0, which reads

∀g ∈ H , F (z) • ∇g(z)θ(z)dz = ∇f * (z) • ∇g(z)θ(z)dz
Moreover, using that J (f * ) = 0, we get that J (0) • g = -2 f * , g for all g ∈ H and then

J (0) = J (f * ) -J (0) • f * -f * 2 = J (f * ) + f * 2 .
As a consequence, for all g ∈ H,

J (g) = J (0) + J (0) • g + g 2 = J (f * ) + f * -g 2 . ( 33
) Proof of Proposition 4. Let f ∈ H and i ∈ 1, d . If J (f ) = inf{J (f + g) , g ∈ Σ i } then the result is correct since 0 ∈ Σ i is a minimizer over Σ i .
Suppose now that 0 is not a minimizer, i.e. that J (f ) > inf{J (f + g) , g ∈ Σ i } and consider a minimizing sequence (g (l) ) l∈N in Σ i such that J f + g (l) converges to inf{J (f + g) , g ∈ Σ i } as l goes to infinity. For l large enough, J f + g (l) < J (f ) so that g (l) = 0, and thus up to an extraction we suppose that g (l) = 0 for all l ∈ N. Moreover the sequence is bounded in H 1 and thus, up to the extraction of a subsequence, we suppose that it weakly converges in H 1 to some g * ∈ H. The function H g → J (f + g) being convex on H, it is weakly lower semi-continuous, so that

J (g * ) inf g∈Σ i J (f + g) .
For all l ∈ N, there exist r

(l) 1 , • • • , r (l)
d ∈ H 1 (T) such that g (l) = d j=1 r

(l) j . Since g (l) = 0, we can normalize the r (l) j 's so that r (l) j L 2 (T) = 1 for all j = i and l ∈ N. As a consequence, up to the extraction of a subsequence, for all 1 ≤ j = i ≤ d, there exists r * j ∈ L 2 (T) such that the sequence (r (l) j ) l∈N weakly converges to r * j in L 2 (T). Now, since the sequence g (l) is bounded in H and

∇g (l) 2 L 2 (T d ) = d j=1 ∂ z j r (l) j 2 L 2 (T) h =j r (l) h 2 L 2 (T) ,
we get that the sequence (r

(l) i ) l∈N is bounded in H 1 (T) (since T r (l)
i = 0 for all l ∈ N). Thus, up to the extraction of another subsequence, the sequence (r

(l) i ) l∈N weakly converges in H 1 (T) to some r * i ∈ H 1 (T) such that T r * i = 0. From [32, Lemma 2], (g (l)
) l∈N converges in the distributional sense to d j=1 r * j . Thus, g * = d j=1 r * j and since g * = 0, this implies that for all 1 ≤ j ≤ d, r * j = 0. Finally, since

∇g * 2 L 2 (T d ) = d j=1 ∂ z j r * j 2 L 2 (T) h =j r * h 2 L 2 (T)
is a finite quantity, this implies that for all 1

≤ j = i ≤ d, ∂ z j r * j L 2 (T) ≤ ∇g * 2 L 2 (T d ) h =j r * h 2 L 2 (T)
< +∞, and thus r * j ∈ H 1 (T). This implies that g * ∈ Σ i and yields the desired result.

Remark 13. The problem would be ill-posed if we were to try and minimize J (f + g -

T d g) over all g ∈ Σ := {r 1 ⊗ • • • ⊗ r d , r j ∈ H 1 (T)
for all 1 ≤ j ≤ d}. This is the reason why we introduced the condition that one of the r j 's has zero mean. Indeed, consider the situation where d = 2, f = 0 and F (z) = (a (z 1 ), b (z 2 )) for some smooth functions a, b : T → R with zero mean. Then, the minimum of J over H is 0, and only attained at f * (z 1 , z 2 ) = a(z 1 ) + b(z 2 ), which is not of the form g -T d g for some g ∈ Σ. Nevertheless, the sequence (g (l) ) l∈N * defined by: for all l ∈ N * , g (l) = r

(l) 1 r (l) 2 with r (l) 1 (z 1 ) = 1 + a(z 1 ) l , r (l) 2 (z 2 ) = l + b(z 2 )
is a minimizing sequence. Indeed, the sequence g (l) -T d g (l) l∈N * weakly converges to f * . In other words, the set g -T d g, g ∈ Σ is not weakly closed in H.

Proposition 4 proves that all the iterations of Algorithm 2 are well-defined. In the following, we consider a sequence (f n ) n∈N given by the latter and g n = f n+1 -f n for n ∈ N. The general idea of the proof of Theorem 5 is that, if the sequence (f n ) n∈N converges to some f ∞ in H, it holds that J (f ∞ ) • g = 0 for all g ∈ ∪ d i=1 Σ i , and a density argument enables to conclude. Nevertheless, remark that Σ i is not a vector space and that its elements all have null integral, so that one should be careful. In the following, we essentially adapt the arguments of [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF].

Lemma 14. For all g ∈ ∪ d i=1 Σ i and n ∈ N, |J (f n ) • g| 6 g d-1 j=0 g n+j . Proof. Let i ∈ 1, d be such that g ∈ Σ i . Let n i ∈ n, n + d -1 be such that g n i ∈ Σ i . We bound, first, |J (f n ) • g| |J (f n i ) • g| + 2 | g, f n i -f n | .
The second term of the right hand side is bounded by 2 g n i -1 j=n g j . To deal with the first one, note that, even though g n i is a minimizer of J (f n i + •) over Σ i , it is not necessarily true that J (f n i + g n i ) • g = 0, since Σ i is not a vector space. We follow the proof of [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF]Proposition 3.3]. By convexity of

t ∈ R → ψ(t) := J (f n i + g + t(g n i -g)) ,
and since t = 0 minimizes ψ(t), we get

ψ (0) ψ(1) -ψ(0) 0, which reads
J (f n i + g) • g J (f n i + g) • g n i .
Hence,

-J (f n i ) • g = -J (f n i + g) • g + 2 g 2 -J (f n i + g) • g n i + 2 g 2 -J (f n i + g n i ) • g n i + 2 g n i , g n i -g + 2 g 2 . Now, 1 being a minimizer over R of t → J (f n i + tg n i ), J (f n i + g n i ) • g n i = 0, so that -J (f n i ) • g 2 g n i 2 + g g n i + g 2 .
When applied to g = ± g n i g/ g , this inequality yields

|J (f n i ) • g| 6 g g n i ,
which concludes the proof.

Proposition 15. Let f * be the unique minimizer of J over H. Then

f n -f * -→ n→∞ 0.
Proof. As in the previous proof,

0 = J (f n + g n ) • g n = 2(1 + λ) ∇g n • (∇f n+1 -F )θ = 0
for all n ∈ N, so that,

J (f n ) -J (f n+1 ) = J (f n+1 -g n ) -J (f n+1 ) = g n 2 . (34)
In particular, since (J (f n )) n 0 is a decreasing sequence bounded from below,

n 0 g n 2 < ∞. (35)
Together with Lemma 14 and the fact J (f * ) = 0, this implies that for all g ∈ Span ∪

d i=1 Σ i , 2 f * -f n , g = J (f * ) • g -J (f n ) • g -→ n→∞ 0.
Now, for all r 1 , • • • , r d ∈ C ∞ (T), denoting by := r i ∈ C ∞ (T d ) (note that we do not have necessarily that

T d h = 0), we can write h - T d h = d j=1 r j - T r j l<j T r l l>j r l , which proves that h - T d h ∈ Span ∪ d i=1 Σ i . As a consequence, (1 + λ) T d θ ∇(f * -f n ) • ∇h = f * -f n , h - T d h -→ n→∞ 0 ,
As a consequence, the limit f ∞ of any convergent (in the weak sense in H) subsequence of (f n ) n∈N necessarily satisfies that

T d θ ∇(f * -f ∞ ) • ∇h = 0,
for any tensor product function h = r i , with r i ∈ C ∞ (T) for all 1 ≤ i ≤ d. By [12, Lemma 2.1], this implies that f ∞ = f * . On the other hand, since ( f n ) n∈N is bounded, all its subsequences admits weak convergent subsequences, and the fact they all have the same limit f * proves that the whole sequence (f n ) n∈N weakly converges in H to f * . In particular

f * -f n , f * -→ n→∞ 0. (36)
Thus, it only remains to prove that ( f * -f n , f n ) n∈N also converges to zero as n goes to infinity to obtain the strong convergence of the sequence (f n

) n∈N to f * in H. From Lemma 14, 2| f * -f n , f n | = |J (f n ) • f n | n k=0 |J (f n ) • g k | 6 n-1 k=0 g k n+d-1 j=n g j 6 nda n ∞ k=0 g k 2 .
with a n = n+d-1 j=n g j 2 . Using [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF], since n∈N a n ≤ d n∈N g n 2 < ∞, there exists an extracted subsequence (n k ) k 1 such that (n k a n k ) k≥1 converges to 0 as k goes to infinity. As a consequence, ( f * -f n k , f n k ) k≥1 goes to zero as k → ∞, and thus so does ( f n k -f * ) k≥1 by [START_REF] Lelièvre | Free energy computations: A mathematical perspective[END_REF]. Finally, from [START_REF] Lelièvre | A general two-scale criteria for logarithmic sobolev inequalities[END_REF], the sequence ( f n -f * ) n∈N is non-increasing, so that the whole sequence goes to zero. Hence the result.

Discussion and variations

For the sake of clarity, the TABF algorithm defined in Section 2.3 has been kept relatively simple, and there is obviously room for many variations or fine-tuning. We list here a few of them.

Extended ABF

Consider general reaction coordinates ξ : T D → M where M is a submanifold of R d or T d . In the Extended ABF (EABF) algorithm introduced in [START_REF] Lelièvre | Computation of free energy profiles with parallel adaptive dynamics[END_REF] (see also [START_REF] Lesage | Smoothed biasing forces yield unbiased free energies with the extended-system adaptive biasing force method[END_REF][START_REF] Fu | Extended adaptive biasing force algorithm. an onthe-fly implementation for accurate free-energy calculations[END_REF]), the state space is extended to T D × M with the addition of auxiliary variables (or fictitious particles) z ∈ M, and the potential V on T D is extended to a potential Ṽ on

T D × M as Ṽ (q, z) = V (x) + 1 2σ 2 (dist M (ξ(q), z)) 2 , ∀(q, z) ∈ T D × M,
for some small parameter σ > 0, where dist M stands for the distance on M. The reaction coordinates on the extended space are then defined by ξ(q, z) = z, which means the framework considered in the present paper is general for the EABF algorithm. If (Q, Z) is distributed according to µ Ṽ ,β , the law of Z is obtained from the law of ξ(Q) through a Gaussian convolution of variance σ 2 /β. There are several practical advantages to EABF:

In the potential Ṽ , in the case where dist M is an Euclidean distance, the z i 's for i ∈ 1, d appear in separate terms of a sum, they are not directly coupled. As a consequence, in the EABF case, µ A,β should be, in some sense, closer to the product of its marginal (namely, at equilibrium, the Z i 's should be closer to be independent) than in the non-extended ABF case. In [START_REF] Fu | Extended adaptive biasing force algorithm. an onthe-fly implementation for accurate free-energy calculations[END_REF], this was a crucial point since the density µ A,β was approximated by a tensor product. But even in our case where the approximation as a sum of tensor product is made at the level of A, we can expect this form of Ṽ to improve the approximation.

After convolution, the so-called mean-force ∇ z A is smoother than the initial mean force in the non-extended ABF case. Since it varies less, its estimation is expected to be easier.

That being said, the tensorized ABF introduced above can also be straightforwardly extended to a general ABF framework, without extended coordinates.

Non-periodic reaction coordinates

In general, M may be different from T d . If it has boundaries, for instance if M = [0, 1] d , the definition of the algorithm is the same except that the diffusion ( 5) is reflected at the boundaries of M. The proof of well-posedness and convergence of the tensor algorithm, i.e. Theorem 5, is unchanged. The proof of the long-time convergence of the idealized algorithm, i.e. Theorem 2, is similar up to technical considerations in particular to take into account boundary conditions in Section 3.2. Moreover, M may not be compact, for instance M = R d , with U satisfying suitable growth conditions at infinity. Since the Lebesgue measure has not a finite mass, a confining biasing potential has to be added to the adaptive biasing potential, see [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF]Section 1.2], in which case the law of ξ(X t ) is not supposed to converge to a uniform law (flat histogram) but to a target unimodal law on R d . Extending our long-time convergence results to this case would lead to additional theoretical difficulties and is out of the scope of the present work.

Real implementation

The algorithm really implemented for the numerical experiments in Section 6 differs from the theoretical Algorithm 3 in the following points:

1. Time and space are discretized. The SDE ( 5) is replaced by an Euler-Maruyama scheme with some time step δt and the time integral in ( 8) is replaced by a discrete sum with a time step ∆t (not necessarily small; it can be of the order of the decorrelation length of the process (Q t , Z t ) t 0 ). The one-dimensional functions in the tensor terms are restricted to be continuous piecewise linear, determined by their value on a discrete grid with q ∈ N * points, so that solving [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF] amount to solve a q × q linear system. In particular, the discrete space interpolation plays a role similar to the regularization kernel K which is no more necessary, hence is discarded.

2. In fact, it is not necessary to re-weight the occupation distribution, namely (8) can be replaced by ν t = 1/t t 0 δ (Qs,Zs) ds. In that case, ν t is expected to converge to µ V -à * ,β for some à * instead of µ V,β but the conditional law of Q given Z = z is the same for these two laws. Since the free energy only depends on these conditional laws, A t is still expected to converge to à * , that should be close to the true free energy in a sense similar to [START_REF] Cancès | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF]. This is clear in the mean-field limit of the algorithm, where no regularization is needed so that à * = A (see [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF]). It is more difficult to establish for the self-interacting ABF process, but in parallel of the present paper it has been done in [START_REF] Benaïm | Analysis of an Adaptive Biasing Force method based on self-interacting dynamics[END_REF]. We tried numerically both cases, and the results were similar. The results presented in Section 6 are obtained with the full (non reweighted) occupation distribution.

3. Instead of a single particle, in practice, several replicas of the process (5) are simulated in parallel. Denoting N the number of replicas and (Q i t , Z i t ) t 0 the i th replica, i ∈ 1, N , the total empirical distribution of the system is

νN,t = 1 N t/∆t N i=1 t/∆t k=1 δ (Q i k∆t ,Z i k∆t ) . (37)
The replicas all use the same bias A t obtained from this empirical distribution by minimizing J νN,t at times t = t k = kT up .

4. It is possible to use the tensor approximation only as a correction of the classical ABF, or more precisely of the Generalized ABF (GABF) algorithm proposed in [START_REF] Zhao | The extended generalized adaptive biasing force algorithm for multidimensional free-energy calculations[END_REF] where the bias is just a sum of one-dimensional functions. Namely, for all time t > 0, for j ∈ 1, d , let

α j,t (z j ) = T p ×T d
∂ y k V (q, y)K(y j , z j )dν t (q, y)

β j,t (z j ) = T p ×T d
K(y j , z j )dν t (q, y) .

These functions can be recorded on d one-dimensional grids and are easily updated on the fly. Denoting γ t,j (z j ) = 1 β j,t (z j )>s α j,t (z j )/β j,t (z j ) for some burn-in time s > 0, let A t,j (z j ) = z j 0 γ t,j (z)dz if the j th reaction coordinate z j lies in R and A t,j (z j ) = z j 0 γ t,j (u)du -z j 1 0 γ t,j (u)du if z j lies in T (so that, in both cases, ∂ z j A t,j is the Helmoltz projection in L 2 (dz j ) of γ t,j ). Then, at time t, in the dynamics [START_REF] Bash | Free energy calculations by computer simulation[END_REF], use the bias

∇ z A t with A t (z) = d j=1 A t k ,j (z j ) + f m (z)
where f m is a tensor approximation obtained through Algorithm 2 of the minimizer of H f → J νt k (f -d j=1 A t k ,j ), where t k = sup{t k < t, k ∈ N} is the last update time.

Some limitations and perspectives

A practical limitation observed in the algorithm is the following. Recall that d is too large to keep in memory the empirical measure on a grid by simply recording how many times each d-dimensional cell has been visited by the process, as in the classical ABF algorithm. Instead, the sequence (Z k∆t , ∇ z V (Q k∆t , Z k∆t )) k∈N is kept in memory for some ∆t > 0, and thus computing an expectation with respect to ν t has a numerical cost proportional to t. Such integrals are computed when solving the one-dimension equations [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF], which have to be solved repeatedly at each addition of a tensor term to the bias. As t grows, the update of the bias gets numerically more expensive. We list here some possible directions to address this question. The analysis of these variations is beyond the reach of the present work.

1. At the beginning of Algorithm 2, a clustering or quantization algorithm (see [START_REF] Pagès | CEM-RACS 2013-modelling and simulation of complex systems: stochastic and deterministic approaches[END_REF]) can be used to reduce the memory (Z k∆t , ∇ z V (Q k∆t , Z k∆t )) k∈ 1,tn/∆t to fewer points.

2. Another way to deal with this problem would be to use a fixed small size for the memory.

For instance, at an update time t k , the empirical measure used to define J νt k could be

νt k = t k t k-l e -βAs(Zs) ds -1 t k t k-l
δ (Qs,Zs) e -βAs(Zs) ds for some small l ∈ N * , say l = 1. In that case, in order to expect a long-time convergence of the bias, following classical stochastic algorithms, we would define the new bias as

A t k = A t k-1 + γ k f k where f k is (a tensor approximation of) a minimizer over H of H f → J νt k (A t k +f
) and (γ k ) k∈N is a positive sequence with γ k → 0 and k l=1 γ l → ∞ as k → ∞.

3.

A third way to deal with the memory management as time increases could be to use a stochastic gradient descent when solving the one-dimensional partial differential equation [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF]. In other words, when optimizing r i for some 1 ≤ i ≤ d, instead of computing averages over all steps l ∈ 1, t k /∆t , only use an approximation of ν t k by picking a random (and comparatively small) set of steps among 1, t k /∆t . Then only an estimation of the gradient of H 1 (T) r i → J t k (f + d j=1 r j ) is computed, which is exactly the settings of the stochastic gradient descent. [START_REF] Barone | Efficient estimation of cardiac conductivities: A proper generalized decomposition approach[END_REF]. Finally, since integrals with respect to the occupation measure are computed at each step of the Alternating Least Square (ALS) method, the cost can be reduced by using more sophisticated tensor approximation algorithms [46,[START_REF] Khoromskij | Tensor Numerical Methods in Scientific Computing[END_REF][START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF].

A second possible limitation is the following. Note that, as the number of reaction coordinates increases, we can expect that, at some point, the biasing scheme becomes unefficient. Indeed, by flattening the energy landscape, we replace the initial sampling problem (that was mainly restricted to low-energy regions, which form a low-dimensional manifold) by the sampling of the uniform measure on some hypercube, which is not so easy. In some sense, following the definitions of [START_REF] Lelièvre | Two mathematical tools to analyze metastable stochastic processes[END_REF], at some point, energy barriers are replaced by entropic ones (which means that, in the exploration of the space, what takes time is not crossing high energy areas but visiting all areas in a relatively high dimensional space). Moreover, the variance of the estimator (3) increases due to the exponential weights. As a consequence, as d increases, a partial biasing with V bias,t = θA t • ξ for some θ ∈ (0, 1) may be more appropriate than the full biasing (i.e. θ = 1). At the biased equilibrium, if A t = A is the true free energy, the marginal law of the reactions coordinates is thus µ (1-θ)A,β , i.e. the temperature is increased. Then the choice of θ such that this measure satisfies a Poincaré inequality with minimal constant (which means the corresponding overdamped Langevin process mixes the fastest) may not be θ = 1, see e.g. [START_REF] Fort | Convergence and efficiency of adaptive importance sampling techniques with partial biasing[END_REF].

Other possible simple variations

1. From the biased trajectory (Q t , Z t ) t 0 provided by the TABF algorithm, in order to compute expectations with respect to the target Gibbs measure µ = µ V,β , an alternative to the reweighting step (3) is the following. Remark that only the Z variable is biased, so that for all z ∈ T d , the conditional expectations T p f (q, z)µ(q, z)dq/ T p µ(q, z)dq can be estimated without re-weighting. On the other hand, the marginal law of Z is estimated by exp(-βA Ttot )/ T d exp(-βA Ttot (z))dz. 2. The bias update period T up and the number m of tensor products added at each update in Algorithm 3, instead of having fixed values, could be adaptively chosen. For instance, the bias could be updated when the histogram of the reaction coordinates have reached some stability, and m could be the lowest integer n ∈ N such that J

ν (A t k + f n-d ) - J ν (A t k + f n ) ε for some threshold ε > 0.
3. A time-dependent weight in the definition (8) of ν t (or, in practice, in [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF] for νN,t ) can be added in such a way that old samples have less influence than new ones since they are more biased toward the initial distribution.

4. The regularization kernel K and parameter λ may depend on time. Indeed, as time goes, the size of the sample increases. Since the problem of minimizing J νt is in practice solved on a finite dimension space, for a time large enough the regularization is actually not necessary anymore and the minimization problem with K(z, y) = δ z (y) and λ = 0 is well-posed.

5. For k ∈ N + , denote by J λ k the function given by [START_REF] Brandt | Machine learning of biomolecular reaction coordinates[END_REF] for some λ > 0 with ν = ν t k . Rather than setting A t k to be the minimizer of J λ k , we can set it to be A t k-1 + f where f is the minimizer of

H f → J 0 k A t k-1 + f + λ T d |∇ z f (z)| 2 dz , (38)
the difference being that A t k-1 does not appear in the last regularization term any more. Note that, when λ = 0, there is no difference. When λ > 0, the theoretical results of Section 4, i.e. the well-posedness of the tensor approximation, can be straightforwardly adapted. The long-time behaviour study of Section 3 should be similar, although a bit more troublesome since A t k would not depend only on the empirical distribution η t k but also on the previous bias A t k-1 . On the other hand, remark that 0 is a minimizer of (38) if and only if A t k is a minimizer of J 0 k . As a consequence, the long-time limit of A t should be the minimizer of J µ with λ = 0 which, in view of ( 12), advocates for this alternative form of cost function.

6. Since the bias is stored in memory in a tensor form, it is possible to use at some times a compression algorithm (see [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]) to reduce the number of tensor terms, if needed.

Numerical experiments

Let us fix some details and parameters that will hold for the different examples below. In this section, the modifications discussed in Section 5.3 are enforced. The one-dimensional functions r n,j are stored for all n ∈ N and 1 ≤ j ≤ d on a discrete grid with M = 30 points, so that the minimization of functions of the form d j=1 r j → J νt (f + d j=1 r j ) is restricted to tensor products of one-dimensional continuous piecewise linear functions on this grid, and the Euler-Lagrange equations [START_REF] Espig | On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations[END_REF] are replaced by M × M linear systems. This discretization replaces the regularization by a kernel K, which is no more necessary. The process (5) is discretized with a time-step δt = 25.10 -5 , while the time integral in the empirical measure ν t defined in (8) is discretized with a time-step ∆t = 20δt. In other words, the reaction coordinates and the associated local mean forces are recorded in memory only every 20 steps of the Euler scheme. Moreover, N independent replicas of the processes are run in parallel and the occupation measure used to defined the bias is νN,t given by [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF]. The update times t k of the bias are fixed at t k = kT , with T a multiple of ∆t and the number of tensor terms g n added at each update time is fixed with value m.

A low dimensional example

We start to test the method on a toy model, with N = 30 replicas, a bias update period of T = 100∆t, a regularization parameter λ = 10 -5 , and m = 8 tensor products added at each update. The reaction coordinates are Euclidean coordinates, more precisely ξ(x) = (x 1 , x 2 ), so that we don't introduce any additional extended coordinate.

The dimensions are D = 3, d = 2, particles start at (0, 0, 0) and V (x 1 , x 2 , x 3 ) = -sin(3x 1 ) sin(x 2 ) cos(x 3 -1) + cos(3x 2 + 2)(0.5 + cos(x 3 -2)) + 2 sin(2x 1 + 0.5) cos(x 3 ) -5 cos(x 1 ) cos(x 2 ) cos(x 3 + 1) .

This potential has the following properties: it is not a tensor product and yields a metastable process but, since V (x 1 , x 2 , x 3 ) = ψ(x 1 , x 2 ) cos(x 3 + ϕ(x 1 , x 2 )) for some functions ψ and ϕ, there is no metastability in the orthogonal space for fixed x 1 , x 2 .

The results are given in Figures 1 and2 (for β = 1) and 3 and 4 (for β = 5). In both cases, the theoretical free energy is successfully computed and the histograms of the reaction coordinates is eventually flat. This is a bit slower with the inverse temperature β = 5, since the initial metastability is very strong. As can be seen in Figure 4, at that temperature and in the same times, a non-biased process is stuck in its initial well.

Polymer ring in solvent

We now consider a system inspired from [START_REF] Alrachid | Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection[END_REF]. The system is constituted of two types of particles, solvent particles, and polymer particles. The polymer particles interact through a potential made precise below to form a ring. The reaction coordinates are the bond lengths between consecutive polymer particles. This gives a large dimensional problem, for which the total dimension and the number of reaction coordinates are easily prescribed, and moreover where the reaction coordinates should exhibit some correlations (if it wasn't the case, then the TABF algorithm would not give different results than GABF [START_REF] Zhao | The extended generalized adaptive biasing force algorithm for multidimensional free-energy calculations[END_REF]).

In a two-dimensional periodic box, we consider D/2 = 100 particles among which d (labeled from 1 to d) form a polymer and the others are solvent particles. The length of the box is L = D/2, to ensure a concentration independent from D. Each pair of particles that involves at least one solvent particle interacts through the purely repulsive WCA pair potential, which is the Lennard-Jones potential truncated at its minimum, namely

V W CA (r) = ε1 r r 0 1 + σ r 12 - σ r 6
where r denotes the distance between the two particles, ε = 1, σ = 0.5 and r 0 = 2 1/6 σ. Each pair of consecutive particles in the polymer ring (where the d th and first particles are considered to be consecutive, closing the loop) interacts through a double well potential

V DW (r) = h 1 - (2r -2r 1 -ω) 2 ω 2 2 ,
where r 1 = r 0 , ω = 1 and h = 3. The minimum of V DW is attained at r = r 1 (compact state) and r = r 1 + w (stretched state). Finally, each triplet of consecutive particles in the polymer also interacts through the angle θ they form with the potential

V A (θ) = 1 2 (cos(θ) -cos(θ d )) 2
with an equilibrium angle θ d = π(1 -2/d) that ensures that the total angular potential is minimized when the polymer particles form a regular polygon. There are d reaction coordinates, which are the distances between two consecutive polymer particles. Following Section 5.1, the interaction between an extended reaction coordinate z and the corresponding distance r in the system is given via the extended potential Figure 5: The slow motions of the system are the transitions of each bond between two consecutive particles of the polymer from its compact state to its stretched state.

The total energy of the (extended) system is thus, for (q, z) ∈ (LT) D × [ξ min , ξ max ] d ,

V (q, z) = D/2 i=d+1 j<i V W CA (|q i -q j |) + d i=1 V E (z i , |q i+1 -q i |) + d i=1 V DW (|q i -qi+1 |) + d-1 i=1 V A arccos qi+1 -q i |q i+1 -q i | • qi+2 -qi-1 |q i+2 -qi+1 | ,
where qi = q i for all i ∈ 1, d and qd+j = q j for j = 1, 2. Initially, the polymer is in a compact state, i.e. the distances between two consecutive of its particles are at distance r 1 , the angles are θ d and all the extended variables (z i ) i∈ 1,d are at 0. For this model, we use the variant described in point 4 of Section 5.3 namely, following the GABF algorithm, we keep in memory one dimensional free energies on a grid and we use the tensor approximation as a correction of this initial guess. There are N = 50 replicas, the update period, regularization parameter, and inverse temperature are respectively T = 10 4 ∆t, λ = 0.05 and β = 1, and at each update, m = 4d tensor products are added.

The free energy is expected to be close to a sum of one-dimensional double well potentials, with minima attained at points close to 0 and 1. Nevertheless the angular force should favor configurations where the consecutive distances in the polymer are close. This fact cannot be grasped by the GABF algorithm alone, for which reaction coordinates are treated independently one from the others.

The results are presented in Figure 6 for d = 3 and Figures 7 and 8 for d = 5. In Figure 6, we see that indeed the one-dimensional free energies recovered by the GABF algorithm have two wells approximately at 0 and 1, and that the non-independent part of the free energy has the following effect: when z 3 = 0, the well (0, 0) is favored, when z 3 = 1 the same goes for (1, 1), and when z 3 is intermediate the landscape is flatter and the two wells {z 1 = z 2 = x} with x ∈ {0, 1} are favored with respect to the wells (0, 1) and (1, 0). The result is similar in Figure 7, even though the quality of the estimation is lower for z 3 = z 4 = z 5 = 0.5, which is to be expected as this lies in a very low probability area (since 0.5 is the saddle point of the two well potential). This shows that the TABF algorithm is able to recover non-trivial dependency structures between reaction coordinates.

For d = 5, the free energy is eventually approximated with 140 tensor terms (7 updates, adding 20 terms each), which means 5 × 140 = 700 one-dimensional functions have been stored, each represented by M = 30 numbers. This is orders of magnitude below the cost 30 5 = 2, 43 × 10 7 required to store a 5-dimensional function on a grid of the same precision. Moreover, since the total number of time steps of the simulation is of order 10 6 , most of the points of the 5-dimensional grid have never been visited during the whole simulation, so we would'nt have any estimation of the free energy with a classical ABF algorithm.

We investigate further the quality of the tensor approximation by computing the evolution of the square error on the free energy as a function of the number of tensor terms. To this aim, first, using the same data as in Figure 7, in Figure 9 we plot A m -A m f 2 2 with m f = 140 where A m is the approximation obtained with the m first tensor terms constructed by the algorithm (the GABF one-dimensional free energies are not taken into account). We see a fast decay at m ∈ {0, 20, 60}, which corresponds to different update times (recall 20 terms are added at each update), and thus to different empirical distributions and different GABF estimations. This makes this figure difficult to interpret. For this reason, second, in order to focus only on the tensor approximation step without the adaptive bias process or the GABF one-dimensional free energies, we consider N = 50 independent long standard (non-biased) trajectories with 5.10 6 time steps followed by a unique tensor approximation of the free energy with Algorithm 2. Denoting again by A m the approximation obtained with m tensor terms, we plot A m -A m f

eee

  -βA * (Zs) e βA * (Zs) ds + o t→+∞ -βAs(Zs) e βA * (Zs) ds + o t→+∞ -βAs(Zs) e βA * (Zs) ds + o t→+∞ βA * •ξ =: κ almost surely. As a consequence,

  ϕ(X s )| 2 e -β(As(Zs)-A * (Zs)) ds + o t→+∞ t |∇R A * ϕ| 2 e βA * •ξ + o t→+∞ A * ϕ| 2 e βA * •ξ . In other words, the asymptotic variance is 2 β µ e βA * •ξ |∇ψ| 2 µ e βA * •ξ

Figure 1 :

 1 Figure 1: For β = 1, left and middle: estimated free energy respectively at t = 1 and t = 30. Right: theoretical free energy.

Figure 2 :

 2 Figure 2: For β = 1, cumulated histograms of the reaction coordinates at t = 30 respectively for the TABF algorithm (left) and a non-biased process (right).

Figure 3 :

 3 Figure 3: For β = 5, left and middle: estimated free energy respectively at t = 10 and t = 100. Right: theoretical free energy.

Figure 4 :

 4 Figure 4: For β = 5, cumulated histograms of the reaction coordinates at t = 100 respectively for the TABF algorithm (left) and a non-biased process (right).

2 for

 2 δ = 0.01. The scaling ensures that the minimum of V E (z, r) + V DW (r) is attained at z = (r -r 1 )/w ∈ {0, 1}. Moreover, in line with Section 5.2, the extended variable is confined in [ξ min , ξ max ] d by orthogonal reflection at the boundary, with ξ min = -0.2 and ξ max = 1.2.

2 2 2 2 / A m f 2 2

 222 with m f = 200. The result, displayed in Figure 10, shows a fast convergence. More quantitatively, we observe that the relative error A m -A m f is below 20% with m = 19, below 10% with m = 53 and below 1% with m = 123.

Figure 7 :

 7 Figure7: For d = 5, up: estimated free energy as a function of (z 1 , z 2 ) when z 3 , z 4 , z 5 are, respectively, (0, 0, 0) (left) (0.5, 0.5, 0.5) (middle) and (1, 1, 1) (right). Bottom: idem but without the independent, one-dimensional parts given by the GABF algorithm.

Figure 8 :

 8 Figure 8: For d = 5, cumulative one-dimensional histograms of the five reaction coordinates at t = 50 for the TABF algorithm (up) and for a non biased process (bottom).

Figure 9 :

 9 Figure9: Square error on the free energy as a function of the number of tensor terms (same simulation as Figure7).

Figure 10 :

 10 Figure 10: Square error on the free energy as a function of the number of tensor terms (long unbiased trajectories followed by Algorithm 2).

  )). Hence, for 1 s t, in principle, it holds that S u S t for u ∈ [t, t + s], so that

	νt+s =	t t + s	νt +	s s + t	1 s	t	t+s	δ Xu du	t t + s	νt +	s s + t	µ.

  applies. It remains to obtain a bound on ∇A t k -∇A t k-1 ∞ . In this proof, to simplify the notation, we write θ k = θ νt k andF k = F νt k . Denoting by

	m :=	t k+1 t k t k δ (Qs,Zs) e βAs(Zs) ds t k+1 e βAs(Zs) ds	and p :=	t k+1 t k t k+1 0	e βAs(Zs) ds e βAs(Zs) ds	,
	it holds that

Acknowledgements

This work was supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492 and under the European Union's Horizon 2020 Research and Innovation Programme, ERC Grant Agreement number 810367, project EMC2. It was also supported by the ANR JCJC project COMODO (ANR-19-CE46-0002) and the ANR Project EFI (ANR-17-CE40-0030) of the French National Research Agency. The authors would like to thank the associate editor and the referees for their work.

Figure 6: For d = 3, up: estimated free energy as a function of (z 1 , z 2 ) when z 3 is, respectively, 0 (left) 0.5 (middle) and 1 (right). Middle: idem but without the independent, one-dimensional parts given by the GABF algorithm. Bottom: one-dimensional potential given by the GABF algorithm for z 1 (left) z 2 (middle) and their sum (right).