Unresolved reliability and generality issues of the Yau-Lu unitary-embedding method for computing all eigenvalues of large sparse symmetric matrices

In 1993 Yau and Lu proposed computing all eigenvalues of a large sparse symmetric matrix by embedding it into a unitary matrix $B=[A+i\sqrt{I-A^2}]^2$ and applying Chebyshev-accelerated Fourier extraction to the trajectory $B^k v_0$, aiming for $O(n)$ matrix-vector products and $O(n)$ memory when only eigenvalues are required. The authors themselves listed the method's unresolved issues: the eigenvalue selection scheme can return fewer than $n$ eigenvalues when eigenvalues are multiple or clustered; accuracy depends on the random starting vector with no systematic guarantee; the parameters $N$ and $\delta$ must be chosen heuristically; the method needs reliable extreme-eigenvalue estimates to scale the spectrum into $[-0.9,0.9]$ and the polynomial approximation of $\sqrt{1-x^2}$ degrades outside this interval; and extensions to generalized and non-symmetric problems were not addressed. These listed limitations have not been resolved. The problem of making this spectral-transform approach a dependable, automatic, and general solver - or proving which of its goals are unattainable - remains open.

Background

Computing the eigenvalues of a very large sparse symmetric matrix is a core task of scientific computing: discretizations of vibration, quantum, and network problems produce matrices of dimension $n$ in the millions whose entries are mostly zero. Dense solvers, which first reduce the matrix to tridiagonal form, cost $O(n^3)$ operations and $O(n^2)$ storage and are impractical, so iterative methods are used that touch the matrix only through sparse matrix-vector multiplications. The Lanczos method builds an orthonormal basis of the Krylov subspace spanned by $\{v_0, Av_0, A^2v_0, \dots\}$ and projects $A$ onto a small tridiagonal matrix; it computes a few extreme eigenvalues very quickly, but interior eigenvalues converge slowly, and in floating-point arithmetic the computed basis vectors lose orthogonality, producing spurious duplicate ('ghost') eigenvalues. Computing all $n$ eigenvalues reliably with Lanczos therefore forces repeated reorthogonalization of the basis, at $O(n^2)$ memory cost.

Yau and Lu (1993) proposed a different route. The matrix is first scaled and shifted so that its spectrum lies in $[-0.9,0.9]$; each eigenvalue is then written as $\lambda_j=\cos\mu_j$. Since the function $\sqrt{1-x^2}$ satisfies $\cos\mu+i\sqrt{1-\cos^2\mu}=e^{i\mu}$, the matrix $B=[A+i\sqrt{I-A^2}]^2$ is a unitary matrix of the form $e^{iF}$ with $A=\cos(F/2)$, and its eigenvalues $e^{2i\mu_j}$ lie on the unit circle: the real spectrum has been converted into a set of rotation angles. On the unit circle there is no growth or decay, so a polynomial in the angle variable can isolate each eigencomponent. The identity $B+B^{-1}=4A^2-2I$ yields the three-term recurrence $v_{k+1}=(4A^2-2I)v_k-v_{k-1}$, which generates the trajectory $v_k=B^k v_0$ using only sparse matrix-vector products; only the first step $v_1=Bv_0$ needs one application of $\sqrt{I-A^2}$, performed once with a tabulated degree-12 minimax polynomial for $\sqrt{1-x^2}$ (coefficients from the Remez exchange algorithm, the classical iterative method for best uniform polynomial approximation), accurate to about $10^{-8}$ on $|x|\le 0.9$ - the endpoints $\pm 1$ are excluded because $\sqrt{1-x^2}$ has infinite derivative there. Eigencomponents are then extracted by a Chebyshev filter: the Chebyshev polynomial $T_m$, defined by $T_m(\cos\theta)=\cos(m\theta)$, yields after rescaling the polynomial $Q_N$ of degree $N-1$ that equals 1 at a target angle and is uniformly smallest on all angles farther than $\delta\pi$, where $\delta$ (typically $O(1/n)$) is tied to the smallest expected eigenvalue spacing. The filtered vectors $u(\mu)=\operatorname{Re}\sum_{k<N}\beta_k e^{-ik\mu}B^k v_0$, where the $\beta_k$ are the real filter weights determined by $Q_N$, are evaluated simultaneously at the $2N$ uniform angles $\mu_l=\pi l/N$, $l=0,\dots,2N-1$, by a fast Fourier transform (specifically a discrete cosine transform, the real-input variant of the Fourier transform, applied to the sequence $\beta_k v_k$ and costing $O(nN\log N)$ instead of $O(nN^2)$). Candidate eigenvalues are obtained from Rayleigh quotients $u^TAu/u^Tu$ of the filtered vectors, and a selection step keeps the local minima of $|\tilde{\mu}_l-\pi l/N|$, where $\tilde{\mu}_l$ is the angle implied by the Rayleigh quotient. When only eigenvalues are needed, memory drops to $O(n)$ by storing just the components of the $v_k$ on a sparsest row of $A$ and forming the Rayleigh-type ratios along that row.

