A Provably Componentwise Backward Stable $O(n^2)$ QR Algorithm for the Diagonalization of Colleague Matrices

The roots of a monic polynomial expressed in a Chebyshev basis are known to be the eigenvalues of the so-called colleague matrix, which is a Hessenberg matrix that is the sum of a symmetric tridiagonal matrix and a rank-1 matrix. The rootfinding problem is thus reformulated as an eigenproblem, making the computation of the eigenvalues of such matrices a subject of significant practical importance. In this manuscript, we describe an $O(n^2)$ explicit structured QR algorithm for colleague matrices and prove that it is componentwise backward stable, in the sense that the backward error in the colleague matrix can be represented as relative perturbations to its components. A recent result of Noferini, Robol, and Vandebril shows that componentwise backward stability implies that the backward error $\delta c$ in the vector $c$ of Chebyshev expansion coefficients of the polynomial has the bound $\lVert \delta c \rVert \lesssim \lVert c \rVert u$, where $u$ is machine precision. Thus, the algorithm we describe has both the optimal backward error in the coefficients and the optimal cost $O(n^2)$. We illustrate the performance of the algorithm with several numerical examples.


Introduction
The problem of finding the roots of the polynomial is one of the oldest and most classical problems in mathematics. Countless methods have been proposed for its solution (see, for example, [28] for a history, and the two volumes [24] and [25] for a detailed account of such methods). In the 1800's, it was observed by Frobenius that the roots of the polynomial are the eigenvalues of a certain matrix called the companion matrix, formed using the polynomial coefficients. A matrix whose eigenvalues are the roots of p(x) is called a linearization of p(x). Given a linearization of p(x), the roots of the polynomial can thus be recovered by computing the eigenvalues of the matrix.
If the roots of the polynomial are found numerically, than the computed roots can be viewed as the exact roots of a perturbed polynomial p(x) + δp(x) with coefficients c i + δc i , where the size of the vector δc is called the backward error in the polynomial coefficients. The best possible bound on the backward error that a linearization method can have for general polynomials is δc c u, (2) or, in other words, that the relative normwise backward error in the polynomial coefficients is bounded by machine precision u (see, for example, [26]). The backward error in the companion matrix method was revealed by the influential paper [19], which analyzed the relationship between perturbations in the companion matrix and perturbations in the polynomial coefficients. There, the authors proved that if the companion matrix C is perturbed by the matrix E, then the matrix C + E is a linearization of the polynomial with perturbed coefficients c i + δc i , and that the perturbation satisfies the normwise bound δc c E u.
If the eigenvalues are computed by a standard QR algorithm, which is known to be backward stable (see, for example, [35]), then the computed eigenvalues are the exact eigenvalues of C + E, where E C u. Since C ≈ c , it follows that the backward error in the polynomial coefficients is bounded by Thus, as c get larger, the relative backward error in the coefficients increases. The companion matrix method, at first glance, would appear then to have two drawbacks: it falls short of the optimal backward error bound (2), and it costs O(n 3 ) operations as a result of using the QR algorthim. The situation improved dramatically in 2007, when Bini, Eidelman, Gemignani, and Gohberg published a paper [10] describing a stable, O(n 2 ) explicit QR method for companion matrices (around the same time, Chandrasekaran, Gu, Xia, and Zhu also discovered an O(n 2 ) method for companion matrices, see [15]). The algorithm is based on the observation that the companion matrix and its QR iterates have a certain structure which allows them to be represented by a collection of O(n) parameters called generators (specifically, the companion matrix is a Hessenberg matrix that is the sum of a unitary matrix and a rank-1 perturbation; matrices of this form are called fellow matrices). In 2010, an implicit version of this algorithm, also stable and O(n 2 ), and also based on generators, was introduced in [9]. Around the same time, Van Barel, Vandebril, Van Dooren, and Frederix discovered in [8] an alternative stable, O(n 2 ) implicit QR algorithm based on representing the unitary part by so-called core transformations, which are rotation matrices acting only on two adjacent rows at a time (see, for example, [2]). The first example of a proof of backward stability for an implicit O(n 2 ) QR algorithm for companion matrices was given by Aurentz, Mach, Vandebril, and Watkins in [4]; this algorithm is again based on core transformations. The backward stability result accompanying this QR algorithm guarantees the sub-optimal bound (2), but has the optimal complexity of O(n 2 ). Amazingly, the authors then discovered that the algorithm they had constructed, with some minor modifications, actually yields the optimal bound (4) in practice. An investigation showed that the reason for this remarkable behavior is that their algorithm is not just backward stable, but is componentwise backward stable, meaning that the backward error in the companion matrix can be decomposed into proportional backward errors in each of its components. They published a proof of the componentwise backward stability of their algorithm, together with a proof that componentwise backward stability guarantees the bound (2), in [3], along with numerical experiments.
Thus, if the coefficients of the polynomial p(x) in the monomial basis are known, then the algorithm of [3] is optimal in both error and time complexity. However, if the coefficients are not known, then the companion matrix cannot be used to the find the roots accurately, since the relationship between the values of the polynomial p(x) and the coefficients of its monomial expansion is highly unstable (this fact has been known for many decades, at least as early as Wilkinson [31]). If the polynomial p(x) is instead expanded in a basis of Chebyshev polynomials p(x) = c 0 + c 1 T 1 (x) + · · · + c n−1 T n−1 (x) + T n (x), (5) where T i (x) is the Chebyshev polynomial of order i, then the relationship between the coefficients and the polynomial is perfectly stable (see, for example, [36]). In fact, this observation is the basis for the Chebfun software package (see [6] and [16]). An analogue of the companion matrix, constructed from the Chebyshev expansion coefficients, was discovered in 1961 by Good [21], who called it the colleague matrix, and independently by Spect in 1957 [33]- [34]. The first O(n 2 ) algorithm for colleague matrices was discovered by Bini, Gemignani, and Pan in 2005 (even before [10] appeared) and is a stable, explicit QR algorithm based on generators [11]. Like the companion matrix, the colleague matrix has a special structure that is preserved over QR iterations (specifically, the colleague matrix is a Hessenberg matrix that is the sum of a Hermitian matrix and a rank-1 perturbation).
In 2008, Eidelman, Gemignani, and Gohberg, in [20], introduced a stable, O(n 2 ) implicit QR algorithm. The relationship between the backward error in the Chebyshev expansion coefficients and perturbations to the colleague matrix was first investigated Nakasukasa and Noferini in [26], where the authors found a lower bound for the backward error in the coefficients, showing that a backward stable QR algorithm can do no better than (4) (around the same time, Lawrence, Van Barel, and Van Dooren published a general analysis in [23], where they also proved a lower bound for colleague matrices). In [30], Perez and Noferini improved on this result and found an upper bound as well, proving that if the perturbation to the colleague matrix is small, then the bound (2) is achieved. The relationship between componentwise perturbations to the colleague matrix and the backward error in the coefficients was described completely in 2019 by Noferini, Robol, and Vandebril in [27].
Recently (see [32] and [14]), it was observed that certain O(n 2 ) structured QR algorithms for colleague matrices are surprisingly stable, attaining the bound (2) in many cases, an observation that mirrors the discovery in [3] for the case of companion matrices. However, unlike in [3], all previously proposed O(n 2 ) structured QR algorithms for colleague matrices have polynomials for which the worst-case bound (4) is attained. Thus, the question of whether or not there exists a structured O(n 2 ) QR algorithm that, when used to find the roots of a colleague matrix, attains the optimal bound (2), has remained open. In this manuscript, we answer this question in the affirmative by presenting a new, explicit O(n 2 ) QR algorithm for colleague matrices (in fact, for all Hessenberg matrices that have a Hermitian plus rank-1 structure), and prove that our algorithm is componentwise backward stable. Combined with the result in [27], this amounts to a proof that the optimal bound (2) is attained for all polynomials p(x). We demonstrate that this is indeed the case with several numerical experiments.
The structure of this manuscript is as follows. Section 2 describes the mathematical and numerical preliminaries. Section 3 describes the algorithm, and explains the significance of each step. In Section 4, we prove rigorously that the algorithm is componentwise backward stable. Section 5 presents the results of several numerical experiments. In Section 6, we discuss possible extensions and generalizations of the algorithm.

