10. Symmetric Matrices and Quadratic Forms

Symmetric matrices are the best-behaved matrices in linear algebra. Their eigenvalues are real, their eigenvectors can be chosen orthonormal, and they are always diagonalizable by a rotation of coordinates. This chapter proves that fact, the spectral theorem, and extends it to Hermitian and normal matrices. It then applies the theorem to quadratic forms, such as energies, variances, and second-order Taylor terms, and finishes with positive definite matrices, the Cholesky factorization, and quadratic optimization. These ideas underpin least squares, statistics, convex optimization, and the SVD of Chapter 11.

<a id="101-spectral-theory"></a>

10.1 Spectral Theory

This subchapter proves that the symmetric matrices are exactly the real matrices with an orthonormal basis of eigenvectors. It extends the result to Hermitian and normal matrices, and rewrites every symmetric matrix as a weighted sum of perpendicular projections.

Properties of Symmetric Matrices

Why it matters. Symmetric matrices appear wherever a matrix records a relationship that does not depend on order (Section 4.2). Examples are covariance matrices, Hessians of smooth functions, stiffness matrices, graph Laplacians, and every Gram matrix $A^TA$. Chapter 8 showed that a general matrix can have complex eigenvalues, too few eigenvectors, or eigenvectors at awkward angles. For a symmetric matrix none of these problems can occur. This topic proves the facts that rule them out. The next topic assembles them into the spectral theorem.

Analogy. Every square matrix splits into a symmetric part and a skew-symmetric part (Section 4.2). Near the identity, the skew-symmetric part acts like an infinitesimal rotation (Section 9.4), and the symmetric part acts like a pure stretch. A symmetric matrix is all stretch and no twist. Think of a rubber sheet pulled along two perpendicular directions: points on those directions stay on their lines, and they are the eigenvectors.

The key identity: symmetric means self-adjoint. A real square matrix $A$ is symmetric if and only if

$$(A\mathbf{x}) \cdot \mathbf{y} = \mathbf{x} \cdot (A\mathbf{y}) \qquad \text{for all } \mathbf{x}, \mathbf{y} \in \mathbb{R}^n.$$

The left side is $\mathbf{x}^TA^T\mathbf{y}$ and the right side is $\mathbf{x}^TA\mathbf{y}$. They agree for all vectors exactly when $A^T = A$; take $\mathbf{x} = \mathbf{e}_i$ and $\mathbf{y} = \mathbf{e}_j$ to compare entries. A matrix with this property is called self-adjoint. In words, $A$ can be moved from one side of the dot product to the other. Every property below follows from this one line.

Property 1: the eigenvalues are real. Let $A\mathbf{z} = \lambda\mathbf{z}$ with $\mathbf{z} \in \mathbb{C}^n$ nonzero, allowing complex values for now. A real symmetric matrix is Hermitian, so $\mathbf{z}^*A\mathbf{z}$ is real (Section 4.2). Then

$$\mathbf{z}^A\mathbf{z} = \lambda,\mathbf{z}^\mathbf{z}, \qquad \mathbf{z}^*\mathbf{z} = \lVert \mathbf{z} \rVert^2 > 0,$$

so $\lambda = \mathbf{z}^A\mathbf{z} / \mathbf{z}^\mathbf{z}$ is a ratio of real numbers. Because $\lambda$ is real, $A - \lambda I$ is a real singular matrix and has a real nonzero null vector. So the eigenvectors can be chosen real as well.

Property 2: eigenvectors for distinct eigenvalues are orthogonal. Let $A\mathbf{v}_1 = \lambda_1\mathbf{v}_1$ and $A\mathbf{v}_2 = \lambda_2\mathbf{v}_2$ with $\lambda_1 \ne \lambda_2$. Moving $A$ across the dot product,

$$\lambda_1(\mathbf{v}_1 \cdot \mathbf{v}_2) = (A\mathbf{v}_1) \cdot \mathbf{v}_2 = \mathbf{v}_1 \cdot (A\mathbf{v}_2) = \lambda_2(\mathbf{v}_1 \cdot \mathbf{v}_2),$$

so $(\lambda_1 - \lambda_2)(\mathbf{v}_1 \cdot \mathbf{v}_2) = 0$, which forces $\mathbf{v}_1 \cdot \mathbf{v}_2 = 0$. For a general matrix, eigenvectors for distinct eigenvalues are only guaranteed to be linearly independent (Section 8.2). Symmetry upgrades "independent" to "perpendicular".

Property 3: orthogonal complements of invariant subspaces are invariant. Suppose a subspace $W$ satisfies $A W \subseteq W$ (Section 8.3). Let $\mathbf{u} \in W^\perp$ and $\mathbf{w} \in W$. Then $(A\mathbf{u}) \cdot \mathbf{w} = \mathbf{u} \cdot (A\mathbf{w}) = 0$, because $A\mathbf{w} \in W$. So $A\mathbf{u} \in W^\perp$. Once one eigenvector has been found, $A$ therefore maps the hyperplane perpendicular to it into itself. Repeating the argument inside that hyperplane drives the proof of the spectral theorem.

Property 4: the column space and the null space are orthogonal complements. For any matrix, $N(A) = C(A^T)^\perp$ (Section 9.2). If $A = A^T$, this reads $N(A) = C(A)^\perp$, so

$$\mathbb{R}^n = C(A) \oplus N(A), \qquad C(A) \perp N(A).$$

For a non-symmetric matrix the two subspaces can overlap or sit at any angle. For example, $\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$ has $C(A) = N(A) = \operatorname{span}(\mathbf{e}_1)$.

Summary of the contrast.

	General real square matrix	Real symmetric matrix
Eigenvalues	may be complex	always real
Eigenvectors for distinct eigenvalues	linearly independent	orthogonal
Diagonalizable	not always (defective matrices)	always, by an orthogonal matrix (next topic)
Geometric vs. algebraic multiplicity	geometric $\le$ algebraic	always equal
$C(A)$ and $N(A)$	unrelated	orthogonal complements
Rank	at least the number of nonzero eigenvalues	exactly the number of nonzero eigenvalues
Sensitivity of eigenvalues to perturbations	can be enormous	at most the size of the perturbation (Section 10.3)

The last three rows follow from the spectral theorem and from Weyl's inequality, which are proved later in this chapter.

Worked example. For $A = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$, the characteristic polynomial is $(2 - \lambda)^2 - 1 = (\lambda - 3)(\lambda - 1)$. The eigenvectors are $(1, 1)$ for $\lambda = 3$ and $(1, -1)$ for $\lambda = 1$, and their dot product is $0$. Compare $B = \begin{bmatrix} 1 & 2 \ 0 & 3 \end{bmatrix}$. It has real eigenvalues $1$ and $3$, but its eigenvectors $(1, 0)$ and $(1, 1)$ meet at $45°$.

import numpy as np

A = np.array([[2., 1.], [1., 2.]])           # symmetric
B = np.array([[1., 2.], [0., 3.]])           # upper triangular, not symmetric
R = np.array([[0., -1.], [1., 0.]])          # rotation by 90 degrees

for name, M in [("A", A), ("B", B), ("R", R)]:
    vals, vecs = np.linalg.eig(M)
    cos = abs(np.vdot(vecs[:, 0], vecs[:, 1]))    # eig returns unit vectors, so this is |cos(angle)|
    print(name, np.real_if_close(vals), round(cos, 4))
# A [3. 1.] 0.0                  real eigenvalues, perpendicular eigenvectors
# B [1. 3.] 0.7071               real eigenvalues, eigenvectors at 45 degrees
# R [0.+1.j 0.-1.j] 0.0          complex eigenvalues (R is normal: see Normal Matrices)

# Self-adjointness: (Ax).y = x.(Ay) holds for every x, y exactly when A is symmetric
rng = np.random.default_rng(0)
x, y = rng.standard_normal((2, 2))
print(np.isclose((A @ x) @ y, x @ (A @ y)), np.isclose((B @ x) @ y, x @ (B @ y)))   # True False

# Singular symmetric matrix: the column space is perpendicular to the null space
S = np.array([[1., 2.], [2., 4.]])
null = np.array([2., -1.])                    # S @ null = 0
print(S @ null, S[:, 0] @ null)               # [0. 0.] 0.0
C = np.array([[1., 2.], [0., 0.]])            # not symmetric, same null space
print(C @ null, C[:, 0] @ null)               # [0. 0.] 2.0: column (1, 0) is not perpendicular to (2, -1)

# Random symmetric matrices always have real spectra; random general matrices usually do not
M = rng.standard_normal((200, 200))
print(np.abs(np.linalg.eigvals(M + M.T).imag).max(), np.abs(np.linalg.eigvals(M).imag).max() > 1)   # 0.0 True

Key takeaways

• A real matrix is symmetric exactly when it is self-adjoint: $(A\mathbf{x}) \cdot \mathbf{y} = \mathbf{x} \cdot (A\mathbf{y})$ for all $\mathbf{x}, \mathbf{y}$.
• The eigenvalues of a real symmetric matrix are real, and its eigenvectors can be chosen real.
• Eigenvectors belonging to distinct eigenvalues are orthogonal, not merely independent.
• If $W$ is invariant under a symmetric $A$, so is $W^\perp$; this is the engine of the spectral theorem.
• For symmetric $A$, the column space and the null space are orthogonal complements.

🧪 Practice

1. Find the eigenvalues and eigenvectors of $\begin{bmatrix} 3 & 4 \ 4 & -3 \end{bmatrix}$ and verify that the eigenvectors are orthogonal.
2. Show directly that $\begin{bmatrix} a & b \ b & c \end{bmatrix}$ has real eigenvalues by computing the discriminant of its characteristic polynomial. When does it have a repeated eigenvalue?
3. Prove that if $A$ is symmetric and $A^2 = 0$, then $A = 0$. Then give a non-symmetric counterexample.
4. Interview: A colleague computes the eigenvalues of a large covariance matrix with np.linalg.eig and gets complex numbers with imaginary parts around 1e-17. Explain where they come from and what they should do instead. Hint: rounding during the computation of the matrix can leave it slightly non-symmetric, and eig is a general solver that cannot assume symmetry; symmetrize with (S + S.T) / 2 and use eigh, which returns real eigenvalues by construction.

Spectral Theorem for Real Symmetric Matrices

Why it matters. Diagonalizing a matrix means finding coordinates in which it acts by independent scalings (Section 8.2). In general, these coordinates may not exist, and when they do, the eigenvector matrix $P$ can be badly conditioned. The spectral theorem says that for symmetric matrices, the best case always happens: the coordinates exist, and the axes are perpendicular. It is the single most used theorem about matrices in applied mathematics. PCA, vibration analysis, quantum mechanics, and spectral graph theory are all applications of it.

Theorem (Spectral theorem, real case). For a real $n \times n$ matrix $A$, the following are equivalent:

1. $A$ is symmetric.
2. $\mathbb{R}^n$ has an orthonormal basis consisting of eigenvectors of $A$.
3. $A = QDQ^T$ for some orthogonal matrix $Q$ ($Q^TQ = I$) and some real diagonal matrix $D$.

The columns of $Q$ are the orthonormal eigenvectors, and the diagonal entries of $D$ are the corresponding eigenvalues.

Proof. (2) $\Leftrightarrow$ (3) is the definition of diagonalization (Section 8.2) with $P = Q$ orthogonal, using $Q^{-1} = Q^T$ (Section 9.4).

(3) $\Rightarrow$ (1): $(QDQ^T)^T = (Q^T)^TD^TQ^T = QDQ^T$, since a diagonal matrix is symmetric.

(1) $\Rightarrow$ (2), by induction on $n$. The case $n = 1$ is trivial. By Property 1 of the previous topic, $A$ has a real eigenvalue $\lambda_1$ with a real unit eigenvector $\mathbf{q}_1$. The line $W = \operatorname{span}(\mathbf{q}_1)$ is invariant, so $W^\perp$ is invariant by Property 3. Choose any orthonormal basis of $W^\perp$ and put it in the columns of an $n \times (n - 1)$ matrix $Q'$. Then $Q_1 = [\mathbf{q}_1 \mid Q']$ is orthogonal, and

$$Q_1^TAQ_1 = \begin{bmatrix} \lambda_1 & \mathbf{0}^T \ \mathbf{0} & A' \end{bmatrix}, \qquad A' = Q'^TAQ'.$$

The off-diagonal blocks are zero because $A\mathbf{q}_1 = \lambda_1\mathbf{q}_1$ is perpendicular to $W^\perp$, and because the matrix is symmetric. $A'$ is symmetric and of size $n - 1$. By induction, it has an orthonormal eigenbasis $\mathbf{y}_2, \dots, \mathbf{y}_n$ of $\mathbb{R}^{n-1}$. The vectors $\mathbf{q}_1, Q'\mathbf{y}_2, \dots, Q'\mathbf{y}_n$ are then an orthonormal eigenbasis of $\mathbb{R}^n$. $\blacksquare$

A second proof, via Schur. Because the eigenvalues of $A$ are real, the real Schur form (Section 8.3) has no $2 \times 2$ blocks: $A = QTQ^T$ with $Q$ orthogonal and $T$ upper triangular. But $T = Q^TAQ$ is symmetric, and a symmetric upper triangular matrix is diagonal. $\blacksquare$

The "only if" direction is a fingerprint. No non-symmetric real matrix has an orthonormal basis of real eigenvectors. If you observe perpendicular eigenvectors, the matrix is symmetric, and if a model's matrix ought to have perpendicular modes, it ought to be symmetric.

Consequences.

• No defective symmetric matrices. Every symmetric matrix is diagonalizable, so algebraic and geometric multiplicities always agree, and the Jordan form is diagonal (Section 8.3).
• Orthogonal eigenspaces. The eigenspaces for distinct eigenvalues are mutually perpendicular, and together they fill $\mathbb{R}^n$: $\mathbb{R}^n = E_{\mu_1} \oplus \cdots \oplus E_{\mu_k}$, an orthogonal direct sum (Section 9.2).
• Rank. $\operatorname{rank}A = \operatorname{rank}D$ is the number of nonzero eigenvalues, counted with multiplicity.
• Trace and determinant. $\operatorname{tr}A = \sum\lambda_i$ and $\det A = \prod\lambda_i$, as for every matrix (Section 8.3).

Geometric meaning. Read $A = QDQ^T$ from right to left as three steps applied to $\mathbf{x}$:

   In the standard basis, A mixes x and y:        In the eigenbasis q1, q2, A only rescales:

   [x']   [2 1] [x]                               [c1']   [3 0] [c1]
   [y'] = [1 2] [y]                               [c2'] = [0 1] [c2]

   A = Q D Q^T, read right to left:

      x  ---Q^T--->  (c1, c2) = (q1.x, q2.x)  ---D--->  (3 c1, 1 c2)  ---Q--->  A x
         rotate into the                      stretch along             rotate back
         eigenbasis                           each axis

   Effect on the unit circle: an ellipse with semi-axis 3 along q1 = (1, 1)/sqrt(2)
   and semi-axis 1 along q2 = (1, -1)/sqrt(2).

Every symmetric matrix is a rotation (or reflection), followed by stretches along perpendicular axes, followed by the inverse rotation. A negative eigenvalue stretches and flips that axis.

In NumPy: eigh, not eig. np.linalg.eigh is the symmetric/Hermitian eigensolver. It returns real eigenvalues in ascending order and an orthonormal matrix of eigenvectors, even when eigenvalues are repeated, and it is faster than eig. It reads only one triangle of the input (the lower one by default), so it silently gives wrong answers for a matrix that is not symmetric.

import numpy as np

A = np.array([[4., 1., 1.],
              [1., 4., 1.],
              [1., 1., 4.]])                  # eigenvalues 6, 3, 3

w, Q = np.linalg.eigh(A)                      # ascending real eigenvalues, orthonormal eigenvectors
print(w)                                      # [3. 3. 6.]
print(np.allclose(Q.T @ Q, np.eye(3)))        # True: Q is orthogonal
print(np.allclose(Q @ np.diag(w) @ Q.T, A))   # True: A = Q D Q^T

# eig does not know A is symmetric: inside the repeated eigenspace its vectors need not be orthogonal
w2, P = np.linalg.eig(A)
print(np.round(np.real_if_close(w2), 6))      # [6. 3. 3.]
print(np.round(np.abs(P.T @ P), 3))           # off-diagonal 0.189 (the exact value depends on the LAPACK build)
# [[1.    0.    0.   ]
#  [0.    1.    0.189]
#  [0.    0.189 1.   ]]

# Pitfall: eigh reads only the lower triangle and silently assumes symmetry
N = np.array([[1., 0.], [5., 2.]])            # lower triangular, eigenvalues 1 and 2
print(np.linalg.eigh(N)[0], np.real_if_close(np.linalg.eigvals(N)))
# [-3.52493781  6.52493781] [2. 1.]          eigh solved the symmetric matrix [[1, 5], [5, 2]] instead

Key takeaways

• A real matrix is orthogonally diagonalizable, $A = QDQ^T$, if and only if it is symmetric.
• The proof peels off one eigenvector at a time, using the invariance of orthogonal complements; alternatively, a symmetric Schur form must be diagonal.
• Symmetric matrices are never defective, their eigenspaces are mutually orthogonal, and their rank equals the number of nonzero eigenvalues.
• Geometrically, a symmetric matrix stretches space along perpendicular axes, the eigenvectors, by factors equal to the eigenvalues.
• Use np.linalg.eigh for symmetric matrices: it is faster, guarantees real output and orthonormal eigenvectors, but trusts its input to be symmetric.

🧪 Practice

