SciPy, оптимизация с условиями
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;
- словарем или списком словарей для методов COBYLA, SLSQP.
План статьи:
1) Рассмотреть применение алгоритма условной оптимизации в доверительной области (method=»trust-constr») с ограничениями, заданными в виде объектов Bounds , LinearConstraint , NonlinearConstraint ;
2) Рассмотреть последовательное программирование методом наименьших квадратов (method=»SLSQP») с ограничениями, заданными в виде словаря ;
3) Разобрать пример оптимизации инвестиционного портфеля с помощью метода SLSQP .
Условная оптимизация 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=, 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=, 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=, 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=, bounds=bounds) print(res.x)Условная оптимизация method=»SLSQP»
Метод SLSQP предназначен для решения задач минимизации функции в виде:
Где и — множества индексов выражений, описывающих ограничения в виде равенств или неравенств. — множества нижних и верхних границ для области определения функции.
Линейные и нелинейные ограничения описываются в виде словарей с ключами type , fun и jac .
Поиск минимума осуществляется следующим образом:
x0 = np.array([0.5, 0]) res = minimize(rosen, x0, method='SLSQP', jac=rosen_der, constraints=[eq_cons, ineq_cons], options=, 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 ]Пример оптимизации
UPD: Вдохновленный статьей Теория инвестиций автора abak, вместо абстрактного производства предлагаю рассмотреть портфельную оптимизацию по Марковицу. Решение этой задачи предполагает поиск оптимальных долей различных инструментов в инвестиционном портфеле с минимальным риском при заданном уровне доходности.
Вектором x обозначим доли различных инструментов в портфеле:
Вектор — ожидаемая доходность каждого инструмента, — ожидаемая дисперсия, — матрица взаимной корреляции между инструментами.
Минимизируется ожидаемая дисперсия, вычисляемая по формуле где S — матрица взаимной корреляции.
Ограничения в виде заданного уровня ожидаемой доходности r :
Доля каждого инструмента в портфеле должна быть положительной:
Сумма долей каждого инструмента в портфеле равна единице:
Обернем решение задачи условной оптимизации scipy.optimize в функцию от r и вычислим оптимальные доли инструментов для заданного уровня доходности r = [4..12]% .
from scipy.optimize import Bounds, minimize import numpy as np import matplotlib.pyplot as plt labels = ['Акции', 'Облигации', 'Недвижимость', 'Золото'] mu = np.array([ 10.9, 5.2, 10.8, 7.0]) var = np.array([15.2, 3.6, 19.2, 15.6]) R = np.array( [[1., 0., 0.59, 0.04], [0., 1.0, 0.19, 0.28], [0.59, 0.19, 1., 0.13], [0.04, 0.28, 0.13, 1.]]) var = np.expand_dims(var, axis=0) S = var.T @ var * R # Initial guess x = np.ones(4) * 0.25 def value(x): return x.T @ S @ x def optimize_portfolio(r): mu_cons = sum_cons = bnds = Bounds (np.zeros_like(x), np.ones_like(x) * np.inf) res = minimize(value, x, method='SLSQP', constraints=[mu_cons, sum_cons], bounds=bnds) return res.x rate = np.arange(4, 12, 0.1) y = np.array(list(map(optimize_portfolio, rate))).T plt.figure(figsize=(16, 6)) plt.stackplot(rate, y, labels=labels) plt.legend(loc='upper left') plt.show()Заключение
В статье изложены основные приемы работы с пакетом scipy.optimize , используемые для решения задач условной минимизации.
Много теории и винрарных примеров можно найти, например, в книге И.Л.Акулича «Математическое программирование в примерах и задачах». Более хардкорное применение scipy.optimize для построения 3D структуры по набору изображений (статья на хабре) можно посмотреть в scipy-cookbook.
Основной источник информации — docs.scipy.org, желающие поконтрибьютить в перевод документации scipy добро пожаловать на GitHub.
Спасибо mephistopheies за участие в подготовке публикации.