凸优化

1 纯凸优化问题、非凸局部优化问题的解法

Created 2021.03.22 by Cong Yu; Last modified: 2024.06.22-v4.3.2

Copyleft! 2024 Cong Yu. Some rights reserved.

1 纯凸优化问题、非凸局部优化问题的解法

refitem: METHODS-FOR-NON-LINEAR-LEAST-SQUARES-PROBLEMS

有代码可查

讨论一类朴素的无约束凸优化问题的解法:

\[\begin{align} & \textrm{Given }F：\mathbb{R} ^{n} \mapsto \mathbb{R} \\ & \textrm{Find } \mathbf{x}^+ = \text {argmin}_\mathbf{x} F(\mathbf x) \\ \end{align}\]

且其中目标函数 $F$ 是凸函数

或者讨论非凸局部优化问题的解法

1.0 解析法

$F(\mathbf x)$ 是可微函数，则最优点 $\mathbf x^+$ 满足充要条件: $\nabla F(\mathbf x^+) = 0$

由该充要条件可得到 ==解析法== 求解上面的问题，比较简单：

求解 $\nabla F(\mathbf x) \stackrel{let}= 0$，解得$x^*$ 即为马鞍点，或者极值点。

但是往往会出现$\nabla F(\mathbf x)$ 比较容易写出，但是$\nabla F(\mathbf x) = 0$ 却无法求解的问题，或者干脆$\nabla F(\mathbf x)$ 也无法写出的问题。故无法使用解析解法，可以采用数值解法，即优化迭代方法。

1.1 下降方法

下降方法是一类基本的优化算法，其核心思想是通过迭代逐步减小目标函数值。主要包括：

1.1.1 精确直线搜索

在每次迭代中，沿着搜索方向找到使目标函数最小的步长：

\[\alpha_k = \text{argmin}_{\alpha > 0} f(x_k + \alpha d_k)\]

其中 $d_k$ 是搜索方向。精确直线搜索虽然能保证每次迭代都取得最大可能的下降，但计算成本较高。

1.1.2 回溯直线搜索

为了避免精确搜索的高计算成本，使用回溯策略：

选择初始步长 $\alpha_0$
如果 $f(x_k + \alpha d_k) > f(x_k) + c_1\alpha \nabla f(x_k)^T d_k$，则减小步长
重复步骤2直到满足条件

其中 $c_1 \in (0,1)$ 是控制参数。回溯搜索通过逐步减小步长来寻找满足条件的点，计算效率更高。

1.1.3 最速下降法

最速下降法选择负梯度方向作为搜索方向：

\[d_k = -\nabla f(x_k)\]

这是最直观的下降方向，因为梯度方向是函数值增长最快的方向，其反方向就是下降最快的方向。

1.2 梯度下降法

梯度下降法是最基本的优化算法之一：

一阶方法：只使用目标函数的一阶导数信息
基本思路：一阶泰勒展开，每次朝梯度下降的方向迭代更新x
优点：实现简单，计算量小
缺点：收敛速度可能较慢，对条件数敏感

原理推导：

\[F(\mathbf x+\Delta \mathbf x) \approx F(\mathbf x) + J \Delta \mathbf x\]

令 $F(\mathbf x+\Delta\mathbf x)< F(\mathbf x)$ 即可，即 $J\Delta\mathbf x < 0$，即 $\Delta \mathbf x$ 与 $J$ 向量的夹角大于90°即可。在这个范围内，下降速度最快的自然就是直接取$J$的反方向：即当 $\Delta \mathbf x=-J^T$ 时，下降最快。

但是如果直接按照梯度更新，步子迈的太大，所以引入参数学习率 $\eta$，通常小于1。

于是最终的迭代公式即

\[\mathbf x:= \mathbf x-\eta J^T\]

学习率的选择对算法性能有重要影响：

太大：可能导致震荡或不收敛
太小：收敛速度慢
自适应：可以根据迭代过程动态调整

1.3 梯度下降衍生算法

普通的梯度下降方法选择的是梯度的反方向。实际计算过程中，由于各参数的尺度不同，可能出现各种收敛缓慢的情况，而且普通梯度下降法中的步长参数$\eta$不可调，可以考虑做自适应的算法。由此衍生出一系列算法:

