最优化原理与算法复习_摩登7注册声学隔声装饰工程

您现在的位置：主页 > 新闻中心 > 公司动态最优化原理与算法复习

最优化原理与算法复习

作者：佚名日期：2024-04-22 阅读( )

线搜索类算法的数学表述为：给定当前迭代点 $x_k$ ，首先通过某种算法选取向量 $d_k$ ，之后确定正数 $α_k$ ，则下一步的迭代点可写作

$x_{k+1}=x_k + α_k d_k$

我们称 $d_k$ 为迭代点 $x_k$ 处的搜索方向， $α_k$ 为相应的步长．这里要求 $d_k$ 是一个下降方向，即

$(d_k)^T \ abla f(x_k)<0$

这个下降性质保证了沿着此方向搜索函数 $f$ 的值会减小．线搜索类算法的关键是如何选取一个好的方向 $d_k \\in R^n$ 以及合适的步长 $α_k$ .

我们讨论如何选取适当的步长.

1.1.1 线搜索准则

下面介绍的线搜索准则都是单调线搜索，还有一些非单调线搜索准则，这里不加介绍

Armijo准则

Armijo准则对步长作了本质要求，即保证每一步迭代充分下降。

定义: 设 $d_k$ 是点 $x_k$ 处的下降方向，若
$f(x_k+\\alpha d_k)\\leq f(x_k) + c_1 \\alpha \ abla f(x_k)^T d_k$
则称步长 $\\alpha$ 满足Armijo准则，其中 $c_1\\in(0,1)$ 是一个常数.

但Armijo准则并不能保证迭代能够顺利进行，例如 $\\alpha=0$ 就是一个满足Armijo准则的步长。因此常常需要其他准则来配合Armijo准则使用。下面列举两种准则。

Goldstein-Armijo准则
为了克服 Armijo 准则的缺陷，我们需要引入其他准则来保证每一步的 $α_k$ 不会太小。
定义: 设 $d_k$ 是点 $x_k$ 处的下降方向，若
$\\begin{align*} f(x_k+\\alpha d_k) &\\leq f(x_k) + c \\alpha \ abla f(x_k)^T d_k \\\\ f(x_k+\\alpha d_k) &\\geq f(x_k) + （1-c） \\alpha \ abla f(x_k)^T d_k \\end{align*}$
则称步长 $\\alpha$ 满足 Goldstein 准则,其中 $c\\in(0,\\frac{1}{2})$ .

# Goldstein-Armijo一维线搜索

Wolfe-Armijo准则

Goldstein准则能够保证函数值充分下降，但是它可能避开了最优函数值，为此我们引入Armijo-Wolfe准则。

定义： 设 $d_k$ 是点 $x_k$ 处的下降方向，若
$\\begin{align*} f(x_k+\\alpha d_k) &\\leq f(x_k) + c_1 \\alpha \ abla f(x_k)^T d_k \\\\ \ abla f(x_k+\\alpha d_k)^T d_k &\\geq c_2 \ abla f(x_k)^T d_k \\end{align*}$
则称步长 $\\alpha$ 满足Wolfe准则，其中 $c_1,c_2 \\in (0,1)$ 为给定的常数且 $c_1<c_2$ .

#Wolfe-Armijo一维搜索

1.1.2 线搜索实现算法

回退法

1.2.1 试探法

黄金分割法

黄金分割法是一种区间收缩法，其主要用于求解单峰函数的极值问题。其利用了单峰函数的性质，排除极值不可能存在的区间，从而达到逐渐缩小目标区间的效果。

求极小值的黄金分割法具体实现流程如下：

给定 $a,b,\\varepsilon >0$ .
计算 $x_1=a + 0.382(b-a),x_2=a + 0.618(b-a)$ .
若 $b-a<\\varepsilon$ ，停止，输出 $x_*=\\frac{a+b}{2}$ ;否则转(4).
计算 $f_1=f(x_1),f_2=f(x_2)$ .
如果 $f_1 > f_2$ ，令 $a=x_1,x_1=x_2,x_2=a+0.618(b-a)$ , 否则令 $b=x_2,x_2=x_1,x_1=a+0.382(b-a)$ ,转(3).

#黄金分割法
def goldenSearch(a,b,fun,epsilon = 0.01):
    x1 = a + 0.382 * (b-a)
    x2 = a + 0.618 * (b-a)
    f1 = fun(x1)
    f2 = fun(x2)
    while b-a > epsilon:
        print('a={},b={},x1={},x2={}'.format(a,b,x1,x2))
        if f1 < f2:
            b = x2
            x2 = x1
            x1 = a + 0.382 * (b - a)
            f1 = fun(x1)
            f2 = fun(x2)
        else:
            a = x1
            x1 = x2
            x2 = a + 0.618 * (b - a)
            f1 = fun(x1)
            f2 = fun(x2)
    return (b+a)/2, fun((b+a)/2)

