最近在课程实践作业中，要求用最速下降法、牛顿法和拟牛顿法三种方法求解高维一致凸二次优化问题的极小值，网上看到的大部分程序都是手动求好了凸二次函数 f 的偏导然后带进去计算，这样的话限制死了维数和次数，也让程序显得比较笨拙，因此就自己用python从零实现了一下，由于要的急也还有很多改进空间吧…这种优化问题其实用 matlab 会比较方便，因此在 python 里想的也是借鉴 matlab中的符号计算体系，所以基于 sympy 库去实现的。
最速下降法
理论部分：
代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import math
import time
def steepest_descent(f, x): #输入函数f，初始点x，根据最速下降法求解无约束最优化问题
n = len(x) #维数n
diff_f = [] #存储梯度向量
sym_x = [] #存储自变量的符号形式
for i in range(n):
sym = 'x' + str(i+1)
sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
diff = sy.diff(f, sym_x[i]) #计算梯度
diff_f.append(diff)
error = 1 #误差
time = 0 #迭代次数
list_time = [] #存储画图所需横坐标
list_error = [] #存储画图所需纵坐标
while(error > 1e-6) and (time < 100) : #设置终止循环条件
d = [] #搜索方向d(k)
sub_d = {}
for i in range(n):
sub_d[sym_x[i]] = x[i]
for i in range(n):
d.append( - diff_f[i].evalf(subs=sub_d)) #相当于把x1等参数换成了具体数值，代入梯度计算
print("当前搜索方向为：", d)
m = sy.Symbol('m') #步长变量alpha
func = f
for i in range(n):
func = func.subs(sym_x[i], x[i] + m * d[i]) # 实现f(x) = f(x + m * d)，相当于把x1等参数换成了含m变量，此时func是关于m的函数
best_m = sy.solve(sy.diff(func)) #解得最优步长m(k)
print("当前最优步长为：", best_m)
x1 = [] #计算下一个点
for i in range(n):
x1.append(x[i] + best_m[0] * d[i]) #x(k+1) = x(k) + m(k) * d(k)
x = x1 #迭代
res = 0 #计算误差
sub_x = {}
for i in range(n):
sub_x[sym_x[i]] = x[i]
for i in range(n):
res += diff_f[i].evalf(subs=sub_x) ** 2 #同样是传入具体数值代替参数进行计算
error = math.sqrt(res)
print("当前误差：", error)
list_error.append(error)
time += 1
list_time.append(time)
print("循环次数为：",time)
plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
plt.xlabel("Number of iterations") #横坐标名字
plt.ylabel("The error") #纵坐标名字
plt.show()
return x
t0 = time.clock()
x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = sy.symbols('x1, x2, x3, x4, x5, x6, x7, x8, x9, x10')
f = x1**2 + 2*x2**2 + 3*x3**2 + 4*x4**2 + 5*x5**2 + 6*x6**2 + 7*x7**2 + 8*x8**2 + 9*x9**2 + 10*x10**2 - x1 - 9*x2 - 8*x3 - 7*x4 - 6*x5 - 2*x6 - 4*x7 - 4*x8 - x9 - 3*x10
x = steepest_descent(f, [100,0,1,0,100,0,1,0,100,0])
print("解为：", x)
print("所用时间为", time.clock()-t0,"秒")
运行结果为：
牛顿法
理论部分：
代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import numpy as np #导入numpy库，矩阵运算所需
import math
import time
import copy
def newton(f, x): #输入函数f，初始点x，根据牛顿法求解无约束最优化问题
n = len(x) #维数n
diff_f = [] #存储梯度向量
sym_x = [] #存储自变量的符号形式
for i in range(n):
sym = 'x' + str(i+1)
sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
diff = sy.diff(f, sym_x[i]) #计算梯度
diff_f.append(diff)
hesse_matrix = [] #存储海森矩阵
for i in range(n):
row = []#当前行的梯度
for j in range(n):
row.append(sy.diff(diff_f[i], sym_x[j])) #计算梯度
hesse_matrix.append(row)
error = 1 #误差
time = 0 #迭代次数
list_time = [0] #存储画图所需横坐标
list_error = [1] #存储画图所需纵坐标
while(error > 1e-6) and (time < 100) : #设置终止循环条件
g = copy.deepcopy(diff_f) #计算梯度向量g(k)
H = copy.deepcopy(hesse_matrix) #计算海森矩阵H(k)
sub_g = {}
for i in range(n):
sub_g[sym_x[i]] = x[i]
for i in range(n):
g[i] = float(diff_f[i].evalf(subs=sub_g))
for j in range(n):
H[i][j] = float(hesse_matrix[i][j].evalf(subs=sub_g))
print("当前g(k)：", g)
print("当前H(k)：", H)
d = - np.dot(np.linalg.inv(np.array(H)), np.array(g)) #牛顿方向d(k) = -H(k)^-1 * g(k),inv()实现求逆
x += d #迭代，x(k+1) = x(k) + d
res = 0 #计算误差
sub_x = {}
for i in range(n):
sub_x[sym_x[i]] = x[i]
for i in range(n):
res += diff_f[i].evalf(subs=sub_x) ** 2
error = math.sqrt(res)
print("当前误差：", error)
list_error.append(error)
time += 1
list_time.append(time)
print("循环次数为：",time)
plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
plt.xlabel("Number of iterations") #横坐标名字
plt.ylabel("The error") #纵坐标名字
plt.show()
return x
拟牛顿法：
理论部分：
代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import numpy as np #导入numpy库，矩阵运算所需
import math
import time
import copy
def quasi_newton(f, x): #输入函数f，初始点x，根据DFP拟牛顿法求解无约束最优化问题
n = len(x) #维数n
diff_f = [] #存储梯度向量
sym_x = [] #存储自变量的符号形式
for i in range(n):
sym = 'x' + str(i+1)
sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
diff = sy.diff(f, sym_x[i]) #计算梯度
diff_f.append(diff)
D = np.zeros((n, n))
for i in range(n):
D[i][i] = 1 #D(0) = I 单位矩阵
error = 1 #误差
time = 0 #迭代次数
list_time = [] #存储画图所需横坐标
list_error = [] #存储画图所需纵坐标
while(error > 1e-6) and (time < 100) : #设置终止循环条件
g = copy.deepcopy(diff_f) #计算梯度向量g(k)
sub_g = {}
for i in range(n):
sub_g[sym_x[i]] = x[i]
for i in range(n):
g[i] = float(diff_f[i].evalf(subs=sub_g))
d = - D@g #d(k) = - D(k) * g(k),@是numpy中另一种矩阵乘法的实现
print("当前搜索方向为：", d)
m = sy.Symbol('m') #步长变量
func = f
for i in range(n):
func = func.subs(sym_x[i], x[i] + m * d[i]) # f(x) = f(x + m * d)
best_m = sy.solve(sy.diff(func)) #解得最优步长m(k)
print("当前最优步长为：", best_m)
s = float(best_m[0]) * d #s(k) = m(k) * d(k)
x = x + s #x(k+1) = x(k) + s(k)
g1 = copy.deepcopy(g) #计算梯度向量g(k+1)
sub_g = {}
for i in range(n):
sub_g[sym_x[i]] = x[i]
for i in range(n):
g1[i] = float(diff_f[i].evalf(subs=sub_g))
y = np.array(g1) - np.array(g) #y(k) = g(k+1) - g(k)
s.shape, y.shape = (n,1), (n,1) #计算一维数组转置需要补充维度
# D(k+1) = D(k) + s(k)*s(k)^T / s(k)^T*y(k) - D(k)*y(k)*y(k)^T*D(k) / y(k)^T*D(k)*y(k)
D += (np.dot(s, s.transpose()) / np.dot(s.transpose(), y)) - (np.dot(np.dot(np.dot(D, y), y.transpose()), D) / np.dot(np.dot(y.transpose(), D), y))
res = 0 #计算误差
for i in range(n):
res += g1[i]** 2
error = math.sqrt(res)
print("当前误差：", error)
list_error.append(error)
time += 1
list_time.append(time)
print("循环次数为：",time)
plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
plt.xlabel("Number of iterations") #横坐标名字
plt.ylabel("The error") #纵坐标名字
plt.show()
return x
存在的问题
要说略显笨拙的地方首先就是符号计算吧，由于第一次像这样大面积使用sympy库，我总感觉同样的思路有更简单的函数什么的，但又不太熟悉；其次就是函数 f 还是得自己敲，如果维数高一点也很麻烦，但输入系数矩阵 A 我感觉也同样麻烦…
另外，程序中的示例函数是十维二次函数，更高维度自然是没有问题，但是更高次数，比如三次、四次函数在
这里会有点问题，这里输出的best_m是个列表，也可能解为虚数，虚数显然没办法进行下一步计算，会报错。
扫一扫
2012
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
