Matrix Diagonalization

Matrix diagonalization turns a finite matrix representation of a Hamiltonian or observable into numerical eigenvalues and eigenvectors. The library call is usually one line. Trustworthy use requires more: choose a solver that respects Hermiticity, understand what degeneracy makes nonunique, validate the returned subspaces, and separate eigensolver error from errors already present in the matrix model.

This page treats dense standard and generalized Hermitian eigenproblems. The algebraic conditions for diagonalizability live in Diagonalization. Large problems for which only part of the spectrum is needed belong to Sparse Eigensolvers.

Problem and Output

For a standard Hermitian eigenproblem,

$$Hv_n=E_nv_n, \qquad H^\dagger=H.$$

A complete solver returns an ordered vector of real eigenvalues and a unitary matrix of eigenvectors,

$$\Lambda = \operatorname{diag}(E_0,\ldots,E_{N-1}), \qquad V=(v_0\ \cdots\ v_{N-1}),$$

such that

$$HV=V\Lambda, \qquad V^\dagger V=I, \qquad H=V\Lambda V^\dagger.$$

Many numerical libraries store eigenvectors as columns of $V$. Verify that convention rather than guessing from the array shape.

For a Hermitian matrix, the mathematical eigenvalues are real. Tiny imaginary parts in quantities reconstructed by floating-point arithmetic should be judged against a scale-aware tolerance, not deleted without diagnosis.

Use a Hermitian Solver

If $H=H^\dagger$, use a Hermitian eigensolver. It exploits the matrix structure and is designed to return:

• real eigenvalues;
• orthonormal eigenvectors;
• better efficiency and storage use than a general nonsymmetric solver;
• error behavior appropriate to a normal matrix.

General eigensolvers solve a harder problem and may return small spurious imaginary parts or less orthogonal vectors. Calling one on a Hermitian Hamiltonian discards valuable information.

Do not silently replace an imperfect input by

$$H_{\mathrm{sym}} = \frac12(H+H^\dagger)$$

unless that operation is part of the mathematical model. First measure

$$\eta_{\mathrm H} = \frac{\lVert H-H^\dagger\rVert} {\max(\lVert H\rVert,\epsilon_{\mathrm{scale}})}.$$

A large $\eta_{\mathrm H}$ may expose an assembly error, a missing complex conjugate, an inconsistent boundary term, or a genuinely non-Hermitian model. Symmetrization can hide each of these.

What Dense Solvers Do

Production routines do not usually find roots of the characteristic polynomial. For a dense Hermitian matrix they typically:

1. reduce $H$ by unitary transformations to a real symmetric or complex Hermitian tridiagonal problem;
2. solve the tridiagonal eigenproblem using a QR-family, divide-and-conquer, bisection/inverse-iteration, or relatively robust representation method;
3. back-transform the eigenvectors when they are requested.

The exact path depends on the library and selected driver. For an $N\times N$ dense matrix, storing the matrix requires $O(N^2)$ memory and a complete eigendecomposition requires $O(N^3)$ arithmetic. Computing eigenvalues without eigenvectors can reduce constants and memory, but it does not change dense asymptotic scaling.

This scaling is the practical boundary of dense diagonalization. A tensor-product Hilbert space can make $N$ grow exponentially even when each local subsystem is small.

Solver Selection

Mathematical task	Appropriate numerical family
All eigenpairs of dense Hermitian $H$	Dense Hermitian eigensolver
Eigenvalues only	Hermitian eigenvalue-only routine
Selected interval of a dense Hermitian spectrum	Hermitian subset-capable driver
$Hc=ESc$ with positive-definite $S$	Generalized Hermitian eigensolver
A few extremal eigenpairs of large sparse $H$	Lanczos or another sparse Hermitian method
Singular values, rank, or least squares	Singular-value decomposition
Genuinely non-Hermitian operator	General or structure-specific non-Hermitian solver

NumPy’s current eigh interface solves real symmetric or complex Hermitian problems, returns eigenvalues in ascending order, and places the associated normalized eigenvectors in columns. SciPy’s eigh additionally supports generalized Hermitian problems and selected eigenvalue subsets. Treat these as interface facts to verify in the installed library documentation, not as universal conventions shared by every language.