Preliminaries
In this section, we describe the mathematical and numerical preliminaries.

Linear Algebra
The following lemma states that if the sum of a Hermitian matrix and a rank-1 update pq * is lower Hessenberg, then the matrix is determined entirely by its diagonal and superdiagonal together with the vectors p and q. Lemma 2.1 (Eidelman, Gemignani, Gohberg [20]). Suppose that A ∈ C n×n is Hermitian, and let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that p, q ∈ C n and that A + pq * is lower Hessenberg. Then where a i,j denotes the (i, j)-th entry of A.
The following lemma states that if the sum of a matrix and a rank-1 update pq * is lower triangular, then the upper Hessenberg part of the matrix is determined entirely by its diagonal and subdiagonal, together with the vectors p and q. Lemma 2.2. Suppose that B ∈ C n×n and let d and γ denote the diagonal and subdiagonal of B, respectively. Suppose that p, q ∈ C n and that B + pq * is lower triangular. Then where b i,j denotes the (i, j)-th entry of B.
The following definition introduces two matrix seminorms that we will need in our error analysis. Definition 2.1. Suppose that A ∈ C n×n and let a i,j denote the (i, j)-th entry of A. We will use the notation · H to denote the square root of the sum of squares of the entries in the upper Hessenberg part of a matrix, so that Likewise, we will use the notation · T to denote the square root of the sum of squares of the entries in the upper triangular part, so that The following is a straightforward lemma stating that if a certain sequence of transformations is applied to a matrix on the right, then the upper triangular part of the result is determined by only the upper Hessenberg part of the original matrix. Lemma 2.3. Suppose that B ∈ C n×n , and let P 2 , P 3 , . . . , P n ∈ C n×n be matrices such that P k only affects the (k − 1, k)-plane of any vector it is applied to. Define P ∈ C n×n by the formula P = P 2 P 3 · · · P n . Then the upper triangular part of BP * is determined entirely by the upper Hessenberg part of B. Furthermore, if P 2 , P 3 , . . . , P n are unitary, then BP * T ≤ B H .

Error Analysis
The following definition introduces the notation used in the error analysis that appears in this manuscript. We follow the notation used in [22] and [3].
Definition 2.2. Evaluation of an expression in floating point arithmetic is denoted by fl (·), and we denote the unit roundoff (or machine epsilon) by u. We assume that where op stands for any of the basic arithmetic operations +, −, * , /. We denote computed quantities by a hat, so that x denotes the computed approximation to x. We use the notation to mean "less than or equal to the right hand side times a modest multiplicative constant depending on n as a low-degree polynomial", where the meaning of n is clear from the context. Whenever a matrix or vector norm appears to the left or right of , we omit the particular choice of norm, since in finite dimensions all norms are equivalent. When u appears in an expression on the right hand side of , we ignore all higher order powers of u.

Floating point computation of complex plane rotations
The following lemma bounds the forward error of the floating point computation of a complex plane rotation (see, for example, §20 of [37]).

Multiplication by complex plane rotations
The following lemma estimates the forward error of applying a plane rotation to a vector (see, for example, §21 of [37]).
Lemma 2.5. Suppose that Q ∈ SU(2) is a complex rotation matrix, and let Q be a floating point approximation to Q satisfying (14). Suppose further that

Colleague Matrices and Polynomial Rootfinding
Suppose that p(x) is a monic polynomial of order n represented by where c j ∈ R, c n = 1, and T j (x) is the Chebyshev polynomial of order j. It turns out that the roots of p(x) are the eigenvalues of the (scaled) n × n colleague matrix where e n is the n-th unit vector (see, for example, [21]). A matrix C whose eigenvalues are the roots of p(x) is called a linearization of p(x). Letting and we see that the colleague matrix C can be written as where A is Hermitian and C is lower Hessenberg. The following beautiful theorem by Noferini, Robol, and Vandebril (see Corollary 5.4 of [27]) bounds the change in the coefficients of the polynomial being linearized by the componentwise perturbations of the colleague matrix C = A + e n q * . Theorem 2.6. Let C = A + e n q * be the linearization (17) of the monic polynomial p(x), expressed in the Chebyshev basis. Consider the perturbations δA ≤ A , δe n ≤ n , and δq ≤ q . Then, the matrix is a linearization of the polynomial where δc n + q + c A .

Stability of Rootfinding Using Linearizations
Suppose that p(x) is a monic polynomial of order n represented in the Chebshev basis where c j ∈ R, c n = 1, and T j (x) is the Chebyshev polynomial of order j, and let the roots of p(x) be denoted by x 1 , x 2 , . . . x n ∈ C. Suppose that a rootfinding algorithm returns the computed roots x 1 , x 2 , . . . , x n ∈ C. If the computed roots are the exact roots of some perturbed polynomial where δc c u, then we say that the rootfinding algorithm is backward stable. In fact, this is the best backward stability bound that can be hoped for, for general polynomials p(x) (see the discussion in Appendix A of [26]).
Remark 2.1. Suppose that p(x) is a polynomial of order n, that is not monic, expressed in the Chebyshev basis where a j ∈ R and T j (x) is the Chebyshev polynomial of order j. Clearly, the roots of p(x) are identical to the roots of p(x)/a n . Let c j = a j /a n , for j = 0, 1, . . . , n. If a backward stable rootfinding algorithm is applied to the monic polynomial p(x)/a n , then, letting δa = δc · a n , the algorithm is also backward stable with respect to the original coefficients a j , since When linearization is used as a rootfinding algorithm, the stability of the computed roots comes from the stability of the eigenvalue algorithm applied to the colleague matrix C. If the eigenvalues of C are computed by an unstructured QR algorithm, then the backward error δC on C is bounded by δC C u. Since the backward error is unstructured, it follows that δA ≈ C u, so, by Theorem 2.6 together with the fact that C ≈ c , the backward error in the coefficients is bounded by δc c 2 u.
Remark 2.2. This backward error can be reduced by partially balancing the matrix C.
Suppose that, instead of computing the eigenvalues of C, we compute the eigenvalues of

