Eigenvalues and Eigenvectors

An eigenvector of a square matrix is a non-zero vector that the matrix only stretches, shrinks or flips, without turning it; the factor is its eigenvalue. Writing a matrix in terms of its eigenvectors and eigenvalues, when that is possible, is its eigendecomposition (also called the eigenvalue decomposition, EVD, or spectral decomposition). Many of the matrices in computer vision are symmetric: covariance matrices of pixels, patches or points, the structure tensor of image gradients, Hessians of images and of cost functions, and graph Laplacians. For them the eigenvectors are orthogonal axes and the eigenvalues say how much happens along each, which is why principal component analysis, corner detection, normal estimation and spectral segmentation are all eigenproblems. Linear algebra gives the definition and the spectral theorem in outline; this article covers the general case, symmetric matrices in more depth, how much noise moves eigenvalues and eigenvectors, how the decomposition is computed and what the libraries return.

Definition

For a square matrix $A \in \mathbb{R}^{n \times n}$, a scalar $\lambda$ and a non-zero vector $\mathbf{u}$ with $A\mathbf{u} = \lambda \mathbf{u}$ are an eigenvalue and an eigenvector of $A$ (Murphy, 2022, §7.4.1). Any non-zero multiple of an eigenvector is again one, so eigenvectors are usually normalized to unit length, which still leaves the sign free; Murphy notes that “we will have to live with this” ambiguity between $\mathbf{u}$ and $-\mathbf{u}$. When several independent eigenvectors share one eigenvalue, only the subspace they span is determined.

Intuition

A linear map turns most vectors. The eigenvectors are the exceptions: the directions it leaves on their own line. For a symmetric matrix there are $n$ such directions and they are perpendicular, so the map is a rotation into those axes, a scaling along each, and a rotation back (Murphy, 2022, §7.4.3). The unit circle becomes an ellipse whose axes are the eigenvectors and whose semi-axes are the absolute eigenvalues.

The structure tensor of an image window is the standard picture in vision. Its dominant eigenvector points along the direction in which intensity changes most and its eigenvalues measure how much the intensity changes along each axis: both small in a flat region, one large on an edge, both large at a corner, as the Harris corner detector uses. When the two eigenvalues are equal, as on a pattern with equal contrast in two perpendicular directions, there is no dominant direction at all, and any estimate of one is decided by the noise. How close two eigenvalues are therefore decides how well their eigenvectors are determined, which is the subject of the example below.

A general matrix need not behave this way. A rotation of the plane turns every real vector, so it has no real eigenvectors; a shear tilts every vector except those along one line, so it has too few.

Formal Definition

The characteristic polynomial and multiplicities

$A\mathbf{u} = \lambda\mathbf{u}$ with $\mathbf{u} \ne 0$ means that $\lambda I - A$ has a non-trivial null space, so it is singular:

$$\det(\lambda I - A) = 0 .$$

This is the characteristic equation; its left side is a polynomial of degree $n$ in $\lambda$, and its $n$ roots, which may be complex, are the eigenvalues (Murphy, 2022, §7.4.1). The sum of the eigenvalues is the trace and their product is the determinant; the eigenvalues of a triangular matrix are its diagonal entries (§7.4.1).

Two counts describe an eigenvalue, in the standard terms. Its algebraic multiplicity is its multiplicity as a root of the characteristic polynomial; its geometric multiplicity is the dimension of the null space of $\lambda I - A$, the number of independent eigenvectors it has. The geometric multiplicity is at least 1 and at most the algebraic one. Stacking the eigenvector equations gives $AU = U\Lambda$, and when the $n$ eigenvectors are independent, $U$ is invertible and

$$A = U \Lambda U^{-1},$$

and $A$ is called diagonalizable (Murphy, 2022, §7.4.2). A matrix with an eigenvalue whose geometric multiplicity is smaller than its algebraic one is defective and has no such factorization. The shear

$$S = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}$$

is the smallest example: it is triangular, so its eigenvalue 1 has algebraic multiplicity 2, but $S - I$ has rank 1, so there is only one eigenvector direction, $(1, 0)$. Distinct eigenvalues always give independent eigenvectors, so a matrix can only be defective when an eigenvalue repeats (a standard result, stated here without a source read for it).

A real matrix can have complex eigenvalues, in conjugate pairs. A 3D rotation by an angle $\varphi$ about a unit axis $\mathbf{a}$ has $R\mathbf{a} = \mathbf{a}$, so the axis is an eigenvector with eigenvalue 1, as Murphy (2022, §7.4.1) notes; the other two eigenvalues are $e^{\pm i\varphi}$ (shown numerically in the example).

Symmetric matrices and the spectral theorem

When $A$ is real and symmetric, all its eigenvalues are real and its eigenvectors can be chosen orthonormal (Murphy, 2022, §7.4.3), so $U$ is orthogonal and

$$A = U \Lambda U^\top = \sum_{i=1}^{n} \lambda_i\, \mathbf{u}_i \mathbf{u}_i^\top .$$

This is the spectral theorem, and every symmetric matrix is diagonalizable, repeated eigenvalues included. Each term $\lambda_i \mathbf{u}_i\mathbf{u}_i^\top$ scales the component along $\mathbf{u}_i$. Inverting is then trivial when no eigenvalue is zero, $A^{-1} = \sum_i \lambda_i^{-1} \mathbf{u}_i\mathbf{u}_i^\top$ (§7.4.3), and so is any function of $A$ that acts on the eigenvalues, such as the inverse square root used for whitening in principal component analysis. The LAPACK Users’ Guide (Anderson et al., 1999) writes the same factorization as $A = Z\Lambda Z^\top$ with the eigenvalues in ascending order, which is the order NumPy and SciPy return.

When an eigenvalue repeats, its eigenspace has an orthonormal basis, but every rotation of that basis is another one. The decomposition is unique only up to signs when the eigenvalues are distinct, and up to a rotation within each eigenspace when they are not.

Quadratic forms and definiteness

With $\mathbf{y} = U^\top\mathbf{x}$,

