1 Introduction
Eigendecomposition of matrices is useful in continuum mechanics [1], quantum mechanics [2]
[3], graph theory [4], and many other areas of science and mathematics.The power iteration algorithm is a well-known method for determining the largest eigenvalue and its associated eigenvector of a general matrix [5]. The following argument is taken from [6]
. Given a diagonalizable matrix
, we may write(1) |
(2) |
where the columns of are the eigenvectors , and is diagonal with the eigenvalues of along the diagonal. Multiplying Equation 2 on the right by ,
(3) |
Randomly choosing a vector , we may write
(4) |
Here, we assume that is the largest eigenvalue of , and is its associated eigenvector. With each factor of , the components along the other eigenvectors are suppressed by a factor of , and thus we expect to become increasingly parallel to .
The power iteration algorithm is usually expressed as below, where is drawn randomly and indicates the final value of :
However, there is a well-known algorithm for the computation of high powers of a matrix: exponentiation through squaring. Without loss of generality, we assume that . again indicates the final value of .
2 Material and Methods
2.1 Proposed algorithms
We propose an algorithm which combines Algorithms 1 and 2 to achieve exponential convergence of the largest eigenvalue:
To generalize Algorithms 1 and 3 to the problem of finding the largest eigenvectors in the case of a self-adjoint (Hermitian) matrix, we may perform a version of Gram-Schmidt orthonormalization carried out on the input matrix , as shown in Algorithm 4.
To demonstrate Algorithm 4, assume we have computed the largest eigenvalue and eigenvector of in line 6 of Algorithm 4 (by using either Algorithm 1 or 3). After normalizing the eigenvector (line 8), Then we can remove that component from the matrix in line 10
(5) |
Thus now has eigenvalue for the matrix . Because we have the additional assumption that is self-adjoint (), its eigendecomposition is guaranteed to be unitary, by the finite-dimensional spectral theorem [7]
. Therefore, the Gram-Schmidt process works since at each step, we do not change the span of the remaining eigenvectors.
2.2 Implementation and computational resources
We implemented these algorithms and visualized results in Python 3.8 [8], used the NumPy 1.19 package for linear algebra as reference [9], and created figures with Matplotlib 3.4 [10]. The code is included in the supplemental data (*.py
). Matrices of various dimensions were generated randomly with entries drawn from the normal distribution (mean 0, standard deviation 1) and symmetrized, to guarantee real eigenvalues and eigenvectors in the case of real matrices. The largest eigenvalue of each matrix was computed using the NumPy linear algebra library (
numpy.linalg.eigvals
). The matrix norms in Algorithms 1 line 4 and Algorithm 3 line 6 are necessary to prevent numerical overflow, but are otherwise arbitrary, so we chose the efficient max norm (maximum matrix element).
We recorded the computation time and the number of iterations of Algorithms 1 and 3 required to reach a numerical tolerance of . Calculations were performed on 12 cores of an AMD Ryzen 9 3900X CPU, with 64 GB of RAM available, allowing for parallelized matrix-matrix and matrix-vector multiplication through the NumPy library.
3 Results
We observed significant speedups in using Algorithm 3 as compared to Algorithm 1 for both real and complex matrices, shown in Tables 1 and 2111Table 1 corresponds to the file —test.py— and Table 2 to *—testcomplex.py—.. We also observed that it took far fewer iterations to reach numerical convergence using Algorithm 3 compared to Algorithm 1. We observed several instances in which it took Algorithm 1 over
iterations to reach the desired tolerance, which would occasionally significantly skewed the timings (
, 300 matrices taking almost 60 s to complete). In contrast, Algorithm 3 consistently reached convergence within 20 iterations, shown in Figure 1.(300 matrices) | (5 matrices) | (1 matrix) | (1 matrix) | |
---|---|---|---|---|
Alg 1 | 35 s (117 ms/matrix) | 29 s (5.8 s/matrix) | 59 s | 110 s |
Alg 3 | 0.54 s (1.8 ms/matrix) | 0.77 s (153 ms/matrix) | 2.8 s | 9.7 s |
speedup | 65 | 36 | 21 | 11 |
(300 matrices) | (5 matrices) | (1 matrix) | (1 matrix) | |
---|---|---|---|---|
Alg 1 | 67 s (223 ms/matrix) | 69 s (14 s/matrix) | 175 s | 359 s |
Alg 3 | 1.38 s (4.6 ms/matrix) | 7.3 s (1.5 s/matrix) | 25 s | 95 s |
speedup | 49 | 9.5 | 7 | 3.8 |
Even though Algorithm 3 was written in Python whereas NumPy code executes highly optimized LAPACK routines for eigendecomposition, we also observed that Algorithm 3 (taken advantage of LAPACK matrix multiplication) is faster than NumPy’s numpy.linalg.eigvals
when we desire only the largest eigenvalue (3 speed up on three hundred matrices). It may be beneficial to write a version of Algorithm 3 in a lower level language like C or Fortran to perform a more direct comparison.
4 Discussion
Algorithm 3 computes a high (normalized) power of , which is used to project to , which as mentioned in the introduction should be close to . Crucially, although naive multiplication of matrices in line 5 of Algorithm 3 takes time, the powers take steps to compute and thus there is an exponential increase in accuracy with each matrix multiplication performed. Overall, Algorithm 3 takes time, using naive matrix multiplication. In contrast, Algorithm 1 takes steps to compute , and thus time. Thus there is a tradeoff between the cost of matrix multiplication and desired accuracy, in which Algorithm 3 appears to be more performant in moderately-sized practical computations ().
Just as in the traditional power iteration algorithm, Algorithm 3 relies on the choice of an appropriate . Random choice of should prevent the degenerate case of orthogonal to except in rare cases. Algorithm 3 may also be more efficient with sparse matrices, as the cost of matrix-matrix multiplication may be greatly reduced, though further benchmarking is required to compare with sparse applications of Algorithm 3. For other applications, the throughput of GPU-accelerated matrix-matrix multiplication may further equalize the performance of matrix-matrix multiplication in Algorithm 3 compared to the matrix-vector multiplication in Algorithm 1.
The major drawback of Algorithm 3 is the same as in traditional power iteration: only the largest eigenvalue and associated eigenvector are computed. Other iterative methods such as the Lanczos algorithm [11, 12] and inverse iteration [13] may compute other eigenvalues. It may also be possible to incorporate Algorithm 2 into these iterative methods to speed up convergence. We also provided a method for computing the largest eigenvalues and associated eigenvectors in Algorithm 4 when applied to self-adjoint matrices, which we observed to work for small values of (see provided code). Algorithm 4222Toy code is shown in —eigenbasis.py—. may be able to compete with other iterative methods for small and large , though benchmarking is needed.
Even with its drawbacks, Algorithm 3 may have practical uses. The largest eigenvalue and associated eigenvector can be highly informative, for example in the graph similarity score (Equation 5) of Zager and Verghese [14]. There may be other specific use cases where Algorithm 3 is preferred over the full eigendecomposition.
To our knowledge, our proposed modification to power iteration has not been reported in the literature. Work by Shi and Petrovic to improve power iteration convergence relies on incorporating information about the second eigenvalue, but does not demonstrate the exponential speed of our Algorithm 3 [15].

