Julia e o método de descida de coordenadas

O método de descida por coordenadas é um dos métodos mais simples de otimização multidimensional e faz um bom trabalho em encontrar um mínimo local de funções com uma topografia relativamente suave; portanto, é melhor começar a se familiarizar com os métodos de otimização a partir dele.

A busca pelo extremo é realizada na direção dos eixos das coordenadas, ou seja, durante o processo de pesquisa, apenas uma coordenada é alterada. Assim, um problema multidimensional se reduz a um problema unidimensional.

Algoritmo

Primeiro ciclo:

• $$$x_1 = var$, $$$x_2 = x_2 ^ 0$, $$$x_3 = x_3 ^ 0$, ..., $$$x_n = x_n ^ 0$.
• Procurando funções extremas $$$f (x_1)$. Coloque o extremo da função no ponto $$$x_1 ^ 1$.
• $$$x_1 = x_1 ^ 1$, $$$x_2 = var$, $$$x_3 = x_3 ^ 0$, ..., $$$x_n = x_n ^ 0$. Função Extremum $$$f (x_2)$é igual a $$$x_2 ^ 1$
• ...
• $$$x_1 = x_1 ^ 1$, $$$x_2 = x_2 ^ 1$, $$$x_3 = x_3 ^ 1$, ..., $$$x_n = var$.
  Como resultado da execução de n etapas, foi encontrado um novo ponto de abordagem ao extremo $$$(x_1 ^ 1, x_2 ^ 1, ..., x_n ^ 1)$. Em seguida, verificamos o critério para o final da conta: se for atendida, a solução será encontrada, caso contrário, realizaremos outro ciclo.

Preparação

Verifique as versões do pacote:

(v1.1) pkg> status Status `C:\Users\User\.julia\environments\v1.1\Project.toml` [336ed68f] CSV v0.4.3 [a93c6f00] DataFrames v0.17.1 [7073ff75] IJulia v1.16.0 [47be7bcc] ORCA v0.2.1 [58dd65bb] Plotly v0.2.0 [f0f68f2c] PlotlyJS v0.12.3 [91a5bcdd] Plots v0.23.0 [ce6b1742] RDatasets v0.6.1 [90137ffa] StaticArrays v0.10.2 [8bb1440f] DelimitedFiles [10745b16] Statistics 

Definimos uma função para desenhar as linhas de superfície ou de nível nas quais seria conveniente ajustar os limites do gráfico:

 using Plots plotly() #   function plotter(plot_fun; low, up) Xs = range(low[1], stop = up[1], length = 80) Ys = range(low[2], stop = up[2], length = 80) Zs = [ fun([xy]) for x in Xs, y in Ys ]; plot_fun(Xs, Ys, Zs) xaxis!( (low[1], up[1]), low[1]:(up[1]-low[1])/5:up[1] ) #   yaxis!( (low[2], up[2]), low[2]:(up[2]-low[2])/5:up[2] ) end 

Como função de modelo, escolhemos um parabolóide elíptico

 parabol(x) = sum(u->u*u, x) fun = parabol plotter(surface, low = [-1 -1], up = [1 1]) 

Descida de coordenadas

Implementamos o método em uma função que leva o nome do método de otimização unidimensional, a dimensão do problema, o erro desejado, a aproximação inicial e as restrições para desenhar linhas de nível. Todos os parâmetros são definidos para os valores padrão.

 function ofDescent(odm; ndimes = 2, ε = 1e-4, fit = [.5 .5], low = [-1 -1], up = [1 1]) k = 1 #      sumx = 0 sumz = 1 plotter(contour, low = low, up = up) #    x = [ fit[1] ] #     y = [ fit[2] ] #      while abs(sumz - sumx) > ε && k<80 fitz = copy(fit) for i = 1:ndimes odm(i, fit, ε) #    end push!(x, fit[1]) push!(y, fit[2]) sumx = sum(abs, fit ) sumz = sum(abs, fitz) #println("$k $fit") k += 1 end #     scatter!(x', y', legend = false, marker=(10, 0.5, :orange) );# , ,  p = title!("Age = $(size(x,1)) f(x,y) = $(fun(fit))") end 

Em seguida, tentaremos vários métodos de otimização unidimensional.

Método de Newton

$$