$$\mathbf{x}^\top A \mathbf{x} = \mathbf{y}^\top \Lambda \mathbf{y} = \sum_{i} \lambda_i y_i^2 ,$$

so a symmetric matrix is positive definite exactly when all its eigenvalues are positive and positive semidefinite when all are non-negative; with eigenvalues of both signs it is indefinite (Murphy, 2022, §7.4.3.1). For a positive definite $A$ the level sets $\mathbf{x}^\top A\mathbf{x} = r$ are ellipsoids whose axes are the eigenvectors, with semi-axes proportional to $\lambda_i^{-1/2}$ (§7.4.4); the Gaussian distribution article describes covariance ellipses this way. Matrices of the form $B^\top B$, such as covariances and structure tensors, are always positive semidefinite.

The Rayleigh quotient

The Rayleigh quotient of a symmetric $A$ at a non-zero vector $\mathbf{x}$ is

$$R(A, \mathbf{x}) = \frac{\mathbf{x}^\top A \mathbf{x}}{\mathbf{x}^\top \mathbf{x}},$$

and at an eigenvector it equals the eigenvalue (Murphy, 2022, §7.4.6). Murphy (§7.4.8) shows with a Lagrange multiplier that a maximizer or minimizer of $\mathbf{x}^\top A\mathbf{x}$ on the unit sphere must be an eigenvector. From the sum above, $R$ is a weighted average of the eigenvalues with weights $y_i^2/\lVert\mathbf{y}\rVert^2$, so (our derivation) it lies between $\lambda_{\min}$ and $\lambda_{\max}$ and reaches them at the corresponding eigenvectors. This is why the leading principal axis is the direction of maximum variance and why normalized cuts relaxes its criterion to a Rayleigh quotient minimized by an eigenvector.

The Rayleigh quotient is also a forgiving eigenvalue estimate. If a unit vector $\mathbf{x}$ makes an angle $\delta$ with the top eigenvector of a $2 \times 2$ symmetric matrix, then $R = \lambda_1\cos^2\delta + \lambda_2\sin^2\delta$, so

$$\lambda_1 - R = (\lambda_1 - \lambda_2) \sin^2\delta$$

(our derivation; in $n$ dimensions the error is likewise of second order in the angle). An eigenvector that is wrong by $10^{-3}$ radians gives an eigenvalue wrong by about $10^{-6}$ of the gap.

The generalized symmetric-definite problem

Several methods in vision solve $A\mathbf{v} = \lambda B\mathbf{v}$ with $A$ symmetric and $B$ symmetric positive definite: Fisher’s discriminant with between- and within-class scatter matrices in linear discriminant analysis and in Fisherfaces (eigenfaces), and normalized cuts with the graph Laplacian and the degree matrix. Factoring $B = LL^\top$ turns it into the ordinary symmetric problem for $L^{-1}AL^{-\top}$ with eigenvectors $L^\top\mathbf{v}$ (our derivation), so the eigenvalues are real and the eigenvectors are orthogonal in the inner product $\mathbf{v}^\top B \mathbf{w}$, the “D-orthogonality” that the normalized cuts article uses. SciPy’s scipy.linalg.eigh(a, b) solves this problem directly, with the gvd driver by default for all eigenvalues and gvx when a subset is requested (SciPy 1.16.3 source).

Properties

How far eigenvalues move

Perturb a symmetric $A$ by a symmetric $E$. By Weyl’s inequality, a standard result stated here in our words, each eigenvalue, sorted, moves by at most the spectral norm of the perturbation:

$$|\lambda_i(A + E) - \lambda_i(A)| \le \lVert E \rVert_2 .$$

Symmetric eigenvalues are therefore always well conditioned in absolute terms. The LAPACK Users’ Guide applies the same idea to rounding: the computed decomposition is nearly the exact decomposition of $A + E$ with $\lVert E\rVert_2 / \lVert A\rVert_2 \le p(n)\varepsilon$, so each computed eigenvalue is within about $p(n)\,\varepsilon\,\lVert A\rVert_2$ of a true one, where $p(n)$ is a modestly growing function of $n$ and $\varepsilon$ the machine precision. Large eigenvalues are accurate in relative terms and small ones may not be (“Further Details: Error Bounds for the Symmetric Eigenproblem”).

Non-symmetric matrices have no such guarantee. Perturbing the shear to $\begin{bmatrix} 1 & 1 \\ \epsilon & 1\end{bmatrix}$ gives eigenvalues $1 \pm \sqrt{\epsilon}$ (our derivation from the characteristic polynomial $(1-\lambda)^2 = \epsilon$), so a change of $10^{-10}$ moves the eigenvalues by $10^{-5}$. Near a defective matrix, eigenvalues are sensitive to the square root of the perturbation.

How far eigenvectors move: the gap

Eigenvectors are a different matter. The LAPACK Users’ Guide gives the approximate bound

$$\theta(\hat{\mathbf{z}}_i, \mathbf{z}_i) \lesssim \frac{p(n)\,\varepsilon\,\lVert A\rVert_2}{\operatorname{gap}_i}, \qquad \operatorname{gap}_i = \min_{j \ne i} |\lambda_i - \lambda_j|,$$

for the angle between a computed and a true unit eigenvector, valid when $p(n)\varepsilon$ is small enough; it concludes that an eigenvector whose eigenvalue is close to others may be inaccurate. The guide adds that a cluster of close eigenvalues far from the rest can still have a well-determined invariant subspace, with the same bound using the gap between the cluster and the other eigenvalues, and that computed eigenvectors are orthogonal to working precision however far each is from the truth (except, occasionally, those computed by inverse iteration for tightly clustered eigenvalues).

The same reasoning with measurement noise in place of rounding error (our reasoning) gives a rule of thumb: the angle by which a perturbation $E$ turns an eigenvector is, to first order, at most $\lVert E\rVert_2$ divided by the gap. This is a first-order statement, approximate as the guide states it, and not a theorem for large perturbations. For a $2\times 2$ symmetric matrix an exact version is easy (our derivation). Write $A = mI + D$ with $D$ traceless; $D$ has eigenvalues $\pm g/2$, where $g$ is the gap, and its dominant direction at angle $\phi$ satisfies

