domain: math info

text:
Vanni Noferini
email: vanni.noferini@maths.manchester.ac.uk

Perez Alvaro
email: javier.perezalvaro@manchester.ac.uk

Javier Pérez

Chebyshev rootfinding via computing eigenvalues of colleague matrices: when is it stable?

Keywords: polynomial, roots, Chebyshev basis, matrix polynomial, colleague matrix, backward stability, polynomial eigenvalue problem, Arnold trasnversality theorem MSC classification: 65H04, 65H17, 65F15, 65G50

Computing the roots of a scalar polynomial, or the eigenvalues of a matrix polynomial, expressed in the Chebyshev basis {T k (x)} is a fundamental problem that arises in many applications. In this work, we analyze the backward stability of the polynomial rootfinding problem solved with colleague matrices. In other words, given a scalar polynomial p(x) or a matrix polynomial P (x) expressed in the Chebyshev basis, the question is to determine whether the whole set of computed eigenvalues of the colleague matrix, obtained with a backward stable algorithm, like the QR algorithm, are the set of roots of a nearby polynomial or not. In order to do so, we derive a first order backward error analysis of the polynomial rootfinding algorithm using colleague matrices adapting the geometric arguments in [A. Edelman and H. Murakami, Polynomial roots for companion matrix eigenvalues, Math. Comp. 210, 763-776, 1995] to the Chebyshev basis. We show that, if the absolute value of the coefficients of p(x) (respectively, the norm of the coefficients of P (x)) are bounded by a moderate number, computing the roots of p(x) (respectively, the eigenvalues of P (x)) via the eigenvalues of its colleague matrix using a backward stable eigenvalue algorithm is backward stable. This backward error analysis also expands on the very recent work [Y. Nakatsukasa and V. Noferini, On the stability of computing polynomial roots via confederate linearizations, To appear in Math. Comp.] that already showed that this algorithm is not backward normwise stable if the coefficients of the polynomial p(x) do not have moderate norms.

Introduction

A popular way to compute the roots of a monic polynomial expressed in the monomial basis is via the eigenvalues of its companion matrix. This is, for instance, the way followed by the MATLAB command roots, that, after balancing the companion matrix, uses the QR-algorithm to get its eigenvalues. The numerical properties of this method for computing roots of polynomials have been extensively studied [START_REF] De Terán | Backward stability of polynomial root-finding using Fiedler companion matrices[END_REF][START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF][START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF], in particular with respect to conditioning and backward errors. It has been shown that, in practice, if the companion matrix is balanced [START_REF] Parlett | Balancing a matrix for calculation of eigenvalues and eigenvectors[END_REF], the rootfinding method using companion matrices is numerically stable, in the sense that the computed roots are the exact roots of a nearby polynomial. However, as it was made famous by Wilkinson [21,[START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF][START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF], roots of polynomials whose roots lie on a real interval can be highly sensitive to perturbations in the coefficients when the monomial basis is used. So, even perturbations in the coefficients of order of the machine precision may produce a catastrophically large forward error. In practice, rootfinding on a real interval is a very frequent and important situation, and one way to circumvent this problem is to use, instead, a polynomial basis such that the roots of a polynomial expressed in that basis are better conditioned functions of its coefficients, like the Chebyshev basis.

Chebyshev polynomials are a family of polynomials, orthogonal with respect to the weight function w(x) = (1 -x 2 ) -1/2 on the interval [-1, 1], which may be computed using the following recurrence relation [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]Chapter 22]:

T 0 (x) = 1,
T 1 (x) = x, and, T k (x) = 2xT k-1 (x) -T k-2 (x), for k ≥ 2.

(1) Moreover, the Chebyshev polynomials T 0 (x), T 1 (x), . . . , T n (x) form a basis for the vector space of polynomials of degree at most n with real coefficients R n [x], so, any real polynomial p(x) ∈ R n [x] can be written uniquely as p(x) = n k=0 a k T k (x). Chebyshev polynomials are widely used in many areas of numerical analysis, and in particular approximation theory [START_REF] Trefethen | Approximation theory and approximation practice[END_REF]. In fact, a common approach, as done in Chebfun [START_REF] Trefethen | Chebfun Version 5. The Chebfun Development Team[END_REF], for computing the roots of a nonlinear smooth function f (x) on an interval is to approximate first f (x) by a polynomial p(x) expressed in the Chebyshev basis via Chebyshev interpolation and then compute the roots of p(x) as the eigenvalues of its colleague matrix [START_REF] Good | The colleague matrix, a Chebyshev analogue of the companion matrix[END_REF].

In this paper, we are interested in the backward stability of the rootfinding problem solved via colleague matrices and a backward stable eigenvalue algorithm. Our work is motivated by [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF], which addresses related issues for confederate matrices (the colleague matrix is a particular example of a confederate matrix [START_REF] Barnett | Polynomials and Linear Control Systems[END_REF][START_REF] Maroulas | Polynomials with respect to a general basis[END_REF]). Also, similar backward error analysis may be found in [START_REF] De Terán | Backward stability of polynomial root-finding using Fiedler companion matrices[END_REF][START_REF] Lawrence | Stability of rootfinding for barycentric Lagrange interpolants[END_REF][START_REF] Lawrence | Structured backward error analysis of polynomial eigenvalue problems solved by linearizations[END_REF]. In [START_REF] De Terán | Backward stability of polynomial root-finding using Fiedler companion matrices[END_REF], the authors study the backward stability of rootfinding methods using Fiedler companion matrices of monic polynomials expressed in the monomial basis; in [START_REF] Lawrence | Stability of rootfinding for barycentric Lagrange interpolants[END_REF], the authors study the backward stability of rootfinding methods using a suitable companion matrix of polynomials expressed in barycentric form; in [START_REF] Lawrence | Structured backward error analysis of polynomial eigenvalue problems solved by linearizations[END_REF], several bases are analyzed at once, for nonstandard linearizations of larger size with respect to the colleague or the companion.

Given a monic scalar polynomial in the Chebyshev basis of degree n p(x) = T n (x) + n-1 k=0 a k T k (x), with a k ∈ R, for k = 0, 1, . . . , n -1, where by monic in the Chebyshev basis we mean that the coefficient of T n (x) is equal to 1, the following matrix

C T = 1 2              -a n-1 -a n-2 + 1 -a n-3 • • • -a 2 -a 1 -a 0 1 0 1 0 • • • • • • 0 0 1 0
• • • • • • 0 1 0 1 0 • • • • • • • • • 0 2 0              ∈ R n×n ,
is known as the colleague matrix of p(x) [START_REF] Good | The colleague matrix, a Chebyshev analogue of the companion matrix[END_REF]. If we write T n (x) = ν n x n + • • • , the colleague matrix has the property that det(xI -C T ) = 1 νn p(x), so the roots of p(x) may be computed as the eigenvalues of C T using, for instance, the QR algorithm. The QR algorithm is a backward stable algorithm, this means than the computed eigenvalues are the exact eigenvalues of a matrix C T + E, where E is a (possibly dense) matrix such that E = O(u) C T , for some matrix norm, where u denotes the machine precision. However, the previous equation does not guarantee that the computed eigenvalues are the roots of a nearby polynomial of p(x) or, in other words, that this rootfinding method is backward stable. In order for the method to be backward stable (in a normwise sense), the computed eigenvalues should be the exact roots of a polynomial p(x

) = T n (x) + n-1 k=0 a k T k (x), such that p -p p = O(u),
for some polynomial norm. The backward stability of the polynomial rootfinding in degree-graded basis using confederate matrices is studied in [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF]. In particular (see [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF]Theorem 4.2]), the authors prove that if C T is the colleague matrix of a polynomial p(x) and E ∈ R n×n is any matrix, then the eigenvalues of C T + E are the exact roots of a polynomial p(x) such that

p(x) -p(x) = n-1 i=0 δ i (p, E)T i (x) + O( E 2 2 ), (2)
where, for i = 0, 1, . . . , n -1, the quantity δ i (p, E) is an affine function of the coefficients of p(x), and, separately, of the entries of E. Equation ( 2) implies that if the roots of p(x) are computed as the eigenvalues of its colleague matrix C T using a backward stable eigenvalue algorithm, then, the computed roots will be the exact roots of a polynomial p(x) such that

p -p p = κ(n)O(u) p ,
for some constant k(n). The previous equation shows, first, that this method is not backward stable if p ≫ 1, and, second, that this method is backward stable if the following two conditions are satisfied: (i) the quantity κ(n) is a low-degree polynomial in n with moderate coefficients; and, (ii) the norm p is moderate. As it is observed in [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF], writing δ i (p, E) = i,j,ℓ β ijℓ a ℓ E ij , since it is not clear what exactly are the constants β ijℓ involved, in principle it could happen that |β ijℓ | >> 1, implying that κ(n) might not be a polynomial in n with moderate coefficients. However, in this work we show that, in fact, |β ijℓ | ≤ 4, and, so,

