SciPy, оптимизация

SciPy (произнася като сай пай) — това е пакет от приложни математически процедури, основан на разширението Numpy на Python. С SciPy интерактивната сесия на Python се превръща в напълно функционална среда за обработка на данни и прототипиране на сложни системи, подобно на MATLAB, IDL, Octave, R-Lab и SciLab. Днес искам накратко да разкажа как да се прилагат някои известни алгоритми за оптимизация в пакета scipy.optimize. По-подробна и актуална справка за използването на функциите винаги може да бъде получена с командата help() или с Shift+Tab.

Въведение

За да спестя на себе си и на читателите търсенето и четенето на първоизточниците, линковете към описанията на методите ще бъдат предимно към Wikipedia. Обикновено тази информация е достатъчна за разбиране на методите в общи линии и условията за тяхното прилагане. За да разберем същността на математическите методи, ще се насочим към линковете към по-авторитетни публикации, които могат да бъдат намерени в края на всяка статия или в любимата ни търсачка.

И така, модулът scipy.optimize включва в себе си реализация на следните процедури:

1. Условна и безусловна минимизация на скалярни функции с множество променливи (minim) с помощта на различни алгоритми (симплекс на Нелдер-Мид, BFGS, сопряжени градиенти на Нютон, COBYLA и SLSQP)
2. Глобална оптимизация (например: basinhopping, diff_evolution)
3. Минимизация на остатъци МНК (least_squares) и алгоритми за подгонване на криви с нелинеен МНК (curve_fit)
4. Минимизация на скалярни функции с една променлива (minim_scalar) и търсене на корени (root_scalar)
5. Многомерни решатели на системи от уравнения (root) с използване на различни алгоритми (гибриден на Пауъл, Левенберг-Марквардт или крупномасштабни методи, такива като Нютон-Крилов).

В тази статия ще разгледаме само първия пункт от целия този списък.

Безусловна минимизация на скалярна функция с множество променливи

Функцията minim от пакета scipy.optimize предоставя общ интерфейс за решаване на задачи за условна и безусловна минимизация на скалярни функции с множество променливи. За да демонстрираме нейната работа, ще ни е необходима подходяща функция с множество променливи, която ще минимизираме по различни начини.

За тези цели идеално подхожда функцията на Розенброк от N променливи, която е с вида:

Въпреки че функцията на Розенброк и нейните Якоби и Хес матрици (първа и втора производна съответно) вече са определени в пакета scipy.optimize, ще я дефинираме сами.

import numpy as np

