SciPy (вимовляється як сай пай) - це пакет прикладних математичних процедур, заснований на розширенні Numpy Python. З SciPy інтерактивний сеанс Python перетворюється на таке ж повноцінне середовище обробки даних та прототипування складних систем, як MATLAB, IDL, Octave, R-Lab та SciLab. Сьогодні хочу коротко розповісти про те, як слід застосовувати деякі відомі алгоритми оптимізації в пакеті scipy.optimize. Докладнішу та актуальну довідку щодо застосування функцій завжди можна отримати за допомогою команди help() або за допомогою Shift+Tab.
Запровадження
Щоб позбавити себе і читачів від пошуку і читання першоджерел, посилання описи методів будуть переважно на вікіпедію. Як правило, цієї інформації достатньо для розуміння методів у загальних рисах та умов їх застосування. Для розуміння суті математичних методів йдемо за посиланнями більш авторитетні публікації, які можна знайти в кінці кожної статті або в улюбленій пошуковій системі.
Отже, модуль scipy.optimize включає реалізацію наступних процедур:
- Умовної та безумовної мінімізації скалярних функцій кількох змінних (minim) за допомогою різних алгоритмів (симплекс Нелдера-Міда, BFGS, сполучених градієнтів Ньютона, и )
- Глобальної оптимізації (наприклад: , )
- Мінімізація залишків (least_squares) та алгоритми припасування кривих нелінійним МНК (curve_fit)
- Мінімізації скалярної функції однієї змінної (minim_scalar) і пошуку коренів (root_scalar)
- Багатовимірні вирішувачі системи рівнянь (root) з використанням різних алгоритмів (гібридний Пауелла, або великомасштабні методи, такі як ).
У цій статті ми розглянемо лише перший пункт із усього цього списку.
Безумовна мінімізація скалярної функції кількох змінних
Функція minim з пакету scipy.optimize надає загальний інтерфейс для вирішення завдань умовної та безумовної мінімізації скалярних функцій кількох змінних. Щоб продемонструвати її роботу, нам знадобиться відповідна функція кількох змінних, яку ми по-різному мінімізуватимемо.
Для цих цілей чудово підійде функція Розенброка від N змінних, яка має вигляд:
Незважаючи на те, що функція Розенброка та її матриці Якобі та Гессе (першої та другої похідної відповідно) вже визначені в пакеті scipy.optimize, визначимо її самостійно.
import numpy as np
def rosen(x):
"""The Rosenbrock function"""
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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 339
Function evaluations: 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 1622
[1. 1. 1. 1. 1.]Алгоритм Бройдена-Флетчера-Голдфарба-Шанно (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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 25
Function evaluations: 30
Gradient evaluations: 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 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})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 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})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 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})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 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.xOptimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 14
Gradient evaluations: 13
Hessian evaluations: 14
array([1., 1., 1., 1., 1.])На цьому, мабуть, зупинимося. У наступній статті постараюся розповісти найцікавіше про умовну мінімізацію, додаток мінімізації у вирішенні завдань апроксимації, мінімізації функції однієї змінної, довільних мінімізаторів та пошук коренів системи рівнянь за допомогою пакету scipy.optimize.
Джерело:
Джерело: habr.com