# Stein-based preconditioners for weak-constraint 4D-var

Davide Palitta<EMAIL_ADDRESS>Jemima M. Tabeart
<EMAIL_ADDRESS>
###### Abstract

Algorithms for data assimilation try to predict the most likely state of a
dynamical system by combining information from observations and prior models.
Variational approaches, such as the weak-constraint four-dimensional
variational data assimilation formulation considered in this paper, can
ultimately be interpreted as a minimization problem. One of the main
challenges of such a formulation is the solution of large  linear systems of
equations which arise within the inner linear step of the adopted nonlinear
solver. Depending on the selected approach, these linear algebraic problems
amount to either a saddle point linear system or a symmetric positive definite
(SPD) one. Both formulations can be solved by means of a Krylov method, like
GMRES or CG, that needs to be preconditioned to ensure fast convergence in
terms of the number of iterations. In this paper we illustrate novel,
efficient preconditioning operators which involve the solution of certain
Stein matrix equations. In addition to achieving better computational
performance, the latter machinery allows us to derive tighter bounds for the
eigenvalue distribution of the preconditioned linear system for certain
problem settings. A panel of diverse numerical results displays the
effectiveness of the proposed methodology compared to current state-of-the-art
approaches.

###### keywords:

4D-Var , data assimilation , preconditioning , Stein equations

[UniBo]organization=Dipartimento di Matematica and AM2, Alma Mater Studiorum -
Università di Bologna, addressline=Piazza di Porta S. Donato, 5, I-40127,
city=Bologna, country=Italy,

[edi]organization=School of Mathematics, University of Edinburgh,
addressline=Peter Guthrie Tait Road, city=Edinburgh, country=UK [ox]
organization=Now at: Mathematical Institute, University of Oxford,
addressline=Andrew Wiles Building, city=Oxford, country=UK

##  1 Introduction

Given a computational model for a dynamical system, data assimilation aims to
merge observational, measured data of that system with prior model information
to obtain a better estimate of the system state at a specified time. The most
mature application of data assimilation is to numerical weather prediction
(NWP), where it is used to obtain the initial conditions for forecasts [7],
but in recent years data assimilation approaches have been studied more
broadly within earth sciences, ecology, and neuroscience; see, e.g., [35, 30,
41]. In particular, observations $y_{i}\in\mathbb{R}^{p_{i}}$ at time
$t_{i}\in[t_{0},t_{N}]$ are combined with prior information
$x_{b}\in\mathbb{R}^{s}$ from a model to compute the most likely state
$x_{i}\in\mathbb{R}^{s}$ of the system at time $t_{i}$. It is typically
assumed that the background state $x_{b}$ can be written as
$x_{b}=x_{0}^{t}+\epsilon^{b}$ where $x_{0}^{t}$ denotes the true initial
state of the system with the error $\epsilon^{b}$ . This error is distributed
according to a normal distribution with error covariance matrix
$B\in\mathbb{R}^{s\times s}$ and zero mean, i.e.,
$\epsilon^{b}\sim\mathcal{N}(0,B)$. Similarly, we write each observation in
terms of the true initial state as
$y_{i}=\mathcal{H}_{i}(x_{i}^{t})+\epsilon_{i}$ with the observation error
$\epsilon_{i}\sim\mathcal{N}(0,R_{i})$ for all $i=0,\ldots,N$. In order to map
between observation and state space, we introduce a possibly nonlinear
operator $\mathcal{H}_{i}:\mathbb{R}^{s}\rightarrow\mathbb{R}^{p_{i}}$ which
connects the true state $x_{i}^{t}$ at time $t_{i}$ and the observational data
$y_{i}$.

One of the notable aspects of the weak-constraint four-dimensional variational
assimilation problem (4D-Var) is the propagation of the computed states. If
$x_{i-1}$ denotes the state variable at time $t_{i-1}$, this is propagated to
the next observation time $t_{i}$ via an imperfect forecast model
$\mathcal{M}_{i}$ such that $x_{i}=\mathcal{M}_{i}(x_{i-1})+\epsilon_{i}^{m}$
where $\epsilon_{i}^{m}\sim\mathcal{N}(0,Q_{i})$ for all $i=1,\ldots,N$. This
is probably the main difference between the weak- and strong-constrained
4D-Var approaches. Indeed, in the latter methodology the forecast model is
supposed to be exact.

The ultimate goal of weak-constrained 4D-Var is then minimizing the following
functional

$\displaystyle J(x)=$
$\displaystyle(x_{0}-x_{b})^{T}B^{-1}(x_{0}-x_{b})+\sum_{i=0}^{N}(y_{i}-\mathcal{H}_{i}(x_{i}))^{T}R_{i}^{-1}(y_{i}-\mathcal{H}_{i}(x_{i}))$
$\displaystyle+\sum_{i=1}^{N}(x_{i}-\mathcal{M}_{i}(x_{i-1}))^{T}Q_{i}^{-1}(x_{i}-\mathcal{M}_{i}(x_{i-1})),$
(1)

where $x=(x_{0}^{T},\ldots,x_{N}^{T})^{T}\in\mathbb{R}^{(N+1)s}$ collects all
the state variables $x_{0},\ldots,x_{N}$.

The objective function (1) is often minimized by means of an incremental
approach [8]. Roughly speaking, this consists of a Gauss-Newton scheme where
at each iteration a linearised problem needs to be solved. It has been shown
that such a linear, inner problem can be reformulated as a large, sparse,
symmetric, but also very structured saddle point linear system; see, e.g.,
[12, 11, 18]. A more traditional approach consists of solving the SPD linear
system stemming from the quadratic optimization problems arising from the
adopted Gauss-Newton procedure; see e.g. [10, 49].

Krylov methods like the Generalized Minimal RESidual method (GMRES) [40] and
the MINimal RESidual method (MINRES) [32] are powerful tools for the solution
of saddle point linear systems. See, e.g., the survey paper [2].  Similarly,
the Conjugate Gradient method (CG) [21] is the most commonly used solver for
SPD linear systems. In both scenarios, it is vital to choose good
preconditioners for the adopted iterative scheme to ensure fast convergence in
terms of both the number of iterations and wallclock times.

A variety of preconditioners for the saddle point and SPD 4D-Var problems have
been proposed in the literature; see, e.g., [14, 17, 47, 9, 50]. While these
operators enjoy some appealing features, e.g., they guarantee parallisability
in the saddle-point context, they also neglect important features of the
original linear system to achieve affordable computational costs. This worsens
the capability of the preconditioners to reduce the overall iteration count.
In the operational NWP setting this is particularly problematic, as in
practice the maximum number of iterations is capped by a very small number
compared to the dimension of the problem [7]. Therefore, any preconditioning
method that reduces the iteration count without dramatically increasing the
computational cost is likely to be highly beneficial.

In this work we propose to fully exploit the inherent block structure of  both
the saddle point and SPD formulations within a matrix-oriented GMRES/CG
approach, see, e.g., [14, 44, 27]. Such machinery naturally leads to the
design of more efficient preconditioning operators with Kronecker structure.
These new preconditioners yield beneficial theoretical properties, thus
achieving faster convergence in terms of number of iterations than state-of-
the-art approaches, while maintaining a low computational cost and an easy-to-
parallelize nature.

The framework proposed in this paper requires moderate values of $p$ and $s$
(e.g. $\mathcal{O}(10^{3})$) to be computationally successful. We must mention
that these restrictions on the problem dimensions may be unrealistic for NWP
applications but can be reasonable for other data assimilation problems such
as parameter estimation tasks for low-dimensional parameter domains see, e.g.,
[36, 23] for some examples in agriculture. The fresh methodology we present in
this paper can also be successfully applied when data assimilation is combined
with model order reduction. This interesting scenario sees a first reduction
step aimed at reducing the state dimension. Then the reduced model is utilized
within the selected data assimilation approach; see, e.g., [26]. Weak-
constraint 4D-Var may be a particularly appropriate choice of data
assimilation scheme to be combined with model order reduction, as it is able
to take the additional model error coming from the reduction step into account
explicitly.  Moreover, the techniques we develop here serve as the initial
step towards novel procedures tackling large-scale problems such as the ones
stemming from NWP. This will be the subject of future work.

We additionally note that the weak-constraint 4D-Var problem requires more
computational resource than the strong-constraint formulation, in addition to
the need to prescribe model error covariance matrices. Alternative approaches
have been proposed which incorporate some model error information within a
strong-constraint 4D-Var approach via an inflated covariance approach, e.g.
[22, 15]. However, as these methods both require the inversion of the inflated
observation error covariance term, users are limited to the use of observation
error covariance matrices that are easy to invert. Our new approach allows
full flexibility for all data assimilation parameters, as well as revealing
and exploiting the Kronecker structure that is obscured in the usual primal
form. Indeed, the inflated covariance approach is expected to perform poorly
for large numbers of observation times, whereas the approach proposed in this
work scales very well with $N$.

Here is a synopsis of the paper. In section 2 we recall the formulation of the
SPD and saddle point linear systems stemming from weak-constraint 4D-Var. We
briefly introduce matrix-oriented GMRES and CG in section 2.1 and in section
2.2 we describe a general preconditioning framework to be embedded in these
routines. In section 3 we address the case of observation-time dependent
$\mathcal{M}_{i}$ and propose an original, efficient preconditioning operator.
The latter is very similar to the original operator with a single exception.
In particular, the original forward operator is approximated by a suitable,
observation-time independent one, namely
$\mathcal{M}_{i}\equiv\widehat{\mathcal{M}}$ for all $i$, in the
preconditioning operator. Thanks to this feature, we are able to show that the
inversion of a certain matrix $\mathbf{L}$, which is the predominant
computational bottleneck of state-of-the-art preconditioning procedures for
4D-Var, is in fact equivalent to solving a Stein matrix equation. In addition
to leading to some insights regarding the selection of a suitable
$\widehat{\mathcal{M}}$, we describe in section 3.2 how the matrix-oriented
perspective allows the efficient incorporation of information from the
observation term within the Schur complement of the saddle point system or,
equivalently, in the preconditioner of the SPD problem, by adapting an
approach proposed in [48]. We note that the observation term has often been
completely neglected in state-of-the-art preconditioners, but can be
incorporated approximately within the Kronecker preconditioning framework at a
moderate computational cost. In section 4, some further considerations are
given in the case that also the original forecast model $\mathcal{M}_{i}$ is
observation-time independent itself, namely $\mathcal{M}_{i}\equiv\mathcal{M}$
for all $i=1,\ldots,N$ in the forward operator. A number of numerical results
showing the potential of our fresh, successful strategy are reported in
section 5. We finish in section 6 by drawing some conclusions and presenting
possible outlooks.

Throughout the paper we adopt the following notation. Capital italic letters
($\mathcal{A}$) denote block matrices whose blocks have a Kronecker structure.
These blocks, and in general matrices having a Kronecker structure, are
denoted by capital bold letters ($\mathbf{A}$) whereas simple capital letters
($A$) are used for general matrices without any Kronecker structure. $I_{N}$
denotes the identity matrix of dimension $N$. The subscript is omitted
whenever the dimension of $I$ is clear from the context. The $i$-th vector of
the canonical basis of $\mathbb{R}^{N}$ is denoted by $e_{i}$. The Kronecker
product is denoted by $\otimes$, whereas $\circ$ represents the Hadamard
product. Given a matrix $X\in\mathbb{R}^{n\times n}$,
$\text{vec}(X)\in\mathbb{R}^{n^{2}}$ is the vector collecting the columns of
$X$ on top of one another. For instance, the variable $x$ in (1) can be
written as $x=\text{vec}([x_{0},\ldots,x_{N}])$. To conclude, $\lambda(A)$
denotes the spectrum of the matrix $A$, with
$\lambda_{\max}(A)=\lambda_{1}(A)\geq\lambda_{2}(A)\geq\ldots\geq\lambda_{N}(A)=\lambda_{\min}(A)$.

##  2 Linear system formulations

As previously mentioned, the vector state $x$ which minimizes (1) can be
computed by an incremental approach [8] where the cost function (1) is
approximated by a quadratic function of the increment $\delta
x^{(\ell)}=x^{(\ell+1)}-x^{(\ell)}$, with $x^{(\ell)}$ being the $\ell$-th
Gauss-Newton iterate. If $\delta x=\text{vec}([\delta x_{0},\ldots,\delta
x_{N}])$, the quadratic objective function is given by

$\displaystyle\delta J^{(\ell)}(\delta x)=$ $\displaystyle(\delta
x_{0}-b_{0}^{(\ell)})^{T}B^{-1}(\delta
x_{0}-b_{0}^{(\ell)})+\sum_{i=0}^{N}(d_{i}^{(\ell)}-H_{i}^{(\ell)}\delta
x_{i})^{T}R_{i}^{-1}(d_{i}^{(\ell)}-H_{i}^{(\ell)}\delta x_{i})$
$\displaystyle+\sum_{i=1}^{N}(\delta x_{i}-M_{i}^{(\ell)}\delta
x_{i-1})-c_{i}^{(\ell)})^{T}Q_{i}^{-1}(\delta x_{i}-M_{i}^{(\ell)}\delta
x_{i-1}-c_{i}^{(\ell)})),$

where $b_{0}^{(\ell)}=x_{b}-x_{0}^{(\ell)}$,
$d_{i}^{(\ell)}=y_{i}-\mathcal{H}_{i}(x_{i}^{(\ell)})$,
$c_{i}^{(\ell)}=\mathcal{M}_{i}(x_{i-1}^{(\ell)})-x_{i}^{(\ell)}$, and
$H_{i}^{(\ell)}$, $M_{i}^{(\ell)}$ are linearizations of $\mathcal{H}_{i}$ and
$\mathcal{M}_{i}$ about $x_{i}^{(\ell)}$, respectively. We note that
${B},{Q}_{i}$ and ${R}_{i}$ are covariance matrices so that they are symmetric
and positive semi-definite by construction. In addition, as inverse covariance
matrices are required in the objective function formulation (1) we assume
these matrices to be strictly positive definite. Therefore, by dropping the
Gauss-Newton index $(\ell)$ for better readability and assuming
$p_{0}=\ldots=p_{N}=p$,  CG can be employed to minimize $\delta J$ by solving
the following linear system

$\displaystyle\underbrace{(\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}+\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H})}_{\mathbf{S}}\delta
x=\mathbf{D}^{-1}b+\mathbf{L}^{T}\mathbf{H}^{T}\mathbf{R}^{-1}d,$ (2)

where $b=\text{vec}([b_{0},c_{1},\ldots,c_{N}])\in\mathbb{R}^{(N+1)s}$,
$d=\text{vec}([d_{0},\ldots,d_{N}])\in\mathbb{R}^{(N+1)p}$, and
$\mathbf{D},\mathbf{L}\in\mathbb{R}^{(N+1)s\times(N+1)s}$,
$\mathbf{R}\in\mathbb{R}^{(N+1)p\times(N+1)p}$,
$\mathbf{H}\in\mathbb{R}^{(N+1)p\times(N+1)s}$ are such that

$\displaystyle\mathbf{D}=\begin{pmatrix}B&&&\\\ &Q_{1}&&\\\ &&\ddots&\\\
&&&Q_{N}\\\ \end{pmatrix},$
$\displaystyle\quad\mathbf{L}=\begin{pmatrix}I&&&\\\ -M_{1}&I&&\\\
&\ddots&\ddots&\\\ &&-M_{N}&I\\\ \end{pmatrix},$
$\displaystyle\mathbf{R}=\begin{pmatrix}R_{0}&&&\\\ &R_{1}&&\\\ &&\ddots&\\\
&&&R_{N}\\\ \end{pmatrix},$
$\displaystyle\quad\mathbf{H}=\begin{pmatrix}H_{0}&&&\\\ &H_{1}&&\\\
&&\ddots&\\\ &&&H_{N}\\\ \end{pmatrix}.$

As an alternative to the quadratic minimization (2), $\delta x$ can be
computed by solving the following saddle point linear system [12]

$\underbrace{\begin{pmatrix}\mathbf{D}&0&\mathbf{L}\\\
0&\mathbf{R}&\mathbf{H}\\\ \mathbf{L}^{T}&\mathbf{H}^{T}&0\\\
\end{pmatrix}}_{=:\mathcal{A}}\begin{pmatrix}\delta\eta\\\ \delta\lambda\\\
\delta x\\\ \end{pmatrix}=\begin{pmatrix}b\\\ d\\\ 0\\\ \end{pmatrix},$ (3)

We note that both (2) and (3) contain a lot of inherent block structure111
Notice that the coefficient matrix in (2) is the Schur complement of
$\mathcal{A}$. This motivates the use of $\mathbf{S}$ for the former.. We
propose to fully exploit this structure by using matrix implementations of
iterative methods and designing preconditioners with explicit Kronecker
structure. We illustrate the main concept by considering a data assimilation
problem where the blocks of $\mathcal{A},\mathbf{S}$ and corresponding
preconditioners $\mathcal{P}$ have Kronecker structure. This could arise
naturally via consistent observation networks, with fixed observation and
model error statistics, at each observation time. In the case that
$Q_{1}=\cdots=Q_{N}\equiv Q$, $R_{0}=\cdots=R_{N}\equiv R$, and
$H_{1}=\cdots=H_{N}\equiv H$, we can write the terms above compactly by using
the inherent Kronecker structure

