损失景观临界点分析：Hessian谱与曲率动力学

深度神经网络训练过程中，优化器在复杂的损失景观中导航，其轨迹会受到**临界点（Critical Points）**的强烈影响。临界点是梯度为零的点，包括局部最小值、鞍点和局部最大值。本文档系统分析神经网络损失景观中临界点的结构、分类方法及其在训练过程中的演化规律。

临界点基础理论

定义

临界点（Critical Point）：满足 $\nabla f({\theta }^{\ast })=0$ 的参数点 ${\theta }^{\ast }$。

分类

设 $H({\theta }^{\ast })={\nabla }^{2}f({\theta }^{\ast })$ 是损失函数在 ${\theta }^{\ast }$ 处的Hessian矩阵。

类型	Hessian特征值	Index	几何特征
局部最小值	全部为正	0	局部”碗”形
鞍点	有正有负	$k>0$	马鞍面
局部最大值	全部为负	$n$	局部”倒碗”形

Index定义：$i({\theta }^{\ast })=|\{{\lambda }_i(H):{\lambda }_i<0\}|$，即负特征值的个数。

High-dimensional视角

在高维神经网络中（参数可达数十亿），临界点的分类具有独特性质：

• 维度诅咒与祝福：在高维空间中，几乎所有critical points都是鞍点
• 局部最小值的稀有性：真正的局部最小值在总critical points中占比极小
• 鞍点的普遍性：这使得优化算法需要能够逃离鞍点

import numpy as np
import torch
 
def classify_critical_point(hessian_eigenvalues, tol=1e-6):
    """
    根据Hessian特征值分类临界点
    
    Args:
        hessian_eigenvalues: Hessian矩阵的特征值
        
    Returns:
        str: 临界点类型
        int: index (负特征值数量)
    """
    n_negative = np.sum(hessian_eigenvalues < -tol)
    n_positive = np.sum(hessian_eigenvalues > tol)
    n_zero = len(hessian_eigenvalues) - n_negative - n_positive
    
    if n_negative == 0 and n_positive > 0:
        return "local_minimum", 0
    elif n_negative > 0 and n_positive > 0:
        return "saddle_point", n_negative
    elif n_negative == len(hessian_eigenvalues):
        return "local_maximum", len(hessian_eigenvalues)
    else:
        return "degenerate_critical_point", n_negative

Hessian谱分析

谱分布理论

深度神经网络损失景观的Hessian谱呈现独特的结构：

1. 批量归一化网络的大特征值

对于带BatchNorm的网络，Hessian的最大特征值通常与以下因素相关：

• Batch size
• 学习率
• 网络深度

2. 特征值尺度分离

$${\lambda }_1\gg {\lambda }_2\ge \dots \ge {\lambda }_n\gg {\lambda }_{\min }$$

Hessian计算方法

精确Hessian（小型网络）

def compute_exact_hessian(model, dataloader, criterion):
    """
    计算精确Hessian（适用于小网络）
    
    H[i,j] = ∂²L/∂θᵢ∂θⱼ
    """
    model.eval()
    n_params = sum(p.numel() for p in model.parameters())
    hessian = torch.zeros(n_params, n_params)
    
    for inputs, targets in dataloader:
        model.zero_grad()
        outputs = model(inputs)
        loss = criterion(outputs, targets)
        loss.backward()
        
        # 收集梯度
        gradients = torch.cat([p.grad.flatten() for p in model.parameters()])
        
        # 计算二阶导数
        for i in range(n_params):
            grad_vec = gradients.clone()
            # 使用Hessian-vector乘积
            hessian[i] = compute_hessian_vector_product(model, grad_vec)
    
    return hessian
 
def compute_hessian_vector_product(model, vec):
    """计算Hessian-向量乘积 H*v"""
    params = [p for p in model.parameters() if p.gradient is not None]
    result = torch.zeros_like(vec)
    idx = 0
    
    for p in params:
        numel = p.numel()
        # 手动设置梯度向量
        p.grad = vec[idx:idx+numel].view(p.shape).clone()
        # 清理之前的梯度
        model.zero_grad()
        
        # 反向传播得到H*v
        if p.grad is not None:
            result[idx:idx+numel] = p.grad.view(-1)
        idx += numel
    
    return result