1. Orthogonally diagonalize $A = \begin{bmatrix} 1 & 2 \ 2 & -2 \end{bmatrix}$: find $Q$ and $D$ with $A = QDQ^T$, and check the trace and determinant.
2. Prove that a symmetric matrix with a single eigenvalue $\lambda$ (of multiplicity $n$) must equal $\lambda I$. Why does this fail for $\begin{bmatrix} \lambda & 1 \ 0 & \lambda \end{bmatrix}$?
3. Give a non-symmetric $2 \times 2$ matrix with two distinct real eigenvalues. Show that its eigenvectors are independent but not orthogonal, and confirm that this is forced by the theorem.
4. Interview: Why does NumPy provide eigh when eig already handles every square matrix? When is using eigh a bug? Hint: think about guarantees (real eigenvalues, orthonormal eigenvectors even for repeated eigenvalues, sorted output), speed (reduction to a symmetric tridiagonal matrix), and the fact that eigh never checks symmetry.

Orthogonal Diagonalization

Why it matters. The spectral theorem says that $Q$ and $D$ exist. This topic shows how to find them by hand and what to do with them afterwards. Once $A = QDQ^T$ is known, powers, inverses, square roots, and exponentials of $A$ reduce to operations on the numbers on the diagonal of $D$. Because $Q^{-1} = Q^T$, nothing has to be inverted, and the computation does not amplify rounding errors.

Procedure. To orthogonally diagonalize a symmetric matrix $A$:

1. Eigenvalues. Find the roots of $\det(A - \lambda I)$, with multiplicities. They are all real.
2. Eigenspaces. For each distinct eigenvalue $\mu$, find a basis of $N(A - \mu I)$. Its dimension equals the multiplicity of $\mu$.
3. Orthonormalize inside each eigenspace. Apply Gram-Schmidt (Section 9.3) to each eigenspace basis separately. Vectors from different eigenspaces are already orthogonal (Property 2), so this step is only needed for repeated eigenvalues. For simple eigenvalues it reduces to normalizing.
4. Assemble. Place the orthonormal vectors in the columns of $Q$, and the matching eigenvalues, in the same order, on the diagonal of $D$. Then $A = QDQ^T$.

A common mistake is to run Gram-Schmidt across all the eigenvectors at once. It produces an orthonormal set, but the mixed vectors are no longer eigenvectors.

Worked example. Let

$$A = \begin{bmatrix} 2 & 1 & 1 \ 1 & 2 & 1 \ 1 & 1 & 2 \end{bmatrix} = I + J,$$

where $J$ is the all-ones matrix. $J$ has rank 1, with $J(1, 1, 1) = 3(1, 1, 1)$ and $J\mathbf{x} = \mathbf{0}$ for every $\mathbf{x}$ perpendicular to $(1, 1, 1)$. Adding $I$ adds 1 to each eigenvalue.

1. Eigenvalues: $4$ (simple) and $1$ (multiplicity 2).
2. $E_4 = \operatorname{span}{(1, 1, 1)}$. $E_1 = N(A - I) = {x + y + z = 0}$, with basis $\mathbf{v}_1 = (1, -1, 0)$, $\mathbf{v}_2 = (1, 0, -1)$.
3. Normalize $(1, 1, 1)$ to $\mathbf{q}_1 = \frac{1}{\sqrt{3}}(1, 1, 1)$. Gram-Schmidt in $E_1$: $\mathbf{q}_2 = \frac{1}{\sqrt{2}}(1, -1, 0)$, and

$$\mathbf{v}_2 - (\mathbf{v}_2 \cdot \mathbf{q}_2)\mathbf{q}_2 = (1, 0, -1) - \tfrac12(1, -1, 0) = \left(\tfrac12, \tfrac12, -1\right), \qquad \mathbf{q}_3 = \tfrac{1}{\sqrt{6}}(1, 1, -2).$$

4. Assemble:

$$Q = \begin{bmatrix} 1/\sqrt{3} & 1/\sqrt{2} & 1/\sqrt{6} \ 1/\sqrt{3} & -1/\sqrt{2} & 1/\sqrt{6} \ 1/\sqrt{3} & 0 & -2/\sqrt{6} \end{bmatrix}, \qquad D = \begin{bmatrix} 4 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{bmatrix}.$$

The choice of $\mathbf{q}_2, \mathbf{q}_3$ is not unique. Any orthonormal basis of the plane $x + y + z = 0$ works, and so does flipping the sign of any column. The matrices $D$ and $\mathbf{q}_2\mathbf{q}_2^T + \mathbf{q}_3\mathbf{q}_3^T$ (the projection onto $E_1$) are unique (see Spectral Decomposition and Rank-One Sums).

Functions of a symmetric matrix. If $A = QDQ^T$, then $A^k = QD^kQ^T$, because the inner factors $Q^TQ = I$ cancel. More generally, for any function $f$ defined on the eigenvalues,

$$f(A) = Q,f(D),Q^T = Q \begin{bmatrix} f(\lambda_1) & & \ & \ddots & \ & & f(\lambda_n) \end{bmatrix} Q^T.$$

Function	Formula	Requirement	Use
Power $A^k$	$QD^kQ^T$	none	Markov-like iterations, $A^{100}$ in $O(n^3)$
Inverse $A^{-1}$	$QD^{-1}Q^T$	all $\lambda_i \ne 0$	analysis; solve systems with Cholesky instead (Section 10.3)
Square root $A^{1/2}$	$QD^{1/2}Q^T$	all $\lambda_i \ge 0$	whitening, matrix means, the polar decomposition
Exponential $e^{tA}$	$Qe^{tD}Q^T$	none	$\mathbf{x}' = A\mathbf{x}$, heat diffusion on graphs

For a positive semidefinite $A$ (all $\lambda_i \ge 0$, Section 10.3), $A^{1/2} = QD^{1/2}Q^T$ is the unique symmetric positive semidefinite square root. Other square roots exist, but they are not symmetric with nonnegative eigenvalues.

Orthogonal versus general diagonalization.

	$A = PDP^{-1}$ (general)	$A = QDQ^T$ (symmetric)
Exists when	$A$ has $n$ independent eigenvectors	always, for symmetric $A$
Inverting the eigenvector matrix	requires solving with $P$	free: $Q^{-1} = Q^T$
Condition number of the eigenvector matrix	can be arbitrarily large	exactly 1
Eigenvector coordinates of $\mathbf{x}$	solve $P\mathbf{c} = \mathbf{x}$	dot products $c_i = \mathbf{q}_i \cdot \mathbf{x}$ (Section 9.2)

import numpy as np

A = np.array([[2., 1., 1.],
              [1., 2., 1.],
              [1., 1., 2.]])

# Step 2: eigenspace bases found by hand
v1 = np.array([1., 1., 1.])                   # lambda = 4
basis_1 = np.array([[1., -1., 0.],            # lambda = 1: the plane x + y + z = 0
                    [1., 0., -1.]]).T

# Step 3: orthonormalize inside each eigenspace (QR performs Gram-Schmidt)
q1 = v1 / np.linalg.norm(v1)
Q1, _ = np.linalg.qr(basis_1)                 # orthonormal basis of the plane (QR may flip signs)
Q = np.column_stack([q1, Q1])
D = np.diag([4., 1., 1.])
print(np.round(Q * np.sqrt(6), 4) + 0.0)      # columns scaled by sqrt(6) to expose the exact entries
# [[ 1.4142 -1.7321 -1.    ]                    column 1: sqrt(2)(1, 1, 1)  = sqrt(6) (1, 1, 1)/sqrt(3)
#  [ 1.4142  1.7321 -1.    ]                    column 2: sqrt(3)(-1, 1, 0) = sqrt(6) (-1, 1, 0)/sqrt(2)
#  [ 1.4142  0.      2.    ]]                   column 3: (-1, -1, 2)       = sqrt(6) (-1, -1, 2)/sqrt(6)
print(np.allclose(Q.T @ Q, np.eye(3)), np.allclose(Q @ D @ Q.T, A))   # True True

# Functions of A through the diagonal: f(A) = Q f(D) Q^T
def sym_fun(A, f):
    w, V = np.linalg.eigh(A)
    return V @ np.diag(f(w)) @ V.T

root = sym_fun(A, np.sqrt)
print(np.allclose(root @ root, A), np.allclose(root, root.T))             # True True
print(np.allclose(sym_fun(A, lambda w: 1 / w), np.linalg.inv(A)))         # True
print(np.allclose(sym_fun(A, lambda w: w**10), np.linalg.matrix_power(A, 10)))   # True
print(np.round(sym_fun(A, lambda w: w**10)[0], 1))   # [349526. 349525. 349525.] = ((4^10 + 2)/3, (4^10 - 1)/3, ...)

Key takeaways

• Orthogonal diagonalization: find the eigenvalues, a basis of each eigenspace, orthonormalize within each eigenspace, and assemble $Q$ and $D$.
• Eigenvectors from different eigenspaces are automatically orthogonal; Gram-Schmidt is needed only inside repeated eigenspaces.
• $Q$ is never unique (signs, and bases of repeated eigenspaces), but $D$ and the eigenspaces are.
• $f(A) = Qf(D)Q^T$ turns powers, inverses, square roots, and exponentials into scalar operations on the eigenvalues.
• A positive semidefinite matrix has a unique positive semidefinite square root $QD^{1/2}Q^T$.

🧪 Practice

1. Orthogonally diagonalize $A = \begin{bmatrix} 0 & 2 & 2 \ 2 & 0 & 2 \ 2 & 2 & 0 \end{bmatrix}$. (Write $A = 2(J - I)$ to guess the eigenvalues.)
2. Compute the positive semidefinite square root of $\begin{bmatrix} 5 & 4 \ 4 & 5 \end{bmatrix}$ and verify your answer by squaring it.
3. Show that if $A = QDQ^T$ is invertible, then $A^{-1}$ is orthogonally diagonalized by the same $Q$. Then show that $e^A$ is symmetric with positive eigenvalues for every symmetric $A$.
4. Interview: Data with covariance matrix $\Sigma$ must be whitened: transformed by a matrix $W$ so that the new covariance $W\Sigma W^T$ is the identity. Give one choice of $W$ based on orthogonal diagonalization, and describe the family of all valid choices. Hint: $\Sigma = QDQ^T$ gives $W = \Sigma^{-1/2} = QD^{-1/2}Q^T$ (ZCA whitening); any $W$ with $W^TW = \Sigma^{-1}$ works, for example $D^{-1/2}Q^T$ (PCA whitening) or the inverse of a Cholesky factor (Section 10.3).

Spectral Theorem for Hermitian Matrices

Why it matters. Quantum mechanics, signal processing, and the analysis of real matrices with complex eigenvalues all happen in $\mathbb{C}^n$. There, the role of symmetric matrices is played by Hermitian matrices ($A^* = A$, Section 4.2), and the role of orthogonal matrices by unitary ones ($U^*U = I$, Section 9.4). The spectral theorem carries over almost word for word. Its physical content is that observables have real measured values and perfectly distinguishable outcomes.

Theorem (Spectral theorem, Hermitian case). A complex $n \times n$ matrix $A$ is Hermitian if and only if

$$A = UDU^*$$

for some unitary $U$ and some real diagonal $D$. Equivalently, $\mathbb{C}^n$ has an orthonormal basis (for the inner product $\langle \mathbf{u}, \mathbf{v} \rangle = \mathbf{v}^*\mathbf{u}$, Section 9.1) of eigenvectors of $A$, and all eigenvalues are real.

Proof. By Schur's theorem (Section 8.3), $A = UTU^*$ with $U$ unitary and $T$ upper triangular. If $A$ is Hermitian, then

$$T^* = (U^AU)^ = U^*A^*U = U^*AU = T.$$

$T$ is upper triangular, and $T^$ is lower triangular, so both are diagonal. The diagonal of $T = T^$ satisfies $t_{ii} = \overline{t_{ii}}$, so it is real. Conversely, $(UDU^)^ = U\overline{D}U^* = UDU^*$ when $D$ is real. $\blacksquare$

The properties from the real case carry over with conjugates inserted. $A$ is self-adjoint, $\langle A\mathbf{u}, \mathbf{v} \rangle = \langle \mathbf{u}, A\mathbf{v} \rangle$. Its eigenvalues are real, because $\lambda = \mathbf{z}^A\mathbf{z} / \mathbf{z}^\mathbf{z}$. Eigenvectors for distinct eigenvalues are orthogonal in the complex inner product. The real theorem is the special case in which all entries are real.

	Real symmetric	Complex Hermitian
Condition	$A^T = A$	$A^* = A$
Eigenvalues	real	real
Eigenvector matrix	orthogonal, $Q^TQ = I$	unitary, $U^*U = I$
Factorization	$A = QDQ^T$	$A = UDU^*$
Quadratic form	$\mathbf{x}^TA\mathbf{x} \in \mathbb{R}$	$\mathbf{z}^*A\mathbf{z} \in \mathbb{R}$
NumPy	np.linalg.eigh	np.linalg.eigh (same function)

Worked example. Let $H = \begin{bmatrix} 2 & 1 - i \ 1 + i & 3 \end{bmatrix}$. Then $\operatorname{tr}H = 5$ and $\det H = 6 - (1 - i)(1 + i) = 6 - 2 = 4$, so the eigenvalues are $1$ and $4$, both real. For $\lambda = 4$, $(H - 4I)\mathbf{z} = \mathbf{0}$ gives $-2z_1 + (1 - i)z_2 = 0$, so $\mathbf{u}_2 = \frac{1}{\sqrt{6}}(1 - i, 2)$. For $\lambda = 1$, $\mathbf{u}_1 = \frac{1}{\sqrt{3}}(1 - i, -1)$. Check: $\langle \mathbf{u}_1, \mathbf{u}_2 \rangle \propto (1 - i)\overline{(1 - i)} + (-1)(2) = 2 - 2 = 0$.

Quantum measurement. In quantum mechanics a state is a unit vector $\boldsymbol{\psi} \in \mathbb{C}^n$, and an observable is a Hermitian matrix $H = \sum_k \lambda_k\mathbf{u}_k\mathbf{u}_k^*$. A measurement returns one of the eigenvalues $\lambda_k$, with probability $p_k = \lvert \langle \boldsymbol{\psi}, \mathbf{u}_k \rangle \rvert^2$. The probabilities add up to 1 because the eigenvectors form an orthonormal basis (Parseval, Section 9.2). The expected value is $\sum p_k\lambda_k = \boldsymbol{\psi}^*H\boldsymbol{\psi}$. Real eigenvalues give real outcomes, and orthogonal eigenvectors make distinct outcomes mutually exclusive. That is why observables must be Hermitian.

Relatives and a pitfall.

• Skew-Hermitian matrices ($K^* = -K$) satisfy $K = iH$ with $H = -iK$ Hermitian. So they are unitarily diagonalizable with purely imaginary eigenvalues. A real skew-symmetric matrix is a special case: its eigenvalues are $0$ and pairs $\pm i\theta$, and it can be diagonalized over $\mathbb{C}$ but not over $\mathbb{R}$.
• Complex symmetric matrices ($A^T = A$ with complex entries) are not Hermitian, and none of the theory applies. The matrix $C = \begin{bmatrix} 1 & i \ i & -1 \end{bmatrix}$ is complex symmetric and satisfies $C^2 = 0$. So it is a nonzero nilpotent matrix, and it is defective. Over $\mathbb{C}$, the useful symmetry is $A^* = A$, not $A^T = A$.

import numpy as np

H = np.array([[2., 1 - 1j],
              [1 + 1j, 3.]])                  # Hermitian: real diagonal, conjugate off-diagonal pair
w, U = np.linalg.eigh(H)
print(w)                                      # [1. 4.]: real, ascending
print(np.allclose(U.conj().T @ U, np.eye(2)), np.allclose(U @ np.diag(w) @ U.conj().T, H))   # True True

# Quantum measurement: outcomes are eigenvalues, probabilities |<psi, u_k>|^2
psi = np.array([1., 1j]) / np.sqrt(2)         # a normalized state
probs = np.abs(U.conj().T @ psi) ** 2         # squared coordinates in the eigenbasis
print(np.round(probs, 4), np.isclose(probs.sum(), 1))       # [0.1667 0.8333] True
print(np.isclose(probs @ w, (psi.conj() @ H @ psi).real))   # True: expected value = sum of p_k * lambda_k

# A complex symmetric (not Hermitian) matrix can be defective
C = np.array([[1, 1j], [1j, -1]])
print(np.allclose(C, C.T), np.allclose(C, C.conj().T))       # True False: symmetric, but not Hermitian
print(np.allclose(C @ C, 0))                                 # True: C^2 = 0 although C != 0, so C is nilpotent

# Real skew-symmetric matrices have purely imaginary eigenvalues
K = np.array([[0., -3.], [3., 0.]])
print(np.linalg.eigvals(K), np.linalg.eigvalsh(1j * K))      # [0.+3.j 0.-3.j] [-3.  3.]: iK is Hermitian

Key takeaways

• A complex matrix is Hermitian if and only if it is unitarily diagonalizable with real eigenvalues: $A = UDU^*$.
• The proof is one line from Schur: a Hermitian triangular matrix is diagonal with a real diagonal.
• Eigenvectors of a Hermitian matrix for distinct eigenvalues are orthogonal in the complex inner product $\mathbf{v}^*\mathbf{u}$.
• Quantum observables are Hermitian because measured values must be real and distinct outcomes must be orthogonal.
• Skew-Hermitian matrices are $i$ times Hermitian ones; complex symmetric matrices ($A^T = A$) have no spectral theorem.

🧪 Practice

1. Unitarily diagonalize the Pauli matrix $\sigma_y = \begin{bmatrix} 0 & -i \ i & 0 \end{bmatrix}$.
2. Prove that the eigenvalues of a skew-Hermitian matrix are purely imaginary, and that a skew-Hermitian matrix is unitarily diagonalizable.
3. Write a Hermitian matrix as $H = S + iK$ with $S$, $K$ real. Show that $S$ is symmetric, $K$ is skew-symmetric, and the real $2n \times 2n$ matrix $\begin{bmatrix} S & -K \ K & S \end{bmatrix}$ is symmetric with the same eigenvalues as $H$, each appearing twice.
4. Interview: Why does quantum mechanics represent observables by Hermitian matrices rather than by arbitrary diagonalizable matrices with real eigenvalues? Hint: a diagonalizable matrix with real eigenvalues gives real outcomes, but the measurement probabilities $\lvert \langle \boldsymbol{\psi}, \mathbf{u}_k \rangle \rvert^2$ only sum to 1 for every state when the eigenvectors are orthonormal.

