# Eigenvalue Bounds for Double Saddle-Point Systems

Susanne Bradley and Chen Greif  Susanne Bradley Department of Computer
Science, The University of British Columbia, Vancouver, Canada V6T 1Z4
<EMAIL_ADDRESS>greif@cs.ubc.ca). Chen Greif22footnotemark: 2

(October 25, 2021)

###### Abstract

We derive bounds on the eigenvalues of a generic form of double saddle-point
matrices. The bounds are expressed in terms of extremal eigenvalues and
singular values of the associated block matrices. Inertia and algebraic
multiplicity of eigenvalues are considered as well. The analysis includes
bounds for preconditioned matrices based on block diagonal preconditioners
using Schur complements, and it is shown that in this case the eigenvalues are
clustered within a few intervals bounded away from zero. Analysis for
approximations of Schur complements is included. Some numerical experiments
validate our analytical findings.

##  1 Introduction

Given positive integer dimensions $n\geq m\geq p$, consider the
$(n+m+p)\times(n+m+p)$ double saddle-point system

$\mathcal{K}u=b,$ (1)

where

$\mathcal{K}=\begin{bmatrix}A&B^{T}&0\\\ B&-D&C^{T}\\\
0&C&E\end{bmatrix};\quad u=\begin{bmatrix}x\\\ y\\\ z\end{bmatrix};\quad
b=\begin{bmatrix}p\\\ q\\\ r\end{bmatrix}.$ (2)

In (2), $A\in\mathbb{R}^{n\times n}$ is assumed symmetric positive definite,
$D\in\mathbb{R}^{m\times m}$ and $E\in\mathbb{R}^{p\times p}$ are positive
semidefinite, and $B\in\mathbb{R}^{m\times n}$, $C\in\mathbb{R}^{p\times m}.$

Under the assumptions stated above, $\mathcal{K}$ is symmetric and indefinite,
and solving (1) presents several numerical challenges. When the linear system
is too large for direct solvers to work effectively, iterative methods [35]
are preferred; this is the scenario that is of interest to us in this work. A
minimum residual Krylov subspace solver such as MINRES [23] is a popular
choice due to its optimality and short recurrence properties. Other solvers
may be highly effective as well.

Linear systems of the form (1)–(2) appear frequently in multiphysics problems,
and their numerical solution is of increasing importance and interest. There
is a large number of relevant applications [6, 18, 37]; liquid crystal
problems [2, 28], Darcy-Stokes equations [7, 15], coupled poromechanical
equations [9], magma-mantle dynamics [31], and PDE-constrained optimization
problems [26, 30] are just a small subset of linear systems that fit into this
framework.

The matrices $D$ and $E$ are often associated with regularization (or
stabilization). We refer to the special case of $D=E=0$ as an unregularized
form of $\mathcal{K}$, and denote it by $\mathcal{K}_{0}$:

$\mathcal{K}_{0}=\begin{bmatrix}A&B^{T}&0\\\ B&0&C^{T}\\\ 0&C&0\end{bmatrix}.$
(3)

This simpler form is of much potential interest, as it may be considered a
direct generalization of the standard saddle-point form for $2\times 2$ block
matrices:

$\tilde{\mathcal{K}}_{0}=\begin{bmatrix}A&B^{T}\\\ B&0\end{bmatrix}.$ (4)

Matrices of the form (4) have been extensively studied, and their analytical
and numerical properties are well understood; see [5, 32] for excellent
surveys. Some properties of $\mathcal{K}$ follow from appropriately reordering
and partitioning the block matrix and then using known results for
block-$2\times 2$ saddle-point matrices (as given in, for example, [11, 33,
34, 36]). Specifically, $\mathcal{K}$ can be reordered and partitioned into a
$2\times 2$ block matrix

$\tilde{\mathcal{K}}=\left[\begin{array}[]{c c | c}A&0&B^{T}\\\ 0&E&C\\\ \hline\cr B&C^{T}&-D\end{array}\right].$ (5)

While this approach is often effective, we may benefit from considering the
block-$3\times 3$ formulation $\mathcal{K}$ directly, without resorting to
(5). When $E$ is rank-deficient or zero (as occurs, for example, in [10, 19]),
the leading $2\times 2$ block of $\tilde{\mathcal{K}}$ is singular, even if
$A$ is full rank. It is then more challenging to develop preconditioners for
$\tilde{\mathcal{K}}$ or to obtain bounds on its eigenvalues, as we are
restricted to methods that can handle singular leading blocks. Additionally,
by considering the block-$3\times 3$ formulation in deriving eigenvalue
bounds, we will see in this paper that we can derive effective bounds using
the singular values of $B$ and $C$ (along with the eigenvalues of the diagonal
blocks). Analyzing $\tilde{\mathcal{K}}$ using established results for
block-$2\times 2$ matrices (such as those given by Rusten and Winther [34])
instead requires singular values of the larger off-diagonal block
$\begin{bmatrix}B&C^{T}\end{bmatrix}$, which may be more difficult to obtain
than singular values of $B$ and $C$. (We can estimate the singular values of
$\begin{bmatrix}B&C^{T}\end{bmatrix}$ in terms of those of $B$ and $C$, but if
these estimates are loose it will result in loose eigenvalue bounds.)

There has been a recent surge of interest in the iterative solution of
multiple saddle-point problems, and our work adds to the increasing body of
literature that considers these problems. Recent papers that provide
interesting analysis are, for example, [2, 3, 6, 25, 37].

The distribution of eigenvalues of $\mathcal{K}$ plays a central role in
determining the efficiency of iterative solvers. It is therefore useful to
gain an understanding of the spectral structure as part of the selection and
employment of solvers. Effective preconditioners are instrumental in
accelerating convergence for sparse and large linear systems; see [4, 24, 35,
38] for general overviews, and [5, 27, 32] for a useful overview of solvers
and preconditioners for saddle-point problems. For multiphysics problems it
has been demonstrated that exploiting the properties of the underlying
discrete differential operators and other characteristics of the problem at
hand is beneficial in the development of robust and fast solvers; see, e.g.,
[9].

There are several potential approaches to the development of preconditioners
for the specific class of problem that we are considering. Monolithic
preconditioners, which work on the entire matrix, have been recently shown to
be extremely effective. Recent work such as [1] has pushed the envelope
towards scalable solvers based on this methodology. Another increasingly
important approach is operator preconditioning based on continuous spaces; see
[14, 20]. This approach relies on the properties of the underlying continuous
differential operators, and uses tools such as Riesz representation and
natural norm considerations to derive block diagonal preconditioners. Block
diagonal preconditioners can also be derived directly by linear algebra
considerations accompanied by properties of discretized PDEs; see, for
example, [8, 27].

The preconditioner we consider for our analysis is:

$\mathcal{M}:=\begin{bmatrix}A&0&0\\\ 0&S_{1}&0\\\ 0&0&S_{2}\end{bmatrix},$
(6)

where

$S_{1}=D+BA^{-1}B^{T};\quad S_{2}=E+CS_{1}^{-1}C^{T}.$ (7)

We assume that $S_{1}$ and $S_{2}$ are both positive definite.

The preconditioner $\mathcal{M}$ is based on Schur complements. It has been
considered in, for example, [6, 37], and is a natural extension of [17, 22]
for block-$2\times 2$ matrices of the form (4). In practice the Schur
complements $S_{1}$ and $S_{2}$ defined in (7) are too expensive to form and
invert exactly. It is therefore useful to consider approximations to those
matrices when a practical preconditioner is to be developed, and we include an
analysis of that scenario.

The recent papers of Sogn and Zulehner [37] and Cai et al. [6] both analyze
the performance of a block-$n\times n$ block diagonal preconditioner analogous
to $\mathcal{M}$ defined in (6), though the former focus on the spectral
properties of the continuous (rather than discretized) preconditioned
operator, and Cai et al. focus their analyses of this preconditioner on the
case where all diagonal blocks except $A$ are zero. Existing analyses of the
eigenvalues of unpreconditioned block-$3\times 3$ matrices have often been
restricted to specific problems, such as interior-point methods in constrained
optimization; see [13, 21].

Our goal in this paper is to provide a general framework for analysis of
eigenvalue bounds for the double saddle-point case, with a rather minimal set
of assumptions on the matrices involved. To our knowledge, ours is the first
work to provide general spectral bounds for the unpreconditioned matrix
$\mathcal{K}$ and for cases where inexact Schur complements are used.

In Section 2 we derive bounds on the eigenvalues of $\mathcal{K}$. In Section
3 we turn our attention to preconditioned matrices of the form
$\mathcal{M}^{-1}\mathcal{K}$. In Section 4 we allow the preconditioners to
contain approximations of Schur complements and an approximate inversion of
the leading block, and show how the bounds are affected as a result. In
Section 5 we validate some of our analytical observations with numerical
experiments. Finally, in Section 6 we draw some conclusions.

Notation. The $n+m+p$ eigenvalues of the unpreconditioned and preconditioned
double saddle-point matrices will be denoted by $\lambda$. We will use $\mu$
with appropriate matrix superscripts and subscripts to denote eigenvalues of
the matrix blocks, and $\sigma$ will accordingly signify singular values. The
eigenvalues of $A$, for example, will be denoted by

$\mu_{i}^{A},\quad i=1,\dots,n,$

and in terms of ordering we will assume that

$\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}>0.$

To increase clarity, we will use $\mu_{\max}^{A}$ to denote $\mu_{1}^{A}$ and
$\mu_{\min}^{A}$ to denote $\mu_{n}^{A}$, and so on.

Based on the above conventions, we use the following notation for eigenvalues
and singular values of the matrices comprising $\mathcal{K}$:

matrix | size | type | number | notation | $\max$ | $\min$
---|---|---|---|---|---|---
$A$ | $n\times n$ | eigenvalues | $n$ | $\mu_{i}^{A},\ i=1,\dots,n$ | $\mu_{\max}^{A}$ | $\mu_{\min}^{A}$
$B$ | $m\times n$ | singular values | $m$ | $\sigma_{i}^{B},\ i=1,\dots,m$ | $\sigma_{\max}^{B}$ | $\sigma_{\min}^{B}$
$C$ | $p\times m$ | singular values | $p$ | $\sigma_{i}^{C},\ i=1,\dots,p$ | $\sigma_{\max}^{C}$ | $\sigma_{\min}^{C}$
$D$ | $m\times m$ | eigenvalues | $m$ | $\mu_{i}^{D},\ i=1,\dots,m$ | $\mu_{\max}^{D}$ | $\mu_{\min}^{D}$
$E$ | $p\times p$ | eigenvalues | $p$ | $\mu_{i}^{D},\ i=1,\dots,p$ | $\mu_{\max}^{E}$ | $\mu_{\min}^{E}$

##  2 Eigenvalue bounds for $\mathcal{K}$

###  2.1 Inertia and solvability conditions

We first discuss the inertia and conditions for nonsingularity of
$\mathcal{K}$ defined in (2). Recall that the inertia of a matrix is the
triplet denoting the number of its positive, negative, and zero eigenvalues
[16, Definition 4.5.6].

######  Proposition 1.

The following conditions are necessary for $\mathcal{K}$ to be invertible:

  1. (i)

$\ker(A)\cap\ker(B)=\emptyset$;

  2. (ii)

$\ker(B^{T})\cap\ker(D)\cap\ker(C)=\emptyset$;

  3. (iii)

$\ker(C^{T})\cap\ker(E)=\emptyset$.

A sufficient condition for $\mathcal{K}$ to be invertible is that $A$,
$S_{1}$, and $S_{2}$ are invertible.

###### Proof.

We begin with the proof of statement (i) by assuming to the contrary that the
intersection of the kernels is not empty – namely, there exists a nonzero
vector $x$ such that $Ax=Bx=0$. This would mean that the block vector
$\begin{bmatrix}x^{T}&0&0\end{bmatrix}^{T}$ was a null vector of
$\mathcal{K}$, which would imply that $\mathcal{K}$ was singular. Similar
reasoning proves (ii) and (iii).

For the sufficient condition, we observe that when $A$, $S_{1}$, and $S_{2}$
are invertible we can write a block-$LDL^{T}$ factorization of $\mathcal{K}$:

