非线性规划——无约束最优化问题精讲

最优化问题的研究历史可以追溯到17世纪的变分法，随着数学、物理学、经济学和计算科学的不断发展，最优化问题逐渐成为一个独立的学科。对于无约束最优化问题的求解，从最早的最速下降法，到后来的牛顿法和共轭梯度法，再到现代的变尺度法和智能算法，发展历程反映了科学技术进步的轨迹。无约束最优化问题也是非线性规划的重要内容，与有约束最优化问题相比，无约束最优化问题的特点是不受约束条件的限制，因此可以在更为宽广的搜索空间内寻找最优解，它在众多实际应用中具有广泛的应用，比如经济学、工程设计、机器学习等领域，这里对常用的无约束解析方法进行介绍。

一、无约束最优化问题

最优化问题的一般形式为：

$$\begin{aligned} & \min f(x) \\ & \text { s.t.} \quad x \in X \end{aligned}$$

其中 $x$ 为决策变量，$f(x)$ 是目标函数，$X$ 为约束集或可行域。特别地，如果 $X=R^n$ ，则最 优化问题成为无约束最优化问题。
最优化方法通常采用迭代法求它的最优解，其基本思想是：给定一个初始点 $x_0$ ，按照某一迭代规则产品一个点列 $\left\{x_n\right\}$ ，使得当 $\left\{x_n\right\}$ 是有穷点列时，其最后一个点是最优化模型问题的最优解。迭代规则由迭代公式决定，迭代公式的基本表示形式如下:

$$x_{k+1}=x_k+\alpha_k d_k$$

式中， $\alpha_k$ 为步长因子， $d_k$ 为搜索方向。在最优化算法中，搜索方向 $d_k$ 是 $f$ 在 $x_k$ 点处的下降方 向, 即:

$$f\left(x_k+\alpha_k d_k\right)<f\left(x_k\right)$$

最优化方法的基本结构如下:

• 给定初始点 $x_0$ ；
• 确定搜索方向 $d_k$ ，即按照一定规则，构造 $f$ 在 $x_k$ 点处的下降方向作为搜索方向；
• 确定步长因子 $\alpha_k$ ，使目标函数有某种意义的下降。

二、最速下降法

最速下降法(Gradient Descent Method)的基本思想是：因为连续函数沿着负梯度方向的下降速度是最快的（这一结论由梯度和方向导数的定义可以推出），所以每次迭代我们都从当前点出发，沿着负梯度方向前进一个最优步长，可以期望能较快逼近函数的极值。

最速下降法流程	函数曲面

最速下降法仅有三个步骤：
设置初始值。设置迭代起点$x_0\in R^n$,允许误差$\epsilon>0$和迭代变量初值$k\leftarrow0$。
检查终止条件。如果$||\nabla f(x_k)||<\epsilon$，停止迭代输出$x_k$作为近似最优解；否则转步骤3。
迭代，通过一维搜索求下一个迭代点。取搜索方向为负梯度方向$d_k=-\nabla f(x)$，求$\lambda_k$使得 $$f(x_k+\lambda_k d_k)=\min_{\lambda\geq0} f(x_k+\lambda d_k)$$ 再令 $$x_{k+1}=x_k+\lambda d_k$$ 转步骤２。

步骤中隐含了条件：函数$f(x)$必须可微，也就是说函数$f(x)$的梯度必须存在，其中最关键的步骤是求解问题
$$\begin{equation} \min_{\lambda\geq0}f(x_k+\lambda d_k) \end{equation}$$
给出它的最优性必要条件
$$\begin{equation} \frac{d f(x_k+\lambda d_k)}{d\lambda}=\nabla f(x_k+\lambda d_k)^T d_k=0 \end{equation}$$
当$f(x)$的形式确定，我们可以通过求解这个一元方程来获得迭代步长$\lambda$。当此方程形式复杂，解析解不存在，我们就需要使用“一维搜索”来求解$\lambda$了。一维搜索是一些数值方法，有0.618法、Fibonacci法、抛物线法等等，这里不详细解释了。

