最优化算法实践

来源：互联网发布：惊讶猫走红网络编辑：程序博客网时间：2024/05/24 05:17

学习资料：《最优化基础理论与方法》复旦大学出版社

一维搜索

在xk+1=xk+ak∗dk,k=0,1,2,...中，假定在xk处的搜索方向dk已经确定，怎样寻找沿dk方向合适的步长ak，以确定xk+1=xk+ak∗dk，且满足f(xk+1)<f(xk)

精确线性搜索

略

搜索区间与单峰函数

给定一个初始值，找到一个搜索区间（要求函数在搜索区间上是单峰函数）

def searchInterval(func, start, h, gamma):    '''    func: 确定步长的一维函数    start：初始值    h：确定搜索区间的进退法的步长    gamma：进退法中的一个常数，用来增大步长进行搜索，大于1    '''    a0 = start    v0 = func(a0)    a1 = a0 + h    v1 = func(a1)    if v0 > v1:  # v0 > v1 则向右进，直到找到大于v1的为止        a2 = a1  # 初始化a2，方便迭代        k = 1  # 用于计数和作为gamma的指数        while True:            a2 += gamma**k*h             v2 = func(a2)            if v2 > v1:  # 如果找到大于v1的值，则确定了一个搜索区间                return (a0, a2)            a0 = a1  # 没有找到则将a1置为a0，a2置为a1, v2置为v1            a1 = a2            v1 = v2            k += 1    else:  # v0 <= v1 则交换a0, a1的值，向左退，直到找到大于v1的为止        temp = a0; a0 = a1; a1 = temp        temp = v0; v0 = v1; v1 = temp        a2 = a1        k = 1        while True:            a2 -= gamma**k*h  # 向左退            v2 = func(a2)            if v2 > v1:                return (a2, a0)            a0 = a1            a1 = a2            v1 = v2            k += 1

可以看到有很多重复代码，且有无线循环的可能，接下来优化时，加入一个迭代上限itr

def searchInterval2(func, start, h, gamma, itr = 100):    '''    func: 确定步长的一维函数    start：初始值    h：确定搜索区间的进退法的步长    gamma：进退法中的一个常数，用来增大步长进行搜索，大于1    '''    a0 = start    v0 = func(a0)    a1 = a0 + h    v1 = func(a1)    if v0 <= v1:        temp = a0; a0 = a1; a1 = temp        temp = v0; v0 = v1; v1 = temp        sgn = -1    a2 = a1    k = 1    while k < itr:        a2 += sgn*(gamma**k)*h        v2 = func(a2)        if v2 > v1:            return (a0, a2) if sgn == 1 else (a2, a0)        a0 = a1        a1 = a2        v1 = func(a1)        k += 1    print "Can not find a valid interval."

测试：无约束规划问题为：
minimize{f(x)=(x1−2)2+(x1−2x2)2}
则在初始点x1=(0,3)T时，由最速下降法得到的下降方向为d1=−▽f(x1)=(16,−24)T. 所以第一个一维搜索问题就是寻找 α1 使得 f(x1+α1∗d1) 最小。首先要找到α1 的搜索范围，这样才能用直接搜索法得到较精确值。

import numpy as npfrom numpy.linalg import *# 定义目标函数及其梯度def objFunc(x):    x = np.array(x)    A = np.array([[2., -2.],                 [-2., 4.]])    b = np.array([-4., 0.])    return x.T.dot(A).dot(x)+b.T.dot(x)+4def grad(x):    x = np.array(x)    A = np.array([[2., -2.],                 [-2., 4.]])    b = np.array([-4., 0.])    return 2.*A.dot(x)+b.Tleft, right = searchInterval(lambda a: objFunc(np.array([0., 3.])-a*grad([0., 3.])), 0., 0.2, 2.)

结果为(−0.4,0.2)

直接搜索法–0.618法

直接搜索法中最简单的0.618法，不断的计算区间中左右两点的值，以缩小区间到给定精度

from math import exp, sqrt, logdef directSearch(func, a, b, e, al=None, ar=None, lvalue=None, rvalue=None):    TAO = (sqrt(5)-1)/2    if al is None:        al = a + (1-TAO)*(b-a)    if ar is None:        ar = a + TAO * (b-a)    if lvalue is None:        lvalue = func(al)    if rvalue is None:        rvalue = func(ar)    if abs(b - a) <= e:        return al if lvalue < rvalue else ar    if lvalue < rvalue:        b = ar        ar = al        rvalue = lvalue        return directSearch(func, a, b, e, ar=ar, rvalue=rvalue)    else:        a = al        al = ar        lvalue = rvalue        return directSearch(func, a, b, e, al=al, lvalue=lvalue)

测试：对上面给出的区间运用0.618法找到较为精确的解

