Advanced Optimization Methods

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

import numpy as np
from typing import Tuple
 
def bfgs_update(
    H: np.ndarray,
    s: np.ndarray,
    y: np.ndarray,
    eps: float = 1e-10
) -> np.ndarray:
    """
    Compute BFGS update of inverse Hessian approximation.
    
    H_{k+1} = (I - ρ·s·yᵀ) H_k (I - ρ·y·sᵀ) + ρ·s·sᵀ
    
    Parameters:
    -----------
    H : np.ndarray
        Current inverse Hessian approximation (n, n)
    s : np.ndarray
        Step vector: s = x_{k+1} - x_k (n,)
    y : np.ndarray
        Gradient difference: y = ∇f_{k+1} - ∇f_k (n,)
    eps : float
        Small constant for numerical stability
        
    Returns:
    --------
    H_new : np.ndarray
        Updated inverse Hessian approximation (n, n)
    """
    ys = np.dot(y, s)
    
    # Curvature condition check
    if ys < eps:
        # Skip update if curvature condition violated
        return H
    
    rho = 1.0 / ys
    n = len(s)
    I = np.eye(n)
    
    # BFGS formula: H = (I - ρsyᵀ) H (I - ρysᵀ) + ρssᵀ
    V = I - rho * np.outer(s, y)
    H_new = V @ H @ V.T + rho * np.outer(s, s)
    
    return H_new
 
 
def bfgs_direct_hessian_update(
    B: np.ndarray,
    s: np.ndarray,
    y: np.ndarray,
    eps: float = 1e-10
) -> np.ndarray:
    """
    Direct BFGS update for Hessian approximation B.
    
    B_{k+1} = B_k - (B_k s sᵀ B_k)/(sᵀ B_k s) + (y yᵀ)/(yᵀ s)
    
    This is the Sherman-Morrison-Woodbury dual of the inverse update.
    """
    ys = np.dot(y, s)
    
    if ys < eps:
        return B
    
    Bs = B @ s
    sBs = np.dot(s, Bs)
    
    if abs(sBs) < eps:
        return B
    
    B_new = B - np.outer(Bs, Bs) / sBs + np.outer(y, y) / ys
    
    return B_new

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

import numpy as np
from typing import Callable, Tuple, Optional
 
def wolfe_line_search(
    f: Callable[[np.ndarray], float],
    grad_f: Callable[[np.ndarray], np.ndarray],
    x: np.ndarray,
    d: np.ndarray,
    grad: np.ndarray,
    c1: float = 1e-4,
    c2: float = 0.9,
    alpha_init: float = 1.0,
    max_iter: int = 25
) -> Tuple[float, np.ndarray]:
    """
    Line search satisfying strong Wolfe conditions.
    
    Returns step size α and gradient at new point.
    Critical for BFGS to maintain positive definiteness.
    """
    f_0 = f(x)
    g_0 = np.dot(grad, d)  # Directional derivative at α=0
    
    alpha_lo, alpha_hi = 0.0, np.inf
    alpha = alpha_init
    
    for _ in range(max_iter):
        x_new = x + alpha * d
        f_new = f(x_new)
        grad_new = grad_f(x_new)
        g_new = np.dot(grad_new, d)
        
        # Check Armijo condition
        if f_new > f_0 + c1 * alpha * g_0:
            alpha_hi = alpha
            alpha = 0.5 * (alpha_lo + alpha_hi)
            continue
            
        # Check curvature condition (strong Wolfe)
        if abs(g_new) > c2 * abs(g_0):
            if g_new * g_0 < 0:
                alpha_hi = alpha
            else:
                alpha_lo = alpha
            
            if alpha_hi < np.inf:
                alpha = 0.5 * (alpha_lo + alpha_hi)
            else:
                alpha = 2 * alpha
            continue
        
        # Both conditions satisfied
        return alpha, grad_new
    
    return alpha, grad_f(x + alpha * d)
 
 
