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.

Solving Steady-State Heat Conduction with Energy Minimization

August 11, 2026

Many physical systems reach a steady (equilibrium) state not because something is being “solved” in the traditional sense, but because nature is minimizing a quantity — usually an energy functional. Steady-state heat conduction is a textbook example: instead of solving the heat equation directly, we can find the temperature distribution that minimizes a scalar energy functional. This is the Dirichlet principle, and it turns out the exact same mathematical machinery describes irrotational, incompressible fluid flow (potential flow), since both are governed by the Laplace/Poisson equation.

In this article, we’ll take a concrete example — a 2D metal plate heated from the left and containing an internal heat source — and find its steady-state temperature field by minimizing energy directly, rather than solving a linear PDE the “usual” way.

The Physics: From Energy to Equation

For steady-state heat conduction with conductivity $k$ and a heat source density $q(x,y)$, the temperature field $T(x,y)$ that the system settles into is the one that minimizes the functional:

$$
E[T] = \int_\Omega \left[ \frac{k}{2} , |\nabla T|^2 ;-; q(x,y), T \right] , dx, dy
$$

The first term is analogous to elastic/kinetic energy stored in the temperature gradient; the second term is the work done by the heat source. Taking the variational derivative and setting it to zero (Euler–Lagrange equation) gives:

$$
\frac{\delta E}{\delta T} = -k \nabla^2 T - q = 0 \quad \Longrightarrow \quad k,\nabla^2 T + q = 0
$$

This is exactly the steady-state heat equation. So minimizing $E[T]$ and solving $k\nabla^2 T + q = 0$ are mathematically equivalent — this is the essence of the Dirichlet principle.

The same functional form, with $T$ replaced by a velocity potential $\phi$ and $q=0$, is the kinetic energy of an incompressible, irrotational fluid:

$$
E[\phi] = \int_\Omega \frac{1}{2},|\nabla \phi|^2 , dx, dy, \qquad \nabla^2 \phi = 0
$$

which is why the code below applies equally well to potential-flow problems.

Discretizing the Problem

On a finite-difference grid, the continuous functional becomes a quadratic form in the vector of unknown (interior) temperatures $\mathbf{T}$:

$$
E(\mathbf{T}) = \frac{1}{2},\mathbf{T}^\top A, \mathbf{T} - \mathbf{f}^\top \mathbf{T}
$$

where $A$ is the discrete (negative) Laplacian operator (positive-definite, sparse) and $\mathbf{f}$ bundles the heat source and boundary-condition contributions. Its gradient is:

$$
\nabla E(\mathbf{T}) = A\mathbf{T} - \mathbf{f}
$$

Setting this to zero recovers the linear system $A\mathbf{T} = \mathbf{f}$ — but instead of solving it as a linear system, we’ll hand $E(\mathbf{T})$ and $\nabla E(\mathbf{T})$ to an optimizer and let it find the minimum directly, which is the whole point of this exercise.

The Example Problem

A square plate, $1\text{m} \times 1\text{m}$:

Left wall ($x=0$) held at $T = 100°C$ (a hot wall)
Top, bottom, and right walls held at $T = 0°C$
A localized internal heat source (a Gaussian “heater”) embedded inside the plate near $(0.65, 0.5)$

We’ll compute the steady-state temperature field two different ways — a general-purpose optimizer (L-BFGS-B) and a specialized, much faster energy-minimizing algorithm (Conjugate Gradient) — and cross-check both against a direct sparse linear solve.

Python Source Code

# ==========================================================
# Steady-State Heat Conduction via Energy Minimization
# ==========================================================
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time

# ---------- 1. Grid and physical parameters ----------
Nx, Ny = 51, 51
Lx, Ly = 1.0, 1.0
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y, indexing='ij')

k_cond = 1.0      # thermal conductivity
T_hot  = 100.0    # left wall temperature
T_cold = 0.0      # other walls temperature

# internal heat source (Gaussian "heater")
q0, xs, ys, sigma = 3000.0, 0.65, 0.5, 0.06
Q = q0 * np.exp(-((X - xs)**2 + (Y - ys)**2) / (2 * sigma**2))

# ---------- 2. Build the discrete energy operator A and load vector f ----------
mx, my = Nx - 2, Ny - 2   # number of interior nodes in each direction

main_x = -2.0 * np.ones(mx)
off_x  = np.ones(mx - 1)
D2x = sp.diags([off_x, main_x, off_x], [-1, 0, 1], format='csr') / dx**2

main_y = -2.0 * np.ones(my)
off_y  = np.ones(my - 1)
D2y = sp.diags([off_y, main_y, off_y], [-1, 0, 1], format='csr') / dy**2

Ix = sp.identity(mx, format='csr')
Iy = sp.identity(my, format='csr')

L = sp.kron(D2x, Iy, format='csr') + sp.kron(Ix, D2y, format='csr')
A = (-L).tocsr()  # positive-definite discrete energy operator

Q_int = Q[1:-1, 1:-1]
boundary_term = np.zeros((mx, my))
boundary_term[0, :] += T_hot / dx**2   # contribution from the hot left wall

f = (Q_int / k_cond + boundary_term).flatten()

# ---------- 3. Discrete energy functional and its gradient ----------
def energy(t):
    return 0.5 * t @ (A @ t) - f @ t

def grad_energy(t):
    return A @ t - f

# ---------- 4. Reference solution (direct sparse solve) ----------
t0 = time.time()
T_ref = spla.spsolve(A.tocsc(), f)
t_direct = time.time() - t0

# ---------- 5. Energy minimization: general-purpose optimizer (L-BFGS-B) ----------
energy_history_bfgs = []
def callback_bfgs(t):
    energy_history_bfgs.append(energy(t))

t0 = time.time()
res_bfgs = minimize(energy, np.zeros(mx * my), jac=grad_energy,
                     method='L-BFGS-B', callback=callback_bfgs,
                     options={'maxiter': 800, 'ftol': 1e-14, 'gtol': 1e-8})
