{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 18.06 Pset 2 - Solutions\n",
    "\n",
    "Due Wed. 9/19 at 10:55am by electronic submission."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1 (10 points)\n",
    "\n",
    "For the matrix\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix}\n",
    "1 & 2 & 1 & -1 & 1 \\\\\n",
    "1 & 0 & 1 & -1 & 0 \\\\\n",
    "1 & 2 & 1 & 0 & 0 \\\\\n",
    "1 & 1 & 0 &  0 & 0 \\\\\n",
    "-1 & 0 & 0 &  0 & 0 \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "and the vector $b = \\begin{pmatrix} 0\\\\0\\\\0\\\\0\\\\1 \\end{pmatrix}$:\n",
    "\n",
    "**(a)** Show *hand calculations* to compute $x = A^{-1} b$.  (Hint: remember what I said in class about computing inverses of matrices.)\n",
    "\n",
    "**(b)** Compute $A^{-1}$ in Julia with `inv(A)` and compare it to your solution $x$.  Explain why your $x$ appears in $A^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "**(a)** To calculate $x=A^{-1}b$, we just need to find the solution of $Ax = b$ (we don't want to explicitly calculate the inverse matrix if we can help it!)  Although the matrix $A$ is neither upper nor lower triangular, we could readily switch rows to put it in lower triangular form. However, we can already solve the linear system $Ax = b$ by forward substitution in its current form. Let $x=\\begin{pmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\\\ x_5\\end{pmatrix}$. Then:\n",
    "\\begin{align}\n",
    "-x_1 = 1 &\\implies x_1 = -1\\\\\n",
    "x_1 + x_2 = 0 &\\implies x_2 = 1\\\\\n",
    "x_1 + 2x_2 + x_3 = 0 &\\implies x_3 = -1\\\\\n",
    "x_1 + x_3 - x_4 = 0 &\\implies x_4 = -2\\\\\n",
    "x_1 +2x_2 + x_3 -x_4 + x_5 = 0 &\\implies x_5 = -2.\n",
    "\\end{align}\n",
    "So that $x = \\begin{pmatrix} -1 \\\\ 1 \\\\ -1 \\\\ -2 \\\\ -2 \\end{pmatrix}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** We can use Julia to explicitly calculate $A^{-1}$ (remembering that in general this is a bad idea):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Int64,2}:\n",
       "  1  2  1  -1  1\n",
       "  1  0  1  -1  0\n",
       "  1  2  1   0  0\n",
       "  1  1  0   0  0\n",
       " -1  0  0   0  0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ 1 2 1 -1 1\n",
    "      1 0 1 -1 0\n",
    "      1 2 1  0 0\n",
    "      1 1 0  0 0\n",
    "     -1 0 0  0 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.0   0.0  0.0   0.0  -1.0\n",
       " 0.0   0.0  0.0   1.0   1.0\n",
       " 0.0   0.0  1.0  -2.0  -1.0\n",
       " 0.0  -1.0  1.0  -2.0  -2.0\n",
       " 1.0  -1.0  0.0  -2.0  -2.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ainv = inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can double check that our $x$ is correct by then explicitly calculating $x = A^{-1}b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = [0\n",
    "     0\n",
    "     0\n",
    "     0\n",
    "     1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " -1.0\n",
       "  1.0\n",
       " -1.0\n",
       " -2.0\n",
       " -2.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Ainv*b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is exactly what we found in part **(a)**. However, notice that $x$ corresponds to the entries in the fifth column of $A^{-1}$. This should not surprise us. Recall from the end of lecture 3 that we can compute $A^{-1}$ by repeatedly solving $Ax = b$, where $b$ is each of the columns of the identity matrix $I$. Each of the corresponding solutions will then be the corresponding column of $A^{-1}$. We solved $Ax=b$ for $b$ being the fifth column of the identity matrix, so we knew already that our solution $x$ would be the fifth column of $A^{-1}$.\n",
    "\n",
    "More generally, if you multiply *any* matrix by a column of $I$ then you get the corresponding column of that matrix.  Here, if you multiply any $5\\times 5$ matrix by this $b$ then you get its 5th column, and hence $A^{-1}b$ gives the fifth column of $A^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2 (15 points)\n",
    "\n",
    "Consider the matrix $$A = \\begin{pmatrix} 1 & 4 & 1 \\\\ 1 & 2 & -1 \\\\ 3 & 14 & 6 \\end{pmatrix}$$ from lecture 4.\n",
    "\n",
    "In this problem, we will consider transforming this matrix by a sequence of invertible linear **column operations** — multiplying *columns* by scalars and adding/subtracting them.\n",
    "\n",
    "**(a)** Come up with column operations that change the first row from $\\begin{pmatrix} 1 & 4 & 1\\end{pmatrix}$ to $\\begin{pmatrix} 1 & 0 & 0\\end{pmatrix}$, i.e. that put zeros to the *right* of the first pivot.   Express these operations in terms of multipling $A$ on the left or right by some matrix $E$.   Give $E$ and $E^{-1}$.\n",
    "\n",
    "(Note that your operations must be invertible, like Gaussian-elimination steps — no fair just multipling the second and third columns by zero, since this would lead to a non-invertible $E$.)\n",
    "\n",
    "**(b)** If we do a *sequence* of column operations that transform $A$ into $I$, and then do the *same* column operations to the $3\\times 3$ identity matrix $I$, what is the resulting matrix?\n",
    "\n",
    "**(c)** Carry out the process from (b): perform column operations on $A$ that transform it into $I$, and perform the *same* column operations on $I$.   Do them both at once by \"augmenting\" $A$ in a certain way (maybe not the usual way).   Check your result against the lecture-4 notes or with Julia to see that you got the expected result from (b)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "**(a)** To turn the first row of $A$ from $\\begin{pmatrix} 1 & 4 & 1\\end{pmatrix}$ to $\\begin{pmatrix} 1 & 0 & 0\\end{pmatrix}$, we can subtract 4 times the first column from the second column, and then subtract the first column from the third column. Recall that column operations can be described by multipling on the right by a matrix $E$, so that $AE$ has the correct first row. The matrix $E$ that corresponds to these column operations is\n",
    "\\begin{align}\n",
    "E = \\begin{pmatrix} 1 & -4 & -1 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{pmatrix}\n",
    "\\end{align}\n",
    "The inverse of these column operations is to add 4 times the first column to the second column, and add the first column to the third column. This allows us to construct $E^{-1}$:\n",
    "\\begin{align}\n",
    "E^{-1} = \\begin{pmatrix} 1 & 4 & 1 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** Suppose we do a sequence of row operations to $A$ to transform it into the identity. We can express this as multiplication on the right by a matrix $R$, so that $AR = I$. But notice that this matrix must be the inverse of $A$, i.e. $R=A^{-1}$!  Applying the same column operations to $I$ is then equivalent to multiplying $I$ on the right by $A^{-1}$, to give $IA^{-1} = A^{-1}$. So applying the column operations to $I$ will give us the the inverse of $A$. This process will work for any non-singular square matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** We can carry out column operations on $A$ that transform it into $I$, and perform the *same* column operations on $I$ simultaneously by augmenting $A$ into the following $6\\times 3$ matrix:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 4 & 1 \\\\ 1 & 2 & -1 \\\\ 3 & 14 & 6 \\\\ 1 & 0 & 0\\\\ 0 & 1 & 0\\\\ 0 & 0 & 1\\end{pmatrix}.\n",
    "\\end{align}\n",
    "applying the column operations from part **(a)** yields:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 1 & -2 & -2 \\\\ 3 & 2 & 3 \\\\ 1 & -4 & -1\\\\ 0 & 1 & 0\\\\ 0 & 0 & 1\\end{pmatrix}.\n",
    "\\end{align}\n",
    "We can then add $1/2$ times the second column to the first column, and subtract the second column from the third column, to yield:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & -2 & 0 \\\\ 4 & 2 & 1 \\\\ -1 & -4 & 3\\\\ 1/2 & 1 & -1\\\\ 0 & 0 & 1\\end{pmatrix}.\n",
    "\\end{align}\n",
    "We can then divide the second column by -2:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 4 & -1 & 1 \\\\ -1 & 2 & 3\\\\ 1/2 & -1/2 & -1\\\\ 0 & 0 & 1\\end{pmatrix}.\n",
    "\\end{align}\n",
    "Finally we can subtract 4 times column 3 from column 1, and add column 3 to column 2, to obtain:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\\\ -13 & 5 & 3\\\\ 9/2 & -3/2 & -1\\\\ -4 & 1 & 1\\end{pmatrix}.\n",
    "\\end{align}\n",
    "We can check our result by letting $R = \\begin{pmatrix}  -13 & 5 & 3\\\\ 9/2 & -3/2 & -1\\\\ -4 & 1 & 1\\end{pmatrix}$. We can then calculate $AR$ in Julia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [1  4  1\n",
    "     1  2 -1\n",
    "     3 14  6];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = [ -13  5    3\n",
    "      4.5 -1.5 -1\n",
    "      -4   1    1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 1.0  0.0  0.0\n",
       " 0.0  1.0  0.0\n",
       " 0.0  0.0  1.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So indeed applying the same column operations to $I$ yields the matrix $R$, which we can multiply $A$ on the right by to yield $I$, i.e. the inverse of $A$. We can also use to Julia to calculate $A^{-1}$ explicity, and we find that it is indeed the same as $R$. Note that $A^{-1}$ was also given in the lecture notes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -13.0   5.0   3.0\n",
       "   4.5  -1.5  -1.0\n",
       "  -4.0   1.0   1.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3 (10 points)\n",
    "\n",
    "From Strang, section 2.6: Compute $L$ and $U$ for the following symmetric matrix $A$:\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix} a & a & a & a \\\\ a & b & b & b \\\\ a & b & c & c \\\\ a & b & c & d \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "and find four conditions on the numbers $a,b,c,d$ to get $A=LU$ with four pivots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "We firstly perform Gaussian elimination on $A$ to obtain $U$. We subtract row 1 from rows 2, 3 & 4 to obtain:\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} a & a & a & a \\\\ 0 & b-a & b-a & b-a \\\\ 0 & b-a & c-a & c-a \\\\ 0 & b-a & c-a & d-a \\end{pmatrix}\n",
    "\\end{align}\n",
    "We then subtract row 2 from rows 3 and 4 to obtain\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} a & a & a & a \\\\ 0 & b-a & b-a & b-a \\\\ 0 & 0 & c-b & c-b \\\\ 0 & 0 & c-b & d-b \\end{pmatrix}\n",
    "\\end{align}\n",
    "Finally we subtract row 3 from row 4 to obtain \n",
    "\\begin{align}\n",
    "\\begin{pmatrix} a & a & a & a\\\\ 0 & b-a & b-a & b-a \\\\ 0 & 0 & c-b & c-b \\\\ 0 & 0 & 0 & d-c \\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Hence $U = \\begin{pmatrix} a & a & a & a \\\\ 0 & b-a & b-a & b-a \\\\ 0 & 0 & c-b & c-b \\\\ 0 & 0 & 0 & d-c \\end{pmatrix}$\n",
    "\n",
    "One way to compute $L$ would be to find the elimination matrices that correspond to each of these steps, invert them, and then multiply them together. But this is tedious and uneccessary! Remember from lecture 5 (see 'LU factorization for real') that to find $L$ we just need to keep track of the multipliers at each step of the Gaussian elimination (they are all $1$). We can then immediately write down that: $L = \\begin{pmatrix} 1 & 0 & 0 & 0\\\\ 1 & 1 & 0 & 0 \\\\ 1 & 1 & 1 & 0 \\\\ 1 & 1 & 1 & 1 \\end{pmatrix}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have found $A=LU$, where \n",
    "\\begin{align}\n",
    "L = \\begin{pmatrix} 1 & 0 & 0 & 0\\\\ 1 & 1 & 0 & 0 \\\\ 1 & 1 & 1 & 0 \\\\ 1 & 1 & 1 & 1 \\end{pmatrix} \\; \\text{and} \\;\n",
    "U = \\begin{pmatrix} a & a & a & a \\\\ 0 & b-a & b-a & b-a \\\\ 0 & 0 & c-b & c-b \\\\ 0 & 0 & 0 & d-c \\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "The conditions to have four pivots are that $a\\neq 0, a \\neq b, b \\neq c$ and $c\\neq d$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4 (10 points)\n",
    "\n",
    "From Strang, section 2.6.  Consider $$L = \\begin{pmatrix} 1 & 0 & 0 \\\\ a & 1 & 0 \\\\ b & c & 1 \\end{pmatrix}$$ for some numbers $a,b,c$.\n",
    "\n",
    "**(a)** When you perform the usual Gaussian elimination steps to $L$, what matrix will you obtain?\n",
    "\n",
    "**(b)** If you apply the *same* row operations to $I$, what matrix will you get?\n",
    "\n",
    "**(c)** If you apply the *same* row operations to $LB$ for some $3\\times n$ matrix $B$, what will you get?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "**(a)** The usual Gaussian elimination method would suggest we firstly subtract $a$ times row 1 from row 2, and then subtract $b$ times row 1 from row 3 to give:\n",
    "$$\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0\\\\ 0 & c & 1\\end{pmatrix}$$\n",
    "We can then subtract $c$ times row 2 from row 3 to get \n",
    "$$\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0\\\\ 0 & 0 & 1\\end{pmatrix}$$\n",
    "So Gaussian elimination on $L$ yields the identity matrix.\n",
    "\n",
    "**(b)** Applying the same sequence of row operations to the identity matrix yields\n",
    "\\begin{align}\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0\\\\ 0 & 0 & 1\\end{pmatrix} \\rightarrow \n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ -a & 1 & 0\\\\ -b & 0 & 1\\end{pmatrix} \\rightarrow \n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ -a & 1 & 0\\\\ ac-b & -c & 1\\end{pmatrix}\n",
    "\\end{align}\n",
    "Remeber that if you do row operations to $L$ that give $I$, then those row operations are exactly equivalent to multiplying by $L^{-1}$ on the left.  So when you do these row operations to $I$ you will get $L^{-1}$, i.e. $L^{-1} = \\begin{pmatrix} 1 & 0 & 0 \\\\ -a & 1 & 0\\\\ ac-b & -c & 1\\end{pmatrix}$\n",
    "\n",
    "**(c)** Applying the same row operations to $LB$ is equivalent to multiplying on the left by the matrix $L^{-1}$, i.e. $L^{-1}(LB)$. But we can use the associativity of matrix multiplication: $L^{-1}(LB) = (L^{-1}L)B$. But $L^{-1}L=I$, and so applying these row operations to the product $LB$ will just yield the matrix $B$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5 (10 points)\n",
    "\n",
    "Consider the following matrices:\n",
    "\n",
    "$$\n",
    "U = \\begin{pmatrix} 1 & 1 & -1 \\\\ 0 & 1 & 2 \\\\ 0 & 0 & 1 \\end{pmatrix}, \\;\n",
    "L = \\begin{pmatrix} 1 & 0 & 0 \\\\ -1 & 1 & 0 \\\\ -2 & 1 & 1 \\end{pmatrix}, \\;\n",
    "B = \\begin{pmatrix} 1 & 2 & 3 \\\\ 3 & 2 & 1 \\\\ 1 & 0 & 1 \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Let $A = U B^{-1} L$.\n",
    "\n",
    "**(a)** Compute the *second column* of $A^{-1}$.  (If you think about it, you can do it *without inverting any matrices*.)\n",
    "\n",
    "**(b)** Check your answer by explicitly computing $A^{-1}$ in Julia."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "**(a)** Recall that the second column of $A^{-1}$ will be the same as the column vector $x$ which solves $Ax = b$ for $b=\\begin{pmatrix} 0 \\\\ 1 \\\\ 0 \\end{pmatrix}$.\n",
    "\n",
    "We want to find $A^{-1}b$, and $A^{-1} = L^{-1} B U^{-1}$. So, we just need to multiply $U^{-1}b$ via backsubstitution, then multiply by $B$, then multiply by $L^{-1}$ by forward-substition.\n",
    "\n",
    "$Ax=b$ is the same as $UB^{-1}Lx = b$. To solve this system, let's write $y = B^{-1}Lx$ and firstly solve $Uy = b$ by back substitution. If $ y = \\begin{pmatrix} y_1 \\\\ y_2 \\\\ y_3 \\end{pmatrix}$, then\n",
    "\\begin{align}\n",
    "y_3 = 0 &\\implies y_3 = 0\\\\\n",
    "y_2 + 2y_3 = 1 &\\implies y_2 = 1\\\\\n",
    "y_1 + y_2 -y_3 = 0 &\\implies y_1 = -1\n",
    "\\end{align}\n",
    "\n",
    "So now we just need to solve $B^{-1}Lx = y$, or equivalently, $Lx = By$. \n",
    "\n",
    "We can calculate $By = \\begin{pmatrix} 1 & 2 & 3 \\\\ 3 & 2 & 1 \\\\ 1 & 0 & 1 \\end{pmatrix}\\begin{pmatrix} -1 \\\\ 1 \\\\ 0\\end{pmatrix} = \\begin{pmatrix} 1 \\\\ -1 \\\\-1\\end{pmatrix}$.\n",
    "\n",
    "So we now solve $Lx = \\begin{pmatrix} 1 \\\\ -1 \\\\-1\\end{pmatrix}$ by forward substitution. If $x = \\begin{pmatrix} x_1 \\\\ x_2 \\\\ x_3 \\end{pmatrix}$, then \n",
    "\\begin{align}\n",
    "x_1 = 1 &\\implies x_1 = 1\\\\\n",
    "-x_1 + x_2 = -1 &\\implies x_2 = 0\\\\\n",
    "-2x_1 + x_2 + x_3 = -1 &\\implies x_3 =1\n",
    "\\end{align}\n",
    "\n",
    "So the solution of $Ax = b$, i.e. the second column of $A^{-1}$, is $x = \\begin{pmatrix} 1 \\\\ 0 \\\\ 1 \\end{pmatrix}$.\n",
    "\n",
    "**(b)** We can check this calculation by explicitly calculating $A^{-1}$ in Julia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  1.0   1.0  2.0\n",
       "  4.0  -0.0  8.0\n",
       " -1.0   1.0  0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = [1 1 -1\n",
    "     0 1 2\n",
    "     0 0 1];\n",
    "L = [1 0 0\n",
    "    -1 1 0\n",
    "    -2 1 1];\n",
    "B = [1 2 3\n",
    "     3 2 1\n",
    "     1 0 1];\n",
    "A = U*(B\\L);\n",
    "inv(A) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 6 (10 points)\n",
    "\n",
    "Suppose you have a 5×5 matrix of the form \n",
    "\n",
    "$$A = \\begin{pmatrix}\n",
    "\\star & \\star & 0 & 0 & 0 \\\\\n",
    "\\star & \\star & \\star & 0 & 0 \\\\\n",
    "0 & \\star & \\star & \\star & 0 \\\\\n",
    "0 & 0 & \\star & \\star & \\star \\\\\n",
    "0 & 0 & 0 & \\star & \\star\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "where \"$\\star$\" denotes nonzero entries (**not necessarily all the same numbers**).   This is called a [tridiagonal matrix](https://en.wikipedia.org/wiki/Tridiagonal_matrix).\n",
    "\n",
    "The following code constructs a random tridiagonal matrix in Julia, using a special `Tridiagonal` type that allows Julia to take advantage of this structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Tridiagonal{Float64}:\n",
       " -1.51486  -1.22025     ⋅           ⋅          ⋅      \n",
       " -1.90274  -1.69823    0.0605162    ⋅          ⋅      \n",
       "   ⋅       -0.875176  -1.40127     0.526583    ⋅      \n",
       "   ⋅         ⋅        -0.0994365  -1.48919   -0.701886\n",
       "   ⋅         ⋅          ⋅          0.143034   1.15643 "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Tridiagonal(randn(4), randn(5), randn(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *inverse* of a tridiagonal matrix has no special pattern of nonzero entries in general:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " -5.65997     3.98062     0.167448     0.0628758   0.0381619\n",
       "  6.207      -4.9417     -0.207877    -0.0780566  -0.0473757\n",
       " -3.77604     3.00629    -0.568656    -0.213527   -0.129598 \n",
       "  0.267743   -0.213163    0.040321    -0.697935   -0.423605 \n",
       " -0.0331159   0.0263652  -0.00498712   0.0863245   0.917122 "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the *LU factorization* of a tridiagonal matrix is very special.\n",
    "\n",
    "**(a)** Assuming no row swaps are required and the matrix is nonsingular, what pattern of nonzero entries do you generically expect to see in the $L$ and $U$ factors of a matrix $A$ of this form, and why?\n",
    "\n",
    "Check your answer against your random 5×5 $A$ in Julia.\n",
    "\n",
    "**(b)** If $A$ is an $m \\times m$ tridiagonal matrix (i.e. same pattern of zeros, but extended to an arbitrary size), how does the count of scalar arithmetic operations (±,×,÷) to compute the L, U factors (i.e. to perform elimination) scale with $m$?  You need not give an exact count, just say whether it is roughtly proportional to $m$, $m^2$, $m^3$, etcetera for large $m$.)\n",
    "\n",
    "You need not count operations on numbers that are *guaranteed* to be zero for *any* tridiagonal matrix.  The computer can store *just* the nonzero entries of the matrix and only operate on those."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "**(a)** Suppose we have a tridiagonal matrix $A$. To obtain the $LU$ factorization, we must perform Gaussian elimination. To transform $A$ into upper triangular form, we only need to eliminate the entries on the lower diagonal. This can be done by subtracting some multiple of the first row from the second row, then subtracting some multiple of the second row from the third row, and so. This means that the matrix $U$ will be tridiagonal, but the lower diagonal will all be zeros. Since the row operations required to perform this Gaussian elimination only require us at each step to to perform row operations on adjacent rows, the matrix $L$, which describes the inverse of the row operations required to perform Gaussian elimination will also be tridiagonal, but now the upper diagonal will all be zeros. \n",
    "\n",
    "We say that $L$ and $U$ are *bidiagonal* matrices.\n",
    "\n",
    "This can be confirmed by calculating $L$ and $U$ for a random $5\\times 5$ tridiagonal matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 1.0      0.0      0.0         0.0        0.0\n",
       " 1.25605  1.0      0.0         0.0        0.0\n",
       " 0.0      5.28665  1.0         0.0        0.0\n",
       " 0.0      0.0      0.0577718   1.0        0.0\n",
       " 0.0      0.0      0.0        -0.0941255  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " -1.51486  -1.22025    0.0         0.0        0.0     \n",
       "  0.0      -0.165545   0.0605162   0.0        0.0     \n",
       "  0.0       0.0       -1.72119     0.526583   0.0     \n",
       "  0.0       0.0        0.0        -1.51961   -0.701886\n",
       "  0.0       0.0        0.0         0.0        1.09037 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L, U = lu(A, Val{false})\n",
    "display(L)\n",
    "display(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** To find the $L, U$ factorization we have to perform Gaussian elimination as described in part **(a)**. Each step of the Guassian elimination will require us to multiply a row by an appropriate factor, which will require 2 multiplications on nonzero entries. We will then subtract this row from the row below, which will require 2 non trivial subtractions. We have to do this (m-1) times to transform the matrix into upper triangular form. The total number of nontrivial operations to perform Gaussian elimination is therefore 4*(m-1). \n",
    "\n",
    "Gaussian elimination therefore requires a number of operations that is roughly proportional to $m$. So obtaining the LU factorization of a tridiagonal matrix has *linear complexity*. \n",
    "\n",
    "(Note that performing Gaussian elimination witohut row exchanges for a general matrix $A$ will require a number of operations that is roughly proportional to $m^3$)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}