Valeria Simoncini
email: valeria.simoncini@unibo.it

Daniele Toni
email: daniele.toni2@studio.unibo.it

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

Keywords: Matrix eigenproblem. Lyapunov equation. Structural properties. Eigenvalue distribution AMS subject classifications. 15A22, 15A24, 65F15, 65F45

We are interested in generalized matrix eigenvalue problems of the type AX + XA = λM XM and AX + XA = λ(M X + 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.

1. Introduction. We consider the following matrix eigenvalue problems (1.1) AX + XA = λM XM and (1.2)

AX + XA = λ(M X + XM ),
with A, M ∈ R n×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 [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF], [START_REF] Meerbergen | A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line[END_REF], 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 [START_REF] Elman | Lyapunov inverse iteration for computing a few rightmost eigenvalues of large generalized eigenvalue problems[END_REF], [START_REF] Meerbergen | A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line[END_REF] and classically studied in the context of multiparameter eigenproblems; see, e.g., [START_REF] Atkinson | Multiparameter Eigenvalue Problems[END_REF]Ch. 7], [START_REF] Michiel | A Jacobi-Davidson type method for a right definite two-parameter eigenvalue problem[END_REF]. 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) (A ⊗ I n ⊗ I n + I n ⊗ A ⊗ I n + I n ⊗ I n ⊗ A)x = λ(M ⊗ M ⊗ M )x,
and

(A ⊗ I n ⊗ I n + I n ⊗ A ⊗ I n + I n ⊗ I n ⊗ A)x = (1.4) λ(M ⊗ I n ⊗ I n + I n ⊗ M ⊗ I n + I n ⊗ I n ⊗ M )x,
where I n is the identity matrix of size n, and ⊗ stands for the Kronecker product, which for matrices A ∈ R n A ×m A , A = (a ij ) i=1,...,n A ,j=1,...,m A , and B ∈ R n B ×m B is defined as ( [START_REF] Horn | Topics in Matrix Analysis[END_REF])

A ⊗ B =      a 11 B a 12 B • • • a 1m H B a 21 B a 22 B • • • a 2m H B . . . . . . a n H 1 B a n H 2 B • • • a n H m H B      ∈ R n A n B ×m A m B . (1.5)
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) AX + XA = λX.

Letting (θ i , z i ), i = 1, . . . , n be the eigenpairs of A, then the eigenpairs of (1.6) are obtained as λ = θ i + θ j , X = z i z T j for all possible i, j = 1, . . . , n [START_REF] Horn | Topics in Matrix Analysis[END_REF]Th. 4.4.5]. Analogously, the following tensor Lyapunov eigenvalue problem

(1.7) (A ⊗ I n ⊗ I n + I n ⊗ A ⊗ I n + I n ⊗ I n ⊗ A)x = λx
is equipped with eigenvalues λ = θ i + θ j + θ l , for all possible i, j, l = 1, . . . , n, and corresponding eigenvectors x = z i ⊗ z j ⊗ z l ; see, e.g., [START_REF] Lim | Tensors in computations[END_REF]. 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 [START_REF] Faber | On the theory of equivalent operators and application to the numerical solution of ununiform elliptic partial differential equations[END_REF] 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 X ∈ R n1×n2×n3 . We let 1 n be the vector of all ones with n components; the subscript will be omitted when clear from the context. We define 0 n similarly. The matrix A = diag(a) is a diagonal matrix with the components of the vector a on its diagonal, and A = blkdiag(B, C) is a block diagonal matrix having the square matrices B, C as diagonal blocks.

We also recall that for x = vec(X) -where "vec" stacks all matrix columns one below the other -and conforming dimensions, it holds that x T (I ⊗ M )x = trace(X T M X) and x T (M ⊗ I)x = trace(XM 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 Kronecker operator, the problems (1.1) and (1.2) can be written as vector generalized eigenproblems, namely

(A ⊗ I + I ⊗ A)x = λ(M ⊗ M )x, (2.1) (A ⊗ I + I ⊗ A)x = λ(M ⊗ I + I ⊗ M )x, (2.2)
with x = 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 = M xλ with A, M ∈ R n d ×n d , d = 2, 3, symmetric and positive definite. Hence, a full set of eigenvectors {x 1 , . . . , x n d } can be determined, satisfying x T i M x j = 0 for i ̸ = 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 [START_REF] Elman | Lyapunov inverse iteration for computing a few rightmost eigenvalues of large generalized eigenvalue problems[END_REF], [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF]. 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 eigenproblem associated with (2.1), thus possibly representing the major similarity between the two problems. Indeed, multiplying by (M

1 2 ⊗M 1 2
) -1 from the left both equations, and letting

A = M -1 2 AM -1 2 we obtain ( A ⊗ M -1 + M -1 ⊗ A)y = λy, (2.3) ( A ⊗ M -1 + M -1 ⊗ A)y = λ(I ⊗ M -1 + M -1 ⊗ I)y, (2.4) with y = (M 1 2 ⊗ M 1 2 )
x. This connection and in particular the form in (2.4) will lead to new developments in section 5.

On the Kronecker pair

(A ⊗ I + I ⊗ A, M ⊗ M ).
With appropriate scaling, the generalized eigenvectors in (2.1) satisfy:

(3.1) x T j (M ⊗ M )x i = δ i,j , x T j (A ⊗ I + I ⊗ A)x i = λ i δ i,j ,
where δ 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 ⟨X, Y ⟩ = 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, . . . , (n 2 + n)/2, while the remaining (n 2 -n)/2 matrices are skew-symmetric;

ii) Orthogonality: ⟨(M

1 2 X i M 1 2 ), (M 1 2 X j M 1 2 )⟩ = δ ij and λ i δ ij = 2 ⟨X i , AX j ⟩, i, j = 1, . . . , n 2 .
Proof. We observe that if X i satisfies (1.1) then also X T i satisfies the equation for the same λ: it suffices to transpose the whole equation and use the symmetry of A and M . If it holds that X T i = ±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) follows. If

