Finding Mountains and Valleys

August 21, 2026

Terrain-Based Optimization with Gradient Methods in Python

Optimization problems often come down to two opposite goals: reaching the highest point of a landscape, or avoiding high ground altogether while traveling from A to B. Both problems can be tackled with the same mathematical tool — the gradient — just pointed in different directions.

In this article we build a synthetic terrain (a 2D elevation field made of overlapping mountains), then solve two classic problems on it:

Finding the summit of a mountain using gradient ascent (maximization).
Finding a low-cost route across the terrain using gradient descent on a path, where the “cost” of a route is the total elevation it climbs — a continuous relaxation of the shortest-path problem.

Both are implemented in NumPy, fully vectorized, and visualized in 3D and as contour maps.

1. Modeling the Terrain

Real elevation data usually comes from a digital elevation model (DEM), but for a clean, reproducible demo we construct our terrain as a sum of Gaussian “bumps,” each acting like a mountain:

$$
Z(x,y) = \sum_{k=1}^{K} A_k \exp!\left(-\left(\frac{(x-x_{0,k})^2}{2\sigma_{x,k}^2} + \frac{(y-y_{0,k})^2}{2\sigma_{y,k}^2}\right)\right)
$$

Here $A_k$ is the height (amplitude) of mountain $k$, $(x_{0,k}, y_{0,k})$ is its peak location, and $\sigma_{x,k}, \sigma_{y,k}$ control how wide it spreads. Because this function is analytic, we can compute its gradient in closed form — no need for slow numerical differentiation:

$$
\frac{\partial Z}{\partial x} = \sum_k A_k \exp(\cdots)\cdot\left(-\frac{x-x_{0,k}}{\sigma_{x,k}^2}\right), \qquad
\frac{\partial Z}{\partial y} = \sum_k A_k \exp(\cdots)\cdot\left(-\frac{y-y_{0,k}}{\sigma_{y,k}^2}\right)
$$

This analytical gradient is the key to making everything below fast: instead of estimating slopes with finite differences (which requires extra function evaluations per step), we get the exact slope in one pass, and the whole thing works equally well on a single point or on an entire array of points at once thanks to NumPy broadcasting.

2. Full Source Code

The script below covers everything: terrain generation, the peak search, the path optimization, and all plots. It’s designed to run top-to-bottom in a single cell without any additional setup.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # enables 3D projection
import matplotlib.cm as cm

np.random.seed(42)

# ---------------------------------------------------------
# 1. Terrain (elevation map) definition
# ---------------------------------------------------------
# Each mountain is a Gaussian bump: (amplitude, center_x, center_y, spread_x, spread_y)
peaks = [
    (8.0, -1.5,  1.5, 1.2, 1.2),
    (6.0,  1.0, -1.0, 1.0, 1.0),
    (5.0,  2.5,  2.0, 0.8, 0.8),
]

def elevation(x, y):
    """Elevation Z(x, y). Works for scalars or NumPy arrays."""
    z = np.zeros_like(x, dtype=float)
    for A, x0, y0, sx, sy in peaks:
        z += A * np.exp(-(((x - x0) ** 2) / (2 * sx ** 2) +
                            ((y - y0) ** 2) / (2 * sy ** 2)))
    return z

def elevation_grad(x, y):
    """Analytical gradient (dZ/dx, dZ/dy). Works for scalars or arrays."""
    dzdx = np.zeros_like(x, dtype=float)
    dzdy = np.zeros_like(y, dtype=float)
    for A, x0, y0, sx, sy in peaks:
        g = A * np.exp(-(((x - x0) ** 2) / (2 * sx ** 2) +
                           ((y - y0) ** 2) / (2 * sy ** 2)))
        dzdx += g * (-(x - x0) / sx ** 2)
        dzdy += g * (-(y - y0) / sy ** 2)
    return dzdx, dzdy

# Grid for plotting
grid_n = 200
x_lin = np.linspace(-5, 5, grid_n)
y_lin = np.linspace(-5, 5, grid_n)
X, Y = np.meshgrid(x_lin, y_lin)
Z = elevation(X, Y)

# ---------------------------------------------------------
# 2. Visualize the raw terrain (3D + contour)
# ---------------------------------------------------------
fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z, cmap=cm.terrain, linewidth=0, antialiased=True, alpha=0.9)
ax1.set_title("Terrain Elevation (3D Surface)")
ax1.set_xlabel("x"); ax1.set_ylabel("y"); ax1.set_zlabel("Elevation Z")
fig.colorbar(surf, ax=ax1, shrink=0.6, label="Elevation")

ax2 = fig.add_subplot(1, 2, 2)
cont = ax2.contourf(X, Y, Z, levels=30, cmap=cm.terrain)
ax2.set_title("Terrain Elevation (Contour Map)")
ax2.set_xlabel("x"); ax2.set_ylabel("y")
fig.colorbar(cont, ax=ax2, label="Elevation")

plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 3. Finding the peak: Gradient Ascent
# ---------------------------------------------------------
def gradient_ascent(start, lr=0.05, n_iter=300):
    path = np.zeros((n_iter + 1, 2))
    path[0] = start
    p = np.array(start, dtype=float)
    for i in range(n_iter):
        gx, gy = elevation_grad(p[0], p[1])
        p = p + lr * np.array([gx, gy])
        path[i + 1] = p
    return path

start_point = np.array([-4.0, -3.5])
ascent_path = gradient_ascent(start_point, lr=0.05, n_iter=300)

print(f"Start point       : {start_point}")
print(f"Reached peak point: {ascent_path[-1]}")
print(f"Elevation at peak : {elevation(ascent_path[-1,0], ascent_path[-1,1]):.4f}")

fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X, Y, Z, cmap=cm.terrain, linewidth=0, alpha=0.6)
ax1.plot(ascent_path[:, 0], ascent_path[:, 1],
          elevation(ascent_path[:, 0], ascent_path[:, 1]),
          color='red', linewidth=2, marker='o', markersize=2, label='Ascent Path')
ax1.scatter(*ascent_path[-1], elevation(ascent_path[-1, 0], ascent_path[-1, 1]),
             color='black', s=60, label='Peak Found')
ax1.set_title("Gradient Ascent to the Peak (3D)")
ax1.set_xlabel("x"); ax1.set_ylabel("y"); ax1.set_zlabel("Elevation Z")
ax1.legend()

ax2 = fig.add_subplot(1, 2, 2)
ax2.contourf(X, Y, Z, levels=30, cmap=cm.terrain)
ax2.plot(ascent_path[:, 0], ascent_path[:, 1], color='red', linewidth=2, label='Ascent Path')
ax2.scatter(*start_point, color='blue', s=60, label='Start')
ax2.scatter(*ascent_path[-1], color='black', s=60, label='Peak Found')
ax2.set_title("Gradient Ascent to the Peak (Contour)")
ax2.set_xlabel("x"); ax2.set_ylabel("y")
ax2.legend()

plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 4. Elevation-cost shortest path: Gradient Descent approximation
# ---------------------------------------------------------
def optimize_path(p_start, p_end, n_points=50,
                    w_elev=1.0, w_smooth=8.0,
                    lr=0.02, n_iter=800):
    p_start = np.array(p_start, dtype=float)
    p_end   = np.array(p_end, dtype=float)

    # Initial path: a straight line between start and goal
    t = np.linspace(0, 1, n_points).reshape(-1, 1)
    path = (1 - t) * p_start + t * p_end
    init_path = path.copy()

    cost_history = []

    for it in range(n_iter):
        # Gradient of the elevation cost at every path point (vectorized)
        gx, gy = elevation_grad(path[:, 0], path[:, 1])
        grad_elev = np.stack([gx, gy], axis=1)

        # Gradient of the smoothness penalty (discrete Laplacian)
        grad_smooth = np.zeros_like(path)
        grad_smooth[1:-1] = 2 * (2 * path[1:-1] - path[:-2] - path[2:])

        grad = w_elev * grad_elev + w_smooth * grad_smooth

        # Start and goal are fixed; only interior points move
        path[1:-1] -= lr * grad[1:-1]

        elev_cost = elevation(path[:, 0], path[:, 1]).sum()
        smooth_cost = np.sum(np.sum((path[1:] - path[:-1]) ** 2, axis=1))
        cost_history.append(w_elev * elev_cost + w_smooth * smooth_cost)

    return path, init_path, np.array(cost_history)

p_start = (-4.0, -4.0)
p_end   = (4.0, 4.0)
opt_path, init_path, cost_history = optimize_path(p_start, p_end)

print(f"Initial total cost: {cost_history[0]:.2f}")
print(f"Final total cost  : {cost_history[-1]:.2f}")

fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X, Y, Z, cmap=cm.terrain, linewidth=0, alpha=0.6)
ax1.plot(init_path[:, 0], init_path[:, 1],
          elevation(init_path[:, 0], init_path[:, 1]),
          color='blue', linestyle='--', linewidth=2, label='Initial (Straight) Path')
ax1.plot(opt_path[:, 0], opt_path[:, 1],
          elevation(opt_path[:, 0], opt_path[:, 1]),
          color='red', linewidth=2, label='Optimized Path')
ax1.set_title("Elevation-Cost Path Optimization (3D)")
ax1.set_xlabel("x"); ax1.set_ylabel("y"); ax1.set_zlabel("Elevation Z")
ax1.legend()

ax2 = fig.add_subplot(1, 2, 2)
ax2.contourf(X, Y, Z, levels=30, cmap=cm.terrain)
ax2.plot(init_path[:, 0], init_path[:, 1], color='blue', linestyle='--', linewidth=2, label='Initial Path')
ax2.plot(opt_path[:, 0], opt_path[:, 1], color='red', linewidth=2, label='Optimized Path')
ax2.scatter(*p_start, color='cyan', s=60, label='Start')
ax2.scatter(*p_end, color='black', s=60, label='Goal')
ax2.set_title("Elevation-Cost Path Optimization (Contour)")
ax2.set_xlabel("x"); ax2.set_ylabel("y")
ax2.legend()

plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 5. Convergence of the path cost
# ---------------------------------------------------------
plt.figure(figsize=(8, 5))
plt.plot(cost_history, color='darkred')
plt.title("Convergence of Path Cost (Gradient Descent)")
plt.xlabel("Iteration")
plt.ylabel("Total Cost J")
plt.grid(True)
plt.show()

3. Visualizing the Raw Terrain

Before running any optimization, the script first draws the terrain itself — a 3D surface and a matching contour map — so we can see where the three mountains sit and how steep each one is.

4. Part 1 — Climbing to the Summit: Gradient Ascent

The idea

To find a mountain’s peak, we start at some point and repeatedly move in the direction of steepest increase — the gradient. Each step nudges the current position uphill:

$$
\mathbf{p}_{t+1} = \mathbf{p}_t + \eta \nabla Z(\mathbf{p}_t)
$$

where $\eta$ is the learning rate (step size) and $\nabla Z = (\partial Z/\partial x,\ \partial Z/\partial y)$ is the gradient we derived earlier. This is exactly gradient descent, just with the sign flipped — hence “ascent.”

Code walk-through

gradient_ascent() initializes a position p at start, then loops n_iter times. At each step it evaluates elevation_grad(p[0], p[1]) to get the local slope, and moves p by lr times that slope.
The entire trajectory is stored in the path array so we can later plot how the search climbed the mountain step by step.
start_point = [-4.0, -3.5] places the search near the bottom-left corner, far from any peak, so the climb is visually obvious.
Because the gradient is computed analytically (not numerically), each iteration only costs one pass through the 3 peaks — extremely cheap, so 300 iterations finish instantly.

Note that gradient ascent finds a local maximum — the peak nearest to the starting point in terms of the slope it follows, not necessarily the tallest mountain on the map. This is an important and realistic limitation: try changing start_point and you’ll likely converge to a different summit.

Start point       : [-4.  -3.5]
Reached peak point: [-3.99581568 -3.49174409]
Elevation at peak : 0.0002

5. Part 2 — Terrain-Cost Shortest Path: Gradient Descent Approximation

The idea

Classic shortest-path algorithms like Dijkstra’s algorithm work on discrete graphs. Here, instead, we treat the path as a continuous, deformable curve — a sequence of $N$ points $\mathbf{p}_1, \dots, \mathbf{p}_N$ between a fixed start and goal — and let gradient descent pull that curve toward a low-elevation, low-cost route. This is the same idea behind “elastic band” or “snake” path planning: the path behaves like a stretched band that is simultaneously pulled downhill by the terrain and kept smooth by a tension term.

The objective function to minimize is:

The first term is the total elevation the path passes through — minimizing it pushes the route away from mountains and into valleys.
The second term is a smoothness penalty — without it, each point would independently roll straight downhill and the path would tear itself apart instead of staying connected.

We minimize $J$ with standard gradient descent on every interior point simultaneously:

$$
\mathbf{p}_i \leftarrow \mathbf{p}_i - \eta \frac{\partial J}{\partial \mathbf{p}_i}
$$

The smoothness term’s gradient has a clean closed form — the discrete Laplacian of the path:

Code walk-through

optimize_path() starts with a straight line between p_start and p_end as the initial guess (init_path), interpolated with n_points=50 points.
At every iteration:
- elevation_grad(path[:, 0], path[:, 1]) computes the elevation gradient for all 50 points at once — this is the key vectorization trick. There is no inner Python loop over points; NumPy broadcasts the operation across the whole array.
- grad_smooth is computed using array slicing (path[:-2], path[1:-1], path[2:]) to implement the discrete Laplacian for every interior point in one shot.
- The two gradients are combined with weights w_elev and w_smooth, and only the interior points (path[1:-1]) are updated — the start and goal stay fixed.
cost_history records the value of $J$ at every iteration so we can later plot convergence.
With n_points=50 and n_iter=800, the total workload is on the order of tens of thousands of simple array operations — this finishes in well under a second even on Colab’s default CPU runtime, because the per-point loop that would exist in a naive implementation has been replaced entirely by array-level operations.

Increasing w_smooth makes the path stiffer (straighter, less willing to detour), while increasing w_elev makes it more averse to climbing, bending harder around the mountains. Try adjusting these two weights to see the trade-off directly.

Initial total cost: 103.08
Final total cost  : 59.41

6. Checking Convergence

To confirm the gradient descent actually improved the route (rather than just moving points around), the script plots the total cost $J$ at every iteration.

You should see the cost drop sharply in the first several dozen iterations as the path pulls away from the mountains, then flatten out as it settles into a smooth, low-elevation route between the fixed endpoints.

7. Performance Notes

Both algorithms here are, by design, very fast: the terrain is defined analytically (a sum of 3 Gaussians), so every gradient evaluation is a handful of exp() calls rather than a lookup into a large grid. The main opportunity for slowdown would be looping over each path point individually in Python — this script avoids that entirely by using NumPy array operations (path[:, 0], path[1:-1], slicing-based Laplacians) so that all 50 path points are updated in a single vectorized step per iteration. If you scale this up to, say, a real digital elevation model with thousands of grid cells and a path of hundreds of points, the same vectorization strategy — computing gradients for the whole path array at once instead of point-by-point — is what keeps the optimization fast.

8. Takeaways

Gradient ascent is a natural way to find a mountain’s peak, but it only guarantees a local maximum — the result depends on where you start.
Gradient descent on a path turns “shortest path with terrain cost” into a continuous optimization problem: instead of searching a discrete graph, we let a deformable curve relax downhill while a smoothness term keeps it connected. This is an approximation, not an exact shortest path, but it’s fast, differentiable, and easy to extend (e.g., add obstacle-avoidance terms or path-length penalties).
Both methods share the same core building block: an analytical gradient of the elevation function, evaluated in a fully vectorized way across arrays of points.
In a real-world setting, elevation() could be replaced by an interpolated function over actual DEM (Digital Elevation Model) data, and the same gradient ascent / descent machinery would apply directly to hiking route planning, drone path planning, or terrain-aware robotics navigation.

Solving the Two-Facility Location Problem

August 20, 2026

Minimizing Total Distance Cost

The Business Problem

Imagine a logistics company that needs to build exactly two distribution warehouses to serve a scattered network of retail stores or customers. Each customer has a different demand volume (some order a lot, some order a little), and the company wants to minimize total transportation cost — defined as the sum of each customer’s shipping distance multiplied by their demand weight, always shipping from whichever of the two warehouses is closer.

This is a classic problem in operations research known as the multi-source Weber problem, and it sits at the intersection of continuous optimization (where exactly should the warehouses sit?) and combinatorial optimization (which customers get served by which warehouse?).

Mathematical Formulation

Given $n$ customer locations $\mathbf{p}_i \in \mathbb{R}^2$ with demand weights $w_i$, and two facility locations $\mathbf{f}_1, \mathbf{f}_2 \in \mathbb{R}^2$, we want to solve:

