# A parallel in time algorithm based ParaExp for optimal control problems

Felix Kwok, Djahou N Tognon

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

## To cite this version:

Felix Kwok, Djahou N Tognon. A parallel in time algorithm based ParaExp for optimal control problems. 2024. hal-04612018

HAL Id: hal-04612018

https://inria.hal.science/hal-04612018

Preprint submitted on 14 Jun 2024

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
<figure><img src="image_1.png" /></figure>

HAL
open science

# A Parallel in Time Algorithm Based on ParaExp for Optimal Control Problems

Felix Kwok, Djahou N Tognon

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

## To cite this version:

Felix Kwok, Djahou N Tognon. A Parallel in Time Algorithm Based on ParaExp for Optimal Control Problems. 2024 IEEE 63rd Conference on Decision and Control (CDC 2024), Dec 2024, Milan, Italy. hal-04612018v4

HAL Id: hal-04612018

https://inria.hal.science/hal-04612018v4

Submitted on 24 Sep 2025

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A PARALLEL IN TIME ALGORITHM BASED PARAEXP FOR
OPTIMAL CONTROL PROBLEMS

Felix Kwok

Département de mathématiques et
de statistique, Université Laval
Canada

felix.kwok@mat.ulaval.ca

Djahou N. Tognon

INRIA Paris, ANGE Project-Team,
Sorbonne Université, (LJLL),
France

djahou-norbert.tognon@inria.fr

ABSTRACT

We propose a new parallel-in-time algorithm for solving optimal control problems constrained by
partial differential equations. Our approach, which is based on a deeper understanding of ParaExp,
considers an overlapping time-domain decomposition in which we combine the solution of homoge-
neous problems using exponential propagation with the local solutions of inhomogeneous problems.
The algorithm yields a linear system whose matrix-vector product can be fully performed in parallel.
We then propose a preconditioner to speed up the convergence of GMRES in the special cases of
the heat and wave equations. Numerical experiments are provided to illustrate the efficiency of our
preconditioners.

Keywords ParaExp Algorithm · Optimality systems · Preconditioning

# 1 Introduction

The efficient solution of optimal control problems constrained by time-dependent partial differential equations (PDEs) has received much interest in recent decades, given its importance in applications and the increasing availability of powerful modern computing clusters. The main computational challenges lie in the large amounts of data that must be processed, the intensive computation required to simulate the underlying PDEs, and the optimization procedure needed to determine the best control function. Parallel numerical algorithms, specifically those of the domain decomposition type, are well suited for tackling these huge problems by subdividing them into manageable chunks that can be handled simultaneously by many processors. Parallelism also has the obvious advantage of reducing the wall-clock time required for finding the optimal solution, which is especially important for real-time applications.

There exists a rich literature on spatial domain decomposition methods for elliptic and initial value problems, see
[17, 7, 12, 13] and the references therein. There is also a growing literature in parallel-in-time methods, in multiple
processors solve the initial value problem on different parts of the time interval simultaneously; see [5, 16] for a survey
of these methods for initial value problems. Time parallelization can also be applied to optimal control problems, either
by decomposing the optimization problem directly [10, 14, 1, 8], or by using a parallel-in-time method to solve the
forward and adjoint problems required by gradient computations [19].

In this paper, we focus on the efficient solution of linear-quadratic optimal control problems when a cheap exponential integrator is available, i.e., when one is able to evaluate $\exp(t\mathcal{L})$ quickly for a linear operator $\mathcal{L}$. For initial value problems, the authors of [6] have used such integrators to construct a parallel-in-time method known as ParaExp. In practice, many different methods can be used to evaluate $\exp(t\mathcal{L})$, and the choice depends on the particular properties of the operator $\mathcal{L}$. In this paper, we show how similar ideas can be used to solve optimal control problems; we will focus on the particular cases of the heat and wave equations, where specific approximations are used to construct the preconditioner for solving for the optimal adjoint state.

The rest of the paper is organized in follows. We present our model optimal control problem in Section 2. The parallel-in-time algorithm is defined in Section 3, where we show how computations from different parts of the time interval can be performed in parallel. The algorithm yields a linear system whose matrix-vector product can be performed in
# A Parallel in Time Algorithm Based on ParaExp for Optimal Control Problems

Felix Kwok and Djahou N. Tognon

**Abstract**— We propose a new parallel-in-time algorithm for solving optimal control problems constrained by discretized partial differential equations. Our approach, which is based on a deeper understanding of ParaExp, considers an overlapping time-domain decomposition in which we combine the solution of homogeneous problems using exponential propagation with the local solutions of inhomogeneous problems. The algorithm yields a linear system whose matrix-vector product can be fully performed in parallel. We then propose a preconditioner to speed up the convergence of GMRES in the special cases of the heat and wave equations. Numerical experiments are provided to illustrate the efficiency of our preconditioners.

## I. INTRODUCTION

The efficient solution of optimal control problems constrained by time-dependent partial differential equations (PDEs) has received much interest in recent decades, given its importance in applications and the increasing availability of powerful modern computing clusters. The main computational challenges lie in the large amounts of data that must be processed, the intensive computation required to simulate the underlying PDEs, and the optimization procedure needed to determine the best control function. Parallel numerical algorithms, specifically those of the domain decomposition type, are well suited for tackling these huge problems by subdividing them into manageable chunks that can be handled simultaneously by many processors. Parallelism also has the obvious advantage of reducing the wall-clock time required to find the optimal solution, which is especially important for real-time applications.

There exists a rich literature on spatial domain decomposition methods for elliptic and initial value problems, see [18], [7], [12], [13] and the references therein. There is also a growing literature on parallel-in-time methods, in which multiple processors solve the initial value problem on different parts of the time interval simultaneously; see [5], [16] for a survey of these methods. Time parallelization can also be applied to optimal control problems, either by decomposing the optimization problem directly, see [10], [14], [1], [8], or by using a parallel-in-time method to solve the forward and adjoint problems required by gradient computations [20].

In this paper, we focus on the efficient solution of linear-quadratic optimal control problems when a cheap exponential integrator is available, i.e., when one is able to evaluate

$\exp(t\mathcal{L})$ quickly for a matrix $\mathcal{L}$ arising from spatial discretization. For initial value problems, the authors of [6] have used such integrators to construct a parallel-in-time method known as ParaExp. In practice, many different methods can be used to evaluate $\exp(t\mathcal{L})$, and the choice depends on the particular properties of the operator $\mathcal{L}$. In this paper, we show how similar ideas can be used to solve optimal control problems; we will focus on the particular cases of the heat and wave equations, where specific approximations are used to construct the preconditioner for solving for the optimal adjoint state. Numerical experiments are provided to illustrate the behavior of these preconditioners.

## II. LINEAR OPTIMAL CONTROL PROBLEM

We consider the optimal control problem constrained by a system of ordinary differential equations (ODEs)

$$ \begin{array}{ll} \min_{v} & \frac{1}{2} \|y(T) - y_{tg}\|^2 + \frac{\alpha}{2} \int_{0}^{T} \|v(t)\|^2 dt, \\ \text{subject to} & \dot{y}(t) = \mathcal{L}y(t) + v(t), \quad y(0) = y_{in}, \quad t \in (0, T), \end{array} $$

where the constraint arises from a semi-discretization in space of a time-dependent PDE. Here, $\alpha$ is a regularization parameter and $y_{in}, y_{tg}$ are the initial and target states respectively. The state function $y$ and the control $v$ are mappings from $(0, T) \to \mathbb{R}^r$, $r$ a positive integer. The operator $\mathcal{L} \in \mathbb{R}^{r \times r}$ is independent of the time variable. Using the Lagrange multiplier approach, one can characterize the optimal trajectory $y(t)$ via the forward and adjoint problems

$$ \begin{cases} \dot{y} = \mathcal{L}y + v & \text{on } (0, T), \\ y(0) = y_{in}, \end{cases} \quad \begin{cases} \dot{\lambda} = -\mathcal{L}^T \lambda & \text{on } (0, T) \\ \lambda(T) = y(T) - y_{tg}, \end{cases} $$

where the control $v$ and adjoint state $\lambda$ are related by the algebraic equation $\lambda(t) = \alpha v(t)$ for all $t \in (0, T)$. Eliminating $v$, the above system can also be written as

$$ \dot{y} = \mathcal{L}y - \frac{1}{\alpha}\lambda, \quad \dot{\lambda} = -\mathcal{L}^T \lambda \quad \text{on } (0, T) \quad (1) $$

with the initial and final condition $y(0) = y_{in}$ and $\lambda(T) = y(T) - y_{tg}$ respectively.

## III. PARALLEL-IN-TIME ALGORITHM

For a positive integer $L$, we consider the non-overlapping sub-intervals $(T_{l-1}, T_l), l = 1, \dots, L$ of $(0, T)$ with $T_l = l\Delta T$ and $\Delta T = T/L$. We define two sets of intermediate states $(Y_l)_{l=1, \dots, L}$ and $(\Lambda_l)_{l=1, \dots, L}$, corresponding to the state $y$ and the adjoint $\lambda$ at times $T_1, \dots, T_L$. We now introduce a new parallelization idea for solving (1) inspired by the ParaExp algorithm [6], namely to decompose the ODEs into

<footnote>F. Kwok, Université Laval, Département de mathématiques et de statistique, Canada felix.kwok@mat.ulaval.ca</footnote>

<footnote>D. N. Tognon, INRIA Paris : ANGE Project-Team, and Sorbonne Université: Laboratoire Jacques-Louis Lions (LJLL), France djahou-norbert.tognon@inria.fr</footnote>
parallel, but the system is generally ill-conditioned; a preconditioner is therefore needed to solve it efficiently using iterative methods. Section 4 is devoted to the construction of the preconditioner for two model problems, namely the heat and wave equations. Numerical experiments are provided to illustrate the behaviour of these preconditioners. Finally, we summarize and outline some future work in Section 5.

## 2 Linear optimal control problem

We consider the optimal control problem given by

$$ 
\begin{array}{ll}
\displaystyle \min_{\nu} & \frac{1}{2} \|y(T) - y_{tg}\|^2 + \frac{\alpha}{2} \int_0^T \|\nu(t)\|^2 dt, \\
\text{subject to} & \dot{y}(t) = \mathcal{L}y(t) + \nu(t), \quad y(0) = y_{in}, \quad t \in (0, T),
\end{array}
$$

where $\alpha$ is a regularization parameter and $y_{in}, y_{tg}$ are the initial and target states respectively. The state function $y$ and the control $\nu$ are defined from $(0, T)$ to $\mathbb{R}^r, r \in \mathbb{N}^*$. The operator $\mathcal{L} \in \mathbb{R}^{r \times r}$ is independent of the time variable. Using the Lagrange multiplier approach (see e.g. the historical review [3]), one can derive the forward and adjoint problems to be

$$ 
\begin{cases}
\dot{y} = \mathcal{L}y + \nu & \text{on } (0, T), \\
y(0) = y_{in},
\end{cases}
\quad
\begin{cases}
\dot{\lambda} = -\mathcal{L}^T \lambda & \text{on } (0, T) \\
\lambda(T) = y(T) - y_{tg},
\end{cases}
$$

where the control $\nu$ and adjoint state $\lambda$ are related by the algebraic equation $\lambda(t) = \alpha\nu(t)$ for all $t \in (0, T)$. Eliminating $\nu$, the above system can also be written as

$$ \dot{y} = \mathcal{L}y - \frac{1}{\alpha}\lambda, \quad \dot{\lambda} = -\mathcal{L}^T \lambda \quad \text{on } (0, T) \qquad (1) $$

with the initial and final condition $y(0) = y_{in}$ and $\lambda(T) = y(T) - y_{tg}$ respectively.

## 3 Parallel-in-time algorithm

For a given $L \in \mathbb{N}^*$, we consider the non-overlapping sub-intervals $(T_{\ell-1}, T_{\ell})$, $\ell = 1, \dots, L$ of $(0, T)$ with $T_{\ell} = \ell\Delta T$ and $\Delta T = T/L$. We define two sets of intermediate states $(Y_{\ell})_{\ell=1,\dots,L}$ and $(\Lambda_{\ell})_{\ell=1,\dots,L}$, corresponding to approximations of the state $y$ and the adjoint state $\lambda$ at times $T_1, \dots, T_L$. We now introduce a new parallelization idea for solving (1) inspired by the ParaExp algorithm [6], namely to decompose the ODEs into homogeneous and inhomogeneous parts. First, we consider the homogeneous sub-problems on $\lambda$, which require the backward propagation of $\lambda$ over each $(T_{\ell-1}, T_L)$ according to

