Newton's Method Second-order Optimization Techniques

Line Search

In the previous page we used a fixed learning rate (or step size) \(\eta\) chosen by experiment. A principled alternative is the line search: given a current iterate \(\boldsymbol{\theta}_t\) and a descent direction \(\mathbf{d}_t\) (i.e., one with \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) < 0\)), select the step that minimizes the objective along the ray \(\{\boldsymbol{\theta}_t + \eta\, \mathbf{d}_t : \eta > 0\}\): \[ \eta_t \in \operatorname*{arg\,min}_{\eta > 0}\, \mathcal{L}(\boldsymbol{\theta}_t + \eta\, \mathbf{d}_t). \]

Exact Line Search on a Quadratic

On a strictly convex quadratic \[ \mathcal{L}(\boldsymbol{\theta}) = \tfrac{1}{2}\, \boldsymbol{\theta}^\top \mathbf{A} \boldsymbol{\theta} + \mathbf{b}^\top \boldsymbol{\theta} + c, \qquad \mathbf{A} \succ 0, \] the one-dimensional problem admits a closed-form solution. Differentiating with respect to \(\eta\): \[ \frac{d}{d\eta}\, \mathcal{L}(\boldsymbol{\theta} + \eta\, \mathbf{d}) = \mathbf{d}^\top (\mathbf{A} \boldsymbol{\theta} + \mathbf{b}) + \eta\, \mathbf{d}^\top \mathbf{A}\, \mathbf{d}. \] Setting this to zero and using \(\mathbf{A} \boldsymbol{\theta} + \mathbf{b} = \nabla \mathcal{L}(\boldsymbol{\theta})\): \[ \eta^*_t = -\,\frac{\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t)}{\mathbf{d}_t^\top \mathbf{A}\, \mathbf{d}_t}. \] Positivity of \(\mathbf{A}\) makes the denominator positive, and the descent-direction condition makes the numerator negative, so \(\eta^*_t > 0\) as required.

Inexact Line Search

Solving the subproblem exactly at every outer iteration is wasteful. Inexact line searches instead seek any step \(\eta\) that delivers "sufficient" progress along \(\mathbf{d}_t\). The most common such rule is the Armijo condition.

The standard backtracking implementation starts from an initial trial step \(\eta\) (often \(\eta = 1\) for Newton-type methods) and shrinks it by a factor \(\tau \in (0, 1)\) (commonly \(\tau = 0.5\)) until (1) holds. By continuity of \(\mathcal{L}\) and the negativity of \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t)\), (1) is eventually satisfied for sufficiently small \(\eta\).

Newton's Method

Gradient descent uses a first-order linear model of \(\mathcal{L}\); Newton's method uses the second-order quadratic model. Instead of following the negative gradient, at each step it jumps to the minimizer of the local quadratic approximation, exploiting curvature information encoded in the Hessian.

The second-order Taylor expansion around the current iterate \(\boldsymbol{\theta}_t\) reads \[ \mathcal{L}(\boldsymbol{\theta}_t + \Delta\boldsymbol{\theta}) \;\approx\; \mathcal{L}(\boldsymbol{\theta}_t) + \mathbf{g}_t^\top \Delta\boldsymbol{\theta} + \tfrac{1}{2}\, \Delta\boldsymbol{\theta}^\top \mathbf{H}_t\, \Delta\boldsymbol{\theta}, \] where \(\mathbf{g}_t = \nabla \mathcal{L}(\boldsymbol{\theta}_t)\) and \(\mathbf{H}_t = \nabla^2 \mathcal{L}(\boldsymbol{\theta}_t)\). Minimizing the right-hand side over \(\Delta\boldsymbol{\theta}\) leads to the Newton step.

