0618法、最速下降法、牛顿法解方程最优解的python实现

大家好，欢迎来到IT知识分享网。

最优化理论（二）0618法、最速下降法、牛顿法（python实现）

一、0618法

0618法代码如下：

# @Time : 2023/5/4 # @title: 0618法 # 写传入的函数 def hanshu(x): return 2*x2-x-1 def golden_section(f,a,b,epsilon): lamuda = a+0.382*(b-a) mu = a+ 0.618 *(b-a) k = 0 result = 0 while True: if hanshu(lamuda) < hanshu(mu): if b - lamuda <= epsilon: result = mu break else: a = lamuda lamuda = mu mu = a + 0.618 * (b-a) else: if mu - a < epsilon: result = lamuda break else: b = mu mu = lamuda lamuda = a + 0.382 * (b-a) print(a,b) k+=1 return result if __name__ == "__main__": x = golden_section(hanshu,-1,1,0.16) print("%3.2f" % x) 

二、最速下降法（二元方程）

代码如下

from sympy import * x_1 = symbols('x_1') x_2 = symbols('x_2') # 函数求解 fun = 2*x_12 + x_22 print(fun) # 求导 grad_1 = diff(fun, x_1) grad_2 = diff(fun, x_2) # 精度0.1 MaxIter = 100 epsilon = 0.1 # 定义起始点 x_1_value = 1 x_2_value = 1 iter_cnt = 0 current_step_size = 10000 grad_1_value = (float)(grad_1.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf()) grad_2_value = (float)(grad_2.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf()) current_obj = fun.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf() print( '迭代次数: %2d 当前点 (%3.2f, %3.2f) 当前数值: %5.4f 梯度x1: %5.4f 梯度x2 : %5.4f 步长 : %5.4f' % ( iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size)) while (abs(grad_1_value) + abs(grad_2_value) >= epsilon): iter_cnt += 1 # find the step size t = symbols('t') x_1_updated = x_1_value + grad_1_value * t x_2_updated = x_2_value + grad_2_value * t Fun_updated = fun.subs({ 
   x_1: x_1_updated, x_2: x_2_updated}) grad_t = diff(Fun_updated, t) t_value = solve(grad_t, t)[0] # solve grad_t == 0 # 更新点的位置 grad_1_value = (float)(grad_1.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf()) grad_2_value = (float)(grad_2.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf()) x_1_value = (float)(x_1_value + t_value * grad_1_value) x_2_value = (float)(x_2_value + t_value * grad_2_value) current_obj = fun.subs({ 
   x_1: x_1_value, x_2: x_2_value}).evalf() current_step_size = t_value print( '迭代次数: %2d 当前点 (%3.2f, %3.2f) 当前数值: %5.4f 梯度x1: %5.4f 梯度x2 : %5.4f 步长 : %5.4f' % ( iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size)) print("实际上，问题的最优解x*=(%3.2f, %3.2f)^T" % (int(x_1_value), int(x_2_value))) 

运行结果：

2*x_12 + x_22 迭代次数: 0 当前点 (1.00, 1.00) 当前数值: 3.0000 梯度x1: 4.0000 梯度x2 : 2.0000 步长 : 10000.0000 迭代次数: 1 当前点 (-0.11, 0.44) 当前数值: 0.2222 梯度x1: 4.0000 梯度x2 : 2.0000 步长 : -0.2778 迭代次数: 2 当前点 (-0.11, 0.44) 当前数值: 0.2222 梯度x1: -0.4444 梯度x2 : 0.8889 步长 : 0.0000 迭代次数: 3 当前点 (0.07, 0.07) 当前数值: 0.0165 梯度x1: -0.4444 梯度x2 : 0.8889 步长 : -0.4167 迭代次数: 4 当前点 (0.07, 0.07) 当前数值: 0.0165 梯度x1: 0.2963 梯度x2 : 0.1481 步长 : 0.0000 迭代次数: 5 当前点 (-0.01, 0.03) 当前数值: 0.0012 梯度x1: 0.2963 梯度x2 : 0.1481 步长 : -0.2778 迭代次数: 6 当前点 (-0.01, 0.03) 当前数值: 0.0012 梯度x1: -0.0329 梯度x2 : 0.0658 步长 : 0.0000 实际上，问题的最优解x*=(0.00, 0.00)^T 

三、牛顿法（二元方程）

代码如下：

针对上面的例题，由于求导后的数值在程序的Hesse matrix中无法表示为浮点数，而是未知数。所以导致不能对hesse矩阵进行取逆，这里我直接输入了一个求导后的数值。

from sympy import * import numpy as np # 定义符号 x1,x2 = symbols('x1, x2') # 定义所求函数 f = (x1-1)4 + x22 # 求解梯度值 def get_grad(f, X): # 计算一阶导数 f1 = diff(f, x1) f2 = diff(f, x2) # 代入具体数值计算 grad = np.array([[f1.subs([(x1, X[0]), (x2, X[1])])], [f2.subs([(x1, X[0]), (x2, X[1])])]]) return grad # 求解Hession矩阵 def get_hess(f, X): # 计算二次偏导 f1 = diff(f, x1) f2 = diff(f, x2) # f11 = diff(f,x1,2) f22 = diff(f,x2,2) f12 = diff(f1,x2) f21 = diff(f2,x1) # 计算具体数值计算 hess = np.array([[12, f12.subs([(x1,X[0]), (x2,X[1])])], [f21.subs([(x1,X[0]), (x2,X[1])]), f22.subs([(x1,X[0]), (x2,X[1])])]]) # 转换数值类型为了后续求逆矩阵 hess = np.array(hess, dtype = 'float') return hess # 牛顿迭代 def newton_iter(X0, err): # 记录迭代次数 count = 0 X1 = np.array([[0],[0]]) while True: # 得到两次迭代结果差值 X2 = X0 - X1 if sqrt(X2[0]2 + X2[1]2) <= err: break else: hess = get_hess(f, X0) # 得到Hession矩阵的逆 hess_inv = np.linalg.inv(hess) grad = get_grad(f, X0) X1 = X0 #保存上次迭代结果 # 迭代公式 X0 = X1 - np.dot(hess_inv, grad) count += 1 print('第',count,'次迭代:','x1=',X0[0],'x2=',X0[1]) print('迭代次数为:',count) print('迭代结果为',X0) # 设置初始点 X0 = np.array([[1],[1]]) err = 1e-10 newton_iter(X0, err) 

运行结果

第 1 次迭代: x1= [1] x2= [0] 第 2 次迭代: x1= [1] x2= [0] 迭代次数为: 2 迭代结果为 [[1] [0]]