23、优化算法全解析

优化算法全解析

1. 背景

优化问题在各个科学领域中普遍存在。文献中存在各种优化算法，这些算法的选择取决于待优化函数的一阶和/或二阶导数是否可用，也依赖于待优化的函数以及约束条件的类型。由于最小化与最大化相反，这里主要关注最小化问题。一般来说，有四种类型的目标函数（和约束条件）：
- 线性
- 二次
- 平滑非线性
- 非平滑

无约束的最小化问题是无约束最小化问题。当目标函数是线性且约束条件为线性不等式时，就形成了线性规划；当目标函数和约束条件都是非线性时，则得到非线性规划。最小化算法也会根据问题的性质而有所不同。例如，在无约束情况下，当提供目标函数的梯度和二阶导数时，在多变量情形中可以使用牛顿法；当仅提供一阶导数时，可以应用拟牛顿法；当问题的维度较大时，可以使用共轭梯度算法。

在存在非线性约束时，可以使用一类梯度方法，如缩减和投影梯度方法。当约束条件为线性或非线性不等式时，可以使用凸规划方法。在大多数情况下，可以使用拉格朗日乘数将有约束问题转化为无约束问题。

一维和多维问题通常有很大差异。单变量和双变量最小化问题通常不难解决，因为函数可以绘制和可视化，特别是当函数平滑时。当一阶导数不可用时，可以使用黄金分割等方法。对于多变量且存在多个最小值的情况，问题会变得更加困难。实际上，多变量情况下最小化的主要障碍是局部最小值问题。例如，当函数为二次函数时可以达到全局最小值：
[f (x) = \frac{1}{2}x^T Ax - b^T x + c]
其中(A)是对称矩阵。该二次方程是一个值得关注的典型例子。(f(.))的梯度是(Ax - b)，最优性的必要条件是(\nabla f (x) = 0)。然而，这个线性方程的解只是部分答案。为了得到完整答案，必须在必要条件的解处计算二阶导数，从而得到海森矩阵：
[H = (h_{ij}) = \begin{bmatrix} \frac{\partial^2f}{\partial x_i\partial x_j} \end{bmatrix} = A]
显然，如果(A)是半正定的，必要条件的解就是全局最小值。但一般来说，待最小化的函数是非二次的，需要应用更高级的工具。

2. 单变量优化

一般来说，一维情况是最简单的最小化问题，特别是当目标函数平滑时。

2.1 直接搜索

简单的直接搜索方法基于连续的函数评估，旨在缩小包含最小值的区间长度。最广泛使用的方法有：
- 二分搜索 ：基于连续对包含最小值的区间进行二等分。经过(n)次迭代后，包含最小值的区间长度(I_n)为：
[I_n = \frac{1}{2^{n/2}} I_0]
其中，如果([x_1, x_2])是初始区间，(I_0 = x_2 - x_1)。
- 黄金分割法 ：基于使用两个额外的点(x_3)和(x_4)将初始区间([x_1, x_2])细分为三个子区间，其中(x_1 < x_3 < x_4 < x_2)。例如，如果(f (x_3) \leq f (x_4))，则预计最小值位于([x_1, x_4])内；否则，位于([x_3, x_2])内。迭代过程如下：
[x^{(i)}_3 = \frac{\tau - 1}{\tau} (x^{(i)}_2 - x^{(i)}_1) + x^{(i)}_1]
[x^{(i)}_4 = \frac{1}{\tau} (x^{(i)}_2 - x^{(i)}_1) + x^{(i)}_1]
其中(\tau)是黄金数(\frac{1 + \sqrt{5}}{2})。

2.2 导数方法

当一阶和可能的二阶导数可用时，已知以下两个条件：
[\frac{d}{dx} f (x^ ) = 0]
[\frac{d^2}{dx^2} f (x^ ) > 0]
是(x^*)为(f(.))最小值的充分条件。在这种情况下，最广泛使用的方法基于牛顿算法，也称为牛顿 - 拉夫逊法，旨在基于(\frac{df (x)}{dx})处的切线计算(\frac{df (x)}{dx})的零点。当(\frac{d^2f (x_k)}{dx^2} \neq 0)时，算法为：
[x_{k + 1} = x_k - \frac{df (x_k)/dx}{d^2f (x_k)/dx^2}]
当二阶导数未提供时，可以使用有限差分方案近似式（E.4）的分母：
[x_{k + 1} = x_k - \frac{x_k - x_{k - 1}}{df (x_k) - df (x_{k - 1})} \frac{df (x_k)}{dx_k}]

