{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import numpy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# http://fgiesen.wordpress.com/2013/06/02/modified-gram-schmidt-orthogonalization/\n",
      "\n",
      "np.set_printoptions(precision=13)\n",
      "\n",
      "if 0:\n",
      "    eps = 1e-8\n",
      "\n",
      "    A = np.array([\n",
      "        [1,  1,  1],\n",
      "        [eps,eps,0],\n",
      "        [eps,0,  eps]\n",
      "        ])\n",
      "else:\n",
      "    A = np.random.randn(3, 3)\n",
      "\n",
      "A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "array([[ 0.5253995478494,  0.0728870856098, -0.1239374131572],\n",
        "       [ 1.4521672418306,  0.2884259189745, -0.9646423363317],\n",
        "       [ 1.9552560268667, -0.4250653013151,  0.4126165423083]])"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def test_orthogonality(orthogonal_vectors):\n",
      "    Q = np.array(orthogonal_vectors)\n",
      "\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"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "orthogonal_vectors = []\n",
      "\n",
      "for k in range(len(A.T)):\n",
      "    q0 = q = A[:, k]\n",
      "    for prev_q in orthogonal_vectors:\n",
      "        q = q - np.dot(prev_q, q0)*prev_q\n",
      "    q = q/la.norm(q)\n",
      "    \n",
      "    orthogonal_vectors.append(q)\n",
      "\n",
      "test_orthogonality(orthogonal_vectors)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Q:\n",
        "[[ 0.2108719041425  0.5828350493858  0.7847524101592]\n",
        " [ 0.2104899290795  0.7568976003929 -0.6187083418506]\n",
        " [ 0.95458212313   -0.2956506853143 -0.0369275300254]]\n",
        "Q^T Q:\n",
        "[[ 1.  0.  0.]\n",
        " [ 0.  1.  0.]\n",
        " [ 0.  0.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "This code:\n",
      "q = q - np.dot(prev_q, q0)*prev_q\n",
      "will fail for the second example above.\n",
      "-->\n",
      "<!--\n",
      "Modified Gram-Schmidt uses:\n",
      "q = q - np.dot(prev_q, q)*prev_q\n",
      "which works.\n",
      "-->\n",
      "(Edit this cell to see instructions.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* Mathematically: fine.\n",
      "* Issue with computed solution?\n",
      "* Pros/cons of `q0`?\n",
      "* Fix?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "my_A = A.copy()\n",
      "orthogonal_vectors = []\n",
      "\n",
      "for k in range(len(A.T)):\n",
      "    q = my_A[:, k]\n",
      "    q = q/la.norm(q)\n",
      "    \n",
      "    for j in range(k+1, len(A.T)):\n",
      "        ...\n",
      "        \n",
      "    orthogonal_vectors.append(q)\n",
      "\n",
      "test_orthogonality(orthogonal_vectors)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Q:\n",
        "[[ 0.2108719041425  0.5828350493858  0.7847524101592]\n",
        " [ 0.2104899290795  0.7568976003929 -0.6187083418506]\n",
        " [ 0.95458212313   -0.2956506853143 -0.0369275300254]]\n",
        "Q^T Q:\n",
        "[[ 1.  0.  0.]\n",
        " [ 0.  1.  0.]\n",
        " [ 0.  0.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "my_A[:, j] -= np.dot(my_A[:, j], q)*q\n",
      "-->\n",
      "(Edit this cell to see solution.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}