Eigenvalues & Positive Definite Matrices

The idea in one sentence

Almost every direction gets turned when a matrix acts. A few come back pointing exactly the way they went in, only longer or shorter, and those few decide everything that happens when you apply the matrix again and again.

The directions that survive

Lecture 2 left you with a matrix as a machine: a vector goes in, a moved vector comes out. Now stop watching the whole plane and watch one arrow at a time. Take

and feed it $(-1, 1)$. Out comes $(-3, 1)$. The arrow got longer, and it also swung around toward the horizontal. Output and input point different ways. That is the ordinary case, and it is what almost every arrow does.

Now feed the same machine $(1, 1)$. Out comes $(5, 5)$: the same direction, five times as long. Feed it $(1, -2)$. Out comes $(2, -4)$: the same direction, twice as long. Those two arrows did not turn at all. The machine that bends everything else leaves these two lines exactly where they are and slides points along them.

An arrow that survives like this is an eigenvector of $A$. The number it gets multiplied by is its eigenvalue, written $\lambda$. The word eigen is German for own, so these are the matrix’s own directions. One line holds the whole subject:

The condition $x \neq 0$ is not fussiness. The zero vector satisfies $A0 = \lambda 0$ for every number $\lambda$, so allowing it would make the definition say nothing. Also, if $x$ is an eigenvector then so is $2x$ and so is $-x$, with the same $\lambda$. What survives is a line through the origin, and $x$ is just a name for a point on it.

The eigenvalue is allowed to be zero. If $Ax = 0$ for some nonzero $x$, then $x$ is an eigenvector with $\lambda = 0$, and $x$ sits in the null space $N(A)$. You have already met one. The ring road machine from Lecture 2 sends $(1, 1, 1)$ to $(0, 0, 0)$, because raising all three cities by a meter changes no climb. That blind direction is an eigenvector with eigenvalue zero. Singular and “has eigenvalue zero” are the same news said twice.

How to find them

Guessing works for small matrices and stops working fast. Here is the procedure, and it comes straight out of Lecture 14.

Move everything to one side of $Ax = \lambda x$. You cannot subtract the number $\lambda$ from the matrix $A$, so put it on the diagonal first, as $\lambda I$:

Read that as a question about one matrix. We want a nonzero $x$ in the null space of $A - \lambda I$. A square matrix has a nonzero null vector exactly when it is singular, and Lecture 14 gave a one-number test for singular. So $\lambda$ is an eigenvalue of $A$ precisely when

In words: the eigenvalues of $A$ are the numbers you can subtract from its diagonal to make it singular. That one sentence is the whole method.

Run it on the matrix from Figure 1:

The roots are $\lambda = 5$ and $\lambda = 2$, which is what Figure 1 showed. Now put each root back and solve. For $\lambda = 5$,

whose null space is spanned by $(1, 1)$. For $\lambda = 2$ the matrix is $\begin{bmatrix} 2 & 1 \\ 2 & 1 \end{bmatrix}$ and the null space is spanned by $(1, -2)$. Two lines, two numbers, done.

Look again at the polynomial $\lambda^2 - 7\lambda + 10$. The 7 is $4 + 3$, the sum along the diagonal of $A$. The 10 is $12 - 2$, the determinant of $A$. That is not a coincidence of this example. For any $n$ by $n$ matrix,

where $\mathrm{tr}(A)$ is the sum of the diagonal entries, the trace. Both checks cost seconds and both catch arithmetic slips. Here $5 + 2 = 7$ and $5 \times 2 = 10$.

One more free result. If $A$ is triangular then $A - \lambda I$ is triangular too, and the determinant of a triangular matrix is the product of its diagonal entries. So $\det(A - \lambda I) = (a_{11} - \lambda)\cdots(a_{nn} - \lambda)$, and the eigenvalues of a triangular matrix are exactly the numbers sitting on its diagonal. Nothing to compute.

A machine that settles down

You can find eigenvalues now. Here is why anyone bothers. Two towns swap residents every year. Of the people in town 1, eighty percent stay and twenty percent move to town 2. Of the people in town 2, seventy percent stay and thirty percent move to town 1. Write the fractions as columns:

Each column sums to 1, because everybody ends up somewhere. A matrix like this is called a Markov matrix. If $u_k$ holds the two population fractions in year $k$, then $u_{k+1} = A u_k$, and the population in year 100 is $A^{100} u_0$. Nobody wants to multiply a hundred matrices.

