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

HAL
open science

# On some structural properties of generalized Lyapunov eigenproblems and application to operator preconditioning

Valeria Simoncini, Daniele Toni

► To cite this version:

Valeria Simoncini, Daniele Toni. On some structural properties of generalized Lyapunov eigenproblems and application to operator preconditioning. 2023. hal-03942108

HAL Id: hal-03942108

https://hal.science/hal-03942108

Preprint submitted on 16 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# ON SOME STRUCTURAL PROPERTIES OF GENERALIZED LYAPUNOV EIGENPROBLEMS AND APPLICATION TO OPERATOR PRECONDITIONING*

VALERIA SIMONCINI†† AND DANIELE TONI‡

**Abstract.** We are interested in generalized matrix eigenvalue problems of the type $AX+XA = \lambda MXM$ and $AX+XA = \lambda(MX+XM)$ with $A$ and $M$ both symmetric and positive definite, and in their tensor counterparts. We collect several structural properties, some of which are known, together with some new spectral results. We also analyze in detail the case where the second problem stems from the discretization of linear elliptic partial differential equations by finite differences. In particular, we derive spectral properties that can be used in the numerical solution of the resulting algebraic linear system.

**Key words.** Matrix eigenproblem. Lyapunov equation. Structural properties. Eigenvalue distribution.

**AMS subject classifications.** 15A22, 15A24, 65F15, 65F45.

## 1. Introduction.
We consider the following matrix eigenvalue problems

$$ (1.1) \qquad AX + XA = \lambda MXM $$

and

$$ (1.2) \qquad AX + XA = \lambda(MX + XM), $$

with $A, M \in \mathbb{R}^{n \times n}$ symmetric and positive definite. Both equations can be viewed as generalized Lyapunov eigenvalue problems, where the right-hand side has one or two terms. With some abuse of notation, we will refer to $X$ as an “eigenvector matrix” or “eigenmatrix”.

Problem (1.2) occurs in different applications, such as in the detection of Hopf bifurcations [25],[26], and in the analysis of preconditioned solvers for discretized elliptic equations. More generally and with $A, M$ not necessarily symmetric, the two problems arise in different numerical strategies to characterize the spectral distance from the imaginary axis [9],[26] and classically studied in the context of multiparameter eigenproblems; see, e.g., [1, Ch. 7],[16]. The two problems share some spectral properties, however the right-hand side term provides a quite different setting in the two cases.

In the following we will collect some of the scarce results available in the literature, and provide new ones, associated with some specific settings. In several occasions we will also generalize our findings to the (three-dimensional, or order-3) tensor setting, namely to the equations,

$$ (1.3) \qquad (A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A)\boldsymbol{x} = \lambda(M \otimes M \otimes M)\boldsymbol{x}, $$

<footnote>*Version of January 11, 2023.</footnote>

<footnote>†IMATI-CNR, Via Ferrata 5/A, Pavia, Italia</footnote>

<footnote>‡Dipartimento di Matematica and (AM)², Alma Mater Studiorum - Università di Bologna, Piazza di Porta San Donato 5, 40126 Bologna, Italy. Emails: valeria.simoncini@unibo.it, daniele.toni2@studio.unibo.it</footnote>
and

$$ (1.4) \quad (A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A)\boldsymbol{x} = \lambda(M \otimes I_n \otimes I_n + I_n \otimes M \otimes I_n + I_n \otimes I_n \otimes M)\boldsymbol{x}, $$

where $I_n$ is the identity matrix of size $n$, and $\otimes$ stands for the Kronecker product, which for matrices $A \in \mathbb{R}^{n_A \times m_A}$, $A = (a_{ij})_{i=1,...,n_A, j=1,...,m_A}$, and $B \in \mathbb{R}^{n_B \times m_B}$ is defined as ([17])

$$ (1.5) \quad A \otimes B = \begin{bmatrix} a_{11}B & a_{12}B & \cdots & a_{1m_H}B \\ a_{21}B & a_{22}B & \cdots & a_{2m_H}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{n_H 1}B & a_{n_H 2}B & \cdots & a_{n_H m_H}B \end{bmatrix} \in \mathbb{R}^{n_A n_B \times m_A m_B}. $$

We refrain from further extending the analysis to $d$-dimensional tensors, as the three-dimensional case already provides a quite good indication of the type of properties obtained in the multiarray setting, and of the tools used for proving them.

The most exercised setting is probably the case corresponding to $M = I$ in (1.1), that is the following standard matrix Lyapunov eigenvalue problem

$$ (1.6) \qquad AX + XA = \lambda X. $$

Letting $(\theta_i, z_i)$, $i = 1, ..., n$ be the eigenpairs of $A$, then the eigenpairs of (1.6) are obtained as $\lambda = \theta_i + \theta_j$, $X = z_i z_j^T$ for all possible $i, j = 1, ..., n$ [17, Th. 4.4.5].

Analogously, the following tensor Lyapunov eigenvalue problem

$$ (1.7) \qquad (A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A)\boldsymbol{x} = \lambda\boldsymbol{x} $$

is equipped with eigenvalues $\lambda = \theta_i + \theta_j + \theta_l$, for all possible $i, j, l = 1, ..., n$, and corresponding eigenvectors $\boldsymbol{x} = z_i \otimes z_j \otimes z_l$; see, e.g., [21].

These spectral constructions do not straightforwardly carry over to the generalized case considered here, that it, the case where matrices arise on both sides of the equality. Nonetheless, in section 3 and section 4 we show that the eigenmatrices and eigentensors do maintain certain structures such as symmetry or low rank characterizations. These properties would not be captured in the vectorized formulation. In section 5 we further specialize the setting (1.2) to the case where $A$ and $M$ are related. This may occur when so-called operator preconditioning strategies are used to solve the linear system obtained by the finite difference discretization of linear elliptic selfadjoint differential problems in two or three space dimensions; see [10] and the references in section 5.

**1.1. Notation.** Capital roman letters denote matrices, usually of size $n$, while bold capital roman letters are used for matrices stemming from Kronecker products and sums, and bold small roman letters refer to vectors of corresponding size. Bold face is also used for three-dimensional (or three-modes) tensors $\boldsymbol{X} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$. We let $\mathbf{1}_n$ be the vector of all ones with $n$ components; the subscript will be omitted when clear from the context. We define $\mathbf{0}_n$ similarly. The matrix $A = \text{diag}(a)$ is a diagonal matrix with the components of the vector $a$ on its diagonal, and $A = \text{blkdiag}(B, C)$ is a block diagonal matrix having the square matrices $B, C$ as diagonal blocks.

We also recall that for $\boldsymbol{x} = \text{vec}(X)$ – where “vec” stacks all matrix columns one below the other – and conforming dimensions, it holds that $\boldsymbol{x}^T(I \otimes M)\boldsymbol{x} = \text{trace}(X^T M X)$ and $\boldsymbol{x}^T(M \otimes I)\boldsymbol{x} = \text{trace}(X M X^T)$.

When working with tensors we shall mainly work with their Kronecker formulation, while we refer, for instance, to [12, section 12.4] for the adopted notation.
**2. Spectral decomposition of linear tensor operators.** By using the Kron- ecker operator, the problems (1.1) and (1.2) can be written as vector generalized eigenproblems, namely

$$ (2.1) \qquad (A \otimes I + I \otimes A)x = \lambda(M \otimes M)x, $$

$$ (2.2) \qquad (A \otimes I + I \otimes A)x = \lambda(M \otimes I + I \otimes M)x, $$

with $x = \text{vec}(X)$. Note that the tensor generalizations (1.3) and (1.4) have precisely this structure. Within the vector framework, all these equations have the form $Ax = Mx\lambda$ with $A, M \in \mathbb{R}^{n^d \times n^d}$, $d = 2, 3$, symmetric and positive definite. Hence, a full set of eigenvectors $\{x_1, \dots, x_{n^d}\}$ can be determined, satisfying $x_i^T M x_j = 0$ for $i \neq j$, associated with corresponding real and positive eigenvalues, thus completely characterizing the eigenvectors, from the vector space point of view. On the other hand, in the vector formulations (2.1) and (2.2), for instance, structural properties such as symmetry or low rank of the eigenmatrices $X$ are lost. These properties may have important implications in applications, but also on the obtainable accuracy of computational methods [9],[25]. Moreover, we will see in the following that taking into account the structure allows one to naturally recover certain optimality properties of acceleration procedures associated with the discretization of partial differential equations (PDEs).

Finally, we observe that problem (2.2) may be viewed as the generalized eigen- problem associated with (2.1), thus possibly representing the major similarity between the two problems. Indeed, multiplying by $(M^{\frac{1}{2}} \otimes M^{\frac{1}{2}})^{-1}$ from the left both equations, and letting $\hat{A} = M^{-\frac{1}{2}}A M^{-\frac{1}{2}}$ we obtain

$$ (2.3) \qquad (\hat{A} \otimes M^{-1} + M^{-1} \otimes \hat{A})y = \lambda y, $$

$$ (2.4) \qquad (\hat{A} \otimes M^{-1} + M^{-1} \otimes \hat{A})y = \lambda(I \otimes M^{-1} + M^{-1} \otimes I)y, $$

with $y = (M^{\frac{1}{2}} \otimes M^{\frac{1}{2}})x$. This connection and in particular the form in (2.4) will lead to new developments in section 5.

**3. On the Kronecker pair** $(A \otimes I + I \otimes A, M \otimes M)$. With appropriate scaling, the generalized eigenvectors in (2.1) satisfy:

$$ (3.1) \qquad x_j^T (M \otimes M) x_i = \delta_{i,j}, \qquad x_j^T (A \otimes I + I \otimes A) x_i = \lambda_i \delta_{i,j}, $$

