# Import necessary libraries
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from IPython.display import display, Math

Introduction to Optimization

Optimization is the process of finding the best solution to a problem from a set of possible solutions. In mathematical terms, we seek to find the minimum or maximum of a function $f: \mathbb{R}^n \rightarrow \mathbb{R}$.

Key Concepts:

Objective Function: The function $f(x)$ we want to minimize or maximize
Decision Variables: The variables $x = (x_1, x_2, …, x_n)$ we can control
Optimal Solution: The point $x^{\ast}$ where $f(x^{\ast})$ is minimal (or maximal)

We can differentiate between:

Unconstrained Optimization: Find $\min_{x \in \mathbb{R}^n} f(x)$
Constrained Optimization: Find $\min_{x \in \mathbb{R}^n} f(x)$ subject to constraints $g_i(x) \leq 0$, $h_j(x) = 0$

Mathematical Foundations

Gradient

The gradient of a function $f: \mathbb{R}^n \rightarrow \mathbb{R}$ at a point $x$ is the vector of partial derivatives:

$$\nabla f(x) = \left[\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, …, \frac{\partial f}{\partial x_n}\right]^T$$

Intuition: The gradient points in the direction of steepest increase of the function. Its negative points in the direction of steepest decrease.

Hessian Matrix

The Hessian matrix is the matrix of second-order partial derivatives:

$$H(x) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{bmatrix}$$

Intuition: The Hessian describes the curvature of the function, helping us determine whether a point is a minimum, maximum, or saddle point.

Stationary Point Analysis

A stationary point (or critical point) is a point where the gradient equals zero:

$$\nabla f(x^*) = 0$$

Classification of Stationary Points

Using the Hessian matrix $H(x^*)$ at a stationary point:

Local Minimum: $H(x^*)$ is positive definite (all eigenvalues > 0)
Local Maximum: $H(x^*)$ is negative definite (all eigenvalues < 0)
Saddle Point: $H(x^*)$ has both positive and negative eigenvalues
Inconclusive: $H(x^*)$ is singular or has zero eigenvalues

Implementation

# Define symbols for symbolic computation
x, y = sp.symbols('x y')

# Example function: f(x,y) = x^2 + y^2
f_symbolic = x**2 + y**2
print("Function: f(x,y) =", f_symbolic)

# Compute gradient symbolically
gradient = [sp.diff(f_symbolic, x), sp.diff(f_symbolic, y)]
print("\nGradient: ∇f =", gradient)

# Solve for stationary points
stationary_points = sp.solve(gradient, (x, y), dict=True)
print("\nStationary Points:", stationary_points)

# Compute and display Hessian matrix
Hessian = [[sp.diff(f_symbolic, x, x), sp.diff(f_symbolic, x, y)],
           [sp.diff(f_symbolic, y, x), sp.diff(f_symbolic, y, y)]]
print("\nHessian Matrix:")
display(Math(sp.latex(sp.Matrix(Hessian))))

# Classify each stationary point
for point in stationary_points:
    # Evaluate Hessian at the point
    H_eval = [[expr.subs(point) for expr in row] for row in Hessian]
    H_matrix = sp.Matrix(H_eval)
    
    # Calculate eigenvalues
    eigenvalues = list(H_matrix.eigenvals().keys())
    eigenvalues_float = [float(sp.N(ev, chop=True)) for ev in eigenvalues]
    
    # Classify based on eigenvalues
    if all(ev > 0 for ev in eigenvalues_float):
        classification = "Local Minimum (convex)"
    elif all(ev < 0 for ev in eigenvalues_float):
        classification = "Local Maximum (concave)"
    elif any(ev > 0 for ev in eigenvalues_float) and any(ev < 0 for ev in eigenvalues_float):
        classification = "Saddle Point"
    else:
        classification = "Inconclusive"
    
    print(f"\nPoint {point}:")
    print(f"Classification: {classification}")
    print(f"Eigenvalues: {eigenvalues_float}")

Function: f(x,y) = x**2 + y**2

Gradient: ∇f = [2*x, 2*y]

Stationary Points: [{x: 0, y: 0}]

Hessian Matrix:

$\displaystyle \left[\begin{matrix}2 & 0\0 & 2\end{matrix}\right]$

Point {x: 0, y: 0}:
Classification: Local Minimum (convex)
Eigenvalues: [2.0]

Gradient Descent Algorithm

Gradient Descent is an iterative optimization algorithm for finding a local minimum of a differentiable function.

Formal Definition

Let $f: \mathbb{R}^n \rightarrow \mathbb{R}$ be a continuously differentiable function. The gradient descent algorithm generates a sequence ${x_k}_{k=0}^{\infty}$ according to the iteration:

$$x_{k+1} = x_k - \alpha_k \nabla f(x_k)$$

where:

$x_0 \in \mathbb{R}^n$ is the initial point
$\alpha_k > 0$ is the step size (or learning rate) at iteration $k$
$\nabla f(x_k) \in \mathbb{R}^n$ is the gradient of $f$ evaluated at $x_k$

Motivation

