Autor Tema: Dogleg Trust Region Method Fail

0 Usuarios y 1 Visitante están viendo este tema.

14 Octubre, 2024, 11:08 pm
Leído 507 veces

Suarezemanuelema

  • $$\Large \color{#6a84c0}\pi$$
  • Mensajes: 1
  • País: ar
  • Karma: +0/-0
I'm implementing in Julia the Dogleg Trust Region Method.
I've got this:

function trust_dogleg(f::Function,xₖ,δ,δₖ,η;ϵ=10e-8,iter_max=1000000)

   iter=0
   gₖ=gradiente(f,xₖ)
   Bₖ=hessian(f,xₖ)
   while norm(gₖ)>ϵ && iter<iter_max
      m(t)=f(xₖ)+gₖ'*t+0.5*t'*Bₖ*t
      #Resuelvo 4.3 con Metodo Dogleg
      p_B=Bₖ\(-gₖ)
      p_U= -(((gₖ'*gₖ)/(gₖ'*Bₖ*gₖ))*gₖ)
      if norm(p_B) <= δₖ
         d=p_B
      elseif norm(p_U)>=δₖ
         τ=δₖ/norm(p_U)
         d=τ*p_U
      else
         #p=(p_U + (τ-1)*(p_B - p_U))'*(p_U + (τ-1)*(p_B - p_U))
         #a=norm(p_B-p_U)^2
         a=p_B'*p_B-2*p_B'*p_U+p_U'*p_U
         #b=2*(p_B-p_U)'*p_U
         b=2*p_B'*p_U-2*p_U'*p_U
         #c=norm(p_U)^2-δₖ^2
         c=p_U'*p_U-δₖ^2
         τ=resolvente(a,b,c)+1
         d=p_U + (τ - 1)*(p_B - p_U)
      end
      cero=zeros(length(d))
      #
      ρ=(f(xₖ)-f(xₖ+d))/(m(cero)-m(d))
      if ρ<1/4
         δₖ=(1/4)*δₖ
      else
         if ρ>3/4 && norm(d)==δₖ
            if 2*δₖ<δ
               δₖ=2*δₖ
            else
               δₖ=δ
            end
         #else δₖ=δₖ
         end
      if ρ>η
         xₖ=xₖ+d
         gₖ=gradiente(f,xₖ)
         Bₖ=hessian(f,xₖ)
      #else xₖ=xₖ
      end
      end
      iter=iter+1
   end
   return xₖ,iter

When i use this with this function:

function valle(v::Vector{Float64})

      x, y = v[1], v[2]
       return (x^2 + y^2) * exp(-((x^2 + y^2) / 20))
   end

And call:

iniciodd_2=[2.0,1.0]

trust_dogleg(valle,iniciodd_2,1.5,1,0.2)

I noticed that the method is growing. And the key is that m(cero)-m(d) is being negative (when must be positive, because d is the minimum in a region where there is the 0). On the first step, is calculating d with this block:

elseif norm(p_U)>=δₖ
         τ=δₖ/norm(p_U)
         d=τ*p_U


Also, i'm calculating the gradient with:
function gradiente(fun,x::Vector{Float64})

    nv = length(x)
    h = 0.00001

    xh1 = x.+h
    xh0 = x.-h
    hh  = xh1 - xh0  # 2h
   
    fderiv = zeros(nv)

    if nv == 1
        f1 = fun(xh1)
        f0 = fun(xh0)
        fderiv = (f1-f0)/(2*h)
    else
        for j=1:nv
            xx = copy(x)
            xx[j] = xh1[j]
            f1 = fun(xx)
            xx[j] = xh0[j]
            f0 = fun(xx)
            fderiv[j] = (f1-f0)/(2*h)
        end
    end
   
    return fderiv
end

And the hessian with the Julia function

Any idea of where is failing the algorithm?
Thank you!