This cost profile was the point of the method: all $n$ eigenvalues with $O(n)$ sparse matrix-vector products and $O(n)$ storage, versus $O(n^2)$ storage for reorthogonalized Lanczos. In the authors' numerical tests (a tridiagonal matrix with the exactly known spectrum $\lambda_j=0.9\cos(j\pi/(n+1))$, $n=128$, and a randomly perturbed diagonal matrix, $n=256$, both with $N=8n$) most eigenvalues were computed to about $10^{-12}$, but the extreme eigenvalues of the tridiagonal, whose spacing is $O(1/n^2)$, were not resolved to high accuracy. The authors listed the method's remaining limitations: (1) the selection scheme can return fewer than $n$ eigenvalues when eigenvalues are multiple or closer than the angular grid $\pi/N$; (2) the accuracy depends on the random initial vector $v_0$ and no systematic way to guarantee accuracy is provided; (3) the parameters $N$ and $\delta$ must be chosen from the eigenvalue distribution (for example the smallest eigenvalue spacing), but no automatic or adaptive procedure is given; (4) reliable estimates of the extreme eigenvalues are needed to scale $A$ into $[-0.9,0.9]$, and if the true spectrum exceeds this interval the polynomial approximation of $\sqrt{1-x^2}$ degrades, the restriction to $|x|\le 0.9$ being tied to the infinite derivative of $\sqrt{1-x^2}$ at $x=\pm 1$; (5) the method is designed for the standard symmetric problem only, and extensions to generalized problems $Ax=\lambda Mx$ or to non-symmetric matrices are not addressed; and (6) as positioning, the method is aimed at computing all eigenvalues - for a few extreme ones, Lanczos is more efficient. The companion journal paper of the same year develops the dense-matrix, matrix-multiplication-based variant of the same idea.

Problem Statement

The problem is the one the source authors left open: turn the Yau-Lu spectral-transform approach into a dependable method for computing all eigenvalues (and, where claimed, eigenvectors) of a large sparse symmetric matrix within its intended cost profile - $O(n)$ sparse matrix-vector products and $O(n)$ memory in the eigenvalues-only mode (the mode that stores only the trajectory components needed to form eigenvalue estimates, not the full vector sequence), $n$ being the matrix dimension - by resolving the limitations the authors themselves listed. These components serve one common objective, a trustworthy all-eigenvalue solver of this type; below, 'this class of schemes' denotes spectral-transform eigensolvers that form eigenvalue estimates from a Chebyshev-filtered trajectory $B^kv_0$ evaluated on the uniform angular grid $\mu_l=\pi l/N$, the Yau-Lu construction being the prototype. Components 1-4 are independently standing open issues of the same program, each admitting either a resolution of the goal it names or a proof that the goal is unattainable for this class of schemes; component 5 is secondary.

1. Completeness of eigenvalue selection. Replace or repair the local-minimum selection of $|\tilde{\mu}_l-\pi l/N|$ (where $\mu_l=\pi l/N$, $l=0,\dots,2N-1$, is the uniform angular grid and $\tilde{\mu}_l$ is the angle implied by the Rayleigh quotient of the filtered vector at $\mu_l$) so that the method provably returns all $n$ eigenvalues with their multiplicities, including eigenvalues whose angular separation is below the grid spacing $\pi/N$ (multiple or clustered eigenvalues), or prove that this completeness goal - returning all $n$ eigenvalues with their multiplicities, including eigenvalues separated below the grid spacing $\pi/N$ - is impossible for this class of schemes.
2. Guaranteed accuracy. Provide a systematic accuracy guarantee - an a priori bound or a certified a posteriori error bound independent of the random starting vector $v_0$ - including a rounding-error analysis of the unreorthogonalized three-term recurrence $v_{k+1}=(4A^2-2I)v_k-v_{k-1}$ (which generates the trajectory $v_k=B^kv_0$ from the starting vector $v_0$, $A$ being the scaled input matrix) over $N=O(n)$ steps and of the effect of the polynomial-approximation error in $\sqrt{1-x^2}$.
3. Automatic parameters. An adaptive or automatic rule for choosing $N$ (the number of filter terms, which fixes the angular grid $\pi l/N$) and $\delta$ (the angular gap parameter of the filter) that does not require prior knowledge of the eigenvalue distribution, with a termination or failure criterion.
4. Robust spectral scaling. Remove or soften the dependence on accurate extreme-eigenvalue estimates and the $|x|\le 0.9$ restriction without abandoning the cost profile.
5. Generality (a secondary component of the same program). Extend the construction to generalized symmetric problems $Ax=\lambda Mx$ (the matrix-pencil problem with a second symmetric matrix $M$) and to non-symmetric matrices, or identify the obstruction.

