Semi-implicit Integration Factor Methods on Sparse Grids for High-Dimensional Systems

Numerical methods for partial differential equations in high-dimensional spaces are often limited by the curse of dimensionality. Though the sparse grid technique, based on a one-dimensional hierarchical basis through tensor products, is popular for handling challenges such as those associated with spatial discretization, the stability conditions on time step size due to temporal discretization, such as those associated with high-order derivatives in space and stiff reactions, remain. Here, we incorporate the sparse grids with the implicit integration factor method (IIF) that is advantageous in terms of stability conditions for systems containing stiff reactions and diffusions. We combine IIF, in which the reaction is treated implicitly and the diffusion is treated explicitly and exactly, with various sparse grid techniques based on the finite element and finite difference methods and a multi-level combination approach. The overall method is found to be efficient in terms of both storage and computational time for solving a wide range of PDEs in high dimensions. In particular, the IIF with the sparse grid combination technique is flexible and effective in solving systems that may include cross-derivatives and non-constant diffusion coefficients. Extensive numerical simulations in both linear and nonlinear systems in high dimensions, along with applications of diffusive logistic equations and Fokker-Planck equations, demonstrate the accuracy, efficiency, and robustness of the new methods, indicating potential broad applications of the sparse grid-based integration factor method.

where P is an elliptic differential operator with respect to the spatial variable x. This equation has been extensively studied because of its wide application in certain fields. For instance, the formation of the morphogen gradient during the development of the embryo is modeled using reaction-diffusion systems [1], where P denotes the Laplacian operator with respect to x. The stochastic behavior of a gene network can be described using the Fokker-Planck equation [2], also known as the backward Kolmogorov equation, where P is a second-order differential operator containing cross-derivatives. In finance, the Black-Scholes equation adopts a similar form when used to estimate the price of options under several risk factors [3]. In population genetics, the site-frequency spectrum can be modeled using such equations as well [4].
From the numerical perspective, solving Eq. (1) in high dimensions can be extremely challenging. Due to the "curse of dimensionality", achieving good accuracy of (for example, α = 2 by a second-order central difference formula) usually requires an number of points in uniform grids. The storage and operations on this large number of grid points can be prohibitively expensive when d is large. In addition, spatial discretization due to high-order spatial derivatives and stiff reactions leads to severe stability constraints on the time step for temporal integration. For the explicit methods, such as Runge-Kutta or Euler methods, a severe constraint is placed on the time step, whereas for implicit methods, such as the Crank-Nicolson method, large nonlinear systems are solved at each time step, leading to excessive computational costs.
The sparse grid technique has been shown to be an efficient approach for dealing with highorder spatial dimensions [5]. The discretization for this approach involves an O(N · (log N) d−1 ) number of points along with an accuracy of O((log N) d−1 /N 2 ) when the onedimensional piecewise linear hierarchical basis is applied. Such an approach may be extended to the general piecewise d-linear hierarchal basis by a tensor product construction in d-dimensional spaces [6,5]. Other hierarchical bases, e.g., high-order polynomials, interpolets, wavelets, have also been developed for a higher order of accuracy [7,8,9]. Moreover, sparse grids have been recently applied to stochastic simulations and parameter estimations [10,11,12].
A popular approach for addressing the problem associated with dimensionality is the sparse grid finite element method using the Galerkin technique [13,14] constructed on piecewise linear element using weak formulas. Another approach is the sparse grid finite difference method [15], which employs the regular second-order central difference approximation on the diffusion terms. The finite volume method [16] and the spectral method [17,18] can also be integrated with the sparse grid technique. Another interesting approach is the so-called sparse grid combination technique [19] that uses multi-level regular uniform grids such that the final solution is constructed using a linear combination of the intermediate solutions at the uniform grids, leading to a straightforward implementation, similar to the standard uniform grid approach.
For temporal integration, the integration factor (IF) and exponential time differencing (ETD) methods are effective ways to deal with the temporal stability constraints arising from highorder spatial derivatives on uniform meshes [20,21,22]. The IF and ETD methods usually treat linear operators of the highest-order derivatives exactly, and hence, they provide good temporal stability by allowing larger sizes of time step in temporal updates [23,24,20]. For addressing the stiffness in reactions, a class of semi-implicit integration factor (IIF) methods, which integrate the differential operators exactly like the IF schemes while treating the reaction terms implicitly, have been developed [25]. In IIF, the calculation of the diffusion and implicit treatment of the reaction is decoupled such that the size of the nonlinear system that needs to be solved at each time step is the same as that of the original continuous PDEs. This property results in good efficiency, in addition to excellent stability conditions (e.g., the second-order IIF is linearly unconditionally stable). Moreover, the IIF method can handle reaction-convection-diffusion equations through an operator splitting technique [26] and can be incorporated with the adaptive meshes and general curvilinear coordinates [27]. Because the exact treatment of the diffusion terms requires computing exponentials of matrices resulting from the discretization of the linear differential operators, a compact representation of IIF (cIIF) [28,27] and an array representation of IIF (AcIIF) [29] for systems with high spatial dimensions have been introduced. In the compact or array representations, the discretized functions in high spatial dimensions can be represented using multi-dimensional arrays rather than vectors or matrices so that the cost and storage associated with the calculation of the exponentials are significantly reduced.
In this paper, we integrate the sparse grids with the IIF methods to take the advantages provided by both methods to solve temporal PDEs (e.g., Eq. (1)) in high spatial dimensions.
In particular, we combine the two temporal schemes, the IIF and AcIIF methods, with the three different sparse grid discretization techniques: finite element, finite difference, and sparse grid combination technique. The combination technique is found to be especially effective in terms of incorporating various implicit integration factor methods, particularly when dealing with systems that include cross-derivatives and non-constant diffusion coefficients.
The paper is organized as follows. In Section 2, we construct the IIF method on sparse grids based on the finite element and Galerkin technique. In Section 3, we construct the AcIIF method on sparse grids using the finite difference approximation. In Section 4, we apply the AcIIF method to the sparse grid combination technique. In Section 5, we describe numerical tests to demonstrate the accuracy, efficiency, and applications of these methods.