$$x_ {n + 1} = x_n - \ frac {f (x_n)} {f '(x_n)}$$

A ideia do método é simples, assim como a implementação

 #  dfun = (x, i) -> i == 1 ? 2x[1] + x[2]*x[2] : 2x[2] + x[1]*x[1] function newton2(i, fit, ϵ) k = 1 oldfit = Inf while ( abs(fit[i] - oldfit) > ϵ && k<50 ) oldfit = fit[i] tryfit = fun(fit) / dfun(fit, i) fit[i] -= tryfit println(" $k $fit") k+=1 end end ofDescent(newton2, fit = [9.1 9.1], low = [-4 -4], up = [13 13]) 

Newton é bastante exigente na aproximação inicial e, sem limitação nos degraus, pode facilmente partir para distâncias desconhecidas. O cálculo da derivada é desejável, mas pequenas variações podem ser dispensadas. Modificamos nossa função:

 function newton(i, fit, ϵ) k = 1 oldfit = Inf x = [] y = [] push!(x, fit[1]) push!(y, fit[2]) while ( abs(oldfit - fit[i]) > ϵ && k<50 ) fx = fun(fit) oldfit = fit[i] fit[i] += 0.01 dfx = fun(fit) fit[i] -= 0.01 tryfit = fx*0.01 / (dfx-fx) #        if( abs(tryfit) > abs(fit[i])*10 ) push!(x, fit[1]) push!(y, fit[2]) break end fit[i] -= tryfit #println(" $k $fit") push!(x, fit[1]) push!(y, fit[2]) k+=1 end #  - -   plot!(x, y, w = 3, legend = false, marker = :rect ) end ofDescent(newton, fit = [9.1 9.1], low = [-4 -4], up = [13 13]) 

Interpolação parabólica inversa

Um método que não requer conhecimento da derivada e tem boa convergência

 function ipi(i, fit, ϵ) # inverse parabolic interpolation n = 0 xn2 = copy(fit) xn1 = copy(fit) xn = copy(fit) xnp = zeros( length(fit) ) xy = copy(fit) xn2[i] = fit[i] - 0.15 xn[i] = fit[i] + 0.15 fn2 = fun(xn2) fn1 = fun(xn1) while abs(xn[i] - xn1[i]) > ϵ && n<80 fn = fun(xn) a = fn1*fn / ( (fn2-fn1)*(fn2-fn ) ) b = fn2*fn / ( (fn1-fn2)*(fn1-fn ) ) c = fn2*fn1 / ( (fn -fn2)*(fn -fn1) ) xnp[i] = a*xn2[i] + b*xn1[i] + c*xn[i] xn2[i] = xn1[i] xn1[i] = xn[i] xn[i] = xnp[i] fn2 = fn1 fn1 = fn n += 1 println(" $n $xn $xn1") xy = [xy; xn] end fit[i] = xnp[i] plot!(xy[:,1], xy[:,2], w = 3, legend = false, marker = :rect ) end ofDescent(ipi, fit = [0.1 0.1], low = [-.1 -.1], up = [.4 .4]) 

Se piorarmos a aproximação inicial, o método começará a exigir um imenso número de etapas para cada era de descida de coordenadas. Nesse sentido, seu irmão ganha

Interpolação Parabólica Sequencial

Também requer três pontos de partida , mas em muitas funções de teste mostra resultados mais satisfatórios.

 function spi(i, fit, ϵ) # sequential parabolic interpolation n = 0 xn2 = copy(fit) xn1 = copy(fit) xn = copy(fit) xnp = zeros( length(fit) ) xy = copy(fit) xn2[i] = fit[i] - 0.01 xn[i] = fit[i] + 0.01 fn2 = fun(xn2) fn1 = fun(xn1) while abs(xn[i] - xn1[i]) > ϵ && n<200 fn = fun(xn) v0 = xn1[i] - xn2[i] v1 = xn[i] - xn2[i] v2 = xn[i] - xn1[i] s0 = xn1[i]*xn1[i] - xn2[i]*xn2[i] s1 = xn[i] *xn[i] - xn2[i]*xn2[i] s2 = xn[i] *xn[i] - xn1[i]*xn1[i] xnp[i] = 0.5(fn2*s2 - fn1*s1 + fn*s0) / (fn2*v2 - fn1*v1 + fn*v0) xn2[i] = xn1[i] xn1[i] = xn[i] xn[i] = xnp[i] fn2 = fn1 fn1 = fn n += 1 println(" $n $xn $xn1") xy = [xy; xn] end fit[i] = xnp[i] plot!(xy[:,1], xy[:,2], w = 3, legend = false, marker = :rect ) end ofDescent(spi, fit = [16.1 16.1], low = [-.1 -.1], up = [.4 .4]) 

Sair de um ponto de partida muito ruim foi feito em três etapas! Bom ... Mas todos os métodos têm uma desvantagem - eles convergem para um mínimo local. Agora vamos assumir a função de Ackley para pesquisa

$$

$$f (x, y) = - 20 \ exp \ left [-0.2 {\ sqrt {0.5 \ left (x ^ {2} + y ^ {2} \ right)}} \ right] - \ exp \ left [ 0,5 \ esquerda (\ cos 2 \ pi x + \ cos 2 \ pi y \ right) \ direita] + e + 20$$

 ekly(x) = -20exp(-0.2sqrt(0.5(x[1]*x[1]+x[2]*x[2]))) - exp(0.5(cospi(2x[1])+cospi(2x[2]))) + 20 + ℯ # f(0,0) = 0, x_i ∈ [-5,5] fun = ekly plotter(surface, low = [-5 -5], up = [5 5]) 

 ofDescent(spi, fit = [1.4 1.4], low = [-.1 -.1], up = [2.4 2.4]) 

Método da proporção áurea

Teoria Embora a implementação seja complicada, o método às vezes se mostra bem pulando os mínimos locais

 function interval(i, fit, st) d = 0. ab = zeros(2) fitc = copy(fit) ab[1] = fitc[i] Fa = fun(fitc) fitc[i] -= st Fdx = fun(fitc) fitc[i] += st if Fdx < Fa st = -st end fitc[i] += st ab[2] = fitc[i] Fb = fun(fitc) while Fb < Fa d = ab[1] ab[1] = ab[2] Fa = Fb fitc[i] += st ab[2] = fitc[i] Fb = fun(fitc) # println("----", Fb, " ", Fa) end if st < 0 c = ab[2] ab[2] = d d = c end ab[1] = d ab end function goldsection(i, fit, ϵ) ϕ = Base.MathConstants.golden ab = interval(i, fit, 0.01) α = ϕ*ab[1] + (1-ϕ)*ab[2] β = ϕ*ab[2] + (1-ϕ)*ab[1] fit[i] = α Fa = fun(fit) fit[i] = β Fb = fun(fit) while abs(ab[2] - ab[1]) > ϕ if Fa < Fb ab[2] = β β = α Fb = Fa α = ϕ*ab[1] + (1-ϕ)*ab[2] fit[i] = α Fa = fun(fit) else ab[1] = α α = β Fa = Fb β = ϕ*ab[2] + (1-ϕ)*ab[1] fit[i] = β Fb = fun(fit) end println(">>>", i, ab) plot!(ab, w = 1, legend = false, marker = :rect ) end fit[i] = 0.5(ab[1] + ab[2]) end ofDescent(goldsection, fit = [1.4 1.4], low = [-.1 -.1], up = [1. 1.]) 

Isso é tudo com descida coordenada. Os algoritmos dos métodos apresentados são bastante simples, portanto, implementá-los no seu idioma preferido não é difícil. No futuro, você poderá considerar as ferramentas integradas da linguagem Julia , mas, por enquanto, deseja tocar em tudo, por assim dizer, com suas mãos, considerar os métodos mais complicados e mais eficientes, para ir para a otimização global, comparando simultaneamente com a implementação em outro idioma.

Literatura

1. Zaitsev V.V. Métodos numéricos para físicos. Equações não lineares e otimização: um tutorial. - Samara, 2005 - 86s.
2. Ivanov A.V. Métodos de computador para otimizar sistemas ópticos. Guia de estudo. –SPb: SPbSU ITMO, 2010 - 114s.
3. Popova T. M. Métodos de otimização multidimensional: diretrizes e tarefas para trabalhos de laboratório na disciplina "Métodos de otimização" para estudantes da direção "Matemática Aplicada" / comp. T. M. Popova. - Khabarovsk: Editora do Pacífico. estado Universidade, 2012 - 44 p.