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

Начало работы с Optimization.jl

Страница в процессе перевода.

В этом руководстве мы познакомим вас с основами Optimization.jl, показав, как легко комбинировать локальные и глобальные оптимизаторы для уравнения Розенброка.

Уравнение Розенброка определяется следующим образом:

Это параметризованная задача оптимизации, в которой мы хотим найти вектор u такой, что u минимизирует f. Простейший код, который можно скопировать и вставить, использующий квазиньютоновский метод (LBFGS) для решения задачи Розенброка, выглядит следующим образом:

# Импортируйте пакет и сформулируйте задачу для оптимизации.
using Optimization, OptimizationLBFGSB, Zygote
rosenbrock(u, p) = (p[1] - u[1])^2 + p[2] * (u[2] - u[1]^2)^2
u0 = zeros(2)
p = [1.0, 100.0]

optf = OptimizationFunction(rosenbrock, AutoZygote())
prob = OptimizationProblem(optf, u0, p)

sol = solve(prob, OptimizationLBFGSB.LBFGSB())
retcode: Success
u: 2-element Vector{Float64}:
 0.9999997057368228
 0.999999398151528
sol.u
2-element Vector{Float64}:
 0.9999997057368228
 0.999999398151528
sol.objective
1.0433892998247468e-13

Вот и всё! Так и нужно делать. Теперь давайте подробнее разберёмся, что означает каждая часть и как всё это настроить под ваши нужды.

Понимание объекта решения

Объект решения представляет собой SciMLBase.AbstractNoTimeSolution, и, следовательно, он соответствует Интерфейсу решений SciMLBase для объектов, не являющихся временными рядами и описан на странице типа решения. Однако для простоты давайте немного покажем его в действии.

Оптимизационное решение имеет интерфейс массива, так что оно действует как массив, для которого оно решается. Этот синтаксис массива — это сокращенная запись для простого получения решения u. Например:

sol[1] == sol.u[1]
true
Array(sol) == sol.u
true

sol.objective возвращает итоговую стоимость оптимизации. Мы можем проверить это, подставив её в нашу функцию:

rosenbrock(sol.u, p)
1.0433892998247468e-13
sol.objective
1.0433892998247468e-13

Параметр sol.retcode предоставляет нам дополнительную информацию о процессе решения.

sol.retcode
ReturnCode.Success = 1

Здесь указано ReturnCode.Success, что означает успешное решение. Подробнее о различных кодах возврата можно узнать в разделе ReturnCode документации SciMLBase по адресу: https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes.

Если вас интересует статистика процесса решения, например, для выбора более эффективного решателя, вы можете изучить sol.stats.

sol.stats
SciMLBase.OptimizationStats
Number of iterations:                              21
Time in seconds:                                   0.098339
Number of function evaluations:                    25
Number of gradient evaluations:                    25
Number of hessian evaluations:                     0

Это лишь малая часть того, что есть, посмотрите другие страницы для получения более подробной информации, а теперь перейдем к настройке.

Импорт другого пакета решателя и решение задачи

OptimizationOptimJL — это обертка для Optim.jl, а OptimizationBBO — это обертка для BlackBoxOptim.jl.

Сначала давайте воспользуемся решателем NelderMead, не использующим производные, из пакета Optim.jl:

using OptimizationOptimJL
sol = solve(prob, Optim.NelderMead())
retcode: Success
u: 2-element Vector{Float64}:
 0.9999634355313174
 0.9999315506115275

BlackBoxOptim.jl предлагает решатели глобальной оптимизации без использования производных, которые требуют установки границ с помощью параметров lb и ub в OptimizationProblem. Давайте воспользуемся решателем BBO_adaptive_de_rand_1_bin_radiuslimited():

using OptimizationBBO
prob = OptimizationProblem(rosenbrock, u0, p, lb = [-1.0, -1.0], ub = [1.0, 1.0])
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited())
retcode: MaxIters
u: 2-element Vector{Float64}:
 1.0
 0.9999999999999999

Решение, полученное с помощью исходного решателя, всегда можно получить через original:

