CSS 420 Lecture 3 Convexity and Hessians.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lecture 3\n",
    "===========================================================================\n",
    "\n",
    "Convexity and Hessians\n",
    "------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "init_printing(use_latex='mathjax')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets find the minima of a few functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = Symbol('x', real=True)\n",
    "functions = (('a', 3*x**2),\n",
    "            ('b', 2*x),\n",
    "            ('c', -5*x**2),\n",
    "            ('d', 2*x**2 - x**3))\n",
    "for (c, f) in functions:\n",
    "    derivative = diff(f, x, 2)\n",
    "    positivex = solve(derivative>0, x)\n",
    "    print c, str(derivative), positivex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also evaluate symbolic expressions of multiple variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = var('x1 x2')\n",
    "f = 2*x1**3 - 3*x1*x2 + 2*x2**3\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example f(1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f.subs([(x1,1),(x2,2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And take the partial differential w.r.t. each variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = [f.diff(m, 1) for m in [x1,x2]]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second order differential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df2 = [f.diff(m, 2) for m in [x1,x2]]\n",
    "df2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Partial differential w.r.t. one variable then the next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "varx = [x1,x2]\n",
    "dfx1x2 = [df[i].diff(varx[1-i], 1) for i in range(0,2)]\n",
    "dfx1x2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember the Hessian is the matrix of repeated partial differentials:\n",
    "\n",
    "$$\\begin{bmatrix} \\frac{\\partial^2 f}{\\partial x_{1}^2} & \\frac{\\partial^2 f}{\\partial x_{1} \\partial x_{2}} \\\\ \\frac{\\partial^2 f}{\\partial x_{2} \\partial x_{1}} & \\frac{\\partial^2 f}{\\partial x_{2}^2} \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H = hessian(f, [x1,x2])\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can still evaluate these expressions at a specific point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H.subs([(x1,1),(x2,2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = var('x1 x2')\n",
    "f2 = 2*x1**2 + 2*x1*x2 + 1.5*x2**2 + 7*x1 + 8*x2 + 24\n",
    "H2 = hessian(f2, [x1, x2])\n",
    "H2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Positive definiteness is the matrix analogue to positiveness in scalars. Observe these definitions:\n",
    "\n",
    "$x^THx \\geq 0 \\forall x \\iff H~\\text{positive semidefinite}$ similar to $h \\geq 0 \\iff h~\\text{positive}$\n",
    "\n",
    "$x^THx > 0 \\forall x \\neq 0 \\iff H~\\text{positive definite}$ similar to $h > 0 \\iff h~\\text{strictly positive}$\n",
    "\n",
    "$x^THx < 0 \\forall x \\neq 0 \\iff H~\\text{negative definite}$ similar to $h < 0 \\iff h~\\text{strictly negative}$\n",
    "\n",
    "$x^THx \\leq 0 \\forall x \\iff H~\\text{positive semidefinite}$ similar to $h \\leq 0 \\iff h~\\text{negative}$\n",
    "\n",
    "There are several ways to determine whether a matrix is positive definite.\n",
    "\n",
    "The first involves evaluating the principal minors. These are the determinants of matrices obtained by removing a row (i) and a column (j), the principal minors are the cases where i = j."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def find_minor(a, n):\n",
    "    a.row_del(n)\n",
    "    a.col_del(n)\n",
    "    return a\n",
    "\n",
    "def det_minors(b):\n",
    "    msize = np.shape(b)[0]\n",
    "    tmvals = np.arange(msize)\n",
    "    mtmp = b[:,:]\n",
    "    for i in range(0,msize):\n",
    "        minor = find_minor(b,i)\n",
    "        b = mtmp[:,:]\n",
    "        tmvals[i] = minor.det()\n",
    "    mvals = tmvals.tolist()\n",
    "    return mvals\n",
    "\n",
    "tmp = H2[:,:]\n",
    "det_minors(tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def posdef_minor(H):\n",
    "    minors = det_minors(H)\n",
    "    return all(minors[i]>=0 for i in range(0,np.size(minors)))\n",
    "\n",
    "tmp = H2[:,:]\n",
    "posdef_minor(tmp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second method is to check if all the eigenvalues are positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "EV = H2.eigenvals()\n",
    "EVs = [N(i,3) for i in EV]\n",
    "EVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def posdef_eigen(H):\n",
    "    EV = H2.eigenvals()\n",
    "    return all(EVs[i]>=0 for i in range(0,np.size(EVs)))\n",
    "\n",
    "posdef_eigen(H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another convenient way to check if the matrix is positive definite is to attempt a Cholesky decomposition, if it succeeds the matrix is positive definite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H2.cholesky()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what's the difference between these methods? They will all work, but some take longer than others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timeit('posdef_minors(H2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timeit('posdef_eigen(H2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timeit(H2.cholesky())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Calculating the minors and the eigenvalues require approximately the same amounts of time, however the Cholesky decomposition takes significantly longer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}