8 Eigenvalues and Eigenvectors

A matrix, read as a linear map, generally moves a vector to somewhere unrelated: a different direction and a different length. But for almost every matrix there are a few special directions where the map does nothing but stretch.

Those directions are the eigenvectors, the stretch factors are the eigenvalues, and together they are the closest thing linear algebra has to a matrix’s DNA. PCA, covariance structure, the stability of a dynamical system, the convergence rate of gradient descent, and Google’s original PageRank are all the same computation.

8.1 The defining equation

A non-zero vector $\mathbf{v}$ is an eigenvector of a square matrix $\mathbf{A}$, with eigenvalue $\lambda$, when

\[ \mathbf{A}\mathbf{v} = \lambda\mathbf{v} \tag{8.1}\]

Read it carefully: applying the matrix has the same effect as multiplying by a single number. The direction survives; only the length changes.

$\mathbf{v}$ must be non-zero — otherwise every $\lambda$ would work and the definition would say nothing. $\lambda$ may be zero, and a zero eigenvalue is informative: it means $\mathbf{A}\mathbf{v} = \mathbf{0}$, so $\mathbf{v}$ is in the null space and the matrix is singular.

Note also that eigenvectors come in families: if $\mathbf{v}$ is an eigenvector then so is $c\mathbf{v}$ for any $c \neq 0$, with the same eigenvalue. What is really special is the direction, so software returns unit-length representatives.

Worked example. Take

\[ \mathbf{A} = \begin{bmatrix}2 & 1 \\ 1 & 2\end{bmatrix} \]

and try $\mathbf{v} = (1,1)$:

\[ \begin{bmatrix}2 & 1 \\ 1 & 2\end{bmatrix} \begin{bmatrix}1\\1\end{bmatrix} = \begin{bmatrix}3\\3\end{bmatrix} = 3\begin{bmatrix}1\\1\end{bmatrix} \]

An eigenvector with eigenvalue 3. Now try $\mathbf{v} = (1,-1)$:

\[ \begin{bmatrix}2 & 1 \\ 1 & 2\end{bmatrix} \begin{bmatrix}1\\-1\end{bmatrix} = \begin{bmatrix}1\\-1\end{bmatrix} = 1\begin{bmatrix}1\\-1\end{bmatrix} \]

Another one, eigenvalue 1. But $(1,0)$ is not an eigenvector:

A <- matrix(c(2, 1, 1, 2), nrow = 2, byrow = TRUE)
drop(A %*% c(1, 1))

[1] 3 3

drop(A %*% c(1, -1))

[1]  1 -1

drop(A %*% c(1, 0)) # not a multiple of (1, 0)

[1] 2 1

8.2 Geometric meaning

The definition is easier to believe once you see it.

draw_plane(
  vectors = list(
    v = c(1, 1),
    Av = drop(A %*% c(1, 1)),
    u = c(1, 0),
    Au = drop(A %*% c(1, 0))
  ),
  color = c(
    amds_gray, amds_colors[1],
    amds_gray, amds_colors[2]
  ),
  line_type = c("dashed", "solid", "dashed", "solid"),
  label_at = c("mid", "tip", "mid", "tip")
)

Figure 8.1: $\mathbf{A}$ stretches the eigenvector $(1,1)$ by 3 without turning it. The non-eigenvector $(1,0)$ is rotated off its own line — no scalar could reproduce that.

Apply $\mathbf{A}$ to every unit vector at once and the picture becomes clearer still. The unit circle becomes an ellipse, and the eigenvectors are the directions where the ellipse’s axes point.

th <- seq(0, 2 * pi, length.out = 400)
circle <- cbind(cos(th), sin(th))
ellipse <- t(A %*% t(circle))
q1 <- c(1, 1) / sqrt(2)
q2 <- c(1, -1) / sqrt(2)
draw_plane(
  curves = list(circle, ellipse),
  curve_color = c(amds_gray, amds_colors[1]),
  vectors = list(`3q1` = 3 * q1, q2 = q2),
  color = c(amds_colors[2], amds_colors[3]),
  label_at = c("mid", "mid")
)