例1： 用最速下降法计算函数 $f\left(x_1, x_2\right)=\left(x_1-2\right)^2+\left(x_2-2\right)^2$ 的极小点和极小值, 设初始点为 $x^0=(0,0)^T$, 迭代误差为 $\varepsilon=0.05$ 。
解: 因为函数 $f\left(x_1, x_2\right)$ 的梯度等于 $\nabla f\left(x_1, x_2\right)=\left(\begin{array}{l}2\left(x_1-2\right) \\ 2\left(x_2-1\right)\end{array}\right)$, 那么在点 $x^0=\left(\begin{array}{l}0 \\ 0\end{array}\right)$ 处的梯 度值等于 $\nabla f\left(x^0\right)=\left(\begin{array}{l}-4 \\ -4\end{array}\right)$
令 $d^0=-\nabla f(0,0)=(4,4)^T$, 则有: $x^0+\alpha d^0=\left(\begin{array}{l}4 \alpha \\ 4 \alpha\end{array}\right), f\left(x^0+\alpha d^0\right)=4(2 \alpha-1)^2$ 。我们需 要求解 $\alpha_0$, 使其满足 $f\left(x^0+\alpha_0 d^0\right)=\min _{\alpha \geq 0}\left\{4(2 \alpha-1)^2\right\}$ 。
令 $\frac{d}{d \alpha} 4(2 \alpha-1)^2=8(2 \alpha-1)=0$, 求得 $\alpha_0=\frac{1}{2}$, 所以下一个迭代点为:

$$x^1=x^0+\alpha_0 d^0=\left(\begin{array}{l} 0 \\ 0 \end{array}\right)+\frac{1}{2}\left(\begin{array}{l} 4 \\ 4 \end{array}\right)=\left(\begin{array}{l} 2 \\ 2 \end{array}\right) \text { 。 }$$

函数 $f\left(x_1, x_2\right)$ 在点 $x^1$ 处的梯度为: $\nabla f\left(x^1\right)=\left(\begin{array}{l}0 \\ 0\end{array}\right)$, 由于 $\left\|\nabla f\left(x^1\right)\right\|=0 \leq 0.05$, 那么 $x^1$ 是局部极小点, 迭代结束。因为 $f\left(x_1, x_2\right)$ 是凸函数, 那么 $x^1$ 也是全局极小点, $f\left(x^1\right)=0$ 是全局极小值。

三、牛顿法

牛顿法是一种利用目标函数的二阶导数信息来快速求解无约束最优化问题的算法。与最速下降法仅依赖梯度不同，牛顿法通过考虑目标函数的曲率（即Hessian矩阵）来调整步长和方向，因此具有更快的收敛速度。牛顿法在目标函数接近二次函数时表现尤其出色，能够以二次收敛的速度逼近最优解。

3.1 牛顿法的基本原理

对于一个目标函数 $f(x)$，假设它在当前点 $x_k$ 处是可微的，并且存在二阶导数。我们通过二次近似来代替原始的目标函数，即在 $x_k$ 处将目标函数用二阶泰勒展开式近似表示为：

$$f(x) \approx f(x_k) + \nabla f(x_k)^T (x - x_k) + \frac{1}{2} (x - x_k)^T \nabla^2 f(x_k) (x - x_k)$$

其中，$\nabla f(x_k)$ 是目标函数在点 $x_k$ 处的梯度，$\nabla^2 f(x_k)$ 是Hessian矩阵，即二阶导数矩阵。

通过最小化上述二次近似表达式，牛顿法的更新公式可以推导为：

$$x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$$

其中，$[\nabla^2 f(x_k)]^{-1}$ 是目标函数在点 $x_k$ 处的Hessian矩阵的逆。

3.2 牛顿法的算法步骤

牛顿法的算法步骤可以总结为以下几步：

