Matrices, Moments and Quadrature: Applications to Time-Dependent Partial Differential Equations

The numerical solution of a time-dependent PDE generally involves the solution of a stiff system of ODEs arising from spatial discretization of the PDE. There are many methods in the literature for solving such systems, such as exponential propagation iterative (EPI) methods, that rely on Krylov projection to compute matrix function-vector products. Unfortunately, as spatial resolution increases, these products require an increasing number of Krylov projection steps, thus drastically increasing computational expense. This paper describes a modification of EPI methods that uses Krylov subspace spectral (KSS) methods, to compute these matrix function-vector products. KSS methods represent a balance between the efficiency of explicit methods and the stability of implicit methods. This balance is achieved by approximating the matrix exponential with different polynomials for each Fourier coefficient of the solution. These polynomials arise from techniques due to Golub and Meurant for computing bilinear forms involving matrix functions by treating them as Riemann-Stieltjes integrals, which are then approximated using Gaussian quadrature rules. This paper describes how the nodes for the quadrature rules required by KSS methods can be estimated very rapidly through asymptotic analysis of block Lanczos iteration, thus drastically reducing computational expense without sacrificing accuracy. Numeri‐ cal experiments demonstrate that this modification causes the number of Krylov projection steps to become bounded independently of the grid size, thus dramatically improving efficiency and scalability.