$$ \dot{\lambda}_{\ell}(t) = -\mathcal{L}^T \lambda_{\ell}(t), \quad \lambda_{\ell}(T_L) = \Lambda_L, \quad \text{on } (T_{\ell-1}, T_L), \qquad (2) $$

Next, we define the inhomogeneous part of the problem on $y$, given by

$$ \dot{w}_{\ell}(t) = \mathcal{L}w_{\ell}(t) - \frac{1}{\alpha}\lambda_{\ell}(t), \quad w_{\ell}(T_{\ell-1}) = 0, \quad \text{on } (T_{\ell-1}, T_{\ell}), \qquad (3) $$

for $\ell = 1, \dots, L$. Note that the $w_\ell$ are solved with no initial conditions; they are instead carried by the homogeneous sub-problems $u_\ell$ defined below, which perform the forward propagation of $y$ over $(T_{\ell-1}, T_\ell)$:

$$ 
\begin{aligned}
\dot{u}_{\ell}(t) &= \mathcal{L}u_{\ell}(t), \quad u_{\ell}(T_0) = y_{in}, \quad \text{on } (T_0, T_L) \\
\dot{u}_{\ell}(t) &= \mathcal{L}u_{\ell}(t), \quad u_{\ell}(T_{\ell}) = w_{\ell-1}(T_{\ell-1}), \quad \text{on } (T_{\ell-1}, T_L),
\end{aligned}
\qquad (4) $$

for $\ell = 2, \dots, L$. Finally, $Y_\ell \approx y_\ell(T_\ell)$ is directly obtained by superposition principle

$$ Y_{\ell} = w_{\ell}(T_{\ell}) + \sum_{j=1}^{\ell} u_{j}(T_{\ell}) \quad \text{on } (T_{\ell-1}, T_{\ell}), \quad \ell = 1, \dots, L. \qquad (5) $$

Thus, for a given $\Lambda_L$, we have the following data dependency: the $\lambda_\ell$ is independent of everything, $w_\ell$ depends only on $\lambda_\ell$ as a source term, and $u_{\ell+1}$ depends on $w_\ell$ for initial conditions. Moreover, the computation for different indices $\ell$ are independent. Thus, for $L$ available processors, it is natural to assign $\lambda_\ell$, $w_\ell$ and $u_{\ell+1}$ to the same processor for $\ell = 1, \dots, L-1$ and $\lambda_L, w_L$ and $u_1$ to the Lth processor, as shown in Figure 1. As we only need $y_\ell$ at the interfaces to recreate the optimal trajectory (5), the only communication needed between processors are the values of $u$ and $w$ at the interfaces $T_\ell$.
homogeneous and inhomogeneous parts. First, we consider the homogeneous sub-problems on $\lambda$, which require the backward propagation of $\lambda$ over each $(T_{\ell-1}, T_L)$:

$$ \dot{\lambda}_{\ell}(t) = -\mathcal{L}^{T} \lambda_{\ell}(t), \quad \lambda_{\ell}(T_L) = \Lambda_L, \quad \text{on } (T_{\ell-1}, T_L), \quad (2) $$

Next, we define the inhomogeneous part of the problem on $y$, given by

$$ \dot{w}_{\ell}(t) = \mathcal{L} w_{\ell}(t) - \frac{1}{\alpha} \lambda_{\ell}(t), \quad w_{\ell}(T_{\ell-1}) = 0, \quad \text{on } (T_{\ell-1}, T_{\ell}), \quad (3) $$

for $\ell = 1, \dots, L$. Note that the $w_\ell$ are solved with no initial conditions; they are instead carried by the homogeneous sub-problems $u_\ell$ defined below, which perform the forward propagation of $y$ over $(T_{\ell-1}, T_\ell)$:

$$ \dot{u}_{1}(t) = \mathcal{L} u_{1}(t), \quad u_{1}(T_{0}) = y_{in}, \quad \text{on } (T_{0}, T_{L}) \quad (4) $$

$$ \dot{u}_{\ell}(t) = \mathcal{L} u_{\ell}(t), \quad u_{\ell}(T_{\ell-1}) = w_{\ell-1}(T_{\ell-1}), \quad \text{on } (T_{\ell-1}, T_{\ell}), $$

for $\ell = 2, \dots, L$. Finally, $y(T_\ell)$ is obtained by superposition:

$$ y(T_\ell) = w_\ell(T_\ell) + \sum_{j=1}^{\ell} u_j(T_\ell), \quad \ell = 1, \dots, L. \quad (5) $$

Figure 1 summarizes the dependency between $\lambda_\ell$, $w_\ell$ and $u_\ell$. Thus, for $L$ available processors, it is natural to assign $\lambda_\ell$, $w_\ell$ and $u_{\ell+1}$ to the same processor for $\ell = 1, \dots, L-1$ and $\lambda_L$, $w_L$ and $u_1$ to the $L^{th}$ processor. Note that to calculate $y(T_\ell)$, the processors only need to exchange the values of $u$ and $w$ at the interfaces $T_\ell$, rather than whole trajectories.

Once the $Y_\ell = y(T_\ell)$ and $\Lambda_{\ell+1} = \lambda(T_{\ell+1})$ are found, the trajectories $y(t)$ and $\lambda(t)$ on $t \in [T_\ell, T_{\ell+1}]$ can easily be constructed in parallel using (1).

Let $\mathcal{P}_\ell$ and $\mathcal{Q}_\ell$ be the solution operators for $u_\ell$ and $\lambda_\ell$ defined by the homogeneous sub-problems (4) and (2). (Even though $\mathcal{L}$ is identical on each interval, we write the index $\ell$ explicitly to indicate the sub-interval on which the propagation is performed.) In other words, the solution of (2) is given by

$$ \lambda_{\ell}(t) = \mathcal{Q}_{\ell}(t) \cdot \Lambda_{L}, \quad \text{on } (T_{\ell-1}, T_{L}), \quad (6) $$

and the one of (4) is given by

$$ u_1(t) = \mathcal{P}_1(t) \cdot y_{in}, \quad \text{on } (T_0, T_L) \quad (7) $$

$$ u_{\ell}(t) = \mathcal{P}_{\ell}(t) \cdot w_{\ell-1}(T_{\ell-1}), \quad \text{on } (T_{\ell-1}, T_{\ell}). $$

Thus, substituting (6) into (3) on $(T_{\ell-1}, T_\ell)$, we get

$$ \dot{w}_{\ell}(t) = \mathcal{L} w_{\ell}(t) - \frac{1}{\alpha} \mathcal{Q}_{\ell}(t) \cdot \Lambda_{L}, \quad w_{\ell}(T_{\ell-1}) = 0. \quad (8) $$

Since $w_\ell(T_{\ell-1}) = 0$, the above implies

$$ w_{\ell}(t) = -\frac{1}{\alpha} \left[ \int_{T_{\ell-1}}^{t} e^{(t-s)\mathcal{L}} \mathcal{Q}_{\ell}(s) ds \right] \cdot \Lambda_{L}. $$

Let $\mathcal{R}_\ell$ denote the solution operator of $w_\ell$ on $(T_{\ell-1}, T_\ell)$, such that

$$ w_{\ell}(t) = -\frac{1}{\alpha} \mathcal{R}_{\ell}(t) \cdot \Lambda_{L}. \quad (9) $$

Therefore, substituting (9) into (7), we obtain

$$ u_1(t) = \mathcal{P}_1(t) \cdot y_{in} \quad \text{on } (T_0, T_L), \quad (10) $$

$$ u_{\ell}(t) = -\frac{1}{\alpha} \mathcal{P}_{\ell}(t) \cdot \mathcal{R}_{\ell-1}(T_{\ell-1}) \cdot \Lambda_{L}, \quad \text{on } (T_{\ell-1}, T_{\ell}), $$

<figure><img src="image_4.png" /><figcaption>Fig. 1. Data dependency and processor assignment for the parallel computation of $\mathcal{M} \cdot \Lambda_L$. Note that $u_1$ is independent of $\Lambda_L$ and is calculated only once when forming the right-hand side $\mathbf{b}$ in (13).</figcaption></figure>

and inserting (9) and (10) into (5) yields

$$ y(T_\ell) = \mathcal{P}_1(T_\ell) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_\ell(T_\ell) \cdot \Lambda_L \\ - \frac{1}{\alpha} \sum_{j=2}^{\ell} \mathcal{P}_j(T_\ell) \cdot \mathcal{R}_{j-1}(T_{j-1}) \cdot \Lambda_L, \quad \text{on } (T_{\ell-1}, T_\ell). $$

Hence, the optimality system (1) in terms of $Y_\ell = y(T_\ell)$ and $\Lambda_\ell = \lambda(T_\ell)$ becomes

$$ Y_{\ell} = \mathcal{P}_{1}(T_{\ell}) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_{\ell}(T_{\ell}) \cdot \Lambda_{L} \\ - \frac{1}{\alpha} \sum_{j=2}^{\ell} \mathcal{P}_{j}(T_{\ell}) \cdot \mathcal{R}_{j-1}(T_{j-1}) \cdot \Lambda_{L}, \quad \ell = 1, \dots, L, \quad (11) $$

and $\Lambda_\ell = \mathcal{Q}_\ell(T_\ell) \cdot \Lambda_L, \ell = 1, \dots, L-1$. Now that we expressed the $Y_\ell$ and $\Lambda_\ell$ in terms of a single unknown $\Lambda_L$, we only need to impose one last equation to close the system: we use the final condition in (1), given by

$$ \Lambda_L - Y_L + y_{tg} = 0. \quad (12) $$

Substituting (11) for $\ell = L$ into the above, we obtain

$$ \mathcal{M} \cdot \Lambda_L + \mathbf{b} = 0, \quad (13) $$

where $\mathbf{b} := y_{tg} - \mathcal{P}_1(T_L) \cdot y_{in}$ and

$$ \mathcal{M} := I + \frac{1}{\alpha} \mathcal{R}_L(T_L) + \frac{1}{\alpha} \sum_{j=2}^{L} \mathcal{P}_j(T_L) \cdot \mathcal{R}_{j-1}(T_{j-1}). $$

with $I$ being the identity matrix. Just like for (5), the matrix-vector product $\mathcal{M} \cdot \Lambda_L$ can be computed fully in parallel: each term of the form $\mathcal{P}_j(T_L) \cdot \mathcal{R}_{j-1}(T_{j-1}) \cdot \Lambda_L$ in the sum above can be performed on the $(j-1)^{th}$ processor for $j = 2, \dots, L$, while the term $\mathcal{R}_L(T_L) \cdot \Lambda_L$ can be handled by the $L^{th}$ processor.

In practice, the solution operators $\mathcal{P}_\ell$ and $\mathcal{R}_\ell$ are matrix exponentials that need to be approximated numerically. Computing this approximation efficiently is a non-trivial task: for small and medium-sized dense matrices ($r \le 1000$), various methods can be found in [15], [11] and [9]. For large sparse matrices, however, the task frequently requires an understanding of the spectral properties of the underlying operator $\mathcal{L}$. One possibility proposed by [6] is projection-based methods. There, the authors claim that “we clearly
Let $\mathcal{P}_\ell$ and $\mathcal{Q}_\ell$ be the exponential propagators for $u_\ell$ and $\lambda_\ell$, as defined by the homogeneous sub-problems (4) and (2). (Even though $\mathcal{L}$ is identical on each interval, we write the index $\ell$ explicitly to indicate the sub-interval on which the propagation is performed.) Then the solution of (2) becomes

$$ \lambda_{\ell}(t) = \mathcal{Q}_{\ell}(t) \cdot \Lambda_{L}, \quad \text{on } (T_{\ell-1}, T_{\ell}). \qquad (6) $$

Substituting (6) into (3) on $(T_{\ell-1}, T_\ell)$ gives

$$ \dot{w}_{\ell}(t) = \mathcal{L}w_{\ell}(t) - \frac{1}{\alpha}\mathcal{Q}_{\ell}(t) \cdot \Lambda_{L}, \quad w_{\ell}(T_{\ell-1}) = 0. \qquad (7) $$

