{ "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": [ "import flame\n", "import laff\n", "\n", "def LU_unb_var5(A):\n", "\n", " ATL, ATR, \\\n", " ABL, ABR = flame.part_2x2(A, \\\n", " 0, 0, 'TL')\n", "\n", " while ATL.shape[0] < A.shape[0]:\n", "\n", " A00, a01, A02, \\\n", " a10t, alpha11, a12t, \\\n", " A20, a21, A22 = flame.repart_2x2_to_3x3(ATL, ATR, \\\n", " ABL, ABR, \\\n", " 1, 1, 'BR')\n", "\n", " #------------------------------------------------------------#\n", "\n", " laff.invscal( alpha11, a21 ) # a21 := a21 / alpha11\n", " laff.ger( -1.0, a21, a12t, A22 ) # A22 := A22 - a21 * a12t\n", "\n", " #------------------------------------------------------------#\n", "\n", " ATL, ATR, \\\n", " ABL, ABR = flame.cont_with_3x3_to_2x2(A00, a01, A02, \\\n", " a10t, alpha11, a12t, \\\n", " A20, a21, A22, \\\n", " 'TL')\n", "\n", " flame.merge_2x2(ATL, ATR, \\\n", " ABL, ABR, A)" ], "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": [ "laff.trsm('Lower triangular', 'No transpose', 'Unit diagonal', C, V)\n", "laff.trsm('Upper triangular', 'No transpose', 'Nonunit diagonal', C, V)" ], "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", " m,n = B.shape # How many rows and columns does B have?\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", " # C = A^T A\n", " C = np.transpose( A ) * A \n", " # W = A^T B\n", " W = np.transpose( A ) * B\n", " # Overwrite C with its LU factorization\n", " LU_unb_var5( C )\n", " \n", " # Solve L(UX) = W, overwriting W with X\n", " laff.trsm('Lower triangular', 'No transpose', 'Unit diagonal', C, W)\n", " laff.trsm('Upper triangular', 'No transpose', 'Nonunit diagonal', C, W)\n", " \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": {} } ] }