Simple exponential acceleration of the power iteration algorithm

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

A \in C^{n \times n}

, we may write

A = U Λ U^{- 1}

(1)

A^{k} = (U Λ U^{- 1}) (U Λ U^{- 1}) \dots (U Λ U^{- 1}) = U Λ^{k} U^{- 1}

(2)

where the columns of $U$ are the eigenvectors $v_{i}$ , and $Λ$ is diagonal with the eigenvalues $λ_{i}$ of $A$ along the diagonal. Multiplying Equation 2 on the right by $U$ ,

A^{k} U = U Λ^{k}

(3)

Randomly choosing a vector $x = U y \in C^{n}$ , we may write

A^{k} x = U Λ^{k} y = λ_{1}^{k} (n \sum i = 1 v_{i} {(\frac{λ_{i}}{λ_{1}})}^{k} y_{i})

(4)

Here, we assume that $λ_{1}$ is the largest eigenvalue of $A$ , and $v_{1}$ is its associated eigenvector. With each factor of $A$ , the components along the other eigenvectors are suppressed by a factor of $| \frac{λ_{i}}{λ_{1}} | < 1$ , and thus we expect $A^{k} x$ to become increasingly parallel to $v_{1}$ .

The power iteration algorithm is usually expressed as below, where $x_{0}$ is drawn randomly and $f$ indicates the final value of $k$ :

Input :

A \in R^{n \times n}

, tolerance

t o l

x_{0} \in R^{n}

k \leftarrow 0

3 do

x_{k + 1} \leftarrow (A x_{k}) / (| A x_{k z} |)

k \leftarrow k + 1

6while $| x_{k + 1} - x_{k} | > t o l$ ;

7return

(x_{f}^{†} A x_{f}) / (x_{f}^{†} x_{f})

and

x_{f}

Output :

Rayleigh quotient estimate of the largest eigenvalue

λ

, estimated eigenvector

x

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 $k = 2^{j}, j \in N$ . $f$ again indicates the final value of $i$ .

Input :

A \in R^{n \times n}

, desired power

k = 2^{j}

A_{0} \leftarrow A

i \leftarrow 0

3 do

A_{i + 1} \leftarrow A_{i} \cdot A_{i}

i \leftarrow i + 1

6while $i \leq j$ ;

return

A_{i}

Output :

A_{f} = A^{k}

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 :

A \in R^{n \times n}

, tolerance

t o l

x_{0} \in R^{n}

A_{0} \leftarrow A

i \leftarrow 0

4 do

A_{i + 1} \leftarrow A_{i} \cdot A_{i}

A_{i + 1} \leftarrow A_{i + 1} / | A_{i + 1} |

i \leftarrow i + 1

8while $| A_{i + 1} - A_{i} | > t o l$ ;

x_{f} \leftarrow A_{f} x_{0}

10 return

(x_{f}^{†} A_{0} x_{f}) / (x_{f}^{†} x_{f})

and

x_{f}

Output : Rayleigh quotient estimate of the largest eigenvalue

λ

, estimated eigenvector

x

Algorithm 3 Power iteration with exponentiation

To generalize Algorithms 1 and 3 to the problem of finding the largest $k$ 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 $A$ , as shown in Algorithm 4.

Input :

A \in R^{n \times n}

, desired number of largest eigenvalues

k

A_{0} \leftarrow A

i \leftarrow 0

B \leftarrow {e_{i}}

Λ \leftarrow []

5 do

λ, v \leftarrow PowerIteration (A_{i})

7 Append

λ

Λ

e_{i} \leftarrow v / | v |

j \leftarrow 0

A_{i + 1} \leftarrow A_{i} - λ e_{i} e_{i}^{†}

11while $i \leq k$ ;

12return

B

Λ

Output : The largest

k

eigenvectors and eigenvalues of

A

Algorithm 4 Gram-Schmidt Orthonormalization

To demonstrate Algorithm 4, assume we have computed the largest eigenvalue and eigenvector of $A$ 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 $A$ in line 10