• ①初始值设定。选择一个初始点 $x_0$ 和一个收敛精度 $\epsilon$。
• ②计算梯度和Hessian矩阵。在当前点 $x_k$ 处，计算目标函数的梯度 $\nabla f(x_k)$ 和Hessian矩阵 $\nabla^2 f(x_k)$。
• ③判断停止条件。如果梯度的范数$\|\nabla f(x_k)\|$小于预设的精度$\epsilon$，则停止迭代，输出当前点$x_k$作为最优解。
• ④更新点。根据牛顿法的更新公式$x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$计算下一个迭代点。
• ⑤重复迭代。将$x_{k+1}$作为新的当前点，返回步骤②，直到满足停止条件。

3.3 练习

例2：考虑下面二元二次函数 $$f(x_1, x_2) = 3x_1^2 + 2x_1x_2 + x_2^2 - 4x_1 - 2x_2$$ ，给出其最优解。

梯度和Hessian矩阵
我们需要计算目标函数的梯度和Hessian矩阵。

• 梯度 $\nabla f(x_1, x_2)$ 是一个向量，表示函数在 $x_1$ 和 $x_2$ 方向上的导数：

$$\nabla f(x_1, x_2) = \begin{pmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \end{pmatrix} = \begin{pmatrix} 6x_1 + 2x_2 - 4 \\ 2x_1 + 2x_2 - 2 \end{pmatrix}$$

• Hessian矩阵 $\nabla^2 f(x_1, x_2)$ 是一个 2x2 矩阵，表示二阶偏导数：

$$\nabla^2 f(x_1, x_2) = \begin{pmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} \end{pmatrix} = \begin{pmatrix} 6 & 2 \\ 2 & 2 \end{pmatrix}$$

算法步骤

• 初始值设定：设定初始点 $(x_1^0, x_2^0)$，例如 $(x_1^0, x_2^0) = (0, 0)$，并设定收敛精度 $\epsilon = 10^{-6}$。
• 计算梯度和Hessian矩阵：在初始点 $(x_1^0, x_2^0)$ 处计算梯度和Hessian矩阵。
• 判断停止条件：如果梯度的范数 $|\nabla f(x_1^k, x_2^k)|$ 小于预设精度 $\epsilon$，则停止迭代，输出最优解。
• 更新点：通过牛顿法更新公式：

$$\begin{pmatrix} x_1^{k+1} \\ x_2^{k+1} \end{pmatrix} = \begin{pmatrix} x_1^k \\ x_2^k \end{pmatrix} - \left[\nabla^2 f(x_1^k, x_2^k)\right]^{-1} \nabla f(x_1^k, x_2^k)$$

• 重复迭代：返回步骤2，直到满足停止条件。

下面以 $(x_1^0, x_2^0) = (0, 0)$ 为初始点，手动计算前两次迭代：

迭代1：

• 初始点：$(x_1^0, x_2^0) = (0, 0)$
• 梯度 $\nabla f(0, 0) = \begin{pmatrix} -4 \ -2 \end{pmatrix}$
• Hessian矩阵 $\nabla^2 f(0, 0) = \begin{pmatrix} 6 & 2 \\ 2 & 2 \end{pmatrix}$
• 计算Hessian矩阵的逆：

$$\left[\nabla^2 f(0, 0)\right]^{-1} = \frac{1}{8} \begin{pmatrix} 2 & -2 \\ -2 & 6 \end{pmatrix}$$

• 更新点 $(x_1^1, x_2^1)$：

$$\begin{pmatrix} x_1^1 \\ x_2^1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} - \frac{1}{8} \begin{pmatrix} 2 & -2 \\ -2 & 6 \end{pmatrix} \begin{pmatrix} -4 \\ -2 \end{pmatrix} = \begin{pmatrix} 0.5 \\0.5 \end{pmatrix}$$

迭代2：

• 当前点：$(x_1^1, x_2^1) = (0.5, 0.5)$
• 梯度 $\nabla f(0.5, 0.5) = \begin{pmatrix} 0 \ 0 \end{pmatrix}$
  梯度为零，满足停止条件，得出最优解 $(x_1^*, x_2^*) = (0.5, 0.5)$。

import numpy as np