$\mathbf{D}=e_{1}e_{1}^{T}\otimes B+(I_{N+1}-e_{1}e_{1}^{T})\otimes
Q,\quad\mathbf{R}=I_{N+1}\otimes R,\quad\mathbf{H}=I_{N+1}\otimes H.$

In the more general setting where the covariance matrices and linearised
observation operator differ at each time, preconditioners with Kronecker
structure can be used within the same setting. We expect the strategy
presented in section 3 to be effective also in the case where we relax the
Kronecker assumptions on $\mathbf{R},$ $\mathbf{H}$ and $\mathbf{D}$. This
will be the subject of future work.

###  2.1 Matrix-oriented GMRES and CG

The Kronecker form of $\mathbf{S}$ and the blocks of the coefficient matrix
$\mathcal{A}$ naturally suggests the use of matrix-oriented Krylov subspace
methods to solve the linear systems (2) and (3). Depending on the adopted
preconditioning operator (see section 2.2), the most popular solution schemes
for solving (3) is GMRES, or MINRES if symmetry is preserved. Similarly, CG is
employed for (2). It is well-known that the original vector form of such
methods can be easily transformed in order to obtain a matrix-oriented
formulation of these routines. These implementations can be obtained by
exploiting the properties of the Kronecker product [51] and the equivalence
between the 2-norm of vectors and the matrix inner product, namely
$\text{vec}(A)^{T}\text{vec}(B)=\text{trace}((AB)^{T}AB)$. See, e.g., [14, 44,
34, 27] and Appendix B.

We would like to point out that none of the Krylov routines used to obtain the
results in section 5 are equipped with low-rank truncations as suggested in
[14, 44, 27]. These truncations steps are essential to obtain a feasibile
storage demand when very large dimensional problems are considered. Here we
suppose $p$, $s$, and $N$ to be moderate, say $\mathcal{O}(10^{3})$, so that
issues related to the memory consumption in our solvers do not occur in
general. Avoiding the employment of any low-rank truncation will be also
crucial to obtain the results stated in Proposition 2 and section 4.1.

###  2.2 Preconditioning operators

It is well-known that Krylov subspace techniques require effective
preconditioning operators to obtain fast convergence in terms of the number of
iterations.

In the 4D-Var context, many authors considered preconditioners for (2) of the
form

$\mathbf{\widehat{S}}:=\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L},$ (4)

which neglect the second term $\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$ in the
definition of $\mathbf{S}$. See, e.g., [14, 17, 47]. This leads to an easier-
to-invert preconditioning operator222 Notice that, due to its structure,
$\mathbf{L}$ is always nonsingular, regardless of the $M_{i}$’s. as
$\mathbf{\widehat{S}}^{-1}=\mathbf{L}^{-1}\mathbf{D}\mathbf{L}^{-T}$. A key
limitation of this preconditioner is that computation of the inverse operators
$\mathbf{L}^{-1}$ and $\mathbf{L}^{-T}$ requires many serial matrix products
and is thus not parallelisable. One of the main strategies to overcome this
issue is the introduction of a further layer of approximation related to
employing an operator $\mathbf{\widehat{L}}\approx\mathbf{L}$ in the
definition of $\mathbf{\widehat{S}}$, such that multiplication of a vector by
$\mathbf{\widehat{S}}^{-1}=\mathbf{\widehat{L}}^{-1}\mathbf{D}\mathbf{\widehat{L}}^{-T}$
can be distributed over multiple processors. Different options for the
selection of $\mathbf{\widehat{L}}$ can be found in, e.g., [12, 17, 14, 47].

For saddle-point linear systems of the form (3), some of the most commonly-
used preconditioners are the _block diagonal_ and _block triangular_
preconditioners. See, e.g., [3, 29, 17, 14]. In particular, the block diagonal
preconditioner is defined as follows

$\mathcal{P}_{\mathcal{D}}:=\begin{pmatrix}\mathbf{D}&&\\\ &\mathbf{R}&\\\
&&\mathbf{\widehat{S}}\\\ \end{pmatrix},$ (5)

where $\mathbf{\widehat{S}}$ is again such that
$\mathbf{\widehat{S}}\approx\mathbf{S}=\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}+\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$
is an approximation to the Schur complement $\mathbf{S}$ of $\mathcal{A}$,
and it is often of the form (4).

Similarly, the block triangular preconditioner is defined as follows

$\mathcal{P}_{\mathcal{T}}:=\begin{pmatrix}\mathbf{D}&0&{\mathbf{L}}\\\
&\mathbf{R}&\mathbf{H}\\\ &&\mathbf{\widehat{S}}\\\ \end{pmatrix}.$ (6)

A different class of preconditioning operators for data assimilation problems
is given by the _inexact constraint_ preconditioner

$\mathcal{P}_{\mathcal{C}}:=\begin{pmatrix}\mathbf{D}&0&\mathbf{\widehat{L}}\\\
0&\mathbf{R}&0\\\ \mathbf{\widehat{L}}^{T}&0&0\\\ \end{pmatrix},$ (7)

which does not involve the inexact Schur complement $\mathbf{\widehat{S}}$.

Clearly, the effectiveness of the preconditioning operators
$\widehat{\mathbf{S}}$, $\mathcal{P}_{\mathcal{D}}$,
$\mathcal{P}_{\mathcal{T}}$, and $\mathcal{P}_{\mathcal{C}}$ significantly
depends on the adopted approximations $\mathbf{\widehat{L}}$ and
$\mathbf{\widehat{S}}$. In this paper we introduce novel tools which allow for
more successful selections of $\mathbf{\widehat{L}}$ and
$\mathbf{\widehat{S}}$.  In particular, in section 3.1 we propose a novel
approximation $\mathbf{\widehat{L}}\approx\mathbf{L}$ which amounts to a Stein
operator. We will show that the inversion of such $\mathbf{\widehat{L}}$ is
still computationally affordable by exploiting its matrix equation structure,
while it leads to a dramatic decrease in the iteration count. Moreover, we
extend this matrix-oriented method to preconditioning operators
$\mathbf{\widehat{S}}$ which explicitly take into account information from the
observation term $\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$ of the Schur
complement $\mathbf{S}$ (see section 3.2) by adapting a low-rank correction
approach that was proposed in [48]. The original techniques proposed in this
paper lead to the design of preconditioning operators with better theoretical
properties (section 4.1) and more competitive computational records (section
5).

We conclude this section by presenting a novel result related to the use of
$\mathcal{P}_{\mathcal{D}}$ in our setting. We report the proof of the
following theorem in Appendix A.

######  Theorem 1.

If $\mathbf{b}=(b^{T},d^{T},0)^{T}$ denotes the right-hand side in (3), then
the orthonormal basis vectors $\\{v_{1},\ldots,v_{m}\\}$ of the Krylov
subspace
$K_{m}(\mathcal{AP}_{\mathcal{D}}^{-1},\mathbf{b})=\text{span}\\{\mathbf{b},\mathcal{AP}_{\mathcal{D}}^{-1}\mathbf{b},\ldots,(\mathcal{AP}_{\mathcal{D}}^{-1})^{m-1}\mathbf{b}\\}$
computed by GMRES are such that

$v_{2k-1}=\begin{bmatrix}u_{2k-1}\\\ w_{2k-1}\\\
0\end{bmatrix},\quad\text{and}\quad v_{2k}=\begin{bmatrix}0\\\ 0\\\
z_{2k}\end{bmatrix},\quad\text{for any}\;k\geq 1.$

The zero block structure of the basis vectors illustrated in Theorem 1 can be
exploited to design more efficient implementations of the preconditioning step
involving $\mathcal{P}_{\mathcal{D}}$. For instance, we can invert the
(inexact) Schur complement $\widehat{\mathbf{S}}$ only for alternate
iterations. Similarly, the linear systems with $\mathbf{D}$ and $\mathbf{R}$
play a role only in case of an odd iteration index. The GMRES
orthogonalization step can also benefit from Theorem 1, as there is no need to
explicitly perform the orthonormalization of the blocks which necessarly have
to be zero in the current iteration.

We take advantage of these observations to obtain all the results related to
the performance achieved by $\mathcal{P}_{\mathcal{D}}$, which are reported in
section 5.

##  3 A new preconditioning operator

In this section we present the main contribution of this paper. In particular,
we propose to use the following operator

$\widehat{\mathbf{L}}=I_{N+1}\otimes
I_{s}-\Sigma\otimes\widehat{M}=\begin{pmatrix}I&&&\\\ -\widehat{M}&I&&\\\
&\ddots&\ddots&\\\ &&-\widehat{M}&I\\\ \end{pmatrix},$ (8)

in place of $\mathbf{L}$ within the selected preconditioning framework. The
matrix $\widehat{M}$ in (8) is chosen to be some representative value of the
$M_{i}$’s defining $\mathbf{L}$.

In the numerical experiments in section 5 we consider a number of options for
$\widehat{M}$ including, one of the $M_{i}$’s (e.g. the smallest/largest in
norm or condition number, first/last in the sequence), possibly cycling on the
index $i$, and the Karcher matrix mean [5] when the $M_{i}$’s are all
SPD333Alternative matrix means can be used depending on the problem at hand..

The employment of the operator $\widehat{\mathbf{L}}$ described in (8) in the
definition of $\widehat{\mathbf{S}}$, $\mathcal{P}_{\mathcal{D}}$,
$\mathcal{P}_{\mathcal{T}}$, and $\mathcal{P}_{\mathcal{C}}$ leads to novel
preconditioning operators for (2) and (3) that can significantly outperform
other state-of-the-art approaches. In the next proposition we provide some
indications of when we might expect using our fresh approach to be
particularly effective. To this end, we present theoretical bounds on the
eigenvalues of
$\mathbf{\widehat{L}}^{-T}\mathbf{L}^{T}\mathbf{L}\mathbf{\widehat{L}}^{-1}$.

######  Proposition 2.

Let $D_{i}=\widehat{M}-M_{i}$ and $\widehat{\mathbf{L}}$ as in (8). The
eigenvalues of
$\widehat{\mathbf{L}}^{-T}\mathbf{L}^{T}\mathbf{L}\widehat{\mathbf{L}}^{-1}$
can be bounded above by

$1+\frac{N}{2}\left(\rho_{N}+\sqrt{\rho_{N}^{2}+4\rho_{N}}\right),$ (9)

where

$\rho_{N}=\left\\{\begin{array}[]{ll}N\cdot\max_{m=1,\dots,N}\lambda_{max}(D_{m}^{T}D_{m}),&\text{if}\;\lambda_{\max}(\widehat{M}^{T}\widehat{M})=1,\\\
&\\\
\frac{1-\lambda_{\max}^{N}(\widehat{M}^{T}\widehat{M})}{1-\lambda_{\max}(\widehat{M}^{T}\widehat{M})}\cdot\max_{m=1,\dots,N}\lambda_{max}(D_{m}^{T}D_{m}),&\text{otherwise}.\end{array}\right.$
(10)

Due to the multiple levels of approximation used to obtain the result of
Proposition 2, the bounds are likely to be loose in practice. However, the
qualitative information encoded in (9) may provide a way to select a ‘good’
choice of $\widehat{M}$, and an indication of when the preconditioner (8) is
likely to be effective.

In particular, Proposition 2 indicates that the best results are likely to be
obtained when $\max_{i,j}\|M_{i}-M_{j}\|$ is small. If the difference between
the linearised model operators is large then the maximum eigenvalue of the
difference terms cannot all be kept small. Similarly, the spectral norm of
$\widehat{M}$ itself must also be small. If not, then the sum in (32) will
blow up rapidly even for moderate values of $N$. Both of these observations
provide insight into a heuristic way to select $\widehat{M}$: begin by
choosing $M_{i}$ with smallest norm. If all the values of $||M_{i}||$ are
similar, then it is likely that the $D_{i}$ term becomes more important – we
can then choose $\widehat{M}$ to be the value of $M_{i}$ that minimises the
average value of $\|D_{i}D_{i}^{T}\|$.  See section 5 for a panel of diverse
numerical experiments displaying such trends.

To obtain a successful preconditioning strategy, operating with
$\widehat{\mathbf{L}}$ in (8) must not be computationally demanding. In
particular, the application of the preconditioners $\widehat{\mathbf{S}}$,
$\mathcal{P}_{\mathcal{D}}$, $\mathcal{P}_{\mathcal{T}}$, and
$\mathcal{P}_{\mathcal{C}}$ always requires the inversion of
$\widehat{\mathbf{L}}$. The efficient computation of
$\widehat{\mathbf{L}}^{-1}$ will be the subject of the next section.

###  3.1 On the inversion of the Stein operator $\widehat{\mathbf{L}}$

Thanks to the properties of the Kronecker product, see, e.g., [42], the action
of $\widehat{\mathbf{L}}$ in (8) on a vector $z=\text{vec}(Z)$ can be written
as follows

$\widehat{\mathbf{L}}z=\text{vec}(Z-\widehat{M}Z\Sigma^{T}).$

A linear operator of the form

$\begin{array}[]{rrll}\mathfrak{L}:&\mathbb{R}^{s\times(N+1)}&\rightarrow&\mathbb{R}^{s\times(N+1)}\\\
&Z&\mapsto&Z-\widehat{M}Z\Sigma^{T}\\\ \end{array}$

is called a Stein operator in the matrix equation literature; see, e.g., [42].
Therefore, $\widehat{\mathbf{L}}z=\text{vec}(\mathfrak{L}(Z))$. Due to this
relation, hereafter, with abuse of notation, we say that also
$\widehat{\mathbf{L}}$ amounts to a Stein operator. The inversion of
$\widehat{\mathbf{L}}$ is thus equivalent to inverting $\mathfrak{L}$, and
hence to solving a so-called Stein matrix equation

$\text{vec}(Z)=\widehat{\mathbf{L}}^{-1}\text{vec}(V)\quad\Longleftrightarrow\quad
Z-\widehat{M}Z\Sigma^{T}=V.$ (11)

Similarly,

$\text{vec}(Z)=\widehat{\mathbf{L}}^{-T}\text{vec}(V)\quad\Longleftrightarrow\quad
Z-\widehat{M}^{T}Z\Sigma=V.$ (12)

Different numerical methods have been proposed in the literature for the
efficient solution of Stein matrix equations. See, e.g., [1, 25] and [42,
Section 6].

In our setting, we need to solve equations (11) and (12) several times.
Indeed, depending on the adopted preconditioning scheme, a couple of Stein
equations have to be solved at each GMRES/CG iteration. By fully exploiting
the structure of the coefficient matrices defining the Stein equations, we
illustrate a novel solution procedure that remarkably reduces the
computational cost of the preconditioning steps. Our original scheme requires
some minor precomputation that can be performed once prior to the start of the
adopted Krylov iterative scheme.

We first notice that we can write

$\Sigma=C-e_{1}e_{N+1}^{T},\qquad C=\begin{pmatrix}0&&&1\\\ 1&0&&\\\
&\ddots&\ddots&\\\ &&1&0\\\ \end{pmatrix}.$ (13)

Thanks to its circulant structure, $C$ can be cheaply diagonalized by the fast
Fourier transform (FFT), namely $C=F^{-1}\Pi F$ where

$\Pi=\text{diag}(\pi_{1},\ldots,\pi_{N+1}),\quad(\pi_{1},\ldots,\pi_{N+1})^{T}=FCe_{1},$

and $F$ denotes the discrete Fourier matrix. The observation in (13) leads to
the following result.

######  Proposition 3.

Let $\widehat{M}=T\Lambda T^{-1}$ be the eigendecomposition of $\widehat{M}$
with $\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{s})$ and
$P\in\mathbb{C}^{s\times(N+1)}$ be such that
$P_{i,j}=1/(1-\lambda_{i}\pi_{j})$. Moreover, let
${\color[rgb]{0,0,0}U}=I+\text{diag}((P(\Lambda Fe_{1}\circ F^{-T}e_{N+1})))$
where $\circ$ denotes the Hadamard product. Then the solution $Z$ to the Stein
equation in (11) can be written as

$Z=T(Y-W)F^{-T},$ (14)

where

$Y=P\circ(T^{-1}VF^{T}),\quad\text{and}\quad
W=P\circ({\color[rgb]{0,0,0}U}^{-1}(\Lambda YF^{-1}e_{N+1})e_{1}^{T}F^{T}).$

Similarly, the solution $Z$ to (12) is such that

${\color[rgb]{0,0,0}Z}=T^{-T}(G-H)F,$ (15)

where

$G=P\circ(T^{T}{\color[rgb]{0,0,0}V}F^{-1}),\quad\text{and}\quad
H=P\circ({\color[rgb]{0,0,0}U}^{-1}(\Lambda GF^{T}e_{1})e_{N+1}^{T}F^{-1}).$