$\begin{bmatrix}A&B^{T}&0\\\ B&-D&C^{T}\\\
0&C&E\end{bmatrix}=\begin{bmatrix}I&0&0\\\ BA^{-1}&I&0\\\
0&-CS_{1}^{-1}&I\end{bmatrix}\underbrace{\begin{bmatrix}A&0&0\\\ 0&-S_{1}&0\\\
0&0&S_{2}\end{bmatrix}}_{=:\mathcal{D}}\begin{bmatrix}I&A^{-1}B^{T}&0\\\
0&I&-S_{1}^{-1}C^{T}\\\ 0&0&I\end{bmatrix}.$ (8)

When $A$ and $S_{1}$ are invertible, $S_{2}$ is well-defined and $\mathcal{K}$
is invertible if and only if $\mathcal{D}$ is invertible. The stated result
follows. ∎

Throughout the rest of this paper, we assume that the sufficient condition
holds: namely, that $A$, $S_{1}$, $S_{2}$ are invertible. Given our
assumptions that $D$ and $E$ are semidefinite, this is equivalent to $A$,
$S_{1}$, and $S_{2}$ being positive definite. This allows us to obtain the
following result on the inertia of $\mathcal{K}$, which will be useful in
deriving our bounds.

######  Lemma 1 (Inertia of a double saddle-point matrix).

The matrix $\mathcal{K}$ has $n+p$ positive eigenvalues and $m$ negative
eigenvalues.

###### Proof.

The matrices $A$, $S_{1}$, and $S_{2}$ are symmetric positive definite, by our
assumptions. Sylvester’s Law of Inertia tells us the inertia of $\mathcal{K}$
is the same as that of $\mathcal{D}$ defined in (8); the stated result
follows. ∎

######  Corollary 1.

Let $a,b,c,d,$ and $e$ be scalars with $a>0$, $d,e\geq 0$,
$s_{1}:=d+\frac{b^{2}}{a}>0$ and $s_{2}:=e+\frac{c^{2}}{s_{1}}>0$. Any cubic
polynomial of the form

$p(\lambda)=\lambda^{3}+(d-a-e)\lambda^{2}+(ae-ad-
de-b^{2}-c^{2})\lambda+(ade+ac^{2}+b^{2}e)$

has two positive real roots and one negative real root.

###### Proof.

Consider the $3\times 3$ matrix

$P=\begin{bmatrix}a&b&0\\\ b&-d&c\\\ 0&c&e\end{bmatrix}.$

Using the characteristic polynomial of $P$, it is straightforward to confirm
that

$\displaystyle\det(\lambda I-P)$
$\displaystyle=\lambda^{3}-\operatorname{Tr}(P)\lambda^{2}-\frac{1}{2}\left(\operatorname{Tr}(P^{2})-\operatorname{Tr}^{2}(P)\right)\lambda-\det(P)$
$\displaystyle=p(\lambda).$

Because $P$ is symmetric its eigenvalues are real, and because $P$ is a double
saddle-point matrix with $n=m=p=1$, the two positive and one negative root
follow by Lemma 1. See Figure 1 for a graphical illustration. ∎

Figure 1: Plot of a cubic polynomial $p(\lambda)$ of the form described in
Corollary 1, with two positive roots and one negative root.

###  2.2 Derivation of bounds

Let us define three cubic polynomials, as follows:

$\displaystyle
p(\lambda)=\lambda^{3}+(\mu_{\max}^{D}-\mu_{\min}^{A})\lambda^{2}-\left(\mu_{\min}^{A}\mu_{\max}^{D}+(\sigma_{\max}^{B})^{2}+(\sigma_{\min}^{C})^{2}\right)\lambda+\mu_{\min}^{A}(\sigma_{\min}^{C})^{2};$
(9a)
$\displaystyle\begin{split}q(\lambda)=&\lambda^{3}+(\mu^{D}_{\min}-\mu^{A}_{\max}-\mu^{E}_{\max})\lambda^{2}\\\
&+\left(\mu^{A}_{\max}\mu^{E}_{\max}-\mu^{A}_{\max}\mu^{D}_{\min}-\mu^{D}_{\min}\mu^{E}_{\max}-(\sigma^{B}_{\max})^{2}-(\sigma^{C}_{\max})^{2}\right)\lambda\\\
&+\left(\mu^{A}_{\max}\mu^{D}_{\min}\mu^{E}_{\max}+\mu_{\max}^{A}(\sigma^{C}_{\max})^{2}+(\sigma^{B}_{\max})^{2}\mu^{E}_{\max}\right);\end{split}$
(9b)

$\displaystyle\begin{split}r(\lambda)=&\lambda^{3}+(\mu^{D}_{\max}-\mu^{A}_{\min}-\mu^{E}_{\min})\lambda^{2}\\\
&+\left(\mu^{A}_{\min}\mu^{E}_{\min}-\mu^{A}_{\min}\mu^{D}_{\max}-\mu^{D}_{\max}\mu^{E}_{\min}-(\sigma^{B}_{\max})^{2}-(\sigma^{C}_{\max})^{2}\right)\lambda\\\
&+(\mu^{A}_{\min}\mu^{D}_{\max}\mu^{E}_{\min}+\mu_{\min}^{A}(\sigma^{C}_{\max})^{2}+(\sigma^{B}_{\max})^{2}\mu_{\min}^{E}).\end{split}$
(9c)

All three of these polynomials are of the form described in Corollary 1. Thus,
all roots are real and each polynomial has two positive roots and one negative
root. For notational convenience, we will denote the negative root of
$p(\lambda)$, for example, by $p^{-}$, and use subscripts max and min to
distinguish between the two positive roots. For example, $p^{+}_{\max}$ will
denote the largest positive root and $p^{+}_{\min}$ will denote the smallest
positive root. The same notational rules apply to $q(\lambda)$ and
$r(\lambda)$.

######  Theorem 1 (Eigenvalue bounds, matrix $\mathcal{K}$).

Using the notation established in (2.2), the eigenvalues of $\mathcal{K}$ are
bounded within the intervals

$\left[r^{-},\frac{\mu^{A}_{\max}-\sqrt{(\mu^{A}_{\max})^{2}+4(\sigma^{B}_{\min})^{2}}}{2}\right]\
\bigcup\ \left[p^{+}_{\min},q^{+}_{\max}\right].$ (10)

###### Proof.

Upper bound on positive eigenvalues. We let
$v=\begin{bmatrix}x^{T}&y^{T}&z^{T}\end{bmatrix}^{T}$ be a vector with
$x\in\mathbb{R}^{n}$, $y\in\mathbb{R}^{m}$, and $z\in\mathbb{R}^{p}$. Because
$\mathcal{K}$ is symmetric, we can derive an upper bound on the eigenvalues of
$\mathcal{K}$ by bounding the value of $\frac{v^{T}\mathcal{K}v}{v^{T}v}$. We
can write

$v^{T}\mathcal{K}v=x^{T}Ax-y^{T}Dy+z^{T}Ez+2x^{T}B^{T}y+2y^{T}C^{T}z.$ (11)

We use Cauchy-Schwarz to bound the mixed bilinear forms – such as
$x^{T}B^{T}y$ – by, for example,

$-||x||\cdot||B^{T}y||\leq x^{T}B^{T}y\leq||x||\cdot||B^{T}y||.$

Using this and the eigenvalues/singular values of the block of $\mathcal{K}$,
we can bound (11) from above by

$\displaystyle v^{T}\mathcal{K}v$
$\displaystyle\leq\mu^{A}_{\max}||x||^{2}-\mu^{D}_{\min}||y||^{2}+\mu^{E}_{\max}||z||^{2}+2\sigma^{B}_{\max}||x||\cdot||y||+2\sigma^{C}_{\max}||y||\cdot||z||$
$\displaystyle=\begin{bmatrix}||x||&||y||&||z||\end{bmatrix}\underbrace{\begin{bmatrix}\mu^{A}_{\max}&\sigma^{B}_{\max}&0\\\
\sigma^{B}_{\max}&-\mu^{D}_{\min}&\sigma^{C}_{\max}\\\
0&\sigma^{C}_{\max}&\mu^{E}_{\max}\end{bmatrix}}_{=:R}\begin{bmatrix}||x||\\\
||y||\\\ ||z||\end{bmatrix}.$

An upper bound on $\frac{v^{T}\mathcal{K}v}{v^{T}v}$ is therefore given by the
maximal eigenvalue of $R$. The largest positive eigenvalue of $\mathcal{K}$ is
therefore less than or equal to the largest root of the characteristic
polynomial $\det{(\lambda I-R)}$, which yields the desired result.

Lower bound on negative eigenvalues. The proof is similar to that for the
upper bound on the positive eigenvalues. Using Cauchy-Schwarz and the
eigenvalues/singular values of the blocks, we bound (11) from below by writing

$\displaystyle v^{T}\mathcal{K}v$
$\displaystyle\geq\mu^{A}_{\min}||x||^{2}-2\sigma^{B}_{\max}||x||\cdot||y||-\mu^{D}_{\max}||y||^{2}-2\sigma^{C}_{\max}||y||\cdot||z||+\mu^{E}_{\min}||z||^{2}$
$\displaystyle=\begin{bmatrix}||x||&||y||&||z||\end{bmatrix}\underbrace{\begin{bmatrix}\mu^{A}_{\min}&-\sigma^{B}_{\max}&0\\\
-\sigma^{B}_{\max}&-\mu^{D}_{\max}&-\sigma^{C}_{\max}\\\
0&-\sigma^{C}_{\max}&\mu^{E}_{\min}\end{bmatrix}}_{=:R}\begin{bmatrix}||x||\\\
||y||\\\ ||z||\end{bmatrix}.$

A lower bound on $\frac{v^{T}\mathcal{K}v}{v^{T}v}$ is therefore given by the
smallest eigenvalue of $R$. Taking the characteristic polynomial of $R$ yields
the stated result.

Upper bound on negative eigenvalues. We derive an upper bound on the negative
eigenvalues of $\mathcal{K}$ by finding a lower bound on the negative
eigenvalues of $\mathcal{K}^{-1}$. We begin by partitioning $\mathcal{K}$ as

$\mathcal{K}=\left[\begin{array}[]{c c | c }A&B^{T}&0\\\ B&-D&C^{T}\\\ \hline\cr 0&C&E\end{array}\right]=\begin{bmatrix}\mathcal{K}_{2}&\bar{C}^{T}\\\ \bar{C}&E\end{bmatrix},$ (12)

where $\mathcal{K}_{2}=\begin{bmatrix}A&B^{T}\\\ B&-D\end{bmatrix}$ and
$\bar{C}=\begin{bmatrix}0&C\end{bmatrix}$. By [5, Equation (3.4)],

$\mathcal{K}^{-1}=\begin{bmatrix}\mathcal{K}_{2}^{-1}+\mathcal{K}_{2}^{-1}\bar{C}^{T}S_{2}^{-1}\bar{C}\mathcal{K}_{2}^{-1}&-\mathcal{K}_{2}^{-1}\bar{C}S_{2}^{-1}\\\
-S_{2}^{-1}\bar{C}\mathcal{K}_{2}^{-1}&S_{2}^{-1}\end{bmatrix},$

where $S_{2}=E+\bar{C}\mathcal{K}_{2}^{-1}\bar{C}^{T}=E+CS_{1}^{-1}C^{T}$.
Notice that

$\mathcal{K}^{-1}=\begin{bmatrix}\mathcal{K}_{2}^{-1}&0\\\
0&0\end{bmatrix}+\begin{bmatrix}\mathcal{K}_{2}^{-1}\bar{C}\\\
-I\end{bmatrix}S_{2}^{-1}\begin{bmatrix}\bar{C}^{T}\mathcal{K}_{2}^{-1}&-I\end{bmatrix}.$