def rosen(x):
    """Функцията на Розенброк"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

За по-голяма яснота ще изобразим в 3D стойностите на функцията на Розенброк от две променливи.

Код за визуализация

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Настройка на 3D графика
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Задаване на ъгъла на виждане
ax.view_init(45, 30)

# Създаване на данни за графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))

# Изчертаване на повърхността
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

Знаейки предварително, че минимумът е 0 при , нека разгледаме примери за това как да определим минималната стойност на функцията на Розенброк с помощта на различни процедури от scipy.optimize.

Методът на Нелдера-Мида (Nelder-Mead)

Нека имаме начална точка x0 в 5-мерното пространство. Ще намерим най-близката до нея точка на минимума на функцията на Розенброк с помощта на алгоритъма на симплекса Нелдер-Мид (алгоритъмът е зададен като стойност на параметъра method):

from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизацията успешно приключи.
         Текуща стойност на функцията: 0.000000
         Итерации: 339
         Оценки на функцията: 571
[1. 1. 1. 1. 1.]

Методът на симплекса е най-простият начин за минимизиране на явно дефинирана и доста гладка функция. Не изисква изчисляване на производни на функцията, достатъчно е да зададете само нейните стойности. Методът на Нелдер-Мид е добро решение за прости задачи по минимизиране. Въпреки това, тъй като не използва оценки на градиента, за намирането на минимума може да е нужно повече време.

Методът на Пауъл

Друг алгоритъм за оптимизация, който изчислява само стойностите на функциите, е методът на Пауъл. За да го използвате, трябва да зададете method = ‘powell’ във функцията minim.

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизацията успешно приключи.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функцията: 1622
[1. 1. 1. 1. 1.]

Алгоритъмът на Бройдена-Флетчера-Голдфарба-Шанно (BFGS)

За постигане на по-бърза конвергенция към решението, процедурата BFGS използва градиент на целевата функция. Градиентът може да се зададе под формата на функция или да се изчисли с помощта на разлики от първи ред. Във всеки случай методът BFGS обикновено изисква по-малко извиквания на функции в сравнение с симплекс метода.

Нека намерим производната на функцията на Розенброк в аналитичен вид:

Това изразение важи за производните на всички променливи, освен за първата и последната, които се определят като:

Нека разгледаме функция на Python, която изчислява този градиент:

def rosen_der (x):
    xm = x [1: -1]
    xm_m1 = x [: - 2]
    xm_p1 = x [2:]
    der = np.zeros_like (x)
    der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
    der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
    der [-1] = 200 * (x [-1] -x [-2] ** 2)
    return der

Функцията за изчисляване на градиента се указва като стойност на параметъра jac на функцията minim, както е показано по-долу.

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Оптимизацията приключи успешно.
         Текуща стойност на функцията: 0.000000
         Итерации: 25
         Оценки на функцията: 30
         Оценки на градиента: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

Алгориъм на спрягатни градиенти (Нютон)

Алгоритъм съпрегатни градиенти на Нютон е модифициран метод на Нютон.
Методът на Нютон е основан на апоксимацията на функцията в локална област посредством полином от втори порядък:

където е матрица на вторите производни (матрица на Хессе, хесиан).
Ако хесианът е положително определен, то локалният минимум на тази функция може да се намери, като се приравни нулевият градиент на квадратичната форма към нула. В резултат ще получим следното изражение:

Обратният хесиан се изчислява с помощта на метода на спрягатни градиенти. Пример за използване на този метод за минимизиране на функцията на Розенброк е показан по-долу. За да използвате метода на Newton-CG, е необходимо да зададете функция, която изчислява хесиана.
Хесианът на функцията на Розенброк в аналитичен вид е равен на:

където и , определят матрицата .

Останалите ненулеви елементи на матрицата са равни на:

Например, в петизмерно пространство N = 5, матрицата на Хессе за функцията на Розенброк има лентов вид:

Кодът, който изчислява този хесиан, заедно с кода за минимизиране на функцията на Розенброк с помощта на метода на спрягатни градиенти (Нютон):

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

res = minimize(rosen, x0, method='Newton-CG', 
               jac=rosen_der, hess=rosen_hess,
               options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация завершилась успешно.
         Текущая функция: 0.000000
         Итерации: 24
         Оценка функции: 33
         Оценки градиента: 56
         Оценки Гессе: 24
[1.         1.         1.         0.99999999 0.99999999]

Пример определения функции произведения гессиана и произвольного вектора

В реальных задачах вычисление и хранение всей матрицы Гессе может требовать значительных ресурсов времени и памяти. При этом фактически нет необходимости задавать саму матрицу Гессе, так как для процедуры минимизации нужен только вектор, равный произведению гессиана с другим произвольным вектором. Таким образом, с вычислительной точки зрения гораздо предпочтительнее сразу определить функцию, которая возвращает результат произведения гессиана с произвольным вектором.

Рассмотрим функцию hess, которая принимает вектор минимизации в качестве первого аргумента, а произвольный вектор — в качестве второго аргумента (вместе с другими аргументами минимизируемой функции). В нашем случае вычислить произведение гессиана функции Розенброка с произвольным вектором не очень сложно. Если p — произвольный вектор, то произведение имеет вид:

Функция, вычисляющая произведение гессиана и произвольного вектора, передается как значение аргумента hessp функции minimize:

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] 
    -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

Оптимизация завершилась успешно.
         Текущая функция: 0.000000
         Итерации: 24
         Оценка функции: 33
         Оценки градиента: 56
         Оценки Гессе: 66

Алгоритм доверительной области (trust region) сопряженных градиентов (Ньютона)

Плохая обусловленность матрицы Гессе и неверные направления поиска могут привести к тому, что алгоритм сопряженных градиентов Ньютона может быть неэффективным. В таких случаях предпочтение отдается методу доверительной области (trust-region) сопряженных градиентов Ньютона.

Пример с определением матрицы Гессе:

res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация успешно завершена.
         Текущая функция: 0.000000
         Итерации: 20
         Оценки функций: 21
         Оценки градиента: 20
         Оценки гессиана: 19
[1. 1. 1. 1. 1.]

Пример с функцией произведения гессиана и произвольного вектора:

res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация успешно завершена.
         Текущая функция: 0.000000
         Итерации: 20
         Оценки функций: 21
         Оценки градиента: 20
         Оценки гессиана: 0
[1. 1. 1. 1. 1.]

Методи на основата на Крилов

Как и метод trust-ncg, методы на основе Крилова са добре подходящи за решаване на мащабни задачи, тъй като използват само матрични-векторни произведения. Същината им е в решаването на задачата в доверителна област, ограничена от отсеченото подпространство на Крилов. За неопределени задачи е по-добре да се използва този метод, тъй като изисква по-малко нелинейни итерации заради по-малкия брой матрични-векторни произведения на една подзадача, в сравнение с метода trust-ncg. Освен това решението на квадратичната подзадача е по-точно, отколкото при метода trust-ncg.
Пример с определением матрицы Гессе:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Оптимизацията е успешно приключила.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функции: 20
         Оценки на градиента: 20
         Оценки на гессиана: 18

print(res.x)

    [1. 1. 1. 1. 1.]

Пример с функцией произведения гессиана и произвольного вектора:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Оптимизацията е успешно приключила.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функции: 20
         Оценки на градиента: 20
         Оценки на гессиана: 0

print(res.x)

    [1. 1. 1. 1. 1.]

Алгоритъм за приближено решение в доверителна област

Всички методи (Newton-CG, trust-ncg и trust-krylov) са подходящи за решаване на мащабни задачи (с хиляди променливи). Това се дължи на факта, че алгоритъмът на методите на основата на свързаните градиенти предвижда приближено намиране на обратната матрица на гессиана. Решението се намира итеративно, без явното разлагане на гессиана. Тъй като е необходимо да се определи само функцията за произведението на гессиана и произволния вектор, този алгоритъм е особено полезен при работа с разредени (лентови диагонални) матрици. Това осигурява ниски разходи за памет и значителна икономия на време.

При задачах със среден размер разходите за съхранение и факторизация на Хесиевата матрица не играят решаваща роля. Това означава, че можем да постигнем решение с по-малко итерации, разрешавайки подсъществата в областта на доверието почти точно. За целта някои нелинейни уравнения се решават итеративно за всяка квадратична подзадача. Такова решение обикновено изисква 3 или 4 разлагания на Холецкия матрицата на Хес. В резултат методът конвергира с по-малко итерации и изисква по-малко изчисления на целевата функция в сравнение с други реализирани методи за област на доверие. Този алгоритъм предполага само определянето на пълната Хесиева матрица и не поддържа възможността за използване на функция за произведение на Хес и произволен вектор.

Пример с минимизация на функцията на Розенброк:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Оптимизацията завърши успешно.
         Текуща стойност на функцията: 0.000000
         Итерации: 13
         Оценки на функцията: 14
         Оценки на градиента: 13
         Оценки на Хес: 14

array([1., 1., 1., 1., 1.])

На това ще спрем. В следващата статия ще се опитам да разкажа най-интересното за условната минимизация, приложението на минимизацията в решаването на задачи по апроксимация, минимизацията на функции с една променлива, произволни минимизатори и търсенето на корени на системи от уравнения с помощта на пакета scipy.optimize.

Източник: https://docs.scipy.org/doc/scipy/reference/

Източник: habr.com