Front. Artif. Intell. Frontiers in Artificial Intelligence Front. Artif. Intell. 2624-8212 Frontiers Media S.A. 668215 10.3389/frai.2021.668215 Artificial Intelligence Original Research Solution of the Fokker–Planck Equation by Cross Approximation Method in the Tensor Train Format Chertkov and Oseledets Tensor Train for Fokker–Planck Equation Chertkov Andrei * Oseledets Ivan Skolkovo Institute of Science and Technology, Moscow, Russia

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

 *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

 02 08 2021 2021 4 668215 15 02 2021 29 06 2021 Copyright © 2021 Chertkov and Oseledets. 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.

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.

 fokker-planck equation probability density function tensor train format cross approximation chebyshev polynomial ornstein-uhlenbeck process dumbbell model 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 (Grathwohl et al., 2018) for neural ordinary differential equation (ODE) (Chen et al., 2018; Chen and Duvenaud, 2019), generative models (Kidger et al., 2021), etc.

Consider a stochastic dynamical system which is described by stochastic differential equation (SDE) of the form 1 d x = f ( x , t )   d t + S ( x , t )   d β ,   d β   d β ⊤ = Q ( t )   d t ,   x = x ( t ) ∈ ℝ d , where d β is a q-dimensional space-time white noise, f is a known d-dimensional vector-function and S ∈ ℝ d × q , Q ∈ ℝ q × q are known matrices. The FPE for the corresponding probability density function (PDF) ρ ( x , t ) of the spatial variable x has the form ∂ ρ ( x , t ) ∂ t = ∑ i = 1 d ∑ j = 1 d ∂ ∂ x i ∂ ∂ x j [ D i j ( x , t ) ρ ( x , t ) ] − ∑ i = 1 d ∂ ∂ x i [ f i ( x , t ) ρ ( x , t ) ] , where D ( x , t ) = 1 2 S ( x , t ) Q ( t ) S ⊤ ( x , t ) 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 (Dolgov et al., 2012). 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) (Oseledets, 2011). 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 S = 0 , i.e. in the deterministic case. For this case the evolution of the PDF along the trajectory is given by the formula ∂ ρ ( x , t ) ∂ t = − Tr ( ∂ f ( x , t ) ∂ x )   ρ ( x , t ) , where Tr (   ⋅   ) is a trace operation for the matrix. Hence, to compute the value of ρ ( x , t ) at the specific point x = x ^ , it is sufficient to find a preimage x ^ 0 such that if it is used as an initial condition for eq. 1, then we arrive to x ^ . 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 ρ ( x , t ) at any x ^ , we can use the cross approximation method (CAM) (Oseledets and Tyrtyshnikov, 2010; Savostyanov and Oseledets, 2011; Dolgov and Savostyanov, 2020) 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 S ≠ 0 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 t + h .

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 code 2 , 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 β ∈ ℝ q has the same dimension as the spatial variable x ∈ ℝ d ( q = d ), and the matrices in eq. 1 and eq. 2 have the form 3 Q ( t ) ≡ I d ,   S ( x , t ) ≡ 2 χ I d ,   D ( x , t ) ≡ χ I d , where χ ≥ 0 is a scalar diffusion coefficient. Then eqs 1, 2 can be rewritten in a more compact form d x = f ( x , t )   d t + 2 χ d β ,   d β   d β ⊤ = I d   d t , ∂ ρ ∂ t = χ Δ ρ − div [ f ( x , t ) ρ ] , where d-dimensional spatial variable x = x ( t ) ∈ Ω ⊂ ℝ d has the corresponding PDF ρ ( x , t ) with initial conditions x ( 0 ) = x 0 ∼ ρ ( x , 0 ) ,   ρ ( x , 0 ) = ρ 0 ( x ) . To construct the PDF at some moment τ ( τ > 0 ) for the known initial distribution ρ 0 ( x ) , we discretize eqs 5, 6 on the uniform time grid with M ( M ≥ 2 ) points t m = m h ,   h = τ M − 1 ,   m = 0,1 , … , M − 1 , and introduce the notation x m = x ( t m ) for value of the spatial variable at the moment t m and ρ m (   ⋅   ) = ρ (   ⋅   , t m ) for values of the PDF at the same moment.

2.1 Splitting Scheme

