Keywords: , Conditioning and Backward Error, Computation of Invariant Pairs and Matrix Solvents

I am also grateful with the University of Costa Rica, for its excellent academical preparation, for the financial support during my studies in Europe and, also, for reserved me a place of work for my return.

Finally, I want to thank my family, specially to my father Eugenio Segura, for all their effort and support over the years. And the last but

Transpose of a matrix/vector A A *

Conjugate transpose of a matrix A A -1 Inverse of a matrix A A +

Pseudo-inverse of a matrix A σ(A) Spectrum of a square matrix A rank(A) Rank of a matrix A det(A) Determinant of a matrix A tr(A), trace(A)

Trace of a matrix A A ⊗ B Kronecker product of the matrices A and B diag (A 1 , A 2 . . . , A n ) The (block) diagonal matrix

     A 1 0 • • • 0 0 A 2 . . . . . . . . . . . . . . . 0 0 • • • 0 A n      a ij , A(i, j), [A] ij
The (i, j)th entry of a matrix A A(i, :)

The ith row of a matrix A A(:, j)

The jth column of a matrix A x Norm of a vector x (see A.3.1)

Introduction

Computation of Invariant Pairs and Matrix Solvents

Invariant pairs, introduced and analyzed in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], [START_REF] Beyn | Continuation of invariant subspaces for parameterized quadratic eigenvalue problems[END_REF], [START_REF] Effenberger | Robust solution methods for nonlinear eigenvalue problems[END_REF] and [START_REF] Szyld | Several properties of invariant pairs of nonlinear algebraic eigenvalue problems[END_REF], are a generalization of eigenpairs for matrix polynomials. Let P (λ) = ℓ j=0 A j λ j be an n × n matrix polynomial: the polynomial eigenvalue problem consists in computing a scalar λ and a nonzero vector x such that P (λ)x = 0. Now, choose a positive integer k. The matrices X, S of sizes n × k and k × k, respectively, are said to form an invariant pair of size k for P (λ) if X = 0 and: P (X, S) := ℓ j=0 A j XS j = 0.

(1)

Note that the eigenvalues of the matrix S are also eigenvalues of P (λ).

Invariant pairs offer a unified theoretical perspective on the problem of computing several eigenvalue-eigenvector pairs for a given matrix polynomial. From a numerical point of view, moreover, the computation of an invariant pair tends to be more stable than the computation of single eigenpairs, particularly in the case of multiple or tightly clustered eigenvalues. The notion of invariant pairs can also be applied to more general nonlinear problems, although here we will limit our presentation to matrix polynomials.

How to compute invariant pairs? Beyn and Thümmler ( [START_REF] Beyn | Continuation of invariant subspaces for parameterized quadratic eigenvalue problems[END_REF]) adopt a continuation method of predictor-corrector type. Betcke and Kressner ([17]), on the other hand, establish a correspondence between invariant pairs of a given matrix polynomial and of its linearizations. Invariant pairs for P (λ) are extracted from invariant pairs of a linearized form and then refined via Newton's method.

The approach we take in this work to compute invariant pairs is based on contour integrals. Being able to specify the contour Γ allows us to select invariant pairs that have eigenvalues in a prescribed part of the complex plane. Contour integrals play an important role in the definition and computation of moments, which form a Hankel matrix pencil yielding the eigenvalues of the given matrix polynomial that belong to the prescribed contour. The use of Hankel pencils of moment matrices is widespread in several applications such as control theory, signal processing or shape reconstruction, but nonlinear eigenvalue-eigenvector problems can also be tackled through this approach, as suggested for instance in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] and [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF]. E. Polizzi's FEAST algorithm [START_REF] Polizzi | Density-matrix-based algorithm for solving eigenvalue problems[END_REF] is also an interesting example of contour-integral based eigensolver applied to large-scale electronic structure computations.

In this work, we adapt such methods to the computation of invariant pairs. We study, in particular, the scalar moment method and its relation with the multiplicity structure of the eigenvalues, but we also explore the behavior of the block version.

These results on invariant pairs can be applied to the particular case of matrix solvents, that is, to the matrix equation:

P (S) := ℓ j=0
A j S j = 0.

The matrix solvent problem has received remarkable attention in the literature since Sylvester's work [START_REF] Sylvester | On Hamilton's quadratic equation and the general unilateral equation in matrices[END_REF] in the 1880s. The relation between the Riccati and the quadratic matrix equation is highlighted in [START_REF] Bini | Transforming algebraic Riccati equations into unilateral quadratic matrix equations[END_REF] whereas a study on the existence of solvents can be found in [START_REF] Dennis | The algebraic theory of matrix polynomials[END_REF]. Several works address the problem of computing a numerical approximation for the solution of the quadratic matrix equation: an approach to compute, when possible, the dominant solvent is proposed in [START_REF] Dennis | Algorithms for solvents of matrix polynomials[END_REF]. Newton's method and some variations are also used to approximate solvents numerically: see, for example, [START_REF] Davis | Numerical solution of a quadratic matrix equation[END_REF], [START_REF] Kim | Numerical methods for solving a quadratic matrix equation[END_REF], [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF], [START_REF] Long | Improved Newton's method with exact line searches to solve quadratic matrix equation[END_REF]. The work in [START_REF] Hashemi | Efficient computation of enclosures for the exact solvents of a quadratic matrix equation[END_REF] uses interval arithmetic to compute an interval matrix containing the exact solution to the quadratic matrix equation. For the case of the general matrix solvent problem, we can also cite [START_REF] Brás | A spectral approach to polynomial matrices solvents[END_REF], [START_REF] Pereira | On solvents of matrix polynomials[END_REF] and [START_REF] Kratz | Numerical solution of matrix polynomial equations by Newton's method[END_REF].

We exhibit computable formulations for the condition number and backward error of the general matrix solvent problem, thus generalizing previous work on the quadratic matrix equation. Moreover, we propose an adaptation of the moment method to the computation of solvents. Finally, we build on existing work on triangularization of matrix polynomials (see [START_REF] Tisseur | Triangularizing quadratic matrix polynomials[END_REF] and [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF]) and explore the relationship between solvents of matrix polynomials in general and in triangularized form.

This thesis is organized as follows. Chapter 1 introduces preliminary notions, definitions and notation concerning matrix polynomials and their applications, the relation with systems of ordinary differential equations, the generalized eigenvalue problem and some solution methods. The last section discusses eigenvalue shifting for matrix polynomials and presents a generalization of the shifting technique proposed in [START_REF] Meini | A "shift-and-deflate" technique for quadratic matrix polynomials[END_REF]. Eigenvalue shifting consists in moving one or several eigenvalues of a given matrix polynomial to prescribed positions in the complex plane (or to infinity). It can be useful, for instance, as a preliminary modification of the polynomial before applying methods that require a particular eigenvalue distribution.

Chapter 2 introduces the general theory of the invariant pair problem, along with an alternative formulation based on the contour integral. The original contribution of this chapter consists in new formulations for the condition number and the backward error of the invariant pair problem. These formulas are obtained using the definition (1) recalled above.

Chapter 3 is devoted to matrix solvents. A quick review of the general theory, of applications, and of an alternative formulation based on the contour integral are given. Next, we specialize the results of Chapter 2 to solvents, and we compute new characterizations for the condition number and the backward error of the matrix solvent problem. Also, motivated by the results in [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF] and [START_REF] Tisseur | Triangularizing quadratic matrix polynomials[END_REF], we study in Section 3.4 the relation between solvents of general and triangularized matrix polynomials.

Chapter 4 focuses on the computation of eigenvalues and invariant pairs through moments and Hankel pencils. Our starting point here is the Sakurai-Sugiura method ( [START_REF] Sakurai | A projection method for generalized eigenvalue problems using numerical integration[END_REF], [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF]) for computing eigenvalues of a matrix polynomial P (λ). The main idea consists in defining a complex function f (z) whose poles are the eigenvalues of P (λ) inside a given contour Γ, and in computing a few moments of f (z) via contour integrals. The moments are then arranged to form a Hankel matrix pencil, whose generalized eigenvalues are the eigenvalues of P (λ) that belong to the interior of Γ. We generalize and adapt this approach to the computation of invariant pairs with eigenvalues belonging to a prescribed region of the complex plane. We also discuss the effectiveness of the scalar and block versions of the method in presence of multiple eigenvalues. In particular, we show that the scalar method cannot capture some eigenvalue multiplicity structures, for which the block version is needed. Our main results here consist in Theorem 14, Corollary 2 and Theorem 16. Moreover, Section 4.4 addresses some questions on the choice of the contour Γ; its content is not new, but it may offer a useful complement to the topics of this thesis, particularly, regarding the estimation of the number of eigenvalues in a given contour. Finally, Section 4.5.1 presents a theoretical and experimental error analysis for the trapezoid rule applied to our quadrature problems.

The techniques presented in Chapter 4 -as well as other direct approaches to the computation of invariant pairs -can be used either alone or in combination with iterative refinement methods. Motivated by the work in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], in Chapter 5 we propose and compare some variants of Newton's method applied to the numerical refinement of invariant pairs. In particular, we experiment with line search strategies and with Šamanskii's acceleration technique [START_REF] Šamanskii | On a modification of the Newton method[END_REF].

Finally, Chapter 6 presents some conclusions and ideas for future work. The Maple and MATLAB implementations of the symbolic and numeric methods presented in this thesis are available online at the URL http://www.unilim.fr/pages perso/esteban.segura/software.html Chapter 1 :

Matrix Polynomials and the

Eigenvalue Problem

Matrix Polynomials and the Eigenvalue Problem

Let K be an arbitrary field. An n × n matrix polynomial is an n × n matrix whose entries are polynomials in K[λ]. Any matrix polynomial P (λ) can be written in the following form:

Definition 1. An n × n matrix polynomial P (λ) is defined as:

P (λ) = A 0 + A 1 λ + A 2 λ 2 + • • • + A ℓ λ ℓ = ℓ i=0 A i λ i , (1.1)
where ℓ ∈ N is the degree of the matrix polynomial and A 0 , A 1 , ..., A ℓ ∈ K n×n .

A ℓ is called the leading coefficient matrix and A 0 is called the trailing coefficient matrix of P (λ). When A ℓ = I n , the matrix polynomial is said to be monic.

The rank of a matrix polynomial P (λ) is defined as: rank(P (λ)) = max{rank(P (λ 0 )) : λ 0 ∈ K}.

Definition 2. An n × n matrix polynomial P (λ) is said to be regular if rank(P (λ)) = n, or, equivalently, if its determinant det(P (λ)) does not vanish identically. Otherwise, it is said to be singular.

In this work, we assume that P (λ) is regular.

Definition 3. The reversal of the matrix polynomial P (λ) is:

rev(P (λ)) := λ ℓ P (1/λ) = λ ℓ A 0 + λ ℓ-1 A 1 + λ ℓ-2 A 2 + • • • + A ℓ (1.2)
In this work, K denotes either the field of complex numbers C or the field of real numbers R. • The homogeneous formulation of equation (1.1):

P (α, β) = β ℓ A 0 + αβ ℓ-1 A 1 + • • • + α ℓ A ℓ = ℓ i=0 α i β ℓ-i A i ,
for λ = α/β and (α, β) = (0, 0), allows the simultaneous treatment of finite and infinite eigenvalues (see [START_REF] Ahmad | Backward errors and pseudospectra for structured nonlinear eigenvalue problems[END_REF], [START_REF] Higham | Backward error of polynomial eigenproblems solved by linearization[END_REF], [START_REF] Higham | The conditioning of linearizations of matrix polynomials[END_REF]).

• Infinite eigenvalues can still be covered using the reversal matrix polynomial (1.2).

Infinite eigenvalues of the matrix polynomial P (λ) are zero eigenvalues of the reversal matrix polynomial rev(P (λ)).

A particular case of special interest is the quadratic eigenvalue problem (QEP), where ℓ = 2:

Q(λ)x = (A 0 + A 1 λ + A 2 λ 2 )x = 0, x = 0. (1.3)
Typical applications of the QEP include the vibration analysis of buildings, machines and vehicles (see [START_REF] Gohberg | Matrix polynomials[END_REF], [START_REF] Lancaster | Lambda-matrices and vibrating systems[END_REF], [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]). A considerable amount of work has been done on the theoretical and computational study of the QEP: see, for instance, [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF].

Remark 2. The case when ℓ = 1 for the PEP corresponds to the generalized eigenvalue problem (GEP) for matrix pencils:

Ax = λBx, (1.4)
and if, moreover, we have that B = I, we obtain the standard eigenvalue problem:

Ax = λx.

Systems of Ordinary Differential Equations

Matrix polynomials play an important role in the study of ordinary differential equations (ODEs). A system of ODEs of order ℓ > 1, with constant coefficients, takes the form (see, e.g., [START_REF] Gohberg | Matrix polynomials[END_REF]):

ℓ i=0 A i d dt i u(t) = 0. (1.5)
Suppose we are looking for solutions of the form u(t) = x 0 e λ 0 t , where x 0 and λ 0 are independent of t, then (1.5) leads to the polynomial eigenvalue problem P (λ 0 )x 0 = 0, where

P (λ) = ℓ i=0 A i λ i .

Computation of Invariant Pairs and Matrix Solvents

More generally, the function:

u(t) = t k k! x 0 + • • • + t 1! x k-1 + x k e λ 0 t
is a solution of the differential equation if and only if the set of vectors x 0 , x 1 , . . . , x k , with x 0 = 0, satisfies the relations:

j i=0 1 i! P (i) (λ 0 )x j-i = 0, j = 0, 1, . . . , k.
This set of vectors x 0 , x 1 , . . . , x k is called a Jordan chain of length k + 1 associated with the eigenvalue λ 0 and the eigenvector x 0 .

In particular, when ℓ = 2, consider the linear second order differential equation:

M q(t) + C q(t) + Kq(t) = f (t), (1.6)
where M, C, K are n × n matrices and f (t) is a vector. Such equations arise in many engineering applications, for instance, when studying mechanical and electrical oscillation (see [START_REF] Lancaster | Lambda-matrices and vibrating systems[END_REF], [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]).

The homogeneous equation of (1.6) leads to the QEP:

(λ 2 M + λC + K)x = 0
for the solutions of the form q(t) = e λt x.

Several examples of real life problems related to (1.6) have been studied in the literature. For instance, in [START_REF] Abrams | Two coupled oscillator models: the Millennium Bridge and the chimera state[END_REF], [START_REF] Jammaz | Jacobi-Davidson method for polynomial eigenvalue problems[END_REF], [START_REF] Datta | Finite element model updating and partial eigenvalue assignment in structural dynamics: recent developments on computational methods[END_REF], [START_REF] Eckhardt | Modeling walker synchronization on the Millennium Bridge[END_REF], [START_REF] Macdonald | Lateral excitation of bridges by balancing pedestrians[END_REF], [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF] the wobbling of the Millennium bridge over the river Thames in London was discussed. On its opening in 2000, this footbridge started to wobble. It has generally been thought that the Millennium Bridge wobble was due to pedestrians synchronizing their footsteps with the bridge motion. However, this is not supported by measurements of the phenomenon on other bridges. In [START_REF] Macdonald | Lateral excitation of bridges by balancing pedestrians[END_REF], a simple model of human balance strategy for normal walking on a stationary surface was considered. This model showed that pedestrian can effectively act as a negative (or positive) damper to the bridge motion and hence inadvertently feed energy into bridge oscillations.

Two days after the opening, the bridge was closed for almost two years while modifications were made to eliminate the wobble entirely. It reopened in 2002. Another interesting real life example is the project of the company SFE GmbH in Berlin, which investigated the noise in rail traffic that is caused by high speed trains (see [START_REF] Hilliges | On the solution of palindromic eigenvalue problems[END_REF], [START_REF] Ipsen | Accurate eigenvalues for fast trains[END_REF], [START_REF] Mackey | Palindromic polynomial eigenvalue problems: Good vibrations from good linearizations[END_REF]).

Smith Normal Form

Definition 5. A matrix polynomial E(λ) is unimodular if its determinant is a nonzero constant, independent of λ. Definition 6. Two matrix polynomials A(λ) and B(λ) of the same size are called equivalent (A(λ) ∼ B(λ)) if there exist two unimodular matrix polynomials E(λ) and F (λ), such that:

A(λ) = E(λ)B(λ)F (λ).
A crucial and well-known property of matrix polynomials is the existence of the Smith form (see, e.g., [START_REF] Gohberg | Matrix polynomials[END_REF]). Theorem 1. [START_REF] Gohberg | Matrix polynomials[END_REF] Every n×n matrix polynomial P (λ) is equivalent to a matrix polynomial D(λ) of the form

D(λ) = diag(d 1 (λ), . . . , d n (λ)), (1.7)
where the d i 's are monic polynomials such that d i (λ) is divisible by d i-1 (λ).

Definition 7. The matrix polynomial D(λ), given by (1.7), is called the Smith normal form of P (λ), and the monic scalar polynomials d i (λ) are called the invariant polynomials of P (λ).

Solving PEPs via Linearization

Linearization is a well established method for computing eigenvalues and eigenvectors of polynomial eigenvalue problems of moderate size. Given a polynomial eigenvalue problem (1.1) of degree ℓ ≥ 2, the linearization method approach consists in converting P (λ) to a linear ℓn × ℓn pencil L(λ) = A + λB, having the same spectrum as P (λ). This linear eigenvalue problem can then be solved by standard methods, e.g., the QZ algorithm (see A.5).

Software libraries as LAPACK provide routines, among others, for solving systems of linear equations (see [START_REF] Anderson | [END_REF]). On the other hand, the library ARPACK uses the implicitly restarted Arnoldi method or, in the case of symmetric matrices, the corresponding variant of the Lanczos algorithm to compute a few eigenvalues and corresponding eigenvectors of large sparse or structured matrices (see [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]).

We introduce the notion of linearization following [START_REF] Gohberg | Matrix polynomials[END_REF] and [START_REF] Mackey | Vector spaces of linearizations for matrix polynomials[END_REF].

Definition 8. Let P (λ) be an n × n matrix polynomial of degree ℓ ≥ 1. A pencil L(λ) = A + λB with A, B ∈ C ℓn×ℓn is called linearization of P (λ) if there exist unimodular matrix polynomials E(λ), F (λ) ∈ C ℓn×ℓn , such that:

F (λ)L(λ)E(λ) = P (λ) 0 0 I (ℓ-1)n . (1.8)
Remark 3. The linearization (1.8) is not unique (see [START_REF] Mackey | Vector spaces of linearizations for matrix polynomials[END_REF], [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]).

Most of the linearizations used in practice are of the first companion form (see [START_REF] Gohberg | Matrix polynomials[END_REF], [START_REF] Mackey | Palindromic polynomial eigenvalue problems: Good vibrations from good linearizations[END_REF]):

C 1 (λ) =       A ℓ-1 A ℓ-2 • • • A 0 -I n 0 • • • 0 . . . . . . . . . . . . 0 • • • -I n 0       + λ        A ℓ 0 • • • 0 0 I n . . . . . . . . . . . . . . . 0 0 • • • 0 I n       
or the second companion form:

C 2 (λ) =        A ℓ-1 -I n • • • 0 A ℓ-2 0 • • • . . . . . . . . . . . . -I n A 0 0 • • • 0        + λ        A ℓ 0 • • • 0 0 I n . . . . . . . . . . . . . . . 0 0 • • • 0 I n        .
It is well known that using this approach to solve the PEP via linearization can have some drawbacks. For instance:

• The linearization approach transforms the original n×n matrix polynomial of degree ℓ into a larger ℓn × ℓn linear eigenvalue problem.

• The conditioning of the linearized problem can depend on the type of linearization used and may be significantly worse than the conditioning of the original problem (see [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF], [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]).

• If there is a special structure in the matrix coefficients A i , such as symmetry, sparsity pattern, palindromicity, the linearization may modify it. In that case, special linearizations can be chosen to exploit the structure (see [START_REF] Adhikari | On backward errors of structured polynomial eigenproblems solved by structure preserving linearizations[END_REF], [START_REF] Adhikari | Structured eigenvalue condition numbers and linearizations for matrix polynomials[END_REF], [START_REF] Bennerý | Solving large-scale quadratic eigenvalue problems with Hamiltonian eigenstructure using a structure-preserving Krylov subspace method[END_REF], [START_REF] Byers | Trimmed linearizations for structured matrix polynomials[END_REF], [START_REF] De Terán | Linearizations of matrix polynomials: Sharp lower bounds for the dimension and structures[END_REF], [START_REF] De Terán | Palindromic companion forms for matrix polynomials of odd degree[END_REF], [START_REF] Gemignani | The Ehrlich-Aberth method for palindromic matrix polynomials represented in the Dickson basis[END_REF], [START_REF] Higham | Symmetric linearizations for matrix polynomials[END_REF], [START_REF] Higham | The conditioning of linearizations of matrix polynomials[END_REF], [START_REF] Higham | Definite matrix polynomials and their linearization by definite pencils[END_REF], [START_REF] Mackey | Structured linearizations for matrix polynomials[END_REF], [START_REF] Mackey | Structured polynomial eigenvalue problems: Good vibrations from good linearizations[END_REF], [START_REF] Mackey | Vector spaces of linearizations for matrix polynomials[END_REF], [START_REF] Mackey | Smith forms of palindromic matrix polynomials[END_REF], [START_REF] Mackey | Skew-symmetric matrix polynomials and their smith forms[END_REF], [START_REF] Mackey | Möbius transformations of matrix polynomials[END_REF], [START_REF] Mackey | Linearizations of matrix polynomials in Bernstein basis[END_REF], [START_REF] Mehrmann | Structure-preserving methods for computing eigenpairs of large sparse skew-Hamiltonian/Hamiltonian pencils[END_REF], [START_REF] Mehrmann | Polynomial eigenvalue problems with Hamiltonian structure[END_REF]).

Direct Methods to Solve PEPs

In the previous section, we have listed some drawbacks of the linearization method for solving PEPs. We pointed out that this approach not only increases the dimension of the problem, but may also be very ill-conditioned. We recall now some approaches which offer the possibility to handle the PEPs directly without linearization.