Figure 8.2: $\mathbf{A}$ maps the unit circle to an ellipse. The long axis lies along the eigenvector with eigenvalue 3, the short axis along the eigenvector with eigenvalue 1.

The eigenvalues are the amounts of stretch along those axes: 3 in one direction, 1 in the other. Every other direction gets some mixture of the two, which is why it comes out rotated.

Going deeper: matrices with no real eigenvectors

A rotation moves every direction, so it should have no eigenvectors at all — and over the real numbers it does not.

Rot <- matrix(c(0, -1, 1, 0), nrow = 2, byrow = TRUE)
eigen(Rot)$values

[1] 0+1i 0-1i

The eigenvalues are $\pm i$. Allowing complex numbers restores the guarantee that an $n \times n$ matrix always has $n$ eigenvalues, and the complex pair encodes the rotation angle. Symmetric matrices — the ones this book mostly cares about — never have this problem (Section 8.7).

8.3 The characteristic polynomial

How do you find the eigenvalues? Rearrange Equation 8.1:

\[ \mathbf{A}\mathbf{v} = \lambda\mathbf{v} \;\Longleftrightarrow\; \mathbf{A}\mathbf{v} - \lambda\mathbf{v} = \mathbf{0} \;\Longleftrightarrow\; (\mathbf{A} - \lambda\mathbf{I})\mathbf{v} = \mathbf{0} \]

We need a non-zero $\mathbf{v}$ solving this. From Section 7.4, a homogeneous system has a non-zero solution exactly when the matrix is singular — and from Section 5.13, singular means zero determinant:

\[ \det(\mathbf{A} - \lambda\mathbf{I}) = 0 \tag{8.2}\]

This is the characteristic equation. Expanding the determinant gives a degree-$n$ polynomial in $\lambda$, the characteristic polynomial, whose roots are the eigenvalues.

Two identities fall out of it, both worth memorizing as sanity checks:

\[ \sum_i \lambda_i = \operatorname{tr}(\mathbf{A}) \qquad \prod_i \lambda_i = \det(\mathbf{A}) \tag{8.3}\]

8.4 Computing eigenvalues by hand

Worked example. For our $\mathbf{A}$:

\[ \mathbf{A} - \lambda\mathbf{I} = \begin{bmatrix}2-\lambda & 1 \\ 1 & 2-\lambda\end{bmatrix} \]

\[ \det(\mathbf{A} - \lambda\mathbf{I}) = (2-\lambda)^2 - 1 = \lambda^2 - 4\lambda + 3 = (\lambda - 3)(\lambda - 1) \]

So $\lambda_1 = 3$ and $\lambda_2 = 1$. Check against Equation 8.3: the trace is 4 and $3 + 1 = 4$; the determinant is 3 and $3 \times 1 = 3$. Both agree.

c(trace = sum(diag(A)), det = det(A))

trace   det 
    4     3

Now find each eigenvector by solving $(\mathbf{A} - \lambda\mathbf{I})\mathbf{v} = \mathbf{0}$.

For $\lambda = 3$:

\[ \begin{bmatrix}-1 & 1 \\ 1 & -1\end{bmatrix} \begin{bmatrix}v_1\\v_2\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix} \]

Both rows say $v_1 = v_2$, so $\mathbf{v} = (1,1)$ — or any multiple. That the two rows say the same thing is not a coincidence: $\mathbf{A} - \lambda\mathbf{I}$ is singular by construction, which is exactly why a non-zero solution exists.

For $\lambda = 1$:

\[ \begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix} \begin{bmatrix}v_1\\v_2\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix} \]

giving $v_1 = -v_2$, so $\mathbf{v} = (1,-1)$.

Watch out

