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

Оптимизация

Открыть пример в Engee

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

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

In [ ]: 
# Загрузим необходимые библиотеки
using Plots
plotlyjs()
Out[0]:
Plots.PlotlyJSBackend()

Оптимизация функции нескольких переменных

http://julianlsolvers.github.io/Optim.jl/v0.9.3/user/minimization/

Найдите минимум функции $f(x) = (1.0 - x_1)^2 + 100.0 * (x_2 - x_1^2)^2$ при помощи библиотеки Optim.

Укажите начальную точку $(0,0)$.

In [ ]: 
using Optim
f_optim(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
results_auto = Optim.optimize(f_optim, [0.0, 0.0])
print(results_auto)

 * Status: success

 * Candidate solution
    Final objective value:     3.525527e-09

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    60
    f(x) calls:    117

Реализация градиентного спуска

https://nextjournal.com/DeepLearningNotes/Ch06StochasticGradientDescent

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

In [ ]: 
using Plots, LinearAlgebra, StatsBase

Входные данные для этой задачи создаются в следующей ячейке.

In [ ]: 
regX = rand(100)
regY = 50 .+ 100 * regX + 2 * randn(100);
scatter( regX, regY, fmt = :png, legend=:bottomright, label="data" )
Out[0]:

Точное решение этой задачи можно выразить следующим образом.

In [ ]: 
ϕ = hcat(regX.^0, regX.^1);
δ = 0.01;
θ = inv(ϕ'ϕ + δ^2*I)*ϕ'regY;

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

In [ ]: 
function loss(X,y,θ)
    cost = (y .- X*θ)'*(y .- X*θ);
    grad = - 2 * X'y + 2*X'X*θ;
    return (cost, grad)
end
Out[0]:
loss (generic function with 1 method)
In [ ]: 
function sgd(loss, X, y, θ_start, η, n, num_iters)
    #=
       Перечень аргументов:
       loss -- функция, которую мы оптимизируем, 
               она принимает датасет (X, y) и
               текущий набор параметров θ;
               функция должна возвращать
               значение потерь и градиент при наборе параметров θ
       theta0 -- исходная оценка вектора параметров θ
       η -- шаг обучения (learning rate)
       n -- размер батча
       num_iters -- макс.количество выполнения SGD
       Возвращает:
       theta -- вектор параметров по завершению SGD
       path -- оценки оптимального θ на каждом шаге обучения
    =#
    data = [(X[i,:],y[i]) for i in 1:length(y)]
    θ = θ_start
    path = θ
    for iter in 1:num_iters
        s = StatsBase.sample(data, n)
        Xs = hcat(first.(s)...)'
        ys = last.(s)
        cost, grad = loss(Xs,ys,θ)
        θ = θ .- ( ( η * num_iters/(num_iters+iter) * grad) / n )
        path = hcat(path,θ)
    end
    return (θ,path)
end
Out[0]:
sgd (generic function with 1 method)

Вывод графика с ходом оптимизации:

In [ ]: 
# Поверхность функции потерь
xr = range(40, stop=60, length=10);
yr = range(90, stop=110, length=10);
X = hcat(ones(length(regX)), regX);
f(x,y) = loss(X,regY,[x y]')[1][1];
plot(xr, yr, f, st=:wireframe, fmt = :png, label="Функция потерь", color="blue")

# Пути оптимизации
X = hcat(ones(length(regX)), regX);
y = regY;
(t,p) = sgd(loss, X, y, [40 90]', 0.1, 1, 500);
scatter!( p[1,:], p[2,:], f.(p[1,:],p[2,:]),
    label="Путь оптимизации", ms=3, markerstrokecolor = "white", camera=(40,30) )
Out[0]:

Выведем результаты оптимизации:

In [ ]: 
println("Координаты найденной точки, в которой функция принимает минимальное значение:")
print( p[:,end] )

Координаты найденной точки, в которой функция принимает минимальное значение:
[50.44203686443982, 99.09098067162014]

Поиск оптимального решения с ограничениями

https://www.juliabloggers.com/optimization-with-jump-in-julia/

Предположим, перед нами стоит такая задача оптимизации:

$$ \begin{aligned} & \max 5 x+4 y \\ & \text { условия } \\ & x+y \leq 5 \\ & 10 x+6 y \leq 45 \\ & x, y \geq 0 \end{aligned} $$

Область допустимых значений показаны на рисунке (x и y – непрерывные, независимые переменные):

edu_optim_demo1.png

Нам заданы целевая задача оптимизации (красная стрелка) и направление (максимизация), зелёная точка – оптимальное решение, её координаты: x = 3.75 и y = 1.,25. Оптимальное значение целевой функции равно 23.75.

Решим эту задачу при помощи JuMP – библиотеки Julia для оптимизации.

In [ ]: 
using JuMP, GLPK

# Зададим модель
m = Model(); # m будет значить "модель"

# Выберем оптимизатор
set_optimizer(m, GLPK.Optimizer); # выберем решатель GLPK (GNU Linear Programming Kit)

# Определим
@variable(m, x>=0);
@variable(m, y>=0);

# Определим ограничения
@constraint(m, x+y<=5);
@constraint(m, 10x+6y<=45);

# Определим целевую функцию
@objective(m, Max, 5x+4y);

# Запустим решатель
optimize!(m);

# Вывод
println(objective_value(m)) # Значение целевой функции в оптимальной точке
println("x = ", value.(x), "\n","y = ",value.(y)) # Координаты оптимальных x и y

[ Info: Precompiling JuMP [4076af6c-e467-56ae-b986-b466b2749572]
[ Info: Precompiling GLPK [60bf3e95-4087-53dc-ae20-288a0d20c6a6]

23.75
x = 3.75
y = 1.2500000000000002