12-Matrix-Inverse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d9d98a73",
   "metadata": {},
   "source": [
    "# 12 - Matrix Inverse\n",
    "- 12.1 Concepts and applications\n",
    "- 12.2 Inverse of a diagonal matrix\n",
    "- 12.3 Inverse of a $2 \\times 2$ matrix\n",
    "- 12.4 MCA Algorithm\n",
    "- 12.5 Inverse via row reduction\n",
    "- 12.6 Left inverse for tall matrices\n",
    "- 12.7 Right inverse for wide matrices\n",
    "- 12.8 Pseudoinverse, part 1\n",
    "- 12.9 Exercises\n",
    "- 12.10 Answers\n",
    "- 12.11 Coding challenges\n",
    "- 12.12 Code solutions\n",
    "\n",
    "Notes, code snippets, and the end of chapter exercises from the book _Linear Algebra: Theory, Intuition, Code_ by Mike X Cohen. \n",
    "\n",
    "Find more information about the book on [github](https://github.com/mikexcohen/LinAlgBook) and [amazon](https://www.amazon.com/Linear-Algebra-Theory-Intuition-Code/dp/9083136604)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "807516f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sym"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88e6bce2",
   "metadata": {},
   "source": [
    "## 12.1 Concepts and applications\n",
    "Inverse is the matrix analog of the multiplicative inverse (division).\n",
    "\n",
    "$$\n",
    "A A^{-1} = A^{-1} A = I\n",
    "$$\n",
    "\n",
    "Inverse of the inverse is the original matrix.\n",
    "\n",
    "$$\n",
    "(A^{-1})^{-1} = A\n",
    "$$\n",
    "\n",
    "Transpose of the inverse is equal to the inverse of the transpose. Notation for this is $A^{-T}$.\n",
    "\n",
    "$$\n",
    "(A^{-1})^{T} = (A^{T})^{-1} = A^{-T}\n",
    "$$\n",
    "\n",
    "Note\n",
    "- Matrix inverse exists if and only if matrix is square and full rank.\n",
    "- Inverse of a square matrix, if it exists, is the same whether the inverse pre-multiplies the matrix on the left or post-multiplies the matrix on the right.\n",
    "- Inverse of a matrix is unique.\n",
    "- Similar to determinants (_section 11.9 Theory vs. Practice_) explicitly computing the inverse is not computationally useful.\n",
    "    - Most problems that require the inverse in theory are solved without explicitly computing the inverse.\n",
    "\n",
    "### Inverse of Symmetric Matrices\n",
    "Inverse of a symmetric matrix is itself symmetric. If $A = A^T$ then $A^{-1} = A^{-T}$.\n",
    "\n",
    "### Methods for Computing Inverse\n",
    "1. MCA aka Cramer's Rule (_section 12.4_)\n",
    "2. RREF, row-reduced echelon form reduction (_section 12.5_)\n",
    "3. SVD, singular value decomposition (_chapter 16_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dbafcd8",
   "metadata": {},
   "source": [
    "## 12.2 Inverse of a diagonal matrix\n",
    "Inverse of a diagonal matrix is equal to the reciprocal of the diagonal elements.\n",
    "\n",
    "Example\n",
    "$$\n",
    "A =\n",
    "\\begin{bmatrix}\n",
    "d_1 & 0 & 0 \\\\\n",
    "0 & d_2 & 0 \\\\\n",
    "0 & 0 & d_3 \\\\\n",
    "\\end{bmatrix} \\quad\n",
    "A^{-1} = \n",
    "\\begin{bmatrix}\n",
    "\\frac{1}{d_1} & 0 & 0 \\\\\n",
    "0 & \\frac{1}{d_2} & 0 \\\\\n",
    "0 & 0 & \\frac{1}{d_3} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Note\n",
    "- If any of the diagonal elements are 0, then the matrix must be singular and inverse does not exist."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57d86466",
   "metadata": {},
   "source": [
    "## 12.3 Inverse of a $2 \\times 2$ matrix\n",
    "Inverse of a $2 \\times 2$ matrix is given as a special case of Cramer's Rule.\n",
    "\n",
    "$$\n",
    "A = \n",
    "\\begin{bmatrix}\n",
    "a & b \\\\\n",
    "c & d \\\\\n",
    "\\end{bmatrix} \\quad\n",
    "A^{-1} = \n",
    "\\frac{1}{ad - bc}\n",
    "\\begin{bmatrix}\n",
    "d & -b \\\\\n",
    "-c & a \\\\\n",
    "\\end{bmatrix} \\quad\n",
    "$$\n",
    "\n",
    "Note\n",
    "- If the $\\text{det}(A) = 0$, then the matrix must be singular and inverse does not exist."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8a7b4f6",
   "metadata": {},
   "source": [
    "## 12.4 MCA Algorithm\n",
    "MCA algorithm stands for minors, cofactors, and adjoint is used to find the inverse of a matrix and is also known as _Cramer's Rule_.\n",
    "\n",
    "The adjoint matrix of $A$ is the inverse $A^{-1}$.\n",
    "\n",
    "1. Compute the minors matrix $M$.\n",
    "- Element $m_{i,j}$ of $M$ equals the determinant of the submatrix formed by excluding the _ith_ row and _jth_ column of $A$.\n",
    "\n",
    "2. Compute the cofactors matrix $C$.\n",
    "- Cofactors matrix $C$ equal to the Hadamard product of chessboard matrix $G$ and minors matrix $M$.\n",
    "$$\n",
    "C = G \\odot M\n",
    "$$\n",
    "- Chessboard matrix $G$ is a matrix of alternating $\\pm 1$ with $g_{i,j} = (-1)^{i+j}$ with upper left hand corner starting from $+1$.\n",
    "\n",
    "3. Compute the adjoint matrix $A^{-1}$.\n",
    "- Adjoint matrix $A^{-1}$ is equal to the transpose of the cofactors matrix multiplied by the reciprocal of the determinant of the original matrix.\n",
    "$$\n",
    "A^{-1} = \\frac{1}{\\text{det}(A)} C^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06d36929",
   "metadata": {},
   "source": [
    "## 12.5 Inverse via row reduction\n",
    "\n",
    "### RREF Inverse\n",
    "Inverse of a matrix $A$ can be obtained via RREF in a manner analogous to Gauss-Jordan elimination.\n",
    "\n",
    "1. Augment $A$ with $I$ to produce the matrix $[ A | I ]$\n",
    "2. Use elementary row operations to transform augmented matrix $\\text{RREF}([ A | I ]) = [ I | A^{-1} ]$\n",
    "\n",
    "Notes\n",
    "- By comparison Gauss-Jordan is used to solve $Ax = b$ for $x$ and the approach to find the inverse solves $A A^{-1} = I$ for $A^{-1}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2373a2b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5\n",
    "A = np.random.random((n, n))\n",
    "I = np.eye(n)\n",
    "\n",
    "# Augment A and I to produce a single matrix [A|I].\n",
    "AI = np.hstack((A, I))\n",
    "\n",
    "# Compute RREF(A|I) = [I|A^{-1}].\n",
    "rrefAI, _ = sym.Matrix(AI).rref()\n",
    "Ainv = np.array(rrefAI[:,n:])\n",
    "\n",
    "# Verify A A^{-1} = I.\n",
    "np.testing.assert_almost_equal(A @ Ainv, np.eye(n))\n",
    "\n",
    "# Verify the inverse computed via RREF equals the library for returning the inverse.\n",
    "np.testing.assert_almost_equal(Ainv, np.linalg.inv(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f6c307f",
   "metadata": {},
   "source": [
    "## 12.6 Left inverse for tall matrices\n",
    "Tall matrices $T$ do not have a two-sided inverse, but they have a one-sided left inverse if $\\text{rank}(T) = n$.\n",
    "\n",
    "Left inverse of a tall matrix $T$ is given by the product of 3 matrices.\n",
    "\n",
    "$$\n",
    "(T^T T)^{-1} T^T T = I\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Tall matrices are rectangular, non-square, and have more rows than columns.\n",
    "- Product $T^T T$ is square $n \\times n$ and invertible if $\\text{rank}(T) = n$ where $n$ is the number of columns.\n",
    "- Left inverse is used in overdetermined systems eg systems with more equations than unknowns aka _least squares_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4d259aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify that the product of the left inverse of a \n",
    "# tall full-column-rank matrix and itself is the identity matrix.\n",
    "m, n = 7, 5\n",
    "T = np.random.random((m, n))\n",
    "\n",
    "Tleftinv = np.linalg.inv(T.T @ T) @ T.T\n",
    "\n",
    "np.testing.assert_almost_equal(Tleftinv @ T, np.eye(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e93820bc",
   "metadata": {},
   "source": [
    "## 12.7 Right inverse for wide matrices\n",
    "Wide matrices $W$ do not have a two-sided inverse, but they can have a one-sided right inverse if $\\text{rank}(W) = m$.\n",
    "\n",
    "Right inverse of a wide matrix $W$ is given by the product of 3 matrices.\n",
    "\n",
    "$$\n",
    "W W^T (W W^T)^{-1} = I\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Wide matrices are rectangular, non-square, and have more columns than rows.\n",
    "- Product $W W^T$ is square $m \\times m$ and invertible if $\\text{rank}(W) = m$ where $m$ is the number of rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f64c41b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify that the product of the right inverse of a \n",
    "# wide full-row-rank matrix and itself is the identity matrix.\n",
    "m, n = 5, 7\n",
    "W = np.random.random((m, n))\n",
    "\n",
    "Wrightinv = W.T @ np.linalg.inv(W @ W.T)\n",
    "\n",
    "np.testing.assert_almost_equal(W @ Wrightinv, np.eye(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e61c090",
   "metadata": {},
   "source": [
    "## 12.8 Pseudoinverse, part 1\n",
    "\n",
    "Notes\n",
    "- Pseudoinverse exists even when the matrix is singular.\n",
    "- Pseudoinverse exists even when the matrix is not square.\n",
    "- Pseudoinverse is the same as the inverse for full-rank square matrices.\n",
    "- Notation used for pseudoinverse varies $A^\\dagger, A^*, A^+$ are in common use.\n",
    "- For tall full-column-rank matrices $A^+ = (T^T T)^{-1} T^T$ (_left-sided inverse_).\n",
    "- For wide full-row-rank matrices $A^+ = W^T (W W^T)^{-1}$ (_right-sided inverse_)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8838115",
   "metadata": {},
   "source": [
    "## 12.11 Coding challenges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3645b4a5",
   "metadata": {},
   "source": [
    "> Implement the MCA algorithm to compute the matrix inverse. Test your algorithm by computing the inverse of random matrices and compare against the library function for returning the inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "91c5f6a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def minors(A):\n",
    "    \"\"\"\n",
    "    minors returns the minors matrix of A\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :return: numpy.ndarray   Matrix M\n",
    "    \"\"\"\n",
    "    assert A.shape[0] == A.shape[1], \"A is not square\"\n",
    "\n",
    "    M, n = np.zeros_like(A), A.shape[0]\n",
    "    for i in range(n):\n",
    "        idx_i = np.r_[:i:,i+1:n:]\n",
    "        for j in range(n):\n",
    "            idx_j = np.r_[:j:,j+1:n:]\n",
    "            # Compute the determinant of the submatrix formed by \n",
    "            # excluding row i and column j of the matrix A.\n",
    "            M[i,j] = np.linalg.det(A[np.ix_(idx_i, idx_j)])\n",
    "    return M\n",
    "\n",
    "\n",
    "def chessboard(m, n):\n",
    "    \"\"\"\n",
    "    chessboard returns the matrix G of +/- 1 values according to $g_{i,j} = (-1)^{i+j}$\n",
    "\n",
    "    :param m: int            Number of rows.\n",
    "    :param n: int            Number of columns.\n",
    "    :return: numpy.ndarray   Matrix G\n",
    "    \"\"\"\n",
    "    G = np.zeros((m, n))\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            G[i,j] = (-1)**(i+j)\n",
    "    return G\n",
    "\n",
    "\n",
    "def adjoint(A):\n",
    "    \"\"\"\n",
    "    adjoint returns the inverse of A using MCA algorithm aka Cramer's Rule\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :return: numpy.ndarray   Matrix Ainv\n",
    "    \"\"\"\n",
    "    assert A.shape[0] == A.shape[1], \"A is not square\"\n",
    "\n",
    "    # Compute the determinant of A.\n",
    "    detA = np.linalg.det(A)\n",
    "    assert detA != 0., \"A is singular\"\n",
    "\n",
    "    # Compute the cofactors matrix.\n",
    "    C = chessboard(A.shape[0], A.shape[0]) * minors(A)\n",
    "\n",
    "    return 1./detA * C.T\n",
    "\n",
    "\n",
    "n = 5\n",
    "A = np.random.random((n, n))\n",
    "Ainv = adjoint(A)\n",
    "\n",
    "# Verify A A^{-1} = I.\n",
    "np.testing.assert_almost_equal(A @ Ainv, np.eye(n))\n",
    "\n",
    "# Verify the adjoint equals the inverse.\n",
    "np.testing.assert_almost_equal(Ainv, np.linalg.inv(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d883f9",
   "metadata": {},
   "source": [
    "> Show that (1) the pseudoinverse is the same as the real inverse for a full-rank square matrix; and (2) the pseudoinverse is the same as the left inverse for a full-column-rank tall matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "80c10b19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify the p-inv is the same as the real inverse for a full-rank square matrix.\n",
    "n = 5\n",
    "A = np.random.random((n, n))\n",
    "\n",
    "Ainv = np.linalg.inv(A)\n",
    "Apinv = np.linalg.pinv(A)\n",
    "\n",
    "np.testing.assert_almost_equal(Ainv, Apinv)\n",
    "\n",
    "# Verify the p-inv is the same as the left inverse for a full-column-rank tall matrix.\n",
    "m, n = 7, 5\n",
    "T = np.random.random((m, n))\n",
    "\n",
    "Tleftinv = np.linalg.inv(T.T @ T) @ T.T\n",
    "Tpinv = np.linalg.pinv(T)\n",
    "\n",
    "np.testing.assert_almost_equal(Tleftinv, Tpinv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dc97a3a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}