Reliable Workflow

A reproducible dense calculation has four stages.

Assemble. Fix basis ordering, units, dtype, boundary conditions, and every normalization factor before diagonalization.

Inspect. Check shape, finiteness, Hermiticity, expected sparsity or block structure, and a few known matrix elements.

Solve. Use the structure-aware routine. Request only eigenvalues when vectors are unnecessary.

Validate and interpret. Scale residuals, test orthogonality and reconstruction, handle degenerate subspaces, and compare against analytic or convergence benchmarks.

The eigensolver only answers the finite matrix problem it receives. A tiny residual cannot prove that a basis truncation or spatial discretization represents the intended continuum operator.

Worked Two-Level Benchmark

Consider

$$H = \begin{pmatrix} \Delta & g\\ g^* & -\Delta \end{pmatrix}.$$

The characteristic equation gives

$$E_\pm = \mathord\pm\sqrt{\Delta^2+\lvert g\rvert^2}.$$

This model is a useful unit test because it exercises complex Hermitian input while retaining an exact answer. Choose $\Delta=0.6$ and $g=0.8i$. Then $E_\pm=\mathord\pm1$.

import numpy as np

H = np.array(
    [[0.6, 0.8j],
     [-0.8j, -0.6]],
    dtype=np.complex128,
)

energies, vectors = np.linalg.eigh(H)
Lambda = np.diag(energies)

residual = H @ vectors - vectors @ Lambda
orthogonality = vectors.conj().T @ vectors - np.eye(2)
reconstruction = H - vectors @ Lambda @ vectors.conj().T

print(energies)
print(np.linalg.norm(residual, ord="fro"))
print(np.linalg.norm(orthogonality, ord="fro"))
print(np.linalg.norm(reconstruction, ord="fro"))

The expected eigenvalue array is $(-1,1)$ up to rounding. The three norms should be small compared with the scale of $H$. The individual eigenvector phases may differ across libraries or runs without changing any physical prediction.

Residuals and Backward Error

For an approximate eigenpair $(\widetilde E,\widetilde v)$, the residual is

$$r = H\widetilde v-\widetilde E\widetilde v.$$

An absolute residual has units and changes if $H$ is rescaled. A dimensionless normwise backward-error indicator is

$$\eta = \frac{\lVert r\rVert_2} {\lVert H\rVert_2\lVert\widetilde v\rVert_2 +\lvert\widetilde E\rvert \lVert\widetilde v\rVert_2}.$$

Small $\eta$ means the returned pair is an exact eigenpair of a nearby matrix problem. It does not, by itself, guarantee a small forward error in the eigenvector. Forward sensitivity also depends on spectral gaps.

For a complete eigensystem, useful matrix-level diagnostics are

$$\begin{aligned} R&=HV-V\Lambda,\\ O&=V^\dagger V-I,\\ D&=H-V\Lambda V^\dagger. \end{aligned}$$

Report scaled norms such as

$$\frac{\lVert R\rVert_{\mathrm F}} {\lVert H\rVert_{\mathrm F}}, \qquad \lVert O\rVert_{\mathrm F}, \qquad \frac{\lVert D\rVert_{\mathrm F}} {\lVert H\rVert_{\mathrm F}},$$

with a stated treatment for the exceptional case $\lVert H\rVert_{\mathrm F}=0$.

Rayleigh Quotient and Energy Variance

For any nonzero trial vector $v$, the Rayleigh quotient is

$$\rho(v) = \frac{v^\dagger Hv}{v^\dagger v}.$$

For normalized $v$ and Hermitian $H$, the residual formed with $\rho(v)$ obeys

$$\begin{aligned} \lVert(H-\rho) v\rVert_2^2 &=\langle H^2\rangle_v-\langle H\rangle_v^2\\ &=(\Delta_v H)^2. \end{aligned}$$