• Ehrlich-Aberth iteration: In [START_REF] Bini | Solving polynomial eigenvalue problems by means of the Ehrlich-Aberth method[END_REF], the authors present an effective approach based on the Ehrlich-Aberth iteration to approximate eigenvalues of matrix polynomials. The main points addressed in this work are the choice of the starting approximations, the computation of the Newton correction, the halting criterion and the case of (multiple) eigenvalues at zero and at infinity.

• Jacobi-Davidson method: This method computes selected eigenvalues and associated eigenvectors of a matrix polynomial. It iteratively constructs approximations of certain eigenvectors by solving a projected problem. The method finds the approximate eigenvector as "best" approximation in some search subspace (see [START_REF] Sleijpen | Jacobidavidson type methods for generalized eigenproblems and polynomial eigenproblems[END_REF], [START_REF] Sleijpen | Quadratic eigenproblems are no problem[END_REF]). This approach has been used for the efficient solution of quadratic eigenproblems associated with acoustic problems with damping (see [START_REF] Van Gijzen | The parallel computation of the smallest eigenpair of an acoustic problem with damping[END_REF]).

• A second-order Arnoldi method for the solution of the quadratic eigenvalue problem (SOAR): This method for solving large-scale QEPs generates an orthonormal basis and then applies the standard Rayleigh-Ritz orthogonal projection technique to approximate eigenvalues and eigenvectors (see [START_REF] Bai | SOAR: A second-order Arnoldi method for the solution of the quadratic eigenvalue problem[END_REF]).

• Arnoldi and Lanczos-type methods: These processes are developed to construct projections of the QEP. The convergence of these methods is usually slower than a Krylov subspace method applied to the equivalent linear eigenvalue problem (see [START_REF] Hoffnung | Krylov type subspace methods for matrix polynomials[END_REF]).

• A subspace approximation method: In [START_REF] Holz | A subspace approximation method for the quadratic eigenvalue problem[END_REF], the authors use perturbation subspaces for block eigenvector matrices to reduce a modified problem to a sequence of problems of smaller dimension. They show that this method converges at least as fast as the corresponding Taylor series, and that Rayleigh quotient iteration can be used for acceleration.

• Contour integral based methods: In [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF], [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF] and [START_REF] Sakurai | Efficient parameter estimation and implementation of a contour integral-based eigensolver[END_REF], the authors use contour integral formulations to find all the eigenvalues of PEPs, which are inside a closed contour in the complex plane.

Scaling of Generalized and Polynomial Eigenvalue Problems

The performance of an algorithm may depend crucially on how the problem is formulated.

Balancing is a preprocessing technique that aims to avoid large differences in magnitude among matrix entries, which may cause a poor numerical performance. A matrix with a norm that is several orders of magnitude larger than the modulus of its eigenvalues typically has eigenvalues that are sensitive to perturbations in the entries.

The process of balancing produces a matrix, which is diagonally similar to the given matrix and reduces the matrix norm. As a consequence, the eigenvalues do not change, but their sensitivity can significantly be reduced. Such a diagonal scaling is, therefore, typically used before running any eigenvalue algorithm.

Balancing Technique for Matrices

In the case of a n × n matrix A, balancing consists in finding a diagonal matrix D such that DAD -1 is a well-scaled matrix (see [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF], [START_REF] Osborne | On pre-conditioning of matrices[END_REF], [START_REF] Parlett | Balancing a matrix for calculation of eigenvalues and eigenvectors[END_REF]). For instance, consider the badly scaled matrix A defined by:

A =    0 2 -20 2 -5 2 20 1 0 2 6 2 -15 0    .
If we balance A using the diagonal matrix D = diag(2 -20 , 1, 2 -15 ), we obtain:

D -1 AD =    0 1 1 1 1 0 2 1 0    .
Let us compute the eigenvalues of A and of D -1 AD using the MATLAB command eig eig(A,'nobalance') eig(D -1 AD) Remark 4. There are also cases in which balancing can lead to a catastrophic increase of the errors in the computed eigenvalues (see [START_REF] Watkins | A case where balancing is harmful[END_REF]).

Balancing Technique for Generalized and PEPs

In the case of the generalized eigenvalue problems Ax = λBx, the authors of [START_REF] Ward | Balancing the generalized eigenvalue problem[END_REF] introduce a balancing technique which aims to find diagonal matrices D 1 and D 2 such that the elements of D 1 AD 2 and D 1 BD 2 are scaled as equal in magnitude as possible (see also [START_REF] Hodel | Computation of system zeros with balancing[END_REF]).

A different approach for the scaling of GEPs is proposed in [START_REF] Lemonnier | Balancing regular matrix pencils[END_REF]. A linearly convergent iteration provides matrices D 1 and D 2 consisting of powers of 2 that approximately satisfy:

D 1 AD 2 e j 2 2 + D 1 BD 2 e j 2 2 = e * i D 1 AD 2 2 2 + e * i D 1 BD 2 2 2 = 1, i, j = 1, . . . , n.
In [START_REF] Betcke | Optimal scaling of generalized and polynomial eigenvalue problems[END_REF], besides studying optimal balancing of GEPs and polynomial eigenvalue problems, it is noted that this iteration can easily be extended to weighted scaling of matrix polynomials P (λ) by:

ℓ k=0 ω 2k D 1 A k D 2 e i 2 2 = 1, ℓ k=0 ω 2k e * j D 1 A k D 2 2 2 = 1, i, j = 1, . . . , n
for some ω > 0 that is chosen to be close in magnitude to the desired eigenvalues.

In the specific case of the quadratic matrix polynomial Q(λ), Fan, Lin and Van Dooren [START_REF] Fan | Normwise scaling of second order polynomial matrices[END_REF] suggest that a good scaling strategy for the QEP is to scale the coefficients A 2 , A 1 and A 0 so that their 2-norms are all close to 1. They consider modifying Q(λ) to:

Q(µ) ≡ βQ(λ) = µ 2 Ã2 + µ Ã1 + Ã0 ,
where:

λ = αµ, Ã2 = α 2 βA 2 , Ã1 = αβA 1 , Ã0 = βA 0 , α = a 0 a 2 , β = 2 (a 0 + a 1 α) and a 2 = A 2 2 , a 1 = A 1 2 , a 0 = A 0 2 .
Note that the eigenvalues of Q(λ) can be recovered from those of Q(µ) by λ = αµ. Moreover, this scaling approach does not affect any sparsity of A 2 , A 1 and A 0 .

In [START_REF] Hammarling | An algorithm for the complete solution of quadratic eigenvalue problems[END_REF], it is presented an eigensolver for the complete solution of QEPs (function quadeig in MATLAB), which uses this scaling technique.

A different approach to balance matrix polynomials can be found in [START_REF] Corless | On a generalized companion matrix pencil for matrix polynomials expressed in the Lagrange basis[END_REF].

Shifting Technique for P (λ)

In [START_REF] Meini | A "shift-and-deflate" technique for quadratic matrix polynomials[END_REF], B. Meini describes a technique to shift two eigenvalues λ 1 and λ 2 of an n × n quadratic matrix polynomial P (λ). This method requires the knowledge of a right and of a left eigenvector associated with λ 1 and λ 2 , respectively. The idea is to shift λ 1 to 0 and λ 2 to ∞ and then deflate those values obtaining as a result a new (n -1) × (n -1) quadratic matrix polynomial, which shares with P (λ) all its eigenvalues, except for λ 1 and λ 2 .

As we pointed out before, this approach requires the knowledge of a right and of a left eigenvector associated with the eigenvalues λ 1 and λ 2 . Here, we present a formulation that allow us to shift to infinity several eigenvalues at the same time, assuming only the knowledge of a right eigenvector associated with those eigenvalues to shift.

Let us first recall some results from [START_REF] Meini | A "shift-and-deflate" technique for quadratic matrix polynomials[END_REF].

Meini's Shifting Formulation

Given a n × n matrix polynomial P (λ) of degree ℓ ≥ 1, we have:

Theorem 2. [Thm. 1, [101]] Let λ 1 ∈ C and v ∈ C n , v = 0, such that P (λ 1 )v = 0.
Then for any η ∈ C and for any vector x such that x * v = 1, the following properties hold:

1. The function

P (λ) = P (λ) I + λ 1 -η λ -λ 1 vx *
is a matrix polynomial of degree ℓ;

2. det P (λ) = λ-η λ-λ 1 det P (λ);

3. P (η)v = 0, i.e., the value η, where we moved the original eigenvalue λ 1 , is an eigenvalue for the shifted matrix polynomial. Note that the eigenvector v is the same for λ 1 and η.

4. If σ ∈ {λ 1 , η} is such that P (σ)w = 0, then P (σ) w = 0, where

w = I - λ 1 -η σ -η vx * w.
5. Moreover, by setting:

P (λ) = ℓ i=0 λ i Ãi , one has Ãℓ = A ℓ , Ãi = A i + (λ 1 -η) ℓ j=i+1 λ j-i-1 A j vx * , i = 0, . . . , ℓ -1.
(1.9)

Equations (1.9) give us an easy way to compute the coefficients of the matrix polynomial whose eigenvalue λ has been shifted.

Proposition 1. [Prop. 2, [101]] Let λ 1 ∈ C, v ∈ C n , v = 0, such that P (λ 1 )v = 0 and let x ∈ C n such that x * v = 1.
Then the following properties hold:

1. The matrix polynomial

P (λ) = P (λ) I + λ 1 λ -λ 1 vx *
is such that det P (λ) = det P (λ) λ λ-λ 1 and P (0)v = 0, i.e., the eigenvalue λ 1 is shifted to zero, while keeping the same right eigenvector v;

2. if λ 1 = 0, the matrix polynomial

P (λ) = P (λ) I + λ λ 1 -λ vx *
is such that det P (λ) = det P (λ) λ 1 λ 1 -λ and P (∞)v = 0, i.e., the eigenvalue λ 1 is shifted to infinity, while keeping the same right eigenvector v.

Generalization of Shifting Technique

The results presented in [START_REF] Meini | A "shift-and-deflate" technique for quadratic matrix polynomials[END_REF] allow us to shift a single eigenvalue of a matrix polynomial. Our goal in this section is to present a formulation to shift several eigenvalues at the same time. We start by shifting the eigenvalue λ k to infinity.

Proposition 2. Let P (λ) be a matrix polynomial, λ k = 0 ∈ C and v = 0 ∈ C n , such that P (λ k )v = 0 and let x ∈ C n be such that x * v = 1. Consider the matrix polynomial P (λ) = ℓ i=0 λ i Âi with coefficients:

Â0 = A 0 , Âi = A i + λ -1 k i-1 j=0 λ -j k A i-j-1 vx * , i = 1, . . . , ℓ. (1.10)
Then the following properties hold:

1. P (∞)v = 0, i.e., the eigenvalue λ k is shifted to infinity, while keeping the same right eigenvector v.

2. If λ r ∈ {λ k , ∞} and w = 0 are such that P (λ r )w = 0, then ŵ = 0 such that P (λ r ) ŵ = 0 can be computed by:

ŵ = I - λ r λ k vx * w. (1.11)
Proof. The shift of λ k to ∞ can be proved by applying Theorem 2 to the reversed matrix polynomial rev(P (λ)), by shifting λ -1 k to 0. The matrix polynomial P (λ) is the reversed polynomial of

P (λ) = rev(P (λ)) I + λ -1 k λ -λ -1 k vx * .
Then, we have:

P (λ) = ℓ i=0 λ i A ℓ-i I + λ -1 k λ -λ -1 k vx * .
Adapting (1.9) to last P (λ), we obtain:

Ãℓ = A ℓ-ℓ = A 0 , Ãi = A ℓ-i + λ -1 k ℓ j=i+1 1 λ j-i-1 k A ℓ-j vx * , i = 0, . . . , ℓ -1.
As P (λ) = rev( P (λ)), we get:

Â0 = A 0 , Âi = A i + λ -1 k i-1 j=0 λ -j k A i-j-1 vx * , i = 1, . . . , ℓ.
Now, to prove the second part of the Theorem note that:

P (λ r ) ŵ = P (λ r ) I + λ r λ k -λ r vx * I - λ r λ k vx * w = P (λ r ) I + λ r λ k vx * + λ r λ k -λ r vx * - λ 2 r λ k (λ k -λ r ) vx * w = P (λ r ) I + -λ r λ k + λ 2 r + λ r λ k -λ 2 r λ k (λ k -λ r ) vx * w = P (λ r )w = 0.
Remark 5. The new matrix polynomial P (λ) has the same eigenvalues as P (λ), except for λ k that is shifted to ∞.

Example 1. Consider the following quadratic matrix polynomial discussed in [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]:

P (λ) = λ 2 A 2 + λA 1 + A 0 = λ 2    1 0 0 0 1 0 0 0 0    + λ    -2 0 1 0 0 0 0 0 0    +    1 0 0 0 -1 0 0 0 1    ,
where the eigenvalues of P (λ) are:

λ 1 = 1, λ 2 = 1, λ 3 = 1, λ 4 = -1, λ 5 = ∞ and λ 6 = -∞.
Let's shift to infinity the eigenvalue λ 1 .

Note that if v = [1; 0; 0] is an eigenvector associated with λ 1 , then the vector

x = v is such that x * v = 1.
To compute the coefficients of the shifted matrix polynomial P (λ), we use (1.10):

Â0 = A 0 , Â1 = A 1 + λ -1 1 0 j=0 λ -j A -j vx * = A 1 + λ -1 1 A 0 vx * , Â2 = A 2 + λ -1 1 1 j=0 λ -j 1 A 1-j vx * = A 2 + λ -1 1 A 1 + λ -1 1 A 0 vx * .
Then we obtain:

P (λ) = λ 2 Â2 + λ Â1 + Â0 = λ 2    0 0 0 0 1 0 0 0 0    + λ    -1 0 1 0 0 0 0 0 0    +    1 0 0 0 -1 0 0 0 1    ,
where the eigenvalues of P (λ) are: λ1 = ∞, λ2 = 1, λ3 = -1, λ4 = 1, λ5 = ∞ and λ6 = -∞. The associated eigenvector to the shifted eigenvalue remains equal, i.e., P (∞)v = 0, where v = [1; 0; 0]. Moreover, suppose we want to compute the eigenvector v4 such that P (-1)v 4 = 0. From formula (1.11), we have:

v4 = I - λ 4 λ 1 vx * v 4 = (I + vx * ) v 4 = [0; -1; 0],
where v 4 = [0; -1; 0] is such that P (-1)v 4 = 0.

In general, to shift m eigenvalues to ∞, we have the following more general result.

Theorem 3. Let P (λ) be an n × n matrix polynomial and its eigenvalues λ i = 0 ∈ C with associated eigenvectors v i = 0 ∈ C n , for i = 1, . . . , m. Consider the n × n matrix polynomial P (λ) = ℓ j=0 λ j m A j , with coefficients:

m A 0 = m-1 A 0 = • • • = 0 A 0 = A 0 , m A i = m-1 A i + λ -1 m i-1 j=0 λ -j m m-1 A i-j-1 m-1 v m m-1 x * m , i = 1, . . . , ℓ, (1.12)
where 0 A i = A i . Then the following properties hold:

• P (∞)v i = 0, i.e., the eigenvalues λ i are shifted to infinity, while keeping the same right eigenvectors v i , for i = 1, . . . , m.

• If P (λ k+1 )v k+1 = 0, then the vector

k v k+1 = 0 ∈ C n such that P (λ k+1 ) k v k+1 = 0, can be computed by: k v k+1 = k-1 j=0 I - λ k+1 λ k-j k-j-1 v k-j k-j-1 x * k-j v k+1 , k = 1, . . . , m -1, (1.13)
where the vectors

k-j-1 x * k-j ∈ C n satisfy k-j-1 x * k-j k-j-1 v k-j = 1 and 0 v 1 = v 1 , 0 x 1 = x 1 .
Proof. To get the coefficients m A i , for i = 0, . . . , ℓ, note that we just must apply m times the Proposition 2, i.e., the matrix polynomial P (λ) is:

P (λ) = P (λ) I + λ λ 1 -λ v 1 x * 1 I + λ λ 2 -λ v 2 x * 2 • • • I + λ λ m -λ v m x * m w.
Note that:

2 P (λ) = P (λ) I + λ λ 1 -λ v 1 x * 1 I + λ λ 2 -λ v 2 x * 2 = 1 P (λ) I + λ λ 2 -λ v 2 x * 2 ,
where the coefficients of 1 P (λ) are:

1 A 0 = A 0 , 1 A i = A i + λ -1 1 i-1 j=0 λ -j 1 A i-j-1 v 1 x * 1 , i = 1, . . . , ℓ.
Then, the coefficients for 2 P (λ) are given by:

2 A 0 = 1 A 0 = A 0 , 2 A i = 1 A i + λ -1 2 i-1 j=0 λ -j 2 1 A i-j-1 1 v 2 1
x * 2 , i = 1, . . . , ℓ.

Continuing this m times we obtain that the coefficients for P (λ) are given by:

m A 0 = m-1 A 0 = • • • = 0 A 0 = A 0 , m A i = m-1 A i + λ -1 m i-1 j=0 λ -j m m-1 A i-j-1 m-1 v m m-1
x * m , i = 1, . . . , ℓ,

for k v k+1 = k-1 j=0 I - λ k+1 λ k-j k-j-1 v k-j k-j-1 x * k-j v k+1 , k = 1, . . . , m -1,
where the vectors

k-j-1 x * k-j ∈ C n satisfy k-j-1 x * k-j k-j-1 v k-j = 1 and 0 v 1 = v 1 , 0 x 1 = x 1 .
Example 2. Consider the quadratic matrix polynomial P (λ) of Example 1 with eigenvalues:

λ 1 = 1, λ 2 = -1, λ 3 = 1, λ 4 = 1, λ 5 = ∞ and λ 6 = -∞.
Let's shift to infinity the eigenvalues λ 1 = 1 and λ 2 = -1 with corresponding associated eigenvectors v 1 = [1; 0; 0] and v 2 = [0; -1; 0]. Consider the vectors:

x 1 = v 1 and x 2 = v 2 such that x * i v i = 1, for i = 1, 2.
To compute the coefficients of the shifted matrix polynomial P (λ), we use (1.12):

2 A 0 = 1 A 0 = 0 A 0 = A 0 , 2 A 1 = 1 A 1 + λ -1 2 1 A 0 1 v 2 1 x * 2 = 0 A 1 + λ -1 1 0 A 0 0 v 1 0 x * 1 + λ -1 2 A 0 1 v 2 1 x * 2 = A 1 + λ -1 1 A 0 v 1 x * 1 + λ -1 2 A 0 1 v 2 1 x * 2 , 2 A 2 = 1 A 2 + λ -1 2 1 j=0 λ -j 2 1 A 1-j 1 v 2 1 x * 2 = 1 A 2 + λ -1 2 1 A 1 + λ -1 2 1 A 0 1 v 2 1 x * 2 = 0 A 2 + λ -1 1 1 j=0 λ -j 1 0 A 1-j 0 v 1 0 x * 1 + λ -1 2 1 A 1 + λ -1 2 1 A 0 1 v 2 1 x * 2 = A 2 + λ -1 1 A 1 + λ -1 1 A 0 v 1 x * 1 + λ -1 2 1 A 1 + λ -1 2 A 0 1 v 2 1 x * 2 = A 2 + λ -1 1 A 1 + λ -1 1 A 0 v 1 x * 1 + λ -1 2 A 1 + λ -1 1 A 0 v 1 x * 1 + λ -1 2 A 0 1 v 2 1 x * 2 .
Now, we need to compute

1 v 2 ∈ C n .
Using (1.13), we obtain:

1 v 2 = I - λ 2 λ 1 v 1 x * 1 v 2 =    0 -1 0    .
Taking

1 x 2 = 1 v 2 , satisfying 1 x * 2 1 v 2 = 1
, we obtain:

P (λ) = λ 2 2 A 2 + λ 2 A 1 + 2 A 0 = λ 2    0 0 0 0 0 0 0 0 0    + λ    -1 0 1 0 1 0 0 0 0    +    1 0 0 0 -1 0 0 0 1    .
Note that P (λ) has the eigenvalues:

λ1 = ∞, λ2 = ∞, λ3 = 1, λ4 = 1, λ5 = ∞ and λ6 = -∞.
The following result is a generalization of Theorem 2. It shows how to shift m eigenvalues λ i to the values η 1 , . . . , η m , with η i = ∞, for i = 1, . . . , m.

Theorem 4. Let P (λ) be an n × n matrix polynomial and its eigenvalues

λ i = 0 ∈ C with associated eigenvectors v i = 0 ∈ C n , for i = 1, . . . , m. Consider the n × n matrix polynomial P (λ) = ℓ j=0 λ j m A j , with coefficients: m A ℓ = m-1 A ℓ = • • • = 0 A ℓ = A ℓ , m A i = m-1 A i + (λ m -η m ) ℓ j=i+1 λ j-i-1 m m-1 A j m-1 v m m-1 x * m , i = 0, . . . , ℓ -1, (1.14)
where 0 A i = A i . Then the following properties hold:

• P (η i )v i = 0, i.e., the eigenvalues λ i are shifted to η i , while keeping the same right eigenvectors v i , for i = 1, . . . , m.

• If P (λ k+1 )v k+1 = 0, then the vector k v k+1 = 0 ∈ C n such that P (λ k+1 ) k v k+1 = 0, can be computed by: k v k+1 = k-1 j=0 I - λ k-j -η k-j λ k+1 -η k-j k-j-1 v k-j k-j-1 x * k-j v k+1 , k = 1, . . . , m -1, (1.15)
where the vectors

k-j-1 x * k-j ∈ C n satisfy k-j-1 x * k-j k-j-1 v k-j = 1 and 0 v 1 = v 1 , 0 x 1 = x 1 .

Computation of Invariant Pairs and Matrix Solvents

Proof. For the proof of this Theorem just note that we can apply Theorem 3 to the reversed matrix polynomial of P (λ) and by shifting λ i to η i .

Example 3. Consider again the quadratic matrix polynomial P (λ) of Example 1 with eigenvalues:

