11_LA_QR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Note:  This material largely follows the text \"Numerical Linear Algebra\" by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "init_cell": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# QR Factorizations and Least Squares"
   ]
  },
  {
   "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 $v$ onto a lower dimensional subspace.  For example, suppose that $v$ lies completely within the subspace, i.e. $\\vec{v} \\in \\text{range}(P)$. If that is the case then $P \\vec{v}$ should not change, or $P\\vec{v} = \\vec{v}$.  This motivates the definition above.\n",
    "\n",
    "If $P\\vec{v} = \\vec{v}'$ then $\\vec{v}'$ should be in the sub-space and subsequent applications of $P$ should result in no change, i.e.\n",
    "$$\n",
    "    P( P \\vec{v} ) = P\\vec{v}' = \\vec{v}'.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As another example, take a vector $\\vec{x} \\notin \\text{range}(P)$ and project it onto the subspace $P\\vec{x} = \\vec{v}$.  If we apply the projection again to $\\vec{v}$ we now have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    P\\vec{x} &= \\vec{v}, \\\\\n",
    "    P^2 \\vec{x} & = P \\vec{v} = \\vec{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 $\\vec{x} \\notin \\text{range}(P)$, if we look at the difference between the projection and the original vector $P\\vec{x} - \\vec{x}$ and apply the projection again we have\n",
    "$$\n",
    "    P(P\\vec{x} - \\vec{x}) = P^2 \\vec{x} - P\\vec{x} = 0\n",
    "$$\n",
    "which means the difference between the projected vector $P\\vec{x} = \\vec{v}$ lies in the null space of $P$, $\\vec{v} \\in \\text{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",
    "    \\vec{x} \\in \\text{null}(P),\n",
    "$$ \n",
    "then \n",
    "$$\n",
    "    (I - P) \\vec{x} = \\vec{x} - P \\vec{x} = \\vec{x}\n",
    "$$ \n",
    "since $P \\vec{x} = 0$ implying that $\\vec{x} \\in \\text{range}(I - P)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We also know that \n",
    "$$\n",
    "    (I - P) \\vec{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",
    "    \\text{range}(I - P) \\cap \\text{range}(P) = \\{0\\}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    \\text{null}(P) \\cap \\text{range}(P) = \\{0\\}\n",
    "$$\n",
    "These two spaces are called **complementary subspaces**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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 $\\vec{s}\\in S = \\text{range}(P)$, and $\\vec{v} \\in V = \\text{null}(P)$.  If we have $\\vec{x} \\in \\mathbb C^{m \\times m}$ that we can split the vector $\\vec{x}$ into components in $S$ and $V$ by using the projections\n",
    "$$\\begin{aligned}\n",
    "    P \\vec{x} = \\vec{x}_S& &\\vec{x}_s \\in S \\\\\n",
    "    (I - P) \\vec{x} = \\vec{x}_V& &\\vec{x}_V \\in V\n",
    "\\end{aligned}$$\n",
    "which we can also observe adds to the original vector as\n",
    "$$\n",
    "    \\vec{x}_S + \\vec{x}_V = P \\vec{x} + (I - P) \\vec{x} = \\vec{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": "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!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A **hermitian** matrix is one whose complex conjugate is itself, i.e.\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": "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",
    "$$\n",
    "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:  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": [
    "$$\n",
    "    Q = \\begin{bmatrix} 1 & 0 \\\\ \n",
    "                        0 & 1 \\\\ \n",
    "                        0 & 0 \n",
    "        \\end{bmatrix}, \\quad \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": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Q = numpy.array([[1, 0],[0, 1],[0, 0]])\n",
    "P = numpy.dot(Q, numpy.conjugate(numpy.transpose(Q)))\n",
    "I = numpy.identity(3)\n",
    "\n",
    "x = numpy.array([3, 4, 5])\n",
    "x_S = numpy.dot(P, x)\n",
    "x_V = numpy.dot(I - P, x)\n",
    "print(x)\n",
    "print(x_S)\n",
    "print(x_V)\n",
    "print(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 $\\vec{q}$.\n",
    "\n",
    "Form the projector $P = \\vec{q} \\vec{q}^\\ast \\in \\mathbb C^{m \\times m}$.  The complement $I - P$ will then include everything **BUT** that direction.  If $||\\vec{q}|| = 1$ we can then simply use $I - \\vec{q} \\vec{q}^\\ast$.  If not we can write the projector in terms of the arbitrary vector $\\vec{a}$ as\n",
    "$$\n",
    "    I - \\frac{\\vec{a} \\vec{a}^\\ast}{||\\vec{a}||^2} = I - \\frac{\\vec{a} \\vec{a}^\\ast}{\\vec{a}^\\ast \\vec{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}(\\vec{a} \\vec{a}^\\ast) = 1$.\n",
    "\n",
    "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": "skip"
    }
   },
   "outputs": [],
   "source": [
    "q = numpy.array([0, 0, 1])\n",
    "P = numpy.outer(q, q.conjugate())\n",
    "P_comp = numpy.identity(3) - P\n",
    "\n",
    "x = numpy.array([3, 4, 5])\n",
    "print(numpy.dot(P, x))\n",
    "print(numpy.dot(P_comp, x))\n",
    "\n",
    "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": [
    "## 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": "subslide"
    }
   },
   "source": [
    "The basic idea is that we want to break up $A$ into its successive spaces spanned by the columns of $A$.  If we have\n",
    "$$\n",
    "    A = \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\vec{a}_1 & \\cdots & \\vec{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix}\n",
    "$$\n",
    "where the columns of $A$ are linearly independent then we want to construct the sequence\n",
    "$$\n",
    "    \\text{span}(\\vec{a}_1) \\subseteq \\text{span}(\\vec{a}_1, \\vec{a}_2) \\subseteq \\text{span}(\\vec{a}_1, \\vec{a}_2, \\vec{a}_3) \\subseteq \\cdots \\subseteq \\text{span}(\\vec{a}_1, \\vec{a}_2, \\ldots , \\vec{a}_n)\n",
    "$$\n",
    "where here $\\text{span}(\\vec{v}_1,\\vec{v}_2,\\ldots,\\vec{v}_m)$ indicates the subspace spanned by the vectors $\\vec{v}_1, \\vec{v}_2, \\ldots, \\vec{v}_m$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "QR factorization attempts to construct a set of orthonormal vectors $\\vec{q}_i$ that span each of the subspaces, i.e. \n",
    "$$\n",
    "    \\text{span}(\\vec{a}_1, \\vec{a}_2, \\ldots, \\vec{a}_j) = \\text{span}(\\vec{q}_1, \\vec{q}_2, \\ldots, \\vec{q}_j).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Before specifying a general procedure we go through the first few steps of the process to get a feel for what needs to be done. In particular we are given a sequence of vectors, $\\vec{a_1}$, $\\vec{a_2}$, $\\ldots$, $\\vec{a_n}$, and generate a sequence of orthonormal vectors  $\\vec{q}_1$, $\\vec{q}_2$, $\\ldots$, and $\\vec{q}_n$ whose spans are the same subspace.\n",
    "\n",
    "First, for $\\text{span}(\\vec{a}_1)$ we can directly use $\\vec{a}_1$ but normalize the vector such that\n",
    "$$\n",
    "    \\vec{q}_1 = \\frac{\\vec{a}_1}{||\\vec{a}_1||}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "2. For $\\text{span}(\\vec{a}_1, \\vec{a}_2)$ we already have $\\vec{q}_1$ so we need to have a vector $\\vec{q}_2$ that is orthogonal to $\\vec{q}_1$, i.e.\n",
    "$$\n",
    "    \\langle \\vec{q}_1, \\vec{q}_2 \\rangle = \\vec{q}_1^\\ast \\vec{q}_2 = 0\n",
    "$$\n",
    "and is again normalized.  We can accomplish this by modifying $\\vec{a}_2$ such that\n",
    "$$\n",
    "    \\vec{q}_2 = \\frac{\\vec{a}_2 - \\langle \\vec{q}_1, \\vec{a}_2\\rangle \\vec{q}_1}{||\\vec{a}_2 - \\langle \\vec{q}_1, \\vec{a}_2\\rangle\\vec{q}_1||}.\n",
    "$$\n",
    "which we can show is orthogonal to $\\vec{q}_1$ as\n",
    "$$\\begin{aligned}\n",
    "    \\langle \\vec{q}_1, \\vec{q}_2 \\rangle &= \\left(\\langle  \\vec{q}_1, \\vec{a}_2 - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle\\vec{q}_1 \\rangle \\right) \\frac{1}{||\\vec{a}_2 - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle \\vec{q}_1||} \\\\\n",
    "    &= \\left(\\langle  \\vec{q}_1, \\vec{a}_2 \\rangle - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle\\langle\\vec{q}_1,\\vec{q}_1\\rangle\\right) \\frac{1}{||\\vec{a}_2 - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle \\vec{q}_1||} \\\\\n",
    "    &= \\left(\\langle  \\vec{q}_1, \\vec{a}_2 \\rangle - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle\\right) \\frac{1}{||\\vec{a}_2 - \\langle  \\vec{q}_1, \\vec{a}_2 \\rangle \\vec{q}_1||} \\\\\n",
    "    &= 0.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "These results suggest then that we may have a matrix factorization that has the following form:\n",
    "$$\n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\vec{a}_1 & \\cdots & \\vec{a}_n \\\\  &  &  \\\\  &  &  \\end{bmatrix} = \n",
    "    \\begin{bmatrix}  &  &  \\\\  &  &  \\\\ \\vec{q}_1 & \\cdots & \\vec{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",
    "    \\vec{a}_1 &= r_{11} \\vec{q}_1 \\\\\n",
    "    \\vec{a}_2 &= r_{22} \\vec{q}_2 + r_{12} \\vec{q}_1 \\\\\n",
    "    \\vec{a}_3 &= r_{33} \\vec{q}_3 + r_{23} \\vec{q}_2 + r_{13} \\vec{q}_1 \\\\\n",
    "    &\\vdots\n",
    "\\end{aligned}$$\n",
    "We can also identify at least the first couple of values of $\\vec{r}$ as\n",
    "$$\\begin{aligned}\n",
    "    r_{11} &= \\left \\Vert \\vec{a}_1 \\right \\Vert \\\\\n",
    "    r_{12} &= \\langle \\vec{q}_1, \\vec{a}_2 \\rangle \\\\\n",
    "    r_{22} &= \\left \\Vert \\vec{a}_2 - r_{12} \\vec{q}_1 \\right \\Vert\n",
    "\\end{aligned}$$\n",
    "It turns out we can generalize this as Gram-Schmidt orthogonalization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gram-Schmidt Orthogonalization\n",
    "\n",
    "As may have been suggestive we can directly construct the arrays $Q$ and $R$ via a process of successive orthogonalization.  We have already shown the first two iterations so lets now consider the $j$th iteration.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want to subtract off the components of the vector $\\vec{a}_j$ in the direction of the $\\vec{q}_i$ vectors where $i=1,\\ldots,j-1$.  This suggests that we define a vector $\\vec{v}_j$ such that\n",
    "$$\\begin{aligned}\n",
    "    \\vec{v}_j &= \\vec{a}_j - \\langle \\vec{q}_1, \\vec{a}_j \\rangle \\vec{q}_1 - \\langle \\vec{q}_2, \\vec{a}_j \\rangle \\vec{q}_2 - \\quad \\cdots \\quad- \\langle \\vec{q}_{j-1}, \\vec{a}_j \\rangle \\vec{q}_{i-1} \\\\\n",
    "    &= \\vec{a}_j - \\sum^{j-1}_{i=1} \\langle \\vec{q}_i, \\vec{a}_j \\rangle \\vec{q}_i.\n",
    "\\end{aligned}$$\n",
    "\n",
    "We also need to normalize $\\vec{v}_j$ which allows to define the $j$th column of $Q$ as\n",
    "$$\n",
    "    \\vec{q}_j = \\frac{\\vec{a}_j - \\sum^{j-1}_{i=1} \\langle \\vec{q}_i, \\vec{a}_j \\rangle \\vec{q}_i}{\\left \\Vert \\vec{a}_j - \\sum^{j-1}_{i=1} \\langle \\vec{q}_i, \\vec{a}_j \\rangle \\vec{q}_i \\right \\Vert}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also discern what the entries of $R$ are as we can write the matrix multiplication as the sequence\n",
    "$$\\begin{aligned}\n",
    "    \\vec{q}_1 &= \\frac{\\vec{a}_1}{r_{11}} \\\\\n",
    "    \\vec{q}_2 &= \\frac{\\vec{a}_2 - r_{12} \\vec{q}_1}{r_{22}} \\\\\n",
    "    \\vec{q}_3 &= \\frac{\\vec{a}_3 - r_{13} \\vec{q}_1 - r_{23} \\vec{q}_2}{r_{33}} \\\\\n",
    "    &\\vdots \\\\\n",
    "    \\vec{q}_n &= \\frac{\\vec{a}_n - \\sum^{n-1}_{i=1} r_{in} \\vec{q}_i}{r_{nn}}\n",
    "\\end{aligned}$$\n",
    "leading us to define\n",
    "$$\n",
    "    r_{ij} = \\left \\{ \\begin{aligned}\n",
    "        &\\langle \\vec{q}_i, \\vec{a}_j \\rangle & &i \\neq j \\\\\n",
    "        &\\left \\Vert \\vec{a}_j - \\sum^{j-1}_{i=1} r_{ij} \\vec{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": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Classical Gram-Schmidt Iteration\n",
    "def classic_GS(A):\n",
    "    m = A.shape[0]\n",
    "    n = A.shape[1]\n",
    "    Q = numpy.empty((m, n))\n",
    "    R = numpy.zeros((n, n))\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\n",
    "\n",
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "Q, R = classic_GS(A)\n",
    "print(A)\n",
    "print(Q)\n",
    "print(numpy.dot(Q.transpose(), Q))\n",
    "print(R)\n",
    "print(numpy.dot(Q, R) - A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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": "subslide"
    }
   },
   "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": [
    "#### Gram-Schmidt in Terms of Projections\n",
    "\n",
    "To start lets rewrite classical Gram-Schmidt as a series of projections:\n",
    "$$\n",
    "    \\vec{q}_1 = \\frac{P_1 \\vec{a}_1}{||P_1 \\vec{a}_1||}, \\quad \\vec{q}_2 = \\frac{P_2 \\vec{a}_2}{||P_2 \\vec{a}_2||}, \\quad \\cdots \\quad \\vec{q}_n = \\frac{P_n \\vec{a}_n}{||P_n \\vec{a}_n||}\n",
    "$$\n",
    "where the $P_i$ are orthogonal projectors onto the $\\vec{q}_1, \\vec{q}_2, \\ldots, \\vec{q}_{i-1}$, in other words the complement of $\\text{span}(\\vec{a}_1, \\vec{a}_2, \\ldots, \\vec{a}_{i-1})$.\n",
    "\n",
    "How should we construct these projectors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We saw before that we can easily construct an orthogonal projector onto the complement of a space by first constructing the projector onto the space itself via\n",
    "$$\n",
    "    \\hat{\\!Q}_{i-1} = \\begin{bmatrix}\n",
    "           &    &       &        \\\\\n",
    "           &    &       &        \\\\\n",
    "        \\vec{q}_1 & \\vec{q}_2 & \\cdots & \\vec{q}_{i-1} \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "Constructing the projection onto the space spanned by $\\vec{q}_1$ through $\\vec{q}_{i-1}$ is then\n",
    "$$\n",
    "    \\hat{\\!P}_{i-1} = \\hat{\\!Q}_{i-1} \\hat{\\!Q}^\\ast_{i-1}\n",
    "$$\n",
    "and therefore the projections in Gram-Schmidt orthogonalization is\n",
    "$$\n",
    "    P_{i} = I - \\hat{\\!P}_{i-1} = I - \\hat{\\!Q}_{i-1} \\hat{\\!Q}^\\ast_{i-1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Modified Gram-Schmidt\n",
    "\n",
    "One problem with the original Gram-Schmidt algorithm is it is not stable numerically.  Instead we can derive a modified method that is more numerically stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Recall that the basic piece of the original algorithm was to take the inner product of $\\vec{a}_j$ and all the relevant $\\vec{q}_i$.  Using the rewritten version of Gram-Schmidt in terms of projections we then have\n",
    "\n",
    "$$\n",
    "    \\vec{v}_i = P_i \\vec{a}_i.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This projection is of rank $m - (i - 1)$ as we know that the resulting $\\vec{v}_i$ 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_i = \\hat{\\!P}_{\\vec{q}_{i-1}} \\hat{\\!P}_{\\vec{q}_{i-2}} \\cdots \\hat{\\!P}_{\\vec{q}_{2}} \\hat{\\!P}_{\\vec{q}_{1}}\n",
    "$$\n",
    "\n",
    "where $\\hat{\\!P}_{\\vec{q}_{i}}$ projects onto the complement of the space spanned by $\\vec{q}_i$.  Note that this performs mathematically the same job as $P_i \\vec{a}_i$ however each of these projectors are of rank $m - 1$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This leads to the following set of calculations:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    1.\\quad \\vec{v}^{(1)}_i &= \\vec{a}_i  \\\\\n",
    "    2.\\quad \\vec{v}^{(2)}_i &= \\hat{\\!P}_{\\vec{q}_1} \\vec{v}_i^{(1)} = \\vec{v}^{(1)}_i - \\vec{q}_1 q_1^\\ast \\vec{v}^{(1)}_i \\\\\n",
    "    3.\\quad \\vec{v}^{(3)}_i &= \\hat{\\!P}_{\\vec{q}_2} \\vec{v}_i^{(2)} = \\vec{v}^{(2)}_i - \\vec{q}_2 \\vec{q}_2^\\ast \\vec{v}^{(2)}_i \\\\\n",
    "    & \\text{  } \\vdots & &\\\\\n",
    "    i.\\quad \\vec{v}^{(i)}_i &= \\hat{\\!P}_{\\vec{q}_{i-1}} \\vec{v}_i^{(i-1)} =  \\vec{v}_i^{(i-1)} - \\vec{q}_{i-1} \\vec{q}_{i-1}^\\ast \\vec{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": "subslide"
    }
   },
   "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": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Modified Gram-Schmidt Iteration\n",
    "def mod_GS(A):\n",
    "    m = A.shape[0]\n",
    "    n = A.shape[1]\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\n",
    "\n",
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "print(A)\n",
    "Q, R = mod_GS(A)\n",
    "print(R)\n",
    "print(numpy.dot(Q.transpose(), Q))\n",
    "print(\"Modified = \")\n",
    "print(numpy.dot(Q, R) - A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Householder Triangularization\n",
    "\n",
    "One way to also interpret Gram-Schmidt orthogonalization is as a series of 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",
    "        \\vec{v}_1 & \\vec{v}_2 & \\vec{v}_3 &  \\cdots & \\vec{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",
    "        \\vec{q}_1 & \\vec{v}_2 &  \\vec{v}_3 & \\cdots & \\vec{v}_n \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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",
    "        \\vec{v}_1 & \\vec{v}_2 &  \\vec{v}_3 & \\cdots & \\vec{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",
    "        \\vec{q}_1 & \\vec{v}_2^{(2)} & \\vec{v}_3^{(2)} & \\cdots & \\vec{v}_n^{(2)} \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The next step can then be placed into the second row:\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "           &    &       &  \\\\\n",
    "           &    &       &  \\\\\n",
    "        \\vec{v}_1 & \\vec{v}_2 & \\vec{v}_3 & \\cdots & \\vec{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",
    "        \\vec{q}_1 & \\vec{q}_2 & \\vec{v}_3^{(3)} & \\cdots & \\vec{v}_n^{(3)} \\\\ \n",
    "           &    &       &        \\\\\n",
    "           &    &       &       \n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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": "subslide"
    }
   },
   "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 it's inverse it's 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": "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{x} & \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x} \\\\\n",
    "    \\end{bmatrix} \\overset{Q_2}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x} \\\\\n",
    "    \\end{bmatrix} \\overset{Q_3}{\\rightarrow}\n",
    "    \\begin{bmatrix}\n",
    "        \\text{x} & \\text{x} & \\text{x} \\\\\n",
    "        0 & \\text{x} & \\text{x} \\\\\n",
    "        0 & 0 & \\text{x} \\\\\n",
    "        0 & 0 & 0\\\\\n",
    "    \\end{bmatrix}\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 matrix of the form\n",
    "\n",
    "$$\n",
    "    Q_k = \\begin{bmatrix} I & 0 \\\\ 0 & F \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $I \\in \\mathbb C^{k-1 \\times k-1}$ identity matrix and $F \\in \\mathbb C^{m - (k - 1) \\times m - (k-1)}$ unitary matrix.  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 $F$. Note that the definition of $Q_k$ implies that\n",
    "$$\n",
    "\\begin{align}\n",
    "Q_k \\vec{x} &= \\begin{bmatrix} I & 0 \\\\ 0 & F \\end{bmatrix} \\vec{x} \n",
    "     &= \\begin{bmatrix} I & 0 \\\\ 0 & F \\end{bmatrix} \n",
    "     \\begin{bmatrix} \\vec{x}_{top} \\\\ \\vec{x}_{bottom} \\end{bmatrix} \\\\\n",
    "     &= \\begin{bmatrix} \\vec{x}_{top} \\\\ F\\vec{x}_{bottom}  \\end{bmatrix}. \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "The vectors $\\vec{x}_{top}$ and $\\vec{x}_{bottom}$ are defined to be\n",
    "$$\n",
    "\\begin{align}\n",
    "     \\vec{x}_{top} &= \\begin{bmatrix} x_1  \\\\  x_2 \\\\ \\vdots \\\\ x_{k-1}\\end{bmatrix} \\\\\n",
    "     \\vec{x}_{bottom} &= \\begin{bmatrix} x_k \\\\ x_{k+1} \\\\ \\vdots \\\\ x_n \\end{bmatrix} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "The task is to construct the matrix $F$ so that after multiplying, $F\\vec{x}_{bottom}$  results in a vector that has some nonzero number in the top position and the rest of the numbers are all zeroes. \n",
    "\n",
    "One way to think of this operation is to consider the transformation that reflects the vector $x$ over the plane $H$ so that $F x = v = ||x|| \\hat{e}_1$:\n",
    "![Householder reflection](./images/householder.png)\n",
    "or mathematically\n",
    "$$\n",
    "    \\vec{x} = \\begin{bmatrix}\n",
    "        x_1 \\\\\n",
    "        x_2 \\\\\n",
    "        \\vdots \\\\\n",
    "        x_{n} \n",
    "     \\end{bmatrix} \\overset{F}{\\rightarrow}\n",
    "    F\\vec{x} = \\begin{bmatrix}\n",
    "        ||\\vec{x}|| \\\\\n",
    "        0 \\\\\n",
    "        \\vdots \\\\\n",
    "        0\n",
    "    \\end{bmatrix} = ||\\vec{x}|| \\hat{e}_1.\n",
    "$$\n",
    "The unit vector $\\hat{e}_1$ is defined to be \n",
    "$$\n",
    "\\hat{x} = \\begin{bmatrix}\n",
    "        1 \\\\\n",
    "        0 \\\\\n",
    "        \\vdots \\\\\n",
    "        0\n",
    "    \\end{bmatrix}.\n",
    "$$   "
   ]
  },
  {
   "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 $v = ||\\vec{x}|| \\hat{e}_1 - \\vec{x}$.  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 - \\frac{\\vec{v} \\vec{v}^\\ast}{\\vec{v}^\\ast \\vec{v}}\\right) \\vec{x} = \\vec{x} - \\vec{x} \\left(\\frac{\\vec{v}^\\ast \\vec{x}}{\\vec{v}^\\ast \\vec{v}} \\right),\n",
    "$$\n",
    "\n",
    "the complement of a projection in the direction of the vector $\\vec{v}$, in other words in the direction of $H$ above.  Since we actually want to transform $\\vec{x}$ to lie in the direction of $\\hat{e}_1$ we need to go twice as far as just the projection onto $H$.  This allows us to identify the matrix $F$ as\n",
    "\n",
    "$$\n",
    "    F = I - 2 \\frac{\\vec{v} \\vec{v}^\\ast}{\\vec{v}^\\ast \\vec{v}}.\n",
    "$$"
   ]
  },
  {
   "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 $\\vec{x}$.  This comes back to having difficulties numerically when the vector $\\vec{x}$ is nearly aligned with $\\hat{e}_1$ and therefore one of the $H$ specification.  By convention the $\\vec{v}$ chosen is defined by\n",
    "\n",
    "$$\n",
    "    v = \\text{sign}(\\vec{x}_1)||\\vec{x}|| \\hat{e}_1 + \\vec{x}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Implementation of Householder QR Factorization\n",
    "def householder_QR(A, verbose=False):\n",
    "    R = A.copy()\n",
    "    v = numpy.empty(A.shape)\n",
    "    m, n = A.shape\n",
    "    for k in range(n):\n",
    "        x = R[k:, k]\n",
    "        e1 = numpy.zeros(x.shape)\n",
    "        e1[0] = 1.0\n",
    "        v[k:, k] = numpy.sign(x[0]) * numpy.linalg.norm(x, ord=2) * e1 + x\n",
    "        v[k:, k] = v[k:, k] / numpy.linalg.norm(v[k:, k], ord=2)\n",
    "        R[k:, k:] -= 2.0 * numpy.dot(numpy.outer(v[k:, k], v[k:, k]), R[k:, k:])\n",
    "\n",
    "    # Form Q\n",
    "    m, n = A.shape\n",
    "    Q = numpy.zeros(A.shape)\n",
    "    for i in range(n):\n",
    "        en = numpy.zeros(m)\n",
    "        en[i] = 1.0\n",
    "        for j in range(n - 1, -1, -1):\n",
    "            en[j:m] -= 2.0 * numpy.dot(numpy.outer(v[j:, j], v[j:, j]), en[j:m])\n",
    "        Q[:, i] = en\n",
    "        \n",
    "    return Q, R\n",
    "\n",
    "A = numpy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)\n",
    "print(\"Matrix A = \")\n",
    "print(A)\n",
    "m, n = A.shape\n",
    "Q, R = householder_QR(A, verbose=False)\n",
    "print(\"Householder (reduced) Q = \")\n",
    "print(Q)\n",
    "print(\"Householder (full) R = \")\n",
    "print(R)\n",
    "\n",
    "print(\"Check to see if factorization worked...\")\n",
    "print(numpy.abs(A - numpy.dot(Q, R[:n, :n])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the above algorithm we do not need to explicitly form the matrix $Q$ to save memory and computation.  If we wanted to for instance solve $A \\vec{x} = \\vec{b}$ we again have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    A \\vec{x} &= \\vec{b} \\\\\n",
    "    Q R \\vec{x} &= \\vec{b} \\\\\n",
    "    R\\vec{x} &= Q^\\ast \\vec{b}.\n",
    "\\end{aligned}$$\n",
    "\n",
    "This requires then the multiplication $Q^\\ast \\vec{b}$.  To do this we just need to recognize that the vectors $\\vec{v}$ we save specify the $Q_k$ so we can form the matrices $F$ and multiply the vector directly with the $\\vec{v}_k$s:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 1:   Random Matrix QR\n",
    "\n",
    "Consider a 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 othorgonal directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "N = 80\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()\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",
    "\n",
    "axes.set_xlabel(\"Index\")\n",
    "axes.set_ylabel(\"$R_{ii}$\")\n",
    "axes.legend(loc=3)\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\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "        0.70000 & 0.70711 \\\\ 0.70001 & 0.70711\n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "Check that the matrix $Q$ is really unitary given this matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%precision 16\n",
    "\n",
    "A = numpy.array([[0.7, 0.70711], [0.70001, 0.70711]])\n",
    "\n",
    "Q, R = classic_GS(A)\n",
    "print(\"Classic: \", numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2)))\n",
    "\n",
    "Q, R = mod_GS(A)\n",
    "print(\"Modified: \", numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2)))\n",
    "\n",
    "Q, R = householder_QR(A)\n",
    "print(\"Housholder:\", numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2)))\n",
    "\n",
    "Q, R = numpy.linalg.qr(A)\n",
    "print(\"Numpy: \", numpy.linalg.norm(numpy.dot(Q.transpose(), Q) - numpy.eye(2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Applications of QR\n",
    "\n",
    "#### Solving $A\\vec{x} = \\vec{b}$ with QR\n",
    "\n",
    "Suppose we want to solve the system $A\\vec{x} = \\vec{b}$ where $A \\in \\mathbb C^{m \\times m}$.  See if you can figure out how to use a QR factorization to help with this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Say we have found the QR factorization of $A$, then\n",
    "$$\\begin{aligned}\n",
    "    A \\vec{x} &= \\vec{b} \\\\\n",
    "    QR \\vec{x} & = \\vec{b} \\\\\n",
    "    Q^\\ast Q R \\vec{x} &= Q^\\ast \\vec{b} \\\\\n",
    "    R \\vec{x} &= Q^\\ast \\vec{b}.\n",
    "\\end{aligned}$$\n",
    "Given that $R$ is upper triangular we can use back-substitution to find $\\vec{b}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Applications: Least Squares Problems\n",
    "\n",
    "Least squares problems have already been introduced but lets consider how our QR factorizations might help us.  As before the least squares problem is characterized by wanting to find the $\\vec{b}$ such that $||\\vec{b} - A\\vec{x}||_2$ is minimized over $\\vec{x} \\in \\mathbb C^n$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we are using the $\\ell_2$ norm and know this is equivalent to the Euclidean norm we know that there is a geometric interpretation to this goal, find the vector $\\vec{x}$ that gives the minimum distance between the vector $\\vec{b}$ and $A \\vec{x}$.  This can be interpreted as a projection:\n",
    "![Least-Squares Projection](./images/lsq_projection.png)\n",
    "where\n",
    "$$\n",
    "    \\vec{r} = \\vec{b} - A\\vec{x}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "    \\vec{y} = A\\vec{x} = P\\vec{b}.\n",
    "$$\n",
    "The vector $\\vec{r}$ is called the residual (and the thing we are trying to minimize).  $P$ represents the orthogonal projector onto the $\\text{range}(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "QR factorization plays a role similar to the ideas we saw from Householder triangularization.  Define the orthogonal projector $P = Q Q^\\ast$ based on the reduced QR factorization of $A$.  We know then that the projector projects onto the span of column space of $A$, $\\text{span}(A)$.  Using this QR factorization we know that the least-squares formulation then becomes\n",
    "$$\\begin{aligned}\n",
    "    A^\\ast A \\vec{x} &= A^\\ast \\vec{b} \\\\\n",
    "    R^\\ast Q^\\ast Q R \\vec{x} &= R^\\ast Q^\\ast \\vec{b} \\\\\n",
    "    R^\\ast R \\vec{x} & = R^\\ast Q^\\ast \\vec{b} \\\\\n",
    "    R \\vec{x} & = Q^\\ast \\vec{b}\n",
    "\\end{aligned}$$\n",
    "reducing the least-squares calculation to one of finding the QR factorization and backwards substitution."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}