本征正交分量（EOC）方法

对于大型网络，使用**本征正交分量（Eigen Orthogonal Components, EOC）**近似Hessian谱¹：

class EigenOrthogonalComponents:
    """
    EOC方法：高效的Hessian谱分析
    
    不需要显式计算整个Hessian矩阵
    """
    def __init__(self, model, data_loader, criterion, n_components=100):
        self.model = model
        self.n_components = n_components
        self.eigenvalues = None
        self.eigenvectors = None
        
        # 计算Hessian的power method迭代
        self._compute_spectrum(data_loader, criterion)
    
    def _compute_spectrum(self, data_loader, criterion):
        """使用随机幂迭代计算特征值分布"""
        n_params = sum(p.numel() for p in self.model.parameters())
        
        # 初始化随机向量
        v = torch.randn(n_params, device=next(self.model.parameters()).device)
        v = v / v.norm()
        
        # Power method迭代
        for _ in range(self.n_components):
            # 计算Hv
            Hv = self._hessian_vector_product(v, data_loader, criterion)
            
            # 正交化
            for j in range(len(self.eigenvectors)):
                Hv = Hv - torch.dot(Hv, self.eigenvectors[j]) * self.eigenvectors[j]
            
            v = Hv / Hv.norm()
            
            # 记录特征值估计
            eigenvalue = torch.dot(v, Hv).item()
            self.eigenvalues.append(eigenvalue)
            self.eigenvectors.append(v.clone())
    
    def _hessian_vector_product(self, vec, data_loader, criterion):
        """计算Hessian-向量乘积"""
        model = self.model
        model.zero_grad()
        
        for inputs, targets in data_loader:
            outputs = model(inputs)
            loss = criterion(outputs, targets)
            # 使用自动微分计算H*v
            grads = torch.autograd.grad(
                loss, model.parameters(), 
                create_graph=True
            )
            
            grad_vec = torch.cat([g.view(-1) for g in grads])
            grad_grad = torch.autograd.grad(
                grad_vec.dot(vec), model.parameters()
            )
            
            Hv = torch.cat([g.view(-1) for g in grad_grad])
        
        return Hv
    
    def plot_spectrum(self):
        """绘制特征值分布"""
        import matplotlib.pyplot as plt
        
        plt.figure(figsize=(10, 6))
        plt.plot(range(len(self.eigenvalues)), sorted(self.eigenvalues, reverse=True))
        plt.xlabel('Index')
        plt.ylabel('Eigenvalue')
        plt.title('Hessian Eigenvalue Spectrum')
        plt.yscale('symlog' if any(e < 0 for e in self.eigenvalues) else 'log')
        plt.grid(True)
        plt.show()

谱的典型结构

训练好的神经网络Hessian谱呈现三区域结构：

特征值密度
    │
    │        ████
    │       ██████
    │      ████████
    │     █████████
    │    ███████████
    │   ████████████
    │  █████████████
    │ ██████████████
    │████████████████
────┴──────────────────────────────────────→ 特征值
    │              │              │
  负特征值        零附近        正特征值
  (鞍点方向)    (平坦方向)    (最小值方向)

def analyze_spectrum_regions(eigenvalues, tol_small=1e-3, tol_large=1.0):
    """
    分析Hessian谱的三个区域
    
    Returns:
        dict: 各区域的统计信息
    """
    eigenvalues = np.array(eigenvalues)
    
    negative_mask = eigenvalues < -tol_large
    small_mask = (eigenvalues >= -tol_large) & (eigenvalues < tol_large)
    positive_mask = eigenvalues >= tol_large
    
    return {
        'n_negative': np.sum(negative_mask),
        'n_small': np.sum(small_mask),
        'n_positive': np.sum(positive_mask),
        'fraction_negative': np.sum(negative_mask) / len(eigenvalues),
        'fraction_small': np.sum(small_mask) / len(eigenvalues),
        'fraction_positive': np.sum(positive_mask) / len(eigenvalues),
        'max_eigenvalue': np.max(eigenvalues),
        'min_eigenvalue': np.min(eigenvalues),
        'mean_eigenvalue': np.mean(eigenvalues),
        'condition_number': np.max(eigenvalues) / max(np.min(eigenvalues), 1e-10)
    }