λ 1 = 1, λ 2 = -1, λ 3 = 1, λ 4 = 1, λ 5 = ∞ and λ 6 = -∞.
Let's shift the eigenvalues λ 1 , λ 2 and λ 3 to the values: η 1 = -5, η 2 = 0 and η 3 = 10, respectively. The corresponding associated eigenvectors with the λ i 's are

v 1 = [1; 0; 0], v 2 = [0; -1; 0] and v 3 = [0; 1; 0].
Consider the vectors:

x 1 = v 1 and x 2 = v 2 such that x * i v i = 1, for i = 1, 2, 3.
Then, using (1.14) and (1.15), we find the coefficients of the shifted matrix polynomial P ( λ):

P ( λ) = λ2 Ã2 + λ Ã1 + Ã0 = λ2    1 0 0 0 1 0 0 0 0    + λ    4 0 1 0 -10 0 0 0 0    +    -5 0 0 0 0 0 0 0 1    .
The eigenvalues of P (λ) are: λ1 = -5, λ2 = 0, λ3 = 10, λ4 = 1, λ5 = ∞ and λ6 = -∞.

Chapter 2 :

Invariant Pairs: Theory, Conditioning and Backward Error

Introduction

Invariant pairs, first introduced in [START_REF] Beyn | Continuation of invariant subspaces for parameterized quadratic eigenvalue problems[END_REF] and further developed in [START_REF] Kressner | A block newton method for nonlinear eigenvalue problems[END_REF] and [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], are a generalization of the notion of eigenpair for matrix polynomials. The notion of invariant pair offers a unified theoretical perspective on the problem of computing several eigenvalueeigenvector pairs for a given matrix polynomial and extends the well-known concepts of standard pair [START_REF] Gohberg | Matrix polynomials[END_REF] and null pair [START_REF] Ball | Interpolation of rational matrix functions[END_REF] (or Jordan pair [START_REF] Gohberg | Matrix polynomials[END_REF]). As noted in the next section, the notion of invariant pairs can also be applied to more general nonlinear problems, but in this thesis, we will limit our presentation to matrix polynomials.

In the following section, we present a brief description of the existing theory on this topic found in some of the references cited above.

Previous Work on Invariant Pairs

The notion of invariant pairs for a quadratic matrix polynomial Q(λ) is introduced in [START_REF] Beyn | Continuation of invariant subspaces for parameterized quadratic eigenvalue problems[END_REF]. In this work, W.-J. Beyn and V. Thümmler consider the quadratic matrix polynomials that depend on the parameter s ∈ R:

Q(λ, s) = A 2 (s)λ 2 + A 1 (s)λ + A 0 (s), λ ∈ C
where it is assumed that A 2 (s), A 1 (s) and A 0 (s) are real square matrices depending smoothly on s. The authors investigate under which conditions smooth solution branches (X, Λ) = (X(s), Λ(s)) of the equation:

Q(Λ, s)X = A 2 (s)XΛ 2 + A 1 (s)XΛ + A 0 (s)X = 0 (2.1)
exist and how to compute them in an efficient way.

For the computation of (X(s), Λ(s)), the work presents a continuation method of predictor-corrector type (see [START_REF] Deuflhard | Newton methods for nonlinear problems: affine invariance and adaptive algorithms[END_REF], [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF]) applied directly on (2.1), thus avoiding the linearization of the problem. Moreover, in the correction step of the method, a Newton-like process to generate the sequence (X i , Λ i , s i ) is studied. Finally, the method is demonstrated on several numerical examples such as: a homotopy between random matrices, a model of fluid conveying pipe problem and a traveling wave of a damped wave equation.

In [START_REF] Kressner | A block newton method for nonlinear eigenvalue problems[END_REF], D. Kressner, inspired by the work in [START_REF] Beyn | Continuation of invariant subspaces for parameterized quadratic eigenvalue problems[END_REF], studies a generalization of the notion of invariant pairs. To this end, consider the nonlinear eigenvalue problem:

(f 1 (λ)A 1 + f 2 (λ)A 2 + • • • + f m (λ)A m )x = 0 (2.2) for holomorphic functions f 1 , . . . , f m : Ω → C (where Ω ⊆ C is an open set) and constant matrices A 1 , . . . , A m ∈ C n×n .
In this case, an invariant pair (X, S) ∈ C n×k × C k×k of (2.2) is defined by:

A 1 Xf 1 (S) + A 2 Xf 2 (S) + • • • + A m Xf m (S) = 0,
where the eigenvalues of S are contained in Ω. Moreover, the work presents a block Newton method to compute such invariant pairs. An inverse iteration approach is used for finding the initial pair (X 0 , S 0 ). It is explained that this method inherits the disadvantages of similar methods for solving linear eigenvalue problems: for instance, its global convergence may be erratic and a single slowly converging eigenvalue contained in S will hinder the convergence of the entire pair.

Finally, T. Betcke and D. Kressner focus on the invariant pair problem for matrix polynomials in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF]. They study the behavior of invariant pairs under perturbations of the matrix polynomial, and a first-order perturbation expansion is given. From a computational point of view, different ways to extract invariant pairs from a linearization of the matrix polynomial are analyzed. Moreover, the authors describe a refinement procedure based on Newton's method and applied directly on the polynomial formulation. This computational approach is tested in some numerical experiments.

Later work on invariant pairs can be found in [START_REF] Effenberger | Robust solution methods for nonlinear eigenvalue problems[END_REF], [START_REF] Effenberger | Robust successive computation of eigenpairs for nonlinear eigenvalue problems[END_REF], [START_REF] Effenberger | Chebyshev interpolation for nonlinear eigenvalue problems[END_REF], [START_REF] Jarlebring | Computing a partial Schur factorization of nonlinear eigenvalue problems using the infinite Arnoldi method[END_REF], [START_REF] Kressner | Memory-efficient Arnoldi algorithms for linearizations of matrix polynomials in Chebyshev basis[END_REF], [START_REF] Mackey | Möbius transformations of matrix polynomials[END_REF], [START_REF] Szyld | Several properties of invariant pairs of nonlinear algebraic eigenvalue problems[END_REF].

Definition and Theory

Definition 9. Let P (λ) be an n × n matrix polynomial. A pair (X, S) ∈ C n×k × C k×k , X = 0, is called an invariant pair if it satisfies the relation:

P (X, S) := A ℓ XS ℓ + • • • + A 2 XS 2 + A 1 XS + A 0 X = 0, (2.3)
where A i ∈ C n×n , i = 0, . . . , ℓ, and k is an integer between 1 and nℓ.

Example 4. Consider again the matrix polynomial of Example 4:

P (λ) = λ 2    1 0 0 0 1 0 0 0 0    + λ    -2 0 1 0 0 0 0 0 0    +    1 0 0 0 -1 0 0 0 1    .
Then (X, S), with:

X =    0 1 0 1 0 1 0 0 0    , S =    1 0 0 0 1 1 0 0 1    ,
is an invariant pair for P (λ).

The eigenvalues of P (λ) are: λ 0 = 1 with multiplicity 3, λ 1 = -1, λ 2 = ∞ and λ 3 = -∞. Note that the eigenvalues of S are λ S = 1 with multiplicity 3.

Remark 6.

• Infinite eigenvalues can still be covered by defining invariant pairs for the reversal polynomial (1.2). If a polynomial has zero and infinite eigenvalues, they have to be handled by separate invariant pairs, one for the original and one for the reverse polynomial.

• The definition of an invariant pair is independent of the choice of basis. Indeed, let T ∈ C k×k be an invertible matrix and consider X = XT , S = T -1 ST . Then, by multiplying (2.3) by T from the right, we obtain:

A 0 XT + A 1 XST + A 2 XS 2 T + • • • + A ℓ XS ℓ T = A 0 X + A 1 XT -1 ST + A 2 XT -1 S 2 T + • • • + A ℓ XT -1 S ℓ T = A 0 X + A 1 X S + A 2 X S2 + • • • + A ℓ X Sℓ = 0. (2.4)
Hence, ( X, S) is also an invariant pair.

Note that if S is diagonalizable, then T can be chosen such that: S = T -1 ST = diag(λ 1 , . . . , λ k ). More generally, S can always be chosen in Schur (upper triangular) form.

• The relation (2.4) implies that the columns x1 , x2 , . . . , xk of the transformed basis X are eigenvectors of P (λ): P (λ i )x i = 0, xi = 0. In particular, S can be transformed into Jordan form: it is then easy to see that the eigenvalues of S form a subset of the eigenvalues of P (λ).

The following definitions proposed in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF] and [START_REF] Gohberg | Matrix polynomials[END_REF] will be helpful for our work, for instance, to allow for rank deficiency in X.

Definition 10. A pair (X, S) ∈ C n×k × C k×k is called minimal if there is m ∈ N * such that: V m (X, S) :=       XS m-1 . . . XS X      
has full rank. The smallest such m is called minimality index of (X, S).

Definition 11. An invariant pair (X, S) for a regular matrix polynomial P (λ) of degree ℓ is called simple if (X, S) is minimal and the algebraic multiplicities of the eigenvalues of S are identical to the algebraic multiplicities of the corresponding eigenvalues of P (λ).

It is well known that eigenvectors associated with a multiple eigenvalue are unstable under perturbations, meaning that an arbitrarily small change in the matrix may cause some of the eigenvectors disappear. In contrast, the notion of invariant pairs offers a theoretical perspective and a numerically more stable approach to the task of computing several eigenpairs of a matrix polynomial.

In particular, simple invariant pairs play an important role when using a linearization approach as in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], and ensure local quadratic convergence of Newton's method, as shown in [START_REF] Kressner | A block newton method for nonlinear eigenvalue problems[END_REF]; see also [START_REF] Szyld | Several properties of invariant pairs of nonlinear algebraic eigenvalue problems[END_REF].

Particular Case: Jordan Pairs

Invariant pairs are closely related to the theory of standard pairs presented in [START_REF] Gohberg | Matrix polynomials[END_REF], and, in particular, to Jordan pairs. If (X, S) is a simple invariant pair and S is in Jordan form, then (X, S) is a Jordan pair.

As an example, consider the quadratic matrix polynomial of Example 4:

P (λ) = λ 2    1 0 0 0 1 0 0 0 0    + λ    -2 0 1 0 0 0 0 0 0    +    1 0 0 0 -1 0 0 0 1    ,
and its eigenvalues: λ 0 = 1, with algebraic multiplicity 3, and λ 1 = -1, with algebraic multiplicity 1. A corresponding Jordan pair (X, J) is given by:

X =    0 0 1 0 1 1 0 1 0 0 0 0    , J = diag -1, 1, 1 1 0 1 .

Formulation of the Invariant Pair Problem Using the Contour Integral

Polynomial eigenpairs and invariant pairs can also be defined in terms of a contour integral. Indeed, an equivalent representation for (2.3) is the following.

Proposition 3. [42]

A pair (X, S) ∈ C n×k × C k×k is an invariant pair if and only if satisfies the relation:

P (X, S) := 1 2πı Γ P (λ)X(λI -S) -1 dλ = 0, (2.5)
where Γ ⊆ C is a contour with the spectrum of S in its interior.

This formulation allows us to choose the contour Γ to compute S with eigenvalues lying in a particular region of the complex plane. See [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF], [START_REF] Bart | A state space approach to canonical factorization with applications[END_REF], [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF], [START_REF] Beyn | Continuation of eigenvalues and invariant pairs for parameterized nonlinear eigenvalue problems[END_REF], [START_REF] Effenberger | Robust successive computation of eigenpairs for nonlinear eigenvalue problems[END_REF], [START_REF] Gohberg | General theory of regular matrix polynomials and band Toeplitz operators[END_REF] for applications of the contour integral formulation.

Linearized Matrix Equation and Fréchet Derivative

In this chapter, we analyze the condition number for P (X, S). To this end, we need to study the matrix equation:

(P + ∆P )(X + ∆X, S + ∆S) = 0, (2.6)
where ∆P (λ) is the perturbation of the matrix polynomial P (λ) defined as:

∆P (λ) = ∆A 0 + ∆A 1 λ + ∆A 2 λ 2 + • • • + ∆A ℓ λ ℓ .
Consider ∆X < ǫ, ∆S < ǫ, ∆A i < ǫ for some sufficiently small ǫ > 0. Omitting terms of order O(ǫ 2 ) as ǫ → 0, we have the linearized system:

L P (∆X, ∆S) = -∆P (X, S) (2.7)
with:

L P : (∆X, ∆S) → P (∆X, S) + ℓ j=1 A j XDS j (∆S),
where DS j denotes the Fréchet derivative of the map S → S j :

DS j : ∆S → j-1 i=0 S i ∆SS j-i-1 (2.8)
For example, for ℓ = 2, we have:

L P (∆X, ∆S) = A 0 ∆X + A 1 ∆XS + A 2 ∆XS 2 + A 1 X∆S + A 2 X(∆SS + S∆S).

Condition Number and Backward Error for the Invariant Pair Problem

In the following sections, we present formulations of the backward error and condition number for an invariant pair (X, S) of the matrix equation (2.3). We follow the ideas presented in the articles [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF] and [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF], which give expressions for backward errors and condition numbers for the polynomial eigenvalue problem and for a solvent of the quadratic matrix equation.

Condition Number for P (X, S)

A normwise condition number of the invariant pair (X, S) can be defined as:

κ(X, S) = lim sup ǫ→0            1 ǫ ∆X ∆S F X S F : (P + ∆P )(X + ∆X, S + ∆S) = 0, ∆A i F ≤ ǫα i , i = 0, . . . , ℓ (2.9)
The α i are nonnegative weights that provide flexibility in how the perturbations are measured. A common choice is α i = A i F ; however, if some coefficients are to be left unperturbed, ∆A i can be forced to zero by setting α i = 0.

Theorem 5. The normwise condition number of the simple invariant pair (X, S) is given by:

κ(X, S) = B X B S + B A 2 X S F , (2.10)
where

B X = ℓ j=0 (S j ) T ⊗ A j , B S = ℓ j=1 j-1 i=0 (S j-i-1 ) T ⊗ A j XS i , B A = α ℓ (XS ℓ ) T ⊗ I n • • • α 0 X T ⊗ I n .
Proof. By expanding the first constraint in (2.9) and keeping only the first order terms, we get:

ℓ j=0 ∆A j XS j + ℓ j=0 A j ∆XS j + ℓ j=1 A j XDS j (∆S) = O(ǫ 2 ), (2.11)
where DS j denotes the Fréchet derivative (2.8). Using on equation (2.11) the vec operator (see A.2), we obtain:

• vec (∆P (X, S)) = vec   ℓ j=0 ∆A j XS j   = ℓ j=0 vec ∆A j XS j = ℓ j=0 (XS j ) T ⊗ I n vec (∆A j ) = = α ℓ (XS ℓ ) T ⊗ I n • • • α 0 X T ⊗ I n     vec (∆A ℓ )/α ℓ . . . vec (∆A 0 )/α 0     =: B A vec (∆A), • vec (P (∆X, S)) = vec   ℓ j=0 A j ∆XS j   = ℓ j=0 vec A j ∆XS j = ℓ j=0 (S j ) T ⊗ A j vec (∆X) = =: B X vec (∆X), • vec   ℓ j=1 A j XDS j (∆S)   = vec   ℓ j=1 A j X j-1 i=0 S i ∆SS j-i-1   = ℓ j=1 j-1 i=0 vec A j XS i ∆SS j-i-1 = = ℓ j=1 j-1 i=0 (S j-i-1 ) T ⊗ A j XS i vec (∆S) =: B S vec (∆S).
Then, we have:

B X B S y = -B A x + O(ǫ 2 ),
where

y = vec (∆X) vec (∆S) , and x =     vec (∆A ℓ )/α ℓ . . . vec (∆A 0 )/α 0    
and therefore

y 2 = vec (∆X) vec (∆S) 2 = ∆X ∆S F .
So we have that the definition (2.9) is equivalent to the following

lim sup ǫ→0            1 ǫ y 2 X S F : B X B S y = -B A x + O(ǫ 2 ), x 2 ≤ ǫ            = B X B S + B A 2 X S F ,
where the matrix B X B S has full rank if the invariant pair (X, S) is simple (see [Thm. 7, [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF]]).

Case k=1

In order to better illustrate Theorem 5, let us consider the particular case k = 1. When k = 1, invariant pairs (X, S) coincide with eigenpairs (x, λ). In this case, the matrices B X , B S and B A in (2.10) are:

B X = ℓ j=0 (λ j ) T ⊗ A j = ℓ j=0 λ j A j = P (λ), B S = ℓ j=1 j-1 i=0 (λ j-i-1 ) T ⊗ A j xλ i = ℓ j=1 j-1 i=0 λ j-1 A j x = P ′ (λ)x, B A = α ℓ λ ℓ x T ⊗ I n α ℓ-1 λ ℓ-1 x T ⊗ I n • • • α 0 x T ⊗ I n
Note that:

B A x = α ℓ λ ℓ x T ⊗ I n • • • α 0 x T ⊗ I n     vec (∆A ℓ )/α ℓ . . . vec (∆A 0 )/α 0     = = vec (λ ℓ ∆A ℓ x + • • • + ∆A 0 x) = vec (∆P (λ)x) = ∆P (λ)x.
Therefore, we obtain:

B X B S y = -B A x + O(ǫ 2 ) ⇔ P (λ) P ′ (λ)x ∆x ∆λ = -∆P (λ)x + O(ǫ 2 ) ⇔ P (λ)∆x + P ′ (λ)x∆λ + ∆P (λ)x = O(ǫ 2 ).
The last equation is consistent with the first part of the computation of the condition number for a nonzero simple eigenvalue λ of P (λ) presented in [Thm. 5, [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF]]. The second part differs, because here we are estimating ∆x ∆λ F , whereas classical condition numbers for eigenvalue problems typically take into account angles between left and right eigenvectors. Of course, it would also be interesting to formalize a similar approach for invariant pairs, based on angles between suitable matrix manifolds (such as partially developed in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF]).

Remark 7. The nonnegative parameters α i are scaling parameters, which can be chosen accordingly to the problem. In practice, we choose them as: α i = A i F , for i = 0, . . . , ℓ.

Numerical Examples

Example 5. Consider Example 4 and the corresponding invariant pair given by:

X =    0 1 0 1 0 1 0 0 0    , S =    1 0 0 0 1 1 0 0 1    .
Using the values α 2 = A 2 F , α 1 = A 1 F and α 0 = A 0 F , the condition number for this problem is: κ(X, S) = 3.8057.

Example 6. In [START_REF] Kressner | A block newton method for nonlinear eigenvalue problems[END_REF], the following matrix polynomial was discussed:

T (λ) = λ 2 1 0 0 1 + λ -1 -6 2 -9 + 0 12 -2 14 . (2.
12)

It has eigenvalues λ 1 = 1, λ 2 = 2, λ 3 = 3 and λ 4 = 4.
An invariant pair is given by:

X = 1 1 1 1 , S = 3 0 0 4 .
Using the values α 2 = A 2 F , α 1 = A 1 F and α 0 = A 0 F , the condition number for this problem is: κ(X, S) = 49.1339.

Example 7. Consider the quadratic matrix polynomial Q(λ) = λ 2 A 2 + λA 1 + A 0 with coefficients (see [START_REF] Kim | Numerical methods for solving a quadratic matrix equation[END_REF]):

A 2 = 1 -1 1 -1 , A 1 = A 2 -τ 1 + τ 0 0 1 , A 0 = -(A 2 + A 1 ). (2.13)
Note that if τ = 0 and S = I 2 , we have:

A 2 XS 2 + A 1 XS + A 0 X = A 2 X + A 2 X + -2A 2 X = 0.
Then, any choice of X ∈ C 2×2 will form an invariant pair (X, I 2 ). Note also that the matrix P = B X B S in the formula for the condition number (2.10) is:

P = ((S 2 ) T ⊗ A 2 + S T ⊗ A 1 + I k ⊗ A 0 ) (I k ⊗ A 2 XS + S T ⊗ A 2 X + I k ⊗ A 1 X) = = (I 2 ⊗ A 2 + I 2 ⊗ A 2 + I 2 ⊗ -(A 2 + A 2 )) (I 2 ⊗ A 2 X + I 2 ⊗ A 2 X + I 2 ⊗ A 2 X) = = 0 3I 2 ⊗ A 2 X ∈ C 4×8 .
Here, the matrix P has rank 2, i.e., it is rank deficient and thus if τ = 0, then we have: We show that κ(X, S) tends to ∞ as τ → 0.

κ(X, S) = ∞. Now, we choose X = I 2 , α 2 = A 2 F , α 1 = A 1 F , α 0 = A 0 F and

Backward Error for P (X, S)

Let α i , for i = 0, . . . , ℓ, be nonnegative weights as in Section 2.6.1. The backward error of a computed solution ( X, S) ∈ C n×k × C k×k to (2.3) can be defined as:

η( X, S) = min{ǫ : (P + ∆P )( X, S) = 0, ∆A i F ≤ ǫα i , i = 0, . . . , ℓ} (2.14)
By expanding the first constraint in (2.14) we get:

-P ( X, S) = ∆A ℓ X Sℓ + • • • + ∆A 0 X. (2.15)
Then, we have

-P ( X, S) = α -1 ℓ ∆A ℓ . . . α -1 1 ∆A 1 α -1 0 ∆A 0        α ℓ X Sℓ . . . α 1 X S α 0 X       
Taking the Frobenius norm, we obtain the lower bound for the backward error:

η( X, S) ≥ P ( X, S) F (α 2 ℓ X Sℓ 2 F + • • • + α 2 1 X S 2 F + α 2 0 X 2 F ) 1/2 .
Consider again equation (2.15). Using the vec operator (see A.2), we obtain:

-vec (P ( X, S)) = (( X Sℓ )

T ⊗ I n ) vec (∆A ℓ ) + • • • + ( XT ⊗ I n ) vec (∆A 0 ) = α ℓ ( X Sℓ ) T ⊗ I n . . . α 0 XT ⊗ I n     vec (∆A ℓ )/α ℓ . . . vec (∆A 0 )/α 0     ,
which can be written as:

Hz = r, H ∈ C nk×(ℓ+1)n 2 (2.16)
Here we assume that H is full rank , to guarantee that (2.16) has a solution (backward error is finite). Then the backward error is the minimum 2-norm solution to:

η( X, S) = H + r 2 , (2.17)
where H + denotes the Moore-Penrose pseudoinverse of H + . Eq. (2.17) yields an upper bound for η( X, S):

η( X, S) ≤ H + 2 r 2 = r 2 σ min (H) ,
where σ min denotes the smallest singular value, which is nonzero by assumption. Note that:

σ min (H) 2 = λ min (HH * ) = λ min (α 2 ℓ ( X Sℓ ) T X Sℓ ⊗ I n + • • • + α 2 0 XT X ⊗ I n ) ≥ ≥ α 2 ℓ σ min ( X Sℓ ) 2 + • • • + α 2 1 σ min ( X S) 2 + α 2 0 σ min ( X) 2 .
Thus we obtain the upper bound for η( X, S):

η( X, S) ≤ P ( X, S) F (α 2 ℓ σ min ( X Sℓ ) 2 + • • • + α 2 1 σ min ( X S) 2 + α 2 0 σ min ( X) 2 ) 1/2
.

Case k=1:

In the particular case k = 1, the approximate invariant pair ( X, S) coincides with an approximate eigenpair (x, λ). In this case, the definition (2.14) becomes:

η(x, λ) = min{ǫ : (P + ∆P )(x, λ) = 0, ∆A i F ≤ ǫα i , i = 0, . . . , ℓ},
which is the definition of the normwise backward error of an approximate eigenpair (x, λ) for P (λ)x = 0, presented in [(2.2), [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF]].

Example 8. Let us consider the power plant problem presented in [START_REF] Betcke | NLEVP: A collection of nonlinear eigenvalue problems[END_REF] and in [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. This is a real symmetric QEP, with P (λ) of size 8 × 8, which describes the dynamic behaviour of a nuclear power plant simplified into an eight-degrees-of-freedom system. The problem is ill-conditioned due to the bad scaling of the matrix coefficients. The maximum condition number for the eigenvalues of P (λ), computed by the MAT-LAB function polyeig, is:

κ max = max λ∈Λ condeig λ = 1.0086e+008.
Using the method that will be presented in Section 4.3.1 and Section 4.5, we compute an invariant pair (X, S) associated with the 11 eigenvalues with largest condition number inside the contour Γ = γ + ρe iθ (γ = 80 + 10i, ρ = 170). The condition number and backward error for (X, S) are κ(X, S) = 565.6746 and η(X, S) = 4.4548e -017.

Observe that κ(X, S) is significantly smaller than κ max .

Matrix Solvents: Definition and Theory

In this section we study the matrix solvent problem as a particular case of the invariant pair problem, and we apply to solvents some results we have obtained for invariant pairs. Definition 12. Let P (λ) be an n × n matrix polynomial. A matrix S ∈ C n×n is called a (right) solvent for P (S) if satisfies the relation:

P (S) := A ℓ S ℓ + • • • + A 2 S 2 + A 1 S + A 0 = 0. (3.1)
A special case is, for ℓ = 2, the quadratic matrix equation:

Q(S) := A 2 S 2 + A 1 S + A 0 = 0,
which has received considerable attention in the literature. For instance, in [START_REF] Higham | Numerical analysis of a quadratic matrix equation[END_REF] and [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF] the authors find formulations for the condition number and the backward error. They also propose functional iteration approaches based on Bernoulli's method and Newton's method with line search to compute the solution numerically. The relation between eigenvalues of P (λ) and solvents is highlighted in [START_REF] Lancaster | Lambda-matrices and vibrating systems[END_REF]: a corollary of the generalized Bézout theorem states that if S is a solvent of P (S), then:

P (λ) = L(λ)(λI -S),
where L(λ) is a matrix polynomial of degree ℓ -1. Then any eigenpair of the solvent S is an eigenpair of P (λ).

An equivalent representation for (3.1) that uses the contour integral is as follows.

Proposition 4. S ∈ C n×n is a solvent if and only if

P (S) := 1 2πı Γ P (λ)(λI -S) -1 dλ = 0, (3.2)
for any closed contour Γ ⊆ C with the spectrum of S in its interior.

As for invariant pairs, this formulation allows us to choose the contour Γ to compute S with specific eigenvalues lying in a particular region of the complex plane.

Existence of Solvents

Let us recall some results that will be needed later. The next result is a generalization of a theorem presented in [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF] which gives information about the number of solvents of P (S). Theorem 6. Suppose P (λ) has p distinct eigenvalues {λ i } p i=1 , with n ≤ p ≤ ℓn, and that the corresponding set of p eigenvectors {v i } p i=1 satisfies the Haar condition (every subset of n of them is linearly independent). Then there are at least p n different solvents of P (λ), and exactly this many if p = ℓn, which are given by

S = W diag(µ i )W -1 , W = w 1 • • • w n ,
where the eigenpairs (µ i , w i ) n i=1 are chosen among the eigenpairs (λ i , v i ) p i=1 of P .

Note that if we have that p = n in Theorem 6, the distinctness of the eigenvalues is not needed, and then we have a sufficient condition for the existence of a solvent.

Corollary 1. If P (λ) has n linearly independent eigenvectors v 1 , v 2 , . . . , v n , then P (S) has a solvent.

An example which illustrates this last result is the following. Consider the quadratic matrix solvent problem (see [START_REF] Dennis | The algebraic theory of matrix polynomials[END_REF], [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF])

Q(S) = S 2 + -1 -6 2 -9 S + 0 12 - 2 14 .
Q(λ) has eigenpairs (λ i , v i ):

(λ 1 , v 1 ) = 1, 1 0 , (λ 2 , v 2 ) = 2, 0 1 , (λ 3 , v 3 ) = 3, 1 1 , (λ 4 , v 4 ) = 4, 1 1 .
Consider the subsets of eigenvectors:

{v 1 , v 2 }, {v 1 , v 3 }, {v 1 , v 4 }, {v 2 , v 3 } and {v 2 , v 4 }.
Each subset consists of vectors that are linearly independent. Therefore, the complete set of solvents is:

1 0 0 2 , 1 2 0 3 , 3 0 1 2 , 1 3 0 4 and 4 0 2 2 .
Note that we cannot construct a solvent whose eigenvalues are 3 and 4 because the associated eigenvectors are linearly dependent.

Condition Number and Backward Error for the

Matrix Solvent Problem

An analysis and a computable formulation for the condition number and backward error of a solvent can be found in [START_REF] Higham | Numerical analysis of a quadratic matrix equation[END_REF] and [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF]. These results are given for the case of the quadratic solvent problem Q(S). This section presents computable expressions for the condition number and backward error for the general matrix solvent problem (3.1). We follow the ideas presented in [START_REF] Higham | Numerical analysis of a quadratic matrix equation[END_REF], [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF] and [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF].

Condition Number of P (S)

We perform here a similar analysis as we did in Section 2.6.1. For self-consistency of the chapter, we will add all the details of the calculations.

A normwise condition number of the solvent S can be defined by:

κ(S) = lim sup ǫ→0 1 ǫ ||∆S|| F ||S|| F : (P + ∆P )(S + ∆S) = 0, ∆A i F ≤ ǫα i , i = 0 : ℓ , (3.3)
where ∆P (λ) = ℓ i=0 λ i ∆A i . The α i are nonnegative weights; in particular, ∆A i can be forced to zero by setting α i = 0.

Theorem 7. The normwise condition number of the solvent S is given by:

κ(S) = B-1 S BA 2 ||S|| F , (3.4)
where

BS = ℓ j=1 j-1 i=0 (S j-i-1 ) T ⊗ A j S i ,
and

BA = α ℓ (S ℓ ) T ⊗ I n α ℓ-1 (S ℓ-1 ) T ⊗ I n • • • α 0 I n 2 .
Proof. By expanding the first constraint in (3.3) and keeping only the first order terms, we get:

ℓ j=0 ∆A j S j + ℓ j=1 A j DS j (∆S) = O(ǫ 2 ), (3.5)
where DS j denotes the Fréchet derivative (2.8). Using on equation (3.5) the vec operator (see A.2), we obtain:

• vec (∆P (S)) = vec   ℓ j=0 ∆A j S j   = ℓ j=0 vec ∆A j S j = ℓ j=0 (S j ) T ⊗ I n vec (∆A j ) = = α ℓ (S ℓ ) T ⊗ I n α ℓ-1 (S ℓ-1 ) T ⊗ I n • • • α 0 I n 2       vec (∆A ℓ )/α ℓ vec (∆A ℓ-1 )/α ℓ-1 . . . vec (∆A 0 )/α 0       = BA vec (∆A),
• vec

  ℓ j=1 A j DS j (∆S)   = vec   ℓ j=1 A j j-1 i=0 S i ∆SS j-i-1   = ℓ j=1 j-1 i=0 vec A j S i ∆SS j-i-1 = = ℓ j=1 j-1 i=0 (S j-i-1 ) T ⊗ A j S i vec (∆S) = BS vec (∆S).
Then, we have that:

BS vec (∆S) = -BA x + O(ǫ 2 ),
where

x =     vec (∆A ℓ ) . . . vec (∆A 0 )    
Using that ||vec (∆S)|| 2 = ||∆S|| F , we have that the definition (3.3) is equivalent to the following lim sup

ǫ→0 1 ǫ ||∆S|| F ||S|| F : BS vec (∆S) = -BA x + O(ǫ 2 ), x 2 ≤ ǫ = B-1 S BA 2 ||S|| F .
Example 9. Consider the quadratic matrix polynomial (see [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF]):

Q(λ) = λ 2 A 2 + λA 1 + A 0 = λ 2 1 0 0 1 + λ 0 0 1 0 + -1 0 -1 0 ,
with eigenvalues: -1, 0, 0, 1. Note that there are three solvents for Q(S):

S 1 = 1 -1 0 -1 , S 2 = 1 0 0 0 , S 3 = -1 0 -2 0 .
Using the parameters α 2 = A 2 F , α 1 = A 1 F and α 0 = A 0 F , the solvent S 1 has condition number: κ(S 1 ) = 3.63971.

Note that this condition number κ(S 1 ) is equal to the one found in [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF].

Remark 8. For the solvents S 2 and S 3 , the matrix BS in Theorem 7 is singular. Then we have that κ(S 2 ) = ∞ and κ(S 3 ) = ∞.

Backward Error of P (S)

We proceed as in Section 2.6.2 and we obtain bounds for the backward error of P (S). Let α i , for i = 0, . . . , ℓ, be nonnegative weights as in Section 2.6.1. The backward error of a computed solution S ∈ C n×n to (3.1) can be defined as:

η( S) = min{ǫ : (P + ∆P )( S) = 0, ∆A i F ≤ ǫα i , i = 0, . . . , ℓ} (3.6)
By expanding the first constraint in (3.6) we get:

-

P ( S) = ∆A ℓ Sℓ + • • • + ∆A 0 . (3.7)
Then, we have

-P ( S) = α -1 ℓ ∆A ℓ α -1 ℓ-1 ∆A ℓ-1 . . . α -1 1 ∆A 1 α -1 0 ∆A 0          α ℓ Sℓ α ℓ-1 Sℓ-1 . . . α 1 S α 0         
Taking the Frobenius norm, we obtain the lower bound for the backward error:

η( S) ≥ P ( S) F (α 2 ℓ Sℓ 2 F + α 2 ℓ-1 Sℓ-1 2 F + • • • + α 2 1 S 2 F + α 2 0 ) 1/2
. Consider again equation (3.7). Using the vec operator (see A.2), we obtain:

-vec (P ( S)) = (( Sℓ ) T ⊗ I n ) vec (∆A ℓ ) + • • • + ( ST ⊗ I n ) vec (∆A 1 ) + I n 2 vec (∆A 0 ) = α ℓ ( Sℓ ) T ⊗ I n . . . α 1 ST ⊗ I n α 0 I n 2       vec (∆A ℓ )/α ℓ . . . vec (∆A 1 )/α 1 vec (∆A 0 )/α 0       ,
which can be written as:

Hz = r, H ∈ C n 2 ×(ℓ+1)n 2 (3.8)
Here, we assume that H has full rank to guarantee that (3.8) has a solution (backward error is finite). Then, the backward error is the minimum 2-norm solution to the least square problem defined from (3.8), that is:

η(T ) = H + r 2 . (3.9)
where H + denotes the Moore-Penrose pseudoinverse of H + . Eq. (3.9) yields an upper bound for η( S):

η( S) ≤ H + 2 r 2 = r 2 σ min (H)
,

where σ min denotes the smallest singular value, which is nonzero by assumption. Note that:

σ min (H) 2 = λ min (HH * ) = λ min (α 2 ℓ ( Sℓ ) T Sℓ ⊗ I n + • • • + α 2 1 ( S) T S ⊗ I n + α 2 0 I n 2 ) = λ min (α 2 ℓ ( Sℓ ) * Sℓ ⊗ I n + • • • + α 2 1 ( S) * S ⊗ I n + α 2 0 I n 2 ) ≥ α 2 ℓ σ min ( Sℓ ) 2 + • • • + α 2 1 σ min ( S) 2 + α 2 0 .
Thus we obtain the upper bound for η( S):

η( S) ≤ P ( S) F (α 2 ℓ σ min ( Sℓ ) 2 + • • • + α 2 1 σ min ( S) 2 + α 2 0 ) 1/2 .

Computation of Solvents

The question of designing symbolic algorithms for computing solvents remains relatively unexplored. Attempts have been made to formulate the problem as a system of polynomial equations, which can be solved via standard methods. However, this approach becomes cumbersome for problems of large size (see [START_REF] Higham | Numerical analysis of a quadratic matrix equation[END_REF]). Motivated by applications to differential equations [START_REF] Brugnano | Solving Differential Equations by Multistep Initial and Boundary Value Methods[END_REF], we study an approach to the symbolic or symbolic-numeric computation of solvents based on the moment method, by specializing the results presented in Section 4.1.

A Computational Approach

Our approach to compute matrix solvents is based on the relation between the matrix solvent problem (3.1) and the invariant pair problem (2.3). We state this in the following result. Proof. As (Y, T ) ∈ C n×n × C n×n is an invariant pair of P (λ), we have:

A ℓ Y T ℓ + • • • + A 2 Y T 2 + A 1 Y T + A 0 Y = 0.
Since Y is invertible, we can post-multiply by Y -1 . Then we get:

A ℓ Y T ℓ Y -1 + • • • + A 2 Y T 2 Y -1 + A 1 Y T Y -1 + A 0 = 0 ⇔ A ℓ S ℓ + • • • + A 2 S 2 + A 1 S + A 0 = 0.
Therefore, S is a matrix solvent of P (λ).

Example 10. Consider again the quadratic matrix polynomial of previous section:

Q(λ) = λ 2 1 0 0 1 + λ -1 -6 2 -9 + 0 12 -2 14 .
Invariant pairs (X i , T i ) for P (λ) are given by:

X 1 = 1 1 1 1 , T 1 = 3 0 0 4 ; X 2 = 1 1 1 2 , T 2 = 0 -2 1 3
;

X 3 = 1 3 2 5 , T 3 = 0 -6 1 5
;

X 4 = 3 11 4 12 , T 4 = 0 -3 1 4
;

X 5 = -1 2 1 4 , T 5 = 0 -4 1 5 ; X 6 = 1 4 2 6 , T 6 = 0 -8 1 6
.

Taking S i = X i T i Y -1 i , for i = 1, . . . , 6
, the S i satisfying (3.1) are:

S 2 = 1 0 0 2 , S 3 = 3 0 1 2 , S 4 = 1 2 0 3 , S 5 = 1 3 0 4 and S 6 = 4 0 2 2 .
Note that the matrix X 1 is singular and then we can't construct the solvent S 1 associated to the eigenvalues 3 and 4 (the associated eigenvectors are linearly dependent).

Matrix p-th Root

Consider the matrix equation

X p -A = 0, (3.10)
where A ∈ C n×n . A matrix X satisfying (3.10) is called a p-th root of A (see, e.g., [START_REF] Guo | On Newton's method and Halley's method for the principal pth root of a matrix[END_REF], [START_REF] Guo | A Schur-Newton method for the matrix p-th root and its inverse[END_REF], [START_REF] Higham | On pth roots of stochastic matrices[END_REF], [START_REF] Iannazzo | On the Newton method for the matrix pth root[END_REF], [START_REF] Smith | A Schur algorithm for computing matrix pth roots[END_REF]). One application of p-th roots is in the computation of the matrix logarithm through the relation (see [START_REF] Cheng | Approximating the logarithm of a matrix to specified accuracy[END_REF])

logA = p logA 1 p .
When p = 2 we have the matrix equation:

X 2 -A = 0. (3.11)
A solution X of (3.11) is called a matrix square root of A (see, e.g., [START_REF] Higham | Newton's method for the matrix square root[END_REF], [START_REF] Higham | Computing real square roots of a real matrix[END_REF], [START_REF] Higham | Stable iterations for the matrix square root[END_REF], [START_REF] Iannazzo | A note on computing the matrix square root[END_REF], [START_REF] Meini | The matrix square root from a new functional perspective: theoretical results and computational issues[END_REF]). In this case, if A is singular, the existence of a square root depends on the Jordan structure of the zero eigenvalues (see [START_REF] Cross | Square roots of complex matrices[END_REF]). For instance, the matrix:

A = 0 1 0 0
has no square root, while the matrix

A =    0 1 0 0 0 0 0 0 0   
does have a square root, specifically:

X =    0 0 1 0 0 0 0 1 0    .
If A is nonsingular, it always has a p-th root. The number of square roots varies from two (for a nonsingular Jordan block) to infinity (any involutory matrix is a square root of the identity matrix). Moreover, when A has no non-positive real eigenvalues, one can define the notion of principal root, denoted by A 1 2 . The method of reference for computing matrix square roots is to compute a Schur decomposition, compute a square root of the triangular factor, and then transform back (see [START_REF] Björck | A Schur method for the square root of a matrix[END_REF] and [START_REF] Higham | Computing real square roots of a real matrix[END_REF]).

Iterative methods are alternatives for computing p-th roots. In [START_REF] Higham | Newton's method for the matrix square root[END_REF], [START_REF] Higham | Stable iterations for the matrix square root[END_REF], [START_REF] Iannazzo | A note on computing the matrix square root[END_REF], [START_REF] Iannazzo | On the Newton method for the matrix pth root[END_REF], [START_REF] Meini | The matrix square root from a new functional perspective: theoretical results and computational issues[END_REF], for instance, several variations of the Newton's method to approximate the (square roots) p-th roots are presented.

Note that equation (3.11) is a particular case of the monic matrix solvent problem of degree p:

P (S) := S p + A p-1 S p-1 + • • • + A 2 S 2 + A 1 S + A 0 = 0,
when the matrices A i are zero for i = 1, . . . , p -1. Therefore, we can use our (block) moment method to construct a solvent S satisfying (3.10) (see Section 4.3.2).

Example 11. Consider the matrix equation (see [START_REF] Smith | A Schur algorithm for computing matrix pth roots[END_REF]):

X 4 = A,
where:

A =       1 -1 -1 -1 0 1.3 -1 -1 0 0 1.7 -1 0 0 0 2       .
Consider the circle: Γ = 1 + 1 4 e ıt . Applying our moment method, which will be presented in Chapter 4, we construct the matrix solvent X:

X =       1 -0.2259665745747 -0.2609342676468 -0.3057660919094 0 1.067789972372 -0.1851709326575 -0.212512633424 0 0 1.141858345435 -0.157829231891 0 0 0 1.189207115003       .
The residual is:

res( X) = X4 -A 2 A 2 = 3.1024e -13.

Solvents and Triangularized Matrix Polynomials

Motivated by the results in [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF] and [START_REF] Tisseur | Triangularizing quadratic matrix polynomials[END_REF], where the authors analyze a method for triangularizing the matrix polynomial P (λ), we aim here to study the relation between solvents of general and triangularized matrix polynomials.

Triangularizing Matrix Polynomials

For any algebraically closed field F, any matrix polynomial P (λ) ∈ F[λ] n×m , with n ≤ m, can be reduced to triangular form via unimodular transformations. In other words, there exist matrix polynomials E(λ) and F (λ) with nonzero constant determinant such that:

T (λ) = E(λ)P (λ)F (λ),
where T (λ) is monic triangular and preserves the degree and the finite and infinite elementary divisors of P (λ) [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF], [START_REF] Tisseur | Triangularizing quadratic matrix polynomials[END_REF].

Theorem 9. [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF] For an algebraically closed field F, any P (λ) ∈ F[λ] n×m with n ≤ m is triangularizable.

A continuation, we recall the procedure to triangularize a quadratic matrix polynomial Q(λ) (see [START_REF] Taslaman | Triangularizing matrix polynomials[END_REF]).

3.4.2 Procedure to Triangularize a Quadratic Matrix Polynomial Theorem 1.7 in [START_REF] Gohberg | Matrix polynomials[END_REF] shows that any complex matrix polynomial of degree ℓ, with nonsingular leading coefficient, is equivalent to a monic triangular matrix polynomial of the same degree. In this section, we recall the procedure for ℓ = 2, i.e., for the quadratic matrix polynomial Q(λ).

Compute the invariant factors

: α 1 (λ) | α 2 (λ) | • • • | α n (λ) of Q(λ) and consider D(λ) = diag(α 1 (λ), . . . , α n (λ)). 2. If deg(α 1 ) = 2, then deg(α i ) = 2, i = 2 : n. Then D(λ)
is a monic triangular quadratic matrix polynomial equivalent to Q(λ) and the construction is done.

Otherwise, go to Step 3.

3.

If ℓ 1 = deg(α 1 ) < 2, then ℓ n = deg(α n ) > 2 and there is a monic polynomial s(λ) of degree ℓ n -2 such that α j-1 (λ) | α 1 (λ)s(λ) | α j , for some index j, 1 < j ≤ n.

The s(λ) is the product of α j-1 (λ)/α 1 (λ) and some of the factors in the prime factorization of α j (λ)/α j-1 (λ).

4. Perform the following elementary transformations on D(λ):

(i) Add to column n the first column multiplied by s(λ). Then add to row n the first row multiplied by -α n (λ)/(α 1 (λ)s(λ)) and permute columns 1 and n.

(ii) Successively interchange rows k and k + 1 for k = 1, 2, . . . , j -2, so that rows 1, 2, . . . , j -2, j -1 of the new matrix are rows 2, 3, . . . , j -1 and 1, respectively, of the former one.

(iii) Permute columns 1 to j -1 in the same way as the rows in (ii).

The resulting matrix polynomial has the form:

T 1 (λ) =                  α 2 (λ) . . . α j-1 (λ) α 1 (λ)s(λ) • • • • • • α 1 (λ) α j (λ) . . . . . . α n-1 (λ) . . . -α n (λ)/s(λ)                  . 5. Let D 1 (λ) = diag(α 2 (λ), . . . , α j-1 (λ), α 1 (λ)s(λ), . . . , α n-1 (λ)). If D 1 (λ) is diagonal of degree 2
, the construction is done; otherwise, use D 1 (λ) to choose the new index j and polynomials s(λ) as in Step 3 and perform the elementary transformations described in Step 4 on the whole matrix T 1 (λ). Repeat these steps until all diagonal entries have degree 2.

What can we say about solvents for a given matrix polynomial P (λ) and for the associated triangularized polynomial? A partial answer will be given in Theorem 11.

Theorem 10. For any ℓn × ℓn monic linearization λI -A of P (λ) ∈ C[λ] n×n with nonsingular leading coefficient, there exists U ∈ C n×ℓn with orthogonal columns such that

M =       U U A . . . U A ℓ-1      
is nonsingular and λI -M AM -1 is a linearization for the polynomial

T (λ) = λ ℓ I + λ ℓ-1 T ℓ-1 + • • • + λ 2 T 2 + λT 1 + T 0 ,
which is upper triangular and equivalent to P (λ).

Theorem 10 is a straightforward generalization of a result found in [START_REF] Tisseur | Triangularizing quadratic matrix polynomials[END_REF]. Note that, for the time being, we assume that the leading coefficient A ℓ is nonsingular.

Theorem 11. Let P (λ) be a n × n matrix polynomial and consider the linearization:

A =          0 I n 0 • • • 0 0 0 I n • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n -A 0 -A 1 -A 2 • • • -A ℓ-1          .
Let M be as in Theorem 10 and let Y 1 be the first n × n block of the matrix

M -1       I n S t . . . S ℓ-1 t      
. If Y 1 is nonsingular and S t is a solvent for the triangularized problem, i.e.,

T (S t ) = 0, then S = Y 1 S t Y -1 1
is a solvent for P (S).

Proof. Note that:

       0 0 . . . 0        =        S t -S t S 2 t -S 2 t . . . -T 0 -T 1 S t -T 2 S 2 t -• • • -S ℓ t        = =          0 I n 0 • • • 0 0 0 I n • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n -T 0 -T 1 -T 2 • • • -T ℓ-1                 I n S t . . . S ℓ-1 t        -        S t S 2 t . . . S ℓ t        = = M -1          0 I n 0 • • • 0 0 0 I n • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n -T 0 -T 1 -T 2 • • • -T ℓ-1          M M -1        I n S t . . . S ℓ-1 t        -M -1        I n S t . . . S ℓ-1 t        S t (iii) = =          0 I n 0 • • • 0 0 0 I n • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n -A 0 -A 1 -A 2 • • • -A ℓ-1          M -1        I n S t . . . S ℓ-1 t        -M -1        I n S t . . . S ℓ-1 t        S t . Since M -1       I n S t . . . S ℓ-1 t       has size ℓn × n, let us partition it as       Y 1 Y 2 . . . Y ℓ      
, where Y i ∈ C n×n for i = 1, . . . , ℓ. Then:

         0 I n 0 • • • 0 0 0 I n • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n -A 0 -A 1 -A 2 • • • -A ℓ-1                Y 1 Y 2 . . . Y ℓ       -       Y 1 Y 2 . . . Y ℓ       S t = 0 0 .
Then we have:

Y i = Y 1 S i-1 t , for i = 2, . . . , ℓ; (3.12) -A 0 Y 1 -A 1 Y 2 -• • • -A ℓ-1 Y ℓ -Y ℓ S t = 0. (3.13)
Substituting equations (3.12) in (3.13) we obtain:

0 = Y 1 S ℓ t + A ℓ-1 Y 1 S ℓ-1 t + • • • + A 1 Y 1 S t + A 0 Y 1 .
If Y 1 is invertible we have:

0 = Y 1 S ℓ t Y -1 1 + A ℓ-1 Y 1 S ℓ-1 t Y -1 1 + • • • + A 1 Y 1 S t Y -1 1 + A 0 . Taking S = Y 1 S t Y -1 1
we have:

0 = S ℓ + A ℓ-1 S ℓ-1 + • • • + A 2 S 2 + A 1 S + A 0 := P (S). Then S = Y 1 S t Y -1 1
is a solvent for P (S).

Example 12. Consider the matrix polynomial:

P (λ) = λ 2 1 0 0 1 + λ -1 -6 2 -9 + 0 12 -2 14 .
Let us triangularize P (λ) using the procedure in Section 3.4.2. One possible result is:

T (λ) = (λ -1)(λ -2) 1 0 (λ -3)(λ -4) = λ 2 -3λ + 2 1 0 λ 2 -7λ + 12 = = λ 2 1 0 0 1 + λ -3 0 0 -7 + 2 1 0 12 = λ 2 I 2 + λT 1 + T 0 .
Now, suppose that the solvent S t ∈ C 2×2 of the triangularized problem is in upper trian-gular form, i.e.:

S t = a b 0 d .
Then, we have:

T (S t ) = S 2 t + T 1 S t + T 0 = a 2 -3a + 2 ab -3b + bd + 1 0 d 2 -7d + 12 .
In the task of solving the problem T (S t ) = 0, we see that:

     a 2 -3a + 2 = 0 ab -3b + bd + 1 = 0 d 2 -7d + 12 = 0 =⇒ a = 1 or a = 2 d = 3 or d = 4
Then, we have four cases:

Case 1: If a = 1 and d = 3, then b = -1. We have:

S t 1 = 1 -1 0 3 .
Case 2: If a = 1 and d = 4, then b = -1 2 . We have:

S t 2 = 1 -1 2 0 4 .
Case 3: If a = 2 and d = 3, then b = -1 2 . We have:

S t 3 = 2 -1 2 0 3 .
Case 4: If a = 2 and d = 4, then b = -1 3 . We have:

S t 4 = 2 -1 3 0 4 .
Having the solvents S t i for the triangularized problem, we would like to find the (associated) solvents S i for the original problem P (λ).

Consider the linearization of P (λ):

A = 0 I 2 -A 0 -A 1 =       0 0 1 0 0 0 0 1 0 -12 1 6 2 -14 -2 9       ,
and the matrix M :

M = U U A =       -20 15 8 -7 -2 0 2 0 -14 2 2 0 0 -24 0 12       . Now, compute M -1 I 2 S t i for i = 1, . . . , 4.
Taking Y i as the first 2 × 2 block of this result and computing S i = Y i S t i Y -1 i , we have:

S 1 = 1 2 0 3 , S 2 = 1 3 0 4 , S 3 = 3 0 1 2 and S 4 = 4 0 2 2 .
By Theorem 11, the S i are solvents for P (S), i.e., they satisfy P (S i ) := S 2 i +A 1 S i +A 0 = 0 for i = 1, . . . , 4.

Note that in the description of the procedure to triangularize a quadratic matrix polynomial in Section 3.4.2, we have the freedom when choosing the polynomial s(λ). Therefore, the final triangularization for the matrix polynomial will change depending on our choice. Moreover, the number of solvents may also change.

To demonstrate this, consider the following example.

Example 13. Consider again the matrix polynomial P (λ) of Example 12, but now with its different triangularization:

T (λ) = (λ -1)(λ -3) 1 0 (λ -2)(λ -4) = λ 2 -4λ + 3 1 0 λ 2 -6λ + 8 = = λ 2 1 0 0 1 + λ -4 0 0 -6 + 3 1 0 8 = λ 2 I 2 + λT 1 + T 0 .
Suppose that the solvent S t ∈ C 2×2 of the triangularized problem is in upper triangular form, i.e.:

S t = a b 0 d .
Then, we have:

T (S t ) = S 2 t + T 1 S t + T 0 = a 2 -4a + 3 ab + bd -4b + 1 0 d 2 -6d + 8 .
When solving the problem T (S t ) = 0, we see that:

     a 2 -4a + 3 = 0 ab + bd -4b + 1 = 0 d 2 -6d + 8 = 0 =⇒ a = 1 or a = 3 d = 2 or d = 4
Then, we have four cases:

Case 1: If a = 1 and d = 2, then b = 1. We have:

S t 1 = 1 1 0 2 .
Case 2: If a = 1 and d = 4, then b = -1. We have:

S t 2 = 1 -1 0 4 .
Case 3: If a = 3 and d = 2, then b = -1. We have:

S t 3 = 3 -1 0 2 .
Case 4: If a = 3 and d = 4, then b = -1 3 . We have:

S t 4 = 3 -1 3 0 4 .
Having the solvents S t i for the triangularized problem, we would like to find the (associated) solvents S i for the original problem P (λ).

Consider the linearization of P (λ):

A = 0 I 2 -A 0 -A 1 =       0 0 1 0 0 0 0 1 0 -12 1 6 2 -14 -2 9       ,
and the matrix M : Now, compute M -1 I 2 S t i for i = 1, . . . , 4. Taking Y i as the first 2 × 2 block of this result and computing S i = Y i S t i Y -1 i , we have:

M = U U A =      
S 1 = 1 0 0 2 , S 2 = 1 3 0 4 , S 3 = 3 0 1 2 .
Note that the matrix Y 4 is equal to:

Y 4 = 1 4 3 1 4 3 .
This matrix is singular and then we can't construct the solvent S 4 associated to the eigenvalues 3 and 4 (the associated eigenvectors are linearly dependent). By Theorem 11, the S i are solvents for P (S), i.e., they satisfy P (S i ) := S 2 i + A 1 S i + A 0 = 0 for i = 1, 2, 3.

A Problem with an Infinite Number of Solvents

What happens to the ideas outlined above when working on problems with an infinite number of solvents? Here is an example taken from [START_REF] Pereira | On solvents of matrix polynomials[END_REF].

Consider the matrix polynomial: One possible triangularization for P (λ) is:

P (λ) = λ 2 I + λ    -7 -2 -2 3
T (λ) =    (λ -3)(λ -4) (λ -3) 0 0 (λ -3) 2 1 0 0 (λ -4) 2    = λ 2 I 2 + λT 1 + T 0 = = λ 2 I + λ    -7 1 0 0 -6 0 0 0 -8    +    12 -3 0 0 9 1 0 0 16    .
Now, suppose that the solvent S t ∈ C 3×3 of the triangularized problem is in upper triangular form, i.e.:

S t =    x 11 x 12 x 13 0 x 22 x 23 0 0 x 33    .
Then, we have:

T (St) =S 2 t + T 1 St + T 0 = =    (x 11 -3)(x 11 -4) x 22 -7x 12 + x 11 x 12 + x 12 x 22 -3 x 23 -7x 13 + x 11 x 13 + x 12 x 23 + x 13 x 33 0 (x 22 -3) 2 x 22 x 23 -6x 23 + x 23 x 33 + 1 0 0 (x 33 -4) 2    .
In the task of solving the problem T (S t ) = 0, we see that: Then we find that x 23 = -1, x 12 = 0 and x 12 = -1, which is a contradiction. In this case, there is no solution and then we can't construct a solvent.

II. If x 11 = 4, x 22 = 3 and x 33 = 4: Then we find that x 23 = -1 and x 13 = x 12 + 1. In this case, the solvent S t has the form:

S t =    4 x 12 x 12 + 1 0 3 -1 0 0 4    =    4 0 1 0 3 -1 0 0 4    + x 12    0 1 1 0 0 0 0 0 0    , for x 12 ∈ C.
Thus T (λ) has an infinite number of solvents and the same holds for P (λ).

Chapter 4 :

Moments, Hankel Pencils and Computation of Invariant Pairs

Introduction

Numerical methods based on contour integrals for the computation of eigenvalues of matrix polynomials and analytic matrix-valued functions have recently met with growing interest. Such techniques are related to the well-known method of moments, where the moments are computed by numerical quadrature (see [START_REF] Golub | Matrices, moments and quadrature with applications[END_REF]).

In this section, we explore a similar approach for computing invariant pairs. Our main reference is the Sakurai-Sugiura method (see [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] and [START_REF] Sakurai | A projection method for generalized eigenvalue problems using numerical integration[END_REF]) as well as the presentation given in [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF].

Toeplitz and Hankel Matrices

A Toeplitz matrix is a matrix where each descending diagonal from left to right is constant (see [START_REF] Golub | Matrix computations[END_REF]). Formally: Definition 13. An n × n matrix T is Toeplitz if there exist scalars r -n+1 , . . . , r 0 , . . . , r n-1 such that [T ] ij = r i-j for 1 ≤ i, j ≤ n.

Then, an n × n Toeplitz matrix has the form:

T =             r 0 r -1 r -2 . . . . . . r -n+1
. . . . . . r 1 r 0 r -1 r n-1 . . . . . . r 2 r 1 r 0             .
The Toeplitz structure can be exploited to design fast algorithms for fundamental problems such as solving linear systems. For instance, when T is symmetric and positive definite, the Toeplitz system T x = b can be solved in O(n 2 ) flops. This is due to the structure of the inverse T -1 . Definition 14. An n × n matrix H is Hankel if there exist scalars h 0 , . . . , h 2n-2 such that [H] i+1,j+1 = h i+j for i, j = 0, 1, . . . , n -1.

Then, an n × n Hankel matrix has the form:

H =       h 0 h 1 • • • h n-1 h 1 h 2 • • • h n . . . . . . . . . h n-1 h n • • • h 2n-2       .
Define the "flip matrix":

E n =          0 0 • • • 0 1 0 0 • • • 1 0 . . . . . . . . . . . . 0 1 • • • 0 0 1 0 • • • 0 0          .
Then

The Moment Method and Eigenvalues

Let us begin by briefly recalling a few basic facts about the Sakurai-Sugiura moment method. Here, we essentially follow the presentation given in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF].

Let P (λ) be an n × n matrix polynomial and let Γ be a closed contour in the complex plane and let u and v be arbitrarily given vectors in C n . Define the function:

f (λ) := u H P (λ) -1 v.
In the following, it will be understood that no eigenvalue of P (λ) should lie exactly on the contour Γ: each eigenvalue should be either inside or outside the contour.

The next theorem, which can be found in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF], gives a representation for f (λ) that will be useful later on. Theorem 12. [Thm. 3.1, [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF]] Let D(λ) = diag (d 1 (λ), . . . , d n (λ)) be the Smith form of P (λ), and let E(λ) and F (λ) be as in