Remark 2.3.
Let bal(C) denote the matrix C after complete balancing. Remarkably, in some situations, bal(C) ≈ 1 even when c is large. Thus, complete balancing can completely eliminate large entries in C, at the expense of destroying its symmetric tridiagonal plus rank-1 structure. See Remark 5.1, as well as the paper [29], for a more detailed discussion.
Remark 2.4. While unstructured QR, applied to the colleague matrix, is known to achieve only the backward error bound δc c 2 u, the QZ algorithm, applied to an appropriately scaled matrix pencil, does result in a backward stable rootfinder with the bound δc c u (see, for example, [26]). This is because the eigenvalue problem for the colleague matrix can be written as a matrix pencil A − λB, where both A and B are small, and a backward stable QZ algorithm applied to the pencil computes the exact eigenvalues of a perturbed pencil (A + δA) − λ(B + δB), where δA A u and δB B u. Unfortunately, it appears to be very difficult to construct structured, O(n 2 ) QZ algorithms for colleague matrices that retain the nice stability properties of the unstructured, O(n 3 ) QZ algorithm.
The following theorem, stated in a slightly different form in [27], says that if the eigenvalues of the colleague matrix C = A + e n q * are computed using a componentwise backward stable algorithm, then linearization is backward stable as a rootfinding algorithm. It follows immediately from Theorem 2.6.
Theorem 2.7. Suppose that the eigenvalues of the colleague matrix C = A + e n q * are computed by a componentwise backward stable algorithm, in the sense that the computed eigenvalues are the exact eigenvalues of the matrix where δA A u ≈ u, δe n e n u ≈ u, and δq q u. Then, linearization is backward stable as a rootfinding algorithm, with δc c u.

Conventions
It was pointed out to the authors that, while the Hessenberg matrices in this manuscript are all lower Hessenberg, the standard convention in numerical linear algebra is to study the transpose of the problem, and consider only upper Hessenberg matrices (see [18]). The upper Hessenberg form is much better notationally since, in upper Hessenberg form, the first elimination step eliminates the entry in the (2, 1)-position, while, in lower Hessenberg form, the entry in the (n − 1, n)-position is eliminated first. Furthermore, the upper Hessenberg form is more convenient when representing polynomials in the Lagrange basis (see, for example, [17]). Unfortunately, at the time that this was all pointed out, most of the writing and numerical codes were complete, and had been written in lower Hessenberg form because of a historical fluke related to the structure of old explicit QR codes that were used as a template for our algorithm.

The Algorithm
In this section, we give an overview of our algorithm. We begin by describing the class of matrices our algorithm can be applied to. Let F n ⊂ C n×n be the set of lower Hessenberg matrices of the form where A ∈ C n×n is Hermitian and p, q * ∈ C n . Eidelman, Gemignani, and Gohberg observed in [20] that the matrix A is determined entirely by: 1. The diagonal entries d i = a i,i , for i = 1, 2, . . . , n; 2. The superdiagonal entries β i = a i,i+1 for i = 1, 2, . . . , n − 1; 3. The vectors p and q (see Lemma 2.1). Following [20], we call these four vectors the basic elements or generators of A. In [20], the authors construct an implicit QR algorithm that takes advantage of this structure to achieve a cost of O(n 2 ). They also prove that their algorithm is backward stable, in the sense that, when the algorithm is used to compute the eigenvalues of a matrix C ∈ F n , the computed eigenvalues are the exact eigenvalues of C + δC, where δC C u. This is the same backward stability bound that is provided by an unstructured QR algorithm.
In this manuscript, we describe a new explicit QR algorithm for matrices A+pq * ∈ F n , that also has the cost O(n 2 ), and prove that our algorithm is componentwise backward stable, in the sense that the computed eigenvalues are the exact eigenvalues of (A + δA) + (p + δp)(q + δq) * , where δA A u, δp p u, and δq q u. To motivate our algorithm, consider first the naive unshifted QR algorithm in exact arithmetic, applied to the matrix C = A + pq * . Let the matrix U n ∈ C n×n be the unitary matrix that rotates the (n − 1, n)-plane so that eliminating the superdiagonal in the (n − 1, n)-th position. Likewise, let U n−1 ∈ C n×n denote the unitary matrix rotating the (n − 2, n − 1)-plane so that eliminating the superdiagonal in the (n − 2, n − 1)-th position. Continuing in this fashion, let U n−2 , U n−3 , . . . , U 2 be the unitary matrices eliminating the superdiagonal entries in the (n−3, n−2), (n−4, n−3), . . . , (1, 2) positions of the matrices (U n−1 U n C), (U n−2 U n−1 U n C), . . . , (U 3 · · · U n−1 U n C), respectively. Letting U = U 2 U 3 · · · U n , we have that U C is lower triangular. This matrix has the form where the upper Hessenberg part of the matrix B = U A is determined entirely by: 1. The diagonal entries d i = b i,i , for i = 1, 2, . . . , n; 2. The subdiagonal entries γ i = b i+1,i , for i = 1, 2, . . . , n − 1; 3. The vectors p = U p and q (see Lemma 2.2). Like with the matrix A, we call these four vectors the basic elements or generators of (the upper Hessenberg part of) B. Next, the matrix is multiplied by U * on the right; clearly, so It's easy to show that, since U C is lower triangular and U * = U * n U * n−1 · · · U * 2 , where U k rotates the (k − 1, k)-plane, the matrix U CU * is lower Hessenberg. Thus, U CU * ∈ F n , and the matrix A = U AU * is determined entirely by its diagonal and superdiagonal, together with p = U p and q = U q. Furthermore, the upper triangular part of U AU * is determined entirely by the upper Hessenberg part of B (see Lemma 2.3), and since U AU * is Hermitian, it follows that the whole of the matrix U AU * is determined entirely by the upper Hessenberg part of B.
In our algorithm, we use only the basic elements of A and B to represent our matrices. This results in a single iteration of our QR algorithm requiring O(n) operations. Furthermore, we prove that the matrix A and vectors p and q, computed by a single iteration of our QR algorithm, have the componentwise forward error bounds A − A A u, p − p p u, and q − q p u. We then show that these componentwise forward error bounds result in componentwise backward stability.
To eliminate the superdiagonal element in the (k − 1, k) position of B (k+1) + p (k+1) q * , we first compute the rotation matrix Q k ∈ SU(2) that eliminates it by a rotation in the (k − 1, k)-plane (see Line 4 of Algorithm 1). Next, we apply the rotation matrix separately to the generators of B (k+1) and to the vector p (k+1) . Since we are only interested in computing the upper Hessenberg part of B (k) , we need to update the subdiagonal element in the (k − 1, k − 2) position of B (k+1) , represented by γ (k+1) k−2 (see Figure 1). However, this calculation requires the sub-subdiagonal entry in the (k, k − 2) position of B (k+1) , which is unknown to us since only the upper Hessenberg part of B (k+1) is available. Fortunately, it can be recovered by the following trick. Since is also Hermitian. Thus, Furthermore, since right-multiplication by U * j only affects columns j and j − 1 (see Figure 1), we have that right-multiplication by We know that the entries in the (k − 2, k), (k − 2, k + 1), . . . , (k − 2, n) positions of B (k+1) are inferred from p (k+1) and q by the formula for = k, k + 1, . . . , n. Combining (37) and (38), we thus have that the sub-subdiagonal entry in the (k, k − 2) position of B (k+1) can be recovered by the formula Defining By computing the vector q , respectively, are updated in Line 9, and the vector p (k+1) is rotated in Line 10.
Since we've eliminated the superdiagonal element in the (k − 1, k) position of B (k+1) + p (k+1) q * , we have that the (k − 1, k) element of the matrix B (k) is inferred from p (k) and q by the formula Now, we would like the upper Hessenberg part of B (k) to have a small componentwise error, so that However, consider the following scenario. Suppose that the norm of ( p u, which will be much larger than A u. In this situation then, even if To remedy this, we must apply a correction to p By Lemma 2.5, we have that and (see Line 12 of Algorithm 1). With this correction to p , then the correction is not neccessary, since in this case the error in p k−1 . This process of eliminating the superdiagonal elements can be repeated, until the upper Hessenberg part of the matrix B, approximating B = U 2 U 3 · · · U n A, is obtained, together with p, approximating p = U 2 U 3 · · · U n p (see Algorithm 1). In Section 4.1, Lemma 4.1, we prove that the forward errors in the upper Hessenberg part of B and in the vector p are proportional to A u and p u, respectively. Algorithm 1 (A single elimination of the superdiagonal) Inputs: This algorithm accepts as inputs two vectors d and β representing the diagonal and superdiagonal, respectively, of an n × n Hermitian matrix A, as well as two vectors p and q of length n, where A + pq * is lower Hessenberg. Outputs: It returns as its outputs the rotation matrices Q 2 , Q 3 , . . . , Q n ∈ C 2×2 so that, letting U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k −1, k)-plane by Q k , U 2 U 3 · · · U n (A+pq * ) is lower triangular. It also returns the vectors d, γ, and p, where d and γ represent the diagonal and subdiagonal, respectively, of the matrix U 2 U 3 · · · U n A, and p = U 2 U 3 · · · U n p.
2: Make a copy of q, setting q ← q.
Rotate the subdiagonal and the sub-subdiagonal: Rotate the diagonal and the subdiagonal:

9:
Rotate the superdiagonal and the diagonal: 10: Rotate p: Correct the vector p, Rotate q:

Rotating Back to Hessenberg Form
In this section, we describe how our algorithm rotates the triangular matrix produced by an elimination of the superdiagonal back to lower Hessenberg form (see Algorithm 2). Suppose that B is an n × n matrix and that p and q are vectors such that B + pq * is lower triangular. Notice that this condition is satisfied by the matrix B and the vectors p and q from the preceding section, produced by an elimination of the superdiagonal. Let γ denote the subdiagonal of B. Suppose that we have already applied the rotation matrices U * n , U * n−1 , · · · , U * k+1 to the right of B and q * , and let q (k+1) = U k+1 U k+2 · · · U n q and A (k+1) = BU * n U * n−1 · · · U * k+1 . Suppose that q (k+1) and A (k+1) are the computed approximations to q (k+1) and A (k+1) , respectively, and that the upper triangular part of A (k+1) is represented by its generators: 3. The vectors p and q (k+1) , from which the remaining elements in the upper triangular part are inferred.
To apply the matrix U * k to A (k+1) + p q (k+1) on the right, we apply the matrix Q * k ∈ SU(2) separately to the generators of A (k+1) and to the vector q (k+1) . We start by rotating the diagonal and superdiagonal elements in the (k − 1, k − 1) and , respectively, in Line 2 of Algorithm 2, saving the superdiagonal element in β (k) k−1 (see Figure 2). Next, we rotate the elements in the (k, k − 1) and (k, k) positions, represented by γ k−1 (the (k − 1)-st element of the subdiagonal of B) and d (k+1) k , respectively, in a straightforward way in Line 3; since we are only interested in computing the upper triangular part of A (k) , we only update the diagonal entry. Finally, we rotate the vector q (k+1) in Line 4.
The process of applying the rotation matrices on the right can be repeated, until the upper triangular part of the matrix A, approximating A = BU * n U * n−1 · · · U * 2 , is obtained, together with q, approximating q = U 2 U 3 · · · U n q (see Algorithm 2). In Section 4.1, Lemma 4.2, we prove that the forward errors in the upper triangular part of A and in the vector q are proportional to B H u and q u, respectively.

The QR Algorithms
The elimination of the superdiagonal described in Algorithm 1, followed by the rotation back to Hessenberg form described in Algorithm 2, can be iterated to find the eigenvalues of A + pq * . Our unshifted explicit QR algorithm, based on this iteration, is described in Algorithm 3. This unshifted QR algorithm can be accelerated by the introduction of shifts; our explicit shifted QR algorithm, with Wilkinson shifts, is described in Algorithm 4.
In Section 4.1, we show that the forward error of one iteration of our QR algorithm (Algorithm 1 followed by Algorithm 2) satisfies componentwise forward error bounds. In Section 4.2, we use this result to prove that both our explicit unshifted QR algorithm (Algorithm 3) and our shifted QR algorithm (Algorithm 4) are componentwise backward stable.
Algorithm 2 (Rotating the matrix back to Hessenberg form) Inputs: This algorithm accepts as inputs n − 1 rotation matrices Q 2 , Q 2 , . . . , Q n ∈ C n×n , two vectors d and γ representing the diagonal and subdiagonal, respectively, of an n × n complex matrix B, and two vectors p and q of length n, where B + pq * is lower triangular. Outputs: Letting U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k − 1, k)-plane by Q k , this algorithm returns as its outputs the vectors d, β, and q, where d and β represent the diagonal and superdiagonal, respectively, of the matrix BU * n U * n−1 · · · U * 2 , and q = U 2 U 3 · · · U n q. Rotate the diagonal and the superdiagonal:

3:
Rotate the subdiagonal and the diagonal: Rotate q: Algorithm 3 (Unshifted explicit QR) Inputs: This algorithm accepts as inputs two vectors d and β representing the diagonal and superdiagonal, respectively, of an n × n Hermitian matrix A, as well as two vectors p and q of length n, where A + pq * is lower Hessenberg. It also accepts a tolerance > 0, which determines the accuracy the eigenvalues are computed to. Outputs: It returns as its output the vector λ of length n containing the eigenvalues of the matrix A + pq * .
Perform one iteration of QR (one step of Algorithm 1 followed by one step of Algorithm 2) on the submatrix (A + pq * ) i:n,i:n defined by the vectors d i:n , β i:n−1 , p i:n , and q i:n . Algorithm 4 (Shifted explicit QR) Inputs: This algorithm accepts as inputs two vectors d and β representing the diagonal and superdiagonal, respectively, of an n × n Hermitian matrix A, as well as two vectors p and q of length n, where A + pq * is lower Hessenberg. It also accepts a tolerance > 0, which determines the accuracy the eigenvalues are computed to. Outputs: It returns as its output the vector λ of length n containing the eigenvalues of the matrix A + pq * .

5:
Set µ to whichever of µ 1 and µ 2 is closest to d i + p i q * i .

8:
Perform one iteration of QR (one step of Algorithm 1 followed by one step of Algorithm 2) on the submatrix (A + pq * ) i:n,i:n defined by the vectors d i:n , β i:n−1 , p i:n , and q i:n . 9: end while 10:

Componentwise Backward Stability
The principal results of this section are Theorems 4.6 and 4.7, which state that our unshifted and shifted QR algorithms, respectively, are componentwise backward stable. In Section 4.1, we prove that the forward error of a single sweep of our unshifted QR algorithm satisfies componentwise bounds. In Section 4.2, we use these bounds show componentwise backward stability of our QR algorithms.

Forward Error Analysis of a Single Sweep of QR
Suppose that A is Hermitian and A + pq * is lower Hessenberg. In this section, we prove in Theorem 4.3 that the forward errors in A, p, and q of single sweep of our explicit QR algorithm are proportional to A u, p u, and q u, respectively. The following lemma bounds the forward error of Algorithm 1 (the elimination of the superdiagonal).
Lemma 4.1. Suppose that A ∈ C n×n is a Hermitian matrix, and that p, q ∈ C n .
Suppose further that A + pq * is lower Hessenberg, and let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that Algorithm 1 is carried out in floating point arithmetic with d, β, p, and q as inputs, and let Q 2 , Q 3 , . . . , Q n ∈ SU(2) be the unitary matrices generated by an exact step of Line 4 of Algorithm 1 applied to the computed vectors at that step. Let U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k − 1, k)-plane by Q k , and define U ∈ C n×n by the formula U = U 2 U 3 · · · U n . Suppose finally that d, γ, and p are the outputs generated by Algorithm 1, and define the upper Hessenberg part of the matrix B ∈ C n×n by the formula where b i,j denotes the (i, j)-th entry of B. Let B = U A and p = U p. Then and where · H denotes the square root of the sum of squares of the entries in the upper Hessenberg part of its argument (see Definition 2.1).
denote the computed vectors in Algorithm 1 after the elimination of the superdiagonal elements in the positions (n − 1, n), (n − 2, n − 1), . . . , (k − 1, k). Suppose further that the upper Hessenberg part of the matrix B (k) ∈ C n×n is defined by the formula where b (k) i,j denotes the (i, j)-th entry of B (k) . Clearly, d = d (2) , γ = γ (2) , p = p (2) , and B = B (2) . Let B (k) = U k U k+1 · · · U n A and p (k) = U k U k+1 · · · U n p. We will prove that p u, for each k = n, n − 1, . . . , 2. We begin by proving this statement for k = n. From Line 4, we have that the matrix Q n ∈ SU(2) satisfies with the computed matrix Q n satisfying Q n − Q n u by Lemma 2.4. In Line 6, we have At this stage q is still equal to q and, according to Lemma 2.1, a n,n−2 = −q n p * n−2 . By definition, a n−1,n−2 = γ n−2 . Therefore, by Lemma 2.5, we have that | γ Since a n−1,n−1 = d n−1 and a n,n−1 = γ n−1 , by Lemma 2.5, we have that | d Since, by definition, a n−1,n = β n−1 and a n,n = d n , it follows from Lemma 2.5 that n,n | |β n−1 | 2 + |d n | 2 u ≤ A u. In Line 10, we have where p (n) † n−1 is a temporary value that will be corrected later. Once again, Lemma 2.5 tells us that | p In Line 12, we apply a correction to p (n) † n−1 , so that From this, we see that, if |p n−1 q * n | 2 + |p n q * n | 2 > |β n−1 | 2 + |d n | 2 , then where the second equality is due to (56). Furthermore, If, on the other hand, |p n−1 q * n | 2 + |p n q * n | 2 ≤ |β n−1 | 2 + |d n | 2 , then Moreover, where the second equality follows from (56). This completes the proof that B (n) − B (n) with the computed matrix Q k satisfying Q k − Q k u by Lemma 2.4. In Line 6, we have We must first show that We begin by observing that since right-multiplication by U * j only affects columns j and j − 1. We now observe that is Hermitian, so from (65) we have that By the induction hypothesis, and for all = k, k + 1, . . . , n. Thus, Combining (70) with (67), where q (k+1) = U k+1 U k+2 · · · U n q. From Line 14 we have and, by repeated application of Lemma 2.5, Thus, where the first equality follows by (68) and the second inequality by (69). Combining (71) and (74), or, equivalently, Finally, since, by the induction hypothesis, we use (76) and (77) and apply Lemma 2.5 to (63) to find that | γ By the induction hypothesis, | d Thus, another application of Lemma 2.5 shows that | d By the induction hypothesis, | β (k+1) In Line 10, we then have where p (k) † k−1 is a temporary value that will be corrected later. By the induction hypothesis, and definep Clearly,β by the definition of Q k (see (62)). Also, by Lemma 2.5, we have that |β In Line 12, we apply a correction to p We then observe that and −β Finally, combining (85), (86), and (87), we find that | p If, conversely, | p Next, we observe that Since −p and −p we combine (90), (91), and (92) to see that |− p
The following lemma bounds the forward error of Algorithm 2 (the rotation back to Hessenberg form).
Lemma 4.2. Suppose that B ∈ C n×n and p, q ∈ C n . Suppose further that B + pq * is lower triangular, and let d and γ denote the diagonal and subdiagonal of B, respectively. Suppose that Q 2 , Q 3 , . . . , Q n ∈ SU(2), and suppose that Algorithm 2 is carried out in floating point arithmetic, using d, γ, p, q, and Q 2 , Q 3 , . . . , Q n as inputs. Suppose finally that d, β, and q are the outputs generated by Algorithm 2, and define the upper triangular part of the matrix A ∈ C n×n by the formula where a i,j denotes the (i, j)-th entry of A. Let U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k − 1, k)-plane by Q k . Define U ∈ C n×n by the formula U = U 2 U 3 · · · U n , and let A = BU * and q = U q. Then and where · T denotes the square root of the sum of squares of the entries in the upper triangular part of its argument and · H denotes the square root of the sum of squares of the upper Hessenberg part (see Definition 2.1).

an application of Lemma 2.5 gives us
We first observe that since right multiplication by U * j only affects columns j and j − 1. Since, by the induction hypothesis, | d By the induction hypothesis, | q All that's left now is to prove that |−p q ,k | B H u, for all for all = 1, 2, . . . , k − 2. By the induction hypothesis, and for all = 1, 2, . . . , k − 2. Defineq Multiplying (112) by p * , we have which, combined with (110) and (111) and the fact that A (k+1) and −p q for all = 1, 2, . . . , k − 2. By Lemma 2.5, from which it follows that and for all = 1, 2, . . . , k − 2, where the second inequality follows from (110) and (111) ,k | B H u, for all for all = 1, 2, . . . , k − 2, and we are done.
The following theorem bounds the forward errors of a full sweep of our QR algorithm, and is the principal result of this subsection. Theorem 4.3. Suppose that A ∈ C n×n is a Hermitian matrix, that p, q ∈ C n , and that A + pq * is lower Hessenberg. Let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that Algorithm 1 is carried out in floating point arithmetic, and let Q 2 , Q 3 , . . . , Q n ∈ SU(2) be the unitary matrices generated by an exact step of Line 4 of Algorithm 1 applied to the computed vectors at that step. Let U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k − 1, k)-plane by Q k , and define U ∈ C n×n by the formula U = U 2 U 3 · · · U n . Suppose that Algorithm 2 is then carried out in floating point arithmetic, using the outputs of Algorithm 1 as inputs. Suppose finally that p is an output of Algorithm 1 and q, d, and β are all outputs of Algorithm 2, and define the matrix A by the formula where a i,j denotes the (i, j)-th entry of A. Let A = U AU * , p = U p, and q = U q. Then and q − q q u. (123) Proof. Suppose that B (defined by (47)), p, and Q 2 , Q 3 , . . . , Q n ∈ C n×n are outputs of Algorithm 1. Let B = U A and p = U p. By Lemma 4.1, and where · H denotes the square root of the sum of squares of the entries in the upper Hessenberg part of its argument (see Definition 2.1). Now suppose that B, p, q, and Q 2 , Q 3 , . . . , Q n ∈ C n×n are used as inputs to Algorithm 2. Let U k ∈ C n×n , k = 2, 3, . . . , n, denote the matrices that rotate the (k − 1, k)-plane by Q k , and define U ∈ C n×n by the formula U = U 2 U 3 · · · U n . Let A = B U * and q = U q, and observe that the upper triangular part of A is well-defined due to Lemma 2.3. By Lemma 4.2 we have that and where · T denotes the square root of the sum of squares of the entries in the upper triangular part of its argument and · H denotes the square root of the sum of squares of the upper Hessenberg part (see Definition 2.1). Let A = BU * = U AU * and let q = U q. We observe that where the last inequality follows from (124) and the fact that U − U u. Since, clearly, B H u A u, we combine (126) and (128) to get Now we observe that, since both A = U AU * and A are Hermitian, Next, we observe that so, combining (127) and (131), we have and we are done.

