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!