3. 直接多变量搜索

与一维情况一样，有多变量直接搜索方法和基于梯度的算法。最广泛使用的直接搜索方法包括：

3.1 下山单纯形法

该方法由Nelder和Mead在1965年提出，最初由Spendley等人在1962年描述。该方法基于单纯形，通常顶点相互等距，通过将目标函数值最大的顶点通过对面（即由其余(m)个点形成的超平面）反射到函数值较小的“较低”顶点，从而形成新的单纯形。该方法可用于快速搜索，但在高维情况下可能效率低下。

3.2 共轭方向/鲍威尔法

大多数多变量最小化算法试图找到函数可以沿其最小化的最佳搜索方向。共轭方向法基于最小化二次函数，并且具有二次收敛性。考虑二次函数：
[f (x) = \frac{1}{2}x^T Gx + b^T x + c]
如果(u^T Gv = 0)，则称方向(u)和(v)（关于(G)）是共轭的。该方法基于找到一组相互共轭的搜索方向，沿着这些方向进行最小化。鲍威尔在1964年表明，如果在第(i)次迭代时的一组搜索方向(u^{(i)}_1, \cdots, u^{(i)}_n)被归一化，使得(u^{(i)}_k Gu^{(i)}_k = 1)，(k = 1, \cdots, n)，那么只有当向量相互共轭（线性独立）时，(\det \begin{bmatrix} u^{(i)}_1, \cdots, u^{(i)}_n \end{bmatrix})才会最大化。这提供了一种找到新搜索方向的方法。

鲍威尔算法步骤如下：
1. 初始化(u_i = e_i)，即规范基向量，(i = 1, \cdots, m)。
2. 初始化(x = x_0)。
3. 最小化(f (x_{i - 1} + \lambda u_i))，(x_i = x_0 + \lambda u_i)，(i = 1, \cdots, m)。
4. 设置(u_{i + 1} = u_i)，(i = 1, \cdots, m)，(u_m = x_m - x_0)。
5. 最小化(f (x_m + \lambda u_m))，(x_0 = x_m + \lambda u_m)，然后回到步骤2。

3.3 模拟退火算法

该算法基于统计力学的概念，利用热力学系统在平衡状态下的玻尔兹曼能量分布概率。该方法使用蒙特卡罗模拟生成移动，特别适用于逃离局部最小值。该算法可应用于连续和离散问题。

4. 多变量基于梯度的方法

与直接搜索方法不同，基于梯度的方法使用目标函数的梯度。假设（平滑）目标函数可以近似为：
[f (x + \delta x) = f (x) + g(x)^T \delta x + \frac{1}{2}\delta x^T H\delta x + o(|\delta x|^2)]
其中(g(x) = \nabla f (x))，(H = \begin{bmatrix} \frac{\partial}{\partial x_i\partial x_j} f (x) \end{bmatrix})分别是(f (x))的梯度向量和海森矩阵。梯度方法也属于下降算法类，在不同迭代中对所需最小值的近似以加法方式进行扰动：
[x_{m + 1} = x_m + \lambda u]
下降算法的区别在于搜索方向(u)的选择方式。大多数梯度方法使用梯度作为搜索方向，因为梯度(\nabla f (x))指向函数增长最快的方向。

4.1 最速下降法

最速下降法使用(u = -\frac{g}{|g|})，并选择(\lambda)来最小化单变量目标函数：
[h(\lambda) = f (x_m + \lambda u)]
迭代(m + 1)的解为：
[x_{m + 1} = x_m - \lambda \frac{\nabla f (x_m)}{|\nabla f (x_m)|}]
当使用式（E.7）时，式（E.9）是二次的，此时解为：
[\lambda = \frac{|\nabla f (x_m)|^3}{\nabla f (x_m)^T H\nabla f (x_m)}]
由于每一步都进行一维最小化，该方法计算成本较高。一些作者使用递减步长选择(\lambda = \alpha^k)，((0 < \alpha < 1))，直到找到第一个使(f(.))下降的(k)。

4.2 牛顿 - 拉夫逊法

这是一维牛顿 - 拉夫逊法的推广，基于最小化二次形式。搜索方向为：
[u = -H^{-1}\nabla f (x)]
在第((m + 1))次迭代时，最小值的近似为：
[x_{m + 1} = x_m - H^{-1}(x_m)\nabla f (x_m)]
也可以选择(x_{m + 1} = x_m - \lambda H^{-1}\nabla f (x_m))，其中(\lambda)可以通过一维最小化找到。牛顿法在每次迭代时需要计算海森矩阵的逆，对于大型问题计算成本很高。此外，该方法只有在海森矩阵正定的情况下才能保证收敛到最小值。同样，最速下降法也存在收敛速度慢的问题。这些缺点促使了更先进和改进算法的发展，主要包括共轭梯度法和拟牛顿法。