(1.7). Let χ j (λ) = u H q j (λ)p j (λ) H v, 1 ≤ j ≤ n. Then f (λ) = n j=1 χ j (λ) d j (λ) , (4.1)
where q j (λ) and p j (λ) are the column vectors of E(λ) and F (λ) H , respectively.

Definition 15. Let k ∈ N. The k-th moment of f (z) is:

µ k = 1 2πı Γ z k f (z)dz. (4.2)
For a positive integer m, define the Hankel matrices H 0 , H 1 ∈ C m×m as follows:

H 0 =       µ 0 µ 1 • • • µ m-1 µ 1 µ 2 • • • µ m . . . . . . . . . µ m-1 µ m • • • µ 2m-2       , H 1 =       µ 1 µ 2 • • • µ m µ 2 µ 3 • • • µ m+1 . . . . . . . . . µ m µ m+1 • • • µ 2m-1       . (4.
3)

The eigenvalue algorithm presented in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] relies on the following result:

Theorem 13. [Thm. 3.4, [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF]] Suppose that the polynomial P (λ) has exactly m eigenvalues λ 1 , . . . , λ m in the interior of Γ, and that these eigenvalues are distinct, simple and non degenerated. If χ n (λ ℓ ) = 0 for 1 ≤ ℓ ≤ m, then the eigenvalues of the pencil H 1 -λH 0 are given by λ 1 , . . . , λ m .

So, in order to approximate λ 1 , . . . , λ m , it suffices to compute by quadrature the first 2m moments of f (λ) and then apply an eigensolver, such as the QZ method (see A.5), to the resulting Hankel pencil. Block versions of this approach have also been proposed; we will say more on this later.

We now wish to investigate the behavior of the above method when the hypothesis that the λ i 's are distinct is removed. In particular, we aim to generalize Theorem 13, which is based on the Vandermonde factorization of H 0 and H 1 . Theorem 14. Suppose that P (λ) has exactly m eigenvalues in the interior of Γ, namely, distinct eigenvalues λ 0 , . . . , λ s with algebraic multiplicities m 0 , . . ., m s , respectively, such that m = m 0 + • • • + m s . Moreover, assume that no eigenvalue of P (λ) lies exactly on the contour Γ. If the geometric multiplicity of the λ i 's, for i = 0, . . . , s, is equal to one, then the matrix H 0 is nonsingular and eigenvalues of the pencil H 1 -λH 0 are given by λ 0 , . . . , λ s with algebraic multiplicities m 0 , . . . , m s .

Proof. Suppose that, in the Smith form (1.7) of P (λ), the matrix D(λ) is in the form

D(λ) = diag(d 1 (λ), . . . , d n-1 (λ), d n (λ)), where d n (λ) = (λ -λ 0 ) m 0 (λ -λ 1 ) m 1 • • • (λ -λ s ) ms r i=s+1 (λ -λ i ) m i ,
and λ s+1 , . . . , λ r are the eigenvalues of P (λ) located outside the contour Γ, with algebraic multiplicities m s+1 , . . . , m r . Moreover, define dn (λ) = s i=0 (λ -λ i ) m i , i.e., dn (λ) is the factor of d n (λ) whose roots are the eigenvalues of P (λ) located inside Γ.

Since the geometric multiplicities of the λ i 's are all equal to one, the factors (λ -λ i ), for i = 0, . . . , s, do not appear in the monic scalar polynomials d 0 (λ), . . . , d n-1 (λ).

Applying Theorem 12, we have:

µ k = 1 2πı Γ z k f (z)dz = 1 2πı Γ n j=1 χ j (z) d j (z) z k dz = = 1 2πı Γ ϕ(z) d n (z) z k dz,
where

ϕ(z) = n j=1 χ j (z)h j (z), with h j (z) = d n (z) d j (z) .
We can introduce partial fraction decompositions and write

µ k = 1 2πı Γ ϕ(z) d n (z) z k dz = = 1 2πı Γ m 0 i=1 c 0,i z k (z -λ 0 ) i + • • • + ms i=1 c s,i z k (z -λ s ) i dz = = s j=0 m j i=1 1 2πı Γ c j,i z k (z -λ j ) i dz,
where c j,i ∈ C, for j = 0, . . . , s and i = 1, . . . , m j . Classical results on residues then yield

µ k = s j=0 m j i=1 c j,i Res z k (z -λ j ) i , λ j = = s j=0 m j i=1 c j,i 1 (i -1)! lim z→λ j d i-1 dz i-1 (z -λ j ) i z k (z -λ j ) i = = s j=0 m j i=1 c j,i 1 (i -1)! lim z→λ j d i-1 dz i-1 z k = = s j=0 m j i=1 ν j,i λ k-i+1 j , (4.4)
where

ν j,i = c j,i (i-1)! (k -i + 2)(k -i + 3) • • • k if k ≥ i -1, 0 otherwise.
Now, consider the pencil H 1 -λH 0 with H 0 and H 1 defined as in (4.3). Because of (4.4), and of the fact that λ 0 , . . . , λ s are roots of dn (λ), the moments µ k satisfy a linear recurrence equation of the form:

µ k = a m-1 µ k-1 + a m-2 µ k-2 + • • • + a 0 µ k-m . (4.5)
Moreover, dn (λ) is the polynomial of smallest degree that has roots λ 0 , . . . , λ s with the prescribed multiplicities m 0 , . . . , m s , so the recurrence (4.5) has the shortest possible length; also, note that the a i 's in (4.5) are actually the coefficients of dn (λ). Therefore, the matrices H 0 and H 1 have full rank. The same argument shows that H 0 and H 1 are rank-deficient if taken of size larger than m × m (see [START_REF] Meurant | Matrices, moments and quadrature with applications[END_REF], [ [START_REF] Gantmacher | Theory of matrices I[END_REF], Vol. 2, pp. 205)] and [START_REF] Boley | A general Vandermonde factorization of a Hankel matrix[END_REF]). As a consequence of the shifted Hankel form of H 0 and H 1 , we have

H 0 C = H 1 , where C is a matrix in companion form C =          0 0 • • • 0 x 0 1 0 • • • 0 x 1 0 1 • • • 0 x 2 . . . . . . . . . . . . . . . 0 0 • • • 1 x m-1          ,
and its last column is given by the solution of the linear system

H 0       x 0 x 1 . . . x m-1       =       µ m µ m+1 . . . µ 2m-1       . (4.6)
The polynomial of degree m:

p(λ) = λ m -x m-1 λ m-1 -• • • -x 0
is a scalar multiple of dn (λ), and its roots are the λ i 's generating the entries of the pencil H 1 -λH 0 . So we also have that the µ i 's satisfy the recurrence (4.5):

µ k = x m-1 µ k-1 + x m-2 µ k-2 + • • • + x 0 µ k-m ,
where k = m, m + 1, . . . , 2m. Consider now the Jordan matrix

J =     J 1 . . . J s     ,
where each block J i , of dimension m i , is a square matrix of the form

J i =        λ i 1 λ i . . . . . . 1 λ i       
, and define the confluent Vandermonde matrix

V = v J T v • • • (J T ) r-1 v , where v T = e [m 1 ]T 1 • • • e [ms]T 1
is partitioned conformally with J and e

[m ℓ ]T 1 = 1 0 • • • 0 T
is the m ℓ -dimensional unit coordinate vector. Then we have

V C = v J T v • • • (J T ) r-1 v C = J T v • • • (J T ) r-1 v -(x 0 I + x 1 J + • • • + x m-1 J r-1 ) T v = J T v • • • (J T ) r-1 v (J T ) r v = J T V,
where we have used the property p(J) = 0, that is, the Cayley-Hamilton theorem.

