11.3.7 Implementing the QR Factorization - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "11.3.7 Implementing the QR Factorization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With this notebook, you will implement the QR factorization of a matrix $ A $ with linearly independent columns, producing matrices $ Q $ and $ R $ such that $ A = Q R $.  The algorithm is equivalent to performing Gram-Schmidt orthogonalization on the vectors that are the columns of $ A $.\n",
      "\n",
      "<font color=red> Be sure to make a copy!!!! </font>\n",
      "\n",
      "<h2> First, let's create matrices $ A $, $ Q $, and $ R $.  </h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "A = np.matrix( ' 1., -1.,  2;\\\n",
      "                 2.,  1., -3;\\\n",
      "                -1.,  3.,  2;\\\n",
      "                 0., -2., -1' )\n",
      "\n",
      "Q = np.matrix( np.zeros( (4,3) ) )\n",
      "\n",
      "R = np.matrix( np.zeros( (3,3) ) )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Implement the QR factorization algorithm from 11.3.7 </h2>\n",
      "\n",
      "Here is the algorithm:\n",
      "\n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/11_3_7_QR_factorization.png\" alt=\"QR factorization algorithm\" width=75%>\n",
      "    \n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> QR_Gram_Schmidt_unb</code>\n",
      "with the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> for the algorithm\n",
      "\n",
      "Hints:\n",
      "<ul>\n",
      "<li>\n",
      "You will want to store $ a_1^\\perp $ in <code> q1 </code>.  This will mean first copying <code> a1 </code> to <code> q1 </code>.\n",
      "</li>\n",
      "<li>\n",
      "The following routines will be useful:\n",
      "<ul>\n",
      "<li>\n",
      "<code> laff.copy ( x, y ) </code>\n",
      "</li>\n",
      "<li>\n",
      "<code> laff.gemv ( trans, alpha, A, x, beta, y ) </code>\n",
      "</li>\n",
      "<li>\n",
      "<code> laff.norm2 ( x ) </code> <br>\n",
      "The annoying thing is that this returns the value.  You will want do to <br>\n",
      "<code> rho11[:,:] = laff.norm2 </code>.\n",
      "</li>\n",
      "<li>\n",
      "<code> laff.invscal ( alpha, x ) </code>\n",
      "</li>\n",
      "</ul>\n",
      "</li>\n",
      "</ul>\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff as laff\n",
      "\n",
      "def QR_Gram_Schmidt_unb(A, Q, R):\n",
      "\n",
      "    AL, AR = flame.part_1x2(A, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    QL, QR = flame.part_1x2(Q, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    RTL, RTR, \\\n",
      "    RBL, RBR  = flame.part_2x2(R, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    while AL.shape[1] < A.shape[1]:\n",
      "\n",
      "        A0, a1, A2 = flame.repart_1x2_to_1x3(AL, AR, \\\n",
      "                                             1, 'RIGHT')\n",
      "\n",
      "        Q0, q1, Q2 = flame.repart_1x2_to_1x3(QL, QR, \\\n",
      "                                             1, 'RIGHT')\n",
      "\n",
      "        R00,  r01,   R02,  \\\n",
      "        r10t, rho11, r12t, \\\n",
      "        R20,  r21,   R22   = flame.repart_2x2_to_3x3(RTL, RTR, \\\n",
      "                                                     RBL, RBR, \\\n",
      "                                                     1, 1, 'BR')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.gemv( 'Transpose', 1.0, Q0, a1, 0.0, r01 )\n",
      "        laff.copy( a1, q1 )\n",
      "        laff.gemv( 'No transpose', -1.0, Q0, r01, 1.0, q1 )\n",
      "        rho11[:,:] = laff.norm2( q1 )\n",
      "        laff.invscal( rho11, q1 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        AL, AR = flame.cont_with_1x3_to_1x2(A0, a1, A2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "        QL, QR = flame.cont_with_1x3_to_1x2(Q0, q1, Q2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "        RTL, RTR, \\\n",
      "        RBL, RBR  = flame.cont_with_3x3_to_2x2(R00,  r01,   R02,  \\\n",
      "                                               r10t, rho11, r12t, \\\n",
      "                                               R20,  r21,   R22,  \\\n",
      "                                               'TL')\n",
      "\n",
      "    flame.merge_1x2(QL, QR, Q)\n",
      "\n",
      "    flame.merge_2x2(RTL, RTR, \\\n",
      "                    RBL, RBR, R)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> Test the routine </h3>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "QR_Gram_Schmidt_unb( A, Q, R )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n",
      "\n",
      "print( 'Q = ' )\n",
      "print( Q )\n",
      "\n",
      "print( 'R = ' )\n",
      "print( R )\n",
      "\n",
      "print( 'Q * R - A:' )\n",
      "print( Q * R - A )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The result should be a 4 x 3 (approximately) zero matrix.\n",
      "\n",
      "Now, let's check if the columns of $ Q $ are mutually orthogonal:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( np.transpose( Q ) * Q )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above should approximately be a 3 x 3 identity matrix."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Related to the enrichment."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you read the enrichment on the Gram-Schmidt method, you will find the following interesting..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "epsilon = 1.0e-7\n",
      "\n",
      "A = np.matrix( ' 1.,       1.,      1.;\\\n",
      "                 0,        0.,      0.;\\\n",
      "                 0.,       0,       0.;\\\n",
      "                 0.,       0.,      0.' )\n",
      "\n",
      "A[ 1,0 ] = epsilon\n",
      "A[ 2,1 ] = epsilon\n",
      "A[ 3,2 ] = epsilon\n",
      "\n",
      "#  This creates the matrix\n",
      "# A = np.matrix( ' 1.,       1.,      1.;\\\n",
      "#                  epsilon,  0.,      0.;\\\n",
      "#                  0.,       epsilon, 0.;\\\n",
      "#                  0.,       0.,      epsilon' )\n",
      "\n",
      "Q = np.matrix( np.zeros( (4,3) ) )\n",
      "\n",
      "R = np.matrix( np.zeros( (3,3) ) )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "QR_Gram_Schmidt_unb( A, Q, R )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n",
      "\n",
      "print( 'Q = ' )\n",
      "print( Q )\n",
      "\n",
      "print( 'R = ' )\n",
      "print( R )\n",
      "\n",
      "print( 'Q * R - A:' )\n",
      "print( Q * R - A )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This should equal, approximately, the zero matrix.\n",
      "\n",
      "Now, let's check if the columns of Q are mutually orthogonal:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( np.transpose( Q ) * Q )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You will notice that the resulting 3 x 3 matrix, which should (approximately) equal the identity matrix, has off-diagonal entries that do not equal zero at all.  \n",
      "\n",
      "One of them even equals 0.5 (when I tried)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}