SD-Newton-Enhanced SD
Comparation Steepest Decent Method, Newton Method and Enhanced Steepest Decent Method¶
Enhanced Steepest Decent
By Mohamed Kamel Riahi on "A new approach to improve ill-conditioned parabolic optimal control problem via time domain decomposition"
Numer Algor (2016) 72:635-666
DOI 10.1007/s11075-015-0060-0
https://www.researchgate.net/publication/278829050_A_new_approach_to_improve_ill-conditioned_parabolic_optimal_control_problem_via_time_domain_decomposition
Example 1¶
Consider a problem minimize : $$f(\bar{x})=x_1^2+2x_2^2+3x_3^2+x_4^2+2x_5^2+4x_6^2-x_1-x_2-x_3-x_4-x_5-x_6$$
The problem will be solved using:
- Steepest Decent Method
- Newton Method
- Enhanced Steepest Decent:
- $\hat{n}=1$, equivalent to Steepest Decent Method
- $\hat{n}=2$
- $\hat{n}=3$
- $\hat{n}=6$, equivalent to Newton Method
In [1]:
import autograd.numpy as np
from autograd import grad, jacobian, hessian
def f(x):
return x[0]**2 +2*x[1]**2+3*x[2]**2+1*x[3]**2+2*x[4]**2+4*x[5]**2-x[0]-x[1]-x[2]-x[3]-x[4]-x[5]
In [2]:
def nablaf(x,f):
j=jacobian(f)
return j(x)
def nabla2f(x,f):
j=hessian(f)
return j(x)
def Anablaf(x,f):
return nabla2f(x,f)@nablaf(x,f)
def theta(x,f):
return -(nablaf(x,f)@nablaf(x,f))/(nablaf(x,f)@Anablaf(x,f))
def steepest_decent(x,f):
stop=False
iter=0; max_iter=250
error=0.1**8
print("Initial x= ",x)
while stop==False:
iter+=1
r=nablaf(x,f)
Ar=Anablaf(x,f)
t=theta(x,f)
x=x+t*r
if iter%10==1:
print('Iter =',iter) ; print('x =',x)
if iter>=max_iter or np.linalg.norm(nablaf(x,f))<=error: stop = True
if iter!=1: print('Iter =',iter) ; print('x =',x)
return x,iter
def newton(x,f):
stop=False
iter=0; max_iter=250
error=0.1**8
print("Initial x= ",x)
while stop==False:
iter+=1
r=nablaf(x,f)
Ar=Anablaf(x,f)
t=theta(x,f)
x=x-np.linalg.inv(nabla2f(x,f))@r
print('Iter =',iter); print('x =',x)
if iter>=max_iter or np.linalg.norm(nablaf(x,f))<=error: stop = True
return x,iter
In [3]:
print('\n=======Newton Method========')
x=np.array([1,1,1,1,1,1], dtype=float)
x,iter=newton(x,f)
print('Total iteration number',iter); print('x =',x)
=======Newton Method======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ] Total iteration number 1 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ]
In [4]:
print('\n=======Steepest Decent Method========')
x=np.array([1,1,1,1,1,1], dtype=float)
x,iter=steepest_decent(x,f)
print('Total iteration number',iter); print('x =',x)
=======Steepest Decent Method======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [ 0.84789644 0.54368932 0.2394822 0.84789644 0.54368932 -0.06472492] Iter = 11 x = [0.50222675 0.25000003 0.16666667 0.50222675 0.25000003 0.12430865] Iter = 21 x = [0.50001276 0.25 0.16666667 0.50001276 0.25 0.12499604] Iter = 31 x = [0.50000007 0.25 0.16666667 0.50000007 0.25 0.12499998] Iter = 38 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ] Total iteration number 38 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ]
In [5]:
def pin(x,n,n_hat):
xtemp= np.zeros(len(x))
m=len(x)
a=int((n-1)*(m/n_hat)+1)
b=int(n*(m/n_hat))
for i in range(a,b+1):
xtemp[i-1]=x[i-1]
return xtemp
def phi_nhat(r):
phi=np.zeros(len(r))
for i in range(len(r)):
phi[i]=r[i]@r[i]
#print('phi',phi)
return phi
def phi2_nhat(r,x,f):
phi2=np.zeros((len(r),len(r)))
nn=nabla2f(x,f)
#print('nn',nn)
for i in range(len(r)):
for j in range(len(r)):
phi2[i][j]=r[i]@(nn@r[j])
return phi2
def Theta_nhat(r,x,f):
try:
return -np.linalg.inv(phi2_nhat(r,x,f))@phi_nhat(r)
except:
print('======================The Phi^2 matrix is singular===========================')
print('Current x=',x)
print('Current Phi^2=',phi2_nhat(r,x,f))
def get_rn(x,f,n_hat):
x_tilda=[]
rn=[]
for n in range(n_hat):
x_tilda.append(pin(x,n+1,n_hat))
nn=nablaf(x,f)
rn.append(pin(nn,n+1,n_hat))
return rn
def update_x(x,rn,n_hat):
s=np.zeros(len(x))
for i in range(n_hat):
s+=Theta_nhat(rn,x,f)[i]*rn[i]
x=x+s
return x
def get_r(rn):
r= 0
for i in range(len(rn)):
r+=np.linalg.norm(rn[i])
#print(r)
return r
def enhanced_sd(x,f,n_hat):
stop=False
iter=0; max_iter=250
error=0.1**3
print("Initial x= ",x)
while stop==False:
iter+=1
rn=get_rn(x,f,n_hat)
try:
x=update_x(x,rn,n_hat)
except:
break
#if iter%10==1:
print('Iter =',iter) ; print('x =',x)
if iter>=max_iter or get_r(get_rn(x,f,n_hat))<=error: stop = True
return x,iter
In [6]:
print('\n=======Enhanced Steepest Decent Method with n_hat=1 ========')
x=np.array([1,1,1,1,1,1], dtype=float) #Initial Value
n_hat=1
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=1 ======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [ 0.84789644 0.54368932 0.2394822 0.84789644 0.54368932 -0.06472492] Iter = 2 x = [0.71488162 0.31911081 0.15596127 0.71488162 0.31911081 0.225433 ] Iter = 3 x = [0.6378638 0.26956945 0.16747235 0.6378638 0.26956945 0.08144431] Iter = 4 x = [0.57838254 0.25268299 0.16642952 0.57838254 0.25268299 0.15661262] Iter = 5 x = [0.54929779 0.25069188 0.16669351 0.54929779 0.25069188 0.1096917 ] Iter = 6 x = [0.52792256 0.25009189 0.16665859 0.52792256 0.25009189 0.13624205] Iter = 7 x = [0.51755714 0.25002367 0.16666758 0.51755714 0.25002367 0.11954894] Iter = 8 x = [0.5099441 0.25000314 0.16666639 0.5099441 0.25000314 0.12900358] Iter = 9 x = [0.50625263 0.25000081 0.1666667 0.50625263 0.25000081 0.12305871] Iter = 10 x = [0.5035414 0.25000011 0.16666666 0.5035414 0.25000011 0.1264258 ] Iter = 11 x = [0.50222675 0.25000003 0.16666667 0.50222675 0.25000003 0.12430865] Iter = 12 x = [0.5012612 0.25 0.16666667 0.5012612 0.25 0.12550777] Iter = 13 x = [0.50079301 0.25 0.16666667 0.50079301 0.25 0.12475379] Iter = 14 x = [0.50044915 0.25 0.16666667 0.50044915 0.25 0.12518083] Iter = 15 x = [0.50028242 0.25 0.16666667 0.50028242 0.25 0.12491232] Iter = 16 x = [0.50015996 0.25 0.16666667 0.50015996 0.25 0.1250644 ] Total iteration number 16 x = [0.50015996 0.25 0.16666667 0.50015996 0.25 0.1250644 ]
In [7]:
print('\n=======Enhanced Steepest Decent Method with n_hat=2 ========')
x=np.array([1,1,1,1,1,1], dtype=float) #Initial Value
n_hat=2
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=2 ======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [0.81382979 0.44148936 0.06914894 0.8627907 0.58837209 0.03953488] Iter = 2 x = [0.65379357 0.24619093 0.21833524 0.69362202 0.27280745 0.19894383] Iter = 3 x = [0.57719708 0.24998512 0.14113505 0.63079932 0.25800722 0.10297642] Iter = 4 x = [0.53844681 0.25000006 0.17958303 0.56326655 0.24973883 0.14846024] Iter = 5 x = [0.51929838 0.25 0.16028405 0.54260807 0.24990939 0.11781827] Iter = 6 x = [0.50961126 0.25 0.16989561 0.5206171 0.25000292 0.13264487] Iter = 7 x = [0.50482437 0.25 0.16507109 0.51388493 0.25000101 0.12265965] Iter = 8 x = [0.5024027 0.25 0.16747387 0.50671861 0.24999997 0.12749128] Iter = 9 x = [0.50120604 0.25 0.16626779 0.50452476 0.24999999 0.12423734] Iter = 10 x = [0.50060065 0.25 0.16686846 0.50218943 0.25 0.12581185] Iter = 11 x = [0.5003015 0.25 0.16656695 0.50147451 0.25 0.12475147] Iter = 12 x = [0.50015016 0.25 0.16671711 0.50071348 0.25 0.12526456] Iter = 13 x = [0.50007537 0.25 0.16664174 0.50048051 0.25 0.12491901] Iter = 14 x = [0.50003754 0.25 0.16667928 0.50023251 0.25 0.12508621] Total iteration number 14 x = [0.50003754 0.25 0.16667928 0.50023251 0.25 0.12508621]
In [8]:
print('\n=======Enhanced Steepest Decent Method with n_hat=3 ========')
x=np.array([1,1,1,1,1,1], dtype=float) #Initial Value
n_hat=3
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=3 ======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [0.73684211 0.21052632 0.14473684 0.82894737 0.59345794 0.05140187] Iter = 2 x = [0.5215311 0.28229665 0.20582707 0.52349624 0.29613614 0.1788255 ] Iter = 3 x = [0.51019894 0.24830018 0.16563613 0.51545805 0.27112777 0.12047262] Iter = 4 x = [0.50092718 0.25139076 0.16850691 0.50110415 0.25283806 0.12831107] Iter = 5 x = [0.50043919 0.2499268 0.16661824 0.50072641 0.25129967 0.1247215 ] Iter = 6 x = [0.50003993 0.25005989 0.16675314 0.50005189 0.25017458 0.12520368] Iter = 7 x = [0.50001891 0.24999685 0.16666439 0.50003414 0.25007995 0.12498287] Total iteration number 7 x = [0.50001891 0.24999685 0.16666439 0.50003414 0.25007995 0.12498287]
In [9]:
print('\n=======Enhanced Steepest Decent Method with n_hat=6 ========')
x=np.array([1,1,1,1,1,1], dtype=float) #Initial Value
n_hat=6
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=6 ======== Initial x= [1. 1. 1. 1. 1. 1.] Iter = 1 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ] Total iteration number 1 x = [0.5 0.25 0.16666667 0.5 0.25 0.125 ]
The result of Enhanced Steepest Decent Method with $\hat{n}=1$ same as Steepest Descent Method and Enhanced Steepest Decent Method with $\hat{n}=6$ same as Newton Method
Change the initial point¶
In [10]:
print('\n=======Enhanced Steepest Decent Method with n_hat=3 ========')
x=np.array([1,1,1,1/2,1/2,1/2], dtype=float) #Initial Value
n_hat=3
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=3 ======== Initial x= [1. 1. 1. 0.5 0.5 0.5] Iter = 1 x = [0.73684211 0.21052632 0.16666667 0.5 0.36842105 0.10526316] Iter = 2 x = [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833] ======================The Phi^2 matrix is singular=========================== Current x= [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833] Current Phi^2= [[0.07046542 0. 0. ] [0. 0. 0. ] [0. 0. 0.14093084]] Total iteration number 3 x = [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833]
The enhanced steepest descent fail to solve the problem, because the $\Phi''_n{\hat{n}}(0)$ is singular
In [11]:
print('\n=======Enhanced Steepest Decent Method with n_hat=6 ========')
x=np.array([1,1,1,1/2,1/2,1/2], dtype=float) #Initial Value
n_hat=6
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=6 ======== Initial x= [1. 1. 1. 0.5 0.5 0.5] ======================The Phi^2 matrix is singular=========================== Current x= [1. 1. 1. 0.5 0.5 0.5] Current Phi^2= [[ 2. 0. 0. 0. 0. 0.] [ 0. 36. 0. 0. 0. 0.] [ 0. 0. 150. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 4. 0.] [ 0. 0. 0. 0. 0. 72.]] Total iteration number 1 x = [1. 1. 1. 0.5 0.5 0.5]
The enhanced steepest descent fail to solve the problem, because the $\Phi''_n{\hat{n}}(0)$ is singular
Example 2¶
Consider a problem minimize : $$f(\bar{x})=x_1^2+2x_2^2+3x_3^2+x_4^2+2x_5^2+4x_6^2-x_1-x_2-x_3-x_4-x_5-x_6$$
The problem will be solved using:
- Enhanced Steepest Decent:
- $\hat{n}=3$
In [13]:
def f(x):
return x[0]**2 +2*x[1]**2+3*x[2]**2+1*x[3]**2+2*x[4]**2+4*x[5]**2-x[0]-x[1]-x[2]-x[3]-x[4]-x[5]
print('\n=======Enhanced Steepest Decent Method with n_hat=6 ========')
x=np.array([1,1,1,1/2,1/2,1/2], dtype=float) #Initial Value
n_hat=3
x,iter=enhanced_sd(x,f,n_hat)
print('Total iteration number',iter); print('x =',x)
=======Enhanced Steepest Decent Method with n_hat=6 ======== Initial x= [1. 1. 1. 0.5 0.5 0.5] Iter = 1 x = [0.73684211 0.21052632 0.16666667 0.5 0.36842105 0.10526316] Iter = 2 x = [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833] ======================The Phi^2 matrix is singular=========================== Current x= [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833] Current Phi^2= [[0.07046542 0. 0. ] [0. 0. 0. ] [0. 0. 0.14093084]] Total iteration number 3 x = [0.5215311 0.28229665 0.16666667 0.5 0.26076555 0.14114833]
The enhanced steepest descent fail to solve the problem, because the $\Phi''_n{\hat{n}}(0)$ is singular
In example 1, the result of
Enhanced Steepest Decent Methodwith $\hat{n}=1$ same asSteepest Descent MethodandEnhanced Steepest Decent Methodwith $\hat{n}=6$ same asNewton Method
However, if the initial point change Enhanced Steepest Decent Method might be fail to solve the problem.
The problem is that we dont know when the matrix $\Phi''_n{\hat{n}}(0)$ become singular.
Apabila ada masukkan silahkan sampaikan
Terima kasih
Comments
Post a Comment