{'fraction_non_alphanumeric': 0.05522393662950032, 'fraction_numerical': 0.014491121372526936, 'mean_word_length': 3.827473722363175, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Implicit QR for Companion-like Pencils

A fast implicit QR algorithm for eigenvalue computation of low rank corrections of unitary matrices is adjusted to work with matrix pencils arising from polynomial zerofinding problems . The modified QZ algorithm computes the generalized eigenvalues of certain NxN rank structured matrix pencils using O(N^2) ops and O(N) memory storage. Numerical experiments and comparisons confirm the effectiveness and the stability of the proposed method.


Introduction.
Computing the roots of a univariate polynomial is a fundamental problem that arises in many applications. One way of numerically computing the roots of a polynomial is to form its companion matrix (pencil) and compute the (generalized) eigenvalues.
In a paper on polynomial root-finding [13], Jónsson and Vavasis present a comparative analysis of the accuracy of different matrix algorithms. The conclusion is that computing the roots of a polynomial by first forming the associated companion pencil A − λB and then solving Ax = λBx using the QZ algorithm provides better backward error bounds than computing the eigenvalues of the associated companion matrix by means of the QR algorithm. The analysis does not take into account the possible use of balancing which would have the effect of allineating the accuracy of the two approaches [16,2]. However, it is worth noting that this use is not allowed if the customary QR and QZ algorithms are modified at the aim of reducing their complexity of one order of magnitude by exploiting the rank structure of the initial companion matrix (pencil). Thus, the use of the QZ algorithm applied to the companion pencil achieves the potential best score in terms of accuracy and efficiency among matrix methods for polynomial root-finding.
Another interesting appearance of the pencil approach is in the design of numerical routines for computing the determinant hence the zeros of polynomial matrices [15]. A classical technique is the interpolation of the determinant of the polynomial matrix A(λ) at the roots of unity via FFT. This gives a representation of p(λ) : = det A(λ) in the basis of Lagrange polynomials generated from the set of nodes. The main stumbling step in this process is the presence of undesirable infinite zeros or, equivalently, zero leading coefficients of p(λ). This makes the companion matrix approach fully unusable. On the contrary, a companion pencil can still be formed whose generalized eigenvalues are the finite roots of the polynomial. If the pencil is A − λB directly constructed from the coefficients of the polynomial expressed in the Lagrange basis then A is a unitary plus rank-one matrix and B is a rank-one perturbation of the identity matrix. By a preliminary reduction using the algorithm in [6] the matrices A and B can be transformed in upper Hessenberg and triangular form, respectively. This paper undertakes the task of developing an efficient variant of the QZ algorithm applied to a generalized companion pair (A, B), where A is upper Hessenberg, B is upper triangular and A and B are rank-one modifications of unitary matrices. These assumptions imply suitable rank structures in both A and B. We show that each matrix pair generated by the QZ process is also rank-structured. By exploiting this property a novel, fast adaptation of the QZ eigenvalue algorithm is obtained which require O(n) flops per iteration and O(n) memory space only. Numerical experiments confirm that the algorithm is stable.
The paper is organized as follows. In Sect. 2 we introduce the computational problem and describe the matrix structures involved. Sect. 3 and Sect. 4 deal with basic issues and concepts concerning the condensed representation and the invariance of these structures under the QZ process, whereas the main algorithm is presented in Sect. 5. Sect. 6 gives a backward error analysis for the QZ method applied to the companion pencil. In Sect. 7 we present an implementation of our fast variant of the QZ algorithm together with the results of extensive numerical experiments. Finally, the conclusion and a discussion are the subjects of Sect. 8.
2. The Problem Statement. Companion pencils and generalized companion pencils expressed in the Lagrange basis at the roots of unity are specific instances of the following general class.
Definition 2.1. The matrix pair (A, B), A, B ∈ C N ×N , belongs to the class P N ⊂ C N ×N × C N ×N of generalized companion pencils iff: There exist two vectors z ∈ C N and w ∈ C N and a unitary matrix V ∈ C N ×N such that 4. There exist two vectors p ∈ C N and q ∈ C N and a unitary matrix U ∈ C N ×N such that In order to characterize the individual properties of the matrices A and B we give some additional definitions.
Definition 2.2. We denote by T N the class of upper triangular matrices B ∈ C N ×N which are rank one perturbations of unitary matrices, i.e., such that (2.2) holds for a suitable unitary matrix U and vectors p, q.
Since B is upper triangular the strictly lower triangular part of the unitary matrix U in (2.2) coincides with the corresponding part of the rank one matrix pq * , i.e., where {p(i)} i=1,...,N and {q(j)} j=1,...,N are the entries of p and q, respectively. Definition 2.3. We denote by U N the class of unitary matrices U ∈ C N ×N which satisfy the condition (2.3), i.e., for which there exist vectors p, q such that the matrix B = U − pq * is an upper triangular matrix.
From the nullity theorem [11], see also [9, p.142] it follows that the same property also holds in the strictly upper triangular part, namely, Definition 2.4. We denote by H N the class of upper Hessenberg matrices A ∈ C N ×N which are rank one perturbations of unitary matrices, i.e., such that (2.1) holds for a suitable unitary matrix V and vectors z, w.
Definition 2.5. We denote by V N the class of unitary matrices V ∈ C N ×N for which there exist vectors z, w such that the matrix A = V − zw * is an upper Hessenberg matrix.
We find that Again from the nullity theorem it follows that a similar property also holds in the upper triangular part, namely, The QZ algorithm is the customary method for solving generalized eigenvalue problems numerically by means of unitary transformations (see e.g. [12] and [18]). Recall that the Hessenberg/triangular form is preserved under the QZ iteration; an easy computation then yields Indeed if Q and Z are unitary then from (2.1) and (2.2) it follows that the matrices A 1 = Q * AZ and B 1 = Q * BZ satisfy the relations with the unitary matrices V 1 = Q * V Z, U 1 = Q * U Z and the vectors z 1 = Q * z, w 1 = Z * w, p 1 = Q * p, q 1 = Z * q. Moreover one can choose the unitary matrices Q and Z such that the matrix A 1 is upper Hessenberg and the matrix B 1 is upper triangular. Thus, one can in principle think of designing a structured QZ iteration that, given in input a condensed representation of the matrix pencil (A, B) ∈ P N , returns as output a condensed representation of (A 1 , B 1 ) ∈ P N generated by one step of the classical QZ algorithm applied to (A, B). In the next sections we first introduce an eligible representation of a rank-structured matrix pencil (A, B) ∈ P N and then discuss the modification of this representation under the QZ process.
3. Quasiseparable Representations. In this section we exploit the properties of quasiseparable representations of rank-structured matrices [8], [9,Chapters 4,5]. First we recall some general results and definitions. Subsequently, we describe their adaptations for the representation of the matrices involved in the structured QZ iteration applied to an input matrix pencil (A, B) ∈ P N .

Preliminaries.
is (r L , r U )-quasiseparable, with r L , r U positive integers, if, using MATLAB * notation, Roughly speaking, this means that every submatrix extracted from the lower triangular part of M has rank at most r L , and every submatrix extracted from the upper triangular part of M has rank at most r U . Under this hypothesis, M can be represented using O((r L + r U )N ) parameters. In this subsection we present such a representation.
The quasiseparable representation of a rank-structured matrix consists of a set of vectors and matrices used to generate its entries. For the sake of notational simplicity, generating matrices and vectors are denoted by a lower-case letter.
In this representation, the entries of M take the form  (2.4) it follows that any matrix from the class U N has upper quasiseparable generators with orders equal one.
The quasiseparable representation can be generalized to the case where M is a block matrix, and to the case where the generators do not all have the same size, provided that their product is well defined. Each block M ij of size m i × n j is represented as in (3.1), except that the sizes of the generators now depend on m i and n j , and possibly on the index of a and b. More precisely: p(i), q(j), a(k) are matrices of sizes m i ×r L i−1 , r L j ×n j , r L k ×r L k−1 , respectively; d(i) (i = 1, . . . , N ) are m i × n i matrices, g(i), h(j), b(k) are matrices of sizes m i ×r U i , r U j−1 ×n j , r U k−1 ×r U k , respectively. The numbers r L k , r U k (k = 1, . . . , N − 1) are called the orders of these generators. It is worth noting that lower and upper quasiseparable generators of a matrix are not uniquely defined. A set of generators with minimal orders can be determined according to the ranks of maximal submatrices located in the lower and upper triangular parts of the matrix.
One advantage of the block representation for the purposes of the present paper consists in the fact that N × N upper Hessenberg matrices can be treated as (N + 1 × N + 1) block upper triangular ones by choosing block sizes as m 1 = · · · = m N = 1, m N +1 = 0, n 1 = 0, n 2 = · · · = n N +1 = 1. (3.2) Such a treatment allows also to consider quasiseparable representations which include the main diagonals of matrices. Assume that C is an N × N scalar matrix with the entries in the upper triangular part represented in the form are called upper triangular generators of the matrix C with orders r k (k = 1, . . . , N ). From (2.5) it follows that any matrix from the class V N has upper triangular generators with orders not greater than two. If we treat a matrix C as a block one with entries if sizes (3.2) we conclude that the elements g(i) Matrix operations involving zero-dimensional arrays (empty matrices) are defined according to the rules used in MATLAB and described in [3]. In particular, the product of a m × 0 matrix by a 0 × m matrix is a m × m matrix with all entries equal to 0. Empty matrices may be used in assignment statements as a convenient way to add and/or delete rows or columns of matrices.

3.2.
Representations of matrix pairs from the class P N . Let (A, B) be a matrix pair from the class P N . The corresponding matrix A from the class H N is completely defined by the following parameters: . From (2.5) it follows that the matrix V ∈ V N has upper triangular generators with orders not greater than two.
The corresponding matrix B from the class T N is completely defined by the following parameters: (2.4) it follows that the matrix U ∈ U N has upper quasiseparable generators with orders equal one.
All the given parameters define completely the matrix pair (A, B) from the class P N .
Each step of structured QZ should update these parameters while maintaining the minimal orders of generators. However, structured algorithms for the multiplication of quasiseparable matrices may output redundant generators for the product matrix. For this reason, we will need an algorithm that compresses generator to minimal order.

The compression of generators.
In this section we present an algorithm that takes as input the possibly redundant quasiseparable generators of a unitary matrix and outputs generators of minimal order.
Compression algorithm Input: lower quasiseparable generators p(j), q(j), a(j) and upper quasiseparable generators g(j), h(j), b(j) of possibly redundant orders for the matrix U .
Output: upper quasiseparable generators g s (j), h s (j), b s (j) of minimal order for U .
1. Set X 0 , Y 0 , z 0 to be the 0 × 0 empty matrices and p (1) Observe that the submatrix Θ k ∆ k has orthonormal columns.

Determine the (s
7. Compute the matrices Y k of the size s k × r U k and z k of the size s k × ρ k by the formulas This algorithm is justified in a similar way as the compression algorithm [9, Theorem 7.5] but with computations in the forward direction. For a matrix from the class U N we have m i = n i = 1, i = 1, . . . , N and r L k = 1, k = 1, . . . , N − 1 which by virtue of (4.1) implies ρ k = s k = 1, k = 1, . . . , N − 1.
For a matrix from the class V N we determine the sizes of blocks via (3.2) and the orders of lower generators r L k = 1, k = 1, . . . , N . Hence using (4.1) we obtain Remark 4.1. With the above hypotheses, the matrix U admits the factorization where W is a block lower triangular unitary matrix with block entries of sizes m i × ν j (i, j = 1, . . . , N ) and F is a block upper triangular unitary matrix with block entries of sizes ν i × n j (i, j = 1, . . . , N ). Moreover one can choose the matrix W in the form and, in addition, some initial conditions are satisfied. In the case of the single-shift implicit QZ step the initial condition is where α ∈ C is the shift. In this case one obtains the unitary Hessenberg matrices Q and Z in the form and Q i , Z i are complex Givens rotation matrices. The condition (5.2) means In the sequel of this section we present a fast adaptation of the implicit singleshift QZ algorithm for an input matrix pair (A, B) ∈ P N . The modified algorithm works on the generators of the two matrices and this explains why it is referred to as a structured implicit QZ iteration with single shift. More specifically, the input pair (A, B) = (A 0 , B 0 ) is first represented by means of a linear set of generators as explained in Subsection 3.2. Then a high-level description of structured QZ iteration goes as follows.
1. Given the number α = α 0 ∈ C perform one step of the implicit single-shift QZ algorithm (2.6) by computing subdiagonal entries of the matrix A 1 , diagonal entries of the matrix B 1 , vectors of perturbation z 1 , w 1 , p 1 , q 1 as well as upper triangular generators of the matrix V 1 and upper quasiseparable generators of the matrix U 1 (with redundant orders). 2. Compress the representations of V 1 and U 1 by using the compression algorithm in Section 4.
We now describe how the QZ computation at step 1 above is efficiently performed by working on the set of generators. Let Then the unitary 2 × 2 matrices Q k , Z k (k = 1, . . . , N − 1) generated from one step of the implicit QZ algorithm with single shift α ∈ C, as well as subdiagonal entries σ A1 k (k = 1, . . . , N − 1) of the matrix A 1 , upper triangular generators g (1) V (k) (k = 1, . . . , N − 1) with orders not greater than two of the matrix V 1 , diagonal entries d U (k) (k = 2, . . . , N − 1) of the matrix U 1 with orders not greater than one and vectors of perturbation z 1 = col(z (1) are determined via the following algorithm.
Structured QZ algorithm Input: subdiagonal entries σ A k for the matrix A, upper triangular generators g V (i), h V (i), b V (k) for the matrix V , diagonal entries d B (k) for the matrix B, upper quasiseparable generators g U (i), h U (i), b U (k) for the matrix U , and perturbation vectors p, q, z and w.
Output: subdiagonal entries σ A1 k for the matrix A 1 , upper quasiseparable generators g (1) V (k) for the matrix V 1 , diagonal entries d U (k) for the matrix U 1 , perturbation vectors p 1 , q 1 , z 1 , w 1 , and, if needed, the matrices Q k and Z k .

Set
and for k = 3, . . . , N perform the following. Determine the 2 × 2 unitary matrix W k and the number X k such that The submatrix Θ k ∆ k has orthonormal columns and one can determine the 3 × 3 unitary matrix F k from the condition Determine the matrices h F (k), d F (k), b s (k), g s (k) of sizes 2×1, 1×1, 2×2, 1×2 from the partition U (1) = 1 (5.63) and for k = 2, . . . , N perform the following. Determine the 2 × 2 unitary matrix W k and the number X k such that and determine the numbers h ′ k , h ′′ k , ∆ k , Ω k from the partition The 2-dimensional column Θ k ∆ k has unit norm and one can determine U (k) from the partition The proof of this algorithm is a combination of arguments in [10, Theorem 31.4] with changing the product Q * AQ by Q * AZ, and in [10, Theorem 36.4].
Here η denotes the computational cost required to compute and apply a 2 × 2 Givens matrix.
6. Backward error analysis. The backward stability of the QZ method implies that QZ applied to a companion pencil (A, B) computes the exact eigenvalues of a slightly perturbed pencil (A + E, B + G), where the matrices E and G have small norm. Does this imply that the characteristic polynomial of (A + E, B + G) is close to the characteristic polynomial of (A, B)? In other words, is QZ applied to a companion pencil backward stable as a rootfinding method?
The corresponding problem for the companion matrix is examined in [5], with a positive answer (Theorem 2.1). Given a monic polynomial p(x) = n j=0 a j x j , let A be its companion matrix and E the perturbation matrix that measures the backward error of the QR method applied to A. Edelman and Murakami show that the coefficients of s k−1 in the polynomial det(sI n − A − E) − det(sI n − A) are given at first order by the expression with a n = 1. We seek to extend this result to companion pencils, at least for the case of nonsingular B: here a n is no longer necessarily equal to 1, although different from 0. First, observe that (6.1) can be generalized to the computation of det(s(I n + G) − (A + E)). Indeed, at first order we have and therefore det(s(I n + G) − (A + E)) = det(I n + G)det(sI n − A − E + GA), (6.2) where (again at first order) det(I n + G) = 1 + tr(G) (6.3) and det(sI n − A − E + GA) can be computed by applying (6.1) to the pencil sI n − A with the perturbation matrix E − GA. In order to fix the notation, let us write at first order Now, recall that B can be seen as a rank-one perturbation of the identity matrix: we can write B = I n + (a n − 1)e n e T n , with the usual notation e n = [0, . . . , 0, 1] T . So, the perturbed companion pencil is s(B + G) − (A + E) = s(I n + G) − (A + E) + s(a n − 1)e n e T n . We know how to compute the determinant in the right-hand side of (6.6) thanks to the discussion above. Now we want to compute the second factor. What we need is the (n, n) entry of the matrix [s(I n + G) − (A + E)] −1 , which can be written as where q(x) = x n + n−1 j=0 a j x j . Finally, from det(s(B+G)−(A+E)) we subtract the quantity det(sB−A) = p(s), and we obtain the following result: Proposition 6.1. With the above notation, the following equality is correct at first order: (6.8)

A geometric approach.
A first-order estimate of det(A + E − s(B + G)) − det(A − sB) can also be obtained via a geometric approach by generalizing the ideas present in [5]; see also [13] and [17]. Note that this approach is related to the Leverrier rootfinding method: see [14] for an adaptation to the matrix pencil case.
We summarize the main idea. The ingredients here are: • A polynomial p(z) = n j=0 a j z j .
• The associated companion pencil A − λB of size n × n, with B invertible. We take here • The Sylvester affine space S of companion pencils through A − λB: it is the set of pencils W − λT where W is the zero matrix except for the last column and T is the zero matrix except for the (n, n)-entry. Clearly dim S = n + 1.
, where E and G are unstructured perturbation matrices. We want to compute at first order the coefficients (δa) 0 , . . . , (δa) n of the polynomialp(z) − p(z). Given full, unstructured perturbation matrices E and G, the goal is to characterize their components E S and G S along the Sylvester space.
A basis for the normal space N is given by {U k − λV k } k=1,...,n , where . and the matrices {M k } k=1,...,n are defined by the recursion M k = CM k+1 + a k an I n , M n = I n . From the special form of M k , and therefore of U k and V k , we deduce that E S (k, n) a n + a k−1 a n G S (n, n). (6.9) Also, observe that dim(T ∩ S) = 1: indeed, the equation tr(U H k E S + V H k G S ) = 0 has a one-dimensional space of solutions. So, if we fix G S (n, n), we have at first order (δa) k−1 = −E S (k, n) = −a n tr(U H k E + V H k G) + a k−1 G S (n, n). If we choose G S (n, n) = a n tr(B −1 G), which is the same perturbation on a n given by (6.8), then (6.9) yields E S (k, n)−a k−1 a n tr(B −1 G) = −a n tr(B −1 M k E)−a k−1 a n tr(B −1 G)+a n tr(B −1 M k−1 G), that is, at first order: (δa) k−1 = −a n tr(B −1 M k E) + a n tr(B −1 M k−1 G).
Using the structure of M k and therefore of a n B −1 M k we obtain at first order a m E 1,m + (a n − 1)E 1,n , (6.10) G i,i+m−k+1 − (a n − 1)G k−1,n (6.11) for k = 2, . . . , n, and (δa) n = a n tr(B −1 G) as specified above. Note that, if the matrices E and G are chosen as in Example 7.3, then the backward errors predicted by (6.8) and (6.11) are the same.
7. Numerical results. In this section we test the performance of the fast QZ method presented above and compare it to classical unstructured QZ applied to the companion pencil and to classical QR applied to the companion matrix. The experiments are performed in Matlab: we use the commands roots for classical QR and eig for classical QZ, whereas fast QZ is applied using our Matlab implementation of the algorithm described above. This implementation can be obtained from the authors upon request.
In particular, we report absolute forward and backward errors, measured in ∞norm. For each polynomial p(x) of degree N , the forward error is computed as where the λ j 's are the roots of p(x) as determined by the eigensolver that is being studied, and the α k 's are the "exact" roots of p(x), either known in advance or computed in high precision using Matlab's Symbolic Toolbox, unless otherwise specified. The backward error is computed as where the p k 's are the exact coefficients of p(x), either known in advance or computed in high precision from the known exact roots, and thep k 's are the coefficients computed in high precision from the λ j 's.
Our aim here is to provide experimental evidence pointing to the stability of our structured QZ method and to the better accuracy of QZ versus QR on some classes of polynomials. We do not report timings for the various methods, since the running times for our Matlab implementation cannot be compared to the running times of a built-in function such as eig or roots. Recall, however, that an arithmetic complexity estimate is given in Section 5.1.
Observe that the normalization of the polynomials is a crucial step for the proper functioning of QZ. Unlike the companion matrix, which is necessarily associated with a monic polynomial, the companion pencil allows for an arbitrary scaling of the polynomial. Unless otherwise specified, we normalize w.r.t. the 2-norm of the vector of coefficients. The polynomials obtained from computed roots, used for computation of backward errors, are also normalized in 2-norm. For the purpose of computing backward errors, we have also experimented with normalization in a least-squares sense, as suggested in [13], but in our examples we have generally found little difference between the 2-norm and the least-squares approach.
Example 7.1. Random polynomials. We apply our structured QZ method to polynomials whose coefficients are random complex numbers with real and imaginary parts uniformly chosen in [−1, 1]. Here N denotes the degree. Such test polynomials are generally well conditioned and they are useful to study the behavior of our method as the polynomial degree grows larger.
This is another set of wellconditioned polynomials for which we consider degrees up to 500. 1. p(x) = x 20 + x 19 + · · · + x + 1, 2. the polynomial with roots equally spaced in the interval [−2.1, 1.9], 3. the Chebyshev polynomial of first kind of degree 20, 4. the Bernoulli polynomial of degree 20. Table 7.3 shows forward errors for fast and classical QZ (after normalization of the polynomial) and for classical balanced QR (before normalization), as well as the maximum eigenvalue condition number.
Following [5], we also test our backward error analysis, given in Section 6, on these examples. Table 7.4 shows the logarithm in base 10, rounded to an entire number, of the computed and of the predicted backward errors for each coefficient of each polynomial. The predicted error is computed by applying (6.8).
As for the choice of E and G, we apply the backward error analysis given in [7], Theorem 4.1. This analysis implies that the backward error matrix for the QR method applied to a small rank perturbation of a Hermitian N × N matrix M is, roughly speaking, bounded by a small multiple of εN 2 M F , where ε is the machine epsilon. In view of this result, each nonzero entry of E and G can be taken as a small multiple of N ε, where N is the degree of the polynomial. For the sparsity pattern of E and G, we follow the ideas in [5]. So, E and G are chosen here as E=10*N*eps*triu(ones(20),-2) and G=10*N*eps*triu(ones(20),-1) in order to model the backward error introduced by the QZ iterations. Note that, under these assumptions, the experimental results confirm the theoretical analysis. Example 7.4. QZ vs. QR. In this example with heavily unbalanced coefficients, the (classical or structured) QZ method applied to the companion pencil computes the roots with better accuracy than QR applied to the companion matrix. We take here the polynomial p = 20 k=0 p k x k , where p k = 10 6(−1) (k+1) −3 for k = 0, . . . , 20. Table  7.5 shows absolute forward and backward errors for several eigenvalue methods, with or without normalizing the polynomial before the computation of its roots. See Table  7

Conclusions.
In this work we have presented and tested a new structured version of the single-shift, implicit QZ method. Our algorithm is designed for the fast computation of eigenvalues of matrix pencils belonging to a class P N that includes companion and Lagrange pencils. The quasiseparable structure of such pencils allows us to achieve quadratic complexity.
Future developments for the present work include writing a Fortran implementation of the proposed algorithm to support our cost analysis by timings comparisons. Our numerical experience says that the fast structured implementation of the QZ algorithm applied to a companion pencil provides a fast and provably backward stable root-finding method.