Engee 文档
Notebook

优化设计

在这个例子中,我们将演示几个简单优化问题的解决方案:使用包查找双参数函数的最小值 Optim,梯度下降的显式实现和线性规划问题中具有约束的最优解的搜索。

让我们做一个加载库的准备单元。

In [ ]:
Pkg.add(["GLPK", "LinearAlgebra", "StatsBase", "JuMP", "Optim"])
In [ ]:
# Загрузим необходимые библиотеки
using Plots
plotlyjs()
Out[0]:
Plots.PlotlyJSBackend()

几个变量函数的优化

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

求函数的最小值 使用图书馆 Optim.

指定起点 .

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]

用约束条件找到最优解

假设我们面临这样的优化问题:

可接受值的范围如图所示(x和y是连续的,自变量):

edu_optim_demo1.png

我们给出了一个优化目标(红色箭头)和一个方向(最大化),绿点是最优解,其坐标为x=3.75和y=1.25。 目标函数的最优值为23.75。

让我们使用JuMPJulia进行优化来解决这个问题。

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