Python 工具包 最优化算法 python求解最优化问题

最近在课程实践作业中，要求用最速下降法、牛顿法和拟牛顿法三种方法求解高维一致凸二次优化问题的极小值，网上看到的大部分程序都是手动求好了凸二次函数 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是个列表，也可能解为虚数，虚数显然没办法进行下一步计算，会报错。