Activity 6.1

This page is an ipython notebook, meaning you can download it and run the executable cells. It directly follows Activity 6.1. Reference answers appear after the questions, challenges, predictions, and reflections. They show one reasonable solution or interpretation, but other correct approaches may exist.

Root Finding, Newton--Raphson, Gradient Descent & Statistical Optimization¶

Tools: Python, NumPy, SciPy, matplotlib, Jupyter/VS Code

Learning objectives¶

By the end of this tutorial, you should be able to:

• formulate a root-finding problem;

• implement bisection;

• implement Newton--Raphson;

• diagnose convergence failures;

• define objective functions;

• compute analytical and numerical derivatives;

• implement gradient descent;

• explore learning-rate behavior;

• fit OLS via optimization;

• formulate a negative log-likelihood;

• estimate parameters using SciPy;

• validate custom implementations.

Part 0 --- Set up¶

Create:

notebooks/week12_optimization.ipynb

python import numpy as np import pandas as pd import matplotlib.pyplot as plt

from scipy.optimize import ( root_scalar, minimize )

python rng = np.random.default_rng(556)

Part 1 --- Define a root problem¶

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

python f(1) f(2)

Part 2 --- Visualize the root¶

python x_grid = np.linspace( 0, 3, 300 )

plt.plot( x_grid, f(x_grid) )

plt.axhline(0) plt.xlabel(“x”) plt.ylabel(“f(x)”) plt.show()

Reference answer¶

The endpoint values have opposite signs, so continuity brackets a root. For x**2 - 2, the exact positive root is sqrt(2).

Part 3 --- One bisection step¶

Start:

python a = 1.0 b = 2.0

midpoint = ( a + b ) / 2

python f(a) f(midpoint) f(b)

Part 4 --- Implement bisection¶

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

for iteration in range(
    max_iter
):
    midpoint = (
        a + b
    ) / 2

    value = f(midpoint)

    if abs(value) < tol:
        return (
            midpoint,
            iteration + 1
        )

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

return midpoint, max_iter

python np.sqrt(2)

Reference answer¶

Bisection keeps the half-interval whose endpoints still have opposite signs. Each iteration halves the bracket and guarantees convergence for a continuous sign-changing interval.

Part 5 --- Invalid brackets¶

Try:

python bisection( f, 2, 3 )

Part 6 --- Validate bisection with SciPy¶

python result = root_scalar( f, bracket=[1, 2], method=“bisect” )

python result.root result.converged result.iterations

Reference answer¶

Failing immediately when the bracket lacks a sign change prevents the algorithm from returning a meaningless result when its assumptions are violated.

Part 7 --- Newton derivative¶

Define:

python def f_prime(x): return 2 * x

python x = 1.5

Part 8 --- Implement Newton--Raphson¶

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

for iteration 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
    )

    history.append(
        x_new
    )

    if abs(x_new - x) < tol:
        return (
            x_new,
            iteration + 1,
            history
        )

    x = x_new

return x, max_iter, history

Part 9 --- Visualize Newton convergence¶

Plot the sequence of iterates against sqrt(2).

Compare the number of iterations with bisection.

Reference answer¶

Newton’s method often converges in fewer iterations than bisection near a well-behaved root, but it relies on derivatives and a suitable starting value.

Part 10 --- Validate Newton with SciPy¶

python result = root_scalar( f, x0=1.5, fprime=f_prime, method=“newton” )

Part 11 --- Poor starting values¶

Define:

python def cubic(x): return ( x ** 3 - 2 * x + 2 )

python def cubic_prime(x): return ( 3 * x ** 2 - 2 )

text x0 root/result iterations converged?

Part 12 --- Numerical derivatives¶

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

Part 13 --- Finite-difference step size¶

Try:

text h = 1e-1 1e-3 1e-5 1e-8 1e-12

Reference answer¶

The numerical derivative should be close to the analytic derivative. Too-large steps create truncation error; extremely small steps suffer floating-point cancellation, so smaller is not always better.

Part 14 --- Simple optimization¶

Define:

python def objective(x): return ( x - 3 ) ** 2

python def gradient(x): return ( 2 * ( x - 3 ) )

Reference answer¶

Optimization searches for parameter values minimizing/maximizing an objective. Step size must balance progress with stability; overly large steps can overshoot or diverge.

Part 15 --- One gradient-descent update¶

Start:

python x = 8.0 learning_rate = 0.1

Part 16 --- Implement gradient descent¶

python def gradient_descent( objective, gradient, x0, learning_rate=0.1, tol=1e-8, max_iter=1000 ): x = x0 history = [ ( x, objective(x) ) ]

for iteration in range(
    max_iter
):
    x_new = (
        x
        - learning_rate
        * gradient(x)
    )

    history.append(
        (
            x_new,
            objective(x_new)
        )
    )

    if abs(x_new - x) < tol:
        return (
            x_new,
            iteration + 1,
            history
        )

    x = x_new

return x, max_iter, history

Part 17 --- Learning-rate experiment¶

Run with:

text 0.01 0.1 0.5 1.1

text learning rate iterations final x final loss

Part 18 --- Multivariable gradient descent¶

Define:

python def objective_2d( theta ): return ( ( theta[0] - 2 ) ** 2 + 2 * ( theta[1] + 1 ) ** 2 )

python def gradient_2d( theta ): return np.array([ 2 * ( theta[0] - 2 ), 4 * ( theta[1] + 1 ) ])

python theta = np.array([ 8.0, 5.0 ])

text [2, -1]

Part 19 --- Simulate regression data¶

python n = 200

x = rng.normal( size=n )

X = np.column_stack([ np.ones(n), x ])

beta_true = np.array([ 5.0, 2.0 ])

y = ( X @ beta_true + rng.normal( 0, 1.5, size=n ) )

Part 20 --- MSE loss¶

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

return np.mean(
    residuals ** 2
)

Part 21 --- MSE gradient¶

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

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

Part 22 --- Numerical gradient checker¶

Write:

python def numerical_gradient( objective, theta, h=1e-5 ): grad = np.zeros_like( theta, dtype=float )

for j in range(
    len(theta)
):
    theta_plus = theta.copy()
    theta_minus = theta.copy()

    theta_plus[j] += h
    theta_minus[j] -= h

    grad[j] = (
        objective(theta_plus)
        - objective(theta_minus)
    ) / (
        2 * h
    )

return grad

Reference answer¶

A numerical gradient checker perturbs one parameter at a time, approximates derivatives with finite differences, and compares them with the analytic gradient.

Part 23 --- Fit OLS with gradient descent¶

Start:

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

Reference answer¶

Numerical and analytic gradients should agree within tolerance. np.allclose() is appropriate because finite differences and floating-point arithmetic are approximate.

Part 24 --- Compare with least squares¶

python beta_lstsq, *_ = ( np.linalg.lstsq( X, y, rcond=None ) )

python beta beta_lstsq beta_true

Reference answer¶

Gradient-descent OLS estimates should converge close to the np.linalg.lstsq() solution when learning rate and stopping conditions are appropriate.

Part 25 --- Visualize regression convergence¶

Plot:

text iteration vs. MSE

Reference answer¶

Convergence appears as a decreasing objective and stabilizing parameters/gradient near a minimum.

Part 26 --- Feature scaling¶

Add a predictor on a much larger scale:

python x2 = rng.normal( 0, 1000, size=n )

Reference answer¶

Feature scaling affects gradient-descent geometry; very different predictor scales can force a tiny learning rate and slow or destabilize convergence.

Part 27 --- Normal likelihood¶

Simulate:

python data = rng.normal( loc=7, scale=2, size=200 )

python def normal_nll_mu( mu, data, sigma ): residuals = ( data - mu )

return (
    len(data)
    * np.log(sigma)
    + np.sum(
        residuals ** 2
    ) / (
        2 * sigma ** 2
    )
)

Reference answer¶

Standardization usually makes a common learning rate work better and improves convergence without changing the underlying fitted relationship once coefficients are interpreted on the proper scale.

Part 28 --- Plot the negative log-likelihood¶

Evaluate the NLL over a grid of mu values.

Plot the curve.

Where does the minimum occur relative to the sample mean?

Reference answer¶

For a normal model, minimizing negative log-likelihood gives the maximum-likelihood estimate; the surface shows how the objective changes around plausible parameter values.

Part 29 --- Estimate the mean with SciPy¶

python result = minimize( lambda theta: normal_nll_mu( theta[0], data, sigma=2 ), x0=np.array([ 0.0 ]) )

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

python data.mean()

Part 30 --- Estimate mean and SD jointly¶

Parameterize:

text theta = [mu, log_sigma]

python def normal_nll( theta, data ): mu = theta[0]

sigma = np.exp(
    theta[1]
)

residuals = (
    data - mu
)

return (
    len(data)
    * np.log(sigma)
    + np.sum(
        residuals ** 2
    ) / (
        2 * sigma ** 2
    )
)

Reference answer¶

Estimate mu and an unconstrained parameter such as log_sigma; recover sigma = exp(log_sigma) to guarantee positivity.

Part 31 --- Optimize both parameters¶

python result = minimize( normal_nll, x0=np.array([ 0.0, 0.0 ]), args=( data, ) )

python mu_hat = result.x[0]

sigma_hat = np.exp( result.x[1] )

Part 32 --- Why log SD?¶

Explain why optimizing unrestricted sigma can be problematic.

How does:

text sigma = exp(log_sigma)

Reference answer¶

Directly optimizing sigma permits invalid nonpositive trial values. Optimizing log_sigma is unconstrained while guaranteeing a positive standard deviation.

Part 33 --- Starting-value sensitivity¶

Repeat the two-parameter optimization from:

text [0, 0] [10, 0] [-20, 2] [100, -2]

Part 34 --- Nonconvex objective¶

Define:

python def nonconvex(x): return ( np.sin( 3 * x ) + 0.1 * x ** 2 )

Part 35 --- Root finding or optimization?¶

Choose the best framing.

A¶

text x³ - 5x + 1 = 0

Part 36 --- Mini optimization study¶

Simulate:

python n = 300

x1 = rng.normal( size=n )

x2 = rng.normal( size=n )

X = np.column_stack([ np.ones(n), x1, x2 ])

beta_true = np.array([ 3.0, 1.5, -2.0 ])

y = ( X @ beta_true + rng.normal( 0, 2, size=n ) )

Reference answer¶

Compare optimizers by convergence, final objective, parameter estimates, sensitivity to starts/scaling, and agreement with a trusted solution—not speed alone.

Part 37 --- Optimization audit¶

Before trusting an optimization result, answer:

text What objective is being optimized? Are we minimizing or maximizing? Are derivatives correct? Are parameters constrained? What are the starting values? What is the stopping rule? Did the algorithm converge? Is the result sensitive to initialization? Can it be checked against an analytical/library result?

Reference answer¶

An optimization audit should record the objective, parameterization, starting values, scaling, derivative method, algorithm, stopping rule, convergence diagnostics, and trusted comparison.

Part 38 --- Git checkpoint¶

bash git status git add . git commit -m “Complete Week 12 numerical optimization exercises” git push

Part 39 --- Final reflection¶

1. Root finding¶

What does it mean to find a root?

2. Bisection¶

Why is it robust?

3. Newton--Raphson¶

Where does the update come from?

4. Starting values¶

Why do they matter?

5. Gradient¶

What information does it contain?

6. Gradient descent¶

What role does the learning rate play?

7. OLS¶

How is regression an optimization problem?

8. MLE¶

Why minimize the negative log-likelihood?

9. Constraints¶

Why optimize log(sigma) instead of sigma?

10. Validation¶

What should be checked before trusting optimizer output?

Completion checklist¶

• Created Week 12 notebook

• Defined and visualized a root problem

• Implemented bisection

• Validated bisection with SciPy

• Implemented Newton--Raphson

• Visualized Newton convergence

• Validated Newton with SciPy

• Explored poor starting values

• Implemented numerical differentiation

• Investigated finite-difference step size

• Defined an objective

• Implemented gradient descent

• Compared learning rates

• Implemented multivariable gradient descent

• Simulated regression data

• Defined MSE loss

• Derived the MSE gradient

• Numerically checked the gradient

• Fit OLS with gradient descent

• Compared with np.linalg.lstsq

• Investigated feature scaling

• Defined a negative log-likelihood

• Estimated a mean using minimize

• Estimated mean and SD jointly

• Used a log transformation for a positive parameter

• Explored starting-value sensitivity

• Explored a nonconvex objective

• Checked SciPy convergence diagnostics

• Completed the root-vs-optimization challenge

• Completed the mini optimization study

• Audited an optimization workflow

• Committed Week 12 work to Git

• Pushed work to GitHub

What you should now understand¶

text Statistical problem ↓ Equation / loss / likelihood ↓ Root or optimum ↓ Derivative information ↓ Numerical iteration ↓ Convergence diagnostics ↓ Parameter estimate ↓ Validation

Next week we will focus on **debugging, testing, defensive programming,
and exception handling**.

Reference answer¶

Root finding solves f(x)=0; optimization minimizes/maximizes an objective. Bisection is robust because it maintains a bracket, Newton can be faster but is more assumption-sensitive, and learning rate controls gradient-descent step size.