# Time-parallel data assimilation using a Paraexp-based observer strategy

Lucas Perrin, Julien Salomon

► To cite this version:

Lucas Perrin, Julien Salomon. Time-parallel data assimilation using a Paraexp-based observer strategy. 2026. hal-05522254

HAL Id: hal-05522254

https://hal.science/hal-05522254v1

Preprint submitted on 22 Feb 2026

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# TIME-PARALLEL DATA ASSIMILATION USING A PARAEXP-BASED OBSERVER STRATEGY

LUCAS PERRIN* AND JULIEN SALOMON†

**Abstract.** In this paper, we propose a time-parallel strategy for data assimilation based on the coupling of the Paraexp algorithm with a Luenberger observer. The observer dynamics, being linear, can be reformulated as a linear initial-value problem, making it particularly suitable for the Paraexp approach. To address assimilation over an unbounded time interval, we introduce a windowing strategy together with an adaptive refinement of the time-parallel accuracy across windows. We show that, under a suitable accuracy condition, the exponential convergence rate of the classical observer is preserved. The proposed method is illustrated on several linear systems arising from discretized partial differential equations, including hyperbolic models, for which classical iterative parallel-in-time methods may encounter difficulties. Numerical experiments demonstrate significant efficiency gains and, in some cases, superlinear speedups due to the adaptive accuracy strategy.

**Key words.** Data assimilation, Luenberger observer, Paraexp algorithm, parallel-in-time methods, exponential convergence, hyperbolic systems, matrix exponential.

**AMS subject classifications.** Primary 93B07, Secondary 65M22, 65F60, 93C20

## 1. Introduction.
The efficient solving of large-scale dynamical systems is a central challenge, particularly in the context of real-time applications such as weather forecasting, fluid dynamics, and other high-dimensional systems. These systems often require robust techniques for both integrating over time and assimilating data. This paper focuses on combining two distinct but complementary approaches: the Paraexp algorithm for time-parallelization and the Luenberger observer for data assimilation.

The Luenberger observer [29] is a classical tool in control theory, designed to estimate the state of a system based on partial measurements. It operates by applying a corrective feedback that is linearly proportional to the error between observed and predicted outputs. This correction stabilizes the system, ensuring the convergence of the estimated state to the true state. The observer is especially useful for systems modelled by partial differential equations (PDEs), in the case where complete knowledge of the initial state is unavailable. The convergence of the observer is highly dependent on the observability of the system and the appropriate choice of the gain matrix, $L$, which ensures that the error decays to zero over time. However, when dealing with high-dimensional systems, the computational cost of integrating such systems over long periods can become prohibitive, requiring more efficient time-stepping methods.

On the other hand, the Paraexp algorithm [11] offers a time-parallel approach specifically designed for linear systems, such as those governed by the Luenberger observer. Unlike the widely-used Parareal method, which can face limitations when applied to hyperbolic systems, Paraexp leverages the linearity of the problem to directly decompose the integration task across multiple processors, theoretically achieving near-optimal efficiency. The method capitalizes on the matrix exponential's exact solution to linear systems, allowing for highly efficient parallel computation, particularly suited for hyperbolic systems where other methods struggle [13, 10]. In contrast to iterative methods like Parareal, Paraexp is a direct method, often requiring fewer

<footnote>*Department of Mathematics and Statistics, University of Konstanz, 78464 Konstanz, Germany previously at INRIA Paris, ANGE Project-Team, 75589 Paris Cedex 12, France (lucas.perrin@uni-konstanz.de).</footnote>

<footnote>†INRIA Paris, ANGE Project-Team, 75589 Paris Cedex 12, France and Sorbonne Université, CNRS, Laboratoire Jacques-Louis Lions, 75005 Paris, France</footnote>
iterations and yielding faster convergence.

The combination of these two methods forms the core of this article. By coupling the Luenberger observer's state estimation framework with the time-parallel capabilities of Paraexp, we propose a strategy that not only improves the speed of data assimilation processes, but also enhances their scalability. This combined approach addresses both the computational demands of time integration and the accuracy requirements of state estimation in real-time applications.

Specifically, the application of the Paraexp algorithm to observer dynamics allows efficient parallelization, overcoming the limitations imposed by serial time-stepping methods.

In this work, we focus on two main aspects. First, we adapt the Paraexp time-parallel algorithm to the dynamics of a Luenberger observer over an unbounded time horizon by introducing a windowing strategy. Second, we propose an adaptive accuracy refinement across time windows that ensures preservation of the observer's exponential convergence rate while reducing the overall computational cost. These ideas are illustrated in several linear systems arising from discretized partial differential equations.

This paper builds on the work of [34], which proposes a coupling between data assimilation and time-parallel computing based on the Parareal algorithm. Here, we investigate an alternative strategy based on Paraexp, which is well suited for linear observer dynamics and for hyperbolic-type problems, where iterative parallel-in-time methods can encounter difficulties.

We begin by revisiting the theory behind the Luenberger observer and its role in data assimilation in Section 2, followed by a discussion of the Paraexp algorithm and its advantages in solving linear initial-value problems in Section 3. We then present the proposed coupling strategy and discuss its convergence properties in Section 4, before demonstrating its effectiveness through numerical experiments on hyperbolic systems in Section 5.

**2. Data assimilation.** We introduce the Luenberger observer, the data assimilation technique that is our focus. Throughout the paper, continuous-time quantities are denoted by standard symbols, while their time-discrete approximations are introduced explicitly when needed.

**2.1. The Luenberger observer setting.** Based on a simple concept, the Luenberger observer is widely used, as least as a first approximation [28]. It consists of applying a penalty on the dynamical system that is linearly proportional to the difference between what is observed and what is expected. We therefore consider a situation where only partial information is known about the state of the system under study with the following system:

$$ (2.1) \qquad \begin{cases} \frac{d}{dt}x(t) = Ax(t) + Bu(t), & x(0) = x_0, \quad t \ge 0, \\ y(t) = Cx(t), & \end{cases} $$

where $A \in \mathcal{M}_{m \times m}(\mathbb{R})$, the matrix of the corresponding linear system, $B \in \mathcal{M}_{m \times p}(\mathbb{R})$, the matrix associated with the control, and $C \in \mathcal{M}_{q \times m}(\mathbb{R})$, the matrix associated with the observation, are assumed to be known. We also consider $p, q, m \in \mathbb{N}^*$, with $p < m$ and generally $q \ll m$. The vector $x(t) \in \mathbb{R}^m$ is called the state vector, $y(t) \in \mathbb{R}^q$ is the measured output, $u(t) \in \mathbb{R}^p$ the control, and $t \in (0, +\infty)$. The initial condition $x(0) = x_0$ is unknown.

Throughout this work, we restrict ourselves to linear time-invariant systems with constant matrices A, B and C. This assumption is essential for the application of
the Paraexp algorithm and will be maintained in all subsequent developments. As a result, we do not know the trajectory $x$ associated with this unknown initial condition. We then want to find $x$ with only knowledge of its dynamics, i.e. the matrices $A$ and $B$ as well as the vector $u$, and the observations $y$. To tackle this problem, Luenberger proposed an observer $\hat{x}(t)$ of $x(t)$ [29], which follows the same dynamics, modified by a linear penalty term on the known variables, i.e., $y$ [36, 26, 9, 5]. In this way, the observer is defined by:

$$ (2.2) \qquad \begin{cases} \frac{d}{dt}\hat{x}(t) = A\hat{x}(t) + Bu(t) + L[y(t) - \hat{y}(t)], & \hat{x}(0) = \hat{x}_0, \quad t \ge 0, \\ \hat{y}(t) = C\hat{x}(t), & \end{cases} $$

