Bayesian Nonparametric Inference in McKean-Vlasov models
Richard Nickl Grigorios A. Pavliotis Kolyan Ray



	

Bayesian Nonparametric Inference in McKean-Vlasov models
    Richard Nickl  [Department of Pure Mathematics and Mathematical Statistics, University of Cambridge, email: nickl@maths.cam.ac.uk]    Grigorios A. Pavliotis [Department of Mathematics, Imperial College London, email: g.pavliotis@imperial.ac.uk, kolyan.ray@imperial.ac.uk]    Kolyan Ray ^† 
    
=====================================================================================================================================================================================================================================================================================================




We consider nonparametric statistical inference on a periodic  interaction potential W from noisy discrete space-time measurements of solutions ρ=ρ_W of the nonlinear McKean-Vlasov equation, describing the probability density of the mean field limit of an interacting particle system. We show how Gaussian process priors assigned to W give rise to posterior mean estimators that exhibit fast convergence rates for the implied estimated densities ρ̅ towards ρ_W. We further show that if the initial condition ϕ is not too smooth and satisfies a standard deconvolvability condition, then one can consistently infer the potential W itself at convergence rates N^-θ for appropriate θ>0, where N is the number of measurements. The exponent θ can be taken to approach 1/2 as the regularity of W increases corresponding to `near-parametric' models.





§ INTRODUCTION


We investigate the problem of conducting statistical inference in a class of dynamical systems whose state at time t ∈ [0,T] is described by the probability density ρ=ρ_W(t, ·) solving the nonlinear partial differential equation (PDE)

    ∂ρ/∂ t   = Δρ + ∇· (ρ∇ (W ∗ρ)) 
    ρ(0, ·)   = ϕ,

known as the McKean-Vlasov equation, where Δ, ∇·, ∇ denote the Laplacian, divergence and gradient operators, respectively. Further ϕ is the probability distribution of the initial state and W is an interaction potential with Lipschitz gradient. In the present paper we consider this equation on the d-dimensional torus ^d, with periodic boundary conditions.

This PDE is a non-linear and non-local Fokker-Planck type equation and plays a fundamental role in a variety of scientific application areas  ranging from opinion dynamics <cit.> and other models in the social sciences <cit.> to mathematical biology 

<cit.>, sampling and optimization <cit.>, plasma physics <cit.>, fluid mechanics <cit.>, as well as statistical physics <cit.> and kinetic theory <cit.>. 

The problem we shall study here is how to recover the interaction potential W from a discrete statistical measurement of the evolution of the macrosopic state ρ of the system over time. We shall take a non-parametric approach where a possibly high- or infinite-dimensional model is postulated for the function W, arising from a Gaussian process prior Π, which then is updated via Bayes rule from `regression type' data (Y_i, t_i, X_i)_i=1^N of the form

    Y_i = ρ_W (t_i, X_i) + ε_i,  ε_i ∼^iid N(0,1),   i=1, …, N,

where the t_i, X_i are sampled discretely from the the time-space cylinder 𝒳 = [0,T] ×^d. This approach is frequently taken in applications, e.g., <cit.>, and fits into the general paradigm of Bayesian non-linear inverse problems with PDEs <cit.>,  which is amenable to an algorithmic treatment via MCMC and numerical PDE methods even in high- and infinite-dimensional models <cit.> – see after (<ref>) below for more details. It also allows us to take advantage of recent theoretical developments in the field, e.g., <cit.> and especially <cit.>, but before we do so let us shed some more light on how such measurements may arise in concrete physical situations, explaining also the interpretation of the function W as driving interactions.


A basic interacting particle model postulates the simultaneous evolution in time t ∈ [0,T] of n particles X_t^i in a d-dimensional state space ^d, solving the coupled system of n stochastic differential equations (SDEs)

    dX_t^i = -1/n∑_i ≠ j∇ W(X_t^i-X_t^j)dt + √(2) dB_t^i,    i=1, …, n,

where the B^i_t are independent d-dimensional Brownian motions, and where chaotic initial conditions are assumed, i.e., the X_0^i are all started independently, drawn randomly from the given initial distribution ϕ. This setup can be considered as well in ^d, in particular if we add a drift -∇ U(X^i_t)dt arising from a (known) confining potential U to the  SDE (<ref>). There has been a sequence of interesting recent papers trying to address the problem of inferring W from measurements in models such as (<ref>), see <cit.> – with earlier references being <cit.>. Key statistical aspects of this task, such as dealing with the multi-dimensional setting d>1 and clarifying when the potential W is indeed statistically identifiable, still remain broadly open, however. The important reference <cit.> shows how the particle densities can be estimated empirically from the data, but as our results indicate, one can obtain faster convergence rates than in <cit.> in important physically relevant settings (see below for more discussion). Furthermore, the hypotheses imposed in all the preceding references to identify W are typically implicit and can generally not be verified for natural parameter spaces. Note also that recently developed proof techniques for non-parametric inference in multi-dimensional diffusion models (see <cit.> and references therein) do not apply to the model (<ref>) when the number n of particles is large as the dimension nd of the underlying state space then diverges rapidly.




To exploit the information present in measurements of such systems in the mean field limit (n →∞), one ought to take advantage of the rich mathematical theory that describes the macroscopic behaviour of particles interacting according to (<ref>). In particular, it is well known that under mild assumptions on the interaction potential, and with chaotic initial data, the propagation of chaos property holds– we refer to <cit.> for a classical reference, to <cit.> for particularly clear proofs, and <cit.> for important recent developments. Specifically, it is shown in <cit.> that if one considers the marginal distributions P_t^k,n of any fixed number of k=o(n) particles, then in the mean field limit n →∞ the laws P_t^k,n are approximated in relative entropy by a product measure P^k_t whose coordinate distributions have  probability density ρ_W solving the PDE (<ref>). Moreover, ρ=ρ_W is the density of the law Law (X_t)  where X_t is the solution of the nonlinear (in the sense of McKean) SDE

    dX_t = -∇ W ∗ρ(t, X_t)dt + √(2) dB_t,  t ≥ 0,   X_0 ∼ϕ.

The existence of the solution to the last equation is not difficult to prove, in particular when the state space is the torus, by a fixed point argument in path space, see, e.g., Ch.I in <cit.>. The intuition is that as n increases, the distribution of any fixed particle is determined by the local dynamics of a Brownian motion with drift term reacting to the local density ρ of nearby particles by convolution with the gradient of the interaction potential. 

Moreover, we can regard the laws of each distinct (fixed) pair chosen among k=o(n) particles as approximately independent. One can then record the `relative frequencies'  of the proportion of these k particles spent in each of N `bins' of a dissection of the time-space cylinder [0,T] ×^d. This allows one to obtain approximate (`histogram' type) measurements of the particle density ρ_W, and by the approximate independence of the k marginals we can justify the measurement model (<ref>). We believe that the results obtained here also strongly inform the non-iid models (<ref>), (<ref>), possibly even by using the non-asymptotic approximations in <cit.> combined with results in <cit.> to bound the (one-sided) Le Cam distance for equivalence of statistical experiments. These questions are of independent interest and will be investigated elsewhere – in the present article we content ourselves with aggregating the information in the data to a measurement (<ref>) of the particle densities. 



Once we have adopted the paradigm of Bayesian nonlinear regression with data (<ref>), the question of whether consistent statistical inference on W, ρ_W is feasible can be determined by studying analytical properties of the parameter to solution map W ↦ρ_W of the underlying PDE – this was first shown in <cit.> and a general theory is laid out in <cit.>. For the present non-linear Fokker-Planck equation (<ref>) we exploit ideas from <cit.>  that allow to realise ρ_W as a limiting fixed point of sequences of suitable linear parabolic problems, as well as a `Girsanov type' estimate in relative entropy for solutions to linear and nonlinear Fokker-Planck equations


<cit.>; this bound permits us to replace the  arguments for linear PDEs from Sec. 2.1.1 in <cit.>. Our key analytical results consist, first (Lemma <ref>), of proving Lipschitz continuity of the non-linear forward map W ↦ρ_W on bounded subsets of C^2(^d);

    ρ_W_1 - ρ_W_2_L^2(𝒳)≲∇ (W_1-W_2) ∗ρ_W_2_L^2(𝒳)≲W_1-W_2_H^-β
 where β>0 is 

 a measure of the regularity of the initial condition ϕ∈ H^β. Using the theory from <cit.> this implies that Gaussian process priors for interaction potentials W give rise to sensible non-parametric models for the relevant regression functions ρ_W. In particular, this enables us to obtain fast convergence rates for recovery of ρ_W by Bayesian plug-in estimates ρ_W̅_N arising from the posterior mean W̅_N = E^Π[W|(Y_i, t_i, X_i)_i=1^N]. We then tackle the inverse problem and prove a partial converse to the first inequality in the last display, specifically a stability estimate

    (W_1-W_2)∗ρ_W_2_L^2(𝒳)≲ρ_W_1- ρ_W_2_+

for a norm ·_+ that involves partial derivatives in the space and time variables, see (<ref>) for details. When W_1 has a bound on the number of non-zero Fourier modes, this further implies novel injectivity (or `inverse stabilty') estimate for the entire map W ↦ρ_W as soon as deconvolving the `factor' ρ_W_2 is possible. This permits us to provide concrete verifiable hypotheses on only the initial condition ϕ under which full inference on the interaction potential becomes feasible at polynomial convergence rates N^-θ, θ>0, as soon as W is Sobolev smooth. The Bayesian posterior mean estimators that attain these convergence rate can be implemented by MCMC techniques, and unlike frequentist methods studied in earlier references, our approach does not require or postulate an explicit identification formula for W. Rather, the aforementioned stability estimates are sufficient to guarantee recovery, which is one of the main reason why Bayes methods are attractive in nonlinear inverse problems arising with PDEs and data assimiliation tasks, see <cit.> and references therein.




