Eigenvalue Bounds for Symmetric, Multiple Saddle-Point Matrices with SPD Preconditioners

1. Introduction

We consider the iterative solution of a block tridiagonal multiple saddle-point linear system $\mathcal{A}x=b$, where

$$\mathcal{A}=\left[\begin{array}{ccccc}{A}_0& B_1^{\top }& 0& \dots & 0\\ B_1& -A_1& B_2^{\top }& \ddots & ⋮\\ 0& B_2& A_2& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & B_N^{\top }\\ 0& \dots & 0& B_N& {(-1)}^N{A}_N\end{array}\right].$$

We assume that $A_0\in {\mathbb{R}}^{n_0\times n_0}$ is symmetric positive definite, all other square block matrices $A_k\in {\mathbb{R}}^{n_k\times n_k}$ are symmetric positive semi-definite and $B_k\in {\mathbb{R}}^{n_k\times n_{k-1}}$ have full rank (for $k=1,\dots ,N$). We assume also that $n_k\le n_{k-1}$ for all k. These conditions are sufficient to ensure the invertibility of $\mathcal{A}$.

Linear systems involving matrix $\mathcal{A}$ arise mostly with $N=2$ (double saddle–point systems) in many scientific applications including magma–mantle dynamics [3], liquid crystal director modeling [4] or in the coupled Stokes–Darcy problem [5,6,7,8], and the preconditioning of such linear systems has been considered in [9,10,11]. In particular, block diagonal preconditioners for matrix $\mathcal{A}$ have been studied in [12,13,14,15,16].

Multiple saddle-point linear systems with $N>2$ have recently attracted the attention of a number of researchers. Such systems often arise from modeling multiphysics processes, i.e. the simultaneous simulation of different aspects of physical systems and the interactions among them. Preconditioning of the $4\times 4$ saddle-point linear system ($N=3$), although with a different block structure, has been addressed in [17] to solve a class of optimal control problems. A practical preconditioning strategy for multiple saddle–point linear systems, based on sparse approximate inverses of the diagonal blocks of the block diagonal Schur complement preconditioning matrix, is proposed in [18], for the solution of coupled poromechanical models and the mechanics of fractured media. Theoretical analysis of such preconditioners has been carried on in [19]. In [17], two $5\times 5$ multiple saddle-point systems arising from optimal control problems constrained with either the heat or the wave equation are addressed, and iteratively solved by robust block diagonal preconditioners.

In this work we consider the preconditioner proposed in [1] (PP preconditioner in short) to develop a spectral analysis of the preconditioned matrix for an arbitrary number of blocks and in presence of inexact Schur complements, in the line of [15] in which the block diagonal preconditioner with exact Schur complements is considered. We will define a sequence of polynomials by a three-term recurrence, whose extremal roots characterize the eigenvalues of the preconditioned matrix in terms of the nonnegative eigenvalues of the preconditioned Schur complements. Using a similar technique, we will extend results, obtained in [2] for a double saddle-point, to multiple saddle-point linear system.

Arguably the most prominent Krylov subspace methods for solving the linear system $\mathcal{A}\mathit{x}=\mathit{b}$ are preconditioned variants of MINRES [20] and GMRES [21]. In contrast to GMRES, the previously-discovered MINRES algorithm can explicitly exploit the symmetry of $\mathcal{A}$. As a consequence, MINRES features a three-term recurrence relation, which is beneficial for its implementation (low memory requirements because subspace bases need not be stored) and its purely eigenvalue-based convergence analysis (via the famous connection to orthogonal polynomials; see [22,23]). Specifically, if the eigenvalues of the preconditioned matrix are contained within $[{\rho }_l^{-},{\rho }_{\mathit{u}}^{-}]\cup [{\rho }_l^{+},{\rho }_{\mathit{u}}^{+}]$, for ${\rho }_l^{-}<{\rho }_{\mathit{u}}^{-}<0<{\rho }_l^{+}<{\rho }_{\mathit{u}}^{+}$ such that ${\rho }_{\mathit{u}}^{+}-{\rho }_l^{+}={\rho }_{\mathit{u}}^{-}-{\rho }_l^{-}$, then at iteration k the Euclidean norm of the preconditioned residual $r_k$ satisfies the bound

$${\displaystyle \frac{\parallel {\mathit{r}}_k\parallel }{\parallel {\mathit{r}}_0\parallel }}\le 2{\left({\displaystyle \frac{\sqrt{\left|{\rho }_l^{-}{\rho }_{\mathit{u}}^{+}\right|}-\sqrt{\left|{\rho }_{\mathit{u}}^{-}{\rho }_l^{+}\right|}}{\sqrt{\left|{\rho }_l^{-}{\rho }_{\mathit{u}}^{+}\right|}+\sqrt{\left|{\rho }_{\mathit{u}}^{-}{\rho }_l^{+}\right|}}}\right)}^{\lfloor k/2\rfloor }.$$

The paper is structured as follows: In Section 2 we describe the preconditioner and write the eigenvalue problem for the preconditioned matrix. In Section 3, a sequence of polynomials is defined by recurrence, and the connection of these with the eigenvalues of the preconditioned matrix is stated. Then, in Section 4, we characterize the extremal roots of such polynomials and the dependence of such roots with two vectors of parameters, containing information on the inexact Schur complements. Section 5 is devoted to comparing the bounds with the exactly computed eigenvalues on a large set of synthetic experiments for the number of blocks ranging from 3 to 5. In Section 6 we present as a realistic test case, the 3D Mixed Finite Element discretization of the Biot model, which gives raise to a double saddle-point linear system, to again validate the bounds and to show the good performance of the PP preconditioner, also in comparison with the most known block diagonal preconditioner. Section 7 draws some conclusions.

2. The Preconditioner

Consider the exact block factorization of matrix $\mathcal{A}=U^{\top }{D}^{-1}U$, where

$$D=\left[\begin{array}{ccccc}{A}_0& & & & \\ & -S_1& & & \\ & & S_2& & \\ & & & \ddots & \\ & & & & {(-1)}^N{S}_N\end{array}\right],U=\left[\begin{array}{ccccc}{A}_0& B_1^{\top }& & & \\ & -S_1& B_2^{\top }& & \\ & & S_2& \ddots & \\ & & & \ddots & B_N^{\top }\\ & & & & {(-1)}^N{S}_N\end{array}\right],$$

(1)

with $S_k=A_k+B_k{S}_{k-1}^{-1}{B}_k^{\top },k=1,\dots ,N$. The preconditioner is initially obtained by removing the minus signs in D, thus obtaining a symmetric positive definite matrix $P=U^{\top }\left|D^{-1}\right|U$, and then, in view of its practical application, by approximating the Schur complements with their inexact counterparts. The final expression of the preconditioner is therefore $\mathcal{P}={\mathcal{P}}_U^{\top }{\mathcal{P}}_D^{-1}{P}_U,$ where

$${\mathcal{P}}_D=\left[\begin{array}{ccccc}{\widehat{S}}_0& & & & \\ & {\widehat{S}}_1& & & \\ & & {\widehat{S}}_2& & \\ & & & \ddots & \\ & & & & {\widehat{S}}_N\end{array}\right],{\mathcal{P}}_U=\left[\begin{array}{ccccc}{\widehat{S}}_0& B_1^{\top }& & & \\ & -{\widehat{S}}_1& B_2^{\top }& & \\ & & {\widehat{S}}_2& \ddots & \\ & & & \ddots & B_N^{\top }\\ & & & & {(-1)}^N{\widehat{S}}_N\end{array}\right],$$

(2)

with

$$\begin{array}{ccc}& & {\widehat{S}}_0\approx A_0\hfill \\ {\tilde{S}}_k\hfill & =A_k+B_k{\widehat{S}}_{k-1}^{-1}{B}_k^{\top }& {\widehat{S}}_k\approx {\tilde{S}}_k.\hfill \end{array}$$

This preconditioner has been proposed by Pearson and Potschka in [1], in the framework of PDE-constrained optimization. Since $\mathcal{P}$ is symmetric and positive definite, this allows the use of the MINRES iterative method.

Finding the eigenvalues of ${\mathcal{P}}^{-1}\mathcal{A}$ is equivalent to solving

$${\mathcal{P}}_D^{-1/2}\mathcal{A}{\mathcal{P}}_D^{-1/2}\mathit{u}=\lambda {\mathcal{P}}_D^{-1/2}{P}_L^{\top }{\mathcal{P}}_D^{-1/2}{\mathcal{P}}_D^{-1/2}{P}_U{\mathcal{P}}_D^{-1/2}\mathit{u},\mathit{u}=\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\\ \dots \\ {\mathit{u}}_{N+1}\end{array}\right].$$

Exploiting the block components of this generalized eigenvalue problem, we obtain

$$\begin{array}{ccc}\hfill & & \left[\begin{array}{cccccc}{E}_0& R_1^{\top }& & & & \\ R_1& -E_1& R_2^{\top }& & & \\ 0& R_2& E_2& R_3^{\top }& & \\ & \dots & \dots & \dots \\ & & R_{N-1}& {(-1)}^{N-1}{E}_{N-1}& R_N^{\top }\\ & & & R_N& {(-1)}^N{E}_N\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\\ {\mathit{u}}_3\\ \dots \\ {\mathit{u}}_N\\ {\mathit{u}}_{N+1}\end{array}\right]=\hfill \\ \hfill & =& \lambda \left[\begin{array}{cccccc}I& R_1^{\top }& & & & \\ R_1& I+R_1{R}_1^{\top }& -R_2^{\top }& & & \\ 0& -R_2& I+R_2{R}_2^{\top }& R_3^{\top }& & \\ & \dots & \dots & \dots \\ & & {(-1)}^N{R}_{N-1}& I+R_{N-1}{R}_{N-1}^{\top }& {(-1)}^{N+1}{R}_N^{\top }\\ & & & {(-1)}^{N+1}{R}_N& I+R_N{R}_N^{\top }\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\\ {\mathit{u}}_3\\ \dots \\ {\mathit{u}}_N\\ {\mathit{u}}_{N+1}\end{array}\right]\hfill \end{array}$$

