04.3.1 Matrix vector multiply with transpose matrix - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Matrix-vector multiplication multiply with transpose matrix"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook walks you through how to implement $ y := A^T x + y $ <i> without explicitly transposing the matrix</i>."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Getting started"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will use some functions that are part of our laff library (of which this function will become a part) as well as some routines from the FLAME API (Application Programming Interface) that allows us to write code that closely resembles how we typeset algorithms using the FLAME notation.  These functions are imported with the \"import laff\" and \"import flame\" statements."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Algorithm (via dot products)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/431Mvmult_n_unb_var1.png\" alt=\"Matrix-vector multiplication via dot products algorithm\" width=\"48%\"><image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/431Mvmult_t_unb_var1.png\" alt=\"Matrix-vector multiplication via dot products algorithm\" width=\"48%\">\n",
      "\n",
      "\n",
      "Above, you find two algorithms.  The one on the left computes $ y := A x + y $ via dot products and the one on the right computes $ y := A^T x + y $ via dot products.  You are to implement the one on the right."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Mvmult_t_unb_var1( A, x, y ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given $ A \\in \\mathbb{R}^{m \\times m} $, $ x \\in \\mathbb{R}^m $, and $ y \\in \\mathbb{R}^m $, computes $ y := A^T x + y $.  The \"_t_\" in the name of the routine indicates this is the \"transpose\" matrix-vector multiplication.  \n",
      "\n",
      "The specific laff functions we will use are \n",
      "<ul>\n",
      "<li> <code> laff.dots( x, y, alpha ) </code> which computes $ \\alpha := x^T y + \\alpha $.  </li>\n",
      "</ul>\n",
      "\n",
      "Use the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> to generate a code skeleton.  (Make sure you adjust the name of the routine.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff as laff\n",
      "\n",
      "def Mvmult_t_unb_var1(A, x, y):\n",
      "\n",
      "    AL, AR = flame.part_1x2(A, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    yT, \\\n",
      "    yB  = flame.part_2x1(y, \\\n",
      "                         0, 'TOP')\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",
      "        y0,   \\\n",
      "        psi1, \\\n",
      "        y2    = flame.repart_2x1_to_3x1(yT, \\\n",
      "                                        yB, \\\n",
      "                                        1, 'BOTTOM')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.dots( a1, x, psi1 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        AL, AR = flame.cont_with_1x3_to_1x2(A0, a1, A2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "        yT, \\\n",
      "        yB  = flame.cont_with_3x1_to_2x1(y0,   \\\n",
      "                                         psi1, \\\n",
      "                                         y2,   \\\n",
      "                                         'TOP')\n",
      "\n",
      "    flame.merge_2x1(yT, \\\n",
      "                    yB, y)\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Testing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's quickly test the routine by creating a 3 x 4 matrix and related vectors, performing the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import random\n",
      "from numpy import matrix\n",
      "\n",
      "A = matrix( random.rand( 3,4 ) )\n",
      "x = matrix( random.rand( 3,1 ) )\n",
      "y = matrix( random.rand( 4,1 ) )\n",
      "yold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "print( 'A before =' )\n",
      "print( A )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )\n",
      "\n",
      "print( 'y before =' )\n",
      "print( y )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "A before =\n",
        "[[ 0.45946038  0.36711289  0.62663084  0.69366725]\n",
        " [ 0.61839962  0.32329523  0.54882629  0.52064797]\n",
        " [ 0.92298778  0.09047515  0.40856054  0.23082752]]\n",
        "x before =\n",
        "[[ 0.47016   ]\n",
        " [ 0.77363531]\n",
        " [ 0.30106805]]\n",
        "y before =\n",
        "[[ 0.86834438]\n",
        " [ 0.30196564]\n",
        " [ 0.06375715]\n",
        " [ 0.70927731]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import transpose\n",
      "\n",
      "laff.copy( y, yold )   # save the original vector y\n",
      "\n",
      "Mvmult_t_unb_var1( A, x, y )\n",
      "\n",
      "print( 'y after =' )\n",
      "print( y )\n",
      "\n",
      "print( 'y - ( transpose( A ) * x + yold ) = ' )\n",
      "print( y - ( transpose( A ) * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "y after =\n",
        "[[ 1.84066219]\n",
        " [ 0.75191922]\n",
        " [ 0.90596983]\n",
        " [ 1.50769835]]\n",
        "y - ( transpose( A ) * x + yold ) = \n",
        "[[  0.00000000e+00]\n",
        " [ -1.11022302e-16]\n",
        " [  0.00000000e+00]\n",
        " [  0.00000000e+00]]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bingo, it seems to work!  (Notice that we are doing floating point computations, which means that due to rounding you may not get an exact \"0\".)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Watch your code in action!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Copy and paste the code into <a href=\"http://edx-org-utaustinx.s3.amazonaws.com/UT501x/PictureFlame/PictureFLAME.html\"> PictureFLAME </a>, a webpage where you can watch your routine in action.  Just cut and paste into the box.  \n",
      "\n",
      "Disclaimer: we implemented a VERY simple interpreter.  If you do something wrong, we cannot guarantee the results.  But if you do it right, you are in for a treat.\n",
      "\n",
      "If you want to reset the problem, just click in the box into which you pasted the code and hit \"next\" again."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Algorithm (via <code>axpy</code>s)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/431Mvmult_n_unb_var2.png\" alt=\"Matrix-vector multiplication via axpys algorithm\" width=\"55%\"><image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/431Mvmult_t_unb_var2.png\" alt=\"Matrix-vector multiplication via dot axpys algorithm\" width=\"45%\">\n",
      "\n",
      "Above, you find two algorithms.  The one on the left computes $ y := A x + y $ via <code> axpy</code>s and the one on the right computes $ y := A^T x + y $ via <code> axpy</code>s.  You are to implement the one on the right."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Mvmult_t_unb_var2( A, x, y ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given $ A \\in \\mathbb{R}^{m \\times m} $, $ x \\in \\mathbb{R}^m $, and $ y \\in \\mathbb{R}^m $, computes $ y := A^T x + y $.  The \"_t_\" in the name of the routine indicates this is the \"transpose\" matrix-vector multiplication.  \n",
      "\n",
      "The specific laff functions we will use are \n",
      "<ul>\n",
      "<li> <code> laff.axpy( alpha, x, y ) </code> which computes $ y := \\alpha x + y $.  </li>\n",
      "</ul>\n",
      "\n",
      "Use the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> to generate a code skeleton.  (Make sure you adjust the name of the routine.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff as laff\n",
      "\n",
      "def Mvmult_t_unb_var1(A, x, y):\n",
      "\n",
      "    AT, \\\n",
      "    AB  = flame.part_2x1(A, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    xT, \\\n",
      "    xB  = flame.part_2x1(x, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    while AT.shape[0] < A.shape[0]:\n",
      "\n",
      "        A0,  \\\n",
      "        a1t, \\\n",
      "        A2   = flame.repart_2x1_to_3x1(AT, \\\n",
      "                                       AB, \\\n",
      "                                       1, 'BOTTOM')\n",
      "\n",
      "        x0,   \\\n",
      "        chi1, \\\n",
      "        x2    = flame.repart_2x1_to_3x1(xT, \\\n",
      "                                        xB, \\\n",
      "                                        1, 'BOTTOM')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.axpy( chi1, a1t, y )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        AT, \\\n",
      "        AB  = flame.cont_with_3x1_to_2x1(A0,  \\\n",
      "                                         a1t, \\\n",
      "                                         A2,  \\\n",
      "                                         'TOP')\n",
      "\n",
      "        xT, \\\n",
      "        xB  = flame.cont_with_3x1_to_2x1(x0,   \\\n",
      "                                         chi1, \\\n",
      "                                         x2,   \\\n",
      "                                         'TOP')\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Testing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's quickly test the routine by creating a 3 x 4 matrix and related vectors, performing the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import random\n",
      "from numpy import matrix\n",
      "\n",
      "A = matrix( random.rand( 3,4 ) )\n",
      "x = matrix( random.rand( 3,1 ) )\n",
      "y = matrix( random.rand( 4,1 ) )\n",
      "yold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "print( 'A before =' )\n",
      "print( y )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )\n",
      "\n",
      "print( 'y before =' )\n",
      "print( y )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "A before =\n",
        "[[ 0.59304972]\n",
        " [ 0.18665298]\n",
        " [ 0.09586729]\n",
        " [ 0.71568854]]\n",
        "x before =\n",
        "[[ 0.2797929 ]\n",
        " [ 0.70916909]\n",
        " [ 0.73616771]]\n",
        "y before =\n",
        "[[ 0.59304972]\n",
        " [ 0.18665298]\n",
        " [ 0.09586729]\n",
        " [ 0.71568854]]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import transpose\n",
      "\n",
      "laff.copy( y, yold )   # save the original vector y\n",
      "\n",
      "Mvmult_t_unb_var2( A, x, y )\n",
      "\n",
      "print( 'y after =' )\n",
      "print( y )\n",
      "\n",
      "print( 'y - ( transpose( A ) * x + yold ) = ' )\n",
      "print( y - ( transpose( A ) * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "y after =\n",
        "[[ 1.03703548]\n",
        " [ 1.19297584]\n",
        " [ 0.95849002]\n",
        " [ 1.58112475]]\n",
        "y - ( transpose( A ) * x + yold ) = \n",
        "[[ 0.]\n",
        " [ 0.]\n",
        " [ 0.]\n",
        " [ 0.]]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bingo, it seems to work!  (Notice that we are doing floating point computations, which means that due to rounding you may not get an exact \"0\".)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Watch your code in action!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Copy and paste the code into <a href=\"http://edx-org-utaustinx.s3.amazonaws.com/UT501x/PictureFlame/PictureFLAME.html\"> PictureFLAME </a>, a webpage where you can watch your routine in action.  Just cut and paste into the box.  \n",
      "\n",
      "Disclaimer: we implemented a VERY simple interpreter.  If you do something wrong, we cannot guarantee the results.  But if you do it right, you are in for a treat.\n",
      "\n",
      "If you want to reset the problem, just click in the box into which you pasted the code and hit \"next\" again."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}