Our proofs exploit information in the observations available at finite (in fact, short) time horizons, and are not based on any kind of ergodicity properties of the particle system.  
 Several of the recent papers on inference for the McKean SDE cited earlier develop inference methodologies based on observing particles over a long time interval, and rely on the ergodic properties of the dynamics in the limit as T →∞. This approach is not necessarily well-suited for the mean field dynamics. Indeed, it is well known <cit.> that the McKean-Vlasov dynamics can exhibit phase transitions, i.e. multiple stationary states can exist. 

In particular, for the dynamics on the torus in the absence of a confining potential, the uniform measure is always a stationary state for the mean field dynamics, and clearly no information about the interaction potential can be deduced from it. For example, the spectral theoretic approach that was proposed in <cit.> and that is based on linearising the McKean-Vlasov operator around the invariant measure, does not work in this case, since the linearised Mckean-Vlasov operator is simply the Laplacian. But even when the McKean SDE has a unique informative invariant measure, long term horizon asymptotics may not be useful as the invariant measures are generally smooth and hence not easily deconvolvable.













 







§ MAIN RESULTS




 §.§ Notation


We will denote by C(^d) the space of (bounded) continuous functions on ^d while L^2(^d) denotes the usual Hilbert space of square integrable functions for Lebesgue (probability) measure dx on ^d. The  spaces H^α(^d), C^α(^d) consist of all functions that have partial (in the former case, weak) derivatives up to order α∈ℕ defining elements of L^2(^d), C(^d), respectively, and we set C^∞(^d) = ∩_α>0 C^α(^d), C^0(^d)=C(^d) by convention. For regularity estimates for the PDE (<ref>), it is also convenient to introduce the space W^α, ∞(^d) of functions in C(^d) whose weak partial derivatives up to order α lie in the space L^∞(^d) of bounded functions. We have C^α(^d) ↪ W^α, ∞(^d) where ↪ means a norm-continuous imbedding. We can further define Sobolev spaces H^α(^d)=(H^-α(^d))^* also for α<0 as the topological dual space. The convolution ∗ of functions and measures is defined as usual, for instance if f,g ∈ C(^d) then f ∗ g(x) = ∫_^d f(x-y)g(y)dy, x ∈^d. 

We further define function spaces on the time-space cylinder 𝒳 =[0,T] ×^d, such as the Hilbert space L^2(𝒳, λ) where λ=λ_T is the uniform probability measure on . For B a normed space we also use standard parabolic PDE notation for the function spaces L^p([0,T], B) of maps H:𝒳→ whose norms H(t,·)_B lie in L^p([0,T]). The spaces H^m([0,T], B) then denote those functions H ∈ L^2([0,T], B) whose (weak) time derivatives (∂^j/∂ t) H  for all 0 ≤ j ≤ m lie in L^2([0,T], B), with corresponding Hilbert space norm.  We also write 𝒫=𝒫(^d) for the set of all (Borel-) probability measures on ^d.



 §.§ Gaussian process regression in the McKean-Vlasov model


For an initial condition ϕ∈ H^β, consider the time marginal densities ρ=ρ_W providing the unique periodic solutions to the non-linear parabolic PDE (<ref>) on the time-space cylinder 𝒳  = [0,T] ×^d, see Theorem <ref>. Consider N independent and identically distributed (iid) observations (t_i,X_i,Y_i)_i=1^N arising from the model (<ref>), where the (t_i, X_i) are drawn iid from the uniform distribution λ on [0,T] ×^d, independently of the noise ε_i. We denote by D_N = ((t_i,X_i,Y_i))_i=1^N  the observation vector, with law P_W^N. As explained in the introduction, such regression-type forward models are widely used in (Bayesian) inverse problems <cit.> and are approximately justified in interacting particle models (<ref>) by the propagation of chaos phenomenon.

We devise a Bayesian model for the solutions ρ_W of the PDE (<ref>) in the non-linear regression framework (<ref>), placing a prior on the interaction potential W, our quantity of interest, while treating the initial condition ϕ as given. Another common approach in data assimilation with such time evolution equations would be to assign a prior to the initial condition ϕ – see for instance the recent article <cit.> in  the 2D Navier-Stokes model.



We employ Gaussian process priors <cit.> for the interaction potential – these will be supported on a separable normed linear space (𝒲, ·_𝒲) satisfying the imbedding

    𝒲↪ C^2(^d) ∩{W: ∫_^dW(x)dx =0}.

Note that the potential W is identifiable only up to constants as the gradient ∇ W relevant in models (<ref>),(<ref>),(<ref>) annihilates constants, so fixing ∫ W=0 is natural for identifiability reasons. The separability is required to apply Gaussian process techniques from <cit.> – we can think of 𝒲 to equal a Sobolev space H^α+1(^d) for α>1+d/2 or, to obtain sharper results, we can also take periodic Besov spaces 𝒲 = B^α_∞ 1(^d), α≥ 2 (<cit.>, p.370). Concrete examples will be discussed below, but for now we let Π be any Borel probability measure on some space 𝒲 satisfying (<ref>). 

The posterior distribution obtained from such a prior Π and data (<ref>) is given by Bayes' formula via standard arguments <cit.>;

    dΠ(W|D_N) ∝ e^ℓ_N(W) dΠ(W);  ℓ_N(W) = -1/2∑_i=1^N |Y_i - ρ_W(t_i, X_i)|^2,  W ∈𝒲.

Posterior draws w ∼Π(·|D_N) can be approximately calculated from Markov chain Monte Carlo (MCMC) techniques, for instance by the pCN, ULA or MALA algorithm, following the now well established paradigm of Bayesian inversion and uncertainty quantification in PDE models, see <cit.> and also Section 1.2.4 in <cit.> for an overview. Each MCMC step requires one numerical solution of the non-linear Fokker-Planck equation (<ref>), for which a variety of methods exists, see <cit.> for recent references. We can then form ergodic averages (1/K) ∑_k=1^K w_k of these MCMC draws w_k to approximately compute the posterior mean estimate W̅ = E^Π[W|D_N] of the interaction potential W, which in turn delivers estimates also of the particle densities ρ_W̅ by one further numerical solution of the Fokker-Planck equation (<ref>) with interaction potential W̅. For a concrete example and further discussion on applying MCMC method for instance to the Keller-Segel model for chemotaxis (which is a non-linear, non-local PDE of the form (<ref>), <cit.>), see for instance <cit.>.