We can now introduce the Vandermonde decomposition of the Hankel matrices H 0 and H 1 . From the results presented in [START_REF] Boley | A general Vandermonde factorization of a Hankel matrix[END_REF] it follows that there exist block matrices B 0 = diag(D

s ), partitioned conformally with J, satisfying the conditions B 0 J T = JB 0 and B 1 J T = JB 1 , so that

H i = V T B i V, for i = 0, 1.
Moreover, we can prove that JB 0 = B 1 :

H 0 C = H 1 ⇔ (V T B 0 V )C = V T B 1 V ⇔ V T B 0 J T V = V T B 1 V ⇔ V T JB 0 V = V T B 1 V ⇔ JB 0 = B 1 ,
where we used the properties V C = J T V and B i J T = JB i . Therefore, we have:

H 1 -λH 0 = V T B 1 V -λV T B 0 V = V T JB 0 V -λV T B 0 V = V T (J -λI)B 0 V.
So the eigenvalues of H 1 -λH 0 are λ 0 , . . . , λ s with respective multiplicities m 0 , . . . , m s . Example 14. Consider the matrix polynomial:

P (λ) = λ 2 I + λ       -2 0 0 0 0 -1 0 0 0 0 -6 0 0 0 0 -5       +       1 0 0 1 0 1 4 0 0 0 0 9 0 0 0 0 6       . P (λ) has eigenvalues: λ 1 = 1 2 , λ 2 = 1, λ 3 = 2, λ 4 = 3 with algebraic multiplicities: m 1 = 2, m 2 = 2, m 3 = 1, m 4 = 3.
The associated Smith form is:

D(λ) = diag ( 1, 1, 1, (λ - 1 2 ) 2 (λ -1) 2 (λ -2)(λ -3) 3 .
Choosing vectors u = 2 -2 1 -1

T and v = 0 1 0 2 T
, we find that:

H 0 =       -3 -7 -9 -21 2 -7 -9 -21 2 -12 -9 -21 2 -12 -109 8 -21 2 -12 -109 8 -123 8       , H 1 =       -7 -9 -21 2 -12 -9 -21 2 -12 -109 8 -21 2 -12 -109 8 -123 8 -12 -109 8 -123 8
-551 32       .
Then, we have:

C = H -1 0 H 1 =       0 0 0 -1 4 1 0 0 3 2 0 1 0 -13 4 0 0 1 3       .
Note that the eigenvalues of C are 1 2 , 1 2 , 1, 1. Moreover, the companion matrix C is associated with the monic polynomial:

p(λ) = λ 4 -3λ 3 + 13 4 λ 2 - 3 2 λ + 1 4 ,
whose roots are indeed 1 2 , 1 2 , 1, 1. In fact, the Vandermonde factorization for H 0 and H 1 is

H i = V T B i V , i = 0, 1, where V =       1 1 2 1 4 1 8 0 1 1 3 4 1 1 1 1 0 1 2 3       , B 0 =       0 -2 0 0 -2 0 0 0 0 0 -3 -2 0 0 -2 0       , B 1 =       -2 -1 0 0 -1 0 0 0 0 0 -5 -2 0 0 -2 0       ,
and JB 0 = B 1 , with

J =       1 2 1 0 0 0 1 2 0 0 0 0 1 1 0 0 0 1       .
Example 15. Consider the quadratic matrix polynomial:

P (λ) = λ 2    1 0 0 2 1 0 -1 1 -2    + λ    0 0 0 -4 -2 0 2 -2 4    +    -2 1 -2 2 1 0 -1 1 -2    ,
with associated Smith form:

D(λ) = diag ((d 1 (λ), d 2 (λ), d 3 (λ)) =    1 0 0 0 (λ -1) 2 0 0 0 (λ -1) 3 (λ + 1)    .
The Jordan matrix associated with the linearized form of P (λ) is:

J = J 1 0 0 J 2 ,
where each block J i is:

J 1 = 1 1 0 1 , J 2 =       -1 0 0 0 0 1 1 0 0 0 1 1 0 0 0 1       .
Note that we have the eigenvalue λ = 1 in different Jordan blocks.

Choose Γ as the circle ϕ(t) = 1 + 1 10 e ıt , which contains 5 eigenvalues λ = 1, and consider vectors u = 3 1 -2

T and v = 3 -1 -2 T
. Theorem 14 implies that the moment method yields the eigenvalues of P (λ) inside the contour, which are roots of d 3 (λ), i.e., λ 0 = 1 with multiplicity m 0 = 3. Let us compute the matrices H 0 and H 1 of size 3 × 3: The matrix H 0 is nonsingular. Then, we have:

H 0 =    7 3 3 3 3 7 3 7 15    , H 1 =    3
C = H -1 0 H 1 =    0 0 1 1 0 -3 0 1 3    .
Note that eig(C) = 1, 1, 1.

As the contour contains 5 eigenvalues, we might ask what happens when taking the Hankel matrix Ĥ0 of size 5 × 5:

Ĥ0 =        
        .
This matrix is singular, therefore we will not be able to find all the 5 eigenvalues inside the contour. The method miss the additional multiplicities associated with the polynomial d 2 (λ).

Remark 10. We conclude that the scalar moment method can be used to compute the (possibly multiple) eigenvalues of P (λ) that belong to the interior of Γ and that are roots of the polynomial d n (λ). The method misses the additional multiplicities associated with the polynomials d 1 (λ), . . . , d n-1 (λ).

The above remark is consistent with the fact that the Jordan form of a companion matrix only contains one Jordan block for each eigenvalue: it is not possible to capture multiple eigenvalues associated with several Jordan blocks.

In order to "see" the additional eigenvalues that are roots of d 1 (λ), . . . , d n-1 (λ), the block version of the moment method may be useful (see Section 4.3.2).

Computing Invariant Pairs via Moment Pencils

Let Γ be a closed contour, λ 1 , . . . , λ m all the eigenvalues of P (λ) in the interior of Γ and the matrices H 0 and H 1 defined as in (4.3). For k = 0, 1, . . . , m -1 and a nonzero vector v ∈ C n , consider the vectors:

s k = 1 2πı Γ z k P (z) -1 vdz, (4.7)
The method proposed in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] for the computation of the eigenvectors of P (λ) is based on the following result.

Theorem 15. [Thm. 3.5, [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF]] Let (λ i , w i ), i = 1, . . . , m be eigenpairs for the matrix pencil H 1 -λH 0 , where the simple, distinct, nondegenerate eigenvalues λ i belong to the interior of a given closed contour Γ.

Let S = [s 0 , • • • , s m-1 ].
Then, for i = 1, . . . , m, the vector y i = Sw i is an eigenvector of P (λ) corresponding to the eigenvalue λ i .

Theorem 15 is readily applied to invariant pairs. Corollary 2. With the hypotheses of Theorem 15, S = s 0 , s 1 , . . . , s m-1 and C = H -1 0 H 1 . Then the pair (S, C) satisfies P (S, C) = 0, i.e., (S, C) is a simple invariant pair for P (λ).

Proof. Note that the pair (Y, Λ), where Λ = diag(λ 1 , . . . , λ m ) and Y = [y 1 , . . . , y m ], is clearly an invariant pair for P (λ), that is,

P (Y, Λ) = ℓ j=0 A j Y Λ j = 0. Moreover, we know that C = H -1 0 H 1 = V -1
ΛV , where V is the classical Vandermonde matrix associated with λ 1 , . . . , λ m , and that the columns of V -1 are eigenvectors of H 1 -λH 0 . So we have

0 = ℓ j=0 A j Y Λ j = ℓ j=0 A j Y Λ j V = = ℓ j=0 A j Y V V -1 Λ j V = ℓ j=0 A j SC j = P (S, C),
that is, (S, C) is also an invariant pair of P (λ).

What can we say about more general cases, where some of the hypotheses of Theorem 15 are removed? If we remove the hypothesis that the λ i 's are distinct, we can prove the following.

Theorem 16. With the hypotheses of Theorem 14, let S = s 0 , s 1 , . . . , s m-1 and C = H -1 0 H 1 . Then the pair (S, C) satisfies P (S, C) = 0, i.e., (S, C) is a simple invariant pair for P (λ).

Proof. Consider again the columns q 1 (λ), . . . , q n (λ) of the matrix F (λ) in the Smith form (1.7) and the definition of s k given in (4.7). A similar computation to (4.4) shows that

S = [s 0 , . . . , s m-1 ] = QV, where Q = [Q 0 , . . . , Q s ],
Q j = [γ 0,j q n (λ j ), γ 1,j q ′ n (λ j ), . . . , γ m j -1,j q (m j -1) n (λ j )], for j = 0, . . . , s, the γ i,j 's are complex coefficients and V is the confluent Vandermonde matrix defined above.

It is shown in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] (Lemma 2.4) that, if a complex number ζ is a root of d j (λ) for some index 1 ≤ j ≤ n, then P (ζ)q j (ζ) = 0. In our case, this implies that the vector q n (λ) is a root polynomial of P (λ) corresponding to the eigenvalue λ j , for each j = 0, . . . , s; see [START_REF] Gohberg | Matrix polynomials[END_REF], section 1.5, for the definition and properties of root polynomials. It follows that [q n (λ j ), q ′ n (λ j ), . . . , q

(m j -1)(λ j ) n
] forms a Jordan chain for the eigenvalue λ j . So we have that (Q, J) is an invariant pair for P (λ). Moreover, if C = H -1 0 H 1 as usual, we have

0 = ℓ j=0 A j QJ j = ℓ j=0 A j QJ j V = = ℓ j=0 A j QV V -1 J j V = ℓ j=0 A j SC j = P (S, C),
Therefore, (S, C) is a simple invariant pair for P (λ).

Example 16. Consider the matrix polynomial:

P (λ) = λ 2 1 0 0 1 + λ -2 0 2 -1 + 1 0 0 0 ,
which has eigenvalues λ 1 = 0, with algebraic multiplicity 1 and, λ 2 = 1, with algebraic multiplicity 3. Suppose we are interested in the eigenvalues λ 2 . Then, we can choose a contour

Γ(t) = z 0 + Re ıt , t ∈ [0, 2π], where z 0 = 1 and R = 1 2 . Choosing the vectors u = 1 -1 T and v = -1 1 T
, we find:

H 0 =    -1 -2 -5 -2 -5 -10 -5 -10 -17    , H 1 =    -2 -5 -10 -5 -10 -17 -10 -17 -26    .
Then, we have that the pair (S, C) given by Theorem 16

S = H -1 0 H 1 =    0 0 1 1 0 -3 0 1 3    and C = 0 -1 -2 1 1 3
is an invariant pair, i.e., it satisfies P (S, C) = 0. Note that the companion matrix C is associated with the monic polynomial:

p(λ) = λ 3 -3λ 2 + 3λ -1,
which has a triple root equal to 1.

The Block Moment Method

Instead of the scalar version of the moment method, we can consider a Hankel pencil constructed by block moments M k ∈ C ξ×ξ , for a suitable positive integer ξ.

Definition 16. Let k be a positive integer and U, V ∈ C n×ξ nonzero matrices with linearly independent columns. For k = 0, 1, . . ., define the block moment M k ∈ C ξ×ξ as:

M k = 1 2πı Γ z k U H P (z) -1 V dz.
Then the block Hankel matrices H ξ0 , H ξ1 ∈ C mξ× mξ are defined as:

H ξ0 =       M 0 M 1 • • • M m-1 M 1 M 2 • • • M m . . . . . . . . . M m-1 M m • • • M 2 m-2       , H ξ1 =       M 1 M 2 • • • M m M 2 M 3 • • • M m+1 . . . . . . . . . M m M m+1 • • • M 2 m-1      
Polynomial eigenvalue computation via the eigenvalues of the pencil H ξ1 -λH ξ0 is discussed in [START_REF] Asakura | A numerical method for polynomial eigenvalue problems using contour integral[END_REF] and [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF]. See also [START_REF] Leblanc | Solving acoustic nonlinear eigenvalue problems with a contour integral method[END_REF] for an application to acoustic nonlinear eigenvalue problems.

Invariant pairs can be computed from block moments by applying an approach that is similar to the one described in the previous section for the scalar version. For k = 0, 1, . . . , m -1, consider the matrices S k ∈ C n×ξ defined as:

S k = 1 2πı Γ z k P (z) -1 V dz.
Then, we have the following result.

Computation of Invariant Pairs and Matrix Solvents

Proposition 5. Let Γ be a closed contour, let the block Hankel matrix H ξ0 ∈ C mξ× mξ be nonsingular and m be the number of eigenvalues inside of Γ. If mξ = m and Y = [S 0 , . . . , S m-1 ], T = H -1 ξ0 H ξ1 , then the pair (Y, T ) satisfies P (Y, T ) = 0, i.e., (Y, T ) is a simple invariant pair for P (λ).

Example 17. Consider again the matrix polynomial of Example 15:

P (λ) = λ 2    1 0 0 2 1 0 -1 1 -2    + λ    0 0 0 -4 -2 0 2 -2 4    +    -2 1 -2 2 1 0 -1 1 -2    .
with associated Smith form:

D(λ) = diag ((d 1 (λ), d 2 (λ), d 3 (λ)) =    1 0 0 0 (λ -1) 2 0 0 0 (λ -1) 3 (λ + 1)    .
In Example 15, we found that the scalar moment method, i.e. when ξ = 1, missed the additional multiplicities associated with the polynomial d 2 (λ).

Consider now ξ = 2 as the size of the block moments M k , the contour ϕ(t) = 1 + 1 10 e ıt , containing 5 eigenvalues λ = 1, as before, and the matrices:

U =    1 0 5 -3 2 -4    , V =    1 3 0 1 -2 4    .
We find the block moments M k : Then, we have the Hankel matrix H L0 :

M 0 = -9 -12 9 12 , M 1 = -1 -22 -1 27 , M 2 = -5 -8 1 18 ,
H ξ0 =    M 0 M 1 M 2 M 1 M 2 M 3 M 2 M 3 M 4    =            -9 -12 -1 -22 -5 -8 9
           .
The matrix H ξ0 is singular. This happens because there are just 5 eigenvalues inside the contour and H ξ0 has size 6×6. Then, we have to reduce the matrices H ξ0 and H ξ1 to match the number of eigenvalues in the contour. Therefore, we get the truncated matrices: Then, we obtain:

Ĥξ0 =         -9 -12 -1 -22 -5 9 12 -1 27 1 -1 -22 -5 -8 -21 -1 27 1 18 15 -5 -8 -21 30 -49         , Ĥξ1 =         -1 -22 -5 -8 -21 -1 27
T = Ĥ-1 ξ0 Ĥξ1 =         0 0 0 -2 1 0 0 0 -1 0 1 0 0 4 -3 0 1 0 2 0 0 0 1 -2 3         .
The eigenvalues of the matrix T are 1, 1, 1, 1, 1, which are all the eigenvalues inside the contour. Moreover, computing the matrix Y = [S 0 , S 1 , S 2 ], using we get:

Ŷ =    0 1 1 2 0 0 -2 -2 0 0 0 -3 2 -7 2 -3 -4    .
Then, ( Ŷ , T ) is an invariant pair for P (λ).

Experimentally, we noted that the block method allows us to better "capture" the multiplicity structure of eigenvalues, when there are several Jordan blocks per eigenvalue. Further investigation of this approach will be the topic of future work. It should be pointed out that the results in [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF], and particularly Theorem 3.3, provide useful insight into a generalized block moment method and into the (good) behavior of the method in presence of multiple eigenvalues.

With the condition that the size of the block Hankel matrix H ξ0 ∈ C mξ× mξ is equal to the number of eigenvalues inside of Γ, i.e., if mξ = m, we get:

T = H -1 ξ0 H ξ1 =          0 0 • • • 0 -X 0 I 0 • • • 0 -X 1 0 I • • • 0 -X 2 . . . . . . . . . . . . . . . 0 0 • • • I -X m-1          , where       -X 0 -X 1 . . . -X m-1       = H -1 ξ0       M m M m+1 . . . M 2m-1       .
Consequently, since T has a block companion form, the problem of finding the eigenvalues λ 1 , . . . , λ m of is equivalent to the problem of finding the eigenvalues of the matrix polynomial:

L(λ) := λ ℓ + X ℓ-1 λ ℓ-1 + • • • + X 1 λ + X 0 = 0.
Example 18. Consider the matrix polynomial:

P (λ) = λ 4 I + λ 3    -12 0 0 0 -12 0 0 0 -12    + λ 2    41 0 0 0 41 0 0 0 41    + λ    -30 0 0 0 -30 0 0 0 -30    +    0 1 0 0 0 0    ,
which has eigenvalues: 0, 0, 0,

Suppose that in this example we only wish to compute the eigenvalues 0 and 1. We choose a contour Γ enclosing just those eigenvalues; for example, we can take the circle Γ(t) = z 0 + Re ıt , where z 0 = 1 2 , R = 1. For L = 3, we find

T = H -1 L0 H L1 =            0 0 0 1
0 0 1 0 0 1            .
The companion matrix T is associated with the monic quadratic matrix polynomial:

H(λ) = λ 2 I + λX 1 + X 0 ,
where

X 0 = -    1 30 1 30 1 30 -1 30 -1 30 -1 30 0 0 0    , X 1 = -    61 60 1 60 1 60 -1 60 59 60 -1 60 0 0 1    ,
which has as eigenvalues: 0,0,0,1,1,1 as sought.

In [START_REF] Beckermann | On the numerical condition of a generalized Hankel eigenvalue problem[END_REF], it is shown that the conditioning of the generalized Hankel eigenvalue problem (H 1 -λH 0 )y = 0 grows exponentially in terms of a quantity that depends on the largest distance between eigenvalues (for conditioning on Hankel matrices see, e.g., [START_REF] Beckermann | On the numerical condition of polynomial bases: estimates for the condition number of Vandermonde, Krylov and Hankel matrices[END_REF], [START_REF] Beckermann | The condition number of real Vandermonde, Krylov and positive definite Hankel matrices[END_REF]). Therefore, our computation of invariant pairs and solvents using the (block) moment method can be affected by ill conditioning. One possible solution to this problem consists in breaking the contour Γ, which contains all the eigenvalues, into smaller contours Γ i . The small contours should be chosen so that each of them contains clustered eigenvalues. This improves the conditioning for each subproblem, and therefore for the whole problem. This strategy is consistent with the fact that invariant pairs typically present a particular interest for clustered groups of eigenvalues. For instance, consider Figure 4.1 and suppose that the largest distance between eigenvalues λ i is large. Then, knowing the locations of the eigenvalues, we can group them into smaller non-overlapping contours Γ i as shown in Figure 4.2.

After computing the simple invariant pairs (X i , S i ) associated with each contour, using the (block) moment method, we can construct a simple invariant pair (X, S) as follows:

X = X 1 X 2 • • • , S = diag(S 1 , S 2 , • • • ).

Choosing the Contour

The choice of Γ is, of course, a crucial step when applying the contour integral formulation of the eigenvalue or invariant pair problem. If some information about the localization of the eigenvalues of P (λ) is available, one can choose the contour accordingly. Some works have addressed this localization problem; for instance: in [START_REF] Bini | Locating the eigenvalues of matrix polynomials[END_REF] and [START_REF] Melman | Generalization and variations of Pellet's theorem for matrix polynomials[END_REF] are presented generalizations of matrix version of Pellet's theorem to the case of matrix polynomials.

In [START_REF] Higham | Bounds for eigenvalues of matrix polynomials[END_REF], N. Higham and F. Tisseur found upper and lower bounds to the moduli of the eigenvalues of a matrix polynomial. On the other hand, in [START_REF] Noferini | Tropical roots as approximations to eigenvalues of matrix polynomials[END_REF], the authors used tropical roots as approximations to eigenvalues of matrix polynomials.

A related question is: how many eigenvalues of P (λ) live inside a given contour? Even an approximate estimate can be useful to choose Γ and k consistently. An answer to this problem is provided in [START_REF] Futamura | Stochastic estimation method of eigenvalue density for nonlinear eigenvalue problem on the complex plane[END_REF] and [START_REF] Futamura | Parallel stochastic estimation method of eigenvalue distribution[END_REF] and it is based on the following result.

Theorem 17. [Thm. 2, [START_REF] Futamura | Stochastic estimation method of eigenvalue density for nonlinear eigenvalue problem on the complex plane[END_REF]] Let F (z) be an n × n regular analytic matrix function and let tr(F (z)) be the matrix trace of F (z). In addition, let m be the number of eigenvalues, counting multiplicity, inside closed curves Γ on the complex plane for the nonlinear eigenvalue problem: F (λ)x = 0. Then, we have:

m = Γ tr F (z) -1 dF (z) dz dz, (4.8)
where det(F (z)) = 0.

The equation (4.8) can be approximated by an N -point quadrature rule by:

m ≈ N -1 j=0 ω j tr(F (z j ) -1 F ′ (z j )), (4.9)
where z j is a quadrature point and ω j is a weight.

In the case when we use the trapezoidal rule on a circle with center C and radius R, quadrature points and weights are defined by:

ω j = R N exp 2πi N j + 1 2 , z j = C + R exp 2πi N j + 1 2 .
To avoid the matrix inversion in (4.9), an estimate for the trace with an unbiased estimation can be used:

tr(F (z j ) -1 F ′ (z j )) ≈ 1 L L i=1 v T i F (z j ) -1 F ′ (z j )v i ,
where v i are the sample vectors, with entries equal to 1 or -1 with equal probability and L is the number of sample vectors.

The number m of eigenvalues inside the contour can then be estimated as:

m ≈ 1 L N -1 j=0 ω j L i=1 v T i F (z j ) -1 F ′ (z j )v i . (4.10)
Remark 11. The choice of Γ can be combined with shifting techniques for the eigenvalues of P (λ): see for instance [START_REF] Meini | A "shift-and-deflate" technique for quadratic matrix polynomials[END_REF].

Numerical Approximation: Trapezoid Rule for

Moments

Consider the equation (4.2). Assuming that Γ has a 2π-periodic smooth parametrization:

ϕ ∈ C 1 (R, C), ϕ(t + 2π) = ϕ(t) ∀t ∈ R.
Then, for k = 0, . . . , 2m -1, we have:

µ k = 1 2πı Γ z k f (z)dz = 1 2πı 2π 0 ϕ(t) k f (ϕ(t))ϕ ′ (t)dt. (4.11)
Taking equidistant nodes t j = 2jπ N , j = 0, . . . , N , and using the trapezoid rule (A.6), we obtain the approximation:

µ k ≈ 1 ıN N -1 j=0 ϕ(t j ) k f (ϕ(t j ))ϕ ′ (t j ). (4.12)
Note that for the special case when ϕ(t) is the circle: ϕ(t) = C + Re ıt , we obtain the formula for (4.12):

µ k ≈ R N N -1 j=0 exp 2πıj N ϕ(t j ) k f (ϕ(t j )) = (4.13) = R N N -1 j=0 exp 2πıj N C + R exp 2πıj N k f C + R exp 2πıj N .

Error Analysis of the Trapezoid Rule for Moments

In this section, we present a explicit formulation for the error of the trapezoid rule for the moments µ k , in the case when Γ is a closed contour in the complex plane. We follow a similar approach as W.-J. Beyn in [START_REF] Beyn | An integral method for solving nonlinear eigenvalue problems[END_REF].

Theorem 18. Let P (λ) be an n × n matrix polynomial and h ∈ H(A(a -, a + ), C) be holomorphic on the annulus:

A(a -, a + ) = z ∈ C : R a - < |z| < Ra + , a ± > 1,
for some R > 0.

Then, the approximation error for the trapezoidal quadrature for the moments µ k :

E N (f ) = 1 2πı |z|=R z k f (z)dz - 1 N N -1 j=0 f (Re 2ıjπ N )(Re 2ıjπ N ) k+1
satisfies for all 1 < ρ -< a -, 1 < ρ + < a + :

|E N (f )| ≤ (ρ + R) k+1 M 1 ρ N + -1 + R ρ - k+1 M 2 ρ N --1 ,
for some M 1 and M 2 such that:

u H P (ρ + Re ıt ) -1 v ≤ M 1 , u H P R ρ - e ıt -1 v ≤ M 2 .
Proof. Consider equations (4.11) and let z = Re ıt . Then we have:

µ k = 1 2πı Γ z k f (z)dz = 1 2π 2π 0 (Re ıt ) k f (Re ıt )(Re ıt )dt = 1 2π 2π 0 f (Re ıt )(Re ıt ) k+1 dt.
(4.14) Now, define the function h(z) as:

h(z) = f (z)z k+1 .
Using the Laurent expansion of h, we have:

h(z) = ∞ ℓ=-∞ h ℓ z ℓ , (4.15)
with coefficients:

h ℓ = 1 2πı |z|=R h(z)z -ℓ-1 dz = 1 2πı 2π 0 h(Re ıt )(Re ıt ) -ℓ-1 (ıRe ıt )dt = = R -ℓ 2π 2π 0 h(Re ıt )e -ıℓt dt. (4.16)
From (4.14) and (4.16), we have:

I := µ k = h 0 . (4.17)
For any positive integer N , we define the trapezoidal rule approximation to µ k in (4.14) by:

µ k ≈ 1 N N -1 j=0 f (z j )z k+1 j = 1 N N -1 j=0 h(z j ) := I N , (4.18)
where z j = Re ıt j for equidistant nodes t j = 2jπ N , j = 0, . . . , N -1. From (4.18) and (4.15), we have:

I N = 1 N N -1 j=0 h(z j ) = 1 N ∞ ℓ=-∞ h ℓ N -1 j=0 z ℓ j = = 1 N ∞ ℓ=-∞ h ℓ N -1 j=0 (Re 2ıπj
N ) ℓ = 1 N ∞ ℓ=-∞ R ℓ h ℓ N -1 j=0 e 2ıπjℓ N .
If ℓ is multiple of N , the numbers e

2ıπjℓ N are all equal to 1, and then, the second sum in previous formula is equal to N . On the other hand, if j is not multiple of N , then that sum is equal to 0.

Then, we have:

I N = ∞ ℓ=-∞ R ℓN h ℓN (4.19)
Therefore, using (4.17) and (4.19), we obtain:

I N -I = ∞ ℓ=-∞ R ℓN h ℓN -h 0 = ∞ ℓ=1 R ℓN h ℓN + R -ℓN h -ℓN .
Now, to estimate the error, we need a bound on the coefficients h jN . Using Cauchy's Theorem, we can shift the contour from |z| = R to |z| = ρ + R (i.e. taking z = ρ + Re ıt ). Then, we obtain:

|h ℓ | = 1 2πı |z|=ρ + R h(z)z -ℓ-1 dz = (ρ + R) -ℓ 2π 2π 0 h(ρ + Re ıt )e -ıℓt dt ≤ (ρ + R) -ℓ 2π 2π 0 h(ρ + Re ıt ) dt,
where:

h(ρ + Re ıt ) = f (ρ + Re ıt )(ρ + Re ıt ) k+1 ≤ (ρ + R) k+1 f (ρ + Re ıt ) = = (ρ + R) k+1 u H P (ρ + Re ıt ) -1 v ≤ (ρ + R) k+1 M 1 ,
for some M 1 such that u H P (ρ

+ Re ıt ) -1 v ≤ M 1 .
Then, we have:

|h ℓ | ≤ (ρ + R) -ℓ 2π 2π 0 h(ρ + Re ıt ) dt ≤ (ρ + R) k+1-ℓ M 1 .
Similarly, using Cauchy's Therorem, we can shift the contour from |z| = R to |z| = R ρ - (i.e. taking z = R ρe ıt ). Then we find that:

|h -ℓ | ≤ R ρ - ℓ 1 2π 2π 0 h R ρ - e ıt dt ≤ R ρ - k+1+ℓ M 2 ,
for some M 2 such that u

H P R ρ -e ıt -1 v ≤ M 2 .
Using this, we find that:

|I N -I| = ∞ ℓ=1 R ℓN h ℓN + R -ℓN h -ℓN = ∞ ℓ=1 R ℓN h ℓN + R -ℓN h -ℓN ≤ ∞ ℓ=1 R ℓN |h ℓN | + ∞ ℓ=1 R -ℓN |h -ℓN | ≤ ∞ ℓ=1 R ℓN (ρ + R) k+1-ℓN M 1 + ∞ ℓ=1 R -ℓN R ρ - k+1+ℓN M 2 = (ρ + R) k+1 M 1 ∞ ℓ=1 ρ -ℓN + + R ρ - k+1 M 2 ∞ ℓ=1 ρ -ℓN - = (ρ + R) k+1 M 1 ρ -N + 1 -ρ -N + + R ρ - k+1 M 2 ρ -N - 1 -ρ -N - = (ρ + R) k+1 M 1 1 ρ N + -1 + R ρ - k+1 M 2 1 ρ N --1 .
Then, the error for the trapezoid sum for the moment µ k is bounded by:

|I N -I| ≤ (ρ + R) k+1 M 1 ρ N + -1 + R ρ - k+1 M 2 ρ N --1 .
Example 19. Consider the matrix polynomial (2.12) and the annulus:

A(a -, a + ) = z ∈ C : R a - < |z| < a + R , a ± > 1,
for some R > 0. Suppose that 1 < ρ -< a -and 1 < ρ + < a + . Now, consider in our example the circle of center C = 0 and radius R = . Then, we obtain the annulus:

A = z ∈ C : 1 10 < |z| < 4 , a ± > 1.

Note that

• As we have that ρ -< a -, then we can take ρ -= 31 2 = 15.5 (⇒ R ρ -= 1 5 ).

• Note that: 3.1 < Rρ + < Ra + = 4. We will take several values of Rρ + to study the behavior of the error.

Matrix Solvents

Once an invariant pair has been numerically computed or approximated, it can be refined using an iterative method such as Newton: this is, for instance, the strategy proposed in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF]. Here, we experiment with some modifications to the classical Newton's method applied to the equation P (X, S) = 0.

