A central question in statistical mechanics is the determination of the equilibrium distribution [1] from a stochastic di erential equation (SDE). Standard approaches include applying the Monte Carlo method to simulate an SDE or variational inference method that minimizes an entropic regularized energy functional [2,3]. While these approaches have various successes, a less explored route to study the equilibrium is solving the partial di erential equation (PDE) associated with the SDE, which characterizes the evolution of the distribution. Obtaining the equilibrium distribution can be done by looking at the time-independent version of such a PDE. However, for a d-dimensional SDE, one has to solve a d-dimensional PDE, which results in the curse-ofdimensionality. While one can solve the resulting PDE via assuming various low-complexity ansa such as a neural network or tensor-network for the solution, obtaining the solution requires optimizing a non-convex objective function over high-dimensional spaces [4][5][6][7]. In this note, we take a di erent approach to it by using a convexrelaxation strategy to obtain the moments of the distribution. While our method could potentially be used to study the transient behaviour of general SDEs, we demonstrate the success of our method for the case of overdamped Langevin dynamics [8], where we solve for the equilibrium distribution via solving the Fokker-Planck equation (FPE).
In the rest of the paper, we assume we are given a con ning potential function (e.g. [9, De nition 4.2]) V ∶ ℝ d → ℝ that has a pairwise structure, i.e.
(1.1)
The pairwise potential function is ubiquitous in a wide range of physical, chemical systems and molecular dynamics simulations [10][11][12].
The goal is to study the equilibrium distribution ⋆ (x) ∝ exp(-V(x)∕T) of the following overdamped Langevin processes:
where x t ∈ ⊂ ℝ d is the state of stochastic system, = 1∕T where T is the temperature and W t is a d-dimensional Wiener process. Standard approaches (for example, in molecular dynamics) use a long-time Markov chain Monte Carlo (MCMC) simulation to study the equilibrium distribution [13]. For V with a complicated landscape and high dimensionality, the mixing time of an MCMC can be long.
Deterministic approaches include the solution to the following variational problem:
where the minimizer is, in fact, the equilibrium distribution ⋆ . Since ∶ ℝ d → ℝ is a highdimensional function, various low-complexity ansa such as mean-eld [2], Bethe approximation [14], generative models [15] and tensor networks [16] have been used to represent . These approaches lead to either a non-convex optimization domain or approximations to the entropy term [17].
On the other hand, the equilibrium distribution satis es the d-dimensional stationary FPE, L (x) = -1 (x) -∇ ⋅ ( (x)∇V(x)) = 0.
(1.4)
There are many approaches available to solve the FPE. The traditional nite di erence and nite element methods work well for low-dimensional problems, but they both scale exponentially with the number of dimensions. To circumvent the curse of dimensionality, researchers propose to pose various low-complexity ansa on the solution of FPE to control the growth of parameters. For example, [6,7,18] propose to parameterize the unknown PDE solution with deep neural networks and optimize their underlying stochastic di erential equations or variational problems instead. [19][20][21][22] approximate the di erential operators with data-sparse hierarchical matrices. [23] parameterizes the PDE solution using tensor networks [24][25][26], where a similar approach [4] is proposed for the adjoint equation. These parametric models e ectively control the complexity of the problem, but again obtaining them requires the use of non-convex optimization.
royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ...........................................................................................................
In this article, we propose a novel approach to obtain the moments of the equilibrium distribution by solving for the FPE in terms of the low-order marginals or cluster moments, which circumvent the curse of dimensionality. The optimization problems we proposed are fully convex, hence avoiding the di culty of optimizing a globally non-convex objective in neural network and tensor network approaches. While this has similarities to the convex hierarchies proposed to solve (1.3) due to the linearity of the FPE, we completely bypass the need to approximate the entropy term in terms of the moments, which often leads to a nonconvex cost (as in the case of Bethe approximation). We show theoretically that a mean-eld approximation can be recovered when being used to obtain the 1-marginals of the equilibrium density. Hence, our approach can be regarded as a convex-relaxation approach for obtaining a mean-eld solution to the mean-eld FPE (or McKean-Vlasov type equations), which is nonlinear [27].
We view this article as extending methods described in [28] for solving one-dimensional PDE to high-dimensional cases. We show via analysis and numerical experiments that convex hierarchies (cluster moment and marginal relaxations) in similar spirit with cluster expansion [29] or marginal relaxation [30][31][32] can be used to handle high-dimensional systems that are weakly correlated. We note that this work is substantially di erent from [33] where the non-negative density is parameterized as a sum-of-squares polynomial. For us, only the moments of the distribution are solved, which allows us to e ciently characterize a degenerate distribution.
The rest of the note is organized as follows. In §2, we propose various convex hierarchies to solve (1.4). In §3, we extend the proposed method to a time-dependent problem. In §4, we analyse the recovery property of recovering the mean-eld approximation to an equilibrium density from the proposed convex hierarchies. In §5, we study the proposed method in a few numerical experiments. We then conclude in §6.
In this section, we present two formulations: cluster moment relaxation ( §2a) and marginal relaxation ( §2b) for solving the FPE.
Let ⊂ ℝ be the computational domain for a single variable. Let ∶= {t t ∶ d → ℝ} be a nitedimensional function space (we allow the functions to be complex valued to include Fourier basis). Furthermore, let ∶= {q q ∶ d → ℝ, q ≥ 0} be a nite dimensional space of non-negative functions.
Suppose we want to solve
In the remainder of the paper, we will omit ℝ
d and use ∫ to denote integrations of function over ℝ d . While this is a linear PDE in , to discretize , the number of basis elements needs to scale exponentially with respect to d. To overcome the dimensionality of the problem, we solve a relaxed version of the PDE problem via enforcing necessary conditions for L = 0 and ≥ 0 via royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... test functions t ∈ and q ∈ . More speci cally, we let
Here, the inner product is de ned as ⟨q, ⟩ = ∫ d q(x) (x)dx. The adjoint operator L * applied to t ∈ is
The adjoint operator is derived using the fact that we impose the boundary condition
This is a reasonable assumption especially when V is con ning. We now discuss choices of space for and that lead to a tractable convex program.
Let { j ∶ ⊂ ℝ → ℝ} n j=1 be a set of single variable basis with n elements. Let ℱ K be the space of one body to K-body functions where
For our purpose, we choose { j } n j=1 to be the monomial or Fourier basis and 1 ∝ 1. In what follows, we often use the notations x ℐ ∶= (x i 1 , … , x i K ) and (x ℐ ) ∶= j 1 (x i 1 ) … j K (x i K ). We also let the set of sum-of-squares functions be
We choose = ℱ K and furthermore = K . In this case, applying the adjoint operator L * to t ∈ gives
where K ≥ 1. In this way, (2.2) can be wri en in terms of the following matrix variable:
which scales as
, respectively. More precisely, the main optimization problem derived from (2.2) becomes
) royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ...........................................................................................................
where l ℐ (M) = 0 is the reformulation of
in terms of the moment matrix M, and the positive semide nite constraints M ⪰ 0 comes from the fact that
for any A ⪰ 0.
(ii) Symmetric cluster basis
For interacting particle systems, often one needs to deal with a potential V that is symmetric, i.e.
for any in the symmetric group Sym(d). In this case, one can let = ℱ S,K , where ℱ S,K = {
, which consists of symmetric functions in ℱ K . This allows us to span ℱ S,K with n K basis instead of n K n K basis, hence achieving a signi cant complexity gain. Similarly, one can let the set of non-negative functions = S,K where
the symmetric sum-of-squares functions.
Another possible convex relaxation to solve (1.4) is via discretized low-order marginals. This can be useful when V exhibits singularities (for example, when pairwise interactions V ij is coulombic). Let ℐ be the ℐ-th marginal of ,
where ℐ ∈ d K . A reduced order PDE in terms of the marginal can be derived by partial integrations of the FPE:
which gives
This gives equations for all ℐ and { ℐ∪j } d j=1 . We further impose necessary conditions that these marginals are derived from the same density , i.e.
Now to solve (2.11) along with (2.12), we let ̃ ℐ be a n K point discretization of ℐ for each ℐ. Then
2 ℐ (x ℐ ) x 2 i k and ℐ∪j x i k are discretized by a central and forward di erence, respectively, in terms of ̃ ℐ and ̃ ℐ∪j , and we get a discretized version of (2.11), Lℐ ( ̃ ℐ , { ̃ ℐ∪j } j ) = 0. (2.13) Along with a quadrature scheme for enforcing (2.12) one has a Sherali-Adams type linear programming hierarchy [34] with d K+1 n K+1 number of variables (K ≥ 1). For example, when royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... discretizing the domain d as a uniform grid X ⊂ ℝ d , one can use the quadrature scheme
where X ℐ denotes the slice of grid points X in ℐth dimension, and w x ℐ 's are the quadrature weights.
As in [30,31,35], one can further add semide nite constraints to strengthen the convex relaxation,
where G has size
We remark that as in the case of cluster moment relaxation, symmetrization can also be done with marginals, leading to variables with size n K , as in [35].
As we show later in §4, suppose the underlying equilibrium density ⋆ admits a mean-eld or block mean-eld-type solution, i.e.
the proposed convex programs can stably determine each ℐ . However, in general, one cannot hope to recover higher order marginals ℐ∪j or moments M ℐ ,ℐ in (2.6) (moments of ℐ∪ℐ ) since the total number of equations scales as d K = O(d K ) and the number of these higher order marginals or moments scales as
It is possible to introduce penalties to further regularize the solution such that it admits a mean-eld-type behaviour. While we do not pursue these approaches in the current paper, we remark that one can add the following regularizations to the convex programs in §2a,b: (1) Minimize the nuclear norm
for marginal relaxation. For ℐ∩j with marginals being ℐ and j , the maximum entropy solution is ℐ∩j = ℐ j . In a similar vein, interior point methods [36] for convex programming always have a log barrier function, which serves as an entropic-like regularization for ̃ ℐ∩j .
While this is not the main focus of this article, in this section, we discuss possible extensions of our framework to study non-equilibrium dynamics with time-dependent FPE
subject to initial on (x, 0) and extra boundary conditions on (x, ). For example, when studying rst passage time [37,38], one would have the boundary conditions
where particle starts at = 0 at a ∈ ℝ d and gets absorbed when hi ing region A. We let A be a semialgebraic set de ned by
for example, a ball or a half-space
Then the support of can be characterized by {x h(x) ≥ 0}.
To apply convex relaxation to such a problem, we perform a direct discretization in time by approximating (x, t) as 1 (x) ∶= (x, t 1 ), … , m (x) ∶= (x, t m ), where t l = (l-1)T m-1
. We again have an equation
where its convex relaxation can be obtained via testing against t ∈ , which gives ⟨
In this section, we present a di erent way to solve the time-dependent problem via marginal relaxations. We denote the ℐth marginal (ℐ
and its discretization on the ℐth slice X ℐ ⊂ ℐ of X ⊂ d as ̃ l ℐ . We again look at the time discretized equation (3.5); however, we integrate the equation in (3.5) with ∫ ⋅ dx [d]⧵ℐ and solve for l ℐ in terms of the discretized ̃ l ℐ while imposing constraints such as (2.14) on ̃ l ℐ for each time l. As opposed to the discretized stationary equation (2.13), we have
where Lℐ is de ned as in (2.13). As for the initial and boundary conditions in (3.2), we approximate them as
where f 0 ∶ ℝ ℐ → ℝ is a locally supported function that approximates the ℐ dimensional delta function (⋅ -a ℐ ), and a ℐ , A ℐ are the ℐth slice of the coordinate a ∈ d and set A ⊂ d , respectively.
In this section, we analyse the recovery property of the continuous marginal relaxation ( §2b).
We show that the proposed convex program can determine a mean-field approximation to the equilibrium density ⋆ . This is in stark contrast to determining the mean-eld approximation
royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... from a nonlinear Mckean-Vlasov type equation [27,39], since the equation we are solving is completely linear in terms of variables on a convex domain. We believe the proof technique is generalizable to moment-based relaxation ( §2a) and discretized marginal relaxation (2.14) by suitably incorporating the error of approximating low-order marginals by moments or interpolations.
For simplicity, we assume that each variable is on a periodic domain
Since the goal is to determine the marginals of ⋆ ∈ C ∞ ( d ), we further assume that for i, i
To facilitate the discussion, we de ne
and
One can think of F i as some 'e ective force' exerted on the ith dimension. We now make a few assumptions needed for the proof. For the Fokker-Planck operator L, it can be shown that for some function s,
and furthermore,
(see, for example, [40]). We now make assumptions on the bilinear form D (also called the Dirichlet's form).
Assumption 1. The bilinear form D satis es
for any s that is not a constant ⋆ -almost everywhere.
This assumption on the lower bound of the bilinear form is a standard 'spectral gap' assumption for showing the convergence of a time-dependent FPE [41]. The second assumption concerns the existence of a mean-eld approximation. Assumption 2. We assume that the probability density ⋆ ∝ exp(-V) can be approximated by a mean eld, i.e.
.7) royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ...........................................................................................................
This assumption is rather natural since we rst need the equilibrium density to be well approximated by a mean eld, before trying to determine the components of the mean eld computationally. We stress that we do not assume that each pair of potentials exhibits a separable structure. We only assume that, collectively, every pair gives rise to a nearly separable equilibrium distribution. Each individual pair of interacting potentials may well be highly inseparable, as shown in the example in §4c. To some extent, an assumption along the lines of assumption 2 cannot be completely avoided. For the most general potentials, it is unlikely that one will ever be able to recover the 1-marginal of the equilibrium distribution based on the lowest-level convex hierarchy. This is not unlike the manner the McKean-Vlasov equation approximates the FPE; the McKean-Vlasov's solution cannot be expected to approximate the ground truth 1-marginal distribution.
The goal is to show the solution of the following convex problem:
recovers the mean elds
Here, the set is de ned to be
By abuse of notation, in this section, we treat a single variable function s i acting on the ith dimensional variable x i as a C ∞ ( ) and C ∞ ( d ) function interchangeably. In the la er case, we think of s i (x i ), x i ∈ as s i (x) = s i (x i ), x ∈ d . Although it seems like the third set of constraints has an in nite number of constraints, it is equivalent to the positive semide nite constraints in (2.15) (with a suitable discretization using basis in ) and therefore can be enforced with semide nite programming.
In order to prove our main theorem, we make another assumption, speci cally pertaining to the convex relaxation (4.8). This assumption is relatively unusual in the literature of the FPE. Assumption 3. We assume that F i de ned in (4.3) satis es
for any
i≠i satisfying the constraints in (4.8).
This assumption essentially says the 'e ective force' computed from any two marginals
⧵i is not that di erent from treating each 2-marginal as independent, i.e. ii = i i . This depends on the nature of V and also on the constraints in (4.8). In some sense, one can think of 2 as a 'relaxation error'. In §4c, we give a non-trivial example with strong pairwise interactions to show how this assumption can hold without separability.
We are now ready to present the main theorem.
royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... Theorem 1. Under assumptions 1, 2, 3, the solution to the marginal relaxation problem (4.8) satisfies √
)
⧵i satisfying the constraints in (4.8), we have
.
(4.12)
The rst equality follows from the rst constraints of (4.8). The second equality follows from integration by parts and the de nitions in (4.2) and (4.3). The fourth and fth equality is due to the property of the Fokker-Planck operator (4.4) and (4.5). The sixth equality is based on le ing ∏d
The rst inequality is due to Cauchy-Schwarz, and the last inequality is due to assumption 3. Now, let
and from (4.12), we get
.
(4.15)
At rst sight, it seems like the only pairwise potential that satis es assumption 3 with a small enough 2 is that each
.e. a separable potential. However, in this section, we argue that it is possible to have non-separable
and further
where each a ii ∈ ℝ n×n , and { j } n j=1 ⊂ L 2 ( ⋆ ) is a set of orthonormal single variable basis. We let a ii (k i , k i ) ∼ (0, 2 ). In what follows, we show that one can have as large o( √ d) and still 1), giving rise to a meaningful upper-bound in theorem 1. In this case, the term ii completely dominates E i , E i when d is large. We remark that this choice of potential is motivated by the ubiquitous Rastrigin function [42] for benchmarking global optimization algorithms, where a separable quadratic potential well is buried in oscillatory functions. It is also motivated by a class of well-known optimization problems in theoretical computer science known as spike model/planted models (e.g. planted clique [43], stochastic block model [44] and synchronization problem [45]), where a clean signal to be recovered is planted in random noise. For these problems, it has been shown that with a well-designed algorithm like semide nite programming [46][47][48]), one can provably obtain the global minimizer in polynomial time. We therefore want to demonstrate that our adaptation of semide nite programming to a PDE se ing retains its ability on the most common potentials in global optimization.
To calculate 2 in assumption 3 for this se ing, we rst notice that
where G ii ∈ ℝ n×n , ti ∈ ℝ n×n are formed by discretizing with { i } n i=1 . In particular,
.19) One can think of G = [G ii ] i,i ∈[d] as a size nd × nd matrix. It has several properties, in particular, G ⪰ 0 due to the constraint in (4.8). Furthermore, tr(G ii ) ≤ n. We can factorize G = YY T where Y = [Y T 1 , … , Y T d ] T and each Y i ∈ ℝ n×n , and further
≤ constant (equation (4.9)). The second inequality follows from standard spectral norm bounds on the nd × nd random Gaussian matrix [49], and the third inequality follows from ‖Y i ‖ F ≤ √ n. In the following, we treat n as a constant. Similarly, by replacing Y in the above derivation by
we can also derive
This gives an upper bound for (4.17), and by le ing
We now compare 2 with the eigenvalue 1 in assumption 1. Let and and be the second lowest eigenvalue and its corresponding right and left eigenfunctions for the following equation:
where MF ∝ exp(-E). Fixing = (d -1) -1 , the operator L becomes
This choice of = (d -1) -1 such that
Alternatively, one can start with a ∑d i≠j=1 V ij that is re-scaled in the aforementioned way and let = O(1). These choices ensure that the transport term arising from pairwise interacting potential scales in a similar way as the di usion term in the FPE.
We now want to show this operator L is not that di erent from (d -1)
The second smallest eigenvalue of the rst term (d -1)
, stemming from the fact that
by choosing 1 (x) = (x 1 ). We now want to study the e ect of perturbation on the eigenvalue caused by ii 's. Let 1 (x) = (x 1 ). First order perturbation theory [50] says the eigenvalue of L,
.26) royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... and further ∑ i∈[d]
where the rst inequality follows from Var ⟨
due to the fact that ⟨
, ⟩ is a sub-Gaussian random variable [49]. The second inequality follows from (4.9) and 1 = . By rst order perturbation theory of the eigenvalue, we get 1 ≈ (d -1) + O( √ d). In this case,
for = o( √ d), leading to a meaningful upper-bound for theorem 1.
A similar perturbation analysis based on rst-order perturbation theory gives eigenfunction perturbation 1 in assumption 2 where
In this subsection, we demonstrate the numerical performance of solving the stationary FPE via cluster moment relaxation. We consider two systems in this subsection, one with the separable double-well potential, which is intrinsically a one-dimensional potential so we can easily visualize, and the other one with the Ginzburg-Landau potential, for which we can only try to compare the moments or marginals since the true density is exponentially sized.
We consider the following double-well potential function:
which essentially only contains one-body terms and can be easily separable. The equilibrium distribution for this is simply
where Z is the partition function. The separable structure gives us a convenient way to visualize the distribution via the 1-marginals. In this example, we take = 2.5 and a d = 10 dimensional system. We use monomials as test functions for this example. More speci cally, we let the single variable basis set to be monomials, i.e. { j (⋅) = (⋅) j-1 } 9 j=1 , and we take the test functional space to be the induced 1-body functional space = ℱ 1 . The potential function V can be exactly expanded by monomials up to order 4. For convenience, we denote the moment matrix M in (2.8) as second-order moments with notation M 2 and introduce the rst-order moment vector
) to be a subset of M 2 . We solve the feasibility semide nite programming problem (2.8) with CVX [51,52]. We compare the estimated moments M1 and M2 from the proposed semide nite relaxation with the ground truth moments. We measure the relative errors in rst-and second-order moments de ned by
respectively. In our example, we achieve E 1 = 5.30 × 10 -3 and E 2 = 7.84 × 10 -3 error for the rstand second-order moments. Next, we can estimate the 1-marginals from the estimated moments by solving the following maximum entropy problem:
where i is the i-th marginal, which is a standard method for density estimation from moments [53]. We show the estimated 1-marginal by solving a discretized program of (5.4) and the true 1-marginal in gure 1.
The Ginzburg-Landau theory was developed to provide a phenomenological description of many-body systems [54]. In this numerical example, we consider a periodic discretized Ginzburg-Landau model,
) royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... ρ 5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 True 1-marginal Reconstruct 1-marginal -3 -2 -1 0 1 2 3 x 5 where h = 1∕(d + 1) and U d+1 = U 1 , U 0 = U d based on periodic condition. We x d = 10, = 0.03 and the temperature = 1∕16. We let the computational domain be U ∈ [-1.6, 1.6] d .
In this example, we demonstrate the performance of the proposed relaxation with the complex Fourier basis. We note that the performance using the monomial basis is similar. More speci cally, we take { j (⋅) = exp(i 3.2 ⋅) j } 10 j=-10 as univariate basis functions for each dimension which induce the test function space = ℱ 1 . The potential function V can be expanded with the same set of 21 complex Fourier basis functions in each dimension with relative error 5.31 × 10 -3 .
Using the test function space discussed above, we achieve E 1 = 1.55 × 10 -2 and E 2 = 5.11 × 10 -1 relative error in rst-and second-order moments. The proposed method can identify the rstorder moment accurately, but the estimated second-order moment has a relatively large error. We may need to use third-order moments if we want to identify the second-order moments with high accuracy due to the deviation of the Ginzburg-Landau model from a mean-eld model.
Similarly, we can solve the maximum entropy problem (5.4) to reconstruct the estimated marginals from the moments for any dimension. Here, we compare the reconstructed marginal (red solid line) and ground truth marginal (blue dashed line) in gure 2.
When having a strongly interacting system, it is possible that assumptions 2 and 3 can fail, rendering the analysis provided in §4 not applicable. In this section, we therefore investigate numerically the proposed method when having a strong repulsive potential. More precisely, we study d particles in one dimension interacting through pairwise soft coulombic potentials, with a quartic con nement potential:
4 . (5.6) It is tricky to discretize such a potential with a global polynomial basis set due to its singular behaviour. Therefore, in what follows, instead of using the cluster moment relaxation, we resort royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ........................................................................................................... to the marginal relaxation in §2b, where local interpolants can be used to discretize the resulting equations.
Notably, the equilibrium distribution for this potential has a permutation symmetry, which gives rise to a further computational cost reduction since all K-th marginals are the same (K ∈ [d]). More precisely, when working with (2.11) where K = 1, i.e. the equations that determine the 1-marginals, one can solve a single equation
where we use the fact that 1 and 12 are the same as all other one marginals and two marginals, respectively. Although, at this point, one can further discretize (2.11) by a nite di erence scheme and solve the resulting nite-dimensional convex program, a further simpli cation can be made to (5.7) that allows us to discretize the equation with fewer quadrature points. In particular, we can integrate (5.7) once more,
which results in
where we use a Dirichlet's boundary condition to remove the boundary contributions. In this case, we only need to discretize the rst derivative of the density, unlike in (5.7). This can reduce the number of points needed to discretize the equation as the equilibrium distribution becomes steep due to the repulsive potential. In (5.9), we discretize 1 and 12 as a vector ̃ 1 and a matrix ̃ 12 , respectively. Then, we discretize 1
x 1 by a central di erence and the integral in (5.9) by a trapezoidal rule. We now need to impose constraints on ̃ 1 and ̃ 12 as in (2.14), which are necessary conditions for ̃ 1 and ̃ 12 being the (discretized) marginals of a d-dimensional distribution. As mentioned previously in §2b, one can leverage the fact that these are the marginals of a symmetric distribution and obtain a compact set of constraints as in [35],
where h is the length between two grid points and 1 is the all-one vector.
In this experiment, we report the results for d = 4, i.e. having four particles, and we discretize each dimension with n = 100 points between [-4.5, 4.5]. We enforce a Dirichlet boundary condition as well, as the solution decays to zero on the boundary due to the quartic potential in (5.6). In this case, we can directly compare our solution with the discretized ground truth. We compute the 1-marginals with = 4 and = 8, and the relative 1-marginal error is 1.4 × 10 -1 and 2.1 × 10 -1 , respectively. Visual comparisons are provided in gure 3. The error for the coulombic case is larger than the error in previous experiments. Such error is probably inevitable if we insist on using the lowest level in the proposed convex hierarchy to deal with a strongly interacting system. The lowest-level convex relaxation with only one and two marginals has commi ed a model error to the FPE. This is analogous to how the McKean-Vlasov equation is already an approximation of the FPE, and one will not expect the McKean-Vlasov's solution to always resemble the ground truth 1-marginal distribution. We expect to see improvements when we use a higher cluster marginals to perform the convex relaxation.
In this subsection, we demonstrate the numerical performance of solving the time-dependent FPE via marginal relaxation. More speci cally, we consider the extension of our framework to the time-dependent FPE as discussed in §3. We consider a double-well system (5.1) with d = 2 for be er visualization purposes. For all spatial dimensions, we restrict the support of our marginals in hypercube [-2, 2] 2 and use 100 evenly distributed grid points to discretize each dimension. This means we consider the density to be negligible outside the hypercube [-2, 2] 2 , which is a property of systems with con ning potential. The time dimension is discretized by 100 grid points as well with a gap t = 2 × 10 -4 . We solve the time-dependent optimization problem (3.8) with marginals of order at most 2.
We use the initial and boundary conditions in (3.2). We assume the particles are concentrated around (-1, 0, 0, 0, 0), one of the local minima of the double-well potential, at the beginning of evolution. More speci cally, we let (x, 0) = 1 (x 1 , 0) ∏d i=2 i (x i , 0), where 1 (x 1 , 0) = Uniform[-1, 2, 0.8], i (x i , 0) ∝ exp(-0.3 x i 2 ), i ≠ 1. We test two scenarios for absorbing states: (1) there is an absorbing state at 0.2 and we require (0.2, ) = 0 for all ∈ [0, 0.02]; and (2) there is an absorbing state at 1.0, and we require (1.0, ) = 0 for all ∈ [0, 0.02]. Since the double-well potential is separable, the ground truth 1-marginals can be obtained by solving the time-dependent FPE for the rst dimension only,
where x a = 0.2, 1.0 corresponds to scenarios (1) and (2), respectively. We compare the estimated 1-marginals with the ground truth marginals by solving (5.11). The average elementwise relative error for all 1-marginals over time is (1) 1.02 × 10 -1 and (2) 1.39 × 10 -1 for the two scenarios mentioned above. We further visualize the evolution of 1-marginals in gures 4 and 5.
Although traditional numerical methods such as nite elements have been wildly successful in solving low-dimensional PDEs to high accuracy, they are prohibitively expensive to use in a high-dimensional se ing due to their exponential scaling with respect to dimensionality. In this article, we propose a novel convex optimization method to solve the high-dimensional Fokker-Planck PDE by solving for the low-order marginals/moments of the distribution. The resulting optimization problem is fully deterministic and convex, and can be solved with existing convex program solvers. By solving the low-order marginals/moments instead, we reduce the total number of unknown parameters to estimate and, as a result, break the curse of dimensionality. The complexity of the algorithm depends on the number of dimensions and the number of marginals/moments to solve for. Moreover, our algorithm can potentially provide a coarse description of the equilibrium distribution by characterizing its low-order statistics and can be used to initialize other local optimization algorithms (e.g. MCMC or tensor networks) for further re nement. Finally, we also discuss the generalization to the time-dependent FPE, which could be used to study non-equilibrium physics.
At this point, we have been using 1-cluster moment and 1-marginal relaxation to perform theoretical analyses and numerical experiments. A natural future direction is to investigate how the use of K > 1 cluster moments/marginals relaxation can lead to a be er performance. Since we approximate the FPE with low-order cluster moments/cluster marginals, we cannot expect the same type of accuracy as traditional numerical PDE methods, which use exponentially large computational resources. However, it is important to understand for which distributions one can approximate the FPE with low-order cluster moments/cluster marginals. While we discuss the possibility of using the proposed method for the time-dependent FPE, we have only investigated its use in relatively simple se ings. All these aspects deserve detailed investigations and individual treatments in future works.
royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 481: 20240001 ...........................................................................................................