# 定义目标函数的梯度
def grad_f(x1, x2):
    return np.array([6*x1 + 2*x2 - 4, 2*x1 + 2*x2 - 2])

# 定义目标函数的Hessian矩阵
def hessian_f(x1, x2):
    return np.array([[6, 2], [2, 2]])

# 牛顿法实现
def newton_method_2d(x1_0, x2_0, epsilon=1e-6, max_iter=100):
    x1, x2 = x1_0, x2_0
    for i in range(max_iter):
        grad = grad_f(x1, x2)
        hessian = hessian_f(x1, x2)
        
        # 如果梯度的范数小于设定的精度，则停止迭代
        if np.linalg.norm(grad) < epsilon:
            print(f"迭代收敛于第 {i} 次，最优解 (x_1, x_2) = ({x1:.2f}, {x2:.2f})")
            return x1, x2
        
        # 更新点
        hessian_inv = np.linalg.inv(hessian)
        update = np.dot(hessian_inv, grad)
        x1, x2 = np.array([x1, x2]) - update
    
    print("迭代未能收敛")
    return x1, x2

# 初始点
x1_0, x2_0 = 0, 0

# 调用牛顿法
optimal_x1, optimal_x2 = newton_method_2d(x1_0, x2_0)
print(f"最优解：(x_1, x_2) = ({optimal_x1:.2f}, {optimal_x2:.2f})")

四、共轭梯度法

共轭梯度法是一种用于求解大规模线性方程组$Ax = b$ 和无约束最优化问题的有效算法。它特别适合于高维稀疏矩阵的情况，且不需要计算Hessian矩阵的逆。在无约束最优化中，它主要应用于二次型目标函数的求解。

