Exercises: Gradients, Hessians, and Convexity

Chapters 01 & 02: Mathematical Foundations & Unconstrained Optimization

Self-assessment guide: Work each problem on paper before expanding the answer. If you can solve 80% without peeking, move on to 02-least-squares-ceres.md.

Section A — Conceptual Questions

A1. Given $f(x, y) = x^2 y + \sin(xy)$, compute the gradient $\nabla f$ and the Hessian $H_f$. At what point(s) is the gradient zero?

A2. Prove that $f(x) = \|Ax - b\|^2$ is convex for any matrix $A$ and vector $b$. Under what condition is it strictly convex?

A3. A colleague says: "Gradient descent always converges to the global minimum if you set the learning rate small enough." Is this correct? Provide a counterexample.

A4. Why does Newton's method converge in exactly 1 step for a quadratic function? What does this tell you about its behavior near any minimum?

Section B — Computation Problems

B1. Implement gradient descent and Newton's method for the Rosenbrock function:

$$f(x, y) = (1 - x)^2 + 100(y - x^2)^2$$

Compare: - Number of iterations to converge ($\|\nabla f\| < 10^{-6}$) - Path taken (plot the trajectory) - Effect of step size on GD convergence

import numpy as np

def rosenbrock(x):
    """f(x,y) = (1-x)^2 + 100*(y-x^2)^2"""
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

def rosenbrock_grad(x):
    """Gradient of Rosenbrock function"""
    # TODO: compute and return [df/dx, df/dy]
    pass

def rosenbrock_hess(x):
    """Hessian of Rosenbrock function"""
    # TODO: compute and return 2x2 Hessian matrix
    pass

def gradient_descent(f, grad, x0, lr=0.001, tol=1e-6, max_iter=100000):
    """
    Returns: (x_final, trajectory, num_iterations)
    """
    # TODO: implement
    pass

def newton_method(f, grad, hess, x0, tol=1e-6, max_iter=100):
    """
    Returns: (x_final, trajectory, num_iterations)
    """
    # TODO: implement with line search for safety
    pass

# Start from (-1, 1)
x0 = np.array([-1.0, 1.0])
# Compare both methods. Expected: GD takes ~10,000+ iterations, Newton ~20.

import numpy as np

def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

def rosenbrock_grad(x):
    dx = -2*(1 - x[0]) + 200*(x[1] - x[0]**2)*(-2*x[0])
    dy = 200*(x[1] - x[0]**2)
    return np.array([dx, dy])

def rosenbrock_hess(x):
    dxx = 2 - 400*x[1] + 1200*x[0]**2
    dxy = -400*x[0]
    dyy = 200.0
    return np.array([[dxx, dxy], [dxy, dyy]])

def gradient_descent(f, grad, x0, lr=0.001, tol=1e-6, max_iter=100000):
    x = x0.copy()
    trajectory = [x.copy()]
    for i in range(max_iter):
        g = grad(x)
        if np.linalg.norm(g) < tol:
            return x, trajectory, i
        x = x - lr * g
        trajectory.append(x.copy())
    return x, trajectory, max_iter

def newton_method(f, grad, hess, x0, tol=1e-6, max_iter=100):
    x = x0.copy()
    trajectory = [x.copy()]
    for i in range(max_iter):
        g = grad(x)
        if np.linalg.norm(g) < tol:
            return x, trajectory, i
        H = hess(x)
        # Add damping if H is not positive definite
        try:
            dx = np.linalg.solve(H, -g)
        except np.linalg.LinAlgError:
            dx = -g  # fallback to gradient step

        # Backtracking line search
        alpha = 1.0
        while f(x + alpha * dx) > f(x) + 1e-4 * alpha * g.dot(dx):
            alpha *= 0.5
            if alpha < 1e-16:
                break

        x = x + alpha * dx
        trajectory.append(x.copy())
    return x, trajectory, max_iter

x0 = np.array([-1.0, 1.0])

x_gd, traj_gd, n_gd = gradient_descent(rosenbrock, rosenbrock_grad, x0, lr=0.001)
print(f"GD: {n_gd} iterations, x = {x_gd}")

x_nt, traj_nt, n_nt = newton_method(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0)
print(f"Newton: {n_nt} iterations, x = {x_nt}")

B2. For the function $f(x) = x_1^2 + 10 x_2^2$ (an elongated bowl):

1. What is the condition number of the Hessian?
2. How many iterations does gradient descent take vs. Newton?
3. What happens to GD if you change it to $f(x) = x_1^2 + 1000 x_2^2$?

Section C — OKS Connection

C1. The OKS navigation estimator's EKF state includes $(x, y, \theta, v, \omega)$. The position variables are in metres ($\sim 0$ to $100$), while the heading $\theta$ is in radians ($\sim 0$ to $6.28$).

1. Why might this cause numerical issues if you were solving an optimization problem with these variables?
2. How would you rescale? What scaling factors would you choose?
3. Why does the EKF avoid this issue (mostly)?