<figure><img src="image_3.png" /><figcaption>Figure 1: Parallel distribution on the processors</figcaption></figure>

Moreover, because $w_\ell(T_{\ell-1}) = 0$ in (7), a continuous solution of (7) is given by

$$ w_{\ell}(t) = -\frac{1}{\alpha} \left[ \int_{T_{\ell-1}}^{t} e^{(t-s)\mathcal{L}} \mathcal{Q}_{\ell}(s) ds \right] \cdot \Lambda_{L}. $$

Let $\mathcal{R}_\ell$ be the solution operator of $w_\ell$ on $(T_{\ell-1}, T_\ell)$, such that

$$ w_{\ell}(t) = -\frac{1}{\alpha} \mathcal{R}_{\ell}(t) \cdot \Lambda_{L}. \qquad (8) $$

Therefore, using (8), (4) become

$$ \begin{aligned} u_1(t) &= \mathcal{P}_1(t) \cdot y_{in} \quad \text{on } (T_0, T_L), \\ u_\ell(t) &= -\frac{1}{\alpha} \mathcal{P}_\ell(t) \cdot \mathcal{R}_{\ell-1}(T_{\ell-1}) \cdot \Lambda_L, \quad \text{on } (T_{\ell-1}, T_L). \end{aligned} \qquad (9) $$

Now, replacing (8) and (9) into (5) leads to

$$ y_{\ell}(t) = \mathcal{P}_{1}(t) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_{\ell}(T_{\ell}) \cdot \Lambda_{L} - \frac{1}{\alpha} \sum_{j=2}^{\ell} \mathcal{P}_{j}(T_{\ell}) \cdot \mathcal{R}_{j-1}(T_{j-1}) \cdot \Lambda_{L}, \quad \text{on } (T_{\ell-1}, T_{\ell}). $$

Hence, the optimality system (1) in terms of $Y_\ell$ and $\Lambda_\ell$ becomes

$$ Y_{\ell} = \mathcal{P}_{1}(T_{\ell}) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_{\ell}(T_{\ell}) \cdot \Lambda_{L} - \frac{1}{\alpha} \sum_{j=2}^{\ell} \mathcal{P}_{j}(T_{\ell}) \cdot \mathcal{R}_{j-1}(T_{j-1}) \cdot \Lambda_{L}, \quad \ell = 1, \dots, L, \qquad (10) $$

with

$$ \Lambda_{\ell} = \mathcal{Q}_{\ell}(T_{\ell}) \cdot \Lambda_{L}, \quad \ell = 1, \dots, L-1. $$

We remark that each $Y_\ell$ depends on the same unknown $\Lambda_L$. Moreover, $Y_L$ and $\Lambda_L$ satisfy the discrete final condition of (1) given by

$$ \Lambda_L - Y_L + y_{tg} = 0. \qquad (11) $$

It follows that, for $\ell = L$, substituting (10) into (11) yields

$$ M \cdot \Lambda_L + \mathbf{b} = 0, \qquad (12) $$

where $\mathbf{b} := y_{tg} - \mathcal{P}_1(T_L) \cdot y_{in}$ and

$$ M := I + \frac{1}{\alpha} \mathcal{R}_L(T_L) + \frac{1}{\alpha} \sum_{j=2}^{L} \mathcal{P}_j(T_L) \cdot \mathcal{R}_{j-1}(T_{j-1}). $$
find that it is beneficial to use the Arnoldi method rather than a time-stepping method for the propagation of the linear homogeneous problems. Note that the difference between the projection and time-stepping methods gets large for high accuracy,” meaning that the Arnoldi method can produce accurate approximations efficiently.

## IV. PRECONDITIONER DESIGN

It remains to solve the linear system (13) using an iterative method, e.g. GMRES. Since the matrix $\mathcal{M}$ is generally ill-conditioned, a preconditioner is needed, which we now construct for the semi-discretized heat and wave equations.

### A. Heat equation

In this section, we propose a preconditioner $\tilde{\mathcal{M}}^{-1}$ for the heat equation to speed up the convergence of (13) in GMRES. Our approach is inspired by the work in the Master's thesis [19], where the author derived a preconditioner based on the behavior of $\mathcal{L}$ at low and high frequencies. Let us consider the heat equation given by

$$ \dot{y} = \Delta y + \mathbf{v}, \text{ on } \Omega \times (0, T), \quad (14) $$

with $y(x, 0) = y_{in}(x)$ on $\Omega$. We consider the 1D case where $\Omega = (0, 1)$ with Dirichlet boundary conditions $y(0, t) = y(1, t) = 0, t \in (0, T)$ with $T \ge 1$. Then, a semi-discretization in space of (14) using second-order centered finite differences leads to the following ODE system

$$ \dot{y}(t) = \mathcal{L}y(t) + \mathbf{v}(t), \quad y(0) = y_{in}, \ t \in (0, T), \quad (15) $$

where $y(t), y_0, \mathbf{v}(t) \in \mathbb{R}^r$, with $r$ being the number of unknowns in space and $\mathcal{L}$ being the discretization of $\Delta$ in space. Since $\mathcal{L} = \mathcal{L}^T$, the optimality system (1) becomes

$$ \dot{y} = \mathcal{L}y - \frac{1}{\alpha}\lambda, \quad \dot{\lambda} = -\mathcal{L}\lambda, \quad (16) $$

with $y(0) = y_{in}$ and $\lambda(T) = y(T) - y_{tg}$.

To obtain an effective preconditioner, we need to study the eigenvalue properties of the matrix $\mathcal{M}$ in (13). To do so, we consider the eigenvalue decomposition $\mathcal{L} = \mathcal{U}\mathcal{D}\mathcal{U}^T$, where $\mathcal{D}$ is a diagonal matrix with the eigenvalues of $\mathcal{L}$ and $\mathcal{U}$ a unitary matrix. The transformation $y(t) \mapsto \mathcal{U}^T \cdot y(t)$ and $\lambda(t) \mapsto \mathcal{U}^T \cdot \lambda(t)$ allows us to diagonalize (16) and obtain scalar equations of the form

$$ \dot{y}(t) = \sigma y(t) - \frac{1}{\alpha}\lambda(t), \quad \dot{\lambda}(t) = -\sigma\lambda(t), \quad t \in (0, T), \quad (17) $$

with $y(t), \lambda(t) \in \mathbb{R}$ and $\sigma \le 0$ an eigenvalue of $\mathcal{L}$ (see [3]). From (17), it is easy to see that $\lambda(t) = e^{\sigma(T-t)}\lambda(T)$ and

$$ y(t) = e^{\sigma t} y_{in} - \frac{1}{\alpha} \int_{0}^{t} e^{\sigma(t-s)} \lambda(s) ds. $$

Substituting $\lambda(t)$ into $y(t)$, we obtain

$$ y(t) = e^{\sigma t} y_{in} - \frac{1}{\alpha} e^{\sigma(T+t)} \left( \int_{0}^{t} e^{-2\sigma s} ds \right) \lambda(T). $$

Evaluating $y$ at $t = T$ and inserting the result into the final condition $\lambda(T) = y(T) - y_{tg}$ of (17) (where $\Lambda_L = \lambda(T)$), we get

$$ f(\sigma)\Lambda_L = e^{\sigma T} y_{in} - y_{tg}, \quad (18) $$

where $f(\sigma) := 1 + (e^{2\sigma T} - 1)/2\alpha\sigma$.

Note that $f(\sigma)$ is an eigenvalue of $\mathcal{M}$ whenever $\sigma$ is an eigenvalue of $\mathcal{L}$. Thus, if $\text{Sp}(\mathcal{L}) = \{\sigma_j, j = 1, \dots, r\}$ is the discrete approximation of the continuous eigenvalues of $\Delta$ given by $-j^2\pi^2$, then $f(\sigma_j) \approx 1 - 1/2\alpha\sigma_j, j = 1, \dots, r$, so that $\mathcal{U}f(\mathcal{D})\mathcal{U}^T \approx (I - \frac{1}{2\alpha}\mathcal{L}^{-1})$. Hence, we propose the following preconditioner

$$ \tilde{\mathcal{M}}^{-1} = \mathcal{L}(\mathcal{L} - \frac{1}{2\alpha}I)^{-1}. \quad (19) $$

Note that there is no longer any mention of the eigenvalue decomposition in the definition of $\tilde{\mathcal{M}}$, meaning that no such computation is needed when applying the preconditioner. Each matrix-vector multiplication of the form $\tilde{\mathcal{M}}^{-1} \cdot \Lambda_L$ only requires a multiplication by $\mathcal{L}$ and the solution of an elliptic problem of the form $(\mathcal{L} - \frac{1}{2\alpha}I)\mathbf{v} = \mathbf{f}$, which can be done cheaply using e.g. algebraic multigrid [17].

The following result shows that for small $\delta t$, the eigenvalues of the preconditioned matrix is clustered around 1. Thus, preconditioned GMRES is expected to converge in a small number of iterations [2].

**Theorem 4.1:** Let a positive integer $N$ be given and $\mathcal{R}_\ell$ be approximated using the Implicit Euler method with $N$ fine sub-intervals over each $(T_{\ell-1}, T_\ell)$. Then any eigenvalue $\mu$ of $\mathcal{M}\tilde{\mathcal{M}}^{-1}$ satisfies

$$ 1 < \mu < 1 + \frac{\delta t}{\alpha}, \quad (20) $$

where $\delta t = T/LN$.

*Proof:* By definition, $\mathcal{R}_\ell(T_\ell)$ for (17) is given by

$$ \mathcal{R}_\ell(T_\ell) = e^{(T+T_\ell)\sigma} \int_{T_\ell-1}^{T_\ell} e^{-2\sigma s} ds. $$

Then, the discrete form of $\mathcal{R}_\ell(T_\ell)$ using Implicit Euler reads $\mathcal{R}_\ell(T_\ell) = \delta t e^{(T+T_\ell-2T_{\ell-1})\sigma} \sum_{n=1}^{N} e^{-2n\sigma\delta t}$, so that,

$$ \mathcal{R}_\ell(T_\ell) = \delta t e^{2(\Delta T - \delta t)\sigma} e^{(L-\ell)\Delta T\sigma} \left( \frac{e^{-2\Delta T\sigma} - 1}{e^{-2\delta t\sigma} - 1} \right). \quad (21) $$

Next, from (11), we have

$$ Y_L = \mathcal{P}_1(T) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_L(T) \Lambda_L - \frac{1}{\alpha} \sum_{j=2}^{L} \mathcal{P}_j(T) \mathcal{R}_{j-1}(T_{j-1}) \Lambda_L. \quad (22) $$

Since $\mathcal{P}_\ell(T) = e^{(T-T_{\ell-1})\sigma} = e^{(L-\ell+1)\Delta T\sigma}$, replacing (21) into (22), we obtain

$$ \begin{aligned} Y_L &= e^{L\Delta T\sigma} \cdot y_{in} \\ &\quad - \frac{\delta t}{\alpha} e^{2(\Delta T - \delta t)\sigma} \left( \frac{e^{-2\Delta T\sigma} - 1}{e^{-2\delta t\sigma} - 1} \right) \left( \frac{e^{2\Delta T\sigma} - 1}{e^{2\Delta T\sigma} - 1} \right) \Lambda_L \\ &= e^{L\Delta T\sigma} \cdot y_{in} - \frac{\delta t}{\alpha} \left( \frac{1 - e^{2\Delta T\sigma}}{1 - e^{2\delta t\sigma}} \right) \Lambda_L. \end{aligned} $$