(3)

where

$$\begin{array}{ccc}\hfill R_k& =& {\widehat{S}}_k^{-1/2}{B}_k{\widehat{S}}_{k-1}^{-1/2},k=1,\dots ,N;E_k={\widehat{S}}_k^{-1/2}{A}_k{\widehat{S}}_k^{1/2},k=0,\dots ,N\hfill \\ \hfill R_k^{\top }{R}_k+E_k& =& {\widehat{S}}_k^{-1/2}{\tilde{S}}_k{\widehat{S}}_k^{-1/2}\equiv {\overline{S}}_k.\hfill \end{array}$$

Componentwise, we write the eigenvalue problem (3) as

$$\begin{array}{ccc}\hfill \mathbf{0}& =& (E_{\mathbf{0}}-\lambda I){\mathit{u}}_1+(1-\lambda )R_1^{\top }{\mathit{u}}_2,\hfill \\ \hfill \mathbf{0}& =& (1-\lambda )R_1{\mathit{u}}_1+(-E_1-\lambda (I+R_1{R}_1^{\top })){\mathit{u}}_2+(1+\lambda )R_2^{\top }{\mathit{u}}_3,\hfill \\ \hfill \mathbf{0}& =& (1+\lambda )R_2{\mathit{u}}_2+(E_2-\lambda (I+R_2{R}_2^{\top })){\mathit{u}}_3+(1-\lambda )R_3^{\top }{\mathit{u}}_4,\hfill \\ \hfill & & \dots \dots \dots \hfill \\ \hfill \mathbf{0}& =& (1+{(-1)}^{N-1}\lambda )R_{N-1}{\mathit{u}}_{N-1}+({(-1)}^{N-1}{E}_{N-1}-\lambda (I-R_{N-1}{R}_{N-1}^{\top })){\mathit{u}}_N\hfill \\ \hfill & & +(1+{(-1)}^N\lambda )R_N^{\top }{\mathit{u}}_{N+1},\hfill \\ \hfill \mathbf{0}& =& (1+{(-1)}^N\lambda )R_N{\mathit{u}}_N+({(-1)}^N{E}_N-\lambda (I+R_N{R}_N^{\top })){\mathit{u}}_{N+1}.\hfill \end{array}$$

(4)

The matrices $R_k{R}_k^{\top }$ are all symmetric positive definite. We define two indicators ${\gamma }_E^{\left(k\right)}$ and ${\gamma }_R^{\left(k\right)}$ using the Rayleigh quotient

$$\begin{array}{ccc}{\alpha }_E^{\left(k\right)}\equiv {\lambda }_{min}\left(E_k\right),\hfill & {\beta }_E^{\left(k\right)}\equiv {\lambda }_{max}\left(E_k\right),\hfill & {\gamma }_E^{\left(k\right)}\left(w\right)={\displaystyle \frac{w^T{E}_{k}w}{w^{T}w}}\in \left[{\alpha }_E^{\left(k\right)},{\beta }_E^{\left(k\right)}\right]\equiv {\mathcal{I}}_{E_k},\hfill \\ \hfill & \hfill & i=0,\dots ,N\hfill \\ {\alpha }_R^{\left(k\right)}\equiv {\lambda }_{min}\left(R_k{R}_k^T\right),\hfill & {\beta }_R^{\left(k\right)}\equiv {\lambda }_{max}\left(R_k{R}_k^T\right),\hfill & {\gamma }_R^{\left(k\right)}\left(w\right)={\displaystyle \frac{w^T{R}_k{R}_k^{T}w}{w^{T}w}}\in \left[{\alpha }_R^{\left(k\right)},{\beta }_R^{\left(k\right)}\right]\equiv {\mathcal{I}}_{R_k},\hfill \\ & \hfill & i=1,\dots ,N.\hfill \end{array}$$

(5)

A third set of indicators, ${\gamma }_S^{\left(k\right)}$, linked to the previous ones as

$${\gamma }_S^{\left(k\right)}\left(w\right)={\gamma }_E^{\left(k\right)}\left(w\right)+{\gamma }_R^{\left(k\right)}\left(w\right),\forall w\in {\mathbb{R}}^{n_k},$$

(6)

describes the field-of-values of each preconditioned Schur complement ${\overline{S}}_k$:

$${\alpha }_S^{\left(k\right)}\equiv {\lambda }_{min}\left({\overline{S}}_k\right),{\beta }_S^{\left(k\right)}\equiv {\lambda }_{max}\left({\overline{S}}_k\right),{\gamma }_S^{\left(k\right)}\left(w\right)={\displaystyle \frac{w^T\overline{S}w}{w^{T}w}}\in \left[{\alpha }_S^{\left(k\right)},{\beta }_S^{\left(k\right)}\right]\equiv {\mathcal{I}}_{S_k},i=1,\dots ,N.$$

(7)

If the Schur complements are exactly inverted, then $E_0=I$ and $R_k{R}_k^{\top }+E_k=I$. In such a case, denoting ${\gamma }_R^{\left(0\right)}=0$, we have ${\gamma }_S^{\left(k\right)}\left(w\right)={\gamma }_E^{\left(k\right)}\left(w\right)+{\gamma }_R^{\left(k\right)}\left(w\right)\equiv 1,j=0,\dots ,N,\forall w$, and the eigenvalues of (3) are either $-1$ or 1, with multiplicity, respectively, $n_1+n_3+\dots$ and $n_0+n_2+\dots$ (see [1], Theorem 2.1). In the following, we will often remove the argument w from the previously defined indicators. We define the two vectors of parameters:

$${\gamma }_E=\left[{\gamma }_E^{\left(0\right)},\dots ,{\gamma }_E^{\left(N\right)}\right],{\gamma }_R=\left[{\gamma }_R^{\left(1\right)},\dots ,{\gamma }_R^{\left(N\right)}\right].$$

We will now make some very mild assumptions on two of these indicators:

The value 1 is strictly included in both ${\mathcal{I}}_{E_0}$ and ${\mathcal{I}}_{S_1}$, namely

1. 

${\alpha }_E^{\left(0\right)}<1<{\beta }_E^{\left(0\right)}$,

2. 

${\alpha }_S^{\left(1\right)}<1<{\beta }_S^{\left(1\right)}$.

These assumptions are very commonly satisfied in practice, meaning that 1 is in both the spectra of the preconditioned $(1,1)$ block and of the preconditioned Schur complement ${\overline{S}}_1$.

3. Characterization of the Eigenvalues of the Preconditioned Matrix

We now recursively define a sequence of polynomials, with ${\gamma }_E$ and ${\gamma }_R$ as parameters, whose roots are strictly related with eigenvalues of ${\mathcal{P}}^{-1}\mathcal{A}$.

Definition 1.

$$\begin{array}{ccc}\hfill U_0(\lambda ,{\gamma }_R,{\gamma }_E)& =& 1,\hfill \\ \hfill U_1(\lambda ,{\gamma }_R,{\gamma }_E)& =& \lambda -{\gamma }_E^{\left(0\right)},\hfill \\ \hfill U_{k+1}(\lambda ,{\gamma }_R,{\gamma }_E)& =& (\lambda (1+{\gamma }_R^{\left(k\right)})+{(-1)}^{k+1}{\gamma }_E^{\left(k\right)})U_k(\lambda ,{\gamma }_R,{\gamma }_E)\hfill \\ \hfill & & -{\gamma }_R^{\left(k\right)}{(\lambda +{(-1)}^k)}^2{U}_{k-1}(\lambda ,{\gamma }_R,{\gamma }_E),k\ge 1.\hfill \end{array}$$

(8)

Then a sequence of matrix valued function is also defined by recurrence:

Definition 2.

$$\begin{array}{cc}\hfill Y_1\left(\lambda \right)& =\lambda I-E_0,\hfill \\ \hfill Y_{k+1}\left(\lambda \right)& ={(-1)}^{k+1}{E}_k+\lambda \left(I+R_k{R}_k^{\top }\right)-{(\lambda +{(-1)}^k)}^2{R}_k{Y}_k{\left(\lambda \right)}^{-1}{R}_k^{\top },k\ge 1,\hfill \\ & \mathit{and}\lambda \mathrm{s}.\mathrm{t}.0\notin \sigma \left(Y_k\left(\lambda \right)\right).\hfill \end{array}$$

Notation. We use the notation ${\mathcal{I}}_k$ to denote the union of intervals bounding the roots of polynomials of the form $U_k$ over the valid range of ${\gamma }_E,{\gamma }_R$.

We recall a technical lemma, based on an idea in [24], whose proof can be found in [2].

Lemma 1.

Let Y be a symmetric matrix valued function defined in $F\subset \mathbb{R}$ and

$$0\notin [min\{\sigma \left(Y\right(\xi \left)\right)\},max\{\sigma \left(Y\right(\xi \left)\right)\left\}\right]\mathit{for} \mathit{all}\xi \in F.$$

Then, for arbitrary $s\ne \mathbf{0}$, there exists a vector $v\ne \mathbf{0}$ such that

