DFP算法_python

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

1.代码

代码如下（示例）：

https://blog.csdn.net/qq_39422642/article/details/78876143

# 创作者：眰恦
# 地  点：桂林
# 时  间：2021/10/31  21:58


import numpy as np

# 函数表达式
fun = lambda x:100*(x[0]**2 - x[1])**2 +(x[0] - 1)**2

# 梯度函数
gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1),-200*(x[0]**2-x[1])])

# Hessian矩阵
hess = lambda x:np.array([[1200*x[0]**2-400*x[1]+2,-400*x[0]],[-400*x[0],200]])
'''
    功能：用DFP算法求解无约束问题
    输入：x0为初始点、
    输出：x，val分别是最优点，最优解，k是迭代次数
'''
def nwton(fun,gfun,hess,x0):
    maxk = 1e5
    rho = 0.05  # [0,1]
    sigma = 0.4 # [0,0.5]
    epsilon = 1e-5  # 迭代停止条件
    k = 0
    n = np.shape(x0)[0]
    # print(n)
    # Hk = np.linalg.inv(hess(x0))    # 将hess矩阵初始化为单位矩阵
    Hk = np.eye(2)

    while k < maxk:
        gk = gfun(x0)
        if np.linalg.norm(gk) < epsilon:
            break
        dk = -1.0*np.dot(Hk,gk)
        m = 0
        mk = 0
        while m < 20:
            if fun(x0+rho**m*dk) <= fun(x0) + sigma*rho**m*np.dot(gk,dk):
                mk = m
                break
            m += 1

        # DFP校正
        x = x0+rho**mk*dk
        print("第"+str(k)+"次迭代结果为："+str(x))
        sk = x-x0
        yk = gfun(x)-gk

        if np.dot(sk,yk) > 0:
            Hy = np.dot(Hk,yk)
            sy = np.dot(sk,yk)
            yHy = np.dot(np.dot(yk,Hk),yk)
            Hk = Hk-1.0*Hy.reshape((n,1))*Hy/yHy+1.0*sk.reshape((n,1))*sk/sy
            # print(Hy.reshape((n,1)))
        k += 1
        x0 = x
    return x0,fun(x0),k

x0,fun0,k=nwton(fun,gfun,hess,np.array([0,0]))
print(x0,fun0,k)

2.更新

代码如下（示例）：

---

总结

> 初始对称正定矩阵取单位矩阵
> Hk = np.eye(2)
> 第59次迭代结果为：[0.99279292 0.98550795]
> 第60次迭代结果为：[0.99893754 0.99786139]
> 第61次迭代结果为：[0.99997475 0.99994266]
> 第62次迭代结果为：[0.99999988 0.99999889]
> 第63次迭代结果为：[1.00000006 1.00000012]
[1.00000006 1.00000012] 4.698292481029117e-15 64

---

> 初始对称正定矩阵取Hesse的逆矩阵
> Hk = np.linalg.inv(hess(x0))
> 第59次迭代结果为：[0.98564979 0.9713031 ]
> 第60次迭代结果为：[0.99800673 0.99583304]
> 第61次迭代结果为：[0.99970966 0.99937232]
> 第62次迭代结果为：[0.99997932 0.99995891]
> 第63次迭代结果为：[1.0000007  1.00000134]
> 第64次迭代结果为：[1. 1.]
[1. 1.] 4.699709099980616e-18 65