{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"A = \n",
"[[ 1. -1. 2.]\n",
" [ 2. 1. -3.]\n",
" [-1. 3. 2.]\n",
" [ 0. -2. -1.]]\n"
]
}
],
"prompt_number": 14
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Implement the QR factorization algorithm from 11.3.7
\n",
"\n",
"Here is the algorithm:\n",
"\n",
"![\"QR](\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/11_3_7_QR_factorization.png\")
\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",
"- \n",
"You will want to store $ a_1^\\perp $ in
q1 . This will mean first copying a1 to q1 .\n",
" \n",
"- \n",
"The following routines will be useful:\n",
"
\n",
"- \n",
"
laff.copy ( x, y ) \n",
" \n",
"- \n",
"
laff.gemv ( trans, alpha, A, x, beta, y ) \n",
" \n",
"- \n",
"
laff.norm2 ( x )
\n",
"The annoying thing is that this returns the value. You will want do to
\n",
" rho11[:,:] = laff.norm2 .\n",
" \n",
"- \n",
"
laff.invscal ( alpha, x ) \n",
" \n",
"
\n",
" \n",
"
\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Insert code here\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 15
},
{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"A = \n",
"[[ 1. -1. 2.]\n",
" [ 2. 1. -3.]\n",
" [-1. 3. 2.]\n",
" [ 0. -2. -1.]]\n",
"Q = \n",
"[[ 0.40824829 -0.17609018 0.8820199 ]\n",
" [ 0.81649658 0.44022545 -0.32318287]\n",
" [-0.40824829 0.70436073 0.23565417]\n",
" [ 0. -0.52827054 -0.24912013]]\n",
"R = \n",
"[[ 2.44948974 -0.81649658 -2.44948974]\n",
" [ 0. 3.7859389 0.26413527]\n",
" [ 0. 0. 3.45401687]]\n",
"Q * R - A:\n",
"[[ 0.00000000e+00 0.00000000e+00 -2.22044605e-16]\n",
" [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n"
]
}
],
"prompt_number": 16
},
{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[[ 1.00000000e+00 1.11022302e-16 1.66533454e-16]\n",
" [ 1.11022302e-16 1.00000000e+00 2.77555756e-17]\n",
" [ 1.52655666e-16 5.55111512e-17 1.00000000e+00]]\n"
]
}
],
"prompt_number": 17
},
{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"A = \n",
"[[ 1.00000000e+00 1.00000000e+00 1.00000000e+00]\n",
" [ 1.00000000e-07 0.00000000e+00 0.00000000e+00]\n",
" [ 0.00000000e+00 1.00000000e-07 0.00000000e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 1.00000000e-07]]\n"
]
}
],
"prompt_number": 38
},
{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"A = \n",
"[[ 1.00000000e+00 1.00000000e+00 1.00000000e+00]\n",
" [ 1.00000000e-07 0.00000000e+00 0.00000000e+00]\n",
" [ 0.00000000e+00 1.00000000e-07 0.00000000e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 1.00000000e-07]]\n",
"Q = \n",
"[[ 1.00000000e+00 6.90840682e-08 4.07886015e-08]\n",
" [ 1.00000000e-07 -7.07106781e-01 -4.17602698e-01]\n",
" [ 0.00000000e+00 7.07106781e-01 -3.98821881e-01]\n",
" [ 0.00000000e+00 0.00000000e+00 8.16424579e-01]]\n",
"R = \n",
"[[ 1.00000000e+00 1.00000000e+00 1.00000000e+00]\n",
" [ 0.00000000e+00 1.41421356e-07 6.90840682e-08]\n",
" [ 0.00000000e+00 0.00000000e+00 1.22485288e-07]]\n",
"Q * R - A:\n",
"[[ 0. 0. 0.]\n",
" [ 0. 0. 0.]\n",
" [ 0. 0. 0.]\n",
" [ 0. 0. 0.]]\n"
]
}
],
"prompt_number": 39
},
{
"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": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[[ 1.00000000e+00 -1.62660994e-09 -9.71668373e-10]\n",
" [ -1.62660994e-09 1.00000000e+00 1.32800433e-02]\n",
" [ -9.71668373e-10 1.32800433e-02 1.00000000e+00]]\n"
]
}
],
"prompt_number": 40
},
{
"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": {}
}
]
}