Newton's Method

Newton's method, also called the Newton-Raphson method, is a root-finding algorithm that uses the first few terms of the Taylor series of a function f (x) in the vicinity of a suspected root (see [START_REF] Newton | Methodus fluxionum et serierum infinitarum[END_REF], [START_REF] Raphson | Analysis aequationum universalis[END_REF]). Consider the Taylor series of f (x) expanded about the point x = x 0 + ǫ:

f (x 0 + ǫ) = f (x 0 ) + f ′ (x 0 )ǫ + 1 2 f ′′ (x 0 )ǫ 2 + • • • .
Keeping terms only to first order, we have:

f (x 0 + ǫ) ≈ f (x 0 ) + f ′ (x 0 )ǫ.
Setting f (x 0 + ǫ) = 0 and solving last equation for ǫ = ǫ 0 , gives:

ǫ 0 = - f (x 0 ) f ′ (x 0 ) ,
which is the first-order adjustment to the root's position. Taking x 1 = x 0 + ǫ 0 and calculating a new ǫ 1 , and so on, the process can be repeated until it converges to a fixed point (which is precisely a root) using:

ǫ n = - f (x n ) f ′ (x n ) .
Then, given a good initial choice of the root, the algorithm is iteratively given by:

x k+1 = x k - f (x k ) f ′ (x k ) , for n = 1, 2, . . . .
In the case of invariant pairs, Newton's method defines (∆X, ∆S) by P (X, S) + DP (X,S) (∆X, ∆S) = 0.

(5.1)

And for the case of matrix solvents, Newton's method defines ∆S by P (S) + DP S (∆S) = 0.

(5.2)

As explained in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], the definition of a simple invariant pair (X, S), i.e., a pair such that P (X, S) = 0, is not sufficient to characterize (X, S). Then, we must add the condition

W H V m (X, S) = I k ,
where m ≤ ℓ is not smaller than the minimality index of (X, S) and the columns of

W = [W H m-1 , . . . , W H 0 ]
H form an orthonormal basis of span(V m (X, S)). In other words, the idea of the Newton's method is the following: given an initial approximation (X 0 , S 0 ) to a simple invariant pair (X, S) ∈ C n×k × C k×k our goal is to compute a correction (∆X k , ∆S k ), which brings(X 0 , S 0 ) closer to (X, S). This is, the Newton's method applied to the equations:

P (X, S) = 0, V(X, S) = 0, where V(X, S) = W H V m (X, S)-I k , for some normalization matrix W = [W H m-1 , . . . , W H 0 ] H ∈ C k×mn , takes the form: (X k+1 , S k+1 ) = (X k , S k ) -L -1 k (P (X k , S k ), V(X k , S k )) ,
where L k is the Jacobian at the current iterate (X k , S k ):

L k (∆X, ∆S) = P (∆X, S k ) + ℓ j=1 A j X k DS j k (∆S), m-1 j=1 W H j (∆XS j k + XDS j k (∆S)) .
The methods for solving the correction equations will be discussed in Section 5.5. For now, we will describe the global algorithm and its variants.

Algorithm: Newton's Method

Given an initial approximation (X 0 , S 0 ) to the solution of (2.5) (resp. S 0 to the solution of (3.2)), a tolerance ǫ and a contour Γ ∈ C (with the spectrum of S 0 in its interior), we have:

Algorithm 1. STEP 1: Set k = 0 STEP 2: If error k < ǫ: STOP STEP 3: Solve for (∆X k , ∆S k ) (resp. for ∆S k ) the equation: DP (X,S) (∆X k , ∆S k ) = -(P (X k , S k ), 0) (resp. DP S (∆S k ) = -(P (S k ), 0)) STEP 4: Update: • X k+1 = X k + ∆X k , S k+1 = S k + ∆S k (resp. S k+1 = S k + ∆S k ). • k = k + 1 • go to STEP 2.
Remark 12. In the case of invariant pairs, we define error k as: error

k = P (X k ,S k ) F X k F
and in the case of matrix solvents: error k = P (S k ) F S k F .

Incorporating Line Search into Newton's Method

Line searches are relatively inexpensive and improve the global convergence properties of Newton's method. Each iteration of a line search method computes a search direction d k and then decides how far to move along that direction. The iteration is given by:

x k+1 = x k + t k d k ,
where the positive scalar t k is the step length. The success of a line search method depends on effective choices of both the direction d k and the step length t k (see [START_REF] Wright | Numerical optimization[END_REF]). A value t k = 1 gives the original Newton iteration.

In this section, we show how to incorporate exact line searches into Newton's method to approximate invariant pairs (X, S) and solvents for the general matrix solvent problem P (S) = 0. For that, we use the contour integral formulations (2.5) and (3.2).

As pointed out before, the use of Newton's method incorporating line searches to find solvents is not new. For instance, in [START_REF] Higham | Solving a quadratic matrix equation by Newton's method with exact line searches[END_REF] and [START_REF] Long | Improved Newton's method with exact line searches to solve quadratic matrix equation[END_REF] this approach is used to approximate solvents for the quadratic matrix equation.

Case of Invariant Pairs

In the specific problem (2.5), the direction d k is given by the solution (∆X k , ∆S k ) of the correction equation (5.5) (see Section 5.5). The step length t k on each iteration is given by the solution of the minimization problem: p(t) = P (X + t∆X, S + t∆S) 2 F .

where:

A = 1 2πı Γ P (λ) ∆X + X(λI -S) -1 ∆S (λI -S) -1 ∆S(λI -S) -1 dλ, B = 1 2πı Γ P (λ)∆X(λI -S) -1 ∆S(λI -S) -1 ∆S(λI -S) -1 dλ, α = P (X, S) 2 F , θ = A 2 F , ϕ = B 2 F , η = trace(A * B + B * A
), β = trace(P (X, S) * A + A * P (X, S)), γ = trace(P (X, S) * B + B * P (X, S)).

Case of Matrix Solvents

For the problem (3.2) we have the minimization problem:

p(t) = P (S + t∆S) 2 F .
Using the formula for the derivative of the matrix inverse:

dA -1 dt = -A -1 dA dt A -1 ,
we write the formula for the total derivative of P at S in direction ∆S as follows:

DP S (∆S) = 1 2πı Γ P (λ)(λI -S) -1 ∆S(λI -S) -1 dλ.

(5.4)

Algorithm: Newton's Method with Line Search

Given an initial approximation (X 0 , S 0 ) to the solution of (2.5) (resp. S 0 to the solution of (3.2)), a tolerance ǫ and a contour Γ ∈ C (with the spectrum of S 0 in its interior), we have: P (S + t∆S) 2 F ) STEP 5: Update:

• X k+1 = X k + t∆X k , S k+1 = S k + t∆S k (resp. S k+1 = S k + t∆S k ). • k = k + 1
• go to STEP 2.

Šamanskii's Technique

Another variation on Newton's method is Šamanskii's technique, which accelerates the quadratic convergence to cubic. Therefore, we must ensure we have quadratic convergence in order to use this technique (see [START_REF] Šamanskii | On a modification of the Newton method[END_REF]).

Algorithm: Newton's Method with Line Search and Samanskii Technique

Given an initial approximation (X 0 , S 0 ) to the solution of (2.5) (resp. S 0 to the solution of (3.2)), tolerances ǫ and ǫ 0 and a contour Γ ∈ C (with the spectrum of S 0 in its interior), we have:

P (S + t∆S) 2 F ) STEP 6: Update • X k+1 = X k + t∆X k , S k+1 = S k + t∆S k (resp. S k+1 = S k + t∆S k ). • k = k + 1 and go to STEP 2. STEP 7: Update X k,1 = X k + ∆X k and S k,1 = S k + ∆S k (resp. S k,1 = S k + ∆S k ). STEP 8: Solve for (∆X k,1 , ∆S k,1 ) (resp. ∆S k,1
) the equation:

DP (X,S) (∆X k,1 , ∆S k,1 ) = -(P (X k,1 , S k,1 ), 0) (resp. DP S ∆S k,1 = -(P (S k,1 ), 0)) STEP 9: Update • X k+1 = X k,1 + ∆X k,1 , S k+1 = S k,1 + ∆S k,1 (resp. S k+1 = S k,1 + ∆S k,1
).

• k = k + 1 and go to STEP 2.

Solution of the Correction Equation

In this section, we use Newton's Method and its variations to approximate a solution for equation (2.5) (resp. for equation(3.2)). We must find the correction (∆X, ∆S) (resp. ∆S) for the equations:

DP (X,S) (∆X k , ∆S k ) = -(P (X k , S k ), 0) (5.5)
and respectively for:

DP S (∆S) = -(P (S k ), 0).
In next sections, we will show two ways to find this correction. The first one uses a Kronecker product approach and the second uses forward substitution. We will focus in the case of the correction equation for invariant pairs; the corresponding equation for matrix solvents can be used in a similar way.

Using the Kronecker Product

Consider the equation (5.5). Vectorizing its left part and using (A.1), we obtain:

vec (DP (∆X, ∆S)) = = vec 1 2πı Γ P (λ) ∆X + X(λI -S) -1 ∆S (λI -S) -1 dλ = = 1 2πı Γ vec P (λ) ∆X + X(λI -S) -1 ∆S (λI -S) -1 dλ = = 1 2πı Γ vec P (λ)∆X(λI -S) -1 dλ+ + 1 2πı Γ vec P (λ)X(λI -S) -1 ∆S(λI -S) -1 dλ = = 1 2πı Γ (λI -S) -T ⊗ P (λ) dλ vec (∆X)+ + 1 2πı Γ (λI -S) -T ⊗ P (λ)X(λI -S) -1 dλ vec (∆S).
Then, we can rewrite the equation DP (∆X, ∆S) = (P (X p , S p ), 0) as:

K 11 K 12 K 21 K 22 vec (∆X) vec (∆S) = vec (P (X p , S p )) 0 , (5.6)
where

K 11 = 1 2πı Γ (λI -S) -T ⊗ P (λ) dλ, K 12 = 1 2πı Γ (λI -S) -T ⊗ P (λ)X(λI -S) -1 dλ, K 21 = m-1 j=0 ((S j ) T ⊗ W H j ), K 22 = m-1 j=1 (I k ⊗ W H j X)K S j ,
where K S j denotes the Kronecker product formulation of the Fréchet derivative, i.e.:

K S j = j-1 i=0
((S j-i-1 ) T ⊗ S i .

Using Forward Substitution

Using the Schur decomposition of S p and an appropriate unitary transformation of the pair (X p , S p ), we can assume that S p is in upper triangular form. Using this special structure of S p , we can easily find the columns of ∆X and ∆S successively using a forward substitution process.

Note that using the upper triangular structure of S p in the equation:

DP (∆X, ∆S) = (P (X p , S p ), 0) (5.7)
the computation of the first columns ∆x 1 and ∆s 1 of ∆X and ∆S, respectively, is simple. Assuming that S is an upper triangular matrix, then (λI -S) -1 is also upper triangular. Its diagonal entries are (λI -s 11 ) -1 , (λI -s 22 ) -1 , . . . , (λI -s nn ) -1 , where s 11 , s 22 , . . . , s nn are the diagonal entries of S. In order to compute ∆x 1 and ∆s 1 , we can write:

DP (X, S) = 1 2πı Γ P (λ) ∆x 1 + X(λI -S) -1 ∆s 1 (λ -s 11 ) -1 dλ = = 1 2πı Γ P (λ)∆x 1 (λ -s 11 ) -1 dλ + 1 2πı Γ P (λ)X(λI -S) -1 ∆s 1 (λ -s 11 ) -1 dλ.
Then, we have:

B 11 B 12 B 21 B 22 ∆x 1 ∆s 1 = r 1 0 , (5.8)
where r 1 denotes the first column of Res = P (X p , S p ) and:

B 11 = 1 2πı Γ P (λ)(λ -s 11 ) -1 dλ, B 12 = 1 2πı Γ P (λ)X(λI -S) -1 (λ -s 11 ) -1 dλ, B 21 = m-1 j=0 s j 11 W H j , B 22 = m-1 j=1 W H j X DS j 11 ,
Multiplying last equation by 0 I k-1 to consider just the last k -1 columns, we have: obtaining for the pair (∆X 2 , ∆S 2 ) ∈ C n×(k-1) × C n×(k-1) the linear matrix equation: Taking r 2 and q 2 as the first columns of Res 2 and Ort 2 , respectively, we can find the second columns ∆x 2 and ∆s 2 of ∆X and ∆S, respectively, solving the linear system: W H j X DS j 22 .

1 2πı Γ P (λ) ∆x 1 (λ -s 11 ) -1 ∆x 1 (λ -s 11 ) -1 s 12 (λI -S 22 ) -1 0 I k-1 dλ+ + 1 2πı Γ P (λ) 0 ∆X 2 (λI -S 22 ) -1 0 I k-1 dλ+ + 1 2πı Γ P (λ)X(λI -S) -1 ∆s 1 (λ -s 11 ) -1 ∆s 1 (λ -s 11 ) -1 s 12 (λI -S 22 ) -1 0 I k-1 dλ+ + 1 2πı Γ P (λ)X(λI -S) -1 0 ∆S 2 (λI -S 22 ) -1 0 I k-1 dλ = 1 2πı Γ P (λ)∆x 1 (λ -s 11 ) -
B

Numerical Results

In this section we compare two methods to refine approximate invariant pairs (X, S) ∈ C n×k × C k×k : Newton's method (N.M.) presented in [START_REF] Betcke | Perturbation, extraction and refinement of invariant pairs for matrix polynomials[END_REF], Newton's method with line search (N.M.L.S.), explained in Section 5.3 and Newton's method with line search and Samanskii technique (N.M.L.S.S.), explained in Section 5.4.

We have implemented Newton's method with the above variants in MATLAB and applied it to several problems taken from the NLEVP collection (see [START_REF] Betcke | NLEVP: A collection of nonlinear eigenvalue problems[END_REF]). The results are shown in Table 5.1. For each problem, an initial invariant pair (X 0 , S 0 ) has first been approximated using the (block) moment method of Section 4.3.1 and approximating the moments µ i in (4.2) via the trapezoid rule discussed in Section 4.5, with N = 20 integration nodes. Moreover, Γ is chosen for each problem as the contour enclosing the k eigenvalues with largest condition number (computed using the MATLAB function polyeig).

Table 5.1 shows that line search is generally effective in reducing the number of iterations and the overall computation time.

Figure 5.1 shows the convergence of the Newton's method with line search, for the Dirac problem presented in Table 5.1. Here we use as contour the circle of center C = -0.1 and radius R = 1.14, which contains the 6 eigenvalues with largest condition number.

Functions Implemented in MATLAB and Maple

In this section, we describe the Maple and MATLAB implementations used in this thesis. These implementations are available online at the URL http://www.unilim.fr/pages perso/esteban.segura/software.html

MATLAB

Let P (λ) be an n × n matrix polynomial of degree ℓ and Γ a contour (we limit our implementations to the case when Γ is a circle).

• conditionNumInvPair.m: This function approximates the condition number for an invariant pair (X, S) using equation (2.10). For this computation we use the values α i = A i , for i = 0, . . . , ℓ. The input of this function is: In this work, we have explored several aspects related to invariant pairs and solvents of matrix polynomials.

After recalling some basics about matrix polynomials and their applications, and the relation with systems of ordinary differential equations, we presented a generalization of a shifting method for matrix polynomials, which we believe can be useful, for instance, when applying methods that require a particular eigenvalue distribution.

We presented some key points in the general theory and applications of the invariant pair and matrix solvent problems. Additionally, we computed new formulations for the condition number and the backward error of invariant pairs and matrix solvents. In the case of solvents, this computation generalized the existent previous work on the quadratic matrix equation. Furthermore, we explored the relationship between solvents of matrix polynomials in general and in triangularized form. This could be useful when computing matrix solvents.

In order to compute invariant pairs and matrix solvents, we studied a moment Hankel pencil method in its scalar and block versions: the method lends itself to a numeric or a hybrid symbolic-numeric application. We showed that the scalar method cannot capture some eigenvalue multiplicity structures, for which the block version is needed. Moreover, we analyzed the error for the trapezoid rule applied to the approximation of the moments via contour integrals.

When computing invariant pairs and matrix solvents, the proposed moment Hankel pencil methods, and some others direct approaches as well, may need to be refined numerically using an iterative method: here, we studied and compared two variants of Newton's method, namely line-search and Šamanskii. We tested the effectiveness and robustness of these approaches on several problems, many of them taken from the NLEVP collection. Some points in this analysis may deserve to be explored in further detail. For instance, we would like to achieve a deeper understanding of the behavior of the block moment method presented in Chapter 4: at this time, we do not have results that tell us in advance what should be the size of the blocks, depending on the structure of multiplicity of the eigenvalues.

On the other hand, we would like to present convergence theorems for Newton's method and its variants analyzed in Chapter 5. Formal convergence results would help us decide, for instance, whether the results given by direct methods are good starting points for the subsequent refinement.

Another possible generalization involves the use of numerical quadrature to approximate the contour integrals. Here we limited ourselves for simplicity to the case when the contour is a circle, but it would be interesting to extend our implementations to different types of contour (e.g. ellipses, lines, etc). This might be more advantageous for some problems with particular eigenvalue distributions. Résumé : Cette thèse porte sur certains aspects symboliques-numériques du problème des paires invariantes pour les polynômes de matrices. Les paires invariantes géneralisent la définition de valeur propre / vecteur propre et correspondent à la notion de sousespaces invariants pour le cas nonlinéaire. Elles trouvent leurs applications dans le calcul numérique de plusieurs valeurs propres d'un polynôme de matrices ; elles présentent aussi un intérêt dans le contexte des systèmes différentiels.

En utilisant une approche basée sur les intégrales de contour, nous déterminons des expressions du nombre de conditionnement et de l'erreur rétrograde pour le problème du calcul des paires invariantes. Ensuite, nous adaptons la méthode des moments de Sakurai-Sugiura au calcul des paires invariantes et nous étudions le comportement de la version scalaire et par blocs de la méthode en présence de valeurs propres multiples. Le résultats obtenus à l'aide des approches directes peuvent éventuellement être améliorés numériquement grâce à une méthode itérative : nous proposons ici une comparaison de deux variantes de la méthode de Newton appliquée aux paires invariantes.

Le problème des solvants de matrices est très proche de celui des paires invariants. Le résultats présentés ci-dessus sont donc appliqués au cas des solvants pour obtenir des expressions du nombre de conditionnement et de l'erreur, et un algorithme de calcul basé sur la méthode des moments. De plus, nous étudions le lien entre le problème des solvants et la transformation des polynômes de matrices en forme triangulaire.

Mots clés : Polynômes de matrices, paires invariantes, solvants, intégrale de contour, nombre de conditionnement, erreur rétrograde, faisceaux de matrices de Hankel, méthode de Newton.

Abstract:

In this thesis, we study some symbolic-numeric aspects of the invariant pair problem for matrix polynomials. Invariant pairs extend the notion of eigenvalueeigenvector pairs, providing a counterpart of invariant subspaces for the nonlinear case. They have applications in the numeric computation of several eigenvalues of a matrix polynomial; they also present an interest in the context of differential systems.

Here, a contour integral formulation is applied to compute condition numbers and backward errors for invariant pairs. We then adapt the Sakurai-Sugiura moment method to the computation of invariant pairs, including some classes of problems that have multiple eigenvalues, and we analyze the behavior of the scalar and block versions of the method in presence of different multiplicity patterns. Results obtained via direct approaches may need to be refined numerically using an iterative method: here we study and compare two variants of Newton's method applied to the invariant pair problem.

The matrix solvent problem is closely related to invariant pairs. Therefore, we specialize our results on invariant pairs to the case of matrix solvents, thus obtaining formulations for the condition number and backward errors, and a moment-based computational approach. Furthermore, we investigate the relation between the matrix solvent problem and the triangularization of matrix polynomials.

Definition 4 .

 4 The polynomial eigenvalue problem (PEP) consists in determining right eigenvalue-eigenvector pairs (λ, x) ∈ C × C n , with x = 0, such that P (λ)x = 0, or left eigenvalue-eigenvector pairs (λ, y) ∈ C × C n , with y = 0, such that y * P (λ) = 0. Remark 1.

Figure 2 . 1 :

 21 Figure 2.1: Location of the eigenvalues of the power plant problem

Figure 2 . 2 :

 22 Figure 2.2: Eigenvalues (blue crosses) inside contour Γ = γ + ρe iθ (γ = 80 + 10i, ρ = 170)

Theorem 8 .

 8 Let P (λ) be a n × n matrix polynomial and consider an invariant pair (Y, T ) ∈ C n×k × C k×k of P (λ). If the matrix Y has size n × n, i.e. k = n, and is invertible, then S = Y T Y -1 satisfies equation (3.1), i.e., S is a matrix solvent of P (λ).

Remark 9 .

 9 The computation of the matrix M in the Examples 12 and 13 was done following the proof of Theorem 5.2 in[START_REF] Taslaman | Triangularizing matrix polynomials[END_REF].

x 11 = 3 or x 11 = 4 , x 22 = 3 and x 33 = 4 .

 114334 Then we have two cases: I. If x 11 = 3, x 22 = 3 and x 33 = 4:

Figure 4 . 1 :

 41 Figure 4.1: Eigenvalues inside a circle Figure 4.2: Grouping clustered eigenvalues

31 10 = 3 . 1 ,

 31 i.e. z = 31 10 e ıπ . Let a + = 40 31 ≈ 1.2903, then a + R = 4; and let a -= 31, then R a -= 1 10

Figure 4 . 5 : 5 Chapter 5 :

 4555 Figure 4.5: log|E N (f )| for µ 4 and µ 5

Algorithm 2 . 1 : 2 F

 212 STEP Set k = 0 STEP 2: If error k < ǫ: STOP STEP 3: Solve for (∆X k , ∆S k ) (resp. for ∆S k ) the equation:DP (X,S) (∆X k , ∆S k ) = -(P (X k , S k ), 0) (resp. DP S (∆S) = -(P (S k ), 0))STEP 4: Find by exact line searches a value of t that minimizes the function: min t∈[0,2] P (X + t∆X, S + t∆S)

Algorithm 3 . 1 :STEP 2 :STEP 4 :STEP 5 : 2 F

 312452 STEP Set k = 0 If error k < ǫ: STOP STEP 3: Solve for (∆X k , ∆S k ) (resp. for ∆S k ) the equation:DP (X,S) (∆X k , ∆S k ) = -(P (X k , S k ),0) (resp. DP S (∆S) = -(P (S k ), 0)) If error k < ǫ 0 : go to STEP 7. Find by exact line searches a value of t that minimizes the function: min t∈[0,2] P (X + t∆X, S + t∆S) (resp. min t∈[0,2]

∆X 2 S j 22 + 1 = 2 ,j 12 +

 2221212 ) ∆X 2 + X(λI -S) -1 ∆S 2 (λI -S 22 ) -1 dλ = Res 2 XDS j ([0, ∆S 2 ]) 0 I k-Ortwhere the updated right hand sides are:Res 2 := Res 2 -1 2πı Γ P (λ) ∆x 1 + X(λI -S) -1 ∆s 1 (λ -s 11 ) -1 s 12 (λI -S 22 ) -1 dλ Ort 2 := -XDS j ([∆s 1 , 0]) 0 I k-1 .

Table 5 . 1 :

 51 Comparison of results for classical Newton, Newton with line search and Newton with line search and Šamanskii's technique.

Figure 5 . 1 :

 51 Figure 5.1: Convergence of Dirac problem using Newton's method with line search.This is a log-10 plot of the relative residual P (X,S) F X F

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: Eigenvalues inside the circle Figure 5.3: The eigenvalues that should be excluded from the computations are contained in the smaller circles

Keywords:

  Matrix polynomials, invariant pairs, solvents, contour integral, condition number, backward error, Hankel pencils, Newton's method. LABO -UMR CNRS n • 7252 123, avenue Albert Thomas -87060 LIMOGES

  2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Previous Work on Invariant Pairs . . . . . . . . . . . . . . . . . . . . . . . 2.3 Definition and Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Particular Case: Jordan Pairs . . . . . . . . . . . . . . . . . . . . . 2.4 Formulation of the Invariant Pair Problem Using the Contour Integral . . . 2.5 Linearized Matrix Equation and Fréchet Derivative . . . . . . . . . . . . . 2.6 Condition Number and Backward Error for the Invariant Pair Problem . . Condition Number and Backward Error for the Matrix Solvent Problem . . 3.2.1 Condition Number of P (S) . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Backward Error of P (S) . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Computation of Solvents . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 A Computational Approach . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Matrix p-th Root . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Solvents and Triangularized Matrix Polynomials . . . . . . . . . . . . . . . 3.4.1 Triangularizing Matrix Polynomials . . . . . . . . . . . . . . . . . . 3.4.2 Procedure to Triangularize a Quadratic Matrix Polynomial . . . . . 3.4.3 A Problem with an Infinite Number of Solvents . . . . . . . . . . .

	Notation
	N	Set of nonnegative integers
	N *	Set of positive integers
	Z	Ring of integers
	R, R n , R n×n	Set of real numbers, n-dimensional (column) real vectors
		and n × n real matrices
	C, C n , C n×n	Set of complex numbers, n-dimensional (column) com-
		plex vectors and n × n complex matrices
	K[λ]	Ring of polynomials in λ over a ring K
	A m×n	Additive group of m × n matrices with entries in a ring
		A
	A n×n	Ring of n × n matrices with entries in a ring A
	A n resp. A 1×n	Additive group of n-dimensional (column) vectors,
		resp. row vectors, with entries in a ring A
	deg(p)	Degree of the polynomial p
	0 n	Square zero matrix of size n
	I n	Identity matrix of size n
	A T

2.6.1 Condition Number for P (X, S) . . . . . . . . . . . . . . . . . . . . 2.6.2 Backward Error for P (X, S) . . . . . . . . . . . . . . . . . . . . . . Chapter 3 : Matrix Solvents . . . . . . . . . . . . . . . . . . 3.1 Matrix Solvents: Definition and Theory . . . . . . . . . . . . . . . . . . . . 3.1.1 Existence of Solvents . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Chapter 4 : Moments, Hankel Pencils and Computation of Invariant Pairs . . . . . . . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  If we compute the condition number for the eigenvalues of A and of D -1 AD using MATLAB command condeig, we obtain:

	2.170086486635809 2.170086486626034
	-1.481194304285063 -1.481194304092016
	0.311107817428706 0.311107817465982
	condeig(A) condeig(D -1 AD)
	335275.9044	1.1585
	237384.3761	1.0392
	98402.2850	1.1197

Table 2 .

 2 1: Condition numbers for matrix polynomial (2.13) using different τ

	let us compute

  , an n × n matrix H is Hankel if and only if it can be written as H = T E n , where T is a suitable Toeplitz matrix. Conversely, an n × n matrix T is Toeplitz if and only if T = H E n , where T is a suitable Hankel matrix.

Table 4 .

 4 the number of nodes N i such that the trapezoid rule error is less than 10 -6 . Then, we obtain the results presented in Tables4.1, 4.2 and 4.3, where we fix M 2 and vary M 1 to study the behavior of the problem. Moreover, in Figures 4.3, 4.4 and 4.5 we show a log-10 plot of |E N (f )| for the moments µ k , for k = 0, . . . , 5. 1: Trapezoid Error Example: µ k for k = 0, 1.

	Consider u = , and compute Rρ + -4 -7 and v = -10 11 M 1 M 2 |E N (f )| k=0 N 0 3.2 3195 71.1 10224 1.03226 N -1 + 14.22 ( 31 2 ) N -1 726 3.3 2297 71.1 7580.1 1.06452 N -1 + 14.22 ( 31 2 ) N -1 364 3.4 1905 71.1 6477 1.09677 N -1 + 14.22 ( 31 2 ) N -1 245 3.5 1739 71.1 6086.5 1.12903 N -1 + 14.22 ( 31 2 ) N -1 186 3.6 1728 71.1 6220.8 1.1613 N -1 + 14.22 ( 31 2 ) N -1 151 22394.88002 |E N (f )| k=1 32716.8 1.03226 N -1 + 2.844 ( 31 N -1 2 ) 25014.33 1.06452 N -1 + 2.844 ( 31 N -1 2 ) 22021.8 1.09677 N -1 + 2.844 ( 31 N -1 2 ) 21302.75 1.12903 N -1 + 2.844 ( 31 N -1 2 ) 1.1613 N -1 + 2.844 ( 31 N -1 2 ) 3.7 1889 71.1 6989.3 1.19355 N -1 + 14.22 ( 31 2 ) N -1 129 25860.41 1.19355 N -1 + 2.844 ( 31 N -1 2 ) 3.8 2378 71.1 9036.4 1.2258 N -1 + 14.22 ( 31 2 ) N -1 113 34338.32002 1.2258 N -1 + 2.844 ( 31 N -1 2 ) 3.9 4062 71.1 15841.8 1.2581 N -1 + 14.22 ( 31 2 ) N -1 103 61783.02 1.2581 N -1 + 2.844 ( 31 N -1 2 ) 3.98 18093 71.1 72010.14 1.2839 N -1 + 14.22 ( 31 2 ) N -1 101 286600.3574 1.2839 N -1 + 2.844 ( 31 2 ) N -1	N 1 763 383 258 196 160 136 120 109 106

N -1

Table 4 .

 4 2: Trapezoid Error Example: µ k for k = 2, 3.

	Figure 4.4: log|E N (f )| for µ 2 and µ 3

N -1

Table 4 . 3

 43

: Trapezoid Error Example: µ k for k = 4, 5.

  1 s 12 (λI -S 22 ) -1 dλ + 1 2πı Γ P (λ)∆X 2 (λI -S 22 ) -1 dλ+ + 1 2πı Γ P (λ)X(λI -S) -1 ∆s 1 (λ -s 11 ) -1 s 12 (λI -S 22 ) -1 dλ+

	+	1 2πı

Γ P (λ)X(λI -S) -1 ∆S 2 (λI -S 22 ) -1 dλ,

  11 B 12 B 21 B 22

	where:
	B 11 =	1 2πı Γ	P (λ)(λ -s 22 ) -1 dλ,
	B12 =	1 2πı Γ	P (λ)X(λI -S) -1 (λ -s 22 ) -1 dλ,
		m-1
	B 21 =	s j 22 W H j ,
		j=0
		m-1
	B 22 =
		j=1
			∆x 2 ∆s 2	=	r 2 q 2	,

  N.M.N.M.L.S. N.M.L.S.S.

	Problem	Deg P Size X	Ite Time	Ite Time	Ite Time
	bicycle	2	2 × 2	23	0.082	16	0.112		0.05
	butterfly	4	64 × 5	67	3.719	22	1.567		1.57
	cd player	2	60 × 6	500	N.C.	19	1.021		1.15
	closed loop	2	2 × 2	8	0.016	7	0.02	6	0.02
	damped beam	2	200 × 6	28	6.109	4	0.6037	3	0.80
	dirac	2	80 × 6	500	N.C.	41	1.965		2.23
	hospital	2	24 × 24	53	5.65	51	6.21		6.20
	metal strip	2	9 × 9	500	N.C.	28	0.589		0.48
	mobile manipulator	2	5 × 2	8	0.014	7	0.030	6	0.03
	pdde stability	2	225 × 6	29	9.644	16	5.622		8.17
	planar waveguide	4	129 × 6	72	11.148	19	3.682		5.34
	plasma drift	3	128 × 6	69	13.059	26	5.596		6.33
	power plant	2	8 × 8	15	0.34	13	0.39		0.50
	railtrack	2	1005 × 3	32 199.365 28 209.471 216.53

Acknowledgements

First of all, I would like to thank my supervisors, Paola Boito and Moulay Barkatou, for their friendship, support, guidance and patience throughout my years as a PhD student at the University of Limoges. Their outstanding efforts have made this thesis We would also like to work more on the design and development of effective symbolicnumeric eigenvalue algorithms based on the moment method, along with extensive numerical tests.

Chapter 3 : Matrix Solvents

Computation of Invariant Pairs and Matrix Solvents

Using the formula for the total derivative of P at (X, S) in direction (∆X, ∆S): DP (X,S) (∆X, ∆S) = 1 2πı Γ P (λ) ∆X + X(λI -S) -1 ∆S (λI -S) -1 dλ, (5.3) we have, at second order in ∆X and ∆S :

P (X + t∆X, S + t∆S) = 1 2πı Γ P (λ)(X + t∆X)(λI -S -t∆S) -1 dλ = = 1 2πı Γ P (λ)(X + t∆X) (λI -S) -1 + (λI -S) -1 t∆S(λI -S) -1 +(λI -S) -1 t∆S(λI -S) -1 t∆S(λI -S)

+t 2 (λI -S) -1 ∆S(λI -S) -1 ∆S(λI -S) -1 dλ = = 1 2πı Γ P (λ)X(λI -S) -1 dλ

Recalling that Newton's method defines (∆X, ∆S) by (5.1), we have:

Thus

Then, we have:

Recalling that Newton's method defines ∆S by (5.2), we have:

then, we have:

Thus

where:

where [DS j ] 11 denotes the Fréchet derivative of the first column of S j with respect to the first column of S:

For the second columns of ∆X and ∆S, we can find an equation of the form of (5.8), but first we must update the right side of the equation (5.7). We will now describe how to do this update. First, consider the partitions:

S = s 11 s 12 0 S 22 and S j s j 11 [S j ] 12 0 S j

22

.

Note also that using the formula for the inverse of a block upper triangular matrix [START_REF] Meyer | Generalized inverses of block triangular matrices[END_REF], we have:

Inserting this in equation ( 5.3), we have:

Chapter 6 :

Conclusions and Future Work

Computation of Invariant Pairs and Matrix Solvents

Chapter A :

Appendix

Computation of Invariant Pairs and Matrix Solvents

A.1 Kronecker Product

If A is an m × n matrix and B is a p × q matrix, then the Kronecker product A ⊗ B is the mp × nq block matrix (see [START_REF] Golub | Matrix computations[END_REF], pp. 27):

A.1.1 Kronecker Product Properties

, for matrices A and B nonsingular.

A.2 Vectorization

The vectorization of an m × n matrix A, denoted by vec (A), is the mn × 1 column vector obtained by stacking the columns of the matrix A on top of one another (see [START_REF] Golub | Matrix computations[END_REF], pp.

A.2.1 Compatibility with Kronecker Products

The vectorization is frequently used together with the Kronecker product to express matrix multiplication as a linear transformation on matrices. In particular:

A.3 Vector and Matrix Norm

A.3.1 Vector Norm

A general vector norm on K n is a function f : K n → R that satisfies the following properties (see [START_REF] Golub | Matrix computations[END_REF], pp. 68):

We denote such function with a double bar notation: f (x) = x . A useful class of vector norms are the p -norms defined by:

The 1-, 2-and ∞-norms are the most important:

A.3.2 Matrix Norm

f : K m×n → R is a matrix norm if the following properties hold (see [START_REF] Golub | Matrix computations[END_REF], pp. 71):

• f (A) ≥ 0, A ∈ K m×n , (f (A) = 0, iff A = 0).

We denote such function with a double bar notation: f (A) = A . The most frequently used matrix norms are the Frobenius norm:

and the p-norms:

In the case of p = 1 and p = ∞, we have:

|a ij |, which is the maximum absolute column sum of the matrix.

|a ij |, which is the maximum absolute row sum of the matrix.

Computation of Invariant Pairs and Matrix Solvents

A.4 Pseudo-inverse of a Matrix

Given an m × n matrix A, the Moore-Penrose generalized matrix inverse is a unique n × m matrix pseudo-inverse A + (see [START_REF] Golub | Matrix computations[END_REF], [START_REF] Penrose | A generalized inverse for matrices[END_REF]). The Moore-Penrose inverse satisfies:

A.5 The Generalized Schur Decomposition

The QZ algorithm of Moler and Stewart ([105], [START_REF] Golub | Matrix computations[END_REF]) solves the generalized eigenvalue problem of finding x and λ such that Ax = λBx, for A, B ∈ K n×n .

A.6 Composite Trapezoidal Rule

An intuitive method of finding the area under the curve y = f (x) over [a, b] is by approximating that area with a series of trapezoids that lie above the intervals {[x k , x k+1 ]} (see [START_REF] Burden | Numerical analysis[END_REF]).

Index

Computation of Invariant Pairs and Matrix Solvents