$$\phi = \tfrac{1}{2}\operatorname{atan2}\!\left(2A_{12},\, A_{11} - A_{22}\right), \qquad \lambda_{1,2} = m \pm \sqrt{\left(\tfrac{A_{11} - A_{22}}{2}\right)^2 + A_{12}^2} .$$

Split the perturbation the same way, $E = tI + F$. The isotropic part $tI$ shifts both eigenvalues by $t$ and turns no eigenvector. The traceless part $F$ adds a vector of length $\lVert F\rVert_2$ to one of length $g/2$ in the $(\cos 2\phi, \sin 2\phi)$ plane, so the eigenvector turns by an angle $\Delta$ with

$$\sin 2\Delta \le \frac{2\lVert F\rVert_2}{g} \qquad \text{when } 2\lVert F\rVert_2 < g,$$

with equality possible. For a $2\times 2$ matrix $\lVert E\rVert_2 = |t| + \lVert F\rVert_2$, so the first-order bound with $\lVert E\rVert_2$ can be very loose when most of the perturbation is isotropic, and, with $t = 0$, it can be exceeded at second order, since $\tfrac12\arcsin(2x) > \arcsin x$ for $0 < x \le \tfrac12$.

At a repeated eigenvalue the gap is zero and the eigenvector is not determined at all. Any vector in the eigenspace is an eigenvector, a library returns some basis, and an arbitrarily small perturbation decides which direction the computed “dominant” eigenvector takes.

Power iteration and its rate

Repeated multiplication by $A$ amplifies the component along the eigenvector of largest absolute eigenvalue. Starting from $\mathbf{v}_0$ with a non-zero component along $\mathbf{u}_1$ and normalizing at each step, $\mathbf{v}_t \propto A\mathbf{v}_{t-1}$ converges to $\mathbf{u}_1$, with the error reduced by approximately $|\lambda_2/\lambda_1|$ per iteration (Murphy, 2022, §7.4.6). The Rayleigh quotient of the iterate then estimates $\lambda_1$, with the error of second order in the angle shown above. Deflation, $A^{(2)} = A - \lambda_1\mathbf{u}_1\mathbf{u}_1^\top$, removes the found eigenpair so that the power method finds the next one (§7.4.7). For a $2\times 2$ symmetric positive semidefinite matrix the rate is exact (our derivation): in the eigenbasis the two components are multiplied by $\lambda_1^k$ and $\lambda_2^k$, so $\tan\delta_k = (\lambda_2/\lambda_1)^k \tan\delta_0$. A small gap relative to $\lambda_1$ means slow convergence. The spectral normalization of GAN discriminators estimates the largest singular value of each weight matrix with one step of power iteration.

Computing Eigenvalues

Dense symmetric matrices

LAPACK, which NumPy and SciPy call, solves the dense symmetric problem in two stages (“Symmetric Eigenproblems”). First, orthogonal Householder transformations reduce $A$ to a symmetric tridiagonal matrix $T$, $A = QTQ^\top$ (routine xSYTRD). Second, the eigenvalues and eigenvectors of $T$ are computed, and those of $A$ are $Z = QS$. LAPACK offers several algorithms for the second stage: the implicitly shifted QR algorithm (xSTEQR, used by the -EV drivers), Cuppen’s divide and conquer (xSTEDC, used by the -EVD drivers, which the guide says can be many times faster than QR for large matrices but needs more workspace), the relatively robust representation (RRR) algorithm (xSTEGR, which the guide describes as faster than all the other routines except in a few cases and as using the least workspace), and bisection with inverse iteration for a subset of eigenvalues (xSTEBZ and xSTEIN, used by the -EVX drivers). Non-symmetric matrices are first reduced to upper Hessenberg form instead (xGEHRD), the condensed form for the non-symmetric problem (“Eigenvalue Problems”).

What the libraries return

• np.linalg.eigh and np.linalg.eigvalsh call LAPACK’s _syevd (divide and conquer) and return eigenvalues in ascending order with eigenvectors as columns. They read only one triangle of the input, the lower by default (UPLO='L'), so a non-symmetric input raises no error and gives the decomposition of a different matrix, the symmetric one that the lower triangle defines (NumPy 2.3.5 source).
• scipy.linalg.eigh also reads the lower triangle by default (lower=True) and, as its notes say, does not check the input for symmetry, so that a matrix can be passed with only one triangle filled. Its default driver for the standard problem is evr (LAPACK syevr), which the notes call the optimal choice for most general cases, with ev (symmetric QR, the slowest and most robust), evd and evx as alternatives; subset_by_index and subset_by_value return only some eigenpairs (SciPy 1.16.3 source). NumPy and SciPy therefore use different drivers by default, and drivers are free to return different signs, and different bases for repeated eigenvalues (our reasoning: the decomposition fixes neither).
• np.linalg.eig calls _geev for general matrices. Its eigenvalues are not necessarily ordered, are complex unless all imaginary parts are zero, and come in conjugate pairs for real input; the eigenvector matrix “may not be of maximum rank”, that is, for a defective matrix some columns are dependent, although rounding may hide it (NumPy 2.3.5 source).
• scipy.sparse.linalg.eigsh wraps ARPACK’s implicitly restarted Lanczos method and returns a few eigenpairs of a large sparse or implicit symmetric matrix, selected by which (largest magnitude by default, or largest or smallest algebraic) or near a shift sigma (SciPy 1.16.3 source). The normalized cuts example uses it for the few eigenvectors it needs.
• skimage.feature.structure_tensor returns the three distinct entries of the windowed structure tensor in 'rc' (row, column) order by default, Arr, Arc, Acc, from SciPy’s unnormalized Sobel derivatives, with zero padding (mode='constant') by default. structure_tensor_eigenvalues returns the eigenvalues in decreasing order; for 2D tensors it uses the explicit formula $m \pm \sqrt{A_{rc}^2 + ((A_{rr} - A_{cc})/2)^2}$, and only for higher dimensions does it call eigvalsh and reverse the result (scikit-image 0.25.2, feature/corner.py).