Do the eigenvalue computation instead:

The trace check gives $1 + 0.5 = 1.5$ and the determinant check gives $1 \times 0.5 = 0.5$. Both agree. For $\lambda = 1$, the matrix $A - I$ has rows $(-0.2, 0.3)$ and $(0.2, -0.3)$, so the null space is spanned by $(3, 2)$, which we rescale to $x_1 = (0.6, 0.4)$ so the entries sum to 1. For $\lambda = 0.5$ the null space of $A - 0.5I$ is spanned by $x_2 = (1, -1)$.

Those two numbers say everything. Start from $u_0 = (1, 0)$, everyone in town 1, and split it along the eigenvectors:

Each year, $A$ leaves the first piece alone and halves the second, because that is what eigenvalues 1 and 0.5 do. So

and the second piece dies out. The population settles at 60 percent and 40 percent no matter where it starts. The eigenvector for $\lambda = 1$ is the steady state, and the eigenvector for $\lambda = 0.5$ is the part that decays.

The code below runs the whole example. It checks the two eigenvectors by hand, checks the trace and the determinant, and then computes $A^{100}$ twice: once by a hundred honest multiplications, once from the eigenvalues alone.

import numpy as np

A = np.array([[0.8, 0.3], [0.2, 0.7]])

x1 = np.array([0.6, 0.4])       # null space of A - I
x2 = np.array([1.0, -1.0])      # null space of A - 0.5 I
print(A @ x1)                   # [0.6 0.4]   unchanged, so lambda = 1
print(A @ x2)                   # [ 0.5 -0.5] halved,    so lambda = 0.5

print(np.trace(A), np.linalg.det(A))   # 1.5 0.5  = 1 + 0.5 and 1 * 0.5

def power(M, k):                # the expensive way
    P = M
    for _ in range(k - 1):
        P = P @ M
    return P

print(power(A, 100))            # [[0.6 0.6]
                                #  [0.4 0.4]]

X = np.column_stack([x1, x2])   # eigenvectors in the columns
L = np.diag([1.0 ** 100, 0.5 ** 100])
print(X @ L @ np.linalg.inv(X)) # [[0.6 0.6]
                                #  [0.4 0.4]]   same answer, two numbers of work

import torch

A = torch.tensor([[0.8, 0.3], [0.2, 0.7]])

x1 = torch.tensor([0.6, 0.4])   # null space of A - I
x2 = torch.tensor([1.0, -1.0])  # null space of A - 0.5 I
print(A @ x1)                   # tensor([0.6000, 0.4000])  lambda = 1
print(A @ x2)                   # tensor([ 0.5000, -0.5000]) lambda = 0.5

print(torch.trace(A), torch.linalg.det(A))   # 1.5 0.5

def power(M, k):
    P = M
    for _ in range(k - 1):
        P = P @ M
    return P

print(power(A, 100))            # [[0.6 0.6]
                                #  [0.4 0.4]]

X = torch.stack([x1, x2], dim=1)   # torch spells column_stack as stack(dim=1)
L = torch.diag(torch.tensor([1.0 ** 100, 0.5 ** 100]))
print(X @ L @ torch.linalg.inv(X)) # [[0.6 0.6]
                                   #  [0.4 0.4]]

import jax.numpy as jnp

A = jnp.array([[0.8, 0.3], [0.2, 0.7]])

x1 = jnp.array([0.6, 0.4])      # null space of A - I
x2 = jnp.array([1.0, -1.0])     # null space of A - 0.5 I
print(A @ x1)                   # [0.6 0.4]   unchanged, so lambda = 1
print(A @ x2)                   # [ 0.5 -0.5] halved,    so lambda = 0.5

print(jnp.trace(A), jnp.linalg.det(A))   # 1.5 0.5

def power(M, k):
    P = M
    for _ in range(k - 1):
        P = P @ M
    return P

print(power(A, 100))            # [[0.6 0.6]
                                #  [0.4 0.4]]

X = jnp.column_stack([x1, x2])
L = jnp.diag(jnp.array([1.0 ** 100, 0.5 ** 100]))
print(X @ L @ jnp.linalg.inv(X)) # [[0.6 0.6]
                                 #  [0.4 0.4]]