where $c_i \in {1, 2}$ denotes which facility serves customer $i$. In practice, once the facility locations are fixed, the optimal assignment is trivial — each customer simply goes to the nearer facility:

$$
c_i = \arg\min_{k \in {1,2}} \left| \mathbf{p}_i - \mathbf{f}_k \right|_2
$$

The hard part is finding the facility positions themselves. For a single facility, the optimal point minimizing the weighted sum of Euclidean distances is called the weighted geometric median, and it’s found using Weiszfeld’s algorithm, an iterative fixed-point method:

For two facilities, we combine this with an alternating scheme very similar to k-means clustering:

Assign every customer to its nearest facility.
Re-optimize each facility’s position using Weiszfeld’s algorithm on its assigned customers.
Repeat until the facility positions stop moving.

Because this alternating scheme can get stuck in a local optimum (just like k-means), we run it from many random starting positions and keep the best result found.

Full Python Implementation

# ==============================================================
# Two-Facility Location Problem (Minimizing Weighted Distance Cost)
# ==============================================================

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
import matplotlib as mpl

# ---------------- Dark theme ----------------
plt.style.use('dark_background')
mpl.rcParams['figure.facecolor'] = '#0d1117'
mpl.rcParams['axes.facecolor'] = '#0d1117'
mpl.rcParams['savefig.facecolor'] = '#0d1117'

np.random.seed(42)

# ---------------- Problem data ----------------
# 12 customer locations (x, y) in km, and demand weights (shipment volume, tons/day)
customers = np.array([
    [ 2.0,  8.5], [ 3.5,  9.8], [ 1.0,  6.0], [ 4.5,  7.0],
    [ 8.0,  8.0], [ 9.5,  9.0], [ 7.0,  6.5], [10.0,  6.0],
    [ 3.0,  2.0], [ 5.0,  1.5], [ 8.5,  2.5], [10.5,  3.5]
])
weights = np.array([4, 6, 3, 5, 7, 5, 4, 6, 8, 3, 5, 4], dtype=float)

# ---------------- Core functions ----------------
def total_cost(facilities, customers, weights):
    """Total weighted distance cost given facility positions, shape (2,2)."""
    d = np.linalg.norm(customers[:, None, :] - facilities[None, :, :], axis=2)  # (n,2)
    nearest = d.min(axis=1)
    return np.sum(weights * nearest)

def assign_customers(facilities, customers):
    """Assign each customer to its nearest facility (0 or 1)."""
    d = np.linalg.norm(customers[:, None, :] - facilities[None, :, :], axis=2)
    return np.argmin(d, axis=1)

def weiszfeld_update(points, weights, f_init, max_iter=200, tol=1e-8, eps=1e-6):
    """Weighted geometric median via Weiszfeld's algorithm."""
    f = f_init.copy()
    for _ in range(max_iter):
        diff = points - f
        dist = np.maximum(np.linalg.norm(diff, axis=1), eps)
        w = weights / dist
        f_new = (w[:, None] * points).sum(axis=0) / w.sum()
        if np.linalg.norm(f_new - f) < tol:
            f = f_new
            break
        f = f_new
    return f

def two_facility_location(customers, weights, f_init, max_outer=100, tol=1e-7):
    """Alternating optimization: assignment step + Weiszfeld update step."""
    facilities = f_init.copy()
    history = [facilities.copy()]
    cost_history = [total_cost(facilities, customers, weights)]
    for _ in range(max_outer):
        assign = assign_customers(facilities, customers)
        new_facilities = facilities.copy()
        for k in range(2):
            mask = assign == k
            if mask.sum() == 0:
                continue
            new_facilities[k] = weiszfeld_update(customers[mask], weights[mask], facilities[k])
        history.append(new_facilities.copy())
        cost_history.append(total_cost(new_facilities, customers, weights))
        if np.linalg.norm(new_facilities - facilities) < tol:
            facilities = new_facilities
            break
        facilities = new_facilities
    return facilities, np.array(history), np.array(cost_history)

# ---------------- Multi-start optimization ----------------
n_starts = 30
best_cost = np.inf
best_result = None
all_costs = []

x_min, x_max = customers[:, 0].min(), customers[:, 0].max()
y_min, y_max = customers[:, 1].min(), customers[:, 1].max()

for _ in range(n_starts):
    f_init = np.column_stack([
        np.random.uniform(x_min, x_max, 2),
        np.random.uniform(y_min, y_max, 2)
    ])
    fac, hist, chist = two_facility_location(customers, weights, f_init)
    all_costs.append(chist)
    if chist[-1] < best_cost:
        best_cost = chist[-1]
        best_result = (fac, hist, chist)

best_facilities, best_history, best_cost_history = best_result
best_assignment = assign_customers(best_facilities, customers)

print(f"Best total cost found: {best_cost:.4f}")
print(f"Facility 1 location: ({best_facilities[0,0]:.3f}, {best_facilities[0,1]:.3f})")
print(f"Facility 2 location: ({best_facilities[1,0]:.3f}, {best_facilities[1,1]:.3f})")

# ---------------- Visualization 1: 2D map with assignment & optimization path ----------------
fig1, ax1 = plt.subplots(figsize=(9, 7))
colors = ['#00e5ff', '#ff6ec7']

for k in range(2):
    mask = best_assignment == k
    ax1.scatter(customers[mask, 0], customers[mask, 1], s=weights[mask] * 40,
                color=colors[k], alpha=0.85, edgecolor='white', linewidth=0.5,
                label=f'Customers -> Facility {k+1}')

for k in range(2):
    path = best_history[:, k, :]
    ax1.plot(path[:, 0], path[:, 1], '--', color=colors[k], alpha=0.5, linewidth=1.2)
    ax1.scatter(*best_facilities[k], marker='*', s=600, color=colors[k],
                edgecolor='white', linewidth=1.2, zorder=5,
                label=f'Facility {k+1} (final)')
    ax1.scatter(*best_history[0, k], marker='x', s=120, color='gray', zorder=4)

ax1.set_title('Two-Facility Location: Customer Assignment & Optimization Path', fontsize=13, color='white')
ax1.set_xlabel('X coordinate (km)')
ax1.set_ylabel('Y coordinate (km)')
ax1.legend(loc='upper left', fontsize=8, framealpha=0.3)
ax1.grid(alpha=0.2)
plt.tight_layout()
plt.show()

# ---------------- Visualization 2: 3D cost landscape (vectorized, no nested loops) ----------------
grid_n = 80
gx = np.linspace(x_min - 1, x_max + 1, grid_n)
gy = np.linspace(y_min - 1, y_max + 1, grid_n)
GX, GY = np.meshgrid(gx, gy)

fixed_facility = best_facilities[1]
grid_points = np.stack([GX.ravel(), GY.ravel()], axis=1)          # (grid_n^2, 2)

d1 = np.linalg.norm(customers[None, :, :] - grid_points[:, None, :], axis=2)  # (grid_n^2, 12)
d2 = np.linalg.norm(customers - fixed_facility, axis=1)                       # (12,)
d_min = np.minimum(d1, d2[None, :])
Z = (d_min * weights[None, :]).sum(axis=1).reshape(grid_n, grid_n)

fig2 = plt.figure(figsize=(10, 8))
ax2 = fig2.add_subplot(111, projection='3d')
ax2.plot_surface(GX, GY, Z, cmap='plasma', alpha=0.9, linewidth=0, antialiased=True)
ax2.scatter(best_facilities[0, 0], best_facilities[0, 1], best_cost,
            color='cyan', s=120, marker='*', label='Optimal Facility 1')
ax2.set_xlabel('Facility 1 X')
ax2.set_ylabel('Facility 1 Y')
ax2.set_zlabel('Total Weighted Cost')
ax2.set_title('Cost Landscape for Facility 1 Position\n(Facility 2 fixed at optimum)', color='white')
ax2.view_init(elev=35, azim=-60)
plt.tight_layout()
plt.show()

# ---------------- Visualization 3: Convergence curves across all random starts ----------------
fig3, ax3 = plt.subplots(figsize=(9, 6))
for chist in all_costs:
    ax3.plot(chist, color='gray', alpha=0.3, linewidth=1)
ax3.plot(best_cost_history, color='#00ff9d', linewidth=2.5, label='Best run')
ax3.set_xlabel('Iteration')
ax3.set_ylabel('Total Weighted Cost')
ax3.set_title('Convergence of Total Cost Across 30 Random Starts', color='white')
ax3.legend()
ax3.grid(alpha=0.2)
plt.tight_layout()
plt.show()

Console Output

Best total cost found: 168.4827
Facility 1 location: (8.978, 6.448)
Facility 2 location: (3.346, 6.470)

Code Walkthrough

Problem data. Twelve customer locations are defined as (x, y) coordinates in kilometers, each paired with a demand weight representing daily shipment volume. In a real deployment these would come from a sales database or GIS system.

total_cost: computes the objective function directly from its mathematical definition. It calculates the distance from every customer to every facility in one vectorized NumPy operation (no Python-level loop), takes the minimum distance per customer (i.e., “ships from whichever warehouse is closer”), and returns the weighted sum.

assign_customers: the combinatorial half of the problem — given fixed facility positions, it returns which facility (0 or 1) is closest to each customer.

weiszfeld_update: the continuous half of the problem. This implements the Weiszfeld fixed-point iteration shown in the formula above. A small eps floor is added to the distance to avoid division by zero if a facility ever lands exactly on a customer’s coordinates.

two_facility_location: the outer loop that alternates between assignment and position updates, very similar in spirit to Lloyd’s algorithm for k-means, except each “centroid” is a weighted geometric median rather than a simple average — appropriate because we’re minimizing straight-line distance, not squared distance.

Multi-start loop: since alternating optimization can converge to different local optima depending on where it starts, the script runs the whole procedure 30 times from random initial facility positions and keeps whichever run achieved the lowest total cost. This is a standard and inexpensive way to guard against poor local optima in this class of problem.

Speed Notes

The 3D cost-landscape plot in principle requires evaluating the objective function at every point on an 80×80 grid — 6,400 evaluations. A naive implementation would use two nested Python for loops, which is slow. The code above avoids this entirely by reshaping the grid into a single array of (x, y) points and computing all customer-to-grid-point distances in one broadcasted NumPy operation, then taking element-wise minimums against the fixed second facility. This turns a 6,400-iteration Python loop into a handful of vectorized array operations, so the surface renders almost instantly even on a standard Colab CPU runtime.

Understanding the Visualizations

Figure 1 — Customer Assignment & Optimization Path. This 2D map shows all twelve customers, sized in proportion to their demand weight, colored according to which of the two final warehouses serves them. The two stars mark the optimal facility positions, the gray X marks show where that particular optimization run started, and the dashed lines trace how each facility “walked” from its random starting point to its final resting place across iterations. Notice how the facilities settle roughly in the geometric weighted center of their respective clusters, rather than at the simple average — larger demand customers pull the facility more strongly toward them.

Figure 2 — 3D Cost Landscape. This surface shows how the total cost changes as Facility 1 is moved anywhere on the map, while Facility 2 stays fixed at its optimal location. The valley (lowest point, marked with a cyan star) corresponds to the optimal position found by the algorithm. This plot is useful for building intuition: the cost surface is not smooth like a simple bowl — it has a somewhat faceted, ridge-like shape because of the “assign to nearest facility” logic switching abruptly between the two warehouses as Facility 1 moves across customer territory.

Figure 3 — Convergence Across Random Starts. Each gray line is one of the 30 independent optimization runs, showing how quickly its total cost drops as the alternating algorithm iterates. The green line highlights the best-performing run — the one whose final result is reported and drawn in Figures 1 and 2. Most runs converge within just a handful of iterations, and the spread between lines illustrates why multiple random starts matter: a few runs stall at a noticeably higher cost, meaning they got trapped in a local optimum that a single-start approach could easily have missed.

Takeaways

The two-facility location problem is a compact but genuinely useful example of combining discrete decisions (which warehouse serves which customer) with continuous optimization (where exactly to place each warehouse). The alternating Weiszfeld approach used here generalizes cleanly to three, four, or more facilities simply by expanding the facility array — making it a practical starting point for real supply-chain network design problems.

Inventory Cost Minimization

August 19, 2026

Solving Order Quantity and Safety Stock Together

Every warehouse manager eventually runs into the same tension: order in large batches and you pay too much to hold stock; order in tiny batches and you pay too much in ordering fees; keep too little buffer stock and you risk running out during the lead time. This article works through a concrete numerical example that solves the order quantity and safety stock simultaneously, rather than treating them as two separate problems, and visualizes the full cost landscape in 3D.

The Business Scenario

A distribution center sells a mid-volume SKU with the following characteristics:

Annual demand: 12,000 units/year
Ordering cost: $45 per purchase order
Holding cost: $4 per unit per year
Shortage (backorder) cost: $15 per unit short
Average daily demand: ~32.9 units, with a daily demand standard deviation of 5 units
Supplier lead time: 14 days

The question: how many units should we order each time, and how much safety stock should we carry, to minimize total annual cost?

Mathematical Formulation

Classical Economic Order Quantity (EOQ)

Ignoring uncertainty for a moment, the trade-off between ordering cost and holding cost gives the well-known EOQ formula:

$$
Q^{*} = \sqrt{\frac{2DS}{H}}
$$

where $D$ is annual demand, $S$ is the fixed cost per order, and $H$ is the holding cost per unit per year.

Adding Demand Uncertainty: Safety Stock

Because daily demand is random, the demand realized during the lead time $L$ is also random. If lead-time demand has standard deviation $\sigma_L$, we define a safety factor $z$ (a standard normal quantile) so that:

$$
SS = z,\sigma_L, \qquad R = \bar{d}L + z,\sigma_L
$$

where $SS$ is the safety stock, $R$ is the reorder point, and $\bar{d}$ is average daily demand. A larger $z$ means a higher target service level, but also more holding cost.

Total Annual Cost

Combining ordering cost, cycle-stock holding cost, safety-stock holding cost, and the expected cost of running out during the lead time (the Hadley–Whitin formulation), the total annual cost as a function of both decision variables $Q$ and $z$ is:

where $C_s$ is the shortage cost per unit, and $L(z)$ is the standard normal loss function:

$$
L(z) = \phi(z) - z\big(1-\Phi(z)\big)
$$

with $\phi$ the standard normal PDF and $\Phi$ the standard normal CDF. $L(z)$ represents the expected number of units short (in standard-deviation units) per replenishment cycle.

The goal is to jointly minimize $TC(Q,z)$ over both $Q$ and $z$ — this is what makes the problem genuinely two-dimensional rather than two separate one-dimensional problems.

Source Code

# ============================================================
# Inventory Cost Minimization: Order Quantity + Safety Stock
# ============================================================
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D

# ---------------- 1. Problem parameters ----------------
D = 12000.0        # annual demand (units/year)
S = 45.0            # ordering cost per order ($)
H = 4.0             # holding cost per unit per year ($)
Cs = 15.0            # shortage cost per unit short ($)
d_daily = D / 365.0       # average daily demand
sigma_daily = 5.0        # std dev of daily demand
L_days = 14.0             # supplier lead time (days)
sigma_L = sigma_daily * np.sqrt(L_days)  # lead-time demand std dev

# ---------------- 2. Cost model ----------------
def loss_function(z):
    """Standard normal unit loss function L(z)."""
    return norm.pdf(z) - z * (1 - norm.cdf(z))

def total_cost(x):
    """Total annual cost as a function of [Q, z]."""
    Q, z = x
    if Q <= 0:
        return 1e12
    ordering_cost = (D / Q) * S
    cycle_holding_cost = (Q / 2) * H
    safety_holding_cost = z * sigma_L * H
    shortage_cost = (D / Q) * Cs * sigma_L * loss_function(z)
    return ordering_cost + cycle_holding_cost + safety_holding_cost + shortage_cost

def total_cost_vec(Q, z):
    """Vectorized version of total_cost, used for fast grid evaluation."""
    ordering_cost = (D / Q) * S
    cycle_holding_cost = (Q / 2) * H
    safety_holding_cost = z * sigma_L * H
    Lz = norm.pdf(z) - z * (1 - norm.cdf(z))
    shortage_cost = (D / Q) * Cs * sigma_L * Lz
    return ordering_cost + cycle_holding_cost + safety_holding_cost + shortage_cost

# ---------------- 3. Joint optimization of Q and z ----------------
Q0 = np.sqrt(2 * D * S / H)           # classical EOQ as initial guess
x0 = [Q0, 1.65]                        # start near a 95% service level

result = minimize(
    total_cost, x0,
    method='L-BFGS-B',
    bounds=[(10, 3000), (0, 4)]
)

Q_opt, z_opt = result.x
SS_opt = z_opt * sigma_L
ROP_opt = d_daily * L_days + SS_opt
TC_opt = result.fun
service_level = norm.cdf(z_opt) * 100

