# -*- coding: utf-8 -*-
# author: 'boliang'
# date: 2018/5/1 20:29
import numpy as np
'''
最速下降法
'''
class Grad(object):
def __init__(self, fun, gfun, x0):
self.__fun = fun
self.__gfun = gfun
self.__x0 = x0
def grad(self):
epsilon = 1e-5
xk = self.__x0
step = 5000
sigma = 0.4
beta = 0.5
while step > 0:
gk = self.__gfun(xk)
if np.linalg.norm(gk) <= epsilon:
self.__step = 5000-step
break
dk = -gk.reshape(len(xk))
m = 0
mk = 0
while m < 20:
tmp_a = self.__fun(xk + (beta ** m) * dk)
tmp_b = self.__fun(xk) + (sigma * (beta ** m) * self.__gfun(xk).T).dot(dk.T)
if tmp_a <= tmp_b:
mk = m
break
m += 1
xk = xk + beta ** mk*dk
step -= 1
self.__val = self.__fun(xk)
self.__x = xk
def show(self):
print('迭代次数: {0}'.format(self.__step))
for ind, x in enumerate(self.__x):
print('x{:}: {:.3f}'.format(ind+1, x))
print('目标函数值: {0}'.format(self.__val))
'''
牛顿法
'''
class Newton(object):
__max_step = 5000
def __init__(self, fun, gfun, hesse, x0):
self.__x0 = x0
self.__fun = fun
self.__gfun = gfun
self.__hesse = hesse
'''
基本牛顿法
'''
def common_solve(self):
epsilon = 1e-5
xk = self.__x0
step = self.__max_step
while step > 0:
gk = self.__gfun(xk)
if np.linalg.norm(gk) <= epsilon:
self.__step = self.__max_step-step
break
dk = np.linalg.solve(self.__hesse(xk), -gk).reshape(len(xk))
xk = xk + dk
step -= 1
self.__val = self.__fun(xk)
self.__x = xk
'''
阻尼牛顿法
'''
def dampnm(self):
epsilon = 1e-5
beta = 0.6
sigma = 0.4
xk = self.__x0
step = self.__max_step
while step > 0:
gk = self.__gfun(xk)
if np.linalg.norm(gk) <= epsilon:
break
dk = np.linalg.solve(self.__hesse(xk), -gk).reshape(len(xk))
m = 0
mk = 0
while m < 20:
tmp_a = self.__fun(xk + (beta ** m) * dk)
tmp_b = self.__fun(xk) + (sigma * (beta ** m) * self.__gfun(xk).T).dot(dk.T)
if tmp_a <= tmp_b:
mk = m
break
m += 1
xk = xk + beta ** mk * dk
step -= 1
self.__step = self.__max_step - step
self.__val = self.__fun(xk)
self.__x = xk
def show(self):
print('迭代次数: {0}'.format(self.__step))
for ind, x in enumerate(self.__x):
print('x{:}: {:.3f}'.format(ind+1, x))
print('目标函数值: {0}'.format(self.__val))
class Frcg(object):
__max_step = 5000
def __init__(self, fun, gfun, x0):
self.__fun = fun
self.__gfun = gfun
self.__x0 = x0
def solve(self):
epsilon = 1e-4
sigma = 0.4
beta = 0.5
step = self.__max_step
xk = self.__x0
while step > 0:
gk = self.__gfun(xk)
if np.linalg.norm(gk) <= epsilon:
break
if step == self.__max_step:
dk = -gk.reshape(len(xk))
else:
betak = gk.T.dot(gk) / g0.T.dot(g0)
dk = -gk.reshape(len(xk)) + betak*d0
dk = dk.reshape(len(xk))
gd = gk.T.dot(dk)
if gd >= 0.0:
dk = -gk.reshape(len(xk))
m = 0
mk = 0
while m < 20:
tmp_a = self.__fun(xk + (beta ** m) * dk)
tmp_b = self.__fun(xk) + (sigma * (beta ** m) * self.__gfun(xk).T).dot(dk.T)
if tmp_a <= tmp_b:
mk = m
break
m += 1
xk = xk + beta ** mk * dk
g0 = gk
d0 = dk
step -= 1
self.__step = self.__max_step - step
self.__val = self.__fun(xk)
self.__x = xk
def show(self):
print('迭代次数: {0}'.format(self.__step))
for ind, x in enumerate(self.__x):
print('x{:}: {:.3f}'.format(ind+1, x))
print('目标函数值: {0}'.format(self.__val))
def fun(x):
return 100*(x[0]**2-x[1])**2 + (x[0] - 1)**2
def gfun(x):
a = 400*x[0]*(x[0]**2-x[1]) + 2*(x[0]-1)
b = -200*(x[0]**2 - x[1])
return np.array([[a], [b]])
def hesse(x):
a = 1200*x[0]**2 - 400*x[1] + 2
b = -400*x[0]
c = b
d = 200
return np.array([[a, b], [c, d]])
def problem1():
def fun(x):
return 3*x[0]**2 + 2*x[1]**2 - 4*x[0] - 6*x[1]
def gfun(x):
a = 6*x[0]-4
b = 4*x[1]-6
return np.array([[a], [b]])
sol = Grad(fun, gfun, np.array([0, 0]))
sol.grad()
sol.show()
def problem2():
def fun(x):
return 4*x[0]**2 + x[1]**2 - 8*x[0] - 4*x[1]
def gfun(x):
a = 8*x[0]-8
b = 2*x[1]-4
return np.array([[a], [b]])
def hesse(x):
a = 8
b = 0
c = 0
d = 2
return np.array([[a, b], [c, d]])
print('牛顿法............................')
sol = Newton(fun, gfun, hesse, np.array([0, 0]))
sol.common_solve()
sol.show()
print('\n阻尼牛顿法.......................')
sol.dampnm()
sol.show()
def problem3():
def fun(x):
return (x[0]-2)**4 + (x[0]-2*x[1])**2
def gfun(x):
a = 4*(x[0]-2)**3 + 2*(x[0]-2*x[1])
b = -4*(x[0]-2*x[1])
return np.array([[a], [b]])
sol = Grad(fun, gfun, np.array([0, 3]))
sol.grad()
sol.show()
def problem4():
def fun(x):
return 100*(x[1]-x[0]**2)**2 + (x[0]-1)**2
def gfun(x):
a = -400*x[0]*(x[1]-x[0]**2) + 2*(x[0]-1)
b = 200*(x[1]-x[0]**2)
return np.array([[a], [b]])
def hesse(x):
a = -400*(x[1]-x[0]**2) + 800*x[0]**2 + 2
b = -400*x[0]
c = -400*x[0]
d = 200
return np.array([[a, b], [c, d]])
sol = Newton(fun, gfun, hesse, np.array([-1.2, 1]))
sol.common_solve()
sol.show()
def problem5():
def fun(x):
return 4*x[0]**2 + 4*x[1]**2 - 4*x[0]*x[1] - 12*x[1]
def gfun(x):
a = 8 * x[0] - 4 * x[1]
b = 8 * x[1] - 4 * x[0] - 12
return np.array([[a], [b]])
sol = Frcg(fun, gfun, np.array([-0.5, 1]))
sol.solve()
sol.show()
def problem6():
def fun(x):
return x[0]**2 - 2*x[0]*x[1] + 2*x[1]**2 + x[2]**2 - x[1]*x[2] + x[0] + 3*x[1] - x[2]
def gfun(x):
a = 2*x[0] - 2*x[1] + 1
b = -2*x[0] + 4*x[1] - x[2] + 3
c = 2*x[2] - x[1] - 1
return np.array([[a], [b], [c]])
sol = Frcg(fun, gfun, np.array([0, 0, 0]))
sol.solve()
sol.show()
if __name__ == '__main__':
# sol = Grad(fun, gfun, np.array([-1.2, 1]))
# sol.grad()
# sol.show()
#
# sol = Newton(fun, gfun, hesse, np.array([0, 0]))
# sol.dampnm()
# sol.show()
#
# sol = Frcg(fun, gfun, np.array([-1.2, 1]))
# sol.solve()
# sol.show()
problem1()
print()
problem2()
print()
problem3()
print()
problem4()
print()
problem5()
print()
problem6()
测试结果如下
Problem1 test 最速下降法
迭代次数: 10
x1: 0.667
x2: 1.500
目标函数值: -5.833333333328483
Problem2 test 牛顿法............................
迭代次数: 1
x1: 1.000
x2: 2.000
目标函数值: -8.0
Problem2 test 阻尼牛顿法.......................
迭代次数: 1
x1: 1.000
x2: 2.000
目标函数值: -8.0
Problem3 test 最速下降法
迭代次数: 2111
x1: 2.014
x2: 1.007
目标函数值: 3.7684525980724333e-08
Problem4 test 牛顿法
迭代次数: 5
x1: 1.000
x2: 1.000
目标函数值: 1.8527397132921126e-11
Problem5 test 共轭梯度下降法
迭代次数: 9
x1: 1.000
x2: 2.000
目标函数值: -11.999999999642554
Problem6 test 共轭梯度下降法
迭代次数: 23
x1: -2.833
x2: -2.333
x3: -0.667
目标函数值: -4.583333332921872
Comments