Semi-implicit integration factor method with finite element method on sparse grids
where x = (x 1 , x 2 , …, x d ) and Δ x is the Laplacian operator with respect to x. For simplicity of presentation, we assume the spatial domain to be (0, 1) d and zero Dirichlet boundary conditions at the boundaries for u. We only focus on the piecewise d-linear hierarchical basis on the sparse grids.

Weak formula under sparse grid finite element method
Let The unknown u can be approximated using sparse grid, by (3) where k = (k 1 , …, k d ) and j = (j 1 , …, j d ). Define the inner product of two functions ψ 1 (x) and ψ 2 (x): (4) Substitute u in Eq. (2) by Eq. (3), and apply the inner product on both sides using ϕ k,j̄( x), for all possible k̄ and j. This leads to the following weak formula on the sparse grids: (5) where f k,j (t) is the hierarchical value of function f(u(x, t)): (6) Note that both < ϕ k,j , ϕ k,j̄ > and < Δ x ϕ k,j , ϕ k,j̄ > can be exactly calculated in advance because the piecewise d-linear functions are used. We define the following two matrices: (7) and let the vectors V(t) = (v k,j (t)) and F(t) = (f k,j (t)) denote the hierarchical values of the unknown function and the reaction, respectively. Then, Eq. (5) becomes the following timedependent problem: Instead of solving u, we only need to update the coefficient vector V at each time step according to Eq. (8). By directly applying the IIF method to Eq. (8), we obtain the following second-order (in time) IIF method based on the sparse grid finite element scheme (IIF2-SG-FEM): (9) Note that because the nodal value is usually not equal to the hierarchical value, the above nonlinear system arising from the implicit treatment of the reaction term on the right-hand side cannot be solved locally at each grid point. However, sparse grid discretization allows us to re-order the grid points such that only a localized nonlinear system at each point needs to be solved. The procedure is introduced below.

