# A new ParaDiag time-parallel time integration method

Martin J. Gander Section de Mathèmatiques, University of Geneva, Switzerland
<EMAIL_ADDRESS>Davide Palitta Dipartimento di Matematica and AM2,
Alma Mater Studiorum - Università di Bologna, Piazza di Porta S. Donato, 5,
I-40127 Bologna, Italy<EMAIL_ADDRESS>
###### Abstract

Time-parallel time integration has received a lot of attention in the high
performance computing community over the past two decades. Indeed, it has been
shown that parallel-in-time techniques have the potential to remedy one of the
main computational drawbacks of parallel-in-space solvers. In particular, it
is well-known that for large-scale evolution problems space parallelization
saturates long before all processing cores are effectively used on today’s
large scale parallel computers. Among the many approaches for time-parallel
time integration, ParaDiag schemes have proven to be a very effective
approach. In this framework, the time stepping matrix or an approximation
thereof is diagonalized by Fourier techniques, so that computations taking
place at different time steps can be indeed carried out in parallel. We
propose here a new ParaDiag algorithm combining the Sherman-Morrison-Woodbury
formula and Krylov techniques. A panel of diverse numerical examples
illustrates the potential of our new solver. In particular, we show that it
performs very well compared to different ParaDiag algorithms recently proposed
in the literature.

##  1 Introduction

Time-parallel time integration is currently a very active field of research
within the high performance computing community. Research interest was
relaunched over twenty years ago with the introduction of the Parareal
algorithm [22], which is a two-level, non-intrusive method that allows
existing codes to be parallelized, and works well on parabolic problems; see,
e.g., [14, 8] for detailed convergence analyses. Nowadays, many other
different techniques for time parallelization of evolution problems can be
found on the market: methods based on multiple shooting (like Parareal),
methods based on Domain Decomposition and Waveform relaxation, Multigrid type
methods, and even direct time-parallel solvers; see, e.g., [7] for a thorough
review of such schemes. The boundaries between these different solvers have
become less and less strict over the years, with tools designed for a specific
method being fully exploited in others. The ParaDiag family of time-parallel
time integrators is a typical example of such permeability. Originally,
ParaDiag methods were designed as direct time-parallel solvers [24]. To
overcome the nondiagonalizability of the time stepping matrix which is a
Jordan block for a constant time step, it was first proposed to use different
time steps. However, this trick can lead to very ill-conditioned eigenvector
matrices which can potentially pollute the entire solution process, especially
for fine time grids. A thoughtful tradeoff between having similar time steps
for accuracy and different ones for diagonalizability is thus crucial to make
this first ParaDiag method successful; see, e.g., [11, 10] for a detailed
analysis. Due to the practical limitation this tradeoff imposes on the number
of time steps that can be parallelized, a recent approach consists in
considering an approximate problem where the time stepping matrix is
periodic111A anonymous reviewer pointed out that this corresponds to
approximating a partial fraction expansion of the rational approximation for
the matrix exponential defined by the numerical scheme (which is equivalent to
a Weierstrass normal form [5], here essentially the Jordan Canonical Form that
was discovered independently [18]), by imposing periodicity to obtain a normal
form which is diagonal.. Thanks to this approximation, the ParaDiag approach
can be well combined with other time-parallel methods as, e.g., Parareal [31,
16] or MGRIT [33]. See also [15] for time-periodic waveform relaxation for
initial value problems and [12] for Parareal algorithms for truly time-
periodic problems.

A different approach consists in employing ParaDiag techniques within
iterative methods. For instance,  one approach approximates the time stepping
matrix by a circulant matrix, with the latter being used to define a
preconditioner for Krylov methods. A ParaDiag scheme can then be used for a
more efficient application of the preconditioning operator; see, e.g., the
iterative time parallelization in [25] for parabolic problems, and [2] for
hyperbolic problems. Also for non-linear problems, ParaDiag algorithms
necessarily become iterative; see, e.g., [9].

ParaDiag techniques have been developed for optimal control problems [32] as
well. In this setting, very good performance is obtained by using
$\alpha$-circulant modifications of the time stepping matrix [23, 34]. See
also [13] for an overview and implementation details about these schemes.

In contrast to the first attempts of time parallelization by diagonalization
using different time steps, these modern ParaDiag methods are very successful
in solving both parabolic and hyperbolic evolution problems, and new ideas in
this direction are currently being developed; see, e.g., [19] where
interpolation and low-rank techniques are proposed and studied.

We present here a new ParaDiag algorithm for solving in a time-parallel
fashion the evolution problem

(1) $\begin{array}[]{rll}u_{t}&=&\mathfrak{L}(u)+f,\quad\text{in
}\Omega\times(0,T],\\\ u&=&g,\quad\text{on }\partial\Omega,\\\
u(0)&=&u_{0},\\\ \end{array}$

where the spatial domain $\Omega$ is such that $\Omega\subset\mathbb{R}^{d}$,
$d=1,2,3$, and $\mathfrak{L}$ is a linear differential operator involving only
spatial derivatives. If we discretize (1) in space with a finite element or
finite difference method with $\bar{n}$ degrees of freedom, and use a backward
Euler scheme with $\ell$ time steps, the _all-at-once_ discretization of (1)
can be written in matrix form as

(2) $(I+\tau
K)U-U\Sigma_{1}^{T}=[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}],$

where $K\in\mathbb{R}^{\bar{n}\times\bar{n}}$ is the stiffness matrix stemming
from the spatial discretization, $\Sigma_{1}\in\mathbb{R}^{\ell\times\ell}$ is
a zero matrix having ones only on the first subdiagonal, the $j$-th column of
$U=[\mathbf{u}_{1},\ldots,\mathbf{u}_{\ell}]\in\mathbb{R}^{\bar{n}\times\ell}$
represents the approximation to the solution $u$ at time $t_{j}$,
$j=1,\ldots,\ell$, $\tau=T/\ell$ is the time step, and $\mathbf{u_{0}}$ and
$\mathbf{f}_{j}$ gather the nodal values of $u_{0}$ and $f(t_{j})$ along with
the boundary conditions. See [27].

Recently, the matrix equation formulation (2) has been used to design new
solution procedures. In particular, low-rank solvers can be very successful in
solving (2) whenever the right-hand side
$[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]$ has low
rank; see, e.g., [27]. On the other hand, the performance of such methods
significantly worsens for right-hand sides with a sizable rank. In this paper,
we  address the performance issue when the right-hand side is possibly full
rank. In particular, we propose a new ParaDiag algorithm which is able to
fully take advantage of the circulant-plus-low-rank structure of $\Sigma_{1}$
so as to design an efficient parallel-in-time algorithm for solving (2). A
significant advantage of this approach over low-rank space-time schemes is
that no assumption on the (numerical) rank of the right hand side
$[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]$ is needed
in (2).

Here is an outline of the paper. In section 2 we show how to use the
circulant-plus-low-rank structure of $\Sigma_{1}$ by combining (2) with the
matrix-oriented Sherman-Morrison-Woodbury formula. Such an approach sees as an
intermediate step the solution of a linear system whose coefficient matrix has
a rather involved structure. An ad-hoc projection technique for the solution
of this _inner_ linear system is proposed in section 3 where we also study
some of the properties of the coefficient matrix in case of symmetric positive
definite (SPD) stiffness matrices $K$. In section 4 we propose a very
successful variant of the ParaDiag scheme which makes use of
$\alpha$-circulant matrices. In section 5 we generalize our approach to the
case of higher-order time discretization schemes. The potential of our new
approach is illustrated in section 6 where several numerical results are
shown. In section 7 we draw our conclusions.

Throughout the paper we adopt the following notation. Capital letters $(A)$
denote matrices, bold, lower-case letters ($\mathbf{a}$) vectors, and plain,
lower-case letters ($a$) scalars.
$I_{n}=[\mathbf{e}_{1},\ldots,\mathbf{e}_{n}]$ denotes the identity matrix of
order $n$ and the subscript is omitted whenever the dimensions of $I$ are
clear from the context; $\otimes$ is the Kronecker product whereas
$\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the minimum and the maximum
eigenvalue of $A$, respectively. Given a matrix
$X=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]\in\mathbb{R}^{m\times n}$,
$\text{vec}(X)\in\mathbb{R}^{mn}$ denotes the vector obtained by stacking the
columns of $X$ on top of each other, namely
$\text{vec}(X)=[\mathbf{x}_{1}^{T},\ldots,\mathbf{x}_{n}^{T}]^{T}$.

##  2 The new ParaDiag algorithm

Even though the discrete backward Euler operator $\Sigma_{1}$ cannot be
diagonalized as it is a Jordan block, its circulant-plus-rank-one structure
can be exploited to design efficient solvers for (2). Indeed, we can write

(3) $\Sigma_{1}=\begin{bmatrix}0&&&\\\ 1&\ddots&&\\\ &\ddots&\ddots&\\\
&&1&0\\\ \end{bmatrix}=C_{1}-\mathbf{e}_{1}\mathbf{e}_{\ell}^{T},\quad
C_{1}=\begin{bmatrix}0&&&1\\\ 1&\ddots&&\\\ &\ddots&\ddots&\\\ &&1&0\\\
\end{bmatrix},$

where $C_{1}$ is a circulant matrix.

The relation above has been used in, e.g., [25] to derive preconditioning
operators for the linear system counterpart of (2). In particular, such
preconditioners were obtained by dropping the rank-1 term
$\mathbf{e}_{1}\mathbf{e}_{\ell}^{T}$ in (3).

In our setting, by inserting (3) into (2), we get

(4) $(I+\tau
K)U-UC_{1}^{T}+U\mathbf{e}_{\ell}\mathbf{e}_{1}^{T}=[\mathbf{u_{0}}+\mathbf{f}_{1},\ldots,\mathbf{f}_{\ell}].$

Since $C_{1}$ is a circulant matrix, it can be diagonalized by the Fast
Fourier Transform (FFT), namely
$\Pi_{1}={F}C_{1}{F}^{-1}=\text{diag}({F}C_{1}\mathbf{e}_{1})=\text{diag}(\pi_{1},\ldots,\pi_{\ell})$,
where ${F}$ denotes the discrete Fourier matrix.

By postmultiplying (4) by ${F}^{T}$ we get

(5) $(I_{\bar{n}}+\tau
K)\widetilde{U}-\widetilde{U}\Pi_{1}+\widetilde{U}{F}^{-T}\mathbf{e}_{\ell}\mathbf{e}_{1}^{T}{F}^{T}=[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T},\quad\widetilde{U}=U{F}^{T},$

which can also be written in Kronecker form as

$\left(I_{\ell}\otimes(I_{\bar{n}}+\tau K)-\Pi_{1}\otimes
I_{\bar{n}}+{F}\mathbf{e}_{1}\mathbf{e}_{\ell}^{T}{F}^{-1}\otimes
I_{\bar{n}}\right)\text{vec}(\widetilde{U})=\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T}).$

If we look at the system matrix above, this can be viewed as the sum of two
components. A main part $P:=I_{\ell}\otimes(I_{\bar{n}}+\tau K)-\Pi_{1}\otimes
I_{\bar{n}}$, and a low-rank modification term $MN^{T}$ where
$M:={F}\mathbf{e}_{1}\otimes I_{\bar{n}}$ and
$N:={F}^{-T}\mathbf{e}_{\ell}\otimes I_{\bar{n}}$. The Sherman-Morrison-
Woodbury formula thus implies that