Because the second term is positive semidefinite, we conclude that the
eigenvalues of $\mathcal{K}^{-1}$ are greater than or equal to the eigenvalues
of $\begin{bmatrix}\mathcal{K}_{2}^{-1}&0\\\ 0&0\end{bmatrix}$. Thus, a lower
bound on the negative eigenvalues of $\mathcal{K}_{2}^{-1}$ is also a lower
bound on the negative eigenvalues of $\mathcal{K}^{-1}$. This means that an
upper bound on the negative eigenvalues of $\mathcal{K}_{2}$ is also an upper
bound on the negative eigenvalues of $\mathcal{K}$. The desired result now
follows from Silvester and Wathen [36, Lemma 2.2].

Lower bound on positive eigenvalues. We begin by noting that, because $E$ is
positive semidefinite, the eigenvalues of $\mathcal{K}$ are greater than or
equal to those of

$\mathcal{K}_{E=0}=\begin{bmatrix}A&B^{T}&0\\\ B&-D&C^{T}\\\
0&C&0\end{bmatrix}.$

Moreover, $\mathcal{K}$ and $\mathcal{K}_{E=0}$ have the same inertia, by
Theorem 2. Thus, the smallest positive eigenvalue of $\mathcal{K}$ is greater
than or equal to the smallest positive eigenvalue of $\mathcal{K}_{E=0}$. We
therefore obtain a lower bound on the positive eigenvalues of $\mathcal{K}$ by
using energy estimates with $\mathcal{K}_{E=0}$. The eigenvalue problem
associated with $\mathcal{K}_{E=0}$ is

$\begin{bmatrix}A&B^{T}&0\\\ B&-D&C^{T}\\\
0&C&0\end{bmatrix}\begin{bmatrix}x\\\ y\\\
z\end{bmatrix}=\lambda\begin{bmatrix}x\\\ y\\\ z\end{bmatrix}.$ (13)

We pre-multiply the first block row by $x^{T}$, the second by $y^{T}$, and the
third by $z^{T}$ to obtain

$\displaystyle x^{T}Ax+x^{T}B^{T}y$ $\displaystyle=\lambda x^{T}x;$ (14a)
$\displaystyle y^{T}Bx-y^{T}Dy+y^{T}C^{T}z$ $\displaystyle=\lambda y^{T}y;$
(14b) $\displaystyle z^{T}Cy$ $\displaystyle=\lambda z^{T}z.$ (14c)

The third block row of (13) gives $z=\frac{1}{\lambda}Cy$. The first block row
gives

$B^{T}y=(\lambda I-A)x.$

We now consider two cases based on the value of $\lambda$.

#### Case I: $\lambda<\mu_{\min}^{A}$:

If $\lambda<\mu_{\min}^{A}$, then $(\lambda I-A)$ is negative definite, so we
can write

$x=(\lambda I-A)^{-1}B^{T}y.$

Substituting these into (14b) and rearranging gives

$\displaystyle\lambda y^{T}y$ $\displaystyle=y^{T}B(\lambda
I-A)^{-1}B^{T}y-y^{T}Dy+\frac{1}{\lambda}y^{T}C^{T}Cy$ (15a)
$\displaystyle\geq\frac{(\sigma_{\max}^{B})^{2}}{\lambda-\mu_{\min}^{A}}y^{T}y-\mu_{\max}^{D}y^{T}y+\frac{(\sigma_{\min}^{C})^{2}}{\lambda}y^{T}y.$
(15b)

Dividing by $y^{T}y$, using the fact that $\lambda>0$ and
$\lambda-\mu_{\min}^{A}<0$, and rearranging yields

$\underbrace{\lambda^{3}+(\mu_{\max}^{D}-\mu_{\min}^{A})\lambda^{2}-\left(\mu_{\min}^{A}\mu_{\max}^{D}+(\sigma_{\max}^{B})^{2}+(\sigma_{\min}^{C})^{2}\right)\lambda+\mu_{\min}^{A}(\sigma_{\min}^{C})^{2}}_{=p(\lambda)}\leq
0.$

By applying Corollary 1 with $a=\mu_{\min}^{A}$, $b=\sigma_{\max}^{B}$,
$c=\sigma_{\min}^{C}$, and $d=\mu_{\max}^{D}$, we know that this polynomial
has two positive roots. Moreover, $p(\lambda)$ is negative between these two
roots: this follows from the fact that $p(0)>0$ and
$\lim_{\lambda\rightarrow\infty}p(\lambda)=\infty$. Therefore, we conclude
that in this case $\lambda$ is greater than or equal to the smaller positive
root of $p(\lambda)$, namely $p^{+}_{\min}$.

#### Case II: $\lambda\geq\mu_{\min}^{A}$:

If $\lambda\geq\mu_{\min}^{A}$, then $(\lambda I-A)$ may be indefinite (and
possibly singular). However, we can still obtain a lower bound for $\lambda$
by observing that $\mu_{\min}^{A}$ is greater than $p^{+}_{\min}$. As stated
earlier, $p(\lambda)$ has two positive roots and the value of $p(\lambda)$ is
negative between those roots. Moreover, these are the only positive values of
$\lambda$ for which $p(\lambda)<0$. We observe, after simplification, that

$p(\mu_{\min}^{A})=-(\sigma_{\max}^{B})^{2}\mu_{\min}^{A}<0.$

Therefore, the bound $\lambda\geq p^{+}_{\min}$ also holds in this case, which
completes the proof. ∎

#### Remark 1.

From Theorem 1 we see that when $B$ and $C$ are rank deficient, the internal
bounds are zero. Under mild conditions on the ranks and kernels of $D$ and
$E$, the statement of the theorem and its proof may be revised to obtain
nonzero internal bounds in this case. However, doing so is rather technical,
and since the case of full rank $B$ and $C$ is common, further details on this
end case are omitted.

The matrix $\mathcal{K}_{0}$ is a special case of $\mathcal{K}$, and bounds on
its eigenvalues can be obtained as a direct consequence of Theorem 1.

######  Corollary 2 (Eigenvalue bounds, matrix $\mathcal{K}_{0}$).

Define the following three cubic polynomials as special cases of $p,q$ and $r$
defined in (2.2) with $D=E=0$:

$\displaystyle\hat{p}(\lambda)$
$\displaystyle=\lambda^{3}-\mu_{\min}^{A}\lambda^{2}-\left((\sigma_{\max}^{B})^{2}+(\sigma_{\min}^{C})^{2}\right)\lambda+\mu_{\min}^{A}(\sigma_{\min}^{C})^{2};$
(16a) $\displaystyle\hat{q}(\lambda)$
$\displaystyle=\lambda^{3}-\mu^{A}_{\max}\lambda^{2}-\left((\sigma^{B}_{\max})^{2}+(\sigma^{C}_{\max})^{2}\right)\lambda+\mu_{A}^{\max}(\sigma^{C}_{\max})^{2};$
(16b) $\displaystyle\hat{r}(\lambda)$
$\displaystyle=\lambda^{3}-\mu^{A}_{\min}\lambda^{2}-\left((\sigma^{B}_{\max})^{2}+(\sigma^{C}_{\max})^{2}\right)\lambda+\mu^{A}_{\min}(\sigma^{C}_{\max})^{2}.$
(16c)

Using (16), the eigenvalues of $\mathcal{K}_{0}$ are bounded within the
intervals

$\left[\hat{r}^{-},\frac{\mu^{A}_{\max}-\sqrt{(\mu^{A}_{\max})^{2}+4(\sigma^{B}_{\min})^{2}}}{2}\right]\
\bigcup\ \left[\hat{p}^{+}_{\min},\hat{q}^{+}_{\max}\right].$ (17)

The proof of Corollary 2 is omitted; it is similar to the proof of Theorem 1,
with simplifications arising from setting $D=E=0.$

###  2.3 Tightness of the bounds

While the tightness of the bounds we have obtained depends on the problem at
hand, all bounds presented in this section are attainable, as we will
demonstrate with small examples. For the extremal bounds (upper positive and
lower negative), we note that both bounds hold for all double saddle-point
matrices with $n=m=p=1$, as the characteristic polynomial of the $3\times 3$
matrix is the same as the polynomial given in the upper positive and lower
negative bounds of Theorem 1.

For the upper bound on negative eigenvalues, we consider the example (with
$n=m=2$, and $p=1$):

$\mathcal{K}=\left[\begin{array}[]{c c | c c | c }\mu_{\max}^{A}&0&\sigma_{\min}^{B}&0&0\\\ 0&\mu_{\max}^{A}&0&\sigma_{\min}^{B}&0\\\ \hline\cr\sigma_{\min}^{B}&0&0&0&0\\\ 0&\sigma_{\min}^{B}&0&-\mu^{D}&\sigma^{C}\\\ \hline\cr 0&0&0&\sigma^{C}&\mu^{E}\end{array}.\right]$

We can permute the rows and columns of $\mathcal{K}$ to obtain a block
diagonal matrix:

$\mathcal{K}_{D}=\left[\begin{array}[]{c c | c c c }\mu_{\max}^{A}&\sigma_{\min}^{B}&0&0&0\\\ \sigma_{\min}^{B}&0&0&0&0\\\ \hline\cr 0&0&\mu_{\max}^{A}&\sigma_{\min}^{B}&0\\\ 0&0&\sigma_{\min}^{B}&-\mu^{D}&\sigma^{C}\\\ 0&0&0&\sigma^{C}&\mu^{E}\end{array}\right].$

The upper left block of $\mathcal{K}_{D}$ has as an eigenvalue
$\frac{\mu^{A}_{\max}-\sqrt{(\mu^{A}_{\max})^{2}+4(\sigma^{B}_{\min})^{2}}}{2}$,
which is the bound given by Theorem 1.

Finally, for the lower positive bound, we consider the matrix (with $n=m=p=2$)

$\mathcal{K}=\left[\begin{array}[]{ c c | c c | c c}\mu_{\min}^{A}&0&\sigma_{\max}^{B}&0&0&0\\\ 0&\mu_{\min}^{A}&0&\sigma_{\max}^{B}&0&0\\\ \hline\cr\sigma_{\max}^{B}&0&-\mu_{\max}^{D}&0&\sigma_{\min}^{C}&0\\\ 0&\sigma_{\max}^{B}&0&-\mu_{\max}^{D}&0&\sigma_{\min}^{C}\\\ \hline\cr 0&0&\sigma_{\min}^{C}&0&0&0\\\ 0&0&0&\sigma_{\min}^{C}&0&\mu^{E}\end{array}\right].$

As in the previous example, we can permute the rows and columns to obtain a
block diagonal matrix:

$\mathcal{K}_{D}=\left[\begin{array}[]{c c c | c c c}\mu_{\min}^{A}&\sigma_{\max}^{B}&0&0&0&0\\\ \sigma_{\max}^{B}&-\mu_{\max}^{D}&\sigma_{\min}^{C}&0&0&0\\\ 0&\sigma_{\min}^{C}&0&0&0&0\\\ \hline\cr 0&0&0&\mu_{\min}^{A}&\sigma_{\max}^{B}&0\\\ 0&0&0&\sigma_{\max}^{B}&-\mu_{\max}^{D}&\sigma_{\min}^{C}\\\ 0&0&0&0&\sigma_{\min}^{C}&\mu^{E}\end{array}\right].$

The characteristic polynomial of the upper left block is precisely
$p(\lambda)$ defined in (9a); thus, the bound of Theorem 1 is obtained.

##  3 Eigenvalue bounds for $\mathcal{M}^{-1}\mathcal{K}$

We now derive bounds for the preconditioned matrix
$\mathcal{M}^{-1}\mathcal{K}$, with $\mathcal{M}$ defined in (6).

###  3.1 Inertia and eigenvalue multiplicity

We begin with some observations on inertia and eigenvalue multiplicity. To
simplify the presentation and proof of this result, we restrict ourselves to
the case that $B$ and $C$ have full row rank. In some later results, we will
lift this restriction to allow for rank-deficient $B$ and $C$.

######  Theorem 2 (Inertia and algebraic multiplicity, matrix
$\mathcal{M}^{-1}\mathcal{K}$).

Let $\mathcal{K}$ be defined as in (2) and $\mathcal{M}$ as in (6), and
suppose that $B$ and $C$ have full row rank. The preconditioned matrix
$\mathcal{M}^{-1}\mathcal{K}$ has:

  1. (i)