sol.original
BlackBoxOptim.OptimizationResults(\
adaptive_de_rand_1_bin_radiuslimited\", \"Max number of steps (10000) reached\", 10001, 1.771283167061867e9, 0.008346080780029297, BlackBoxOptim.ParamsDictChain[BlackBoxOptim.ParamsDictChain[Dict{Symbol, Any}(:RngSeed => 110509, :SearchRange => [(-1.0, 1.0), (-1.0, 1.0)], :TraceMode => :silent, :Method => :adaptive_de_rand_1_bin_radiuslimited, :MaxSteps => 10000),Dict{Symbol, Any}()],Dict{Symbol, Any}(:CallbackInterval => -1.0, :TargetFitness => nothing, :TraceMode => :compact, :FitnessScheme => BlackBoxOptim.ScalarFitnessScheme{true}(), :MinDeltaFitnessTolerance => 1.0e-50, :NumDimensions => :NotSpecified, :FitnessTolerance => 1.0e-8, :TraceInterval => 0.5, :MaxStepsWithoutProgress => 10000, :MaxSteps => 10000…)], 10123, BlackBoxOptim.ScalarFitnessScheme{true}(), BlackBoxOptim.TopListArchiveOutput{Float64, Vector{Float64}}(1.232595164407831e-30, [1.0, 0.9999999999999999]), BlackBoxOptim.PopulationOptimizerOutput{BlackBoxOptim.FitPopulation{Float64}}(BlackBoxOptim.FitPopulation{Float64}([0.9999999999999999 0.9999999999999993 … 0.9999999999998591 0.9999999999999493; 0.9999999999999999 0.9999999999999951 … 0.9999999999997199 0.999999999999901], NaN, [1.2449211160519093e-30, 1.2626211826128057e-27, 2.609403963051378e-27, 1.9811378818010626e-27, 1.232595164407831e-30, 6.842986248291311e-27, 2.6437933681383566e-28, 1.1252410559685665e-26, 7.13968423031592e-28, 1.2795601216596803e-25  …  8.019441633341023e-26, 3.2024991738804806e-26, 7.435014031708036e-27, 1.2390500680686085e-24, 1.6884362321444483e-24, 1.3050866744765007e-25, 3.922271887541269e-26, 5.1150023358823055e-25, 2.0126565727501352e-26, 3.170838734487501e-27], 0, BlackBoxOptim.Candidate{Float64}[BlackBoxOptim.Candidate{Float64}([0.9999999999995454, 0.9999999999990971], 34, 2.106988929575013e-25, BlackBoxOptim.AdaptiveDiffEvoRandBin{3}(BlackBoxOptim.AdaptiveDiffEvoParameters(BlackBoxOptim.BimodalCauchy(Distributions.Cauchy{Float64}(μ=0.65, σ=0.1), Distributions.Cauchy{Float64}(μ=1.0, σ=0.1), 0.5, false, true), BlackBoxOptim.BimodalCauchy(Distributions.Cauchy{Float64}(μ=0.1, σ=0.1), Distributions.Cauchy{Float64}(μ=0.95, σ=0.1), 0.5, false, true), [0.24442092887615585, 0.6392413926342186, 0.7200072945856475, 0.574071962668566, 0.671210303023849, 0.8492928641112012, 1.0, 0.6734779426152385, 0.6610270839541772, 1.0  …  0.11061016201125584, 1.0, 0.6540961114128819, 1.0, 0.5390726083867172, 1.0, 0.39497618633170134, 1.0, 0.6004812648321092, 0.8903051110198892], [0.8250017743361556, 0.8654992885751132, 0.9988005853932467, 0.8756522408845833, 0.11739514186337224, 0.973598273818434, 0.1038237373338912, 0.060834820318409685, 0.2331410683532765, 1.0  …  1.0, 0.9444848826061091, 1.0, 0.868155498857225, 1.0, 1.0, 0.0878024832073188, 0.8218291819391484, 0.7927084972576752, 1.0])), 0), BlackBoxOptim.Candidate{Float64}([0.9999999999995454, 0.9999999999996151], 34, 2.770186049042427e-23, BlackBoxOptim.AdaptiveDiffEvoRandBin{3}(BlackBoxOptim.AdaptiveDiffEvoParameters(BlackBoxOptim.BimodalCauchy(Distributions.Cauchy{Float64}(μ=0.65, σ=0.1), Distributions.Cauchy{Float64}(μ=1.0, σ=0.1), 0.5, false, true), BlackBoxOptim.BimodalCauchy(Distributions.Cauchy{Float64}(μ=0.1, σ=0.1), Distributions.Cauchy{Float64}(μ=0.95, σ=0.1), 0.5, false, true), [0.24442092887615585, 0.6392413926342186, 0.7200072945856475, 0.574071962668566, 0.671210303023849, 0.8492928641112012, 1.0, 0.6734779426152385, 0.6610270839541772, 1.0  …  0.11061016201125584, 1.0, 0.6540961114128819, 1.0, 0.5390726083867172, 1.0, 0.39497618633170134, 1.0, 0.6004812648321092, 0.8903051110198892], [0.8250017743361556, 0.8654992885751132, 0.9988005853932467, 0.8756522408845833, 0.11739514186337224, 0.973598273818434, 0.1038237373338912, 0.060834820318409685, 0.2331410683532765, 1.0  …  1.0, 0.9444848826061091, 1.0, 0.868155498857225, 1.0, 1.0, 0.0878024832073188, 0.8218291819391484, 0.7927084972576752, 1.0])), 0)], Base.Threads.SpinLock(0))))"

Определение целевой функции

Optimization.jl предполагает, что ваша целевая функция принимает два аргумента: objective(x, p)

  1. Переменные оптимизации x.

2Другие параметры p, такие как гиперпараметры функции стоимости. Если у вас нет «других параметров», вы можете смело игнорировать этот аргумент. Если ваша целевая функция определена кем-то другим, вы можете создать анонимную функцию, которая просто отбрасывает дополнительные параметры, следующим образом:

obj = (x, p) -> objective(x) # Передайте эту функцию в функцию OptimizationFunction

Управление вычислениями градиентов (автоматическое дифференцирование)

Обратите внимание, что оба описанных выше метода не требовали использования производных, поэтому для оптимизации градиенты не были необходимы. Однако часто оптимизация первого порядка (т.е. с использованием градиентов) гораздо эффективнее. Определение градиентов можно выполнить двумя способами. Один из способов — вручную указать определение градиента в конструкторе OptimizationFunction. Однако более удобный способ получения градиентов — указать тип бэкенда автоматического дифференцирования.

Например, давайте теперь используем метод BFGS библиотеки OptimizationOptimJL для решения той же задачи. Мы импортируем библиотеку автоматического дифференцирования в прямом режиме (using ForwardDiff), а затем укажем в OptimizationFunction, чтобы автоматически строить функции производных с помощью ForwardDiff.jl. Это выглядит так:

using ForwardDiff
optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())
prob = OptimizationProblem(optf, u0, p)
sol = solve(prob, OptimizationOptimJL.BFGS())
retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999373603
 0.99999999986862

