01-Fundamentals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vectors, matrices and norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The notebook demonstrate the computation and use of some important concepts in linear algebra. NumPy is used for the numerical computations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vector norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $l_{p}$-norm,of a vector $\\boldsymbol{x} \\in \\mathbb{C}^{n}$ is \n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{x} \\|_{p} = \\left( \\sum_{i=1}^{n} |x_{i}|^{p} \\right)^{1/p}\n",
    "$$\n",
    "\n",
    "Recall that when $p = \\infty$, we have have the maxiumum norm:\n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{x} \\|_{\\infty} = \\max(|x_{1}|, \\ldots , |x_{n}|)\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "NumPy can compute $l_{p}$ norms of vectors. To see how, we first import NumPy and create a random vectors of length 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.4359949 +0.62113383j 0.02592623+0.52914209j 0.54966248+0.13457995j\n",
      " 0.43532239+0.51357812j 0.4203678 +0.18443987j 0.33033482+0.78533515j\n",
      " 0.20464863+0.85397529j 0.61927097+0.49423684j 0.29965467+0.84656149j\n",
      " 0.26682728+0.07964548j]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.seed(2)\n",
    "\n",
    "x = np.random.rand(10) + 1j*np.random.rand(10)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compute a number of $l_{p}$ norms of $\\boldsymbol{x}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The l_1 norm of x is: 6.685801095190599\n",
      "The l_2 norm of x is: 2.202178795407456\n",
      "The l_3 norm of x is: 1.5477495848471308\n",
      "The l_4 norm of x is: 1.3088747768586462\n"
     ]
    }
   ],
   "source": [
    "for p in range(1, 5):\n",
    "    x_norm = np.linalg.norm(x, p)\n",
    "    print(\"The l_{} norm of x is: {}\".format(p, x_norm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the $l_{\\infty}$ norm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The max norm of x is: 0.8980307744980833\n"
     ]
    }
   ],
   "source": [
    " x_inf = np.linalg.norm(x, np.inf)\n",
    "print(\"The max norm of x is: {}\".format(x_inf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Norms of matrices can also be computed. The more interesting (and abstract) norms are *operator* norms. These are also known as *induced* norms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Operator norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an $n \\times n$ matrix $\\boldsymbol{A}$, the norm of the matrix is a measure of the 'maximum change' in relative length it can induce when applied to a vector. If we consider:  \n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{A}  \\boldsymbol{x} \\| \\le C \\| \\boldsymbol{x}\\| \\quad \\forall \\boldsymbol{x} \\in \\mathbb{C}^{d},\n",
    "$$\n",
    "\n",
    "then the smallest possible $C$ is the norm of $\\boldsymbol{A}$. The norm of $\\boldsymbol{A}$ is denoted by $\\|\\boldsymbol{A}\\|$:\n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{A}  \\boldsymbol{x} \\| \\le \\| \\boldsymbol{A}\\| \\| \\boldsymbol{x}\\| \\quad \\forall \\boldsymbol{x} \\in \\mathbb{C}^{d},\n",
    "$$\n",
    "\n",
    "This can be rearranged to provide the usual definition of a matrix norm:\n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{A} \\| = \\max_{\\boldsymbol{x} \\in \\mathbb{C}^{n} \\backslash \\boldsymbol{0}}\n",
    "    \\frac{\\| \\boldsymbol{A} \\boldsymbol{x}\\|}{\\|\\boldsymbol{x}\\| }\n",
    "$$\n",
    "\n",
    "To compute actual norms of a matrix, we need to choose how we measure the length of a vector, i.e. which norm to use. If we choose the $l_{2}$-norm, then:\n",
    "\n",
    "$$\n",
    "\\| \\boldsymbol{A} \\|_{2} = \\max_{\\boldsymbol{x} \\in \\mathbb{C}^{n} \\backslash \\boldsymbol{0}}\n",
    "    \\frac{\\| \\boldsymbol{A} \\boldsymbol{x}\\|_{2}}{\\|\\boldsymbol{x}\\|_{2} }\n",
    "$$\n",
    "\n",
    "As discussed in the lectures, some norms are relatively inexpensive to compute for large matrices, and others are expensive. We can again use NumPy to compute some matrix norms. We first create a matrix filled with random numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.50524609+0.42754596j 0.0652865 +0.43674726j 0.42812233+0.77655918j\n",
      "  0.09653092+0.53560417j 0.12715997+0.95374223j]\n",
      " [0.59674531+0.54420816j 0.226012  +0.08209492j 0.10694568+0.3663424j\n",
      "  0.22030621+0.8508505j  0.34982629+0.40627504j]\n",
      " [0.46778748+0.02720237j 0.20174323+0.24717724j 0.64040673+0.06714437j\n",
      "  0.48306984+0.99385201j 0.50523672+0.97058031j]\n",
      " [0.38689265+0.80025835j 0.79363745+0.60181712j 0.58000418+0.76495986j\n",
      "  0.1622986 +0.16922545j 0.70075235+0.29302323j]\n",
      " [0.96455108+0.52406688j 0.50000836+0.35662428j 0.88952006+0.04567897j\n",
      "  0.34161365+0.98315345j 0.56714413+0.44135492j]]\n"
     ]
    }
   ],
   "source": [
    "A = np.random.rand(5, 5) + 1j*np.random.rand(5, 5)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then compute some norms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The 1-norm of A is: 4.0707141285209385\n",
      "The 2-norm of A is: 3.5641387076264954\n",
      "The max-norm of A is: 4.362031308991134\n"
     ]
    }
   ],
   "source": [
    "print(\"The 1-norm of A is: {}\".format(np.linalg.norm(A, 1)))\n",
    "print(\"The 2-norm of A is: {}\".format(np.linalg.norm(A, 2)))\n",
    "print(\"The max-norm of A is: {}\".format(np.linalg.norm(A, np.inf)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vector-like norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It sometimes convenient to work with matrix norms that are similar to vector norms. A commonly used matrix norm is the Frobenius norm. It is analogous to the $l_{2}$ norm of a vector, and is defined by:\n",
    "\n",
    "$$\n",
    "\\|\\boldsymbol{A} \\|_{F} = \\left( \\sum_{i}\\sum_{i} a_{ij}^{2} \\right)^{1/2}.\n",
    "$$\n",
    "\n",
    "To compute the Frobenius norm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Frobenius norm of A is: 3.875861005158402\n"
     ]
    }
   ],
   "source": [
    "A_frobenius = np.linalg.norm(A, 'fro')\n",
    "print(\"The Frobenius norm of A is: {}\".format(A_frobenius))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Condition number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condition number of a matrix is important when working with matrices numerically because is tells us something about how stable algorithms will be with respect to round-off errors, and how fast some iterative techniques will converge. Recall that the condition number $\\kappa$ of a matrix $\\boldsymbol{A}$ is defined as:\n",
    "\n",
    "$$\n",
    "\\kappa(\\boldsymbol{A}) = \\| \\boldsymbol{A} \\| \\|\\boldsymbol{A}^{-1}\\|\n",
    "$$\n",
    "\n",
    "If we use the 2-norm, it was shown that:\n",
    "\n",
    "$$\n",
    "\\kappa_{2}(\\boldsymbol{A}) = \\frac{\\sqrt{\\lambda_{\\max}(\\boldsymbol{A}^{T}\\boldsymbol{A})}}{\\sqrt{\\lambda_{\\min}(\\boldsymbol{A}^{T}\\boldsymbol{A})}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Effect of poor conditioning on errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was shown in lectures that when solving $\\boldsymbol{A} \\boldsymbol{x} = \\boldsymbol{b}$, if the condition number of $\\boldsymbol{A}$ is large then small errors in $\\boldsymbol{b}$ can manifest themselves as large errors in the solution, $\\boldsymbol{x}$. We explore this now for the notoriously ill-conditioned *Hilbert matrix*. Entries of the Hilbert matrix $\\boldsymbol{H}$ are given by\n",
    "\n",
    "$$\n",
    "H_{ij} = \\frac{1}{i + j - 1}.\n",
    "$$\n",
    "\n",
    "We can use a SciPy function to create a $n \\times n$ Hilbert matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         0.5        0.33333333 0.25       0.2       ]\n",
      " [0.5        0.33333333 0.25       0.2        0.16666667]\n",
      " [0.33333333 0.25       0.2        0.16666667 0.14285714]\n",
      " [0.25       0.2        0.16666667 0.14285714 0.125     ]\n",
      " [0.2        0.16666667 0.14285714 0.125      0.11111111]]\n",
      "Condition number is: 16024897439077.97\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg as la\n",
    "\n",
    "H = la.hilbert(10)\n",
    "print(H[:5,:5])\n",
    "print(\"Condition number is: {}\".format(np.linalg.cond(H, 2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for this small Hilbert matrix, the condition number is large.\n",
    "\n",
    "We now experiment with solving $\\boldsymbol{A} (\\boldsymbol{x}+ \\delta \\boldsymbol{x}) = \\boldsymbol{b} + \\delta \\boldsymbol{b}$, and compare the error $\\|\\delta{\\boldsymbol{x}}\\| / \\|\\boldsymbol{x}\\|$ to $\\|\\delta{\\boldsymbol{b}}\\| / \\|\\boldsymbol{b}\\|$. We will presume that the NumPy linear solvers can cope with the exact system $\\boldsymbol{A}\\boldsymbol{x} =\\boldsymbol{b}$ (in practice this will be an issue).\n",
    "\n",
    "We first construct $\\boldsymbol{b}$, $\\delta\\boldsymbol{b}$ and $\\boldsymbol{b} + \\delta\\boldsymbol{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.ones(H.shape[0])\n",
    "b_delta = 1.0e-6*np.random.rand(H.shape[0])\n",
    "\n",
    "# Perturbed RHS\n",
    "b1  = b + b_delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve for $\\boldsymbol{A} \\boldsymbol{x}=  \\boldsymbol{b}$ and $\\boldsymbol{A} (\\boldsymbol{x}+ \\delta \\boldsymbol{x}) = \\boldsymbol{b} + \\delta \\boldsymbol{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x  = np.linalg.solve(H, b)\n",
    "x1 = np.linalg.solve(H, b1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compare $\\|\\delta{\\boldsymbol{x}}\\| / \\|\\boldsymbol{x}\\|$ and $\\|\\delta{\\boldsymbol{b}}\\| / \\|\\boldsymbol{b}\\|$ using the $l_{2}$-norm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error in x and b: 0.23126319010332205, 5.160707276540258e-07\n"
     ]
    }
   ],
   "source": [
    "error_x = np.linalg.norm(x - x1, 2)/np.linalg.norm(x, 2)\n",
    "error_b = np.linalg.norm(b_delta, 2)/np.linalg.norm(b, 2)\n",
    "print(\"Relative error in x and b: {}, {}\".format(error_x, error_b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for this small Hilbert matrix, a small error in $\\boldsymbol{b}$  leads to a much larger error in $\\boldsymbol{x}$. This will get worse with problem size.  We'll now put the test in side a loop to test larger matrix sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- For 20 x 20 matrix, the condition number is: 1.0843011915344271e+18\n",
      "    Relative error in x and b: 724.9726498675508, 5.333713914529744e-06\n",
      "- For 100 x 100 matrix, the condition number is: 7.2384591535819874e+19\n",
      "    Relative error in x and b: 1170.8226982642284, 5.52552393262059e-06\n",
      "- For 1000 x 1000 matrix, the condition number is: 2.7714433370477366e+20\n",
      "    Relative error in x and b: 1020.3591058646597, 5.688311109647494e-06\n"
     ]
    }
   ],
   "source": [
    "for n in (20, 100, 1000):\n",
    "    H = la.hilbert(n)\n",
    "    print(\"- For {} x {} matrix, the condition number is: {}\".format(n, n, np.linalg.cond(H, 2)))\n",
    "    \n",
    "    b       = np.ones(H.shape[0])\n",
    "    b_delta = 1.0e-5*np.random.rand(H.shape[0])\n",
    "    b1  = b + b_delta\n",
    "    \n",
    "    x  = np.linalg.solve(H, b)\n",
    "    x1 = np.linalg.solve(H, b1)\n",
    "    \n",
    "    error_x = np.linalg.norm(x - x1, 2)/np.linalg.norm(x, 2)\n",
    "    error_b = np.linalg.norm(b_delta, 2)/np.linalg.norm(b, 2)\n",
    "    print(\"    Relative error in x and b: {}, {}\".format(error_x, error_b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Condition number versus determinant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was discussed in lectures that the condition number of a matrix and its determinant are not necessarily related. Some small examples were presented. Here we consider some larger problems.\n",
    "\n",
    "We consider an $n \\times n$ upper triangular matrix filled with two, and one on the diagonal:\n",
    "\n",
    "$$\n",
    "\\boldsymbol{A} = \\begin{bmatrix} \n",
    "1 & 2 & \\ldots & 2 \n",
    "\\\\\n",
    "  & 1 & \\ddots & \\vdots \n",
    "\\\\\n",
    "  &   & \\ddots & 2\n",
    "\\\\\n",
    "   &  &        & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This matrix has a determinant of one, and a condition number that grows with $n$. We can explore this with NumPy for increasing $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- Matrix size: 2 x 2\n",
      "   * l_2 condition number is: 5.828427124746189\n",
      "   * determinant is: 1.0\n",
      "- Matrix size: 10 x 10\n",
      "   * l_2 condition number is: 161.44763879758855\n",
      "   * determinant is: 1.0\n",
      "- Matrix size: 100 x 100\n",
      "   * l_2 condition number is: 16210.722720219874\n",
      "   * determinant is: 1.0\n",
      "- Matrix size: 500 x 500\n",
      "   * l_2 condition number is: 405284.0679028582\n",
      "   * determinant is: 1.0\n"
     ]
    }
   ],
   "source": [
    "def test_matrix(n):\n",
    "    A = np.zeros((n, n))\n",
    "    A[np.triu_indices(n)] = 2.0\n",
    "    np.fill_diagonal(A, 1.0)\n",
    "    return A\n",
    "\n",
    "for n in (2, 10, 100, 500):\n",
    "    A = test_matrix(n)\n",
    "    print(\"- Matrix size: {} x {}\".format(n, n))\n",
    "    print(\"   * l_2 condition number is: {}\".format(np.linalg.cond(A, 2)))\n",
    "    print(\"   * determinant is: {}\".format(np.linalg.det(A)))"
   ]
  }
 ],
 "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}