Therefore, inserting $Y_L$ into $\Lambda_L = Y_L - y_{tg}$, we get $f_{\delta t}(\sigma)\Lambda_L = e^{\sigma T}y_{in} - y_{tg}$, where $f_{\delta t}(\sigma) := 1 + \frac{\delta t}{\alpha} \left( \frac{1 - e^{2\Delta T\sigma}}{1 - e^{2\delta t\sigma}} \right)$. Note that for $\sigma \in \text{Sp}(\mathcal{L})$, $f_{\delta t}(\sigma)$ is eigenvalue of the discrete $\mathcal{M}$.
with $I$ being the identity matrix. Note that $M\Lambda_L$ contains a sum of terms, each one of which can be computed in parallel over the sub-intervals $[T_{\ell-1}, T_{\ell}]$. We remark that the problem of approximating $\mathcal{P}$ and $\mathcal{R}$ is non-trivial: for small and medium sized dense matrices ($r \le 1000$), various methods can be found in [15, 11] and [9]. For large sparse matrices, however, the task is more challenging and frequently requires an understanding of the spectral properties of the underlying operator $\mathcal{L}$. One possibility proposed by [6] is projection-based methods: there, the authors claim that “we clearly find that it is beneficial to use the Arnoldi method rather than a time-stepping method for the propagation of the linear homogeneous problems. Note that the difference between the projection and time-stepping methods gets large for high accuracy,” meaning that the Arnoldi method can produce accurate approximations efficiently.

It remains to solve the linear system (12) using an iterative method, e.g. GMRES. Since the matrix $M$ is generally ill-conditioned, an appropriate preconditioner is needed, which we will construct in the next section for the heat and wave equations.

# 4 Preconditioner design

## 4.1 Heat equation case

In this section, we propose a preconditioner $P^{-1}$ for the heat equation to speed up the convergence of (12) in GMRES. Our approach is inspired by the work in the Master’s thesis [18], where the author derived a preconditioner based the behaviour of $\mathcal{L}$ at low and high frequencies. Let us consider the heat equation given by

$$ \dot{y} = \Delta y + \nu, \text{ on } \Omega \times (0, T), \quad (13) $$

with $y(x, 0) = y_{in}(x)$ on $\Omega$. We consider the 1D case where $\Omega = (0, 1)$ with Dirichlet boundary conditions $y(0, t) = y(1, t) = 0, t \in (0, T)$ with $T \ge 1$. Then, a semi-discretization in space of (13) using finite differences leads the following ODE system

$$ \dot{y}(t) = \mathcal{L}y(t) + \nu(t), \quad y(0) = y_{in}, \quad t \in (0, T), \quad (14) $$

where $y(t), y_0, \nu(t) \in \mathbb{R}^r$ with $r$ being the number of unknowns in space and $\mathcal{L}$ being the discretization of $\Delta$ in space. Since $\mathcal{L}$ is symmetric ($\mathcal{L} = \mathcal{L}^T$), the optimality system (1) becomes

$$ \dot{y} = \mathcal{L}y - \frac{1}{\alpha}\lambda, \quad \dot{\lambda} = -\mathcal{L}\lambda, \quad (15) $$

with $y(0) = y_{in}$ and $\lambda(T) = y(T) - y_{tg}$.

We now consider the eigenvalue decomposition $\mathcal{L} = QDQ^T$, where $D$ is a diagonal matrix with the eigenvalues of $\mathcal{L}$ and $Q$ a unitary matrix. The transformation $y(t) \mapsto Q^T \cdot y(t)$ and $\lambda(t) \mapsto Q^T \cdot \lambda(t)$ allows us to diagonalize (15) to obtain scalar equations of the form

$$ \dot{y}(t) = \sigma y(t) - \frac{1}{\alpha}\lambda(t), \quad \dot{\lambda}(t) = -\sigma\lambda(t), \quad t \in (0, T), \quad (16) $$

with $y(t), \lambda(t) \in \mathbb{R}$ and $\sigma$ an eigenvalue of $\mathcal{L}$, which is real and negative. From (16), it is easy to see that $\lambda(t) = e^{\sigma(T-t)}\lambda(T)$ and

$$ y(t) = e^{\sigma t}y_{in} - \frac{1}{\alpha}\int_{0}^{t}e^{\sigma(t-s)}\lambda(s)ds. $$

Substituting $\lambda(t)$ into $y(t)$ gives

$$ y(t) = e^{\sigma t}y_{in} - \frac{1}{\alpha}e^{\sigma(T+t)}\left(\int_{0}^{t}e^{-2\sigma s}ds\right)\lambda(T). $$

Evaluating $y$ at $t = T$ and inserting the result into the final condition $\lambda(T) = y(T) - y_{tg}$ of (16) (where $\Lambda_L = \lambda(T)$), we get

$$ f(\sigma)\Lambda_L = e^{\sigma T}y_{in} - y_{tg}, \quad (17) $$

where $f(\sigma) := 1 + (e^{2\sigma T} - 1)/2\alpha\sigma$.

Note that $f(\sigma)$ is an eigenvalue of $M$ whenever $\sigma$ is an eigenvalue of $\mathcal{L}$. Thus, if $\text{Sp}(\mathcal{L}) = \{\sigma_j, j = 1, \dots, r\}$ is the discrete approximation of the continuous eigenvalues of $\Delta$ given by $-j^2\pi^2$, then $f(\sigma_j) \approx 1 - 1/2\alpha\sigma_j, j = 1, \dots, r$, so that $Qf(D)Q^T \approx (I - \frac{1}{2\alpha}\mathcal{L}^{-1})$. Hence, we propose the following preconditioner

$$ \widehat{M}^{-1} = \mathcal{L}(\mathcal{L} - \frac{1}{2\alpha}I)^{-1}, \quad (18) $$
The largest eigenvalue of $\mathcal{L}$ as function of $r$ is given by
$\sigma_1(r) = -4(r+1)^2 \sin^2(\frac{\pi}{2(r+1)}), r \ge 1$ (see [3]). Then,

$$
\sigma'_{1}(r) = -8(r+1)\sin\left(\frac{\pi}{2(r+1)}\right) \times \left[\sin\left(\frac{\pi}{2(r+1)}\right) - \frac{\pi}{2(r+1)}\cos\left(\frac{\pi}{2(r+1)}\right)\right].
$$

Since $r \ge 1$, $\frac{\pi}{2(r+1)} \in (0, \frac{\pi}{4}]$. Using the fact that $\tan \theta - \theta > 0$, i.e., $\sin \theta - \theta \cos \theta > 0$ for $0 \le \theta \le \frac{\pi}{4}$, we conclude that $\sigma'_1(r) < 0$. Thus, $\sigma_1(r)$ decreases on $(1, +\infty)$ and its maximum is $\sigma_1(1) = -8$.

Now, let $\psi_0(\sigma) := f_{\delta t}(\sigma)f^{-1}(\sigma)$ for $\sigma \le -8$. Then, $\psi_0(x) - 1 = \varphi(\sigma)/(2\alpha\sigma - 1)(1 - e^{2\delta t\sigma})$, where $\varphi(\sigma) := 2\sigma\delta t(1 - e^{2T\sigma}) + (1 - e^{2\delta t\sigma})$. We have $\varphi'(\sigma) = 2\delta t(1 - e^{2\delta\sigma}) - 2\delta t(2\sigma T + 1)e^{2T\sigma} > 0$, since $(2\sigma T + 1) < 0$. But $\varphi(-8) = -16\delta t(1 - e^{-16T}) + (1 - e^{-16\delta t}) < 0$ for $\delta t > 0$, which means $\varphi(\sigma) < 0$ for $\sigma \le -8$. Thus, $\psi_0(\sigma) - 1 > 0$, i.e., $\psi_0(\sigma) > 1$.

Next, we consider $f_{\delta t}^0(\sigma) := 1 + \frac{\delta t}{\alpha(1-e^{2\delta t}\sigma)} \ge f_{\delta t}(\sigma)$ and $\psi(\sigma) := f_{\delta t}^0(\sigma)f^{-1}(\sigma)$ for $\sigma \in (-\infty, 0)$, whose derivative is $\psi'(\sigma) = \psi_1(\sigma)/\alpha(2\delta t\alpha - 1)^2(1 - e^{2\delta t\sigma})^2$, where $\psi_1(\sigma) := -2\alpha(1 - e^{2\delta t\sigma})^2 - \delta t(1 - e^{2\delta t\sigma}) + 4\delta t^2\sigma(2\alpha\sigma - 1)e^{2\delta t\sigma}$. We have $\psi_1'(\sigma) = 8\delta t e^{2\delta t\sigma} \psi_2(\sigma)$, where $\psi_2(\sigma) := \alpha(1 - e^{2\delta t\sigma}) + \sigma\delta t^2(2\alpha\sigma - 1) + 2\alpha\delta t\sigma$. Since $\psi_2''(\sigma) = 4\alpha\delta t^2(1 - e^{2\delta t\sigma}) > 0$ and $\psi_2'(0) = -\delta t^2 < 0$, we get $\psi_2'(\sigma) < 0$. Also, since $\psi_2(0) = 0$, $\psi_2(\sigma) \ge 0$, which implies $\psi_1'(\sigma) > 0$. Finally, since $\psi_1(0) = 0$, $\psi_1(\sigma) \le 0$ and $\psi'(\sigma) \le 0$. Hence, $\psi$ is decreases on $(-\infty, 0)$ and tends to $1 + \delta t/\alpha$ as $\sigma$ tends to $-\infty$, such that $\psi(\sigma) < 1 + \delta t/\alpha$.

Now, the quantity $\psi_0(\sigma)$ is the eigenvalue of $\mathcal{M}\tilde{\mathcal{M}}^{-1}$ associated with the fixed $\sigma \in \text{Sp}(\mathcal{L})$. Then we have $1 < \psi_0(\sigma) \le \psi(\sigma) \le 1 + \frac{\delta t}{\alpha}$, which concludes the proof. ■

Remark 1: Note that $\tilde{\mathcal{M}}^{-1}$ is derived from the continuous form of $\mathcal{M}$. Therefore, we expect the preconditioner to be more efficient for a high-order approximation of $\mathcal{R}_\ell$ than for a lower-order one.

Analogously, we can prove that (19) yields an effective preconditioner for (13) when the boundary condition in (14) is Dirichlet-Neumann in 1D or Dirichlet in 2D.

In our numerical test, we use the example studied in [4]: we set $T = 1$, $L = 10$, $r = 100$, $y_{in}(x) = \exp(-100(x-1/2)^2)$ and $y_{tg}(x) = \frac{1}{2}\exp(-100(x-1/4)^2) + \frac{1}{2}(-100(x-3/4)^2)$ with Dirichlet boundary conditions. We use Matlab R2021b as our numerical environment. Moreover, to get $\mathcal{R}_\ell$, we use two quadrature formulas, namely Implicit Euler and the Singly Diagonal Implicit Runge-Kutta (SDIRK) method of order 3, whose stages and weights are $c = [1/2 + \sqrt{3}/6, 1/2 - \sqrt{3}/6]$ and $d = [1/2, 1/2]$ respectively.

First, we investigate the numerical behavior of the upper bound in (20) and the maximal eigenvalue $\sigma_{max}$ of $\mathcal{M}\tilde{\mathcal{M}}^{-1}$ for various $N = 100k$, $k = 1, \dots, 10$, i.e., $\delta t = 1/N$ and fixed $\alpha = 10^{-4}$ with only Implicit Euler for getting $\mathcal{R}_\ell$. The result is shown in Fig. 2, where we observe that the upper bound in (20) follows closely the evolution of $\sigma_{max}$ for various $\delta t$. Thus, (20) provides an efficient estimation of the upper bound of the eigenvalues of $\mathcal{M}\tilde{\mathcal{M}}^{-1}$.

Next, we investigate the convergence speed of the linear system (13) in GMRES for $N = 1000$ and $\alpha = \{10^{-4}, 10^{-6}\}$. We use for this instance the function `gmres` in Matlab with `restart=1`, `tol=1e-8` and `maxit=500`.

For $\alpha = 10^{-4}$, the result is shown in Table I. We observe that when solving (13), unpreconditioned GMRES reaches 500 iterations without satisfying the tolerance for either quadrature formula. Moreover, the eigenvalues of $\mathcal{M}$ for both quadrature formulas lie in (1,500). In contrast, for the preconditioned case, when SDIRK is used, GMRES reaches the tolerance with only 2 iterations, and 9 iterations for the Implicit Euler case. We observe a similar behavior for $\alpha = 10^{-6}$, with poor convergence without preconditioning, fairly good convergence for Implicit Euler (44 iterations) and excellent convergence for SDIRK (4 iterations). The results are shown in Table II.

<figure><img src="image_5.png" /><figcaption>Fig. 2. The maximal eigenvalue of $\mathcal{M}\tilde{\mathcal{M}}^{-1}$ and the upper bound in (20) for various $\delta t$ with $\alpha = 10^{-4}$.</figcaption></figure>