Мы можем просмотреть исходный код, чтобы увидеть статистику по количеству необходимых шагов и вычисленным градиентам:

sol.original
 * Status: success

 * Candidate solution
    Final objective value:     7.645553e-21

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 3.48e-07 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.48e-07 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.91e-14 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 9.03e+06 ≰ 0.0e+00
    |g(x)|                 = 2.31e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    16
    f(x) calls:    53
    ∇f(x) calls:   53

Действительно, это намного меньше, чем при использовании методов без производных!

Однако вычислительная стоимость прямого автоматического дифференцирования масштабируется в зависимости от количества входных данных, и, следовательно, по мере роста нашей задачи оптимизации она замедляется. Чтобы компенсировать это, для более крупных задач оптимизации (>100 переменных состояния) обычно рекомендуется использовать обратное автоматическое дифференцирование. Одним из распространенных вариантов обратного автоматического дифференцирования является Zygote.jl. Мы можем продемонстрировать это следующим образом:

using Zygote
optf = OptimizationFunction(rosenbrock, Optimization.AutoZygote())
prob = OptimizationProblem(optf, u0, p)
sol = solve(prob, OptimizationOptimJL.BFGS())
retcode: Success
u: 2-element Vector{Float64}:
 0.9999999999373603
 0.99999999986862

Установка ограничений на ограничивающую рамку

Во многих случаях потенциальные границы значений решения известны. В файле Optimization.jl их можно задать в качестве аргументов lb и ub для нижней и верхней границ соответственно, предоставив вектор значений, по одному на каждую переменную состояния. Давайте теперь выполним нашу градиентную оптимизацию с ограничениями на ограничивающую рамку, перестроив задачу OptimizationProblem:

prob = OptimizationProblem(optf, u0, p, lb = [-1.0, -1.0], ub = [1.0, 1.0])
sol = solve(prob, OptimizationOptimJL.BFGS())
retcode: Success
u: 2-element Vector{Float64}:
 0.9999999992272249
 0.9999999984496138

Для получения дополнительной информации об обработке ограничений, в частности ограничений равенства и неравенства, ознакомьтесь с руководством по ограничениям (@ref constraints).