3-Otimização.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Otimização\n",
    "\n",
    "Nesta aula estudaremos alguns métodos de minimização de funções de uma variável. Dada uma função $f(x)$, nosso objetivo é encontrar $x$ que produza o menor valor possível para $y = f(x)$. Este tipo de problema aparece em diversas áreas da matemática e engenharia, ainda que normalmente na versão multivariada que veremos em outra ocasião.\n",
    "\n",
    "Considere como exemplo a função $f(x) = -x e^{-x}$. Vemos no gráfico que ela claramente possui um mínimo próximo de $x = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return -x * np.exp(-x)\n",
    "\n",
    "X = np.linspace(0, 10, 100)\n",
    "Y = f(X)\n",
    "plt.plot(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos mostrar analiticamente que o mínimo está de fato no ponto $x = 1$, basta derivar e igualar a zero:\n",
    "\n",
    "$$\n",
    "\\begin{align}f'(x) &= -e^{-x} + x e^{-x} = 0 \\\\ \n",
    "                   & \\Rightarrow x = 1\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Dependendo da função, o cálculo da derivada pode ser bastante complicado ou a etapa de isolar $x$ quando a derivada é igualada à zero pode ser inviável de se fazer analiticamente. Para isto, vamos mostrar alguns métodos numéricos que podem ser aplicados sem utilizar estes passos complicados."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Busca ternária\n",
    "\n",
    "Vamos supor que não saibamos que o mínimo esteja exatamente no ponto $x = 1$, mas temos uma idéia razoável que ele esteja dentro do intervalo $x\\in[0,10]$ (como se pode ver claramente pelo gráfico).\n",
    "\n",
    "Iniciamos nossa busca fixando estes dois pontos avaliando dois pontos dentro deste intervalo. A estratégia lembra um pouco o método da bisseção, com a diferença que guardamos sempre os valores dos extremos dos intervalos junto com um valor no centro menor que ambos os extremos (se o valor do centro for maior que algum dos extremos, significa que a função possui mais de um mínimo e não pode ser tratada por este método).\n",
    "\n",
    "A figura mostra o primeiro passo do algoritmo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Prepara pontos no interior do intervalo\n",
    "x1 = 0\n",
    "x2 = 10\n",
    "xa = x1 + (x2 - x1)/3\n",
    "xb = x2 - (x2 - x1)/3\n",
    "\n",
    "# Gráfico\n",
    "plt.plot(X, Y)\n",
    "plt.xticks([x1, xa, xb, x2], ['$x_1$','$x_a$', '$x_b$', '$x_2$'])\n",
    "xmin, xmax, ymin, ymax = plt.axis()\n",
    "for x in [x1, xa, xb, x2]:\n",
    "    plt.plot([x, x], [ymin, ymax], 'k--')\n",
    "    plt.plot([x], [f(x)], 'ko')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que dos dois novos pontos, $x_a$ possui um valor menor. Assim, atualizamos nosso intervalo para que $x_b$\n",
    "seja promovido ao extremo direito e $x_a$ continue dentro do intervalo. Desta forma, no passo seguinte do algoritmo temos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Prepara pontos no interior do intervalo\n",
    "x1 = 0\n",
    "x2 = xb                # x2 vira o antigo xb!\n",
    "xa = x1 + (x2 - x1)/3  # atualizamos os novos pontos dentro do intervalo\n",
    "xb = x2 - (x2 - x1)/3\n",
    "\n",
    "# Gráfico\n",
    "plt.plot(X, Y)\n",
    "plt.xticks([x1, xa, xb, x2], ['$x_1$','$x_a$', '$x_b$', '$x_2$'])\n",
    "xmin, xmax, ymin, ymax = plt.axis()\n",
    "for i, x in enumerate([x1, xa, xb, x2]):\n",
    "    plt.plot([x, x], [ymin, ymax], ('k--' if i in [1, 2] else 'k-'))\n",
    "    plt.plot([x], [f(x)], 'ko')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O algoritmo segue assim por mais alguns passos. Devemos sempre reduzir o intervalo escolhendo o *maior* ponto entre os dois pontos escolhidos em cada etapa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "for i in range(7):\n",
    "    xa = x1 + (x2 - x1)/3\n",
    "    xb = x2 - (x2 - x1)/3\n",
    "\n",
    "    # Primeiro fazemos o gráfico\n",
    "    plt.plot(X, Y)\n",
    "    plt.xticks([x1, xa, xb, x2], ['$x_1$','$x_a$', '$x_b$', '$x_2$'])\n",
    "    xmin, xmax, ymin, ymax = plt.axis()\n",
    "    for i, x in enumerate([x1, xa, xb, x2]):\n",
    "        plt.plot([x, x], [ymin, ymax], ('k--' if i in [1, 2] else 'k-'))\n",
    "        plt.plot([x], [f(x)], 'ko')\n",
    "    plt.show()\n",
    "        \n",
    "    # Agora atualizamos o resultado\n",
    "    fa = f(xa)\n",
    "    fb = f(xb)\n",
    "    if fa > fb:\n",
    "        x1 = xa\n",
    "    else:\n",
    "        x2 = xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Veja como o intervalo vai gradualmente encolhendo em torno do mínimo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Implementação\n",
    "\n",
    "Vamos agora implementar a busca ternária mais formalmente como uma função chamada `ternaria()` que recebe uma função objetivo e um intervalo e retorna o mínimo da função. A estratégia é continuar a iteração até que o intervalo seja suficientemente pequeno."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def ternaria(f, x1, x2, xtol=1e-6):\n",
    "    \"\"\"\n",
    "    Realiza a busca ternária pelo mínimo da função f(x) data dentro do intervalo [x1, x2].\n",
    "    \n",
    "    Args:\n",
    "        f: função objetivo f(x).\n",
    "        x1, x2: início e fim do intervalo de busca.\n",
    "        xtol: tolerância no intervalo de busca.\n",
    "        \n",
    "    Returns:\n",
    "        Retorna o mínimo quando o tamanho do intervalo de busca for menor que xtol.\n",
    "    \"\"\"\n",
    "    \n",
    "    while True:\n",
    "        xa, xb = x1 + (x2 - x1)/3, x1 + (x2 - x1) * 2 / 3\n",
    "        fa, fb = f(xa), f(xb)\n",
    "        \n",
    "        if fa > fb:\n",
    "            x1 = xa\n",
    "        else:\n",
    "            x2 = xb\n",
    "        \n",
    "        if abs(x1 - x2) < xtol:\n",
    "            break\n",
    "    \n",
    "    return xa if fa < fb else xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora testamos com a nossa função exemplo e algumas outras funções conhecidas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ternaria(f, 0, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ternaria(np.cos, 0, 2 * np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Busca pela razão áurea\n",
    "\n",
    "Na busca ternária, dividimos o intervalo em três partes iguais. Ainda que isto pareça intuitivo, não é a melhor opção. Um problema da busca ternária é que a cada iteração precisamos avaliar a função em dois pontos. Se fizermos a escolha adequada, é possível avaliar a função em apenas um ponto e reciclar o ponto anterior. Para que isto ocorra, é necessário dividir o intervalo em partes proporcionais à razão áurea $phi$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "phi = (1 + np.sqrt(5)) / 2\n",
    "x1, x2 = 0, 10\n",
    "\n",
    "for i in range(7):\n",
    "    xa = x2 - (x2 - x1)/phi\n",
    "    xb = x1 + (x2 - x1)/phi\n",
    "    \n",
    "    # Primeiro fazemos o gráfico\n",
    "    plt.plot(X, Y)\n",
    "    plt.xticks([x1, xa, xb, x2], ['$x_1$','$x_a$', '$x_b$', '$x_2$'])\n",
    "    xmin, xmax, ymin, ymax = plt.axis()\n",
    "    for i, x in enumerate([x1, xa, xb, x2]):\n",
    "        plt.plot([x, x], [ymin, ymax], ('k--' if i in [1, 2] else 'k-'))\n",
    "        plt.plot([x], [f(x)], 'ko')\n",
    "    plt.show()\n",
    "        \n",
    "    # Agora atualizamos o resultado\n",
    "    fa = f(xa)\n",
    "    fb = f(xb)\n",
    "    if fa > fb:\n",
    "        x1 = xa\n",
    "    else:\n",
    "        x2 = xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Note que uma das linhas em $x_a$ ou $x_b$ sempre se alinha com um dos $x_a$ e $x_b$ da iteração posterior. Isto é proposital e permite que avaliemos a função apenas em um dos pontos interiores ao novo intervalo de iteração. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "phi = (1 + np.sqrt(5)) / 2\n",
    "\n",
    "def razao_aurea(f, x1, x2, xtol=1e-6):\n",
    "    \"\"\"\n",
    "    Retorna o mínimo da função f(x) data dentro do intervalo [x1, x2] utilizando a busca pela razão áurea.\n",
    "    \n",
    "    Args:\n",
    "        f: função objetivo f(x).\n",
    "        x1, x2: início e fim do intervalo de busca.\n",
    "        xtol: tolerância no intervalo de busca.\n",
    "        \n",
    "    Returns:\n",
    "        Retorna o valor de x que minimiza f(x) quando o tamanho do intervalo de busca for menor que xtol.\n",
    "    \"\"\"\n",
    "    \n",
    "    xa = x2 - (x2 - x1)/phi\n",
    "    xb = x1 + (x2 - x1)/phi\n",
    "    fa, fb = f(xa), f(xb)\n",
    "    update = None\n",
    "    \n",
    "    while True:\n",
    "        if fa > fb:\n",
    "            x1 = xa\n",
    "            xa, fa = xb, fb\n",
    "            xb = x1 + (x2 - x1)/phi\n",
    "            fb = f(xb)\n",
    "        else:\n",
    "            x2 = xb\n",
    "            xb, fb = xa, fa\n",
    "            xa = x2 - (x2 - x1)/phi\n",
    "            fa = f(xa)\n",
    "        \n",
    "        if abs(x1 - x2) < xtol:\n",
    "            break\n",
    "            \n",
    "    return xa if fa < fb else xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora testamos alguns resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "razao_aurea(f, 0, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "razao_aurea(np.cos, 0, 2*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Finalmente, comparamos a performance de cada método."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 ternaria(f, 0, 10)\n",
    "%timeit -n1000 razao_aurea(f, 0, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 ternaria(np.cos, 0, 2*np.pi)\n",
    "%timeit -n1000 razao_aurea(np.cos, 0, 2*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que nos dois casos a razão áurea oferece quase o dobro da performance da busca ternária. Isto ocorre porque reaproveitamos uma avaliação da função por iteração."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Descendo o gradiente\n",
    "\n",
    "Uma outra estratégia de minimização consiste em olhar a derivada da função e decidir, a partir desta informação,\n",
    "para qual direção devemos procurar o mínimo. Como a derivada fornece a direção de crescimento, devemos andar sempre para a esquerda se a derivada for positiva (função crescendo no ponto) ou para a direita se a derivada for negativa.\n",
    "\n",
    "Vamos decidir sobre um passo e realizar algumas iterações desta estratégia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def f_linha(x):\n",
    "    return (x  - 1) * np.exp(-x)\n",
    "\n",
    "# Calculamos vários valores de x obtidos com a estratégia de andar na direção oposta à derivada\n",
    "x = 5\n",
    "passo = 0.3\n",
    "X = [x]\n",
    "for i in range(30):\n",
    "    deriv = f_linha(x)\n",
    "    if deriv > 0:\n",
    "        x -= passo\n",
    "    else:\n",
    "        x += passo\n",
    "    X.append(x)\n",
    "    \n",
    "# Fazemos um gráfico dos resultados\n",
    "plt.plot(X, 'ko--', label='aproximação via derivada')\n",
    "plt.plot([1] * len(X), 'r-', lw=2, label='valor correto')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Temos um problema! O algoritmo anda na direção correta até chegar próximo ao mínimo. A partir daí, começa a oscilar entre dois pontos em torno do mínimo sem nunca se aproximar verdadeiramente do mínimo. O motivo para que isto aconteça é que o passo é fixo. Talvez fosse melhor que o passo fosse proporcional ao próprio valor da derivada. Desta forma, quando a derivada fosse menor que zero, daríamos passos menores para aproximar mais lentamente do mínimo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x = 5\n",
    "alpha = 3.0\n",
    "X = [x]\n",
    "for i in range(30):\n",
    "    deriv = f_linha(x)\n",
    "    passo = -alpha * deriv\n",
    "    x += passo\n",
    "    X.append(x)\n",
    "    \n",
    "# Fazemos um gráfico dos resultados\n",
    "plt.plot(X, 'ko--', label='aproximação via derivada')\n",
    "plt.plot([1] * len(X), 'r-', lw=2, label='valor correto')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Bem melhor!\n",
    "\n",
    "Agora queremos diminuir o número de iterações calibrando melhor o tamanho de cada passo. Vemos que é possível dar passos bem maiores que utilizamos para construir a figura acima. Isto pode ser feito controlando a variável alpha. No entanto, o ideal é que o próprio alpha evolua para que seja possível utilizar valores altos quando a função é mais previsível e valores menores quando uma maior precisão for necessária. \n",
    "\n",
    "Nossa estratégia é aumentar o valor de alpha a cada iteração e reduzí-lo sempre que o passo tente mandar para um novo ponto que aumenta o valor da função ao invés de diminuí-lo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x = 5\n",
    "alpha = 3.0\n",
    "r_mais = 1.5\n",
    "r_menos = 0.5\n",
    "fmin = f(x)\n",
    "\n",
    "X = [x]\n",
    "for i in range(30):\n",
    "    deriv = f_linha(x)\n",
    "    passo = -alpha * deriv\n",
    "    \n",
    "    # Limita a 10 iterações\n",
    "    for i in range(10):\n",
    "        fnew = f(x + passo)\n",
    "        if fnew < fmin:\n",
    "            fmin = fnew\n",
    "            break\n",
    "        else:\n",
    "            alpha *= r_menos\n",
    "            passo = -alpha * deriv\n",
    "        \n",
    "    x += passo\n",
    "    alpha *= r_mais\n",
    "    \n",
    "    X.append(x)\n",
    "    Y.append(f(x))\n",
    "    \n",
    "# Fazemos um gráfico dos resultados\n",
    "plt.plot(X, 'ko--', label='aproximação via derivada')\n",
    "plt.plot([1] * len(X), 'r-', lw=2, label='valor correto')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n",
    "print('xmin:', x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Enquanto a primeira versão pedia 20 iterações esta se aproxima bastante do resultado com apenas 10 iterações.\n",
    "Vamos à implementação."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def gradiente(f, f_linha, x0, xtol=1e-6, dtol=1e-6, alpha=5, r_mais=1.5, r_menos=0.25):\n",
    "    \"\"\"\n",
    "    Retorna o mínimo da função f(x) com derivada f_linha fazendo uma descida de gradiente.\n",
    "    \n",
    "    Args:\n",
    "        f: função objetivo f(x).\n",
    "        f_linha: derivada de f(x).\n",
    "        x1, x2: início e fim do intervalo de busca.\n",
    "        xtol: tolerância no intervalo de busca.\n",
    "        dtol: tolerância na derivada.\n",
    "        alpha: constante que define tamanho inicial do passo.\n",
    "        r_mais: fator de expansão do passo a cada iteração (r_mais >= 1).\n",
    "        r_menos: fator de redução do passo a cada iteração mal sucedida (r_menos <= 1).\n",
    "        \n",
    "    Returns:\n",
    "        Retorna o valor de x que minimiza f(x).\n",
    "    \"\"\"\n",
    "    x = x0\n",
    "    fmin = f(x)\n",
    "    while True:\n",
    "        deriv = f_linha(x)\n",
    "        if abs(deriv) < dtol:\n",
    "            return x\n",
    "        passo = -alpha * deriv\n",
    "\n",
    "        # Limita a 10 iterações de redução do alpha\n",
    "        for i in range(10):\n",
    "            fnew = f(x + passo)\n",
    "            if fnew < fmin:\n",
    "                fmin = fnew\n",
    "                break\n",
    "            else:\n",
    "                alpha *= r_menos\n",
    "                passo = -alpha * deriv\n",
    "\n",
    "        x += passo\n",
    "        alpha *= r_mais\n",
    "        if abs(passo) < xtol:\n",
    "            return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Testamos para funções conhecidas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "gradiente(f, f_linha, 5, dtol=0), gradiente(np.cos, lambda x: -np.sin(x), 5, dtol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "E comparamos a performance..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 razao_aurea(f, 0, 10)\n",
    "%timeit -n1000 gradiente(f, f_linha, 5, dtol=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 razao_aurea(f, 0, 10, xtol=1e-10)\n",
    "%timeit -n1000 gradiente(f, f_linha, 5, xtol=1e-10, dtol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que o método da descida do gradiente não supera a razão áurea em todos exemplos.\n",
    "\n",
    "O principal problema talvez seja que os métodos de descida do gradiente tendem a se aproximar rapidamente de uma vizinhança do mínimo mas normalmente são relativamente lentos em finalmente atingir este valor. Isto está relacionado ao fato que cada passo é proporcional ao valor do gradiente, sendo que este tende a zero próximo do mínimo. \n",
    "\n",
    "Uma solução para este dilema é mudar de método: vamos estudar o método das aproximações parabólicas e posteriormente implementar uma estratégia que tente extrair o melhor de cada um."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Método de Newton\n",
    "\n",
    "Uma estratégia muito bem sucedida consiste em aproximar a função por parábolas e estimar o mínimo a cada iteração pelo mínimo desta parábola. Muitas vezes fazemos esta estimativa utilizando a derivada segunda da função (também chamada de hessiana).\n",
    "\n",
    "Utilizamos para isto a expansão da função em série de Taylor:\n",
    "\n",
    "$$f(x) \\simeq f(x_0) + f'(x_0)(x - x_0) + \\frac12 f''(x_0)(x - x_0)^2$$\n",
    "\n",
    "Tirando a derivada, calculamos o mínimo desta parábola\n",
    "\n",
    "$$f'(x_0) + f''(x_0)(x - x_0) = 0$$\n",
    "\n",
    "Ou seja:\n",
    "\n",
    "$$x = x_0 - \\frac{f'(x_0)}{f''(x_0)}$$\n",
    "\n",
    "Os gráficos abaixo ilustram a otimização com base na Hessiana"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def hess(x):\n",
    "    return (2 - x) * np.exp(-x)\n",
    "\n",
    "x0 = 1.4\n",
    "X = np.linspace(0, 10, 100)\n",
    "for i in range(5):\n",
    "    f0 = f(x0)\n",
    "    f1 = f_linha(x0)\n",
    "    f2 = hess(x0)\n",
    "    \n",
    "    plt.plot(X, f(X), 'k-', label='função')\n",
    "    plt.plot(X, f0 + f1 * (X - x0) + 0.5 * f2 * (X - x0)**2, label='aproximação quadrática')\n",
    "    plt.plot([x0], [f0], 'ko')\n",
    "    plt.axis([0, 10, -0.5, 0.25])\n",
    "    plt.axis([0, 10, -0.5, 0.25])\n",
    "    plt.show()\n",
    "    \n",
    "    x0 -= f1 / f2\n",
    "    print('xmin:', x0)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Observe que uma vez que o método se aproxima do mínimo, a convergência é muito rápida. Apesar disto, trata-se de um método frágil. Se não estivermos próximos do mínimo e a função não for muito bem descrita por uma parábola, o método gera valores absurdos, como pode-se experimentar escolhendo um valor de $x$ inicial longe de $x=1$.\n",
    "\n",
    "Outro problema é a necessidade de fornecer explicitamente a Hessiana. Isto deixa o método mais trabalhoso de implementar, ainda que a performance próxima do mínimo seja muito boa.\n",
    "\n",
    "A boa notícia é que o método de Newton e o método de descida do gradiente são complementares: o primeiro consegue se aproximar rapidamente do mínimo enquanto que o outro é muito eficiente mas falha se não for inicializado o suficientemente próximo do mínimo. Vamos tentar uma abordagem mista: fazemos a descida do gradiente normalmente, mas em alguns passos tentamos utilizar o método de Newton aproximando a Hessiana por duas avaliações seguidas da derivada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def gradiente_newton(f, f_linha, x0, xtol=1e-6, dtol=1e-6, alpha=5, r_mais=1.5, r_menos=0.25):\n",
    "    \"\"\"\n",
    "    Retorna o mínimo da função f(x) com derivada f_linha fazendo uma descida de gradiente.\n",
    "    \n",
    "    Args:\n",
    "        f: função objetivo f(x).\n",
    "        f_linha: derivada de f(x).\n",
    "        x1, x2: início e fim do intervalo de busca.\n",
    "        xtol: tolerância no intervalo de busca.\n",
    "        dtol: tolerância na derivada.\n",
    "        alpha: constante que define tamanho inicial do passo.\n",
    "        r_mais: fator de expansão do passo a cada iteração (r_mais >= 1).\n",
    "        r_menos: fator de redução do passo a cada iteração mal sucedida (r_menos <= 1).\n",
    "        \n",
    "    Returns:\n",
    "        Retorna o valor de x que minimiza f(x).\n",
    "    \"\"\"\n",
    "    x = x0\n",
    "    fmin = f(x)\n",
    "    while True:\n",
    "        deriv = f_linha(x)\n",
    "        if abs(deriv) < dtol:\n",
    "            return x\n",
    "        passo = -alpha * deriv\n",
    "\n",
    "        # Limita a 10 iterações de redução do alpha\n",
    "        for i in range(10):\n",
    "            fnew = f(x + passo)\n",
    "            if fnew < fmin:\n",
    "                fmin = fnew\n",
    "                \n",
    "                # Aqui chegamos a um valor possivelmente próximo do zero.\n",
    "                # Vamos tentar aproximar a Hessiana pela formula abaixo\n",
    "                hess = (f_linha(x + passo) - deriv) / passo\n",
    "                \n",
    "                # Se a hessiana for negativa, temos um máximo, caso contrário, \n",
    "                # temos um mínimo\n",
    "                if hess > 0:\n",
    "                    xnew = x - deriv/hess\n",
    "                    f_newton = f(xnew)\n",
    "                    if f_newton < fmin:\n",
    "                        fmin = f_newton \n",
    "                        x = xnew - passo  # tiramos um passo pois ele será \n",
    "                                          # acrescentado posteriormente\n",
    "                break\n",
    "            else:\n",
    "                alpha *= r_menos\n",
    "                passo = -alpha * deriv\n",
    "\n",
    "        x += passo\n",
    "        alpha *= r_mais\n",
    "        if abs(passo) < xtol:\n",
    "            return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Testamos o método"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "gradiente_newton(f, f_linha, 5, dtol=0), gradiente_newton(np.cos, lambda x: -np.sin(x), 5, dtol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora vamos analizar a performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 razao_aurea(f, 0, 10)\n",
    "%timeit -n1000 gradiente(f, f_linha, 5, dtol=0)\n",
    "%timeit -n1000 gradiente_newton(f, f_linha, 5, dtol=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 razao_aurea(f, 0, 10, xtol=1e-10)\n",
    "%timeit -n1000 gradiente(f, f_linha, 5, xtol=1e-10, dtol=0)\n",
    "%timeit -n1000 gradiente_newton(f, f_linha, 5, xtol=1e-10, dtol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos então que o método torna-se mais eficaz quando exige-se uma maior precisão."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}