<table><caption>TABLE I
GMRES FOR SOLVING (13) FOR THE HEAT EQUATION, WITH $\alpha = 10^{-4}$.</caption>
  <caption>Matrix M</caption>
  <thead>
    <tr>
      <th></th>
      <th>σ<sub>min</sub></th>
      <th>σ<sub>max</sub></th>
      <th># Iters</th>
      <th>Res</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Euler</th>
      <td>2.0</td>
      <td>4.9e2</td>
      <td>500</td>
      <td>1.3e-8</td>
    </tr>
    <tr>
      <th>SDIRK</th>
      <td>1.08</td>
      <td>4.9e2</td>
      <td>500</td>
      <td>4.08e-7</td>
    </tr>
  </tbody>
</table>

<table><caption>TABLE II
GMRES FOR SOLVING (13) FOR THE HEAT EQUATION, WITH $\alpha = 10^{-6}$.</caption>
  <caption>Preconditioned matrix M̃M<sup>-1</sup></caption>
  <thead>
    <tr>
      <th></th>
      <th>σ<sub>min</sub></th>
      <th>σ<sub>max</sub></th>
      <th># Iters</th>
      <th>Res</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Euler</th>
      <td>1.0</td>
      <td>1.78</td>
      <td>9</td>
      <td>7.7e-9</td>
    </tr>
    <tr>
      <th>SDIRK</th>
      <td>0.97</td>
      <td>1.0</td>
      <td>2</td>
      <td>9.64e-9</td>
    </tr>
  </tbody>
</table>

<table>
  <caption>Matrix M</caption>
  <thead>
    <tr>
      <th></th>
      <th>σ<sub>min</sub></th>
      <th>σ<sub>max</sub></th>
      <th># Iters</th>
      <th>Res</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Euler</th>
      <td>1e2</td>
      <td>4.97e4</td>
      <td>500</td>
      <td>1.08e-6</td>
    </tr>
    <tr>
      <th>SDIRK</th>
      <td>9.68</td>
      <td>4.96e4</td>
      <td>500</td>
      <td>4.2e-5</td>
    </tr>
  </tbody>
</table>

<table>
  <caption>Preconditioned matrix M̃M<sup>-1</sup></caption>
  <thead>
    <tr>
      <th></th>
      <th>σ<sub>min</sub></th>
      <th>σ<sub>max</sub></th>
      <th># Iters</th>
      <th>Res</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Euler</th>
      <td>1.0</td>
      <td>7.76</td>
      <td>44</td>
      <td>8.37e-9</td>
    </tr>
    <tr>
      <th>SDIRK</th>
      <td>0.74</td>
      <td>1.0</td>
      <td>4</td>
      <td>7.31e-9</td>
    </tr>
  </tbody>
</table>

In both Tables I and II, we observe that preconditioned GMRES converges more quickly for SDIRK than for Implicit
where I is the identity matrix. Each application of $\tilde{\mathcal{M}}^{-1}$ thus requires the solution of the elliptic problem $(\mathcal{L} - \frac{1}{\alpha}I)v = f$, which is cheaper in terms of numerical computing cost than the inverse of the full matrix $M$.

To understand the behaviour of GMRES applied to the preconditioned system, we prove the following result.

**Theorem 4.1.** Let $N$ be given and $\mathcal{R}_\ell$ be approximated using implicit Euler with $N$ fine sub-intervals over each $[T_{\ell-1}, T_\ell]$. Then any eigenvalue $\mu$ of $\tilde{\mathcal{M}}^{-1}$ satisfies

$$1 < \mu < 1 + \frac{\delta t}{\alpha}, \qquad (19)$$

where $\delta t = T/LN$.

*Proof.* By definition, $\mathcal{R}_\ell(T_\ell)$ for (16) is given by

$$\mathcal{R}_\ell(T_\ell) = e^{(T+T_\ell)} \int_{T_{\ell-1}}^{T_\ell} e^{-2\sigma s} ds.$$

Then, the discrete form of $\mathcal{R}_\ell(T_\ell)$ using implicit Euler reads

$$\mathcal{R}_\ell(T_\ell) = \delta t e^{(T+T_\ell-2T_{\ell-1})\sigma} \sum_{n=1}^{N} e^{-2n\sigma\delta t},$$

so that,

$$\mathcal{R}_\ell(T_\ell) = \delta t e^{2(\Delta T - \delta t)\sigma} \cdot e^{(L-\ell)\Delta T\sigma} \left( \frac{e^{-2\Delta T\sigma} - 1}{e^{-2\delta t\sigma} - 1} \right). \qquad (20)$$

Next, from (5), we have

$$Y_L = \mathcal{P}_1(T) \cdot y_{in} - \frac{1}{\alpha} \mathcal{R}_L(T) \Lambda_L - \frac{1}{\alpha} \sum_{j=2}^{L} \mathcal{P}_j(T) \mathcal{R}_{j-1}(T_{j-1}) \Lambda_L. \qquad (21)$$

Since $\mathcal{P}_\ell(T) = e^{(T-T_{\ell-1})\sigma} = e^{(L-\ell+1)\Delta T\sigma}$, we obtain using (20) that

$$\begin{aligned} Y_L &= e^{L\Delta T\sigma} \cdot y_{in} - \frac{\delta t}{\alpha} e^{2(\Delta T - \delta t)\sigma} \left( \frac{e^{-2\Delta T\sigma} - 1}{e^{-2\delta t\sigma} - 1} \right) \left( \frac{e^{2T\sigma} - 1}{e^{2\Delta T\sigma} - 1} \right) \Lambda_L \\ &= e^{L\Delta T\sigma} \cdot y_{in} - \frac{\delta t}{\alpha} \left( \frac{1 - e^{2T\sigma}}{1 - e^{2\delta t\sigma}} \right) \Lambda_L. \end{aligned}$$

Inserting $Y_L$ into $\Lambda_L = Y_L - y_{tg}$ yields $f_{\delta t}(\sigma)\Lambda_L = e^{T\sigma}y_{in} - y_{tg}$, where $f_{\delta t}(\sigma) := 1 + \frac{\delta t}{\alpha} \left( \frac{1-e^{2T\sigma}}{1-e^{2\delta t\sigma}} \right)$. Note that for $\sigma \in \text{Sp}(\mathcal{L})$, $f_{\delta t}(\sigma)$ is eigenvalue of $M$.

The largest eigenvalue of $\mathcal{L}$ as function of $r$ is given by $\sigma_1(r) = -4(r+1)^2 \sin^2(\frac{\pi}{2(r+1)})$, $r \ge 1$ (see [2]). Then, we have

$$\sigma'_1(r) = -8(r+1) \sin \left( \frac{\pi}{2(r+1)} \right) \left[ \sin \left( \frac{\pi}{2(r+1)} \right) - \frac{\pi}{2(r+1)} \cos \left( \frac{\pi}{2(r+1)} \right) \right].$$

Note that $\frac{\pi}{2(r+1)} \in (0, \frac{\pi}{4}]$ since $r \ge 1$. Then, using the fact that $\tan\theta - \theta > 0$, i.e., $\sin\theta - \theta \cos\theta > 0$ for $0 \le \theta \le \frac{\pi}{4}$, we conclude that $\sigma'_1(r) < 0$. Hence, $\sigma_1(r)$ is decreasing as $r$ increases and the maximum of $\sigma_1(r)$ for $r \ge 1$ is $-8$.

Now, let consider $\psi_0(\sigma) = f_{\delta t}(\sigma)f^{-1}(\sigma)$ for $\sigma \le -8$. Then, we have $\psi_0(x) - 1 = \varphi(\sigma)/(2\alpha\sigma - 1)(1 - e^{2\delta t\sigma})$, where $\varphi(\sigma) := 2\sigma\delta t(1 - e^{2T\sigma}) + (1 - e^{2\delta t\sigma})$. Next, $\varphi'(\sigma) = 2\delta t(1 - e^{2\delta\sigma}) - 2\delta t(2\sigma T + 1)e^{2T\sigma} > 0$, since $(2\sigma T + 1) < 0$. Moreover, $\varphi(-8) = -16\delta t(1 - e^{-16T}) + (1 - e^{-16\delta t}) < 0$ for $\delta t > 0$ then, $\varphi(\sigma) < 0$. It follows that, $\psi_0(\sigma) - 1 > 0$, i.e., $\psi_0(\sigma) > 1$.

Next, since $f_{\delta t}(\sigma) > 0$, let $f_{\delta t}(\sigma) \le f_{\delta t}^0(\sigma) := 1 + \frac{\delta t}{\alpha(1-e^{2\delta t\sigma})}$ and $\psi(\sigma) := f_{\delta t}^0(\sigma)f^{-1}(\sigma)$ for $\sigma \in (-\infty, 0)$. The derivative of $\psi$ is given by $\psi'(\sigma) = \psi_1(\sigma)/(2\delta t\alpha - 1)^2(1 - e^{2\delta t\sigma})^2$, where

$$\psi_1(\sigma) := -2\alpha(1 - e^{2\delta t\sigma})^2 - \delta t(1 - e^{2\delta t\sigma}) + 4\delta t^2\sigma(2\alpha\sigma - 1)e^{2\delta t\sigma}.$$

Moreover, $\psi'_1(\sigma) = 8\delta t e^{2\delta t\sigma}\psi_2(\sigma)$, where $\psi_2(\sigma) := \alpha(1 - e^{2\delta t\sigma}) + \sigma\delta t^2(2\alpha\sigma - 1) + 2\alpha\delta t\sigma$. Also, $\psi''_2(\sigma) = 4\alpha\delta t^2(1 - e^{2\delta t\sigma}) > 0$ and $\psi'_2(0) = -\delta t^2 < 0$ then $\psi'_2(\sigma) < 0$. Since $\psi_2(0) = 0$, $\psi_2(\sigma) \ge 0$ which implies that
Euler. This is consistent with Remark 1, since SDIRK is of order 3, which is higher than the order 1 of Implicit Euler.

Theorem 4.1 asserts that the condition number $cond(M\tilde{M}^{-1})$ is bounded by $1 + \frac{\delta_l}{\alpha}$ when Implicit Euler is used to compute $\mathcal{R}_l$. Since $M\tilde{M}^{-1}$ is symmetric (as $M$ and $\tilde{M}$ commute and are both symmetric), we deduce that the number of preconditioned GMRES iterations remains bounded as $r \to \infty$, which is confirmed by the results in Table III. We also observe that GMRES converges faster for a higher order method like SDIRK than for Implicit Euler.

<table><caption>TABLE III
GMRES FOR SOLVING (13) WITH PRECONDITIONING FOR THE HEAT EQUATION, WITH $\alpha = 10^{-4}$ AND FOR VARIOUS $r$.</caption><thead><tr><th rowspan="2">r</th><th colspan="2"># Iters(Euler)</th><th colspan="2"># Iters(SDIRK)</th></tr><tr><th>L = 10<sup>3</sup></th><th>L = 3 · 10<sup>3</sup></th><th>L = 10<sup>3</sup></th><th>L = 3 · 10<sup>3</sup></th></tr></thead><tbody><tr><td>100</td><td>9</td><td>6</td><td>3</td><td>2</td></tr><tr><td>200</td><td>10</td><td>7</td><td>3</td><td>3</td></tr><tr><td>250</td><td>11</td><td>7</td><td>3</td><td>3</td></tr><tr><td>600</td><td>11</td><td>7</td><td>3</td><td>3</td></tr></tbody></table>

Finally, one can use an explicit method as a numerical integrator to approximate $\mathcal{R}_l$, as long as a CFL condition is satisfied.

## B. Wave equation

We now consider the optimal control problem constrained by the 1D wave equation given by

$$ \partial_{tt}u = \Delta u + v, \quad \text{on } (0,1) \times (0,T), \quad (23) $$

with $u(x, 0) = u_0(x)$, $\partial_t u(x, 0) = 0$, $x \in (0, 1)$. For simplicity, we will consider homogeneous Dirichlet boundary conditions in space. Discretizing in space using second-order centered finite differences leads to the ODE system

$$ \dot{y} = \mathcal{L}y + \mathcal{B}v, \quad \text{with } \mathcal{L} = \begin{bmatrix} 0 & I \\ \Delta_h & 0 \end{bmatrix}, \quad \mathcal{B} = \begin{bmatrix} 0 \\ I \end{bmatrix}, \quad (24) $$