The gradient descent algorithm is based on the first-order Taylor approximation:

$$f(x + \Delta x) \approx f(x) + \nabla f(x)^T \Delta x$$

To minimize $f$, we want to choose $\Delta x$ such that $f(x + \Delta x) < f(x)$. The direction of steepest descent is:

$$\Delta x = -\alpha \nabla f(x)$$

for some $\alpha > 0$. This gives us the maximum decrease in $f$ for small steps.

Finite Difference Approximation

In practice, we often cannot compute $\nabla f$ analytically. We use the finite difference method to approximate partial derivatives:

Forward Difference:

$$\frac{\partial f}{\partial x_i} \approx \frac{f(x + h e_i) - f(x)}{h}$$

Central Difference (more accurate):

$$\frac{\partial f}{\partial x_i} \approx \frac{f(x + h e_i) - f(x - h e_i)}{2h}$$

where:

$h > 0$ is a small step size (typically $h \in [10^{-8}, 10^{-4}]$)
$e_i$ is the $i$-th standard basis vector

Learning Rate Properties

The choice of learning rate $\alpha_k$ is crucial:

Constant Learning Rate: $\alpha_k = \alpha$ for all $k$
- Simple but may not converge or converge slowly
Diminishing Learning Rate: $\alpha_k \rightarrow 0$ as $k \rightarrow \infty$
- Ensures convergence under certain conditions
- Example: $\alpha_k = \frac{\alpha_0}{1 + \beta k}$
Adaptive Learning Rate: Based on local properties of $f$
- Line search methods
- Trust region methods

Stopping Criteria

Common criteria for terminating the algorithm:

Gradient norm: $|\nabla f(x_k)| < \varepsilon$
Relative change in function value: $\frac{|f(x_{k+1}) - f(x_k)|}{|f(x_k)|} < \varepsilon$
Relative change in parameters: $\frac{|x_{k+1} - x_k|}{|x_k|} < \varepsilon$
Maximum iterations: $k \geq k_{\max}$

Challenges and Limitations

Local Minima: Gradient descent only guarantees convergence to local minima
Saddle Points: Can get stuck or slow down significantly
Ill-Conditioning: Poor performance when Hessian has large condition number
Step Size Selection: Too large → divergence; Too small → slow convergence

# Numerical functions for optimization
def numerical_gradient(f, x, h=1e-4):
    # Compute numerical gradient of f at point x
    grad = np.zeros_like(x)                             
    for i in range(len(x)):                                   
        x_plus, x_minus = x.copy(), x.copy()            
        x_plus[i] += h                            
        x_minus[i] -= h 
        grad[i] = (f(x_plus) - f(x_minus)) / (2 * h)    
    return grad

def gradient_descent(f, x0, alpha=0.1, max_iter=100, tol=1e-6): 
    # Perform gradient descent optimization
    x = x0.copy()
    history = [x.copy()]
    
    for i in range(max_iter):
        grad = numerical_gradient(f, x)
        if np.linalg.norm(grad) < tol:
            print(f"Converged in {i} iterations")
            break
        x -= alpha * grad
        history.append(x.copy())
    
    return x, history

# Example: Quadratic function f(x,y) = x² + y²
f = lambda p: p[0]**2 + p[1]**2
x0 = np.array([4.0, 3.0])

# Run gradient descent from one starting point
x_opt, history = gradient_descent(f, x0, alpha=0.1, max_iter=50)
print(f"Starting point: {x0}")
print(f"Optimal point: {np.round(x_opt, 6)}")
print(f"Optimal value: {f(x_opt):.6e}\n")

# Create visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Create meshgrid
x_vals = np.linspace(-5, 5, 100)
y_vals = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

# Left: Contour plot with gradient descent path
contours = ax1.contour(X, Y, Z, levels=15, cmap='viridis', linewidths=0.8, alpha=0.7)
ax1.clabel(contours, inline=True, fontsize=8, fmt='%g')

history = np.array(history)
ax1.plot(history[:, 0], history[:, 1], 'r-', linewidth=1.5, alpha=0.8)
ax1.plot(history[0, 0], history[0, 1], 'ro', markersize=7, label='Start')
ax1.plot(history[-1, 0], history[-1, 1], 'r^', markersize=7, label='Optimum')

ax1.set_xlabel('$x$', fontsize=11)
ax1.set_ylabel('$y$', fontsize=11)
ax1.set_title('Gradient Descent Path', fontsize=12)
ax1.legend(frameon=False)
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Right: 3D surface
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7, edgecolor='none', antialiased=True)
ax2.plot(history[:, 0], history[:, 1], 
         [f(np.array([x, y])) for x, y in zip(history[:, 0], history[:, 1])], 
         'r-', linewidth=2, alpha=0.9)
ax2.scatter(0, 0, 0, color='red', s=60, marker='^')

ax2.set_xlabel('$x$', fontsize=10)
ax2.set_ylabel('$y$', fontsize=10)
ax2.set_zlabel('$f(x,y)$', fontsize=10)
ax2.set_title('3D Surface', fontsize=12)
ax2.view_init(elev=20, azim=45)