The computational cost of the solution of the Stein equations (11)–(12) by
(14)–(15) amounts to $\mathcal{O}(s^{3}(N+1)\log(N+1))$ floating point
operations: the cubic term $s^{3}$ arises from the eigendecomposition of
$\widehat{M}$,  while the use of the FFT, namely computing the action of the
matrices $F$ and $F^{-1}$ in (14)–(15), leads to the polylogarithmic term in
$N+1$. Even though the eigendecomposition of $\widehat{M}$ can be computed
once, prior to the start of the Krylov routine, the approach presented in
Proposition 3 requires the matrix $\widehat{M}$ to be of moderate size. On the
other hand, by fully exploiting the circulant-plus-low-rank structure of
$\Sigma$, we can afford sizable values of $N$. See [33, Section 5] for
constructions similar to the ones stated in Proposition 3 derived for the
solution of certain Sylvester equations.

The procedures for solving this Stein equation and its transpose are
summarized in Algorithm 1 and 2, respectively, and they rely on the results
presented in Proposition 3. Notice that only matrix-matrix multiplications and
entry-wise operations are performed in Algorithm 1 and 2 making the
preconditioning step easy to parallelize. See, e.g., [19].

input :

$P\in\mathbb{C}^{s\times(N+1)}$ and
$\Lambda,T,{\color[rgb]{0,0,0}U}\in\mathbb{C}^{s\times s}$ as in Proposition
3, $V\in\mathbb{R}^{s\times(N+1)}$.

output :

$Z\in\mathbb{R}^{s\times(N+1)}$ solution to $Z-\widehat{M}Z\Sigma^{T}=V$.

Compute $Y=P\circ(T^{-1}VF^{T})$

Compute $W=P\circ({\color[rgb]{0,0,0}U}^{-1}(\Lambda
YF^{-1}e_{N+1})e_{1}^{T}F^{T})$

Set $Z=T^{-1}(Y-W)F^{-T}$

Algorithm 1 Solution of the Stein equation $Z-\widehat{M}Z\Sigma^{T}=V$.

input :

$P\in\mathbb{C}^{s\times(N+1)}$ and
$\Lambda,T,{\color[rgb]{0,0,0}U}\in\mathbb{C}^{s\times s}$ as in Proposition
3, ${\color[rgb]{0,0,0}V\in\mathbb{R}^{s\times(N+1)}}$.

output :

$Z\in\mathbb{R}^{s\times(N+1)}$ solution to $Z-\widehat{M}^{T}Z\Sigma=V$.

Compute $G=P\circ(T^{T}VF^{-1})$

Compute $H=P\circ({\color[rgb]{0,0,0}U}^{-1}(\Lambda
GF^{T}e_{1})e_{N+1}^{T}F^{-1})$

Set $Z=T^{-T}(G-H)F$

Algorithm 2 Solution of the Stein equation $X-\widehat{M}^{T}X\Sigma=U$.

###  3.2 Influence of Schur complement approximations

The quality of the Schur complement approximation $\widehat{\mathbf{S}}$ plays
an important role in determining the effectiveness of the preconditioning
operators for (2) and (3). In this work we make use of two choices of
$\widehat{\mathbf{S}}$. We briefly consider classic approximations of the form
$\widehat{\mathbf{S}}=\widehat{\mathbf{L}}^{\top}\mathbf{D}^{-1}\widehat{\mathbf{L}}$
as studied in [14, 17] but involving the new approach for
$\widehat{\mathbf{L}}$ illustrated in the previous section. The second option
is motivated by a low-rank approximation proposed in [48] and includes
information from the observation term explicitly. If

$\widehat{\mathbf{S}}={\mathbf{L}}^{\top}\mathbf{D}^{-1}{\mathbf{L}},$ (16)

then the eigenvalues of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ are given by
$(N+1)(s-p)$ unit eigenvalues, and $p(N+1)$ eigenvalues given by
$1+\lambda({\mathbf{L}}^{-1}\mathbf{D{L}}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H})$.
We note that
${\mathbf{L}}^{-1}\mathbf{D{L}}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$
is of rank $p(N+1)$ with non-negative eigenvalues, meaning that the minimum
eigenvalue of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ is 1. However, the
remaining non-unit eigenvalues can be large, for example in the case that
$\mathbf{R}$ is ill-conditioned; see, e.g., [46]. An alternative choice of
$\widehat{\mathbf{S}}$ which allows any extreme eigenvalues arising from the
observation term to be accounted for within the preconditioner, comes from
considering a low-rank update to (16) of the form

$\widehat{\mathbf{S}}={\mathbf{L}}^{\top}\mathbf{D}^{-1}{\mathbf{L}}+\mathbf{K}_{r}\mathbf{K}_{r}^{T},$
(17)

where $\mathbf{K}_{r}=I_{N+1}\otimes
V_{r}\Upsilon_{r}^{1/2}\in\mathbb{R}^{(N+1)s\times(N+1)r}$ and
$V_{r}\Upsilon_{r}^{1/2}\in\mathbb{R}^{s\times r}$ is constructed from the
leading $r$ terms of the eigendecomposition $V\Upsilon V^{T}=H^{T}R^{-1}H$.
This approach has been studied theoretically in [48], where it was proved that
in addition to increasing the number of unit eigenvalues, increasing $r$
reduces the largest eigenvalues of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$.

A similar low-rank update approach can also be considered when using an
approximation $\widehat{\mathbf{L}}$ to $\mathbf{L}$. In this setting the
smallest eigenvalue of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ can now be
smaller than one, and including more information from the observation term is
not guaranteed to reduce bounds on the largest eigenvalue of the
preconditioned system. However, this approach has been found to perform well
for a number of problems, particularly where $\widehat{\mathbf{L}}$ is a
spectrally good approximation to $\mathbf{L}$. In what follows we apply the
low-rank update to an approximate first term.

In a true low-rank approach ($r<<p$) computational efficiency is ensured by
applying the Woodbury identity. This avoids applying the inverse of
$\mathbf{D}$, which is expensive and allows the re-use of parallelisable or
inexpensive approximations to $\mathbf{\widehat{L}}^{-1}$. As the setting
considered in this paper requires blocks that are not too large
($\mathcal{O}(10^{3})$), it is not unreasonable to compute a full
decomposition of $H^{T}R^{-1}H$. We therefore propose using the low-rank
approach and Woodbury implementation with large values of $r\leq p$, i.e.,

$\widehat{\mathbf{S}}^{-1}=\widehat{\mathbf{L}}^{-1}\mathbf{D}\widehat{\mathbf{L}}^{-T}-\widehat{\mathbf{L}}^{-1}\mathbf{D}\widehat{\mathbf{L}}^{-T}\mathbf{K}_{r}\left(I_{r(N+1)}+\mathbf{K}_{r}^{T}\widehat{\mathbf{L}}^{-1}\mathbf{D}\widehat{\mathbf{L}}^{-T}\mathbf{K}_{r}\right)^{-1}\mathbf{K}_{r}^{T}\widehat{\mathbf{L}}^{-1}\mathbf{D}\widehat{\mathbf{L}}^{-T}.$
(18)

####  3.2.1 Algorithmic considerations

While the use of  the preconditioner $\widehat{\mathbf{S}}^{-1}$ in (18) leads
to great gains in the convergence properties of the selected preconditioned
iterative scheme, especially for $r\approx p$ – see section 5 – it also poses
some computational challenges. We now demonstrate how to implement (18) in a
feasible way by exploiting the Kronecker structure of the new preconditioner.

One benefit of using the Woodbury formulation is that the efficient
implementations of $\widehat{\mathbf{L}}^{-1}\text{vec}(V)$ and
$\widehat{\mathbf{L}}^{-T}\text{vec}(Z)$ that were introduced in section 3.1
can be reused to apply $\hat{\mathbf{L}}^{-1}\mathbf{D}\hat{\mathbf{L}}^{-T}$.
Similarly, we can exploit the Kronecker structure of
$\mathbf{K}_{r}=I_{N+1}\otimes V_{r}\Upsilon_{r}^{1/2}$ to cheaply apply this
operator and its transpose.

Therefore the main computational bottleneck of the relation (18) is the
solution of the $r(N+1)\times r(N+1)$ linear system with
$I_{r(N+1)}+\mathbf{K}_{r}^{T}\widehat{\mathbf{L}}^{-1}\mathbf{D}\widehat{\mathbf{L}}^{-T}\mathbf{K}_{r}$.
We obtain computational gains by first transforming this problem into an
equivalent one, and then solving the transformed problem iteratively using an
inner matrix-oriented CG problem.

Solving a transformed problem:
We can make considerable computational savings by solving a transformed
problem that exploits the identity plus Kronecker structure of
$\hat{\mathbf{L}}$. Writing  $\widehat{M}=T\Lambda T^{-1}$, we define

$\widetilde{\mathbf{L}}=I_{N+1}\otimes
I_{s}-\Sigma\otimes\Lambda,\quad\widetilde{\mathbf{S}}=\widetilde{\mathbf{L}}^{\top}\widetilde{\mathbf{D}}^{-1}\widetilde{\mathbf{L}}+\widetilde{\mathbf{K}}_{r}\widetilde{\mathbf{K}}_{r}^{T},$

where $\mathbf{\widetilde{D}}=e_{1}e_{1}^{T}\otimes
T^{-1}BT^{-T}+(I_{N+1}-e_{1}e_{1}^{T})\otimes T^{-1}QT^{-T}$, and
$\mathbf{\widetilde{K}}_{r}=I_{N+1}\otimes V_{r}\Upsilon_{r}^{1/2}$ with
$V_{r}$, $\Upsilon_{r}$ coming now from the eigendecomposition of
$\widetilde{H}^{T}R^{-1}\widetilde{H}$, $\widetilde{H}=HT$. The computation of
$\widetilde{\mathbf{S}}^{-1}$ now involves the solution of a linear system
with
$I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\mathbf{\widetilde{L}}^{-1}\mathbf{\widetilde{D}}\mathbf{\widetilde{L}}^{-T}\mathbf{\widetilde{K}}_{r}$
whose action is performed by following Algorithm 3. Notice that the cost of
the latter algorithm is now linear in $s$ and polylogarithmic in $N+1$ thanks
to the semidiagonalization of $\widehat{\mathbf{L}}$.

For the SPD problem (2), we apply the preconditioner
$\widehat{\mathbf{S}}=(I_{N+1}\otimes
T^{-T})\widetilde{\mathbf{S}}(I_{N+1}\otimes T^{-1})$ as

$\widehat{\mathbf{S}}^{-1}\text{vec}(V)=(I_{N+1}\otimes
T)(\widetilde{\mathbf{S}}^{-1}\text{vec}(T^{T}V)).$ (19)

We note the equality here, and that the only assumption required is that the
full eigendecomposition of $\widehat{M}$ is available.

Similarly, for the saddle-point linear system (3), we still write
$\mathcal{P}_{\mathcal{D}}=\mathcal{T}\mathcal{\widetilde{P}}_{\mathcal{D}}\mathcal{T}^{T}$
and
$\mathcal{P}_{\mathcal{T}}=\mathcal{T}\mathcal{\widetilde{P}}_{\mathcal{T}}\mathcal{T}^{T}$
where

$\mathcal{\widetilde{P}}_{\mathcal{D}}=\begin{pmatrix}\mathbf{\widetilde{D}}&&\\\
&\mathbf{R}&\\\ &&\mathbf{\widetilde{S}}\\\
\end{pmatrix},\qquad\mathcal{\widetilde{P}}_{\mathcal{T}}=\begin{pmatrix}\mathbf{\widetilde{D}}&0&\mathbf{\widetilde{L}}\\\
&\mathbf{R}&\mathbf{\widetilde{H}}\\\ &&\mathbf{\widetilde{S}}\\\
\end{pmatrix},$

and

$\mathcal{T}=\begin{pmatrix}I_{N+1}\otimes T&&\\\ &I_{N+1}\otimes I_{p}&\\\
&&I_{N+1}\otimes T^{-T}\\\ \end{pmatrix}.$ (20)

At the $j$th GMRES iteration we perform

$\mathcal{P}_{\mathcal{D}}^{-1}\text{vec}(V)=\mathcal{T}^{-T}(\mathcal{\widetilde{P}}_{\mathcal{D}}^{-1}(\mathcal{T}^{-T}\text{vec}(V))),$
(21)

and $\mathcal{\widetilde{P}}_{\mathcal{D}}^{-1}$ is computed by following the
strategy presented in the previous sections. Notice that $\mathcal{T}$ is
block diagonal with blocks having a Kronecker form. This rich structure can be
exploited to cheaply perform the transformations involving $\mathcal{T}$
itself. The same approach is adopted for the block triangular preconditioner
$\mathcal{P}_{\mathcal{T}}$. The inexact constraint preconditioner
$\mathcal{P}_{\mathcal{C}}$ would not benefit from the semi-diagonalization of
$\mathbf{\widehat{L}}$ as its definition does not include the approximate
Schur complement $\mathbf{\widehat{S}}$. We note at this stage that the
transformed preconditioner is equivalent to (18).

Inner matrix-oriented CG:
It now remains to solve
$(I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\mathbf{\widetilde{L}}^{-1}\mathbf{\widetilde{D}}\mathbf{\widetilde{L}}^{-T}\mathbf{\widetilde{K}}_{r})^{-1}x$
efficiently. A naive strategy would consist of assembling the coefficient
matrix first, by possibly exploiting the Kronecker structure of the involved
factors $\mathbf{\widetilde{K}}_{r}$, $\widetilde{\mathbf{L}}$, and
$\mathbf{\widetilde{D}}$. To this end, the scheme proposed in [20] could be
employed with straightforward modifications. Even though this step can be
carried out prior to the GMRES/CG iterations, the construction of the dense
matrix
$\mathbf{\widetilde{K}}_{r}^{T}{\mathbf{\widetilde{L}}}^{-1}\mathbf{\widetilde{D}}{\mathbf{\widetilde{L}}}^{-T}\mathbf{\widetilde{K}}_{r}$
requires the solution of $r(N+1)$ matrix equations making this task
computationally unaffordable.

We thus pursue a different path. Since $\mathbf{D}$ and $\mathbf{R}$ are SPD
by construction
$I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\widetilde{\mathbf{L}}^{-1}\mathbf{\widetilde{D}}\widetilde{\mathbf{L}}^{-T}\mathbf{\widetilde{K}}_{r}$
is SPD as well. Moreover, the Kronecker structure of the latter matrix can be
exploited to cheaply compute its action as mentioned above. We therefore
propose using an iterative method to approximate the solution of
$(I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\mathbf{\widetilde{L}}^{-1}\mathbf{\widetilde{D}}\mathbf{\widetilde{L}}^{-T}\mathbf{\widetilde{K}}_{r})^{-1}x$
within the preconditioner by means of a matrix-oriented CG method. This
iterative method only requires only applications of the operator
$I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\widetilde{\mathbf{L}}^{-1}\mathbf{\widetilde{D}}\widetilde{\mathbf{L}}^{-T}\mathbf{\widetilde{K}}_{r}$.
Therefore, the overall scheme for solving (2) and (3) can be seen as an inner-
outer iteration [43] whenever the approximation (18) with $r>0$ is adopted
within the selected preconditioning framework. In particular, the outer Krylov
routine (GMRES/CG) is preconditioned with a scheme involving a second, inner
Krylov method (CG). Notice that the use of CG within the preconditioning step
requires the employment of a flexible variant of the outer Krylov method as we
are using different approximate preconditioners for each outer iteration; see
[43]. The matrix-oriented implementation of flexible GMRES and CG can be
easily obtained from their standard form [38, 31].

We stress once again that combining the inner matrix-oriented CG method with
the semi-diagonalised approach to solve for $\widetilde{\mathbf{S}}^{-1}$
significantly lowers the computational cost of the preconditioning step,
especially for the case $r>0$. Indeed, in this case, the cost of the CG
iterations involved in the computation of $\widetilde{\mathbf{S}}^{-1}$ is
linear in $s$ and polylogarithmic in $N$ thanks to the semi-diagonalization of
$\mathbf{\widehat{L}}$; see Algorithm 3. We note that by using an inner
iterative solver we obtain an approximation to the ‘true’ preconditioner
$\widetilde{\mathbf{S}}$. However, our numerical experiments in section 5
reveal that we can obtain near optimal performance using this nested approach
for a reasonable choice of tolerance within the inner matrix-oriented CG
problem.

input :

$Z\in\mathbb{C}^{r\times(N+1)}$, $V_{r}$, $\Upsilon_{r}\in\mathbb{R}^{s\times
r}$, and $\Lambda$, $\widetilde{B}$, $\widetilde{Q}\in\mathbb{R}^{s\times s}$.

output :

$X\in\mathbb{R}^{r\times(N+1)}$ such that
$\text{vec}(X)=(I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\mathbf{\widetilde{L}}^{-1}\mathbf{\widetilde{D}}\mathbf{\widetilde{L}}^{-T}\mathbf{\widetilde{K}}_{r})\text{vec}(Z)$.

Solve

$Y-\Lambda Y\Sigma=V_{r}\Upsilon_{r}Z$

by means of Algorithm 2 with $T=I$

Compute
$W=\widetilde{B}Ye_{1}e_{1}^{T}+\widetilde{Q}Y(I_{N+1}-e_{1}e_{1}^{T})$

Solve

$U-\Lambda U\Sigma^{T}=W$

by means of Algorithm 1 with $T=I$

Set $X=Z+\Upsilon_{r}V_{r}^{T}U$