Introduction
Consider an autonomous system of ordinary differential equations (ODEs) of the form 0 0 ' = ( ), ( ) = , F t y y y y (1) such as one that would naturally arise from spatial discretization of a partial differential equation (PDE). Such systems are generally stiff when the underlying PDE has significant diffusive or advective terms, and this stiffness is exacerbated by increasing the spatial resolution. Stiffness leads to significantly increased computational expense for both explicit and implicit time-stepping methods. For explicit methods, the time step is severely restricted by the CFL condition, while implicit methods require the solution of ill-conditioned systems of linear equations during each time step. Such systems are generally sparse and therefore best suited for iterative methods, but for these stiff systems of ODEs, iterative methods require many iterations or a specially developed preconditioner [1].
Exponential propagation iterative (EPI) methods, introduced by Tokman et al. [1,2], are Runge-Kutta like time-stepping methods that are designed to minimize the number of matrix function-vector products of the form w = φ(Aτ)b, where φ is a smooth function, A is a square matrix, τ is a parameter determined by the time step, and b is a column vector. The approach used by EPI methods to compute w for a given symmetric matrix A is to apply the Lanczos algorithm to A with the initial vector b, until we obtain a matrix X j with orthonormal columns and a tridiagonal matrix T j such that X j T AX j = T j . Then, we can compute the approximation where e 1 = 1 0 ⋯ 0 T . Since each column x k of the matrix X j is of the form x k = p k − 1 (A)b, where p n (A) is a polynomial of degree n in A, w j is the product of a polynomial in A of degree j − 1 and b. Since the matrix A arises from a stiff PDE, the eigenvalues of A are distributed over several orders of magnitude. As it is generally not possible to accurately approximate φ(λ) on such a large interval with a low-degree polynomial, a large number of Lanczos iterations is generally necessary in order to obtain a sufficiently accurate of w.
The difficulty that time-stepping methods have with stiffness stems from the fact that low-and high-frequency components of the solution, which change at widely varying speeds, are coupled and therefore must be evolved using a common time step. While the low-frequency components generally make the most significant contribution to the solution, the highfrequency components change most rapidly, and therefore constrain the time step. The greater the spatial resolution, the greater the bandwidth of the solution, thus constraining the time step even further.
This coupling, however, is not the only cause of the greater computational expense incurred by time-stepping methods applied to stiff systems. A key contributing factor is the use of the same polynomial or rational function to approximate all of these components of φ(Aτ)b, when such a function cannot effectively approximate φ(λτ) on a large interval except at high degree, which increases computational expense. An alternative is to use Krylov subspace spectral (KSS) methods [3,4], which employ a component-wise approach to these matrix functionvector products. In KSS methods, which are explicit, each Fourier coefficient of the solution is computed using an interpolating polynomial with frequency-dependent interpolation points. This individualized approach to computing each Fourier component yields high-order accuracy and stability like that of implicit methods, even though KSS methods themselves are explicit.
To date, KSS methods have been applied mostly to linear PDEs on d-dimensional boxes, for d = 1, 2, 3, with either periodic or homogeneous Dirichlet or Neumann boundary conditions. A first-order KSS method was applied to nonlinear diffusion equations for image processing by Guidotti et al. [5], but for that problem, a straightforward linearization was used during each time step. In order to compute solutions of nonlinear PDEs with higher-order accuracy in time, it is necessary to treat the nonlinearity more carefully. This can be accomplished by combining KSS methods with another high-order time-stepping method, such as EPI methods.
In this paper, this combination is presented, for the purpose of solving systems of ODEs of the form (1) that are obtained through spatial discretization of nonlinear PDEs defined on rectangular domains with periodic, homogeneous Dirichlet, or homogeneous Neumann boundary conditions. The proposed algorithm, first described in [6,7], uses the componentwise approach of KSS methods to more efficiently compute matrix function-vector products of the form y = φ(Aτ)b. The following features distinguish the combined approach from previous work on either EPI or KSS methods: • Instead of applying Krylov projection (e.g., see [8][9][10]) with initial vector b to compute y, it is applied only to a low-frequency projection of b, in order to avoid the larger number of iterations that Krylov projection typically incurs at higher spatial resolution. Furthermore, for advection-dominated problems, denoising is applied to the Krylov subspace basis vectors, to eliminate spurious high-frequency oscillations in these vectors which slow convergence.
• For the high-frequency portion of b, a KSS method, as described in [11], is used to apply φ(Aτ). In this latest version of KSS methods, each Fourier component of the output vector y is still approximated using its own block quadrature rule, as before. However, now the quadrature nodes (i.e., frequency-dependent interpolation points for φ) are obtained through high-frequency analysis of block Lanczos iteration, which yields simple formulas for the nodes, instead of having to obtain them by solving an eigenvalue problem [3]. This greatly improves efficiency over the approach used in [3,4], in which the nodes are obtained by explicitly applying block Lanczos iteration to A for each Fourier component, without sacrificing accuracy [11].
• Once the frequency-dependent interpolation points are obtained, it is necessary to construct and apply frequency-dependent interpolating polynomial approximations p ω (λ) of φ to A, for each wave number ω, and then the Fourier coefficients of the solution are inner products of the form e^ω H p ω (Aτ)b, where ê ω is a discretization of an appropriate Fourier basis function.
This paper provides implementation details for this task, and explains how it can be accomplished with as few Fourier transforms as possible [6,7].
The outline of the paper is as follows. Section 2 gives an overview of KSS methods. Section 3 reviews their acceleration based on asymptotic analysis of recursion coefficients, first presented in [11], and discusses extension to non-self-adjoint operators. Section 4 provides a brief overview of EPI methods and demonstrates the spurious high-frequency oscillations that can occur when using standard Krylov projection within an EPI method. Section 5 describes how KSS and EPI methods are combined. Numerical results are presented in Section 6. Concluding remarks and directions for future work are given in Section 7.

Matrices, moments, quadrature, and PDE
To provide context for the latest advancements with KSS methods, we begin with an overview of KSS methods as described in [3], applied to the 1-D parabolic PDE = 0, 0 < < 2 , > 0, = ( ( ) ) ( ) , where p(x) > 0 and q(x) ≥ 0 on [0, 2π], and q(x) is not identically zero. In KSS methods, each Fourier coefficient of the solution ũ(x, t n + 1 ) is obtained by applying the exact solution operator of the PDE to ũ(x, t n ). For a given wave number ω, such a Fourier coefficient is given by where Applied Linear Algebra in Action

Matrices, moments, and quadrature
The spatial discretization of (6) yields the bilinear form Because the Sturm-Liouville operator L is self-adjoint and positive definite, it follows that the matrix A is symmetric positive definite. As such, it has real eigenvalues b = λ 1 ≥ λ 2 ≥ ⋯ ≥ λ N = a > 0, and corresponding orthonormal eigenvectors q j , j = 1, …, N. From the spectral decomposition of A, we obtain the following representation of (6): As described by Golub and Meurant in [12], (8) can be viewed as a Riemann-Stieltjes integral where the measure α(λ) is defined by if , This allows approximation of (6) using Gaussian quadrature rules, where the nodes and weights are obtained by applying the Lanczos algorithm to A with initial vectors u and v [12].   (8) where the matrix A represents a spectral discretization of a 1-D, second-order differential operator with smooth leading coefficient (top plot) and discontinuous leading coefficient (bottom plot), where u = v is a discretization of e 2ix (solid curve) or e 64ix (dashed curve).

