Efficient Sampling of Stochastic Differential Equations with Positive Semi-Definite Models
Authors: Anant Raj
Umut Simsekli
Alessandro Rudi
Venue: NeurIPS 2023 poster (2023)
Keywords: Kernel Methods
Sampling
Fokker-Planck Equation
Fractional Fokker-Planck Equation
Stochastic Differential Equations
Partial Differential Equations
TL;DR: Sampling of Stochastic Differential Equations with PSD Models

Abstract: This paper deals with the problem of efficient sampling from a stochastic differential equation, given the drift function and the diffusion matrix. The proposed approach leverages a recent model for probabilities (Rudi and Ciliberto, 2021) (the positive semi-definite -- PSD model) from which it is possible to obtain independent and identically distributed (i.i.d.) samples at precision $\varepsilon$ with a cost that is $m^2 d \log(1/\varepsilon)$ where $m$ is the dimension of the model, $d$ the dimension of the space. The proposed approach consists in: first, computing the PSD model that satisfies the Fokker-Planck equation (or its fractional variant) associated with the SDE, up to error $\varepsilon$, and then sampling from the resulting PSD model. Assuming some regularity of the Fokker-Planck solution (i.e. $\beta$-times differentiability plus some geometric condition on its zeros) We obtain an algorithm that: (a) in the preparatory phase obtains a PSD model with L2 distance $\varepsilon$ from the solution of the equation, with a model of dimension $m = \varepsilon^{-(d+1)/(\beta-2s)} (\log(1/\varepsilon))^{d+1}$ where $1/2\leq s\leq1$ is the fractional power to the Laplacian, and total computational complexity of $O(m^{3.5} \log(1/\varepsilon))$ and then (b) for Fokker-Planck equation, it is able to produce i.i.d.\ samples with error $\varepsilon$ in Wasserstein-1 distance, with a cost that is $O(d \varepsilon^{-2(d+1)/\beta-2} \log(1/\varepsilon)^{2d+3})$ per sample. This means that, if the probability associated with the SDE is somewhat regular, i.e. $\beta \geq 4d+2$, then the algorithm requires $O(\varepsilon^{-0.88} \log(1/\varepsilon)^{4.5d})$ in the preparatory phase, and $O(\varepsilon^{-1/2}\log(1/\varepsilon)^{2d+2})$ for each sample. Our results suggest that as the true solution gets smoother, we can circumvent the curse of dimensionality without requiring any sort of convexity.