4.3 共轭梯度法

下降方向(-g = -\nabla f (x))与到最小值的方向可能接近正交，这可以解释最速下降法收敛速度慢的原因。对于二次函数，例如，最佳搜索方向与上一步的方向共轭。这是共轭梯度法的基本思想，新的搜索方向被构造为与上一步的梯度共轭。该方法可以看作是共轭性与最速下降法的结合，也称为弗莱彻 - 里夫斯（或投影）方法。从一组共轭梯度(-g_k)，(k = 1, \cdots, m)，通过线性组合形成一组新的共轭方向：
[u_k = -g_{k - 1} + \sum_{j = 1}^{k - 1} \alpha_{jk}u_j]
其中(\alpha_{jk} = -\frac{g_{k - 1}^T H u_j}{u_j^T H u_j})，(j = 1, \cdots, k - 1)，(g_k = \nabla f (x_k))。在二次形式中，经过推导可以简化为：
[u_k = -g_{k - 1} + \frac{g_{k - 1}^T g_{k - 1}}{g_{k - 2}^T g_{k - 2}} u_{k - 1}]
对于二次函数，该算法最多在(n)次迭代内收敛，其中(n)是问题的维度。对于一般函数，每次迭代都可以使用式（E.16）更新搜索方向，并且实际上每(n)次迭代后将(u_k)重置为(-g_{k - 1})。

4.4 拟牛顿法

牛顿 - 拉夫逊方向(-H^{-1}g)可以看作是对最速下降方向(-g = -\nabla f (x))的改进。拟牛顿方法试图利用最速下降法和基本二阶牛顿法的二次收敛速度。它基于近似海森矩阵(H)的逆。最广泛使用的拟牛顿方法是戴维森 - 弗莱彻 - 鲍威尔法，有时也称为变度量法，通过迭代过程近似(H^{-1})：
[x_{k + 1} = x_k - \lambda_k S_k g_k]
其中(S_k)是收敛到(H^{-1})的序列，由下式给出：
[S_{k + 1} = S_k - \frac{S_k \delta g_k \delta g_k^T S_k}{\delta g_k^T S_k \delta g_k} + \frac{\delta x_k \delta x_k^T}{\delta x_k^T \delta g_k}]
其中(\delta g_k = g_{k + 1} - g_k)，(\delta x_k = x_{k + 1} - x_k = -\lambda_k S_k \delta g_k)。

4.5 基于常微分方程的方法

基于求解常微分方程组的优化技术已经被提出并使用了一段时间，旨在通过跟踪常微分方程组的轨迹来找到(\min_x F(x))的解。例如，如果(x^ )是式（E.21）的解，那么(\nabla F(x^ ) = 0)。因此，通过从合适的初始条件开始积分动态系统：
[\frac{dx}{dt} = -\nabla F(x)]
原则上应该收敛到(x^*)。该方法可以看作是最速下降算法的连续版本。实际上，当(\frac{dx}{dt})用简单有限差分(\frac{x_{t + h} - x_t}{h})近似时，式（E.22）等价于最速下降算法。式（E.22）也可以用连续牛顿方程代替：
[\frac{dx}{dt} = -H^{-1} (x) \nabla F(x)]
其中(H)是(F(.))在(x)处的海森矩阵。

下面用mermaid流程图展示多变量基于梯度的方法的选择流程：

graph TD;
    A[开始] --> B{是否有梯度信息};
    B -- 是 --> C{问题维度大吗};
    C -- 是 --> D[共轭梯度法];
    C -- 否 --> E{二阶导数信息可用吗};
    E -- 是 --> F[牛顿 - 拉夫逊法];
    E -- 否 --> G[拟牛顿法];
    B -- 否 --> H[直接多变量搜索方法];
    D --> I[迭代求解];
    F --> I;
    G --> I;
    H --> I;
    I --> J[结束];

5. 约束最小化

5.1 背景