Pure Newton's method takes \(\boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t + \mathbf{d}_t\) (unit step). This scheme enjoys local quadratic convergence near a strict local minimum, but can diverge when started far from the optimum or when the Hessian is not positive definite. To restore global convergence, the Damped Newton's Method scales the Newton direction by a line-search step size: \[ \boldsymbol{\theta}_{t+1} \;=\; \boldsymbol{\theta}_t + \eta_t\, \mathbf{d}_t, \qquad \eta_t \text{ chosen to satisfy the } \href{#D-armijo_condition}{\textbf{Armijo condition}}. \]

BFGS Method

To avoid the \(O(n^3)\) cost of Hessian computation and inversion, quasi-Newton methods maintain a symmetric positive-definite matrix \(\mathbf{B}_t \approx \mathbf{H}_t\) (or its inverse \(\mathbf{C}_t \approx \mathbf{H}_t^{-1}\)) that is updated cheaply at each iteration using only gradient information. The BFGS method of Broyden–Fletcher–Goldfarb–Shanno is the quasi-Newton method most widely used in practice.

Secant Condition

Let \(\mathbf{s}_t = \boldsymbol{\theta}_{t+1} - \boldsymbol{\theta}_t\) denote the parameter step and \(\mathbf{y}_t = \mathbf{g}_{t+1} - \mathbf{g}_t\) the gradient change. A first-order Taylor expansion of \(\nabla \mathcal{L}\) around \(\boldsymbol{\theta}_t\) gives \[ \mathbf{g}_{t+1} \approx \mathbf{g}_t + \mathbf{H}_t\, \mathbf{s}_t \quad\Longrightarrow\quad \mathbf{y}_t \approx \mathbf{H}_t\, \mathbf{s}_t. \] This motivates the quasi-Newton curvature-matching requirement on the updated approximation:

Derivation of the BFGS Update

The secant condition alone under-determines \(\mathbf{B}_{t+1}\): in \(n\) dimensions there are \(n(n+1)/2\) free parameters in a symmetric matrix but only \(n\) linear constraints. BFGS makes the update unique by requiring (i) the update be a symmetric rank-2 correction of \(\mathbf{B}_t\), and (ii) \(\mathbf{B}_{t+1} \succ 0\) be preserved. Writing \(\mathbf{B}_{t+1} = \mathbf{B}_t + \mathbf{U}\), the secant condition forces \(\mathbf{U}\, \mathbf{s}_t = \mathbf{y}_t - \mathbf{B}_t\, \mathbf{s}_t\), and the rank-2 ansatz \[ \mathbf{U} = \alpha\, \mathbf{y}_t \mathbf{y}_t^\top + \beta\, (\mathbf{B}_t \mathbf{s}_t)(\mathbf{B}_t \mathbf{s}_t)^\top \] is chosen so that the first term injects curvature along \(\mathbf{y}_t\) while the second removes the now-stale curvature along \(\mathbf{B}_t \mathbf{s}_t\). Enforcing \(\mathbf{U}\, \mathbf{s}_t = \mathbf{y}_t - \mathbf{B}_t \mathbf{s}_t\) pins the scalars to \(\alpha = 1/(\mathbf{y}_t^\top \mathbf{s}_t)\) and \(\beta = -1/(\mathbf{s}_t^\top \mathbf{B}_t \mathbf{s}_t)\), yielding the BFGS update formula.

The curvature condition \(\mathbf{y}_t^\top \mathbf{s}_t > 0\) is not automatic — it fails on non-convex objectives unless the step size \(\eta_t\) is chosen carefully. It is guaranteed by the Wolfe conditions introduced below, which is why Wolfe (not plain Armijo) is the standard line search for BFGS.

Wolfe Conditions

BFGS's positive-definiteness preservation and fast convergence both depend on the step \(\eta_t\) satisfying conditions stronger than Armijo alone. The Wolfe conditions pair the sufficient-decrease condition with a curvature condition that prevents steps that are too short.

The curvature condition (W2) is precisely what guarantees the BFGS prerequisite \(\mathbf{y}_t^\top \mathbf{s}_t > 0\). Indeed, for a descent direction \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) < 0\), (W2) gives \[ \mathbf{y}_t^\top \mathbf{s}_t = \eta_t\, \mathbf{d}_t^\top \bigl(\nabla \mathcal{L}(\boldsymbol{\theta}_{t+1}) - \nabla \mathcal{L}(\boldsymbol{\theta}_t)\bigr) \geq \eta_t\, (c_2 - 1)\, \mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) > 0, \] since \(\eta_t > 0\), \(c_2 < 1\), and \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) < 0\). This is the reason BFGS line searches are implemented with (strong) Wolfe rather than plain Armijo.