$${\displaystyle \frac{{\mathit{s}}^{\top }Y{\left(\xi \right)}^{-1}\mathit{s}}{{\mathit{s}}^{\top }\mathit{s}}}={\displaystyle \frac{1}{{\gamma }_Z}}\mathit{with}{\gamma }_Z={\displaystyle \frac{v^{\top }Y\left(\xi \right)v}{v^{\top }v}}.$$

The next lemma, which will be used in the proof of the subsequent Theorem 1, links together the two sequences previously defined.

Lemma 2.

For every $\mathit{u}\ne 0$, there is a choice of γ for which

$${\displaystyle \frac{{\mathit{u}}^{\top }{Y}_{k+1}\left(\lambda \right)\mathit{u}}{{\mathit{u}}^{\top }\mathit{u}}}={\displaystyle \frac{U_{k+1}\left(\lambda \right)}{U_k\left(\lambda \right)}}\mathit{for} \mathit{all}\lambda \notin \bigcup _{j=1}^k{\mathcal{I}}_j.$$

Proof. 

This is shown by induction. We first define

$${\eta }_k\left(\lambda \right)={(-1)}^{k+1}{\gamma }_E^{\left(k\right)}+\lambda (1+{\gamma }_R^{\left(k\right)}),{\overline{\mu }}_k\left(\lambda \right)={(\lambda +{(-1)}^k)}^2.$$

(9)

For $k=0$ we have ${\displaystyle \frac{{\mathit{u}}^{\top }{Y}_1\left(\lambda \right)\mathit{u}}{{\mathit{u}}^{\top }\mathit{u}}=\lambda -{\gamma }_E^{\left(0\right)}=\frac{U_1\left(\lambda \right)}{U_0\left(\lambda \right)}}$ for all $\lambda \in \mathbb{R}$. If $k\ge 1$, the condition $\lambda \notin {\mathcal{I}}_k$ , together with the inductive hypothesis ${\displaystyle \frac{{\mathit{u}}^{\top }{Y}_k\left(\lambda \right)\mathit{u}}{{\mathit{u}}^{\top }\mathit{u}}=\frac{U_k\left(\lambda \right)}{U_{k-1}\left(\lambda \right)}}$, implies invertibility of $Y_k\left(\lambda \right)$. Moreover, this is equivalent to the condition $0\notin [min\{\sigma \left(Y\right(\xi \left)\right)\},max\{\sigma \left(Y\right(\xi \left)\right)\left\}\right]$ that guarantees the applicability of Lemma 1. Therefore, we can write

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mathit{u}}^{\top }{Y}_{k+1}\left(\lambda \right)\mathit{u}}{{\mathit{u}}^{\top }\mathit{u}}}& ={\eta }_k\left(\lambda \right)-{\overline{\mu }}_k\left(\lambda \right){\displaystyle \frac{{\mathit{u}}^{\top }{R}_k{Y}_k{\left(\lambda \right)}^{-1}{R}_k^{\top }\mathit{u}}{{\mathit{u}}^{\top }\mathit{u}}}\underset{w=R_k^{\top }\mathit{u}}{=}{\eta }_k\left(\lambda \right)-{\overline{\mu }}_k\left(\lambda \right){\displaystyle \frac{w^{\top }{Y}_k{\left(\lambda \right)}^{-1}w}{w^{\top }w}}{\gamma }_R^{\left(k\right)}\hfill \end{array}$$

(10)

We then apply Lemma 1 and the inductive hypothesis to write

$${\displaystyle \frac{w^{\top }{Y}_k{\left(\lambda \right)}^{-1}w}{w^{\top }w}}={\displaystyle \frac{U_{k-1}\left(\lambda \right)}{U_k\left(\lambda \right)}}.$$

Substituting into (10) and using relation (8) we obtain

$${\eta }_k\left(\lambda \right)-{\overline{\mu }}_k\left(\lambda \right){\displaystyle \frac{w^{\top }{Y}_k{\left(\lambda \right)}^{-1}w}{w^{\top }w}}{\gamma }_R^{\left(k\right)}={\displaystyle \frac{{\eta }_k\left(\lambda \right)U_k\left(\lambda \right)-{\overline{\mu }}_k\left(\lambda \right){\gamma }_R^{\left(k\right)}{U}_{k-1}\left(\lambda \right)}{U_k\left(\lambda \right)}}={\displaystyle \frac{U_{k+1}\left(\lambda \right)}{U_k\left(\lambda \right)}}.$$

   □

We now state the main results of this Section:

Theorem 1.

The eigenvalues of ${\mathcal{P}}^{-1}\mathcal{A}$ are located in ${\displaystyle \mathcal{I}=\bigcup _{k=1}^{N+1}{\mathcal{I}}_k}$.

Proof. 

The proof is carried out through induction on k, that is for every $k\le N+1$ either

$$\left(\mathrm{i}\right)\lambda \in {\mathcal{I}}_k\mathrm{or}\left(\mathrm{ii}\right){\mathit{u}}_k=(1+(-1^k)\lambda )Y_k{\left(\lambda \right)}^{-1}{R}_k^{\top }{\mathit{u}}_{k+1},$$

and for $k=N+1$ only condition (i) can hold. Let $\mathit{u}={\left[{\mathit{u}}_1^{\top },\dots ,{\mathit{u}}_{k+1}^{\top }\right]}^{\top }$ be an eigenvector of (4). Assume that $\lambda \notin {\mathcal{I}}_{E_0}$ (for $k=0$), then $Y_1\left(\lambda \right)$ is invertible. From the first equation of (4) we obtain

$$(\lambda I-E_0){\mathit{u}}_1=(1-\lambda )R_1^{\top }{\mathit{u}}_2\Rightarrow Y_1\left(\lambda \right){\mathit{u}}_1=(1-\lambda )R_1^{\top }{\mathit{u}}_2,$$

(11)

whereupon inserting (11) into the second row of (4) yields

$$\underset{Y_2\left(\lambda \right)}{\underbrace{\left(E_1+\lambda (I+R_1{R}_1^{\top })-{(\lambda -1)}^2{R}_1{(\lambda I-E_0)}^{-1}{R}_1^{\top }\right)}}{\mathit{u}}_2=(1+\lambda )R_2^{\top }{\mathit{u}}_3.$$

(12)

Pre-multiplying the left hand side of (12) by ${\displaystyle \frac{{\mathit{u}}_2^{\top }}{{\mathit{u}}_2^{\top }{\mathit{u}}_2}}$, we obtain that

$${\displaystyle \frac{{\mathit{u}}_2^{\top }{Y}_2\left(\lambda \right){\mathit{u}}_2}{{\mathit{u}}_2^{\top }{\mathit{u}}_2}}={\displaystyle \frac{U_2\left(\lambda \right)}{U_1\left(\lambda \right)}}$$

If $\lambda$ is a zero of $U_2$ then $\lambda \in {\mathcal{I}}_2$, Otherwise $Y_2\left(\lambda \right)$ is invertible and definite, and we can write ${\mathit{u}}_2=(1+\lambda )Y_2{\left(\lambda \right)}^{-1}{\mathit{u}}_3$. Assume now the inductive hypothesis holds for $k-1$. If $\lambda \notin {\mathcal{I}}_{k-1}$, then $Y_{k-1}\left(\lambda \right)$ is definite and invertible. We can write

$$\begin{array}{ccc}\hfill Y_k\left(\lambda \right){\mathit{u}}_k& \equiv & \left({(-1)}^k{E}_{k-1}+\lambda (I+R_{k-1}{R}_{k-1}^{\top })-{(\lambda +{(-1)}^{k-1})}^2{R}_{k-1}{Y}_{k-1}^{-1}{R}_{k-1}^{\top }\right){\mathit{u}}_k\hfill \\ & =& (1+{(-1)}^k\lambda )R_k^{\top }{\mathit{u}}_{k+1}.\hfill \end{array}$$

Since

$${\displaystyle \frac{{\mathit{u}}_k^{\top }{Y}_k\left(\lambda \right){\mathit{u}}_k}{{\mathit{u}}_k^{\top }{\mathit{u}}_k}}={\displaystyle \frac{U_k\left(\lambda \right)}{U_{k-1}\left(\lambda \right)}},$$

then either $\lambda$ is a zero of $U_k\left(\lambda \right)$ and hence $\lambda \in {\mathcal{I}}_k$ or $Y_k\left(\lambda \right)$ is definite, in such a case we can write ${\mathit{u}}_k=(1+{(-1)}^k\lambda )Y_k{\left(\lambda \right)}^{-1}{\mathit{u}}_{k+1}$. The induction process ends for $k=N+1$. In this case we have that

$$Y_{N+1}\left(\lambda \right){\mathit{u}}_{N+1}=0,$$

and the condition $\lambda \in {\mathcal{I}}_{N+1}$ must hold, noticing that ${\mathit{u}}_{N+1}=0$ would imply that also ${\mathit{u}}_N=\dots ={\mathit{u}}_1=0$ contradicting the definition of an eigenvector.    □

4. Bounds on the Roots of $\left\{{\mathit{U}}_{\mathit{k}}\right\}$

In this Section we will first characterize the set $\mathcal{I}=\bigcup {\mathcal{I}}_k$. Then we will specify the dependence of the zeros $U_k\left(\lambda \right)$ on the parameters ${\gamma }_E,{\gamma }_R$.

Proposition 1.

Since the polynomials of the sequence (8) are monic

$$\underset{\lambda \to +\infty }{lim}{U}_k\left(\lambda \right)=+\infty \underset{\lambda \to -\infty }{lim}{U}_k\left(\lambda \right)={(-1)}^k·\infty$$

