{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lesson 16 - One Eigenpair from Iteration\n", "\n", "**Learning goals.** By the end of this lecture, students should be able to:\n", "\n", "- define eigenvalues and eigenvectors using the equation $\\mathbf{A}\\mathbf{x}=\\lambda\\mathbf{x}$,\n", "- explain how repeated multiplication by a matrix can reveal a dominant eigenvector,\n", "- implement power iteration and monitor convergence with eigenvalue changes or residual norms,\n", "- use normalization to keep an eigenvector iteration numerically well scaled,\n", "- describe how inverse iteration shifts the problem from repeated matrix-vector products to repeated linear solves.\n", "\n", "**Key terms.** eigenvalue, eigenvector, power iteration, Rayleigh quotient, inverse iteration.\n", "\n", "The matrix-vector products, norms, residuals, and sparse-iteration vocabulary from Lesson 14 are the computational prerequisites for this chapter. The normalization and geometric language from Lesson 15 help us interpret eigenvectors as directions rather than as uniquely scaled vectors.\n", "\n", "## Starting from Richardson iteration\n", "\n", "Suppose we want to solve the linear system\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{x}=\\mathbf{b}.\n", "$$\n", "\n", "A simple fixed-point idea is Richardson iteration, which, for a given initial $\\mathbf{x}^{(0)}$, produces the sequence\n", "\n", "$$\n", " \\mathbf{x}^{(i+1)}=(\\mathbf{I}-\\mathbf{A})\\mathbf{x}^{(i)}+\\mathbf{b} \\, , \\qquad i = 0, 1, 2, \\ldots\n", "$$\n", "\n", "If $\\mathbf{b}=\\mathbf{0}$, then the update reduces to the homogeneous iteration\n", "\n", "$$\n", " \\mathbf{x}^{(i+1)}=(\\mathbf{I}-\\mathbf{A})\\mathbf{x}^{(i)}.\n", "\\label{eq:homoiter}\n", "$$\n", "\n", "The homogeneous iteration equation looks like a special case, but it is more important than it first appears. If $\\mathbf{x}^{\\ast}$ is the exact solution of $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$, then the error\n", "\n", "$$\n", " \\mathbf{e}^{(i)}=\\mathbf{x}^{(i)}-\\mathbf{x}^{\\ast}\n", "$$\n", "\n", "satisfies\n", "\n", "$$\n", " \\mathbf{e}^{(i+1)}=(\\mathbf{I}-\\mathbf{A})\\mathbf{e}^{(i)}.\n", "$$\n", "\n", "So even when $\\mathbf{b}\\neq \\mathbf{0}$, the behavior of the iteration is governed by repeated multiplication by a matrix.\n", "\n", "**Definition.**\n", "\n", "A nonzero vector $\\mathbf{v}$ is an *eigenvector* of $\\mathbf{A}$ if there exists a scalar $\\lambda$ such that\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}.\n", "$$\n", "\n", "The scalar $\\lambda$ is the corresponding *eigenvalue*. (The German word {\\it eigen} translates to {\\it own} in English, and its use in our present context suggests the eigenpair is \"owned\" by (i.e., is characteristic of) the matrix $\\mathbf{A}$.)\n", "\n", "The key point is that an eigenvector is a direction that the matrix does not rotate away from itself. (Recall that the angle $\\theta$ made between parallel vectors is 0 (or $\\pi$, if antiparallel); equivalently, $\\cos(\\theta)$ has a magnitude of one.) It may stretch, shrink, or flip sign, but it remains on the same line. In that sense, eigenvectors are the *fixed directions* or *modes* of the matrix.\n", "\n", "If $\\mathbf{v}$ is an eigenvector of $\\mathbf{A}$, note that\n", "\n", "$$\n", " (\\mathbf{I}-\\mathbf{A})\\mathbf{v}=\\mathbf{v}-\\lambda\\mathbf{v}=(1-\\lambda)\\mathbf{v}.\n", "$$\n", "\n", "Therefore, $\\mathbf{A}$ and $\\mathbf{I}-\\mathbf{A}$ share the same eigenvectors. The only change is the growth factor: for $\\mathbf{A}$ it is $\\lambda$, and for $\\mathbf{I}-\\mathbf{A}$ it is $1-\\lambda$.\n", "\n", "**Takeaway.**\n", "\n", "In a fixed direction $\\mathbf{v}$, the matrix $\\mathbf{A}$ multiplies the magnitude by $\\lambda$. In the corresponding error iteration, the same direction is multiplied by $1-\\lambda$. That is why eigenvalues control both *mode growth* and *iterative convergence*.\n", "\n", "## Why repeated multiplication reveals a preferred direction\n", "\n", "Consider the repeated multiplication by $\\mathbf{A}$ of a given initial vector $\\mathbf{x}^{(0)}$, or\n", "\n", "$$\n", " \\mathbf{x}^{(k)}=\\mathbf{A}\\mathbf{x}^{(k-1)} \\, , \\qquad k = 1, 2, \\ldots\n", "$$\n", "\n", "Assume for now that $\\mathbf{A}$ has a basis of eigenvectors $\\mathbf{v}_1,\\ldots,\\mathbf{v}_n$, with eigenvalues $\\lambda_1,\\ldots,\\lambda_n$, and that one eigenvalue is dominant in magnitude, i.e.,\n", "\n", "$$\n", " |\\lambda_1|>|\\lambda_2|\\geq \\cdots \\geq |\\lambda_n|.\n", "$$\n", "\n", "Suppose we can write the initial vector as\n", "\n", "$$\n", " \\mathbf{x}^{(0)}=c_1\\mathbf{v}_1+c_2\\mathbf{v}_2+\\cdots+c_n\\mathbf{v}_n \\, .\n", "$$\n", "\n", "After $k$ matrix multiplications, we have\n", "\n", "$$\n", "\\mathbf{x}^{(k)}=\\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.\n", "$$\n", "\n", "By factoring out the dominant term, the result is\n", "\n", "$$\n", " \\mathbf{x}^{(k)}=\\lambda_1^k\\left(c_1\\mathbf{v}_1+\\sum_{j=2}^{n} c_j\\left(\\frac{\\lambda_j}{\\lambda_1}\\right)^k\\mathbf{v}_j\\right).\n", "$$\n", "\n", "Because $|\\lambda_j/\\lambda_1|<1$ for $j\\geq 2$, the ratios decay to zero as $k$ grows. Therefore, provided $c_1\\neq 0$, the direction of $\\mathbf{x}^{(k)}$ approaches the direction of $\\mathbf{v}_1$. One simple way to ensure $c_1 \\neq 0$ almost always is to start with a vector of random numbers, possibly all drawn from $U(0, 1)$.\n", "\n", "This is the motivation for the eigenvalue equation. Iteration tells us that there may be a special direction that remains unchanged except for a scale factor. That observation leads naturally to\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}.\n", "$$\n", "\n", "The vector $\\mathbf{v}$ is the special direction; the scalar $\\lambda$ is the amount of growth or shrinkage induced by the matrix in that direction.\n", "\n", "If the dominant eigenvalue is negative, the iterate flips sign from step to step while still approaching a fixed line. If two or more eigenvalues have the same largest magnitude, or if the dominant pair is complex, the simple picture above must be modified. The one-direction story is cleanest when a single, real dominant eigenvalue exists.\n", "\n", "## Power iteration\n", "\n", "Because repeated multiplication reveals the dominant direction, the simplest numerical procedure for the dominant eigenpair is *power iteration*. To prevent overflow or underflow, we normalize after each multiply, or\n", "\n", "$$\n", " \\mathbf{y}^{(k+1)}=\\mathbf{A}\\mathbf{x}^{(k)} \\, ,\n", "$$\n", "\n", "$$\n", " \\mathbf{x}^{(k+1)}=\\frac{\\mathbf{y}^{(k+1)}}{\\left\\lVert \\mathbf{y}^{(k+1)} \\right\\rVert_2} \\, .\n", "$$\n", "\n", "Once the direction has stabilized, the corresponding eigenvalue can be estimated from the {\\it Rayleigh quotient}\n", "\n", "$$\n", " \\lambda^{(k)}=\\frac{\\mathbf{x}^{(k)\\mathsf{T}}\\mathbf{A}\\mathbf{x}^{(k)}}{\\mathbf{x}^{(k)\\mathsf{T}}\\mathbf{x}^{(k)}} \\, .\n", "$$\n", "\n", "If the iterates are normalized so that $\\left\\lVert \\mathbf{x}^{(k)} \\right\\rVert_2=1$, then this reduces to\n", "\n", "$$\n", " \\lambda^{(k)}=\\mathbf{x}^{(k)\\mathsf{T}}\\mathbf{A}\\mathbf{x}^{(k)}.\n", "$$\n", "\n", "A useful stopping criterion is the residual norm\n", "\n", "$$\n", " r^{(k)}=\\left\\lVert \\mathbf{A}\\mathbf{x}^{(k)}-\\lambda^{(k)}\\mathbf{x}^{(k)} \\right\\rVert_2.\n", "$$\n", "\n", "When $r^{(k)}$ is small, the pair $(\\lambda^{(k)},\\mathbf{x}^{(k)})$ is close to a true eigenpair.\n", "\n", "**Algorithm.**\n", "\n", "**Algorithm: Power iteration for the dominant eigenpair**\n", "\n", "1. Choose a nonzero starting vector $\\mathbf{x}^{(0)}$ and normalize it.\n", "1. For $k=0,1,2,\\ldots$, compute $\\mathbf{y}^{(k+1)}=\\mathbf{A}\\mathbf{x}^{(k)}$.\n", "1. Normalize: $\\mathbf{x}^{(k+1)}=\\mathbf{y}^{(k+1)}/\\left\\lVert \\mathbf{y}^{(k+1)} \\right\\rVert$.\n", "1. Estimate the eigenvalue using the Rayleigh quotient.\n", "1. Stop when the residual $\\left\\lVert \\mathbf{A}\\mathbf{x}^{(k)}-\\lambda^{(k)}\\mathbf{x}^{(k)} \\right\\rVert$ or the change in $\\lambda^{(k)}$ is sufficiently small.\n", "\n", "Under the dominant-eigenvalue assumption, the error in the direction typically decays at a rate proportional to\n", "\n", "$$\n", " \\left|\\frac{\\lambda_2}{\\lambda_1}\\right|^k.\n", "$$\n", "\n", "That is why power iteration is fast when the dominant eigenvalue is well separated from the rest of the spectrum, and slow when the spectral gap is small.\n", "\n", "**Example.**\n", "\n", "Consider\n", "\n", "$$\n", " \\mathbf{A}=\\begin{bmatrix}2 & 1 \\\\ 1 & 2\\end{bmatrix}.\n", "$$\n", "\n", "Its eigenpairs are\n", "\n", "$$\n", " \\lambda_1=3, \\qquad \\mathbf{v}_1=\\begin{bmatrix}1 \\\\ 1\\end{bmatrix},\n", "$$\n", "\n", "$$\n", " \\lambda_2=1, \\qquad \\mathbf{v}_2=\\begin{bmatrix}1 \\\\ -1\\end{bmatrix}.\n", "$$\n", "\n", "If we start from $\\mathbf{x}^{(0)}=[1\\;0]\\mathsf{T}$, then\n", "\n", "$$\n", " \\mathbf{x}^{(0)}=\\frac{1}{2}\\mathbf{v}_1+\\frac{1}{2}\\mathbf{v}_2.\n", "$$\n", "\n", "After $k$ multiplications,\n", "\n", "$$\n", " \\mathbf{A}^k\\mathbf{x}^{(0)}=\\frac{1}{2}3^k\\mathbf{v}_1+\\frac{1}{2}1^k\\mathbf{v}_2.\n", "$$\n", "\n", "As $k$ grows, the $3^k\\mathbf{v}_1$ term dominates, so the iterate aligns with $\\mathbf{v}_1$, and the magnitude grows by a factor of roughly $3$ each step.\n", "\n", "## What the magnitude tells us\n", "\n", "Once the direction has settled onto an eigenvector $\\mathbf{v}$, the action of the matrix becomes especially simple:\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}.\n", "$$\n", "\n", "Applying the matrix again gives\n", "\n", "$$\n", " \\mathbf{A}^2\\mathbf{v}=\\lambda^2\\mathbf{v},\n", "$$\n", "\n", "and after $k$ applications,\n", "\n", "$$\n", " \\mathbf{A}^k\\mathbf{v}=\\lambda^k\\mathbf{v}.\n", "$$\n", "\n", "This is the clearest interpretation of an eigenvalue: it tells us how fast the magnitude changes under repeated application of the matrix, once the direction has become an eigendirection.\n", "\n", "For a simple eigenvalue, the eigendirection is unique up to scaling. Multiplying an eigenvector by any nonzero scalar gives the same direction, so in practice we normalize eigenvectors to a convenient length, often\n", "\n", "$$\n", " \\left\\lVert \\mathbf{v} \\right\\rVert_2=1.\n", "$$\n", "\n", "## Smallest eigenvalue, inverse iteration, and shifted inverse iteration\n", "\n", "### Recovering the smallest magnitude eigenvalue\n", "\n", "If $\\mathbf{A}$ is invertible and\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{v}_j=\\lambda_j\\mathbf{v}_j,\n", "$$\n", "\n", "then multiplying by $\\mathbf{A}^{-1}$ gives\n", "\n", "$$\n", " \\mathbf{A}^{-1}\\mathbf{v}_j=\\frac{1}{\\lambda_j}\\mathbf{v}_j.\n", "$$\n", "\n", "Therefore, the dominant eigenvalue of $\\mathbf{A}^{-1}$ corresponds to the eigenvalue of $\\mathbf{A}$ with smallest magnitude. This means that power iteration applied to $\\mathbf{A}^{-1}$ identifies the eigenvector associated with the smallest $|\\lambda|$.\n", "\n", "In practice, we almost never form $\\mathbf{A}^{-1}$ explicitly. Instead, each iteration solves the linear system\n", "\n", "$$\n", " \\mathbf{A}\\mathbf{y}^{(k+1)}=\\mathbf{x}^{(k)} \\,\n", "$$\n", "\n", "and follows with the normalization\n", "\n", "$$\n", " \\mathbf{x}^{(k+1)}=\\frac{\\mathbf{y}^{(k+1)}}{\\left\\lVert \\mathbf{y}^{(k+1)} \\right\\rVert}.\n", "$$\n", "\n", "### Shifted inverse iteration\n", "\n", "Suppose instead that we want the eigenpair nearest a target value $\\mu$. Consider the shifted matrix\n", "\n", "$$\n", " \\mathbf{A}-\\mu\\mathbf{I}.\n", "$$\n", "\n", "If $\\mathbf{A}\\mathbf{v}_j=\\lambda_j\\mathbf{v}_j$, then\n", "\n", "$$\n", " (\\mathbf{A}-\\mu\\mathbf{I})\\mathbf{v}_j=(\\lambda_j-\\mu)\\mathbf{v}_j,\n", "$$\n", "\n", "and therefore\n", "\n", "$$\n", " (\\mathbf{A}-\\mu\\mathbf{I})^{-1}\\mathbf{v}_j=\\frac{1}{\\lambda_j-\\mu}\\mathbf{v}_j.\n", "$$\n", "\n", "The eigenvalue nearest $\\mu$ produces the largest magnitude reciprocal, so it becomes dominant for the inverse-shifted iteration.\n", "\n", "The resulting algorithm is\n", "\n", "$$\n", " (\\mathbf{A}-\\mu\\mathbf{I})\\mathbf{y}^{(k+1)}=\\mathbf{x}^{(k)},\n", "$$\n", "\n", "$$\n", " \\mathbf{x}^{(k+1)}=\\frac{\\mathbf{y}^{(k+1)}}{\\left\\lVert \\mathbf{y}^{(k+1)} \\right\\rVert}.\n", "$$\n", "\n", "A common refinement is to update $\\mu$ using the latest Rayleigh quotient. This often accelerates convergence substantially.\n", "\n", "**Takeaway.**\n", "\n", "Power iteration finds the largest eigenvalue in *magnitude*. Inverse iteration finds the smallest eigenvalue in *magnitude*, assuming invertibility. Shifted inverse iteration finds the eigenpair *closest to the shift* $\\mu$.\n", "\n", "## Applications of the dominant eigenpair\n", "\n", "### Stability of an iterative scheme\n", "\n", "Return to the error iteration\n", "\n", "$$\n", " \\mathbf{e}^{(i+1)}=(\\mathbf{I}-\\mathbf{A})\\mathbf{e}^{(i)}.\n", "$$\n", "\n", "If $\\mathbf{v}_j$ is an eigenvector of $\\mathbf{A}$ with eigenvalue $\\lambda_j$, then along that mode the error behaves as\n", "\n", "$$\n", " \\mathbf{e}^{(i+1)}=(1-\\lambda_j)\\mathbf{e}^{(i)}.\n", "$$\n", "\n", "Thus, the iteration converges only if every modal factor shrinks:\n", "\n", "$$\n", " |1-\\lambda_j|<1 \\qquad \\text{for all } j.\n", "$$\n", "\n", "Equivalently,\n", "\n", "$$\n", " \\rho(\\mathbf{I}-\\mathbf{A})<1,\n", "$$\n", "\n", "where $\\rho(\\cdot)$ denotes the spectral radius. In other words, eigenvalues predict whether a fixed-point iteration damps or amplifies each error mode.\n", "\n", "### Two-group infinite-medium problem\n", "\n", "A common reactor-physics eigenvalue problem has the form\n", "\n", "$$\n", " \\mathbf{M}\\boldsymbol{\\phi}=\\frac{1}{k_{\\infty}}\\mathbf{F}\\boldsymbol{\\phi},\n", "$$\n", "\n", "where\n", "\n", "$$\n", " \\boldsymbol{\\phi}=\\begin{bmatrix}\\phi_1 \\\\ \\phi_2\\end{bmatrix}\n", "$$\n", "\n", "is the two-group flux shape, $\\mathbf{M}$ contains removal and intergroup terms, and $\\mathbf{F}$ contains production by fission. Multiplying by $\\mathbf{M}^{-1}$ gives\n", "\n", "$$\n", " \\mathbf{M}^{-1}\\mathbf{F}\\boldsymbol{\\phi}=k_{\\infty}\\boldsymbol{\\phi}.\n", "$$\n", "\n", "This is an eigenvalue problem in standard form.\n", "\n", "The interpretation is exactly the same as before:\n", "\n", "- the dominant eigenvector $\\boldsymbol{\\phi}$ gives the asymptotic group-to-group flux shape,\n", "- the dominant eigenvalue $k_{\\infty}$ gives the multiplication factor.\n", "\n", "If $k_{\\infty}>1$, production exceeds losses; if $k_{\\infty}<1$, losses dominate. The dominant eigenpair is therefore the physically important one.\n", "\n", "### Asymptotic mode dominance\n", "\n", "More broadly, whenever a linear model is applied repeatedly,\n", "\n", "$$\n", " \\mathbf{x}^{(k+1)}=\\mathbf{A}\\mathbf{x}^{(k)},\n", "$$\n", "\n", "the long-time state is controlled by the eigenvalue of largest magnitude, provided the assumptions behind power iteration are satisfied. This is why dominant eigenpairs appear in neutron multiplication, repeated transport sweeps, linear stability analysis, and many other settings.\n", "\n", "## Exercises\n", "\n", "1. Verify that $\\mathbf{v}=[1,1]^T$ is an eigenvector of\n", "$$\n", " \\mathbf{A}=\n", " \\begin{bmatrix}\n", " 2 & 1\\\\\n", " 1 & 2\n", " \\end{bmatrix}.\n", "$$\n", "\n", "What is the associated eigenvalue?\n", "\n", "1. Starting from $\\mathbf{x}^{(0)}=[1,0]^T$, decompose the vector into the eigenvector directions $[1,1]^T$ and $[1,-1]^T$. Which component becomes dominant under repeated multiplication by $\\mathbf{A}$?\n", "1. Compute the Rayleigh quotient for the matrix in the first exercise and the normalized vector $\\mathbf{x}=(1/\\sqrt{2})[1,1]^T$.\n", "1. Explain why inverse iteration tends to find the eigenvalue of smallest magnitude.\n", "\n", "**Answers.**\n", "\n", "1. Multiplication gives\n", "$$\n", " \\mathbf{A}\\mathbf{v}=\n", " \\begin{bmatrix}\n", " 3\\\\\n", " 3\n", " \\end{bmatrix}\n", " =\n", " 3\\mathbf{v},\n", "$$\n", "\n", "so the eigenvalue is $\\lambda=3$.\n", "\n", "1. The starting vector can be written as\n", "$$\n", " \\begin{bmatrix}\n", " 1\\\\\n", " 0\n", " \\end{bmatrix}\n", " =\n", " \\frac{1}{2}\n", " \\begin{bmatrix}\n", " 1\\\\\n", " 1\n", " \\end{bmatrix}\n", " +\n", " \\frac{1}{2}\n", " \\begin{bmatrix}\n", " 1\\\\\n", " -1\n", " \\end{bmatrix}.\n", "$$\n", "\n", "The $[1,1]^T$ component becomes dominant because its eigenvalue has magnitude $3$, while the other component has eigenvalue $1$.\n", "\n", "1. The Rayleigh quotient is $3$, matching the exact eigenvalue for that eigenvector.\n", "1. The eigenvalues of $\\mathbf{A}^{-1}$ are $1/\\lambda_j$. The largest reciprocal magnitude corresponds to the smallest $|\\lambda_j|$, so power iteration applied through solves with $\\mathbf{A}$ targets the smallest-magnitude eigenvalue.\n", "\n", "## End-of-lecture summary\n", "\n", "**Takeaway.**\n", "\n", "- Eigenvectors are invariant directions of a matrix.\n", "- Eigenvalues are the corresponding growth, shrinkage, or sign-flip factors.\n", "- The dominant eigenpair can be discovered by iteration because repeated multiplication amplifies the largest-magnitude eigenmode relative to all others.\n", "- Inverse and shifted inverse iteration use linear solves to target the smallest-magnitude eigenvalue or the eigenvalue nearest a chosen shift.\n", "- In reactor physics, the dominant eigenpair carries direct physical meaning: multiplication factor and flux shape.\n" ], "id": "l16-oneeigenpairfromiteration" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" } }, "nbformat": 4, "nbformat_minor": 5 }