{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 3x3 Householder QR Demo"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This demo constructs a $3\\times 3$ QR factorization using Householder reflectors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.06427586,  1.03785412,  1.19643376],\n",
              "       [-1.06405176,  0.48693156, -0.83335032],\n",
              "       [-0.20227028, -1.10857722,  0.27986689]])"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "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(n, n)\n",
        "A"
      ]
    },
    {
      "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",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = A[:, 0]\n",
        "v = a-la.norm(a)*e1\n",
        "\n",
        "H1 = np.eye(3) - 2*np.outer(v, v)/(v@v)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.51848691e+00,   5.33870203e-01,   1.38523070e+00],\n",
              "       [  3.40401588e-16,  -6.93719959e-01,  -3.91067578e-01],\n",
              "       [  9.99443798e-17,  -1.33301245e+00,   3.63942360e-01]])"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "A1 = H1 @ A\n",
        "A1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "NB: Never build full Householder matrices in actual code! (Why? How?)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = A1[:, 1].copy()\n",
        "a[0] = 0\n",
        "v = a-la.norm(a)*e2\n",
        "\n",
        "H2 = np.eye(3) - 2*np.outer(v, v)/(v@v)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.51848691e+00,   5.33870203e-01,   1.38523070e+00],\n",
              "       [ -2.45801147e-16,   1.50272073e+00,  -1.42307423e-01],\n",
              "       [ -2.55820086e-16,   1.61523465e-16,   5.14914060e-01]])"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "R = H2 @ A1\n",
        "R"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "collapsed": false,
        "scrolled": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "8.5458331069102901e-16"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Q = np.dot(H2, H1).T\n",
        "la.norm(np.dot(Q, R) - A)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}