p -p p = O(u) p ,
holds. The previous equation implies that computing the roots of p(x) via the eigenvalues of its colleague matrix using a backward stable eigenvalue algorithm is a backward stable rootfinding algorithm, provided that p 1.

Computing the eigenvalues of matrix polynomials in the Chebyshev basis is becoming an important problem [START_REF] Effenberger | Chebyshev interpolation for nonlinear eigenvalue problems[END_REF]. Given a p × p monic matrix polynomial in the Chebyshev basis of degree n

P (x) = I p T n (x) + n-1 k=0 A k T k (x), with A k ∈ R p×p , for k = 0, 1, . . . , n -1, (3)
where by monic in the Chebyshev basis we mean that the coefficient of T n (x) is equal to I p (the p × p identity matrix), the polynomial eigenvalue problem consists of finding the eigenvalues of P (x), that is, finding the roots of the scalar polynomial det (P (x)) (note that the monicity of P (x) implies its regularity), that is, det (P (x)) is not identically zero). A common approach to solve the polynomial eigenvalue problem for P (x) is to use the block colleague matrix

C T = 1 2              -A n-1 -A n-2 + I p -A n-3 • • • -A 2 -A 1 -A 0 I p 0 I p 0 • • • • • • 0 0 I p 0 I p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • 0 I p 0 I p 0 • • • • • • • • • 0 2I p 0              ∈ R np×np , (4)
since it is known (see [START_REF] Amiraslami | Linearization of matrix polynomials expressed in polynomial bases[END_REF]) that the eigenvalues of (4) coincide with the eigenvalues of P (x).

The backward stability question arises also in the matrix polynomial case, since the eigenvalues of C T may be computed using a backward stable eigenvalue algorithm, like the QR algorithm. As in the scalar polynomial case, the backward stability of the QR algorithm implies that the computed eigenvalues are the exact eigenvalues of the perturbed matrix C T + E, with E = O(u) C T . However, again, the backward stability of the QR algorithm does not guarantee that the computed eigenvalues are the eigenvalues of a nearby matrix polynomial of P (x) or, in other words, that the method is backward stable. In order for the method to be backward stable (in a normwise sense), the computed eigenvalues should be the exact eigenvalues of a matrix polynomial

P (x) = I p T n (x) + n-1 k=0 A k T k (x), with A k ∈ R p×p , such that P -P P = O(u),
for some matrix polynomial norm. Nevertheless, we will prove using some arguments inspired by [START_REF] Arnold | On matrices depending on parameters[END_REF][START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF][START_REF] Lemmonier | Optimal scaling of block companion pencils[END_REF] that if the eigenvalues of a matrix polynomial P (x) are computed as the eigenvalues of its colleague matrix using a backward stable eigenvalue algorithm, then, the computed eigenvalues are the exact eigenvalues of a monic matrix polynomial in the Chebyshev basis P (x) such that

P -P P = O(u) P .
The previous equation implies that this method is backward stable if P is moderate.

The paper is organized as follows. In Section 2 we present Arnold transversality theorem for colleague matrices, which will be the main tool to study the polynomial backward stability of the rootfinding method using colleague matrices. In Section 3, to give a flavor of our algebraic approach based on interpreting the companion or colleague linearizations as multiplication-by-x operators in certain quotient modules (see also [START_REF] Nakatsukasa | Vector spaces of linearizations for matrix polynomials: a bivariate polynomial approach[END_REF]Sec. 5]), we first review the backward error analysis of the rootfinding method using companion matrices in [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF]. Then, in Section 4 we prove Arnold transversality theorem for colleague matrices, and we use this theorem to study the backward stability of the rootfinding method using colleague matrices.

Throughout this paper, for a p × p matrix polynomial P (x) = n k=0 A k T k (x), non necessarily monic, P F is the norm on the vector space of p × p matrix polynomials of degree less than or equal to n defined as

P F = n k=0 A k 2 F .
Notice that, since we are going to deal with monic polynomials in the Chebyshev basis, A n = I p . Also notice that for a scalar polynomial p(x) = n k=0 a k T k (x), that is, for p = 1, this norm reduces to the usual 2-norm:

p F = p 2 = n-1 k=0 |a k | 2 .

Arnold transversality theorem for colleague matrices

Arnold transversality theorem will be the main tool in Section 4 to study what kind of polynomial backward stability is provided by matrix backward stability when the roots of scalar polynomials or the eigenvalues of matrix polynomials are computed as the eigenvalues of its colleague matrix with a backward stable eigenvalue algorithm. This theorem was first stated in [START_REF] Arnold | On matrices depending on parameters[END_REF] for companion matrices, and later generalized in [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF] to confederate matrices of scalar polynomials. Following [START_REF] Arnold | On matrices depending on parameters[END_REF][START_REF] De Terán | Backward stability of polynomial root-finding using Fiedler companion matrices[END_REF][START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF][START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF], we consider the Euclidian matrix space R n×n with the usual Frobenius inner product < A, B >:= tr (AB T ), where M T denotes the transpose of M ∈ R n×n . In this space, the set of matrices similar to a given matrix A ∈ R n×n is a differentiable manifold in R n×n . This manifold is called the orbit of A under the action of similarity: O(A) := {SAS -1 : S ∈ R n×n and det(S) = 0}.

A first-order expansion shows that the tangent space of O(A) at A is the set

T A O(A) := {AX -XA for some X ∈ R n×n }.
We also consider the vector subspace of "first block row matrices", denoted by BFR n,p ⊂ R np×np , which is defined as those n × n block matrices [X ij ], with X ij ∈ R p×p , whose block rows are all zero except (possibly) the first:

BFR n,p := X = I P 0 • • • 0 T X 1 X 2 • • • X n , for some X 1 , X 2 , . . . , X n ∈ R p×p ⊂ R np×np .
Note that taking p = 1 the space BFR n,p reduces to the vector subspace F R n of "first row matrices" introduced in [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF].

Arnold transversality theorem for a block colleague matrix C T of a monic matrix polynomial P (x) in the Chebyshev basis states that any matrix E ∈ R np×np may be decomposed as

E = F 0 + T,
where F 0 ∈ BFR n,p is a first block row matrix and T ∈ T CT O(C T ). Notice that taking p = 1 in the previous decomposition, this "block" version of Arnold transversality theorem reduces to a special case of [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF]Theorem 4.1].

In Section 4.2 we present a proof of Arnold transversality theorem, different to the one in [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF], extending (for the important case of the Chebyshev basis) [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF]Theorem 4.1] to the more complicated case of matrix polynomials. The new approach allows us to compute explicitly the matrix F 0 . Then, using this explicit expression, we study the polynomial backward stability of the rootfinding method using colleague matrices. However, before doing that, to give a flavor of our approach, we first treat the much easier case of the monomial basis (see [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF]).

Arnold transversality theorem for companion matrices and backward error analysis

The backward stability of the polynomial rootfinding using companion matrices is studied in [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF]. In this section we review their results in the matrix polynomial case, since results for scalar polynomials may be seen as corollaries of them. Given a p × p monic matrix polynomial expressed in the monomial basis

P (x) = I p x n + n-1 k=0 B k x k , with B k ∈ R p×p , for k = 0, 1, . . . , n -1,
its eigenvalues may be computed as the eigenvalues of the block companion matrix (see [START_REF] Lancaster | The Theory of Matrices[END_REF]Chapter 14])

C 1 :=         -B n-1 -B n-2 • • • -B 1 -B 0 I p 0 • • • 0 0 0 I p . . . 0 0 . . . . . . . . . . . . . . . 0 • • • 0 I p 0        
, using a backward stable eigenvalue algorithm.

In [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF], it is shown how a dense perturbation E to the block companion matrix leads to first order perturbations in the coefficients B k . To do so, they use Arnold transversality theorem for companion matrices, that is, they show how any matrix E ∈ R np×np may be decomposed as

E = E 0 + T = E 0 + C 1 X -XC 1 , (5)
where E 0 ∈ BFR n,p is a first block row matrix and

T = XC 1 -C 1 X ∈ T C1 O(C 1
) is a matrix in the tangent space of the orbit of C 1 . As a consequence of the previous decomposition we have that if E is a small perturbation of C 1 , then, to first order in E, the matrix C 1 + E is similar to C 1 + E 0 . In other words,

C 1 + E = (I + X) -1 (C 1 + E 0 )(I + X) + O( E 2 2 ),
for some matrix X ∈ R np×np . Noting that C 1 + E 0 is in turn a block companion matrix of another matrix polynomial, we see that a small perturbation of the block companion matrix of P (x) is similar, to first order in the norm of the perturbation, to a block companion matrix of a perturbed polynomial P (x).

In order to compute E 0 in ( 5), define the matrices

