10_LA_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Note:  This material largely follows the text \"Numerical Linear Algebra\" by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "init_cell": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%precision 3\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Linear Algebra\n",
    "\n",
    "Numerical methods for linear algebra problems lies at the heart of many numerical approaches and is something we will spend some time on.  Roughly we can break down problems that we would like to solve into two general problems, solving a system of equations\n",
    "\n",
    "$$A \\vec{x} = \\vec{b}$$\n",
    "\n",
    "and solving the eigenvalue problem\n",
    "\n",
    "$$A \\vec{v} = \\lambda \\vec{v}.$$\n",
    "\n",
    "We examine each of these problems separately and will evaluate some of the fundamental properties and methods for solving these problems. We will be careful in deciding how to evaluate the results of our calculations and try to gain some understanding of when and how they fail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## General Problem Specification\n",
    "\n",
    "The number and power of the different tools made available from the study of linear algebra makes it an invaluable field of study. Before we dive in to numerical approximations we first consider some of the pivotal problems that numerical methods for linear algebra are used to address.\n",
    "\n",
    "For this discussion we will be using the common notation $m \\times n$ to denote the dimensions of a matrix $A$.  The $m$ refers to the number of rows and $n$ the number of columns.  If a matrix is square, i.e. $m = n$, then we will use the notation that $A$ is $m \\times m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Systems of Equations\n",
    "\n",
    "The first type of problem is to find the solution to a linear system of equations.  If we have $m$ equations for $m$ unknowns it can be written in matrix/vector form,\n",
    "\n",
    "$$A \\vec{x} = \\vec{b}.$$\n",
    "\n",
    "For this example $A$ is an $m \\times m$ matrix, denoted as being in $\\mathbb{R}^{m\\times m}$, and $\\vec{x}$ and $\\vec{b}$ are column vectors with $m$ entries, denoted as $\\mathbb{R}^m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Vandermonde Matrix\n",
    "\n",
    "We have data $(x_i, y_i), ~~ i = 1, 2, \\ldots, m$ that we want to fit a polynomial of order $m-1$.  Solving the linear system $A p = y$ does this for us where\n",
    "\n",
    "$$A = \\begin{bmatrix}\n",
    "    1 & x_1 & x_1^2 & \\cdots & x_1^{m-1} \\\\\n",
    "    1 & x_2 & x_2^2 & \\cdots & x_2^{m-1} \\\\\n",
    "    \\vdots & \\vdots & \\vdots & & \\vdots \\\\\n",
    "    1 & x_m & x_m^2 & \\cdots & x_m^{m-1}\n",
    "\\end{bmatrix} \\quad \\quad y = \\begin{bmatrix}\n",
    "y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_m\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "and $p$ are the coefficients of the interpolating polynomial $\\mathcal{P}_N(x) = p_0 + p_1 x + p_2 x^2 + \\cdots + p_m x^{m-1}$. The solution to this system satisfies $\\mathcal{P}_N(x_i)=y_i$ for $i=1, 2, \\ldots, m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Linear least squares 1\n",
    "\n",
    "In a similar case as above, say we want to fit a particular function (could be a polynomial) to a given number of data points except in this case we have more data points than free parameters.  In the case of polynomials this could be the same as saying we have $m$ data points but only want to fit a $n - 1$ order polynomial through the data where $n - 1 \\leq m$.  One of the common approaches to this problem is to minimize the \"least-squares\" error between the data and the resulting function:\n",
    "$$\n",
    "    E = \\left( \\sum^m_{i=1} |y_i - f(x_i)|^2 \\right )^{1/2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But how do we do this if our matrix $A$ is now $m \\times n$ and looks like\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "        1 & x_1 & x_1^2 & \\cdots & x_1^{n-1} \\\\\n",
    "        1 & x_2 & x_2^2 & \\cdots & x_2^{n-1} \\\\\n",
    "        \\vdots & \\vdots & \\vdots & & \\vdots \\\\\n",
    "        1 & x_m & x_m^2 & \\cdots & x_m^{n-1}\n",
    "    \\end{bmatrix}?\n",
    "$$\n",
    "\n",
    "Turns out if we solve the system\n",
    "\n",
    "$$A^T A x = A^T b$$\n",
    "\n",
    "we can gaurantee that the error is minimized in the least-squares sense[<sup>1</sup>](#footnoteRegression).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Linear least squares 2\n",
    "\n",
    "Fitting a line through data that has random noise added to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Linear Least Squares Problem\n",
    "\n",
    "# First define the independent and dependent variables.\n",
    "N = 20\n",
    "x = numpy.linspace(-1.0, 1.0, N)\n",
    "y = x + numpy.random.random((N)) \n",
    "\n",
    "# Define the Vandermonde matrix based on our x-values\n",
    "A = numpy.ones((x.shape[0], 2))\n",
    "A[:, 1] = x\n",
    "\n",
    "# Determine the coefficients of the polynomial that will\n",
    "# result in the smallest sum of the squares of the residual.\n",
    "p = numpy.linalg.solve(numpy.dot(A.transpose(), A), numpy.dot(A.transpose(), y))\n",
    "print(\"Error in slope = %s, y-intercept = %s\" % (numpy.abs(p[1] - 1.0), numpy.abs(p[0] - 0.5)))\n",
    "\n",
    "# Plot it out, cuz pictures are fun!\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, y, 'ko')\n",
    "axes.plot(x, p[0] + p[1] * x, 'r')\n",
    "axes.set_title(\"Least Squares Fit to Data\")\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"$f(x)$ and $y_i$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Eigenproblems\n",
    "\n",
    "Eigenproblems come up in a variety of contexts and often are integral to many problem of scientific and engineering interest. It is such a powerful idea that it is not uncommon for us to take a problem and convert it into an eigenproblem. Here we introduce the idea and give some examples.\n",
    "\n",
    "As a review, if $A \\in \\mathbb{C}^{m\\times m}$ (a square matrix with complex values), a non-zero vector $\\vec{v}\\in\\mathbb{C}^m$ is an **eigenvector** of $A$ with a corresponding **eigenvalue** $\\lambda \\in \\mathbb{C}$ if \n",
    "\n",
    "$$A \\vec{v} = \\lambda \\vec{v}.$$\n",
    "\n",
    "One way to interpret the eigenproblem is that we are attempting to ascertain the \"action\" of the matrix $A$ on some subspace of $\\mathbb{C}^m$ where this action acts like scalar multiplication.  This subspace is called an **eigenspace**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Compute the eigenspace of the matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "        1 & 2 \\\\\n",
    "        2 & 1\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Recall that we can find the eigenvalues of a matrix by computing $\\det(A - \\lambda I) = 0$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this case we have\n",
    "$$\\begin{aligned}\n",
    "    A - \\lambda I &= \\begin{bmatrix}\n",
    "        1 & 2 \\\\\n",
    "        2 & 1\n",
    "    \\end{bmatrix} - \\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix} \\lambda\\\\\n",
    "    &= \\begin{bmatrix}\n",
    "        1 - \\lambda & 2 \\\\\n",
    "        2 & 1 - \\lambda\n",
    "    \\end{bmatrix}.\n",
    "\\end{aligned}$$\n",
    "The determinant of the matrix is\n",
    "$$\\begin{aligned}\n",
    "\\begin{vmatrix}\n",
    "    1 - \\lambda & 2 \\\\\n",
    "    2 & 1 - \\lambda\n",
    "\\end{vmatrix} &= (1 - \\lambda) (1 - \\lambda) - 2 \\cdot 2 \\\\\n",
    "&= 1 - 2 \\lambda + \\lambda^2 - 4 \\\\\n",
    "&= \\lambda^2 - 2 \\lambda - 3.\n",
    "\\end{aligned}$$\n",
    "This result is sometimes referred to as the characteristic equation of the matrix, $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Setting the determinant equal to zero we can find the eigenvalues as\n",
    "$$\\begin{aligned}\n",
    "    & \\\\\n",
    "    \\lambda &= \\frac{2 \\pm \\sqrt{4 - 4 \\cdot 1 \\cdot (-3)}}{2} \\\\\n",
    "    &= 1 \\pm 2 \\\\\n",
    "    &= -1 \\mathrm{~and~} 3\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The eigenvalues are used to determine the eigenvectors. The eigenvectors are found by going back to the equation $(A - \\lambda I) \\vec{v}_i = 0$ and solving for each vector.  A trick that works some of the time is to normalize each vector such that the first entry is 1 ($\\vec{v}_1 = 1$):\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "        1 - \\lambda & 2 \\\\\n",
    "        2 & 1 - \\lambda\n",
    "    \\end{bmatrix} \\begin{bmatrix} 1 \\\\ v_2 \\end{bmatrix} = 0\n",
    "$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    1 - \\lambda + 2 v_2 &= 0 \\\\\n",
    "    v_2 &= \\frac{\\lambda - 1}{2}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can check this by\n",
    "$$\\begin{aligned}\n",
    "    2 + \\left(1- \\lambda \\frac{\\lambda - 1}{2}\\right) & = 0\\\\\n",
    "    (\\lambda - 1)^2 - 4 &=0\n",
    "\\end{aligned}$$\n",
    "\n",
    "which by design is satisfied by our eigenvalues.  Another sometimes easier approach is to plug-in the eigenvalues to find each corresponding eigenvector.  The eigenvectors are therefore\n",
    "\n",
    "$$\\vec{v} = \\begin{bmatrix}1 \\\\ -1 \\end{bmatrix}, \\begin{bmatrix}1 \\\\ 1 \\end{bmatrix}.$$\n",
    "\n",
    "Note that these are linearly independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fundamentals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Matrix-Vector Multiplication\n",
    "\n",
    "One of the most basic operations we can perform with matrices is to multiply them be a vector.  This matrix-vector product $A \\vec{x} = \\vec{b}$ is defined as\n",
    "$$\n",
    "    b_i = \\sum^n_{j=1} a_{ij} x_j \\quad \\text{where}\\quad i = 1, \\ldots, m\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Writing the matrix-vector product this way we see that one interpretation of this product is that each column of $A$ is weighted by the value $x_j$, or in other words $\\vec{b}$ is a linear combination of the columns of $A$ where each column's weighting is $x_j$.\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\vec{b} &= A \\vec{x}, \\\\\n",
    "    \\vec{b} &= \n",
    "    \\begin{bmatrix}  &  &  &  \\\\  &  &  &  \\\\ \\vec{a}_1 & \\vec{a}_2 & \\cdots & \\vec{a}_n \\\\  &  &  &   \\\\  &  &  &  \\end{bmatrix}\n",
    "    \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n \\end{bmatrix}, \\\\\n",
    "    \\vec{b} &= x_1 \\vec{a}_1 + x_2 \\vec{a}_2 + \\cdots + x_n \\vec{a}_n.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This view will be useful later when we are trying to interpret various types of matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One important property of the matrix-vector product is that is a **linear** operation, also known as a **linear operator**.  This means that the for any $\\vec{x}, \\vec{y} \\in \\mathbb{C}^n$ and any $c \\in \\mathbb{C}$ we know that\n",
    "\n",
    "1. $A (\\vec{x} + \\vec{y}) = A\\vec{x} + A\\vec{y}$\n",
    "1. $A\\cdot (c\\vec{x}) = c A \\vec{x}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Vandermonde Matrix\n",
    "\n",
    "In the case where we have $m$ data points and want $m - 1$ order polynomial interpolant the matrix $A$ is a square, $m \\times m$, matrix as before.  Using the above interpretation the polynomial coefficients $p$ are the weights for each of the monomials that give exactly the $y$ values of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Numerical matrix-vector multiply\n",
    "\n",
    "Write a matrix-vector multiply function and check it with the appropriate `numpy` routine.  Also verify the linearity of the matrix-vector multiply."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "#A x = b\n",
    "#(m x n) (n x 1) = (m x 1)\n",
    "def matrix_vector_product(A, x):\n",
    "    m, n = A.shape\n",
    "    b = numpy.zeros(m)\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            b[i] += A[i, j] * x[j]\n",
    "    return b\n",
    "\n",
    "m = 4\n",
    "n = 3\n",
    "A = numpy.random.uniform(size=(m,n))\n",
    "x = numpy.random.uniform(size=(n))\n",
    "y = numpy.random.uniform(size=(n))\n",
    "c = numpy.random.uniform()\n",
    "b = matrix_vector_product(A, x)\n",
    "print(numpy.allclose(b, numpy.dot(A, x)))\n",
    "print(numpy.allclose(matrix_vector_product(A, (x + y)), matrix_vector_product(A, x) + matrix_vector_product(A, y)))\n",
    "print(numpy.allclose(matrix_vector_product(A, c * x), c*matrix_vector_product(A, x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Matrix-Matrix Multiplication\n",
    "\n",
    "The matrix product with another matrix $A C = B$ is defined as\n",
    "$$\n",
    "    b_{ij} = \\sum^m_{k=1} a_{ik} c_{kj}.\n",
    "$$\n",
    "Again, a useful interpretation of this operation is that the product result $B$ is the a linear combination of the columns of $A$.\n",
    "\n",
    "What are the dimensions of $A$ and $C$ so that the multiplication works?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Outer Product\n",
    "\n",
    "The product of two vectors $\\vec{u} \\in \\mathbb{C}^m$ and $\\vec{v} \\in \\mathbb{C}^n$ is a $m \\times n$ matrix where the columns are the vector $u$ multiplied by the corresponding value of $v$:\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\vec{u} \\vec{v}^T &= \n",
    "    \\begin{bmatrix}  u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_n  \\end{bmatrix}\n",
    "    \\begin{bmatrix} v_1 & v_2 & \\cdots & v_n \\end{bmatrix}, \\\\\n",
    "    & = \\begin{bmatrix} v_1u_1 & \\cdots & v_n u_1 \\\\ \\vdots &  & \\vdots \\\\ v_1 u_m & \\cdots & v_n u_m \\end{bmatrix}.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It is useful to think of these as operations on the column vectors, and an equivalent way to express this relationship is \n",
    "$$\n",
    "\\begin{align}\n",
    "    \\vec{u} \\vec{v}^T &=\n",
    "    \\begin{bmatrix}  \\\\ \\vec{u} \\\\ \\\\  \\end{bmatrix}\n",
    "    \\begin{bmatrix} v_1 & v_2 & \\cdots & v_n \\end{bmatrix}, \\\\\n",
    "    &=\n",
    "    \\begin{bmatrix}  &  &  &  \\\\  &  &  &  \\\\ \\vec{u}v_1  & \\vec{u} v_2  & \\cdots & \\vec{u} v_n  \\\\  &  &  & \\\\  &  &  &  \\end{bmatrix}, \\\\\n",
    "    & = \\begin{bmatrix} v_1u_1 & \\cdots & v_n u_1 \\\\ \\vdots &  & \\vdots \\\\ v_1 u_m & \\cdots & v_n u_m \\end{bmatrix}.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Upper Triangular Multiplication\n",
    "\n",
    "Consider the multiplication of a matrix $A \\in \\mathbb{C}^{m\\times n}$ and the **upper-triangular** matrix $R$ defined as the $n \\times n$ matrix with entries $r_{ij} = 1$ for $i \\leq j$ and $r_{ij} = 0$ for $i > j$.  The product can be written as\n",
    "$$\n",
    "    \\begin{bmatrix}  \\\\  \\\\ \\vec{b}_1 & \\cdots & \\vec{b}_n \\\\ \\\\ \\\\ \\end{bmatrix} = \\begin{bmatrix} \\\\ \\\\  \\vec{a}_1 & \\cdots & \\vec{a}_n \\\\ \\\\ \\\\ \\end{bmatrix} \\begin{bmatrix} 1 & \\cdots & 1 \\\\  & \\ddots & \\vdots \\\\  &  & 1 \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "The columns of $B$ are then\n",
    "$$\n",
    "    \\vec{b}_j = A \\vec{r}_j = \\sum^j_{k=1} \\vec{a}_k\n",
    "$$\n",
    "so that $\\vec{b}_j$ is the sum of the first $j$ columns of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Write Matrix-Matrix Multiplication\n",
    "\n",
    "Write a function that computes matrix-matrix multiplication and demonstrate the following properties:\n",
    "1. $A (B + C) = AB + AC$ (for square matrices))\n",
    "1. $A (cB) = c AB$ where $c \\in \\mathbb{C}$\n",
    "1. $AB \\neq BA$ in general"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def matrix_matrix_product(A, B):\n",
    "    C = numpy.zeros((A.shape[0], B.shape[1]))\n",
    "    for i in range(A.shape[0]):\n",
    "        for j in range(B.shape[1]):\n",
    "            for k in range(A.shape[1]):\n",
    "                C[i, j] += A[i, k] * B[k, j]\n",
    "    return C\n",
    "\n",
    "\n",
    "m = 4\n",
    "n = 4\n",
    "p = 4\n",
    "A = numpy.random.uniform(size=(m, n))\n",
    "B = numpy.random.uniform(size=(n, p))\n",
    "C = numpy.random.uniform(size=(m, p))\n",
    "c = numpy.random.uniform()\n",
    "print(numpy.allclose(matrix_matrix_product(A, B), numpy.dot(A, B)))\n",
    "print(numpy.allclose(matrix_matrix_product(A, (B + C)), matrix_matrix_product(A, B) + matrix_matrix_product(A, C)))\n",
    "print(numpy.allclose(matrix_matrix_product(A, c * B), c*matrix_matrix_product(A, B)))\n",
    "print(numpy.allclose(matrix_matrix_product(A, B), matrix_matrix_product(B, A)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### NumPy Products\n",
    "\n",
    "NumPy and SciPy contain routines that ware optimized to perform matrix-vector and matrix-matrix multiplication.  Given two `ndarray`s you can take their product by using the `dot` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "n = 10\n",
    "m = 5\n",
    "# Matrix vector with identity\n",
    "A = numpy.identity(n)\n",
    "x = numpy.random.random(n)\n",
    "print(numpy.allclose(x, numpy.dot(A, x)))\n",
    "\n",
    "# Matrix vector product\n",
    "A = numpy.random.random((m, n))\n",
    "print(numpy.dot(A, x))\n",
    "\n",
    "# Matrix matrix product\n",
    "B = numpy.random.random((n, m))\n",
    "print(numpy.dot(A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Range and Null-Space\n",
    "\n",
    "#### Range\n",
    "- The **range** of a matrix $A \\in \\mathbb R^{m \\times n}$ (similar to any function), denoted as $\\text{range}(A)$, is the set of vectors that can be expressed as $A x$ for $x \\in \\mathbb R^n$.  \n",
    "- We can also then say that that $\\text{range}(A)$ is the space **spanned** by the columns of $A$.  In other words the columns of $A$ provide a basis for $\\text{range}(A)$, also called the **column space** of the matrix $A$.  \n",
    "\n",
    "#### Null-Space\n",
    "- Similarly the **null-space** of a matrix $A$, denoted $\\text{null}(A)$ is the set of vectors $x$ that satisfy $A x = 0$.\n",
    "- A similar concept is the **rank** of the matrix $A$, denoted as $\\text{rank}(A)$, is the dimension of the column space.  A matrix $A$ is said to have **full-rank** if $\\text{rank}(A) = \\min(m, n)$.  This property also implies that the matrix mapping is **one-to-one**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inverse\n",
    "\n",
    "A **non-singular** or **invertible** matrix is characterized as a matrix with full-rank.  This is related to why we know that the matrix is one-to-one, we can use it to transform a vector $x$ and using the inverse, denoted $A^{-1}$, we can map it back to the original matrix.  The familiar definition of this is\n",
    "\\begin{align*}\n",
    "    A \\vec{x} &= \\vec{b}, \\\\\n",
    "    A^{-1} A \\vec{x} & = A^{-1} \\vec{b}, \\\\\n",
    "    x &=A^{-1} \\vec{b}.\n",
    "\\end{align*}\n",
    "Since $A$ has full rank, its columns form a basis for $\\mathbb{R}^m$ and the vector $\\vec{b}$ must be in the column space of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are a number of important properties of a non-singular matrix A.  Here we list them as the following equivalent statements\n",
    "1. $A$ has an inverse $A^{-1}$\n",
    "1. $\\text{rank}(A) = m$\n",
    "1. $\\text{range}(A) = \\mathbb{C}^m$\n",
    "1. $\\text{null}(A) = {0}$\n",
    "1. 0 is not an eigenvalue of $A$\n",
    "1. $\\text{det}(A) \\neq 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Properties of invertible matrices\n",
    "\n",
    "Show that given an invertible matrix that the rest of the properties hold.  Make sure to search the `numpy` packages for relevant functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "m = 3\n",
    "for n in range(100):\n",
    "    A = numpy.random.uniform(size=(m, m))\n",
    "    if numpy.linalg.det(A) != 0:\n",
    "        break\n",
    "        \n",
    "print(numpy.dot(numpy.linalg.inv(A), A))\n",
    "print(numpy.linalg.matrix_rank(A))\n",
    "print(\"range\")\n",
    "print(numpy.linalg.solve(A, numpy.zeros(m)))\n",
    "print(numpy.linalg.eigvals(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Orthogonal Vectors and Matrices\n",
    "\n",
    "Orthogonality is a very important concept in linear algebra that forms the basis of many of the modern methods used in numerical computations.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Two vectors are said to be orthogonal if their **inner-product** or **dot-product** defined as\n",
    "$$\n",
    "    < \\vec{x}, \\vec{y} > \\equiv (\\vec{x}, \\vec{y}) \\equiv \\vec{x}^T\\vec{y} \\equiv \\vec{x} \\cdot \\vec{y} = \\sum^m_{i=1} x_i y_i\n",
    "$$\n",
    "Here we have shown the various notations you may run into (the inner-product is in-fact a general term for a similar operation for mathematical objects such as functions).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $\\langle \\vec{x},\\vec{y} \\rangle = 0$ then we say $\\vec{x}$ and $\\vec{y}$ are orthogonal.  The reason we use this terminology is that the inner-product of two vectors can also be written in terms of the angle between them where\n",
    "$$\n",
    "    \\cos \\theta = \\frac{\\langle \\vec{x}, \\vec{y} \\rangle}{||\\vec{x}||_2~||\\vec{y}||_2}\n",
    "$$\n",
    "and $||\\vec{x}||_2$ is the Euclidean ($\\ell^2$) norm of the vector $\\vec{x}$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can write this in terms of the inner-product as well as\n",
    "$$\n",
    "    ||\\vec{x}||_2^2 = \\langle \\vec{x}, \\vec{x} \\rangle = \\vec{x}^T\\vec{x} = \\sum^m_{i=1} |x_i|^2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The generalization of the inner-product to complex spaces is defined as\n",
    "$$\n",
    "    \\langle x, y \\rangle = \\sum^m_{i=1} x_i^* y_i\n",
    "$$\n",
    "where $x_i^*$ is the complex-conjugate of the value $x_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Orthonormality\n",
    "\n",
    "Taking this idea one step further we can say a set of vectors $\\vec{x} \\in X$ are orthogonal to $\\vec{y} \\in Y$ if $\\forall \\vec{x},\\vec{y}$ $< \\vec{x}, \\vec{y} > = 0$.  If $\\forall \\vec{x},\\vec{y}$ $||\\vec{x}|| = 1$ and $||\\vec{y}|| = 1$ then they are also called orthonormal. Note that we dropped the 2 as a subscript to the notation for the norm of a vector. Later we will explore other ways to define a norm of a vector other than the Euclidean norm defined above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Another concept that is related to orthogonality is linear-independence.  A set of vectors $\\vec{x} \\in X$ are **linearly independent** if $\\forall \\vec{x} \\in X$ that each $\\vec{x}$ cannot be written as a linear combination of the other vectors in the set $X$.  \n",
    "\n",
    "An equivalent statement is that there does not exist a set of scalars $c_i$ such that\n",
    "\n",
    "$$\n",
    "    \\vec{x}_k = \\sum^n_{i=1, i \\neq k} c_i \\vec{x}_i.\n",
    "$$\n",
    "\n",
    "Another way to write this is that $\\vec{x}_k \\in X$ is orthogonal to all the rest of the vectors in the set $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This can be related directly through the idea of projection.  If we have a set of vectors $\\vec{x} \\in X$ we can project another vector $\\vec{v}$ onto the vectors in $X$ by using the inner-product.  This is especially powerful if we have a set of linearly-independent vectors $X$, which are said to **span** a space (or provide a **basis** for a space), s.t. any vector in the space spanned by $X$ can be expressed as a linear combination of the basis vectors $X$\n",
    "$$\n",
    "    \\vec{v} = \\sum^n_{i=1} \\, \\langle \\vec{v}, \\vec{x}_i \\rangle \\, \\vec{x}_i.\n",
    "$$\n",
    "Note if $\\vec{v} \\in X$ that \n",
    "$$\n",
    "    \\langle \\vec{v}, \\vec{x}_i \\rangle = 0 \\quad \\forall \\vec{x}_i \\in X \\setminus \\vec{v}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Looping back to matrices, the column space of a matrix is spanned by its linearly independent columns.  Any vector $v$ in the column space can therefore be expressed via the equation above.  A special class of matrices are called **unitary** matrices when complex-valued and **orthogonal** when purely real-valued if the columns of the matrix are orthonormal to each other.  Importantly this implies that for a unitary matrix $Q$ we know the following\n",
    "\n",
    "1. $Q^* = Q^{-1}$\n",
    "1. $Q^*Q = I$\n",
    "\n",
    "where $Q^*$ is called the **adjoint** of $Q$.  The adjoint is defined as the transpose of the original matrix with the entries being the complex conjugate of each entry as the notation implies.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As an example if we have the matrix\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    Q &= \\begin{bmatrix} q_{11} & q_{12} \\\\ q_{21} & q_{22} \\\\ q_{31} & q_{32} \\end{bmatrix} \\quad \\text{then} \\\\\n",
    "    Q^* &= \\begin{bmatrix} q^*_{11} & q^*_{21} & q^*_{31} \\\\ q^*_{12} & q^*_{22} & q^*_{32} \\end{bmatrix}\n",
    "\\end{aligned}\n",
    "$$\n",
    "The important part of being an unitary matrix is that the projection onto the column space of the matrix $Q$ preserves geometry in an Euclidean sense, i.e. preserves the Cartesian distance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Vector Norms\n",
    "\n",
    "Norms (and also measures) provide a means for measure the \"size\" or distance in a space.  In general a norm is a function, denoted by $||\\cdot||$, that maps $\\mathbb{C}^m \\rightarrow \\mathbb{R}$.  In other words we stick in a multi-valued object and get a single, real-valued number out the other end.  All norms satisfy the properties:\n",
    "\n",
    "1. $||\\vec{x}|| \\geq 0$, and $||\\vec{x}|| = 0$ only if $\\vec{x} = \\vec{0}$\n",
    "1. $||\\vec{x} + \\vec{y}|| \\leq ||\\vec{x}|| + ||\\vec{y}||$ (triangle inequality)\n",
    "1. $||c \\vec{x}|| = |c| ~ ||\\vec{x}||$ where $c \\in \\mathbb{C}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are a number of relevant norms that we can define beyond the Euclidean norm, also know as the 2-norm or $\\ell_2$ norm:\n",
    "\n",
    "1. $\\ell_1$ norm:\n",
    "$$\n",
    "    ||\\vec{x}||_1 = \\sum^m_{i=1} |x_i|,\n",
    "$$\n",
    "1. $\\ell_2$ norm:\n",
    "$$\n",
    "    ||\\vec{x}||_2 = \\left( \\sum^m_{i=1} |x_i|^2 \\right)^{1/2},\n",
    "$$\n",
    "1. $\\ell_p$ norm:\n",
    "$$\n",
    "    ||\\vec{x}||_p = \\left( \\sum^m_{i=1} |x_i|^p \\right)^{1/p}, \\quad \\quad 1 \\leq p < \\infty,\n",
    "$$\n",
    "1. $\\ell_\\infty$ norm:\n",
    "$$\n",
    "    ||\\vec{x}||_\\infty = \\max_{1\\leq i \\leq m} |x_i|,\n",
    "$$\n",
    "1. weighted $\\ell_p$ norm:\n",
    "$$\n",
    "    ||\\vec{x}||_{W_p} = \\left( \\sum^m_{i=1} |w_i x_i|^p \\right)^{1/p}, \\quad \\quad 1 \\leq p < \\infty,\n",
    "$$\n",
    "\n",
    "These are also related to other norms denoted by capital letters ($L_2$ for instance).  In this case we use the lower-case notation to denote finite or discrete versions of the infinite dimensional counterparts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Comparisons Between Norms\n",
    "\n",
    "Compute the norms given some vector $\\vec{x}$ and compare their values.  Verify the properties of the norm for one of the norms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "m = 10\n",
    "p = 4\n",
    "x = numpy.random.uniform(size=m)\n",
    "\n",
    "ell_1 = 0.0\n",
    "for i in range(m):\n",
    "    ell_1 += numpy.abs(x[i])\n",
    "    \n",
    "ell_2 = 0.0\n",
    "for i in range(m):\n",
    "    ell_2 += numpy.abs(x[i])**2\n",
    "ell_2 = numpy.sqrt(ell_2)\n",
    "\n",
    "ell_p = 0.0\n",
    "for i in range(m):\n",
    "    ell_p += numpy.abs(x[i])**p\n",
    "ell_p = (ell_2)**(1.0 / p)\n",
    "\n",
    "ell_infty = numpy.max(numpy.abs(x))\n",
    "\n",
    "print(\"L_1 = %s, L_2 = %s, L_%s = %s, L_infty = %s\" % (ell_1, ell_2, p, ell_p, ell_infty))\n",
    "\n",
    "y = numpy.random.uniform(size=m)\n",
    "print()\n",
    "print(\"Properties of norms:\")\n",
    "print(numpy.max(numpy.abs(x + y)), numpy.max(numpy.abs(x)) + numpy.max(numpy.abs(y)))\n",
    "print(numpy.max(numpy.abs(0.1 * x)), 0.1 * numpy.max(numpy.abs(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Matrix Norms\n",
    "\n",
    "The most direct way to consider a matrix norm is those induced by a vector-norm.  Given a vector norm, we can define a matrix norm as the smallest number $C$ that satisfies the inequality\n",
    "$$\n",
    "    ||A \\vec{x}||_{m} \\leq C ||\\vec{x}||_{n}.\n",
    "$$\n",
    "or as the supremum of the ratios so that\n",
    "$$\n",
    "    C = \\sup_{\\vec{x}\\in\\mathbb{C}^n ~ \\vec{x}\\neq\\vec{0}} \\frac{||A \\vec{x}||_{m}}{||\\vec{x}||_n}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Noting that $||A \\vec{x}||$ lives in the column space and $||\\vec{x}||$ on the domain we can think of the matrix norm as the \"size\" of the matrix that maps the domain to the range.  Also noting that if $||\\vec{x}||_n = 1$ we also satisfy the condition we can write the induced matrix norm as\n",
    "$$\n",
    "    ||A||_{(m,n)} = \\sup_{\\vec{x} \\in \\mathbb{C}^n ~ ||\\vec{x}||_{n} = 1} ||A \\vec{x}||_{m}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Induced Matrix Norms\n",
    "\n",
    "Consider the matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix} 1 & 2 \\\\ 0 & 2 \\end{bmatrix}.\n",
    "$$\n",
    "Compute the induced-matrix norm of $A$ for the vector norms $\\ell_2$ and $\\ell_\\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$\\ell^2$: For both of the requested norms the unit-length vectors $[1, 0]$ and $[0, 1]$ can be used to give an idea of what the norm might be and provide a lower bound.  \n",
    "\n",
    "$$\n",
    "    ||A||_2 = \\sup_{x \\in \\mathbb{R}^n} \\left( ||A \\cdot [1, 0]^T||_2, ||A \\cdot [0, 1]^T||_2 \\right )\n",
    "$$\n",
    "\n",
    "computing each of the norms we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\begin{bmatrix} 1 & 2 \\\\ 0 & 2 \\end{bmatrix} \\cdot \\begin{bmatrix} 1 \\\\ 0 \\end{bmatrix} &= \\begin{bmatrix} 1 \\\\ 0 \\end{bmatrix} \\\\\n",
    "    \\begin{bmatrix} 1 & 2 \\\\ 0 & 2 \\end{bmatrix} \\cdot \\begin{bmatrix} 0 \\\\ 1 \\end{bmatrix} &= \\begin{bmatrix} 2 \\\\ 2 \\end{bmatrix}\n",
    "\\end{aligned}$$\n",
    "\n",
    "which translates into the norms $||A \\cdot [1, 0]^T||_2 = 1$ and $||A \\cdot [0, 1]^T||_2 = 2 \\sqrt{2}$.  This implies that the $\\ell_2$ induced matrix norm of $A$ is at least $||A||_{2} = 2 \\sqrt{2} \\approx 2.828427125$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The exact value of $||A||_2$ can be computed using the spectral radius defined as\n",
    "$$\n",
    "    \\rho(A) = \\max_{i} |\\lambda_i|,\n",
    "$$\n",
    "where $\\lambda_i$ are the eigenvalues of $A$.  With this we can compute the $\\ell_2$ norm of $A$ as\n",
    "$$\n",
    "    ||A||_2 = \\sqrt{\\rho(A^\\ast A)}\n",
    "$$\n",
    "\n",
    "Computing the norm again here we find\n",
    "$$\n",
    "    A^\\ast A = \\begin{bmatrix} 1 & 0 \\\\ 2 & 2 \\end{bmatrix} \\begin{bmatrix} 1 & 2 \\\\ 0 & 2 \\end{bmatrix} = \\begin{bmatrix} 1 & 2 \\\\ 2 & 8 \\end{bmatrix}\n",
    "$$\n",
    "which has eigenvalues \n",
    "$$\n",
    "    \\lambda = \\frac{1}{2}\\left(9 \\pm \\sqrt{65}\\right )\n",
    "$$\n",
    "so $||A||_2 \\approx 2.9208096$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$\\ell^\\infty$: We can again bound $||A||_\\infty$ by looking at the unit vectors which give us the matrix lower bound of 2.  To compute it turns out $||A||_{\\infty} = \\max_{1 \\leq i \\leq m} ||a^\\ast_i||_1$ where $a^\\ast_i$ is the $i$th row of $A$.  This represents then the maximum of the row sums of $A$.  Therefore $||A||_\\infty = 3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([[1, 2], [0, 2]])\n",
    "print(numpy.linalg.norm(A, ord=2))\n",
    "print(numpy.linalg.norm(A, ord=numpy.infty))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: General  Norms of a Matrix\n",
    "\n",
    "Compute a bound on the induced norm of the $m \\times n$ dimensional matrix $A$ using $\\ell_1$ and $\\ell_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One of the most useful ways to think about matrix norms is as a transformation of a unit-ball to an ellipse.  Depending on the norm in question, the norm will be some combination of the resulting ellipse.  For the above cases we have some nice relations based on these ideas.\n",
    "\n",
    "1. $||A \\vec{x}||_1 = || \\sum^n_{j=1} x_j \\vec{a}_j ||_1 \\leq \\sum^n_{j=1} |x_j| ||\\vec{a}_j||_1 \\leq \\max_{1\\leq j\\leq n} ||\\vec{a}_j||_1$\n",
    "1. $||A \\vec{x}||_\\infty = || \\sum^n_{j=1} x_j \\vec{a_j} ||_\\infty \\leq \\sum^n_{j=1} |x_j| ||\\vec{a}_j||_\\infty \\leq \\max_{1 \\leq i \\leq m} ||a^*_i||_1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note: that this code is a bit fragile to angles that go beyond pi\n",
    "# due to the use of arccos.\n",
    "\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "A = numpy.array([[1, 2], [0, 2]])\n",
    "\n",
    "def draw_unit_vectors(axes, A, head_width=0.1):\n",
    "    head_length = 1.5 * head_width\n",
    "    image_e = numpy.empty(A.shape)\n",
    "    angle = numpy.empty(A.shape[0])\n",
    "    image_e[:, 0] = numpy.dot(A, numpy.array((1.0, 0.0)))\n",
    "    image_e[:, 1] = numpy.dot(A, numpy.array((0.0, 1.0)))\n",
    "    for i in range(A.shape[0]):\n",
    "        angle[i] = numpy.arccos(image_e[0, i] / numpy.linalg.norm(image_e[:, i], ord=2))\n",
    "        axes.arrow(0.0, 0.0, image_e[0, i] - head_length * numpy.cos(angle[i]), \n",
    "                             image_e[1, i] - head_length * numpy.sin(angle[i]), \n",
    "                             head_width=head_width, color='b', alpha=0.5)\n",
    "    \n",
    "head_width = 0.2\n",
    "head_length = 1.5 * head_width\n",
    "# ============\n",
    "# 1-norm\n",
    "# Unit-ball\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.suptitle(\"1-Norm\")\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1, aspect='equal')\n",
    "axes.plot((1.0, 0.0, -1.0, 0.0, 1.0), (0.0, 1.0, 0.0, -1.0, 0.0), 'r')\n",
    "draw_unit_vectors(axes, numpy.eye(2))\n",
    "axes.set_title(\"Unit Ball\")\n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid(True)\n",
    "\n",
    "# Image\n",
    "axes = fig.add_subplot(1, 2, 2, aspect='equal')\n",
    "axes.plot((1.0, 2.0, -1.0, -2.0, 1.0), (0.0, 2.0, 0.0, -2.0, 0.0), 'r')\n",
    "draw_unit_vectors(axes, A, head_width=0.2)\n",
    "\n",
    "axes.set_title(\"Images Under A\")\n",
    "axes.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# ============\n",
    "# 2-norm\n",
    "\n",
    "# Unit-ball\n",
    "fig = plt.figure()\n",
    "fig.suptitle(\"2-Norm\")\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1, aspect='equal')\n",
    "axes.add_artist(plt.Circle((0.0, 0.0), 1.0, edgecolor='r', facecolor='none'))\n",
    "draw_unit_vectors(axes, numpy.eye(2))\n",
    "axes.set_title(\"Unit Ball\")\n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid(True)\n",
    "\n",
    "# Image\n",
    "# Compute some geometry\n",
    "u, s, v = numpy.linalg.svd(A)\n",
    "theta = numpy.empty(A.shape[0])\n",
    "ellipse_axes = numpy.empty(A.shape)\n",
    "theta[0] = numpy.arccos(u[0][0]) / numpy.linalg.norm(u[0], ord=2)\n",
    "theta[1] = theta[0] - numpy.pi / 2.0\n",
    "for i in range(theta.shape[0]):\n",
    "    ellipse_axes[0, i] = s[i] * numpy.cos(theta[i])\n",
    "    ellipse_axes[1, i] = s[i] * numpy.sin(theta[i])\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2, aspect='equal')\n",
    "axes.add_artist(patches.Ellipse((0.0, 0.0), 2 * s[0], 2 * s[1], theta[0] * 180.0 / numpy.pi,\n",
    "                                edgecolor='r', facecolor='none'))\n",
    "for i in range(A.shape[0]):\n",
    "    axes.arrow(0.0, 0.0, ellipse_axes[0, i] - head_length * numpy.cos(theta[i]), \n",
    "                         ellipse_axes[1, i] - head_length * numpy.sin(theta[i]), \n",
    "                         head_width=head_width, color='k')\n",
    "draw_unit_vectors(axes, A, head_width=0.2)\n",
    "axes.set_title(\"Images Under A\")\n",
    "axes.set_xlim((-s[0] + 0.1, s[0] + 0.1))\n",
    "axes.set_ylim((-s[0] + 0.1, s[0] + 0.1))\n",
    "axes.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# ============\n",
    "# infty-norm\n",
    "# Unit-ball\n",
    "fig = plt.figure()\n",
    "fig.suptitle(\"$\\infty$-Norm\")\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1, aspect='equal')\n",
    "axes.plot((1.0, -1.0, -1.0, 1.0, 1.0), (1.0, 1.0, -1.0, -1.0, 1.0), 'r')\n",
    "draw_unit_vectors(axes, numpy.eye(2))\n",
    "axes.set_title(\"Unit Ball\")\n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid(True)\n",
    "\n",
    "# Image\n",
    "# Geometry - Corners are A * ((1, 1), (1, -1), (-1, 1), (-1, -1))\n",
    "# Symmetry implies we only need two.  Here we just plot two\n",
    "u = numpy.empty(A.shape)\n",
    "u[:, 0] = numpy.dot(A, numpy.array((1.0, 1.0)))\n",
    "u[:, 1] = numpy.dot(A, numpy.array((-1.0, 1.0)))\n",
    "theta[0] = numpy.arccos(u[0, 0] / numpy.linalg.norm(u[:, 0], ord=2))\n",
    "theta[1] = numpy.arccos(u[0, 1] / numpy.linalg.norm(u[:, 1], ord=2))\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2, aspect='equal')\n",
    "axes.plot((3, 1, -3, -1, 3), (2, 2, -2, -2, 2), 'r')\n",
    "for i in range(A.shape[0]):\n",
    "    axes.arrow(0.0, 0.0, u[0, i] - head_length * numpy.cos(theta[i]), \n",
    "                         u[1, i] - head_length * numpy.sin(theta[i]), \n",
    "                         head_width=head_width, color='k')\n",
    "\n",
    "draw_unit_vectors(axes, A, head_width=0.2)\n",
    "axes.set_title(\"Images Under A\")\n",
    "axes.set_xlim((-4.1, 4.1))\n",
    "axes.set_ylim((-3.1, 3.1))\n",
    "axes.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Cauchy-Schwarz and Hölder Inequalities\n",
    "\n",
    "Computing matrix norms where $p \\neq 1$ or $\\infty$ is more difficult unfortunately.  We have a couple of tools that can be useful however.  \n",
    "\n",
    " - **Cauchy-Schwarz Inequality**:  For the special case where $p=q=2$, for any vectors $\\vec{x}$ and $\\vec{y}$\n",
    "$$\n",
    "    |\\vec{x}^*\\vec{y}| \\leq ||\\vec{x}||_2 ||\\vec{y}||_2\n",
    "$$\n",
    " - **Hölder's Inequality**:  Turns out this holds in general if given a $p$ and $q$ that satisfy $1/p + 1/q = 1$ with $1 \\leq p, q \\leq \\infty$\n",
    "\n",
    "$$\n",
    "    |\\vec{x}^*\\vec{y}| \\leq ||\\vec{x}||_p ||\\vec{y}||_q.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### General Matrix Norms (induced and  non-induced)\n",
    "\n",
    "In general matrix-norms have the following properties whether they are induced from a vector-norm or not:\n",
    "1. $||A|| \\geq 0$ and $||A|| = 0$ only if $A = 0$\n",
    "1. $||A + B|| \\leq ||A|| + ||B||$ (Triangle Inequality)\n",
    "1. $||c A|| = |c| ||A||$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The most widely used matrix norm not induced by a vector norm is the **Frobenius norm** defined by\n",
    "$$\n",
    "    ||A||_F = \\left( \\sum^m_{i=1} \\sum^n_{j=1} |A_{ij}|^2 \\right)^{1/2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Invariance under unitary multiplication\n",
    "\n",
    "One important property of the matrix 2-norm (and Forbenius norm) is that multiplication by a unitary matrix does not change the product (kind of like multiplication by 1).  In general for any $A \\in \\mathbb{C}^{m\\times n}$ and unitary matrix $Q \\in \\mathbb{C}^{m \\times m}$ we have\n",
    "\\begin{align*}\n",
    "    ||Q A||_2 &= ||A||_2 \\\\ ||Q A||_F &= ||A||_F.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<sup>1</sup><span id=\"footnoteRegression\"> http://www.utstat.toronto.edu/~brunner/books/LinearModelsInStatistics.pdf</span>"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.9"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}