Matrix Diagonalization

Matrix diagonalization turns a represented Hamiltonian into energies and stationary vectors. Its result is exact only for the supplied finite matrix, up to roundoff and solver error: agreement between two eigensolvers says nothing by itself about a spatial grid or basis cutoff. This page treats finite Hermitian matrices; Rayleigh–Ritz connects them to continuum operators.

Helpful background. The practical spectral theorem supplies orthogonal eigenspaces and spectral projectors.

Hermitian eigensystems and backward error

Let $H=H^\dagger$ be an $n\times n$ matrix in an orthonormal basis. We seek

$$HC=C\Lambda,\qquad C^\dagger C=I,$$

with real diagonal $\Lambda$. A structure-aware dense solver usually reduces $H$ by unitary transformations to real symmetric tridiagonal form, solves that smaller-structure problem, and transforms vectors back. It needs $O(n^2)$ storage and typically $O(n^3)$ work for the full eigensystem. A general nonsymmetric solver discards useful structure and can return small imaginary parts that distract from the real error analysis.

Request eigenvalues alone when vectors are unnecessary; use a subset-capable Hermitian driver for a selected interval. NumPy’s eigh returns ascending eigenvalues with their eigenvectors in columns. SciPy’s eigh also supports positive-definite generalized problems and selected subsets. Check the installed interface rather than assuming every library uses those conventions. Singular values and rank require an SVD, not a Hermitian eigensolve applied indiscriminately to the input matrix.

Before solving, measure $\|H-H^\dagger\|/\|H\|$. Do not silently replace a badly assembled matrix by $(H+H^\dagger)/2$: that can conceal an incorrect boundary term or conjugation. Some libraries read only one triangle, so the test must inspect the matrix actually constructed.

For a normalized candidate $c$ with real Rayleigh quotient $\theta=c^\dagger Hc$, form the residual $r=Hc-\theta c$. Expand $c=\sum_j a_j u_j$ in exact eigenvectors. Then

$$\|r\|^2=\sum_j |a_j|^2(E_j-\theta)^2, \qquad \min_j|E_j-\theta|\leq\|r\|.$$

This proves an absolute spectral-distance bound, without assuming an isolated eigenvalue. The scaled residual $\|r\|/(\|H\|+|\theta|)$ is a useful backward-error diagnostic, but it must be accompanied by an absolute tolerance when a physically important energy is close to zero. The LAPACK guide explains why backward stability and accuracy of individual vectors are different questions.

For a complete eigensystem, also check

$$R=HC-C\Lambda,\qquad O=C^\dagger C-I,\qquad D=H-C\Lambda C^\dagger.$$

Report $\|R\|_{\mathrm F}/\|H\|_{\mathrm F}$, $\|O\|_{\mathrm F}$, and $\|D\|_{\mathrm F}/\|H\|_{\mathrm F}$. If $H=0$, use absolute residual and reconstruction tolerances instead of dividing by zero. Reconstruction uses all eigenpairs; it is not a full-matrix accuracy test for a selected subset.

For normalized $c$, the Rayleigh-quotient residual has the physical meaning

$$\|(H-\theta I)c\|^2 =\langle H^2\rangle_c-\langle H\rangle_c^2 =(\Delta_c H)^2.$$

Thus residual norm is energy standard deviation. Near an eigenstate, evaluate the norm directly rather than subtracting two nearly equal expectation values. A small variance still does not select a unique eigenvector inside a tight spectral cluster.

Degeneracy and observable evaluation

If a target eigenvalue is separated from the rest of the spectrum by a gap $g$, eigenvector error is controlled in scale by $\|r\|/g$. Near degeneracy, tiny perturbations can rotate individual vectors substantially. Compare the projector $P=C_D C_D^\dagger$ onto the whole cluster, or its principal angles, rather than comparing vector entries.