X T i ̸ = ±X i , then we let V 1 = (X i + X T i ) symmetric and V 2 = (X i -X T i ) skew-symmetric. Then range([vec(V 1 ), vec(V 2 )]) = range([vec(X i ), vec(X T i )]
, that is vec(V 1 ), vec(V 2 ) are again eigenvectors and they span the same (vector) eigenspace as vec(X i ), vec(X T i ). The result in (i) thus also follows. The results in (ii) directly follow from the matricization of the orthogonality properties of the eigenvectors x i 's.

The proof of the previous result shows that if X i is neither symmetric nor skewsymmetric, then the corresponding eigenvalue must have multiplicity (at least) two, and that the symmetric, skew-symmetric matrices (X i + X T i ), (X i -X T i ) span the same eigenspace. In particular, simple eigenvalues only admit symmetric or skewsymmetric 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 {λ i }, {X i } i=1,...,n 2 be the eigenvalues and eigenmatrices associated with the pair (A⊗I +I ⊗A, M ⊗M ). The solution to AY +Y A+M Y M = C with C symmetric can be written as

Y = n 2 i=1 X i sym α i X i , with α i = trace(X T i C) 1 + λ i . Proof. Since {x i } with x i = vec(X i ) are linearly independent, we can write Y = n 2 i=1 α i X i . Hence, C = n 2 i=1 α i (AX i + X i A + M X i M ) = n 2 i=1 α i (1 + λ i )M X i M.
Since trace(X T j (M X i M )) = δ i,j , it follows that trace(X T i C) = α i (1 + λ i ), from which the result follows. Finally, for C symmetric, we can write

C = n k=1 c k η k c T k , so that trace(X T i C) = n k=1 η k trace(c T k X T i 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(θ 1 , . . . , θ n ) and Γ = diag(γ 1 , . . . , γ n ), respectively. Then the eigenpairs have the form (λ k , QE k Q T ) with λ k = (θ i + θ j )/(γ i γ j ) and E k = e i e T j , for i, j = 1, . . . , n, k = (i -1)n + j. This characterization immediately follows from substituting A = QΘQ T and M = QΓQ T in (1.1).

Remark 3.3. The related problem

(3.2) AXM + M XA = 2M XM λ, X ∈ R n×n
has been thoroughly analyzed in [START_REF] Atkinson | Multiparameter Eigenvalue Problems[END_REF] and more recently, e.g., in [START_REF] Meerbergen | A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line[END_REF]; see also [START_REF] Elman | Lyapunov inverse iteration for computing a few rightmost eigenvalues of large generalized eigenvalue problems[END_REF] for a direct application to stability analysis, where A = I is used. In [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF], [START_REF] Meerbergen | A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line[END_REF] 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 AX + X A T = 2λX, A = M -1 A.

Since both M and A are symmetric and positive definite, 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 {θ i } i=1,...,n . Simple eigenvalues have X = z i y T j as corresponding eigenmatrices with eigenvalues λ = (θ i +θ j )/2, for i, j = 1, . . . , n. Due to the structure of 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 [START_REF] Meerbergen | A reflection on the implicitly restarted arnoldi method for computing eigenvalues near a vertical line[END_REF] 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 {x i } i=1,...,n 3 , and their tensor form

{X i } i=1,...,n 3 . Let the corresponding mode-m (with m ∈ {1, 2, 3}) matricization be [X (i) 1 , X (i) 2 , . . . , X (i) n ]. Then the matrices X (i) k satisfy: i) Symmetry: X (i)
k is either symmetric or skew-symmetric; ii) Orthogonality:

n k=1 n ℓ=1 M k,ℓ trace(M X (j) k M X (i) ℓ ) = δ i,j .
where trace(M

X (j) k M X (i) ℓ ) = ⟨(M 1 2 X (j) k M 1 2 ), (M 1 2 X (i) ℓ M 1 2 )⟩.
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 ⊗ I n + I n ⊗ A) = λM X(M ⊗ M ).
For each block of columns X j we can write

