Lecture6-Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 6: Optimization\n",
    "\n",
    "## Topics\n",
    "\n",
    "* Introduction\n",
    "* Unconstrained Optimization\n",
    "* Constrained Optimization\n",
    "    \n",
    "## References\n",
    "\n",
    "* D. Bindel and J. Goodman: Principles of Scientific Computing, Chapter 6, 2009.\n",
    "* S. Boyd and L. Vandenberghe: Convex Optimization, Cambridge University Press, Chapters 5, 9, 11, 2004.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "lecture = 6\n",
    "\n",
    "import sys\n",
    "sys.path.append(\"lib\")\n",
    "import fmt\n",
    "import sympy as sp\n",
    "from IPython.display import display\n",
    "\n",
    "# assert sp.__version__ == \"0.7.5\", \"Need sympy version 0.7.5 to render properly\"\n",
    "sp.init_printing(use_latex = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "# Introduction\n",
    "\n",
    "* Optimization is an indispensible tool in computational finance\n",
    "* Problems such as: asset allocation, risk management, model calibration, etc all use optimization in one way or another\n",
    "* Optimization itself is a vast subject, more suitably covered in a full semester course\n",
    "* We will cover some of the basic theories and go into a little bit of details in a couple of examples.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "## Formulation\n",
    "\n",
    "An optimization problem usually involves three elements:\n",
    "* **Objective function**: the function or quantity to be optimized: profit, loss, risk, etc\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* **Variables**: numbers of shares in each stock, amount of capital to be invested in each sector, etc\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "* **Constraints**: sometimes the allowable variables are restricted, e.g. total risk cannot exceed certain criteria, total amount of capital is limited, etc\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Mathematically,\n",
    "\n",
    "$$\n",
    "\\renewcommand{bx}{\\boldsymbol x}\n",
    "\\Large{ \\min_{\\bx\\in \\Psi, \\; \\Psi\\subset \\mathbb{R}^n}} f(\\bx)\n",
    "$$\n",
    "\n",
    "where $ f(\\bx): \\mathbb{R}^n \\to \\mathbb{R} $ is a scalar **objective function**, $\\Psi$ is a subset of $\\mathbb{R}^n$ and called the ** feasible region**.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* The minimization and maximization problems are equivalent\n",
    "$$\n",
    "\\min_{\\bx\\in \\Psi} f(\\bx) = -\\max_{\\bx\\in \\Psi} (-f(\\bx))\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Types of Optimization Problems\n",
    "\n",
    "\n",
    "* If $\\Psi = \\mathbb{R}^n$, the problem is **unconstrained**, otherwise, it is **constrained**\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* If $f(\\bx)$ is linear and $\\Psi$ is a polyhedron, then **Linear Programming**, otherwise, **Nonlinear Programming**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* If $f(\\bx)$ is quadratic and $\\Psi$ is a polyhedron, then **Quadratic Programming**\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* If $f(\\bx)$ and $\\Psi$ are convex, then **Convex Optimization** (of which linear and quadratic programming are special cases)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "* If $\\Psi$ contains discrete variables, then **Discrete Optimization**\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* If only integer variables are allowed, then **Integer Programming**, **Mixed Integer Programming** involves problems in which only some of the variables are constrained to be integers.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* If the specifications of $f(\\bx)$ and $\\Psi$ are NOT deterministic, then **Stochastic Programming** \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* Another name you often hear is **Dynamic Programming**, this does not refer to a particular type of optimization problem, rather it is a method for solving an optimization problem by breaking it down to a collection of simpler subproblems, using Bellman's Principle of Optimality.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Basic Concepts\n",
    "\n",
    "* A point $\\bx^*$ is called a **local minimum**, if $\\exists \\epsilon > 0$ \n",
    "$$\n",
    "f(\\bx^*) \\leq f(\\bx), \\forall \\bx\\in\\Psi, s.t. \\|  \\bx  - \\bx^* \\| \\lt \\epsilon.\n",
    "$$\n",
    "\n",
    "* A point $\\bx^*$ is a **global minimum**, if  \n",
    "$$\n",
    "f(\\bx^*) \\leq f(\\bx), \\forall \\bx\\in\\Psi.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "* Finding a global minimum is considerably harder than a local minimum --- imagine the difference in the difficulties of getting to the top of your neighborhood hill top vs getting to the top of Mt. Everest.\n",
    "\n",
    "* With one exception: for convex problems, the local minimum is the also global solution.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* If $f(\\bx)$ is differentiable (i.e. first derivatives exist), then a **necessary condition** for a local minimizer is that $\\bx^*$ is a **critical point**\n",
    "\n",
    "$$\n",
    "\\renewcommand{bs}{\\boldsymbol}{g} (\\bx^*) = {\\bs\\nabla} _\\bx f(\\bx^*)  = 0\n",
    "$$\n",
    "\n",
    "* Clearly, finding the critical points is equivalent to root searching problems you have encountered earlier.\n",
    "\n",
    "* On the other hand, observing that solving $f(\\bx) = \\bs{0}$ is equivalent to\n",
    "\n",
    "$$\n",
    "\\min_{\\bx}\\left[ f(\\bx)^T f(\\bx)\\right],\n",
    "$$\n",
    "\n",
    "(although this is rarely advised --- doing anything through square or higher power is generally a bad idea --- it  makes the problem harder to solve).\n",
    "\n",
    "* But this shows the root searching problem and the optimization problem are closely related.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* Critical point is only a necessary condition, it's not sufficient,\n",
    "\n",
    "<center><img src = \"img/saddle_point_func.png\" height=450 width=450></center>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* If $f(\\bx)$ is twice-differentiable (i.e. second derivatives exist), then a **sufficent condition** for a local minimizer is, in addition to being a critical point,  the **Hessian** at $\\bx^*$ is positive definite\n",
    "\n",
    "$$\n",
    "{\\bs H} (\\bx^*) = {\\bs\\nabla^2} _\\bx f(\\bx^*)  \\succ 0\n",
    "$$\n",
    "\n",
    "*by J:\n",
    "For the scalar case $f(x)=0$ (root finding), we can use\n",
    "1. Bisection\n",
    "2. Secant method\n",
    "3. Newton's method\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\n",
    "* These are all  natural conclusions from what you are already very familiar with (hopefully!): properties of a quadratic function and the Taylor series for a general (twice-differentiable) function\n",
    "$$\n",
    "f(\\bx^* + \\bs\\delta x) = f(\\bx^*) \n",
    "                       + {\\bs\\nabla} f(\\bx^*)^T {\\bs\\delta x} \n",
    "                       + \\frac{1}{2} {\\bs\\delta x}^T{\\bs\\nabla^2} f(\\bx^*)  {\\bs\\delta x}\n",
    "                       + {\\bs O}(||{\\bs\\delta x}||^3)\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Optimization Examples in Finance\n",
    "\n",
    "### 1. Portfolio Optimization\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\min_{\\bx} & &  \\frac{1}{2} \\lambda\\; \\bx^T \\Sigma \\bx - \\mu^T \\bx \n",
    "\\\\\n",
    "s.t. & & \\Sigma x_i = 1\n",
    "\\end{array}\n",
    "$$\n",
    "where $\\lambda$ is the risk-aversion coefficient, $\\mu$ is the expected asset return vector and $\\Sigma$ is the covariance matrix. This is a quadratic programming problem.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "### 2. (Static) Asset-Liability Management\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\min_{\\bx} & &  \\Sigma_j x_j P_j\n",
    "\\\\\n",
    "s.t. & & \\Sigma_j x_j C_j(t) \\geq L(t) \\;\\;  \\forall t\n",
    "\\\\\n",
    "     & & x_j \\geq 0 \\;\\; \\forall j\n",
    "\\end{array}\n",
    "$$\n",
    "where $x_j, P_j, C_j(t)$ are the amount, price and cashflow at time $t$ of asset $j$. $L(t)$ is the liability payment at time $t$. This is a linear programming problem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "### 3. Volatility Surface fitting\n",
    "\n",
    "$$\n",
    "\\min_{\\sigma(S,t)} {\\Large\\Sigma_j^n} (C(\\sigma(S,t),K_j,T_j) - C_j)^2 \n",
    "$$\n",
    "where: $\\sigma(S,t) > 0$ is the volatility value at the surface point $(S,t)$; $C(\\sigma(S,t),K_j,T_j)$ is the standard Black-Scholes formule for European call options; and $C_j$ is the market quoted price of the option. This a non-linear optimization problem.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Unconstrained Optimization\n",
    "\n",
    "* Unconstrained means:\n",
    "\n",
    "$$\n",
    "\\min_{\\bx\\in  \\mathbb{R}^n} f(\\bx)\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "* Will be focusing on smooth functions (typically at least twice differentable). The main purpose here is to show the various types of algorithms for solving the unconstrained optimization problems.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "## Optimality Characterization\n",
    "\n",
    "* The **necessary condition**: the solution must be a **critical point**:\n",
    "\n",
    "$$\n",
    "{\\bs g} (\\bx^*) = {\\bs\\nabla} _\\bx f(\\bx^*)  = 0\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "* The **sufficient condition**: the **Hessian** at the optimal point must be positive definite,\n",
    "\n",
    "$$\n",
    "{\\bs H} (\\bx^*) = {\\bs\\nabla^2} _\\bx f(\\bx^*)  \\succ 0\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Solution Methods\n",
    "\n",
    "* In practice, optimization problems are often solved using an **iterative algorithm**, which searches for a sequence of points,\n",
    "\n",
    "$$\n",
    "\\bx^0, \\bx^1, \\bx^2, \\cdots, \\bx^n, \\cdots\n",
    "$$\n",
    "with \n",
    "$$\n",
    "f(\\bx^{k+1}) < f(\\bx^k).\n",
    "$$\n",
    "\n",
    "* The algorithm typically stops at $||{\\bs\\nabla} f(\\bx) || < \\epsilon $ for some $\\epsilon$ small.\n",
    "\n",
    "* No guarantee for finding the global minimum.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "### 1. Direct Search Method\n",
    "\n",
    "* Similar in spirit to the Bisection method in one dimemsion. Requires only function evaluations.\n",
    "\n",
    "* Quoting M. Wright: \"A direct serach method does not 'in its heart' develop an approximate gradient\".\n",
    "\n",
    "\n",
    "* Representative: Nelder-Mead Method (or Simplex Search method)\n",
    "    * Searches through the simplex vertices (polytope of N+1 vertices in N dimensions)\n",
    "    * Techniques: reflection - expansion - contraction - reduction\n",
    "\n",
    "* Scipy example:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.000000\n",
      "         Iterations: 650\n",
      "         Function evaluations: 1031\n",
      "The solution from Nelder-Mead:\n",
      "[ 1.  1.  1.  1.  1.  1.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def rosen(x):\n",
    "     return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)\n",
    "    \n",
    "x0 = np.array([1.3, 0.7, 0.8, 2.2, 1.2, 2.1])\n",
    "res = minimize(rosen, x0, method='nelder-mead',\n",
    "                options={'xtol': 1e-8, 'disp': True})\n",
    "\n",
    "print \"The solution from Nelder-Mead:\"\n",
    "print  (res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Advantages: simple, only function evaulation needed.\n",
    "\n",
    "* Deficiencies: slow, may fail to converge in higher dimensions\n",
    "\n",
    "* Suffers from the \"curse of dimensionality\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### 2. Descent Methods\n",
    "\n",
    "** Algorithm for general descent method **\n",
    "1. Given a starting point $ \\bx^0$\n",
    "2. Repeat\n",
    "    1. Determine a descent direction $\\delta \\bx$;\n",
    "    2. Line search. Choose a step size $t > 0$;\n",
    "    3. Update. $\\bx^{k+1} = \\bx^k + t \\delta \\bx $\n",
    "\n",
    "3. Until stopping criterion is satisfied\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "The algorithm alternates between two main decisions: determine a descent direction  $\\delta \\bx$ and choose a step size $t$.\n",
    "* Different ways of choosing the descent direction giving rise to different descent method and convergence rate\n",
    "\n",
    "* The line search method falls into two categories: exact line search and backtracking line search. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "#### 2.1 Steepest Descent Method\n",
    "\n",
    "* If the objective function is differentiable, we have\n",
    "$$\n",
    "f(\\bx^k + t \\bs\\delta x) \\approx f(\\bx^k) \n",
    "                       +t [ {\\bs\\nabla} f(\\bx^k)^T {\\bs\\delta x} ]                        \n",
    "$$\n",
    "\n",
    "* This means choosing **gradient** direction\n",
    "$$\n",
    "\\bs\\delta x = - {\\bs\\nabla} f(\\bx^k)\n",
    "$$\n",
    "will lead to the **steepest** descent at points sufficiently close to $\\bx^k$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* The line search step size can be done in one-dimensional minimization:\n",
    "$$\n",
    "t^k = \\min_t f(\\bx^k + t \\bs\\delta x) \\triangleq \\min_t \\phi(t).\n",
    "$$\n",
    "\n",
    "\n",
    "* Exact line search (choosing the minimizing $t$ above) leads to zig-zag path towards the minimum: which means slow convergence\n",
    "$$\n",
    "\\phi'(t) = 0 =  [{\\bs\\nabla} f(\\bx^k + t \\bs\\delta x)]^T \\bs\\delta x ,\n",
    "$$\n",
    "(notice the two consecutive search directions will be perpendicular to each other, we've met this problem before, what's the strategy?).\n",
    "\n",
    "* Convergence: the steepest descent method converges linearly, and it will behave badly if the condition number of the **Hessian** (the second order derivative matrix) is large.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "#### 2.2 Newton's Method\n",
    "\n",
    "* If the objective function is twice differentiable, we have with more accuracy\n",
    "$$\n",
    "f(\\bx^k + \\bs\\delta x) \\approx f(\\bx^k) \n",
    "                       + {\\bs\\nabla} f(\\bx^k)^T {\\bs\\delta x}                             \n",
    "                       + \\frac{1}{2} {\\bs\\delta x}^T{\\bs\\nabla^2} f(\\bx^k)  {\\bs\\delta x}\n",
    "$$\n",
    "\n",
    "* RHS is a quadratic function in ${\\bs\\delta x}$, so the minimum is achieved at\n",
    "$$\n",
    "\\bs\\delta x = - [{\\bs\\nabla^2} f(\\bx^k)]^{-1} {\\bs\\nabla} f(\\bx^k).\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "* Convergence of the Newton's method is rapid in general, and quadratic once entering into the pure Newton phase.\n",
    "\n",
    "\n",
    "* Disadvantages of Newton's method:\n",
    "    * The cost of computing and storing the Hessian can be very high, if not outright prohibitive\n",
    "    * The cost of solving the set of linear equation at the Newton step\n",
    "\n",
    "\n",
    "* There are various ways to compute an approximation of the Hessian to substantially reduce the cost of computing the Newton step. This leads to a family of algorithms called **Quasi-Newton methods**.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Constrained Optimization\n",
    "\n",
    "* Now add  constraints.\n",
    "\n",
    "* Constrainted problems are much harder: even a seemingly simple Integer Programming problem is **NP-Complete** (i.e. can't be solved in polynomial time).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "## General Framework\n",
    "\n",
    "\n",
    "* Constrained Optimization Problem\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\min_{\\bx\\in \\mathbb{R}^n } &  f(\\bx) \n",
    "\\\\\n",
    "s.t. &\\bs{h}(\\bx)  = \\bs{0}\n",
    "\\\\\n",
    "  &\\bs{g}(\\bx) \\leq \\bs{0}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "\n",
    "* Will be focusing on smooth functions (typically at least twice differentable). \n",
    "\n",
    "* The goal is to find a local minimum satisfying the constriants.\n",
    "\n",
    "* And we will denote the feasible region as domain \n",
    "$\\mathcal{D} = \\{\\bx \\in \\mathbb{R}^n | \\; \\bs{h}(\\bx)  = \\bs{0}, \\; \\bs{g}(\\bx) \\leq \\bs{0}\\}$.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "## Duality\n",
    "\n",
    "* Define the **Lagrangian ** as,\n",
    "$$\n",
    "\\renewcommand{ml}{\\mathbb{\\mathcal L}}\n",
    "\\renewcommand{bmu}{\\boldsymbol{ \\mu}}\n",
    "\\renewcommand{bld}{\\boldsymbol{ \\lambda}}\n",
    "\\ml(\\bx, \\bmu, \\bld) = f(\\bx) + \\bmu^T \\bs{h}(\\bx) + \\bld^T \\bs{g}(\\bx).\n",
    "$$\n",
    "\n",
    "   Here the vectors $\\bmu, \\bld$ are called the **dual variables** or **Lagrange Multipliers**.\n",
    "\n",
    "\n",
    "* Further, define the **Lagrangian dual function** as,\n",
    "\n",
    "$$\n",
    "\\renewcommand{mD}{\\mathbb{\\mathcal D}}\n",
    "\\renewcommand{df}{\\hat{f}}\n",
    "\\df(\\bmu, \\bld) = \\inf_{\\bx\\in \\mD}\\; \\ml(\\bx, \\bmu, \\bld) =\\inf_{\\bx\\in\\mD} \\left( f(\\bx) + \\bmu^T \\bs{h}(\\bx) + \\bld^T \\bs{g}(\\bx) \\right).\n",
    "$$\n",
    "\n",
    "\n",
    "* The dual function takes the pointwise infimum of a family of affine functions of $(\\bmu, \\bld)$, it is a concave function.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* If $\\bx^*$ is a solution to the original optimization problem (**Primal Problem**), then $\\forall \\bld \\succeq 0$ and any $\\bmu$,\n",
    "\n",
    "$$\n",
    "\\df(\\bmu, \\bld) \\leq f(\\bx^*).\n",
    "$$\n",
    "     \n",
    "   \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "* Which leads to the optimization problem (**Dual Problem**)\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\max_{(\\bmu, \\bld)} &  \\df(\\bmu, \\bld) \n",
    "\\\\\n",
    "s.t. &\\bld \\succeq 0\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "* If $(\\bmu^*, \\bld^*)$ is a solution to the  **Dual Problem**, it's clear that **the weak duality** holds,\n",
    "\n",
    "$$\n",
    "\\df(\\bmu^*, \\bld^*) \\leq f(\\bx^*).\n",
    "$$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "* However, if the **Primal problem** satisfies certain **constraint qualifications** (such as, convexity or Slater's condition), then the **strong duality** holds,\n",
    "$$\n",
    "\\df(\\bmu^*, \\bld^*) = f(\\bx^*),\n",
    "$$\n",
    "\n",
    "    which implies the primal and the dual problems are equivalent.\n",
    "\n",
    "\n",
    "* In the **Maximum Entropy Method** lecture later, this will be explored in greater details.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "## Optimality Conditions\n",
    "\n",
    "\n",
    "\n",
    "* From the duality principle, \n",
    "\n",
    "    1. $ \\bs{\\nabla}_{\\bx}\\ml(\\bx^*, \\bmu^*, \\bld^*) = \\bs{0} $, stationality \n",
    "    2. $ \\bs{h}(\\bx^*)  = \\bs{0}; \\bs{g}(\\bx^*) \\leq \\bs{0} $, feasibility  \n",
    "    3. $ \\bld^* \\succeq \\bs{0} $, dual feasibility  (component-wise) \n",
    "    4. $ \\bld^* \\circ \\bs{g}(\\bx^*) = \\bs{0} $, complementary slackness (component-wise multiplication: convolution, *by J)\n",
    "    5. $ \\bs{\\nabla}^2\\ml(\\bx^*, \\bmu^*, \\bld^*) \\succ 0 $, Hessian positive definite constraints\n",
    "\n",
    "\n",
    "* These conditions are called **KKT conditions** (Karush-Kuhn-Tucker) --- the necessary and sufficient conditions for $\\bx^*$ to be a local minimizer.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "### Exploring the KKT conditions further\n",
    "\n",
    "* For the **unconstrained** case, the conditions 2, 3 and 4 drop out, what's left are:\n",
    "\n",
    "    1. $ \\bs{\\nabla}_{\\bx}f(\\bx^*) = \\bs{0} $, stationality \n",
    "    5. $ \\bs{\\nabla}^2f(\\bx^*) \\succ 0 $, Hessian positive definite constraints\n",
    "\n",
    "\n",
    "* which we are familiar with---the necessary and sufficient condition for $\\bx^*$ to be a (local) minimizer.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* For the **equality constrained** case, let's take a look at the Lagrange function\n",
    "\n",
    "$$\n",
    "\\ml(\\bx, \\bmu) = f(\\bx) + \\bmu^T \\bs{h}(\\bx)\n",
    "$$\n",
    "\n",
    "\n",
    "* If we simply consider this as an unconstrained problem with $(\\bx, \\bmu)$ as the new unknown vector and apply the two conditions in previous slide\n",
    "\n",
    "    1. $ \\bs{\\nabla}_{\\bx}\\ml(\\bx^*, \\bmu^*) = \\bs{0} $, stationality \n",
    "    2. $ \\bs{\\nabla}_{\\bmu}\\ml(\\bx^*, \\bmu^*) = \\bs{0} $, stationality \n",
    "    5. $ \\bs{\\nabla}^2_{\\bx\\bx}\\ml(\\bx^*, \\bmu^*) \\succ 0 $, Hessian positive definite constraints\n",
    "\n",
    "\n",
    "* Notice the second condition above is simply $\\bs{h}(\\bx^*)  = \\bs{0}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* For the **inequality(only) constrained** problems, the answer can be analyzed in two scenarios based on how the local minimizer $\\bx^*$ is situated in the feasible region.\n",
    "\n",
    "\n",
    "* **Scenario 1**: $\\bs{g}(\\bx^*) < \\bs{0}$,  in this case the point $\\bx^*$ is an interior point of the feasible domain and the constraints is called **inactive** and the case simply reduces to the unconstrained case:\n",
    "\n",
    "    1. $ \\bs{\\nabla}_{\\bx}f(\\bx^*) = \\bs{0} $, stationality \n",
    "    5. $ \\bs{\\nabla}^2_{\\bx\\bx}f(\\bx^*) \\succ 0 $, Hessian positive definite constraints\n",
    "    2. and of course $ \\bs{g}(\\bx^*) < \\bs{0} $,\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "* **Scenario 2**: $\\bs{g}(\\bx^*) = \\bs{0}$,  in this case the point $\\bx^*$ is right on the boundary of the feasible domain and the constraints is called **active** and the case reduces to the equality constrained case:\n",
    "\n",
    "    1. $ \\bs{\\nabla}_{\\bx}f(\\bx^*) + \\bld^T \\bs{\\nabla}_{\\bx}\\bs{g}(\\bx^*) = \\bs{0} $, stationality \n",
    "    5. $ \\bs{\\nabla}^2_{\\bx\\bx}f(\\bx^*) \\succ 0 $, Hessian positive definite constraints\n",
    "    2. $ \\bs{g}(\\bx^*) = \\bs{0} $\n",
    "    3. and $ \\bld > 0 $ component wise\n",
    "\n",
    "\n",
    "* Note the last condition ensures that the decreasing direction of $f(\\bx)$ is  strictly pointing outwards of the feasible region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Solution Methods\n",
    "\n",
    "* The solution method for the constrained optimization problems are generally an extension of those for unconstrained version. We are not going into great details, but will use the Sequential Quadratic Programming (**SQP**) method as an example.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* The SQP is a natural extension of the Newton's Descent method introduced earlier for the unconstrained case: recall that the Newton's method does a quadratic optimization in each iteration step. For SQP, at each step, the original optimization problem is approximated with\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\min_{\\bx\\in \\mathbb{R}^n} & {\\bs\\nabla} f(\\bx^k)^T {(\\bx - \\bx^k)}                             \n",
    "                       + \\frac{1}{2} {(\\bx - \\bx^k)}^T{\\bs\\nabla^2} f(\\bx^k)  {(\\bx - \\bx^k)}\n",
    "\\\\\n",
    "s.t. & {\\bs\\nabla} \\bs{h}(\\bx^k)^T {(\\bx - \\bx^k)} + \\bs{h}(\\bx^k)  = \\bs{0}\n",
    "\\\\\n",
    "  & {\\bs\\nabla} \\bs{g}(\\bx^k)^T {(\\bx - \\bx^k)} + \\bs{g}(\\bx^k)  \\leq \\bs{0}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "* This is a constrained quadratic programming problem, which is \"slightly\" easier to deal with than the original problem.\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}