### 最优化理论与算法——梯度法与阻尼牛顿法
#### 一、梯度法
**1、算法原理**
梯度法是一种基于一阶导数（或梯度）来寻找函数最小值的方法。其基本思想是沿着负梯度的方向进行搜索，因为负梯度方向是指函数下降最快的方向。
**2、算法步骤**
- **初始化**: 选择一个初始点 \( x_0 \) 和一个允许误差 \( \epsilon > 0 \)。
- **计算梯度**: 在当前点 \( x_k \) 处计算梯度 \( g_k = \nabla f(x_k) \)。
- **确定搜索方向**: 搜索方向通常定义为负梯度方向，即 \( d_k = -g_k \)。
- **线搜索**: 使用某种线搜索方法确定步长 \( \alpha_k \)，使得 \( x_{k+1} = x_k + \alpha_k d_k \)。
- **停止条件**: 如果 \( \| g_k \| < \epsilon \)，则停止迭代，返回 \( x_k \) 作为近似解；否则，令 \( k = k + 1 \)，并重复步骤2至4。
**3、代码实现**
下面给出了梯度法的一个简单MATLAB实现：
```matlab
% 定义目标函数 fun
% 定义梯度计算函数 gfun
% 梯度法主程序
function [k, x, val] = grad(fun, gfun, x0, epsilon)
maxk = 5000; % 最大迭代次数
beta = 0.5; sigma = 0.4; % 参数设置
k = 0;
while (k < maxk)
gk = feval(gfun, x0); % 计算梯度
dk = -gk; % 计算搜索方向
if (norm(gk) < epsilon) % 终止准则
break;
end
m = 0; mk = 0;
while (m < 20) % Armijo搜索求步长
if (feval(fun, x0 + beta^m * dk) <= feval(fun, x0) + sigma * beta^m * gk' * dk)
mk = m; break;
end
m = m + 1;
end
x0 = x0 + beta^mk * dk; % 更新点
k = k + 1;
end
x = x0; % 近似最优点
val = feval(fun, x0); % 最优值
end
```
**4、结果**
假设我们采用上述MATLAB代码，并使用一个具体的测试函数，例如Rosenbrock函数，我们可以观察到梯度法的收敛过程。
**5、分析**
梯度法的优点在于简单且易于实现。然而，它的收敛速度可能比较慢，尤其是在接近极小点时，由于梯度变得非常小。此外，如果目标函数具有鞍点，则梯度法可能会被这些点所吸引，导致不收敛或收敛速度非常慢。
#### 二、阻尼牛顿法
**1、算法原理**
阻尼牛顿法是在牛顿法的基础上加入阻尼因子的一种改进算法。牛顿法利用了目标函数的二阶导数信息（即Hessian矩阵），以期望更快地收敛到极小点。但在某些情况下，Hessian矩阵可能不是正定的，这会导致迭代不稳定甚至发散。为了克服这一问题，阻尼牛顿法通过引入步长因子 \( \alpha_k \) 来控制每次迭代的步长大小。
**2、算法步骤**
- **初始化**: 选择一个初始点 \( x_0 \) 和允许误差 \( \epsilon > 0 \)。
- **计算梯度和Hessian**: 在当前点 \( x_k \) 处计算梯度 \( g_k = \nabla f(x_k) \) 和Hessian矩阵 \( H_k = \nabla^2 f(x_k) \)。
- **确定搜索方向**: 搜索方向定义为 \( d_k = -H_k^{-1}g_k \)。
- **线搜索**: 使用某种线搜索方法确定步长 \( \alpha_k \)，使得 \( x_{k+1} = x_k + \alpha_k d_k \)。
- **停止条件**: 如果 \( \| g_k \| < \epsilon \)，则停止迭代，返回 \( x_k \) 作为近似解；否则，令 \( k = k + 1 \)，并重复步骤2至4。
**3、代码实现**
阻尼牛顿法的关键在于如何计算Hessian矩阵以及它的逆。下面是一个简单的MATLAB实现示例：
```matlab
% 理查德外推法求海瑟矩阵
function He = compute_Hess(fun, x0, h)
n = length(x0);
He = zeros(n, n);
temp = x0;
for i = 1:n
temp(i) = temp(i) + 0.5*h;
He(i, i) = 16 * feval(fun, temp);
temp(i) = temp(i) - h;
He(i, i) = He(i, i) + 16 * feval(fun, temp);
temp(i) = temp(i) + 3*h/2;
He(i, i) = He(i, i) - feval(fun, temp);
temp(i) = temp(i) - 2*h;
He(i, i) = He(i, i) - feval(fun, temp);
He(i, i) = He(i, i) - 30 * feval(fun, x0);
He(i, i) = He(i, i) / (3*h*h);
temp = x0;
end
% 以下代码用于计算非对角线元素
for i = 1:n
for j = i+1:n
temp(i) = temp(i) + h/2;
temp(j) = temp(j) + h/2;
He(i, j) = 4 * feval(fun, temp);
temp(j) = temp(j) - h;
He(i, j) = He(i, j) - 4 * feval(fun, temp);
temp(j) = temp(j) + h;
temp(i) = temp(i) - h;
He(i, j) = He(i, j) - feval(fun, temp);
temp(i) = temp(i) + h;
He(i, j) = He(i, j) / (h*h);
He(j, i) = He(i, j); % 因为Hessian是对称的
temp = x0;
end
end
end
```
这段代码实现了海瑟矩阵的数值估计。接下来，可以结合线搜索方法来确定步长，并完成阻尼牛顿法的完整实现。
**总结**
- **梯度法** 是一种基于梯度下降方向的简单迭代方法，适用于大规模问题，但由于仅使用了一阶导数信息，其收敛速度可能较慢。
- **阻尼牛顿法** 利用了二阶导数信息来加速收敛，但在实际应用中需要注意Hessian矩阵的正定性问题，并通过引入阻尼因子来确保算法的稳定性。
通过以上分析可以看出，两种方法各有优势和局限性，选择哪种方法取决于具体的应用场景和目标函数的特性。