Block Lanczos algorithm
In the case where u ≠ v, the presence of a negative weight would destabilize the quadrature rule [14]. Alternatively, we consider the approximation of the 2 × 2 matrix integral We use the most general K-node quadrature formula, as described in [12], to get an approximation for (9) of the form  (11) where, for each j, λ j is a scalar and v j is a 2-vector. Each node λ j is an eigenvalue of the matrix 1 1 which is a block-tridiagonal matrix of order 2K. The vector v j consists of the first two elements of the corresponding normalized eigenvector. The matrices M j and B j are computed using the block Lanczos algorithm [15], described below.

Krylov subspace spectral methods
The block KSS method [3,4] for (3) begins by defining where ê ω is a discretization of 1 2π e iωx and u n is a discretization of the approximate solution u(x, t) at time t n = nΔt. Next, we compute the QR factorization of R 0 , and where Then we apply the block Lanczos algorithm [15] to the matrix L N , which comes from the discretization of L, with initial block X 1 (ω). This produces a block tridiagonal matrix T K of the form (12), where every entry of T K is a function of ω. Then, at time t n + 1 , each Fourier coefficient of the solution is [ ] An inverse FFT applied to the vector of Fourier coefficients û n + 1 yields the vector û n + 1 , which consists of the values of the solution u(x, t n + 1 ) at the grid points.
This algorithm has temporal accuracy O(Δt 2K − 1 ) for parabolic problems [3]. Even higher-order accuracy, O(Δt 4K − 2 ), is obtained for the second-order wave equation [4], due to the secondorder time derivative. Furthermore, under appropriate assumptions on the coefficients of the PDE, the 1-node KSS method is unconditionally stable [3,4]. More generally, the temporal order of accuracy is O(Δt (2K − 1)d ), where d is the highest order of a time derivative in the PDE; this order of accuracy has been observed with the Schrödinger equation [16] and Maxwell's equations [18].

Asymptotic analysis of block Lanczos iteration
For each Fourier coefficient of the solution, KSS methods use a different quadrature rule, because the measure α(λ) in (8) is defined in terms of the frequency. It follows that if S(L; Δt) is the solution operator of the PDE (e.g. S(L; Δt) = e − LΔt for the PDE (3)), then the function S(λ; Δt) is approximated by a polynomial of degree 2K that interpolates S(λ; Δt) at different points for each frequency, which are the block Gaussian quadrature nodes λ j in (11). Taking all Fourier coefficients together, the computed solution at time t n + 1 , u n + 1 can be expressed as where L N is a discretization of L on an N-point grid and D j (Δt) is a matrix that is diagonal in discrete Fourier space. The diagonal entries are the coefficients of these frequency-dependent interpolating polynomials in power form, with each row corresponding to a different frequency.
In the block KSS method described in [3,4] and reviewed in the previous section, these interpolation points are the eigenvalues of the block tridiagonal matrix T K from (12) that is produced by block Lanczos iteration. Although each such matrix is small, computing the eigenvalues and eigenvectors for each frequency is computationally expensive. Therefore, in this section, we describe a much faster approach to obtaining estimates of these nodes, at least for high frequencies. This approach was first presented in [11] and generalized in [6,7]. The basic idea is to examine the entries of T K as |ω| → ∞, where ω is the wave number.