t_bfgs = time.time() - t0
T_bfgs = res_bfgs.x

# ---------- 6. Energy minimization: fast specialized version (Conjugate Gradient) ----------
energy_history_cg = []
def callback_cg(t):
    energy_history_cg.append(energy(t))

t0 = time.time()
T_cg, info = spla.cg(A, f, rtol=1e-12, maxiter=1000, callback=callback_cg)
t_cg = time.time() - t0

print(f"Direct sparse solve : {t_direct*1000:8.2f} ms")
print(f"L-BFGS-B minimizer  : {t_bfgs*1000:8.2f} ms  ({len(energy_history_bfgs)} iterations)")
print(f"Conjugate Gradient  : {t_cg*1000:8.2f} ms  ({len(energy_history_cg)} iterations)")
print(f"max|T_bfgs - T_ref| = {np.max(np.abs(T_bfgs - T_ref)):.3e}")
print(f"max|T_cg   - T_ref| = {np.max(np.abs(T_cg   - T_ref)):.3e}")

# ---------- 7. Assemble the full temperature field (with boundary values) ----------
def assemble_full(t_int):
    T_full = np.zeros((Nx, Ny))
    T_full[-1, :] = T_cold
    T_full[:, 0]  = T_cold
    T_full[:, -1] = T_cold
    T_full[0, :]  = T_hot
    T_full[1:-1, 1:-1] = t_int.reshape(mx, my)
    return T_full

T_field = assemble_full(T_cg)

# ---------- 8. Heat flux field  q = -k*grad(T)  (a fluid-like transport field) ----------
dTdx, dTdy = np.gradient(T_field, dx, dy)
qx, qy = -k_cond * dTdx, -k_cond * dTdy

# ---------- 9. Visualization ----------
fig1, ax1 = plt.subplots(figsize=(7, 6))
c = ax1.contourf(X, Y, T_field, levels=40, cmap='inferno')
skip = 4
ax1.quiver(X[::skip, ::skip], Y[::skip, ::skip],
           qx[::skip, ::skip], qy[::skip, ::skip],
           color='white', scale=4000, width=0.003)
ax1.set_xlabel('x [m]')
ax1.set_ylabel('y [m]')
ax1.set_title('Steady-State Temperature Field and Heat Flux')
fig1.colorbar(c, ax=ax1, label='Temperature [°C]')
plt.tight_layout()
plt.show()

fig2, ax2 = plt.subplots(figsize=(7, 5))
ax2.plot(energy_history_bfgs, 'o-', label='L-BFGS-B', markersize=3)
ax2.plot(energy_history_cg, 's-', label='Conjugate Gradient', markersize=3)
ax2.set_xlabel('Iteration')
ax2.set_ylabel('Energy E[T]')
ax2.set_title('Convergence of the Energy Functional')
ax2.legend()
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()

fig3 = plt.figure(figsize=(8, 6))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(X, Y, T_field, cmap='inferno', linewidth=0, antialiased=True)
ax3.set_xlabel('x [m]')
ax3.set_ylabel('y [m]')
ax3.set_zlabel('T [°C]')
ax3.set_title('3D Temperature Distribution')
fig3.colorbar(surf, ax=ax3, shrink=0.6, label='Temperature [°C]')
plt.tight_layout()
plt.show()

Direct sparse solve :    36.16 ms
L-BFGS-B minimizer  :  1754.44 ms  (222 iterations)
Conjugate Gradient  :    86.39 ms  (169 iterations)
max|T_bfgs - T_ref| = 8.850e-05
max|T_cg   - T_ref| = 8.624e-11

Code Walkthrough

Section 1 — Grid setup. We define a $51 \times 51$ grid over a unit square. X, Y = np.meshgrid(..., indexing='ij') keeps the array’s first axis aligned with $x$ and the second with $y$, which matters later when we flatten arrays into vectors.

Section 2 — Building $A$ and $\mathbf{f}$. This is the heart of the discretization. D2x and D2y are standard tridiagonal second-derivative operators (the 1D Laplacian stencil $[1, -2, 1]/h^2$) built only over the interior grid points — boundary points are not unknowns, they’re known data. We combine the two 1D operators into a 2D Laplacian using a Kronecker sum: sp.kron(D2x, Iy) + sp.kron(Ix, D2y). This is the standard trick for turning a 2D finite-difference stencil into a single sparse matrix without writing nested loops. We negate the Laplacian (A = -L) so that $A$ is positive-definite, matching the sign convention required for an energy-minimization (rather than energy-maximization) problem.

Because the hot boundary ($T=100°C$ at $x=0$) is not one of our unknowns, its influence has to be folded into the right-hand side vector $\mathbf{f}$ manually — this is the boundary_term array, which only affects the row of interior nodes immediately adjacent to the hot wall.

Section 3 — Energy and gradient functions. energy(t) directly implements $E(\mathbf{T}) = \frac12 \mathbf{T}^\top A \mathbf{T} - \mathbf{f}^\top \mathbf{T}$ and grad_energy(t) implements its exact analytic gradient $A\mathbf{T}-\mathbf{f}$. Supplying the exact gradient (rather than letting the optimizer estimate it by finite differences) is the single biggest performance factor here — without it, scipy.optimize.minimize would need $O(n)$ extra function evaluations per iteration just to approximate the gradient.

Section 4 — Reference solution. spla.spsolve solves $A\mathbf{T}=\mathbf{f}$ directly using sparse LU factorization. This is our ground truth to validate the optimization-based approaches against.

Sections 5 & 6 — Two flavors of energy minimization.

