Работа с постоянными параметрами

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

Рассмотрим функцию потерь с квадратичной ошибкой, которая зависит от некоторых данных x и y, а также параметров betas. Что касается решателя, у функции, которую мы хотим минимизировать, должен быть только один входной аргумент, назовем его sqerror.

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

x = [1.0, 2.0, 3.0]
y = 1.0 .+ 2.0 .* x .+ [-0.3, 0.3, -0.1]

Затем мы просто определяем функцию для трех переменных

function sqerror(betas, X, Y)
    err = 0.0
    for i in 1:length(X)
        pred_i = betas[1] + betas[2] * X[i]
        err += (Y[i] - pred_i)^2
    end
    return err
end

а затем оптимизируем следующую анонимную функцию

res = optimize(b -> sqerror(b, x, y), [0.0, 0.0])

Или можно определить закрытие sqerror(betas), которому известны только что определенные переменные

function sqerror(betas)
    err = 0.0
    for i in 1:length(x)
        pred_i = betas[1] + betas[2] * x[i]
        err += (y[i] - pred_i)^2
    end
    return err
end

Затем можно оптимизировать функцию sqerror так же, как и любую другую функцию

res = optimize(sqerror, [0.0, 0.0])

Избегайте повторяющихся вычислений

Допустим, вы оптимизируете функцию

f(x) = x[1]^2+x[2]^2
g!(storage, x) = copyto!(storage, [2x[1], 2x[2]])

В этой ситуации никакие вычисления из f не могут быть повторно использованы в g!. Однако иногда между целевой функцией и градиентом есть существенное сходство, и некоторые вычисления можно использовать повторно.

Чтобы избежать повторных вычислений, определите функции fg! или fgh!, которые одновременно вычисляют целевую функцию, градиент и гессиан (при необходимости). Эти функции можно написать так, чтобы не повторять общие вычисления.

Например, здесь мы определяем функцию fg! для вычисления целевой функции и градиента, если это необходимо:

function fg!(F, G, x)
  # здесь выполняются обычные вычисления
  # ...
  if G !== nothing
    # здесь приводится код для вычисления градиента
    # запись результата в вектор G
    # G .= ...
  end
  if F !== nothing
    # value = ... код для вычисления целевой функции
    return value
  end
end

Optim будет вызывать эту функцию только с аргументом G, который имеет значение nothing (если градиент не требуется) или Vector, который должен быть заполнен (на месте) градиентом. Такая гибкость удобна для алгоритмов, которые используют градиент только в одних итерациях, но не в других.

Теперь мы вызываем optimize со следующим синтаксисом:

Optim.optimize(Optim.only_fg!(fg!), [0., 0.], Optim.LBFGS())

Аналогично для вычислений, требующих гессиана, мы можем написать следующее:

function fgh!(F, G, H, x)
  G === nothing || # вычисление градиента и сохранение в G
  H === nothing || # вычисление гессиана и сохранение в H
  F === nothing || return f(x)
  nothing
end

Optim.optimize(Optim.only_fgh!(fgh!), [0., 0.], Optim.Newton())

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

Как уже говорилось во введении, передача аналитических градиентов может повлиять на производительность. Чтобы показать соответствующий пример, рассмотрим сепарабельное расширение функции Розенброка в размерности 5000. См. SROSENBR в CUTEst.

Ниже мы используем градиенты и целевые функции из mastsif через CUTEst.jl. Мы показываем только первые пять итераций попытки минимизировать функцию с помощью градиентного спуска.

julia> @time optimize(f, initial_x, GradientDescent(),
                      Optim.Options(show_trace=true, iterations = 5))
Iter     Function value   Gradient norm
     0     4.850000e+04     2.116000e+02
     1     1.018734e+03     2.704951e+01
     2     3.468449e+00     5.721261e-01
     3     2.966899e+00     2.638790e-02
     4     2.511859e+00     5.237768e-01
     5     2.107853e+00     1.020287e-01
 21.731129 seconds (1.61 M allocations: 63.434 MB, 0.03% gc time)
