EigenValuesAndWhereToFindThem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Eigenvalues and where to find them\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This ipython notebook is mostly a copy of the work done by Marc Khoury (https://marckhoury.github.io/numerical-algorithms-for-computing-eigenvectors/)\n",
    "Most parts of it are direct reproduction, all credit belongs to him.\n",
    "\n",
    "The eigenvalues and eigenvectors of a matrix are essential in many applications across the sciences. Despite their utility, students often leave their linear algebra courses with very little intuition for eigenvectors. In this notebook we describe several surprisingly simple algorithms for computing the eigenvalues and eigenvectors of a matrix, while attempting to convey as much geometric intuition as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The  nice world of positive definite matrices\n",
    "\n",
    "### The  nice world of positive definite matrices\n",
    "\n",
    "#### Basic definition\n",
    "In linear algebra, a symmetric $ n\\times n$ real matrix $M$ is said to be positive definite if the scalar $ z^{T}Mz$ is strictly positive for every non-zero column vector $z \\in \\mathbb{R}^n \\setminus \\vec{0}$. When interpreting $Mz$ as the output of an operator, $M$, acting on an input, $z$, the property of positive definiteness implies that the output always has a positive inner product with the input, as often observed in physical processes.\n",
    "\n",
    "More generally, a complex $n\\times n$ Hermitian matrix $M$ is said to be positive definite if the scalar $z^{*}Mz$ is strictly positive for every non-zero column vector $z \\in \\mathbb{C}^n \\setminus \\vec{0}$.\n",
    "Here $z^{*}$ denotes the conjugate transpose of $z$ and we recall that an hermitian matrix $M$ is such that $M_{i,j} = \\overline{M_{j,i}}$ (that implies that its diagonal is real). Note that $z^{*}Mz$ is automatically real since $M$ is Hermitian.\n",
    "\n",
    "Positive semi-definite matrices are defined similarly, except that the above scalars $z^TMz$ or $z^{*}Mz$ must be positive or zero (i.e. non-negative). Negative definite and negative semi-definite matrices are defined analogously. A matrix that is not positive semi-definite and not negative semi-definite is called indefinite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Alternative definition through bilinear form / inner product\n",
    "\n",
    "The matrix $M$ is positive definite if and only if the bilinear form $\\langle z,w\\rangle =z^{T}Mw$ is positive definite and similarly for a positive definite sesquilinear form in the complex case.\n",
    "\n",
    "##### Bilinear form\n",
    "A bilinear form on a vector space $V$ is a bilinear map $V \\times V \\to K$, where $K$ is the field of scalars. In other words, a bilinear form is a function $B : V \\times V \\to K$ that is linear in each argument separately:\n",
    "\n",
    "\\begin{align*}\n",
    "  & B(u + v, w) = B(u, w) + B(v, w) \\qquad \\text{and} \\qquad B(\\lambda u, v) &= \\lambda B(u, v) \\\\\n",
    "  & B(u, v + w) = B(u, v) + B(u, w) \\qquad \\text{and} \\qquad B(u, \\lambda v) &= \\lambda B(u, v)\n",
    "\\end{align*}\n",
    "\n",
    "##### Sesquilinear form\n",
    "The definition of sesquilinear extends this property to the complex case\n",
    "\n",
    "\\begin{align*}\n",
    "  & \\phi(x+y,z+w) & = \\phi(x,z) + \\phi (x,w) + \\phi(y,z) + \\phi (y,w) \\\\\n",
    "  & \\phi(ax,by) & = \\overline{a} b \\, \\phi(x,y)\n",
    "\\end{align*}\n",
    "\n",
    "##### Inner product\n",
    "\n",
    "A bilinear/sesquilinear form is called an inner product iff it has the following properties:\n",
    "* Conjugate symmetry: $\\langle x,y\\rangle =\\overline{\\langle y,x\\rangle }$\n",
    "* Linearity in the first argument: \n",
    "\\begin{align*}\n",
    "  \\langle ax,y\\rangle &= a \\langle x,y\\rangle \\\\\n",
    "  \\langle x+y,z\\rangle &=\\langle x,z\\rangle +\\langle y,z\\rangle\n",
    "\\end{align*}\n",
    "* Positive-definite: $\\langle x,x\\rangle >0,\\quad x\\in V\\setminus \\{\\mathbf{\\vec{0}} \\}$\n",
    "\n",
    "Inner product can be use to provide a vector space with a lot of useful properties, indeed inner product space is a specific case of more generic structures that are less constrained:\n",
    "![caption](data/250px-Mathematical_Spaces.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eigenvalues and eigenvectors\n",
    "\n",
    "Given a vector space $V$ over the field $F$ and a linear transformation $M : F \\to F$ let's now define what is an \"eigenvalue\".\n",
    "\\begin{align*}\n",
    "  \\lambda \\in F \\text{ is an eigenvalue of } M \\iff $\\exists v \\in F \\setminus \\vec{0} \\text{ s.t } M v = \\lambda v\n",
    "\\end{align*}\n",
    "\n",
    "$v$ is called the eigenvector corresponding the eigenvalue $\\lambda$. There can be multiple eigenvector associated with one eigenvalue, in this case we talk about the multiplicity of this eigenvalue.\n",
    "\n",
    "#### Characteristic polynomial\n",
    "From what we have seen in the previous definition, the eigenvalue $\\lambda$ and the eigenvector $v$ can be seen a solution to the following problem:\n",
    "\n",
    "\\begin{align}\n",
    "  (\\lambda Id - M) v = 0\n",
    "\\end{align}\n",
    "\n",
    "Since $v$ is non-zero, this means that the matrix $\\lambda I − M$ is singular (non-invertible), which in turn means that its determinant is 0. Thus the roots of the function $det(\\lambda I − M)$ are the eigenvalues of M, and it is clear that this determinant is a polynomial in $\\lambda$.\n",
    "This \"trick\" of finding the proper determinant for a matrix to be singular allows to separate the search of eigenvalue from the search of its corresponding eigenvectors.\n",
    "\n",
    "An interesting property of the characteristic polynomial is that, thanks to Cayley–Hamilton theorem, every square matrix over a commutative ring ($\\mathbb{R} or \\mathbb{C} for instance$) satisfies its own characteristic equation.\n",
    "We slightly touched upon this topic in the notebook called \"KrylovBasedMethods\"\n",
    "\n",
    "Unfortunately, finding roots of a polynomial of arbitrary degree can be a hard task. In particular, Abel–Ruffini theorem proves that there are polynomial of degree 5 whose roots cannot be expressed with radicals, and later on, in 1830, Évariste Galois proved that most equations of degree higher than four cannot be solved by radicals, hence no closed form is easily available."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The fundamental theorem of algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eigenvalues of \n",
    "Let's now try to imagine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "Lets recall what posit\n",
    "\n",
    "Let \n",
    "A\n",
    " be a symmetric positive definite matrix. Since \n",
    "A\n",
    " is symmetric all of the eigenvalues of \n",
    "A\n",
    " are real and \n",
    "A\n",
    " has a full set of orthogonal eigenvectors. Let \n",
    "λ\n",
    "1\n",
    "≥\n",
    "λ\n",
    "2\n",
    "≥\n",
    "…\n",
    "≥\n",
    "λ\n",
    "n\n",
    ">\n",
    "0\n",
    "denote the eigenvalues of \n",
    "A\n",
    " and let \n",
    "u\n",
    "1\n",
    ",\n",
    "…\n",
    "u\n",
    "n\n",
    " denote their corresponding eigenvectors. The fact that \n",
    "A\n",
    " is positive definite means that \n",
    "λ\n",
    "i\n",
    ">\n",
    "0\n",
    " for all \n",
    "i\n",
    ". This condition isn’t strictly necessary for the algorithms described below; I’m assuming it so that I can refer to the largest eigenvalue as opposed to the largest in magnitude eigenvalue.\n",
    "\n",
    "All of my intuition for positive definite matrices comes from the geometry of the quadratic form \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    ". Figure 1 plots \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " in \n",
    "R\n",
    "3\n",
    " for several \n",
    "2\n",
    "×\n",
    "2\n",
    " matrices. When \n",
    "A\n",
    " is positive definite, the quadratic form \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " is shaped like a bowl. More rigorously it has positive curvature in every direction and the curvature at the origin in the direction of each eigenvector is proportional to the eigenvalue of that eigenvector. In \n",
    "R\n",
    "3\n",
    ", the two eigenvectors give the directions of the maximum and minimum curvature at the origin. These are also known as principal directions in differential geometry, and the curvatures in these directions are known as principal curvatures. I often shorten this intuition by simply stating that positive definite matrices are bowls, because this is always the picture I have in my head when discussing them.\n",
    "\n",
    "\n",
    "Figure 1: The geometry of the quadratic form \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " for, from left to right, a positive definite matrix, a positive semi-definite matrix, an indefinite matrix, and a negative definite matrix. When \n",
    "A\n",
    " is positive definite it has positive curvature in every direction and is shaped like a bowl. The curvature at the origin in the direction of an eigenvector is proportional to the eigenvalue. A positive semi-definite matrix may have one or more eigenvalues equal to 0. This creates a flat (zero curvature) subspace of dimension equal to the number of eigenvalues with value equal to 0. An indefinite matrix has both positive and negative eigenvalues, and so has some directions with positive curvature and some with negative curvature, creating a saddle. A negative definite matrix has all negative eigenvalues and so the curvature in every direction is negative at every point.\n",
    "Now suppose we wanted to compute a single eigenvector of \n",
    "A\n",
    ". This problem comes up more often than you’d think and it’s a crime that undergraduate linear algebra courses don’t often make this clear. The first algorithm that one generally learns, and the only algorithm in this post that I knew as an undergraduate, is an incredibly simple algorithm called Power Iteration. Starting from a random unit vector \n",
    "v\n",
    " we simply compute \n",
    "A\n",
    "t\n",
    "v\n",
    " iteratively. For sufficiently large \n",
    "t\n",
    ", \n",
    "A\n",
    "t\n",
    "v\n",
    " converges to the eigenvector corresponding to the largest eigenvalue of \n",
    "A\n",
    ", hereafter referred to as the “top eigenvector”.\n",
    "\n",
    "def PowerIteration(A, max_iter):\n",
    "  v = np.random.randn(A.shape[0])\n",
    "  v /= np.linalg.norm(v) #generate a uniformly random unit vector\n",
    "  for t in range(max_iter):\n",
    "    v = np.dot(A, v) #compute Av\n",
    "    v /= np.linalg.norm(v)\n",
    "  return v\n",
    "To see why Power Iteration converges to the top eigenvector of \n",
    "A\n",
    " it helps to write \n",
    "v\n",
    " in the eigenbasis of \n",
    "A\n",
    " as \n",
    "v\n",
    "=\n",
    "∑\n",
    "n\n",
    "i\n",
    "=\n",
    "1\n",
    "β\n",
    "i\n",
    "u\n",
    "i\n",
    "for some coefficients \n",
    "β\n",
    "i\n",
    ". Then we have that\n",
    "\n",
    "A\n",
    "t\n",
    "v\n",
    "=\n",
    "A\n",
    "t\n",
    "(\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "1\n",
    " \n",
    "β\n",
    "i\n",
    "u\n",
    "i\n",
    ")\n",
    "=\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "1\n",
    " \n",
    "β\n",
    "i\n",
    "A\n",
    "t\n",
    "u\n",
    "i\n",
    "=\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "1\n",
    " \n",
    "β\n",
    "i\n",
    "λ\n",
    "t\n",
    "i\n",
    "u\n",
    "i\n",
    "=\n",
    "λ\n",
    "t\n",
    "1\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "1\n",
    " \n",
    "β\n",
    "i\n",
    "(\n",
    "λ\n",
    "i\n",
    "λ\n",
    "1\n",
    ")\n",
    "t\n",
    "u\n",
    "i\n",
    "=\n",
    "λ\n",
    "t\n",
    "1\n",
    "(\n",
    "β\n",
    "1\n",
    "u\n",
    "1\n",
    "+\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "2\n",
    " \n",
    "β\n",
    "i\n",
    "(\n",
    "λ\n",
    "i\n",
    "λ\n",
    "1\n",
    ")\n",
    "t\n",
    "u\n",
    "i\n",
    ")\n",
    ".\n",
    "Since \n",
    "λ\n",
    "1\n",
    " is the largest eigenvalue, the fractions \n",
    "(\n",
    "λ\n",
    "i\n",
    "λ\n",
    "1\n",
    ")\n",
    "t\n",
    " go to 0 as \n",
    "t\n",
    "→\n",
    "∞\n",
    ", for all \n",
    "i\n",
    "≠\n",
    "1\n",
    ". Thus the only component of \n",
    "A\n",
    "t\n",
    "v\n",
    " that has any weight is that of \n",
    "u\n",
    "1\n",
    ". How quickly each of those terms goes to 0 depends on the ratio \n",
    "λ\n",
    "2\n",
    "λ\n",
    "1\n",
    ". If this term is close to 1 then it may take many iterations to disambiguate between the top two (or more) eigenvectors. We say that the Power Iteration algorithm converges at a rate of \n",
    "O\n",
    "(\n",
    "(\n",
    "λ\n",
    "2\n",
    "λ\n",
    "1\n",
    ")\n",
    "t\n",
    ")\n",
    ", which for some unfortunate historical reason is referred to as “linear convergence”.\n",
    "\n",
    "\n",
    "Figure 2: An illustration of the Power Iteration algorithm. The \n",
    "i\n",
    "th bar represents the component of the current iterate on the \n",
    "i\n",
    "th eigenvector, in order of decreasing eigenvalue. Notice that the components corresponding to the smallest eigenvalues decrease most rapidly, whereas the components on the largest eigenvalues take longer to converge. This animation represents 50 iterations of Power Iteration.\n",
    "Power Iteration will give us an estimate of the top eigenvector \n",
    "u\n",
    "1\n",
    ", but what about the other extreme? What if instead we wanted to compute \n",
    "u\n",
    "n\n",
    ", the eigenvector corresponding to the smallest eigenvalue? It turns out there is a simple modification to the standard Power Iteration algorithm that computes \n",
    "u\n",
    "n\n",
    ". Instead of multiplying by \n",
    "A\n",
    " at each iteration, multiply by \n",
    "A\n",
    "−\n",
    "1\n",
    ". This works because the eigenvalues of \n",
    "A\n",
    "−\n",
    "1\n",
    " are \n",
    "1\n",
    "λ\n",
    "i\n",
    ", and thus the smallest eigenvalue of \n",
    "A\n",
    ", \n",
    "λ\n",
    "n\n",
    ", corresponds to the largest eigenvalue of \n",
    "A\n",
    "−\n",
    "1\n",
    ", \n",
    "1\n",
    "λ\n",
    "n\n",
    ". Furthermore the eigenvectors of \n",
    "A\n",
    "−\n",
    "1\n",
    " are unchanged. This slight modification is called Inverse Iteration, and it exhibits the same convergence as Power Iteration, by the same analysis.\n",
    "\n",
    "def InverseIteration(A, max_iter):\n",
    "  v = np.random.randn(A.shape[0])\n",
    "  v /= np.linalg.norm(v) #generate a uniformly random unit vector\n",
    "  lu, piv = scipy.linalg.lu_factor(A) # compute LU factorization of A\n",
    "  for t in range(max_iter):\n",
    "    v = scipy.linalg.lu_solve((lu, piv), v) #compute A^(-1)v\n",
    "    v /= np.linalg.norm(v)\n",
    "  return v\n",
    "Note that we don’t actually compute \n",
    "A\n",
    "−\n",
    "1\n",
    " explicitly. Instead we compute an LU factorization of \n",
    "A\n",
    " and solve the system \n",
    "L\n",
    "U\n",
    "v\n",
    "t\n",
    "+\n",
    "1\n",
    "=\n",
    "v\n",
    "t\n",
    ". The matrix that we’re multiplying by does not change at each iteration, so we can compute the LU factorization once and quickly solve a linear system to compute \n",
    "A\n",
    "−\n",
    "1\n",
    "v\n",
    " at each iteration.\n",
    "\n",
    "\n",
    "Figure 3: The Inverse Iteration algorithm. Notice that in this case the algorithm converges to the eigenvector corresponding to the smallest eigenvalue.\n",
    "Power Iteration and Inverse Iteration find the eigenvectors at the extremes of the spectrum of \n",
    "A\n",
    ", but sometimes we may want to compute a specific eigenvector corresponding to a specific eigenvalue. Suppose that we have an estimate \n",
    "μ\n",
    " of an eigenvalue. We can find the eigenvector corresponding to the eigenvalue of \n",
    "A\n",
    " closest to \n",
    "μ\n",
    " by a simple modification to Inverse Iteration. Instead of multiplying by \n",
    "A\n",
    "−\n",
    "1\n",
    " at each iteration, multiply by \n",
    "(\n",
    "μ\n",
    "I\n",
    "n\n",
    "−\n",
    "A\n",
    ")\n",
    "−\n",
    "1\n",
    " where \n",
    "I\n",
    "n\n",
    " is the identity matrix. The eigenvalues of \n",
    "(\n",
    "μ\n",
    "I\n",
    "n\n",
    "−\n",
    "A\n",
    ")\n",
    "−\n",
    "1\n",
    " are \n",
    "1\n",
    "μ\n",
    "−\n",
    "λ\n",
    "i\n",
    ". Thus the largest eigenvalue of \n",
    "(\n",
    "μ\n",
    "I\n",
    "n\n",
    "−\n",
    "A\n",
    ")\n",
    "−\n",
    "1\n",
    " corresponds to the eigenvalue of \n",
    "A\n",
    " whose value is closest to \n",
    "μ\n",
    ". By the same analysis as Power Iteration, Shifted Inverse Iteration also exhibits linear convergence. However the better the estimate \n",
    "μ\n",
    " the larger \n",
    "1\n",
    "μ\n",
    "−\n",
    "λ\n",
    "i\n",
    " and, consequently, the faster the convergence.\n",
    "\n",
    "def ShiftedInverseIteration(A, mu, max_iter):\n",
    "  I = np.identity(A.shape[0])\n",
    "  v = np.random.randn(A.shape[0])\n",
    "  v /= np.linalg.norm(v) #generate a uniformly random unit vector\n",
    "  lu, piv = scipy.linalg.lu_factor(mu*I - A) # compute LU factorization of (mu*I - A)\n",
    "  for t in range(max_iter):\n",
    "    v = scipy.linalg.lu_solve((lu, piv), v) #compute (mu*I - A)^(-1)v\n",
    "    v /= np.linalg.norm(v)\n",
    "  return v\n",
    "\n",
    "Figure 4: The Shifted Inverse Iteration algorithm. In this case we converge to the eigenvector corresponding to the eigenvalue nearest \n",
    "μ\n",
    ".\n",
    "Shifted Inverse Iteration converges quickly if a good estimate of the target eigenvalue is available. However if \n",
    "μ\n",
    " is a poor approximation of the desired eigenvalue, Shifted Inverse Iteration may take a long time to converge. In fact all of the algorithms we’ve presented so far have exactly the same convergence rate; they all converge linearly. If instead we could improve on the eigenvalue estimate at each iteration we could potentially develop an algorithm with a faster convergence rate. This is the main idea behind Rayleigh Quotient Iteration.\n",
    "\n",
    "The Rayleigh quotient is defined as \n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "=\n",
    "v\n",
    "⊤\n",
    "A\n",
    "v\n",
    "v\n",
    "⊤\n",
    "v\n",
    " for any vector \n",
    "v\n",
    ". There are many different ways in which we can understand the Rayleigh quotient. Some intuition that is often given is that the Rayleigh quotient is the scalar value that behaves most like an “eigenvalue” for \n",
    "v\n",
    ", even though \n",
    "v\n",
    " may not be an eigenvector. What is meant is that the Rayleigh quotient is the minimum to the optimization problem \n",
    "min\n",
    "λ\n",
    "∈\n",
    "R\n",
    "|\n",
    "|\n",
    "A\n",
    "v\n",
    "−\n",
    "λ\n",
    "v\n",
    "|\n",
    "|\n",
    "2\n",
    ". This intuition is hardly satisfying.\n",
    "\n",
    "Let’s return to the geometry of the quadratic forms \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " and \n",
    "x\n",
    "⊤\n",
    "x\n",
    " which comprise the Rayleigh quotient, drawn in orange and blue respectively in Figure 5. Without loss of generality we can assume that \n",
    "A\n",
    " is a diagonal matrix. (This is without loss of generality because we’re merely rotating the surface so that the eigenvectors align with the \n",
    "x\n",
    " and \n",
    "y\n",
    " axes, which does not affect the geometry of the surface. This is a common trick in the numerical algorithms literature.) In this coordinate system, the quadratic form \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    "=\n",
    "λ\n",
    "1\n",
    "x\n",
    "2\n",
    "1\n",
    "+\n",
    "λ\n",
    "2\n",
    "x\n",
    "2\n",
    "2\n",
    ", where \n",
    "λ\n",
    "1\n",
    " and \n",
    "λ\n",
    "2\n",
    " are the diagonal entries, and thus the eigenvalues, of \n",
    "A\n",
    ".\n",
    "\n",
    "Consider any vector \n",
    "v\n",
    " and let \n",
    "h\n",
    "=\n",
    "span\n",
    "{\n",
    "v\n",
    ",\n",
    "(\n",
    "0\n",
    ",\n",
    "0\n",
    ",\n",
    "1\n",
    ")\n",
    "}\n",
    " be the plane spanned by \n",
    "v\n",
    " and the vector \n",
    "(\n",
    "0\n",
    ",\n",
    "0\n",
    ",\n",
    "1\n",
    ")\n",
    ". The intersection of \n",
    "h\n",
    " with the quadratic forms \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " and \n",
    "x\n",
    "⊤\n",
    "x\n",
    " is comprised of two parabolas, also shown in Figure 5. (This is a common trick in the geometric algorithms literature.) If \n",
    "v\n",
    " is aligned with the \n",
    "x\n",
    "-axis, then, within the coordinate system defined by \n",
    "h\n",
    ", \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " can be parameterized by \n",
    "y\n",
    "=\n",
    "λ\n",
    "1\n",
    "x\n",
    "2\n",
    " and \n",
    "x\n",
    "⊤\n",
    "x\n",
    " can be parameterized by \n",
    "y\n",
    "=\n",
    "x\n",
    "2\n",
    ". (Note that here \n",
    "y\n",
    " and \n",
    "x\n",
    " refer to local coordinates within \n",
    "h\n",
    " and are distinct from the vector \n",
    "x\n",
    " used in \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    ".) Similarly if \n",
    "v\n",
    "is aligned with the \n",
    "y\n",
    "-axis, then \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " can be parameterized by \n",
    "y\n",
    "=\n",
    "λ\n",
    "2\n",
    "x\n",
    "2\n",
    ". (If \n",
    "v\n",
    " is any other vector then \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " can be parameterized by \n",
    "y\n",
    "=\n",
    "κ\n",
    "x\n",
    "2\n",
    " for some \n",
    "κ\n",
    " dependent upon \n",
    "v\n",
    ".) The Rayleigh quotient at \n",
    "v\n",
    " is \n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "=\n",
    "λ\n",
    "1\n",
    "x\n",
    "2\n",
    "x\n",
    "2\n",
    "=\n",
    "λ\n",
    "1\n",
    ". The curvature of the parabola \n",
    "y\n",
    "=\n",
    "λ\n",
    "1\n",
    "x\n",
    "2\n",
    " at the origin is \n",
    "2\n",
    "λ\n",
    "1\n",
    ". Thus the Rayleigh quotient is proportional to the the curvature of \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " in the direction \n",
    "v\n",
    "!\n",
    "\n",
    "\n",
    "Figure 5: The quadratic form \n",
    "x\n",
    "⊤\n",
    "A\n",
    "x\n",
    " is shown in orange and \n",
    "x\n",
    "⊤\n",
    "x\n",
    " is shown in blue. Intersecting both surfaces with a plane \n",
    "h\n",
    " gives two parabola. Within the plane \n",
    "h\n",
    " we can define a local coordinate system and parameterize both parabola as \n",
    "κ\n",
    "x\n",
    "2\n",
    " and \n",
    "x\n",
    "2\n",
    ". The Rayleigh quotient is equal to the ratio of the heights of the parabolas at any point, which is always equal to \n",
    "κ\n",
    ".\n",
    "From this intuition it is clear that the value of the Rayleigh quotient is identical along any ray starting at, but not including, the origin. The length of \n",
    "v\n",
    " corresponds to the value of \n",
    "x\n",
    " in the coordinate system defined by \n",
    "h\n",
    ", which does not affect the Rayleigh quotient. We can also see this algebraically, by choosing a unit vector \n",
    "v\n",
    " and parameterizing a ray in the direction \n",
    "v\n",
    " as \n",
    "α\n",
    "v\n",
    " for \n",
    "α\n",
    "∈\n",
    "R\n",
    " and \n",
    "α\n",
    ">\n",
    "0\n",
    ". Then we have that\n",
    "\n",
    "λ\n",
    "R\n",
    "(\n",
    "α\n",
    "v\n",
    ")\n",
    "=\n",
    "(\n",
    "α\n",
    "v\n",
    "⊤\n",
    ")\n",
    "A\n",
    "(\n",
    "α\n",
    "v\n",
    ")\n",
    "α\n",
    "2\n",
    "v\n",
    "⊤\n",
    "v\n",
    "=\n",
    "v\n",
    "⊤\n",
    "A\n",
    "v\n",
    "v\n",
    "⊤\n",
    "v\n",
    "=\n",
    "v\n",
    "⊤\n",
    "A\n",
    "v\n",
    ".\n",
    "Thus it is sufficient to consider the values of the Rayleigh quotient on the unit sphere.\n",
    "\n",
    "For a unit vector \n",
    "v\n",
    " the value of the Rayleigh quotient can be written in the eigenbasis as \n",
    "v\n",
    "⊤\n",
    "A\n",
    "v\n",
    "=\n",
    "n\n",
    "∑\n",
    "i\n",
    "=\n",
    "1\n",
    " \n",
    "λ\n",
    "i\n",
    "⟨\n",
    "v\n",
    ",\n",
    "u\n",
    "i\n",
    "⟩\n",
    "2\n",
    " where \n",
    "∑\n",
    "n\n",
    "i\n",
    "=\n",
    "1\n",
    "⟨\n",
    "v\n",
    ",\n",
    "u\n",
    "i\n",
    "⟩\n",
    "2\n",
    "=\n",
    "1\n",
    ". Thus the Rayleigh quotient is a convex combination of the eigenvalues of \n",
    "A\n",
    " and so its value is bounded by the minimum and maximum eigenvalues \n",
    "λ\n",
    "n\n",
    "≤\n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "≤\n",
    "λ\n",
    "1\n",
    "for all \n",
    "v\n",
    ". This fact is also easily seen from the geometric picture above, as the curvature at the origin is bounded by twice the minimum and maximum eigenvalues. It can be readily seen by either direct calculation or by the coefficients of the convex combination, that if \n",
    "v\n",
    " is an eigenvector, then \n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    " is the corresponding eigenvalue of \n",
    "v\n",
    ".\n",
    "\n",
    "Recall that a critical point of a function is a point where the derivative is equal to 0. It should come as no surprise that the eigenvalues are the critical values of the Rayleigh quotient and the eigenvectors are the critical points. What is less obvious is the special geometric structure of the critical points.\n",
    "\n",
    "The gradient of the Rayleigh quotient is \n",
    "2\n",
    "v\n",
    "⊤\n",
    "v\n",
    "(\n",
    "A\n",
    "v\n",
    "−\n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "v\n",
    ")\n",
    ", from which it is easy to see that every eigenvector is a critical point of \n",
    "λ\n",
    "R\n",
    ". The type of critical point is determined by the Hessian matrix, which at the critical point \n",
    "u\n",
    "i\n",
    " is \n",
    "2\n",
    "(\n",
    "A\n",
    "−\n",
    "λ\n",
    "i\n",
    "I\n",
    ")\n",
    ". The eigenvalues of the Hessian are \n",
    "λ\n",
    "j\n",
    "−\n",
    "λ\n",
    "i\n",
    " for \n",
    "j\n",
    "∈\n",
    "[\n",
    "1\n",
    ",\n",
    "n\n",
    "]\n",
    ". Assuming for a moment that the eigenvalues are all distinct, the matrix \n",
    "2\n",
    "(\n",
    "A\n",
    "−\n",
    "λ\n",
    "i\n",
    "I\n",
    ")\n",
    " has \n",
    "i\n",
    "−\n",
    "1\n",
    " eigenvectors that are positive, one eigenvalue that is 0, and \n",
    "n\n",
    "−\n",
    "i\n",
    " eigenvalues that are negative. The 0 eigenvalue represents the fact that the value of the Rayleigh quotient is unchanged along the ray \n",
    "α\n",
    "u\n",
    "i\n",
    ". The other eigenvalues represent the fact that at \n",
    "u\n",
    "i\n",
    ", along the unit sphere, there are \n",
    "i\n",
    "−\n",
    "1\n",
    " directions in which we can walk to increase the value of the Rayleigh quotient, and \n",
    "n\n",
    "−\n",
    "i\n",
    " directions that decrease the Rayleigh quotient. Thus each eigenvector gives rise to a different type of saddle, and there are exactly two critical points of each type on the unit sphere.\n",
    "\n",
    "\n",
    "Figure 6: Contours of the Rayleigh quotient on the unit sphere and the gradient of the Rayleigh quotient at each point. We clearly see one minimum in blue corresponding to the minimum eigenvalue, one saddle point, and one maximum in bright yellow corresponding to the maximum eigenvalue.\n",
    "Finally we come to the crown jewel of the algorithms in this post. The Rayleigh Quotient Iteration algorithm simply updates the estimate \n",
    "μ\n",
    " at each iteration with the Rayleigh quotient. Other than this slight modification, the algorithm is exactly like Shifted Inverse iteration.\n",
    "\n",
    "def RayleighQuotientIteration(A, max_iter):\n",
    "  I = np.identity(A.shape[0])\n",
    "  v = np.random.randn(A.shape[0])\n",
    "  v /= np.linalg.norm(v) #generate a uniformly random unit vector\n",
    "  mu = np.dot(v, np.dot(A, v))\n",
    "  for t in range(max_iter):\n",
    "    v = np.linalg.solve(mu * I - A, v) #compute (mu*I - A)^(-1)v\n",
    "    v /= np.linalg.norm(v)\n",
    "    mu = np.dot(v, np.dot(A, v)) #compute Rayleigh quotient\n",
    "  return (v, mu)\n",
    "This slight modification drastically improves the convergence rate. Unlike the other algorithms in this post which converge linearly, Rayleigh quotient iteration exhibits local cubic convergence! This means that, assuming \n",
    "∥\n",
    "v\n",
    "t\n",
    "−\n",
    "u\n",
    "i\n",
    "∥\n",
    "≤\n",
    "ϵ\n",
    " for some \n",
    "u\n",
    "i\n",
    ", on the next iteration we will have that \n",
    "∥\n",
    "v\n",
    "t\n",
    "+\n",
    "1\n",
    "−\n",
    "u\n",
    "i\n",
    "∥\n",
    "≤\n",
    "ϵ\n",
    "3\n",
    ". In practice this means that you should expect triple the number of correct digits at each iteration. It’s hard to understate how crazy fast cubic convergence is, and, to the best of the author’s knowledge, algorithms that exhibit cubic convergence are rare in the numerical algorithms literature.\n",
    "\n",
    "\n",
    "Figure 7: The Rayleigh Quotient Iteration algorithm. After only 6 iterations the eigenvalue estimate \n",
    "μ\n",
    "t\n",
    " is so accurate that the resulting matrix \n",
    "(\n",
    "μ\n",
    "t\n",
    "I\n",
    "n\n",
    "−\n",
    "A\n",
    ")\n",
    " is singular up-to machine precision and we can no longer solve the system for an inverse. Note that every other figure in this post shows 50 iterations.\n",
    "Intuitively, the reason that Rayleigh Quotient Iteration exhibits cubic convergence is because, while the Shifted Inverse Iteration step converges linearly, the Rayleigh quotient is a quadratically good estimate of an eigenvalue near an eigenvector. To see this consider the Taylor series expansion of \n",
    "λ\n",
    "R\n",
    " near an eigenvector \n",
    "u\n",
    "i\n",
    ".\n",
    "\n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "=\n",
    "λ\n",
    "R\n",
    "(\n",
    "u\n",
    "i\n",
    ")\n",
    "+\n",
    "(\n",
    "v\n",
    "−\n",
    "u\n",
    "i\n",
    ")\n",
    "⊤\n",
    "∇\n",
    "λ\n",
    "R\n",
    "(\n",
    "u\n",
    "i\n",
    ")\n",
    "+\n",
    "O\n",
    "(\n",
    "|\n",
    "|\n",
    "v\n",
    "−\n",
    "u\n",
    "i\n",
    "|\n",
    "|\n",
    "2\n",
    ")\n",
    "=\n",
    "λ\n",
    "R\n",
    "(\n",
    "u\n",
    "i\n",
    ")\n",
    "+\n",
    "O\n",
    "(\n",
    "|\n",
    "|\n",
    "v\n",
    "−\n",
    "u\n",
    "i\n",
    "|\n",
    "|\n",
    "2\n",
    ")\n",
    "λ\n",
    "R\n",
    "(\n",
    "v\n",
    ")\n",
    "−\n",
    "λ\n",
    "R\n",
    "(\n",
    "u\n",
    "i\n",
    ")\n",
    "=\n",
    "O\n",
    "(\n",
    "|\n",
    "|\n",
    "v\n",
    "−\n",
    "u\n",
    "i\n",
    "|\n",
    "|\n",
    "2\n",
    ")\n",
    "The second step follows from the fact that \n",
    "u\n",
    "i\n",
    " is a critical point of \n",
    "λ\n",
    "R\n",
    " and so \n",
    "∇\n",
    "λ\n",
    "R\n",
    "(\n",
    "u\n",
    "i\n",
    ")\n",
    "=\n",
    "0\n",
    ".\n",
    "\n",
    "While Rayleigh Quotient Iteration exhibits very fast convergence, it’s not without its drawbacks. First, notice that the system \n",
    "(\n",
    "μ\n",
    "t\n",
    "I\n",
    "−\n",
    "A\n",
    ")\n",
    "−\n",
    "1\n",
    " changes at each iteration. Thus we cannot precompute a factorization of this matrix and quickly solve the system using forward and backward substitution at each iteration, like we did in the Shifted Inverse Iteration algorithm. We need to solve a different linear system at each iteration, which is much more expensive. Second, Rayleigh Quotient Iteration gives no control over to which eigenvector it converges. The eigenvector it converges to depends on which basin of attraction the initial random vector \n",
    "v\n",
    "0\n",
    " falls into. Thus cubic convergence comes at a steep cost. This balance between an improved convergence rate and solving a different linear system at each iteration feels like mathematical poetic justice. The price to pay for cubic convergence is steep.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}