本文的代码与课本算法的流程框架完全一致,不再重复说明理论,只是汇总一下代码实现。
文章目录
拟牛顿法-BFGS
# -*- coding: utf-8 -*-#
# Author: xhc
# Date: 2021-05-28 16:01
# project: 1_zyh
# Name: BFGS.py
import numpy as np
# 定义求解方程
fun = lambda x: 0.5*x[0]**2+x[1]**2-x[0]*x[1]-2*x[0]
# 定义函数 用于求梯度数值 数据放入mat矩阵
def gfun(x):
x = np.array(x)
part1 = x[0][0]-x[1][0]-2
part2 = -1.0*x[0][0]+2*x[1][0]
result = np.mat([part1,part2])
# print(gf) 打印测试
return result
# 传入参数有三个 fun为函数 gfun为梯度 x0为初始值
def BGFS(fun, gfun, x0):
maxk = 5000 # 迭代次数
rho = 0.45 # 步长
sigma = 0.3 # 常数 在 0到1/2之间
epsilon = 1e-6 # 终止条件
k = 0 # 第几次迭代
Bk = np.mat([[1.0,0.0],[0.0,1.0]]) # 初始值为单位矩阵
while k < maxk:
gk = gfun(x0) # (1,2)
if np.linalg.norm(gk) < epsilon: # 范数小于epsilon 终止
break
dk = -Bk.I.dot(gk.T) # Bk*d+梯度=0
m = 0
mk = 0
while m < 30: # Arijo算法求步长
if fun(x0 + (rho**m)*dk) < fun(x0) + sigma * (rho**m)*gk.dot(dk):
mk = m
break
m = m + 1
x = x0 + (rho**mk)*dk # (2,1) # 迭代公式
# 以下就是BFGS修正公式
sk = x - x0 # (2,1)
yk = gfun(x) - gk # (1,2)
if yk * sk > 0:
Bk = Bk - (Bk * sk * sk.T * Bk)/(sk.T*Bk*sk) + (yk.T*yk)/(yk*sk)
k = k + 1
x0 = x
print("--------------------------")
print(