N k := n j=k C n-j 1 (I n ⊗ B j ), for k = 1, 2, . . . , n, (6)
where B n := I p and ⊗ denotes the Kronecker product, which satisfy, for k = 1, 2, . . . , n -1,

N k = N k+1 C 1 + I n ⊗ B k , with N n = I n ⊗ I p . (7)
Using the recurrence relation [START_REF] De Terán | Backward stability of polynomial root-finding using Fiedler companion matrices[END_REF], the matrix N k can be computed explicitly [8, Lemma 4.1]:

N k =             0 -B k-1 • • • -B 1 -B 0 . . . . . . . . . . . . . . . 0 -B k-1 • • • -B 1 -B 0 I p B n-1 • • • B k 0 I p B n-1 • • • B k . . . . . . . . . . . . . . . 0 I p B n-1 • • • B k             , for k = 1, 2, . . . , n, (8)
where the first block-row contains n -k block rows, and the second block row contains k block rows.

Writing the matrix E and the matrix N k as n×n block matrices E = [E ij ] and

N k = [(N k ) ij ], respectively, with E ij , (N k ) ij ∈ R p×p ,
the authors in [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF] prove that the matrix E 0 in ( 5) is given by

E 0 =      n i,j=1 E ij (N n ) ji • • • n i,j=1 E ij (N 2 ) ji n i,j=1 E ij (N 1 ) ji 0 • • • 0 0 . . . . . . . . . 0 • • • 0 0      .
Therefore, the matrix C 1 + E is similar, to first order in the norm of E, to a block companion matrix associated with the monic matrix polynomial

P (x) = I p x n + n-1 k=0 B k x k = I p x n + n-1 k=0 (B k - n ij=1 E ij (N k+1 ) ji ))x k .
From the previous equation, together with (8), we get that the eigenvalues of C 1 + E are the eigenvalues of a monic matrix polynomial

P (x) = I p x n + n-1 k=0 B k x k such that, to first order in the norm of E, B k -B k F = n s=k+1   n j=n-k B s E j-s+k+1,j   - k s=0   n-k-1 j=1 B s E j-s+k+1,j   F , (9)
for k = 0, 1, . . . , n -1. Then, if the eigenvalues of the block companion matrix C 1 are computed with a backward stable eigenvalue algorithm, it may be proved from ( 9) that, to first order in E, the computed eigenvalues are the exact eigenvalues of a polynomial P (x) such that

P -P F P F = O(u) P F . (10)
Note that [START_REF] Good | The colleague matrix, a Chebyshev analogue of the companion matrix[END_REF] implies that if the norm of the coefficients of P (x) are moderate, then computing its eigenvalues via the eigenvalues of its block companion matrix C 1 is a backward stable method.

In the following, we present a different approach to the matrices N k , for k = 1, 2, . . . , n, in (8) in order to better understand them and the block Toeplitz structure of their two blocks. First, by direct multiplication, it may be easily checked (see also [START_REF] Nakatsukasa | Vector spaces of linearizations for matrix polynomials: a bivariate polynomial approach[END_REF]Section 5] and [START_REF] Noferini | Duality of matrix pencils, Wong chains and linearizations[END_REF]Section 9]) that the block companion matrix C 1 satisfies

x      x n-1 . . . x 1      ⊗ I p = C 1      x n-1 . . . x 1      ⊗ I p + e 1 ⊗ P (x), (11)
where e 1 denotes the first column of the n × n identity matrix. Then, define the degree k Horner shift of P (x) as the matrix polynomial (see [6, Definition 4.1])

P k (x) := I p x k + k-1 j=0 B n-k+j x j . ( 12
)
The Horner shifts of P (x) satisfy, for k = 1, 2, . . . , n -1,

P k (x) = xP k-1 (x) + B n-k , with P 0 (x) = I p . (13)
Using ( 13), together with ( 7) and ( 11), it may be proved that the matrix N k in ( 6) is the unique matrix satisfying     

x n-1 . . .

x 1      ⊗ P n-k (x) = N k      x n-1 . . . x 1      ⊗ I p +      r 1k (x) . . . r n-1,k (x) r nk (x)      ⊗ P (x) (14)
for some scalar polynomials r 1k (x), . . . , r nk (x). In words, the ith block row of N k contains the coefficients of the unique matrix polynomial Q ik (x) of degree less than or equal to n -1 such that

x n-i P n-k (x) = Q ik (x) + r ik (x)P (x), (15)
for some scalar polynomial r ik (x).

Observing that the Horner shifts in (12) satisfy

x k P n-k (x) = P (x) -B k-1 x k-1 -• • • -B 1 x -B 0 , for k = 1, 2, . . . , n -1, (16)
it is immediate to check that the polynomial Q ik (x) in ( 15) is given by

Q ik (x) = n j=k B j x n+j-i-k if i ≥ n -k + 1, and k-1 j=0 -B j x n+j-i-k if i ≤ n -k, (17)
for i, k = 1, 2, . . . , n, where we set B n := I p . Then, from ( 15) and ( 17) we get

         
x n-1 . . .

x k x k-1 . . . 1           ⊗ P n-k (x) =            - k-1 j=0 B j x n-1+j-k . . . -k-1 j=0 B j x j n j=k B j x j-1 . . . n j=l B j x j-k            +           r 1k (x) . . . r n-k,k (x) r n-k+1,k (x)
. . .

r nk (x)           ⊗ P (x) (18)
and comparing ( 14) with [START_REF] Nakatsukasa | Vector spaces of linearizations for matrix polynomials: a bivariate polynomial approach[END_REF], we recover the expression in [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF] for the matrix N k , and the Toeplitz structure of its two blocks is also immediately explained. In Section 4.2 we will use an approach similar to [START_REF] Lemmonier | Optimal scaling of companion pencils for the QZ-algorithm[END_REF] in the case of the Chebyshev basis. However, there are more technicalities to be addressed with respect to the monomial basis. [START_REF] Barnett | Polynomials and Linear Control Systems[END_REF] Arnold transversality theorem for colleague matrices and backward error analysis

Clenshaw shifts and Clenshaw matrices

In this section we introduce some matrix polynomials and some matrices, named here as Clenshaw shifts and Clenshaw matrices, respectively, associated with a monic matrix polynomial in the Chebyshev basis P (x), that will be used through Section 4.2 and will be key in the following developments. Clenshaw shifts and Clenshaw matrices are the generalization of the Horner shifts in [START_REF] Lawrence | Stability of rootfinding for barycentric Lagrange interpolants[END_REF] and the matrices N 1 , N 2 , . . . , N k in ( 6), respectively, when the polynomial P (x) is expressed in the Chebyshev basis. Associated with the p × p monic matrix polynomial in the Chebyshev basis P (x) in (3), we define the following p × p matrix polynomials:

H 0 (x) = 2I p , H 1 (x) = 2xH 0 (x) + 2A n-1 , H k (x) = 2xH k-1 (x) -H k-2 (x) + 2A n-k , for k = 2, 3, . . . , n -2, H n-1 (x) = xH n-2 (x) -H n-3 (x)/2 + A 1 . (19)
We will refer, for k = 1, 2, . . . , n, to the matrix polynomial H k (x) as the degree k Clenshaw shifts of P (x), since for p = 1 they coincide with the well known Clenshaw shifts associated with a scalar polynomial expressed in the Chebyshev basis [START_REF] Clenshaw | A note on the summation of Chebyshev series[END_REF]. Clenshaw shifts are related with the polynomial P (x) through the following equation [START_REF] Clenshaw | A note on the summation of Chebyshev series[END_REF]:

2P (x) = 2xH n-1 (x) -H n-2 (x) + 2A 0 . (20)
In Theorem 4.1, given the Chebyshev polynomial T n-i (x) and the Clenshaw shift H n-k (x), we show how to express T n-i (x)H n-k (x) uniquely as Q ij (x) + r ik (x)P (x), where Q ij (x) is a p × p matrix polynomial of degree less than or equal to n -1 and r ik (x) is a scalar polynomial. This result is the Chebyshev analogue of [START_REF] Lemmonier | Optimal scaling of block companion pencils[END_REF] and [START_REF] Nakatsukasa | On the stability of computing polynomial roots via confederate linearizations[END_REF]. The proof of Theorem 4.1 is elementary but rather technical, so we leave it to the appendix. In order to write down a reasonably simple formula for T n-i (x)H n-k (x), we define the following quantities

Γ 2k+1 = Γ 2k-1 + 2A n-2k-1 , for k = 1, 2, . . . , n 2 -1, with Γ 0 = 2I p , and
Γ 2k = Γ 2(k-1) + 2A n-2k , for k = 1, 2, . . . , n 2 -1, with Γ 1 = 2A n-1 . (21)
Notice that in Γ k only appear coefficients of P (x) with indices of the same parity.

Theorem 4.1.