$\displaystyle\text{vec}(\widetilde{U})=$ $\displaystyle
P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T})$
(6)
$\displaystyle-P^{-1}M(I+N^{T}P^{-1}M)^{-1}N^{T}P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T}).$

Since $P$ is block diagonal,

$P=\begin{bmatrix}(1-\pi_{1})I_{\bar{n}}+\tau K&&\\\ &\ddots&\\\
&&(1-\pi_{\ell})I_{\bar{n}}+\tau K\\\
\end{bmatrix}\in\mathbb{R}^{\bar{n}\ell\times\bar{n}\ell},$

the action of its inverse can be efficiently computed in parallel. In
particular, we can write

$P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T})=\text{vec}(L),$

where

$\displaystyle L:=[((1-\pi_{1})I_{\bar{n}}+\tau
K)^{-1}[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T}\mathbf{e}_{1},\ldots,$
$\displaystyle\hskip 100.00015pt((1-\pi_{\ell})I_{\bar{n}}+\tau
K)^{-1}[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T}\mathbf{e}_{\ell}].$

By using the property of the Kronecker product, the last part in the second
term on the right of (2) is such that

$N^{T}P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T})=N^{T}\text{vec}(L)=L{F}^{-T}\mathbf{e}_{\ell},$

and the linear system with $I+N^{T}P^{-1}M$ we need to solve in (2) can thus
be written as

(7)
$(I+N^{T}P^{-1}M)\mathbf{x}=\mathbf{b},\qquad\mathbf{b}:=L{F}^{-T}\mathbf{e}_{\ell}.$

A naive approach for solving (7) would be to first construct the coefficient
matrix, and then apply one’s favorite linear system solver. However, we notice
that we cannot explicitly compute the coefficient matrix $I+N^{T}P^{-1}M$ as
this  requires applying $P^{-1}$ to all the $\bar{n}$ columns of the matrix
$M$, destroying potential parallel-in-time acceleration.

We propose to use an ad-hoc projection scheme to solve (7). To this end, we
start by deriving an explicit expression for $N^{T}P^{-1}M$. We first notice
that ${F}\mathbf{e}_{1}\in\mathbb{R}^{\ell}$ is the vector of all ones, a
property of the discrete Fourier transform. Therefore,
$M=[I_{\bar{n}},\ldots,I_{\bar{n}}]^{T}$. Moreover, if

${F}^{-T}\mathbf{e}_{\ell}=\begin{bmatrix}\gamma_{1}\\\ \vdots\\\
\gamma_{\ell}\\\ \end{bmatrix}\quad\text{ then }\quad
N=\begin{bmatrix}\gamma_{1}I_{\bar{n}}\\\ \vdots\\\
\gamma_{\ell}I_{\bar{n}}\\\ \end{bmatrix},$

which implies that the term $N^{T}P^{-1}M$ we want to express explicitly can
be written as

$\displaystyle N^{T}P^{-1}M=$
$\displaystyle\begin{bmatrix}\gamma_{1}I_{\bar{n}}&\cdots&\gamma_{\ell}I_{\bar{n}}\\\
\end{bmatrix}\begin{bmatrix}((1-\pi_{1})I_{\bar{n}}+\tau K)^{-1}&&\\\
&\ddots&\\\ &&((1-\pi_{\ell})I_{\bar{n}}+\tau K)^{-1}\\\
\end{bmatrix}\begin{bmatrix}I_{\bar{n}}\\\ \vdots\\\ I_{\bar{n}}\\\
\end{bmatrix}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}.$

Therefore, the linear system (7) can be rewritten in the form

(8)
$\underbrace{\left(I+\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}\right)}_{=:J_{\ell}}\mathbf{x}=\mathbf{b}.$

We propose to apply a projection method for solving problem (8) where the
residual vector is imposed to be orthogonal to a suitable subspace. This is a
very general machinery and its effectiveness depends mainly on the
approximation space one uses. The details of our approach are given in section
3, and for the moment we simply denote by $\mathbf{x}_{m}$ the computed
approximation222The index $m$ denotes the number of iterations performed by
the projection method in Algorithm 3.1 to obtain the sought approximation; see
section 3. to $\mathbf{x}$, solution to (8). Going back to (2), we use
$\mathbf{x}_{m}$ to compute

$W=P^{-1}M\mathbf{x}_{m}=P^{-1}\text{vec}(\mathbf{x}_{m}\mathbf{e}_{1}^{T}{F}^{T})=P^{-1}\text{vec}(\mathbf{x}_{m}\mathbf{1}_{\ell}^{T}),$

where $\mathbf{1}_{\ell}\in\mathbb{R}^{\ell}$ denotes the vector of all ones,
and $P^{-1}$ can be applied column-wise in parallel once again, since it is
block diagonal. We thus get

$W=[((1-\pi_{1})I_{\bar{n}}+\tau
K)^{-1}\mathbf{x}_{m},\ldots,((1-\pi_{\ell})I_{\bar{n}}+\tau
K)^{-1}\mathbf{x}_{m}].$

To conclude, the solution $U$ is then obtained by computing

$U=(L-W){F}^{-T}.$

The pseudocode of our new ParaDiag method is given in Algorithm 2.1.

input : $K\in\mathbb{R}^{\bar{n}\times\bar{n}}$, $\mathbf{u}_{0}$,
$\mathbf{f}_{i}\in\mathbb{R}^{\bar{n}}$, $i=1,\ldots,\ell$,
$\ell\in\mathbb{N}$.

output : $U\in\mathbb{R}^{\bar{n}\times\ell}$ approximate solution to (2).

1 Compute $[\pi_{1},\ldots,\pi_{\ell}]^{T}={F}C_{1}\mathbf{e}_{1}$, and
$[\gamma_{1},\ldots,\gamma_{\ell}]^{T}={F}^{-T}\mathbf{e}_{\ell}$

2

3 parfor $i=1,\ldots,\ell$

4 Set $L\mathbf{e}_{i}=((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}[\mathbf{u}_{0}+\tau\mathbf{f}_{1},\tau\mathbf{f}_{2},\ldots,\tau\mathbf{f}_{\ell}]{F}^{T}\mathbf{e}_{i}$

5

6

7Compute $\mathbf{b}=L{F}^{-T}\mathbf{e}_{\ell}$

8

9Compute $\mathbf{x}_{m}$ by applying Algorithm 3.1 to (8)

10

11 parfor $i=1,\ldots,\ell$

12 Set $W\mathbf{e}_{i}=((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}\mathbf{x}_{m}$

13

14

Set $U=(L-W){F}^{-T}$

Algorithm 2.1 New ParaDiag Algorithm with Backward Euler

##  3 The projection method for the inner linear system

In this section we describe the solution of (8) by a suitable projection
method. In particular, we compute a numerical approximation
$\mathbf{x}_{m}=V_{m}\mathbf{y}_{m}\approx\mathbf{x}$, where the orthonormal
columns of $V_{m}\in\mathbb{R}^{\bar{n}\times m}$ span a suitable subspace
$\mathcal{K}_{m}$, namely $\mathcal{K}_{m}=\text{Range}(V_{m})$, and the
vector $\mathbf{y}_{m}\in\mathbb{R}^{m}$ is computed by imposing an
orthogonality condition on the residual
$\mathbf{r}_{m}=J_{\ell}\mathbf{x}_{m}-\mathbf{b}$, i.e.,
$V_{m}^{T}\mathbf{r}_{m}=\mathbf{0}$.

Some computational aspects and the effectiveness of any projection method
strongly depend on the adopted approximation space $\mathcal{K}_{m}$.  Even
though the structure of $J_{\ell}$ in (8) is rather involved – it is a linear
combination of inverses of shifted and scaled $K$’s – we propose using the
polynomial Krylov subspace

(9)
$\mathcal{K}_{m}={\color[rgb]{0,0,0}\mathcal{K}_{m}(K,\mathbf{b})=}\text{span}\\{\mathbf{b},K\mathbf{b},\ldots,K^{m-1}\mathbf{b}\\},$

so that our projection method can be seen as a FOM-like scheme333Our
projection scheme does not amount to the actual FOM method as the matrix used
to define the subspace $\mathcal{K}_{m}$, namely $K$, does not coincide with
the coefficient matrix of the linear system we want to solve, namely
$J_{\ell}$.. However, other options such as rational Krylov subspaces can be
used instead of (9) as well.  In our numerical experiments, using (9) already
leads to remarkable performance, especially in terms of number of iterations;
see section 6. This means that using more sophisticated spaces, with a more
expensive construction step, may not reap benefits.

We would like to mention here that the FOM-like method we propose for solving
(8) can indeed be seen as a Krylov method for the numerical approximation of
$f(K)\mathbf{b}$, with $f$ being the rational function
$f(K)=\left(I+\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}\right)^{-1}\mathbf{b}$. Further considerations about the connection
between these two points of view will be made at the end of this section when
$K$ is SPD.

The use of the approximation space (9) allows us to compute the vector
$\mathbf{y}_{m}$ and derive a relation for a cheap computation of the residual
norm $\|\mathbf{r}_{m}\|$, using the Arnoldi relation arising from the
construction of the Krylov subspace [29],

(10) $KV_{m}=V_{m}T_{m}+t_{m+1,m}\mathbf{v}_{m+1}\mathbf{e}_{m}^{T},$

where $T_{m}=V_{m}^{T}KV_{m}$, $t_{m+1,m}$ stems from the orthogonalization
procedure, and $\mathbf{v}_{m+1}$ is the $(m+1)$-st basis vector. Since the
relation (10) is shift-invariant, we can shift $K$ using any eigenvalue
$\pi_{i}$ of the circulant matrix, $C_{1}$, giving

$((1-\pi_{i})I_{\bar{n}}+\tau K)V_{m}=V_{m}((1-\pi_{i})I_{m}+\tau
T_{m})+t_{m+1,m}\mathbf{v}_{m+1}\mathbf{e}_{m}^{T},\quad\forall\;i=1,\ldots,\ell.$

Moreover, by premultiplying by $((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}$,
postmultiplying by $((1-\pi_{i}I_{m})+\tau T_{m})^{-1}$, and moving some terms
we get

$\displaystyle((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}V_{m}=$ $\displaystyle
V_{m}((1-\pi_{i})I_{m}+\tau T_{m})^{-1}$ (11) $\displaystyle-
t_{m+1,m}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}\mathbf{v}_{m+1}\mathbf{e}_{m}^{T}((1-\pi_{i})I_{m}+\tau T_{m})^{-1}.$

Looking at the residual vector, if $\beta:=\|\mathbf{b}\|$,
$\mathbf{h}_{i}:=((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}\mathbf{v}_{m+1}$, and
$S_{i}:=((1-\pi_{i})I_{m}+\tau T_{m})^{-1}$, we have

$\displaystyle\mathbf{r}_{m}=$ $\displaystyle
J_{\ell}V_{m}\mathbf{y}_{m}-\mathbf{b}$ $\displaystyle=$ $\displaystyle
V_{m}\left(\left(I_{m}+\sum_{i=1}^{\ell}\gamma_{i}S_{i}-t_{m+1,m}V_{m}^{T}\sum_{i=1}^{\ell}\gamma_{i}\mathbf{h}_{i}\mathbf{e}_{m}^{T}S_{i}\right)\mathbf{y}_{m}-\beta\mathbf{e}_{1}\right)$
$\displaystyle-
t_{m+1,m}(I-V_{m}V_{m}^{T})\sum_{i=1}^{\ell}\gamma_{i}\mathbf{h}_{i}\mathbf{e}_{m}^{T}S_{i}\mathbf{y}_{m}.$

Imposing the Galerkin condition $V_{m}^{T}\mathbf{r}_{m}=0$ is thus equivalent
to computing $\mathbf{y}_{m}$ as the solution of the $m\times m$ linear system

(12)
$\left(I_{m}+\sum_{i=1}^{\ell}\gamma_{i}\left(I_{m}-t_{m+1,m}V_{m}^{T}\mathbf{h}_{i}\mathbf{e}_{m}^{T}\right)S_{i}\right)\mathbf{y}_{m}=\beta\mathbf{e}_{1}.$

Then the residual norm is such that

(13)
$\|\mathbf{r}_{m}\|=|t_{m+1,m}|\cdot\left\|(I-V_{m}V_{m}^{T})\sum_{i=1}^{\ell}\gamma_{i}\mathbf{h}_{i}\mathbf{e}_{m}^{T}S_{i}\mathbf{y}_{m}\right\|.$

Our FOM-like method for (8) is summarized in Algorithm 3.1.

1

input : $K\in\mathbb{R}^{\bar{n}\times\bar{n}}$,
$\mathbf{b}\in\mathbb{R}^{\bar{n}}$, $\pi_{i}$,$\gamma_{i}\in\mathbb{C}$,
$i=1,\ldots,\ell$, $\text{{\rm maxit}}\in\mathbb{N}$, $\epsilon>0$,
${\color[rgb]{0,0,0}q}\geq 1$.

output : $\mathbf{x}_{m}\in\mathbb{R}^{\bar{n}}$ approximate solution to (8).

2 Set $\beta=\|\mathbf{b}\|$, $\mathbf{v}_{1}=\mathbf{b}/\beta$,
$V_{1}=\mathbf{v}_{1}$, $m=0$, $\|\mathbf{r}\|=1$

3

4while  _$\|\mathbf{r}\| >\epsilon\cdot\beta$ and $m\leq\text{{\rm maxit}}$_
do

5 Set $m=m+1$

6

7 Compute $\widetilde{\mathbf{v}}=K\mathbf{v}_{m}$

8

9 for  _$k=1,\ldots,m$_ do

10 Compute $t_{k,m}=\mathbf{v}_{k}^{T}\widetilde{\mathbf{v}}$

11

12 Set $\widetilde{\mathbf{v}}=\widetilde{\mathbf{v}}-t_{k,m}\mathbf{v}_{k}$

13

14 Set $t_{m+1,m}=\|\widetilde{\mathbf{v}}\|$,
$\mathbf{v}_{m+1}=\widetilde{\mathbf{v}}/t_{m+1,m}$, and
$V_{m+1}=[V_{m},\mathbf{v}_{m+1}]$

15

16 if  _$\mathtt{mod}(m,{\color[rgb]{0,0,0}q})=0$_ then

17 parfor $i=1,\ldots,\ell$

18 Set $\mathbf{h}_{i}=((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}\mathbf{v}_{m+1}$
and $S_{i}=((1-\pi_{i})I_{m}+\tau T_{m})^{-1}$

19

20

21 Solve
$(I_{m}+\sum_{i=1}^{\ell}\gamma_{i}(I_{m}-t_{m+1,j}V_{m}^{T}\mathbf{h}_{i}\mathbf{e}_{j}^{T})S_{i})\mathbf{y}_{m}=\beta\mathbf{e}_{1}$

22

23 Compute
$\|\mathbf{r}\|=|t_{m+1,m}|\cdot\|(I-V_{m}V_{m}^{T})\sum_{i=1}^{\ell}\mathbf{h}_{i}\mathbf{e}_{m}^{T}S_{i}\mathbf{y}_{m}\|$

24

25

26

27Set $\mathbf{x}_{m}=V_{m}\mathbf{y}_{m}$

Algorithm 3.1 FOM-like method for (8)

The most computationally demanding step of this Krylov method is the residual
norm computation. In particular, the computation of the vectors
$\mathbf{h}_{i}$ requires the parallel solution of the linear systems with
$(1-\pi_{i})I_{\bar{n}}+\tau K$ for all $i=1,\ldots,\ell$. In order to reduce
the computational cost, we may want to solve (12) and compute (13) only every
${\color[rgb]{0,0,0}q}\geq 1$ iterations, namely the residual norm gets
_frozen_ for ${\color[rgb]{0,0,0}q}$ iterations. In the worst case scenario,
this procedure leads to computing a slightly larger subspace than what would
have been necessary by checking (13) at each iteration. On the other hand, the
overall number of parallel-in-time loops performed by Algorithm 2.1 becomes
$m/{\color[rgb]{0,0,0}q}+2$.

Another option could be using an _inner_ Krylov method to compute the
$\mathbf{h}_{i}$’s. Since $((1-\pi_{i})I_{\bar{n}}+\tau
K)\mathbf{h}_{i}=\mathbf{v}_{m+1}$, by exploiting the fact the the right-hand
side $\mathbf{v}_{m+1}$ does not depend on the shift index $i$ along with the
shift-invariance property of the Krylov subspace, one may want to construct
$\mathcal{K}_{t}(K,\mathbf{v}_{m+1})$ and employ well-established Krylov
routines for shifted linear systems; see, e.g. [3, 30]. On the other hand,
this procedure would compute only approximations
$\mathbf{\widetilde{h}}_{i}\approx\mathbf{h}_{i}$ in general. The impact of
such approximations on the vector $\mathbf{y}_{m}$ and, ultimately, on
$\mathbf{x}_{m}$ may be tricky to assess. Moreover, an underlying assumption
of this paper is that we are able to solve $\ell$ shifted linear systems with
$K$ by using parallelization.

We now consider in more detail the special case where the stiffness matrix $K$
is SPD, and we show that our algorithm for (8) can largely take advantage of
this structure.  See the Appendix for the proof of the following result.

######  Theorem 3.1.

Let the spatial stiffness matrix $K$ be symmetric positive definite. Then the
coefficient matrix
$J_{\ell}:=I+\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}$
in (8) is Hermitian positive definite for any $\ell\geq 1$. Moreover,

(14) $\kappa(J_{\ell})\leq 1+\frac{1}{\tau\lambda_{\min}(K)}.$

Recalling that Algorithm 3.1 imposes a Galerkin condition on the Krylov
subspace $\mathcal{K}_{m}$, we point out that Theorem 3.1 implies that the
solution $\mathbf{x}_{m}$ provided by Algorithm 3.1 minimizes the error in the
$J_{\ell}$-energy norm over $\mathcal{K}_{m}$ whenever $K$ is SPD.

It is well-known that the convergence rate of FOM-like methods for symmetric
positive definite problems is related to the condition number of the
coefficient matrix444Though better insight is provided by the entire
eigenvalue distribution., see, e.g., [21]. The bound (14) displays the
interplay between the spatial and time discretization and how they contribute
to the convergence of Algorithm 3.1. In particular, if $\lambda_{\min}(K)$ is
far from zero and the time grid is rather coarse, i.e. $\tau$ is large, we
expect Algorithm 3.1 to converge fast. On the other hand, we may need many
iterations to reach the desired level of accuracy for problems posed on very
fine time grids with a small $\lambda_{\min}(K)$, depending on the scaling of
$K$, e.g., $1/h^{2}$ or $1/h$.

Since $K$ is SPD, $\lambda_{\min}(K)$ can be cheaply computed by, e.g., the
inverse power method555The implementation of such a method must make use of
the properties of $K$ to be efficient. In particular, its symmetric positive
definite nature and the possible sparsity coming from the adopted
discretization in space must be taken into account.. Therefore, it is easy to
check the magnitude of $\tau\lambda_{\min}(K)$. On the other hand, one may
want to select a time grid such that the latter value is big enough to
guarantee a well-conditioned $J_{\ell}$.

For $K$ SPD, the solution of (12) by Algorithm 3.1 can be seen as the
numerical evaluation of the action of the rational matrix function
$f(K)=(I+\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1})^{-1}$ on the vector $\mathbf{b}$, namely $f(K)\mathbf{b}$, by the
Lanczos method. Therefore, one may want to take advantage of this viewpoint to
derive a-priori upper bounds on the 2- and/or $J_{\ell}$-norm of the error;
see, e.g., [6, 4]. With these bounds at hand, it would be possible to predict
how many iterations are indeed sufficient to get the desired level of accuracy
in Algorithm 3.1, thus avoiding the expensive residual norm computation in
(13). In this scenario, a single parallel-in-time loop is necessary within
Algorithm 3.1. In particular, once a sufficiently large Krylov subspace is
constructed, this parallel-in-time loop is involved in the definition of the
coefficient matrix in (12) and thus in the computation of the vector
$\mathbf{y}_{m}$.

##  4 $\alpha$-acceleration

In the context of preconditioning operators for all-at-once linear systems
stemming from (1), $\alpha$-circulant matrices have been largely used; see,
e.g., [23]. We now explore the impact of this technique also on the scheme we
proposed in the previous sections. In particular, given $\alpha\in(0,1]$, we
write

(15) $\Sigma_{1}=C_{\alpha}-\alpha\mathbf{e}_{1}\mathbf{e}_{\ell}^{T},\quad
C_{\alpha}=\begin{bmatrix}0&&&\alpha\\\ 1&\ddots&&\\\ &\ddots&\ddots&\\\
&&1&0\\\ \end{bmatrix}\in\mathbb{R}^{\ell\times\ell}.$

Now $C_{\alpha}$ is an $\alpha$-circulant matrix, and it can be diagonalized
by the _scaled_ Fast Fourier Transform

$C_{\alpha}=D_{\alpha}^{-1}{F}^{-1}\Pi_{\alpha}{F}D_{\alpha},\qquad
D_{\alpha}=\begin{bmatrix}1&&&\\\ &\alpha^{1/\ell}&&\\\ &&\ddots&\\\
&&&\alpha^{(\ell-1)/\ell}\\\
\end{bmatrix},\quad\Pi_{\alpha}=\alpha^{1/\ell}\Pi_{1},$

see, e.g., [26].

It is well-known that the eigenvector matrix ${F}D_{\alpha}$ can be very ill-
conditioned for small $\alpha$ and sizable values of $\ell$. Indeed,
$\kappa({F}D_{\alpha})=\alpha^{-(\ell-1)/\ell}$. Such ill-conditioning will be
one of the major obstacles in using very small values of $\alpha$.

The use of (15) does not lead to any particular difficulty in the solver we
proposed in section 2. By following the same exact steps as in section 2 and
adopting the same notation, a direct computation shows that (2) translates
into

$\displaystyle\text{vec}(\widetilde{U})=$ $\displaystyle
P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T})$
(16)
$\displaystyle-\alpha^{1/\ell}P^{-1}M(I+\alpha^{1/\ell}N^{T}P^{-1}M)^{-1}N^{T}P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T}).$

Once $\widetilde{U}$ is computed, the original solution can be retrieved by
computing $U=\widetilde{U}{F}^{-T}D_{\alpha}^{-1}$.

From (4), it is clear how the use of $\alpha$-circulant matrices is equivalent
to introducing a _weight_ in our setting. Indeed, the scalar $\alpha\in(0,1]$
determines the contribution of the correction term