约束最小化问题比无约束问题更为复杂。一个典型的（平滑）约束最小化问题形式如下：
[
\begin{align }
\min_x f(x)\
\text{s.t. } g_i(x) = 0, &\quad i = 1, \cdots, r\
h_j(x) \leq 0, &\quad j = 1, \cdots, m
\end{align }
]
当上述方程中涉及的函数为凸函数或多项式时，该问题被称为数学规划。例如，如果(f(.))是二次或凸函数，且约束条件为线性，则存在有效的编程程序进行最小化。

大多数算法尝试将上述问题转化为无约束问题。当约束条件简单时，可通过变量变换轻松实现。示例如下：
- 对于(x \geq 0)形式的约束，变量变换为(x = y^2)。
- 对于(a \leq x \leq b)，可令(x = \frac{a + b}{2} + \frac{b - a}{2} \sin y)。

在许多情况下，不等式(h(x) \leq 0)可通过引入松弛变量(y)处理，得到(h(x) + y^2 = 0)。等式约束通常通过引入拉格朗日乘数来处理。在一些正则条件下，(x^ )成为上述约束局部最小值的必要条件是存在拉格朗日乘数(u^ = (u_1^ , \cdots, u_r^ )^T)和(v^ = (v_1^ , \cdots, v_m^ )^T)，使得：
[
\begin{align }
\nabla f(x^ ) + \sum_{i = 1}^{r} u_i^ \nabla g_i(x^ ) + \sum_{j = 1}^{m} v_j^ \nabla h_j(x^ ) &= 0\
v_j^ h_j(x^ ) &= 0, \quad j = 1, \cdots, m\
v_j^ &\geq 0, \quad j = 1, \cdots, m
\end{align }
]
这些条件被称为库恩 - 塔克最优性条件，表达了拉格朗日函数：
[L(x; u, v) = f(x) + \sum_{i = 1}^{r} u_i g_i(x) + \sum_{j = 1}^{m} v_j h_j(x)]
在(x^ )处对于最优值(u^ )和(v^ )的平稳性。上述第一个向量方程可通过最小化其元素的平方和来求解，即(\min \sum_{k = 1}^{n} \left(\frac{\partial L}{\partial x_k}\right)^2)。在数学规划中，上述系统通常被称为原问题的对偶问题。

5.2 约束最小化的方法

5.2.1 拉格朗日方法

该方法基于在每次迭代时最小化拉格朗日函数：
[L(x; u, v) = f(x) + u_k^T g(x) + v_k^T h(x)]
在下一步迭代(k + 1)得到最小值(x_{k + 1})。乘数(u_{k + 1})和(v_{k + 1})被视为线性化约束的最优乘数：
[
\begin{align }
g(x_{k + 1}) + (x - x_{k + 1})^T \nabla g(x_{k + 1}) &= 0\
h(x_{k + 1}) + (x - x_{k + 1})^T \nabla h(x_{k + 1}) &\leq 0
\end{align }
]
此方法基于在当前点(x_{k + 1})对约束进行线性化。在大多数迭代技术中，可通过最小化(\sum_{j} h_j(x) + \sum_{i} g_i^2(x))获得初始可行点。

5.2.2 惩罚函数法

惩罚函数法的基本思想是，在搜索某个函数的约束最小值时，常遇到约束形式为(g(x) \leq 0)的情况，每次迭代新形成的(x)需满足这些约束。简单的处理方法是形成(g(x))元素的线性组合(u^T g(x))，即惩罚函数，用于考虑(g(x))的正分量。当(g(x))的相应分量不违反约束（即非正）时，(u)的分量为零；否则为大正数。然后需要最小化原始目标函数和惩罚函数的和，即惩罚目标函数。最小化惩罚目标函数可防止搜索算法选择违反约束的方向。

一般而言，惩罚方法基于顺序最小化以下形式的无约束问题：
[F(x) = f(x) + \sum_{j} w_j G(h_j(x), \rho) + \sum_{i} H(g_i(x, \rho))]
其中(w_j)（(j = 1, \cdots, m)）和(\rho)是在最小化过程中可改变值的参数，通常(\rho)随迭代次数增加而减小到零。(G(.))和(H(.))是惩罚函数。例如，函数(G(u, \rho) = \frac{u^2}{\rho})是广泛使用的惩罚函数之一。当存在不等式约束且(\rho)固定时，障碍函数(G(.))在可行区域内部（(h_j(x) \leq 0)，(j = 1, \cdots, m)）非零，在边界上为无穷大。这使迭代点(x_k)保持在可行集内，随着(\rho \to 0)，接近约束最小值。此类障碍函数的例子包括(\log(-h(x)))和(\frac{\rho}{h^2(x)})。以下惩罚函数也用于上述问题：
[\rho^3 \sum_{j} w_j h_j^2(x) + \frac{1}{\rho} \sum_{i} g_i^2(x)]