Examples

Eigenvector accuracy and the gap, on structure tensors of gratings

The example builds images whose structure tensor is known exactly and measures how well its dominant eigenvector, the orientation of the pattern, survives noise as the gap between the two eigenvalues closes. Each $128 \times 128$ image is a sum of two whole-period sinusoidal gratings on a periodic grid: one with wave vector $(12, 5)$ cycles per image, at orientation $\theta = \operatorname{atan2}(5, 12) = 22.62°$, and a perpendicular one, $(-5, 12)$, with amplitude $a \in \{0, 0.25, 0.5, 0.75, 0.9, 1\}$ relative to the first, with random phases. Conventions: $x$ is the column index and $y$ the row index (pointing down), angles are measured from the $x$ axis toward $+y$, and the structure tensor is in $(x, y)$ order. Derivatives are computed spectrally (by the FFT), which is exact for these band-limited images, and the tensor is averaged over the whole image, one window. Over whole periods the mean of $\sin^2$ is $\tfrac12$ and the cross terms vanish, so (our derivation)

$$J_0 = \frac{\omega^2}{2}\left(\mathbf{n}_1\mathbf{n}_1^\top + a^2\,\mathbf{n}_2\mathbf{n}_2^\top\right), \qquad \omega = \frac{2\pi \cdot 13}{128},$$

with $\mathbf{n}_1 = (\cos\theta, \sin\theta)$ the dominant eigenvector, eigenvalues $\omega^2/2$ and $a^2\omega^2/2$, and gap $(\omega^2/2)(1 - a^2)$, which closes at $a = 1$. Gaussian noise of standard deviation $\sigma \in \{0.1, 0.5\}$ is added, with seeds 0–9 and the same noise field for every $a$ and $\sigma$ of a seed, so that differences between amplitudes are not masked by different noise draws. With $J$ the noisy tensor, $E = J - J_0$, and $\theta_{\text{err}}$ is the angle between the dominant eigenvector from np.linalg.eigh and $\mathbf{n}_1$.

The example tests five hypotheses and eight checks, all fixed before any code ran:

• H1: $\sin\theta_{\text{err}} \le \lVert E\rVert_2/\text{gap}$ in all 100 noisy cases with $a < 1$, the first-order bound applied to noise.
• H2: at each $\sigma$ the mean orientation error over the ten seeds increases strictly with $a$ over $a = 0$ to $0.9$; at $a = 1$ the eigenvector is arbitrary, and those rows are printed and not scored.
• H3: power iteration from $(1, 0)$ on every matrix with $a < 1$ (the 100 noisy ones and the five noise-free ones of seed 0) satisfies $\sin\delta_k / \sin\delta_0$ within a factor of 2 of $r^k$, with $r = \lambda_2/\lambda_1$ of that matrix and $\delta_k$ the angle to its own dominant eigenvector, for every $k \le 30$ at which $r^k\sin\delta_0 > 10^{-10}$; below that, double-precision rounding decides the ratio.
• H4: the closed-form angle and eigenvalues above agree with eigh to $10^{-12}$ (radians, and relative to $\lambda_{\max}$) on all 120 noisy cases and the five noise-free cases with $a < 1$.
• H5: $\lambda_1(J_0) - R(J_0, \hat{\mathbf{u}}) = \text{gap}\cdot\sin^2\theta_{\text{err}}$ to a relative $10^{-6}$ on the 100 noisy cases with $a < 1$, with the left side computed as $\hat{\mathbf{u}}^\top(\lambda_1 I - J_0)\hat{\mathbf{u}}$.
• Checks: C1, the spectral derivative against the analytic one and the noise-free tensor against $J_0$ (both to $10^{-12}$); C2, each noise field’s standard deviation within 2% of 1; C3, Weyl’s inequality on every noisy case (ratio at most $1 + 10^{-12}$); C4, the exact $2\times2$ bound $\sin 2\theta_{\text{err}} \le 2\lVert F\rVert_2/\text{gap}$ where it applies (ratio at most $1 + 10^{-9}$); C5, scikit-image’s tensor against 64 times one built from a hand-written normalized Sobel kernel, and its eigenvalues against eigvalsh reversed (both to $10^{-12}$), on one image per $a$ and $\sigma$ with sigma=2; C6, the shear under eig and eigh (eigenvectors less than $10^{-6}$ degrees apart, rank of $S - I$ equal to 1, an eigh residual of at least 0.5, and UPLO='U' giving 0 and 2); C7, the $\sqrt{\epsilon}$ shift of the perturbed shear (to a relative $10^{-6}$); and C8, a rotation by 40° about $(1, 2, 2)/3$ (eigenvalue 1 and the axis to $10^{-12}$, moduli 1 to $10^{-12}$, arguments ±40° to $10^{-10}$ degrees).

The code runs with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 in about half a second on one CPU thread.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage as ndi
from skimage.feature import structure_tensor, structure_tensor_eigenvalues
from skimage.filters import gaussian

N = 128
y, x = np.mgrid[0:N, 0:N].astype(float)          # x = column, y = row
k1, k2 = np.array([12, 5]), np.array([-5, 12])   # wave vectors, cycles per image
omega = 2 * np.pi * np.hypot(*k1) / N            # radians per pixel
theta = np.arctan2(k1[1], k1[0])                 # grating 1 orientation, 22.62 deg
n1 = np.array([np.cos(theta), np.sin(theta)])
n2 = np.array([-np.sin(theta), np.cos(theta)])
AMPS, SIGMAS, SEEDS = [0.0, 0.25, 0.5, 0.75, 0.9, 1.0], [0.1, 0.5], range(10)

def image(a, ph1, ph2):
    return (np.cos(2 * np.pi * (k1[0] * x + k1[1] * y) / N + ph1)
            + a * np.cos(2 * np.pi * (k2[0] * x + k2[1] * y) / N + ph2))