L-BFGS-B is a general-purpose quasi-Newton optimizer. It works for any differentiable energy functional (including nonlinear ones), which is why it’s the natural first choice pedagogically.
Conjugate Gradient (scipy.sparse.linalg.cg) is not just “a solver” — it is literally an algorithm that minimizes a quadratic energy functional of exactly this form, one conjugate direction at a time. Because it’s purpose-built for symmetric positive-definite quadratic problems (like ours) rather than general nonlinear ones, it converges in dramatically fewer, cheaper iterations. This is our “fast version”: same energy-minimization idea, specialized to the structure of the problem for a large speedup, while still being philosophically the same energy-minimization approach as L-BFGS-B.

Section 7 — Reassembling the full grid. The optimizer only ever sees the vector of interior unknowns; assemble_full puts the known boundary values back around it to produce a complete $(N_x, N_y)$ temperature field for plotting.

Section 8 — Heat flux. Fourier’s law says heat flows down the temperature gradient: $\mathbf{q} = -k\nabla T$. This vector field is mathematically identical in form to a fluid velocity field in potential flow, which is why we visualize it with arrows — it makes the heat-conduction/fluid-flow analogy visually concrete.

Why Conjugate Gradient Wins on Speed

L-BFGS-B has to build up an approximation to the inverse Hessian using gradient history, which costs extra bookkeeping every iteration and typically needs many iterations to reach high accuracy on an ill-conditioned quadratic like a discretized Laplacian. Conjugate Gradient, by contrast, exploits the fact that our energy is exactly quadratic with a known, sparse, symmetric positive-definite matrix — each iteration is just one sparse matrix-vector product plus a few dot products, and CG is mathematically guaranteed to reach the exact minimum within $n$ steps (and in practice, far fewer, since the Poisson matrix’s eigenvalues cluster). The console output above should show CG finishing in a small fraction of the time L-BFGS-B needs, while both agree with the direct solve to within numerical precision.

Understanding the Results

Figure 1 — Temperature field and heat flux. The heatmap shows temperature decaying smoothly from the hot left wall (bright) toward the cooler edges, with a distinct hot spot around the embedded heat source. The white arrows show the heat flux vector field $\mathbf{q}=-k\nabla T$ — heat flows from hot to cold, “downhill” on the temperature surface, exactly like a fluid flowing from high to low potential.

Figure 2 — Energy convergence. This plot shows the value of the energy functional $E(\mathbf{T})$ at each iteration for both algorithms. Both curves should descend monotonically toward the same minimum value (since $A$ is positive-definite, the quadratic has a unique global minimum), but the Conjugate Gradient curve should reach convergence in far fewer iterations — visually demonstrating that both algorithms are doing the same conceptual thing (rolling downhill on the energy landscape) at very different speeds.

Figure 3 — 3D temperature surface. This is the same data as Figure 1, but viewed as a literal energy landscape — high near the hot wall and the internal heater, sloping down toward the cold boundaries. Seeing it as a physical surface makes the “minimization” framing intuitive: the true physical steady state is the shape a stretched elastic membrane would settle into if pinned at the given boundary heights and pushed up by the heat source, minimizing its stored elastic energy.

Closing Thoughts

Framing steady-state heat conduction as an energy-minimization problem does more than provide an alternative numerical method — it reveals the same mathematics underlying incompressible potential flow, electrostatics, and membrane mechanics. Once a problem is expressed as minimizing $\frac12\mathbf{T}^\top A\mathbf{T} - \mathbf{f}^\top\mathbf{T}$, any tool from the optimization world — gradient descent, L-BFGS, conjugate gradient, or beyond — becomes a legitimate physics solver, and the choice between them becomes purely a question of computational efficiency rather than modeling correctness.

Minimizing the Weight of a Truss Structure with Python

August 10, 2026

A Two-Bar Truss Case Study

Structural weight minimization is one of the classic entry points into engineering optimization. The idea is simple to state but rich in behavior: given a structure that has to carry a load safely, find the combination of member cross-sectional areas and member lengths (or, equivalently, the geometry that determines those lengths) that uses the least material while still satisfying stress and stiffness requirements.

In this article we work through a concrete, fully worked example — a symmetric two-bar truss — set up the mechanics equations, formulate it as a constrained nonlinear optimization problem in Python, solve it with scipy.optimize, and then visualize the result from several angles, including a 3D view of the design space.

The Structure

Picture two support points fixed on the ground, separated by a fixed half-width $b$, and a single apex node above them where a vertical load $P$ is applied. Two bars connect the apex to each support, forming a symmetric “A-frame” truss. Both bars share the same cross-sectional area $A$.

The two design variables are:

$A$ — the cross-sectional area of each bar
$h$ — the height of the apex above the base

Because the base width $b$ is fixed, the bar length is fully determined by $h$:

$$
L(h) = \sqrt{b^2 + h^2}
$$

So treating $h$ as a design variable is equivalent to treating the member length directly as a design variable — raising the apex lengthens the bars and simultaneously changes the load angle.

The Mechanics

By vertical equilibrium at the apex, each bar carries an axial force:

$$
F(h) = \frac{P , L(h)}{2h}
$$

The axial stress in each bar is:

$$
\sigma(A,h) = \frac{F(h)}{A} = \frac{P , L(h)}{2Ah}
$$

The vertical deflection of the apex node under the load (from the bar’s elastic elongation, projected back through the load angle) is:

$$
\delta(A,h) = \frac{P , L(h)^3}{2AEh^2}
$$

where $E$ is Young’s modulus of the bar material.

The total structural weight is:

$$
W(A,h) = 2 \rho A L(h) = 2 \rho A \sqrt{b^2+h^2}
$$

where $\rho$ is the material density.

The Optimization Problem

We want to minimize weight subject to an allowable stress $\sigma_{allow}$ and an allowable deflection $\delta_{allow}$:

$$
\min_{A,,h} \quad W(A,h) = 2\rho A \sqrt{b^2+h^2}
$$

$$
\text{subject to} \quad \sigma(A,h) \le \sigma_{allow}, \qquad \delta(A,h) \le \delta_{allow}
$$

