05.5.1 Multiplying upper triangular matrices.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Multiplying upper triangular matrices"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook helps you implement the operation $ C := U R $ where $ C, U, R \\in \\mathbb{R}^{n \\times n} $, and $ U $ and $ R $ are upper triangular.  $ U $ and $ R $ are stored in the upper triangular part of numpy matrices <code> U </code> and <code> R </code>. The upper triangular part of matrix <code> C </code> is to be overwritten with the resulting upper triangular matrix.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, we create some matrices."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "n = 5\n",
      "\n",
      "C = np.matrix( np.random.random( (n, n) ) )\n",
      "print( 'C = ' )\n",
      "print( C )\n",
      "\n",
      "Cold = np.matrix( np.zeros( (n,n ) ) )\n",
      "Cold = np.matrix( np.copy( C ) )           # an alternative way of doing a \"hard\" copy, in this case of a matrix\n",
      "    \n",
      "U = np.matrix( np.random.random( (n, n) ) )\n",
      "print( 'U = ' )\n",
      "print( U )\n",
      "\n",
      "R = np.matrix( np.random.random( (n, n) ) )\n",
      "print( 'R = ' )\n",
      "print( R )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "<h2>The algorithm  </h2>\n",
      "\n",
      "(The homework was to propose the updates to $ C $.)\n",
      "\n",
      "<image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/Trtrmm_uu_unb_var1_blank.png\" alt=\"Upper triangular matrix-matrix multiplication\" width=\"80%\">\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> The routine <code> Trtrmm_uu_unb_var1( U, R, C ) </code> </h2>\n",
      "\n",
      "This routine computes $ C := U R + C $.  The \"\\_uu\\_\" means that $ U $ and $ R $ are upper triangular matrices (which means $ C $ is too).  However, the lower triangular parts of numpy matrices <code> U </code>, <code> R </code>, and <code> C </code> are not to be \"touched\".\n",
      "    \n",
      "The specific laff function you will want to use is some subset of\n",
      "<ul>\n",
      "<li> <code> laff.gemv( trans, alpha, A, x, beta, y ) </code> which computes $ y := \\alpha A x + \\beta y $ or $ y := \\alpha A^T x + \\beta y $ depending on whether <code> trans = 'No transpose' </code> or <code> trans = 'Transpose' </code> </li>\n",
      "    <li> <code> laff.ger( alpha, x, y, A ) </code> which computes the rank-1 update (adds a multiple of an outer product to a matrix)\n",
      "$ A := \\alpha x y^T + A $. </li>\n",
      "    <li> <code> laff.axpy( alpha, x, y ) </code></li>\n",
      "    <li>    <code> laff.dots( x, y, alpha ) </code></li>\n",
      "</ul>\n",
      "\n",
      "<font color=red> Hint:</font>  If you multiply with $ U_{00} $, you will want to use <code> np.triu( U00 ) </code> to make sure you don't compute with the nonzeroes below the diagonal.\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": [
      "# insert code here\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.matrix( np.copy( Cold ) )               # restore C \n",
      "\n",
      "Trtrmm_uu_unb_var1( U, R, C )\n",
      "\n",
      "# compute it using numpy *.  This is a little complex, since we want to make sure we\n",
      "#   don't change the lower triangular part of C\n",
      "# Cref = np.tril( Cold, -1 ) + np.triu( U ) * np.triu( R )\n",
      "Cref = np.triu( U ) * np.triu( R ) + np.tril( Cold, -1 )\n",
      "print( 'C - Cref' )\n",
      "print( C - Cref )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In theory, you should get a matrix of (approximately) all zeroes."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Watch the algorithm at work!"
     ]
    },
    {
     "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."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}