Lemma 3.

Given the sequence of polynomials (8), for all $k\ge 1$

$$\mathrm{sgn}\left(U_{2k-1}\left(0\right)\right)=\mathrm{sgn}\left(U_{2k}\left(0\right)\right)={(-1)}^k.$$

In other words, the signs of $U_k\left(0\right)$ follow the repeating pattern $--,++,--,++,\dots$ starting from $k=1$. Equivalently, ${\displaystyle sgn\left(U_k\left(0\right)\right)={(-1)}^{\lceil \frac{k}{2}\rceil }}$.

Proof. 

We prove this by induction of k. The claim holds for $k=1$:

$$U_1\left(0\right)=-{\gamma }_E^{\left(0\right)}<0,U_2\left(0\right)=-{\gamma }_E^{\left(1\right)}{\gamma }_E^{\left(0\right)}-{\gamma }_R^{\left(1\right)}<0.$$

Assume that the claim holds for some $k\ge 1$, that is the sign of $U_{2k-1}\left(0\right)$ and $U_{2k}\left(0\right)$ is $s={(-1)}^k$, then

$$\begin{array}{ccc}\hfill \mathrm{sgn}\left(U_{2k+1}\left(0\right)\right)& =& {(-1)}^{2k+1}\mathrm{sgn}\left(U_{2k}\left(0\right)\right)-\mathrm{sgn}\left(U_{2k-1}\left(0\right)\right)=-s={(-1)}^{k+1}.\hfill \\ \hfill \mathrm{sgn}\left(U_{2k+2}\left(0\right)\right)& =& {(-1)}^{2k+2}\mathrm{sgn}\left(U_{2k+1}\left(0\right)\right)-\mathrm{sgn}\left(U_{2k}\left(0\right)\right)=-s={(-1)}^{k+1}.\hfill \end{array}$$

   □

We now show that, based on Assumptions 1, both $-1$ and 1 are strictly included in $\mathcal{I}$, in fact ${\alpha }_E^{\left(0\right)}<1<{\beta }_E^{\left(0\right)}$ implies 1 is strictly included in ${\mathcal{I}}_1$, and ${\alpha }_S^{\left(1\right)}<1<{\beta }_S^{\left(1\right)}$ implies that $-1$ is strictly included in ${\mathcal{I}}_2.$ In fact

$$U_2(-1;{\gamma }_E^{\left(0\right)},{\gamma }_E^{\left(1\right)},{\gamma }_R^{\left(1\right)})=1+{\gamma }_E^{\left(0\right)}-3{\gamma }_R^{\left(1\right)}+{\gamma }_E^{\left(0\right)}({\gamma }_R^{\left(1\right)}-{\gamma }_E^{\left(1\right)})-{\gamma }_E^{\left(1\right)}.$$

Using ${\gamma }_E^{\left(0\right)}=1,$ we have $U_2(-1;1,{\gamma }_E^{\left(1\right)},{\gamma }_R^{\left(1\right)})=2(1-{\gamma }_E^{\left(1\right)}-{\gamma }_R^{\left(1\right)})=2(1-{\gamma }_S^{\left(1\right)})\equiv U_2(-1;1,{\gamma }_S^{\left(1\right)})$. That $-1$ strictly belongs to ${\mathcal{I}}_2$ comes from $U_2(-1;1,{\alpha }_S^{\left(1\right)})>0$ and $U_2(-1;1,{\beta }_S^{\left(1\right)})<0$.

We now state a results regarding the extremal roots of the polynomials $\left\{U_k\right\}$.

Theorem 2.

Assume that the polynomial $U_k\left(\lambda \right)$ has $s_k$ distinct roots, and let us denote as ${\xi }_1^{\left(k\right)}<{\xi }_2^{\left(k\right)}<\dots {\xi }_{s_k}^{\left(k\right)}$ the roots of $U_k\left(\lambda \right)$. If ${\xi }_1^{\left(k\right)}<-1$ and ${\xi }_k^{\left(k\right)}>1,\forall k$, then the extremal roots of the polynomials $U_k$ satisfy

$${\xi }_1^{\left(k\right)}<{\xi }_1^{(k-1)}<\dots <{\xi }_1^{\left(2\right)},\mathit{and}{\xi }_1^{\left(1\right)}<\dots <{\xi }_{s_{k-1}}^{(k-1)}<{\xi }_{s_k}^{\left(k\right)}.$$

Proof. 

We prove the claim by induction, starting from $l=2$. Consider first the positive roots. The basis of the induction is

$$1<{\xi }_1^{\left(1\right)}<{\xi }_2^{\left(2\right)},$$

which is readily proved by taking ${\xi }_1^{\left(1\right)}={\gamma }_E^{\left(0\right)}>1$ and observing that $U_2\left({\gamma }_E^{\left(0\right)}\right)<0.$ Assume now that the claim holds for $l-1$, that is, ${\xi }_{s_{l-1}}^{(l-1)}<{\xi }_{s_l}^{\left(l\right)}$. This implies that $U_{l-1}\left({\xi }_{s_l}^{\left(l\right)}\right)>0$. Then, from the recursion (8)

$$U_{l+1}\left({\xi }_{s_l}^{\left(l\right)}\right)=-{\gamma }_R{({\xi }_{s_l}^{\left(l\right)}+{(-1)}^l)}^2{U}_{l-1}\left({\xi }_{s_l}^{\left(l\right)}\right)<0,$$

which implies that ${\xi }_{s_l}^{\left(l\right)}<{\xi }_{s_{l+1}}^{(l+1)}$.

Consider now the negative roots. If ${\gamma }_E^{\left(0\right)}<1$ and ${\gamma }_S^{\left(1\right)}>1$, direct computation shows that $U_2(-1)<0,$ providing ${\xi }_1^{\left(2\right)}<-1.$ Moreover, $U_3\left({\xi }_1^{\left(2\right)}\right)=-{\gamma }_R^{\left(2\right)}{({\xi }_1^{\left(2\right)}+1)}^2({\xi }_1^{\left(2\right)}-{\gamma }_E^{\left(0\right)})>0$, which implies ${\xi }_1^{\left(3\right)}<{\xi }_1^{\left(2\right)}$. Assume now that the claim holds for $l-1$, that is ${\xi }_1^{(l-1)}>{\xi }_1^{\left(l\right)}$ which implies that $\mathrm{sgn}\left(U_{l-1}\left({\xi }_1^{\left(l\right)}\right)\right)={(-1)}^{l-1}$. Then, from the recursion (8)

$$\mathrm{sgn}\left(U_{l+1}\left({\xi }_1^{\left(l\right)}\right)\right)=-\mathrm{sgn}\left(U_{l-1}\left({\xi }_1^{\left(l\right)}\right))\right)=-{(-1)}^{l-1}={(-1)}^l,$$

which finally shows that ${\xi }_1^{(l+1)}<{\xi }_1^{\left(l\right)}$, as desired.    □

The results of this section allow us to characterize the set $\mathcal{I}$ in Theorem 1 like

$$\mathcal{I}=[{\lambda }_{-}^{LB},{\lambda }_{-}^{UB}]\cup [{\lambda }_{+}^{LB},{\lambda }_{+}^{UB}],$$

where

$$\begin{array}{ccc}\hfill {\lambda }_{-}^{LB}& =& {\xi }_1^{(N+1)}\hfill \\ \hfill {\lambda }_{-}^{UB}& =& \underset{j}{max}\{max\{{\xi }_l^{\left(j\right)},\mathrm{s}.\mathrm{t}.{\xi }_l^{\left(j\right)}<0\}\}\hfill \\ \hfill {\lambda }_{+}^{LB}& =& \underset{j}{min}\{min\{{\xi }_l^{\left(j\right)},\mathrm{s}.\mathrm{t}.{\xi }_l^{\left(j\right)}>0\}\}\hfill \\ \hfill {\lambda }_{+}^{UB}& =& {\xi }_{s_{N+1}}^{(N+1)}.\hfill \end{array}$$

4.1. How Zeros of $\left\{U_k\right\}$ Move Depending on ${\gamma }_E^{\left(j\right)},{\gamma }_R^{\left(j\right)}$

Now the question is: for which values of the parameters ${\gamma }_E^{\left(i\right)},{\gamma }_R^{\left(i\right)}$ the extremal values of the roots is attained? The following results will allow us to state that the extremal values of the zeros of $U_{k+1}(\lambda ,{\gamma }_E,{\gamma }_R)$ are obtained at the extremal values of ${\gamma }_E$ and ${\gamma }_R$, namely $\{{\alpha }_E^{\left(i\right)},{\beta }_E^{\left(i\right)};{\alpha }_R^{\left(i\right)},{\beta }_R^{\left(i\right)}\}$.

Let $\xi$ an extremal zero of the polynomial $U_{k+1}$ (smallest/largest negative, smallest/largest positive) for a particular combination of the parameters. Then taking separately one of the parameters, ${\gamma }_{*}$ we can write ${\gamma }_{*}$ as a convex combination of its extremal values, ${\gamma }_{*}=\alpha {\gamma }_{*}^{min}+(1-\alpha ){\gamma }_{*}^{max}$ and write