$\displaystyle\text{vec}(U_{2})=$
$\displaystyle\alpha^{1/\ell}(D_{\alpha}^{-1}{F}^{-1}\otimes
I)P^{-1}M(I+\alpha^{1/\ell}N^{T}P^{-1}M)^{-1}N^{T}P^{-1}$
$\displaystyle\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T}),$

to the point that, for sufficiently small $\alpha$, the solution can be often
well approximated even when neglecting this term; see section 6.

The inner FOM-like method can benefit from the introduction of the parameter
$\alpha$ as well. Indeed, it holds that

(17)
$I+\alpha^{1/\ell}N^{T}P^{-1}M=I+\alpha^{1/\ell}\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}.$

Therefore, we need to solve a linear system whose coefficient matrix can be
seen as a small perturbation of the identity for $\alpha\ll 1$. Using small
values of $\alpha$ has thus the potential of remarkably reducing the number of
iterations performed by the FOM-like method to attain a prescribed accuracy.
This motivates the name _$\alpha$ -acceleration_ whenever $\alpha$-circulant
matrices of the form (15) are used in our context.

For the sake of completeness we report in Algorithm 4.1 the overall procedure
when using the $\alpha$-circulant matrix in (15). In line 4.1 we compute the
residual norm provided by the first term of the solution, namely
$\text{vec}(U_{1})=(D_{\alpha}^{-1}{F}^{-1}\otimes
I)P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T})$.
If this residual norm is sufficiently small, we stop the algorithm and set
$U=U_{1}$ thus saving a lot of computational effort.

1

input : $K\in\mathbb{R}^{\bar{n}\times\bar{n}}$, $\mathbf{u}_{0}$,
$\mathbf{f}_{i}\in\mathbb{R}^{\bar{n}}$, $i=1,\ldots,\ell$,
$\ell\in\mathbb{N}$, $\alpha\in(0,1]$, $\tau,\epsilon>0$.

output : $U\in\mathbb{R}^{\bar{n}\times\ell}$ approximate solution to (2).

2 Compute
$[\pi_{1},\ldots,\pi_{\ell}]^{T}=\alpha^{1/\ell}{F}C_{1}\mathbf{e}_{1}$,
$[\gamma_{1},\ldots,\gamma_{\ell}]^{T}={F}^{-T}\mathbf{e}_{\ell}$,
$D_{\alpha}=\text{diag}(1,\alpha^{1/\ell},\ldots,\alpha^{(\ell-1)/\ell})$

3

4 parfor $i=1,\ldots,\ell$

5 Set $L\mathbf{e}_{i}=((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}[\mathbf{u}_{0}+\tau\mathbf{f}_{1},\tau\mathbf{f}_{2},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T}\mathbf{e}_{i}$

6

7

8Set $U_{1}=L{F}^{-T}D_{\alpha}^{-1}$

9

10if  _$\|(I+\tau
K)U_{1}-U_{1}\Sigma_{1}^{T}-[\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]\|_{F}\leq\epsilon\cdot\|U_{1}\|_{F}$_
then

11

12 Set $U=U_{1}$ and return

13

14Compute $\mathbf{b}=L{F}^{-T}\mathbf{e}_{\ell}$

15

16Compute $\mathbf{x}_{m}$ by applying a variant of Algorithm 3.1 to
$\left(I+\alpha^{1/\ell}\sum_{i=1}^{\ell}\gamma_{i}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}\right)\mathbf{x}=\mathbf{b}$

17 parfor $i=1,\ldots,\ell$