This is how to understand eigenvalues, not how to compute them. Root-finding for polynomials is badly conditioned, and beyond degree 4 there is no formula at all.

Real software never forms the characteristic polynomial. It runs iterative algorithms — the QR algorithm, or variants of Section 8.10 — directly on the matrix. For $n$ beyond about 3, use eigen().

8.5 Eigendecomposition

Collect the eigenvectors as the columns of a matrix $\mathbf{V}$ and the eigenvalues along the diagonal of $\boldsymbol{\Lambda}$. Then Equation 8.1, written for all eigenvectors at once, becomes

\[ \mathbf{A}\mathbf{V} = \mathbf{V}\boldsymbol{\Lambda} \]

and if $\mathbf{V}$ is invertible — which it is when the eigenvectors are independent —

\[ \mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \tag{8.4}\]

This is the eigendecomposition. Read right to left, it says every such matrix does the same three things: change to the eigenvector basis, scale each coordinate independently, change back.

e <- eigen(A)
e$values

[1] 3 1

round(e$vectors, 4)

       [,1]    [,2]
[1,] 0.7071 -0.7071
[2,] 0.7071  0.7071

V <- e$vectors
L <- diag(e$values)
round(V %*% L %*% solve(V), 10)

     [,1] [,2]
[1,]    2    1
[2,]    1    2

R returns unit-length eigenvectors, so $(1,1)$ appears as $(0.707, 0.707)$. Signs are arbitrary — $-\mathbf{v}$ is as valid an eigenvector as $\mathbf{v}$ — so do not be surprised if they flip between machines or versions.

Not every matrix has one. A matrix without $n$ independent eigenvectors is defective and cannot be diagonalized; $\begin{bmatrix}1&1\\0&1\end{bmatrix}$ is the standard example, with eigenvalue 1 repeated but only one eigenvector direction. The SVD (Section 9.6) exists for every matrix, which is one reason it is the more useful tool in practice.

8.6 Diagonalization

The payoff of Equation 8.4 is that powers become trivial. Because the inner factors cancel,

\[ \mathbf{A}^2 = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} = \mathbf{V}\boldsymbol{\Lambda}^2\mathbf{V}^{-1} \]

and in general

\[ \mathbf{A}^k = \mathbf{V}\boldsymbol{\Lambda}^k\mathbf{V}^{-1} \tag{8.5}\]

Raising a diagonal matrix to a power just raises each diagonal entry, so $\mathbf{A}^k$ costs one decomposition and $n$ scalar powers instead of $k$ matrix multiplications.

A10 <- V %*% diag(e$values^10) %*% solve(V)
round(A10)

      [,1]  [,2]
[1,] 29525 29524
[2,] 29524 29525

Same answer as multiplying it out:

Apow <- diag(2)
for (i in 1:10) Apow <- Apow %*% A
Apow

      [,1]  [,2]
[1,] 29525 29524
[2,] 29524 29525

Equation 8.5 also explains long-run behavior, which is what people usually want it for. As $k$ grows, $\lambda^k$ explodes if $|\lambda| > 1$ and vanishes if $|\lambda| < 1$. So the largest eigenvalue eventually dominates everything, and the direction of the system settles onto its eigenvector. Here $3^{10} = 59049$ against $1^{10} = 1$ — after ten steps the second eigenvector contributes almost nothing, which is why every entry of $\mathbf{A}^{10}$ is close to $59049/2$.

In machine learning

Eigenvalues govern stability. A Markov chain converges at a rate set by its second-largest eigenvalue. A recurrent network whose weight matrix has eigenvalues above 1 suffers exploding gradients and below 1 vanishing ones — the same $\lambda^k$, seen from a different direction. Gradient descent’s convergence rate on a quadratic is governed by the ratio of the largest to smallest eigenvalue of the Hessian (Section 16.4).

8.7 Symmetric matrices and the spectral theorem

Symmetric matrices are the well-behaved case, and it is worth knowing exactly how well behaved. The spectral theorem says that if $\mathbf{A} = \mathbf{A}^\top$ then:

All eigenvalues are real — no complex pairs.
Eigenvectors for distinct eigenvalues are orthogonal.
$\mathbf{A}$ is always diagonalizable, never defective.

So the eigenvectors can be chosen orthonormal, and $\mathbf{V}$ becomes an orthogonal matrix $\mathbf{Q}$ with $\mathbf{Q}^{-1} = \mathbf{Q}^\top$ (Section 7.11). The decomposition simplifies to

\[ \mathbf{A} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top \tag{8.6}\]

An inverse replaced by a transpose, and a guarantee that it always exists.

round(crossprod(e$vectors), 10) # Q'Q = I

     [,1] [,2]
[1,]    1    0
[2,]    0    1

round(V %*% L %*% t(V), 10) # QLQ' recovers A

     [,1] [,2]
[1,]    2    1
[2,]    1    2

Our eigenvectors $(1,1)$ and $(1,-1)$ are indeed perpendicular. Compare a non-symmetric matrix, where they are not:

B <- matrix(c(3, 1, 0, 2), nrow = 2, byrow = TRUE)
eb <- eigen(B)
eb$values

[1] 3 2

round(eb$vectors, 4)

     [,1]    [,2]
[1,]    1 -0.7071
[2,]    0  0.7071

round(crossprod(eb$vectors), 4) # off-diagonal is not 0

        [,1]    [,2]
[1,]  1.0000 -0.7071
[2,] -0.7071  1.0000

This matters constantly, because the matrices that arise in statistics are symmetric almost by construction: covariance matrices, correlation matrices, Gram matrices $\mathbf{X}^\top\mathbf{X}$, kernel matrices, Hessians.

In machine learning

PCA is Equation 8.6 applied to a covariance matrix. The eigenvectors are the principal components — orthogonal by the spectral theorem, which is why they give uncorrelated scores — and the eigenvalues are the variance captured along each. “The first PC explains 60% of the variance” means $\lambda_1 / \sum_i \lambda_i = 0.6$.

8.8 Quadratic forms

A quadratic form is a function built from a symmetric matrix:

\[ f(\mathbf{x}) = \mathbf{x}^\top\mathbf{A}\mathbf{x} = \sum_i\sum_j a_{ij}x_ix_j \]

It takes a vector and returns a scalar, and every term is degree 2. For our $\mathbf{A}$:

\[ f(\mathbf{x}) = 2x_1^2 + 2x_1x_2 + 2x_2^2 \]

qform <- function(x, A) drop(crossprod(x, A %*% x))
c(qform(c(1, 1), A), qform(c(1, -1), A), qform(c(1, 0), A))

[1] 6 2 2

Eigenvectors make quadratic forms transparent. Substituting Equation 8.6 and letting $\mathbf{y} = \mathbf{Q}^\top\mathbf{x}$ be the coordinates in the eigenvector basis:

\[ f(\mathbf{x}) = \mathbf{x}^\top\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top\mathbf{x} = \mathbf{y}^\top\boldsymbol{\Lambda}\mathbf{y} = \sum_i \lambda_i y_i^2 \tag{8.7}\]

In the right basis, every cross term disappears. The form is just a weighted sum of squares, weighted by the eigenvalues.

That makes the level sets easy to describe. The set $f(\mathbf{x}) = 1$ is an ellipse whose axes lie along the eigenvectors, with semi-axis length $1/\sqrt{\lambda_i}$.

lev <- (cos(th) / sqrt(3)) %o% q1 + (sin(th) / sqrt(1)) %o% q2
draw_plane(
  curves = list(lev),
  curve_color = amds_colors[1],
  vectors = list(q1 = q1 / sqrt(3), q2 = q2),
  color = c(amds_colors[2], amds_colors[3]),
  label_at = c("tip", "tip")
)