$$\begin{array}{cc}\hfill U_{k+1}(\lambda ,{\gamma }_{*})& =s_1\left(\lambda \right)+{\gamma }_{*}{s}_2\left(\lambda \right)=s_1\left(\lambda \right)+(\alpha {\gamma }_{*}^{min}+(1-\alpha ){\gamma }_{*}^{max})s_2\left(\lambda \right)\hfill \\ \hfill & =\alpha (s_1\left(\lambda \right)+{\gamma }_{*}^{min}{s}_2\left(\lambda \right))+(1-\alpha )(s_1\left(\lambda \right)+{\gamma }_{*}^{max}{s}_2\left(\lambda \right))\hfill \\ \hfill & =\alpha U_{k+1}(\lambda ,{\gamma }_{*}^{min})+(1-\alpha )U_{k+1}(\lambda ,{\gamma }_{*}^{max}).\hfill \end{array}$$

Hence,

$$0=U_{k+1}(\xi ,{\gamma }_{*})=\alpha U_{k+1}(\xi ,{\gamma }_{*}^{min})+(1-\alpha )U_{k+1}(\xi ,{\gamma }_{*}^{max}),$$

with ${\rho }_1=U_{k+1}(\xi ,{\gamma }_{*}^{min})$ and ${\rho }_2=U_{k+1}(\xi ,{\gamma }_{*}^{max})$ therefore either taking opposite signs, or being both zero. In the first case it is clear that one between ${\rho }_1$ and ${\rho }_2$ improves the sought root, in particular we select $\gamma \in \{{\gamma }_{*}^{min},{\gamma }_{*}^{max}\}$ corresponding to the ${\rho }_{*}\in \{{\rho }_1,{\rho }_2\}$ such that

1.

(Smallest negative root): The sign of ${\rho }_{*}$ is opposite to the sign of $U_{k+1}$ at $-\infty$.

2.

(Largest negative root): The sign of ${\rho }_{*}$ is opposite to the sign of $U_{k+1}$ at 0.

3.

(Smallest positive root): The sign of ${\rho }_{*}$ is opposite to the sign of $U_{k+1}$ at 0.

4.

(Largest positive root): The sign of ${\rho }_{*}$ is opposite to the sign of $U_{k+1}$ at $+\infty$.

If, finally, both ${\rho }_1$ and ${\rho }_2$ are zero, for linearity, it means that the root is independent of this parameter and we can choose one of its extremal values.

We will now prove (in Theorem 3, and in Section 4.2) that we can predict whether ${\gamma }_E^{\left(i\right)}={\alpha }_E^{\left(i\right)}$ or ${\gamma }_E^{\left(i\right)}={\beta }_E^{\left(i\right)}$ in order to obtain the intervals bounding the roots of $U_k$.

Lemma 4.

Given a fixed k, there exist polynomials $W_l$ for any $0\le l\le k$ such that

$$U_{k+1}=U_j{W}_{k-j+1}-U_{j-1}{W}_{k-j}{\mu }_j,$$

with ${\mu }_j={\gamma }_R^{\left(j\right)}{(\lambda +{(-1)}^j)}^2.$

Proof. 

The statement holds for $j=k$, by setting $W_1={\eta }_k$ and $W_0=1$. Let the statement be true for $j\ge 2$, then,

$$\begin{array}{ccc}\hfill U_{k+1}& =& U_j{W}_{k-j+1}-U_{j-1}{W}_{k-j}{\mu }_j=\hfill \\ & =& W_{k-j+1}({\eta }_{j-1}{U}_{j-1}-{\mu }_{j-1}{U}_{j-2})-U_{j-1}{W}_{k-j}=\hfill \\ & =& \underset{W_{k-j+2}}{\underbrace{(W_{k-j+1}{\eta }_{j-1}-W_{k-j})}}{U}_{j-1}-W_{k-j+1}{\mu }_{j-1}{U}_{j-2}\hfill \end{array}$$

It is also clear by induction that $W_{k-j}$ depends only on the parameters with index greater than j.    □

We denote with $e_j\in \{{\alpha }_E^{\left(j\right)},{\beta }_E^{\left(j\right)}\}$ one choice of the extremal values of ${\gamma }_E^{\left(j\right)}$ and, analogously, $r_j\in \{{\alpha }_R^{\left(j\right)},{\beta }_R^{\left(j\right)}\}$; moreover we denote by $e_j^{*}$ and $r_j^{*}$ the other extremal value of the parameter. We now consider the polynomial $U_k$ with ${\gamma }_E^{\left(j\right)}\equiv e_j,{\gamma }_R^{\left(j\right)}\equiv r_j$ and also define ${\phantom{U}}^E{U}_k^j$ as the polynomial $U_k$ with ${\gamma }_E^{\left(j\right)}\equiv e_j^{*}$. Similarly we define ${\phantom{U}}^R{U}_k^j$ as the polynomial with ${\gamma }_R^{\left(j\right)}=r_j^{*}$. Then it is immediate from the recursion of $U_k$ that

$$\begin{array}{ccc}\hfill {\phantom{U}}^E{U}_{j+1}^j& =& U_{j+1}+{(-1)}^{j+1}(e_j^{*}-e_j)U_j\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\phantom{U}}^R{U}_{j+1}^j& =& U_{j+1}+[\lambda U_j-{(\lambda +{(-1)}^j)}^2{U}_{j-1}](r_j^{*}-r_j).\hfill \end{array}$$

(14)

Setting $z_j\left(\lambda \right)=\lambda U_j-{(\lambda +{(-1)}^j)}^2{U}_{j-1}$, the second equality can be rewritten simply as

$${\phantom{U}}^R{U}_{j+1}^j=U_{j+1}+z_j(r_j^{*}-r_j).$$

Lemma 5.

We have the following relations

$$\begin{array}{ccc}\hfill {\phantom{U}}^E{U}_{k+1}^j& =& U_{k+1}+{(-1)}^{j+1}{W}_{k-j}{U}_j(e_j^{*}-e_j)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\phantom{U}}^R{U}_{k+1}^j& =& U_{k+1}+W_{k-j}{z}_j(r_j^{*}-r_j).\hfill \end{array}$$

(16)

Proof. 

We write

$$U_{k+1}=U_{j+1}{W}_{k-j}-U_j{W}_{k-j-1}{\mu }_{j+1},$$

and notice that the only term in the right hand side that depends on ${\gamma }_E^{\left(j\right)}$ is $U_{j+1}$. Thus, using (13),

$$\begin{array}{cc}\hfill {\phantom{U}}^E{U}_{k+1}^j-U_{k+1}& =W_{k-j}({\phantom{U}}^E{U}_{j+1}^j-U_{j+1})={(-1)}^{j+1}{W}_{k-j}{U}_j(e_j^{*}-e_j).\hfill \end{array}$$

This proves the first formula. The second one is proved in the same way, using (14).    □

Theorem 3.

Let ξ be the largest (resp. smallest) positive or negative root of $U_{k+1}$ and let $e_j,r_j$ be the parameters that realize this root. Then ξ assumes the largest (resp. smallest) value when the ${\gamma }_E^{\left(j\right)}$ are alternatively the maximum or minimum.

Proof. 

It is enough to show that $(e_j^{*}-e_j)(e_{j-1}^{*}-e_{j-1})$ is negative. First we observe that if $W_{k-j}\left(\xi \right)$ is zero, then $\xi$ does not depend on ${\gamma }_E^{\left(j\right)}$. Thus we can assume $W_{k-j}\left(\xi \right),W_{k-j-1}\left(\xi \right)\ne 0$. We also assume, without loss of generality, that $U_{j+1}\left(\xi \right),U_j\left(\xi \right)\ne 0$. With these assumptions, and recalling that ${\mu }_j>0$, we have

$$\mathrm{sgn}\left(U_j{W}_{k-j}\right)=\mathrm{sgn}\left(U_{j-1}{W}_{k-j-1}\right).$$

Now we observe that ${\phantom{U}}^E{U}_{k+1}^j\left(\xi \right)$ must have a fixed sign as j varies to guarantee that the root $\xi$ is either maximal or minimal (this sign depends on the type of of the root – smallest/largest, positive/negative – and on k but cannot depend on j). We call this sign s. Now we have, using (15) and recalling that $U_{k+1}\left(\xi \right)=0$,