Algorithm 3 Computing the action of
$I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\mathbf{\widetilde{L}}^{-1}\mathbf{\widetilde{D}}\mathbf{\widetilde{L}}^{-T}\mathbf{\widetilde{K}}_{r}$.

As previously mentioned, matrix CG is often equipped with some low-rank
truncations to reduce the overall memory demand; see, e.g., [27]. However,
storage will not be an issue in our context thanks to the modest problem
dimensions we consider. Moreover, the introduction of any low-rank truncation
would worsen the performance of the preconditioning step in general.
Therefore, we perform no low-rank truncations within the inner matrix CG.

We would like to mention that, similarly to the outer Krylov routine, the
inner matrix CG can also be preconditioned to achieve a faster convergence in
terms of the number of iterations. However, we were not able to design an
effective preconditioning operator for
$I_{r(N+1)}+\mathbf{\widetilde{K}}_{r}^{T}\widetilde{\mathbf{L}}^{-1}\mathbf{\widetilde{D}}\widetilde{\mathbf{L}}^{-T}\mathbf{\widetilde{K}}_{r}$
with a reasonable computational cost. A number of natural preconditioning
strategies did not yield improvement in terms of computational speed compared
to a plain, unpreconditioned matrix CG implementation. We therefore present an
unpreconditioned inner matrix CG in all the numerical experiments reported in
section 5.

One may also consider performing the FFT transformations involved in the
solution of the Stein equations outside the GMRES/CG iteration, thinking that
this would further decrease the computational cost of the preconditioning
steps. However, this would also introduce complex arithmetic in the GMRES/CG
iteration, increasing the cost of the overall scheme. Moreover, the
application of the FFT can be cheaply performed without forming the discrete
Fourier matrix $F$. In particular, in all our numerical tests we employed the
Matlab fft and ifft functions. In light of these considerations, we confine
the use of FFT to the preconditioning step only.

##  4 Observation-time independent $\mathcal{M}$

If the model forecast $\mathcal{M}$ takes a constant value between each
observation time, the same is true for its linearization. In this case, all
the matrices $M_{i}$ in the definition of $\mathbf{L}$ are the same, namely
$M=M_{i}$ for all $i=1,\ldots,N+1$. In this case, the operator

$\mathbf{L}=\begin{pmatrix}I&&&\\\ -M&I&&\\\ &\ddots&\ddots&\\\ &&-M&I\\\
\end{pmatrix}=I_{N+1}\otimes I_{s}-\Sigma\otimes M,$ (22)

is a Stein operator itself and we can thus use
$\widehat{\mathbf{L}}=\mathbf{L}$ in our preconditioning strategy. In this
easier setting, the latter choice leads to narrow eigenvalue distributions of
the preconditioned coefficient matrices – see section 4.1 – while maintaining
high computational efficiency.

In contrast to what happens in the more general case discussed at the end of
section 3.2.1, we can now perform the transformation based on the eigenvector
matrix $T$ once before the Krylov routine starts, and not every time the
preconditioner is applied. For instance, if $M=T\Lambda T^{-1}$, (2) can be
written as

$\widetilde{\mathbf{S}}\widetilde{\delta x}=f,$ (23)

where
$\widetilde{\mathbf{S}}=\widetilde{\mathbf{L}}^{T}\widetilde{\mathbf{D}}^{-1}\widetilde{\mathbf{L}}+\widetilde{\mathbf{H}}^{T}\mathbf{R}\widetilde{\mathbf{H}}$,
$\widetilde{\mathbf{D}}$, $\widetilde{\mathbf{H}})$, and
$\widetilde{\mathbf{L}}$ are as in section 3.2.1, and $f=(I\otimes
T^{T})(\mathbf{D}^{-1}b+\mathbf{L}^{T}\mathbf{H}^{T}\mathbf{R}^{-1}d)$. Once
$\widetilde{\delta x}$ is computed, we retrieve the actual solution by
performing $\delta x=(I_{N+1}\otimes T)\widetilde{\delta x}$.

The same approach can be followed for the saddle point linear system (3). We
can write

$\displaystyle\mathcal{A}=$ $\displaystyle\begin{pmatrix}e_{1}e_{1}^{T}\otimes
B+(I_{N+1}-e_{1}e_{1}^{T})\otimes Q&&I_{N+1}\otimes I_{s}-\Sigma\otimes M\\\
0&I_{N+1}\otimes R&I_{N+1}\otimes H\\\ I_{N+1}\otimes I_{s}-\Sigma^{T}\otimes
M^{T}&I_{N+1}\otimes H^{T}&0\end{pmatrix}$ $\displaystyle=$
$\displaystyle\mathcal{T}\underbrace{\begin{pmatrix}e_{1}e_{1}^{T}\otimes\widetilde{B}+(I_{N+1}-e_{1}e_{1}^{T})\otimes\widetilde{Q}&&I_{N+1}\otimes
I_{s}-\Sigma\otimes\Lambda\\\ 0&I_{N+1}\otimes
R&I_{N+1}\otimes\widetilde{H}\\\ I_{N+1}\otimes
I_{s}-\Sigma^{T}\otimes\Lambda&I_{N+1}\otimes\widetilde{H}^{T}&0\\\
\end{pmatrix}}_{\mathcal{\widetilde{A}}}\mathcal{T}^{T},$

where $\widetilde{B}$, $\widetilde{Q}$, $\widetilde{H}$, and $\mathcal{T}$ are
as in section 3.2.1. In place of (3) we can thus solve the transformed system

$\mathcal{\widetilde{A}}\begin{pmatrix}\widetilde{\delta\eta}\\\
\widetilde{\delta\lambda}\\\ \widetilde{\delta x}\\\
\end{pmatrix}=\begin{pmatrix}\widetilde{b}\\\ \widetilde{d}\\\ 0\\\
\end{pmatrix},$ (24)

where

$\begin{pmatrix}\widetilde{\delta\eta}\\\ \widetilde{\delta\lambda}\\\
\widetilde{\delta x}\\\
\end{pmatrix}=\mathcal{T}^{T}\begin{pmatrix}\delta\eta\\\ \delta\lambda\\\
\delta x\\\ \end{pmatrix},\quad\text{and}\quad\begin{pmatrix}\widetilde{b}\\\
\widetilde{d}\\\ 0\\\ \end{pmatrix}=\mathcal{T}^{-1}\begin{pmatrix}b\\\ d\\\
0\\\ \end{pmatrix}.$

Once $(\widetilde{\delta\eta},\widetilde{\delta\lambda},\widetilde{\delta
x})^{T}$ is computed, we retrieve the original solution by
$(\delta\eta,\delta\lambda,\delta
x)^{T}=\mathcal{T}^{-T}(\widetilde{\delta\eta},\widetilde{\delta\lambda},\widetilde{\delta
x})^{T}$.

The preconditioning operators for (23) and (24) can be obtained by mimicking
what we presented in the previous sections. The major difference is the
cheaper inversion of the Stein operator $\widetilde{\mathbf{L}}=I_{N+1}\otimes
I_{s}-\Sigma\otimes\Lambda$ which is now “semi”-diagonalized. The cost of
computing $\widetilde{\mathbf{L}}^{-1}$ is thus linear in $s$ and
polylogarithmic in $N+1$.

###  4.1 Spectral results

In this section we present bounds on the eigenvalues of the preconditioned
systems using our new approach when $\widehat{\mathbf{L}}=\mathbf{L}$.

The spectral properties of linearised data assimilation problem (2) have been
studied in [45] for the unpreconditioned 3D-Var formulation, and in [46] for
strong-constraint 4D-Var preconditioned with the exact first term. Bounds on
the spectrum of the (preconditioned) Hessian $\mathbf{S}$ for the weak-
constraint problem can be obtained using the same theoretical approaches and
replacing $\mathbf{B}$ in those bounds with
$\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}$.

######  Proposition 4.

Let $\widehat{\mathbf{S}}=\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}$. Then
$\widehat{\mathbf{S}}^{-1}\mathbf{S}=\mathbf{I}_{s(N+1)}+\mathbf{L}^{-1}\mathbf{D}\mathbf{L}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$
and the eigenvalues of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ are bounded as
follows

$\lambda(\widehat{\mathbf{S}}^{-1}\mathbf{S})\in\left[1,\frac{\lambda_{\max}(\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H})}{\lambda_{\min}(\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L})}\right].$
(25)

We note that $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ has $(s-p)(N+1)$ unit
eigenvalues.

###### Proof.

Apply [46, Theorem 4] replacing $\mathbf{B}^{-1}$ with
$\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}$. ∎

In the case where the low-rank update to the Schur complement preconditioner
presented in section 3.2 is applied with $\widehat{\mathbf{L}}=\mathbf{L}$, a
bound on the maximum eigenvalue is controlled by the largest neglected
eigenvalue that is not included in the approximation $\mathbf{K}_{r}$.

######  Proposition 5 ([48]).

Let
$\widehat{\mathbf{S}}=\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}+\mathbf{K}_{r}\mathbf{K}_{r}^{T}$
with $\mathbf{K}_{r}$ defined as in section 3.2. Then
${\mathbf{S}}^{-1}\mathbf{S}=\mathbf{I}_{s(N+1)}+(\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}+\mathbf{K}_{r}\mathbf{K}_{r}^{T})^{-1}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$
and the eigenvalues of $\widehat{\mathbf{S}}^{-1}\mathbf{S}$ are bounded
between

$\lambda(\widehat{\mathbf{S}}^{-1}\mathbf{S})\in\left[1,\frac{\lambda_{r+1}}{\lambda_{\min}(\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L})}\right].$
(26)

where $\lambda_{r+1}$ is the $(r+1)$th largest eigenvalue of
$\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$, i.e. the largest eigenvalue that is
neglected by the low-rank approximation $\mathbf{K}_{r}\mathbf{K}_{r}^{T}$ to
$\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}$. We note that
$\widehat{\mathbf{S}}^{-1}\mathbf{S}$ has $(s+r-p)(N+1)$ unit eigenvalues.

######  Corollary 6.

If $r=p$ then $\widehat{\mathbf{S}}^{-1}\mathbf{S}=\mathbf{I}_{n}$ and the
eigenvalues of the preconditioned system are all units.

Approximations of $\widehat{\mathbf{S}}$, either by using randomised
approximations of $\mathbf{K}_{r}$ as proposed in [48], or inner iterative
methods to compute (18) as proposed in this paper may lead to eigenvalues of
the preconditioned system that are smaller than $1$ or larger than the
theoretical upper bound given by Proposition 5. We will study the performance
of these approximate preconditioners in section 5.

The spectral properties of the preconditioned saddle point problem (3) have
been studied in [18, 14, 48, 11], although typically by considering
approximations $\widehat{\mathbf{R}},\widehat{\mathbf{D}}$ and
$\widehat{\mathbf{L}}$ rather than using the exact forward model matrices. In
this work we instead propose using the exact covariance matrices within the
preconditioner, i.e. $\widehat{\mathbf{R}}=\mathbf{R}$ and
$\widehat{\mathbf{D}}=\mathbf{D}$. This is possible due to the exploitation of
the Kronecker structure and the use of matrix iterative methods.

Let
$\lambda(\widehat{\mathbf{S}}^{-1}\mathbf{S})\in[\lambda_{\mathbf{S}},\Lambda_{\mathbf{S}}]$.
In what follows we consider how each of the three preconditioners for the
saddle point problem (3) are affected by the approximation of
$\widehat{\mathbf{S}}$ to $\mathbf{S}$. We now state bounds on the eigenvalues
of the preconditioned saddle point system using each of the preconditioners
introduced in section 2.2 with $\mathbf{\widehat{L}}=\mathbf{L}$.

######  Proposition 7.

With the definitions as stated above, the eigenvalues of
$\mathcal{P}_{\mathcal{D}}^{-1}\mathcal{A}$ are real, and satisfy:

$\displaystyle\lambda(\mathcal{P}_{D}^{-1}\mathcal{A})\in$
$\displaystyle\left[\frac{1-\sqrt{1+4\Lambda_{\mathbf{S}}}}{2},\frac{1-\sqrt{1+{4\lambda_{\mathbf{S}}}}}{2}\right]\cup\\{1\\}\cup\left[\frac{1+\sqrt{1+{4\lambda_{\mathbf{S}}}}}{2},\frac{1+\sqrt{1+4\Lambda_{\mathbf{S}}}}{2}\right].$

###### Proof.

The result directly comes from [37, Theorem 4.2.1]. ∎

We can see that obtaining an improved estimate of the Schur complement (in a
spectral sense) will lead to tighter bounds on the eigenvalues of the
preconditioned system when using $\mathcal{P}_{D}^{-1}\mathcal{A}$.

######  Corollary 8.

If $\widehat{\mathbf{S}}=\mathbf{S}$

$\lambda(\mathcal{P}_{D}^{-1}\mathcal{A})\in\left\\{\frac{1-\sqrt{5}}{2},1,\frac{1+\sqrt{5}}{2}\right\\}.$

The quality of the Schur complement approximation also affects the bounds on
the eigenvalues of the block triangular preconditioner.

######  Proposition 9.

With the definitions as stated above, the eigenvalues of
$\mathcal{P}_{\mathcal{T}}^{-1}\mathcal{A}$ are given by $(s+p)(N+1)$ units,
and the remaining $s(N+1)$ eigenvalues are given by the eigenvalues of
$\mathbf{S}\widehat{\mathbf{S}}^{-1}$.

###### Proof.

We consider the product

$\begin{split}\mathcal{A}\mathcal{P}_{\mathcal{T}}^{-1}&=\begin{pmatrix}\mathbf{D}&0&\mathbf{L}\\\
0&\mathbf{R}&\mathbf{H}\\\
\mathbf{L}^{T}&\mathbf{H}^{T}&0\end{pmatrix}\begin{pmatrix}\mathbf{D}^{-1}&0&\mathbf{D}^{-1}\mathbf{L}\mathbf{\widehat{S}}^{-1}\\\
0&\mathbf{R}^{-1}&\mathbf{R}^{-1}\mathbf{H}\mathbf{\widehat{S}}^{-1}\\\
0&0&-\mathbf{S}^{-1}\end{pmatrix}\\\ &=\begin{pmatrix}I&0&0\\\ 0&I&0\\\
\mathbf{L}^{T}\mathbf{D}^{-1}&\mathbf{H}^{T}\mathbf{R}^{-1}&\mathbf{S}\mathbf{\widehat{S}}^{-1}\end{pmatrix}.\end{split}$
(27)

The eigenvalues of $\mathcal{A}\mathcal{P}_{\mathcal{T}}^{-1}$, and by
similarity $\mathcal{P}_{\mathcal{T}}^{-1}\mathcal{A}$ are therefore given by
$1$ and the eigenvalues of $\mathbf{S}\widehat{\mathbf{S}}^{-1}$. ∎

######  Corollary 10.

If $\widehat{\mathbf{S}}=\mathbf{S}$, then
$\lambda(\mathcal{P}_{\mathcal{T}}^{-1}\mathcal{A})=1$.

######  Proposition 11.

With the definitions as stated above, the eigenvalues of
$\mathcal{P}_{\mathcal{C}}^{-1}\mathcal{A}$ consist of $(2s-p)(N+1)$ unit
eigenvalues, with the remaining $2p(N+1)$ eigenvalues given by

$\lambda(\mathcal{P}_{\mathcal{C}}^{-1}\mathcal{A})=1\pm\sqrt{\lambda_{i}(\mathbf{R}^{-1}\mathbf{HL}^{-1}\mathbf{DL}^{-T}\mathbf{H}^{T})}i$
(28)

###### Proof.

The proof comes from [13, Appendix A]. ∎

##  5 Numerical experiments

###  5.1 Experimental framework

In this section we display some numerical results achieved by the novel
preconditioning framework we presented in this paper for the two problems of
interest (2) and (3). Our matrix-oriented strategy is compared to state-of-
the-art vector-oriented approaches designed for linear systems stemming from
data assimilation problems. In particular, we consider the scheme from [47]
where a user-specified parameter $k$ defines the approximation
$\mathbf{\widehat{L}}^{-1}$. This is used within the Hessian/Schur complement
approximation
$\mathbf{\widehat{S}}=\mathbf{\widehat{L}}^{T}\mathbf{D}^{-1}\mathbf{\widehat{L}}$;
see [47, Section 4] for further details444Notice that choosing
$\mathbf{\widehat{L}}=\mathbf{L}$ is equivalent to setting $k=N+1$ in [47]..
We note that the approach of [47] is designed to increase parallelisability of
the application of $\mathbf{\widehat{L}}$ and hence a preconditioner. All
experiments presented in this section are performed in serial; by exploiting
parallel architectures we expect to see large decreases in wallclock times for
the approach of [47] when $k<<N+1$, which is not the case for our new
strategy.

We also compare matrix-oriented CG with the improved Schur complement (as
presented in section 3.2) against the limited memory preconditioner (LMP)
approach of [9], where an alternative identity plus low-rank preconditioner is
applied as a second level preconditioner. This method can only be implemented
in the vectorised setting, and requires the use of
$\widehat{\mathbf{L}}\equiv\mathbf{L}$. Hence, in its current formulation, the
LMP approach cannot be easily parallelised, making comparison of wallclock
times with the matrix-oriented approach more meaningful than for the approach
of [47].