Let V ^ and W ^ be diffusion and convection operators from the eq. 6 V ^ v ≡ χ Δ v ,   W ^ w ≡ − div [ f ( x , t ) w ] , then on each time step m ( m = 0,1 , … , M − 2 ) we can integrate equation ∂ ρ ∂ t = ( V ^ + W ^ ) ρ ,   ρ ( ⋅ , t m ) = ρ m ( ⋅ ) , on the interval ( t m , t m + h ) , to find ρ m + 1 for the known value ρ m from the previous time step. Its solution can be represented in the form of the product of an initial solution with the matrix exponential ρ m + 1 = e h ( V ^ + W ^ ) ρ m , and if we apply the standard second order operator splitting technique (Glowinski et al., 2017), then ρ m + 1 ≈ e   h 2 V ^ e h W ^ e h 2 V ^ ρ m , which is equivalent to the sequential solution of the following equations ∂ v ( 1 ) ∂ t = χ Δ v ( 1 ) ,   v ( 1 ) (   ⋅   , t m ) = ρ m (   ⋅   ) , ∂ w ∂ t = − div [ f ( x , t ) w ] ,   w (   ⋅   , t m ) = v ( 1 ) (   ⋅   , t m + h 2 ) , ∂ v ( 2 ) ∂ t = χ Δ v ( 2 ) ,   v ( 2 ) (   ⋅   , t m ) = w (   ⋅   , t m + h ) , with the final approximation of the solution ρ m + 1 (   ⋅   ) = v ( 2 ) (   ⋅   , t m + h 2 ) .

 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 (Trefethen, 2000).

