In [1]:
function newtmin( obj, x0; maxIts=100, optTol=1e-6)
# Newton method with Armijo backtracking and some treatment with Hessian, treatment 2
# Minimize a function f using Newton’s method.
# obj: a function that evaluates the objective value,
# gradient, and Hessian at a point x, i.e.,
# (f, g, H) = obj(x)
# x0: starting point.
# maxIts (optional): maximum number of iterations.
# optTol (optional): optimality tolerance based on
# ||grad(x)|| <= optTol*||grad(x0)||
mu=1e-4
epsilon=0.01
x=x0
status = 0
its = 0
(f0,g0,H0)=obj(x0)
(f,g,H)=obj(x)
(V,S)=eig(H)
bar=ones(size(H,1))*epsilon
bar=Diagonal(bar)
bH=S*max(bar,abs(Diagonal(V)))*S'
while status != 1
alpha = 1
xnew = x-alpha*inv(bH)*g
(fnew,gnew,Hnew)=obj(xnew)
sts =-fnew+f-alpha*mu*g'*inv(bH)*g
while sts[1]<0
alpha=alpha/2
xnew = x-alpha*inv(bH)*g
(fnew,gnew,Hnew)=obj(xnew)
sts=-fnew+f-alpha*mu*g'*inv(bH)*g
end
x = x-alpha*inv(bH)*g
(f,g,H)=obj(x)
(V,S)=eig(H)
bar=ones(size(H,1))*epsilon
bar=Diagonal(bar)
bH=S*max(bar,abs(Diagonal(V)))*S'
its = its+1
if norm(g)<= optTol*norm(g0)
status = 1
end
if its>maxIts
status = 1
end
end
return (x, its)
end
Out[1]:
In [2]:
using ForwardDiff
using ReverseDiffSource
In [16]:
x0=[-1.2;1]
f(x)=(1-x[1])^2
c(x)=10*(x[2]-x[1]^2)
g(x)=ForwardDiff.gradient(x -> f(x))(x)
H(x)=ForwardDiff.hessian(x -> f(x))(x)
J(x)=(ForwardDiff.gradient(x -> c(x))(x))'
Out[16]:
In [17]:
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=ForwardDiff.hessian(z -> d1(z))(z)
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=ForwardDiff.hessian(z -> d2(z))(z)
result=H(x)+dh1(x)+dh2(x)
return result
end
Out[17]:
In [22]:
function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou) # objective value at x
lagg = nxLrou(x,y,rou) # gradient at x
lagH = nxxLrou(x,y,rou) # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
y=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
sts=1
end
x=xnew
end
return x
end
Out[22]:
In [27]:
auglag( x0; maxIts=100, optTol=1e-6)
Out[27]:
In [53]:
x0=[2;2]
f(x)=log(1+x[1]^2)-x[2]
c(x)=(1+x[1]^2)^2+x[2]^2-4
J(x)=(ForwardDiff.gradient(x -> c(x))(x))'
g(x)=[2*x[1]/(1+x[1]^2);-1]
H(x)=[(2-2*x[1]^2)/(1+x[1]^2)^2 0;0 0]
Out[53]:
In [54]:
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(x)=-y*[4+12*x[1]^2 0;0 2]
d2(z)=rou/2*c(z)'*c(z)
dh2(x)=rou*[4+12*x[1]^2 0;0 2]*c(x)+rou*J(x)'*J(x)
result=H(x)+dh1(x)+dh2(x)
return result
end
Out[54]:
In [55]:
function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou) # objective value at x
lagg = nxLrou(x,y,rou) # gradient at x
lagH = nxxLrou(x,y,rou) # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
ynew=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
if abs(abs(J(xnew)*g(xnew)/(norm(J(xnew))*norm(g(xnew))))-1)[1]<optTol
sts=1
end
end
x=xnew
y=ynew
end
return x
end
auglag( x0; maxIts=100, optTol=1e-6)
Out[55]:
In [6]:
x0=[2;1]
f(x)=-1
c(x)=[x[1]^2+x[2]^2-25; x[1]*x[2]-9]
g(x)=[0;0]
H(x)=[0 0 ;0 0 ]
J(x)=([2*x[1] x[2];2*x[2] x[1]])'
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=-y[1]*[2 0;0 2]-y[2]*[0 1; 1 0]
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=rou*[2 0;0 2]*c(x)[1]+rou*[0 1; 1 0]*c(x)[2]+rou*J(x)'*J(x)
result=H(x)+dh1(x)+dh2(x)
return result
end
Out[6]:
In [54]:
function auglag( x0; maxIts=100, optTol=1e-6)
y0=[1;1]
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou) # objective value at x
lagg = nxLrou(x,y,rou) # gradient at x
lagH = nxxLrou(x,y,rou) # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
ynew=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
sts=1
end
x=xnew
y=ynew
end
return x
end
auglag( x0; maxIts=100, optTol=1e-6)
Out[54]:
In [119]:
x0=[0;0]
f(x)=sin(pi*x[1]/12)*cos(pi*x[2]/16)
c(x)=4*x[1]-3*x[2]
g(x)=[pi/12*cos(pi*x[1]/12)*cos(pi*x[2]/16); -pi/16*sin(pi*x[1]/12)*sin(pi*x[2]/16)]
H(x)=[-pi^2/(12^2)*sin(pi*x[1]/12)*cos(pi*x[2]/16) -pi^2/(12*16)*cos(pi*x[1]/12)*sin(pi*x[2]/16);-pi^2/(16*12)*cos(pi*x[1]/12)*sin(pi*x[2]/16) -pi^2/(16^2)*sin(pi*x[1]/12)*cos(pi*x[2]/16)]
J(x)=(ForwardDiff.gradient(x -> c(x))(x))'
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=ForwardDiff.hessian(z -> d1(z))(z)
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=ForwardDiff.hessian(z -> d2(z))(z)
result=H(x)+dh1(x)+dh2(x)
return result
end
Out[119]:
In [120]:
function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou) # objective value at x
lagg = nxLrou(x,y,rou) # gradient at x
lagH = nxxLrou(x,y,rou) # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
y=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
if abs(abs(J(xnew)*g(xnew)/(norm(J(xnew))*norm(g(xnew))))-1)[1]<optTol
sts=1
end
end
x=xnew
end
return x
end
Out[120]:
In [121]:
auglag( x0; maxIts=100, optTol=1e-6)
Out[121]: