{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"11.2.5 Rank-k Approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Be sure to make a copy!!!! \n",
"\n",
"This notebook walks you through the operations required to compute a low rank approximation of a matrix $ B $. We will create a matrix $A$ whose column space will be used in the approximation of $B$.\n",
"\n",
"We start by creating a random $ m \\times n $ matrix $ B $. We then take $ k $ columns of $ B $ to be matrix $ A $, whose columns will be used in the approximation $ B \\approx A V $.\n",
"(In the text and the videos, we talk about computing $ W $ so that $ B \\approx A W^T $. Here we find it more convenient to compute the transpose of that matrix instead. We call it $ V $ to distinguish it from $ W $. So, $ W = V^T $.)\n",
"\n",
"$ V $ is computed as $ ( A^T A )^{-1} A^T B $. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"import laff\n",
"import flame\n",
"\n",
"m = 8\n",
"n = 8\n",
"k = 3\n",
"\n",
"# Random matrix of size mxn\n",
"B = np.matrix( np.random.rand( m, n ) )\n",
"\n",
"# A is k columns of B taken at even intervals\n",
"if 2*k <= n: #k is less than half of n\n",
" interval = np.ceil( n/k ) \n",
" A = B[ :, ::interval ] # These are slices in python. \n",
" # This says take all rows of B, and columns \n",
" # from 0 to the end at interval steps\n",
"else:\n",
" A = B[ :, :k] #If k is greater than half of n, then just take the first k columns\n",
"\n",
"print( 'A = ' )\n",
"print( A )\n",
"\n",
"print( 'B = ' )\n",
"print( B )\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start the process of computing $( A^T A )^{-1} A^T B$ by computing $ A^T A $ and storing the result in a matrix, $C$. In this implementation, we ignore symmetry."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"C = np.transpose( A ) * A \n",
"\n",
"print( 'C = ' )\n",
"print( C )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, form $ V = A^T B $, notice that we are not done forming $V$ after this step."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"V = np.transpose( A ) * B"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instead of computing $ C^{-1} = ( A^T A )^{-1} $ explicitly, we notice that we can instead store the $ L $ and $ U $ factorization of $C$ in $ C $ and then just solve $ L ( U X ) = V $. First we will overwrite $ V $ with the \n",
"result of solving $ L Z = V $, and then we will overwrite $ V $ with the result of solving $ U X = V $."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Copy your `LU_unb_var5` routine from *Notebook 6.3: Solving A x b via LU Factorization and Triangular Solves*"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Insert your LU_unb_var5 method here"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now run `LU_unb_var5` on the matrix $C$ to store $L$ and $U$ in it."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"LU_unb_var5( C )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve $L ( U X ) = V$, overwriting $V$ where $U$ and $L$ are stored in the upper and the strictly lower portions of $C$ respectively."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Insert Code Here...\n",
"# Use laff.trsm(uplo, transpose, diag, A, B) to solve A X = B overwriting B with X"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The $ j $th column of $ A V $ now equals the projection of the $ j $th column of $ B $ onto the column space of $ A $, $ {\\cal C}( A ) $. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A couple of notes:\n",
" \n",
"- The matrix $ A^T A $ is symmetric positive definite. As a result, one does not need to pivot when performing LU factorization. (The reason for this is beyond the scope of this course.)\n",
"\n",
"- One could use what is called a \"Symmetric rank-k update\" operation to compute only the lower (or upper) triangular part of $ A^T A $. This would (approximately) halve the number of floating point operations that are required.\n",
"\n",
"- In one of the enrichments, 8.4.2, we discussed the Cholesky factorization of a symmetric positive definite matrix. One should ideally use that here, since it also takes advantage of symmetry.\n",
"\n",
"- This would then leave us with $ L $, a lower triangular matrix, such that $ C = A^T A = L L^T $. Computing $ V $ would then require the steps\n",
" - $ V = A^T B $.\n",
" - Solve $ L Z = V $ overwriting $ V $ with $ Z $.\n",
" - Solve $ L^T X = V $ overwriting $ V $ with $ X $.\n",
" "
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"A routine"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above computation should be implemented as the routine RankKApprox( B, k )
\n",
"where $ B $ is the $ m \\times n $ matrix to be approximated, and $k$ is the rank of the eventual approximation that will be returned by the method. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def RankKApprox( B, k ):\n",
" \n",
" m,n = B.shape\n",
" \n",
" # A is k columns of B taken at even intervals\n",
" if 2*k <= n: #k is less than half of n\n",
" interval = np.ceil( n/k ) \n",
" A = B[ :, ::interval ] # These are slices in python. \n",
" # This says take all rows of B, and columns \n",
" # from 0 to the end at interval steps\n",
" else:\n",
" A = B[ :, :k] #If k is greater than half of n, then just take the first k columns\n",
"\n",
" #------------------------------------------------------------#\n",
" \n",
" # Replace the comments below with their respective operations from the notebook\n",
" \n",
" # C = A^T A\n",
" # W = A^T B\n",
" # Overwrite C with its LU factorization\n",
" # Solve L(UX) = W, overwriting W with X\n",
" \n",
" #------------------------------------------------------------#\n",
" \n",
" # Return A * W\n",
" return A * W"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"An Application: Rank k Image Approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have implemented routines to create low rank approximations to matrices we will explore what a rank k approximation to an image looks like. Each pixel in an image can be thought of as a value in a matrix. For a grayscale image, this value corresponds to how black or white it is on a relative scale.\n",
"\n",
"We will use two techniques for these approximations. First, the normal approximation developed above and second, the SVD which is a very useful matrix decomposition that guarantees the best approximation given $k$ columns. The SVD might take a while to compute, so don't panic if one of the code blocks takes a bit to complete.\n",
"\n",
"Try experimenting with the number of columns below by changing the `numCols` variable.\n",
"\n",
"If you want to use your own images, make sure that they are in the `png` format and just place them in your notebooks directoy. Then change the `filename` variable to reflect the file name of your image."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
" %pylab inline"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from im_approx import *\n",
"\n",
"# Try varying the number of columns used for the approximation \n",
"numCols = 40\n",
"\n",
"# Load an image\n",
"filename = 'building.png'\n",
"\n",
"img = np.matrix(read_image( filename ))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Make the approximations using normal equations and the SVD\n",
"# This step might take a while as it is a lot of computation.\n",
"normalApprox, SVDApprox = create_approximations( img, k=numCols, approximator=RankKApprox )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Now plot the approximations that we have created\n",
"\n",
"# Note that we're having some issues with our \n",
"# approximations being somewhat darker than the \n",
"# real image and are investigating.\n",
"plot_approximations( img, normalApprox, SVDApprox, k=numCols )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Hints for implementing RankKApprox"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" # C = A^T A: This comment corresponds to computing $C = A^T A$. Try using numpy's built in transpose method by calling `np.transpse(A)`. \n",
"\n",
"# W = A^T B: This comment corresponds to computing $W = A^T B$. See above for a hint.\n",
"\n",
"# Overwrite C with its LU factorization: Use your implementation of LU_unb_var5 from an earlier notebook for this.\n",
" \n",
"# Solve L(UX) = W, overwriting W with X: Use `laff.trsm` to do this. Recall that \"trsm\" means triangular solve with multiple right hand sides."
]
}
],
"metadata": {}
}
]
}