# Spectral analysis of implicit 2 stage block Runge-Kutta preconditioners

Martin J Gander, Michal Outrata

<figure><img src="image_1.png" /></figure>

To cite this version:

Martin J Gander, Michal Outrata. Spectral analysis of implicit 2 stage block Runge-Kutta preconditioners. 2023. hal-04016604

HAL Id: hal-04016604

https://hal.science/hal-04016604

Preprint submitted on 6 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spectral analysis of implicit 2 stage block Runge-Kutta preconditioners

Martin J. Gander*, Michal Outrata†

## Abstract

We analyze the recently introduced family of preconditioners in [13] for the stage equations of implicit Runge-Kutta methods for two stage methods. We give explicit formulas for the eigenvalues and eigenvectors of the preconditioned systems for a general method and use these to give explicit convergence estimates of preconditioned GMRES for some common choices of the implicit Runge-Kutta methods. This analysis also allows us to qualitatively predict and explain the main observed features of the GMRES convergence behavior, not only bound it. We illustrate our analysis with numerical experiments. We also consider the direction of numerical optimization for improving the preconditioners performance, as suggested in [13]. We consider two different ways – both distinct to the one introduced in [13] – and numerically optimize these, using the explicit bounds obtained beforehand.

**Keywords:** implicit Runge-Kutta methods, two stages, preconditioners, GMRES, bounds

**Classification:** 65F08, 65F10

## 1 Introduction

Runge-Kutta methods are a well-established family of one-step solvers for systems of ordinary differential equations (ODEs; see [20, 19] for an overview and further references). For implicit methods (IRK), their efficiency relies on having an efficient solver for the so-called *stage equations* – in general a system of *ms* non-linear equations, where *m* is the number of scalar ODEs in the system and *s* is the number of stages of the Runge-Kutta method. An important application arises from the space discretization of time-dependent partial differential equations (PDEs), resulting in a system of ODEs with *very* large *m*. If the spatial operator is *linear*, then the stage equations also become a system of linear algebraic equations, which are often solved by an iterative solver, e.g., a Krylov method. In [13], the authors introduced a family of preconditioners for GMRES for the stage equations, numerically showing that these preconditioners give an *outstanding* performance, especially under refinement of the spatial mesh, i.e., as *m* grows. Recently, there has been also other contributions in the direction of preconditioning the *fully implicit* Runge-Kutta stage equations for PDEs, see [16, 15] but also [10] and [2], introducing new ideas and testing these numerically on a variety of test problems.

<footnote>*Section de Mathématiques, Université de Genève; this work was partially supported by the SNF grant number 178752.</footnote>

<footnote>†Section de Mathématiques, Université de Genève; this work was partially supported by the FCS Swiss Excellence PhD Fellowship program of the Swiss Federation.</footnote>
We focus on the setting considered by Howle et al. in [13] and analyze the convergence of the preconditioned GMRES method, giving a theoretical background for the performance observed for two stage methods. The general s-stage case will be treated in a follow-up manuscript. In Section 2 we summarize some preliminary knowledge and introduce the problem and the preconditioners followed by some general results for Kronecker-like matrices in Section 3. We then give the analysis of each of the block type of the preconditioners in Section 4.1 – 4.3 and also comment on the possibility of numerical optimization to improve the solution process in Section 4.4.

## 2 Model problem and preliminaries

As our model problem we consider the heat equation on the unit square and a time interval $(0, T_{end})$, i.e.,

$$ 
\begin{gathered}
\frac{\partial u}{\partial t} = \Delta u + f \quad \text{in } \Omega \times (0, T_{\text{end}}), \\
u = g \quad \text{on } \partial\Omega \times (0, T_{\text{end}}) \quad \text{and} \quad u = u_0 \quad \text{in } \Omega \times \{0\},
\end{gathered}
\tag{1}
$$

where $\Delta$ is the Laplace operator, $f, g, u_0$ are given functions and $\Omega$ is the unit square $\Omega := (0, 1) \times (0, 1)$. We discretize in space using finite difference scheme on an equidistant grid with $N + 1$ rows and columns and with the mesh size $h = 1/N$ as in Figure 1. The values at the interior grid points become unknown functions of time, which are governed by the system of ODEs,

$$ \frac{\partial}{\partial t} u_i(t) = \frac{u_{i-N}(t) + u_{i-1}(t) - 4u_i(t) + u_{i+1}(t) + u_{i+N}(t)}{h^2} + b_i^{(ST)}(t), \tag{2} $$

for $i = N + 1, \dots, N(N - 1) - 1$, where $b_i^{(ST)}(t)$ collects the known values from the source terms, given by $g$ and $f$, at the given point. Combining the unknowns in each grid column into one vector denoted by $\mathbf{u}_k(t)$, i.e.,

$$ \mathbf{u}_k(t) := [u_{Nk+2} \quad u_{Nk+3} \quad \cdots \quad u_{N(k+1)-1}]^T(t), \quad \mathbf{u}(t) := [\mathbf{u}_1(t) \quad \cdots \quad \mathbf{u}_{N-1}(t)]^T, $$

and also analogously for $\mathbf{b}_k(t)$ and $\mathbf{b}(t)$, we rewrite (2) as

$$ \frac{\partial}{\partial t} \mathbf{u}(t) = \frac{1}{h^2} L \mathbf{u}(t) + \mathbf{b}^{(\mathrm{ST})}(t), \tag{3} $$

with

$$ L = \begin{bmatrix} T & I \\ I & \ddots & \ddots \\ \ddots & \ddots & I \\ I & T \end{bmatrix}, \quad T = \begin{bmatrix} -4 & 1 & & \\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & -4 \end{bmatrix}, \quad I = \begin{bmatrix} 1 & & & \\ & \ddots & & \\ & & \ddots & 1 \\ & & & T \end{bmatrix}, \tag{4} $$

where $L$ is of dimension $n := (N - 1)^2$, the blocks $T, I$ are of dimension $N - 1$. We discretize $[0, T_{end}]$ with $M_{T_{end}} + 1$ equidistant time points with time step $\tau = T_{end}/M_{T_{end}}$, i.e.,

$$ \{0 = t_0 < t_1 < \cdots < t_{M_{T_{\text{end}}} - 1} < t_{M_{T_{\text{end}}}} = T_{\text{end}}\}, \quad \tau = \frac{T_{\text{end}}}{M_{T_{\text{end}}}} \quad \text{and} \quad t_m = \tau \cdot m, \ m = 0, \dots, M_{T_{\text{end}}}. $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Left: grid points for $N + 1 = 4$; right: lexicographical ordering of the unknowns for $N + 1 = 4$.</figcaption></figure>

Having a *Butcher tableau*

$$ \begin{array}{c|c} \mathbf{c} & A \\ \hline \mathbf{b} & := \begin{array}{c|c} c_1 & \begin{array}{cccc} a_{1,1} & \dots & a_{1,s} \end{array} \\ \vdots & \vdots \vdots \vdots \\ c_s & \begin{array}{cccc} a_{s,1} & \dots & a_{s,s} \end{array} \\ \hline b_1 & \dots \\ \end{array} \end{array}, \qquad (5) $$

the corresponding IRK method applied to (3) at the m-th time step gives the approximation $\mathbf{u}^m \approx \mathbf{u}(t_m)$ as

$$ \mathbf{u}^m = \mathbf{u}^{m-1} + \tau \sum_{i=1}^{s} b_i \mathbf{k}_i^m, \qquad (6) $$