import tensorflow as tf

A = tf.constant([[0.8, 0.3], [0.2, 0.7]])

x1 = tf.constant([0.6, 0.4])    # null space of A - I
x2 = tf.constant([1.0, -1.0])   # null space of A - 0.5 I
print(tf.linalg.matvec(A, x1))  # [0.6 0.4]   unchanged, so lambda = 1
print(tf.linalg.matvec(A, x2))  # [ 0.5 -0.5] halved,    so lambda = 0.5

print(tf.linalg.trace(A), tf.linalg.det(A))   # 1.5 0.5

def power(M, k):
    P = M
    for _ in range(k - 1):
        P = P @ M
    return P

print(power(A, 100))            # [[0.6 0.6]
                                #  [0.4 0.4]]

X = tf.stack([x1, x2], axis=1)  # tensorflow has no column_stack
L = tf.linalg.diag(tf.constant([1.0 ** 100, 0.5 ** 100]))
print(X @ L @ tf.linalg.inv(X)) # [[0.6 0.6]
                                #  [0.4 0.4]]

Both columns of $A^{100}$ are the steady state. That is the promise of eigenvalues: a hundred matrix multiplications collapse into raising two numbers to the hundredth power, and one of those numbers is so small that it vanishes.

Matrices whose eigenvalues you can read off

Some machines wear their answers on the outside. Four of them are worth memorizing, because they show up inside everything else.

Start with the identity. $Ix = x$ for every $x$, so every vector is an eigenvector and the only eigenvalue is 1. Nothing turns because nothing moves.

Next a projection. Lecture 11 built $P = aa\T / a\T a$, the machine that drops every vector onto the line through $a$. Take $a = (1, 1)$, so

A vector already on the line, like $(1, 1)$, is left alone: $\lambda = 1$. A vector perpendicular to the line, like $(1, -1)$, is crushed to zero: $\lambda = 0$. There is a reason those are the only two options. Projecting twice is the same as projecting once, so $P^2 = P$, and applying that to an eigenvector gives $\lambda^2 = \lambda$. The only numbers that survive are 0 and 1. The column space of $P$ is the eigenvectors for $\lambda = 1$, and the null space is the eigenvectors for $\lambda = 0$.

Double a projection and subtract the identity and you get a reflection, a mirror, $R = 2P - I = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$. Vectors on the mirror line are fixed, so $\lambda = 1$. Vectors perpendicular to it are sent to their opposites, so $\lambda = -1$. Reflecting twice does nothing, so $R^2 = I$ and $\lambda^2 = 1$. The same argument as before, one letter different.

Last, a rotation. Turn the plane a quarter circle, $Q = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$. Every direction moves. No line survives, and the algebra agrees: $\det(Q - \lambda I) = \lambda^2 + 1$, which no real number solves.

Two shortcuts follow from the definition and cost nothing. Apply $A$ twice to an eigenvector and you get $A^2 x = \lambda^2 x$, so squaring a matrix squares its eigenvalues and keeps its eigenvectors. Shift the matrix by $cI$ and every eigenvalue shifts by $c$, since $(A + cI)x = (\lambda + c)x$. The eigenvectors do not move at all.

The shortcut you might expect next is false. Eigenvalues of $A + B$ are not the sums of the eigenvalues, and eigenvalues of $AB$ are not the products, because $A$ and $B$ generally have different eigenvectors. Take $A = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}$ and $B = \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix}$. Every eigenvalue of each is 0. Yet $A + B = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$ has eigenvalues $1$ and $-1$, and $AB = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}$ has eigenvalues 1 and 0. Sums and products of eigenvalues only behave when the two matrices share eigenvectors, which happens exactly when $AB = BA$.

Diagonalization

So far the eigenvectors have been handled one at a time. Take them all at once and the matrix comes apart into numbers.

Suppose an $n$ by $n$ matrix has $n$ independent eigenvectors. Put them in the columns of a matrix $X$ and put the eigenvalues on the diagonal of $\Lambda$. Then look at what $AX$ is, column by column: the $k$-th column of $AX$ is $Ax_k$, which is $\lambda_k x_k$. That is the $k$-th column of $X$ scaled, which is exactly the $k$-th column of $X\Lambda$. So

$X$ is invertible because its columns are independent, which is the only thing we assumed.