Thus the energy variance is exactly the squared residual norm for a normalized trial state evaluated at its Rayleigh quotient. This links a standard numerical diagnostic to a physical statement: an exact energy eigenstate has zero energy variance.

A small variance confirms proximity to some spectral subspace. If several eigenvalues are clustered, it need not identify a unique eigenvector inside that subspace.

Degeneracy and Nonuniqueness

If an eigenvalue $E$ has multiplicity $d$, every orthonormal basis of its eigenspace is valid. If $V_D$ contains one numerical basis for that subspace and $W\in U(d)$, then

$$V_D \longmapsto V_DW$$

changes the returned vectors but not the eigenspace. Comparing individual vectors in an exactly or nearly degenerate cluster is therefore unreliable.

Compare the spectral projector instead:

$$P_D=V_DV_D^\dagger.$$

Two computations agree on the subspace when their projectors agree. Principal angles provide a more detailed comparison: the singular values of $V_D^\dagger \widetilde V_D$ are the cosines of the principal angles between the two subspaces.

If an additional Hermitian symmetry $Q$ commutes with $H$, diagonalizing the restriction of $Q$ inside the degenerate eigenspace can select physically meaningful labels. This is a basis choice supplied by extra structure, not by $H$ alone.

Near Degeneracy and Spectral Gaps

Let an isolated eigenvalue or cluster be separated from the rest of the spectrum by a gap $g$. A perturbation $\delta H$ can rotate its invariant subspace by an amount controlled schematically by

$$\sin\Theta \lesssim \frac{\lVert\delta H\rVert}{g}.$$

When $g$ is small, eigenvectors can change substantially even though eigenvalues and the combined clustered subspace remain accurate. This is why a tiny residual and a visually different eigenvector are not contradictory near degeneracy.

Do not decide degeneracy from a fixed number of decimal places. Compare the observed splitting with relevant scales:

• matrix norm and machine precision;
• estimated assembly or discretization error;
• symmetry expectations;
• convergence under basis or grid refinement.

Phases, Signs, and State Tracking

An isolated normalized eigenvector is defined only up to phase:

$$v_n\longmapsto e^{i\phi_n}v_n.$$

Real symmetric solvers show the same freedom as an arbitrary sign. Component-by-component comparison can therefore fail even for identical physical states.

For a parameter-dependent Hamiltonian $H(s)$, a simple phase alignment for an isolated state is

$$v_n(s+\delta s) \longmapsto e^{-i\arg\langle v_n(s),v_n(s+\delta s)\rangle} v_n(s+\delta s).$$

Near crossings, sorting by eigenvalue alone can swap state labels. Match states by overlaps, symmetry quantum numbers, or continuity of spectral projectors. For a degenerate block, align whole subspaces with a unitary Procrustes or singular-value-decomposition step rather than aligning columns independently.

This local gauge fixing aids plotting and differentiation. It does not remove global geometric effects such as Berry phase.

Generalized Hermitian Eigenproblems

A nonorthogonal basis produces

$$Hc=ESc,$$

where the overlap matrix satisfies

$$S=S^\dagger>0.$$

The eigenvectors are normalized in the $S$ metric:

$$c_m^\dagger S c_n=\delta_{mn}.$$

If $S=LL^\dagger$ is a Cholesky factorization, define $y=L^\dagger c$. Then

$$\left(L^{-1}HL^{-\dagger}\right)y=Ey, \qquad c=L^{-\dagger}y.$$

In software, use a generalized Hermitian driver or triangular solves; do not explicitly form matrix inverses. If $S$ is not positive definite or is extremely ill-conditioned, the basis may contain exact or near linear dependencies. That is a modeling and conditioning problem, not merely an eigensolver inconvenience.

Validation must use the generalized residual and metric:

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

Symmetry Blocks Before Diagonalization

If a Hermitian operator $Q$ commutes with $H$,

$$[H,Q]=0,$$

the Hilbert space can often be decomposed into invariant sectors before solving. Block diagonalization:

• reduces time and memory;
• prevents mixing between distinct symmetry labels;
• makes degeneracies easier to interpret;
• supplies more stable state labels across parameter sweeps.