where the vectors $\mathbf{k}_1^m, \dots, \mathbf{k}_s^m \in \mathbb{R}^n$ are the solutions of the linear system

$$ \underbrace{\left( \begin{bmatrix} I & & & \\ & \ddots & & \\ & & a_{1,s}L & \\ & & & a_{s,1}L \\ & & & a_{s,s}L \end{bmatrix} - \frac{\tau}{h^2} \begin{bmatrix} a_{1,1}L & \dots & a_{1,s}L \\ \vdots & \ddots & \vdots \\ a_{s,1}L & \dots & a_{s,s}L \end{bmatrix} \right)}_{\equiv I_s \otimes I_n - \frac{\tau}{h^2}(A \otimes L) := M} \mathbf{k}^m = \begin{bmatrix} \frac{1}{h^2} L \mathbf{u}^{m-1} + \mathbf{b}^{(\text{ST})} (t_{m-1} + c_1 \tau) \\ \vdots \\ \frac{1}{h^2} L \mathbf{u}^{m-1} + \mathbf{b}^{(\text{ST})} (t_{m-1} + c_s \tau) \end{bmatrix}, \quad (7) $$

with

$$ \mathbf{k}^m := [\mathbf{k}_1^m \quad \dots \quad \mathbf{k}_s^m]^T \in \mathbb{R}^{ns}. $$

The symbol $\otimes$ stands for the Kronecker product (see [18] and references therein) and we would like to note here that (7) can be reformulated into a *matrix equation*, which is in general better suited for using a Krylov solver (see [12]). Here we focus on the analysis of the results in [13] and thus we do not address this any further but a study of the preconditioners from [13] in the matrix equations setting seems worthwhile. Having $p \le 2s$ as the order of convergence of the IRK method we assume that it is balanced with the spatial discretization error, i.e., that $h^2 = C_e \tau^p$ for some $C_e > 0$.

The problem (7) with the system matrix $M$ sparse can be very large for $h$ (and $\tau$) small, suggesting an iterative solver such as GMRES, BiCG or GCR should be used which in turn usually requires
a preconditioner to become truly efficient. In [13], the authors introduce the block preconditioners

$$P^d = I_s \otimes I_n - \frac{\tau}{h^2} \mathrm{diag}(A) \otimes L, \qquad (8)$$

$$P^u = I_s \otimes I_n - \frac{\tau}{h^2} D_A U_A \otimes L \quad \text{and} \quad P^l = I_s \otimes I_n - \frac{\tau}{h^2} L_A D_A \otimes L,$$

where $L_A, D_A, U_A$ are the LDU factors of the Butcher tableau matrix A. In addition, the authors also consider the block triangular preconditioners

$$P^{\mathrm{GSL}} = I_s \otimes I_n - \frac{\tau}{h^2} A_L \otimes L \quad \text{and} \quad P^{\mathrm{GSU}} = I_s \otimes I_n - \frac{\tau}{h^2} A_U \otimes L, \qquad (9)$$

where GSL/GSU stands for Gauss-Seidel lower/upper, and $A_{L,U}$ is the lower/upper triangular part of A, i.e.,

$$(A_L)_{ij} = \begin{cases} a_{ij} & \text{if } i \ge j \\ 0 & \text{otherwise} \end{cases}, \quad (A_U)_{ij} = \begin{cases} a_{ij} & \text{if } i \le j \\ 0 & \text{otherwise} \end{cases}.$$

Some of these were considered already in [17]. Notice that if $a_{ii} > 0$ for all $i = 1, \dots, s$, then the preconditioners are invertible as L is symmetric negative-definite. More general conditions for non-singularity of the preconditioners can be also derived analogously to [16, Lemma 1].

Using GMRES for a linear system $Cx = f$ with C being diagonalizable, i.e., $C = S\Lambda S^{-1}$ and $\Lambda = \mathrm{diag}(\lambda_1, \dots, \lambda_d)$, a standard convergence bound for the residuals $r_\ell$ reads

$$\frac{\|r_\ell\|}{\|r_0\|} \le \kappa(S) \min_{\substack{\varphi(0)=1 \\ \deg(\varphi) \le \ell}} \max_{1 \le i \le d} |\varphi(\lambda_i)|, \qquad (10)$$

where $\kappa(S)$ is the 2-norm condition number of the matrix S, see, e.g., [8, Section 5.7.2]. We would like to highlight some aspects of this bound that is often used to study GMRES convergence behavior.

**Remark 1.** As indicated above, the spectral information of the system matrix in GMRES (in our case of the preconditioned system) does not generally govern the convergence (see [5], [6] and [1] and also [8, Chapter 2 and 5.7] and the references therein). If the system matrix is normal, i.e., it is diagonalizable with S unitary, then the spectral information is enough to evaluate the ideal GMRES bound. However, if C is non-normal, then a convincing argument needs to be put forward to validate linking spectral information with the convergence behavior of GMRES as the authors in [8, p. 303, Remark 1] point out.

Moreover, particular knowledge of the interaction of S and the initial residual $r_0$ can lead to a qualitative and quantitative improvement on (10), see, e.g., [7]. However, studying GMRES behavior with the bound (10), this interaction is completely lost.

## 3 Analysis of the block preconditioners

We start by transforming the calculations into the eigenbasis of the spatial operator. Denoting the eigenpairs of L by $(\lambda_k, v_k)$, we organize the eigenvectors into an n-by-n matrix V and define the block transformation matrix Q,

$$V := [v_1, \dots, v_n], \quad \text{and} \quad Q := \begin{bmatrix} V & \cdots & V \end{bmatrix} \in \mathbb{R}^{s n \times s n}. \qquad (11)$$
Transforming $M$ blockwise into the $V$ basis gives $\tilde{M} := QMQ^T$,

$$ \tilde{M} = \begin{bmatrix} I & & \\ & \ddots & \\ & & I \end{bmatrix} - \frac{\tau}{h^2} \begin{bmatrix} a_{1,1}\Lambda & \dots & a_{1,s}\Lambda \\ \vdots & \ddots & \vdots \\ a_{s,1}\Lambda & \dots & a_{s,s}\Lambda \end{bmatrix}, \qquad (12) $$

with $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_n)$. With the preconditioners proposed in (8-9) we write the spectrum of the preconditioned system as

$$ \text{sp}(MP^{-1}) = \text{sp}(Q^T MP^{-1}Q) = \text{sp}(Q^T MQQ^T P^{-1}Q) = \text{sp}(\tilde{M}\tilde{P}^{-1}), $$

where $\tilde{P} := Q^T P Q$ stands for one of the right-preconditioners $P^{\text{d,GSU,u}}$ and an analogous formulation follows also for the left-preconditioners $P^{\text{l,GSL}}$. As the preconditioners are defined blockwise as scalar multiplications of $L$ and $I$, their blockwise transformation into the eigenbasis of $L$ is a straight-forward calculation - replacing $L$ with $\Lambda$ (and keeping $I$). Next, such matrices - block matrices with each block being a square, diagonal matrix - can be permuted into classical block-diagonal matrices as the following lemma shows.

**Lemma 1.** Let $C \in \mathbb{R}^{ns \times ns}$ be a real matrix with block structure such that every block is a square diagonal matrix, i.e.,

