{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-1 fw=all -->\n",
    "\n",
    "# Eigendecompositions\n",
    "\n",
    "In that section we saw a matrix as a\n",
    "*geometric distortion* of space: it skews, rotates, and rescales the grid, and\n",
    "the determinant records the net volume factor. Eigenvalues refine that picture.\n",
    "For a well-behaved square matrix there is a special set of directions, the\n",
    "*eigenvectors*, along which the distortion is a pure stretch, and an\n",
    "*eigenbasis* in which the whole map decouples into independent one-dimensional\n",
    "stretches. This single observation is what makes stability analysis, PCA, and\n",
    "the curvature story of optimization tractable. The goal of this section is to\n",
    "convey, with both pictures and proofs, why eigenvalues are so central. Our\n",
    "running example is an *iterated map*: repeatedly applying the same matrix, as a\n",
    "deep *linear* network (one with the nonlinearities stripped away) does, makes\n",
    "the role of the largest eigenvalue unmistakable.\n",
    "\n",
    "The numerical checks in this section use small dense matrices."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-imports"
   },
   "source": [
    "import numpy as np"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-2 fw=all -->\n",
    "\n",
    "## Eigenvalues and Eigenvectors\n",
    "\n",
    "### Definition and Geometry\n",
    "\n",
    "Suppose that we have a matrix $\\mathbf{A}$ with the following entries:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "2 & 0 \\\\\n",
    "0 & -1\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "If we apply $\\mathbf{A}$ to any vector $\\mathbf{v} = [x, y]^\\top$,\n",
    "we obtain a vector $\\mathbf{A}\\mathbf{v} = [2x, -y]^\\top$.\n",
    "This has an intuitive interpretation:\n",
    "stretch the vector to be twice as wide in the $x$-direction,\n",
    "and then flip it in the $y$-direction.\n",
    "\n",
    "However, there are *some* vectors for which the direction remains unchanged.\n",
    "Namely $[1, 0]^\\top$ gets sent to $[2, 0]^\\top$\n",
    "and $[0, 1]^\\top$ gets sent to $[0, -1]^\\top$.\n",
    "These vectors are still on the same line through the origin;\n",
    "the only modification is that $\\mathbf{A}$ scales $[1,0]^\\top$ by $2$\n",
    "and $[0,1]^\\top$ by $-1$.\n",
    "We call such vectors *eigenvectors*\n",
    "and the factor they are stretched by their *eigenvalues*.\n",
    "\n",
    "In general, if we can find a number $\\lambda$\n",
    "and a *nonzero* vector $\\mathbf{v}$ such that\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\mathbf{v} = \\lambda \\mathbf{v},\n",
    "$$\n",
    "\n",
    "we say that $\\mathbf{v}$ is an *eigenvector* of $\\mathbf{A}$ and $\\lambda$ its\n",
    "*eigenvalue*. (We insist $\\mathbf{v}\\neq\\mathbf 0$: the zero vector trivially\n",
    "satisfies the equation for every $\\lambda$ and carries no information.\n",
    "The scalar $\\lambda$ *is* allowed to be zero, and we will see that $\\lambda=0$\n",
    "signals a non-invertible matrix.) For a fixed eigenvalue $\\lambda$, the set of\n",
    "all vectors satisfying $\\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}$, together with\n",
    "$\\mathbf 0$, forms a subspace, the *eigenspace* of $\\lambda$.\n",
    "\n",
    "To picture eigenvectors, ask what a matrix does to the\n",
    "*unit circle*, the set of all unit vectors. A full-rank matrix maps it to an\n",
    "ellipse (a singular matrix flattens the ellipse further, into a segment or a\n",
    "point), and for a *symmetric* matrix the axes of that ellipse lie exactly along\n",
    "the eigenvectors, with half-lengths $|\\lambda_i|$ (an axis flips when\n",
    "$\\lambda_i<0$). the figure draws this for our\n",
    "$\\operatorname{diag}(2,-1)$ above and for the symmetric\n",
    "$[[2,1],[1,2]]$ that we revisit in the exercises. This is the same \"circle becomes an ellipse\"\n",
    "picture that the singular value decomposition will generalize to *every* matrix\n",
    "in that section; here the special feature is that one set of\n",
    "axes does the whole job.\n",
    "\n",
    "![The unit circle maps to an ellipse. For a symmetric matrix the ellipse axes lie along the eigenvectors (green), with images scaled by the eigenvalues (red). Left: $\\operatorname{diag}(2,-1)$, axes along the coordinate directions. Right: the symmetric $\\left(\\begin{smallmatrix}2&1\\\\1&2\\end{smallmatrix}\\right)$, axes along the diagonal directions.](../img/mdl-la-eig-ellipse.svg)\n",
    "\n",
    "The green arrows are the unit eigenvectors; the red arrows are their images\n",
    "$\\mathbf{A}\\mathbf{v}_i=\\lambda_i\\mathbf{v}_i$, which fall exactly on the\n",
    "ellipse's axes. For $\\operatorname{diag}(2,-1)$ the axes are the coordinate\n",
    "directions; for the symmetric $[[2,1],[1,2]]$ they are the diagonal directions\n",
    "$[1,1]^\\top/\\sqrt2$ (scaled by $3$) and $[1,-1]^\\top/\\sqrt2$ (scaled by $1$).\n",
    "\n",
    "### Finding Eigenvalues\n",
    "Let's figure out how to find them. By subtracting off the $\\lambda \\mathbf{v}$ from both sides,\n",
    "and then factoring out the vector,\n",
    "we see the above is equivalent to:\n",
    "\n",
    "$$(\\mathbf{A} - \\lambda \\mathbf{I})\\mathbf{v} = 0.$$\n",
    "\n",
    "For the equation to happen, we see that $(\\mathbf{A} - \\lambda \\mathbf{I})$\n",
    "must compress some direction down to zero,\n",
    "hence it is not invertible, and thus the determinant is zero.\n",
    "Thus, we can find the *eigenvalues*\n",
    "by finding for what $\\lambda$ is $\\det(\\mathbf{A}-\\lambda \\mathbf{I}) = 0$.\n",
    "Once we find the eigenvalues, we can solve\n",
    "$\\mathbf{A}\\mathbf{v} = \\lambda \\mathbf{v}$\n",
    "to find the associated *eigenvector(s)*.\n",
    "\n",
    "Let's see this with a more challenging matrix\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "2 & 1\\\\\n",
    "2 & 3\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "If we consider $\\det(\\mathbf{A}-\\lambda \\mathbf{I}) = 0$,\n",
    "we see this is equivalent to a polynomial equation in $\\lambda$: the\n",
    "*characteristic polynomial* $p(\\lambda) = \\det(\\mathbf{A}-\\lambda\\mathbf{I})$,\n",
    "here $(2-\\lambda)(3-\\lambda)-2 = (4-\\lambda)(1-\\lambda)$.\n",
    "Thus the two eigenvalues are $\\lambda = 1$ and $\\lambda = 4$.\n",
    "To find the associated eigenvectors, we solve\n",
    "$(\\mathbf{A}-\\mathbf{I})\\mathbf{v}=\\mathbf 0$ for $\\lambda = 1$ and\n",
    "$(\\mathbf{A}-4\\mathbf{I})\\mathbf{v}=\\mathbf 0$ for $\\lambda = 4$, i.e.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "2 & 1\\\\\n",
    "2 & 3\n",
    "\\end{bmatrix}\\begin{bmatrix}x \\\\ y\\end{bmatrix} = \\begin{bmatrix}x \\\\ y\\end{bmatrix}  \\; \\textrm{and} \\;\n",
    "\\begin{bmatrix}\n",
    "2 & 1\\\\\n",
    "2 & 3\n",
    "\\end{bmatrix}\\begin{bmatrix}x \\\\ y\\end{bmatrix}  = \\begin{bmatrix}4x \\\\ 4y\\end{bmatrix} .\n",
    "$$\n",
    "\n",
    "These are solved (up to scale) by $[1, -1]^\\top$ for $\\lambda = 1$\n",
    "and $[1, 2]^\\top$ for $\\lambda = 4$, respectively.\n",
    "(Each eigenvector is determined only up to a nonzero scalar multiple: if\n",
    "$\\mathbf v$ is an eigenvector, so is $c\\mathbf v$ for any $c\\neq0$, since both\n",
    "sides of the equation scale by $c$.)\n",
    "\n",
    "We can check this in code using the built-in `eig` routine."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-an-example"
   },
   "source": [
    "import numpy as onp\n",
    "onp.linalg.eig(onp.array([[2, 1], [2, 3]]))"
   ],
   "execution_count": 2,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "EigResult(eigenvalues=array([1., 4.]), eigenvectors=array([[-0.70710678, -0.4472136 ],\n       [ 0.70710678, -0.89442719]]))"
     },
     "metadata": {},
     "execution_count": 2
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-3 fw=mxnet -->\n",
    "\n",
    "MXNet arrays have no complex dtype, so `mxnet.np.linalg.eig` fails whenever the\n",
    "spectrum is complex, as it is for some matrices later in this section; we\n",
    "therefore use plain NumPy for eigenvalue computations here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-4 fw=all -->\n",
    "\n",
    "Note that the library normalizes the eigenvectors to be of length one,\n",
    "whereas we took ours to be of arbitrary length.\n",
    "Additionally, the choice of sign is arbitrary.\n",
    "However, the vectors computed are parallel\n",
    "to the ones we found by hand with the same eigenvalues.\n",
    "\n",
    "What did `eig` actually do? For $n\\ge5$ there is no\n",
    "formula in radicals for the roots of the characteristic polynomial, and\n",
    "practical libraries never form that polynomial anyway. They instead reduce\n",
    "$\\mathbf{A}$ to a nearly-triangular *Hessenberg* form and run the shifted *QR\n",
    "algorithm*, an iteration that drives the matrix toward triangular form while\n",
    "preserving its eigenvalues [@Golub.Van-Loan.1996]. The iteration targets\n",
    "triangular form because a triangular matrix displays its eigenvalues on its\n",
    "diagonal: $\\mathbf{T}-\\lambda\\mathbf{I}$ is again triangular, so by the\n",
    "triangular determinant rule of that section,\n",
    "$\\det(\\mathbf{T}-\\lambda\\mathbf{I})=\\prod_i(t_{ii}-\\lambda)$, which vanishes\n",
    "exactly at the diagonal values. The *power iteration* we\n",
    "develop later in this section matters in the complementary regime, when\n",
    "$\\mathbf{A}$ is too large to factor and only matrix--vector products are\n",
    "affordable.\n",
    "\n",
    "### Eigendecomposition and What It Computes\n",
    "\n",
    "Let's continue the previous example one step further.  Let\n",
    "\n",
    "$$\n",
    "\\mathbf{W} = \\begin{bmatrix}\n",
    "1 & 1 \\\\\n",
    "-1 & 2\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "be the matrix where the columns are the eigenvectors of the matrix $\\mathbf{A}$. Let\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\Lambda} = \\begin{bmatrix}\n",
    "1 & 0 \\\\\n",
    "0 & 4\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "be the matrix with the associated eigenvalues on the diagonal.\n",
    "Then the definition of eigenvalues and eigenvectors tells us that\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\mathbf{W} =\\mathbf{W} \\boldsymbol{\\Lambda} .\n",
    "$$\n",
    "\n",
    "The matrix $\\mathbf{W}$ is invertible, so we may multiply both sides by $\\mathbf{W}^{-1}$ on the right,\n",
    "we see that we may write\n",
    "\n",
    "$$\\mathbf{A} = \\mathbf{W} \\boldsymbol{\\Lambda} \\mathbf{W}^{-1}.$$\n",
    "\n",
    "Below we will see some nice consequences of this,\n",
    "but for now we need only know that such a decomposition\n",
    "will exist as long as we can find a full collection\n",
    "of linearly independent eigenvectors (so that $\\mathbf{W}$ is invertible).\n",
    "A matrix that admits such a decomposition is called *diagonalizable*.\n",
    "\n",
    "The eigendecomposition the equation turns many matrix operations\n",
    "into scalar ones. As a first example, consider:\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^n = \\overbrace{\\mathbf{A}\\cdots \\mathbf{A}}^{\\textrm{$n$ times}} = \\overbrace{(\\mathbf{W}\\boldsymbol{\\Lambda} \\mathbf{W}^{-1})\\cdots(\\mathbf{W}\\boldsymbol{\\Lambda} \\mathbf{W}^{-1})}^{\\textrm{$n$ times}} =  \\mathbf{W}\\overbrace{\\boldsymbol{\\Lambda}\\cdots\\boldsymbol{\\Lambda}}^{\\textrm{$n$ times}}\\mathbf{W}^{-1} = \\mathbf{W}\\boldsymbol{\\Lambda}^n \\mathbf{W}^{-1}.\n",
    "$$\n",
    "\n",
    "This tells us that for any positive power of a matrix,\n",
    "the eigendecomposition is obtained by just raising the eigenvalues to the same power.\n",
    "The same can be shown for negative powers,\n",
    "so if we want to invert a matrix we need only consider\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^{-1} = \\mathbf{W}\\boldsymbol{\\Lambda}^{-1} \\mathbf{W}^{-1},\n",
    "$$\n",
    "\n",
    "or in other words, just invert each eigenvalue.\n",
    "This will work as long as each eigenvalue is non-zero,\n",
    "so we see that invertible is the same as having no zero eigenvalues.\n",
    "The same recipe defines the *matrix exponential*\n",
    "$e^{\\mathbf{A}t}=\\mathbf{W}e^{\\boldsymbol{\\Lambda}t}\\mathbf{W}^{-1}$, which\n",
    "solves the linear differential equation $\\dot{\\mathbf{x}}=\\mathbf{A}\\mathbf{x}$;\n",
    "there the *real parts* of the eigenvalues decide stability, as we show in\n",
    "that section.\n",
    "\n",
    "**Determinant and trace from the eigenvalues.** Two quantities we met in\n",
    "that section read off the spectrum directly.\n",
    "For a diagonalizable matrix, both follow in two lines from identities we have\n",
    "already proved. The determinant is multiplicative\n",
    "the equation, and applying multiplicativity to\n",
    "$\\mathbf{W}\\mathbf{W}^{-1}=\\mathbf{I}$ gives\n",
    "$\\det(\\mathbf{W}^{-1})=1/\\det(\\mathbf{W})$, so\n",
    "\n",
    "$$\n",
    "\\det(\\mathbf{A}) = \\det(\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^{-1})\n",
    "  = \\det(\\mathbf{W})\\det(\\boldsymbol{\\Lambda})\\det(\\mathbf{W}^{-1})\n",
    "  = \\det(\\boldsymbol{\\Lambda})\n",
    "  = \\lambda_1 \\lambda_2 \\cdots \\lambda_n,\n",
    "$$\n",
    "\n",
    "the product of all the eigenvalues. This matches the geometric picture: whatever\n",
    "stretching $\\mathbf{W}$ does, $\\mathbf{W}^{-1}$ undoes, so the only net volume\n",
    "change is by the diagonal $\\boldsymbol{\\Lambda}$, which scales volumes by the\n",
    "product of its entries. The trace is cyclic, $\\operatorname{tr}(\\mathbf{A}\\mathbf{B})\n",
    "=\\operatorname{tr}(\\mathbf{B}\\mathbf{A})$ the equation, so the\n",
    "same conjugation collapses:\n",
    "\n",
    "$$\n",
    "\\operatorname{tr}(\\mathbf{A})\n",
    "  = \\operatorname{tr}(\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^{-1})\n",
    "  = \\operatorname{tr}(\\boldsymbol{\\Lambda}\\mathbf{W}^{-1}\\mathbf{W})\n",
    "  = \\operatorname{tr}(\\boldsymbol{\\Lambda})\n",
    "  = \\lambda_1+\\cdots+\\lambda_n,\n",
    "$$\n",
    "\n",
    "the *sum* of the eigenvalues. So determinant and trace are the product and sum\n",
    "of the spectrum; both are *similarity invariants* (unchanged by\n",
    "$\\mathbf{A}\\mapsto\\mathbf{W}\\mathbf{A}\\mathbf{W}^{-1}$). The trace identity\n",
    "returns when we estimate $\\log\\det$ Jacobians with the Hutchinson trace\n",
    "estimator in that section.\n",
    "\n",
    "Both identities in fact hold for *every* square matrix, diagonalizable or not,\n",
    "provided the eigenvalues are counted with algebraic multiplicity over the\n",
    "complex numbers. The characteristic polynomial\n",
    "$p(\\lambda)=\\det(\\mathbf{A}-\\lambda\\mathbf{I})$ has degree $n$, and the\n",
    "*fundamental theorem of algebra*, which we borrow here as an unproved fact,\n",
    "says that it factors completely over $\\mathbb{C}$:\n",
    "$p(\\lambda)=\\prod_{i=1}^{n}(\\lambda_i-\\lambda)$. Setting $\\lambda=0$ recovers\n",
    "$\\det(\\mathbf{A})=\\prod_i\\lambda_i$, and comparing the coefficients of\n",
    "$\\lambda^{n-1}$ on both sides recovers\n",
    "$\\operatorname{tr}(\\mathbf{A})=\\sum_i\\lambda_i$: on the right the coefficient is\n",
    "$(-1)^{n-1}\\sum_i\\lambda_i$, on the left it is $(-1)^{n-1}\\sum_i a_{ii}$ (this\n",
    "uses the permutation expansion of that section; only\n",
    "the all-diagonal term contributes at that degree). Complex eigenvalues of a real\n",
    "matrix enter in conjugate pairs, so both the product and the sum come out real.\n",
    "We meet these complex pairs concretely in that section.\n",
    "\n",
    "Finally, recall that the rank was the maximum number of linearly independent\n",
    "columns of a matrix. For a *diagonalizable* matrix, the rank equals the number\n",
    "of non-zero eigenvalues: multiplying by the invertible $\\mathbf{W}$ and\n",
    "$\\mathbf{W}^{-1}$ changes no dimensions, so\n",
    "$\\operatorname{rank}\\mathbf{A}=\\operatorname{rank}\\boldsymbol{\\Lambda}$, the\n",
    "number of non-zero diagonal entries. (For a general matrix, the right invariant\n",
    "is the number of non-zero *singular values*; that section.)\n",
    "\n",
    "#### When Does an Eigenbasis Exist? Multiplicity and Diagonalizability\n",
    "\n",
    "Whether we can build an invertible $\\mathbf{W}$ comes down to *counting*\n",
    "eigenvectors, and the precise bookkeeping uses two notions of multiplicity,\n",
    "both read off the characteristic polynomial of\n",
    "that section. The number of times $\\lambda$\n",
    "appears as a root of $p(\\lambda)$ is its *algebraic multiplicity*. The dimension\n",
    "of its eigenspace (the number of independent eigenvectors it contributes) is its\n",
    "*geometric multiplicity*. These always satisfy\n",
    "\n",
    "$$\n",
    "1 \\le \\textrm{geometric mult.}(\\lambda) \\le \\textrm{algebraic mult.}(\\lambda),\n",
    "$$\n",
    "\n",
    "and they need not be equal. The matrix is diagonalizable precisely when its\n",
    "eigenvectors span $\\mathbb{R}^n$, which happens **if and only if geometric\n",
    "multiplicity equals algebraic multiplicity for every eigenvalue**. One\n",
    "sufficient condition covers the generic case.\n",
    "\n",
    "**Proposition (distinct eigenvalues $\\Rightarrow$ diagonalizable).**\n",
    "*Eigenvectors belonging to distinct eigenvalues are linearly independent. In\n",
    "particular, an $n\\times n$ matrix with $n$ distinct eigenvalues is\n",
    "diagonalizable.*\n",
    "\n",
    "**Proof.** Suppose, for contradiction, that some eigenvectors\n",
    "$\\mathbf{w}_1,\\ldots,\\mathbf{w}_m$ with distinct eigenvalues\n",
    "$\\lambda_1,\\ldots,\\lambda_m$ are dependent, and take a dependence relation\n",
    "$\\sum_{i=1}^{m} c_i\\mathbf{w}_i = \\mathbf 0$ with the *fewest* nonzero\n",
    "coefficients; relabel so $c_m\\neq0$. Any such relation has at least *two*\n",
    "nonzero coefficients: a single-term relation $c_m\\mathbf{w}_m=\\mathbf 0$ would\n",
    "force $\\mathbf{w}_m=\\mathbf 0$, and eigenvectors are nonzero by definition.\n",
    "Apply $\\mathbf{A}$ and subtract $\\lambda_m$\n",
    "times the relation:\n",
    "\n",
    "$$\n",
    "\\mathbf 0 = \\mathbf{A}\\!\\sum_i c_i\\mathbf{w}_i - \\lambda_m\\!\\sum_i c_i\\mathbf{w}_i\n",
    "   = \\sum_{i=1}^{m-1} c_i(\\lambda_i-\\lambda_m)\\,\\mathbf{w}_i .\n",
    "$$\n",
    "\n",
    "The $\\mathbf{w}_m$ term cancels, so we obtain a *shorter* relation. It is still\n",
    "nontrivial: the original relation had a second nonzero coefficient $c_i$ with\n",
    "$i<m$, and it survives as $c_i(\\lambda_i-\\lambda_m)\\neq0$ because the\n",
    "eigenvalues are distinct, contradicting minimality. Hence the eigenvectors are independent.\n",
    "With $n$ distinct eigenvalues we get $n$ independent eigenvectors, which form a\n",
    "basis of $\\mathbb{R}^n$, so $\\mathbf{W}$ is invertible. $\\blacksquare$\n",
    "\n",
    "When eigenvalues are *repeated*, diagonalizability can fail: an eigenvalue's\n",
    "eigenspace may be too small to fill out its algebraic multiplicity. The\n",
    "simplest example is the *shear* matrix\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "1 & 1 \\\\\n",
    "0 & 1\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Its characteristic polynomial is\n",
    "$\\det(\\mathbf{A}-\\lambda\\mathbf{I})=(1-\\lambda)^2$, so $\\lambda=1$ has algebraic\n",
    "multiplicity $2$. But solving $(\\mathbf{A}-\\mathbf{I})\\mathbf{v}=\\mathbf 0$ gives\n",
    "$\\bigl[\\begin{smallmatrix}0&1\\\\0&0\\end{smallmatrix}\\bigr]\\mathbf v=\\mathbf 0$,\n",
    "i.e. $v_2=0$: the eigenspace is the *one-dimensional* line spanned by\n",
    "$[1,0]^\\top$. The geometric multiplicity is $1$ while the algebraic multiplicity\n",
    "is $2$, so there is no eigenbasis and $\\mathbf{A}$ is not diagonalizable; a\n",
    "matrix with this shortfall is called *defective*.\n",
    "the figure makes the shortage visible: the shear\n",
    "slides every horizontal layer of the grid sideways, and the horizontal axis is\n",
    "the *only* line through the origin left pointing where it started. A\n",
    "diagonalizable $2\\times2$ matrix would exhibit two independent such lines; the\n",
    "shear has just one, and no change of basis can conjure a second.\n",
    "\n",
    "![The defective shear acting on the plane. Every horizontal line slides rigidly to the right by an amount proportional to its height, so the grid tilts, but the $x$-axis (green) is carried to itself: $\\mathbf{A}\\mathbf{e}_1=\\mathbf{e}_1$, the sole eigendirection, with $\\lambda=1$. The other basis vector picks up a horizontal component, $\\mathbf{A}\\mathbf{e}_2=(1,1)^\\top$, and *no* second direction is preserved: geometric multiplicity $1$ against algebraic multiplicity $2$.](../img/mdl-la-defective-shear.svg)\n",
    "\n",
    "#### Beyond Diagonalization: The Jordan Normal Form\n",
    "\n",
    "When the geometric multiplicity falls short of the algebraic multiplicity, no\n",
    "eigenbasis exists. What structure remains? The answer is completely catalogued.\n",
    "Over the complex numbers, every square matrix is similar to a block-diagonal\n",
    "matrix whose diagonal blocks are *Jordan blocks*\n",
    "\n",
    "$$\n",
    "\\mathbf{J}_k(\\lambda) =\n",
    "\\begin{bmatrix}\n",
    "\\lambda & 1      &        &        \\\\\n",
    "        & \\lambda& \\ddots &        \\\\\n",
    "        &        & \\ddots & 1      \\\\\n",
    "        &        &        & \\lambda\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "$k\\times k$ blocks with a single eigenvalue $\\lambda$ repeated down the\n",
    "diagonal and $1$s just above it. This is the *Jordan normal form*\n",
    "[@Horn.Johnson.2012], which we use without proof.\n",
    "\n",
    "The form makes both multiplicities visible at a glance, because each block\n",
    "contributes exactly *one* eigenvector: within a block,\n",
    "$(\\mathbf{J}_k(\\lambda)-\\lambda\\mathbf{I})\\mathbf{v}=\\mathbf 0$ forces every\n",
    "coordinate of $\\mathbf{v}$ except the first to vanish. So the geometric\n",
    "multiplicity of $\\lambda$ is the *number* of blocks carrying $\\lambda$, the\n",
    "algebraic multiplicity is the *total size* of those blocks, and a matrix is\n",
    "diagonalizable exactly when every block is $1\\times1$. The shear\n",
    "the equation is already in Jordan form: it is the single\n",
    "block $\\mathbf{J}_2(1)$, one block of size two, hence one eigenvector for an\n",
    "eigenvalue of algebraic multiplicity two, exactly the shortage that\n",
    "the figure displays.\n",
    "\n",
    "What does a Jordan block do under iteration? Write\n",
    "$\\mathbf{J}_k(\\lambda)=\\lambda\\mathbf{I}+\\mathbf{N}$, where $\\mathbf{N}$ has\n",
    "$1$s just above the diagonal and zeros elsewhere. $\\mathbf{N}$ is a shift: it\n",
    "sends each basis vector to its predecessor, so $\\mathbf{N}^k=\\mathbf 0$ (a\n",
    "matrix some power of which vanishes is called *nilpotent*). Because\n",
    "$\\mathbf{I}$ and $\\mathbf{N}$ commute, the binomial theorem applies to\n",
    "$(\\lambda\\mathbf{I}+\\mathbf{N})^n$, and for $k=2$ every term beyond the first\n",
    "two contains $\\mathbf{N}^2=\\mathbf 0$:\n",
    "\n",
    "$$\n",
    "\\mathbf{J}_2(\\lambda)^n = (\\lambda\\mathbf{I}+\\mathbf{N})^n\n",
    "  = \\lambda^n\\mathbf{I} + n\\lambda^{n-1}\\mathbf{N}.\n",
    "$$\n",
    "\n",
    "The eigenvalue still decides growth versus decay, but multiplied by a\n",
    "*polynomial* in $n$ rather than acting alone. A defective matrix with\n",
    "$|\\lambda|=1$, like the shear, drifts linearly instead of staying put:\n",
    "$\\mathbf{J}_2(1)^n$ has an $n$ in its corner. And with $|\\lambda|$ slightly\n",
    "below $1$, the off-diagonal entry $n\\lambda^{n-1}$ first *grows*, peaks near\n",
    "$n\\approx-1/\\log\\lambda$, and only then decays: the iterates rise transiently\n",
    "before the geometric factor wins. We can watch this happen for\n",
    "$\\lambda=0.99$, where the formula puts the peak near $n\\approx100$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-eigendecomposition-jordan-power"
   },
   "source": [
    "J = np.array([[0.99, 1.], [0., 0.99]])   # the Jordan block J_2(0.99)\n",
    "for n in [1, 50, 100, 200, 800]:\n",
    "    Jn = np.linalg.matrix_power(J, n)\n",
    "    pred = n * 0.99 ** (n - 1)   # predicted corner entry, n * lambda^(n-1)\n",
    "    print(f'n = {n:3d}   ||J^n|| = {np.linalg.norm(Jn, 2):8.3f}   '\n",
    "          f'n * 0.99^(n-1) = {pred:8.3f}')"
   ],
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "n =   1   ||J^n|| =    1.609   n * 0.99^(n-1) =    1.000\nn =  50   ||J^n|| =   30.568   n * 0.99^(n-1) =   30.556\nn = 100   ||J^n|| =   36.977   n * 0.99^(n-1) =   36.973\nn = 200   ||J^n|| =   27.067   n * 0.99^(n-1) =   27.067\nn = 800   ||J^n|| =    0.260   n * 0.99^(n-1) =    0.260\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-5 fw=all -->\n",
    "\n",
    "Every eigenvalue of $\\mathbf{J}$ is $0.99$, safely below $1$, yet the norm\n",
    "climbs to about $37$ before the decay takes over. This *transient growth* is a\n",
    "signature of defective (more generally, non-normal) matrices, and it is\n",
    "invisible to the eigenvalues alone; the singular values of\n",
    "that section measure it directly.\n",
    "\n",
    "Jordan blocks are the first of two ways a real matrix can fall short of a real\n",
    "eigenbasis: too few eigenvectors. The second way is different in kind: a real\n",
    "matrix can have no real eigenvalue at all in some plane, because it *rotates*\n",
    "that plane.\n",
    "\n",
    "#### Complex Eigenvalues Are Rotations\n",
    "\n",
    "Consider the planar rotation by angle $\\theta$,\n",
    "\n",
    "$$\n",
    "\\mathbf{R} = \\begin{bmatrix}\\cos\\theta & -\\sin\\theta\\\\\n",
    "                            \\sin\\theta & \\phantom{-}\\cos\\theta\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "If $0<\\theta<\\pi$, $\\mathbf{R}$ sends *no* nonzero real vector to a multiple of\n",
    "itself (a true rotation has no fixed direction), so it has no real eigenvector.\n",
    "Its characteristic polynomial $\\lambda^2-2\\cos\\theta\\,\\lambda+1$ has the\n",
    "complex-conjugate roots $\\lambda=\\cos\\theta\\pm i\\sin\\theta=e^{\\pm i\\theta}$. The\n",
    "*modulus* $|\\lambda|=1$ records that $\\mathbf{R}$ preserves length (it is\n",
    "orthogonal), while the *argument* $\\pm\\theta$ records the rotation angle."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-complex-rotation"
   },
   "source": [
    "theta = np.pi / 6\n",
    "R = np.array([[np.cos(theta), -np.sin(theta)],\n",
    "              [np.sin(theta),  np.cos(theta)]])\n",
    "print('eigenvalues of R:', np.round(np.linalg.eigvals(R), 4))\n",
    "print('expected e^{+/- i*theta}:', np.round(np.exp([1j * theta, -1j * theta]), 4))\n",
    "print('moduli (should be 1):', np.round(np.abs(np.linalg.eigvals(R)), 4))"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "eigenvalues of R: [0.866+0.5j 0.866-0.5j]\nexpected e^{+/- i*theta}: [0.866+0.5j 0.866-0.5j]\nmoduli (should be 1): [1. 1.]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-6 fw=all -->\n",
    "\n",
    "This is the general rule. Complex eigenvalues of a real matrix come in\n",
    "conjugate pairs $a\\pm ib=re^{\\pm i\\theta}$, and each pair encodes a *scale by\n",
    "$r$ and rotate by $\\theta$* on some two-dimensional invariant plane. But where\n",
    "is that plane, and how do we pass between the complex pair and the real\n",
    "rotation? The dictionary is a two-line computation. Let\n",
    "$\\mathbf{A}\\mathbf{w}=\\lambda\\mathbf{w}$ with $\\lambda=a+ib$, $b\\neq0$, and\n",
    "split the eigenvector into real and imaginary parts,\n",
    "$\\mathbf{w}=\\mathbf{u}+i\\mathbf{v}$. Expanding\n",
    "$\\mathbf{A}(\\mathbf{u}+i\\mathbf{v})=(a+ib)(\\mathbf{u}+i\\mathbf{v})$ and\n",
    "matching real and imaginary parts gives\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\mathbf{u} = a\\mathbf{u} - b\\mathbf{v},\n",
    "\\qquad\n",
    "\\mathbf{A}\\mathbf{v} = b\\mathbf{u} + a\\mathbf{v}.\n",
    "$$\n",
    "\n",
    "Both images land back in $\\textrm{span}\\{\\mathbf{u},\\mathbf{v}\\}$, so this span\n",
    "is the invariant plane. (It really is a plane: if $\\mathbf{v}$ were a multiple\n",
    "of $\\mathbf{u}$, then $\\mathbf{u}$ would be a *real* eigenvector of\n",
    "$\\mathbf{A}$, forcing its eigenvalue to be real and contradicting $b\\neq0$.)\n",
    "Now collect the two identities into columns. With\n",
    "$\\mathbf{P}=[\\mathbf{v}\\;\\;\\mathbf{u}]$, they say exactly that\n",
    "$\\mathbf{A}\\mathbf{P}=\\mathbf{P}\\left(\\begin{smallmatrix}a & -b\\\\ b & \\phantom{-}a\\end{smallmatrix}\\right)$,\n",
    "i.e.\n",
    "\n",
    "$$\n",
    "\\mathbf{P}^{-1}\\mathbf{A}\\mathbf{P}\n",
    "= \\begin{bmatrix} a & -b\\\\ b & \\phantom{-}a \\end{bmatrix}\n",
    "= \\underbrace{\\sqrt{a^2+b^2}}_{r}\\,\n",
    "  \\begin{bmatrix} \\cos\\theta & -\\sin\\theta\\\\ \\sin\\theta & \\phantom{-}\\cos\\theta \\end{bmatrix},\n",
    "\\qquad \\theta=\\arg(a+ib).\n",
    "$$\n",
    "\n",
    "In the basis $(\\mathbf{v},\\mathbf{u})$ of its invariant plane, $\\mathbf{A}$\n",
    "*is* the rotation--scaling block: scale by the modulus, rotate by the\n",
    "argument. The dictionary reads just as easily right to left. Identify\n",
    "$(x,y)^\\top\\in\\mathbb{R}^2$ with the complex number $x+iy$; then the block acts\n",
    "as $(x,y)^\\top\\mapsto(ax-by,\\,bx+ay)^\\top$, which is precisely multiplication\n",
    "by $a+ib$, and its eigenvalues are $a\\pm ib$ with eigenvectors\n",
    "$(1,\\mp i)^\\top$. Complex pair and real block carry the same two numbers in\n",
    "different coordinates, polar $(r,\\theta)$ versus Cartesian $(a,b)$. Two\n",
    "conventions to keep straight: choosing the conjugate eigenvalue $a-ib$ instead\n",
    "flips the sign of $b$, giving the same plane traversed in the opposite\n",
    "orientation; and replacing $\\mathbf{w}$ by any nonzero complex multiple changes\n",
    "$\\mathbf{P}$ but not the block, since the derivation used only\n",
    "$\\mathbf{A}\\mathbf{w}=\\lambda\\mathbf{w}$. So the block is well defined even\n",
    "though `eig` normalizes eigenvectors with an arbitrary phase. Let's verify the\n",
    "dictionary on a matrix with no visible rotation in it."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-eigendecomposition-complex-eigenvalues-are-rotations-1"
   },
   "source": [
    "A = np.array([[0.5, -1.2], [0.9, 0.3]])\n",
    "lam, W = np.linalg.eig(A)\n",
    "k = np.argmax(lam.imag)                  # the eigenvalue a + ib with b > 0\n",
    "r, theta = np.abs(lam[k]), np.angle(lam[k])\n",
    "u, v = W[:, k].real, W[:, k].imag\n",
    "P = np.column_stack([v, u])              # basis (v, u) of the invariant plane\n",
    "print(f'lambda = {lam[k]:.4f}   r = {r:.4f}   theta = {theta:.4f}')\n",
    "print('P^-1 A P =\\n', np.round(np.linalg.inv(P) @ A @ P, 4))\n",
    "print('r * R(theta) =\\n', np.round(r * np.array(\n",
    "    [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]), 4))"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "lambda = 0.4000+1.0344j   r = 1.1091   theta = 1.2018\nP^-1 A P =\n [[ 0.4    -1.0344]\n [ 1.0344  0.4   ]]\nr * R(theta) =\n [[ 0.4    -1.0344]\n [ 1.0344  0.4   ]]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-7 fw=all -->\n",
    "\n",
    "Nothing about the entries of $\\mathbf{A}$ advertises a rotation, yet in the\n",
    "basis read off from its eigenvector the map is exactly scale-by-$r$,\n",
    "rotate-by-$\\theta$. Allowing such rotation--scaling blocks alongside\n",
    "the Jordan blocks of that section gives the *real Jordan form*\n",
    "[@Horn.Johnson.2012], and with it a complete inventory of what a linear\n",
    "map can do over the reals: every square matrix decomposes into pure stretch\n",
    "directions (real eigenvalues with a full set of eigenvectors), shear-like\n",
    "Jordan blocks (repeated eigenvalues with too few eigenvectors), and\n",
    "rotation--scaling planes (complex-conjugate pairs). The inventory also predicts\n",
    "how iterating a matrix can misbehave: a dominant rotating pair keeps the\n",
    "iterate spinning forever instead of settling on a direction, and a dominant\n",
    "defective block converges with a polynomial drag. Both reappear when we develop\n",
    "*power iteration* later in this section. The dictionary itself returns in\n",
    "that section, where it converts eigenvalues into the\n",
    "phase portraits of continuous-time flows $\\dot{\\mathbf{x}}=\\mathbf{A}\\mathbf{x}$:\n",
    "a complex pair $a\\pm ib$ produces a spiral, with the real part $a$ setting\n",
    "growth or decay and the imaginary part $b$ the angular velocity.\n",
    "\n",
    "## Non-Normal Matrices and Transient Amplification\n",
    "\n",
    "The Jordan example showed that eigenvalues can miss a large temporary excursion.\n",
    "Defectiveness is the sharpest version of the problem, but it is not the whole\n",
    "problem: a matrix can have distinct eigenvalues, be diagonalizable, and still\n",
    "amplify some inputs dramatically before its asymptotic decay becomes visible.\n",
    "The property separating the clean case from this behavior is **normality**.\n",
    "\n",
    "### Normality: When Eigenvalues Control Norms\n",
    "\n",
    "Over the complex numbers, a matrix is **normal** when it commutes with its\n",
    "conjugate transpose,\n",
    "\n",
    "$$\n",
    "\\mathbf A^*\\mathbf A=\\mathbf A\\mathbf A^*.\n",
    "$$\n",
    "\n",
    "For a real matrix, $\\mathbf A^*$ is simply $\\mathbf A^\\top$. Symmetric,\n",
    "orthogonal, and unitary matrices are normal, as are planar rotation--scaling\n",
    "blocks. The complex spectral theorem says that a matrix is normal if and only\n",
    "if it has an orthonormal eigenbasis over $\\mathbb C$:\n",
    "$\\mathbf A=\\mathbf U\\boldsymbol\\Lambda\\mathbf U^*$ with $\\mathbf U$ unitary.\n",
    "Consequently,\n",
    "\n",
    "$$\n",
    "\\|\\mathbf A^k\\|_2\n",
    "=\\|\\mathbf U\\boldsymbol\\Lambda^k\\mathbf U^*\\|_2\n",
    "=\\max_i|\\lambda_i|^k\n",
    "=\\rho(\\mathbf A)^k.\n",
    "$$\n",
    "\n",
    "For normal matrices the eigenvalues therefore control the largest possible\n",
    "amplification at *every* iteration, not merely its asymptotic exponential\n",
    "rate. This is why the orthonormal eigenbasis in the next section is more than\n",
    "an aesthetic convenience: the change of coordinates neither stretches nor\n",
    "nearly collapses any direction.\n",
    "\n",
    "A diagonalizable non-normal matrix instead has\n",
    "$\\mathbf A=\\mathbf W\\boldsymbol\\Lambda\\mathbf W^{-1}$ with a generally\n",
    "non-orthogonal eigenvector matrix. The elementary bound\n",
    "\n",
    "$$\n",
    "\\|\\mathbf A^k\\|_2\n",
    "\\le \\underbrace{\\|\\mathbf W\\|_2\\|\\mathbf W^{-1}\\|_2}_{\\kappa_2(\\mathbf W)}\n",
    "   \\rho(\\mathbf A)^k\n",
    "$$\n",
    "\n",
    "exposes the missing factor: nearly parallel eigenvectors make\n",
    "$\\kappa_2(\\mathbf W)$ large. Individual eigenmodes then carry large\n",
    "coefficients that can reinforce before their eventual cancellation. The bound\n",
    "need not be tight, but it explains why eigenvalues alone can be a poor\n",
    "finite-time diagnostic.\n",
    "\n",
    "### A Stable Matrix That First Grows\n",
    "\n",
    "Consider the upper-triangular matrix\n",
    "\n",
    "$$\n",
    "\\mathbf A=\\begin{bmatrix}0.9&4\\\\0&0.8\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Its distinct eigenvalues $0.9$ and $0.8$ make it diagonalizable, and both lie\n",
    "strictly inside the unit circle. Nevertheless, its eigenvectors are nearly\n",
    "parallel and its first few powers are large. The off-diagonal entry has the\n",
    "closed form\n",
    "\n",
    "$$\n",
    "(\\mathbf A^k)_{12}\n",
    "=4\\,\\frac{0.9^k-0.8^k}{0.9-0.8}\n",
    "=40(0.9^k-0.8^k),\n",
    "$$\n",
    "\n",
    "which initially grows even though both geometric factors decay. Singular\n",
    "values, rather than eigenvalue moduli, measure this one-step and finite-time\n",
    "stretch. The NumPy cell compares the non-normal matrix with the diagonal\n",
    "matrix having exactly the same eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-eigendecomposition-nonnormal-transient"
   },
   "source": [
    "A = np.array([[0.9, 4.0], [0.0, 0.8]])\n",
    "D = np.diag([0.9, 0.8])                         # same eigenvalues, normal\n",
    "ks = np.arange(0, 61)\n",
    "norm_A = np.array([np.linalg.norm(np.linalg.matrix_power(A, k), 2) for k in ks])\n",
    "norm_D = np.array([np.linalg.norm(np.linalg.matrix_power(D, k), 2) for k in ks])\n",
    "k_peak = int(norm_A.argmax())\n",
    "W = np.linalg.eig(A)[1]\n",
    "z = 1.0\n",
    "res_A = np.linalg.norm(np.linalg.inv(z * np.eye(2) - A), 2)\n",
    "res_D = np.linalg.norm(np.linalg.inv(z * np.eye(2) - D), 2)\n",
    "print('eigenvalues:', np.linalg.eigvals(A))\n",
    "print(f'normal matrix: max ||D^k|| = {norm_D.max():.3f}')\n",
    "print(f'non-normal:    max ||A^k|| = {norm_A[k_peak]:.3f} at k = {k_peak}')\n",
    "print(f'eigenvector condition number = {np.linalg.cond(W):.1f}')\n",
    "print(f'at z=1: ||(zI-A)^-1|| = {res_A:.1f},  normal comparison = {res_D:.1f}')"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "eigenvalues: [0.9 0.8]\nnormal matrix: max ||D^k|| = 1.000\nnon-normal:    max ||A^k|| = 10.788 at k = 6\neigenvector condition number = 80.0\nat z=1: ||(zI-A)^-1|| = 200.3,  normal comparison = 10.0\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-8 fw=all -->\n",
    "\n",
    "The normal matrix contracts from the first step. The non-normal matrix has the\n",
    "same spectral radius but amplifies a suitably chosen input by more than an\n",
    "order of magnitude before decay wins. This is **transient amplification**:\n",
    "stability as $k\\to\\infty$ coexisting with substantial finite-time growth.\n",
    "the figure draws the whole story: on a log scale\n",
    "the two norm curves share the same asymptotic slope — that is the common\n",
    "spectral radius — but only after the non-normal matrix has climbed through a\n",
    "ten-fold excursion that eigenvalues alone gave no hint of.\n",
    "\n",
    "![Two matrices with identical eigenvalues, radically different finite-time behavior. Left: the spectral norm of the powers of the non-normal matrix climbs by an order of magnitude before the asymptotic decay — the same log-slope as its normal twin — takes over. Right: the reason, seen through pseudospectra; perturbations of size $\\varepsilon$ can move the non-normal matrix's eigenvalues far outside the unit circle, while the normal matrix's stay in tight disks.](../img/mdl-la-transient-growth.svg)\n",
    "\n",
    "### Pseudospectra: Stability under Perturbation\n",
    "\n",
    "A second symptom of non-normality is eigenvalue sensitivity. For\n",
    "$\\varepsilon>0$, the **$\\varepsilon$-pseudospectrum** is\n",
    "\n",
    "$$\n",
    "\\Lambda_\\varepsilon(\\mathbf A)\n",
    "=\\{z\\in\\mathbb C:\\sigma_{\\min}(z\\mathbf I-\\mathbf A)\\le\\varepsilon\\}\n",
    "=\\{z:z\\text{ is an eigenvalue of }\\mathbf A+\\mathbf E\n",
    "       \\text{ for some }\\|\\mathbf E\\|_2\\le\\varepsilon\\}.\n",
    "$$\n",
    "\n",
    "Away from an eigenvalue this can equivalently be read through the\n",
    "**resolvent** $(z\\mathbf I-\\mathbf A)^{-1}$: membership means its norm is at\n",
    "least $1/\\varepsilon$. For a normal matrix the resolvent norm is exactly the\n",
    "reciprocal distance from $z$ to the nearest eigenvalue. For a non-normal\n",
    "matrix it can be much larger, so small perturbations can move eigenvalues far\n",
    "and a forced linear system can respond strongly even when $z$ is not close to\n",
    "the displayed spectrum. The final line of the cell measures this effect at\n",
    "$z=1$: the two matrices have the same eigenvalues, yet the non-normal\n",
    "resolvent is roughly twenty times larger. The right panel of\n",
    "the figure draws the comparison whole — around the\n",
    "normal matrix's eigenvalues the pseudospectral level sets are tight\n",
    "concentric disks, while around the non-normal matrix's they balloon far\n",
    "beyond the unit circle at the same $\\varepsilon$.\n",
    "\n",
    "Pseudospectra are primarily a diagnostic rather than another decomposition to\n",
    "memorize. In numerical linear algebra they warn that computed eigenvalues may\n",
    "be sensitive; in dynamical systems they warn about transient responses; and\n",
    "in deep learning they explain why a spectral-radius constraint alone does not\n",
    "control a layer's Lipschitz constant or its gradient amplification.\n",
    "\n",
    "### Products of Jacobians\n",
    "\n",
    "A recurrent or very deep network adds one complication: its backward map is a\n",
    "product of generally different Jacobians,\n",
    "\n",
    "$$\n",
    "\\mathbf J_T\\mathbf J_{T-1}\\cdots\\mathbf J_1.\n",
    "$$\n",
    "\n",
    "Even if every factor has eigenvalues inside the unit circle, the factors need\n",
    "not share eigenvectors and their expanding singular directions can rotate into\n",
    "one another. There is therefore no rule that obtains the product's norm by\n",
    "multiplying the factors' spectral radii. Singular values give finite-horizon\n",
    "bounds,\n",
    "$\\|\\mathbf J_T\\cdots\\mathbf J_1\\|_2\\le\\prod_t\\sigma_{\\max}(\\mathbf J_t)$;\n",
    "long-run average logarithmic growth is summarized by **Lyapunov exponents**.\n",
    "The practical hierarchy is now precise: eigenvalues describe asymptotic powers\n",
    "of one fixed map, singular values describe finite-time stretch, pseudospectra\n",
    "describe sensitivity, and Lyapunov exponents describe long products of\n",
    "changing maps.\n",
    "\n",
    "## Symmetric Matrices and Positive Definiteness\n",
    "\n",
    "### The Spectral Theorem\n",
    "\n",
    "Shears and rotations tell us what can go wrong. We now restrict attention to\n",
    "the family of matrices for which a full set of orthonormal eigenvectors is\n",
    "*guaranteed*. The most commonly encountered such family are the\n",
    "*symmetric matrices*, those with $\\mathbf{A}=\\mathbf{A}^\\top$. Their guarantee has\n",
    "a name.\n",
    "\n",
    "**Theorem (Spectral Theorem).** *Every real symmetric matrix $\\mathbf{A}$ has a\n",
    "full orthonormal eigenbasis with real eigenvalues. Equivalently,\n",
    "$\\mathbf{A}=\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^\\top$ where\n",
    "$\\boldsymbol{\\Lambda}=\\operatorname{diag}(\\lambda_1,\\ldots,\\lambda_n)$ is real and\n",
    "$\\mathbf{W}$ is orthogonal* ($\\mathbf{W}^\\top\\mathbf{W}=\\mathbf{I}$, so\n",
    "$\\mathbf{W}^{-1}=\\mathbf{W}^\\top$).\n",
    "\n",
    "We prove the three claims in turn (real eigenvalues, orthogonal eigenvectors,\n",
    "and a complete basis); each proof is short and instructive. The inner product\n",
    "of complex vectors here is $\\langle\\mathbf x,\\mathbf y\\rangle=\\mathbf x^*\\mathbf y$\n",
    "with $\\mathbf x^*$ the conjugate transpose; on real vectors it is the ordinary\n",
    "dot product. The one fact we borrow is that a real polynomial, here the\n",
    "characteristic polynomial, has at least one (possibly complex) root, a special\n",
    "case of the fundamental theorem of algebra we invoked above, so at least\n",
    "one eigenpair exists over $\\mathbb{C}$.\n",
    "\n",
    "**Proof, part (i): the eigenvalues are real.** Let\n",
    "$\\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}$ with $\\mathbf{v}\\neq\\mathbf 0$, allowing\n",
    "complex entries. Since $\\mathbf{A}$ is real and symmetric, $\\mathbf A^*=\\mathbf A^\\top=\\mathbf A$, so\n",
    "\n",
    "$$\n",
    "\\lambda\\,\\mathbf{v}^*\\mathbf{v}\n",
    "   = \\mathbf{v}^*\\mathbf{A}\\mathbf{v}\n",
    "   = (\\mathbf{A}\\mathbf{v})^*\\mathbf{v}\n",
    "   = \\bar\\lambda\\,\\mathbf{v}^*\\mathbf{v}.\n",
    "$$\n",
    "\n",
    "Because $\\mathbf{v}^*\\mathbf{v}=\\|\\mathbf v\\|^2>0$, we may divide it out to get\n",
    "$\\lambda=\\bar\\lambda$, i.e. $\\lambda$ is real. The eigenvector can then be taken\n",
    "real as well (it solves the real system $(\\mathbf A-\\lambda\\mathbf I)\\mathbf v=\\mathbf0$). $\\blacksquare$\n",
    "\n",
    "**Proof, part (ii): eigenvectors for distinct eigenvalues are orthogonal.** Let\n",
    "$\\mathbf{A}\\mathbf{u}=\\lambda\\mathbf{u}$ and\n",
    "$\\mathbf{A}\\mathbf{v}=\\mu\\mathbf{v}$ with $\\lambda\\neq\\mu$. Using\n",
    "$\\mathbf{A}=\\mathbf{A}^\\top$ to slide $\\mathbf{A}$ across the inner product,\n",
    "\n",
    "$$\n",
    "\\lambda\\langle\\mathbf{u},\\mathbf{v}\\rangle\n",
    "  = \\langle\\mathbf{A}\\mathbf{u},\\mathbf{v}\\rangle\n",
    "  = \\langle\\mathbf{u},\\mathbf{A}\\mathbf{v}\\rangle\n",
    "  = \\mu\\langle\\mathbf{u},\\mathbf{v}\\rangle,\n",
    "$$\n",
    "\n",
    "so $(\\lambda-\\mu)\\langle\\mathbf{u},\\mathbf{v}\\rangle=0$. Since\n",
    "$\\lambda\\neq\\mu$, we conclude $\\langle\\mathbf{u},\\mathbf{v}\\rangle=0$: the\n",
    "eigenvectors are orthogonal. $\\blacksquare$\n",
    "\n",
    "**Proof, part (iii): a full orthonormal basis (induction on the orthogonal complement).** We induct on\n",
    "$n$. For $n=1$ the claim is trivial. For $n>1$, part (i) gives a real eigenpair\n",
    "$(\\lambda_1,\\mathbf{w}_1)$ with $\\|\\mathbf{w}_1\\|=1$. Let\n",
    "$U=\\mathbf{w}_1^{\\perp}$ be its orthogonal complement (dimension $n-1$). This\n",
    "subspace is *$\\mathbf{A}$-invariant*: for any $\\mathbf{x}\\in U$,\n",
    "\n",
    "$$\n",
    "\\langle\\mathbf{A}\\mathbf{x},\\mathbf{w}_1\\rangle\n",
    "  = \\langle\\mathbf{x},\\mathbf{A}\\mathbf{w}_1\\rangle\n",
    "  = \\lambda_1\\langle\\mathbf{x},\\mathbf{w}_1\\rangle = 0,\n",
    "$$\n",
    "\n",
    "so $\\mathbf{A}\\mathbf{x}\\in U$ too. The restriction $\\mathbf{A}|_U$ is again\n",
    "symmetric on an $(n-1)$-dimensional space, so by the inductive hypothesis it has\n",
    "an orthonormal eigenbasis $\\mathbf{w}_2,\\ldots,\\mathbf{w}_n$ of $U$. Together with\n",
    "$\\mathbf{w}_1$ these form an orthonormal eigenbasis of $\\mathbb{R}^n$. Collecting\n",
    "them as the columns of $\\mathbf{W}$ gives the stated factorization. $\\blacksquare$\n",
    "\n",
    "In this special case, then, we can write the eigendecomposition\n",
    "the equation as\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^\\top .\n",
    "$$\n",
    "\n",
    "This factorization is what makes PCA (the variance-maximizing projections of\n",
    "that section), covariance analysis, and the Hessian test for minima\n",
    "(that section) computable: all three concern symmetric matrices. It\n",
    "also reads geometrically: $\\mathbf{W}^\\top$ rotates the orthonormal eigenvectors onto the\n",
    "coordinate axes, $\\boldsymbol{\\Lambda}$ stretches independently along each axis,\n",
    "and $\\mathbf{W}$ rotates back. This is precisely the \"circle becomes an ellipse\n",
    "with axes along the eigenvectors\" picture we drew above, now justified.\n",
    "\n",
    "The defective shear the equation shows that the\n",
    "eigendecomposition is *picky*: it needs a square matrix, and to be fully\n",
    "well-behaved it really wants a symmetric one. The fix that works for *every*\n",
    "matrix $\\mathbf{A}$ rests on a single observation we will use repeatedly: the\n",
    "matrix $\\mathbf{A}^\\top\\mathbf{A}$ is always **symmetric** ($(\\mathbf A^\\top\\mathbf\n",
    "A)^\\top=\\mathbf A^\\top\\mathbf A$) and **positive semidefinite**, because\n",
    "\n",
    "$$\n",
    "\\mathbf{x}^\\top(\\mathbf{A}^\\top\\mathbf{A})\\mathbf{x} = \\|\\mathbf{A}\\mathbf{x}\\|^2 \\ge 0 .\n",
    "$$\n",
    "\n",
    "The spectral theorem therefore applies to $\\mathbf{A}^\\top\\mathbf{A}$, and feeding\n",
    "its orthonormal eigenvectors and (non-negative) eigenvalues through $\\mathbf{A}$\n",
    "manufactures the singular value decomposition\n",
    "$\\mathbf{A}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$: in short, **the SVD is\n",
    "the eigendecomposition of $\\mathbf{A}^\\top\\mathbf{A}$ in disguise**, with singular\n",
    "values $\\sigma_i=\\sqrt{\\lambda_i(\\mathbf{A}^\\top\\mathbf{A})}$. We carry out the\n",
    "construction, and give the defective shear its SVD, in\n",
    "that section. Positive semidefiniteness, the property doing\n",
    "the work here, is what we take up next.\n",
    "\n",
    "### Positive (Semi)Definiteness\n",
    "\n",
    "A symmetric matrix's eigenvalues answer a *sign* question: is the quadratic form\n",
    "$\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}$ ever negative? This one property is what\n",
    "makes covariance matrices, Gram matrices, and \"is this a minimum?\" Hessian tests\n",
    "work.\n",
    "\n",
    "A symmetric matrix $\\mathbf{A}$ is *positive semidefinite* (PSD), written\n",
    "$\\mathbf{A}\\succeq0$, if\n",
    "\n",
    "$$\n",
    "\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x} \\ge 0 \\quad\\textrm{for all } \\mathbf{x},\n",
    "$$\n",
    "\n",
    "and *positive definite* (PD), written $\\mathbf{A}\\succ0$, if the inequality is\n",
    "strict for every $\\mathbf{x}\\neq\\mathbf 0$. The eigenvalues decide the question\n",
    "completely.\n",
    "\n",
    "**Proposition (PSD/PD via eigenvalues).** *For symmetric $\\mathbf{A}$,*\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\succeq0 \\iff \\lambda_i\\ge0\\ \\forall i,\n",
    "\\qquad\n",
    "\\mathbf{A}\\succ0 \\iff \\lambda_i>0\\ \\forall i .\n",
    "$$\n",
    "\n",
    "**Proof.** Use the spectral theorem\n",
    "$\\mathbf{A}=\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^\\top$ and the change of\n",
    "variables $\\mathbf{y}=\\mathbf{W}^\\top\\mathbf{x}$. Since $\\mathbf{W}$ is\n",
    "orthogonal this is just a rotation, and\n",
    "\n",
    "$$\n",
    "\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}\n",
    "  = \\mathbf{x}^\\top\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^\\top\\mathbf{x}\n",
    "  = \\mathbf{y}^\\top\\boldsymbol{\\Lambda}\\mathbf{y}\n",
    "  = \\sum_{i=1}^{n}\\lambda_i\\,y_i^2\n",
    "  = \\sum_{i=1}^{n}\\lambda_i\\,(\\mathbf{w}_i^\\top\\mathbf{x})^2 .\n",
    "$$\n",
    "\n",
    "The quadratic form is a *weighted sum of squares* with the eigenvalues as\n",
    "weights. If every $\\lambda_i\\ge0$ the sum is $\\ge0$, so $\\mathbf{A}\\succeq0$.\n",
    "Conversely, choosing $\\mathbf{x}=\\mathbf{w}_j$ (so $\\mathbf y=\\mathbf e_j$) gives\n",
    "$\\mathbf{w}_j^\\top\\mathbf{A}\\mathbf{w}_j=\\lambda_j$, which forces\n",
    "$\\lambda_j\\ge0$. The strict (PD) statement is identical with $\\ge$ replaced by\n",
    "$>$. $\\blacksquare$\n",
    "\n",
    "Three consequences follow immediately, and together they cover most of the places\n",
    "PSD matrices appear in deep learning.\n",
    "\n",
    "* **PD $\\Rightarrow$ invertible.** If $\\mathbf{A}\\succ0$ and\n",
    "  $\\mathbf{A}\\mathbf{x}=\\mathbf 0$, then\n",
    "  $\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}=0$, which forces $\\mathbf{x}=\\mathbf 0$; the\n",
    "  null space is trivial. (Equivalently, all $\\lambda_i>0$, so\n",
    "  $\\det\\mathbf A=\\prod_i\\lambda_i>0$ by the equation.)\n",
    "* **Gram and covariance matrices are PSD.** For any matrix $\\mathbf{X}$, the Gram\n",
    "  matrix $\\mathbf{G}=\\mathbf{X}^\\top\\mathbf{X}$ satisfies\n",
    "  $\\mathbf{x}^\\top\\mathbf{G}\\mathbf{x}=\\|\\mathbf{X}\\mathbf{x}\\|^2\\ge0$, so\n",
    "  $\\mathbf{G}\\succeq0$; it is PD exactly when $\\mathbf{X}$ has full column rank\n",
    "  (so that $\\mathbf{X}\\mathbf{x}=\\mathbf 0$ implies $\\mathbf{x}=\\mathbf 0$).\n",
    "  Covariance matrices, being Gram matrices of centered data divided by $n$,\n",
    "  inherit this: they are always PSD.\n",
    "* **The Hessian and minima.** A twice-differentiable function has a local minimum\n",
    "  at a critical point when its *Hessian*, the symmetric matrix $\\mathbf{H}$ of\n",
    "  second partial derivatives $h_{ij}=\\partial^2 f/\\partial x_i\\partial x_j$\n",
    "  (developed in the calculus chapter), is PD there: by\n",
    "  the equation the second-order Taylor term\n",
    "  $\\tfrac12\\mathbf{x}^\\top\\mathbf{H}\\mathbf{x}$ curves *upward* in every direction.\n",
    "  This is the second-order optimality test we use throughout optimization.\n",
    "\n",
    "Positive definiteness also has a purely *algorithmic* certificate, the\n",
    "*Cholesky factorization*: $\\mathbf{A}\\succ0$ if and only if\n",
    "$\\mathbf{A}=\\mathbf{L}\\mathbf{L}^\\top$ for a lower-triangular $\\mathbf{L}$ with\n",
    "positive diagonal (one direction is immediate: such an $\\mathbf{L}$ is\n",
    "invertible, so\n",
    "$\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}=\\|\\mathbf{L}^\\top\\mathbf{x}\\|^2>0$ for\n",
    "$\\mathbf{x}\\neq\\mathbf 0$). Attempting the factorization, far cheaper than an\n",
    "eigendecomposition, is the standard numerical test for definiteness. The same\n",
    "factorization is how multivariate Gaussians are sampled, an application we\n",
    "develop in that section.\n",
    "\n",
    "Geometrically, the sign pattern of the eigenvalues is the *shape* of the surface\n",
    "$z=\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}$, viewed through the eigenvector axes:\n",
    "all-positive eigenvalues give an upward **bowl** (PD, a strict minimum at the\n",
    "origin), a zero eigenvalue gives a flat-bottomed **trough** along that\n",
    "eigenvector (PSD, not strict), and a mix of signs gives a **saddle**\n",
    "(indefinite). the figure shows all three.\n",
    "\n",
    "![The quadratic-form surface $z=\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}$ takes its shape from the eigenvalue signs: a bowl when all $\\lambda_i>0$ (positive definite), a flat-bottomed trough when one $\\lambda_i=0$ (positive semidefinite), and a saddle for mixed signs (indefinite).](../img/mdl-la-psd.svg)\n",
    "\n",
    "Let us verify the eigenvalue test against the definition numerically. We build a\n",
    "Gram matrix $\\mathbf{G}=\\mathbf{X}^\\top\\mathbf{X}$ from a tall $\\mathbf{X}$ with\n",
    "full column rank (so it should be PD, all eigenvalues positive), and a second one\n",
    "whose columns are dependent (so it should be PSD with a zero eigenvalue)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-psd-gram"
   },
   "source": [
    "X_full = np.array([[1., 0.], [1., 1.], [0., 1.], [2., 1.]])  # rank 2\n",
    "X_dep = np.array([[1., 2.], [1., 2.], [0., 0.], [2., 4.]])   # col 2 = 2*col 1\n",
    "for X, name in [(X_full, 'full column rank'), (X_dep, 'dependent columns')]:\n",
    "    eigvals = np.linalg.eigvalsh(X.T @ X).round(4)\n",
    "    eigvals[np.abs(eigvals) < 1e-9] = 0.0  # clean up tiny round-off / -0.0\n",
    "    print(f'{name}: eigenvalues of X^T X = {eigvals}')"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "full column rank: eigenvalues of X^T X = [1.1459 7.8541]\ndependent columns: eigenvalues of X^T X = [ 0. 30.]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-9 fw=all -->\n",
    "\n",
    "The full-rank case has both eigenvalues strictly positive (PD); the dependent\n",
    "case has a zero eigenvalue (PSD but not PD), exactly as the rank condition\n",
    "predicts.\n",
    "\n",
    "### The Rayleigh Quotient: Eigenvalues as Extreme Stretches\n",
    "\n",
    "Identity the equation does more than test signs: it pins down the\n",
    "*largest* and *smallest* eigenvalues as optimization problems. Define the\n",
    "*Rayleigh quotient* of a symmetric $\\mathbf{A}$,\n",
    "\n",
    "$$\n",
    "R(\\mathbf{x}) = \\frac{\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}}{\\mathbf{x}^\\top\\mathbf{x}},\n",
    "\\qquad \\mathbf{x}\\neq\\mathbf 0 .\n",
    "$$\n",
    "\n",
    "**Proposition (Rayleigh).** *For symmetric $\\mathbf{A}$ with eigenvalues\n",
    "$\\lambda_1\\ge\\cdots\\ge\\lambda_n$,*\n",
    "\n",
    "$$\n",
    "\\lambda_1 = \\max_{\\mathbf{x}\\neq\\mathbf 0} R(\\mathbf{x}),\n",
    "\\qquad\n",
    "\\lambda_n = \\min_{\\mathbf{x}\\neq\\mathbf 0} R(\\mathbf{x}),\n",
    "$$\n",
    "\n",
    "*attained at the corresponding eigenvectors.*\n",
    "\n",
    "**Proof.** Restrict to unit vectors (scaling $\\mathbf{x}$ leaves $R$ unchanged).\n",
    "With $\\mathbf{y}=\\mathbf{W}^\\top\\mathbf{x}$ we have $\\|\\mathbf y\\|=\\|\\mathbf x\\|=1$\n",
    "and, from the equation,\n",
    "$\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}=\\sum_i\\lambda_i y_i^2$, a weighted average of\n",
    "the eigenvalues with non-negative weights $y_i^2$ summing to $1$. Any weighted\n",
    "average lies between the smallest and largest values, so\n",
    "$\\lambda_n\\le R(\\mathbf x)\\le\\lambda_1$; the bounds are achieved by putting all\n",
    "weight on a single coordinate, i.e. $\\mathbf x=\\mathbf w_1$ or\n",
    "$\\mathbf x=\\mathbf w_n$. $\\blacksquare$\n",
    "\n",
    "The proposition is easy to *watch*: sweep a unit vector once around the circle,\n",
    "record $R(\\mathbf{x})$ at each angle, and compare the extremes of the sweep to\n",
    "the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-eigendecomposition-the-rayleigh-quotient-eigenvalues-as-extreme-stretches"
   },
   "source": [
    "A = np.array([[3., 1.], [1., 2.]])\n",
    "thetas = np.linspace(0, np.pi, 721)   # R has period pi: x and -x agree\n",
    "xs = np.stack([np.cos(thetas), np.sin(thetas)])   # unit vectors, all angles\n",
    "R = np.sum(xs * (A @ xs), axis=0)     # R(x) = x^T A x  (denominator is 1)\n",
    "print(f'sweep of R:  min = {R.min():.6f}, max = {R.max():.6f}')\n",
    "print(f'eigvalsh:    {np.linalg.eigvalsh(A).round(6)}')"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "sweep of R:  min = 1.381967, max = 3.618033\neigvalsh:    [1.381966 3.618034]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-10 fw=all -->\n",
    "\n",
    "The sweep never escapes the eigenvalue interval, and its extremes reproduce\n",
    "$\\lambda_{\\min}$ and $\\lambda_{\\max}$ to the resolution of the angular grid;\n",
    "the two angles where they occur are precisely the eigenvector directions. Every\n",
    "direction in between mixes the two eigenvalues in proportion to its squared\n",
    "eigen-coordinates, exactly as the weighted-average proof says.\n",
    "\n",
    "Maximizing over subspaces of fixed dimension recovers *every* eigenvalue in\n",
    "between (the *Courant--Fischer min-max theorem*,\n",
    "$\\lambda_k=\\max_{\\dim S=k}\\min_{\\mathbf 0\\neq\\mathbf x\\in S}R(\\mathbf x)$),\n",
    "though we will not need the general form. The Rayleigh quotient reappears in\n",
    "that section as\n",
    "$\\sigma_1=\\max_{\\|\\mathbf x\\|=1}\\|\\mathbf{A}\\mathbf x\\|$ (the largest singular\n",
    "value) and as the variational characterization of PCA's top component; in\n",
    "optimization, $\\lambda_{\\max}(\\mathbf H)$ is the worst-case curvature, the\n",
    "constant $L$ that bounds a safe gradient-descent step size.\n",
    "\n",
    "A large spread of Hessian eigenvalues means the PD bowl above is *elongated*,\n",
    "and gradient descent is slow on such a bowl. The ratio\n",
    "$\\kappa=\\lambda_{\\max}/\\lambda_{\\min}$ is the *condition number*, defined in\n",
    "general and analyzed in that section.\n",
    "\n",
    "## Localizing and Computing Eigenvalues\n",
    "\n",
    "### Gershgorin Discs\n",
    "\n",
    "Eigenvalues are often difficult to reason with intuitively.\n",
    "If presented an arbitrary matrix, there is little that can be said\n",
    "about what the eigenvalues are without computing them.\n",
    "There is, however, one theorem that lets us *localize* them cheaply\n",
    "when the largest values are on the diagonal.\n",
    "\n",
    "**Theorem (Gershgorin** [@Gershgorin.1931]**).** *Let $\\mathbf{A}=(a_{ij})$ be any $n\\times n$ matrix\n",
    "and let $r_i=\\sum_{j\\neq i}|a_{ij}|$ be the off-diagonal absolute row sum. Let\n",
    "$\\mathcal{D}_i$ be the disc in the complex plane centered at $a_{ii}$ with radius\n",
    "$r_i$. Then every eigenvalue of $\\mathbf{A}$ lies in the union\n",
    "$\\bigcup_i\\mathcal{D}_i$.*\n",
    "\n",
    "**Proof (one line, once you spot the trick).** Let\n",
    "$\\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}$ with $\\mathbf{v}\\neq\\mathbf 0$, and pick\n",
    "the index $i$ of the entry of *largest magnitude*, so $|v_i|\\ge|v_j|$ for all\n",
    "$j$ (and $v_i\\neq0$). Reading off row $i$ of $\\mathbf{A}\\mathbf{v}=\\lambda\\mathbf{v}$,\n",
    "\n",
    "$$\n",
    "\\sum_j a_{ij}v_j = \\lambda v_i\n",
    "\\quad\\Longrightarrow\\quad\n",
    "(\\lambda - a_{ii})\\,v_i = \\sum_{j\\neq i} a_{ij} v_j .\n",
    "$$\n",
    "\n",
    "Divide by $v_i$ and take absolute values; since each $|v_j/v_i|\\le1$,\n",
    "\n",
    "$$\n",
    "|\\lambda - a_{ii}|\n",
    "  = \\Bigl|\\sum_{j\\neq i} a_{ij}\\,\\tfrac{v_j}{v_i}\\Bigr|\n",
    "  \\le \\sum_{j\\neq i} |a_{ij}| = r_i ,\n",
    "$$\n",
    "\n",
    "so $\\lambda\\in\\mathcal{D}_i$. $\\blacksquare$\n",
    "\n",
    "A useful corollary ties this section to the previous one. If $\\mathbf{A}$ is\n",
    "*strictly diagonally dominant* with positive diagonal ($a_{ii}>r_i$ for every\n",
    "$i$), then every disc $\\mathcal D_i$ lives strictly in the right half-plane, so\n",
    "$0$ is excluded from $\\bigcup_i\\mathcal D_i$ and no eigenvalue can be zero. For a\n",
    "symmetric such matrix the eigenvalues are real and positive, so it is positive\n",
    "definite, hence invertible (that section). Localization thus gives a\n",
    "no-computation proof of invertibility, which is exactly the property one needs to\n",
    "guarantee numerical solvers will succeed.\n",
    "\n",
    "Let's look at an example.\n",
    "Consider the matrix:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "1.0 & 0.1 & 0.1 & 0.1 \\\\\n",
    "0.1 & 3.0 & 0.2 & 0.3 \\\\\n",
    "0.1 & 0.2 & 5.0 & 0.5 \\\\\n",
    "0.1 & 0.3 & 0.5 & 9.0\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "We have $r_1 = 0.3$, $r_2 = 0.6$, $r_3 = 0.8$ and $r_4 = 0.9$.\n",
    "The matrix is symmetric, so all eigenvalues are real.\n",
    "This means that all of our eigenvalues will be in one of the ranges of\n",
    "\n",
    "$$[a_{11}-r_1, a_{11}+r_1] = [0.7, 1.3], $$\n",
    "\n",
    "$$[a_{22}-r_2, a_{22}+r_2] = [2.4, 3.6], $$\n",
    "\n",
    "$$[a_{33}-r_3, a_{33}+r_3] = [4.2, 5.8], $$\n",
    "\n",
    "$$[a_{44}-r_4, a_{44}+r_4] = [8.1, 9.9]. $$\n",
    "\n",
    "\n",
    "Performing the numerical computation shows\n",
    "that the eigenvalues are approximately $0.99$, $2.97$, $4.95$, $9.08$,\n",
    "all comfortably inside the ranges provided."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-gershgorin-circle-theorem"
   },
   "source": [
    "import numpy as onp\n",
    "A = onp.array([[1.0, 0.1, 0.1, 0.1],\n",
    "              [0.1, 3.0, 0.2, 0.3],\n",
    "              [0.1, 0.2, 5.0, 0.5],\n",
    "              [0.1, 0.3, 0.5, 9.0]])\n",
    "\n",
    "v, _ = onp.linalg.eig(A)\n",
    "v"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([9.08033648, 0.99228545, 4.95394089, 2.97343718])"
     },
     "metadata": {},
     "execution_count": 9
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-11 fw=all -->\n",
    "\n",
    "The theorem is most vivid in a picture. the figure draws the\n",
    "four discs in the complex plane and overlays the true eigenvalues, each one\n",
    "landing inside its disc.\n",
    "\n",
    "![The four Gershgorin discs for the matrix above, centered at the diagonal entries with radii the off-diagonal row sums. The true eigenvalues (crosses) each land inside a disc, exactly as the theorem guarantees.](../img/mdl-la-gershgorin.svg)\n",
    "\n",
    "In this way, eigenvalues can be approximated,\n",
    "and the approximations will be fairly accurate\n",
    "when the diagonal is\n",
    "significantly larger than all the other elements.\n",
    "\n",
    "### Power Iteration\n",
    "\n",
    "Now that we understand what eigenvectors are in principle,\n",
    "let's see how they can be used to provide a deep understanding\n",
    "of a problem central to neural network behavior: proper weight initialization.\n",
    "\n",
    "The full mathematical investigation of the initialization\n",
    "of deep neural networks is beyond the scope of the text,\n",
    "but we can see a toy version here to understand\n",
    "how eigenvalues can help us see how these models work.\n",
    "As we know, neural networks operate by interspersing layers\n",
    "of linear transformations with non-linear operations.\n",
    "For simplicity here, we will assume that there is no non-linearity,\n",
    "and that the transformation is a single repeated matrix operation $A$,\n",
    "so that the output of our model is\n",
    "\n",
    "$$\n",
    "\\mathbf{v}_{out} = \\mathbf{A}\\cdot \\mathbf{A}\\cdots \\mathbf{A} \\mathbf{v}_{in} = \\mathbf{A}^N \\mathbf{v}_{in}.\n",
    "$$\n",
    "\n",
    "For context, let's think of a generic ML problem, where we turn input data, like\n",
    "an image, into a prediction, like the probability the image is a picture of a\n",
    "cat. If repeated application of $\\mathbf{A}$ stretches a random vector out to be\n",
    "very long, then small changes in input will be amplified into large changes in\n",
    "output: tiny modifications of the input image would lead to vastly different\n",
    "predictions. On the flip side, if $\\mathbf{A}$ shrinks random vectors, then after\n",
    "many layers the vector shrinks to nothing and the output stops depending on the\n",
    "input. We need to walk the narrow line between growth and decay. What governs\n",
    "which happens? Repeatedly multiplying a random vector by $\\mathbf{A}$ and tracking\n",
    "its norm reveals that the *ratio of consecutive norms* stabilizes, and the\n",
    "eigendecomposition tells us exactly what it stabilizes to.\n",
    "\n",
    "Suppose $\\mathbf{A}$ is diagonalizable with a *strictly dominant* eigenvalue,\n",
    "$|\\lambda_1|>|\\lambda_2|\\ge\\cdots\\ge|\\lambda_n|$, and eigenvectors\n",
    "$\\mathbf{w}_1,\\ldots,\\mathbf{w}_n$. We call $\\lambda_1$ the *principal eigenvalue*\n",
    "and $\\mathbf{w}_1$ the *principal eigenvector*. Expand the input in the\n",
    "eigenbasis, $\\mathbf{v}_0=\\sum_i c_i\\mathbf{w}_i$, and apply $\\mathbf{A}$ a total\n",
    "of $k$ times. Because $\\mathbf{A}^k\\mathbf{w}_i=\\lambda_i^k\\mathbf{w}_i$,\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^k\\mathbf{v}_0\n",
    "  = \\sum_{i=1}^{n} c_i\\lambda_i^k\\mathbf{w}_i\n",
    "  = \\lambda_1^k\\Bigl(c_1\\mathbf{w}_1\n",
    "       + \\sum_{i\\ge2} c_i\\Bigl(\\tfrac{\\lambda_i}{\\lambda_1}\\Bigr)^{k}\\mathbf{w}_i\\Bigr).\n",
    "$$\n",
    "\n",
    "Every ratio $|\\lambda_i/\\lambda_1|<1$, so as $k\\to\\infty$ the bracket\n",
    "*converges to* $c_1\\mathbf{w}_1$: its tail decays to zero at the geometric rate\n",
    "$|\\lambda_2/\\lambda_1|$, set by the slowest-decaying term. Provided\n",
    "$c_1\\neq0$ (the input has *some* component\n",
    "along $\\mathbf{w}_1$), the iterate aligns with the principal eigenvector and its\n",
    "norm grows like $|\\lambda_1|^k$, so consecutive norms have ratio\n",
    "$\\to|\\lambda_1|$. This *is* the classical *power iteration* for the dominant\n",
    "eigenpair [@Golub.Van-Loan.1996]. the figure makes the\n",
    "convergence concrete for the small symmetric matrix $\\mathbf{B}=[[3,1],[1,2]]$,\n",
    "which has a genuine strictly dominant real eigenvalue\n",
    "$\\lambda_1=(5+\\sqrt5)/2\\approx3.618$: the renormalized iterates swing onto\n",
    "$\\mathbf{w}_1$, their *direction* gap closing at the rate\n",
    "$|\\lambda_2/\\lambda_1|=\\tfrac{3-\\sqrt5}{2}\\approx0.382$, while the norm ratio\n",
    "flattens to $\\lambda_1$ even faster, at the squared rate\n",
    "$(\\lambda_2/\\lambda_1)^2\\approx0.146$, because for a symmetric matrix the\n",
    "eigenvectors are orthogonal, so the norm is insensitive, to first order, to the\n",
    "remaining misalignment.\n",
    "\n",
    "![Power iteration on the symmetric $\\mathbf{B}=\\left(\\begin{smallmatrix}3&1\\\\1&2\\end{smallmatrix}\\right)$. Left: the renormalized iterates (lightening arrows) swing onto the principal eigenvector $\\mathbf{w}_1$ (orange), the direction gap closing at rate $|\\lambda_2/\\lambda_1|\\approx0.382$. Right: the ratio of consecutive norms converges to $\\lambda_1\\approx3.618$, its gap closing at the squared rate $(\\lambda_2/\\lambda_1)^2\\approx0.146$.](../img/mdl-la-power-iter.svg)\n",
    "\n",
    "Two caveats, both places where the short version above would mislead. First,\n",
    "**genericity**: the argument needs $c_1\\neq0$, but a\n",
    "random Gaussian input has $c_1=0$ with probability zero, so the method works\n",
    "almost surely. Second, **strict dominance**: when the two largest eigenvalues are\n",
    "a *complex-conjugate pair* of equal modulus, the dominant contribution *rotates*\n",
    "(see that section), and for a non-normal matrix the\n",
    "consecutive-norm ratio $\\|\\mathbf{A}^{k+1}\\mathbf{v}\\|/\\|\\mathbf{A}^{k}\\mathbf{v}\\|$\n",
    "can then oscillate *forever*, never converging. What always converges (for\n",
    "generic $\\mathbf{v}$) is the averaged growth rate\n",
    "$\\|\\mathbf{A}^{k}\\mathbf{v}\\|^{1/k}\\to\\max_i|\\lambda_i|$. The random matrix in\n",
    "the demonstration below happens to draw a real, strictly dominant eigenvalue, so\n",
    "even the simple ratio settles.\n",
    "\n",
    "The payoff is a one-line numerical fact: run a couple hundred iterations on a random\n",
    "$5\\times5$ matrix, read off the stabilized norm ratio, and compare it to the\n",
    "largest eigenvalue modulus. They agree."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "eigendecomposition-power-iteration"
   },
   "source": [
    "import numpy as onp\n",
    "onp.random.seed(8675309)\n",
    "A = onp.random.randn(5, 5)\n",
    "v = onp.random.randn(5, 1)\n",
    "for _ in range(200):\n",
    "    Av = A @ v\n",
    "    ratio = onp.linalg.norm(Av) / onp.linalg.norm(v)\n",
    "    v = Av / onp.linalg.norm(Av)\n",
    "rho = max(abs(onp.linalg.eigvals(A)))\n",
    "print(f'stabilized norm ratio = {ratio:.10f}   max|eigenvalue| = {rho:.10f}')"
   ],
   "execution_count": 10,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "stabilized norm ratio = 1.9744593215   max|eigenvalue| = 1.9744593215\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-12 fw=all -->\n",
    "\n",
    "The stabilized stretching factor matches the largest eigenvalue modulus to ten\n",
    "decimal places. One caveat: the random draw here happens\n",
    "to produce a *real*, strictly dominant eigenvalue, which is why the plain\n",
    "consecutive-norm ratio converges. Had the dominant eigenvalues been a\n",
    "complex pair, we would have read off the averaged rate\n",
    "$\\|\\mathbf{A}^k\\mathbf{v}\\|^{1/k}$ instead, per the caveat above. Either way,\n",
    "the quantity being measured has a name: the *spectral radius*\n",
    "$\\rho(\\mathbf A)=\\max_i|\\lambda_i|$, the largest\n",
    "eigenvalue modulus, which we examine next.\n",
    "\n",
    "#### Aside: PageRank and the Perron--Frobenius Theorem\n",
    "\n",
    "For an important class of matrices, power iteration is also the algorithm of\n",
    "choice for *computing* the dominant eigenvector, and that eigenvector is often\n",
    "the whole point. The **Perron--Frobenius theorem** [@Horn.Johnson.2012]\n",
    "says that a matrix with non-negative entries (and a mild connectivity\n",
    "condition, *irreducibility*: from every state you can reach every other) has a\n",
    "real, positive dominant eigenvalue whose eigenvector can be chosen with all\n",
    "entries non-negative, so it is interpretable as a distribution. Google's\n",
    "original **PageRank** [@Brin.Page.1998] is exactly this eigenvector.\n",
    "Model a random web surfer who at each step follows a random outgoing link.\n",
    "Before forming a transition matrix, repair every **dangling** page (one with no\n",
    "outgoing links), for example by replacing its empty column by the uniform\n",
    "distribution. The resulting click probabilities form a column-stochastic\n",
    "matrix $\\mathbf{P}$ (every column sums to $1$, and every entry is\n",
    "non-negative). Two lines show\n",
    "that $\\rho(\\mathbf{P})=1$. Because the columns sum to one,\n",
    "$\\mathbf{1}^\\top\\mathbf{P}=\\mathbf{1}^\\top$, so $1$ is an eigenvalue of\n",
    "$\\mathbf{P}^\\top$ and hence of $\\mathbf{P}$ (a matrix and its transpose share\n",
    "the characteristic polynomial, since $\\det(\\mathbf{M}^\\top)=\\det(\\mathbf{M})$).\n",
    "And no eigenvalue can exceed modulus $1$: applying the Gershgorin theorem of\n",
    "that section to\n",
    "$\\mathbf{P}^\\top$, whose rows are the columns of $\\mathbf{P}$, puts every\n",
    "eigenvalue in a disc centered at $p_{ii}\\ge0$ with radius $1-p_{ii}$, hence\n",
    "within modulus $p_{ii}+(1-p_{ii})=1$. The fraction of time the\n",
    "surfer spends on each page is the stationary distribution $\\boldsymbol\\pi$ with\n",
    "$\\mathbf{P}\\boldsymbol\\pi=\\boldsymbol\\pi$, the dominant eigenvector. One wrinkle:\n",
    "the repaired link matrix can still fail irreducibility or be periodic because\n",
    "the web contains disconnected pockets and cycles. The separate fix is\n",
    "*damping*: iterate the blended matrix\n",
    "$\\alpha\\mathbf{P}+(1-\\alpha)\\tfrac1n\\mathbf{1}\\mathbf{1}^\\top$ with\n",
    "$\\alpha\\approx0.85$, i.e. with probability $0.15$ the surfer teleports to a\n",
    "uniformly random page. Every entry of the blend is positive, so\n",
    "Perron--Frobenius applies *by construction*, and it guarantees the gap\n",
    "$|\\lambda_2|\\le\\alpha$ [@Haveliwala.Kamvar.2003], so power iteration\n",
    "converges fast. Power\n",
    "iteration started from the uniform distribution converges to $\\boldsymbol\\pi$, and\n",
    "that is precisely how PageRank is computed at web scale: forming $\\mathbf{P}$\n",
    "explicitly is unthinkable, but multiplying by it is cheap (and the teleportation\n",
    "term is rank one, so the damped multiply costs no more). The same\n",
    "dominant-eigenvector-by-power-iteration idea reappears across machine learning, in\n",
    "the stationary distribution of a Markov chain and in spectral clustering.\n",
    "\n",
    "## Spectral Radius, Stability, and Deep Networks\n",
    "\n",
    "The *spectral radius*\n",
    "$\\rho(\\mathbf A)=\\max_i|\\lambda_i|$ controls the **asymptotic exponential\n",
    "rate** of powers of one fixed matrix:\n",
    "\n",
    "$$\n",
    "\\lim_{k\\to\\infty}\\|\\mathbf A^k\\|^{1/k}=\\rho(\\mathbf A).\n",
    "$$\n",
    "\n",
    "This statement is weaker than saying that $\\|\\mathbf A^k\\|$ is approximately\n",
    "$\\rho(\\mathbf A)^k$ at every finite $k$. A defective Jordan block contributes\n",
    "polynomial factors, and a non-normal matrix can amplify some vectors greatly\n",
    "before its eventual asymptotic behavior takes over. For example,\n",
    "$\\mathbf A=\\left[\\begin{smallmatrix}1&1\\\\0&1\\end{smallmatrix}\\right]$ has\n",
    "$\\rho(\\mathbf A)=1$ but\n",
    "$\\mathbf A^k=\\left[\\begin{smallmatrix}1&k\\\\0&1\\end{smallmatrix}\\right]$.\n",
    "Thus rescaling to $\\rho(\\mathbf A)=1$ controls the exponential rate but does\n",
    "not guarantee bounded iterates. This distinction matters for initialization\n",
    "and raises the question of how large $\\rho$ typically is to begin with.\n",
    "\n",
    "### What Random Matrices Look Like\n",
    "\n",
    "Random matrix theory answers it. The *real Ginibre ensemble*\n",
    "[@Ginibre.1965] is the $n\\times n$ matrix whose entries are independent\n",
    "with mean zero and variance one. For this ensemble the rescaled eigenvalues\n",
    "$\\lambda_i/\\sqrt n$ fill the *unit* disk in the complex plane roughly uniformly\n",
    "as $n\\to\\infty$: the *circular law* [@Tao.Vu.2010]. Consequently\n",
    "$\\rho(\\mathbf A)/\\sqrt n\\to1$: the spectral radius grows like $\\sqrt n$, with\n",
    "the largest modulus sitting just outside the bulk of the disk. At finite $n$\n",
    "this is only a tendency (for the small $5\\times5$ example above we should not\n",
    "expect to land exactly on $\\sqrt5\\approx2.24$, since finite-size fluctuations\n",
    "are substantial), but the *scaling* with $\\sqrt n$ is what matters for\n",
    "initialization: it is why sensible schemes scale random weights by $1/\\sqrt n$\n",
    "(equivalently, by $1/\\sqrt{\\textrm{fan-in}}$) to keep $\\rho$ near $1$.\n",
    "\n",
    "A word of caution: the spectral radius is *not* the same as the largest\n",
    "*singular value* of such a matrix, which by the Marchenko--Pastur law sits near\n",
    "$2\\sqrt{n}$ at finite $n$. We return to singular values, and to the deeper\n",
    "connections between random-matrix spectra and initialization\n",
    "[@Pennington.Schoenholz.Ganguli.2017], in that section.\n",
    "\n",
    "Symmetric random matrices have their own law, and it is the one a data analyst\n",
    "meets first. Take $n$ samples of a $d$-dimensional vector with independent,\n",
    "mean-zero, unit-variance coordinates (pure noise, true covariance\n",
    "$\\mathbf{I}$) and form the sample covariance\n",
    "$\\hat{\\boldsymbol\\Sigma}=\\tfrac1n\\mathbf{X}^\\top\\mathbf{X}$. The relevant\n",
    "parameter is the *aspect ratio* $\\gamma=d/n$, the number of dimensions each\n",
    "sample must account for. One might hope the eigenvalues all sit near $1$. They\n",
    "spread instead: for $\\gamma$ held fixed, they cover the **Marchenko--Pastur\n",
    "bulk** $[(1-\\sqrt\\gamma)^2,\\,(1+\\sqrt\\gamma)^2]$ [@Marchenko.Pastur.1967];\n",
    "with $d=n/4$ the \"noise\" eigenvalues range across $[0.25, 2.25]$. This bulk is\n",
    "the *null hypothesis* of spectral data analysis: an eigenvalue inside it is\n",
    "indistinguishable from sampling noise, and only eigenvalues that escape *above*\n",
    "the bulk edge testify to real structure (the calibration behind the\n",
    "singular-value threshold of that section). It is also the\n",
    "baseline against which the heavy-tailed spectra of *trained* weight matrices,\n",
    "which spill far outside any Marchenko--Pastur bulk, are read as evidence of\n",
    "learned correlation in that section.\n",
    "\n",
    "The iterated map is exactly what a recurrent network runs. It computes a hidden state\n",
    "$\\mathbf{h}_t=\\phi(\\mathbf{W}\\mathbf{h}_{t-1}+\\cdots)$ by applying (almost) the\n",
    "same transformation at every time step, and backpropagation through time\n",
    "multiplies the per-step Jacobians:\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{h}_0}\\mathcal{L}\n",
    "   = \\mathbf{J}_1^\\top\\mathbf{J}_2^\\top\\cdots\\mathbf{J}_T^\\top\n",
    "   \\,\\nabla_{\\mathbf{h}_T}\\mathcal{L},\n",
    "\\qquad \\mathbf{J}_t=\\frac{\\partial\\mathbf{h}_t}{\\partial\\mathbf{h}_{t-1}} .\n",
    "$$\n",
    "\n",
    "This is a product of matrices applied to a vector. In the idealized linear case\n",
    "where every Jacobian is the same diagonalizable matrix $\\mathbf J$, its\n",
    "asymptotic exponential rate is governed by $\\rho(\\mathbf J)$. Actual recurrent\n",
    "networks multiply **different**, generally non-normal Jacobians. Their gradient\n",
    "norm is controlled directly by singular values of the product\n",
    "$\\mathbf J_T\\cdots\\mathbf J_1$, or asymptotically by Lyapunov exponents; the\n",
    "spectral radii of the individual factors do not determine it. Large singular\n",
    "values can produce exploding gradients and small ones vanishing gradients\n",
    "[@bengio1994learning; @Pascanu.Mikolov.Bengio.2013].\n",
    "\n",
    "This refined view explains why orthogonal/identity initialization\n",
    "[@Saxe.McClelland.Ganguli.2014], gradient clipping, and the gating\n",
    "mechanisms of LSTMs and GRUs help: they control finite-time amplification and\n",
    "information flow, not merely eigenvalue moduli. Singular values and condition\n",
    "numbers, introduced in that section, supply the right\n",
    "single-layer language; products require the same ideas applied along the whole\n",
    "trajectory.\n",
    "\n",
    "## Summary\n",
    "* Eigenvectors are directions a matrix stretches without rotating; eigenvalues\n",
    "  are the stretch factors. The unit circle maps to an ellipse whose axes, for a\n",
    "  symmetric matrix, lie along the eigenvectors with half-lengths $|\\lambda_i|$.\n",
    "* A matrix is *diagonalizable*\n",
    "  ($\\mathbf{A}=\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^{-1}$) exactly when its\n",
    "  eigenvectors span the space, i.e. geometric multiplicity equals algebraic\n",
    "  multiplicity for every eigenvalue; $n$ distinct eigenvalues guarantee it.\n",
    "* The eigenvalues reduce many operations to scalar ones: $\\mathbf{A}^n$ raises\n",
    "  them to the $n$, and because determinant is multiplicative and trace is\n",
    "  cyclic, similarity invariance gives $\\det\\mathbf{A}=\\prod_i\\lambda_i$ and\n",
    "  $\\operatorname{tr}\\mathbf{A}=\\sum_i\\lambda_i$.\n",
    "* Over the reals, every square matrix decomposes into pure stretch directions,\n",
    "  shear-like Jordan blocks (defective eigenvalues), and rotation--scaling\n",
    "  planes (complex-conjugate pairs). A Jordan block grows like $\\lambda^n$\n",
    "  times a polynomial in $n$, so a defective matrix can grow transiently even\n",
    "  when every $|\\lambda_i|<1$.\n",
    "* The *spectral theorem* guarantees real symmetric matrices an orthonormal\n",
    "  eigenbasis with real eigenvalues, $\\mathbf{A}=\\mathbf{W}\\boldsymbol{\\Lambda}\\mathbf{W}^\\top$;\n",
    "  applied to $\\mathbf{A}^\\top\\mathbf{A}$ it builds the SVD for *every* matrix.\n",
    "* *Positive (semi)definiteness* is decided by the sign of the eigenvalues, via\n",
    "  $\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}=\\sum_i\\lambda_i(\\mathbf{w}_i^\\top\\mathbf{x})^2$;\n",
    "  this governs Gram/covariance matrices and the Hessian minimum test. The Rayleigh\n",
    "  quotient identifies $\\lambda_{\\max}$ and $\\lambda_{\\min}$ as extreme stretches.\n",
    "* The Gershgorin Circle Theorem localizes the eigenvalues in discs around the\n",
    "  diagonal, giving cheap bounds and a no-computation invertibility test.\n",
    "* The *spectral radius* gives the asymptotic exponential rate of powers of one\n",
    "  fixed matrix, but Jordan blocks and non-normality can cause polynomial or\n",
    "  transient growth. Recurrent gradients involve products of changing\n",
    "  Jacobians, so their singular values and Lyapunov exponents, rather than the\n",
    "  factors' spectral radii alone, control exploding and vanishing gradients.\n",
    "\n",
    "## Exercises\n",
    "1. For the non-normal matrix in the equation, derive the\n",
    "   formula for $(\\mathbf A^k)_{12}$ by diagonalization or induction. Find the\n",
    "   integer $k$ at which this entry is largest, and compare it with the NumPy\n",
    "   result for $\\|\\mathbf A^k\\|_2$.\n",
    "1. Prove the equation from a unitary diagonalization.\n",
    "   Then construct a diagonalizable $2\\times2$ matrix with spectral radius\n",
    "   $1/2$ but one-step norm larger than an arbitrary prescribed $M>0$.\n",
    "1. What are the eigenvalues and eigenvectors of\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "2 & 1 \\\\\n",
    "1 & 2\n",
    "\\end{bmatrix}?\n",
    "$$\n",
    "1.  What are the eigenvalues and eigenvectors of the following matrix, and what is strange about this example compared to the previous one?\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "2 & 1 \\\\\n",
    "0 & 2\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "1. Without computing the eigenvalues, is it possible that the smallest eigenvalue of the following matrix is less than $0.5$? *Note*: this problem can be done in your head.\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix}\n",
    "3.0 & 0.1 & 0.3 & 1.0 \\\\\n",
    "0.1 & 1.0 & 0.1 & 0.2 \\\\\n",
    "0.3 & 0.1 & 5.0 & 0.0 \\\\\n",
    "1.0 & 0.2 & 0.0 & 1.8\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "1. Classify each matrix as positive definite, positive semidefinite, or\n",
    "   indefinite *from its eigenvalues*:\n",
    "$\\begin{bmatrix}2&1\\\\1&2\\end{bmatrix}$,\n",
    "$\\begin{bmatrix}1&2\\\\2&1\\end{bmatrix}$, and\n",
    "$\\begin{bmatrix}1&1\\\\1&1\\end{bmatrix}$.\n",
    "1. Prove that for any matrix $\\mathbf{X}$, the Gram matrix\n",
    "   $\\mathbf{X}^\\top\\mathbf{X}$ is positive semidefinite, and is positive definite\n",
    "   if and only if $\\mathbf{X}$ has full column rank. Conclude that the sum of two\n",
    "   positive semidefinite matrices is again positive semidefinite.\n",
    "1. Using $\\det\\mathbf{A}=\\prod_i\\lambda_i$ and\n",
    "   $\\operatorname{tr}\\mathbf{A}=\\sum_i\\lambda_i$, find both eigenvalues of\n",
    "   $\\begin{bmatrix}4&1\\\\1&4\\end{bmatrix}$ without forming the characteristic\n",
    "   polynomial directly. (*Hint:* you know their sum and product.)\n",
    "1. Verify the Rayleigh-quotient bounds for\n",
    "   $\\mathbf{A}=\\begin{bmatrix}3&1\\\\1&2\\end{bmatrix}$: compute\n",
    "   $\\mathbf{x}^\\top\\mathbf{A}\\mathbf{x}/\\mathbf{x}^\\top\\mathbf{x}$ for several unit\n",
    "   vectors and confirm the values stay between $\\lambda_{\\min}$ and\n",
    "   $\\lambda_{\\max}$, hitting the endpoints at the eigenvectors.\n",
    "1. The power-iteration tail in the equation decays at rate\n",
    "   $|\\lambda_2/\\lambda_1|$. For $\\mathbf{B}=\\begin{bmatrix}3&1\\\\1&2\\end{bmatrix}$,\n",
    "   compute this rate and predict roughly how many iterations are needed to reduce\n",
    "   the misalignment by a factor of $100$.\n",
    "1. Compute $\\mathbf{J}_2(\\lambda)^n$ by hand by writing\n",
    "   $\\mathbf{J}_2(\\lambda)=\\lambda\\mathbf{I}+\\mathbf{N}$ and using\n",
    "   $\\mathbf{N}^2=\\mathbf{0}$. For $\\lambda=0.99$, at which $n$ does\n",
    "   $\\|\\mathbf{J}_2(\\lambda)^n\\|$ peak? Check your prediction against the code\n",
    "   in that section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-eigendecomposition-md-13 fw=mxnet -->\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/411)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-mxnet",
   "display_name": "d2l (MXNet)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