if __name__ == '__main__':
    fun = lambda x: 2*x**2 - x - 1
    print(goldenSearch(-10,10,fun,0.01))

斐波那契法

斐波那契法是分割法求解一维极小值问题的最优策略，而黄金分割法只是近似最优。

求极小值的斐波那契法实现流程如下。

#斐波那契法
import math
class Fibs:
    def __init__(self):
        self.a = 0
        self.b = 1
    def __next__(self):
        self.a, self.b = self.b, self.a + self.b
        return self.a
    def __iter__(self):
        return self
def FibsSearch(a,b,fun,epsilon = 0.01,delta = 0.001):
    ctrl = math.ceil((b-a)/epsilon)
    i = 0
    fibs = Fibs()
    F = []
    for f in fibs:
        if f > ctrl:
            break
        else:
            F.append(f)
            i += 1
    F.insert(0,0)
    n = len(F) - 1

    x1 = a + F[n-2]/F[n] * (b-a)
    x2 = a + F[n-1]/F[n] * (b-a)
    f1 = fun(x1)
    f2 = fun(x2)
    print(a, '\	', b, '\	', fun(x1), '\	', fun(x2))

    for k in range(1,n-1):
        if k == n-2:
            i = k
            break
        if f1 > f2:
            a = x1
            x1 = x2
            x2 = a + F[n-k-1]/F[n-k] * (b-a)
            f1 = fun(x1)
            f2 = fun(x2)
        else:
            b = x2
            x2 = x1
            x1 = a + F[n-k-2]/F[n-k] * (b-a)
            f1 = fun(x1)
            f2 = fun(x2)
        print(a,'\	', b,'\	', fun(x1),'\	', fun(x2))
    k = i
    x1 = x1
    x2 = x1 + delta
    f1 = fun(x1)
    f2 = fun(x2)
    if f1 > f2:
        a = x1
    else:
        b = x2
    print(a,'\	',b,'\	',fun(x1),'\	',fun(x2))
    return a,b,fun(x1),fun(x2)
if __name__ == '__main__':
    fun = lambda x: 2 * x ** 2 - x - 1
    print(FibsSearch(-1,1,fun,0.16,0.01))

1.2.2 函数逼近法