Numerically check the commutator relative to matrix scales. A purported symmetry that fails this test may have been broken by boundary conditions, truncation, or an assembly error.

Sorting and Labeling Eigenpairs

Never sort eigenvalues without applying the same permutation to eigenvector columns. If $p$ is the sorting permutation,

$$E_n\longmapsto E_{p(n)}, \qquad v_n\longmapsto v_{p(n)}.$$

Ascending energy is often useful, but it is not always the best label. Near exact crossings, states can be tracked by symmetry sectors. Near avoided crossings, overlap continuation may better preserve identity. When a reported state is defined by an observable $A$, record

$$\langle A\rangle_n=v_n^\dagger A v_n$$

alongside its energy rather than relying on array position alone.

Error Budget

Four distinct error sources should be separated:

Source	Typical diagnostic
Matrix assembly	Hermiticity, units, known entries, symmetry commutators
Basis truncation or discretization	Convergence under basis size, grid spacing, or domain size
Floating-point conditioning	Scale changes, precision changes, gap estimates
Eigensolver	Residual, orthogonality, reconstruction

A backward-stable eigensolver can solve the wrong discretized Hamiltonian extremely accurately. Conversely, physically converged low-energy observables may remain useful even when high-energy grid states are poor.

Harmonic-Oscillator Benchmark

On an interior uniform grid with spacing $\Delta x$ and Dirichlet endpoints, the centered finite-difference oscillator Hamiltonian has

$$\begin{aligned} H_{jj} &=\frac{\hbar^2}{m(\Delta x)^2} +\frac12m\omega^2x_j^2,\\ H_{j,j\mathord\pm1} &=-\frac{\hbar^2}{2m(\Delta x)^2}. \end{aligned}$$

Its low eigenvalues should approach

$$E_n=\hbar\omega\left(n+\frac12\right).$$

A credible convergence study varies both the grid spacing and the finite domain size. Decreasing $\Delta x$ at fixed box size reduces discretization error but not boundary truncation error. Increasing the box at fixed number of points can make the grid coarser. Coordinate the two limits and monitor several low-lying levels.

The construction belongs to Finite Difference Methods. Imaginary-Time Projection Notebook uses the resulting eigensystem as an independent benchmark for an iterative ground-state calculation.

Reproducibility Checklist

Record enough information to rerun and audit the calculation:

• basis and basis ordering;
• physical units and nondimensionalization;
• matrix dimension and dtype;
• symmetry sectors and truncation rules;
• library and backend versions;
• routine and nondefault solver options;
• whether eigenvectors were requested;
• residual, orthogonality, and reconstruction tolerances;
• convergence data for the physical discretization;
• the rule used to sort, phase-align, or track states.

Bitwise-identical eigenvectors are not a reasonable portability requirement, especially in degenerate spaces. Reproducible eigenvalues, residuals, projectors, and observables are stronger scientific targets.

Common Mistakes

• Using a general eigensolver for a Hermitian matrix.
• Trusting the input triangle without checking how the library interprets it.
• Silently symmetrizing a matrix before diagnosing the anti-Hermitian part.
• Comparing degenerate eigenvectors column by column instead of comparing subspaces.
• Forgetting that eigenvectors are arbitrary up to phase or sign.
• Sorting eigenvalues without permuting the corresponding eigenvectors.
• Treating a small residual as proof of continuum or basis convergence.
• Using a fixed absolute tolerance across differently scaled Hamiltonians.
• Forming $A^\dagger A$ to replace an SVD and thereby squaring the condition number.
• Forming explicit inverses in a generalized eigenproblem.
• Applying dense diagonalization when only a few eigenpairs of a large sparse matrix are needed.

Boundary of This Page

This page owns the workflow for dense Hermitian matrix eigenproblems. It does not duplicate the mathematical spectral theorem, the construction of a particular Hamiltonian, or large-scale iterative methods. See Spectral Decomposition for the algebra, Discretization for continuum-to-matrix modeling, and Sparse Eigensolvers for selected eigenpairs at scale.