The block case
As in the previous section, let u n be a discretization of the approximate solution u(x, t) at time t n = nΔt on a uniform N-point grid. Then, KSS methods use the initial block R 0 = ê ω u n , for each ω = − N/2 + 1, …, N/2. We start the first iteration of the block Lanczos algorithm by finding the QR-factorization of R 0 : with u ω n is defined as in (16). We note that if the solution u is continuous, then as |ω| → ∞, | û n (ω)| → 0, so that in the limit B 0 is diagonal.
The next step is to compute where the matrix L N is a spectral discretization of the operator L defined by Lu = pu xx + q(x)u, with p being a constant. Substituting the value of X 1 from (18) into (19) yields Rayleigh quotient of L N and u ω n . As |ω| increases, the Fourier coefficients of a continuous function decay to zero; therefore, as long as the solution is sufficiently regular, the non-diagonal entries of M 1 become negligible, that is, Proceeding with the iteration, and neglecting any terms that are Fourier coefficients or are of lower order in ω, we obtain where q is a vector consisting of the value of q(x) at the grid points, q = q − q , and multiplication of vectors is component-wise.
To obtain X 2 , we perform the QR-factorization R 1 = X 2 B 1 . We note that the (1, 2) entry of B 1 , modulo lower-order terms, is the Fourier It follows that if the coefficient q(x) and solution u(x, t) are sufficiently regular, then B 1 approaches a diagonal matrix as |ω| → ∞, just as B 0 does. Continuing this process, it can be shown that given sufficient regularity of the solution and coefficients of L, each block M j or B j of T K from (12) becomes approximately diagonal at high frequencies.
Because T K ij ≈ 0 when i + j is odd in the high-frequency case, it follows that if the rows and columns of T K are permuted so that odd-numbered and even-numbered rows and columns are grouped together and in order, then the eigenvalue problem for T K approximately decouples. Therefore, approximate eigenvalues of T K can be obtained by computing the eigenvalues of the tridiagonal matrices obtained by performing "non-block" Lanczos on L N with initial vectors equal to the two columns of R 0 separately, rather than applying block Lanczos with these two columns together in the initial block. In [11], it is shown that this decoupling also takes place if the leading coefficient p(x) of L is not constant.

The non-block case
The decoupling observed in the preceding discussion reveals that we can obtain approximations of half of the block Gaussian quadrature nodes for all Fourier coefficients by applying "non-block" Lanczos iteration to the matrix L N with initial vector u n , the computed solution, as is done in Krylov projection methods such as those described in [8][9][10]. These nodes, denoted by λ 1 , …, λ m where m is the number of iterations, will be referred to as frequency-independent nodes. Because this iteration does not depend on the wave number ω, the frequency-independent nodes need only be computed once for each vector u for which an expression of the form φ(L N τ)u n is required. The other half of the nodes can be estimated through an asymptotic analysis of Lanczos iteration applied to L N with initial vector ê ω [11]; these are called frequencydependent nodes and will be denoted by λ 1,ω , …, λ m,ω .
As an example, we consider the case where the matrix A comes from a spectral discretization of the operator Lu = − pu xx + q(x)u, where p is a constant, and assuming periodic boundary conditions [11]. Carrying out three iterations, which corresponds to a fifth-order accurate KSS method for a parabolic PDE, yields the following recursion coefficients as functions of the wave number ω, after neglecting lower-order terms: It follows that the frequency-dependent nodes can easily be estimated as whereas for a third-order KSS method, the frequency-dependent nodes are In [6,7], similar formulas for the nodes are derived for a PDE with homogeneous Neumann boundary conditions, and for a 2-D PDE with periodic boundary conditions.
When the matrix A is a finite-difference representation of the underlying differential operator, the block Gaussian quadrature nodes can be represented more accurately if formulas for the eigenvalues of symmetric Toeplitz matrices are used for the leading-order terms in the nodes. For example, in (20), pω 2 is replaced by 2p(N/π) 2 (1 − cos(πω/N)), where N is the number of grid points.

Non-self-adjoint operators
The theory developed in [12] applies to symmetric positive definite matrices, but this property is not essential [17,18]. That said, care must be taken with nonsymmetric matrices, as a straightforward use of unsymmetric Lanczos to treat bilinear forms as Riemann-Stiejtjes integrals, as described in [19], can suffer from "serious breakdown" [20]. Therefore, a more robust approach for applying KSS methods to linear PDE with non-self-adjoint spatial differential operators is to use Arnoldi iteration to approximate φ(Aτ)b. The algorithm for Arnoldi iteration, applied to a matrix A and initial vector z 0 , is as follows: The output of Arnoldi iteration is an upper Hessenberg matrix H m , and a matrix V m with orthonormal columns, such that By analogy with (17), to approximate u n + 1 = φ(Aτ)u n , we can compute each Fourier coefficient [û n + 1 ] ω of u n + 1 by applying block Arnoldi iteration [21] to A, with initial block R 0 (ω) as defined in (13). After m iterations, this iteration produces a block upper Hessenberg matrix H m (ω), from which we obtain where B 0 (ω) and E 12 are as defined in (15) and (17), respectively. As shown in [6,7], like with block Lanczos, the eigenvalue problem for H m (ω) approximately decouples for high frequen-cies, due to the decay of the Fourier coefficients of u n . It follows that we can easily estimate the frequency-dependent eigenvalues of H m (ω), which are used as interpolation points for a polynomial approximation of φ(λ), by applying "non-block" Arnoldi iteration, as described in the above algorithm, with an initial vector z 0 chosen to be a discretization of an appropriate Fourier basis function. Additional details can be found in [6,7].