训练过程中的曲率演化

Critical Learning Rate

关键发现：存在一个临界学习率（Critical Learning Rate），超过这个值训练会变得不稳定²。

理论推导：

对于梯度下降更新 ${\theta }_{t+1}={\theta }_t-\eta \nabla f({\theta }_t)$，考虑在局部最小值附近的稳定性：

在Hessian为 $H$ 的点，梯度下降的更新可以近似为：

$${\theta }_{t+1}\approx (I-\eta H){\theta }_t$$

稳定性条件：$|1-\eta {\lambda }_i|<1$ for all $i$

这给出：

$$0<\eta <\frac{2}{{\lambda }_{\max }(H)}$$

临界学习率：

$${\eta }_{\text{crit}}=\frac{2}{{\lambda }_{\max }}$$

def compute_critical_learning_rate(hessian_eigenvalues):
    """
    计算临界学习率
    
    η_crit = 2 / λ_max
    """
    lambda_max = np.max(hessian_eigenvalues)
    return 2.0 / lambda_max
 
def analyze_lr_stability(lr, hessian_eigenvalues):
    """
    分析给定学习率的稳定性
    """
    stable_mask = np.abs(1 - lr * hessian_eigenvalues) < 1
    unstable_mask = ~stable_mask
    
    return {
        'stable_fraction': np.mean(stable_mask),
        'unstable_fraction': np.mean(unstable_mask),
        'max_spectral_radius': np.max(np.abs(1 - lr * hessian_eigenvalues)),
        'is_stable': np.all(unstable_mask == 0)
    }

Edge of Stability现象

实验发现：当使用接近或超过 ${\eta }_{\text{crit}}$ 的学习率时，会出现Edge of Stability现象：

1. 瞬态不稳定期：训练初期loss可能上升
2. 自适应稳定化：网络自动调整使得 $\eta {\lambda }_{\max }\approx 2$
3. 长期稳定：最终进入稳定状态

class EdgeOfStabilityTracker:
    """
    追踪Edge of Stability现象
    """
    def __init__(self):
        self.lr_history = []
        self.loss_history = []
        self.sharpness_history = []  # λ_max的追踪
        
    def update(self, lr, loss, hessian_sharpness):
        self.lr_history.append(lr)
        self.loss_history.append(loss)
        self.sharpness_history.append(hessian_sharpness)
    
    def compute_effective_lr(self):
        """
        计算有效学习率与sharpness的乘积
        """
        return [lr * sharp for lr, sharp in 
                zip(self.lr_history, self.sharpness_history)]
    
    def plot_eos_dynamics(self):
        """可视化Edge of Stability动态"""
        import matplotlib.pyplot as plt
        
        fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 10))
        
        steps = range(len(self.loss_history))
        
        # Loss曲线
        ax1.plot(steps, self.loss_history)
        ax1.set_ylabel('Loss')
        ax1.set_title('Training Loss')
        ax1.grid(True)
        
        # Sharpness曲线
        ax2.plot(steps, self.sharpness_history)
        ax2.set_ylabel('λ_max (Sharpness)')
        ax2.set_title('Hessian Sharpness')
        ax2.grid(True)
        
        # 有效学习率
        effective_lr = self.compute_effective_lr()
        ax3.plot(steps, effective_lr, label='η * λ_max')
        ax3.axhline(y=2.0, color='r', linestyle='--', label='Stability boundary')
        ax3.set_xlabel('Training Step')
        ax3.set_ylabel('η * λ_max')
        ax3.set_title('Edge of Stability')
        ax3.legend()
        ax3.grid(True)
        
        plt.tight_layout()
        plt.show()

训练曲线的曲率演化

典型的训练过程中，Hessian谱的演化规律：

训练阶段	Loss水平	λ_max范围	曲率特征
初始阶段	高	大	负曲率方向多
快速下降期	快速下降	先增后减	接近临界
收敛期	低	稳定	主要正曲率