Normal Matrices

Why it matters. Hermitian matrices are not the only ones that can be diagonalized by a unitary matrix. Rotations, unitary matrices, skew-symmetric matrices, and circulant matrices also have orthonormal eigenbases, even though their eigenvalues are complex. Normal matrices are exactly the matrices with an orthonormal eigenbasis. For them, eigenvalues tell the whole story: powers, norms, and sensitivity to perturbations are all controlled by the eigenvalues alone. For non-normal matrices this is false, and the difference matters in numerical work and dynamics.

Definition. A complex square matrix $A$ is normal if it commutes with its conjugate transpose:

$$AA^* = A^*A.$$

For real matrices this reads $AA^T = A^TA$.

Normal	Not normal
Hermitian ($A^* = A$) and real symmetric	$\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$ (a shear, defective)
Skew-Hermitian and real skew-symmetric	$\begin{bmatrix} 1 & 2 \ 0 & 3 \end{bmatrix}$ (diagonalizable, but not unitarily)
Unitary and real orthogonal	most upper triangular matrices
Diagonal matrices	nonzero nilpotent matrices
Circulant matrices

Theorem (Spectral theorem for normal matrices). A complex square matrix $A$ is unitarily diagonalizable, $A = UDU^*$ with $U$ unitary and $D$ diagonal (possibly complex), if and only if $A$ is normal.

Proof. ($\Rightarrow$) If $A = UDU^$, then $AA^ = UD\overline{D}U^*$ and $A^A = U\overline{D}DU^$. These are equal because diagonal matrices commute.

($\Leftarrow$) By Schur, $A = UTU^$ with $T$ upper triangular. Normality transfers to $T$: $TT^ = U^*AA^*U = U^*A^*AU = T^T$. A normal upper triangular matrix is diagonal (Section 8.3). Compare the $(1, 1)$ entries: $(TT^){11} = \sum_j \lvert t{1j} \rvert^2$, the squared length of the first row, while $(T^*T){11} = \lvert t{11} \rvert^2$, since the first column has one nonzero entry. So $t_{12} = \cdots = t_{1n} = 0$, and the same argument applies to the remaining block. $\blacksquare$

Equivalent characterizations. For a square matrix $A$, each of the following is equivalent to normality.

• $\lVert A\mathbf{x} \rVert = \lVert A^\mathbf{x} \rVert$ for every $\mathbf{x}$. (The Hermitian matrix $A^A - AA^$ has $\mathbf{x}^(A^A - AA^)\mathbf{x} = 0$ for all $\mathbf{x}$, which forces it to be zero.)
• Every eigenvector of $A$ is an eigenvector of $A^$, with the conjugate eigenvalue: $A\mathbf{x} = \lambda\mathbf{x} \Rightarrow A^\mathbf{x} = \bar\lambda\mathbf{x}$.
• $A = H_1 + iH_2$ with $H_1$, $H_2$ Hermitian and $H_1H_2 = H_2H_1$.

Where the eigenvalues lie decides the type. Within the normal matrices, $A = UDU^*$ inherits each property of $D$:

Normal $A$ with eigenvalues ...	is exactly ...
real	Hermitian
purely imaginary	skew-Hermitian
on the unit circle, $\lvert \lambda \rvert = 1$	unitary
real and $\ge 0$ ($> 0$)	Hermitian positive semidefinite (definite), Section 10.3

For example, $(UDU^)^ = U\overline{D}U^*$ equals $A$ exactly when $\overline{D} = D$. Outside the normal class this breaks down: $\begin{bmatrix} 1 & 2 \ 0 & 3 \end{bmatrix}$ has real eigenvalues but is not symmetric.

Every unitary matrix is $e^{iH}$. A unitary matrix is normal with eigenvalues $e^{i\varphi_k}$, so $U = V\operatorname{diag}(e^{i\varphi_1}, \dots, e^{i\varphi_n})V^$ with $V$ unitary. Set $H = V\operatorname{diag}(\varphi_1, \dots, \varphi_n)V^$, which is Hermitian because the $\varphi_k$ are real. Then $e^{iH} = V\operatorname{diag}(e^{i\varphi_k})V^* = U$. This completes the claim made in Section 9.4. The angles $\varphi_k$ are determined only up to multiples of $2\pi$, so $H$ is not unique.

Real normal matrices. A real normal matrix may have complex eigenvalues, so it cannot always be diagonalized over $\mathbb{R}$. It can be orthogonally block-diagonalized: $A = Q,\operatorname{diag}(B_1, \dots, B_k),Q^T$ with $Q$ real orthogonal, where each block is either a real $1 \times 1$ eigenvalue or a $2 \times 2$ scaled rotation $\begin{bmatrix} a & -b \ b & a \end{bmatrix}$ carrying the pair $a \pm ib$. For a rotation of $\mathbb{R}^3$ this gives one block $1$ (the axis) and one rotation block (the plane of rotation).

Why normality matters numerically. For a normal matrix, perturbing $A$ by $E$ moves each eigenvalue by at most $\lVert E \rVert_2$ (the Bauer-Fike theorem with a perfectly conditioned eigenvector matrix; Chapter 11). For a non-normal matrix, the same perturbation can move eigenvalues by far more. Perturbing the $2 \times 2$ Jordan block by $10^{-10}$ moves its eigenvalues by $10^{-5}$. Likewise, for normal $A$, $\lVert A^k \rVert_2 = \rho(A)^k$ exactly, while non-normal matrices can show large transient growth of $A^k$ even when all $\lvert \lambda \rvert < 1$.

import numpy as np
from scipy.linalg import schur, expm

def is_normal(A):
    return np.allclose(A @ A.conj().T, A.conj().T @ A)

theta = 0.7
examples = {
    "symmetric":      np.array([[2., 1.], [1., 2.]]),
    "skew-symmetric": np.array([[0., -2.], [2., 0.]]),
    "rotation":       np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]),
    "circulant":      np.array([[1., 2., 3.], [3., 1., 2.], [2., 3., 1.]]),
    "shear":          np.array([[1., 1.], [0., 1.]]),
    "triangular":     np.array([[1., 2.], [0., 3.]]),
}
for name, A in examples.items():
    print(f"{name:15s}", is_normal(A))       # True for the first four, False for shear and triangular

# Normal => the Schur form is diagonal, so the Schur vectors are an orthonormal eigenbasis
T, Z = schur(examples["circulant"], output="complex")
print(np.allclose(T, np.diag(np.diag(T))), np.round(np.diag(T), 4) + 0)
# True [ 6. +0.j    -1.5+0.866j -1.5-0.866j]

# Every unitary matrix is exp(iH) with H Hermitian
rng = np.random.default_rng(1)
U, _ = np.linalg.qr(rng.standard_normal((3, 3)) + 1j * rng.standard_normal((3, 3)))   # a random unitary
T, Z = schur(U, output="complex")             # diagonal, entries e^{i phi_k} on the unit circle
print(np.round(np.abs(np.diag(T)), 12))       # [1. 1. 1.]
H = Z @ np.diag(np.angle(np.diag(T))) @ Z.conj().T
print(np.allclose(H, H.conj().T), np.allclose(expm(1j * H), U))   # True True

# Eigenvalue sensitivity: the same perturbation of size 1e-10 applied to a normal and a non-normal matrix
E = np.array([[0., 0.], [1e-10, 0.]])
N = np.diag([1., 2.])                          # normal
J = np.array([[1., 1.], [0., 1.]])             # Jordan block: non-normal
print(np.abs(np.sort(np.linalg.eigvals(N + E).real) - [1., 2.]).max() <= 1e-10)   # True: moved by <= 1e-10
print(np.linalg.eigvals(J + E).real - 1)       # [ 1.e-05 -1.e-05]: moved 100,000 times farther

Key takeaways

• $A$ is normal when $AA^* = A^*A$; this includes Hermitian, skew-Hermitian, unitary, orthogonal, diagonal, and circulant matrices.
• A matrix is unitarily diagonalizable if and only if it is normal (Schur plus "normal triangular implies diagonal").
• Among normal matrices, the location of the eigenvalues determines the type: real means Hermitian, imaginary means skew-Hermitian, unit modulus means unitary.
• Every unitary matrix equals $e^{iH}$ for some Hermitian $H$.
• Eigenvalues of normal matrices are insensitive to perturbations; for non-normal matrices they can be extremely sensitive.

🧪 Practice

1. Decide which matrices are normal: $\begin{bmatrix} 1 & i \ i & 1 \end{bmatrix}$, $\begin{bmatrix} 1 & 1 \ -1 & 1 \end{bmatrix}$, $\begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix}$. For each normal one, find its eigenvalues and say whether it is Hermitian, unitary, or neither.
2. Prove that if $A$ is normal and $A\mathbf{x} = \lambda\mathbf{x}$, then $A^\mathbf{x} = \bar\lambda\mathbf{x}$. Use the fact that $B = A - \lambda I$ is normal and $\lVert B\mathbf{x} \rVert = \lVert B^\mathbf{x} \rVert$.
3. Show that a normal matrix with real eigenvalues is Hermitian, and that a normal matrix whose eigenvalues all satisfy $\lvert \lambda \rvert = 1$ is unitary.
4. Interview: A discrete-time system $\mathbf{x}_{k+1} = A\mathbf{x}_k$ has all eigenvalues inside the unit disk, yet simulations show $\lVert \mathbf{x}_k \rVert$ growing by a factor of 1000 before eventually decaying. How is that possible, and why can it never happen when $A$ is normal? Hint: for normal $A$, $\lVert A^k \rVert_2 = \rho(A)^k$; for non-normal $A$, nearly parallel eigenvectors allow large cancellations in $\mathbf{x}_0$ that unwind over time (transient growth).

Spectral Decomposition and Rank-One Sums

Why it matters. $A = QDQ^T$ is a product of three matrices. Multiplying it out gives a sum: a symmetric matrix is a weighted sum of perpendicular projections, with the eigenvalues as weights. This view explains what $A$ does to a vector, gives functions of $A$ without coordinates, and shows how to compress $A$ by keeping only its largest terms. Low-rank approximation, PCA, and the SVD of Chapter 11 grow from this idea.

The rank-one expansion. Expanding $QDQ^T$ column by row (Section 4.1),

$$A = \begin{bmatrix} \mathbf{q}_1 & \cdots & \mathbf{q}_n \end{bmatrix}\begin{bmatrix} \lambda_1 & & \ & \ddots & \ & & \lambda_n \end{bmatrix}\begin{bmatrix} \mathbf{q}_1^T \ \vdots \ \mathbf{q}n^T \end{bmatrix} = \sum{i=1}^n \lambda_i,\mathbf{q}_i\mathbf{q}_i^T.$$

Each $\mathbf{q}_i\mathbf{q}_i^T$ is the orthogonal projection onto the line through $\mathbf{q}_i$ (Section 9.2), a symmetric rank-one matrix. So

$$A\mathbf{x} = \sum_{i=1}^n \lambda_i,(\mathbf{q}_i \cdot \mathbf{x}),\mathbf{q}_i:$$

                                  x
                                  |
         +------------------------+------------------------+
         |                        |                        |
   (q1 . x) q1              (q2 . x) q2       ...    (qn . x) qn      1. split x along the
         |                        |                        |             orthonormal eigenvectors
     x lambda1                x lambda2       ...     x lambdan       2. scale each piece by
         |                        |                        |             its eigenvalue
         +------------------------+------------------------+
                                  |
                                 A x                                  3. add the pieces back up

Grouping by eigenvalue: the spectral decomposition. Let $\mu_1, \dots, \mu_k$ be the distinct eigenvalues, and let $P_j$ be the orthogonal projection onto the eigenspace $E_{\mu_j}$. If $\mathbf{q}i$ runs over an orthonormal basis of $E{\mu_j}$, then $P_j = \sum \mathbf{q}_i\mathbf{q}_i^T$. Then

$$A = \sum_{j=1}^k \mu_j P_j, \qquad P_j^2 = P_j = P_j^T, \qquad P_iP_j = 0 \ (i \ne j), \qquad \sum_{j=1}^k P_j = I.$$

This is the decomposition of Section 9.2 that the spectral theorem promised. The last identity, called a resolution of the identity, says that the eigenspaces split $\mathbb{R}^n$ into perpendicular pieces. The individual terms $\mathbf{q}_i\mathbf{q}_i^T$ depend on which eigenvectors were chosen, but the projectors $P_j$ do not. So the spectral decomposition is unique up to the order of the terms.

Functional calculus. Because $P_iP_j = 0$ and $P_j^2 = P_j$, squaring the sum gives $A^2 = \sum \mu_j^2P_j$, and in general

$$f(A) = \sum_{j=1}^k f(\mu_j),P_j,$$

which covers $A^{-1}$, $A^{1/2}$, $e^{tA}$, and so on. The projectors are themselves polynomials in $A$ (Lagrange interpolation):

$$P_j = \prod_{i \ne j} \frac{A - \mu_iI}{\mu_j - \mu_i}.$$

This formula finds $P_j$ without computing any eigenvectors, which is handy for $2 \times 2$ and $3 \times 3$ matrices by hand.

Worked example. For $A = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$ with eigenvalues $3$ and $1$:

$$P_3 = \frac{A - I}{3 - 1} = \frac12\begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}, \qquad P_1 = \frac{A - 3I}{1 - 3} = \frac12\begin{bmatrix} 1 & -1 \ -1 & 1 \end{bmatrix},$$

$$A = 3P_3 + 1P_1, \qquad A^{-1} = \tfrac13P_3 + P_1 = \frac13\begin{bmatrix} 2 & -1 \ -1 & 2 \end{bmatrix}.$$

Truncation: the best low-rank approximation. Order the terms by $\lvert \lambda_i \rvert$, largest first, and keep the first $k$:

$$A_k = \sum_{i=1}^k \lambda_i,\mathbf{q}i\mathbf{q}i^T, \qquad \lVert A - A_k \rVert_2 = \lvert \lambda{k+1} \rvert, \qquad \lVert A - A_k \rVert_F = \sqrt{\textstyle\sum{i > k}\lambda_i^2}.$$

No rank-$k$ matrix does better in either norm. This is the Eckart-Young theorem, proved for all matrices with the SVD in Section 11.2. $A_k$ needs only $k(n + 1)$ numbers instead of $n^2$, and $A_k\mathbf{x}$ costs $O(nk)$ instead of $O(n^2)$. For a covariance matrix, the terms are the principal components. $\lambda_i$ is the variance along the direction $\mathbf{q}_i$, and keeping the largest terms keeps most of the variance (Section 11.2).

import numpy as np

A = np.array([[2., 1.], [1., 2.]])
w, Q = np.linalg.eigh(A)                      # w = [1, 3]
terms = [lam * np.outer(q, q) for lam, q in zip(w, Q.T)]   # lambda_i q_i q_i^T
print(np.round(terms[0], 4).tolist(), np.round(terms[1], 4).tolist())
# [[0.5, -0.5], [-0.5, 0.5]] [[1.5, 1.5], [1.5, 1.5]]     1 * P_1 and 3 * P_3
print(np.allclose(sum(terms), A))             # True

# The projector onto an eigenspace does not depend on the basis chosen inside it
B = np.array([[2., 1., 1.], [1., 2., 1.], [1., 1., 2.]])   # eigenvalues 1, 1, 4
w, Q = np.linalg.eigh(B)
P1 = Q[:, :2] @ Q[:, :2].T                    # projector onto the eigenvalue-1 plane
G = np.array([[0.6, -0.8], [0.8, 0.6]])       # rotate the basis inside the plane
Q_alt = Q[:, :2] @ G
print(np.allclose(Q_alt @ Q_alt.T, P1), np.allclose(P1, np.eye(3) - np.ones((3, 3)) / 3))   # True True

# Truncating the sum: keep the k terms with the largest |lambda|
rng = np.random.default_rng(4)
M = rng.standard_normal((6, 6))
S = M + M.T
w, Q = np.linalg.eigh(S)
order = np.argsort(-np.abs(w))                # largest |lambda| first
k = 2
idx = order[:k]
S_k = Q[:, idx] @ np.diag(w[idx]) @ Q[:, idx].T
print(np.isclose(np.linalg.norm(S - S_k, 2), abs(w[order[k]])))                       # True: spectral-norm error
print(np.isclose(np.linalg.norm(S - S_k, 'fro'), np.sqrt(np.sum(w[order[k:]]**2))))   # True: Frobenius error

Key takeaways

• Every symmetric matrix is a sum of rank-one projections weighted by eigenvalues: $A = \sum\lambda_i\mathbf{q}_i\mathbf{q}_i^T$.
• Grouping by distinct eigenvalues gives $A = \sum\mu_jP_j$ with orthogonal projectors that are unique, mutually annihilating, and sum to $I$.
• $f(A) = \sum f(\mu_j)P_j$, and each $P_j$ is a polynomial in $A$.
• Keeping the $k$ terms with the largest $\lvert \lambda \rvert$ gives the best rank-$k$ approximation, with error $\lvert \lambda_{k+1} \rvert$ in the spectral norm.

🧪 Practice

1. Find the spectral decomposition of $A = \begin{bmatrix} 7 & -2 \ -2 & 4 \end{bmatrix}$, and use it to write $A^{-1}$ and $A^3$ as combinations of the projectors.
2. Use the Lagrange formula to find the eigenspace projectors of $\begin{bmatrix} 1 & 2 \ 2 & 1 \end{bmatrix}$ without computing eigenvectors. Verify $P_1 + P_2 = I$ and $P_1P_2 = 0$.
3. For a symmetric matrix with a repeated eigenvalue, show that the rank-one terms $\mathbf{q}_i\mathbf{q}_i^T$ depend on the choice of eigenvectors, but their sum over the eigenspace does not.
4. Interview: You must multiply a dense $10^4 \times 10^4$ symmetric similarity matrix by thousands of vectors, and its eigenvalues decay quickly. How would you speed this up, how would you choose the rank, and how would you compute the factors? Hint: replace $A$ by $A_k$, so that each product costs $O(nk)$; the spectral-norm error is the first discarded $\lvert \lambda \rvert$; compute only the top-$k$ eigenpairs with a Krylov method such as Lanczos (Section 12.3) instead of a full eigh.