with $L \in \mathcal{M}_{m \times q}(\mathbb{R})$ the observer gain, and $\hat{x}(0)$ an arbitrary initial condition. This system can be written as:

$$ (2.3) \qquad \frac{d}{dt}\hat{x}(t) = M\hat{x}(t) + g(t), \quad \hat{x}(0) = \hat{x}_0, \quad t \ge 0, $$

with:

$$ (2.4) \qquad \begin{cases} M = (A - LC), \\ g(t) = Bu(t) + Ly(t). \end{cases} $$

In this construction, the matrix $L$ needs to be specified, and, as pointed out in [34], it plays a crucial role in the estimation of the error $\varepsilon(t) := x(t) - \hat{x}(t)$ as we have:

$$ (2.5) \qquad \varepsilon(t) = e^{(A-LC)t}\varepsilon(0), $$

and therefore, if the eigenvalues of $A - LC$ sit in the left half complex plane (i.e., $(A - LC)$ is Hurwitz), so that the error decays to 0 as $t \to +\infty$.

Let us now detail how to design and compute $L$ in practice to ensure convergence. Such a property relies on classical results from observability and control theory [7]. We briefly recall these notions here for completeness and to fix notation; detailed proofs and further discussions can be found in the cited references. The first definition is the one of a controllable system.

**DEFINITION 2.1.** A system is said to be completely controllable if for any $t_0$, any initial state $x(t_0) = x_0$ and any given a state $x_f$ there exists a finite time $t_1 > t_0$ and a control $u(t)$, $t_0 < t < t_1$ such that $x(t_1) = x_f$.

The well-known Kalman's rank condition [23] provides a criterion to guarantee controllability.

**THEOREM 2.2.** A system is completely controllable if and only if one of the following equivalent conditions holds:

1. rank $[B, AB, \dots, A^{n-1}B] = n$;

2. rank $[B, A - \lambda I] = n$, ($\forall \lambda \in \mathbb{C}$);

3. $\{s^T A = \mu s^T \text{ and } s^T B = 0\} \iff s^T = 0$.

Luenberger theory also makes use of observability condition which is defined as follows.

**DEFINITION 2.3.** A system is completely observable if for any $t_0$ and any initial state $x(t_0) = x_0$, there exists a finite time $t_1 > t_0$ such that knowledge of $u(t)$ and $y(t)$ for $t_0 \le t \le t_1$ suffices to determine $x_0$ uniquely. There is no loss of generality in assuming $u(t) \equiv 0$ throughout the interval.
126 Analogously to the case of controllability, a rigorous criterion ensuring observability
can be stated.
127

128
129
THEOREM 2.4. A system is said to be completely observable if and only if one of the following equivalent conditions holds:

$$
\begin{align*}
1. \text{ rank} \quad & \begin{bmatrix} C \\ CA \\ \vdots \\ CA^{n-1} \end{bmatrix} &&= n; \\
2. \text{ rank} \quad & \begin{bmatrix} A - \lambda I \\ C \end{bmatrix} &&= n \quad (\forall \lambda \in \mathbb{C}); \\
3. \{As = \mu s \quad \text{and} \quad Cs = 0\} &&\iff \\
& s = 0. &&
\end{align*}
$$

131 To put it simply, controllability refers to the ability to bring a system to any desired state, while observability refers to the ability to determine the initial state from 132 subsequent states. The following theorem from [29, p.303] links observability and 133 controllability in our case.
134

135
136
THEOREM 2.5. *Given a completely observable system ((2.1)), an observer of the form (2.2) exists. Moreover, the eigenvalues of A − LC can be selected arbitrarily.*

137 In other words, if our system is observable, a control of the form +L[y − ŷ] can bring
138 our observer ̂x to the desired state x.
139

139
140
141
Given an arbitrary set of eigenvalues of (A − LC), the actual computation of
the corresponding matrix L can be found in [21] and is theoretically justified by 2.5.
Nevertheless, we recall [34, Proposition 2.2].

142
143
PROPOSITION 2.6. Suppose that (2.1) is observable and that the eigenvalues of A − LC are negative and simple. Then, we have:

$$
\frac{\| e^{(A-LC)t} \|}{\| e^{(A-LC)t} \|} \le \eta e^{-\mu t},
$$

145
146
with $\mu := \min_{\nu \in \sigma(A-LC)} |\nu|$, and $\eta := \text{cond}(V) = \|V\| \|V^{-1}\|$ where $V$ is the matrix whose rows are the left-eigenvectors of $A-LC$.

147
The following estimate results from this proposition and (2.5):

$$
\|\epsilon(t)\| \le \eta \|x(0) - \hat{x}(0)\| e^{-\mu t}.
$$

149
150
As we do not know $\|x(0) - \hat{x}(0)\|$, only the convergence rate $\mu$ and $\eta$ are determined
by the choice of $\sigma(A - LC)$.

151
152
153
Finally, let us remark that (2.3)-(2.4) is a linear initial-value problem which can
be solved using a parallel scheme. This is the core idea of this article. We explore
this in Section 4.

154
155
156
157
158
159
160
**2.2. Data assimilation, discretized PDEs, and observability.** We state
here some remarks about observability and controllability in discretized PDEs. While
observer theory is well understood at the continuous level, its application to systems
governed by PDEs raises additional difficulties once space discretization is introduced.
In particular, classical convergence analyses based solely on the spectral properties of
the observer dynamics operator are generally no longer sufficient. Indeed, it is now
well known that systems that are observable at the continuous level may lose this
property after discretization, due, for instance, to spurious high-frequency numerical modes; see, e.g., [37, 20, 33]. As a consequence, classical observer designs based on pole placement do not necessarily yield convergent discrete observers without additional numerical stabilization.

Moreover, the underlying properties of the PDEs give us certain constraints on their observability. Indeed, the heat equation, for instance, is, in 1D and 2D, observable in the sense of Definition 2.3 whatever the observable domain. Whereas, for instance, in the case of the 2D wave equation, the observable domain region $\omega$ needs to catch all the bouncing rays that the border reflects. This leads to a constraint requiring part of the boundary to be observed in order to capture all geometric rays, as characterized by the geometric control condition [2].

In such cases, Theorem 2.5 does not apply, and prevent us arbitrarily choosing the eigenvalues of $A - LC$. We maintain a general analysis of the problem until offering some examples, and for each examples we carefully motivate our choice for the gain matrix $L$. We also explain the construction of the matrix $C$ for each examples.

For the wave equation, [17, 16] reports that a simple matrix $L$ of the form: $L = \gamma C^T$, where the gain, $\gamma > 0$, is a tunable parameter, suffices to have a convergent observer. This follows from [27][theorem 2.3]. This is of course in the case of an observable setting where Bardos-Lebeau condition [2] is met.

For the wave equation, we consider the case where only a part of the velocity field is observed as it is presented in [6]. We justify not using an observation on the position as [6] show us that implementing *Schur Displacement Feedback* yields similar results. Finally, as it is known that spatial discretization induces undesirable numerical artifacts in the stabilization and control of a wave equation, we choose to use the scheme proposed by A. Münch [32], together with a numerical crime, in order to mitigate spurious high-frequency effects induced by spatial discretization.

### 3. Time parallelization.

**3.1. The Paraexp algorithm.** In [34], Kwok, Riffo and Salomon present a time parallelization of the *Luenberger observer* using the well known *Parareal* [30] algorithm. We are here interested in another time-parallel scheme, *Paraexp* [11], which is a direct method based on the decomposition of linear initial-value problems into homogeneous and non-homogeneous subproblems. Two key points motivate the choice to extend [34] work to the *Paraexp* algorithm. First, one of the drawbacks of the *Parareal* scheme is that it is a multiple shooting method and therefore reaches limits when tackling hyperbolic problems, as shown in [13, 10]. Second, the *Luenberger observer dynamics* is linear. The *Paraexp* algorithm is designed to perform extremely well on such problems. Moreover, in a linear setting, for *Parareal* at least two iterations are required, whereas *Paraexp* is a direct method. This means that the efficiency of *Parareal* is at most 50% whereas the one of *Paraexp* is at most 100%. Hence, most of our numerical examples in Section 5 are applications to hyperbolic problems.