$$ C = \begin{bmatrix} \Lambda_{11} & \dots & \Lambda_{1s} \\ \vdots & \ddots & \vdots \\ \Lambda_{s1} & \dots & \Lambda_{ss} \end{bmatrix}, \quad \text{with } \Lambda_{ij} = \text{diag}(\lambda_1^{(ij)}, \dots, \lambda_n^{(ij)}) \quad \forall i,j. \qquad (13) $$

Then there exists a permutation matrix $\Pi \in \mathbb{R}^{ns \times ns}$ such that

$$ \Pi^T C \Pi = \begin{bmatrix} C_1 & & \\ & \ddots & \\ & & C_n \end{bmatrix} \quad \text{with } C_\ell = \begin{bmatrix} \lambda_\ell^{(11)} & \dots & \lambda_\ell^{(1s)} \\ \vdots & \ddots & \vdots \\ \lambda_\ell^{(s1)} & \dots & \lambda_\ell^{(ss)} \end{bmatrix} \in \mathbb{R}^{s \times s}, \qquad (14) $$

for any $\ell = 1, \dots, n$.

Hence, $C$ is diagonalizable if and only if $C_\ell$ is diagonalizable for all $\ell = 1, \dots, n$ and if $C_\ell = V_\ell^{-1} D_\ell V_\ell$ is the eigendecomposition of $C_\ell$ with $D_\ell = \text{diag}(\mu_\ell^{(1)}, \dots, \mu_\ell^{(s)})$, then

$$ \text{sp}(C) = \bigcup_{\ell=1}^{n} \bigcup_{i=1}^{s} \mu_{\ell}^{(i)}, $$

and if $(\mu, \mathbf{v})$ is an eigenpair of some $C_\ell$, then $(\mu, \Pi^T (\mathbf{v} \otimes \mathbf{e}_\ell))$ is an eigenpair of $C$. As a result, if $C$ is diagonalizable with $C = V^{-1}DV$, then

$$ \kappa(V) = \max_{\ell=1,\dots,s} \kappa(V_\ell), \qquad (15) $$

where $\kappa(\cdot)$ is the 2-norm condition number.
Proof. Setting $E_\ell = \text{diag}(0, \dots, 0, 1, 0, \dots, 0) \in \mathbb{R}^{n \times n}$ as the matrix with the only non-zero entry being at the position $(\ell, \ell)$ with value one, we observe that

$$C = \sum_{\ell=1}^{n} C_{\ell} \otimes E_{\ell}.$$

Using the Kronecker product permutation property from [18, Eqn. (1) and below], we take $\Pi$ such that

$$\Pi^T C \Pi = \sum_{\ell=1}^{n} E_{\ell} \otimes C_{\ell} = \begin{bmatrix} C_1 & & & \\ & \ddots & & \\ & & C_n & \end{bmatrix},$$

proving the first part of the statement. The rest follows by a direct calculation and the properties of block-diagonal matrices. $\square$

Using Lemma 1 we can analyse the eigenproperties of $\tilde{M}\tilde{P}^{-1}$ directly and thereby evaluate the GMRES bound.

**Remark 2.** We note that an analogous lemma to Lemma 1 can also be formulated for non-normal matrices (replacing $Q^T$ by $Q^{-1}$). Considering the Jordan canonical (or the Schur decomposition form) of $C_\ell$, Lemma 1 can be reformulated to obtain a block upper bidiagonal (or block upper-triangular) matrix.

To shorten the notation we set

$$\theta_k := \frac{\tau}{h^2} \lambda_k \quad \text{and} \quad \Theta := \frac{\tau}{h^2} \Lambda, \qquad (16)$$

as these quantities appear always together in the computations. By a direct calculation (see [11, Appendix B.8]) we get the limit behavior of $\theta_k$ as $\tau, h \to 0$,

$$\begin{gathered} (\theta_n, \theta_1) \to (-\frac{8}{C_e}, 0), \qquad (\theta_n, \theta_1) \to (-\infty, 0), \\ \underbrace{(\theta_1^{-1}, \theta_n^{-1}) \to (-\infty, -\frac{C_e}{8})}_{(\text{LIM})_{p=1}} \qquad \underbrace{(\theta_1^{-1}, \theta_n^{-1}) \to (-\infty, 0)}_{(\text{LIM})_{p>1}} \end{gathered} \qquad (17)$$

and continue by explicit calculations for the two stage methods.

# 4 Two stage methods

For the case $s = 2$ we have

$$\begin{align} \tilde{P}^d &= \begin{bmatrix} I - a_{11}\Theta & 0 \\ 0 & I - a_{22}\Theta \end{bmatrix}, & \tilde{P}^{\text{GSL}} &= \begin{bmatrix} I - a_{11}\Theta & 0 \\ -a_{21}\Theta & I - a_{22}\Theta \end{bmatrix}, & \tilde{P}^{\text{GSU}} &= \begin{bmatrix} I - a_{11}\Theta & -a_{12}\Theta \\ 0 & I - a_{22}\Theta \end{bmatrix}, \notag \\ \tilde{P}^u &= \begin{bmatrix} I - a_{11}\Theta & -a_{12}\Theta \\ 0 & I - \left(a_{22} - \frac{a_{21}a_{12}}{a_{11}}\right)\Theta \end{bmatrix}, & \tilde{P}^1 &= \begin{bmatrix} I - a_{11}\Theta & 0 \\ -a_{21}\Theta & I - \left(a_{22} - \frac{a_{21}a_{12}}{a_{11}}\right)\Theta \end{bmatrix}, \tag{18} \end{align}$$

and we start with a useful lemma summarizing some direct calculations.
**Lemma 2** ([11, Appendix B.8]). Let $C \in \mathbb{R}^{2 \times 2}$.

(i) The eigenvalues $\mu_{1,2}$ are given by $\mu_{1,2} = \frac{c_{11}+c_{22}}{2} \pm \frac{1}{2}\sqrt{D}$ with $D = (c_{11} - c_{22})^2 + 4c_{12}c_{21}$. In particular, if $D \neq 0$, then C is diagonalizable.

(ii-a) If $c_{12} \neq 0$ and $D \neq 0$, then C is diagonalizable and the eigenvectors $\mathbf{v}_{1,2}$ are given by $\mathbf{v}_{1,2} = \frac{1}{\|\tilde{\mathbf{v}}_{1,2}\|} \tilde{\mathbf{v}}_{1,2}$ with $\tilde{\mathbf{v}}_{1,2} = \begin{bmatrix} 1 \\ \alpha_{1,2} \end{bmatrix}$ and $\alpha_i = \begin{cases} 0 & \text{if } c_{11} = \lambda_i, \\ -\frac{c_{11}-c_{22}\mp\sqrt{D}}{2c_{12}} & \text{if } c_{11} \neq \lambda_{1,2}. \end{cases}$

(ii-b) If $c_{21} \neq 0$ and $D \neq 0$, then C is diagonalizable and the eigenvectors $\mathbf{v}_{1,2}$ are given by $\mathbf{v}_{1,2} = \frac{1}{\|\tilde{\mathbf{v}}_{1,2}\|} \tilde{\mathbf{v}}_{1,2}$ with $\tilde{\mathbf{v}}_{1,2} = \begin{bmatrix} \alpha_{1,2} \\ 1 \end{bmatrix}$ and $\alpha_i = \begin{cases} 0 & \text{if } c_{22} = \lambda_i, \\ -\frac{-c_{11}+c_{22}\mp\sqrt{D}}{2c_{21}} & \text{if } c_{22} \neq \lambda_{1,2}. \end{cases}$