$$\begin{array}{cc}\hfill & \mathrm{sgn}(e_j^{*}-e_j)(e_{j-1}^{*}-e_{j-1})\hfill \\ & =\mathrm{sgn}({\phantom{U}}^E{U}_{k+1}^j\left(\xi \right)\left({\phantom{U}}^E{U}_{k+1}^{j-1}\left(\xi \right)\right){(-1)}^{k-j}{(-1)}^j\mathrm{sgn}\left(W_{k-j-1}{U}_j\right)\mathrm{sgn}\left(W_j{U}_{j-1}\right)\hfill \\ & =s^2{(-1)}^{2j+1}=-1.\hfill \end{array}$$

   □

4.2. Choice of ${\gamma }_E^{\left(k\right)}$

In the previous section we have shown that the ${\gamma }_E^{\left(k\right)}$ assume alternatively the maximum or minimum value. It is enough to fix the value of ${\gamma }_E^{\left(k\right)}$ to determine all the ${\gamma }_E^{\left(j\right)}$. We know, from (13) that $\mathrm{sign}(e_k^{*}-e_k)$ is equal to ${(-1)}^{k+1}\mathrm{sign}(U_k\left(\xi \right){\phantom{U}}^E{U}_{k+1}^k\left(\xi \right)$. We will consider four cases

• Let $\xi ={\lambda }_{+}^{UB}={\xi }_{s_{k+1}}^{(k+1)}$. We know that $\xi$ is larger than all the roots of $U_k$ and ${\phantom{U}}^E{U}_{k+1}^k$ and so $U_k\left(\xi \right),{\phantom{U}}^E{U}_{k+1}^k\left(\xi \right)>0$. Thus

  $$\mathrm{sign}(e_k^{*}-e_k)={(-1)}^{k+1}$$

  and so

  $$\mathrm{sign}(e_j^{*}-e_j)=(e_k^{*}-e_k){(-1)}^{k-j}={(-1)}^{k+1}{(-1)}^{k-j}={(-1)}^{j+1}.$$

  (17)
• Let $\xi ={\lambda }_{-}^{LB}={\xi }_1^{(k+1)}$. The root $\xi$ is negative and smaller than all the roots of $U_k$ and ${\phantom{U}}^E{U}_{k+1}^k$ and so

  $$\mathrm{sign}\left(U_k\left(\xi \right)\right)={(-1)}^k\mathrm{sign}\left({\phantom{U}}^E{U}_{k+1}^k\left(\xi \right)\right)={(-1)}^{k+1}.$$
  Thus

  $$\mathrm{sign}(e_k^{*}-e_k)={(-1)}^{k+1}{(-1)}^k{(-1)}^{k+1}={(-1)}^k$$

  and we have

  $$\mathrm{sign}(e_j^{*}-e_j)=\mathrm{sign}(e_k^{*}-e_k){(-1)}^{k-j}={(-1)}^k{(-1)}^{k-j}={(-1)}^j.$$

  (18)
• Let now $\xi ={\lambda }_{+}^{LB}$ and let m be the (smallest) index such that $U_{m+1}\left(\xi \right)=0$. Then the sign of $U_m\left(\xi \right)$ and ${\phantom{U}}^E{U}_{m+1}^m$ must be the sign that this polynomial assumes on 0. Thus

  $$\begin{array}{cc}\hfill \mathrm{sign}(e_m^{*}-e_m)& ={(-1)}^{m+1}\mathrm{sign}\left(U_m\left(0\right)\right)\mathrm{sign}\left({\phantom{U}}^E{U}_{m+1}^m\left(0\right)\right)\hfill \\ & ={(-1)}^{m+1}{(-1)}^{⌈{\displaystyle \frac{m+1}{2}}⌉}{(-1)}^{⌈{\displaystyle \frac{m}{2}}⌉}\hfill \\ & ={(-1)}^{m+1}{(-1)}^{m+1}=+\hfill \end{array}$$

  and

  $$\mathrm{sign}(e_j^{*}-e_j)=\mathrm{sign}(e_m^{*}-e_m){(-1)}^{m-j}={(-1)}^{m-j}.$$

  (19)
• Let finally $\xi ={\lambda }_{-}^{UB}$ and let m be the (smallest) index such that $U_{m+1}\left(\xi \right)=0$. Reasoning as above shows that

  $$\mathrm{sign}(e_m^{*}-e_m)=+\mathrm{and}\mathrm{sign}(e_j^{*}-e_j)={(-1)}^{m-j}.$$

  (20)

4.3. Sign of $r_j^{*}-r_j$ in Terms of the Sign of $z_j\left(\xi \right)$

Recall that we have

$$\begin{array}{cc}\hfill {\phantom{U}}^E{U}_{k+1}^j\left(\xi \right)& ={(-1)}^{j+1}{W}_{j+1}\left(\xi \right)U_j\left(\xi \right)(e_j^{*}-e_j),\hfill \\ \hfill {\phantom{U}}^R{U}_{k+1}^j\left(\xi \right)& =W_{j+1}\left(\xi \right)z_j\left(\xi \right)(r_j^{*}-r_j).\hfill \end{array}$$

If $W_{j+1}\left(\xi \right)=0$, then the root is independent of ${\gamma }_R^{\left(j\right)}$, so we can assume $W_{j+1}\left(\xi \right)\ne 0$. We also assume, without loss of generality, that $U_j\left(\xi \right)\ne 0$. Then we have

$$\begin{array}{cc}\hfill \mathrm{sign}(r_j^{*}-r_j)& =\mathrm{sign}\left(z_j\left(\xi \right){\phantom{U}}^R{U}_{k+1}^j\left(\xi \right)\right)\mathrm{sign}\left(W_{j+1}\right)\hfill \\ & =\mathrm{sign}\left(z_j\left(\xi \right)\right)\mathrm{sign}\left({\phantom{U}}^R{U}_{k+1}^j\left(\xi \right)\right){(-1)}^{j+1}\mathrm{sign}\left({\phantom{U}}^E{U}_{k+1}^j\left(\xi \right)U_j\left(\xi \right)\right)\mathrm{sign}(e_j^{*}-e_j)\hfill \\ & ={(-1)}^{j+1}\mathrm{sign}\left(U_j\left(\xi \right)\right)\mathrm{sign}(e_j^{*}-e_j)\mathrm{sign}\left(z_j\left(\xi \right)\right),\hfill \end{array}$$

where the last step follows from the fact that ${\phantom{U}}^E{U}_{k+1}^j\left(\xi \right)$ and ${\phantom{U}}^R{U}_{k+1}^j\left(\xi \right)$ must have the same sign. Again we distinguish four cases depending on the type of root we are considering.

• Let $\xi ={\lambda }_{+}^{UB}={\xi }_{s_{k+1}}^{(k+1)}$. Then $U_j\left(\xi \right)>0$ and $\mathrm{sign}(e_j^{*}-e_j)={(-1)}^{j+1}$ so

  $$\mathrm{sign}(r_j^{*}-r_j)={(-1)}^{j+1}{(-1)}^{j+1}\mathrm{sign}\left(z_j\left(\xi \right)\right)=\mathrm{sign}\left(z_j\left(\xi \right)\right).$$

  (21)
• Let $\xi ={\lambda }_{-}^{UB}={\xi }_1^{(k+1)}$. Then

  $$\mathrm{sign}\left(U_j\left(\xi \right)\right)={(-1)}^j\mathrm{and}\mathrm{sign}(e_j^{*}-e_j)={(-1)}^j$$

  so

  $$\begin{array}{cc}\hfill \mathrm{sign}(r_j^{*}-r_j)& ={(-1)}^{j+1}{(-1)}^j{(-1)}^j\mathrm{sign}\left(z_j\left(\xi \right)\right)\hfill \\ \hfill & ={(-1)}^{j+1}\mathrm{sign}\left(z_j\left(\xi \right)\right).\hfill \end{array}$$

  (22)
• Let now $\xi ={\lambda }_{+}^{LB}$ and $m\le k$ be the (smallest) index such that $U_{m+1}\left(\xi \right)=0$. Then

  $$\mathrm{sign}\left(U_j\left(\xi \right)\right)={(-1)}^{⌈{\displaystyle \frac{j}{2}}⌉}\mathrm{and}\mathrm{sign}(e_j^{*}-e_j)={(-1)}^{m-j}$$

  so

  $$\begin{array}{cc}\hfill \mathrm{sign}(r_j^{*}-r_j)& ={(-1)}^{j+1}{(-1)}^{⌈{\displaystyle \frac{j}{2}}⌉}{(-1)}^{m-j}\mathrm{sign}\left(z_j\left(\xi \right)\right)\hfill \\ & ={(-1)}^{m+1}{(-1)}^{⌈{\displaystyle \frac{j}{2}}⌉}\mathrm{sign}\left(z_j\left(\xi \right)\right).\hfill \end{array}$$
• Let finally $\xi ={\lambda }_{-}^{UB}$ and $m\le k$ be the (smallest) index such that $U_{m+1}\left(\xi \right)=0$. Reasoning as above shows that

  $$\mathrm{sign}(r_j^{*}-r_j)={(-1)}^{m+1}{(-1)}^{⌈{\displaystyle \frac{j}{2}}⌉}\mathrm{sign}\left(z_j\left(\xi \right)\right).$$

The sign of $z_j\left(\xi \right)$ cannot be computed in general.

4.4. Bounds for the Eigenvalues of the Preconditioned Matrix

A procedure to compute the intervals bounding the eigenvalues of ${\mathcal{P}}^{-1}\mathcal{A}$ with $N+1$ blocks can be described as follows:

• Upper bound for the positive eigenvalues. In view of Theorem 2, this is given by the largest root of the polynomial $U_{N+1}\left(\lambda \right)$. The sign of (17) is ${(-1)}^{j+1}$. When the sign of $e_j^{*}-e_j$ is positive, this means that the wrong indicator is larger than the correct one, and therefore we must set ${\gamma }_E^{\left(j\right)}={\alpha }_E^{\left(j\right)}$. Hence, as $\mathrm{sgn}(e_0^{*}-e_0)=-1$, we should start with ${\gamma }_E^{\left(0\right)}={\beta }_E^{\left(0\right)}$, and then alternating between the two extremal values: ${\beta }_E^{\left(0\right)},{\alpha }_E^{\left(1\right)},{\beta }_E^{\left(2\right)},\dots$.
• Lower bound for the negative eigenvalues. In view of Theorem 2, this is given by the smallest zero of the polynomial $U_{N+1}\left(\lambda \right)$. The sign of (18) is ${(-1)}^j$. This means that we should start with ${\gamma }_E^{\left(0\right)}={\alpha }_E^{\left(0\right)}$ (since $e_0^{*}-e_0>0$) and then alternating between the two extremal values: ${\alpha }_E^{\left(0\right)},{\beta }_E^{\left(1\right)},{\alpha }_E^{\left(2\right)},\dots$.
• Upper bounds for the negative eigenvalues. We consider evaluating the largest negative root ${\xi }_{-}^{\left(m\right)}$ of all polynomials $U_m,2\le m\le N+1$. The sign in (19) is ${(-1)}^{m+j}$ suggesting that we have to start with ${\alpha }_E^{\left(0\right)}$ if m is even and with ${\beta }_E^{\left(0\right)}$, if m is odd, and proceed as before. Finally we compute

  $${\lambda }_{-}^{UB}=\underset{m\le N+1}{max}\left\{{\xi }_{-}^{\left(m\right)}\right\}.$$
• Lower bounds for the positive eigenvalues. We consider evaluating the smallest positive root ${\xi }_{+}^{\left(m\right)}$ of all polynomials $U_m,1\le m\le N+1$. The sign in (20) is again ${(-1)}^{m+j}$ suggesting that we have to start with ${\alpha }_E^{\left(0\right)}$ if m is even and with ${\beta }_E^{\left(0\right)}$, if m is odd, and proceed as before. Finally we compute

  $${\lambda }_{+}^{LB}=\underset{m\le N+1}{min}\left\{{\xi }_{+}^{\left(m\right)}\right\}.$$

In each of the four cases, the choice of ${\gamma }_R^{\left(j\right)}\in [{\alpha }_R^{\left(j\right)},{\beta }_R^{\left(j\right)}]$ must be performed by taking all the combinations of these values. Some insights on how to select the extremal values of the ${\gamma }_R$ indicators, upon additional assumptions, are given in Lemma 6, whose proof is deferred to Appendix A, and in Theorem 4.

Lemma 6.

Let ρ a positive number and assume that, for every $j=0,\dots ,k$, ${\alpha }_E^{\left(j\right)}\le 2-\rho$ and ${\beta }_E^{\left(j\right)}\ge 2+\rho$. Assume further that $\left|\xi \right|>{\displaystyle \frac{1}{\rho }}$. Then $z_j\left({\xi }_{k+1}^{(k+1)}\right)<0$, and $sgn\left(z_j\left({\xi }_1^{(k+1)}\right)\right)={(-1)}^j$.

Under the assumptions of Lemma 6 we are now able to determine the monotonicity of the extremal roots of $U_{k+1}$ depending on ${\gamma }_R^{\left(j\right)}$.

Theorem 4.

Under the assumptions of Lemma 6, the largest (respectively, the smallest) root of $U_{k+1}$ assumes its largest (respectively, smallest) value, in combination with ${\gamma }_R^{\left(j\right)}={\beta }_R^{\left(j\right)},j=1,\dots ,k$.

Proof. 

Let $\xi$ be the largest positive root. Then the previous lemma shows that $\mathrm{sign}\left(z_j\left(\xi \right)\right)=(-1)$. But then (21) implies

$$\mathrm{sign}(r_j^{*}-r_j)=-1$$

and so $r_j={\beta }_R^{\left(j\right)}$. Now let $\xi$ be the smallest negative root. The previous lemma shows that $\mathrm{sign}\left(z_j\left(\xi \right)\right)={(-1)}^j$ but then (22) implies

$$\begin{array}{c}\hfill \mathrm{sign}(r_j^{*}-r_j)={(-1)}^{j+1}{(-1)}^j=-1\end{array}$$

and so again $r_j={\beta }_R^{\left(j\right)}$.    □

5. Numerical Results. Randomly Generated Matrices

We now undertake numerical tests to validate the theoretical bounds of Theorem 1, and subsequent characterizations. We determine the extremal eigenvalues of ${\mathcal{P}}^{-1}\mathcal{A}$ on randomly-generated linear systems. Specifically, we considered a simplified case with $E_i\equiv 0,i>0$ and run a number of different test cases , combining the values of the extremal eigenvalues (Table 1) of the symmetric positive definite matrices involved. In particular

• Case $N=2$. We have three parameters with 6 different endpoints of the corresponding intervals. Overall we run $3^6=729$ test cases.
• Case $N=3$. We have 4 parameters with 4 different endpoints of the corresponding intervals. Overall we run $4^4=256$ test cases.
• Case $N=4$. We have 5 parameters with 4 different endpoints of the corresponding intervals. Overall we run $5^4=1024$ test cases.

From Figure 1, Figure 2 and Figure 3 we notice that three out of four bounds perfectly capture the eigenvalues of the preconditioned matrix, while the upper bounds for negative eigenvalues are not as tight for $N=2,4$ and lower bounds for positive eigenvalues are less effective for $N=3$. All in all, our technique is able to predict the behavior of the MINRES iterative solver applied to the preconditioned saddle-point linear system.

5.1. Comparisons Against the Block Diagonal Preconditioner

The preconditioner we have studied so far, being SPD, allows the MINRES iteration for the iterative solution of the multiple saddle-point linear systems. Another, most known, preconditioner sharing the same properties is the block diagonal preconditioner ${\mathcal{P}}_D$, based on inexact Schur complements, defined in (1).

We compare the MINRES number of iterations in all the previously described test cases for $N=3$. In Figure 4 we plot the number of iterations with the PP preconditioner and the block diagonal preconditioner, together with the square root of $\kappa$, an indicator of the ill-conditioning of the overall saddle-point system, namely

$$\kappa ={\displaystyle \frac{{\beta }_E^{\left(0\right)}}{{\alpha }_E^{\left(0\right)}}}{\displaystyle \frac{{\beta }_R^{\left(1\right)}{\beta }_R^{\left(2\right)}{\beta }_R^{\left(3\right)}}{{\alpha }_R^{\left(1\right)}{\alpha }_R^{\left(2\right)}{\alpha }_R^{\left(3\right)}}}.$$

From the figure, we notice that the two preconditioners are both affected by the $\kappa$ indicator, with the PP preconditioner more effective than ${\mathcal{P}}_D$ for small to modest values of $\kappa$. When the Schur complements are instead not optimally approximated, the block diagonal preconditioner is to be preferred.

6. A Realistic Example: The 3-D Biot’s Consolidation Model

A numerical experiment, concerning the mixed form of Biot’s poroelasticity equations, is used to investigate the quality of the bounds for the eigenvalues of the preconditioned matrix and the effectiveness of the proposed triangular preconditioners on a realistic problem. For the PDEs and boundary conditions governing this model, we refer to the description in [19]. As a reference test case for the experiments, we considered the porous cantilever beam problem originally introduced in [25] and already used in [26]. The domain is the unit cube with no-flow boundary conditions along all sides, zero displacements along the left edge, and a uniform load applied on top. The material properties are described in ([19], Table 5,1). We considered the Mixed-Hybrid Finite Element discretization of the coupled Biot equations, which gives rise to a double saddle point linear system. The test case refers to a unitary cubic domain, uniformly discretized using a $h=0.05$ mesh size in all the spatial dimensions. The size of the block matrices and the overall number of nonzeros in the double saddle-point system are reported in Table 2.

6.1. Handling the Case $n_2>n_1$

The theory previously developed is based on the assumption $n_0\ge n_1\ge n_2\dots$, which is not verified in this problem, where $n_2=25200>n_1=8000$. The main consequence of this is that $R_2{R}_2^{\top }$ is singular and ${\gamma }_R^{\left(2\right)}=0$. This drawback can be circumvented by attacking the eigenvalue problem (3) in a different way. Let us write (3) for a double saddle-point linear system:

$$\begin{array}{ccc}\hfill (E_0-\lambda I){\mathit{u}}_1+(1-\lambda )R_1^{\top }{\mathit{u}}_2& =& 0,\hfill \\ \hfill (1-\lambda )R_1{\mathit{u}}_1+(-E_1-\lambda (I+R_1{R}_1^{\top })){\mathit{u}}_2+(1+\lambda )R_2^{\top }{\mathit{u}}_3& =& 0,\hfill \\ \hfill (1+\lambda )R_2{\mathit{u}}_2+(E_2-\lambda (I+R_2{R}_2^{\top })){\mathit{u}}_3& =& 0.\hfill \end{array}$$

(23)

As in the proof of Theorem 1, we have, for all $\lambda \notin [{\alpha }_E^{\left(0\right)},{\beta }_E^{\left(0\right)}]$

$$\left(E_1+\lambda (I+R_1{R}_1^{\top })-{(\lambda -1)}^2{R}_1{(\lambda I-E_0)}^{-1}{R}_1^{\top }\right){\mathit{u}}_2=(1+\lambda )R_2^{\top }{\mathit{u}}_3.$$

(24)

Assuming now that $\lambda \notin \left[{\displaystyle \frac{{\alpha }_E^{\left(2\right)}}{1+{\beta }_R^{\left(2\right)}}},{\displaystyle \frac{{\beta }_E^{\left(2\right)}}{1+{\alpha }_R^{\left(2\right)}}}\right]\equiv \overline{\mathcal{I}}$, which ensures that $E_2-\lambda (I+R_2{R}_2^{\top })$ is SPD; we obtain ${\mathit{u}}_3$ from the third of (23)

$${\mathit{u}}_3=-(\lambda +1){(E_2-\lambda (I+R_2{R}_2^{\top }))}^{-1}{R}_2{\mathit{u}}_2$$

Substituting this into (24) yields

$$\begin{array}{ccc}\hfill & & \left(E_1+\lambda (I+R_1{R}_1^{\top })-{(\lambda -1)}^2{R}_1{(\lambda I-E_0)}^{-1}{R}_1^{\top }\right.\hfill \\ & & +\left.{(1+\lambda )}^2{R}_2^{\top }{(E_2-\lambda (I+R_2{R}_2^{\top }))}^{-1}{R}_2\right){\mathit{u}}_2=0\hfill \end{array}$$

(25)

Premultiplying now by ${\mathit{u}}_2^{\top }$ and defining $w=R_1^{\top }{\mathit{u}}_2$ and $z=R_2{\mathit{u}}_2$, yields

$$\begin{array}{ccc}& & {(\lambda -1)}^2{\displaystyle \frac{w^{\top }{(E_0-\lambda I)}^{-1}w}{w^{\top }w}}{\gamma }_R^{\left(1\right)}+\lambda (1+{\gamma }_R^{\left(1\right)})+{\gamma }_E^{\left(1\right)}\hfill \\ & & +{(1+\lambda )}^2{\displaystyle \frac{z^{\top }{(E_2-\lambda (I+R_2^{\top }{R}_2))}^{-1}z}{z^{\top }z}}{\overline{\gamma }}_R^{\left(2\right)}=0.\hfill \end{array}$$

Before proceeding we notice that ${\overline{\gamma }}_R^{\left(2\right)}$, differently from definition (5), is the Rayleigh quotient of $R_2^{\top }{R}_2$. Due to the fact that $R_2$ has full rank we can bound ${\overline{\gamma }}_R^{\left(2\right)}$ as $0<{\alpha }_R^{\left(2\right)}\le {\overline{\gamma }}_R^{\left(2\right)}\le {\beta }_R^{\left(2\right)}.$ Hence, from now on, we will refer to ${\overline{\gamma }}_R^{\left(2\right)}$ as simply ${\gamma }_R^{\left(2\right)}$. Applying Lemma 1 to both $E_0-\lambda I$ and $E_2-\lambda (I+R_2^{\top }{R}_2)$, we rewrite the previous as

$$\begin{array}{ccc}\hfill 0& =& -{\gamma }_R^{\left(1\right)}{\displaystyle \frac{{(\lambda -1)}^2}{\lambda -{\gamma }_E^{\left(0\right)}}}+\lambda (1+{\gamma }_R^{\left(1\right)})+{\gamma }_E^{\left(1\right)}-{(1+\lambda )}^2{\gamma }_R^{\left(2\right)}{\displaystyle \frac{1}{\lambda (1+{\gamma }_R^{\left(2\right)})-{\gamma }_E^{\left(2\right)}}}\hfill \\ & =& {\displaystyle \frac{U_2\left(\lambda \right)}{U_1\left(\lambda \right)}}-{\displaystyle \frac{{(1+\lambda )}^2{\gamma }_R^{\left(2\right)}}{\lambda (1+{\gamma }_R^{\left(2\right)})-{\gamma }_E^{\left(2\right)}}}={\displaystyle \frac{U_3\left(\lambda \right)}{U_1\left(\lambda \right)(\lambda (1+{\gamma }_R^{\left(2\right)})-{\gamma }_E^{\left(2\right)})}},\hfill \end{array}$$

which tells that any eigenvalue of ${\mathcal{P}}^{-1}\mathcal{A}$ outside $[{\alpha }_E^{\left(0\right)},{\beta }_E^{\left(0\right)}]\cup \overline{\mathcal{I}}$ is a zero of $U_3\left(\lambda \right)$, meaning that any eigenvalue is in ${\mathcal{I}}_1\cup \overline{\mathcal{I}}\cup {\mathcal{I}}_3$.

6.2. Verification of the Bounds

To approximate the $(1,1)$ block, we employed the classical incomplete Cholesky factorization (IC) with fill-in based on a drop tolerance $\delta$. The approximations ${\widehat{S}}_1$ and ${\widehat{S}}_2$ of the Schur complements are as in [19], which prove completely scalable with the problem size. Even if the IC preconditioner does not scale with the discretization parameter, it is helpful to conduct a spectral analysis as a function of the drop tolerance.

In Table 3 we report the extremal values of the five indicators for this problem.

We computed the relevant eigenvalues of the preconditioned matrix ${\mathcal{P}}^{-1}\mathcal{A}$, using the Matlab function eigs with a function handle for the application of the preconditioned matrix (which we could not compute explicitly, due to its size and density) as well as the four endpoints of $\left\}\right\}{I}_3$, as predicted by the theory, and reported all of them in Table 4.