def simulate_curvature_evolution(n_steps=1000, lr=0.01):
    """
    模拟训练过程中的曲率演化
    """
    sharpness = []
    loss = []
    
    # 初始高sharpness
    current_sharpness = 50.0
    current_loss = 2.0
    
    for t in range(n_steps):
        loss.append(current_loss)
        sharpness.append(current_sharpness)
        
        # 模拟动态
        # 1. Loss下降
        current_loss *= (1 - lr * 0.5)
        current_loss = max(0.01, current_loss)
        
        # 2. Sharpness动态
        if t < 100:
            # 初期增加
            current_sharpness *= 1.02
        elif t < 500:
            # 接近临界
            current_sharpness *= (0.9995 if current_sharpness > 2.0/lr else 1.001)
        else:
            # 收敛
            current_sharpness *= 0.9999
        
        current_sharpness = max(0.1, current_sharpness)
    
    return loss, sharpness

鞍点逃离机制

鞍点的几何结构

在高维空间中，鞍点具有以下性质：

• 指数多的负曲率方向：存在大量负特征值方向
• 能量壁垒较低：从一个局部最小值到另一个的路径上鞍点壁垒不高
• 维度依赖逃离难度：逃离时间随维度指数增长

逃离机制分析

def analyze_saddle_escape(steps=100, hessian_eigenvalues=None):
    """
    分析在给定Hessian谱下的鞍点逃离动态
    """
    if hessian_eigenvalues is None:
        # 使用典型的神经网络Hessian谱
        hessian_eigenvalues = np.concatenate([
            np.random.randn(10) * 10,  # 10个大特征值
            np.random.randn(90) * 0.1,  # 90个小特征值
        ])
    
    eigenvalues = np.array(hessian_eigenvalues)
    
    # 计算逃离时间估计
    # 对于负曲率方向 -λ，需要 η > 1/|λ| 才能逃离
    negative_eigenvalues = eigenvalues[eigenvalues < 0]
    positive_eigenvalues = eigenvalues[eigenvalues > 0]
    
    if len(negative_eigenvalues) > 0:
        hardest_escape = 1.0 / np.abs(np.min(negative_eigenvalues))
        avg_escape = 1.0 / np.abs(np.mean(negative_eigenvalues))
    else:
        hardest_escape = float('inf')
        avg_escape = float('inf')
    
    return {
        'n_negative_directions': len(negative_eigenvalues),
        'n_positive_directions': len(positive_eigenvalues),
        'hardest_escape_lr': hardest_escape,
        'avg_escape_lr': avg_escape,
        'saddle_index': len(negative_eigenvalues)
    }

随机梯度下降的鞍点逃离

SGD的噪声为逃离鞍点提供了机制：

噪声驱动的逃离：

$${\theta }_{t+1}={\theta }_t-\eta \nabla f({\theta }_t)+\sqrt{\eta }\sigma \cdot {\xi }_t$$

其中 ${\xi }_t\sim \mathcal{N}(0,I)$ 是噪声项。

def simulate_sgd_saddle_escape(
    n_steps=1000, 
    lr=0.01, 
    noise_std=0.01,
    initial_position=0.0
):
    """
    模拟SGD在鞍点附近的逃离
    """
    position = initial_position
    trajectory = [position]
    
    # 鞍点的Hessian: 正交方向正曲率，一个方向负曲率
    # 势能函数: f(x,y) = x² - y²
    def gradient(x, y):
        grad_x = 2 * x  # 正曲率方向
        grad_y = -2 * y  # 负曲率方向
        return grad_x, grad_y
    
    for t in range(n_steps):
        # 计算梯度
        grad_x, grad_y = gradient(position[0], position[1])
        
        # SGD更新
        noise = np.random.randn(2) * noise_std
        position = position - lr * np.array([grad_x, grad_y]) + noise
        
        trajectory.append(position.copy())
    
    return np.array(trajectory)

局部最小值的质量

好最小值 vs 坏最小值

并非所有局部最小值都一样好：

属性	好最小值	坏最小值
泛化能力	强	弱
Sharpness	小（平坦）	大（尖锐）
** Hessian谱**	均匀小特征值	分散大特征值
Mode Connectivity	高	低

Sharpness与泛化