plt.tight_layout()
plt.show()

Starting point: [4. 3.]
Optimal point: [5.7e-05 4.3e-05]
Optimal value: 5.092590e-09

png

Constrained Optimization: Lagrange Multipliers

Real optimization problems are rarely unconstrained. The method of Lagrange multipliers transforms these constrained problems into a form we can solve systematically.

Imagine you’re hiking along a mountain ridge trail (your constraint curve $g(x) = 0$). The mountain has contour lines—curves of constant elevation. At the highest point you can reach while staying on the trail, the trail is tangent to a contour line.

If the trail crossed through the contour line at an angle, you could walk along the trail to reach a higher contour. Only when the trail is perfectly tangent to a contour line—touching it at exactly one point—have you found the highest elevation on that trail.

Now here’s where the gradients come in:

The gradient $\nabla f$ points perpendicular to the contour lines (pointing uphill, in the direction of steepest ascent)
The gradient $\nabla g$ points perpendicular to the constraint curve (the trail)
When the contour line and trail are tangent to each other, their perpendiculars must be parallel

This is why $\nabla f = \lambda \nabla g$ at the optimum.

The Lagrangian Method

For a problem with multiple constraints:

$$x^* = \arg\min_{x \in \mathbb{R}^n} f(x) \quad \text{s.t.} \quad g_i(x) = 0, \quad i = 1, …, m$$

We construct the Lagrangian function: $$\mathcal{L}(x, \lambda) = f(x) + \sum_{i=1}^m \lambda_i g_i(x)$$

The optimality conditions are simply: $$\nabla_x \mathcal{L} = 0 \quad \text{(gradient condition)}, \quad \nabla_\lambda \mathcal{L} = 0 \quad \text{(recovers constraints)}$$

This transforms our constrained problem into solving a system of equations—much more tractable computationally.

Inequality Constraints: Karush-Kuhn-Tucker Conditions

When constraints can be inequalities:

$$x^* = \arg\min_{x \in \mathbb{R}^n} f(x) \quad \text{s.t.} \quad g_i(x) = 0, \quad h_j(x) \leq 0$$

We augment the Lagrangian: $$\mathcal{L}(x, \lambda, \mu) = f(x) + \sum_i \lambda_i g_i(x) + \sum_j \mu_j h_j(x)$$

The key difference with inequalities:

An inequality constraint $h(x) \leq 0$ might be:

Active: $h(x) = 0$ (you’re at the limit)
Inactive: $h(x) < 0$ (you have slack)

The KKT conditions encode this:

Complementary slackness: $\mu_j h_j(x) = 0$
- If constraint is active ($h_j = 0$), then $\mu_j$ can be positive
- If constraint has slack ($h_j < 0$), then $\mu_j = 0$
- Result: At the optimum, $\sum_j \mu_j h_j(x^*) = 0$ always!
- But the individual $\mu_j$ values tell us which constraints matter
Non-negativity: $\mu_j \geq 0$
- Tightening a constraint can only make things worse (or same), never better
- So the sensitivity must be non-negative

Intuition: The Lagrangian “penalizes” constraint violations during optimization. At the solution, complementary slackness ensures no penalty remains, but the multipliers $\mu_j$ reveal which constraints are binding (positive $\mu$) versus irrelevant (zero $\mu$).

Example: Minimize cost subject to production $\geq$ demand:

If you’re producing exactly the demand: constraint active, $\mu > 0$
If you’re overproducing anyway: constraint inactive, $\mu = 0$

Economic Interpretation

The multiplier $\lambda_i$ has a concrete meaning (for a constraint $g_i(x) = b_i$): $$\lambda_i = -\frac{\partial f^*}{\partial b_i}$$

It’s the shadow price—how much the optimal value improves per unit relaxation of constraint $i$. In economics:

Budget constraint: $\lambda$ = marginal utility of income
Resource constraint: $\lambda$ = marginal value of that resource
Production quota: $\lambda$ = marginal cost of production

This duality between optimization and valuation is fundamental to economic theory and pricing mechanisms.

Practical Implementation

Modern solvers handle the mechanics, but understanding helps diagnose issues:

Multiple solutions: Check second-order conditions (bordered Hessian)
No solution: Constraints may be inconsistent
Slow convergence: Problem may be ill-conditioned
Large multipliers: Constraint is very restrictive, consider reformulation

Note on Sufficiency

The KKT conditions are necessary for optimality. For sufficiency, check:

Convex objective function (for minimization)
Convex feasible region
Or verify second-order conditions (bordered Hessian)

Implementation Example

Let’s solve a concrete problem: minimize the distance from the origin to a line.

Problem: $\min_{x,y} f(x,y) = x^2 + y^2$ subject to $g(x,y) = x + y - 1 = 0$

Physical interpretation: Find the point on the line $x + y = 1$ closest to the origin.

We’ll solve this using the Lagrangian approach and verify our geometric intuition.

# Define symbolic variables
x, y, lam = sp.symbols('x y lambda')

# Objective function and constraint
f = x**2 + y**2
g = x + y - 1