AX j + X j A + n k=1 A j,k X k = λM n k=1 M j,k X k M.
By transposing both equation sides, we observe that both X ℓ and X T ℓ , ℓ = 1, . . . , n satisfy the matrix equation. If X ℓ is symmetric or skew-symmetric for all ℓ = 1, . . . , n, then we are done. If X ℓ ̸ = X T ℓ then the matrices V ℓ := X ℓ + X T ℓ and W ℓ := X ℓ -X T ℓ , ℓ = 1, . . . , n also satisfy the matrix equation and their vectorization spans the same space as the vectorization of X ℓ , X T ℓ . ii) We only have to unfold the first vector property in (3.1). Let

X i = [X (i) 1 , X (i) 2 , . . . , X (i) n ] with X (i) k = X i (:, :, k) and X j = [X (j) 1 , X (j) 2 , . . . , X (j) n ] with X (j) k = X j (:, :, k) be the matricizations of the eigenvectors x i , x j , respec- tively, such that x T i (M ⊗ M ⊗ M )x j = δ i,j .
Unfolding the Kronecker product we obtain

δ i,j = x T i (M ⊗ M ⊗ M )x j = trace(M X j (M ⊗ M ) X T i ) = trace(M n k=1 (X (j) k ( n ℓ=1 M k,ℓ M X (i) ℓ ))) = trace(M n k=1 (X (j) k M ( n ℓ=1 M k,ℓ X (i) ℓ ))).
Reordering terms yields the desired relation.

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 [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF]. In there, the problem

(4.1) AXB + BXA = λ(BXN + N XB), X ∈ R n×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 [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF] can be readily employed in our setting. Indeed, let us collect all terms in (1.2) on the left-hand side as

(4.2) (A -λM )X + X(A -λM ) = 0,
and let G(λ) := A-λM . The original eigenproblem can be related to the solvability of the matrix equation G(λ)X + XG(λ) = 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

(ζ i , z i ), (ζ j , z j ) of G(λ k ), with i ̸ = j, such that ζ i = -ζ j and X k = [z i , z j ]Γ[z j , z i ] T = z i γ 1 z T j + z j γ 2 z T i , for any nonzero Γ = diag([γ 1 , γ 2 ]); ii) There exists an eigenpair (ζ, z) of G(λ k ) with simple ζ = 0 such that X k = zz T . Proof. For a given λ k , a solution to G(λ k )X + XG(λ k ) = 0 is determined in the null space of G(λ k ) = I ⊗ G(λ k ) + G(λ k ) ⊗ I, that is, λ k should be such that G(λ k ) is singular.
The eigenvalues of G(λ k ) are given by ζ i + ζ j , where ζ's are the eigenvalues of G(λ k ), for all combinations of i, j ∈ {1, . . . , n}. Hence, G(λ k ) is singular if and only if there exist two eigenvalues such that ζ i + ζ j = 0. The rest of the result follows by substitution.

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 characterization of Proposition 4.1 motivated the authors of [START_REF] Meerbergen | Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large scale problems[END_REF] to use this matrix equation in the so-called bialternate product method for detecting Hopf bifurcations [START_REF] Guckenheimer | Computing Hopf bifurcations I[END_REF].

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

(A ⊗ I n ⊗ I n + I n ⊗ A ⊗ I n + I n ⊗ I n ⊗ A)x = (4.3) λ(M ⊗ I n ⊗ I n + I n ⊗ M ⊗ I n + I n ⊗ I n ⊗ M )x,
by exploiting the same matrix G(λ) defined above, so as to write

(G(λ) ⊗ I n ⊗ I n + I n ⊗ G(λ) ⊗ I n + I n ⊗ I n ⊗ G(λ))x = 0, or G(λ)x = 0 in short.
Proposition 4.2. With the notation above, we have that (λ k , x k ) is an eigenpair of (4.3) if and only if one of the following holds i) There exist three eigenpairs (θ i , z i ), (θ j , z j ), (θ l , z l ) of G(λ k ) with i, j, l not equal, such that θ i + θ j + θ l = 0 and

x k = γ 1 z i ⊗ z j ⊗ z l + γ 2 z i ⊗ z l ⊗ z j + γ 3 z j ⊗ z i ⊗ z l + γ 4 z j ⊗z l ⊗z i +γ 5 z l ⊗z i ⊗z j +γ 6 z l ⊗z j ⊗z i , for any nonzero Γ = diag(γ 1 , γ 2 , γ 3 , γ 4 , γ 5 , γ 6 );
ii) There exists an eigenpair

(θ i , z i ) of G(λ k ) with simple θ i = 0 so that x k = z i ⊗ z i ⊗ z i .
Proof. For a given λ k , a solution x to G(λ k )x = 0 is determined in the null space of the coefficient matrix G(λ k ). The eigenvalues of G(λ k ) are given by θ i + θ j + θ l , where θ's are the eigenvalues of G(λ k ), for all combinations of i, j, l. Hence, G(λ k ) is singular if and only if there exist three eigenvalues such that θ i + θ j + θ l = 0. The rest of the result follows by substitution.

