{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Gram-Schmidt and Modified Gram-Schmidt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = np.random.randn(3, 3)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def test_orthogonality(Q):\n",
        "    print(\"Q:\")\n",
        "    print(Q)\n",
        "    \n",
        "    print(\"Q^T Q:\")\n",
        "    QtQ = np.dot(Q.T, Q)\n",
        "    QtQ[np.abs(QtQ) < 1e-15] = 0\n",
        "    print(QtQ)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Q = np.zeros(A.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let us generalize the process we used for three vectors earlier:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for k in range(A.shape[1]):\n",
        "    avec = A[:, k]\n",
        "    q = avec\n",
        "    for j in range(k):\n",
        "        q = q - np.dot(avec, Q[:,j])*Q[:,j]\n",
        "    \n",
        "    Q[:, k] = q/la.norm(q)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This procedure is called [Gram-Schmidt Orthonormalization](https://en.wikipedia.org/wiki/Gram\u2013Schmidt_process)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Q:\n",
            "[[-0.6932589320501 -0.6855758663147 -0.2222111263183]\n",
            " [-0.7199408381809  0.6447564063972  0.2568547564853]\n",
            " [ 0.032821374928  -0.3380457187077  0.9405572015626]]\n",
            "Q^T Q:\n",
            "[[ 1.  0.  0.]\n",
            " [ 0.  1.  0.]\n",
            " [ 0.  0.  1.]]\n"
          ]
        }
      ],
      "source": [
        "test_orthogonality(Q)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let us try a different example ([Source](http://fgiesen.wordpress.com/2013/06/02/modified-gram-schmidt-orthogonalization/)):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.0000000000000e+00,   1.0000000000000e+00,   1.0000000000000e+00],\n",
              "       [  1.0000000000000e-08,   1.0000000000000e-08,   0.0000000000000e+00],\n",
              "       [  1.0000000000000e-08,   0.0000000000000e+00,   1.0000000000000e-08]])"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "\n",
        "np.set_printoptions(precision=13)\n",
        "\n",
        "eps = 1e-8\n",
        "\n",
        "A = np.array([\n",
        "    [1,  1,  1],\n",
        "    [eps,eps,0],\n",
        "    [eps,0,  eps]\n",
        "    ])\n",
        "\n",
        "A"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Q = np.zeros(A.shape)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[  1.0000000000000e+00   1.0000000000000e-08   1.0000000000000e-08]\n",
            "norm --> [  1.0000000000000e+00   1.0000000000000e-08   1.0000000000000e-08]\n",
            "-------\n",
            "[  1.0000000000000e+00   1.0000000000000e-08   0.0000000000000e+00]\n",
            "[  0.0000000000000e+00   0.0000000000000e+00  -1.0000000000000e-08]\n",
            "norm --> [ 0.  0. -1.]\n",
            "-------\n",
            "[  1.0000000000000e+00   0.0000000000000e+00   1.0000000000000e-08]\n",
            "[  0.0000000000000e+00  -1.0000000000000e-08   0.0000000000000e+00]\n",
            "[  0.0000000000000e+00  -1.0000000000000e-08  -1.0000000000000e-08]\n",
            "norm --> [ 0.              -0.7071067811865 -0.7071067811865]\n",
            "-------\n"
          ]
        }
      ],
      "source": [
        "for k in range(A.shape[1]):\n",
        "    avec = A[:, k]\n",
        "    q = avec\n",
        "    for j in range(k):\n",
        "        print(q)\n",
        "        q = q - np.dot(avec, Q[:,j])*Q[:,j]\n",
        "    \n",
        "    print(q)\n",
        "    q = q/la.norm(q)\n",
        "    Q[:, k] = q\n",
        "    print(\"norm -->\", q)\n",
        "    print(\"-------\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Q:\n",
            "[[  1.0000000000000e+00   0.0000000000000e+00   0.0000000000000e+00]\n",
            " [  1.0000000000000e-08   0.0000000000000e+00  -7.0710678118655e-01]\n",
            " [  1.0000000000000e-08  -1.0000000000000e+00  -7.0710678118655e-01]]\n",
            "Q^T Q:\n",
            "[[  1.0000000000000e+00  -1.0000000000000e-08  -1.4142135623731e-08]\n",
            " [ -1.0000000000000e-08   1.0000000000000e+00   7.0710678118655e-01]\n",
            " [ -1.4142135623731e-08   7.0710678118655e-01   1.0000000000000e+00]]\n"
          ]
        }
      ],
      "source": [
        "test_orthogonality(Q)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Questions:\n",
        "\n",
        "* What happened?\n",
        "* How do we fix it?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 44,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Q = np.zeros(A.shape)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 47,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for k in range(A.shape[1]):\n",
        "    q = A[:, k]\n",
        "    for j in range(k):\n",
        "        q = q - np.dot(q, Q[:,j])*Q[:,j]\n",
        "    \n",
        "    Q[:, k] = q/la.norm(q)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 48,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Q:\n",
            "[[  1.0000000000000e+00   0.0000000000000e+00   0.0000000000000e+00]\n",
            " [  1.0000000000000e-08   0.0000000000000e+00  -1.0000000000000e+00]\n",
            " [  1.0000000000000e-08  -1.0000000000000e+00   0.0000000000000e+00]]\n",
            "Q^T Q:\n",
            "[[  1.0000000000000e+00  -1.0000000000000e-08  -1.0000000000000e-08]\n",
            " [ -1.0000000000000e-08   1.0000000000000e+00   0.0000000000000e+00]\n",
            " [ -1.0000000000000e-08   0.0000000000000e+00   1.0000000000000e+00]]\n"
          ]
        }
      ],
      "source": [
        "test_orthogonality(Q)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This procedure is called *Modified* Gram-Schmidt Orthogonalization.\n",
        "\n",
        "Questions:\n",
        "\n",
        "* Is there a difference mathematically between modified and unmodified?\n",
        "* Why are there $10^{-8}$ values left in $Q^TQ$?"
      ]
    }
  ],
  "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
}