notebooks/Chapter 2 - Basic Matrix Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sy\n",
    "sy.init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\" # display multiple lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "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 between two types of multiplication: _Hadamard multiplication_, denoted as $A \\odot B$ (element-wise multiplication), and _matrix multiplication_, denoted simply as $AB$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 2], [3, 4]])\n",
    "B = np.array([[5, 6], [7, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5, 12],\n",
       "       [21, 32]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*B # this is Hadamard elementwise product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[19, 22],\n",
       "       [43, 50]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@B # this is matrix product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's show explicitly the matrix multiplication rule for each element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "43"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "22"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(A[0,:]*B[:,0]) # element at (1, 1)\n",
    "np.sum(A[1,:]*B[:,0]) # element at (2, 1)\n",
    "np.sum(A[0,:]*B[:,1]) # element at (1, 2)\n",
    "np.sum(A[1,:]*B[:,1]) # element at (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 need to use them repeatedly. With this library, we can perform computations analytically, making it a valuable tool for learning linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "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": 61,
   "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": 61,
     "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": 61,
     "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": 62,
   "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": 62,
     "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": 62,
     "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": 63,
   "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": 63,
     "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": 63,
     "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": 64,
   "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": 64,
     "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": [
    "Matrix multiplication usually does not commute, meaning ${AB} \\neq {BA}$. For instance, consider matrices $A$ and $B$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "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": 65,
     "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": 65,
     "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? Let's try find out analytically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "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": 66,
     "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": 66,
     "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": 67,
   "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": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = A*B - B*A; M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That gives us a system of linear equations\n",
    "$$\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}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve it, we can use the Gauss-Jordan elimination method. If we treat $a, b, c, d$ as coefficients of the system, we can extract an augmented matrix.\n",
    "$$\n",
    "\\begin{align*}\n",
    "b g-c f=0 \\quad &\\Rightarrow-c f+b g+0 e+0 h=0\\\\\n",
    "a f-b e+b h-d f=0 \\quad &\\Rightarrow(a-d) f+0 g-b e+b h=0\\\\\n",
    "-a g+c e-c h+d g=0 \\quad &\\Rightarrow 0 f+(d-a) g+c e-c h=0\\\\\n",
    "-b g+c f=0 \\quad &\\Rightarrow c f-b g+0 e+0 h=0\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "So the augmented matrix takes the form\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left[\\begin{array}{cccc:c}\n",
    "-c & b & 0 & 0 & 0 \\\\\n",
    "a-d & 0 & -b & b & 0 \\\\\n",
    "0 & d-a & c & -c & 0 \\\\\n",
    "c & -b & 0 & 0 & 0\n",
    "\\end{array}\\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- c & b & 0 & 0\\\\a - d & 0 & - b & b\\\\0 & - a + d & c & - c\\\\c & - b & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ -c      b     0   0 ⎤\n",
       "⎢                     ⎥\n",
       "⎢a - d    0     -b  b ⎥\n",
       "⎢                     ⎥\n",
       "⎢  0    -a + d  c   -c⎥\n",
       "⎢                     ⎥\n",
       "⎣  c      -b    0   0 ⎦"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_aug = sy.Matrix([[-c, b, 0, 0], [a-d,0, -b, b], [0, d-a, c, -c], [c, -b, 0, 0]]); A_aug"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform _Gaussian-Jordon_ elimination till row reduced formed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & - \\frac{b}{a - d} & \\frac{b}{a - d}\\\\0 & 1 & \\frac{b c}{- a b + b d} & \\frac{c}{a - d}\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right], \\  \\left( 0, \\  1\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡         -b         b  ⎤        ⎞\n",
       "⎜⎢1  0    ─────     ─────⎥        ⎟\n",
       "⎜⎢        a - d     a - d⎥        ⎟\n",
       "⎜⎢                       ⎥        ⎟\n",
       "⎜⎢         b⋅c        c  ⎥        ⎟\n",
       "⎜⎢0  1  ──────────  ─────⎥, (0, 1)⎟\n",
       "⎜⎢      -a⋅b + b⋅d  a - d⎥        ⎟\n",
       "⎜⎢                       ⎥        ⎟\n",
       "⎜⎢0  0      0         0  ⎥        ⎟\n",
       "⎜⎢                       ⎥        ⎟\n",
       "⎝⎣0  0      0         0  ⎦        ⎠"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_aug.rref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general solution is \n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left(\\begin{array}{l}\n",
    "e \\\\\n",
    "f \\\\\n",
    "g \\\\\n",
    "h\n",
    "\\end{array}\\right)=c_1\\left(\\begin{array}{c}\n",
    "\\frac{b}{a-d} \\\\\n",
    "\\frac{b c}{a b-b d} \\\\\n",
    "1 \\\\\n",
    "0\n",
    "\\end{array}\\right)+c_2\\left(\\begin{array}{c}\n",
    "-\\frac{b}{a-d} \\\\\n",
    "-\\frac{c}{a-d} \\\\\n",
    "0 \\\\\n",
    "1\n",
    "\\end{array}\\right)\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution for B elements:\n",
      "e: h + g*(a - d)/c\n",
      "f: b*g/c\n",
      "\n",
      "Parameterized solution in LaTeX format:\n",
      "e: c_{2} + \\frac{c_{1} \\left(a - d\\right)}{c}\n",
      "f: \\frac{b c_{1}}{c}\n",
      "g: c_{1}\n",
      "h: c_{2}\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "\n",
    "# Define symbolic entries for A (2x2 matrix)\n",
    "a, b, c, d = sp.symbols('a b c d')\n",
    "\n",
    "# Define symbolic entries for B (2x2 matrix)\n",
    "e, f, g, h = sp.symbols('e f g h')\n",
    "\n",
    "# Define matrices A and B\n",
    "A = sp.Matrix([[a, b], [c, d]])\n",
    "B = sp.Matrix([[e, f], [g, h]])\n",
    "\n",
    "# Compute AB and BA\n",
    "AB = A * B\n",
    "BA = B * A\n",
    "\n",
    "# Set up equations AB = BA\n",
    "equations = [\n",
    "    sp.Eq(AB[i, j], BA[i, j]) for i in range(2) for j in range(2)\n",
    "]\n",
    "\n",
    "# Solve the system of equations\n",
    "solution = sp.solve(equations, (e, f, g, h))\n",
    "\n",
    "# Print the general solutions\n",
    "print(\"Solution for B elements:\")\n",
    "for sol in solution:\n",
    "    print(f\"{sol}: {solution[sol]}\")\n",
    "\n",
    "# To express the solution in a parameterized form, let's substitute specific values\n",
    "# Define parameters c1 and c2\n",
    "c1, c2 = sp.symbols('c1 c2')\n",
    "\n",
    "# Define the parameterized solution\n",
    "parametric_solution = {\n",
    "    e: solution[e].subs({g: c1, h: c2}),\n",
    "    f: solution[f].subs({g: c1, h: c2}),\n",
    "    g: c1,\n",
    "    h: c2\n",
    "}\n",
    "\n",
    "# Print the parameterized solution in LaTeX format\n",
    "print(\"\\nParameterized solution in LaTeX format:\")\n",
    "for sol in parametric_solution:\n",
    "    print(f\"{sp.latex(sol)}: {sp.latex(parametric_solution[sol])}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "e: c_{2} + \\frac{c_{1} \\left(a - d\\right)}{c}\\\\\n",
    "f: \\frac{b c_{1}}{c}\\\\\n",
    "g: c_{1}\\\\\n",
    "h: c_{2}\\\\\n",
    "$$"
   ]
  },
  {
   "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": 73,
   "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": 73,
     "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": 73,
     "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": 111,
   "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": 111,
     "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_t = AB.transpose(); AB_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transpose $A$ and $B$, then multiply, meaning $B^TA^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "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": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_t_A_t = B.transpose()*A.transpose()\n",
    "B_t_A_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if they are equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AB_t == B_t_A_t"
   ]
  },
  {
   "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$, only $1$'s on principal diagonal, all rest elements are $0$'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "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": 117,
     "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": 118,
   "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": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = np.eye(5); I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[42., 73., 33., 12., 98.],\n",
       "       [28., 33., 18., 45., 19.],\n",
       "       [47., 92., 58., 54., 15.],\n",
       "       [43., 16., 51., 58., 75.],\n",
       "       [27., 71., 16., 17.,  6.]])"
      ]
     },
     "execution_count": 119,
     "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": [
    "Check if they are equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A@I == I@A).all()"
   ]
  },
  {
   "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[\n",
    "\\begin{array}{ccc}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\begin{array}{c}\n",
    "R_1 \\leftrightarrow R_2 \\\\\n",
    "~ \\\\\n",
    "~\n",
    "\\end{array}\n",
    "\\qquad \\Longrightarrow \\qquad\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    "0 & 1 & 0 \\\\\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where $R_1 \\leftrightarrow R_2$ means exchanging row 1 and row 2, we denote the transformed matrix as $E$. Then, left multiply $E$ onto a matrix $A$ will perform the exact the same row operation. \n",
    "\n",
    "First, generate matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 22 & 39\\\\0 & 63 & 77\\\\22 & 72 & 74\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   22  39⎤\n",
       "⎢          ⎥\n",
       "⎢0   63  77⎥\n",
       "⎢          ⎥\n",
       "⎣22  72  74⎦"
      ]
     },
     "execution_count": 81,
     "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": 82,
   "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": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = sy.Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]]);E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that by left-multiplying $E$ onto $A$, matrix $A$ also switches rows 1 and 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 63 & 77\\\\0 & 22 & 39\\\\22 & 72 & 74\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   63  77⎤\n",
       "⎢          ⎥\n",
       "⎢0   22  39⎥\n",
       "⎢          ⎥\n",
       "⎣22  72  74⎦"
      ]
     },
     "execution_count": 83,
     "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",
    "\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{array}\n",
    "\\right]\\\\\n",
    "\n",
    "\\begin{array}{c}\n",
    "R_3-7R_1\n",
    "\\end{array}\n",
    "\\longrightarrow\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "-7 & 0 & 1\n",
    "\\end{array}\n",
    "\\right]\n",
    "\n",
    "$$\n",
    "\n",
    "Let's verify with SymPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 77\\\\90 & 35 & 0\\\\69 & 4 & 16\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   0   77⎤\n",
       "⎢          ⎥\n",
       "⎢90  35  0 ⎥\n",
       "⎢          ⎥\n",
       "⎣69  4   16⎦"
      ]
     },
     "execution_count": 84,
     "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": 84,
     "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": "markdown",
   "metadata": {},
   "source": [
    "We will see the $R_3-7R_1$ takes places on $A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 77\\\\90 & 35 & 0\\\\69 & 4 & -523\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   0    77 ⎤\n",
       "⎢            ⎥\n",
       "⎢90  35   0  ⎥\n",
       "⎢            ⎥\n",
       "⎣69  4   -523⎦"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also reproduce this by explicit row operation on $ A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 77\\\\90 & 35 & 0\\\\69 & 4 & -523\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   0    77 ⎤\n",
       "⎢            ⎥\n",
       "⎢90  35   0  ⎥\n",
       "⎢            ⎥\n",
       "⎣69  4   -523⎦"
      ]
     },
     "execution_count": 86,
     "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 the next section, we will refresh an important conclusion: 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": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-10., -14.,  16.,   3.,  -7.],\n",
       "       [  1.,  -3.,   9.,   5.,   0.],\n",
       "       [-15.,   2.,  -8.,  17.,  15.],\n",
       "       [ 23.,  -7.,  -6., -21.,  -7.],\n",
       "       [ 18.,   2.,  18.,  -6., -12.]])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.round(10*np.random.randn(5,5)); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.131, -0.442,  0.357,  0.116,  0.302],\n",
       "       [-0.157,  0.346, -0.286, -0.115, -0.198],\n",
       "       [-0.336,  1.136, -0.736, -0.2  , -0.608],\n",
       "       [ 0.483, -1.549,  1.082,  0.267,  0.915],\n",
       "       [-0.575,  1.873, -1.158, -0.278, -1.032]])"
      ]
     },
     "execution_count": 88,
     "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": 89,
   "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": 89,
     "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\"> Gauss-Jordan Elimination Method for Matrix Inversion</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A convenient way to calculate the inverse of a matrix is to construct an augmented matrix $[A \\,|\\, \\mathbf{I}]$. Then, multiply a series of elementary matrices $E$ (representing elementary row operations) until matrix $A$ is in row-reduced form. If $A$ is of full rank, this process will transform $A$ into an identity matrix $\\mathbf{I}$. Consequently, the identity matrix on the right-hand side of the augmented matrix will be converted into $A^{-1}$ automatically.\n",
    "\n",
    "We can demonstrate this using SymPy's `.rref()` function on the augmented matrix $[A \\,|\\, \\mathbf{I}]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-10.0 & -14.0 & 16.0 & 3.0 & -7.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\1.0 & -3.0 & 9.0 & 5.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0\\\\-15.0 & 2.0 & -8.0 & 17.0 & 15.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0\\\\23.0 & -7.0 & -6.0 & -21.0 & -7.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0\\\\18.0 & 2.0 & 18.0 & -6.0 & -12.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-10.0  -14.0  16.0   3.0   -7.0   1.0  0.0  0.0  0.0  0.0⎤\n",
       "⎢                                                         ⎥\n",
       "⎢ 1.0   -3.0   9.0    5.0    0.0   0.0  1.0  0.0  0.0  0.0⎥\n",
       "⎢                                                         ⎥\n",
       "⎢-15.0   2.0   -8.0  17.0   15.0   0.0  0.0  1.0  0.0  0.0⎥\n",
       "⎢                                                         ⎥\n",
       "⎢23.0   -7.0   -6.0  -21.0  -7.0   0.0  0.0  0.0  1.0  0.0⎥\n",
       "⎢                                                         ⎥\n",
       "⎣18.0    2.0   18.0  -6.0   -12.0  0.0  0.0  0.0  0.0  1.0⎦"
      ]
     },
     "execution_count": 90,
     "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": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & 0.130965315503552 & -0.442457166736314 & 0.356957793564563 & 0.116005014625992 & 0.302131216046803\\\\0 & 1 & 0 & 0 & 0 & -0.157459256163811 & 0.34625992478061 & -0.285582950271626 & -0.114584203928124 & -0.19828666945257\\\\0 & 0 & 1 & 0 & 0 & -0.335645633096532 & 1.13618888424572 & -0.73614709569578 & -0.199665691600502 & -0.607918930213122\\\\0 & 0 & 0 & 1 & 0 & 0.48349352277476 & -1.54889260342666 & 1.08232344337652 & 0.26744671959883 & 0.914855829502717\\\\0 & 0 & 0 & 0 & 1 & -0.575010447137485 & 1.87275386544087 & -1.15754283326369 & -0.278311742582533 & -1.03249059757626\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2, \\  3, \\  4\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡1  0  0  0  0  0.130965315503552   -0.442457166736314  0.356957793564563   0 ↪\n",
       "⎜⎢                                                                             ↪\n",
       "⎜⎢0  1  0  0  0  -0.157459256163811   0.34625992478061   -0.285582950271626  - ↪\n",
       "⎜⎢                                                                             ↪\n",
       "⎜⎢0  0  1  0  0  -0.335645633096532   1.13618888424572   -0.73614709569578   - ↪\n",
       "⎜⎢                                                                             ↪\n",
       "⎜⎢0  0  0  1  0   0.48349352277476   -1.54889260342666    1.08232344337652     ↪\n",
       "⎜⎢                                                                             ↪\n",
       "⎝⎣0  0  0  0  1  -0.575010447137485   1.87275386544087   -1.15754283326369   - ↪\n",
       "\n",
       "↪ .116005014625992   0.302131216046803 ⎤                 ⎞\n",
       "↪                                      ⎥                 ⎟\n",
       "↪ 0.114584203928124  -0.19828666945257 ⎥                 ⎟\n",
       "↪                                      ⎥                 ⎟\n",
       "↪ 0.199665691600502  -0.607918930213122⎥, (0, 1, 2, 3, 4)⎟\n",
       "↪                                      ⎥                 ⎟\n",
       "↪ 0.26744671959883   0.914855829502717 ⎥                 ⎟\n",
       "↪                                      ⎥                 ⎟\n",
       "↪ 0.278311742582533  -1.03249059757626 ⎦                 ⎠"
      ]
     },
     "execution_count": 91,
     "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": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.130965315503552 & -0.442457166736314 & 0.356957793564563 & 0.116005014625992 & 0.302131216046803\\\\-0.157459256163811 & 0.34625992478061 & -0.285582950271626 & -0.114584203928124 & -0.19828666945257\\\\-0.335645633096532 & 1.13618888424572 & -0.73614709569578 & -0.199665691600502 & -0.607918930213122\\\\0.48349352277476 & -1.54889260342666 & 1.08232344337652 & 0.26744671959883 & 0.914855829502717\\\\-0.575010447137485 & 1.87275386544087 & -1.15754283326369 & -0.278311742582533 & -1.03249059757626\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0.130965315503552   -0.442457166736314  0.356957793564563   0.116005014625992 ↪\n",
       "⎢                                                                              ↪\n",
       "⎢-0.157459256163811   0.34625992478061   -0.285582950271626  -0.11458420392812 ↪\n",
       "⎢                                                                              ↪\n",
       "⎢-0.335645633096532   1.13618888424572   -0.73614709569578   -0.19966569160050 ↪\n",
       "⎢                                                                              ↪\n",
       "⎢ 0.48349352277476   -1.54889260342666    1.08232344337652    0.26744671959883 ↪\n",
       "⎢                                                                              ↪\n",
       "⎣-0.575010447137485   1.87275386544087   -1.15754283326369   -0.27831174258253 ↪\n",
       "\n",
       "↪    0.302131216046803 ⎤\n",
       "↪                      ⎥\n",
       "↪ 4  -0.19828666945257 ⎥\n",
       "↪                      ⎥\n",
       "↪ 2  -0.607918930213122⎥\n",
       "↪                      ⎥\n",
       "↪    0.914855829502717 ⎥\n",
       "↪                      ⎥\n",
       "↪ 3  -1.03249059757626 ⎦"
      ]
     },
     "execution_count": 92,
     "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, just for sake of readability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.131 & -0.4425 & 0.357 & 0.116 & 0.3021\\\\-0.1575 & 0.3463 & -0.2856 & -0.1146 & -0.1983\\\\-0.3356 & 1.1362 & -0.7361 & -0.1997 & -0.6079\\\\0.4835 & -1.5489 & 1.0823 & 0.2674 & 0.9149\\\\-0.575 & 1.8728 & -1.1575 & -0.2783 & -1.0325\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 0.131   -0.4425   0.357    0.116   0.3021 ⎤\n",
       "⎢                                           ⎥\n",
       "⎢-0.1575  0.3463   -0.2856  -0.1146  -0.1983⎥\n",
       "⎢                                           ⎥\n",
       "⎢-0.3356  1.1362   -0.7361  -0.1997  -0.6079⎥\n",
       "⎢                                           ⎥\n",
       "⎢0.4835   -1.5489  1.0823   0.2674   0.9149 ⎥\n",
       "⎢                                           ⎥\n",
       "⎣-0.575   1.8728   -1.1575  -0.2783  -1.0325⎦"
      ]
     },
     "execution_count": 93,
     "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": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 1.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 & 0.0 & 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1.0  0.0  0.0  0.0  0.0⎤\n",
       "⎢                       ⎥\n",
       "⎢0.0  1.0   0   0.0  0.0⎥\n",
       "⎢                       ⎥\n",
       "⎢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  0.0  1.0⎦"
      ]
     },
     "execution_count": 94,
     "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": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3, lamda, 1],\n",
       "       [2, -1, 6],\n",
       "       [1, 9, 4]], dtype=object)"
      ]
     },
     "execution_count": 123,
     "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": 96,
   "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": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AI = np.hstack((A, I))\n",
    "AI = sy.Matrix(AI); AI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "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 - 1.0⋅λ    2.0⋅λ + 3.0 ⎥           ⎟\n",
       "⎜⎢0  0  1   ──────────    ─────────────  ─────────────⎥           ⎟\n",
       "⎝⎣          -6⋅λ - 465    2.0⋅λ + 155.0  2.0⋅λ + 155.0⎦           ⎠"
      ]
     },
     "execution_count": 97,
     "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 multiple conditions to be satisfied, (in the denominators):\n",
    "$$\n",
    "\\begin{align*}\n",
    "-6\\lambda -465 &\\neq0\\\\\n",
    "4 \\lambda + 310 &\\neq 0\\\\\n",
    "2 \\lambda + 155 &\\neq 0\n",
    "\\end{align*}\n",
    "$$\n",
    "However, they are actually one condition, because they are multiples of each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve for $\\lambda$'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFoAAAAUCAYAAAAN+ioeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAABJ0AAASdAHeZh94AAADU0lEQVR4nO3YXYhVVRgG4GdSSrC6KVCIyuqihEKDLopArSEo+nMs6cayyy5CDIMgim++IOgiyCDKIJgBq7uKLtTqwnGmIAiC/oiiH6cIy7JQmiQqnS72PuPxeM6Z8zenGeyFw9577fWu8653r/3tb30D09PT/sfc44z/WsDpgsWVk8xch7Gqe19GxBX9FrSQkZnn45fqtogYoMroKoxjHw7NgZC7sRarsQrn4JWI2NSEM4mLG9w+GBHLe8EpectwAM9jGEO4FVfhAvyFTzGCkYg4XjPEUWR5fn+1hnpG74uI4QYiu8VjCoOn8ANafWOOYHud9qkec+5UhNPXsREv4EfFm/49lmEDXsItmbkxImY+chFxVPGAKhGiqdFziYcUBn+tWNljzbvP4HAHD78TzhB+xQQGcAd2Va/czHwUH+AuhemvtTJwX42OiBljM7NZ174jM8/FjXg1Io5hb71+EfFTZu7Ak1hnPhrdBc7KzE24CH/gE0yUhvSKcxvOxBst6Pm7PP7TingWTnq3HDsVq2i7YrV9lZlre8gZUjyQd5oJyczFuK+8fKs1+QvD6BEMKoxbqsgAXsQK7MnMVd1yMnMJbsaeiPhzFj1P4Ursjoi3W51Ey6FjlpSpHpqmba0iImqD+Wd4IDOnsM2JNKwbzk042yxhIzO3lPwvcG8782gnRn+D2Z52NQ60I6QD7FBMek0POBsUOfKuRsTMfBDP4nMMRsRv7Yht2eiIGGxn4D6gsgNb2g0nMxfhduyNiCP1SJm5Fc8o3ozBiPi5XbELIUY3wrXl8dsuOWtwngZhIzMfUZj8EW7oxGTmudGZuTIzT1mxmbkCz5WXL3fJGcJxvFmH87ji4/ehYiV3XJboax6dmeuxvrys1Buuy8zR8vxQRDxcRbkH2zJzAt/hd1ymqD8swW48XfM3LXMyc6DU835EHKzRuhlP4BjexZY6m6zJiBitbayHfm9YVmNzTdul5Y/CmGqjx3A5rsb1ith6GO8pcuSd1bWGDjjX4EL1ayKXlMdF2NpgPuMYbXDvJPR7Cz6sLLq02H9cMZl2/qMdTiXFOyU+t6t1NtQzOjIznB716CF8HBH7ezFYvXp0BdVGTzpRS2UO6tHzDRGxssdDVtejT8K/1Lo4zMoWgV8AAAAASUVORK5CYII=",
      "text/latex": [
       "$\\displaystyle \\left[ - \\frac{155}{2}\\right]$"
      ],
      "text/plain": [
       "[-155/2]"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy.solvers.solve(-6*lamb-465, lamb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this is the $\\lambda$ that makes the matrix invertible. Let's test this with the _determinant_. If $|A| = 0$, then the matrix is not invertible, so plug in $\\lambda$ back in $A$ then calculate determinant. Don't worry about determinants; we will revisit this topic later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACEAAAAQCAYAAACYwhZnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAABJ0AAASdAHeZh94AAABZklEQVR4nM3VsWoVQRSH8d+9aGMIEQQrQeUS9QFEA6mCEHwJ+6AQCJYWhyPYCipaC76AtgGxUnwDTSRgl0pUJJZeizuRuYtLdizUaf4we863386wM6PpdOpfj2Pdicw8g7u4jlPYx3NkRHxugQ9ljeqVyMwJ3uA0XuA9rmANO1iNiE8DBQazuivxpDRtRsSjCngfW7iHjSESLaxxx3wdH/G4Awwc4EZmLhz19lbWuHq4VnI7In7MdUV8w2ucwMpREq2sWuJiyd0e8IeSFwZINLFqiaWSX3saD+dPDpBoYo17iv7qqCUO7ZZ+V1jNfxnAbWLVEjsl+/Z8uWTfPtejiVVLvCq5nplz25SZi1jFd7wdINHE+lUQEXvYxjnc6kATC3gWEQcVcJKZlzLzeF3cyuqemDfNjtqHmXkN73DV7L/fxZ1O/UucxXmzg+mPWHNLVb7gMp6WhtuY4AFWht4brazR/3CV/wTkupu37oWhFwAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle 0.0$"
      ],
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 99,
     "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",
    "$$\n",
    "\\begin{align}\n",
    "ABB^{-1}A^{-1}=AIA^{-1}=I=AB(AB)^{-1}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "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.12.0"
  },
  "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
}