Как стать автором
Поиск
Написать публикацию
Обновить
77.1
Рейтинг
Open Data Science
Крупнейшее русскоязычное Data Science сообщество

SciPy, оптимизация с условиями

Время на прочтение7 мин
Количество просмотров48K
Блог компании Open Data SciencePython*Алгоритмы*Математика*Учебный процесс в IT
Туториал


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) Разобрать пример оптимизации инвестиционного портфеля с помощью метода SLSQP.


Условная оптимизация method="trust-constr"


Реализация метода trust-constr основана на EQSQP для задач с ограничениями вида равенства и на TRIP для задач с ограничениями в виде неравенств. Оба метода реализованы алгоритмами поиска локального минимума в доверительной области и хорошо подходят для крупномасштабных задач.


Математическая постановка задачи поиска минимума в общем виде:


$ \min_x f(x) $


$ c^l \leq c(x) \leq c^u ,$


$x^l \leq x \leq x^u$


Для ограничений строгого равенства нижняя граница устанавливается равной верхней $c^l_j = c^u_j$.
Для одностороннего ограничения верхняя или нижняя граница устанавливается np.inf с соответствующим знаком.
Пусть необходимо найти минимум известной функцию Розенброка от двух переменных:


$ \min_{x_0, x_1} 100(x_0 - x_1^2)^2 + (1-x_0)^2 $


При этом заданы следующие ограничения на ее область определения:


$ x_0^2 + x_1 \leq 1 $


$ x_0^2 - x_1 \leq 1 $


$ 2x_0 + x_1 = 1 $


$ x_0 + 2x_1 \leq 1 $


$ 0 \leq x_0 \leq 1 $


$ -0.5 \leq x_1 \leq 2.0 $


В нашем случае имеется единственное решение в точке $[x_0, x_1] = [0.4149, 0.1701]$, для которой справедливы только первое и четвертое ограничения.
Пройдемся по ограничениям снизу вверх и рассмотрим, как можно их записать в scipy.
Ограничения $0 \leq x_0 \leq 1$ и $-0.5 \leq x_1 \leq 2.0$ определим с помощью объекта Bounds.


from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])

Ограничения $x_0 + 2 x_1 \leq 1$ и $2 x_0 + x_1 = 1$ запишем в линейной форме:


$ \begin{bmatrix} - \infty \\ 1 \end{bmatrix} \leq \begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 1 \end{bmatrix}$


Определим эти ограничения в виде объекта LinearConstraint:


import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

И наконец нелинейное ограничение в матричной форме:


$ c(x) = \begin{bmatrix} x_0^2 + x_1 \\ x_0^2 - x_1 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 1 \end{bmatrix} $


Определим матрицу Якоби для этого ограничения и линейную комбинацию матрицы Гессе с произвольным вектором $v$:


$ J(x) = \begin{bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \end{bmatrix} $


$ H(x, v) = \sum\limits_{i=0}^1 v_i \nabla^2 c_i(x) = v_0 \begin{bmatrix} 2 & 0 \\ 2 & 0 \end{bmatrix} + v_1 \begin{bmatrix} 2 & 0 \\ 2 & 0 \end{bmatrix} $


Теперь нелинейное ограничение можем определить как объект 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)

Когда вычисление матрицы Гессе $H (x, v)$ требует больших затрат, можно использовать класс 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 предназначен для решения задач минимизации функции в виде:


$ \min_x f(x) $


$ c_j(x) = 0, j \in \mathcal {E} $


$ c_j(x) \geq 0, j \in \mathcal {I} $


$ lb_i \leq x_i \leq ub_i, i=1,...,N $


Где $\mathcal {E}$ и $\mathcal {I}$ — множества индексов выражений, описывающих ограничения в виде равенств или неравенств. $[lb_i, ub_i]$ — множества нижних и верхних границ для области определения функции.


Линейные и нелинейные ограничения описываются в виде словарей с ключами 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 ]

Пример оптимизации


UPD: Вдохновленный статьей Теория инвестиций автора abak, вместо абстрактного производства предлагаю рассмотреть портфельную оптимизацию по Марковицу. Решение этой задачи предполагает поиск оптимальных долей различных инструментов в инвестиционном портфеле с минимальным риском при заданном уровне доходности.
Вектором x обозначим доли различных инструментов в портфеле:


  • x[0] — Акции
  • x[1] — Облигации
  • x[2] — Недвижимость
  • x[3] — Золото

Вектор $\mu$ — ожидаемая доходность каждого инструмента, $\sigma$ — ожидаемая дисперсия, $R$ — матрица взаимной корреляции между инструментами.


Минимизируется ожидаемая дисперсия, вычисляемая по формуле $x^TSx$ где S — матрица взаимной корреляции.


$ S = \begin{bmatrix} \sigma_1^2 & \rho_{1,2}\sigma_1\sigma_2 & \cdots & \rho_{1,n}\sigma_1\sigma_n \\ \rho_{2,1}\sigma_2\sigma_1 & \sigma_2^2 & \cdots & \rho_{2,n}\sigma_2\sigma_n \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{n,1}\sigma_n\sigma_1 & \rho_{n,2}\sigma_n\sigma_2 & \cdots & \sigma_n^2 \end{bmatrix} $


$ x^TSx \to \min $


Ограничения в виде заданного уровня ожидаемой доходности r:


$ \mu^Tx = r $


Доля каждого инструмента в портфеле должна быть положительной:


$ x \ge 0 $


Сумма долей каждого инструмента в портфеле равна единице:


$ |x| = 1 $


Обернем решение задачи условной оптимизации 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 = {'type': 'eq',
                 'fun': lambda x: np.sum(mu @ x.T) - r
                }
    sum_cons = {'type': 'eq',
                 'fun': lambda x: np.sum(x) - 1
                }
    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 за участие в подготовке публикации.

Теги:
Хабы:
+43
Комментарии5
Open Data Science
Крупнейшее русскоязычное Data Science сообщество
Сайт
30
Карма
0
Рейтинг
Андрей Лабинцев @balezz

Программист исследователь

Отправить сообщение

Публикации

Изменить настройки темы

Информация

Сайт
ods.ai
Дата регистрации
Дата основания
Численность
5 001–10 000 человек
Местоположение
Россия

Истории

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr