Simple exponential acceleration of the power iteration algorithm

09/22/2021
by   Congzhou M Sha, et al.
Penn State University
0

Many real-world problems rely on finding eigenvalues and eigenvectors of a matrix. The power iteration algorithm is a simple method for determining the largest eigenvalue and associated eigenvector of a general matrix. This algorithm relies on the idea that repeated multiplication of a randomly chosen vector x by the matrix A gradually amplifies the component of the vector along the eigenvector of the largest eigenvalue of A while suppressing all other components. Unfortunately, the power iteration algorithm may demonstrate slow convergence. In this report, we demonstrate an exponential speed up in convergence of the power iteration algorithm with only a polynomial increase in computation by taking advantage of the commutativity of matrix multiplication.

READ FULL TEXT VIEW PDF

Authors

page 6

06/24/2021

A sparse approximate inverse for triangular matrices based on Jacobi iteration

In this paper, we propose a sparse approximate inverse for triangular ma...
10/22/2021

Multiplication-Avoiding Variant of Power Iteration with Applications

Power iteration is a fundamental algorithm in data analysis. It extracts...
02/28/2020

Dynamical perturbation theory for eigenvalue problems

Many problems in physics, chemistry and other fields are perturbative in...
12/29/2020

A fast iterative algorithm for near-diagonal eigenvalue problems

We introduce a novel iterative eigenvalue algorithm for near-diagonal ma...
06/17/2020

A simple extrapolation method for clustered eigenvalues

This paper introduces a simple variant of the power method. It is shown ...
12/03/2019

A FEAST variant incorporated with a power iteration

We present a variant of the FEAST matrix eigensolver for solving restric...
06/11/2020

Fast Coherent Point Drift

Nonrigid point set registration is widely applied in the tasks of comput...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Eigendecomposition of matrices is useful in continuum mechanics [1], quantum mechanics [2]

, machine learning

[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 :

Input : , tolerance
1
2
3 do
4      
5      
6while ;
7return and
Output :

 Rayleigh quotient estimate of the largest eigenvalue

, estimated eigenvector
Algorithm 1 Power Iteration

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 .

Input : , desired power
1
2
3 do
4      
5      
6while ;
return Output : 
Algorithm 2 Matrix exponentiation through squaring

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:

Input : , tolerance
1
2
3
4 do
5      
6      
7      
8while ;
9
10 return and
Output : Rayleigh quotient estimate of the largest eigenvalue , estimated eigenvector
Algorithm 3 Power iteration with exponentiation

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.

Input : , desired number of largest eigenvalues
1
2
3
4
5 do
6      
7       Append to
8      
9      
10      
11while ;
12return ,
Output : The largest eigenvectors and eigenvalues of
Algorithm 4 Gram-Schmidt Orthonormalization

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
Table 1: Benchmarking results (real matrices)
(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
Table 2: Benchmarking results (complex matrices)

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].

Figure 1: Histogram of number of iterations needed to reach convergence (, 300 matrices). A: real matrices, Algorithm 1. B: real matrices, Algorithm 3. C: complex matrices, Algorithm 1. D: complex matrices, Algorithm 3. The distribution of iterations for convergence in Algorithm 1 is similar to linear scale for iterations of Algorithm 3. Algorithm 3 reaches convergence within 20 iterations due to the exponential growth of its accuracy. There is long tail to the right in the execution times of both algorithms, which is greatly exaggerated when performing naive power iteration, leading to slow and unpredictable convergence when using Algorithm 1.

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.