We remark that the non-nullity of the γ i s in Proposition 4.2 is related to the multiplicity of the eigenvalues θ's in G(λ k ). In other words, each x k is composed by up to six terms, with the maximum occurring when the θ's are such that {z i , z j , z l } spans a three dimensional subspace.

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) u x (1, y) = 0, y ∈ [0, 1], u y (x, 1) = 0, x ∈ [0, 1], u(x, y) = 0 elsewhere on ∂Ω.

To determine the finite difference discretization of the PDE, let us consider the gridpoints (x i , y j ), i = 1, . . . , n of a uniform discretization of Ω. 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) D 0 :=       1 -1 . . . . . . . . . -1 1       ∈ R n×n , S a,0 := diag(a i ′ ) i=1,...n ∈ R n×n , S b,0 := diag(b i ′ ) i=1,...n ∈ R n×n .
where a i ′ = a(x i-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) A a U + U A b = F, A k = D T 0 S k,0 D 0 , k = a, b,
where F = f (x i , y j ) + b.c. and U ≈ u(x i , y j ). Note that the Poisson equation, corresponding to constant a(x) = 1, b(y) = 1 is given by1

(5.4)

M U + U M = F, M = D T 0 D 0 . Classically, equation (5.
3) is vectorized into the equation Au = f , with u = vec(U ), and f = vec(F ), and A = I ⊗ A a + A b ⊗ I; see, e.g., [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF] 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 → M U + U M on the left-hand side of (5.4), with the same boundary conditions [START_REF] Manteuffel | Preconditioning and boundary conditions[END_REF]. 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., [START_REF] Faber | On the theory of equivalent operators and application to the numerical solution of ununiform elliptic partial differential equations[END_REF], [START_REF] Bjørn Fredrik Nielsen | Preconditioning by inverting the Laplacian: an analysis of the eigenvalues[END_REF], and [START_REF] Manteuffel | Optimal equivalent preconditioners[END_REF] for a spectral analysis of this strategy on convection-diffusion operators and [START_REF] Sun | Finite difference preconditioning for solving orthogonal collocation equations for boundary value problems[END_REF] 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 [START_REF] Axelsson | Mesh independent superlinear PCG rates via compact-equivalent operators[END_REF], [START_REF] Faber | On the theory of equivalent operators and application to the numerical solution of ununiform elliptic partial differential equations[END_REF], [START_REF] Hiptmair | Operator preconditioning[END_REF], [START_REF] Manteuffel | Optimal equivalent preconditioners[END_REF]. In particular, in the recent work [START_REF] Gergelits | Laplacian preconditioning of elliptic PDEs: Localization of the eigenvalues of the discretized operator[END_REF] the authors have proved what they called a pairing between the values of operator coefficient function κ 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.

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 [START_REF] Manteuffel | Preconditioning and boundary conditions[END_REF]. 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 -1 2 ⊗ M -1 2 yields a convenient form for our analysis. In matrix form, this multiplication transforms (2.4) into