where $\delta_{i,j}$ is the Kronecker delta. Some of the results below are matrix rephrasings of these vector relations that highlight structural properties. In the following, we use $\langle X, Y \rangle = \text{trace}(X^T Y)$ for matrices $X, Y$ having the same dimensions.

**PROPOSITION 3.1.** The eigenmatrices $\{X_i\}_{i=1,...,n^2}$ satisfy

i) Symmetry: The matrix indexes can be ordered so that $X_i$ is symmetric for $i = 1, \dots, (n^2 + n)/2$, while the remaining $(n^2 - n)/2$ matrices are skew-symmetric;

ii) Orthogonality: $\langle (M^{\frac{1}{2}} X_i M^{\frac{1}{2}}), (M^{\frac{1}{2}} X_j M^{\frac{1}{2}}) \rangle = \delta_{ij}$ and $\lambda_i \delta_{ij} = 2 \langle X_i, AX_j \rangle$, $i, j = 1, \dots, n^2$.

*Proof.* We observe that if $X_i$ satisfies (1.1) then also $X_i^T$ satisfies the equa- tion for the same $\lambda$: it suffices to transpose the whole equation and use the sym- metry of $A$ and $M$. If it holds that $X_i^T = \pm X_i$, then $X_i$ is either symmetric or skew-symmetric. Since a symmetric matrix has $(n^2 + n)/2$ degrees of freedom, and skew-symmetric ones have $(n^2 - n)/2$ degrees of freedom, the result in (i) fol- lows. If $X_i^T \neq \pm X_i$, then we let $V_1 = (X_i + X_i^T)$ symmetric and $V_2 = (X_i - X_i^T)$
skew-symmetric. Then range([vec(V₁), vec(V₂)]) = range([vec(Xᵢ), vec(Xᵢᵀ)], that is vec(V₁), vec(V₂) are again eigenvectors and they span the same (vector) eigenspace as vec(Xᵢ), vec(Xᵢᵀ). The result in (i) thus also follows.

The results in (ii) directly follow from the matricization of the orthogonality properties of the eigenvectors xᵢ's. □

The proof of the previous result shows that if Xᵢ is neither symmetric nor skew-symmetric, then the corresponding eigenvalue must have multiplicity (at least) two, and that the symmetric, skew-symmetric matrices (Xᵢ + Xᵢᵀ), (Xᵢ - Xᵢᵀ) span the same eigenspace. In particular, simple eigenvalues only admit symmetric or skew-symmetric eigenmatrices.

As a side result highlighting the role of the eigenmatrix structure, we report a characterization for the solution of a related linear matrix equation.

**PROPOSITION 3.2.** Let {λᵢ}, {Xᵢ}ᵢ=1,...,n² be the eigenvalues and eigenmatrices associated with the pair (A⊗I+I⊗A, M⊗M). The solution to AY+YA+MYM = C with C symmetric can be written as

$$
Y = \sum_{i=1}^{n^2} \alpha_i X_i, \quad \text{with} \quad \alpha_i = \frac{\operatorname{trace}(X_i^T C)}{1 + \lambda_i}.
$$

Proof. Since {$\mathbf{x}_i$} with $\mathbf{x}_i = \operatorname{vec}(X_i)$ are linearly independent, we can write $Y = \sum_{i=1}^{n^2} \alpha_i X_i$. Hence,

$$
C = \sum_{i=1}^{n^2} \alpha_i (AX_i + X_i A + MX_i M) = \sum_{i=1}^{n^2} \alpha_i (1 + \lambda_i) MX_i M.
$$

Since trace(X_j^T (M X_i M)) = δ_{i,j}, it follows that trace(X_i^T C) = α_i(1 + λ_i), from which the result follows. Finally, for C symmetric, we can write C = ∑_{k=1}^n c_k η_k c_k^T, so that trace(X_i^T C) = ∑_{k=1}^n η_k trace(c_k^T X_i^T c_k). For X_i skew-symmetric, each of these addends is zero. Hence, only the terms corresponding to the symmetric matrices X_i need be considered in the sum in i. □

The general results of Proposition 3.1 do not indicate special rank properties, and in fact eigenmatrices of full rank are expected. On the other hand, rank structured eigenmatrices arise in special cases, when for instance the two symmetric matrices A, M have further structure. As an example, consider the case when A, M commute, so that they share the same orthogonal eigenbasis Q [18, Th.1.3.12], and they have eigenvalue matrices Θ = diag(θ₁, ..., θₙ) and Γ = diag(γ₁, ..., γₙ), respectively. Then the eigenpairs have the form (λₖ, QEₖQᵀ) with λₖ = (θᵢ + θⱼ) / (γᵢγⱼ) and Eₖ = eᵢeⱼ, for i, j = 1, ..., n, k = (i - 1)n + j. This characterization immediately follows from substituting A = QΘQᵀ and M = QΓQᵀ in (1.1).

REMARK 3.3. The related problem

$$
(3.2) \qquad AXM + MXA = 2MXM\lambda, \quad X \in \mathbb{R}^{n \times n}
$$

has been thoroughly analyzed in [1] and more recently, e.g., in [26]; see also [9] for a direct application to stability analysis, where A = I is used. In [25],[26] it is shown that for each real eigenpair (λ, q) of (A, M), there exists an eigenpair (λ, X) where vec(X) = q ⊗ q. More specifically, under the assumption that M is nonsingular, (3.2) can be rewritten as

$$
\hat{A}X + X\hat{A}^T = 2\lambda X, \quad \hat{A} = M^{-1}A.
$$
Since both M and A are symmetric and positive definite, $\hat{A}$ admits a full set of right eigenvectors $\{z_i\}_{i=1,...,n}$ and left eigenvectors $\{y_i\}_{i=1,...,n}$ (both sets real) associated with real eigenvalues $\{\theta_i\}_{i=1,...,n}$. Simple eigenvalues have $X = z_i y_j^T$ as corresponding eigenmatrices with eigenvalues $\lambda = (\theta_i+\theta_j)/2$, for $i, j = 1, ..., n$. Due to the structure of $\hat{A}$, the $z_i$'s can be taken to be M-orthogonal, and $y_i = M z_i$. For multiple eigenvalues we can proceed as in the proof of Proposition 3.1, with symmetric and skew-symmetric eigenmatrices. We also remark that in [26] the non-Hermitian case is addressed, and the occurrence of complex eigenvalues is also discussed; this is not the case in our setting.

The structural symmetry and orthogonality properties of Proposition 3.1 can be somehow generalized to the corresponding tensor eigenproblem (1.3). Clearly, orthogonality properties as those in (3.1) still hold, though a different interpretation is required.

**PROPOSITION 3.4.** Consider the eigenvectors $\{\mathbf{x}_i\}_{i=1,...,n^3}$, and their tensor form $\{\mathbf{X}_i\}_{i=1,...,n^3}$. Let the corresponding mode-m (with $m \in \{1, 2, 3\}$) matricization be $[X_1^{(i)}, X_2^{(i)}, ..., X_n^{(i)}]$. Then the matrices $X_k^{(i)}$ satisfy:

i) Symmetry: $X_k^{(i)}$ is either symmetric or skew-symmetric;

ii) Orthogonality:

$$ \sum_{k=1}^{n} \sum_{\ell=1}^{n} M_{k,\ell} \text{trace}(M X_k^{(j)} M X_{\ell}^{(i)}) = \delta_{i,j}. $$

where $\text{trace}(M X_k^{(j)} M X_{\ell}^{(i)}) = \langle (M^{\frac{1}{2}} X_k^{(j)} M^{\frac{1}{2}}), (M^{\frac{1}{2}} X_{\ell}^{(i)} M^{\frac{1}{2}}) \rangle$.

*Proof.* Symmetry could be inspected elementwise, but we prefer to proceed with a less technical matricization approach. We observe that the symmetric structure of both left and right operators allows us to work on the matricization in any of the modes. For simplicity, we consider the first mode.

i) We write the tensor equation in matrix form as

$$ AX + X(A \otimes I_n + I_n \otimes A) = \lambda MX(M \otimes M). $$

For each block of columns $X_j$ we can write

$$ AX_j + X_j A + \sum_{k=1}^{n} A_{j,k} X_k = \lambda M \sum_{k=1}^{n} M_{j,k} X_k M. $$

By transposing both equation sides, we observe that both $X_\ell$ and $X_\ell^T$, $\ell = 1, ..., n$ satisfy the matrix equation. If $X_\ell$ is symmetric or skew-symmetric for all $\ell = 1, ..., n$, then we are done. If $X_\ell \neq X_\ell^T$ then the matrices $V_\ell := X_\ell + X_\ell^T$ and $W_\ell := X_\ell - X_\ell^T$, $\ell = 1, ..., n$ also satisfy the matrix equation and their vectorization spans the same space as the vectorization of $X_\ell, X_\ell^T$.

ii) We only have to unfold the first vector property in (3.1).

Let $\hat{X}_i = [X_1^{(i)}, X_2^{(i)}, ..., X_n^{(i)}]$ with $X_k^{(i)} = X_i(\cdot, \cdot, k)$ and $\hat{X}_j = [X_1^{(j)}, X_2^{(j)}, ..., X_n^{(j)}]$ with $X_k^{(j)} = X_j(\cdot, \cdot, k)$ be the matricizations of the eigenvectors $\mathbf{x}_i, \mathbf{x}_j$, respectively, such that $\mathbf{x}_i^T(M \otimes M \otimes M)\mathbf{x}_j = \delta_{i,j}$. Unfolding the Kronecker product we
obtain

