Abstract
We propose the novel numerical scheme for solution of the multidimensional Fokker–Planck equation, which is based on the Chebyshev interpolation and the spectral differentiation techniques as well as low rank tensor approximations, namely, the tensor train decomposition and the multidimensional cross approximation method, which in combination makes it possible to drastically reduce the number of degrees of freedom required to maintain accuracy as dimensionality increases. We demonstrate the effectiveness of the proposed approach on a number of multidimensional problems, including Ornstein-Uhlenbeck process and the dumbbell model. The developed computationally efficient solver can be used in a wide range of practically significant problems, including density estimation in machine learning applications.
1 Introduction
Fokker–Planck equation (FPE) is an important in studying properties of the dynamical systems, and has attracted a lot of attention in different fields. In recent years, FPE has become widespread in the machine learning community in the context of the important problems of density estimation () for neural ordinary differential equation (ODE) (; ), generative models (), etc.
Consider a stochastic dynamical system which is described by stochastic differential equation (SDE) of the form1where is a q-dimensional space-time white noise, f is a known d-dimensional vector-function and , are known matrices. The FPE for the corresponding probability density function (PDF) of the spatial variable has the formwhere is a diffusion tensor.
One of the major complications in solution of the FPE is the high dimensionality of the practically significant computational problems. Complexity of using grid-based representation of the solution grows exponentially with d, thus some low-parametric representations are required. One of the promising directions is the usage of low-rank tensor methods, studied in (). The equation is discretized on a tensor-product grid, such that the solution is represented as a d-dimensional tensor, and this tensor is approximated in the low-rank tensor train format (TT-format) (). Even with such complexity reduction, the computations often take a long time. In this paper we propose another approach of using low-rank tensor methods for the solution of the FPE, based on its intimate connection to the dynamical systems.
The key idea can be illustrated for , i.e. in the deterministic case. For this case the evolution of the PDF along the trajectory is given by the formulawhere is a trace operation for the matrix. Hence, to compute the value of at the specific point , it is sufficient to find a preimage such that if it is used as an initial condition for eq. 1, then we arrive to . To find the preimage, we need to integrate the eq. 1 backwards in time, and then to find the PDF value, we integrate a system of eqs 1, 3. Since we can evaluate the value of at any , we can use the cross approximation method (CAM) (; ; ) in the TT-format to recover a supposedly low-rank tensor from its samples. In this way we do not need to have any compact representation of f, but only numerically solve the corresponding ODE. For the situation is more complicated, but we develop a splitting and multidimensional interpolation schemes that allow us effectively recompute the values of the density from some time moment t to the next step .
To summarize, main contributions of our paper are the following:
• we derive a formula to recompute the values of the PDF on each time step, using the second order operator splitting, Chebyshev interpolation and spectral differentiation techniques;
• we propose to use a TT-format and CAM to approximate the solution of the FPE which makes it possible to drastically reduce the number of degrees of freedom required to maintain accuracy as dimensionality increases;
• we implement FPE solver, based on the proposed approach, as a publicly available python code2, and we test our approach on several examples, including multidimensional Ornstein-Uhlenbeck process and dumbbell model, which demonstrate its efficiency and robustness.
2 Computation of the Probability Density Function
For ease of demonstration of the proposed approach, we suppose that the noise has the same dimension as the spatial variable (), and the matrices in eq. 1 and eq. 2 have the form3where is a scalar diffusion coefficient. Then eqs 1, 2 can be rewritten in a more compact formwhere d-dimensional spatial variable has the corresponding PDF with initial conditionsTo construct the PDF at some moment τ () for the known initial distribution , we discretize eqs 5, 6 on the uniform time grid with M () pointsand introduce the notation for value of the spatial variable at the moment and for values of the PDF at the same moment.
2.1 Splitting Scheme
Let and be diffusion and convection operators from the eq. 6then on each time step m () we can integrate equationon the interval , to find for the known value from the previous time step. Its solution can be represented in the form of the product of an initial solution with the matrix exponentialand if we apply the standard second order operator splitting technique (), thenwhich is equivalent to the sequential solution of the following equationswith the final approximation of the solution .
2.2 Interpolation of the Solution
To efficiently solve the convection eq. 14, we need the ability to calculate the solution of the diffusion eq. 13 at arbitrary spatial points, hence the natural choice for the discretization in the spatial domain are Chebyshev nodes, which makes it possible to interpolate the corresponding function on each time step by the Chebyshev polynomials ().
We introduce the d-dimensional spatial grid as a tensor product of the one-dimensional grids4where () is a number of points along the kth spatial axis (), and the total number of the grid points is . Note that this grid can be also represented in the flatten form as a following matrixwhere , and by we denoted an operation of construction of the multi-index from the flatten long index according to the big-endian conventionSuppose that we calculated PDF on some time step m () at the nodes of the spatial grid [note that for the case , the corresponding values come from the known initial condition ]. These values can be collected as elements of a tensor5 such thatwhere ().
Let us interpolate PDF via the system of orthogonal Chebyshev polynomials of the first kindin the form of the naturally cropped sumwhere is some spatial point and interpolation coefficients are elements of the tensor . For construction of this tensor we should set equality in the interpolation nodes eq. 16for all combinations of ().
Therefore the interpolation process can be represented as a transformation of the tensor to the tensor according to the system of eq. 22. If the Chebyshev polynomials and nodes are used for interpolation, then a good way is to apply a fast Fourier transform (FFT) () for this transformation. However the exponential growth of computational complexity and memory consumption with the growth of the number of spatial dimensions makes it impossible to calculate and store related tensors for the multidimensional case in the dense data format. Hence in the next sections we present an efficient algorithm for construction of the tensor in the low-rank TT-format.
2.3 Solution of the Diffusion Equation
To solve the diffusion eqs 13, 15 on the Chebyshev grid, we discretize Laplace operator using the second order Chebyshev differential matrices [see, for example, ()] such that , where for each spatial dimension with if or and otherwise, and one dimensional grid points defined from eq. 16. Then discretized Laplace operator has the form6Let be the known initial condition for the diffusion equation on the time step m (), then for the solution at the moment we havewhere an operation constructs a vector from the tensor by a standard reshaping procedure like eq. 18. And finally due to the well known property of the matrix exponential, we come toIf we can represent the initial condition in the form of Kronecker product of the one-dimensional tensors (for example, in terms of the TT-format in the form of the Kronecker products of the TT-cores, as will be presented below in this work), then we can efficiently evaluate the formula eq. 26 to obtain the desired approximation for solution .
2.4 Solution of the Convection Equation
Convection eq. 14 can be reformulated in terms of the FPE without diffusion part, when the corresponding ODE has the formIf we consider the differentiation along the trajectory of the particles, as was briefly described in the Introduction, thenwhere we replaced the term by the right hand side of eq. 14 and by the right hand side of the corresponding equation in eq. 27.
Hence equation for w may be rewritten in terms of the trajectory integration of the following systemLet us integrate eq. 29 on a time step m (). If we set any spatial grid point () as initial condition for the spatial variable, then we’ll obtain solution for some point outside the grid (see Figure 1 with the illustration for the two-dimensional case). Hence we should firstly solve eq. 27 backward in time to find the corresponding spatial point that will be transformed to the grid point by the step . If we select this point and the related value as initial conditions for the system eq. 29, then its solution will be related to the point of interest .
FIGURE 1
Note that, according to our splitting scheme, we solve the convection part eq. 14 after the corresponding diffusion eq. 13, and hence the initial condition is already known and defined as a tensor on the Chebyshev spatial grid. Using this tensor, we can perform interpolation according to the formula eq. 22 and calculate the tensor of interpolation coefficients . Then we can evaluate the approximated value at the point as according to eq. 21.
Hence our solution strategy for convection equation is the following. For the given spatial grid point we integrate equationbackward in time to find the corresponding point . Then we find the value of w at this point, using interpolation , and then we solve the system eq. 29 on the time interval with initial condition to obtain the value at the point . The described process should be repeated for each grid point () and, ultimately, we’ll obtain a tensor which is the approximated solution of convection part eq. 14 of the splitting scheme on the Chebyshev spatial grid.
An important contribution of this paper is an indication of the possibility and a practical implementation of the usage of the multidimensional CAM in the TT-format to recover a supposedly low-rank tensor from computations on only a part of specially selected spatial grid points. This scheme will be described in more details later in the work after setting out the fundamentals of the TT-format.
3 Low-Rank Representation
There has been much interest lately in the development of data-sparse tensor formats for high-dimensional problems. A very promising tensor format is provided by the tensor train (TT) approach (; ), which was proposed for compact representation and approximation of high-dimensional tensors. It can be computed via standard decompositions (such as SVD and QR-decomposition) but does not suffer from the curse of dimensionality7.
In many analytical considerations and practical cases a tensor is given implicitly by a procedure enabling us to compute any of its elements, so the tensor appears rather as a black box. For example, to construct the convection part of PDF (i.e., the tensor introduced above), we should compute the corresponding function for all possible sets of indices. This process requires an extremely large number of operations and can be time-consuming, so it may be useful to find some suitable low-parametric approximation of this tensor using only a small portion of all tensor elements. CAM () which is a widely used method for approximation of high-dimensional tensors looks appropriate for this case.
In this section we describe the properties of the TT-format and multidimensional CAM that are necessary for efficient solution of our problem, as well as the specific features of the practical implementation of interpolation by the Chebyshev polynomials in terms of the TT-format and CAM.
3.1 Tensor Train Format
A tensor is said to be in the TT-format (), if its elements are represented by the formulawhere (), three-dimensional tensors are named TT-cores, and integers (with convention ) are named TT-ranks. The latter formula can be also rewritten in a more compact formwhere is an matrix for each fixed (since , the result of matrix multiplications in eq. 32 is a scalar). And a vector form of the TT-decomposition looks likewhere the slices of the TT-cores are vectors of length ().
The benefit of the TT-decomposition is the following. Storage of the TT-cores requires less or equal than memory cells (instead of cells for the uncompressed tensor, where is an average size of the tensor modes), and hence the TT-decomposition is free from the curse of dimensionality if the TT-ranks are bounded.
The detailed description of the TT-format and linear algebra operations in terms of this format8 is given in works (; ). It is important to note that for a given tensor in the full format, the TT-decomposition (compression) can be performed by a stable TT-SVD algorithm. This algorithm constructs an approximation in the TT-format to the given tensor with a prescribed accuracy in the Frobenius norm9but a procedure of the tensor approximation in the full format is too costly, and is even impossible for large dimensions due to the curse of dimensionality. Therefore more efficient algorithms like CAM are needed to quickly construct the tensor in the low rank TT-format.
3.2 Cross Approximation Method
The CAM allows to construct a TT-approximation of the tensor with prescribed accuracy , using only part of the full tensor elements. This method is a multi-dimensional analogue of the simple cross approximation method for the matrices () that allows one to approximate large matrices in time by computing only elements, where is an average size of the matrix modes and R is the rank of the matrix. The CAM and the TT-format can significantly speed up the computation and reduce the amount of consumed memory as will be illustrated in the next sections on the solution of the model equations.
The CAM constructs a TT-approximation to the tensor , given as a function , that returns the th entry of for a given set of indices. This method requires only
3.3 Multidimensional Interpolation
As was discussed in the previous sections, we discretize the FPE on the multidimensional Chebyshev grid and interpolate solution of the first diffusion equation in the splitting scheme eq. 13 by the Chebyshev polynomials to obtain its values on custom spatial points (different from the grid nodes) and then perform efficient trajectory integration of the convection eq. 14.
The desired interpolation may be constructed from solution of the system of eq. 22 in terms of the FFT (), but for the high dimension numbers we have the exponential growth of computational complexity and memory consumption, hence it is very promising to construct tensor of the nodal values and the corresponding interpolation coefficients in the TT-format.
Consider a TT-tensor with the list of TT-cores , which collects PDF values on the nodes of the Chebyshev grid at some time step (the related function is , and this tensor is obtained, for example, by the CAM or according to TT-SVD procedure from the tensor in the full format). Then the corresponding TT-tensor of interpolation coefficients with the TT-cores can be constructed according to the scheme, which is presented in Algorithm 2 [we’ll call it as a function below].
In this Algorithm we use standard linear algebra operations and , which rearrange the axes and change the dimension of the given tensor respectively, function for construction of the one-dimensional FFT for the given vector, and function from the ttpy package, which round the given tensor to the prescribed accuracy ϵ. Note that the inner loop in Algorithm 2 for may be replaced by the vectorized computations of the corresponding two-dimensional FFT.
For the known tensor we can perform a fast computation of the function value at any given spatial point by a matrix product of the convolutions of the TT-cores of with appropriate column vectors of Chebyshev polynomialsWe’ll call the corresponding function as below. This function constructs a list of values for the given set of I points (), using interpolation coefficients and sequentially applying the formula eq. 35 for each spatial point.
4 Detailed Algorithms
In Algorithms 3, 4 and 5 we combine the theoretical details discussed in the previous sections of this work and present the final calculation scheme for solution of the multidimensional FPE in the TT-format, using CAM (function , Algorithm 1)10 and interpolation by the Chebyshev polynomials (function from Algorithm 2 that constructs interpolation coefficients and function that evaluates interpolation result at given points according to the formula eq. 35).
We denote by the standard linear algebra operation that evaluates the Einstein summation convention on the operands (see, for example, the numpy python package). Function stack arrays in sequence vertically, function (where and are initial and final times, is the right hand side of equations, and matrix collects initial conditions) solves a system of ODE with vectorized initial condition by the one step of the fourth order Runge-Kutta method.
5 Numerical Examples
In this section we illustrate the proposed computational scheme, which was presented above, with the numerical experiments. All calculations were carried out in the Google Colab cloud interface11 with the standard configuration (without GPU support).
Firstly we consider an equation with a linear convection term—Ornstein-Uhlenbeck process (OUP) () in one, three and five dimensions. For the one-dimensional case, which is presented for convention, we only solve equation using the dense format (not TT-format), hence the corresponding results are used to verify the general correctness and convergence properties of the proposed algorithm, but not its efficiency. In the case of the multivariate problems we use the proposed tensor based solver, which operates in accordance with the algorithm described above. To check the results of our computations, we use the known analytic stationary solution for the OUP, and for the one-dimensional case we also perform comparison with constructed analytic solution at any time moment.
Then we consider more complicated dumbbell problem () which may be represented as a three-dimensional FPE with a nonlinear convection term. For this case we consider the Kramer expression and compare our computation results with the results from another works for the same problem.
In the numerical experiments we consider the spatial region such that PDF is almost vanish on the boundaries , and the initial condition is selected in the form of the Gaussian functionwhere parameter s is selected as . To estimate the accuracy of the obtained PDF () we use the relative 2-norm of deviation from the exact value ()We compute the value through the given function, using a CAM with an accuracy parameter two orders of magnitude higher, than the one that was used in the solver.
5.1 Numerical Solution of the Ornstein-Uhlenbeck Process
Consider FPE of the form
eq. 6in the
d-dimensional case with
where
is invertible real matrix,
is the long-term mean,
and
(
),
(
). This equation is a well known multivariate OUP with the following properties [see for example (
;
)]:
• mean vector is
• covariance matrix is
and, in our case as noted above
;
• transitional PDF is
• stationary solution is
where matrix
can be found from the following equation
• the (multivariate) OUP at any time is a (multivariate) normal random variable;
• the OUP is mean-reverting (the solution tends to its long-term mean as time t tends to infinity) if all eigenvalues of A are positive (if in the one-dimensional case).
5.1.1 One-Dimensional Process
Let consider the one-dimensional () OUP withWe can calculate the analytic solution in terms of only spatial variable and time via integration of the transitional PDF eq. 41Accurate computations lead to the following formulawhere is defined by eq. 40 and for the one-dimensional case may be represented in the formUsing the formulas eq. 42 and eq. 43 we can represent a stationary solution for the one-dimensional case in the explicit formWe perform computation for spatial points and time points and compare the numerical solution with the known analytic eq. 46 and stationary eq. 48 solution. In the Figure 2 we present the corresponding result. Over time, the error of the numerical solution relative to the analytical solution first increases slightly, and then stabilizes at approximately . At the same time, the numerical solution approaches the stationary one, and the corresponding error at large times also becomes approximately . Note that the time to build the solution was about 5 s.
FIGURE 2
5.1.2 Three-Dimensional Process
Our next example is the three-dimensional () OUP with the following parametersWhen carrying out numerical calculation, we select as the accuracy of the CAM, 100 as a total number of time points and 30 as a number of points along each of the spatial dimensions. The computation result is compared with the stationary solution eq. 42 which was obtained as solution of the related matrix eq. 43 by a standard solver for Lyapunov equation.
The result is shown in Figure 3. As can be seen, the TT-rank12 remains limited, and the accuracy of the solution over time grows, reaching by the time . The time to build the solution was about 26 s.
FIGURE 3
To evaluate the efficiency of the proposed algorithm in the TT-format, we also solve these three-dimensional OUP, using dense format (as for the one-dimensional case, all arrays are presented in its full form). The corresponding calculation took about 376 s, so in this case we have an acceleration of calculations by more than an order of magnitude.
5.1.3 Five-Dimensional Process
This multidimensional case is considered in the same manner as the previous one. We select the following parametersWe select the same values as in the previous example for the CAM accuracy (), the number of time points (100) and the number of spatial points (30), and compare result of the computation with the stationary solution from eq. 42 and eq. 43.
The results are presented on the plots on Figure 4. The TT-rank of the solution remains limited and reaches the value 4.5 at the end time step, and the solution accuracy reaches almost . The time to build the solution was about 100 s.
FIGURE 4
5.2 Numerical Solution of the Dumbbell Problem
Now consider a more complex non-linear example corresponding to the three-dimensional (
) dumbbell model of the form
eq. 6with
13where
Making simple calculations (taking into account the specific form of the matrix
A), we get explicit expressions for the function and the required partial derivatives (
)
Next, we consider the Kramer expression
and as the values of interest (as in the works (
;
)) we select
During the calculations we used the following solver parameters:
• the accuracy of the CAM is ;
• the number of time grid points is 100;
• the number of grid points along each of the spatial dimensions is 60.
The results are presented on the plots on Figure 5. The time to build the solution was about 200 s (also additional time was required to calculate the values and from eq. 56 and eq. 57 respectively). As can be seen, the TT-rank remains limited, and its stationary value is about 8. We compared the obtained stationary values of the and variables:with the corresponding results from ()14, and we get the following values for relative errors
FIGURE 5
6 Related Works
The problem of uncertainty propagation through nonlinear dynamical systems subject to stochastic excitation is given by the FPE, which describes the evolution of the PDF, and has been extensively studied in the literature. A number of numerical methods such as the path integral technique (; ), the finite difference and the finite element method (; ) have been proposed to solve the FPE.
These methods inevitably require mesh or associated transformations, which increase the amount of computation. The problem becomes worse when the system dimension increases. To maintain accuracy in traditional discretization based numerical methods, the number of degrees of freedom of the approximation, i.e. the number of unknowns, grows exponentially as the dimensionality of the underlying state-space increases.
On the other hand, the Monte Carlo method, that is common for such kind of problems (; ), has slow rate of convergence, causing it to become computationally burdensome as the underlying dimensionality increases. Hence, the so-called curse of dimensionality fundamentally limits the use of the FPE for uncertainty quantification in high dimensional systems.
In recent years, low-rank tensor approximations have become especially popular for solving multidimensional problems in various fields of knowledge (). However, for the FPE, this approach is not yet widely used. We note the works (; ; ; ; ) in which the low-rank TT-decomposition was proposed for solution of the multidimensional FPE. In these works, the differential operator and the right-hand side of the system are represented in the form of TT-tensor. Moreover, in paper () the joint discretization of the solution in space-time is considered. The difference of our approach from these works is its more explicit iterative form for time integration, as well as the absence of the need to represent the right hand side of the system in a low-rank format, which allows to use this approach in machine learning applications.
7 Conclusion
In this paper we proposed the novel numerical scheme for solution of the multidimensional Fokker–Planck equation, which is based on the Chebyshev interpolation and spectral differentiation techniques as well as low rank tensor approximations, namely, the tensor train decomposition and cross approximation method, which in combination make it possible to drastically reduce the number of degrees of freedom required to maintain accuracy as dimensionality increases.
The proposed approach can be used for the numerical analysis of uncertainty propagation through nonlinear dynamical systems subject to stochastic excitations, and we demonstrated its effectiveness on a number of multidimensional problems, including Ornstein-Uhlenbeck process and dumbbell model.
As part of the further development of this work, we plan to conduct more rigorous estimates of the convergence of the proposed scheme, as well as formulate a set of heuristics for the optimal choice of number of time and spatial grid points and tensor train rank. Another promising direction for further research is the application of established approaches and developed solver to the problem of density estimation for machine learning models.
Statements
Data availability statement
Program code, input data and calculation results can be found here: https://github.com/AndreiChertkov/fpcross
Author contributions
IO contributed to general formulation of the problem and formulation of the preliminary version of the algorithm. AC contributed to development of the final version of algorithms and program code, carrying out numerical experiments.
Funding
The work was supported by Ministry of Science and Higher Education grant No. 075-10-2021-068.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Footnotes
1.^Vectors and matrices are denoted hereinafter by lower case bold letters () and upper case letters () respectively. We denote the th element of an matrix A as and assume that , . For vectors we use the same notation: is the i-th element of the vector (). In addition, for a compact representation of an i-th (, where ) element of a vector function , we will use the notation , which means .
2.^The code is publicly available from https://github.com/AndreiChertkov/fpcross.
3.^We use notation for the () identity matrix.
4.^We suppose that for each spatial dimension the variable x varies within . In other cases, an appropriate scaling can be easily applied.
5.^By tensors we mean multidimensional arrays with a number of dimensions d (). A two-dimensional tensor () is a matrix, and when it is a vector. For tensors with we use upper case calligraphic letters (). The th entry of a d-dimensional tensor is denoted by , where () and is a size of the k-th mode, and mode-k slice of such tensor is denoted by .
6.^Note that for the case , we have only one matrix which greatly simplifies the computation process.
7.^By the full format tensor representation or uncompressed tensor we mean the case, when one calculates and saves in the memory all tensor elements. The number of elements of an uncompressed tensor (hence, the memory required to store it) and the amount of operations required to perform basic operations with such tensor grows exponentially in the dimensionality, and this problem is called the curse of dimensionality.
8.^All basic operations in the TT-format are implemented in the ttpy python package https://github.com/oseledets/ttpy and its MATLAB version https://github.com/oseledets/TT-Toolbox.
9.^An exact TT-representation exists for the given full tensor , and TT-ranks of such representation are bounded by ranks of the corresponding unfolding matrices (). Nevertheless, in practical applications it is more useful to construct TT-approximation with a prescribed accuracy , and then carry out all operations (summations, products, etc) in the TT-format, maintaining the same accuracy of the result.
10.^Note that in Algorithm 4 we use the solution of the convection part from the previous time step () as an initial guess for CAM on the next time step (k). As it was found empirically It may seem more logical to use the solution of the diffusion part from the same time step (k) as an initial guess, but we found empirically that it leads to higher TT-ranks of the result.
11.^Actual links to the corresponding Colab notebooks are available in our public repository https://github.com/AndreiChertkov/fpcross.
12.^Hereinafter, we present effective TT-rank of the computation result. For TT-tensor with TT-ranks () the effective TT-rank is a solution of quadratic equation The representation with a constant TT-rank (, , ) yields the same total number of parameters as in the original decomposition of the tensor .
13.^This choice of parameters corresponds to the problem of polymer modeling from the work (). In the corresponding model, the molecules of the polymer are represented by beads and interactions are indicated by connecting springs. Accordingly, for the case of only two particles we come to the dumbbell problem, which can be mathematically written in the form of the FPE.
14.^As values for comparison, we used the result of the most accurate calculation from work (), within which , and .
References
1
ChenR. T.DuvenaudD. (2019). Neural Networks with Cheap Differential Operators. New York, NY: arXiv preprint arXiv:1912.03579.
2
ChenT. Q.RubanovaY.BettencourtJ.DuvenaudD. K. (2018). Neural Ordinary Differential Equations. Adv. Neural Inf. Process. Syst., 6571–6583.
3
CichockiA.LeeN.OseledetsI.PhanA.-H.ZhaoQ.MandicD. P. (2016). Tensor Networks for Dimensionality Reduction and Large-Scale Optimization: Part 1 Low-Rank Tensor Decompositions. FNT Machine Learn.9, 249–429. 10.1561/2200000059
4
DolgovS.SavostyanovD. (2020). Parallel Cross Interpolation for High-Precision Calculation of High-Dimensional Integrals. Comp. Phys. Commun.246, 106869. 10.1016/j.cpc.2019.106869
5
DolgovS. V. (2019). A Tensor Decomposition Algorithm for Large Odes with Conservation Laws. Comput. Methods Appl. Math.19, 23–38. 10.1515/cmam-2018-0023
6
DolgovS. V.KhoromskijB. N.OseledetsI. V. (2012). Fast Solution of Parabolic Problems in the Tensor Train/Quantized Tensor Train Format with Initial Application to the Fokker--Planck Equation. SIAM J. Sci. Comput.34, A3016–A3038. 10.1137/120864210
7
FoxC.DolgovS.MorrisonM. E.MoltenoT. C. (2020). Grid Methods for Bayes-Optimal Continuous-Discrete Filtering and Utilizing a Functional Tensor Train Representation. Inverse Probl. Sci. Eng., 1–19. 10.1080/17415977.2020.1862109
8
GlowinskiR.OsherS. J.YinW. (2017). Splitting Methods in Communication, Imaging, Science, and Engineering. Springer.
9
GrathwohlW.ChenR. T.BettencourtJ.SutskeverI.DuvenaudD. (2018). Ffjord: Free-form Continuous Dynamics for Scalable Reversible Generative Models. arXiv preprint arXiv:1810.01367.
10
KidgerP.FosterJ.LiX.OberhauserH.LyonsT. (2021). Neural Sdes as Infinite-Dimensional gans. arXiv preprint arXiv:2102.03657.
11
KikuchiK.YoshidaM.MaekawaT.WatanabeH. (1991). Metropolis Monte Carlo Method as a Numerical Technique to Solve the Fokker-Planck Equation. Chem. Phys. Lett.185, 335–338. 10.1016/s0009-2614(91)85070-d
12
KüchlinS.JennyP. (2017). Parallel Fokker-Planck-DSMC Algorithm for Rarefied Gas Flow Simulation in Complex Domains at All Knudsen Numbers. J. Comput. Phys.328, 258–277. 10.1016/j.jcp.2016.10.018
13
KumarP.NarayananS. (2006). Solution of Fokker-Planck Equation by Finite Element and Finite Difference Methods for Nonlinear Systems. Sadhana31, 445–461. 10.1007/bf02716786
14
OseledetsI.TyrtyshnikovE. (2010). Tt-cross Approximation for Multidimensional Arrays. Linear Algebra its Appl.432, 70–88. 10.1016/j.laa.2009.07.024
15
OseledetsI. V. (2010). Approximation of $2^d\times2^d$ Matrices Using Tensor Decomposition. SIAM J. Matrix Anal. Appl.31, 2130–2145. 10.1137/090757861
16
OseledetsI. V. (2011). Tensor-train Decomposition. SIAM J. Sci. Comput.33, 2295–2317. 10.1137/090752286
17
OseledetsI. V.TyrtyshnikovE. E. (2009). Breaking the Curse of Dimensionality, or How to Use Svd in many Dimensions. SIAM J. Sci. Comput.31, 3744–3759. 10.1137/090748330
18
PichlerL.MasudA.BergmanL. A. (2013). “Numerical Solution of the Fokker-Planck Equation by Finite Difference and Finite Element Methods-A Comparative Study,” in Computational Methods in Stochastic Dynamics (Springer), 69–85. 10.1007/978-94-007-5134-7_5
19
SavostyanovD.OseledetsI. (2011). “Fast Adaptive Interpolation of Multi-Dimensional Arrays in Tensor Train Format,” in The 2011 International Workshop on Multidimensional (nD) Systems (IEEE), 1–8. 10.1109/nds.2011.6076873
20
SinghR.GhoshD.AdhikariR. (2018). Fast Bayesian Inference of the Multivariate Ornstein-Uhlenbeck Process. Phys. Rev. E98, 012136. 10.1103/PhysRevE.98.012136
21
SubramaniamG. M.VedulaP. (2017). A Transformed Path Integral Approach for Solution of the Fokker-Planck Equation. J. Comput. Phys.346, 49–70. 10.1016/j.jcp.2017.06.002
22
SunY.KumarM. (2015). A Numerical Solver for High Dimensional Transient Fokker-Planck Equation in Modeling Polymeric Fluids. J. Comput. Phys.289, 149–168. 10.1016/j.jcp.2015.02.026
23
SunY.KumarM. (2014). Numerical Solution of High Dimensional Stationary Fokker-Planck Equations via Tensor Decomposition and Chebyshev Spectral Differentiation. Comput. Math. Appl.67, 1960–1977. 10.1016/j.camwa.2014.04.017
24
TrefethenL. N. (2000). Spectral Methods in MATLAB, Vol. 10. Philadelphia, PA: Siam.
25
TyrtyshnikovE. (2000). Incomplete Cross Approximation in the Mosaic-Skeleton Method. Computing64, 367–380. 10.1007/s006070070031
26
VatiwutipongP.PhewcheanN. (2019). Alternative Way to Derive the Distribution of the Multivariate Ornstein-Uhlenbeck Process. Adv. Differ. Equ2019, 276. 10.1186/s13662-019-2214-1
27
VenkiteswaranG.JunkM. (2005). A QMC Approach for High Dimensional Fokker-Planck Equations Modelling Polymeric Liquids. Mathematics Comput. Simulation68, 43–56. 10.1016/j.matcom.2004.09.002
28
WehnerM. F.WolferW. G. (1983). Numerical Evaluation of Path-Integral Solutions to Fokker-Planck Equations. Phys. Rev. A.27, 2663–2670. 10.1103/physreva.27.2663
Summary
Keywords
fokker-planck equation, probability density function, tensor train format, cross approximation, chebyshev polynomial, ornstein-uhlenbeck process, dumbbell model
Citation
Chertkov A and Oseledets I (2021) Solution of the Fokker–Planck Equation by Cross Approximation Method in the Tensor Train Format. Front. Artif. Intell. 4:668215. doi: 10.3389/frai.2021.668215
Received
15 February 2021
Accepted
29 June 2021
Published
02 August 2021
Volume
4 - 2021
Edited by
Evangelos Papalexakis, University of California, Riverside, United States
Reviewed by
Devin Matthews, Southern Methodist University, United States
Prakash Vedula, University of Oklahoma, United States
Updates
Copyright
© 2021 Chertkov and Oseledets.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Andrei Chertkov, a.chertkov@skoltech.ru
This article was submitted to Machine Learning and Artificial Intelligence, a section of the journal Frontiers in Artificial Intelligence
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.