EPI methods
In this section, we give an overview of EPI methods, as developed by Tokman et al. in [1,2].
To solve an autonomous system of ODEs of the form (1), a time step from time t n to time t n + 1 = t n + Δt is taken as follows. First, the time derivative F(y) is expressed in terms of its Taylor expansion around y(t n ), which yields where A n is the Jacobian matrix of F(y) evaluated at y(t n ), and R(y(t)) is the Taylor remainder.
Next, we multiply both sides by an integrating factor e −A n t and then integrate from t n to t n + 1 to The integral on the right side is then approximated using a quadrature rule. This requires the evaluation of matrix function-vector products of the form φ(Aτ)b, with A = A n , for various "φfunctions" and various choices of the vectors b and scaling factors τ, which are determined so as to satisfy order conditions.
Any such matrix function-vector product is computed using Krylov projection. Arnoldi iteration is applied to A (or Lanczos iteration, if A is symmetric) with initial vector b. After m iterations, we obtain (21), from which we obtain an approximation where, as before, H m is upper Hessenberg (or tridiagonal of A is symmetric), and the columns of V m form an orthonormal basis for the Krylov subspace The accuracy of this approximation is discussed in [9]. In the case where A is ill-conditioned, the number of iterations m needed for convergence of (25) can be quite large, and this is exacerbated by increasing the spatial resolution in the discretization of the underlying PDE from which (1) arises.
When the number of iterations is large, an additional issue, particularly for advectiondominated problems, is the appearance of spurious high-frequency oscillations in the columns of V m , even if the initial vector b represents a smooth function. As can be seen in Figure 2, even after relatively few iterations, these oscillations can occur and reduce the columns of V m to essentially noise. This can be alleviated by filtering out high-frequency components of the columns of V m after each matrix-vector multiplication. As can be seen in Figure 3, the spurious high-frequency oscillations are reduced, resulting in more meaningful data in V m . Future work will include the automatic, adaptive selection of an appropriate threshold for filtering high-frequency components. The behavior of the unfiltered Krylov vectors is not surprising, as similar behavior is displayed by the unsmoothed Fourier method applied to hyperbolic PDEs [22]. In that work, the proposed remedies were to either use filtering, or increase the number of grid points; the former remedy serves as the motivation for denoising in this context. It is worth noting that in both the unfiltered and filtered cases, no loss of orthogonality was observed in the Krylov vectors; that is, 1 / m ∥ V m T V m − I m ∥ F was negligibly small (that is, on the order of 10 − 10 or smaller) in both cases.

KSS-EPI methods
We now describe the combination of KSS and EPI methods. The EPI method itself is not changed; what is modified is the approach to computing any matrix function-vector product of the form φ(Aτ)b. The main steps in the new approach are as follows: 1. First, it is necessary to determine a cutoff frequency N c . In the numerical experiments presented in this paper, the value of N c has been determined by experimentation; it is demonstrated in [7] that the performance is not unduly sensitive to the choice of N c . In future work, an adaptive approach to choosing N c will be developed.

Using an FFT, b is decomposed into low-frequency and high-frequency components
3. φ(Aτ)b L is computed using Krylov projection, as is normally done in EPI methods.

φ(Aτ)b H
is computed using a KSS method, as described in Section 3.