Figure 8.3: The level set $\mathbf{x}^\top\mathbf{A}\mathbf{x} = 1$. Its axes align with the eigenvectors, and it is narrow where the eigenvalue is large — the semi-axis is $1/\sqrt{\lambda}$.

Note the inversion: a large eigenvalue makes the level set narrow in that direction, because $f$ grows quickly there and reaches 1 sooner. This is exactly the picture of an ill-conditioned optimization problem — a long thin valley — and we return to it in Section 16.3.

If some eigenvalues are negative, the level sets are hyperbolas instead:

s <- seq(-1.6, 1.6, length.out = 200)
draw_plane(
  curves = list(
    cbind(cosh(s), sinh(s)),
    cbind(-cosh(s), sinh(s))
  ),
  curve_color = amds_colors[1]
)

Figure 8.4: An indefinite form, $x_1^2 - x_2^2 = 1$. With eigenvalues of opposite sign the level sets are hyperbolas, and the origin is a saddle point rather than a minimum.

8.9 Positive definite matrices

A symmetric matrix $\mathbf{A}$ is:

Name	Condition	Eigenvalues
positive definite	$\mathbf{x}^\top\mathbf{A}\mathbf{x} > 0$ for all $\mathbf{x} \neq \mathbf{0}$	all $> 0$
positive semi-definite	$\mathbf{x}^\top\mathbf{A}\mathbf{x} \geq 0$	all $\geq 0$
indefinite	takes both signs	mixed signs

The equivalence between the two columns is immediate from Equation 8.7: if $f = \sum_i\lambda_iy_i^2$, it is positive for every non-zero $\mathbf{y}$ exactly when every $\lambda_i$ is positive.

mats <- list(
  PD = matrix(c(2, 1, 1, 2), 2, byrow = TRUE),
  PSD = matrix(c(1, 1, 1, 1), 2, byrow = TRUE),
  indefinite = matrix(c(1, 0, 0, -1), 2, byrow = TRUE)
)
sapply(mats, function(M) eigen(M, symmetric = TRUE)$values)

     PD PSD indefinite
[1,]  3   2          1
[2,]  1   0         -1

Positive definiteness is the matrix version of “positive number”, and it is the condition that makes things work:

A positive definite matrix is invertible (no zero eigenvalue).
It has a Cholesky factorization (Section 9.4), the fastest way to solve with it.
A quadratic form with a positive definite matrix has a unique minimum. This is what makes least squares have one answer, and it is the definition of a convex quadratic (Section 15.3).

$\mathbf{X}^\top\mathbf{X}$ is always positive semi-definite, since $\mathbf{x}^\top\mathbf{X}^\top\mathbf{X}\mathbf{x} = \|\mathbf{X}\mathbf{x}\|^2 \geq 0$. It is positive definite exactly when $\mathbf{X}$ has independent columns — the identifiability condition again, in yet another costume.

8.10 Power iteration

The simplest algorithm that actually finds an eigenvector, and the one whose idea powers the serious ones.

Repeatedly multiply by $\mathbf{A}$ and renormalize:

\[ \mathbf{v}_{k+1} = \frac{\mathbf{A}\mathbf{v}_k}{\|\mathbf{A}\mathbf{v}_k\|} \]

Why it works: write the starting vector in the eigenvector basis. Each multiplication scales component $i$ by $\lambda_i$, so after $k$ steps the components are in proportion $c_i\lambda_i^k$. The largest $|\lambda|$ wins by an ever-widening margin, and renormalizing keeps the vector from overflowing.

power_iterate <- function(A, v, n) {
  out <- matrix(0, nrow = length(v), ncol = n + 1L)
  out[, 1] <- v / sqrt(sum(v^2))
  for (k in seq_len(n)) {
    w <- A %*% out[, k]
    out[, k + 1] <- w / sqrt(sum(w^2))
  }
  out
}
iters <- power_iterate(A, c(1, 0), 5)
round(iters, 4)

     [,1]   [,2]   [,3]   [,4]   [,5]   [,6]
