Stable Symmetric Eigenvalue Computation

Why Symmetry Matters in Eigenvalue Computation

Eigenvalues and eigenvectors of a symmetric matrix summarize important features of multivariate data, including directions of variation and the structure of estimated covariance matrices. In computation, however, an eigenvalue list is not enough to establish that the result is reliable. Rounding can affect both the matrix and the computed eigenvectors, and eigenvectors can be sensitive when eigenvalues are close. A useful assessment therefore combines the eigenpair residuals, the orthogonality of the computed vectors, and the size of the perturbation that would make the result exact.

For a real symmetric matrix \(A\in\mathbb{R}^{p\times p}\), the spectral theorem guarantees an orthogonal matrix \(Q\) and a diagonal matrix \(\Lambda\) such that \(A=Q\Lambda Q^T\). The columns of \(Q\) are unit eigenvectors, and the diagonal entries of \(\Lambda\) are real eigenvalues. A numerical eigensolver aims to produce approximations to this decomposition while respecting the symmetry that makes it well behaved.

How Symmetric Eigensolvers Use Structure

A standard dense symmetric eigensolver first reduces \(A\) to a symmetric tridiagonal matrix by orthogonal transformations. In exact arithmetic, this has the form \(T=U^TAU\), where \(U^TU=I\) and \(T\) is tridiagonal: its entries are zero more than one position from the diagonal. The eigenvalues of \(T\) are those of \(A\), and if \(Tz=\lambda z\), then \(Uz\) is an eigenvector of \(A\). Orthogonal transformations preserve vector lengths and do not amplify errors through repeated changes of scale in the way general, nonorthogonal transformations can.

The tridiagonal problem can then be solved with methods designed for its structure, such as implicit shifted QR iteration or divide-and-conquer methods. These methods differ in speed, memory use, and how they obtain eigenvectors. In practical use, a symmetric eigensolver should be selected for a symmetric input, rather than treating \(A\) as a general matrix. The precise numerical guarantees depend on the implementation and floating-point conditions, but well-designed symmetric algorithms aim for a backward-stable result: the computed decomposition should be exact for a matrix close to the input.

If the supplied matrix is only approximately symmetric because of data processing or rounding, inspect the discrepancy before proceeding. Replacing it by \((A+A^T)/2\) gives a symmetric matrix, but it is a change to the input, not a harmless correction in every application. Record the size of that change and make sure averaging is justified by the model.

What a Residual Can Guarantee

For a symmetric matrix, a small residual has a precise implication for eigenvalue proximity. It does not require that the candidate vector already be an exact eigenvector.

Proof. By the spectral theorem, write \(A=Q\Lambda Q^T\), with orthonormal eigenvectors \(q_j\) and eigenvalues \(\lambda_j\). Expand \(\hat q=\sum_{j=1}^p c_jq_j\), where \(\sum_j c_j^2=1\). Then \[ \|A\hat q-\hat\lambda\hat q\|_2^2 =\sum_{j=1}^p c_j^2(\lambda_j-\hat\lambda)^2. \] The right side is a weighted average of the squared distances from \(\hat\lambda\) to the eigenvalues. It cannot be smaller than the smallest such squared distance. Taking square roots proves the claim. \(\square\)

This result certifies closeness to some eigenvalue. It does not identify which eigenvalue, and it does not guarantee a small eigenvector error when eigenvalues are clustered. A gap condition is needed for that.

Proof. Expand \(\hat q=\sum_j c_jq_j\). By definition, \(\sin^2\theta=\sum_{j\ne k}c_j^2\). The residual expansion gives \[ \|A\hat q-\hat\lambda\hat q\|_2^2 =\sum_j c_j^2(\lambda_j-\hat\lambda)^2 \geq \delta^2\sum_{j\ne k}c_j^2 =\delta^2\sin^2\theta. \] Taking square roots and dividing by \(\delta>0\) proves the bound. \(\square\)

Orthogonality and Backward Error

For a full symmetric eigendecomposition, check both the eigenpair equations and the orthogonality of the reported vectors. Given a computed matrix \(\hat Q\) and diagonal matrix \(\hat\Lambda\), useful diagnostics include the residual norm \(\|A\hat Q-\hat Q\hat\Lambda\|_2\) or Frobenius norm, and the orthogonality error \(\|\hat Q^T\hat Q-I\|_2\) or Frobenius norm. Relative versions should scale these quantities by the size of the input where appropriate.

Proof. Since \(\hat Q\) is orthogonal, \(\hat Q^T\hat Q=I\), and \[ \hat Q^TR=\hat Q^TA\hat Q-\hat\Lambda. \] Both terms on the right are symmetric, so \(\hat Q^TR\) is symmetric. Also \(R\hat Q^T=\hat Q(\hat Q^TR)\hat Q^T\), which is symmetric; hence \(\Delta A\) is symmetric. From the definition of \(R\), \(A\hat Q=\hat Q\hat\Lambda+R\). Therefore \[ (A+\Delta A)\hat Q =A\hat Q-R\hat Q^T\hat Q =\hat Q\hat\Lambda+R-R =\hat Q\hat\Lambda. \] Finally, multiplication by an orthogonal matrix preserves the spectral norm, so \(\|\Delta A\|_2=\|R\hat Q^T\|_2=\|R\|_2\). \(\square\)

The theorem explains the backward-error interpretation of a small residual when the computed vectors are orthogonal. If \(\hat Q^T\hat Q\) is not close to \(I\), the theorem's premise fails: a small residual alone does not certify an accurate orthogonal eigendecomposition. In that case, report both diagnostics and investigate the computation, including whether the input was genuinely symmetric.

Small Gaps and a Common Mistake

A common mistake is to judge an eigensolver by whether its reported eigenvalues match another calculation to many decimal places, or to conclude that a small absolute residual guarantees accurate eigenvectors. Eigenvalue differences can be misleading when methods use different rounding, and eigenvectors associated with a tight cluster can change substantially while the residual remains small. Prefer residuals and orthogonality checks, scaled to the problem, and interpret eigenvector accuracy in light of the relevant spectral gap. Within a clustered eigenspace, the subspace may be more stable and meaningful than any particular choice of basis vectors.