Sharpness定义：

$$\text{Sharpness}({\theta }^{\ast })={\lambda }_{\max }(H({\theta }^{\ast }))$$

Sharpness-泛化关系：

经验发现：泛化误差与Sharpness正相关

def compute sharpness泛化_proxy(model, train_loader, test_loader, criterion):
    """
    计算Sharpness与泛化能力的代理指标
    """
    # 1. 计算训练和测试损失
    train_loss = evaluate_loss(model, train_loader, criterion)
    test_loss = evaluate_loss(model, test_loader, criterion)
    generalization_gap = test_loss - train_loss
    
    # 2. 估计Sharpness
    sharpness = estimate_sharpness(model, train_loader)
    
    # 3. 计算平坦度（small eigenvalue count）
    n_small_eigenvalues = count_small_eigenvalues(model, train_loader)
    
    return {
        'train_loss': train_loss,
        'test_loss': test_loss,
        'generalization_gap': generalization_gap,
        'sharpness': sharpness,
        'n_flat_directions': n_small_eigenvalues
    }

实用技术

1. 曲率感知的优化

class SharpnessAwareOptimizer(torch.optim.Optimizer):
    """
    曲率感知的优化器
    
    根据Hessian谱动态调整学习率
    """
    def __init__(self, params, lr=1e-3, sharpness_target=1.5):
        defaults = dict(lr=lr, sharpness_target=sharpness_target)
        super().__init__(params, defaults)
        self.sharpness_estimator = CurvatureEstimator()
    
    def step(self, closure=None):
        loss = None
        if closure is not None:
            loss = closure()
        
        # 估计当前sharpness
        sharpness = self.sharpness_estimator.estimate(self.param_groups)
        
        for group in self.param_groups:
            lr = group['lr']
            
            # 动态调整学习率
            if sharpness > group['sharpness_target']:
                adjusted_lr = lr * 0.9  # 减小学习率
            else:
                adjusted_lr = lr * 1.01  # 适当增大学习率
            
            # 执行更新
            for p in group['params']:
                if p.grad is not None:
                    p.data.add_(p.grad, alpha=-adjusted_lr)
        
        return loss

2. Fisher信息矩阵近似

使用Fisher信息矩阵 $F\approx \mathbb{E}[\nabla \log p(y|x)\nabla \log p(y|x{)}^T]$ 作为Hessian的替代：

class FisherCurvature:
    """
    Fisher信息矩阵曲率估计
    """
    def __init__(self, model, dataloader):
        self.model = model
        self.dataloader = dataloader
        self.fisher = None
    
    def compute_fisher(self, n_samples=1000):
        """计算对角Fisher信息"""
        self.fisher = {}
        
        for name, param in self.model.named_parameters():
            self.fisher[name] = torch.zeros_like(param)
        
        self.model.eval()
        for i, (inputs, targets) in enumerate(self.dataloader):
            if i >= n_samples:
                break
            
            self.model.zero_grad()
            outputs = self.model(inputs)
            loss = F.cross_entropy(outputs, targets)
            loss.backward()
            
            for name, param in self.model.named_parameters():
                if param.grad is not None:
                    self.fisher[name] += param.grad ** 2
        
        # 归一化
        for name in self.fisher:
            self.fisher[name] /= n_samples
    
    def estimate_natural_gradient(self, gradient):
        """计算自然梯度近似"""
        natural_gradient = {}
        idx = 0
        
        for name, grad in gradient.items():
            if name in self.fisher:
                # F⁻¹∇ = ∇ / F
                natural_gradient[name] = grad / (self.fisher[name] + 1e-8)
        
        return natural_gradient

参考

---

相关阅读

• 损失景观的层次结构 — 嵌入原理
• Sharp vs Flat Minima — 局部极小的曲率与泛化
• 模式连接理论 — 极小值之间的连通路径
• 训练不稳定性与平坦性偏差 — Edge of Stability深入分析

Footnotes

1. Ghorbani et al., “An Investigation into the Neural Tangent Kernel”, ICML 2020 ↩

2. Cohen et al., “Gradient Descent on Neural Networks Typically Occurs at the Edge of Stability”, ICLR 2021 ↩