import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
mpl.rcParams['font.sans-serif'] = 'SimHei'
mpl.rcParams['axes.unicode_minus'] = False
def num_grad(x, h):
df = np.zeros(x.size)
for i in range(x.size):
x1, x2 = x.copy(), x.copy()
x1[i] = x[i] - h
x2[i] = x[i] + h
y1, y2 = f(x1), f(x2)
df[i] = (y2 - y1) / (2 * h)
return df
def num_hess(x, h):
hess = np.zeros((x.size, x.size))
for i in range(x.size):
x1 = x.copy()
x1[i] = x[i] - h
df1 = num_grad(x1, h)
x2 = x.copy()
x2[i] = x[i] + h
df2 = num_grad(x2, h)
d2f = (df2 - df1) / (2 * h)
hess[i] = d2f
return hess
def linesearch(x, pk:float,f):
'''
黄金分割法
'''
eps = 0.001
r = 500
l = -500
while r - l > eps:
t1 = l + 0.382 * (r - l)
t2 = l + 0.618 * (r - l)
x_1 = x+ t1 * pk
x_2 = x + t2 * pk
f1 = f(x_1)
f2 = f(x_2)
if f1 < f2:
r = t2
else:
l = t1
if t2 - t1 < eps:
return t2
def steepest(x,f,epsilon=0.01,h=10**-5,maxiter=10**4):
'''
最速下降法
:param x: 初始值
:param f: 函数
:param epsilon: eps
:param h: 求导的delta x
:param maxiter: 最大迭代次数
:return: 答案的list用于画图
'''
x_list=[]
ans_list=[]
x_list.append(x)
ans_list.append(f(x))
for iter1 in range(maxiter):
grad = num_grad(x, h)
pk = -grad
ak = linesearch(x, pk,f)
x = x + ak * pk
print('grad:', grad, 'x_{i+1}:', x, 't:', ak,'ans:',f(x))
x_list.append(x)
ans_list.append(f(x))
if np.linalg.norm(grad) < epsilon:
return x_list,ans_list
return x_list,ans_list
def newTonFuction(x,f,epsilon=0.01,h1=10**-5,h2=10**-5,maxiter=10**4):
x_list=[]
ans_list=[]
x_list.append(x)
ans_list.append(f(x))
for iter1 in range(maxiter):
grad = num_grad(x, h1)
hess = num_hess(x, h2)
pk = -np.dot((np.linalg.inv(hess)), grad)
x = x + pk
x_list.append(x)
ans_list.append(f(x))
print('grad:', grad, 'x_{i+1}:', x, 'ans:', f(x))
if np.linalg.norm(grad) < epsilon:
return x_list,ans_list
return x_list,ans_list
def BFGS(x,f,epsilon=0.01,h=10**-5,maxiter=10**4):
'''
拟牛顿法
:param x: 初始点
:param f: 求解函数
:param epsilon: eps
:param h: 求导的delta x
:param maxiter: 最大迭代数
:return: 答案的list便于画图
'''
Bk = np.eye(x.size)
x_list=[]
ans_list=[]
x_list.append(x)
ans_list.append(f(x))
for i in range(maxiter):
grad = num_grad(x, h)
pk = -np.dot((np.linalg.inv(Bk)), grad)
ak = linesearch(x, pk,f)
x = x + pk * ak
yk = num_grad(x, h) - grad
sk = ak * pk
if np.dot(yk.reshape(1, grad.shape[0]), sk) > 0:
Bk = Bk - np.dot(np.dot(np.dot(Bk, sk).reshape(sk.shape[0], 1), sk.reshape(1, sk.shape[0])), Bk)\
/ np.dot(np.dot(sk.reshape(1, sk.shape[0]), Bk), sk) + \
np.dot(yk.reshape(yk.shape[0], 1),yk.reshape(1, yk.shape[0])) / \
np.dot(yk.reshape(1, yk.shape[0]), sk)
print('grad{%d}:' % (i + 1), np.round(num_grad(x, h), 3),
'p{%d}:' % (i + 1), np.round(pk, 3),
'x_{%d}:' % (i + 1), np.round(x, 3), f(x))
x_list.append(x)
ans_list.append(f(x))
if np.linalg.norm(grad) < epsilon:
return x_list,ans_list
return x_list,ans_list
def modify_newton(x,f,epsilon=0.01,h1=10**-5,h2=10**-5,maxiter=10**4):
ans_list=[]
x_list=[]
x_list.append(x)
ans_list.append(f(x))
for i in range(maxiter):
grad=num_grad(x,h1)
hess = num_hess(x, h2)
pk = -np.dot((np.linalg.inv(hess)), grad)
t=linesearch(x,pk,f)
x = x + pk*t
print('grad:',np.round(grad,3),'pk:',np.round(pk,3),'t:',np.round(t,3),'x_{%d}:'%i,np.round(x,3),f(x))
x_list.append(x)
ans_list.append(f(x))
if np.linalg.norm(grad)<epsilon:
return x_list,ans_list
return x_list,ans_list
def conjugate_gradient(x,f,epsilon=0.01,h1=10**-5,h2=10**-5,maxiter=10**4):
ans_list = []
x_list = []
x_list.append(x)
ans_list.append(f(x))
pk=-num_grad(x,h1)
if np.linalg.norm(pk)<epsilon:
return x
for i in range(maxiter):
grad=num_grad(x,h1)
t=np.round(linesearch(x,pk,f),2)
used_x=x.copy()
x = x + pk*t
pk=-num_grad(x,h1)+(np.linalg.norm(num_grad(x,h1))**2/np.linalg.norm(num_grad(used_x,h1))**2)*pk
print('grad{%d}:'%(i+1),np.round(num_grad(x,h1),3),'p{%d}:'%(i+1),np.round(pk,3),'t_{%d}:'%(i+1),np.round(t,3),'x_{%d}:'%(i+1),np.round(x,3),f(x))
x_list.append(x)
ans_list.append(f(x))
if np.linalg.norm(grad)<epsilon:
return x_list,ans_list
return x_list,ans_list
def simplex_method(x:list,f,maxiter:int=100,alpha:float=1.,beta:float=1.,epsilon:float=0.01):
temp_x_list=x
x_list=[]
ans_list=[]
for i in range(maxiter):
l=[]
for j in temp_x_list:
l.append(f(j))
l=np.array(l)
x_min = temp_x_list[l.argsort()[0]]
x_mid = temp_x_list[l.argsort()[1]]
x_max = temp_x_list[l.argsort()[2]]
loss=sum([(f(j)-f(x_min))**2 for j in x])
tempx=(np.sum(x,axis=0)-x_max)/ (len(x)-1)
reflect_x=2*tempx-x_max
print('reflect_x:',reflect_x)
if f(x_min)>f(reflect_x):
inflation_x=tempx+alpha*(reflect_x-tempx)
if f(inflation_x)<f(reflect_x):
temp_x_list[l.argsort()[2]]=inflation_x
else:
temp_x_list[l.argsort()[2]]=reflect_x
print('f(x_min):%.3f>f(reflect):%.3f'%(f(x_min),f(reflect_x)))
print('inflation_x:',inflation_x)
if f(x_min)<=f(reflect_x) and f(reflect_x)<f(x_mid):
temp_x_list[l.argsort()[2]]=reflect_x
print('f(x_min):%.3f<=f(reflect):%.3f<=f(x_mid):%.3f' % (f(x_min), f(reflect_x),f(x_mid)))
if f(x_mid)<=f(reflect_x) and f(reflect_x)<f(x_max):
shrink_x=tempx+beta*(reflect_x-tempx)
temp_x_list[l.argsort()[2]]=shrink_x
print('f(x_mid):%.3f<=f (reflect):%.3f <f(x_max):%.3f' % (f(x_mid),f(reflect_x), f(x_max)))
print('shrink_x:',shrink_x)
if f(reflect_x)>=f(x_max):
shrink_x=tempx+beta*(x_max-tempx)
if f(shrink_x)>f(x_max):
for j in x:
j=(x_min+j)/2
temp_x_list=x
else:
temp_x_list[l.argsort()[2]]=shrink_x
print('f(reflect):%.3f>f(x_max):%.3f' % (f(reflect_x),f(reflect_x)))
print('shrink_x:',shrink_x)
print('temp_x_list:',temp_x_list)
print('loss:%.3f' % loss)
if loss < epsilon:
break
x_list.append(x_min)
ans_list.append(f(x_min))
print()
return x_list,ans_list
def cyclic_coordinate_method(x:list,f,epsilon:float=0.01,h1:float=1e-5,maxiter:int=100):
ans_list = []
x_list = []
x_list.append(x)
ans_list.append(f(x))
I=np.eye(2)
for i in range(maxiter):
used_x=x.copy()
for pk in I:
t = np.round(linesearch(x, pk,f),2)
x = x + pk * t
print('p{%d}:' % (i + 1), np.round(pk, 3),
't_{%d}:' % (i + 1), np.round(t, 2), 'x_{%d}:' % (i + 1), np.round(x, 3), f(x),'eps:',np.linalg.norm(used_x-x))
x_list.append(x)
ans_list.append(f(x))
if np.linalg.norm(used_x-x)< epsilon:
return x_list,ans_list
return x_list,ans_list
def draw(x:list=None,y:list=None,ans:list=None):
fig=plt.figure()
ax = fig.gca(projection='3d')
x1 = np.linspace(min(x),max(x), 100)
x2 = np.linspace(min(x), max(x), 100)
tempx, tempy = np.meshgrid(x1, x2)
surf = ax.plot_surface(tempx, tempy, f([tempx,tempy]), cmap=cm.jet, zorder=10)
scat=ax.scatter(x,y,ans,zorder=1000,marker='x')
line=ax.plot(x,y,ans,zorder=100,linestyle='--')
contour=ax.contour(tempx,tempy,f([tempx,tempy]),50,zorder=100000)
plt.legend()
plt.show()
return
def f(x):
x1 = x[0]
x2 = x[1]
y=x1**2+x2**2+(x2-1)**2
return y
if __name__ == '__main__':
x0 = np.array([1., 1.])
x_list,ans_list=BFGS(x0,f,maxiter=100)
draw([i[0] for i in x_list],[i[1] for i in x_list],ans_list)
|
