Single Precision Calculation of Iterative Refinement of Pairs of aReal Symmetric-Definite Generalized Eigenproblem by Using a
Filter Composed of a Single Resolvent
Hiroshi Murakami

Tokyo Metropolitan Universitymurakami@tmu.ac.jp

KEYWORDS
Eigenproblem, Filter, Iterative Refinement, Single Precision, Resolvent
ACM Reference Format:
Hiroshi Murakami. 2019. Single Precision Calculation of Iterative Refinement of Pairs of a Real Symmetric-Definite Generalized Eigenproblem byUsing a Filter Composed of a Single Resolvent. In Proceedings of 48th International Conference on Parallel Processing (ICPP2019). ACM, New York, NY,USA, 2 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn

1INTRODUCTION

For a given generalized eigenproblem Av = λBv whose matrices
A and B are real symmetric and B is positive definite, we solve all
those approximate pairs whose eigenvalues are in the specified
real interval [a, b] by using a filter. There are many references ofmethods which can be classified as filter diagonalization methods[1, 2, 4–7, 10].
The filter is usually composed of several resolvents R(ρ i ) ≡(A − ρ i B)−1 B whose shift ρ i are complex numbers.
For a given vector x, an application of the resolvent y ← R(ρ)xis to solve a system of linear equations C(ρ)y = Bx for y, hereC(ρ) ≡ A − ρB is the shifted matrix corresponds to both matricesof the generalized eigenproblem. In this study, for the solution ofthis kind of system, we assume to use some direct method whichuses matrix factorization.
When the shift ρ is a real number, the matrix C(ρ) is real symmetric. When the shift is a real number less than the minimum
eigenvalue λ min of the eigenproblem, the matrix is real symmetricpositive definite. When the shift is an imaginary number, the matrix
is complex symmetric and non-singular. For a symmetric matrix
either real or complex, the modified Cholesky method can be used
to solve the system of linear equations by a matrix decomposition
and forward and backward substitutions (The modified Cholesky
method for complex symmetric matrices is derived from the method
for real symmetric matrices by replacing numbers and arithmetic
expressions from real to complex).
Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACM
must be honored. Abstracting with credit is permitted. To copy otherwise, or republish,
to post on servers or to redistribute to lists, requires prior specific permission and/or afee. Request permissions from permissions@acm.org.
ICPP2019, Aug 5–8, Kyoto Research Park, Kyoto, Japan© 2019 Association for Computing Machinery.
ACM ISBN 978-x-xxxx-xxxx-x/YY/MM. . . $15.00
https://doi.org/10.1145/nnnnnnn.nnnnnnn

When computing resources are limited but the problem size islarge, the filtering calculation tends to be restricted by both amountsof computation to factor shifted matrices and especially storage tohold factors of matrices. Both are proportional to the number ofresolvents to construct the filter. So it is desirable to reduce the number of resolvents. There are two kinds of filters composed of only asingle resolvent: 1) The filter which is a real polynomial of a resolvent with a real shift, 2) The filter which is a real polynomial of theimaginary part of a resolvent with an imaginary shift. In this study,in order to make the filter design simple, a Chebyshev polynomial
is used to express the "real polynomial". When the interval is located at the lower-end of the eigenvalue distribution, we use a filter
F = дs Tn (2 γ R(ρ) − I ) which is an n-th degree Chebyshev polynomial of a single resolvent R(ρ) whose shift ρ is real and less thanthe minimum eigenvalue. When the interval is in the middle of theeigenvalue distribution, we use a filter F = дs Tn (2 γ ′ Im R(ρ ′ ) − I )which is an n-th degree Chebyshev polynomial of the imaginary
part of a single resolvent R(ρ ′ ) whose shift ρ ′ is imaginary. Here,дs is the tight upper-bound of the transfer function magnitude of
the filter in the stop-band, γ and γ ′ are real constants, and I is theidentity operator.
However, properties of these simple filters are not very well,because they are composed of a single resolvent rather than many,
and also their real polynomials are expressed by just using Chebyshev polynomials. For example, their transfer functions cannot
have steep changes of values, thus µ − 1 the geometrical ratio of
the width of transition-bands to the width of the pass-band cannot
be made very small. Also, if the value of дs is set very small, whichis the upper-bound of the transfer-rates in stop-bands, the valueof 1/дp will be large, which is the max-min ratio of the transfer
function of the filter in the passband λ ∈ [a, b]. When this max-min
ratio is very large, the contained rates of required eigenvectors inthe set of vectors after a filtration tend to have different orders ofmagnitudes. Therefore, within a vector, suppressed by large values of those eigenvectors whose transfer-rates are larger, smaller
values of those eigenvectors whose transfer-rates are smaller loseaccuracy by rounding errors. By this reason, those eigenvectors
whose transfer-rates are smaller, which are extracted from a setof filtered vectors, tend to have lower accuracy. Therefore, someapproximate pairs may not attain the level of required accuracy.
In the above explanation about the filtering method, we assumed
to apply the filter only once to a set of random initial vectors. The
following procedure shows how to calculate approximate pairs witha single application of the filter.