Solving the nonlinear system Eq. (9)
To solve Eq. (9), we first calculate the hierarchical value of k = (0, …, 0) and j = (1, …, 1). Note that at this point, the nodal value equals the hierarchical value, i.e., (10) so that we can directly solve the following equation: (11) where L k,j is the computed value of at the corresponding grid point. Then, we solve the hierarchical values for k such that k 1 + k 2 + … + k d = 1. Notice that there are two ways to compute the nodal value of f at these points: (i) from the hierarchical value of f (left side of Eq. (12)) and (ii) from the nodal value of u (right side of Eq. (12)). (12) Because v (0,…,0),(1,…,1) (t n+1 ) and f (0,…,0),(1,…,1) (t n+1 ) are known, substituting (12) into Eq. (9) for F(t n+1 ), we obtain (13) In general, with all v k,j (and the corresponding f k,j ) for which k 1 + … + k d < K obtained, the equality Eq. (14) follows and can be used to compute the nodal value of f at point x k,j̄ for k̄ satisfying , (14) where (15) Putting Eq. (14) into Eq. (9), we derive the local nonlinear system for v k,j̄ with k̄ satisfying . Recursively, we solve all unknowns in Eq. (9).

Error estimation, computational cost and stability
The error of the IIF2-SG-FEM method mainly comes from two sources: (i) spatial discretization by the sparse grids and (ii) the second-order IIF method for the temporal integration. Therefore, the estimated overall error is (16) The major cost of the IIF2-SG-FEM method associated with the temporal integration arises from calculating the exponential of M −1 DΔt. Although this calculation only needs to be performed once before the temporal update, it can be very expensive for high dimensions (d) and fine spatial resolution (N x ) because of the large size of the matrix of order O(N x (log N x )d −1 ). In addition, at each time step, one must solve the local nonlinear system obtained by substituting Eq. (14) into Eq. (9), leading to higher costs. Later, we discuss other methods that can reduce the size of the exponential matrix to reduce the computational costs.
Because M is the mass matrix and D is associated with the Laplacian operator, M −1 DΔt only contains non-negative eigenvalues, and the stability of IIF2-SG-FEM can be studied using a scalar case [25]. When the reaction term is linear, showing that IIF2-SG-FEM is A-stable is straightforward.

Array-representation semi-implicit integration factor method (AcIIF) with the sparse grids finite difference method
In this section, to reduce the size of the required exponential matrix and associated computational time, we construct the AcIIF method with the sparse grid finite difference scheme to solve Eq. (2), which has the same spatial domain and boundary conditions as that described in Section 2.

An introduction to the finite difference scheme on sparse grids
Following the methods previously studied in [15], let V(t) = (v k,j (t)) be the hierarchical value of the unknown u and U(t) = (u(x k,j , t)) be the nodal value. The nodal-hierarchical transformations, denoted by H̄ and M, are obtained using a dimensional splitting scheme [15]. Let H̄i be the (nodal to) hierarchical basis transformation along direction x i . It is equivalent to the one-dimensional (nodal to) hierarchical basis transformation acting on U k i j i(t) for all possible k i and j i , where (17) and U k i j i(t) is a sub-vector of U(t) obtained by fixing k s , j s , s ≠ i and only varying k i and j i . The (hierarchical to) nodal basis transformation along direction x i , M̄i, is defined similarly, except it acts on the hierarchical value V(t) by a one-dimensional (hierarchical to) nodal transformation. The nodal-hierarchical transformation is the composition of H̄i or M̄i: (18) Define H̄I /{i} := H̄1…H̄i −1 H̄i +1 …H̄d as the (nodal to) hierarchical basis transformation along all directions except for the direction x i . We also define the finite difference operator D̄i acting on any vector U(t) for all possible k i and j i to be the regular one-dimensional central Then, the finite difference scheme on sparse grids for is , with an estimated error [15] as follows: (19) With this approximation, Eq. (2) can be written into the following ODE system: (20) By applying the second-order IIF method, we obtain the following second-order (in time) IIF with sparse grid finite difference method (IIF2-SG-FD): (21) where . In the IIF2-SG-FD scheme, we need to compute and store the exponential of an When N x and d are large, the computation is intractable.

