Solution of the Fokker–Planck Equation by Cross Approximation Method in the Tensor Train Format

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.


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 and Duvenaud, 2019), generative models (Kidger et al., 2021), etc.
Consider a stochastic dynamical system which is described by stochastic differential equation where dβ is a q-dimensional space-time white noise, f is a known d-dimensional vector-function and S ∈ R d×q , Q ∈ R q×q are known matrices. The FPE for the corresponding probability density function (PDF) ρ(x, t) of the spatial variable x has the form where D(x, t) 1 2 S(x, t)Q(t)S u (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 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 lowrank 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.

COMPUTATION OF THE PROBABILITY DENSITY FUNCTION
For ease of demonstration of the proposed approach, we suppose that the noise β ∈ R q has the same dimension as the spatial variable x ∈ R d (q d), and the matrices in eq. 1 and eq. 2 have the form 3 where χ ≥ 0 is a scalar diffusion coefficient. Then eqs 1, 2 can be rewritten in a more compact form where d-dimensional spatial variable x x(t) ∈ Ω ⊂ R d has the corresponding PDF ρ(x, t) with initial conditions 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 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.

Splitting Scheme
Let V and W be diffusion and convection operators from the eq. 6 then on each time step m (m 0, 1, . . . , M − 2) we can integrate equation 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 and if we apply the standard second order operator splitting technique (Glowinski et al., 2017), then which is equivalent to the sequential solution of the following equations with the final approximation of the solution ρ m+1 ( · ) v (2) · , t m + h 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 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 where n 1, 2, . . . , N, k 1, 2, . . . , d and by mind(n) [n 1 , n 2 , . . . , n u d we denoted an operation of construction of the multi-index from the flatten long index according to the bigendian convention 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 R m ∈ R N1×N2×...×N d such that 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 in the form of the naturally cropped sum where x (x 1 , x 2 , . . . , x d ) is some spatial point and interpolation coefficients are elements of the tensor A m ∈ R N 1 ×N 2 ×...×N d . For construction of this tensor we should set equality in the interpolation nodes eq. 16 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 R 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.

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) with c i 2 if i 1 or i N k and c i 1 otherwise, and one dimensional grid points x (g) k defined from eq. 16. Then discretized Laplace operator has the form 6 4 We suppose that for each spatial dimension the variable x varies within [−1, 1]. In other cases, an appropriate scaling can be easily applied. 5 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, B, C, . . .). The (n 1 , n 2 , . . . , n d ) th entry of a d-dimensional tensor A ∈ R N1 ×N2×...×Nd 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 ].
6 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 ∈ R N0×N0 which greatly simplifies the computation process.
Frontiers in Artificial Intelligence | www.frontiersin.org August 2021 | Volume 4 | Article 668215 Let V m ∈ R N1×N2×...×N d be the known initial condition for the diffusion equation on the time step m (t m mh), then for the solution V m+ 1 2 at the moment t m + h 2 we have 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 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 .

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 If we consider the differentiation along the trajectory of the particles, as was briefly described in the Introduction, then where we replaced the term zw zt by the right hand side of eq. 14 and zx k zt 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 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 * . 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 ∈ R 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 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 ∈ R N1×N2×...×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.
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 timeconsuming, so it may be useful to find some suitable lowparametric 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 TTformat 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.

Tensor Train Format
A tensor R ∈ R N 1 ×N 2 ×...×N d is said to be in the TT-format (Oseledets, 2011), if its elements are represented by the formula where n k 1, 2, . . . , N k (k 1, 2, . . . , d), three-dimensional tensors G k ∈ R 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 TTranks. The latter formula can be also rewritten in a more compact form R[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 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 × 0 cells for the uncompressed tensor, where N 0 is an average size of the tensor modes), and hence the TTdecomposition 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 R in the full format, the TT-decomposition (compression) can be performed by a stable TT-SVD algorithm. This algorithm constructs an approximation R in the TT-format to the given tensor R with a prescribed accuracy ϵ TT in the Frobenius norm 9 FIGURE 1 | 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.
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.

Cross Approximation Method
The CAM allows to construct a TT-approximation of the tensor with prescribed accuracy ϵ CA , 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 R to the tensor R, given as a function f (n 1 , n 2 , . . . , n d ), that returns the (n 1 , n 2 , . . . , n d ) th entry of R for a given set of indices. This method requires only ) operations for the construction of the approximation with a prescribed accuracy ϵ CA , where R 0 , R 1 , . . . , R d (R 0 R d 1) are TT-ranks of the tensor R [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 ϵ CA , but for a wide class of practically interesting tasks the TT-ranks are bounded or depend polylogarithmically on ϵ CA [ (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 cross 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.

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 R ∈ R N1×N2×...×N d with the list of TTcores [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 ∈ R 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 twodimensional 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 u d ] by a matrix product of the convolutions of the TT-cores of A with appropriate column vectors of Chebyshev polynomials 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 ∈ R d×I (I ≥ 1), using interpolation coefficients A and sequentially applying the formula eq. 35 for each spatial point.

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 cross 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.

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)| zΩ ≈ 0, and the initial condition is selected in the form of the Gaussian function 10 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. 11 Actual links to the corresponding Colab notebooks are available in our public repository https://github.com/AndreiChertkov/fpcross.
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 (ρ exact ) We compute the value ρ exact 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.

Numerical Solution of the Ornstein-Uhlenbeck Process
Consider FPE of the form eq. 6 in the d-dimensional case with where A ∈ R d×d is invertible real matrix, μ ∈ R d is the long-term mean, x min ∈ R and x max ∈ R (x min < x max ), τ ∈ R (τ > 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 • covariance matrix is and, in our case as noted above S 2χ I d ; • transitional PDF is • stationary solution is where matrix W ∈ R d×d 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 longterm mean μ as time t tends to infinity) if all eigenvalues of A are positive (if A > 0 in the one-dimensional case).
We can calculate the analytic solution in terms of only spatial variable and time via integration of the transitional PDF eq. 41 Accurate computations lead to the following formula where Σ(t) is defined by eq. 40 and for the one-dimensional case may be represented in the form Using the formulas eq. 42 and eq. 43 we can represent a stationary solution for the one-dimensional case in the explicit form 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.

Three-Dimensional Process
Our next example is the three-dimensional (d 3) OUP with the following parameters 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.
To evaluate the efficiency of the proposed algorithm in the TTformat, 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.

Five-Dimensional Process
This multidimensional case is considered in the same manner as the previous one. We select the following parameters 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 TTrank 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.

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  12 Hereinafter, we present effective TT-rank of the computation result. For TTtensor X ∈ R N1 ×N2×...×Nd 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 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 . 13 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.

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.