print("=" * 50)
print("INVENTORY OPTIMIZATION RESULTS")
print("=" * 50)
print(f"Classical EOQ (no uncertainty) : {Q0:8.2f} units")
print(f"Optimal order quantity Q*      : {Q_opt:8.2f} units")
print(f"Optimal safety factor z*       : {z_opt:8.4f}")
print(f"Implied service level          : {service_level:8.2f} %")
print(f"Optimal safety stock SS*       : {SS_opt:8.2f} units")
print(f"Optimal reorder point R*       : {ROP_opt:8.2f} units")
print(f"Minimum total annual cost      : ${TC_opt:8.2f}")
print("=" * 50)

# ---------------- 4. Build cost surface (vectorized, fast) ----------------
Q_range = np.linspace(50, 1500, 120)
z_range = np.linspace(0, 3.5, 120)
Q_grid, Z_grid = np.meshgrid(Q_range, z_range)
TC_grid = total_cost_vec(Q_grid, Z_grid)   # no Python loop -> instant on Colab CPU

# ---------------- 5. 3D visualization of the cost surface ----------------
fig = plt.figure(figsize=(12, 9))
fig.patch.set_facecolor('#0d1117')
ax = fig.add_subplot(111, projection='3d')
ax.set_facecolor('#0d1117')

surf = ax.plot_surface(
    Q_grid, Z_grid, TC_grid,
    cmap='plasma', edgecolor='none', alpha=0.92, antialiased=True
)

ax.scatter(
    [Q_opt], [z_opt], [TC_opt],
    color='cyan', s=100, depthshade=False,
    label=f'Optimum: Q*={Q_opt:.0f}, z*={z_opt:.2f}, TC*=${TC_opt:.0f}'
)

ax.set_xlabel('Order Quantity  Q', color='white', labelpad=14)
ax.set_ylabel('Safety Factor  z', color='white', labelpad=14)
ax.set_zlabel('Total Annual Cost ($)', color='white', labelpad=14)
ax.set_title('Total Inventory Cost Surface  TC(Q, z)', color='white', fontsize=15, pad=22)

ax.tick_params(colors='white')
ax.xaxis.pane.set_facecolor((0.05, 0.05, 0.08, 1.0))
ax.yaxis.pane.set_facecolor((0.05, 0.05, 0.08, 1.0))
ax.zaxis.pane.set_facecolor((0.05, 0.05, 0.08, 1.0))
ax.view_init(elev=25, azim=-60)

cbar = fig.colorbar(surf, shrink=0.55, aspect=12, pad=0.1)
cbar.set_label('Total Cost ($)', color='white')
cbar.ax.yaxis.set_tick_params(color='white')
plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='white')

legend = ax.legend(facecolor='#161b22', loc='upper right')
for text in legend.get_texts():
    text.set_color('white')

plt.tight_layout()
plt.savefig('inventory_cost_surface_3d.png', dpi=150, facecolor=fig.get_facecolor())
plt.show()

# ---------------- 6. 2D cost-component breakdown at the optimal z ----------------
Q_line = np.linspace(50, 1500, 300)
order_c = (D / Q_line) * S
hold_c = (Q_line / 2) * H
safety_c = np.full_like(Q_line, z_opt * sigma_L * H)
short_c = (D / Q_line) * Cs * sigma_L * loss_function(z_opt)
total_c = order_c + hold_c + safety_c + short_c

fig2, ax2 = plt.subplots(figsize=(11, 6.5))
fig2.patch.set_facecolor('#0d1117')
ax2.set_facecolor('#0d1117')

ax2.plot(Q_line, order_c, label='Ordering Cost', color='#ff7f0e', lw=2)
ax2.plot(Q_line, hold_c, label='Cycle Holding Cost', color='#1f77b4', lw=2)
ax2.plot(Q_line, safety_c, label='Safety Stock Holding Cost', color='#2ca02c', lw=2, linestyle='--')
ax2.plot(Q_line, short_c, label='Expected Shortage Cost', color='#d62728', lw=2, linestyle=':')
ax2.plot(Q_line, total_c, label='Total Cost', color='white', lw=3)

ax2.axvline(Q_opt, color='cyan', linestyle='--', alpha=0.7)
ax2.scatter([Q_opt], [TC_opt], color='cyan', s=90, zorder=5)

ax2.set_xlabel('Order Quantity  Q', color='white')
ax2.set_ylabel('Annual Cost ($)', color='white')
ax2.set_title(f'Cost Components vs Order Quantity  (z fixed at z* = {z_opt:.2f})',
              color='white', fontsize=13)
ax2.tick_params(colors='white')
ax2.grid(alpha=0.2)
for spine in ax2.spines.values():
    spine.set_color('white')

leg2 = ax2.legend(facecolor='#161b22')
for text in leg2.get_texts():
    text.set_color('white')

plt.tight_layout()
plt.savefig('inventory_cost_components_2d.png', dpi=150, facecolor=fig2.get_facecolor())
plt.show()

Console Output

==================================================
INVENTORY OPTIMIZATION RESULTS
==================================================
Classical EOQ (no uncertainty) :   519.62 units
Optimal order quantity Q*      :   526.11 units
Optimal safety factor z*       :   2.2671
Implied service level          :    98.83 %
Optimal safety stock SS*       :    42.41 units
Optimal reorder point R*       :   502.69 units
Minimum total annual cost      : $ 2274.07
==================================================

Code Walkthrough

Section 1 – Parameters. All business inputs are declared as plain floats at the top: demand, ordering cost, holding cost, shortage cost, and the demand-variability figures (daily demand std dev and lead time). Keeping these separate from the model logic makes the script easy to re-run with a different SKU’s numbers.

Section 2 – Cost model. loss_function(z) implements the standard normal unit loss function $L(z)$, which converts a safety factor into an expected number of units short per cycle. total_cost(x) is the scalar objective function that scipy.optimize.minimize calls; it unpacks x = [Q, z] and sums the four cost terms described in the math section above. A guard clause returns a very large penalty if Q is non-positive, which keeps the optimizer away from degenerate values. total_cost_vec(Q, z) is a NumPy-vectorized twin of the same function — it accepts full arrays for Q and z and returns an array of costs with no Python-level loop, which is what makes the 3D surface below build instantly instead of looping over thousands of grid points one at a time.

Section 3 – Joint optimization. Rather than solving order quantity and safety stock as two independent problems, minimize searches over both Q and z at once, using the classical EOQ value and a 95%-service-level guess (z = 1.65) as the starting point. L-BFGS-B is used because it supports box constraints (bounds), which keeps Q and z in economically sensible ranges during the search. The results are then translated into safety stock and reorder point using the formulas from the math section.

Section 4 – Cost surface. A 120×120 grid of (Q, z) combinations is built with np.meshgrid, and the entire cost surface is evaluated in one vectorized call to total_cost_vec. This avoids a nested double loop (120×120 = 14,400 evaluations) and keeps the whole computation well under a second.

Section 5 – 3D plot. plot_surface draws the cost landscape, and the optimum found in Section 3 is marked as a single cyan point. The surface makes the trade-off visually obvious: moving along the $Q$-axis away from the optimum increases cost because of the ordering/holding trade-off, while moving along the $z$-axis away from the optimum increases cost because of the safety-stock/shortage trade-off. The dark theme (background, panes, tick colors) matches a typical technical-blog dark layout.

Section 6 – 2D breakdown. Fixing $z$ at its optimal value, this chart decomposes total cost into its four components as a function of $Q$ alone — this is the classic “U-shaped EOQ curve” but now shown alongside the safety-stock and shortage-cost lines so the reader can see how much of the total cost each component actually contributes at the optimum.

The 3D surface is bowl-shaped, with a single global minimum in the interior of the plotted region — that’s the optimum reported by the optimizer. Moving toward small $Q$ makes the surface rise steeply because ordering cost explodes as $D/Q$; moving toward small $z$ makes it rise because expected shortage cost grows. The optimum sits at a moderate order quantity with a fairly high safety factor, since in this example the shortage cost is significant relative to holding cost.

In the 2D breakdown, the white “Total Cost” curve is clearly the sum of the other four curves, and its minimum lines up exactly with the cyan marker. Ordering cost falls steadily as $Q$ grows, cycle holding cost rises linearly, and the safety-stock holding cost stays flat (since $z$ is fixed in this chart) while shortage cost falls as larger, more frequent… actually less frequent orders reduce the number of cycles per year in which a shortage can occur.

Interpretation of the Numbers

For this example, the model finds an optimal order quantity of roughly 526 units per order — close to, but not identical to, the classical EOQ of about 520 units, because uncertainty slightly shifts the ideal batch size. The optimal safety factor comes out to about $z \approx 2.27$, which corresponds to an implied service level of roughly 98.8%. That translates into a safety stock of about 42 units and a reorder point of about 503 units, at a minimum total annual cost of roughly $2,274.

The relatively high service level here is a direct consequence of the numbers chosen: a shortage cost of $15 per unit is considerably higher than the $4 holding cost per unit, so the optimizer prefers to carry a bit more safety stock rather than risk running short. Changing that ratio — for example, lowering the shortage cost or raising the holding cost — would pull the optimal $z$ down and reduce the safety stock accordingly. That sensitivity is exactly what the 3D surface makes visible at a glance: the cost landscape’s shape along the $z$-axis directly reflects how expensive shortages are relative to holding inventory.

Profit Maximization with a Cobb-Douglas Production Function

August 18, 2026

Optimizing Labor and Capital in Python

Every firm faces the same fundamental question: how much labor and how much capital should it employ to maximize profit? When production follows a Cobb-Douglas technology, this question turns into a clean, well-behaved optimization problem — one that’s perfect for illustrating both the economics and the numerical methods behind it. In this article, we’ll set up a two-variable profit maximization problem, solve it with a closed-form derivation and a numerical solver, and visualize the profit landscape in 3D.

The Economic Setup

Consider a firm that produces output $Q$ using labor $L$ and capital $K$ according to a Cobb-Douglas production function:

$$
Q(L, K) = A , L^{\alpha} K^{\beta}
$$

where $A$ is total factor productivity, and $\alpha, \beta \in (0,1)$ are the output elasticities of labor and capital. When $\alpha + \beta < 1$, the technology exhibits decreasing returns to scale, which guarantees a well-defined, interior profit-maximizing point (rather than a corner solution or an unbounded profit).

The firm sells output at price $p$, pays wage $w$ per unit of labor, and pays rental rate $r$ per unit of capital. Profit is:

$$
\pi(L, K) = p , A , L^{\alpha} K^{\beta} - wL - rK
$$

The firm’s problem is:

$$
\max_{L > 0,, K > 0} ; \pi(L, K)
$$

First-Order Conditions

Taking partial derivatives and setting them to zero gives the classic marginal-revenue-product-equals-input-price conditions:

$$
\frac{\partial \pi}{\partial L} = p A \alpha L^{\alpha - 1} K^{\beta} - w = 0
$$

$$
\frac{\partial \pi}{\partial K} = p A \beta L^{\alpha} K^{\beta - 1} - r = 0
$$

Dividing the first equation by the second eliminates $p$ and $A$, yielding a simple ratio between the optimal capital-labor ratio and the relative input prices:

$$
\frac{K^*}{L^*} = \frac{\beta w}{\alpha r}
$$

Substituting this back into the first FOC and solving for $L^*$ gives a closed-form solution:

$$
L^* = \left[ \frac{w}{pA\alpha \left(\dfrac{\beta w}{\alpha r}\right)^{\beta}} \right]^{\frac{1}{\alpha + \beta - 1}}, \qquad K^* = \frac{\beta w}{\alpha r} , L^*
$$

This is exactly what we’ll implement in Python, then cross-check numerically with scipy.optimize.

Concrete Numerical Example

We’ll use the following parameters:

$A = 8$ (productivity)
$\alpha = 0.35$ (labor elasticity)
$\beta = 0.25$ (capital elasticity)
$p = 10$ (output price)
$w = 6$ (wage rate)
$r = 4$ (capital rental rate)

Since $\alpha + \beta = 0.6 < 1$, the profit function is strictly concave in $(L, K)$, so it has a unique interior maximum.

Source Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# ---------- Dark theme for Colab plots ----------
plt.rcParams['figure.facecolor']  = '#1e1e1e'
plt.rcParams['axes.facecolor']    = '#1e1e1e'
plt.rcParams['savefig.facecolor'] = '#1e1e1e'
plt.rcParams['text.color']        = 'white'
plt.rcParams['axes.labelcolor']   = 'white'
plt.rcParams['xtick.color']       = 'white'
plt.rcParams['ytick.color']       = 'white'
plt.rcParams['axes.edgecolor']    = 'white'
plt.rcParams['grid.color']        = '#555555'

# ---------- Model parameters ----------
A = 8.0      # total factor productivity
alpha = 0.35 # output elasticity of labor
beta = 0.25  # output elasticity of capital
p = 10.0     # output price
w = 6.0      # wage rate
r = 4.0      # capital rental rate

# ---------- Core functions ----------
def production(L, K):
    """Cobb-Douglas production function."""
    return A * np.power(L, alpha) * np.power(K, beta)

def profit(L, K):
    """Profit as a function of labor and capital."""
    return p * production(L, K) - w * L - r * K

def neg_profit(x):
    """Negative profit for minimization (scipy minimizes by default)."""
    L, K = x
    if L <= 0 or K <= 0:
        return 1e12
    return -profit(L, K)

# ---------- Closed-form analytical solution ----------
def analytical_solution():
    ratio = (beta * w) / (alpha * r)                 # optimal K*/L*
    inner = p * A * alpha * (ratio ** beta)
    exponent = alpha + beta - 1.0
    L_star = (w / inner) ** (1.0 / exponent)
    K_star = ratio * L_star
    return L_star, K_star

L_analytical, K_analytical = analytical_solution()

# ---------- Numerical cross-check with scipy ----------
result = minimize(
    neg_profit,
    x0=[1.0, 1.0],
    method='L-BFGS-B',
    bounds=[(1e-6, None), (1e-6, None)]
)
L_numeric, K_numeric = result.x
Q_numeric = production(L_numeric, K_numeric)
profit_numeric = -result.fun

# ---------- Print results ----------
print("=== Profit Maximization: Cobb-Douglas Production ===")
print(f"Parameters: A={A}, alpha={alpha}, beta={beta}, p={p}, w={w}, r={r}")
print()
print("Closed-form analytical solution:")
print(f"  L* = {L_analytical:.6f}")
print(f"  K* = {K_analytical:.6f}")
print()
print("Numerical solution (scipy L-BFGS-B):")
print(f"  L* = {L_numeric:.6f}")
print(f"  K* = {K_numeric:.6f}")
print(f"  Q*  = {Q_numeric:.6f}")
print(f"  Profit* = {profit_numeric:.6f}")

# ============================================================
# Figure 1: 3D profit surface + contour map
# ============================================================
L_range = np.linspace(0.5, L_numeric * 3, 120)
K_range = np.linspace(0.5, K_numeric * 3, 120)
L_grid, K_grid = np.meshgrid(L_range, K_range)
Profit_grid = profit(L_grid, K_grid)

fig = plt.figure(figsize=(16, 7))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(L_grid, K_grid, Profit_grid, cmap='plasma',
                         alpha=0.9, edgecolor='none', antialiased=True)
ax1.scatter([L_numeric], [K_numeric], [profit_numeric],
            color='cyan', s=120, marker='o', depthshade=False,
            label='Optimal point (L*, K*)')
ax1.set_xlabel('Labor  L')
ax1.set_ylabel('Capital  K')
ax1.set_zlabel('Profit  π')
ax1.set_title('Profit Surface over Labor and Capital', color='white')
ax1.view_init(elev=28, azim=-135)
ax1.legend(loc='upper left')
fig.colorbar(surf, ax=ax1, shrink=0.5, pad=0.12, label='Profit')

ax2 = fig.add_subplot(1, 2, 2)
cf = ax2.contourf(L_grid, K_grid, Profit_grid, levels=40, cmap='plasma')
ax2.contour(L_grid, K_grid, Profit_grid, levels=15, colors='white',
            linewidths=0.4, alpha=0.6)
ax2.plot(L_numeric, K_numeric, 'o', color='cyan', markersize=11,
         label='Optimal point (L*, K*)')
ax2.set_xlabel('Labor  L')
ax2.set_ylabel('Capital  K')
ax2.set_title('Profit Contour Map', color='white')
ax2.legend()
fig.colorbar(cf, ax=ax2, label='Profit')

plt.tight_layout()
plt.show()

# ============================================================
# Figure 2: Profit slices through the optimum
# ============================================================
fig2, axes = plt.subplots(1, 2, figsize=(14, 5))

L_slice = np.linspace(0.5, L_numeric * 3, 200)
profit_L_slice = profit(L_slice, K_numeric)
axes[0].plot(L_slice, profit_L_slice, color='cyan', linewidth=2.2)
axes[0].axvline(L_numeric, color='magenta', linestyle='--',
                 label=f'L* = {L_numeric:.2f}')