AcIIF method with the sparse grid finite difference method
Here, we apply the array representation to U(t) and V(t) to decompose them into smaller vectors [29]. Starting from the hierarchical value V(t), for a fixed direction x i and k i , j i , varying k i and j i gives a vector of T(k i , j i ) elements. Denote it by V k i ,j i(t), where (22) Let k i and j i go through all possible values; the resulting smaller vectors form the original vector V(t). We denote {(k i , j i ) : Σ r≠i k r ≤ K, j r = 1, 2, …, 2 k r , r ≠ i} by and use the notation ⊗ for such a representation: (23) Note that this representation is slightly different from the original array representation, where all of the subarrays have the same length [29]. In this case, the different subvectors , and the sizes are smaller than 2N x − 1. With the array representation, we can use smaller matrices to represent the (hierarchical to) nodal transformation along direction x i : (24) and the central difference scheme D̄i: (25) Here, M k i ,j i denotes the one-dimensional (hierarchical to) nodal transformation acting on V k i ,j i(t), and D k i ,j i is the tridiagonal matrix: (26) Both matrices have a size of

Define
, and apply the fact that ; under the array representation, the exponential of acting on U(t) becomes equivalent to (27) which is then multiplied by H̄− 1 to obtain the nodal value. Compared with the IIF method in which the exponential matrix is Through the array representation, i.e., utilizing multiple matrices of much smaller sizes instead of a single large exponential matrix, we obtain the array-representation semi-implicit integration factor method (AcIIF2-SG-FD), which can be summarized in the following steps: i. Transform the nodal value of to the hierarchical value W 1 (t), (28) ii. Adopt array representation to compute the exponential matrices and vector multiplication dimension by dimension: (29) iii. Update the nodal value by solving the nonlinear system:

Error estimation, computational cost, stability and high-order (in time) method
Naturally, the errors for the IIF2-SG-FD and AcIIF2-SG-FD methods take forms similar to that described in Eq. (16) because the sources of the discretization errors are similar.
Computing large matrix exponentials are required in IIF2-SG-FD and IIF2-SG-FEM. In contrast, in the AcIIF2-SG-FD method, the computational cost of the matrix exponential is significantly reduced, by a factor of at least (logN x ) d−1 /2. Even though some additional hierarchical-nodal transformations are introduced at each time step in AcIIF2-SG-FD, the cost is negligible relative to the cost of calculating the matrix exponentials because the fast computation of these transformations is straightforward [18]. Similar to IIF2-SG-FEM, the stability of both the IIF2-SG-FD and AcIIF2-SG-FD methods can be analyzed through the scalar case, and the methods are linearly A-stable [25].
To improve the order of accuracy for temporal integration, one can apply the third-order IIF (IIF3) to Eq. (20) in a straightforward way, leading to the IIF3-SG-FD method: (31) with error: (32)

