Abstract
We design a fast implicit real QZ algorithm for eigenvalue computation of structured companion pencils arising from linearizations of polynomial rootfinding problems. The modified QZ algorithm computes the generalized eigenvalues of an \(N\times N\) structured matrix pencil using O(N) flops per iteration and O(N) memory storage. Numerical experiments and comparisons confirm the effectiveness and the stability of the proposed method.
Similar content being viewed by others
Notes
MATLAB is a registered trademark of The Mathworks, Inc.
Both implementations are available for download at http://www.unilim.fr/pages_perso/paola.boito/software.html.
References
Aurentz, J.L., Mach, T., Robol, L., Vandebril, R., Watkins, D.S.: Roots of polynomials: on twisted QR methods for companion matrices and pencils, arXiv:1611.02435 [math.NA] (2016)
Austin, A.P., Kravanja, P., Trefethen, L.N.: Numerical algorithms based on analytic function values at roots of unity. SIAM J. Numer. Anal. 52, 1795–1821 (2014)
Betcke, T., Higham, N.J., Mehrmann, V., Schröder, C., Tisseur, F.: NLEVP: a collection of nonlinear eigenvalue problems. ACM Trans. Math. Softw. 39, 7–28 (2013)
Boito, P., Eidelman, Y., Gemignani, L.: Implicit QR for companion-like pencils. Math. Comput. 85, 1753–1774 (2016)
Corless, R.: Generalized companion matrices for the Lagrange basis. In: Proceedings EACA (2004)
de Boor, C.: An empty exercise. ACM SIGNUM Newslett. 25(4), 2–6 (1990)
Effenberger, C.: Robust solution methods for nonlinear eigenvalue problems, Ph.D. thesis, EPFL (2013)
Eidelman, Y., Gohberg, I.: On a new class of structured matrices. Integral Equ. Oper. Theory 34, 293–324 (1999)
Eidelman, Y., Gohberg, I., Haimovici, I.: Separable Type Representations of Matrices and Fast Algorithms. Volume 1. Basics. Completion Problems. Multiplication and Inversion Algorithms, Operator Theory: Advances and Applications. Birkhäuser, Basel (2013)
Eidelman, Y., Gohberg, I., Haimovici, I.: Separable Type Representations of Matrices and Fast Algorithms. Volume 2. Eigenvalue method, Operator Theory: Advances and Applications. Birkhäuser, Basel (2013)
Fiedler, M., Markham, T.L.: Completing a matrix when certain entries of its inverse are specified. Linear Algebra Appl. 74, 225–237 (1986)
Golub, G.H., Van Loan, C.F.: Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences, 3rd edn. Johns Hopkins University Press, Baltimore, MD (1996)
Jenkins, M.A., Traub, J.F.: Principles for testing polynomial zerofinding programs. ACM Trans. Math. Softw. 1, 26–34 (1975)
Lawrence, P.W.: Eigenvalues methods for interpolation bases. University of West Ontario, Electronic Thesis and Dissertation Repository, paper 1359 (2013)
Lawrence, P.W.: Fast reduction of generalized companion matrix pairs for barycentric Lagrange interpolants. SIAM J. Matrix Anal. Appl. 34(3), 1277–1300 (2013)
Lawrence, P.W., Corless, R.M.: Stability of rootfinding for barycentric Lagrange interpolants. Numer. Algorithms 65(3), 447–464 (2014)
Solov’ëv, S.I.: Preconditioned iterative methods for a class of nonlinear eigenvalue problems. Linear Algebra Appl. 415, 210–229 (2006)
Toh, K.-C., Trefethen, L.N.: Pseudozeros of polynomials and pseudospectra of companion matrices. Numer. Math. 68, 403–425 (1994)
Watkins, D.S.: Fundamentals of Matrix Computations. Pure and Applied Mathematics (New York), 2nd edn. Wiley-Interscience [John Wiley & Sons], New York (2002)
Watkins, D.S.: The Matrix Eigenvalue Problem: GR and Krylov Subspace Methods. SIAM, Philadelphia (2007)
Acknowledgements
Thanks to Thomas Mach for useful suggestions concerning the Fortran implementation of Givens transformations.
Author information
Authors and Affiliations
Corresponding author
Additional information
This work was partially supported by GNCS-INDAM and University of Pisa.
Appendix
Appendix
In this appendix we give a formal proof of the correctness of the algorithm stated in Sect. 4. Specifically, we prove the following:
Theorem 1
Let \((A,B)\in {\mathcal {P}}_N\) be a matrix pair with an upper Hessenberg matrix \(A=V-{{\varvec{z}}}{{\varvec{w}}}^*\) from the class \({\mathcal {H}}_N\) and an upper triangular matrix \(B=U-{{\varvec{p}}} {{\varvec{q}}}^*\) from the class \({{\mathcal {T}}}_{N}\) with the unitary matrices \(V\in {{\mathcal {V}}}_{N},U\in {{\mathcal {U}}}_{N}\) and the vectors \({{\varvec{z}}},{{\varvec{w}}},{{\varvec{p}}},{{\varvec{q}}}\in {\mathbb {R}}^N\). Let \(p(z)=\alpha + \beta z +\gamma z^2 \in {\mathbb {R}}[z]\) be a polynomial of degree at most 2. Let Q, Z be unitary matrices defined as in (2.8), (2.9) where the matrices \(Q_i\) and \(Z_i\), \(1\le i\le N-1\), are generated by the algorithm in Sect. 4. Then \(A_1=Q^** A Z\) and \(B_1=Q^* B Z\) are upper Hessenberg and upper triangular, respectively, and, moreover, \(Q^* p(A B^{-1}) {{\varvec{e}}}_1=\alpha {{\varvec{e}}}_1\) for a suitable scalar \(\alpha \in {\mathbb {R}}\).
Proof
The property \(Q^* p(A B^{-1}) {{\varvec{e}}}_1=\alpha {{\varvec{e}}}_1\) easily follows by construction. The proof of the remaining properties is constructive by showing that \(A_1\) is upper Hessenberg and \(B_1\) is upper triangular and then providing structured representations of the entries of their unitary components \(V_1=Q^* V Z\) and \(U_1=Q^* U Z\). We restrict ourselves to consider \(A_1\) and \(V_1\) since the computation of \(B_1\) and \(U_1\) and of the perturbation vectors can be treated in a similar way.
We treat A and V as block matrices with entries of sizes \(m^A_i\times n^A_j,\;i,j=1,\dots ,N+3\), where
Relative to this partition the matrix V has diagonal entries
upper quasiseparable generators
and lower quasiseparable generators
Relative to the partition (6.32) the matrix A is a block upper triangular matrix with diagonal entries
Moreover using (2.1) we obtain upper quasiseparable of the matrix A relative to the partition (6.32) with orders
by the formulas
Using (2.8) and setting
we get
We have
with
and
with
We treat the lower Hessenberg matrix \(Q^*\) as a block matrix with entries of sizes \(\tau ^A_i\times m^A_j,\;i,j=1,\dots ,N+3\), where
The matrix \(Q^*\) has the representation considered in Lemma 31.1 in [10] with the matrices \(S_k\;(k=1,\dots ,N+2)\) of sizes \((\tau ^A_1+r^S_1) \times m^A_1,\; (\tau ^A_k+r^S_k)\times (m^A_k+r^S_{k-1})\;(k=2,\dots ,N+1),\;\tau ^A_{N+2}\times (m^A_{N+2}+r^S_{N+1})\), where
We treat the upper Hessenberg matrix Z as a block matrix with entries of sizes \(n^A_i\times \nu ^A_j,\;i,j=1,\dots ,N+2\), where
The matrix Z has the representation considered in Lemma 31.1 in [10] with the matrices \(T^A_k\; (k=1,\dots ,N+2)\) of sizes \(n^A_1\times (\nu ^A_1+r^*_1),\;(n^A_k+r^*_{k-1})\times (\nu ^A_k+r^*_k)\; (k=2,\dots ,N+1),\;(n^A_{N+2}+r^*_{N+1})\times \nu ^A_{N+2}\), where
Now we apply the structured multiplication algorithm for quasiseparable representations stated in Corollary 31.2 in [10] in order to determine diagonal entries \({\tilde{d}}_A(k)\;(k=3,\dots ,N+2)\) and quasiseparable generators \({\tilde{q}}(j)\;(j=3,\dots ,N+1);\;{\tilde{g}}_A(i)\;(i=3,\dots ,N+2)\) of the matrix \(A_1=Q^*AZ\) as well as auxiliary variables \({\beta }_k^A,f_k^A,\phi _k^A, \varphi _k^A\). The matrix \(A_1\) is obtained as a block one with entries of sizes \(\tau ^A_i\times \nu ^A_j,\;i,j=1,\dots ,N+3\).
For the variables \({\beta }_k={\beta }^A_k,f_k=f^A_k,\phi _k=\phi ^A_k\) used in Corollary 31.2 we use the partitions
with the matrices \(f^A_k,\varphi ^A_k,\chi _k\) of sizes \(2\times 2,2\times r^V_k, 2\times 1\). For \(k=1,\dots ,N-2\) combining Corollary 31.2 with (6.42), (6.43), (6.36) and (6.38),(6.39) we get
Using (6.46) and (6.39) we get
and
with \(\epsilon ^A_{k+2}\) as in (4.14).
Inserting (6.48), (6.49) in (6.47) and using (6.46), (6.38) we obtain
From (6.50) using (4.15) we obtain the relations
and
From (6.51) using (4.16) we have
and
The formulas (6.44) and (6.45) mean that \((\sigma ^A_k){(1)},\;k=1,\dots ,N-2\) are subdiagonal entries of the matrix \(A_1\) (treated as an usual scalar matrix). The equalities (6.54) imply that \(A_1\) is an upper Hessenberg matrix.
Next we apply the structured multiplication algorithm stated in Lemma 31.1 in [10] to compute (block) upper quasiseparable generators \({\tilde{g}}_V(i)\;(i=1,\dots ,N+2),\;{\tilde{h}}_V(j)\;(j=2,\dots ,N+3), \;{\tilde{b}}_V(k)\;(k=2,\dots ,N+2)\) with orders
and diagonal entries \({\tilde{d}}_{V_1}(k)\;(k=1,\dots ,N+3)\) of the matrix \(V_1=Q^*VZ\). The matrix \(V_1\) is obtained as a block one with entries of sizes \(\tau ^A_i\times \nu ^A_j,\;i,j=1,\dots ,N+3\), where the numbers \(\tau ^a_i,\nu ^A_j\) are defined in (6.44), (6.45).
Using Lemma 31.1 and (6.42), (6.43) we obtain that
together with the relation (4.26).
From (6.55) we find that the auxiliary matrices \({\alpha }_k\;(k=3,\dots ,N)\) satisfy the relations
Comparing this with (4.4), (4.21) we get
Thus using (6.56) and (6.55) we obtain (4.24), (4.25).
Next we show that the auxiliary variables \(f^A_k,\varphi ^A_k\;(k=3,\dots ,N)\) may be determined via relations (4.5), (4.28). Take \(\gamma _2\) as in (4.4) and assume that for some k with \(1\le k\le N-2\) the relations
hold. By (4.15) and (6.52) we have
Using (4.23) and (4.14) we get
Thus combining (6.58) and (6.59) together and using (4.24), (4.25) and (4.27) we obtain (4.28). \(\square \)
Rights and permissions
About this article
Cite this article
Boito, P., Eidelman, Y. & Gemignani, L. A real QZ algorithm for structured companion pencils. Calcolo 54, 1305–1338 (2017). https://doi.org/10.1007/s10092-017-0231-6
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10092-017-0231-6
Keywords
- Rank-structured matrix
- Quasiseparable matrix
- Real QZ algorithm
- Lagrange approximation
- Eigenvalue computation
- Complexity