where $y = \begin{bmatrix} u_h \\ \partial_t u_h \end{bmatrix}$ and $\Delta_h$ is the discrete Laplacian.

To derive a preconditioner for the system (13), we again look to the continuous problem. Since $\mathcal{L}$ is no longer symmetric, the continuous analogue of $M$ takes the form

$$ \mathcal{M}_c := I + \frac{1}{\alpha} \int_0^T e^{s\mathcal{L}} \mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} ds. $$

Let us evaluate this integral. First note that

$$ e^{s\mathcal{L}} = I + s\mathcal{L} + \frac{s^2}{2!}\mathcal{L}^2 + \frac{s^3}{3!}\mathcal{L}^3 + \dots $$

Introducing the notation $A := -\Delta_h$, we see that

$$ \mathcal{L}^2 = \begin{bmatrix} \Delta_h & 0 \\ 0 & \Delta_h \end{bmatrix} = - \begin{bmatrix} A & 0 \\ 0 & A \end{bmatrix} := -\mathcal{A} $$

is a symmetric negative definite matrix, so we can write

$$ \begin{aligned} e^{s\mathcal{L}} &= \left( I - \frac{s^2}{2!}\mathcal{A} + \frac{s^4}{4!}\mathcal{A}^2 + \dots \right) \\ &\quad + s\mathcal{L} \left( I - \frac{s^2}{3!}\mathcal{A} + \frac{s^4}{5!}\mathcal{A}^2 + \dots \right) \\ &= \cos(s\mathcal{A}^{1/2}) + s\mathcal{L}\operatorname{sinc}(s\mathcal{A}^{1/2}), \end{aligned} $$

where $\operatorname{sinc}(x) = \sin(x)/x$ is the sinc function and $\mathcal{A}^{1/2}$ is the symmetric square root of $\mathcal{A}$ (i.e., a symmetric matrix with positive eigenvalues whose square is $\mathcal{A}$). This leads to

$$ \begin{aligned} e^{s\mathcal{L}} \mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} = & [\cos(s\mathcal{A}^{1/2}) + s\mathcal{L}\operatorname{sinc}(s\mathcal{A}^{1/2})] \cdot \mathcal{B} \\ & \cdot \mathcal{B}^T \cdot [\cos(s\mathcal{A}^{1/2}) + s\operatorname{sinc}(s\mathcal{A}^{1/2})\mathcal{L}^T]. \end{aligned} $$

By direct calculation and using $x \cdot \operatorname{sinc}(xA^{1/2}) = A^{-1/2}\sin(xA^{1/2})$, we obtain

$$ e^{s\mathcal{L}} \mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} = \begin{bmatrix} \mathcal{C}_{11}(s) & \mathcal{C}_{12}(s) \\ \mathcal{C}_{21}(s) & \mathcal{C}_{22}(s) \end{bmatrix}, $$

where $\mathcal{C}_{11}(s) := A^{-1}\sin^2(sA^{1/2})$, $\mathcal{C}_{22}(s) := \cos^2(sA^{1/2})$, $\mathcal{C}_{12}(s) := A^{-1/2}\cos(sA^{1/2})\sin(sA^{1/2})$ and $\mathcal{C}_{21}(s) := \mathcal{C}_{12}(s)$. We can now integrate these coefficients to obtain

$$ \mathcal{M}_c = \begin{bmatrix} \mathcal{M}_{11} & \mathcal{M}_{12} \\ \mathcal{M}_{21} & \mathcal{M}_{22} \end{bmatrix}, $$

where $\mathcal{M}_{12} := \frac{1}{2\alpha}A^{-1}(I - \cos(2TA^{1/2}))$, $\mathcal{M}_{21} := \mathcal{M}_{12}$,

$$ \mathcal{M}_{11} := I + \frac{T}{2\alpha}A^{-1} - \frac{1}{4\alpha}A^{-3/2}\sin(2TA^{1/2}), $$

$$ \mathcal{M}_{22} := \left(1 + \frac{T}{2\alpha}\right)I + \frac{1}{4\alpha}A^{-1/2}\sin(2TA^{1/2}). $$

To find a good preconditioner for $\mathcal{M}_c$, let us analyze $\mathcal{M}_c^{-1}$. Since all $\mathcal{M}_{ij}$ are square and they all commute, we have

$$ \mathcal{M}_c^{-1} = \begin{bmatrix} \mathcal{H}^{-1}\mathcal{M}_{22} & -\mathcal{H}^{-1}\mathcal{M}_{12} \\ -\mathcal{H}^{-1}\mathcal{M}_{21} & \mathcal{H}^{-1}\mathcal{M}_{11} \end{bmatrix}, $$

with $\mathcal{H} := \mathcal{M}_{11}\mathcal{M}_{22} - \mathcal{M}_{12}\mathcal{M}_{21}$. Just like for the heat equation, we replace $A$ by a scalar $\sigma > 0$ and study the eigenvalues of $\mathcal{M}_c^{-1}$ when $\sigma \approx k^2\pi^2$ approximates the eigenvalues of the negative Laplacian. We observe for large $k$ that $\mathcal{H} \approx \frac{1}{4\alpha^2}(2\alpha + T)A^{-1}(TI + 2\alpha A)$, $\mathcal{H}^{-1}\mathcal{M}_{12} = \mathcal{H}^{-1}\mathcal{M}_{21} \approx 0$ and

$$ \bullet \mathcal{H}^{-1}\mathcal{M}_{11} \approx I - (TI + 2\alpha A)^{-1}, $$

$$ \bullet \mathcal{H}^{-1}\mathcal{M}_{22} \approx \left(1 - \frac{1}{(2\alpha + T)}\right)I. $$

These observations lead us to approximate $\mathcal{M}_c^{-1}$ by a block diagonal matrix $\tilde{\mathcal{M}}^{-1}$ of the form

$$ \tilde{\mathcal{M}}^{-1} := I - \begin{bmatrix} (aI + bA)^{-1} & 0 \\ 0 & cI \end{bmatrix}, \quad (25) $$

and use it as a preconditioner for (13), where $a = 1, b = 2\alpha/T$ and finally $c = T/(2\alpha + T)^1$. As before, one need not compute the eigenvalue decomposition in order to apply the preconditioner $\tilde{\mathcal{M}}$: each application only requires solving an elliptic problem of the form $(aI + bA)v = f$, which again can be done cheaply using e.g. multigrid.

