13_LA_eigen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Note:  This material largely follows the text \"Numerical Linear Algebra\" by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "init_cell": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Eigenproblems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "We will now consider eigenproblems of the form\n",
    "$$\n",
    "    A x = \\lambda x\n",
    "$$\n",
    "where $A \\in \\mathbb C^{m \\times m}$, $x \\in \\mathbb C^m$ and $\\lambda \\in \\mathbb C$.  The vector $x$ is known as the **eigenvector** and $\\lambda$ the **eigenvalue**.  The set of all eigenvalues is called the **spectrum** of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Characteristic Polynomial\n",
    "\n",
    "The **characteristic polynomial** of a matrix $A$ is the polynomial $\\mathcal{P}_A(z)$ defined by\n",
    "$$\n",
    "    \\mathcal{P}_A(z) = \\det(z I - A).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can then state the following theorem regarding the zeros of $\\mathcal{P}_A$ and the eigenvalues of $A$:\n",
    "\n",
    "*Theorem:* $\\lambda$ is an eigenvalue of $A$ if and only if $\\mathcal{P}_A(\\lambda) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Proof:* \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{If } \\lambda \\text{ is an eigenvalue of } A &\\Leftrightarrow \\text{ there is a non-zero vector } x \\text{ s.t. } \\lambda x - A x = 0 \\\\\n",
    "    &\\Leftrightarrow \\lambda I A \\text{ is singular (since }x\\text{ is a non-trivial vector in the null space of } \\lambda I - A) \\\\\n",
    "    &\\Leftrightarrow \\det(\\lambda I - A) = 0\n",
    "\\end{aligned}$$\n",
    "\n",
    "Note that this theorem implies that even though $A \\in \\mathbb R^{m \\times m}$ that $\\lambda \\in \\mathbb C$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Eigenvalue Decomposition\n",
    "Similar to QR factorization, an eigendecomposition is possible such that $A$ can be written as\n",
    "\n",
    "$$\n",
    "    A = X \\Lambda X^{-1}\n",
    "$$\n",
    "\n",
    "where $X$ is the matrix formed by the eigenvectors $x$ as its columns and $\\Lambda$ is a diagonal matrix with the eigenvalues along its diagonal.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This equation comes from the similar equation $A X = X \\Lambda$ which is of course related to the original problem statement.  This latter equation can be written out as\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   & A &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &  \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "        x_1 & x_2 & \\cdots & x_{m-1} & x_m \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &  \n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "        x_1 & x_2 & \\cdots & x_{m-1} & x_m \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &  \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix} \n",
    "        \\lambda_1 &   &   &   &   \\\\\n",
    "          & \\lambda_2 &   &   &   \\\\\n",
    "          &   & \\ddots &   &   \\\\\n",
    "          &   &   & \\lambda_{m-1} &   \\\\\n",
    "          &   &   &   & \\lambda_m\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Here we note that the eigenpair $(x_j, \\lambda_j)$ are matched as the $j$th column of $X$ and the $j$th element of $\\Lambda$ on the diagonal.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Algebraic multiplicity** is the number of times overall an eigenvalue repeats itself.\n",
    "\n",
    "**Geometric multiplicity** is defined as the number of linearly independent eigenvectors that belong to each eigenvalue.\n",
    "\n",
    "If the algebraic multiplicity is equal to the geometric multiplicity for all $\\lambda$ then we can say that there is a full eigenspace."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Computing Multiplicities\n",
    "\n",
    "Compute the geometric and algebraic multiplicities for the following matrices.  What is the relationship between the algebraic and geometric multiplicities?\n",
    "\n",
    "$$A = \\begin{bmatrix} \n",
    "    2 &   &  \\\\\n",
    "      & 2 &  \\\\\n",
    "      &   & 2 \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "$$B = \\begin{bmatrix} 2\n",
    "      & 1 &   \\\\\n",
    "      & 2 & 1 \\\\\n",
    "      &   & 2 \n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1. The characteristic polynomial of $A$ is\n",
    "  \n",
    "  $$\n",
    "      \\mathcal{P}_A(z) = (2 - z)(2 - z)(2 - z) = (2 - z)^3\n",
    "  $$\n",
    "  \n",
    "  so the eigenvalues are all $\\lambda = 2$ so we know the algebraic multiplicity is 3 of this eigenvalue.  The geometric multiplicity is determined by the number of linearly independent eigenvectors.  For this matrix we have three eigenvectors that are all linearly independent which happen to be the unit vectors in each direction (check!).  This means that the geometric multiplicity is also 3.\n",
    "\n",
    "1. The characteristic polynomial of $B$ is the same as $A$ so again we know $\\lambda = 2$ but now we need to be a bit careful about the eigenvectors.  In this case the only eigenvector is a scalar multiple of $e_1$ so the geometric multiplicity is 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Interpretations of the Eigenspace\n",
    "\n",
    "One way to interpret the eigenproblem is that of one that tries to find the subspaces of $\\mathbb C^m$ which act like scalar multiplication by $\\lambda$.  The eigenvectors associated with one eigenvalue then form a subspace of $S \\subseteq \\mathbb C^m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "When an eigenvalue has algebraic multiplicity that equals its geometric then it is called non-defective and otherwise defective.  This property is also inherited to the matrix so in the above example $A$ and $B$ are non-defective and defective matrices respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Determinant and Trace\n",
    "\n",
    "Two important properties of matrices have important relationships with their eigenvalues, namely the determinant and trace.  The determinant we have seen, the **trace** is defined as the sum of the elements on the diagonal of a matrix, in other words\n",
    "$$\n",
    "    \\text{tr}(A) = \\sum^m_{i=1} A_{ii}.\n",
    "$$\n",
    "\n",
    "The relationship between the determinant and the eigenvalues is not difficult to guess due to the nature of the characteristic polynomial.  The trace of a diagonal matrix is clear and provides another suggestion to the relationship.\n",
    "\n",
    "**Theorem:** The determinant $\\det(A)$ and trace $\\text{trace}(A)$ are equal to the product and sum of the eigenvalues of $A$ respectively counting algebraic multiplicity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Similarity Transformations\n",
    "\n",
    "The relationship between a matrix's eigenvalues and its determinant and trace are due to the special relationship between the eigenvalue decomposition and what is called similarity transformations.  A **similarity transformation** is defined as a transformation that takes A and maps it to $X^{-1} A X$ (assuming $X$ is non-singular).  Two matrices are said to be **similar** if there is a similarity transformation between them.  \n",
    "\n",
    "The most important property of similar matrices is that they have the same characteristic polynomial, eigenvalues, and multiplicities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This allows us to relate geometric and algebraic multiplicity as \n",
    "\n",
    "**Theorem:** The algebraic multiplicity of an eigenvalue $\\lambda$ is at least as great as its geometric multiplicity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Schur Factorization\n",
    "\n",
    "A **Schur factorization** of a matrix $A$ is defined as\n",
    "\n",
    "$$\n",
    "    A = Q T Q^\\ast\n",
    "$$\n",
    "\n",
    "where $Q$ is unitary and $T$ is upper-triangular.  In particular note that due do the structure of the resulting characteristic polynomial that $A$ and $T$ have identical eigenvalues. \n",
    "\n",
    " - Why is this?  \n",
    " - What relationship does this have to similarity transformations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " - The characteristic polynomials of $A$ and $T$ are identical.\n",
    " - The Schur factorization is a similarity transformation due to $Q$ being unitary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem:** Every matrix $A \\in \\mathbb C^{m \\times m}$ has a Schur factorization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that the above results imply the following\n",
    " - An eigen-decomposition $A = X \\Lambda X^{-1}$ exists if and only if $A$ is non-defective (it has a complete set of eigenvectors)\n",
    " - A unitary transformation $A = Q \\Lambda Q^\\ast$ exists if and only if $A$ is normal ($A^\\ast A = A A^\\ast$)\n",
    " - A Schur factorization always exists\n",
    " \n",
    "Note that each of these lead to a means for isolating the eigenvalues of a matrix and will be useful when considering algorithms for finding them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Condition Number of a Simple Eigenvalue\n",
    "\n",
    "Before we discuss a number of approaches to computing eigenvalues it good to consider what the condition number of a given eigenproblem is.  \n",
    "\n",
    "Let \n",
    "$$\n",
    "    Ax = \\lambda x\n",
    "$$ \n",
    "define the eigenvalue problem in question.  Here we will introduce a related problem \n",
    "$$\n",
    "    y^\\ast A = \\lambda y^\\ast\n",
    "$$ \n",
    "where $y$ is the **left eigenvector** and from before $x$ is the **right eigenvector**.  These vectors also can be shown to have the relationship $y^\\ast x \\neq 0$ for a simple eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now consider the perturbed problem\n",
    "$$\n",
    "    (A + \\delta A) (x + \\delta x) = (\\lambda + \\delta \\lambda) (x + \\delta x).\n",
    "$$\n",
    "Expanding this and throwing out quadratic terms and removing the eigenproblem we have\n",
    "$$\n",
    "    \\delta A x + A \\delta x = \\delta \\lambda x + \\lambda \\delta x.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Multiple both sides of the above by the left eigenvector and use $y^\\ast x \\neq 0$ to find\n",
    "$$\\begin{aligned}\n",
    "    y^\\ast \\delta A x + y^\\ast A \\delta x &= y^\\ast \\delta \\lambda x + y^\\ast \\lambda \\delta x \\\\\n",
    "    y^\\ast \\delta A x &= y^\\ast \\delta \\lambda x\n",
    "\\end{aligned}$$\n",
    "where we again use the slightly different definition of the eigenproblem.  We can then solve for $\\delta \\lambda$ to find\n",
    "$$\n",
    "    \\delta \\lambda = \\frac{y^\\ast \\delta A x}{y^\\ast x}\n",
    "$$\n",
    "meaning that the ratio between the dot-product of the left and right eigenvectors and the conjugate dot-product of the matrix $\\delta A$ then form a form of bound on the expected error in the simple eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computing Eigenvalues\n",
    "\n",
    "The most obvious approach to computing eigenvalues is a direct computation of the roots of the characteristic polynomial.  Unfortunately the following theorem suggests this is not a good way to compute eigenvalues:\n",
    "\n",
    "**Theorem:** For an $m \\geq 5$ there is a polynomial $\\mathcal{P}(z)$ of degree $m$ with rational coefficients that has a real root $\\mathcal{P}(z_0) = 0$ with the property that $z_0$ cannot be written using any expression involving rational numbers, addition, subtraction, multiplication, division, and $k$th roots.\n",
    "\n",
    "Not all is lost however, we just cannot use any direct methods to solve for the eigenvalues.  Instead we must use an iterative approach, in other words we want to construct a sequence that converges to the eigenvalues.  How does this relate to how we found roots previously?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Almost all approaches to computing eigenvalues do so through the computation of the Schur factorization.  The Schur factorization, as we have seen, will preserve the eigenvalues.  The steps to compute the Schur factorization are usually broken down into two steps\n",
    "1. Directly transform $A$ into a **Hessenberg** matrix, a matrix that contains zeros below its first sub-diagonal, directly using Householder reflections.\n",
    "1. Use an iterative method to change the sub-diagonal into all zeros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Hessenberg and Tridiagonal form\n",
    "\n",
    "What we want to do is construct a sequence of unitary matrices that turns $A$ into a Hessenberg matrix to start.  We can use Householder reflections to do this with the important distinction that we only want to remove zeros below the first sub-diagonal.  The sequence would look something like\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} & \\text{x} & \\text{x}\n",
    "    \\end{bmatrix} \\overset{Q_1}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x}& \\text{x} & \\text{x}\n",
    "    \\end{bmatrix} \\overset{Q_2}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x}& \\text{x} & \\text{x} \n",
    "    \\end{bmatrix} \\overset{Q_3}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x}& \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & 0 & \\text{x} & \\text{x}\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "so we have the sequence $Q^\\ast_1 A Q_1$.  Note we need both to preserve the entries of the first column that are not being transformed to zeros."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One important special case of this sequence of transformations is that if the matrix $A$ is hermitian (the matrix is its own conjugate transpose, $A = A^\\ast$, or symmetric in the real case) then the Hessenberg matrix is tridiagonal.\n",
    "\n",
    "We now will focus on how to formulate the iteration step of the eigenproblem.  We will also restrict our attention to symmetric, real matrices.  This implies that all eigenvalues will be real and have a complete set of orthogonal eigenvectors.  Generalizations can be made of many of the following algorithms but is beyond the scope of this class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Rayleigh Quotient and Inverse Iteration\n",
    "\n",
    "There are a number of classical approaches to computing the iterative step above which we will review here.  Inverse power iteration in particular is today still the dominant means of finding the eigenvectors once the eigenvalues are known."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Rayleigh Quotient\n",
    "\n",
    "The **Rayleigh quotient** of a vector $x \\in \\mathbb R^m$ is the scalar\n",
    "$$\n",
    "    r(x) = \\frac{x^T A x}{x^T x}.\n",
    "$$\n",
    "The importance of the Rayleigh quotient is made clear when we evaluate $r(x)$ at an eigenvector.  When this is the case the quotient evaluates to the corresponding eigenvalue.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The Rayleigh quotient can be motivated by asking the question, given an eigenvector $x$, what value $\\alpha$ acts most like an eigenvalue in an $\\ell_2$ sense:\n",
    "$$\n",
    "    \\min_\\alpha ||A x - \\alpha x||_2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This can be reformulated as a least-squares problem noting that $x$ is the \"matrix\", $\\alpha$ is the unknown vector (scalar) and $Ax$ is the right-hand side so we have\n",
    "$$\n",
    "    (x^T x) \\alpha = x^T (A x)\n",
    "$$\n",
    "which can be solved so that\n",
    "$$\n",
    "    \\alpha = r(x) = \\frac{x^T A x}{x^T x}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Power Iteration\n",
    "\n",
    "Power iteration is a straight forward approach to finding the eigenvector of the largest eigenvalue of $A$.  The basic idea is that the sequence\n",
    "$$\n",
    "    \\frac{x}{||x||}, \\frac{Ax}{||Ax||}, \\frac{A^2x}{||A^2x||}, \\frac{A^3x}{||A^3x||}, \\ldots\n",
    "$$\n",
    "will converge (although very slowly) to the desired eigenvector.\n",
    "\n",
    "We implement this method by initializing the algorithm with some vector $v$ with $||v|| = 1$.  We then apply the sequence of multiplications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The reason why this works can be seen by considering the initial vector $v$ as a linear combination of the orthonormal eigenvectors (which we have assumed exist) such that\n",
    "\n",
    "$$\n",
    "    v^{(0)} = a_1 q_1 + a_2 q_2 + \\cdots + a_m q_m.\n",
    "$$\n",
    "\n",
    "Multiplying $v^{(0)}$ by $A$ then leads to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    Av^{(0)} = v^{(1)} &= a_1 A q_1 + a_2 A q_2 + \\cdots + a_m A q_m \\\\\n",
    "    &= c_1 (a_1 \\lambda_1 q_1 + a_2 \\lambda_2 q_2 + \\cdots + a_m \\lambda_m q_m) \\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "where $c_1$ is some constant due to the fact the eigenvectors are not uniquely specified.  Repeating this $k$ times we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    Av^{(k-1)} = v^{(k)} &= a_1 A^k q_1 + a_2 A^k q_2 + \\cdots + a_m A^k q_m \\\\\n",
    "    &= c_k (a_1 \\lambda_1^k q_1 + a_2 \\lambda_2^k q_2 + \\cdots + a_m \\lambda_m^k q_m) \\\\\n",
    "    &= c_k \\lambda_1^k \\left(a_1 q_1 + a_2 \\frac{\\lambda_2^k}{\\lambda_1^k} q_2 + \\cdots + a_m \\frac{\\lambda_m^k}{\\lambda_1^k} q_m \\right)\n",
    "\\end{aligned}$$\n",
    "\n",
    "Since $\\lambda_1 > \\lambda_i$ for all $i \\neq 1$ then in the limit the terms $\\lambda_2^k / \\lambda_1^k$ will approach zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Inverse Iteration\n",
    "\n",
    "Inverse iteration uses a similar approach with the difference being that we can use it to find any of the eigenvectors for the matrix $A$.  \n",
    "\n",
    "Consider the matrix \n",
    "\n",
    "$$\n",
    "    (A - \\mu I)^{-1},\n",
    "$$ \n",
    "\n",
    "the eigenvectors of this matrix are the same as $A$ with the eigenvalues \n",
    "\n",
    "$$\n",
    "    (\\lambda_j - \\mu)^{-1}\n",
    "$$ \n",
    "\n",
    "where $\\lambda_j$ are the eigenvalues of $A$.  \n",
    "\n",
    "If $\\mu$ is close to a particular $\\lambda_j$, say $\\lambda_J$, then \n",
    "\n",
    "$$\n",
    "    (\\lambda_J - \\mu)^{-1}\n",
    "$$ \n",
    "\n",
    "will be larger than any of the other $(\\lambda_j - \\mu)^{-1}$.  In this way we effectively have picked out the eigenvalue we want to consider in the power iteration!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Rayleigh Quotient Iteration\n",
    "\n",
    "By themselves the above approaches are not particularly useful but combining them we can iterate back and forth to find the eigenvalue, eigenvector pair:\n",
    "1. Compute the Rayleigh quotient and find an estimate for $\\lambda_j$\n",
    "1. Compute one step of inverse iteration to approximate $x_j$\n",
    "1. Repeat..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "m = 3\n",
    "A = numpy.array([[2, 1, 1], [1, 3, 1], [1, 1, 4]])\n",
    "\n",
    "num_steps = 10\n",
    "v = numpy.empty((num_steps, m))\n",
    "lam = numpy.empty(num_steps)\n",
    "\n",
    "v[0, :] = numpy.array([1, 1, 1])\n",
    "v[0, :] = v[0, :] / numpy.linalg.norm(v[0, :], ord=2)\n",
    "lam[0] = numpy.dot(v[0,:], numpy.dot(A, v[0, :]))\n",
    "for k in range(1, num_steps):\n",
    "    w = numpy.linalg.solve(A - lam[k-1] * numpy.identity(m), v[k-1, :])\n",
    "    v[k, :] = w / numpy.linalg.norm(w, ord=2)\n",
    "    lam[k] = numpy.dot(v[k,:], numpy.dot(A, v[k, :]))\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(range(10), numpy.abs(lam - numpy.linalg.eigvals(A)[0]), 'o')\n",
    "\n",
    "axes.set_title(\"Convergence of eigenvalue\")\n",
    "axes.set_xlabel(\"Step\")\n",
    "axes.set_ylabel(\"Error\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR Algorithm\n",
    "\n",
    "The most basic use of a $QR$ factorization to find eigenvalues is to iteratively compute the factorization and multiply the resulting $Q$ and $R$ in the reverse order.  This sequence will eventually converge to the Schur decomposition of the matrix $A$.\n",
    "\n",
    "Code this up and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%precision 6\n",
    "m = 3\n",
    "A = numpy.array([[2, 1, 1], [1, 3, 1], [1, 1, 4]])\n",
    "MAX_STEPS = 10\n",
    "\n",
    "for i in range(MAX_STEPS):\n",
    "    Q, R = numpy.linalg.qr(A)\n",
    "    A = numpy.dot(R, Q)\n",
    "    print()\n",
    "    print(\"A(%s) =\" % (i))\n",
    "    print(A)\n",
    "\n",
    "print()\n",
    "print(\"True eigenvalues: \")\n",
    "print(numpy.linalg.eigvals(A))\n",
    "print()\n",
    "print(\"Computed eigenvalues: \")\n",
    "for i in range(m):\n",
    "    print(A[i, i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So why does this work?  The first step is to find the $QR$ factorization of $A^{(k-1)}$ which is equivalent to finding\n",
    "\n",
    "$$\n",
    "    (Q^{(k)})^T A^{(k-1)} = R^{(k)}\n",
    "$$\n",
    "\n",
    "and multiplying on the right leads to\n",
    "\n",
    "$$\n",
    "    (Q^{(k)})^T A^{(k-1)} Q^{(k)} = R^{(k)} Q^{(k)}.\n",
    "$$\n",
    "\n",
    "In this way we can see that this is a similarity transformation of the matrix $A^{(k-1)}$ since the $Q^{(k)}$ is an orthogonal matrix ($Q^{-1} = Q^T$). This of course is not a great idea to do directly but works great in this case as we iterate to find the upper triangular matrix $R^{(k)}$ which is exactly where the eigenvalues appear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In practice this basic algorithm is modified to include a few additions:\n",
    "\n",
    "1. Before starting the iteration $A$ is reduced to tridiagonal form.\n",
    "1. Motivated by the inverse power iteration we observed we instead consider a shifted matrix $A^{(k)} - \\mu^{(k)} I$ for factoring.  The $\\mu$ picked is related to the estimate given by the Rayleigh quotient.  Here we have\n",
    "\n",
    "$$\n",
    "    \\mu^{(k)} = \\frac{(q_m^{(k)})^T A q_m^{(k)}}{(q_m^{(k)})^T q_m^{(k)}} = (q_m^{(k)})^T A q_m^{(k)}.\n",
    "$$\n",
    "\n",
    "1. Deflation is used to reduce the matrix $A^{(k)}$ into smaller matrices once (or when we are close to) finding an eigenvalue to simplify the problem.\n",
    "\n",
    "This has been the standard approach until recently for finding eigenvalues of a matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Alternatives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Jacobi\n",
    "\n",
    "Jacobi iteration employs the idea that we know the eigenvalues of a matrix of size equal to or less than 4 (we know the roots of the characteristic polynomial directly).  Jacobi iteration therefore attempts to break the matrix down into at most 4 by 4 matrices along the diagonal via a series of similarity transformations based on only diagonalizing sub-matrices 4 by 4 or smaller."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bisection \n",
    "\n",
    "It turns out if you do not want all of the eigenvalues of a matrix that using a bisection method to find some subset of the eigenvalues is often the most efficient way to get these.  This avoids the pitfall of trying to find the eigenvalues via other root-finding approaches by only needing evaluations of the function and if a suitable initial guess is provided can find the eigenvalue quickly that is closest to the initial bracket provided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Divide-and-conquer\n",
    "\n",
    "This algorithm is actually the one used most often used if both eigenvalues and eigenvectors are needed and performs up to twice as fast as the $QR$ approach.  The basic idea is to split the matrix into two pieces at every iteration by introducing zeros on the appropriate off-diagonals which neatly divides the problem into two pieces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Arnoldi and Lanczos Iteration\n",
    "\n",
    "Krylov subspace methods (which we will unfortunately not cover) are another approach to finding eigenvalues of a matrix.  These methods generally use some piece of the $QR$ approach outlined above and are extremely effective at finding the \"extreme\" eigenvalues of the matrix."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.0"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}