5.
The results of the preceding two steps are added to obtain an approximation of φ(Aτ)b.
We now provide details on how step 4 can be performed efficiently, by minimizing the number of FFTs. To simplify the exposition, we consider the 1-D case, with periodic boundary conditions. For each wave number ω, we obtain the frequency-dependent nodes {λ 1,ω , λ 2,ω , …, λ K,ω }, as described in Section 3, by approximating the entries of the tridiagonal matrix T K (ω) produced by Lanczos iteration with initial vector ê ω . As demonstrated in Section 3, this is readily accomplished using the coefficients of the spatial differential operator. Next, we use Krylov projection with initial vector b H , as in Section 4, and compute the eigenvalues of the resulting tridiagonal (or Hessenberg, if A is nonsymmetric) matrix. These are the frequency-independent nodes {λ 1 , λ 2 , …, λ K }.
Arranging the interpolation points in the order indicated above allows us to reduce the number of FFTs needed. Using the relation from Lanczos iteration, = ,  Finally, using F to denote the discrete Fourier transform, we have and it can easily be seen that the solution at each time step requires K FFTs and one inverse FFT. It is worth noting that once the Krylov subspace vectors A j w, j = 0, 1, …, K − 1 are computed, the K FFTs can be performed in parallel.

Numerical results
In this section, we compare several versions of EPI methods, as applied to two test problems; additional test problems are featured in [6,7]. The versions differ in the way in which they compute matrix function-vector products of the form φ(Aτ)b: • standard Krylov projection, as in (25), hereafter referred to as "Krylov-EPI", either with or without denoising as in Section 4; • using the KSS approach, as described in Section 5, hereafter referred to as "KSS-EPI"; • Newton interpolation using Leja points [23], hereafter referred to as "LEJA"; and • adaptive Krylov projection [24], hereafter referred to as "AKP".
All of these approaches are used in the context of two EPI methods. The first is a third-order, two-stage EPI method [1] where a 11 = 9/4 and b 1 = 32/81, and For this method,  å and the coefficients g ij , a ij , b j , and p ij are obtained from the description of the EPIRK5s3 method in [25].
For all methods used to compute matrix function-vector products, efficiency will be measured in two ways: (1) total execution time required to integrate over the entire time interval, and (2) the average number of matrix-vector products (plus the average number of FFTs, in the case of KSS-EPI) performed for each evaluation of a matrix function-vector product of the form φ(Aτ)b. The phrase "number of iterations" is used throughout this section to refer to this quantity. Throughout this section, for the purpose of discussing performance as a function of spatial resolution, N refers to the number of grid points per dimension.

Diffusive problem
For our first test problem, we choose a diffusion-dominated PDE: the 2-D Allen-Cahn equation, with α = 0.1. We impose homogeneous Neumann boundary conditions, and the initial condition  (2 ). u x y x y p p + The Laplacian is discretized using a centered finite difference. For KSS-EPI, the low-frequency portion b L consists of all components with wave numbers ω i ≤ 7, i = 1, 2. That is, for this problem, the value of N c , as defined in the previous section, is 7. The low value of this threshold is due to the smoothness of the initial data.  The results are shown in Figures 4 and 5 for the third-order EPI method (29), and in Figures 6 and 7 for the firth-order EPI method (30). It can be seen from Figures 4 and 6 that for both EPI methods, the accuracy of all four approaches to computing matrix function-vector products are comparable. However, the efficiency and scalability of these four approaches are very different. In particular, these figures show that for KSS-EPI methods, the growth in execution time as a function of the number of grid points per dimension is much less. For example, note that for the coarsest grid, with N = 50, the speed of KSS-EPI is similar to that of Krylov-EPI, but for N = 150 and N = 300, KSS-EPI is much faster. Furthermore, as N increases, this gap in performance grows. This is due to the fact that while both methods use Krylov projection, KSS-EPI applies it to an initial vector with only low-frequency components, thus significantly reducing the number of iterations needed for convergence. The difference in scalability is more clearly illustrated in Figures 5 and 7. It can be seen that for KSS-EPI, the number of overall iterations (matrix-vector multiplications + FFTs) shows almost no sensitivity to the grid size, compared to Krylov-EPI, AKP and Leja interpolation, all of which exhibit substantial growth as the number of grid points increases. It can also be seen that denoising results in a slightly higher number of overall iterations, so it is not beneficial for this problem. This is not surprising, as this is a diffusive problem that has a smooth solution and is therefore less susceptible to high-frequency oscillations during Lanczos iteration. Although AKP is slightly more accurate than KSS-EPI for fifth order, this is more than offset by the superior efficiency of KSS-EPI with denoising, particularly at larger grid sizes and larger time steps.

