{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import numpy as np\n",
      "import numpy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 3\n",
      "\n",
      "e1 = np.array([1,0,0])\n",
      "e2 = np.array([0,1,0])\n",
      "e3 = np.array([0,0,1])\n",
      "\n",
      "A = np.random.randn(3, 3)\n",
      "A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "array([[ 0.47748905, -0.39916192,  0.40139911],\n",
        "       [ 0.43535359,  0.24187108, -0.66839074],\n",
        "       [ 0.30276693, -1.2088603 , -0.87829135]])"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Householder reflector:\n",
      "$$I-2\\frac{vv^T}{v^Tv}$$\n",
      "\n",
      "Choose $v=a-\\|a\\|e_1$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = A[:, 0]\n",
      "v = a-la.norm(a)*e1\n",
      "\n",
      "H1 = ...\n",
      "A1 = np.dot(H1, A)\n",
      "A1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "array([[  7.13579951e-01,  -6.32443387e-01,  -5.11842021e-01],\n",
        "       [  5.55111512e-17,   6.72044039e-01,   1.01563347e+00],\n",
        "       [  8.32667268e-17,  -9.09696239e-01,   2.92864361e-01]])"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "H1 = np.eye(3) - 2*np.outer(v, v)/np.dot(v, v)\n",
      "-->\n",
      "(Edit this cell for solution.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "NB: Never build full Householder matrices in actual code! (Why? How?)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = A1[:, 1].copy()\n",
      "a[0] = 0\n",
      "v = a-la.norm(a)*e2\n",
      "\n",
      "H2 = ...\n",
      "R = np.dot(H2, A1)\n",
      "R"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "array([[  7.13579951e-01,  -6.32443387e-01,  -5.11842021e-01],\n",
        "       [ -3.39885479e-17,   1.13101301e+00,   3.67929286e-01],\n",
        "       [ -9.41255243e-17,  -1.11022302e-16,  -9.90913173e-01]])"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "H2 = np.eye(3) - 2*np.outer(v, v)/np.dot(v, v)\n",
      "-->\n",
      "(Edit this cell for solution.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q = np.dot(H2, H1).T\n",
      "la.norm(np.dot(Q, R) - A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "5.2735593669694936e-16"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}