ICPP2019, Aug 5–8, Kyoto Research Park, Kyoto, Japan
1) Let Y (0) be an initial set of m random column vectors.
2) B-orthonormalize Y (0) to make X (1) ;X (1) is filtered to make Y (1) .
3) Considering properties of the filter, approximate pairs areconstructed from both sets of vectors X (1) and Y (1) .

2

ITERATIVE REFINEMENT OF EIGENPAIRSBY USING A FILTER

When properties of the filter are not good, approximate pairs arealso not good in accuracy which are obtained by an applicationof the filter to a set of initial random vectors. Even in such case,
we can improve approximate pairs with a few applications of thecombination of B-orthonormalization and filtering.
The following procedure shows how to calculate approximate
pairs with applications of the filter IT times.
1) Let Y (0) be an initial set of m random column vectors.
2) Iterate the followings for i = 1, . . ., IT
B-orthonormalize Y (i−1) to make X (i) ;X (i) is filtered to make Y (i) .
3) Considering properties of the filter, approximate pairs areconstructed from both sets of vectors X (IT) and Y (IT) .
(During the iteration, in B-orthonormalization in the above step2, if the effective rank of the set of vectors is found decreased, wedecrease the number m of orthonormalized vectors in the set.)
The orthonormalization prevents eigenvectors of small transferrates from losing their accuracy by numerical rounding errors. It
prevents the set of vectors from being occupied by those eigenvectors whose transfer-rates are larger. The principle to use orthogonalization of vectors in each iteration step is called "orthogonaliteration" and is well known [3, 8, 9].

3CONCLUSION

We made some experiments to solve pairs of a real symmetricdefinite generalized eigenproblem whose eigenvalues are in thespecified interval by using a filter.
In this study we used filters composed of an action of a singleresolvent. We used a real shift for the resolvent to solve pairs whoseeigenvalues are lowest. When we used an imaginary shift, theinterval for eigenvalues may be placed anywhere. The filter weused is a real polynomial of a single resolvent whose shift is real, ora real polynomial of an imaginary part of a single resolvent whoseshift is imaginary. When the degree of the real polynomial is n,
in an application of the filter the resolvent is applied n times. AChebyshev polynomial is used to represent the real polynomial to
make the filter design simple, and an application of the filter canbe calculated by using the three-term recursion.
An application of a resolvent to a vector is to solve a systemof linear equations whose coefficient is the shifted matrix madefrom both matrices of the generalized eigenproblem. In our study,
the system of linear equations is assumed to be solved by somedirect method using decomposition of the coefficient matrix. Since
we use a filter consists of a single resolvent, we need to factor thecoefficient matrix only once, and matrix factors are hold and usedsequentially n times to solve a system of linear equations inside thefiltering, here n is degree of the polynomial of the filter.