Efficient Sampling of Stochastic Differential Equations with Positive Semi-Definite Models
Anant Raj
Umut Şimşekli
Coordinated Science Laboraotry University of Illimois Trbana-Chhappen Inria, Ecole Normale Supérieure PSL Research University, Paris, France. anant. raj fiiriiiia
Inria, CNRS, Ecole Normale Supérieure PSL Research University, Paris, France.
umut. wettmmekkiiinniiii
Alessandro Rudi
Inria, Ecole Normale Supérieure
PSL Research University, Paris, France. alessandro rudi wwui@innia.....i
Abstract
This paper deals with the problem of efficient sampling from stochastic differen- tial equation, given the drift function and the diffusion matrix. The proposed ap- proach leverages recent model for probabilities [RC2]] (the positive semi-definite PSD model) from which is possible to obtain independent and identically dis- tributed (i.i.d.) samples at precision with cost that is m m2d log(/ where is the dimension of the model, the dimension of the space. The proposed ap- proach consists of: first, computing the PSD model that satisfies the Fokker-Plande equation (or its fractional variant) associated with the SDE, up to error s, and then sampling from the resulting PSD model. Assuming some regularity of the Fokker- Planck solution (i.e. ß-times differentiability plus some geometric condition on zeros) We obtain an algorithm that: (a) in the preparatory phase obtains PSD model with L2 distance from the solution of the equation, with model of di- mension m = E (y/ (log(1/e))
where 1/2 < s Σ the fractional power to the Laplacian, and total computational complexity of O( log(1 (1))) and then (b) for Pokker-Planck equation, it is able to produce i.i.d. samples with er- ror E in Wasserstein-1 distance, with a cost that is O(d (//d=-2(++1/// f log( (1/=)2d++)
per sample. This means that, if the probability associated with the SDE is somewhat regular, i.e. 3 2 4d - 2, then the algorithm requires O log(1/e) 4.5d in the preparatory phase, and O 1/2 1og(1/))2d+3 2d+2) for each sample. Our results suggest that as the true solution gets smoother, we can circumvent the curse of dimensionality without requiring any sort of convexity.
Introduction
High dimensional stochastic differential equations (SDE) and their associated partial differential equations (PDEs) arise often in various scientific and engineering applications such control systems, aviation, fluid dynamics, etc (LT00 HF08 Tro10]. Sampling from an SDE an active area of research. Generally, a direct approach like Langevin sampling [DMM19] 19 CDJB20 is commonly employed to sample from SDEs. However, an alternative approach involves approximating the solution of the associated PDE in form that enables easy sampling. In this paper, focus the latter approach.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

The exact solutions of most of the high- sigh-dimensional PDEs are not computable in closed form. Hence, for this reason, finding approximate numerical solutions high-dimensional SDEs and related PDEs has been an important area of research in numerical methods [Ame14] [GN13] T2d12 With the advent of new machine learning methods, the question naturally arises as to whether we can leverage the expressive power of these models to effectively capture the solution of high sigh-dimensional partial differential equations (PDEs). Recent research efforts have explored this direction by employing various machine learning models for different classes of PDEs [BHJK20 LLMD18 2KL+ For instance, works such as [CHS221 Ste21 BQL 22 have harnessed the capabilities of reproducing kernel Hilbert spaces (RKHS) to model PDE solutions. With the recent developments in neural network models and architecture, [LLMDI8 KLY21 KPRS22 CMW consider neural networks to model and upproximate the solution of high dimensional PDEs. Despite impressive algorithmic innovation, the theoretical understanding of these approaches is still limited, leaving a significant class of PDEs beyond the reach of existing machine learning-based -based modeling frameworks with inadequate theoretical insights into the obtained solutions.
In this paper, our goal is to sample from the stochastic differential equations (SDE) that are driven by a Brownian motion or an a-stable process. The time evolution of probability density functions (PDFs) of stochastic differential equations (SDE) that are driven by Brownian motion or a-stable process is given by Fokker-Planch Eguation (FPE) and the fractional Fokker-Plane
Equation (fractional FPE) respectively (see e.g. [DMM19 [Sm HMW21]. Therefore, our focus in this study lies in approximating the solution of the Fokker-Planck Equation (FPE), along with its fractional counterpart, as discussed in [UHK18], and subsequently sampling from the approximated solution. The modeling of solutions for FPEs (FPEs) and fractional FPEs poses significant challenges, encompassing the following key aspects: (i) Solutions of FPEs and fractional FPEs are probability densities that always non-negative and vanish at infinity, (ii) Modeling a probability density is hard because of the normalization property, and (iii) Fractional FPEs involve fractional Laplacian operator that is non-local in nature. In order for an approximate solution of FPEs or fractional FPEs to be useful for sampling from the corresponding stochastic differential equations (SDEs), it is necessary that sampling from the approximate solution of the corresponding partial differential equation (PDE) be straightforward. Unfortunately, comprehensive algorithmic approach that addresses the above- mentioned issues and includes a formal approximation analysis is currently lacking in the literature.
This paper adopts a positive semi-definite (PSD) model-based approach effectively tackle the problem at hand. The recently proposed PSD model [MFBR20] RMFB20 RC21] offers promising solutions to the challenges encountered in probability density modeling, making it an excellent choice for representing the solution of FPEs and their fractional counterparts. In a recent work by (RC21) it was demonstrated that employing a Gaussian kernel in the PSD model allows for the exact determination of the normalization constant and enables the characterization of a broad range of probability densities using such PSD models. Moreover, an algorithm for sampling from the PSD model is also proposed recently with statistical guarantees [MFBR22].
Driven by these insights, explore the utilization of PSD model for approximating the solution of FPEs and fractional FPEs. Consequently, we employ the sampling algorithm presented in MF[MFB222 to generate samples from the approximate solution. Our primary objective in this study is to obtain rigorous guarantees on the approximation error of PSD models when applied to the approximation of solutions of FPEs and fractional FPEs. Under regularity conditions on the solution of these PDEs, we make the following contributions.
Contributions:
We make the following contributions in this paper:
We show that the PSD-based representation of the solution of FPE has good approximation proper- ties under the regularity assumptions on the true solution of the PDE. In particular, we show that we can achieve an approximation error of O(e) while using m 0( (d=))(()))((8))()) (log ! 1)d++1/2) number of samples for PSD representation where d is the dimension of space variables and the true solution is 3 times differentiable (Theorem 1 and 2). In addition, we are able to provide Wasserstein guarantees for the sampling error of PSD model-driven solution when sampling from an SDE driven by a brownian motion (Corollary 1). An important contribution of our approach is that it relies on much weaker assumptions on the solution of associated PDE (33 times differentiable) compared to existing methods that employ dissipativity, log- log-Sobolev, and Poincaré like conditions
to provide sampling error guarantees (DMM19 DK Dall7 DM17 CEL 22
BCE GCE+)2).

---

We demonstrate that utilizing a positive semi-definite (PSD) based representation is natural and effective method to approximate the fractional Laplacian operator when applied to a probability density. In particular, we provide results that illustrate the effectiveness of PSDD-base representa- tions for Gaussian and shift-invariant kernels, which can be represented using Bochner's theorem (Theorem 3).
As a final contribution, we also establish an approximation guarantee for the solution of the fractional FPE under certain regularity assumptions the true solution. This extends our earlier result for the FPE and demonstrates the effectiveness of the PSD>-ased approach in approximating the solution of both FPE and fractional FPE We achieve an approximation error of O(e) while using m = 00s=((41))(C (d+1))((--2)) (log ! )dd+1)/2) number of samples for PSD representation where d is the dimension of space variables and the true solution is /3 times differentiable (Theorem 4 and 5). This result directly implies Wasserstein guarantee on the sampling error when sampling from an SDE driven by an q-stable levy process using the PSD model. To the best of our knowledge, our algorithm stands as the first of its kind, providing a provable guarantee for sampling from an SDE driven by an o-stable levy process as well as for approximating the solution of fractional FPEs under regularity assumption on the solution.
Our results clearly suggest that as the true solution gets smoother, we require less number of samples to approximate the true solution using PSD-based representation, hence circumventing the curse of dimensionality. All the proofs are provided in the appendix.
2
Related Work
Eva22 provides a good introduction on PDEs. GW79 1se09 methods solve PDE length.
Amele have discussed numerical
Modeling the solution of a PDE using machine learning approaches and in particular kernel methods is a well-established technique leading to optimal sample complexity and very efficient algorithms (see for example [BKO 96 SW06 Sch07l Fas07 PR20, CHOS2T (te2II). With the recent advancement of neural networks in machine learning, new approaches have emerged. A deep learning- learning-based approach that can handle general high-dimensional parabolic PDEs was introduced in (H)188). Two popular approaches that were proposed in this context are (i) Physics inspired neural networks (PINN) [RPK19] and (ii) Deep Ritz method 18) However, no statistical optimality was provided in the original papers. Recently, 21) studied the statistical limits of deep learning and Fourier basis techniques for solving elliptic PDEs. The sample complexity to obtain solution via a PINN and Deep Ritz model are shown be O of the solution. Subsequently, in a series of papers [DR222 LL22 LLW2] CLL2I, the theoretical properties of solving different class of PDEs are discussed under restrictive assumptions on the function class (barron space). For an overview, we refer to [BH)20 Solutions of PDEs that involve a fractional operator have also been modeled by function RKHS [Arq AA20 ASA 19 AAS20] However, these works do not provide finite sample approximating error guarantees for the solution obtained. review of kernel-based numerical methods for PDEs can be found in (SS 16 FFTS.
Machine learning-based solution to PDEs:
(n- (n~ J12 logn) and O (n
28-2
a=-2 log respectively where 3 is the order of smoothness 29-2
The (fractional) Fokker-Planck equation: Various efforts have been made to approximate the Yan
solution of FPE and its fractional counterpart by numerical methods [TSIG22 PMB13
LAT04 Bar01] [KRP12] DM14 SS22) utilized radial basis kernel based approach to numerically compute the solution of an FPE. RKHS- based methods have also been used to obtain an approximate solution BH17 SAA20 However, the above mentioned works are numerical in nature and do not discuss the approximation limits of the solution. On the more applied side, various works have been proposed by harnessing the expressive power of deep networks to model the solution of FPEs and fractional FPES fxZL 20, ZDL22 WWLZ22] However, a systematic theoretical study is missing from the literature.
Langevin Monte Carlo: An important line of research considers the probabilistic counterpart of FPE that is directly simulating the SDEs associated with the FPEs by dime-disrettizatinn Under convexity assumptions [ DM17 DM19 proved guarantees for obtaining a single (approximate) sample from the SDE for a given time. These frameworks were later extended in [RRTI7] XCZG18 EH211 so that the convexity assumptions were weakened to more general notions such as uniform dissipativity. [CCF15 introduced a more general SDE framework, which covers a larger class of PDEs. In the case of sampling SDEs with domain constraints [BEL15] BDMP17 HKRC18 SR20
3

---

proposed Langevin- based methods under convexity assumption, which were later weakened in
2L22 [5im17 NSR19 considered SDEs which corresponded the fractional FPE in dimension one and this approach was later extended in [HWW21 7722 An overview of these methods is provided [NF21] Contrary to this line of research, which typically requires convexity or dissipativity, our framework solely requires smoothness conditions.
Lam2
PSD models are an effective way to model non-negative functions and enjoy the nice properties of linear models and more generally RKHS. It was introduced in (MFBR20) and their effectiveness for modeling probability distributions was shown in [RC21]. It has been also effectively used in optimal transport estimation [MVB*21 VMR+21] |
PSD Models:
in finding global minima of non-convex function [MMB22)) in optimal control [BCRB22] In [MFBR22], an algorithm to draw samples from PSD model was proposed.
Notations and Background
Notation: Space input domain is denoted X and time input domain is denoted as T. Combined space X is denoted X. L (X)), (+} and (3) denote respectively the space of essentially bounded, absolutely integrable and square-integrable functions with respect to Lebesgue measure over X. W (Xt denotes the Sobolev space of functions whose weak derivatives up to order / are interalin on X. We denote by R the space vectors in Re with positive entries, Raxd the space of n x matrices, S = S+ (R") the space of positive semi- memindefinie n x n matrices. Given a vector 1) E Rdd, we denote diag(n) E ddd the diagonal matrix associated to
3.1 (Fractional) Fokker-Planch Equation
The FPE describes the time evaluation of particle density if the particles are moving in vector field and are also influenced by random noise like Brownian motion. Let us consider the following SDE driven by the standard Brownian motion (Wiener process) W,:
dX; = (Xx,t))t + odW,
where X₁ E Rd is random variable, p(X,, t) Fd is drift, Bgd and W, is a k-dimensional Wienner process. tx( denotes the ith element of vector p(X,,t). The diffusion tensor" is defined by D = 1000
The probability density function p(z, t) of the random variable X, solving the SDE given by the FPE, given as follows [RR96 Pav14]
Öp(z, t) 0t
Σ
0 (mila, Op(a, ()) + Dos dr,Os, p(x, t), such that p(0, a) = =po(2).
8²
(2)
dx.
It is known that when (xX,, f) is the gradient of potential objective function X - R, i.e.,
= Vf. where I satisfies certain smoothness and growth conditions, the stationary distribution of the stochastic process in equation (I) exists and is described by the so-called Gibbs distribution 14y
Despite of having nice theoretical properties, Brownian-driven SDEs are not suitable in certain models in engineering, physics and biology [STL%, A popular non-Gaussiaa SDE used in such applications is given follows:
dx, t)d, dIg,
where L denotes the (rotationally invariant) a-stable Lévy process in Rd and is defined as follows
for a E (0, 2] [ual55]
1. LG 0 almost surely;
2. For any to f1
tN. the increments L,; LG independent;
3. The difference I T. and L have the same distribution, with the characteristic function
E [e exp(- (t- [[u|(() for t > s and = V-I;
4. LÝ is continuous in probability.
In the SDE described in equation 3 we choose = 1 for the simplicity of the expressions. When
= 2, the process recovers the Brownian motion, i.c., L /2W, The fundamental difference
'oo simplify analysis, will consider diffusion tensor 1) independent from X; and in this work.

---

between L and W when a < 2 is that the increments of this process are heavy-tailed, with the following property; |||| L [|"] Do for all p 2 a and t - 0. This heavy-tailed structure makes the solution paths of (3) discontinuous: (Xt) can exhibit at most countable number of discontinuities. SDEs based on stabee processes have also received some attention in machine learning as well [SSDE20]
H [SSDD] 20). The resulting governing equation for stable-driven SDEs is similar to the traditional Fokker Planck equation except that the order a of the highest derivative fractional when a < 2. The following generalization of the FPE is known the literature as fractional FPE because of the presence of fractional Laplacian operator [UHK18]
Öt Öp(z,t) t) = - dx.
Ø (mi( ()(a,, f)) - (-A)"/2 t), such that 2) = po(2),
(4)
where (-A)* is the non-local Laplacian operator which is also known as the fractional Lapla- cian operator and is induced by the stable process. For < s < 1, the fractional Laplacian of order s, (-A)* can be defined on functions f : I - I as singular integral follows: f(a) = ca. Sgd ||=- f(x)=f(i) dy, where Cd,s = 4) [- A more useful representation of the
fractional Laplacian operator is given as a Fourier multiplier, F()-A) = ||w])>>[[]]( (w), where Fi) Sne ff (z))e= dr. We will utilize the Fourier multiplier- sittiilee-bbeed representation of the fractional Laplacian operator to prove our approximation result for the fractional FPE.
3.2 PSD Models for Probability Representation
Consider feature map representation Φ X - H. from input space X to a Hilbert space H, and linear operator M (30), then PSD models are represented as in [MFBR20)
[(z; M, の) <(x) MΦ(x).
(5)
It is clear that f<; M, d) is a non-negative function and hence PSD models offer general way to parameterize non-negative functions. We consider H. to be the RKHS corresponding to the kernel k such that (xx) k(x,,) )) In particular, we will consider the case where: Φ 6, Rd K,, is the feature map associated with the Gaussian kernel (2, y) = 6,(xx) (T = e p--n|i--p/ with n 0. Let us assume that M lies in the span of feature maps corresponding to data points, for some A E ST. Hence, a Gaussian PSD) model can be defined as, fa;; A, X,,) = Eij Asykin (2t, ( Given two base point matrices X E Raxdd and X' 12mxxd then X,x',, .n denotes the kernel matrix with entries (Kx,x'm) m))s xz)) where z; and z'j are the i-th and j-th rows of X and X' respectively. Gaussian PSD model has lot of interesting properties and the details can be found in [RC2T]. But, we would like to mention below an important aspect of PSD models which enable them to be used as an effective representation for probability density:
X = [I1, } then M- = Eis A150(x;))(zzz)
Sum Rule (Marginalization and Integration): The integral of a PSD model can be computed as,
f f(z; dz C2n Tr(AKx,x.9), where c, / k(0,z) dz.
Similarly, only one variable of a PSD model can also be integrated to result in the sum rule.
Proposition 1 (Marginalization in one variable (RC21)). Let X E Buxd, R**' AES, and n, of E IR+,, then the following integral is a PSD model,
J (xx,v: A, YJ,Y, (1) dx fly B,Y, with B = c2n K Kx
In a recent work, [MFBR22 gave an efficient algorithm to sample from PSD model, with the following result
Proposition 2 (Efficient sampling from PSD models MFBR22 Let E > O. Given a PSD model f(x; A, with A mm for Rd. There exists algorithm (presented in (MFBR22V) that samples i..dd points from probability is such that the I-Wasserstein distance between f(; A, and fe satisfies (ff. fe) S e. Moreover, the cost of the algorithm is 2mm2dlog(dd dlog(d/=)).
4 Approximation of the (Fractional) Fokker-Planck Equation
In this section, we investigate the approximating properties of PSD models for solving FPEs and fractional FPEs. Our focus is on the approximation guarantee achieved by fitting a single PSD model
5

---

in the joint space-time domain. Let's denote the solution of the PDE as p² (I, t), where we consider
the FPE in section 4. and the fractional FPE in section 4.2 We introduce mild assumptions on the domain and the solution probability density. Throughout the paper, we assume the spatial domain X
be (-R, R) and the time domain I to be (0, R). We have the following assumption from [RC21
on the solution p¹ of the PDE,
Assumption 1. Let 3 2, q E N. There exists J2, So E density p² 5 R satisfies, p² p=(,t)= -
(X),),, such that the
The assumption above is quite general and satisfied by a wide family of probabilities, as discussed in the Proposition of [RC21). We reiterate the result here below.
Proposition 3 (Proposition 5, [RC21)). The assumption above is satisfied by
(i) any (2, f) that is B-times differentiable and strictly positive on [-R, RJ x [0, R).
(ii) any exponential model p(x,t) =((,t,, such that E W. ( (
(iii) any mixture model of Gaussians or, more generally, of exponential models from (ii).
(iv) any that is 3 + 2-times differentiable on [-R, R, Rt) x [0, Rt with a finite set of zeros in X
and with positive definite Hessian each zero.
Moreover when p* is B-times differentiable over | R. R L (X)).
X [0, R). then it belongs to W (20 n
We will need to impose some extra mild assumptions on SJ for j E (1,... . to approximate the solution of a fractional Fokker- Planckk equation because of involvement of the fractional laplacian op- erator (non-local operator). We will cite these extra assumptions when describe the approximation results in section 4.2 Apart from that, we also have a boundedness assumption for the coefficients of the PDE in the domain consider.
Assumption 2. We have for all (x,t) E Xt (G(x,)) t) s Rp. (1,5) E (1, (1,...,d).
< RRp, and Dij S for all Oz:
It is also obvious to see that the PSD-based representation of probability density vanishes at infinity.
4.1 Approximating Solution of Fokker-Planck Equation
Let us consider the Fokker-Planck equation given in equation 2 and pz t) is the solution of Fokker-Planne equation in equation (2) which satisfies the assumption Hence, p" p"(x, satisfies the following,
Op (z,t) at
Σ
a (ui(z, (I,,,,) Dr.
Do dr,Ou, ().
0²
For now, we are ignoring the initial value as we only are interested in approximation error and we will have uniform approximation of the solution density p" and its derivatives up to the order 2 in the domain.
Given m base points pair (æ,, [ X for i 21. .m of space and time, let us consider the following approximation model for the solution of FPE,
7(x, () Σ A,5 (1x, 1), (ā,, ()) kx.(((z f), ())
(6)
1
where Æx,r is the joint kernel across space and time and A is a positive semi-definite matrix of dimension Rmxm kx.r can be defined as the product of two Gaussian kernels kx and by which defined space and time domain respectively i.e kx,(((x t), (æ', t')) = kx (2, (t,,, We represent ox ) = kx (x,, .) E Hxx, ΦT(()= = (t, .) Hy and Φx.7(a, () Φx(x)
) 8 dr(().. Here, we assume that xx and ky are both Gaussian kernels and hence the RKHS Hx and Hy correspond to the RKHS of the Gaussian kernel. For simplicity, we keep the kernel bandwidth
clear that the joint kernel  is also Gaussian with the kernel bandwidth parameter equal to 17. Given M Er, Aujxx, (2,, t. (2,,,,)* p(x, t) can also be represented with function f((x, t), A, X, 17) as, T)T M (x, p(z, t) = (x, T
parameter 77 the same for lx and ky.
ola t) Av(2,))
(7)
= ((.A))))
T= A.G.)..

---

where (xxx,) t) E RR, [+(x,t)]; = (x, Rz) (r(t, matrix of space and time for m pair of samples (I, t2) for i {1. m).
and X E mmxd+1 is the concatenated data
We can now easily optimize an associate loss function with the FPE to obtain the unknown PSD
matrix A. Let us assume that there exists a positive semi-definite matrix A for which the PSD model
p(x. t) approximates the true solution p" (a, t) well. In that case, our goal is to obtain an upper bound on the following error objectíve:
0(p' (I, 1)
Σ
Oz (x,(x))))' (x(x,))
d ΣΣ Do 1m1 dr.On,
0
(p' -a) t)
(8)
Ot
(2(())
Our approximation guarantee works in two steps. In the first step, we find an infinite dimen- sional positive operator Mg : H[ @ HGr + H(x @ HCr such that if we denote p(æ, t) = (x,r(x,t) (M,x,,rr, t) then for some Me, p is e-approximation of the true density p'. Clearly, for a set of functions fi, .o E Hx @ Hr, Mg = - is a positive operator from Hx @ Hy + Hxx & Hy. Now, we carefully construct f, for {1, [1,..., .q) in RKHS Hxx & Hr such that approximation guarantee in theorem Dis satisfied. Details of the construction is given in Appendix. We represent (r, t) as a d + dimensional vector and a mollifier function g (details are in the appendix) such that (z) = v (d+1) (ā/U). Then, we construct our f, = fy * 9u for j E (1,2.... .9}. We prove the result in Theorem I with this construction. To prove approx- imation error guarantee for the solution of FPE using PSD model, we assume the derivatives of fi, Sg up order 2 are bounded which is already covered in Assumption
Theorem 1. Let 3 > 2, q E N. Let ft,...,Ja E W2(R)) n (R()) and the function p" Sh. Let s E (0. and let 1) E R Let Φx. be the feature map of the Gaussian kernel with bandwidth n and let Hx @ Hr be the associated RKHS. Then there exists M2 E S+Hx @ HG) with ran S q, such that for the representation p(z,t) = 4x X..x, follloooo holds under assumption on the coesficients of the FPE,
0p(z, t)
8
2z Dis dr, p(a, 1)
(1, 1)p(xx)))
= =O(e).
(9)
Dx.
1=
12(3))
and tr(M,) S 2/2
exp (1+6#
20
( CC. (t))))))))
where |r] det(diag(n)) and Č" andC' depend only on 3, d. (Ray) (eyy
Af Sa; L3(20)
and
f2 0zjmk L3(24)
Even though, from theorem we know that there exists positive operator M Hxx 2 Hy. - Hx & Hy such that 6(x, t) Φx.r t) M. t) approximate the solution p* (r. of FPE from Theorem [ it is not clear how to compute M2 as we do not have access of fy for j 21,2...... ... q) beforehand. Hence, to get the final approximation error bound we need to further approximate p(z, !) with PSD model of finite dimension p(x, 1) which can be computed if we have the access to a finite number of base points (z,, t,) for {1, .m) in the domain X. Let us now construct matrix Am as follows. Consider Am K K (X.T)x,T) ZMZ K - (X,T)XXTT) where ((X.))XXTT) is a kernel matrix in Rmxm such that [K [K(x.Txx.)) = (x, 1.) (x,((1;;t)) tt). Let us define the operator Ž H[ @ HG + Rm such that Zu = (#(21, t1) u, (2 tm) u). and Z'a Σm (x t.)a; We further define projection operator P Z" K (x.mXx Z. Hence, we now define the approximation density p(z, t) as (x,,)) ((n,t)TAmnii... tl (xx :- ((.(a))).m.....)))))
(x, ( (M
-- ((.(a.0)PMMMPPP..
In the main result of this section below, we show that p(x,t) is a good approximation for p(z, f) and hence for the true solution p" (z, Theorem Let p" satisfy Assumptions I [ Let > and 3 > There exists a Gaussian PSD model of dimension E N. i.e., p(x,t) f((2, t), Am, X, 17) with Am E ST and X E mmxxd+1
7

---

and 1m E RQ such that with m = 0(= (d+11//1 (0-1))/((222)))) (log (d+11)/)) we have,
Op(z. t)
0 dz:
(1,
Dis dr, P(a, t)
= O(e). (2(0))
(10)
Theorem 2 indicates that there is a positive semi-definite matrix Am E R such that p(r, t) = (ar, 1)" Am (z, 1) is good approximation of p" (x, t). Next, we would see the usefulness of the above result in Theorem in sampling from the solution of FPE given the drift term μ satisfies the regularity condition in [BRS16 Theorem 1.1] and j W (X). In this case, the drift function j can be approximated by a function f in a Gaussian RKHS, and the approximation guarantees are due to interpolation result in Lemma The following result is the direct consequence of results in the Theorem 2
and Proposition 2
Efficient sampling method Let ja HG, be the approximation of μ such that || A L2(1) O(e),
by Gaussian kernel approximation. Let us apply the construction of Thm 2 j.
Step 1: Find the best model Am by minimizing the l.h.s. of equation (12) over the set of m- dimensional Gaussian PSD models. Note that in that case, since everything is a linear combination of Gaussians (or products of Gaussians), we can compute the L distance in closed form and minimize it over the positive semidefinite cone of m-dimensional matrices, obtaining A, exactly.
Step 2: Denote the resulting PSD model by (x, 1) = f(x, I; A,, 11). At any time / of interest we sample from the probability (.,t) via the sampling algorithm in [MFBR22).
Corollary (Error of sampling). Assume that μ E W. (X: satisfies the regularity condition as in [BRSTO Theorem 1.1]. For any time t, the algorithm described above provides i.i.d. samples from a probability Pe that satisfies in expectation E, W (p',, ( p.)² O(e2), where p == (;,t) the probability at time t of stochastic process with drift μ and diffusion D. Moreover, the cost of the algorithm is (mm² log(d/=)) where m O( 4/3 4/3- (log(()))
This result is proven in the appendix after theorem 2 It turns out that the parameter of the PSD model can be learned by simple semi-definite program [RMFB20] which has computational complexity of O(m 3.5 log(1/=)) where m is the dimension of the model.
Remark. It is evident from the aforementioned result that our method excels in sampling from the Stochastic Differential Equation (SDE) while ensuring the particle density adheres to Assumption 7 within a bounded domain. In contrast, traditional sampling methods such as Langevin struggle to sample from this extensive class due to the absence of dissipativity, tog-sobolev, and Poincaré-like conditions.
4.2 Approximating Solution of Fractional Fokker-Planck Equation
In this section, our aim is to approximate the solution pxx () of fractional Fokker-Planck equation, similar to the previous section. The equation replaces the Laplacian operator with a fractional Laplacian operator while keeping other terms unchanged. The fractional Laplacian operator, being non-local, is generally challenging to compute analytically. Therefore, approximating the fractional operator function has become separate research direction. Below, demonstrate that utilizing a PSD model-based representation is a natural choice for approximating the fractional operator acting on probability density. We present results for two cases: (i) Gaussian kernel with the bandwidth parameter and (ii) a shift-invariant kernel represented using Bochner's 's theorem [Rud17].
Theorem 3. Consider probubility density represented by PSD model as p(x) f(x; A, X,n) Σnj= Aigk(1,, a))((z,,)) for some PSD matrix A, where X E modd is the data matrix whose i-th row represents sample x, Bd for forie, e 1. if the kernel is Gaussian with bandwidth parameter n, we have the
f (1...., m,), then
(-2) C E Anll:11 +(l,=11+11e
csTT=]
i,5=1
where μ.- = 2jn(x, + x.),) Σ == I and C can be computed in the closed form.
2. if the kernel k have the bochner's theorem representation, i.e. k(æ - 1) = Sea g(w)e wn T (x-w) di where q a probability densiry, then
(-a) Σ Aij Elllue wel so (acta) gut (1)-a)),
11
where expectations are over We, ~ g(.).
8

---

We utilize the Fourier multiplier- guuteplier-based definition of the fractional Laplacian operator in Theorem 3 both cases discussed in Theorem B empirical estimation of (-A)"p(z), obtained from finite samples (empirical average), is sufficient for practical purposes. Shifting our focus from the representation of a non-local operator acting on probability density, we now turn to the approximation result for the solution of a fractional Biokker-Planch equation in this section. We assume that the optimal solution p'(x, t) can be expressed as a sum of squares of q functions, as stated in assumption I To establish the approximation results for the fractional ookker-Plancls equation, we need to impose additional mild assumption on fj for j E 1,... ,q in addition to assumption [
In
which we state next.
Assumption& 3. Let F())()) denotes the Fourier transform of a function f. then F(S(() E (20) (ix),, and 3((-A)" E L' (2t for i E (1....q..
Similar to the previous section, it requires two steps to obtain the approximation error bound. By similar construction as in section [4.1] we define an infinite dimensional positive operator M2-"
Hx @ Hr + Hxx @ He and based on this positive operator, we have probability density p as p(x, t) = 6x, t) Mc (z, t) which is (s) approximation of (z, t). The result is stated below in Theorem and proven Appendix
Theorem 4. Let B > 2.9 E N. Let fi...., satisfy assumptions T and 3 and the function p'= E'= f.. Let (0. 1] and l E B4 Let Φx.T be the feature map of the Gaussian kernel with bandwidth y and let Hx Hor be the associated RKHS. Then there exists Me E (Xxx @o HGr with ran S q, such that for the representation p(z, t) = & Mx following holds under
assumption 2 on the coefficients of the fractional FPE.
0p(z,)) 1)
0 (M6 ()(a,, f)) + (-))
Σ
= -(s).
(11)
Ot
Ooz
(2(30)
and tr(M.) S 0/m|2/2 (1 (1+ exp
C" g )ーक
where |n| (ct(ggig())), and ö and C depend only B, d and properties of fi,'.. fq.
Major difficulty in proving the result in Theorem Marises from the challenge of effectively managing the fractional Laplacian term. To tackle this issuc, we utilize the Fourier- based representation to gain control over the error associated with the fractional Laplacian. The details can be found in Lemma 7
presented in Appendix C.T] As a next step, we would follow similar procedure as in section 4. obtain the final approximation bound. In the second part of the proof, we utilize Gagliardo-Nirenber inequality [MP18) for fractional Laplacian to obtain the final approximation guarantee. In the next result, we show that for a particular choice of a PSD matrix Am. p(x,,) = (1x t) ( t) =
(x, 1)T PM PΦx,T 1) is a good approximation for the true solution p (1,, 1) of the fractional Fokker-Planck equation.
Theorem S. Let p¹ satisfy Assumptions 7 and Let > 0 and B > 2. There exists Gaussian PSD model of dimension m E N. i.e., t(z,t) f((x, t), Am. X, OX,T), with Am E St and
X E mmxd and Tm E R4
+. such that with m = (d+11/)))/(8 20) (log 2 (d++1)))) we have
0p(x,t) t) 0t
Σ
0
(mi(a, ()((,, ()) (-)) ((2,
0(E)
(12)
Des
(2(0))
Similar results as in Corollary [ can be obtained for the case fractional Fokker-Planck equations as well. We can directly utilize results from [RZGŞ23 bound the Wasserstein metric between the solution of two fractional FPE whose drift terms are close in 0 metric. However, to find a PSD matrix A, from finite samples does require two separate sampling procedures which would contribute to the estimation error, (i) sampling to estimate the mean in the estimation of the fractional Laplacian and (ii) sampling of a finite number of data points from the data generating distribution.
Remark. As far as we know, Theorem 5 represents the first proof of approximation error bounds for the solution of the Fractional Fokker- Planck Equation when the approximated solution is density function. This significant result allows for the utilization of the algorithu presented in MMFR222| to sample from the approximate solution.
9

---

Conclusion
In this paper, we study the approximation properties of PSD models in modeling the solutions of FPEs and fractional FPEs. For the FPE, we show that a PSD model of size m Ö( (ff=-(441/(322) can
approximate the true solution density up to order in L metric where /3 is the order of smoothness of the true solution. Furthermore, we extend our result to the fractional FPEs to show that the required model size to achieve E-approximation error to the solution density is m = Ö(s (d+1)dd11//--000 (8)))) where is the fractional order of the fractional Laplacian. In the process, we also show that PSD model-based representations for probability densities allow an easy way to approximate the fractional Laplacian operator (non-local operator) acting on probability density. As future research direction, we would like to investigate and obtain a finite sample bound on the estimation error under regularity conditions the drift term μ.
Acknowledgement
Anant Raj is supported by the a Marie Sklodowska-Curie Fellowship (project NN-OVEROPT 101030817) 17). Umut Şimşekli's research is supported by the French government under management of Agence Nationale de la Recherche as part of the "Investissements a'vvenir" program, reference MAR-19-P3IA-0001 (PRAIRIE 3IA Institute) and the European Research Council Starting Grant DYNASTY - 101039676. Alessandro Rudi acknowleges support of the French government under management of Agence Nationale de la Recherche as part of the "Investissements d' d'avenir" program, reference ANR--19-P3 19--333-00-0001 (PRAIRIE 31A Institute)and support of the European Research Council
(grant REAL 947908).
References
[AA20] Omar Abu Arqub. Numerical simulation of fluid flows via reproducing kernel method. International Journal of Numerical Methods for Heat & Fluid Flow, 30(1 1):4711-44733 2020.
partial differential equations arising in
[AAS20] Omar Abu Arqub and Mohammed Al-Smadi. An adaptive numerical approach for the solutions of fractional advection-diifuusui diffussin and dispersion equations in singular case under riesz's 's derivative operator. Physica A: Statistical Mechanics and its Applications, 540: 140:123257, 2020.
[Ame14] William F Ames. Numerical methods for partial differential equations. Academic press, 2014.
[Arq18] Omar Abu Arqub. Numerical solutions for the robin line-fractional partial differential equations of heat and fluid flows based on the reproducing kernel algorithm. International Journal of Numerical Methods for Heat & Fluid Flow, 28(4):828-88 856, 2018
[ASA19] Mohammed Al-Smadi and Omar Ahu Arqub. Computational algorithm for solving fredholm fractfraciiona partial integrodifferential equations of dirichlet functions type with error estimates. Applied Mathematics and Computation, 342:280-294, 2019.
[Bar01] E Barkai. Fractional fokker-planck equation, solution, and application. Physical Review E. 63(4):046118, 2001.
[BCE+22] Krishna Balasubramanian, Sinho Chewi, Murat A Erdogdu, Adil Salim, and Shunshi Zhang. Towards theory of non-log-concave sampling: first-order stationarity guarantees for langevin monte carlo. In Conference on Learning Theory, pages 2896-2923. PMLR, 2022.
[BCRB22] Eloise Berthier. Justin Carpentier, Alessandro Rudi, and Francis Bach. dnfinite-dimensin sums- of-squares for optimal control. 2022 IEEE 61 Conference on Decision and Control (CDC), pages 577-582. JEEE, 2022.
[BDMP17] Nicolas Brosse, Alain Durmus, fric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal langevin monte carlo. In Conference learning theory, pages 319 342. PMLR, 2017.
[BEL15] Sebastien Bubeck, Ronen Eldan, and Joseph Lehec. Finite-time analysis of projected langevin monte carlo. Advances in Neural Information Processing Systems, 28, 2015.
[BH17] Jake Bouvrie and Boumediene Hamzi. Kernel methods for the approximation of some key quantities of nonlinear systems. Journal of Computational Dynamics, 4(1&2):1-19 2017.
10

---

[BIUK20] Christian Beck, Martin Hutzenthaler. Armulf Jentzen, and Benno Kuckuck. An overview on deep learning-ban based approximation methods for partial differential equations. arXiv preprint UrXiv:2012.12344 2020.
[BKO* 96] Ted Belytschko, Yury Krongauz, Daniel Organ, Mark Fleming, and Petr Krysl. Meshless methods: overview and recent developments. Computer methods in applied mechanics and engineering. 139(1-4):3 1996.
[BQ Kaijun Bao, Xu Qian, Ziyuan Liu, Haiteng Wang, and Songhe Song. Function- valued rkhs baseb operator learning for differential equations, arXiv preprint arXiv:2202.09458 2022.
[BRS16] Vladimir Bogachev, Michael Röckner, and Stanislav V Shaposhnikov. Distances between transition probabilities of diffusions and applications nonlinear fokker-planck-kolmogony equations. Journal of Functional Analysis, 271(5):1262 1300, 2016.
[CDJB20] Niladri Chatterji, Jelena Diakonikolas, Michael I Jordan, and Peter Bartlett. Langevin monte carlo without smoothness. In International Conference on Artificial Intelligence and Statistics, pages 1716-1726. PMLR, 2020.
[CEL*21] Sinho Chewi, Murat A Erdogdu, Mufan Bill Li, Ruoqi Shen, and Matthew Zhang Analysis of langevin monte carlo from poincar "e log-soholev. arXiv preprint 12. 2021.
[CHOS21] Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Solving and learning nonlinear pdes with gaussian processes. Journal of Computational Physics, 447:1 10668, 2021.
[CLL21] Ziang Chen, Jianfeng Lu, and Yulong Lu. On the representation of solutions to clliptic pdes in barron spaces, Advances in neural information processing systems, 34:6454-644 6465, 2021.
[CMW*21] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physics- informed neural networks (pinns) for fluid mechanics: A review. Acta Mechanica Sinica, 37(12)11277-1738 2021.
[Dal17] Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Sociery. Series B (Statistical Methodology), pages 651-676,2011 2017.
[DK19] Amak S Dalalyan and Avetik Karagulyan. User- friendly guarantees for the langevin monte carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):5278 5311, 2019.
[DM14] Mehdi Dehghan and Vahid Mohammadi. The numerical solution of fokker-planck equation with radial basis functions (rbfs) based on the meshless technique of kansas approach and galerkin method. Engineering Analysis with Boundary Elements, 47:38-66 63, 2014.
[DM17] Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted langevin algorithm. 2017.
[DM19] Alain Durmus and Eric Moulines. High-dimensional bayesian inference via the unadjusted langevin algorithm. 2019.
[DM119] Alain Durmus, Szymon Majewski, and Blakej Miasojedow. Analysis of langevin monte carlo via convex optimization. The Journal of Machine Learning Research, 20(1):266--211, 2019.
[DRM22] Tim De Ryck and Siddhartha Mishra. Generic bounds the approximation error for physics- informed (and) operator learning. arXiv preprint arXiv:2205. 1193,, 2022.
[Dua15] Jinqiao Duan. An introduction to stochastic dynamics, volume 51, Cambridge University Press, 2015.
[EH21] Murat A Erdogdu and Rasa Hosseinzadeh. On the convergence of langevin monte carlo: The interplay between tail growth and smoothness. In Conference Learning Theory, pages 1776- 1822. PMLR, 2021.
[Eva22] Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Society, 2022.
[Fas07] Gregory B Fasshauer. Meshfree approximation methods with MATLAB, volume 6. World Scientifie, 2007.
[FF15] Bengt Fornberg and Natasha Flyer. Solving pdes with radial basis functions, Acta Numerica, 24:215 258, 2015.

---

[FGNI3] Xiaobing, Feng, Roland Glowinski, and Michael Neilan. Recent developments in numerical
methods for fully nonlinear second order partial differential equations. siam REVIEW, 55(2):205- 267, 2013.
[GW79] Ian Gladwell and R Wait. A survey of numerical methodes for partial differential equations. Oxford: Clarendon Press, 1979.
[HFO8] Isom Herron and Michael Foster. Partial differential equations in fluid dynamics. Cambridge Univ. Press, 2008.
[HJE18] Jiequn Han, Amulf Jentzen, and Weinan E. Solving high- digen-immnniinal partial differential equations using deep learning. Proceedings the National Academy of Sciences, 115(34):8505 85100 2018.
[HKRC18] Ya-Ping Hsich, Ali Kavis, Paul Rolland, and Volkan Cevher. Mirrored langevin dynamics. Advances in Neural Information Processing Systems, 31, 2018.
[HMW21] Lu-Jing Huang, Mateusz B Majka, and Jian Wang. Approximation of heavy-tailed distributions via stable-driven driven sdes. 2021.
[Isc09] Arich Iserles. A first course in the numerical analysis of differential equations. Cambridge university press, 2009.
Number 44,
[KK]+21] George Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang. and Liu Yang. Physics-infornet machine learning. Nature Reviews Physics, 3(6):222 2402224022222 2021.
[KLY21] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying. Solving parametric pde problems with artificial neural networks. Fumpean Journal of Applied Mathematics, 32(3):421-435 2021.
[KPRS22] Gitta Kutyniok, Philipp Petersen, Mones Raslan, and Reinhold Schneider. A theoretical analysis of deep neural networks and parametric pdes. Constructive Approximation, 55(1):73-125, 2022.
[KRP12] S Kazem, JA Rad, and K. Parand. Radial basis functions methods for solving fokker-planck equation. Engineering Analysis with Boundary Elements, 36(2):181-189 2012.
[Lam21] Andrew Lamperski. Projected stochastic gradient langevin algorithms for constrained sampling and non-convex learning. In Conference on Learning Theory, pages 2891-2937. PMLR, 2021.
[LAT04] Fawang Liu, Vo Anh, and Ian Turner. Numerical solution of the space fractional fokker-planck equation. Journal of Computational and Applied Mathematics, 166(1):209 219, 2004.
[LCL+21] Yiping Lu, Haoxuan Chen, Jianfeng Lu, Lexing Ying, and Jose Blanchet. Machine learning for elliptic pdes: fast rate gencralization bound, neural scaling law and minimax optimality. arXiv preprint arXiv: 068067 2021.
[LL22] Jianfeng Lu and Yulong Lu. A priori generalization error analysis of two-layer neural networks for solving high dimensional schrödinger eigenvalue problems. Communications of the American Mathematical Society, 2(01):1 1211, 2022.
[LLMD18] Zichao Long, Yiping La, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In International conference on machine learning, pages 3208-3216. PMLR, 2018
[LLW21] Jianfeng Lu, Yulong Lu, and Min Wang. A priori generalization analysis of the deep ritz method for solving high dimensional elliptic equations. arXiv preprint arXiv:2101.01708, 2021.
[LT00] Irena Lasiecka and Roberto Triggiani. Control theory for partial differential equations: Volume
1, Abstract parabolic systems: Continuous and approximation theories, volume 1. Cambridge University Press, 2000.
[MCF15] Yi-An Ma, Tianqi Chen, and Emily Fox. A complete recipe for stochastic gradient meme. Advances neural information processing systems, 28, 2015.
[MFBR20] Ulysse Martcau-Ferey, Francis Bach, and Alessandro Rudi. Nar-paraeetric models for non- negative functions. Advances in neural information processing systems, 33:12816 11211281612280 11826,, 2020.
[MFBR22] Ulysse Marteau-Feree Fercy, Francis Bach, and Alessandro Rudi. Sampling from arbitrary functions via psd models. In International Conference on Artificial Intelligence and Statistics, pages PMLR, 2022.
[MP18] Carlo Morosi and Livio Pizzocchero. On the constants for some fractional gagliardo gagliardo-nirenbery and sobolev inequalities. Expositiones Mathematicae, 36(1):32-77, -7, 2018.
12

---

[MVB*21] Boris Muzellec, Adrien Vacher, Francis Bach, François-Xave Xaviir Vialard, and Alessandro Rudi. Near-optimal estimation of smooth transport maps with kernel sums-of-squares. arXiv preprint arXiv:2112.01907, 2021.
[NF21] Christopher Nemeth and Paul Fearnhead. Stochastic gradient markov chain monte carlo. Journal of the American Statistical Association, 1 1665333)433-450, 2021.
[]irl Louis Nirenberg. On elliptic partial diffcrential equations. I| principio di minimo e sue applicationi alle equazioni funzionali, pages 1-48. Springer, 2011.
[9SGR19] Thanh Huy Nguyen, Umut Simsekli, Mert Gurbuzbalaban, and Gaël Richard. First exit time analysis of stochastic gradient descent under heavy-tailed gradient noise. Advances in neural information processing systems, 32, 2019.
[NSR19] Than Huy Nguyen, Umut Simsekli, and Gaël Richard. Non-asymptotic analysis of fractional langevin monte carlo for non-convex optimization. In International Conference on Machine Learning, pages 4810-4819. PMLR, 2019.
[Pav14] Grigorios A Pavliotis. Stochastic processes and applications: diffusion processes, the Fokker- Planck and Langevin equations, volume 60. Springer, 2014.
[PMB13] Lukas Pichler, Arif Masud, and Lawrence A Bergman. Numerical solution of the fokker-planck equation by finite difference and finite element methods a comparative study. Computational Methods in Stochastic Dynamics: Volume 2, pages 69 85, 2013.
[PR20] Vivek G Patel and Nikunj V Rachchh. Meshless method-review recent developments, Materials today: proceedings, 26:1598-1603, 2020.
[PRDVR20] Nicolò Pagliana, Alessandro Rudi, Ernesto De Vito, and Lorenzo Rosasco. Interpolation and learning with scale dependent kernels. arXiv preprint arXiv:2006. 09984, 2020.
[RC21] Alessandro Rudi and Carlo Ciliberto. Psd representations for effective probability models, Ad- vances in Neural Information Processing Systems, 34, 2021.
[RMFB20] Alessandro Rudi, Ulysse Marteau
Marteau-Ferey, Ferey, and Francis Bach. Finding global minima via kernel approximations. arXiv preprint arXiv:2012. 19788 2020.
[RPK19] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physicx, 378:686-707, 2019.
[RR96] Hannes Risken and Hannes Risken. Fokker-planck equation. Springer, 1996.
[RRTI7] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient langevin dynamics: nonasymptotic analysis. Conference on Learning Theory, pages 1674-1703. PMLR, 2017.
[Rud17] Walter Rudin. Fourier analysis on groups. Courier Dover Publications, 2017.
[RZ10] Christian Rieger and Barbara Zwicknagl. Sampling incqualitics for infinitely smooth functions, with applications to interpolation and machine lcarning. Advances in Computational Mathematics, 103):103-129, 2010.
[RZG$23] Anant Raj. Lingjiong Zhu, Mert Gürbüzbalaban, and Umut Şimşekli. Algorithmic stability of heavy-tailed sgd with general loss functions. arXiv preprint arXiv:2301.11885, 2023.
[SAA20] Szchr, Dario Azzimonti, and Laura Azzimonti. exact kernel framework for spatio dynamics. arXiv preprint arXiv:201 1.06848, 2020.
[Sch07] Schaback. Kernel-based meshless methods. Lecture Notes for Taught Course Approximation Theory. Georg-August Universität Görtüngen, 2007.
[$im17] Umut Şimşekli. Fractional langevin monte carlo: Exploring lévy driven stochastic differential equations for markov chain monte carlo. In International Conference on Machine Learning, pages 3200-3209, PMLR, 2017.
[SR20] Adil Salim and Peter Richtárik. Primal dual interpretation of the proximal stochastic gradient langevin algorithm. Advances in Neural Information Processing Systems, 33:3786-311 3796, 2020.
13

---

[SS 16] Saburou Saitoh, Yoshihiro Sawano, et al. Theory of reproducing kernels and upplications. Springer. 2016.
[SS22] B Sepehrian and Z Shamohammadi. Solution of the liouville caputo time-and riesz space fractional fokker- planck equation via radial basis functions. Asian- Bura-Eu Journal of Mathematics. :5(15))2550995 2022.
[SSDE20] Umut Simsekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu. Hausdorff dimension, heavy tails, and generalization in neural networks. Advances in Neural Information Processing Systems, 33:5 138-5115,,,, 2020.
[SSG19] Umut Simsekli, Levent Sagun, and Mert Gurbuzbalaban. A tail-index analysis of stochastic gradient noise in deep neural networks. In International Conference on Machine Learning, pages 5827-5837. PMLR, 2019.
[Ste21] George Stepaniants. Learning partial differential equations in reproducing kernel hilbert spaces. arXiv preprint arXiv:2108. 1580, 2021.
[STL9)] Gennady Samorodnitsky. Murad Taqqu. and RW Linde. Stable non con-gaussian random processes: stochastic models with infinite variance. Bulletin the London Mathematical Society, 28(1 28114))554-- 555, 1996.
[SW06] Robert Schaback and Holger Wendland. Kernel techniques: from machine learning to meshless methods. Acta numerica, 15:543-639, 2006.
[Tdd12] Eitan Tadmor. A review of numerical methods for nonlinear partial differential equations. Bulletin of the American Mathematical Society, 49(4):507-554, 2012.
[T6610] Fredi Tröttzsch. Optimal control of partial differential equations: theory, methods, and applica- tions, volume 122. American Mathematical Soc., 2010.
[TSIG22] Armin Tabandeh, Neetesh Sharma. Leandro lannacone, and Paolo Gardoni. Numerical solution of the pokker-planck equation using physics-basee based mixture models. Computer Methods Applied Mechanics and Engineering, 399: 15424, 2022.
[UHK18] Sabir Umarov, Marjorie Hahn, and Kei Kobayashi. Beyond the Triangle: Brownian Motion, Itō Calculus, and Fokker-Planck Equation Equation-Fractional Generalizations. World Scientific, 2018.
[VMR+21] Adrien Vacher, Boris Muzellec, Alessandro Rudi, Francis Bach, and Francois- Trancois-Xavier Vialard. A dimension-free free computational upper-bound for smooth optimal transport estimation. In Conference on Learning Theory, pages 4143-4173. PMLR, 2021.
[Wen04] Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
[WWLZ22] Jia-Li Wei, Guo-Cheng Wu, Bao- Liu, and Zhengang Zhao. New semi-analytical solutions of the fractionaioll fokker-planck equation by the neural network method. Optik, 259: 168896, 2022.
[XCZG18] Pan Xu, Jinghui Chen, Difan Zou, and Quanquan Gu. Global convergence of langevin dynamics based algorithms for nonconvex optimization. Advances in Neural Information Processing Systems, 31, 2018.
[XZL+20] Yong Xu, Hao Zhang, Yongge Li, Kuang Zhou, Qi Liu, and Jürgen Kurths. Solving fokker- planck equation using deep learning. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(1):013133, 2020.
[Y+ 18] Bing Yu al. The deep ritz method: deep learning- based numerical algorithm for solving variational problems. Communications Mathematics and Statistics, 6(1):1-12, 2018
[Yan13] Limei Yan. Numerical solutions of fractional fokker fokker-planck equations using iterative laplace transform method. In Abstract and applied analysis, volume 2013. Hindawi, 2013.
[ZDL22] Jiayu Zhai, Mathhw Dobson, and Yao Li. deep learning method for solving fokker- -lanr-placck equations. In Mathematicul and Scientific Machine Learning, pages 568-597. PMLR, 2022.
[ZFM*20 20] Pan Zhou, Jiashi Feng. Chao Ma, Caiming Xiong, Steven Chu Hong Hoi, et al. Towards theoreti- cally understanding why sgd generalizes better than adam in deep learning. Advances in Neural Information Processing Systems, 33:21285 -22966 2020.
14

---

[Zho08] Ding-Xuan Zhou. Derivative reproducing properties for kernel methods in learning theory. Journal computational and Applied Mathematics, 220 (0000022):6646400 2008.
[71.22] Yuping Zheng and Andrew Lamperski. Constrained langevin algorithms with l-mixing external random variables. arXiv preprint arXiv:2205.44192, 2022.
[ZZ22] Xiaolong Zhang and Xicheng Zhang. Ergodicity of supercritical sdes driven by alpha-stable processes and heavy-tailed sampling, arXiv preprint arXiv:2201.1013 2022.
15

---

Appendix
A Existence of an appropriate PSD matrix Mg (Fokker-Planck Equation)
A.1 Useful Results
Consider the mollifier function g defined as in RC221]
2-d/2 |
g(x)
(13)
Va
where Ja/2 is the Bessel function of the first kind of order d/2 and Va = Spalis de I(d/2+1)
4/3
go(x) = = = (xx/v).
Lemma ([RC2T), Lemma D.3). The function g defined above satisfies L (R(RR)) n (.R)) and (() 1. Moreover, for any w E TRd, we have
[[]])))(() S Flg((W) S [[]
The following proposition provides useful characterization of the space W.2 (RR
Proposition 4 (Characterization of the Sobolev space Wh [Wen04).. Let k N. The norm of the Sobolev space || (lw3 is equivalent to the following norm
|I [ = / (3(S((W)? (1 + [101 du, Af E (T((R))
(14)
and satisfies
(2m I (Bd) $ 11fllwg & 22 12 AS E (RR
(15)
Moreover, when d/2, then W (RR) is reproducing kernel Hilbert space.
Lemma 2. Consider the definition of gu (xx) g=(((x/v) where g is defined equation
13 Let
3 > 2,9 E N. Let Sq E W2 (IRC) n L (Rd) and the function p* = Ei= f2. Let us assume that the Assumption is also satisfied. Let & E (0. and let n E Rd__. Let Φx.r be the feature map the Gaussian kernel with bandwidth 17 and let Hxx @ Hr be the associated RKHS. Then there exists M. E (Hx & Htr with rank(M.)) g. such that for the representation plz, & x. M following holds for and k E (1....d).
@p(z,t) 8p" 1x"(x,t)
I (T((2)
II we (pdy (ffllln(R)))
Øt
Ot
(l(2d)
O.fe
2(26) 8
ll [lfillangee (ad))
I1glle, (2-));
(16)
L2(20)
Op (x,t) Øt
Op(x, t) Øt
$ (S0(22)(0-))
II (2e)) (illooo
1.((d))
Of 0t
2(22) 9
I1 billing(ee (ad))
|191l
(17)
(2-))
Op (x,1)
Op(e, 1) Or
(8)((2))
|I (ee)) W (fillegack)
Oak
(((d))
2(20)
Σ I[fill (R())
af: arr
|9 (x()
(18)
(i((Rd)
Op" (a,
Op(c,1) L)
(8T(2)) (3-1)
[1 (wdy (fillle
(19)
Oek
12(3d)
16

---

0.6
+ 2(2o)² 8
M w. (Ra)
(20)
Ora (R0)
Proof. We have g defined as in RC21
4/2 2
gu (x(x)
(21)
Va
where Ja/2 is the Bessel function of the first kind of order d/2 and Va JpalsI dx r(d/2+1))
2//2
(x( = o =00gg(//v) use result lemma
to prove the statement. For a function f W.2 (IR").
(f x go)
U5. /
Oor (f go)
Dr:
L2(R(")
L2(8d)
= 112 - F( (2(20) = 12TT Dlgo))
1m S.a. wi w [vw)( dw S ?
w (لmاا\1W de 11111
= -17 [ Jenellzs
wf + |
( du
4m W
< sup e||=||21
+1+ || dw
(1 + ||0| 2)2
wu?
4m²
s sup
228 [fll w.f (ee)
S sup
22 1lfll
+||wil21 (1+ ||ا|||| 3)
<||il|2 (1 2
16n 16r²(2n) (2v) 2(8 1)
|1 wg (R.)
(1 + x2 02)0-1
Hence,
0
16, 1 2(8-1) (1 02)8-1 || [1 (Rt)
f
(s * go
(22)
0x:
dz.
(R(R))
The same result will also hold for
&
without loss of generality. Hence,
øt
a δ
f
( * go)
S
16 (20)2(11))))
|| (K+) (1 + 02)8-1
(23)
Øt
at
L2(3ed)
We now bound
(2.t)
ap(2.1) St
and
(2,1) 0z;
ap(z,()
(
8a,
(L(()))
As in RC21 let us denote, fi * 9o as fi.v We consider
M
f.
Hence,
a Σ (x,t).
p(x = f(x,t) = (x, M.xx.((ett) =
i=1
This gives,
i
Of.
Qf..v
9.
fo 0f:
fo Ofe
Of..
fuw Of
p = 2
fi.v dx:
2
+ ft
0x:
dz:
Ox:
Ox,
Ox:
Dx;
Ox.
i=1
1=1
fi Öf:
0f.,
+ (fs - f Ofr.
= 2 Σ
dx,
dz,
im 1
17

---

Hence, we have
0f:
0fc
15..
izz p
Ow. p
2
11 + +
11 (i
Oz:
dv:
dz: L2(34)
(2d))
i=1
(2(R4)
1=1
(24)
Similarly,
0
0 01"
2
4
Of:
0f..e
1[Sallez(R4) + +2_
I1 (emlle2((*))
@fr.
01
L1(Rd)
0t
Øt
(2(8d)
0t
(2(3d))
(25)
Hence,
0p' (x,t) 0p(x,t)
9
0f.
0f
0fer
S 2
Σ
[ +2 22) 2 IIfe - fr.ell2a
01
0t
(R()
01
01
12(0d)
01 (.2(0d))
(0-2)
0fi.c
$ (1 + v2) 8-1
Σ
[I W (e + 2(2v))
II I1 (Rd)
(1((
S 82(22)
flfe w (20) (lfelllllllllaae + +2(22)
Σ [I llfall (144)
05: Jv
L2(R'")
88x(2 (0-1)
[[ W (20) (lfellllllllll + +((20)
[I Ily (2d)
aff
1190 [goll
Öt
(2(RR))
8T((2)
[| W www + 2(20)
Σ
[| w (md)
aff
1 (Rd)
at
La((d))
If choose
==min min (
x
E
=
20+1
(26)
where
C1 El= 11 [3 (R-)) IIf (ll and C2 = Ei= llfal ws (Re))
0/4 LI(Bs)
Hence,
Opt
0p(x,t)
I E.
Ōt
Ot
('(dd)
Similarly,
0p' (x. t)
0p(x,t) t) Øt
0f
afi,
1f..
2
[Ifillow + Σ II fiella Ot
0t
(2(R4)
1=1
0t
Øt
L2(2d)
(o0
Now,
0
af
0f
0fi Ot 1.oo
gu
11 1 =
Øt
(zd)
Øt
(R00
Ot
(Rd)
Hence,
Op" (x, 1)
1p(x,1)
< (8))
[| Wg (R4) (14)) (
0t
0t
(2((2)
1=1
+ 2(2c)
IIfe w (Rd)
01
(
1=1
we choose
(1
1/1
==min min
s
20+1C2
18

---

where
G1 = Ei= ll (Re) ( ((1 and C2 Ei= |l [iflll (34)
L (.9))
Hence,
Op"(2,
tp(z,t)
S.
σ1
01
(2(24)
Similarly,
Öp" (2, t)
tp(r,t)
s 8T (20)
Σ
fi w/ (aoy)lffilll
(27)
Ork
80p
L'(R4)
Ofi Dor Lg(299)
+ 2(22))
| wflllw (Rd)
1l0ll,,
and
Öp (x,
Op(x,t)
S (T() (8-1)
Σ
II ll wy (Re1) (ete)l/llllllmmee
0xk
Bark
(2(RR)
Of: Oxs (=()))
+ +((2v)
|1 I1 (R-)
Lemma 3. Consider the definition of = g(x/o) where is defined in equation
[3]
Let B > 2,q E N. Let fi.....f E W2 (I2) n L (Bd) and the function p" = >'= f.- Let
us assume that the Assumption 2 is also satisfied. Ler E (0. 1] and let n E B Let σx.r be the feature map of the Gaussian kernel with bandwidth 1) and let Hx @ Hy be the associated RKHS. Then there exists Mg E S+ (Hx @ (T) with rank (Me) S such that for the representation p(a,t) t) = =f(xx) = x followinng holds for all i (1,...dd..
0²
8²
16m(2u)(0-) 1) (1 ((9)l))()))
05m Dzy (2(a)))
p" (5,
dr,ün, p(z, t)
Σ II (Bet)
(1+ (1+ 4
(1(3))
m=1
32 32m) (8-2) (20) (1 + 22)
9 |I (Re) 115ml (2(20) + 2(2v)
8 fm
|I 11mlllg (Rd)
Z-2
(((()))
and
0²
8² dr, p(r, t)
(2) (0-1)
18
(5,
(1 (1oll
Σ
flll (R()) (Rd) Der (Ro)))
(1 + 02) _
(2(2))
m=1
32 (1-2)
(20)
9
pm
|1 (Ro) 200 + 2(26)
[m I (xd)
8-2 (1 22) 2
Ox, (xd))
Proof. We have y defined as in IRC21]
2-4/2
(x)
(28)
Va
where Jaj2 is the Bessel function of the first kind of order d/2 and Va = S de = r(d/2+1) (x(x) = (g((wx/)). We use result of lemma to prove the statement. We have, for a function
4/2
E W
8² dr,0x,
0²
8 0x,;dx;
a² 1 - Ox,0xx
1
(f + go) dx,
=
($ * go) 122299
(2(
|| w [[ 1x= ww gell 22
| (Iw - Fla.))ll (22
19

---

16T
i|l - F[o)( de
JRd
S $16T' [
afw31111)(E))
=16m"
wi + | ( [I de
< <16m"
(1 [[IWW1))STT du
e||l||21 (1 + |
wfw
< 16
sup
(1 + |
22 fllligg (Rn)
e||=||21
< 16 sup
(1 + | - 2 ||
fl (ar)
o||w||21
256 (25))(8-2)
|| w
(1 (1+ 12)2-2
As in RC21], let us denote, ft 9o as fe.. We consider
M Σ Stsle
Hence,
((x, = f(x, = (xx 1) (x, = Σ (z,
i=1
Now we have,
8
82
05 0fm Oz, Oz;
0fm,v fm.. dx; Oxz
Ox, (x,
dx,dx, f(z,t) t)
Σ
=
(A(Rd)
m=
fm
8 8p
+fm 0z,02;
fm,e dx,0xz (1(d))
05 05m.v
05m
05m
05m,1
15
< 2
+
Oxi dx:
Oa, (I((RR))
Ox, dx; (R((Rd)
Ox,
(2(Rd))
(.2(8R)
m=1
8 fm
fm.e
02 fm. ax,0x L2(34)]
+ Σ
[lfmll ( + Il.m fm.elli2(RR)
dx,dx, 0x,0x, (2(22))
m=1
Similarly,
82 (Disp"
02
dr,0r, (Di, (Dissa,
dr, On,
LI(R
82
02
I sup
D
p² (I,
0x,0x, f(x,t)
(29)
dr,wn,
( im
(RR)))
R2
Now,
8²
82
p¹ (a, t)
fl, t)
Ox,Ox,
dx,0x
('(Rd)
05m
0fm.v
05m
05m. 0fm.v Ox,
Ø 0f
2
Dr:
dx:
Oay
dx, LP(Rd)) Ox, (P(R))
(2(Rd)
(2(2d))
m=1
20

---

9 0
8 fm.. dz L2(34)
02 fm. Ox, L2(3))
+
Σ
[l (I2((R)) + Ilfm -moolle2(K))
Or,Ūs,
Let's consider all the terms separately.
0m
0f dry (2(4))
(2u( (3-1)
O5 Öer I2((R))
11 Wl I1 (Rd)
Ōr:
(R=)
(1 + u²)
x
d-1
mmm
05m
mf..
05m.v
16 (20) (8-1)
0fm I1 wg
Σ
5
m=1
dz,
(2(30)
Oxy L2(RR)
(1 + v2
0
m=1
ary 12(00))
9
8 fm
0 fm, 0x,0x, 12((R))
327 (0-2)
9
2 >
[l LP(8d))
|l W (24) llfmlleaRR
Ür,0z,
(1 + x2) 2
m=
8 fm, 0r,0zz (P(BR)
02 fm. Ox, (2(R2)
2
1l sm, (2(8d)
5 2(2x)
Σ
fl Il (B4)
(B.)
m=1
m=1
4 =2(20)
02 fm Jr dr,Os, (3(33))
11 Il (
m=1
4
02 Jm 0x,0r, 122(R))
= =2(20
fl ll FWw
Hence,
||
8² 0z, p¹ (x, 0x,0x, f(e, t)
0²
16m(2u) (0-1)
85m Σ |I II ( Oxx, 1,2(2d)
(1 (1ollaea)))))
m122))
(1 12)
32sh (2v) (9-2) 9
82 8m Orrciaa (2(2))
Σ II I/ (Rd) [l ea + 2(2v)
Σ (fmlllg((e)) (md)
1191l
(1 + 02)
Hence, if we choose,
E ((15728-1c)))
e (
# (2842C5
v min
26
where
Ci
(1
([a))))(R))
29h
|l/l (R.)
C2 119llz
D2,
(R())))
Em=n (f 1l. (a then,
and Cs
Em llf Il I (Re)
2² fm fm dadda
(2(d))
82 1" (x,
0²
f(x,
Se.
Ou,
(((R9)
This gives,
82 Dij dz, (x, Dij dx,Ox, $(z, t)
S dd Rae.
((((R))
In similar way, we can also obtain bound for L metric. Now we have,
0² 0x. p² (x, 1) dr,Os, (xx,
82
0fm 05m
0fm. Ofm Ox: 0x,
Ox Ox,
(2(RR))
m==
21

---

80
8 fm.
+ +Im m dr,00,
Sm Ox,
12((2))
0fm
ffm.v
0m
0fm dx,
0fm.y dx,
0 f dzj
2
Or;
dx;
(RR)
L (324)
(R(R))
1000 (Rd)_
4
8 fm
fmm Or,,Or (R((R))
02 fm. Or,00.r
Σ
[l [[mlllZ0 Lo0 + 11 Sm,wll
Or,On
(B4))
Let us now compute the following,
O 0m,"" Ōuz (RU)
05 Ozj (R()
119llz)(R4)
and
@ sm. 0r, 0x, (Bd)
02 fm 0x,0x, (.0)
Similarly,
y dx,0z, (D,p"(Ixt))
0² 0x,0z, (Disf(z)))
j=1
(2((2)
82
82
sup Dis
p² (x, 1)
dr,0a, f(x,
(30)
dr,ūs,
j=1
(2((3))
:=
Combining everything together, we have
0² du,
0²
(67(
9
8fm axz
p² (x,
fxx ()
(1 I1lll
Σ |I I1 (wff(R))) (Rd)
dx,du, (2(30) (1 v2) _
m=1
32sh 32x) (9-2)
9
0 8
Σ I1fllwg [lfml (o) + 2(2)))
Σ llf (f
I|alle.
(1 x2) 2
Oz, (20))
Hence, we get the final result by replacing fxx, t) with with
A.2 Approximation Result by a PSD Model (Proof of Theorem
Proof. Now, just need to collect results from Lemmas 2 and to get the final bound. Let us assume that there exist an optimal density function p" (z, t) which is the solution of fokker- planck equation. Hence, We would like to show that there exists psd operator M, such that if p(z, t) = fzx t) = Φx.T(x; t) M,0x,T t) such that we can obtain bound for following three quantities for 21,
Op" (ar.t)
ap(a,t) Lr(Rd+1) = 0ss).
1.
dt
dt
2.
d θ
(2,(2,t)
Op(a, t) Ox;
(Hi(z, t)p² () ()0x,))
Σ
[i(x, t)
dz,
1=1
(T(R-11)
Opi
(p" (x,t) p(2, p(a,
-((s)
Ox:
(I(R--1)
3.
>d mm
Σd ==
(Disp*(z,
82
(Disk(z)))) Li(Re+1) = 0(s).
Bartay
22

---

If the above three condition holds, then, we can show that for r {1,2},
0p(x,11 1)
0² dz,dzy (D (Dis((e;)) ())
(1,(2) t))((2,))
0()).
Ot
Or,
(=(R0=1)
We have from Lemma 2 and equation
26). if we choose
min 2 16
20 28+1C2
where
C1
Eis Il (fllwffeet)) ( /
and
C2
Eis (1
0/s dt
[19 (R Hence,
L (R
Op" (rr 1)
(f(x;
S.
0t
Ot La(eet))
From RC21 Theorem D.4 (Eq. D.25) J, we have,
|l wy (3 (0ll (,, 1) p(ar, 1) (2((4111 S (2 + l Ia (BO+L))
This proves the first part. Going further,
a (1, t)p"tt) (x, - pa(a, t)((z, t)) 12(R0+1)
d Σ (x)
(x,
Op(az, t)
Σ
dr.
Ox:
Ox:
1ui(z,1) (p" (p"(x,) --(aa; b(o.:))] L2(20t1)
Ox:
(x,))
ap(a,t)
Ops(z,, t) (p' (p'(x,) p(x,,)) (2(20011
Σ [x(2,)) ()
oz.
oz.
12(24)))
dr.
Op (x,
Öp(a, f) Or.
< sup (up t) Σ a.1E((,9)
Ox.
(2(2+)
SRu
(x,
+ d sup 2.1e(2.7) i
sup
||I" (x, Ō(1, (1leaead))))
Or:
SRAP
d
(x, t) 0p(x,t)
< R2 Σ
+ dR [1 (a, Ō(a, (1222
(31)
Oa:
Oe:
i=1
Now,
0
Σ
(M t)ptt)p' (x, ()((x,))(zz) t))
dx,
(2(Re+1)
Op* (x.1) 1)((,,1)
< R
dRRu IIp' (z, ((z)) f) Lqquet)
Ox; dx; L
From RC21 Theorem D.4 (Eq. D.25) ], we have,
9 (x, p(z, ()ll, S (2v) l (R
And previously, we showed in lemma
(x,t) Ork
0p(x,t) t)
5 STT (8-1)
[I [ln (iilee/tiilllllllerrre
(32)
dzk (2(dd11)
23

---

2()2) 8
|| I1 (R2 1)
der (
(33)
If we choose,
& Jt : ((2--1667)))
(
(2stica) (20ticc))
e -mminn
where
Ci
Eisr IIfellly (ifelttnllllll...xrrett ceee |19 (1
C2
Ei=1 IIfillwz( 88
0%
and C3
(1
(1((())()н))) Eis Ilfellug (Rans)) (Rats) f (R(att)) then, we have
d
Σ
(mi(2,, t)p" - ml(x, 0) 03)
S <d(Ry + RugJe.
dz;
122(8011)
We showed in l that
02 Oz,00z p" (xx,t)
0²
dx,dr, p(z, t)
L2(B011)
16T(2c)
0m
(1 (2941))
Σ I1fll (ad+1)
(1+ 02) +
Örz
1x (Q441)
32 (2u) 22 (0-2)
q
9
fm
(lllg I1. ( + +((20)²
11 15mllls (Re+1)
(1 02) 2
ac. Dxr
(xd+))
m=1
m=1
Hence, if we choose,
V min ((
3
F
#
e
s
26
28+2C
where
Ci = (1 + (19ll
Eim= ||ʃ wy (Rd+1)
8fm Jey
C2
L (
Em [I (R4+1) 1l م (pet1) and C33 Em [1mlllg then,
2² fm 0z;d0
(191)
82 dr,Ür, p² (xr,t)
8²
d2,0x, (z,
Se.
12(Be+2)
This gives,
ap(x,t) Øt
t)
d (((x,)))
à
ΣΣ 82 dr,dz, (Disple,t)
+ dz,
< d(Rp +Rp,) + d² Ra).).
1=1
(2(R411)
We can bound the trace of the matrix as follows. We have,
M
Σ f.
1=1
From [C22| Equation D.21J, we have
Tr(M)) = Σ
1 = =c,,28(11 +
(1= 20 2 e 2
f
3
Wil (R++11
1=1
Since, we choose,
E
26 2C
24

---

where C C
C1 c1 (19111911)))11) C2 Em= |I [ifflll:] Iw (R
ox, ( C2
Em= (0llw I1 (R and C3 Em= I15 (w 0² fm fm
40 (1 (2)0 Tr(M)) c,,2
evonr e
Σ
fall (alliwg (R-12)
96 2C 1 1 = =
8u
||
W3 (R++1)
i=1
20 1
c, 28
3
e n (z)
W. (B++1)
E" .
=|11///3
1 e p J
20
S
C/ml1/2
1 exp
exp
-
-
where Č 927 -C, = S9 and and =dd2222 E1m
|1
W.
(B4+1)
max
1, 3 29
A.3 Compression Result
Let X Rd be an open with Lipschitz boundary, contained in the hypercube [-R, RR with R > O. Given ....... Fm X be m points in [-R, RJd. Define the base point matrix X e Ip mxd to be the matrix whose j-th row is the point Tj. The following result holds. We introduce the called fill distance [Wen04]
h = max min || - 2/|. xel- R. =EX
(34)
Now, we state some useful results below.
A.4
Useful Results
Lemma 4 ([RC2]) Lemma C.I ). Let = (-R, and n, E R.... Let u E Hy satisfying u(ā1) u(ām) = O. There exists three constants c, c,C, C² depending only on d (and in such that, when S o/C", then,
particular, independent from R, u.
[|]| S Cap e -
-- log ah ||
with g, = det( + diag(n)) and a = min( RR, -:i
and T/- )
max:=l..... T/-
Lemma 5 (Norm function's 's derivative with scattered zeros). LetT= (-R, R2d and E Ed u(fm) = O. There exists three constants c, C, C' depending
Let E H,, satisfying u(ō1) only on (and in particular, independent from R, 7. u, .......,,,,, such that, when h S a/C', then,
DDuullz=(T) S Ge-co/h log(co (odco//]llullx,,,
with G be some positive constant and o = min(R, R, and n4 = max;=1....dd T3-
Proof. Here, utilize Eogliardo-Nirenbbrd inequality Nirll We have, for 0 < k <
oo and m) + 2k < 1, then
D' S Gluull 1-0()0))
D /|| (²(T))
(35)
where G 22/2a///-11///T(//P ko
Considering the case for j 1 and j = 2. we have,
J12 (
|D ull & G1llu||
25

---

| S G21llull (2(T)
ul (T(
where G 1//2://4-1/2//T/PPP
and
412) 212 (
G2 24/2///4-1/2// (14/2)11
We here have assumed that +
"-"² a+4
sin di44
2
From the defeinition, we have
DD'ullz 1 Gillull D (2())
( 1 G2llullia(7) 1- D (2()))
From the definition,
d+3 D D'ul u 1 (T)) Gi Gi |
wyy (T)))
d+1 (ry
D ( & G2 (((())
|
Now by Theorem 4.3 of [RZ10] we have for some constants Ba and B' depending only on d,
Byk*
||]]I2(TT) k!
h W
when kh < R By Hence, we have.
Blgk
| (' SG1
h² [|
B ht ||
D ( S G2
k!
From equation C.9 in RC21 Proof of Lemma C.I], we have
[|]llwt(r) S ||w|bor, 2/2 (47) )/22 1/2
Hence,
Bakk |D''ulL=>rr S
heth 1ullo,,, Co) 1/lloor,,,/ (42 k/2 /k2 = G (2 (2/n+Bahte Bahk)* )/2
||ullac,,
k!
/k!
Bhk** uullac, 1/2
(4n+) 2/2 6/2 VRl G2 (2 Bahk)**
|D S22 k!!
h
||u||ss,,
vk!
Now, we make similar choice of the parameter f and k as in RC21 If we choose, C3 max(2B.,BB) h < Cs 4+4 min(R, 1/ & [s] and = C. min R, 1/ /n+)), then we have hk S It and
2.TT- Moreover, VRl log & If also assume that G max(Gi, G2), then for
j = (1,2)
ID'ulLLSTT S G
VR! V |ullxx, S G/4 1or///ooo//////llllaaa/a
The final result is obtained by writing s/4 = 00/h with o = min (R, 1/ (T+ and C3/4.
Lemma (Lemma 3, page 28 PRDVR20). Let X Rd with non-zero volume. Let H be reproducing kernel Hilbert space on X, associated to a continuous uniformly bounded featurr map X H. Let A K. be a bounded linear operator. Then,
sup [Ab(x) xEX
sup fllc())
11
In particular, if X C Rd a non-empty open set, then supsex (| suPpsii<s ||A (
S
26

---

Theorem 6 Theorem C.3 (RC21)). Let R > 0, E +++, E N. Let CT be a non-empty open set and let T....., associated to the points w.r.t T (defined in eq. (34). Let P :H, + Hy be the associated projection operator (see definition in eq. (38) ). There exists three constants C, C", such that, when h S a /C',,
......... Emm be a set of distinct points. Let h be the fill distance
sup aex
(1 ()),(())la,, I Cqn e
= log
Here 4n - det( + diag(n)) and a min(R, / 11 max; Tie The constants C', C" depend only on d and, in particular, are independent from R. 17, .i.....,,a Im.
Similar to the results in theorem 6 we will obtain results for derivatives.
Theorem 7. Let R >> 0. - E N. Let X TT= (-R, R2S be a non-empty open set and let .i,.... Fm be set of distinct points. Let h >> 0 be the fill distance associated the points w.r. T (defined in eq. (34)). Let P :H + Hy be the associated projection operator (see definition in eq. 388 There exists three constants c, C, C², such that, when h S o/C",,
sup |I(I (B)DjPPDD (x))b,,,, C1qy e- 3f log 42 ZEX
for some positive constant Ct. Here 9n det(, 1 diag(n)) 1/4 and a min(R. (or nt M = max; TB- The constants c, ,,C',C" depend only on d and, in particular, are independent from R, Ii,...., Imm
Proof. If the function lies in the RKHS, then using the property of the projection operator P, have
((Pf())=D((PPf Di = (o(.))))))))) = D' ('. (f.)())))))))))))))))()()))))( (f. (f.D'ó(z))) ())z).
Hence, PDf-PD/S ((x))(a))>>>> 0 for all E E (1,... .m) By lemmal 4 know that there exist three constants e, C, C" depending only on d such that when h S a/C" we have that the following holds ||u]| (ullix(() I Cgp e of log - for any u E H, such that u(ō1) = ... = u(im) = O. Since, for any E Hyy we have that DĀ f P belongs to Hy and satisfies such property, we can apply lemma 4 with u = (I DS f, obtaining, under the same assumption h.
|I ( fT) S Cqp C Tf log of DDD
Vf E HG,
where we used the fact that II(r P)DP flx,, S S PIll flloe, and || - PI| S 1, since P is a projection operator and so also I - P satisfies this property. Derivative reproducing properties for kernel can be looked at in [Zho08]. The final result is obtained by applying lemma 6 with A = -- P. from which we have
sup [I(I - (DDPP)D''(()lI S sup =EX
[I(I 1(1-P)D'SlL-O
Cgy e log || (DSflos fl
sup [Diflis]
||
Cign e log o
for some positive C1. The same result can also be obtained by the use of Lemma5Jwin condition that f is d/2 + differentiable function which is a slightly weaker condition on the order of smoothness.
A.5 Compression of PSD model for Fokker-Planck
We have (01, (1), 22 ... (Zm, tm) boundary of the cube [-R, Rj x [0, X be the base point matrix whose j-rows are the points Zz. Let us now consider the model p = M(, v) and & = s(: X, v) where
x T which is an open bounded subset with Lipschitz
Am = K(X.TTXXTY ZM K-1 (XIT)XTT)
Žu = tT) Tu, u.... (( T
(36)
27

---

and
2 a Σ (i Wa)s.
(37)
Hence, for any E RR
A Σ
4,==
Let us also define associated projection operator,
p= 2" K (X,T)xX Z.
(38)
Some properties associated with it,
p2 P, p7 2-, 2p=2.
Theorem 8. Let E R4+ and ler M E (3G). Let be an open bounded subset with the Lipschitz
boundary of the cube [-R, RJ [0, RR, R O. Let Ø1,..., Jm E ' and X be the base point matrix m. Consider the model f(:; M, 1/) and the the compressed model p f(:; An: X, 17) with
whose j-rows are the points žj with j = 1,...,m.
Am K-1 8.8.
2MŽ K 8,8
(36) in terms of Xm- Let h be the fill distance (defined in eq. (34 associated to the points Tts.... Jmm The there exist three constants c, C, C' depending only on such that, when S a/C", with a = min(R, 1/ /++); 7 = max;=l,...d TKis Qn det( 1 diag(n)) 1/4, then we have 7t
where 7. HG, - Thm is defined in eq.
1(p(x,1)- ((, 0t
d (M,la,,))(((.z p(x,
d ΣΣ im1jm1
8 dx. dr,dr,
D.,
(p(z,t) p( t))
dr,
S dR + 2dR,, + 4d Ra)l dee - lo + (2 (2+ dR + 2dR + 2d Ra) /IIMilp(xr t) Cgre 50 log 22
+ +2
te [IM M, + (2H,, + Cdu, - +
(1;
+ 2R oe log te
where Fx, (r; M,
ov(r))) Oz;
M
00(z.1) D2,
Üz (x)
(1
P)
00(z.1) 02, 10
faye, (,
Pulz))) M
v(x.)) and ūzz = (x) =0| -T-
0² ((x.))
Proof. We have,
p(x,t) - p(x,t) (xx,t) (PMP)
(39)
Let now consider every term separately.
Öp(x,t) 0p(x,t) 2 Oel(x,)) Mú((,t -2 Ou)(x,t) PMPW(1,
0t Øt
Ot
Ot
Orl Ovl (M (MM))((2))) Øt
(40)
Similarly,
0p(x,1) 1) Or,
0p(r, 1)
Mi(xx,t) - Ax((,, T PM PW(xtt
= dr,
0z:
Or:
28

---

(x))
2
(M - (MM-P)))(z,t))
(41)
dxi
Finally,
pp(x,t) pī(x,1) (2) 2 (x,
22 (2) (2
dr,ez dr,zi
(2. PMPU,( (
(2/ t) (MP
(x, 1aatt))))) (M - PPMP)((x (z, 6) (M TM (MMP)
(42)
We have,
PMP - M- P)M - P) M(I - -P) - PPMM
Since la T AaAl < Aa| 3e |B] and Ja ABa] < ||]||//]| |///||||BB//aa||| for any a in a Hilbert space H.. and for A, B bounded symmetric linear operators with B E S++ (HH). by bounding the terms of the equation above, we have for any E Rdd
(W(,, ((x,t)((PM))) MM t)1 - 21((T ()),(2, 1)b MI ( tll,,,
+ I1 P)b(2,, (x, |M
=20 fxx M, // (x)) ca (x)² 1/2
(43)
where CM =| and we denoted by ūz the quantity (1))) = Il(I P) t)b,,, and noted that ||M t) = (xxx t) Mab(a,, f(z; (x, t)). Now,
(M - PMP))(E 1) S I1( 1---P)),,, 111 || (1-
B)) @u(a,t) t)
Da:
Oc:;
- -11 1-P) 0v((111)
- (1 (1- P)p 1) 1011,, ||Ar]1//2
M1/2
Bv ()
Oz,
Ox
1944
56,
(x)ū (x) + cis f(x;; M, w) 1/2 u.n (z) + 1/2 fax (x;; M, M, 2/) (z)),
(44)
where Tz. (I; M, v))
@v(x.)
M
@0(2.)) and ūzz (x) |( - (p) -
av(xtt)
oz,
Uz
Following similar steps, we get
tv( T
(M PMP)6(0,) f) I {c)((())ū(())) + 1/2 fl; M, (z)
2/2
Ōt
CM fz(:; M, eb) 1/2 (1 (z)),,
1/2
(45)
where M,, (1((a.))) dt
M
(wh().))) dt
and ūx(x)) = (1I P)
(a((a.)) dt
Kr
t) TM CMMTMMMPMMM PMPP (MP))(1, 1) (x) (x) + CM [x; M, bb) 1/2 22 ir ar (x)
1/2
1/2 F fayZ 25 (x; M, 2b) 1/2 (x), +Ca
(46)
where Jazar (2; M,1))
M ( and ū.fo;( (x) = || - p) 8 ((a1))
0zy 1r.
And,
(x(I,t) (M (1) < Caruz, (z) (x) - 1/2
CM fz, (2; M, 16) (x)
+ Ja (x; M, 1/ // (x).
1/2
(47)
Hence, combining everything together, we have:
0p(z,1) 0p(z, 1)
Ous. ()(()))) aar(a, 1)p(I, 1))
d D
82
0²
p(x, t) Di,
2x,0x, p(2,
Øt
Øt
Ox.
Ox.
Fx. Fn,
29

---

1p(æ,))
Öp(z,() ()
0p(z,1)
05(x,1) aon
d
Opi(a, 1)
Our(a, 1) do. p(a, 1)
Mi t)
Me( t)
dr. p(a,t) tt
öt
0t
do:
d Σ
d
82
8²
+
Du
Ozz p(a, t)
Dus drn, (ar,))
1=1 3=1
@p(2z f)
0p(z, t)
Op(a,,
0p((x, ()
<
+ R
+ (ol,,
(0(2x, p(2,
Öt
Øt
Ox.
0x.
0
0²
Ra
p(a,
dr, 5(2,1) 1))
Or,0r,
Ovb())))) T
2
(M - (MMP)), + 22 Σ
(M - PMP))(a, f)
+ ddR,,l((, (MM (M ParP))(ar t)
Us.
+ 2R. (M PPMP)ME,, (i, 1) (M - PMPPOS 01
(M - MMPPMMMGT ()) +2 Σ ()( 0t
(v,))11) T
2
(M PMPP)e(o 1)
+ (d (M- (MrP))(Z, f)
1=1
+ 2Ra oi (z, 1) (M PMP) (1, + 2R (Ta
PM P) ))) (x,
Now, combining everything together, we get
0p(z,t) t) 01
ap(z,t) t)
(0, t)p(x, t))
0( t))
0²
Dis Dr dz, dr,
p(a, 1) - Dis dr,0o, pla, ()
Or:
du:
2cMun (x) 20 1/2
f(z: M, / (x()) + 20 1/2
M f((;; M, 1)) 1/2 M1 (x) + d dR (20 1/2 (x; M, 1b) 1/2 in (x) (x)))
+ 2R
Jevun (1r uz (x) cM f(;; M, vb) 1/2 Üz (x)
CM fm, (z; M, 1b) 1/2 (z))
1/2
+ 2R
[eau ūi (xr (z)) (z) + 1/2 (x; M, xb) 1/2 uzyar (x) + 1/2 Jayar (z; M, ab) 1/2 122 (x)
+ 2Ra
Jew ūzy (z) (z) 2/2 So, (a; M, v)) 1/ u., (x) + chr Te,(2; (a; M, 3) "u, (x)
1=
Now from theorems
6 and 7 we have that when the fill distance h < T with a = -miin R2 f), then | (x)) S Cape
% log K (| S Cape % log % [ 1| (x)) S Ca % log xr (48)
and نمحنت (x) I Ca of log o
for all i E (1.-..,............ dd, E {1, ,d and some constant C, c, C' depending only on d. Hence, we have
A(p(z, () ()x,
(luaa,,((((ez p(x, ())
8²
D,y dr,Or,
(p(z, t) p( f))
dt
dz;
j=1
S + dR, + +2dR,
+ 4d Ra)l 4 M Ehe 200 log x + (2 dR + 2dR + 2d Ra) 2d) /IIMIpp( t) Ogge "f log x
te log n + 2C
VIMIL M, (2R, Adn The of
[I;
+ 2RaC o log so
A.6 Final Approximation Bound (Proof of Theorem 2
Proof. From the previous result in theorem 1 know that there exists M2 E ( & Hy) with rank(M.)) s q. such that for the representation p(z, t) 6x.r following holds under
30

---

assumption
on the coefficients of the fractional FPE,
Öp(z, f)
2 dr dz, ple, t)
(ua ()(,,),) 6)
Dis
= -(e),
(49)
Öt
ax:
im1jm1
12(8))
for an appropriate choice of parameters. We also know that
0p" (x, t)
0² Dis dr,ur p" (a, 1) =0
(pi(a, (x, ())
(50)
dL
dz,
We know that, for given Am- p(x, t) = (xx t) PM,PPuaa t) where P is the projection operator defined earlier. From the results in theorem 8
we have
0(p(r, - ((a,)))
Mal(z, ()((,))) 1) ())
0² (p(2, t) - p(ax, t))
D., dn,dn,
0t
Ox:
(2 + dR, + Rd., 4d Ra) Ra)l | T5e log th (2 dRR + 2dR, 2d Ra) v IMM lIp(2, t) Cqne
+ +C R los ViMellf((a; M, + (2R + CdRa)Cane R. lox (M;
+ RRa o log o
[Mellfaas, Me,
It is clear that some big constant Q, Impo Mc, n) S &IlMe for all i,j E (1,...,, .). Similarly, fz. M.. 11 S M for all i .4) and (r(;; M.. s @ll- Hence, we have, 0(0(x, - -p(a,))
M(ula, (p((x,)) - p(x, ())
8²
D&s Or,Os,
(f(z, t) - p(ar, t))
Øt
0s;
(2 + dR,, + 2dR, 1d Ra)l M. C o2 200 log
dR + RdR. 2d Ra) / M l|p(z, t) Cane log To |Mal] + Rad Ca.e # log ||Mall
T8 log 28
+ +2QCq,e -e- log th |Mall (2QR,,
4Qd Ca
< (2 + dRy + 22dR, 4d Ra) M. &002,e
4 - log 4 + (2 dR + 2dR. 2d Ra) R.) /IMM,Ilp(e,) Cyne log
+ +Q22 + 2R eQd Ra)C
S8 log s ||Mell-
Now, let compute |Me]l. We have
MMall < Tr(M)) = C-(4+1)/2
(1 :5 x2))
E'e- 2
Let us choose,
Hence,
3-2 log( R.
23 d+1 |M S C2(1 + 1R)++R - 8-2
Now,
0(j(x,t) - -5(x, t))
a(u (a, (f(r,)) ((x,
@²
Dij dr,du, (p(x,1) - (le,1))
Øt
de.
12(3))
(2 + 3dR, + 4d R.)] C 20 -200 log p/d+1)/2
Q(2C + 2R d 6Qd (oddRa))Cap Ra) oee log Se ||Me] Rdd+1)/2
+ (2 dR 2dR, 2d Rd) [Mellpp Cg.e 2 log %f
= (2 + 3dR + 4d Ra) [M, C dhe
-20 los v 1//21)/2 Q(2C + 2R, d + 6Qd Ra) Cane o log o0 |Mal] ((4+1)/2
P1
+ 3dR, + 2dl Ra) |l [MMllppa (2200 Og,e" 4 log 4
P=
= =//l|Ma]l02 dhe - loz tे (d-1)//2 + P&Cque T log 5E [MMe] |Me|je(e-1112
+ P "22 d log * Il p(a, (lle.))lee 12
(51)
31

---

Then, note that ||p [2)))) = (1 Hence, using the argument in [RC2] we have
lp 11 - 2R.
Term R comes from the fact that for each time instance p(æ, 1) is a density and 1 E (0 < R). Hence,
D(f(x,t) p(x,,))
((,,(, (0((x,), ;(x,1)))
4
0
dr,da, (f(x, 8) p(z, 6))
01
12(3)
$ Pill MM thz 200 log th Hdd P.Cqre log * |Mee [Me|lall-1112
+ +2 11 Cane * log %
P1 R) 3 C² th
- log th n(d+1)/2 P2Cane T log T. Cs(1 + R) R))= # (1411)22
1+l + +2,3RC, (1 + R) 3-25 20-
Cane o log %
By choosing h = ca/s with max (C', (1 +
d+1 log 1 + (1+ log(1 + R) + log((() +e),
2
for some big enough Č. Since & 2 e, then log 2 1, so
Ce % log %f Ce #log v < Ce S CC( R) d/2 2(n-2).
Hence, we have,
Op(xx,)) 1) p(x,t)) 1))
p(x, 1)))
8²
Dis dr, (f(x, ((z,
da;
12(3))
0(e).
Combining evrything together and by the use of triangles incquality, we have
Op(x, 1)
az, (m, t)p(z,t))
ΣΣ Dos
dz,dz, plz,t)
O(e).
Öt
12(3)
Em is h-covering of T, guarantees that the number of To conclude we recall the fact that 21 centers m in the covering satisfies
m S (1 + 2Rya+1 d+1
h
Then, since h ca (C4 log Cs Rebg((+RR with C4 1 + d/ min 2) and C3 (C)1/C, and since a min (R 1/ 1/VT), then R/o = max RVT) s 1+ 1/18-118-2 (log 1+ 112 so after final calculation, we have
क < 11 11122R0++// + log + Clue x log
1 R
1 R
1/2
for some constants Cs Cg and C1o independent of e.
A.7 Proof of Corollary
Proof of Corollary 7 Let p" be the probability corresponding to the solution of the Fokker-Planck equation with drift f and be the conditional probability corresponding the solution of the Fokker- Planck equation with drift jt. And denote by (x, L) the m dimensional PSD model obtained by applying the construction of Thm. 2 to the Fokker-Planck equation defined terms of μ. Hence, from the direct application of the result in [BRS16] Theorem 1.1] since |||
Ail L. = O(s) we have,
I|p' -1l O(e))
(52)
From, theorem we have,
| - p.y 11111X(Txx) O
Finally due to proposition 2 for any PSD model f there exist algorithm that samples i.i.d. from the probability satisfying
O(E).
32

---

Now note that for any t. the function p. (t, is still a PSD model (see [RC21]). Then, for any t (0, T), we have that there exists probabilities Pr.o such that
Wi (f. )), (.) O(e).
Then, since W a distance and is bounded by the 1.2 distance since we are on bounded domain, there exists C' depending on the diameter and the volume of the domain, such that for any
Wi(p" (4, (.)) x C(IPP" (4, (1, (1,-))222(x)) + (1b((,,) - (2((4-)12((x) + Wi(
Then, by taking the expectation over t (O, we have
E. Wi(p" (t, (t.)), < C(IIp* 00 [[ax(xx + IIø(t,-) -t, Pa)
B Estimation of Fractional Laplacian Operator on a Probability Density
Before go into the details of the estimation, let us first recall bochner's theorem.
Theorem 9 (Rudi7)). A continuous kernel k(z, y) = k(æ - y) on Id is positive definite if and only if ((6) is the Fourier transform of a non-negative measure.
Proof of Theorem
In psd model approximation, we have
p(z) = Ax;
Σ
()
15=1
We will first prove the part (i) of the result. Now, we know that, fractional laplacian operator is defined as,
F((-A)*1({> ((E()=)((A(EE)
=
4==1
= =|||||
Alaaa, (1E) F[k(es,
Now that, & a gaussian kernel with bandwidth n.. Hence,
k(y, x)) ~111~=11
Flk(,, )](E) = it y J. I1El1
ce J. (1&111 4jng y T 4o
2|aا|| ce -elllllll" e ll 2imell
for some e > which can be computed in closed form. Now, we know that convolution of two multivariate gaussian functions are still gaussian. Hence,
Flk( ()()) F(h(2y, ({S) c'e
for some constant e' > which again can be computed in the closed form. Hence, we have
F((-A)*[(E)
2
Σ
AyJ (E{)) Flk(2,, ()())
f,f=1
m
A ||||||
(53)
33

---

Finally we have,
Σ Asj / [16]220 e +1l=s)1)) stl 2in(a,t;;)1f e -58" dg
(-A)")(x)- C"
25=1
Aye Mye / 1 J5 के 8 2==n(()tzz))1 de.
C"
By inspecting the above equation carefully, we can see that after appropriate scaling,
(-A))f()) - Σ
Age nlles ++ N( 22,)=
(54)
1,j=1
where μ 2jn(z; + 22)) and Σ = 4nI. One can sample from N(p, Σ) to estimate the fractional
laplacian operator for PSD model.
Now, we will prove part (ii) of the theorem statement. From the bochner's theorem, have,
k(r ( y) = / q(w)e (a-w) (e-s) dw,
where g is the probability density. In our psd model approximation, we have
m Σ
p(x) = (xx,)
1-
Σ Ais (Jre jus (or soi--al 1) glwsr) duss) (se ejui (as (a)--)) ((ww) dwo
(
( cdul (ec-x) a ( (e, ( dlwselan. ne
Aij
Rd
Now, we know that, fractional laplacian operator is defined as,
F((-A)* f(E) (T(((E))))(E))
Σ Ais I / edint (ac- (o) (eg==)((u)g((u))) dweddur e da 1,j=1
||<|| 2.
m
Ass Iec ( (Jre couT (as -at didl
= ||E|] 2+
dx g(we)g(wx) dwsedor Ais [ I (ʃ estuut antwaas) e jo (wetent)) dx a(we)g(wx) duredon
i-=1
= ||E|| 2x
kJ=1
(fe e-deTwurtuntt (ustontt) dx) g(we) dusedwa
||<||
28
Aij IRd JRd
Glwe x: )ctwers))
(.j
1
Aij
I
cdlu a. a.twae)) 20l120((ur + we + E) g(wz)g((x) disddok-
1
JRd JRd
1,j=1
Hence,
Hleo aitweey) ( 151220(ww + Wik + E)e Te as) glwe) dwrduix
(-A)*f(x)
= Σ
Aij Sge Jzn
=
Aij j 2d
edlut as aitwers) wwwe + wi 2 (we+ūr) g(we ((we))(w))) dwordon
1
Σ
Aij I ( clwi (mi+ awn (es+ + wx ((we)g(wx) dwopdor
Rd zd
34

---

As Bl[lus + wxll22 2 (a, (artt)dda (1stt)))))))
+=
C Existence of an appropriate PSD matrix M₂ (Fractional Fokker-Planck Equation)
Let us redefine the mollifier function defined as in [C221].
2 2-d/2 d/2
g(x)
Va || (a/2(2n) (1p22/2))))))))))2((((((1111)
(55)
where Ja/2 is the Bessel function of the first kind of order d/2 and Va Si dx r(d/2+1) (x) = - (00//(//)) F-) denotes the fourier transform of f.
x4/2
C.1 Useful Results (Fractional Fokker Planck Equation)
Lemma 7. Consider the definition of (x:) =00((//oo where is defined in equation 133 Let
2, >> and let 1j E R44.. Let Ωx. be the feature map of the Gaussian kernel with bandwidth 1j and let Hx @HT be the associated RKHS. Then there exists Mg E S (Xxxxx &HT) with rank(M S q, such that for the representation f(x,t) 2x Me x, following holds for O and f1....d) d}.
E N. Let S1 W.2 (R4) L and the function p² = Ei f?. Let E (0,
I1(-A)'(" p S (20)
Cl fl w (ec+1) + (20)
Σ C fl (eitly
for some positive constants 6 and ez where Clu ( (1TAA pla and c( = =111:11==0(/))))00>==0====
|||||| [f |واا| ER for i f1.....gh ,q).
Proof. As in RC211, let us denote, fi * Io as fi.. We consider
M Σ St.
Hence,
= ((x,() = (x, 1) M, Σ 1,,
I((-A)(pp -))1LA(R ((((-))"(p" D) [1((211))
||| ||||||2001p* 1) |||u|/ee0
Σ
1 12 (w)
L2(R011)
9 E(.
| [|0leag 22g
- f.,
Faw) (w)
1|
F F((S [1,0)(: (w)
12(R211)
(2(2112)
Σ [||W|| I(f - fl. + f(. (W) =
|| ||||o|12 [ fa,e) (w) * + F + fr. (w) (w)l
Let us consider the function.
||||| 25 of [(ft f(,w)) (w) (w) / ||/| [ -y)(()))) de
(w)
:=)(w)
35

---

/ || y + y)n(v) du Sa f IW - 11 22)(()) dy + 02 / [lu1| 2))(()) dy = n(w) + com(w) * ([(w)]=)(ww)}
Hence,
II(-A)(p" -00l"a--lllla S <cllllw/ n(w))l + c2llm(() (wwll~n((w)] (w)))l (R(0)))
Now, we can apply Young's convolution inequality
incuuly on I2 and T2. Hence,
Ti S c
|
2m
m
|
+ |
2 m
I Sei ||| |
26 (2)
[[]T[A]] + [
=cil
| - (2) [[[[[A]]] (R + -
S c:
|||| (0 [[ + ( [T
In the last equation, we applied holder's inequality. mann-Lebesgue lemma, we have
From uniform continuity and the Rie-
| (fgg))ll ( & ( = / (o(vx)) dx [g(g///)) [1(x)) de llal v (((et1)-
Hence,
Ti S ci [ ( =| 100
Similarly,
Ts S (allm (Z)) + +
S <c>llm(w>)l>>>>>>
=c)llm(w)]l<><e) [[]+ L (xe+1)
In the above equation, we can apply holder's 's inequality to get,
To & 22 (2llm(w]]l [[][=|]>([[[((]]]]]]]] (x++1) + |||||| 2 F |1 [x(
$ [a)lm(w))l [|] + ||be||=*"((f))))))=======: (w)) (
Hence,
T2 I cgllm(w))ll
[[]] () (1e+1) ||
The term [|m(w)]l II 2 and [|||||?= (2(3)13) can be bounded using similar technique as before. We will use the result from Lemmal
|||w|| ((W))
[ I|l (ww) du
||||| 48 | [+]( dw
/ J11w12
[|o]|1b| (1 + ||ا|| - + (T)))))))))) du
||:0 as
228 sup x||=||21
(1 + || 12)0 ||
W (Ra11)
< 220 sup
Il
=||w|||1 (1 + || 12)3-22
Wy (Re+1)
2(8-22)) 26)
220
(1 1= llf
W2 (R&-1)
36

---

22 (20)218-23 26)
11bbl
W (1dtt) S (20) (226))(0-2))
[[lill (f (10t2)"
(1 + 02)8-20
Combining everything together, we have
- IK--A)(0--Plleeaaee -P)1 I (2 (8- 26)
Σ Ti | W/ (B441)
(34+1) + 02 (2
Ilfil Wwg (
22
where
chi |||||| (f ( +
1 + || ~ ( (011))1 and CT2
|||ي [A 2 Leretty
C.2 Proof of Theorem
From the previous lemma, we have,
9
(K--))(p* ( S (2v) (8-00)0-000 2n)
Σ
c
Ti
[[
We (2at1) + (20)
Cit [Ifill
w3 (a--1)
==1
for some positive ei and c where c
| and ||||||| [|llo//>==f]([\\@l====>>))l | D2(R111) Now, if we choose, for some constant
|F |
Chi =|I1|1==((()WW)/lrrrr
c28 2s Ei= (CH, ( flll
CT,
w. (Re11))
then
(p(-A)<pp* - 1) Sc.
We can bound the other two terms using lemma We have,
0p" (x,
0f(r;; t)
S (2T) (8-1)
Σ (f (
q + 2(2))
0f.
[| Villwa((eet) (Re+1)
11 (
Öt (Be+1))
If we choose
min ((00 (22-116mC
x
28+ C2
where
C1
Eim II Iilll ( and
C2
E. (R))))
|1 (151,,, Hence,
- (Rd+1)
Op (ax,t Op(a,t)
Se.
(((R0+1) Ø1
Similarly, from lemma 2
(x,
0p(z, t) dx
< 8T (2v)
|| (R (illlo (eal
(Ra11)
Oxy
(((444))
[5 || ||fill (wet-))
+ 2(2)) 8
119ll,, I10l,,,
37

---

From RC21 Theorem D.4 (Eq. D.27) I, we have,
(1,( p(a, t)]l 23) I (2u) + lgg (9(lo)())41)))
II w (0 >
1=1
Hence, we choose,
min (( (28--1677)
7
#
28+1
2
where
C1
Eis1 Ilfe (enkll/e
C2
=
E' [1fillag(ROs))) d/
0%: (e 1191le
and C3
(1
|19|l (R Eis llfillwg (2011)) Il (Re1)) then, we have
2 (1, (x, - ull ((x,t))(z,))) t))
S dRR, RupDE.
Oz,
L2(B011)
Hence, by combining everything together, we have
Op(z,t) t) O1
Σ
(mile 6)) - -(-a)"
= + +((R +Rp)))
dr.
(2(20+3)
We can bound the trace of the matrix as follows. We have,
M = > Si
From [R22| Equation D.21J, we have
q
Tr(M)) Σ
89 Ma c,,220 (1+ ()23
e n0
fall
w (R412)
Since, we choose,
7-28
263 2C
where C depends on fis, F(Ji] and it's derivative.
Tr(Me) = 2, (1+ (g)
llf W2 (1)
1=1
28
1 1 2,22
25
J
80 agu
Σ
f Wil (Re+1)
16T 967228-2
2,
I| Wiy (3d+1)
1 3
e nn (z) ( 3
1=1
( Ö
_ e पंअ ))
2/2
1 g
1 2/2
2 exp
23
170
where
C
96.228-2C.
8
and
C
π -d/2228 D |l
Wi (R
max
1,
3 28
87
38

---

C.3
Approximation Properties of PSD Model for Fractional Laplacian (Proof of Theorem
Proof. From the previous result in theorem 4 we know that there exists Mg E ( @ HG) such that for the representation p(x, t) = x following holds under assumption 2 on the coefficients of the fractional FPE,
with (M)) S
Op(x,()
(M(a,, (mi(2, f)f(,, ()) (-))))(z,))
= =(s),
(56)
0t
dos
12(30)
for an appropriate choice of parameters. We also know that
(x,
Σ
0 Or.
(ms(,, (x, + (- (z, =
(57)
Öt
We know that, for given Am- p(a, t) (xx t)T defined earlier. We will utilize the Gagliardo Nirenberg inequality to control the error compression of fractional laplacian compressed model with respect to the original one. For any fixed 4,
t) where P is the projection operator
Eegliardo-Nirebber
(~ p(z; t)] (-A)"M(<z,)) (PM_ M2) t)
By applying Gagliardo-Nirenben inequality for fractional laplacian IMP18] we have
((-A))B((()))>a())) S (1(((2) (aaa
(1) TmT2
where 0 = 8 1 and A is normal laplacian operator. Let us consider the term Tj and T2 separately. In our case, (x) = (zx (PMM M.) (x, t). In the proof of theorem 8 we have proved the bound on A (z)) and Be (x)). From equation (43). have,
(xx t)T (PM, M,)o(z,t)) t) = 21((1 P) t) lac, [IM& 1/2 /IMM (/z, t)ll,,
|I - ()y,(r, (), |M&ll
2c 1/2 fx;; Me. (6) /2 (x) + cxrū ( (x:)²
where сы = ||M&l and we denoted by ūn (I) the quantity (z = 11(I P)) (z, lll,,, and we noted that |IM|/ (zx B),, t(z,t) (x, f(x; MG, (x, t)). Similarly, we have shown that p'((,t) pp(att)
(x, t)_ (M.
2V((,)) PM. PMes) t). (58)
On, dr,z,
From equations (46 and 47
(x, ( (M. PM₂ P) t) s
(x) + cxl fa: Me.
(x)
1/2 f0zz, 1z, (z: Me, 16) 2/) (x), (59)
+
where Eajar (x; MeM 43
0 ((.())))) dieyas
M
0 ((rt) and uzyan (x) = (1 - Po ((x())
Jayz.
Jzz 36.
And,
1)
PM. cxruz, (z)ū,, (x) + G 2/2 (2; Mes 1/2 (I)
1/2
+0 Ja, (a;; Me, 2/2 (ur (2x)
1/2
(60)
Follwoing similar steps of calculations that is Theorem 8 and Theroem we get a(p(z, 1) p(xx 1)) +
5(2, t)) a(m(z, ()(, - p(a, t))) (- (x (x,e)
0t
Ox:
(2(3)
39

---

Pa IMclII h log % R(d+1)/2 + P.Cage %. los t [M, |
PallM.11' # log # 0le,0)1
(61)
for some constant P1, P2 and P3- Now, let compute ||Mell. We have
|Mel & Tr.M.) Cy(411//2
C' log( )
Let us choose,
Hence,
23 2-
28 |M&ll S C3(1 + R) 25g 22
d+1
For some constant C3. Hence, Then, note that 111/211(2(33) in RC21]) we have
|| (xx) Hence, using the argument
1 ( S 2R.
Term R comes from the fact that for each time instance p(x, 1) is density and / E (0 < R). Now,
0(p(,, x. p(x, t))
a(u (f(rr) ((z, + (-A)" (c. ()(((x v(x, ())
0
De.
1,2(8)
P1Cs(1 R) 8-2 8-20 02
-20 log R R P&Cqne log Cs(1 + R) 2-)-20 8- R(d11//2
qne
PPC R(1 R) s- z(-22) Cane log IVp(e, O1l
all
By choosing h = ca/s with max(C", (1
d+1 200-20
log _ + + () log(1 + R) + +lgg(C) +e).
for some big enough Ć. Since & 2 e, then log 8 2 1, 50
Ce - log 4
logs s Ce < C' Ri d/2 11 x(8-2))
Hence, combining evrything, we have from triangles inequality,
A(p" (a, t) p(x, t))
0(ui(z, () (2,t) p(æ, t)))
+ (- Maa f))
O
Or:
(2(3)
0e)
To conclude we recall the fact that 21, En is h-covering of T, guarantees that the number of centers m in the covering satisfies
(1 2fvd+1 ydtt1
h
Then, since h 2 ca / (Ca log Co log(1+R) with C. 1 + and since o min(II 1/ /T), then R/o max(1, (1, (1,) < 1118- (log R so after final calculation, we have
22), and Cs 1C),
6
1 R mat < 51+ 211 Cy log
1 R 1/2
+ C1oE チーzZ log
E
E
for constants Cg, Cg and Cio independent of
40


Reviews:
[{"review_id": "z3lttDsA8f", "reviewer": "Reviewer_biEz", "summary": "This paper studies the problem of efficiently sampling from a SDE given a drift function and diffusion tensor. It proposes a solution based on computing a PSD model that satisfies the Fokker-Planck solution associated with the SDE, and then sampling from the resulting PSD model.", "questions": "- Can you provide experimental evidence for the theoretical claims?\n- Can you provide simplified theorem statements?\n- Can you provide a more clear comparison with related work?", "limitations": "Yes", "rating": 3, "confidence": 2, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "The paper seems to have a lot of technical content, and it is clear that the authors have put in a lot of effort.", "weaknesses": "This paper is extremely difficult to read. The theorem statements (and the paper overall) need to be drastically simplified to be understandable. The problem doesn't seem to be well motivated, and it's not clear from reading the paper how the solution compares to others in the literature. For instance, the [LCL+21] paper is mentioned in related work, but it's not at all clear how the solution proposed of learning a PSD model compares to their work. There are no experiments to support this idea of fitting a PSD model. Overall, I think this paper needs to be rewritten to be more easily understandable.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for spending time reviewing our work. We agree that our paper is on the technical side; however, we suspect that the reviewer might have found the paper difficult to read if the subject is not close enough to their field of expertise. In such a case, we believe that the strong rejection would be an unfair suggestion, as the paper would be clear and understandable and the contributions would be significant for its intended audience at NeurIPS. \n\nWe would now address the concerns raised by the reviewer here. We also request the reviewer to read our responses to other reviews to get further clarification. We hope that the reviewer could reconsider their decision based on our explanations.  \n\n**Theorem Statements**: Theorem 2 and theorem 5 are about the approximation error guarantees of approximating the solution of the Fokker-Planck and fractional Fokker-Planck equation with the PSD model. More specifically, the statement says that we get $O(\\varepsilon)$ approximation error given enough samples $m=\\tilde{O}(\\varepsilon^{-(d+1)/\\beta-2s})$ where $\\beta$ and $s$ are described in the paper. Corollary 1 provides a guarantee about the sample quality of samples obtained from the PSD model and the error guarantee is provided in Wasserstein-1 distance. Theorem 1 and Theorem 4 are intermediate approximation results for the infinite-dimensional PSD model. These results are useful in obtaining the final approximation results in theorems 2 and 5. Finally, the result of theorem 3 is about approximating the fractional Laplacian operator acting on a probability density represented by a PSD model. Theorem 3 clearly states that while using PSD model-based representation for probability density, one can approximate the fractional Laplacian operator (non-local operator) acting on a probability density for a wide choice of kernel functions. Approximating non-local operators is an independent research problem where we show the effectiveness of PSD-based representation (theorem 3). \n\n**Motivation**: As we described in the introduction of our paper, the main motivation behind this work is to propose an approach that can provably sample from a stochastic differential equation. We take an altogether different approach to address this problem that was not explored in the literature before to the best of our knowledge. We propose a two-step procedure. First, we approximate the solution of a (fractional) Fokker-Planck equation (using the PSD model) that corresponds to a stochastic differential equation.  In the next step, a sample is sampled from the PSD model. Under our proposed approach, we managed to prove the sampling and approximation results under much weaker conditions on the solution density functions ($\\beta$ times differentiable) than what exists in the literature.  \n\n**Advantage and Comparison**: To the best of our knowledge, we are not aware of any algorithm which provably can sample from the solution of a Fokker-Planck in a bounded domain under as much of a weak condition on the solution($\\beta$ times differentiable) as ours. The convergence bound for existing sampling algorithms like Langevin requires dissipativity, log-Sobolev constants, or similar quantities. For example, very recent works ([Lam21] ; [ZL22] ) are able to provide convergence only under dissipativity-like assumptions.  We make no such assumption and still show that there exists an algorithm that can be used to sample from complex solutions of (fractional) Fokker-Planck equation in a ball given the solution is regular enough. If the solution is $d$ times differentiable or more (that is often the case), the dependence on the dimension vanishes up to logarithmic factors.\n\nIt is clear from the results in theorem 2 and 5 that the sample complexity result in our work is better than the one reported in [LCL+ 21]. Also, the solution function obtained using   *Physics-inspired neural networks* (PINN) [RPK19], *Deep Ritz method* [Y+18], and other existing machine learning methods are not bound to stay positive in the domain and are often not integrable. Since the solution function is a probability density in our case and PSD based model satisfies the property of a density function (section 3.2 and [RC21]), hence comparison with results for methods like PINN and Deep Ritz Method [LCL+ 21] is not fair as those methods are not applicable in our case and can not be used in sampling from the solution.  We have tried our best to include existing related results from PDE and sampling in the related work section of our paper. \n\n**Experiment and Contribution**: There is a significant lack of literature on solving the (fractional) Fokker-Planck equation and\nsampling from the solution. Our framework not only addresses this problem but also allows for the approximation of the fractional Laplacian acting on density, which is a distinct research problem. We believe that our work represents a valuable contribution to the field and has the potential to inspire further research in this area. The theoretical results presented in this paper are non-trivial and come with worst-case guarantees, which makes them stand out from methods that rely on stronger assumptions. Our findings are particularly significant in situations where traditional assumptions, such as log-Sobolev and Poincare fail. Despite the absence of real-world experiments, we hope that the reviewer could re-evaluate our theoretical contributions from the perspective of not requiring strict assumptions and the novel contributions we make in multiple research problems.\n\n"}, {"review_id": "4lvljrCM4s", "reviewer": "Reviewer_oceV", "summary": "This paper considers numerical methods for the Fokker-Plank equation and its \"fractional Laplacian\" version. These partial differential equations describe the evolution of densities of certain random processes: in the case of standard FP, the process is a regular diffusion (an SDE driven by Brownian motion), whereas the fractional version corresponds to a process driven by $\\alpha$-stable noise. Accordingly, the authors focus on initial conditions given by probability densities. \n\nThe idea behind the paper is to adapt the fairly recent methodology of positive semidefinite models (PSD) to this PDE problem. PSD models estimate nonnegative functions via feature maps $\\phi:X\\to H$, where $X$ is the domain of the function and $H$ is a Hilbert space. Given $\\phi$, the nonnegative function of interest is estimated as $\\langle \\phi, M\\phi\\rangle$, where the PSD linear operator $M$ is chosen from the data. The whole method is based on the same idea as the well-known kernel trick: there is no need to explicitly deal with the Hilbert space $H$; only the kernel matrix over data points will matter.  \n\nIn a nutshell, the paper gives a method for approximately solving the two aforementioned PDEs via PSD models. The salient features of this approach are as follows. \n\n1) The PSD estimator can be obtained via a semidefinite program (as shown in previous work). \n\n2) The method will output probability densities, and one can sample from them via methods from previous work. \n\n3) The complexity of the problem grows like $\\approx \\epsilon^{-d/\\beta}$, where $\\epsilon$ is the target accuracy, $d$ is the dimension and $\\beta$ measures the degree of smoothness of the solution. The exponent is not quite right, but the point is that sufficient smoothness will avoid the curse of dimensionality.\n\nOne may note that this method has rigorous error guarantees. This sets it apart from other PDE solvers based on machine learning ideas, such as Physics Inspired Neural Nets (PINNs). \n ", "questions": "* Have you attempted experiments against other ML approaches, or other methods from the literature?\n* Methods have been proposed recently to speed up SDPs, including randomized sketches  (https://math.paperswithcode.com/paper/scalable-semidefinite-programming). Could these methods be applicable in the present setting?\n\n", "limitations": "They have not addressed limitations explicitly, but I do not think that that would have been necessary. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 4, "strengths": "The paper gives strong convergence guarantees for a new PDE solving method. Although several ingredients come from previous work, they are combined with new ideas into an impressive package. As mentioned above, this is a ML-based PDE solver with rigorous guarantees. Importantly, these guarantees do not require coercivity or hypercontractivity properties that are usually required from Langevin samples for solutions.  ", "weaknesses": "* No experiments are performed to compare the method with other PDE solvers (though this is understandable in a heavily mathematical paper). \n* The writing is a bit sloppy in places: I found quite a few typos and have several minor comments on writing (I could probably find some more typos if I kept looking). See the list below. \n\nNB: were it not for the sloppiness, I would have given this paper a higher score.\n\n---\n\n_List of typos and comments on writing_\n\n(Line numbers refer to the full paper in the SM.)\n\nLines 146/7: \"It was\" and \"It has\" should be \"They were\" and \"They have\".\n\nLine 180: why define $j=\\sqrt{-1}$ if you could have just written $\\sqrt{-1}$ in the exponent? (You never use this notation again, do you?)\n\nLine 185: Shouldn't the condition for infinite $p$-th moment be $p\\geq \\alpha$?\n\nLine 201 and elsewhere: I find it a bit odd that the transpose notation is used for vectors in Hilbert space. If you are going to use it, at least explain here what is going on.\n\n \nLine 251: why is this sensence needed here?\n\nLine 205: $\\eta$ was treated as a vector before, and will be treated as a vector later in the text. Here, however, it seems to be a scalar. \n\nLine 223: \"The cost of the algorithm is...\" -- per sample point, right?\n\nLine 254: \"Hence, $p^\\star(x,t)$\"\n\nLine 255 should end with \":\".\n \nLine 270: the definition of $\\psi$ should be in a numbered display for easier refererence.\n\nDisplay (8): if $p^*$ is the true solution, why do we have it in here? (Its contribtion vanishes.) Indeed, the definition of the loss later does not involve $p^*$; see eg. (9). (On the other hand, $p^\\star$ shows up again in the display right after line 989.)\n\nLine 294 and elsewhere: when you are mentioned a numbered theorem or assumption, start the corresponding word with a capital letter. \n\nLine 295: $\\mathcal{H}_X\\otimes\\mathcal{H}_T$\n\nLine 318: \"by Gaussian kernel approximation\" -- should this be here?\n\nLine 365: \"Similarly to the previous section, two steps are required to obtain\"\n\nLine 327 and elsewhere: I think $\\mathbb{E}_t$ was not defined anywhere.\n\nLine 352: \"If the kernel satisfies the Bochner...\"\n\nLine 355 and elsewhere: this is a matter of taste, but couldn't you replace \"utilize\" with \"use\". \n\nLine 376: \"The major difficulty\"\n\nLine 391 (taste): \"does require\" could be \"require\".\n\nLine 395: I think \"approximated\" should be \"approximating\".\n\nProposition 4 is also in Evans's book, which was cited previously (I am saying this from memory, so please check. $W_2$ is replaced by $H$ in that book.)\n\nMany spots of the SM: the inequality $\\|f\\star g\\|_p\\leq \\|f\\|_p\\|g\\|_1$ should be mentioned, as it is used in several steps. \n\nLemma 2: it seems that $\\tilde{p}$ is replaced by $f$ in the statement of the Lemma and at several steps of the proof. (Eg. the display right after 665.) This also seems to take place in the next Lemma.\n\nLine 680: Since you are repeating a definition from the preceding proof, please say so. \n\nLine 682: \"We use the result of Lemma 1\"\n\nLine 746: Should \"The following result holds\" be here?\n\nDisplay (34): this \"fill distance\" is the smallest value of $h$ such that $\\tilde{X}$ is an $h$-net of $[-R,R]^d$.\n\nDisplay below line 777: how should I interpret the norm in the RHS?\n\nLines 873/4: why are you spaeking of conditional probabilities here?\n\nLine 891: what is $f(\\delta)$?\n\nLine 910: \"redefine\" should be \"define again\" (they are not the same).\n\nLine 928: capitalize \"Holder\".", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks a lot for assessing our work positively. Here below we address the major concerns raised by the reviewer.  \n\n**Typos and Presentation**: Thank you very much for pointing us to typos. We will fix the typos and improve the presentation. We will also include suggested citations in our paper.  \n\nLine 251: If a probability density function has a limit at infinity, that limit is always zero. We mention this because the PSD model-based representation of the density vanishes at infinity. We will clarify.  \n\nConfusion regarding transpose in RKHS and $\\eta$ will be clarified. \n\nLine 777: For our case, this is just the supremum of $f$ in the set $\\mathcal{X}$. \n\nLine 873: We will fix that. That would be just \"probability\". \n\nLine 891: $f(\\delta)$ is $k(\\delta)$. \n\nOther minor confusion regarding $\\tilde{p}$ and $f$ would be fixed. \n\n\n\n**Contributions and Experiments**: There is a significant lack of literature on solving the (fractional) Fokker-Planck equation and\nsampling from the solution. Our framework not only addresses this problem but also allows for the approximation of the fractional Laplacian acting on density, which is a distinct research problem. We believe that our work represents a valuable contribution to the field and has the potential to inspire further research in this area.\nThe theoretical results presented in this paper are non-trivial and come with worst-case guarantees, which makes them stand out from methods that rely on stronger assumptions. Our findings are particularly significant in situations where traditional assumptions, such as log-Sobolev and Poincare fail. Despite the absence of real-world experiments, we hope that the reviewer could evaluate our theoretical contributions from the perspective of not requiring strict assumptions. Thank you very much for pointing out a reference for faster SDP via sketching. This method can certainly be utilized in solving the SDP program arising in the PSD model faster.\n"}, {"review_id": "cHBmGgAc0G", "reviewer": "Reviewer_T6Bn", "summary": "The authors introduced a new method of sampling solutions of the Fokker--Planck equation (FPE) using a positive semidefinite (PSD) model. While the algorithm used for the PSD model is not new, the authors created a framework to sample from the FPE, which is quite original based on my knowledge. \n\nWhile I believe the results are very neat, I would like to clarify several questions with the authors before providing a confident score. For now I will recommend borderline accept, and will raise the score once my questions are adequately addressed. ", "questions": "1. To start, while I understand the algorithm of [MFBR22] is not part of the contribution for this work, can the authors provide a high level explanation of how the algorithm works, and comment on its computational complexity? I believe all readers can benefit from understanding the core algorithm, and plus I believe readers should be able to implement the algorithm presented in this paper without having to reference other works. \n\n2. Can authors clarify the role that $q$ plays in this problem? In particular, is $q$ given or chosen by the user, and how does the value of $q$ affect the computational complexity of this algorithm? \n\n3. Under the discussion of \"Efficient sampling method\" the authors mentioned a step of finding the best model $\\hat{A}_m$ by minimizing an error function. Suppose I am an applied scientist who wants to sample from a given FPE, how would I compute this model $\\hat{A}_m$ from the PDE? On a high level, I would like to understand how to implement this sampler end to end, and I don't think I fully understand this part. \n\nOnce again, I would be happy to raise my score once my questions are adequately answered. ", "limitations": "N/A ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The approach to drawing samples for the Fokker--Planck is very unique and original. I believe alternative methods that hold promise should be welcomed and more carefully studied, regardless of how practical it is at the moment. ", "weaknesses": "There are a few clarifications I would like to have regarding the implementation of the algorithm. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for assessing our work positively. Here below we address the major concerns raised by the reviewer.\n\n**Sampling Algorithm**: We will add the details of the sampling algorithm in the revision. Here below, we are describing the methodology of [MFBR22] for sampling algorithm for PSD models on a finite hyper-rectangle. The algorithm is an iterative sampling procedure. Given a bounded hyper-rectangle $Q$, $p_Q$ denotes the PSD models-based density in the hyper-rectangle Q. Details are in [MFBR22]. Given the function $f$, the algorithm will take three inputs $(Q,N, \\varepsilon)$: the hyper-rectangle $Q$ (with\nsides parallel to the axes) from which we would like to sample, the number of i.i.d. samples $N$\nwhich we would like to obtain, and a parameter $\\varepsilon$  which defines the quality of the approximation. at first,  $Q$ is cut in half in its longest direction forming two sub-rectangles $Q_1$ and $Q_2$. If $X_Q$ were a random variable following the law of $p_Q$, then $X_Q \\in Q_i$ with $p_i$ (can be computed in closed form) and $X_Q | \\{X_Q \\in Q_i\\}$ follows the law of $p_{Q_i}$ Therefore when looking for a sample from $p_Q$,   one of the two smaller sub-rectangles $Q_i$ are chosen randomly with probability $p_i$ in which we look for the sample and then call the algorithm recursively to get a sample from $p_{Q_i}$. By applying appropriate stopping criteria that depend on $\\varepsilon$, we get a random sample in $Q$. More details about the algorithm are provided in [MFBR22]. The computational and statistical guarantees are provided in Theorems 1 and 2 of [MFBR22].   \n\n**Effect of q**: From Lemma 2 in the appendix, it is clear that the approximation error grows linearly with $q$. Hence, the sample complexity will increase linearly with $q$.  Intuitively, as $q$ grows, the complexity of the target density function increases and hence, it would require more samples to approximate the target density using the PSD model. However, we are not allowed to choose $q$ as it is a property of the target density. For example, if the probability is bounded away from 0 then q=1. (For further details, please look at: ''*Second order conditions to decompose smooth functions as sums of squares*\", Ulysse Marteau-Ferey, Francis Bach, and Alessandro Rudi)\n\n**Algorithm**: The algorithm to sample consists of two parts. In the first part, we get the approximate solution of the corresponding (fractional)-FPE using the PSD model and in the later part, sampling can be done using the method described in [MFBR22]. There are various ways to obtain the PSD model-based approximation for the solution of the (fractional)-FPE. Here below we describe a couple of ways. \\\\\n\n(1).  If the drift $\\mu$ satisfies the regularity condition as in [BRS16, Theorem 1.1], then the drift function $\\mu$ can be approximated by a function $\\hat{\\mu}$ in Gaussian RKHS and the approximation guarantees are due\n to interpolation result in Lemma 4. Once, we have $\\hat{\\mu}$, we can solve the following constraint SDP to obtain the ${A}$. \n\n$\\min_{\\tilde{p}(x,t)} {\\left\\\\| \\sum_{i =1}^d \\sum_{j =1}^d D_{ij}\\frac{\\partial^2}{\\partial x_i \\partial x_j} \\tilde{p} (x,t) -  \\frac{\\partial  \\tilde{p}(x,t)}{\\partial t}  - \\sum_{i =1}^d \\frac{\\partial }{\\partial x_i}(\\hat{\\mu}_i(x,t)  \\tilde{p}(x,t) )   \\right\\\\|}$  such that $\\tilde{p}(x,0) = p_0(x)$.  $p_0(x)$ is the initial density and the norm is $L^{2} (\\tilde{\\mathcal{X}})$ norm. The above optimization problem is clearly a constraint semi-definite program because of the representation of $\\tilde{p}(x,t)$  as a PSD model. Given $m$ base points pair $(\\tilde{x}_i,\\tilde{t}_i) \\in \\tilde{\\mathcal{X}}$ for $i \\in \\{1,2, \\cdots , m\\}$ of space and time (that is selected uniformly over the ball of radius $R$), we have, $\\tilde{p}(x,t)$ given in equation 6 of the paper. Hence, the above optimization problem optimizes for a PSD matrix $A$. Now, the only remaining question to answer is if the norm $L^{2} (\\tilde{\\mathcal{X}})$ in the above optimization problem can be computed in the closed form or not. \n\nSince, the representation of $\\tilde{p}(x,t)$ is in the span of the product of two functions in Gaussian RKHS and $\\hat{\\mu}$ lies in Gaussian RKHS by approximation, hence $\\\\|\\cdot\\\\|_{L^{2} (\\tilde{\\mathcal{X}})}$ can be computed in closed form. Result in Corollary 1 of the paper is presented under this scheme.\n\n(2). As a second procedure, we can follow the approach described in section 3.2 in [RC21] and optimize an empirical version of the problem that is  \n\n$\\min_{A \\geq 0} \\frac{1}{n} \\sum_{k=1}^n {\\left\\\\| \\sum_{i =1}^d \\sum_{j =1}^d D_{ij}\\frac{\\partial^2}{\\partial x_i \\partial x_j} \\tilde{p} (x_k,t_k) -  \\frac{\\partial  \\tilde{p}(x_k,t_k)}{\\partial t}  - \\sum_{i =1}^d \\frac{\\partial }{\\partial x_i}(\\hat{\\mu}_i(x_k,t_k)  \\tilde{p}(x_k,t_k) )   \\right\\\\|}^2$  such that $\\tilde{p}(x,0) = p_0(x)$ and $\\tilde{p}(x,t)$ is given in equation 6 of the paper.\n\nIn the above optimization problem, $(x_k,t_k)$ for $k \\in \\{1,\\cdots , n\\}$ are iid samples from the hyper-rectangle of side $R$ and  $(\\tilde{x}_i,\\tilde{t}_i) \\in \\tilde{\\mathcal{X}}$ for $i \\in \\{1,2, \\cdots , m \\}$ are $m$ base points pair of space and time (that is selected uniformly over the ball of radius $R$). As long as $(x_k,t_k)$ for $k \\in \\{1,\\cdots , n\\}$ remain iid, the sampling procedure of $(x_k,t_k)$ would not change our theoretical result. Hence, we can also choose uniformly at random. A small amount of regularization can be added to control the norm of $A$. However, this procedure would require us to study the error arising due to empirical approximation in the problem (Theorem 7, RC21). Since the focus of this work is to study the approximation properties of the PSD model in approximating the solution of the (fractional) FPE and sampling from it, we do not discuss the empirical approximation method in the paper. \n\nWe hope that we have answered all your concerns and that the scores will be raised further. \n\n"}, {"review_id": "GA03QzNuWj", "reviewer": "Reviewer_jfmV", "summary": "This article investigates the approximation of the Fokker-Planck equation and the fractional Fokker-Planck equation using positive semi-definite (PSD) models. The authors demonstrate that for a given accuracy $\\epsilon$, there exists a PSD model with a dimension not exceeding $\\epsilon^{-(d+1)/(\\beta - 2s)}(\\log 1/\\epsilon)^{d+1}$ that approximates the corresponding partial differential equation (PDE). By utilizing this result in conjunction with the PSD model-based sampling method proposed in [1], the authors are able to provide samples of the stochastic differential equation (SDE) corresponding to the Fokker-Planck equation or fractional Fokker-Planck equation at any given time $t$.\n\n[1]Marteau-Ferey U, Bach F, Rudi A. Sampling from arbitrary functions via psd models[C]//International Conference on Artificial Intelligence and Statistics. PMLR, 2022: 2823-2861.", "questions": "1. In line 180, should it be $E e^{j\\langle u, L_{t-s}^\\alpha\\rangle} = \\exp\\{-(t-s)||u||_2^\\alpha\\}$?\n2. In line 216, should it be $Y \\in \\mathbb{R}^{n\\times d^\\prime}$?", "limitations": "Despite the theoretical guarantees provided by the author, I believe the biggest limitation of this article lies in the lack of clear demonstration of the advantages of the PSD model. It would be beneficial if the author could explain how the PSD model stands out compared to other similar models in terms of approximating FPE and subsequent sampling tasks. This would significantly enhance the persuasiveness of the article.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This article demonstrates clear and logical writing, making it easy to understand the problem formulation considered by the authors as well as the fundamental concepts and properties of the PSD model adopted.\n2. As a newly proposed model, this paper presents a novel application scenario for the PSD model, which has not been previously explored in related works. This contributes to its significant novelty.\n3. From the perspective of sampling based on SDEs, the authors provide a new design approach for sampling by combining existing works, offering a fresh perspective on sampling techniques.", "weaknesses": "1. Although the author mentions it in the paper, the influence of initial conditions on the solution of the PDE can sometimes be significant. It would enhance the completeness of the theoretical analysis in this article if the author could provide an analysis of whether the PSD model can satisfy the initial conditions or, in cases where exact satisfaction is not possible, characterize the approximation of the initial conditions' impact on the final solution's approximation.\n2. Personally, I have limited knowledge of non-log-concave sampling and numerical modeling of PDEs. However, it would be beneficial if the author could compare the effectiveness of the PSD model-based sampling method with other traditional methods (such as Langevin Monte Carlo) in terms of sampling performance. Additionally, comparing the approximation methods of the PSD model for PDE solutions with other existing approaches, either from the perspective of complexity analysis or purely experimental results, would provide valuable insights.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for positively assessing our paper. Here below we address the concerns raised by the reviewer.  \n\n**Approximation Error with Initial Condition**:\nWe bound the $L_2$ error of the difference between the true value and approximated value of the (fractional) Fokker-Planck equation. However, to achieve this, we need to get the approximation error of all three terms individually in the equation and that also includes approximating $p(x,t)$. Please refer to Eqns 39, 40, 41, and 42, and the block of equations after 47 in the appendix where we bound each term individually. This continues in the rest of the proofs. We get the same guarantees for $p(x,t)$ that we get for our PDE which is similar to (Rudi and Ciliberto, 2021, Theorem 6). Hence, the initial distribution will not change any approximation result in the paper. It will only act as a constraint in the optimization problem to obtain the parameter of the PSD model. Also, please see our response to reviewer T6Bn regarding the algorithm. \n\n**Weakness, Comparison, and Limitations**:   Previous papers about Langevin Monte Carlo and projected Langevin monte carlo algorithms relied on assumptions such as dissipativity, log-Sobolev, or Poincare inequality, which are not always satisfied and we do not require any of such assumptions. Furthermore, we are not\naware of an analysis that covers such a broad class of solutions that we cover in this paper. We only assume that the density function is $\\beta$ time differentiable. This assumption is much weaker.   The existing result for fractional Fokker-Planck is even more scarce. We are not aware of an approach that provides approximation guarantees of the solution of the fractional Fokker-Planck equation that we can sample from. Since the assumptions about the density function are quite different in sampling literature and in our paper, it is hard to do one on one comparison of the results. We have provided worst-case guarantees under the much weaker assumption ($\\beta$ times differentiable) and it is not fair to compare them to methods that make stronger assumptions (dissipativity, log-Sobolev, Poincare) to prove convergence results.\n\nOur sample complexity results for approximating the solution of the (fractional) Fokker-Planck Equation in theorems 2 and 5 are better than approximation guarantees in [LCL+21]. On top of that, our method comes with one more advantage over existing machine learning models to approximate the solution of a (fractional) Fokker-Planck equation. Generally, the approximate solution of a PDE obtained by ML models [LCL+21, RPK19, Y+18] is not bound to stay positive and hence can not act as a candidate to approximate the solution of a (fractional) Fokker-Planck equation whose solution is a probability density. In the previous work [RC21], it was shown that the PSD model can be used for the effective representation of probability density. Hence, PSD model-based representation acts as a good candidate to approximate the solution of a (fractional) Fokker-Planck equation. We theoretically show the validity of the PSD model in approximating the solution of the (fractional) Fokker-Planck equation and provide a sample complexity guarantee that was not known before. This guarantee immediately can be used to provide a sample quality guarantee in Corollary 1. Sampling result in such generality was not known before. \n\nWe will take care of the minor comments in the revision. In line 180, $(t-s)^\\alpha$ is correct. We will fix the dimension of Y in line 216.\n\nAlso, please look at our response to the reviewer oceV regarding contributions and experiments. Please look at our response to reviewer biEz for general remarks. \n\n We hope that we have answered all your concerns and that the scores will be raised further.\n"}]

Author rebuttal:

Meta-review: Most reviewers think the paper is interesting.  Please revise the paper and show more benefits of your method over Langevin algorithms.
Decision: Accept (poster)