f = np.fft.fftfreq(N)                            # cycles per pixel
f[N // 2] = 0.0                                  # drop the Nyquist term
def gradient(img):                               # exact for band-limited images
    F = np.fft.fft2(img)
    gx = np.fft.ifft2(F * 2j * np.pi * f[None, :]).real
    gy = np.fft.ifft2(F * 2j * np.pi * f[:, None]).real
    return gx, gy

def tensor(img):                                 # one window: the whole image
    gx, gy = gradient(img)
    return np.array([[np.mean(gx * gx), np.mean(gx * gy)],
                     [np.mean(gx * gy), np.mean(gy * gy)]])

def angle(u, v):                                 # angle between lines, radians
    return np.arctan2(abs(u[0] * v[1] - u[1] * v[0]), abs(u @ v))

def power_iteration(J, u, k):
    for _ in range(k):
        u = J @ u
        u = u / np.linalg.norm(u)
    return u

print(f"theta = {np.degrees(theta):.4f} deg, omega^2/2 = {omega**2 / 2:.6f}")

# C1: noise-free tensor against the closed form; spectral against analytic derivative
I = image(0.5, 0.3, 1.1)
gx, _ = gradient(I)
gx_true = (-(2 * np.pi * k1[0] / N) * np.sin(2 * np.pi * (k1[0] * x + k1[1] * y) / N + 0.3)
           - 0.5 * (2 * np.pi * k2[0] / N) * np.sin(2 * np.pi * (k2[0] * x + k2[1] * y) / N + 1.1))
print(f"C1 spectral vs analytic derivative, max abs diff {np.abs(gx - gx_true).max():.2e}")
c1 = max(np.abs(tensor(image(a, 0.3, 1.1))
                - omega**2 / 2 * (np.outer(n1, n1) + a**2 * np.outer(n2, n2))).max()
         for a in AMPS) / (omega**2 / 2)
print(f"C1 noise-free J vs closed form, max rel diff {c1:.2e}")

rows, h1, h4, h4v, h5, c3, c4, h3 = [], [], [], [], [], [], [], []
noise_sd = []
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    ph1, ph2 = rng.uniform(0, 2 * np.pi, 2)
    z = rng.standard_normal((N, N))              # same noise field for every a, sigma
    noise_sd.append(z.std())
    for a in AMPS:
        J0 = omega**2 / 2 * (np.outer(n1, n1) + a**2 * np.outer(n2, n2))
        lam0 = np.linalg.eigvalsh(J0)            # ascending
        gap = lam0[1] - lam0[0]
        cases = [(0.0, tensor(image(a, ph1, ph2)))] if seed == 0 else []
        cases += [(s, tensor(image(a, ph1, ph2) + s * z)) for s in SIGMAS]
        for s, J in cases:
            w, V = np.linalg.eigh(J)
            u = V[:, 1]                          # dominant eigenvector
            phi = 0.5 * np.arctan2(2 * J[0, 1], J[0, 0] - J[1, 1])
            m, d = (J[0, 0] + J[1, 1]) / 2, np.hypot((J[0, 0] - J[1, 1]) / 2, J[0, 1])
            if s > 0 or a < 1:                   # H4: closed form against eigh
                h4.append(angle(u, np.array([np.cos(phi), np.sin(phi)])))
                h4v.append(np.abs(np.array([m - d, m + d]) - w).max() / w[1])
            if s == 0:
                continue
            E = J - J0
            normE = np.linalg.norm(E, 2)
            err = angle(u, n1)
            c3.append(np.abs(w - lam0).max() / normE)
            rows.append((s, a, seed, np.degrees(err)))
            if a == 1:
                continue
            h1.append((s, a, np.sin(err) / (normE / gap)))
            F = E - np.trace(E) / 2 * np.eye(2)  # traceless part of E
            normF = np.linalg.norm(F, 2)
            if 2 * normF < gap:
                c4.append(np.sin(2 * err) / (2 * normF / gap))
            rq_err = u @ (lam0[1] * np.eye(2) - J0) @ u  # lambda1 - R(J0, u)
            h5.append(abs(rq_err / (gap * np.sin(err) ** 2) - 1))
        for s, J in cases:                       # H3: power iteration on each matrix
            if a == 1:
                continue
            w, V = np.linalg.eigh(J)
            r, u0 = w[0] / w[1], np.array([1.0, 0.0])
            s0 = np.sin(angle(u0, V[:, 1]))
            for k in range(1, 31):
                if r**k * s0 > 1e-10:
                    sk = np.sin(angle(power_iteration(J, u0, k), V[:, 1]))
                    h3.append(sk / s0 / r**k)

print(f"C2 noise field sd over seeds: {min(noise_sd):.4f} to {max(noise_sd):.4f}")
rows = np.array(rows)
print("orientation error of the dominant eigenvector (deg), mean and range over 10 seeds")
for s in SIGMAS:
    for a in AMPS:
        e = rows[(rows[:, 0] == s) & (rows[:, 1] == a), 3]
        gap = omega**2 / 2 * (1 - a**2)
        r1 = max(v for (ss, aa, v) in h1 if ss == s and aa == a) if a < 1 else float("nan")
        print(f"sigma {s:.1f}  a {a:.2f}  gap {gap:.5f}  mean {e.mean():.4e}  "
              f"range {e.min():.4e} to {e.max():.4e}  max sin/(|E|/gap) {r1:.4f}")
print(f"H1 cases {len(h1)}, max ratio {max(v for *_, v in h1):.6f}")
print(f"C3 Weyl: max |lambda_i - lambda0_i| / |E|_2 = {max(c3):.6f}, 1 - max = {1 - max(c3):.2e}")
print(f"C4 sin(2 err) / (2|F|/gap): cases {len(c4)}, max {max(c4):.6f}")
print(f"H5 Rayleigh: max |(lambda1 - R) / (gap sin^2 err) - 1| = {max(h5):.2e}")
print(f"H3 power iteration: {len(h3)} (case, k) pairs, "
      f"ratio range {min(h3):.6f} to {max(h3):.6f}")
print(f"H4 closed form vs eigh: {len(h4)} cases, max angle {max(h4):.2e} rad, "
      f"max eigenvalue diff {max(h4v):.2e}")
J1 = tensor(image(1.0, 0.3, 1.1))
w, V = np.linalg.eigh(J1)
print(f"a = 1 noise-free: eigenvalues {w}, eigh dominant vector {V[:, 1]}")

# C5: scikit-image's structure tensor, scale factor and eigenvalue order
d = np.array([-0.5, 0.0, 0.5])[None, :] * np.array([0.25, 0.5, 0.25])[:, None]
c5a, c5b, c5c = 0.0, 0.0, True
for s in SIGMAS:
    rng = np.random.default_rng(0)
    ph1, ph2 = rng.uniform(0, 2 * np.pi, 2)
    z = rng.standard_normal((N, N))
    for a in AMPS:
        img = image(a, ph1, ph2) + s * z
        Arr, Arc, Acc = structure_tensor(img, sigma=2)        # defaults: 'constant', 'rc'
        dx = ndi.correlate(img, d, mode="constant")           # normalized Sobel
        dy = ndi.correlate(img, d.T, mode="constant")
        mine = [gaussian(p, sigma=2, mode="constant") for p in (dy * dy, dx * dy, dx * dx)]
        scale = max(np.abs(Arr).max(), np.abs(Acc).max())
        c5a = max(c5a, max(np.abs(A - 64 * B).max() for A, B in zip((Arr, Arc, Acc), mine)) / scale)
        ev = structure_tensor_eigenvalues([Arr, Arc, Acc])
        M = np.stack([np.stack([Arr, Arc], -1), np.stack([Arc, Acc], -1)], -2)
        ref = np.linalg.eigvalsh(M)[..., ::-1]
        c5b = max(c5b, np.abs(np.moveaxis(ev, 0, -1) - ref).max() / ev[0].max())
        c5c = c5c and bool(np.all(ev[0] >= ev[1]))
print(f"C5 skimage tensor vs 64 x normalized-Sobel tensor, max rel diff {c5a:.2e}")
print(f"C5 skimage eigenvalues vs eigvalsh reversed, max rel diff {c5b:.2e}, decreasing: {c5c}")

# C6-C8: matrices that are not symmetric
S = np.array([[1.0, 1.0], [0.0, 1.0]])
w, V = np.linalg.eig(S)
print(f"C6 shear eig: eigenvalues {w}, angle between eigenvectors "
      f"{np.degrees(angle(V[:, 0], V[:, 1])):.2e} deg, rank(S - I) = {np.linalg.matrix_rank(S - np.eye(2))}")
for uplo in ("L", "U"):
    w, V = np.linalg.eigh(S, UPLO=uplo)
    res = max(np.linalg.norm(S @ V[:, i] - w[i] * V[:, i]) for i in range(2))
    print(f"C6 eigh(S, UPLO='{uplo}'): eigenvalues {w}, max residual |Sv - lambda v| {res:.4f}")
eps = 1e-10
w = np.linalg.eig(np.array([[1.0, 1.0], [eps, 1.0]])).eigenvalues
print(f"C7 eps {eps:.0e}: eigenvalue shift {np.abs(w - 1).max():.6e}, sqrt(eps) {np.sqrt(eps):.6e}, "
      f"rel diff {abs(np.abs(w - 1).max() / np.sqrt(eps) - 1):.1e}")
axis = np.array([1.0, 2.0, 2.0]) / 3
K = np.array([[0, -axis[2], axis[1]], [axis[2], 0, -axis[0]], [-axis[1], axis[0], 0]])
t = np.radians(40)
R = np.eye(3) + np.sin(t) * K + (1 - np.cos(t)) * K @ K          # Rodrigues' formula
w, V = np.linalg.eig(R)
i = np.argmin(np.abs(w - 1))
v = V[:, i].real
print(f"C8 rotation eig: {np.round(w, 6)}")
print(f"C8 |lambda - 1| {abs(w[i] - 1):.1e}, |sin(eigvec, axis)| {np.linalg.norm(np.cross(v, axis)):.1e}, "
      f"max |modulus - 1| {np.abs(np.abs(np.delete(w, i)) - 1).max():.1e}, "
      f"max ||arg| - 40 deg| {np.abs(np.abs(np.degrees(np.angle(np.delete(w, i)))) - 40).max():.1e}")

Output:

theta = 22.6199 deg, omega^2/2 = 0.203609
C1 spectral vs analytic derivative, max abs diff 1.74e-14
C1 noise-free J vs closed form, max rel diff 4.09e-16
C2 noise field sd over seeds: 0.9949 to 1.0085
orientation error of the dominant eigenvector (deg), mean and range over 10 seeds
sigma 0.1  a 0.00  gap 0.20361  mean 5.3921e-02  range 1.6959e-02 to 1.0285e-01  max sin/(|E|/gap) 0.0113
sigma 0.1  a 0.25  gap 0.19088  mean 5.7525e-02  range 1.8081e-02 to 1.0972e-01  max sin/(|E|/gap) 0.0113
sigma 0.1  a 0.50  gap 0.15271  mean 7.1910e-02  range 2.2594e-02 to 1.3712e-01  max sin/(|E|/gap) 0.0113
sigma 0.1  a 0.75  gap 0.08908  mean 1.2325e-01  range 3.8727e-02 to 2.3467e-01  max sin/(|E|/gap) 0.0112
sigma 0.1  a 0.90  gap 0.03869  mean 2.8362e-01  range 8.9186e-02 to 5.3737e-01  max sin/(|E|/gap) 0.0112
sigma 0.1  a 1.00  gap 0.00000  mean 3.6468e+01  range 2.2165e+00 to 7.4193e+01  max sin/(|E|/gap) nan
sigma 0.5  a 0.00  gap 0.20361  mean 1.3117e+00  range 4.4368e-01 to 2.3529e+00  max sin/(|E|/gap) 0.0103
sigma 0.5  a 0.25  gap 0.19088  mean 1.3982e+00  range 4.7318e-01 to 2.4997e+00  max sin/(|E|/gap) 0.0103
sigma 0.5  a 0.50  gap 0.15271  mean 1.7383e+00  range 5.9618e-01 to 3.0665e+00  max sin/(|E|/gap) 0.0101
sigma 0.5  a 0.75  gap 0.08908  mean 2.9204e+00  range 9.1726e-01 to 5.2380e+00  max sin/(|E|/gap) 0.0100
sigma 0.5  a 0.90  gap 0.03869  mean 6.3378e+00  range 1.5685e+00 to 1.1613e+01  max sin/(|E|/gap) 0.0096
sigma 0.5  a 1.00  gap 0.00000  mean 3.4803e+01  range 3.4481e+00 to 7.8444e+01  max sin/(|E|/gap) nan
H1 cases 100, max ratio 0.011268
C3 Weyl: max |lambda_i - lambda0_i| / |E|_2 = 1.000000, 1 - max = -8.88e-16
C4 sin(2 err) / (2|F|/gap): cases 100, max 0.999867
H5 Rayleigh: max |(lambda1 - R) / (gap sin^2 err) - 1| = 9.18e-10
H3 power iteration: 2629 (case, k) pairs, ratio range 1.001487 to 1.177005
H4 closed form vs eigh: 125 cases, max angle 1.67e-16 rad, max eigenvalue diff 2.21e-16
a = 1 noise-free: eigenvalues [0.20360878 0.20360878], eigh dominant vector [0. 1.]
C5 skimage tensor vs 64 x normalized-Sobel tensor, max rel diff 3.31e-16
C5 skimage eigenvalues vs eigvalsh reversed, max rel diff 2.68e-16, decreasing: True
C6 shear eig: eigenvalues [1. 1.], angle between eigenvectors 1.27e-14 deg, rank(S - I) = 1
C6 eigh(S, UPLO='L'): eigenvalues [1. 1.], max residual |Sv - lambda v| 1.0000
C6 eigh(S, UPLO='U'): eigenvalues [0. 2.], max residual |Sv - lambda v| 0.7071
C7 eps 1e-10: eigenvalue shift 1.000000e-05, sqrt(eps) 1.000000e-05, rel diff 6.6e-12
C8 rotation eig: [0.766044+0.642788j 0.766044-0.642788j 1.      +0.j      ]
C8 |lambda - 1| 1.1e-16, |sin(eigvec, axis)| 1.6e-16, max |modulus - 1| 2.2e-16, max ||arg| - 40 deg| 0.0e+00

The data. C1 held: the spectral derivative matches the analytic one to $1.74 \times 10^{-14}$ and the noise-free tensor matches $J_0$ to $4.09\times10^{-16}$ of $\omega^2/2$, so the ground truth is exact. C2 held: the noise fields’ standard deviations lie between 0.9949 and 1.0085.

Error against the gap. H1 held: the largest ratio $\sin\theta_{\text{err}}/(\lVert E\rVert_2/\text{gap})$ is 0.011268, so the first-order bound is met in all 100 cases with about 90 times to spare. The construction explains the slack (our derivation, which holds for every case): for a $2\times2$ matrix $\lVert E\rVert_2 = |t| + \lVert F\rVert_2$, and the isotropic part $t$ shifts both eigenvalues without turning the eigenvector, so a bound written with $\lVert E\rVert_2$ counts a part of the perturbation that cannot rotate anything. C4 held: the exact $2\times2$ bound, written with the traceless part alone, applies in all 100 cases and its largest ratio is 0.999867, close to equality. H2 held at both noise levels: the mean error grows strictly with $a$, from 0.053921° to 0.28362° at $\sigma = 0.1$ and from 1.3117° to 6.3378° at $\sigma = 0.5$, as the gap falls from 0.20361 to 0.03869. At $a = 1$, printed and not scored, the gap is zero: the mean errors are 36.468° and 34.803°, and the ranges, from 2.2165° to 74.193° and from 3.4481° to 78.444°, run over most of the possible 0° to 90°. The noise-free tensor at $a = 1$ is a multiple of the identity, both eigenvalues are 0.20360878, and eigh returns the basis vector $(0, 1)$ as its “dominant” eigenvector, a direction that has nothing to do with either grating.

Eigenvalues. C3 held: no eigenvalue moved by more than $\lVert E\rVert_2$; the largest ratio exceeds 1 by $8.88\times10^{-16}$, rounding, within the bar of $10^{-12}$. The ratio reaches 1 because at $a = 1$ the noise-free tensor is a multiple of the identity, so the eigenvalues of $J$ are those of $J_0$ shifted by exactly $t \pm \lVert F\rVert_2$ (our derivation). H5 held: the Rayleigh quotient of the true tensor at the noisy eigenvector falls short of $\lambda_1$ by $\text{gap}\cdot\sin^2\theta_{\text{err}}$, to a relative $9.18\times10^{-10}$: the eigenvalue estimate’s error is the square of the eigenvector’s.

Power iteration and the closed form. H3 held: over 2629 (matrix, $k$) pairs, the ratio of $\sin\delta_k/\sin\delta_0$ to $r^k$ lies between 1.001487 and 1.177005. Because the rate is exact for $2\times2$ matrices (our derivation, under Power iteration above), this tests the implementation and where rounding takes over, not an asymptotic rate: the ratio equals $\cos\delta_k/\cos\delta_0$, which lies between 1 and $1/\cos\delta_0$. H4 held: the closed-form angle and eigenvalues agree with eigh to $1.67\times10^{-16}$ radians and $2.21\times10^{-16}$ in all 125 cases, the 20 noisy cases at $a = 1$ included.

Libraries. C5 held: scikit-image’s tensor equals 64 times the tensor from the normalized Sobel kernel ($[-1, 0, 1]/2$ across, $[1, 2, 1]/4$ along) to $3.31\times10^{-16}$, with Arr the $yy$ entry and Acc the $xx$ entry, and its eigenvalues are in decreasing order and equal eigvalsh’s reversed to $2.68\times10^{-16}$. The factor 64 is the square of 8, the gain of SciPy’s unnormalized Sobel kernel. The scale matters whenever a threshold on eigenvalues is carried between implementations; the order matters whenever code written for eigh’s ascending order is pointed at scikit-image’s output.

Matrices that are not symmetric. C6 held. eig returns eigenvalue 1 twice with two eigenvector columns $1.27\times10^{-14}$ degrees apart, numerically the same direction, and $S - I$ has rank 1: the shear is defective. eigh with the default UPLO='L' reads the lower triangle, which is the identity, and silently returns eigenvalues 1 and 1 with the identity’s eigenvectors, one of which misses $S\mathbf{v} = \lambda\mathbf{v}$ by a residual of 1. With UPLO='U' it reads $\begin{bmatrix}1 & 1\\1 & 1\end{bmatrix}$ and returns 0 and 2, eigenvalues that $S$ does not have. C7 held: a perturbation of $10^{-10}$ in the lower-left entry moves the eigenvalues by $1.000000\times10^{-5}$, equal to $\sqrt{\epsilon}$ to a relative $6.6\times10^{-12}$. C8 held: the rotation’s real eigenvalue is 1 to $1.1\times10^{-16}$ with an eigenvector along the axis to $1.6\times10^{-16}$, and the other two are $0.766044 \pm 0.642788i$, of modulus 1 and argument ±40°.

Process note. All five hypotheses and eight checks, with their bars, were written in a file before any script ran. A planning script, identical to the code above except for three printed lines, was then run once to compute the values a correct implementation reaches; no bar was changed after it. Because it printed C3 exactly at its limit (1.000000) and printed C8’s moduli and arguments at NumPy’s default precision, which cannot decide bars of $10^{-12}$ and $10^{-10}$, three print-only lines were added before the pilot: C3’s “1 − max” and C8’s largest modulus and argument errors, and with them C7’s relative difference, which the planning script had shown only as two seven-digit values. No bar, set-up or measure changed with them. The plan for the example listed nine departures from the batch plan’s description, all fixed before the pilot: (1) spectral derivatives and a single window covering the whole image, as the exact route for whole-period gratings; (2) the open parameters (image size, wave vectors, the two noise levels, random phases, seeds 0–9 and one noise field per seed shared across $a$ and $\sigma$); (3) H2 made precise as strictly increasing seed means over $a < 1$, with $a = 1$ printed only; (4) H3’s measure fixed as $\sin\delta_k/\sin\delta_0$ against $r^k$ from $(1, 0)$, scored only above $10^{-10}$, on the noisy and the seed-0 noise-free matrices; (5) H4’s cases set to all noisy cases and the noise-free ones with $a < 1$, with eigenvalues compared too; (6) a bar added for the Rayleigh-quotient measure, H5; (7) checks C1–C4, C7 and C8 added, and UPLO='U' added to C6; (8) the scikit-image comparison made concrete as the scale test against a hand-built kernel and the order test against eigvalsh; (9) the three print-precision changes above. The pilot’s output was identical to the first run’s, and the first run gave identical output when run twice; no setting, measure or printed value changed after the pilot.

Other examples on the site

The aperture problem uses similar data, a grating and a perpendicular plaid, but prints the structure tensor’s eigenvalues to show that one grating leaves it rank-deficient and only the normal flow recoverable; it does not measure eigenvector accuracy. Gabor filter estimates the orientation of noisy gratings from the energies of a filter bank, a different estimator. Eigenproblems in their own settings appear in the examples of principal component analysis (covariance eigenvectors against the SVD, whitening), spectral clustering (three graph Laplacians and the eigengap), Harris corner detector (corner scores from the structure tensor) and linear algebra.

Common Misconceptions

• “Every matrix has a full set of eigenvectors.” Defective matrices, such as a shear, do not, and real matrices such as rotations can have complex eigenvalues. The spectral theorem is a statement about symmetric matrices.
• “eigh checks that its input is symmetric.” NumPy’s and SciPy’s eigh read one triangle and ignore the other. A non-symmetric input gives an answer, with no warning, for a different matrix.
• “The eigenvector is a property of the data.” Only up to sign, and only when its eigenvalue is separated from the others. At a repeated eigenvalue any vector in the eigenspace is an eigenvector, and the one a library returns is a property of the algorithm; with noise it is decided by the noise, as the $a = 1$ rows show. Comparing eigenvectors from two runs, two libraries or two drivers needs a sign convention, and for clusters a comparison of subspaces.
• “Close eigenvalues are hard to compute.” For symmetric matrices the eigenvalues themselves move at most by the size of the perturbation; it is the individual eigenvectors of close eigenvalues that are poorly determined, while the subspace they span can remain well determined.
• “Eigenvalues come sorted from largest to smallest.” eigh returns them ascending, eig in no particular order, scikit-image’s structure tensor eigenvalues descending, and Murphy (2022, §7.4.1) calls it standard to sort them with the largest magnitude first. Check the order before indexing.
• “Singular values and eigenvalues are the same thing.” They coincide for symmetric positive semidefinite matrices. In general the singular values of $A$ are the square roots of the eigenvalues of $A^\top A$, as the singular value decomposition article explains, and they are the better tool for rectangular and non-symmetric matrices.

Where It Is Used

• Principal component analysis and whitening: eigenvectors of a covariance matrix as principal axes, eigenvalues as variances (principal component analysis, eigenfaces).
• The structure tensor: corners, edges and flat regions from its two eigenvalues (Harris corner detector, feature tracking), and the rank-deficient tensor of the aperture problem and Lucas–Kanade flow.
• Hessian eigenvalues: blob and tube shape from their ratio in lesion detection; minima, maxima and saddle points and the condition number of a cost function in optimization and gradient descent.
• Spectral graph methods: the Laplacian eigenvectors of spectral clustering and the generalized eigenproblem of normalized cuts.
• Discriminant directions: the generalized eigenproblem of linear discriminant analysis.
• Geometry: surface normals as the eigenvector of smallest eigenvalue of a neighbourhood covariance in a point cloud; Horn’s rotation as the top eigenvector of a symmetric $4\times4$ matrix in the Kabsch algorithm; the axis of a rotation as its eigenvector for eigenvalue 1; covariance ellipses of the Gaussian distribution.
• Kernels: kernel PCA and the Nyström approximation in kernel methods.