[1,]    1 0.8944 0.7809 0.7328 0.7158 0.7100
[2,]    0 0.4472 0.6247 0.6805 0.6983 0.7042

draw_plane(
  curves = list(circle),
  curve_color = amds_gray,
  curve_line_type = "dashed",
  vectors = list(
    v0 = iters[, 1],
    v1 = iters[, 2],
    v2 = iters[, 3],
    v5 = iters[, 6]
  ),
  label_at = "mid"
)

Figure 8.5: Power iteration. Starting from $(1,0)$, each step swings the vector further toward the dominant eigenvector direction $(1,1)/\sqrt{2}$, with the remaining gap shrinking by a factor of $\lambda_2/\lambda_1$ each time.

Once the direction has converged, the eigenvalue comes from the Rayleigh quotient:

\[ \lambda = \frac{\mathbf{v}^\top\mathbf{A}\mathbf{v}}{\mathbf{v}^\top\mathbf{v}} \]

v_final <- iters[, 6]
qform(v_final, A) / sum(v_final^2)

[1] 2.999966

Convergence is governed by $|\lambda_2 / \lambda_1|$ — here $1/3$, so the error shrinks by two-thirds per step. When the top two eigenvalues are close, power iteration crawls, which is why production code uses subspace and Krylov methods instead.

8.11 Eigen in R

ev <- eigen(A, symmetric = TRUE)
ev

eigen() decomposition
$values
[1] 3 1

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

Notes that matter in practice:

Eigenvalues come back in decreasing order, so $values[1] is the largest and $vectors[, 1] matches it.
Pass symmetric = TRUE when you know the matrix is symmetric. R then uses a different algorithm that is faster and guarantees real output — without it, rounding can hand you complex results with tiny imaginary parts.
Use only.values = TRUE when you do not need the vectors; it is substantially cheaper.
For a few eigenvalues of a large sparse matrix, eigen() is the wrong tool — it computes all of them, densely.

eigen(A, symmetric = TRUE, only.values = TRUE)$values

[1] 3 1

8.12 Summary

Idea	Statement
Definition	$\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$, $\mathbf{v} \neq \mathbf{0}$
Finding them	roots of $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$
Sanity checks	$\sum\lambda_i = \operatorname{tr}(\mathbf{A})$, $\prod\lambda_i = \det(\mathbf{A})$
Decomposition	$\mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}$
Symmetric case	$\mathbf{A} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top$, always exists
Powers	$\mathbf{A}^k = \mathbf{V}\boldsymbol{\Lambda}^k\mathbf{V}^{-1}$
Positive definite	all $\lambda_i > 0$
R	`eigen(A, symmetric = TRUE)`

8.13 Exercises

1. Find the eigenvalues of $\begin{bmatrix}4 & 0\\0 & -2\end{bmatrix}$ and $\begin{bmatrix}5 & 2\\0 & 3\end{bmatrix}$ by inspection. What is the rule?

Solution

$4$ and $-2$; then $5$ and $3$. The eigenvalues of a triangular (or diagonal) matrix are its diagonal entries, because $\det(\mathbf{A} - \lambda\mathbf{I})$ of a triangular matrix is the product of its diagonal entries.

eigen(matrix(c(5, 2, 0, 3), 2, byrow = TRUE))$values

[1] 5 3

This is why triangular and diagonal forms are worth aiming for: they hand you the eigenvalues for free.

2. Verify that $\sum\lambda_i = \operatorname{tr}(\mathbf{A})$ and $\prod\lambda_i = \det(\mathbf{A})$ for $\mathbf{B} = \begin{bmatrix}3&1\\0&2\end{bmatrix}$.

Solution

Eigenvalues 3 and 2. Trace $= 5 = 3 + 2$; determinant $= 6 = 3 \times 2$.

c(
  sum_lambda = sum(eb$values), trace = sum(diag(B)),
  prod_lambda = prod(eb$values), det = det(B)
)

 sum_lambda       trace prod_lambda         det 
          5           5           6           6