Cross-Links

• Eigenvalues and Eigenvectors
• Hermitian Operators
• Diagonalization
• Spectral Decomposition
• Norms and Metrics
• Floating-Point Arithmetic
• Conditioning and Stability
• Discretization
• Finite Difference Methods
• Sparse Matrices
• Sparse Eigensolvers
• Singular Value Decomposition
• Matrix Exponentials Numerically
• Exact Diagonalization Preview

References

• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, 1997.
• G. W. Stewart and J.-G. Sun, Matrix Perturbation Theory, Academic Press, 1990.
• E. Anderson et al., LAPACK Users’ Guide, 3rd ed., SIAM, 1999.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• NumPy documentation: numpy.linalg.eigh.
• SciPy documentation: scipy.linalg.eigh.

Exercises

1. Derive the eigenvalues of the two-level benchmark and check them using trace and determinant.

Solution

The characteristic determinant is

$$\begin{aligned} \det(H-EI) &=(\Delta-E)(-\Delta-E)-\lvert g\rvert^2\\ &=E^2-\Delta^2-\lvert g\rvert^2. \end{aligned}$$

Therefore

$$E_\pm=\mathord\pm\sqrt{\Delta^2+\lvert g\rvert^2}.$$

Their sum is zero, matching $\operatorname{Tr}H=0$, and their product is

$$E_+E_- =-\Delta^2-\lvert g\rvert^2 =\det H.$$

2. Prove that energy variance equals squared residual norm at the Rayleigh quotient.

Solution

Let $v$ be normalized and set $\rho=v^\dagger Hv$. Since $H$ is Hermitian, $\rho$ is real. Then

$$\begin{aligned} \lVert(H-\rho)v\rVert_2^2 &=v^\dagger(H-\rho)^2v\\ &=v^\dagger H^2v -2\rho\,v^\dagger Hv+\rho^2\\ &=\langle H^2\rangle_v-\rho^2\\ &=(\Delta_vH)^2. \end{aligned}$$

3. Show why individual vectors are not invariant inside a degenerate eigenspace.

Solution

Let $HV_D=EV_D$ for a $d$-column orthonormal basis $V_D$, and let $W\in U(d)$. Then

$$H(V_DW) =(HV_D)W =E(V_DW).$$

Also,

$$(V_DW)^\dagger(V_DW) =W^\dagger W=I_d.$$

Thus $V_DW$ is another valid orthonormal eigenbasis. Its projector is unchanged:

$$(V_DW)(V_DW)^\dagger =V_DWW^\dagger V_D^\dagger =V_DV_D^\dagger.$$

4. Reduce a generalized Hermitian eigenproblem to a standard one.

Solution

For $S>0$, take $S=LL^\dagger$ and define $y=L^\dagger c$, so $c=L^{-\dagger}y$. Substituting into $Hc=ESc$ gives

$$HL^{-\dagger}y =ELL^\dagger L^{-\dagger}y =ELy.$$

Multiplying by $L^{-1}$ yields

$$\left(L^{-1}HL^{-\dagger}\right)y=Ey.$$

The transformed matrix is Hermitian. If the $y$ vectors are orthonormal, then

$$c_m^\dagger Sc_n =y_m^\dagger y_n =\delta_{mn}.$$

5. Design a convergence test for the finite-difference harmonic oscillator.

Solution

Choose several box half-widths $L$ and several grid spacings $\Delta x$. For each pair, build the Hermitian tridiagonal matrix, compute a fixed number of low eigenvalues, and record:

$$\frac{\lvert E_n^{\mathrm{num}}-E_n^{\mathrm{exact}}\rvert} {E_n^{\mathrm{exact}}}, \qquad \frac{\lVert Hv_n-E_nv_n\rVert_2} {\lVert H\rVert_2}.$$

At fixed sufficiently large $L$, refine $\Delta x$ to expose discretization convergence. At fixed sufficiently small $\Delta x$, enlarge $L$ to expose boundary truncation. A level is credible only when both studies stabilize it and the eigensolver residual remains much smaller than the physical discretization error.