<a id="102-quadratic-forms"></a>

10.2 Quadratic Forms

This subchapter studies functions of the form $\mathbf{x}^TA\mathbf{x}$. It shows how a rotation of coordinates removes their cross terms, and how the signs of the eigenvalues classify the curves and surfaces they define.

Quadratic Forms and Matrix Representation

Why it matters. Many central quantities in science and engineering are homogeneous quadratic polynomials in several variables:

Setting	Quadratic form	Matrix
Squared length	$\mathbf{x}^T\mathbf{x} = x_1^2 + \cdots + x_n^2$	$I$
Weighted squared length	$\sum w_ix_i^2$	$\operatorname{diag}(w_1, \dots, w_n)$
Elastic energy of a spring network	$\tfrac12\mathbf{u}^TK\mathbf{u}$	stiffness matrix $K$
Variance of a portfolio	$\mathbf{w}^T\Sigma\mathbf{w}$	covariance matrix $\Sigma$
Squared Mahalanobis distance	$(\mathbf{x} - \boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu})$	$\Sigma^{-1}$
Second-order Taylor term	$\tfrac12\mathbf{h}^T\nabla^2 f(\mathbf{x}_0),\mathbf{h}$	Hessian matrix

Writing each one as $\mathbf{x}^TA\mathbf{x}$ makes the spectral theorem available. The eigenvalues and eigenvectors of $A$ then reveal the shape of the function: where it grows fastest, where it is flat, and whether it ever becomes negative.

Definition. A quadratic form on $\mathbb{R}^n$ is a function

$$Q(\mathbf{x}) = \mathbf{x}^TA\mathbf{x} = \sum_{i=1}^n\sum_{j=1}^n a_{ij}x_ix_j,$$

where $A$ is an $n \times n$ real matrix. Every term has total degree exactly 2. There are no linear terms and no constant, so $Q(c\mathbf{x}) = c^2Q(\mathbf{x})$.

Use the symmetric matrix. Many matrices give the same form, because $a_{ij}x_ix_j$ and $a_{ji}x_jx_i$ multiply the same monomial. Only $a_{ij} + a_{ji}$ matters. Equivalently, the skew-symmetric part of any matrix contributes nothing, since $\mathbf{x}^TK\mathbf{x} = 0$ for skew-symmetric $K$ (Section 4.2). So we may replace $A$ by its symmetric part $\tfrac12(A + A^T)$, and there is exactly one symmetric matrix for each quadratic form. By convention, "the matrix of $Q$" means this symmetric one. The form also determines the matrix:

$$a_{ii} = Q(\mathbf{e}i), \qquad a{ij} = \tfrac12\big(Q(\mathbf{e}_i + \mathbf{e}_j) - Q(\mathbf{e}_i) - Q(\mathbf{e}_j)\big).$$

Reading off the matrix. The coefficient of $x_i^2$ goes on the diagonal. The coefficient of the cross term $x_ix_j$ is split in half between positions $(i, j)$ and $(j, i)$:

$$Q(x, y, z) = 3x^2 + 4xy - y^2 + 6yz \qquad\Longrightarrow\qquad A = \begin{bmatrix} 3 & 2 & 0 \ 2 & -1 & 3 \ 0 & 3 & 0 \end{bmatrix}.$$

In two variables, writing the cross coefficient as $2b$ avoids the halves:

$$ax^2 + 2bxy + cy^2 = \begin{bmatrix} x & y \end{bmatrix}\begin{bmatrix} a & b \ b & c \end{bmatrix}\begin{bmatrix} x \ y \end{bmatrix}.$$

Calculus of quadratic forms. For symmetric $A$,

$$\nabla(\mathbf{x}^TA\mathbf{x}) = 2A\mathbf{x}, \qquad \nabla^2(\mathbf{x}^TA\mathbf{x}) = 2A.$$

For a non-symmetric $B$, the gradient is $(B + B^T)\mathbf{x}$, another reason to symmetrize first. Every smooth function looks like a quadratic form near a point, after the constant and linear terms:

$$f(\mathbf{x}_0 + \mathbf{h}) \approx f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0) \cdot \mathbf{h} + \tfrac12\mathbf{h}^T\nabla^2f(\mathbf{x}_0),\mathbf{h}.$$

At a critical point the linear term vanishes, and the quadratic form decides whether $f$ has a minimum, a maximum, or a saddle there (Section 10.3).

Bilinear forms and polarization. A symmetric matrix also defines the symmetric bilinear form $B(\mathbf{x}, \mathbf{y}) = \mathbf{x}^TA\mathbf{y}$, with $Q(\mathbf{x}) = B(\mathbf{x}, \mathbf{x})$. The quadratic form determines the bilinear one through the polarization identity

$$B(\mathbf{x}, \mathbf{y}) = \tfrac14\big(Q(\mathbf{x} + \mathbf{y}) - Q(\mathbf{x} - \mathbf{y})\big).$$

An inner product (Section 9.1) is exactly a symmetric bilinear form whose quadratic form is positive for every nonzero vector. Abstract bilinear forms are treated in Section 13.1. Over $\mathbb{C}$, the analogue is the Hermitian form $\mathbf{z}^*A\mathbf{z}$ with $A$ Hermitian, which is always real (Section 4.2).

import numpy as np

def form_matrix(n, coeffs):
    """Symmetric matrix of the form sum c_ij x_i x_j; coeffs maps (i, j) with i <= j to c_ij."""
    A = np.zeros((n, n))
    for (i, j), c in coeffs.items():
        if i == j:
            A[i, i] += c                       # x_i^2 coefficient goes on the diagonal
        else:
            A[i, j] += c / 2                   # cross-term coefficient is split between (i, j) and (j, i)
            A[j, i] += c / 2
    return A

# Q(x, y, z) = 3x^2 + 4xy - y^2 + 6yz
A = form_matrix(3, {(0, 0): 3, (0, 1): 4, (1, 1): -1, (1, 2): 6})
print(A)
# [[ 3.  2.  0.]
#  [ 2. -1.  3.]
#  [ 0.  3.  0.]]
x = np.array([1., 2., -1.])
print(x @ A @ x, 3*1 + 4*1*2 - 2**2 + 6*2*(-1))           # -5.0 -5

# A non-symmetric matrix defines the same form as its symmetric part
B = np.array([[3., 4., 0.], [0., -1., 6.], [0., 0., 0.]])   # all cross terms above the diagonal
print(x @ B @ x, np.allclose((B + B.T) / 2, A))            # -5.0 True

# The gradient of x^T A x is 2Ax: compare with central finite differences
h = 1e-6
num_grad = np.array([((x + h*e) @ A @ (x + h*e) - (x - h*e) @ A @ (x - h*e)) / (2*h) for e in np.eye(3)])
print(np.round(num_grad, 6), 2 * A @ x)                     # [14. -6. 12.] [14. -6. 12.]

# Polarization: the bilinear form is recovered from Q alone
Q = lambda v: v @ A @ v
y = np.array([0., 1., 3.])
print(np.isclose((Q(x + y) - Q(x - y)) / 4, x @ A @ y))     # True

Key takeaways

• A quadratic form is $Q(\mathbf{x}) = \mathbf{x}^TA\mathbf{x}$, a homogeneous polynomial of degree 2; $Q(c\mathbf{x}) = c^2Q(\mathbf{x})$.
• Each quadratic form has exactly one symmetric matrix: squared-term coefficients on the diagonal, half of each cross-term coefficient in the two off-diagonal positions.
• For symmetric $A$, the gradient is $2A\mathbf{x}$ and the Hessian is $2A$; quadratic forms are the second-order part of every smooth function.
• The polarization identity recovers the bilinear form $\mathbf{x}^TA\mathbf{y}$ from the quadratic form.

🧪 Practice

1. Write the symmetric matrix of $Q(x_1, x_2, x_3) = x_1^2 - 2x_2^2 + 5x_3^2 + 6x_1x_2 - 4x_2x_3 + x_1x_3$, and evaluate $Q(1, -1, 2)$ both directly and as $\mathbf{x}^TA\mathbf{x}$.
2. Expand $\mathbf{x}^TA\mathbf{x}$ for $A = \begin{bmatrix} 1 & 3 \ -1 & 2 \end{bmatrix}$. Find the symmetric matrix that gives the same form, and explain why $\mathbf{x}^TA\mathbf{x}$ cannot distinguish $A$ from $\tfrac12(A + A^T)$.
3. Compute the gradient and the Hessian of $f(\mathbf{x}) = (\mathbf{x} - \boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu})$ for symmetric invertible $\Sigma$.
4. Interview: A risk system computes the variance $\mathbf{w}^T\Sigma\mathbf{w}$ of a portfolio of 5000 assets and its gradient with respect to the weights. How expensive is each, and why must $\Sigma$ satisfy $\mathbf{w}^T\Sigma\mathbf{w} \ge 0$ for every $\mathbf{w}$? Hint: one matrix-vector product $\Sigma\mathbf{w}$ gives both the variance and the gradient $2\Sigma\mathbf{w}$; a variance cannot be negative, and an estimated $\Sigma$ that violates this lets an optimizer build "negative-risk" portfolios.

Change of Variable and Principal Axes Theorem

Why it matters. The cross terms are what make a quadratic form hard to read. Without them, $3u^2 + 7v^2$ is visibly an ellipse aligned with the axes, positive everywhere except at the origin. With them, $5x^2 - 4xy + 5y^2$ hides the same information. The spectral theorem provides a rotation of coordinates that removes every cross term at once. The new axes, the principal axes, are the natural directions of the problem: principal stresses in a material, principal moments of inertia of a rigid body, and principal components of a data set.

Change of variable. Substitute $\mathbf{x} = P\mathbf{y}$ with $P$ invertible:

$$Q(\mathbf{x}) = (P\mathbf{y})^TA(P\mathbf{y}) = \mathbf{y}^T(P^TAP),\mathbf{y}.$$

In the new variables, the form has matrix $P^TAP$, which is again symmetric. The goal is a $P$ that makes $P^TAP$ diagonal. If $P$ is orthogonal, then $P^TAP = P^{-1}AP$ is also a similarity, so the eigenvalues are preserved, and the spectral theorem says exactly how to choose $P$.

Theorem (Principal axes theorem). Let $A$ be an $n \times n$ symmetric matrix with orthogonal diagonalization $A = QDQ^T$. The orthogonal change of variable $\mathbf{x} = Q\mathbf{y}$ transforms $\mathbf{x}^TA\mathbf{x}$ into a form with no cross terms:

$$\mathbf{x}^TA\mathbf{x} = \mathbf{y}^TD\mathbf{y} = \lambda_1y_1^2 + \lambda_2y_2^2 + \cdots + \lambda_ny_n^2.$$

The columns of $Q$, the unit eigenvectors of $A$, are the principal axes of the form. The new coordinate $y_i = \mathbf{q}_i \cdot \mathbf{x}$ measures $\mathbf{x}$ along the $i$-th axis.