18 Set $W\mathbf{e}_{i}=((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}\mathbf{x}_{m}$

19

20

Set $U=U_{1}-\alpha^{1/\ell}W{F}^{-T}D_{\alpha}^{-1}$

Algorithm 4.1 ParaDiag and Backward Euler with $\alpha$-acceleration

Even though the use of tiny values of $\alpha$ would be desirable, this would
lead to an extremely ill-conditioned eigenvector matrix ${F}D_{\alpha}$. For
very small values of $\alpha$, this poor conditioning pollutes the computed
solution $U$. In particular, it can be observed that the error is proportional
to $\kappa({F}D_{\alpha})$. This drawback is well-known also in the case of
the variable time-stepping procedure proposed in [11, 10]. Finding the value
of $\alpha$ providing the best trade-off between acceleration and loss of
accuracy is not an easy task. We refer to our numerical experiments to supply
some insight to this important issue.

##  5 Higher-order BDFs

In the previous sections we assumed that the backward Euler method was used
for time discretization so that $\Sigma_{1}$ was of the form (3). In this
section we show how to generalize our methodology  to higher-order BDF
approximations of time-dependent problems (1).

The backward Euler method belongs to the larger class of backward
differentiation formulas (BDFs) for time-dependent problems (1). In
particular, the backward Euler method is a BDF of order one, and BDFs of order
$s>1$ are well-established time integrator schemes as well.

The matrix formulation of the discrete problem stemming from  an all-at-once
discretization of (1) using a BDF of order $s$ can be written as

(18) $(I+\tau\beta K)U-U\Sigma_{s}^{T}=G,$

where

$G=\left[\sum_{j=1}^{s}\alpha_{j}\mathbf{u}_{1-j}+\tau\beta\mathbf{f}_{1},\sum_{j=2}^{s}\alpha_{j}\mathbf{u}_{2-j}+\tau\beta\mathbf{f}_{2}\ldots,\alpha_{s}\mathbf{u}_{0}+\tau\beta\mathbf{f}_{s},\tau\beta\mathbf{f}_{s+1},\ldots\tau\beta\mathbf{f}_{\ell}\right],$

and

(19) $\Sigma_{s}=\begin{bmatrix}0&&&&&&\\\ \alpha_{1}&\ddots&&&&&\\\
\vdots&\ddots&&&&&\\\ \alpha_{s}&&\ddots&\ddots&&&\\\ 0&\ddots&&\ddots&&&\\\
\vdots&&\ddots&&\ddots&\ddots&\\\ 0&&\cdots&\alpha_{s}&\cdots&\alpha_{1}&0\\\
\end{bmatrix},$

see, e.g., [27].

In the equation above, $\mathbf{u}_{1-s},\ldots,\mathbf{u}_{0}$ denote the $s$
initial values required by the adopted $s$-order BDF. In Table 1 we report the
scalars $\alpha_{j}=\alpha_{j}(s)$ and $\beta=\beta(s)$ defining a BDF of
order $s$; see, e.g., [1, Table 5.3]666Notice that in Table 1 we have changed
the sign of the $\alpha_{j}$’s with respect to the values listed in [1, Table
5.3] in order to have a term of the form $-U\Sigma_{s}^{T}$ also in (18) as it
is done for equation (2)..

Table 1:  BDF coefficients for $s\leq 6$. $s$ | $\beta$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\alpha_{4}$ | $\alpha_{5}$ | $\alpha_{6}$
---|---|---|---|---|---|---|---
1 | 1 | 1 |  |  |  |  |
2 | 2/3 | 4/3 | -1/3 |  |  |  |
3 | 6/11 | 18/11 | -9/11 | 2/11 |  |  |
4 | 12/25 | 48/25 | -36/25 | 16/25 | -3/25 |  |
5 | 60/137 | 300/137 | -300/137 | 200/137 | -75/137 | 12/137 |
6 | 60/147 | 360/147 | -450/147 | 400/147 | -225/147 | 72/147 | -10/147

It has been shown in, e.g., [1, Section 5.2.3] that any BDF of order $s>6$ is
unstable. We thus restrict ourselves to $s\leq 6$.

For a general BDF of order $s$, the time discrete operator $\Sigma_{s}$ can be
written as a circulant plus a matrix of rank $s$, namely

$\Sigma_{s}=C_{s}-[\mathbf{e}_{1},\ldots,\mathbf{e}_{s}]\bm{\alpha}_{s}[\mathbf{e}_{\ell-s+1},\ldots,\mathbf{e}_{\ell}]^{T},\
\bm{\alpha}_{s}=\begin{bmatrix}\alpha_{s}&\ldots&&\ldots&\alpha_{1}\\\
&\alpha_{s}&\ldots&\ldots&\alpha_{1}\\\ &&\ddots&&\vdots\\\
&&&\alpha_{s}&\alpha_{s-1}\\\ &&&&\alpha_{s}\\\
\end{bmatrix}\in\mathbb{R}^{s\times s}.$

This means that the same strategy presented in section 2 can be applied by
setting in (2) $P:=I_{\ell}\otimes(I_{\bar{n}}+\tau\beta K)-\Pi_{s}\otimes
I_{\bar{n}}$, $C_{s}:={F}\Pi_{s}{F}^{-1}$,
$\Pi_{s}=\text{diag}(FC_{s}\mathbf{e}_{1})$,
$M:={F}[\mathbf{e}_{1},\ldots,\mathbf{e}_{s}]\otimes I_{\bar{n}}$, and
$N:={F}^{-T}[\mathbf{e}_{\ell-s+1},\ldots,\mathbf{e}_{\ell}]\bm{\alpha}_{s}^{T}$.
The main difference lies in the linear system (7). In this setting, the right-
hand side, simply denoted by $\mathbf{b}$ in (7), is given by
$\mathbf{b}^{(s)}=\text{vec}(B^{(s)})$ where

$B^{(s)}=[((1-\pi_{1})I_{\bar{n}}+\tau\beta
K)^{-1}G{F}^{T}\mathbf{e}_{1},\ldots,((1-\pi_{\ell})I_{\bar{n}}+\tau\beta
K)^{-1}G{F}^{T}\mathbf{e}_{\ell}]{F}^{-T}[\mathbf{e}_{\ell-s+1},\ldots,\mathbf{e}_{\ell}]\bm{\alpha}_{s}^{T}.$

Moreover, the coefficient matrix $I+N^{T}P^{-1}M\in\mathbb{R}^{s\bar{n}\times
s\bar{n}}$ can be seen as an $s\times s$ block matrix with
$\bar{n}\times\bar{n}$ blocks; see, e.g., [27, Section 5.2]. If

${F}\mathbf{e}_{i}=\begin{bmatrix}\theta_{1}^{(i)}\\\ \vdots\\\
\theta_{\ell}^{(i)}\\\
\end{bmatrix},\quad{F}^{-T}[\mathbf{e}_{\ell-s+1},\ldots,\mathbf{e}_{\ell}]\bm{\alpha}_{s}\mathbf{e}_{i}=\begin{bmatrix}\gamma_{1}^{(i)}\\\
\vdots\\\ \gamma_{\ell}^{(i)}\\\ \end{bmatrix},\quad i=1,\ldots,s,$

so that

$M=\begin{bmatrix}\theta_{1}^{(1)}I_{\bar{n}}&\cdots&\theta_{1}^{(s)}I_{\bar{n}}\\\
\vdots&&\vdots\\\
\theta_{\ell}^{(1)}I_{\bar{n}}&\cdots&\theta_{\ell}^{(s)}I_{\bar{n}}\\\
\end{bmatrix},\quad
N=\begin{bmatrix}\gamma_{1}^{(1)}I_{\bar{n}}&\cdots&\gamma_{1}^{(s)}I_{\bar{n}}\\\
\vdots&&\vdots\\\
\gamma_{\ell}^{(1)}I_{\bar{n}}&\cdots&\gamma_{\ell}^{(s)}I_{\bar{n}}\\\
\end{bmatrix},$

then we have

$N^{T}P^{-1}M=\begin{bmatrix}\displaystyle\sum_{i=1}^{\ell}\gamma_{i}^{(1)}\theta_{i}^{(1)}((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}\par&\cdots&\displaystyle\sum_{i=1}^{\ell}\gamma_{i}^{(1)}\theta_{i}^{(s)}((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}\\\ \vdots&&\vdots\\\
\displaystyle\sum_{i=1}^{\ell}\gamma_{i}^{(s)}\theta_{i}^{(1)}((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}&\cdots&\displaystyle\sum_{i=1}^{\ell}\gamma_{i}^{(s)}\theta_{i}^{(s)}((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}\\\ \end{bmatrix}.$

Denoting the matrix $I+N^{T}P^{-1}M$ by $J_{\ell}^{(s)}$, we thus need to
solve the linear system

(20) $J_{\ell}^{(s)}\mathbf{x}=\mathbf{b}^{(s)}.$

To this end, we adopt a projection scheme similar to the one presented in
section 3. Instead of (9), we use the block Krylov subspace

(21)
$\mathcal{K}_{m}^{\square}=\text{range}\\{[B^{(s)},KB^{(s)},\ldots,K^{m-1}B^{(s)}]\\}.$

If $V_{m}=[\mathcal{V}_{1},\ldots,\mathcal{V}_{m}]\in\mathbb{R}^{n\times sm}$,
$\mathcal{V}_{i}\in\mathbb{R}^{n\times s}$ for all $i=1,\ldots,m$, has
orthonormal columns and is such that
$\text{range}(V_{m})=\mathcal{K}_{m}^{\square}$, then the block counterpart of
the Arnoldi relation (10) holds. In particular,

(22) $KV_{m}=V_{m}T_{m}+\mathcal{V}_{m+1}t_{m+1,m}(\mathbf{e}_{m}^{T}\otimes
I_{s}),$

where $T_{m}=V_{m}^{T}KV_{m}\in\mathbb{R}^{sm\times sm}$, and
$t_{m+1,m}\in\mathbb{R}^{s\times s}$; see, e.g., [17, 20] for further details
on block Krylov methods and efficient schemes for the computation of the basis
$V_{m}$.

In view of (22) we can derive the block counterparts of (3). In particular, by
using the same arguments of section 3, we have

(23) $\displaystyle((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}V_{m}=V_{m}((1-\pi_{i})I_{m}+\tau\beta T_{m})^{-1}$
$\displaystyle\qquad-((1-\pi_{i})I_{\bar{n}}+\tau\beta
K)^{-1}\mathcal{V}_{m+1}t_{m+1,m}(\mathbf{e}_{m}^{T}\otimes
I_{s})((1-\pi_{i})I_{m}+\tau\beta T_{m})^{-1}.$

We propose to solve (20) by imposing a Galerkin condition with respect to the
space spanned by $I_{s}\otimes V_{m}$. In particular, if
$\mathbf{x}_{m}=(I_{s}\otimes V_{m})\mathbf{y}_{m}$, we want the residual
vector $\mathbf{r}_{m}=J_{\ell}^{(s)}\mathbf{x}_{m}-\mathbf{b}^{(s)}$ to be
such that $(I_{s}\otimes V_{m}^{T})\mathbf{r}_{m}=0$.

######  Proposition 5.1.

If $B^{(s)}=V_{1}\beta$, $\beta\in\mathbb{R}^{s\times s}$,
$H_{i}=((1-\pi_{i})I_{\bar{n}}+\tau\beta K)^{-1}\mathcal{V}_{m+1}$, and
$S_{i}=((1-\pi_{i})I_{m}+\tau\beta T_{m})^{-1}$, then the residual vector
$\mathbf{r}_{m}=J_{\ell}^{(s)}\mathbf{x}_{m}-\mathbf{b}^{(s)}=J_{\ell}^{(s)}(I_{s}\otimes
V_{m})\mathbf{y}_{m}-\mathbf{b}^{(s)}$ can be written as

$\displaystyle\mathbf{r}_{m}=$ $\displaystyle(I_{s}\otimes
V_{m})\left(\left(I_{sm}+\mathcal{T}_{m}\right)\mathbf{y}_{m}-(\mathbf{e}_{1}\otimes
I_{s})\beta\right)-(I-I_{s}\otimes
V_{m}V_{m}^{T})\mathcal{S}_{m}\mathbf{y}_{m},$

where $\mathcal{T}_{m},\mathcal{S}_{m}\in\mathbb{R}^{ms\times ms}$ are
$s\times s$ block matrices whose $(k,h)$ blocks $(\mathcal{T}_{m})_{k,h}$,
$(\mathcal{S}_{m})_{k,h}$ are given by

$(\mathcal{T}_{m})_{k,h}=\sum_{i=1}^{\ell}\gamma_{i}^{(k)}\theta_{i}^{(h)}S_{i}-V_{m}^{T}\sum_{i=1}^{\ell}\gamma_{i}^{(k)}\theta_{i}^{(h)}H_{i}t_{m+1,m}(\mathbf{e}_{m}^{T}\otimes
I_{s})S_{i},$
$(\mathcal{S}_{m})_{k,h}=\sum_{i=1}^{\ell}\gamma_{i}^{(k)}\theta_{i}^{(h)}H_{i}t_{m+1,m}(\mathbf{e}_{m}^{T}\otimes
I_{s})S_{i}.$

Moreover, imposing the Galerkin condition
$\mathbf{r}_{m}\perp\mathcal{K}_{m}^{\square}$ is equivalent to computing
$\mathbf{y}_{m}$ as the solution of the linear system

$\left(I_{sm}+\mathcal{T}_{m}\right)\mathbf{y}_{m}=(\mathbf{e}_{1}\otimes
I_{s})\beta,$

and

$\|\mathbf{r}_{m}\|=\|(I-I_{s}\otimes
V_{m}V_{m}^{T})\mathcal{S}_{m}\mathbf{y}_{m}\|.$

######  Proof 5.2.

The results can be shown by applying the exact same arguments of section 3
_block-wise_ noticing that

$I_{s}\otimes V_{m}=\begin{bmatrix}V_{m}&&\\\ &\ddots&\\\ &&V_{m}\\\
\end{bmatrix}.$

##  6 Numerical examples

In this section we show several numerical examples to illustrate different
aspects of our new ParaDiag algorithm. We would like to mention that in the
following we refrain from reporting any running times achieved by the routines
we test. Indeed, all the experiments are carried out on a simple laptop, so
that reporting running times would be only informative in a relative sense.
Moreover, the application of $P$ is still performed sequentially and no real
parallel-in-time implementation is adopted so far. On the other hand, we
report the number of parallel-in-time loops that would be performed by the
tested schemes as we believe this can be a fair measure of the computational
cost.

Unless stated otherwise, the threshold used in Algorithm 3.1 is
$\epsilon=10^{-8}$. Moreover, we check the residual norm in Algorithm 3.1 at
each iteration, i.e. we set ${\color[rgb]{0,0,0}q}=1$.

###  6.1 Heat Equation

We consider the heat equation

(24) $\left\\{\begin{array}[]{rlll}u_{t}&=&\Delta u+f,&\quad\text{in
}\Omega\times(0,1],\\\ u&=&0,&\quad\text{on }\partial\Omega,\\\
u(x,0)&=&u_{0},&\end{array}\right.$

where $\Omega=(0,1)^{2}$. If the $2$-dimensional Laplace operator is
discretized using second order centered finite differences with $n$ nodes in
each direction and $h=1/(n+1)$, the stiffness matrix
$K\in\mathbb{R}^{\bar{n}\times\bar{n}}$, $\bar{n}=n^{2}$, can be written as

${\color[rgb]{0,0,0}K=I_{n}\otimes T+T\otimes I_{n}},\quad
T=\frac{1}{h^{2}}\begin{pmatrix}2&-1&&\\\ -1&\ddots&\ddots&\\\
&\ddots&\ddots&-1\\\ &&-1&2\\\ \end{pmatrix}\in\mathbb{R}^{n\times n}.$

Since the eigenvalues of the SPD matrix $T$ are known in closed-form, by using
the properties of the Kronecker product, it is easy to show that

$\lambda_{\min}(K)=\frac{8}{h^{2}}\sin\left(\frac{\pi}{2(n+1)}\right)^{2}.$

By defining $\tau=1/\ell$ for a given number of time steps $\ell\geq 1$, we
can compute the upper bound (14). In Table 2 we show this upper bound for
common values of $\bar{n}$ and $\ell$.

Table 2: Upper bound  of the coefficient matrix, $J_{\ell}$, (14) for the $2$-dimensional heat equation (24). Different values of the number of spatial degrees of freedom $\bar{n}$ and time steps $\ell$ are used in the discretization. $\bar{n}$ | $\ell$ | $1+\frac{1}{\tau\lambda_{\min}(K)}$
---|---|---
65536 | 256 | 13.969
512 | 26.938
1024 | 52.877
262144 | 256 | 13.969
512 | 26.938
1024 | 52.876
1048576 | 256 | 13.969
512 | 26.938
1024 | 52.876

From the results in Table 2 we first notice that while the value of the upper
bound (14) linearly grows with $\ell$, an increment in the number of spatial
degrees of freedom $\bar{n}$ does not have a significant impact since
$\lambda_{\min}(K)$ only moderately increases with $n$. Moreover,
$1+1/(\tau\lambda_{\min}(K))$ is $\mathcal{O}(10)$ for any value of $\bar{n}$
and $\ell$ we tested. This implies that the matrix $J_{\ell}$ is always very
well-conditioned regardless of the source term $f$ and the initial condition
$u_{0}$ in (24).

We tested Algorithm 2.1 on an instance of the heat equation (24). In
particular, we consider [25, Example 6.1] where $f=0$ and
$u_{0}=xy(x-1)(1-y)$.  For all the values of $\bar{n}$ and $\ell$ listed in
Table 2, a single iteration of Algorithm 3.1 already results in a relative
residual norm smaller than $\epsilon=10^{-8}$. Indeed, a more careful
inspection of the eigenvalue distribution of $J_{\ell}$ shows that the latter
is a very small perturbation of the identity. In particular,

$\lambda_{\min}(J_{\ell})=1+\sum_{i=1}^{\ell}\frac{\gamma_{i}}{1-\pi_{i}+\tau\lambda_{\max}(K)},\quad\lambda_{\max}(J_{\ell})=1+\sum_{i=1}^{\ell}\frac{\gamma_{i}}{1-\pi_{i}+\tau\lambda_{\min}(K)},$

and for all the values of $\bar{n}$ and $\ell$ listed in Table 2 we have
$[\lambda_{\min}(J_{\ell}),\lambda_{\max}(J_{\ell})]\subset[1,1.004]$.

###  6.2 Advection-Diffusion Problem

We now consider the time-dependent advection-diffusion equation

(25) $\left\\{\begin{array}[]{rlll}u_{t}-\nu\Delta u+\mathbf{w}\cdot\nabla
u&=&0,&\text{in }\Omega\times(0,1],\;\Omega\colon=(0,1)^{2},\\\
u&=&g(x,y),&\text{on }\partial\Omega,\\\
u_{0}=u(x,y,0)&=&0&\text{otherwise,}\end{array}\right.$

where $\nu>0$, $\mathbf{w}=(2y(1-x^{2}),-2x(1-y^{2}))$ and
$g(1,y)=g(x,0)=g(x,1)=0$ while $g(0,y)=1$; see, e.g., [25]. To obtain the
stiffness matrix $K\in\mathbb{R}^{\bar{n}\times\bar{n}}$, we used again
centered finite differences, and for the time integration backward Euler. The
initial vector $\mathbf{u}_{0}$ is set to be zero everywhere except the
boundaries, where it satisfies the boundary conditions.

We start by exploring numerically the $\alpha$-acceleration technique
presented in section 4. We show in Figure 1

Figure 1: Advection-diffusion equation: $\|U_{1}\|_{F}$ (solid line with
stars), $\|U_{2}\|_{F}$ (solid line with circles), relative residual norm
achieved by $U$ (dashed-dotted line) and $U_{1}$ (dotted line), and the error
$\|\widehat{U}-U\|_{F}/\|\widehat{U}\|_{F}$ (dashed line) for different values
of $\alpha$, $\bar{n}=65\,536$, $\nu=2^{-5}$, and $\ell=64$. $\widehat{U}$
denotes the solution obtained by running GMRES on the all-at-once system with
a small residual tolerance ($10^{-13}$).

the results obtained by running a single iteration of the FOM-like method,
namely $\text{maxit}=1$ in Algorithm 3.1, for different values of $\alpha$. We
set $\bar{n}=65\,536$, $\nu=2^{-5}$, and use $\ell=64$ time steps, and recall
that we write the computed solution $U$ as $U=U_{1}+U_{2}$ where $U_{1}$ and
$U_{2}$ are such that

$\begin{array}[]{l}\text{vec}(U_{1})=(D_{\alpha}^{-1}{F}^{-1}\otimes
I)\text{vec}(L),\quad\text{vec}(L)=P^{-1}\text{vec}([\mathbf{u_{0}}+\tau\mathbf{f}_{1},\ldots,\tau\mathbf{f}_{\ell}]D_{\alpha}{F}^{T}),\\\
\\\ \text{vec}(U_{2})=-\alpha^{1/\ell}(D_{\alpha}^{-1}{F}^{-1}\otimes
I)P^{-1}M(I+\alpha^{1/\ell}N^{T}P^{-1}M)^{-1}N^{T}\text{vec}(L).\end{array}$

In Figure 1 we show $\|U_{1}\|_{F}$ (solid line with stars), $\|U_{2}\|_{F}$
(solid line with circles), the relative residual norm achieved by the whole
$U$ (dashed-dotted line) and the one attained by using only $U_{1}$ (dotted
line) for different values of $\alpha$. We also computed the _exact_ algebraic
solution $\widehat{U}$ by running GMRES on the all-at-once system with the
very small tolerance $10^{-13}$ on the relative residual norm. In Figure 1 we
thus plot also the trend of $\|\widehat{U}-U\|_{F}/\|\widehat{U}\|_{F}$
(dashed line).

The first thing to notice from Figure 1 is that a change in $\alpha$ does not
really have an impact on $\|U_{1}\|_{F}$; this value remains (almost) constant
for all the $\alpha$’s we tested. On the other hand, the trend of
$\|U_{2}\|_{F}$ closely follows the values of $\alpha$ showing how the
contribution of $U_{2}$ strongly depends on the selected $\alpha$. If a
$\mathcal{O}(10^{-7})$ relative residual norm is good enough777Notice that
having this kind of accuracy in the algebraic problem is often exceeding the
discretization error by orders of magnitude., $U_{2}$ can be completely
neglected if $\alpha=\mathcal{O}(10^{-7})$ in this example. This means that
$U=U_{1}$ can be computed by performing a single parallel-in-time loop.  On
the other hand, $U_{2}$ has an important role in the accuracy of the overall
solution. Indeed, it is interesting to note that the error and the residual
norm achieved by $U=U_{1}+U_{2}$ scale like $\alpha^{2}$ whereas the residual
norm attained by $U_{1}$ depends only linearly on $\alpha$.

Note also that $\alpha$ plays a role in the inner linear system we need to
solve to compute $U_{2}$. Figure 1 shows that the smaller $\alpha$, the more
accurate a solution we get by performing a single iteration of the FOM-like
method. This is due to the fact that the coefficient matrix in (17) becomes a
small perturbation of the identity for small values of $\alpha$.

As already mentioned, however, using too small values of $\alpha$ leads to a
remarkable increase in the condition number of the transformation matrix
${F}D_{\alpha}$. This is clearly visible in Figure 1: for $\alpha\leq 10^{-4}$
the error starts increasing, since $\kappa({F}D_{\alpha})$ is becoming the
dominant factor polluting the quality of the computed solution $U$.

###  6.3 Comparison with other ParaDiag Methods

We now compare our new ParaDiag scheme, namely Algorithm 4.1, with different
state-of-the-art ParaDiag solvers. As already mentioned, in what follows we
report the number of parallel-in-time loops (#PinT) needed by the routines we
test as a measure to assess their computational cost, instead of running time,
since the latter strongly depends on the precise implementations and computing
infrastructures used. Recording fair comparisons in terms of running times
would not be straightforward.

The first ParaDiag technique we compare to is the one proposed in [25], where
preconditioned GMRES is used to iteratively solve the
$\bar{n}\ell\times\bar{n}\ell$ linear system corresponding to the Kronecker
form of (2).  The preconditioning operator we use within GMRES with right
preconditioning888To have fair comparisons, we use a matrix oriented GMRES
implementation, so that the Kronecker form of (2) is never explicitly
computed. See, e.g., [28].  is

$\mathcal{P}:x\mapsto(F\otimes I_{\bar{n}})(I_{\ell}\otimes(I_{\bar{n}}+\tau
K)+\Pi_{1}\otimes I_{\bar{n}})(F^{-1}\otimes I_{\bar{n}})x.$

Therefore, at each iteration we need to apply $\mathcal{P}^{-1}$ to the
current basis vector, namely we need to perform a parallel-in-time loop every
time we apply the preconditioner. This means that, in total, this routine
performs $p+1$ parallel-in-time loops if $p$ denotes the number of GMRES
iterations needed to converge. The threshold on the relative residual norm we
used in GMRES is $\epsilon=10^{-8}$.  Notice that the use of right
preconditioning allows us to have access to the actual, unpreconditioned
residual norm computed by GMRES. The latter quantity is thus comparable with
the residual norm achieved by Algorithm 4.1.

The second ParaDiag technique we compare to is the new interpolation scheme
presented in [19, Section 3] and denoted by Ev-Int in what follows. In this
method, given two parameters $\rho$ and $r$, one needs to perform $r$
parallel-in-time loops involving different coefficient matrices. The quality
of the computed solution depends on $\rho$ and $r$. The authors in [19] do not
comment much on the selection of $\rho$ and $r$. They suggest to use
$\rho=5\cdot 10^{-4}$ and $r=2$, values adopted in most of the experiments
shown in [19].

In Table 3

|  |  | Algorithm 4.1 ($\alpha=10^{-4}$) | GMRES | Ev-Int
---|---|---|---|---|---
$\bar{n}$ | $\ell$ | $\nu$ | #PinT | Rel. Res. | #PinT | Rel. Res. | #PinT | Rel. Res.
$16\,384$ | 32 | $10^{-1}$ | 3 | 8.41e-11 | 5 | 8.02e-13 | 2 | 7.26e-11
$10^{-2}$ | 3 | 6.98e-12 | 5 | 7.28e-12 | 2 | 1.24e-11
$10^{-3}$ | 3 | 2.77e-12 | 5 | 2.13e-10 | 2 | 7.22e-11
64 | $10^{-1}$ | 3 | 1.74e-11 | 5 | 1.12e-13 | 2 | 3.72e-11
$10^{-2}$ | 3 | 1.42e-12 | 5 | 1.31e-12 | 2 | 4.19e-12
$10^{-3}$ | 3 | 7.41e-13 | 5 | 9.42e-12 | 2 | 1.99e-11
128 | $10^{-1}$ | 3 | 1.20e-11 | 5 | 4.51e-14 | 2 | 2.30e-11
$10^{-2}$ | 3 | 1.01e-12 | 5 | 4.24e-13 | 2 | 2.23e-12
$10^{-3}$ | 3 | 3.99e-13 | 5 | 1.88e-12 | 2 | 8.11e-12
$65\,536$ | 32 | $10^{-1}$ | 3 | 3.42e-10 | 5 | 6.67e-13 | 2 | 8.58e-11
$10^{-2}$ | 3 | 2.72e-11 | 5 | 3.91e-12 | 2 | 1.09e-11
$10^{-3}$ | 3 | 3.61e-12 | 5 | 3.37e-11 | 2 | 4.49e-11
64 | $10^{-1}$ | 3 | 7.26e-11 | 5 | 1.31e-13 | 2 | 3.62e-11
$10^{-2}$ | 3 | 5.21e-12 | 5 | 6.73e-13 | 2 | 3.33e-12
$10^{-3}$ | 3 | 7.29e-13 | 5 | 8.71e-13 | 2 | 1.00e-11
128 | $10^{-1}$ | 3 | 4.84e-11 | 5 | 1.26e-13 | 2 | 3.07e-11
$10^{-2}$ | 3 | 3.71e-12 | 5 | 2.11e-13 | 2 | 2.28e-12
$10^{-3}$ | 3 | 5.40e-13 | 5 | 1.50e-13 | 2 | 3.55e-12
Table 3: Advection-Diffusion equation: results for different values of
$\bar{n}$, $\ell$, and $\nu$.

we show the results for different values of $\bar{n}$, $\ell$, and the
viscosity parameter $\nu$. We use Algorithm 4.1 with $\alpha=10^{-4}$, and
recall that this algorithm performs $m+2$ parallel-in-time loops where $m$ is
the number of iterations needed by Algorithm 3.1 to converge. From the results
in Table 3 we see that the number of parallel-in-time loops performed by both
Algorithm 4.1 and GMRES are very robust with respect to $\bar{n}$, $\ell$, and
$\nu$. We also see that the accuracy attained by our solver improves by
decreasing $\nu$, for fixed $\bar{n}$ and $\ell$, whereas GMRES shows the
opposite trend, in general. Ev-Int has a less regular behavior in this regard.

Notice however that, even though Algorithm 4.1 is often more accurate than Ev-
Int, the latter algorithm only performs two parallel-in-time loops, in
contrast to the three loops performed by our new ParaDiag algorithm. By
setting $r=3$, similar results in terms of accuracy can be obtained for Ev-Int
as well.

To make a fair comparison between our procedure and Ev-Int, in Table 4 we show

|  |  | #PinT$=2$ | #PinT$=1$
---|---|---|---|---
$\bar{n}$ | $\ell$ | $\nu$ | Alg. 4.1 ($\alpha=10^{-4}$) | Ev-Int | Alg. 4.1 ($\alpha=10^{-6}$) | Ev-Int
$16\,384$ | 32 | $10^{-1}$ | 8.44e-11 | 7.26e-11 | 1.87e-8 | 8.95e-6
$10^{-2}$ | 6.97e-12 | 1.24e-11 | 3.97e-8 | 1.98e-5
$10^{-3}$ | 3.14e-12 | 7.22e-11 | 5.88e-8 | 2.94e-5
64 | $10^{-1}$ | 1.76e-11 | 3.72e-11 | 1.09e-8 | 5.45e-6
$10^{-2}$ | 1.43e-12 | 4.19e-12 | 2.45e-8 | 1.22e-5
$10^{-3}$ | 8.43e-13 | 1.99e-11 | 3.99e-8 | 1.99e-5
128 | $10^{-1}$ | 1.20e-11 | 2.30e-11 | 8.23e-9 | 3.55e-6
$10^{-2}$ | 1.01e-12 | 2.23e-12 | 1.59e-8 | 7.96e-6
$10^{-3}$ | 4.30e-13 | 8.11e-12 | 2.79e-8 | 1.39e-5
$65\,536$ | 32 | $10^{-1}$ | 3.43e-10 | 8.58e-11 | 2.84e-8 | 8.66e-6
$10^{-2}$ | 2.71e-11 | 1.09e-11 | 3.84e-8 | 1.92e-5
$10^{-3}$ | 3.73e-12 | 4.49e-11 | 5.64e-8 | 2.82e-5
64 | $10^{-1}$ | 7.27e-11 | 3.62e-11 | 1.17e-8 | 5.26e-6
$10^{-2}$ | 5.21e-12 | 3.33e-12 | 2.35e-8 | 1.17e-5
$10^{-3}$ | 7.57e-13 | 1.00e-11 | 3.81e-8 | 1.90e-5
128 | $10^{-1}$ | 4.83e-11 | 3.07e-11 | 1.78e-8 | 3.43e-6
$10^{-2}$ | 3.71e-12 | 2.28e-12 | 1.52e-8 | 7.59e-6
$10^{-3}$ | 5.46e-13 | 3.55e-12 | 2.66e-8 | 1.33e-5
Table 4: Advection-Diffusion equation: results for different values of
$\bar{n}$, $\ell$, and $\nu$.

the relative residual norms obtained by fixing the main computational cost of
the two algorithms, namely we perform the same number of parallel-in-time
loops in both schemes. For Ev-Int we use $r=2$ (same results as in Table 3)
and $r=1$. For Algorithm 4.1 ($\alpha=10^{-4}$) we perform only two parallel-
in-time loops by approximating the matrix in (17) by the identity, namely the
FOM-like method is not performed and we set $\mathbf{x}_{m}=\mathbf{b}$ in
line 16 of Algorithm 4.1. The single parallel-in-time loop scenario is
addressed by reporting the accuracy attained by $U_{1}$ in Algorithm 4.1 for
$\alpha=10^{-6}$.

As shown in Table 4, our new ParaDiag algorithm is in general at least as
accurate as Ev-Int in case of two parallel-in-time loops. On the other hand,
whenever a single parallel-in-time loop is performed, our new Algorithm 4.1
($\alpha=10^{-6}$) is often three orders of magnitude more accurate than Ev-
Int. Therefore, in general, our new ParaDiag algorithm achieves better results
in terms of accuracy than Ev-Int whenever a cap on the computational cost of
the adopted solver is imposed.  We would like to mention, however, that a
careful tuning of the Ev-Int parameters may improve the performance of the
solver.

##  7 Conclusions

We presented a new ParaDiag algorithm which fully exploits the circulant-plus-
low-rank structure of the discrete operator stemming from the discretization
of evolution problems (1) by BDFs, one of the most commonly used family of
implicit time integrators. A clever use of the matrix-oriented Sherman-
Morrison-Woodbury formula along with the design of an ad-hoc projection scheme
make our new strategy very successful.

We studied our algorithm for parabolic problems, and also introduced a variant
using $\alpha-$acceleration, based on $\alpha-$circulant matrices. A
comparison with recent ParaDiag techniques from the literature shows that our
new ParaDiag algorithm is competitive and capable of delivering superior
accuracy for comparable cost.

Our methodology can be easily generalized to other discretization schemes as
long as the discrete time operator can be written as a circulant matrix plus a
low-rank correction.  For instance, in [19, Section 4] it is shown that the
all-at-once discretization of (1) by Runge-Kutta methods leads to discrete
problems similar to (2), with a comparable $\Sigma_{1}$. Therefore, our
ParaDiag scheme can also be used to parallelize these Runge-Kutta methods in
time.

## Acknowledgments

This work was supported by the Swiss National Science Foundation, and part has
been carried out at CIRM in Marseille, France, in the context of the Morlet
chair of the first author. The second author is member of the Italian INdAM
Research group GNCS. His work was partially supported by the research project
“Tecniche avanzate per problemi evolutivi: discretizzazione, algebra lineare
numerica, ottimizzazione” (INdAM \- GNCS Project CUP_E55F22000270001).

We would like to thank the anonymous reviewers for their valuable comments and
remarks.

## Appendix

We show now the proof of Theorem 3.1.

######  Proof 7.1.

If $K$ is symmetric, then by its definition also $J_{\ell}$ is symmetric. The
only complex values involved in the definition of $J_{\ell}$ are the
$\pi_{i}$’s and $\gamma_{i}$’s and they appear as either diagonal elements
(the former ones) or scalar multipliers (the latter ones). The eigenvalues
$\pi_{i}$ of the circulant matrix $C_{1}$ are of the form
$\pi_{i}=\omega^{-(\ell-1)(i-1)}$ for all $i=1,\ldots,\ell$, where
$\omega=e^{-\frac{2\pi\imath}{\ell}}\in\mathbb{C}$. Therefore, $|\pi_{i}|\leq
1$ for all $i=1,\ldots,\ell$ and they come in complex conjugate pairs.
Moreover,

$\begin{bmatrix}\gamma_{1}\\\ \vdots\\\ \gamma_{\ell}\\\
\end{bmatrix}={F}^{-T}\mathbf{e}_{\ell}=\frac{1}{\ell}\mathcal{\bar{F}}\mathbf{e}_{\ell}=\frac{1}{\ell}\begin{bmatrix}1\\\
\omega^{-(\ell-1)}\\\ \omega^{-2(\ell-1)}\\\ \vdots\\\
\omega^{-(\ell-1)^{2}}\\\ \end{bmatrix}=\frac{1}{\ell}\begin{bmatrix}1\\\
\pi_{1}\\\ \pi_{2}\\\ \vdots\\\ \pi_{\ell}\\\ \end{bmatrix},$

so that also the $\gamma_{i}$’s come in complex conjugate pairs, and
$|\gamma_{i}|\leq 1/\ell\leq 1$ for any $i=1,\ldots,\ell$, $\ell\geq 1$.

Since $J_{\ell}$ is symmetric, $|\pi_{i}|\leq 1$ for any $i$, and the
$\pi_{i}$’s and $\gamma_{i}$’s all come in complex conjugate pairs, $J_{\ell}$
is Hermitian.

We now show that it is also positive definite. For any $z\in\mathbb{C}^{\ell}$
satisfying $\|z\|=1$, we have
$\text{Re}\left(z^{*}((1-\pi_{i})I_{\bar{n}}+\tau K)z\right)>0$ since $K$ is
SPD, $\tau>0$ and $|\pi_{i}|\leq 1$ for any $i=1,\ldots,\ell$. This means that
also $((1-\pi_{i})I_{\bar{n}}+\tau K)^{-1}$ is positive definite for any $i$.
Furthermore, recalling that $\gamma_{i}=\pi_{i}/\ell$, we have

$\displaystyle\min_{z\in\mathbb{C}^{\ell},\|z\|=1}\text{Re}(z^{*}J_{\ell}z)=$
$\displaystyle
1+\sum_{i=1}^{\ell}\text{Re}(\gamma_{i})\min_{z\in\mathbb{C}^{\ell},\|z\|=1}\text{Re}(z^{*}((1-\pi_{i})I+\tau
K)^{-1}z)$ (26) $\displaystyle=$ $\displaystyle
1+\frac{1}{\ell}\sum_{i=1}^{\ell}\frac{\text{Re}(\pi_{i})}{1-\text{Re}(\pi_{i})+\tau\lambda_{\max}(K)}.$

By construction we know that $\sum_{i=1}^{\ell}\text{Re}(\pi_{i})=0$.
Moreover, if we denote by $\mathfrak{J}$, $\mathfrak{K}$ the index sets such
that $\text{Re}(\pi_{i})\geq 0$ and $\text{Re}(\pi_{i})<0$, respectively,
then, for an even $\ell$, $\mathfrak{J}$ and $\mathfrak{K}$ have the same
cardinality. In particular, $\\#\mathfrak{J}=\\#\mathfrak{K}=\ell/2$ and we
have

$0=\sum_{i=1}^{\ell}\text{Re}(\pi_{i})=\sum_{j\in\mathfrak{J}}\text{Re}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})=\sum_{j\in\mathfrak{J}}(\text{Re}(\pi_{j})-\text{Re}(\pi_{j})),$

namely for a given $j\in\mathfrak{J}$ there exists an index $k\in\mathfrak{K}$
such that $\text{Re}(\pi_{j})=-\text{Re}(\pi_{k})$ if $\ell$ is even.

On the other hand, if $\ell$ is odd, we have $\\#\mathfrak{J}=(\ell-1)/2+1$,
$\\#\mathfrak{K}=(\ell-1)/2$ and we can write

$0=\sum_{i=1}^{\ell}\text{Re}(\pi_{i})=\sum_{j\in\mathfrak{J}}\text{Re}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})=1+\sum_{j\in\mathfrak{J},\pi_{j}\neq
1}\text{Re}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k}),$

which means that

$0<\sum_{j\in\mathfrak{J},\pi_{j}\neq
1}\text{Re}(\pi_{j})=-1-\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k}).$

However, in (7.1) we have a weighted sum of the real parts of the eigenvalues
$\pi_{i}$. In particular, we can write

$\frac{1}{\ell}\sum_{i=1}^{\ell}\frac{\text{Re}(\pi_{i})}{1-\text{Re}(\pi_{i})+\tau\lambda_{\max}(K)}=\sum_{i=1}^{\ell}\text{Re}(\pi_{i})w(\pi_{i}),$

where

$w(\pi_{i}):=\frac{1}{\ell(1-\text{Re}(\pi_{i})+\tau\lambda_{\max}(K))}>0,\quad\text{for
any }i=1,\ldots,\ell.$

We now focus on the scalars $w(\pi_{i})$ and we define

$w(\pi_{i}):=\left\\{\begin{array}[]{l}\check{w}(\pi_{i}),\quad\text{if
}\;i\in\mathfrak{K}{\color[rgb]{0,0,0},}\\\ \hat{w}(\pi_{i}),\quad\text{if
}\;i\in\mathfrak{J}.\end{array}\right.$

A direct computation shows that $\hat{w}(\pi_{j})\geq\check{w}(\pi_{k})$ for
any $j\in\mathfrak{J},k\in\mathfrak{K}$, since

(27)
$\min_{j\in\mathfrak{J}}\hat{w}(\pi_{j})\geq\frac{1}{\ell\tau\lambda_{\max}(K)}\geq\max_{k\in\mathfrak{K}}\check{w}(\pi_{k}).$

For an even $\ell$, it holds that

$\displaystyle\sum_{i=1}^{\ell}\text{Re}(\pi_{i})w(\pi_{i})=$
$\displaystyle\sum_{j\in\mathfrak{J}}\text{Re}(\pi_{j})\hat{w}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})\check{w}(\pi_{k})$
$\displaystyle=$
$\displaystyle\sum_{j\in\mathfrak{J}}\text{Re}(\pi_{j})(\hat{w}(\pi_{j})-\check{w}(\pi_{k_{j}}))\geq
0,$