One can also attempt to minimise the negative log-posterior over 𝒲 (amounting to a MAP estimate, or Tikhonov regulariser), but we note that the problem is non-convex due to the nonlinearity of W ↦ρ_W, hence MCMC may be more robust to the presence of local optima in the criterion function. Computational guarantees for gradient based MCMC methods can be obtained in principle following ideas from <cit.> combined with some of the stability estimates from the proofs that follow. But this is beyond the scope of the present paper – we shall focus here on analysing the statistical properties of the posterior distribution Π(·|D_N) under the objective (`frequentist') assumption that an actual ground truth potential W_0 generated the interacting particle system and thus data from equation (<ref>). Therefore the statistical analysis is under the law P_W_0^N of the data vector D_N, and →^P_W_0^N will denote convergence in P_W_0^N-probability. 






 §.§ Posterior contraction rates for the particle densities


We first establish contraction rates for posterior distributions (<ref>) arising from rescaled and possibly projected Gaussian process priors over ^d. We start with a mean-zero Gaussian process (V(x): x ∈^d) with reproducing kernel Hilbert space (RKHS)  whose law Π_V satisfies the following assumption. For definitions and background material on Gaussian processes and their associated RKHS, see Chapter 11 of <cit.> or <cit.>. 

α
Let Π_V be a centred Gaussian probability measure on 𝒲 from (<ref>) with RKHS (,·_). Further suppose that for some α>d/2+1, the continuous embedding ↪ H^α+1(^d) holds. 

To deal with various non-linearities in our regression problem, we follow ideas in <cit.> and rescale the `base prior' via

    W = V/√(N)δ_N;        δ_N = N^-α+1+β/2(α+1)+2β+d,

to introduce extra regularisation in the posterior distribution. Here, α+1 models the regularity of the interaction potential W (so α models the regularity of ∇ W), while β >0 is determined by that of the initial condition ρ_W(0, ·) = ϕ∈ H^β. The extra regularisation allows us for instance to deduce space or time regularity estimates for ρ_W corresponding to interaction potentials W drawn from the posterior, see Theorem <ref> below. Note that such priors are special cases of the rescaled Gaussian process priors considered in several `direct' statistical settings <cit.>.




We will further consider L^2(^d)-projections π_𝒲_N(W) of the law of W onto sequences of linear subspaces 𝒲_N ⊆𝒲, including the case 𝒲_N = 𝒲 for all N. The final prior law of such W on 𝒲 or 𝒲_N will be denoted by Π=Π_N. A typical example of such a projection is the truncation of an infinite dimensional Gaussian process onto its finite-dimensional counterpart, as is usually done in practice. We note that less aggressive shrinkage would be permitted in (<ref>) but we opt for the present choice as this allows us to directly use results from <cit.> in the proofs.



The following theorem shows that we can solve the regression problem underlying (<ref>) via the implied prior arising from a Gaussian process model for the interaction potential W. [We restrict to even integers β to facilitate the application of parabolic PDE theory in the proofs.]


Suppose ϕ∈ H^β(^d) for some even integer β≥ 3+d and let Π_V be a Gaussian measure satisfying Condition <ref>. Consider the rescaled Gaussian process prior from (<ref>) projected onto a linear subspace 𝒲_N of 𝒲. Let W_0 ∈ℋ satisfy ∫_^d W_0 =0, and assume further that there exists a sequence W_0,N∈∩𝒲_N such that W_0,N_ = O(1) and W_0 - W_0,N_H^-β = O(δ_N) as N→∞. Then for M,L>0 large enough, we have

    Π(W ∈𝒲_M,N : ρ_W-ρ_W_0_L^2()≤ L δ_N|D_N) →^P^N_W_0 1

as N→∞, where

    _M,N = { W=W_1+W_2 ∈𝒲_N: W_1_H^-β≤ Mδ_N,  W_2_≤ M,  W_C^2≤ M }.



The convergence rate δ_N = N^-α+1+β/2(α+1)+2β+d derived above corresponds to the forward rate of an (α+1)-smooth truth in an inverse problem with `polynomial ill-posedness' β. The self-convolutional term ∇ W ∗ρ_W appearing in models (<ref>) and  (<ref>) suggests the ill-posedness here should be related to the smoothness of ρ_W_0. In the present nonlinear parabolic PDE setting, the regularity of the densities ρ_W_0 is driven by that of the initial condition ϕ rather than that of W_0, see Theorem <ref>, which is reflected in the role of β in the forward rate, see (<ref>) below. Whether this is optimal is a delicate question concerning the optimality of regularity estimates in nonlinear parabolic PDEs. Note that the convergence rate δ_N accelerates to the `parametric' rate 1/√(N) as soon as either the model for W or the initial condition ϕ is very smooth, i.e., as α or β→∞. 

To compare the convergence rates obtained in the preceding theorem to the existing literature, let us consider a general version of the interacting particle system (<ref>) on all of ^d:

    dX_t^i = - ∇ U(X^i_t)dt -1/n∑_i ≠ j∇ W(X_t^i-X_t^j)dt + √(2) dB_t^i,    i=1, …, n,

where U is a suitable confining potential. For (spatially) β-smooth ρ  solving the corresponding non-linear Fokker-Planck equation, Della Maestra and Hoffmann <cit.> show that the minimax rate for pointwise estimation of ρ is the `usual' nonparametric rate n^-β/2β+d. One can check that their minimax lower bound is established by considering W=0, that is, in a diffusion model without interaction term, where regularity of U translates directly into that of ρ. However, when U is zero or much smoother than W, our results show that we can expect significantly faster rates n^-α+1+β/2(α+1) +2β+d for the forward map ρ_W if the deconvolution structure is correctly taken advantage of.

This suggests that methods that are only based on the smoothness of ρ∈ H^β and not (as we do via our Lemma <ref>) on properties of the underlying parameter to solution map (U,W) ↦ρ_U,W may not be able to fully exploit the information available in the data. This highlights a significant advantage of modelling W directly, as the Bayes method does. We also note that by placing a prior on W, and hence ρ_W, one obtains posterior draws ρ_W that lie in the range of the forward map (i.e., they solve the PDE (<ref>)). This is key for applying our stability estimate (<ref>) below for recovery of W. While the Bayes method does this automatically, it is not necessarily the case for other methods that first estimate ρ_W directly by some smoothing method which relaxes the PDE constraint.

In summary, our results demonstrate that the estimation rates for ρ are more complex than the regularity of ρ alone might suggest, and that they depend on a delicate interplay of the interaction term W, the initial condition ϕ, and the confining potential U. 





  §.§.§ Examples of Gaussian priors


We now discuss examples of Gaussian priors satisfying the conditions of Theorem <ref>.

[Periodized Matérn process]Consider the Gaussian process on ^d given by

    V(x) =  (2π)^d/2∑_0 ≠ k∈^d1/(1+4π^2|k|^2)^(α+1)/2 g_k e_k(x),      g_k ∼^iid N(0,1),    x∈^d,

where the e_k ∝ e^2π i k (·) are the L^2(^d)-orthonormal trigonometric polynomials over ^d. This is the L^2(^d)-series expansion resulting from periodising the usual Matérn process on ^d, see Section A.1.1 of <cit.> for details. We have excluded the constant function from the trigonometric basis to incorporate our identifiability condition ∫_^d W=0.

By standard arguments (Ex. 2.6.15 in <cit.>), the RKHS of V equals 

    := H^α+1(^d) ∩{∫_^d W=0},

with equivalent norm ·_≃·_H^α+1. The above process has a version V whose sample paths are in C^α+1-d/2-η(^d) for any η>0 (<cit.>, Proposition I.4). Thus using <cit.>, Lemma I.7, we see that V defines a Gaussian random element of C^r_0(^d) for any r_0<α+1-d/2, and since for α>1+d/2 the space C^r_0(^d) imbeds continuously into the separable Besov space B^r_∞ 1(^d) ↪ C^2(^d), 2<r<r_0, (cf. p.370f. in <cit.>) we can realise the law of V as a Gaussian Borel probability measure on 
    𝒲 = B^r_∞ 1(^d) ∩{∫_^d W =0},  any  r<α + 1 -d/2,
 so that Condition <ref> is verified.

 Passing through a separable space is important to apply techniques from Gaussian process theory <cit.>, Ch.2, but the introduction of Besov spaces could be avoided by increasing the smoothness to α>1+d and taking 𝒲 = H^α-d/2 as in Theorem B.1.3 in <cit.>. 


Now assume W_0 ∈ defined in (<ref>) which also lies in 𝒲 (by what precedes, or abstractly by Corollary 2.6.17 in <cit.>). Applying Theorem <ref> with 𝒲_N =𝒲, W_0,N= W_0, then gives that for M,L>0 large enough,

    Π(W∈𝒲_N:  ρ_W-ρ_W_0_L^2()≤ LN^-α+1+β/2(α+1)+2β+d,  W_C^2≤ M |D_N) →^P^N_W_0 1,

as N→∞ for this prior.



[Truncated Fourier prior] For implementation by MCMC, and also in order to conduct inference on the interaction potential later, it is of interest to remove the higher Fourier frequencies from the Matérn prior for W just constructed. For K_N ∈ℕ, K_N →_N →∞∞, consider the truncated Fourier series prior

    π_𝒲_N(W)(x) =1/√(N)δ_N∑_k∈^d: 0<|k| ≤ K_N1/(1+|k|^2)^(α+1)/2 g_k e_k(x),      g_k ∼^iid N(0,1),

which equals the L^2-projection of W onto 𝒲_N = E_K_N, the span of the trigonometric polynomials up to frequency K_N (excluding constants). Let W_0 ∈, and set W_0,N = ∑_k∈^d: 0<|k| ≤ K_N⟨ W_0 , e_k⟩_L^2  e_k ∈∩𝒲_N to be the Fourier projection of W_0 onto the first K_N frequencies. Then W_0,N_≤W_0_H^α+1 and W_0-W_0,N_H^-β≲ K_N^-α-1-β =O(δ_N) for K_N ≳ N^1/2(α+1)+2β+d. Applying Theorem <ref> gives that for M,L>0 large enough,

    Π(W∈𝒲_M,N:  ρ_W-ρ_W_0_L^2()≤ LN^-α+1+β/2(α+1)+2β+d |D_N) →^P^N_W_0 1,

as N→∞, with regularisation 𝒲_M,N set given in (<ref>). Note the condition K_N ≳ N^1/2(α+1)+2β+d means the prior scaling rather than the truncation drives the prior smoothness.


 In many physical applications, the interaction potential W is symmetric. This can be encoded into the prior by separating the real and imaginary parts of the trigonometric polynomials e_k ∝ e^2π i k (·) into ϕ_2k(x) ∝cos (2π kx) and ϕ_2k+1(x) ∝sin (2π kx), and  setting the prior coefficients of ϕ_2k+1 to be zero. This yields prior and then also posterior draws that are symmetric about zero and all our results apply equally for Sobolev smooth symmetric truths W_0.



 One can also consider smoother Gaussian processes, for instance

    V(x) =  (2π)^d/2∑_0 ≠ k∈^d e^-r|k|_1/2 g_k e_k(x),      g_k ∼^iid N(0,1),    x∈^d,

where |k|_1 = |k_1|+…+|k_d| and r>0. These model infinitely differentiable periodic functions with exponentially decaying Fourier coefficients. For example, as in <cit.>, we consider the rescaled Gaussian process

    W = V/log N,

where less rescaling is needed than above since the sample paths of V are already very regular. The RKHS of V is

    _r = { h = ∑_0 ≠ k∈^d h_k e_k: h_^2:= ∑_0 ≠ k∈^d h_k^2 e^r|k|_1 < ∞},

which embeds continuously into any Sobolev space H^α+1, α>0, meaning we have all the required regularity properties needed above. Our arguments can then be modified to extend to such settings by standard arguments, see Section C of <cit.> for details. For W_0 ∈_r, we thus have for M,L>0 large enough,

    Π(W∈𝒲_N':  ρ_W-ρ_W_0_L^2()≤ Lδ_N',  W_W^2,∞≤ M |D_N) →^P^N_W_0 1

as N→∞, where δ_N' = (log N)^η/√(N) for some η>0 and _N' = {W=W_1+W_2: W_1_H^-β≤ M δ_N', W_2__r≤ M, W_C^2≤ M}. Such priors work for truths with exponentially decaying Fourier coefficients and are thus relevant when the interaction potential W_0 can be well approximated by relatively few Fourier modes. This occurs in several physical applications, for example the Kuramoto-Shinomoto-Sakaguchi (O(2)) model for synchronization or the (wrapped) Gaussian (attractive or repulsive) interaction potential; we refer to <cit.> for further discussion and several examples.






 §.§ Recovering the interaction potential


We now consider conditions under which we can consistently recover the interaction potential W_0 from noisy observations (<ref>) of the McKean-Vlasov PDE (<ref>). 



We notice that the posterior contraction rate δ_N for ρ_W - ρ_W_0_L^2 from the previous theorem becomes close to the parametric rate 1/√(N) as β→∞, even if we keep the regularity α of the prior model for the potential W fixed. This reflects the self-convolving nature of the McKean-Vlasov dynamics: the density ρ_W can be very smooth even when W is not, as long as the initial condition is sufficiently regular in a Sobolev sense. A consequence is that the interacting particle dynamics partially mask information about W unless the initial condition is somehow atypical. Viewed this way, inferring W becomes a statistical inverse problem which compounds a deconvolution step and a nonlinear inverse problem for the Fokker-Planck equation (<ref>). We now show that this inverse problem can be solved, provided that the initial condition ϕ is not too smooth, measured by a lower bound on the decay of its Fourier coefficients. In other words, if we can prepare the system to start in a somewhat irregular initial state, then the observed dynamics of the interacting particle system will reveal the potential W to the Bayesian posterior distribution arising from a prior that concentrates on at most a finite (but growing) number of Fourier modes. 

At the heart of the main statistical result Theorem <ref> is the following stability estimate for the forward map W ↦ρ_W under such hypotheses.



Suppose W_0∈ H^α+1 for α>1+d/2 and ϕ∈ H^β(^d) where β≥ 3+d is an even integer. Assume ϕ is strictly positive ϕ≥ϕ_min>0 on ^d and that the Fourier coefficients ϕ̂ of ϕ satisfy 

    |ϕ̂(k)| ≥ c |k|^-ζ,           k ∈^d,

for some ζ > β + d/2. Let W = ∑_k∈^d: |k| ≤ K W_k e_k ∈ E_K be band-limited by K Fourier frequencies and assume W_W^2,∞ + W_0_W^2,∞≤ M as well as ∫_^d (W-W_0)dx=0. Then there exists a constant C=C(d,T,M,ϕ_H^β, ϕ_min)>0 such that for any 0<t_0 ≤c/2C K^-ζ we have

    W - W_0_L^2^2 ≤ C (  K^-2α-2 + t_0^-1K^2ζ(ρ_W- ρ_W_0^2_H^1([0,t_0], H^-1) + ρ_W- ρ_W_0^2_L^2([0,t_0], H^1 )) ).




The analytical intuition of this stability estimate could be phrased as saying that if we can obtain precise estimates of the space and time derivatives of ρ_W_0, we will be able to identify W_0 ∈

under the hypothesis (<ref>). So both the space and time dynamics of the densities ρ_W are relevant for reconstruction. We can combine this estimate with the truncated Fourier series prior in Example <ref> to consistently recover the underlying potential. 


Suppose W_0 ∈

for α>1+d/2 and ϕ∈ H^β(^d), where β≥ 3+d is an even integer. Further assume that ϕ is lower bounded by ϕ_min>0 and satisfies (<ref>). Consider the truncated Fourier series prior in Example <ref> with K_N ≃ N^1/2(α+1)+2β+d. Then for L>0 large enough,

    Π(W∈ E_K_N: W-W_0_L^2≥ LK_N^3ζ/2δ_N^β-2/β  |D_N) →^P^N_W_0 0,

as N→∞, where δ_N = N^-α+1+β/2(α+1)+2β+d. In particular, the rate is N^-θ for

    θ = (α+1+β)(β-2)/β - 3ζ/2/2(α+1)+2β+d.

Moreover, if W̅_N = E^Π[W|D_N] ∈ E_K_N is the posterior mean, then we also have 
    W̅_N - W_0_L^2 = O_P_W_0^N(N^-θ)  as  N →∞.
 


Theorem <ref> uses the stability estimate (<ref>) with t_0 ≃ K_N^-ζ→ 0 small, which shows that as long as the initial condition ϕ is sufficiently irregular, the short-time dynamics of the system are already enough to consistently recover W as N→∞. This is especially relevant in our setting where ergodic averages may be uninformative for W, since the uniform distribution is always a stationary state for the mean-field dynamics on the torus <cit.>. Whether one can further exploit long-time information under relevant physical scenarios is a subtle question given the existence of phase transitions and, consequently, of multiple stationary states for McKean-Vlasov dynamics. In any case, our result already shows that Bayesian methods can extract relevant short-time information from the data.

The condition ζ >β+d/2 ensures that ϕ can satisfy both ϕ∈ H^β and (<ref>) simultaneously. If the lower bound (<ref>) is too weak, i.e. the initial density is much harder to deconvolve that its Sobolev smoothness β suggests, then our result may not imply consistency. Consistency with polynomial rates is possible as soon as θ>0, which holds for ζ satisfying

    β + d/2 < ζ < 2(α+1+β)(β-2)/3β.

This range is non-empty for α>0 large enough. For semi-parametric applications to Bernstein-von Mises theorems (e.g., as in Sec. 4.1.1 in <cit.>) it is important that these rates can be sufficiently fast in regular models: for instance we can always obtain 3 θ >1 as long as α≥α_0(β,d) for some α_0=α_0(β, d), if (<ref>) holds for ζ=β+d/2+η and all η>0.


One can further check that the rate K_N^3ζ/2δ_N^β-2/β holds in Theorem <ref> if one takes any K_N ≳ N^1/2(α+1)+2β+d, though taking larger K_N here deteriorates the final rate. One can also extend the theorem to K_N ≲ N^1/2(α+1)+2β+d giving deconvolution rate

    K_N^-α-1 + K_N^3ζ/2( ξ^β-2/β + ξ_N^β/β+1),

where now ξ_N = δ_N + K_N^-α-1-β is a possibly slower contraction rate (since the prior truncation may now drive the forward rate).

 We can also consider the truncated version of the prior in Remark <ref>, truncating at level K_N as in (<ref>). For W_0 ∈_r defined in (<ref>) and K_N ≥1/rlog N, we obtain the same near-parametric contraction rate δ_N' = (log N)^η /√(N) for the densities ρ_W as in Remark <ref>. Turning to the deconvolution, we deduce the rate

    W-W_0_L^2^2 ≲ e^-rK_N + K_N^3ζ (δ_N')^2(β-2)/β≲ (log N)^3ζ + 2(β-2)η/β N^-β-2/β.

Thus for W_0 ∈_r, one can also recover the true potential W_0, and thus solve the inverse problem, at the near-parametric rate (log N)^η'/√(N) as the regularity β→∞.




We next give a concrete example of a distribution ϕ for which the above theorem applies and consistent recovery of the potential W_0 with polynomial rates is possible.


[Periodised symmetric multivariate Laplace distribution] Consider the symmetric d-dimensional multivariate Laplace distribution on ^d, usually defined via its characteristic function 1/1+|ξ|^2/2, ξ∈^d. This distribution has a density function of the form (<cit.>, p. 235)

    f_d(y) = 2/(2π)^d/2( |y|^2/2)^2-d/2 K_(2-d)/2(√(2)|y|),     y ∈^d,

where K_λ is the modified Bessel function of the second kind, having integral representation K_λ(x) = ∫_0^∞ e^-x cosh tcosh (λ t) dt.
Consider the periodized version of this density

    f_d^per(x) = ∑_k∈^d f_d(x+k),      x ∈(0,1]^d,

which defines a strictly positive density function on the torus ^d with Fourier series

    f̂_d^per(k) = F_^d[f_d^per](k) = ∫_^d f_d^per(x) e^-2π i k.x dx =  F_^d[f_d](k) = 1/1+2π^2|k|^2,      k ∈^d,

where the relation with the Fourier transform of the unperiodized density f_d follows from Theorem 8.35 of <cit.>. For any m≥ 1, consider the m-fold convolution ϕ_m = f_d^per∗…∗ f_d^per which defines a positive probability density function on ^d with Fourier transform

    ϕ̂_m(k) = 1/(1+2π^2|k|^2)^m.

Then ϕ_m ∈ H^β for any β<2m-d/2 and satisfies the deconvolution condition (<ref>) with ζ = 2m.

Applying Theorem <ref> with W_0 ∈ℋ⊂ H^α+1∩ W^2,∞ and β = 2m-d/2-η for η>0 arbitrarily small, we obtain contraction rate 

    Π(W∈ E_K_N: W-W_0_L^2≥ LN^-θ'  |D_N) →^P_W_0 0

as N→∞, where θ' is any constant satisfying

    θ' < (α+1+2m-d/2)2m-d/2-2/2m-d/2 - 3m/2(α+1)+4m

(which can get arbitrarily close to the upper bound by taking η>0 small enough). The right-hand side is strictly positive if and only if

    α+1 > (2m-d/2)(m+d/2+2)/(2m-d/2-2),

in which case we can take θ'>0. Thus for any fixed m, one can consistently recover W_0 ∈ 

with polynomial rates for α large enough.





§ PROOFS




 §.§ Proof of Theorem <ref>: forward contraction rate


We verify the conditions of Theorem 2.2.2 in <cit.>, namely Condition 2.1.1 of <cit.> with regularization space = W^2,∞ on which the prior concentrates (in view of (<ref>) and the continuous imbedding C^2 ↪ W^2,∞). This consists of showing (i) boundedness of ρ_W on the domain  and (ii) Lipschitz continuity of the forward map W ↦ρ_W: H^-β(^d) → L^2(), uniformly over ·_W^2,∞-balls of radius M>0. While that theorem requires that W_0 ∈∩, it is straightforward to modify the proof to require only W_0 ∈ = W^2,∞ along with a suitably good approximation sequence W_0,N in the RKHS , see (2.22) or Exercise 2.4.3 in <cit.>.

(i) Boundedness: using the Sobolev embedding theorem H^d/2+(^d) ↪ C(^d) for any >0 and the L^∞([0,T] ; H^k+1(^d)) regularity estimate (<ref>) from Theorem <ref> with k=d-1,

    sup_W: W_≤ Mρ_W_L^∞()   ≤sup_W: W_W^2,∞≤ Msup_0≤ t ≤ Tρ_W(t)_H^d/2+
       ≤sup_W: W_W^2,∞≤ Mρ_W_L^∞([0,T] ; H^d(^d))≤ C_M,d,T,ϕ_0_H^d,

since ϕ∈ H^β with β≥ 3+d. This verifies (2.3) in Condition 2.1.1 of <cit.>.


(ii) Lipschitz continuity of the forward map: we control the distance between solutions to the McKean-Vlasov equation (<ref>) for different potentials via the relative entropy/Kullback-Leibler divergence of the marginal densities. We first require the following estimate from <cit.>, adapted on the torus.


Let W_i ∈ W^2,∞ (^d),   i=1,2, ϕ∈ H^3+d(^d) ∩(^d), and consider the solutions ρ_i ,   i=1,2, of the corresponding Fokker-Planck equations 

    ∂ρ_i/∂ t    =  Δρ_i + ∇· (ρ_i (∇ W_i ∗ρ_i) ),
              
    ρ_i(0,·)      =  ϕ.

Then for all t>0,

    H(ρ_1(t) | ρ_2(t))    ≡∫_^dρ_1(t,x) logρ_1(t,x)/ρ_2(t,x) dx  
       ≤1/2∫_0^t ∫_^d |∇ W_1 ∗ρ_1(u) - ∇ W_2 ∗ρ_2(u)|^2 ρ_1(u)   dx du.





The estimate follows from <cit.>, Lemma 3.1, which follows from <cit.>, Theorem 1.1. The proof can be easily adapted so that it holds for Fokker-Planck equations posed on ^d instead of ^d. We now check the conditions (H1)-(H3) of <cit.>[Lem. 3.1]. 

Under our assumptions on the interaction potential and on the initial condition ϕ, from Theorem <ref> below it follows that the Fokker-Planck equations (<ref>) have a unique, strictly positive solution ρ_i(t) ∈(^d) ∩ C^2(^d) for all t >0. 
By standard convolution inequalities, ∇ W_2 ∗ρ_2 is bounded, therefore (H1) is satisfied. This implies, in particular, that the weighted L^p spaces in the proof of <cit.>, Lemma 3.1 can be replaced by the flat spaces L^p(^d).  Furthermore, since W_i ∈ W^2,∞ (^d) and ρ_i is bounded on ^d for all T>0, together with Young's inequality for convolutions, we have that ∇ W_1 ∗ρ_1 ∈ L^2((0,T) ; L^p (^d )) for p > d+2 (in fact, we can take p=∞),  ∇ W_2 ∗ρ_2 ∈ L^2((0,T) ; L^2 (^d)) and ∇ W_1 ∗ρ_1 ∈ L^∞((0,T) ; L^2 (^d)), i.e. (H2)-(H3) are satisfied.





Let W_i ∈ W^2,∞ (^d), ϕ∈ H^3+d(^d) ∩(^d) and consider the corresponding solutions ρ_i = ρ_W_i to (<ref>). If W_1_W^2,∞≤ M, then for all t ∈ [0,T],

    ∫_0^t H(ρ_1(u) | ρ_2(u))   du ≤ C ∫_0^t∇ (W_1 - W_2) ∗ρ_2(t) ^2_L^2(^d),

where C = C(d,T,M,ϕ_H^d).




Arguing as in (<ref>) and using Theorem <ref>,

    ρ_1_L^∞()≤ C(d) ρ_1_L^∞([0,T];H^d(^d))≤ C_0(d,T,M,ϕ_H^d).

Using Lemma <ref>,  Young's convolution inequality and Pinsker's inequality (e.g., Prop. 6.1.7 in <cit.>) we have for all 0 ≤ t ≤ T:

    H(ρ_1(t) | ρ_2(t))     ≤   1/2ρ_1_L^∞()∫_0^t ∫_^d| ∇ W_1 ∗ρ_1(u) - ∇ W_2 ∗ρ_2(u) |^2    dx du 
       ≤    C_0 ∫_0^t ∫_^d| ∇ W_1 ∗ (ρ_1(u) - ρ_2(u)) |^2   dx du 
             + C_0 ∫_0^t ∫_^d| ∇ (W_1-W_2) ∗ρ_2(u) |^2   dx du 
       ≤    C_0 ∇ W_1 ^2_L^2(^d)∫_0^t ρ_1(u) - ρ_2(u) ^2_L^1(^d)  du  
             + C_0 ∇ (W_1  - W_2) ∗ρ_2 _L^2([0,t] ×^d)^2 
       ≤    2 C_0 W_1 ^2_H^1∫_0^t H(ρ_1(u) | ρ_2(u))   du + C_0 ∇ (W_1  - W_2) ∗ρ_2 _L^2([0,t] ×^d)^2.

Setting η (t) = ∫_0^t H(ρ_1(u) | ρ_2(u))   du, the above estimate can be written as

    η̇(t) ≤ 2 C_0 M^2 η(t) + C_0 ∇ (W_1  - W_2) ∗ρ_2 _L^2([0,t] ×^d)^2.

The result now follows using the differential form of Gronwall's inequality <cit.>.



The second part of the following lemma shows that the map W ↦ρ_W is β-smoothing, suggesting a mildly ill-posed inverse problem. We have not attempted to optimise β here but believe it close to being sharp. 



Let W_i ∈ W^2,∞ (^d), ϕ∈ H^3+d(^d) ∩(^d) and consider the corresponding solutions ρ_W_i to (<ref>). Suppose W_i_W^2,∞≤ M for i=1,2.


  (i) Then for some constant C = C(d,T,M,ϕ_H^d),

    ρ_W_1 - ρ_W_2_L^2(, λ)^2 ≤ C ∫_0^T ∇ (W_1-W_2) ∗ρ_W_2(t)_L^2(^d)^2 dt.


  (ii)
If in addition ϕ∈ H^β(^d), where β≥ 3+d is an even integer, then

    ρ_W_1 - ρ_W_2_L^2(, λ)≤ C W_1-W_2_H^-β

for some constant C(d,β,T,M,ϕ_H^β).




(i) For any densities p,q, recall that the Hellinger distance h(p,q) = ∫ (√(p)-√(q))^2 satisfies h^2(p,q) ≤ H(p|q)  (e.g. Lemma B.1 in <cit.>). Moreover, writing (p-q)^2 = (√(p)+√(q))^2(√(p)-√(q))^2 yields p-q_L^2^2 ≤ 4max (p_∞,q_∞)h^2(p,q). Using these facts and that ρ_W_i(t) are probability densities for all t≥ 0,

    1/T∫_0^T ρ_W_1(t,·) - ρ_W_2(t,·)_L^2(^d)^2 dt ≤4/Tmax_i=1,2ρ_W_i_L^∞()∫_0^T H(ρ_W_1(t,·)|ρ_W_2(t,·)) dt.

Using (<ref>) and Lemma <ref> then gives the result.


(ii) Let e_m(x) = e^2π i m.x, m ∈^d, denote the usual trigonometric polynomials. Since ⟨ f(x-·),e_m ⟩_L^2(^d) = - e_-m(x) ⟨ f,e_-m⟩_L^2(^d), Parseval's theorem yields

    | ∂_x_i (W_1-W_2) ∗ρ_W_2(t,·)(x) |^2 
    
       ≤ 4π^2 ( ∑_m∈^d |m_i| |⟨ W_1-W_2,e_m ⟩_L^2| |⟨ρ_W_2(t,·),e_m ⟩_L^2| )^2 
       ≤ 4π^2 ∑_m∈^d |m_i|^2 (1+|m|^2)^-1-β |⟨ W_1-W_2,e_m ⟩_L^2|^2 
           ×∑_m∈^d  (1+|m|^2)^1+β |⟨ρ_W_2(t,·),e_m ⟩_L^2|^2

for any x∈^d. This implies that

    |∇ (W_1-W_2) ∗ρ_W_2(t,·) (x)|^2 ≤ 4π^2 ρ_W_2(t,·)_H^1+β^2 W_1-W_2_H^-β^2.

By the regularity estimate (<ref>) in Theorem <ref>, ρ_W_2_L^2([0,T];H^β+1(^d))≤ C(d,β,T,M,ϕ_H^β)<∞. Integrating the second last display in space and time, and substituting this into the conclusion of (i) then implies

    ρ_W_1-ρ_W_2_L^2()^2    ≤ C ∫_0^T ρ_W_2(t,·)_H^β+1^2 dt  W_1-W_2_H^-β^2  ≤ C' W_1-W_2_H^-β^2.



The estimate given in Lemma <ref>(ii) is uniform over ·_=·_W^2,∞-balls for W, which verifies the Lipschitz estimate (2.4) in Condition 2.1.1 of <cit.> with κ = β. We may thus apply Theorem 2.2.2 of <cit.>, which completes the proof of Theorem <ref>.






 §.§ Proof of Theorems <ref> and <ref>: stability estimate and deconvolution



Let Ł_W,nψ :=Δψ  + ∇· (ψ∇ W ∗ρ_n-1), for smooth ψ, denote the operator defined in the linearized McKean-Vlasov equation (<ref>). We will prove the result for the linearized solutions before passing to the limit n→∞. Further define the linear parabolic operator

𝒫_W,n=(∂/∂ t - ℒ_W,n), which satisfies 𝒫_W,nρ_W,n=0 where ρ_W,n is a solution to (<ref>), see Theorem <ref>.



Step 1: Reduction to deconvolution with the particle density.
Let ρ_W, ρ_W_0 be two solutions of (<ref>) for distinct W,W_0, and denote by ρ_W,n,ρ_W_0,n their n-approximated linearised solutions. As the 𝒫_W,n operators are linear, we can write

    𝒫_W,n (ρ_W,n- ρ_W_0,n)    =0 - ∂ρ/∂ tρ_W_0,n + ℒ_W,nρ_W_0,n
       =  (ℒ_W,n-ℒ_W_0,n)ρ_W_0,n
       = ∇· (ρ_W_0,n∇ (W ∗ρ_W,n-1 - W_0 ∗ρ_W_0,n-1)).

The densities ρ_W_0,n(t,·) are strictly positive throughout [0,T] ×^d in view of (<ref>). Hence by standard theory for periodic elliptic equations, the operator L_0,n,tψ:=∇· (ρ_W_0,n(t) ∇ψ) is invertible on L^2_0(^d):=L^2(^d) ∩{∫_^d f =0} for every t: if we select the integral zero solution u = L_0,n,t^-1h to the elliptic equation L_0,n,t  u = h, we have a uniform in n, t elliptic regularity estimate (proved just as on p.127-128 in <cit.>),

    L_0,n,t^-1h_L^2≤ c h_H^-1,   h ∈ L^2_0(^d),

with c depending only on the (universal) constant in the Poincaré inequality on ^d as well as the lower bound ϕ_min from (<ref>). We also notice that for any t,

    ∫_^d𝒫_W,n (ρ_W,n- ρ_W_0,n)  = ∫_^d (Ł_W, n- Ł_W_0,n) ρ_W_0,n =0

as both Ł_W_n and Ł_W_0,n are in divergence form. Therefore we can apply the inverse of L_0,n,t to (<ref>) and obtain

    W ∗ρ_W,n-1 - W_0 ∗ρ_W_0,n-1 = L_0,n,t^-1[𝒫_W,n (ρ_W,n- ρ_W_0,n)],    on  [0,T] ×^d,

where we note that the left hand side integrates to zero since ∫ W = ∫ W_0, ∫ρ_W_0, n-1=∫ρ_W, n-1 =1. Rearranging we obtain the identity

    (W-W_0) ∗ρ_W_0,n-1 = L_0,n,t^-1[𝒫_W,n (ρ_W,n- ρ_W_0,n)] + W ∗ (ρ_W_0,n-1 - ρ_W,n-1),    on  [0,T] ×^d.

We can now fix any time horizon 0<t_0 ≤ T and estimate

    ∫_0^t_0(W-W_0) ∗ρ_W_0,n-1(t)^2_L^2(^d)dt 
       ≲∫_0^t_0L_0,n,t^-1[𝒫_W,n (ρ_W,n- ρ_W_0,n)(t)]_L^2(^d)^2dt + ∫_0^t_0W ∗ (ρ_W_0,n-1 - ρ_W,n-1)(t)_L^2(^d)^2 dt 
       ≲∫_0^t_0(∂/∂ t - Ł_W,n) (ρ_W,n- ρ_W_0,n)(t)_H^-1(^d)^2dt + ∫_0^t_0ρ_W_0,n-1 - ρ_W,n-1(t)_L^2(^d)^2 dt 
       ≲ρ_W,n- ρ_W_0,n^2_H^1([0,t_0], H^-1) + ρ_W,n- ρ_W_0,n^2_L^2([0,t_0], H^1).

We can take limits in this inequality: this is clear for the left hand side by Lemma <ref> and since convolution with W-W_0 is uniformly Lipschitz on L^2(^d). For the right hand side of (<ref>) we notice that by Theorem <ref> and Rellich's compactness theorem (p.305 in <cit.>), the sequences {ρ_W,n - ρ_W_0,n: n ∈ℕ} converge in H^1([0,t_0], H^-1) and L^2([0,t_0], H^1) along a subsequence to their unique limit, which in view of Lemma <ref> equals ρ_W - ρ_W_0. As the constants implied by ≲ in (<ref>) are n-independent, we deduce for any 0<t_0≤ T that

    ∫_0^t_0(W-W_0) ∗ρ_W_0(t)^2_L^2(^d)dt ≤ C( ρ_W- ρ_W_0^2_H^1([0,t_0], H^-1) + ρ_W- ρ_W_0^2_L^2([0,t_0], H^1))

for a constant C=C(T,d, M, ϕ_min, ϕ_H^β)>0.





Step 2. Deconvolving at small times.
Let W_0,K = ∑_k∈^d: |k| ≤ K⟨ W_0 , e_k⟩_L^2  e_k denote the Fourier projection of W_0∈ H^α+1 onto the first K frequencies. By Parseval's theorem and since ∫ W = ∫ W_0 we have

    W - W_0_L^2^2    ≤W_0_H^α+1^2 K^-2(α+1) + W-W_0,K^2_L^2
       ≲ K^-2α-2 + ∑_0<|k| ≤ K |Ŵ(k)-Ŵ_0(k)|^2 |ρ̂_W_0,n(t,k)|^2/|ρ̂_W_0,n(t, k)|^2
       ≲ K^-2α-2+sup_0<|k|≤ K1/|ρ̂_W_0,n(t, k)|^2(W-W_0) ∗ρ_W_0,n(t, ·)_L^2^2

For t>0, define the stability constant

    ι_K, n(t) ≡sup_0<|k|≤ K1/|ρ̂_W_0,n(t, k)|^2,

which we now upper bound. By Theorem <ref> we know 
    dρ_W_0,n/dt_L^∞([0,T];H^β)≤ C(d,T,β,M,ϕ_H^β),
 and hence t ↦ρ_W_0,n(t) is Lipschitz as a map from [0,T] → L^1(^d) which implies ρ_W_0,n(t) - ϕ_L^1(^d)≤ Ct. But since L^1-norms bound the Fourier coefficients uniformly, for any k∈^d,

    |ρ̂_W_0,n(k)| ≥ |ϕ̂(k)| - |ρ̂_W_0,n(k) - ϕ̂(k)| ≥ c|k|^-ζ - Ct.

Thus for any 0<t≤ t_0 = c2C K^-ζ, the stability constant in (<ref>) satisfies ι_K,n(t) ≤c^2/4 K^2ζ.

Integrating the inequality (<ref>) in time then gives

    ∫_0^t_0W - W_0_L^2^2 dt ≲ t_0 K^-2α-2 + K^2ζ∫_0^t_0(W-W_0) ∗ρ_W_0,n(t, ·)_L^2^2dt

and hence by (<ref>),

    W - W_0_L^2^2 ≲ K^-2α-2 + t_0^-1 K^2ζ(ρ_W- ρ_W_0^2_H^1([0,t_0], H^-1) + ρ_W- ρ_W_0^2_L^2([0,t_0], H^1)),

completing the proof.





The proof follows from combining the contraction rate for the forward map from Example <ref> with the stability estimate from Theorem <ref>. By Example <ref>, the posterior concentrates on the set

    𝒲_N = {W∈ E_K_N:  ρ_W-ρ_W_0_L^2() > Lδ_N ,  W_W^2,∞≤ M},

where δ_N = N^-α+1+β/2(α+1)+2β+d, so it suffices to prove the result for W∈𝒲_N. Using the interpolation inequality (A.5) in <cit.>, u_H^1(^d)≲u_L^2(^d)^1-1/γu_H^γ(^d)^1/γ for any γ>1. By (<ref>) of Theorem <ref>, ρ_W_L^2([0,T];H^β+1)≤ C for any W∈𝒲_N, and hence applying the interpolation with γ=β+1 and using Hölder's inequality yields

    ρ_W - ρ_W_0_L^2([0,T];H^1)^2    ≲( ∫_0^T ρ_W(t) - ρ_W_0(t)_L^2(^d)^2 dt )^γ-1/γ  = ρ_W - ρ_W_0_L^2()^2β/β+1.

Turning to the time regularity, using Proposition 2.1 in Ch. 4 in <cit.>, we have the interpolation inequality 
    v_H^1([0,T];L^2)≲v_L^2([0,T];L^2)^(m-1)/mv_H^m([0,T];L^2)^1/m

for any integer m≥ 1. But by (<ref>) of Theorem <ref>, d^ℓρ_W/dt^ℓ_L^2([0,T];H^1)≤ C for all ℓ =0,1,…,β/2 and hence we may take m=β/2 to obtain

    ρ_W - ρ_W_0_H^1([0,T];H^-1)≲ρ_W - ρ_W_0_L^2([0,T];L^2)^m-1/m = ρ_W - ρ_W_0_L^2()^β-2/β.

Apply Theorem <ref> with t_0 ≃ K_N^-ζ gives that for any W ∈𝒲_N ⊂ E_K_N,

    W-W_0_L^2(^d)^2    ≲ K_N^-2α-2 + K_N^3ζ( δ_N^2(β-2)/β + δ_N^2β/β+1) ≲  K_N^-2α-2 + K_N^3ζδ_N^2(β-2)/β.

Substituting in K_N ≃ N^1/2(α+1)+2β+d and δ_N shows that the second term dominates for all α>0 and β≥ 3+d since ζ > β + d/2. The convergence of the posterior mean vector now follows from arguments as in Theorem 2.3.2 in <cit.> and the details are left to the reader.






 §.§ Regularity properties for the McKean-Vlasov PDE


In this section, we collect properties and regularity estimates for the nonlinear parabolic McKean-Vlasov PDE (<ref>) driving the dynamics of our model, as well as its linearization. Recall that we work on the torus ^d = ^d/^d= (0,1]^d and that = (^d) denotes the set of probability densities (with respect to Lebesgue measure) on ^d. The basic existence and uniqueness of classical solutions to (<ref>) follows from <cit.>; see also <cit.>. 



If W ∈ W^2, ∞(^d) and ϕ∈ H^3+d(^d) ∩(^d), then there exists a classical solution ρ:[0,T]×^d→[0,∞) such that ρ (t, ·) ∈(^d) ∩ C^2(^d) for all t >0.


Similar to <cit.>, we approximate the McKean-Vlasov equation (<ref>) by a sequence of linear parabolic PDEs. For n=1,2,…, consider the sequence of linear equations on ^d

    ∂ρ_n/∂ t   =  Δρ_n + ∇· (ρ_n ∇ W ∗ρ_n-1)
    =: ℒ_W,nρ_n, 
    ρ(0,·)     = ϕ,

where ϕ∈ H^β(^d) ∩(^d) and we take a smooth time-independent initialization ρ_0 ∈ C^∞((^d)) (^d). One then shows as in <cit.>:

 
Let W ∈ W^2, ∞(^d) and ϕ∈ H^3+d(^d) ∩(^d). As k →∞ we have (if necessary along a subsequence)

    sup_0 ≤ t ≤ Tρ_W,k(t) - ρ_W(t)_H^1→ 0.
 



We now establish quantitative regularity estimates which hold uniformly over potentials W of bounded norm, which was crucial in our proofs above. It builds on the results of <cit.>. We remark that the smoothness of ρ_W is driven by the regularity of the initial condition ϕ, while the regularity of W is less crucial as it factors into the estimates only after convolution with ρ_n-1, which is a smoothing operation.


Let β≥ 3+d be an integer and suppose W ∈ W^2,∞(^d) satisfies W_W^2,∞≤ M. Let further ϕ∈ H^β(^d) ∩(^d) be a strictly positive probability density inf_xϕ(x) ≥ϕ_min>0 and let ρ_0 ∈ C^∞(^d) ∩(^d) be strictly positive. Then the sequence of PDEs (<ref>) has a unique solution ρ_n ∈ C^1,2([0,T] ×^d), which is strictly positive

    inf_x ∈^d, t ∈ [0,T]ρ_n(t,x) ≥λ (ϕ_min, M,T)>0,

for a constant λ that depends only on ϕ_min,M,T, and satisfies ∫_^dρ_n(t,x)   dx =1 for all t∈[0,T]. For k=0,1,…,β-1, we have the following regularity estimate in space:

    ρ_n _L^2([0,T] ; H^k+2(^d)) + ρ_n _L^∞([0,T] ; H^k+1(^d))≤ C(d,k,T,M,ϕ_H^k+1).

Furthermore, if β = 2j is an even integer, we have the time regularity

    d^ℓρ_n/dt^ℓ∈ L^2 ([0,T]; H^2j-2ℓ+1(^d)) ∩ L^∞ ([0,T] ; H^2j-2ℓ(^d))     0 ≤ℓ≤ j,

with the estimate

    ∑_ℓ=0^j( d^ℓρ_n/d t^ℓ_L^2([0,T] ; H^2j - 2ℓ +1(^d)) + d^ℓρ_n/d t^ℓ_L^∞([0,T] ; H^2j -2ℓ(^d))) 
               ≤ C(d,j,T,M,ϕ_H^2j).





Moreover, all the above estimates hold with ρ_n replaced by ρ, the solution to the McKean-Vlasov equation (<ref>).


The existence of a unique solution ρ_n to (<ref>) that is a probability density follows from standard linear parabolic PDE theory, e.g. <cit.>. The basic regularity in time and space follow from Theorems 27.2-27.3 (time) and Theorem 27.5 (space) of <cit.>, see also <cit.>. In particular, since ρ_0 is smooth, for n=1  we have a linear parabolic PDE (written in non-divergence form) with smooth coefficients ∇ W ∗ρ_0 ∈ C^∞ and Δ W ∗ρ_0 ∈ C^∞. Furthermore, the compatibility conditions are automatically satisfied since we have periodic boundary conditions. Therefore, from <cit.> it follows that ρ_1 ∈ L^2 ([0,T]; H^β+1(^d)) ∩ H^β +1/2 ([0,T] ; L^2(^d)). From this we deduce that we have the same space/time regularity as ρ_1 for ∇ W ∗ρ_n-1 and of Δ W ∗ρ_n-1, since W ∈ W^2, ∞. Applying again <cit.> we conclude the existence and uniqueness of solutions ρ_n ∈ L^2 ([0,T]; H^β+1(^d)) ∩ H^β +1/2 ([0,T] ; L^2(^d)) to (<ref>), which justifies taking any derivatives we need in what follows. It remains to show the quantitative regularity estimates, which we show hold uniformly in n and over potentials W_W^2,∞≤ M. 












We start with space regularity and follow ideas in <cit.>. We prove (<ref>) by induction on k, starting with k=0. We multiply (<ref>) by ρ_n, integrate in x, do the standard integration by parts, use Hölder and Young inequalities to deduce

    1/2d/d tρ_n ^2_L^2 + ∇ρ_n ^2_L^2   = ⟨ρ_n, d/dtρ_n - Δρ_n ⟩_L^2
        = - ⟨∇ρ_n, ρ_n ∇ W_n-1⟩_L^2
       ≤1/2∇ρ_n_L^2^2 + 1/2ρ_n_L^2^2∇ W_n-1_L^2^2,

where we have used the notation W_n-1 := W ∗ρ_n-1. Since ∇ W_n-1_L^2≤∇ W_L^2ρ_n-1_L^1≤ M by Young's convolution inequality, we deduce

    1/2d/d tρ_n ^2_L^2 + 1/2∇ρ_n ^2_L^2≤M^2/2ρ_n_L^2^2

for all t∈[0,T], where we recall ρ_n = ρ_n(t,·). Using the differential form of Gronwall's inequality (<cit.>, Section B.2), we have

    ρ_n(t,·) ^2_L^2≤ e^M^2Tρ(0,·) ^2_L^2 = e^M^2Tϕ_L^2^2

for all t ∈ [0,T]. Furthermore, integrating (<ref>) in time and using (<ref>),

    sup_t ∈ [0,T]ρ_n(t,·) ^2_L^2 + ∫_0^T ∇ρ_n(t,·) ^2_L^2  dt ≤ (M^2T e^M^2T+1) ϕ^2_L^2.

We proceed in similar manner with the higher regularity estimates. For the H^2-in-space estimate, we multiply (<ref>) by ∇ρ_n, integrate by parts, use Hölder's and Young's inequalities 

and the fact that ρ_n-1_L^1 =1 to deduce that

    1/2d/d t∇ρ_n _L^2^2 + Δρ_n ^2_L^2    =    ∫_^d∇·(ρ_n ∇ W_n-1) Δρ_n   dx 
    
       ≤   1/2Δρ_n^2_L^2 + 1/2∫_^d| ∇·( ρ_n ∇ W_n-1) |^2   dx
    
       ≤   1/2Δρ_n^2_L^2 +  ∇ρ_n·∇ W_n-1_L^2^2 +  ρ_n Δ W_n-1_L^2^2
    
       ≤   1/2Δρ_n^2_L^2 +  ∇ W_n-1_L^∞^2 ∇ρ_n^2_L^2 + Δ W_n-1_L^∞^2 ρ_n _L^2^2.

Using (<ref>)-(<ref>), we similarly obtain 
ddt∇ρ_n _L^2^2 + Δρ_n ^2_L^2≤ C(M,T)ϕ_L^2^2. Integrating in time and using (<ref>) and Poincaré's inequality on the torus,

    sup_t ∈ [0,T]∇ρ_n ^2_L^2 + c ∫_0^T ρ_n ^2_H^2  dt ≤ C(T,M)ϕ_H^1^2.

This establishes (<ref>) for k=0.


For the higher regularity estimates (<ref>) with integers 1≤ k ≤β-1, we proceed by induction on k as in <cit.>. Assume that (<ref>) holds for k-1. We differentiate (<ref>) by ∂_ν for some multi-index ν with |ν| = k, multiply by Δ∂_νρ_n, integrate by parts and use Young's inequality to obtain

    1/2d/d t∇∂_νρ_n ^2_L^2 +  Δ∂_νρ_n ^2_L^2   = ∫_^d|Δ∂_νρ_n ∂_α∇· (ρ_n ∇ W_n-1 )  |   dx 
       ≤1/2Δ∂_νρ_n^2_L^2 + C_d ∑_i=1^d  ρ_n ( ∂_x_i W ∗ρ_n-1)  ^2_H^k+1.









But by Lemma <ref>,

    ρ_n ( ∂_x_i W ∗ρ_n-1) ^2_H^k+1≤ C(d,k,M) ρ_n _H^k+1ρ_n-1_H^k.

Substituting this estimate into (<ref>), summing over all ν = |k|, integrating over [0,T] and using the induction hypothesis,

    sup_t ∈ [0,T]ρ_n(t) ^2_H^k+1 + ∫_0^T ρ_n(t) ^2_H^k+2dt    ≤ C(d,k,M) ∫_0^T ρ_n(t)_H^k+1^2 ρ_n-1(t)_H^k dt  + ϕ_H^k+1^2 
       ≤ C(d,k,M,T,ϕ_H^k+1, ρ_0_H^k+1),

which establishes (<ref>) for k as desired. Now since the ρ_n are uniformly bounded in the relevant Hilbert space norms, the Banach-Alaoglu theorem implies that they converge weakly in these Hilbert spaces to their limits along a subsequence. Combined with the inequality h_H≤lim inf_n h_n_H for any weakly convergent sequence h_n → h in a Hilbert space, we may take limits in the last inequalities, to give the following apriori estimate for the solution to the McKean-Vlasov PDE

    ρ_L^2([0,T] ; H^k+2(^d)) + ρ_L^∞([0,T] ; H^k+1(^d))≤ C(T,ρ_0, M).

We now consider regularity in time. For β =2 j, we have that ρ_n ∈ L^2 ([0,T]; H^2j+1(^d)) ∩ H^j+1/2 ([0,T] ; L^2(^d)) with ∂^j-1_t ρ_n ∈ C([0,T] ; H^1(^d)), see <cit.>.




To get the quantitative estimates, we follow the proof of <cit.>, proceeding by induction on the number of time-derivatives ℓ. For ℓ=0, this follows immediately from (<ref>). Suppose now we have the time regularity estimate (<ref>) for some integer 0 ≤ℓ < j. We differentiate (<ref>) ℓ-times in t to obtain, using the notation ρ^(ℓ):=d^ℓρ/d t^ℓ,

    ρ_n^(ℓ+1)   =   Δρ_n^(ℓ) + ∇·(ρ_n ∇ W ∗ρ_n-1)^(ℓ)
        =    Δρ_n^(ℓ) + ∑_k=0^ℓℓk∇·(ρ_n^(k)∇ W ∗ρ_n-1^(ℓ-k)).

Taking the H^2j - 2 ℓ -1(^d)-norm in space and using Lemma <ref> since 2j - 2ℓ≥ 2,

    ρ_n^(ℓ+1)^2_H^2j - 2 ℓ -1(^d)   ≤    2 ρ_n^(ℓ)^2_H^2j - 2 ℓ +1(^d) + C_ℓ∑_k=0^ℓρ_n^(k)∇ W ∗ρ_n-1^(ℓ-k)^2_H^2j - 2 ℓ(^d)
       ≤    
                     2 ρ_n^(ℓ)^2_H^2j - 2 ℓ +1(^d) + C_j,M∑_k=0^ℓρ_n^(k)^2_H^2j - 2 ℓ(^d)ρ_n-1^(ℓ-k)^2_H^2j - 2ℓ-1(^d),

where C depends only on j and W_W^2,∞(^d)≤ M. Integrating over time t ∈ [0,T] then gives

    ρ_n^(ℓ+1)^2_L^2(0,T ; H^2j - 2 ℓ -1(^d))   ≤    2 ρ_n^(ℓ)^2_L^2(0,T ; H^2j - 2 ℓ +1(^d))
           
                      + C_j,M∑_k=0^ℓρ_n^(k)^2_L^2(0,T ; H^2j - 2 ℓ(^d))ρ_n-1^(ℓ-k)^2_L^∞(0,T ; H^2j - 2 ℓ-1(^d)).

Since 0 ≤ k,ℓ-k ≤ℓ, we may now apply the inductive hypothesis (<ref>) to conclude that

    ρ_n^(ℓ+1)^2_L^2(0,T ; H^2j- 2ℓ -1(^d))≤  C(d,j,T,W_W^2,∞,ϕ_H^2j,ρ_0_H^2j),

which is the first part of the required estimate for ℓ+1. Similarly, taking the H^2j-2ℓ-2-norm in space and arguing as above yields

    ρ_n^(ℓ+1)^2_L^∞(0,T ; H^2j - 2 ℓ -2(^d))≤  C(d,j,T,W_W^2,∞,ϕ_H^2j,ρ_0_H^2j),

which proves the desired estimate for ℓ+1. This proves (<ref>) for ℓ = 0,…,j as required, based on the smoothness of the initial condition ϕ∈ H^2j. As before, we can then pass to the limit n →∞ to obtain the estimates for the mean field PDE.






To obtain the lower bound (<ref>) to the solution to the McKean-Vlasov PDE, we follow the proof of <cit.>[Lem. 5.1] that is based on the argument presented in <cit.>[Thm. 2]. Fix T'∈(0,T], let ρ_t(x)= ρ(t,x) denote the solution to the McKean-Vlasov PDE (<ref>),  and consider the unique strong solution of the linear (in the sense of McKean) SDE

    d Y_t = ∇ W ∗ρ_T'-t (Y_t)   dt + √(2)  dB_t,    Y_0 =x.

Using Itô's formula, the PDE (<ref>) and taking expectations gives for t∈[0,T'],

    ^x ρ_T'-t (Y_t) = ρ_T'(x) - ^x ∫_0^t ρ_T'-s (Y_s) Δ W ∗ρ_T'-s (Y_s)   ds,

where ^x denotes the expectation with respect to the law of the process Y_t starting at x. We use now the fact that ρ_u(x) ≥ 0,   u ∈ [0,T], together with the estimate Δ W ∗ρ_u _L^∞≤W_W^2,∞ρ_u_L^1≤ M to deduce that

    ^x ρ_T'-t (Y_t) ≤ρ_T'(x) + M ∫_0^t ^x ρ_T'-s (Y_s)    ds.

From the integral form of Gronwall's inequality, it follows that for all t∈[0,T'],

    ^x ρ_T'-t(Y_t) ≤ρ_T'(x) (1 + M t e^M t) ≤ρ_T'(x) ( 1+MTe^MT).

Setting now t=T' gives that for any T' ∈ (0,T],

    ρ_T'(x) ≥ϕ_min/1 + MT e^MT,

as required. The same argument also applies to the sequence of PDEs (<ref>). We omit the details.





Let W ∈ W^2,∞(^d), u∈ H^k and v ∈ H^k-1 for integer k≥ 2. Then

    f (∇ W *g) _H^k≤ C(k,d,W_W^2,∞) f _H^kg _H^k-1.




The proof follows from the estimates in <cit.>, pp. 65-66.




Acknowledgements. RN was supported by an ERC Advanced Grant (Horizon Europe UKRI G116786) as well as by EPSRC programme grant EP/V026259. GP is partially supported by an ERC-EPSRC Frontier Research Guarantee through Grant No. EP/X038645, ERC through Advanced Grant No. 247031.




 






plain