5.2.3 梯度投影法

该方法基于将梯度(-g = -\nabla f(x))投影到当前点(x)处与可行集相切的超平面上，以找到搜索方向。将不满足的不等式约束和等式约束在当前点附近线性化，考虑(Kx = 0)，其中(K)是((r + l_1) \times n)矩阵，(l_1)是(h(x) \geq 0)的约束数量。

算法步骤如下：
1. 从可行集中选择(x_0)。
2. 对(g_i())（(i = 1, \cdots, r)）和当前起作用的不等式（即(h_j(x_k) \geq 0)的那些）进行线性化，以计算(K)。
3. 计算投影梯度(u)：
- 由(-g = u + K^T w)和(Kx = 0)，可得(w = -(KK^T)^{-1} Kg)。
- 负投影梯度为(u = -\left(I_n - K^T (KK^T)^{-1} K\right) g)。
4. 最小化单变量函数(f(x_k + \lambda_{k + 1} u))，设(x_{k + 1}^{(1)} = x_k + \lambda_{k + 1} u)。若(x_{k + 1}^{(1)})可行，则设(x_{k + 1}^{(2)} = x_{k + 1}^{(1)})；否则，使用牛顿法的适当版本求解(\frac{1}{\rho} \sum_{j} h_j^2(x))，找到可行区域边界上的点(x_{k + 1}^{(2)})。
5. 若(f(x_{k + 1}^{(2)}) \leq f(x_k))，则设(x_{k + 1} = x_{k + 1}^{(2)})，回到步骤2；否则，回到步骤4，通过生成序列(x_{k + 1}^{(t)} = x_k + \frac{1}{\tau^{t - 2}} \lambda_{k + 1} u)求解(\lambda_{k + 1})，直到(f(x_{k + 1}^{(t)}) \leq f(x_k))满足。
6. 迭代步骤2 - 5，直到获得约束最小值。对应起作用的不等式的最优乘数(v_i)和(u^*)由上述(w)给出。

5.2.4 其他与梯度相关的方法

• 多重梯度求和法 ：搜索方向为(u = -\frac{\nabla f(x_k)}{|\nabla f(x_k)|} \sum_{j} \frac{\nabla h_j(x_k)}{|\nabla h_j(x_k)|})，其中求和是对当前点(x_k)处违反的约束进行的。
• 小步长梯度法 ：搜索方向为(u = -\nabla f(x_k) - \sum_{j = 1}^{m} w_j(x_k) \nabla h_j(x_k))，其中当(h_j(x_k) > 0)时(w_j(x_k) = w)（(w)是适当选择的大常数），否则为零。

基于常微分方程的方法在将问题转化为无约束最小化问题后，也可用于约束最小化。

下面用表格对比几种约束最小化方法的特点：
|方法|优点|缺点|适用场景|
| ---- | ---- | ---- | ---- |
|拉格朗日方法|理论基础完善，可处理多种类型约束|需要线性化约束，计算较复杂|约束条件相对规则，可线性化的问题|
|惩罚函数法|实现简单，可防止违反约束|惩罚参数选择困难，可能影响收敛速度|约束条件简单，对收敛速度要求不高的问题|
|梯度投影法|能有效处理约束，利用梯度信息|计算投影矩阵较复杂，高维问题计算量大|约束条件明确，可线性化的高维问题|
|多重梯度求和法|考虑了违反约束的情况|搜索方向计算依赖违反约束的数量|违反约束较少的问题|
|小步长梯度法|可根据约束情况调整搜索方向|大常数(w)选择困难|约束条件较简单的问题|

下面用mermaid流程图展示约束最小化方法的选择流程：

graph TD;
    A[开始] --> B{约束类型简单吗};
    B -- 是 --> C{可线性化吗};
    C -- 是 --> D[拉格朗日方法];
    C -- 否 --> E[惩罚函数法];
    B -- 否 --> F{是否可利用梯度信息};
    F -- 是 --> G[梯度投影法];
    F -- 否 --> H[其他方法];
    D --> I[迭代求解];
    E --> I;
    G --> I;
    H --> I;
    I --> J[结束];

综上所述，优化算法在不同场景下各有优劣。在实际应用中，需要根据目标函数的性质、约束条件以及问题的维度等因素，综合考虑选择合适的优化算法，以达到高效、准确地求解最小化问题的目的。