Lecture8.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "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 = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#Lecture 8: Optimization\n",
    "\n",
    "## Topics\n",
    "\n",
    "* I. Introduction\n",
    "\n",
    "    * Formulation\n",
    "    * Types of Optimization Problems\n",
    "    * Basic Concepts\n",
    "    * Optimization Examples in Finance\n",
    "    \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* II. Unconstrained Optimization\n",
    "\n",
    "    * Optimality Characterization\n",
    "    * Solution Methods\n",
    "        * Direct search method\n",
    "        * Descent Methods\n",
    "        \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "       \n",
    "\n",
    "* III. Constrained Optimization\n",
    "    * General framework\n",
    "    * Duality\n",
    "    * Optimality conditions\n",
    "        * Why KKT conditions\n",
    "    * Solution Methods\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "# I. Introduction\n",
    "\n",
    "Optimization exists in every field, particularly in finance: most activities in finance revolve around some goals that needed to be optimized, e.g. maximizing profits, minimizing risk, maximizing utility, and minimizing costs etc.\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**: 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**: some choices are restricted, 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": [
    "* If $f(\\bx)$ is linear and $\\Psi$ is a polyhedron, then **Linear Programming**, otherwise, **Nonlinear Programming**. (Polyhedron is defined as intersection of finite number of halfspaces and hyperplanes.)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* If $f(\\bx)$ is quadratic and $\\Psi$ is a polyhedron, then **Quadratic Programming**"
   ]
  },
  {
   "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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Finding the global minimum is considerably harder than a local minimum --- imaging 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, any local minimum is also globally optimal (homework).\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 normally ill-advised --- doing anything through square or higher power when you can avoid it is generally a bad idea --- it makes the problem harder to solve, think ill-conditioning as a starter).\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 **sufficient 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",
    "\n",
    "\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",
    "$$\n",
    "f(\\bx^* + \\bs\\delta x) = f(\\bx^*) + {\\bs\\nabla} f(\\bx^*)^T {\\bs\\delta x} + \\frac{1}{2} {\\bs\\delta x}^T{\\bs\\nabla^2} f(\\bx^*)  {\\bs\\delta x} + {\\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, $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": [
    "# II. 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 differentiable). 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 dimension. Requires only function evaluations.\n",
    "\n",
    "* Quoting M. Wright: \"A direct search 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": {
    "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: [ 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:\", res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Advantages: simple, only function evaluation 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",
    "$$\n",
    "f(\\bx^k + t \\bs\\delta x) \\approx f(\\bx^k) + t [ {\\bs\\nabla} f(\\bx^k)^T {\\bs\\delta x} ]                        \n",
    "$$\n",
    "\n",
    "* This means choosing **gradient** direction\n",
    "\n",
    "$$\n",
    "\\bs\\delta x = - {\\bs\\nabla} f(\\bx^k)\n",
    "$$\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",
    "$$\n",
    "f(\\bx^k + \\bs\\delta x) \\approx f(\\bx^k) + {\\bs\\nabla} f(\\bx^k)^T {\\bs\\delta x} + \\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",
    "$$\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",
    "* 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",
    "* 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# III. Constrained Optimization\n",
    "\n",
    "* Now add  constraints.\n",
    "\n",
    "* Constrained 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 differentiable). \n",
    "\n",
    "* The goal is to find a local minimum satisfying the constraints.\n",
    "\n",
    "* And we will denote the overall domain of definition as\n",
    "$\\mathcal{D} = DOM(f)\\cap DOM(h) \\cap DOM(g)$."
   ]
  },
  {
   "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 for $\\forall \\bld \\succeq 0$ and any $\\bmu$,\n",
    "\n",
    "$$\n",
    "\\df(\\bmu, \\bld) \\leq f(\\bx^*).\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 straightforward (homework) 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",
    "\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)\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": [
    "\n",
    "\n",
    "* Best way to visualize the conditions 1 and 4 above is through the following example, where\n",
    "\n",
    "$$ f(\\bx^*) = (x_1 -2)^2 + (x_2 + 2)^2, \\;\\;\\; g(\\bx^*) = x_1^2 + x_2^2 - 1$$\n",
    "<center><img src = \"img/KKT_cond_4.PNG\" height=700 width=700></center>\n",
    "\n",
    "  * with $ \\bld > 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Solution Methods\n",
    "\n",
    "* There are two main competing approaches to solve the constrained optimization problems:\n",
    "    * Sequential Quadratic Programming (SQP): BFGS, NLOPT, SNOPT\n",
    "    * Interior Point Method (IPT): IPOPT, KNITRO, LOQO\n",
    "* We are not going into great details, but will expand the SQP method slightly (and will get into IPT method later).\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)} + \\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "# References: \n",
    "S. Boyd and L. Vandenberghe (2004), Convex Optimization, Cambridge University Press.\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}