Let us explain the *Paraexp* approach in the case of (2.3). This formulation includes, e.g., space-discretization of partial differential equations (via, e.g., finite difference schemes, or more complex ones like finite volumes or finite element schemes).

The time parallelization relies on two facts. First, the exact solution of the homogeneous initial value problem

$$ \begin{cases} \frac{d}{dt}w(t) = Mw(t), & t \in [0, T], \\ w(0) = x_0, \end{cases} $$
208 is given by the matrix exponential:

209 $$w(t) = \exp(tM)x_0,$$

210 where $\exp(tM) = (I + (tM) + (tM)^2/2 + ...).$ Second, if we integrate the non-homogeneous problem with a null initial value

212 $$ \begin{cases} \frac{d}{dt}v(t) = Mv(t) + g(t), & t \in [0, T], \\ v(0) = 0_{\mathbb{C}^m}. \end{cases} $$

213 The solution of (2.3) is given by:

214 $$x(t) = v(t) + w(t).$$

215 In order to take advantage of N processors for the time parallelization, we consider a decomposition of $[0, T]$ of the form:

217 $$[0, T] = \bigcup_{j} [t_{j-1}, t_j], \quad j = 1, \dots, N, \quad 0 = t_0 \le t_1 \le \dots \le t_N = T,$$

218 and the following subproblems:

219 Type 1: For $j = 1, \dots, N$ solve:

220 $$ \frac{d}{dt}v_j(t) = Mv_j(t) + g(t), \quad v_j(t_{j-1}) = 0_{\mathbb{C}^m}, \quad t \in [t_{j-1}, t_j], $$

221 using some serial integrator.

222 Type 2: For $j = 1, \dots, N$ solve:

223 $$ \frac{d}{dt}w_j(t) = Mw_j(t), \quad w_j(t_{j-1}) = v_{j-1}(t_{j-1}), \quad t \in [t_{j-1}, T], $$

224 setting $v_0(t_0) := x_0$, using an exponential propagator.

225 Type 1 subproblems are independent from each other. as well as the subproblems of Type 2 once $v_{j-1}(t_{j-1})$ is known:

227 $$ w_j(t) = \exp((t - t_{j-1})M)v_{j-1}(t_{j-1}). $$

228 In 3.2 we recall the concept of exponential propagator, as used in [11].

229 The solution is:

230 $$ x(t) = v_k(t) + \sum_{j=1}^{k} w_j(t) \quad \text{with } k \text{ such that } t \in [t_{k-1}, t_k] $$

231 The parallel procedure is given by simultaneously integrating the subproblems of Type

232 1 and Type 2, and finally forming the sum for $x(t)$ at the required time points $t_k$.

233 It is important to note that by choosing

234 $$ \begin{cases} x^S(t) = \sum_{j=1}^{N} \tilde{v}_j(t), \\ x^H(t) = \sum_{j=1}^{N} \tilde{w}_j(t), \end{cases} $$
where $\tilde{v}_j$ is defined on $[0, T]$, it coincides with $v_j(t)$ on $[t_{j-1}, t_j]$ and cancels everywhere else. The same property holds for $\tilde{w}_j$ in the sense that $\tilde{w}_j$ coincides with $\sum_{j=1}^k w_j$ on $[t_{k-1}, t_k]$, and cancels everywhere else. As a consequence:

$$x(t) = x^S(t) + x^H(t).$$

<figure><img src="image_2.png" /><figcaption>FIGURE 3.1. Overlapping time decomposition of an initial-value problem into four non-homogeneous problems with zero initial guess (Type 1, solid red curves) and four homogeneous problems (Type 2, dashed blue curves). Type 2 sub-problems are exponentially propagated. The solution of the original problem is obtained by summation of all these curves. (from [11]).</figcaption></figure>

## 3.2. Matrix exponential computation.
In this Section, we present the concept of exponential propagator as introduced by Gander and Güttel [11]. Recall that solving the homogeneous equation:

$$\begin{cases} \frac{d}{dt}x(t) = Mx(t), \\ x(0) = x_0, \end{cases}$$

yields the solution $x(t) = \exp(tM)x_0$, which yields a matrix exponential computation.

In the original Paraexp paper, Gander and Güttel suggest two methods to obtain a "near-optimal approximation of the matrix exponential", namely, projection based methods, which use rational Krylov methods, and expansion-based methods, which use Chebychev polynomials. The former is preferred when addressing parabolic problems (e.g. Heat equation), the latter when addressing hyperbolic problems (e.g. Wave equation). These choices are mostly due to the eigenstructure of the involved matrices. Indeed, the choice of a single real pole of high multiplicity for the rational Krylov method in [11] is relevant for matrices with real eigenvalues, which is the case for parabolic problems, whereas interpolation by Chebyshev polynomials is more relevant for matrices with eigenvalues on the imaginary axis, which is the case for hyperbolic problems. However, these two methods are not applicable to the case of a Luenberger observer when the observed state follows a hyperbolic dynamics, which is precisely the case we are interested in Section 5.

Indeed, although the state follows a hyperbolic trajectory and the matrix of its dynamics (of the semi-discretized system in space) has purely imaginary eigenvalues, the corresponding observer does not, with different and non-pure imaginary eigenvalues.

It is fairly easy to intuit this, for instance by considering that the observer's standard should decrease at a constant rate, which is not the case for hyperbolic
systems that conserve mass.

We propose two ways to tackle this issue. First, a robust method for computing the exponential of any matrix can be obtained by the means of Padé approximations coupled with the scaling and squaring method [18]. This approach is implemented in Matlab for the command expm. The main drawback of this strategy lies in its computational cost when tackling large problems. Second, one can use a rational Krylov method with optimally chosen poles in order to have an accurate approximation in a Krylov space of small dimension [3, 4]. The drawback of this solution is that it requires to study the matrix of the observer dynamics beforehand which can also reveal to be costly if the problem is of large size.

**4. Coupling Paraexp with the Luenberger observer.** We now present a scheme to parallelize in time an observer $\hat{x}$ of a state $x$.

**4.1. Time parallelisation on unbounded time interval.** In order to accelerate the assimilation and take advantage of a time-parallelization procedure, we divide the time interval $[0, +\infty)$ into windows of a given length $T > 0$ denoted by:

$$W_{\ell} := [T_{\ell-1}, T_{\ell}], \quad \ell \in \mathbb{N}^{*},$$

where $T_{\ell} = \ell \cdot T$. We then introduce a parallel version of (2.2) on each window $W_{\ell}$:

$$ (4.1) \qquad \begin{cases} \frac{d}{dt} \hat{x}_{\parallel}(t) = A \hat{x}_{\parallel}(t) + Bu(t) + L[y(t) - C \hat{x}_{\parallel}(t)], \\ \hat{x}_{\parallel}(T_{\ell-1}) = \hat{\mathbf{X}}_{\ell-1}^{h_{\ell-1}}, \end{cases} \quad t \in [T_{\ell-1}, T_{\ell}]. $$