def bfgs_optimize(
    f: Callable[[np.ndarray], float],
    grad_f: Callable[[np.ndarray], np.ndarray],
    x0: np.ndarray,
    tol: float = 1e-8,
    max_iter: int = 1000,
    verbose: bool = False
) -> Tuple[np.ndarray, list]:
    """
    BFGS quasi-Newton optimization with Wolfe line search.
    
    Parameters:
    -----------
    f : callable
        Objective function f(x) -> scalar
    grad_f : callable
        Gradient function ∇f(x) -> n-vector
    x0 : np.ndarray
        Starting point (n,)
    tol : float
        Gradient norm tolerance for convergence
    max_iter : int
        Maximum number of iterations
    verbose : bool
        Print iteration information
        
    Returns:
    --------
    x : np.ndarray
        Approximate minimizer
    history : list
        List of (x, f(x), ||grad||) at each iteration
    """
    n = len(x0)
    x = x0.copy().astype(float)
    H = np.eye(n)  # Initial inverse Hessian: identity
    
    grad = grad_f(x)
    history = [(x.copy(), f(x), np.linalg.norm(grad))]
    
    for k in range(max_iter):
        grad_norm = np.linalg.norm(grad)
        
        if verbose and k % 10 == 0:
            print(f"Iter {k:4d}: f = {f(x):.6e}, ||∇f|| = {grad_norm:.6e}")
        
        # Check convergence
        if grad_norm < tol:
            if verbose:
                print(f"Converged in {k} iterations")
            break
        
        # Compute search direction
        d = -H @ grad
        
        # Line search (Wolfe conditions)
        alpha, grad_new = wolfe_line_search(f, grad_f, x, d, grad)
        
        # Compute step and gradient difference
        s = alpha * d
        y = grad_new - grad
        
        # Update position and gradient
        x = x + s
        grad = grad_new
        
        # BFGS update
        ys = np.dot(y, s)
        if ys > 1e-10:  # Curvature condition
            rho = 1.0 / ys
            I = np.eye(n)
            V = I - rho * np.outer(s, y)
            H = V @ H @ V.T + rho * np.outer(s, s)
        
        history.append((x.copy(), f(x), np.linalg.norm(grad)))
    
    return x, history
 
 
# Example: Rosenbrock function
def rosenbrock(x):
    """Rosenbrock function (extended to n dimensions)"""
    return sum(100 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2 
               for i in range(len(x) - 1))
 
def rosenbrock_grad(x):
    """Gradient of Rosenbrock function"""
    n = len(x)
    grad = np.zeros(n)
    for i in range(n - 1):
        grad[i] += -400 * x[i] * (x[i+1] - x[i]**2) - 2 * (1 - x[i])
        grad[i+1] += 200 * (x[i+1] - x[i]**2)
    return grad
 
 