$m$ negative eigenvalues;

  2. (ii)

$p$ eigenvalues in $(0,1)$;

  3. (iii)

$n-m$ eigenvalues equal to $1$; and

  4. (iv)

$m$ eigenvalues greater than $1$.

###### Proof.

Because $\mathcal{M}$ is symmetric positive definite, $\mathcal{M}^{1/2}$
exists and is invertible, and the inertias of
$\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}$ and
$\mathcal{M}^{-1}\mathcal{K}$ are equal to the inertia of $\mathcal{K}$. Thus,
Lemma 1 establishes that $\mathcal{M}^{-1}\mathcal{K}$ has $n+p$ positive and
$m$ negative eigenvalues, which proves (i).

For (ii)-(iv) we split the $n+p$ positive eigenvalues into eigenvalues less
than, equal to, or greater than 1. For this we compute the inertia of the
shifted, split-preconditioned matrix

$\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}-I=\begin{bmatrix}0&\tilde{B}^{T}&0\\\
\tilde{B}&-\tilde{D}-I&\tilde{C}^{T}\\\ 0&\tilde{C}&\tilde{E}-I\end{bmatrix},$

where $\tilde{B}=S_{1}^{-1/2}BA^{-1/2}$,
$\tilde{C}=S_{2}^{-1/2}CS_{1}^{-1/2}$, $\tilde{D}=S_{1}^{-1/2}DS_{1}^{-1/2}$,
$\tilde{E}=S_{2}^{-1/2}ES_{2}^{-1/2}$. The number of positive, negative, and
zero eigenvalues of $\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}-I$ will
be equal to the number of eigenvalues of $\mathcal{M}^{-1}\mathcal{K}$ greater
than, less than, or equal to 1, respectively.

Noting that $\tilde{E}+\tilde{C}\tilde{C}^{T}=I$, we write

$\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}-I=\begin{bmatrix}0&\mathcal{B}\\\
\mathcal{B}^{T}&\mathcal{T}\end{bmatrix},$

where

$\mathcal{T}:=\begin{bmatrix}-\tilde{D}-I&\tilde{C}^{T}\\\
\tilde{C}&-\tilde{C}\tilde{C}^{T}\end{bmatrix}\ \ \ \textrm{and}\ \ \
\mathcal{B}=\begin{bmatrix}\tilde{B}^{T}&0\end{bmatrix}.$

The matrix $\mathcal{B}$ is in $\mathbb{R}^{n\times(m+p)}$ and has rank $m$.
We define a matrix

$N:=\begin{bmatrix}\textbf{0}_{m\times p}\\\ I_{p}\end{bmatrix},$

where $\textbf{0}_{m\times p}$ is the $m\times p$ zero matrix and $I_{p}$ is
the $p\times p$ identity matrix. The columns of $N$ form a basis for
$\ker(\mathcal{B})$. Denote the inertia of $M$ by
$\operatorname{In}(M)=(n_{+},n_{-},n_{0})$ . It is well known (see, e.g., [12,
Lemma 3.4]) that

$\operatorname{In}(\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}-I)=\operatorname{In}(N^{T}\mathcal{T}N)+(m,m,n-m).$

Because $N^{T}\mathcal{T}N=-\tilde{C}\tilde{C}^{T}$ is negative definite, this
gives
$\operatorname{In}(\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}-I)=(m,m+p,n-m)$,
which yields (ii)-(iv). ∎

###  3.2 Derivation of bounds

It is possible to obtain eigenvalue bounds on the preconditioned system
$\mathcal{M}^{-1}\mathcal{K}$ by using the results of Theorem 1 on the
(symmetric) preconditioned system
$\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}$; however, some of the
resulting bounds will be loose. The reason for this is that Theorem 1 assumes
no relationships between the blocks of $\mathcal{K}$ and considers each block
individually. In this section, we will derive tight eigenvalue bounds using
energy estimates, which allow us to fully exploit the relationships between
the blocks of the preconditioned system.

We begin by recalling the following result from Horn and Johnson [16, Theorem
7.7.3], which will be useful throughout the analysis that follows.

######  Lemma 2.

Let $M$ and $N$ be symmetric positive semidefinite matrices such that $M+N$ is
positive definite. Then $(M+N)^{-1}M$ is symmetric positive semidefinite with
all eigenvalues in $[0,1]$.

### The case $D=E=0$

When $D$ and $E$ are both zero, $\mathcal{M}^{-1}\mathcal{K}$ has six distinct
eigenvalues given by: $1$, $\frac{1\pm\sqrt{5}}{2}$, and the roots of the
cubic polynomial $\lambda^{3}-\lambda^{2}-2\lambda+1$. We refer to Cai et al.
[6, Theorem 2.2] for proof.

### The case $D=0$ and $E\neq 0$

When $D=0$, it is necessary that $B$ have full row rank in order for $S_{1}$
to be invertible; however, $C$ may be rank-deficient. Suppose that $C^{T}$ has
nullity $k$. The following result holds.

######  Theorem 3 (Eigenvalue bounds, matrix $\mathcal{M}^{-1}\mathcal{K}$,
$D=0$, $E\neq 0$).

When $D=0$ and $E\neq 0$, the eigenvalues of $\mathcal{M}^{-1}\mathcal{K}$ are
given by: $\lambda=1$ with multiplicity $n-m+k$;
$\lambda=\frac{1\pm\sqrt{5}}{2}$, each with multiplicity $m-p+k$; and $p-k$
eigenvalues located in each of the three intervals:
$\left[\zeta^{-},\frac{1-\sqrt{5}}{2}\right)$,
$\left[\zeta_{\min}^{+},1\right)$, and
$\left(\frac{1+\sqrt{5}}{2},\zeta_{\max}^{+}\right]$, where $\zeta^{-}$,
$\zeta_{\min}^{+}$ and $\zeta_{\max}^{+}$ are the roots of the cubic
polynomial $\lambda^{3}-\lambda^{2}-2\lambda+1$, ordered from smallest to
largest. (These are approximately $-1.2470$, $0.4450$, and $1.8019$,
respectively.)

###### Proof.

We write out the (left-)preconditioned operator

$\mathcal{M}^{-1}\mathcal{K}=\begin{bmatrix}I&A^{-1}B^{T}&0\\\
S_{1}^{-1}B&0&S_{1}^{-1}C^{T}\\\ 0&S_{2}^{-1}C&S_{2}^{-1}E\end{bmatrix},$

with corresponding eigenvalue equations

$\displaystyle x+A^{-1}B^{T}y$ $\displaystyle=\lambda x;$ (18a) $\displaystyle
S_{1}^{-1}Bx+S_{1}^{-1}C^{T}z$ $\displaystyle=\lambda y;$ (18b) $\displaystyle
S_{2}^{-1}Cy+S_{2}^{-1}Ez$ $\displaystyle=\lambda z.$ (18c)

We obtain $n-m$ eigenvectors for $\lambda=1$ by choosing $x\in\ker(B)$ and
$y,z=0$. By considering $y\in\ker(C)$ and $z=0$, we obtain eigenvalues
$\lambda=\frac{1\pm\sqrt{5}}{2}$, each with geometric multiplicity $m-p+k$.

For the remaining eigenvalues, we assume that $z\neq 0$ and
$\lambda\notin\\{1,\frac{1\pm\sqrt{5}}{2}\\}$. From (18a) we obtain

$x=\frac{1}{\lambda-1}A^{-1}B^{T}y,$

which we substitute into (18b) and rearrange to get

$y=\frac{\lambda-1}{\lambda^{2}-\lambda-1}S_{1}^{-1}C^{T}z.$

Substituting this into (18c) gives

$\frac{\lambda-1}{\lambda^{2}-\lambda-1}S_{2}^{-1}CS_{1}^{-1}C^{T}z+S_{2}^{-1}Ez=\lambda
z.$ (19)

Because $S_{2}=E+CS_{1}^{-1}C^{T}$, we can write
$S_{2}^{-1}E=I-S_{2}^{-1}CS_{1}^{-1}C^{T}$. We substitute this into (19) and
rearrange to obtain

$\bigg{(}-\lambda^{2}+2\lambda\bigg{)}S_{2}^{-1}CS_{1}^{-1}C^{T}z=\bigg{(}\lambda^{3}-2\lambda^{2}+1\bigg{)}z.$

Let $\\{\delta_{j},v_{j}\\}$, for $1\leq j\leq p$, denote an eigenpair of
$S_{2}^{-1}CS_{1}^{-1}C^{T}$. For $k$ of these eigenpairs corresponding to
$v_{j}\in\ker(C^{T})$, we note that

$S_{2}^{-1}Ev_{j}=(I-S_{2}^{-1}CS_{1}^{-1}C^{T})v_{j}=v_{j},$

and therefore $\begin{bmatrix}0&0&v_{j}^{T}\end{bmatrix}^{T}$ is an
eigenvector of $\mathcal{M}^{-1}\mathcal{K}$ with $\lambda=1$. For the $p-k$
remaining eigenpairs, we have $0<\delta_{j}\leq 1$ (by Lemma 2). We can then
write

$\bigg{(}-\lambda_{j}^{2}+2\lambda_{j}\bigg{)}\delta_{j}z_{j}=\bigg{(}\lambda_{j}^{3}-2\lambda_{j}^{2}+1\bigg{)}z_{j},$

where $\lambda_{j}$ corresponds to an eigenpair $\\{\delta_{j},z_{j}\\}$ of
$S_{2}^{-1}CS_{1}^{-1}C^{T}$. Because $z_{j}\neq 0$ this implies that

$\lambda_{j}^{3}-(2-\delta_{j})\lambda_{j}^{2}-2\delta_{j}\lambda_{j}+1=0.$
(20)

Thus, each of the $p-k$ positive eigenvalues $\delta_{j}$ of
$S_{2}^{-1}CS_{1}^{-1}C^{T}$ yields three distinct corresponding eigenvalues
$\lambda_{j}^{(1)},\lambda_{j}^{(2)}$, and $\lambda_{j}^{(3)}$ of
$\mathcal{M}^{-1}\mathcal{K}$, corresponding to the roots of the cubic
polynomial (20). These $3(p-k)$ eigenvalues, combined with the eigenvalues
described earlier, account for all eigenvalues of
$\mathcal{M}^{-1}\mathcal{K}$. Substituting $\delta_{j}=0$ and $\delta_{j}=1$
into (20) gives us the three intervals for these eigenvalues stated in the
theorem. ∎

### The case $D\neq 0$ and $E=0$

When $D\neq 0$ and $E=0$, the bounds are the same as when $D$ and $E$ are both
nonzero, which are given next.

### The case $D,E\neq 0$

######  Theorem 4 (Eigenvalue bounds, matrix $\mathcal{M}^{-1}\mathcal{K}$,
$D,E\neq 0$).

The eigenvalues of $\mathcal{M}^{-1}\mathcal{K}$ are bounded within the
intervals

$\left[-\frac{1+\sqrt{5}}{2},\frac{1-\sqrt{5}}{2}\right]\cup\left[\zeta_{\min}^{+},\zeta_{\max}^{+}\right],$

where $\zeta_{\min}^{+}$ and $\zeta_{\max}^{+}$ are respectively the smallest
and largest positive roots of the cubic polynomial
$\lambda^{3}-\lambda^{2}-2\lambda+1$ (approximately 0.4450 and 1.8019,
respectively).

###### Proof.

Upper bound on positive eigenvalues. We know from Theorem 2 that the upper
bound on positive eigenvalues is greater than 1, so we assume here that
$\lambda>1$. We begin by writing the eigenvalue equations as

$\displaystyle Ax+B^{T}y$ $\displaystyle=\lambda Ax;$ (21a) $\displaystyle Bx-
Dy+C^{T}z$ $\displaystyle=\lambda S_{1}y;$ (21b) $\displaystyle Cy+Ez$
$\displaystyle=\lambda S_{2}z.$ (21c)

From (21a) we get

$x=\frac{1}{\lambda-1}A^{-1}B^{T}y,$ (22)