Results of Optimization Algorithm
 * Algorithm: Gradient Descent
 * Starting Point: [1.2,1.0, ...]
 * Minimizer: [1.0287767703731154,1.058769439356144, ...]
 * Minimum: 2.107853e+00
 * Iterations: 5
 * Convergence: false
   * |x - x'| < 0.0: false
   * |f(x) - f(x')| / |f(x)| < 0.0: false
   * |g(x)| < 1.0e-08: false
   * Reached Maximum Number of Iterations: true
 * Objective Function Calls: 23
 * Gradient Calls: 23

julia> @time optimize(f, g!, initial_x, GradientDescent(),
                      Optim.Options(show_trace=true, iterations = 5))
Iter     Function value   Gradient norm
     0     4.850000e+04     2.116000e+02
     1     1.018769e+03     2.704998e+01
     2     3.468488e+00     5.721481e-01
     3     2.966900e+00     2.638792e-02
     4     2.511828e+00     5.237919e-01
     5     2.107802e+00     1.020415e-01
  0.009889 seconds (915 allocations: 270.266 KB)
Results of Optimization Algorithm
 * Algorithm: Gradient Descent
 * Starting Point: [1.2,1.0, ...]
 * Minimizer: [1.0287763814102757,1.05876866832087, ...]
 * Minimum: 2.107802e+00
 * Iterations: 5
 * Convergence: false
   * |x - x'| < 0.0: false
   * |f(x) - f(x')| / |f(x)| < 0.0: false
   * |g(x)| < 1.0e-08: false
   * Reached Maximum Number of Iterations: true
 * Objective Function Calls: 23
 * Gradient Calls: 23

Цель получила значение, которое очень схоже в двух выполнениях, но выполнение с аналитическим градиентом намного быстрее. Возможно, код конечных разностей может быть усовершенствован, но в целом оптимизация будет замедлена из-за всех вычислений функций, необходимых для вычисления центральных конечных разностей.

Разделение времени, затраченного на код Optim и заданные пользователем функции

Рассмотрим задачу Розенброка.

using Optim, OptimTestProblems
prob = UnconstrainedProblems.examples["Rosenbrock"];

Допустим, мы оптимизируем эту функцию и затем хотим узнать общее время выполнения оптимизации (optimize) с помощью метода Ньютона с доверительной областью. Это удивительно, но она занимает много времени. Затем мы задаемся вопросом, тратится ли время в собственном коде Optim (например, на решение подзадачи) или на вычисление цели, градиента или гессиана, которые мы предоставили. Здесь может быть очень полезен пакет TimerOutputs.jl. Он позволяет запустить общий таймер для optimize, а также добавить отдельные таймеры для f, g! и h!. Рассмотрим приведенный ниже пример, который принадлежит автору пакета Кристоферу Карлссону (Kristoffer Carlsson).

using TimerOutputs
const to = TimerOutput()

f(x    ) =  @timeit to "f"  prob.f(x)
g!(x, g) =  @timeit to "g!" prob.g!(x, g)
h!(x, h) =  @timeit to "h!" prob.h!(x, h)

begin
reset_timer!(to)
@timeit to "Trust Region" begin
    res = Optim.optimize(f, g!, h!, prob.initial_x, NewtonTrustRegion())
end
show(to; allocations = false)
end

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

Ранняя остановка

Иногда бывает интересно остановить работу оптимизатора раньше времени. Проще всего это сделать, задав в качестве значения ключевого слова iterations в Optim.Options некоторое число. В этом случае счетчик итераций не сможет превысить определенный предел, стандартное значение которого равно 1000. В качестве альтернативы можно установить мягкое ограничение на время выполнения процедуры оптимизации, задав ключевое слово time_limit в конструкторе Optim.Options.

using Optim, OptimTestProblems
problem = UnconstrainedProblems.examples["Rosenbrock"]

f = problem.f
initial_x = problem.initial_x

function slow(x)
    sleep(0.1)
    f(x)
end

start_time = time()

optimize(slow, zeros(2), NelderMead(), Optim.Options(time_limit = 3.0))

Выполнение остановится примерно через три секунды. Если нам важнее остановиться до достижения предела, можно использовать обратный вызов с простой моделью для предсказания значения времени, которое пройдет до окончания следующей итерации. Рассмотрим следующий код.

using Optim, OptimTestProblems
problem = UnconstrainedProblems.examples["Rosenbrock"]

f = problem.f
initial_x = problem.initial_x

function very_slow(x)
    sleep(.5)
    f(x)
end

start_time = time()
time_to_setup = zeros(1)
function advanced_time_control(x)
    println(" * Iteration:       ", x.iteration)
    so_far =  time()-start_time
    println(" * Time so far:     ", so_far)
    if x.iteration == 0
        time_to_setup .= time()-start_time
    else
        expected_next_time = so_far + (time()-start_time-time_to_setup[1])/(x.iteration)
        println(" * Next iteration ≈ ", expected_next_time)
        println()
        return expected_next_time < 13 ? false : true
    end
    println()
    false
end
optimize(very_slow, zeros(2), NelderMead(), Optim.Options(callback = advanced_time_control))

Он попытается предсказать время, прошедшее после завершения следующей итерации, и остановится, если ожидается, что оно превысит предел в 13 секунд. Запустив его, мы получим примерно следующий результат.

julia> optimize(very_slow, zeros(2), NelderMead(), Optim.Options(callback = advanced_time_control))
 * Iteration:       0
 * Time so far:     2.219298839569092

 * Iteration:       1
 * Time so far:     3.4006409645080566
 * Next iteration ≈ 4.5429909229278564

 * Iteration:       2
 * Time so far:     4.403923988342285
 * Next iteration ≈ 5.476739525794983

 * Iteration:       3
 * Time so far:     5.407265901565552
 * Next iteration ≈ 6.4569235642751055

 * Iteration:       4
 * Time so far:     5.909044027328491
 * Next iteration ≈ 6.821732044219971

 * Iteration:       5
 * Time so far:     6.912338972091675
 * Next iteration ≈ 7.843148183822632

 * Iteration:       6
 * Time so far:     7.9156060218811035
 * Next iteration ≈ 8.85849153995514

 * Iteration:       7
 * Time so far:     8.918903827667236
 * Next iteration ≈ 9.870419979095459

 * Iteration:       8
 * Time so far:     9.922197818756104
 * Next iteration ≈ 10.880185931921005

 * Iteration:       9
 * Time so far:     10.925468921661377
 * Next iteration ≈ 11.888488478130764

 * Iteration:       10
 * Time so far:     11.92870283126831
 * Next iteration ≈ 12.895747828483582

 * Iteration:       11
 * Time so far:     12.932114839553833
 * Next iteration ≈ 13.902462200684981

Results of Optimization Algorithm
 * Algorithm: Nelder-Mead
 * Starting Point: [0.0,0.0]
 * Minimizer: [0.23359374999999996,0.042187499999999996, ...]
 * Minimum: 6.291677e-01
 * Iterations: 11
 * Convergence: false
   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: false
   * Reached Maximum Number of Iterations: false
 * Objective Function Calls: 24