Let P (x) = I p T n (x) + n-1 k=0 A k T k (x)
be a p × p monic matrix polynomial in the Chebyshev basis of degree n, let T n-i (x) and H n-k (x) be, respectively, the degree n -i Chebyshev polynomial and the degree n -k Clenshaw shift of P (x), with i, k ∈ {1, 2, . . . , n}. Then, there exist a unique p × p matrix polynomial Q ik (x) of degree less than or equal to n -1 and a unique scalar polynomial r ik (x) such that

T n-i (x)H n-k (x) = Q ik (x) + r ik (x)P (x),
where,

• if i ≥ n -k + 1 and k ≥ 2, Q ik (x) = n-k-1 ℓ=0 Γ ℓ (T 2n-i-k-ℓ (x) + T |k+ℓ-i| (x)) + Γ n-k T n-i (x); (22)
• if i = n and k = 1, Q ik (x) = n-2 ℓ=0 Γ ℓ T n-1-ℓ (x) + Γ n-1 2 T 0 (x); (23)
• if i ≤ n -k and n -1 ≥ k ≥ 2, Q ik (x) = i-2 ℓ=0 Γ ℓ (T i+k-2-ℓ (x) + T |k+ℓ-i| (x)) + Γ i-1 T k-1 (x) - n-k+1-i ℓ=1 k-1+ℓ r=1 2A k-1+ℓ-r T |n-i+1-ℓ-r| (x); (24)
• if i ≤ n -k and k = 1 Q ik (x) = i-2 ℓ=0 Γ ℓ T i-1-ℓ (x) + Γ i-1 2 T 0 (x) - n-i ℓ=1 ℓ r=1 A ℓ-r T |n-i+1-ℓ-r| (x); ( 25
)
where Γ ℓ , for ℓ = 0, 1, 2, . . ., is defined in [START_REF] Peters | Practical problems arising in the solution of polynomial equations[END_REF].

From Theorem 4.1, it is clear that there exists a unique n × n block matrix

M k = [(M k ) ij ], with (M k ) ij ∈ R p×p , such that      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H n-k (x) = M k      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ I p +      r 1k (x) . . . r n-1,k (x) r nk (x)      ⊗ P (x), for k = 1, 2, . . . , n, (26)
where ⊗ denotes the Kronecker product, for some scalar polynomials r 1k (x), . . . , r nk (x). We will refer to the matrix M k in [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF] as the kth Clenshaw matrix of P (x).

By direct multiplication, it may be easily checked that the block colleague matrix

C T satisfies x     
T n-1 (x) . . .

T 1 (x) T 0 (x)      ⊗ I p = C T     
T n-1 (x) . . .

T 1 (x) T 0 (x)      ⊗ I p + e 1 ⊗ P (x), (27)
Then, using ( 26) and ( 27), in Proposition 4.2 we show that the Clenshaw matrices M 1 , M 2 , . . . , M n in (26) satisfy a simple recurrence relation. This recurrence relation is the analogue of (13) when the matrix polynomial P (x) is expressed in the Chebyshev basis instead of the monomial basis.

Proposition 4.2. Let P (x) = I p T n (x)+ n-1 k=0 A k T k (x)
be a p×p monic matrix polynomial in the Chebyshev basis of degree n, let C T be the block colleague matrix of P (x), and let M 1 , M 2 , . . . , M n be the Clenshaw matrices in [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF]. Then,

M n = I n ⊗ 2I p , M n-1 = 2M n C T + I n ⊗ 2A n-1 , M k = 2M k+1 C T -M k+2 + I n ⊗ 2A k , for k = n -2, . . . , 3, 2,
and

M 1 = M 2 C T -M 3 /2 + I n ⊗ A 1 .
(28)

Proof. The proof proceeds by induction on k (backwards from k = n). First, we prove that the result is true for k = n. From (19), we have

     T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H 0 (x) =      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ 2I p =      2I p . . . 2I p 2I p           T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ I p .
Comparing the previous equation with [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF], we deduce that M n = I n ⊗ 2I p .

Second, we prove that the result is true for k = n -1. From ( 19), we have

     T n-1 (x)
. . .

T 1 (x) T 0 (x)      ⊗ H 1 (x) =      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ (2xH 0 (x) + 2A n-1 ) = 2x      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H 0 (x) +      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ 2A n-1 .
Using the inductive hypothesis, together with (27), we get

     T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H 1 (x) = (2M n C T + I n ⊗ 2A n-1 )      T n-1 (x) . . . T 1 (x) T 0 (x)      +      r 1,n-1 (x) . . . r n-1,n-1 (x) r n,n-1 (x)      ⊗ P (x),
for some scalar polynomials r 1,n-1 (x), . . . , r n,n-1 (x). Comparing the previous equation with ( 26), we deduce that

M n-1 = 2M n C T + I n ⊗ 2A n-1 .
Third, assume that the result is true for

M n , M n-1 , . . . , M k+1 , with n -2 ≥ k ≥ 2. From (19), we have      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H n-k (x) =      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ (2xH n-k-1 (x) -H n-k-2 (x) + 2A k ) = 2x      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H n-k-1 (x) -      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H n-k-2 (x) +      T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ 2A k .
Using that the result is true for k + 1 and k + 2, together with (27), we get

     T n-1 (x) . . . T 1 (x) T 0 (x)      ⊗ H n-k (x) = (2C T M k+1 -M k+2 + I n ⊗ 2A k )      T n-1 (x) . . . T 1 (x) T 0 (x)      +      r 1k (x) . . . r n-1,k (x) r n,k (x)      ⊗ P (x),
for some scalar polynomials r 1k (x), . . . , r nk (x). Comparing the previous equation with [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF], we deduce that

M k = 2C T M k+1 -M k+2 + I n ⊗ 2A k .
Finally, the proof of the last case (k = 1) is similar to the proof for the previous cases (n

-2 ≥ k ≥ 2), but using H n-1 (x) = xH n-2 (x) -H n-3 (x)/2 + A 1 , so we omit it.
In contrast with (8), the Clenshaw matrices M 1 , M 2 , . . . , M n have a complicated structure. For example, for n = 6 and k = 3, it is easy to check using (28) that the matrix M k is equal to

        0 -2A 2 -2A 3 -2A 1 2I p -2A 4 -2A 2 -4A 0 -2A 3 -4A 1 -2A 2 -2A 0 0 0 2I p -2A 2 2A 5 -2A 3 -2A 1 2I p -2A 2 -4A 0 -2A 1 0 2I p 2A 5 2I p + 2A 4 -2A 2 2A 5 -2A 1 2I p -2A 0 2I p 2A 5 2I p + 2A 4 2A 5 + 2A 3 4I p + 2A 4 2A 5 0 2I p 2A 5 4I p + 2A 4 4A 5 + 2A 3 2A 4 + 2I p 0 0 4I p 4A 5 4I p + 4A 4 2A 5 + 2A 3         . ( 29
)
Two observations about the matrix (29) are: (i) its first block column is equal to e n-k+1 ⊗ 2I p , where e ℓ denotes the ℓth column of the n × n identity matrix; and, (ii) if we set A n := I p , each block entry has the form n i=0 c i A i , where |c i | ≤ 4. In Theorem 4.3, we show that the two previous observations are true for any n and k. Property (i) will be key to prove Arnold transversality theorem, and property (ii) will be key to study what kind of backward stability of a linearization-based algorithm for the polynomial eigenvalue problem is provided by the backward stability of an eigensolver for the linearized problem. Theorem 4.3. Let P (x) = I p T n (x) + n-1 k=0 A k T k (x) be a p × p monic matrix polynomial in the Chebyshev basis of degree n, and let M k , for k = 1, 2, . . . , n, be the kth Clenshaw matrix in [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF]. Then, the following statements hold:

(a) The first block column of M k is equal to e n-k+1 ⊗ 2I p , where e ℓ denotes the ℓth column of the n × n identity matrix.

(b) For i, j = 1, 2, . . . , n, the (i, j)th block entry of

M k satisfies (M k ) ij = n t=0 α t,
H n-k (x)T n-i (x) = n j=1 (M k ) ij T n-j (x) + r ik (x)P (x).
Therefore, to prove part (a) it is enough to show that

T n-i (x)H n-k (x) = 2I p T n-1 (x) + • • • + r ik (x)P (x), ( 30
) if i = n -k + 1
, and that