We introduce the d-dimensional spatial grid X ( g ) as a tensor product of the one-dimensional grids 4 x k ( g ) ∈ ℝ N k ,   x k ( g ) [ n k ] = cos π ⋅ ( n k − 1 ) N k − 1 ,   n k = 1,2 , … , N k , where N k ( N k ≥ 2 ) is a number of points along the kth spatial axis ( k = 1,2 , … , d ), and the total number of the grid points is N = N 1 ⋅ N 2 ⋅ … ⋅ N d . Note that this grid can be also represented in the flatten form as a following matrix X ( g ) ∈ ℝ d × N ,   X ( g ) [ k , n ] = x k ( g ) [ mind ( n ) [ k ] ] , where n = 1,2 , … , N , k = 1,2 , … , d and by mind ( n ) = [ n 1 , n 2 , … , n d ⊤ we denoted an operation of construction of the multi-index from the flatten long index according to the big-endian convention n = n d + ( n d − 1 − 1 ) N d + … + ( n 1 − 1 ) N 2 N 3 … N d . Suppose that we calculated PDF ρ m on some time step m ( m ≥ 0 ) at the nodes of the spatial grid X ( g ) [note that for the case m = 0 , the corresponding values come from the known initial condition ρ 0 ( x ) ]. These values can be collected as elements of a tensor 5 ℛ m ∈ ℝ N 1 × N 2 × … × N d such that ℛ m [ n 1 , n 2 , … , n d ] = ρ m ( x 1 ( g ) [ n 1 ] , x 2 ( g ) [ n 2 ] , … , x d ( g ) [ n d ] ) , where n k = 1,2 , … , N k ( k = 1,2 , … , d ).

Let us interpolate PDF ρ m via the system of orthogonal Chebyshev polynomials of the first kind T 0 ( x ) = 1 ,   T 1 ( x ) = x ,   T k + 1 ( x ) = 2 x T k ( x ) − T k − 1 ( x )   f o r   k = 1,2 , … , in the form of the naturally cropped sum ρ m ( x ) ≈ ρ m ˜ ( x ) = = ∑ n 1 = 1 N 1 ∑ n 2 = 1 N 2 ⋅ ⋅ ⋅ ∑ n d = 1 N d A m [ n 1 , n 2 , … , n d ]   T n 1 − 1 ( x 1 ) T n 2 − 1 ( x 2 ) … T n d − 1 ( x d ) , where x = ( x 1 , x 2 , … , x d ) is some spatial point and interpolation coefficients are elements of the tensor A m ∈ ℝ N 1 × N 2 × … × N d . For construction of this tensor we should set equality in the interpolation nodes eq. 16 ρ m ˜ ( x 1 ( g ) [ n 1 ] ,     x 2 ( g ) [ n 2 ] ,   … ,   x d ( g ) [ n d ] ) = ρ m ( x 1 ( g ) [ n 1 ] ,   x 2 ( g ) [ n 2 ] ,   … ,   x d ( g ) [ n d ] ) for all combinations of n k = 1,2 , … , N k ( k = 1,2 , … , d ).

Therefore the interpolation process can be represented as a transformation of the tensor ℛ m to the tensor A m 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) (Trefethen, 2000) 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 A m 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, (Trefethen, 2000)] D k ∈ ℝ N k × N k such that D k = D ˜ k D ˜ k , where for each spatial dimension k = 1,2 , … , d D ˜ k [ i , j ] = 2 ( N k − 1 ) 2 + 1 6 ,   i = j = 1 , − x k ( g ) [ j ] 2 ( 1 − ( x k ( g ) [ j ] ) 2 ) ,   i = j = 2,3 , … , N k − 1 , c i c j ( − 1 ) i + j x k ( g ) [ i ] − x k ( g ) [ j ] ,   i ≠ j ,   i , j = 2,3 , … , N k − 1 , − 2 ( N k − 1 ) 2 + 1 6 ,   i = j = N k , with c i = 2 if i = 1 or i = N k and c i = 1 otherwise, and one dimensional grid points x k ( g ) defined from eq. 16. Then discretized Laplace operator has the form 6 Δ = D 1   ⊗   I N 2   ⊗ … ⊗   I N d + I N 1   ⊗   D 2   ⊗ …   ⊗   I N d + … + I N 1 ⊗ I N 2 ⊗ … ⊗ D d . Let V m ∈ ℝ N 1 × N 2 × … × N d be the known initial condition for the diffusion equation on the time step m ( t m = m h ), then for the solution V m + 1 2 at the moment t m + h 2 we have vec ( V m + 1 2 ) = e h 2 χ Δ vec ( V m ) , where an operation vec ( ⋅ ) 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 to vec ( V m   +   1 2 ) = ( e h 2 χ D 1   ⊗   e h 2 χ D 2   ⊗   …   ⊗   e h 2 χ D d ) vec ( V m ) . If we can represent the initial condition V m 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 vec ( V m + 1 2 ) .

 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 form d x = f ( x , t )   d t ,   x = x ( t ) ∈ ℝ d ,   x ∼ ρ ( x , t ) . If we consider the differentiation along the trajectory of the particles, as was briefly described in the Introduction, then ( ∂ w ∂ t ) x = x ( t ) = ∑ k = 1 d ∂ w ∂ x k ∂ x k ∂ t + ∂ w ∂ t = ∑ k = 1 d ∂ w ∂ x k ∂ x k ∂ t − div [ f w ] = = ∑ k = 1 d ∂ w ∂ x k f k − ∑ k = 1 d ∂ f k ∂ x k w − ∑ k = 1 d f k ∂ w ∂ x k = − ∑ k = 1 d ∂ f k ∂ x k w , where we replaced the term ∂ w ∂ t by the right hand side of eq. 14 and ∂ x k ∂ t 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 system { ∂ x ∂ t = f ( x , t ) , ∂ w ∂ t = − Tr ( ∂ f ∂ x ( x , t ) ) w Let us integrate eq. 29 on a time step m ( m = 0,1 , … , M − 2 ). If we set any spatial grid point x * = X ( g ) [ : , n ] ( n = 1,2 , … , N ) as initial condition for the spatial variable, then we’ll obtain solution w ^ m + 1 for some point x ^ m + 1 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 x ^ m that will be transformed to the grid point x * by the step m + 1 . If we select this point x ^ m and the related value w ^ m = w ( x ^ m , t m ) as initial conditions for the system eq. 29, then its solution w m + 1 will be related to the point of interest x * .

Evolution of the spatial variable and the corresponding PDF for two consecutive time steps related to the fixed Chebyshev grid in the case of two dimensions.

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 w m is already known and defined as a tensor W m ∈ ℝ N 1 × N 2 × … × N d 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 A m . Then we can evaluate the approximated value at the point x ^ m as w m ˜ ( x ^ m ) according to eq. 21.

Hence our solution strategy for convection equation is the following. For the given spatial grid point x * = X ( g ) [ : , n ] we integrate equation ∂ x ∂ t = f ( x , t ) ,   x ( t m + 1 ) = x * , backward in time to find the corresponding point x ^ m = x ( t m ) . Then we find the value of w at this point, using interpolation w m ˜ , and then we solve the system eq. 29 on the time interval ( t m , t m + h ) with initial condition ( x ^ m , w m ˜ ( x ^ m ) ) to obtain the value w m + 1 at the point x * . The described process should be repeated for each grid point ( n = 1,2 , … , N ) and, ultimately, we’ll obtain a tensor W m + 1 ∈ ℝ N 1 × N 2 × … × N d 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 W m + 1 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 (Oseledets and Tyrtyshnikov, 2009; Oseledets, 2011), 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 dimensionality 7 .

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 W m 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 (Oseledets and Tyrtyshnikov, 2010) 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 ℛ ∈ ℝ N 1 × N 2 × … × N d is said to be in the TT-format (Oseledets, 2011), if its elements are represented by the formula ℛ [ n 1 , n 2 , … , n d ] = ∑ r 1 = 1 R 1 ∑ r 2 = 1 R 2 ⋅ ⋅ ⋅ ∑ r d − 1 = 1 R d − 1 G 1 [ 1 , n 1 , r 1 ] G 2 [ r 1 , n 2 , r 2 ] …                                                                                                               G d − 1 [ r d − 2 , n d − 1 , r d − 1 ] G d [ r d − 1 , n d , 1 ] where n k = 1,2 , … , N k ( k = 1,2 , … , d ), three-dimensional tensors G k ∈ ℝ R k − 1 × N k × R k are named TT-cores, and integers R 0 , R 1 , … , R d (with convention R 0 = R d = 1 ) are named TT-ranks. The latter formula can be also rewritten in a more compact form ℛ [ n 1 , n 2 , … , n d ] = G 1 ( n 1 ) G 2 ( n 2 ) … G d ( n d ) , where G k ( n k ) = G k [ : , n k , : ] is an R k − 1 × R k matrix for each fixed n k (since R 0 = R d = 1 , the result of matrix multiplications in eq. 32 is a scalar). And a vector form of the TT-decomposition looks like vec ( ℛ ) = ∑ r 1 = 1 R 1 ∑ r 2 = 1 R 2 ⋅ ⋅ ⋅ ∑ r d − 1 = 1 R d − 1 G 1 [ 1 , : , r 1 ] ⊗ G 2 [ r 1 , : , r 2 ] ⊗ … ⊗ G d [ r d − 1 , : , 1 ] , where the slices of the TT-cores G k are vectors of length N k ( k = 1,2 , … , d ).

The benefit of the TT-decomposition is the following. Storage of the TT-cores G 1 , G 2 , … , G d requires less or equal than d × max 1 ≤ k ≤ d ( N k R k 2 ) memory cells (instead of N = N 1 N 2 … N d ∼ N 0 d cells for the uncompressed tensor, where N 0 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 format 8 is given in works (Oseledets and Tyrtyshnikov, 2009; Oseledets, 2011). 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 ϵ T T in the Frobenius norm 9 ‖ ℛ − ℛ ^ ‖ F ≤ ϵ T T ⋅ ‖ ℛ ^ ‖ F , but 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 ϵ C A , using only part of the full tensor elements. This method is a multi-dimensional analogue of the simple cross approximation method for the matrices (Tyrtyshnikov, 2000) that allows one to approximate large matrices in O ( N 0 R 2 ) time by computing only O ( N 0 R ) elements, where N 0 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 f ( n 1 , n 2 , … , n d ) , that returns the ( n 1 , n 2 , … , n d ) th entry of ℛ ^ for a given set of indices. This method requires only

O ( d × max 1 ≤ k ≤ d ( N k R k 3 ) ) operations for the construction of the approximation with a prescribed accuracy ϵ C A , where R 0 , R 1 , … , R d ( R 0 = R d = 1 ) are TT-ranks of the tensor ℛ [see detailed discussion of the CAM in (Oseledets and Tyrtyshnikov, 2010)]. It should be noted that TT-ranks can depend on the value of selected accuracy ϵ C A , but for a wide class of practically interesting tasks the TT-ranks are bounded or depend polylogarithmically on ϵ C A [(Oseledets, 2010; Oseledets, 2011) for more details and examples].

In Algorithm 1 the description of the process of construction of the tensor in the TT-format on the Chebyshev grid by the CAM is presented (we’ll call it as a function cros s X (   ⋅   ) below). We prepare function func, which transforms given indices into the spatial grid points and return an array of the corresponding values of the target r (   ⋅   ) . Then this function is passed as an argument to the standard rank adaptive method tt_rectcross from the ttpy package. The CAM is described in more detail in the original papers (Oseledets and Tyrtyshnikov, 2010; Savostyanov and Oseledets, 2011), as well as in a recent work (Dolgov and Savostyanov, 2020), which formulates a computationally efficient parallel implementation of the algorithm.

 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 (Trefethen, 2000), 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 ℛ ∈ ℝ N 1 × N 2 × … × N d with the list of TT-cores [ G 1 , G 2 , … G d ] , which collects PDF values on the nodes of the Chebyshev grid at some time step (the related function is r ( x ) , 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 A ∈ ℝ N 1 × N 2 × … × N d of interpolation coefficients with the TT-cores [ G ˜ 1 , G ˜ 2 , … G ˜ d ] can be constructed according to the scheme, which is presented in Algorithm 2 [we’ll call it as a function interpolate (   ⋅   ) below].

In this Algorithm we use standard linear algebra operations swapaxes and reshape , which rearrange the axes and change the dimension of the given tensor respectively, function fft for construction of the one-dimensional FFT for the given vector, and function tt_round from the ttpy package, which round the given tensor to the prescribed accuracy ϵ. Note that the inner loop in Algorithm 2 for r * may be replaced by the vectorized computations of the corresponding two-dimensional FFT.

For the known tensor A we can perform a fast computation of the function value at any given spatial point x = [ x 1 , x 2 , … , x d ⊤ ] by a matrix product of the convolutions of the TT-cores of A with appropriate column vectors of Chebyshev polynomials r ( x ) ≈ ∑ r 1 = 1 R 1 ∑ r 2 = 1 R 2 … ∑ r d − 1 = 1 R d − 1 ( ∑ n 1 = 1 N 1 G ˜ 1 [ 1 , n 1 , r 1 ] T n 1 − 1 ( x 1 ) ) ( ∑ n 2 = 1 N 2 G ˜ 2 [ r 1 , n 2 , r 2 ] T n 2 − 1 ( x 2 ) ) … ( ∑ n d = 1 N d G ˜ d [ r d − 1 , n d , 1 ] T n d − 1 ( x d ) ) , We’ll call the corresponding function as inter _ eval ( A , X ) below. This function constructs a list of r (   ⋅   ) values for the given set of I points X ∈ ℝ d × I ( I ≥ 1 ), using interpolation coefficients A 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 cros s X , Algorithm 1) 10 and interpolation by the Chebyshev polynomials (function interpolate from Algorithm 2 that constructs interpolation coefficients and function inter_eval that evaluates interpolation result at given points according to the formula eq. 35).

We denote by einsum the standard linear algebra operation that evaluates the Einstein summation convention on the operands (see, for example, the numpy python package). Function vstack stack arrays in sequence vertically, function ode _ solve ( rhs , t 1 , t 2 , Y 0 ) (where t 1 and t 2 are initial and final times, rhs is the right hand side of equations, and matrix Y 0 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 interface 11 with the standard configuration (without GPU support).

Firstly we consider an equation with a linear convection term—Ornstein-Uhlenbeck process (OUP) (Vatiwutipong and Phewchean, 2019) 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 (Venkiteswaran and Junk, 2005) 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 ρ ( x , t ) | ∂ Ω ≈ 0 , and the initial condition is selected in the form of the Gaussian function ρ ( x , 0 ) = ρ 0 ( x ) = ( 2 π s ) − d 2 exp [ − 1 2 s | | x | | 2 ] ,   s ∈ ℝ ,   s > 0 , where parameter s is selected as s = 1 . To estimate the accuracy of the obtained PDF ( ρ ) we use the relative 2-norm of deviation from the exact value ( ρ e x a c t ) e = | | ρ − ρ e x a c t | | 2 | ρ e x a c t | | 2 . We compute the value ρ e x a c t 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. 6 in the d-dimensional case with f ( x , t ) = A ( μ − x ( t ) ) ,   χ = 1 2 ,   x ∈ Ω = [ x m i n , x m a x ] d ,   t ∈ [ 0 , τ ] , where A ∈ ℝ d × d is invertible real matrix, μ ∈ ℝ d is the long-term mean, x m i n ∈ ℝ and x m a x ∈ ℝ ( x m i n < x m a x ), τ ∈ ℝ ( τ > 0 ). This equation is a well known multivariate OUP with the following properties [see for example (Singh et al., 2018; Vatiwutipong and Phewchean, 2019)]:

• mean vector is

M ( t , x 0 ) = e − A t x 0 + ( I d − e − A t ) μ ;

• covariance matrix is

 Σ ( t ) = ∫ 0 t e A ( s − t ) S S ⊤ e A ⊤ ( s − t ) d   s , and, in our case as noted above S = 2 χ I d ;

• transitional PDF is

 ρ ( x , t , x 0 ) = exp [ − 1 2 ( x − M ( t , x 0 ) ) ⊤ Σ − 1 ( t ) ( x − M ( t , x 0 ) ) ] | 2 π Σ ( t ) | ;

• stationary solution is

 ρ s t ( x ) = exp [ − 1 2 x ⊤ W − 1 x ] ( 2 π ) d d e t ( W ) , where matrix W ∈ ℝ d × d can be found from the following equation A W + W A ⊤ = 2 χ I d ;

• 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 A > 0 in the one-dimensional case).

5.1.1 One-Dimensional Process

Let consider the one-dimensional ( d = 1 ) OUP with A = 1 ,   μ = 0 ,   x m i n = − 5 ,   x m a x = 5 ,   τ = 10. We can calculate the analytic solution in terms of only spatial variable and time via integration of the transitional PDF eq. 41 ρ ( x , t ) = ∫ − ∞ ∞ ρ ( x , t , x 0 ) ρ 0 ( x 0 )   d x 0 . Accurate computations lead to the following formula ρ ( x , t ) = 1 2 π ( Σ ( t ) + s e − 2 A t ) exp [ − x 2 2 ( Σ ( t ) + s e − 2 A t ) ] , where Σ ( t ) is defined by eq. 40 and for the one-dimensional case may be represented in the form Σ ( t ) = 1 − e − 2 A t 2 A . Using the formulas eq. 42 and eq. 43 we can represent a stationary solution for the one-dimensional case in the explicit form ρ s t a t ( x ) = A π e − A x 2 . We perform computation for N 1 = 50 spatial points and M = 1000 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 10 − 5 . At the same time, the numerical solution approaches the stationary one, and the corresponding error at large times also becomes approximately 10 − 5 . Note that the time to build the solution was about 5 s.

Relative error of the calculated solution vs known analytic and stationary solutions for the one-dimensional OUP.

 5.1.2 Three-Dimensional Process

Our next example is the three-dimensional ( d = 3 ) OUP with the following parameters A = [ 1.5 1 0 0 1 0 0.5 0.3 1 ] ,   μ = 0 ,   x m i n = − 5 ,   x m a x = 5 ,   τ = 5. When carrying out numerical calculation, we select 10 − 4 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-rank 12 remains limited, and the accuracy of the solution over time grows, reaching 10 − 3 by the time t = 5 . The time to build the solution was about 26 s.

Relative error of the calculated solution vs known stationary solution (A) and the effective TT-rank (B) for the three-dimensional OUP.

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 parameters A = [ 1.5 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0.5 0.3 0.2 0 1 ] ,   μ = 0 ,   x m i n = − 5 ,   x m a x = 5 ,   τ = 5. We select the same values as in the previous example for the CAM accuracy ( 10 − 4 ), 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 10 − 3 . The time to build the solution was about 100 s.

Relative error of the calculated solution vs known stationary solution (A) and the effective TT-rank (B) for the five-dimensional OUP.

 5.2 Numerical Solution of the Dumbbell Problem

Now consider a more complex non-linear example corresponding to the three-dimensional ( d = 3 ) dumbbell model of the form eq. 6 with 13 f ( x , t ) = A x − 1 2 ∇ ϕ ,   A = β [ 0 1 0 0 0 0 0 0 0 ] ,   ϕ = | | x | | 2 2 + α p 3 e − | | x | | 2 2 p 2 , where χ = 1 2 ,   x ∈ Ω = [ − 10,10 ] 3 ,   t ∈ [ 0,10 ] ,   α = 0.1 ,   β = 1 ,   p = 0.5. 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 ( k = 1,2,3 ) f = β [ x 2 0 0 ] − 1 2 x + α 2 p 5 e − | | x | | 2 2 p 2 x , ∂ f k ∂ x k = − 1 2 + α 2 p 5 e − | | x | | 2 2 p 2 − α 2 p 7 e − | | x | | 2 2 p 2 x k 2 . Next, we consider the Kramer expression τ ( t ) = ∫ ρ ( x , t ) [ x ⊗ ∇ ϕ ] d x , and as the values of interest (as in the works (Venkiteswaran and Junk, 2005; Dolgov et al., 2012)) we select ψ ( t ) = τ 11 ( t ) − τ 22 ( t ) β 2 = 1 β 2 ρ ( x , t ) ( x 1 ∂ ϕ ∂ x 1 − x 2 ∂ ϕ ∂ x 2 ) , η ( t ) = τ 12 ( t ) β = 1 β ρ ( x , t ) x 1 ∂ ϕ ∂ x 2 . During the calculations we used the following solver parameters:

• the accuracy of the CAM is 10 − 5 ;

• 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 ψ ( t ) and η ( t ) 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 ψ ( t ) and η ( t ) variables: ψ ( t = 10 ) = 2.0707 ,   η ( t = 10 ) = 1.0318 , with the corresponding results from (Dolgov et al., 2012) 14 , and we get the following values for relative errors ϵ ψ = 1.9 × 10 − 4 ,   ϵ η = 9.7 × 10 − 4 .

Computed values (A) and the effective TT-rank (B) for the three-dimensional dumbbell problem.

 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 (Wehner and Wolfer, 1983; Subramaniam and Vedula, 2017), the finite difference and the finite element method (Kumar and Narayanan, 2006; Pichler et al., 2013) 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 (Kikuchi et al., 1991; Küchlin and Jenny, 2017), 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 (Cichocki et al., 2016). However, for the FPE, this approach is not yet widely used. We note the works (Dolgov et al., 2012; Sun and Kumar, 2014; Sun and Kumar, 2015; Dolgov, 2019; Fox et al., 2020) 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 (Dolgov et al., 2012) 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.

 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.

Vectors and matrices are denoted hereinafter by lower case bold letters ( a , b , c , … ) and upper case letters ( A , B , C , … ) respectively. We denote the ( i 1 , i 2 ) th element of an N 1 × N 2 matrix A as A [ i 1 , i 2 ] and assume that 1 ≤ i 1 ≤ N 1 , 1 ≤ i 2 ≤ N 2 . For vectors we use the same notation: a [ i ] is the i-th element of the vector a ( i = 1,2 , … , N ). In addition, for a compact representation of an i-th ( i = 1,2 , … , d , where d ≥ 1 ) element of a vector function f =   [ f 1 , f 2 , … , f d ] ⊤ , we will use the notation f i , which means f i (   ⋅   ) = f i (   ⋅   ) [ i ] .

The code is publicly available from https://github.com/AndreiChertkov/fpcross.

We use notation I k for the k × k ( k = 1,2 , … ) identity matrix.

We suppose that for each spatial dimension the variable x varies within [ − 1 ,   1 ] . In other cases, an appropriate scaling can be easily applied.

By tensors we mean multidimensional arrays with a number of dimensions d ( d ≥ 1 ). A two-dimensional tensor ( d = 2 ) is a matrix, and when d = 1 it is a vector. For tensors with d > 2 we use upper case calligraphic letters ( A , ℬ , C , … ). The ( n 1 , n 2 , … , n d ) th entry of a d-dimensional tensor A ∈ ℝ N 1 × N 2 × … × N d is denoted by A [ n 1 , n 2 , … , n d ] , where n k = 1,2 , … , N k ( k = 1,2 , … , d ) and N k is a size of the k-th mode, and mode-k slice of such tensor is denoted by A [ n 1 , … , n k − 1 , : , n k + 1 , … , n d ] .

Note that for the case N 1 = N 2 = … = N d ≡ N 0 , we have only one matrix D 1 = D 2 = … = D d ≡ D 0 ∈ ℝ N 0 × N 0 which greatly simplifies the computation process.

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.

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.

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 (Oseledets, 2011). Nevertheless, in practical applications it is more useful to construct TT-approximation with a prescribed accuracy ϵ T T , and then carry out all operations (summations, products, etc) in the TT-format, maintaining the same accuracy ϵ T T of the result.

Note that in Algorithm 4 we use the solution of the convection part from the previous time step ( k − 1 ) as an initial guess W 0 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.

Actual links to the corresponding Colab notebooks are available in our public repository https://github.com/AndreiChertkov/fpcross.

Hereinafter, we present effective TT-rank of the computation result. For TT-tensor X ∈ ℝ N 1 × N 2 × … × N d with TT-ranks R 0 , R 1 , … , R d ( R 0 = R d = 1 ) the effective TT-rank R ^ is a solution of quadratic equation N 1 R ^ + ∑ α = 2 d − 1 N α R ^ 2 + N d R ^ = ∑ α = 1 d N α R α − 1 R α . The representation with a constant TT-rank R ^ ( R ^ 0 = 1 , R ^ 1 = R ^ 2 … = R ^ d − 1 = R ^ , R ^ d = 1 ) yields the same total number of parameters as in the original decomposition of the tensor X .

This choice of parameters corresponds to the problem of polymer modeling from the work (Venkiteswaran and Junk, 2005). 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.

As values for comparison, we used the result of the most accurate calculation from work (Dolgov et al., 2012), within which ψ ^ ( t = 10 ) = 2.071143 , and η ^ ( t = 10 ) = 1.0328125 .

 References Chen R. T. Duvenaud D. (2019). Neural Networks with Cheap Differential Operators. New York, NY: arXiv preprint arXiv:1912.03579. Chen T. Q. Rubanova Y. Bettencourt J. Duvenaud D. K. (2018). Neural Ordinary Differential Equations. Adv. Neural Inf. Process. Syst., 6571–6583. Cichocki A. Lee N. Oseledets I. Phan A.-H. Zhao Q. Mandic D. 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 Dolgov S. Savostyanov D. (2020). Parallel Cross Interpolation for High-Precision Calculation of High-Dimensional Integrals. Comp. Phys. Commun. 246, 106869. 10.1016/j.cpc.2019.106869 Dolgov S. V. (2019). A Tensor Decomposition Algorithm for Large Odes with Conservation Laws. Comput. Methods Appl. Math. 19, 23–38. 10.1515/cmam-2018-0023 Dolgov S. V. Khoromskij B. N. Oseledets I. 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 Fox C. Dolgov S. Morrison M. E. Molteno T. 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 Glowinski R. Osher S. J. Yin W. (2017). Splitting Methods in Communication, Imaging, Science, and Engineering. Springer. Grathwohl W. Chen R. T. Bettencourt J. Sutskever I. Duvenaud D. (2018). Ffjord: Free-form Continuous Dynamics for Scalable Reversible Generative Models. arXiv preprint arXiv:1810.01367. Kidger P. Foster J. Li X. Oberhauser H. Lyons T. (2021). Neural Sdes as Infinite-Dimensional gans. arXiv preprint arXiv:2102.03657. Kikuchi K. Yoshida M. Maekawa T. Watanabe H. (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 Küchlin S. Jenny P. (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 Kumar P. Narayanan S. (2006). Solution of Fokker-Planck Equation by Finite Element and Finite Difference Methods for Nonlinear Systems. Sadhana 31, 445–461. 10.1007/bf02716786 Oseledets I. Tyrtyshnikov E. (2010). Tt-cross Approximation for Multidimensional Arrays. Linear Algebra its Appl. 432, 70–88. 10.1016/j.laa.2009.07.024 Oseledets I. V. (2010). Approximation of $2^d\times2^d$ Matrices Using Tensor Decomposition. SIAM J. Matrix Anal. Appl. 31, 2130–2145. 10.1137/090757861 Oseledets I. V. (2011). Tensor-train Decomposition. SIAM J. Sci. Comput. 33, 2295–2317. 10.1137/090752286 Oseledets I. V. Tyrtyshnikov E. 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 Pichler L. Masud A. Bergman L. 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 Savostyanov D. Oseledets I. (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 Singh R. Ghosh D. Adhikari R. (2018). Fast Bayesian Inference of the Multivariate Ornstein-Uhlenbeck Process. Phys. Rev. E 98, 012136. 10.1103/PhysRevE.98.012136 Subramaniam G. M. Vedula P. (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 Sun Y. Kumar M. (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 Sun Y. Kumar M. (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 Trefethen L. N. (2000). Spectral Methods in MATLAB, Vol. 10. Philadelphia, PA: Siam. Tyrtyshnikov E. (2000). Incomplete Cross Approximation in the Mosaic-Skeleton Method. Computing 64, 367–380. 10.1007/s006070070031 Vatiwutipong P. Phewchean N. (2019). Alternative Way to Derive the Distribution of the Multivariate Ornstein-Uhlenbeck Process. Adv. Differ. Equ 2019, 276. 10.1186/s13662-019-2214-1 Venkiteswaran G. Junk M. (2005). A QMC Approach for High Dimensional Fokker-Planck Equations Modelling Polymeric Liquids. Mathematics Comput. Simulation 68, 43–56. 10.1016/j.matcom.2004.09.002 Wehner M. F. Wolfer W. G. (1983). Numerical Evaluation of Path-Integral Solutions to Fokker-Planck Equations. Phys. Rev. A. 27, 2663–2670. 10.1103/physreva.27.2663