6.3. Improving the Upper Bounds for the Negative Eigenvalues

From Table 4 we see that the bounds perfectly capture the extremal eigenvalues, except for the upper bound for the negative eigenvalues, which is still loose. However, this can be improved by taking into account the indicator ${\gamma }_S^{\left(k\right)}\left(w\right),$ defined in (7). The endpoints of this indicator are, for this test case:

$${\mathcal{I}}_{S_1}=[0.3187,1.2408],{\mathcal{I}}_{S_2}=[0.9998,1.0002],$$

for every value of $\delta$. In order to use these indicators in the eigenvalue analysis, we perform a change of variable, using (6), in the polynomial recurrence which reads as, for $k\ge 1$,

$$U_{k+1}\left(\lambda \right)=(\lambda (1+{\gamma }_S^{\left(k\right)}-{\gamma }_E^{\left(k\right)})+{(-1)}^{k+1}{\gamma }_E^{\left(k\right)})U_k\left(\lambda \right)+({\gamma }_E^{\left(k\right)}-{\gamma }_S^{\left(k\right)}){(\lambda +{(-1)}^k)}^2{U}_{k-1}\left(\lambda \right),$$

(26)

for which all the results obtained in Section 3 still hold, as well as the observation at the beginning of Section 4.1, stating that the bounds for the roots of $U_{k+1}$ are taken at the endpoints of the parameters. On the contrary, in this case, we cannot select which of the extremal value of ${\gamma }_E^{\left(k\right)}$ and ${\gamma }_S^{\left(k\right)}$ provides the bound. However, in view of the low degree of the polynomial we can compute its roots for every combination of the endpoints of the indicators, and take the minimum/maximum values of the roots, depending on the bounds sought. To do this we also should consider the following