(iii) If $c_{12}, c_{21} \neq 0$ and $D \neq 0$, then the condition number $\kappa(V)$ of the matrix of eigenvectors of C is given by $\kappa(V) = \sqrt{\frac{\|\tilde{\mathbf{v}}_1\|\|\tilde{\mathbf{v}}_2\| + \sqrt{(1+\overline{\alpha_1}\alpha_2)(1+\alpha_1\overline{\alpha_2})}}{\|\tilde{\mathbf{v}}_1\|\|\tilde{\mathbf{v}}_2\| - \sqrt{(1+\overline{\alpha_1}\alpha_2)(1+\alpha_1\overline{\alpha_2})}}}$, where $\alpha_i$ is given as in (ii-a).

(iv-a) If $c_{21} = 0, c_{12} \neq 0$ and $c_{11} \neq c_{22}$, then C is diagonalizable with real eigenvalues $c_{11}, c_{22}$ and eigenvectors $\mathbf{e}_1, \mathbf{v}_2$, and the formula from (iii) simplifies to $\kappa(V) = \sqrt{\frac{\|\tilde{\mathbf{v}}_2\|+1}{\|\tilde{\mathbf{v}}_2\|-1}}$, where $\mathbf{v}_1$ is the vector with non-zero entries from (ii-a), and $\alpha_1$ is its first component.

(iv-b) If $c_{12} = 0, c_{21} \neq 0$ and $c_{11} \neq c_{22}$, then C is diagonalizable with real eigenvalues $c_{11}, c_{22}$ and eigenvectors $\mathbf{v}_1, \mathbf{e}_2$, and the formula from (iii) simplifies to $\kappa(V) = \sqrt{\frac{\|\tilde{\mathbf{v}}_1\|+1}{\|\tilde{\mathbf{v}}_1\|-1}}$, where $\mathbf{v}_2$ is the vector with non-zero entries from (ii-b).

We analyze first the block diagonal preconditioners, and then continue with the block triangular ones. The calculations below give insight into the results presented in [13], e.g., give explicit formulas for the results in Figure 4.1 and 4.3, Table 4.3 and Table 5.1, 5.2, and 5.3 from [13] for $s = 2$.

## 4.1 Block diagonal preconditioner

A direct calculation gives

$$ \tilde{M}(\tilde{P}^{\text{diag}})^{-1} = \begin{bmatrix} I & -a_{12}\Theta(I - a_{22}\Theta)^{-1} \\ -a_{21}\Theta(I - a_{11}\Theta)^{-1} & I \end{bmatrix}, \quad (19) $$

and using Lemma 1 the eigen-information of the preconditioned system can be obtained from the 2-by-2 matrices

$$ X_k^d := \begin{bmatrix} 1 & -\frac{a_{12}\theta_k}{1-a_{22}\theta_k} \\ -\frac{a_{21}\theta_k}{1-a_{11}\theta_k} & 1 \end{bmatrix}. $$

We immediately notice that $X_k^d$ is diagonalizable if and only if<fnref n="1" />

$$ a_{12} = 0 \iff a_{21} = 0. \quad (20) $$

<footnote n="1">¹If $a_{12} = a_{21} = 0$, then A is in fact diagonal and hence $M = P^d$, making this case uninteresting.</footnote>
Assuming $a_{12}, a_{21} \neq 0$, we calculate the characteristic polynomial $p_{X_k^d}(\lambda)$ of $X_k^d$,

$$p_{X_k^d}(\lambda) = \lambda^2 - 2\lambda + 1 - \frac{a_{12}a_{21}\theta_k^2}{(1 - a_{11}\theta_k)(1 - a_{22}\theta_k)},$$

and therefore the eigenvalues $\xi_{1,2}^{(k)}$ of $X_k^d$ are given by

$$\xi_{1,2}^{(k)} = 1 \pm \sqrt{D_k} \quad \text{with} \quad D_k = \frac{a_{12}a_{21}}{(|\theta_k^{-1}| + a_{11})(|\theta_k^{-1}| + a_{22})}. \qquad (21)$$

We write $\xi_{1,2}^{(k)}$ as functions of $|\theta_k|^{-1}$,

$$\xi_1^{(k)} = 1 \pm \sqrt{\phi(|\theta_k|^{-1})} \quad \text{with} \quad \phi(\alpha) = \frac{a_{12}a_{21}}{(\alpha + a_{11})(\alpha + a_{22})},$$

and $\alpha \in (|\theta_1|^{-1}, |\theta_n|^{-1})$ – an interval converging towards the limit interval in (17). If $a_{11}, a_{22} \ge 0$ (e.g., Gauss, Radau or Lobatto methods), then $\xi_{1,2}^{(k)}$ lie on a line segment in $\mathbb{C}$ going through the point 1, which is either a part of the real axis (if $\text{sign}(a_{12}a_{21}) \ge 0$) or on the line $1 + \beta i$, $\beta \in \mathbb{R}$ (otherwise). We have

$$|\xi_{1,2}^{(k)} - 1| = \left| \sqrt{\phi(|\theta_k^{-1}|)} \right| = \sqrt{\left| \frac{a_{12}a_{12}}{(|\theta_k^{-1}| + a_{11})(|\theta_k^{-1}| + a_{22})} \right|},$$

and hence<fnref n="2" /> the maximum of $|\xi_{1,2}^{(k)} - 1|$ as a function of $|\theta_k^{-1}|$ is attained either at one of the endpoints of the interval in (17) or at an interior stationary point. Calculating the derivative, we get

$$(|\phi|)'(\alpha) = -\text{sign}(\phi(\alpha)) a_{12}a_{21} \frac{2\alpha + a_{11} + a_{22}}{(\alpha + a_{11})^2(\alpha + a_{22})^2},$$

and thus the only candidate for a stationary point is $-(a_{11} + a_{22})/2$ assuming it belongs to the domain of $\phi$ (as mentioned above, this is not the case for the commonly used Gauss, Radau or Lobatto methods). Assuming it does not, e.g., because $a_{11}, a_{22} > 0$, the maximum is attained at the left endpoint of the interval in (17), bounded from above by the value at $\alpha = 0$ which gives

$$|\xi_{1,2}^{(k)} - 1| \le \sqrt{\left| \frac{a_{12}a_{21}}{a_{11}a_{22}} \right|}. \qquad (22)$$

Thus making this quantity small will make the eigenvalues cluster tightly around 1. Notice that the bound above is suggesting to make the diagonal entries large compared to the off-diagonal ones, making the matrix diagonal in the limit (and hence making the preconditioner exact). Assuming $a_{ij} \neq 0$ for $i, j = 1, 2$ and $a_{11}, a_{22} \ge 0$ we use Lemma 2, and the condition number of the matrix of eigenvectors $S_k^d$ of $X_k^d$ is given by

$$\kappa(S_k^d) = \sqrt{\frac{1 + \left| \frac{a_{21}(|\theta_k|^{-1} + a_{22})}{a_{12}(|\theta_k|^{-1} + a_{11})} \right| + \left| 1 - \frac{a_{21}(|\theta_k|^{-1} + a_{22})}{a_{12}(|\theta_k|^{-1} + a_{11})} \right|}{1 + \left| \frac{a_{21}(|\theta_k|^{-1} + a_{22})}{a_{12}(|\theta_k|^{-1} + a_{11})} \right| - \left| 1 - \frac{a_{21}(|\theta_k|^{-1} + a_{22})}{a_{12}(|\theta_k|^{-1} + a_{11})} \right|}}. \qquad (23)$$