axes[0].set_xlabel('Labor  L')
axes[0].set_ylabel('Profit  π')
axes[0].set_title('Profit vs Labor (K fixed at K*)', color='white')
axes[0].grid(alpha=0.3)
axes[0].legend()

K_slice = np.linspace(0.5, K_numeric * 3, 200)
profit_K_slice = profit(L_numeric, K_slice)
axes[1].plot(K_slice, profit_K_slice, color='orange', linewidth=2.2)
axes[1].axvline(K_numeric, color='magenta', linestyle='--',
                 label=f'K* = {K_numeric:.2f}')
axes[1].set_xlabel('Capital  K')
axes[1].set_ylabel('Profit  π')
axes[1].set_title('Profit vs Capital (L fixed at L*)', color='white')
axes[1].grid(alpha=0.3)
axes[1].legend()

plt.tight_layout()
plt.show()

# ============================================================
# Figure 3: Comparative statics — how the optimum shifts with price p
# ============================================================
p_range = np.linspace(5.0, 20.0, 60)
L_path = np.zeros_like(p_range)
K_path = np.zeros_like(p_range)
profit_path = np.zeros_like(p_range)

ratio_fixed = (beta * w) / (alpha * r)
exponent = alpha + beta - 1.0

for i, p_val in enumerate(p_range):
    inner = p_val * A * alpha * (ratio_fixed ** beta)
    L_val = (w / inner) ** (1.0 / exponent)
    K_val = ratio_fixed * L_val
    L_path[i] = L_val
    K_path[i] = K_val
    profit_path[i] = p_val * production(L_val, K_val) - w * L_val - r * K_val

fig3, axes3 = plt.subplots(1, 3, figsize=(18, 5))

axes3[0].plot(p_range, L_path, color='cyan', linewidth=2.2)
axes3[0].set_xlabel('Output price  p')
axes3[0].set_ylabel('Optimal Labor  L*')
axes3[0].set_title('L* vs Output Price', color='white')
axes3[0].grid(alpha=0.3)

axes3[1].plot(p_range, K_path, color='orange', linewidth=2.2)
axes3[1].set_xlabel('Output price  p')
axes3[1].set_ylabel('Optimal Capital  K*')
axes3[1].set_title('K* vs Output Price', color='white')
axes3[1].grid(alpha=0.3)

axes3[2].plot(p_range, profit_path, color='lime', linewidth=2.2)
axes3[2].set_xlabel('Output price  p')
axes3[2].set_ylabel('Maximum Profit  π*')
axes3[2].set_title('Optimal Profit vs Output Price', color='white')
axes3[2].grid(alpha=0.3)

plt.tight_layout()
plt.show()

=== Profit Maximization: Cobb-Douglas Production ===
Parameters: A=8.0, alpha=0.35, beta=0.25, p=10.0, w=6.0, r=4.0

Closed-form analytical solution:
  L* = 49.118372
  K* = 52.626828

Numerical solution (scipy L-BFGS-B):
  L* = 49.118528
  K* = 52.626637
  Q*  = 84.202941
  Profit* = 336.811696

Code Walkthrough

Dark theme setup. The plt.rcParams block configures every plot to use a dark background with white text, matching the blog’s visual style, and applies to all three figures automatically.

production(L, K) implements the Cobb-Douglas function $Q = AL^{\alpha}K^{\beta}$ using np.power, which works cleanly whether L and K are scalars or full meshgrid arrays — this lets the same function serve both the optimizer and the plotting code.

profit(L, K) computes $\pi = pQ - wL - rK$ directly from the production function.

neg_profit(x) wraps profit for scipy.optimize.minimize, which only performs minimization. It also guards against non-positive inputs by returning a very large penalty value, keeping the optimizer inside the economically meaningful region $L, K > 0$.

analytical_solution() implements the closed-form formula derived above. It first computes the optimal capital-labor ratio ratio = (β·w)/(α·r), then solves for $L^*$ using the derived exponent formula, and finally recovers $K^*$ from the ratio. This gives an exact answer with no iterative solver involved.

scipy.optimize.minimize cross-checks the analytical result numerically using the L-BFGS-B algorithm, which handles the box constraints ($L, K > 0$) efficiently. Starting from [1.0, 1.0], it converges to the same optimum as the closed-form solution — a good sanity check that both the math and the code agree.

Figure 1 (3D surface + contour) visualizes the entire profit landscape. The 3D surface shows profit as a dome-shaped peak — a direct consequence of decreasing returns to scale making the profit function strictly concave. The contour map on the right shows the same landscape from above, with concentric rings collapsing toward the optimal point marked in cyan.

Figure 2 (profit slices) cuts through the 3D surface along each axis, holding the other input fixed at its optimal value. Both curves are single-peaked, confirming that $L^*$ and $K^*$ are each true local maximizers along their respective directions — a visual confirmation of the first-order conditions.

Figure 3 (comparative statics) explores how the optimum responds to a change in the output price p, holding A, α, β, w, r fixed. Because the capital-labor ratio K*/L* depends only on α, β, w, r (not on p), both L* and K* scale up together as p rises — the firm expands scale but keeps its input mix constant. Maximum profit π* grows even faster than either input, since revenue rises with both price and quantity simultaneously.

Interpreting the Results

The closed-form and numerical solutions should match to several decimal places, confirming the correctness of both derivations. Economically, the result illustrates a core insight of Cobb-Douglas theory: the ratio of capital to labor at the optimum is pinned down entirely by relative factor prices and output elasticities ($\beta w / \alpha r$), while the scale of production is what responds to the output price. This decomposition — between “input mix” and “input scale” — is one of the reasons Cobb-Douglas functions remain a workhorse in microeconomics and production theory.

Wrapping Up

This example shows how a two-input profit maximization problem can be solved two independent ways — analytically via the first-order conditions, and numerically via constrained optimization — with both approaches converging to the same answer. The 3D visualization makes the concavity of the profit function tangible, while the comparative statics plot reveals how the firm’s optimal scale (but not its input mix) responds to market prices. The same framework extends naturally to more complex production technologies, multiple outputs, or additional constraints such as a fixed budget for total input spending.

Maximizing Utility with Two Consumption Goods

August 17, 2026

A Cobb-Douglas Example in Python

Consumer choice theory sits at the heart of microeconomics, and one of its cleanest applications is the classic utility maximization problem with two goods. Given a fixed budget, how should a consumer split spending between two goods to get the most satisfaction possible? In this article we’ll work through a concrete Cobb-Douglas example, solve it both analytically and numerically, and visualize the solution with indifference curves, a 3D utility surface, and a contour map.

The Setup

A consumer chooses quantities of two goods, $x$ and $y$, to maximize a Cobb-Douglas utility function:

$$
U(x, y) = x^{\alpha} y^{\beta}
$$

subject to a linear budget constraint:

$$
p_x x + p_y y = I
$$

where $p_x$ and $p_y$ are the prices of the two goods and $I$ is total income. In our example we’ll use:

$$
\alpha = 0.6, \quad \beta = 0.4, \quad p_x = 4, \quad p_y = 2, \quad I = 100
$$

Solving with Lagrange Multipliers

Setting up the Lagrangian:

$$
\mathcal{L}(x, y, \lambda) = x^{\alpha} y^{\beta} + \lambda (I - p_x x - p_y y)
$$

Taking first-order conditions and eliminating $\lambda$ gives the tangency condition where the marginal rate of substitution equals the price ratio:

$$
\frac{\alpha y}{\beta x} = \frac{p_x}{p_y}
$$

Combining this with the budget constraint yields a closed-form solution:

This closed-form result gives us a perfect benchmark to check against a numerical optimizer.

Python Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D

# ----------------------------
# Dark theme for matplotlib
# ----------------------------
plt.rcParams['figure.facecolor'] = '#1e1e2e'
plt.rcParams['axes.facecolor'] = '#1e1e2e'
plt.rcParams['axes.edgecolor'] = '#cdd6f4'
plt.rcParams['axes.labelcolor'] = '#cdd6f4'
plt.rcParams['text.color'] = '#cdd6f4'
plt.rcParams['xtick.color'] = '#cdd6f4'
plt.rcParams['ytick.color'] = '#cdd6f4'
plt.rcParams['grid.color'] = '#45475a'
plt.rcParams['grid.alpha'] = 0.3

# ----------------------------
# Parameters
# ----------------------------
alpha = 0.6
beta = 0.4
px = 4.0
py = 2.0
I = 100.0

def utility(x, y, alpha=alpha, beta=beta):
    return (x ** alpha) * (y ** beta)

def neg_utility(vars):
    x, y = vars
    if x <= 0 or y <= 0:
        return 1e10
    return -utility(x, y)

def budget_constraint(vars):
    x, y = vars
    return I - (px * x + py * y)

# ----------------------------
# Analytical (closed-form) solution
# ----------------------------
x_analytic = (alpha / (alpha + beta)) * I / px
y_analytic = (beta / (alpha + beta)) * I / py
u_analytic = utility(x_analytic, y_analytic)

# ----------------------------
# Numerical solution via SLSQP
# ----------------------------
x0 = [I / (2 * px), I / (2 * py)]
constraints = [{'type': 'eq', 'fun': budget_constraint}]
bounds = [(1e-6, I / px), (1e-6, I / py)]

result = minimize(
    neg_utility, x0,
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options={'ftol': 1e-12, 'maxiter': 500}
)
x_opt, y_opt = result.x
u_opt = utility(x_opt, y_opt)

print("=== Analytical Solution (Cobb-Douglas closed form) ===")
print(f"x* = {x_analytic:.4f}")
print(f"y* = {y_analytic:.4f}")
print(f"U* = {u_analytic:.4f}")
print()
print("=== Numerical Solution (SLSQP) ===")
print(f"x* = {x_opt:.4f}")
print(f"y* = {y_opt:.4f}")
print(f"U* = {u_opt:.4f}")
print(f"Converged: {result.success}, message: {result.message}")
print()
print("=== Spending Check ===")
print(f"Total spend (analytical): {px*x_analytic + py*y_analytic:.4f} (budget = {I})")
print(f"Total spend (numerical):  {px*x_opt + py*y_opt:.4f} (budget = {I})")

# ============================================================
# Figure 1: Indifference curves + budget line (2D)
# ============================================================
x_max = I / px
y_max = I / py

x_grid = np.linspace(0.01, x_max * 1.1, 400)
y_grid = np.linspace(0.01, y_max * 1.1, 400)
X, Y = np.meshgrid(x_grid, y_grid)
Z = utility(X, Y)

fig1, ax1 = plt.subplots(figsize=(8, 7))

u_levels = [u_opt * f for f in [0.4, 0.6, 0.8, 1.0, 1.2]]
colors = ['#89b4fa', '#74c7ec', '#94e2d5', '#f9e2af', '#fab387']
contour = ax1.contour(X, Y, Z, levels=sorted(u_levels), colors=colors, linewidths=2)
ax1.clabel(contour, inline=True, fontsize=8, fmt=lambda v: f"U={v:.1f}")

x_budget = np.linspace(0, x_max, 200)
y_budget = (I - px * x_budget) / py
ax1.plot(x_budget, y_budget, color='#f38ba8', linewidth=2.5, label='Budget line')

ax1.scatter([x_opt], [y_opt], color='#a6e3a1', s=120, zorder=5,
            edgecolor='white', linewidth=1.5, label='Optimal bundle (x*, y*)')

ax1.set_xlim(0, x_max * 1.1)
ax1.set_ylim(0, y_max * 1.1)
ax1.set_xlabel('Good X')
ax1.set_ylabel('Good Y')
ax1.set_title('Indifference Curves and Budget Line')
ax1.legend(loc='upper right', facecolor='#313244', edgecolor='#45475a')
ax1.grid(True)
plt.tight_layout()
plt.show()

# ============================================================
# Figure 2: 3D utility surface with budget constraint path
# ============================================================
fig2 = plt.figure(figsize=(9, 7))
ax2 = fig2.add_subplot(111, projection='3d')

surf = ax2.plot_surface(X, Y, Z, cmap='mako' if 'mako' in plt.colormaps() else 'viridis',
                         alpha=0.75, linewidth=0, antialiased=True)

u_budget = utility(x_budget, np.clip(y_budget, 1e-6, None))
ax2.plot(x_budget, y_budget, u_budget, color='#f38ba8', linewidth=3, label='Utility along budget line')

ax2.scatter([x_opt], [y_opt], [u_opt], color='#a6e3a1', s=80,
            edgecolor='white', linewidth=1.2, label='Optimal point')

ax2.set_xlabel('Good X')
ax2.set_ylabel('Good Y')
ax2.set_zlabel('Utility U(x, y)')
ax2.set_title('Utility Surface with Budget-Constrained Path')
ax2.view_init(elev=28, azim=-50)
fig2.colorbar(surf, shrink=0.5, aspect=12, pad=0.1)
ax2.legend(loc='upper left', facecolor='#313244', edgecolor='#45475a')
plt.tight_layout()
plt.show()

# ============================================================
# Figure 3: Filled contour map with optimum highlighted
# ============================================================
fig3, ax3 = plt.subplots(figsize=(8, 7))

filled = ax3.contourf(X, Y, Z, levels=30, cmap='mako' if 'mako' in plt.colormaps() else 'viridis')
ax3.plot(x_budget, y_budget, color='#f38ba8', linewidth=2.5, label='Budget line')
ax3.scatter([x_opt], [y_opt], color='#a6e3a1', s=120, zorder=5,
            edgecolor='white', linewidth=1.5, label='Optimal bundle')

ax3.set_xlim(0, x_max * 1.1)
ax3.set_ylim(0, y_max * 1.1)
ax3.set_xlabel('Good X')
ax3.set_ylabel('Good Y')
ax3.set_title('Utility Landscape (Filled Contour)')
fig3.colorbar(filled, ax=ax3, shrink=0.85, label='Utility')
ax3.legend(loc='upper right', facecolor='#313244', edgecolor='#45475a')
plt.tight_layout()
plt.show()

=== Analytical Solution (Cobb-Douglas closed form) ===
x* = 15.0000
y* = 20.0000
U* = 16.8293

=== Numerical Solution (SLSQP) ===
x* = 15.0000
y* = 20.0000
U* = 16.8293
Converged: True, message: Optimization terminated successfully

=== Spending Check ===
Total spend (analytical): 100.0000 (budget = 100.0)
Total spend (numerical):  100.0000 (budget = 100.0)

Code Walkthrough

Utility and constraint functions. utility(x, y) implements the Cobb-Douglas form $x^{\alpha}y^{\beta}$ directly. neg_utility wraps it with a sign flip because scipy.optimize.minimize only minimizes — maximizing utility is equivalent to minimizing its negative. It also guards against non-positive quantities by returning a huge penalty value, which keeps the optimizer away from invalid regions without needing complicated bound logic.

Analytical solution. Because Cobb-Douglas preferences have a well-known closed-form solution, we compute x_analytic and y_analytic directly from the formula derived above. This isn’t just for display — it acts as a ground-truth check against the numerical result.

Numerical solution. We use scipy.optimize.minimize with the SLSQP (Sequential Least Squares Programming) method, which handles equality constraints natively. The budget constraint is passed as a dictionary with 'type': 'eq', and bounds keep both goods within a sensible positive range. Starting from a naive 50/50 budget split (x0), SLSQP converges to the same optimum as the analytical formula, which is a nice sanity check that the numerical approach is correctly specified.

Why this is already fast. This is a small, smooth, twice-differentiable convex optimization problem in two variables — SLSQP converges in a handful of iterations, so there’s no need for any special acceleration here. The computational bottleneck, if any, is in the plotting: the 3D surface and contour plots use a $400 \times 400$ grid, which is dense enough for smooth-looking curves while still rendering instantly.

Visualizing the Solution

Figure 1 — Indifference curves and the budget line. Each colored curve traces a set of $(x, y)$ bundles that give the same utility level. The steepness of these curves at any point reflects the marginal rate of substitution — how much of $y$ the consumer is willing to give up for one more unit of $x$ while staying equally satisfied. The red line is the budget constraint: every point on it costs exactly $I$. The green dot marks the optimum, and geometrically it’s exactly where the budget line is tangent to the highest reachable indifference curve. Any point further out on that indifference curve isn’t affordable, and any affordable point not on that curve leaves utility on the table.

Figure 2 — The 3D utility surface. This lifts the whole picture into three dimensions: height now directly represents utility $U(x, y)$ instead of being encoded as contour lines. The red curve traces the utility value along every affordable bundle on the budget line, and the green marker sits at its peak. Rotating this surface makes it visually obvious that the constrained problem is really about finding the highest point reachable while walking along that one path defined by the budget constraint — the unconstrained peak of the full surface would require unlimited income.

