mp573_lecture_optim_002.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lecture 2: Review of Vectors, Matrices, and Linear Least-Squares Problems\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## MP 573 Part 2, Lecture 2: Review of vectors, matrices, and LLS\n",
    "##\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, Image\n",
    "import matplotlib.image as mpimg\n",
    "from os.path import dirname, join as pjoin\n",
    "import scipy.io as sio\n",
    "import scipy.linalg as slinalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear problems (aka linear systems of equations)\n",
    "\n",
    "In this course, we will encounter multiple instances of linear problems (ie: linear systems of equations):\n",
    "\\begin{equation}\n",
    "    \\mathbf{A} \\mathbf{x} = \\mathbf{y}\n",
    "\\end{equation}\n",
    "Note that these linear systems need not have an exact solution. For instance, does the following linear system (where $\\mathbf{x}$ is a length-1 vector) have a solution?\n",
    "\\begin{equation}\n",
    "      \\left[ {\\begin{array}{c}\n",
    "  1  \\\\\n",
    "  1  \\\\\n",
    "  \\end{array} } \\right] \n",
    "  \\mathbf{x} = \n",
    "        \\left[ {\\begin{array}{c}\n",
    "  1  \\\\\n",
    "  7  \\\\\n",
    "  \\end{array} } \\right] \n",
    "\\end{equation}\n",
    "(Why or why not?)\n",
    "\n",
    "How about the following linear system (where $\\mathbf{x}$ is a length-2 vector): does it have a solution?\n",
    "\\begin{equation}\n",
    "      \\left[ {\\begin{array}{cc}\n",
    "  1 & 1 \\\\\n",
    "  3 & 3 \\\\\n",
    "  \\end{array} } \\right] \n",
    "  \\mathbf{x} = \n",
    "        \\left[ {\\begin{array}{c}\n",
    "  1  \\\\\n",
    "  1  \\\\\n",
    "  \\end{array} } \\right] \n",
    "\\end{equation}\n",
    "(Why or why not?)\n",
    "\n",
    "How about the following linear system (where $\\mathbf{x}$ is a length-2 vector): does it have a solution?\n",
    "\\begin{equation}\n",
    "      \\left[ {\\begin{array}{cc}\n",
    "  1 & 0 \\\\\n",
    "  2 & 1 \\\\\n",
    "  \\end{array} } \\right] \n",
    "  \\mathbf{x} = \n",
    "        \\left[ {\\begin{array}{c}\n",
    "  1  \\\\\n",
    "  1  \\\\\n",
    "  \\end{array} } \\right] \n",
    "\\end{equation}\n",
    "(Why or why not?)\n",
    "\n",
    "In addition, systems of equations oftentimes have multiple solutions. For instance, consider the system of equations (where in this case $\\mathbf{x}$ is a length-3 vector): \n",
    "\\begin{equation}\n",
    "      \\left[ {\\begin{array}{ccc}\n",
    "  1 & 0 & 1 \\\\\n",
    "  2 & 1 & 1 \\\\\n",
    "  \\end{array} } \\right] \n",
    "  \\mathbf{x} = \n",
    "        \\left[ {\\begin{array}{c}\n",
    "  1  \\\\\n",
    "  1  \\\\\n",
    "  \\end{array} } \\right] \n",
    "\\end{equation}\n",
    "Can you give two or more possible solutions $\\mathbf{x}$ to the linear system above? Can you provide a closed form description of *all* possible solutions to this linear system?  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = \n",
      "[[1 0 1]\n",
      " [2 1 1]]\n",
      "\n",
      " y = \n",
      "[[1]\n",
      " [1]]\n",
      "\n",
      " xest = \n",
      "[[ 0.33333333]\n",
      " [-0.33333333]\n",
      " [ 0.66666667]]\n",
      "\n",
      " A xest = \n",
      "[[1.]\n",
      " [1.]]\n"
     ]
    }
   ],
   "source": [
    "# Create our matrix A and vector y\n",
    "# (feel free to replace by one of the matrices/vectors above)\n",
    "A = np.array([[1,0,1], \n",
    "    [2,1,1]])\n",
    "\n",
    "y = np.array([[1],[1]])\n",
    "\n",
    "# Solve the linear system of equations\n",
    "xest = np.linalg.lstsq(A,y,rcond=10)[0]\n",
    "yest = np.matmul(A,xest)\n",
    "\n",
    "\n",
    "\n",
    "print('A = ')\n",
    "print(A)\n",
    "\n",
    "print('\\n y = ')\n",
    "print(y)\n",
    "\n",
    "print('\\n xest = ')\n",
    "print(xest)\n",
    "\n",
    "print('\\n A xest = ')\n",
    "print(yest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "### Question: how does the LLS-solving function in the previous cell pick xest out of all the possible candidate solutions? \n",
    "\n",
    "Let's see if our matrix has a null space, and if so let's modify xest by adding a component within the null space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Basis for null space of A: \n",
      "[[-0.57735027]\n",
      " [ 0.57735027]\n",
      " [ 0.57735027]]\n",
      "\n",
      " xest2: \n",
      "[[-0.09344603]\n",
      " [ 0.09344603]\n",
      " [ 1.09344603]]\n",
      "\n",
      " yest2: \n",
      "[[1.]\n",
      " [1.]]\n"
     ]
    }
   ],
   "source": [
    "# Does our matrix have a null space? \n",
    "ns = slinalg.null_space(A)\n",
    "print('Basis for null space of A: ')\n",
    "print(ns)\n",
    "\n",
    "# Generate alternative solution by adding something \n",
    "# within the nullspace of A to our current solution\n",
    "xest2 = xest + np.matmul(ns,rnd.rand(1,ns.shape[1]))\n",
    "print('\\n xest2: ')\n",
    "print(xest2)\n",
    "\n",
    "# Check that our new solution still fits the data just as well\n",
    "yest2 = np.matmul(A,xest2)\n",
    "print('\\n yest2: ')\n",
    "print(yest2)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q1: A linear least squares problem with a square matrix A\n",
    "\n",
    "Suppose we have a linear least-squares problem defined as $$\\mathbf{x}^* = \\arg \\min_{\\mathbf{x}} \\| \\mathbf{Ax} - \\mathbf{y}\\|^2$$\n",
    "\n",
    "for some matrix $\\mathbf{A} \\in \\mathbb{R}^{N \\times N}$, some data vector $\\mathbf{y} \\in \\mathbb{R}^{N}$, and some unknown vector $\\mathbf{x} \\in \\mathbb{R}^{N}$. \n",
    "\n",
    "Is the following statement correct? Since we know that $\\mathbf{A}$ is a square matrix (ie: there are as many unknowns as there are equations), we know that this LLS problem necessarily has a single unique solution $\\mathbf{x}^*$ that perfectly fits the data, ie: $\\| \\mathbf{A \\, x}^* - \\mathbf{y}\\|^2 = 0$. \n",
    "\n",
    "*Possible answers:*\n",
    "\n",
    "A: Correct\n",
    "\n",
    "B: Not correct. The system may have more than one solution, although we know that any solution will perfectly fit the data, ie: $\\| \\mathbf{A \\, x}^* - \\mathbf{y}\\|^2 = 0$\n",
    "\n",
    "C: Not correct. The system only has one solution, but this solution may not perfectly fit the data, ie: it is possible that $\\| \\mathbf{A \\, x}^* - \\mathbf{y}\\|^2 > 0$\n",
    "\n",
    "D: Not correct. The system may have multiple solutions, and these solutions may or may not fit the data perfectly. \n",
    "\n",
    "E: Not correct. This system has no solution, ever. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q2: LLS solution\n",
    "\n",
    "Using basic calculus, in order to minimize $ \\| \\mathbf{Ax} - \\mathbf{y} \\|_2^2 $ we can look for the choice of $\\mathbf{x}$ (which we call $\\mathbf{x}^*$) such that the gradient of $ \\|  \\mathbf{Ax} - \\mathbf{y} \\|_2^2 $ is zero, ie: \n",
    "\n",
    "\\begin{equation}\n",
    "   \\mathbf{A}^T \\left( \\mathbf{A} \\mathbf{x}^* - \\mathbf{y} \\right) = 0\n",
    "\\end{equation}\n",
    "(let's stick to real-valued vectors and matrices, for simplicity), which can be rewritten as:\n",
    "\\begin{eqnarray}\n",
    "\\mathbf{A}^T \\mathbf{A} \\mathbf{x}^* = \\mathbf{A}^T \\mathbf{y}  \n",
    "\\end{eqnarray}\n",
    "Now, if the matrix $\\mathbf{A}^T \\mathbf{A}$ has an inverse matrix $(\\mathbf{A}^T \\mathbf{A})^{-1}$, we can apply it to both sides of the equation above, ie:\n",
    "\\begin{eqnarray}\n",
    "(\\mathbf{A}^T \\mathbf{A})^{-1} (\\mathbf{A}^T \\mathbf{A}) \\mathbf{x}^* = (\\mathbf{A}^T \\mathbf{A})^{-1}\\mathbf{A}^T \\mathbf{y}  \n",
    "\\end{eqnarray}\n",
    "or in other words,\n",
    "\\begin{equation} \\label{eq:lls_closedsolution}\n",
    "\\mathbf{x}^* = (\\mathbf{A}^T \\mathbf{A})^{-1}\\mathbf{A}^T \\mathbf{y} \n",
    "\\end{equation}\n",
    "And this is our linear least squares solution. \n",
    "\n",
    "*Question:* Suppose we are given one data vector $\\mathbf{y_a}$ and this leads to a solution $\\mathbf{x}_a^*$ as described above. Now suppose we are given a different data vector $\\mathbf{y_b} = 7 \\cdot \\mathbf{y_a}$. What will be the solution to our LLS problem with this new data vector $\\mathbf{y_b}$?\n",
    "\n",
    "A: The same, ie: $\\mathbf{x}_b^* = \\mathbf{x}_a^*$ \n",
    "\n",
    "B: The solution will change linearly with the data ie: $\\mathbf{x}_b^* = 7 \\cdot \\mathbf{x}_a^*$ \n",
    "\n",
    "C: The solution will be all zeros\n",
    "\n",
    "D: We don't have enough information to answer this question - it will depend on further specifics of $\\mathbf{A}$ and $\\mathbf{y_a}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q3: Weighted linear least squares\n",
    "\n",
    "Now we want to elaborate on the previous example, still assuming that the matrix $\\mathbf{A}^T \\mathbf{A}$ has an inverse, as follows: we'd like to minimize $ \\| \\mathbf{W} ( \\mathbf{Ax} - \\mathbf{y}) \\|_2^2 $ for some positive definite matrix $\\mathbf{W}$. For example, let us assume that $ \\mathbf{W}$ is a diagonal matrix where the entries are all strictly positive, ie: $W_{kk} > 0$, but these diagonal entries are generally different from each other. This weighted least squares formulation allows us to include information about which entries in our data vector we really want to fit closely (those entries $y_k$ with high corresponding $W_{kk}$ values), and which ones we do not care about as much (those entries $y_k$ with lower corresponding $W_{kk}$ values). \n",
    "\n",
    "Can you derive a closed form solution to our weighted LLS problem? \n",
    "\n",
    "A: Not in general.\n",
    "\n",
    "B: Yes, we can just call $\\mathbf{z} = \\mathbf{W y}$ our new data vector, $\\mathbf{B} = \\mathbf{W A}$ our new system matrix, and our solution reduces to the LLS solution above, ie: \n",
    "\\begin{equation} \n",
    "\\mathbf{x}^* = (\\mathbf{B}^T \\mathbf{B})^{-1}\\mathbf{B}^T \\mathbf{z} \n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}