<footnote n="2">²We assumed at the beginning of the section that the inverse $(\tilde{P}^d)^{-1}$ exists and hence the denominator of $\phi$ is non-zero. Hence $|\phi(\alpha)|$ is a smooth function.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: We show the GMRES performance together with the evaluated bound (10) for different preconditioners and Butcher tableaus. We set $N = 100$, $\tau = h^{2/p}$ ($p$ is the order of convergence of the IRK method; $p_{\text{Gauss}} = 4$, $p_{\text{RadauIIA}} = 3$, $p_{\text{LobattoIIIC}} = 2$) and use a random right-hand side.</figcaption></figure>

Recalling (17), we see that $\kappa(S_k^d)$ has no singularities in it and the limits at 0 and $+\infty$ are also real and bounded. Therefore the conditioning of the eigenbasis of the preconditioned system will be uniformly bounded with respect to mesh refinement in space and time – just as the clustering diameter.

**Remark 3.** Some authors call such preconditioners, i.e., preconditioners such that the eigenproperties of the preconditioned system can be bounded independently of $h$ and $\tau$, order-optimal and the general way to show order-optimality for these kind of preconditioners has been laid out in [9]. However, we want to emphasize that this independence does not mean, practically speaking, that the given preconditioner is in some sense optimal or even well-performing (and this is even more pronounced if the order-optimality is considered only wrt the spectrum, omitting the conditioning of the eigenbasis of the preconditioned system) – but in particular settings with additional reasoning this might be a useful property as we will demonstrate.

Given a Butcher tableau $A$, the bound (10) can be further approximated using the Joukowsky bound with the Chebyshev polynomials in the complex plane, see [14, Section 6.11 and Corollary 6.33] and also [8, Section 5.7.2] and references therein. A relatively direct and elementary calculation then allows us to evaluate the bound (10), and in Figure 2 we show the GMRES convergence behavior together with these bounds for different Butcher tableaus for $P^d$ on the left (for validation of the above formulas and detailed calculations we refer the interested reader to [11, Section 7.3.1]). Notice that as $h$ is fixed for all methods we obtain different $\tau$ for each of the IRK methods and therefore the scaling factor $\tau/h^2$ becomes $h^{2(1/p-1)}$. Hence, the difference is not only in the choice of $A$ but also in the interval spanned by $\theta_k$ and how close (or far) these are to the limit in (17).

We also want to address the notable stair-like behavior. Since GMRES is invariant (in exact arithmetic) to an orthogonal transformation we can focus on GMRES applied to a problem with
the block-diagonal matrix $X^d := \text{diag}(X_1^d, X_2^d, \dots, X_n^d)$. As noted in [3], the optimal polynomial$^3$ that solves the min-max problem

$$ \frac{\|\mathbf{r}_{\ell}\|}{\|\mathbf{r}_0\|} \le \kappa(S^d) \min_{\substack{\varphi(0)=1 \\ \deg(\varphi) \le \ell}} \max_{1 \le k \le n} \|\varphi(X_k^d)\|, $$

equioscillates over the blocks – but only for even degrees, i.e., for $\ell = 2l$ for some $l \in \mathbb{N}$. Considering our situation, the optimal polynomial is trying to minimize its value over these 2-by-2 blocks, which in our case have eigenvalues $\xi_1^{(k)} \neq \xi_2^{(k)}$. Hence, the optimal polynomial aims to be small (in the maximum norm) on the line segment connecting $1 + i\sqrt{|(a_{12}a_{21})/(a_{11}a_{22})|}$ and $1 - i\sqrt{|(a_{12}a_{21})/(a_{11}a_{22})|}$, symmetrically about 1. However, polynomials with real coefficients of odd degree have at least one real root. This is why we see notably less improvement at odd iterations – the added degree of freedom for the optimization (placement of one additional root) is constrained in a quite inconvenient way for our setting. In contrast to even iterations, on odd iterations we can make the polynomial small only in the vicinity of the real axis, which after a couple of odd iterations decreases the overall maximum norm only very little, e.g., for LobattoIIIC in Figure 4 we see that initially the real root can contribute almost equally to the decrease of the maximum norm of the polynomial, but after seven or nine iterations this is no longer true.

**Remark 4.** This reasoning is similar to the one in [8, Section 5.7.2, p.291], where the GMRES convergence bound (10) is adapted to a symmetric indefinite problem with spectrum in two intervals $I^- \cup I^+ \equiv [-1, \nu] \cup [\nu, 1]$. Then (10) can be simplified to

$$ \frac{\|\mathbf{r}_{\ell}\|}{\|\mathbf{r}_0\|} \le 2 \left( \frac{\kappa - 1}{\kappa + 1} \right)^{[\ell/2]}, \qquad (24) $$

where $\kappa = 1/\nu$ is the condition number of the given system matrix and $[\ell/2]$ is the integer part of $\ell/2$. Notice that this bound gives linear convergence but only over every two iterations of GMRES – precisely the observed behavior in Figure 4.

However adapting this bound to our setting becomes very quickly very complicated and goes beyond the scope of this text. The main difficulty, in our eyes, is in the treatment of the optimal polynomials with real coefficients on two separated line-segments in $\mathbb{C}$. For more details on the real case and the derivation of (24) we refer also to [4, Chapter 3].

## 4.2 Block upper-triangular preconditioner

We consider the preconditioners $P^{\text{GSU}}$, $P^u$ and where necessary we join the quantities corresponding to the preconditioners by writing, e.g., $X_k^{\text{GSU},u}$ instead of $X_k^{\text{GSU}}$ and $X_k^u$. By analogous calculations to Section 4.1, we obtain the formulas

$$ X_k^{\text{GSU}} = \begin{bmatrix} 1 & 0 \\ -\frac{a_{21}\theta_k}{1-a_{11}\theta_k} & \xi_k^{\text{GSU}} \end{bmatrix} \quad \text{and} \quad X_k^u = \begin{bmatrix} 1 & 0 \\ -\frac{a_{21}\theta_k}{1-a_{11}\theta_k} & \xi_k^u \end{bmatrix}, $$

<footnote>$^3$In [3], the authors call these the Chebyshev polynomials of the given matrix – of $X^d$ in this case.</footnote>
hence obtaining the spectrum of the preconditioned systems as the union of $\{1, \xi_k^{\text{GSU,u}}\}_{k=1}^n$ with

$$ \begin{aligned} \xi_k^{\text{GSU}} &= \frac{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22}) + 1}{|\theta|^2 a_{11} a_{22} + |\theta|(a_{11} + a_{22}) + 1}, \\ \xi_k^u &= \frac{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22}) + 1}{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22} - \frac{a_{21}a_{12}}{a_{11}}) + 1}. \end{aligned} \quad (25) $$

Moreover, Lemma 2 shows that for $X_k^{\text{GSU},d} = (S_k^{\text{GSU},d})^{-1}\text{diag}(1, \xi_k^{\text{GSU,u}})S_k^{\text{GSU},d}$ we have

$$ \kappa(S_k^{\text{GSU,u}}) = \sqrt{\frac{\sqrt{1 + (\beta_{\text{GSU,u}})^2} + 1}{\sqrt{1 + (\beta_{\text{GSU,u}})^2} - 1}}, $$

where the scalars $\beta_{\text{GSU,u}}$ are given as functions of $\theta_k$,

