Abstract
The Schur decomposition of a square matrix A is an important intermediate step of state-of-the-art numerical algorithms for addressing eigenvalue problems, matrix functions, and matrix equations. This work is concerned with the following task: Compute a (more) accurate Schur decomposition of A from a given approximate Schur decomposition. This task arises, for example, in the context of parameter-dependent eigenvalue problems and mixed precision computations. We have developed a Newton-like algorithm that requires the solution of a triangular matrix equation and an approximate orthogonalization step in every iteration. We prove local quadratic convergence for matrices with mutually distinct eigenvalues and observe fast convergence in practice. In a mixed low-high precision environment, our algorithm essentially reduces to only four high-precision matrix-matrix multiplications per iteration. When refining double to quadruple precision, it often needs only 3–4 iterations, which reduces the time of computing a quadruple precision Schur decomposition by up to a factor of 10–20.
Similar content being viewed by others
Notes
For the matrix A and indices i1 ≤ i2, j1 ≤ j2, with A(i1 : i2,j1 : j2) we denote the submatrix of A consisting of all elements aij such that i ∈{i1,i1 + 1,…,i2} and j ∈{j1,j1 + 1,…,j2}. If i1 = i2, or j1 = j2, then one index can be dropped, e.g., A(i1,j1 : j2) = A(i1 : i1,j1 : j2) and \(A(i_{1}, j_{1}) = A(i_{1}:i_{1}, j_{1}:j_{1}) = a_{i_{1},j_{1}}\).
Note that Matlab’s Symbolic Toolbox vpa (variable precision arithmetic) supports eigenvalue computations but it does not support the computation of Schur decompositions. Executing eig in vpa takes 33.70 s for a 100 × 100 matrix in quadruple precision, compared to only 0.30 s needed by Advanpix for the complete Schur decomposition.
References
Advanpix: Multiprecision computing toolbox for Matlab. http://www.advanpix.com/. Accessed: 2021-12-15
Abdelfattah, A., Anzt, H., Boman, E. G., Carson, E., Cojean, T., Dongarra, J., Fox, A., Gates, M., Higham, N. J., Li, X. S., Loe, J., Luszczek, P., Pranesh, S., Rajamanickam, S., Ribizel, T., Smith, B. F., Swirydowicz, K., Thomas, S., Tomov, S., Tsai, Y. M., Yang, U. M.: A survey of numerical linear algebra methods utilizing mixed-precision arithmetic. Int J. High Perform Comput. Appl. 35(4), 344–369 (2021)
Ablin, P., Peyré, G.: Fast and accurate optimization on the orthogonal manifold without retraction. arXiv:2102.07432 (2021)
Absil, P.-A., Mahony, R., Sepulchre, R.: Optimization algorithms on matrix manifolds. Princeton University Press, Princeton (2008). With a foreword by Paul Van Dooren
Absil, P. -A., Mahony, R., Sepulchre, R., Van Dooren, P.: A grassmann-Rayleigh quotient iteration for computing invariant subspaces. SIAM Rev. 44(1), 57–73 (2002)
Ascher, U. M., Chin, H., Reich, S.: Stabilization of DAEs and invariant manifolds. Numer. Math. 67(2), 131–149 (1994)
Bai, Z., Day, D., Demmel, J. W., Dongarra, J. J.: A test matrix collection for non-Hermitian eigenvalue problems (release 1.0). Technical Report CS-97-355, Department of Computer Science, University of Tennessee, Knoxville, TN, USA, March 1997. Also available online from http://math.nist.gov/MatrixMarket
Beyn, W. -J., Kleß, W., Thümmler, V.: Continuation of low-dimensional invariant subspaces in dynamical systems of large dimension. In: Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems, pp 47–72. Springer, Berlin (2001)
Bindel, D. S., Demmel, J. W., Friedman, M.: Continuation of invariant subspaces in large bifurcation problems. SIAM J. Sci Comput. 30(2), 637–656 (2008)
Boumal, N.: An introduction to optimization on smooth manifolds To appear with Cambridge University Press (2022)
Braman, K., Byers, R., Mathias, R.: The multishift QR algorithm. I. Maintaining well-focused shifts and level 3 performance. SIAM J. Matrix Anal. Appl. 23(4), 929–947 (2002)
Braman, K., Byers, R., Mathias, R.: The multishift QR algorithm. II. Aggressive early deflation. SIAM J. Matrix Anal. Appl. 23(4), 948–973 (2002)
Chatelin, F.: Simultaneous Newton’s Iteration for the Eigenproblem. In: Defect Correction Methods (Oberwolfach, 1983), Volume 5 of Comput. Suppl., pp 67–74. Springer, Vienna (1984)
Davies, R.O., Modi, J.J.: A direct method for completing eigenproblem solutions on a parallel computer. Linear Algebra Appl. 77, 61–74 (1986)
Demmel, J. W.: Three methods for refining estimates of invariant subspaces. Computing 38, 43–57 (1987)
Eberlein, P. J.: A Jacobi-like method for the automatic computation of eigenvalues and eigenvectors of an arbitrary matrix. J. Soc. Indust. Appl. Math. 10, 74–88 (1962)
Gao, B., Liu, X., Chen, X., Yuan, Y.-X.: A new first-order algorithmic framework for optimization problems with orthogonality constraints. SIAM J. Optim. 28(1), 302–332 (2018)
Golub, G.H., Van Loan, C.F.: Matrix Computations. Johns Hopkins Studies in the Mathematical Sciences, 4th edn. Johns Hopkins University Press, Baltimore (2013)
Greenstadt, J.: A method for finding roots of arbitrary matrices. Math. Tables Aids Comput. 9, 47–52 (1955)
Higham, N. J.: Accuracy and Stability of Numerical Algorithms, 2nd edn. SIAM, Philadelphia (2002)
Higham, N. J.: Functions of Matrices. Society for Industrial and Applied Mathematics. SIAM, Philadelphia (2008)
Higham, N.J., Schreiber, R.S.: Fast polar decomposition of an arbitrary matrix. SIAM J. Sci Statist. Comput. 11(4), 648–655 (1990)
Hüper, K., Van Dooren, P.: New algorithms for the iterative refinement of estimates of invariant subspaces. Journal Future Generation Computer Systems 19, 1231–1242 (2003)
Ipsen, I. C. F.: Computing an eigenvector with inverse iteration. SIAM Rev. 39(2), 254–291 (1997)
Jahn, H. A.: Improvement of an approximate set of latent roots and modal columns of a matrix by methods akin to those of classical perturbation theory. Quart. J. Mech. Appl. Math. 1, 131–144 (1948)
Jonsson, I., Kågström, B.: Recursive blocked algorithm for solving triangular systems. I. one-sided and coupled Sylvester-type matrix equations. ACM Trans. Math Software 28(4), 392–415 (2002)
Jonsson, I., Kågström, B.: Recursive blocked algorithm for solving triangular systems. II. Two-sided and generalized Sylvester and Lyapunov matrix equations. ACM Trans. Math Software 28(4), 416–435 (2002)
Kahan, W.: Refineig: a program to refine eigensystems, 2007. Available from https://people.eecs.berkeley.edu/wkahan/Math128/RefinEig.pdf
Konstantinov, M. M., Petkov, P. Hr., Christov, N. D.: Nonlocal perturbation analysis of the Schur system of a matrix. SIAM J Matrix Anal. Appl. 15(2), 383–392 (1994)
Lange, M., Rump, S.M.: An alternative approach to Ozaki’s scheme for error-free transformation of matrix multiplication, 2022 Presentation at International Workshop on Reliable Computing and Computer-Assisted Proofs (reCAP2022)
Lundström, E.: Adaptive eigenvalue computations using Newton’s method on the Grassmann manifold. SIAM J Matrix Anal. Appl. 23(3), 819–839 (2002)
Mehl, C.: On asymptotic convergence of nonsymmetric Jacobi algorithms. SIAM J. Matrix Anal. Appl. 30(1), 291–311 (2008)
Nocedal, J., Wright, S. J.: Numerical optimization. Springer Series in Operations Research and Financial Engineering, 2nd edn. Springer, New York (2006)
Ogita, T., Aishima, K.: Iterative refinement for symmetric eigenvalue decomposition. Jpn. J. Ind. Appl Math. 35, 1007–1035 (2018)
Ogita, T., Aishima, K.: Iterative refinement for symmetric eigenvalue decomposition II: clustered eigenvalues. Jpn. J. Ind. Appl. Math. 36(2), 435–459 (2019)
O’Leary, D.P., Stewart, G. W.: Data-flow algorithms for parallel matrix computations. Comm. ACM 28(8), 840–853 (1985)
Ozaki, K., Ogita, T., Oishi, S., Rump, S. M.: Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numerical Algorithms 59(1), 95–118 (2012)
Stewart, G. W.: A Jacobi-like algorithm for computing the Schur decomposition of a non-Hermitian matrix. SIAM J. Sci. Statist. Comput. 6(4), 853–864 (1985)
Stewart, G. W., Sun, J. -G.: Matrix Perturbation Theory. Academic Press, New York (1990)
Sun, J.-G.: Perturbation Bounds for the Schur Decomposition Report UMINF-92.20, Department of Computing Science, Umeå University, Umeå, Sweden (1992)
Sun, J. -G.: Perturbation bounds for the generalized Schur decomposition. SIAM J. Matrix Anal. Appl. 16(4), 1328–1340 (1995)
Wielandt, H.: Beiträge Zur Mathematischen Behandlung Komplexer Eigenwertprobleme, Teil V: Bestimmung Höherer Eigenwerte Durch Gebrochene Iteration. Bericht B 44/J/37. Aerodynamische Versuchsanstalt Göttingen, Germany (1944)
Wilkinson, J.H.: The Algebraic Eigenvalue Problem. Clarendon Press, Oxford (1965)
Acknowledgements
The authors thank Takeshi Ogita for providing the Matlab toolbox acc based on [37]. They are also grateful to Nicolas Boumal and Christian Lubich for discussions related to Remark 1.
Funding
The first author was supported by the Croatian Science Foundation under grant IP-2019-04-6268.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare no competing interests.
Additional information
Availability of data and materials
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Zvonimir Bujanović, Daniel Kressner, and Christian Schröder contributed equally to this work.
Rights and permissions
About this article
Cite this article
Bujanović, Z., Kressner, D. & Schröder, C. Iterative refinement of Schur decompositions. Numer Algor 92, 247–267 (2023). https://doi.org/10.1007/s11075-022-01327-6
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11075-022-01327-6