psets/pset2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 18.S096 pset 2\n",
    "\n",
    "Due Wednesday, January 24, 2pm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1: Backsubstitution\n",
    "\n",
    "An [upper-triangular system of equations](https://en.wikipedia.org/wiki/Triangular_matrix) can be solved efficiently ($O(n^2)$ time) by *backsubstitution*: starting at the last row (= one equation in one unknown) and working upwards, eliminating one variable at a time.\n",
    "\n",
    "Here is \"textbook\" code that implements backsubstition for an upper-triangular matrix `U` and a matrix or vector `B` of right-hand sides, i.e. solving $UX=B$ for $X$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "backsolve (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function backsolve(U::UpperTriangular, B::Union{AbstractMatrix,AbstractVector})\n",
    "    m = LinAlg.checksquare(U)\n",
    "    m == size(B,1) || throw(DimensionMismatch(\"matrix sizes don't match\"))\n",
    "    n = size(B, 2)\n",
    "    X = similar(B, typeof(zero(eltype(B))/oneunit(eltype(U)))) # allocate result array\n",
    "    for j = 1:n, i = m:-1:1\n",
    "        s = zero(eltype(X))\n",
    "        for k = i+1:m\n",
    "            s += U[i,k] * X[k,j]\n",
    "        end\n",
    "        X[i,j] = (B[i,j] - s) / U[i,i]\n",
    "    end\n",
    "    return X\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that it gives the same answer as `\\` in a couple of cases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 0.567182  0.851563  0.0675331\n",
       "  ⋅        0.664635  0.463687 \n",
       "  ⋅         ⋅        0.771664 "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = UpperTriangular(rand(3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.0172276\n",
       " 0.094541 \n",
       " 0.15564  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = rand(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.00405888\n",
       " 0.00153189\n",
       " 0.201694  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U \\ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.00405888\n",
       " 0.00153189\n",
       " 0.201694  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "backsolve(U, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = rand(3,4)\n",
    "U \\ B ≈ backsolve(U, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But it is **vastly slower**, especially for a matrix of right-hand sides, for matrices of significant size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /Users/stevenj/.julia/lib/v0.6/BenchmarkTools.ji for module BenchmarkTools.\n",
      "\u001b[39m"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  874.866 μs (1 allocation: 7.94 KiB)\n",
      "  311.113 μs (2 allocations: 7.95 KiB)\n",
      "  969.399 ms (2 allocations: 7.63 MiB)\n",
      "  9.551 ms (3 allocations: 7.63 MiB)\n"
     ]
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "U = UpperTriangular(rand(1000,1000))\n",
    "b = rand(1000)\n",
    "B = rand(1000,1000)\n",
    "@btime backsolve($U,$b)\n",
    "@btime (\\)($U,$b)\n",
    "@btime backsolve($U,$B)\n",
    "@btime (\\)($U,$B);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3x slower for a vector, but 100x slower for an $m\\times m$ matrix of right-hand sides, where backsubstitution is $O(m^3)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assignment: Speed up `backsolve`\n",
    "\n",
    "* Things like `@inbounds` and `@simd` are obviously the first things to try.\n",
    "\n",
    "* Think about a blocked or recursive/cache-oblivious algorithm.  If you divide the matrices into half-size blocks, can you write down a backsubstitution algorithm acting on the blocks?  For example, divide the matrices as $$B = UX = \\begin{pmatrix} U_1 & U_2 \\\\ 0 & U_3 \\end{pmatrix} \\begin{pmatrix} X_1 \\\\ X_2 \\end{pmatrix} = \\begin{pmatrix} B_1 \\\\ B_2 \\end{pmatrix}$$ and see if the blocked algorithm (backsolving recursively on the blocks, with some base case) speeds things up.\n",
    "\n",
    "* Can you parallelize it, especially for a large matrix?  (i.e. solve for different columns of `X` in parallel.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}