Limited-Memory BFGS (L-BFGS)

For large-scale problems, storing the \(n \times n\) matrix \(\mathbf{C}_t\) explicitly is infeasible — \(n \sim 10^6\) already requires terabytes. L-BFGS sidesteps this by never forming \(\mathbf{C}_t\): it retains only the most recent \(m\) pairs \(\{(\mathbf{s}_i, \mathbf{y}_i)\}\) (typically \(m \in \{5, \ldots, 20\}\)) and computes \(\mathbf{C}_t\, \mathbf{g}_t\) implicitly via the two-loop recursion below. A direct expansion of the two-factor form \(\mathbf{C}_{t+1} = (\mathbf{I} - \rho_t \mathbf{s}_t \mathbf{y}_t^\top)\, \mathbf{C}_t\, (\mathbf{I} - \rho_t \mathbf{y}_t \mathbf{s}_t^\top) + \rho_t \mathbf{s}_t \mathbf{s}_t^\top\) over the last \(m\) pairs, applied to a vector \(\mathbf{g}_t\) against an implicit initial \(\mathbf{C}_t^{(0)}\) (often \(\mathbf{I}\)), produces exactly the two-loop recursion, and note that the recursion is also verifiable mechanically against the explicit formula, as done in our verification. The memory cost drops to \(O(mn)\) and each iteration costs \(O(mn)\).

Note: For numerical efficiency, the scalar \(\rho_t = 1/(\mathbf{y}_t^\top \mathbf{s}_t)\) is computed once and stored alongside each \((\mathbf{s}_t, \mathbf{y}_t)\) pair rather than being recomputed whenever it is needed in the two loops.

Python Implementation