and from (21c) we get $z=(\lambda S_{2}-E)^{-1}Cy$ (because $\lambda>1$, we
are guaranteed that $\lambda S_{2}-E=(\lambda-1)E+\lambda CS_{1}^{-1}C^{T}$ is
positive definite). Substituting these values back into (21b) and pre-
multiplying by $y^{T}$ gives

$\left(\frac{1}{\lambda-1}\right)y^{T}BA^{-1}B^{T}y+y^{T}(D+\lambda
S_{1})y-y^{T}C^{T}(\lambda S_{2}-E)^{-1}Cy=0.$

Recalling that $S_{1}=D+BA^{-1}B^{T}$ we can rewrite this as

$\left(\frac{1}{1-\lambda}+\lambda\right)y^{T}S_{1}y+\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy-y^{T}C^{T}(\lambda
S_{2}-E)^{-1}Cy=0.$ (23)

Because $\lambda>1$, we have

$\displaystyle y^{T}C^{T}(\lambda S_{2}-E)^{-1}Cy$ $\displaystyle\leq
y^{T}C^{T}\Big{(}\lambda(S_{2}-E)\Big{)}^{-1}Cy$
$\displaystyle=\frac{1}{\lambda}y^{T}C^{T}\left(CS_{1}^{-1}C^{T}\right)^{-1}Cy.$

Therefore, (23) gives

$\left(\frac{1}{1-\lambda}+\lambda\right)y^{T}S_{1}y+\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy-\frac{1}{\lambda}y^{T}C^{T}(CS_{1}^{-1}C^{T})^{-1}Cy\leq
0.$ (24)

Next, we let $\tilde{y}=S_{1}^{1/2}y$ and rewrite the first and third terms in
(24) in terms of $\tilde{y}$:

$\left(\frac{1}{1-\lambda}+\lambda\right)\tilde{y}^{T}\tilde{y}+\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy-\frac{1}{\lambda}\tilde{y}^{T}\underbrace{S_{1}^{-1/2}C^{T}(CS_{1}^{-1}C^{T})^{-1}CS_{1}^{-1/2}}_{:=P}\tilde{y}\leq
0.$ (25)

Because $P$ defined in (25) is an orthogonal projector, we have
$\tilde{y}^{T}P\tilde{y}\leq\tilde{y}^{T}\tilde{y}$. And because $\lambda>1$,
we have $1-\frac{1}{1-\lambda}>0$, which means that the second term
$\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy$ is non-negative and can be dropped
from the inequality. The inequality (25) therefore becomes

$\left(\frac{1}{1-\lambda}+\lambda\right)\tilde{y}^{T}\tilde{y}-\frac{1}{\lambda}\tilde{y}^{T}{y}\leq
0,$

which we can divide by $\tilde{y}^{T}\tilde{y}$ rearrange to give

$\lambda^{3}-\lambda^{2}-2\lambda+1\leq 0.$

This, combined with the assumption that $\lambda>1$, gives us the stated
result that $\lambda$ is less than or equal to the largest root of
$\lambda^{3}-\lambda^{2}-2\lambda+1$. The polynomial
$\lambda^{3}-\lambda^{2}-2\lambda+1$ is of the form given in Corollary 1 with
$a=b=c=1$, $d=e=0$, so its value is negative between the two positive roots.

Lower bound on negative eigenvalues. We begin from (23) and note that when
$\lambda<0$, by similar reasoning as was shown for the upper bound on positive
eigenvalues,

$\left(\frac{1}{1-\lambda}+\lambda\right)y^{T}S_{1}y+\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy-\frac{1}{\lambda}y^{T}C^{T}(CS_{1}^{-1}C^{T})^{-1}Cy\geq
0.$

Rewriting the inequality in terms of $\tilde{y}=S_{1}^{1/2}y$ gives

$\left(\frac{1}{1-\lambda}+\lambda\right)\tilde{y}^{T}\tilde{y}+\left(1-\frac{1}{1-\lambda}\right)\tilde{y}^{T}S_{1}^{-1/2}DS_{1}^{-1/2}\tilde{y}-\frac{1}{\lambda}\tilde{y}^{T}P\tilde{y}\geq
0,$

where $P$ is the orthogonal projector defined in (25). For the second term,
note that $S_{1}^{-1/2}DS_{1}^{-1/2}$ is similar to $S_{1}^{-1}D$ which has
all eigenvalues between 0 and 1 (by Lemma 2). Thus, both
$\tilde{y}^{T}S_{1}^{-1/2}DS_{1}^{-1/2}\tilde{y}$ and
$\tilde{y}^{T}P\tilde{y}$ are less than or equal to $\tilde{y}^{T}\tilde{y}$
and we can therefore write

$\left(\frac{1}{1-\lambda}+\lambda\right)\tilde{y}^{T}\tilde{y}+\left(1-\frac{1}{1-\lambda}\right)\tilde{y}^{T}\tilde{y}-\frac{1}{\lambda}\tilde{y}^{T}\tilde{y}\geq
0,$

which, after dividing by $\tilde{y}^{T}\tilde{y}$ and simplifying, yields

$\lambda^{2}+\lambda-1\leq 0.$

This along with the assumption that $\lambda<0$ gives the desired bound of
$\lambda\geq-\frac{1+\sqrt{5}}{2}$.

Lower bound on positive eigenvalues. Assume that $0<\lambda<1$ (we know from
Theorem 2 that the lower bound is in this interval). Substituting (22) into
(21b) and solving for $y$ gives

$y=\underbrace{\left(\frac{1}{1-\lambda}BA^{-1}B^{T}+D+\lambda
S_{1}\right)^{-1}}_{=:Q}C^{T}z,$ (26)

When $0<\lambda<1$, the value $\frac{1}{1-\lambda}$ is positive, so we are
guaranteed that $Q$ in (26) is positive definite. If $z\in\ker(C^{T})$, (26)
gives $y=0$ and (22) gives $x=0$, and we can see from the eigenvalue equations
(21a)-(21c) that this eigenvector corresponds to $\lambda=1$ (because
$Ez=S_{2}z$ for $z\in\ker(C^{T})$). This contradicts our assumption that
$\lambda<1$; thus, we assume $z\notin\ker(C^{T})$.

We can then write (21c) as

$C\left(\frac{1}{1-\lambda}BA^{-1}B^{T}+D+\lambda
S_{1}\right)^{-1}C^{T}z+(1-\lambda)Ez-\lambda CS_{1}^{-1}C^{T}z=0.$ (27)

When $0<\lambda<1$, we have $\frac{1}{1-\lambda}>1$. Therefore, if we take the
inner product of $z^{T}$ with (27), replace $D$ by $\frac{1}{1-\lambda}D$, and
drop the non-negative term $(1-\lambda)z^{T}Ez$, we obtain the inequality

$z^{T}C\left(\left(\frac{1}{1-\lambda}+\lambda\right)S_{1}\right)^{-1}C^{T}z-\lambda
z^{T}CS_{1}^{-1}C^{T}z\leq 0.$

After simplifying and dividing by $z^{T}CS_{1}^{-1}C^{T}z$, we obtain

$\lambda^{3}-\lambda^{2}-2\lambda+1\leq 0.$ (28)

This, combined with the assumption that $0<\lambda<1$, gives us that $\lambda$
must be greater than or equal to the smaller positive root of (28), as
required.

Upper bound on negative eigenvalues. Assume that $\lambda<0$. We begin from
Equation (23) and note that $\lambda S_{2}-E$ is negative definite. Therefore,

$\left(\frac{1}{1-\lambda}+\lambda\right)y^{T}S_{1}y+\left(1-\frac{1}{1-\lambda}\right)y^{T}Dy\leq
0.$

Since $1-\frac{1}{1-\lambda}>0$, the $y^{T}Dy$ term is non-negative and can be
dropped while keeping the inequality, giving us

$\left(\frac{1}{1-\lambda}+\lambda\right)y^{T}S_{1}y\leq 0.$

We thus require

$\frac{1}{1-\lambda}+\lambda\leq 0\ \ \textrm{with }\ \lambda<0,$

which leads to the desired bound $\lambda\leq\frac{1-\sqrt{5}}{2}$. ∎

##  4 Bounds for preconditioners with approximations of Schur complements

In practice, it is too expensive to invert $A$, $S_{1}$, and $S_{2}$ exactly.
In this section we examine eigenvalue bounds on matrices of the form
$\tilde{\mathcal{M}}^{-1}\mathcal{K}$, where $\tilde{\mathcal{M}}$ uses
symmetric positive definite and ideally spectrally equivalent approximations
for $A$, $S_{1}$, and $S_{2}$. Specifically, we consider the approximate block
diagonal preconditioner

$\tilde{\mathcal{M}}=\begin{bmatrix}\tilde{A}&0&0\\\ 0&\tilde{S}_{1}&0\\\
0&0&\tilde{S}_{2}\end{bmatrix},$ (29)

with $\tilde{A}$, $\tilde{S}_{1}$ and $\tilde{S}_{2}$ satisfying the
following:

######  Assumption 5.

Let $\Lambda(\cdot)$ denote the spectrum of a matrix. The diagonal blocks of
the approximate preconditioner $\tilde{\mathcal{M}}$, given by $\tilde{A}$,
$\tilde{S}_{1}$, and $\tilde{S}_{2}$, satisfy:

  1. i.

$\Lambda(\tilde{A}^{-1}A)\in[\alpha_{0},\beta_{0}]$;

  2. ii.

$\Lambda(\tilde{S}_{1}^{-1}S_{1})\in[\alpha_{1},\beta_{1}]$;

  3. iii.

$\Lambda(\tilde{S}_{2}^{-1}S_{2})\in[\alpha_{2},\beta_{2}]$,

where $0<\alpha_{i}\leq 1\leq\beta_{i}$.

We note that to obtain spectral equivalence we seek approximations that yield
values of $\alpha_{i}$ independent of the mesh size and bounded uniformly away
from zero. It is also worth noting that the values of $\alpha_{i}$ and
$\beta_{i}$ are typically not explicitly available. We briefly address this at
the end of this section, following our derivation of the bounds.

To simplify our analyses, we define:

$\displaystyle\tilde{A}^{-1}A=:Q_{0};$ (30a)
$\displaystyle\tilde{S}_{1}^{-1}S_{1}=:Q_{1};$ (30b)
$\displaystyle\tilde{S}_{2}^{-1}S_{2}=:Q_{2}.$ (30c)

To derive the eigenvalue bounds, we consider the split preconditioned matrix

$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}=\begin{bmatrix}\tilde{Q}_{0}&\tilde{B}^{T}&0\\\
\tilde{B}&-\tilde{D}&\tilde{C}^{T}\\\ 0&\tilde{C}&\tilde{E}\end{bmatrix},$
(31)

where $\tilde{Q}_{0}=\tilde{A}^{-1/2}A\tilde{A}^{-1/2}$,
$\tilde{B}=\tilde{S}_{1}^{-1/2}B\tilde{A}^{-1/2}$,
$\tilde{C}=\tilde{S}_{2}^{-1/2}C\tilde{S}_{1}^{-1/2}$,
$\tilde{D}=\tilde{S}_{1}^{-1/2}D\tilde{S}_{1}^{-1/2}$, and
$\tilde{E}=\tilde{S}_{2}^{-1/2}E\tilde{S}_{2}^{-1/2}$. We proceed by bounding
the eigenvalues and singular values of the blocks of
$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}$ and then
applying the results for general double saddle-point matrices presented in
Section 2.2. In order to avoid providing internal eigenvalue bounds equal to
zero (see Remark 1) we assume here that $B$ and $C$ have full row rank.

######  Lemma 3.

When $B$ and $C$ have full row rank, bounds on the eigenvalues/singular values
of the blocks of the matrix
$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}$ are as
follows:

  * •

$\tilde{Q}_{0}$: eigenvalues are in $[\alpha_{0},\beta_{0}]$.

  * •

$\tilde{B}$: singular values are in
$\left[\sqrt{\frac{\alpha_{0}\alpha_{1}}{1+\eta_{D}}},\sqrt{\beta_{0}\beta_{1}}\right]$,
where $\eta_{D}$ is the maximal eigenvalue of $(BA^{-1}B^{T})^{-1}D$.

  * •