where the index $k_{j}\in\mathfrak{K}$ is such that
$\text{Re}(\pi_{k_{j}})=-\text{Re}(\pi_{j})$ for a given $j\in\mathfrak{J}$.

On the other hand, if $\ell$ is odd and $j^{*}\in\mathfrak{J}$ is such that
$\pi_{j^{*}}=1$, we can write

$\displaystyle\sum_{i=1}^{\ell}\text{Re}(\pi_{i})w(\pi_{i})=$
$\displaystyle\,\hat{w}(\pi_{j^{*}})+\sum_{j\in\mathfrak{J}\setminus\\{j^{*}\\}}\text{Re}(\pi_{j})\hat{w}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})\check{w}(\pi_{k})$
$\displaystyle\geq$
$\displaystyle\,\hat{w}(\pi_{j^{*}})+\frac{1}{\ell\tau\lambda_{\max}(K)}\sum_{j\in\mathfrak{J}\setminus\\{j^{*}\\}}\text{Re}(\pi_{j})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})\check{w}(\pi_{k})$
$\displaystyle=$
$\displaystyle\,\hat{w}(\pi_{j^{*}})-\frac{1}{\ell\tau\lambda_{\max}(K)}-\frac{1}{\ell\tau\lambda_{\max}(K)}\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})\check{w}(\pi_{k})$
$\displaystyle=$
$\displaystyle\,\hat{w}(\pi_{j^{*}})-\frac{1}{\ell\tau\lambda_{\max}(K)}+\sum_{k\in\mathfrak{K}}\text{Re}(\pi_{k})\left(\check{w}(\pi_{k})-\frac{1}{\ell\tau\lambda_{\max}(K)}\right)\geq
0.$