$$ |\beta_{\text{GSU}}| = \frac{|a_{12}|}{||\theta_k|^{-1} + a_{22}|}, \quad \text{and} \quad |\beta_u| = \frac{|a_{12}/a_{11}|}{|\theta_k| \frac{\det(A)}{a_{11}} + a_{11}|}. $$

Further calculations give

$$ \begin{aligned} |\xi_k^{\text{GSU}} - 1| &= \left| \frac{|\theta|^2 a_{12} a_{21}}{|\theta|^2 a_{11} a_{22} + |\theta|(a_{11} + a_{22}) + 1} \right|, \\ |\xi_k^u - 1| &= \left| \frac{|\theta| a_{12} a_{21} / a_{11}}{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22} - \frac{a_{21} a_{12}}{a_{11}}) + 1} \right|, \end{aligned} $$

and assuming $a_{11}, a_{22} \ge 0$, elementary calculus reveals that the cluster diameter for $P^{\text{GSU}}$ is bounded from above by $|a_{12}a_{21}|/(a_{11}a_{22})$ – the limit as $|\theta_k| \to +\infty$. For $P^u$ analogous calculations show that the cluster diameter is maximized at $|\theta_k| = \det(A)^{-1/2}$ (assuming $\det(A) \in (|\theta_n|, |\theta_1|)$) with the value

$$ |\xi_k^u - 1| \le \left| \frac{a_{21}a_{12}}{a_{11}^2 + \det(A) + \frac{2a_{11}}{\sqrt{\det(A)}}} \right| $$

The conditioning of the eigenbasis $S_k^{\text{GSU,u}}$ of $X_k^{\text{GSU,u}}$ can be treated similarly, first observing that $\kappa(S_k^{\text{GSU,u}})$ is a decreasing function of $|\beta_{\text{GSU,u}}| \in (0, +\infty)$ with a singularity at 0. Hence, as $h, \tau \to 0$, the matrix $X_1^{\text{GSU}}$ (and $X_n^u$) becomes non-diagonalizable$^4$ as $|\theta_1| \to 0$ (and $|\theta_n| \to +\infty$). We show the preconditioned GMRES convergence behavior and the evaluated bounds in Figure 2; for detailed calculations and validation of the above formulas we refer the interested reader to [11, Section 7.3.2]. Also, notice that the stair-like behavior of $P^d$ is not present as the entire spectrum is real and covers reasonably uniformly an interval on a real line and a single point 1. As the authors point out in [8, Section 5.6.2 and 5.7.2], as long as the condition number $\kappa(S^{\text{GSU,u}})$ is not too large, the classical linear bound based on the condition number of the preconditioned system matrix (i.e., the condition number of $M(P^{\text{GSU,u}})^{-1}$) can be quite descriptive for the worst-case GMRES behavior. If the preconditioned system is not far away from being normal, then this explains why we only

<footnote>$^4$We can see this already in (25) since in the limits considered we get $\xi_k^{\text{GSU,u}} = 1$ thus obtaining a 2-by-2 Jordan block.</footnote>
see the linear convergence – the spectrum is populating the interval considered densely enough so that the superlinear convergence argument used in exact arithmetic is not applicable (i.e., we cannot single out any outliers if there are none, see [8, Section 5.6.4 (Figure 5.7 in particular) and also Section 5.6.4]). Hence, linear convergence is, in principle, to be expected. Also, this explains the observation that the number of GMRES iterations does not grow as $h \to 0$ as observed in [13] (an analogous argument applies also to the block-diagonal preconditioner). Notice that the accuracy of the bound (24) is usually supported by the same kind of arguments. By analogy, this also qualitatively explains the linear in the linear-over-two-iterations convergence behavior of the block-diagonal preconditioned system.

## 4.3 Block lower-triangular preconditioner

The results for $P^{GSL}, P^1$ are completely analogous to the ones from Section 4.2 and hence we just present these without much comment; more details can be found in [11, Section 7.3.3]. We have

$$X_k^{GSL} = \begin{bmatrix} 1 & -\frac{a_{12}\theta_k}{1-a_{11}\theta_k} \\ 0 & \xi_k^{GSL} \end{bmatrix} \quad \text{and} \quad X_k^1 = \begin{bmatrix} 1 & -\frac{a_{12}\theta_k}{1-a_{11}\theta_k} \\ 0 & \xi_k^1 \end{bmatrix},$$

hence obtaining the spectrum of the preconditioned systems as the union of $\{1, \xi_k^{GSL,1}\}_{k=1}^n$ with

$$\begin{align}
\xi_k^{GSL} &= \xi_k^{GSU} = \frac{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22}) + 1}{|\theta|^2 a_{11} a_{22} + |\theta|(a_{11} + a_{22}) + 1}, \tag{26} \\
\xi_k^1 &= \xi_k^u = \frac{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22}) + 1}{|\theta|^2 \det(A) + |\theta|(a_{11} + a_{22} - \frac{a_{21}a_{12}}{a_{11}}) + 1}. \notag
\end{align}$$

Moreover, Lemma 2 shows that for $X_k^{GSL,l} = (S_k^{GSL,l})^{-1}\text{diag}(1, \xi_k^{GSL,l})S_k^{GSL,l}$ we have

$$\kappa(S_k^{GSL,l}) = \sqrt{\frac{\sqrt{1 + (\beta_{GSL,l})^2 + 1}}{\sqrt{1 + (\beta_{GSL,l})^2} - 1}},$$

where the scalars $\beta_{GSL,l}$ are given as functions of $\theta_k$,

$$|\beta_{GSL}| = \frac{|a_{21}|}{||\theta_k|^{-1} + a_{22}|}, \quad \text{and} \quad |\beta_l| = \frac{|a_{21}/a_{11}|}{||\theta_k| \frac{\det(A)}{a_{11}} + a_{11}|}.$$

Notice that the results are either identical or very similar to the ones obtained with $P^{GSU}, P^u$ and a comparison of the convergence behavior as well as of the bounds in Figure 2 and Figure 3 reflects this.

## 4.4 Optimized Butcher tableaus

In [13], the authors consider taking a different $A$, let us denote it by $\tilde{A}$, for the construction of the preconditioner, motivated by [17]. $\tilde{A}$ is chosen as a result of an optimization routine seeking to minimize $\kappa(\tilde{A}^{-1}A)$ (or $\kappa(A\tilde{A}^{-1})$ for a right preconditioner) subject to a particular non-zero pattern and $\text{diag}(\tilde{A}) = \text{diag}(A)$. This optimization thus evaluates only quantities corresponding to the s-bys Butcher tableaus, which should be negligible in cost compared to the solution process of the
<figure><img src="image_5.png" /><figcaption>Figure 3: We show the GMRES performance together with the evaluated bound (10) for different preconditioners and Butcher tableaus. We set $N = 100$, $\tau = h^{2/p}$ ($p$ is the order of convergence of the IRK method; $p_{Gauss} = 4$, $p_{RadauIIA} = 3$, $p_{LobattoIIIC} = 2$) and use a random right-hand side.</figcaption></figure>

stage equations. Numerical results in [13, Table 4.3] show that the resulting preconditioner lowers the condition number of the preconditioned system. Having explicit formulas from Section 4.1–4.3 we can move from minimizing the quantity $\kappa(\tilde{A}^{-1}A)$ (or $\kappa(A\tilde{A}^{-1})$) to minimizing quantities in the bound (10) so that a provable bound can be obtained based on the result.