5 Conclusions
Modifying the power iteration algorithm to take advantage of matrix exponentiation by squaring produces a fast, practical, and robust method of determining the largest eigenvalue and its associated eigenvector for matrices. Our Algorithm 3 outperforms traditional power iteration in convergence rate and speed on a variety of matrix sizes. We recommend Algorithm 3 as a drop-in replacement for power iteration due to its favorable performance and simple implementation.
6 Competing interests
The authors have no competing interests to disclose.
7 Funding
We acknowledge support from the National Institutes for Health (R35 GM134864), National Science Foundation (2040667) and the Passan Foundation. This project is also funded, in part, under a grant with the Pennsylvania Department of Health using Tobacco CURE Funds. The Department specifically disclaims responsibility for any analyses, interpretations or conclusions.
References
- [1] L. Landau and E. Lifshitz, Theory of Elasticity (Third Edition). Oxford: Butterworth-Heinemann, third edition ed., 1986.
- [2] J. J. Sakurai, Modern quantum mechanics; rev. ed. Reading, MA: Addison-Wesley, 1994.
-
[3]
I. T. Jolliffe and J. Cadima, “Principal component analysis: a review and recent developments,”
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 374, no. 2065, p. 20150202, 2016. - [4] F. R. K. Chung, Spectral Graph Theory. Providence, RI: American Mathematical Society, 1997.
- [5] R. Mises and H. Pollaczek-Geiringer, “Praktische verfahren der gleichungsauflösung,” Journal of Applied Mathematics and Mechanics, vol. 9, no. 2, pp. 152–164, 1929.
- [6] D. Bindel, “Lecture notes for matrix computations,” October 2016.
- [7] M. Reed and B. Simon, eds., Methods of Modern Mathematical Physics, vol. 1. Academic Press, 1972.
- [8] G. Van Rossum and F. L. Drake, Python 3 Reference Manual. Scotts Valley, CA: CreateSpace, 2009.
- [9] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature, vol. 585, pp. 357–362, Sept. 2020.
- [10] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering, vol. 9, no. 3, pp. 90–95, 2007.
- [11] C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” J. Res. Natl. Bur. Stand. B, vol. 45, pp. 255–282, 1950.
- [12] I. U. Ojalvo and M. Newman, “Vibration modes of large structures by an automatic matrix-reductionmethod,” AIAA Journal, vol. 8, no. 7, pp. 1234–1239, 1970.
- [13] E. Pohlhausen, “Berechnung der eigenschwingungen statisch-bestimmter fachwerke,” ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, vol. 1, no. 1, pp. 28–42, 1921.
- [14] L. A. Zager and G. C. Verghese, “Graph similarity scoring and matching,” Applied Mathematics Letters, vol. 21, no. 1, pp. 86–94, 2008.
- [15] B. Shi and B. Petrovic, “Implementation of the modified power iteration method to two-group monte carlo eigenvalue problems,” Annals of Nuclear Energy, vol. 38, no. 4, pp. 781–787, 2011.