Backward Error Analysis of the QR Algorithms
Suppose that A is Hermitian and A + pq * is lower Hessenberg. In this section, we prove in Theorems 4.6 and 4.7 that the backward errors in A, p, and q of both our explicit unshifted QR algorithm (see Algorithm 3) and explicit shifted QR algorithm (see Algorithm 4) are proportional to A u, p u, and q u, respectively. The following lemma states that a single iteration of our QR algorithm is componentwise backward stable.
Lemma 4.4. Suppose that A ∈ C n×n is a Hermitian matrix, that p, q ∈ C n , and that A + pq * is lower Hessenberg. Let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that a single iteration of our QR algorithm (Algorithm 1 followed by Algorithm 2) is carried out in floating point arithmetic, and let p, q, d, and β denote the outputs of the algorithm. Define the matrix A by the formula (120). Then there exists a unitary matrix U ∈ C n×n , a matrix δA ∈ C n×n , and vectors δp, δq ∈ C n , such that and where δA A u, δp p u, and δq q u.
Proof. Let U ∈ C n×n be the unitary matrix defined in the statement of Theorem 4.3, and let A = U AU * , p = U p, and q = U q. By Theorem 4.3, where δA A u. Thus, where δA = U * δAU . Since U is unitary, clearly δA A u. Likewise, by Theorem 4.3, where δp p u, so where δp = U * δp and δp p u. Similarly, where δq q u.
The following lemma states that repeated iterations of our QR algorithm are componentwise backward stable.
Lemma 4.5. Suppose that A ∈ C n×n is a Hermitian matrix, that p, q ∈ C n , and that A + pq * is lower Hessenberg. Let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that k iterations of our QR algorithm (Algorithm 1 followed by Algorithm 2) are carried out in floating point arithmetic, and let p (k) , q (k) , d (k) , and β (k) denote the outputs of the algorithm. Define the matrix A (k) by the formula (120), making the obvious substitutions. Then there exists a unitary matrix U ∈ C n×n , a matrix δA ∈ C n×n , and vectors δp, δq ∈ C n , such that and where δA A u, δp p u, and δq q u.
The following theorem states that our explicit unshifted QR algorithm is componentwise backward stable. Theorem 4.6 (Explicit unshifted QR). Suppose that A ∈ C n×n is a Hermitian matrix, that p, q ∈ C n , and that A + pq * is lower Hessenberg. Let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that Algorithm 3 is carried out in floating point arithmetic with A u, and let λ 1 , λ 2 , . . . , λ n denote the outputs. Then there exist a matrix δA ∈ C n×n and vectors δp, δq ∈ C n such that λ 1 , λ 2 , . . . , λ n are the exact eigenvalues of the matrix where δA A u, δp p u, and δq q u.
The following theorem states that our explicit shifted QR algorithm is componentwise backward stable, for those eigenvalues for which the shifts are small. Theorem 4.7 (Explicit shifted QR). Suppose that A ∈ C n×n is a Hermitian matrix, that p, q ∈ C n , and that A + pq * is lower Hessenberg. Let d and β denote the diagonal and superdiagonal of A, respectively. Suppose that Algorithm 4 is carried out in floating point arithmetic with A u, and suppose that µ ( ) is the largest total shift encountered at any point during the course of the algorithm from i = 1, 2, . . . , in the outer loop. Let λ 1 , λ 2 , . . . , λ n denote the outputs of the algorithm. Then, for each = 1, 2, . . . , n, there exist a matrix δA ∈ C n×n and vectors δp, δq ∈ C n such that λ 1 , λ 2 , . . . , λ are exact eigenvalues of the matrix where δA ( A + |µ ( ) |)u, δp p u, and δq q u.
Proof. The proof is essentially identical to the proof of Theorem 4.6, and we omit it.
Remark 4.1. Notice that Theorems 4.6 and 4.7 do not make any mention of convergence. What they say is that, if the algorithm converges, then it is componentwise backward stable. We observe that, in practice, Algorithm 4 always converges rapidly, at least quadratically, for ≈ A u.

Remark 4.2.
Notice that the bound on δA in Theorem 4.7 involves µ ( ) , which is the largest total shift encountered at any point during the calculation of λ 1 , λ 2 , . . . , λ . While this bound appears weaker than the corresponding bound in Theorem 4.6, in practice it turns out to be essentially the same, as follows. We can always assume that A is much smaller than p, or q, or both; if this isn't the case, then componentwise stability no longer has any special meaning, since it follows immediately from the usual Bauer-Fike perturbation bounds (see [7]). Furthermore, we tend to be interested in the componentwise backward stability of small eigenvalues λ i , where | λ i | ≈ A . If we perform a few iterations of unshifted QR on the matrix, then the eigenvalues of the top-left 2×2 block will approach the two smallest eigenvalues of the matrix (recalling that our QR algorithm works with lower Hessenberg matrices). If we now use Algorithm 4, we'll find that the initial shift is small and, as a result, all the eigenvalues are computed roughly in order from smallest to largest. This means that |µ (i) | ≈ | λ i | and (approximately) λ 1 < λ 2 < · · · < λ i . For λ i such that | λ i | ≈ A , we have then that the bound δA ( A + |µ (i) |)u becomes δA A u. Finally, we point out that the dependence of the bound on µ ( ) could likely be removed entirely by reformulating our QR algorithm as an implicit method.