T n-i (x)H n-k (x) = (M k ) iν T ν (x) + • • • + r ik (x)P (x), (31)
with ν < n -1, if i = n -k + 1
, where the dots correspond to Chebyshev polynomials with lower indices. First, suppose that i ≥ n-k+1. We will prove that T n-i (x)H n-k (x) is of the form (30) when i = n-k+1 and it is of the form (31) otherwise. We need to distinguish several cases. First, let k = n. From [START_REF] Trefethen | Approximation theory and approximation practice[END_REF] we get that T n-i (x)H 0 (x) = Γ 0 T n-i (x) = 2I p T n-i (x). Since the index n -i is equal to n -1 if and only if i = 1, the result is true in this case. Then, consider the case n -1 ≥ k ≥ 2. There are three kinds of indices of Chebyshev polynomials in [START_REF] Trefethen | Approximation theory and approximation practice[END_REF]. The first is 2n -i -k -ℓ, which is equal to n -1 if and only if ℓ = 0 and i = n -k + 1. This gives a contribution Γ 0 T n-1 (x) = 2I p T n-1 (x) only when i = n -k + 1. The second one is |k + ℓ -i|. Taking into account the possible values that k, ℓ, and i can take in [START_REF] Trefethen | Approximation theory and approximation practice[END_REF], it may be easily checked that this index is smaller than or equal to n -2. The third index is n -i which necessarily is smaller than or equal to n -2, and, hence, the result is true in this case. Finally, consider the case k = 1 and i = n. There are two kinds on indices of Chebyshev polynomials in [START_REF] Trefethen | Chebfun Version 5. The Chebfun Development Team[END_REF]. The first one is n -1 -ℓ, which is equal to n -1 if and only if ℓ = 0. This gives a contribution Γ 0 T n-1 (x) = 2I p T n-1 (x). The second index is 0, which is smaller than n -2. Therefore, the result is also true in this case. Now suppose that i ≤ n -k. We will we prove that T n-i (x)H n-k (x) is of the form (31). Notice that there are four kinds of indices in [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF] when k ≥ 2, namely, i + k -2 -ℓ, |k + ℓ -i|, |n -i + 1 -ℓ -r| and k -1, and three kinds on indices in [START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF] when k = 1, namely, i -1 -ℓ, i -1 and |n -i + 1 -ℓ -r|. Taking into account the possible values that k, ℓ, r, and i can take in (31), in both cases (k ≥ 2 and k = 2), it may be checked that these indices do not exceed n -2. Now, we proceed to prove part (b). Again, we need to distinguish several cases. First, suppose that i ≥ n -k + 1 and also assume that k ≥ 2 (the argument when k = 1 is similar and simpler, so we omit it), and consider the three kinds of indices of Chebyshev polynomials that appear in [START_REF] Trefethen | Approximation theory and approximation practice[END_REF], namely, 2n -i -k -ℓ, |k + ℓ -i|, and n -i. For ℓ = 0, 1, . . . , n -k, a careful look at these indices reveals that if k + ℓ -i ≥ 0, then the three of them are different. Therefore, we can write [START_REF] Trefethen | Approximation theory and approximation practice[END_REF] as

n-1 ℓ=0 B ℓ T ℓ (x) + -1 ℓ=-(1-n) B ℓ T -ℓ (x), ( 32
)
where B ℓ is equal to either 0 or Γ t for some t. It follows that (M k ) ij is equal to either 0, Γ t for some t, or Γ t1 + Γ t2 for some t 1 , t 2 . Finally, recall from (21) that Γ t is equal to 2I

p + 2A n-2 + 2A n-4 + • • • if t is even, or to 2A n-1 + 2A n-3 + • • • if t is odd. Therefore, (M k ) ij = n t=0 α t,ijk A t , with |α t,ijk | ≤ 4.
Then suppose that i ≤ n -k and also assume that k ≥ 2 (again, the argument when k = 1 is similar and simpler, so we omit it). First, consider the three kinds of indices of Chebyshev polynomials that appear in the first summand in [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF], namely, i + k -2 -ℓ, |k + ℓ -i|, k -1. For ℓ = 0, 1, . . . , i -2, again, it may be checked that if k + ℓ -i ≥ 0, then these three indices are different. Therefore, the first summand in [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF] is also of the form (32). Finally, consider the index of the Chebyshev polynomials and the index of the coefficients A i that appear in the second summand in [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF], namely, |n -i + 1 -ℓ -r|, and k + 1 + ℓ -r. If n -i + 1 -ℓ -r ≥ 0, it may be checked that for any two allowed different pairs (ℓ, r) that realize the same value of n -i + 1 -ℓ -r, then the associate indices k + 1 + ℓ -r must be different. Since the same occur when n -i + 1 -ℓ -r < 0, it follows that ( 24) is of the form

n-2 ℓ=0 C ℓ T ℓ (x) + -1 ℓ=-(2-n) C ℓ T -ℓ (x) -2 n-2 ℓ=0 D ℓ T ℓ (x) -2 -1 ℓ=2-n D ℓ T -ℓ (x) + r ik (x)P (x)
where C ℓ is equal to either 0 or Γ t for some t, and D ℓ is equal to

q ℓ t=1 A it , where i t1 = i t2 whenever t 1 = t 2 .
Then, it follows that

(M k ) ij = n ℓ=0 δ ℓ A ℓ - n ℓ=0 ρ ℓ A ℓ ,
where δ ℓ and ρ ℓ are equal to either 4, or 2 or 0, therefore

(M k ) ij = n t=0 α t,ijk A t with |α t,ijk | ≤ 4.
If necessary, explicit expressions of the entries of the Clenshaw matrices M 1 , M 2 , . . . , M n may be obtained from Theorem 4.1. However, since Theorem 4.3 is the only information that we will need about them to prove our main results in the following section, we do not pursue that idea.

Backward error of the Chebyshev rootfinding method using colleague matrices

In this section we prove Arnold transversality theorem for colleague matrices of monic polynomials in the Chebyshev basis. That is, we show that any matrix E ∈ R pn×pn may be decomposed as

