da_and_opt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Оптимизация функций. Символьные вычисления\n",
    "Шестаков А.В. Майнор по анализу данных - 02/02/2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сегодня мы рассмотрим следущие темы:\n",
    "\n",
    "    1) Методы оптимизации функций в Python и графика с ней связанная\n",
    "    2) Символьные вычисления\n",
    "    3) Решение задачи линейной регрессии тремя способами!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "import sympy\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Методы оптимизации функции в Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Выпуклые функции\n",
    "\n",
    "Рассмотрим две функции:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "\n",
    "y1 = lambda(x): 0.5*x**2 + 10*np.sin(x) - 2\n",
    "y2 = lambda(x): x**2 + 0.5**x - 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14,7))\n",
    "ax[0].plot(x, y1(x))\n",
    "ax[0].set_ylabel('$y_1$', fontsize=20)\n",
    "ax[0].set_xlabel('$x$', fontsize=20)\n",
    "\n",
    "ax[1].plot(x, y2(x))\n",
    "ax[1].set_ylabel('$y_2$', fontsize=20)\n",
    "ax[1].set_xlabel('$x$', fontsize=20)\n",
    "\n",
    "t = [0,5]\n",
    "alpha = 0.5\n",
    "f1 = alpha*y1(t[0]) + (1-alpha)*y1(t[1])\n",
    "f2 = y1(alpha*t[0] + (1-alpha)*t[1])\n",
    "\n",
    "ax[0].plot(alpha*t[0] + (1-alpha)*t[1], f1, '*r')\n",
    "ax[0].plot(alpha*t[0] + (1-alpha)*t[1], f2, 'sb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Функция $f(x)$ называется выпуклой (convex), если для любых двух точек $u$ и $v$ и для любого числа $\\alpha \\in [0,1]$ выполняется: \n",
    "$$ f\\left(\\alpha u + [1-\\alpha] v\\right) \\leq \\alpha f(u) + (1-\\alpha) f(v) $$\n",
    "\n",
    "Какая из этих функций является выпуклой? И заодно посмотрим, что это за выражение такое."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Чем хороши выпуклые фукнции с точки хрения поиска оптимального значения?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на выпуклую функцию в 3D!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "y = np.linspace(-5, 5, 100)\n",
    "\n",
    "X, Y = np.meshgrid(x,y)\n",
    "Z = X**2 + 0.25*Y**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(14, 7))\n",
    "ax = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax.view_init(40, 25)\n",
    "ax.plot_surface(X, Y, Z, alpha=0.3,)\n",
    "# ax.plot_(X, Y, Z)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "contour = ax.contour(X, Y, Z)\n",
    "plt.clabel(contour, inline=1, fontsize=10)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Гладкие функции"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гладкие функции (непрерывно дифференцируемые, smooth functions) удобны тем, что в любой точке можно посчитать её производую (градиент)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "\n",
    "y1 = lambda(x): 10*x**(2)\n",
    "y2 = lambda(x): abs(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14,7))\n",
    "ax[0].plot(x, y1(x))\n",
    "ax[0].set_ylabel('$y_1$', fontsize=20)\n",
    "ax[0].set_xlabel('$x$', fontsize=20)\n",
    "\n",
    "ax[1].plot(x, y2(x))\n",
    "ax[0].set_ylabel('$y_2$', fontsize=20)\n",
    "ax[0].set_xlabel('$x$', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Поиск экстемумов в Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В зависимости от того, какой функции вы смотрите в глаза, нужно выбирать свои методы. О описанием различных ситуаций и методов.\n",
    "\n",
    "Для вас есть хорошая новость. Как правило мы будем работать с выпуклыми и гладкими функциями. В общем случае это делается так:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Задаём функцию\n",
    "def f(x):\n",
    "    return x**2 + np.exp(x) - 4\n",
    "\n",
    "# Находим минимум\n",
    "x_min = opt.minimize(f, -3)\n",
    "x_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "plt.plot(x, f(x))\n",
    "plt.plot(x_min['x'], x_min['fun'], 'sb') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Символьные вычисления"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Символьные вычисления в Python выполняются в модуле `sympy`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Числа"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.init_printing(use_latex=True)\n",
    "\n",
    "a = sympy.Rational(1, 2)\n",
    "print a\n",
    "print a + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print sympy.pi\n",
    "print sympy.pi.evalf(n=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выражения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = sympy.Symbol('x')\n",
    "y, z = sympy.symbols('y z')\n",
    "w = (x + y)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.expand(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.factor(x**2 - y**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.simplify(w - 2*x*y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.simplify(sympy.sin(x)/sympy.cos(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Дифференцирование"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.diff(w, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.diff(w, x, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Решение уравнений"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sympy.solve(x**2 - 4, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = sympy.solve([x + 2*y - 2, \n",
    "                   x + 3*y - 6], x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При оптимизации функций, для некоторых методлов необходимо расчитать градиент функций.\n",
    "Используйте `sympy` если испытывайте какие-то сложности, затем, дополнительно проверяйте с `scipy.optimize.check_grad()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sympy.init_printing(use_latex=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Решение задачи линейной регрессии"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим 3 наблюдения:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = np.array([[1, 1],\n",
    "              [1, 2],\n",
    "              [1, 3]])\n",
    "y = np.array([1, 2, 2])\n",
    "\n",
    "plt.scatter(X[:,1], y)\n",
    "plt.xlabel('$x_1$')\n",
    "plt.ylabel('$y$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можем ли мы провести через них прямую?\n",
    "Почему?\n",
    "\n",
    "Как бы вы решали эту задачу, если бы нужно было провести прямую через 2 точки?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normal Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если система $$ A x = b $$ не имеет решения, то решайте $$A^\\top A x = A^\\top b$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Beta = np.linalg.inv((X.T.dot(X))).dot(X.T).dot(y)\n",
    "print Beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_hat = X.dot(Beta)\n",
    "\n",
    "plt.scatter(X[:,1], y)\n",
    "plt.plot(X[:,1], y_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Дифференцирование"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$X$ - признаковое описание наблюдений,<br\\> $y$ - прогнозируемая величина\n",
    "\n",
    "Пусть задана функция ошибки (функция потерь) $L(\\cdot)$. <br\\>\n",
    "Нам надо построить такой функционал $f(X)$, который будет выдавать значение наиболее близкие к $y$, иначе говоря: $$L\\left(f(X) - y\\right) \\rightarrow\\min $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Определим функцию потерь, как сумму квадратов разности выдаваемого ответа функционала и реального значения: \n",
    "$$ L(\\cdot) = \\frac{1}{2n}\\sum_{i=1}^n(f(x^{(i)}) - y^{(i)})^2$$\n",
    "\n",
    "Так как среди всего множества моделей мы выбрали линейную регрессию, то $$f(X) = \\beta_0 + \\beta_1x_1$$\n",
    "Подставляем это выражение в $L(\\cdot)$ и находим $\\beta_0$,\n",
    "$\\beta_1$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Градиентный спуск (gradient descent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Градиентый спукс - это итеративный метод оптимизации функции. Он заключается в постепенном перемещении к точке экспетмума в направлении градиента этой функции в точке.\n",
    "\n",
    "Пусть дана функция $L(\\theta)$, необходимо найти $\\hat{\\theta}:\\ L(\\hat{\\theta}) = \\min L(\\theta)$ <br\\>\n",
    "\n",
    "Шаги алгоритма:\n",
    "\n",
    "    1. Случайным образом фиксируется начальное состояние \n",
    "    2. Пока алгоритм не сойдется выполняется обновление:\n",
    "$$ \\theta = \\theta - \\alpha \\frac{\\partial L}{\\partial \\theta} $$\n",
    "    где $\\alpha$ - скорость спуска"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пример:**\n",
    "$$L(\\theta) = 3\\theta^2 + 2\\theta - 10$$\n",
    "Найдем градиент (в данном случае простую производую, так как у нас одна переменная)\n",
    "$$ \\frac{\\partial L}{\\partial \\theta}  = 6\\theta + 2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 3*x**2 + 2*x - 10\n",
    "\n",
    "def grad(x):\n",
    "    return 6*x + 2\n",
    "\n",
    "iter_max = 500 # максимальное кол-во итераций\n",
    "alpha = 0.01 # скорость спуска\n",
    "old_min = 0\n",
    "temp_min = 4\n",
    "precision = 0.001 # требуемая точность\n",
    "i = 0\n",
    "\n",
    "mins = [temp_min] # значений на каждой итерации алгоритма\n",
    "cost = [] # разность с инстинным решением\n",
    "\n",
    "# Напишите код, реализующий данный алгоритм\n",
    "# \n",
    "while abs(temp_min - old_min) > precision and i<iter_max:\n",
    "    ## ...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "\n",
    "# Рисуем функцию\n",
    "fig, ax = plt.subplots(1, 2, figsize=(14,7))\n",
    "ax[0].plot(x, f(x), 'b', linewidth=2)\n",
    "ax[0].set_ylabel('$ L(\\\\theta) $', fontsize=20)\n",
    "ax[0].set_xlabel('$\\\\theta$', fontsize=20)\n",
    "\n",
    "# Рисуем пусть градиента\n",
    "mins = np.array(mins)\n",
    "idx = np.linspace(0, len(mins)-1, 10).astype(int)\n",
    "\n",
    "ax[0].plot(mins[idx], f(mins[idx]), '-*')\n",
    "\n",
    "\n",
    "# Смотрим как уменьшается ошибка с итерацией\n",
    "ax[1].plot(cost)\n",
    "ax[1].set_ylabel('$Error$', fontsize=20)\n",
    "ax[1].set_xlabel('$Iteration$', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вернемся к линейной регрессии..\n",
    "\n",
    "Теперь, когда мы знаем, градиент функции квадратичной ошибки для линейной регрессии $\\beta_0$, $\\beta_1$ можно использовать метод градиентного спуска и для неё!\n",
    "\n",
    "**Важно!** Веса должны обновляться одновременно, а не по очереди!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Напишите код, реализующий данный алгоритм для линейной регрессии\n",
    "# \n",
    "def gradient_descent(X, y, iters, alpha):\n",
    "    costs = []\n",
    "    n = y.shape[0] \n",
    "    Beta = np.random.rand(2)\n",
    "    history = [Beta] \n",
    "    preds = []\n",
    "    for i in xrange(iters):\n",
    "        #....\n",
    "        \n",
    "    return history, costs, preds, Beta"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}