Approximating eigenvalues and eigenvectors through iteration — the power method, inverse iteration, the QR algorithm, and the sensitivity of computed eigenvalues
Chapter 10 introduces numerical methods for approximating eigenvalues and eigenvectors of matrices. Exact symbolic eigenvalue computation becomes impractical or impossible for large matrices, so numerical methods are essential in applications across structural engineering, data science, physics, economics, and network analysis. Students will study the power method, normalization, convergence, inverse iteration, shifted inverse iteration, the QR method as a preview, eigenvalue sensitivity, and real-world applications.
The chapter connects to linear algebra concepts students already know — eigenvalues, eigenvectors, and the characteristic polynomial — and shows how those concepts extend into computational practice once matrices grow too large for exact methods to remain feasible.
A large social network has millions of users and billions of connections. How would you identify the most influential nodes in that network? One of the most powerful techniques — used by web search engines, recommendation systems, and biological network analysis — relies entirely on computing a single dominant eigenvalue and its corresponding eigenvector. Now suppose the matrix representing this network is one million by one million. Finding eigenvalues by solving a characteristic polynomial is completely out of the question. Symbolic methods fail. How do numerical mathematicians approach this problem? The answer is iteration. Rather than solving for all eigenvalues at once, numerical methods create a sequence of vectors that gradually align with the dominant eigenvector. This chapter develops those ideas carefully, beginning with the mathematical meaning of eigenvalues and building through the power method, convergence analysis, and extensions into more sophisticated numerical eigenvalue algorithms.
Eigenvalue problems appear throughout science, engineering, and computation. They describe natural frequencies of vibrating structures, stability of equilibrium solutions of differential equations, principal directions in data analysis, ranks of web pages, population dynamics in ecology, energy states in quantum mechanics, and modes of stress in materials science.
In Linear Algebra, students learn to find eigenvalues by computing \(\det(\mathbf{A} - \lambda \mathbf{I}) = 0\) and solving the resulting characteristic polynomial. This works well for small matrices — two-by-two or three-by-three. But for a hundred-by-hundred matrix, the characteristic polynomial has degree one hundred. Computing its roots exactly is neither practical nor numerically stable.
Numerical eigenvalue methods replace exact polynomial root-finding with iterative approximation. They find one eigenvalue at a time, or gradually refine all eigenvalues, in a way that is computationally efficient and numerically reliable for large matrices. This chapter focuses on the ideas a student of introductory numerical methods needs to understand: what the power method does and why it works, how convergence depends on the ratio of eigenvalues, how inverse iteration finds small eigenvalues, and why the QR method is the modern workhorse for complete eigenvalue computation.
By the end of this chapter, you should be able to:
An eigenvalue problem is a question about how a matrix acts on special vectors. For a square matrix \(\mathbf{A}\), a nonzero vector \(\mathbf{v}\) is an eigenvector with eigenvalue \(\lambda\) if
\[ \mathbf{A}\mathbf{v} = \lambda \mathbf{v}. \]
This equation says that multiplying the matrix \(\mathbf{A}\) by the vector \(\mathbf{v}\) produces a scalar multiple of \(\mathbf{v}\). The direction of \(\mathbf{v}\) is preserved; only its magnitude is scaled by \(\lambda\).
Eigenvalues carry profound meaning across mathematics and its applications. In structural mechanics, they describe the natural resonant frequencies of a bridge or aircraft wing. If an external vibration matches a natural frequency, resonance can cause catastrophic failure — the Tacoma Narrows Bridge collapse in 1940 is a historical example. Engineers must compute eigenvalues to ensure safety.
In differential equations, eigenvalues of the Jacobian matrix at an equilibrium determine whether small perturbations grow or decay, telling us whether the equilibrium is stable or unstable. In data science, eigenvalues of covariance matrices identify the directions of greatest variance in a dataset, which is the foundation of principal component analysis. In network science, the largest eigenvalue of an adjacency matrix encodes the most influential node in the network.
The mathematical problem is clear. The computational challenge is that for large matrices, exact methods are not feasible. A \(1000 \times 1000\) matrix has a characteristic polynomial of degree 1000. Finding its roots symbolically or by any polynomial root-finding method is numerically unstable and computationally ruinous. Numerical eigenvalue methods solve this problem by replacing exact root-finding with controlled iteration.
Numerical eigenvalue methods do not try to find a formula. They build a sequence of vectors or matrices that gradually reveal the eigenvalues through repeated, well-controlled computation. This is the same philosophy that runs through every numerical method in this book: replace an impossible exact calculation with a feasible approximate one, and understand precisely how good that approximation is.
Before developing numerical methods, it is worth reviewing the mathematical definitions precisely.
Definition. Let \(\mathbf{A}\) be an \(n \times n\) matrix. A scalar \(\lambda\) is an eigenvalue of \(\mathbf{A}\) if there exists a nonzero vector \(\mathbf{v}\) such that
\[ \mathbf{A}\mathbf{v} = \lambda \mathbf{v}. \]
The vector \(\mathbf{v}\) is called an eigenvector corresponding to \(\lambda\). This equation can be rewritten as
\[ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v} = \mathbf{0}. \]
For this system to have a nonzero solution \(\mathbf{v}\), the matrix \(\mathbf{A} - \lambda \mathbf{I}\) must be singular, which means its determinant must be zero:
\[ \det(\mathbf{A} - \lambda \mathbf{I}) = 0. \]
This equation is called the characteristic equation. The polynomial
\[ p(\lambda) = \det(\mathbf{A} - \lambda \mathbf{I}) \]
is the characteristic polynomial of \(\mathbf{A}\). Its roots are the eigenvalues of \(\mathbf{A}\).
For a \(2 \times 2\) matrix
\[ \mathbf{A} = \begin{pmatrix} a & b \\ c & d \end{pmatrix}, \]
the characteristic polynomial is
\[ p(\lambda) = \lambda^2 - (a + d)\lambda + (ad - bc) = 0. \]
For a \(3 \times 3\) matrix, the characteristic polynomial has degree three. For an \(n \times n\) matrix, it has degree \(n\). When \(n\) is large — say, \(n = 500\) — computing and then finding the roots of a degree-500 polynomial is not a viable strategy. This is where numerical eigenvalue methods become indispensable.
A matrix may have repeated eigenvalues, and eigenvalues may be complex even if the matrix is real. The numerical methods in this chapter focus on the common case of real matrices with distinct real eigenvalues, with appropriate notes about more general situations.
Among all eigenvalues of a matrix, the one with the largest absolute value is called the dominant eigenvalue. It governs the long-run behavior of repeated matrix multiplication.
Definition. Let \(\lambda_1, \lambda_2, \ldots, \lambda_n\) be the eigenvalues of \(\mathbf{A}\), ordered so that
\[ |\lambda_1| \geq |\lambda_2| \geq \cdots \geq |\lambda_n|. \]
If \(|\lambda_1| > |\lambda_2|\), then \(\lambda_1\) is called the dominant eigenvalue, and its corresponding eigenvector \(\mathbf{v}_1\) is the dominant eigenvector.
The dominant eigenvalue is often the most important eigenvalue in applications:
The power method exploits the dominance of \(\lambda_1\) by showing that repeated multiplication by \(\mathbf{A}\) causes any starting vector to align with the dominant eigenvector.
The spectral radius of \(\mathbf{A}\) is \(\rho(\mathbf{A}) = \max_i |\lambda_i|\), the magnitude of the dominant eigenvalue.
The power method is the fundamental iterative algorithm for finding the dominant eigenvalue. Its logic is elegant and its convergence is driven by a simple ratio.
The key idea. Suppose \(\mathbf{A}\) has \(n\) linearly independent eigenvectors \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_n\) with corresponding eigenvalues \(\lambda_1, \lambda_2, \ldots, \lambda_n\) where \(|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n|\).
Any starting vector \(\mathbf{x}^{(0)}\) can be written as a linear combination of eigenvectors:
\[ \mathbf{x}^{(0)} = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \cdots + c_n \mathbf{v}_n, \]
provided \(c_1 \neq 0\) (meaning the starting vector has a nonzero component in the dominant eigenvector direction).
Now multiply repeatedly by \(\mathbf{A}\). Since \(\mathbf{A}\mathbf{v}_i = \lambda_i \mathbf{v}_i\):
\[ \mathbf{A}^k \mathbf{x}^{(0)} = c_1 \lambda_1^k \mathbf{v}_1 + c_2 \lambda_2^k \mathbf{v}_2 + \cdots + c_n \lambda_n^k \mathbf{v}_n. \]
Factor out \(\lambda_1^k\):
\[ \mathbf{A}^k \mathbf{x}^{(0)} = \lambda_1^k \left[ c_1 \mathbf{v}_1 + c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^k \mathbf{v}_2 + \cdots + c_n \left(\frac{\lambda_n}{\lambda_1}\right)^k \mathbf{v}_n \right]. \]
Because \(|\lambda_1| > |\lambda_i|\) for all \(i > 1\), each ratio \(\left|\frac{\lambda_i}{\lambda_1}\right| < 1\). As \(k \to \infty\), all terms \(\left(\frac{\lambda_i}{\lambda_1}\right)^k \to 0\). Therefore:
\[ \mathbf{A}^k \mathbf{x}^{(0)} \approx \lambda_1^k \cdot c_1 \mathbf{v}_1. \]
The iterates align with the dominant eigenvector. The dominant eigenvalue can then be extracted as the ratio of corresponding components of consecutive iterates.
The convergence rate depends on the ratio \(r = \left|\frac{\lambda_2}{\lambda_1}\right|\). When \(r\) is small, convergence is rapid. When \(r\) is close to 1 (the two largest eigenvalues are nearly equal in magnitude), convergence is very slow.
Algorithm: Power Method
Purpose: Approximate the dominant eigenvalue and eigenvector of A.
Inputs: n x n matrix A
Initial vector x (nonzero, chosen arbitrarily or randomly)
Tolerance tol
Maximum iterations N
Steps:
1. Set k = 0.
2. Normalize: x = x / ||x|| (to prevent overflow)
3. Repeat:
a. Compute y = A * x
b. Find the index p such that |y[p]| is largest
c. Set lambda = y[p] (estimated dominant eigenvalue)
d. Normalize: x_new = y / ||y||
e. If ||x_new - x|| < tol, stop and return lambda and x_new
f. Set x = x_new, increment k
4. If k reaches N without convergence, report failure
or continue with best estimate.
Stopping criterion: ||x_new - x|| < tol
Output: Approximate dominant eigenvalue lambda,
approximate dominant eigenvector x
Reliability notes:
- Requires |lambda_1| > |lambda_2| for convergence.
- Convergence rate ~ |lambda_2 / lambda_1|^k.
- Fails if c_1 = 0 (rare in practice with random starts).
- Fails if two eigenvalues tie in absolute value.| Iteration | Eigenvalue Estimate | x₁ | x₂ |
|---|---|---|---|
| 0 | — | 1.000 | 0.000 |
| 1 | 4.000 | 0.894 | 0.447 |
| 2 | 4.472 | 0.789 | 0.614 |
| 3 | 4.770 | 0.734 | 0.679 |
| 5 | 4.923 | 0.711 | 0.703 |
| 8 | 4.984 | 0.7075 | 0.7067 |
| 10 | 4.994 | 0.7072 | 0.7070 |
The eigenvalue estimate converges to 5, and the eigenvector converges to \(\frac{1}{\sqrt{2}}\begin{pmatrix}1\\1\end{pmatrix}\) (the normalized form of \(\begin{pmatrix}1\\1\end{pmatrix}\)).
Check. The true dominant eigenvalue is 5, and \(\mathbf{v}_1 = \begin{pmatrix}1\\1\end{pmatrix}\). The method converges correctly.
Interpret. The power method requires about 10 iterations to achieve four decimal places of accuracy, consistent with the convergence ratio \(r = 0.4\): the error is reduced by a factor of roughly 0.4 each iteration.
Normalization is not merely a convenience — it is essential for preventing numerical overflow and for extracting eigenvalue estimates from the iteration.
Why normalize? Without normalization, the iterates \(\mathbf{A}^k \mathbf{x}^{(0)}\) grow without bound (when \(|\lambda_1| > 1\)) or collapse to zero (when \(|\lambda_1| < 1\)). After just 50 iterations with a dominant eigenvalue of 10, the vector magnitude would be on the order of \(10^{50}\), far beyond what floating-point arithmetic can represent accurately.
How normalization is performed. The most common approach normalizes by the infinity norm (the largest absolute component):
\[ \mathbf{x}^{(k+1)} = \frac{\mathbf{A}\mathbf{x}^{(k)}}{\|\mathbf{A}\mathbf{x}^{(k)}\|_\infty}. \]
Under this normalization, the scaling factor applied at each step converges to the dominant eigenvalue \(\lambda_1\).
Alternatively, one can normalize by the Euclidean norm:
\[ \mathbf{x}^{(k+1)} = \frac{\mathbf{A}\mathbf{x}^{(k)}}{\|\mathbf{A}\mathbf{x}^{(k)}\|_2}, \]
keeping eigenvector iterates as unit vectors.
Extracting the eigenvalue estimate. After normalization by the infinity norm, the normalization constant at each step approximates \(\lambda_1\). After normalization by the Euclidean norm, the eigenvalue is estimated using the Rayleigh quotient:
\[ \lambda \approx \frac{(\mathbf{x}^{(k)})^T \mathbf{A} \mathbf{x}^{(k)}}{(\mathbf{x}^{(k)})^T \mathbf{x}^{(k)}}. \]
The Rayleigh quotient gives a very accurate eigenvalue estimate even when the eigenvector approximation is still somewhat rough. It converges faster than the eigenvector itself.
The Rayleigh quotient is a numerically efficient and robust way to estimate eigenvalues. When you have a good approximation to an eigenvector, the Rayleigh quotient gives you an even better approximation to the corresponding eigenvalue. This asymmetry — eigenvector approximation leads eigenvalue approximation — is a useful practical fact.
The convergence of the power method is governed entirely by the ratio of the two largest eigenvalues in absolute value.
Theorem. If \(|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n|\) and the starting vector has a nonzero component in the \(\mathbf{v}_1\) direction, then the error in the eigenvector iterate decreases like
\[ \text{error} \approx C \left|\frac{\lambda_2}{\lambda_1}\right|^k \]
for some constant \(C\). The eigenvalue estimate converges at twice this rate (because the Rayleigh quotient squares the convergence ratio).
This theorem contains the practical wisdom of the power method:
When the power method converges slowly. If the two largest eigenvalues are nearly equal in magnitude, the power method struggles. Techniques like shifts (discussed in Section 10.8) can dramatically accelerate convergence by making the effective ratio much smaller.
When the power method fails outright. Failure modes include:
Always check whether your eigenvalue estimate is stabilizing as you iterate. If consecutive estimates are not converging toward a fixed value, suspect one of the failure modes above. Plotting the eigenvalue estimate versus iteration number is an easy diagnostic.
The power method finds the dominant eigenvalue — the one with largest absolute value. But many applications require other eigenvalues, particularly the smallest one. Inverse iteration solves this problem by transforming the eigenvalue problem.
The key idea. If \(\mathbf{A}\mathbf{v} = \lambda \mathbf{v}\), then
\[ \mathbf{A}^{-1}\mathbf{v} = \frac{1}{\lambda}\mathbf{v}. \]
The eigenvalues of \(\mathbf{A}^{-1}\) are the reciprocals of the eigenvalues of \(\mathbf{A}\). The smallest eigenvalue of \(\mathbf{A}\) (in absolute value) becomes the largest eigenvalue of \(\mathbf{A}^{-1}\).
Therefore: Applying the power method to \(\mathbf{A}^{-1}\) finds the eigenvalue of \(\mathbf{A}\) with smallest absolute value.
Implementation. Rather than computing \(\mathbf{A}^{-1}\) explicitly (which is expensive and potentially inaccurate), each iteration of inverse iteration solves a linear system:
\[ \mathbf{A}\mathbf{y}^{(k+1)} = \mathbf{x}^{(k)}, \]
then normalizes \(\mathbf{y}^{(k+1)}\) to get \(\mathbf{x}^{(k+1)}\). Solving this linear system at each step (using LU decomposition computed once and stored) is much more efficient than computing \(\mathbf{A}^{-1}\) directly.
Algorithm: Inverse Iteration
Purpose: Approximate the eigenvalue of A with smallest
absolute value.
Inputs: n x n matrix A (nonsingular)
Initial vector x (nonzero)
Tolerance tol
Steps:
1. Compute LU decomposition of A (done once).
2. Repeat:
a. Solve A * y = x using stored LU factorization
b. lambda_inv = component of y with largest magnitude
c. x_new = y / ||y||
d. If ||x_new - x|| < tol, stop
e. Set x = x_new
3. Return eigenvalue estimate: lambda = 1 / lambda_inv
Output: Approximate smallest eigenvalue of A,
approximate eigenvector.Convergence. Inverse iteration converges at rate \(|\lambda_{\text{min}}/\lambda_{\text{second smallest}}|\), which is analogous to the power method but inverted. If the smallest eigenvalue is well-separated from the others, convergence is rapid.
A powerful extension of inverse iteration uses a shift to target any eigenvalue, not just the smallest.
The shift idea. If \(\mu\) is a real number (the shift), then the matrix \(\mathbf{A} - \mu\mathbf{I}\) has eigenvalues \(\lambda_i - \mu\). If \(\mu\) is chosen close to some eigenvalue \(\lambda_j\), then \(\lambda_j - \mu\) is the smallest eigenvalue (in absolute value) of \(\mathbf{A} - \mu\mathbf{I}\).
Applying inverse iteration to \(\mathbf{A} - \mu\mathbf{I}\) therefore finds the eigenvalue of \(\mathbf{A}\) nearest to \(\mu\).
Algorithm: Shifted Inverse Iteration
Purpose: Approximate the eigenvalue of A nearest to a
given shift mu.
Inputs: n x n matrix A
Shift mu (estimate of target eigenvalue)
Initial vector x
Tolerance tol
Steps:
1. Compute LU decomposition of (A - mu * I) (done once).
2. Repeat:
a. Solve (A - mu * I) * y = x
b. x_new = y / ||y||
c. Compute Rayleigh quotient: lambda = x_new^T * A * x_new
d. If ||x_new - x|| < tol, stop
e. Set x = x_new
3. Return eigenvalue estimate lambda and eigenvector x.Convergence. The convergence rate of shifted inverse iteration is
\[ \left|\frac{\lambda_j - \mu}{\lambda_{\text{next nearest}} - \mu}\right|. \]
When \(\mu\) is very close to \(\lambda_j\), the numerator is near zero and convergence becomes extremely rapid — sometimes requiring only a few iterations. This makes shifted inverse iteration extraordinarily powerful when a rough estimate of the target eigenvalue is already known.
Shifted inverse iteration is one of the fastest single-eigenvalue methods available. If you already have an approximate eigenvalue — perhaps from the power method or a rough physical estimate — shifted inverse iteration can refine it to high accuracy in just a few steps.
For computing all eigenvalues of a matrix simultaneously, the modern standard algorithm is the QR method (also called QR iteration or QR algorithm). This section introduces the key ideas without full development.
QR Decomposition Review. Any matrix \(\mathbf{A}\) can be factored as
\[ \mathbf{A} = \mathbf{Q}\mathbf{R}, \]
where \(\mathbf{Q}\) is an orthogonal matrix (\(\mathbf{Q}^T\mathbf{Q} = \mathbf{I}\)) and \(\mathbf{R}\) is upper triangular. This factorization is called the QR decomposition.
The QR Iteration. The QR method builds a sequence of matrices:
\[ \mathbf{A}^{(0)} = \mathbf{A} \]
At each step \(k\):
The key fact is that each \(\mathbf{A}^{(k+1)}\) is similar to \(\mathbf{A}^{(k)}\) (they share the same eigenvalues), but the sequence converges to an upper triangular (or quasi-upper triangular) matrix whose diagonal entries are the eigenvalues of \(\mathbf{A}\).
Why does this work? The QR iteration is related to simultaneous power iteration applied to all eigenvectors at once. The orthogonal factors \(\mathbf{Q}^{(k)}\) accumulate the eigenvectors, while the diagonal of \(\mathbf{R}^{(k)}\mathbf{Q}^{(k)}\) converges to the eigenvalues.
Practical QR methods include shifts (analogous to shifted inverse iteration) that dramatically accelerate convergence, and preprocessing steps that reduce \(\mathbf{A}\) to simpler form (Hessenberg form) before iteration. The full practical QR algorithm with shifts is the method used in virtually all numerical software (MATLAB, NumPy, LAPACK) when eig or similar functions are called.
The full development of the QR algorithm, including convergence theory, Householder reductions, and implicit shifts, belongs to a course in numerical linear algebra or numerical analysis. Here, the goal is for students to understand that:
The QR method is the standard approach for computing all eigenvalues simultaneously.
It is iterative, producing a sequence of similar matrices converging to triangular form.
Shifts are essential for rapid convergence.
The diagonal entries of the converged matrix give all eigenvalues.
Not all eigenvalue problems are equally well-conditioned. The sensitivity of eigenvalues to perturbations in the matrix is a critical practical concern.
The question. If the entries of \(\mathbf{A}\) are perturbed by small amounts (as happens due to measurement error, rounding, or data uncertainty), how much do the eigenvalues change?
A simple bound. The Bauer–Fike theorem states that if \(\mathbf{A}\) is diagonalizable with eigenvector matrix \(\mathbf{V}\) (so \(\mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}\)), and \(\mathbf{E}\) is a perturbation, then each eigenvalue \(\hat\lambda\) of \(\mathbf{A} + \mathbf{E}\) satisfies
\[ \min_j |\hat{\lambda} - \lambda_j| \leq \|\mathbf{V}\|_p \|\mathbf{V}^{-1}\|_p \|\mathbf{E}\|_p = \kappa_p(\mathbf{V}) \|\mathbf{E}\|_p, \]
where \(\kappa_p(\mathbf{V})\) is the condition number of the eigenvector matrix.
What this means. When the eigenvectors of \(\mathbf{A}\) are nearly parallel (nearly linearly dependent), the condition number \(\kappa(\mathbf{V})\) is large, and small perturbations in the matrix can cause large changes in the eigenvalues. Such a matrix is called ill-conditioned for eigenvalue computation.
For symmetric matrices, the eigenvectors are always orthogonal and \(\kappa(\mathbf{V}) = 1\). Eigenvalues of symmetric matrices are therefore always well-conditioned — a fact that makes symmetric eigenvalue problems much more stable than general ones. This is why algorithms for symmetric matrices (like the symmetric QR algorithm or Lanczos iteration) are simpler and more reliable than general eigenvalue algorithms.
Individual eigenvalue condition numbers. For a specific eigenvalue \(\lambda_j\) of a diagonalizable matrix, the eigenvalue condition number is
\[ \kappa_j = \frac{1}{|\mathbf{u}_j^T \mathbf{v}_j|}, \]
where \(\mathbf{v}_j\) is the right eigenvector and \(\mathbf{u}_j\) is the left eigenvector (the eigenvector of \(\mathbf{A}^T\) for the same eigenvalue). When right and left eigenvectors are nearly orthogonal, \(\kappa_j\) is large and the eigenvalue is sensitive.
Before computing eigenvalues numerically, ask whether the matrix is symmetric, nearly symmetric, or general. A symmetric positive definite matrix arising from a physical problem is likely to have well-conditioned eigenvalues. A general matrix assembled from uncertain data may have poorly conditioned eigenvalues, and reported eigenvalues may not be meaningful to the full claimed precision.
One of the most important applications of eigenvalue computation is in the analysis of stability and vibrations of physical systems.
Linear stability of differential equations. Consider a system of differential equations
\[ \frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x}. \]
The general solution is built from terms of the form \(e^{\lambda_i t} \mathbf{v}_i\), where \(\lambda_i\) and \(\mathbf{v}_i\) are the eigenvalues and eigenvectors of \(\mathbf{A}\). The equilibrium \(\mathbf{x} = \mathbf{0}\) is:
Determining stability requires computing or bounding the eigenvalues — specifically, their real parts. This is a fundamental reason eigenvalue computation appears throughout engineering and physics.
Mechanical vibrations. For a mechanical system described by
\[ \mathbf{M}\ddot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{0}, \]
where \(\mathbf{M}\) is the mass matrix and \(\mathbf{K}\) is the stiffness matrix, the natural frequencies \(\omega_i\) of vibration satisfy the generalized eigenvalue problem:
\[ \mathbf{K}\mathbf{v} = \omega^2 \mathbf{M}\mathbf{v}. \]
The natural frequencies are the square roots of the eigenvalues of \(\mathbf{M}^{-1}\mathbf{K}\) (or of the equivalent generalized problem). Computing these frequencies is essential in structural design — an engineer must ensure that no natural frequency coincides with typical operating frequencies that could cause resonance.
Suppose a chemical reactor is modeled by a system of differential equations, and linearization at an operating point gives a \(10 \times 10\) Jacobian matrix \(\mathbf{A}\). The reactor is safe if all eigenvalues of \(\mathbf{A}\) have negative real parts. Computing these ten eigenvalues numerically (using the QR method, for instance) tells the engineer whether the operating point is stable.
Eigenvalue methods play a central role in data analysis and network science.
Google PageRank. The World Wide Web can be modeled as a directed graph. Web pages are nodes, and hyperlinks are edges. Google's PageRank algorithm assigns importance scores to pages by finding the dominant eigenvector of a modified adjacency matrix called the Google matrix \(\mathbf{G}\). The entry \(G_{ij}\) represents the probability that a random web surfer follows a link from page \(j\) to page \(i\). The PageRank vector is the dominant eigenvector of \(\mathbf{G}\), corresponding to the dominant eigenvalue \(\lambda_1 = 1\).
The power method, applied to the Google matrix (which has hundreds of billions of rows and columns), converges surprisingly quickly because the convergence ratio is controlled by the second eigenvalue, and the random-walk structure of the Google matrix keeps the second eigenvalue well below 1.
Principal Component Analysis (PCA). In data science, a dataset with \(p\) measured variables can be represented by a \(p \times p\) covariance matrix \(\mathbf{C}\). The eigenvalues of \(\mathbf{C}\) measure the variance in each principal direction, and the eigenvectors give those directions. The dominant eigenvector points in the direction of greatest variance in the data. Projecting data onto the first few eigenvectors (principal components) reduces dimensionality while preserving most of the variance.
PCA requires computing the eigenvalues and eigenvectors of the covariance matrix. For modern high-dimensional datasets, this is done using specialized numerical algorithms — singular value decomposition (SVD) in particular — which are more stable than forming the covariance matrix explicitly.
Spectral graph theory. The eigenvalues of the Laplacian matrix of a graph (a matrix derived from the graph's connectivity structure) reveal important properties of the graph: the number of connected components, the graph's robustness, how quickly information diffuses through the network, and more. The second-smallest eigenvalue of the Laplacian (the algebraic connectivity or Fiedler value) measures how well-connected the graph is. Computing this eigenvalue requires the inverse iteration or Lanczos methods rather than the basic power method, since the power method finds only the dominant eigenvalue.
Chapter 10 completes Part IV's treatment of numerical methods for problems involving matrices. The arc of Part IV is:
Part V will complete the book with numerical methods for optimization (Chapter 11), ordinary differential equations (Chapter 12), partial differential equations (Chapter 13), and a capstone chapter on scientific computing (Chapter 14).
Chapter 11 — Numerical Optimization introduces methods for finding minimum and maximum values of functions. Eigenvalues return in optimization: the eigenvalues of the Hessian matrix of second derivatives determine whether a critical point is a minimum, maximum, or saddle point. The relationship between linear algebra, eigenvalues, and optimization runs deep throughout computational mathematics.
| Method | Finds | Key Formula | Convergence Rate |
|---|---|---|---|
| Power Method | Dominant eigenvalue | \(\mathbf{x}^{(k+1)} = \mathbf{A}\mathbf{x}^{(k)} / \|\mathbf{A}\mathbf{x}^{(k)}\|\) | \(|\lambda_2/\lambda_1|^k\) |
| Rayleigh Quotient | Refined eigenvalue estimate | \(\lambda \approx \mathbf{x}^T\mathbf{A}\mathbf{x} / \mathbf{x}^T\mathbf{x}\) | Squared rate |
| Inverse Iteration | Smallest eigenvalue | \(\mathbf{A}\mathbf{y}^{(k+1)} = \mathbf{x}^{(k)}\) | \(|\lambda_{\min}/\lambda_{2nd\,\min}|^k\) |
| Shifted Inverse Iteration | Eigenvalue nearest \(\mu\) | \((\mathbf{A}-\mu\mathbf{I})\mathbf{y}^{(k+1)} = \mathbf{x}^{(k)}\) | \(\left|\frac{\lambda_j-\mu}{\lambda_{next}-\mu}\right|^k\) |
| QR Method | All eigenvalues | \(\mathbf{A}^{(k)} = \mathbf{Q}^{(k)}\mathbf{R}^{(k)}\), \(\mathbf{A}^{(k+1)} = \mathbf{R}^{(k)}\mathbf{Q}^{(k)}\) | Depends on shifts used |
Problems 1–5: Eigenvalue computation by hand.
Problems 6–10: Power method iterations.
Chapter 10 developed methods for approximating eigenvalues and eigenvectors when exact computation is impossible. Chapter 11 turns to numerical optimization — finding minimum and maximum values of functions when calculus alone cannot produce a closed-form answer. Students will study one-dimensional search methods, Newton's method for optimization, gradient descent, and an introduction to multivariable and constrained optimization. Eigenvalues return here: the eigenvalues of the Hessian matrix at a critical point determine whether that point is a minimum, maximum, or saddle point, and the ratio of largest to smallest Hessian eigenvalue governs how quickly gradient descent converges — directly mirroring the convergence ratio studied in this chapter.