As previously mentioned, we employ a matrix-oriented implementation of CG
(respectively GMRES) whereas the standard _vector_ form of CG (resp. GMRES) is
adopted whenever a preconditioning strategy different from the one introduced
in this paper is considered. This is mainly due to the possibility of using
existing code for the preconditioners in [14, 17, 47]. Indeed, these routines
have been designed for standard GMRES and CG and not for their matrix-oriented
counterpart. Notice, however, that the two GMRES/CG implementations are
equivalent in exact arithmetic as we do not perform any low-rank truncation
within matrix GMRES/CG. On the other hand, the matrix-oriented form of
GMRES/CG may present some computational advantages due to the Kronecker form
of the blocks of the coefficient matrix $\mathcal{A}$ in (3) and $\mathbf{S}$
in (2). See, e.g., Table 5.

In all the results reported here, the algorithms have been stopped as soon as
the iterative method (either in matrix or vector form) relative residual norm
becomes smaller than $10^{-8}$. The same threshold has been used for the inner
CG relative residual norm when we adopt the strategy presented in section
3.2.1.

In what follows, the matrix-oriented implementation of CG (resp. GMRES) will
be denoted by matCG (resp. matGMRES) whereas its standard, vector counterpart
by vecCG (resp vecGMRES).

All the experiments have been run using Matlab (version 2022a) on a machine
with a 1.8GHz Intel quad-core i7 processor with 15GB RAM on an Ubuntu 20.04.2
LTS operating system.

For both the problem settings we considered, we used data assimilation terms
based on those introduced in [47], which we now present briefly. For all
experiments the dimension of the state is $s=1000$ and the number of
observations at each observation time is given by $p=500$. The background
error covariance matrix and model error covariance matrices are produced using
an adapted SOAR correlation function [47, Equation (15)] with parameters
$L_{B}=0.6,L_{Q}=0.75$, $\sigma_{B}=0.5$, $\sigma_{Q}=0.2$, $100$ non-zero
entries per row for $B$ and $120$ for $Q$. The observation error covariance
matrix $R\in\mathbb{R}^{500\times 500}$ is produced using the block approach
of [47]. The observation operator $H\in\mathbb{R}^{500\times 1000}$ has a
single unit entry per row, arranged in ascending column order. Each of these
terms is repeated in a Kronecker structure to obtain
$\mathbf{D}=e_{1}e^{T}_{1}\otimes B+(I_{N+1}-e_{1}e^{T}_{1})\otimes
Q\in\mathbb{R}^{1000(N+1)\times 1000(N+1)},$ $\mathbf{R}=I_{N+1}\otimes
R\in\mathbb{R}^{500(N+1)\times 500(N+1)}$, and $\mathbf{H}=I_{N+1}\otimes
H\in\mathbb{R}^{500(N+1)\times 1000(N+1)}$. We discuss the two classes of
model matrices in the relevant section.

###  5.2 Results for Lorenz96

Our first example is the Lorenz96 problem [28], a nonlinear set of coupled
ODEs that is often used as a data assimilation test problem due to its chaotic
nature.

Consider $s$ equally spaced points on the unit line, e.g. $\Delta
x=\frac{1}{s}$. For $i=1,\dots,s$, we consider

$\frac{dx_{i}}{dt}=(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+8,$

with periodic boundary conditions (i.e. $x_{-1}=x_{s-1},x_{0}=x_{s}$). We
discretize the equations above by using the numerical implementation of [10],
which integrates the model forward in time using a fourth-order Runge-Kutta
scheme. We consider $s=1000$, and unless otherwise mentioned, we use $\Delta
t=10^{-6}$.

As illustrated in Proposition 2, we expect the strategy proposed in section 3
to work well whenever the selected $\widehat{M}$ has small spectral norm and
is such that the matrices $D_{i}=\widehat{M}-M_{i}$ has small spectral norm
for all $i=1,\ldots,N$ as well. In Table 1 we report $\|\widehat{M}\|$ and
$\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ for $N=10$ and different choices of
$\widehat{M}$ varying $\Delta t$. In particular, we consider as $\widehat{M}$
the symmetrised first and the last matrices in the block subdiagonal of
$\mathbf{L}$, namely $Sym(M_{1})=1/2(M_{1}+M_{1}^{T})$ and
$Sym(M_{10})=1/2(M_{10}+M_{10}^{T})$, and the Karcher mean555The Karcher mean
is computed by means of the routine positive_definite_karcher_mean included in
the Matlab toolbox Manopt 6.0 [6]. As optimization procedure, we adopted the
Barzilai-Borwein approach presented in [24] within
positive_definite_karcher_mean. of the symmetric parts of the matrices
$M_{i}$, $i=1,\ldots,10$, namely $Sym(M_{i})=1/2(M_{i}+M_{i}^{T})$, as these
are all SPD. In what follows, we denote the Karcher mean by $Karch(M_{i})$. We
also report the computed upper bound of Proposition 2.

For $\Delta t\leq 1\times 10^{-3}$ the aforementioned selections of
$\widehat{M}$ lead to very similar results. As $\|\widehat{M}\|$ is close to
one and $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ is rather small, the upper
bound on the eigenvalues of
$\widehat{\mathbf{L}}^{-T}\mathbf{L}^{T}\mathbf{L}\widehat{\mathbf{L}}^{-1}$
is also close to $1$. As $\Delta t$ increases, the norm of the linearised
$M_{i}$ operators move further away from $1$, leading to increases in the
upper bound. For the Lorenz96 problem, $\|M_{i}\|$ increases monotonically
with $i$, meaning that for larger choices of $\Delta t$ the norm of
$Sym(M_{10})$ is much larger than $1$. We see that in the case $\Delta
t=10^{-1}$ the choice of $\widehat{M}$ makes a large difference to the upper
bound in Proposition 2. In section 3, we proposed to select $\widehat{M}$ with
the smallest norm in order to minimise the upper bound on the preconditioned
spectrum, or by minimising $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ in the
case that the norms have similar values. This approach is supported by the
results of Table 1, and motivates the selection of $\widehat{M}=Sym(M_{1})$
for this problem for the experiments that follow.

$\Delta t$ | $\widehat{M}$ | $Sym(M_{1})$ | $Sym(M_{10})$ | $Karch(M_{i})$
---|---|---|---|---
$1\times 10^{-6}$ | $\|\widehat{M}\|$ | $1.000$ | $1.0023$ | $1.0011$
| $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ | $1.7641\times 10^{-5}$ | $1.7641\times 10^{-5}$ | $1.3981\times 10^{-5}$
| Upper bound (9) | $1.0012$ | $1.0012$ | $1.0001$
$1\times 10^{-3}$ | $\|\widehat{M}\|$ | $0.9980$ | $1.0110$ | $1.0046$
| $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ | $1.401\times 10^{-2}$ | $1.401\times 10^{-2}$ | $1.009\times 10^{-2}$
| Upper bound (9) | $1.3796$ | $1.3912$ | $1.3261$
$5\times 10^{-2}$ | $\|\widehat{M}\|$ | $0.9980$ | $1.0110$ | $1.0046$
| $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ | $1.401\times 10^{-2}$ | $1.401\times 10^{-2}$ | $1.009\times 10^{-2}$
| Upper bound (9) | $66.9385$ | $11278$ | $320.97$
$1\times 10^{-1}$ | $\|\widehat{M}\|$ | $0.8932$ | $8.2359$ | $1.5188$
| $\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ | $7.9840$ | $10.8191$ | $7.7113$
| Upper bound (9) | $516.7533$ | $2.148\times 10^{10}$ | $8828.74$
Table 1: Example 5.2. $\|\widehat{M}\|$,
$\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})$ and the upper bound (9) of
Proposition 2 for different selection of $\widehat{M}$ and $N=10$.

In Table 2 we report the performance achieved by our matrix-oriented
preconditioning frameworks: $\widehat{\mathbf{S}}$ (with both $r=0$ and $r=p$
in (18)) for (2) and $\mathcal{P}_{\mathcal{D}}$, $\mathcal{P}_{\mathcal{T}}$
(with both $r=0$ and $r=p$ in (18)), and $\mathcal{P}_{\mathcal{C}}$ for (3),
varying $\widehat{M}$. We compare these results with those attained by the
strategy proposed in [47], where $\widehat{\mathbf{L}}$ is chosen as follows:

$\text{the $(i,j)$th block of
}\widehat{\mathbf{L}}=\begin{cases}\mathbf{I},&\quad i=j,\\\ -M_{i},&\quad
i=j\text{ and }i-k\lfloor\frac{i}{k}\rfloor,\\\ 0,&\quad\text{ otherwise.
}\end{cases}$

As this $\widehat{\mathbf{L}}$ does not have the form of a Stein equation, it
is applied using vecCG/vecGMRES, with the parameter $k=3.$

| $\widehat{\mathbf{S}}$ | $\mathcal{P}_{\mathcal{D}}$ | $\mathcal{P}_{\mathcal{T}}$ | $\mathcal{P}_{\mathcal{C}}$
---|---|---|---|---
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=0$ | 25.7 | 45.0 | 29.0 | 46.0
matCG/matGMRES, $\widehat{M}=Sym(M_{10})$, $r=0$ | 25.8 | 45.0 | 29.0 | 46.0
matCG/matGMRES, $\widehat{M}=karch(M_{i})$, $r=0$ | 26.3 | 45.0 | 28.1 | 46.0
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=p$ | 3.8 | 7.0 | 6.0 | -
matCG/matGMRES, $\widehat{M}=Sym(M_{10})$, $r=p$ | 3.2 | 7.0 | 6.0 | -
matCG/matGMRES, $\widehat{M}=karch(M_{i})$, $r=p$ | 3.5 | 7.0 | 6.0 | -
vecCG/vecGMRES [47], $k=3$ | 529.5 | 358.0 | 188.7 | 66.7
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=0$ | 9.3816 | 7.5924 | 4.9876 | 7.0108
matCG/matGMRES, $\widehat{M}=Sym(M_{10})$, $r=0$ | 9.3436 | 7.7803 | 5.0595 | 6.9689
matCG/matGMRES, $\widehat{M}=karch(M_{i})$, $r=0$ | 9.4084 | 7.6916 | 4.8374 | 6.8653
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=p$ | 3.3120 | 3.2902 | 3.6184 | -
matCG/matGMRES, $\widehat{M}=Sym(M_{10})$, $r=p$ | 3.0478 | 3.5083 | 3.7205 | -
matCG/matGMRES, $\widehat{M}=karch(M_{i})$, $r=p$ | 3.4022 | 3.4488 | 3.7182 | -
vecCG/vecGMRES [47], $k=3$ | 222.9978 | 153.9620 | 91.1391 | 19.7812
Table 2: Example 5.2. Iterations (top) and wallclock time (bottom) to
convergence for the Lorenz96 problem with $\Delta t=10^{-6}$ for $N=10$ for
different preconditioners. Values are averaged over 10 realisations. The pre-
computation for the Karcher mean took 28.7684 seconds.

From the results in Table 2 we can see that the use of the
$\widehat{\mathbf{L}}$ proposed in section 3 is very effective in reducing the
overall iteration count for all $\widehat{\mathbf{S}}$,
$\mathcal{P}_{\mathcal{D}}$, $\mathcal{P}_{\mathcal{T}}$, and
$\mathcal{P}_{\mathcal{C}}$. The number of iterations achieved by
$\widehat{\mathbf{S}}$, $\mathcal{P}_{\mathcal{D}}$ and
$\mathcal{P}_{\mathcal{T}}$ with $r=p$ is remarkably small, especially when
compared to the one attained by employing the $\widehat{\mathbf{L}}$ coming
from [47] with $k=3$. These small numbers of iterations impact on the
wallclock time of the overall solution problem too with
$\widehat{\mathbf{S}}$, $r=p$, with $\widehat{M}=Sym(M_{10})$ being the
fastest approach we tested. However, we note that the approach of [47] is
designed to increase parallelisability of preconditioners, and significant
speed up is to be expected for this preconditioner in a parallel setting (all
experiments presented here are performed in serial).

In Table 3 we report the results obtained for $N=100$ for a selection of
parameters. For our new strategy, we document the results achieved using
$\widehat{M}=Sym(M_{1})$. We can see that $\widehat{\mathbf{S}}$,
$\mathcal{P}_{\mathcal{D}}$, and $\mathcal{P}_{\mathcal{T}}$ with $r=p$ lead
to a very small number of matCG/matGMRES iterations also for this problem
setting. Moreover, they scale much better than the strategy from [47] in terms
of computational timing in this serial setting. However, we recall from Table
1 that the norm of the difference between the blocks $M_{i}$ is small for
$\Delta t=10^{-6}$. This is the best scenario for our novel preconditioning
approach.

In Table 4 we consider the performance of different approximations
$\widehat{\mathbf{S}}$ for the SPD problem (2). In addition to the low-rank
approach presented in section 3.2 we consider the limited memory
preconditioning (LMP) approach studied in [9]. The LMP approach approximates
eigenvalues of
$\mathbf{I}+\mathbf{D}^{1/2}\mathbf{L}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}\mathbf{L}^{-1}\mathbf{D}^{1/2}$,
i.e. symmetrically preconditioning with the exact
$\mathbf{L}^{T}\mathbf{D}^{-1}\mathbf{L}$ operator. Spectral information is
typically approximated using randomised numerical linear algebra approaches.
In addition to requiring spectral information of a much large linear system
than the approach considered in section 3.2, LMP requires the use of
$\widehat{\mathbf{L}}=\mathbf{L}$, meaning that it cannot be readily applied
in the matrix-oriented approach. Hence, computing the full spectrum of
$\mathbf{I}+\mathbf{D}^{1/2}\mathbf{L}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}\mathbf{L}^{-1}\mathbf{D}^{1/2}$
is prohibitively costly both in terms of storage and computation. We therefore
compare matCG preconditioned with $\widehat{\mathbf{S}}$ as in (18) with
$r=0,r=p$ with vecCG preconditioned in two different ways. In the first place,
we use $\widehat{\mathbf{S}}$ as in (17) where the inverse of the exact
$\mathbf{L}$ is computed by means of the algorithm in [47] by setting $k=N+1$.
The low-rank term $\mathbf{K}_{r}\mathbf{K}_{r}^{T}$ is constructed by
considering $r=10,$ or $r=p$ eigenpairs of $H^{T}RH$. The second
preconditioning approach for vecCG is given by LMP. Also in this case
$\mathbf{L}^{-1}$ is computed by following [47] with $k=N+1$. In LMP, the rank
of the low-rank approximation to
$\mathbf{D}^{1/2}\mathbf{L}^{-T}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}\mathbf{L}^{-1}\mathbf{D}^{1/2}$
is set to $10(N+1)$. Notice that this means that we are employing the same
number of eigenpairs as in (17) for $r=10$. Indeed, if
$V_{r}\Upsilon_{r}V_{r}^{T}\approx H^{T}RH$, $V_{r}\in\mathbb{R}^{p\times
10}$, then $\mathbf{K}_{r}=I_{N+1}\otimes V_{r}\Upsilon_{r}^{1/2}$ has rank
$10(N+1)$.

| $\widehat{\mathbf{S}}$ | $\mathcal{P}_{\mathcal{D}}$ | $\mathcal{P}_{\mathcal{T}}$ | $\mathcal{P}_{\mathcal{C}}$
---|---|---|---|---
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=0$ | 213 | 239 | 173.5 | 242.5
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=p$ | 9 | 15 | 14 | -
vecCG/vecGMRES [47], $k=3$ | 955.4 | - | 938.1 | 570.9
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=0$ | 827.34 | 499.70 | 356.71 | 504.56
matCG/matGMRES, $\widehat{M}=Sym(M_{1})$, $r=p$ | 94.80 | 105.70 | 116.43 | -
vecCG/vecGMRES [47], $k=3$ | 4757 | 5054 | 5553 | 2243
Table 3: Example 5.2. Iterations (top) and wallclock time (bottom) to convergence for the Lorenz96 problem with $N=100$ for different preconditioners. 10 realisations. vecGMRES with $\mathcal{P}_{D}$ and $\mathbf{\widehat{L}}^{-1}$ computed as in [47] ($k=3$) did not converge in 1000 iterations. Preconditioner | Iterations | Wallclock time
---|---|---
matCG, $\widehat{M}=Sym(M_{1})$, $r=0$ | 212.5 | 845.71
matCG, $\widehat{M}=Sym(M_{1})$, $r=p$ | 9 | 98.98
vecCG, $k=N+1$, $r=0$ | 178.5 | 967.73
vecCG, $k=N+1$, $r=10$ | 126.2 | 958.20
vecCG, $k=N+1$, $r=10(N+1)$, LMP | 28 | 119.22
Table 4: Example 5.2. Iterations (left) and wallclock time (right) to
convergence for the Lorenz96 problem with $N=100$ for different
preconditioners which approximate $\widehat{\mathbf{S}}$. Values are averaged
over 10 realisations.

We see that in terms of iterations the LMP approach results in much larger
reductions than the approach of section 3.2 for the same number of eigenpairs.
However, by exploiting the Kronecker structure of the new preconditioning
approach, we can incorporate many more terms, leading to very small number of
iterations for matCG with $r=p$. We also note the improvement in wallclock
times when using the matrix-oriented approach, with the fastest times occuring
for matCG with $r=p$.