$\tilde{C}$: singular values are in
$\left[\sqrt{\frac{\alpha_{1}\alpha_{2}}{1+\eta_{E}}},\sqrt{\beta_{1}\beta_{2}}\right]$,
where $\eta_{E}$ is the maximal eigenvalue of $(CS_{1}^{-1}C^{T})^{-1}E$.

  * •

$\tilde{D}$: eigenvalues are in $[0,\beta_{1}]$.

  * •

$\tilde{E}$: eigenvalues are in $[0,\beta_{2}]$.

###### Proof.

The eigenvalue bounds on $\tilde{Q}_{0}$ follow from the fact that
$\tilde{Q}_{0}$ is similar to $Q_{0}$. For $\tilde{D}$ and $\tilde{E}$, the
lower bounds follow from the fact that $D$ and $E$ are semidefinite (if $D$
and/or $E$ are definite, this bound will be loose, but we use the zero bound
to simplify some results that we present in this section). For the upper bound
on $\tilde{D}$, we note that
$\tilde{D}=\tilde{S}_{1}^{-1/2}D\tilde{S}_{1}^{-1/2}$, which is similar to
$\tilde{S}_{1}^{-1}D=Q_{1}S_{1}^{-1}D$. The eigenvalues of $S_{1}^{-1}D$ are
less than or equal to 1 by Lemma 2, meaning that those of $Q_{1}S_{1}^{-1}D$
are less than or equal to $\beta_{1}$. Analogous reasoning gives the upper
bound for $\tilde{E}$.

We now present the results for $\tilde{B}$. For the upper bound, note that the
matrix
$\tilde{B}\tilde{B}^{T}=\tilde{S}_{1}^{-1/2}B\tilde{A}^{-1}B^{T}\tilde{S}_{1}^{-1/2}$
is similar to
$\tilde{S}_{1}^{-1}B\tilde{A}^{-1}B^{T}=Q_{1}S_{1}^{-1}B\tilde{A}^{-1}B^{T}$.
Because the eigenvalues of $Q_{1}$ are in $[\alpha_{1},\beta_{1}]$, we need
only bound the eigenvalues of $S_{1}^{-1}B\tilde{A}^{-1}B^{T}$. These are the
same as the nonzero eigenvalues of

$\tilde{A}^{-1}B^{T}S_{1}^{-1}B=Q_{0}A^{-1}B^{T}S_{1}^{-1}B.$

The nonzero eigenvalues of $A^{-1}B^{T}S_{1}^{-1}B$ are the same as those of
$S_{1}^{-1}BA^{-1}B^{T}$, which are all less than or equal to 1 by Lemma 2.
Thus, the eigenvalues of $Q_{0}A^{-1}B^{T}S_{1}^{-1}B$ are less than or equal
to $\beta_{0}$, from which we conclude that the eigenvalues of
$\tilde{B}\tilde{B}^{T}$ are less than or equal to $\beta_{0}\beta_{1}$,
giving an upper singular value bound of $\sqrt{\beta_{0}\beta_{1}}$.
Similarly, a lower bound on the eigenvalues of $\tilde{B}\tilde{B}^{T}$ is
given by $\alpha_{0}\alpha_{1}$ times a lower bound on the eigenvalues of
$S_{1}^{-1}BA^{-1}B^{T}$. Because we have assumed that $B$ is full rank,
$BA^{-1}B^{T}$ is invertible, implying that

$S_{1}^{-1}BA^{-1}B^{T}=\left((BA^{-1}B^{T})^{-1}S_{1}\right)^{-1}=\left(I+(BA^{-1}B^{T})^{-1}D\right)^{-1}.$

The stated result then follows because the eigenvalues of
$S_{1}^{-1}BA^{-1}B^{T}$ are greater than or equal to $\frac{1}{1+\eta_{D}}$,
where $\eta_{D}$ is the maximal eigenvalue of $(BA^{-1}B^{T})^{-1}D$. Thus, a
lower bound on the singular values of $\tilde{B}$ is given by the square root
of this value.

The bounds for $\tilde{C}$ are obtained in the same way as those for
$\tilde{B}$. ∎

We can now present bounds on the eigenvalues of
$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}$. We define
three cubic polynomials:

$\displaystyle u(\lambda)$
$\displaystyle=\lambda^{3}+(\beta_{1}-\alpha_{0})\lambda^{2}-\left(\alpha_{0}\beta_{1}+\frac{\alpha_{1}\alpha_{2}}{1+\eta_{E}}+\beta_{0}\beta_{1}\right)\lambda+\frac{\alpha_{0}\alpha_{1}\alpha_{2}}{1+\eta_{E}};$
(32a) $\displaystyle v(\lambda)$
$\displaystyle=\lambda^{3}-(\beta_{0}+\beta_{2})\lambda^{2}+(\beta_{0}\beta_{2}-\beta_{0}\beta_{1}-\beta_{1}\beta_{2})\lambda+2\beta_{0}\beta_{1}\beta_{2};$
(32b) $\displaystyle w(\lambda)$
$\displaystyle=\lambda^{3}+(\beta_{1}-\alpha_{0})\lambda^{2}-(\alpha_{0}\beta_{1}+\beta_{0}\beta_{1}+\beta_{1}\beta_{2})\lambda+\alpha_{0}\beta_{1}\beta_{2}.$
(32c)

These polynomials all have two positive roots and one negative root, by
Corollary 1. As in Section Section 2, we let $u^{-}$ denote the (single)
negative root of a polynomial $u$, and let $u_{\min}^{+}$ and $u_{\max}^{+}$
respectively denote the smallest and largest positive roots.

######  Theorem 6 (Eigenvalue bounds, matrix
$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}$).

When $B$ and $C$ have full row rank, the eigenvalues of
$\tilde{\mathcal{M}}^{-1/2}\mathcal{K}\tilde{\mathcal{M}}^{-1/2}$ are bounded
within the intervals

$\left[w^{-},\frac{\beta_{0}-\sqrt{\beta_{0}^{2}+\frac{4\alpha_{0}\alpha_{1}}{1+\eta_{D}}}}{2}\right]\
\bigcup\ \left[u^{+}_{\min},v^{+}_{\max}\right].$ (33)

###### Proof.

The stated bounds follow from Lemma 3 and Theorem 1. ∎

#### Remark 2.

Obtaining $\eta_{D}$ and $\eta_{E}$ requires computation of eigenvalues of two
matrices related to the Schur complements $S_{1}$ and $S_{2}$:
$(S_{1}-D)^{-1}D=(BA^{-1}B^{T})^{-1}D$ and
$(S_{2}-E)^{-1}E=(CS_{1}^{-1}C^{T})^{-1}E$, respectively. As we have
previously mentioned, in practice when solving (1) the matrices $S_{1}$ and
$S_{2}$ would typically be approximated by sparse and easier-to-invert
matrices, rather than formed and computed explicitly. Therefore, we cannot
expect to compute $\eta_{D}$ and $\eta_{E}$ exactly. Spectral equivalence
relations may be helpful in providing reasonable approximations here. For
example, in common formulations of the Stokes-Darcy problem [7] the matrix $A$
is a discrete Laplacian and $B$ and $C$ are discrete divergence operators (or
scaled variations thereof). In such a case, for certain finite element
discretizations $BA^{-1}B^{T}=S_{1}-D$ is spectrally equivalent to the mass
matrix [8, Section 5.5]. Denote it by $Q$. Then, estimating the maximal
eigenvalue of $\eta_{D}$ amounts to computing an approximation to the maximal
eigenvalue of $Q^{-1}D$, which is computationally straightforward given the
favorable spectral properties of $Q$. If the above-mentioned Schur complement
is approximated by the mass matrix, then it can be shown that
$CS_{1}^{-1}C^{T}=S_{2}-E$ is strongly related to the scalar Laplacian, and
therefore maximal eigenvalue of $\eta_{E}$ would be relatively easy to
approximate as well.

When $D=0$, $\eta_{D}$ and the maximal eigenvalue of $\tilde{D}$ are zero.
Similarly, when $E=0$, $\eta_{E}$ and the maximal eigenvalue of $\tilde{E}$
are zero. In these cases, we can simplify some of the bounds of Theorem 6.
Some of the cubic polynomials that define the bounds will change in these
cases. We will use a bar (e.g., $\bar{u}$) to denote cubic polynomials where
$D=0$, a hat (e.g., $\hat{u}$) to denote the polynomials where $E=0$, and both
(e.g., $\hat{\bar{u}}$) to denote both $D$ and $E$ being zero. We define the
following cubic polynomials:

$\displaystyle\bar{u}(\lambda)$
$\displaystyle=\lambda^{3}-\alpha_{0}\lambda^{2}-\left(\frac{\alpha_{1}\alpha_{2}}{1+\eta_{E}}+\beta_{0}\beta_{1}\right)\lambda+\frac{\alpha_{0}\alpha_{1}\alpha_{2}}{1+\eta_{E}};$
(34a) $\displaystyle\hat{u}(\lambda)$
$\displaystyle=\lambda^{3}+(\beta_{1}-\alpha_{0})\lambda^{2}-(\alpha_{0}\beta_{1}+\alpha_{1}\alpha_{2}+\beta_{0}\beta_{1})\lambda+\alpha_{0}\alpha_{1}\alpha_{2};$
(34b) $\displaystyle\hat{\bar{u}}(\lambda)$
$\displaystyle=\lambda^{3}-\alpha_{0}\lambda^{2}-(\alpha_{1}\alpha_{2}+\beta_{0}\beta_{1})\lambda+\alpha_{0}\alpha_{1}\alpha_{2};$
(34c) $\displaystyle\hat{v}(\lambda)$
$\displaystyle=\lambda^{3}-\beta_{0}\lambda^{2}-(\beta_{0}\beta_{1}+\beta_{1}\beta_{2})\lambda+\beta_{0}\beta_{1}\beta_{2};$
(34d) $\displaystyle\bar{w}(\lambda)$
$\displaystyle=\lambda^{3}-\alpha_{0}\lambda^{2}-(\beta_{0}\beta_{1}+\beta_{1}\beta_{2})\lambda+\alpha_{0}\beta_{1}\beta_{2}.$
(34e)

The following results then follow directly from Theorem 6; the proof is
omitted.

######  Corollary 3.

In the case $D=0,\ E\neq 0$: the eigenvalues of
$\mathcal{M}^{-1/2}\mathcal{K}\mathcal{M}^{-1/2}$ are bounded within the
intervals

$\left[\bar{w}^{-},\frac{\beta_{0}-\sqrt{\beta_{0}^{2}+4\alpha_{0}\alpha_{1}}}{2}\right]\
\bigcup\ \left[\bar{u}^{+}_{\min},v^{+}_{\max}\right].$ (35)

In the case $D\neq 0,\ E=0$: the eigenvalues are bounded in

$\left[w^{-},\frac{\beta_{0}-\sqrt{\beta_{0}^{2}+\frac{4\alpha_{0}\alpha_{1}}{1+\eta_{D}}}}{2}\right]\
\bigcup\ \left[\hat{u}^{+}_{\min},\hat{v}^{+}_{\max}\right].$ (36)

In the case $D=0,\ E=0$: the eigenvalues are bounded in

$\left[\bar{w}^{-},\frac{\beta_{0}-\sqrt{\beta_{0}^{2}+4\alpha_{0}\alpha_{1}}}{2}\right]\
\bigcup\ \left[\hat{\bar{u}}^{+}_{\min},\hat{v}^{+}_{\max}\right].$ (37)

#### Remark 3.

