t14_matrix_inverse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1ouzTiKEWJbR"
   },
   "source": [
    "# Topic 14: Matrix INverse and Pseudo Inverses\n",
    "Author: Evan Chrisney, echrisney@gmail.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "CX1H5eeAWWPy"
   },
   "source": [
    "# Introduction\n",
    " \n",
    "Matrix inverses are matrices which when multiplied by their oringinal matrices equate to the identity matrix. Only square matrices are invertible, but non square matrices may have a left or right psuedo inverse. A matrix $\\mathbf{A}$ has a left inverse if there exists a matrix $\\mathbf{B}$ such that $\\mathbf{BA} = I$, and a right inverse if there exists a matrix $\\mathbf{C}$ such that $\\mathbf{AC} = I$ (Moon and Stirling, 2000). Square matrices that are not invertible are called singular matrices. A matrix is invertible if it has both a right and left inverse. Left and right inverses are types of pseudo inverses. \n",
    "\n",
    "Matrix inverses are important because they provide solutions to the ubiquitous equation $\\mathbf{A}x = b$. There are many conditions that satisfy the invertibility of a matrix which will be explained in detail below. If a true matrix inverse does not exist for $\\mathbf{A}$, or if $\\mathbf{A}$ is not square, then the left and right pseudo inverses are solutions to $\\mathbf{A}x = b$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "b9u2BgA5Wi1I"
   },
   "source": [
    "# Explanation of the theory\n",
    "\n",
    "$\\textbf{Matrix Inverse (nxn):}$\n",
    "\n",
    "A square (n x n) matrix is invertible if there exists a square (n x n) matrix such that\n",
    "$$\\mathbf{A}\\mathbf{B} = \\mathbf{B}\\mathbf{A} = I,$$\n",
    "where $I$ is the identity matrix of size n x n. $\\mathbf{B}$ is a unique matrix hereby notated as $\\mathbf{A^{-1}}$. \n",
    "\n",
    "A matrix $\\mathbf{A}$ is invertible if and only if it has the following properties: \n",
    "\n",
    "1. The null space of $\\mathbf{A}$ is the null vector {0}. \n",
    "2. $\\mathbf{A}$ is full rank, meaning the rank of $\\mathbf{A}$ is n. \n",
    "3. $\\mathbf{A}x = 0$ implies that $x = 0$. \n",
    "4. The rows and columns of $\\mathbf{A}$ are linearly independent.  \n",
    "5. The determinant of $\\mathbf{A}$ is non zero. \n",
    "6. $\\mathbf{A}$ has no zero eigenvalues. \n",
    "7. $\\mathbf{A^HA}$ is positive definite\n",
    "8. $\\mathbf{A}$ is nonsingular. \n",
    "9. $\\mathbf{A}$ has n pivots. \n",
    "10. The transpose of $\\mathbf{A}$, $\\mathbf{A}^T$, is also invertible. \n",
    "11. $\\mathbf{A}$ has both a left inversve and a right inverse $\\mathbf{B}$ and $\\mathbf{C}$ such that $\\mathbf{B}$ = $\\mathbf{C}$ = $\\mathbf{A}^{-1}$. \n",
    "12. $(\\mathbf{A}^{-1})^{-1} = \\mathbf{A}$\n",
    "13. $(\\mathbf{A}^{T})^{-1} = (\\mathbf{A}^{-1})^{T}$\n",
    "\n",
    "The adjugate of $\\mathbf{A}$ can be used to determine $\\mathbf{A}^{-1}$ as\n",
    "$$ \\mathbf{A}^{-1} = \\frac{1}{det(\\mathbf{A})}adj(\\mathbf{A}), $$\n",
    "where $adj$ denotes the adjugate and $det$ denotes the determinant. \n",
    "\n",
    "$\\textbf{Matrix Pseudo Inverses (nxm):}$\n",
    "A non-square (n x m) matrix may have a left or right pseudo inverse, as explained above. \n",
    "\n",
    "$\\textit{Left Inverse} (nxm), m >= n:\\\\$\n",
    "For a system $\\mathbf{A}x = b$, there exists a unique left inverse if and only if the columns of $\\mathbf{A}$ are linearly independent, or $rank(\\mathbf{A} = n)$, and $m >= n$, ie $\\mathbf{A}$ is a tall matrix. \n",
    "\n",
    "The left inverse, or moore penrose pseudo inverse, is \n",
    "$\\mathbf{B} = (\\mathbf{A}^{H}\\mathbf{A})^{-1}\\mathbf{A}^H$, \n",
    "which is the identical to the least squares solution for $\\mathbf{A}x = b$, or \n",
    "$x = (\\mathbf{A}^{H}\\mathbf{A})^{-1}\\mathbf{A}^{H}b$.\n",
    "\n",
    "$\\textit{Right Inverse} (nxm), n >= m:\\\\$\n",
    "For a system $\\mathbf{A}x = b$, there exists at least one solution for any $b$ if and only if the rows of $\\mathbf{A}$ are linearly independent, or $rank(\\mathbf{A} = m)$, and $n >= m$, ie $\\mathbf{A}$ is a fat matrix. In this case, the solution is a right inverse. \n",
    "\n",
    "The right inverse is \n",
    "$\\mathbf{C} = \\mathbf{A}^H(\\mathbf{A}\\mathbf{A}^{H})^{-1}$, \n",
    "which is the identical to the minimum norm solution for $\\mathbf{A}x = b$, or \n",
    "$x = \\mathbf{A}^H(\\mathbf{A}\\mathbf{A}^{H})^{-1}b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zNIUvElqWwdh"
   },
   "source": [
    "# Simple Numerical Example\n",
    "\n",
    "Here is some simple python code that shows the concepts of matrix inversion for a square matrix, as well as computation of the left and right inverses, which all should be equal for this example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 168
    },
    "colab_type": "code",
    "id": "OU8eOXxwW9ET",
    "outputId": "47d57def-d903-45f5-949b-de80af39e9d7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.26666667 -0.6        -0.13333333]\n",
      " [-0.23333333 -0.4        -0.36666667]\n",
      " [-0.13333333  0.2        -0.06666667]]\n",
      "[[-0.26666667 -0.6        -0.13333333]\n",
      " [-0.23333333 -0.4        -0.36666667]\n",
      " [-0.13333333  0.2        -0.06666667]]\n",
      "[[-0.26666667 -0.6        -0.13333333]\n",
      " [-0.23333333 -0.4        -0.36666667]\n",
      " [-0.13333333  0.2        -0.06666667]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# numpy includes its own inverse, but this function drives home the topics above. \n",
    "######################################################\n",
    "# Compute the matrix inverse of a square matrix using the adjugate method\n",
    "\n",
    "\n",
    "# Declare a matrix that we know is invertible\n",
    "A = np.matrix([[-3, 2, -5], [-1, 0, 2], [3, -4, 1]])\n",
    "\n",
    "# get the determinant\n",
    "m = np.linalg.det(A)\n",
    "\n",
    "# Grab length of A which we know is nxn\n",
    "l = len(A) \n",
    "\n",
    "# Initialize the adjugate matrix to 0's\n",
    "C = np.zeros((l,l))\n",
    "\n",
    "# Loop through A and compute the adjugate\n",
    "for i in range(l):\n",
    "        for j in range(l):\n",
    "                # Single out the rows/cols to compute minor determinants\n",
    "                #\n",
    "                # I realize i and j are swapped, this saves the step of\n",
    "                # transposing C later. \n",
    "                temp = np.delete(A, (j), axis = 0)\n",
    "                temp = np.delete(temp, (i), axis = 1)\n",
    "                # Compute minor determinant\n",
    "                M = np.linalg.det(temp)\n",
    "                # Compute the adjugate\n",
    "                C[i][j] = (-1)**(i+j)*M\n",
    "\n",
    "# Finish by computing the adjugate\n",
    "A_inv = 1/m*C\n",
    "# Print out matrix inverse\n",
    "print(A_inv)\n",
    "\n",
    "#####################################################\n",
    "# Now find the left inverse, which should be the same\n",
    "A_tran = A.getH()\n",
    "Grammian_l = np.linalg.inv(A_tran*A)\n",
    "A_dagger_l = Grammian_l*A_tran\n",
    "print(A_dagger_l)\n",
    "\n",
    "\n",
    "#####################################################\n",
    "# Now find the right inverse which is the same\n",
    "Grammian_r = np.linalg.inv(A*A_tran)\n",
    "A_dagger_r = A_tran*Grammian_r\n",
    "print(A_dagger_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wPb-J4i3XNXo"
   },
   "source": [
    "# Engineering Application: Remote Sensing\n",
    "\n",
    "As written above, the matrix inverses provide a solution to the ubiquitous equation $\\mathbf{A}x = b$. An engineering example used in my research quite frequently is using least squares (LS) to fit a line to some data, where we'll solve for $s$ in $\\mathbf{A}s = x$ using a left pseudo inverse. \n",
    "\n",
    "As a brief introduction to this application, I will fit normalized radar cross section (RCS) $\\sigma^0$ data vs incidence angle from a satellite radar at C band know as the advanced scatteromer, or ASCAT. ASCAT is a fan beam scatterometer that observes the earth surface at a variable range of incidence angles, from about 30 to 60 degrees. The purpose of this application is to determine what the $\\sigma^0$ incidence angle dependence ($s$) is at C band in order to normalize ASCAT measurements to one incidence angle for cross calibration purposes.  \n",
    "\n",
    "Previous studies have shown that $\\sigma^0$ exhibits a log-linear dependence with incidence angle over tropical rainforests over the mid incidence angle range from about 30 to 60 degrees incidence (N. Madsen, BYU Masters Thesis, 2015). Due to the log-linear dependence, the dependence is easily estimated as the slope $s$ of the first order polynomial \n",
    "\n",
    "$s_1\\theta + s_2 = \\sigma^0, $\n",
    "\n",
    "where $s_1$ and $s_2$ are the coefficients of $s$ and $s_1$ is the slope which we're solving for, $\\theta$ is the ASCAT incidence angle data, and $\\sigma^0$ is the ASCAT RCS data. In matrix form, this equation is\n",
    "\n",
    "$ [\\mathbf{\\theta}, \\mathbf{1}] [s_1, s_2]^{T} = \\sigma^0$,\n",
    "\n",
    "which easily seen of the form $\\mathbf{A}s = x$. \n",
    "\n",
    "To determine the ASCAT $\\sigma^0$ incidence angle dependence, we create the LS equation\n",
    "\n",
    "$ s = \\mathbf{A}^{\\dagger}\\mathbf{x}, $\n",
    "\n",
    "where $s$ is the dB/degree dependence, $\\mathbf{A}^{\\dagger}$ is the left psuedo inverse of the ASCAT incidence angle data, and $x$ is $\\sigma^0$ RCS data. Again, only the first coefficient of $s$ is of interest, as it gives the slope of the line.\n",
    "\n",
    "In our application, we will only use a small portion of the actual ASCAT data, since there are millions of measurements in a given day. To simplify this, I truncated millions of measurements for a day into 21 measurements to give a feel for what the approximate incidence angle dependence $s$ is. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 50
    },
    "colab_type": "code",
    "id": "FkE9JDD1Xc_B",
    "outputId": "8fcce972-a4bc-4420-81da-29c86ac366e1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.07582702]\n",
      " [-5.07314663]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "##############################################################################################\n",
    "# Evan Chrisney\n",
    "# 671 Application code for Pseudo Inverses\n",
    "# This script will use a left inverse to solve a 1st order poylnomial equation as described above\n",
    "\n",
    "# Create A matrix as described above using actual ASCAT incdience angle values\n",
    "A = np.matrix([[59, 48.63, 50.14, 54.22, 39.23, 52.36, 52.6, 44.27, 55.07, 59.06\\\n",
    ",58.32, 58.5, 37.49, 45.61, 52.13, 42.04, 54.74, 58.53, 38.19, 45.74, 59.31], \\\n",
    "[1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])\n",
    "\n",
    "# Create sigma_0 row vector using actual ASCAT RCS values\n",
    "sigma_0 = np.array([-9.5103, -8.2828, -8.0002, -8.2415, -8.4277, -9.0573, \\\n",
    "-8.6784, -8.5909, -9.0481, -9.6975, -8.8720, -9.7875, -8.3998, -10.0404, \\\n",
    "-10.0735, -6.8197, -9.2742, -9.5769, -7.5610, -8.9871, -10.3787])\n",
    "\n",
    "# Get length of sigma_0 for use in reshape the row vector to a column vec\n",
    "l = len(sigma_0)\n",
    "\n",
    "# Transpose A into a 21x2 matrix instead of 2x21\n",
    "A = A.transpose()\n",
    "# Transpose sigma_0 into a column vector\n",
    "sigma_0 = np.reshape(sigma_0, (l, 1))\n",
    "\n",
    "# Get A hermitian, although transpose would also work since these are reals\n",
    "A_tran = A.getH()\n",
    "# Create the left inverse using the Grammian as above\n",
    "Grammian_l = np.linalg.inv(A_tran*A)\n",
    "# Compute the left inverse of A as above\n",
    "A_dagger_l = Grammian_l*A_tran\n",
    "# Solve for s using the left inverse\n",
    "s = A_dagger_l*sigma_0\n",
    "# Print out s, where the first coefficient is the approximate sigma^0 incidence angle dependence for ASCAT, \n",
    "# which is -0.0758 dB/degree \n",
    "print(s)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Untitled0.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}