6.1 Numerical Optimization

Why this week matters¶

Many statistical procedures can be written as optimization problems.

Examples include ordinary least squares, maximum likelihood estimation, logistic regression, regularized regression, neural-network training, nonlinear curve fitting, and numerical equation solving.

The central question for Week 12 is:

How do we turn a statistical estimation problem into an objective function that a numerical algorithm can solve?

This week focuses on root finding, bisection, Newton--Raphson, objective and loss functions, gradients and Hessians, gradient descent, learning rates, convergence criteria, ordinary least squares as optimization, maximum likelihood as optimization, and validation against SciPy.

1. Root finding vs. optimization¶

Two common numerical tasks are:

Root finding:
find x such that f(x) = 0

Optimization:
find x* that minimizes f(x)

These are related because a differentiable optimum often satisfies:

f'(x*) = 0

2. A root-finding example¶

Suppose:

f(x) = x² - 2

We want:

x² - 2 = 0

Define:

def f(x):
    return x ** 2 - 2

The positive solution is sqrt(2), but we will approximate it numerically.

3. Bracketing¶

If f(a) and f(b) have opposite signs and f is continuous, at least one root lies between a and b.

For example:

f(1)
f(2)

give opposite signs.

This provides the bracket:

[1, 2]

4. Bisection¶

Bisection repeatedly halves a valid bracket.

[a,b]
  ↓
m = (a+b)/2
  ↓
keep half containing sign change
  ↓
repeat

A simple implementation:

def bisection(
    f,
    a,
    b,
    tol=1e-8,
    max_iter=100
):
    if f(a) * f(b) > 0:
        raise ValueError(
            "Root is not bracketed."
        )

    for _ in range(max_iter):
        midpoint = (a + b) / 2
        value = f(midpoint)

        if abs(value) < tol:
            return midpoint

        if f(a) * value < 0:
            b = midpoint
        else:
            a = midpoint

    return midpoint

Bisection is slow but robust when its assumptions hold.

5. Newton--Raphson¶

Newton’s method uses derivative information:

x_(k+1)
=
x_k - f(x_k)/f'(x_k)

Geometrically, the tangent line at the current point predicts where the function crosses zero.

6. Deriving Newton’s method¶

Use a first-order Taylor approximation:

f(x)
≈
f(x_k) + f'(x_k)(x-x_k)

Set this equal to zero and solve for x:

x
=
x_k - f(x_k)/f'(x_k)

This becomes the next iterate.

7. Newton--Raphson in Python¶

def newton(
    f,
    f_prime,
    x0,
    tol=1e-8,
    max_iter=100
):
    x = x0

    for _ in range(max_iter):
        derivative = f_prime(x)

        if abs(derivative) < 1e-12:
            raise ValueError(
                "Derivative is too close to zero."
            )

        x_new = (
            x
            - f(x) / derivative
        )

        if abs(x_new - x) < tol:
            return x_new

        x = x_new

    return x

8. Why Newton can be fast¶

Near a well-behaved simple root, Newton’s method can converge approximately quadratically.

But it may fail because of poor starting values, near-zero derivatives, oscillation, divergence, multiple roots, or invalid parameter regions.

Fast convergence does not imply guaranteed convergence.

9. Convergence criteria¶

Possible stopping rules include:

|x_(k+1) - x_k| < tolerance

and:

|f(x_k)| < tolerance

A robust implementation should also have a maximum number of iterations.

10. Validate with SciPy¶

SciPy provides:

from scipy.optimize import root_scalar

For bisection:

result = root_scalar(
    f,
    bracket=[1, 2],
    method="bisect"
)

For Newton:

result = root_scalar(
    f,
    x0=1.5,
    fprime=f_prime,
    method="newton"
)

Inspect:

result.root
result.converged
result.iterations

11. Optimization objectives¶

An optimization problem has the form:

minimize f(theta)

The function may be called:

objective
loss
cost
negative log-likelihood
residual sum of squares

The optimizer searches parameter space for values that reduce the objective.

12. Gradient¶

For:

f(theta_1,...,theta_p)

the gradient is:

∇f(theta)
=
[
∂f/∂theta_1,
...,
∂f/∂theta_p
]^T

The gradient points in the direction of steepest local increase.

Therefore -∇f points toward local decrease.

13. Gradient descent¶

The update is:

theta_(k+1)
=
theta_k - alpha ∇f(theta_k)

where alpha is the learning rate.

Conceptually:

current parameters
      ↓
calculate gradient
      ↓
move opposite gradient
      ↓
new parameters
      ↓
repeat

14. Learning rate¶

If the learning rate is too small:

convergence is slow

If it is too large:

overshooting
oscillation
divergence

The step size is part of the numerical method.

15. Hessian¶

For a multivariable objective, the Hessian is the matrix of second derivatives.

It describes local curvature and appears in Newton optimization and uncertainty approximations.

For multiple dimensions, Newton’s method can be written:

theta_(k+1)
=
theta_k - H(theta_k)^(-1) ∇f(theta_k)

In practice, solve the linear system instead of explicitly forming an inverse.

16. OLS as optimization¶

Linear regression:

y = X beta + epsilon

Ordinary least squares minimizes:

RSS(beta)
=
||y - X beta||²

Example loss:

def mse_loss(
    beta,
    X,
    y
):
    residuals = (
        y - X @ beta
    )

    return np.mean(
        residuals ** 2
    )