(A - λ_{1} v_{1} v_{1}^{†}) v_{1} = λ_{1} v_{1} - λ_{1} v_{1} (v_{1}^{†} v_{1}) = 0

(5)

Thus $v_{1}$ now has eigenvalue $0$ for the matrix $A - λ_{1} v_{1} v_{1}^{T}$ . Because we have the additional assumption that $A$ is self-adjoint ( $A = A^{†}$ ), 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 $10^{- 10}$ . 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 2¹¹1Table 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 $10^{6}$

iterations to reach the desired tolerance, which would occasionally significantly skewed the timings (

n = 100

, 300 matrices taking almost 60 s to complete). In contrast, Algorithm 3 consistently reached convergence within 20 iterations, shown in Figure 1.

	$n = 100$ (300 matrices)	$n = 1000$ (5 matrices)	$n = 3000$ (1 matrix)	$n = 5000$ (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 $\times$	36 $\times$	21 $\times$	11 $\times$

Table 1: Benchmarking results (real matrices)

	$n = 100$ (300 matrices)	$n = 1000$ (5 matrices)	$n = 3000$ (1 matrix)	$n = 5000$ (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 $\times$	9.5 $\times$	7 $\times$	3.8 $\times$

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 $\times$ speed up on three hundred $100 \times 100$ 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 $A$ , which is used to project $x_{0}$ to $x_{f}$ , which as mentioned in the introduction should be close to $v_{1}$ . Crucially, although naive multiplication of $n \times n$ matrices in line 5 of Algorithm 3 takes $O (n^{3})$ time, the powers ${(\frac{λ_{i}}{λ_{1}})}^{k}$ take $O (lg k)$ steps to compute and thus there is an exponential increase in accuracy with each matrix multiplication performed. Overall, Algorithm 3 takes $O (n^{3} lg k)$ time, using naive matrix multiplication. In contrast, Algorithm 1 takes $O (k)$ steps to compute ${(\frac{λ_{i}}{λ_{1}})}^{k}$ , and thus $O (n^{2} k)$ 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 ( $n = 100$ ).

Just as in the traditional power iteration algorithm, Algorithm 3 relies on the choice of an appropriate $x_{0}$ . Random choice of $x_{0}$ should prevent the degenerate case of $x_{0}$ orthogonal to $v_{1}$ except in rare cases. Algorithm 3 may also be more efficient with sparse matrices, as the cost of $O (n^{3})$ 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 $k$ largest eigenvalues and associated eigenvectors in Algorithm 4 when applied to self-adjoint matrices, which we observed to work for small values of $k$ (see provided code). Algorithm 4²²2Toy code is shown in —eigenbasis.py—. may be able to compete with other iterative methods for small $k$ and large $n$ , 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 ( $n = 100$ , 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 ${log}_{2}$ 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.

Simple exponential acceleration of the power iteration algorithm

Authors

A sparse approximate inverse for triangular matrices based on Jacobi iteration

Multiplication-Avoiding Variant of Power Iteration with Applications

Dynamical perturbation theory for eigenvalue problems

A fast iterative algorithm for near-diagonal eigenvalue problems

A simple extrapolation method for clustered eigenvalues

A FEAST variant incorporated with a power iteration

Fast Coherent Point Drift

1 Introduction

2 Material and Methods

2.1 Proposed algorithms

2.2 Implementation and computational resources

3 Results

4 Discussion

5 Conclusions

6 Competing interests

7 Funding

References

Simple exponential acceleration of the power iteration algorithm

Authors

Related Research

A sparse approximate inverse for triangular matrices based on Jacobi iteration

Multiplication-Avoiding Variant of Power Iteration with Applications

Dynamical perturbation theory for eigenvalue problems

A fast iterative algorithm for near-diagonal eigenvalue problems

A simple extrapolation method for clustered eigenvalues

A FEAST variant incorporated with a power iteration

Fast Coherent Point Drift

This week in AI

1 Introduction

2 Material and Methods

2.1 Proposed algorithms

2.2 Implementation and computational resources

3 Results

4 Discussion

5 Conclusions

6 Competing interests

7 Funding

References