The nonegativity of the quantity above follows from (27). Indeed,
$\hat{w}(\pi_{j^{*}})\geq\frac{1}{\ell\tau\lambda_{\max}(K)}$ whereas
$\check{w}(\pi_{k})-\frac{1}{\ell\tau\lambda_{\max}(K)}\leq 0$ for any
$k\in\mathfrak{K}$. On the other hand, $\text{Re}(\pi_{k})<0$ for any
$k\in\mathfrak{K}$ so that
$\text{Re}(\pi_{k})\left(\check{w}(\pi_{k})-\frac{1}{\ell\tau\lambda_{\max}(K)}\right)\geq
0$ for any $k\in\mathfrak{K}$.

Therefore, for any $\ell$, we have

$1+\frac{1}{\ell}\sum_{i=1}^{\ell}\frac{\text{Re}(\pi_{i})}{1-\text{Re}(\pi_{i})+\tau\lambda_{\max}(K)}>0,$

and this shows the positive definiteness of $J_{\ell}$.

To conclude, we show the upper bound (14). For any $\ell$, we have

$\displaystyle\kappa(J_{\ell})=$
$\displaystyle\,\frac{\max_{i}|\lambda_{i}(J_{\ell})|}{\min_{i}|\lambda_{i}(J_{\ell})|}=\frac{\max_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}J_{\ell}z|}{\min_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}J_{\ell}z|}$
$\displaystyle=$
$\displaystyle\,\frac{1+\max_{z\in\mathbb{C}^{\ell},\|z\|=1}|\sum_{i=1}^{\ell}\gamma_{i}z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}z|}{1+\min_{z\in\mathbb{C}^{\ell},\|z\|=1}|\sum_{i=1}^{\ell}\gamma_{i}z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}z|}.$