where $\hat{x}_{\parallel}(T_{\ell}) = \hat{\mathbf{X}}_{\ell}^{h_{\ell}}$ is the approximation of $\hat{x}(T_{\ell})$ obtained by the parallel solver over the window $\ell$ with $h_{\ell}$ a auxiliary variable which stands for the desired precision of our parallel solver over this window. We adopt the convention that $\hat{x}_{\parallel}(0) = \hat{\mathbf{X}}_{0}^{h_{0}} = \hat{x}_{0}$.

Since the reconstruction error is potentially significant at the beginning of the observation, it is unnecessary to solve (4.1) with very high precision in the first windows. We then start with a rough approximation, and refine it from one window to another in such a way that the rate of convergence of the original (non-parallelized) observer is preserved. Hence, we will apply Paraexp with increasingly higher precision as the windows are processed.

**4.2. Convergence of the parallel observer.** We analyze the convergence of the time-parallel observer and identify the origin of the additional error introduced by the parallel-in-time procedure.

Let $\hat{x}$ denote the solution of the continuous observer (2.2), and let $\hat{x}_{\parallel}$ denote the approximation produced by the parallel solver over a given window $W_{\ell} = [T_{\ell-1}, T_{\ell}]$.

Both $\hat{x}$ and $\hat{x}_{\parallel}$ satisfy linear dynamics of the form

$$ \frac{d}{dt} z = Mz + g(t), $$

with identical data, but $\hat{x}_{\parallel}$ is affected by the numerical error induced by the time-parallel integration.

Throughout this analysis, we assume that the homogeneous subproblems (Type 2) are solved either exactly or with an accuracy that is negligible compared to the discretization error of the non-homogeneous subproblems. This assumption is satisfied, for instance, when using exact diagonalization, scaling-and-squaring with Padé approximants, or sufficiently accurate Krylov-based exponential propagators.
We denote the parallel error as:

$$ \varepsilon_{\parallel}(t) := x(t) - \hat{x}_{\parallel}(t), $$

and the Paraexp decompositions at the end of window $W_{\ell}$:

$$ \hat{x}_{\parallel}(T_{\ell}) = \hat{x}_{\parallel}^{S}(T_{\ell}) + \hat{x}_{\parallel}^{H}(T_{\ell}), $$

where the superscripts S and H denote the non-homogeneous (Type 1) and homogeneous (Type 2) contributions, respectively.

Since the homogeneous subproblems are solved exactly up to the accuracy of the exponential propagator, the difference between $\hat{x}$ and $\hat{x}_{\parallel}$ originates solely from the numerical approximation of the Type 1 subproblems. Consequently, we obtain the decomposition:

$$ \varepsilon_{\parallel}(T_{\ell}) = \varepsilon(T_{\ell}) + (\hat{x}^{S}(T_{\ell}) - \hat{x}_{\parallel}^{S}(T_{\ell})), $$

and therefore the estimate:

$$ (4.2) \qquad \|\varepsilon_{\parallel}(T_{\ell})\| \le \|\varepsilon(T_{\ell})\| + \|\hat{x}^{S}(T_{\ell}) - \hat{x}_{\parallel}^{S}(T_{\ell})\|. $$

The first term on the right-hand side corresponds to the classical observer error, which converges exponentially under the assumptions of Proposition 2.6. The second term represents the numerical error induced by the parallel-in-time discretization of the Type 1 subproblems. The convergence of the parallel observer therefore follows provided this additional error is controlled appropriately.

**4.3. Preservation of the convergence rate.** We now derive a sufficient condition ensuring that the parallel observer preserves the exponential convergence rate of the classical observer.

From Proposition 2.6, the serial observer satisfies

$$ \|\varepsilon(T_{\ell})\| \le \eta \|\varepsilon(0)\| e^{-\mu T_{\ell}}, \quad \ell \ge 1. $$

Combining this estimate with (4.2), it is sufficient to require that the numerical error introduced by the Type 1 subproblems decays at least at the same rate, namely

$$ (4.3) \qquad \|\hat{x}_{\parallel}^{S}(T_{\ell+1}) - \hat{x}^{S}(T_{\ell+1})\| \le e^{-\mu T} \|\hat{x}_{\parallel}^{S}(T_{\ell}) - \hat{x}^{S}(T_{\ell})\|, \quad \ell \ge 1. $$

Assume that the Type 1 subproblems are discretized by a Runge-Kutta scheme of order k with timestep $h_{\ell}$ on window $W_{\ell}$. Then there exists a constant $C > 0$, independent of $\ell$, such that

$$ \|\hat{x}_{\parallel}^{S}(T_{\ell}) - \hat{x}^{S}(T_{\ell})\| \le C h_{\ell}^{k}. $$

A sufficient condition for (4.3) is therefore

$$ (4.4) \qquad h_{\ell+1} \le h_{\ell} e^{-\mu T_{\ell}/k}, \quad \ell \ge 1. $$

Condition (4.4) provides a practical guideline for refining the accuracy of the time-parallel solver from one window to the next so as to preserve the exponential convergence rate of the observer. Under this condition, the parallel observer converges with the same asymptotic rate $\mu$ as the serial observer.

In practice, the state $x(t)$ is not available and the stopping criterion is formulated in terms of the observable output. Accordingly, the algorithm is terminated when $\|y(t) - \hat{y}(t)\| \le \text{tol}$. When the convergence rate $\mu$ is not known a priori, it may be estimated numerically, for instance by power-method-based techniques [31, 25].
**4.4. Proposed strategy and its optimisation.** As we want to set a non in-trusive method, we present a technique where the Paraexp algorithm can be used as a black box. As a matter of a fact, our approach does not depend on the time paral-lelization method which is used in an inner loop on each window, and an algorithm like Parareal also fits this setting perfectly. Nevertheless, we remark that the values of the states (observed quantities and observer state) are only required at time $T_\ell$. Hence, some optimization techniques can be applied.

First, we do not need to compute the solution of Type 2 sub-problems over the entire window $W_\ell$, but only at time $T_{\ell+1}$. This means that only $N$ matrix exponentials need to be computed, instead of $N \frac{T}{\Delta t}$. Second, as all the windows $W_\ell$ have the same length, and the matrix exponential is the same for all windows, the $N$ matrix exponential computations need to be done only once, and can be used again from one window to another. Note that this is the case only when the computation method of the matrix exponential does not rely on the initial condition. This argument is therefore valid for scaling and squaring with Padé approximations or Chebyshev polynomials, but is not relevant in the case of rational Krylov. Those properties leads to a very low computational cost when solving Type 2 sub-problems. In practical cases, this part can be neglected so that our method is fully efficient. We therefore design an updated version of the Paraexp algorithm in Algorithm 1. In this updated version, lines (6) and (7) are virtually free in computation costs compare to the rest of the procedure.

<caption>**Algorithm 1:** updated Paraexp: ParaExp_end</caption>

Input : $N, M, g(t), \Delta t, x_0, t_0, T$
Output: $x_\parallel(T)$
1 for $j$ in $\{1, \dots, N\}$: $t_j \leftarrow t_0 + j \frac{T}{N}$;
2 for $j$ in $\{1, \dots, N\}$: precompute: $\exp\{M(T + t_0 - t_{j-1})\}$;
3 $x_{t_0}^S \leftarrow x_0$;
4 for $j = 1, \dots, N$ do in parallel
5 $x_{t_j}^S \leftarrow \text{solve}(\frac{d}{dt}x = Mx + g(t); x_0 = 0; [t_{j-1}, t_j]; \Delta t) [\text{end}]$;
6 for $j = 1, \dots, N$ do in parallel
7 $x_{t_{j-1}}^H = \exp\{M(T + t_0 - t_{j-1})\} x_{t_{j-1}}^S$ (precomputed);
8 $x_\parallel(T) \leftarrow x_{t_N}^S + \sum_{j=1}^N x_{t_{j-1}}^H$;