AcIIF with sparse grids combination technique
In the previous sections, we incorporated the IIF and AcIIF methods with sparse grids in the finite element or finite difference methods. In both approaches, the formulations are most effective in dealing with the Laplacian differential operators in Eq. (1). Here, we incorporate a sparse grid combination technique with the AcIIF approach for more general cases for solving Eq. (1), where the spatial partial differential operator P may contain crossderivatives or non-constant diffusion coefficients.
In the sparse grid combination technique [19], instead of using non-uniform grids in a typical sparse grid discretization, the combination technique solves the PDEs on different levels of uniform grids and then combines the solutions for the final solution [19,30]. The sparse grid combination technique typically consists of two major steps: Step 1: Choose the spatial resolution N x = 2 K , and solve Eq. (1) on different uniform grids defined as below: (33) where 2 k i is the grid number along direction x i . The result on each grid is denoted by w k 1 ,…,k d (x, t).
Step 2: Combine all the solutions at different uniform grid levels as follows: (34) If each individual solution w k 1 ,…,k d (x, t) to the exact solution u(x, t) satisfies (35) then the overall error of w K (x, t) is [19] Wang et al.

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript (36) In the sparse grid combination technique, each subproblem contains 2 m grid points for K ≤ m ≤ K + d − 1, and the overall approximated solution has an error similar to the error of a typical sparse grid method. Each subproblem is independent of the other; thus, the implementation is much simpler than that for the finite difference or element approaches, and distributed (or parallel) calculations of each subproblem become straightforward.
If the explicit methods (e.g., RK2) are used to solve each subproblem, due to the stability constraints, the time step must be . Using IIF or AcIIF, the stability constraint for each subproblem is relaxed; in particular, the second-order IIF or AcIIF leads to no constraint on the temporal stability. If we set the uniform time step for each subproblem, e.g., Δt ~ 1/N x in AcIIF, the overall error in both space and time for each subproblem takes the form of Eq. (35).
In particular, on a uniform grid, after spatial discretization, the solution of Eq. (1), U, can be expressed as a d-dimensional array. The array representation for , the partial differential operator along the x i and x j directions, is denoted by . Thus, we obtain (37) where N r is the number of grid points in the x r direction. The exponential of is computed by the array representation: (38) Applying the implicit integration factor method yields the following second-order AcIIF method for Eq. (1) (39) where V(t n ) = U(t n ) + Δt/2F(U(t n )). The derivation of Eq. (39) is based on direct application of compact implicit integration factor method if 's are commutative. Otherwise, a splitting similar to Strang splitting method [31] needs to be used as shown in our previous work [29] to decompose the exponential matrices of Eq. (38) to obtain Eq. (39) for a second order method. The method is linearly absolute stable regardless of the commutative property [29]. Overall, the total computational cost for AcIIF along with the sparse grid combination technique is significantly reduced upon utilizing both advantages.

Numerical simulations
In this section, we investigate the three different sparse grid methods, all of which are based on the IIF or AcIIF temporal integration approach by applying them to five different PDE systems, where some have exact solutions. In particular, we examine the convergence and computational efficiency of the new methods and compare them with other existing methods for several cases.

A nonlinear reaction-diffusion system
We first implement the IIF2-SG-FEM method on the following nonlinear reaction-diffusion equation: (40) with zero Dirichlet boundary conditions and the initial condition: (41) The exact solution for this system is (42) Direct implementation of the weak formula Eq. (5) using the implicit methods (e.g., Crank-Nicolson) requires solving large nonlinear systems at each time step because of the presence of the nonlinear reaction u(x, y, t) 2 . Instead, we only make a comparison between the IIF2-SG-FEM method and RK2 for this case. In the RK2 method, if one chooses the time step , it is convergent, whereas a larger time step Δt = 1/N x leads to instability. In contrast, numerical results show that the IIF2-SG-FEM, in which the time step is chosen to be Δt = 1/N x , is convergent with the expected order of accuracy. The error of the IIF2-SG-FEM method is plotted in Figure 1. Because the discretization error is analytically, the errors measured under the log-log scale on different spatial resolutions should follow a straight line with a slope of one. This is shown in Figure 1, where the line of errors is almost parallel to the red line when the spatial resolution increases. In general, IIF2 can use a much larger time step in O(1/N x ) instead of for the RK2; thus, the efficiency of IIF2 is better than that of RK2.

