{"nbformat": 3, "nbformat_minor": 0, "metadata": {"name": "", "signature": "sha256:f3474079a26bfc08e0d1a45b1d4b17bb9eba2d3f592496c73cd9f102814618cc"}, "worksheets": [{"metadata": {}, "cells": [{"metadata": {}, "cell_type": "markdown", "source": ["# Gram-Schmidt and Modified Gram-Schmidt"]}, {"language": "python", "outputs": [], "input": ["import numpy as np\n", "import numpy.linalg as la"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 1}, {"language": "python", "outputs": [], "input": ["A = np.random.randn(3, 3)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 35}, {"language": "python", "outputs": [], "input": ["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", "collapsed": false, "metadata": {}, "prompt_number": 36}, {"language": "python", "outputs": [], "input": ["Q = np.zeros(A.shape)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 37}, {"metadata": {}, "cell_type": "markdown", "source": ["Now let us generalize the process we used for three vectors earlier:"]}, {"language": "python", "outputs": [], "input": ["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": "code", "collapsed": false, "metadata": {}, "prompt_number": 38}, {"metadata": {}, "cell_type": "markdown", "source": ["This procedure is called [Gram-Schmidt Orthonormalization](https://en.wikipedia.org/wiki/Gram\u2013Schmidt_process)."]}, {"language": "python", "outputs": [{"stream": "stdout", "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"], "output_type": "stream"}], "input": ["test_orthogonality(Q)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 39}, {"metadata": {}, "cell_type": "markdown", "source": ["Now let us try a different example ([Source](http://fgiesen.wordpress.com/2013/06/02/modified-gram-schmidt-orthogonalization/)):"]}, {"language": "python", "outputs": [{"prompt_number": 2, "output_type": "pyout", "text": ["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]])"], "metadata": {}}], "input": ["\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", "collapsed": false, "metadata": {}, "prompt_number": 2}, {"language": "python", "outputs": [], "input": ["Q = np.zeros(A.shape)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 3}, {"language": "python", "outputs": [{"stream": "stdout", "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"], "output_type": "stream"}], "input": ["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", "collapsed": false, "metadata": {}, "prompt_number": 17}, {"language": "python", "outputs": [{"stream": "stdout", "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"], "output_type": "stream"}], "input": ["test_orthogonality(Q)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 43}, {"metadata": {}, "cell_type": "markdown", "source": ["Questions:\n", "\n", "* What happened?\n", "* How do we fix it?"]}, {"language": "python", "outputs": [], "input": ["Q = np.zeros(A.shape)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 44}, {"language": "python", "outputs": [], "input": ["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", "collapsed": false, "metadata": {}, "prompt_number": 47}, {"language": "python", "outputs": [{"stream": "stdout", "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"], "output_type": "stream"}], "input": ["test_orthogonality(Q)"], "cell_type": "code", "collapsed": false, "metadata": {}, "prompt_number": 48}, {"metadata": {}, "cell_type": "markdown", "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$?"]}]}]}