# Construct Lagrangian: L(x, y, λ) = f(x, y) + λ·g(x, y)
L = f + lam * g

print("Objective function f(x,y):", f)
print("Constraint g(x,y) = 0:", g, "= 0")
print("\nLagrangian L(x,y,λ):")
display(Math(r'\mathcal{L}(x,y,\lambda) = ' + sp.latex(L)))

# Compute first-order conditions: ∇L = 0
dL_dx = sp.diff(L, x)
dL_dy = sp.diff(L, y)
dL_dlam = sp.diff(L, lam)

print("\nFirst-order conditions (∇L = 0):")
print(f"∂L/∂x = {dL_dx} = 0")
print(f"∂L/∂y = {dL_dy} = 0")
print(f"∂L/∂λ = {dL_dlam} = 0")

# Solve system of equations
solution = sp.solve([dL_dx, dL_dy, dL_dlam], (x, y, lam), dict=True)
print("\nSolution:")
for sol in solution:
    print(f"x* = {sol[x]}, y* = {sol[y]}, λ* = {sol[lam]}")
    print(f"f(x*, y*) = {float(f.subs(sol)):.6f}")
    print(f"\nVerification:")
    print(f"  Constraint satisfied: g(x*, y*) = {float(g.subs(sol)):.6f}")
    print(f"  Distance from origin: {float(sp.sqrt(f.subs(sol))):.6f}")
    
    # Check that gradients are parallel
    grad_f = [sp.diff(f, x).subs(sol), sp.diff(f, y).subs(sol)]
    grad_g = [sp.diff(g, x).subs(sol), sp.diff(g, y).subs(sol)]
    grad_f_float = [float(v) for v in grad_f]
    grad_g_float = [float(v) for v in grad_g]
    lambda_grad_g = [float(sol[lam] * v) for v in grad_g]
    
    print(f"  ∇f = {grad_f_float}")
    print(f"  ∇g = {grad_g_float}")
    print(f"  λ·∇g = {lambda_grad_g}")
    print(f"  Gradients parallel? {np.allclose(grad_f_float, lambda_grad_g)}")

Objective function f(x,y): x**2 + y**2
Constraint g(x,y) = 0: x + y - 1 = 0

Lagrangian L(x,y,λ):

$\displaystyle \mathcal{L}(x,y,\lambda) = \lambda \left(x + y - 1\right) + x^{2} + y^{2}$

First-order conditions (∇L = 0):
∂L/∂x = lambda + 2*x = 0
∂L/∂y = lambda + 2*y = 0
∂L/∂λ = x + y - 1 = 0

Solution:
x* = 1/2, y* = 1/2, λ* = -1
f(x*, y*) = 0.500000

Verification:
  Constraint satisfied: g(x*, y*) = 0.000000
  Distance from origin: 0.707107
  ∇f = [1.0, 1.0]
  ∇g = [1.0, 1.0]
  λ·∇g = [-1.0, -1.0]
  Gradients parallel? False

# Visualization of Lagrange multiplier solution
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Create meshgrid
x_range = np.linspace(-0.5, 2, 100)
y_range = np.linspace(-0.5, 2, 100)
X, Y = np.meshgrid(x_range, y_range)
Z = X**2 + Y**2

# Left: Contour plot with constraint
contours = ax1.contour(X, Y, Z, levels=15, cmap='viridis', linewidths=0.8, alpha=0.7)
ax1.clabel(contours, inline=True, fontsize=8, fmt='%g')

# Constraint line: x + y = 1
x_line = np.linspace(-0.5, 2, 100)
y_line = 1 - x_line
ax1.plot(x_line, y_line, 'r-', linewidth=2.5, label='Constraint: $x+y=1$')

# Optimal point
x_opt, y_opt = 0.5, 0.5
ax1.plot(x_opt, y_opt, 'ro', markersize=10, label='Optimum', zorder=5)

# Gradient vectors
grad_f = np.array([2*x_opt, 2*y_opt])
grad_g = np.array([1, 1])
scale = 0.2

ax1.arrow(x_opt, y_opt, scale*grad_f[0], scale*grad_f[1], 
          head_width=0.06, head_length=0.04, fc='blue', ec='blue', linewidth=1.5)
ax1.text(x_opt + scale*grad_f[0] + 0.08, y_opt + scale*grad_f[1], 
         r'$\nabla f$', fontsize=11, color='blue')

ax1.arrow(x_opt, y_opt, scale*grad_g[0], scale*grad_g[1], 
          head_width=0.06, head_length=0.04, fc='orange', ec='orange', linewidth=1.5)
ax1.text(x_opt + scale*grad_g[0] - 0.15, y_opt + scale*grad_g[1] + 0.05, 
         r'$\nabla g$', fontsize=11, color='orange')

ax1.plot([0, x_opt], [0, y_opt], 'k--', linewidth=1, alpha=0.4)
ax1.plot(0, 0, 'k^', markersize=6)

