Gram-Schmidt.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:5c17b9d84f37fe41e17c63ca91fd1034a55477548bff302c05560cfc216cab70"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "%matplotlib inline\n",
      "np.set_printoptions(precision=2) # So printed matrices aren't huge\n",
      "np.set_printoptions(suppress=True) # Don't use scientific notation"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "QR Factorization by the Gram-Schmidt algorithm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's look at the accuracy of the classical and modified Gram-Schmidt algorithms (see, e.g. Trefethen & Bau lecture 7).  Here is an implementation of the classical Gram-Schmidt algorithm; notice that it uses more storage than necessary."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def clgs(A):\n",
      "    \"QR factorization by the classical Gram-Schmidt algorithm.\"\n",
      "    n = A.shape[1] # Number of columns of A\n",
      "    R = np.zeros([n,n])\n",
      "    V = np.zeros(A.shape)\n",
      "    Q = np.zeros(A.shape)\n",
      "    \n",
      "    for j in range(n):      # loop over columns of R\n",
      "        V[:,j] = A[:,j]\n",
      "        for i in range(j):  # Make jth column of V orthogonal to ith column of Q\n",
      "            R[i,j] = np.dot(Q[:,i].T,A[:,j])\n",
      "            V[:,j] = V[:,j] - R[i,j]*Q[:,i]\n",
      "        \n",
      "        R[j,j] = np.linalg.norm(V[:,j],2)\n",
      "        Q[:,j] = V[:,j]/R[j,j]             # Normalize jth column of V to get jth column of Q\n",
      "    \n",
      "    return Q, R"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we'll generate a random $4 \\times 4$ matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 4\n",
      "A = np.random.rand(n,n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's factor it:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q, R = clgs(A)\n",
      "print 'Q = \\n', Q\n",
      "print 'R = \\n', R"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Q = \n",
        "[[ 0.54 -0.75  0.38  0.04]\n",
        " [ 0.58  0.07 -0.73  0.35]\n",
        " [ 0.38  0.59  0.57  0.42]\n",
        " [ 0.47  0.29 -0.   -0.83]]\n",
        "R = \n",
        "[[ 1.28  1.09  1.18  0.77]\n",
        " [ 0.    0.44  0.07  0.18]\n",
        " [ 0.    0.    0.37  0.13]\n",
        " [ 0.    0.    0.    0.32]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Clearly $R$ is upper triangular.  Is $Q$ unitary?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.max(np.abs(np.dot(Q.T,Q)-np.eye(n)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.22514845491e-15\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not bad.  We can also check our factorization:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.max(np.abs(np.dot(Q,R)-A))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.11022302463e-16\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That looks pretty good. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Instability of classical Gram-Schmidt"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's try it on a different matrix.  The function below sets up a Hilbert matrix; these matrices are notoriously ill-conditioned (i.e., close to singular)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Hilbert(n):\n",
      "    \"\"\"Return the n x n Hilbert matrix.\"\"\"\n",
      "    A = np.zeros([n,n])\n",
      "    for i in range(n):\n",
      "        for j in range(n):\n",
      "            A[i,j] = 1./(i+j+1)\n",
      "    return A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's set up and factor a modest-sized Hilbert matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 10\n",
      "H = Hilbert(n)\n",
      "print H"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.    0.5   0.33  0.25  0.2   0.17  0.14  0.12  0.11  0.1 ]\n",
        " [ 0.5   0.33  0.25  0.2   0.17  0.14  0.12  0.11  0.1   0.09]\n",
        " [ 0.33  0.25  0.2   0.17  0.14  0.12  0.11  0.1   0.09  0.08]\n",
        " [ 0.25  0.2   0.17  0.14  0.12  0.11  0.1   0.09  0.08  0.08]\n",
        " [ 0.2   0.17  0.14  0.12  0.11  0.1   0.09  0.08  0.08  0.07]\n",
        " [ 0.17  0.14  0.12  0.11  0.1   0.09  0.08  0.08  0.07  0.07]\n",
        " [ 0.14  0.12  0.11  0.1   0.09  0.08  0.08  0.07  0.07  0.06]\n",
        " [ 0.12  0.11  0.1   0.09  0.08  0.08  0.07  0.07  0.06  0.06]\n",
        " [ 0.11  0.1   0.09  0.08  0.08  0.07  0.07  0.06  0.06  0.06]\n",
        " [ 0.1   0.09  0.08  0.08  0.07  0.07  0.06  0.06  0.06  0.05]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q, R = clgs(H)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How good is our factorization?  Let's check."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.max(np.abs(np.dot(Q,R)-H))\n",
      "print np.max(np.abs(np.dot(Q.T,Q)-np.eye(n)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "5.55111512313e-17\n",
        "0.999965967367\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Whoa!  We still have an accurate factorization of $A$, but $Q$ isn't even close to orthogonal:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.dot(Q.T,Q)-np.eye(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.    0.   -0.    0.   -0.    0.    0.    0.    0.    0.  ]\n",
        " [ 0.    0.   -0.    0.   -0.    0.   -0.   -0.   -0.   -0.  ]\n",
        " [-0.   -0.    0.    0.   -0.    0.   -0.   -0.   -0.   -0.  ]\n",
        " [ 0.    0.    0.   -0.   -0.    0.   -0.   -0.   -0.   -0.  ]\n",
        " [-0.   -0.   -0.   -0.    0.    0.   -0.   -0.   -0.   -0.  ]\n",
        " [ 0.    0.    0.    0.    0.   -0.   -0.05 -0.04 -0.03 -0.03]\n",
        " [ 0.   -0.   -0.   -0.   -0.   -0.05  0.    1.    1.    1.  ]\n",
        " [ 0.   -0.   -0.   -0.   -0.   -0.04  1.   -0.    1.    1.  ]\n",
        " [ 0.   -0.   -0.   -0.   -0.   -0.03  1.    1.    0.    1.  ]\n",
        " [ 0.   -0.   -0.   -0.   -0.   -0.03  1.    1.    1.    0.  ]]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Modified Gram-Schmidt (to the rescue!)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's try the modified algorithm on the same Hilbert matrix.  Here is the modified algorithm (again, using more storage than necessary):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mgs(A):\n",
      "    \"QR factorization by the modified Gram-Schmidt algorithm.\"\n",
      "    n = A.shape[1]\n",
      "    R = np.zeros([n,n])\n",
      "    V = np.zeros(A.shape)\n",
      "    Q = np.zeros(A.shape)\n",
      "    for i in range(n):\n",
      "        V[:,i] = A[:,i]\n",
      "    for i in range(n):\n",
      "        R[i,i] = np.linalg.norm(V[:,i],2)\n",
      "        Q[:,i] = V[:,i]/R[i,i]\n",
      "        for j in range(i,n):\n",
      "            R[i,j] = np.dot(Q[:,i].T,V[:,j])\n",
      "            V[:,j] = V[:,j] - R[i,j]*Q[:,i]\n",
      "    return Q, R"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q2, R2 = mgs(H)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How good is our factorization?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.max(abs(np.dot(Q2,R2)-H))\n",
      "print np.max(abs(np.dot(Q2.T,Q2)-np.eye(n)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.11022302463e-16\n",
        "0.000124230592673\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The orthogonality of $Q$ is far from perfect, but much better than what we got from the classical algorithm."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can also see the inaccuracy of cGS if we look at the diagonal entries of $R$, which approximate the singular values of $A$ (see Trefethen & Bau lecture 9):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "plt.figure(figsize=(8,4))\n",
      "plt.semilogy(np.diag(R),'r',linewidth=3)\n",
      "plt.hold(True)\n",
      "plt.semilogy(np.diag(R2),'ok')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10b57c850>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAEDCAYAAAAY+lwJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lPW99/H3lyCJYF3A2ipih0axUG2VqmhbZFrxSWiq\neFxKcanLUbuZ0N3l9FGec55Wq20VYq0bcEARLNW6xSIuZ6ytVcFWixWspMay41qLmMjyPX/8JmQy\nmUAmmczMnfm8ritXuZe57+9cFT753fdvMXdHREREil+/QhcgIiIiXaPQFhERiQiFtoiISEQotEVE\nRCJCoS0iIhIRCm0REZGIUGiLiIhEhEJbREQkIvIa2mY23MxuNbMF+byviIhIX5DX0Hb3V9z9/Hze\nU0REpK/ocWib2UwzW29mS9P2V5vZcjN72cwu7ul9RERESl0uWtqzgOrUHWZWBlyf3D8KmGxmI3Nw\nLxERkZLV49B29yeAt9J2HwWscPcmd98MzAcmmtlgM7sROEytbxERkez076XrDgVWpmyvAsa4+5vA\n13b0QTPTsmMiIlJy3N12dk5vdURT8IqIiORYb4X2amBYyvYwQmu7y9w98j9XXHFFwWvQd+g736Mv\nfAd9j+L66QvfoS98j2z0VmgvAQ4ys5iZDQAmAff10r1ERERKQi6GfM0DngRGmNlKMzvX3bcAFwEP\nAS8Cd7r7sp7eS0REpJT1uCOau0/uZP9vgd/29PpRFo/HC11Cj/WF7wB943v0he8A+h7FpC98B+g7\n36MrLNvn6b2ttfd4sdUlIiLSG8xCp3EvYO9xERERyTGFtoiISEQotEVERCJCoS0iIhIRCm0REZGI\n6K25xzMys0HADUALkHD3O/J5fxERkSjLd0v7ZOBX7n4hcOKOTqyqqqKhoSE/VYmIiERAj1vaZjYT\nqAE2uPuhKfurgeuAMuBWd/8JYfWv55OnbN3RdRctWkRjYyMANTU1PS1TREQk8nLR0p4FVKfuMLMy\n4Prk/lHAZDMbSVg0pHUhkZ3eu7Gxkfof/hDefTcHZYqIiERbj0Pb3Z8A3krbfRSwwt2b3H0zMB+Y\nCNwNnGJmN9DFBUSan3sO9t4bTjgBZsyADRt6WrKIiEgk9VZHtKHAypTtVcAYd98EnJfNhSoAmpvh\ngQfCjxl8+tNw0kkwcSIcdFDuqhYRESlivRXaOZk4fLdddmHwrruSeOcd4tuv7PCHP4Sf738fRo0K\n4X3SSXDEEdBPo9hERKS4JRIJEokETU1NWX0uJwuGmFkMuL+1I5qZHQ1Mdffq5PalwLZkZ7SdXcsh\n9B6vra0NndBWrIB77w0/f/gDbNuW+cP77QcnnhgC/HOfgwEDevzdREREelM2C4b0Vmj3B14CjgPW\nAM8Ak7uypvZOV/l67bXwmPyee2DRovDoPJMPfAC+8IUQ4BMmwB57ZPu1REREel1eQ9vM5gHjgCHA\nBuByd59lZhNoG/I1w92v7OL1ur4057vvwsMPhwC//354883M5+2yS2h5T5wYfoYO7UopIiIivS7v\nLe1c6vZ62lu2hEfn99wTHqO/8krn5x55ZNt78FGjQuc2ERGRAijN0E7lDkuXhvC+5x740586P7ey\nMoT3SSfBMcdAWVn37ysiIpIlhXa6f/wD7rsvBHgiAVs7mYztgx8M48EnToTjj4ddd81dDSIiIhko\ntHfkrbfgwQdDgC9cCBs3Zj5v4ECoqgoB/sUvwpAhvVOPiIiUNIV2VzU3w2OPhQC/7z5Yvz7zef36\nwdixbRO6DB/e+7WJiEhJUGh3x7Zt8PTTbe/BX3qp83M/8Ym2AD/8cHVkExGRblNo58Ly5W0B/tRT\nnZ93wAFtQ8mOPTYMLxMREemiog5tMxsO/Aewh7ufluF4cYR2qrVrwzjwe++FRx6B99/PfN6ee0JN\nDQ3DhjH9mWdo2bqV8vJy6urqtLyoiIhkVNShvf3GZgsiE9qp/vWv0IHt3nvDzGz//Ge7ww3AFKAx\nZV9lZSXTpk1TcIuISAfZhHa3V9cws5lmtt7Mlqbtrzaz5Wb2spld3N3rF60PfABOOw1uvz1MqfrI\nI3DRRbD//gBMp31gQ3Jd8G98A559NowhFxER6YaeLIk1C6hO3WFmZcD1yf2jgMlmNtLMzjKza81s\nvx7cr/jssgscdxzU14ex4M8+S8sBB2Q8tfkf/wirkB16KFxzDaxZk+diRUQk6rod2u7+BPBW2u6j\ngBXu3uTum4H5wER3v83dv+3ua8xssJndCBzWp1riZjB6NOUf+1jGwxWtf/jrX+EHP4Bhw6C6GubN\ng02b8lamiIhEV64Xnx4KrEzZXpXct527v+nuX3P3g7qyVGfU1NXVUVlZ2W5f5b77UhuPhwlbWm3b\nBg89BKefDvvuC+efD088ocfnIiLSqf45vl7OEicejxOLxYjFYsTjceLxeK4u3ataO5vV19fT3NxM\nRUVF27rg//oX3HUXzJ4dplNt9c47MGNG+Bk+HL7ylfDz0Y8W5kuIiEivSiQSJBIJmpqasvpcj3qP\nZ1hH+2hgqrtXJ7cvBbZl06Iu+t7jufLqq3DbbSHAV6zIfM5nPwtnnx06vmk9cBGRPilvQ74yhHZ/\n4CXgOGAN8Aww2d2XZXHN0gjtVu7wxz/CnDkwf36HIWQAVFSEGdjOPhvGj4f+uX5AIiIihZKX0Daz\necA4YAiwAbjc3WeZ2QTgOqAMmOHuV2Z53dIK7VTNzWESl9mzw1jwTKuR7bsvnHFGCPBDDsl/jSIi\nklORmFylMyUd2qnWrQs9y2fPhuefz3zO6NHh3ffpp4dlRUVEJHIU2n3N88+Hx+dz52Zeiax/f5gw\nIbS+v/hFKC/Pf40iItItCu2+asuWMExszpwwjWpLS8dz9toLvvzlEOBHHaUVyEREipxCuxS89RYs\nWBAenz/5ZOZzDj44PD4/66wwmYuIiBQdhXapefnlMHxszpwwlCydGXz+8yHATz4Zdtst/zWKiEhG\nCu1StW0b/O53IbwXLICNGzueM2gQnHpqCPB4HPrlelI8ERHJhkJb4N134Te/CQH+yCOZp0cdNiw8\nOj/7bBgxIv81ioiIQlvSrFoVep7Png3LOpnn5uijQ+t70iQYPDi/9YmIlLCiDm0zmwjUALsTJl95\nOO24Qru3uMOSJaH1fccd8OabHc8ZMABOPJGGkSOZ/tRTtLz/PuXl5dTV1W2fV11ERHKnqEN7+43N\n9gR+6u7np+1XaOfD++9DQ0MI8AceCMPJkhqAKUBjyumVlZVMmzZNwS0ikmP5msZ0JqHFvKF17vHk\n/mrapjG9tbPFQszsp8Dt7v5c2n6Fdr69/nqYfW3OHFiyhCpgUYbTqg49lIVPPw277prvCkVE+qxs\nQrsnXYdnAdVpNy4Drk/uHwVMNrORZnaWmV1rZvtZ8BPgt+mBLQWy995QWwuLF8MLL9DSyZju5qVL\nYehQ+Na34MUX81ykiIh0O7Td/QngrbTdRwEr3L3J3TcD84GJ7n6bu3/b3dcAtYRVwE41s6929/7S\nSz7+ccpHjsx4qALCpC7TpsHHPw5jx4bx4e+9l9cSRURKVa7XeBwKrEzZXgWMST3B3acD03d2oXg8\nTiwWIxaLEY/HicfjOS1UOldXV0djYyONjW1vtSsHD6a2f3/YsKHtxN//PvxMmRJ6nl94IYwaVYCK\nRUSiJZFIkEgkaGpqyupzuV5P+xSg2t0vSG6fCYxx99osrql32kWgoaGB+vp6mpubqaiooLa2lpoJ\nE+DRR+Gmm8Lc5ymd17b77GdDeJ96qt59i4h0Qd56j2cI7aOBqe5endy+FNjWWWe0Tq6p0I6Cdetg\n1iy45RZ45ZWOx/faK0zacsEFan2LiOxAIUO7P/AS4Z31GuAZYLK7dzKjR8ZrKrSjZNu2MOPazTfv\nuPX91a/CKaeo9S0ikiZfQ77mAeOAIcAG4HJ3n2VmE2gb8jXD3a/M8roK7ajqauv7wguhk85uIiKl\nJhKTq3RGod0HdKX1PXZs27vvior81ygiUiQU2lI81q6F//5vtb5FRDqh0Jbi09r6vukmuO8+tb5F\nRJIU2lLcdtb6Hjy4bdy3Wt8i0scptCUaUlvf994LW7d2PGfs2Lae52p9i0gfpNCW6Fm7tq3neaYZ\ngtT6FpE+SqEt0bVtGzz8cFvP80yt72OPDeGt1reI9AFFHdpm9jHCcs1DgIfcfUbacYW2BF1pfbfO\nuqbWt4hEVFGH9vYbm/UD5rv7l9L2K7SlPbW+RaQPy8t62mY208zWm9nStP3VZrbczF42s4s7+ewJ\nQANh6U6RHevXD6qq4K67YOVK+NGPIBZrf87vfgdnnhnW+/7Od2D58oKUKiLSm3oyjelYYCMwJ2Xu\n8TLC3OPjgdXAYmAycAQwGrgmuaZ26zXudfeJaddVS1t2rrX13Truu5PWd8ORRzL9+edp2byZ8vJy\n6urqqKmpyX+9IiKdKOSCIccAV6Ss8nVJspCrUj4zDjgZqACWuft1addUaEt21q6FmTPDu+9XX92+\nu4HQeaIx5dTKykqmTZum4BaRolHI0D4VqMrFetrjxo0jFosRi8WIx+PE4/Fu1yklIq31XbV1K4sy\nnFZ1xBEsfOYZsJ3+/RAR6RWJRIJEIkFTUxOzZ88GChPapwDVuQhttbSlR9auJf7pT/N4hl7n44DE\nyJFw0UVh7Pduu+W9PBGRVnnpiNaJ1cCwlO1hwKoc30Nk5/bdl/IRIzIeqgBYtgy++c3Qce3b34YV\nK/JanohId+Q6tJcAB5lZzMwGAJOA+3J8D5Euqauro7Kyst2+yt13p3bXXdt2vPMOXHcdjBgBNTWw\ncGF4zC4iUoR60nt8HuFJ4xBgA3C5u88yswnAdUAZMMPdr8zyuno8LjnT0NBAfX09zc3NVFRUUFtb\nS83YsTB7Nlx/Pfztbx0/NGJEeHR+9tmw++75L1pESkokJlfpjEJb8qa141p9PTz4IKT/N7fbbnDO\nOSHADz64ICWKSN+n0BbJ1ooV8ItfhKFj77zT8XhVFdTWwoQJYbIXEZEcUWiLdNfGjXDbbaH1vWxZ\nx+MHHhg6sJ17LuyxR/7rE5E+R6Et0lPu8OijIbzvv7/jo/NBg8JwsYsuglGjClOjiPQJCm2RXPr7\n3+GGG2DGDHj77Y7Hx48Pj85raqCsLP/1iUikKbRFesO778LcuTB9Ovz1rx2PDx8eHp2fdx7stVf+\n6xORSFJoi/Qmd0gkwqPze+/tOK574MCw4lhtLRxySEFKFJHoKPrQNrNBQAKY6u4NaccU2hIdr74a\nHp3feiu8+WbH45/7XAjvE06A/v3zX5+IFL0ohPb/A/5FWOVLoS3Rt2kTzJsXHp3/5S8djx9wQHh0\n/u//DkOG5L8+ESlaeZl73Mxmmtl6M1uatr/azJab2ctmdnGGzx0PvAi81t17ixSdgQNDID/3HDz+\nOJx6avtOaf/4B1x8Mey/P5x/Pjz/fOFqFZHI6sk0pmOBjcCclFW+yoCXgPGExUMWA5OBI4DRwDXA\nN4BBwCjgPeDfPKUItbSlz1i5Em68EW6+GV5/vePxY48Nj85POkmPzkVKWCHX0z4GuMLdq5PblyQL\nuSrDZ88GXnP3B9P2K7Slb2luhvnzw6PzP/+54/H994evfx0uuAA++MH81yciBVXI0D4VqMrFetrj\nxo0jFosRi8WIx+PE4/Fu1ylSFNzhySdDr/O77oItW9ofLy+HyZND63v06MLUKCJ5kUgkSCQSNDU1\nMXv2bKAwoX0KUJ2L0FZLW/q01avhppvCz4YNHY9/5jMhvE8+GXbZJf/1iUjeFLKlfTRhGFfr4/FL\ngW3u/pMsrqnQltLR0gK/+lVofS9e3PH4fvvR8LnPMX3VKlqA8vJy6urqqKmpyXupItI7Chna/Qkd\n0Y4D1gDPAJPdPcPKC51eU6Etpenpp8N77wULYPNmABqAKUBjymmVlZVMmzZNwS3SR+QltM1sHjAO\nGAJsAC5391lmNgG4DigDZrj7lVleV6EtpW3duvDY/MYbqVq3jkUZTqn61KdYuHgx2E7/jotIkSv6\nyVV2RKEtkvT++8Q/+UkeX768w6FxQOLjH4fvfhdOPz10YhORSMrL5Coi0ssGDKD8gAMyHqqAsGjJ\needBLAZXXglvvZXP6kSkABTaIkWsrq6OysrKdvsq99iD2tSW9bp1cNllMGwYTJkCr7yS5ypFJF/0\neFykyDU0NFBfX09zczMVFRXU1tZS8+lPh5nWpk+HNWvaf6BfvzCN6ve+B0ceWZiiRaTL9E5bpFS8\n/36Ybe2nP4WlSzseP/bYEN41NSHMRaToKLRFSo07PPxwCO+HH+54/OCDQ6e1s86Cior81ycinVJo\ni5Sy55+Hn/0sLBWaPlXqPvvARReFuc733rsw9YlIO0Ud2mYWB/4LeAGY7+6Ppx1XaIvkwqpV4Z33\nTTfBO++0P7brrnDuufDtb8OBBxamPhEBin/I1zbgX0A5sKoA9xcpDfvvD1dfHZYI/dnPQu/yVu+9\nBzfcACNGhPnNn3yycHWKSJf1ZEa0mUANsKF1GtPk/mraZkS7NX3ecTMzd3cz2wf4ubufmXZcLW2R\n3rB5c5gi9ac/zbxE6DHHhE5rEydCWVn+6xMpUflqac8CqtNuXAZcn9w/CphsZiPN7Cwzu9bM9vO2\nNH6b0NoWkXzYZZcwe9qzz8Kjj8KECe2P//GPcMopodPaDTfApk2FqVNEOpXrBUOOAa5IWeXrEgB3\nvyrlM/8GVAF7Aje4++/SrqmWtki+vPAC/PzncPvt2xcp2W7IEPjGN+Cb34QPfagw9YmUgEK+0x4K\nrEzZXpXct527/8bdv+buX04PbBHJs0MOgZkz4dVX4dJLYc8924698Qb813/BRz4CF14IGeZAF5H8\nynVo56x5HI/HOeecc5g6dSqJRCJXlxWRTPbdF37849Bpbdq0MJ95q5YWuOUWGDkSTjwRfve7MC5c\nRLotkUgwdepUzjnnnKw+l+vH40cDU1Mej18KbEvvjLaTa+rxuEihbdkCd98dOq0tXtzx+JFHhk5r\nJ58M/fvnvz6RPqSQj8eXAAeZWczMBgCTgPtyfA8R6W39+8OXvgRPPx1a1iee2P744sUwaRIcdFAY\nC75xY2HqFCkxPRnyNY+wrO8QYANwubvPMrMJtA35muHuV2Z5XbW0RYrR8uVw7bUwe3Z4ZJ5qzz3D\nLGu1teFRu4h0WVHPiLYzCm2RIrdhA/ziF+HnjTfaH9tlFzjzTPjOd0InNxHZKYW2iPS+TZtCq/vn\nP4cVKzoer64O770//3mwnf5bJFKyFNoikj9bt8J994VOa5mmQz3sMPje92jYbTem33ADLS0tlJeX\nU1dXR01NTf7rFSkyCm0RKYwnnwzznP/mN+2GhTUAU8rKaNy6dfu+yspKpk2bpuCWvsE9TFDU3Jz1\nj02ZkryEQltECmHFitBpbdYseO89qoBFGU6rqqpi4cKF+a5O+qIehGbOfrrbsXv7V1Boi0ghvf46\n/PKXxP/zP3k8fW1vYNzQoSSeeQb2268AxUnRe/99WLcO1q6FNWvaflK3166Fd9/tUWgWmkJbRIpK\n1fHHs+iRRzruBxYOGBDW9r74Yhg+PP/FSf5t3gzr13cM4/RAfu21Qleanf79oaIi6x+7/nqgSEPb\nwhv3/w98AFji7nPSjiu0RfqYhoYGpkyZQmNj4/Z9lcA0wvq+QFgO9PTTwxzoI0cWoErpsS1bwpDA\nnYXxhg290yruZmjm5Ke8vNuzAxZ1R7TkKl8TgdeBB939sbTjCm2RPqihoYH6+nqam5upqKig9uij\nqXnoIXjqqfYnmoXpUf/jP+DwwwtTrLS3dWto9e4sjNevh23bcndfs7DC3H77tf3su2/H7d1371Fo\nFlpeQtvMZhJ+Sd7QOvd4cn81bTOi3Zo+77iZXQy86e63mNkCdz8t7bhCW6RUuMP//A/86Efw2GMd\nj0+YEML7M5/Jf22lYNu20O9gZ2G8bl0I7lzaZ5/Og7h134c+FNkgzka+QnsssBGYk7JgSBnwEjAe\nWA0sBiYDRwCjgWuAzwHvu/sCM7vT3SelXVehLVKKnnoqhPcDD3Q8Nm5cCO/x4zVRS0+tWgXz5sH8\n+fCXv4RH2rm09947D+MPfzjMnidAHh+PZ1jl6xjgipRVvi5JFnJVymd2BeqBTcAyd/9l2jUV2iKl\n7PnnwzKhCxZ0fO955JEhvE84Afrler2jPuztt+Guu2DuXEgkuvc+efDgroVxeXnOy+/rChnapwJV\n7n5BcvtMYIy712ZxTYW2iMBLL8FVV8Htt3dsDR5yCFx2WViJrKysMPUVu+ZmePDBENQPPBCGT2Wy\n556dB3Hr/n33DZ2tpFdkE9q5flmQs6SNx+PEYjFisRjxeJx4PJ6rS4tIFBx8cJicZepUuPpqmDGj\nbXWxF14IPc0vvxwuuQTOOgsGDChouUVh2zZ4/PEQ1L/+Nfzznx3P6dcvzAd/xhlw0kkhtCXvEokE\niUSCpqamrD6X65b20cDUlMfjlwLb0juj7eSaammLSEdr14bFSX75yzCZRqphw+D734fzz4dddy1M\nfYXiHl4pzJ0b3lWvXp35vNGjQ1B/+cuazKbIFPLxeH9CR7TjgDXAM8Bkd1+WxTUV2iLSuTfegOnT\nw8/bb7c/ts8+YVnQr389DAPqy159Fe64I4T1X/+a+ZyPfjQ8kTjjDPjYx/Jbn3RZvnqPzwPGAUOA\nDcDl7j7LzCbQNuRrhrtfmeV1FdoisnPvvBNa3T//eZisI9Wee0JdXfgZMqQw9fWGN94IHfTmzoXf\n/z7zOXvvDZMmhaA++mj1to+Aop5cZWcU2iKSlU2bwvvuq68Ow5lSDRoUWt3f/W7o2RxFmzbB/feH\noF64MEwBmm7gwPB++owz4PjjNZwqYhTaIlJ63n8fbrst9DhfsaL9sfLy8L77+9+Hj3ykMPVlY8uW\nMNnM3Llw992wcWPHc8rKQkCfeSZMnAi77Zb/OiUnFNoiUrq2bAmPkH/849DLPFX//iHkLr0URowo\nTH2dcYdnnw1BPX9+mIUskzFjQot60qTwDl8iT6EtIrJtW3is/KMfweLF7Y+ZhTHel10Gn/hEYepr\n1dgYgnruXPjb3zKfM2JECOrTT4cDD8xvfdLrFNoiIq3c4ZFHQng//njH4yecEGZZGzMmfzVt2AB3\n3hmC+umnM5/zoQ/B5MkhrD/1KXUo68MU2iIimfzhDyG8f/vbjseOOy6EdzzeOwG5cSPce2+Y4e3h\nhzMvwLHbbmGFszPOCBOglMBiGaLQFhHZsT/9KbzzvvvujvNwH3NMCO8vfKHn4b15cwjouXPhnntC\nT/B0/fuH1czOOCO0+gcO7Nk9JXKKOrTN7LPAGYQpVEe5+2fSjiu0RSQ/Xnwx9Da/446OLd/DDgvv\nvE8+Obv5zd3DimVz58KvfhXWoc7ks58NQX3aaX1rLLlkrahDe/uNzSYC+7j7LWn7Fdoikl9//3sY\n5z1rVseFNQ4+OPQ2P/30HY9/Xr48BPUdd4TrZTJqVFuHslgsZ+VLtOVrRrSZQA2woXUa0+T+atpm\nRLu1s3nHzexO4Dx3fzdtv0JbRApj9Wr42c/gxhvhvffaH4vF4Ac/oOHDH2b6jTfS0tJCOVB34IHU\nPPdcGK6VydChbR3KPvlJdSiTDvIV2mOBjcCclLnHywhzj48HVgOLgcnAEcBo4Bp3X2NmBwA/dPcL\nM1xXoS0ihfXaa3DddXD99WG61KQGYEpZGY0pj9IrgWmEFsx2e+wBp54agvrYY7V8qOxQIRcMOQa4\nImWVr0uShVyV9rmpwEJ3fyrDNRXaIlIc/vlP+MUv4Npr4fXXqQIWZTitClg4YADU1ITJW77wBa0/\nLV2WTWj3y/G9hwIrU7ZXJfe14+5TMwW2iEhR2WOP0BmtqQmuvZaWTt5pN48YEWYwu/vu0HFNgS29\nJNeDAHPWPI7H48RiMWKxGPF4nHg8nqtLi4hkZ9Ag+Na3KG9oCBO1pKkYPhz22qsAhUlUJRIJEokE\nTU1NWX0u16G9GhiWsj2M0NrOWiKRyEU9IiI5U/etb9H4yis0NjZu31dZWUltbW0Bq5IoSm2Mzp49\nu8ufy3VoLwEOSr7rXgNMInREExGJvJqa0N2svr6e5uZmKioqqK2t3b5fpLf1pPf4PGAcMATYAFzu\n7rPMbAJtQ75muPuVWV5XHdFERKRkRGJylc4otEVEpJQUsve4iIiI9BKFtoiISEQotEVERCJCoS0i\nIhIRCm0REZGIUGiLiIhEhEJbREQkInI9I9pOmdn+wHTgLeBvna23LSIiIu0VoqV9KHCXu/87cHgB\n7i8iIhJJ3Q5tM5tpZuvNbGna/mozW25mL5vZxRk++iRwoZk9Cizs7v1FRERKTU/mHh8LbATmuPuh\nyX1lwEvAeMKKX4sJC4YcAYwGrgG+BDzr7k+Y2QJ3Py3tuprGVERESkY205h2+512MnRjabuPAla4\ne1OykPnARHe/Crgtue8x4HIzOx14pbv3FxERKTW57og2FFiZsr0KGJN6grv/BTh1ZxeKx+PEYjFi\nsVi7dUdFRESiLpFIkEgkaGpqyupzPVrlK9nSvj/l8fgpQLW7X5DcPhMY4+5dXiFej8dFRKSUFHKV\nr9XAsJTtYYTWtoiIiPRQrkN7CXCQmcXMbAAwCbgvx/cQEREpST0Z8jWPMHxrhJmtNLNz3X0LcBHw\nEPAicKe7L8tNqSIiIqWtR++0e4PeaYuI9K6GhgamT59OS0sL5eXl1NXVUVNTU+iySlZehnyJiEj0\nNDQ0MGXKFBobG7fva/2zgrv4acEQEZESMn369HaBDSG06+vrC1SRZEOhLSJSQlpaWjLub25uznMl\n0h0KbRGRElJeXp5xf0VFRZ4rke5QaIuIlJC6ujoqKyvb7ausrKS2tstzYEkBqSOaiEgJae1sVl9f\nT3NzMxUVFdTW1qoTWkTkfciXmY0CrgDeAB5197vSjmvIl4iIlIxCTmPaFdVAvbt/A/hKAe4vIiIS\nST2ZEW2mma03s6Vp+6vNbLmZvWxmF2f46G3Al83samBId+8vIiJSarr9eNzMxgIbgTkpq3yVAS8B\n4wmLhyx5mccYAAAKcklEQVQGJgNHAKOBa9x9Tcq5d7n7SWnX1eNxEREpGXmZEc3dn0guzZnqKGCF\nuzclC5kPTHT3qwgtbMzsI8BlwCDg6u7eX0REpNTkuvf4UGBlyvYqYEzqCe7+KvDVnV0oHo8Ti8WI\nxWLE43Hi8XhOCxURESmURCJBIpGgqakpq8/1qPd4sqV9f8rj8VOAane/ILl9JjDG3bs8AFCPx0VE\npJQUsvf4amBYyvYwQmtbREREeijXob0EOMjMYmY2AJgE3Jfje4iIiJSkngz5mgc8CYwws5Vmdq67\nbwEuAh4CXgTudPdluSlVRESktOV9RrSd0TttEREpJcU+I5qIiIh0g0JbREQkIhTaIiIiEaHQFhER\niQiFtoiISEQotEVERCKiV0PbzIab2a1mtiC5PcjMZpvZzWZ2em/eW0REpK/p1dB291fc/fyUXScD\nv3L3C4ETe/PeIiIifU2XQtvMZprZejNbmra/2syWm9nLZnZxFy6VugrY1ixrFRERKWldbWnPAqpT\nd5hZGXB9cv8oYLKZjTSzs8zsWjPbL8N1VtG2oIjep4uISLc1NDRQVVVFPB6nqqqKhoaGQpfU67q0\nnra7P5FchjPVUcAKd28CMLP5wER3vwq4LblvMPBj4PBkS7weuN7MatBCIiIi0k0NDQ1MmTKFxsbG\n7fta/1xTU1Oosnpdl0K7E6mPuiG0oseknuDubwJfS/vceV25eDweJxaLEYvFiMfjxOPxHpQqIiJ9\nyfTp09sFNoTQrq+vj0RoJxIJEokETU1NWX2uJ6Hdqyt6JBKJ3ry8iIhEWEtLS8b9zc3Nea6ke1Ib\no7Nnz+7y53ryXnk1be+nSf55VQ+uJyIi0iXl5eUZ91dUVOS5kvzqSWgvAQ4ys5iZDQAmoffUIiKS\nB3V1dVRWVrbbV1lZSW1tbYEqyo8uPR43s3nAOGCIma0ELnf3WWZ2EfAQUAbMcPdlvVeqiIhI0Pre\nur6+nubmZioqKqitrY3E++yeMPdefTWdNTNzgGKrS0REpDeYGQDubjs7V2OlRUREIkKhLSIiEhEK\nbRERkYhQaIuIiESEQltERCQiFNoiIiIR0auhbWbDzexWM1uQaVtERES6rldD291fcffzO9sWERGR\nrutSaJvZTDNbb2ZL0/ZXm9lyM3s5ufSmiIiI9JKutrRnAdWpO8ysDLg+uX8UMNnMRprZWWZ2rZnt\nl9tSRURESluXQtvdnwDeStt9FLDC3ZvcfTMwH5jo7re5+7fdfY2ZDTazG4HDzOzi9O2cfhMREZE+\nrifraQ8FVqZsrwLGpJ7g7m8CX0v7XPq2iIiIdEFPQrtXV/SIx+PEYjFisVi7xcJFRESiLpFIkEgk\naGpqyupzXV7ly8xiwP3ufmhy+2hgqrtXJ7cvBba5+0+yqqDjfbTKl4iIlIx8rfK1BDjIzGJmNgCY\nBNzXg+uJiIjIDnR1yNc84ElghJmtNLNz3X0LcBHwEPAicKe7L+u9UkVEREpblx+P54sej4uISCnJ\n1+NxERERySOFtoiISEQotEVERCJCoS0iIhIRCm0REZGIUGiLiIgUSENDQ1bn93pom9lwM7vVzBYk\ntyea2c1mNt/Mju/t+4uIiBSjhoYGpkyZktVn8jZO28wWuPtpKdt7Aj919/PTztM4bRER6fOqqqpY\ntGjR9u2cjtM2s5lmtt7Mlqbtrzaz5Wb2cpbLbf6QsB53n5VIJApdQo/1he8AfeN79IXvAPoexaQv\nfAeI7vdoaWnJ+jPZPB6fBVSn7jCzMkLwVgOjgMlmNtLMzjKza81sv/SLWPAT4Lfu/lzWFUdIVP9D\nStUXvgP0je/RF74D6HsUk77wHSC636O8vDzrz3R5aU53fyK50leqo4AV7t4EYGbzgYnufhVwW3Lf\nYODHwGFmdgnwLnAcsLuZHejuN2VdtYiISMTV1dXR2NhIY2Njlz+T1TvtDMtzngpUufsFye0zgTHu\nXptF3en30MtsEREpOfmYe1wBKyIikiddfjzeidXAsJTtYcCqnlywK79piIiIlKKetrSXAAeZWczM\nBgCTgPt6XpaIiIiky2bI1zzgSWCEma00s3PdfQtwEfAQ8CJwp7sv651SRURESlveJlfpCjOrBq4D\nyoBb3f0nBS4pa2Y2E6gBNrR22IsiMxsGzAH2IfRduNndpxe2quyYWQXwOFAODADudfdLC1tV9yWH\nWC4BVrn7CYWupzvMrAl4B9gKbHb3owpbUfaSE0PdCnyc8HfjPHd/qrBVZcfMDgbmp+z6KPB/I/h3\n/FLgTGAbsBQ4192zH/xcYGY2BTgfMOAWd5/W6bnFEtrJf5BeAsYT3pUvBiZHreVuZmOBjcCciIf2\nh4EPu/tzZrYb8CxwUgT//xjo7pvMrD/we+B77v77QtfVHWb2HeBTwAfc/cRC19MdZvYK8Cl3f7PQ\ntXSXmc0GHnf3mcn/rga5+z8LXVd3mVk/wr+5R7n7ykLX01XJ0UyPASPdvcXM7gQedPfZBS0sS2Z2\nCDAPOBLYDCwEvubuGceBFdOCIdvHfLv7ZsJvgRMLXFPW3P0J4K1C19FT7r6udfIbd98ILAM6TJZT\n7Nx9U/KPAwhPcCIZFma2P/AFQgsv6p01I1u/me0BjHX3mQDuviXKgZ00HmiMUmAnvUMIuYHJX54G\nEn75iJqPAU+7e7O7byU8HTy5s5OLKbSHAqn/0axK7pMCS/5GezjwdGEryZ6Z9TOz54D1wP+4+4uF\nrqmbrgW+T3gMGGUOPGJmS8zsgkIX0w3DgdfMbJaZ/cnMbjGzgYUuqoe+DNxR6CKylXxa8zPgH8Aa\n4G13f6SwVXXLC8BYMxuc/G+pBti/s5OLKbSL4zm9tJN8NP5rYEqyxR0p7r7N3Q8j/CU41sziBS4p\na2b2RUIfiT8T4VZq0mfc/XBgAvDN5OukKOkPjAZucPfRhBkeLylsSd2XHPVzArCg0LVky8wqgW8B\nMcJTwN3M7IyCFtUN7r4c+AmwCPgt8Gd28Mt5MYV2zsd8S8+Y2S7AXcDt7n5PoevpieQjzAbgiELX\n0g2fBk5Mvg+eB3zezOYUuKZucfe1yf99DfgN4bVYlKwidARcnNz+NSHEo2oC8Gzy/4+oOQJ40t3f\nSI5kupvwdyVy3H2mux/h7uOAtwn9uzIqptDWmO8iYmYGzABedPfrCl1Pd5jZ3smevpjZrsDxhN9i\nI8XdL3P3Ye4+nPAo8zF3/0qh68qWmQ00sw8k/zwI+D+EHr+R4e7rgJVmNiK5azzw1wKW1FOTCb8I\nRtFy4Ggz2zX579V4wtDjyDGzfZL/ewDwb+zgdUVPZ0TLGXffYmatY77LgBlR66kM28ezjwOGmNlK\n4HJ3n1XgsrrjM4ShFH8xs9agu9TdFxawpmztC8xO9o7tB9zm7o8WuKZciOqrpA8Bvwn/vtIfmOvu\ni3b8kaJUC8xNNi4agXMLXE+3JH9xGg9EsW8B7v588onTEsLj5D8BNxe2qm77tZkNIXSs+4a7v9PZ\niUUz5EtERER2rJgej4uIiMgOKLRFREQiQqEtIiISEQptERGRiFBoi4iIRIRCW0REJCIU2iIiIhGh\n0BYREYmI/wU/aGBhEoofTwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10b4b8c90>"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Householder transformations are even better"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q3, R3 = np.linalg.qr(H)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.max(np.abs(np.dot(Q3,R3)-H))\n",
      "print np.max(np.abs(np.dot(Q3.T,Q3)-np.eye(n)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.4408920985e-16\n",
        "4.4408920985e-16\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}