Why insist on an orthogonal $P$? An orthogonal change of variable preserves lengths and angles (Section 9.4), so it moves the curve $Q(\mathbf{x}) = c$ rigidly without distorting it. The move is a rotation, or a reflection, which can be avoided by flipping the sign of one column of $Q$ to make $\det Q = +1$. A non-orthogonal $P$ can also remove the cross terms (completing the square does this, see Congruence and Sylvester's Law of Inertia), but it stretches the picture. That is fine for deciding signs, but wrong for measuring axis lengths.

Geometry of $\mathbf{x}^TA\mathbf{x} = 1$. In principal coordinates the level set is $\sum\lambda_iy_i^2 = 1$. If all $\lambda_i > 0$, this is an ellipse or an ellipsoid whose $i$-th semi-axis points along $\mathbf{q}_i$ and has length $1/\sqrt{\lambda_i}$. Large eigenvalues give short axes: the form grows fastest in that direction, so it reaches the value 1 soonest.

Worked example. Identify the curve $5x^2 - 4xy + 5y^2 = 48$.

1. Matrix: $A = \begin{bmatrix} 5 & -2 \ -2 & 5 \end{bmatrix}$.
2. Eigenvalues: $\operatorname{tr}A = 10$, $\det A = 21$, so $\lambda = 3, 7$.
3. Unit eigenvectors: $A(1, 1) = (3, 3)$ gives $\mathbf{q}_1 = \frac{1}{\sqrt2}(1, 1)$ for $\lambda = 3$. $A(-1, 1) = (-7, 7)$ gives $\mathbf{q}_2 = \frac{1}{\sqrt2}(-1, 1)$ for $\lambda = 7$. $Q = [\mathbf{q}_1\ \mathbf{q}_2]$ has determinant $+1$: it is the rotation by $45°$.
4. New form: $3u^2 + 7v^2 = 48$, that is,

$$\frac{u^2}{16} + \frac{v^2}{48/7} = 1.$$

The curve is an ellipse with semi-axis $4$ along $(1, 1)$ and semi-axis $\sqrt{48/7} \approx 2.619$ along $(-1, 1)$.

  original coordinates (x, y)                           principal-axis coordinates (u, v)

  5x^2 - 4xy + 5y^2 = 48            [x]   [u]           3u^2 + 7v^2 = 48
  tilted ellipse, cross term        [y] = Q [v]         u^2/16 + v^2/(48/7) = 1
                                                        axis-aligned ellipse
                     1   [ 1  -1 ]
          with  Q = ---  [ 1   1 ]   (rotation by +45 degrees)
                    sqrt2

  u-axis = q1 = (1, 1)/sqrt2  : lambda = 3, semi-axis sqrt(48/3) = 4      (long axis: small lambda)
  v-axis = q2 = (-1, 1)/sqrt2 : lambda = 7, semi-axis sqrt(48/7) ~ 2.62   (short axis: large lambda)

In two dimensions, the angle directly. For $ax^2 + bxy + cy^2$, rotating the axes by an angle $\theta$ with $\cot 2\theta = (a - c)/b$ removes the cross term. For $a = c$, as in the example, $\theta = 45°$. This is the classical formula from analytic geometry. The eigenvector method gives the same result and works in any dimension.

Linear terms: shift, then rotate. For a full quadratic $\mathbf{x}^TA\mathbf{x} + \mathbf{b}^T\mathbf{x} + c = 0$ with $A$ invertible, substitute $\mathbf{x} = \mathbf{x}_0 + \mathbf{z}$ with center $\mathbf{x}_0 = -\tfrac12A^{-1}\mathbf{b}$. The linear terms cancel and

$$\mathbf{z}^TA\mathbf{z} = \mathbf{x}_0^TA\mathbf{x}_0 - c.$$

This is completing the square in matrix form. For example, $5x^2 - 4xy + 5y^2 - 18x + 24y - 15 = 0$ has center $\mathbf{x}_0 = (1, -2)$ and $\mathbf{x}_0^TA\mathbf{x}_0 - c = 33 + 15 = 48$. It is the ellipse of the example, translated to $(1, -2)$.

import numpy as np

# 5x^2 - 4xy + 5y^2 = 48
A = np.array([[5., -2.], [-2., 5.]])
w, Q = np.linalg.eigh(A)
Q[:, 0] *= np.sign(Q[0, 0])                  # eigenvector signs are arbitrary: pick one
if np.linalg.det(Q) < 0:
    Q[:, 1] *= -1                             # make Q a rotation (det +1), not a reflection
print(w)                                      # [3. 7.]: coefficients of u^2 and v^2
print(np.round(Q * np.sqrt(2), 4))            # columns: principal axis directions (times sqrt 2)
# [[ 1. -1.]
#  [ 1.  1.]]
print(np.round(Q.T @ A @ Q, 10) + 0.0)        # diagonal: no cross term in (u, v)
# [[3. 0.]
#  [0. 7.]]
print(np.round(np.sqrt(48 / w), 4))           # [4.     2.6186]: semi-axis lengths
print(round(np.degrees(np.arctan2(Q[1, 0], Q[0, 0])), 4))   # 45.0: rotation angle of the first axis

# Points built in (u, v) coordinates land on the original curve after x = Q y
t = np.linspace(0, 2 * np.pi, 7)
uv = np.vstack([4 * np.cos(t), np.sqrt(48 / 7) * np.sin(t)])
xy = Q @ uv
print(np.allclose(np.einsum('ik,ij,jk->k', xy, A, xy), 48))   # True: x^T A x = 48 for each column

# With linear terms: x^T A x + b^T x + c = 0 is centered at x0 = -A^{-1} b / 2
b, c = np.array([-18., 24.]), -15.
x0 = -np.linalg.solve(A, b) / 2
print(x0, x0 @ A @ x0 - c)                    # [ 1. -2.] 48.0: center, and the right-hand side of Q(x - x0) = 48

Key takeaways

• Under $\mathbf{x} = P\mathbf{y}$, the matrix of a quadratic form changes to $P^TAP$.
• Principal axes theorem: with $P = Q$ from $A = QDQ^T$, the form becomes $\sum\lambda_iy_i^2$, with no cross terms.
• The principal axes are the eigenvectors; on the level set $\mathbf{x}^TA\mathbf{x} = 1$ the semi-axes have lengths $1/\sqrt{\lambda_i}$, so large eigenvalues mean short axes.
• Orthogonal substitutions preserve shape; choose $\det Q = +1$ to get a pure rotation.
• Linear terms are removed by shifting to the center $\mathbf{x}_0 = -\tfrac12A^{-1}\mathbf{b}$ before rotating.

🧪 Practice

1. Find the principal axes and the standard form of $3x^2 + 2xy + 3y^2 = 8$, and give the lengths and directions of the semi-axes.
2. Show that $3x^2 + 2xy + 3y^2 - 8x - 8y = 0$ is the curve of Exercise 1 translated, by finding its center with $\mathbf{x}_0 = -\tfrac12A^{-1}\mathbf{b}$.
3. Using the orthogonal diagonalization of $\begin{bmatrix} 2 & 1 & 1 \ 1 & 2 & 1 \ 1 & 1 & 2 \end{bmatrix}$ from earlier in this chapter, find an orthogonal substitution that removes all cross terms from $2x_1^2 + 2x_2^2 + 2x_3^2 + 2x_1x_2 + 2x_1x_3 + 2x_2x_3$, and describe the surface where this form equals 4.
4. Interview: You fit a 2D Gaussian with mean $\boldsymbol{\mu}$ and covariance $\Sigma$ to data, and must draw its 95% confidence ellipse. How do you get the ellipse's center, orientation, and axis lengths? Hint: the ellipse is $(\mathbf{x} - \boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu}) = \chi^2_2(0.95) \approx 5.991$; $\Sigma^{-1}$ has the same eigenvectors as $\Sigma$ and reciprocal eigenvalues, so the semi-axes are $\sqrt{5.991,\lambda_i(\Sigma)}$ along the eigenvectors of $\Sigma$.

Classification of Conics and Quadric Surfaces

Why it matters. Every second-degree equation in two variables describes an ellipse, a hyperbola, a parabola, or a degenerate case such as a pair of lines. In three variables there are about a dozen surface types. Knowing which one an equation describes matters in orbit mechanics (conic orbits), optics (parabolic mirrors), computer graphics and CAD (quadric primitives), and in recognizing the level sets of an objective function. The principal axes theorem reduces the whole classification to the signs of the eigenvalues.

Conics. A general conic is

$$ax^2 + bxy + cy^2 + dx + ey + f = 0, \qquad A_2 = \begin{bmatrix} a & b/2 \ b/2 & c \end{bmatrix},$$

where $A_2$ is the matrix of the quadratic part. Rotate to principal axes to get $\lambda_1u^2 + \lambda_2v^2 + d'u + e'v + f = 0$, then complete the squares. The outcome depends on the signs of $\lambda_1$ and $\lambda_2$, or equivalently on $\det A_2 = \lambda_1\lambda_2 = -\tfrac14(b^2 - 4ac)$:

Eigenvalues of $A_2$	$\det A_2$	$b^2 - 4ac$	Non-degenerate conic	Standard form
same sign	$> 0$	$< 0$	ellipse (circle if $\lambda_1 = \lambda_2$)	$\dfrac{u^2}{p^2} + \dfrac{v^2}{q^2} = 1$
opposite signs	$< 0$	$> 0$	hyperbola	$\dfrac{u^2}{p^2} - \dfrac{v^2}{q^2} = 1$
one is zero	$= 0$	$= 0$	parabola	$v^2 = 4pu$

   lambda1, lambda2 same sign       lambda1, lambda2 opposite signs     one lambda zero
          ellipse                          hyperbola                       parabola

           .-----.                    \               /                  \           /
         /         \                   \             /                    \         /
        |     +     |                   )     +     (                      \       /
         \         /                   /             \                      '.___.'
           '-----'                    /               \
   bounded, closed curve             two branches, two asymptotes       one branch, no asymptotes

The discriminant $b^2 - 4ac$ does not change under rotations, because it equals $-4\det A_2$ and similar matrices share determinants. So it classifies the conic without rotating anything.

Degenerate conics. The full $3 \times 3$ symmetric matrix

$$M = \begin{bmatrix} a & b/2 & d/2 \ b/2 & c & e/2 \ d/2 & e/2 & f \end{bmatrix}, \qquad \text{conic: } \begin{bmatrix} x & y & 1 \end{bmatrix} M \begin{bmatrix} x \ y \ 1 \end{bmatrix} = 0,$$

detects degeneracy: the conic is degenerate exactly when $\det M = 0$. Then it is a single point ($x^2 + y^2 = 0$), two crossing lines ($x^2 - y^2 = 0$), two parallel lines ($x^2 = 1$), one line ($x^2 = 0$), or empty ($x^2 = -1$). When $\det A_2 > 0$ and $\det M \ne 0$, the ellipse is real when $\det M / (a + c) < 0$ and empty otherwise. An example of an empty one is $x^2 + y^2 + 1 = 0$, called an imaginary ellipse.

Quadric surfaces. In $\mathbb{R}^3$, after rotating to principal axes and shifting, an equation $\mathbf{x}^TA\mathbf{x} + \mathbf{b}^T\mathbf{x} + c = 0$ becomes one of the following standard forms ($p, q, r > 0$):

Signs of $(\lambda_1, \lambda_2, \lambda_3)$	Standard form	Surface
$+, +, +$	$\dfrac{u^2}{p^2} + \dfrac{v^2}{q^2} + \dfrac{w^2}{r^2} = 1$	ellipsoid
$+, +, -$	$\dfrac{u^2}{p^2} + \dfrac{v^2}{q^2} - \dfrac{w^2}{r^2} = 1$	hyperboloid of one sheet
$+, -, -$	$\dfrac{u^2}{p^2} - \dfrac{v^2}{q^2} - \dfrac{w^2}{r^2} = 1$	hyperboloid of two sheets
$+, +, -$	$\dfrac{u^2}{p^2} + \dfrac{v^2}{q^2} - \dfrac{w^2}{r^2} = 0$	elliptic cone
$+, +, 0$	$w = \dfrac{u^2}{p^2} + \dfrac{v^2}{q^2}$	elliptic paraboloid (a bowl)
$+, -, 0$	$w = \dfrac{u^2}{p^2} - \dfrac{v^2}{q^2}$	hyperbolic paraboloid (a saddle)
$+, +, 0$	$\dfrac{u^2}{p^2} + \dfrac{v^2}{q^2} = 1$	elliptic cylinder
$+, -, 0$	$\dfrac{u^2}{p^2} - \dfrac{v^2}{q^2} = 1$	hyperbolic cylinder
$+, 0, 0$	$u^2 = 4pw$	parabolic cylinder

Multiplying an equation by $-1$ flips all signs, including the sign of the right-hand side, so the table lists each pattern only once. A zero eigenvalue produces a paraboloid if a linear term survives in that direction, and a cylinder if it does not. One sheet or two depends on how many signs agree with the right-hand side. For example, $u^2 + v^2 - w^2 = 1$ meets every horizontal plane in a circle (one sheet), while $u^2 - v^2 - w^2 = 1$ forces $u^2 \ge 1$ (two sheets).

Worked examples.

• $x^2 + 4xy + y^2 = 3$: $A_2 = \begin{bmatrix} 1 & 2 \ 2 & 1 \end{bmatrix}$ has eigenvalues $3$ and $-1$, with opposite signs. In principal coordinates, $3u^2 - v^2 = 3$, or $u^2 - v^2/3 = 1$: a hyperbola whose $u$-axis points along $(1, 1)/\sqrt2$, with asymptotes $v = \pm\sqrt3,u$.
• $2xy + 2yz + 2zx = 1$: $A = J - I$ (with $J$ the all-ones matrix) has eigenvalues $2, -1, -1$. In principal coordinates, $2u^2 - v^2 - w^2 = 1$: a hyperboloid of two sheets around the axis $(1, 1, 1)/\sqrt3$.

import numpy as np

def classify_conic(a, b, c, d, e, f, tol=1e-12):
    """Classify a x^2 + b xy + c y^2 + d x + e y + f = 0 (absolute tol: fine for modest integer inputs)."""
    A2 = np.array([[a, b/2], [b/2, c]])                       # quadratic part
    M = np.array([[a, b/2, d/2], [b/2, c, e/2], [d/2, e/2, f]])   # full 3 x 3 matrix
    det2, det3 = np.linalg.det(A2), np.linalg.det(M)
    if abs(det3) < tol:
        kind = {1: "point", -1: "two crossing lines", 0: "parallel lines, one line, or empty"}
        return "degenerate: " + kind[int(np.sign(np.round(det2, 12)))]
    if det2 > tol:                                            # eigenvalues of A2 have the same sign
        return "ellipse" if det3 / (a + c) < 0 else "empty (imaginary ellipse)"
    if det2 < -tol:                                           # opposite signs
        return "hyperbola"
    return "parabola"                                         # one zero eigenvalue

tests = {
    "5x^2 - 4xy + 5y^2 - 48":      (5, -4, 5, 0, 0, -48),
    "x^2 + 4xy + y^2 - 3":         (1, 4, 1, 0, 0, -3),
    "x^2 - 2xy + y^2 - 4x":        (1, -2, 1, -4, 0, 0),
    "x^2 + y^2 + 1":               (1, 0, 1, 0, 0, 1),
    "x^2 - y^2":                   (1, 0, -1, 0, 0, 0),
    "x^2 + 2xy + y^2 - 1":         (1, 2, 1, 0, 0, -1),
}
for eq, coeffs in tests.items():
    print(f"{eq:28s} -> {classify_conic(*coeffs)}")
# 5x^2 - 4xy + 5y^2 - 48       -> ellipse
# x^2 + 4xy + y^2 - 3          -> hyperbola
# x^2 - 2xy + y^2 - 4x         -> parabola
# x^2 + y^2 + 1                -> empty (imaginary ellipse)
# x^2 - y^2                    -> degenerate: two crossing lines
# x^2 + 2xy + y^2 - 1          -> degenerate: parallel lines, one line, or empty   ((x + y)^2 = 1: parallel lines)

# Quadric surface 2xy + 2yz + 2zx = 1: the eigenvalue signs decide the shape
A3 = np.array([[0., 1., 1.], [1., 0., 1.], [1., 1., 0.]])
print(np.round(np.linalg.eigvalsh(A3), 10))                  # [-1. -1.  2.]: one positive, two negative -> two sheets

Key takeaways

• Rotating to principal axes and completing squares reduces every conic and quadric to a standard form.
• For conics: eigenvalues of the quadratic part with the same sign give an ellipse, opposite signs a hyperbola, and one zero eigenvalue a parabola; the discriminant $b^2 - 4ac = -4\det A_2$ decides this without rotating.
• The determinant of the $3 \times 3$ matrix $M$ detects degenerate conics (points, line pairs, empty sets).
• For quadric surfaces, the sign pattern of the three eigenvalues, the right-hand side, and any surviving linear term pick out ellipsoids, hyperboloids, cones, paraboloids, and cylinders.

🧪 Practice

1. Classify each conic: (a) $2x^2 + 3xy - 2y^2 = 5$; (b) $x^2 - 6xy + 9y^2 + 4x = 0$; (c) $4x^2 + 4xy + 5y^2 = 0$.
2. Identify the quadric surface $x^2 + y^2 - z^2 + 2xy = 1$ by finding the eigenvalues of its matrix, and give its principal axes.
3. Show that the discriminant $b^2 - 4ac$ is unchanged when the coordinate axes are rotated, using the relation $b^2 - 4ac = -4\det A_2$.
4. Interview: A robotics pipeline fits a general quadric $\mathbf{x}^TA\mathbf{x} + \mathbf{b}^T\mathbf{x} + c = 0$ to a 3D point cloud and must decide programmatically whether the fit is an ellipsoid, and if so report its center, axes, and radii. Outline the computation and its numerical pitfalls. Hint: center $\mathbf{x}_0 = -\tfrac12A^{-1}\mathbf{b}$, right-hand side $k = \mathbf{x}_0^TA\mathbf{x}_0 - c$; ellipsoid if all eigenvalues of $A$ share the sign of $k$; radii $\sqrt{k/\lambda_i}$ along the eigenvectors; compare eigenvalues against a tolerance relative to $\max\lvert\lambda_i\rvert$, not against zero.

Congruence and Sylvester's Law of Inertia

Why it matters. Many questions about a quadratic form depend only on the signs of its eigenvalues. Is this critical point a minimum? Is this conic an ellipse or a hyperbola? In how many independent directions does this energy decrease? Computing eigenvalues exactly is hard, since by hand it means solving a polynomial equation. Sylvester's law of inertia says that any invertible change of variables that removes the cross terms reveals the same sign counts. So completing the square, or Gaussian elimination, which are finite and exact, is enough.

Congruence. Two symmetric matrices $A$ and $B$ are congruent if $B = P^TAP$ for some invertible $P$. This is how the matrix of a quadratic form changes under $\mathbf{x} = P\mathbf{y}$. Compare it with similarity, which is how the matrix of a linear map changes under a change of basis:

	Similarity $P^{-1}AP$	Congruence $P^TAP$
Arises from	changing basis for a linear map	changing variables in a quadratic form
Preserves	eigenvalues, trace, determinant, rank	symmetry, rank, inertia (eigenvalue signs)
Determinant	unchanged	multiplied by $(\det P)^2 > 0$, so its sign is unchanged
Normal form for symmetric $A$	$\operatorname{diag}(\lambda_1, \dots, \lambda_n)$	$\operatorname{diag}(I_{n_+}, -I_{n_-}, 0_{n_0})$

The two notions coincide when $P$ is orthogonal, which is the setting of the principal axes theorem.

Inertia. The inertia of a symmetric matrix $A$ is the triple

$$\operatorname{In}(A) = (n_+, n_-, n_0),$$

the numbers of positive, negative, and zero eigenvalues. Then $\operatorname{rank}A = n_+ + n_-$, and the signature is $n_+ - n_-$ (some texts call the pair $(n_+, n_-)$ the signature).

Theorem (Sylvester's law of inertia). Congruent symmetric matrices have the same inertia. Equivalently, if a quadratic form is written as $\sum d_iy_i^2$ by any invertible change of variables, then the numbers of positive, negative, and zero coefficients $d_i$ do not depend on the change of variables. In particular, every symmetric matrix is congruent to $\operatorname{diag}(I_{n_+}, -I_{n_-}, 0)$: take the principal axes and rescale each $y_i$ by $1/\sqrt{\lvert \lambda_i \rvert}$.

Proof sketch. The number $n_0 = n - \operatorname{rank}$ is preserved, because $P$ is invertible. Suppose the same form equals both $y_1^2 + \cdots + y_p^2 - (\text{squares})$ and $z_1^2 + \cdots + z_{p'}^2 - (\text{squares})$ in two coordinate systems, with $p < p'$. On the subspace $U = {y_1 = \cdots = y_p = 0}$, of dimension $n - p$, the form is $\le 0$. On the subspace $W$ where only $z_1, \dots, z_{p'}$ can be nonzero, of dimension $p'$, the form is $> 0$ except at $\mathbf{0}$. Since $(n - p) + p' > n$, the two subspaces share a nonzero vector. The form would be both $\le 0$ and $> 0$ there, a contradiction. $\blacksquare$

The proof also shows that $n_+$ is the largest dimension of a subspace on which the form is positive definite. This characterization does not mention any coordinates.

Method 1: completing the square (Lagrange's method).

$$x^2 + 4xy + y^2 = (x + 2y)^2 - 3y^2.$$

With $u = x + 2y$ and $v = y$, an invertible but non-orthogonal substitution, the form is $u^2 - 3v^2$, so the inertia is $(1, 1, 0)$. The eigenvalues, $3$ and $-1$, have the same signs but different sizes, as the theorem predicts.

Method 2: symmetric elimination. Gaussian elimination without row exchanges on a symmetric matrix gives $A = LDL^T$ (Section 4.4), with $L$ unit lower triangular and $D$ holding the pivots. This is a congruence, $D = L^{-1}AL^{-T}$, so the signs of the pivots are the signs of the eigenvalues. For

$$A = \begin{bmatrix} 1 & 2 & 0 \ 2 & 1 & 1 \ 0 & 1 & 3 \end{bmatrix},$$

the pivots are $1$, $1 - 4 = -3$, and $3 - (1)(1)/(-3) = 10/3$, so $\operatorname{In}(A) = (2, 1, 0)$. The eigenvalues are approximately $-1.125$, $2.363$, and $3.762$. The pivots are not the eigenvalues, but their signs agree, and so do their products: $1 \cdot (-3) \cdot \tfrac{10}{3} = -10 = \det A$.

Method 3: counting eigenvalues by slicing. The eigenvalues of $A - \sigma I$ are $\lambda_i - \sigma$. So $n_-(A - \sigma I)$ is the number of eigenvalues of $A$ below $\sigma$. One $LDL^T$ factorization counts them without computing any eigenvalue. Repeating this with different shifts $\sigma$ and bisecting locates individual eigenvalues. This is the basis of the bisection method for symmetric tridiagonal matrices, and of the checks that sparse eigensolvers use to make sure no eigenvalue was missed.

Where inertia appears.

• Critical points. At a critical point with an invertible Hessian, the number of negative eigenvalues of the Hessian is the Morse index: $0$ means a local minimum, $n$ a local maximum, and anything in between a saddle (Section 10.3).
• Conics and quadrics. The classification tables of the previous topic depend only on inertia.
• Special relativity. The Minkowski form $c^2t^2 - x^2 - y^2 - z^2$ has inertia $(1, 3, 0)$ in every inertial frame, because Lorentz transformations are congruences that preserve it.

import numpy as np
from scipy.linalg import ldl

def inertia(A, tol=1e-10):
    w = np.linalg.eigvalsh(A)
    return int((w > tol).sum()), int((w < -tol).sum()), int((np.abs(w) <= tol).sum())

A = np.array([[1., 2., 0.],
              [2., 1., 1.],
              [0., 1., 3.]])
print(np.round(np.linalg.eigvalsh(A), 4), inertia(A))   # [-1.1249  2.3633  3.7616] (2, 1, 0)

# Congruence by a random invertible P changes the eigenvalues but not their signs
rng = np.random.default_rng(0)
for _ in range(3):
    P = rng.standard_normal((3, 3))
    B = P.T @ A @ P
    print(np.round(np.linalg.eigvalsh(B), 2), inertia(B))
# [-0.13  1.73  8.57] (2, 1, 0)
# [-0.79  0.22 28.91] (2, 1, 0)
# [-2.09  1.06  9.99] (2, 1, 0)

# Symmetric elimination A = L D L^T is a congruence: the signs of the pivots give the inertia
def ldl_no_pivot(A):
    U = A.astype(float).copy()
    n = len(U)
    L = np.eye(n)
    for k in range(n - 1):
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]       # assumes nonzero pivots
            U[i] -= L[i, k] * U[k]
    return L, np.diag(U)

L, d = ldl_no_pivot(A)
print(np.round(d, 4), np.allclose(L @ np.diag(d) @ L.T, A))   # [ 1.     -3.      3.3333] True

# Counting eigenvalues below a shift sigma via the inertia of A - sigma I.
# scipy's ldl pivots for stability, so D may contain 2 x 2 blocks (eigvalsh(D) is used here for brevity;
# real code reads the signs off the 1 x 1 and 2 x 2 blocks directly in O(n))
def count_below(A, sigma):
    _, D, _ = ldl(A - sigma * np.eye(len(A)))
    return inertia(D)[1]

T = 2 * np.eye(50) - np.eye(50, k=1) - np.eye(50, k=-1)      # second-difference matrix
print(count_below(T, 1.0), int((np.linalg.eigvalsh(T) < 1.0).sum()))   # 16 16

Key takeaways

• $A$ and $P^TAP$ ($P$ invertible) are congruent: they represent the same quadratic form in different variables.
• Congruence changes eigenvalues but preserves the inertia $(n_+, n_-, n_0)$ (Sylvester's law of inertia).
• Inertia can be computed exactly without eigenvalues, by completing squares or from the signs of the pivots of $A = LDL^T$.
• $n_-(A - \sigma I)$ counts the eigenvalues of $A$ below $\sigma$, which allows eigenvalues to be located by bisection.
• Every symmetric matrix is congruent to $\operatorname{diag}(I_{n_+}, -I_{n_-}, 0)$.

🧪 Practice

1. Complete the square to find the inertia of $Q = x^2 + 2xy + 3y^2 + 4yz + z^2$. Confirm with the pivots of its matrix, and check that their product equals the determinant.
2. Find an invertible $P$ with $P^T\begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}P = \begin{bmatrix} 1 & 0 \ 0 & -1 \end{bmatrix}$. (The form is $2xy$; try $x = u + v$, $y = u - v$, then rescale.)
3. Prove that $n_+(A)$ equals the largest dimension of a subspace $W$ such that $\mathbf{x}^TA\mathbf{x} > 0$ for every nonzero $\mathbf{x} \in W$.
4. Interview: How would you count the eigenvalues of a large sparse symmetric matrix that lie in an interval $[a, b)$, without computing any of them? Hint: by Sylvester's law, the count is $n_-(A - bI) - n_-(A - aI)$, and each term comes from the pivots of a sparse $LDL^T$ factorization.

<a id="103-definiteness-and-optimization"></a>

10.3 Definiteness and Optimization

This subchapter singles out the quadratic forms that are positive in every direction, gives practical tests for recognizing them, and factors them with Cholesky. It then characterizes eigenvalues as the answers to quadratic optimization problems.

Positive Definite and Semidefinite Matrices

Why it matters. A positive definite matrix describes a quantity that grows in every direction away from the origin. Examples are an elastic energy that rises whenever the structure is deformed, a variance that is positive for every non-trivial combination of features, and a squared distance. These matrices define inner products (Section 9.1), make least squares and Newton's method well posed, certify that a critical point is a minimum, and are exactly the covariance matrices and kernel matrices of statistics and machine learning. Positive definite matrices play the role among matrices that positive numbers play among numbers.

Analogy. Picture the graph $z = \mathbf{x}^TA\mathbf{x}$ over the plane. A positive definite form is a bowl: every path away from the bottom goes uphill. A semidefinite form is a trough: uphill in some directions, flat along the bottom of the valley. An indefinite form is a saddle (a riding saddle, or a Pringles chip): uphill one way, downhill another.

Definitions. Let $A$ be a real symmetric $n \times n$ matrix and $Q(\mathbf{x}) = \mathbf{x}^TA\mathbf{x}$.

Class	Condition for all $\mathbf{x} \ne \mathbf{0}$	Eigenvalues	Graph of $Q$ ($n = 2$)	Notation
positive definite (PD)	$Q(\mathbf{x}) > 0$	all $> 0$	bowl	$A \succ 0$
positive semidefinite (PSD)	$Q(\mathbf{x}) \ge 0$	all $\ge 0$	trough	$A \succeq 0$
negative definite	$Q(\mathbf{x}) < 0$	all $< 0$	upside-down bowl	$A \prec 0$
negative semidefinite	$Q(\mathbf{x}) \le 0$	all $\le 0$	upside-down trough	$A \preceq 0$
indefinite	$Q$ takes both signs	both signs occur	saddle

The eigenvalue column follows from the principal axes theorem. With $\mathbf{x} = Q\mathbf{y}$, $\mathbf{x}^TA\mathbf{x} = \sum\lambda_iy_i^2$, which is positive for every $\mathbf{y} \ne \mathbf{0}$ exactly when every $\lambda_i > 0$. If some $\lambda_i \le 0$, then $\mathbf{x} = \mathbf{q}_i$ gives $Q(\mathbf{q}_i) = \lambda_i \le 0$. For Hermitian matrices the definitions are the same with $\mathbf{z}^*A\mathbf{z}$.

   positive definite                positive semidefinite              indefinite
   level sets: nested ellipses      level sets: parallel lines         level sets: hyperbolas
   (bowl, unique minimum at 0)      (trough: flat along null space)    (saddle: up in x, down in y)

        .-----------.                   |   |   | |   |   |              \    |    /
      /   .-------.   \                 |   |   | |   |   |            .  \   |   /  .
     |   |    0    |   |                |   |   | |   |   |              '. \  |  / .'
      \   '-------'   /                 |   |   | |   |   |            ---  ) 0 (  ---
        '-----------'                   |   |   | |   |   |              .' /  |  \ '.
                                        |   |   | |   |   |            '  /   |   \  '
   x^2 + 2y^2                        x^2 (zero along the y-axis)           x^2 - y^2

Properties. Let $A$ and $B$ be symmetric.

• Diagonal entries. If $A \succ 0$, then $a_{ii} = \mathbf{e}i^TA\mathbf{e}i > 0$. If $A \succeq 0$, then $a{ii} \ge 0$, and $a{ii} = 0$ forces the whole $i$-th row and column to vanish.
• Invertibility. $A \succ 0$ implies $\det A = \prod\lambda_i > 0$, so $A$ is invertible. Moreover, $A^{-1} \succ 0$, because its eigenvalues are $1/\lambda_i$.
• Cone. If $A, B \succeq 0$ and $s, t \ge 0$, then $sA + tB \succeq 0$, and it is $\succ 0$ if one of them is $\succ 0$ with a positive coefficient. The PSD matrices form a convex cone.
• Congruence. If $A \succ 0$ and $B$ is $n \times k$, then $B^TAB \succeq 0$, and $B^TAB \succ 0$ when $B$ has independent columns, because $\mathbf{x}^TB^TAB\mathbf{x} = (B\mathbf{x})^TA(B\mathbf{x})$.
• Gram matrices. With $A = I$: $X^TX \succeq 0$ for every matrix $X$, since $\mathbf{x}^TX^TX\mathbf{x} = \lVert X\mathbf{x} \rVert^2$. It is $\succ 0$ exactly when the columns of $X$ are independent. Conversely, every PSD matrix is a Gram matrix, $A = B^TB$, for example with $B = A^{1/2}$ or with a Cholesky factor (next topics).
• Covariance matrices. For a random vector $\mathbf{X}$ with covariance $\Sigma$, $\mathbf{w}^T\Sigma\mathbf{w} = \operatorname{Var}(\mathbf{w}^T\mathbf{X}) \ge 0$. So $\Sigma \succeq 0$, and $\Sigma \succ 0$ unless some non-trivial combination of the features is constant.
• Principal submatrices. Every principal submatrix of a PD matrix is PD: restrict $\mathbf{x}$ to vectors that are zero outside the chosen indices.
• Loewner order. $A \succeq B$ means $A - B \succeq 0$. This is a partial order. For example, $\operatorname{diag}(1, 0)$ and $\operatorname{diag}(0, 1)$ are not comparable.
• Products. $AB$ is usually not symmetric, even when $A, B \succ 0$. Its eigenvalues are still positive, because $AB$ is similar to $A^{1/2}BA^{1/2} \succ 0$.

Definiteness is about symmetric matrices. For a non-symmetric $A$, $\mathbf{x}^TA\mathbf{x}$ only sees the symmetric part $\tfrac12(A + A^T)$ (Section 10.2). The matrix $\begin{bmatrix} 1 & 4 \ -4 & 1 \end{bmatrix}$ has $\mathbf{x}^TA\mathbf{x} = x_1^2 + x_2^2 > 0$, yet its eigenvalues are $1 \pm 4i$. So this guide, like most texts, applies the words "positive definite" only to symmetric (or Hermitian) matrices.

import numpy as np

def definiteness(A, tol=1e-10):
    w = np.linalg.eigvalsh(A)                  # assumes A is symmetric
    if np.all(w > tol):   return "positive definite"
    if np.all(w >= -tol): return "positive semidefinite"
    if np.all(w < -tol):  return "negative definite"
    if np.all(w <= tol):  return "negative semidefinite"
    return "indefinite"

for M in ([[2, -1], [-1, 2]], [[1, 1], [1, 1]], [[1, 3], [3, 1]], [[-2, 0], [0, -5]]):
    M = np.array(M, float)
    print(M.tolist(), np.linalg.eigvalsh(M), definiteness(M))
# [[2.0, -1.0], [-1.0, 2.0]] [1. 3.] positive definite
# [[1.0, 1.0], [1.0, 1.0]] [0. 2.] positive semidefinite
# [[1.0, 3.0], [3.0, 1.0]] [-2.  4.] indefinite
# [[-2.0, 0.0], [0.0, -5.0]] [-5. -2.] negative definite

# Gram matrices X^T X are always PSD; PD exactly when the columns of X are independent
rng = np.random.default_rng(5)
X = rng.standard_normal((100, 3))
print(definiteness(X.T @ X))                  # positive definite
X_dep = np.column_stack([X, X[:, 0] - 2 * X[:, 1]])          # 4th column depends on the others
print(definiteness(X_dep.T @ X_dep), np.round(np.linalg.eigvalsh(X_dep.T @ X_dep)[0], 8) + 0.0)
# positive semidefinite 0.0                  the zero eigenvalue's eigenvector is (1, -2, 0, -1)/sqrt(6)

# Covariance matrix: w^T Sigma w is the variance of the combined feature X w, so it is never negative
Sigma = np.cov(X, rowvar=False)
wvec = np.array([1., -2., 0.5])
print(np.isclose(wvec @ Sigma @ wvec, np.var(X @ wvec, ddof=1)))   # True

# Definiteness is about symmetric matrices: x^T A x only "sees" the symmetric part of A
A = np.array([[1., 4.], [-4., 1.]])           # x^T A x = x1^2 + x2^2 > 0, yet eigenvalues are 1 +- 4i
print(np.linalg.eigvals(A), definiteness((A + A.T) / 2))      # [1.+4.j 1.-4.j] positive definite

Key takeaways

• A symmetric $A$ is positive definite when $\mathbf{x}^TA\mathbf{x} > 0$ for all $\mathbf{x} \ne \mathbf{0}$, which is equivalent to all eigenvalues being positive; semidefinite allows zero.
• PD matrices have positive diagonals, positive determinants, PD inverses, and PD principal submatrices; PSD matrices form a convex cone.
• $X^TX$ is always PSD and is PD exactly when $X$ has independent columns; every PSD matrix is such a Gram matrix.
• Covariance matrices are PSD because variances are nonnegative.
• Definiteness is defined for symmetric (Hermitian) matrices; $\mathbf{x}^TA\mathbf{x}$ ignores the skew-symmetric part of $A$.

🧪 Practice

1. Classify $\begin{bmatrix} 3 & 1 \ 1 & 3 \end{bmatrix}$, $\begin{bmatrix} 1 & -2 \ -2 & 4 \end{bmatrix}$, $\begin{bmatrix} 2 & 3 \ 3 & 2 \end{bmatrix}$, and $\begin{bmatrix} -1 & 1 \ 1 & -1 \end{bmatrix}$ as PD, PSD, ND, NSD, or indefinite.
2. Prove that if $A \succ 0$ then $A^{-1} \succ 0$, and that $A, B \succ 0$ implies $A + B \succ 0$. Give an example where $A, B \succ 0$ but $AB$ is not symmetric.
3. Show that $X^TX$ is PD exactly when the columns of $X$ are linearly independent. What does this say about when the normal equations of least squares (Section 9.5) have a unique solution?
4. Interview: A sample covariance matrix estimated from 50 daily returns of 200 assets is passed to a portfolio optimizer that inverts it. What goes wrong, and what are standard remedies? Hint: the rank is at most 49, so the matrix is PSD but singular; its null directions look like riskless portfolios; remedies include shrinkage $\Sigma + \delta I$ (Ledoit-Wolf), factor models, and more data.

Tests for Definiteness

Why it matters. Computing all eigenvalues is the most informative way to decide definiteness, but often not the cheapest, and impractical by hand. Several other tests are equivalent for symmetric matrices. Each fits a different situation: a quick check on a $2 \times 2$ Hessian, a proof, or a production numerical routine. Knowing why they agree also ties together elimination, determinants, and eigenvalues.

Theorem (Tests for positive definiteness). For a real symmetric $n \times n$ matrix $A$, the following are equivalent:

Test	Condition	Cost / use
Energy	$\mathbf{x}^TA\mathbf{x} > 0$ for all $\mathbf{x} \ne \mathbf{0}$	the definition; used in proofs
Eigenvalues	all $\lambda_i > 0$	$O(n^3)$ with eigvalsh; shows the margin $\lambda_{\min}$
Pivots	elimination without row exchanges produces $n$ positive pivots	$\approx n^3/3$ flops
Sylvester's criterion	all leading principal minors $D_k = \det A_{[1..k, 1..k]}$ are $> 0$	convenient by hand for $n \le 3$
Cholesky	$A = LL^T$ with $L$ lower triangular with positive diagonal	$\approx n^3/3$ flops; the standard numerical test
Gram	$A = B^TB$ for some $B$ with independent columns	used in proofs and modeling

Why the tests agree. Energy $\Leftrightarrow$ eigenvalues is the principal axes theorem. For eigenvalues $\Leftrightarrow$ pivots: elimination without row exchanges gives $A = LDL^T$, a congruence, so by Sylvester's law of inertia the pivots have the same signs as the eigenvalues. Pivots $\Leftrightarrow$ minors: elimination on the top-left $k \times k$ block never touches rows below it, so $D_k = d_1d_2\cdots d_k$ and

$$d_k = \frac{D_k}{D_{k-1}} \qquad (D_0 = 1).$$

So all pivots are positive exactly when all leading minors are. If a pivot is zero, elimination stops, but then $D_k = 0$ for that $k$ and $A$ is not PD either. Pivots $\Leftrightarrow$ Cholesky: $LDL^T = (LD^{1/2})(LD^{1/2})^T$ needs $d_i > 0$ to take square roots (next topic). Gram $\Leftrightarrow$ energy: $\mathbf{x}^TB^TB\mathbf{x} = \lVert B\mathbf{x} \rVert^2$.

The $2 \times 2$ case. $\begin{bmatrix} a & b \ b & c \end{bmatrix} \succ 0$ if and only if $a > 0$ and $ac - b^2 > 0$. The pivots are $a$ and $c - b^2/a$.

Worked example. Take

$$A = \begin{bmatrix} 2 & -1 & 0 \ -1 & 2 & -1 \ 0 & -1 & 2 \end{bmatrix}.$$

• Leading minors: $D_1 = 2$, $D_2 = 4 - 1 = 3$, $D_3 = 2 \cdot 3 - (-1)(-2) = 4$. All positive.
• Pivots: $2$, $3/2$, $4/3$, the ratios $D_k/D_{k-1}$.
• Eigenvalues: $2 - \sqrt2$, $2$, $2 + \sqrt2$. All positive.
• Energy: $\mathbf{x}^TA\mathbf{x} = x_1^2 + (x_1 - x_2)^2 + (x_2 - x_3)^2 + x_3^2$. This sum of squares is zero only when $x_1 = x_2 = x_3 = 0$.

The last form shows where this matrix comes from: it is the energy of a chain of springs with both ends fixed.

Semidefinite and negative definite variants.

• Trap: leading minors $\ge 0$ do not imply PSD. The matrix $\begin{bmatrix} 0 & 0 \ 0 & -1 \end{bmatrix}$ has leading minors $0, 0$ but is negative semidefinite. A correct PSD criterion requires all principal minors (every choice of the same rows and columns, $2^n - 1$ of them) to be $\ge 0$. In practice, check $\lambda_{\min} \ge -\text{tol}$ instead.
• $A \prec 0$ if and only if $-A \succ 0$, which in terms of minors means that the leading minors alternate in sign, starting negative: $(-1)^kD_k > 0$.

Quick necessary and sufficient conditions.

• Necessary: a PD matrix has $a_{ii} > 0$ and $a_{ij}^2 < a_{ii}a_{jj}$ (each $2 \times 2$ principal minor is positive), so its largest entry in absolute value lies on the diagonal. Any violation rules out PD immediately.
• Sufficient: a symmetric matrix with positive diagonal that is strictly diagonally dominant, $a_{ii} > \sum_{j \ne i}\lvert a_{ij} \rvert$ for all $i$, is PD (by Gershgorin's theorem, every eigenvalue lies in a disk centered at some $a_{ii}$ with radius smaller than $a_{ii}$).

The second-derivative test. At a critical point $\mathbf{x}_0$ of a smooth $f$ ($\nabla f(\mathbf{x}_0) = \mathbf{0}$), the Taylor expansion of Section 10.2 gives:

Hessian $\nabla^2f(\mathbf{x}_0)$	Conclusion
positive definite	strict local minimum
negative definite	strict local maximum
indefinite	saddle point
semidefinite and singular	inconclusive: higher-order terms decide

import numpy as np
from itertools import combinations

def leading_minors(A):
    return [np.linalg.det(A[:k, :k]) for k in range(1, len(A) + 1)]

def pivots(A):
    U = A.astype(float).copy()
    for k in range(len(U) - 1):
        U[k+1:] -= np.outer(U[k+1:, k] / U[k, k], U[k])     # eliminate below pivot k (no row exchanges)
    return np.diag(U)

def is_pd_cholesky(A):
    try:
        np.linalg.cholesky(A)
        return True
    except np.linalg.LinAlgError:
        return False

A = np.array([[2., -1., 0.], [-1., 2., -1.], [0., -1., 2.]])
print(np.round(np.linalg.eigvalsh(A), 4))     # [0.5858 2.     3.4142]: 2 - sqrt 2, 2, 2 + sqrt 2
print(np.round(leading_minors(A), 4))         # [2. 3. 4.]
print(np.round(pivots(A), 4))                 # [2.     1.5    1.3333]: ratios of consecutive minors
print(is_pd_cholesky(A))                      # True

# Semidefinite trap: leading minors >= 0 do not imply PSD; all principal minors are needed
B = np.array([[0., 0.], [0., -1.]])
print(np.round(leading_minors(B), 4) + 0.0, np.linalg.eigvalsh(B))   # [0. 0.] [-1.  0.]
def all_principal_minors(A):
    n = len(A)
    return [np.linalg.det(A[np.ix_(S, S)]) for k in range(1, n + 1) for S in combinations(range(n), k)]
print(np.round(all_principal_minors(B), 4) + 0.0)   # [ 0. -1.  0.]: the minor -1 exposes it

# Second-derivative test for f(x, y) = x^3 - 3x + y^2 - 2xy at its critical points
def hessian(x, y):
    return np.array([[6 * x, -2.], [-2., 2.]])
# grad f = (3x^2 - 3 - 2y, 2y - 2x) = 0  ->  y = x and 3x^2 - 2x - 3 = 0
for x in np.roots([3, -2, -3]):
    w = np.linalg.eigvalsh(hessian(x, x))
    kind = "local min" if w.min() > 0 else "local max" if w.max() < 0 else "saddle"
    print(round(x, 4), np.round(w, 4), kind)
# 1.3874 [1.4206 8.9039] local min
# -0.7208 [-4.9039  2.5794] saddle

Key takeaways

• For symmetric $A$: positive energy, positive eigenvalues, positive pivots, positive leading principal minors, an existing Cholesky factor, and $A = B^TB$ with independent columns are all equivalent.
• Pivots are ratios of consecutive leading minors, $d_k = D_k/D_{k-1}$, and have the same signs as the eigenvalues.
• For semidefiniteness, leading minors are not enough; all principal minors (or the smallest eigenvalue) must be checked.
• In numerical code, attempting a Cholesky factorization is the standard PD test; eigenvalues are used when the margin $\lambda_{\min}$ matters.
• The second-derivative test is definiteness of the Hessian at a critical point.

🧪 Practice

1. Use Sylvester's criterion to find all $k$ for which $\begin{bmatrix} k & 1 & 0 \ 1 & k & 1 \ 0 & 1 & k \end{bmatrix}$ is positive definite, and confirm with its eigenvalues $k$ and $k \pm \sqrt2$.
2. Decide whether $\begin{bmatrix} 1 & 2 & 3 \ 2 & 5 & 4 \ 3 & 4 & 9 \end{bmatrix}$ is positive definite by computing its pivots. Which pivot fails, and how does it relate to the leading minors?
3. Find the critical points of $f(x, y) = x^2 + xy + y^2 - 3x$ and of $g(x, y) = x^2 + 4xy + y^2$, and classify them with the Hessian.
4. Interview: You are asked to write is_positive_definite(A) for a production library. Which test do you use, and which edge cases do you handle? Hint: check symmetry first (within a tolerance); attempt a Cholesky factorization and catch the failure (about $n^3/3$ flops, and it stops early); avoid determinants and minors, which overflow, underflow, and lose accuracy; decide what to do near the boundary, where rounding can flip the answer for a matrix whose smallest eigenvalue is about $10^{-16}\lVert A \rVert$; handle NaN and infinite entries.

Cholesky Decomposition

Why it matters. Symmetric positive definite systems $A\mathbf{x} = \mathbf{b}$ are everywhere: normal equations, finite element stiffness systems, Newton steps for convex functions, Gaussian process regression, and Kalman filters. For them, the Cholesky factorization $A = LL^T$ costs half as much as LU, needs no pivoting, and is numerically stable. The same factor also gives log-determinants, samples from Gaussian distributions, and the cheapest reliable test for definiteness.

Analogy. Positive numbers have real square roots: $a = (\sqrt a)^2$. Positive definite matrices have a triangular "square root": $A = LL^T$. Just as $\sqrt a$ fails to exist for negative $a$, the Cholesky algorithm fails, at an attempted square root of a non-positive number, exactly when $A$ is not positive definite.

Theorem (Cholesky factorization). A real symmetric matrix $A$ is positive definite if and only if

$$A = LL^T$$

with $L$ lower triangular with positive diagonal entries. This $L$ is unique.

Proof. If $A \succ 0$, elimination gives $A = L_1DL_1^T$ with positive pivots on the diagonal of $D$ (previous topic). Then $L = L_1D^{1/2}$ works. Conversely, if $A = LL^T$ with $L$ invertible, then $\mathbf{x}^TA\mathbf{x} = \lVert L^T\mathbf{x} \rVert^2 > 0$ for $\mathbf{x} \ne \mathbf{0}$. For uniqueness: if $L_1L_1^T = L_2L_2^T$, then $L_2^{-1}L_1 = L_2^TL_1^{-T}$ is both lower and upper triangular, hence diagonal. Call it $\Delta$. Then $\Delta = \Delta^{-T}$ gives $\Delta^2 = I$, and positive diagonals force $\Delta = I$. $\blacksquare$

The upper triangular form $A = R^TR$ with $R = L^T$ is the same factorization.

The algorithm. Compare entries of $A = LL^T$ column by column. Because $L$ is lower triangular, $a_{ij} = \sum_{k \le j} l_{ik}l_{jk}$ for $i \ge j$, which can be solved for the newest unknowns:

$$l_{jj} = \sqrt{a_{jj} - \sum_{k < j} l_{jk}^2}, \qquad l_{ij} = \frac{1}{l_{jj}}\Big(a_{ij} - \sum_{k < j} l_{ik}l_{jk}\Big) \quad (i > j).$$

If the quantity under a square root is $\le 0$, then $A$ is not positive definite, and the algorithm stops. The cost is about $n^3/3$ flops, half of LU's $2n^3/3$, because symmetry means only half the matrix is computed. No pivoting is needed: $a_{ii} = \sum_k l_{ik}^2$ gives $\lvert l_{ik} \rvert \le \sqrt{a_{ii}}$, so the entries of $L$ cannot grow, and that bounded growth is what makes the algorithm stable.

Worked example. Factor

$$A = \begin{bmatrix} 4 & 2 & -2 \ 2 & 10 & 2 \ -2 & 2 & 6 \end{bmatrix}.$$

• Column 1: $l_{11} = \sqrt4 = 2$, $l_{21} = 2/2 = 1$, $l_{31} = -2/2 = -1$.
• Column 2: $l_{22} = \sqrt{10 - 1^2} = 3$, $l_{32} = (2 - (-1)(1))/3 = 1$.
• Column 3: $l_{33} = \sqrt{6 - (-1)^2 - 1^2} = 2$.

$$L = \begin{bmatrix} 2 & 0 & 0 \ 1 & 3 & 0 \ -1 & 1 & 2 \end{bmatrix}.$$

To solve $A\mathbf{x} = \mathbf{b}$ with $\mathbf{b} = (6, 24, 8)$, first solve $L\mathbf{y} = \mathbf{b}$ by forward substitution: $\mathbf{y} = (3, 7, 2)$. Then solve $L^T\mathbf{x} = \mathbf{y}$ by back substitution: $\mathbf{x} = (1, 2, 1)$.

What the factor is used for.

Task	How $L$ is used
Solve $A\mathbf{x} = \mathbf{b}$ for many right-hand sides	factor once ($n^3/3$), then two triangular solves per $\mathbf{b}$ ($2n^2$)
$\log\det A$	$2\sum_i\log l_{ii}$, with no overflow, since $\det A = \prod l_{ii}^2$
Sample $\mathbf{x} \sim \mathcal{N}(\boldsymbol\mu, \Sigma)$	$\mathbf{x} = \boldsymbol\mu + L\mathbf{z}$ with $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, I)$, since $\operatorname{Cov}(L\mathbf{z}) = LL^T = \Sigma$
Whitening	$L^{-1}(\mathbf{x} - \boldsymbol\mu)$ has identity covariance
Test positive definiteness	the factorization succeeds
Generalized eigenproblem $A\mathbf{v} = \lambda B\mathbf{v}$, $B \succ 0$	$B = LL^T$ reduces it to the symmetric problem $L^{-1}AL^{-T}\mathbf{y} = \lambda\mathbf{y}$
Least squares	$A^TA = R^TR$, with $R$ from QR (Section 9.3); forming $A^TA$ squares the condition number, so prefer QR

Variants. $LDL^T$ avoids square roots. With symmetric pivoting (Bunch-Kaufman, which scipy.linalg.ldl uses), it also handles indefinite matrices. Pivoted Cholesky handles semidefinite and nearly low-rank matrices. Sparse Cholesky reorders the unknowns to limit fill-in (Section 12.1). Gaussian-process codes often add a small "jitter" $\varepsilon I$ before factoring, because kernel matrices of nearby points are positive definite in theory but numerically singular.

Library conventions. np.linalg.cholesky returns the lower factor $L$. scipy.linalg.cholesky returns the upper factor $R = L^T$ unless you pass lower=True. scipy.linalg.cho_factor and cho_solve factor once and solve many times.

import numpy as np
from scipy.linalg import cho_factor, cho_solve, cholesky

def cholesky_lower(A):
    """Return lower triangular L with A = L L^T (column by column)."""
    n = len(A)
    L = np.zeros_like(A, dtype=float)
    for j in range(n):
        s = A[j, j] - L[j, :j] @ L[j, :j]     # what remains of the diagonal entry
        if s <= 0:
            raise ValueError("matrix is not positive definite")
        L[j, j] = np.sqrt(s)
        L[j+1:, j] = (A[j+1:, j] - L[j+1:, :j] @ L[j, :j]) / L[j, j]
    return L

A = np.array([[4., 2., -2.], [2., 10., 2.], [-2., 2., 6.]])
L = cholesky_lower(A)
print(L)
# [[ 2.  0.  0.]
#  [ 1.  3.  0.]
#  [-1.  1.  2.]]
print(np.allclose(L, np.linalg.cholesky(A)))  # True: NumPy returns the lower factor
print(np.allclose(cholesky(A), L.T))          # True: SciPy returns the UPPER factor by default

# Solve A x = b with two triangular solves, reusing the factorization
b = np.array([6., 24., 8.])
c, low = cho_factor(A)
x = cho_solve((c, low), b)
print(x, np.allclose(A @ x, b))               # [1. 2. 1.] True

# log det A = 2 * sum(log diag L): no overflow even when det A itself overflows
print(np.isclose(2 * np.log(np.diag(L)).sum(), np.log(np.linalg.det(A))))   # True
big = np.diag(np.full(400, 50.0))             # det = 50^400 overflows float64
print(np.linalg.det(big), round(2 * np.log(np.diag(np.linalg.cholesky(big))).sum(), 4))
# inf 1564.8092                              (det also emits an overflow warning); 400 * log(50) = 1564.8092

# Sampling N(0, A): x = L z with z ~ N(0, I) has covariance L L^T = A
rng = np.random.default_rng(0)
Z = rng.standard_normal((3, 200_000))
samples = L @ Z
print(np.round(np.cov(samples), 1))           # recovers A up to sampling error
# [[ 4.  2. -2.]
#  [ 2. 10.  2.]
#  [-2.  2.  6.]]

# Failure is the cheapest reliable definiteness test
try:
    np.linalg.cholesky(np.array([[1., 2.], [2., 1.]]))
except np.linalg.LinAlgError as err:
    print("not positive definite:", err)      # not positive definite: Matrix is not positive definite

Key takeaways

• A symmetric matrix is positive definite exactly when it has a Cholesky factorization $A = LL^T$ with positive diagonal; $L$ is unique.
• The algorithm computes $L$ column by column at about $n^3/3$ flops, half the cost of LU, with no pivoting and guaranteed stability.
• A non-positive number under a square root proves that $A$ is not positive definite.
• Uses: SPD solves, $\log\det A = 2\sum\log l_{ii}$, Gaussian sampling $\boldsymbol\mu + L\mathbf{z}$, whitening, and reducing generalized eigenproblems.
• NumPy returns the lower factor, SciPy's cholesky the upper one by default.

🧪 Practice

1. Compute by hand the Cholesky factor of $A = \begin{bmatrix} 9 & 3 & 3 \ 3 & 5 & 1 \ 3 & 1 & 10 \end{bmatrix}$, and use it to solve $A\mathbf{x} = (15, 9, 14)$.
2. Run the algorithm on $\begin{bmatrix} 1 & 2 \ 2 & 1 \end{bmatrix}$ and identify the step where it fails. Relate the failed quantity to a pivot and to a leading minor.
3. Prove that the Cholesky factor of a PD matrix satisfies $\lvert l_{ij} \rvert \le \sqrt{a_{ii}}$. Then explain why this bound makes pivoting unnecessary, in contrast with LU (Section 4.4).
4. Interview: A Gaussian-process regression code computes $K^{-1}\mathbf{y}$ and $\log\det K$ for 5000 training points using np.linalg.inv and np.linalg.det, and it intermittently crashes with "singular matrix" or returns -inf. How would you rewrite it? Hint: factor once with cho_factor, use cho_solve for $K^{-1}\mathbf{y}$ and $2\sum\log l_{ii}$ for $\log\det K$; add a small jitter $\varepsilon I$ (or the noise variance) to the diagonal, because kernel matrices of nearby points are numerically singular.

Constrained Optimization of Quadratic Forms

Why it matters. "Which direction maximizes the variance of the data?", "Along which axis is this material most stressed?", and "Which input does this symmetric system amplify most?" are all questions of the form: maximize $\mathbf{x}^TA\mathbf{x}$ over directions $\mathbf{x}$. Because $Q(c\mathbf{x}) = c^2Q(\mathbf{x})$, the question only makes sense once the length of $\mathbf{x}$ is fixed. The answer is one of the most useful facts in applied linear algebra: the extreme values are the extreme eigenvalues, attained at the corresponding eigenvectors.

Theorem (Extreme values on the unit sphere). Let $A$ be symmetric, with eigenvalues $\lambda_1 \le \lambda_2 \le \cdots \le \lambda_n$ and orthonormal eigenvectors $\mathbf{q}_1, \dots, \mathbf{q}_n$. Then

$$\min_{\lVert \mathbf{x} \rVert = 1}\mathbf{x}^TA\mathbf{x} = \lambda_1 \ \ (\text{at } \mathbf{x} = \pm\mathbf{q}1), \qquad \max{\lVert \mathbf{x} \rVert = 1}\mathbf{x}^TA\mathbf{x} = \lambda_n \ \ (\text{at } \mathbf{x} = \pm\mathbf{q}_n).$$

Moreover, the maximum of $\mathbf{x}^TA\mathbf{x}$ over unit vectors with $\mathbf{x} \perp \mathbf{q}n$ is $\lambda{n-1}$, attained at $\mathbf{q}_{n-1}$. Continuing this way, adding one orthogonality constraint at a time, produces every eigenvalue.

Proof. Substitute $\mathbf{x} = Q\mathbf{y}$. Then $\lVert \mathbf{y} \rVert = \lVert \mathbf{x} \rVert = 1$ and

$$\mathbf{x}^TA\mathbf{x} = \lambda_1y_1^2 + \cdots + \lambda_ny_n^2, \qquad y_1^2 + \cdots + y_n^2 = 1.$$

This is a weighted average of the eigenvalues, with weights $y_i^2 \ge 0$ summing to 1. So it lies between $\lambda_1$ and $\lambda_n$, with equality for $\mathbf{y} = \mathbf{e}_1$ and $\mathbf{y} = \mathbf{e}_n$. The constraint $\mathbf{x} \perp \mathbf{q}_n$ means $y_n = 0$, which removes $\lambda_n$ from the average. $\blacksquare$

The Lagrange multiplier is the eigenvalue. With $\mathcal{L}(\mathbf{x}, \lambda) = \mathbf{x}^TA\mathbf{x} - \lambda(\mathbf{x}^T\mathbf{x} - 1)$, setting $\nabla_{\mathbf{x}}\mathcal{L} = 2A\mathbf{x} - 2\lambda\mathbf{x} = \mathbf{0}$ gives exactly $A\mathbf{x} = \lambda\mathbf{x}$. So the critical points of a quadratic form on the sphere are the unit eigenvectors, and the value there is $\mathbf{x}^TA\mathbf{x} = \lambda\mathbf{x}^T\mathbf{x} = \lambda$. Geometrically, a level curve $Q(\mathbf{x}) = c$ is tangent to the unit circle exactly where the gradient $2A\mathbf{x}$ is parallel to the circle's normal $\mathbf{x}$.

Worked example. For $Q(x, y) = 5x^2 - 4xy + 5y^2$ (eigenvalues $3$ and $7$), parametrize the unit circle as $\mathbf{x} = (\cos\theta, \sin\theta)$:

$$Q(\cos\theta, \sin\theta) = 5 - 4\cos\theta\sin\theta = 5 - 2\sin2\theta.$$

$\theta$	$0°$	$45°$	$90°$	$135°$	$180°$
$Q$	$5$	$3$ (min, along $\mathbf{q}_1 = (1, 1)/\sqrt2$)	$5$	$7$ (max, along $\mathbf{q}_2 = (-1, 1)/\sqrt2$)	$5$

The extreme values are exactly $3 = \lambda_{\min}$ and $7 = \lambda_{\max}$, at the principal axes.

Other constraints.

• Ellipsoid constraint $\mathbf{x}^TB\mathbf{x} = 1$ with $B \succ 0$: Lagrange gives $A\mathbf{x} = \lambda B\mathbf{x}$, the generalized symmetric eigenproblem. Mechanical vibrations ($K\mathbf{v} = \omega^2M\mathbf{v}$ with stiffness $K$ and mass $M$) and Fisher's linear discriminant are of this form. It is solved by reduction with $B = LL^T$ (previous topic), which scipy.linalg.eigh(A, B) does.
• No constraint, with a linear term: $f(\mathbf{x}) = \tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$ has gradient $A\mathbf{x} - \mathbf{b}$ and Hessian $A$. If $A \succ 0$, $f$ is strictly convex with the unique minimizer $A\mathbf{x}^* = \mathbf{b}$ and minimum value $-\tfrac12\mathbf{b}^TA^{-1}\mathbf{b}$. If $A$ is indefinite, $f$ is unbounded below. If $A \succeq 0$ is singular, a minimum exists exactly when $\mathbf{b} \in C(A)$. So solving an SPD system is the same as minimizing a convex quadratic. That is the idea behind the conjugate gradient method (Section 12.2), and least squares is the case $A = X^TX$ (Section 9.5).

Applications. The first principal component of centered data $X$ maximizes the projected variance $\mathbf{w}^T\Sigma\mathbf{w}$ over unit vectors $\mathbf{w}$, so it is the top eigenvector of the covariance matrix (Section 11.2). The principal stresses of a stress tensor are its extreme eigenvalues. And $\max_{\lVert \mathbf{x} \rVert = 1}\lVert A\mathbf{x} \rVert^2 = \lambda_{\max}(A^TA)$ is the squared spectral norm of any matrix, which is where the singular values of Chapter 11 come from.

import numpy as np
from scipy.linalg import eigh

A = np.array([[5., -2.], [-2., 5.]])          # Q(x) = 5x^2 - 4xy + 5y^2

# Brute force over the unit circle vs. the eigenvalue answer
t = np.linspace(0, 2 * np.pi, 100_001)
X = np.vstack([np.cos(t), np.sin(t)])
q = np.einsum('ik,ij,jk->k', X, A, X)          # Q at every sampled point
print(round(q.min(), 6), round(q.max(), 6), np.linalg.eigvalsh(A))   # 3.0 7.0 [3. 7.]
print(np.round(X[:, q.argmax()], 4))           # [ 0.7071 -0.7071]: the maximizer (1, -1)/sqrt(2), up to sign

# The second eigenvalue is the max of Q on the unit sphere, orthogonal to the top eigenvector
S = np.array([[3., 1., 0.], [1., 3., 1.], [0., 1., 3.]])
w, V = np.linalg.eigh(S)
rng = np.random.default_rng(2)
Y = rng.standard_normal((3, 200_000))
Y -= np.outer(V[:, 2], V[:, 2] @ Y)            # remove the component along q_max
Y /= np.linalg.norm(Y, axis=0)
print(np.round(w, 4), round(np.einsum('ik,ij,jk->k', Y, S, Y).max(), 4))   # [1.5858 3.     4.4142] 3.0

# Constraint x^T B x = 1 (an ellipse): generalized eigenproblem A v = lambda B v
B = np.array([[4., 0.], [0., 1.]])
lam, Vg = eigh(A, B)                           # SciPy reduces it with a Cholesky factorization of B
print(np.round(lam, 4), np.round(np.diag(Vg.T @ B @ Vg), 6))   # [1.   5.25] [1. 1.]: v^T B v = 1

# Unconstrained: f(x) = 1/2 x^T A x - b^T x with A positive definite is minimized where A x = b
b = np.array([3., 3.])
x_star = np.linalg.solve(A, b)
f = lambda x: 0.5 * x @ A @ x - b @ x
print(x_star, all(f(x_star) < f(x_star + 1e-3 * d) for d in rng.standard_normal((100, 2))))   # [1. 1.] True

Key takeaways

• On the unit sphere, $\mathbf{x}^TA\mathbf{x}$ ranges exactly over $[\lambda_{\min}, \lambda_{\max}]$, with the extremes at the corresponding eigenvectors.
• In eigen-coordinates, $\mathbf{x}^TA\mathbf{x}$ is a weighted average of the eigenvalues; constraints $\mathbf{x} \perp \mathbf{q}n, \mathbf{q}{n-1}, \dots$ reveal the other eigenvalues in turn.
• The Lagrange multiplier of the constraint $\lVert \mathbf{x} \rVert = 1$ is the eigenvalue itself.
• The constraint $\mathbf{x}^TB\mathbf{x} = 1$ leads to the generalized eigenproblem $A\mathbf{x} = \lambda B\mathbf{x}$.
• Minimizing $\tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$ with $A \succ 0$ is equivalent to solving $A\mathbf{x} = \mathbf{b}$.

🧪 Practice

1. Find the maximum and minimum of $2x^2 + 2y^2 + 2z^2 + 2xy + 2xz + 2yz$ on the unit sphere, and all points where they are attained.
2. For $A = \begin{bmatrix} 3 & 1 & 0 \ 1 & 3 & 1 \ 0 & 1 & 3 \end{bmatrix}$, find the maximum of $\mathbf{x}^TA\mathbf{x}$ over unit vectors orthogonal to the top eigenvector, and a maximizer. (The eigenvalues are $3 - \sqrt2$, $3$, $3 + \sqrt2$.)
3. Write $f(x_1, x_2) = x_1^2 + x_1x_2 + x_2^2 - 3x_1$ as $\tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$, check that $A \succ 0$, and find the minimizer by solving a linear system.
4. Interview: Starting from "find the unit direction along which the projected data has the largest variance", derive why the first principal component of a centered data matrix $X$ is the top eigenvector of $X^TX$. Hint: the projections are $X\mathbf{w}$, their sample variance is $\mathbf{w}^TX^TX\mathbf{w}/(m - 1)$, and maximizing a quadratic form over unit vectors gives $\lambda_{\max}$.

Rayleigh Quotient and Courant-Fischer Theorem

Why it matters. The previous topic characterized the largest and smallest eigenvalues by optimization. The Rayleigh quotient extends this to every nonzero vector, and the Courant-Fischer theorem extends it to every eigenvalue, without mentioning eigenvectors. That shift is powerful. It gives the best eigenvalue estimate from an approximate eigenvector, which drives iterative eigensolvers. It also proves that the eigenvalues of symmetric matrices barely move under perturbations, and it yields interlacing inequalities used in graph theory, vibration analysis, and quantum mechanics.

Definition. For symmetric $A$ and $\mathbf{x} \ne \mathbf{0}$, the Rayleigh quotient is

$$R_A(\mathbf{x}) = \frac{\mathbf{x}^TA\mathbf{x}}{\mathbf{x}^T\mathbf{x}}.$$

It is scale invariant, $R_A(c\mathbf{x}) = R_A(\mathbf{x})$, so it is the quadratic form on the unit sphere extended to all nonzero vectors. By the previous topic, $\lambda_{\min} \le R_A(\mathbf{x}) \le \lambda_{\max}$, and $R_A(\mathbf{q}_i) = \lambda_i$.

Properties.

• Stationary points are eigenvectors. $\nabla R_A(\mathbf{x}) = \dfrac{2}{\mathbf{x}^T\mathbf{x}}\big(A\mathbf{x} - R_A(\mathbf{x}),\mathbf{x}\big)$. This vanishes exactly when $A\mathbf{x} = R_A(\mathbf{x})\mathbf{x}$.
• Best eigenvalue estimate. For a given $\mathbf{x}$, the scalar $\mu$ minimizing $\lVert A\mathbf{x} - \mu\mathbf{x} \rVert$ is $\mu = R_A(\mathbf{x})$. It is the least squares solution of the $n \times 1$ system $\mathbf{x}\mu \approx A\mathbf{x}$ (Section 9.5).
• Quadratic accuracy. Let $\mathbf{x} = \mathbf{q} + \varepsilon\mathbf{d}$, where $A\mathbf{q} = \lambda\mathbf{q}$, $\lVert \mathbf{q} \rVert = \lVert \mathbf{d} \rVert = 1$, and $\mathbf{d} \perp \mathbf{q}$. The cross terms vanish because $\mathbf{q}^TA\mathbf{d} = \lambda\mathbf{q}^T\mathbf{d} = 0$, so

$$R_A(\mathbf{x}) - \lambda = \frac{\varepsilon^2,(\mathbf{d}^TA\mathbf{d} - \lambda)}{1 + \varepsilon^2} = O(\varepsilon^2).$$

An eigenvector that is accurate to 1% gives an eigenvalue accurate to about 0.01%. Rayleigh quotient iteration, which uses $R_A(\mathbf{x})$ as the shift for inverse iteration, therefore converges cubically for symmetric matrices (Section 12.3).

Rayleigh's principle. Since $\lambda_{\min} \le R_A(\mathbf{x})$ for every $\mathbf{x}$, any trial vector gives an upper bound on the smallest eigenvalue. Engineers estimate the fundamental vibration frequency of a structure by plugging in a plausible deflection shape. In quantum mechanics, the same idea is the variational principle: the ground-state energy is at most $\boldsymbol\psi^H\boldsymbol\psi / \boldsymbol\psi^\boldsymbol\psi$ for any trial state.

Theorem (Courant-Fischer min-max). Let $A$ be symmetric with eigenvalues $\lambda_1 \le \lambda_2 \le \cdots \le \lambda_n$. Then, for each $k$,

$$\lambda_k = \min_{\dim S = k}\ \max_{\mathbf{x} \in S,\ \mathbf{x} \ne \mathbf{0}} R_A(\mathbf{x}) = \max_{\dim S = n - k + 1}\ \min_{\mathbf{x} \in S,\ \mathbf{x} \ne \mathbf{0}} R_A(\mathbf{x}),$$

where $S$ ranges over subspaces of $\mathbb{R}^n$ of the given dimension.

In words: on every $k$-dimensional subspace, the Rayleigh quotient reaches at least $\lambda_k$ somewhere, and the span of $\mathbf{q}_1, \dots, \mathbf{q}_k$ is a subspace on which it never exceeds $\lambda_k$. For $k = 1$ and $k = n$, the theorem reduces to the extreme-value theorem of the previous topic.

Proof of the first formula. Let $S$ be any $k$-dimensional subspace, and let $T = \operatorname{span}(\mathbf{q}_k, \dots, \mathbf{q}_n)$, of dimension $n - k + 1$. Since $\dim S + \dim T = n + 1 > n$, they share a nonzero vector $\mathbf{x}$. On $T$, $R_A$ is a weighted average of $\lambda_k, \dots, \lambda_n$, so $R_A(\mathbf{x}) \ge \lambda_k$. Hence $\max_S R_A \ge \lambda_k$ for every $S$. For $S = \operatorname{span}(\mathbf{q}_1, \dots, \mathbf{q}_k)$, $R_A$ is an average of $\lambda_1, \dots, \lambda_k$, so $\max_S R_A = \lambda_k$. The minimum over $S$ is therefore exactly $\lambda_k$. The second formula follows by applying the first to $-A$. $\blacksquare$

Consequences. Each follows by comparing Rayleigh quotients and then applying min-max.

• Weyl's inequality. For symmetric $A$ and $E$, $\lambda_k(A) + \lambda_{\min}(E) \le \lambda_k(A + E) \le \lambda_k(A) + \lambda_{\max}(E)$, so

$$\lvert \lambda_k(A + E) - \lambda_k(A) \rvert \le \lVert E \rVert_2.$$

The eigenvalues of symmetric matrices are perfectly conditioned. This is the property promised in the table of Properties of Symmetric Matrices, and it contrasts with the Jordan block example in Normal Matrices.

• Monotonicity. If $A \preceq B$ (Loewner order), then $\lambda_k(A) \le \lambda_k(B)$ for every $k$. Stiffening a structure raises all of its vibration frequencies.
• Cauchy interlacing. If $B$ is $A$ with its last row and column deleted, with eigenvalues $\mu_1 \le \cdots \le \mu_{n-1}$, then $\lambda_k \le \mu_k \le \lambda_{k+1}$. ($B$ is $A$ restricted to the subspace $x_n = 0$.)

   eigenvalues of A (n = 4):     lambda1        lambda2              lambda3     lambda4
                          ----------*--------------*--------------------*-----------*------->
   eigenvalues of B (3 x 3):               mu1               mu2                mu3
                          -----------------o-----------------o------------------o----------->

                   lambda1 <= mu1 <= lambda2 <= mu2 <= lambda3 <= mu3 <= lambda4

• Rank-one updates. The eigenvalues of $A + \mathbf{v}\mathbf{v}^T$ interlace with those of $A$, every eigenvalue moves up, and the total increase is exactly $\lVert \mathbf{v} \rVert^2 = \operatorname{tr}(\mathbf{v}\mathbf{v}^T)$.

Applications. For a graph Laplacian $L$ (Section 13.2), $\lambda_2 = \min{R_L(\mathbf{x}) : \mathbf{x} \perp \mathbf{1}}$ measures how well connected the graph is, and its minimizer, the Fiedler vector, drives spectral clustering. Courant-Fischer also explains why removing a degree of freedom from a mechanical system (a principal submatrix) cannot push any frequency outside the original interlacing bounds.

import numpy as np

def rayleigh(A, x):
    return (x @ A @ x) / (x @ x)

n = 100
T = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)       # discrete vibrating string
w, V = np.linalg.eigh(T)

# Quadratic accuracy: a vector error of size eps gives an eigenvalue error of size about eps^2
rng = np.random.default_rng(3)
for eps in (1e-2, 1e-3, 1e-4):
    d = rng.standard_normal(n)
    d -= (V[:, 0] @ d) * V[:, 0]               # perturb perpendicular to the eigenvector
    x = V[:, 0] + eps * d / np.linalg.norm(d)
    print(eps, f"{rayleigh(T, x) - w[0]:.1e}")
# 0.01 2.3e-04
# 0.001 1.9e-06
# 0.0001 2.2e-08                              each 10x better vector gives a 100x better eigenvalue

# Rayleigh's principle: any trial shape gives an upper bound on the lowest eigenvalue
j = np.arange(1, n + 1)
trial = j * (n + 1 - j)                        # a parabola that vanishes at both ends
print(f"{w[0]:.6e} {rayleigh(T, trial):.6e}")  # 9.674354e-04 9.802000e-04: within 1.4%, from above

# Weyl: a symmetric perturbation E moves every eigenvalue by at most ||E||_2
E = rng.standard_normal((n, n)) * 1e-3
E = (E + E.T) / 2
shift = np.abs(np.linalg.eigvalsh(T + E) - w).max()
print(shift <= np.linalg.norm(E, 2), f"{shift:.2e} {np.linalg.norm(E, 2):.2e}")   # True 2.28e-03 1.38e-02

# Cauchy interlacing: eigenvalues of a principal submatrix sit between those of the full matrix
S = rng.standard_normal((6, 6)); S = S + S.T
mu = np.linalg.eigvalsh(S[:5, :5])
lam = np.linalg.eigvalsh(S)
print(np.all(lam[:-1] <= mu) and np.all(mu <= lam[1:]))   # True

Key takeaways

• The Rayleigh quotient $R_A(\mathbf{x}) = \mathbf{x}^TA\mathbf{x}/\mathbf{x}^T\mathbf{x}$ lies in $[\lambda_{\min}, \lambda_{\max}]$; its stationary points are the eigenvectors, and its stationary values the eigenvalues.
• $R_A(\mathbf{x})$ is the best eigenvalue estimate for a given $\mathbf{x}$, with error quadratic in the eigenvector error.
• Any trial vector gives an upper bound on $\lambda_{\min}$ (Rayleigh's principle, the variational principle).
• Courant-Fischer: $\lambda_k = \min_{\dim S = k}\max_{\mathbf{x} \in S}R_A(\mathbf{x})$ characterizes every eigenvalue without eigenvectors.
• Consequences: Weyl's inequality $\lvert \lambda_k(A + E) - \lambda_k(A) \rvert \le \lVert E \rVert_2$, monotonicity in the Loewner order, and Cauchy interlacing.

🧪 Practice

1. For $A = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$ and $\mathbf{x} = (1, 0.9)$, compute $R_A(\mathbf{x})$ and compare it with the nearest eigenvalue $3$. Then compare the eigenvalue error with the angle between $\mathbf{x}$ and $(1, 1)$.
2. Use the Courant-Fischer theorem to prove that if $B - A$ is positive semidefinite, then $\lambda_k(A) \le \lambda_k(B)$ for every $k$.
3. The eigenvalues of $\begin{bmatrix} 2 & -1 & 0 \ -1 & 2 & -1 \ 0 & -1 & 2 \end{bmatrix}$ are $2 - \sqrt2$, $2$, $2 + \sqrt2$. Without computing, bound the eigenvalues of its top-left $2 \times 2$ block, and then check your bounds.
4. Interview: A covariance matrix is stored with rounding noise of size about $10^{-12}$ in each entry. Why can you trust its eigenvalues to roughly that accuracy, when for some non-symmetric matrices the same noise shifts eigenvalues by $10^{-6}$? Hint: Weyl's inequality bounds eigenvalue shifts of symmetric matrices by $\lVert E \rVert_2$; for an $n \times n$ Jordan block, a perturbation of size $\varepsilon$ can move eigenvalues by about $\varepsilon^{1/n}$.

<a id="11-singular-value-decomposition-and-matrix-analysis"></a>