A linear d-dimensional reaction-diffusion equation
We implement the AcIIF2-SG-FD method on the equation (43) with the initial condition (44) and zero Dirichlet boundary condition for d = 2, 3, 4, 6 respectively. The exact solution for the equation is: The L ∞ error versus for different spatial mesh sizes is plotted in Figure 2 to show the order of accuracy. The time step is set as Δt = 1/N x , and all simulations are terminated at t = 0.5 except that for d = 6. Compared with the red reference line with a slope of 1, in different dimensions, the numerical errors always fall on a straight line that is almost parallel to the red line as the spatial resolution increases, showing the expected order of accuracy.
One major advantage of the AcIIF method is the storage. For higher dimensions, such as d = 6, 8, 10, both IIF and RK2 methods rehire much more storage space. For example, when d = 8 and K = 6, IIF2-SG-FD and RK2 methods run out of memory for a computer of 16 GB memory whereas the AcIIF2-SG-FD method only requires matrix exponentials with a size no larger than 128 × 128. Nevertheless, when d = 8, 10, the CPU time taking for the standard sparse grid (mostly due to the transformations on different grid levels) is still very costly, in particular, for large K.
Another advantage of AcIIF is the efficiency. To demonstrate it, we compare it with IIF2 and RK2 based on the sparse grid finite difference method. Constrained by the stability condition, the time step for RK2 must be . For the other two methods, the time step is 1/N x . The CPU times consumed by different methods are listed in Table 1 for d = 2, 3, 4. As shown in the table, although RK2 requires less computational time on a coarse mesh, the IIF2 method is more efficient when the mesh is refined in a low-dimensional space. Overall, AcIIF2 performs most efficiently on the fine mesh across different dimensions. A detailed analysis shows that for the IIF2-SG-FD method, the most timeconsuming part is calculating the matrix exponential, consuming around 90 percent of the overall CPU time. In the AcIIF2-SG-FD method, the corresponding time spent on this calculation is significantly decreased. We also compare the CPU time by these methods when the error is comparable on different dimensions. For example, the error at N = 32 (or N = 64) for d = 3 is comparable to the error level at N = 64 (or N = 128) for the case of d = 4 in Tables 1 whereas the ratio between the computational time of AcIIF2 and RK2 changes significantly less than the ratio between d = 3 and d = 4 at the resolution N = 64 (or N = 128).

A diffusive logistic equation
The diffusive logistic equation in the following form: (46) describes the population density evolution, in which the unknown u denotes the population density of a species at time t and at location x [32]. The smooth function g defines the birth rate of the species, which may be positive or negative. Various numerical methods have been developed to solve this type of equation [33,34]. We use the AcIIF2-SG-FD method to solve the diffusive logistic equation in different dimensions. To study the accuracy of the method, we set up the equation for which an exact solution can be constructed.
and the zero Dirichlet boundary condition is imposed. The exact solution for this equation is (48) The implementation is the same as that for Section 5.2, except that the Newton method with a tolerance of 10 −6 is adopted to solve the small nonlinear system at each grid point, similar to the standard IIF method [25]. The log-log error plot for d = 2, 3, 4, 6 is displayed in Figure 3. From these plots, we can see that in different dimensions, the errors always follow a line with a slope close to 1 as the spatial resolution increases, consistent with the expected order of accuracy.
We next compare the CPU time used for AcIIF2-SG-FD and the RK2 method coupled with sparse grid finite difference for d = 2, 3, 4. The time step for RK2 is to guarantee that the stability condition will be satisfied. The CPU times for both methods are provided in Table 2. It is clear that regardless of dimensions of the system, AcIIF2 reduces the computational time significantly, particularly with fine meshes.

A three-dimensional reaction-diffusion system with cross-derivatives
To deal with cross-derivatives in the differential operators, we implement the AcIIF2 method with the sparse grid combination technique on this three-dimensional PDE: where x, y, z ∈ (0, 2π), with a periodic boundary condition and the initial condition: (50) The exact solution for the system is (51) The L ∞ error is plotted in Figure 4, which exhibits the correct and expected order of accuracy, i.e., parallel to the line of reference when the spatial resolution becomes large enough. Compared with the IIF method, AcIIF is clearly more effective in dealing with the matrix exponentials, as seen in Table 3, in which the CPU times for calculating the matrix exponentials, recognized to be the dominant cost for both methods, are listed.