$$
A_{min} \le A \le A_{max}, \qquad h_{min} \le h \le h_{max}
$$

This is a genuinely interesting problem because increasing $h$ makes the load angle steeper (reducing the axial force, and therefore the area needed to satisfy the stress limit), but it also increases $L$, which pushes the deflection up as $L^3$. There is no free lunch: some intermediate height minimizes the total weight, and the optimizer has to find it.

Python Implementation

The code below is a single, self-contained cell. It sets up the mechanics, runs a gradient-based constrained optimizer (SLSQP), cross-checks the numerical result against a fast analytical reduction of the problem, and produces four plots.

# ==========================================================
# Weight Minimization of a Two-Bar Truss
# Design Variables: Cross-sectional Area (A) and Height (h)
# ==========================================================

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

# ----------------------------------------------------------
# 1. Fixed Physical Parameters
# ----------------------------------------------------------
P           = 50000.0     # Applied vertical load at apex [N]
b           = 1.0         # Fixed half base width [m]
rho         = 7850.0      # Density of steel [kg/m^3]
E           = 200e9       # Young's modulus of steel [Pa]
sigma_allow = 165e6       # Allowable axial stress [Pa]
delta_allow = 5e-3        # Allowable vertical deflection at apex [m]

# ----------------------------------------------------------
# 2. Geometry / Mechanics helper functions
# ----------------------------------------------------------
def bar_length(h):
    """Length of one bar as a function of apex height h."""
    return np.sqrt(b**2 + h**2)

def axial_force(h):
    """Axial force carried by one bar (symmetric 2-bar truss)."""
    L = bar_length(h)
    return P * L / (2.0 * h)

def axial_stress(A, h):
    """Axial stress in one bar."""
    return axial_force(h) / A

def apex_deflection(A, h):
    """Vertical deflection at the apex node."""
    L = bar_length(h)
    return (P * L**3) / (2.0 * A * E * h**2)

def total_weight(x):
    """Total weight of the two bars: x = [A, h]."""
    A, h = x
    L = bar_length(h)
    return 2.0 * rho * A * L

# ----------------------------------------------------------
# 3. Constraint functions (SLSQP requires g(x) >= 0 form)
# ----------------------------------------------------------
def stress_constraint(x):
    A, h = x
    return sigma_allow - axial_stress(A, h)

def deflection_constraint(x):
    A, h = x
    return delta_allow - apex_deflection(A, h)

# ----------------------------------------------------------
# 4. Optimization set-up
# ----------------------------------------------------------
A_min, A_max = 1.0e-5, 5.0e-3     # [m^2]  (0.1 cm^2 -- 50 cm^2)
h_min, h_max = 0.5, 4.0           # [m]

bounds = [(A_min, A_max), (h_min, h_max)]
x0     = np.array([5.0e-4, 1.5])   # initial guess [A, h]

history = [x0.copy()]
def record(xk):
    history.append(xk.copy())

result = minimize(
    total_weight,
    x0,
    method='SLSQP',
    bounds=bounds,
    constraints=[
        {'type': 'ineq', 'fun': stress_constraint},
        {'type': 'ineq', 'fun': deflection_constraint},
    ],
    callback=record,
    options={'maxiter': 200, 'ftol': 1e-10}
)

A_opt, h_opt = result.x
W_opt        = result.fun
history      = np.array(history)

# ----------------------------------------------------------
# 5. Report
# ----------------------------------------------------------
print("=" * 55)
print(" OPTIMIZATION RESULT ")
print("=" * 55)
print(f"Success              : {result.success}")
print(f"Iterations           : {result.nit}")
print(f"Optimal area  A*     : {A_opt*1e4:.4f} cm^2")
print(f"Optimal height h*    : {h_opt:.4f} m")
print(f"Optimal bar length L*: {bar_length(h_opt):.4f} m")
print(f"Minimum weight W*    : {W_opt:.4f} kg")
print(f"Stress at optimum    : {axial_stress(A_opt,h_opt)/1e6:.3f} MPa "
      f"(allowable {sigma_allow/1e6:.1f} MPa)")
print(f"Deflection at optimum: {apex_deflection(A_opt,h_opt)*1e3:.4f} mm "
      f"(allowable {delta_allow*1e3:.1f} mm)")
print("=" * 55)

# ----------------------------------------------------------
# 6. Analytical cross-check: reduce to a 1-D problem in h
#    (since weight increases monotonically with A, the
#     optimal A for a given h is exactly the tighter of the
#     two constraints)
# ----------------------------------------------------------
h_fine = np.linspace(h_min, h_max, 2000)
L_fine = bar_length(h_fine)

A_stress_curve = P * L_fine / (2.0 * h_fine * sigma_allow)
A_defl_curve   = (P * L_fine**3) / (2.0 * E * h_fine**2 * delta_allow)
A_active_curve = np.maximum(A_stress_curve, A_defl_curve)
W_active_curve = 2.0 * rho * A_active_curve * L_fine

idx_min   = np.argmin(W_active_curve)
h_star_1d = h_fine[idx_min]
W_star_1d = W_active_curve[idx_min]

print(f"1-D cross-check   h* : {h_star_1d:.4f} m,  W*: {W_star_1d:.4f} kg")

# ----------------------------------------------------------
# 7. Figure 1: 3-D weight surface with optimization path
# ----------------------------------------------------------
A_grid = np.linspace(A_min, A_max, 120)
h_grid = np.linspace(h_min, h_max, 120)
AG, HG = np.meshgrid(A_grid, h_grid)
WG = 2.0 * rho * AG * bar_length(HG)

fig1 = plt.figure(figsize=(10, 7))
ax1 = fig1.add_subplot(111, projection='3d')
surf = ax1.plot_surface(AG*1e4, HG, WG, cmap='viridis',
                         alpha=0.75, linewidth=0, antialiased=True)