Figure 1: Example 5.2. Iterations to reach convergence with changing time
discretization, $\Delta t$, for the Lorenz96 problem for different choices of
preconditioner for $N=10$. Panel (a) shows $\widehat{\mathbf{S}}$, (b) shows
$\mathcal{P}_{\mathcal{D}}$, (c) shows $\mathcal{P}_{\mathcal{T}}$ and (d)
shows $\mathcal{P}_{\mathcal{C}}$. For all panels the red solid line
represents matGMRES/matCG with $r=0$ for $\widehat{M}=Sym(M_{1})$, blue dashed
line represents matGMRES/matCG with $r=p$, and black dot-dashed line
vecGMRES/vecCG for $k=N+1$. For panel (a) the cyan solid line with cross
markers shows vecCG for $k=N+1$ and $r=10(N+1)$ using the LMP approach. We
report averaged behaviour over 10 realisations.

Figure 1 shows how the Kronecker preconditioners perform as $\Delta t$
increases, and the difference between linearised model operators increases. We
compare against the approach of [47] using $k=N+1$, i.e.
$\widehat{\mathbf{L}}\equiv\mathbf{L}$. For the SPD problem, we also plot
convergence for the vectorised approach with LMP. For $\Delta
t\in[10^{-6},10^{-2}]$ both matCG/matGMRES and vecCG/vecGMRES behave
similarly, with only a small increase in iterations with $\Delta t$. However,
for larger values of $\Delta t$ the number of iterations required to reach
convergence increases for all methods, with the largest impact seen for
matCG/matGMRES. The use of $r=p$ within matCG/matGMRES is much more sensitive
to the choice of $\Delta t$, with a steady increase in the number of
iterations required to reach convergence as $\Delta t$ increases. For $\Delta
t=5\times 10^{-2}$ there is negligible benefit in terms of iterations to the
inclusion of observation information within the Schur complement. For large
values of $\Delta t$, the difference between
$\widehat{\mathbf{L}}=I_{N+1}\otimes I_{s}-\Sigma\otimes\widehat{M}$ and
$\mathbf{L}$ increases, meaning that including $r>0$ factors coming from the
observation term is obtaining an improved estimate of the wrong
preconditioner. In the future it might be possible to design alternative
additional terms that can correct for this discrepancy. We note that for LMP
as $\widehat{\mathbf{L}}\equiv\mathbf{L}$ the ‘correct’ low-rank update is
used for all choices of $\Delta t$.

###  5.3 Results for heat equation

We now present an example with $M_{i}=M$ for all $i$. This simpler setting
allows us to validate the theoretical properties of our new approach, and
consider computational aspects such as the cost of the approach proposed in
section 3.2 and scaling of our methods with increasing observation times. Our
second numerical example comes from [47, Section 6.1]. For this example,
$\mathbf{L}$ is a Stein operator of the form (22) where the matrix $M$ amounts
to the discrete operator stemming from the discretization of the one-
dimensional heat equation on the unit line

$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial x^{2}}$

with Dirichlent boundary conditions $u(0,t)=u_{0}$, $u(1,t)=u_{1}$ for all
$t\in(0,1]$. By discretising the equation above by means of the forward Euler
method in time and second-order central differences in space, $M$ can be
written as follows

$M=\begin{bmatrix}0&0&0&0&\cdots&0&0\\\ 0&1-2r&r&0&\cdots&0&0\\\
0&r&\ddots&\ddots&&\vdots&\vdots\\\ 0&0&\ddots&\ddots&\ddots&0&0\\\
\vdots&&\ddots&\ddots&\ddots&r&0\\\ 0&\cdots&0&0&r&1-2r&0\\\
0&\cdots&0&0&0&0&0\\\ \end{bmatrix},\quad r=\frac{\Delta t}{(\Delta x)^{2}}.$

For this example we use $\Delta x=10^{-3}$ and $\Delta t=4\times 10^{-7}$.

We begin by comparing the performance of the matrix and vector oriented
approaches when using the exact $\widehat{\mathbf{L}}\equiv\mathbf{L}$ in all
four choices of preconditioner proposed in this paper, and $r=0$ in (18). For
vecCG/vecGMRES, $\mathbf{L}^{-1}$ is computed by using the procedure coming
from [47] for $k=N+1$. We remind the reader that the two versions of CG/GMRES
are equivalent in exact arithmetic. However, the matrix oriented approaches
present some computational advantages for this example. Table 5 demonstrates
the remarkable gain in efficiency that occurs for $\widehat{\mathbf{S}}$,
$\mathcal{P}_{\mathcal{D}}$, and $\mathcal{P}_{\mathcal{T}}$ for $N=10$. We
notice that the strategy based on $\mathcal{P}_{\mathcal{C}}$ is faster than
the ones related to the block diagonal and block triangular preconditioners,
in spite of the fact that it requires a larger number of iterations. This is
due to the small $N$ selected in this example which makes the inversion of
$\mathbf{D}$ the most expensive step in the preconditioners that involve
$\mathbf{\widehat{S}}$. We recall that $\mathcal{P}_{\mathcal{C}}$ avoids the
application of $\mathbf{D}^{-1}$. Moreover, for this choice of $N$, vecGMRES
is faster than matGMRES when $\mathcal{P}_{\mathcal{C}}$ is adopted as
preconditioning operator. This is related to the cost of the
eigendecomposition of $M$ described at the end of section 3.1. This step
cubically depends on $s$; taking about 0.2s for this problem, namely about
$1/2$ of the overall running time achieved by matGMRES. Nevertheless, we
anticipate the cost of this eigendecomposition to be amortised for larger
choices of $N$ (see e.g. Figure 2).

| $\widehat{\mathbf{S}}$ | $\mathcal{P}_{\mathcal{D}}$ | $\mathcal{P}_{\mathcal{T}}$ | $\mathcal{P}_{\mathcal{C}}$
---|---|---|---|---
Iterations – vecCG/vecGMRES | 18.4 | 36.6 | 18.8 | 37.6
Iterations – matCG/matGMRES | 18.4 | 36.6 | 18.8 | 37.6
Wallclock time – vecCG/vecGMRES | 6.2475 | 7.7263 | 4.1889 | 0.4013
Wallclock time – matCG/matGMRES | 1.3356 | 0.8472 | 0.8493 | 0.5713
Table 5: Example 5.3. Iterations and wallclock time to convergence for $N=10$
and $\widehat{\mathbf{L}}\equiv\mathbf{L}$ for the objective function
formulation (2) (column 1) and the saddle point formulation (3) with different
preconditioners (columns 2-4) averaged over 10 realisations.

We now focus on the performance achieved by the novel Schur complement
approximation (17) and its implementation illustrated in section 3.2.1. To
this end, we consider only $\widehat{\mathbf{S}}$ and
$\mathcal{P}_{\mathcal{D}}$. Similar results have also been obtained for
$\mathcal{P}_{\mathcal{T}}$ but are not presented here. In Table 6 (left), for
different values of $N$, we report the iteration count and the overall running
time of matCG/matGMRES when the approximate Schur complement
$\widehat{\mathbf{S}}$ is adopted for $\mathbf{S}$ and
$\mathcal{P}_{\mathcal{D}}$ for $r=0$ and $r=p$. We notice that choosing $r=p$
in (17) leads to a remarkable decrease in the number of iterations needed to
converge. Moreover, the number of CG/GMRES iterations we perform turns out to
be $N$-independent for both problem.

We notice that the number of iterations required to solve the matCG
formulation is smaller than for matGMRES, leading to smaller wallclock times.
We note that for $r=p$ we are approximating the inverse of $\mathbf{S}$ to a
small tolerance, and hence obtain convergence in a single iteration of matCG.
Wallclock times for the case $r=p$ are comparable for both problems. Even
though the use of inner CG introduces some inexactness in the preconditioning
step, we notice that the number of iterations performed by matCG/matGMRES with
$\widehat{\mathbf{S}}$ and $r=p$ is independent of $N$ and equal to the number
of iterations expected in case of an exact computation of
$\widehat{\mathbf{S}}^{-1}$.

The large reduction in the iteration count also leads to a significant speed-
up of the overall solution process, which is not obvious in general. Indeed,
the use of (18) for large values of $r$ can be computationally demanding due
to the need to solve the linear system with
$I_{r(N+1)}+\mathbf{K}_{r}^{T}\mathbf{L}^{-1}\mathbf{D}\mathbf{L}^{-T}\mathbf{K}_{r}$.
However, thanks to the matrix CG strategy presented in section 3.2.1, which
takes full advantage of the semi-diagonalization of $\mathbf{L}$, dealing with
$\widehat{\mathbf{S}}$ by (18) turns out to be computationally affordable for
$r=p$. In Table 6 (far right), we report the number of CG iterations and the
related running time needed to approximately invert $\widehat{\mathbf{S}}$
within the matGMRES iteration. We remind the reader that, in light of Theorem
1, $\widehat{\mathbf{S}}^{-1}$ has to be computed only every other GMRES
iteration. From the results in Table 6 (right), we can see that the CG steps
correspond to a small proportion of the overall GMRES running time for small
$N$. However, as $N$ increases, the cost of the inner CG iteration becomes a
larger proportion of the overall wallclock time. The number of CG iterations
increases with $N$, leading to a more demanding preconditioning step for
larger numbers of observation times. In this scenario, equipping the inner CG
solve with effective preconditioning operators may be largely beneficial.
However, as we previously mentioned, natural preconditioning candidates were
not able to reduce the CG iteration count without significantly increasing its
computational cost per iteration. We will explore this challenging topic of
designing bespoke preconditioners for the inner CG solver in the future.

| matCG | matGMRES | Inner CG
---|---|---|---
|  |  |  |  | 2nd GMRES it.
$N$ | Its. | Time | Its. | Time | Its. | Time
10 ($r=0$) | 18.3 | 1.1732 | 36.4 | 1.4561 |  |
($r=p$) | 1 | 0.6683 | 3 | 0.6332 | 20 | 0.0598
20 ($r=0$) | 26.1 | 1.8153 | 51.8 | 2.1692 |  |
($r=p$) | 1 | 0.8008 | 3 | 0.7577 | 30.3 | 0.1420
30 ($r=0$) | 34.2 | 2.3684 | 63.0 | 2.9455 |  |
($r=p$) | 1 | 0.9045 | 3 | 0.9111 | 40 | 0.2674
40 ($r=0$) | 42.1 | 2.9996 | 74.2 | 4.0273 |  |
($r=p$) | 1 | 1.0495 | 3 | 1.0598 | 49.3 | 0.4227
50 ($r=0$) | 49.5 | 3.5876 | 83.0 | 5.1370 |  |
($r=p$) | 1 | 1.2258 | 3 | 1.2289 | 59 | 0.5731
60 ($r=0$) | 57.1 | 4.2990 | 92.6 | 6.7560 |  |
($r=p$) | 1 | 1.5402 | 3 | 1.5267 | 67 | 0.8732
Table 6: Example 5.3. Iterations and wallclock time to convergence for
$\widehat{\mathbf{S}}$ (first two columns) $\mathcal{P}_{\mathcal{D}}$
(columns 3-6) with the two different choices $r=0$ and $r=p$ in (17) as $N$
varies. In columns 5 and 6 we report the iterations and wallclock time needed
by inner CG method to solve the linear system with
$I_{r(N+1)}+\mathbf{K}_{r}^{T}\mathbf{L}^{-1}\mathbf{D}\mathbf{L}^{-T}\mathbf{K}_{r}$
involved in (18) when $r>0$ in $\mathcal{P}_{\mathcal{D}}$.

We conclude the heat equation example by comparing the novel preconditioning
operators developed in this paper with state-of-the-art approaches. In
particular, we consider matCG equipped with $\widehat{\mathbf{S}}$ and
matGMRES equipped with $\mathcal{P}_{\mathcal{D}}$,
$\mathcal{P}_{\mathcal{T}}$ (with both $r=0$ and $r=p$ in (17)), and
$\mathcal{P}_{\mathcal{C}}$. Using $k=N+1$ within the strategy of [47] becomes
infeasible for large values of $N$. We therefore use the strategy coming from
[47] with $k=3$ to approximate $\widehat{\mathbf{L}}^{-1}$ as a comparison.
This is implemented with vecCG/vecGMRES and is otherwise equipped with the
same preconditioning frameworks as our novel approach.

Figure 2: Example 5.3. Iterations (top) and wallclock time (bottom) to reach
convergence with increasing problem size (number of observation times) for the
heat equation problem for different choices of preconditioner. Red solid lines
denote matCG/matGMRES with $r=0$, blue dashed lines denote matCG/matGMRES with
$r=p$, and black dot-dashed lines denote vecCG/vecGMRES with $k=3$. Crosses
denote $\widehat{\mathbf{S}}$, pluses denote $\mathcal{P}_{\mathcal{D}}$,
circles denote $\mathcal{P}_{\mathcal{T}}$ and triangle denote
$\mathcal{P}_{\mathcal{C}}$. We report averaged behaviour over 10
realisations.

Figure 2 (top) shows the number of iterations to reach convergence with an
increasing number of observation times $N$. We can observe that
$\widehat{\mathbf{S}}$, $\mathcal{P}_{\mathcal{D}}$ and
$\mathcal{P}_{\mathcal{T}}$ require a very small number of iterations whenever
we select $r=p$ in (17). The performance of these operators is optimal as the
number of performed CG/GMRES iterations is constant with increasing $N$. We
recall that such optimality is not guaranteed, as the linear system with
$I_{r(N+1)}+\mathbf{K}_{r}^{T}\mathbf{L}^{-1}\mathbf{D}\mathbf{L}^{-T}\mathbf{K}_{r}$
involved in (18) is solved iteratively to a relative residual tolerance of
$10^{-8}$. As previously mentioned, such _inexactness_ does not allow us to
claim that we are working within the scenarios depicted in Corollary 8–10 –
for which we would be able to guarantee an $N$-independent number of CG/GMRES
iterations – even though $r=p$ in (17). Nevertheless, for this example the
overall solution process does demonstrate $N$ independence of iterations to
reach convergence.

Competitive performance is attained also when $r=0$ in the
$\widehat{\mathbf{S}}$, $\mathcal{P}_{\mathcal{D}}$, and
$\mathcal{P}_{\mathcal{T}}$ preconditioning frameworks compared to the
$\mathbf{L}^{-1}$ approximation with $k=3$. The large difference in iterations
for the strategy coming from [47] also leads to much longer wallclock times.

We notice that the use of $\mathcal{P}_{\mathcal{C}}$ leads to a number of
GMRES iterations which is always very similar to the one achieved by
$\mathcal{P}_{\mathcal{D}}$ with $r=0$. The performance of
$\widehat{\mathbf{S}}$ is also similar to the performance of
$\mathcal{P}_{\mathcal{T}}$ in the $r=0$ setting.

In Figure 2 (bottom) we report the computational time of the overall CG/GMRES
solution process for all the preconditioning operators we mentioned above. We
can see that, except for $\mathcal{P}_{\mathcal{C}}$ with small $N$, vecGMRES
is much slower than matGMRES equipped with our novel preconditioning
strategies. In particular, from the results depicted in Figure 2 (bottom) we
see that selecting $r=p$ in $\widehat{\mathbf{S}}$,
$\mathcal{P}_{\mathcal{D}}$ and $\mathcal{P}_{\mathcal{T}}$ is a favourable
choice over $r=0$ also in terms of running time, with competitive scaling with
$N$ when using the improved Schur complement approximation.

##  6 Conclusions and outlook

To fully exploit the rich Kronecker structure of the matrices stemming from
weak-constraint 4D-Var problems, matrix-oriented Krylov methods can be
employed to solve both (2) and (3). The use of such machinery naturally leads
to the design of new preconditioning approaches. In particular, by selecting a
fresh option for the operator $\widehat{\mathbf{L}}$ whose inversion can be
recast in terms of the solution of a Stein matrix equation, we designed
improved preconditioners able to drastically reduce the Krylov iteration count
for certain problems. Our new approach also allows for the efficient inclusion
of information from the observation term of the Schur complement $\mathbf{S}$,
leading to more accurate approximations $\widehat{\mathbf{S}}$.

In the case of observation-time independent forecast models $\mathcal{M}$, our
new preconditioning frameworks achieve optimal performance in terms of the
number of iterations, remarkably without increasing the computational cost of
the overall solution process.

The implementation presented in this paper requires a number of assumptions on
the structure of the data assimilation system, which we hope to relax in
future work. Firstly, the machinery developed here relies on having a moderate
spatial dimension $s$. This assumption is crucial for the Stein equation
solution scheme presented in section 3.1. We plan to extend the
preconditioning framework presented in this paper to the case of sizable $s$
in near future. This can be achieved, e.g., by using projection-based methods
for large-scale Stein equations. The approach currently also requires that a
number of other components of the assimilation problem have a strict Kronecker
structure, meaning that the model error and the observing system are constant
for all observation times. As reported in section 3, we could approximate each
term at the preconditioning level by means of some Kronecker forms. However,
the selection of such approximations may be cumbersome. These aspects will be
investigated elsewhere.

## Acknowledgements

