auto_math_trust_regions_and_optimization_3

Snippet from Math-Trust-Regions-and-Optimization.md

def dogleg_step(g: np.ndarray, B: np.ndarray, Delta: float) -> np.ndarray:
    """
    Compute dogleg trust region step.
    
    Args:
        g: Gradient (d,).
        B: Hessian approximation (d x d).
        Delta: Trust radius.
    
    Returns:
        Dogleg step p.
    """
    # Cauchy point (steepest descent)
    g_norm_sq = g @ g
    g_B_g = g @ B @ g
    
    if g_B_g <= 0:
        tau = 1.0
    else:
        # Nocedal & Wright (2006) eq. 4.12: ||g||^3 == (g@g)**1.5, NOT (g@g)**3
        tau = min(g_norm_sq**1.5 / (Delta * g_B_g), 1.0)
    
    p_C = -tau * (Delta / np.sqrt(g_norm_sq)) * g
    
    # Newton point
    try:
        p_N = -np.linalg.solve(B, g)
    except np.linalg.LinAlgError:
        return p_C  # Fallback to Cauchy
    
    # Check if Newton within trust region
    if np.linalg.norm(p_N) <= Delta:
        return p_N
    
    # Dogleg interpolation
    # Find tau in [0, 2] where ||p_C + (tau - 1)(p_N - p_C)|| = Delta
    diff = p_N - p_C
    a = diff @ diff
    b = 2 * p_C @ diff
    c = p_C @ p_C - Delta**2
    
    discriminant = b**2 - 4*a*c
    if discriminant < 0:
        return p_C
    
    tau = (-b + np.sqrt(discriminant)) / (2*a)
    tau = np.clip(tau, 0, 1)
    
    return p_C + tau * diff