path_W = np.array([total_weight(pt) for pt in history])
ax1.plot(history[:, 0]*1e4, history[:, 1], path_W,
         color='red', marker='o', markersize=4, linewidth=2,
         label='Optimization path')
ax1.scatter([A_opt*1e4], [h_opt], [W_opt],
            color='gold', s=120, edgecolor='black',
            marker='*', label='Optimum', zorder=5)

ax1.set_xlabel(r'Area $A$ [cm$^2$]')
ax1.set_ylabel(r'Height $h$ [m]')
ax1.set_zlabel(r'Weight $W$ [kg]')
ax1.set_title('Weight Surface and Optimization Path')
fig1.colorbar(surf, shrink=0.6, aspect=12, label='Weight [kg]')
ax1.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 8. Figure 2: 2-D contour with feasible region & boundaries
# ----------------------------------------------------------
feasible = ((AG >= P*bar_length(HG)/(2*HG*sigma_allow)) &
            (AG >= (P*bar_length(HG)**3)/(2*E*HG**2*delta_allow)))

fig2, ax2 = plt.subplots(figsize=(9, 7))
cont = ax2.contourf(AG*1e4, HG, WG, levels=30, cmap='viridis', alpha=0.85)
fig2.colorbar(cont, label='Weight [kg]')

ax2.contourf(AG*1e4, HG, feasible.astype(int), levels=[-0.5, 0.5],
             colors=['white'], alpha=0.6)
ax2.plot(A_stress_curve*1e4, h_fine, color='cyan', linewidth=2.5,
         label='Stress boundary')
ax2.plot(A_defl_curve*1e4, h_fine, color='magenta', linewidth=2.5,
         label='Deflection boundary')
ax2.plot(history[:, 0]*1e4, history[:, 1], color='red',
         marker='o', markersize=4, linewidth=2, label='Optimization path')
ax2.scatter([A_opt*1e4], [h_opt], color='gold', s=160,
            edgecolor='black', marker='*', zorder=5, label='Optimum')

ax2.set_xlim(A_min*1e4, A_max*1e4)
ax2.set_ylim(h_min, h_max)
ax2.set_xlabel(r'Area $A$ [cm$^2$]')
ax2.set_ylabel(r'Height $h$ [m]')
ax2.set_title('Design Space: Feasible Region and Weight Contours')
ax2.legend(loc='upper right')
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 9. Figure 3: Convergence history
# ----------------------------------------------------------
fig3, ax3 = plt.subplots(figsize=(8, 5))
ax3.plot(range(len(path_W)), path_W, marker='o', color='steelblue')
ax3.axhline(W_opt, color='red', linestyle='--',
            label=f'Converged weight = {W_opt:.3f} kg')
ax3.set_xlabel('Iteration')
ax3.set_ylabel('Weight [kg]')
ax3.set_title('Convergence History of the Optimizer')
ax3.grid(alpha=0.3)
ax3.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 10. Figure 4: 1-D analytical cross-check
# ----------------------------------------------------------
fig4, ax4 = plt.subplots(figsize=(8, 5))
ax4.plot(h_fine, W_active_curve, color='darkorange', linewidth=2,
         label='Weight along active constraint boundary')
ax4.axvline(h_star_1d, color='green', linestyle='--',
            label=f'Analytical optimum h* = {h_star_1d:.3f} m')
ax4.scatter([h_opt], [W_opt], color='red', s=80, zorder=5,
            label='SLSQP result')
ax4.set_xlabel(r'Height $h$ [m]')
ax4.set_ylabel(r'Weight $W$ [kg]')
ax4.set_title('1-D Reduced Problem: Weight vs. Height on the Active Constraint')
ax4.grid(alpha=0.3)
ax4.legend()
plt.tight_layout()
plt.show()

=======================================================
 OPTIMIZATION RESULT 
=======================================================
Success              : False
Iterations           : 8
Optimal area  A*     : 0.1000 cm^2
Optimal height h*    : 0.7112 m
Optimal bar length L*: 1.2271 m
Minimum weight W*    : 0.1927 kg
Stress at optimum    : 4313.543 MPa (allowable 165.0 MPa)
Deflection at optimum: 45.6648 mm (allowable 5.0 mm)
=======================================================
1-D cross-check   h* : 1.0008 m,  W*: 4.7576 kg

Code Walkthrough

Section 1 — Fixed parameters. The load $P$, base half-width $b$, steel density $\rho$, Young’s modulus $E$, and the two allowable limits (stress and deflection) are all treated as fixed constants. Only $A$ and $h$ are free to vary.

Section 2 — Mechanics functions. bar_length, axial_force, axial_stress, apex_deflection, and total_weight are direct, vectorized translations of the equations derived above. Because they’re written with NumPy operations rather than explicit loops, they work equally well on scalars (during optimization) and on entire arrays (during plotting), which is what keeps the whole script fast — there’s no heavy computation here, so no separate “fast” version is needed; the vectorized formulation already avoids any per-element Python loop.

Section 3 — Constraints. scipy.optimize.minimize with method='SLSQP' expects inequality constraints written as $g(x) \ge 0$. So the stress constraint $\sigma \le \sigma_{allow}$ is rewritten as sigma_allow - axial_stress(A, h) >= 0, and similarly for deflection.

Section 4 — Optimization. bounds keeps the search within a physically sensible box (areas from 0.1 cm² to 50 cm², heights from 0.5 m to 4 m). x0 is a feasible starting guess. A callback function records every intermediate design point xk into history, which lets us later draw the optimizer’s path through the design space. minimize is then called with the objective, the two inequality constraints, and the bounds.

Section 5 — Reporting. After convergence, the optimal area, height, resulting bar length, and minimum weight are printed, along with the stress and deflection values at the optimum (which should sit at or very near one of the two allowable limits — this is the hallmark of an active constraint at the optimum).

