13_LA_eigen_2020.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\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Eigenproblems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "We will now consider eigenproblems of the form\n",
    "\n",
    "$$\n",
    "    A \\mathbf{x} = \\lambda \\mathbf{x}\n",
    "$$\n",
    "\n",
    "where $A \\in \\mathbb C^{m \\times m}$, $\\mathbf{x} \\in \\mathbb C^m$ and $\\lambda \\in \\mathbb C$.  The vector $\\mathbf{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": [
    "### The basics\n",
    "\n",
    "The eigenproblem\n",
    "$$\n",
    "    A \\mathbf{x} = \\lambda \\mathbf{x}\n",
    "$$\n",
    "can be rewritten as\n",
    "\n",
    "$$\n",
    "    ( A - \\lambda I)\\mathbf{x} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "which implies that the eigenvectors are in the Null space of $A-\\lambda I$. \n",
    "\n",
    "However for this matrix to have a non-trivial Null space, requires that $A-\\lambda I$ is singular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Characteristic Polynomial\n",
    "\n",
    "If $A-\\lambda I$ is singular, it follows that\n",
    "\n",
    "$$\n",
    "    \\det( A - \\lambda I) = {\\cal P}_A(\\lambda) = 0\n",
    "$$\n",
    "\n",
    "where ${\\cal P}_A(\\lambda)$ can be shown to be a $m$th order polynomial in $\\lambda$ known as  the **characteristic polynomial** of a matrix $A$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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$.\n",
    "\n",
    "i.e. the eigenvalues are the roots  of ${\\cal P}_A(\\lambda)$, and therefore there are exactly $m$ eigenvalues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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": [
    "## Computing Eigenvalues\n",
    "\n",
    "In basic linear algebra classes we usually find the eigenvalues by directly calculating the roots of ${\\cal P}_A(\\lambda)$ which can work for low-degree polynomials.    Unfortunately the following theorem (due to Galois) 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",
    "I.e., there is no way to find the roots of a polynomial of degree $>4$ in a deterministic, fixed number of steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Not all is lost however!\n",
    "\n",
    "We just must use an iterative approach where we construct a sequence that converges to the eigenvalues.  \n",
    "\n",
    "**Some Questions**\n",
    "* How does this relate to how we found roots previously?\n",
    "* Why will it still be difficult to use our rootfinding routines to find Eigenvalues?\n",
    "\n",
    "We will return to how we actually find Eigenvalues (and roots of polynomials) after a bit more review"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Eigenvalue Factorization and Diagonalization\n",
    "\n",
    "Given that there are exactly $m$ (possibly repeated) eigenvalues for a $m\\times m$ system,  the eigenproblem is really more correctly written as \n",
    "\n",
    "$$\n",
    "    A\\mathbf{x}_i = \\lambda_i \\mathbf{x}_i, \\quad i=1,2,\\ldots,m\n",
    "$$\n",
    "\n",
    "Or in Matrix form as\n",
    "$$\n",
    "    AX = X\\Lambda\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": [
    "Expanded, $A X = X \\Lambda$ looks like:\n",
    "$$\n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   & A &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &  \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "        \\mathbf{x}_1 & \\mathbf{x}_2 & \\cdots  & \\mathbf{x}_m \\\\\n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &  \n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix} \n",
    "          &   &   &   &   \\\\\n",
    "          &   &   &   &   \\\\\n",
    "        \\mathbf{x}_1 & \\mathbf{x}_2 & \\cdots & \\mathbf{x}_{m-1} & \\mathbf{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 $(\\mathbf{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": "slide"
    }
   },
   "source": [
    "#### Diagonalization\n",
    "\n",
    "Eigenproblems can always be written as \n",
    "$$\n",
    "    AX = X\\Lambda\n",
    "$$\n",
    "\n",
    "However, if there are a linearly independent set of eigenvectors then the matrix $X$ is invertible (why?).  Under this condition we can transform $A$ into the diagonal matrix $\\Lambda$ by\n",
    "\n",
    "$$\n",
    "    \\Lambda = X^{-1}A X\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Factorization \n",
    "\n",
    "or more usefully,  we can rewrite $A$ as a product of more useful matrices\n",
    "\n",
    "$$\n",
    "    A = X\\Lambda X^{-1}\n",
    "$$\n",
    "\n",
    "And use it in similar ways to how we use other factorization such as $A=QR$ and $A=LU$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Rules for diagonalizability\n",
    "\n",
    "Unfortunately, not all matrices can be diagonalized (i.e. don't have a full linearly independent set of eigenvectors that span $\\mathbb{R}^m$ or $\\mathbb{C}^m$).  Here are the rules.\n",
    "\n",
    "A matrix can be factored as $A = X\\Lambda X^{-1}$ if\n",
    "\n",
    "1. all eigenvalues are distinct (never repeat). These are known as simple eigenvalues\n",
    "\n",
    "2. Eigenvalues repeat, but for every repeated eigenvalue there can be found a linearly independent eigenvector\n",
    "\n",
    "3. The matrix is Hermitian ($A^\\ast = A$, or if real, then symmetric $A^T = A$).  In this case the eigenvalues are always real and the eigenvectors can always be chosen orthonormal. In this special case $X=Q$ and\n",
    "\n",
    "$$\n",
    "    A = Q\\Lambda Q^\\ast\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The issue of repeated eigenvalues can be made a bit more precise by defining\n",
    "\n",
    "*   **Algebraic multiplicity**: the number of times  an eigenvalue is repeated\n",
    "\n",
    "*   **Geometric multiplicity**: the number of linearly independent eigenvectors corresponding  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 and the matrix is diagonalizable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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",
    "Generally, we say any two matrices $A$ and $B$ are **similar** if they can be related through an invertible matrix $M$ as\n",
    "\n",
    "$$\n",
    "    A = M^{-1} B M\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Example**\n",
    "\n",
    "a diagonalizable matrix $A$ is similar to the diagonal matrix $\\Lambda$ through the invertible matrix $X$\n",
    "\n",
    "$$\n",
    "    A = X\\Lambda X^{-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem**:  If $A$ and $B$ are similar matrices, they have the same eigenvalues and their eigenvectors are related through an invertible matrix $M$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Proof**: Let\n",
    "\n",
    "$$\n",
    "    B = M A M^{-1}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    BM = MA\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "if $A\\mathbf{x} = \\lambda\\mathbf{x}$  then\n",
    "$$\n",
    "    BM\\mathbf{x} = M A\\mathbf{x} = \\lambda M\\mathbf{x}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    B\\mathbf{y} = \\lambda\\mathbf{y}\n",
    "$$\n",
    "\n",
    "which shows that $\\lambda$ is also an eigenvalue of $B$ with corresponding eigenvector $\\mathbf{y} = M\\mathbf{x}$"
   ]
  },
  {
   "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.  Because $Q^\\ast=Q^{-1}$ (for square unitary matrices). It follows directly that $A$ and $T$ are similar.  \n",
    "\n",
    "*  Good News!  $T$ is upper triangular so its eigenvalues can just be read of the diagonal\n",
    "*  Bad News! There is no deterministic way to calculate $T$ as that would violate Galois theory of polynomials"
   ]
  },
  {
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "subslide"
    }
   },
   "source": [
    "## Computing Eigenvalues\n",
    "\n",
    "Almost all useful 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. This is as close to triangular that you can get by direct similarity transformations.\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 similarity transformations matrices that turns $A$ into a Hessenberg matrix with the same eigenvalues as $A$.  We use Householder reflections to do this with the important distinction that we only want to remove zeros below the first sub-diagonal. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\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{H_1^\\ast A_0 H_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{H_2^\\ast A_1H_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{H_3^\\ast A_2H_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 $H = Q^\\ast A Q$ which has the same eigenvalues as the original matrix $A$.  \n",
    "\n",
    "**Question**?  Why can't we just use Householder to take $A\\rightarrow T$ like we did for the $QR$?"
   ]
  },
  {
   "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 straightforward 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": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Demo\n",
    "# generate a random symmetric matrix\n",
    "A = numpy.random.rand(3,3)\n",
    "A = 0.5*( A + A.T)\n",
    "print('A=\\n{}'.format(A))\n",
    "\n",
    "lams = numpy.linalg.eigvals(A)\n",
    "print('\\nLambda = {}'.format(lams))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#the Rayleigh Quotient\n",
    "def rayleighq(A, x):\n",
    "    return numpy.dot(x.T, A.dot(x))/numpy.dot(x.T, x)\n",
    "\n",
    "def power_iteration(A, tol=1.e-6):\n",
    "    \"\"\" power_iteration to find the largest eigenvector and corresponding eigenvalue \n",
    "    \n",
    "    parameters: \n",
    "    -----------\n",
    "    \n",
    "    A: ndarray (square)\n",
    "        m x m matrix\n",
    "        \n",
    "    tol: float\n",
    "        stopping criteria for iteration.  \n",
    "        iteration will cease when ||x_{i+1} - x_{i}|| < tol or\n",
    "        MAX_ITS exceeded\n",
    "        \n",
    "    returns:\n",
    "    --------\n",
    "    x: ndarray \n",
    "        array of iterates of the eigenvector\n",
    "    r: ndarray\n",
    "    \"\"\"\n",
    "    \n",
    "    MAX_ITS = 100\n",
    "    m = A.shape[0]\n",
    "    \n",
    "    x = numpy.empty((MAX_ITS,m))\n",
    "    lam = numpy.empty(MAX_ITS)\n",
    "    res = numpy.empty(MAX_ITS)\n",
    "    \n",
    "    # generate a random unit vector\n",
    "    x0 = numpy.random.rand(A.shape[0])\n",
    "    x[0,:] = x0/numpy.linalg.norm(x0, ord=2)\n",
    "    lam[0] = rayleighq(A, x0)\n",
    "    \n",
    "    for i in range(1,MAX_ITS+1):\n",
    "        xi = A.dot(x[i-1,:])\n",
    "        x[i,:] = xi/numpy.linalg.norm(xi, ord=2)\n",
    "        lam[i] = rayleighq(A, x[i,:])\n",
    "        res[i-1] = numpy.abs(lam[i] - lam[i - 1])/numpy.abs(lam[i])\n",
    "        if res[i-1] < tol:\n",
    "            break\n",
    "    \n",
    "    if i == MAX_ITS:\n",
    "        warnings.warn('Maximum iterations exceeded')\n",
    "    \n",
    "    x.resize(i+1,m)\n",
    "    lam.resize(i+1)\n",
    "    res.resize(i)\n",
    "    return x, lam, res\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x, r, res = power_iteration(A, tol=1.e-8)\n",
    "print('{} Iterations'.format(len(r)))\n",
    "print('x = {}'.format(x[-1]))\n",
    "print('Eigenvalue = {}'.format(r[-1]))\n",
    "print('eigs(A) = {}'.format(numpy.linalg.eigvals(A)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.plot(r,'o',markersize=8)\n",
    "lam_max = numpy.max(numpy.linalg.eigvals(A))\n",
    "axes.plot(lam_max*numpy.ones(r.shape),'k--')\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iteration', fontsize=16)\n",
    "axes.set_ylabel('Rayleigh Quotient',fontsize=16)\n",
    "\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "axes.semilogy(res,'o-',markersize=8)\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iteration', fontsize=16)\n",
    "axes.set_ylabel('Residual',fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "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",
    "here $c_1$ is some constant due to the fact the eigenvectors are not uniquely specified.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "wRepeating 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 and on normalization $v^{(k)}/||v^{(k)}||\\rightarrow \\mathbf{q}_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Inverse Iteration with shifts\n",
    "\n",
    "Inverse iteration with shifts 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",
    "**Some Preliminaries**:  inverse and shift rules of Eigenvalues\n",
    "\n",
    "Show that if $\\mathbf{x}$ is an eigenvector of $A$ with eigenvalue $\\lambda$, then\n",
    "\n",
    "* $\\mathbf{x}$ is an eigenvector of $A^{-1}$ with eigenvalue $1/\\lambda$\n",
    "* $\\mathbf{x}$ is an eigenvector of $A -\\sigma I$ with eigenvalue $\\lambda - \\sigma$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So...\n",
    "\n",
    "If we want to find the smallest eigenvalue we can consider the power method on $A^{-1}$,\n",
    "\n",
    "But we really don't want to find $A^{-1}$ which is expensive,  instead we can  do the equivalent iteration\n",
    "\n",
    "```python\n",
    "x[0] = x0\n",
    "for i in range(MAX_ITS):\n",
    "    solve A w[i] = x[i]\n",
    "        x[i+1] = w[i]/norm(w[i]) \n",
    "```\n",
    "\n",
    "as $\\mathbf{w}_i = A^{-1}\\mathbf{x}_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and if we want to find the eigenvalue closest to some number $\\mu$ we can apply the power method to \n",
    "\n",
    "$$\n",
    "    (A - \\mu I)^{-1},\n",
    "$$ \n",
    "\n",
    "the eigenvectors of this matrix are \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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def rayleigh_quotient_iteration(A, tol=1.e10):\n",
    "    \"\"\" rayleigh quotient iteration to find eigenvalues\n",
    "    parameters: \n",
    "    -----------\n",
    "    \n",
    "    A: ndarray (square)\n",
    "        m x m matrix\n",
    "        \n",
    "    tol: float\n",
    "        stopping criteria for iteration.  \n",
    "        iteration will cease when ||x_{i+1} - x_{i}|| < tol or\n",
    "        MAX_ITS exceeded\n",
    "        \n",
    "    returns:\n",
    "    --------\n",
    "    x: ndarray \n",
    "        array of iterates of the eigenvector\n",
    "    r: ndarray\n",
    "    \"\"\"\n",
    "    \n",
    "    MAX_ITS = 100\n",
    "    m = A.shape[0]\n",
    "    \n",
    "    x = numpy.empty((MAX_ITS,m))\n",
    "    lam = numpy.empty(MAX_ITS)\n",
    "    res = numpy.empty(MAX_ITS)\n",
    "    \n",
    "    I = numpy.eye(m)\n",
    "    \n",
    "    # generate a random unit vector\n",
    "    x0 = numpy.random.rand(A.shape[0])\n",
    "    x[0,:] = x0/numpy.linalg.norm(x0, ord=2)\n",
    "    lam[0] = rayleighq(A, x0)\n",
    "    \n",
    "    for i in range(1,MAX_ITS+1):\n",
    "        # this is the only different line\n",
    "        w = numpy.linalg.solve(A - lam[i-1]*I, x[i-1,:])\n",
    "        x[i,:] = w/numpy.linalg.norm(w, ord=2)\n",
    "        lam[i] = rayleighq(A, x[i,:])\n",
    "        res[i-1] = numpy.abs(lam[i] - lam[i-1])/numpy.abs(lam[i])\n",
    "        if res[i-1] < tol:\n",
    "            break\n",
    "    if i == MAX_ITS:\n",
    "        warnings.warn('Maximum iterations exceeded')\n",
    "    \n",
    "    x.resize(i+1,m)\n",
    "    lam.resize(i+1)\n",
    "    res.resize(i)\n",
    "    return x, lam, res\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x, r, res = rayleigh_quotient_iteration(A, tol=1.e-10)\n",
    "print('{} Iterations'.format(len(r)))\n",
    "print('x = {}'.format(x[-1]))\n",
    "print('Rayleighquotient = {}'.format(r[-1]))\n",
    "print('eigs(A) = {}'.format(numpy.linalg.eigvals(A)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.plot(r,'o',markersize=8)\n",
    "lam_max = numpy.max(numpy.linalg.eigvals(A))\n",
    "axes.plot(lam_max*numpy.ones(r.shape),'k--')\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iteration', fontsize=16)\n",
    "axes.set_ylabel('Rayleigh Quotient',fontsize=16)\n",
    "\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "axes.semilogy(res,'o-',markersize=8)\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iteration', fontsize=16)\n",
    "axes.set_ylabel('Residual',fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR/RQ Algorithm\n",
    "\n",
    "All of the above methods pick out at most a few eigenvalues at a time.  However it turns out we can use the $QR$ algorithm, to iterate towards the Schur factorization and find all the eigenvalues simultaneously. \n",
    "\n",
    "The simplest algorithm just iterates \n",
    "```python\n",
    "    while not converged:\n",
    "        Q, R = numpy.linalg.qr(A)\n",
    "        A = R.dot(Q)        \n",
    "```\n",
    "calculating the $QR$ factorization of $A$, then forming a new $A=RQ$,  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": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%precision 6\n",
    "m = 3\n",
    "A = numpy.array([[2, 1, 1], [1, 3, 1], [1, 1, 4]])\n",
    "print(numpy.linalg.eigvals(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import hessenberg\n",
    "A = hessenberg(A)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "MAX_STEPS = 10\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print()\n",
    "print(\"True eigenvalues: \")\n",
    "print(numpy.linalg.eigvals(A))\n",
    "print()\n",
    "print(\"Computed eigenvalues: \")\n",
    "print(numpy.diag(A))"
   ]
  },
  {
   "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)}$ or\n",
    "\n",
    "$$\n",
    "    A^{(k-1)} = Q^{(k)}R^{(k)}\n",
    "$$\n",
    "\n",
    "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 Hessenberg 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",
    "    \\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",
    "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": "subslide"
    }
   },
   "source": [
    "### Application:  Finding the roots of a polynomial\n",
    "\n",
    "Numpy has a nice function called roots which returns the $n$ roots of a $n$th degree polynomial\n",
    "\n",
    "$$\n",
    "   p(x) = c_0 + c_1 x + c_2 x^2 + \\ldots + c_n x^n\n",
    "$$\n",
    "\n",
    "described by a $n+1$ vector of coefficients $\\mathbf{c}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "c = numpy.array([ -1, -1, 1])\n",
    "r = numpy.roots(c)\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c = numpy.random.rand(6)\n",
    "r = numpy.roots(c)\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This routine, does not try and actually find the roots of a high-order polynomial,  instead it actually calculates the eigenvalues of a **companion matrix** $C$ whose characteristic polynomial $P_C(\\lambda)$ is the **monic** polynomial $p_m(x) = p(x)/c_n$.\n",
    "\n",
    "It can be shown that this matrix can be constructed as ([see e.g.](https://en.wikipedia.org/wiki/Companion_matrix))\n",
    "$$\n",
    "C(p)=\\begin{bmatrix}\n",
    "0 & 0 & \\dots & 0 & -c_0 \\\\\n",
    "1 & 0 & \\dots & 0 & -c_1 \\\\\n",
    "0 & 1 & \\dots & 0 & -c_2 \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
    "0 & 0 & \\dots & 1 & -c_{n-1}\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def myroots(p, verbose=True):\n",
    "    ''' Calculate the roots of a polynomial described by coefficient vector \n",
    "    p(x) = p_0 + p_1 x + p_2 x^2 + \\ldots + p_n x^n    \n",
    "    by finding the eigenvalues of the companion matrix\n",
    "    '''\n",
    "    \n",
    "    # construct the companion matrix of the coefficient vector c\n",
    "    # make p monic and drop the last coefficient\n",
    "    c = p/p[-1]\n",
    "    if verbose:\n",
    "        print(c)\n",
    "    m = len(c) - 1\n",
    "    C = numpy.zeros((m,m))\n",
    "    C[:,-1] = -c[:-1]\n",
    "    C[1:,:-1] = numpy.eye(m-1)\n",
    "    if verbose:\n",
    "        print('C = \\n{}'.format(C))\n",
    "        \n",
    "    return numpy.linalg.eigvals(C)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "c = numpy.array([ -1, -1, 1])\n",
    "r = numpy.roots(c)\n",
    "print(r)\n",
    "mr = myroots(c) \n",
    "print\n",
    "print(mr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c = numpy.random.rand(5)\n",
    "r = numpy.roots(c)\n",
    "print(r)\n",
    "mr = myroots(c) \n",
    "print\n",
    "print(mr)"
   ]
  },
  {
   "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.8.8"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}