给定一个目标函数： $$f(x) = \frac{1}{2} x^T A x - b^T x$$
其中$A$ 是对称正定矩阵，(b$ 是常数向量。目标是找到使$f(x)$ 最小化的$x$。

共轭梯度法的主要步骤如下：

• 初始化：设定初始点$x_0$，计算初始梯度$r_0 = -\nabla f(x_0) = Ax_0 - b$，设定初始搜索方向$p_0 = r_0$，并设定最大迭代次数和收敛精度。
• ①计算步长$\alpha_k$： $$\alpha_k = \frac{r_k^T r_k}{p_k^T A p_k}$$
• ②更新解： $$x_{k+1} = x_k + \alpha_k p_k$$
• ③计算新的残差： $$r_{k+1} = r_k - \alpha_k A p_k$$
• ④判断收敛条件：如果$\|r_{k+1}\| < \epsilon$，则停止迭代，输出最优解。
• ⑤计算新的搜索方向$\beta_k$： $$\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$$

$$p_{k+1} = r_{k+1} + \beta_k p_k$$

• ⑥重复迭代：返回步骤②，直到满足停止条件或达到最大迭代次数。

以一个具体的二次函数为例：

$$f(x_1, x_2) = 3x_1^2 + 2x_1x_2 + x_2^2 - 4x_1 - 2x_2$$

可以将其形式转换为矩阵表达：

$$A = \begin{pmatrix} 6 & 2 \\ 2 & 2 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 2 \end{pmatrix}$$

初始点设定为$x_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$。

import numpy as np

# 定义Hessian矩阵和常数向量
A = np.array([[6, 2], [2, 2]])
b = np.array([4, 2])

# 定义目标函数的梯度
def grad_f(x):
    return A @ x - b  # 梯度为 Ax - b

# 解析求解最优解
def solve_optimal_solution(A, b):
    # 解线性方程 Ax = b
    x_optimal = np.linalg.solve(A, b)
    return x_optimal

# 计算目标函数
def objective_function(x):
    return 3*x[0]**2 + 2*x[0]*x[1] + x[1]**2 - 4*x[0] - 2*x[1]

# 调用求解函数
optimal_x = solve_optimal_solution(A, b)
print(f"最优解：(x_1, x_2) = ({optimal_x[0]:.2f}, {optimal_x[1]:.2f})")

# 计算最优解的目标函数值
optimal_value = objective_function(optimal_x)
print(f"最优解的目标函数值 = {optimal_value:.2f}")

五、变尺度法

变尺度法（Variable Metric Method）是一种优化算法，旨在通过在每次迭代中更新一个可变的尺度（或度量）来改进最优解的搜索过程。该方法尤其适合无约束优化问题，并通过引入二次型近似来提高收敛速度。

算法背景
变尺度法利用目标函数的梯度信息，构造一个二次型模型以近似目标函数，并使用该模型来指导搜索过程。与牛顿法不同，变尺度法不需要计算Hessian矩阵的逆，而是利用梯度的线性组合来调整搜索方向。

算法步骤

• ①初始化：设定初始点 $x_0$，初始搜索方向 $p_0$，以及初始步长 $\alpha_0$。
• ②迭代过程：
  • 计算当前点的梯度 $g_k = \nabla f(x_k)$。
  • 计算步长 $ \alpha_k$：$$\alpha_k = \frac{-g_k^T p_k}{p_k^T H_k p_k}$$其中$ H_k$ 为可变的近似Hessian矩阵。
  • 更新解： $$x_{k+1} = x_k + \alpha_k p_k$$
  • 更新梯度： $$g_{k+1} = \nabla f(x_{k+1})$$
  • 更新搜索方向 $ p_{k+1}$：$$p_{k+1} = -g_{k+1} + \beta_k p_k$$其中$ \beta_k$ 是通过某种方式计算得到的（如Polak-Ribiere方法）。
  • 判断收敛条件：如果 $|g_{k+1}| < \epsilon$，则停止迭代。
• ③重复迭代：返回步骤②，直到满足停止条件或达到最大迭代次数。

算例示范
以一个具体的目标函数为例： $$f(x_1, x_2) = (x_1 - 1)^2 + (x_2 - 2)^2 + (x_1 - x_2)^2$$
初始点设定为 $x_0 = \begin{pmatrix} 0 \ 0 \end{pmatrix}$。

import numpy as np
from scipy.optimize import minimize

# 定义新的目标函数
def objective_function(x):
    return (x[0] - 1)**2 + (x[1] - 2)**2 + (x[0] - x[1])**2

# 初始点
x0 = np.array([0.0, 0.0])

# 调用优化函数
result = minimize(objective_function, x0)

# 输出结果，保留两位小数
print("最优解: [ {:.2f}, {:.2f} ]".format(result.x[0], result.x[1]))
print("目标函数值: {:.2f}".format(result.fun))

总结

无约束最优化问题作为优化理论中的一个重要分支，拥有深厚的理论基础和广泛的实际应用。最速下降法、牛顿法、共轭梯度法和变尺度法等算法各具特色，适用于不同类型的优化问题。随着计算技术的发展和优化算法的不断进步，解决大规模、复杂无约束最优化问题的方法将会越来越多样化，应用领域也将更加广泛。未来的研究可能会更加注重算法的全局收敛性、计算效率以及应用的广泛性，为实际问题的优化求解提供更加有效的工具。约束最优化问题在实际应用中有着广泛的用途，特别是在机器学习、信号处理、工程优化、经济预测等领域。例如，在机器学习中的模型训练过程中，常常需要通过无约束最优化算法来最小化损失函数。随着计算能力的提升和算法的改进，无约束最优化算法在解决更大规模、更复杂的应用问题时展现出了巨大的潜力。近年来，智能优化算法如遗传算法、粒子群算法等逐渐崭露头角，成为无约束最优化问题求解的有力补充。尽管这些算法不依赖于目标函数的梯度信息，但它们在面对高度非线性、不可微和多峰函数时表现出了较强的鲁棒性和全局搜索能力。然而，这类智能算法的收敛速度通常不如传统的梯度类算法，因此在实际应用中常常需要与经典的无约束最优化算法相结合。

参考资料

1. 无约束最优化方法
2. Scipy优化使用教程】一、Scipy中无约束优化的六大算法