The sample implementation below tests L-BFGS on the Rosenbrock function \[ f(\mathbf{x}) = \sum_{i=1}^{n-1} \bigl[a\, (x_{i+1} - x_i^2)^2 + (b - x_i)^2\bigr], \qquad a = 100,\; b = 1, \] a standard optimization benchmark with global minimum \(f(\mathbf{1}) = 0\). Despite the simple appearance, the Rosenbrock function is notoriously difficult for optimization algorithms because its global minimum lies at the bottom of a narrow, curved, "banana-shaped" valley. The surrounding landscape is ill-conditioned (strong anisotropy of the Hessian) and contains broad flat regions where gradients are small, so naive gradient descent crawls while Newton-type methods must handle curvature carefully.

                            import numpy as np

                            # Line search with Wolfe conditions 
                            def line_search(f, grad_f, theta, p, c1 = 1e-4, c2 = 0.9, max_iter = 100):
                                eta = 1.0
                                eta_low = 0.0
                                eta_high = None
                            
                                phi_0 = f(theta)
                                grad_phi_0 = np.dot(grad_f(theta), p)
                            
                                for _ in range(max_iter):
                                    phi_eta = f(theta + eta * p)
                            
                                    # Check Armijo condition
                                    if phi_eta > phi_0 + c1 * eta * grad_phi_0:
                                        eta_high = eta
                                    else:
                                        # Check Curvature condition
                                        grad_phi_eta = np.dot(grad_f(theta + eta * p), p)
                                        if grad_phi_eta < c2 * grad_phi_0:
                                            eta_low = eta
                                        else:
                                            return eta
                            
                                    # Update step size using bisection method
                                    if eta_high is not None:
                                        eta = (eta_low + eta_high) / 2.0
                                    else:
                                        eta *= 2.0
                            
                                return eta
                            
                            # Limited memory BFGS 
                            def limited_bfgs(f, grad_f, theta0, m = 10, tol = 1e-6, max_iter = 2000):
                                
                                theta = theta0.copy()
                                g = grad_f(theta)
                                s_list = []
                                y_list = []
                                rho_list = [] # We introduce rho =  1/s^\top y instead of directly using 1/s^\top y for efficiency & stability. 
                            
                                for _ in range(max_iter):
                                    
                                    if np.linalg.norm(g) < tol:
                                        break
                            
                                    q = g.copy()
                                    alpha_list = []  # Need this for the second loop to avoid computing alpha again 
                            
                                    # Loop backward through stored (s, y) pairs
                                    for s, y, rho in reversed(list(zip(s_list, y_list, rho_list))):
                                        alpha = rho * np.dot(s, q)
                                        alpha_list.append(alpha)
                                        q -= alpha * y
                            
                                    # Initial Hessian approximation is identity: H0 = I, so H0 * q = q
                                    # Note: In the first iteration (when s_list is empty), the loops do nothing,
                                    # so r = q = g, and thus p = -r = -g (steepest descent direction)
                                    r = q
                            
                                    # Loop forward through stored (s, y) pairs
                                    for (s, y, rho), alpha in zip(zip(s_list, y_list, rho_list), reversed(alpha_list)):
                                        beta = rho * np.dot(y, r)
                                        r += s * (alpha - beta)
                            
                                    # Search direction
                                    p = -r
                            
                                    # Compute the step size by Line search satisfying Wolfe conditions
                                    eta = line_search(f, grad_f, theta, p)
                            
                                    # Update parameters
                                    theta += eta * p
                                    grad_next = grad_f(theta)
                            
                                    # Update memory for (s, y) pairs
                                    s = eta * p
                                    y = grad_next - g
                                    if np.dot(s, y) > 1e-10 * np.linalg.norm(s) * np.linalg.norm(y): 
                                        if len(s_list) == m:
                                            s_list.pop(0)
                                            y_list.pop(0)
                                            rho_list.pop(0)
                                        s_list.append(s)
                                        y_list.append(y)
                                        rho_list.append(1.0 / np.dot(y, s))
                                        
                                    # Update gradient
                                    g = grad_next
                            
                                return theta
                            
                            # Objective function and its gradient: 
                            # The Rosenbrock function is commonly used for testing optimization algorithms.
                            def rosenbrock(x):
                                return np.sum(100 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)
                            
                            def grad_rosenbrock(x):
                                grad = np.zeros_like(x)
                                grad[:-1] = -400 * x[:-1] * (x[1:] - x[:-1]**2) - 2 * (1 - x[:-1]) # # x_1 to x_n-1
                                grad[1:] += 200 * (x[1:] - x[:-1]**2) # x_2 to x_n
                                return grad
                            
                            # Finite difference gradient to check grad_rosenbrock().
                            def finite_difference_gradient(f, x, epsilon=1e-6):
                                grad = np.zeros_like(x)
                                for i in range(len(x)):
                                    x_forward = x.copy()
                                    x_backward = x.copy()
                                    x_forward[i] += epsilon
                                    x_backward[i] -= epsilon
                                    grad[i] = (f(x_forward) - f(x_backward)) / (2 * epsilon)
                                return grad
                            
                            if __name__ == "__main__":
                                
                                n = 50 # Dimensionality
                                
                                # Randomly generate an initial point x0:
                                # You could try : x0 = np.ones(n) + 0.3 * np.random.randn(n), which represents adding small 
                                # perturbation around the global minimum x* = [1, ... , 1]^\top 
                                x0 =  np.random.randn(n)  
                                numeric_opt = limited_bfgs(rosenbrock, grad_rosenbrock, x0)
                                print("Initial point: \n", x0.tolist())
                                print("\n Numerical optimum: \n", numeric_opt.tolist())
                            
                                '''
                                # If you are not sure about the gradient of the objective, always you can compare it with finite difference.
                                grad_analytic = grad_rosenbrock(x0)
                                grad_numeric = finite_difference_gradient(rosenbrock, x0)
                                relative_error_grad = np.linalg.norm(grad_analytic - grad_numeric) / np.linalg.norm(grad_analytic)
                                print("Relative error:", relative_error_grad) 
                                ''' 
                                # Relative error between the numerical optimum and the global optimum
                                global_opt = np.ones(n)  # The actual global optimum of Rosenbrock function is x* = [1, ... , 1]^\top
                                relative_error_optimal = np.linalg.norm(numeric_opt - global_opt) / np.linalg.norm(global_opt)
                                print(f"\n Relative error to the global minimum: {relative_error_optimal*100:.8f}%")