E = F 0 + T, (33)
where F 0 ∈ BFR n,p is a first block row matrix and T ∈ T CT O(C T ), constructing the matrix F 0 explicitly. As in the case of the monomial basis, generically, dim (T CT O(C T )) + dim (BFR n,p ) = n 2 p 2 -np + np 2 ≥ n 2 p 2 (see [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of block companion pencils[END_REF]). In words, the tangent space T CT O(C T ) and the vector space of first block row matrices BFR n,p may have a nontrivial intersection when p > 1. For this reason, following [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF][START_REF] Lemmonier | Optimal scaling of block companion pencils[END_REF] we choose a particular subspace of the tangent space that will give a unique decomposition (33). This subspace is denoted by Sub T CT O(C T ) and it is given by

Sub T CT O(C T ) = {X ∈ T CT O(C T ) such that X has 0 first block column}.
In order to get the decomposition (33) with T ∈ Sub T CT O(C T ) we will make use of the Clenshaw matrices M 1 , M 2 , . . . , M n ∈ R np×np , defined in [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF], of the matrix polynomial P (x) in [START_REF] Arnold | On matrices depending on parameters[END_REF]. Recall that the Clenshaw matrices satisfy the following recurrence relation (see Proposition 4.2):

M n = I n ⊗ 2I p , M n-1 = 2M n C T + I n ⊗ 2A n-1 , M k = 2M k-1 C T -M k-2 + I n ⊗ 2A k , for k = n -2, . . . , 3, 2, and M 1 = M 2 C T -M 3 /2 + I n ⊗ A 1 . ( 34
)
We illustrate these matrices with a small example (n = 4). From (34), it is easy to check that

M 4 =     2I n 0 0 0 0 2I n 0 0 0 0 2I n 0 0 0 0 2I n     , M 3 =     0 2I n -2A 2 -2A 1 -2A 0 2I n 2A 3 2I n 0 0 2I n 2A 3 2I n 0 0 4I n 2A 3     , M 2 =     0 -2A 1 2I n -2A 2 -4A 0 -2A 1 0 2I n 2A 3 -2A 1 2I n -2A 0 2I n 2A 3 4I n + 2A 2 2A 3 0 4I n 4A 3 2I n + 2A 2     , and
M 1 =     0 -2A 0 -2A 1 I n -A 2 -A 0 0 0 2I n -2A 0 A 3 -A 1 0 2I n 2A 3 I n + A 2 -A 0 2I n 2A 3 2I n + 2A 2 A 3 + A 1     .
The only information that we need about Clenshaw matrices are those stated in Theorem 4.3.

Following [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF], we also define the block trace of a np Proof. From part (a) in Theorem 4.3, the (n -k, 1) block entry of M k+1 is equal to 2I p , and the (i, 1) block entry of M k+1 , with i = n -k, is equal to 0. Therefore, Z 1,n-k , for k = 0, 1, . . . , n -1, is uniquely determined from (35). Also, if X has 0 first block column, it may be easily checked that the map from X to the last n -1 block rows of C T X -XC T has a trivial nullspace. Thus, Z is uniquely determined by (36).

× np block matrix Z = [Z ij ], with Z ij ∈ R p×p ,
To finish the proof we need to prove that (36) implies (35). That is, we need to show that tr p (M k+1 (C T X-XC T )) = 0 for any block matrix X with 0 first block column. In order to do this, first, we show that if X has 0 first block column, then tr p (M k+1 XC T ) = tr p (XC T M k+1 ). The proof of the previous equation is not completely immediate when p > 1 since, in this situation, tr p (AB) = tr p (BA) does not hold in general. So, consider a block matrix Y that has 0 first block column. Then,

tr p (C T Y ) = p-2 i=1 Y i,i+1 2 + Y i+2,i+1 2 + Y p-1,p = tr p (Y C T ).
Therefore, if X has 0 first block column, then tr p (M k+1 XC T ) = tr p (C T M k+1 X).

Then, we show that tr p (C T M k+1 X) = tr p (M k+1 C T X). To do this, note that the Clenshaw matrix

M k+1 is of the form 2 n-k C n-k-1 T + n-k-1 t=1 (I n ⊗ B k )C n-k-1-t T , for some B 1 , B 2 , . . . , B n-k-1 ∈ R p×p (
this can be verified by induction using (34)). So, we only need to show that tr p (C T (I n ⊗B)C j T X) = tr p ((I ⊗B)C j T C T X). Indeed, since the matrix C T (I n ⊗ B) -(I n ⊗ B)C T is 0 except the first block row, and since C T X has 0 first block column, it follows that tr p (C

T (I n ⊗ B)C j T X -(I ⊗ B)C j T C T X) = 0. Therefore, tr p (C T M k+1 X) = tr p (M k+1 C T X). Thus, we conclude that tr p (M k+1 XC T ) = tr p (C T M k+1 X) = tr p (M k+1 C T X).
In Theorem 4.5 we present the proof of Arnold transversality theorem for block colleague matrices. Part (a) in Theorem 4.3 will be key here.

Theorem 4.5. Let p(x) = I p T n (x) + n-1 k=0 A k T k (x)
be a p × p monic matrix polynomial in the Chebyshev basis of degree n, and let C T be its block colleague matrix. Then, any matrix E ∈ R np×np can be expressed as

E = F 0 + T, (37)
where F 0 ∈ BFR n,p is a first block row matrix, and T ∈ Sub T CT O(C T ). Moreover, if the first block row of

F 0 is written as F (n-1) 0 • • • F (1) 0 F (0) 0
, then

F (k) 0 = 1 2 tr p (EM k+1 ), for k = 0, 1, . . . , n -1, ( 38
)
where the matrix M k+1 is Clenshaw matrix defind in [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF].

Proof. Define F (k) 0 = 1 2 tr p (EM k+1
), for k = 0, 1, . . . , n -1, and let F 0 ∈ BFR n,p be a first block row matrix such that its first block row is

F (n-1) 0 • • • F (1) 0 F (0) 0
. We may write the matrix T := E -F 0 . Then, we have to check that T ∈ Sub T CT O(C T ). From Theorem 4.4, we see that it is sufficient to show that tr p (T M k+1 ) = 0, for k = 0, 1, . . . , n -1. Indeed, using part (a) in Theorem 4.3,

tr p (T M k+1 ) = tr p (EM k+1 ) -tr p (F 0 M k+1 ) = tr p (EM k+1 ) 2F (k) 0 = tr p (EM k+1 ) -tr p (EM k+1 ) = 0, for k = 0, 1, . . . , n -1. So, we conclude that T ∈ Sub T CT O(C T ).
An important consequence of the decomposition in Theorem 4.5 is that if E is a small perturbation of the block colleague matrix C T , then

C T + E = C T + F 0 + T = C T + F 0 + (C T X -XC T ) = (I + X) -1 (C T + F 0 )(I + X) + O( E 2 2 ),
where we have used that T can be written as C T X -XC T , for some X ∈ R np×np with 0 first block column.

Noticing that C T + F 0 is in turn a block colleague matrix of another matrix polynomial, we deduce that a small perturbation of the block colleague matrix of P (x) is similar, to first order in the norm of the perturbation, to a block colleague matrix of a perturbed polynomial P (x). This observation allows us to formulate the next corollary.

Corollary 4.6. Let P (x) = I p T n (x) + n-1 k=0 A k T k (x) be a p × p monic matrix polynomial in the Chebyshev basis of degree n, and let C T be its block colleague matrix. Assume that the eigenvalues of P (x) are computed as the eigenvalues of C T with a backward stable algorithm, i.e., an algorithm that computes the exact eigenvalues of some matrix C T + E, with E F = O(u) C T F , where u is the machine precision. Then, to first order in u, the computed roots are the exact roots of a polynomial P (x) such that

P -P F P F = O(u) P F .
Proof. If a backward stable eigensolver is given C T as an input, the computed eigenvalues are the exact eigenvalues of a matrix C T + E, for some E with E F = ǫ C T F , where ǫ = uh(n), for some low degree polynomial h with moderate coefficients. In other words, the computed eigenvalues are the exact roots of the polynomial det(xI -C T -E). Using Theorem 4.5, we can write E = F 0 + T , where T ∈ Sub T CT O(C T ) and F 0 is a first block row matrix with first block row as in (38). Therefore,

C T + E = C T + F 0 + C T X -XC T = (I + X) -1 (C T + F 0 )(I + X) + O(u 2 ),
for some matrix X ∈ R np×np with 0 first block column, that is, the matrix C T + E is similar, to first order in u, to the colleague matrix C T + F 0 . Writing E as a np × np block matrix E = [E ij ], with E ij ∈ R p×p , and noticing that C T + F 0 is the colleague matrix of the matrix polynomial

P (x) = I p T n (x) + n-1 k=0 (A k - F (k) 0 )T k (x)
, we have that, to first order in u, the computed eigenvalues are the exact eigenvalues of a matrix polynomial P

(x) = I p T n (x) + n-1 k=0 A k T k (x), with A k -A k F = F (k) 0 F = tr p (EM k+1 ) F = n i,j=1 E ij (M k+1 ) ji F . Therefore, for k = 0, 1, . . . , n -1, we have A k -A k F ≤ n i,j=1 E ij F (M k+1 ) ji F ≤ n i,j=1 E ij 2 F n i,j=1 (M k+1 ) ij 2 F = E F M k+1 F ,
where we have used the Cauchy-Schwarz inequality. Then, using part (b) of Theorem 4.3, we have

M k+1 F = n i,j=1 (M k+1 ) ij 2 F = n i,j=1 n t=0 α t,ij,k+1 A t 2 F ≤ n i,j=1 n t=0 α t,ij,k+1 A t F 2 ≤ 4 n i,j=1 n t=0 A t F 2 ≤ 4n 3/2 P F ,
where we have used n t=0 A t 2 ≤ √ n P 2 . Finally, using that E F = ǫ C T F , we get that

E F = ǫ C T F ≤ ǫ 2         -A n-1 -A n-1 • • • -A 0 0 0 • • • 0 . . . . . . . . . . . . . . . . . . 0 0 • • • 0         F + ǫ 2        I p I p I p . . . . . . I p 0 I p 2I p        F ≤ ǫ 2   n-1 t=0 A t 2 F + 2np   ≤ ǫ √ 2n 2   n-1 t=0 A t 2 F + √ p   ≤ ǫ √ 2n P F .
Thus, the computed eigenvalues, to first order in u, are the exact eigenvalues of a monic matrix polynomial in the Chebyshev basis P (x) such that,

P -P F = n-1 k=0 A k -A k 2 F ≤ n-1 k=0 A k -A k F ≤ n-1 k=0 M k+1 F E F ≤ 4n 5/2 P F E F ≤ ǫ P 2 F ,
where ǫ = u h(n), for some low degree polynomial h with moderate coefficients.

Conclusions

In this paper, we have analyzed the backward stability of a Chebyshev-basis polynomial rootfinder (or matrix polynomial eigensolver) based on the solution of the standard eigenvalue problem for the corresponding colleague matrix. More precisely, given a monic scalar polynomial in the Chebyshev basis p(x), we have proved that if the roots of p(x) are computed as the eigenvalues of a colleague matrix using a backward stable eigenvalue algorithm, like the QR algorithm, then the computed roots are the exact roots of a monic polynomial in the Chebyshev basis p(x) such that

p -p 2 p 2 = O(u) p 2 ,
Similarly, if the eigenvalues of a monic matrix polynomial in the Chebyshev basis are computed as the eigenvalues of a block colleague matrix using a backward stable eigenvalue algorithm, then the computed eigenvalues are the exact eigenvalues of a monic matrix polynomial in the Chebyshev basis P (x) such that

P -P F P F = O(u) P F ,
These backward error analysis show that these methods are backward stable when the norms p 2 and P F are moderate.

A Proof of Theorem 4.1

In this section we present the proof of Theorem 4.1, that is, given the Clenshaw shift H n-k (x) associated with the matrix polynomial P (x) in (3), and the Chebyshev polynomial T n-i (x), we show that

T n-i (x)H n-k (x) = Q ik (x) + r ik (x)P (x), (39)
for some scalar polynomial r ik (x), where Q ik (x) is the matrix polynomial of degree less than or equal to n -1 in ( 22)- [START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF]. Moreover, we show that the decomposition (39) is unique.

Along the proof, quite often products of two of Chebyshev polynomials will occur. For this reason, the following formula [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]Chapter 22] is of fundamental importance here:

2T m (x)T n (x) = T m+n (x) + T |m-n| (x). ( 40
)
The first step is to expand the Clenshaw shifts H k (x), for k = 0, 1, . . . , n -1, in the Chebyshev basis. We will prove

H k (x) = k-1 ℓ=0 2Γ ℓ T k-ℓ (x) + Γ k T 0 (x), for k = 0, 1, . . . , n -2, and (41)
H n-1 (x) = n-2 ℓ=0 Γ ℓ T n-1-ℓ (x) + 1 2 Γ n-1 T 0 (x), (42)
where Γ ℓ is defined in [START_REF] Peters | Practical problems arising in the solution of polynomial equations[END_REF]. The proof proceeds by induction on k. From [START_REF] Noferini | Duality of matrix pencils, Wong chains and linearizations[END_REF] we get H 0 (x) = 2I p = Γ 0 T 0 (x) and

H 1 (x) = 4I p x + 2A n-1 = 2Γ 0 T 1 (x) + Γ 1 T 0 (x)
, so the result is true for k = 0 and k = 1. Then, assume that the result is true for H 0 (x), H 1 (x), . . . , H k-1 (x), with 2 ≤ k ≤ n -2. Using the induction hypothesis, together with [START_REF] Noferini | Duality of matrix pencils, Wong chains and linearizations[END_REF], we have

H k (x) = 2xH k-1 (x) -H k-2 (x) + 2A n-k = 2x k-2 ℓ=0 2Γ l T k-1-ℓ (x) + Γ k-1 T 0 (x) - k-3 ℓ=0 2Γ ℓ T k-2-ℓ (x) -Γ k-2 T 0 (x) + 2A n-k .
Using T 0 (x) = 1, T 1 (x) = x, and (40) with m = 1 and n = k, from the previous equation we get

H k (x) = k-2 ℓ=0 2Γ ℓ (T k-ℓ (x) + T k-2-ℓ (x)) + 2Γ k-1 T 1 (x) - k-3 ℓ=0 2Γ ℓ T k-2-ℓ (x) -Γ k-2 T 0 (x) + 2A n-k T 0 (x) = k-2 ℓ=0 2Γ ℓ T k-ℓ (x) + 2Γ k-2 T 0 (x) + 2Γ k-1 T 1 (x) -Γ k-2 T 0 (x) + 2A n-k T 0 (x) = k-1 ℓ=0 2Γ ℓ T k-ℓ (x) + (Γ k-2 + 2A n-k )T 0 (x) = k-1 ℓ=0 2Γ ℓ T k-ℓ (x) + Γ k T 0 (x),
where in the last equality we have used Γ k-2 + 2A n-k = Γ k . Therefore, the result is also true for H k (x). Finally, the proof that (42) holds is similar to the previous one, but starting with H n-1 (x) = xH n-2 (x) -H n-3 (x)/2 + A 1 , so we omit the details. Now we proceed to show that (39) holds with Q ik (x) as in ( 22)- [START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF]. In order to do that, we will proceed in certain order. To help the reader to follow the steps, we depict all the possible products T n-i (x)H n-k (x) for n = 10 in the following 10 × 10 grid.

From ( 1) and ( 19), we have

H n-2 (x)T 2 (x) = H n-2 (x) (2xT 1 (x) -T 0 (x)) = 2xT 1 (x)H n-2 (x) -T 0 (x)H n-2 (x) = T 1 (x) (2H n-1 (x) + H n-3 (x) -2A 1 ) -T 0 (x)H n-2 (x) = 2T 1 (x)H n-1 (x) + T 1 (x)H n-3 (x) -T 0 (x)H n-2 (x) -2A 1 T 1 (x).
As can be seen from the previous equation, the product H n-2 (x)T 2 (x) may be computed from products represented by two triangular vertices: T 1 (x)H n-3 (x) and T 0 (x)H n-2 (x), and the product T 1 (x)H n-1 (x). Then, using (41), (42), and the result previously obtained for T 1 (x)H n-1 (x), we get

T 2 (x)H n-2 (x) = n-4 ℓ=0 Γ ℓ (T n-2-ℓ (x) + T |ℓ+4-n| (x)) + Γ n-3 T 1 (x) -2A 0 T 0 (x) -2A 1 T 1 (x) + • • • ,
where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. The previous equation corresponds to [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF] with i = n -2 and k = 2. Finally, we consider the white circular vertices labeled with 3, that is, the products T k (x)H n-k (x), for k = 3, 4, . . . , n. From ( 1) and ( 19), we have

T k (x)H n-k (x) = (2xT k-1 (x) -T k-2 (x))H n-k (x) = 2xT k-1 (x)H n-k (x) -T k-2 (x)H n-k (x) = T k-1 (x)(H n-k+1 (x) + H n-k-1 (x) -2A k-1 ) -T k-2 (x)H n-k (x) = T k-1 (x)H n-k+1 (x) + T k-1 (x)H n-k-1 (x) -T k-2 (x)H n-k (x) -2A k-1 T k-1 (x).
As can be seen from the previous equation,

T k (x)H n-k (x) may be computed from T k-1 (x)H n-k-1 (x) and T k-2 (x)H n-k (x)
, represented in the grid by triangular vertices, and T k-1 (x)H n-k+1 (x), represented in the grid by the white circular vertex above the white circular vertex corresponding to T k (x)H n-k (x). Since we have previously seen that Theorem 4.1 holds for T 1 (x)H n-1 (x) and T 2 (x)H n-2 (x), and for products represented by triangular vertices, this shows how to prove inductively (from top to bottom) that Theorem 4.1 holds for products represented by white circular vertices labeled with 3. Indeed, assuming that the result holds for T k-1 (x)H n-k+1 (x) and using [START_REF] Trefethen | Approximation theory and approximation practice[END_REF], we get

T k (x)H n-k (x) = T k-1 (x)H n-k-1 + k r=2 (-2A k-r )T k-r (x) -2A k-1 T k-1 (x) + • • • = n-k-2 ℓ=0 Γ ℓ (T n-2-ℓ (x) + T |ℓ+2k-n| ) + Γ n-k-1 T k-1 (x) + k r=1 (-2A k-r )T k-r (x) + • • • ,
where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. It is immediate to check that the previous equation corresponds to [START_REF] Toh | Pseudozeros of polynomials and pseudospectra of companion matrices[END_REF] when i = n -k. The second step is to consider the products T k+1 (x)H n-k (x), for k = 2, 3, . . . , n -2, that is, the diagonal with black circular vertices in the grid. This step is very similar to the previous one, so we will only sketch the main ideas. We have to distinguish the cases k = 2, k = 3 and k > 3. When k = 2, using (1), ( 19) and ( 20), it may be proved

T 2 (x)H n-1 (x) = T 1 (x)H n-2 (x) -T 0 (x)H n-1 (x) -2A 0 T 1 (x) + • • • ,
where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. The previous equation shows that T 2 (x)H n-1 (x) may be computed from two products represented by triangular vertices in the grid: T 1 (x)H n-2 (x) and T 0 (x)H n-1 (x). Since we have seen that Theorem 4.1 holds for products represented by triangular vertices, it may be proved that (25) holds for T 2 (x)H n-1 (x).

Then, from (1) and ( 19), it may be proved that, when k = 2,

T 3 (x)H n-2 (x) = 2T 2 (x)H n-1 (x) + T 2 (x)H n-3 (x) -2A 1 T 2 (x),
and, when k > 3,

T k+1 (x)H n-k (x) = T k (x)H n-k+1 (x) + T k (x)H n-k-1 (x) -T k-1 (x)H n-k (x) -2A k-1 T k (x).
These two equations show that T k+1 (x)H n-k (x) may be computed from two products represented by triangular vertices, and the product represented by the black circular vertex above the black circular vertex corresponding to T k+1 (x)H n-k (x). Assuming that Theorem 4.1 holds for T 2 (x)H n-1 (x), the previous observation shows how to prove inductively (from top to bottom) that Theorem 4.1 holds for products corresponding to black circular vertices labeled with 2 and 3. Now, we address the products represented by circular vertices colored with different shades of grey, that is, the products T k+r-1 (x)H n-k (x), for r = 3, 4, . . . , n -2 and k = 1, 2, . . . , n -1 -r. We will show that Theorem 4.1 holds for products represented by vertices in the same diagonal (same shade of grey) assuming that it holds for products represented by (non-triangular) vertices in the diagonal on its left. Since we have previously proved that Theorem 4.1 holds for products represented by the white and black diagonals, this will imply that Theorem 4.1 holds for all products represented by grey vertices. For each grey diagonal, we have to distinguish the products represented by vertices labeled with 1, 2, and 3.