牛顿法
我们首先介绍解方程的牛顿迭代法，即常说的切线迭代法，其选取一个可微函数上的初始点 $x_1$ ,在其处作切线，在切线与 $x$ 轴交点处作垂线与函数得到一个新的交点作为 $x_2$ ，如此不断迭代。
容易得到牛顿迭代法的迭代公式为
$x_{k+1}=x_k - \\frac{f(x_k)}{f'(x_k)}$
而极值优化中的牛顿法则可以看作是对 $f'$ 通过切线法求解零点，从而求解函数的驻点，这要求函数是二阶可微的。
很明显牛顿迭代法是一个局部算法，初始点的选取对算法的收敛性和收敛速度有很大影响
牛顿迭代法的计算步骤：

取初始点 $x_0$ ,允许误差 $\\varepsilon$ ,计数器 $k=0$ .
若 $\\parallel f'(x_k) \\parallel < \\varepsilon$ ，停止迭代，得到 $x_*=x_k$ .
否则，根据迭代公式
$x_{k+1}=x_k - \\frac{f'(x_k)}{f''(x_k)}$
计算 $x_{k+1}$ ,置 $k=k+1$ ，回到（2）.

#牛顿法
def NewtonSearch(x0,fun,Dfun,D2fun,epsilon = 0.01):
    while abs(Dfun(x0)) >= epsilon:
        x0 = x0 - Dfun(x0)/D2fun(x0)
        print(x0)
    return x0,fun(x0)
if __name__ == '__main__':
    fun = lambda x: 2* x**4 - x - 1
    Dfun = lambda x:8* x**3 - 1
    D2fun = lambda x:24* x**2
    print(NewtonSearch(6,fun,Dfun,D2fun,0.01))

割线法

和牛顿法类似，割线法不断用割线逼近目标函数的导数曲线，并把割线的零点作为目标函数零点的近似。

割线法的实现流程如下：

给定初始点 $x_0,x_1$ ,允许误差 $\\varepsilon >0$ ,计数器 $k=1$ .
若 $\\parallel \ abla f(x_k) \\parallel < \\varepsilon$ ,则停止迭代，得到 $x_*=x_k$ .
根据迭代公式

$x_{k+1}=x_k - \\frac{x_k - x_{k-1}}{f'(x_k)-f'(x_{x-1})}\ abla f(x_k)$

计算新的迭代点,置 $k=k+1$ ,转(2).

抛物线法

基本思想为在极小值点附近，选取三个初始点 $x_1,x_2,x_3$ ，满足 $f(x_1)>f(x_2)$ , $f(x_2)<f(x_3)$ ，采用二次插值法进行插值，并求插值函数极小值，然后再根据算法进行迭代。

在插值过程中，我们可以采用Lagrange插值法，直接得到插值多项式

$\\begin{align*}\\varphi(x)&=\\sum_{i=1}^3 f(x_i)(\\prod_{j=1\\\\j\ ot=i}^{3}\\frac{x-x_j}{x_i-x_j}) \\end{align*}$

对其求导

$\\begin{align*}\\varphi'(x)&=\\sum_{i=1}^3 f(x_i)(\\prod_{j=1\\\\j\ ot=i}^{3}\\frac{x-x_j}{x_i-x_j})' \\\\&=\\sum_{i=1}^3 f(x_i)(\\sum_{1\\leq j,k\\leq 3\\\\j\ ot=k\\\\j,k\ ot=i}\\frac{x-x_j}{(x_i-x_j)(x_i-x_k)})\\end{align*}$

得到 $\\varphi(x)$ 的极小值点

$x_k=\\frac{(x_2^2-x_3^2)f(x_1)+(x_3^2-x_1^2)f(x_2)+(x_1^2-x_2^2)f(x_3)}{2[(x_2-x_3)f(x_1)+(x_3-x_1)f(x_2)+(x_1-x_2)f(x_3)]}$

置

$\\begin{align*}& x_2=\\hat{x}:\\hat{x}=x_2 or \\hat{x}=x_k,f(\\hat{x})=\\min\\{f(x_2),f(x_k)\\}\\\\& x_1=left\\quad x_2\\\\&x_3=right\\quad x_2\\end{align*}$

即 $x_1,x_3$ 分别为 $x_2$ 左右第一点，从而达到缩小搜索区间的目的；同样的，当误差 $\\parallel f'(x_k-x_{k-1})\\parallel < \\varepsilon or \\parallel x_k - x_{k-1}\\parallel <\\delta$ ,终止迭代并输出 $x_*=x_k$ .

import numpy as np
import matplotlib.pyplot as plt
import math


def phi(x):
    '''
        测试函数1
    :param x:
    :return:
    '''
    return x * x - 2 * x + 1


def complicated_func(x):
    '''
        测试函数2
    :param x:
    :return:
    '''
    return x * x * x + 5 * math.sin(2 * x)


def parabolic_search(f, a, b, epsilon=1e-1):
    '''
        抛物线法，迭代函数
    :param f: 目标函数
    :param a:   起始点
    :param b:   终止点
    :param epsilon: 阈值
    :return:
    '''
    h = (b - a) / 2
    s0 = a
    s1 = a + h
    s2 = b
    f0 = f(s0)
    f1 = f(s1)
    f2 = f(s2)
    h_mean = (4 * f1 - 3 * f0 - f2) / (2 * (2 * f1 - f0 - f2)) * h
    s_mean = s0 + h_mean
    f_mean = f(s_mean)
    # 调试
    k = 0
    while s2 - s0 > epsilon:
        h = (s2 - s0) / 2
        h_mean = (4 * f1 - 3 * f0 - f2) / (2 * (2 * f1 - f0 - f2)) * h
        s_mean = s0 + h_mean
        f_mean = f(s_mean)
        if f1 <= f_mean:
            if s1 < s_mean:
                s2 = s_mean
                f2 = f_mean
                # 重新计算一次，书上并没有写，所以导致一直循环
                s1 = (s2 + s0)/2
                f1 = f(s1)
            else:
                s0 = s_mean
                f0 = f_mean
                s1 = (s2 + s0)/2
                f1 = f(s1)
        else:
            if s1 > s_mean:
                s2 = s1
                s1 = s_mean
                f2 = f1
                f1 = f_mean
            else:
                s0 = s1
                s1 = s_mean
                f0 = f1
                f1 = f_mean
        # print([k, (s2 - s0), f_mean, s_mean])
        print(k)
        k += 1
        plt.scatter(s_mean, f_mean)
    return s_mean, f_mean


if __name__ == '__main__':
    x = np.linspace(1, 3, 200)
    y = []
    index = 0
    for i in x:
        y.append(complicated_func(x[index]))
        index += 1
    plt.plot(x, y)

    result = parabolic_search(complicated_func, 1.0, 3.0)
    print(result)
    plt.scatter(result[0],result[1])
    plt.show()

    # x=np.linspace(0, 2, 200)
    # plt.plot(x, phi(x))
    # plt.show()
    # result=parabolic_search(phi, 0, 2.0)
    # print(result)
————————————————
版权声明：本文为CSDN博主「shawn_zhu1」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_26479655/article/details/82889160

对搜索步长的讨论告一段落，下面来讨论如何确定搜索方向

最速下降法的基本思想是选取当前搜索点 $x_k$ 的函数负梯度方向 $-\ abla f(x_k)$ 作为搜索方向 $d_k$ ，在该方向上进行精确一维搜索

$\\alpha_k=\\min_{\\alpha}f(x_k+\\alpha d_k)$

确定步长进行迭代 $x_{k+1}=x_k+a_k d_k$ ，直至 $\\parallel d_k \\parallel <\\varepsilon$ ，这里 $\\varepsilon$ 是提前给定的误差.

特别的，对于二次凸函数

$f(x)=\\frac{1}{2}x^TAx+bx+c，A=A^T,A=B^TB,\\det{B}\ ot=0$

的极小值优化问题,其负梯度方向 $-\ abla f(x_k)=-Ax-b$ ,并且求解一元可微函数 $f(x_k+\\alpha d_k)$ 极值,

$\\frac{df(x_k+\\alpha d_k)}{d\\alpha}=\\frac{d(\\frac{1}{2}x_k^TAx_k+\\alpha x_k^TAd_k+\\frac{1}{2}\\alpha^2 d_k^T A d_k+bx_k+\\alpha b d_k+c)}{d\\alpha}=0$

即求解

$x_k^TAd_k+\\alpha d_k^T A d_k+bd_k=0$

得到

$\\alpha_k=\\frac{d_k^Td_k}{d_k^T A d_k}$

但是由于最速下降法依赖于矩阵的条件数 $r$ ，即最大特征值与最小特征值的比值 $r=\\frac{\\lambda_{max}}{\\lambda_{min}}$ ,条件数越小最速下降法收敛速度越快，反之，当条件数很大时，收敛速度非常慢，且越靠近精确解越慢(会发生锯齿效应)，因而在实际应用中效果并不太良好。

首先需要指出，对于二次凸函数的极小值问题，共轭梯度法本质上是一个求精确值的直接方法，扩张子空间定理保证了它可以在有限步内求得优化问题的精确解，但由于舍入误差的存在，在实际计算中无法直接得到很好的精确解

共轭方向法是在一组关于对称正定矩阵 $A$ 的共轭方向 $d_1,d_2,\\cdots,d_n$ 上进行搜索从而得到最优解的方法。这里对共轭进行一点说明，如果在 $\\mathbb{R}^n$ 上定义内积 $(d_i,d_j)=d_i^T A d_j$ ，两个方向在这个意义下的共轭实际上就是广义内积空间中的正交，从而这样的 $n$ 个共轭向量组成的向量组实际上构成了全空间的一组基。

不加证明的引入扩张子空间定理

扩张子空间定理： 如果 $d_1,d_2,\\cdots,d_n$ 是关于对称正定矩阵 $A$ 的一组共轭方向，且有二次凸函数
$f(x)=\\frac{1}{2}x^TAx+bx+c$
那么以任何初始点 $x_0$ 沿方向 $d_1,d_2,\\cdots,d_k$ 进行 $k$ 次搜索 ,得到的点 $x_k$ 是 $f$ 在线性流形 $x_0+\\mathcal{L}\\{d_1,d_2,\\cdots,d_k\\}$ 上的唯一极小值点。特别的，经过 $n$ 次搜索后，得到的点 $x_n$ 是函数的全局唯一极小值.

基于这一定理，我们只需要找到一组共轭方向，便可以进行最优化求解。根据格莱姆施密特(非标准)正交化方法，这只需要找到一组线性无关的方向，即可进行正交化处理得到一组 $A$ 共轭方向。

共轭梯度法作为一种特殊的共轭方向法，它选取最速下降方向 $d_1=g_1=-\ abla f(x_0)=Ax_0+b$ 为初始方向，并进行迭代

$x_{1}=x_0+\\alpha_1 d_1$

这里步长采用最速下降法公式

$\\alpha_1=\\frac{g_1^Td_1}{d_1^T A d_1}$

从而得到 $x_1$ 处的最速下降方向 $g_2=-\ abla f(x-1)=-Ax_1-b$ ,根据格莱姆施密特方法，对其作 $A$ 共轭处理

$d_2=g_2-\\frac{(g_2,d_1)}{(d_1,d_1)}d_1$

!!请注意这里的内积并非标准内积，而是定义为 $(d_i,d_j)=d_i^T A d_j$

便得到了新的方向，且该方向与 $d_1$ 共轭。

更一般的可以得到迭代公式

$\\begin{align*}& x_{k}=x_{k-1}+\\frac{g_k^T d_k}{(d_k^T,d_k)}d_k\\\\& d_{k}=g_{k}-\\frac{(g_k,d_{k-1})}{(d_{k-1},d_{k-1})}d_{k-1}\\\\\\end{align*}$

#共轭梯度算法
import numpy as np
import matplotlib.pyplot as plt
def testFun(x, y):
    t = 3.0*x - 2.0*y
    t1 = x - 1.0
    z = np.power(t, 2) + np.power(t1, 4)
    return z

def gradTestFun(x, y):
    '''
    求函数的梯度
    '''
    delta_x = 1e-6      #x方向差分小量
    delta_y = 1e-6      #y方向差分小量
    grad_x = (testFun(x+delta_x, y)-testFun(x-delta_x, y))/(2.0*delta_x)
    grad_y = (testFun(x, y+delta_y)-testFun(x, y-delta_y))/(2.0*delta_y)
    grad_xy = np.array([grad_x, grad_y])
    return grad_xy

def getStepLengthByNewton(array_xy, array_d):
    '''
    采用牛顿法，精确线性搜索确定移动步长
    '''
    a0 = 1.0           #初始猜测值
    e0 = 1e-6          #退出搜索循环的条件
    delta_a = 1e-6     #对a作差分的小量
    while(1):
        new_a = array_xy + a0*array_d
        new_a_l = array_xy + (a0-delta_a)*array_d
        new_a_h = array_xy + (a0+delta_a)*array_d
        diff_a0 = (testFun(new_a_h[0], new_a_h[1]) - testFun(new_a_l[0], new_a_l[1]))/(2.0*delta_a)
        if np.abs(diff_a0) < e0:
            break
        ddiff_a0 = (testFun(new_a_h[0], new_a_h[1]) + testFun(new_a_l[0], new_a_l[1]) - 2.0*testFun(new_a[0], new_a[1]))/(delta_a*delta_a)
        a0 = a0 - diff_a0/ddiff_a0
    return a0

def plotResult(array_xy_history):
    x = np.linspace(-1.0, 4.0, 100)
    y = np.linspace(-4.0, 8.0, 100)
    X, Y = np.meshgrid(x, y)
    Z = testFun(X, Y)
    plt.figure()
    plt.xlim(-1.0, 4.0)
    plt.ylim(-4.0, 8.0)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.pcolormesh(X, Y, Z)
    plt.scatter(array_xy_history[:,0], array_xy_history[:,1])
    xy_count = array_xy_history.shape[0]
    plt.show()

def mainFRCG():
    '''
    使用CG算法优化，用FR公式计算组合系数
    '''
    ls_xy_history = []                           #存储初始坐标的迭代结果
    xy0 = np.array([4.0, -2.0])                   #初始点
    grad_xy = gradTestFun(xy0[0], xy0[1])
    d = -1.0*grad_xy                             #初始搜索方向
    e0 = 1e-6                                    #迭代退出条件
    xy = xy0
    while(1):
        ls_xy_history.append(xy)
        grad_xy = gradTestFun(xy[0], xy[1])
        tag_reach = np.abs(grad_xy) < e0
        if tag_reach.all():
            break
        step_length = getStepLengthByNewton(xy, d)
        xy_new = xy + step_length*d
        grad_xy_new = gradTestFun(xy_new[0], xy_new[1])
        b = np.dot(grad_xy_new, grad_xy_new)/np.dot(grad_xy, grad_xy)       #根据FR公式计算组合系数
        d = b*d - grad_xy_new
        xy = xy_new
    array_xy_history = np.array(ls_xy_history)
    plotResult(array_xy_history)
    return array_xy_history
print(mainFRCG())

本文使用 Zhihu On VSCode 创作并发布

上一篇：保健品能在抖音小店卖吗？需要哪些资质？下一篇：交叉熵损失函数和Adam优化器

某某声学隔音工程有限公司

解决方案

摩登7登录

经典案例

新闻中心

关于我们

联系我们

关注我们