ax1.set_xlabel('$x$', fontsize=11)
ax1.set_ylabel('$y$', fontsize=11)
ax1.set_title('Constrained Optimization', fontsize=12)
ax1.legend(frameon=False)
ax1.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)
ax1.set_aspect('equal')
ax1.set_xlim(-0.3, 1.8)
ax1.set_ylim(-0.3, 1.8)

# Right: 3D view
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.6, edgecolor='none', antialiased=True)

# Constraint in 3D
z_line = x_line**2 + y_line**2
ax2.plot(x_line, y_line, z_line, 'r-', linewidth=3)

# Optimal point
z_opt = x_opt**2 + y_opt**2
ax2.scatter([x_opt], [y_opt], [z_opt], color='red', s=80, marker='o')
ax2.plot([x_opt, x_opt], [y_opt, y_opt], [0, z_opt], 'r--', linewidth=1.5, alpha=0.6)

ax2.set_xlabel('$x$', fontsize=10)
ax2.set_ylabel('$y$', fontsize=10)
ax2.set_zlabel('$f(x,y)$', fontsize=10)
ax2.set_title('3D View', fontsize=12)
ax2.view_init(elev=20, azim=45)

plt.tight_layout()
plt.show()

print("\nKey Insight:")
print("• The gradients ∇f and ∇g are parallel at the optimum")
print("• The optimal point is where the smallest circle touches the line")

png

Key Insight:
• The gradients ∇f and ∇g are parallel at the optimum
• The optimal point is where the smallest circle touches the line

Newton’s Method

While gradient descent follows the slope downhill, Newton’s method uses second-order information (curvature) to take smarter steps toward the minimum.

The Core Idea

Formally, at each point $x_k$, we approximate $f$ with a quadratic (second-order Taylor expansion):

$$f(x) \approx f(x_k) + \nabla f(x_k)^T (x - x_k) + \frac{1}{2}(x - x_k)^T H(x_k) (x - x_k)$$

The minimum of this quadratic occurs where its gradient equals zero. Solving for $x$ gives the Newton update:

$$x_{k+1} = x_k - H(x_k)^{-1} \nabla f(x_k)$$

where $H(x_k)$ is the Hessian matrix at $x_k$.

Intuition: Gradient Descent vs Newton

Consider a simple 1D function near a minimum:

Gradient Descent: Takes steps proportional to the slope: $x_{k+1} = x_k - \alpha f’(x_k)$
- In a narrow valley: overshoots, oscillates, zigzags
- Needs many small steps
Newton’s Method: Accounts for curvature: $x_{k+1} = x_k - \frac{f’(x_k)}{f’'(x_k)}$
- In a narrow valley: recognizes the curvature, takes appropriately sized step
- Can reach the minimum in one step (for quadratic functions!)

The key difference: Dividing by $f’'$ (or $H^{-1}$ in higher dimensions) rescales the step based on local curvature. Steep curvature → smaller steps. Flat curvature → larger steps.

Convergence Properties

Quadratic Convergence: Near the minimum, the number of correct digits roughly doubles each iteration. If you’re 0.01 away, next iteration you’re 0.0001 away, then 0.00000001 away.

This is dramatically faster than gradient descent’s linear convergence (constant factor reduction each step).

Requirements for convergence:

Starting point $x_0$ must be sufficiently close to the minimum
Hessian must be positive definite (convex region)
Hessian must be non-singular (invertible)

Computational Cost

The trade-off: Newton’s method needs:

$\mathcal{O}(n^2)$ storage for the Hessian
$\mathcal{O}(n^3)$ operations to solve the linear system $H \Delta x = -\nabla f$
Second derivatives (analytically or via finite differences)

Gradient descent needs only:

$\mathcal{O}(n)$ storage for the gradient
$\mathcal{O}(n)$ operations per iteration
First derivatives

When to use Newton:

Small to medium problems (n < 1000): Newton’s fast convergence outweighs cost
High accuracy needed: Quadratic convergence reaches precision quickly
Hessian available: Analytical derivatives or structured problems
Good initial guess: Near the solution already

When NOT to use Newton:

Large-scale problems (n > 10,000): Memory and computation prohibitive
Far from solution: May diverge or take poor steps
Non-convex regions: Hessian may not be positive definite

def numerical_hessian(f, x, h=1e-5):
    # Compute numerical Hessian of f at point x
    n = len(x)
    H = np.zeros((n, n))
    
    for i in range(n):
        for j in range(n):
            ei, ej = np.zeros(n), np.zeros(n)
            ei[i] = h
            ej[j] = h
            
            if i == j:
                H[i,j] = (f(x + ei) - 2*f(x) + f(x - ei)) / h**2
            else:
                H[i,j] = (f(x + ei + ej) - f(x + ei) - f(x + ej) + f(x)) / h**2
    
    return H

def newton_method(f, x0, max_iter=20, tol=1e-6):
    """Newton's method for optimization."""
    x = x0.copy()
    history = [x.copy()]
    grad_norms = []
    
    for i in range(max_iter):
        grad = numerical_gradient(f, x)
        grad_norm = np.linalg.norm(grad)
        grad_norms.append(grad_norm)
        
        if grad_norm < tol:
            print(f"Converged in {i} iterations (||grad|| = {grad_norm:.2e})")
            break
        
        H = numerical_hessian(f, x)
        
        try:
            delta = np.linalg.solve(H, -grad)
        except np.linalg.LinAlgError:
            print(f"Singular Hessian at iteration {i}")
            break
        
        x += delta
        history.append(x.copy())
    
    return x, history, grad_norms

# Compare Newton's Method with Gradient Descent
print("Rosenbrock function: f(x,y) = (1-x)² + 100(y-x²)²")

# Rosenbrock function
rosenbrock = lambda p: (1 - p[0])**2 + 100*(p[1] - p[0]**2)**2

x0 = np.array([0.5, 0.5])
print(f"Starting point: {x0}\n")

# Run both methods
print("Gradient Descent:")
gd_result, gd_hist = gradient_descent(rosenbrock, x0, alpha=0.001, max_iter=50, tol=1e-6)
print(f"Final: {np.round(gd_result, 4)}, f = {rosenbrock(gd_result):.2e}\n")

print("Newton's Method:")
newton_result, newton_hist, newton_grads = newton_method(rosenbrock, x0, max_iter=50, tol=1e-6)
print(f"Final: {np.round(newton_result, 4)}, f = {rosenbrock(newton_result):.2e}\n")

# Visualization
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4.5))