Advective problem
For our second test problem, we choose an advection-dominated PDE: a 1-D Burgers' equation with ν = 0.03. We impose homogeneous Dirichlet boundary conditions, and the initial condition  For KSS-EPI, the low-frequency portion b L consists of all components with wave numbers ω ≤ N c = 40. The higher threshold for this problem, compared to the Allen-Cahn equation (31), is due to the fact that the initial data is less smooth.
The results are shown in Figures 8 and 9 for the third-order EPI method (29), and in Figures 10 and 11 for the fifth-order EPI method (30). From Figures 8 and 10, it can be seen that Krylov-EPI and KSS-EPI have comparable accuracy, but even for the coarsest grid with N = 500 points, KSS-EPI is much faster, and as with the Allen-Cahn equation, this gap in performance only grows with N. As shown in Figures 9 and 11, this is again due to the increasing number of Krylov projection steps needed for Krylov-EPI. While Leja interpolation is more accurate than KSS-EPI for both the third-and fifth-order EPI methods, and AKP is also more accurate in the fifth-order case, both of these methods are still significantly slower than KSS-EPI, and a similar scalability gap can also be observed.  The difference in scalability among the four approaches to computing matrix function-vector products is more apparent in Figures 9 and 11. Denoising applied to KSS-EPI is advantageous for this problem, unlike with the Allen-Cahn equation. As can be seen in these figures, for KSS-EPI without denoising, the number of iterations is increasing with the number of grid points, though not nearly as rapidly as with the other methods. However, with denoising included, the insensitivity of the number of iterations to the grid size is restored.

Discussion of efficiency
The major components of the computational cost of KSS-EPI stem from Krylov projection that is applied to low-frequency parts, and FFTs that are applied to both the low-and highfrequency parts. Specifically, suppose the EPI method is of order p, and q Krylov projection steps are needed for convergence of the low-frequency part. Then, the task of evaluating φ(Aτ)b requires q + p matrix-vector multiplications, (p + 3)/2 FFTs and 2 inverse FFTs if denoising is not used, and q + p matrix-vector multiplications, q + (p + 3)/2 FFTs and q + 2 inverse FFTs if denoising is used. Denoising, therefore, is only worthwhile if the value of q can be substantially reduced, as in the case of Burgers' equation.

Conclusions and future work
We have demonstrated that when solving stiff systems of nonlinear ODEs derived from PDEs, the growth in the computational cost that results from an increase in the number of grid points can be significantly reduced by performing a relatively low and grid-insensitive number of Krylov projection steps and FFTs on low-and high-frequency portions of the solution separately, instead of a number of Krylov projection steps on the entire solution that grows substantially with the number of grid points. The component-wise approach employed by KSS methods, in which each Fourier coefficient of a matrix function-vector product φ(τA)b is computed using a frequency-dependent interpolant of φ, allows the Krylov subspace dimension to be bounded independently of the grid size and instead determined by the desired temporal order of accuracy (for the high-frequency part) and by the number of Krylov projection steps required on a coarse grid (for the low-frequency part).
Future work on KSS-EPI methods will focus on the computation of low-frequency Fourier components of the solution, as this step accounts for the most significant part of the running time. This work requires an efficient approach to automatically selecting the threshold N c that is used to distinguish low-frequency from high-frequency components. Such adaptivity can be based on the smoothness of the solution and the performance of Krylov projection during previous time steps. In addition, computing low-frequency components with methods other than Krylov projection, including Leja interpolation and adaptive Krylov projection, will be investigated. Such combination requires examination of error estimation and stopping criteria for these methods, to determine whether their convergence can be accelerated if it is known that the initial vector represents a bandlimited function.
It is important to note that the decomposition of the block Gaussian quadrature nodes into frequency-dependent and frequency-independent nodes is not limited to cases in which Fourier decomposition is used. In [26], the same idea is used for PDEs defined on a disk, in which a Legendre polynomial expansion is used for the radial part of the solution. The decoupling of the eigenvalue problem for T K first observed in [11], and further exploited in [6,7], occurs due to the rapid decay of the coefficients in whatever series expansion is used to represent the solution. Future work will include taking advantage of this behavior in order to generalize KSS methods to PDEs defined on non-rectangular domains.