There is some looseness in the bounds of Theorem 6 when applied to the exact
preconditioning case (i.e., $\alpha_{i}=\beta_{i}=1$). This is a consequence
of the fact that Theorem 6 is based on the bounds for unpreconditioned
matrices, which consider each matrix block individually, as opposed to the
energy estimates approach in Section 3, which fully exploits the relationships
between the blocks of the preconditioned matrix. For the extremal (lower
negative and upper positive) bounds, the looseness is minor: Theorem 6 gives a
bound of 2 for the positive eigenvalues and approximately $-1.9$ on the
negative eigenvalues, while we know from Theorem 4 that tight bounds are
approximately 1.8 and $-1.6$, respectively. For the interior bounds, however,
we note that the bounds may be quite loose if $\eta_{D}$ or $\eta_{E}$ is very
large. Fortunately, this is often not a concern in practical settings as
having a large $\eta_{D}$ or $\eta_{E}$ generally implies that the spectral
norm of $D$ or $E$ is large relative to that of $BA^{-1}B^{T}$ or
$CS_{1}^{-1}C^{T}$, respectively. Often $D$ and/or $E$ are regularization
terms, which tend to have a fairly small norm. Nonetheless, it should be
acknowledged that the interior bounds may be pessimistic for some problems.

#### Remark 4.

The values of $\alpha_{i}$ and $\beta_{i}$ are rarely available in practical
applications; they are typically coercivity constants or related quantities,
which are proven to be independent of the mesh size but are not known
explicitly. The values of $\beta_{i}$ should not typically generate a
difficulty, and approximating them by $1$ or a value close to 1 should provide
a reasonable approximation for the bounds.

For the $\alpha_{i}$, let us provide some (partial) observations. For the
solutions of the quadratic equations in Corollary 3, we can use the Taylor
approximation $\sqrt{1+x}\lessapprox 1+\frac{x}{2}$ for $0<x\ll 1$ to conclude
that

$\frac{\beta_{0}-\sqrt{\beta_{0}^{2}+4\alpha_{0}\alpha_{1}}}{2}\gtrapprox-\frac{\alpha_{0}\alpha_{1}}{\beta_{0}}.$

This means that if $\alpha_{0}$ and $\alpha_{1}$ are small, then the above
displayed expression would be a valid (albeit slightly less tight) upper
negative bound, and if the $\alpha_{i}$ are uniformly bounded away from zero,
then so is the bound.

##  5 Numerical experiments

In this section we consider two slightly different variants of a Poisson
control problem as in [30]. First, we consider a distributed control problem
with Dirichlet boundary conditions:

$\displaystyle\min_{u,f}\frac{1}{2}$
$\displaystyle||u-\hat{u}||_{L_{2}(\Omega)}^{2}+\frac{\beta}{2}||f||_{L_{2}(\Omega)}^{2}$
(38a) $\displaystyle\textrm{s.t.}\ \ -\nabla^{2}u$ $\displaystyle=f\ \
\textrm{in }\Omega,$ (38b) $\displaystyle u$ $\displaystyle=g\ \ \textrm{on
}\partial\Omega,$ (38c)

where $u$ is the state, $\hat{u}$ is the desired state, $0<\beta\ll 1$ is a
regularization parameter, $f$ is the control, and $\Omega$ is the domain with
boundary $\partial\Omega$. Upon discretization, we obtain the system

$\underbrace{\begin{bmatrix}M&K&0\\\ K&0&-M\\\ 0&-M&\beta
M\end{bmatrix}}_{=:\mathcal{K}}\begin{bmatrix}u_{h}\\\ \lambda\\\
f_{h}\end{bmatrix}=\begin{bmatrix}b\\\ d\\\ 0\end{bmatrix},$ (39)

where $M$ is a symmetric positive definite mass matrix and $K$ is a symmetric
positive definite discrete Laplacian. All blocks of $\mathcal{K}$ are square
(i.e., $n=m=p$).

As a second experiment we consider a boundary control problem:

$\displaystyle\min_{u,f}\frac{1}{2}$
$\displaystyle||u-\hat{u}||_{L_{2}(\Omega)}^{2}+\frac{\beta}{2}||g||_{L_{2}(\Omega)}^{2}$
(40a) $\displaystyle\textrm{s.t.}\ \ -\nabla^{2}u$ $\displaystyle=0\ \
\textrm{in }\Omega,$ (40b) $\displaystyle\frac{\partial u}{\partial n}$
$\displaystyle=g\ \ \textrm{on }\partial\Omega,$ (40c)

which after discretization yields the linear system

$\underbrace{\begin{bmatrix}M&K&0\\\ K&0&-E\\\ 0&-E^{T}&\beta
M_{b}\end{bmatrix}}_{=:\mathcal{K}_{\partial}}\begin{bmatrix}u_{h}\\\
\lambda\\\ g_{h}\end{bmatrix}=\begin{bmatrix}b_{\partial}\\\ d_{\partial}\\\
0\end{bmatrix},$ (41)

where $M_{b}\in\mathbb{R}^{n_{b}\times n_{b}}$ (with $n_{b}<n$) is a boundary
mass matrix. Thus, in the distributed control problem the mass matrix in the
(2,3)/(3,2)-block is square and in the boundary control problem we consider a
rectangular version of it.

In all experiments that follow we set $\Omega$ to be the unit square. We use
uniform Q1 finite elements and set $\beta=10^{-3}$. The MATLAB code of Rees
[29] was used to generate the linear systems.

###  5.1 Eigenvalues of unpreconditioned matrices

Here we compare the eigenvalues of $\mathcal{K}$ (39) and
$\mathcal{K}_{\partial}$ (41) to the eigenvalue bounds predicted by Theorem 1.
We use MATLAB’s eigs/svds functions to compute the minimum and maximum
eigenvalues/singular values of the matrix blocks $M,K,M_{b}$, and $E$.

We note that for the distributed control matrix $\mathcal{K}$, all blocks are
square and the (3,3)-block is positive definite; therefore, we can also use
our results to obtain bounds on the re-ordered matrix

$\mathcal{K}_{\textrm{flip}}=\begin{bmatrix}\beta M&-M&0\\\ -M&0&K\\\
0&K&M\end{bmatrix}.$ (42)

Both orderings generate the same extremal bounds but different interior
bounds.

(a) Positive eigenvalues of $\mathcal{K}$

(b) Negative eigenvalues of $\mathcal{K}$

(c) Positive eigenvalues of $\mathcal{K}_{\partial}$

(d) Negative eigenvalues of $\mathcal{K}_{\partial}$

Figure 2: Largest and smallest positive (left) and negative (right)
eigenvalues of the distributed control matrix $\mathcal{K}$ (top) and the
boundary control matrix $\mathcal{K}_{\partial}$ (bottom). Blue circles
indicate the eigenvalues, and black lines the bounds given by Theorem 1. For
$\mathcal{K}$, the dashed lines indicate the bounds obtained by applying
Theorem 1 to the reordered matrix $\mathcal{K}_{\textrm{flip}}$.

Comparisons of the predicted eigenvalue bounds to the actual eigenvalues are
shown in Section 2. In the distributed control case, we show the bounds
obtained for both the original matrix $\mathcal{K}$ and those for the
reordered matrix $\mathcal{K}_{\textrm{flip}}$. In all cases, the bounds for
the extremal eigenvalues are quite tight: this is because the $K$ block has
the largest eigenvalues ($O(1)$ compared to $O(h^{2})$ for the others – see
[8, Proposition 1.29 and Theorem 1.32]). Therefore, the largest positive and
smallest negative eigenvalue of both $\mathcal{K}$ and
$\mathcal{K}_{\partial}$ tend towards $\mu_{\max}^{K}$ and $-\mu_{\max}^{K}$,
respectively. Referring to the proofs of the extremal bounds in Theorem 1, the
$3\times 3$ matrices $R$ will contain a $\mu_{\max}^{K}$ (or
$-\mu_{\max}^{K}$) term, with other lower-order terms. This means that the
extremal eigenvalues of $R$ (and therefore the predicted eigenvalue bounds)
will also be close to $\pm\mu_{\max}^{K}$.

It is more difficult to capture the interior bounds. In the distributed
control case, each ordering (either the original ordering with $M$ in the
leading block or the flipped ordering with $\beta M$ in the leading block)
gives one bound that is quite tight and one that is loose. In the boundary
control case, both interior bounds are loose.

###  5.2 Eigenvalues of preconditioned matrices

We now consider preconditioning strategies for PDE-constrained optimization.
We examine eigenvalue bounds with both exact and approximate Schur
complements.

For the distributed control problem, we work with the reordered matrix
$\mathcal{K}_{\textrm{flip}}$ (42). The Schur complement preconditioner for
$\mathcal{K}_{\textrm{flip}}$ is

$\mathcal{M}=\begin{bmatrix}\beta M&0&0\\\ 0&\frac{1}{\beta}M&0\\\ 0&0&M+\beta
KM^{-1}K\end{bmatrix}.$ (43)

The first and second blocks are mass matrices, which are cheap to invert, so
we leave these terms as they are. For the second Schur complement
$S_{2}=M+\beta KM^{-1}K$, we use the approximation
$\tilde{S}_{2}=\left(M+\sqrt{\beta}K\right)M^{-1}\left(M+\sqrt{\beta}K\right)$
proposed by Pearson and Wathen [26] to obtain the preconditioner:

$\tilde{\mathcal{M}}=\begin{bmatrix}\beta M&0&0\\\ 0&\frac{1}{\beta}M&0\\\
0&0&\left(M+\sqrt{\beta}K\right)M^{-1}\left(M+\sqrt{\beta}K\right)\end{bmatrix}.$
(44)

Per [26, Theorem 4], the eigenvalues of $\tilde{S}_{2}^{-1}S_{2}$ satisfy
$\Lambda(\tilde{S}_{2}^{-1}S_{2})\in\left[\frac{1}{2},1\right]$. Thus
$\tilde{\mathcal{M}}$ satisfies Theorem 5 with
$\alpha_{0}=\beta_{0}=\alpha_{1}=\beta_{1}=1$, $\alpha_{2}=\frac{1}{2}$, and
$\beta_{2}=1$.

(a) $\mathcal{M}^{-1}\mathcal{K}_{\textrm{flip}}$

(b) $\tilde{\mathcal{M}}^{-1}\mathcal{K}_{\textrm{flip}}$

Figure 3: Eigenvalue plots for reordered distributed Poisson control matrix
$\mathcal{K}_{\textrm{flip}}$, with exact preconditioner $\mathcal{M}$ (left)
and approximate preconditioner $\tilde{\mathcal{M}}$ (right). Eigenvalues are
shown by the blue circles; eigenvalue bounds (from Theorem 3 on left and
Corollary 3 on right) are shown by lines.

Plots of the preconditioned eigenvalues for $h=2^{-4}$ are shown in Section 3.
The value $\eta_{E}$, defined as the maximal eigenvalue of

$(CS_{1}^{-1}C^{T})^{-1}E=\beta(KM^{-1}K)^{-1}M,$

is approximately $2.6\times 10^{-7}$. We note that for 2D problems with
uniform Q1 finite element discretizations the value $\eta_{E}$ is $O(\beta
h^{4})$, and will thus be small in general.

Comparing Figures 3(a) and 3(b), we notice that the bounds on the negative
eigenvalue bounds do not change when we use the approximate Schur complement,
but the lower positive bound becomes smaller (from 0.4450 with the exact Schur
complement to 0.2929 for the approximate Schur complement) and the upper
positive bound becomes larger (1.8019 in the exact case and 2 in the
approximate case). We note that the eigenvalues appear to be very close to the
predicted bounds except for the upper positive eigenvalues of
$\tilde{\mathcal{M}}^{-1}\mathcal{K}$. As discussed in Remark 3, this kind of
minor looseness in the upper bound can happen when we have highly accurate
Schur complement approximations (as in this case, where we are exactly
inverting the two $M$ blocks).

For the matrix $\mathcal{K}_{\partial}$ arising from the boundary control
problem, the Schur complement preconditioner is

$\mathcal{M}_{\partial}=\begin{bmatrix}M&0&0\\\ 0&KM^{-1}K&0\\\ 0&0&\beta
M_{b}+E^{T}\left(KM^{-1}K\right)^{-1}E\end{bmatrix}.$ (45)