The algorithm corresponding to this strategy is presented in Algorithm 2, where $A, B, u(t), y(t), C$ and $\hat{x}_0$ come from (2.2), $T$ is the size of the windows, $h_1$ the desired initial precision, and $\mathcal{L} = \{\lambda_1, \lambda_2, \dots, \lambda_n\}$ the desired eigenvalues for the pole placement. The notation Para-in-time refers to a generic parallel in time procedure
over $N$ computers, with a precision given by $h_\ell$, solving $\frac{d}{dt}x(t) = Mx(t) + g(t)$ at time $t = T_{\ell-1} + T = T_\ell$, with $\hat{\mathbf{X}}_{\ell-1}^{h_\ell-1}$ as an initial condition at time $t = T_{\ell-1}$. The previously exposed ParaExp-end algorithm can be used in this instance instead.

<caption>**Algorithm 2:** Strategy, with a generic Parallel-in-time procedure</caption>

Input : $A, B, u(t), y(t), C, \hat{x}_0, T, N, h_1, (\mathcal{L})$
Output: $\{\hat{\mathbf{X}}_\ell^{h_\ell}\}_{\ell \ge 0}$
1 if PDE = true then
2 | $L \leftarrow \gamma \cdot C^T$;
3 else
4 | $\lfloor (V, L) \leftarrow \text{place}(A^T, C^T, \mathcal{L}) \rfloor$;
5 | $M \leftarrow (A - LC)$;
6 | $g(t) \leftarrow Bu(t) + Ly(t)$;
7 | $\hat{\mathbf{X}}_0^{h_0} \leftarrow \hat{x}_0$;
8 | $\ell \leftarrow 1$;
9 while $\|C\hat{\mathbf{X}}_\ell^{h_\ell} - \mathbf{y}(T_\ell)\| > \text{tol}$ do
10 | $\quad (T_{\ell-1}, T_\ell) \leftarrow ((\ell-1) \cdot T, \ell \cdot T)$;
11 | Determine $h_\ell$ such that (4.4) holds;
12 | $\quad \hat{\mathbf{X}}_\ell^{h_\ell} \leftarrow \text{Para-in-time}(N, M, g(t), h_\ell, \hat{\mathbf{X}}_{\ell-1}^{h_\ell-1}, T_{\ell-1}, T)$;
13 | $\quad \ell \leftarrow \ell + 1$;

5. Numerical experiments. Numerical experiments are obtained using Matlab [19]. The efficiency is computed as follow:

$$375 \qquad (5.1) \qquad \text{Efficiency} := \frac{\tau_s}{N \times \tau_p^{\mu,N,T}},$$

with, as stated before, $N$ the number of computers, $T$ the size of windows, and $\mu$ the rate of convergence. The CPU time required to achieve $\|\varepsilon\|(t) \le \text{tol}$ using a parallel solver is $\tau_p^{\mu,N,T}$ and the CPU time required to achieve $\|\varepsilon\|(t) \le \text{tol}$ using a sequential solver is $\tau_s$. The computational times for the parallel solver are simulated parallel computation times. They are defined by enforcing Matlab to use only one processor with input `maxNumCompThreads(1)`, and taking the reference time for performing a parallel `for` loop as the longest time to calculate one of the `for` instructions. The Speedup can also be computed by:

$$384 \qquad (5.2) \qquad \text{Speedup} := \frac{\tau_s}{\tau_p^{\mu,N,T}}.$$

Due to the design of our method, the time $t$ at which the parallel solver matches the stopping criteria can be greater than the time $t$ at which the serial solver matches the stopping criteria. We call those respectively $t = t_p^{\text{stop}}$ and $t = t_s^{\text{stop}}$. Indeed, once the window size $T$ is set, if the convergence criteria is not met at the end of the window, the procedures continues onto another window, even if it was close to convergence. For this reason, we introduce what we call the overshoot of our method:

$$391 \qquad (5.3) \qquad \text{Overshoot} := \frac{t_p^{\text{stop}}}{t_s^{\text{stop}}}.$$
5.1. A small system. We consider here an academic small numerical example.

We refer to [14] for the working physics. An autopilot is to be designed for a subma-
rine to maintain a constant depth under severe wave disturbances and yields of the
following system:

$$
\begin{align*}
\frac{d}{dt} \begin{pmatrix} w \\ q \\ z \\ \theta \end{pmatrix} &= A \begin{pmatrix} w \\ q \\ z \\ \theta \end{pmatrix} + B \begin{pmatrix} \delta_B \\ \delta_S \end{pmatrix} \\
(5.4) \qquad \begin{pmatrix} z \\ \theta \end{pmatrix} &= C \begin{pmatrix} w \\ q \\ z \\ \theta \end{pmatrix}
\end{align*}
$$

where:

$$
A = \begin{pmatrix} -0.038 & 0.896 & 0 & 0.0015 \\ 0.0017 & -0.092 & 0 & -0.0056 \\ 1 & 0 & 0 & -3.086 \\ 0 & 1 & 0 & 0 \end{pmatrix}, \quad B = \begin{pmatrix} -0.0075 & -0.023 \\ 0.0017 & -0.0022 \\ 0 & 0 \\ 0 & 0 \end{pmatrix},
$$

and $C = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}$ so that $x = (w, q, z, \theta)^T$ and $y = (z, \theta)^T$. This system is completely controllable and observable as per Theorem 2.2 and Theorem 2.4. We use the following control term: $u(t) = (\delta_B(t), \delta_S(t))^T = 100 \cdot (\sin(t), \exp(-t))^T$. We construct an observer $\hat{x}$ for the system (5.4), as in (2.2). The matrix $L$ is chosen with a pole placement procedure using the place function in Matlab. For more details about this algorithm that allows us to choose the eigenvalues of the matrix $(A - LC)$ we refer to [21]. This procedure gives us accordingly a matrix $L$ (as well as the matrix $V$ whose columns are the eigenvectors of the matrix $(A - LC)$ for the desired eigenvalues that we use thereafter). Here, solving Type 1 subproblem is done with the usual RK4 scheme, and solving the Type 2 subproblem is direct and exact. Indeed it is well known that the exponential of a diagonalised matrix is easy to compute:

$$
\exp(t(A - LC)) = \exp(tM) = \exp(tVDV^{-1}) = V \exp(tD)V^{-1},
$$

where $D$ is a diagonal matrix whose entries are the chosen eigenvalues of $(A - LC)$,
so that $\exp(tD)$ is a diagonal matrix whose entries are the exponential at time $t$ of
the chosen eigenvalues of $(A - LC)$.

We perform 500 runs of the problem, using a random initial value $\hat{x}_0 \in (0, 5)^4$.
We iterate the solve over different values of $T$ and $\mu$. Results are presented in 5.1.
We see that our efficiency is close to 100% as $\mu$ is small (meaning that the rate of
convergence is small), and gets lower as $\mu$ increases up to a point where the algorithm
becomes inefficient. This comes from the fact that the parallel observer needs at least
2 windows to converge, in other words it can not have converged before $t = 2T$. If the
sequential observer has a stopping time way smaller than $2T$, then it is more efficient.
This situation happens when $T$ and/or $\mu$ are large. This numerical example allows
us to apply the pole placement procedure described in 2.1. We are therefore entitled
to choose a very large $\mu$ for extremely fast convergence of the observer, making the
proposed strategy obsolete. But as mentionned in 2.2 we do not have complete control
over $\mu$ on some systems arising from discretized PDEs. Nonetheless this example
<figure><img src="image_3.png" /><figcaption>FIGURE 5.1. Results for the submarine example.</figcaption></figure>

426 validates our technique for some choices of T and µ as we were able to reach a 100% accuracy. 427

