scipy.optimize.minimize ('SLSQP') слишком медленно при заданной переменной dim 2000

На мой взгляд, scipy.minimze предоставляет интуитивно понятный интерфейс для оптимизации. Я обнаружил, что ускорение функции стоимости (и, в конечном итоге, функции градиента) может дать вам хорошее ускорение.

Например, возьмем N-мерную функцию Розенброка:

import numpy as np
from scipy.optimize import minimize


def rosenbrock(x, N):
    out = 0.0
    for i in range(N-1):
        out += 100.0 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2
    return out

# slow optimize
N = 20
x_0 = - np.ones(N)
%timeit minimize(rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
res = minimize(rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
print(res.message)

Сроки оптимизации доходности

102 ms ± 1.86 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Optimization terminated successfully.

Теперь вы можете ускорить целевую функцию с помощью numba и предоставить простую функцию для вычисления градиента следующим образом:

from numba import jit, float64, int64

@jit(float64(float64[:], int64), nopython=True, parallel=True)
def fast_rosenbrock(x, N):
    out = 0.0
    for i in range(N-1):
        out += 100.0 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2
    return out


@jit(float64[:](float64[:], int64), nopython=True, parallel=True)
def fast_jac(x, N):
    h = 1e-9
    jac = np.zeros_like(x)
    f_0 = fast_rosenbrock(x, N)
    for i in range(N):
        x_d = np.copy(x)
        x_d[i] += h
        f_d = fast_rosenbrock(x_d, N)
        jac[i] = (f_d - f_0) / h
    return jac

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

print('with fast jacobian')
%timeit minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4}, jac=fast_jac)
print('without fast jacobian')
%timeit minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
res = minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4}, jac=fast_jac)
print(res.message)

Испытание того и другого, с функцией быстрого якобиана и без нее. Результатом этого является:

with fast jacobian
9.67 ms ± 488 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
without fast jacobian
27.2 ms ± 2.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Optimization terminated successfully.

Это примерно в 10 раз ускорение без особых усилий. Улучшение, которого вы можете достичь с помощью этого, во многом зависит от неэффективности вашей функции затрат. У меня была функция стоимости с несколькими вычислениями, и я смог получить ускорение примерно на 10 ^ 2 - 10 ^ 3.

Преимущество этого подхода в том, что он не требует больших усилий, и вы можете остаться с scipy и его приятным интерфейсом.

Мы мало что знаем о модели, но вот несколько примечаний:

1. SLSQP действительно разработан для небольших (плотных), хорошо масштабируемых моделей.
2. SLSQP - это локальный решатель. Он принимает невыпуклые проблемы, но предоставляет только локальные решения.
3. Я сомневаюсь, что для SLSQP есть такая сложность. Во всяком случае, это мало что говорит о производительности той или иной задачи.
4. IPOPT - это крупномасштабный решатель разреженных внутренних точек. Он найдет локальные решения. Он может решать действительно большие модели.
5. Глобальные решатели, такие как BARON, ANTIGONE и COUENNE, находят глобальные решения (или ограничивают качество решения, если вы нетерпеливы и преждевременно останавливаетесь). Эти решатели (большую часть времени) медленнее, чем локальные решатели. Мне не известны прямые ссылки на Python.
6. Если у вас есть хорошая отправная точка, локальный решатель может быть именно тем, что вам нужно. Использование многозадачной стратегии может помочь найти лучшие решения (не доказано, что оптимальность в глобальном масштабе не доказана, но вы получаете некоторую уверенность в том, что не нашли действительно плохого локального оптимума).
7. Фреймворк Python PYOMO предлагает доступ ко многим решателям. Однако вам нужно будет переписать модель. PYOMO имеет автоматическую дифференциацию: нет необходимости предоставлять градиенты.
8. Для тестирования вы можете попробовать расшифровать модель в AMPL или GAMS и решить онлайн через NEOS. Это позволит вам опробовать несколько решателей. И AMPL, и GAMS имеют автоматическое различие.

Удивительно, но я нашел относительно нормальное решение, используя оптимизатор для платформы глубокого обучения, Tensorflow, с использованием базового градиентного спуска (на самом деле RMSProp, градиентного спуска с импульсом) после того, как я изменил функцию стоимости, включив ограничение неравенства и ограничивающие ограничения в качестве штрафов. (Полагаю, это то же самое, что и метод лагранжа). Он очень быстро тренируется и быстро сходится с правильными лямбда-параметрами для штрафных ограничений. Мне даже не пришлось переписывать якобианов, поскольку TF позаботится об этом без особого влияния на скорость, очевидно.

До этого мне удалось заставить работать NLOPT, и он намного быстрее, чем scipy / SLSQP, но все же медленнее на более высоких измерениях. Также NLOPT / AUGLANG очень быстрый, но плохо сходится.

Тем не менее, при 20k переменных он все еще медленный. Частично из-за подкачки памяти и функции стоимости, составляющей не менее O (n ^ 2) от попарного евклидова расстояния (я использую (x-x.t) ^ 2 + (y-y.t) ^ 2 с широковещательной передачей). Так что все еще не оптимально.