**4. EPI methods**

**3.3. Non-self-adjoint operators**

12 Applied Linear Algebra in Action

**Arnoldi iteration** *v*1 = **z**0/∥**z**0 ∥ <sup>2</sup> **for** *j* = 1, 2, … **z***<sup>j</sup>* = *A***v***<sup>j</sup>* **for** *k* = 1, 2, …, *j*

*hkj* =**v***<sup>k</sup>*

 = **z***<sup>j</sup>* − *hkj* **v***<sup>k</sup>*

 = ∥ **z***<sup>j</sup>* ∥ <sup>2</sup>

which we obtain

/*hj* + 1,*<sup>j</sup>*

thonormal columns, such that

1

w

**z***<sup>j</sup>*

**v***j* + 1 = **z***<sup>j</sup>*

**end**

 **end** *hj* + 1,*<sup>j</sup>* *<sup>H</sup>* **<sup>z</sup>** *<sup>j</sup>*

The theory developed in [12] applies to symmetric positive definite matrices, but this prop‐ erty 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 in‐ tegrals, as described in [19], can suffer from "serious breakdown" [20]. Therefore, a more ro‐ bust 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

The output of Arnoldi iteration is an upper Hessenberg matrix *Hm*, and a matrix *Vm* with or‐

= . 1, 1

By analogy with (17), to approximate **u***<sup>n</sup>* + 1 = *φ*(*Aτ*)**u***<sup>n</sup>*, we can compute each Fourier coefficient [**û***<sup>n</sup>* + 1]*ω* of **u***<sup>n</sup>* + 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 *Hm*(*ω*), from

> wt

where *B*0(*ω*) and *E*<sup>12</sup> are as defined in (15) and (17), respectively. As shown in [6, 7], like with block Lanczos, the eigenvalue problem for *Hm*(*ω*) approximately decouples for high frequen‐

*m*

0 12 12 0 <sup>12</sup> [ ] = ( ) ( ( )) ( ) , <sup>ˆ</sup> *<sup>n</sup> H H*

<sup>+</sup> é ù

*B E H EB*

 j

w

*<sup>H</sup> AV V H h m m m m mm m* + + + **z e** (21)

 w

ë û **u** (22)

Arnoldi iteration, applied to a matrix *A* and initial vector **z**0, is as follows:

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 *tn* to time *tn* + 1 = *tn* + *Δt* is taken as follows. First, the time derivative *F*(**y**) is expressed in terms of its Taylor expansion around **y**(*tn*), which yields

$$\frac{d\mathbf{y}}{dt} = F(\mathbf{y}(t\_u)) + A\_\imath(\mathbf{y}(t) - \mathbf{y}(t\_u)) + R(\mathbf{y}(t)),\tag{23}$$

where *An* is the Jacobian matrix of **F**(**y**) evaluated at **y**(*tn*), and *R*(**y**(*t*)) is the Taylor remainder. Next, we multiply both sides by an integrating factor *e* −*Ant* and then integrate from *tn* to *tn* + 1 to obtain

$$\mathbf{y}(t\_{n+1}) = \mathbf{y}(t\_n) + [e^{t\_n \mathbf{u}} - I]A\_n^{-1}F(\mathbf{y}(t\_n)) + \int\_{t\_n}^{t\_{n+1}} e^{A\_n(t\_{n+1}-\tau)}R(\mathbf{y}(\tau))d\tau. \tag{24}$$

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* = *An*, 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

$$\|\varphi(A\tau)\mathbf{b}\approx\|\|\mathbf{b}\|\|\_{2}\,V\_{m}\varphi(H\_{m}\tau)\mathbf{e}\_{1},\tag{25}$$

where, as before, *Hm* is upper Hessenberg (or tridiagonal of *A* is symmetric), and the columns of *Vm* form an orthonormal basis for the Krylov subspace

$$\mathcal{K}(A, \mathbf{b}, m) = \text{span}\left\{ \mathbf{b}, Ab, A^2 \mathbf{b}, \dots, A^{m-1} \mathbf{b} \right\}.$$

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 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 *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 essentially noise.

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

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

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

The behavior of the unfiltered Krylov vectors is not surprising, as similar behavior is dis‐ played 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*∥*Vm TVm* <sup>−</sup> *Im*∥*F* was negligibly small (that is, on the order of 10<sup>−</sup><sup>10</sup> or smaller) in both cases.
