Fast methods for resumming matrix polynomials and Chebyshev matrix polynomials

Abstract

Fast and effective algorithms are discussed for resumming matrix polynomials and Chebyshev matrix polynomials. These algorithms lead to a significant speed-up in computer time by reducing the number of matrix multiplications required to roughly twice the square root of the degree of the polynomial. A few numerical tests are presented, showing that evaluation of matrix functions via polynomial expansions can be preferable when the matrix is sparse and these fast resummation algorithms are employed.

Introduction

Functions of a matrix can be defined by their series expansions, such as the matrix exponential:

$$\text{e}{}^{\mathrm{<mtext>X</mtext>}}\text{=1+}\text{X}\text{+}\text{1}\text{2}\text{X}{}^2\text{+}\text{1}\text{3!}\text{X}{}^3\text{+⋯.}$$

It is relatively unusual for a matrix function to be evaluated this way because typically it is more efficient and precise to employ the eigenvectors and eigenvalues of X, if it is diagonalizable. For example, if X is Hermitian, i.e., $\text{X}{}^{\mathrm{<mtext>T</mtext>}}\text{=}\text{X}$, we can evaluate the matrix function $\text{e}{}^{\mathrm{<mtext>X</mtext>}}$ as

$$\text{e}{}^{\mathrm{<mtext>X</mtext>}}\text{=}\text{U}\text{e}{}^{\mathrm{<mtext>x</mtext>}}\text{U}{}^{\mathrm{†}}\text{,}$$

where U are the eigenvectors and x the diagonal matrix of eigenvalues of matrix X.

However, some situations arise where evaluation of the series is actually preferable. The particular case we have in mind is where the dimension of X is very large, and in addition (or as a consequence) X is very sparse. Then sparse matrix multiplication methods can be used to evaluate the products entering the matrix polynomial. The resulting computational effort can then in principle increase only linearly with the dimension of the matrix, if its bandwidth is constant. By contrast, explicitly obtaining the eigenvalues and eigenvectors cannot scale linearly – indeed it is cubic scaling if standard linear algebra routines are employed. For example, this situation arises in formulating linear scaling algorithms to perform tight-binding or Kohn–Sham density functional theory calculations on very large molecules [1], [2], [3], [4], [5], [6], [7], in evaluation of time-evolution operator e^−iHt [8], and the calculation of Boltzmann factors and partition functions of very large molecules in computational physics [9], [10], [11].

In this paper, we report fast algorithms for summing matrix and Chebyshev matrix polynomials. The problem is to minimize the computer time required to perform the sum to a given power, N_pol, the degree of the polynomial. To an excellent approximation this is equivalent to minimizing the number of matrix multiplications, since the cost of matrix multiplications completely dominates the cost of matrix additions. The other degree of freedom that is available is to store multiple intermediate quantities to permit reduction of the matrix multiplication effort.

Simple term-by-term evaluation of a matrix polynomial of degree N_pol requires N_pol−1 matrix multiplications. However it is possible to do substantially better, as has been shown in a number of papers where fast algorithms have been reported for resumming matrix polynomials with fewer matrix multiplications [12], [13], [14], [15], [16]. Perhaps the most effective algorithm was suggested by Paterson and Stockmeyer [14], which requires only $\text{O}\text{(}\text{N}{}_{\mathrm{<mtext>pol</mtext>}}\text{)}$ nonscalar multiplications to evaluate polynomials of degree N_pol. These authors also presented proofs that at least $\text{N}{}_{\mathrm{<mtext>pol</mtext>}}$ nonscalar multiplications are required to resum a polynomial of degree N_pol. The algorithm has practical application in the evaluation of matrix polynomials with scalar coefficients and offers huge improvement over simple term-by-term evaluation, at the expense of requiring more stored matrices. Storage of $\text{O}\text{(}\text{N}{}_{\mathrm{<mtext>pol</mtext>}}\text{)}$ matrices is required in this algorithm, where N is the dimension of the matrix X. Van Loan [15] later showed how this procedure could be implemented with greatly reduced storage, although significantly more matrix multiplies were required.

In the following section (Section 2), we review the algorithm of Paterson and Stockmeyer [14] and identify other two algorithms for matrix polynomial evaluation that also scale with $\text{N}{}_{\mathrm{<mtext>pol</mtext>}}$, the square root of the maximum degree of the polynomial. The coefficient is approximately 2. These two relatively distinct algorithms turn out to yield very similar matrix multiplication counts but require fewer stored matrices than Paterson and Stockmeyer’s algorithm.

Another purpose of the paper is to generalize these algorithms for summing simple matrix polynomials to accelerate the summation matrix series based on Chebyshev matrix polynomials of the form

$\text{T}{}_{\mathrm{i}}\text{(}\text{X}\text{)}$ is the Chebyshev matrix polynomial of degree i, which is recursively defined in terms of T_i−1 and T_i−2, and a₀,a₁,…,a_N are the scalar coefficients. If the eigenvalues of X lie in the range [−1,1], then Chebyshev approximations are numerically very stable and they are less susceptible to round-off errors due to limited machine precision than the equivalent power series [17]. In the context of linear scaling electronic structure methods, Chebyshev matrix polynomials have therefore been proposed as a more stable and efficient alternative to simple matrix polynomials [3], [18], [19]. These generalizations are described in Section 3. They yield algorithms that are broadly similar to the ones described for matrix polynomials, with differences arising from the recursion relations associated with Chebyshev matrix polynomials. Detailed application of this approach to linear scaling electronic structure calculations based on Chebyshev polynomials is described elsewhere [20].

In Section 4, we present three types of data to characterize the performance of the 3 fast summation algorithms. First, the number of matrix multiplies, as well as the number of stored intermediates, are reported for polynomials of various degrees. Second, computer time for the evaluation of fast summations is compared with the time for conventional term-by-term summation on some test matrices. Third, the computer time for evaluating matrix polynomials by fast summation is compared with evaluation by explicit diagonalization for the case where the matrix is large and sparse. Finally we conclude in Section 5.