17. Gradient of OLS loss¶

For:

L(beta)
=
(1/n)||y-Xbeta||²

the gradient is:

∇L(beta)
=
(2/n)X^T(Xbeta-y)

In NumPy:

def mse_gradient(
    beta,
    X,
    y
):
    n = len(y)

    return (
        2 / n
        * X.T
        @ (
            X @ beta - y
        )
    )

18. Gradient descent for regression¶

Initialize:

beta = np.zeros(
    X.shape[1]
)

Then repeatedly update:

beta = (
    beta
    - learning_rate
    * mse_gradient(
        beta,
        X,
        y
    )
)

Compare the result with:

np.linalg.lstsq(
    X,
    y,
    rcond=None
)

Two correct computational routes should agree closely.

19. Feature scaling¶

Gradient descent can behave poorly when predictor scales differ dramatically.

Standardization:

means = X.mean(axis=0)
sds = X.std(axis=0)

Z = (
    X - means
) / sds

can make the optimization geometry much easier.

20. Maximum likelihood estimation¶

For density:

f(x | theta)

the likelihood is:

L(theta)
=
product f(x_i | theta)

The log-likelihood is:

ell(theta)
=
sum log f(x_i | theta)

Maximum likelihood chooses:

theta_hat
=
argmax ell(theta)

21. Negative log-likelihood¶

Most optimization libraries minimize.

Therefore:

maximize log-likelihood

is equivalent to:

minimize negative log-likelihood

This is one of the central patterns in statistical computing.

22. Normal mean MLE¶

Suppose:

X_i ~ Normal(mu, sigma²)

with known sigma.

Ignoring constants, minimizing the negative log-likelihood in mu is equivalent to minimizing:

sum (x_i - mu)²

The numerical optimizer should recover approximately the sample mean.

23. Constrained parameters¶

If sigma is unknown, it must satisfy:

sigma > 0

One useful parameterization is:

eta = log(sigma)

so:

sigma = exp(eta)

The optimizer can search over any real eta while always producing a valid positive sigma.

24. SciPy `minimize`¶

Typical use:

from scipy.optimize import minimize

result = minimize(
    objective,
    x0=initial_parameters
)

Always inspect:

result.x
result.fun
result.success
result.message
result.nit

Do not trust result.x without checking convergence.

25. Local vs. global minima¶

A nonconvex objective can have several local minima.

Different initial values may lead to different solutions.

Possible strategies include multiple starting values, better initialization, visualization, or global optimization methods.

26. Convexity¶

For a convex objective:

every local minimum is a global minimum

OLS loss is convex.

This makes least-squares optimization easier to reason about than many complex machine-learning objectives.

27. Numerical differentiation¶

If an analytical derivative is inconvenient, approximate it:

f'(x)
≈
[f(x+h)-f(x-h)]/(2h)

Implementation:

def numerical_derivative(
    f,
    x,
    h=1e-5
):
    return (
        f(x + h)
        - f(x - h)
    ) / (
        2 * h
    )

28. Finite-difference step size¶

If h is too large, approximation error can be substantial.

If h is too small, floating-point cancellation can dominate.

This illustrates a common theme in numerical computing:

More precision in an input setting does not always produce more accurate computation.

29. Gradient checking¶

Analytical gradients should be checked against numerical approximations.

This is extremely useful when implementing statistical or machine-learning algorithms from scratch.

A wrong gradient may still produce code that runs.

30. Optimization diagnostics¶

Always inspect:

objective history
gradient magnitude
parameter changes
iteration count
convergence flag
final objective
sensitivity to starting values

An optimizer returning numbers does not prove it solved the intended problem.

31. Common mistakes¶

maximizing when the library minimizes;
incorrect derivative signs;
learning rate too large;
no maximum iteration limit;
ignored parameter constraints;
ignored convergence messages;
explicit matrix inverses when a solve is preferable;
failure to check sensitivity to starting values.

32. Statistical optimization workflow¶

Statistical model
      ↓
Define equation / objective
      ↓
Derive or approximate derivatives
      ↓
Choose starting values
      ↓
Choose numerical method
      ↓
Iterate
      ↓
Check convergence
      ↓
Validate
      ↓
Interpret

33. Key ideas¶

By the end of Week 12, you should be able to explain:

Root finding vs. optimization.
Bracketing and bisection.
Newton--Raphson updates.
Why Newton’s method can converge quickly or fail.
Why starting values matter.
Objective and loss functions.
Gradients and Hessians.
Gradient descent.
Learning-rate behavior.
OLS as optimization.
The OLS gradient.
Why feature scaling matters.
MLE as optimization.
Negative log-likelihood.
Constrained parameter transformations.
Local vs. global minima.
Numerical differentiation.
Gradient checking.
Convergence diagnostics.
Validation against trusted routines.

34. Recommended reading¶

35. YouTube recommendations¶

1. MathAndPhysics --- Newton-Raphson Algorithm in Python¶

A focused implementation of Newton--Raphson in Python, including the update rule and practical considerations.

Recommended use: Watch before implementing Newton’s method from scratch.

Watch on YouTube

2. Dataquest --- Gradient Descent From Scratch In Python¶

A practical implementation of gradient descent for linear regression. It covers loss, partial derivatives, parameter updates, learning rates, and the training loop.

Recommended use: Watch alongside the OLS and gradient-descent sections.