(注：本小节$\theta$ 即前文所述被优化的变量 $\mathbf{x}$)

refitem:

1.3.1 Momentum 动量法

动量法通过引入动量项来加速收敛，类似于物理中的动量概念：

\[\begin{align} \upsilon_t &= \gamma \upsilon_{t−1} + \eta \nabla_{\theta}J(\theta) \\ \theta &:= \theta− \upsilon_t \end{align}\]

$\gamma$ 通常设置为0.9或更小，用于控制历史梯度的影响程度。动量法的主要优点：

加速收敛
减少震荡
帮助跳出局部最小值

1.3.2 Adagrad

Adagrad为每个参数提供自适应学习率，根据历史梯度信息动态调整：

\[θ_{t+1,i}=θ_{t,i}−\frac η {\sqrt {G_{t,ii}+ϵ}}∇_θJ(θ_{t,i})\]

其中 $G_{t,ii}$ 是参数 $\theta_i$ 的历史梯度平方和。特点：

自动调整学习率
适合处理稀疏数据
学习率会随时间单调递减

1.3.3 RMSProp

RMSProp（Root Mean Square Propagation）是 Adagrad 的改进版本，解决了 Adagrad 学习率单调递减的问题：

\[\begin{align} E[g^2]_t &= \gamma E[g^2]_{t-1} + (1-\gamma)g_t^2 \\ \theta_{t+1} &= \theta_t - \frac{\eta}{\sqrt{E[g^2]_t + \epsilon}} g_t \end{align}\]

其中 $\gamma$ 通常设置为 0.9，$\epsilon$ 是一个很小的数（如 1e-8）用于防止除零。优点：

学习率不会单调递减
对非平稳目标函数表现更好
适合在线学习

1.3.4 AdaDelta

AdaDelta 是 RMSProp 的扩展，不需要设置初始学习率：

\[\begin{align} E[g^2]_t &= \gamma E[g^2]_{t-1} + (1-\gamma)g_t^2 \\ \Delta\theta_t &= -\frac{\sqrt{E[\Delta\theta^2]_{t-1} + \epsilon}}{\sqrt{E[g^2]_t + \epsilon}} g_t \\ E[\Delta\theta^2]_t &= \gamma E[\Delta\theta^2]_{t-1} + (1-\gamma)\Delta\theta_t^2 \\ \theta_{t+1} &= \theta_t + \Delta\theta_t \end{align}\]

特点：

不需要手动设置学习率
对超参数不敏感
计算效率高

1.3.5 Adam

Adam（Adaptive Moment Estimation）结合了动量法和 RMSProp 的优点：

\[\begin{align} m_t &= \beta_1 m_{t-1} + (1-\beta_1)g_t \\ v_t &= \beta_2 v_{t-1} + (1-\beta_2)g_t^2 \\ \hat{m}_t &= \frac{m_t}{1-\beta_1^t} \\ \hat{v}_t &= \frac{v_t}{1-\beta_2^t} \\ \theta_{t+1} &= \theta_t - \frac{\eta}{\sqrt{\hat{v}_t} + \epsilon} \hat{m}_t \end{align}\]

其中 $\beta_1$ 和 $\beta_2$ 通常分别设置为 0.9 和 0.999。优点：

结合了动量和自适应学习率的优点
对超参数不敏感
计算效率高
适合大规模数据和参数

1.3.6 AdaMax

AdaMax 是 Adam 的变体，使用无穷范数来更新参数：

\[\begin{align} m_t &= \beta_1 m_{t-1} + (1-\beta_1)g_t \\ u_t &= \max(\beta_2 u_{t-1}, |g_t|) \\ \theta_{t+1} &= \theta_t - \frac{\eta}{u_t} m_t \end{align}\]

其中 $u_t$ 依赖于参数梯度的无穷范数，这使其在某些情况下比 Adam 更稳定。特点：

对梯度范围更敏感
在某些问题上表现更好
计算量略小于Adam

所有方法的对比：

pytorch

https://pytorch.org/docs/stable/optim.html#

pytorch中提供的优化方法

1.4 牛顿法

牛顿法是一种二阶优化方法，利用目标函数的二阶导数信息：