First, adapting the eigenpair formulas in Section 4.1 – 4.3 for the case $A \neq \tilde{A}$ we could optimize the eigenproperties of the preconditioned system and thus obtain a theoretical bound on the resulting preconditioner. Using Lemma 1, the matrices $X_k$ are easily accessible with any $A$ fixed.

Second, having $A = \tilde{A}$ as in Section 4.1 – 4.3, we can still consider the same optimization problem – but naturally with extra constraints so that the resulting IRK method has some desired order of convergence and stability properties.

We numerically test these approaches, defining the objective function to be minimized as a weighted sum of the cluster diameter<fnref n="5" /> $|\xi_k - 1|$ and the conditioning of the eigenbasis of the preconditioned system matrix $\kappa(S)$,

$$f_{\text{obj}}(A, \tilde{A}) := \max_k |\xi_k - 1| + \omega\kappa(S),$$

with some positive small weight $\omega > 0$. Minimizing $f_{\text{obj}}$ then aims to minimize the min-max part of (10) while keeping the eigenbasis conditioning under control. In case that the clustering is very tight we expect to get good GMRES convergence. Before showing the results we comment on the usefulness of the optimization approach for the case $A = \tilde{A}$.

**Remark 5.** In order to optimize over $A = \tilde{A}$ we have to keep in mind that we are at the same time changing the IRK method we are using to solve the system of ODEs. This introduces some

<footnote n="5">⁵The clustering point is for all of our preconditioners naturally equal to 1.</footnote>
optimization constraints to ensure good order of convergence (equality constraints on A, b, c) and stability (equality and inequality constraints on A, b, c). In our experience, in order to obtain a notable improvement in minimization of $f_{obj}$ over choosing $A = A_{Gauss, RadauIIA,}$..., we need to relax some of the qualities of the IRK methods, e.g., decrease the order of convergence by one or give up some stability. This seems quite unappealing as the maximal order of convergence $p_{Gauss} = 4$ is already not too large for practical use.

In our eyes, this makes the optimization for $A = \tilde{A}$ better suited for larger $s$ (e.g., $s \ge 5$), where decreasing the order of convergence by one or two might not be a significant drawback as the order of convergence of the standard IRK methods is much higher; see [11, Section 7.5].

Note that even though Remark 5 suggests that the optimization with $A = \tilde{A}$ is not really of interest here, it is still reasonable to start studying it for $s = 2$, to get some insight in the case where the formulas are explicitly available. We show the numerical results in Figure 4 and the resulting matrices in Table 1.

We see that for the setting $A = \tilde{A}$ in the first row of Figure 4 the optimized $A$ are such that we converge after two to five iterations, basically turning GMRES into a direct solver$^6$. This seems natural looking at Table 1 – we observe that the Butcher tableau matrix $A$ adapts the non-zero structure of the preconditioner, e.g., $A$ becomes close to diagonal for $P^d$, so that the preconditioner then becomes almost identical to $M$, reinforcing the point of Remark 5. Notably, this is achieved while not exploding the condition number of the eigenbasis of the preconditioned system.

For the second row, i.e., $A \neq \tilde{A}$ and optimizing $\tilde{A}$ that is used to construct the preconditioner, we see that for the block-diagonal preconditioner we can still obtain a considerable speed-up, as well as for the preconditioners $P^{GSL,GSU}$. This is not the case, however, for the preconditioners $P^{l,u}$. The reason is that with $A = \tilde{A}$, the spectrum of $M(P^{l,u})^{-1}$ was real but this property is lost in the general case $A \neq \tilde{A}$. Hence, even though we have tightened the clustering of the eigenvalues by optimizing $\tilde{A}$ this wasn't significant enough to off-set the introduction of the complex eigenvalues of the preconditioned system.

We would like to note that we obtained similar results to the first row when considering $p = 3$ and $A$-stability (i.e., giving up the more restrictive property of $L$-stability) but for $p = 3$ and $L$-stable methods there seemed to be next to no gains from the extra optimization (see [11, Section 7.5]).

<table><caption>Table 1: The resulting matrices $A = \tilde{A}$ (first row) or $\tilde{A}$ (second row) corresponding to the results presented in Figure 4.</caption><thead><tr><th></th><th>block diagonal</th><th colspan="2">block lower-triangular</th><th colspan="2">block upper-triangular</th></tr><tr><th></th><th><i>P</i><sup>d</sup></th><th><i>P</i><sup>GSL</sup></th><th><i>P</i><sup>l</sup></th><th><i>P</i><sup>GSU</sup></th><th><i>P</i><sup>u</sup></th></tr></thead><tbody><tr><th><i>A</i> = <i>A</i>&#773;</th><td><table><tbody><tr><td>1.8 &middot; 10<sup>-1</sup></td><td>-3.8 &middot; 10<sup>-6</sup></td></tr><tr><td>2.3 &middot; 10<sup>-6</sup></td><td>3.9 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>3.5 &middot; 10<sup>-1</sup></td><td>-4.5 &middot; 10<sup>-9</sup></td></tr><tr><td>2.5 &middot; 10<sup>-6</sup></td><td>2.3 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>3.2 &middot; 10<sup>-1</sup></td><td>-3.0 &middot; 10<sup>-8</sup></td></tr><tr><td>9.6 &middot; 10<sup>-1</sup></td><td>2.7 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>3.5 &middot; 10<sup>-1</sup></td><td>4.2 &middot; 10<sup>-1</sup></td></tr><tr><td>-5.8 &middot; 10<sup>-9</sup></td><td>2.3 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>4.9 &middot; 10<sup>-1</sup></td><td>-5.2 &middot; 10<sup>-2</sup></td></tr><tr><td>6.3 &middot; 10<sup>-8</sup></td><td>1.7 &middot; 10<sup>-2</sup></td></tr></tbody></table></th></tr><tr><th><i>A</i> &ne; <i>A</i>&#773;</th><td><table><tbody><tr><td>9.3 &middot; 10<sup>-1</sup></td><td>4.8 &middot; 10<sup>-2</sup></td></tr><tr><td>4.0 &middot; 10<sup>-1</sup></td><td>1.3 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>5.1 &middot; 10<sup>-1</sup></td><td>-6.6 &middot; 10<sup>-2</sup></td></tr><tr><td>8.0 &middot; 10<sup>-1</sup></td><td>3.0 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>4.1 &middot; 10<sup>-1</sup></td><td>-8.4 &middot; 10<sup>-2</sup></td></tr><tr><td>7.4 &middot; 10<sup>-1</sup></td><td>2.2 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>4.1 &middot; 10<sup>-1</sup></td><td>-8.2 &middot; 10<sup>-2</sup></td></tr><tr><td>7.0 &middot; 10<sup>-1</sup></td><td>3.7 &middot; 10<sup>-1</sup></td></tr></tbody></table></th><td><table><tbody><tr><td>4.1 &middot; 10<sup>-1</sup></td><td>-8.1 &middot; 10<sup>-2</sup></td></tr><tr><td>7.2 &middot; 10<sup>-1</sup></td><td>2.3 &middot; 10<sup>-1</sup></td></tr></tbody></table></th></tr></tbody></table>

Last but not least we comment on the computational costs. In order to evaluate the complete eigenproperties of the preconditioned system $M(P^*)^{-1}$ (or $(P^*)^{-1}M$) we first calculate the eigen-decomposition of $L$ (which is prohibitively costly) and then we calculate eigendecompositions of $n$