# Create meshgrid
x_vals = np.linspace(-0.5, 1.5, 200)
y_vals = np.linspace(-0.5, 1.5, 200)
X, Y = np.meshgrid(x_vals, y_vals)
Z = (1 - X)**2 + 100*(Y - X**2)**2
Z_log = np.log10(Z + 1)

# Gradient Descent
contours = ax1.contour(X, Y, Z_log, levels=20, cmap='viridis', linewidths=0.8, alpha=0.7)
gd_hist = np.array(gd_hist)
ax1.plot(gd_hist[:, 0], gd_hist[:, 1], 'r-', linewidth=1.5, alpha=0.8)
ax1.plot(gd_hist[0, 0], gd_hist[0, 1], 'ro', markersize=6, label='Start')
ax1.plot(gd_hist[-1, 0], gd_hist[-1, 1], 'rs', markersize=6, label='End')
ax1.plot(1, 1, 'r^', markersize=7, label='True min')

ax1.set_xlabel('$x$', fontsize=11)
ax1.set_ylabel('$y$', fontsize=11)
ax1.set_title('Gradient Descent', fontsize=12)
ax1.legend(frameon=False, fontsize=9)
ax1.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Newton's Method
ax2.contour(X, Y, Z_log, levels=20, cmap='viridis', linewidths=0.8, alpha=0.7)
newton_hist = np.array(newton_hist)
ax2.plot(newton_hist[:, 0], newton_hist[:, 1], 'b-', linewidth=1.5, alpha=0.8)
ax2.plot(newton_hist[0, 0], newton_hist[0, 1], 'bo', markersize=6, label='Start')
ax2.plot(newton_hist[-1, 0], newton_hist[-1, 1], 'bs', markersize=6, label='End')
ax2.plot(1, 1, 'b^', markersize=7, label='True min')

ax2.set_xlabel('$x$', fontsize=11)
ax2.set_ylabel('$y$', fontsize=11)
ax2.set_title('Newton Method', fontsize=12)
ax2.legend(frameon=False, fontsize=9)
ax2.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Convergence comparison
gd_values = [rosenbrock(p) for p in gd_hist]
newton_values = [rosenbrock(p) for p in newton_hist]

ax3.semilogy(gd_values, 'r-', linewidth=1.5, alpha=0.7, label='Gradient Descent')
ax3.semilogy(newton_values, 'b--', linewidth=1.5, alpha=0.7, label='Newton')

ax3.set_xlabel('Iteration', fontsize=11)
ax3.set_ylabel('$f(x)$ (log scale)', fontsize=11)
ax3.set_title('Convergence', fontsize=12)
ax3.legend(frameon=False, fontsize=10)
ax3.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

plt.tight_layout()
plt.show()

Rosenbrock function: f(x,y) = (1-x)² + 100(y-x²)²
Starting point: [0.5 0.5]

Gradient Descent:
Final: [0.6377 0.4048], f = 1.32e-01

Newton's Method:
Converged in 6 iterations (||grad|| = 2.46e-09)
Final: [1. 1.], f = 4.01e-12

png

Quasi-Newton Methods: BFGS

Newton’s method is powerful but expensive. Quasi-Newton methods offer a clever compromise: approximate the Hessian using only gradient information, building up curvature knowledge as you optimize.

The Core Problem

Newton’s method requires computing and inverting the Hessian at each step:

Computing $H$: $\mathcal{O}(n^2)$ second derivatives
Inverting $H$: $\mathcal{O}(n^3)$ operations

Can we get Newton-like convergence without explicitly computing $H$?

The BFGS Idea

Instead of computing $H^{-1}$, we:

Start with an initial guess: $B_0 = I$ (identity matrix)
At each step, observe how the gradient changes
Update our approximation $B_k \approx H^{-1}$ to match this observation

Key insight: If we move from $x_k$ to $x_{k+1}$, we observe:

Position change: $s_k = x_{k+1} - x_k$
Gradient change: $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$

For a quadratic function, $H$ relates these: $y_k = H s_k$, or equivalently $s_k = H^{-1} y_k$.

BFGS constructs $B_{k+1}$ to satisfy this secant equation while staying close to $B_k$ and maintaining symmetry and positive definiteness.

The BFGS Update Formula

The update that achieves all these properties is:

$$B_{k+1} = \left(I - \rho_k s_k y_k^T\right) B_k \left(I - \rho_k y_k s_k^T\right) + \rho_k s_k s_k^T$$

where $\rho_k = \dfrac{1}{y_k^T s_k}$.

Equivalent expanded form (used in implementation):

$$B_{k+1} = B_k + \left(1 + \rho_k y_k^T B_k y_k\right) \rho_k s_k s_k^T - \rho_k (B_k y_k s_k^T + s_k y_k^T B_k)$$

where:

$s_k = x_{k+1} - x_k$ (step taken)
$y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$ (gradient change)

Intuition:

The update maintains symmetry and positive definiteness
It incorporates new curvature information from the gradient change
The secant equation $B_{k+1} y_k = s_k$ is satisfied
Together, these properties approximate the inverse Hessian

The algorithm then takes a Newton-like step: $x_{k+1} = x_k - \alpha_k B_k \nabla f(x_k)$

where $\alpha_k$ is found via line search (ensuring sufficient decrease).

Why BFGS Works

The secant equation $B y = s$ is a first-order approximation to $H^{-1}$. By enforcing this at each step, BFGS accumulates curvature information from the gradient history.

After $n$ steps in $n$ dimensions, BFGS has enough information to reconstruct $H^{-1}$ exactly for a quadratic function. This gives it superlinear convergence—faster than linear (gradient descent) but not quite quadratic (Newton).

Practical Advantages

No Hessian needed: Only requires gradients (first derivatives)
Efficient: $\mathcal{O}(n^2)$ storage, $\mathcal{O}(n^2)$ operations per iteration
Robust: Maintains positive definiteness automatically
Fast convergence: Superlinear near the minimum
Line search compatible: Step size $\alpha_k$ can be tuned for stability

The Line Search Component

BFGS typically uses a line search to find $\alpha_k$ satisfying the Wolfe conditions:

Sufficient decrease: $f(x_k + \alpha_k d_k) \leq f(x_k) + c_1 \alpha_k \nabla f(x_k)^T d_k$
- Ensures we make progress
Curvature condition: $|\nabla f(x_k + \alpha_k d_k)^T d_k| \leq c_2 |\nabla f(x_k)^T d_k|$
- Ensures we don’t stop too early (allows BFGS update to work properly)

where $d_k = -B_k \nabla f(x_k)$ is the search direction.

When to Use BFGS

Use BFGS for:

Medium-scale problems (100 to 10,000 variables)
When Hessian is unavailable or expensive
When gradient descent is too slow
When Newton’s method is too expensive

BFGS vs alternatives:

vs Gradient Descent: Much faster convergence, worth the extra cost per iteration
vs Newton: Cheaper per iteration, no need for second derivatives, nearly as fast
vs L-BFGS: Full BFGS for medium problems; L-BFGS for very large (>10,000 variables)

def line_search(f, x, grad, direction, alpha=1.0, rho=0.5, c=1e-4):
    # Backtracking line search with Armijo condition
    f_x = f(x)
    grad_dot_dir = np.dot(grad, direction)
    
    while f(x + alpha * direction) > f_x + c * alpha * grad_dot_dir:
        alpha *= rho
        if alpha < 1e-10:
            break
    
    return alpha

def bfgs_method(f, x0, max_iter=100, tol=1e-6):
    # BFGS optimization method
    n = len(x0)
    x = x0.copy()
    B = np.eye(n)  # Initial inverse Hessian approximation
    history = [x.copy()]
    grad_norms = []
    
    for i in range(max_iter):
        grad = numerical_gradient(f, x)
        grad_norm = np.linalg.norm(grad)
        grad_norms.append(grad_norm)
        
        if grad_norm < tol:
            print(f"Converged in {i} iterations (||grad|| = {grad_norm:.2e})")
            break
        
        # Search direction
        direction = -B @ grad
        
        # Line search
        alpha = line_search(f, x, grad, direction)
        
        # Update position
        s = alpha * direction
        x_new = x + s
        grad_new = numerical_gradient(f, x_new)
        y = grad_new - grad
        
        # BFGS update
        y_dot_s = np.dot(y, s)
        if y_dot_s > 1e-10:
            rho = 1.0 / y_dot_s
            I = np.eye(n)
            B_y = B @ y
            B = B + (1 + rho * np.dot(y, B_y)) * rho * np.outer(s, s) \
                  - rho * (np.outer(B_y, s) + np.outer(s, B_y))
        
        x = x_new
        history.append(x.copy())
    
    return x, history, grad_norms