The source question concerns multiple and clustered eigenvalues and the absence of accuracy guarantees. The authors' positioning remark - that the method targets all eigenvalues while Lanczos is more efficient for a few extreme eigenvalues - is not part of the question.

References

1. Shing-Tung Yau and Ya Yan Lu, "A New Approach to Sparse Matrix Eigenvalues", Proceedings of the First IEEE Regional Conference on Aerospace Control Systems, 1993, pp. 132-137. DOI: 10.1109/aerocs.1993.720909, https://doi.org/10.1109/aerocs.1993.720909 (source paper)
2. Shing-Tung Yau and Ya Yan Lu, "Reducing the Symmetric Matrix Eigenvalue Problem to Matrix Multiplications", SIAM Journal on Scientific Computing 14(1):121-136, 1993. DOI: 10.1137/0914008, https://doi.org/10.1137/0914008
3. Françoise Tisseur, "Parallel Implementation of the Yau and Lu Method for Eigenvalue Computation", The International Journal of Supercomputer Applications and High Performance Computing 11(3):197-204, 1997. DOI: 10.1177/109434209701100302, https://doi.org/10.1177/109434209701100302
4. Stéphane Domas and Françoise Tisseur, "Parallel implementation of a symmetric eigensolver based on the Yau and Lu method", in Vector and Parallel Processing - VECPAR'96, Lecture Notes in Computer Science 1215, Springer, pp. 140-153, 1997. DOI: 10.1007/3-540-62828-2_117, https://doi.org/10.1007/3-540-62828-2_117
5. Françoise Tisseur, "Méthodes numériques pour le calcul d'éléments spectraux : étude de la précision, la stabilité et la parallélisation", doctoral thesis, Saint-Étienne, France, 1997 (no DOI; theses.fr record 1997STET4006). https://theses.fr/1997STET4006
6. Ya Yan Lu, "A Padé Approximation Method for Square Roots of Symmetric Positive Definite Matrices", SIAM Journal on Matrix Analysis and Applications 19(3):833-845, 1998. DOI: 10.1137/S089547989731631X, https://doi.org/10.1137/S089547989731631X
7. Ya Yan Lu, "A complex coefficient rational approximation of 1+x", Applied Numerical Mathematics 27(2):141-154, 1998. DOI: 10.1016/S0168-9274(98)00009-9, https://doi.org/10.1016/S0168-9274(98)00009-9
8. Ruipeng Li, Yuanzhe Xi, Eugene Vecharynski, Chao Yang, and Yousef Saad, "A Thick-Restart Lanczos Algorithm with Polynomial Filtering for Hermitian Eigenvalue Problems", SIAM Journal on Scientific Computing 38(4):A2512-A2534, 2016. DOI: 10.1137/15M1054493, https://doi.org/10.1137/15M1054493
9. Ruipeng Li, Yuanzhe Xi, Lucas Erlandson, and Yousef Saad, "The Eigenvalues Slicing Library (EVSL): Algorithms, Implementation, and Software", SIAM Journal on Scientific Computing 41(4):C393-C415, 2019 (preprint arXiv:1802.05215, 2018). DOI: 10.1137/18M1170935, https://doi.org/10.1137/18M1170935; preprint: https://arxiv.org/abs/1802.05215
10. Andreas Pieper, Moritz Kreutzer, Andreas Alvermann, Martin Galgon, Holger Fehske, Georg Hager, Bruno Lang, and Gerhard Wellein, "High-performance implementation of Chebyshev filter diagonalization for interior eigenvalue computations", Journal of Computational Physics 325:226-243, 2016. DOI: 10.1016/j.jcp.2016.08.027, https://doi.org/10.1016/j.jcp.2016.08.027