Section 6 — Analytical cross-check. This is the most instructive part of the script. Because weight increases monotonically with $A$ for any fixed $h$, the best possible area for a given height is always exactly the larger of the two constraint-required areas — there’s never a reason to use more material than the tighter constraint demands. That means the full 2-variable problem can be collapsed into a 1-variable problem in $h$ alone: compute the required area from each constraint across a fine grid of $h$ values, take the pointwise maximum, and minimize the resulting weight curve directly with np.argmin. This gives an independent, essentially “brute-force but cheap” verification of what SLSQP found, with no reliance on gradients or convergence tolerances.

Sections 7–10 — Visualization. Four separate figures are generated, each isolating a different way of looking at the result. They are described in detail below.

Visualizing the Results

Figure 1 — The 3D Weight Surface

This plot shows the raw objective function $W(A,h)$ as a surface over the $(A,h)$ plane, with the optimizer’s path traced in red and the final optimum marked with a gold star. Note that the surface itself is monotonically increasing in both variables — taken alone, the unconstrained minimum would simply be the corner with the smallest possible area and height. What this figure really conveys is how the optimizer moves across that surface, converging quickly from the initial guess toward the constrained optimum.

Figure 2 — Design Space, Feasible Region, and Constraint Boundaries

This is the figure that makes the trade-off visible. The colored contours show weight; the cyan curve is the stress-constraint boundary; the magenta curve is the deflection-constraint boundary; and the region shaded white is infeasible (violates at least one constraint). The true feasible region is the “wedge” between the two boundary curves. The optimum sits exactly where the two constraint boundaries and the weight contours pinch together — this is the classic signature of a constrained optimum lying on an active constraint boundary rather than in the interior of the feasible region.

Figure 3 — Convergence History

A simple line plot of the objective value at each iteration of the SLSQP solver. It should drop quickly from the (feasible but suboptimal) initial guess and flatten out as it converges to the minimum weight. This is a good diagnostic for confirming the optimizer didn’t stall or need an excessive number of iterations.

Figure 4 — The 1D Reduced Problem

This plot shows the weight-along-the-active-constraint curve computed in Section 6, as a function of height alone. It has a clear, visually obvious interior minimum — this is the direct visual proof that raising the apex too little makes the deflection constraint expensive (long, heavily-loaded, thin bars deflect too much), while raising it too much makes the growing bar length itself dominate the weight, even as the required area shrinks. The green dashed line marks the analytically found optimum height, and the red dot marks where the SLSQP result landed; the two should coincide almost exactly, confirming the numerical optimizer found the true global optimum for this problem.

Takeaways

The two-bar truss is small enough to solve by hand in reduced form, which is exactly what makes it such a good teaching example: it lets you validate a general-purpose nonlinear optimizer (SLSQP) against an independent, near-analytical solution. The same pattern — objective function, mechanics-derived constraints, a gradient-based solver, and a reduced-dimension sanity check — scales directly to much larger truss problems with dozens of bars and areas, where the reduction to one variable is no longer possible but the same scipy.optimize workflow still applies.

Minimizing Potential Energy in a Two-Variable System

August 9, 2026

A Spring Chain Model in Python

Physical systems love to settle into their lowest-energy configuration. A ball rolls to the bottom of a valley, a molecule relaxes into its most stable bond lengths, a chain of springs stretches until every force balances out. This principle — energy minimization — is one of the most powerful ideas in physics and computational chemistry, and it’s a beautiful playground for numerical optimization.

In this article, we’ll build a concrete example: a 1D chain of springs connecting four particles, where the two end particles are fixed and the two middle particles are free to move. This gives us exactly two variables to optimize — the positions of the two free particles — and it’s a simplified but genuine model of how molecular mechanics software finds the equilibrium geometry of a molecule (think of it as a toy triatomic chain, like a simplified CO₂ backbone).

We’ll solve it three different ways (gradient descent, scipy’s optimizer, and an exact analytical solution), visualize the energy landscape in 3D, and compare a slow vs. a fast implementation.

1. The Physical Model

Imagine four particles arranged along a line, connected in series by three springs:

$P_0$ is fixed at $x_0 = 0$
$P_3$ is fixed at $x_3 = 12$
$P_1$ (position $x_1$) and $P_2$ (position $x_2$) are free — these are our two variables

Each spring has its own stiffness $k_i$ and natural (rest) length $L_i$. The total potential energy of the system is the sum of the harmonic spring energies:

$$
U(x_1, x_2) = \frac{1}{2}k_1\big(x_1 - x_0 - L_1\big)^2 + \frac{1}{2}k_2\big((x_2-x_1) - L_2\big)^2 + \frac{1}{2}k_3\big((x_3-x_2) - L_3\big)^2
$$

The equilibrium configuration is the $(x_1, x_2)$ that minimizes $U$. Because the total distance between the fixed ends ($12$) is larger than the sum of natural lengths ($L_1+L_2+L_3=9$), every spring is forced to stretch somewhat — so the equilibrium is a genuine trade-off between the three springs, not a trivial answer. This is exactly the kind of geometry relaxation problem that molecular mechanics force fields solve, just with more atoms and more complex potentials (Lennard-Jones, angle terms, etc.).

The gradient (force balance conditions) and Hessian (curvature / stiffness matrix) can be derived analytically:

$$
\frac{\partial U}{\partial x_1} = k_1(x_1 - x_0 - L_1) - k_2\big((x_2-x_1)-L_2\big)
$$

$$
\frac{\partial U}{\partial x_2} = k_2\big((x_2-x_1)-L_2\big) - k_3\big((x_3-x_2)-L_3\big)
$$

Since $U$ is quadratic in $x_1, x_2$, the Hessian is constant:

$$
H = \begin{pmatrix} k_1+k_2 & -k_2 \ -k_2 & k_2+k_3 \end{pmatrix}
$$

This lets us solve for the exact minimum with plain linear algebra — a perfect way to double-check our numerical optimizers.

