{
"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$ - признаковое описание наблюдений,
$y$ - прогнозируемая величина\n",
"\n",
"Пусть задана функция ошибки (функция потерь) $L(\\cdot)$.
\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)$
\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