Efficient spectral sparse grid approximations for solving multi-dimensional forward backward SDEs

This is the second part in a series of papers on multi-step schemes for solving coupled forward backward stochastic differential equations (FBSDEs). We extend the basic idea in our former paper [W. Zhao, Y. Fu and T. Zhou, SIAM J. Sci. Comput., 36 (2014), pp. A1731-A1751] to solve high-dimensional FBSDEs, by using the spectral sparse grid approximations. The main issue for solving high dimensional FBSDEs is to build an efficient spatial discretization, and deal with the related high dimensional conditional expectations and interpolations. In this work, we propose the sparse grid spatial discretization. We use the sparse grid Gaussian-Hermite quadrature rule to approximate the conditional expectations. And for the associated high dimensional interpolations, we adopt an spectral expansion of functions in polynomial spaces with respect to the spatial variables, and use the sparse grid approximations to recover the expansion coefficients. The FFT algorithm is used to speed up the recovery procedure, and the entire algorithm admits efficient and high accurate approximations in high-dimensions, provided that the solutions are sufficiently smooth. Several numerical examples are presented to demonstrate the efficiency of the proposed methods.

1. Introduction.Backward stochastic differential equation (BSDE) in the linear sense was first introduced by J.M.Bismut in 1973 [3].Then, in 1990, Pardoux and Peng showed the existence and uniqueness of the adapted solution for nonlinear BSDE for the first time [21].Since then, forward backward stochastic differential equations(FBSDEs) have been extensively studied, and have been shown disperse applications in different fields, such as stochastic optimal control, nonlinear filtering, nonlinear expectations, ect.The general high-dimensional FBSDEs defined on a complete probability space (Ω, F , F, P) take the following form (1.1) where t ∈ [0, T ] with T being the fixed time horizon; F = (F t ) 0≤t≤T is the natural filtration of the standard d-dimensional Brownian motion W = (W t ) 0≤t≤T ; X 0 ∈ F 0 and ξ ∈ F T are the initial and terminal conditions for the forward stochastic differential equation(SDE) and BSDE respectively; b : are the unknowns.It is worth to note that b(•, x, y, z), σ(•, x, y, z) and f (•, x, y, z) are all F t -adapted for fixed x, y and z, and that the two stochastic integrals with respect to W s are of the Itô type.A triple (X t , Y t , Z t ) is called an L 2 -adapted solution for FBSDEs (1.1) if it is F t -adapted, square integrable and satisfies the FBSDEs (1.1).FBSDEs (1.1) is called decoupled when b and σ are both independent of Y and Z.In this paper, we consider the numerical solution for FBSDEs (1.1) with ξ = ϕ(X T ), b(t, X t , Y t , Z t ), σ(t, X t , Y t , Z t ) and f (t, X t , Y t , Z t ) being deterministic functions.Due to the complex solution structure, solutions of FBSDEs in closed form can seldom be constructed.However, for decoupled FBSDEs, Peng [23] introduced the following nonlinear Feynman-Kac formula, which established a deep relationship between parabolic PDEs and FBSDEs (1.1): consider the following parabolic PDE (1.2) u t (t, x) + Lu(t, x) + f (t, x, u(t, x), u x σ(t, x)) = 0, t ∈ [0, T ), x ∈ R q u(T, x) = ϕ(x), where σσ ⊤ i,j ∂ xixj u(t, x).
If the above PDE (1.2) has a classical solution u(t, x) ∈ C 1,2 , then the triple X t , u(t, X t ), u x (t, X t )σ(t, X t ) solves the decoupled FBSDEs (1.1).
From the numerical point of view, one can use the above connections between PDEs and FBSDEs to design the so called probabilistic numerical methods for PDEs, by solving the equivalent FBSDEs.While there are a lot of works dealing with numerical schemes for BSDEs [4,2,5,7,15,33,36,34,24], however, there are only a few work on numerical methods for FBSDEs [8,18,16,17,31,35,37] and the second order FBSDEs (which are related to fully non-linear PDEs) [9,10,13,14].
The key issue in numerical methods is to balance the accuracy and computational complexity.Typically, the computational complexity increases dramatically as the dimension increases.Some of the above mentioned works are designed with high order accuracy that can however only be used to deal with low dimensional FBSDEs.While some of them are low order numerical methods that are suitable for solving high dimensional problems.In particular, we highlight the work [15], where the constructed numerical methods can deal with very high dimensional BSDEs, however, the convergence rate is only 1/2.We also mention the work [10], where a numerical example for a 12-dimensional coupled FBSDE is reported, and it is shown by numerical test that the numerical method converges with order 1.
In this work, we aim to design high order numerical schemes for multi-dimensional FBSDEs, by extending our previous work in [W.Zhao, Y. Fu and T. Zhou, SIAM J. Sci.Comput., 36 (2014), pp.A1731-A1751].The main difficulty for solving high dimensional FBSDEs is to efficiently evaluate the high dimensional conditional expectations and design the corresponding high dimensional interpolations.Here we shall use the sparse grid approximation technique to deal with these issues, and the FFT algorithm will be used to speed up the interpolation procedure.Several multi-dimensional examples with dimension up to 6 will be presented, and high order convergence rates up to 3 will be shown.
The rest of the paper is organized as follows.In Section 2, we present some preliminaries, and review the multi-step schemes introduced in our previous work [35].The sparse grid approximation will be discussed in Section 3, and this is followed by our fully discrete numerical schemes for high dimensional FBSDEs in Section 4. In Section 5, we shall present several numerical examples.Finally in Section 6 we give some concluding remarks.Now we introduce some notations to be used.
1.A ⊤ : the transpose of a vector/matrix A.
where α = (α 1 , α 2 , . . ., α q ), D ⊂ R q is a bounded domain and 2. Time discretization for FBSDEs.In this section, we shall briefly review the time-discrete schemes proposed in [35].To this end, let us consider the uniform time partition over [0, T ], with ∆t = T N and t i = i∆t, i = 0, 1, . . ., N .We set ∆W n,k = W t n+k − W tn , ∆W tn,t = W t − W tn , ∆t n,k = t n+k − t k and ∆t tn,t = t − t n for convenience.Then under certain regularity assumptions, for x ∈ R q the following two reference equations can be derived from the FBSDEs (1.1) (for details, one can refer to [35]): The multi-step schemes in [35] rely on efficiently approximating the derivatives in (2.1) and (2.2).The following classical approximations are used for functions One can easily check that du dt (t n ) admits the following representation where α k,i , i = 0, 1, . . ., k solve the following linear system (2.5) . Now by inserting the similar approach to equations (2.1) and (2.2) one gets To propose the multi-step schemes in [35], one needs the following local property of the generator of diffusion processes (see also in [35] for details) Theorem 2.1.Let t 0 < t be a fixed time, and x 0 ∈ R q be a fixed space point.If with bs = b(s, Xs ; t 0 , x 0 ), σs = σ(s, Xs ; t 0 , x 0 ) being smooth functions of (s, Xs ) with parameters (t 0 , x 0 ) that satisfy b(t 0 , Xt0 ; t 0 , x 0 ) = b(t 0 , x 0 ), σ(t 0 , Xt0 ; t 0 , x 0 ) = σ(t 0 , x 0 ).
By combining this local property and equations (2.6)-(2.7), the following multi-step numerical schemes are proposed in [35] Algorithm 1 Multi-step semi-discrete schemes for coupled FBSDEs Assume that Y N −i and Z N −i (i = 0, 1, . . ., k − 1) are known.For n = N − k, . . ., 0, solve X n,j (j = 1, 2, . . ., k), Y n = Y n (X n ) and 2. Set l = 0 and let ǫ 0 be a given tolerance.Solve Y n,l+1 (X n ) and Z n,l+1 (X n ) by the following steps, where Ȳ n+j are the values of Y n+j at the space point X n,j .
The above schemes are semi-discrete schemes proposed in [35], and it was shown that the k-th step scheme admits a k order convergence rate, provided that 1 ≤ k ≤ 6.To efficiently solve the FBSDEs, however, we also have to introduce a spatial discretizition, and should guarantee a high quality spatial approximation, e.g., the approximation of conditional expectations, to balance the entire numerical error.In fact, this is the main purpose of this work, and we will build efficient algorithms to deal with this issue in the following sections.Meanwhile, it is noted that an iterative procedure is used in the above schemes, this is due to the couple property of FBSDEs.For Decoupled FBSDEs, such a procedure may be omitted.

