11.2.5 Rank-k Approximation - Answer.ipynb

{
 "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": [
      "<font color=red> Be sure to make a copy!!!! </font>\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 <code> RankKApprox( B, k ) </code>\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": [
      "<font color=blue> # C = A^T A: </font> This comment corresponds to computing $C = A^T A$. Try using numpy's built in transpose method by calling `np.transpse(A)`. \n",
      "\n",
      "<font color=blue># W = A^T B: </font> This comment corresponds to computing $W = A^T B$. See above for a hint.\n",
      "\n",
      "<font color=blue># Overwrite C with its LU factorization: </font> Use your implementation of LU_unb_var5 from an earlier notebook for this.\n",
      "    \n",
      "<font color=blue># Solve L(UX) = W, overwriting W with X: </font> Use `laff.trsm` to do this. Recall that \"trsm\" means triangular solve with multiple right hand sides."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}