$$
\begin{align*}
\delta_{i,j} &= \boldsymbol{x}_i^T (M \otimes M \otimes M) \boldsymbol{x}_j = \operatorname{trace}(M \hat{X}_j (M \otimes M) \hat{X}_i^T) \\
&= \operatorname{trace}(M \sum_{k=1}^n (X_k^{(j)} (\sum_{\ell=1}^n M_{k,\ell} M X_{\ell}^{(i)}))) \\
&= \operatorname{trace}(M \sum_{k=1}^n (X_k^{(j)} M (\sum_{\ell=1}^n M_{k,\ell} X_{\ell}^{(i)}))).
\end{align*}
$$

Reordering terms yields the desired relation.

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

**4. On the Kronecker pair** (*A* ⊗ *I* + *I* ⊗ *A*, *M* ⊗ *I* + *I* ⊗ *M*). Although (1.2) can be written as the more familiar generalized vector problem in (2.2), the original matrix form is particularly convenient in certain PDE contexts, as we will see later, in addition to the already cited Hopf bifurcation analysis [25]. In there, the problem

$$
(4.1) \quad AXB + BXA = \lambda(BXN + NXB), \quad X \in \mathbb{R}^{n \times n}
$$

was thoroughly investigated; this corresponds to (1.2) for $B = I$, $N = M$ and $A$
symmetric and positive definite, and to (2.4) for $B = M^{-1}$ and $N = I$; hence certain
results in [25] can be readily employed in our setting. Indeed, let us collect all terms
in (1.2) on the left-hand side as

$$
(4.2) \qquad (A - \lambda M)X + X(A - \lambda M) = 0,
$$

and let $G(\lambda) := A - \lambda M$. The original eigenproblem can be related to the solvability of the matrix equation $G(\lambda)X + XG(\lambda) = 0$. The following result is just a reformulation of [25, Th.2.2] under our hypotheses.

PROPOSITION 4.1. With the notation above, we have that (λ_k, X_k) is an eigenpair of (1.2) if and only if one of the following holds

i) There exist two eigenpairs $(\zeta_i, z_i)$, $(\zeta_j, z_j)$ of $G(\lambda_k)$, with $i \neq j$, such that $\zeta_i = -\zeta_j$ and $X_k = [z_i, z_j]\Gamma[z_j, z_i]^T = z_i\gamma_1 z_j^T + z_j\gamma_2 z_i^T$, for any nonzero $\Gamma = \operatorname{diag}([\gamma_1, \gamma_2])$;

ii) There exists an eigenpair $(\zeta, z)$ of $G(\lambda_k)$ with simple $\zeta = 0$ such that $X_k = \bar{z}\bar{z}^T$.

Proof. For a given $\lambda_k$, a solution to $G(\lambda_k)X + XG(\lambda_k) = 0$ is determined in the null space of $G(\lambda_k) = I \otimes G(\lambda_k) + G(\lambda_k) \otimes I$, that is, $\lambda_k$ should be such that $G(\lambda_k)$ is singular. The eigenvalues of $G(\lambda_k)$ are given by $\zeta_i + \zeta_j$, where $\zeta$'s are the eigenvalues of $G(\lambda_k)$, for all combinations of $i, j \in \{1, \dots, n\}$. Hence, $G(\lambda_k)$ is singular if and only if there exist two eigenvalues such that $\zeta_i + \zeta_j = 0$. The rest of the result follows by substitution.
$\square$

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

A remarkable consequence of Proposition 4.1 is that the eigenmatrix $X_k$ has
always rank at most two. We also notice that in the nonsymmetric case, the charac-
terization of Proposition 4.1 motivated the authors of [25] to use this matrix equation
in the so-called bialternate product method for detecting Hopf bifurcations [14].

We can naturally extend this result to the tensorial setting, with the generalized
eigentensor problem

$$
(4.3) \quad (A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A)x = \\
\qquad \lambda(M \otimes I_n \otimes I_n + I_n \otimes M \otimes I_n + I_n \otimes I_n \otimes M)x,
$$

by exploiting the same matrix $G(\lambda)$ defined above, so as to write

$$
(G(\lambda) \otimes I_n \otimes I_n + I_n \otimes G(\lambda) \otimes I_n + I_n \otimes I_n \otimes G(\lambda)) \boldsymbol{x} = 0,
$$
or $G(\lambda)x = 0$ in short.

**PROPOSITION 4.2.** With the notation above, we have that $(\lambda_k, x_k)$ is an eigenpair of (4.3) if and only if one of the following holds

i) There exist three eigenpairs $(\theta_i, z_i)$, $(\theta_j, z_j)$, $(\theta_l, z_l)$ of $G(\lambda_k)$ with $i, j, l$ not equal, such that $\theta_i + \theta_j + \theta_l = 0$ and $x_k = \gamma_1 z_i \otimes z_j \otimes z_l + \gamma_2 z_i \otimes z_l \otimes z_j + \gamma_3 z_j \otimes z_i \otimes z_l + \gamma_4 z_j \otimes z_l \otimes z_i + \gamma_5 z_l \otimes z_i \otimes z_j + \gamma_6 z_l \otimes z_j \otimes z_i$, for any nonzero $\Gamma = \text{diag}(\gamma_1, \gamma_2, \gamma_3, \gamma_4, \gamma_5, \gamma_6)$;

ii) There exists an eigenpair $(\theta_i, z_i)$ of $G(\lambda_k)$ with simple $\theta_i = 0$ so that $x_k = z_i \otimes z_i \otimes z_i$.

*Proof.* For a given $\lambda_k$, a solution $x$ to $G(\lambda_k)x = 0$ is determined in the null space of the coefficient matrix $G(\lambda_k)$. The eigenvalues of $G(\lambda_k)$ are given by $\theta_i + \theta_j + \theta_l$, where $\theta$'s are the eigenvalues of $G(\lambda_k)$, for all combinations of $i, j, l$. Hence, $G(\lambda_k)$ is singular if and only if there exist three eigenvalues such that $\theta_i + \theta_j + \theta_l = 0$. The rest of the result follows by substitution. $\square$

We remark that the non-nullity of the $\gamma_i$s in Proposition 4.2 is related to the multiplicity of the eigenvalues $\theta$'s in $G(\lambda_k)$. In other words, each $x_k$ is composed by up to six terms, with the maximum occurring when the $\theta$'s are such that $\{z_i, z_j, z_l\}$ spans a three dimensional subspace.

**5. Operator preconditioning.** In this section we consider in detail an application of the problem (1.2) to the spectral analysis of preconditioned coefficient matrices in the solution of large linear systems stemming from the discretization of linear PDEs. After a short description of the discretization procedure, tailored to our setting, we report on several eigenvalue properties whose derivation takes advantage of the structure analyzed so far.

**5.1. The discretized problem.** We consider the following equation

$$-\nabla \cdot (\kappa(x, y)\nabla u) = f, \quad (x, y) \in \Omega = (0, 1)^2,$$

with either Dirichlet or mixed (Dirichlet plus Neumann) boundary conditions, and $\kappa(x, y) = \text{diag}(a(x), b(y))$. Consider the following (mixed) conditions:

$$ (5.1) \quad u_x(1, y) = 0, y \in [0, 1], \quad u_y(x, 1) = 0, x \in [0, 1], \quad u(x, y) = 0 \text{ elsewhere on } \partial\Omega. $$

To determine the finite difference discretization of the PDE, let us consider the grid-points $(x_i, y_j)$, $i = 1, \dots, n$ of a uniform discretization of $\Omega$. Here we focus on mixed conditions, so the boundary nodes associated with Neumann conditions are included; see Appendix A for pure Dirichlet conditions. We define

$$ (5.2) \quad D_0 := \begin{bmatrix} 1 \\ -1 & \ddots \\ \vdots & \ddots \\ \vdots & \ddots \\ -1 & 1 \end{bmatrix} \in \mathbb{R}^{n \times n}, \quad \begin{array}{l} S_{a,0} := \text{diag}(a_{i'})_{i=1,\dots,n} \in \mathbb{R}^{n \times n}, \\ S_{b,0} := \text{diag}(b_{i'})_{i=1,\dots,n} \in \mathbb{R}^{n \times n}. \end{array} $$