We thank Adam El-Said for his code for the Lorenz 96 weak constraint 4D-Var
assimilation problem. We also thank Ieva Daužickaitė for providing us with her
code for the LMP implementation of [9].

The first author is member of the Italian INdAM Research group GNCS. His work
was partially supported by the research project “Tecniche avanzate per
problemi evolutivi: discretizzazione, algebra lineare numerica,
ottimizzazione” (INdAM - GNCS Project CUP_E55F22000270001).

The second author gratefully acknowledges funding from the Engineering and
Physical Sciences Research Council (EPSRC) grant EP/S027785/1.

## Appendix A

Here we report the proof of Theorem 1.

###### Proof.

We first write

$\mathcal{AP}_{\mathcal{D}}^{-1}=\begin{bmatrix}I&0&\mathbf{L}\widehat{\mathbf{S}}^{-1}\\\
0&I&\mathbf{H}\widehat{\mathbf{S}}^{-1}\\\
\mathbf{L}^{T}\mathbf{D}^{-1}&\mathbf{H}^{T}\mathbf{R}^{-1}&0\\\
\end{bmatrix}.$

We show the statement by induction on $k\geq 1$.

For $k=1$, we define

$v_{1}=\frac{1}{\sqrt{\|b\|^{2}+\|d\|^{2}}}\begin{bmatrix}b\\\ d\\\
0\end{bmatrix}.$

Then,

$\widetilde{v}_{2}=\mathcal{AP}_{\mathcal{D}}^{-1}v_{1}=\frac{1}{\sqrt{\|b\|^{2}+\|d\|^{2}}}\begin{bmatrix}b\\\
d\\\
\mathbf{L}^{T}\mathbf{D}^{-1}b+\mathbf{H}^{T}\mathbf{R}^{-1}d\end{bmatrix},$

and the latter vector needs to be orthogonalized with respect to $v_{1}$. A
direct computation shows that the outcome of this orthogonalization is

$\widehat{v}_{2}=\frac{1}{\sqrt{\|b\|^{2}+\|d\|^{2}}}\begin{bmatrix}0\\\ 0\\\
\mathbf{L}^{T}\mathbf{D}^{-1}b+\mathbf{H}^{T}\mathbf{R}^{-1}d\end{bmatrix},$

and $v_{2}=\widehat{v}_{2}/\|\widehat{v}_{2}\|$.

We now assume that the result has been shown for a certain $\bar{k}>1$ and we
prove the inductive step for $\bar{k}+1$.

It holds

$\widetilde{v}_{2(\bar{k}+1)-1}=\widetilde{v}_{2\bar{k}+1}=\mathcal{AP}_{\mathcal{D}}^{-1}v_{2\bar{k}}=\begin{bmatrix}\mathbf{L}\widehat{\mathbf{S}}^{-1}z_{2\bar{k}}\\\
\mathbf{H}\widehat{\mathbf{S}}^{-1}z_{2\bar{k}}\\\ 0\end{bmatrix}.$

Then the orthogonalization step is such that

$\widehat{v}_{2(\bar{k}+1)-1}=\widetilde{v}_{2(\bar{k}+1)-1}-\sum_{j=1}^{\bar{k}}\alpha_{j}v_{2j-1}-\sum_{j=1}^{\bar{k}}\beta_{j}v_{2j},$

and the only term that may potentially contribute to the third block of
$\widehat{v}_{2(\bar{k}+1)-1}$ is

$\sum_{j=1}^{\bar{k}}\beta_{j}v_{2j}=\begin{bmatrix}0\\\ 0\\\
\sum_{j=1}^{\bar{k}}\beta_{j}z_{2j}\end{bmatrix}.$

However, all the scalars $\beta_{j}$’s are zero since

$\beta_{j}=\widetilde{v}_{2(\bar{k}+1)-1}^{T}v_{2j}=\begin{bmatrix}(\mathbf{L}\widehat{\mathbf{S}}^{-1}z_{2\bar{k}})^{T},(\mathbf{H}\widehat{\mathbf{S}}^{-1}z_{2\bar{k}})^{T},0\end{bmatrix}\begin{bmatrix}0\\\
0\\\ z_{2j}\end{bmatrix}=0.$

Therefore,
$v_{2(\bar{k}+1)-1}=\widehat{v}_{2(\bar{k}+1)-1}/\|\widehat{v}_{2(\bar{k}+1)-1}\|$
has a third zero block.

To conclude, if
$v_{2(\bar{k}+1)-1}=[u_{2(\bar{k}+1)-1}^{T},w_{2(\bar{k}+1)-1}^{T},0]^{T}$,
then

$\widetilde{v}_{2(\bar{k}+1)}=\mathcal{AP}_{\mathcal{D}}^{-1}v_{2(\bar{k}+1)-1}=\begin{bmatrix}u_{2(\bar{k}+1)-1}\\\
w_{2(\bar{k}+1)-1}\\\
\mathbf{L}^{T}\mathbf{D}^{-1}u_{2(\bar{k}+1)-1}+\mathbf{H}^{T}\mathbf{R}^{-1}w_{2(\bar{k}+1)-1}\end{bmatrix},$

and orthonormalizing such a vector with respect to the computed basis, and in
particular $v_{2(\bar{k}+1)-1}$, leads to a $v_{2(\bar{k}+1)}$ whose third
block is the only nonzero block. Hence the result of Theorem 2.1 holds by
induction. ∎

Here we report the proof of Proposition 2.

###### Proof.

We begin by observing that
$\widehat{\mathbf{L}}^{-T}\mathbf{L}^{T}\mathbf{L}\widehat{\mathbf{L}}^{-1}=I+A(\widehat{M})$
where the $(i,j)$th block of $A(\widehat{M})$ is given by

$\begin{cases}\sum_{k=i}^{N}\widehat{M}^{(k-i)T}D_{k}^{T}D_{k}\widehat{M}^{k-i}&\emph{
if }i=j,\\\
D_{i-1}\widehat{M}^{i-j-1}+\sum_{k=i}^{N}\widehat{M}^{(k-i)T}D_{k}^{T}D_{k}\widehat{M}^{k-j}&\emph{
if }i>j,\\\
\widehat{M}^{(j-i-1)T}D_{j-1}^{T}+\sum_{k=j}^{N}\widehat{M}^{(k-i)T}D_{k}^{T}D_{k}\widehat{M}^{k-j}&\emph{
if }j>i.\\\ \end{cases}$ (29)

We then write $A(\widehat{M})=\sum_{m=1}^{N}A_{m}$ where $A_{m}$ contains
terms which depend only on $D_{m}$ and powers of $\widehat{M}$ and their
transposes only. We bound the eigenvalues of $A(\widehat{M})$ above by
applying [4, Equation 5.12.2] to obtain

$\lambda_{\max}(A(\widehat{M}))\leq 1+\sum_{k=1}^{N}\lambda_{\max}({A_{k}}).$
(30)

We now bound the eigenvalues of $A_{m}$. For $m=1,\dots,N$, the $(i,j)$th
block of $A_{m}$ is given by

$\begin{cases}\widehat{M}^{(m-i)T}D_{m}^{T}D_{m}\widehat{M}^{m-j}&\emph{ if
}i,j\leq m,\\\ \widehat{M}^{(m-i)T}D_{m}^{T}&\emph{ if }j=m+1>i,\\\
D_{m}\widehat{M}^{m-i}&\emph{ if }i=m+1>j.\end{cases}$ (31)

For all choices of $m$, $A_{m}$ has $m^{2}-1$ non-zero blocks, and has rank
$2s$. If $0_{\ell}\in\mathbb{R}^{\ell}$ denotes the zero vector of length
$\ell$, the $(m-2)s$ eigenvectors corresponding to zero take the form
$(e_{i},-\widehat{M}e_{i},0_{(m-2)s})$, or
$(0_{s},e_{i},-\widehat{M}e_{i},0_{(m-1)s},\dots,(0_{(m-1)s},e_{i},-\widehat{M}e_{i},0_{s})$.

The non-zero eigenvalues of $A_{m}$ can be found by solving the $s\times s$
system

$\left(D_{m}\left(\sum_{k=0}^{N-1}\widehat{M}^{k}\widehat{M}^{kT}\right)D_{m}^{T}\right)v=\frac{\mu^{2}}{\mu+1}v,$

i.e.

$\mu=0.5(\rho\pm\sqrt{\rho^{2}+4\rho}),$ (32)

where $\rho$ are the eigenvalues of
$(D_{m}(\sum_{k=0}^{m-1}\widehat{M}^{k}\widehat{M}^{kT})D_{m}^{T})$.

By the monotonicity of (32), the largest value of $\mu$ occurs for the largest
value of $\rho$ with the positive option, and the smallest value of $\mu$
occurs for the largest value of $\rho$ taking the negative option. Therefore
an upper bound for $\rho$ provides us with an upper bound for $\mu$, and hence
$\lambda_{\max}(A(\widehat{M}))$.

By similarity

$\displaystyle\max(\rho_{m})$
$\displaystyle=\lambda_{\max}(D_{m}^{T}D_{m}(\sum_{k=0}^{m-1}\widehat{M}^{k}\widehat{M}^{kT}))$
$\displaystyle\leq\lambda_{\max}(D_{m}^{T}D_{m})\lambda_{\max}(\sum_{k=0}^{m-1}\widehat{M}^{k}\widehat{M}^{kT}))$
$\displaystyle\leq\lambda_{\max}(D_{m}^{T}D_{m})\sum_{k=0}^{m-1}\lambda_{\max}(\widehat{M}^{T}\widehat{M})^{k}.$

A loose upper bound can be obtained by defining

$\rho_{N}=\max_{m}\lambda_{\max}(D_{m}^{T}D_{m})\sum_{k=0}^{N-1}\lambda_{\max}(\widehat{M}^{T}\widehat{M})^{k}.$

Moreover,

