{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Keeping track of coefficients in Gram-Schmidt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = np.random.randn(3, 3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start from regular old (modified) Gram-Schmidt:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "\n",
        "Q = np.zeros(A.shape)\n",
        "\n",
        "q = A[:, 0]\n",
        "Q[:, 0] = q/la.norm(q)\n",
        "\n",
        "# -----------\n",
        "\n",
        "q = A[:, 1]\n",
        "coeff = np.dot(Q[:, 0], q)\n",
        "q = q - coeff*Q[:, 0]\n",
        "Q[:, 1] = q/la.norm(q)\n",
        "\n",
        "# -----------\n",
        "\n",
        "q = A[:, 2]\n",
        "coeff = np.dot(Q[:, 0], q)\n",
        "q = q - coeff*Q[:, 0]\n",
        "coeff = np.dot(Q[:, 1], q)\n",
        "q = q - coeff*Q[:, 1]\n",
        "Q[:, 2] = q/la.norm(q)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.00000000e+00,   6.15868752e-17,  -5.86239841e-16],\n",
              "       [  6.15868752e-17,   1.00000000e+00,  -3.18779032e-16],\n",
              "       [ -5.86239841e-16,  -3.18779032e-16,   1.00000000e+00]])"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Q.dot(Q.T)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we want to keep track of what vector got added to what other vector, in the style of an elimination matrix.\n",
        "\n",
        "Let's call that matrix $R$.\n",
        "\n",
        "* Would it be $A=QR$ or $A=RQ$? Why?\n",
        "* Where are $R$'s nonzeros?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "R = np.zeros((A.shape[0], A.shape[0]))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Q = np.zeros(A.shape)\n",
        "\n",
        "q = A[:, 0]\n",
        "Q[:, 0] = q/la.norm(q)\n",
        "\n",
        "R[0,0] = la.norm(q)\n",
        "\n",
        "# -----------\n",
        "\n",
        "q = A[:, 1]\n",
        "coeff = np.dot(Q[:, 0], q)\n",
        "R[0,1] = coeff\n",
        "q = q - coeff*Q[:, 0]\n",
        "Q[:, 1] = q/la.norm(q)\n",
        "\n",
        "R[1,1] = la.norm(q)\n",
        "\n",
        "# -----------\n",
        "\n",
        "q = A[:, 2]\n",
        "coeff = np.dot(Q[:, 0], q)\n",
        "R[0,2] = coeff\n",
        "q = q - coeff*Q[:, 0]\n",
        "coeff = np.dot(Q[:, 1], q)\n",
        "R[1,2] = coeff\n",
        "q = q- coeff*Q[:, 1]\n",
        "Q[:, 2] = q/la.norm(q)\n",
        "\n",
        "R[2,2] = la.norm(q)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.37598334,  1.41035995, -2.41234028],\n",
              "       [ 0.        ,  0.79661434,  1.04638607],\n",
              "       [ 0.        ,  0.        ,  0.54718387]])"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "R"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "5.5511151231257827e-17"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(Q@R - A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is called [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----------\n",
        "* When does it break?\n",
        "* Does it need something like pivoting?\n",
        "* Can we use it for something?"
      ]
    }
  ],
  "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
}