Remark 1.

The definition of ${\gamma }_S^{\left(k\right)}\left(w\right)$, constrains this indicator to be larger than the corresponding ${\gamma }_E^{\left(k\right)}\left(w\right)$. Hence in correspondence of ${\gamma }_E^{\left(k\right)}={\beta }_E^{\left(k\right)}$, the lower bound for ${\gamma }_S^{\left(k\right)}$ is $min\{{\beta }_E^{\left(k\right)},{\alpha }_S^{\left(k\right)}\}$.

Taking this into account, we compute the 3 roots of

$$U_3(\lambda ;{\gamma }_E^{\left(0\right)},{\gamma }_E^{\left(1\right)},{\gamma }_E^{\left(2\right)},{\gamma }_S^{\left(1\right)},{\gamma }_S^{\left(2\right)}),$$

for the $2^5=32$ combinations of the extremal values of the indicators, and call ${\xi }_{-}^{LB},{\xi }_{-}^{UB},{\xi }_{+}^{LB},{\xi }_{+}^{UB}$ the extremal values of the roots. In view of the discussion in Section 6.3, we finally compute as the eigenvalue bounds the following expressions:

$$\begin{array}{ccc}\hfill {\lambda }_{-}^{\mathrm{LB}}={\xi }_{-}^{LB},& & {\lambda }_{-}^{\mathrm{UB}}={\xi }_{-}^{UB}\hfill \\ \hfill {\lambda }_{+}^{\mathrm{LB}}=min\left\{{\alpha }_E^{\left(0\right)},{\displaystyle \frac{{\alpha }_E^{\left(2\right)}}{1+{\beta }_R^{\left(2\right)}}},{\xi }_{+}^{LB}\right\},& & {\lambda }_{+}^{\mathrm{UB}}=max\left\{{\beta }_E^{\left(0\right)},{\displaystyle \frac{{\beta }_E^{\left(2\right)}}{1+{\alpha }_R^{\left(2\right)}}},{\xi }_{+}^{UB}\right\}\hfill \end{array}$$

which we report in Table 5, together with the exact eigenvalues.

The adherence of the bounds to the extremal eigenvalues is now impressive.

6.4. Comparisons with the Block Diagonal Preconditioner for the Biot Problem

We finally compare the block diagonal with the PP preconditioners in terms of number of iterations, for different values of the parameter $\delta$, and report the obtained results in Table 6, showing that the preconditioner analyzed in this work represents a valid alternative to the block diagonal preconditioner, especially when the Schur complements are all well approximated. Figure 5 compare the convergence profile of MINRES with both preconditioners for $\delta ={10}^{-5},{10}^{-6}$.

A final remark is in order: In this study we have not been concerned with the cost of the application of the two preconditioners, at each MINRES iteration, which is crucial to evaluate the performance of both. This comparison is outside the scope of this paper, which is mainly concerned with eigenvalue distribution and number of iterations. However, even taking into account that the application of the PP preconditioner is from 1.5 to (slightly less than) 2 times more expensive than that of the block diagonal preconditioner, our results show that the PP preconditioner can be a viable alternative to accelerate the MINRES solution of multiple saddle-point linear systems.

7. Conclusions

We have developed an eigenvalue analysis of the preconditioned multiple saddle-point linear systems, with the SPD preconditioner (PP in short) proposed by Pearson and Potschka in [1]. This analysis is based on the (approximate) knowledge of the smallest and the largest eigenvalues of a number of SPD matrices, related to the preconditioned blocks of the original system. The obtained bounds reveal very close to the exact extremal eigenvalues of the preconditioned matrix, and therefore, this study represents a valuable tool to predict with some accuracy the number of MINRES iterations to solve the multiple saddle-point linear system. In a number of both synthetic and realistic test cases, the PP preconditioner reveals more efficient than the block diagonal preconditioner when the Schur complement matrices are well approximated.