SciPy (произносится как сай пай) — это основанный на numpy математический пакет, включающий в себя также библиотеки на C и Fortran. С SciPy интерактивный сеанс Python превращается в такую же полноценную среду обработки данных, как MATLAB, IDL, Octave, R или SciLab.
В этой статье рассмотрим основные приемы математического программирования — решения задач условной оптимизации для скалярной функции нескольких переменных с помощью пакета scipy.optimize. Алгоритмы безусловной оптимизации уже рассмотрены в прошлой статье. Более подробную и актуальную справку по функциям scipy всегда можно получить с помощью команды help(), Shift+Tab или в официальной документации.
Введение
Общий интерфейс для решения задач как условной, так и безусловной оптимизации в пакете scipy.optimize предоставляется функцией minimize()
. Однако известно, что универсального способа для решения всех задач не существует, поэтому выбор адекватного метода как всегда ложится на плечи исследователя.
Подходящий алгоритм оптимизации задается с помощью аргумента функции minimize(..., method="")
.
Для условной оптимизации функции нескольких переменных доступны реализации следующих методов:
trust-constr
— поиск локального минимума в доверительной области. Статья на wiki, статья на хабре;SLSQP
— последовательное квадратичное программирование с ограничениями, ньютоновский метод решения системы Лагранжа. Статья на вики.TNC
— Truncated Newton Constrained, ограниченное число итераций, хорош для нелинейных функций с большим числом независимых переменных. Статья на wiki.L-BFGS-B
— метод от четверки Broyden–Fletcher–Goldfarb–Shanno, реализованный с уменьшенным потреблением памяти за счет частичной загрузки векторов из матрицы Гессе. Статья на wiki, статья на хабре.COBYLA
—КОБЫЛАConstrained Optimization By Linear Approximation, ограниченная оптимизация с линейной аппроксимацией (без вычисления градиента). Статья на wiki.
В зависимости от выбранного метода, по-разному задаются условия и ограничения для решения задачи:
- объектом класса
Bounds
для методов L-BFGS-B, TNC, SLSQP, trust-constr; - списком
(min, max)
для этих же методов L-BFGS-B, TNC, SLSQP, trust-constr; - объектом или списком объектов
LinearConstraint
,NonlinearConstraint
для методов COBYLA, SLSQP, trust-constr; - словарем или списком словарей
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}
для методов COBYLA, SLSQP.
План статьи:
1) Рассмотреть применение алгоритма условной оптимизации в доверительной области (method="trust-constr") с ограничениями, заданными в виде объектов Bounds
, LinearConstraint
, NonlinearConstraint
;
2) Рассмотреть последовательное программирование методом наименьших квадратов (method="SLSQP") с ограничениями, заданными в виде словаря {'type', 'fun', 'jac', 'args'}
;
3) Разобрать пример оптимизации выпускаемой продукции на примере веб-студии.
Условная оптимизация method="trust-constr"
Реализация метода trust-constr
основана на EQSQP для задач с ограничениями вида равенства и на TRIP для задач с ограничениями в виде неравенств. Оба метода реализованы алгоритмами поиска локального минимума в доверительной области и хорошо подходят для крупномасштабных задач.
Математическая постановка задачи поиска минимума в общем виде:
Для ограничений строгого равенства нижняя граница устанавливается равной верхней .
Для одностороннего ограничения верхняя или нижняя граница устанавливается np.inf
с соответствующим знаком.
Пусть необходимо найти минимум известной функцию Розенброка от двух переменных:
При этом заданы следующие ограничения на ее область определения:
В нашем случае имеется единственное решение в точке , для которой справедливы только первое и четвертое ограничения.
Пройдемся по ограничениям снизу вверх и рассмотрим, как можно их записать в scipy.
Ограничения и определим с помощью объекта Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])
Ограничения и запишем в линейной форме:
Определим эти ограничения в виде объекта LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])
И наконец нелинейное ограничение в матричной форме:
Определим матрицу Якоби для этого ограничения и линейную комбинацию матрицы Гессе с произвольным вектором :
Теперь нелинейное ограничение можем определить как объект NonlinearConstraint
:
from scipy.optimize import NonlinearConstraint
def cons_f(x):
return [x[0]**2 + x[1], x[0]**2 - x[1]]
def cons_J(x):
return [[2*x[0], 1], [2*x[0], -1]]
def cons_H(x, v):
return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
Если размер велик, матрицы можно задавать и в разреженном виде:
from scipy.sparse import csc_matrix
def cons_H_sparse(x, v):
return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_sparse)
или как объект LinearOperator
:
from scipy.sparse.linalg import LinearOperator
def cons_H_linear_operator(x, v):
def matvec(p):
return np.array([p[0]*2*(v[0]+v[1]), 0])
return LinearOperator((2, 2), matvec=matvec)
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_linear_operator)
Когда вычисление матрицы Гессе требует больших затрат, можно использовать класс HessianUpdateStrategy
. Доступны следующие стратегии: BFGS
и SR1
.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())
Гессиан также может быть вычислен с помощью конечных разностей:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')
Матрицу Якоби для ограничений также можно вычислить с помощью конечных разностей. Однако, в этом случае матрицу Гессе с помощью конечных разностей уже не вычислить. Гессиан должен быть определен в виде функции или с помощью класса HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())
Решение оптимизационной задачи выглядит следующим образом:
from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod
x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]
При необходимости, функцию вычисления гессиана можно определить с помощью класса LinearOperator
def rosen_hess_linop(x):
def matvec(p):
return rosen_hess_prod(x, p)
return LinearOperator((2, 2), matvec=matvec)
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)
или произведение Гессиана и произвольного вектора через параметр hessp
:
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)
Альтернативно, первая и вторая производные оптимизируемой функции могут быть вычислены приближенно. Например, гессиан может быть аппроксимирован с помощью функции SR1
(квази-Ньютоновского приближения). Градиент может быть аппроксимирован конечными разностями.
from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(),
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)
Условная оптимизация method="SLSQP"
Метод SLSQP предназначен для решения задач минимизации функции в виде:
Где и — множества индексов выражений, описывающих ограничения в виде равенств или неравенств. — множества нижних и верхних границ для области определения функции.
Линейные и нелинейные ограничения описываются в виде словарей с ключами type
, fun
и jac
.
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
1 - x [0] ** 2 - x [1],
1 - x [0] ** 2 + x [1]]),
'jac': lambda x: np.array ([[- 1.0, -2.0],
[-2 * x [0], -1.0],
[-2 * x [0], 1.0]])
}
eq_cons = {'type': 'eq',
'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
'jac': lambda x: np.array ([2.0, 1.0])
}
Поиск минимума осуществляется следующим образом:
x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
bounds=bounds)
print(res.x)
Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]
Пример оптимизации
В связи с переходом к пятому технологическому укладу, рассмотрим оптимизацию производства на примере веб-студии, которая приносит нам небольшой, но стабильный доход. Представим себя директором галеры, на которой производится три вида продукции:
- x0 — продающие лэндинги, от 10 т.р.
- x1 — корпоративные сайты, от 20 т.р.
- x2 — интернет магазины, от 30 т.р.
Наш дружный рабочий коллектив включает в себя четырех джунов, двух мидлов и одного сеньора. Фонд их рабочего времени на месяц:
- джуны:
4 * 150 = 600 чел * час
, - мидлы:
2 * 150 = 300 чел * час
, - сеньор:
150 чел * час
.
Пусть на разработку и деплой одного сайта типа (x0, x1, x2) первый попавшийся джуниор должен потратить (10, 20, 30) часов, мидл — (7, 15, 20), сеньор — (5, 10, 15) часов лучшего времени своей жизни.
Как любому нормальному директору, нам хочется максимизировать ежемесячную прибыль. Первый шаг к успеху — записываем целевую функцию value
как сумму доходов от произведенной за месяц продукции:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]
Это не ошибка, при поиске максимума целевая функция минимизируется с обратным знаком.
Следующий шаг — запрещаем перерабатывать своим сотрудникам и вводим ограничения на фонд рабочего времени:
Что эквивалентно:
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
300 - 7 * x [0] - 15 * x [1] - 20 * x[2],
150 - 5 * x [0] - 10 * x [1] - 15 * x[2]])
}
Формальное ограничение — выпуск продукции должен быть только положительным:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])
И наконец самое радужное допущение — из-за низкой цены и высокого качества к нам постоянно выстраивается очередь из довольных клиентов. Мы можем сами выбирать ежемесячные объемы производства продукции, исходя из решения задачи условной оптимизации с scipy.optimize
:
x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)
[7.85714286 5.71428571 3.57142857]
Нестрого округлим до целых и посчитаем месячную загрузку гребцов при оптимальном раскладе продукции x = (8, 6, 3)
:
- джуны:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час
; - мидлы:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час
; - сеньор:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час
.
Вывод: чтобы директор получал свой заслуженный максимум, оптимально делать в месяц по 8 лэндингов, 6 средних сайтов и 3 магазина. Сеньор при этом должен пахать не отрываясь от станка, загрузка мидлов составит примерно 2/3, джунов меньше половины.
Заключение
В статье изложены основные приемы работы с пакетом scipy.optimize
, используемые для решения задач условной минимизации. Лично я использую scipy
чисто в академических целях, поэтому приведенный пример носит такой шуточный характер.
Много теории и винрарных примеров можно найти, например, в книге И.Л.Акулича "Математическое программирование в примерах и задачах". Более хардкорное применение scipy.optimize
для построения 3D структуры по набору изображений (статья на хабре) можно посмотреть в scipy-cookbook.
Основной источник информации — docs.scipy.org, желающие поконтрибьютить в перевод этого и других разделов scipy
добро пожаловать на GitHub.
Спасибо mephistopheies за участие в подготовке публикации.
Комментарии (5)
balezz Автор
17.04.2019 22:03Задачи у нас в основном учебные, поэтому про оптимизацию зашумленных функций вряд ли что смогу рассказать. А откуда берется шум? Вы его моделируете или это в результате измерений получается?
akhkmed
18.04.2019 09:55+1В моём случае шум происходил от округления, и градиенты было не посчитать. Если делать оптимизацию параметров для реальных систем, а они несовершенны, то источников шума будет много.
Вы верно заметили, что при моделировании шум и ошибки возможны. Например, в методе Монте-Карло.
akhkmed
Спасибо за руководство.
Как я понимаю, для зашумлённых функций невозможно вычислить градиент, а значит, будут работать не все методы. Как пишут на scipy-lectures.org, рекомендуется Nelder-Mead или Powell. Экспериментально установлено, что COBYLA тоже не вычисляет градиент. То есть, получается, это единственный метод без градиента, в котором можно задать ограничения. С другой стороны, COBYLA делает попытки вычисления даже за ограничениями, это тоже нужно учитывать.
Поделитесь, пожалуйста, опытом, как для Ваших задач работает COBYLA и использовали ли вы его для зашумлённых функций?
iroln
Посмотрите этот пакет
https://github.com/facebookresearch/nevergrad
А также, возможно, будет интересно вот это:
https://github.com/uqfoundation/mystic
malkovsky
А я бы еще порекомендовал
www.cvxpy.org
UPD. Пардон, я не посмотрел суть обсуждения, пакет выше — это вообще говоря DSL для описания задач выпуклой оптимизации, но в общем предполагается использование градиентных методов