directSearch(lambda a: objFunc(np.array([0., 3.])-a*grad([0., 3.])), left, right, 1e-4)

结果为0.0956

非精确一维搜索方法

略

无约束规划方法

最速下降法

搜索方向 dk=−▽f(xk) 时，在 xk 处下降的幅度最大。
算法代码如下，其中利用一维搜索的方法获得步长：

def gd(func, grad, start, e):    start = np.array(start)    if sum(grad(start)**2) < e:        return start    d = -grad(start)    left, right = searchInterval(lambda a: func(start+a*d), 0., 0.1, 2)  # h=0.1和gamma=2可以变化    a = directSearch(lambda a: func(start+a*d), left, right, 1e-4)    return np.round(gd(func, grad, start+a*d, e),2)  # 只返回两位数的坐标

测试：现在用可视化的方式观察我们的例子，先画出目标函数的3d图，以及它的等高线图，针对某一初始点，我们用最速下降法得到一个解，并在等高线图中画出迭代轨迹。
在这之前需要保存迭代轨迹，因此，稍微修改最速下降法的函数，再里面添加一个 pts.append(start) ， pts 为全局变量，保存迭代过的点

3d图的绘制：

import numpy as npfrom matplotlib import cmimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Ddelta = 0.2x = np.arange(-10, 10, delta)y = np.arange(-10, 10, delta)X, Y = np.meshgrid(x, y)x=X.flatten()y=Y.flatten()z=np.array([objFunc(z) for z in zip(x, y)])  # 直接由外部objFunc函数定义目标函数Z = z.reshape(np.shape(X))  # 便于绘制等高线图fig = plt.figure()ax = fig.gca(projection='3d')ax.plot_trisurf(x, y, z, cmap=cm.jet, linewidth=0.01)plt.show()

结果：
目标函数
进行最速下降迭代

pts = []gd(objFunc, grad, [0.， 3.], 1e-4)

结果为array([ 2., 1.]) 正是目标函数取得最小时的点
pts 中包含迭代轨迹，则可以画出等高线和迭代轨迹：

plt.contour(X, Y, Z)plt.plot(np.array(pts)[:,0],np.array(pts)[:,1],'o-' )plt.show()

结果：
最速下降法
可以看到，每次下降的方向与上一方向是正交的。

但最速下降法容易造成zig-zag迭代路径，比如，当取初始值为(-10, -10)时，等高线和迭代图为：
最速下降法2

固定步长的问题
如果不用一维搜索决定步长，而采用固定步长，则可能产生如下问题：
- 步长过长：本题中用0.2以上的步长都会导致不收敛
- 步长过短：迭代太慢（下图为步长0.05）

如果能找到合适的步长，也能得到正确的结果，但精度可能受损（下图步长为0.15，最后迭代结果为(1.99, 1.)）
0.15步长
而0.1的步长可能更适合，结果为（2., 1.）
0.1步长

牛顿法

迭代公式为：
xk+1=xk+dk,
其中 dk 是线性方程组 ▽2f(xk)d=−▽f(xk) 的解。算法代码如下：

from numpy.linalg import *def newton(func, grad, hesse, start, e):    start = np.array(start)    if sum(grad(start)**2) < e:        return start    d = solve(hesse(start), -grad(start))    return newton(func, grad, hesse, start+d, e)

对于之前那个题，以(0, 3) 为初始值（不论什么为初始值），可以一步得到最优解，因为目标函数是二次多项式。
牛顿法
如果是非凸函数，则结果一般不会得到精确的最优解。取决于初始值，最后结果可能是极小点，鞍点，或者hesse矩阵是奇异的。

共轭梯度法

设 f(x)=12xTGx+rTx+δ，其中G 是正定矩阵。由扩展子空间定理，若d1,d2,…,dn 为 n 个 G 共轭方向，那么从任意的初始点 x1 出发，至多作 n 次精确一维搜索，就可以得到目标函数唯一的极小点。
PRP算法 代码只针对正定二次函数

from numpy.linalg import *import numpy as npdef cgm(func, grad, G, start, e):    x = np.array(start, dtype='f')    d0 = np.zeros_like(start)    beta = 0.    while sum(grad(x)**2) > e:        gradx = grad(x)        d = -gradx+ beta*d0  # 保证是G的共轭方向        a = -d.T.dot(gradx)/(d.T.dot(G).dot(d))  # 精确一维搜索的解        x = x + a*d        beta = grad(x).T.dot(grad(x))/(gradx.T.dot(gradx))        d0 = d    return x

对于一般的可微函数，共轭方向也是下降方向。但beta如何计算，一维搜索怎么得到，还需要具体问题具体分析。如果n步后，没有终止，则搜索方向应重新开始，成为n步重新开始策略。

0 0