Figure 3 — Filled contour map. This is a top-down view of the same surface, using color intensity to encode utility instead of height. It makes the “climbing” intuition especially clear: the optimal bundle sits at the point on the budget line where the colors are most intense, i.e. the warmest reachable region.

Interpreting the Result

With $\alpha = 0.6$ and $\beta = 0.4$, the consumer values good $X$ relatively more, so the closed-form solution allocates a larger budget share to $X$: specifically a share of $\frac{\alpha}{\alpha+\beta} = 0.6$ of income goes to $X$ and $\frac{\beta}{\alpha+\beta} = 0.4$ goes to $Y$. This is a defining feature of Cobb-Douglas preferences — the optimal expenditure shares depend only on the exponents $\alpha$ and $\beta$, not on prices or income at all. That’s why both the analytical and numerical methods land on the same answer regardless of how the starting guess for the optimizer is chosen.

This framework generalizes readily: swapping in a CES or quasi-linear utility function, adding more goods, or introducing non-linear budget constraints (like quantity discounts) all fit naturally into the same Lagrangian and scipy.optimize machinery used here.

Minimizing Hinge Loss with L2 Regularization

August 16, 2026

A Hands-On Look at the Soft-Margin SVM

Support Vector Machines are often introduced through the lens of the quadratic programming dual problem, but the primal formulation tells a much more intuitive story: it’s just an unconstrained (well, almost) optimization problem where we’re minimizing a hinge loss term balanced against an L2 regularization term. Today we’ll build this from scratch, watch gradient descent carve out an optimal decision boundary, and visualize the loss landscape in 3D.

The Objective Function

The soft-margin SVM primal objective for a binary classification problem with labels $y_i \in {-1, +1}$ is:

$$
\mathcal{L}(w, b) = \frac{1}{2}|w|^2 + C \sum_{i=1}^{n} \max\left(0,\ 1 - y_i(w^\top x_i + b)\right)
$$

Here:

The first term, $\frac{1}{2}|w|^2$, is the L2 regularizer — it penalizes large weight vectors, which corresponds to maximizing the margin $\frac{2}{|w|}$ between the two classes.
The second term is the hinge loss, which only penalizes points that are either misclassified or sitting inside the margin. Points correctly classified with enough margin contribute zero loss.
$C$ controls the trade-off between margin width and classification error tolerance. Large $C$ pushes toward fewer margin violations (harder margin); small $C$ favors a wider margin at the cost of some misclassifications.

Since the hinge loss $\max(0, 1 - z)$ is not differentiable at $z = 1$, we can’t use plain gradient descent — but we can use subgradient descent, which picks a valid subgradient wherever the function is non-smooth. The subgradient of the objective with respect to $w$ and $b$ is:

$$
\nabla_w \mathcal{L} = w - C \sum_{i \in \mathcal{V}} y_i x_i, \qquad
\nabla_b \mathcal{L} = -C \sum_{i \in \mathcal{V}} y_i
$$

where $\mathcal{V} = { i : y_i(w^\top x_i + b) < 1 }$ is the set of points violating the margin.

The Example Problem

We’ll generate a 2D dataset of two overlapping Gaussian blobs, which forces the optimizer to genuinely trade off margin width against a handful of unavoidable violations — a much more interesting case than perfectly separable data.

Source Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets import make_blobs

# ---------------------------------------------------------
# Reproducibility
# ---------------------------------------------------------
np.random.seed(42)

# ---------------------------------------------------------
# 1. Generate a 2D binary classification dataset
#    Two overlapping blobs -> not perfectly separable
# ---------------------------------------------------------
X, y_raw = make_blobs(n_samples=200, centers=2, cluster_std=1.8, random_state=7)
y = np.where(y_raw == 0, -1, 1).astype(np.float64)  # labels in {-1, +1}

# Standardize features (helps gradient descent convergence)
X_mean, X_std = X.mean(axis=0), X.std(axis=0)
X = (X - X_mean) / X_std

n_samples, n_features = X.shape

# ---------------------------------------------------------
# 2. Vectorized hinge loss + L2 regularization
# ---------------------------------------------------------
def compute_loss(w, b, X, y, C):
    margins = y * (X @ w + b)
    hinge = np.maximum(0.0, 1.0 - margins)
    reg_term = 0.5 * np.dot(w, w)
    return reg_term + C * np.sum(hinge)

def compute_subgradient(w, b, X, y, C):
    margins = y * (X @ w + b)
    violating = margins < 1.0          # boolean mask, fully vectorized
    # Sum over violating points without any explicit Python loop
    grad_w = w - C * (y[violating, None] * X[violating]).sum(axis=0)
    grad_b = -C * y[violating].sum()
    return grad_w, grad_b

# ---------------------------------------------------------
# 3. Subgradient descent with momentum for faster convergence
# ---------------------------------------------------------
def train_svm(X, y, C=1.0, lr=0.05, momentum=0.9, n_iters=2000):
    w = np.zeros(X.shape[1])
    b = 0.0
    v_w = np.zeros_like(w)
    v_b = 0.0

    loss_history = np.zeros(n_iters)
    w_history = np.zeros((n_iters, 2))  # for 3D trajectory plot

    for t in range(n_iters):
        grad_w, grad_b = compute_subgradient(w, b, X, y, C)

        # Momentum-accelerated update (much faster than vanilla GD)
        v_w = momentum * v_w - lr * grad_w
        v_b = momentum * v_b - lr * grad_b
        w = w + v_w
        b = b + v_b

        loss_history[t] = compute_loss(w, b, X, y, C)
        w_history[t] = w

    return w, b, loss_history, w_history

C_value = 1.0
w_opt, b_opt, loss_history, w_history = train_svm(
    X, y, C=C_value, lr=0.05, momentum=0.9, n_iters=2000
)

final_loss = loss_history[-1]
margins_final = y * (X @ w_opt + b_opt)
n_violations = int(np.sum(margins_final < 1.0))
n_misclassified = int(np.sum(margins_final < 0.0))
margin_width = 2.0 / np.linalg.norm(w_opt)

print(f"Optimal w: {w_opt}")
print(f"Optimal b: {b_opt:.4f}")
print(f"Final objective value: {final_loss:.4f}")
print(f"Margin width (2/||w||): {margin_width:.4f}")
print(f"Number of margin violations: {n_violations} / {n_samples}")
print(f"Number of misclassified points: {n_misclassified} / {n_samples}")

# ▼コンソール出力ここに挿入▼

# ---------------------------------------------------------
# 4. Plot 1: Decision boundary + margins on the dataset
# ---------------------------------------------------------
plt.style.use('dark_background')
fig1, ax1 = plt.subplots(figsize=(8, 7))

x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max, 300),
                        np.linspace(x2_min, x2_max, 300))
grid_points = np.c_[xx1.ravel(), xx2.ravel()]
decision_vals = (grid_points @ w_opt + b_opt).reshape(xx1.shape)

ax1.contourf(xx1, xx2, decision_vals, levels=np.linspace(decision_vals.min(), decision_vals.max(), 25),
             cmap='coolwarm', alpha=0.35)
ax1.contour(xx1, xx2, decision_vals, levels=[-1, 0, 1],
            colors=['#00e5ff', '#ffffff', '#ff4081'], linewidths=[1.5, 2.5, 1.5],
            linestyles=['dashed', 'solid', 'dashed'])

colors = np.where(y == 1, '#ff4081', '#00e5ff')
ax1.scatter(X[:, 0], X[:, 1], c=colors, edgecolors='white', linewidths=0.5, s=45, zorder=3)

violating_mask = margins_final < 1.0
ax1.scatter(X[violating_mask, 0], X[violating_mask, 1],
            facecolors='none', edgecolors='yellow', s=140, linewidths=1.5,
            label='Margin violations', zorder=4)

ax1.set_title(f'Soft-Margin SVM Decision Boundary (C={C_value})', fontsize=14, color='white')
ax1.set_xlabel('Feature 1 (standardized)')
ax1.set_ylabel('Feature 2 (standardized)')
ax1.legend(loc='upper right', facecolor='#222222', edgecolor='gray')
plt.tight_layout()
plt.show()

# ▼決定境界と余白の可視化グラフをここに挿入▼

# ---------------------------------------------------------
# 5. Plot 2: Loss convergence curve
# ---------------------------------------------------------
fig2, ax2 = plt.subplots(figsize=(8, 5))
ax2.plot(loss_history, color='#00e5ff', linewidth=2)
ax2.set_title('Objective Value During Subgradient Descent', fontsize=14, color='white')
ax2.set_xlabel('Iteration')
ax2.set_ylabel(r'$\mathcal{L}(w,b)$')
ax2.grid(alpha=0.2)
plt.tight_layout()
plt.show()

# ▼損失関数の収束グラフをここに挿入▼

# ---------------------------------------------------------
# 6. Plot 3: 3D loss surface over (w1, w2) with descent trajectory
#    b is fixed at its optimal value for this visualization
# ---------------------------------------------------------
w1_range = np.linspace(w_opt[0] - 3, w_opt[0] + 3, 80)
w2_range = np.linspace(w_opt[1] - 3, w_opt[1] + 3, 80)
W1, W2 = np.meshgrid(w1_range, w2_range)

# Vectorized surface computation (no nested Python loops)
margins_grid = y[None, None, :] * (
    X[None, None, :, 0] * W1[:, :, None] + X[None, None, :, 1] * W2[:, :, None] + b_opt
)
hinge_grid = np.maximum(0.0, 1.0 - margins_grid)
Loss_surface = 0.5 * (W1**2 + W2**2) + C_value * hinge_grid.sum(axis=2)

fig3 = plt.figure(figsize=(9, 7))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(W1, W2, Loss_surface, cmap='plasma', alpha=0.85,
                         linewidth=0, antialiased=True)

# Overlay the trajectory taken by subgradient descent
traj_loss = np.array([compute_loss(w_history[t], b_opt, X, y, C_value)
                       for t in range(0, len(w_history), 20)])
traj_w1 = w_history[::20, 0]
traj_w2 = w_history[::20, 1]
ax3.plot(traj_w1, traj_w2, traj_loss, color='#00ff88', linewidth=2.5,
          marker='o', markersize=2, label='Descent path')
ax3.scatter([w_opt[0]], [w_opt[1]], [final_loss], color='white', s=80,
            edgecolors='#00ff88', linewidths=2, label='Optimum', zorder=5)

ax3.set_title('Loss Landscape over (w1, w2) with Descent Trajectory', fontsize=13, color='white')
ax3.set_xlabel('w1')
ax3.set_ylabel('w2')
ax3.set_zlabel(r'$\mathcal{L}(w,b)$')
fig3.colorbar(surf, shrink=0.5, aspect=12, pad=0.1)
ax3.legend(loc='upper left', facecolor='#222222', edgecolor='gray')
plt.tight_layout()
plt.show()

Console Output

Optimal w: [ 2.32509315 -0.46864248]
Optimal b: -0.0653
Final objective value: 17.9965
Margin width (2/||w||): 0.8432
Number of margin violations: 22 / 200
Number of misclassified points: 5 / 200

Code Walkthrough

Data generation and preprocessing. We use make_blobs with a fairly large cluster_std=1.8 so the two classes overlap somewhat — this is what makes soft-margin behavior (as opposed to hard-margin) actually necessary. Standardizing the features to zero mean and unit variance is important here: gradient-based optimization of the SVM objective converges much faster and more reliably on standardized inputs, since the regularization term $\frac{1}{2}|w|^2$ and the hinge term operate on comparable scales.

compute_loss. This directly implements $\mathcal{L}(w,b)$ from the formula above. The margin for every sample, $y_i(w^\top x_i + b)$, is computed in a single matrix-vector product X @ w, then combined with the hinge via np.maximum(0.0, 1.0 - margins) — no loop over individual samples.

compute_subgradient. This is the core of the optimization. Rather than looping over each sample and checking if margin < 1, we build a boolean mask violating over the entire array at once. Then X[violating] selects only the rows corresponding to margin-violating points, and (y[violating, None] * X[violating]).sum(axis=0) computes $\sum_{i \in \mathcal{V}} y_i x_i$ as a single reduction. This vectorization is what keeps the whole training loop fast even though we run 2000 iterations — there’s no per-sample Python-level iteration anywhere in the hot path.

train_svm. We use subgradient descent with momentum (momentum=0.9), which accumulates a velocity vector v_w/v_b instead of stepping directly along the raw subgradient. This significantly speeds up convergence compared to vanilla subgradient descent, because it dampens the oscillation that hinge-loss subgradients tend to cause near the optimum (since the active set $\mathcal{V}$ can flip abruptly from iteration to iteration). We record loss_history and w_history at every step so we can later visualize both convergence and the optimization path in weight space.

Post-training diagnostics. After training, we recompute the final margins to count how many points are margin violators (margin < 1) versus actually misclassified (margin < 0) — these are different things: a point can sit inside the margin and still be correctly classified. We also report the margin width $2/|w|$, which is the geometric quantity the regularization term is implicitly maximizing.

Visualizing the Results

Plot 1 — Decision boundary. The solid white line is the decision boundary $w^\top x + b = 0$; the dashed cyan and pink lines are the margin boundaries $w^\top x + b = \pm 1$. Points circled in yellow are margin violators — some of these are still correctly classified but sit too close to the boundary, while others have crossed to the wrong side entirely. The background shading shows the signed distance field.

Plot 2 — Convergence curve. This tracks $\mathcal{L}(w,b)$ over all 2000 iterations. Because momentum is used, the curve typically drops sharply in the first few hundred iterations and then flattens as the active constraint set stabilizes.

Plot 3 — 3D loss landscape. This is the most illuminating visualization: we fix $b$ at its converged value and sweep $w_1, w_2$ over a grid to render $\mathcal{L}(w_1, w_2, b^*)$ as a 3D surface. Because hinge loss is piecewise-linear and the regularizer is quadratic, the surface is convex but has visible creases where individual hinge terms switch from active to inactive. The green trajectory shows the actual path subgradient descent took through this landscape, sampled every 20 iterations, ending at the white marker — the converged optimum.

Interpreting the Trade-off

Try re-running the training with different values of C_value. Increasing C (e.g. to 10.0) will shrink the margin width and reduce the number of violations, since misclassification becomes more costly relative to margin size. Decreasing C (e.g. to 0.1) will widen the margin substantially and tolerate more violations — the regularization term starts to dominate. This single hyperparameter is the entire story of the bias-variance trade-off for SVMs, and watching the 3D surface’s minimum shift as C changes is a great way to build intuition for it.

Minimizing the Negative Log-Likelihood of a Bivariate Gaussian Mixture Model

August 15, 2026

A Full Walkthrough in Python

Gaussian Mixture Models (GMMs) are one of the most elegant tools in probabilistic machine learning: they let us describe a complex, multi-modal cloud of data as a weighted sum of simple Gaussian “blobs.” Fitting a GMM comes down to one core task — minimizing the negative log-likelihood (NLL) of the data under the mixture model. In this post, we’ll build a complete, runnable example in Python (Google Colaboratory) that generates synthetic two-dimensional data from a known mixture, fits a GMM by directly minimizing the NLL, and visualizes the result with both 2D contour plots and a 3D density surface.

1. The Math Behind a Gaussian Mixture Model

A bivariate ($D=2$) Gaussian Mixture Model with $K$ components describes each data point $\mathbf{x} \in \mathbb{R}^2$ as being drawn from a weighted sum of Gaussian densities:

$$
p(\mathbf{x} \mid \Theta) = \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \Sigma_k)
$$

where $\pi_k$ are the mixing weights ($\sum_k \pi_k = 1$, $\pi_k \geq 0$), $\boldsymbol{\mu}_k \in \mathbb{R}^2$ is the mean of component $k$, and $\Sigma_k$ is its $2 \times 2$ covariance matrix. The multivariate normal density itself is:

$$
\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \Sigma) = \frac{1}{2\pi \sqrt{|\Sigma|}} \exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu})\right)
$$

Given $N$ i.i.d. data points ${\mathbf{x}_1, \dots, \mathbf{x}_N}$, the log-likelihood of the entire dataset is:

$$
\log \mathcal{L}(\Theta) = \sum_{i=1}^{N} \log \left( \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \Sigma_k) \right)
$$

Fitting the model means finding $\Theta = {\pi_k, \boldsymbol{\mu}_k, \Sigma_k}$ that minimizes the negative log-likelihood:

$$
\text{NLL}(\Theta) = -\sum_{i=1}^{N} \log \left( \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \Sigma_k) \right)
$$

This is usually solved with the EM algorithm, but it can also be solved as a direct numerical optimization problem — which is exactly what we’ll do here, since it makes the “NLL minimization” framing explicit and lets us plot its convergence curve.

The tricky part of direct optimization is that $\pi_k$ must sum to 1 and $\Sigma_k$ must be positive-definite. We handle this with two standard reparameterization tricks:

Weights: parametrize $K-1$ free logits and pass them through a softmax, guaranteeing $\pi_k \geq 0$ and $\sum_k \pi_k = 1$.
Covariances: parametrize each $\Sigma_k$ via its Cholesky factor $L_k$ (a lower-triangular matrix with positive diagonal, enforced via $\exp(\cdot)$), so that $\Sigma_k = L_k L_k^\top$ is automatically positive-definite.

2. The Example We’ll Solve

We generate 600 synthetic points from a true 3-component bivariate Gaussian mixture with known weights, means, and covariances, then pretend we don’t know the true parameters and recover them by minimizing the NLL with scipy.optimize.minimize.

3. Full Python Source Code (Google Colaboratory)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)
from scipy.optimize import minimize
from scipy.stats import multivariate_normal

# ----------------------------------------------------------------------
# 1. Reproducibility
# ----------------------------------------------------------------------
np.random.seed(42)

# ----------------------------------------------------------------------
# 2. Generate synthetic 2D data from a true 3-component Gaussian mixture
# ----------------------------------------------------------------------
true_weights = np.array([0.35, 0.25, 0.40])
true_means = np.array([
    [0.0, 0.0],
    [5.0, 5.0],
    [0.0, 6.0],
])
true_covs = np.array([
    [[1.0, 0.3], [0.3, 0.8]],
    [[1.2, -0.4], [-0.4, 1.0]],
    [[0.6, 0.0], [0.0, 1.5]],
])

n_samples = 600
K_true = len(true_weights)
component_choice = np.random.choice(K_true, size=n_samples, p=true_weights)
X = np.zeros((n_samples, 2))
for k in range(K_true):
    idx = component_choice == k
    n_k = idx.sum()
    X[idx] = np.random.multivariate_normal(true_means[k], true_covs[k], size=n_k)

# ----------------------------------------------------------------------
# 3. Parametrize the GMM so the optimizer only sees unconstrained reals
# ----------------------------------------------------------------------
K = 3   # number of components we fit
D = 2   # dimensionality

def unpack_params(theta, K=K, D=D):
    idx = 0

    # --- mixture weights via softmax of (K-1) free logits ---
    logits = np.concatenate([theta[idx:idx + (K - 1)], [0.0]])
    idx += (K - 1)
    weights = np.exp(logits - logits.max())
    weights /= weights.sum()

    # --- means ---
    means = theta[idx:idx + K * D].reshape(K, D)
    idx += K * D

    # --- covariances via Cholesky factors (guarantees PD covariance) ---
    covs = np.zeros((K, D, D))
    for k in range(K):
        l11 = np.exp(theta[idx]);  idx += 1
        l21 = theta[idx];          idx += 1
        l22 = np.exp(theta[idx]);  idx += 1
        L = np.array([[l11, 0.0], [l21, l22]])
        covs[k] = L @ L.T

    return weights, means, covs

n_params = (K - 1) + K * D + K * 3

def negative_log_likelihood(theta, X):
    weights, means, covs = unpack_params(theta)
    n = X.shape[0]
    component_pdf = np.zeros((n, K))
    for k in range(K):
        component_pdf[:, k] = weights[k] * multivariate_normal.pdf(
            X, mean=means[k], cov=covs[k]
        )
    mixture_pdf = component_pdf.sum(axis=1)
    mixture_pdf = np.clip(mixture_pdf, 1e-300, None)
    return -np.sum(np.log(mixture_pdf))

# ----------------------------------------------------------------------
# 4. Initialize parameters and run the optimizer
# ----------------------------------------------------------------------
rng = np.random.default_rng(0)
theta0 = np.zeros(n_params)

idx = 0
theta0[idx:idx + (K - 1)] = 0.0
idx += (K - 1)

init_means = X[rng.choice(n_samples, size=K, replace=False)]
theta0[idx:idx + K * D] = init_means.flatten()
idx += K * D

data_std = X.std(axis=0).mean()
for k in range(K):
    theta0[idx] = np.log(data_std); idx += 1
    theta0[idx] = 0.0;              idx += 1
    theta0[idx] = np.log(data_std); idx += 1

nll_history = []
def callback(theta):
    nll_history.append(negative_log_likelihood(theta, X))

result = minimize(
    negative_log_likelihood,
    theta0,
    args=(X,),
    method="BFGS",
    callback=callback,
    options={"maxiter": 500, "disp": False},
)

fitted_weights, fitted_means, fitted_covs = unpack_params(result.x)

print("Converged:", result.success)
print("Final negative log-likelihood:", result.fun)
print("Fitted weights:", np.round(fitted_weights, 3))
print("Fitted means:\n", np.round(fitted_means, 3))

# ----------------------------------------------------------------------
# 5. Assign each point to its most likely component (responsibilities)
# ----------------------------------------------------------------------
resp = np.zeros((n_samples, K))
for k in range(K):
    resp[:, k] = fitted_weights[k] * multivariate_normal.pdf(
        X, mean=fitted_means[k], cov=fitted_covs[k]
    )
resp /= resp.sum(axis=1, keepdims=True)
labels = resp.argmax(axis=1)

# ----------------------------------------------------------------------
# 6. Build a grid for contour / 3D surface plotting
# ----------------------------------------------------------------------
x_min, x_max = X[:, 0].min() - 2, X[:, 0].max() + 2
y_min, y_max = X[:, 1].min() - 2, X[:, 1].max() + 2
xx, yy = np.meshgrid(
    np.linspace(x_min, x_max, 150),
    np.linspace(y_min, y_max, 150),
)
grid_points = np.column_stack([xx.ravel(), yy.ravel()])

density = np.zeros(grid_points.shape[0])
for k in range(K):
    density += fitted_weights[k] * multivariate_normal.pdf(
        grid_points, mean=fitted_means[k], cov=fitted_covs[k]
    )
density = density.reshape(xx.shape)

# ----------------------------------------------------------------------
# 7. Plot 1: data scatter + fitted contours
# ----------------------------------------------------------------------
fig1, ax1 = plt.subplots(figsize=(7, 6))
ax1.scatter(X[:, 0], X[:, 1], c=labels, cmap="viridis", s=15, alpha=0.7)
ax1.contour(xx, yy, density, levels=10, cmap="Reds")
ax1.scatter(fitted_means[:, 0], fitted_means[:, 1], c="black", marker="x",
            s=120, linewidths=3, label="Fitted centers")
ax1.set_xlabel("x1")
ax1.set_ylabel("x2")
ax1.set_title("Fitted 2D Gaussian Mixture Model (contours) over data")
ax1.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------------------
# 8. Plot 2: 3D surface of the fitted density
# ----------------------------------------------------------------------
fig2 = plt.figure(figsize=(8, 6))
ax2 = fig2.add_subplot(111, projection="3d")
ax2.plot_surface(xx, yy, density, cmap="viridis", linewidth=0, antialiased=True)
ax2.set_xlabel("x1")
ax2.set_ylabel("x2")
ax2.set_zlabel("Probability density")
ax2.set_title("3D Surface of the Fitted GMM Density")
plt.tight_layout()
plt.show()

# ----------------------------------------------------------------------
# 9. Plot 3: negative log-likelihood convergence
# ----------------------------------------------------------------------
fig3, ax3 = plt.subplots(figsize=(7, 5))
ax3.plot(nll_history, marker="o", markersize=3)
ax3.set_xlabel("Optimizer iteration")
ax3.set_ylabel("Negative log-likelihood")
ax3.set_title("Convergence of the Negative Log-Likelihood")
ax3.grid(alpha=0.3)
plt.tight_layout()
plt.show()

Console Output

Converged: False
Final negative log-likelihood: 2282.7653343388647
Fitted weights: [0.415 0.233 0.352]
Fitted means:
 [[0.078 6.027]
 [4.726 5.235]
 [0.03  0.024]]

4. Code Walkthrough

Sections 1–2 (data generation): We fix a random seed for reproducibility, then define three “ground truth” bivariate Gaussians with different means, covariances, and mixing weights. np.random.choice decides which component each of the 600 points belongs to, and np.random.multivariate_normal draws the actual samples. This is our synthetic dataset — in a real project, X would simply be your observed 2D data.

Section 3 (parametrization): This is the heart of the “NLL minimization” approach. Instead of optimizing $\pi_k$, $\boldsymbol{\mu}_k$, $\Sigma_k$ directly (which have constraints), we optimize a single unconstrained vector theta. unpack_params converts theta back into valid weights, means, and covariances every time it’s called:

The softmax trick turns $K-1$ free numbers into $K$ probabilities that sum to 1.
The Cholesky trick turns 3 free numbers per component into a valid $2\times 2$ positive-definite covariance matrix, since any matrix of the form $LL^\top$ is guaranteed to be positive semi-definite when $L$ has positive diagonal entries.

negative_log_likelihood implements the NLL formula from Section 1: for each component we compute the weighted density at every data point via scipy.stats.multivariate_normal.pdf, sum across components to get the mixture density per point, then sum the log of that (with a small clip to avoid log(0)).

Section 4 (optimization): We initialize the means by randomly picking 3 actual data points (a common, effective heuristic) and initialize each covariance to be a scaled identity matrix based on the data’s overall spread. scipy.optimize.minimize with the BFGS method then searches for the theta that minimizes the NLL, using a callback to record the NLL after every iteration for later plotting.

Section 5 (responsibilities): Once fitted, we compute the posterior probability that each point belongs to each component (its “responsibility”), and assign each point to its most probable component via argmax. This gives us cluster labels for coloring the scatter plot.

Sections 6–9 (plotting): We build a fine grid over the data range, evaluate the fitted mixture density on that grid, and produce three plots: a 2D contour plot over the scattered data, a 3D surface of the density function, and the NLL convergence curve.

5. Why Vectorization Matters Here

A naive implementation of negative_log_likelihood might loop over every data point and every component with plain Python for loops, calling the Gaussian density formula one point at a time. For $N=600$ points, $K=3$ components, and an optimizer that might call the objective function hundreds of times (once per BFGS iteration, plus extra calls for numerical gradient estimation), that naive approach would execute the density formula on the order of hundreds of thousands to millions of times in pure Python — noticeably slow, and potentially the difference between a cell that finishes instantly and one that hangs for a long time.

The code above avoids this entirely by calling multivariate_normal.pdf(X, mean=..., cov=...) once per component, letting SciPy evaluate the density for all $N$ points simultaneously using vectorized, compiled NumPy operations under the hood. This reduces the inner loop from $N \times K$ scalar Python operations to just $K$ vectorized calls, which is what makes this example fast and reliable even inside an iterative optimizer.

6. Visualizing the Results

The first plot overlays the raw data (colored by which fitted component most likely generated each point) with red contour lines showing the fitted mixture density, and black X markers at the three fitted component centers. This is the most direct way to see whether the GMM has correctly identified the three underlying blobs.

The second plot renders the same fitted density as a full 3D surface, where height represents probability density. The three “peaks” correspond to the three Gaussian components, and you can visually inspect how their shapes (steepness, orientation, spread) reflect the fitted covariance matrices — an elongated, tilted peak indicates correlation between the two variables, while a symmetric peak indicates roughly independent variables.

The third plot tracks the negative log-likelihood at every optimizer iteration. Because we are minimizing the NLL, this curve should decrease monotonically (or nearly so) and flatten out as the optimizer converges — a flattening curve is a good visual confirmation that BFGS successfully found a local minimum of the NLL surface.

Conclusion

We’ve walked through the full pipeline of fitting a bivariate Gaussian Mixture Model by directly minimizing its negative log-likelihood: reparameterizing constrained parameters into an unconstrained optimization space, vectorizing the likelihood computation for speed, running a quasi-Newton optimizer with convergence tracking, and visualizing the fitted density both in 2D and 3D. This direct-optimization approach is a great complement to the more commonly taught EM algorithm, and it generalizes naturally to more components, higher dimensions, or custom priors on the parameters.

Maximizing the Log-Likelihood in Logistic Regression

August 14, 2026

A Hands-On Example with Gradient Ascent vs. Newton-Raphson

Logistic regression is one of the most widely used models for binary classification, and at its core lies a beautiful optimization problem: finding the parameter vector that maximizes the log-likelihood of the observed data. In this article, we’ll build a concrete example from scratch — a synthetic medical diagnosis dataset — and solve the maximum likelihood estimation (MLE) problem using two different optimization strategies: plain gradient ascent and the much faster Newton-Raphson (IRLS) method.

The Problem Setup

Suppose we have $n$ observations, each with a feature vector $\mathbf{x}_i \in \mathbb{R}^p$ and a binary label $y_i \in {0, 1}$. Logistic regression models the probability of the positive class as:

$$P(y_i = 1 \mid \mathbf{x}_i) = \sigma(\mathbf{x}_i^\top \boldsymbol{\beta}) = \frac{1}{1 + e^{-\mathbf{x}_i^\top \boldsymbol{\beta}}}$$

where $\boldsymbol{\beta}$ is the parameter vector we want to estimate (including an intercept term).

The Log-Likelihood Function

Assuming the observations are independent, the likelihood of the entire dataset is:

$$L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \sigma(\mathbf{x}_i^\top \boldsymbol{\beta})^{y_i} \left(1 - \sigma(\mathbf{x}_i^\top \boldsymbol{\beta})\right)^{1 - y_i}$$

Taking the logarithm turns this product into a sum, which is far easier to optimize:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ y_i \log \sigma(\mathbf{x}_i^\top \boldsymbol{\beta}) + (1 - y_i) \log\left(1 - \sigma(\mathbf{x}_i^\top \boldsymbol{\beta})\right) \right]$$

Our goal is:

$$\boldsymbol{\beta}^{*} = \underset{\boldsymbol{\beta}}{\arg\max}\ \ell(\boldsymbol{\beta})$$

Since $\ell(\boldsymbol{\beta})$ is concave in $\boldsymbol{\beta}$, this problem has a unique global maximum, which makes it a perfect candidate for gradient-based optimization.

Gradient of the Log-Likelihood

Differentiating $\ell(\boldsymbol{\beta})$ with respect to $\boldsymbol{\beta}$ gives a remarkably clean expression:

$$\nabla \ell(\boldsymbol{\beta}) = \mathbf{X}^\top (\mathbf{y} - \boldsymbol{\sigma})$$

where $\boldsymbol{\sigma} = \sigma(\mathbf{X}\boldsymbol{\beta})$ is the vector of predicted probabilities. This gradient tells us how to move $\boldsymbol{\beta}$ to increase the log-likelihood, and it’s the basis of gradient ascent:

$$\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + \eta \nabla \ell(\boldsymbol{\beta}^{(t)})$$

Gradient ascent is simple, but it can take hundreds or thousands of iterations to converge, especially when features are on different scales or the log-likelihood surface is elongated (ill-conditioned).

Speeding Things Up: Newton-Raphson / IRLS

To converge dramatically faster, we can use second-order information — the Hessian of the log-likelihood:

$$H(\boldsymbol{\beta}) = -\mathbf{X}^\top \mathbf{W} \mathbf{X}, \qquad \mathbf{W} = \mathrm{diag}\big(\sigma_i(1-\sigma_i)\big)$$

The Newton-Raphson update, also known as Iteratively Reweighted Least Squares (IRLS) in the context of logistic regression, is:

$$\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - H(\boldsymbol{\beta}^{(t)})^{-1} \nabla \ell(\boldsymbol{\beta}^{(t)})$$

Because it uses curvature information, Newton-Raphson typically converges in fewer than 10 iterations, compared to hundreds for plain gradient ascent — a huge speedup when the dataset or the number of features grows.

The Example: A Synthetic Tumor Diagnosis Dataset

We’ll generate a synthetic dataset with two features — “tumor size” and “cell irregularity score” — and a binary label indicating malignant (1) or benign (0). We’ll then fit logistic regression using both gradient ascent and Newton-Raphson, compare their convergence speed, and visualize the log-likelihood landscape and the resulting decision boundary in 3D.

Full Source Code

# ==========================================================
# Logistic Regression via Log-Likelihood Maximization
# Gradient Ascent vs. Newton-Raphson (IRLS)
# ==========================================================

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# ----------------------------------------------------------
# 0. Global settings
# ----------------------------------------------------------
np.random.seed(42)
plt.style.use('dark_background')

DARK_BG = '#0d1117'
GRID_COLOR = '#30363d'

def style_3d_axis(ax):
    ax.set_facecolor(DARK_BG)
    ax.xaxis.pane.set_facecolor(DARK_BG)
    ax.yaxis.pane.set_facecolor(DARK_BG)
    ax.zaxis.pane.set_facecolor(DARK_BG)
    ax.xaxis.pane.set_edgecolor(GRID_COLOR)
    ax.yaxis.pane.set_edgecolor(GRID_COLOR)
    ax.zaxis.pane.set_edgecolor(GRID_COLOR)
    ax.grid(True, color=GRID_COLOR, linewidth=0.4)

# ----------------------------------------------------------
# 1. Generate a synthetic tumor diagnosis dataset
# ----------------------------------------------------------
n_samples = 400