Read the middle formula from right to left and it becomes a sentence. $X^{-1}$ rewrites a vector in eigen-coordinates, saying how much of each eigenvector it contains. $\Lambda$ scales each of those amounts by its own eigenvalue. $X$ writes the result back in ordinary coordinates. Change the language, do something trivial, change the language back.

Watch it on the matrix from Figure 1, whose eigenvectors are $x_1 = (1, 1)$ with $\lambda_1 = 5$ and $x_2 = (1, -2)$ with $\lambda_2 = 2$. Take $u = (3, 0)$. In eigen-coordinates, $u = 2x_1 + 1x_2$, since $2(1,1) + (1,-2) = (3, 0)$. Now apply $A$ by scaling each amount:

which is what $A(3,0)$ gives directly. The matrix never had to be multiplied.

Powers come out of this for free. Two copies give $A^2 = X\Lambda X^{-1} X \Lambda X^{-1} = X \Lambda^2 X^{-1}$, since the inner pair cancels, and the same cancellation happens all the way up:

Raising a diagonal matrix to the hundredth power means raising $n$ numbers to the hundredth power. That is the entire trick behind the Markov example, and it also answers the stability question at a glance: $A^k$ goes to the zero matrix exactly when every $|\lambda_i| < 1$.

It does not always work. Diagonalization needs $n$ independent eigenvectors, and some matrices do not have them. Take

It is triangular, so the eigenvalues are 1 and 1. The root $\lambda = 1$ appears twice in the characteristic polynomial, which we call its algebraic multiplicity. But $A - I = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}$ has rank 1, so its null space is only one dimensional: there is a single eigendirection, $(1, 0)$, and no second one. The number of independent eigenvectors for a given $\lambda$ is its geometric multiplicity, and it is never larger than the algebraic multiplicity. When it comes up short there is no eigenvector basis and no diagonalization. Distinct eigenvalues always give independent eigenvectors, so the trouble can only start when a root repeats.

Fibonacci and the golden ratio

Now a matrix where everything works, and where the answer turns out to be famous.

The numbers $0, 1, 1, 2, 3, 5, 8, 13, \dots$ obey $F_{k+2} = F_{k+1} + F_k$. That is a rule for getting the next number from the previous two, which sounds like nothing to do with matrices until you stack two consecutive terms into a vector. Let $u_k = (F_{k+1}, F_k)$. Then

The top row adds the two entries, the bottom row remembers the old top. So $u_k = A^k u_0$ with $u_0 = (1, 0)$, and the Fibonacci numbers are the powers of one small matrix. Now diagonalize it:

The larger root is the golden ratio $\varphi \approx 1.618$, and the smaller is $\psi \approx -0.618$. The trace check gives $\varphi + \psi = 1$ and the determinant check gives $\varphi \psi = -1$, both of which you can read straight off the polynomial. For the eigenvectors, the second row of $A - \lambda I$ says $x_1 = \lambda x_2$, so the eigenvector is $(\lambda, 1)$.

Split $u_0 = (1, 0)$ along them. Solving $c_1(\varphi, 1) + c_2(\psi, 1) = (1, 0)$ gives $c_2 = -c_1$ and $c_1(\varphi - \psi) = 1$, and since $\varphi - \psi = \sqrt{5}$ we get $c_1 = 1/\sqrt5$. Multiply by $A$ a hundred times, which now means multiplying each piece by its own eigenvalue a hundred times, and read the bottom entry:

A sequence defined by addition has a closed form built from square roots. And because $|\psi| < 1$, the second term shrinks to nothing: $F_k$ is the nearest integer to $\varphi^k / \sqrt5$. It also explains a fact people notice without knowing why. The ratio $F_{k+1}/F_k$ marches toward $\varphi$, because the $\varphi$ piece grows and the $\psi$ piece dies, so after a while the vector $u_k$ is pointing almost exactly along the eigenvector $(\varphi, 1)$.

import numpy as np

phi = (1 + 5 ** 0.5) / 2
psi = (1 - 5 ** 0.5) / 2
print(phi, psi)            # 1.618033988749895 -0.6180339887498949
print(phi + psi, phi * psi)  # 1.0 -1.0000000000000002  the trace and the determinant

def fib_closed(n):
    return (phi ** n - psi ** n) / 5 ** 0.5