AXM -1 + M -1 X A = λ(XM -1 + M -1 X) with A = M -1 2 AM -1 2 .
If M and A are spectrally equivalent, then A is spectrally equivalent to the identity matrix, which leads to conclude that the eigenvalues λs do not depend on the problem size, but only on the constants yielding spectral equivalence of A and M .

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 [START_REF] Faber | On the theory of equivalent operators and application to the numerical solution of ununiform elliptic partial differential equations[END_REF], [START_REF] Joubert | Preconditioning second-order elliptic operators: experiment and theory[END_REF], [START_REF] Manteuffel | Optimal equivalent preconditioners[END_REF], [START_REF] Manteuffel | Preconditioning and boundary conditions[END_REF].

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 (κ 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., [START_REF] Howard | Finite Elements and Fast Iterative Solvers, with applications in incompressible fluid[END_REF], [START_REF] Mardal | Preconditioning discretizations of systems of partial differential equations[END_REF], [START_REF] Hiptmair | Operator preconditioning[END_REF]. Hence, in the one-dimensional case for instance,

∇ • κ∇ → D T 0 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 ̸ = a.

In the one-dimensional case with mixed boundary conditions, the spectrum of the pair (D T 0 S a,0 D 0 , D T 0 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 α ≤ a(x) ≤ β for all x ∈ [0, 1]. Then A and M are spectrally equivalent, and it holds that

α ≤ x T Ax x T M x ≤ β, ∀0 ̸ = x ∈ R n 2 .
Proof. We first note that from a(x) ∈ [α, β] ∀x ∈ [0, 1], it follows that spec(S 0 ) ⊆ [α, β]. Using the spectral properties of the Kronecker product, it holds that

M = D T 0 ⊗ I n , I n ⊗ D T 0 D 0 ⊗ I n I n ⊗ D 0 =: D T 0 D 0 (5.5) A = D T 0 ⊗ I n , I n ⊗ D T 0 S 0 ⊗ I n I n ⊗ S 0 D 0 ⊗ I n I n ⊗ D 0 =: D T 0 S 0 D 0 (5.6)
with S 0 ∈ R n 2 ×n 2 a diagonal matrix. Using the spectral properties of the Kronecker sum, it also follows that spec(S 0 ) ⊆ [α, β]. For 0 ̸ = x ∈ R n 2 we have

x T Ax x T M x = x T D T 0 S 0 D 0 x x T D T 0 D 0 x = 0̸ =y=D0x y T S 0 y y T y ∈ [α, β]. □
We remark that completely analogous formulations can be obtained in the tensor case [START_REF] Khoromskij | Tensor Numerical Methods in Scientific Computing[END_REF], 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 [α, β], 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 ) ̸ = 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

λ of A a X + XA b = λ(M X + XM ) satisfy λ ∈ 2 min x∈(0,1) {a(x), b(x)}, 2 max x∈(0,1) {a(x), b(x)} .
Piecewise constant coefficients. An interesting characterization can be derived for a(x) piecewise constant in [α, β], showing that the eigenvalues do not distribute uniformly in [α, β], as, for instance, empty gaps arise. Though the problem with piecewise constant values could be addressed with more powerful techniques (see, e.g., [START_REF] Serra-Capizzano | Asymptotic results on the spectra of block Toeplitz preconditioned matrices[END_REF]), our approach would also allow one to naturally predict the behavior of continuous coefficients with steep value jumps by means of a perturbation analysiswe 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 [START_REF] Gergelits | Laplacian preconditioning of elliptic PDEs: Localization of the eigenvalues of the discretized operator[END_REF]. 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., [START_REF] Cai | An analysis of a preconditioner for the discretized pressure equation arising in reservoir simulation[END_REF], [START_REF] Dolgov | A reciprocal preconditioner for structured matrices arising from elliptic problems with jump coefficients[END_REF] and their references; in general, algebraic multigrid and domain decomposition preconditioners seem to be robust with respect to jumps, see, e.g., [START_REF] Dryja | Multilevel Schwarz methods for elliptic problems with discontinuous coefficients in three dimensions[END_REF], [START_REF] Blatt | A Parallel Algebraic Multigrid Method for Elliptic Problems with Highly Discontinuous Coefficients[END_REF], [START_REF] Graham | Unstructured additive Schwarz-conjugate gradient method for elliptic problems with highly discontinuous coefficients[END_REF], and [START_REF] Zhu | Domain decomposition preconditioners for elliptic with jump coefficients[END_REF] 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 ∈ (0, 1), that is a(x) = α for x ≤ x 0 and a(x) = β for x > x 0 , with α ̸ = β, so that after discretization, S 0 := blkdiag(αI n1 , βI n2 ), 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 2 1 eigenpairs (α, X), k = 1, . . . , n 2 1 such that (α, z i ) and (α, z j ) are an eigenpair of (A, M ), and X = [z i , z j ]Γ[z j , z i ] T for all possible combinations of i and j, and Γ symmetric, either positive definite or indefinite. In particular, for a fixed eigenvalue λ k , there are two independent eigenvectors, which are

x 1 = vec([z 1 , z 2 ]Γ 1 [z 2 , z 1 ] T ) and x 2 = vec([z 1 , z 2 ]Γ 2 [z 2 , z 1 ] T )
with the two Γs having different signature. Analogously, the pencil has n 2 2 eigenvalues equal to β, with corresponding eigenvectors.

The following theorem shows that all other eigenvalues are a convex linear combination of α and β. To this end, let n 1 , n 2 ∈ N with n 1 + n 2 = n, and define

P (0) 1 = blkdiag(P 1 , 0 n2 ) ∈ R n×n , P 1 = D T 0 D 0 ∈ R n1×n1 and P (0) 2 = blkdiag(0 n1-1 , P 2 ) ∈ R n×n , P 2 =        1 -1 -1 2 -1 . . . . . . . . . -1 2 -1 -1 1        ∈ R (n2+1)×(n2+1) .
We note that both matrices P the two nonzero blocks overlap in the (n 1 , n 1 ) element. Then we define

P 1 = I n ⊗ P (0) 1 + P (0) 1 ⊗ I n , P 2 = I n ⊗ P (0) 2 + P (0) 2 ⊗ 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 θ j , j = 1, . . . , k be the non-negative real eigenvalues of the pencil (P 2 , P 1 ), that is θ j ̸ = ∞. Then the eigenvalues λ of (A, M ) satisfy

λ ∈ β, βθ 1 + α θ 1 + 1 , . . . , βθ k + α θ k + 1 .
Proof. We can write A = αP 1 + βP 2 , M = D T 0 D 0 = P 1 + P 2 . The problem Ax = λM x can be written as (α -λ)P 1 x + (β -λ)P 2 x = 0, that is, for λ ̸ = β,

P 2 x = θP 1 x, with θ = - α -λ β -λ .
Writing down λ in terms of θ yields λ = βθ+α θ+1 . The result also includes the case λ = α for θ = 0, and shows that the sought after linear combination is λ = θ θ+1 β + 1 θ+1 α. Subtracting (α + β)/2 on both sides and rearranging terms, we obtain

λ -α+β 2 β-α 2 = θ -1 θ + 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 θ ≥ 0. The distribution of θ around the value 1 provides information on the distribution of λ around the middle value (α + β)/2.

To proceed further, it is thus useful to obtain a refined analysis on the location of the eigenvalues θ of (P 2 , P 1 ). Proposition 5.5. Let λ max (P i ) be the largest eigenvalue of P i , i = 1, 2, and λ min (P i ) be the smallest nonzero eigenvalue of P i , i = 1, 2. Then the finite nonzero eigenvalues θ of (P 2 , P 1 ) satisfy

1 2 λ min (P 2 ) λ max (P 1 )
≤ θ ≤ 2 λ max (P 2 ) λ 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.

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)