<footnote>$^6$Notice that in practice the solve accuracy also needs to be balanced with the discretization errors, hence even five iterations to obtain a relative residual smaller than machine precision might be considered a direct solver (i.e., converging after one or two steps) in some applications.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4: GMRES convergence (and bounds) for block preconditioners $P^d$ (left column), $P^{GSL}$, $P^l$ (middle column) and $P^{GSU}$, $P^u$ (right column) for the two stage RadauIIA Butcher tableau (filled markers, same results as in Figure 2 and 3) and also the optimized setting considered above (half-filled markers). In the second row the optimization process changes only the matrix $\tilde{A}$ used to construct the preconditioner (and not $M$) in contrast to the first row where the optimized is used for construction of both the system matrix and the preconditioner (in this case we imposed $A$ to be such that the resulting IRK method is at least of order two and is $L$-stable, the constraints can be found in [11, Section 7.3.4]). We used $\omega = 10^{-5}$, $N = 100$ and $\tau = h^{2/p}$ as above.</figcaption></figure>

matrices $X_k^*$, each of dimension $s$, which can be done in parallel. Even without any parallelization the cost is linear in $n$ with the constant corresponding to $s^3$. In practice, we often have (or can reasonably cheaply obtain) some estimate of the eigeninformation of $L$, e.g., estimates $\mu_1, \mu_n$ of the true eigenvalues $\lambda_1, \lambda_n$ or on the conditioning of $V$. The eigendecomposition of $L$ can then be replaced by, e.g., considering only $q$ “fake” $\lambda_k$ sampled from $(\mu_{\min}, \mu_{\max})$ (assuming we know the spectrum of $L$ is real).

We show a numerical illustration, considering the same setting as for Figure 4 and sample only $q$ distinct $\lambda_k$ (and thus $\theta_k$) in the interval $\mu_{\min} := \theta_{\min}, \mu_{\max} := \theta_{\max}$. We show the evolution of the number of preconditioned GMRES with the optimized preconditioners depending on $q$ in Figure 5. We see that for the spatial dimension $n = 99^2$ we already get comparable results to Figure 4 by sampling *only very few* “fake” $\theta_k$ compared to optimizing over the actual roughly ten thousands of them. Also, the bounds stay still quite descriptive – only for $q = 40$ and $P^{GSU}$ have we found that running the optimization as a black-box does not give a useful descriptive bound$^7$. This observation is key for either of the optimization approaches to be viable and further study of $q$ for more involved settings (and its dependency on $n$) and/or sampling strategies for complex spectra is necessary. However, as the main interest is in $s$ large, this will be presented in the follow-up work.

<footnote>$^7$In fact for that particular setting the optimization routine sets off in a direction of *inadmissible* $A$ (singular) and does not find its way back. Fine tuning the parameters of the optimization routine does, however, fix this issue.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5</figcaption></figure>

# 5 Generalizations and conclusion

We have shown that for two stage IRK methods the preconditioners from [13] can be analyzed explicitly using spectral techniques and the GMRES convergence behavior can be reasonably predicted using the worst-case GMRES bound. This bound can be evaluated directly before the computations and gives a theoretical background to the results observed in [13]. These results are given for a simple test problem but the analysis clearly extends to any diagonalizable spatial operator $L$ – not just the Laplacian. The same is true also for the discretization scheme used – a finite elements scheme analysis can be done completely analogously (see [11, Section 7.7]). Some of the above can be generalized to $s$-stage IRK methods but due to the space restrictions these results will be presented in an upcoming manuscript.

We would also like to mention that both in [10] and [16], the authors multiplied $M$ with $A^{-1} \otimes I_n$ from the left. The analysis above for this case is analogous and reveals that for some preconditioners we obtain much better properties and for others the performance deteriorates. This is naturally very dependent on the choice of $A$ as well; for more details see [11, Section 7.6].

# References

[1] M. Arioli, V. Pták, and Z. Strakoš. Krylov sequences of maximal length and convergence of GMRES. *BIT Numerical Mathematics*, 38(4):636–643, 1998.

[2] M. R. Clines, V. E. Howle, and K. R. Long. Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs. arXiv: https://arxiv.org/abs/2206.08991, 2022.
[3] V. Faber, J. Liesen, and P. Tichý. On Chebyshev polynomials of matrices. *SIAM Journal on Matrix Analysis and Applications*, 31(4):2205–2221, 2010.

[4] B. Fischer. *Polynomial Based Iteration Methods for Symmetric Linear Systems*. Springer, Wiesbaden, 1996.

[5] A. Greenbaum and Z. Strakoš. *Matrices that generate the same Krylov residual spaces*. Springer, 1994.

[6] A. Greenbaum, V. Pták, and Z. Strakoš. Any nonincreasing convergence curve is possible for GMRES. *SIAM Journal on Matrix Analysis and Applications*, 17(3):465–469, 1996.

[7] J. Liesen and Z. Strakoš. GMRES convergence analysis for a convection-diffusion model problem. *SIAM Journal on Scientific Computing*, 26(6):1989–2009, 2005.

[8] J. Liesen and Z. Strakoš. *Krylov Subspace Methods: Principles and Analysis*. Oxford University Press, Oxford, 2013.

[9] K.-A. Mardal, T. K. Nilssen, and G. A. Staff. Order-optimal preconditioners for implicit Runge-Kutta schemes applied to parabolic PDEs. *SIAM Journal on Scientific Computing*, 29 (1):361–375, 2007.

[10] M. Neytcheva and O. Axelsson. Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order. In P. Frolkovič, K. Mikula, and D. Ševčovič, editors, *Proceedings of the Conference Algoritmy 2020*. Slovak University of Technology in Bratislava, Vydavateľstvo SPEKTRUM, 2020.

[11] M. Outrata. *Schwarz methods, Schur complements, preconditioning and numerical linear algebra*. PhD thesis, University of Geneva, Math Department, 2022.

[12] D. Palitta and V. Simoncini. Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations. *Vietnam Journal of Mathematics*, 48(4):791–807, 2020.

[13] M. M. Rana, V. E. Howle, K. Long, A. Meek, and W. Milestone. A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems. *SIAM Journal on Scientific Computing*, 43(5):S475–S495, 2021.

[14] Y. Saad. *Iterative Methods for Sparse Linear Systems*. Other Titles in Applied Mathematics. SIAM, Philadelphia, Second edition, 2003. ISBN 978-0-89871-534-7.

[15] B. S. Southworth, O. Krzysik, and W. Pazner. Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs. *SIAM Journal on Scientific Computing*, 44(2):636–663, 2022.

[16] B. S. Southworth, O. Krzysik, W. Pazner, and H. De Sterck. Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting. *SIAM Journal on Scientific Computing*, 44(1):416–443, 2022.

[17] G. A. Staff, K.-A. Mardal, and T. K. Nilssen. Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs. *Modeling, Identification and Control*, 27(2):109–123, 2006.
[18] C. F. Van Loan. The ubiquitous Kronecker product. *Journal of Computational and Applied Mathematics*, 123(1-2):85-100, 2000.

[19] G. Wanner and E. Hairer. *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer Berlin, Heidelberg, 1996.

[20] G. Wanner, S. P. Nørsett, and E. Hairer. *Solving Ordinary Differential Equations I: Non-Stiff Problems*. Springer Berlin, Heidelberg, 1987.