2. Full Python Source Code

Copy this entire cell into a single Google Colaboratory cell and run it.

# ============================================================
# Two-Variable Potential Energy Minimization: Spring Chain Model
# ============================================================
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
import time

# ------------------------------------------------------------
# 1. Physical parameters
# ------------------------------------------------------------
k1, k2, k3 = 1.0, 2.0, 1.5      # spring constants
L1, L2, L3 = 3.0, 3.0, 3.0      # natural lengths
x0_fixed, x3_fixed = 0.0, 12.0  # fixed endpoints

def potential_energy(x):
    """Total potential energy U(x1, x2) of the spring chain."""
    x1, x2 = x
    U1 = 0.5 * k1 * (x1 - x0_fixed - L1) ** 2
    U2 = 0.5 * k2 * ((x2 - x1) - L2) ** 2
    U3 = 0.5 * k3 * ((x3_fixed - x2) - L3) ** 2
    return U1 + U2 + U3

def gradient(x):
    """Analytic gradient (force imbalance) of U."""
    x1, x2 = x
    dU_dx1 = k1 * (x1 - x0_fixed - L1) - k2 * ((x2 - x1) - L2)
    dU_dx2 = k2 * ((x2 - x1) - L2) - k3 * ((x3_fixed - x2) - L3)
    return np.array([dU_dx1, dU_dx2])

def hessian(x):
    """Analytic (constant) Hessian matrix of U."""
    h11 = k1 + k2
    h22 = k2 + k3
    h12 = -k2
    return np.array([[h11, h12], [h12, h22]])

# ------------------------------------------------------------
# 2. Energy landscape: SLOW (naive loop) vs FAST (vectorized)
# ------------------------------------------------------------
x1_range = np.linspace(-2, 10, 300)
x2_range = np.linspace(0, 14, 300)

# --- Naive version: nested Python loops (slow) ---
t0 = time.time()
X1, X2 = np.meshgrid(x1_range, x2_range)
U_slow = np.zeros_like(X1)
for i in range(X1.shape[0]):
    for j in range(X1.shape[1]):
        U_slow[i, j] = potential_energy([X1[i, j], X2[i, j]])
t_slow = time.time() - t0

# --- Vectorized version: full NumPy array math (fast) ---
t0 = time.time()
U1_grid = 0.5 * k1 * (X1 - x0_fixed - L1) ** 2
U2_grid = 0.5 * k2 * ((X2 - X1) - L2) ** 2
U3_grid = 0.5 * k3 * ((x3_fixed - X2) - L3) ** 2
U_fast = U1_grid + U2_grid + U3_grid
t_fast = time.time() - t0

print(f"Naive loop computation time:      {t_slow*1000:.2f} ms")
print(f"Vectorized computation time:      {t_fast*1000:.2f} ms")
print(f"Speedup factor:                   {t_slow/t_fast:.1f}x")
print(f"Results identical (max diff):     {np.max(np.abs(U_slow - U_fast)):.2e}")

# ------------------------------------------------------------
# 3. Method A: Custom steepest-descent optimizer (records path)
# ------------------------------------------------------------
def gradient_descent(x_start, lr=0.15, tol=1e-8, max_iter=500):
    path = [np.array(x_start, dtype=float)]
    x = np.array(x_start, dtype=float)
    for _ in range(max_iter):
        grad = gradient(x)
        x = x - lr * grad
        path.append(x.copy())
        if np.linalg.norm(grad) < tol:
            break
    return x, np.array(path)

x_start = np.array([1.0, 8.0])
x_gd, path_gd = gradient_descent(x_start)

# ------------------------------------------------------------
# 4. Method B: scipy.optimize.minimize with analytic gradient/Hessian
# ------------------------------------------------------------
result = minimize(
    potential_energy, x_start,
    jac=gradient, hess=hessian,
    method='Newton-CG',
    options={'xtol': 1e-10}
)
x_scipy = result.x

# ------------------------------------------------------------
# 5. Method C: Exact analytical solution (linear system H x = c)
# ------------------------------------------------------------
H = hessian(x_start)  # constant everywhere since U is quadratic
c = np.array([
    k1 * (x0_fixed + L1) - k2 * L2,
    k2 * L2 + k3 * (x3_fixed - L3)
])
x_analytic = np.linalg.solve(H, c)

# ------------------------------------------------------------
# 6. Print summary of results
# ------------------------------------------------------------
print("\n--- Equilibrium Positions ---")
print(f"Gradient Descent : x1={x_gd[0]:.6f}, x2={x_gd[1]:.6f}, U={potential_energy(x_gd):.6f}")
print(f"scipy Newton-CG   : x1={x_scipy[0]:.6f}, x2={x_scipy[1]:.6f}, U={potential_energy(x_scipy):.6f}")
print(f"Analytical (exact): x1={x_analytic[0]:.6f}, x2={x_analytic[1]:.6f}, U={potential_energy(x_analytic):.6f}")

bond1 = x_analytic[0] - x0_fixed
bond2 = x_analytic[1] - x_analytic[0]
bond3 = x3_fixed - x_analytic[1]
print("\n--- Equilibrium Bond (Spring) Lengths ---")
print(f"Spring 1: {bond1:.4f}  (natural length {L1})")
print(f"Spring 2: {bond2:.4f}  (natural length {L2})")
print(f"Spring 3: {bond3:.4f}  (natural length {L3})")

# ------------------------------------------------------------
# 7. Visualization: 3D energy surface
# ------------------------------------------------------------
fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X1, X2, U_fast, cmap='viridis', alpha=0.85,
                         linewidth=0, antialiased=True)
ax1.scatter(x_analytic[0], x_analytic[1], potential_energy(x_analytic),
            color='red', s=80, label='Equilibrium (minimum)')