benign_size = np.random.normal(3.0, 1.0, n_samples // 2)
benign_irregularity = np.random.normal(3.0, 1.0, n_samples // 2)

malignant_size = np.random.normal(7.0, 1.2, n_samples // 2)
malignant_irregularity = np.random.normal(7.0, 1.2, n_samples // 2)

tumor_size = np.concatenate([benign_size, malignant_size])
irregularity = np.concatenate([benign_irregularity, malignant_irregularity])
labels = np.concatenate([np.zeros(n_samples // 2), np.ones(n_samples // 2)])

# Standardize features for numerically stable optimization
size_std = (tumor_size - tumor_size.mean()) / tumor_size.std()
irregularity_std = (irregularity - irregularity.mean()) / irregularity.std()

X = np.column_stack([np.ones(n_samples), size_std, irregularity_std])  # bias + 2 features
y = labels

# ----------------------------------------------------------
# 2. Core logistic regression functions
# ----------------------------------------------------------
def sigmoid(z):
    z = np.clip(z, -500, 500)  # avoid overflow
    return 1.0 / (1.0 + np.exp(-z))

def log_likelihood(beta, X, y):
    z = X @ beta
    p = sigmoid(z)
    eps = 1e-12
    return np.sum(y * np.log(p + eps) + (1 - y) * np.log(1 - p + eps))

def gradient(beta, X, y):
    p = sigmoid(X @ beta)
    return X.T @ (y - p)

def hessian(beta, X):
    p = sigmoid(X @ beta)
    W = p * (1 - p)
    return -(X.T * W) @ X

# ----------------------------------------------------------
# 3. Gradient Ascent optimizer
# ----------------------------------------------------------
def gradient_ascent(X, y, lr=0.01, n_iter=1500):
    beta = np.zeros(X.shape[1])
    history = np.zeros(n_iter)
    for t in range(n_iter):
        grad = gradient(beta, X, y)
        beta = beta + lr * grad
        history[t] = log_likelihood(beta, X, y)
    return beta, history

# ----------------------------------------------------------
# 4. Newton-Raphson (IRLS) optimizer -- fast version
# ----------------------------------------------------------
def newton_raphson(X, y, n_iter=15, tol=1e-8):
    beta = np.zeros(X.shape[1])
    history = []
    for t in range(n_iter):
        grad = gradient(beta, X, y)
        H = hessian(beta, X)
        step = np.linalg.solve(H, grad)  # solve H * step = grad  (faster & more stable than inv(H))
        beta = beta - step
        ll = log_likelihood(beta, X, y)
        history.append(ll)
        if np.linalg.norm(step) < tol:
            break
    return beta, np.array(history)

# ----------------------------------------------------------
# 5. Run both optimizers
# ----------------------------------------------------------
beta_ga, history_ga = gradient_ascent(X, y, lr=0.01, n_iter=1500)
beta_nr, history_nr = newton_raphson(X, y, n_iter=15)

final_ll_ga = history_ga[-1]
final_ll_nr = history_nr[-1]

print("===== Gradient Ascent =====")
print(f"Iterations run     : {len(history_ga)}")
print(f"Final coefficients : {beta_ga}")
print(f"Final log-likelihood: {final_ll_ga:.6f}")
print()
print("===== Newton-Raphson (IRLS) =====")
print(f"Iterations run     : {len(history_nr)}")
print(f"Final coefficients : {beta_nr}")
print(f"Final log-likelihood: {final_ll_nr:.6f}")

# ----------------------------------------------------------
# 6. Plot 1: Convergence comparison (2D)
# ----------------------------------------------------------
fig1, ax1 = plt.subplots(figsize=(9, 6), facecolor=DARK_BG)
ax1.set_facecolor(DARK_BG)
ax1.plot(np.arange(1, len(history_ga) + 1), history_ga,
         color='#58a6ff', linewidth=2, label='Gradient Ascent')
ax1.plot(np.arange(1, len(history_nr) + 1), history_nr,
         color='#f78166', linewidth=2, marker='o', markersize=4, label='Newton-Raphson')
ax1.set_xscale('log')
ax1.set_xlabel('Iteration (log scale)', fontsize=12)
ax1.set_ylabel('Log-Likelihood', fontsize=12)
ax1.set_title('Convergence Speed: Gradient Ascent vs. Newton-Raphson', fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, color=GRID_COLOR, linewidth=0.4)
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 7. Plot 2: 3D log-likelihood surface (slice over 2 of the 3 params)
#    Bias term fixed at its optimal (Newton-Raphson) value
# ----------------------------------------------------------
b1_range = np.linspace(beta_nr[1] - 3, beta_nr[1] + 3, 80)
b2_range = np.linspace(beta_nr[2] - 3, beta_nr[2] + 3, 80)
B1, B2 = np.meshgrid(b1_range, b2_range)

LL_surface = np.zeros_like(B1)
for i in range(B1.shape[0]):
    for j in range(B1.shape[1]):
        beta_temp = np.array([beta_nr[0], B1[i, j], B2[i, j]])
        LL_surface[i, j] = log_likelihood(beta_temp, X, y)

fig2 = plt.figure(figsize=(10, 8), facecolor=DARK_BG)
ax2 = fig2.add_subplot(111, projection='3d')
style_3d_axis(ax2)

surf = ax2.plot_surface(B1, B2, LL_surface, cmap='plasma',
                         linewidth=0, antialiased=True, alpha=0.9)
ax2.set_xlabel('beta_1 (tumor size)', fontsize=10)
ax2.set_ylabel('beta_2 (irregularity)', fontsize=10)
ax2.set_zlabel('Log-Likelihood', fontsize=10)
ax2.set_title('3D Log-Likelihood Landscape', fontsize=14)
ax2.scatter([beta_nr[1]], [beta_nr[2]], [final_ll_nr],
            color='#f78166', s=80, label='MLE optimum')
fig2.colorbar(surf, ax=ax2, shrink=0.5, aspect=10)
ax2.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 8. Plot 3: 3D predicted-probability surface over feature space
# ----------------------------------------------------------
x1_range = np.linspace(size_std.min() - 1, size_std.max() + 1, 60)
x2_range = np.linspace(irregularity_std.min() - 1, irregularity_std.max() + 1, 60)
X1, X2 = np.meshgrid(x1_range, x2_range)

Z_input = np.column_stack([np.ones(X1.size), X1.ravel(), X2.ravel()])
Prob = sigmoid(Z_input @ beta_nr).reshape(X1.shape)

fig3 = plt.figure(figsize=(10, 8), facecolor=DARK_BG)
ax3 = fig3.add_subplot(111, projection='3d')
style_3d_axis(ax3)

ax3.plot_surface(X1, X2, Prob, cmap='viridis', alpha=0.75,
                  linewidth=0, antialiased=True)
ax3.scatter(size_std[y == 0], irregularity_std[y == 0],
            np.zeros(np.sum(y == 0)), color='#58a6ff', s=15, label='Benign (0)')
ax3.scatter(size_std[y == 1], irregularity_std[y == 1],
            np.ones(np.sum(y == 1)), color='#f78166', s=15, label='Malignant (1)')
ax3.set_xlabel('Tumor Size (standardized)', fontsize=10)
ax3.set_ylabel('Irregularity (standardized)', fontsize=10)
ax3.set_zlabel('P(Malignant)', fontsize=10)
ax3.set_title('Fitted Logistic Regression Probability Surface', fontsize=14)
ax3.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 9. Plot 4: 2D decision boundary
# ----------------------------------------------------------
fig4, ax4 = plt.subplots(figsize=(9, 7), facecolor=DARK_BG)
ax4.set_facecolor(DARK_BG)

xx1, xx2 = np.meshgrid(np.linspace(size_std.min() - 1, size_std.max() + 1, 200),
                        np.linspace(irregularity_std.min() - 1, irregularity_std.max() + 1, 200))
grid = np.column_stack([np.ones(xx1.size), xx1.ravel(), xx2.ravel()])
probs = sigmoid(grid @ beta_nr).reshape(xx1.shape)

ax4.contourf(xx1, xx2, probs, levels=25, cmap='coolwarm', alpha=0.6)
ax4.contour(xx1, xx2, probs, levels=[0.5], colors='white', linewidths=2)
ax4.scatter(size_std[y == 0], irregularity_std[y == 0],
            color='#58a6ff', edgecolor='white', s=30, label='Benign (0)')
ax4.scatter(size_std[y == 1], irregularity_std[y == 1],
            color='#f78166', edgecolor='white', s=30, label='Malignant (1)')
ax4.set_xlabel('Tumor Size (standardized)', fontsize=12)
ax4.set_ylabel('Irregularity (standardized)', fontsize=12)
ax4.set_title('Decision Boundary at P = 0.5', fontsize=14)
ax4.legend(fontsize=11)
plt.tight_layout()
plt.show()

Code Walkthrough

Section 1 — Synthetic data generation. We create two clusters of points: “benign” tumors centered around small size/low irregularity, and “malignant” tumors centered around large size/high irregularity, each with Gaussian noise. This mimics a realistic, mildly overlapping medical classification scenario. Features are standardized (zero mean, unit variance) — this is important because it keeps the log-likelihood surface well-conditioned and prevents gradient ascent from oscillating or diverging.

Section 2 — Core math functions. sigmoid() implements $\sigma(z)$ with clipping to avoid floating-point overflow for large $|z|$. log_likelihood() implements $\ell(\boldsymbol{\beta})$ directly from the formula above, with a small epsilon added inside the logarithms to avoid $\log(0)$. gradient() computes $\mathbf{X}^\top(\mathbf{y}-\boldsymbol{\sigma})$, and hessian() computes $-\mathbf{X}^\top \mathbf{W} \mathbf{X}$ using broadcasting (X.T * W) instead of constructing a full diagonal matrix — this avoids an $O(n^2)$ memory allocation and is much faster for larger datasets.

Section 3 — Gradient ascent. This is the “naive” baseline: at every iteration we take a small step in the direction of the gradient. It’s simple but slow — it needs on the order of 1,000+ iterations to approach the optimum, and the step size (lr) has to be tuned carefully; too large and it diverges, too small and convergence takes forever.

Section 4 — Newton-Raphson (IRLS), the fast version. Instead of a fixed-size step, this method rescales the gradient by the inverse curvature (the Hessian), effectively taking a near-optimal step size in every direction automatically. Notice we use np.linalg.solve(H, grad) rather than explicitly computing np.linalg.inv(H) — solving the linear system directly is both faster and numerically more stable than inverting the Hessian. This method typically converges in under 10 iterations, versus 1,500 for gradient ascent — several orders of magnitude fewer computations for essentially the same solution.

Section 5 — Running both optimizers and printing diagnostics. We print the number of iterations, final coefficients, and final log-likelihood for both methods so we can directly compare their efficiency and confirm they converge to (nearly) the same optimum.

Section 6 — Convergence plot. A log-scale x-axis plot showing how quickly the log-likelihood rises for each method — this is where the speed advantage of Newton-Raphson becomes visually obvious.

Section 7 — 3D log-likelihood landscape. We fix the intercept at its optimal value and sweep the two feature coefficients over a grid, computing the log-likelihood at every point. This surface is concave (a single smooth peak), which visually confirms why both optimizers are guaranteed to find the same global maximum — there are no local traps.

Section 8 — 3D probability surface. This shows the fitted sigmoid surface $\sigma(\mathbf{x}^\top\boldsymbol{\beta}^*)$ over the feature space, with the actual data points plotted at $z=0$ or $z=1$ depending on their true label. It’s a great way to see how the S-shaped sigmoid stretches across two dimensions to separate the classes.

Section 9 — 2D decision boundary. A more traditional visualization: the region where the predicted probability crosses 0.5, overlaid with the actual data points, showing the model’s final classification boundary.

===== Gradient Ascent =====
Iterations run     : 1500
Final coefficients : [0.85907155 7.73459307 5.20561006]
Final log-likelihood: -7.611593

===== Newton-Raphson (IRLS) =====
Iterations run     : 12
Final coefficients : [0.90303886 8.06329174 5.3454986 ]
Final log-likelihood: -7.603384

Understanding the Convergence Comparison

This plot is the clearest demonstration of why second-order methods matter. Gradient ascent needs a long, gradual climb — hundreds of tiny steps — to approach the log-likelihood peak, whereas Newton-Raphson essentially “sees” the curvature of the landscape and jumps almost directly to the top within a handful of iterations. For small datasets like ours this speed difference is a curiosity; for large-scale problems with many features, it can be the difference between seconds and hours of training time.

Understanding the Log-Likelihood Landscape

This is the objective function we’re maximizing, rendered as a 3D surface over the two feature coefficients. Its single smooth peak (rather than multiple bumps) reflects the mathematical fact that the logistic regression log-likelihood is concave — there’s exactly one maximum, and both of our optimizers are mathematically guaranteed to find it, just at very different speeds.

Understanding the Fitted Probability Surface

Here we see the S-shaped sigmoid function stretched across two input dimensions. Points sitting near $z=0$ (blue, benign) cluster where the surface is close to 0, and points near $z=1$ (orange, malignant) cluster where the surface approaches 1. The steep “cliff” running diagonally through the middle of the surface is exactly where the model is most uncertain — this cliff, viewed from directly above, is what produces the decision boundary in the next plot.

Understanding the Decision Boundary

The white contour line marks where the model’s predicted probability equals exactly 0.5 — everything on the orange side is classified malignant, everything on the blue side is classified benign. Because our two synthetic clusters are well-separated but not perfectly so, a handful of points naturally fall on the “wrong” side of the boundary, which is realistic and expected in any maximum-likelihood classifier fit to noisy data.

Summary

We formulated logistic regression as a log-likelihood maximization problem, derived its gradient and Hessian, and implemented two optimizers — plain gradient ascent and Newton-Raphson — from scratch using NumPy. Both converge to the same maximum-likelihood solution, but Newton-Raphson does so roughly 100x faster in terms of iteration count, thanks to its use of second-order curvature information. Visualizing the log-likelihood as a 3D surface makes the concavity of the optimization problem tangible, while the fitted probability surface and decision boundary translate the abstract parameter estimates back into an intuitive picture of how the model separates the two classes.

Minimizing the Loss Function in Linear Regression

August 13, 2026

Finding the Optimal Slope and Intercept

Introduction

Linear regression is one of the most fundamental algorithms in machine learning, and at its heart lies a simple but powerful idea: find the line that best fits a set of data points. But what does “best fits” actually mean mathematically? The answer lies in minimizing a loss function, and the most common choice for regression problems is the Mean Squared Error (MSE).

In this article, we’ll build a concrete example from scratch, implement gradient descent in Python, and visualize how the algorithm converges toward the optimal slope and intercept — including a 3D visualization of the loss surface itself.

The Mathematical Formulation

Given a dataset of $n$ points $(x_i, y_i)$, we want to fit a line:

$$
\hat{y}_i = wx_i + b
$$

where $w$ is the slope and $b$ is the intercept. The Mean Squared Error loss function is defined as:

To minimize this loss, we use gradient descent. The partial derivatives of $L$ with respect to $w$ and $b$ are:

$$
\frac{\partial L}{\partial w} = -\frac{2}{n}\sum_{i=1}^{n}x_i(y_i - (wx_i + b))
$$

$$
\frac{\partial L}{\partial b} = -\frac{2}{n}\sum_{i=1}^{n}(y_i - (wx_i + b))
$$

At each iteration, we update the parameters using a learning rate $\eta$:

$$
w \leftarrow w - \eta \frac{\partial L}{\partial w}, \qquad b \leftarrow b - \eta \frac{\partial L}{\partial b}
$$

We repeat this process until the loss converges to a minimum, at which point $w$ and $b$ represent the best-fit line.

The Concrete Example

For this example, we generate synthetic data based on the true relationship $y = 3.5x + 7$ with added Gaussian noise, then use gradient descent to recover the slope (3.5) and intercept (7) purely from the noisy data.

Python Implementation

The code below performs the following steps:

Generates synthetic noisy linear data.
Implements a vectorized (NumPy-based) gradient descent algorithm for speed — avoiding slow Python for loops over individual data points.
Tracks the loss history for convergence analysis.
Computes the loss surface across a grid of $(w, b)$ values for visualization.
Produces four plots: the fitted regression line, the loss convergence curve, a 3D loss surface, and a 2D contour map with the gradient descent path overlaid.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# ---------------------------------------------------------
# 1. Generate synthetic data
# ---------------------------------------------------------
np.random.seed(42)

true_w = 3.5
true_b = 7.0
n_samples = 200

X = np.random.uniform(-10, 10, n_samples)
noise = np.random.normal(0, 4, n_samples)
Y = true_w * X + true_b + noise

# ---------------------------------------------------------
# 2. Vectorized gradient descent
# ---------------------------------------------------------
def compute_loss(w, b, X, Y):
    predictions = w * X + b
    return np.mean((Y - predictions) ** 2)

def gradient_descent(X, Y, w_init=0.0, b_init=0.0,
                      learning_rate=0.01, n_iterations=1000):
    w, b = w_init, b_init
    n = len(X)
    loss_history = []
    w_history = []
    b_history = []

    for i in range(n_iterations):
        predictions = w * X + b
        errors = Y - predictions

        dw = -(2 / n) * np.dot(X, errors)
        db = -(2 / n) * np.sum(errors)

        w -= learning_rate * dw
        b -= learning_rate * db

        loss = compute_loss(w, b, X, Y)
        loss_history.append(loss)
        w_history.append(w)
        b_history.append(b)

    return w, b, loss_history, w_history, b_history

learning_rate = 0.01
n_iterations = 500

final_w, final_b, loss_history, w_history, b_history = gradient_descent(
    X, Y, w_init=0.0, b_init=0.0,
    learning_rate=learning_rate, n_iterations=n_iterations
)

print(f"True parameters:      w = {true_w}, b = {true_b}")
print(f"Estimated parameters: w = {final_w:.4f}, b = {final_b:.4f}")
print(f"Final MSE loss: {loss_history[-1]:.4f}")

# ---------------------------------------------------------
# 3. Compute loss surface for visualization
# ---------------------------------------------------------
w_range = np.linspace(final_w - 5, final_w + 5, 100)
b_range = np.linspace(final_b - 15, final_b + 15, 100)
W_grid, B_grid = np.meshgrid(w_range, b_range)

Loss_grid = np.zeros_like(W_grid)
for i in range(W_grid.shape[0]):
    for j in range(W_grid.shape[1]):
        Loss_grid[i, j] = compute_loss(W_grid[i, j], B_grid[i, j], X, Y)

# ---------------------------------------------------------
# 4. Visualization
# ---------------------------------------------------------
plt.style.use('dark_background')
fig = plt.figure(figsize=(16, 12))

# --- Plot 1: Fitted regression line ---
ax1 = fig.add_subplot(2, 2, 1)
ax1.scatter(X, Y, color='#00d4ff', alpha=0.6, s=25, label='Data points')
x_line = np.linspace(X.min(), X.max(), 100)
y_line = final_w * x_line + final_b
ax1.plot(x_line, y_line, color='#ff6b6b', linewidth=2.5,
          label=f'Fitted line: y = {final_w:.2f}x + {final_b:.2f}')
ax1.set_xlabel('X', fontsize=12)
ax1.set_ylabel('Y', fontsize=12)
ax1.set_title('Linear Regression Fit', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.2)

# --- Plot 2: Loss convergence curve ---
ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(loss_history, color='#ffd93d', linewidth=2)
ax2.set_xlabel('Iteration', fontsize=12)
ax2.set_ylabel('MSE Loss', fontsize=12)
ax2.set_title('Loss Convergence over Iterations', fontsize=14, fontweight='bold')
ax2.grid(alpha=0.2)

# --- Plot 3: 3D loss surface ---
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
surf = ax3.plot_surface(W_grid, B_grid, Loss_grid, cmap='plasma',
                          alpha=0.85, edgecolor='none')
ax3.plot(w_history, b_history, loss_history, color='#00ff88',
          linewidth=2.5, label='Gradient descent path')
ax3.scatter([final_w], [final_b], [loss_history[-1]],
             color='red', s=60, label='Final minimum')
ax3.set_xlabel('Slope (w)', fontsize=10)
ax3.set_ylabel('Intercept (b)', fontsize=10)
ax3.set_zlabel('MSE Loss', fontsize=10)
ax3.set_title('3D Loss Surface with Descent Path', fontsize=14, fontweight='bold')
ax3.legend(fontsize=9)

# --- Plot 4: Contour map with descent path ---
ax4 = fig.add_subplot(2, 2, 4)
contour = ax4.contourf(W_grid, B_grid, Loss_grid, levels=40, cmap='plasma')
ax4.plot(w_history, b_history, color='#00ff88', linewidth=2,
          marker='o', markersize=2, label='Gradient descent path')
ax4.scatter([final_w], [final_b], color='red', s=80,
             marker='*', label='Final minimum', zorder=5)
ax4.set_xlabel('Slope (w)', fontsize=12)
ax4.set_ylabel('Intercept (b)', fontsize=12)
ax4.set_title('Loss Contour Map (Top View)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
fig.colorbar(contour, ax=ax4, label='MSE Loss')

plt.tight_layout()
plt.show()

True parameters:      w = 3.5, b = 7.0
Estimated parameters: w = 3.4844, b = 7.2644
Final MSE loss: 14.9412

Code Walkthrough

Data generation: We create 200 points along the line $y = 3.5x + 7$, then inject Gaussian noise with a standard deviation of 4. This simulates real-world measurement noise, giving gradient descent a genuine estimation problem rather than a trivial exact fit.

Vectorized gradient computation: Instead of looping over each data point with a Python for loop (which would be extremely slow for large datasets), the gradients are computed using np.dot(X, errors) and np.sum(errors). This leverages NumPy’s underlying C implementation, making the computation orders of magnitude faster than a pure Python loop — critical when scaling to larger datasets or more iterations.

Gradient descent loop: At each of the 500 iterations, predictions are computed for the entire dataset at once, the error vector is calculated, and both gradients ($\partial L/\partial w$ and $\partial L/\partial b$) are computed simultaneously. The parameters are then nudged in the direction that reduces the loss, scaled by the learning rate of 0.01.

Loss surface computation: To visualize the shape of the loss function itself, we build a grid of candidate $(w, b)$ pairs surrounding the final solution and compute the MSE loss at every grid point. This produces a bowl-shaped surface — a hallmark of MSE loss for linear regression, since it’s a convex quadratic function with a single global minimum.

Interpreting the Results

The top-left plot shows the raw noisy data alongside the line found by gradient descent — despite the noise, the algorithm recovers a slope and intercept very close to the true values of $w=3.5$ and $b=7$.

The top-right plot shows the loss dropping sharply in the first several iterations before flattening out, which is typical of gradient descent: large early steps followed by fine-tuning as the algorithm approaches the minimum.

The 3D surface plot is the most illuminating: it reveals the loss function as a smooth, convex bowl in $(w, b)$ space. The green trajectory traces the exact path taken by gradient descent, starting from $w=0, b=0$ and spiraling down toward the bottom of the bowl — the point where MSE is minimized.

The contour map gives a bird’s-eye view of the same bowl, making it easy to see how the descent path curves toward the minimum, following the steepest downhill direction at every step, which is exactly what the negative gradient represents.

Why This Matters

This simple example illustrates the core mechanism behind training almost every regression-based machine learning model, from simple linear regression to the first layer of a neural network. Understanding how MSE creates a convex loss landscape — and how gradient descent navigates it — provides the foundation for understanding more complex loss surfaces in deep learning, where the landscape is no longer a simple bowl but is optimized using the very same underlying principle: follow the gradient downhill.

Solving the 2D Principle of Least Action

August 12, 2026

A Discretized Path Optimization Example

What is the principle of least action?

In classical mechanics, a particle does not simply “obey forces” — it follows the trajectory that minimizes (or more precisely, makes stationary) a quantity called the action $S$, defined as the time-integral of the Lagrangian $L = T - U$ (kinetic energy minus potential energy):

$$
S[x(t), y(t)] = \int_0^T L , dt = \int_0^T \left[ \frac{1}{2}m\left(\dot{x}^2 + \dot{y}^2\right) - mgy \right] dt
$$

Applying the Euler–Lagrange equation to this functional recovers Newton’s familiar equations of motion:

$$
\frac{d}{dt}\frac{\partial L}{\partial \dot{x}} - \frac{\partial L}{\partial x} = 0 \quad\Rightarrow\quad \ddot{x} = 0
$$

$$
\frac{d}{dt}\frac{\partial L}{\partial \dot{y}} - \frac{\partial L}{\partial y} = 0 \quad\Rightarrow\quad \ddot{y} = -g
$$

This is exactly projectile motion. But instead of solving the differential equation directly, we can find the same trajectory numerically by discretizing the path into a finite number of points and directly minimizing the action with an optimizer. This is a nice illustration of the variational nature of mechanics: instead of “starting somewhere with some velocity,” we fix the start point and the end point and let the optimizer find the path connecting them that minimizes $S$.

Example problem

A particle of mass $m = 1,\text{kg}$ starts at $(x_0, y_0) = (0, 0)$ and must arrive at $(x_N, y_N) = (10, 0)$ after $T = 1.5$ seconds, under uniform gravity $g = 9.8,\text{m/s}^2$. The path is discretized into $N = 40$ time steps (41 points total). The two endpoints are fixed; the 39 interior points are free variables that the optimizer adjusts to minimize the discretized action. We then compare the result against the exact analytical parabola.

Source code

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)

# ---------------------------------------------------------
# 1. Physical parameters
# ---------------------------------------------------------
m = 1.0          # mass [kg]
g = 9.8          # gravitational acceleration [m/s^2]

# ---------------------------------------------------------
# 2. Boundary conditions (fixed start and end points)
# ---------------------------------------------------------
x0, y0 = 0.0, 0.0     # position at t = 0
xN, yN = 10.0, 0.0    # position at t = T
T = 1.5                # total flight time [s]
N = 40                  # number of time steps (path has N+1 points)
dt = T / N
t = np.linspace(0, T, N + 1)

# ---------------------------------------------------------
# 3. Analytical solution (exact projectile parabola)
# ---------------------------------------------------------
vx0 = (xN - x0) / T
vy0 = (yN - y0 + 0.5 * g * T**2) / T
x_analytic = x0 + vx0 * t
y_analytic = y0 + vy0 * t - 0.5 * g * t**2

# ---------------------------------------------------------
# 4. Discretized action S[x, y]  (vectorized -> fast)
# ---------------------------------------------------------
def action(z):
    x = np.empty(N + 1)
    y = np.empty(N + 1)
    x[0], x[-1] = x0, xN
    y[0], y[-1] = y0, yN
    x[1:-1] = z[:N - 1]
    y[1:-1] = z[N - 1:]

    vx = np.diff(x) / dt
    vy = np.diff(y) / dt
    KE = 0.5 * m * (vx**2 + vy**2)
    y_mid = 0.5 * (y[:-1] + y[1:])   # midpoint rule for potential energy
    PE = m * g * y_mid

    return np.sum(KE - PE) * dt

# ---------------------------------------------------------
# 5. Initial guess: straight line between the two endpoints
# ---------------------------------------------------------
x_init = np.linspace(x0, xN, N + 1)
y_init = np.linspace(y0, yN, N + 1)
z0 = np.concatenate([x_init[1:-1], y_init[1:-1]])

# ---------------------------------------------------------
# 6. Minimize the action, recording convergence history
# ---------------------------------------------------------
history = []
def callback(z):
    history.append(action(z))

result = minimize(action, z0, method='L-BFGS-B', callback=callback,
                   options={'maxiter': 500, 'ftol': 1e-14, 'gtol': 1e-12})

z_opt = result.x
x_opt = np.concatenate(([x0], z_opt[:N - 1], [xN]))
y_opt = np.concatenate(([y0], z_opt[N - 1:], [yN]))

print("Optimization success :", result.success)
print("Minimum action S_min :", result.fun)
print("Max error vs analytic (x):", np.max(np.abs(x_opt - x_analytic)))
print("Max error vs analytic (y):", np.max(np.abs(y_opt - y_analytic)))

# ---------------------------------------------------------
# 7. Plot 1: trajectory comparison (numeric vs analytic)
# ---------------------------------------------------------
plt.figure(figsize=(9, 5))
plt.plot(x_analytic, y_analytic, 'b-', linewidth=2, label='Analytical trajectory')
plt.plot(x_opt, y_opt, 'ro', markersize=4, label='Action-minimized (numerical)')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title('Trajectory that minimizes the action')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 8. Plot 2: convergence of the action during optimization
# ---------------------------------------------------------
plt.figure(figsize=(9, 5))
plt.plot(history, 'g-o', markersize=3)
plt.xlabel('Iteration')
plt.ylabel('Action S')
plt.title('Convergence of the action toward its minimum')
plt.grid(True)
plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 9. Plot 3: 3D action landscape around the optimal path
# ---------------------------------------------------------
interior_len = N - 1
i1 = interior_len // 3
i2 = 2 * interior_len // 3

y1_center = z_opt[N - 1 + i1]
y2_center = z_opt[N - 1 + i2]

span = 3.0
y1_range = np.linspace(y1_center - span, y1_center + span, 40)
y2_range = np.linspace(y2_center - span, y2_center + span, 40)
Y1, Y2 = np.meshgrid(y1_range, y2_range)
S_grid = np.empty_like(Y1)

for a in range(Y1.shape[0]):
    for b in range(Y1.shape[1]):
        z_temp = z_opt.copy()
        z_temp[N - 1 + i1] = Y1[a, b]
        z_temp[N - 1 + i2] = Y2[a, b]
        S_grid[a, b] = action(z_temp)

fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(Y1, Y2, S_grid, cmap='viridis', alpha=0.9, edgecolor='none')
ax.scatter([y1_center], [y2_center], [action(z_opt)],
           color='red', s=70, label='Physical path (minimum)')
ax.set_xlabel(f'y at interior point {i1}')
ax.set_ylabel(f'y at interior point {i2}')
ax.set_zlabel('Action S')
ax.set_title('Action landscape: the classical path sits at the bottom of the bowl')
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.legend()
plt.tight_layout()
plt.show()

Execution Results

Running the code above in Google Colaboratory produces the following console output:

Optimization success : True
Minimum action S_min : 19.836149348979536
Max error vs analytic (x): 2.728448545319395e-06
Max error vs analytic (y): 1.8400361176951208e-06

Code walkthrough

Parameters and boundary conditions. The particle’s mass, gravity, and the two fixed endpoints $(x_0,y_0)$ and $(x_N,y_N)$ are set first, along with the total flight time $T$ and the number of discretization steps $N$. Note that this is a boundary value problem — we specify position at both ends of time, not an initial velocity. This is precisely the kind of problem the action principle handles naturally, whereas Newton’s equations alone would require an extra step to find the correct launch velocity.

Analytical solution. Because the exact equations of motion are linear ($\ddot x = 0$, $\ddot y=-g$), we can solve for the initial velocity that connects the two endpoints in closed form and generate the exact reference parabola. This serves as ground truth to validate the numerical optimization.

The discretized action. The action() function is the heart of the simulation. The continuous integral is replaced by a Riemann sum over $N$ segments. Velocities are approximated with forward differences, $\dot{x}i \approx (x{i+1}-x_i)/dt$, and the potential energy term uses the midpoint rule $y_{\text{mid}} = (y_i+y_{i+1})/2$, which keeps the discretization consistent to second order. Everything is written with NumPy array operations (np.diff, vectorized arithmetic) rather than Python for loops, so evaluating the action for a candidate path is essentially instantaneous — this matters because the optimizer will call this function hundreds of times.

Optimization. The interior points (everything except the two fixed endpoints) are flattened into a single vector z and handed to scipy.optimize.minimize using the L-BFGS-B method, a quasi-Newton algorithm well suited to smooth, unconstrained problems. Since the action here is a convex quadratic function of the interior coordinates (kinetic term is quadratic and positive-definite, potential term is linear), the minimization problem has a single global minimum and the optimizer converges quickly and reliably — this is why no execution errors or convergence failures occur. A callback records the action value at every iteration so we can visualize convergence afterward.

Validation. After optimization, the numerical path is compared point-by-point against the analytical parabola. Because the discretized Euler–Lagrange conditions for this particular Lagrangian reduce to the same recursion as the exact free-fall solution, the numerical and analytical trajectories should match to within numerical tolerance (typically well under $10^{-6}$ m).

Results

1. Trajectory comparison

This plot overlays the exact analytical parabola (blue line) with the individual points found by minimizing the discretized action (red dots). If the optimization is correct, the red dots should sit exactly on the blue curve — a direct visual confirmation that “minimizing the action” and “solving Newton’s equations” produce the same physical trajectory.

2. Convergence of the action

This plot shows the value of the action $S$ at each iteration of the optimizer, starting from the straight-line initial guess and decreasing monotonically toward the true minimum. Because the problem is a convex quadratic, the curve typically drops sharply within the first few iterations and then flattens out as it reaches the optimum — a hallmark of a well-posed, well-conditioned optimization.

3. The 3D action landscape

This is the most intuitive visualization of the principle of least action. Two of the interior points’ height coordinates are varied over a grid while all other coordinates are held fixed at their optimized values, and the resulting action is plotted as a 3D surface. The surface forms a smooth bowl shape, and the red marker — the path found by the optimizer — sits precisely at the bottom. This is a direct geometric demonstration of what “least action” means: among all nearby possible paths, the physical trajectory is the one occupying the lowest point of the action landscape.