Документация Engee

Оптимизация без ограничений

Чтобы показать, как можно использовать пакет Optim, мы минимизируем функцию Розенброка — классическую тестовую задачу для численной оптимизации. Будем считать, что вы уже установили пакет Optim с помощью диспетчера пакетов Julia. Сначала загрузим Optim и определим функцию Розенброка:

using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

После определения этой функции можно найти минимизатор (входное значение, которое минимизирует цель) и минимум (значение цели в минимизаторе) с помощью любого предпочитаемого алгоритма оптимизации. Определив функцию, мы просто указываем начальную точку x и вызываем optimize с отправной точкой x0:

x0 = [0.0, 0.0]
optimize(f, x0)

Примечание. Важно передавать initial_x в виде массива. Если ваша задача одномерная, ее придется заключить в массив. Это можно легко сделать, написав optimize(x->f(first(x)), [initial_x]), чтобы входные данные были массивом, но анонимная функция автоматически передавала первый (и единственный) элемент в вашу функцию f.

Optim по умолчанию будет использовать метод Нелдера-Мида в многомерном случае, поскольку мы не задали градиент. Это также можно указать явно, используя следующее.

optimize(f, x0, NelderMead())

Также доступны другие решатели. Ниже мы используем L-BFGS, квазиньютоновский метод, требующий градиента. Если мы передадим только f, Optim построит приблизительный градиент с помощью центрального конечного дифференцирования:

optimize(f, x0, LBFGS())

Для повышения производительности и точности можно передать собственную функцию градиента. Если цель написана полностью в коде Julia без специальных обращений к внешним (то есть не принадлежащим Julia) библиотекам, вы также можете применять автоматическое дифференцирование, используя ключевое слово autodiff и задав ему значение :forward:

optimize(f, x0, LBFGS(); autodiff = :forward)

Для примера Розенброка можно показать, что аналитический градиент имеет следующий вид:

function g!(G, x)
    G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    G[2] = 200.0 * (x[2] - x[1]^2)
end

Обратите внимание, что функции, которые мы используем для вычисления градиента (а затем и гессиана h!) функции Розенброка, изменяют массив хранения фиксированного размера, который в этих примерах передается в качестве дополнительного аргумента G (или H для гессиана). Изменяя один массив в течение многих итераций, этот стиль определения функции устраняет иногда значительные затраты, связанные с выделением нового массива при каждом вызове функций g! или h!. Если вы предпочитаете, чтобы ваши градиенты просто принимали x, вы все равно можете использовать optimize, задав ключевому слову inplace значение false:

optimize(f, g, x0; inplace = false)

где g является функцией только x.

В нашей версии на месте вы просто передаете g! вместе с f из предыдущей версии, чтобы использовать градиент:

optimize(f, g!, x0, LBFGS())

Для некоторых методов, таких как имитация отжига, градиент игнорируется:

optimize(f, g!, x0, SimulatedAnnealing())

Помимо указания градиентов можно также предоставить функцию гессиана h!. В нашем случае она выглядит так:

function h!(H, x)
    H[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
    H[1, 2] = -400.0 * x[1]
    H[2, 1] = -400.0 * x[1]
    H[2, 2] = 200.0
end

Теперь мы можем использовать метод Ньютона для оптимизации, выполнив следующее.

optimize(f, g!, h!, x0)

По умолчанию используется Newton(), поскольку была задана функция гессиана. Как и градиенты, функция гессиана будет игнорироваться, если вы используете метод, которому она не требуется:

optimize(f, g!, h!, x0, LBFGS())

Обратите внимание, что Optim не будет генерировать приближенные гессианы с помощью конечного дифференцирования из-за потенциально низкой точности приближений к гессианам. Кроме метода Ньютона, ни один из алгоритмов, предоставляемых пакетом Optim, не использует точные гессианы.

Оптимизация с блочными ограничениями

Доступен прямой алгоритм внутренней точки для простых блочных ограничений (нижние и верхние границы). Еще раз возьмем пример Розенброка, приведенный выше; блочная минимизация выполняется следующим образом.

lower = [1.25, -2.1]
upper = [Inf, Inf]
initial_x = [2.0, 2.0]
inner_optimizer = GradientDescent()
results = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))

Выполняется оптимизация с барьерным штрафом, последовательно уменьшая барьерный коэффициент и используя выбранный inner_optimizer (GradientDescent() выше) для сходимости на каждом шаге. Чтобы изменить специфические параметры алгоритма, например алгоритм линейного поиска, укажите его непосредственно в конструкторе inner_optimizer:

lower = [1.25, -2.1]
upper = [Inf, Inf]
initial_x = [2.0, 2.0]
# требует использования линейного поиска
inner_optimizer = GradientDescent(linesearch=LineSearches.BackTracking(order=3))
results = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))

Этот алгоритм использует диагональное предобусловливание для повышения точности и, следовательно, является хорошим примером использования ConjugateGradient или LBFGS с предобусловливанием. Другие методы в настоящее время не используют предобусловливание. Используются только блочные ограничения. Если вы можете аналитически вычислить диагональ гессиана целевой функции, возможно, вам стоит подумать о написании собственного предобусловливателя.

Существует два параметра итераций: параметр внешних итераций, используемый для управления Fminbox, и параметр внутренних итераций, используемый для управления внутренним оптимизатором. Например, ниже приведено ограничение оптимизации до двух основных итераций.

results = optimize(f, g!, lower, upper, initial_x, Fminbox(GradientDescent()), Optim.Options(outer_iterations = 2))

