{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left; width: 300px;\" src=\"images/contourplot.png\">\n",
    "***\n",
    "# <font color=\"grey\">    Lecture 2 - Minima of Convex Functions</font>\n",
    "***\n",
    "$\\newcommand{\\vct}[1]{\\mathbf{#1}}$\n",
    "$\\newcommand{\\mtx}[1]{\\mathbf{#1}}$\n",
    "$\\newcommand{\\e}{\\varepsilon}$\n",
    "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
    "$\\newcommand{\\minimize}{\\text{minimize}\\quad}$\n",
    "$\\newcommand{\\maximize}{\\text{maximize}\\quad}$\n",
    "$\\newcommand{\\subjto}{\\quad\\text{subject to}\\quad}$\n",
    "$\\newcommand{\\R}{\\mathbb{R}}$\n",
    "$\\newcommand{\\trans}{T}$\n",
    "$\\newcommand{\\ip}[2]{\\langle {#1}, {#2} \\rangle}$\n",
    "$\\newcommand{\\zerovct}{\\vct{0}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lecture we will study the unconstrained problem\n",
    "\n",
    "[1]\n",
    "\\begin{equation*}\n",
    " \\minimize\\ f(\\vct{x}),\n",
    "\\end{equation*}\n",
    "\n",
    "where $\\vct{x}\\in \\R^n$. **Optimality conditions** aim to identify properties that potential minimizers need to satisfy in relation to $f(\\vct{x})$. We will review the well known local optimality conditions for differentiable functions from calculus. We then introduce convex functions and discuss some of their properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## <font color=\"grey\">Unconstrained optimization</font>\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solutions to [1] come in different flavours, as in the following definition.\n",
    "\n",
    "**Definition.** A point $\\vct{x}^*\\in \\R^n$ is a\n",
    "\n",
    " * **global minimizer** of [1] if for all $\\vct{x}\\in U$, $f(\\vct{x}^*)\\leq f(\\vct{x})$;\n",
    " * a **local minimizer**, if there is an open neighbourhood $U$ of $\\vct{x}$ such that $f(\\vct{x}^*)\\leq f(\\vct{x})$ for all $\\vct{x}\\in U$;\n",
    " * a **strict local minimizer**, if there is an open neighbourhood $U$ of $\\vct{x}$ such that $f(\\vct{x}^*)<f(\\vct{x})$ for all $\\vct{x}\\in U$;\n",
    " * an **isolated minimizer** if there is an open neighbourhood $U$ of $\\vct{x}^*$ such that $\\vct{x}^*$ is the only local minimizer in $U$.\n",
    "\n",
    "Without any further assumptions on $f$, finding a minimizer is a hopeless task: we simply can't examine the function at *all* points in $\\R^n$. \n",
    "The situation becomes more tractable if we assume some *smoothness* conditions. Recall that $C^k(U)$ denotes the set of functions that are $k$ times continuously differentiable on some set $U$. The following **first-order necessary condition** for optimality is well known. We write $\\nabla f(\\vct{x})$ for the gradient of $f$ at $\\vct{x}$, i.e., the vector \n",
    "\n",
    "\\begin{equation*}\n",
    " \\nabla f(\\vct{x}) = \\left(\\frac{\\partial f}{\\partial x_1}(\\vct{x}),\\dots,\\frac{\\partial f}{\\partial x_n}(\\vct{x})\\right)^{\\trans}\n",
    "\\end{equation*}\n",
    "\n",
    "**Theorem.** Let $\\vct{x}^*$ be a local minimizer of $f$ and assume that $f\\in C^1(U)$ for a neighbourhood of $U$ of $\\vct{x}^*$. Then $\\nabla f(\\vct{x}^*) = \\zerovct$. \n",
    "\n",
    "There are simple examples that show that this is not a sufficient condition: maxima and saddle points will also have a vanishing gradient. If we have access to *second-order information*, in form of the second derivative, or Hessian, of $f$, then we can say more. Recall that the Hessian of $f$ at $\\vct{x}$, $\\nabla^2f(\\vct{x})$, is the $n\\times n$ symmetric matrix given by the second derivatives,\n",
    "\n",
    "\\begin{equation*}\n",
    " \\nabla^2f(\\vct{x}) = \\left(\\frac{\\partial^2 f}{\\partial x_i \\partial x_j}\\right)_{1\\leq i,j\\leq n}.\n",
    "\\end{equation*}\n",
    "\n",
    "In the one-variable case we know that if $x^*$ is a local minimizer of $f\\in C^2([a,b])$, then $f'(x^*)=0$ *and* $f''(x^*)\\geq 0$. Moreover, the conditions $f'(x^*)=0$ *and* $f''(x^*)>0$ guarantee that we have a local minimizer. These conditions generalise to higher dimension, but first we need to know what $f''(x)>0$ means when we have more than one variable.\n",
    "\n",
    "Recall that a matrix $\\mtx{A}$ is **positive semidefinite**, written $\\mtx{A}\\succeq \\zerovct$, if for every $\\vct{x}\\in \\R^n$, $\\vct{x}^{\\top}\\mtx{A}\\vct{x}\\geq 0$, and **positive definite**, written $\\mtx{A}\\succ \\zerovct$, if $\\vct{x}^{\\top}\\mtx{A}\\vct{x}>0$. The property that the Hessian matrix is positive semidefinite is a multivariate generalization of the property that the second derivative is nonnegative. The known conditions for a minimizer involving the second derivative generalise accordingly.\n",
    "\n",
    "**Theorem.** Let $f\\in C^2(U)$ for some open set $U$ and $\\vct{x}^*\\in U$. \n",
    " If $\\vct{x}^*$ is a local minimizer, then $\\nabla f(\\vct{x}^*)=0$ and $\\nabla^2f(\\vct{x}^*)$  is positive semidefinite. Conversely, if $\\nabla f(\\vct{x}^*)=\\zerovct$ and $\\nabla^2f(\\vct{x}^*)$ is positive definite, then $\\vct{x}^*$ is a strict local minimizer. \n",
    "\n",
    "Unfortunately, the above criteria are not able to identify global minimizers, as differentiability is a local property. If, however, the function is **convex**, then we can say a lot more!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## <font color=\"grey\">Convex sets and functions</font>\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now come to the central notion of this course. \n",
    "\n",
    "**Definition.** A set $C\\subseteq \\R^n$ is **convex** if for all $\\vct{x},\\vct{y}\\in C$ and $\\lambda \\in [0,1]$, the line $\\lambda \\vct{x}+(1-\\lambda)\\vct{y}\\in C$. A **convex body** is a convex set that is closed and bounded.\n",
    "\n",
    "**Definition.**\n",
    "Let $S\\subseteq \\R^n$. A function $f\\colon S\\to \\R$ is called **convex** if $S$ is convex and for all $\\vct{x},\\vct{y}\\in \\R^n$ and $\\lambda\\in [0,1]$,\n",
    "\n",
    "\\begin{equation*}\n",
    " f(\\lambda \\vct{x}+(1-\\lambda)\\vct{y})\\leq \\lambda f(\\vct{x})+(1-\\lambda)f(\\vct{y}).\n",
    "\\end{equation*}\n",
    "\n",
    "The function $f$ is called **strictly convex** if\n",
    "\n",
    "\\begin{equation*}\n",
    " f(\\lambda \\vct{x}+(1-\\lambda)\\vct{y})< \\lambda f(\\vct{x})+(1-\\lambda)f(\\vct{y}).\n",
    "\\end{equation*}\n",
    "\n",
    "A function $f$ is called **concave**, if $-f$ is convex. \n",
    "\n",
    "The next figure illustrates how a convex set in $\\R^2$ and a function of one variable looks like. The graph of the function lies below any line connecting two points on it.\n",
    "\n",
    "![Convex function](images/convset.png)\n",
    "\n",
    "Convex function have pleasant properties, while at the same time covering many of the functions that arise in applications. Perhaps the most important property is that local minima are global minima.\n",
    "\n",
    "**Theorem.** Let $f\\colon \\R^n\\to \\R$ be a convex function. Then any local minimizer of $f$ is a global minimizer.\n",
    "\n",
    "**Proof.** Let $\\vct{x}^*$ be a local minimizer and assume that it is not a global minimizer. Then there exists a vector $\\vct{y}\\in \\R^n$ such that $f(\\vct{y})<f(\\vct{x}^*)$. Since $f$ is convex, for any $\\lambda\\in [0,1]$ and $\\vct{x}=\\lambda \\vct{y}+(1-\\lambda)\\vct{x}^*$ we have\n",
    "\n",
    " \\begin{equation*}\n",
    "  f(\\vct{x}) \\leq \\lambda f(\\vct{y})+(1-\\lambda) f(\\vct{x}^*) < \\lambda f(\\vct{x}^*)+(1-\\lambda)f(\\vct{x}^*) = f(\\vct{x}^*).\n",
    " \\end{equation*}\n",
    "\n",
    "This holds for all $\\vct{x}$ on the line segment connecting $\\vct{y}$ and $\\vct{x}^*$. Since every open neighbourhood $U$ of $\\vct{x}^*$ contains a bit of this line segment, this means that every open neighbourhood $U$ of $\\vct{x}^*$ contains an $\\vct{x}\\neq \\vct{x}^*$ such that $f(\\vct{x})\\leq f(\\vct{x}^*)$, in contradiction to the assumption that $\\vct{x}^*$ is a local minimizer. It follows that $\\vct{x}^*$ has to be a global minimizer.\n",
    "\n",
    "** Remark.** Note that in the above theorem we made no assumptions about the differentiability of the function $f$! In fact, while a convex function is always *continuous*, it need not be differentiable. The function $f(x) = |x|$ is a typical example: it is convex, but not differentiable at $x=0$.\n",
    "\n",
    "**Example.** Affine functions $f(\\vct{x})=\\ip{\\vct{x}}{\\vct{a}}+b$ and the exponential function $e^x$ are examples of convex functions. \n",
    "\n",
    "**Example.** In optimization we will often work with functions of matrices, where an $m\\times n$ matrix is considered as a vector in $\\R^{m\\times n}\\cong \\R^{mn}$. If the matrix is symmetric, that is, if $\\mtx{A}^{\\trans}=\\mtx{A}$, then we only care about the upper diagonal entries, and we consider the space $\\mathcal{S}^n$ of symmetric matrices as a vector space of dimension $n(n+1)/2$ (the number of entries on and above the main diagonal). Important functions on symmetric matrices that are convex are the operator norm $\\norm{\\mtx{A}}_2$, defined as\n",
    " \\begin{equation*}\n",
    "  \\norm{\\mtx{A}}_2 := \\max_{\\vct{x}\\colon \\norm{\\vct{x}}\\leq 1} \\frac{\\norm{\\mtx{A}\\vct{x}}_2}{\\norm{\\vct{x}}_2},\n",
    " \\end{equation*}\n",
    "or the function $\\log \\det(\\mtx{X})$, defined on the set of *positive semidefinite* symmetric matrices $\\mathcal{S}_+^n$.\n",
    "\n",
    "There are useful ways of characterising convexity using differentiability.\n",
    "\n",
    "**Theorem.**\n",
    "* Let $f\\in C^1(\\R^n)$. Then $f$ is convex if and only if for all $\\vct{x}, \\vct{y}\\in \\R^n$,\n",
    "\n",
    " \\begin{equation*}\n",
    "  f(\\vct{y})\\geq f(\\vct{x})+\\nabla f(\\vct{x})^{\\trans} (\\vct{y}-\\vct{x}).\n",
    " \\end{equation*}\n",
    "\n",
    "* Let $f\\in C^2(\\R^n)$. Then $f$ is convex if and only if $\\nabla^2 f(\\vct{x})$ is positive semidefinite. If $\\nabla^2f(\\vct{x})$ is positive definite, then $f$ is strictly convex.\n",
    "\n",
    "**Example.** Consider a quadratic function of the form\n",
    " \n",
    " \\begin{equation*}\n",
    "f(\\vct{x}) = \\frac{1}{2}\\vct{x}^{\\trans}\\mtx{A}\\vct{x}+\\vct{b}^{\\trans}\\vct{x}+c, \n",
    " \\end{equation*}\n",
    " \n",
    " where $\\mtx{A}\\in \\R^{n\\times n}$ is symmetric. Writing out the product, we get\n",
    " \n",
    " \\begin{equation*}\n",
    "   \\mtx{x}^{T}\\mtx{A}\\vct{x} = \\begin{pmatrix} x_1 & \\cdots & x_n\n",
    "   \\end{pmatrix}\n",
    "   \\begin{pmatrix}\n",
    "   a_{11} & \\cdots & a_{1n}\\\\\n",
    "   \\vdots & \\ddots & \\vdots\\\\\n",
    "   a_{n1} & \\cdots & a_{nn}\n",
    "   \\end{pmatrix}\n",
    "   \\begin{pmatrix}\n",
    "   x_1\\\\ \\vdots \\\\ x_n\n",
    "   \\end{pmatrix} = \\begin{pmatrix}\n",
    "   x_1 & \\cdots & x_n\n",
    "   \\end{pmatrix}\n",
    "   \\begin{pmatrix}\n",
    "   a_{11}x_1+\\cdots+a_{1n}x_n\\\\\n",
    "   \\vdots \\\\\n",
    "   a_{n1}x_1+\\cdots+a_{nn}x_n\n",
    "   \\end{pmatrix} = \\sum_{i=1}^n \\sum_{j=1}^n a_{ij} x_i x_j.\n",
    " \\end{equation*} \n",
    " \n",
    " Because $\\mtx{A}$ is symmetric, we have $a_{ij}=a_{ji}$, and the above product simplifies to\n",
    " \n",
    " \\begin{equation*}\n",
    " \\vct{x}^{T}\\mtx{A}\\vct{x} = \\sum_{i=1}^n a_{ii} x_i^2 + 2\\sum_{1\\leq i<j\\leq n} a_{ij}x_i x_j.\n",
    " \\end{equation*}\n",
    " \n",
    " This is a quadratic function, because it involves products of the $x_i$. The gradient and the Hessian of $f(\\vct{x})$ are found by computing the first and second order partial derivatives of $f$:\n",
    " \n",
    " \\begin{equation*}\n",
    " \\frac{\\partial f}{\\partial x_i} = \\sum_{j=1}^n a_{ij} x_j + b_i, \\quad \\frac{\\partial^2 f}{\\partial x_i \\partial x_j} = a_{ij}.\n",
    " \\end{equation*}\n",
    " \n",
    " In summary, we have\n",
    " \\begin{equation*}\n",
    "   \\nabla f(\\vct{x}) = \\mtx{A}\\vct{x} + \\vct{b}, \\quad \\nabla^2f(\\vct{x}) = \\mtx{A}.\n",
    " \\end{equation*}\n",
    " \n",
    " Using the previous theorem, we see that $f$ is convex **if and only if** $\\mtx{A}$ is positive semidefinite. A typical example for such a function is the least squares function encountered in Lecture 1,\n",
    " \n",
    " \\begin{equation*}\n",
    "   f(\\vct{x}) = \\|\\mtx{X}\\vct{\\beta}-\\vct{y}\\|^2 = (\\mtx{X}\\vct{\\beta}-\\vct{y})^{T}(\\mtx{X}\\vct{\\beta}-\\vct{y}) = \\vct{\\beta}^{T}\\mtx{X}^{T}\\mtx{X}\\vct{\\beta} -2\\vct{y}^{T}\\mtx{X}\\vct{\\beta}+\\vct{y}^T\\vct{y}.\n",
    " \\end{equation*}\n",
    " \n",
    " The matrix $\\mtx{A}=\\mtx{X}^T\\mtx{X}$ is always symmetric and positive semidefinite (why?) so that the function $f$ is convex. In the example below, we plot the mean squared error of a random least squares problem (see Lecture 1 for an explanation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A convenient way to visualise a function $f\\colon \\R^2\\to \\R$ is through **contour plots**. A **level set** of the function $f$ is a set of the form\n",
    "\n",
    "\\begin{equation*}\n",
    "  \\{\\vct{x}\\mid f(\\vct{x}) = c\\},\n",
    "\\end{equation*}\n",
    "\n",
    "where $c$ is the **level**. Each such level set is a curve in $\\R^2$, and a contour plot is a plot of a collection of such curves for various $c$. If one colours the areas between adjacent curves, one gets a plot as in the following figure. A *convex function* is has the property that there is only one *sink* in the contour plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAD8CAYAAAC4uSVNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFetJREFUeJzt3X+MZWV9x/H3BxYSjLpgTLDZrS6spa5IXHBhqSiroHWD\ncU2NTfxRjJpgY1WIGoICEUJMBbSxqG2CUTAaiVE0FFtUlkjGtroLC2yzwlh+CArEH11TbAwRQb79\nY+7Mzs7cuXPuPefc53nO+bwSkntnzj3zsLPz3u8899wZRQRmZlaWQ1IvwMzMxud4m5kVyPE2MyuQ\n421mViDH28ysQI63mVmBGom3pLWSviFpVtLdkrY2cV4zMxtuTUPnuQq4KSL+WtIa4BlLD5DkC8rN\nzCYQEVr6NtV9kY6kZwN3RcTGVY6LO844beS5Nm4/EYBP3Lybj/5lc8P74Se/vLFzzfv4l67n4ne+\neeH+L/70lMY/xmrufuwFqx5z3dWX8ba//Vir67hztv6/y7d+6zJe/aZ21znK3l0PVTpu9rbPsOmU\nc9tdTAO8zma1vc7Np25Y8X2Xnn340Hg3sW1yDLBf0rWS7pT0eUlHjHuS+XA3rY1wL5VruKehiXCn\nVjXcZm0YFe5Rmoj3GuAk4J8i4iTgceAj45zA4R6Pw92MvbsecritWE3seT8CPBwRewb3rwcuGHbg\n1T/9+cLtlx21li1HrR0a7ldsXNfAstp1+uYXFxHuE162rZV1NB3uDZvaWedKJo32c9eV8Vy819ms\nttY5bOp+cHaGh2ZnVn1s7T1vAEkzwDkRca+kS4BnRMQFS45Ztufd1sQN7U/dJYS7LV2YuM1yUGXL\nZKU976auNjkX+Kqkw4CfAu9q6LwTmcZ2ybQ53M1wuK0rGol3RPwXcPI4j/E+d3UOdzMcbsvJpE9U\nzkvyCkuHuzqHuxkOt3XN1OPtcFfncDfD4bYummq8He7qHO5mONyWo7pbJtCBH0zlcLfH4TbLV1NX\nm3RWX8NdMkfb+qDoybtr13LnFO5Sp26H23LXxJYJFBzvroU7Jw63Wf6KjLdfhNMeh9usPU1N3VBg\nvP0EZXscbrNyFBVvh7s9DrdZWYqJt8PdHofbrH1NbplAQfFum8NdFofb+q6IeHftyhKHux6H20rT\n9NQNBcTb4W6Hw21Wtqzj7XC3w+E2m542pm7ION5dvJY7Bw63WTdkGW9fWdIOh9usO7KLt8PdDofb\nbPra2jKBDOPdNoe7HA632coai7ekQyTdKenGSc/hJyhtnsNtpWtz6oZmJ+/zgHsmfbDD3Y4Sp26H\n22x1jcRb0nrgLOALkzze4W6Hw22WRttTNzQ3eX8aOB8YuxZduyTQ4Z6cw21WXe14S3o98KuI2Ato\n8F8lXbyyJAcOt1k605i6oZnfYXkasEPSWcARwLMkfTki3rH0wE/cvHvh9is2ruPMkxv46CP0cbvE\n4TYr24OzMzw0O7PqcYpo7otd0jbgwxGxY8j74rEr379w3/vczXO4zdJqY+q+9OzDiYhlOxpJrvN2\nuJvncJv1SxPbJgsiYgYYOe873M1zuM36Z6qTt8PdPIfbLA/TeqJyXu9eHt8Uh3syDrdZMzoT72lO\n3TmEu0QOt3XVtKdu6Ei8fS13/hxus2YVH2/vc+fP4bYuSzF1Q+Hxdrjz53CbtaPYeDvc+XO4retS\nTd1QcLynyeEen8Nt1q4i4923K0scbrP8pJy6ocB4O9x5c7itD1KHGwqLdx8vCSyJw202PcXE209Q\n5s3htr7IYeqGQuLtcOfN4TabviLiPU0Ot5mtJJepGwqIt5+gzJunbrM0so63w503h9v6JKepGzKO\nd9+uLHG4zWwcWca7j09QlsThtr7JbeqGTOM9TTmEu6Sp2+E2y0PteEtaL+n7ku6WtE/SuXXO533u\nfDnc1kc5Tt3QzC8gfgr4UETslfRM4A5JN0fET8Y9kcOdL4fb+ijXcEMDk3dE/DIi9g5u/w6YBdaN\nex6HO18Ot1l+Gt3zlrQB2AzsHudxvrIkXw639VXOUzc0GO/Blsn1wHmDCTxLqaduh9vMmtDEnjeS\n1jAX7q9ExL+sdNzHv3T9wu3TN7+Y0ze/uHfbJWaWv5RT94OzMzw0O7PqcYqoPwlK+jKwPyI+NOKY\nePzW6w56W9/C7anbrAw5bZlcevbhRISWvr2JSwVPA94OnCHpLkl3Stq+2uMc7nw53NZnOYV7lNrb\nJhHxn8Ch4zzGT1Dmy+G2Pisl3JDgFZZ9e+m7w21mbej0y+Md7uocbuu7kqZumHK8+7bPXQqH26w8\nnZy8cwh3KVO3w21W3tQNHY13aqWE28zK1bl4p566Swq3p26zMqdu6Fi8He7qHG6zcsMNHYp36nCX\nxOE2K18n4p1DuEuZuh1uszklT93QkXin5nCblaX0cEMH4p166na4zSyFouPtcJvZuEqauk/atOyH\nCS4oNt4Od3Weus3KMyrcUGi8U4e7JA632QGlTN2rhRsKjXdqpUzdDrfZAaWEu6ri4p166na4zaxN\nVaZuKCzeDnc1DrfZwUqZuquGGwqKt8NtZpMoJdzjKiLeqcNdEk/dZmUaZ+qGhuItabukn0i6V9IF\nTZwzJ6VM3Q632cFKmbrHDTc089vjDwE+B7wOOB54q6QX1T3vvNRTt8NtZjlqYvI+BbgvIn4WEU8C\nXwPe2MB5He6KHG6z5bo8dUMz8V4HPLzo/iODt9WSOtxmVq6uhxsKecIyBU/dZpazNQ2c41Hg+Yvu\nrx+8bZmrrrpq4fbWrVs59dRTh54w9dTtcJuVq/Spe9+eGfbdMbPq4xVRL1SSDgX+GzgT+AVwG/DW\niJhdclw8cP/9q57P4a7G4TZbrvRwD7Njy2FExLIH1J68I+KPkt4P3MzcNswXl4a7Koe7GofbzJrY\nNiEivgv8eZ1zpA63mZWti1P3KH7CcsBTt1m5+hZuyCTeqaduh9vM2tZkuCGDeDvc1TjcZsOVMnU3\nLWm8U4e7FA632XClhLvpqRsymLxTKmXqNrPl+hxuSBjv1FN3KeH21G1mwySJt8NdjcNtNlzfp25I\nEO/U4S6Fw202nMM9p3d73qVM3WZmo0w13qmn7lLC7anbbDhP3Qf0ZvJ2uM3K5nAfrDfxLoHDbVa2\naYUbehLvUqZuMxuulKl7mjof71LC7anbbLhSwj3NqRs6Hm+H26xsDvfKOh3vEjjcZjaJzsa7lKnb\nzIbz1D1aJ+NdSrg9dZsN53CvrnPxdrjNbBpShhtqxlvSlZJmJe2V9E1Jz25qYV3mcJutrJSpO7W6\nk/fNwPERsRm4D/ho/SVNrpSp28yGKyXcqaduqBnviLglIp4e3N0FrK+/pMmUEm5P3WbDOdzjaXLP\n+93Adxo8X2UOt5lNQy7hBliz2gGSdgJHL34TEMBFEfHtwTEXAU9GxHWjznXd1Zct3D7hZds4Ycu2\nSdZsZh1TytQ9Dfv2zLDvjplVj1NEvalV0juBc4AzIuKJEcfFjXuerPWxhvHUbVa2UsKdaureseUw\nImLZB1918h5F0nbgfOD0UeFui8NtVjaHe3J197w/CzwT2CnpTkn/3MCaOsXhNitbjuGGmpN3RPxZ\nUwsZVylTt5kNV8rUnasiX2FZSrg9dZsNV0q4U0/dxx/5sxXfV1y8HW6zsjnc1YwKNxQYbzOztqUO\ndxVFxdtTt1nZSpm6U1tt6oaC4u1wm5WtlHCnnrqrhBsKincJHG6z4RzuaqqGGwqJdylTt5mVK3W4\nx5V9vEsJt6dus+FKmbpTG2fqhgLiXQKH22y4UsKdeuoeN9yQebxLmbrNbDmHu5pJwg0Zx7uUcHvq\nNlvO4W5flvF2uM2sDyaduiHTeJtZuTx1V1Mn3JBhvD11m5XL4a6mbrghs3g73GblcriraSLckFm8\nzaxMDvf0ZRNvT91m1nVNTd2QSbwdbrNyeequpslwQybxNrMyOdzVNB1uaCjekj4s6WlJzxn3sZ66\nzcrkcKdVO96S1gOvBcb+p6WUcJvZwUoJdw7amLqhmcn708D5DZwnW566zcqUeupuK9xQM96SdgAP\nR8S+cR9bytTtcJsdrJSpu8vhBliz2gGSdgJHL34TEMDFwIXMbZksft+Krrv6MgB+sT/YsGkbx2za\nNu56zSwhh7uaOuHetWsXu3fvXvU4RUw2AUt6CXAL8Dhz0V4PPAqcEhG/HnJ83LjnScBTt1mJHO7q\nmpy6N77whUTEsv+pVSfvlUTEj4Hnzd+X9CBwUkT876jHOdxm5Skl3Dloe7tkXpPXeQerbJuUEm4z\nO6CkcKeeuqcVbqgxeS8VEcc2da7UPHWbladP4Qa/wnIZh9vsgFKm7r6FGxxvM1uBw503x3sRT91m\nc0oJdw5STN3geJvZEiWFO/XUnSrc4Hgv8NRt5nCPI2W4wfEGHG4zcLjHkTrc4HibWWFShzsXvY+3\np26zcqbuHMI9zan7Tx6+bcX39T7eZn1XSrhzkEu4oefx9tRtfVdSuFNP3Tnscy/W23g73NZ3Dnd1\n0w73alM39DjeZn3mcOerSrihp/H21G195nCPJ6d97sV6GW+zvnK4x5NruKGH8fbUbWZV5PYE5VK9\ni7dZX3nqri7HJyiX6lW8PXVbXznc+Zok3NCjeDvc1lcO93hy3uderDfxNusjh3s8pYQbGoi3pA9I\nmpW0T9Lldc/XBk/d1kcO93hKCjfU/AXEkl4FvAE4ISKekvTc2isys9pKCncOcr+yZJi6k/d7gcsj\n4imAiNhff0nN8tRtfVNauFNP3SVcWTJM3XgfB5wuaZekWyVtaWJRZjYZhztvTYUbKmybSNoJHL34\nTUAAFw8ef1REnCrpZODrwLErnevWb122cHvDpm0cs2nbhMuuxlO39YnDPb4c97l/sPcefrD3nlWP\nU0RMvBhJNwFXRMTM4P79wNaI+M2QY+PSr/xh4o81Ccfb+qSkeDvc1T3j1W8jIpb9gdXdNrkBOANA\n0nHAYcPCnYLDbX3icI+nxCcol6p1tQlwLXCNpH3AE8A76i+pPofb+sThHk+pT1AuVSveEfEkcHZD\nazGzMTnc4+lKuKGDr7D01G19UVK4+6jNcEMH423WB6WFu29Td9vhho7F21O39YHDPb6uhRs6Fm+z\nrnO4x9eFK0uG6Uy8PXVb1znc4+vSE5RLdSbeZl3mcI+vy+GGjsTbU7d1mcM9vq6HGzoSb7Oucrjz\nlyLc0IF4e+q2riot3Lno4pUlwxQfb7MuKjHcOUzdXbuy5A+3/3DF9xUdb0/d1kUO92S6ts89KtxQ\neLzNusbhnkzfwg0Fx9tTt3WNwz2ZPoYbCo63WZc43JPpa7ih0Hh76rYucbjLkPLKkmGKjLdZVzjc\nk+vaJYHjTN3geJsl43BPrk+XBK6kuHh7y8S6wOGeXJ/3uRerFW9JL5X0I0l3SbpN0pY65zPrgxLD\nnQuH+4C6k/eVwCURcSJwCfDJmucbyVO3la7UcOcwdTvcB6sb76eBtYPbRwKP1jyfWWc53JPrWrib\nUOu3xwMfBL4n6R8AAS+vv6ThPHVbyRzuyXUx3HWnbqgQb0k7gaMXvwkI4CLgNcB5EXGDpDcD1wCv\nrb0qsw5xuMtRSrgBFBGTP1h6LCKOXHT/txGxdoVjY9tfXbxwf8OmbRyzaVulj+Op20rlcNfTtam7\nSrj//YFH+I8HDuxAX3HL7UTEsk9I3W2TRyVti4gZSWcC9446+NVv+ljND2dWDoe7nj6GG+CVG9fz\nyo3rF+5fccvtQ4+rG+9zgM9IOhT4PfCemucz6wSHu56+hnscteIdET8EWr2221smVhqHux6Hu5ri\nXmFpljOHu56uhbtNWcfbU7eVxOGup4vhbmvqhszjbVYKh7seh3t82cbbU7eVwuGup2s/IRDaDzdk\nHG+zEjjc9aQId6lPUC7leJtNyOEuT1fCDZnG21smlrtSw52Tru1zTzPckGm8zXJWcrhzmbq7Fu4U\n6r7C0qw3So42ONxtamvqfuC7d634vuwmb2+ZWI4c7mY43NWNCjdkGG+z3DjczXC4q1st3JBZvD11\nW24c7mY43NVVCTdkFm+znDjczfCLcKqrGm5wvM2Gcrib4RfhtCebeHvLxHLhcDeji+Fu0zhTN2QU\nb7McONzN6Gq4c9gumZdlvP/n0d2pl7CqEtYIXuc4qoT7wdmZ9hcyocXh3rcn3TrHCfeuXbsa+Zht\nh/sHe+/JKtyQSbyXbpnsz+ALeTUlrBG8zqqqTtwPZRrvpRP3vjvyXOdSu3fX/7xPY+K+9d++18p5\nJw03+BWW1nOlb5NAPlsl4EsCx1En3JDJ5G2WgsPdLId7uhQR0/lA0nQ+kJlZx0TEsn+lpxZvMzNr\njrdNzMwK5HibmRUoy3hLeqmkH0m6S9JtkrakXtNKJH1A0qykfZIuT72eUSR9WNLTkp6Tei3DSLpy\n8Ge5V9I3JT079ZrmSdou6SeS7pV0Qer1DCNpvaTvS7p78Pfx3NRrGkXSIZLulHRj6rWsRNJaSd8Y\n/L28W9LW1Gual2W8gSuBSyLiROAS4JOJ1zOUpFcBbwBOiIgTgE+lXdHKJK0HXgvk/FOCbgaOj4jN\nwH3ARxOvB5iLDPA54HXA8cBbJb0o7aqGegr4UEQcD/wF8L5M1znvPOCe1ItYxVXATRGxCXgpMJt4\nPQtyjffTwNrB7SOBRxOuZZT3ApdHxFMAEbE/8XpG+TRwfupFjBIRt0TE04O7u4D1KdezyCnAfRHx\ns4h4Evga8MbEa1omIn4ZEXsHt3/HXGjWpV3VcINh4izgC6nXspLBd36vjIhrASLiqYj4v8TLWpBr\nvD8IfErSz5mbwrOYwIY4Djhd0i5Jt+a6vSNpB/BwROxLvZYxvBv4TupFDKwDHl50/xEyjeI8SRuA\nzUCuL7GdHyZyvtztGGC/pGsH2zufl3RE6kXNS/YKS0k7gaMXv4m5T+RFwGuA8yLiBklvBq5h7lv+\nqRuxzouZ+/M7KiJOlXQy8HXg2OmvctV1XsjBf37JXtkx6vMeEd8eHHMR8GREXJdgicWT9Ezgeua+\nhn6Xej1LSXo98KuI2DvYesznlUYHWwOcBLwvIvZI+kfgI8xt5SaX5XXekh6LiCMX3f9tRKwd9ZgU\nJN0EXBERM4P79wNbI+I3aVd2gKSXALcAjzP3RbKeuW2oUyLi1ynXNoykdwLnAGdExBOJlwOApFOB\nSyNi++D+R4CIiCvSrmw5SWuAfwW+ExFXpV7PMJL+Hvgb5vbojwCeBXwrIt6RdGFLSDoa+FFEHDu4\n/wrggoh4Q9qVzcl12+RRSdsAJJ0J3Jt4PSu5ATgDQNJxwGE5hRsgIn4cEc+LiGMj4hjmvuU/MdNw\nb2fuW+kduYR74HbghZJeIOlw4C1ArldIXAPck2u4ASLiwoh4/iCKbwG+n1u4ASLiV8DDg69tgDPJ\n6AnWXH8w1TnAZyQdCvweeE/i9azkWuAaSfuAJ4Ds/gIOEeT7bepngcOBnZIAdkXE36VdEkTEHyW9\nn7mrYQ4BvhgR2Vx1ME/SacDbgX2S7mLuc31hRHw37cqKdi7wVUmHAT8F3pV4PQuy3DYxM7PRct02\nMTOzERxvM7MCOd5mZgVyvM3MCuR4m5kVyPE2MyuQ421mViDH28ysQP8Pa5Xt5v4bpc0AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb861bb3550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "% matplotlib inline\n",
    "\n",
    "# Create random data: we use the randn function\n",
    "X = np.random.randn(3,2)\n",
    "y = np.random.randn(3)\n",
    "\n",
    "# Solve least squares problem minimize \\|X\\beta-y\\|^2\n",
    "# the index 0 says that we get the first component of the solution\n",
    "# (the function lstsq give more output than just the beta vector)\n",
    "beta = la.lstsq(X,y)[0]\n",
    "\n",
    "# Create function and plot the contours\n",
    "def f(a,b):\n",
    "    return sum((a*X[:,0]+b*X[:,1]-y)**2)\n",
    "\n",
    "# Find the \"right\" boundaries around the minimum\n",
    "xx = np.linspace(beta[0]-8,beta[0]+8,100)\n",
    "yy = np.linspace(beta[1]-8,beta[1]+8,100)\n",
    "\n",
    "# A mesh grid gives a pair of matrix XX and YY,\n",
    "# such that for each (i,j), the pairs (XX[i,j],YY[i,j])\n",
    "# cover all the points on the square defined by xx and yy\n",
    "XX, YY = np.meshgrid(xx,yy)\n",
    "\n",
    "# Compute the Z values corresponding to the meshgrid\n",
    "Z = np.zeros(XX.shape)\n",
    "for i in range(Z.shape[0]):\n",
    "    for j in range(Z.shape[1]):\n",
    "        Z[i,j] = f(XX[i,j],YY[i,j])\n",
    "\n",
    "# Before plotting the contour map, choose a nice colormap\n",
    "cmap = plt.cm.get_cmap(\"coolwarm\")\n",
    "plt.contourf(XX,YY,Z, cmap = cmap)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [py27]",
   "language": "python",
   "name": "Python [py27]"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}