Numerical Results
In this section, we demonstrate the componentwise backward stability of our shifted QR algorithm (see Algorithm 4) by illustrating its stability when it is used as a rootfinding algorithm (see Sections 2.3 and 2.4). Consider a polynomial p(x) of order n, not necessarily monic, expressed in a Chebyshev polynomial basis where a j ∈ R and T j (x) is the Chebyshev polynomial of order j. By Theorem 2.7 and Remark 2.1, we have that if the eigenvalues of the linearization (17), where c j = a j /a n , j = 0, 1, . . . , n, are computed by a componentwise backward stable algorithm, then the computed roots x 1 , x 2 , . . . , x n are the exact roots of the perturbed polynomial where δa a u.
By applying our QR algorithm to linearizations of various polynomials p(x), we demonstrate our algorithm's componentwise backward stability by showing that the bound (179) always holds. We estimate the size of the backward error δa in the coefficients by using the following observation (see the discussion accompanying Table 1 in [26]). By the definition of p(x) + δp(x), we have that (p + δp)( x i ) = 0 for i = 1, 2, . . . , n. From (177) and (178), it follows that Since −1 ≤ T j (x) ≤ 1 for all j when x ∈ [−1, 1], we have whenever x i ∈ C is not too far from the interval [−1, 1].
Even though x i is already a floating point number, the value p( x i ) cannot be computed exactly in floating point arithmetic. Letting p( x i ) denote the approximation to p( x i ) computed in floating point arithmetic, we know that where is the absolute condition number of p(x) at x = x i . When κ(p; x i ) is large, the error in evaluating p( x i ) dominates, while when κ(p; x i ) is of modest size, we have p( x i ) ≈ p( x i ).
In this section, we investigate the quantity for various polynomials p(x). When κ(p; x i ) ≥ a , we have When κ(p; x i ) ≤ a , Thus, if our QR algorithm is indeed componentwise backward stable and (179) is satisfied, then, by (181), we expect to find that for all polynomials p(x).
For p( x i ) to be a good approximation to δa (see (181)), we stated that x i ∈ C should be "not too far from the interval [−1, 1]." We make this notion precise as follows. Let z i , i = 1, 2, . . . , n denote the exact roots of the order-n polynomial p(x), and let z i denote the computed roots. We select roots close to the interval [−1, 1] by choosing some δ > 0 (for example, δ = 10 −3 ), and letting x i ∈ R denote the real part of all roots z i that are inside the rectangle If the polynomial p(x) has the real root z i , then taking the real part of z i will not result in any additional error. The number of real roots inside the region (188) will often be less than the order n, and we denote the number of such roots by n roots .
In our numerical experiments, we compute the eigenvalues of the colleague matrix using three different algorithms: our Algorithm 4; MATLAB's eig function; and MATLAB's eig function with balancing turned off (using the option 'nobalance'), which we call eig nb. For our experiments in extended (quadruple) precision, we use the Advanpix Multiprecision Computing Toolbox and its implementation of eig (see [1]). Since the Advanpix Multiprecision Computing Toolbox's eig function always balances the matrix, and does not support the 'nobalance' option, we omit the test of eig nb in extended precision.
For each example, we report the degree of the underlying polynomial, the order n of the Chebyshev expansion used to approximate it, the size of the vector c in the Euclidean norm, the Frobenius norm of the completely balanced colleague matrix, which we denote by bal(C), the number n roots of computed roots inside the region (188) for the given value of δ > 0, the size max i |z i | of the largest complex root of the colleague matrix, and the value of max i |η(p; x i )|, where the maximum is taken over all of the real parts of the computed roots inside the region (188). We

p rand (x): Polynomials with Random Coefficients
Following [14], we construct polynomials p rand (x) by sampling Chebyshev expansion coefficients a i independently from a standard normal distribution, so that a i ∼ N (0, 1), for i = 0, 1, . . . , n − 1. Then, we choose the desired value of c by setting a n = a / c , so that the vector of coefficients c appearing in the colleague matrix (17), where c i = a i /a n for i = 1, 2, . . . , n, has the specified norm. For this example, we choose n = 30 and set δ = 10 −5 to extract the real roots (see formula (188)).
We report the results in Figure 3. We see that our algorithm shows the expected backward stability over the entire range of c , while MATLAB's eig, both balanced and unbalanced, shows the expected growth with c (see the discussion in Section 2.4). Interestingly, for this example, balancing appears to only improve the error by an order of magnitude or two, while leaving the growth in the error with respect to c unchanged. This turns out be completely consistent the following explanation. Balancing the colleague matrix can reduce the magnitude of the all elements from c to c 1 2 , except for the element in the (n, n)-position, which balancing cannot change. In this example, all of the elements of the vector c are around the same size as c , so there are n large elements in the colleague matrix. Thus, balancing reduces the number of large elements of size c from n to 1, resulting in an n-fold reduction in the norm of the matrix. In this example, n = 30, which corresponds well with the approximately 30-fold reduction in error due to balancing that we observe in Figure 3.
We found that the colleague matrix has a single large eigenvalue of the size c , and the rest of the eigenvalues are small. This is not surprising, since there is an entry of size c is the (n, n)-position of the matrix (from which it follows that Ce n ≈ c e n ). Thus, for all three algorithms, max i | z i | ≈ c .

p wilk (x): Wilkinson's Polynomial
Here we consider the famous Wilkinson polynomial, normalized so that all of its roots are inside the interval [−1, 1], defined by the formula We construct an order-n Chebyshev expansion of this degree-m polynomial, sampling it at n Chebshev points and applying a linear transformation to obtain the expansion coefficients (see [36]). We then compute the eigenvalues of the colleague matrix, and set δ = 10 −3 to extract the real roots (see formula (188)). The results of our numerical experiment are shown in Tables 1 and 2. We observe that our algorithm is backward stable, while eig loses accuracy whenever the roots of the colleague matrix are large. Plots of the real and complex roots of the order-100 Chebyshev expansion are shown for various degrees of p wilk (x) in Figure 4. When the order of the Wilkinson polynomial becomes large, spurious real roots begin appearing in the middle of the interval [−1, 1]. It turns out that the roots computed by our algorithm are still backward stable, even in this situation; the function is so small near the middle of the interval that a small relative perturbation in the Chebyshev coefficients causes additional roots to appear.
Remark 5.1. The remarkable stability of eig for many of the examples in Tables 1 and 2 is explained by the following observation. The colleague matrix is the sum of a tridiagonal matrix and a matrix that is all zeros except for the last row, which is essentially equal to the coefficient vector c (see formula (17)). When there are large elements of c near the tail of the vector, the corresponding large entries in the colleague matrix cannot be balanced away, since they are very close to the diagonal of the matrix. On the other hand, when all of the elements of c near the tail of the vector are relatively small, and the large elements of c appear near the head of the vector, these large elements can be easily balanced away, since they are far from the diagonal, and the corresponding elements on the other side of the diagonal are all zero. The coefficient vector c is usually large only because the last coefficient of the corresponding non-monic Chebyshev expansion is small. If the function being approximated by this non-monic Chebyshev expansion has been adequately represented, then taking additional terms in the expansion will result in corresponding expansion coefficients which are all machine epsilon in size. Thus, adding terms to the Chebyshev expansion has the effect of adding elements of size approximately one to the tail of the coefficient vector c; if enough such elements are added, then all the large elements of c will be closer to the head of the vector, and can be balanced away. We also observe that, not unexpectedly, the size of the largest eigenvalue of the colleague matrix is approximately the same size as the norm of the colleague matrix after balancing. Thus, if enough terms are taken in a Chebyshev expansion, all of the  24 24 0.95 · 10 4 0.65 · 10 1 0.92 · 10 0 24 0.15 · 10 −14 24 0.32 · 10 −14 25 0.66 · 10 15 0.35 · 10 11 0.35 · 10 11 14 † 0.11 · 10 −4 24 0.19 · 10 −14 26 0.22 · 10 15 0.11 · 10 6 0.76 · 10 5 24 0.51 · 10 −10 24 0.24 · 10 −14 27 0.40 · 10 16 0.66 · 10 4 0.38 · 10 4 24 0.90 · 10 −12 24 0.19 · 10 −14 28 0.12 · 10 15 0.37 · 10 3 0.17 · 10 3 24 0.58 · 10 −13 24 0.14 · 10 −14 100 0.38 · 10 14 0.24 · 10 2 0.11 · 10 1 24 0.10 · 10 −13 24 0.24 · 10 −14 34 100 0.34 · 10 14 0.17 · 10 2 0.11 · 10 1 34 0.14 · 10 −13 34 0.12 · 10 −13 44 100 0.32 · 10 14 0.16 · 10 2 0.11 · 10 1 44 0.53 · 10 −13 44 0.41 · 10 −14 54 100 0.30 · 10 14 0.16 · 10 2 0.11 · 10 1 60 0.15 · 10 −13 60 0.28 · 10 −13 Table 1: The results of computing the roots of the Wilkinson polynomial p wilk (x), using our algorithm and eig, with δ = 10 −3 . † The error was so large here that some roots were outside of the region (188).  Table 2: The results of computing the roots of the Wilkinson polynomial p wilk (x) in extended precision, using our algorithm and eig, with δ = 10 −3 . † The error was so large here that some roots were outside of the region (188).
eigenvalues of the colleague matrix will eventually be small. All of this indicates that, provided enough terms are taken, a dense eigensolver combined with balancing can result in a backward stable rootfinding algorithm. Of course, we note that balancing the colleague matrix destroys the Hermitian plus rank-1 structure, which bars the use of structured QR algorithms depending on this property.