В отличие от этого, в следующем примере максимальное количество итераций для каждой оптимизации GradientDescent() равно двум.

results = optimize(f, g!, lower, upper, initial_x, Fminbox(GradientDescent()), Optim.Options(iterations = 2))

Использование информации второго порядка

При наличии гессиана целевой функции можно использовать прямо-двойственный алгоритм, реализованный в IPNewton. Интерфейс аналогичен:

results = optimize(f, lower, upper, initial_x, IPNewton())
results = optimize(f, g!, lower, upper, initial_x, IPNewton())
results = optimize(f, g!, h!, lower, upper, initial_x, IPNewton())

Минимизация одномерной функции в ограниченном интервале

Минимизация одномерных функций без производных выполняется через интерфейс optimize:

optimize(f, lower, upper, method; kwargs...)

Обратите внимание на отсутствие начального x. В качестве конкретного примера можно привести следующую квадратичную функцию.

julia> f_univariate(x) = 2x^2+3x+1
f_univariate (generic function with 1 method)

julia> optimize(f_univariate, -2.0, 1.0)
Results of Optimization Algorithm
 * Algorithm: Brent's Method
 * Search Interval: [-2.000000, 1.000000]
 * Minimizer: -7.500000e-01
 * Minimum: -1.250000e-01
 * Iterations: 7
 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
 * Objective Function Calls: 8

Результат показывает, что мы задали начальные нижнюю и верхнюю границы, что есть конечный минимизатор и минимум, и что для этого потребовалось семь основных итераций. Важно отметить, что мы также видим, что была объявлена сходимость. По умолчанию используется метод Брента, который является одним из двух доступных методов:

  • Метод Брента, используемый по умолчанию (может быть выбран явным образом с помощью Brent()).

  • Поиск золотого сечения; доступен при использовании GoldenSection().

Чтобы указать этот метод вручную, используется обычный синтаксис, как для многомерной оптимизации.

    optimize(f, lower, upper, Brent(); kwargs...)
    optimize(f, lower, upper, GoldenSection(); kwargs...)

Для настройки параметров этого особого типа оптимизации применяются ключевые слова. Помимо параметров iterations, store_trace, show_trace, show_warnings и extended_trace, также доступны следующие.

  • rel_tol: относительный допуск, используемый для определения сходимости. Значение по умолчанию — sqrt(eps(T)).

  • abs_tol: абсолютный допуск, используемый для определения сходимости. Значение по умолчанию — eps(T).

Получение результатов

После того как мы получили результаты в res, можно использовать API для получения результатов оптимизации. Он состоит из набора функций. Они не экспортируются, поэтому их нужно предварять префиксом Optim.. Допустим, мы выполняем следующую оптимизацию.

res = optimize(x->dot(x,[1 0. 0; 0 3 0; 0 0 1]*x), zeros(3))

Если мы не можем вспомнить, какой метод использовали, просто применяем

summary(res)

и в результате возвращается "Nelder Mead". Полезной для нас информацией являются минимизатор и минимум целевых функций, которые можно найти следующим образом.

julia> Optim.minimizer(res)
3-element Array{Float64,1}:
 -0.499921
 -0.3333
 -1.49994

julia> Optim.minimum(res)
 -2.8333333205768865

Полный список функций

Полный список функций можно найти ниже.

Определены для всех методов:

  • summary(res)

  • minimizer(res)

  • minimum(res)

  • iterations(res)

  • iteration_limit_reached(res)

  • trace(res)

  • x_trace(res)

  • f_trace(res)

  • f_calls(res)

  • converged(res)

Определены для одномерной оптимизации:

  • lower_bound(res)

  • upper_bound(res)

  • x_lower_trace(res)

  • x_upper_trace(res)

  • rel_tol(res)

  • abs_tol(res)

Определены для многомерной оптимизации:

  • g_norm_trace(res)

  • g_calls(res)

  • x_converged(res)

  • f_converged(res)

  • g_converged(res)

  • initial_state(res)

Определены для NelderMead с параметром trace_simplex=true:

  • centroid_trace(res) (с extended_trace=true)

  • simplex_trace(res)

  • simplex_values_trace(res)

Входные типы

Большинство пользователей вводят Vector в качестве initial_x и получают Optim.minimizer(res), который также является вектором. Для методов нулевого и первого порядка также можно передавать матрицы или даже массивы более высокой размерности. Единственное ограничение, налагаемое при наличии случая с Vector, заключается в том, что больше нельзя использовать аппроксимацию конечными разностями или автоматическое дифференцирование. Методы второго порядка (варианты метода Ньютона) не поддерживают этот более общий тип входных данных.

Примечания о флагах и проверках сходимости

В настоящее время можно получить доступ к минимизатору, используя Optim.minimizer(result), даже если все флаги сходимости имеют значение false. Это означает, что пользователю следует с осторожностью использовать выходные данные решателей. Рекомендуется включить проверки на сходимость, если минимизатор или минимум используются для дальнейших вычислений.

Следует отметить, что методы первого и второго порядков проверяют сходимость градиента перед входом в цикл оптимизации. Это делается для того, чтобы избежать ошибок при линейном поиске, если initial_x является стационарной точкой. Обратите внимание, что это только проверка первого порядка. Если initial_x представляет собой любой тип стационарной точки, g_converged будет иметь значение true. Сюда входят локальные минимумы, седловые точки и локальные максимумы. Если iterations имеет значение 0, а g_converged имеет значение true, пользователь должен помнить об этом.