**5. KSS-EPI methods**

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 exa‐ cerbated by increasing the spatial resolution in the discretization of the underlying PDE

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 *Vm*, 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 *Vm* to

**Figure 2.** Columns of *Vm* from (25) generated by Arnoldi iteration applied to the matrix from Burgers' equation (see

This can be alleviated by filtering out high-frequency components of the columns of *Vm* 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 *Vm*. Future work will include the automatic, adaptive selection of an appropriate threshold for filtering high-frequency

**Figure 3.** Columns of *Vm* from (25) generated by Arnoldi iteration, with denoising, applied to the matrix from Burgers'

from which (1) arises.

14 Applied Linear Algebra in Action

essentially noise.

Section 6.2).

components.

equation (see Section 6.2).

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:


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,*<sup>ω</sup>*, *λ*2,*<sup>ω</sup>*, …, *λ<sup>K</sup>*,*<sup>ω</sup>*}, as described in Section 3, by approximating the entries of the tridiagonal matrix T*<sup>K</sup>* (*ω*) 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*}.

The Fourier coefficient of *φ*(*Aτ*)**b***H* corresponding to the wave number *ω* is obtained by computing the same Fourier coefficient of *p*2*<sup>K</sup>* <sup>−</sup> 1,*<sup>ω</sup>*(*Aτ*)**b***H*, where *p*2*<sup>K</sup>* <sup>−</sup> 1,*<sup>ω</sup>* is the polynomial interpolant of *φ*(*λ*) with interpolation points {*λ<sup>i</sup>* , *λi*,*ω*}*i*=1 *<sup>K</sup>* . Expressing this interpolant in Newton form, we have

$$\begin{split} p\_{2K-1,\alpha}(\boldsymbol{\lambda}) &= \sum\_{j=1}^{K} g[\boldsymbol{\lambda}\_{1}, \ldots, \boldsymbol{\lambda}\_{j}] \prod\_{i=1}^{j-1} (\boldsymbol{\lambda} - \boldsymbol{\lambda}\_{i}) + \\ &\sum\_{j=1}^{K} g[\boldsymbol{\lambda}\_{1}, \ldots, \boldsymbol{\lambda}\_{K}, \boldsymbol{\lambda}\_{1,\alpha}, \ldots, \boldsymbol{\lambda}\_{j,\alpha}] \left( \prod\_{i=1}^{j-1} (\boldsymbol{\lambda} - \boldsymbol{\lambda}\_{i,\alpha}) \right) \prod\_{k=1}^{K} (\boldsymbol{\lambda} - \boldsymbol{\lambda}\_{k}). \end{split} \tag{26}$$

Arranging the interpolation points in the order indicated above allows us to reduce the number of FFTs needed. Using the relation from Lanczos iteration,

$$AX\_{\kappa} = X\_{\kappa} T\_{\kappa} + \mathbf{r}\_{\kappa} \mathbf{e}\_{K}^{\mathbb{T}},\tag{27}$$

we define

$$\begin{split} \mathbf{v} &= p\_{K-1}(A)\mathbf{b}\_{H} = \left\{\boldsymbol{\varrho}[\hat{\lambda}\_{1}] + \boldsymbol{\varrho}[\hat{\lambda}\_{1}, \hat{\lambda}\_{2}](A - \hat{\lambda}\_{1}I) + \cdots \right\} \\ &+ \boldsymbol{\varrho}[\hat{\lambda}\_{1}, \hat{\lambda}\_{2}, \dots, \hat{\lambda}\_{K}](A - \hat{\lambda}\_{1}I)\cdots(A - \hat{\lambda}\_{K-1}I) \right\} \mathbf{b}\_{H} \\ &= \| (\mathbf{b}\_{H} \parallel\_{\mathbb{L}} X\_{K} p\_{K-1}(A) \mathbf{e}\_{1}) \\ \mathbf{w} &= q\_{K}(A)\mathbf{b}\_{H} = (A - \hat{\lambda}\_{1}I)\cdots(A - \hat{\lambda}\_{K}I)\mathbf{b}\_{H} \\ &= \beta\_{1}\beta\_{2}\dots\beta\_{K-1}\mathbf{r}\_{K} \\ &= \beta\_{1}\dots\beta\_{K}X\_{K+1}\mathbf{e}\_{K+1}, \end{split}$$

and

$$
\begin{split}
\tilde{p}\_{K-1,\alpha}(\mathcal{\lambda}) &= \varphi[\mathcal{\lambda}\_{1}, \dots, \mathcal{\lambda}\_{K}, \mathcal{\lambda}\_{1,\alpha}] + \varphi[\mathcal{\lambda}\_{1}, \dots, \mathcal{\lambda}\_{2,\alpha}](\mathcal{\lambda} - \mathcal{\lambda}\_{1,\alpha}) + \dots \\ &+ \varphi[\mathcal{\lambda}\_{1}, \dots, \mathcal{\lambda}\_{K,\alpha}](\mathcal{\lambda} - \mathcal{\lambda}\_{1,\alpha}) \cdots (\mathcal{\lambda} - \mathcal{\lambda}\_{K-1,\alpha}) \\ &= C\_{K-1}^{\alpha} \mathcal{\lambda}^{K-1} + \dots + C\_{1}^{\alpha} \mathcal{\lambda} + C\_{0}^{\alpha},
\end{split}
$$

where *Cj <sup>ω</sup>*, for *j* = 0, 1, …, *K* − 1, are the coefficients of *p*˜ *<sup>K</sup>* <sup>−</sup>1,*<sup>ω</sup>* in power form, which can easily be computed by repeated application of nested multiplication to the last *K* terms of the Newton form of *p*2*K* − 1,*<sup>ω</sup>*.

Finally, using F to denote the discrete Fourier transform, we have

$$\rho(A\mathbf{\tau})\mathbf{b}\_H \approx p\_{2K-1,a}(A\mathbf{\tau})\mathbf{b}\_H = \mathbf{v} + \tilde{p}\_{K-1,a}(A)\mathbf{w} = \mathbf{v} + \mathbb{F}^1 \sum\_{j=0}^{K-1} [C\_j^a] \mathcal{F} A^j \mathbf{w},\tag{28}$$

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 *Aj* **w**, *j* = 0, 1, …, *K* − 1 are computed, the *K* FFTs can be performed in parallel.