if __name__ == "__main__":
    # 10-dimensional Rosenbrock
    x0 = np.zeros(10)
    
    x_opt, history = bfgs_optimize(
        rosenbrock, rosenbrock_grad, x0,
        verbose=True
    )
    
    print(f"\nOptimum: {x_opt}")
    print(f"Optimal value: {rosenbrock(x_opt):.2e}")
    print(f"Total iterations: {len(history)}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190

import numpy as np
from collections import deque
from typing import Callable, Tuple, List
 
def lbfgs_two_loop_recursion(
    grad: np.ndarray,
    s_history: List[np.ndarray],
    y_history: List[np.ndarray],
    rho_history: List[float],
    gamma: float = 1.0
) -> np.ndarray:
    """
    L-BFGS two-loop recursion to compute search direction.
    
    Computes H_k @ grad without forming H_k explicitly.
    Uses only the m most recent (s, y, ρ) triplets.
    
    Parameters:
    -----------
    grad : np.ndarray
        Current gradient (n,)
    s_history : list
        Recent s_i = x_{i+1} - x_i vectors
    y_history : list
        Recent y_i = ∇f_{i+1} - ∇f_i vectors  
    rho_history : list
        Recent ρ_i = 1/(y_i·s_i) scalars
    gamma : float
        Scaling for initial H_0 = gamma * I
        
    Returns:
    --------
    d : np.ndarray
        Search direction -H_k @ grad (n,)
    """
    q = grad.copy()
    m = len(s_history)
    
    if m == 0:
        # No history: steepest descent
        return -gamma * grad
    
    # First loop: right-to-left
    alpha = np.zeros(m)
    for i in range(m - 1, -1, -1):
        alpha[i] = rho_history[i] * np.dot(s_history[i], q)
        q = q - alpha[i] * y_history[i]
    
    # Initial direction: H_0 @ q where H_0 = gamma * I
    r = gamma * q
    
    # Second loop: left-to-right
    for i in range(m):
        beta = rho_history[i] * np.dot(y_history[i], r)
        r = r + (alpha[i] - beta) * s_history[i]
    
    return -r
 
 
def lbfgs_optimize(
    f: Callable[[np.ndarray], float],
    grad_f: Callable[[np.ndarray], np.ndarray],
    x0: np.ndarray,
    m: int = 10,
    tol: float = 1e-8,
    max_iter: int = 1000,
    verbose: bool = False
) -> Tuple[np.ndarray, list]:
    """
    L-BFGS optimization with limited memory.
    
    Parameters:
    -----------
    f : callable
        Objective function f(x) -> scalar
    grad_f : callable
        Gradient function ∇f(x) -> n-vector
    x0 : np.ndarray
        Starting point (n,)
    m : int
        Number of (s, y) pairs to store (memory parameter)
    tol : float
        Gradient norm tolerance for convergence
    max_iter : int
        Maximum number of iterations
    verbose : bool
        Print iteration information
        
    Returns:
    --------
    x : np.ndarray
        Approximate minimizer
    history : list
        List of (x, f(x), ||grad||) at each iteration
    """
    x = x0.copy().astype(float)
    
    # History storage (using deque for efficient pop-left)
    s_history = deque(maxlen=m)
    y_history = deque(maxlen=m)
    rho_history = deque(maxlen=m)
    
    grad = grad_f(x)
    history = [(x.copy(), f(x), np.linalg.norm(grad))]
    
    for k in range(max_iter):
        grad_norm = np.linalg.norm(grad)
        
        if verbose and k % 20 == 0:
            print(f"Iter {k:4d}: f = {f(x):.6e}, ||∇f|| = {grad_norm:.6e}")
        
        # Check convergence
        if grad_norm < tol:
            if verbose:
                print(f"Converged in {k} iterations")
            break
        
        # Compute initial H_0 scaling
        if len(y_history) > 0:
            # Barzilai-Borwein scaling
            gamma = (np.dot(s_history[-1], y_history[-1]) / 
                     np.dot(y_history[-1], y_history[-1]))
        else:
            gamma = 1.0 / grad_norm  # Scale by inverse gradient norm initially
        
        # Compute search direction using two-loop recursion
        d = lbfgs_two_loop_recursion(
            grad, list(s_history), list(y_history), 
            list(rho_history), gamma
        )
        
        # Line search (simplified backtracking for brevity)
        alpha = 1.0
        c1 = 1e-4
        f_x = f(x)
        grad_d = np.dot(grad, d)
        
        for _ in range(20):
            if f(x + alpha * d) <= f_x + c1 * alpha * grad_d:
                break
            alpha *= 0.5
        
        # Update
        x_new = x + alpha * d
        grad_new = grad_f(x_new)
        
        s = x_new - x
        y = grad_new - grad
        ys = np.dot(y, s)
        
        # Store new (s, y, ρ) if curvature condition satisfied
        if ys > 1e-10:
            s_history.append(s)
            y_history.append(y)
            rho_history.append(1.0 / ys)
        
        x = x_new
        grad = grad_new
        
        history.append((x.copy(), f(x), np.linalg.norm(grad)))
    
    return x, history
 
 
if __name__ == "__main__":
    # Test on high-dimensional Rosenbrock
    n = 1000
    
    def rosenbrock_nd(x):
        return sum(100 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2 
                   for i in range(len(x) - 1))
    
    def rosenbrock_grad_nd(x):
        n = len(x)
        grad = np.zeros(n)
        for i in range(n - 1):
            grad[i] += -400 * x[i] * (x[i+1] - x[i]**2) - 2 * (1 - x[i])
            grad[i+1] += 200 * (x[i+1] - x[i]**2)
        return grad
    
    x0 = np.zeros(n)
    
    x_opt, history = lbfgs_optimize(
        rosenbrock_nd, rosenbrock_grad_nd, x0,
        m=20, verbose=True, max_iter=500
    )
    
    print(f"\nFinal gradient norm: {np.linalg.norm(rosenbrock_grad_nd(x_opt)):.2e}")
    print(f"Optimal value: {rosenbrock_nd(x_opt):.2e}")
    print(f"Total iterations: {len(history)}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

import numpy as np
from scipy.optimize import minimize
 
# Example 1: Logistic Regression with L-BFGS
def logistic_loss_and_grad(w, X, y, reg=0.01):
    """
    Compute L2-regularized logistic loss and gradient.
    
    Loss = Σ log(1 + exp(-y_i * wᵀx_i)) + (reg/2)||w||²
    """
    N = len(y)
    z = y * (X @ w)
    
    # Numerically stable log-sum-exp
    loss = np.sum(np.logaddexp(0, -z)) / N + 0.5 * reg * np.dot(w, w)
    
    # Gradient
    sigmoid_neg = 1 / (1 + np.exp(z))  # σ(-z) = 1 - σ(z)
    grad = -X.T @ (y * sigmoid_neg) / N + reg * w
    
    return loss, grad
 
 
def train_logistic_lbfgs(X, y, reg=0.01):
    """Train logistic regression using scipy's L-BFGS-B."""
    n_features = X.shape[1]
    w0 = np.zeros(n_features)
    
    result = minimize(
        logistic_loss_and_grad,
        w0,
        args=(X, y, reg),
        method='L-BFGS-B',
        jac=True,  # Function returns (loss, gradient)
        options={'maxiter': 1000, 'disp': True}
    )
    
    return result.x
 
 
# Example 2: Neural Network training with L-BFGS
def mlp_forward_backward(params, X, y, hidden_size=100):
    """
    Two-layer MLP: X -> hidden -> softmax output
    Returns loss and flattened gradient.
    """
    n_features = X.shape[1]
    n_classes = np.max(y) + 1
    
    # Unpack parameters
    W1_size = n_features * hidden_size
    b1_size = hidden_size
    W2_size = hidden_size * n_classes
    
    W1 = params[:W1_size].reshape(n_features, hidden_size)
    b1 = params[W1_size:W1_size + b1_size]
    W2 = params[W1_size + b1_size:W1_size + b1_size + W2_size].reshape(hidden_size, n_classes)
    b2 = params[W1_size + b1_size + W2_size:]
    
    # Forward pass
    h = np.maximum(0, X @ W1 + b1)  # ReLU
    logits = h @ W2 + b2
    
    # Softmax and cross-entropy
    logits_max = logits.max(axis=1, keepdims=True)
    exp_logits = np.exp(logits - logits_max)
    probs = exp_logits / exp_logits.sum(axis=1, keepdims=True)
    
    N = len(y)
    loss = -np.log(probs[np.arange(N), y] + 1e-8).mean()
    
    # Backward pass
    dlogits = probs.copy()
    dlogits[np.arange(N), y] -= 1
    dlogits /= N
    
    dW2 = h.T @ dlogits
    db2 = dlogits.sum(axis=0)
    
    dh = dlogits @ W2.T
    dh[h <= 0] = 0  # ReLU gradient
    
    dW1 = X.T @ dh
    db1 = dh.sum(axis=0)
    
    # Pack gradients
    grad = np.concatenate([dW1.ravel(), db1, dW2.ravel(), db2])
    
    return loss, grad
 
 
# Example 3: Using PyTorch's L-BFGS
"""
import torch
from torch.optim import LBFGS
 
model = MyModel()
optimizer = LBFGS(
    model.parameters(),
    lr=1.0,
    max_iter=20,
    history_size=10,
    line_search_fn='strong_wolfe'
)
 
def closure():
    optimizer.zero_grad()
    output = model(input)
    loss = criterion(output, target)
    loss.backward()
    return loss
 
# L-BFGS requires calling optimizer.step(closure) with a closure
for epoch in range(num_epochs):
    optimizer.step(closure)
"""