λ min (P 2 ) λ max (P 1 ) = 1 -cos π 2n2+1 1 -cos (2(n-n2)-1)π 2(n-n2)+1 , λ max (P 2 ) λ min (P 1 ) = 1 -cos (2n2-1)π 2n2+1 1 -cos π 2(n-n2)+1
, thus allowing a better localization of the eigenvalues λ distinct from α, β of the pair (A, M ), in the interval [α, β]. In particular, for n 2 ≪ n, a gap can be observed between λ = α (corresponding to θ = 0) and the next distinct eigenvalue, as shown by the lower bound for θ. For n 2 ≈ n, a gap can be observed between λ = β (corresponding to an infinite θ) and the previous distinct eigenvalue, monitored by the upper bound for θ. Finally, our computational experiments have shown that for n 2 ≈ (n -n 2 ) ≈ n/2 a cluster of eigenvalues λ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 [λ min , λ max ] instead of [α, β], where λ min , λ max are the eigenvalues of (A, M ) closest to α, β, respectively, but distinct from them.

n 2 =2 n 2 =4 n 2 =6 n 2 =8 n 2 =10 n 2 =12 n 2 =14 n 2 =16 n 2 =18 n 2 =20 n 2 =22 n 2 =24 n 2 =26 n 2 =28 n 2 =30 n 2 =32 n 2 =34 n 2 =36 n 2 =38 n 2 =40 n 2 =42 n 2 =44 n 2 =46 n 2 =48

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 ) = β in the last n 2 discretization nodes and to α elsewhere. Varying n 2 corresponds to a smaller or larger portion of the domain where the value β 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 α = 4, β = 10 and n = 50 nodes in each direction. Figure 5.1 reports the eigendistribution for n 2 = 2, 4, 6, . . . , 48. Each set of eigenvalues is depicted as a monotonically increasing curve with values from α to β, with increasing values of n 2 from left to right. The eigenvalue gap near β for n 2 ≪ n is clearly visible among the leftmost curves, while the rightmost curves display the gap near α. We also notice that for values n 2 ≈ 50 = n/2 more eigenvalues take values around 7 = (α + β)/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 λ = θ θ+1 β + 1 θ+1 α not equal to β, as n 2 varies. The bound is obtained by using the upper estimate of θ in (5.7) together with Proposition 5.5. The bound appears to be quite sharp. Corresponding similar values can be obtained for α. Example 5.7. We consider again the data from Example 5.6, but with different values of α and β, namely α = 1, β = 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 × n 2 . The shown curves refer to the energy-norm of the error. We also report (dashed line) the expected asymptotic convergence rate ρ k as the iteration k progresses, with

ρ = cond ef f -1 cond ef f + 1 , cond ef f := λ max λ min ;
here cond ef f denotes the "effective" condition number of the preconditioned matrix M -1 A, obtained by excluding the eigenvalues α, β. 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 β/α, much larger than cond ef f , associated with the whole spectral interval. We refer, e.g., to [START_REF] Axelsson | On the rate of convergence of the preconditioned conjugate gradient method[END_REF], [START_REF] Van Der Sluis | The rate of convergence of conjugate gradients[END_REF], [START_REF] Axelsson | Mesh independent superlinear PCG rates via compact-equivalent operators[END_REF], [START_REF] Gergelits | Laplacian preconditioning of elliptic PDEs: Localization of the eigenvalues of the discretized operator[END_REF] 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 ≈ 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. We then compare the eigenvalue distribution for the pair (A, M ) with that obtained with the continuous coefficients

a(x) =    10 4
x < 0.5 10 4 exp(25(0.5 -x)) + 1.1 x ∈ [0.5, 1.1) 1

