computeGradient 梯度计算方法见文章
理查德外推法计算偏导数近似值-python实现
梯度法
梯度下降法的优化思想是用当前位置负梯度方向作为搜索方向,因为该方向为当前位置的最快下降方向,所以也被称为是”最速下降法“。最速下降法越接近目标值,步长越小,前进越慢
流程
实现
from matplotlib import pyplot
from computeGradient import computeGrad
def fun(x):
return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
def gradient(x):
grad = computeGrad(fun, x, 1e-3)
return grad
def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
"""
非精确线搜索Armijo准则
:param beta:(0,1)
:param sigma:(0,0.5)
:param maxm:最大m值(m 非负整数0,1,2,3,...)
:return:mk
"""
m = 0
mk = 0
while m < maxm:
if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
mk = m
break
m += 1
return mk
def Gradient(fun, gfun, x0, epsilon):
"""
Gradient 梯度法
:param fun:目标函数
:param gfun:梯度
:param x0:初始点
:param epsilon:容许误差
:return:k 迭代次数,x 近似最优点,val 近似最优值
"""
maxIteration = 5000
beta = 0.5
sigma = 0.4
k = 0
graphArr = []
graphArr.append(x0)
while k < maxIteration:
gk = gfun(x0)
dk = -1.0 * gk
if np.linalg.norm(gk) < epsilon:
break
mk = armijo(fun,x0,dk,gk,beta,sigma,maxm=20)
x0 = x0 + beta ** mk * dk
graphArr.append(x0)
k += 1
x = x0
val = fun(x)
return k, x, val, graphArr
k, x, val, graphArr = Gradient(fun, gradient, np.array([-1.2, 1.0]), 1e-5)
print "minPoint:", x, "minValue:", int(val), "iteration:", k
arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()
结果
minPoint: [ 0.99998906 0.99997808] minValue: 0 iteration: 1435
阻尼牛顿法
基础牛顿法初始点需要靠近极小点,否则算法可能不收敛。为了克服初始点选取困难的问题,阻尼牛顿法引入了线搜索方法以得到大范围收敛。
流程
实现
import numpy as np
from matplotlib import pyplot
from computeGradient import computeGrad
from computeHessian import computeHes
def fun(x):
return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
def gradient(x):
grad = computeGrad(fun, x, 1e-3)
return grad
def hessianMatrix(x):
Hes = computeHes(fun, x, 1e-3)
return Hes
def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
"""
非精确线搜索Armijo准则
:param beta:(0,1)
:param sigma:(0,0.5)
:param maxm: 最大m值(m 非负整数0,1,2,3,...)
:return:mk
"""
m = 0
mk = 0
while m < maxm:
if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
mk = m
break
m += 1
return mk
def dampingNewton(fun, gfun, Hess, x0, epsilon):
maxIteration = 5000
beta = 0.5
sigma = 0.4
k = 0
graphArr = []
graphArr.append(x0)
while k < maxIteration:
gk = gfun(x0)
Gk = Hess(x0)
dk = -1.0 * np.linalg.solve(Gk, gk)
if np.linalg.norm(gk) < epsilon:
break
mk = armijo(fun,x0,dk,gk,beta,sigma,maxm=20)
x0 = x0 + beta ** mk * dk
k = k + 1
graphArr.append(x0)
x = x0
val = fun(x)
return k, x, val, graphArr
k, x, val, graphArr = dampingNewton(fun, gradient, hessianMatrix, np.array([-1.2, 1.0]), 1e-5)
print "minPoint:", x, "minValue:", int(val), "iteration:", k
arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()
结果
minPoint: [ 0.9994168 0.99883101] minValue: 0 iteration: 5000
拟牛顿法
解决无约束优化问题,在计算过程中不需要计算目标函数的Hesse阵,却在某种意义下具有使用Hesse阵时的功效。其基本思想是,用迭代点处的一阶导数和拟牛顿条件来构造目标函数的曲率近似,然后把极小值作为新的迭代点,并重复这一过程,直至求得满足精度的近似极小值点。
流程
实现: 对称秩1
import numpy as np
from matplotlib import pyplot
from computeGradient import computeGrad
def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
"""
非精确线搜索Armijo准则
:param beta:(0,1)
:param sigma:(0,0.5)
:param maxm:最大m值(m 非负整数0,1,2,3,...)
:return:mk
"""
m = 0
mk = 0
while m < maxm:
if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
mk = m
break
m += 1
return mk
def sr1Newton(fun, gfun, x0, epsilon, maxIteration):
"""
基于Armijo的对称秩1算法
:param fun:目标函数
:param gfun:梯度
:param x0:初始点
:param epsilon:容许误差
:param maxIteration:最大迭代次数
:return:k 迭代次数,x 近似最优点,val 近似最优值
"""
beta = 0.55
sigma = 0.4
k = 0
n = len(x0)
Hk = np.eye(n)
graphArr = []
graphArr.append(x0)
while k < maxIteration:
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
dk = -1.0 * np.dot(Hk, gk)
mk = armijo(fun, x0, dk, gk, beta, sigma, maxm=20)
x = x0 + beta ** mk * dk
sk = x - x0
yk = gfun(x0) - gk
sHy = sk - np.dot(Hk, yk)
Hk = Hk + 1.0 * sHy * sHy.reshape((n, 1)) / 1.0 * sHy.reshape((n, 1)) * yk
k = k + 1
x0 = x
graphArr.append(x0)
x = x0
val = fun(x0)
return k, x, val, graphArr
k, x, val, graphArr = sr1Newton(fun, gradient, np.array([-1.0, 1.0]), 1e-5, 10000)
print "minPoint:", x, "minValue:", int(val), "iteration:", k
arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()
结果
minPoint: [ 0.99998982 0.99997959] minValue: 0 iteration: 8121
实现: BFGS
# coding=utf-8
import numpy as np
from matplotlib import pyplot
from computeGradient import computeGrad
def fun(x):
return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
def gradient(x):
# grad = np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])
grad = computeGrad(fun, x, 1e-5)
return grad
def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
"""
非精确线搜索Armijo准则
:param beta:(0,1)
:param sigma:(0,0.5)
:param maxm:最大m值(m 非负整数0,1,2,3,...)
:return:mk
"""
gamma = 0.7
m = 0
mk = 0
while m < maxm:
# gk1 = gfun(x0 + beta ** m * dk)
# if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk) and np.dot(gk1.T,dk) >= gamma * np.dot(gk.T, dk):
if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
mk = m
break
m += 1
return mk
def bfgs(fun, gfun, x0, epsilon, maxIteration=1000):
beta = 0.55
sigma = 0.4
gamma = 0.7
k = 0
n = np.shape(x0)[0]
Bk = np.eye(n)
graphArr = []
graphArr.append(x0)
while k < maxIteration:
gk = gfun(x0) # 计算梯度
if np.linalg.norm(gk) < epsilon: # 检验终止准则
break
# 通过解方程组 Bk * dk = -gk ,计算搜索方向 dk = -Bk/gk
dk = -1.0 * np.linalg.solve(Bk, gk)
mk = armijo(fun, x0, dk, gk, beta, sigma, maxm=20)
# BFGS
x = x0 + beta ** mk * dk
# print "Iteration" + str(k) + ":" + str(x)
sk = x - x0
yk = gfun(x) - gk
if np.dot(sk, yk) > 0:
Bs = np.dot(Bk, sk)
ys = np.dot(yk, sk)
sBs = np.dot(np.dot(sk, Bk), sk)
Bk = Bk - 1.0 * Bs.reshape((n, 1)) * Bs / sBs + 1.0 * yk.reshape((n, 1)) * yk / ys
k += 1
x0 = x
graphArr.append(x0)
return x0, fun(x0), k, graphArr
x0, fun0, k, graphArr = bfgs(fun, gradient, np.array([-1.2, 1.0]), 1e-5, 10000)
print "minPoint:", x0, "minValue:", int(fun0), "iteration:", k
arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()
结果
minPoint: [ 0.99999999 0.99999999] minValue: 0 iteration: 32