First, we consider the product T r (x)H n-1 (x), with r ≥ 3, represented by a grey vertex labeled with 1. From ( 1) and ( 20), we get T r (x)H n-1 (x) = (2xT r-1 (x) -T r-2 (x))H n-1 (x) = 2xT r-1 (x)H n-1 (x) -T r-2 (x)H n-1 (x) = T r-1 (x)H n-2 (x) -T r-2 (x)H n-1 (x) -2A 0 T r-1 (x) + • • • , where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. The previous equation shows that T r (x)H n-1 (x) may be computed from two products represented by vertices in the diagonal on its left: T r-1 (x)H n-2 (x) and T r-2 (x)H n-1 (x). Assuming that ( 24) and ( 25) hold for those products, we have where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. As can be checked, the previous equation corresponds to [START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF] with k = 1 and i = n -r.

The proof that Theorem 4.1 holds for products represented by grey vertices labeled with 2 is very similar to the previous one, so we omit it.

Finally, consider a product T n-i (x)H n-k (x) represented by a grey vertex labeled with 3. From ( 1) and ( 19), we have

T n-i (x)H n-k (x) =(2xT n-i-1 (x) -T n-i-2 (x))H n-k (x) =T n-i-1 (x)H n-k+1 (x) + T n-i-1 (x)H n-k-1 (x) -T n-i-2 (x)H n-k (x) -2A k-1 T n-i-1 (x)
The previous equation shows that T n-i (x)H n-k (x) may be computed from two products represented by (nontriangular) vertices in the diagonal on its left: T n-i-1 (x)H n-k-1 (x) and T n-i-2 (x)H n-k (x), and a product represented by a vertex in the same diagonal, above the vertex corresponding to T n-i (x)H n-k (x), that is, the product T n-i-1 (x)H n-k+1 (x). This observation shows how to prove inductively (from top to bottom) that Theorem 4.1 holds for the grey vertices labeled with 3 in the same diagonal. Assuming that ( 24) holds for T n-i-1 (x)H n-k-1 (x), T n-i-2 (x)H n-k (x) and T n-i-2 (x)H n-k (x), and using Γ i+1 -Γ i-1 = 2A n-i-1 , where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial, which shows that (24) holds also for T n-i (x)H n-k (x).

The final step of the proof consists in proving the uniqueness of r ik (x) and Q ik (x) in (39). For this purpose, assume that there exist two scalar polynomials r ik (x) and r ik (x), and two matrix polynomials Q ik (x) and Q ik (x) of degree at most n -1 such that T n-i (x)H n-k (x) = Q ik (x) + r ik (x)P (x) = Q ik (x) + r ik (x)P (x). Then, Q ik (x) -Q ik (x) = ( r ik (x) -r ik (x))P (x) is a matrix polynomial of degree at most n -1, but, if r ik (x) = r ik (x), the matrix polynomial ( r ik (x) -r ik (x))P (x) has degree larger than or equal to n, hence r ik (x) = r ik (x) and Q ik (x) = Q ik (x).

Theorem 4 . 4 .

 44 as the p × p matrix tr p (Z) := n i=1 Z ii . The block trace is used in Theorem 4.4, which provides a characterization of the subspace Sub T CT O(C T ), and is a generalization of [8, Theorem 4.1] when the matrix polynomial P (x) is expressed in the Chebyshev basis. For any Z ∈ R pn×pn , tr p (M k+1 Z) = 0, for k = 0, 1, . . . , n -1, (35) if and only if Z = C T X -XC T for some X ∈ R np×np with 0 first block column. (36)Moreover, either condition determines the first block row of Z uniquely given the remaining block rows.

T r- 1 AAAAAA

 1 (x)H n-2 (x) = n-r-1 ℓ=0 Γ ℓ T n-r+1-ℓ (x) + T |ℓ-n+r+1| (x) + Γ n-r T 1 (x)s T |r-ℓ-s| (x) + • • • , and T r-2 (x)H n-1 (x) = ℓ-s T |r-ℓ-s-1| (x) + • • • ,where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. Usingn-r-1 ℓ=0 Γ ℓ T n-r+1-ℓ (x) + T |ℓ-n+r+1| (x) + Γ n-r T 1 (x) -n-r ℓ=0 Γ ℓ T n-r+1-ℓ (x) -1 2 Γ n-r+1 T 0 (x) = n-r-2 ℓ=0 Γ ℓ T n-r-1-ℓ (x) + 1 2 Γ n-r-1 T 0 (x) -A r-1 T 0 (x),where we have used (Γ n-r+1 -Γ n-r-1 )/2 = A r-1 , andℓ-s T |r-ℓ-s-1| (x) = -ℓ+1-s T |r-ℓ-s| (x) -ℓ+1-s T |r-ℓ-s| (x) -r s=1 A r-s T |s-1| (x) + A r-1 T 0 (x) + A 0 T r-1 ℓ+1-s T |r-ℓ-s| (x) + A r-1 T 0 (x) + 2A 0 T r-1 ℓ-s T |r+1-ℓ-s| (x) + A r-1 T 0 (x) + 2A 0 T r-1 (x), we get T r-1 (x)H n-2 (x) = n-r-2 ℓ=0 Γ ℓ T n-r-1-ℓ (x) + 1 2 Γ n-r-1 T 0 (x) -r ℓ=1 ℓ s=1 2A ℓ-s T |r+1-ℓ-s| (x) + • • • ,

i- 1 ℓ=0Γ- 2 ℓ=0Γ- 2 ℓ=0Γ

 122 ℓ (T i+k-2-ℓ (x) + T |k+ℓ-i-2| (x)) + Γ i T k-2 (x) + i-1 ℓ=0 Γ ℓ (T i+k-ℓ (x) + T |k+ℓ-i| (x)) + Γ i T k (x) -i ℓ=0 Γ ℓ (T i+k-ℓ (x) + T |k+ℓ-i-2| (x)) -Γ i+1 T k-1 (x) = iℓ (T i+k-2-ℓ (x) + T |k+ℓ-i| (x)) + Γ i-1 T k-1 (x) -2A n-i-1 T k-1 , 2+ℓ-r T |n-i-ℓ-r| (x) -1+ℓ-r T |n-i-1-ℓ-r| (x) = -1+ℓ-r T |n-i+1-ℓ-r| (x) + 2A k-1 T n-i-1 (x) + 2A n-i-1 T k-1 (x). we get T n-i (x)H n-k (x) = T n-i-1 (x)H n-k+1 (x) + T n-i-1 (x)H n-k-1 (x) -T n-i-2 (x)H n-k (x) -2A k-1 T n-i-1 (x) = iℓ (T i+k-2-ℓ (x) + T |k+ℓ-i| (x)) + Γ i-1 T k-1 (x) -1+ℓ-r T |n-i+1-ℓ-r| (x) + • • • ,

  ijk A t with |α t,ijk | ≤ 4, where we set A n := I p .

	Proof. From Theorem 4.1 together with (26), we have

The vertices with triangular shape in the previous grid represent the cases in which the degree of T n-i (x)H n-k (x) does not exceed n -1, that is, when i ≥ n -k + 1. In this case, the polynomial Q ik (x) coincide with T n-i (x)H n-k (x), so we just need to expand T n-i (x)H n-k (x) in the Chebyshev basis. Indeed, when i = n and k = 1, from (42), we have

and when n - 40) and (41), we have

As can be checked, the two previous equations correspond to [START_REF] Trefethen | Approximation theory and approximation practice[END_REF] and [START_REF] Trefethen | Chebfun Version 5. The Chebfun Development Team[END_REF], respectively. Next, we consider the products T n-i (x)H n-k (x) with i < n -k + 1, represented in the grid by vertices with circular shape. This case is much more involved, since the degree of T n-i (x)H n-k (x) is larger than or equal to n. We will prove that (39) holds, with Q ik (x) as in ( 22)-( 25), each diagonal in the grid at a time (from left to right), showing that each product T n-i (x)H n-k (x) can be computed using, at most, a product represented by a vertex in the same diagonal and two products represented by vertices in the diagonal on its left.

The first step is to consider the products T k (x)H n-k (x), for k = 1, 2, . . . , n -1, that is, products represented by the diagonal with white circular vertices in the grid. We show that Theorem 4.1 holds for those products from top to bottom. We start with the white circular vertex labeled with 1 in the grid, that is, with the product T 1 (x)H n-1 (x). From ( 20) and (41), together with T 1 (x) = x, we have

where the dots correspond to something of the form r(x)P (x), with r(x) a scalar polynomial. As can be easily checked, the previous equation corresponds to [START_REF] Wilkinson | Rounding Errors in Algebraic Processes[END_REF] with i = n -1. Then, we consider the white circular vertex labeled with 2 in the grid, that is, the product T 2 (x)H n-2 (x).