11_LA_QR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": [
    "# Orthogonalization:\n",
    "# The QR Factorization, Projections,Least Squares Problems and Applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Projections\n",
    "\n",
    "A **projector** is a square matrix $P$ that satisfies\n",
    "$$\n",
    "    P^2 = P.\n",
    "$$\n",
    "\n",
    "Why does this definition make sense?  Why do we require it to be square?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A projector comes from the idea that we want to project a vector $\\mathbf{v}$ onto a lower dimensional subspace.  For example, suppose that $\\mathbf{v}$ lies completely within the subspace, i.e. $\\mathbf{v} \\in \\text{range}(P)$. If that is the case then $P \\mathbf{v}$ should not change, or $P\\mathbf{v} = \\mathbf{v}$.  This motivates the definition above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "i.e.  if\n",
    "$$\n",
    "    P\\mathbf{v} = \\mathbf{v} \n",
    "$$\n",
    "then\n",
    "$$\n",
    "    P( P \\mathbf{v} ) = P\\mathbf{v} = \\mathbf{v}.\n",
    "$$\n",
    "or $$P^2 = P$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As another example, take a vector $\\mathbf{x} \\notin \\text{range}(P)$ and project it onto the subspace $P\\mathbf{x} = \\mathbf{v}$.  If we apply the projection again to $\\mathbf{v}$ we now have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    P\\mathbf{x} &= \\mathbf{v}, \\\\\n",
    "    P^2 \\mathbf{x} & = P \\mathbf{v} = \\mathbf{v} \\\\\n",
    "    \\Rightarrow P^2 &= P.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It is also important to keep in mind the following, given again $\\mathbf{x} \\notin \\text{range}(P)$, if we look at the difference between the projection and the original vector $P\\mathbf{x} - \\mathbf{x}$ and apply the projection again we have \n",
    "\n",
    "$$\n",
    "    P(P\\mathbf{x} - \\mathbf{x}) = P^2 \\mathbf{x} - P\\mathbf{x} = 0\n",
    "$$\n",
    "\n",
    "which means the difference between $\\mathbf{x}$ and the projected vector $P\\mathbf{x} = \\mathbf{v}$ lies in the null space of $P$, i.e.\n",
    "\n",
    "$$\n",
    "(\\mathbf{x}-\\mathbf{v}) \\in \\text{null}(P).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "But\n",
    "$$\n",
    "    \\mathbf{x} - \\mathbf{v} = \\mathbf{x} - P\\mathbf{x} = (I - P)\\mathbf{x}\n",
    "$$\n",
    "\n",
    "Thus $I-P$ is a matrix that projects onto $\\mathrm{null}(P)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Complementary Projectors\n",
    "\n",
    "A projector also has a complement defined to be $I - P$.  \n",
    "\n",
    "Show that this complement is also a projector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can show that this a projector by examining a repeated application of $(I-P)$:\n",
    "$$\\begin{aligned}\n",
    "    (I - P)^2 &= I - IP - IP + P^2 \\\\\n",
    "    &= I - 2 P + P^2 \\\\\n",
    "    &= I - 2P + P \\\\\n",
    "    &= I - P.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out that the complement projects exactly onto $\\text{null}(P)$.  \n",
    "\n",
    "Take \n",
    "$$\n",
    "    \\mathbf{x} \\in \\text{null}(P),\n",
    "$$ \n",
    "\n",
    "then \n",
    "$$\n",
    "    (I - P) \\mathbf{x} = \\mathbf{x} - P \\mathbf{x} = \\mathbf{x}\n",
    "$$ \n",
    "\n",
    "since $P \\mathbf{x} = 0$ implying that $\\mathbf{x} \\in \\text{range}(I - P)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We also know that \n",
    "$$\n",
    "    (I - P) \\mathbf{x }\\in \\text{null}(P)\n",
    "$$ \n",
    "as well.  \n",
    "\n",
    "This shows that the \n",
    "$$\n",
    "    \\text{range}(I - P) \\subseteq \\text{null}(P)\n",
    "$$ \n",
    "and \n",
    "$$\n",
    "    \\text{range}(I - P) \\supseteq \\text{null}(P)\n",
    "$$ \n",
    "implying that \n",
    "$$\n",
    "    \\text{range}(I - P) = \\text{null}(P)\n",
    "$$ \n",
    "exactly.  \n",
    "\n",
    "Reflect on these subspaces and convince yourself that this all makes sense."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This result provides an important property of a projector and its complement, namely that they divide a space into two subspaces whose intersection is \n",
    "\n",
    "$$\n",
    "    \\text{range}(I - P) \\cap \\text{range}(P) = \\{0\\}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "    \\text{null}(P) \\cap \\text{range}(P) = \\{0\\}\n",
    "$$\n",
    "\n",
    "These two spaces are called **complementary subspaces**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given this property we can take any $P \\in \\mathbb C^{m \\times m}$ which will split $\\mathbb C^{m \\times m}$ into two subspaces $S$ and $V$, assume that $\\mathbf{s}\\in S = \\text{range}(P)$, and $\\mathbf{v} \\in V = \\text{null}(P)$.  If we have $\\mathbf{x} \\in \\mathbb C^{m \\times m}$ that we can split the vector $\\mathbf{x}$ into components in $S$ and $V$ by using the projections\n",
    "$$\\begin{aligned}\n",
    "    P \\mathbf{x} = \\mathbf{x}_S& &\\mathbf{x}_s \\in S \\\\\n",
    "    (I - P) \\mathbf{x} = \\mathbf{x}_V& &\\mathbf{x}_V \\in V\n",
    "\\end{aligned}$$\n",
    "which we can also observe adds to the original vector as\n",
    "$$\n",
    "    \\mathbf{x}_S + \\mathbf{x}_V = P \\mathbf{x} + (I - P) \\mathbf{x} = \\mathbf{x}.\n",
    "$$\n",
    "\n",
    "Try constructing a projection matrix so that $P \\in \\mathbb R^3$ that projects a vector into one of the coordinate directions ($\\mathbb R$).  \n",
    " - What is the complementary projector?\n",
    " - What is the complementary subspace?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Example:  A non-orthogonal non-linear projector\n",
    "\n",
    "Given a vector of mols of $N$ chemical components\n",
    "\n",
    "$$\n",
    "\\mathbf{n} = \\begin{bmatrix} n_1 \\\\ n_2 \\\\ \\vdots \\\\ n_N\\end{bmatrix}\n",
    "$$\n",
    "where (e.g. $n_1$ is the number of moles of component $1$) and $n_i\\geq 0$\n",
    "\n",
    "Then we can define the mol fraction  of component $i$ as\n",
    "$$\n",
    "    x_i = \\frac{n_i}{\\mathbf{n}^T\\mathbf{1}}\n",
    "$$  \n",
    "and the \"vector\" of mole fractions \n",
    "$$\n",
    "\\mathbf{x} = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_N\\end{bmatrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the homework you will show that\n",
    "\n",
    "* $\\mathbf{x}^T\\mathbf{1} = 1$ (the sum of the mole fractions add to 1)\n",
    "* mole fractions do not form a vector space or subspace\n",
    "* There exists a non Orthogonal projector $f$ such that $f(\\mathbf{n})=\\mathbf{x}$, $f^2=f$\n",
    "* $P$ is singular (like all projection matrices) such that if you know the mole-fractions you don't know how many moles you have.\n",
    "* Find $N(P)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Orthogonal Projectors\n",
    "\n",
    "An **orthogonal projector** is one that projects onto a subspace $S$ that is orthogonal to the complementary subspace $V$ (this is also phrased that $S$ projects along a space $V$).  Note that we are only talking about the subspaces (and their basis), not the projectors!  i.e. for orthogonal **subspaces** all the vectors in $S$ are orthogonal to all the vectors in $V$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "A **hermitian** matrix is one whose complex conjugate transposed is itself, i.e.\n",
    "\n",
    "$$\n",
    "    P = P^\\ast.\n",
    "$$\n",
    "\n",
    "With this definition we can then say:  *A projector $P$ is orthogonal if and only if $P$ is hermitian.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Quick Proof**\n",
    "\n",
    "Show that if $P^2 = P$ and $P^\\ast = P$, then \n",
    "$$\n",
    "\\langle P\\mathbf{x}, (I-P)\\mathbf{x}\\rangle=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Projection with an Orthonormal Basis\n",
    "\n",
    "We can also directly construct a projector that uses an orthonormal basis on the subspace $S$.  If we define another matrix $Q \\in \\mathbb C^{m \\times n}$ which is unitary (its columns are orthonormal) we can construct an orthogonal projector as\n",
    "$$\n",
    "    P = Q Q^*.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Check that $P^* = P$ and $P^2 = P$  (Remember $Q^*Q=I$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Note that the resulting matrix $P$ is in $\\mathbb C^{m \\times m}$ as we require.  This means also that the dimension of the subspace $S$ is $n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Orthonormal projection and Least-Squares Problems... A review\n",
    "\n",
    "Consider the overdetermined problem $A\\mathbf{x}=\\mathbf{b}$ where $A\\in\\mathbb{R}^{3\\times2}$ and $\\mathbf{b}\\in\\mathbb{R}^3$ i.e. \n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix} | & | \\\\\n",
    "                    \\mathbf{a}_1 & \\mathbf{a}_2\\\\\n",
    "                     | & | \\\\\n",
    "     \\end{bmatrix} \\begin{bmatrix} x_1 \\\\ x_2\\\\ \\end{bmatrix} \n",
    "         = \\begin{bmatrix} |   \\\\\n",
    "                    \\mathbf{b} \\\\\n",
    "                     |  \\\\ \n",
    "                     \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and $\\mathbf{a}_1$, $\\mathbf{a}_2$ are linearly independent vectors that span a two-dimensional subspace of $\\mathbb{R}^3$.                   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Geometry\n",
    "\n",
    "Geometrically this problem looks like\n",
    "\n",
    "<img align=center, src=\"./images/least_squares_geometry.jpg\" alt=\"Drawing\" width=600/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If $\\mathbf{b}\\notin C(A)$, then there is clearly no solution to $A\\mathbf{x}=\\mathbf{b}$.  However, we can find the point $\\mathbf{p}=A\\hat{\\mathbf{x}}\\in C(A)$ that minimizes the length of the the error $\\mathbf{e}=\\mathbf{b}-\\mathbf{p}$.  While we could resort to calculus to find the values of $\\hat{\\mathbf{x}}$ that minimizes $||\\mathbf{e}||_2$.  It should be clear from the figure that the shortest error (in the $\\ell_2$ norm) is the one that is perpendicular to every vector in $C(A)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But the sub-space of vectors orthogonal to $C(A)$ is the left-Null Space $N(A^T)$,  and therefore we simply seek solutions of \n",
    "\n",
    "\\begin{align}\n",
    "    0 &= A^T\\mathbf{e} \\\\\n",
    "    &= A^T(\\mathbf{b}-\\mathbf{p})\\\\\n",
    "    &= A^T(\\mathbf{b} - A\\hat{\\mathbf{x}})\\\\\n",
    "\\end{align}\n",
    "\n",
    "or we just need to solve the \"Normal Equations\"  $A^T A\\hat{\\mathbf{x}} = A^T\\mathbf{b}$ for the least-squares solution\n",
    "$$\n",
    "    \\hat{\\mathbf{x}} = \\left(A^T A\\right)^{-1}A^T\\mathbf{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "if we're actually interested in $\\mathbf{p}$ which is the orthogonal projection of $\\mathbf{b}$ onto $C(A)$ we get\n",
    "\n",
    "$$\n",
    "    \\mathbf{p}= A\\hat{\\mathbf{x}} = A \\left(A^T A\\right)^{-1}A^T\\mathbf{b} = P\\mathbf{b}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "    P = A \\left(A^T A\\right)^{-1}A^T\n",
    "$$ \n",
    "\n",
    "is an orthogonal projection matrix (verify that $P^2 = P$ and $(I - P)\\mathbf{b}\\perp P\\mathbf{b})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For a general matrix $A$,  this form of the projection matrix is rather horrid to find,  however,  if the columns of $A$ formed an orthonormal basis for $C(A)$,  i.e. $A=Q$,  then the form of the projection matrix is much simpler as $Q^T Q=I$, therefore\n",
    "\n",
    "$$\n",
    "    P = QQ^T\n",
    "$$\n",
    "\n",
    "This is actually quite general. Given any orthonormal basis for a vector space $S=\\mathrm{span}\\langle \\mathbf{q}_1,\\mathbf{q}_2,\\ldots,\\mathbf{q}_N\\rangle$.  If these vectors form the columns of $Q$,  then the orthogonal projector onto $S$ is always $QQ^T$ and the complement is always $I-QQ^T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Example:  Construction of an orthonormal projector**\n",
    "\n",
    "Take $\\mathbb R^3$ and derive a projector that projects onto the x-y plane and is an orthogonal projector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "the simplest Orthonormal basis for the $x-y$ plane are the columns of \n",
    "$$\n",
    "    Q = \\begin{bmatrix} 1 & 0 \\\\ \n",
    "                        0 & 1 \\\\ \n",
    "                        0 & 0 \n",
    "        \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and an orthogonal projector onto the $x-y$ plane is simply\n",
    "\n",
    "$$ \n",
    "    Q Q^\\ast = \\begin{bmatrix} 1 & 0 \\\\ \n",
    "                               0 & 1 \\\\ \n",
    "                               0 & 0 \n",
    "               \\end{bmatrix} \n",
    "               \\begin{bmatrix} 1 & 0 & 0 \\\\\n",
    "                               0 & 1 & 0 \n",
    "               \\end{bmatrix} = \\begin{bmatrix} \n",
    "                               1 & 0 & 0 \\\\\n",
    "                               0 & 1 & 0 \\\\ \n",
    "                               0 & 0 & 0 \n",
    "               \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Q = numpy.array([[1., 0.],[0., 1.],[0., 0.]])\n",
    "P = numpy.dot(Q, Q.T)\n",
    "I = numpy.identity(3)\n",
    "\n",
    "x = numpy.array([3., 4., 5.])\n",
    "\n",
    "print('x = {}'.format(x))\n",
    "print('\\nQ = \\n{}'.format(Q))\n",
    "print('\\nP = \\n{}'.format(P))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x_S = numpy.dot(P, x)\n",
    "x_V = numpy.dot(I - P, x)\n",
    "print('x_S = {}'.format(x_S))\n",
    "print('x_V = {}\\n'.format(x_V))\n",
    "print('x_S + x_V = {}'.format(x_S+x_V))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A numerically more sensible approach\n",
    "\n",
    "The previous problem calculated the projection matrix first, $P=QQ^T$ then calculated the projection and its complement by matrix vector multiplication, i.e.\n",
    "\n",
    "$$\n",
    "    \\mathbf{x}_S = P\\mathbf{x} = (QQ^T)\\mathbf{x}, \\quad \\mathbf{x}_V = (I - P)\\mathbf{x}\n",
    "$$\n",
    "\n",
    "A mathematically equivalent, but numerically more efficient method is to calculate the following\n",
    "$$\n",
    "    \\mathbf{x}_S = Q(Q^T\\mathbf{x}),\\quad \\mathbf{x}_V = \\mathbf{x} - \\mathbf{x}_S\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "How much more efficient is the latter than the former?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Check "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Q = numpy.array([[1, 0],[0, 1],[0, 0]])\n",
    "x = numpy.array([3., 4., 5.])\n",
    "\n",
    "print('x = {}'.format(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x_S = Q.dot(Q.T.dot(x))\n",
    "x_V = x - x_S\n",
    "print('x_S = {}'.format(x_S))\n",
    "print('x_V = {}\\n'.format(x_V))\n",
    "print('x_S + x_V = {}'.format(x_S+x_V))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: Construction of a projector that eliminates a direction\n",
    "\n",
    "Goal:  Eliminate the component of a vector in the direction $\\mathbf{q}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "<img align=center, src=\"./images/least_squares_geometry.jpg\" alt=\"Drawing\" width=600/>\n",
    "\n",
    "e.g. We can calculate the projection $\\mathbf{p}$ in two equivalent ways\n",
    "\n",
    "1. find the matrix $P=QQ^T$ that projects $\\mathbf{b}$ onto $C(A)$\n",
    "2. find the matrix $P'=\\mathbf{q}\\mathbf{q}^T$ that projects onto the unit vector $\\mathbf{q}$ normal to the plane (i.e. parallel to $\\mathbf{e}$), and then take its Complement projection $\\mathbf{p} = (I - P')\\mathbf{b} = \\mathbf{b} - \\mathbf{q}(\\mathbf{q}^T\\mathbf{b})$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Form the projector $P = \\mathbf{q} \\mathbf{q}^\\ast \\in \\mathbb C^{m \\times m}$.  The complement $I - P$ will then include everything **BUT** that direction. \n",
    "\n",
    "If $||\\mathbf{q}|| = 1$ we can then simply use $I - \\mathbf{q} \\mathbf{q}^\\ast$.  If not we can write the projector in terms of the arbitrary vector $\\mathbf{a}$ as\n",
    "$$\n",
    "    I - \\frac{\\mathbf{a} \\mathbf{a}^\\ast}{||\\mathbf{a}||^2} = I - \\frac{\\mathbf{a} \\mathbf{a}^\\ast}{\\mathbf{a}^\\ast \\mathbf{a}}\n",
    "$$\n",
    "Note that differences in the resulting dimensions between the two values in the fraction.  Also note that as we saw with the outer product, the resulting $\\text{rank}(\\mathbf{a} \\mathbf{a}^\\ast) = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now again try to construct a projector in $\\mathbb R^3$ that projects onto the $x$-$y$ plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "q = numpy.array([0, 0, 1])\n",
    "P = numpy.outer(q, q.conjugate())\n",
    "P_comp = numpy.identity(3) - P\n",
    "print('P = \\n{}\\n'.format(P))\n",
    "print('I - P = \\n{}'.format(P_comp))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.array([3, 4, 5])\n",
    "print(numpy.dot(P, x), q*(q.dot(x)))\n",
    "print(numpy.dot(P_comp, x), x - q*(q.dot(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "a = numpy.array([0, 0, 3])\n",
    "P = numpy.outer(a, a.conjugate()) / (numpy.dot(a, a.conjugate()))\n",
    "P_comp = numpy.identity(3) - P\n",
    "print(numpy.dot(P, x))\n",
    "print(numpy.dot(P_comp, x)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quick Review\n",
    "\n",
    "* A projection matrix is any square matrix $P$ such that $P^2=P$\n",
    "* $P$ projects onto a subspace $S=\\mathrm{range}(P)$\n",
    "* The complementary projection matrix $I-P$ projects onto $V=\\mathrm{null}(P)$\n",
    "* An *orthogonal* projector can always be constructed as $P=QQ^T$ where the columns of $Q$ form an *orthonormal* basis for $S$ and $S\\perp V$\n",
    "* Solutions of Least-squares problems $A\\mathbf{x}=\\mathbf{b}$ are essentially projection problems where we seek to solve $A\\mathbf{x}=\\mathbf{b}$ where $\\mathbf{b}\\notin C(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution of Least Squares problems by the Normal Equations\n",
    "\n",
    "given $A\\mathbf{x}=\\mathbf{b}$ where $\\mathbf{b}\\notin C(A)$ we can always solve them using the Normal Equations\n",
    "\n",
    "$$\n",
    "A^T A\\hat{\\mathbf{x}} = A^T\\mathbf{b}\n",
    "$$ \n",
    "\n",
    "which actually solves $A\\hat{\\mathbf{x}} =\\mathbf{p}$  where $\\mathbf{p}$ is the orthogonal projection of $\\mathbf{b}$ onto $C(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "However...as you will show in the homework,  this can be numerically innaccurate because the condition number\n",
    "\n",
    "$$ \\kappa(A^T A) = \\kappa(A)^2 $$\n",
    "\n",
    "but there is a better way"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution of Least Squares problems by the QR factorization\n",
    "\n",
    "given any matrix $A$ that is full column rank, we will show that we can always factor it as \n",
    "\n",
    "$$\n",
    "    A=QR\n",
    "$$\n",
    "\n",
    "where $Q$ is a unitary matrix whose columns form an *orthonormal* basis for $C(A)$,  and $R$ is an upper triangular matrix that says how to reconstruct the columns of $A$ from the columns of $Q$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{a}_1 & \\cdots & \\mathbf{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix} = \n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{q}_1 & \\cdots & \\mathbf{q}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix}\n",
    "    \\begin{bmatrix} r_{11} & r_{12} & \\cdots & r_{1n} \\\\  & r_{22} &  &  \\\\  &  & \\ddots & \\vdots \\\\  &  &  & r_{nn} \\end{bmatrix}.\n",
    "$$\n",
    "If we write this out as a matrix multiplication we have\n",
    "$$\\begin{aligned}\n",
    "    \\mathbf{a}_1 &= r_{11} \\mathbf{q}_1 \\\\\n",
    "    \\mathbf{a}_2 &= r_{22} \\mathbf{q}_2 + r_{12} \\mathbf{q}_1 \\\\\n",
    "    \\mathbf{a}_3 &= r_{33} \\mathbf{q}_3 + r_{23} \\mathbf{q}_2 + r_{13} \\mathbf{q}_1 \\\\\n",
    "    &\\vdots\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "i.e we can construct the columns of $A$ from linear combinations of the columns of $Q$. (And those specific combinations come from $R=Q^TA$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "given $A=QR$,  then solving\n",
    "\n",
    "$$\n",
    "    A\\mathbf{x} = \\mathbf{b}\n",
    "$$\n",
    "\n",
    "becomes\n",
    "\n",
    "$$\n",
    "    QR\\mathbf{x} = \\mathbf{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or since $Q^T Q=I$\n",
    "\n",
    "$$\n",
    "    R\\mathbf{x} = Q^T\\mathbf{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which can be solve quickly by back-substitution as it is a triangular matrix.\n",
    "\n",
    "Moreover,  this problem is much better conditioned as it can be shown that $\\kappa(R)=\\kappa(A)$ not $\\kappa(A)^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Multiplying both sides by $Q$ again shows that this problem is equivalent to\n",
    "\n",
    "$$\n",
    "    QR\\mathbf{x} = QQ^T\\mathbf{b}\n",
    "$$\n",
    "\n",
    "or $A\\mathbf{x} = QQ^T\\mathbf{b}$ i.e. we are just solving $A\\mathbf{x}=\\mathbf{p}$ where $\\mathbf{p}$ is the orthogonal projection of  $\\mathbf{b}$ onto $C(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Question**... how to find $QR$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR Factorization\n",
    "\n",
    "One of the most important ideas in linear algebra is the concept of factorizing an original matrix into different constituents that may have useful properties.  These properties can help us understand the matrix better and lead to numerical methods.  In numerical linear algebra one of the most important factorizations is the **QR factorization**.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "There are actually multiple algorithm's that accomplish the same factorization but with very different methods and numerical stability.  Here we will discuss three of the algorithms\n",
    "\n",
    "1. **Classical Gram-Schmidt Orthogonalization**: Transform $A\\rightarrow Q$ by a succession of projections and calculate $R$ as a by-product of the algorithm.  (Unfortunately, prone to numerical floating point error)\n",
    "2. **Modified Gram-Schmidt Orthogonalization**: Transform $A\\rightarrow Q$ by a different set of projections.  Yields the same $QR$ but more numerically stable.\n",
    "3. **Householder Triangularization**:  Transform $A\\rightarrow R$ by a series of Unitary transformations, can solve least squares problems directly without accumulating $Q$, or can build $Q$ on the fly.\n",
    "\n",
    "And want to find a sequence of orthonormal vectors $\\mathbf{q}_j$ that span the sequence of subspaces\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Classical Gram-Schmidt\n",
    "\n",
    "We begin with a matrix $A\\in\\mathbb C^{m\\times n} $ with $n$ linearly independent columns. \n",
    "$$\n",
    "    A = \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{a}_1 & \\cdots & \\mathbf{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix}\n",
    "$$\n",
    "And want to find a sequence of orthonormal vectors $\\mathbf{q}_j$ that span the sequence of subspaces\n",
    "\n",
    "$$\n",
    "    \\text{span}(\\mathbf{a}_1) \\subseteq \\text{span}(\\mathbf{a}_1, \\mathbf{a}_2) \\subseteq \\text{span}(\\mathbf{a}_1, \\mathbf{a}_2, \\mathbf{a}_3) \\subseteq \\cdots \\subseteq \\text{span}(\\mathbf{a}_1, \\mathbf{a}_2, \\ldots , \\mathbf{a}_n)\n",
    "$$\n",
    "\n",
    "where here $\\text{span}(\\mathbf{v}_1,\\mathbf{v}_2,\\ldots,\\mathbf{v}_m)$ indicates the subspace spanned by the vectors $\\mathbf{v}_1, \\mathbf{v}_2, \\ldots, \\mathbf{v}_m$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The individual $\\mathbf{q}_j$ will form the columns of the matrix $Q$ such that $C(A_k)=C(Q_k)$ where the $A_k$ is a matrix with the first $k$ columns of A (or Q)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Starting with the first vector $\\mathbf{a}_1$, this forms the basis for  a 1-dimensional subspace of $\\mathbb{R}^m$ (i.e. a line).  Thus $\\mathbf{q}_1$ is a unit vector in that line or"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "    \\mathbf{q}_1 = \\frac{\\mathbf{a}_1}{||\\mathbf{a}_1||}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For $\\text{span}(\\mathbf{a}_1, \\mathbf{a}_2)$ (which is a plane in $\\mathbb{R}^m$ we already have $\\mathbf{q}_1$ so we need to find a vector $\\mathbf{q}_2$ that is orthogonal to $\\mathbf{q}_1$, but still in the plane.\n",
    "\n",
    "An obvious option is to find the component of $\\mathbf{a}_2$ that is orthogonal to $\\mathbf{q}_1$ but this is just"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "    \\mathbf{v}_2 &= (I - \\mathbf{q}_1\\mathbf{q}_1^T)\\mathbf{a}_2 \\\\\n",
    "    &= \\mathbf{a}_2 - \\mathbf{q}_1(\\mathbf{q}_1^T\\mathbf{a}_2)\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "By construction it's easy to show that $\\mathbf{v_2}$ is orthogonal to $\\mathbf{q}_1$ but not necessarily a unit vector but we can find $\\mathbf{q}_2$ by normalizing\n",
    "\n",
    "$$\n",
    "    \\mathbf{q}_2 = \\frac{\\mathbf{v}_2}{||\\mathbf{v}_2||}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Classical Gram-Schmidt as a series of projections\n",
    "If we define the unitary matrix \n",
    "$$\n",
    "    Q_k =     \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{q}_1 & \\cdots & \\mathbf{q}_k \\\\  &  &  \\\\  &  &  \\end{bmatrix}\n",
    "$$ \n",
    "as the first $k$ columns of Q,  then we could also rewrite the first two steps  of Classical Gram-Schmidt as\n",
    "\n",
    "\n",
    "set $\\mathbf{v}_1 = \\mathbf{a}_1$, then\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbf{q}_1 &= \\frac{\\mathbf{v}_1}{||\\mathbf{v}_1||}\\\\\n",
    "\\mathbf{v_2} &= (I - Q_1Q_1^T)\\mathbf{a}_2 = \\mathbf{a}_2 - Q_1(Q_1^T\\mathbf{a}_2)\\\\\n",
    "\\mathbf{q}_2 &= \\frac{\\mathbf{v}_2}{||\\mathbf{v}_2||}\\\\\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and we can continue\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbf{v_3} &= (I - Q_2Q_2^T)\\mathbf{a}_3 = \\mathbf{a}_3 - Q_2(Q_2^T\\mathbf{a}_3)\\\\\n",
    "\\mathbf{q}_3 &= \\frac{\\mathbf{v}_3}{||\\mathbf{v}_3||}\\\\\n",
    "& \\vdots \\\\\n",
    "\\mathbf{v_k} &= (I - Q_{k-1}Q_{k-1}^T)\\mathbf{a}_k = \\mathbf{a}_k - Q_{k-1}(Q_{k-1}^T\\mathbf{a}_k)\\\\\n",
    "\\mathbf{q}_k &= \\frac{\\mathbf{v}_k}{||\\mathbf{v}_k||}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "With each step finding the component of $a_k$ that is orthogonal to all the other vectors before it, and normalizing it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A picture is probably useful here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### But what about $R$?\n",
    "\n",
    "The above algorithm appears to transform $A$ directly to $Q$ without calculating $R$ (which we need for least-squares problems).\n",
    "\n",
    "But actually that's not true as $R$ is hiding in the algorithm (much like $L$ hides in the $LU$ factorizations $A\\rightarrow U$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If we consider the full factorization $A = QR$\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{a}_1 & \\cdots & \\mathbf{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix} = \n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{q}_1 & \\cdots & \\mathbf{q}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix}\n",
    "    \\begin{bmatrix} r_{11} & r_{12} & \\cdots & r_{1n} \\\\  & r_{22} &  &  \\\\  &  & \\ddots & \\vdots \\\\  &  &  & r_{nn} \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "it follows that $R=Q^T A$ or\n",
    "$$ \\begin{bmatrix} r_{11} & r_{12} & \\cdots & r_{1n} \\\\  & r_{22} &  &  \\\\  &  & \\ddots & \\vdots \\\\  &  &  & r_{nn}\\end{bmatrix} =\n",
    "   \\begin{bmatrix}  & -\\mathbf{q}_1^T- &  \\\\  & -\\mathbf{q}_2^T- &  \\\\  & \\vdots &  \\\\  & -\\mathbf{q}_n^T-  &  \\\\ \\end{bmatrix}\n",
    "\\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\mathbf{a}_1 & \\cdots & \\mathbf{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix}  \n",
    "    .\n",
    "$$\n",
    "Which with a little bit of work, you can show that the $j$th column of $R$\n",
    "$$\n",
    "    \\mathbf{r}_j = Q_j^T\\mathbf{a}_j\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So the full Classical Gram-Schmidt algorithm looks something like\n",
    "\n",
    "Set $\\mathbf{v}_1 = \\mathbf{a}_1$, $R_{11}=||\\mathbf{v}_1||$, $\\mathbf{q}_1 = \\mathbf{v}_1/R_{11}$\n",
    "\n",
    "Loop over columns for $j=2,\\ldots,n$\n",
    "\n",
    "* Find $\\mathbf{r}_j = Q^T_{j-1}\\mathbf{a_j}$\n",
    "* $\\mathbf{v}_j = \\mathbf{a}_j - Q_{j-1}\\mathbf{r}_j$\n",
    "* $R_{jj} = ||\\mathbf{v}_j||_2$\n",
    "* $\\mathbf{q}_j = \\mathbf{v}_j/R_{jj}$\n",
    "\n",
    "Which builds up both $Q$ and $R$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Or unrolled\n",
    "$$\\begin{aligned}\n",
    "    \\mathbf{q}_1 &= \\frac{\\mathbf{a}_1}{r_{11}} \\\\\n",
    "    \\mathbf{q}_2 &= \\frac{\\mathbf{a}_2 - r_{12} \\mathbf{q}_1}{r_{22}} \\\\\n",
    "    \\mathbf{q}_3 &= \\frac{\\mathbf{a}_3 - r_{13} \\mathbf{q}_1 - r_{23} \\mathbf{q}_2}{r_{33}} \\\\\n",
    "    &\\vdots \\\\\n",
    "    \\mathbf{q}_n &= \\frac{\\mathbf{a}_n - \\sum^{n-1}_{i=1} r_{in} \\mathbf{q}_i}{r_{nn}}\n",
    "\\end{aligned}$$\n",
    "leading us to define\n",
    "$$\n",
    "    r_{ij} = \\left \\{ \\begin{aligned}\n",
    "        &\\langle \\mathbf{q}_i, \\mathbf{a}_j \\rangle & &i \\neq j \\\\\n",
    "        &\\left \\Vert \\mathbf{a}_j - \\sum^{j-1}_{i=1} r_{ij} \\mathbf{q}_i \\right \\Vert & &i = j\n",
    "    \\end{aligned} \\right .\n",
    "$$\n",
    "\n",
    "This is called the **classical Gram-Schmidt** iteration.  Turns out that the procedure above is unstable because of rounding errors introduced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Which can be easily coded up in python as "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Classical Gram-Schmidt Iteration\n",
    "def classic_GS(A):\n",
    "    m, n = A.shape\n",
    "    Q = numpy.empty((m, n))\n",
    "    R = numpy.zeros((n, n))\n",
    "    # loop over columns\n",
    "    for j in range(n):\n",
    "        v = A[:, j]\n",
    "        for i in range(j):\n",
    "            R[i, j] = numpy.dot(Q[:, i].conjugate(), A[:, j])\n",
    "            v = v - R[i, j] * Q[:, i]\n",
    "        R[j, j] = numpy.linalg.norm(v, ord=2)\n",
    "        Q[:, j] = v / R[j, j]\n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### And check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Q, R = classic_GS(A)\n",
    "print('Q=\\n{}\\n'.format(Q))\n",
    "print('Q^TQ=\\n{}'.format(numpy.dot(Q.transpose(), Q)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('R=\\n{}\\n'.format(R))\n",
    "print('QR - A=\\n{}'.format(numpy.dot(Q, R) - A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### Full vs. Reduced QR\n",
    "\n",
    "If the original matrix $A \\in \\mathbb C^{m \\times n}$ where $m \\ge n$ then we can still define a QR factorization, called the **full QR factorization**, which appends columns full of zeros to $R$ to reproduce the full matrix.\n",
    "$$\n",
    "    A = Q R = \\begin{bmatrix} Q_1 & Q_2 \\end{bmatrix} \n",
    "              \\begin{bmatrix} R_1 \\\\\n",
    "                              0 \n",
    "              \\end{bmatrix} = Q_1 R_1\n",
    "$$\n",
    "The factorization $Q_1 R_1$ is called the **reduced** or **thin QR factorization** of $A$.\n",
    "\n",
    "We require that the additional columns added $Q_2$ are an orthonormal basis that is orthogonal itself to $\\text{range}(A)$.  If $A$ is full ranked then $Q_1$ and $Q_2$ provide a basis for $\\text{range}(A)$ and $\\text{null}(A^\\ast)$ respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### QR Existence and Uniqueness\n",
    "Two important theorems exist regarding this algorithm which we state without proof:\n",
    "\n",
    "*Every $A \\in \\mathbb C^{m \\times n}$ with $m \\geq n$ has a full QR factorization and therefore a reduced QR factorization.*\n",
    "\n",
    "*Each $A \\in \\mathbb C^{m \\times n}$ with $m \\geq n$ of full rank has a unique reduced QR factorization $A = QR$ with $r_{jj} > 0$.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Modified Gram-Schmidt\n",
    "\n",
    "Unfortunately the classical Gram-Schmidt algorithm is is not stable numerically.  Instead we can derive a modified method that is more numerically stable but calculates the same $Q$ and $R$ with just a different order of projections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Recall that the basic piece of the original algorithm was to take the inner product of $\\mathbf{a}_j$ and all the relevant $\\mathbf{q}_i$.  Using the rewritten version of Gram-Schmidt in terms of projections we then have\n",
    "\n",
    "$$\n",
    "    \\mathbf{v}_j = (I - Q_{j-1}Q^T_{j-1}) \\mathbf{a}_j.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Each of these projections is of different rank $m - (j - 1)$ although we know that the resulting $\\mathbf{v}_j$ are linearly independent by construction.  The modified version of Gram-Schmidt instead uses projections that are all of rank $m-1$.  To construct this projection remember that we can again construct the complement to a projection and perform the following sequence of projections\n",
    "\n",
    "$$\n",
    "    P_j = \\hat{\\!P}_{\\mathbf{q}_{j-1}} \\hat{\\!P}_{\\mathbf{q}_{j-2}} \\cdots \\hat{\\!P}_{\\mathbf{q}_{2}} \\hat{\\!P}_{\\mathbf{q}_{1}}\n",
    "$$\n",
    "\n",
    "where \n",
    "$$\\hat{\\!P}_{\\mathbf{q}_{j}} = I - \\mathbf{q}_i\\mathbf{q_i}^T$$\n",
    "\n",
    "which projects onto the complementary space orthogonal to $\\mathbf{q}_i$. \n",
    "\n",
    "Note that this performs mathematically the same job as $P_i \\mathbf{a}_i$ however each of these projectors are of rank $m - 1$.  The reason why this approach is more stable is that we are not projecting with a possibly arbitrarily low-rank projector, instead we only take projectors that are high-rank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "again...a picture is probably worth a lot here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This leads to the following set of calculations:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    1.\\quad \\mathbf{v}^{(1)}_i &= \\mathbf{a}_i  \\\\\n",
    "    2.\\quad \\mathbf{v}^{(2)}_i &= \\hat{\\!P}_{\\mathbf{q}_1} \\mathbf{v}_i^{(1)} = \\mathbf{v}^{(1)}_i - \\mathbf{q}_1 q_1^\\ast \\mathbf{v}^{(1)}_i \\\\\n",
    "    3.\\quad \\mathbf{v}^{(3)}_i &= \\hat{\\!P}_{\\mathbf{q}_2} \\mathbf{v}_i^{(2)} = \\mathbf{v}^{(2)}_i - \\mathbf{q}_2 \\mathbf{q}_2^\\ast \\mathbf{v}^{(2)}_i \\\\\n",
    "    & \\text{  } \\vdots & &\\\\\n",
    "    i.\\quad \\mathbf{v}^{(i)}_i &= \\hat{\\!P}_{\\mathbf{q}_{i-1}} \\mathbf{v}_i^{(i-1)} =  \\mathbf{v}_i^{(i-1)} - \\mathbf{q}_{i-1} \\mathbf{q}_{i-1}^\\ast \\mathbf{v}^{(i-1)}_i\n",
    "\\end{aligned}$$\n",
    "\n",
    "The reason why this approach is more stable is that we are not projecting with a possibly arbitrarily low-rank projector, instead we only take projectors that are high-rank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Example: Implementation of modified Gram-Schmidt**\n",
    "Implement the modified Gram-Schmidt algorithm checking to make sure the resulting factorization has the required properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Modified Gram-Schmidt Iteration\n",
    "def mod_GS(A):\n",
    "    m, n  = A.shape\n",
    "    Q = numpy.empty((m, n))\n",
    "    R = numpy.zeros((n, n))\n",
    "    v = A.copy()\n",
    "    for i in range(n):\n",
    "        R[i, i] = numpy.linalg.norm(v[:, i], ord=2)\n",
    "        Q[:, i] = v[:, i] / R[i, i]\n",
    "        for j in range(i + 1, n):\n",
    "            R[i, j] = numpy.dot(Q[:, i].conjugate(), v[:, j])\n",
    "            v[:, j] -= R[i, j] * Q[:, i]\n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### And check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Q, R = mod_GS(A)\n",
    "print('Q=\\n{}\\n'.format(Q))\n",
    "print('Q^TQ=\\n{}'.format(numpy.dot(Q.transpose(), Q)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print('R=\\n{}\\n'.format(R))\n",
    "print('QR - A=\\n{}'.format(numpy.dot(Q, R) - A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Householder Triangularization\n",
    "\n",
    "One way to also interpret Gram-Schmidt orthogonalization is as a series of right multiplications by upper triangular matrices of the matrix A.  For instance the first step in performing the modified algorithm is to divide through by the norm $r_{11} = ||v_1||$ to give $q_1$:\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "           &    &       &  \\\\\n",
    "           &    &       &  \\\\\n",
    "        \\mathbf{v}_1 & \\mathbf{v}_2 & \\mathbf{v}_3 &  \\cdots & \\mathbf{v}_n \\\\\n",
    "           &    &       &  \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "        \\frac{1}{r_{11}}   &  &\\cdots &       \\\\\n",
    "           & 1  &          \\\\\n",
    "         &  & \\ddots &  \\\\ \n",
    "           &    &  & 1 \n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "           &    &       &        \\\\\n",
    "           &    &       &        \\\\\n",
    "        \\mathbf{q}_1 & \\mathbf{v}_2 &  \\mathbf{v}_3 & \\cdots & \\mathbf{v}_n \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can also perform all the step (2) evaluations by also combining the step that projects onto the complement of $q_1$ by add the appropriate values to the entire first row:\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "           &    &       &  \\\\\n",
    "           &    &       &  \\\\\n",
    "        \\mathbf{v}_1 & \\mathbf{v}_2 &  \\mathbf{v}_3 & \\cdots & \\mathbf{v}_n \\\\\n",
    "           &    &       &  \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "        \\frac{1}{r_{11}}   & -\\frac{r_{12}}{r_{11}} & -\\frac{r_{13}}{r_{11}} & \\cdots      \\\\\n",
    "           & 1   &         \\\\\n",
    "         &  & \\ddots &  \\\\ \n",
    "           &    &  & 1 \n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "           &    &       &        \\\\\n",
    "           &    &       &        \\\\\n",
    "        \\mathbf{q}_1 & \\mathbf{v}_2^{(2)} & \\mathbf{v}_3^{(2)} & \\cdots & \\mathbf{v}_n^{(2)} \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The next step can then be placed into the second row:\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "           &    &       &  \\\\\n",
    "           &    &       &  \\\\\n",
    "        \\mathbf{v}_1 & \\mathbf{v}_2 & \\mathbf{v}_3 & \\cdots & \\mathbf{v}_n \\\\\n",
    "           &    &       &  \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "    \\cdot R_1 \\cdot\n",
    "    \\begin{bmatrix}\n",
    "        1 &  &  &  & \\\\\n",
    "         & \\frac{1}{r_{22}}   & -\\frac{r_{23}}{r_{22}} & -\\frac{r_{25}}{r_{22}} & \\cdots      \\\\\n",
    "         &   & 1   &         \\\\\n",
    "         &  &  & \\ddots &  \\\\ \n",
    "         &    &  & & 1 \n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "           &    &       &        \\\\\n",
    "           &    &       &        \\\\\n",
    "        \\mathbf{q}_1 & \\mathbf{q}_2 & \\mathbf{v}_3^{(3)} & \\cdots & \\mathbf{v}_n^{(3)} \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we identify the matrices as $R_1$ for the first case, $R_2$ for the second case and so on we can write the algorithm as\n",
    "\n",
    "$$\n",
    "    A \\underbrace{R_1R_2 \\quad \\cdots \\quad R_n}_{\\hat{R}^{-1}} = \\hat{\\!Q}.\n",
    "$$\n",
    "\n",
    "This view of Gram-Schmidt is called Gram-Schmidt triangularization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Householder triangularization is similar in spirit.  Instead of multiplying $A$ on the right Householder multiplies $A$ on the left by unitary matrices $Q_k$.  Remember that a unitary matrix (or an orthogonal matrix when strictly real) has as its inverse its adjoint (transpose when real) $Q^* = Q^{-1}$ so that $Q^* Q = I$.  We therefore have\n",
    "\n",
    "$$\n",
    "    Q_n Q_{n-1} \\quad \\cdots \\quad Q_2 Q_1 A = R\n",
    "$$\n",
    "\n",
    "which if we identify $Q_n Q_{n-1} \\text{  } \\cdots \\text{  } Q_2 Q_1 = Q^*$ and note that if $Q = Q^\\ast_n Q^\\ast_{n-1} \\text{  } \\cdots \\text{  } Q^\\ast_2 Q^\\ast_1$ then $Q$ is also unitary.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Householder Triangularization\n",
    "\n",
    "While Gram-Schmidt and its variants transform $A\\rightarrow Q$ and calculate $R$ on the way.  Householder triangularization is a rather different algorithm that transforms $A\\rightarrow R$ by multiplying $A$ by a series of *square* Unitary matrices that systematically put zeros in the subdiagonal (similar to the LU factorization) \n",
    "\n",
    "$$\n",
    "    Q_n Q_{n-1} \\quad \\cdots \\quad Q_2 Q_1 A = R\n",
    "$$\n",
    "\n",
    "which if we identify $Q^* = Q_n Q_{n-1} \\text{  } \\cdots \\text{  } Q_2 Q_1$ then we can define  $Q = Q^\\ast_1 Q^\\ast_{2} \\text{  } \\cdots \\text{  } Q^\\ast_{n-1} Q^\\ast_n$ which is also unitary and \n",
    "\n",
    "$$\n",
    "    A = QR\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can then write this as\n",
    "$$\\begin{aligned}\n",
    "    Q_n Q_{n-1} \\text{  } \\cdots \\text{  } Q_2 Q_1 A &= R \\\\\n",
    "    Q_{n-1} \\cdots Q_2 Q_1 A &= Q^\\ast_n R \\\\\n",
    "    & \\text{  } \\vdots \\\\\n",
    "    A &= Q^\\ast_1 Q^\\ast_2 \\text{  } \\cdots \\text{  } Q^\\ast_{n-1} Q^\\ast_n R \\\\\n",
    "    A &= Q R\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This was we can think of Householder triangularization as one of introducing zeros into $A$ via orthogonal matrices.\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "    \\end{bmatrix} \\overset{Q_1}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{+} & \\text{+} & \\text{+} \\\\\n",
    "         & \\text{+} & \\text{+} \\\\\n",
    "         & \\text{+} & \\text{+} \\\\\n",
    "         & \\text{+} & \\text{+} \\\\\n",
    "    \\end{bmatrix} \\overset{Q_2}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{+} & \\text{+} & \\text{+} \\\\\n",
    "         & \\text{-} & \\text{-} \\\\\n",
    "         &  & \\text{-} \\\\\n",
    "         &  & \\text{-} \\\\\n",
    "    \\end{bmatrix} \\overset{Q_3}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{+} & \\text{+} & \\text{+} \\\\\n",
    "         & \\text{-} & \\text{-} \\\\\n",
    "         &  & \\text{*} \\\\\n",
    "         &  & \\\\\n",
    "    \\end{bmatrix} = R\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now the question is how do we construct the $Q_k$.  The construction is usually broken down into a block matrix of the form\n",
    "\n",
    "$$\n",
    "    Q_k = \\begin{bmatrix} I & 0 \\\\ 0 & H \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $I \\in \\mathbb C^{k-1 \\times k-1}$ is the identity matrix that preserves the top $k-1$ rows\n",
    "\n",
    "and $H \\in \\mathbb C^{m - (k - 1) \\times m - (k-1)}$ is a unitary matrix that just modifies the lower $m-k-1$ rows.\n",
    "\n",
    "Note that this will leave the rows and columns we have already worked on alone and be unitary.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We now turn to the task of constructing the matrix $H$. Note that the definition of $Q_k$ implies that\n",
    "$$\n",
    "\\begin{align}\n",
    "Q_k \\mathbf{x} &= \\begin{bmatrix} I & 0 \\\\ 0 & H \\end{bmatrix} \\mathbf{x} \n",
    "     &= \\begin{bmatrix} I & 0 \\\\ 0 & H \\end{bmatrix} \n",
    "     \\begin{bmatrix} \\mathbf{x}_{top} \\\\ \\mathbf{x}_{bottom} \\end{bmatrix} \\\\\n",
    "     &= \\begin{bmatrix} \\mathbf{x}_{top} \\\\ H\\mathbf{x}_{bottom}  \\end{bmatrix}. \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "The vectors $\\mathbf{x}_{top}$ and $\\mathbf{x}_{bottom}$ are defined to be\n",
    "$$\n",
    "\\begin{align}\n",
    "     \\mathbf{x}_{top} &= \\begin{bmatrix} x_1  \\\\  x_2 \\\\ \\vdots \\\\ x_{k-1}\\end{bmatrix} \\\\\n",
    "     \\mathbf{x}_{bottom} &= \\begin{bmatrix} x_k \\\\ x_{k+1} \\\\ \\vdots \\\\ x_m \\end{bmatrix} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "The task is to construct the matrix $H$ so that after multiplying, $H\\mathbf{x}_{bottom}$  results in a vector that has some nonzero number in the top position and the rest of the numbers are all zeroes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or  \n",
    "\n",
    "$$H\\mathbf{x}_{bottom} = \\begin{bmatrix} \\alpha \\\\ 0 \\\\ \\vdots \\\\ 0 \\end{bmatrix}$$\n",
    "\n",
    "But if $H$ is unitary, it can't change the length of a vector which implies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$ \\alpha = ||\\mathbf{x}_{bottom}||$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One unitary operation that does this is the *Householder Reflection Matrix*  that reflects the vector $\\mathbf{x}$ over the hyper plane $H$ that is normal to the vector $\\mathbf{q}$ so that $H \\mathbf{x} = ||x|| \\hat{\\mathbf{e}}_1$:\n",
    "<table>\n",
    "    <tr align='center'>\n",
    "<img align=center, src=\"./images/householder.png\" alt=\"Drawing\" width=1000/>\n",
    "    <tr\\>\n",
    "<table\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or mathematically\n",
    "$$\n",
    "    \\mathbf{x} = \\begin{bmatrix}\n",
    "        x_1 \\\\\n",
    "        x_2 \\\\\n",
    "        \\vdots \\\\\n",
    "        x_{n} \n",
    "     \\end{bmatrix}, \\quad\n",
    "    H\\mathbf{x} = \\begin{bmatrix}\n",
    "        ||\\mathbf{x}|| \\\\\n",
    "        0 \\\\\n",
    "        \\vdots \\\\\n",
    "        0\n",
    "    \\end{bmatrix} = ||\\mathbf{x}|| \\mathbf{e}_1.\\quad\\text{where}\\quad\\mathbf{e}_1 = \\begin{bmatrix}\n",
    "        1 \\\\\n",
    "        0 \\\\\n",
    "        \\vdots \\\\\n",
    "        0\n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "is the first unit vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This is of course the effect on only one vector.  Any other vector will be reflected across $H$ (technically a hyperplane) which is orthogonal to \n",
    "\n",
    "$$\\mathbf{v} = ||\\mathbf{x}|| \\hat{e}_1 - \\mathbf{x}.$$\n",
    "\n",
    "This has a similar construction as to the projector complements we were working with before.  Consider the projector defined as\n",
    "\n",
    "$$\n",
    "    P x = \\left (I - \\mathbf{q} \\mathbf{q}^T\\right)\\mathbf{x} = \\mathbf{x} - \\mathbf{q}(\\mathbf{q}^T\\mathbf{x}) \n",
    "$$\n",
    "\n",
    "where $$\\mathbf{q} = \\frac{\\mathbf{v}}{||\\mathbf{v}||}.$$\n",
    "\n",
    "This vector $P\\mathbf{x}$ is now orthogonal to $\\mathbf{v}$. i.e. lies in the plane $H$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we actually want to transform $\\mathbf{x}$ to lie in the direction of $\\hat{e}_1$ we need to go twice as far as which allows us to identify the matrix $H$ as\n",
    "\n",
    "$$\n",
    "    H = I - 2 \\mathbf{q} \\mathbf{q}^T.\n",
    "$$\n",
    "\n",
    "<table>\n",
    "    <tr align='center'>\n",
    "<img align=center, src=\"./images/householder.png\" alt=\"Drawing\" width=1000/>\n",
    "    <tr\\>\n",
    "<table\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There is actually a non-uniqueness to which direction we reflect over since another definition of $\\hat{H}$ which is orthogonal to the one we originally choose is available.  For numerical stability purposes we will choose the reflector that is the most different from $\\mathbf{x}$.  This comes back to having difficulties numerically when the vector $\\mathbf{x}$ is nearly aligned with $\\hat{e}_1$ and therefore one of the $H$ specification.  By convention the $\\mathbf{v}$ chosen is defined by\n",
    "\n",
    "$$\n",
    "    \\mathbf{v} = \\text{sign}(x_1)||\\mathbf{x}|| \\hat{e}_1 + \\mathbf{x}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Implementation of Householder QR Factorization\n",
    "def householder_QR(A, verbose=False):\n",
    "    R = A.copy()\n",
    "    m, n = A.shape\n",
    "    QT = numpy.eye(m)\n",
    "    for k in range(n):\n",
    "        x = numpy.zeros(m)\n",
    "        e = numpy.zeros(m)\n",
    "        x[k:] = R[k:, k]\n",
    "        e[k] = 1.0\n",
    "        # simplest version v = ||x||e - x\n",
    "        #v =  numpy.linalg.norm(x, ord=2) * e - x\n",
    "        # alternate version\n",
    "        v = numpy.sign(x[k]) * numpy.linalg.norm(x, ord=2) * e + x\n",
    "        v = v / numpy.linalg.norm(v, ord=2)\n",
    "        R  -= 2.0 * numpy.outer(v,numpy.dot(v.T,R))\n",
    "        QT  -= 2.0 * numpy.outer(v,numpy.dot(v.T,QT))\n",
    "    Q = QT.T[:,:n]   \n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "print(\"Matrix A = \")\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%precision 6\n",
    "Q, R = householder_QR(A, verbose=False)\n",
    "print(\"Householder (reduced) Q =\\n{}\\n\".format(Q))\n",
    "print(\"Householder (full) R = \")\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "m, n = A.shape\n",
    "print(\"Check to see if factorization worked...||A -QR|| = {}\".format(numpy.linalg.norm(A - numpy.dot(Q, R[:n, :n]))))\n",
    "print(A - numpy.dot(Q, R[:n, :n]))\n",
    "print(\"\\nCheck if Q is unitary...||Q^TQ -I|| = {}\".format(numpy.linalg.norm(Q.T.dot(Q)-numpy.eye(n))))\n",
    "print(Q.T.dot(Q))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Comparison of accuracy of the different algorithms\n",
    "\n",
    "As it turns out, not all $QR$ algorithms produce the same quality of orthogonalization.  Here we provide a few examples that compare the behavior of the 3 different algorithms and `numpy.linalg.qr`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### Example 1:   Random Matrix QR\n",
    "\n",
    "Here we construct a large matrix $A$ with a random eigenspace and widely varying eigenvalues.  The values along the diagonal of $R$ gives us some idea of the size of the projections as we go, i.e. the larger the values the less effective we are in constructing orthogonal directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "N = 80\n",
    "# construct a random matrix with known singular values\n",
    "U, X = numpy.linalg.qr(numpy.random.random((N, N)))\n",
    "V, X = numpy.linalg.qr(numpy.random.random((N, N)))\n",
    "S = numpy.diag(2.0**numpy.arange(-1.0, -(N + 1), -1.0))\n",
    "A = numpy.dot(U, numpy.dot(S, V))\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "Q, R = classic_GS(A)\n",
    "axes.semilogy(numpy.diag(R), 'bo', label=\"Classic\")\n",
    "Q, R = mod_GS(A)\n",
    "axes.semilogy(numpy.diag(R), 'ro', label=\"Modified\")\n",
    "Q, R = householder_QR(A)\n",
    "axes.semilogy(numpy.diag(R), 'ko', label=\"Householder\")\n",
    "Q, R = numpy.linalg.qr(A)\n",
    "axes.semilogy(numpy.diag(R), 'go', label=\"numpy\")\n",
    "\n",
    "axes.set_xlabel(\"Index\", fontsize=16)\n",
    "axes.set_ylabel(\"$R_{ii}$\", fontsize=16)\n",
    "axes.legend(loc=3, fontsize=14)\n",
    "axes.plot(numpy.arange(0, N), numpy.ones(N) * numpy.sqrt(numpy.finfo(float).eps), 'k--')\n",
    "axes.plot(numpy.arange(0, N), numpy.ones(N) * numpy.finfo(float).eps, 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 2:  Comparing Orthogonality\n",
    "\n",
    "Consider the $QR$ factorization of the ill-conditioned matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "        0.70000 & 0.70711 \\\\ 0.70001 & 0.70711\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%precision 16\n",
    "\n",
    "A = numpy.array([[0.7, 0.70711], [0.70001, 0.70711]])\n",
    "print('\\ncond(A) = {}'.format(numpy.linalg.cond(A)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To check that the matrix $Q$ is really unitary, we compute $A=QR$ with the different algorithms and compare\n",
    "\n",
    "$$\n",
    "    ||Q^TQ - I||\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "Q_c, R = classic_GS(A)\n",
    "r_c = numpy.linalg.norm(numpy.dot(Q_c.transpose(), Q_c) - numpy.eye(2))\n",
    "print(\"Classic: \", r_c )\n",
    "\n",
    "Q, R = mod_GS(A)\n",
    "print(\"Modified: \", numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2)))\n",
    "\n",
    "Q_h, R = householder_QR(A)\n",
    "r_h = numpy.linalg.norm(numpy.dot(Q_h.transpose(), Q_h) - numpy.eye(2))\n",
    "print(\"Householder:\", r_h)\n",
    "\n",
    "Q, R = numpy.linalg.qr(A)\n",
    "r = numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2))\n",
    "print(\"Numpy: \", r)\n",
    "\n",
    "print('\\ncond(A) = {}'.format(numpy.linalg.cond(A)))\n",
    "print('r_classic/r_householder = {}'.format(r_c/r_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Applications of the QR\n",
    "\n",
    "The $QR$ factorization and unitary transformation such as Householder reflections,  play important roles in a wide range of algorithms for Numerical Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Least-Squares problems: Solving $A\\mathbf{x} = \\mathbf{b}$ with QR\n",
    "\n",
    "We have already discussed solving overdetermined least-squares problems using the $QR$ factorization.  In general,  if we seek least-squares solutions to  $A\\mathbf{x}=\\mathbf{b}$ and $A=QR$, the problem reduces to solving the triangular problem \n",
    "$$\n",
    "R\\mathbf{x} = Q^T\\mathbf{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we use Householder triangularization to transform $A\\rightarrow R$ directly,  we do not need to explicitly form the matrix $Q$ and can  save memory and computation.  If we  consider the augmented system $\\begin{bmatrix} A & \\mathbf{b}\\end{bmatrix}\n",
    "$ and Apply the same sequence of unitary, rank 1 tranformations we used in the Householder algorithm, the sequence becomes\n",
    "\n",
    "$$\n",
    "Q_n\\ldots Q_2Q_1\\begin{bmatrix} A & \\mathbf{b}\\end{bmatrix} = \\begin{bmatrix} Q^TA & Q^T\\mathbf{b}\\end{bmatrix} =\\begin{bmatrix} R & \\mathbf{c}\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Thus we just need to solve\n",
    "$$ \n",
    "    R\\mathbf{x} = \\mathbf{c}\n",
    "$$ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Alternatively,  during the $QR$ factorization by Householder, we can just store the vectors $\\mathbf{q}_1\\ldots\\mathbf{q}_n$ and reconstruct $Q^T\\mathbf{b}$ by rank-1 updates as neccessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Finding Eigenvalues\n",
    "\n",
    "As it turns out,  the $QR$ factorization and Householder transformation are also extremely useful for finding eigenvalues of a matrix $A$.  There is an entire notebook `13_LA_eigen.ipynb` that develops these ideas in detail. Here we will just discuss the parts that relate to the $QR$ and orthogonalization algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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": "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": "skip"
    }
   },
   "source": [
    "**Example**\n",
    "\n",
    "The general Eigen problem\n",
    "$$\n",
    "A\\mathbf{x}=\\lambda\\mathbf{x}\n",
    "$$ \n",
    "\n",
    "is really $n$ problems for each eigenvalue, eigenvector pair i.e.\n",
    "\n",
    "$$\n",
    "    A\\mathbf{x}_i = \\lambda_i\\mathbf{x}_i\\quad\\text{for}\\quad i=1,\\ldots,n \n",
    "$$\n",
    "\n",
    "which can be written concisely in matrix form as\n",
    "\n",
    "$$\n",
    "AX =X\\Lambda\n",
    "$$\n",
    "\n",
    "where $X$ is a matrix whose columns contain the eigenvectors and $\\Lambda$ is a diagonal matrix of corresponding eigenvalues.  This form is always true"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example Diagonalizable matrices\n",
    "\n",
    "If a matrix $A\\in\\mathbb{R}^{n\\times n}$ has $n$ linearly independent eigenvectors, then  we say $A$ is *diagonalizable* and we can factor it as\n",
    "\n",
    "$$\n",
    "    A = X\\Lambda X^{-1}\n",
    "$$\n",
    "Where $X$ is a matrix of eigenvectors, and $\\Lambda$ is a diagonal matrix of corresponding Eigenvalues.\n",
    "\n",
    "which says $A$ is similar to $\\Lambda$ with similarity transform $X$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Theorem:\n",
    "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": "fragment"
    }
   },
   "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.\n",
    "\n",
    "**Partial Proof** for a diagonalizable matrix.  If $A$ is diagonalizable, $A=X\\Lambda X^{-1}$.  But we know we can always factor $X=QR$ and substitute to show\n",
    "\n",
    "$$\n",
    "    A = Q(R\\Lambda R^{-1})Q^T\n",
    "$$\n",
    "\n",
    "and it is not hard to show that the product $R\\Lambda R^{-1}$ is also an upper triangular matrix (exercise left to the reader).  \n",
    "\n",
    "(For a non-diagonalizable matrix the proof requires showing the existence of the Jordan form $A=MJM^{-1}$)"
   ]
  },
  {
   "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 Hermitian ($A^\\ast = A$)\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": "subslide"
    }
   },
   "source": [
    "### Hessenberg form\n",
    "\n",
    "The first step to finding the Schur factorization is to try and get $A$ as close to triangular as possible without changing its eigenvalues. This requires a series of similarity transformations. \n",
    "\n",
    "As it turns out,  the closest we can do is to reduce it to Hessenberg form which is upper triangular with one extra subdiagonal.  Which can be done with $n$ explicit similarity transformations using Householder Transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "slide"
    }
   },
   "source": [
    "## QR/RQ Algorithm\n",
    "\n",
    "Given a matrix in Hessenberg form, it turns out we can use repeated $QR$ factorizations to reduce the size of the subdiagonal and iterate towards the Schur factorization  to 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",
    "A0 = numpy.array([[2, 1, 1], [1, 3, 1], [1, 1, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "MAX_STEPS = 10\n",
    "\n",
    "A=A0\n",
    "print('A=')\n",
    "print(A)\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({}) =\".format(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.sort(numpy.linalg.eigvals(A0)))\n",
    "print()\n",
    "print(\"Computed eigenvalues: \")\n",
    "print(numpy.sort(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  \n",
    "\n",
    "$$A^{(k-1)} = Q^{(k)}R^{(k)}$$ \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which is equivalent to finding\n",
    "\n",
    "$$\n",
    "    (Q^{(k)})^T A^{(k-1)} = R^{(k)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and multiplying on the right leads to\n",
    "\n",
    "$$\n",
    "    (Q^{(k)})^T A^{(k-1)} Q^{(k)} = R^{(k)} Q^{(k)} = A^{(k)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "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 or 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",
    "$$\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": "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) = p_0 x^n + p_1 x^{n-1} + p_2 x^{n-2} + \\ldots + p_n\n",
    "$$\n",
    "\n",
    "described by a $n+1$ vector of coefficients $\\mathbf{p}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "p = numpy.array([ 1, 1, -1])\n",
    "r = numpy.roots(p)\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p = numpy.random.rand(8)\n",
    "r = numpy.roots(p)\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 \n",
    "\n",
    "$$c(x) = c_0 + c_1 x + c_2 x^2 + \\ldots + c_{n-1} x^{n-1} + x^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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def myroots(p, verbose=False):\n",
    "    ''' Calculate the roots of a polynomial described by coefficient vector \n",
    "    in numpy.roots order\n",
    "    p(x) = p_0 x^n + p_1 x^{n-1} + p_2 x^{n-2} + \\ldots + p_n   \n",
    "    by finding the eigenvalues of the companion matrix\n",
    "    \n",
    "    returns:\n",
    "    --------\n",
    "    eigenvalues sorted by |\\lambda|\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # construct the companion matrix of the coefficient vector c\n",
    "    # make p monic and reverse the order for this definition of the companion matrix\n",
    "    c = numpy.flip(p/p[0])\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",
    "    # calculate the eigenvalues of the companion matrix, then sort by |lambda|\n",
    "    eigs = numpy.linalg.eigvals(C)\n",
    "    index = numpy.flip(numpy.argsort(numpy.abs(eigs)))\n",
    "    return eigs[index]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "p = numpy.array([ 1, 1, -1])\n",
    "r = numpy.roots(p)\n",
    "print(r)\n",
    "mr = myroots(p) \n",
    "print\n",
    "print(mr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%precision 4\n",
    "p = numpy.random.rand(5)\n",
    "r = numpy.roots(p)\n",
    "print('nproots: {}'.format(r))\n",
    "mr = myroots(p) \n",
    "print()\n",
    "print('myroots: {}'.format(mr))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}