def fib_loop(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

for n in [10, 30, 50]:
    print(n, fib_loop(n), fib_closed(n))
# 10 55 55.000000000000014
# 30 832040 832040.0000000008
# 50 12586269025 12586269025.00002

A = np.array([[1.0, 1.0], [1.0, 0.0]])
u = np.array([1.0, 0.0])
for _ in range(10):
    u = A @ u
print(u)                     # [89. 55.]  which is (F11, F10)
print(fib_loop(101) / fib_loop(100))   # 1.618033988749895

import torch

phi = (1 + 5 ** 0.5) / 2
psi = (1 - 5 ** 0.5) / 2
print(phi, psi)            # 1.618033988749895 -0.6180339887498949
print(phi + psi, phi * psi)  # 1.0 -1.0000000000000002  the trace and the determinant

def fib_closed(n):
    return (phi ** n - psi ** n) / 5 ** 0.5

def fib_loop(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

for n in [10, 30, 50]:
    print(n, fib_loop(n), fib_closed(n))
# 10 55 55.000000000000014
# 30 832040 832040.0000000008
# 50 12586269025 12586269025.00002

A = torch.tensor([[1.0, 1.0], [1.0, 0.0]], dtype=torch.float64)
u = torch.tensor([1.0, 0.0], dtype=torch.float64)
for _ in range(10):
    u = A @ u
print(u)                     # tensor([89., 55.])  which is (F11, F10)
print(fib_loop(101) / fib_loop(100))   # 1.618033988749895

import jax.numpy as jnp
from jax import config
config.update("jax_enable_x64", True)   # jax defaults to float32, so ask for doubles

phi = (1 + 5 ** 0.5) / 2
psi = (1 - 5 ** 0.5) / 2
print(phi, psi)            # 1.618033988749895 -0.6180339887498949
print(phi + psi, phi * psi)  # 1.0 -1.0000000000000002  the trace and the determinant

def fib_closed(n):
    return (phi ** n - psi ** n) / 5 ** 0.5

def fib_loop(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

for n in [10, 30, 50]:
    print(n, fib_loop(n), fib_closed(n))
# 10 55 55.000000000000014
# 30 832040 832040.0000000008
# 50 12586269025 12586269025.00002

A = jnp.array([[1.0, 1.0], [1.0, 0.0]])
u = jnp.array([1.0, 0.0])
for _ in range(10):
    u = A @ u
print(u)                     # [89. 55.]  which is (F11, F10)
print(fib_loop(101) / fib_loop(100))   # 1.618033988749895

import tensorflow as tf

phi = (1 + 5 ** 0.5) / 2
psi = (1 - 5 ** 0.5) / 2
print(phi, psi)            # 1.618033988749895 -0.6180339887498949
print(phi + psi, phi * psi)  # 1.0 -1.0000000000000002  the trace and the determinant

def fib_closed(n):
    return (phi ** n - psi ** n) / 5 ** 0.5

def fib_loop(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

for n in [10, 30, 50]:
    print(n, fib_loop(n), fib_closed(n))
# 10 55 55.000000000000014
# 30 832040 832040.0000000008
# 50 12586269025 12586269025.00002

A = tf.constant([[1.0, 1.0], [1.0, 0.0]], dtype=tf.float64)
u = tf.constant([1.0, 0.0], dtype=tf.float64)
for _ in range(10):
    u = tf.linalg.matvec(A, u)
print(u)                     # [89. 55.]  which is (F11, F10)
print(fib_loop(101) / fib_loop(100))   # 1.618033988749895

The closed form drifts from the integers in the sixteenth significant digit, which is a double doing its best, and the loop stays exact. Both are worth having. The formula tells you the answer grows like $\varphi^k$; the loop tells you the answer.

That last observation, that repeated multiplication drags any starting vector toward the largest eigenvalue’s direction, is an algorithm. It is called power iteration, and it is how the biggest eigenvalue of a huge matrix gets computed when nobody can afford a characteristic polynomial. Press the button and watch it happen.

Symmetric matrices are the good case

Everything so far has worked for any matrix with enough eigenvectors. One family does much better than that, and it happens to be the family your data arrives in.

A matrix is symmetric when $S = S\T$, so the entry in row $i$ column $j$ equals the entry in row $j$ column $i$. This is not a rare accident. Covariance matrices are symmetric. So is every Gram matrix $A\T A$, and so is the second-derivative matrix of any smooth function. Symmetric matrices come with two gifts.

The first gift is that the eigenvalues are real. No rotations hide inside a symmetric matrix, so the complex case never arises.

The second is that the eigenvectors are perpendicular. Here is the whole proof for two different eigenvalues. Suppose $Sx_1 = \lambda x_1$ and $Sx_2 = \mu x_2$ with $\lambda \neq \mu$. Then

where the middle step is where symmetry is used, since $(Sx_1)\T x_2 = x_1\T S\T x_2$ and $S\T = S$. So $(\lambda - \mu)(x_1 \cdot x_2) = 0$, and since $\lambda \neq \mu$ the dot product must be zero. Perpendicular.

Scale each eigenvector to length 1 and call them $q_1, \dots, q_n$. Put them in the columns of $Q$. Perpendicular unit columns mean $Q\T Q = I$, so $Q^{-1} = Q\T$, and diagonalization loses its inverse:

This is the spectral theorem, and it holds for every symmetric matrix, even when eigenvalues repeat. Read the three factors right to left: $Q\T$ rotates space so the eigenvectors become the coordinate axes, $\Lambda$ stretches along those axes, $Q$ rotates back. Every symmetric matrix is a stretch in disguise, and the disguise is a rotation.

Take $S = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$. The trace is 4 and the determinant is 3, so the eigenvalues are 3 and 1. The eigenvectors are $q_1 = (1,1)/\sqrt2$ and $q_2 = (1,-1)/\sqrt2$, which are perpendicular as promised. Feed $S$ the unit circle and it comes out as an ellipse with half-axes 3 and 1, pointing along $q_1$ and $q_2$.

There is one more way to read $S = Q\Lambda Q\T$, and it is the one that will matter later. Multiply it out as a sum of columns times rows:

Every piece $q_i q_i\T$ is a familiar object. It is Lecture 11’s projection matrix $aa\T / a\T a$ with $a = q_i$, and since $q_i$ has length 1 the denominator is 1. So each piece projects onto one eigendirection, and $S$ scales that projection by $\lambda_i$ and adds the results. A symmetric matrix is a weighted sum of projections onto perpendicular lines. For the example, $q_1 q_1\T = \tfrac12\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$ and $q_2 q_2\T = \tfrac12\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$, and $3$ times the first plus $1$ times the second rebuilds $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ exactly.

Positive definite

Some symmetric matrices have all their eigenvalues on the positive side. That single condition turns out to be the difference between a surface with a bottom and a surface without one, which is why this section is really about optimization.

A symmetric matrix $S$ is positive definite when the number

is greater than zero for every nonzero $x$. That number is called the energy. For a 2 by 2 matrix with rows $(a, b)$ and $(b, c)$ the energy works out to $a x_1^2 + 2b x_1 x_2 + c x_2^2$. So the condition says: this quadratic is positive everywhere except at the origin.

Why is that the same as “all eigenvalues positive”? Because of the spectral theorem. Write $x$ in the eigenvector basis, $x = c_1 q_1 + \dots + c_n q_n$. The $q_i$ are perpendicular unit vectors, so the cross terms vanish and

The energy is a weighted sum of squares with the eigenvalues as weights. It is positive for every nonzero $x$ exactly when every weight is positive. Take one $\lambda_i$ negative and the energy goes negative along $q_i$.

For a symmetric matrix, four tests all say the same thing, and each is worth having because each is cheap in a different situation:

• every eigenvalue is positive;
• every pivot from elimination is positive;
• every leading top-left determinant is positive, the $1\times1$, the $2\times2$, and so on;
• $S = A\T A$ for some matrix $A$ with independent columns.

The last one is the reason positive definite matrices are everywhere in machine learning. If $A$ holds your data, then $x\T A\T A x = (Ax)\T(Ax) = \norm{Ax}^2$, which is a squared length and so never negative, and it is zero only when $Ax = 0$, which independent columns forbid. Least squares and covariance both produce $A\T A$, so both produce positive definite matrices whenever the data has no redundant column.

Try all four on the tridiagonal matrix from Lecture 5:

Its pivots are $2, \tfrac32, \tfrac43$, all positive. Its leading determinants are $2, 3, 4$, all positive. Its eigenvalues are $2 - \sqrt2 \approx 0.586$, then $2$, then $2 + \sqrt2 \approx 3.414$, all positive, and they multiply to 4, which is the determinant, and sum to 6, which is the trace. And Cholesky factors it as $K = LL\T$, which is the fourth test with $A = L\T$. Four tests, one verdict.

Positive semidefinite relaxes every strict inequality to an inclusive one. The energy is allowed to reach zero: $x\T S x \geq 0$ for every $x$, some eigenvalue may be 0, some pivot may be 0, and the determinant may be 0. The two names get swapped constantly, so keep the distinction sharp: definite means the bowl has exactly one lowest point, semidefinite means the bowl may have a flat valley floor running through it.

$N = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$ is the small example. Its energy is $x_1^2 + 2x_1x_2 + x_2^2 = (x_1 + x_2)^2$, which is never negative and is zero all along the line $x_1 + x_2 = 0$. Its eigenvalues are 2 and 0. Semidefinite, not definite. And $\begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}$, with eigenvalues 1 and $-1$, has energy $x_1^2 - x_2^2$, which is positive in one direction and negative in another. That one is called indefinite, and its graph is a saddle.

You have met this test before. In one variable you look at $f''$: positive means the curve holds water and you are sitting at the bottom of it. The Hessian is that same number with more directions to account for, and you read the verdict off the signs exactly as you always did.

import numpy as np

K = np.array([[2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0]])

print(np.linalg.eigvalsh(K))    # [0.58578644 2. 3.41421356]   2 - sqrt2, 2, 2 + sqrt2

L = np.linalg.cholesky(K)       # succeeds only for positive definite input
print(np.diag(L) ** 2)          # [2. 1.5 1.33333333]  the pivots

print(np.round([np.linalg.det(K[:1, :1]),
                np.linalg.det(K[:2, :2]),
                np.linalg.det(K)], 6))            # [2. 3. 4.]  the leading determinants

for x in [np.array([1.0, 2.0, 3.0]), np.array([1.0, 1.0, 1.0])]:
    print(x @ K @ x)            # 12.0 then 2.0   energy, positive away from the origin

lam = np.linalg.eigvalsh(K)
Q = np.linalg.eigh(K)[1]
S = sum(lam[i] * np.outer(Q[:, i], Q[:, i]) for i in range(3))
print(np.round(S, 10))          # rebuilds K from rank-one projections

N = np.array([[1.0, 1.0], [1.0, 1.0]])
print(np.linalg.eigvalsh(N))          # [0. 2.]   semidefinite, not definite
d = np.array([1.0, -1.0])
print(d @ N @ d)                      # 0.0   the flat direction

import torch

K = torch.tensor([[2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0]], dtype=torch.float64)

print(torch.linalg.eigvalsh(K))  # [0.58578644 2. 3.41421356]  2 - sqrt2, 2, 2 + sqrt2

L = torch.linalg.cholesky(K)     # raises only for non positive definite input
print(torch.diag(L) ** 2)        # [2. 1.5 1.33333333]  the pivots

print(torch.round(torch.stack([torch.linalg.det(K[:1, :1]),
                               torch.linalg.det(K[:2, :2]),
                               torch.linalg.det(K)]), decimals=6))   # [2. 3. 4.]

for x in [torch.tensor([1.0, 2.0, 3.0], dtype=torch.float64),
          torch.tensor([1.0, 1.0, 1.0], dtype=torch.float64)]:
    print(x @ K @ x)             # 12.0 then 2.0   energy, positive away from the origin

lam, Q = torch.linalg.eigh(K)
S = sum(lam[i] * torch.outer(Q[:, i], Q[:, i]) for i in range(3))
print(torch.round(S, decimals=10))    # rebuilds K from rank-one projections

N = torch.tensor([[1.0, 1.0], [1.0, 1.0]], dtype=torch.float64)
print(torch.linalg.eigvalsh(N))       # [0. 2.]   semidefinite, not definite
d = torch.tensor([1.0, -1.0], dtype=torch.float64)
print(d @ N @ d)                      # 0.0   the flat direction

import jax.numpy as jnp
from jax import config
config.update("jax_enable_x64", True)

K = jnp.array([[2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0]])

print(jnp.linalg.eigvalsh(K))   # [0.58578644 2. 3.41421356]  2 - sqrt2, 2, 2 + sqrt2

L = jnp.linalg.cholesky(K)      # returns nan entries for non positive definite input
print(jnp.diag(L) ** 2)         # [2. 1.5 1.33333333]  the pivots

print(jnp.round(jnp.array([jnp.linalg.det(K[:1, :1]),
                           jnp.linalg.det(K[:2, :2]),
                           jnp.linalg.det(K)]), 6))       # [2. 3. 4.]

for x in [jnp.array([1.0, 2.0, 3.0]), jnp.array([1.0, 1.0, 1.0])]:
    print(x @ K @ x)            # 12.0 then 2.0   energy, positive away from the origin

lam, Q = jnp.linalg.eigh(K)
S = sum(lam[i] * jnp.outer(Q[:, i], Q[:, i]) for i in range(3))
print(jnp.round(S, 10))         # rebuilds K from rank-one projections

N = jnp.array([[1.0, 1.0], [1.0, 1.0]])
print(jnp.linalg.eigvalsh(N))         # [0. 2.]   semidefinite, not definite
d = jnp.array([1.0, -1.0])
print(d @ N @ d)                      # 0.0   the flat direction

import tensorflow as tf

K = tf.constant([[2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0]], dtype=tf.float64)

print(tf.linalg.eigvalsh(K))    # [0.58578644 2. 3.41421356]  2 - sqrt2, 2, 2 + sqrt2

L = tf.linalg.cholesky(K)       # fills with nan for non positive definite input
print(tf.linalg.diag_part(L) ** 2)   # [2. 1.5 1.33333333]  the pivots
                                     # tensorflow splits diag into diag and diag_part

print(tf.round(tf.stack([tf.linalg.det(K[:1, :1]),
                         tf.linalg.det(K[:2, :2]),
                         tf.linalg.det(K)]) * 1e6) / 1e6)   # [2. 3. 4.]

for x in [tf.constant([1.0, 2.0, 3.0], dtype=tf.float64),
          tf.constant([1.0, 1.0, 1.0], dtype=tf.float64)]:
    print(tf.tensordot(x, tf.linalg.matvec(K, x), 1))   # 12.0 then 2.0   the energy

lam, Q = tf.linalg.eigh(K)
S = sum(lam[i] * tf.tensordot(Q[:, i], Q[:, i], axes=0) for i in range(3))
print(S)                        # rebuilds K from rank-one projections
                                # tensorflow has no outer, so tensordot with axes=0

N = tf.constant([[1.0, 1.0], [1.0, 1.0]], dtype=tf.float64)
print(tf.linalg.eigvalsh(N))          # [0. 2.]   semidefinite, not definite
d = tf.constant([1.0, -1.0], dtype=tf.float64)
print(tf.tensordot(d, tf.linalg.matvec(N, d), 1))   # 0.0   the flat direction

Why gradient descent zigzags

Now the payoff. Training a model means walking downhill on a loss surface. Near a minimum, every smooth loss looks like a quadratic, and the quadratic is written with a symmetric matrix, the matrix of second derivatives called the Hessian. So study the simplest possible loss, with $S$ symmetric and positive definite:

Its lowest point is $x = 0$ and its gradient is $Sx$. The level sets $f(x) = c$ are ellipses, which Figure 7 already showed you: axes along the eigenvectors, and along $q_i$ the ellipse reaches out to $\sqrt{2c/\lambda_i}$. A big eigenvalue means a tight, steeply curved direction. A small eigenvalue means a long, nearly flat one.

Gradient descent takes the step $x \leftarrow x - \alpha S x$. Write $x$ in the eigenvector basis again and the update falls apart into $n$ separate scalar updates:

One number per direction, no interaction. That one line explains everything people say about learning rates. For the method to converge you need $|1 - \alpha\lambda_i| < 1$ in every direction, so $\alpha < 2/\lambda_{\max}$: the steepest direction sets the speed limit. But the error along the flattest direction shrinks by $|1 - \alpha\lambda_{\min}|$, which sits close to 1 when $\lambda_{\min}$ is small. The steep direction caps your step; the flat direction decides how many steps you take. The ratio

is the condition number, and it is the number that matters. With the best possible constant step $\alpha = 2/(\lambda_{\min} + \lambda_{\max})$, the error shrinks by a factor $(\kappa - 1)/(\kappa + 1)$ each step, which is close to 1 whenever $\kappa$ is large.