f sin (x): A Smooth Function
Here we construct an order-n Chebyshev expansion of the smooth function f sin (x) = sin(2 + 20(x + 0.222) 2 ).
Since f sin (x) is analytic, its expansion coefficients decay exponentially. When i ≥ 80, the coefficients a i are around 10 −14 in size (see Figure 5). If the coefficients are computed in extended precision, then, when i ≥ 125, the coefficients a i are around 10 −34 in size.
Since the function is approximated accurately by a Chebyshev expansion, its roots can be computed from the corresponding colleague matrix (see, for example, [12] for a nice discussion). The results of our numerical experiment are shown in Tables 3 and 4. Plots of the real and complex roots of the order-100 Chebyshev expansion are shown in Figure 6.

p mult (x): A Polynomial with Multiple Roots
Here we construct an order n Chebyshev expansion of the degree m polynomial This polynomial has four simple roots on the interval [−1, 1], and a root of multiplicity (m − 4) at the point 1 − 10 −3 (see Figure 7). The results of our numerical experiments are shown in Tables 5 and 6. We observe that, in double precision, when the multiplicity of the root is greater than or equal to 5, not all real roots are found. This is because the error in these roots is approximately equal to 1 5 , and when ≈ 10 −14 , we have that 1 5 ≈ 1.6 · 10 −3 ; when δ = 10 −3 , this means that some of these roots can be outside the region (188). Likewise, since when ≈ 10 −34 , 1 12 ≈ 1.4 · 10 −3 , it follows that in extended precision, real roots are missed when their multiplicity is greater than or equal to approximately 12. See the excellent discussion in [13] for more details.
bottom right corner of the corresponding colleague matrix is around 10 15 in size. This polynomial has seven real roots on the interval [−1, 1], and a single large imaginary root. We report the results of our numerical experiment in Table 7. Clearly, eig struggles to produce any accuracy at all, while our algorithm returns all the roots to machine precision.
at least 1430 to resolve it. Our numerical experiments are shown in Table 8. Plots of the real and complex roots of the order-1600 Chebyshev expansion are shown in Figure 9.

CPU Times
The CPU times of our algorithm are compared to the times of MATLAB's eig in Figure 10.
These timing experiments were performed on polynomials with random, independent, normally distributed Chebyshev expansion coefficients, with the last coefficient chosen so that the vector c has the desired norm (see Section 5.1). We found that the CPU times do not depend on c , so we report the results only for c = 2. We observe that our algorithm is strictly faster than eig, even for small inputs, except perhaps for n = 7, for which our algorithm and eig cost about the same. The growth in CPU times taken by our algorithm agrees nicely with the expected asymptotic cost of O(n 2 ), while eig shows a growth of O(n 3 ).

Conclusions and Generalizations
In this manuscript, we describe an explicit, O(n 2 ) structured QR algorithm for colleague matrices (more generally, for Hessenberg matrices that are the sum of a Hermitian matrix and a rank-1 matrix), and prove that it is componentwise backward stable. These results can be generalized in several directions, of which we describe four. First, the algorithm can be modified in a fairly straightforward way to work on Hessenberg matrices that are the sum of a Hermitian matrix and a rank-k perturbation (as opposed to a rank-1 perturbation). Like in the rank-1 case, most of the entries in the Hermitian part are inferred from the low rank part, except that they are inferred from a rank-k matrix instead of a rank-1 matrix. The QR iteration proceeds similarly, the main difference being that the correction in Line 12 of Algorithm 1 becomes a correction to a row of an n × k matrix. Second, the extension of this algorithm to an implicit, O(n 2 ) structured QR algorithm that is also componentwise backward stable is fairly straightforward. The key observation of this manuscript (that, to maintain componentwise error bounds, a correction must be applied to the rank-1 part whenever an entry of the matrix is eliminated) can be applied to a bulge-chasing algorithm where the matrix is similarly represented by generators. Third, we observe that this algorithm can be used to accelerate the calculation of eigenvalues of general matrices, not necessarily in Hessenberg form, that are representable as the sum of a Hermitian matrix and a rank-1 (or rank-k) matrix. Such matrices can be quickly reduced to Hessenberg form in O(n 3 ) operations, and once they are in Hessenberg form, our O(n 2 ) algorithm can be used to compute the eigenvalues. Thus, the cost of the algorithm is dominated by the reduction to Hessenberg form, which will have a much smaller constant than the standard algorithm for the evaluation of the eigenvalues of the original dense matrix. Furthermore, if the reduction to Hessenberg form can be done in a componentwise backward stable fashion, then this scheme results in a componentwise backward stable eigensolver for general matrices of the form Hermitian plus rank-1 (or rank-k).
Fourth, we observe that our algorithm can be used to find the roots of polynomials expressed in other bases besides Chebyshev polynomials. It was observed in [5] that, given any orthogonal polynomial basis that satisfies a three-term recurrence relation, and given a polynomial expressed in that basis, it is possible to construct an analogue of the colleague matrix from the expansion coefficients. This matrix is a Hessenberg matrix that is the sum of a (not necessarily symmetric) tridiagonal matrix and a rank-1 matrix; matrices of this form are called comrade matrices. For all classical orthogonal polynomials, the tridiagonal part can made symmetric by balancing, without making any entries of the matrix much larger or much smaller. Our algorithm can then be applied to this new matrix, which is a Hessenberg matrix that is the sum of a symmetric tridiagonal matrix and a rank-1 matrix.