Recalling that $|\gamma_{i}|\leq 1/\ell$ and $0\leq|1-\pi_{i}|\leq 2$ for any
$i$ we get

$\displaystyle\kappa(J_{\ell})\leq$
$\displaystyle\,\frac{1+\max_{i}\max_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}z|}{1+\ell\min_{i}\min_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)^{-1}z|}$ $\displaystyle=$
$\displaystyle\,\frac{1+\max_{i}\frac{1}{\min_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)z|}}{1+\ell\min_{i}\frac{1}{\max_{z\in\mathbb{C}^{\ell},\|z\|=1}|z^{*}((1-\pi_{i})I_{\bar{n}}+\tau
K)z|}}$ $\displaystyle\leq$
$\displaystyle\,\frac{1+\frac{1}{\tau\lambda_{\min}(K)}}{1+\frac{\ell}{2+\tau\lambda_{\max}(K)}}=\frac{\tau\lambda_{\min}(K)+1}{2+\tau\lambda_{\max}(K)+\ell}\cdot\frac{2+\tau\lambda_{\max}(K)}{\tau\lambda_{\min}(K)}$
$\displaystyle=$
$\displaystyle\,\frac{\tau\lambda_{\min}(K)}{\tau\lambda_{\max}(K)}\frac{1+\frac{1}{\tau\lambda_{\min}(K)}}{1+\frac{2+\ell}{\tau\lambda_{\max}(K)}}\cdot\frac{\tau\lambda_{\max}(K)}{\tau\lambda_{\min}(K)}\left(1+\frac{2}{\tau\lambda_{\max}(K)}\right)$
$\displaystyle\leq$ $\displaystyle\,1+\frac{1}{\tau\lambda_{\min}(K)}.$

## References

  * [1] U. M. Ascher and L. R. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1998.
  * [2] F. Danieli and A. J Wathen, All-at-once solution of linear wave equations, Numerical Linear Algebra with Applications, 28 (2021), p. e2386.
  * [3] A. Frommer and U. Glässner, Restarted GMRES for shifted linear systems, SIAM J. Sci. Comput., 19 (1998), pp. 15–26.  Special issue on iterative methods (Copper Mountain, CO, 1996).
  * [4] A. Frommer, K. Kahl, Th. Lippert, and H. Rittich, 2-norm error bounds and estimates for Lanczos approximations to linear systems and rational matrix functions, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 1046–1065.
  * [5] A. Frommer, K. Kahl, and M. Tsolakis, Matrix functions via linear systems built from continued fractions, arXiv preprint arXiv:2109.03527, (2021).
  * [6] A. Frommer and V. Simoncini, Stopping criteria for rational matrix functions of Hermitian and symmetric matrices, SIAM J. Sci. Comput., 30 (2008), pp. 1387–1412.
  * [7] M. J. Gander, 50 years of time parallel time integration, in Multiple Shooting and Time Domain Decomposition Methods, T. Carraro, M. Geiger, S. Körkel, and R. Rannacher, eds., Springer, 2015, pp. 69–114.
  * [8] M. J. Gander and E. Hairer, Nonlinear convergence analysis for the Parareal algorithm, in Domain Decomposition Methods in Science and Engineering, Lecture Notes in Computational Science and Engineering, O. B. Widlund and D. E. Keyes, eds., vol. 60, Springer, 2008, pp. 45–56.
  * [9] M. J. Gander and L. Halpern, Time parallelization for nonlinear problems based on diagonalization, in Domain decomposition methods in science and engineering XXIII, Springer, 2017, pp. 163–170.
  * [10] M. J. Gander, L. Halpern, J. Rannou, and J. Ryan, A direct time parallel solver by diagonalization for the wave equation, SIAM J. Sci. Comput., 41 (2019), pp. A220–A245.
  * [11] M. J. Gander, L. Halpern, J. Ryan, and T. T. B. Tran, A direct solver for time parallelization, in Domain decomposition methods in science and engineering XXII, vol. 104 of Lect. Notes Comput. Sci. Eng., Springer, Cham, 2016, pp. 491–499.
  * [12] M. J. Gander, Y.-L. Jiang, B. Song, and H. Zhang, Analysis of two parareal algorithms for time-periodic problems, SIAM Journal on Scientific Computing, 35 (2013), pp. A2393–A2415.
  * [13] M. J. Gander, J. Liu, S.-L. Wu, X. Yue, and T. Zhou, Paradiag: Parallel-in-time algorithms based on the diagonalization technique, arXiv preprint arXiv:2005.09158, (2020).
  * [14] M. J. Gander and S. Vandewalle, Analysis of the Parareal time-parallel time-integration method, SIAM j. Sci. Comput, 29 (2007), pp. 556–578.
  * [15] M. J. Gander and S.-L. Wu, Convergence analysis of a periodic-like waveform relaxation method for initial-value problems via the diagonalization technique, Numerische Mathematik, 143 (2019), pp. 489–527.
  * [16] , A diagonalization-based parareal algorithm for dissipative and wave propagation problems, SIAM Journal on Numerical Analysis, 58 (2020), pp. 2981–3009.
  * [17] M. H. Gutknecht, Krylov subspace algorithms for systems with multiple right hand sides: an introduction, in Modern mathematical models, methods and algorithms for real world systems, A. H. Siddiqi, I. Duff, and O. Christensen, eds., Anshan Ltd, 2007.  Available at http://www.sam.math.ethz.ch/~mhg/pub/delhipap.pdf.
  * [18] T. Hawkins, Weierstrass and the theory of matrices, Archive for history of exact sciences, 17 (1977), pp. 119–163.
  * [19] D. Kressner, S. Massei, and J. Zhu, Improved parallel-in-time integration via low-rank updates and interpolation, ArXiv Preprint: 2204.03073, (2022).
  * [20] M. Kubínová and K. M. Soodhalter, Admissible and Attainable Convergence Behavior of Block Arnoldi and GMRES, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 464–486.
  * [21] J. Liesen and Z. Strakoš, Krylov subspace methods, Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2013\.
  * [22] J.-L. Lions, Y. Maday, and G. Turinici, A ”Parareal” in time discretization of PDE’s, Compte Rendu de l’Academie de Sciences, 332 (2001), pp. 661–668.
  * [23] J. Liu and S.-L. Wu, A fast block $\alpha$-circulant preconditoner for all-at-once systems from wave equations, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 1912–1943.
  * [24] Y. Maday and E. M. Rønquist, Parallelization in time through tensor-product space-time solvers, C. R. Math. Acad. Sci. Paris, 346 (2008), pp. 113–118.
  * [25] E. McDonald, J. Pestana, and A. Wathen, Preconditioning and iterative solution of all-at-once systems for evolutionary partial differential equations, SIAM J. Sci. Comput., 40 (2018), pp. A1012–A1033.
  * [26] S. Noschese and L. Reichel, Generalized circulant Strang-type preconditioners, Numer. Linear Algebra Appl., 19 (2012), pp. 3–17.
  * [27] D. Palitta, Matrix equation techniques for certain evolutionary partial differential equations, J. Sci. Comput., 87 (2021).
  * [28] D. Palitta and P. Kürschner, On the convergence of Krylov methods with low-rank truncations, Numer. Algorithms, 88 (2021), pp. 1383–1417.
  * [29] Y. Saad, Iterative methods for sparse linear systems, SIAM, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2nd ed., 2003\.
  * [30] V. Simoncini, Restarted full orthogonalization method for shifted linear systems, BIT Numerical Mathematics, 43 (2003), pp. 459–466.
  * [31] S.-L. Wu, Toward parallel coarse grid correction for the parareal algorithm, SIAM Journal on Scientific Computing, 40 (2018), pp. A1446–A1472.
  * [32] S.-L. Wu and J. Liu, A parallel-in-time block-circulant preconditioner for optimal control of wave equations, SIAM Journal on Scientific Computing, 42 (2020), pp. A1510–A1540.
  * [33] S.-L. Wu and T. Zhou, Acceleration of the two-level MGRIT algorithm via the diagonalization technique, SIAM Journal on Scientific Computing, 41 (2019), pp. A3421–A3448.
  * [34] , Parallel implementation for the two-stage SDIRK methods via diagonalization, Journal of Computational Physics, 428 (2021).