For two orthonormal cluster bases $C_D$ and $\widetilde C_D$ of the same dimension, the singular values of $C_D^\dagger\widetilde C_D$ are the cosines of their principal angles. They compare the subspaces independently of the chosen basis within each cluster.

An expectation value is $c_j^\dagger Ac_j$. In a degenerate multiplet, an arbitrary solver basis need not diagonalize another commuting observable. Diagonalize $C_D^\dagger AC_D$ within that multiplet if definite quantum numbers are required. Sort eigenvalues and vector columns together; never sort one without the other.

Exploit known symmetry sectors before solving when possible. Check the symmetry commutator on the assembled matrix: a continuum symmetry can be broken by the chosen boundary conditions or truncation.

Along a parameter sweep, even an isolated eigenvector can acquire an arbitrary phase. If $c(s)^\dagger c(s+\delta s)\ne0$, align neighboring vectors by

$$c(s+\delta s)\longmapsto e^{-i\arg[c(s)^\dagger c(s+\delta s)]}c(s+\delta s).$$

This fixes a local plotting convention; it does not remove Berry holonomy. A small overlap warns that the parameter step or state assignment needs attention. Near a crossing, track overlaps, symmetry labels, or cluster projectors instead of energy order alone. Align a degenerate block as a whole, using the singular values and unitary factors of its overlap matrix, rather than aligning its columns independently.

Nonorthogonal representations

A nonorthogonal basis gives $Hc=ESc$, with $S>0$. Solve it as a Hermitian definite generalized problem or orthogonalize the basis. Rayleigh–Ritz derives the Cholesky reduction and explains how near-dependent basis vectors amplify assembly noise. Computing $S^{-1}H$ explicitly creates a generally non-Hermitian matrix in the Euclidean metric and worsens the numerical problem.

The diagnostics must respect the same overlap metric:

$$R=HC-SC\Lambda,\qquad C^\dagger SC\approx I.$$

For an individual vector, a scale-aware algebraic residual is

$$\eta_S=\frac{\|Hc-ESc\|} {(\|H\|+|E|\|S\|)\|c\|}.$$

It measures the supplied matrix problem, not the conditioning of the basis or continuum accuracy. For an absolute spectral-distance bound, the residual in the dual overlap metric uses $\sqrt{r^\dagger S^{-1}r/(c^\dagger Sc)}$, where $r=Hc-ESc$. In a truncated trial space this bounds distance to the projected spectrum; the full Hilbert-space residual additionally contains the component outside that space. A negative overlap eigenvalue inconsistent with rounding calls for checking the basis integrals; it is not repaired by choosing a different eigensolver.

Two-state verification examples

For

$$H=\begin{pmatrix}1&\epsilon\\\epsilon&1\end{pmatrix}, \qquad \epsilon>0,$$

the normalized eigenvectors are $(1,-1)^T/\sqrt2$ and $(1,1)^T/\sqrt2$, with energies $1-\epsilon$ and $1+\epsilon$. This checks ordering, normalization, and off-diagonal signs. At $\epsilon=0$ every orthonormal basis of the two-dimensional space is valid; a componentwise test against the previous vectors becomes invalid, while the projector remains $I$.

Check the largest residual, $\|C^\dagger C-I\|$, the Hermiticity defect, and the difference between $\operatorname{tr}H$ and $\sum_jE_j$. The trace check can reveal missing or duplicated eigenvalues, but matching sums alone cannot establish completeness: different erroneous lists can have the same sum. For large problems requiring only a few states, move to sparse eigensolvers.

Complex input tests a different possible failure: missing conjugation. For real $\Delta$ and complex $g$,

$$H=\begin{pmatrix}\Delta&g\\g^*&-\Delta\end{pmatrix}, \qquad E_\pm=\pm\sqrt{\Delta^2+|g|^2}.$$

With $\Delta=0.6$ and $g=0.8i$, the energies are $(-1,1)$. This short check tests the returned columns and full reconstruction:

import numpy as np

H = np.array([[0.6, 0.8j], [-0.8j, -0.6]], dtype=np.complex128)
energies, C = np.linalg.eigh(H)
Lambda = np.diag(energies)
scale = np.linalg.norm(H, ord="fro")
print(energies)
print(np.linalg.norm(H @ C - C @ Lambda, ord="fro") / scale)
print(np.linalg.norm(C.conj().T @ C - np.eye(2), ord="fro"))
print(np.linalg.norm(H - C @ Lambda @ C.conj().T, ord="fro") / scale)

The three errors should be near binary64 rounding scale for this well-conditioned small matrix. Eigenvector phases can differ across implementations. The exact trace is zero and determinant is $-1$, providing independent checks on the eigenvalue pair.

Separating solver error from the model error

Record basis ordering, units, dtype, symmetry sector, solver and backend, tolerances, and any state-tracking convention. Then distinguish four questions: whether the matrix was assembled correctly, whether its basis or grid is adequate, whether precision resolves the observable, and whether the algebraic solver has converged. A small residual answers only the last question.

For the finite-difference harmonic oscillator, vary grid spacing and box size independently. Refining the grid at fixed box does not remove boundary error; enlarging a box at fixed point count can make the grid coarser. Compare several low energies with $\hbar\omega(n+1/2)$ and require algebraic residuals below the discretization error. Reproducible projectors and observables are more informative than bitwise-identical eigenvector columns.

Exercises

For $H=\operatorname{diag}(0,2)$ and $c=(\cos\alpha,\sin\alpha)^T$, compute the energy error and residual near the ground state.

Solution

The Rayleigh quotient is $\theta=2\sin^2\alpha$, so the energy error is quadratic in $\alpha$. The residual norm is $2|\sin\alpha\cos\alpha|$, linear for small $\alpha$. A rapidly stabilizing energy can therefore conceal a less accurate vector and its observables.

Show that energy variance equals the squared Rayleigh-quotient residual, and explain why this does not prove that a trial state is the ground state.

Solution

For $c^\dagger c=1$ and $\theta=c^\dagger Hc$,

$$\begin{aligned} \|(H-\theta I)c\|^2 &=c^\dagger H^2c-2\theta c^\dagger Hc+\theta^2\\ &=\langle H^2\rangle_c-\theta^2. \end{aligned}$$

Every exact eigenstate has zero variance, including excited states. Identifying the ground state additionally requires spectral ordering, sector checks, or a suitable variational argument.

Design a two-parameter convergence test for oscillator eigenvalues obtained on a finite spatial grid. What would a tiny residual with an unchanged but inaccurate energy mean?

Solution

At each box half-width $L$, refine the spacing $h$ while tracking a fixed set of low states. Then increase $L$ at sufficiently small fixed $h$. Record errors against $\hbar\omega(n+1/2)$ and the original matrix residuals separately. A tiny residual with an inaccurate energy means the finite matrix has been solved accurately; a fixed box, unresolved spatial structure, or an assembly error can still limit the physical answer.

References

• Anderson, E., et al., LAPACK Users’ Guide, 3rd ed., SIAM, 1999, online guide — Hermitian eigensolvers, error bounds, and generalized problems.
• Bai, Z., J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, eds., Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, 2000, online text — solver selection and residual diagnostics.
• Golub, G. H., and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — dense reductions, generalized eigenproblems, and subspace comparisons.
• NumPy Developers, “numpy.linalg.eigh”, NumPy reference documentation — Hermitian input, triangle selection, ordering, and eigenvector-column conventions.
• SciPy Developers, “scipy.linalg.eigh”, SciPy reference documentation — definite generalized problems and subset-capable drivers.
• Stewart, G. W., and J.-G. Sun, Matrix Perturbation Theory, Academic Press, 1990 — spectral gaps, invariant-subspace sensitivity, and principal angles.