where $a_{i'} = a(x_{i-\frac{1}{2}})$, are the collocation values of the coefficients, see, e.g., [28, formula (2.16)]; the same indexing is used for $b(y)$. The discretized equation is given by

$$ (5.3) \quad A_a U + U A_b = F, \quad A_k = D_0^T S_{k,0} D_0, \quad k = a, b, $$
where $F = f(x_i, y_j) + b.c.$ and $U \approx u(x_i, y_j)$. Note that the Poisson equation, corresponding to constant $a(x) = 1, b(y) = 1$ is given by<fnref n="1" />

$$ (5.4) \qquad MU + UM = F, \quad M = D_0^T D_0. $$

Classically, equation (5.3) is vectorized into the equation

$$ Au = f, \text{ with } u = \text{vec}(U), \text{ and } f = \text{vec}(F), $$

and $A = I \otimes A_a + A_b \otimes I$; see, e.g., [30] for some early references. A largely studied preconditioner for this linear system is the matrix obtained by discretization using, e.g., finite elements or finite differences, of the constant version of the operator, which often just corresponds to the operator $P: U \mapsto MU + UM$ on the left-hand side of (5.4), with the same boundary conditions [23]. Hence, the system actually solved is $P^{-1}Ax = P^{-1}f$, where $P$ is the vectorized version of $P$. These classical procedures have been well studied, see, e.g., [10],[27], and [22] for a spectral analysis of this strategy on convection-diffusion operators and [31] for an example of combined discretization techniques. The whole discussion can be generalized to more than two space variables; the multiarray derivation in Kronecker form can be found, for instance, in [20, section 5.6.4].

Towards an accurate assessment of the preconditioning effectiveness, a spectral analysis of the preconditioned coefficient matrix $P^{-1}A$ is crucial. In particular, not only the spectral interval is of interest, but also the distribution of the eigenvalues themselves in that interval. The problem has been analyzed in the literature by making fundamental connections with the continuous problem and the employed discretization [2],[10],[15],[22]. In particular, in the recent work [11] the authors have proved what they called a *pairing* between the values of operator coefficient function $\kappa$ and the eigenvalues of the discretized algebraic eigenvalue problem.

For discretization by finite differences on a rectangular domain or a parallelepiped, this problem now falls into our framework, with $P$ playing the role of $M$, as the eigenvalue problem of interest is given by (2.2), (1.4), respectively. In the following we will thus replace $P$ with $M$, for notational consistency with respect to the previous sections.

**5.2. Spectral analysis of the preconditioned algebraic problem.** In the following we present results on the discretized problem after the same mixed boundary conditions have been imposed to both the coefficient matrix and the preconditioner [23]. The slightly technically more involved case of Dirichlet boundary conditions is postponed to Appendix A.

**REMARK 5.1.** The vector form obtained in (2.4) after multiplication by $M^{-\frac{1}{2}} \otimes M^{-\frac{1}{2}}$ yields a convenient form for our analysis. In matrix form, this multiplication transforms (2.4) into

$$ \hat{A}XM^{-1} + M^{-1}X\hat{A} = \lambda(XM^{-1} + M^{-1}X) $$

with $\hat{A} = M^{-\frac{1}{2}}AM^{-\frac{1}{2}}$. If $M$ and $A$ are spectrally equivalent, then $\hat{A}$ is spectrally equivalent to the identity matrix, which leads to conclude that the eigenvalues $\lambda$s do not depend on the problem size, but only on the constants yielding spectral equivalence of $A$ and $M$.

<footnote n="1">¹The matrix $M$ corresponds to minus the standard one-dimensional second order derivative matrix tridiag($-1, 2, -1$), except for the last diagonal element, equal to one, that takes into account the Neumann boundary condition at the right end of the interval.</footnote>
The facts highlighted in the remark above are now very well known, and they have been derived for different discretizations; see, e.g., the very thorough analyses in [10],[19],[22],[23].

In the following we will derive again this invariance property using the convenient matrix-oriented form derived from finite differences. In practice, this formulation allows us to visually distinguish at the linear algebra level between the operator coefficients ($\kappa$ dependence) and the contribution from the approximate derivatives ($h$ dependence). This way we can parallel a standard strategy in the discretization of weak formulations of elliptic problems, where preconditioners are first represented in the continuous setting, and then appropriately discretized, typically by finite element methods; see, e.g., [8],[24],[15]. Hence, in the one-dimensional case for instance,

$$ \nabla \cdot \kappa \nabla \rightarrow D_0^T S_{a,0} D_0. $$

To lighten the presentation, here we illustrate the result for $b(y_j) = a(x_j)$; see Remark 5.3 for the case $b \neq a$.

In the one-dimensional case with mixed boundary conditions, the spectrum of the pair $(D_0^T S_{a,0} D_0, D_0^T D_0)$ readily corresponds to the diagonal elements of $S_{a,0}$, that is the values of $a(x)$ at the staggered nodes. In the two-dimensional case the evaluation is less accurate, though the general picture is preserved, as expected by the related literature.

**THEOREM 5.2.** Let $a$ be the coefficient function in each derivative and assume that $\alpha \le a(x) \le \beta$ for all $x \in [0, 1]$. Then $\mathbf{A}$ and $\mathbf{M}$ are spectrally equivalent, and it holds that

$$ \alpha \le \frac{\mathbf{x}^T \mathbf{A} \mathbf{x}}{\mathbf{x}^T \mathbf{M} \mathbf{x}} \le \beta, \quad \forall 0 \ne \mathbf{x} \in \mathbb{R}^{n^2}. $$

*Proof.* We first note that from $a(x) \in [\alpha, \beta] \ \forall x \in [0, 1]$, it follows that $\text{spec}(S_0) \subseteq [\alpha, \beta]$. Using the spectral properties of the Kronecker product, it holds that

$$ (5.5) \qquad \mathbf{M} = [D_0^T \otimes I_n, I_n \otimes D_0^T] \begin{bmatrix} D_0 \otimes I_n \\ I_n \otimes D_0 \end{bmatrix} := \mathbf{D}_0^T \mathbf{D}_0 $$

$$ (5.6) \qquad \mathbf{A} = [D_0^T \otimes I_n, I_n \otimes D_0^T] \begin{bmatrix} S_0 \otimes I_n & \\ & I_n \otimes S_0 \end{bmatrix} \begin{bmatrix} D_0 \otimes I_n \\ I_n \otimes D_0 \end{bmatrix} := \mathbf{D}_0^T \mathbf{S}_0 \mathbf{D}_0 $$

with $S_0 \in \mathbb{R}^{n^2 \times n^2}$ a diagonal matrix. Using the spectral properties of the Kronecker sum, it also follows that $\text{spec}(S_0) \subseteq [\alpha, \beta]$. For $0 \ne \mathbf{x} \in \mathbb{R}^{n^2}$ we have

$$ \frac{\mathbf{x}^T \mathbf{A} \mathbf{x}}{\mathbf{x}^T \mathbf{M} \mathbf{x}} = \frac{\mathbf{x}^T \mathbf{D}_0^T \mathbf{S}_0 \mathbf{D}_0 \mathbf{x}}{\mathbf{x}^T \mathbf{D}_0^T \mathbf{D}_0 \mathbf{x}} \underset{0 \ne \mathbf{y}=\mathbf{D}_0 \mathbf{x}}{=} \frac{\mathbf{y}^T \mathbf{S}_0 \mathbf{y}}{\mathbf{y}^T \mathbf{y}} \in [\alpha, \beta]. \quad \square $$

We remark that completely analogous formulations can be obtained in the tensor case [20], though the general multiarray case does not allow for a ready-to-use analysis of the remaining parameter dependence. Indeed, the spectral properties of the preconditioned operator, and thus its effectiveness, still depend on the width of the interval $[\alpha, \beta]$, and on how the eigenvalues distribute within this interval (see section 5.4 for the three-dimensional case).

**REMARK 5.3.** Theorem 5.2 can be generalized to different coefficient functions, that is $a(x_j) \ne b(x_j)$, because the only information used is the decomposition together with the diagonal form of $S_0$, which remains true in this case. As a consequence, with the notation in (5.3), the eigenvalues $\lambda$ of $A_a X + X A_b = \lambda(MX + XM)$ satisfy

$$ \lambda \in \left( 2 \min_{x \in (0,1)} \{a(x), b(x)\}, 2 \max_{x \in (0,1)} \{a(x), b(x)\} \right). $$
**Piecewise constant coefficients.** An interesting characterization can be derived for $a(x)$ piecewise constant in $[\alpha, \beta]$, showing that the eigenvalues do not distribute uniformly in $[\alpha, \beta]$, as, for instance, empty gaps arise. Though the problem with piecewise constant values could be addressed with more powerful techniques (see, e.g., [29]), our approach would also allow one to naturally predict the behavior of continuous coefficients with steep value jumps by means of a perturbation analysis – we do not report such analysis here. Our results complement several studies in the literature, and in particular a recent analysis for the matrices stemming from finite element discretization presented in [11]. The use of a finite difference grid provides us with a somewhat simpler setting.

The role of jumps in the PDE coefficients in the numerical solution of preconditioned linear systems has been highly regarded, and preconditioners tailored to attack possible misbehaviors have been proposed, see, e.g., [5],[6] and their references; in general, algebraic multigrid and domain decomposition preconditioners seem to be robust with respect to jumps, see, e.g., [7],[4],[13], and [33] for further references and for additional numerical linear algebra considerations. Here we are interested in characterizing how these jumps influence the eigenvalue distribution in a simple setting.

Let $a(x)$ be piecewise constant, with a single jump at $x = x_0 \in (0,1)$, that is $a(x) = \alpha$ for $x \le x_0$ and $a(x) = \beta$ for $x > x_0$, with $\alpha \neq \beta$, so that after discretization, $S_0 := \text{blkdiag}(\alpha I_{n_1}, \beta I_{n_2})$, where $n_1 + n_2 = n$, and the actual values of $n_1, n_2$ depend on the location of $x_0$ in the given interval. Thanks to Proposition 4.1, we can say that the pencil $(A, M)$ has $n_1^2$ eigenpairs $(\alpha, X)$, $k = 1, \dots, n_1^2$ such that $(\alpha, z_i)$ and $(\alpha, z_j)$ are an eigenpair of $(A, M)$, and $X = [z_i, z_j] \Gamma[z_j, z_i]^T$ for all possible combinations of $i$ and $j$, and $\Gamma$ symmetric, either positive definite or indefinite. In particular, for a fixed eigenvalue $\lambda_k$, there are two independent eigenvectors, which are $x_1 = \text{vec}([z_1, z_2] \Gamma_1[z_2, z_1]^T)$ and $x_2 = \text{vec}([z_1, z_2] \Gamma_2[z_2, z_1]^T)$ with the two $\Gamma$s having different signature. Analogously, the pencil has $n_2^2$ eigenvalues equal to $\beta$, with corresponding eigenvectors.

The following theorem shows that all other eigenvalues are a convex linear combination of $\alpha$ and $\beta$. To this end, let $n_1, n_2 \in \mathbb{N}$ with $n_1 + n_2 = n$, and define

$$P_1^{(0)} = \text{blkdiag}(P_1, 0_{n_2}) \in \mathbb{R}^{n \times n}, \quad P_1 = D_0^T D_0 \in \mathbb{R}^{n_1 \times n_1}$$

$$\text{and } P_2^{(0)} = \text{blkdiag}(0_{n_1-1}, P_2) \in \mathbb{R}^{n \times n},$$

$$P_2 = \begin{bmatrix} 1 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & -1 & 1 & \end{bmatrix} \in \mathbb{R}^{(n_2+1) \times (n_2+1)}.$$

We note that both matrices $P_1^{(0)}, P_2^{(0)}$ are symmetric and positive semidefinite and that in the sum $P_1^{(0)} + P_2^{(0)}$ the two nonzero blocks overlap in the $(n_1, n_1)$ element. Then we define

$$P_1 = I_n \otimes P_1^{(0)} + P_1^{(0)} \otimes I_n, \quad P_2 = I_n \otimes P_2^{(0)} + P_2^{(0)} \otimes I_n.$$

We observe that $(P_2, P_1)$ is a regular pair, since $P_2 + P_1 = M$ is nonsingular, with non-negative or infinite eigenvalues. Moreover, the distribution of its eigenvalues
depends only on the splitting of the domain given by the coefficient discontinuity, and not on the actual values of the operator $a(x)$.

**THEOREM 5.4.** Let $\theta_j, j = 1, \dots, k$ be the non-negative real eigenvalues of the pencil $(\mathbf{P}_2, \mathbf{P}_1)$, that is $\theta_j \neq \infty$. Then the eigenvalues $\lambda$ of $(\mathbf{A}, \mathbf{M})$ satisfy

$$ \lambda \in \left\{ \beta, \frac{\beta\theta_1 + \alpha}{\theta_1 + 1}, \dots, \frac{\beta\theta_k + \alpha}{\theta_k + 1} \right\}. $$

*Proof.* We can write $\mathbf{A} = \alpha \mathbf{P}_1 + \beta \mathbf{P}_2$, $\mathbf{M} = \mathbf{D}_0^T \mathbf{D}_0 = \mathbf{P}_1 + \mathbf{P}_2$. The problem $\mathbf{A}x = \lambda \mathbf{M}x$ can be written as $(\alpha - \lambda)\mathbf{P}_1x + (\beta - \lambda)\mathbf{P}_2x = 0$, that is, for $\lambda \neq \beta$,

$$ \mathbf{P}_2x = \theta \mathbf{P}_1x, \quad \text{with} \quad \theta = -\frac{\alpha - \lambda}{\beta - \lambda}. $$

Writing down $\lambda$ in terms of $\theta$ yields $\lambda = \frac{\beta\theta+\alpha}{\theta+1}$. $\square$

The result also includes the case $\lambda = \alpha$ for $\theta = 0$, and shows that the sought after linear combination is $\lambda = \frac{\theta}{\theta+1}\beta + \frac{1}{\theta+1}\alpha$. Subtracting $(\alpha + \beta)/2$ on both sides and rearranging terms, we obtain

$$ \frac{\lambda - \frac{\alpha + \beta}{2}}{\frac{\beta - \alpha}{2}} = \frac{\theta - 1}{\theta + 1}. $$

The quantity on the left-hand side may be viewed as a *standardized* eigenvalue around the interval mean, taking values in $[-1, 1]$ for $\theta \ge 0$. The distribution of $\theta$ around the value 1 provides information on the distribution of $\lambda$ around the middle value $(\alpha + \beta)/2$.

To proceed further, it is thus useful to obtain a refined analysis on the location of the eigenvalues $\theta$ of $(\mathbf{P}_2, \mathbf{P}_1)$.

**PROPOSITION 5.5.** Let $\lambda_{\max}(P_i)$ be the largest eigenvalue of $P_i$, $i = 1, 2$, and $\lambda_{\min}(P_i)$ be the smallest nonzero eigenvalue of $P_i$, $i = 1, 2$. Then the finite nonzero eigenvalues $\theta$ of $(\mathbf{P}_2, \mathbf{P}_1)$ satisfy

$$ \frac{1}{2} \frac{\lambda_{\min}(P_2)}{\lambda_{\max}(P_1)} \le \theta \le 2 \frac{\lambda_{\max}(P_2)}{\lambda_{\min}(P_1)}. $$

*Proof.* The proof mainly dwells with the singularity of both matrices in the pencil, and it otherwise follows from standard results. It is postponed to Appendix B. $\square$

We notice that the extreme eigenvalues of $P_1, P_2$ are known analytically, so that in the bounds of Proposition 5.5 we can write

$$ (5.7) \qquad \frac{\lambda_{\min}(P_2)}{\lambda_{\max}(P_1)} = \frac{1 - \cos \frac{\pi}{2n_2+1}}{1 - \cos \frac{(2(n-n_2)-1)\pi}{2(n-n_2)+1}}, \qquad \frac{\lambda_{\max}(P_2)}{\lambda_{\min}(P_1)} = \frac{1 - \cos \frac{(2n_2-1)\pi}{2n_2+1}}{1 - \cos \frac{\pi}{2(n-n_2)+1}}, $$

thus allowing a better localization of the eigenvalues $\lambda$ distinct from $\alpha, \beta$ of the pair $(\mathbf{A}, \mathbf{M})$, in the interval $[\alpha, \beta]$. In particular, for $n_2 \ll n$, a gap can be observed between $\lambda = \alpha$ (corresponding to $\theta = 0$) and the next distinct eigenvalue, as shown by the lower bound for $\theta$. For $n_2 \approx n$, a gap can be observed between $\lambda = \beta$ (corresponding to an infinite $\theta$) and the previous distinct eigenvalue, monitored by the upper bound for $\theta$. Finally, our computational experiments have shown that for
$n_2 \approx (n - n_2) \approx n/2$ a cluster of eigenvalues $\lambda$s can be observed around the interval midpoint; see section 5.3.

The obtained gap is important in analyzing the convergence rate of the Preconditioned Conjugate Gradient (PCG) method. Indeed, after the first iterations, the algorithm behaves as if the spectral interval were reduced to the “effective” spectral interval $[\lambda_{\min}, \lambda_{\max}]$ instead of $[\alpha, \beta]$, where $\lambda_{\min}, \lambda_{\max}$ are the eigenvalues of $(A, M)$ closest to $\alpha, \beta$, respectively, but distinct from them.

<figure><img src="image_4.png" /><figcaption>FIG. 5.1. Example 5.6, for $n = 50$, $n_2 = 2:2:48$, $\alpha = 4$, $\beta = 10$.</figcaption></figure>

**5.3. Numerical experiments.** The next examples help us illustrate our findings in the previous section.

**EXAMPLE 5.6.** We consider the preconditioned algebraic problem as if it came from the PDE $-(a(x)u_x)_x - (a(y)u_y)_y = f$ with mixed boundary conditions. While discretizing using $n$ nodes in each direction, the coefficients are set to $a(x_j) = \beta$ in the last $n_2$ discretization nodes and to $\alpha$ elsewhere. Varying $n_2$ corresponds to a smaller or larger portion of the domain where the value $\beta$ is taken, that is, the ratio $n_2/n$ changes. Clearly, in general this variable setting is not representative of a physical problem, however it serves well to our argumentation purposes. We consider $\alpha = 4$, $\beta = 10$ and $n = 50$ nodes in each direction. Figure 5.1 reports the eigendistribution for $n_2 = 2, 4, 6, \dots, 48$. Each set of eigenvalues is depicted as a monotonically increasing curve with values from $\alpha$ to $\beta$, with increasing values of $n_2$ from left to right. The eigenvalue gap near $\beta$ for $n_2 \ll n$ is clearly visible among the leftmost curves, while the rightmost curves display the gap near $\alpha$. We also notice that for values $n_2 \approx 50 = n/2$ more eigenvalues take values around $7 = (\alpha + \beta)/2$ (depicted by the curve bending in the figure). This particular behavior will be analyzed in greater detail elsewhere.

Table 5.1 shows the upper bound for the largest eigenvalue $\lambda = \frac{\theta}{\theta+1}\beta + \frac{1}{\theta+1}\alpha$ not equal to $\beta$, as $n_2$ varies. The bound is obtained by using the upper estimate of $\theta$ in (5.7) together with Proposition 5.5. The bound appears to be quite sharp. Corresponding similar values can be obtained for $\alpha$.
<table><caption>TABLE 5.1</caption><thead><tr><th>n<sub>2</sub></th><th>&lambda;<sub>max</sub></th><th>&lambda;<sub>bound</sub></th><th>n<sub>2</sub></th><th>&lambda;<sub>max</sub></th><th>&lambda;<sub>bound</sub></th></tr></thead><tbody><tr><td>2</td><td>9.9976</td><td>9.9988</td><td>⋮</td><td>⋮</td><td>⋮</td></tr><tr><td>6</td><td>9.9980</td><td>9.9990</td><td>34</td><td>9.9868</td><td>9.9932</td></tr><tr><td>10</td><td>9.9977</td><td>9.9988</td><td>38</td><td>9.9772</td><td>9.9882</td></tr><tr><td>14</td><td>9.9972</td><td>9.9986</td><td>42</td><td>9.9516</td><td>9.9745</td></tr><tr><td>18</td><td>9.9965</td><td>9.9982</td><td>46</td><td>9.8384</td><td>9.9108</td></tr><tr><td>⋮</td><td>⋮</td><td>⋮</td><td>48</td><td>9.5400</td><td>9.7263</td></tr></tbody></table>

Example 5.6. Largest eigenvalue distinct from $\beta$ and its upper bound, for $n = 50$.

EXAMPLE 5.7. We consider again the data from Example 5.6, but with different values of $\alpha$ and $\beta$, namely $\alpha = 1$, $\beta = 10^4$. We set $f$ to have random values (scaled to have unit norm). If Figure 5.2 we report the convergence of the preconditioned conjugate gradient method applied to $Ax = f$ with $M$ as preconditioner for different values of $n$, so that the problem to be solved has dimensions $n^2 \times n^2$. The shown curves refer to the energy-norm of the error. We also report (dashed line) the expected asymptotic convergence rate $\rho^k$ as the iteration $k$ progresses, with

$$ \rho = \frac{\sqrt{\mathrm{cond}_{\mathrm{eff}}} - 1}{\sqrt{\mathrm{cond}_{\mathrm{eff}}} + 1}, \quad \mathrm{cond}_{\mathrm{eff}} := \frac{\lambda_{\max}}{\lambda_{\min}}; $$

here $\mathrm{cond}_{\mathrm{eff}}$ denotes the “effective” condition number of the preconditioned matrix $M^{-1}A$, obtained by excluding the eigenvalues $\alpha, \beta$. The eigenvalues of the preconditioned matrix were computed with the finest grid, though the theory ensures that mesh independence is preserved, as it is also clear from the actual PCG convergence curves. The plot illustrates that the asymptotic ratio computed with the interior spectral interval is very descriptive of most of the convergence history. The initial almost stagnating phase instead conforms with the behavior expected by a condition number of the order of $\beta/\alpha$, much larger than $\mathrm{cond}_{\mathrm{eff}}$, associated with the whole spectral interval. We refer, e.g., to [3],[32],[2],[11] for an analysis of superlinear convergence of the conjugate gradient method and its relevance in discretized elliptic equations. The left plot of Figure 5.2 refers to $n_2 = 4$, whereas the right plot considers $n_2 = 10$.

It is also interesting to observe the different behavior for $n_2 \approx n/2$ in Figure 5.3. According to the plot in Figure 5.1, the central part of the spectral interval is more populated, or, in other words, the two interval ends contain fewer eigenvalues, though no gaps occur. We have observed that these more isolated eigenvalues are gradually identified by the space generated by PCG, accelerating convergence as the iterations proceed. This determines a progressive steeper curve (sequentially increasing superlinear convergence), illustrated by the solid curve in the plot.

EXAMPLE 5.8. We consider the PDE $-(a(x)u_x)_x - (a(y)u_y)_y = f$ with $f(x, y) = \sin(x\pi)\cos(y\pi)$, $(x, y) \in [0, \pi]^2$ and Dirichlet boundary conditions, with the coefficients choice

$$ a(x) = \begin{cases} 10^4 & x < 0.5 \\ 1 & x \ge 0.5. \end{cases} $$

Figure 5.4(left) displays the approximate solution obtained for $n = 60$ nodes in each direction, showing how the coefficients' jumps affect the obtained approximate solution.
<figure><img src="image_5.png" /><figcaption>FIG. 5.2. Example 5.6. Convergence of PCG. Left: $n_2 = 4$. Right: $n_2 = 10$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5.3. Example 5.6. Convergence of PCG for $n_2 \ll n$ (dotted line) and for $n_2 \approx n/2$ (solid line).</figcaption></figure>

We then compare the eigenvalue distribution for the pair $(\mathbf{A}, \mathbf{M})$ with that obtained with the continuous coefficients

$$ \hat{a}(x) = \begin{cases} 10^4 & x < 0.5 \\ 10^4 \exp(25(0.5 - x)) + 1.1 & x \in [0.5, 1.1) \\ 1 & x > 1.1. \end{cases} $$

Figure 5.4(center) displays the graphs in $[0, \pi]^2$ for a and $\hat{a}$ as defined above. The dashed curve in Figure 5.4(right) shows the distribution of the eigenvalues of the pair $(\mathbf{A}, \mathbf{M})$ for the piecewise constant coefficients $a(x)$, together with that obtained with $\hat{a}(x)$, $x \in [0, \pi]$. The spectral behavior for the settings with piecewise and continuous functions is extremely similar, including the clustering around the interval center.

**5.4. The three-dimensional case.** Consider the elliptic problem $-(au_x)_x - (bu_y)_y - (cu_z)_z = f$ in $\Omega = (0, 1)^3$, equipped with properly chosen Neumann-Dirichlet boundary conditions. Though our reasoning can be extended to purely Dirichlet boundary conditions, we focus here on the former case. We show that the two-dimensional study can be extended to the three-dimensional setting. Once again, to simplify the presentation, we consider all equal coefficients. The general case will then
<figure><img src="image_7.png" /><figcaption>FIG. 5.4. Example 5.8. Left: solution U from $\mathbf{u} = \mathbf{A}^{-1}\mathbf{f}$, $\mathbf{u} = \text{vec}(U)$. Center: graphs of $a(x)$ and $\hat{a}(x)$. Right: eigenvalues of the pair $(\mathbf{A}, \mathbf{M})$ for $\mathbf{A}$ defined using $a(x)$ or $\hat{a}(x)$.</figcaption></figure>

follow the path of Remark 5.3. Let $A = D_0^T S_0 D_0 \in \mathbb{R}^{n \times n}$ and

$$ \begin{aligned} \mathbf{A} &= A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A \in \mathbb{R}^{n^3 \times n^3} \\ \mathbf{M} &= M \otimes I_n \otimes I_n + I_n \otimes M \otimes I_n + I_n \otimes I_n \otimes M \in \mathbb{R}^{n^3 \times n^3} \end{aligned} $$

**THEOREM 5.9.** Let $a = a(x)$ be the coefficient function and suppose $a(x) \in [\alpha, \beta]$ for $x \in (0, 1)$. Then $\mathbf{A}$ and $\mathbf{M}$ are spectrally equivalent, and it holds

$$ \alpha \le \frac{\mathbf{x}^T \mathbf{A} \mathbf{x}}{\mathbf{x}^T \mathbf{M}_0 \mathbf{x}} \le \beta, \quad \forall \mathbf{x} \in \mathbb{R}^{n^3}, \mathbf{x} \ne 0. $$

*Proof.* We define

$$ \mathcal{D}_0 := \begin{bmatrix} D_0 \otimes I_n \otimes I_n \\ I_n \otimes D_0 \otimes I_n \\ I_n \otimes I_n \otimes D_0 \end{bmatrix}, \quad \mathcal{S} := \begin{bmatrix} S_0 \otimes I_n \otimes I_n & & \\ & I_n \otimes S_0 \otimes I_n & \\ & & I_n \otimes I_n \otimes S_0 \end{bmatrix}, $$

with $\mathcal{D}_0 \in \mathbb{R}^{3n^3 \times n^3}$ and $\mathcal{S} \in \mathbb{R}^{3n^3 \times 3n^3}$. Observe that $\mathbf{A} = \mathcal{D}_0^T \mathcal{S} \mathcal{D}_0$ and $\mathbf{M} = \mathcal{D}_0^T \mathcal{D}_0$. Hence the proof follows the same lines as that of Theorem 5.2, since $\text{spec}(\mathcal{S})$ is contained in the interval including all values of the coefficient function. □

Once again, the extension of Theorem 5.4 is immediate, after having defined

$$ \mathbf{P}_i = I_n \otimes I \otimes P_i^{(0)} + I_n \otimes P_i^{(0)} \otimes I_n + P_i^{(0)} \otimes I_n \otimes I_n, \quad i = 1, 2. $$

**THEOREM 5.10.** Let $\theta_j$, $j = 1, \dots, k$ be the nonnegative real eigenvalues of the pencil $(\mathbf{P}_2, \mathbf{P}_1)$, that is $\theta_j \ne \infty$. Then the eigenvalues $\lambda$ of $(\mathbf{A}, \mathbf{M})$ satisfy

$$ \lambda \in \left\{ \beta, \frac{\beta\theta_1 + \alpha}{\theta_1 + 1}, \dots, \frac{\beta\theta_k + \alpha}{\theta_k + 1} \right\}. $$

These results can be formulated in any dimension $d$, with $d \ge 2$, whenever the discretization leads to the same Kronecker structure.

**6. Conclusions.** Generalized matrix eigenvalue problems provide a rich source for structural properties, which would be hardly uncovered by the apparently more accessible vectorized form. We have described several of these properties, and illustrated in detail the occurrence of this type of problems in the well established spectral analysis associated with operator preconditioning of elliptic problems on rectangular and
parallelepipedal domains, when finite differences are used. In fact, other discretization strategies relying on tensorial approximation spaces such as IGA or spectral methods may lead to similar frameworks. Possible generalizations of our applied analysis include exploring these methodologies, together with the adaptation to non-self adjoint elliptic operators.

**Acknowledgments.** The first author is a member of Indam-GNCS. Its support is gratefully acknowledged.

REFERENCES

[1] F. Atkinson. *Multiparameter Eigenvalue Problems*. Academic Press, New York, 1972.

[2] O. Axelsson and J. Karátson. Mesh independent superlinear PCG rates via compact-equivalent operators. *SIAM J. Numer. Anal.*, 45(4):1945–1516, 2007.

[3] O. Axelsson and G. Lindskog. On the rate of convergence of the preconditioned conjugate gradient method. *Numer. Math.*, 48:499–523, 1986.

[4] Markus Blatt. *A Parallel Algebraic Multigrid Method for Elliptic Problems with Highly Discontinuous Coefficients*. PhD thesis, Heidelberg Universität (D), 2 July 2010.

[5] Xing Cai, Bjørn Fredrik Nielsen, and Aslak Tveito. An analysis of a preconditioner for the discretized pressure equation arising in reservoir simulation. *IMA Journal of Numerical Analysis*, 19:291–316, 1999.

[6] S. V. Dolgov, B. N. Khoromskij, I. V. Oseledets, and E. Tyrtyshnikov. A reciprocal preconditioner for structured matrices arising from elliptic problems with jump coefficients. *Lin. Alg. Appl.*, 436:2980–3007, 2012.

[7] M. Dryja, M. Sarkis, and Widlund O. Multilevel Schwarz methods for elliptic problems with discontinuous coefficients in three dimensions. *Numer. Math.*, 72:313–348, 1996.

[8] Howard C. Elman, David J. Silvester, and Andrew J. Wathen. *Finite Elements and Fast Iterative Solvers, with applications in incompressible fluid dynamics*, volume 21 of *Numerical Mathematics and Scientific Computation*. Oxford University Press, ii edition, 2014.

[9] Howard C. Elman and Minghao Wu. Lyapunov inverse iteration for computing a few rightmost eigenvalues of large generalized eigenvalue problems. *SIAM J. Matrix Anal. Appl.*, 34(4):1685–1707, 2013.

[10] V. Faber, T. A. Manteuffel, and S. V. Parter. On the theory of equivalent operators and application to the numerical solution of ununiform elliptic partial differential equations. *Advances in Applied Mathematics*, 11:109–163, 1990.

[11] Tomáš Gergelits, Kent-André Mardal, Bjørn Fredrik Nielsen, and Zdeněk Strakoš. Laplacian preconditioning of elliptic PDEs: Localization of the eigenvalues of the discretized operator. *SIAM Journal on Numerical Analysis*, 57(3):1369–1394, 2019.

[12] G. Golub and C. F. Van Loan. *Matrix Computations*. The Johns Hopkins University Press, Baltimore, 4th edition, 2013.

[13] I. G. Graham and M. J. Hagger. Unstructured additive Schwarz-conjugate gradient method for elliptic problems with highly discontinuous coefficients. *SIAM Journal on Scientific Computing*, 20(6):2041–2066, 1999.

[14] J. Guckenheimer, M. Myers, and B. Sturmfels. Computing Hopf bifurcations. *I. SIAM Journal on Numerical Analysis*, 34(1):1–21, 1997.

[15] R. Hiptmair. Operator preconditioning. *J. Computers & Mathematics with Applications*, 52(5):699–706, 2006.

[16] Michiel E. Hochstenbach and Bor Plestenjak. A Jacobi–Davidson type method for a right definite two-parameter eigenvalue problem. *SIAM J. Matrix Anal. Appl.*, 24(2):392–410, 2002.

[17] R. A. Horn and C. R. Johnson. *Topics in Matrix Analysis*. Cambridge University Press, Cambridge, 1991.

[18] R. A. Horn and C. R. Johnson. *Matrix Analysis*. Cambridge University Press, Cambridge, II edition, 2013.

[19] W. Joubert, Th. Manteuffel, S. Parter, and S.-P. Wong. Preconditioning second-order elliptic operators: experiment and theory. *SIAM J. Sci. Stat. Comput.*, 13(1):259–288, 1992.

[20] Boris N. Khoromskij. *Tensor Numerical Methods in Scientific Computing*. De Gruyter, Berlin, Boston, 2018.

[21] L.-H. Lim. Tensors in computations. In *Acta Numerica*, pages 555–764. Oxford Univ. Press, 2021.
[22] T. Manteuffel and J. Otto. Optimal equivalent preconditioners. SIAM J. Numer. Anal., 30(3):790-812, June 1993.

[23] T. A. Manteuffel and S. V. Parter. Preconditioning and boundary conditions. SIAM J. Numer. Anal., 27(3):656-694, 1990.

[24] Kent-Andre Mardal and Ragnar Winther. Preconditioning discretizations of systems of partial differential equations. Numerical Linear Algebra with Applications, 18(1):1-40, 2011.

[25] Karl Meerbergen and Alastair Spence. Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems. SIAM Journal on Matrix Analysis and Applications, 31(4):1982-1999, May 2010.

[26] Karl Meerbergen and Raf Vandebril. A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line. Lin. Alg. Appl., 436:2828-2844, 2012.

[27] Bjørn Fredrik Nielsen, Aslak Tveito, and Wolfgang Hackbusch. Preconditioning by inverting the Laplacian: an analysis of the eigenvalues. IMA Journal of Numerical Analysis, 29:24-42, 2009.

[28] Y. Saad. Iterative methods for sparse linear systems. SIAM, Society for Industrial and Applied Mathematics, 2nd edition, 2003.

[29] S. Serra-Capizzano. Asymptotic results on the spectra of block Toeplitz preconditioned matrices. SIAM J. Matrix Anal. Appl., 20(1):31-44, 1998.

[30] V. Simoncini. Computational methods for linear matrix equations. SIAM Rev., 58(3):377-441, Sept 2016.

[31] Weiwei Sun, Weizhang Huang, and Robert Russell. Finite difference preconditioning for solving orthogonal collocation equations for boundary value problems. SIAM Journal on Numerical Analysis, 33(6):2268-2285, 1996.

[32] A. van der Sluis and H. A. van der Vorst. The rate of convergence of conjugate gradients. Numer. Math., 48:543-560, 1986.

[33] Yunrong Zhu. Domain decomposition preconditioners for elliptic equations with jump coefficients. Num. Lin. Alg. with Appl., 15(2-3):271-289, 2008.

**Appendix A.** In this section we derive mesh independence of the pair $(A, M)$ when the problem is equipped with all Dirichlet boundary conditions.

In the one-dimensional problem, finite differences applied to the Laplace operator yield the tridiagonal matrix $M = \text{tridiag}(-1, 2, -1)$. The algebraic linear system associated to the differential problem (5.3) is $\frac{1}{h^2}Au = f$ or more precisely,

$$ (6.1) \quad \frac{1}{h^2} \begin{bmatrix} a_1' + a_2' & -a_2' & & & \\ -a_2' & \ddots & \ddots & & \\ & \ddots & \ddots & -a_{n'} & \\ & & -a_{n'} & a_{n'} + a_{(n+1)'} \end{bmatrix} \begin{bmatrix} u_1 \\ \vdots \\ \vdots \\ u_n \end{bmatrix} = \begin{bmatrix} f_1 + a_{1'} \frac{u_\beta}{h^2} \\ f_2 \\ \vdots \\ f_{n-1} \\ f_n + a_{(n+1)'} \frac{u_\beta}{h^2} \end{bmatrix}, $$

where we recall that $a_{i'} = a(x_{i-\frac{1}{2}})$. We next prove spectral equivalence in the one-dimensional case. To this end, recalling the definition of $D_0$ in (5.2), we first define

$$ D := \begin{bmatrix} 1 \\ -1 \\ \vdots \\ \vdots \\ 1 \\ -1 \end{bmatrix} \in \mathbb{R}^{(n+1) \times n}, \quad \text{with} \quad D = \begin{bmatrix} D_0 \\ -e_n^T \end{bmatrix} = \begin{bmatrix} I_n \\ -\mathbf{1}^T \end{bmatrix} D_0. $$

**PROPOSITION 6.1.** Let $a(x) \in [\alpha, \beta]$, $x \in [0, 1]$ and $S := \text{diag}(a_1', \dots, a_{(n+1)'}) \in \mathbb{R}^{(n+1) \times (n+1)}$. Let $M, A$ be the matrices defined above. Then $M = D^T D$ and $A = D^T S D$, moreover, $A$ and $M$ are spectrally equivalent, that is

$$ \alpha \le \frac{x^T A x}{x^T M x} \le \beta, \quad \forall x \in \mathbb{R}^n, x \ne 0. $$
Proof. The first equalities for M and A can be proved by examining the matrices elementwise. For the Rayleigh quotient we have

$$ \frac{x^T A x}{x^T M x} = \frac{x^T D^T S D x}{x^T D^T D x} \underset{y=Dx}{=} \frac{y^T S y}{y^T y}, \quad y \in \mathbb{R}^{n+1}. $$

The two bounds follow from recalling that $\alpha y^T y \le y^T S y \le \beta y^T y$ for all $y \in \mathbb{R}^{n+1}$, as $a(x) \in [\alpha, \beta]$ for $x \in [0, 1]$. □

Setting $S = \text{blkdiag}(S_0, a_{(n+1)'})$, we also notice that the eigenproblem $D^T S D x = \lambda D^T D x$ can be written in the following simplified form

$$ D_0^T [I_n, -\mathbf{1}] \begin{bmatrix} S_0 & \\ & a_{(n+1)'} \end{bmatrix} \underbrace{[\mathbf{1}_n \quad -\mathbf{1}^T]}_z = \lambda D_0^T [I_n, -\mathbf{1}] \underbrace{[\mathbf{1}_n \quad -\mathbf{1}^T]}_z D_0 x, $$

that is

$$ (6.2) \qquad (S_0 + a_{(n+1)'} \mathbf{1} \mathbf{1}^T) z = \lambda (I + \mathbf{1} \mathbf{1}^T) z. $$

This form allows us to get more insight into the spectral distribution of the pair $(A, M)$ when the function $a(x)$ is, e.g., constant or piecewise constant in most of the interval interior.

**PROPOSITION 6.2.** i) If $S_0 := \alpha I_n$, then $\text{spec}(A, M) = \{\alpha, (\alpha+na_{(n+1)'})/(1+n)\}$; ii) If $S_0 := \text{blkdiag}(\alpha I_{n_1}, \beta I_{n_2}) \in \mathbb{R}^{n \times n}$, then $\text{spec}(A, M) = \{\alpha, \lambda_1, \lambda_2, \beta\}$, with $\lambda_{1,2} = \frac{-\eta_1 \pm \sqrt{\eta_1^2 - 4\eta_2}}{2}$, where

$$ \eta_1 = - \frac{n a_{(n+1)'} + (1+n_1)\beta + (1+n_2)\alpha}{n+1}, \qquad \eta_2 = \frac{\alpha\beta + a_{(n+1)'}\beta n_1 + a_{(n+1)'}\alpha n_2}{n+1} $$

*Proof.* i) Using the equation (6.2) we obtain $(\alpha I + a_{(n+1)'} \mathbf{1} \mathbf{1}^T) z = \lambda (I + \mathbf{1} \mathbf{1}^T) z$. There exist $n-1$ linearly independent vectors $z_i \perp \mathbf{1}$, so that $(\alpha, z_i)_{i=1,...,n-1}$ are eigenpairs of $(A, M)$.

Moreover, $z = \mathbf{1}$ is an eigenvector with eigenvalue $\lambda = (\alpha + na_{(n+1)'})/(1+n)$.

ii) Let $z = [z^{(1)}; z^{(2)}]$. Using (6.2) we obtain

$$ \begin{bmatrix} \alpha I_{n_1} + a_{(n+1)'} \mathbf{1}_{n_1} \mathbf{1}_{n_1}^T & a_{(n+1)'} \mathbf{1}_{n_1} \mathbf{1}_{n_2}^T \\ a_{(n+1)'} \mathbf{1}_{n_2} \mathbf{1}_{n_1}^T & \beta I_{n_2} + a_{(n+1)'} \mathbf{1}_{n_2} \mathbf{1}_{n_2}^T \end{bmatrix} \begin{bmatrix} z^{(1)} \\ z^{(2)} \end{bmatrix} = \lambda \begin{bmatrix} I_{n_1} + \mathbf{1}_{n_1} \mathbf{1}_{n_1}^T & \mathbf{1}_{n_1} \mathbf{1}_{n_2}^T \\ \mathbf{1}_{n_2} \mathbf{1}_{n_1}^T & I_{n_2} + \mathbf{1}_{n_2} \mathbf{1}_{n_2}^T \end{bmatrix} \begin{bmatrix} z^{(1)} \\ z^{(2)} \end{bmatrix}. $$

For $z^{(1)} \perp \mathbf{1}_{n_1}$, and $z^{(2)} = \mathbf{0}_{n_2}$, we are back to case (i) and $(\alpha, z_i)_{i=1,...,n_1-1}$ are $n_1-1$ eigenpairs. The same can be done for $\beta$, taking $z^{(2)} \perp \mathbf{1}_{n_2}$, $z^{(1)} = \mathbf{0}_{n_1}$, yielding $n_2-1$ more eigenpairs $(\beta, z_i)_{i=1,...,n_1-1}$. The two missing eigenpairs can be obtained as $z^{(1)} = \gamma_1 \mathbf{1}_{n_1}$ and $z^{(2)} = \gamma_2 \mathbf{1}_{n_2}$ as follows. Explicit rewriting allows one to express the eigenproblem in $2 \times 2$ form as follows

$$ \begin{bmatrix} \alpha + n_1 a_{(n+1)'} & n_2 a_{(n+1)'} \\ n_1 a_{(n+1)'} & \beta + n_2 a_{(n+1)'} \end{bmatrix} \begin{bmatrix} \gamma_1 \\ \gamma_2 \end{bmatrix} = \lambda \begin{bmatrix} 1 + n_1 & n_2 \\ n_1 & 1 + n_2 \end{bmatrix} \begin{bmatrix} \gamma_1 \\ \gamma_2 \end{bmatrix}. $$

Computing the roots of the associated characteristic polynomial of degree two, yields $\lambda_{1,2} = \frac{1}{2}(-\eta_1 \pm \sqrt{\eta_1^2 - 4\eta_2})$. □

For $a_{(n+1)'} = \beta$, the quantities $\eta_1, \eta_2$ nicely simplify, showing linear combinations of $\alpha$ and $\beta$, that is $\lambda_1 = \beta, \lambda_2 = \beta n_1/(n+1) + \alpha(n_2+1)/(n+1)$.
In the two-dimensional case, we have the two matrices $\mathbf{A} = I \otimes A + A \otimes I$ and $\mathbf{M} = I \otimes M + M \otimes I$ with $A$ and $M$ as defined above. The problem $\mathbf{A}\mathbf{x} = \lambda\mathbf{M}\mathbf{x}$ can then be written as

$$ (6.3) \qquad \begin{aligned} & (D_0^T S_0 D_0 + a_{(n+1)'}) e_n e_n^T) X + X (D_0^T S_0 D_0 + a_{(n+1)'}) e_n e_n^T) \\ & = \lambda ((D_0^T D_0 + e_n e_n^T) X + X (D_0^T D_0 + e_n e_n^T)). \end{aligned} $$

We notice that vectors $u \perp e_n$ generate eigenvectors in the form $X = uu^T$ that are also eigenvectors of the problem with mixed boundary conditions. Hence, in case of piecewise constant coefficients the same properties of the eigenpairs of this latter problem can be established. We can extend Theorem 5.4 to the problem (6.3). We first define

$$ P_3^{(0)} = \text{blkdiag}(0_{n_1-1}, P_3) \in \mathbb{R}^{n \times n}, \quad P_3 = D_0 D_0^T \in \mathbb{R}^{(n_2+1) \times (n_2+1)} $$

(as opposed to $P_2$, the matrix $P_3$ takes the value one in the (1,1) position, and the value two in the $(n,n)$ position); we note that $P_3^{(0)}$ is symmetric and positive semidefinite. Then we define

$$ \mathbf{P}_3 = I_n \otimes P_3^{(0)} + P_3^{(0)} \otimes I_n, $$

**THEOREM 6.3.** Let $\theta_j, j = 1, \dots, k$ be the nonnegative real eigenvalues of the pencil $(\mathbf{P}_3, \mathbf{P}_1)$, that is $\theta_j \neq \infty$. Then the eigenvalues $\lambda$ of $(\mathbf{A}, \mathbf{M})$ satisfy

$$ \lambda \in \left\{ \beta, \frac{\beta\theta_1 + \alpha}{\theta_1 + 1}, \dots, \frac{\beta\theta_k + \alpha}{\theta_k + 1} \right\}. $$

The proof follows the same steps as that for Theorem 5.4, hence it is omitted.

**Appendix B.** We report the proof of Proposition 5.5. The proof could be derived using general arguments for the regular pair $(\mathbf{P}_2, \mathbf{P}_1)$, when both matrices are singular. We prefer to report a proof that explicitly uses the structure derived in Proposition 4.1.

*Proof.* Let $(\theta, \mathbf{x})$ with $\|\mathbf{x}\| = 1$ and $\theta \in \mathbb{R}$, be an eigenpair of $(\mathbf{P}_2, \mathbf{P}_1)$, so that $\mathbf{P}_2\mathbf{x} = \theta\mathbf{P}_1\mathbf{x}$. Thanks to Proposition 4.1, $\mathbf{x} = \text{vec}(X)$ where either $X = z_i\gamma_1 z_j^T + z_j\gamma_2 z_i^T$ with $z_i, z_j$ unit norm eigenvectors of $(P_2^{(0)} - \theta P_1^{(0)})$ orthogonal to each other, or $X = \mathit{zz}^T$ with $(\theta, z)$ an eigenpair of $(P_2^{(0)}, P_1^{(0)})$ with $\theta \neq 0$. In the latter case, the bound readily follows. Let us thus focus on $X = z_i\gamma_1 z_j^T + z_j\gamma_2 z_i^T$, with $\gamma_1^2 + \gamma_2^2 = 1$. Multiplying by $\mathbf{x}^T$ from the left the eigenequation, we obtain $\mathbf{x}^T \mathbf{P}_2 \mathbf{x} = \theta \mathbf{x}^T \mathbf{P}_1 \mathbf{x}$, which in matrix form reads

$$ (6.4) \qquad \text{trace}(X^T P_2^{(0)} X + X P_2^{(0)} X^T) = \theta \text{trace}(X^T P_1^{(0)} X + X P_1^{(0)} X^T). $$

For $k = 1, 2$, using $z_i^T z_j = 0$, we have

$$ \text{trace}(X^T P_k^{(0)} X) = \gamma_1^2 \|z_j\|^2 z_i^T P_k^{(0)} z_i + \gamma_2^2 \|z_i\|^2 z_j^T P_k^{(0)} z_j, $$

$$ \text{trace}(X P_k^{(0)} X^T) = \gamma_1^2 \|z_i\|^2 z_j^T P_k^{(0)} z_j + \gamma_2^2 \|z_j\|^2 z_i^T P_k^{(0)} z_i. $$

Hence,

$$ (6.5) \qquad \begin{aligned} \text{trace}(X^T P_k^{(0)} X) + \text{trace}(X P_k^{(0)} X^T) &= \|z_j\|^2 z_i^T P_k^{(0)} z_i + \|z_i\|^2 z_j^T P_k^{(0)} z_j \\ &\leq 2 \|z_i\|^2 \|z_j\|^2 \lambda_{\max}(P_k^{(0)}). \end{aligned} $$
For the lower bound, we recall that $z_i$, say, could be the eigenvector associated with the zero eigenvalue of $P_k^{(0)}$, so that $P_k^{(0)}z_i = 0$, $k = 1, 2$. Hence,

$$ \mathrm{trace}(X^T P_k^{(0)} X) + \mathrm{trace}(X P_k^{(0)} X^T) \geq \|z_i\|^2 \|z_j\|^2 \lambda_{\min}(P_k^{(0)}), $$

Summarizing, by plugging the bounds in (6.4),

$$ \lambda_{\min}(P_2^{(0)}) \|z_i\|^2 \|z_j\|^2 \leq \theta 2 \lambda_{\max}(P_1^{(0)}) \|z_i\|^2 \|z_j\|^2 $$

and

$$ 2 \lambda_{\max}(P_2^{(0)}) \|z_i\|^2 \|z_j\|^2 \geq \theta \lambda_{\min}(P_1^{(0)}) \|z_i\|^2 \|z_j\|^2 $$

from which the result follows. □