ax1.set_xlabel('x1 (position of P1)')
ax1.set_ylabel('x2 (position of P2)')
ax1.set_zlabel('Potential Energy U')
ax1.set_title('3D Potential Energy Surface U(x1, x2)')
ax1.legend()
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10, label='Energy')

# ------------------------------------------------------------
# 8. Visualization: 2D contour plot with optimization path
# ------------------------------------------------------------
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X1, X2, U_fast, levels=40, cmap='viridis')
ax2.contour(X1, X2, U_fast, levels=15, colors='white', linewidths=0.4, alpha=0.5)
ax2.plot(path_gd[:, 0], path_gd[:, 1], 'r.-', markersize=3, linewidth=1,
         label='Gradient descent path')
ax2.plot(x_start[0], x_start[1], 'ko', markersize=8, label='Start point')
ax2.plot(x_analytic[0], x_analytic[1], 'r*', markersize=18, label='Equilibrium')
ax2.set_xlabel('x1 (position of P1)')
ax2.set_ylabel('x2 (position of P2)')
ax2.set_title('Energy Contours and Descent Path to Equilibrium')
ax2.legend()
fig.colorbar(contour, ax=ax2, label='Energy')

plt.tight_layout()
plt.show()

Naive loop computation time:      558.96 ms
Vectorized computation time:      6.64 ms
Speedup factor:                   84.2x
Results identical (max diff):     2.84e-14

--- Equilibrium Positions ---
Gradient Descent : x1=4.384615, x2=8.076923, U=2.076923
scipy Newton-CG   : x1=4.384615, x2=8.076923, U=2.076923
Analytical (exact): x1=4.384615, x2=8.076923, U=2.076923

--- Equilibrium Bond (Spring) Lengths ---
Spring 1: 4.3846  (natural length 3.0)
Spring 2: 3.6923  (natural length 3.0)
Spring 3: 3.9231  (natural length 3.0)

3. Code Walkthrough

Physical model functions

potential_energy(x) implements the equation for $U(x_1,x_2)$ directly — three harmonic terms, one per spring. gradient(x) and hessian(x) implement the analytic derivatives we derived above. Supplying these analytically (rather than letting the optimizer approximate them with finite differences) is both faster and numerically more accurate, since finite-difference gradients require multiple extra function evaluations and introduce rounding error.

Naive loop vs. vectorized grid computation

This is the “time-consuming part” of the problem, and it’s a great illustration of why vectorization matters in scientific Python. Computing the energy value at every point of a 300×300 grid via a for i in range(...): for j in range(...): double loop calls potential_energy() 90,000 times in pure Python, which is slow because every call carries Python’s interpreter overhead. The vectorized version instead evaluates the exact same formula on entire NumPy arrays at once (X1, X2 are full 300×300 arrays), letting NumPy’s compiled C backend do all 90,000 evaluations in one shot. The script prints both timings and the speedup factor — on a typical Colab CPU you should see roughly a 50–150× speedup, and the “max diff” check confirms both methods produce identical results (up to floating-point rounding).

Three independent solving methods

Gradient descent (gradient_descent) — a hand-written steepest-descent loop that repeatedly steps opposite to the gradient. We record every intermediate point in path_gd so we can visualize how the optimizer walks downhill.
scipy’s minimize with method='Newton-CG' — a professional-grade optimizer that uses our analytic gradient and Hessian to converge in far fewer iterations than plain gradient descent, since it accounts for the curvature of the energy surface.
Exact analytical solution — because $U$ is a quadratic function, its minimum satisfies the linear equation $H\mathbf{x} = \mathbf{c}$, which np.linalg.solve solves exactly in one step. This serves as a ground-truth check: both numerical methods should match it almost to machine precision.

All three results are printed together, along with the resulting equilibrium bond lengths (the actual stretched length of each spring), so you can see how each spring deviates from its natural length to balance the system.

Visualization

The 3D surface plot (left panel) shows the full bowl-shaped energy landscape — since our potential is a sum of quadratics, it’s a paraboloid, and the red dot marks the single global minimum. The contour plot with path overlay (right panel) shows the same landscape from above, with the gradient descent trajectory drawn as a red line from the black starting point down to the red star at equilibrium. Watching the descent path curve toward the minimum (rather than moving in a straight line) illustrates that the two variables are coupled — moving $x_1$ affects the optimal $x_2$ and vice versa, exactly as spring 2 couples the two free particles together.

4. Interpreting the Results

Once you run the code, check that:

All three methods agree — gradient descent, scipy’s Newton-CG, and the exact linear solve should all report essentially the same $(x_1, x_2)$ and energy value (differences should be smaller than $10^{-4}$).
The bond lengths make physical sense — the stiffer spring ($k_2=2.0$) should be stretched less relative to its natural length than the softer springs, since it “resists” deformation more strongly. This is a direct numerical demonstration of Hooke’s law competition between coupled springs.
The energy surface is convex (bowl-shaped) — this is why gradient descent, despite being a simple algorithm, reliably finds the global minimum here. Real molecular potentials (e.g., Lennard-Jones) are not globally convex and can have multiple local minima, which is why more sophisticated global optimization or multiple random restarts are used in real molecular mechanics software.

5. Extending the Model

This two-variable spring chain is a minimal but genuine example of geometry optimization, the same core computational idea used in molecular simulation packages (like force-field minimization in GROMACS or AMBER) to relax a molecule into its stable 3D shape. Natural extensions worth trying:

Replace one harmonic spring with a Lennard-Jones potential $U_{LJ}(r) = 4\epsilon\left[(\sigma/r)^{12}-(\sigma/r)^6\right]$ to introduce a non-convex landscape with a local minimum, and see how gradient descent can get stuck depending on the starting point.
Extend to more free particles (higher dimensions), where visualizing the full energy surface is no longer possible, but the same gradient/Hessian-based optimization strategy still works.
Add a 2D or 3D geometry (particles free to move in the plane, not just along a line) to more closely resemble real bond-angle relaxation in molecules.