Sparse discretization and corresponding function approximations.
In the last section, we have introduced the high order semi-discrete schemes for coupled FBSDEs.However, to make the multi-step schemes more efficient, one should design efficient numerical methods for evaluating the conditional expectations and high dimensional interpolations.In [35], one uses the uniform tensor spatial meshes, and the tensorized Gaussian Hermite quadrature rule was used to evaluate the conditional expectations, moreover, the Lagrange interpolation method is used to compute the non-grid information.It was shown that such a combination is less efficient for high dimensional FBSDEs, as the required computational work increases exponentially as the dimension increases.To this end, we shall introduce the sparse grid approximation method, which is introduced originally for approximating high dimensional integral [1,28], to deal with this issue.We remark that the sparse grid approximation has been used in many different research topics, see e.g., [20,30,26,27,32] and references therein.
3.1.The multi-dimensional sparse grids.In this section, we follow closely the idea and notations in [28] to given a basic introduction for constructing high dimensional sparse grids.To begin, let us consider a sequence of grids The total number N i of points is usually chosen to be N i = 2 i + 1.Based on such an one dimensional sequence, one can build the q-dimensional sparse grids via where p ≥ q is an integer, and i = (i 1 , ..., i q ) is a multi-index.We call the onedimensional sequence nested if it satisfies In such cases, one can rearrange the sequence into the following hierarchical order: The nested structure admits many advantages in constructing high dimensional sparse grids.From the computational cost point of view, one needs only count each different point once, and then (3.1) can be written as In this paper, we shall adopt two papular types of one dimensional grids to construct the high dimensional sparse grid approximations.The first one is the sparse grid based on the one dimensional Chebyshev-Gauss-Lobatto(CGL) grids {C i } ∞ i=1 that is defined by It is easy to see that the CGL sequence is nested.Using such a nested property, in our framework, the sparse grid based on CGL points will be used to build high dimensional function approximations and related high dimensional interpolations.
The other type of sparse grid we shall use is the one based on the one dimensional Gauss-Hermite(GH) grid {G i } ∞ i=1 which is defined as where {x i j , j = 1, . . ., 2 i − 1} are roots of the Hermite polynomial of order 2 i − 1. Obviously, the GH sequence {G i } ∞ i=1 is not nested.Note that the above definitions are for the one dimensional case, and one can construct the associated high dimensional sparse grids (G p q and C p q ) based on these one-dimensional sets and the formula (3.1).One can refer to Fig. 1 to have a first glance at the two types of sparse grids in two dimensions.In this work, we shall use the sparse grids quadrature rule based on the GH sequence to approximate the associated high dimensional integrals (conditional expectations) in our multi-step schemes for solving FBSDEs.3.2.Function approximations on sparse grids.In this section, we discuss function approximations on sparse grids, and this will play an important role in our numerical schemes for solving FBSDEs.To this end, let us begin with the onedimensional case, and consider a bounded interval I ⊂ R in which we have a set of N i points χ i = {x k , k ∈ I i }.Let ω(x) > 0 be a weight function in I and {φ k } k∈I i be a set of basis functions (usually orthogonal polynomials) in L 2 ω (I).Then for f ∈ C(I), if the values of f at x ∈ χ i are known, we can construct an interpolation approximation of f via where the coefficients {b i k } are obtained by solving the following linear system Then, one can define the multi-variate interpolation approximation for R q -functions by the following sparse interpolation operator: or equivalently, (3.4) where p ≥ q is a positive integer, and with the basis functions φ k = q j=1 φ kj , k j ∈ I ij , and k = (k 1 , ..., k q ).To further reduce the computational complexity, we shall use hierarchical bases in this work.We recall the following definition: Definition 3.1.For nested grids, a set of basis functions By the above definition, it is easy to see that the expansion coefficients {b i k } k∈I i under the hierarchical bases will independent of the level index i.That is, for any f (x) ∈ C(I), we can write where Ĩi = I i \ I i−1 , and the coefficients {b k } k∈I i can be determined by Thus, in case a nested sparse grid and the corresponding hierarchical bases are used, the interpolation procedure (3.3) can be simplified as (3.5) where x = (x 1 , . . ., x q ), φk = q i=1 φki , and Then, the expansion coefficients {b k , k ∈ I p q } can be solved by Thus, I q d defines a unique interpolation operator onto the finite space In [27], the authors developed a fast transform to solve the equation (3.6).To illustrate the solving procedure, let's take I p q = I 4 2 as an example, i.e., To make it more clearly, first of all we consider how to compute the values of b k , k ∈ I 1 × I 1 by the following equation, The coefficients can be computed in two steps.
1. Perform the following transform on {f (x j ), j ∈ I 1 × I 1 } along the first dimension, and get {b ′ j , j ∈ I where matrix T is the inverse of ( φk (x j )) k,j .2. Perform the transform on {b ′ j , j ∈ I 1 × I 1 } along the second dimension, and get the coefficients {b j , j ∈ I 1 × I 1 }.
where matrix T is the inverse of ( φk (x j )) k,j .Based on this fact, we can compute the coefficients {b k , k ∈ I 4  2 } in the following way.Firstly, we perform one-dimensional transforms on {f (x j ), j ∈ I 3 × I 1 }, {f (x j ), j ∈ I 2 × Ĩ2 } and {f (x j ), j ∈ I 1 × Ĩ3 } along the first dimension, respectively.We get {b ′ j , j ∈ I 4 2 }.Then apply one-dimensional transforms on {b ′ j , j ∈ I 1 × I 3 }, {b ′ j , j ∈ Ĩ2 × I 2 } and {b ′ j , j ∈ Ĩ3 × I 1 } along the second dimension.We obtain all the values of {b j , j ∈ I 4  2 }.Now we extend this procedure to the q-dimensional case, and introduce the following algorithm.For more details about the fast transform, readers may refer to the work [27].
Algorithm 2 Fast transform on sparse grid.
Input: q, p, χ p q , {f (x j ), j ∈ I p q } and {φ k (x), k ∈ I p q } Output: {b k , k ∈ I p q |f (x j ) = k∈I p q b k φk (x j )} function FastTran(q,p,{f (x j ), j ∈ I p q }) b j ← f (x j ), for all j ∈ I p q for q ′ = 1 → q do for all i ′ ∈ {i ′ = (i 1 , . . ., i ..,jq T k,j q ′ end for end for return {b k , k ∈ I p q } end function 3.3.The approximation of conditional expectations.Recall our multistep schemes in Algorithm 1 for solving FBSDEs, to eventually solve the FBSDEs, one needs to do spatial discretizations.In this work, we aim at solving the solution pair (Y (x), Z(x)) on a given sparse grid x ∈ χ p q , on each time level.Note that other required information of (Y, Z) can be obtained by using the interpolation procedure based on the sparse grid information.Although the FBSDEs are essentially defined on the unbounded domain, however, we are usually interested in the values of (Y 0 , Z 0 ) in a bounded domain of x, which we denote by [a 0 , b 0 ].Then for each time level t n , we only need to approximate the solution of (Y tn , Z tn ) in a bounded domain [a n , b n ] so that the values of (Y tn , Z tn ) outside will not influence our computation of (Y 0 , Z 0 ) in [a 0 , b 0 ].More precisely, for the time level t = t n , we aim at solving (Y tn (x), Z tn (x)) for x ∈ χ pn q ⊂ [a n , b n ].To do this, one has to approximate the associated conditional expectations in Algorithm 1.As mentioned before, we shall adopt the sparse grid GH quadrature rule to approximate the conditional expectations.To this end, we propose the following choice for the bounded domain: for any x ∈ χ pn q , there holds (3.7) x + b(t n , x, y, z) for every j = 1, 2, . . ., k, where M is defined by where G p q is the sparse grid GH points that are used to approximate the conditional expectations, and The motivation for the above choice is that one needs to include all information that is used in the sparse grid GH quadrature rule.More precisely, consider our multi-step schemes in Algorithm 1, and suppose that (for j = 1, 2, . . ., k) the values of Y n+j (x) on x ∈ χ pn+j q are known.That is, we consider the k-step scheme.Then, we can compute the value of Y n+j (x) for every x ∈ [a n+j , b n+j ] by where the coefficients {β i } i∈I p n+j q can be obtained by Algorithm 2 based on the sparse grid information x ∈ χ pn+j q , and I pn+j q is the associated interpolation operator.Now, to solve (Y n , Z n ) for x ∈ χ pn q , we need to approximate the associated conditional expectations, e.g., E x tn Ȳ n+j .Note that we have In [35], we have proposed the tensor grid of GH points to approximation the above integral, namely, where ξ i and ω i are the GH quadrature points and the corresponding weights, respectively.We have defined the associated quadrature operator by Q p,T q .Then, it is clear to see that one needs the following function information Thus, we can see that the choice of (3.7) is reasonable because it makes every point z stays inside the interval [a n+j , b n+j ].
However, when the dimension of the Brownian motion grows higher, the computational requirements of the above tensor quadrature rule increase exponentially with respect to the dimension.Hence we shall resort to the sparse grid GH quadrature rule, in this work, to approximate the conditional expectations.We denote by Q p q the sparse grid Gauss-Hermite quadrature operator, namely, where the quadrature operator Q k denotes the one-dimensional Gauss-Hermite quadrature based on the grid G k .Then, the conditional expectations involved in Algorithm 1 can be approximated by where the two terms R y,n+j ), the error of the onedimensional Clenshaw-Curtis quadrature rule is given by, where Q stands for the quadrature using Clenshaw-Curtis points, N Q is the number of quadrature points and the constant C relies on the upper bound of the k-th derivative of f .From this one-dimensional estimates, the authors obtain the following result of the sparse grid quadrature for functions defined on a high-dimensional cube where N Q is the number of spare grid quadrature points and the constant C q,k depends on q and the upper bound of the k-th derivative of f .Analogously for functions f ∈ F k 1 (R), we have the following error estimates in [25], Therefore, by conducting a similar procedure in [1,19], we can obtain the following theorem.Theorem 3.2.For any functions f ∈ F k q (R q ), the error of sparse grid quadrature rule where N Q = N Q (q, p) is the number of the sparse grid quadrature points used by Q p q and the constant C q,k depends only on q and the upper bound of the k-th derivative of f .Note that the above quadrature rule uses non-grids information (i.e., not all grid points x belong to χ pn+j q . As mentioned, these information will be obtained via the sparse grid interpolation, i.e., the procedure (3.8).To this end, we propose the fully discrete schemes for approximating the conditional expectations: where pn+j , p n+j ≥ q are positive integers, and the interpolation errors are defined by We note that both the interpolation error and the quadrature error can be well controlled, provided that the functions admit certain regularities, and meanwhile, suitable sparse grids are used.For detailed error estimates of sparse grid quadratures & interpolations, one can refer to [1,20,26].Remark 3.3.In the above discussions, we have used the notations Q pn+j q and I pn+j q to stand for the quadrature and interpolation operators.In both notations, we use pn+j and p n+j to specify the different levels of sparse grids that are used for each time level t = t n+j .We remark that one can of course use uniform sparse grid in each time level, however, one would benefit if different level can be used according to the different time level, from the view of computational cost.In particular, for the time level with a larger domain [a k , b k ], one may need a high level sparse grids to obtain a good accuracy.
In what follows, We are aimed to show how to find the effective computational domains [a n , b n ] for each time level t = t n , and this is also a key issue for solving FBSDEs.In the beginning of this section, we have supposed that the following holds (3.13) for any x ∈ χ pn q , where χ pn q is the set of sparse grids for t = t n .For simple cases, the coefficients of the forward SDE are bounded functions, i.e., However, if b and σ are unbounded functions, the determination of these intervals becomes much more complex, as the values of b and σ keep changing during the iterative procedures.Nevertheless, one can always find a large computational domain to fix this issue (yet with huge computational cost).Another solution is to adopt an efficient approximation method for the whole space R q instead of using bounded domains approximations.And this is our ongoing project.
4. Fully discrete multi-step schemes for multi-dimensional FBSDEs.We summarize in this section the entire multi-step schemes for solving multi-dimensional FBSDEs.Let us consider the k-step scheme, and assume that we have obtained the initial values, i.e., Y N −k (x) and Z N −k (x), i = 0, 1, . . ., k − 1, are known for , we do the following steps: • We choose a suitable computational domain [a n , b n ] for t = t n , and construct the corresponding sparse grids C pn q ⊂ [a n , b n ].In our setting, for each time level, we construct the sparse grids by transforming the standard CGL sparse grids from [−1, 1] q to [a n , b n ].
• For each x ∈ C pn q we use Algorithm 1 to solve Y n (x), Z n (x) with X n = x.In this procedure, we will use the approximation methods in (3.11)-(3.12) to deal with the high dimensional conditional expectations.In particular, we shall use the following hierarchical basis functions { Tk (x)} k∈I pn q that is introduced in [26]: i.e, T k (x) is the transformed classical k-th order one dimensional Chebyshev polynomial.We note that by using this type of bases and the CGL sparse grid, the FFT algorithm can be used to speed up the recovery procedure in Algorithm 2, see also in [27].• One do the above procedure until the t = 0.A detailed description for our multi-step schemes is also shown in Algorithm 3.
Remark 4.1.Note that in the above algorithm, one needs to use some initial values, i.e., Y N −i (x), Z N −i (x) for i = 0, 1, . . ., k − 1.This can be obtained by running some standard low order numerical schemes with small time steps, or use the Runge-Kutta methods [5] or DC methods [31] to do the initialization.In our numerical examples, we shall directly set these values to be known to avoid the initialization error.

Numerical experiments.
In this section, we provide several constructive numerical examples to show the efficiency of our numerical algorithm for solving high-dimensional FBSDEs.We shall also present a numerical comparison between the spectral sparse grid(SSG) method and our previous approach in [35] (where standard Lagrange interpolations on tensor grids (LTG) are used).To specify the main differences between the two approaches, we list the main techniques used in both methods in Table 1.Algorithm 3 k-step scheme for solving high-dimensional FBSDEs {β n j , j ∈ I pn q } ← FastTran(p n , q, {Y n (x j ), j ∈ I pn q }).{γ n j , j ∈ I pn q } ← FastTran(p n , q, {Z n (x j ), j ∈ I pn q }).

end if end for end for
In what follows, we will denote by CR and RT the convergence rate and the running time, respectively.In all our numerical tests, the numerical results, which include numerical errors, convergence rates, and running times, are obtained on a computer with 16 Intel Xeon E5620 CPUs (2.40 GHz), and 3.0 GB free RAM, coding in FORTRAN 95.
Example 1: We first consider a two dimensional example, and we shall also report the numerical comparison between the SSG and LTG methods.In this two dimensional example, we set the components of b as and we set σ to be a diagonal matrix with diagonal components Furthermore, we choose a one dimensional generator function f as + sin(4(t + x 2 )) cos 2 (4(t + x 1 ))(sin(4(t + x 1 )) − 1) + sin(4(t + x 1 )) cos 2 (4(t + x 2 ))(sin(4(t + x 2 )) − 1), so that the exact solution is Due to the periodic property of this example, we can solve the FBSDEs in a fixed space domain [−π, π] 2 .For the SSG method, we shall use the standard Chebyshev-Gauss-Lobatto sparse grid C 7 2 for the spatial discretization, and the Gauss-Hermite sparse grid G 3 2 will be used for the high dimensional quadratures (conditional expectations).In the LTG method, as mentioned before, we shall use the tensor grid of uniform mesh, and the degree of Lagrangian interpolation polynomials n are decided by the following formula to balance the time discretization error and the interpolation error.For more detailed explanations of the LTG method, one can refer to [35].We now solve example 1 with the SSG method and the LTG method, and the errors, running times and convergence rates are shown in Table 2. From the above two tables, we immediately learn that both methods admit high order convergence rates, more precisely, the k-step schemes admits a k-order conver-gence rate (we only listed the numerical results for 1 ≤ k ≤ 3).However, the computational complexity exhibits big differences between the two proposed methods.For example, for 3-step methods with N = 128 (the 7th column), the running time is 197s for the SSG methods vases 4930s for the LTG methods, and this is obvious due to the efficient sparse grid discretization and the efficient sparse grid interpolations.
Example 2: our next example is q-dimensional decoupled FBSDEs.More precisely, we set the components of b as b i (t, x, y, z) = 1 q x i exp (−x 2 i ), i = 1, ..., q, and again, we set σ to be a diagonal matrix with diagonal components The generator function f is chosen to be It is easy to show that the exact solution takes the following form We solve the above FBSDEs with q = 3, 4, 5, 6, and the numerical results are listed in Table 3. Again, the proposed multi-step schemes admit high order convergence rates even for the 6-dimensional problem.To show the computational complexity of the schemes with respect to the dimension, we show in Fig. 2 the growth of running time against the dimension, and it is seems that the running time grows in certain polynomial level (non-exponential).
Example 3: we next consider the q-dimensional coupled FBSDEs.We set And the diagonal matrix σ is chosen with components σ i,i = t 2 sin 2 (y + x i ).The generator function is chosen as It can be checked that the exact solution is X 2 t,j (X t,j+1 + t) + X t,i   .
We solve this example for q = 2, 3, 4, 5, and the errors, running time and convergence rate are shown in Table 4. Again, high order convergence rates are obtained for this coupled example.6. Conclusions.In this work, we have extended our previous work [35] of multistep schemes to solve high-dimensional FBSDEs, by combining the sparse grid spatial discretizations and the sparse grid quadrature & interpolations, and the entire algorithm admits efficient and high accurate approximations in high-dimensions.It is shown that the proposed numerical techniques in this work are more efficient than our previous work in [35].However, we remark that how to build efficient high dimensional approximations is a long-term open question, and more efforts are still needed.Possible extensions along this direction include the adaptive sparse grid approaches, the anisotropic constructions of spatial approximations.Also, a rigorous error analysis will also be part of our future studies.
2. |•|: the Euclidean norm in the Euclidean space R, R q and R q×d .3. F t,x s : the σ-algebra generated by the diffusion process {X r , t ≤ r ≤ s, X t = x}.The symbols in bold denote the corresponding vector or multi-index.For a multi-index i = (i 1 , i 2 , . . ., i d ), |i| 1

Table 1 :
Main techniques used in the SSG and LTG methods for solving FBSDEs.TP: tensor product.SG: sparse grid.GH: Gaussian-Hermite.

Table 2 :
Numerical Errors and convergence rates for Example 1 by SSG (Top) and LTG (Bottom).

Table 3 :
Errors and convergence rates for Example 2.

Table 4 :
Errors and convergence rates for Example 3.