$\sum_{k=0}^{N-1}\lambda_{\max}(\widehat{M}^{T}\widehat{M})^{k}=\left\\{\begin{array}[]{ll}N,&\text{if}\;\lambda_{\max}(\widehat{M}^{T}\widehat{M})=1,\\\
&\\\
\frac{1-\lambda_{\max}^{N}(\widehat{M}^{T}\widehat{M})}{1-\lambda_{\max}(\widehat{M}^{T}\widehat{M})},&\text{otherwise}.\end{array}\right.$

For every choice of $m$ it holds $\mu_{m}\leq
0.5(\rho_{N}+\sqrt{\rho_{N}^{2}+4\rho_{N}})$, therefore

$\lambda_{max}(A(\widehat{M}))\leq
1+\frac{N}{2}(\rho_{N}+\sqrt{\rho_{N}^{2}+4\rho_{N}}).$

∎

Here we report the proof of Proposition 3.

###### Proof.

We show (14). The proof for (15) is analagous.

By plugging (13) into (11) we get

$Z-\widehat{M}ZC^{T}+\widehat{M}Ze_{N+1}e_{1}^{T}=V.$

Premultiplying by $T^{-1}$ and postmultiplying by $F^{T}$ yields

$\widetilde{Z}-\Lambda\widetilde{Z}\Pi+\Lambda\widetilde{Z}(F^{-T}e_{N+1})(e_{1}^{T}F^{T})=T^{-1}VF^{T},\qquad\widetilde{Z}:=T^{-1}ZF^{T},$

whose Kronecker form is given by

$(I_{N+1}\otimes
I_{s}-\Pi\otimes\Lambda+(Fe_{1})(e_{N+1}^{T}F^{-1})\otimes\Lambda)\text{vec}(\widetilde{Z})=\text{vec}(T^{-1}VF^{T}).$

If $\mathbf{G}:=I_{N+1}\otimes
I_{s}-\Pi\otimes\Lambda\in\mathbb{C}^{(N+1)s\times(N+1)s}$,
$\mathbf{M}:=Fe_{1}\otimes\Lambda\in\mathbb{C}^{(N+1)s\times s}$, and
$\mathbf{N}:=F^{-T}e_{N+1}\otimes I_{s}\in\mathbb{C}^{(N+1)s\times s}$, the
Sherman-Morrison-Woodbury formula [16, Equation (2.1.4)] shows that

$\text{vec}(\widetilde{Z})=\mathbf{G}^{-1}\text{vec}(T^{-1}VF^{T})-\mathbf{G}^{-1}\mathbf{M}(I+\mathbf{N}^{T}\mathbf{G}^{-1}\mathbf{M})^{-1}\mathbf{N}^{T}\mathbf{G}^{-1}\text{vec}(T^{-1}VF^{T}).$

Once $\widetilde{Z}$ is computed, we retrieve $Z$ by performing
$Z=T\widetilde{Z}F^{-T}$.

We now derive a cheap procedure for the computation of $\widetilde{Z}$ which
does not involve the construction of any large matrix.

We first notice that, since $\mathbf{G}$ is diagonal it holds

$\text{vec}(Y)=\mathbf{G}^{-1}\text{vec}(T^{-1}VF^{T})\quad\Longleftrightarrow\quad
Y=P\circ(T^{-1}VF^{T}),$

where $P\in\mathbb{C}^{s\times(N+1)}$ is such that
$P_{i,j}=1/(1-\lambda_{i}\pi_{j})$.

Moreover, by exploiting the Kronecker structure of $\mathbf{N}$, we have

$\mathbf{N}^{T}\mathbf{G}^{-1}\text{vec}(T^{-1}VF^{T})=((e_{N+1}^{T}F^{-1})\otimes
I_{s})\text{vec}(Y)=YF^{-T}e_{N+1}.$

We now focus on the computation of $\mathbf{N}^{T}\mathbf{G}^{-1}\mathbf{M}$.
We remind the reader that

$\mathbf{G}=\begin{pmatrix}I_{s}-\pi_{1}\Lambda&&\\\ &\ddots&\\\
&&I_{s}-\pi_{N+1}\Lambda\\\
\end{pmatrix},\;\mathbf{N}=\begin{pmatrix}(F^{-T}e_{N+1})_{1}I_{s}\\\
\vdots\\\ (F^{-T}e_{N+1})_{N+1}I_{s}\\\
\end{pmatrix},\,\mathbf{M}=\begin{pmatrix}(Fe_{1})_{1}\Lambda\\\ \vdots\\\
(Fe_{1})_{N+1}\Lambda\\\ \end{pmatrix},$

and a direct computation shows that

$\mathbf{N}^{T}\mathbf{G}^{-1}\mathbf{M}=\sum_{j=1}^{N+1}(I_{s}-\pi_{j}\Lambda)^{-1}\Lambda(Fe_{1})_{j}(F^{-T}e_{N+1})_{j}=\text{diag}(P(\Lambda
Fe_{1}\circ F^{-T}e_{N+1})).$

The formulation above provides a cheap expression for the construction of
$\mathbf{N}^{T}\mathbf{G}^{-1}\mathbf{M}$ along with illustrating its diagonal
structure, hence solving the linear system with
${\color[rgb]{0,0,0}U}=I+\text{diag}(P(\Lambda Fe_{1}\circ F^{-T}e_{N+1}))$
does not significantly increase the cost of computing $Z$.

Returning to the computation of $\widetilde{Z}$, we have

$\displaystyle\mathbf{G}^{-1}\mathbf{M}(I+\mathbf{N}^{T}\mathbf{G}^{-1}\mathbf{M})^{-1}\mathbf{N}^{T}\mathbf{G}^{-1}\text{vec}(T^{-1}VF^{T})=$
$\displaystyle\mathbf{G}^{-1}\mathbf{M}{\color[rgb]{0,0,0}U}^{-1}YF^{-T}e_{N+1}$
$\displaystyle=$ $\displaystyle
P\circ(\Lambda{\color[rgb]{0,0,0}U}^{-1}YF^{-T}e_{N+1}e_{1}^{T}F^{T})$
$\displaystyle=$ $\displaystyle W.$

Combining the steps above yields the statement in (14). ∎

## Appendix B

For the sake of completeness, in Algorithm 4 we report the pseudocode of the
matrix-oriented GMRES method applied to (3). The $m$-th basis vector of the
Krylov subspace $K_{m}(\mathcal{A},\mathbf{b})$ is represented in terms of the
matrices $\mathcal{V}_{1,m}\in\mathbb{R}^{s\times(N+1)}$,
$\mathcal{V}_{2,m}\in\mathbb{R}^{p\times(N+1)}$, and
$\mathcal{V}_{3,m}\in\mathbb{R}^{s\times(N+1)}$, namely

$v_{m}=\text{vec}\begin{bmatrix}\mathcal{V}_{1,m}\\\ \mathcal{V}_{2,m}\\\
\mathcal{V}_{3,m}\end{bmatrix}.$

The residual norm in line 4 can be cheaply computed by following, e.g., the
classic Given rotation approach presented in [39, Section 6.5.3].

Similarly, in Algorithm 5 we report the pseudocode of the matrix-oriented CG
method applied to (2).

In what follows, $(A)_{i,j}$ will denote the $(i,j)$-th entry of the matrix
$A$.

input :

$B,Q_{1},\ldots,Q_{N}\in\mathbb{R}^{s\times s},$
$R_{0},\ldots,R_{N}\in\mathbb{R}^{p\times p}$,
$H_{0},\ldots,H_{N}\in\mathbb{R}^{p\times s}$,
$M_{1},\ldots,M_{N}\in\mathbb{R}^{s\times s}$
$[b_{0},c_{1},\ldots,c_{N}]\in\mathbb{R}^{s\times(N+1)}$,
$[d_{0},\ldots,d_{N}]\in\mathbb{R}^{p\times(N+1)}$, $m_{\max}$,
$\varepsilon>0$

output :

$\delta\Theta_{m},\delta X_{m}\in\mathbb{R}^{s\times(N+1)}$,
$\delta\Lambda_{m}^{T}\in\mathbb{R}^{p\times(N+1)}$ such that
$\text{vec}(\delta\Theta_{m}^{T},\delta\Lambda_{m}^{T},\delta X_{m}^{T})^{T}$
is an approximate solution to (3).

Compute
$\beta=\sqrt{\|[b_{0},c_{1},\ldots,c_{N}]\|_{F}^{2}+\|[d_{0},\ldots,d_{N}]\|_{F}^{2}}$
and set $\mathcal{V}_{1,1}=[b_{0},c_{1},\ldots,c_{N}]/\beta$,
$\mathcal{V}_{2,1}=[d_{0},\ldots,d_{N}]/\beta$, and $\mathcal{V}_{3,1}=0$

for  _$m=1,2,\dots,$ till $m_{\max}$_ do

2 Set
$\widehat{\mathcal{V}}_{1,m+1}=[B\mathcal{V}_{1,m}e_{1},Q_{1}\mathcal{V}_{1,m}e_{2},\ldots,Q_{N}\mathcal{V}_{1,m}e_{N+1}]+[\mathcal{V}_{3,m}e_{1},\mathcal{V}_{3,m}e_{2}-M_{1}\mathcal{V}_{3,m}e_{1},\ldots,\mathcal{V}_{3,m}e_{N+1}-M_{N}\mathcal{V}_{3,m}e_{N}]$
$\widehat{\mathcal{V}}_{2,m+1}=[R_{0}\mathcal{V}_{2,m}e_{1},\ldots,R_{N}\mathcal{V}_{2,m}e_{N+1}]+[H_{0}\mathcal{V}_{3,m}e_{1},\ldots,H_{N}\mathcal{V}_{3,m}e_{N+1}]$
$\widehat{\mathcal{V}}_{3,m+1}=[\mathcal{V}_{1,m}e_{1}-M_{1}^{T}\mathcal{V}_{1,m}e_{2},\ldots,\mathcal{V}_{1,m}e_{N}-M_{N}^{T}\mathcal{V}_{1,m}e_{N+1},\mathcal{V}_{1,m}e_{N}]+[H_{0}^{T}\mathcal{V}_{2,m}e_{1},\ldots,H_{N}^{T}\mathcal{V}_{2,m}e_{N+1}]$
Set $(\underline{T}_{m})_{j,m}=0$ for $j=1,\ldots,m+1$ for  _$\ell=1,2$_ do

3 Compute
$(\underline{T}_{m})_{j,m}=(\underline{T}_{m})_{j,m}+\sqrt{\text{trace}(\widehat{\mathcal{V}}_{1,m+1}^{T}\mathcal{V}_{1,j})^{2}+\text{trace}(\widehat{\mathcal{V}}_{2,m+1}^{T}\mathcal{V}_{2,j})^{2}+\text{trace}(\widehat{\mathcal{V}}_{3,m+1}^{T}\mathcal{V}_{3,j})^{2}},\quad
j=1,\ldots,m$ Set
$\widehat{\mathcal{V}}_{i,m+1}=\widehat{\mathcal{V}}_{1,m+1}-\sum_{j=1}^{m}(\underline{T}_{m})_{j,m}\mathcal{V}_{i,j}$,
for $i=1,2,3$

4 Compute
$(\underline{T}_{m})_{m+1,m}=\sqrt{\|\widehat{\mathcal{V}}_{1,m+1}\|_{F}^{2}+\|\widehat{\mathcal{V}}_{2,m+1}\|_{F}^{2}+\|\widehat{\mathcal{V}}_{3,m+1}\|_{F}^{2}}$
Set
$\mathcal{V}_{i,m+1}=\widehat{\mathcal{V}}_{i,m+1}/(\underline{T}_{m})_{m+1,m}$,
$i=1,2,3$ Solve
$y_{m}=\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\|\underline{T}_{m}y-\beta
e_{1}\|$ Compute the residual norm $\|r_{m}\|$ if
_$\|r_{m}\|\leq\varepsilon\beta$_ then

5 Go to 4

6

Set $\delta\Theta_{m}=\sum_{j=1}^{m}\mathcal{V}_{1,j}(e_{i}^{T}y_{m})$,
$\delta\Lambda_{m}=\sum_{j=1}^{m}\mathcal{V}_{2,j}(e_{i}^{T}y_{m})$, and
$\delta X_{m}=\sum_{j=1}^{m}\mathcal{V}_{3,j}(e_{i}^{T}y_{m})$

Algorithm 4 Matrix-oriented GMRES for (3)

input :

$B,Q_{1},\ldots,Q_{N}\in\mathbb{R}^{s\times s},$
$R_{0},\ldots,R_{N}\in\mathbb{R}^{p\times p}$,
$H_{0},\ldots,H_{N}\in\mathbb{R}^{p\times s}$,
$M_{1},\ldots,M_{N}\in\mathbb{R}^{s\times s}$
$[b_{0},c_{1},\ldots,c_{N}]\in\mathbb{R}^{s\times(N+1)}$,
$[d_{0},\ldots,d_{N}]\in\mathbb{R}^{p\times(N+1)}$, $m_{\max}$,
$\varepsilon>0$

output :

$\delta X_{m}\in\mathbb{R}^{s\times(N+1)}$ approximate solution to (2)

Set
$\mathcal{R}_{0}=W_{0}=[Bb_{0},Q_{1}c_{1},\ldots,Q_{N}c_{N}]+[H_{0}^{T}R_{0}^{-1}d_{0}-M_{1}^{T}H_{1}^{T}R_{1}^{-1}d_{1},\ldots,H_{N-1}^{T}R_{N-1}^{-1}d_{N-1}-M_{N}^{T}H_{N}^{T}R_{N}^{-1}d_{N},H_{N}^{T}R_{N}^{-1}d_{N}]$,
$\delta X_{0}=0$, and compute $\rho_{0}=\|\mathcal{R}_{0}\|_{F}^{2}$

for  _$m=1,2,\dots,$ till $m_{\max}$_ do

2 Set  $\displaystyle\mathcal{W}_{m}=$
$\displaystyle[B^{-1}W_{m-1}e_{1}-M_{1}^{T}Q_{1}^{-1}(W_{m-1}e_{2}-M_{1}W_{m-1}e_{1}),$
$\displaystyle
Q_{1}^{-1}(W_{m-1}e_{2}-M_{1}W_{m-1}e_{1})-M_{2}^{T}Q_{2}^{-1}(W_{m-1}e_{3}-M_{2}W_{m-1}e_{2}),\ldots,Q_{N}^{-1}(W_{m-1}e_{N+1}-M_{N}W_{m-1}e_{N})]$
$\displaystyle+[H_{0}^{T}R_{0}^{-1}H_{0}W_{m-1}e_{1},\ldots,H_{N}^{T}R_{N}^{-1}H_{N}W_{m-1}e_{N+1}]$
$\alpha_{m}=\rho_{m-1}/\text{trace}(\mathcal{W}_{m}^{T}W_{m-1})$ $\delta
X_{m}=\delta X_{m-1}+\alpha_{m}W_{m-1}$
$\mathcal{R}_{m}=\mathcal{R}_{m-1}-\alpha_{k}\mathcal{W}_{m}$
$\rho_{m}=\|\mathcal{R}_{m}\|_{F}^{2}$ if
_$\sqrt{\rho_{m}}\leq\varepsilon\rho_{0}$_ then

3 Return $\delta X_{m}$

4 $\beta_{m}=\rho_{m}/\rho_{m-1}$ $W_{m}=\mathcal{R}_{m}+\beta_{m}W_{m-1}$

Algorithm 5 Matrix-oriented CG for (2)

## References

  * [1] A. Barraud, A numerical algorithm to solve $A^{T}XA-X=Q$, IEEE Transactions on Automatic Control, 22 (1977), pp. 883–885.
  * [2] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems, Acta Numerica, 14 (2005), p. 1–137.
  * [3] M. Benzi and A. J. Wathen, Some preconditioning techniques for saddle point problems, in Model order reduction: theory, research aspects and applications, vol. 13 of Math. Ind., Springer, Berlin, 2008, pp. 195–211.
  * [4] D. S. Bernstein, Matrix mathematics, Princeton university press, 2009\.
  * [5] D. A. Bini and B. Iannazzo, Computing the Karcher mean of symmetric positive definite matrices, Linear Algebra and its Applications, 438 (2013), pp. 1700–1710.
  * [6] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre, Manopt, a Matlab toolbox for optimization on manifolds, Journal of Machine Learning Research, 15 (2014), pp. 1455–1459.
  * [7] A. Carrassi, M. Bocquet, L. Bertino, and G. Evensen, Data assimilation in the geosciences: An overview of methods, issues, and perspectives, Wiley Interdisciplinary Reviews: Climate Change, 9 (2018), p. e535.
  * [8] P. Courtier, J.-N. Thépaut, and A. Hollingsworth, A strategy for operational implementation of 4D-Var, using an incremental approach, Quarterly Journal of the Royal Meteorological Society, 120 (1994), pp. 1367–1387.
  * [9] I. Daužickaitė, A. S. Lawless, J. A. Scott, and P. J. van Leeuwen, Randomised preconditioning for the forcing formulation of weak constraint 4d-var, arXiv preprint arXiv:2101.07249, (2021).
  * [10] A. El-Said, Conditioning of the weak-constraint variational data assimilation problem for numerical weather prediction, PhD thesis, University of Reading, 2015.
  * [11] M. Fisher, S. Gratton, S. Gürol, Y. Trémolet, and X. Vasseur, Low rank updates in preconditioning the saddle point systems arising from data assimilation problems, Optimization Methods and Software, 33 (2018), pp. 45–69.
  * [12] M. Fisher and S. Gürol, Parallelization in the time dimension of four-dimensional variational data assimilation, Quarterly Journal of the Royal Meteorological Society, 143 (2017), pp. 1136–1147.
  * [13] M. Fisher, Y. Tremolet, H. Auvinen, D. Tan, and P. Poli, Weak-constraint and long-window 4D-Var, ECMWF Reading, UK, 2012.
  * [14] M. A. Freitag and D. L. H. Green, A low-rank approach to the solution of weak constraint variational data assimilation problems, J. Comput. Phys., 357 (2018), pp. 263–281.
  * [15] I. Gejadze, H. Oubanas, and V. Shutyaev, Implicit treatment of model error using inflated observation-error covariance, Quarterly Journal of the Royal Meteorological Society, 143 (2017), pp. 2496–2508.
  * [16] G. H. Golub and C. F. Van Loan, Matrix computations, Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, Baltimore, MD, fourth ed., 2013.
  * [17] S. Gratton, S. Gürol, E. Simon, and P. L. Toint, Guaranteeing the convergence of the saddle formulation for weakly constrained 4D-Var data assimilation, Quarterly Journal of the Royal Meteorological Society, 144 (2018), pp. 2592–2602.
  * [18] D. Green, Model order reduction for large-scale data assimilation problems, PhD thesis, University of Bath, 2019.
  * [19] A. Gupta and V. Kumar, Scalability of Parallel Algorithms for Matrix Multiplication, in 1993 International Conference on Parallel Processing - ICPP’93, vol. 3, IEEE, 1993, pp. 115–123.
  * [20] Y. Hao and V. Simoncini, The Sherman–Morrison–Woodbury formula for generalized linear matrix equations and applications, Numerical Linear Algebra with Applications, 28 (2021), p. e2384.
  * [21] M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436.
  * [22] K. Howes, A. M. Fowler, and A. Lawless, Accounting for model error in strong-constraint 4d-var data assimilation, Quarterly Journal of the Royal Meteorological Society, 143 (2017), pp. 1227–1240.
  * [23] J. Huang, J. L. Gómez-Dans, H. Huang, H. Ma, Q. Wu, P. E. Lewis, S. Liang, Z. Chen, J.-H. Xue, Y. Wu, F. Zhao, J. Wang, and X. Xie, Assimilation of remote sensing into crop growth models: Current status and perspectives, Agricultural and Forest Meteorology, 276-277 (2019), p. 107609.
  * [24] B. Iannazzo and M. Porcelli, The Riemannian Barzilai–Borwein method with nonmonotone line search and the matrix geometric mean computation, IMA Journal of Numerical Analysis, 38 (2017), pp. 495–517.
  * [25] K. Jbilou and A. Messaoudi, A computational method for symmetric Stein matrix equations, in Numerical linear algebra in signals, systems and control, vol. 80 of Lect. Notes Electr. Eng., Springer, Dordrecht, 2011, pp. 295–311.
  * [26] M. Kärcher, S. Boyaval, M. A. Grepl, and K. Veroy, Reduced basis approximation and a posteriori error bounds for 4d-var data assimilation, Optim Eng, 19 (2018), p. 663–695.
  * [27] D. Kressner, M. Plešinger, and C. Tobler, A preconditioned low-rank CG method for parameter-dependent Lyapunov matrix equations, Numerical Linear Algebra with Applications, 21 (2014), pp. 666–684.
  * [28] E. N. Lorenz, Predictability: a problem partly solved, in Seminar on Predictability, 4–8 September 1995, vol. 1, Reading, 1995, ECMWF, pp. 1–18.
  * [29] M. F. Murphy, G. H. Golub, and A. J. Wathen, A note on preconditioning for indefinite linear systems, SIAM J. Sci. Comput., 21 (2000), pp. 1969–1972.
  * [30] G. Nakamura and R. Potthast, Inverse modeling, IOP Publishing, 2015\.
  * [31] Y. Notay, Flexible conjugate gradients, SIAM J. Sci. Comput., 22 (2000), pp. 1444–1460.
  * [32] C. C. Paige and M. A. Saunders, Solutions of sparse indefinite systems of linear equations, SIAM J. Numer. Anal., 12 (1975), pp. 617–629.
  * [33] D. Palitta, Matrix equation techniques for certain evolutionary partial differential equations, J Sci Comput, 87 (2021).
  * [34] D. Palitta and P. Kürschner, On the convergence of Krylov methods with low-rank truncations, Numer. Algorithms, 88 (2021), pp. 1383–1417.
  * [35] E. Pinnington, T. Quaife, and E. Black, Impact of remotely sensed soil moisture and precipitation on soil moisture prediction in a data assimilation system with the jules land surface model, Hydrology and Earth System Sciences, 22 (2018), pp. 2575–2588.
  * [36] X. Quan, B. He, X. Li, and Z. Tang, Estimation of grassland live fuel moisture content from ratio of canopy water content and foliage dry biomass, IEEE Geoscience and Remote Sensing Letters, 12 (2015), pp. 1903–1907.
  * [37] T. Rees, Preconditioning iterative methods for PDE constrained optimization, PhD thesis, Citeseer, 2010.
  * [38] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Comput., 14 (1993), pp. 461–469.
  * [39] , Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, second ed., 2003.
  * [40] Y. Saad and M. H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856–869.
  * [41] S. J. Schiff, Neural control engineering: the emerging intersection between control theory and neuroscience, MIT Press, 2011.
  * [42] V. Simoncini, Computational methods for linear matrix equations, SIAM Rev., 58 (2016), pp. 377–441.
  * [43] V. Simoncini and D. B. Szyld, Flexible Inner-Outer Krylov Subspace Methods, SIAM Journal on Numerical Analysis, 40 (2002), pp. 2219–2239.
  * [44] M. Stoll and T. Breiten, A low-rank in time approach to PDE-constrained optimization, SIAM J. Sci. Comput., 37 (2015), pp. B1–B29.
  * [45] J. M. Tabeart, S. L. Dance, S. A. Haben, A. S. Lawless, N. K. Nichols, and J. A. Waller, The conditioning of least-squares problems in variational data assimilation, Numerical Linear Algebra with Applications, 25 (2018), p. e2165.
  * [46] J. M. Tabeart, S. L. Dance, A. S. Lawless, N. K. Nichols, and J. A. Waller, New bounds on the condition number of the hessian of the preconditioned variational data assimilation problem, Numerical Linear Algebra with Applications, 29 (2022), p. e2405.
  * [47] J. M. Tabeart and J. W. Pearson, Saddle point preconditioners for weak-constraint 4d-var, arXiv preprint arXiv:2105.06975, (2021).
  * [48] J. M. Tabeart and J. W. Pearson, Using low-rank observation information to precondition weak-constraint 4D-Var, In prep, (2021).
  * [49] Y. Tr’emolet, Accounting for an imperfect model in 4d-var, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 132 (2006), pp. 2483–2504.
  * [50] J. Tshimanga, S. Gratton, A. T. Weaver, and A. Sartenaer, Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 134 (2008), pp. 751–769.
  * [51] C. F. Van Loan, The ubiquitous kronecker product, Journal of Computational and Applied Mathematics, 123 (2000), pp. 85–100.  Numerical Analysis 2000. Vol. III: Linear Algebra.