5.2. A 1-dimensional wave equation. We take here an academical partial differential equation example. We consider the following wave equation:

$$
(5.5) \qquad \partial_t^2 u(x,t) = \Delta u(x,t), \quad x \in [0, 2\pi],
$$

with initial condition $u(x,0) = \sin(x)$, and Dirichlet boundary condition $u(0,t) = 0$,
$u(2\pi, t) = 0$. The solution is known as: $u_{\text{true}}(x,t) = \sin(x)\cos(t)$, as well as its time
derivatives $(u_{\text{true}})_t(x,t) = -\sin(x)\sin(t)$ and $(u_{\text{true}})_{tt}(x,t) = -\sin(x)\cos(t)$. As our
aim here is to investigate the performance of our proposed method, we are interested
in an observer that must converge up to machine precision. We therefore choose a
slightly different setting, namely a state that is semi-discretized in space:

$$
(5.6) \qquad \left\{
\begin{array}{ll}
\displaystyle \partial_t^2 u(x,t) = Du(x,t) + g(x,t), & x \in [0, 2\pi], \\
\\
\displaystyle g(x,t) = (u_{\text{true}})_{tt}(x,t) - Du_{\text{true}}(x,t),
\end{array}
\right.
$$

438 in order to be able have convergence of the joint observer up to machine precision.
439 Here, the matrix $D$ represents the standard 3-points finite differnce Laplacian, The
440 source term $g$ is chosen in order to correct the inevitable numerical error introduced
441 by the spatial discretization. This allows the observer defined aftewards to converge
442 up to machine precision, and not just up to spatial discretization precision. Moreover,
443 we choose to solve (5.5) with a scheme proposed by A. Munch in [32]. Equation (5.6)
444 is then written in the form of a system of two first order PDE and semidiscretized in
445 space with $m_x = 20$ spatial points:

$$
\partial_t \begin{pmatrix} u(x,t) \\ u_2(x,t) \end{pmatrix} = (L_m)^{-1} \underbrace{\begin{pmatrix} 0 & I \\ D & 0 \end{pmatrix}}_{A} \begin{pmatrix} u(x,t) \\ u_2(x,t) \end{pmatrix} \\
+ (L_m)^{-1} \begin{pmatrix} 0 \\ \Theta(u_{\text{true}})_{tt}(x,t) - Du_{\text{true}}(x,t) \end{pmatrix},
$$
where $L_m$ and $\Theta$ are the matrices suggested by [32], i.e.:

$$ \Theta = \begin{pmatrix} 1-2\theta & \theta & 0 & \cdots & 0 \\ \theta & 1-2\theta & \theta & \cdots & \vdots \\ 0 & \theta & 1-2\theta & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \theta \\ 0 & \cdots & 0 & \theta & 1-2\theta \end{pmatrix}, \quad L_m = \begin{pmatrix} I & 0 \\ 0 & \Theta \end{pmatrix}, $$

with parameter $\theta = 1/4$. In this framework, (2.2) reads:

$$ \partial_t \begin{pmatrix} \hat{u}(x,t) \\ \hat{u}_2(x,t) \end{pmatrix} = (L_m)^{-1} \underbrace{\begin{pmatrix} 0 & I \\ D & -\gamma C^T C \end{pmatrix}}_{M} \begin{pmatrix} \hat{u}(x,t) \\ \hat{u}_2(x,t) \end{pmatrix} + (L_m)^{-1} \begin{pmatrix} 0 \\ \Theta(u_{\text{true}})_{tt}(x,t) - Du_{\text{true}}(x,t) + \gamma C^T C u_{\text{true}}(x,t) \end{pmatrix}, $$

where we fix $\gamma = 8$, and with $C$ the discretization of our observation operator, namely $\mathbf{1}_{[\pi/2,3\pi/2]}$.

In order to investigate the performance of the proposed method, we set ourselves to compare it to the solve of (5.7) with a classical time-integration scheme RK4, both up to a certain stopping criteria. The nested grid on which our method choose the appropriate time step for solving Type 1 according to (5.7) is $\Delta t = \frac{1}{2^4}$, $i \in [9, 14]$. The convergence rate $\mu$ is computed in advance for this matter with an eigenvalue analysis of the matrix $M$. Time step $\Delta t = \frac{1}{2^9}$ is the largest timestep that allows convergence of the RK4 scheme, and $\Delta t = \frac{1}{2^{14}}$ is the largest timestep that brings our observer error to computer precision. Convergence plots of the observer run with a RK4 scheme and the various timesteps of the nested-grid is shown on 5.2. We see that

<figure><img src="image_4.png" /><figcaption>FIGURE 5.2. Convergence of (5.7) on the proposed nested grid.</figcaption></figure>
we have quite the flexibility over 6 different time steps to be refined as the observer converges.

Regarding matrix exponential computation, the unusual eigenstructure of the matrix of the dynamics of the observer (complex eigenvalues sitting on the left half part of the complex plane) prevent us from using easily classical methods such as explained in 3.2. The test were therefore carried out with Matlab `expm` function. In the future, we want to execute this test with a rational Krylov method, using S. Güttel toolbox [15] in order to pre-compute optimal poles for the matrix.

The two variables of interest are $T \in \{1, 2, 5, 12, 28, 65, 150\}$, the size of the $W_\ell$ windows, and the number of processors available to compute the Paraexp part of our scheme, $N \in \{2, 4, 8, 16, 32, 64, 128\}$. Stopping criteria tol is set at $10^{-12}$. Results are shown in 5.3. Those results are very encouraging. Recall the formula for computing

<figure><img src="image_5.png" /><figcaption>FIGURE 5.3. Results for the 1D wave equation.</figcaption></figure>

the efficiency (5.1). For a simple Paraexp solve, over the same time interval, and with the same solver as a classical integration scheme, the maximal theoretical possible efficiency is 1. Here our method allows us to exceed this threshold up to an efficiency of almost 3. This is thanks to the refinement of the timestep to solve Type 1. Those results are found for a size of windows $T = 28$. Note that the most performing size of windows $T = 28$ is equivalent to approximately $0.18t_s^{\text{stop}}$, with here $t_s^{\text{stop}} \approx 150$. It is also interesting to look at the overshoot of the method. We clearly see large overshoot for large values $T$, pushing the scheme to compute a solution up to time $t = 450$ when the classical scheme finds a convergence around time $t = 165$. This reduces a lot the efficiency of the method but the refinement of the timestep approach allows us to still get interesting results. As a general rule, we observe that our efficiency is globally very good, at least equal to 1.

These numerical results, while very encouraging, should be viewed with certain
487 reservations. First, the size of the problem is still quite small, with only 40 unknowns
488 in the resulting system of ODEs. This is due to two main limitations: the CFL
489 condition and the cost of running this test. Indeed, if the space grid is refined, so
490 is the CFL condition, and therefore the number of nested grids available to us. For
491 instance, choosing a $m_x = 100$ points discretization, resulting in a 200 unknowns
492 system of ODEs, leads to a nested grid with only two levels: $\Delta t = \frac{1}{2^{13}}$ and $\Delta t = \frac{1}{2^{14}}$.
493 The argument could be made that we can still run our method with de facto the
494 finest grid for the RK4 solve in the parallel scheme. This would result in a maximal
495 efficiency of 1. Moreover, the time span over which we are sure that the observer
496 attains full convergence up to machine precision is very large with respect to the time
497 step. Taking the previous example of $m_x = 100$ points discretization in space with
498 the largest time step available $\Delta t = \frac{1}{2^{13}}$ that respects CFL condition means more than
499 $1.2 \cdot 10^6$ iteration of the RK4 scheme with a matrix of size $(2m_x, 2m_x)$ which is not
500 feasible on a desktop computer in a reasonable time (it is for a single run, but not for
501 50 runs in order to compare the methods). This limitation is also the explanation as
502 to why no 2D tests have been done at this time.

503 **5.3. A linearized water waves system.** We are also interested in addressing
504 a linearized water waves system in the setting of [35]. The setting of this numerical
505 example has a lot of similarities with the one of 5.2. We are then interested is the
506 following system of equations, with $x \in [0, \pi]$ and periodic boundary conditions:

$$
507 \quad (5.8) \qquad \left\{
\begin{aligned}
\partial_t \phi(x,t) &= -g\eta(x,t), \\
\partial_t \eta(x,t) &= \mathcal{G}(\phi(x,t)), \\
y(x,t) &= \eta.
\end{aligned}
\right.
$$

508 where $\phi$ represents the velocity potential, $\eta$ represents the free surface elevation, $y$ the
509 observation of the free surface and $g$ the gravity constant $g = 9.81$. As we set our sevles
510 in a deep water regime, $\mathcal{G}$ is defined as the Fourier multiplier $\mathcal{G}(f) : \xi \mapsto \xi \mathcal{F}(f)(\xi)$ for
511 a given frequency $\xi \in \mathbb{R}$, where $\mathcal{F}$ is the Fourier transform. A solution of (5.8) can
512 be defined by:

$$
513 \quad \left\{
\begin{aligned}
\phi_{\text{true}}(x,t) &= \sum_{n=0}^{N_{\text{fr}}} \frac{a_n g}{o_n} \cos(o_n t - k_n x), \\
\eta_{\text{true}}(x,t) &= \sum_{n=0}^{N_{\text{fr}}} a_n \sin(o_n t - k_n x),
\end{aligned}
\right.
$$

514 with its derivatives:

$$
515 \quad \left\{
\begin{aligned}
(\phi_{\text{true}})_t(x,t) &= \sum_{n=0}^{N_{\text{fr}}} -a_n g \sin(o_n t - k_n x), \\
(\eta_{\text{true}})_t(x,t) &= \sum_{n=0}^{N_{\text{fr}}} a_n o_n \cos(o_n t - k_n x),
\end{aligned}
\right.
$$

516 where $a_n = \frac{N_{fr}-n}{4N_{fr}^2}$ represent wave amplitudes, $o_n = \sqrt{k_n/g}$ the angular frequency,
517 and $k_n = 2n$ the angular wave vector. Those are classical real solution of the water
518 waves system in this setting [8]. Initial conditions are $\phi(x,0) = \phi_{\text{true}}(x, T_{\text{start}})$,
519 $\eta(x,0) = \eta_{\text{true}}(x, T_{\text{start}})$ with $T_{\text{start}} = 100$. We show a representation of the free
520 surface $\eta$ between time $t = T_{\text{start}}$ and $t = T_{\text{start}} + 5$ in 5.4.

521 We therefore construct the following observer:

$$
522 \quad \left\{
\begin{aligned}
\partial_t \hat{\phi}(x,t) &= -g\hat{\eta}(x,t), \\
\partial_t \hat{\eta}(x,t) &= \mathcal{G}(\hat{\phi}(x,t)) + \gamma(\eta_{\text{true}}(x,t) - \hat{\eta}(x,t)).
\end{aligned}
\right.
$$
<figure><img src="image_6.png" /><figcaption>FIGURE 5.4. Representation of the free surface $\eta$ between time $t = T_{start}$ and $t = T_{start} + 5$.</figcaption></figure>

As shown in [35], this system is observable and controllable under the considered setting.

As for the previous example in 5.2, we semi-discretize in space this equation with $m_x = 64$ points, with the addition of a source term in order to be able to converge up to machine precision, and not be limited by the space discretization scheme. We would therefore have the semi-discretized in space state dynamic:

$$ (5.9) \qquad \frac{d}{dt} \begin{pmatrix} \phi \\ \eta \end{pmatrix} = \begin{pmatrix} O & -gI \\ G & O \end{pmatrix} \begin{pmatrix} \phi \\ \eta \end{pmatrix} + \begin{pmatrix} (\phi_{\text{true}})_t + g\eta_{\text{true}} \\ (\eta_{\text{true}})_t - G\phi_{\text{true}} \end{pmatrix}, $$

with $G$ the matrix corresponding to discretized operator $G$, $I$ the identity matrix and $O$ the null matrix. The observer that goes along with (5.9) is then:

$$ (5.10) \qquad \frac{d}{dt} \begin{pmatrix} \hat{\phi} \\ \hat{\eta} \end{pmatrix} = \underbrace{\begin{pmatrix} O & -gI \\ G & -\gamma \end{pmatrix}}_{M} \begin{pmatrix} \hat{\phi} \\ \hat{\eta} \end{pmatrix} + \begin{pmatrix} (\phi_{\text{true}})_t + g\eta_{\text{true}} \\ (\eta_{\text{true}})_t - G\phi_{\text{true}} + \eta_{\text{true}} \end{pmatrix}. $$

This setting being very similar as the one of 5.2, the nested grid on which our method choose the appropriate time step for solving Type 1 according to (5.10) is $\Delta t = \frac{1}{2\pi}$, $i \in [4, 13]$. The convergence rate $\mu$ is computed in advance for this matter with an eigenvalue analysis of the matrix $M$. Convergence plots of the observer run with a RK4 scheme and the various timesteps of the nested-grid is shown on 5.5.

Our two variable of interest are the same as the above numerical example, with $T = \{1, 2, 3, 5, 9, 16, 28\}$ and $N = \{2, 4, 8, 16, 32, 64, 128\}$. Stopping criteria tol is set at $10^{-12}$. Results are shown in 5.6. They are very similar to the one of 5.2. We also observe globally very encouraging results with a minimum efficiency around 1, and some great results of efficiency around 3 for a size of window $T = 5$. For large values of $T$, the overshoot becomes large. Note that the most performing size of windows $T = 5$ is equivalent to approximately $0.17t_s^{\text{stop}}$, with here $t_s^{\text{stop}} \approx 28$. This is very similar to what was found in 5.2. Overshoot also behaves as it should, increasing as $T$ becomes very large. The same reservations as for the previous example are also made for this one.
<figure><img src="image_7.png" /><figcaption>FIGURE 5.5. Convergence of (5.10) on the proposed nested grid.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 5.6. Results for the linearized water waves system. Top: (left) Efficiency, (right) overshoot. Bottom: Speedup.</figcaption></figure>

**6. Perspectives.** The results presented in this article showcase the potential of combining the Paraexp algorithm with the Luenberger observer for efficient data assimilation. Moving forward, several directions for future research and numerical experiments can be envisioned, each aimed at extending the applicability and performance of the methods discussed here.
First, one of our primary goals is to extend the numerical experiments to larger-scale problems. While the current study focuses on small or moderate sized systems, applying the Paraexp strategy to higher-dimensional models will allow us to assess its robustness and scalability. These larger-scale tests would further highlight the computational efficiency gains achieved by the time-parallelization scheme and demonstrate its potential for practical use in fields such as fluid dynamics.

Additionally, we see an opportunity to apply this strategy to non-linear PDEs, leveraging the non-linear extension of the Paraexp algorithm [12]. This extension, while conceptually similar to Parareal, is more suited to systems where the non-linearity is not overly dominant. In these cases, the non-linear Paraexp algorithm can still retain some of the computational advantages of its linear counterpart while accommodating the complexities introduced by non-linearities. A targeted study on systems with mild non-linearities would be a natural next step to evaluate this approach and understand the trade-off in performance.

Beyond numerical testing, there are also potential adaptations of the current strategy. One such adaptation involves applying our strategy the back and forth nudging algorithm [1]. Our concept would need to be adapted, given that this algorithm operates over a finite time horizon. This would allow us to extend our approach to a modern data assimilation method that has proved its effectiveness.

Another possible extension involves dynamically adjusting the size of the time windows, $T$ in our strategy. In its current form, the window size remains fixed, but by allowing it to vary over time based on error estimates, we could optimize the computational load and improve its overall performance.

Finally, an intriguing avenue of research would be to adapt the current approach to non-deterministic settings, such as those encountered with Kalman filters [22, 24]. The Luenberger observer, as discussed in this paper, operates in a deterministic framework. However, extending the methodology to a stochastic setting, where measurements are uncertain or noisy, could open new doors for its application. The challenge would be to integrate the Paraexp strategy within the recursive update framework of Kalman filtering.

In conclusion, the work presented in this paper provides a foundation for future research in both theoretical and applied domains. By expanding the numerical tests to larger scales, exploring applications to non-linear systems, adapting the strategy to a broader class of nudging algorithms, and incorporating stochastic elements, we aim to continue exploring the possibilities of coupling time-parallelization and data assimilation methods.

## REFERENCES

[1] D. Auroux and J. Blum. Back and forth nudging algorithm for data assimilation problems. *Comptes Rendus Mathematique*, 340(12):873-878, 2005.

[2] C. Bardos, G. Lebeau, and J. Rauch. Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary. *SIAM Journal on Control and Optimization*, 30(5):1024-1065, 1992.

[3] M. Berljafa and S. Güttel. Generalized rational krylov decompositions with an application to rational approximation. *SIAM Journal on Matrix Analysis and Applications*, 36(2):894-916, 2015.

[4] M. Berljafa and S. Güttel. The rkfit algorithm for nonlinear rational approximation. *SIAM Journal on Scientific Computing*, 39(5):A2049-A2071, 2017.

[5] P. Bernard and V. Andrieu. Luenberger observers for non autonomous nonlinear systems. *IEEE Transactions on Automatic Control*, 64(1):270-281, Sept. 2018.

[6] D. Chapelle, N. Cindea, M. Buhan, and P. Moireau. Exponential convergence of an observer
603 based on partial field measurements for the wave equation. *Mathematical Problems in Engineering*, 2012, 10 2012.

605 [7] J. Coron. *Control and Nonlinearity*. Mathematical surveys and monographs. American Mathematical Society, 2007.

607 [8] V. Duchêne. *Many Models for Water Waves*. Habilitation à diriger des recherches, Université de Rennes 1, July 2021.

609 [9] G. Freudenthaler, F. Göttsch, and T. Meurer. Backstepping-based extended luenberger observer design for a burgers-type pde for multi-agent deployment. *IFAC-PapersOnLine*, 50(1):6780–6785, 2017. 20th IFAC World Congress.

612 [10] M. Gander. Talk by Martin Gander (University of Geneva) at ENLA seminar. https://www.youtube.com/watch?v=DFY_w6KnhfU, 6 2020.

614 [11] M. Gander and S. Güttel. Paraexp: a parallel integrator for linear initial-value problems. *SIAM Journal on Scientific Computing*, 35, 03 2013.

616 [12] M. Gander, S. Güttel, and M. Petcu. *A Nonlinear ParaExp Algorithm*, pages 261–270. 01 2018.

618 [13] M. Gander and S. Vandewalle. Analysis of the parareal time-parallel time-integration method. *SIAM Journal on Scientific Computing*, 29(2):556–578, 2007.

620 [14] M. Grimble, G. Van der Molen, and E. Liceaga-Castro. Submarine depth and pitch control. In *Proceedings of IEEE International Conference on Control and Applications*, pages 953–958 vol.2, 1993.

623 [15] S. Güttel. Rktoolbox guide, 7 2020.

624 [16] G. Haine. Recovering the observable part of the initial data of an infinite-dimensional linear system with skew-adjoint generator. *Mathematics of Control Signals and Systems*, 26:435–462, 01 2014.

627 [17] G. Haine and K. Ramdani. Reconstructing initial data using observers: Error analysis of the semi-discrete and fully discrete approximations. *Numerische Mathematik*, 120:307–343, 07 2011.

630 [18] N. Higham. The scaling and squaring method for the matrix exponential revisited. *SIAM Rev.*, 51(4):747–764, 2009.

632 [19] T. M. Inc. Matlab version: 24.1.0.2628055 (r2024a) update 4, 2024.

633 [20] J. A. Infante and E. Zuazua. Boundary observability for the space semi-discretizations of the 1-d wave equation. *ESAIM: Modélisation mathématique et analyse numérique*, 33(2):407–438, 1999.

634 [21] N. K. N. J. Kautsky and P. V. Dooren. Robust pole assignment in linear state feedback. *International Journal of Control*, 41(5):1129–1155, 1985.

638 [22] R. Kalman. A New Approach to Linear Filtering and Prediction Problems. *Journal of Basic Engineering*, 82(1):35–45, 03 1960.

639 [23] R. Kalman. On the general theory of control systems. *IFAC Proceedings Volumes*, 1(1):491–502, 1960. 1st International IFAC Congress on Automatic and Remote Control, Moscow, USSR, 1960.

643 [24] R. Kalman and R. Bucy. New Results in Linear Filtering and Prediction Theory. *Journal of Basic Engineering*, 83(1):95–108, 03 1961.

645 [25] S. Kamvar, T. Haveliwala, C. Manning, and G. Golub. Extrapolation methods for accelerating pagerank computations. In *Proceedings of the 12th International Conference on World Wide Web*, WWW ’03, page 261–270, New York, NY, USA, 2003. Association for Computing Machinery.

649 [26] P. Ligarius and J.-F. Couchouron. Asymptotic observers for a class of evolution equations. a nonlinear approach. *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics*, 324(3):355–360, 1997.

652 [27] K. Liu. Locally distributed control and damping for the conservative systems. *Siam Journal on Control and Optimization - SIAM J CONTR OPTIMIZAT*, 35, 09 1997.

654 [28] D. Luenberger. Observers for multivariable systems. *Automatic Control, IEEE Transactions on*, AC-11:190 – 197, 05 1966.

656 [29] D. Luenberger. *Introduction to Dynamic Systems: Theory, Models, and Applications*. John Wiley & Sons, New York, 1979.

658 [30] Y. Maday and G. Turinici. A parareal in time procedure for the control of partial differential equations. *Comptes Rendus Mathematique*, 335(4):387–392, 2002.

660 [31] R. Mises and H. Pollaczek-Geiringer. Praktische verfahren der gleichungsauflösung. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik*, 9(2):152–164, 1929.

663 [32] A. Münch. Famille de schémas implicites uniformément contrôlables pour l'équation des ondes 1-d. *Comptes Rendus Mathematique*, 339(10):733–738, 2004.
665

666

667

668

669

670

671

672

673

674

675

[33] K. Ramdani, T. Takahashi, and M. Tucsnak. Uniformly exponentially stable approximations for a class of second order evolution equations. Application to LQR optimization problems. *ESAIM: Control, Optimisation and Calculus of Variations*, 13(3):503-527, 2007.

[34] S. Riffo, F. Kwok, and J. Salomon. Time-parallelization of sequential data assimilation problems. preprint, Dec. 2022.

[35] Y. Yu, H.-L. Pei, and C.-Z. Xu. Estimation of velocity potential of water waves using a luenberger-like observer. *Science China Information Sciences*, 63, 12 2020.

[36] M. Zeitz. The extended luenberger observer for nonlinear systems. *Systems & Control Letters*, 9(2):149-156, 1987.

[37] E. Zuazua. Propagation, observation, and control of waves approximated by finite difference methods. *SIAM Review*, 47(2):197-243, 2005.