x > 1.1. x -(bu y ) y -(cu z ) z = f in Ω = (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 case. We show that the twodimensional 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

A = A ⊗ I n ⊗ I n + I n ⊗ A ⊗ I n + I n ⊗ I n ⊗ A ∈ R n 3 ×n 3 M = M ⊗ I n ⊗ I n + I n ⊗ M ⊗ I n + I n ⊗ I n ⊗ M ∈ R n 3 ×n 3
Theorem 5.9. Let a = a(x) be the coefficient function and suppose a(x) ∈ [α, β] for x ∈ (0, 1). Then A e M are spectrally equivalent, and it holds

α ≤ x T Ax x T M 0 x ≤ β, ∀x ∈ R n 3 , x ̸ = 0.
Proof. We define

D 0 :=   D 0 ⊗ In ⊗ In In ⊗ D 0 ⊗ In In ⊗ In ⊗ D 0  , S :=   S 0 ⊗ In ⊗ In In ⊗ S 0 ⊗ In In ⊗ In ⊗ S 0   ,
with D 0 ∈ R 3n 3 ×n 3 and C ∈ R 3n 3 ×3n 3 . Observe that A = D T 0 SD 0 and M = D T 0 D 0 . Hence the proof follows the same lines as that of Theorem 5.2, since spec(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

P i = I n ⊗ I ⊗ P (0) i + I n ⊗ P (0) i ⊗ I n + P (0) i ⊗ I n ⊗ I n , i = 1, 2.
Theorem 5.10. Let θ j , j = 1, . . . , k be the nonnegative real eigenvalues of the pencil (P 2 , P 1 ), that is θ j ̸ = ∞. Then the eigenvalues λ of (A, M ) satisfy

λ ∈ β, βθ 1 + α θ 1 + 1 , . . . , βθ k + α θ k + 1 .
These results can be formulated in any dimension d, with d ≥ 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.

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 = tridiag(-1, 2, -1). The algebraic linear system associated to the differential problem (5.3) is

u n       =        f 1 + a 1 ′ uα h 2 f 2 . . . f n-1 f n + a (n+1) ′ u β h 2       
, where we recall that a i ′ = a(x i- 1 2 ). We next prove spectral equivalence in the onedimensional case. To this end, recalling the definition of D 0 in (5. The two bounds follow from recalling that αy T y ≤ y T Sy ≤ βy T y for all y ∈ R n+1 , as a(x) ∈ [α, β] for x ∈ [0, 1].

Setting S = blkdiag(S 0 , a (n+1) ′ ), we also notice that the eigenproblem D T SDx = λD T Dx can be written in the following simplified form

D T 0 [I n , -1] S 0 a (n+1) ′ I n -1 T D 0 x z = λD T 0 I n -1 I n -1 T D 0 x z ,
that is (6.2) (S 0 + a (n+1) ′ 11 T )z = λ(I + 11 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. Proof. i) Using the equation (6.2) we obtain αI + a (n+1) ′ 11 T z = λ I + 11 T z. There exist n -1 linearly independent vectors z i ⊥ 1, so that (α, z i ) i=1,...,n-1 are eigenpairs of (A, M ).

Moreover, z = 1 is an eigenvector with eigenvalue λ = (α + na (n+1) ′ )/(1 + n).

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

αIn 1 + a (n+1) ′ 1n 1 1 T n 1 a (n+1) ′ 1n 1 1 T n 2 a (n+1) ′ 1n 2 1 T n 1 βIn 2 + a (n+1) ′ 1n 2 1 T n 2 z (1) z (2) = λ In 1 + 1n 1 1 T n 1 1n 1 1 T n 2 1n 2 1 T n 1 In 2 + 1n 2 1 T n 2 z (1)
z (2) .

For z (1) ⊥ 1 n1 , and z (2) = 0 n2 , we are back to case (i) and (α, z i ) i=1,...,n1-1 are n 1 -1 eigenpairs. The same can be done for β, taking z (2) ⊥ 1 n2 , z (1) = 0 n1 , yielding n 2 -1 more eigenpairs (β, z i ) i=1,...,n1-1 . The two missing eigenpairs can be obtained as z (1) = γ 1 1 n1 and z (2) = γ 2 1 n2 as follows. Explicit rewriting allows one to express the eigenproblem in 2 × 2 form as follows

α + n 1 a (n+1) ′ n 2 a (n+1) ′ n 1 a (n+1) ′ β + n 2 a (n+1) ′ γ 1 γ 2 = λ 1 + n 1 n 2 n 1 1 + n 2 γ 1 γ 2 .
Computing the roots of the associated a characteristic polynomial of degree two, yields λ 1,2 = 1 2 (-η 1 ± η 2 1 -4η 2 ). For a (n+1) ′ = β, the quantities η 1 , η 2 nicely simplify, showing linear combinations of α and β, that is λ 1 = β, λ 2 = βn 1 /(n + 1) + α(n 2 + 1)/(n + 1).

5. 1 .

 1 The discretized problem. We consider the following equation-∇ • (κ(x, y)∇u) = f, (x, y) ∈ Ω = (0, 1) 2 ,with either Dirichlet or mixed (Dirichlet plus Neumann) boundary conditions, and κ(x, y) = diag(a(x), b(y)). Consider the following (mixed) conditions:

Fig. 5 . 1 .

 51 Fig. 5.1. Example 5.6, for n = 50, n 2 = 2 : 2 : 48, α = 4, β = 10.

Example 5 . 8 .

 58 We consider the PDE -(a(x)u x ) x -(a(y)u y ) y = f with f (x, y) = sin(xπ) cos(yπ), (x, y) ∈ [0, π] 2 and Dirichlet boundary conditions, with the coefficients choice a(x) = 10 4 x < 0.5 1 x ≥ 0.5.

Figure 5 . 4 (

 54 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.

Fig. 5 . 2 .Fig. 5 . 3 .

 5253 Fig. 5.2. Example 5.6. Convergence of PCG. Left: n 2 = 4. Right: n 2 = 10.

Figure 5 . 4 ( 5 . 4 .

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

Fig. 5 . 4 .

 54 Fig. 5.4. Example 5.8. Left: solution U from u = A -1 f , u = vec(U ). Center: graphs of a(x) and a(x). Right: eigenvalues of the pair (A, M ) for A defined using a(x) or a(x).

Proposition 6 . 1 .

 61 Let a(x) ∈ [α, β], x ∈ [0, 1] and S := diag(a 1 ′ , . . . , a (n+1) ′ ) ∈ R (n+1)×(n+1) . Let M , A be the matrices defined above. Then M = D T D and A = D T SD, moreover, A e M are spectrally equivalent, that isα ≤ x T Ax x T M x ≤ β, ∀x ∈ R n , x ̸ = 0.Proof. The first equalities for M and A can be proved by examining the matrices elementwise. For the Rayleigh quotient we havex T Ax x T M x = x T D T SDx x T D T Dx = y=Dx y T Sy y T y , y ∈ R n+1 .

Proposition 6 . 2 . 1 1 , η 2 =

 62112 i) If S 0 := αI n , then spec(A, M ) = {α, (α+na (n+1) ′ )/(1+n)}; ii) If S 0 := blkdiag(αI n1 , βI n2 ) ∈ R n×n , then spec(A, M ) = {α, λ 1 , λ 2 , β}, with λ na (n+1) ′ + (1 + n 1 )β + (1 + n 2 )α n + αβ + a (n+1) ′ βn 1 + a (n+1) ′ αn 2 n + 1

Table 5 . 1

 51 Example 5.6. Largest eigenvalue distinct from β and its upper bound, for n = 50.

	n 2	λ max λ bound n 2	λ max λ bound
	2 9.9976 9.9988	. . .	. . .	. . .
	6 9.9980 9.9990 34 9.9868 9.9932
	10 9.9977 9.9988 38 9.9772 9.9882
	14 9.9972 9.9986 42 9.9516 9.9745
	18 9.9965 9.9982 46 9.8384 9.9108
	. . .	. . .	. . . 48 9.5400 9.7263

  1 h 2 Au = f or more precisely, ′ -a n ′ a n ′ + a (n+1) ′

	(6.1)	1 h 2	      a 1 ′ + a 2 ′ -a 2 ′ -a 2 ′ . . . . . .	. . . . . .	-a n      	     	u 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.

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

In the two-dimensional case, we have the two matrices A = I ⊗ A + A ⊗ I and M = I ⊗ M + M ⊗ I with A and M as defined above. The problem Ax = λM x can then be written as

We notice that vectors u ⊥ 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

(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 (0) 3 is symmetric and positive semidefinite. Then we define

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 (P 2 , 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 (θ, x) with ∥x∥ = 1 and θ ∈ R, be an eigenpair of (P 2 , P 1 ), so that P 2 x = θP 1 x. Thanks to Proposition 4.1, x = vec(X) where either X = z i γ 1 z T j + z j γ 2 z T i with z i , z j unit norm eigenvectors of (P

1 ) orthogonal to each other, or X = zz T with (θ, z) an eigenpair of (P (0) 2 , P (0) 1 ) with θ ̸ = 0. In the latter case, the bound readily follows. Let us thus focus on X = z i γ 1 z T j + z j γ 2 z T i , with γ 2 1 + γ 2 2 = 1. Multiplying by x T from the left the eigenequation, we obtain x T P 2 x = θx T P 1 x, which in matrix form reads trace(X T P (0)

Hence, trace(X T P (0)

k z j (6.5)

For the lower bound, we recall that z i , say, could be the eigenvector associated with the zero eigenvalue of P (0) k , so that P (0)

Summarizing, by plugging the bounds in (6.4), λ min (P (0)

2 )∥z i ∥ 2 ∥z j ∥ 2 ≤ θ 2 λ max (P (0) 1 )∥z i ∥ 2 ∥z j ∥ 2 and 2 λ max (P (0)

2 )∥z i ∥ 2 ∥z j ∥ 2 ≥ θ λ min (P (0) 1 )∥z i ∥ 2 ∥z j ∥ 2 from which the result follows.