By the use of a single resolvent for the filter instead of many
resolvents, we reduced both costs to factor matrix and especiallyto store matrix factors. However, properties of those filters which
are composed of a single resolvent are not good compared fromones composed of many resolvents, especially when the precisionof numbers and arithmetics used in computation is low.
The set of initial vectors generated from random numbers is
B-orthonormalized and then filtered to give another set of vectors,
to which we analyze and try to extract approximate pairs. If filter’s
properties are not good, approximate pairs obtained are inaccurate or some of them are lost, especially when the precision ofcomputation is low.
However, in the similar way as "orthogonal iteration" [3, 8, 9]
which is a well known method, we generate an initial set of vectors from random numbers, and to the set the combination of orthonormalization and filtering are applied a few times. The orthonormalization prevents the tendency to linear dependence ofthe set of vectors, and the filtering decrease content rates of thoseeigenvectors to be removed. By this refinement, the set of vectors
spans better approximation of the invariant-subspace spanned byrequired eigenvectors. From the set of refined vectors, the basisof approximate invariant-subspace is constructed, and to the basis
the Rayleigh-Ritz procedure is applied to obtain approximate pairs
required.
We made some experiments for a banded real symmetric-definite
generalized eigenproblem whose size of matrices is 210000 withlower-bandwidth 3051, which is a FEM discretization of the Laplacian eigenproblem in a cube with zero-Dirichlet boundary condition(Even banded matrices are sparse inside their bands, they are treatedas if their bands are dense). From experiments which used only
single precision for computations, even we used a filter whose
properties were not good since it was composed of only a single
resolvent to reduce requirements for computer resources, we foundthis approach of iterative refinement of eigenpairs worked well.

REFERENCES[1] Anthony P. Austin and Lloyd N. Trefethen. 2015. Computing eigenvalues of real
symmetric matrices with rational filters in real arithmetic. SIAM J. Sci. Comput37, 3 (2015), A1365–A1387.
[2] Martin Galgon, Lukas Krämer, and Brunno Lang. 2011. The FEAST algorithm forlarge eigenvalue problems. PAMM · Proc. Appl. Math. Mech. 11, 1 (2011), 747–748.
[3] Gene H. Golub and Charles F. van Loan. 2013. Matrix Computations (4 ed.). TheJohn Hopkins Univ. Press. §8.2.4:‘Orthogonal Iteration’.
[4] Stefan Güttel, Eric Polizzi, Ping Tak Peter Tang, and Gautier Viaud. 2015.
Zolotarev quadrature rules and load balancing for the FEAST eigensolver. SIAM
J. Sci. Comput 37, 4 (2015), A2100–A2122.
[5] Hiroshi Murakami. 2018. Filter diagonalization method by using a polynomialof a resolvent as the filter for a real symmetric-definite generalized eigenproblem. In Eigenvalue Problems: Algorithms, Software and Applications in PetascaleComputing, T. Sakurai, S.L. Zhang, T. Imamura, Y. Yamamoto, Y. Kuramashi, andT. Hoshi (Eds.). LNCSE, Vol. 117. Springer, 368–394.
[6] Hiroshi Murakami. 2019. Filters consists of a few resolvents to solve real
symmetric-definite generalized eigenproblems. JJIAM 36, 2 (2019). (EPASA2018Special Issue) to appear.
[7] Eric Polizzi. 2009. A density matrix-based algorithm for solving eigenvalue
problems. Phys. Rev. B 79, 1 (2009), 115112–115117.
[8] Heinz Rutishauser. 1969. Computational aspects of F. L. Bauer’s simultaneous
iteration method. Numer. Math. 13, 1 (1969), 4–13.
[9] Heinz Rutishauser. 1970. Simultaneous iteration method for symmetric matrices.
Numer. Math. 16, 3 (1970), 205–223.
[10] Sivan Toledo and Eran Rabani. 2002. Very large electronic structure calculations
using an out-of-core filter-diagonalization method. J. Comput. Phys. 180, 1 (2002),
256–269.