A four-dimensional Fokker-Planck equation
In biology, the Fokker-Planck equation (FPE) is often used to describe and analyze the temporal evolution of the probability density functions for species in biochemical networks [2]. The generalized FPE takes the form: (52) In the FPE, R denotes the total number of chemical reactions taking place in the biochemical network, and d denotes the number of different species participating in the reactions. x = (x 1 , …, x d ) denotes a particular biochemical state, where each component x i denotes the copy number of the i-th reactant. The unknown function, p(x, t), is the probability density function for each state x at time t. n ri represents the change in the copy number of reactant i when the r-th reaction occurs once. In addition, one can define (53) where w r (x, t) is the reaction propensity function for the r-th reaction at state x.
FPE is a d-dimensional differential equation with non-constant diffusion coefficients and second-order cross-derivatives. Previously, the AcIIF scheme coupled with the finite difference method has been applied to the FPE on uniform grids, for which the AcIIF method has good efficiency [29]. Here, we implement the AcIIF method with a sparse grid combination technique to study the following biochemical reaction system [35,36,29]. In the metabolite-enzyme system, there are two metabolites A, B and two enzymes E A , E B : (54) where we set k A = k B = 0.3s −1 , k 2 = 0.001s −1 , K I = 60, μ = 0.002s −1 , k E A = k E B = 0.02s −1 and K R = 30. We start with a Gaussian distribution centered at (30,40,15,12) with a standard deviation of as the initial condition for p. There are nine reactions in the system, and for the AcIIF method, we group these reactions if their associated cross-derivative term is present. In this case, the first to fifth reactions are grouped in a one array-representation operator, the sixth and seventh are grouped into one, and the eighth and ninth are grouped into another. In Figure 5, we plot the density distribution at t = 45 solved using the sparse grid method. We also plot the trace of the ODE system, which is consistent with the FPE solution, suggesting the correctness of our approach.
To use the sparse grid combination technique, we first set up a spatial resolution N x = 2 K and then solve the FPE for each of its subproblems. To meet the accuracy requirement, we choose Δt = 1/max i N x i for the AcIIF approach. We use multiple cores to run the simulation and then combine all the results. We also implement RK2 to solve each problem, and we observe that, to maintain convergence, for a certain subproblem, e.g., N x 1 = N x 2 = …. = N x d −1 = 2 and N x d = 2 K+1−d , a very small time step must be used. As a result, RK2 requires significantly more time for these problems than the AcIIF approach, suggesting that AcIIF is a better approach.

Conclusions and Discussions
The sparse grid method is one of the major approaches used to reduce the computational cost associated with spatial discretization in solving PDEs of high dimensions. For a timedependent problem, the size of the time step, dictated by the choice of temporal integrator, becomes critical, particularly for systems with strong stiffness in reactions or large diffusion coefficients. The implicit Integration Factor (IIF), which treats the reaction implicitly and handles diffusion by exact integration, is clearly a natural choice for solving stiff PDEs in high dimensions.
In this paper, we have combined the sparse grid approach with IIF methods for a new class of methods to solve reaction-diffusion equations. We have integrated IIF methods with three different sparse grid approaches: finite element, finite difference, and a multi-level combination technique. In addition, we have employed an array compact representation of IIF (AcIIF) previously developed for uniform and regular grids to reduce the cost associated with storage and computing the exponentials of matrices in sparse grids. After studying various combinations of different sparse grid techniques and IIF methods, we have found that the approach based on AcIIF and the sparse grid combination technique is most effective for solving reaction-diffusion equations, especially for the system containing crossderivatives and non-constant coefficients.
Because the sparse grid combination technique has intrinsic distributed structures, parallel implementation becomes straightforward, in comparison with the finite difference or the finite element methods in combination with IIF. Therefore, a GPU implementation can be adopted to significantly improve the efficiency of the sparse grid combination technique integrated with IIF. So far, most of the systems we deal with have a Dirichlet boundary condition on a rectangular domain. It is possible to extend the sparse grid-based IIF method to problems with complicated boundary conditions as discussed in [37] or irregular domains. The order of accuracy in space can be improved by introducing higher-order polynomials as the hierarchical bases, particularly for problems with high-order or mixed derivatives. The integration factor methods on sparse grids can likely be used in broad applications involving PDEs in high dimensions such as Fokker-Planck equations in biology and Black-Scholes equations in finance.      The CPU time (in seconds) of the AcIIF2-SG-FD, IIF2-SG-FD and RK2-SG-FD methods for different spatial mesh sizes N  Table 2 The CPU time (in seconds) for AcIIF2-SG-FD and RK2-SG-FD for different spatial mesh sizes N x = 2 3 , 2 4 , 2 5 , 2 6 , 2 7 for d = 2, 3, 4 respectively. For AcIIF2-SG-FD, the time step is Δt = 1/N x . In RK2-SG-FD, the time step is . All simulations end at t = 0.5.