These are the cheapest error check available on a hand computation.

3. Show that if $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ then $\mathbf{A}^2\mathbf{v} = \lambda^2\mathbf{v}$. What are the eigenvalues of $\mathbf{A}^{-1}$?

Solution

$\mathbf{A}^2\mathbf{v} = \mathbf{A}(\mathbf{A}\mathbf{v}) = \mathbf{A}(\lambda\mathbf{v}) = \lambda(\mathbf{A}\mathbf{v}) = \lambda^2\mathbf{v}$.

For the inverse, multiply $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ by $\mathbf{A}^{-1}$ and divide by $\lambda$ (non-zero, since an invertible matrix has no zero eigenvalue): $\mathbf{A}^{-1}\mathbf{v} = \lambda^{-1}\mathbf{v}$.

rbind(
  A = eigen(A)$values,
  A_squared = eigen(A %*% A)$values,
  # rev() because eigen() sorts descending, and
  # inverting flips which eigenvalue is largest
  A_inverse = rev(eigen(solve(A))$values)
)

               [,1] [,2]
A         3.0000000    1
A_squared 9.0000000    1
A_inverse 0.3333333    1

Eigenvectors are unchanged; eigenvalues follow the function. This generalizes: any polynomial or power series in $\mathbf{A}$ acts on eigenvalues one at a time, which is how $e^{\mathbf{A}}$ is defined.

4. Is $\begin{bmatrix}2&3\\3&2\end{bmatrix}$ positive definite? Answer from the eigenvalues, then find an $\mathbf{x}$ making the quadratic form negative.

Solution

M <- matrix(c(2, 3, 3, 2), 2, byrow = TRUE)
eigen(M, symmetric = TRUE)$values

[1]  5 -1

qform(c(1, -1), M)

[1] -2

Eigenvalues $5$ and $-1$: indefinite, not positive definite. The eigenvector for the negative eigenvalue is $(1,-1)$, and it gives $f = -1 < 0$.

To make a quadratic form negative, look along the eigenvector with the most negative eigenvalue — Equation 8.7 says that is where the form is smallest.

5. Run power iteration on $\mathbf{B} = \begin{bmatrix}3&1\\0&2\end{bmatrix}$ from a starting vector of your choice. Which eigenvector does it find, and how fast?

Solution

itb <- power_iterate(B, c(1, 1), 12)
round(itb[, c(1, 2, 4, 8, 13)], 5)

        [,1]    [,2]    [,3]    [,4]    [,5]
[1,] 0.70711 0.89443 0.98521 0.99955 0.99999
[2,] 0.70711 0.44721 0.17134 0.03013 0.00387

It converges to $(1, 0)$, the eigenvector for the dominant eigenvalue 3. The ratio $|\lambda_2/\lambda_1| = 2/3$ governs the rate, so convergence is noticeably slower than the $1/3$ in Section 8.10 — after 12 steps there is still a visible second component.

Power iteration only ever finds the dominant eigenvector, and it does so slowly when the top two eigenvalues are close.

6. Build a covariance matrix from data, take its eigendecomposition, and confirm that the eigenvalues sum to the total variance.

Solution

set.seed(42)
n <- 200
z1 <- rnorm(n)
z2 <- 0.8 * z1 + 0.6 * rnorm(n)
X <- cbind(z1, z2)
S <- cov(X)
es <- eigen(S, symmetric = TRUE)
round(es$values, 4)

[1] 1.6215 0.1879

c(sum_eigen = sum(es$values), total_var = sum(diag(S)))

sum_eigen total_var 
 1.809394  1.809394

round(es$values / sum(es$values), 4) # proportion explained

[1] 0.8961 0.1039

The eigenvalues sum to the trace, which is the total variance — so “proportion of variance explained” is just $\lambda_i / \sum_j\lambda_j$. That single line is the whole of PCA’s variance accounting.