二阶方法：使用目标函数的一阶和二阶导数信息
二阶导矩阵 $H$ 计算量比较大
基本思路：2阶泰勒展开，使用2次曲线拟合；然后求近似后函数的导数，令导数等于0，即求出了当前近似二次函数的最小值；然后在此处继续进行2阶泰勒展开，重复上述步骤。

优点：

收敛速度快（二次收敛）
对条件数不敏感
不需要手动设置学习率

缺点：

计算Hessian矩阵及其逆的计算量大
需要目标函数二阶可导
对初始点敏感

原理推导：

目标函数展开：

对于当前参数$\mathbf x$，我们希望找到增量 $\Delta \mathbf x$ 使得 $F(\mathbf x + \Delta \mathbf x)$ 最小。将$F$在$\mathbf x$位置二阶泰勒展开 $F(\mathbf x+\Delta \mathbf x) \approx F(\mathbf x) + J \Delta \mathbf x + \frac12\Delta\mathbf x^T H \Delta\mathbf x \equiv G$ 其中 $J$ 是雅可比矩阵，$H$ 是海森矩阵。

增量方程：将二阶展开后的近似式子定义为G，G对$\Delta \mathbf x$求导，得

\[\frac {\partial G }{\partial \Delta \mathbf x} = J^T + H \Delta \mathbf x\]

令其=0，得增量方程

\[H \Delta x = -J^T\]

解得

\[\Delta \mathbf x = - H^{-1} J^T\]

得出迭代公式：

\[\mathbf x:= \mathbf x-\eta H^{-1}J^T\]

1.5 混合方法

混合方法根据Hessian矩阵的性质采用不同的策略：

当Hessian矩阵正定时：
- 使用完整的牛顿步长
- 利用二次收敛性质
- 适用于局部凸区域
当Hessian矩阵接近奇异时：
- 使用梯度下降步长
- 避免数值不稳定性
- 适用于非凸区域
当Hessian矩阵负定时：
- 使用负梯度方向
- 避免向最大值方向移动
- 适用于局部凹区域

混合方法的优点：

自适应性强
数值稳定性好
收敛性能好

1.6 阻尼方法

阻尼方法通过引入正则化项来改善牛顿法的数值稳定性：

\[\Delta x = -(H + \lambda I)^{-1}J^T\]

其中 $\lambda$ 是阻尼参数，$I$ 是单位矩阵。特点：

改善Hessian矩阵的条件数
控制步长大小
提高算法鲁棒性

阻尼参数的选择策略：

自适应调整
基于信赖域方法
基于线搜索方法

1.7 L-BFGS

L-BFGS（Limited-memory BFGS）是一种拟牛顿法，适用于大规模优化问题：

基本思想：
- 近似Hessian矩阵的逆
- 只存储最近的m个向量对
- 通过递归公式计算搜索方向
数学推导：

2.1 BFGS更新公式

首先，BFGS算法的Hessian矩阵更新公式为：
\[H_{k+1} = (I - \rho_k s_k y_k^T)H_k(I - \rho_k y_k s_k^T) + \rho_k s_k s_k^T\]
其中：
- $s_k = x_{k+1} - x_k$ 是参数更新向量
- $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$ 是梯度更新向量
- $\rho_k = \frac{1}{y_k^T s_k}$
2.2 L-BFGS的递归公式

L-BFGS通过递归方式计算搜索方向，避免存储完整的Hessian矩阵。对于给定的初始Hessian近似$H_0$，搜索方向$d_k$可以通过以下递归公式计算：
\[d_k = -H_k \nabla f(x_k)\]
其中$H_k$可以通过以下递归公式计算：
\[H_k = (V_{k-1}^T \cdots V_{k-m}^T)H_0(V_{k-m} \cdots V_{k-1}) + \rho_{k-m}(V_{k-1}^T \cdots V_{k-m+1}^T)s_{k-m}s_{k-m}^T(V_{k-m+1} \cdots V_{k-1}) + \rho_{k-m+1}(V_{k-1}^T \cdots V_{k-m+2}^T)s_{k-m+1}s_{k-m+1}^T(V_{k-m+2} \cdots V_{k-1}) + \cdots + \rho_{k-1}s_{k-1}s_{k-1}^T\]
其中：
- $V_k = I - \rho_k y_k s_k^T$
- $H_0$通常取为单位矩阵的缩放版本
2.3 双循环递归算法