In practice, the first Schur complement $KM^{-1}K$ can be inverted
approximately with, for example, a multigrid method. For the second Schur
complement, we note that for the 2D Poisson control problem on a uniform Q1
grid, the eigenvalues of $K$ are between $O(h^{2})$ and $O(1)$, while those of
the mass matrices are all $O(h^{2})$. Therefore, the term $\beta M_{b}$ will
dominate $E^{T}\left(KM^{-1}K\right)^{-1}E$ for all but very small values of
$\beta$ (which are not commonly used in practice). Therefore, an approximate
preconditioner for $\mathcal{K}_{\partial}$ is given by

$\tilde{\mathcal{M}}_{\partial}=\begin{bmatrix}M&0&0\\\ 0&KM^{-1}K&0\\\
0&0&\beta M_{b}\end{bmatrix}.$ (46)

We note that this preconditioner is presented in [30], though it is derived
there from the block-$2\times 2$ formulation of $\mathcal{K}$.

(a) $\mathcal{M}_{\partial}^{-1}\mathcal{K}_{\partial}$

(b) $\mathcal{M}_{\partial}^{-1}\mathcal{K}_{\partial}$ (blue);
$\tilde{\mathcal{M}}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$ (red)

Figure 4: Eigenvalue plots for the boundary control problem
$\mathcal{K}_{\partial}$, with exact preconditioner $\mathcal{M}$ and
approximate preconditioner $\tilde{\mathcal{M}}$. On left, eigenvalue bounds
(from Theorem 3) are shown by horizontal lines. Right: all but the single
smallest and largest eigenvalues for both preconditioners; values for
$\mathcal{M}$ are in blue, and those for $\tilde{\mathcal{M}}$ are in red.

Plots of the preconditioned eigenvalues for $h=2^{-4}$ are shown in Figure 4.
For $\tilde{\mathcal{M}}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$, we note
that there is a single large-magnitude positive eigenvalue and a large-
magnitude negative eigenvalue, whose absolute values are nearly 200, and they
make it difficult to see how the other eigenvalues compare to to those of
$\mathcal{M}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$. Thus in Figure 4(b),
we omit the largest and smallest eigenvalues of
$\tilde{\mathcal{M}}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$ and overlay
the others on the corresponding eigenvalues of
$\mathcal{M}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$. We notice that most
other eigenvalues of
$\tilde{\mathcal{M}}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$ remain close
to those of $\mathcal{M}_{\partial}^{-1}\mathcal{K}_{{}_{\partial}}$. We also
note that the performance of this preconditioner depends on $\beta$: in
particular, the performance of the Schur complement approximation deteriorates
for small $\beta$. We refer to [26, 30] for further discussion of this.

It is evident from Figure 4(a) that our bounds are tight and effective. Unlike
in the boundary control example, the Schur complement approximation
$\tilde{S}_{2}$ used here does not have $\beta$\- and $h$-independent
constants $\alpha_{2},\beta_{2}$ such that
$\Lambda(\tilde{S}_{2}^{-1}S_{2})\in[\alpha_{2},\beta_{2}]$, so our analyses
on Schur complement approximations in Section Section 4 are difficult to
apply. Nonetheless, we observe from Figure 4(b) that most of the eigenvalues
of $\tilde{\mathcal{M}}_{\partial}^{-1}\mathcal{K}_{\partial}$ remain very
close to those of $\mathcal{M}_{\partial}^{-1}\mathcal{K}_{\partial}$. Thus,
we see that the eigenvalue bounds for the “ideal” Schur complement
preconditioner may still be of use, provided that we have an effective
approximation of the Schur complement.

##  6 Conclusions

The increasing importance of double saddle-point systems requires attention to
spectral properties of the matrices involved. We have shown that energy
estimates are an effective tool for obtaining eigenvalue bounds in this case.
The assumptions we make are rather general, and the analysis covers a large
class of problems.

There are several directions for potential future work. Following Remark
LABEL:rem:looseness, specific assumptions on the magnitudes of the norms of
the matrices $D$ and $E$ may yield additional results and insights on the
bounds. The rank structure of the blocks may also have a significant effect,
and it may be useful to eliminate the positive definiteness requirement of $A$
and/or consider rank-deficient $B$ and $C$. It may also be useful to consider
nonsymmetric double saddle-point systems, and block triangular rather than
block diagonal preconditioners, although when symmetry is lost eigenvalue
bounds may not be a sufficient tool for predicting convergence rates of Krylov
subspace solvers.

The tightness of our bounds indicates that spectral analysis is beneficial in
capturing the properties of the matrices involved. This, in turn, makes it
possible to effectively predict the convergence rate of Krylov subspace
solvers for this important class of problems.

## References

  * [1] J. H. Adler, T. R. Benson, E. C. Cyr, P. E. Farrell, S. P. MacLachlan, and R. S. Tuminaro.  Monolithic multigrid methods for magnetohydrodynamics.  SIAM Journal on Scientific Computing, 43(5):S70–S91, 2021.
  * [2] F. Ali Beik and M. Benzi.  Iterative methods for double saddle point systems.  SIAM Journal on Matrix Analysis and Applications, 39(2):902–921, 2018.
  * [3] A. Beigl, J. Sogn, and W. Zulehner.  Robust preconditioners for multiple saddle point problems and applications to optimal control problems.  SIAM Journal on Matrix Analysis and Applications, 41(4):1590–1615, 2020.
  * [4] M. Benzi.  Preconditioning techniques for large linear systems: a survey.  J. Comput. Phys., 182(2):418–477, 2002.
  * [5] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numerica, 14:1–137, 2005.
  * [6] M. Cai, G. Ju, and J. Li.  Schur complement based preconditioners for twofold and block tridiagonal saddle point problems, 2021.
  * [7] M. Cai, M. Mu, and J. Xu.  Preconditioning techniques for a mixed Stokes/Darcy model in porous media applications.  Journal of Computational and Applied Mathematics, 233(2):346 – 355, 2009.
  * [8] H. C. Elman, D. J. Silvester, and A. J. Wathen.  Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics.  Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2005.
  * [9] M. Ferronato, A. Franceschini, C. Janna, N. Castelletto, and H. A. Tchelepi.  A general preconditioning framework for coupled multiphysics problems with application to contact- and poro-mechanics.  Journal of Computational Physics, 398:108887, 2019.
  * [10] G. Gatica and N. Heuer.  A dual-dual formulation for the coupling of mixed-FEM and BEM in hyperelasticity.  SIAM Journal on Numerical Analysis, 38:380–400, July 2000.
  * [11] N. Gould and V. Simoncini.  Spectral analysis of saddle point matrices with indefinite leading blocks.  SIAM J. Matrix Analysis Applications, 31:1152–1171, 01 2009.
  * [12] N. I. Gould.  On practical conditions for the existence and uniqueness of solutions to the general equality quadratic programming problem.  Mathematical Programming, 32:90–99, 1985.
  * [13] C. Greif, E. Moulding, and D. Orban.  Bounds on eigenvalues of matrices arising from interior-point methods.  SIAM Journal on Optimization, 24(1):49–83, 2014.
  * [14] R. Hiptmair.  Operator preconditioning.  Comput. Math. Appl., 52(5):699–706, 2006.
  * [15] K. E. Holter, M. Kuchta, and K.-A. Mardal.  Robust preconditioning for coupled Stokes-Darcy problems with the Darcy problem in primal form.  Comput. Math. Appl., 91:53–66, 2021.
  * [16] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Cambridge University Press, 1 edition, 1985.
  * [17] I. C. F. Ipsen.  A note on preconditioning nonsymmetric matrices.  SIAM J. Sci. Comput., 23(3):1050–1051, 2001.
  * [18] D. E. Keyes, L. C. McInnes, C. S. Woodward, W. Gropp, E. Myra, M. Pernice, J. Bell, J. Brown, A. Clo, J. M. Connors, E. M. Constantinescu, D. J. Estep, K. Evans, C. Farhat, A. Hakim, G. Hammond, G. Hansen, J. Hill, T. Isaac, X. Jiao, K. E. Jordan, D. K. Kaushik, E. Kaxiras, A. E. Koniges, K. Lee, A. Lott, Q. Lu, J. Magerlein, R. Maxwell, M. McCourt, M. Mehl, R. P. Pawlowski, A. P. Randles, D. R. Reynolds, B. Rivière, U. Rüde, T. D. Scheibe, J. N. Shadid, B. Sheehan, M. S. Shephard, A. R. Siegel, B. Smith, X. Tang, C. Wilson, and B. I. Wohlmuth.  Multiphysics simulations: Challenges and opportunities.  Int. J. High Perform. Comput. Appl., 27(1):4–83, 2013.
  * [19] U. Langer, G. Of, O. Steinbach, and W. Zulehner.  Inexact data-sparse boundary element tearing and interconnecting methods.  SIAM Journal on Scientific Computing, 29(1):290–314, 2007.
  * [20] K.-A. Mardal and R. Winther.  Preconditioning discretizations of systems of partial differential equations.  Numer. Linear Algebra Appl., 18(1):1–40, 2011.
  * [21] B. Morini, V. Simoncini, and M. Tani.  Spectral estimates for unreduced symmetric KKT systems arising from interior point methods.  Numer. Linear Algebra Appl., 23:776–800, 2016.
  * [22] M. F. Murphy, G. H. Golub, and A. J. Wathen.  A note on preconditioning for indefinite linear systems.  SIAM J. Sci. Comput., 21(6):1969–1972, 2000.
  * [23] C. C. Paige and M. A. Saunders.  Solutions of sparse indefinite systems of linear equations.  SIAM J. Numer. Anal., 12(4):617–629, 1975.
  * [24] J. W. Pearson and J. Pestana.  Preconditioners for Krylov subspace methods: an overview.  GAMM-Mitt., 43(4):e202000015, 35, 2020.
  * [25] J. W. Pearson and A. Potschka.  A note on symmetric positive definite preconditioners for multiple saddle-point systems, 2021.
  * [26] J. W. Pearson and A. J. Wathen.  A new approximation of the Schur complement in preconditioners for PDE-constrained optimization.  Numerical Linear Algebra with Applications, 19(5):816–829, 2012\.
  * [27] J. Pestana and A. J. Wathen.  Natural preconditioning and iterative methods for saddle point systems.  SIAM Rev., 57(1):71–91, 2015.
  * [28] A. Ramage and E. C. Gartland, Jr.  A preconditioned nullspace method for liquid crystal director modeling.  SIAM J. Sci. Comput., 35(1):B226–B247, 2013.
  * [29] T. Rees.  Github - tyronerees/poisson-control.  https://github.com/tyronerees/poisson-control, 2010.  Accessed: 2021-10-22.
  * [30] T. Rees, H. Dollar, and A. Wathen.  Optimal solvers for PDE-constrained optimization.  SIAM Journal on Scientific Computing, 32(1):271–298, 2010.
  * [31] S. Rhebergen, G. Wells, A. Wathen, and R. Katz.  Three-field block preconditioners for models of coupled magma/mantle dynamics.  SIAM Journal on Scientific Computing, 37(5):A2270–A2294, 2015.
  * [32] M. Rozložník.  Saddle-point problems and their iterative solution.  Nečas Center Series. Birkhäuser/Springer, Cham, 2018.
  * [33] D. Ruiz, A. Sartenaer, and C. Tannier.  Refining the lower bound on the positive eigenvalues of saddle point matrices with insights on the interactions between the blocks.  SIAM Journal on Matrix Analysis and Applications, 39(2):712–736, 2018.
  * [34] T. Rusten and R. Winther.  A preconditioned iterative method for saddlepoint problems.  SIAM J. Matrix Anal. Appl., 13:887–904, 1992.
  * [35] Y. Saad.  Iterative methods for sparse linear systems.  Society for Industrial and Applied Mathematics, Philadelphia, PA, second edition, 2003.
  * [36] D. Silvester and A. Wathen.  Fast iterative solution of stabilised Stokes systems part II: Using general block preconditioners.  SIAM Journal on Numerical Analysis, 31(5):1352–1367, 1994.
  * [37] J. Sogn and W. Zulehner.  Schur complement preconditioners for multiple saddle point problems of block tridiagonal form with application to optimization problems.  IMA Journal of Numerical Analysis, 39(3):1328–1359, 05 2018.
  * [38] A. J. Wathen.  Preconditioning.  Acta Numer., 24:329–376, 2015.