In our numerical experiment, we take $T = 2.3$, $u_0(x) = \exp(-100(x - 1/2)^2)$, so that $y_{in} = (u_0(x), 0)$ and $y_{tg}(x) = [(\frac{1}{2}\exp(-100(x - 1/4)^2) + \frac{1}{2}(-100(x - 3/4)^2), 0]$. We set $L = 10, r = 100$, and we use the same Matlab functions for computing $\mathcal{R}_l$ and $\mathcal{P}_l$ as in Section IV-A. We first investigate the convergence speed of the linear system (13)

<footnote n="1">¹The parameters $a, b$, and $c$ must be as specified in this paper for the general case, rather than in the IEEE conference paper.</footnote>
$\psi_1'(\sigma) > 0$. And, since $\psi_1(0) = 0$, $\psi_1(\sigma) \le 0$ and then $\psi'(\sigma) \le 0$. Hence, $\psi$ is decreasing in $(-\infty, 0)$ and tends to $1 + \delta t/\alpha$ as $\sigma$ tends to $-\infty$, such that $\psi(\sigma) < 1 + \delta t/\alpha$.

Now, the quantity $\psi_0(\sigma)$ is the eigenvalue of $M\tilde{M}^{-1}$ associated to the fixed $\sigma \in \text{Sp}(\mathcal{L})$. Then we have $1 < \psi_0(\sigma) \le \psi(\sigma) \le 1 + \frac{\delta t}{\alpha}$, which concludes the proof. $\square$

**Remark 1.** Note that $\tilde{M}^{-1}$ is derived from the continuous form of $M$. Therefore, we expect the preconditioner to be more efficient for a high-order approximation of $\mathcal{R}_\ell$ than for a lower order one.

Analogously, we can also easily prove that (18) is efficient for (12) when the boundary condition of the Laplacian in (13) is Dirichlet-Neumann in 1D or Dirichlet-Dirichlet boundary condition in 2D.

In our numerical test, we use the example studied in [4]: we set $T = 1$, $L = 10$, $r = 100$, $y_{in}(x) = \exp(-100(x - 1/2)^2)$ and $y_{tg}(x) = \frac{1}{2}\exp(-100(x-1/4)^2) + \frac{1}{2}(-100(x-3/4)^2)$ with Dirichlet boundary conditions. We use Matlab R2021b as our numerical environment. The operator $\mathcal{P}_\ell$ is computed by the Matlab function `expm`. As $\mathcal{L}$ is full rank, `expm(\mathcal{L})$ computes the approximations $\tilde{Q}$ of $Q$ and $\tilde{D}$ of $D$ then approximates $\exp(\mathcal{L})$ by $\tilde{Q}\exp(\tilde{D})\tilde{Q}^T$. Moreover, to get $\mathcal{R}_\ell$, we use two quadrature formulas. The first one is Implicit Euler and the second is the Singly Diagonal Implicit Runge-Kutta (SDIRK) method of order 3, whose stages and weights are $c = [1/2 + \sqrt{3}/6, 1/2 - \sqrt{3}/6]$ and $d = [1/2, 1/2]$ respectively.

First, we investigate the numerical behavior of the upper bound in (19) and the maximal eigenvalue $\sigma_{max}$ of $MP^{-1}$ for various $N = 1000$, $k = 1, \dots, 10$, i.e., $\delta t = 1/N$ and fixed $\alpha = 10^{-4}$ with only Implicit Euler for getting $\mathcal{R}_\ell$. The result is shown in Fig. 2, wherein the upper bound in (19) follows efficiently the evolution of $\sigma_{max}$ for various $\delta t$. So, (19) provides an efficient estimation of the upper bound of the eigenvalues of $MP^{-1}$.

Next, we investigate the convergence speed of the linear system (12) in GMRES for $N = 1000$ and $\alpha = \{10^{-4}, 10^{-6}\}$. We use for this instance the function `gmres` in Matlab with `restart=1`, `tol=1e-8` and `maxit=500`.

For $\alpha = 10^{-4}$, the result is shown in Table 1. We observe that for solving (12), GMRES without preconditioner reaches 500 iterations without satisfying the tolerance for both quadrature formulas. Moreover, the eigenvalues of $M$ for both quadrature formulas lie in (1, 500). But, for the preconditioned case, when SDIRK is used, GMRES reaches the tolerance with only 2 iterations, and 9 iterations for Implicit Euler case.

<figure><img src="image_4.png" /><figcaption>Figure 2: The maximal eigenvalue of $M\tilde{M}^{-1}$ and the upper bound in (19) for various $\delta t$ when $\alpha = 10^{-4}$.</figcaption></figure>

For $\alpha = 10^{-6}$, the result is shown in Table 2. We observe that GMRES without preconditioner reaches 500 iterations without satisfying the tolerance for both quadrature formulas. Moreover, the eigenvalues of $M$ increase compare to the ones in $\alpha = 10^{-4}$ case. Besides, for the preconditioned system when SDIRK is used GMRES convergence in 4 iterations whereas for Implicit Euler, 44 iterations.

Furthermore, for both Tables 1 and 2, we observe that GMRES for the preconditioned system for SDIRK case converges faster than the Implicit Euler case. This follows from Remark 1 because the order of SDIRK is 3 higher than the order 1 of Implicit Euler. We also observe that the maximal eigenvalue of $M\tilde{M}^{-1}$ for Implicit Euler method is less than for SDIRK, meaning that the upper bound in (19) can be lower for high order methods than for Implicit Euler.
in GMRES for $N = 1000$ and $\alpha = 10^{-6}$. We use for this instance the function gmres in Matlab with restart=[, tol=1e-8 and maxit=size($M$, 1). The results are shown in Table IV, where we observe that unpreconditioned GMRES applied to (13) converges in 84 iterations for both quadrature formulas, whereas the preconditioned iteration converges in only 4 iterations.

<table><caption>TABLE IV
GMRES FOR SOLVING (13) FOR THE WAVE EQUATION, WITH $\alpha = 10^{-6}$.</caption><thead><tr><th colspan="3">Matrix M</th></tr><tr><th></th><th>cond(M)</th><th># Iters</th><th>Res</th></tr></thead><tbody><tr><td>Euler</td><td>3.8e4</td><td>84</td><td>8.74e-9</td></tr><tr><td>SDIRK</td><td>3.8e4</td><td>84</td><td>8.74e-9</td></tr></tbody></table><table><thead><tr><th colspan="3">Preconditioned matrix M(M<sup>-1</sup>)</th></tr><tr><th></th><th>cond(M(M<sup>-1</sup>)</th><th># Iters</th><th>Res</th></tr></thead><tbody><tr><td>Euler</td><td>2.59</td><td>4</td><td>1.58e-9</td></tr><tr><td>SDIRK</td><td>2.59</td><td>4</td><td>1.58e-9</td></tr></tbody></table>

Unlike for the heat equation, the higher order of SDIRK does not lead to better convergence compared to Implicit Euler. Nonetheless, for both methods, we have observed that the number of GMRES iterations does not change as we increase the number of time steps $N$ from 100 to 1000, as long as preconditioning is used.

These results show that for fixed $r$, the discretization parameters in time do not change the behavior of the convergence of (13) in GMRES. We also observe that when $r$ increases, the number of iterations of the unpreconditioned system increases, whereas the number of iterations of the preconditioned system decreases. The results are shown in Tables V and VI.

<table><caption>TABLE V
GMRES FOR SOLVING (13) FOR THE WAVE EQUATION WITH IMPLICIT EULER WITH $\alpha = 10^{-6}$ FOR VARIOUS $r$.</caption><thead><tr><th>r</th><th>cond(M)</th><th># Iters</th><th>cond(M(M<sup>-1</sup>)</th><th># Iters</th></tr></thead><tbody><tr><td>10</td><td>4.83e2</td><td>10</td><td>2.47</td><td>5</td></tr><tr><td>150</td><td>7.75e4</td><td>76</td><td>2.81</td><td>3</td></tr><tr><td>350</td><td>2.48e5</td><td>104</td><td>2.49</td><td>3</td></tr></tbody></table>

<table><caption>TABLE VI
GMRES FOR SOLVING (13) FOR THE WAVE EQUATION WITH IMPLICIT EULER FOR VARIOUS $\alpha$.</caption><thead><tr><th>α</th><th>cond(M)</th><th># Iters</th><th>cond(M(M<sup>-1</sup>)</th><th># Iters</th></tr></thead><tbody><tr><td>1e-5</td><td>2.26e4</td><td>52</td><td>2.41</td><td>3</td></tr><tr><td>1e-3</td><td>4.98e2</td><td>20</td><td>1.09</td><td>3</td></tr><tr><td>1e-1</td><td>6.03</td><td>9</td><td>1.01</td><td>3</td></tr><tr><td>1e1</td><td>1.05</td><td>4</td><td>1.0</td><td>2</td></tr></tbody></table>

V. ONGOING WORK

We are currently studying the behavior of our preconditioners when $M$ is obtained from an explicit method while respecting the CFL condition, in order to better understand the impact of using explicit versus implicit methods. We are also working on more general convergence properties of the

algorithm, its error analysis, and on understanding its performance compared to existing parallel-in-time algorithms.

VI. ACKNOWLEDGMENTS

The work described in this paper is partially supported by a grant from the ANR/RGC joint research scheme sponsored by the Research Grants Council of the Hong Kong Special Administrative Region, China (Project no. A-CityU203/19), by the French National Research Agency (Project ALLOWAPP, grant ANR-19-CE46-0013-01 and Projet DEEPNUM, grant ANR-21-CE23-0017) and by Mitacs-Globalink Canada and INRIA France (grant IT34057, 2023).

REFERENCES

[1] Andrew T Barker and Martin Stoll. Domain decomposition in time for PDE-constrained optimization. *Computer Physics Communications*, 197:136-143, 2015.

[2] Michele Benzi. Preconditioning techniques for large linear systems: a survey. *Journal of Computational Physics*, 182(2):418-477, 2002.

[3] Albrecht Böttcher and Sergei M. Grudsky. *Spectral Properties of Banded Toeplitz Matrices*. Society for Industrial and Applied Mathematics, 2005.

[4] M. J. Gander, F. Kwok, and J. Salomon. ParaOpt: A parareal algorithm for optimality systems. *SIAM Journal on Scientific Computing*, 42(5):A2773-A2802, 2020.

[5] Martin J Gander. 50 years of time parallel time integration. In *Multiple Shooting and Time Domain Decomposition Methods: MuSTDD*, Heidelberg, May 6-8, 2013, pages 69-113. Springer, 2015.

[6] Martin J. Gander and Stefan Güttel. Paraexp: A parallel integrator for linear initial-value problems. *SIAM Journal on Scientific Computing*, 35(2):C123-C142, 2013.

[7] Martin J Gander and Laurence Halpern. Optimized schwarz waveform relaxation methods for advection reaction diffusion problems. *SIAM Journal on Numerical Analysis*, 45(2):666-697, 2007.

[8] Martin J Gander and Felix Kwok. Schwarz methods for the time-parallel solution of parabolic control problems. In *Domain decomposition methods in science and engineering XXII*, pages 207-216. Springer, 2016.

[9] G.H. Golub and C.F. Van Loan. *Matrix computations*. The John Hopkins University Press, Baltimore, 2013.

[10] Matthias Heinkenschloss. A time-domain decomposition iterative method for the solution of distributed linear quadratic optimal control problems. *Journal of Computational and Applied Mathematics*, 173(1):169-198, 2005.

[11] Nicholas J. Higham. *Functions of Matrices*. Society for Industrial and Applied Mathematics, 2008.

[12] Thi-Thao-Phuong Hoang, Jérôme Jaffré, Caroline Japhet, Michel Kern, and Jean E Roberts. Space-time domain decomposition methods for diffusion problems in mixed formulations. *SIAM Journal on Numerical Analysis*, 51(6):3532-3559, 2013.

[13] Shishun Li and Xiao-Chuan Cai. Convergence analysis of two-level space-time additive Schwarz method for parabolic equations. *SIAM Journal on Numerical Analysis*, 53(6):2727-2751, 2015.

[14] Yvon Maday, Julien Salomon, and Gabriel Turinici. Monotonic parareal control for quantum systems. *SIAM Journal on Numerical Analysis*, 45(6):2468-2482, 2007.

[15] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. *SIAM Review*, 45(1):3-49, 2003.

[16] Benjamin W Ong and Jacob B Schroder. Applications of time parallelization. *Computing and Visualization in Science*, 23:1-15, 2020.

[17] John W Ruge and Klaus Stüben. Algebraic multigrid. In *Multigrid methods*, pages 73-130. SIAM, 1987.

[18] Andrea Toselli and Olof Widlund. *Domain Decomposition Methods – Algorithms and Theory*. Springer-Verlag, Berlin, 2005.

[19] Colin S.C. Tsang. Preconditioners for linear parabolic optimal control problems. Master's thesis, Hong Kong Baptist University, Hong Kong, China, August 2017.
<table><caption>Table 1: GMRES for solving (12) for $\alpha = 10^{-4}$.</caption><thead><tr><th colspan="5">Matrix M</th></tr><tr><th>Schemes</th><th>σ<sub>min</sub></th><th>σ<sub>max</sub></th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>2.0</td><td>4.9e2</td><td>500</td><td>1.3e-8</td></tr><tr><td>SDIRK</td><td>1.08</td><td>4.9e2</td><td>500</td><td>4.08e-7</td></tr></tbody></table>

<table><caption>Table 2: GMRES for solving (12) for $\alpha = 10^{-6}$.</caption><thead><tr><th colspan="5">Preconditioned matrix MP<sup>-1</sup></th></tr><tr><th>Schemes</th><th>σ<sub>min</sub></th><th>σ<sub>max</sub></th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>1.0</td><td>1.78</td><td>9</td><td>7.7e-9</td></tr><tr><td>SDIRK</td><td>0.97</td><td>1.0</td><td>2</td><td>9.64e-9</td></tr></tbody></table>

<table><thead><tr><th colspan="5">Matrix M</th></tr><tr><th>Schemes</th><th>σ<sub>min</sub></th><th>σ<sub>max</sub></th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>1e2</td><td>4.97e4</td><td>500</td><td>1.08e-6</td></tr><tr><td>SDIRK</td><td>9.68</td><td>4.96e4</td><td>500</td><td>4.2e-5</td></tr></tbody></table>

<table><thead><tr><th colspan="5">Preconditioned matrix MP<sup>-1</sup></th></tr><tr><th>Schemes</th><th>σ<sub>min</sub></th><th>σ<sub>max</sub></th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>1.0</td><td>7.76</td><td>44</td><td>8.37e-9</td></tr><tr><td>SDIRK</td><td>0.74</td><td>1.0</td><td>4</td><td>7.31e-9</td></tr></tbody></table>

Finally, one can use an explicit method as quadrature formula to approximate $R_\ell$ as long as a CFL condition is satisfied, just like in the implicit method case.

## 4.2 Wave equation case

We now consider the optimal control problem constrained by the 1D wave equation given by

$$u_{tt} = \Delta u + \nu, \text{ on } (0, 1) \times (0, T), \qquad (22)$$

with $u(x, 0) = u_0(x)$, $u_t(x, 0) = 0$, $x \in (0, 1)$. We can rewrite (22) as the ODE system

$$\dot{y} = \mathcal{L}y + \mathcal{B}\nu, \qquad (23)$$

with

$$y = \begin{bmatrix} u \\ u_t \end{bmatrix}, \quad \mathcal{L} = \begin{bmatrix} 0 & I \\ \Delta & 0 \end{bmatrix} \quad \text{and} \quad \mathcal{B} = \begin{bmatrix} 0 \\ I \end{bmatrix}.$$

For simplicity, we will consider homogeneous Dirichlet boundary conditions in space. In order to derive a preconditioner for the system (12), we will again look to the continuous problem. Since $\mathcal{L}$ is no longer symmetric, the continuous analogue of $M$ takes the form

$$\mathcal{M} := I + \frac{1}{\alpha} \int_{0}^{T} \exp(s\mathcal{L})\mathcal{B}\mathcal{B}^{T} \exp(s\mathcal{L}^{T}) ds.$$

Let us evaluate this integral. First note that

$$\exp(s\mathcal{L}) = I + s\mathcal{L} + \frac{s^2}{2!}\mathcal{L}^2 + \frac{s^3}{3!}\mathcal{L}^3 + \dots$$

Introducing the notation $A = -\Delta$, we see that

$$\mathcal{L}^2 = \begin{bmatrix} \Delta & 0 \\ 0 & \Delta \end{bmatrix} = - \begin{bmatrix} A & 0 \\ 0 & A \end{bmatrix} =: -\mathcal{A}$$

is a symmetric negative definite matrix, so we can write

$$\begin{aligned} e^{s\mathcal{L}} &= \left( I - \frac{s^2}{2!}\mathcal{A} + \frac{s^4}{4!}\mathcal{A}^2 + \dots \right) + s\mathcal{L} \left( I - \frac{s^2}{3!}\mathcal{A} + \frac{s^4}{5!}\mathcal{A}^2 + \dots \right) \\ &= \cos(s\mathcal{A}^{1/2}) + s\mathcal{L}\sin(s\mathcal{A}^{1/2}), \end{aligned}$$
[20] Stefan Ulbrich. Generalized SQP methods with “parareal” time-domain decomposition for time-dependent pde-constrained optimization. In *Real-time PDE-constrained optimization*, pages 145–168. SIAM, 2007.
where $\mathrm{sinc}(x) = \sin(x)/x$ is the sinc function and $\mathcal{A}^{1/2}$ is the symmetric square root of $\mathcal{A}$ (i.e., a symmetric matrix with positive eigenvalues whose square is $\mathcal{A}$). This leads to

$$e^{s\mathcal{L}}\mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} = [\cos(s\mathcal{A}^{1/2}) + s\mathcal{L}\mathrm{sinc}(s\mathcal{A}^{1/2})] \cdot \mathcal{B} \cdot \mathcal{B}^T \cdot [\cos(s\mathcal{A}^{1/2}) + s \mathrm{sinc}(s\mathcal{A}^{1/2})\mathcal{L}^T].$$

By direct calculation and using $x\mathrm{sinc}(xA^{1/2}) = A^{-1/2}\sin(xA^{1/2})$, we obtain

$$e^{s\mathcal{L}}\mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} = \begin{bmatrix} C_{11}(s) & C_{12}(s) \\ C_{21}(s) & C_{22}(s) \end{bmatrix},$$

where

$$C_{11}(s) = A^{-1} \sin^2(sA^{1/2}),$$

$$C_{12}(s) = A^{-1/2} \cos(sA^{1/2}) \sin(sA^{1/2}),$$

$$C_{21}(s) = C_{12}(s),$$

$$C_{22}(s) = \cos^2(sA^{1/2}).$$

We can now integrate these coefficients to obtain

$$I + \frac{1}{\alpha} \int_{0}^{T} e^{s\mathcal{L}}\mathcal{B}\mathcal{B}^T e^{s\mathcal{L}^T} ds = \mathcal{M} = \begin{bmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{bmatrix},$$

where

$$M_{11} = I + \frac{T}{2\alpha}A^{-1} - \frac{1}{4\alpha}A^{-3/2}\sin(2TA^{1/2}),$$

$$M_{12} = \frac{1}{4\alpha}A^{-1}(I - \cos(2TA^{1/2})),$$

$$M_{21} = M_{12},$$

$$M_{22} = \left(1 + \frac{T}{2\alpha}\right)I + \frac{1}{4\alpha}A^{-1/2}\sin(2TA^{1/2}).$$

To use $\mathcal{M}$ as a preconditioner, we must calculate $\mathcal{M}^{-1}$. Since all the $M_{ij}$ are square and they all commute, we have

$$\mathcal{M}^{-1} = \begin{bmatrix} D^{-1}M_{22} & -D^{-1}M_{12} \\ -D^{-1}M_{21} & D^{-1}M_{11} \end{bmatrix},$$

with $D = M_{11}M_{22} - M_{12}M_{21}$. Of course we do not want to compute square roots and trigonometric functions of $\mathcal{A}$; to approximate the blocks of $\mathcal{M}^{-1}$, we again study their symbols, i.e., the functions obtained by replacing $\mathcal{A}$ by a scalar $\sigma > 0$ (this is equivalent to diagonalizing the blocks and studying their eigenvalues as a function of the eigenvalues of $\mathcal{A}$). Since $\sigma$ is an approximation of $k^2\pi^2$, $k = 1, 2, \dots$, we observe that for highest frequencies $\sigma$,

$$\bullet D \approx \frac{1}{4\alpha^2}(2\alpha + T)A^{-1}(TI + 2\alpha A);$$

$$\bullet D^{-1}M_{11} \approx \left[1 - \frac{1}{(2\alpha + T)}\right] I;$$

$$\bullet D^{-1}M_{12} = D^{-1}M_{21} \approx 0;$$

$$\bullet D^{-1}M_{22} \approx I - (TI + 2\alpha A)^{-1}.$$

These observations lead us to approximate $\mathcal{M}^{-1}$ by a block diagonal matrix $\widehat{\mathcal{M}}^{-1}$ of the form

$$\widehat{\mathcal{M}}^{-1} = I - \begin{bmatrix} (aI + bA)^{-1} & 0 \\ 0 & cI \end{bmatrix}, \qquad (24)$$

where $a = T, b = 2\alpha$ and finally $c = 1/(2\alpha + T)$. Each application of the preconditioner therefore involves only one solution of the elliptic problem $(aI + bA)v = f$, just like for the heat equation.

In our numerical experiment, we take $T = 2.3$, $u_0(x) = \exp(-100(x - 1/2)^2)$, so that $y_{in} = (u_0(x), 0)$ and $y_{tg}(x) = [(\frac{1}{2}\exp(-100(x - 1/4)^2) + \frac{1}{2}(-100(x - 3/4)^2), 0]$. We set $L = 10, r = 100$, and we use the same Matlab functions for computing $\mathcal{R}_\ell$ and $\mathcal{P}_\ell$.

Let us investigate the convergence speed of the linear system (12) in GMRES for $N = 1000$. We use for this instance the function gmres in Matlab with restart=(), tol=1e-8 and maxit=size(M,1). The results are shown in Table
<table><caption>Table 3: GMRES for solving (12) involved (23) for $\alpha = 10^{-6}$.</caption><thead><tr><th colspan="4">Matrix M</th></tr><tr><th>Schemes</th><th>cond(M)</th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>3.8e4</td><td>60</td><td>9.85e-10</td></tr><tr><td>SDIRK</td><td>3.8e4</td><td>60</td><td>9.85e-10</td></tr></tbody><thead><tr><th colspan="4">Preconditioned matrix MP<sup>-1</sup></th></tr><tr><th>Schemes</th><th>cond(MP<sup>-1</sup>)</th><th>Nbr Iters</th><th>Residuals</th></tr></thead><tbody><tr><td>Euler</td><td>2.59</td><td>4</td><td>3.33e-14</td></tr><tr><td>SDIRK</td><td>2.59</td><td>4</td><td>3.31e-14</td></tr></tbody></table>

3, where we observe that (12) converges in 84 iterations for both quadrature formulas and $\alpha = 10^{-6}$ whereas the preconditioned one converges in only 4 iterations.

Unlike for the heat equation, the high order of SDIRK does not improve the convergence speed for preconditioned and unpreconditioned system compare to Implicit Euler. Nonetheless, even for $N \in \{100, ..., 1000\}$ the number of iterations in GMRES for the preconditioned system remains the same whether for SDIRK or Implicit Euler.

These results show that for fixed $r$, the discretization parameters in time do not really change the behaviour of the convergence of (12) in GMRES. We also observe that when $r$ increases, the number of iterations of the unpreconditioned

<table><caption>Table 4: GMRES for solving (12) involved (23) with Implicit Euler with $\alpha = 10^{-6}$ for various $r$.</caption><thead><tr><th>r</th><th>cond(M)</th><th>Nbr Iters</th><th>cond(MP<sup>-1</sup>)</th><th>Nbr Iters</th></tr></thead><tbody><tr><td>10</td><td>4.83e2</td><td>10</td><td>2.47</td><td>5</td></tr><tr><td>50</td><td>1.03e4</td><td>41</td><td>2.59</td><td>4</td></tr><tr><td>150</td><td>7.75e4</td><td>76</td><td>2.81</td><td>3</td></tr><tr><td>250</td><td>1.68e5</td><td>92</td><td>2.53</td><td>3</td></tr><tr><td>350</td><td>2.48e5</td><td>104</td><td>2.49</td><td>3</td></tr></tbody></table>

system increases whereas the number of iterations of the preconditioned system decreases. The results are shown in Table 4.

<table><caption>Table 5: GMRES for solving (12) involved (23) with Implicit Euler for various $\alpha$.</caption><thead><tr><th>α</th><th>cond(M)</th><th>Nbr Iters</th><th>cond(MP<sup>-1</sup>)</th><th>Nbr Iters</th></tr></thead><tbody><tr><td>1e-5</td><td>2.26e4</td><td>52</td><td>2.41</td><td>3</td></tr><tr><td>1e-4</td><td>4.48e3</td><td>34</td><td>1.64</td><td>3</td></tr><tr><td>1e-3</td><td>4.98e2</td><td>20</td><td>1.09</td><td>3</td></tr><tr><td>1e-2</td><td>5.12e1</td><td>13</td><td>1.01</td><td>3</td></tr><tr><td>1e-1</td><td>6.03</td><td>9</td><td>1.01</td><td>3</td></tr><tr><td>1e0</td><td>1.5</td><td>6</td><td>1.0</td><td>3</td></tr><tr><td>1e1</td><td>1.05</td><td>4</td><td>1.0</td><td>2</td></tr></tbody></table>

Finally, in Table 5, we observe that when $\alpha$ increases, our preconditioner remains efficient.

# 5 Conclusion

We introduced a new time parallel algorithm for time dependent linear quadratic optimal control problem when a cheap exponential integrator is available. We proposed two preconditioners to solve efficiently the linear equation that comes from the algorithm on the particular cases of heat and wave equations respectively.

We are currently working the behavior of the preconditioner when an explicit method is used to obtain $M$, while respecting the CFL condition to better understand the impact of the use of explicit and implicit methods for solving (23). Finally, we are studying the optimal control problem involving the wave equation with boundary control in one and several dimensions.
# Acknowledgements

The work described in this paper is partially supported by a grant from the ANR/RGC joint research scheme sponsored by the Research Grants Council of the Hong Kong Special Administrative Region, China (Project no. A-CityU203/19), by the French National Research Agency (Project ALLOWAPP, grant ANR-19-CE46-0013-01 and Projet DEEPNUM, grant ANR-21-CE23-0017) and by Mitacs-Globalink Canada and INRIA France (grant IT34057, 2023).

# References

[1] Andrew T Barker and Martin Stoll. Domain decomposition in time for pde-constrained optimization. *Computer Physics Communications*, 197:136-143, 2015.

[2] Albrecht Böttcher and Sergei M. Grudsky. *Spectral Properties of Banded Toeplitz Matrices*. Society for Industrial and Applied Mathematics, 2005.

[3] M. Gander, F. Kwok, and G. Wanner. Constrained optimization: From Lagrangian mechanics to optimal control and PDE constraints. *in Optimization with PDE Constraints*, Springer, pages 151-202, 2014.

[4] M. J. Gander, F. Kwok, and J. Salomon. ParaOpt: A parareal algorithm for optimality systems. *SIAM Journal on Scientific Computing*, 42(5):A2773-A2802, 2020.

[5] Martin J Gander. 50 years of time parallel time integration. In *Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD*, Heidelberg, May 6-8, 2013, pages 69-113. Springer, 2015.

[6] Martin J. Gander and Stefan Güttel. Paraexp: A parallel integrator for linear initial-value problems. *SIAM Journal on Scientific Computing*, 35(2):C123-C142, 2013.

[7] Martin J Gander and Laurence Halpern. Optimized schwarz waveform relaxation methods for advection reaction diffusion problems. *SIAM Journal on Numerical Analysis*, 45(2):666-697, 2007.

[8] Martin J Gander and Felix Kwok. Schwarz methods for the time-parallel solution of parabolic control problems. In *Domain decomposition methods in science and engineering XXII*, pages 207-216. Springer, 2016.

[9] G.H. Golub and C.F. Van Loan. *Matrix computations*. The John Hopkins University Press, Baltimore, 2013.

[10] Matthias Heinkenschloss. A time-domain decomposition iterative method for the solution of distributed linear quadratic optimal control problems. *Journal of Computational and Applied Mathematics*, 173(1):169-198, 2005.

[11] Nicholas J. Higham. *Functions of Matrices*. Society for Industrial and Applied Mathematics, 2008.

[12] Thi-Thao-Phuong Hoang, Jérôme Jaffré, Caroline Japhet, Michel Kern, and Jean E Roberts. Space-time domain decomposition methods for diffusion problems in mixed formulations. *SIAM Journal on Numerical Analysis*, 51(6):3532-3559, 2013.

[13] Shishun Li and Xiao-Chuan Cai. Convergence analysis of two-level space-time additive schwarz method for parabolic equations. *SIAM Journal on Numerical Analysis*, 53(6):2727-2751, 2015.

[14] Yvon Maday, Julien Salomon, and Gabriel Turinici. Monotonic parareal control for quantum systems. *SIAM Journal on Numerical Analysis*, 45(6):2468-2482, 2007.

[15] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. *SIAM Review*, 45(1):3-49, 2003.

[16] Benjamin W Ong and Jacob B Schroder. Applications of time parallelization. *Computing and Visualization in Science*, 23:1-15, 2020.

[17] Andrea Toselli and Olof Widlund. *Domain Decomposition Methods – Algorithms and Theory*. Springer-Verlag, Berlin, 2005.

[18] Colin S.C. Tsang. Preconditioners for linear parabolic optimal control problems. Master's thesis, Hong Kong Baptist University, Hong Kong, China, August 2017.

[19] Stefan Ulbrich. Generalized sqp methods with “parareal” time-domain decomposition for time-dependent pde-constrained optimization. In *Real-time PDE-constrained optimization*, pages 145-168. SIAM, 2007.