{ "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", " Be sure to make a copy!!!! \n", "\n", "

First, let's create matrices $ A $, $ Q $, and $ R $.

" ] }, { "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": [ "

Implement the QR factorization algorithm from 11.3.7

\n", "\n", "Here is the algorithm:\n", "\n", "\"QR\n", " \n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " QR_Gram_Schmidt_unb\n", "with the Spark webpage for the algorithm\n", "\n", "Hints:\n", "\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": [ "

Test the routine

\n", "\n", " 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. " ] }, { "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": {} } ] }