Newton's Method Second-order Optimization Techniques

Line Search

Our treatment of gradient descent 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\) (that is, one with \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) \lt 0\)), we select a step that minimizes the objective along the ray \(\{\boldsymbol{\theta}_t + \eta\, \mathbf{d}_t : \eta \gt 0\}\), assuming such a step exists: \[ \eta_t \in \operatorname*{arg\,min}_{\eta \gt 0}\, \mathcal{L}(\boldsymbol{\theta}_t + \eta\, \mathbf{d}_t). \]

Exact Line Search on a Quadratic

On a strictly convex quadratic \[ \begin{align*} \mathcal{L}(\boldsymbol{\theta}) &= \tfrac{1}{2}\, \boldsymbol{\theta}^\top \mathbf{A} \boldsymbol{\theta} + \mathbf{b}^\top \boldsymbol{\theta} + c, \\\\ \mathbf{A} &\succ 0, \end{align*} \] 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 \gt 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. Since \(\mathcal{L}\) is differentiable at \(\boldsymbol{\theta}_t\), \[ \mathcal{L}(\boldsymbol{\theta}_t + \eta\, \mathbf{d}_t) = \mathcal{L}(\boldsymbol{\theta}_t) + \eta\, \mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) + o(\eta), \] and because \(c_1 \lt 1\) and \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) \lt 0\), the leftover \((1 - c_1)\, \eta\, \mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t)\) dominates the \(o(\eta)\) term, so (1) holds for all 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 \[ \begin{align*} \mathcal{L}(\boldsymbol{\theta}_t + \Delta\boldsymbol{\theta}) &\approx \mathcal{L}(\boldsymbol{\theta}_t) + \mathbf{g}_t^\top \Delta\boldsymbol{\theta} \\\\ &\quad + \tfrac{1}{2}\, \Delta\boldsymbol{\theta}^\top \mathbf{H}_t\, \Delta\boldsymbol{\theta}, \end{align*} \] 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). Near a minimizer at which the Hessian is positive definite and Lipschitz continuous, this scheme converges quadratically (local quadratic convergence), a classical result we state without proof on this page, but can diverge when started far from the optimum or when the Hessian is not positive definite.

To recover global convergence (under standard assumptions that we do not prove here), 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, \quad \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 \[ \begin{align*} \mathbf{g}_{t+1} &\approx \mathbf{g}_t + \mathbf{H}_t\, \mathbf{s}_t \\\\ &\Longrightarrow \mathbf{y}_t \approx \mathbf{H}_t\, \mathbf{s}_t. \end{align*} \] The quasi-Newton family therefore imposes a 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 singles out one update 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 specific rank-2 ansatz built from \(\mathbf{y}_t\) and \(\mathbf{B}_t \mathbf{s}_t\), the two vectors this equation involves, \[ \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, by comparing coefficients when \(\mathbf{y}_t\) and \(\mathbf{B}_t \mathbf{s}_t\) are linearly independent (the formula is then adopted in all cases), 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 \gt 0\) is not automatic. It fails on non-convex objectives unless the step size \(\eta_t\) is chosen carefully. The Wolfe conditions introduced below guarantee it, 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 \gt 0\). Indeed, for a descent direction \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) \lt 0\), (W2) gives \[ \begin{align*} \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) \\\\ &\gt 0, \end{align*} \] since \(\eta_t \gt 0\), \(c_2 \lt 1\), and \(\mathbf{d}_t^\top \nabla \mathcal{L}(\boldsymbol{\theta}_t) \lt 0\).

That chain of inequalities 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. At \(n \sim 10^6\) the matrix 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.

Unrolling the two-factor recurrence 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 that recursion. The memory cost drops to \(O(mn)\) and each iteration costs \(O(mn)\).

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 \[ \begin{align*} f(\mathbf{x}) &= \sum_{i=1}^{n-1} \bigl[a\, (x_{i+1} - x_i^2)^2 + (b - x_i)^2\bigr], \\\\ a &= 100, \quad b = 1, \end{align*} \] 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 along the valley floor the gradient is 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}%")