Chapter 2 - Basic Matrix Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import scipy as sp\n",
    "import sympy as sy\n",
    "sy.init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\" # display multiple results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_expr(expr, num_digits):\n",
    "    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(sy.Number)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Matrix Operations</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix _addition_ operations are straightforward:\n",
    "1. $A+ B= B+ A$\n",
    "2. $(A+B)+ C=A+(B+C)$\n",
    "3. $c(A+B)=cA+cB$\n",
    "4. $(c+d)A=cA+c{D}$\n",
    "5. $c(dA)=(cd)A$\n",
    "6. $A+{0}=A$, where ${0}$ is the zero matrix\n",
    "7. For any $A$, there exists an $- A$, such that $ A+(- A)=0$.\n",
    "\n",
    "They are as obvious as it looks, so no proofs are provided. And the matrix _multiplication_ properties are:\n",
    "1. $ A({BC})=({AB}) C$\n",
    "2. $c({AB})=(cA)B=A(cB)$\n",
    "3. $A(B+ C)={AB}+{AC}$\n",
    "4. $(B+C)A={BA}+{CA}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we need to differentiate two kinds of multiplication, **Hadamard multiplication** (element-wise multiplication) and **matrix multiplication**: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 2], [3, 4]])\n",
    "B = np.array([[5, 6], [7, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5, 12],\n",
       "       [21, 32]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*B # this is Hadamard elementwise product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[19, 22],\n",
       "       [43, 50]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@B # this is matrix product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We show explicitly the matrix multiplication rule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "43"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "22"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(A[0,:]*B[:,0]) # (1, 1)\n",
    "np.sum(A[1,:]*B[:,0]) # (2, 1)\n",
    "np.sum(A[0,:]*B[:,1]) # (1, 2)\n",
    "np.sum(A[1,:]*B[:,1]) # (2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> SymPy Demonstration: Addition </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define all the letters as symbols in case we might use them repetitively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z = sy.symbols('a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z', real = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 a & 2 b & 2 c\\\\2 d & 2 e & 2 f\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡2⋅a  2⋅b  2⋅c⎤\n",
       "⎢             ⎥\n",
       "⎣2⋅d  2⋅e  2⋅f⎦"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\\\0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0  0  0⎤\n",
       "⎢       ⎥\n",
       "⎣0  0  0⎦"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[a, b, c], [d, e, f]])\n",
    "A + A\n",
    "A - A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a + g & b + h & c + i\\\\d + j & e + k & f + l\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a + g  b + h  c + i⎤\n",
       "⎢                   ⎥\n",
       "⎣d + j  e + k  f + l⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a - g & b - h & c - i\\\\d - j & e - k & f - l\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a - g  b - h  c - i⎤\n",
       "⎢                   ⎥\n",
       "⎣d - j  e - k  f - l⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = sy.Matrix([[g, h, i], [j, k, l]])\n",
    "A + B\n",
    "A - B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> SymPy Demonstration: Multiplication </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix multiplication rules can be clearly understood by using symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a & b & c\\\\d & e & f\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a  b  c⎤\n",
       "⎢       ⎥\n",
       "⎣d  e  f⎦"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}g & h & i\\\\j & k & l\\\\m & n & o\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡g  h  i⎤\n",
       "⎢       ⎥\n",
       "⎢j  k  l⎥\n",
       "⎢       ⎥\n",
       "⎣m  n  o⎦"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[a, b, c], [d, e, f]])\n",
    "B = sy.Matrix([[g, h, i], [j, k, l], [m, n, o]])\n",
    "A\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a g + b j + c m & a h + b k + c n & a i + b l + c o\\\\d g + e j + f m & d h + e k + f n & d i + e l + f o\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a⋅g + b⋅j + c⋅m  a⋅h + b⋅k + c⋅n  a⋅i + b⋅l + c⋅o⎤\n",
       "⎢                                                 ⎥\n",
       "⎣d⋅g + e⋅j + f⋅m  d⋅h + e⋅k + f⋅n  d⋅i + e⋅l + f⋅o⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AB = A*B; AB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Commutability </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix multiplication usually do not commute, such that ${AB} \\neq {BA}$. For instance, consider $ A$ and $ B$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}23 & 13\\\\51 & 29\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡23  13⎤\n",
       "⎢      ⎥\n",
       "⎣51  29⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}36 & 44\\\\13 & 16\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡36  44⎤\n",
       "⎢      ⎥\n",
       "⎣13  16⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[3, 4], [7, 8]])\n",
    "B = sy.Matrix([[5, 3], [2, 1]])\n",
    "A*B\n",
    "B*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we find a commutable matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a e + b g & a f + b h\\\\c e + d g & c f + d h\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a⋅e + b⋅g  a⋅f + b⋅h⎤\n",
       "⎢                    ⎥\n",
       "⎣c⋅e + d⋅g  c⋅f + d⋅h⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a e + c f & b e + d f\\\\a g + c h & b g + d h\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a⋅e + c⋅f  b⋅e + d⋅f⎤\n",
       "⎢                    ⎥\n",
       "⎣a⋅g + c⋅h  b⋅g + d⋅h⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[a, b], [c, d]])\n",
    "B = sy.Matrix([[e, f], [g, h]])\n",
    "A*B\n",
    "B*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To show ${AB} = {BA}$, we need to prove ${AB} - {BA} = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}b g - c f & a f - b e + b h - d f\\\\- a g + c e - c h + d g & - b g + c f\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡      b⋅g - c⋅f         a⋅f - b⋅e + b⋅h - d⋅f⎤\n",
       "⎢                                             ⎥\n",
       "⎣-a⋅g + c⋅e - c⋅h + d⋅g       -b⋅g + c⋅f      ⎦"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = A*B - B*A; M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That means\n",
    "\\begin{align}\n",
    "b g - c f&=0 \\\\\n",
    " a f - b e + b h - d f&=0\\\\\n",
    "- a g + c e - c h + d g&=0 \\\\\n",
    "- b g + c f&=0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we treat $a, b, c, d$ as coefficients of the system, we and extract an augmented matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & - c & b & 0\\\\- b & a - d & 0 & b\\\\c & 0 & - a + d & - c\\\\0 & c & - b & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0    -c      b     0 ⎤\n",
       "⎢                     ⎥\n",
       "⎢-b  a - d    0     b ⎥\n",
       "⎢                     ⎥\n",
       "⎢c     0    -a + d  -c⎥\n",
       "⎢                     ⎥\n",
       "⎣0     c      -b    0 ⎦"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_aug = sy.Matrix([[0, -c, b, 0], [-b, a-d, 0, b], [c, 0, d -a, -c], [0, c, -b, 0]]); A_aug"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform **Gaussian-Jordon** elimination till row reduced formed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & - \\frac{a}{c} + \\frac{d}{c} & -1\\\\0 & 1 & - \\frac{b}{c} & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right], \\  \\left( 0, \\  1\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡        a   d    ⎤        ⎞\n",
       "⎜⎢1  0  - ─ + ─  -1⎥        ⎟\n",
       "⎜⎢        c   c    ⎥        ⎟\n",
       "⎜⎢                 ⎥        ⎟\n",
       "⎜⎢        -b       ⎥        ⎟\n",
       "⎜⎢0  1    ───    0 ⎥, (0, 1)⎟\n",
       "⎜⎢         c       ⎥        ⎟\n",
       "⎜⎢                 ⎥        ⎟\n",
       "⎜⎢0  0     0     0 ⎥        ⎟\n",
       "⎜⎢                 ⎥        ⎟\n",
       "⎝⎣0  0     0     0 ⎦        ⎠"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_aug.rref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general solution is \n",
    "\\begin{align}\n",
    "e - \\frac{a-d}cg - h &=0\\\\\n",
    "f - \\frac{b}{c} & =0\\\\\n",
    "g &= free\\\\\n",
    "h & =free\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "if we set coefficients $a = 10, b = 12, c = 20, d = 8$, or $ A = \\left[\\begin{matrix}10 & 12\\\\20 & 8\\end{matrix}\\right]$ then general solution becomes\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "e - .1g - h &=0\\\\\n",
    "f - .6 & =0\\\\\n",
    "g &= free\\\\\n",
    "h & =free\n",
    "\\end{align}\n",
    "Then try a special solution when $g = h = 1$\n",
    "\\begin{align}\n",
    "e  &=1.1\\\\\n",
    "f & =.6\\\\\n",
    "g &=1 \\\\\n",
    "h & =1\n",
    "\\end{align}\n",
    "And this is a **commutable matrix of $A$**, we denote $ C$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.1 & 0.6\\\\1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1.1  0.6⎤\n",
       "⎢        ⎥\n",
       "⎣ 1    1 ⎦"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = sy.Matrix([[1.1, .6], [1, 1]]);C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see that ${AC}={CA}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}23.0 & 18.0\\\\30.0 & 20.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡23.0  18.0⎤\n",
       "⎢          ⎥\n",
       "⎣30.0  20.0⎦"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}23.0 & 18.0\\\\30 & 20\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡23.0  18.0⎤\n",
       "⎢          ⎥\n",
       "⎣ 30    20 ⎦"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[10, 12], [20, 8]])\n",
    "A*C\n",
    "C*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Transpose of Matrices </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix $A_{n\\times m}$ and its transpose is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[1, 4],\n",
       "       [2, 5],\n",
       "       [3, 6]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3], [4, 5, 6]]); A\n",
    "A.T # transpose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 2 & 3\\\\4 & 5 & 6\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  2  3⎤\n",
       "⎢       ⎥\n",
       "⎣4  5  6⎦"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4\\\\2 & 5\\\\3 & 6\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  4⎤\n",
       "⎢    ⎥\n",
       "⎢2  5⎥\n",
       "⎢    ⎥\n",
       "⎣3  6⎦"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[1, 2, 3], [4, 5, 6]]); A\n",
    "A.transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The properties of transpose are\n",
    "1. $(A^T)^T$\n",
    "2. $(A+B)^T=A^T+B^T$\n",
    "3. $(cA)^T=cA^T$\n",
    "4. $(AB)^T=B^TA^T$\n",
    "\n",
    "We can show why the last property holds with SymPy, define $A$ and $B$, multiply them, then transpose, that means $(AB)^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a g + b j & c g + d j & e g + f j\\\\a h + b k & c h + d k & e h + f k\\\\a i + b l & c i + d l & e i + f l\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a⋅g + b⋅j  c⋅g + d⋅j  e⋅g + f⋅j⎤\n",
       "⎢                               ⎥\n",
       "⎢a⋅h + b⋅k  c⋅h + d⋅k  e⋅h + f⋅k⎥\n",
       "⎢                               ⎥\n",
       "⎣a⋅i + b⋅l  c⋅i + d⋅l  e⋅i + f⋅l⎦"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[a, b], [c, d], [e, f]])\n",
    "B = sy.Matrix([[g, h, i], [j, k, l]])\n",
    "AB = A*B\n",
    "AB_tr = AB.transpose(); AB_tr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transpose each of them, then multiply, that means $B^TA^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a g + b j & c g + d j & e g + f j\\\\a h + b k & c h + d k & e h + f k\\\\a i + b l & c i + d l & e i + f l\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a⋅g + b⋅j  c⋅g + d⋅j  e⋅g + f⋅j⎤\n",
       "⎢                               ⎥\n",
       "⎢a⋅h + b⋅k  c⋅h + d⋅k  e⋅h + f⋅k⎥\n",
       "⎢                               ⎥\n",
       "⎣a⋅i + b⋅l  c⋅i + d⋅l  e⋅i + f⋅l⎦"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_tr_A_tr = B.transpose()*A.transpose()\n",
    "B_tr_A_tr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see if they are equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\\\0 & 0 & 0\\\\0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0  0  0⎤\n",
       "⎢       ⎥\n",
       "⎢0  0  0⎥\n",
       "⎢       ⎥\n",
       "⎣0  0  0⎦"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AB_tr - B_tr_A_tr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Identity Matrices </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an identity matrix $I_5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0\\\\0 & 1 & 0 & 0 & 0\\\\0 & 0 & 1 & 0 & 0\\\\0 & 0 & 0 & 1 & 0\\\\0 & 0 & 0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  0  0  0  0⎤\n",
       "⎢             ⎥\n",
       "⎢0  1  0  0  0⎥\n",
       "⎢             ⎥\n",
       "⎢0  0  1  0  0⎥\n",
       "⎢             ⎥\n",
       "⎢0  0  0  1  0⎥\n",
       "⎢             ⎥\n",
       "⎣0  0  0  0  1⎦"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy.eye(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identity matrix properties:\n",
    "\n",
    "$$\n",
    "AI=IA = A\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate $ I$ and $ A$ and show if it holds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., 0., 0.],\n",
       "       [0., 1., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0.],\n",
       "       [0., 0., 0., 1., 0.],\n",
       "       [0., 0., 0., 0., 1.]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = np.eye(5); I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[59., 75., 69., 73., 34.],\n",
       "       [94.,  8., 73.,  6., 10.],\n",
       "       [73., 23., 55., 27., 67.],\n",
       "       [ 1., 36., 25., 68., 75.],\n",
       "       [35., 10., 33., 31., 79.]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.around(np.random.rand(5, 5)*100); A # generate a random matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently it holds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@I - I@A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Elementary Matrix </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An **elementary matrix** is a matrix that can be obtained from a single elementary row operation on an identity matrix. Such as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{matrix}1 & 0 & 0\\cr 0 & 1 & 0\\cr 0 & 0 & 1\\end{matrix}\\right]\\ \\matrix{R_1\\leftrightarrow R_2\\cr ~\\cr ~}\\qquad\\Longrightarrow\\qquad \\left[\\begin{matrix}0 & 1 & 0\\cr 1 & 0 & 0\\cr 0 & 0 & 1\\end{matrix}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $R_1\\leftrightarrow R_2$ means exchanging row $1$ and row $2$, and we denote the transformed matrix as ${E}$, then left multiply $  E$ onto a matrix $  A$. Generate $ A$ first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}34 & 5 & 0\\\\34 & 67 & 4\\\\14 & 32 & 44\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡34  5   0 ⎤\n",
       "⎢          ⎥\n",
       "⎢34  67  4 ⎥\n",
       "⎢          ⎥\n",
       "⎣14  32  44⎦"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.randMatrix(3, percent = 80); A # generate a random matrix with 80% of entries being nonzero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an elementary matrix with $R_1\\leftrightarrow R_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 1 & 0\\\\1 & 0 & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0  1  0⎤\n",
       "⎢       ⎥\n",
       "⎢1  0  0⎥\n",
       "⎢       ⎥\n",
       "⎣0  0  1⎦"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = sy.Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]]);E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that by left-multiplying $ E$ onto $ A$, $ A$ also switches the row $1$ and $2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}34 & 67 & 4\\\\34 & 5 & 0\\\\14 & 32 & 44\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡34  67  4 ⎤\n",
       "⎢          ⎥\n",
       "⎢34  5   0 ⎥\n",
       "⎢          ⎥\n",
       "⎣14  32  44⎦"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding a multiple of a row onto another row in the identity matrix also gives us an elementary matrix.\n",
    "\n",
    "$$\n",
    "\\left[\\begin{matrix}1 & 0 & 0\\cr 0 & 1 & 0\\cr 0 & 0 & 1\\end{matrix}\\right]\\ \\matrix{~\\cr ~\\cr R_3-7R_1}\\qquad\\longrightarrow\\left[\\begin{matrix}1 & 0 & 0\\cr 0 & 1 & 0\\cr -7 & 0 & 1\\end{matrix}\\right]\n",
    "$$\n",
    "\n",
    "Let's verify with SymPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}72 & 31 & 34\\\\25 & 12 & 7\\\\54 & 8 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡72  31  34⎤\n",
       "⎢          ⎥\n",
       "⎢25  12  7 ⎥\n",
       "⎢          ⎥\n",
       "⎣54  8   0 ⎦"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0\\\\0 & 1 & 0\\\\-7 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1   0  0⎤\n",
       "⎢        ⎥\n",
       "⎢0   1  0⎥\n",
       "⎢        ⎥\n",
       "⎣-7  0  1⎦"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.randMatrix(3, percent = 80); A\n",
    "E = sy.Matrix([[1, 0, 0], [0, 1, 0], [-7, 0, 1]]); E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}72 & 31 & 34\\\\25 & 12 & 7\\\\-450 & -209 & -238\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 72    31    34 ⎤\n",
       "⎢                ⎥\n",
       "⎢ 25    12    7  ⎥\n",
       "⎢                ⎥\n",
       "⎣-450  -209  -238⎦"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also show this by explicit row operation on $ A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}72 & 31 & 34\\\\25 & 12 & 7\\\\-450 & -209 & -238\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 72    31    34 ⎤\n",
       "⎢                ⎥\n",
       "⎢ 25    12    7  ⎥\n",
       "⎢                ⎥\n",
       "⎣-450  -209  -238⎦"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EA = sy.matrices.MatrixBase.copy(A)\n",
    "EA[2,:]=-7*EA[0,:]+EA[2,:]\n",
    "EA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In next section, we will see an important conclusion of elementary matrices multiplication is that an **invertible matrix** is a product of a series of elementary matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Inverse Matrices </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If ${AB}={BA}=\\mathbf{I}$, $ B$ is called the inverse of matrix $  A$, denoted as $ B=  A^{-1}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy has convenient function ```np.linalg.inv()``` for computing inverse matrices. Generate $ A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.,  -8.,  -6.,  -2., -16.],\n",
       "       [ -1.,  -2.,   8., -13.,   4.],\n",
       "       [  1.,  11., -13.,   2.,   3.],\n",
       "       [ -1., -13.,   6.,   8.,   9.],\n",
       "       [-12.,   9.,  -4.,  -1.,   7.]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.round(10*np.random.randn(5,5)); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.026,  0.013,  0.042, -0.005, -0.078],\n",
       "       [-0.048, -0.038, -0.034, -0.064,  0.01 ],\n",
       "       [-0.05 , -0.034, -0.096, -0.043,  0.003],\n",
       "       [-0.026, -0.08 , -0.04 , -0.002,  0.006],\n",
       "       [-0.015,  0.041,  0.055,  0.049, -0.002]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ainv = np.linalg.inv(A); Ainv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify if they are truly inverse of each other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0., -0.,  0., -0.],\n",
       "       [-0.,  1.,  0., -0.,  0.],\n",
       "       [-0.,  0.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  1., -0.],\n",
       "       [ 0., -0., -0.,  0.,  1.]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@Ainv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ```-0.``` means there are more digits after point, but omitted here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> $[A\\,|\\,I]\\sim [I\\,|\\,A^{-1}]$ Algorithm </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A convenient way of calculating inverse is that we can construct an augmented matrix $[ A\\,|\\,\\mathbf{I}]$, then multiply a series of $ E$'s which are elementary row operations till matrix $A$ becomes row reduced form, if full rank it means $ A \\rightarrow \\mathbf{I}$. Then $\\mathbf{I}$ on the RHS of augmented matrix will be converted into $ A^{-1}$ automatically. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can show with SymPy's ```.rref()``` function on the augmented matrix $[A\\,|\\,I]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1.0 & -8.0 & -6.0 & -2.0 & -16.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\-1.0 & -2.0 & 8.0 & -13.0 & 4.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0\\\\1.0 & 11.0 & -13.0 & 2.0 & 3.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0\\\\-1.0 & -13.0 & 6.0 & 8.0 & 9.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0\\\\-12.0 & 9.0 & -4.0 & -1.0 & 7.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1.0   -8.0   -6.0   -2.0   -16.0  1.0  0.0  0.0  0.0  0.0⎤\n",
       "⎢                                                          ⎥\n",
       "⎢-1.0   -2.0    8.0   -13.0   4.0   0.0  1.0  0.0  0.0  0.0⎥\n",
       "⎢                                                          ⎥\n",
       "⎢ 1.0   11.0   -13.0   2.0    3.0   0.0  0.0  1.0  0.0  0.0⎥\n",
       "⎢                                                          ⎥\n",
       "⎢-1.0   -13.0   6.0    8.0    9.0   0.0  0.0  0.0  1.0  0.0⎥\n",
       "⎢                                                          ⎥\n",
       "⎣-12.0   9.0   -4.0   -1.0    7.0   0.0  0.0  0.0  0.0  1.0⎦"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AI = np.hstack((A, I)) # stack the matrix A and I horizontally\n",
    "AI = sy.Matrix(AI); AI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & -0.0257937689879073 & 0.0134627985941502 & 0.0420563531303985 & -0.00488473223327575 & -0.0783939953535474\\\\0 & 1 & 0 & 0 & 0 & -0.0475017432957576 & -0.0379022282648688 & -0.0338111774166985 & -0.0638116329512683 & 0.00961703559126635\\\\0 & 0 & 1 & 0 & 0 & -0.0496953174550335 & -0.0339723665721724 & -0.0957148213428458 & -0.0433738992707943 & 0.00261056349626286\\\\0 & 0 & 0 & 1 & 0 & -0.0259829910399854 & -0.0804474050298025 & -0.0399608941092372 & -0.00152428875285147 & 0.00566614922611685\\\\0 & 0 & 0 & 0 & 1 & -0.0152533998647413 & 0.04090525231359 & 0.0551652364049212 & 0.0486668605608682 & -0.00159612305039964\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2, \\  3, \\  4\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡1  0  0  0  0  -0.0257937689879073  0.0134627985941502   0.0420563531303985 \n",
       "⎜⎢                                                                            \n",
       "⎜⎢0  1  0  0  0  -0.0475017432957576  -0.0379022282648688  -0.0338111774166985\n",
       "⎜⎢                                                                            \n",
       "⎜⎢0  0  1  0  0  -0.0496953174550335  -0.0339723665721724  -0.0957148213428458\n",
       "⎜⎢                                                                            \n",
       "⎜⎢0  0  0  1  0  -0.0259829910399854  -0.0804474050298025  -0.0399608941092372\n",
       "⎜⎢                                                                            \n",
       "⎝⎣0  0  0  0  1  -0.0152533998647413   0.04090525231359    0.0551652364049212 \n",
       "\n",
       "  -0.00488473223327575  -0.0783939953535474 ⎤                 ⎞\n",
       "                                            ⎥                 ⎟\n",
       "  -0.0638116329512683   0.00961703559126635 ⎥                 ⎟\n",
       "                                            ⎥                 ⎟\n",
       "  -0.0433738992707943   0.00261056349626286 ⎥, (0, 1, 2, 3, 4)⎟\n",
       "                                            ⎥                 ⎟\n",
       "  -0.00152428875285147  0.00566614922611685 ⎥                 ⎟\n",
       "                                            ⎥                 ⎟\n",
       "   0.0486668605608682   -0.00159612305039964⎦                 ⎠"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AI_rref = AI.rref(); AI_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the RHS block, this is the $A^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-0.0257937689879073 & 0.0134627985941502 & 0.0420563531303985 & -0.00488473223327575 & -0.0783939953535474\\\\-0.0475017432957576 & -0.0379022282648688 & -0.0338111774166985 & -0.0638116329512683 & 0.00961703559126635\\\\-0.0496953174550335 & -0.0339723665721724 & -0.0957148213428458 & -0.0433738992707943 & 0.00261056349626286\\\\-0.0259829910399854 & -0.0804474050298025 & -0.0399608941092372 & -0.00152428875285147 & 0.00566614922611685\\\\-0.0152533998647413 & 0.04090525231359 & 0.0551652364049212 & 0.0486668605608682 & -0.00159612305039964\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-0.0257937689879073  0.0134627985941502   0.0420563531303985   -0.00488473223\n",
       "⎢                                                                             \n",
       "⎢-0.0475017432957576  -0.0379022282648688  -0.0338111774166985  -0.06381163295\n",
       "⎢                                                                             \n",
       "⎢-0.0496953174550335  -0.0339723665721724  -0.0957148213428458  -0.04337389927\n",
       "⎢                                                                             \n",
       "⎢-0.0259829910399854  -0.0804474050298025  -0.0399608941092372  -0.00152428875\n",
       "⎢                                                                             \n",
       "⎣-0.0152533998647413   0.04090525231359    0.0551652364049212    0.04866686056\n",
       "\n",
       "327575  -0.0783939953535474 ⎤\n",
       "                            ⎥\n",
       "12683   0.00961703559126635 ⎥\n",
       "                            ⎥\n",
       "07943   0.00261056349626286 ⎥\n",
       "                            ⎥\n",
       "285147  0.00566614922611685 ⎥\n",
       "                            ⎥\n",
       "08682   -0.00159612305039964⎦"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ainv = AI_rref[0][:,5:];Ainv # extract the RHS block"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I wrote a function to round the float numbers to the $4$th digits, on the top of this file, but this is not absolutely necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-0.0258 & 0.0135 & 0.0421 & -0.0049 & -0.0784\\\\-0.0475 & -0.0379 & -0.0338 & -0.0638 & 0.0096\\\\-0.0497 & -0.034 & -0.0957 & -0.0434 & 0.0026\\\\-0.026 & -0.0804 & -0.04 & -0.0015 & 0.0057\\\\-0.0153 & 0.0409 & 0.0552 & 0.0487 & -0.0016\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-0.0258  0.0135   0.0421   -0.0049  -0.0784⎤\n",
       "⎢                                           ⎥\n",
       "⎢-0.0475  -0.0379  -0.0338  -0.0638  0.0096 ⎥\n",
       "⎢                                           ⎥\n",
       "⎢-0.0497  -0.034   -0.0957  -0.0434  0.0026 ⎥\n",
       "⎢                                           ⎥\n",
       "⎢-0.026   -0.0804   -0.04   -0.0015  0.0057 ⎥\n",
       "⎢                                           ⎥\n",
       "⎣-0.0153  0.0409   0.0552   0.0487   -0.0016⎦"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round_expr(Ainv, 4) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify if $AA^{-1}=\\mathbf{I}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.0 & 0 & 0 & 0 & 0.0\\\\0.0 & 1.0 & 0.0 & 0 & 0.0\\\\0.0 & 0.0 & 1.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 1.0 & 0.0\\\\0.0 & 0 & 0 & 0.0 & 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1.0   0    0    0   0.0⎤\n",
       "⎢                       ⎥\n",
       "⎢0.0  1.0  0.0   0   0.0⎥\n",
       "⎢                       ⎥\n",
       "⎢0.0  0.0  1.0  0.0  0.0⎥\n",
       "⎢                       ⎥\n",
       "⎢0.0  0.0  0.0  1.0  0.0⎥\n",
       "⎢                       ⎥\n",
       "⎣0.0   0    0   0.0  1.0⎦"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix(A)\n",
    "round_expr(A*Ainv, 4) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got $\\mathbf{I}$, which means the RHS block is indeed $A^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> An Example of Existence of Inverse </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the values of $\\lambda$ such that the matrix\n",
    "$$A=\\left[ \\begin{matrix}3 &\\lambda &1\\cr 2 & -1 & 6\\cr 1 & 9 & 4\\end{matrix}\\right]$$\n",
    "is not invertible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still,we are using SymPy to solve the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3, lamda, 1],\n",
       "       [2, -1, 6],\n",
       "       [1, 9, 4]], dtype=object)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lamb = sy.symbols('lamda') # SymPy will automatically render into LaTeX greek letters\n",
    "A = np.array([[3, lamb, 1], [2, -1, 6], [1, 9, 4]])\n",
    "I = np.eye(3); A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Form the augmented matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}3 & \\lambda & 1 & 1.0 & 0.0 & 0.0\\\\2 & -1 & 6 & 0.0 & 1.0 & 0.0\\\\1 & 9 & 4 & 0.0 & 0.0 & 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡3  λ   1  1.0  0.0  0.0⎤\n",
       "⎢                       ⎥\n",
       "⎢2  -1  6  0.0  1.0  0.0⎥\n",
       "⎢                       ⎥\n",
       "⎣1  9   4  0.0  0.0  1.0⎦"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AI = np.hstack((A, I))\n",
    "AI = sy.Matrix(AI); AI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & \\frac{116.0}{4.0 \\lambda + 310.0} & \\frac{8.0 \\lambda - 18.0}{4.0 \\lambda + 310.0} & \\frac{- 12.0 \\lambda - 2.0}{4.0 \\lambda + 310.0}\\\\0 & 1 & 0 & \\frac{4.0}{4.0 \\lambda + 310.0} & - \\frac{22.0}{4.0 \\lambda + 310.0} & \\frac{32.0}{4.0 \\lambda + 310.0}\\\\0 & 0 & 1 & \\frac{57.0}{- 6 \\lambda - 465} & \\frac{27.0 - 1.0 \\lambda}{2.0 \\lambda + 155.0} & \\frac{2.0 \\lambda + 3.0}{2.0 \\lambda + 155.0}\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡             116.0       8.0⋅λ - 18.0  -12.0⋅λ - 2.0⎤           ⎞\n",
       "⎜⎢1  0  0  ─────────────  ─────────────  ─────────────⎥           ⎟\n",
       "⎜⎢         4.0⋅λ + 310.0  4.0⋅λ + 310.0  4.0⋅λ + 310.0⎥           ⎟\n",
       "⎜⎢                                                    ⎥           ⎟\n",
       "⎜⎢              4.0           -22.0           32.0    ⎥           ⎟\n",
       "⎜⎢0  1  0  ─────────────  ─────────────  ─────────────⎥, (0, 1, 2)⎟\n",
       "⎜⎢         4.0⋅λ + 310.0  4.0⋅λ + 310.0  4.0⋅λ + 310.0⎥           ⎟\n",
       "⎜⎢                                                    ⎥           ⎟\n",
       "⎜⎢             57.0          27.0 - λ     2.0⋅λ + 3.0 ⎥           ⎟\n",
       "⎜⎢0  0  1   ──────────    ─────────────  ─────────────⎥           ⎟\n",
       "⎝⎣          -6⋅λ - 465    2.0⋅λ + 155.0  2.0⋅λ + 155.0⎦           ⎠"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AI_rref = AI.rref()\n",
    "AI_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the matrix $A$ invertible we notice that is one conditions to be satisfied (in the denominators):\n",
    "\\begin{align}\n",
    "-6\\lambda -465 &\\neq0\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve for $\\lambda$'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFoAAAAVCAYAAADGpvm7AAADF0lEQVR4nO3YXYhVVRQH8J8fmDJBYIE9RGo9lJhkDz34omMDkUTlNX1LJewhEERJKKToTk89SNkHqCAk2GuJSF8ETpPigyBoRAhpjhFmYdKYmfg1Paxz88ydc+eec+fMdUb7w+Zszt7/fdb+733WWntPqFar/sfoY2Kq3omBVDl2Mwwa57jHYA0Hag0TMzr3ohsfjoIhy/EB9uN8YsjHTTh9hhpfK2dK5MAMXMP7uBsvYTeO4x/04wDWyNbuotCuG6fSDZMzOn+D6jDGjASv41FcwC94OCevH1sy3l8omfOcEHA3VmArfkUPfhYLsQw7sCTpM5DiX3RDu07MrDVkCT2a2CAEPo5FYgJ58Kfii98Kp4I/8C0m4Fl8huupPptwCM8L0T/JM3DW9h9N9OBHg3fBWMFdeAJ7hfvYl9Sv1/U7g21JvTPv4O3e0a3iDryA+/E3vhO77lqJnKcxBZ/msOdK8ryaoy/Gj9D3Ylfdu5N4UQTvMjgVsSBfN7FlMlYl9S+b9P0P7XYdreAjdAnhOjAP2zELX4jgOlLOVDyVtF1qYs/beASf46u8kygidJ/GKVNWaZa25UW38Je/iaj+PV7GO5gmO+AV5TyJO0W2MRzW4RVxxlhZZBJFXMcJzVc7jdNFDGkB28SkF5bAqeCyyDAaYS3eww/ibzlX4LuFhO4qMnAb8Hvy7BghZxKeEX9AfwPeerwr/oyu1Di5MR58dCMsSJ4/jZCzUJwCG7mNV4XIR7BYCyIz9oWei+kZ72e6cUVQHwuKcpaJXHlPBucNEfwOi518NpfVGWh3erc0KURGQOyynUn9LDam+q/Aa+KgcxJ/4UGR804VkX9z3TeKcpbioAicaazGWyLv3i8CYT36UrYPi3YLPV9MII0HkkJcxKSF7sFDeEwsSIc4Wh8QOfIuQ0+ZRTiP4z7hGuoxO3lOEj46C73GqNBVxe4fejU+kJTBqSTPLP9cVeLlWpaPftPtcx9dwVHhYspA+j56UbohvaP7RKJfQ8uOfxxhTsnj1e6jh6Be6GrJH77dkL6PHoSxnt7dMvgX2QreoXFYmrMAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ - \\frac{155}{2}\\right]$"
      ],
      "text/plain": [
       "[-155/2]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy.solvers.solve(-6*lamb-465, lamb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this is one of $\\lambda$ that cause the matrix invertible. Let's test with **determinant**. If $| A|=0$, then the matrix is not invertible. Don't worry about determinants, we will get back to this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACEAAAASCAYAAADVCrdsAAABZklEQVR4nM3Vv0pcQRTH8c8N2iToapMqoHEx6wMEE2OlhVj6CqkFLWws1yK1YmEZQt4gbZpUBt8gf2RBCJjKQoKIjVrsCLPjXndmISSn+V3mzvneH2funFO1223/Oh71WXuG9zjFFU6wh8kh+FmskSSpia94ik/4jnlsYhWLOMs0kM1KK3EQkjawhm0sYxctvMs0UMSqon9iBh3dkjVxHQHH8BtVAF8MMFDEiiuxHPRzkgR/cIjHeD3AQDErNtEK+rMGfBz0RYaJIlZsohH0vCbxbn0iw0QRq98VrYsq6E1BThYrNnHnrqF/jCf7HooiVmziR9C6M58NWnfOcRSxYhNfgq64f0xjus3lEkcZJopY8YaO7pWaxnqSuIMn+Ki3RzQxh9FkfxGrSgZY2mq/4RWWdEv3Rm/bPsEUnofnoVhpqTp4iQ8hYSvA9rEgf24UsdIBBr/wNvND0wPeZ7FK+sRfi//CxC1ZQWbEM4QvuwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 0.0$"
      ],
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[3, -155/2, 1], [2, -1, 6], [1, 9, 4]])\n",
    "np.linalg.det(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $| A|$ is $0$. \n",
    "\n",
    "So we found that one condition, as long as $\\lambda \\neq -\\frac{155}{2}$, the matrix $A$ is invertible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Properties of Inverse Matrices </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. If $A$ and $B$ are both invertible, then $(AB)^{-1}=B^{-1}A^{-1}$.\n",
    "2. If $A$ is invertible, then $(A^T)^{-1}=(A^{-1})^T$.\n",
    "3. If $A$ and $B$ are both invertible and symmetric such that $AB=BA$, then $A^{-1}B$ is symmetric."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The _first property_ is straightforward\n",
    "\\begin{align}\n",
    "ABB^{-1}A^{-1}=AIA^{-1}=I=AB(AB)^{-1}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trick of _second property_ is to show that\n",
    "$$\n",
    "A^T(A^{-1})^T = I\n",
    "$$\n",
    "We can use the property of transpose\n",
    "$$\n",
    "A^T(A^{-1})^T=(A^{-1}A)^T = I^T = I\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The _third property_ is to show\n",
    "$$\n",
    "A^{-1}B = (A^{-1}B)^T\n",
    "$$\n",
    "Again use the property of transpose\n",
    "$$\n",
    "(A^{-1}B)^{T}=B^T(A^{-1})^T=B(A^T)^{-1}=BA^{-1}\n",
    "$$\n",
    "We use the $AB = BA$ condition to proceed\n",
    "\\begin{align}\n",
    "AB&=BA\\\\\n",
    "A^{-1}ABA^{-1}&=A^{-1}BAA^{-1}\\\\\n",
    "BA^{-1}&=A^{-1}B\n",
    "\\end{align}\n",
    "The plug in the previous equation, we have\n",
    "$$\n",
    "(A^{-1}B)^{T}=BA^{-1}=A^{-1}B\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}