Method Comparison

We compare three fundamental optimization algorithms:

Gradient Descent: First-order, cheap per iteration, slow convergence
Newton’s Method: Second-order, expensive per iteration, fast convergence
BFGS: Quasi-Newton, moderate cost, fast convergence

The comparison reveals the fundamental trade-off: cost per iteration vs. number of iterations needed.

# Compare all three methods
print("Comparison: Gradient Descent vs Newton vs BFGS")
print("Rosenbrock function: f(x,y) = (1-x)² + 100(y-x²)²\n")

x0 = np.array([0.5, 0.5])

# Run all methods
print("Gradient Descent:")
gd_result, gd_hist = gradient_descent(rosenbrock, x0, alpha=0.001, max_iter=100, tol=1e-6)
print(f"Final: {np.round(gd_result, 4)}, iterations: {len(gd_hist)-1}\n")

print("Newton's Method:")
newton_result, newton_hist, newton_grads = newton_method(rosenbrock, x0, max_iter=100, tol=1e-6)
print(f"Final: {np.round(newton_result, 4)}, iterations: {len(newton_hist)-1}\n")

print("BFGS:")
bfgs_result, bfgs_hist, bfgs_grads = bfgs_method(rosenbrock, x0, max_iter=100, tol=1e-6)
print(f"Final: {np.round(bfgs_result, 4)}, iterations: {len(bfgs_hist)-1}\n")

# Visualization
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 11))

# Create meshgrid
x_vals = np.linspace(-0.2, 1.3, 200)
y_vals = np.linspace(-0.2, 1.3, 200)
X, Y = np.meshgrid(x_vals, y_vals)
Z = (1 - X)**2 + 100*(Y - X**2)**2
Z_log = np.log10(Z + 1)

gd_hist = np.array(gd_hist)
newton_hist = np.array(newton_hist)
bfgs_hist = np.array(bfgs_hist)

# Plot 1: Gradient Descent
ax1.contour(X, Y, Z_log, levels=20, cmap='viridis', linewidths=0.8, alpha=0.7)
ax1.plot(gd_hist[:, 0], gd_hist[:, 1], 'r-', linewidth=1.5, alpha=0.8)
ax1.plot(gd_hist[0, 0], gd_hist[0, 1], 'ro', markersize=6)
ax1.plot(1, 1, 'r^', markersize=7)
ax1.set_xlabel('$x$', fontsize=11)
ax1.set_ylabel('$y$', fontsize=11)
ax1.set_title('Gradient Descent', fontsize=12)
ax1.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Plot 2: Newton
ax2.contour(X, Y, Z_log, levels=20, cmap='viridis', linewidths=0.8, alpha=0.7)
ax2.plot(newton_hist[:, 0], newton_hist[:, 1], 'b-', linewidth=1.5, alpha=0.8)
ax2.plot(newton_hist[0, 0], newton_hist[0, 1], 'bo', markersize=6)
ax2.plot(1, 1, 'b^', markersize=7)
ax2.set_xlabel('$x$', fontsize=11)
ax2.set_ylabel('$y$', fontsize=11)
ax2.set_title('Newton', fontsize=12)
ax2.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Plot 3: BFGS
ax3.contour(X, Y, Z_log, levels=20, cmap='viridis', linewidths=0.8, alpha=0.7)
ax3.plot(bfgs_hist[:, 0], bfgs_hist[:, 1], 'g-', linewidth=1.5, alpha=0.8)
ax3.plot(bfgs_hist[0, 0], bfgs_hist[0, 1], 'go', markersize=6)
ax3.plot(1, 1, 'g^', markersize=7)
ax3.set_xlabel('$x$', fontsize=11)
ax3.set_ylabel('$y$', fontsize=11)
ax3.set_title('BFGS', fontsize=12)
ax3.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

# Plot 4: Gradient norm comparison
ax4.semilogy(newton_grads, 'b--', linewidth=2, alpha=0.8, label='Newton', marker='o', markersize=4)
ax4.semilogy(bfgs_grads, 'g:', linewidth=2, alpha=0.8, label='BFGS', marker='s', markersize=4)
ax4.axhline(1e-6, color='red', linestyle='--', linewidth=1, alpha=0.5, label='Tolerance')

ax4.set_xlabel('Iteration', fontsize=11)
ax4.set_ylabel('$||\u2207f||$ (log scale)', fontsize=11)
ax4.set_title('Gradient Norm', fontsize=12)
ax4.legend(frameon=False, fontsize=10)
ax4.grid(True, alpha=0.2, linestyle='--', linewidth=0.5)

plt.tight_layout()
plt.show()

Comparison: Gradient Descent vs Newton vs BFGS
Rosenbrock function: f(x,y) = (1-x)² + 100(y-x²)²

Gradient Descent:
Final: [0.651  0.4221], iterations: 100

Newton's Method:
Converged in 6 iterations (||grad|| = 2.46e-09)
Final: [1. 1.], iterations: 6

BFGS:
Final: [1. 1.], iterations: 100

png