L-BFGS使用双循环递归算法来高效计算搜索方向：
1. 内循环：计算$q = \nabla f(x_k)$
2. 外循环：通过递归公式计算$H_k \nabla f(x_k)$
算法步骤：
- 初始化搜索方向
- 更新历史信息
- 计算新的搜索方向
- 执行线搜索
优点：
- 内存需求低
- 计算效率高
- 适合大规模问题
实现细节：
- 通常m取5-20
- 需要定期重启
- 可以使用不同的线搜索策略
优化策略：

a. 初始Hessian近似：
- 通常使用$H_0 = \gamma I$，其中$\gamma = \frac{y_k^T s_k}{y_k^T y_k}$
- 这个选择可以改善算法的收敛性
b. 线搜索策略：
- 使用Wolfe条件确保充分下降
- 可以使用回溯线搜索提高效率
c. 内存管理：
- 定期重启算法
- 动态调整内存大小m
d. 数值稳定性：
- 使用缩放技术
- 处理病态问题
优缺点分析：

优点：
- 内存需求低：只需要存储最近的m个向量对
- 计算效率高：通过递归公式避免存储完整Hessian矩阵
- 收敛性能好：在大多数问题上表现出超线性收敛
- 适合大规模问题：特别适合高维优化问题
缺点：
- 对初始点敏感
- 需要定期重启
- 在非凸问题上可能陷入局部最小值

参考资源：

https://aria42.com/blog/2014/12/understanding-lbfgs
https://www.chokkan.org/software/liblbfgs/
https://github.com/ZJU-FAST-Lab/LBFGS-Lite/

1.9 代码实现

下面是一个简单的Python实现示例，展示了不同优化方法的使用：

import numpy as np
from scipy.optimize import minimize

def gradient_descent(f, grad_f, x0, learning_rate=0.01, max_iter=1000, tol=1e-6):
    """
    梯度下降法实现
    
    参数:
    f: 目标函数
    grad_f: 梯度函数
    x0: 初始点
    learning_rate: 学习率
    max_iter: 最大迭代次数
    tol: 收敛容差
    
    返回:
    x: 最优解
    history: 优化历史
    """
    x = x0.copy()
    history = []
    
    for i in range(max_iter):
        grad = grad_f(x)
        if np.linalg.norm(grad) < tol:
            break
            
        x = x - learning_rate * grad
        history.append(f(x))
        
    return x, history

def newton_method(f, grad_f, hess_f, x0, max_iter=1000, tol=1e-6):
    """
    牛顿法实现
    
    参数:
    f: 目标函数
    grad_f: 梯度函数
    hess_f: Hessian矩阵函数
    x0: 初始点
    max_iter: 最大迭代次数
    tol: 收敛容差
    
    返回:
    x: 最优解
    history: 优化历史
    """
    x = x0.copy()
    history = []
    
    for i in range(max_iter):
        grad = grad_f(x)
        if np.linalg.norm(grad) < tol:
            break
            
        hess = hess_f(x)
        try:
            delta = np.linalg.solve(hess, -grad)
        except np.linalg.LinAlgError:
            # 如果Hessian矩阵奇异，使用伪逆
            delta = -np.linalg.pinv(hess) @ grad
            
        x = x + delta
        history.append(f(x))
        
    return x, history

# 示例：优化Rosenbrock函数
def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

def rosenbrock_grad(x):
    return np.array([
        -2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2),
        200 * (x[1] - x[0]**2)
    ])

def rosenbrock_hess(x):
    return np.array([
        [2 - 400 * x[1] + 1200 * x[0]**2, -400 * x[0]],
        [-400 * x[0], 200]
    ])

# 使用示例
x0 = np.array([-1.5, 1.5])
x_gd, hist_gd = gradient_descent(rosenbrock, rosenbrock_grad, x0)
x_newton, hist_newton = newton_method(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0)

详细解读凸优化问题，以及凸优化求解方法，如：梯度下降，牛顿法，混合方法，LBFGS方法等