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.

Minimizing the Goldstein-Price Function with Python

August 8, 2026

A Global Optimization Walkthrough

Global optimization is one of those problems that looks simple on paper but gets nasty fast once your objective function has multiple local minima. A classic benchmark for testing optimization algorithms is the Goldstein-Price function, a two-dimensional function riddled with local minima that only reveals its single global minimum after careful searching.

In this article, we’ll define the function, minimize it with Python and SciPy, visualize the results in 3D, and break down every part of the code so you understand exactly what’s happening under the hood.

What Is the Goldstein-Price Function?

The Goldstein-Price function is defined as:

$$
f(x, y) = \left[1 + (x + y + 1)^2 \left(19 - 14x + 3x^2 - 14y + 6xy + 3y^2\right)\right] \times \left[30 + (2x - 3y)^2 \left(18 - 32x + 12x^2 + 48y - 36xy + 27y^2\right)\right]
$$

It is typically evaluated over the domain:

$$
x, y \in [-2, 2]
$$

The function has a known global minimum of:

$$
f(0, -1) = 3
$$

What makes this function a good stress test is that it’s not smooth and simple — it has several local minima and a huge dynamic range (values swing from 3 up to over a million within the search domain), which means naive gradient-based methods starting from a bad initial guess can easily get stuck.

To reliably find the global minimum, we’ll use a two-step approach:

Differential Evolution — a population-based global optimization algorithm from SciPy that doesn’t require gradients and is good at escaping local minima.
Nelder-Mead — a local simplex-based method, used here to confirm/refine a result from a specific starting point, illustrating how a local method can get “close enough” but benefits from a good starting guess.

The Source Code

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

# --- Define the Goldstein-Price function (scalar version, for optimizers) ---
def goldstein_price(v):
    x, y = v
    term1 = 1 + (x + y + 1)**2 * (19 - 14*x + 3*x**2 - 14*y + 6*x*y + 3*y**2)
    term2 = 30 + (2*x - 3*y)**2 * (18 - 32*x + 12*x**2 + 48*y - 36*x*y + 27*y**2)
    return term1 * term2

# --- Vectorized version (for fast grid evaluation used in plotting) ---
def goldstein_price_vec(X, Y):
    term1 = 1 + (X + Y + 1)**2 * (19 - 14*X + 3*X**2 - 14*Y + 6*X*Y + 3*Y**2)
    term2 = 30 + (2*X - 3*Y)**2 * (18 - 32*X + 12*X**2 + 48*Y - 36*X*Y + 27*Y**2)
    return term1 * term2

# --- Search domain ---
bounds = [(-2, 2), (-2, 2)]

# --- Step 1: Global search using Differential Evolution ---
start_time = time.time()
result_de = differential_evolution(
    goldstein_price,
    bounds,
    tol=1e-12,
    seed=42
)
elapsed = time.time() - start_time

# --- Step 2: Local refinement from an arbitrary starting point ---
result_local = minimize(
    goldstein_price,
    x0=[0.5, -0.5],
    method='Nelder-Mead'
)

# --- Report results ---
print("=== Differential Evolution (Global Search) ===")
print(f"  x = {result_de.x[0]:.6f}, y = {result_de.x[1]:.6f}")
print(f"  f(x, y) = {result_de.fun:.6f}")
print(f"  Elapsed time: {elapsed:.4f} sec")

print("\n=== Nelder-Mead (Local Search from [0.5, -0.5]) ===")
print(f"  x = {result_local.x[0]:.6f}, y = {result_local.x[1]:.6f}")
print(f"  f(x, y) = {result_local.fun:.6f}")

=== Differential Evolution (Global Search) ===
  x = -0.000000, y = -1.000000
  f(x, y) = 3.000000
  Elapsed time: 0.2694 sec

=== Nelder-Mead (Local Search from [0.5, -0.5]) ===
  x = 0.000019, y = -0.999989
  f(x, y) = 3.000000

Code Walkthrough

The two function definitions. goldstein_price(v) takes a single vector v = [x, y] and returns a scalar — this is the signature SciPy’s optimizers expect. goldstein_price_vec(X, Y) does the exact same math but accepts NumPy arrays (meshgrids) directly, using array broadcasting instead of loops. This second version is purely for speed when we evaluate the function over thousands of grid points for plotting — looping point-by-point in Python would be dramatically slower, so vectorizing with NumPy’s broadcasted arithmetic keeps the whole grid evaluation to a fraction of a second.

differential_evolution is a stochastic, population-based global optimizer. Instead of following a gradient, it maintains a population of candidate solutions and evolves them generation by generation using mutation and crossover, which makes it well-suited to functions like Goldstein-Price that have multiple local minima capable of trapping gradient-based methods. We pass tol=1e-12 to tighten the convergence criterion and seed=42 to make the result reproducible.

minimize(..., method='Nelder-Mead') performs a local simplex search starting from [0.5, -0.5]. Because Nelder-Mead doesn’t require derivatives, it works even on non-smooth objective functions, but it only guarantees convergence to whichever minimum is closest to the starting point — not necessarily the global one. Comparing its output to the differential evolution result is a good way to see how starting-point-dependent local optimization can be.

When you run this, both methods should converge to essentially the same point, (x, y) ≈ (0, -1) with f(x, y) ≈ 3, which matches the function’s known global minimum. The differential evolution run also typically finishes in well under a second, since the search domain is small and two-dimensional.

Visualizing the Function

Because Goldstein-Price spans values from 3 to over a million across the domain, plotting it on a linear scale would flatten the interesting structure near the minimum into an indistinguishable blob. The fix is to plot the base-10 logarithm of the function instead, which compresses the huge dynamic range into something visually readable while preserving the location of the minimum.

# --- Build a grid over the search domain ---
x = np.linspace(-2, 2, 200)
y = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(x, y)
Z = goldstein_price_vec(X, Y)
Z_log = np.log10(Z)  # compress the huge dynamic range for visualization

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

# --- 3D surface plot (log scale) ---
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z_log, cmap='viridis', alpha=0.9, edgecolor='none')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('log10(f(x, y))')
ax1.set_title('Goldstein-Price Function (log scale)')
ax1.scatter(
    [result_de.x[0]], [result_de.x[1]], [np.log10(result_de.fun)],
    color='red', s=100, marker='*', label='Global Minimum'
)
ax1.legend()
fig.colorbar(surf, ax=ax1, shrink=0.5)

# --- 2D contour plot (log scale) ---
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X, Y, Z_log, levels=50, cmap='viridis')
ax2.contour(X, Y, Z_log, levels=20, colors='white', linewidths=0.3, alpha=0.5)
ax2.plot(result_de.x[0], result_de.x[1], 'r*', markersize=20, label='Global Minimum')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('Contour Plot (log scale)')
ax2.legend()
fig.colorbar(contour, ax=ax2)

plt.tight_layout()
plt.show()

Graph Explanation

Left panel — 3D surface plot. This shows the shape of log10(f(x, y)) across the search domain. You can clearly see a steep, funnel-shaped basin near (0, -1), marked with a red star, which is the global minimum. Away from that basin, the surface rises sharply into several ridges and bumps — these are the local minima and saddle regions that make this function a challenging optimization benchmark. The logarithmic z-axis is what makes both the deep basin and the surrounding terrain visible in the same plot; without it, the basin would be invisible next to the much larger values elsewhere in the domain.

Right panel — 2D contour plot. This is essentially a bird’s-eye view of the same log-scaled surface, where color represents function value (darker purple = lower, i.e. closer to the minimum) and the white contour lines trace constant-value bands. The red star again marks the global minimum at (0, -1). Notice how tightly the contour lines bunch up around that point — that steepness is exactly why gradient-based optimizers can converge quickly once they’re near it, but can also get misled by the other local minima elsewhere on the map if they start too far away.

Together, these two plots make it easy to visually confirm what the numerical optimizers already told us: the true minimum sits at (0, -1) with a function value of 3, nestled at the bottom of a narrow, steep-sided basin surrounded by a much bumpier landscape.

Taming the Beale Function

August 7, 2026

A Journey Through a Deceptively Simple Optimization Landscape

Introduction

Among the classic benchmark functions used to test optimization algorithms, the Beale function holds a special place. It looks innocent enough on paper — just a sum of three squared terms — but its landscape hides sharp, narrow valleys that can trap naive gradient-based solvers. In this article, we’ll dissect the Beale function mathematically, implement a robust two-stage optimization pipeline in Python, and visualize the entire search process in 3D.

The Mathematics of the Beale Function

The Beale function is defined as:

$$
f(x, y) = (1.5 - x + xy)^2 + (2.25 - x + xy^2)^2 + (2.625 - x + xy^3)^2
$$

It is evaluated over the domain $x, y \in [-4.5, 4.5]$, and its global minimum is known analytically:

$$
f(3, ; 0.5) = 0
$$

What makes this function tricky is the steep, curved valley leading toward the minimum, combined with flat plateaus and extremely large function values near the corners of the domain (values can exceed $10^5$). This asymmetry means a single fixed-step gradient method can easily overshoot or stall, which is exactly why a hybrid global-then-local strategy is the right tool for the job.

Optimization Strategy

Rather than relying purely on a local method (which is highly sensitive to the starting point) or purely on a slow global method, we combine two techniques:

Differential Evolution (DE) — a population-based global optimizer that explores the whole domain without needing gradient information, giving us a good approximate basin.
L-BFGS-B — a fast, gradient-based local refinement step that polishes the DE result down to near machine precision.

This two-stage approach is both fast and reliable: DE avoids getting stuck in the wrong region, while L-BFGS-B converges to high accuracy in a handful of iterations once we’re already close.

Source Code

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

plt.style.use('dark_background')

# -----------------------------
# Beale function definition
# -----------------------------
def beale(params):
    x, y = params
    term1 = (1.5 - x + x * y) ** 2
    term2 = (2.25 - x + x * y ** 2) ** 2
    term3 = (2.625 - x + x * y ** 3) ** 2
    return term1 + term2 + term3

# -----------------------------
# History tracking for visualization
# -----------------------------
convergence_history = []
positions_history = []

def de_callback(xk, convergence):
    positions_history.append(np.array(xk, dtype=float))
    convergence_history.append(beale(xk))

bounds = [(-4.5, 4.5), (-4.5, 4.5)]

# -----------------------------
# Stage 1: Global search with Differential Evolution
# -----------------------------
start_time = time.time()
result_de = differential_evolution(
    beale,
    bounds,
    seed=42,
    maxiter=300,
    tol=1e-12,
    mutation=(0.5, 1.5),
    recombination=0.7,
    callback=de_callback,
    polish=False
)
de_elapsed = time.time() - start_time

# -----------------------------
# Stage 2: Local refinement with L-BFGS-B
# -----------------------------
start_time = time.time()
result_local = minimize(
    beale,
    result_de.x,
    method='L-BFGS-B',
    bounds=bounds,
    tol=1e-15
)
local_elapsed = time.time() - start_time

true_min = np.array([3.0, 0.5])
distance_to_true = np.linalg.norm(result_local.x - true_min)

print("=" * 60)
print("Beale Function Minimization Results")
print("=" * 60)
print(f"[Stage 1] Differential Evolution")
print(f"  x = {result_de.x[0]:.10f}")
print(f"  y = {result_de.x[1]:.10f}")
print(f"  f(x,y) = {result_de.fun:.15f}")
print(f"  Generations: {result_de.nit}, Time: {de_elapsed:.4f}s")
print("-" * 60)
print(f"[Stage 2] L-BFGS-B Local Refinement")
print(f"  x = {result_local.x[0]:.10f}")
print(f"  y = {result_local.x[1]:.10f}")
print(f"  f(x,y) = {result_local.fun:.2e}")
print(f"  Time: {local_elapsed:.4f}s")
print("-" * 60)
print(f"Known global minimum: (3.0, 0.5), f = 0.0")
print(f"Distance to true minimum: {distance_to_true:.2e}")
print(f"Total optimization time: {de_elapsed + local_elapsed:.4f}s")
print("=" * 60)

============================================================
Beale Function Minimization Results
============================================================
[Stage 1] Differential Evolution
  x = 3.0000000000
  y = 0.5000000000
  f(x,y) = 0.000000000000000
  Generations: 146, Time: 0.6139s
------------------------------------------------------------
[Stage 2] L-BFGS-B Local Refinement
  x = 3.0000000000
  y = 0.5000000000
  f(x,y) = 3.20e-31
  Time: 0.0116s
------------------------------------------------------------
Known global minimum: (3.0, 0.5), f = 0.0
Distance to true minimum: 4.97e-16
Total optimization time: 0.6255s
============================================================

Code Walkthrough

The objective function. beale() implements the three-term formula exactly as written mathematically. Because scipy.optimize passes parameters as a single array, we unpack x, y = params at the top of the function.

History tracking. The de_callback function is invoked after every generation of the differential evolution algorithm. We store both the candidate position (positions_history) and its function value (convergence_history), which lets us later draw the search trajectory over the contour map and plot the convergence curve.

Stage 1 — global search. differential_evolution maintains a population of candidate solutions and evolves them via mutation and crossover, requiring no gradient information. We disable the built-in polish step (polish=False) because we handle refinement ourselves in Stage 2, giving us explicit control and separate timing for each stage.

Stage 2 — local refinement. L-BFGS-B takes the best point found by DE and rapidly converges toward the true minimum using quasi-Newton updates with box constraints, typically needing only a few iterations since the starting point is already close to the optimum.

Why this is already fast. For a 2D problem like this, DE with 300 generations and a small population converges in well under a second, and L-BFGS-B refinement adds negligible overhead. No further acceleration (e.g., multiprocessing or vectorized batch evaluation) is necessary here — the bottleneck for 2D benchmark functions is never raw computation time.

Visualizing the Optimization Landscape

Numbers alone don’t do justice to how treacherous the Beale function’s landscape really is. Let’s render it in three complementary views: a 3D surface, a contour map with the search trajectory overlaid, and the convergence curve.

from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(16, 12))
fig.patch.set_facecolor('#0d1117')
gs = GridSpec(2, 2, width_ratios=[1.3, 1], height_ratios=[1, 1],
              hspace=0.35, wspace=0.3, figure=fig)

# -----------------------------
# Prepare grid data
# -----------------------------
x_range = np.linspace(-4.5, 4.5, 220)
y_range = np.linspace(-4.5, 4.5, 220)
X, Y = np.meshgrid(x_range, y_range)
Z = (1.5 - X + X * Y) ** 2 + (2.25 - X + X * Y ** 2) ** 2 + (2.625 - X + X * Y ** 3) ** 2
Z_log = np.log1p(Z)  # log-scale to compress the huge dynamic range

# -----------------------------
# Panel 1: 3D Surface
# -----------------------------
ax1 = fig.add_subplot(gs[:, 0], projection='3d')
ax1.set_facecolor('#0d1117')
surf = ax1.plot_surface(X, Y, Z_log, cmap='inferno', edgecolor='none',
                         alpha=0.92, antialiased=True, rcount=150, ccount=150)
ax1.scatter([3.0], [0.5], [np.log1p(beale([3.0, 0.5]))],
            color='cyan', s=140, marker='*', edgecolor='white',
            linewidth=0.5, label='Global Minimum (3, 0.5)', depthshade=False)
ax1.scatter([result_local.x[0]], [result_local.x[1]], [np.log1p(result_local.fun)],
            color='lime', s=80, marker='o', edgecolor='white',
            linewidth=0.5, label='Found Minimum', depthshade=False)
ax1.set_xlabel('x', color='white', labelpad=10)
ax1.set_ylabel('y', color='white', labelpad=10)
ax1.set_zlabel('log(1 + f(x,y))', color='white', labelpad=10)
ax1.set_title('Beale Function Surface (log scale)', color='white', fontsize=13, pad=15)
ax1.tick_params(colors='white')
ax1.legend(loc='upper left', facecolor='#161b22', edgecolor='gray', labelcolor='white', fontsize=9)
ax1.view_init(elev=28, azim=-55)
fig.colorbar(surf, ax=ax1, shrink=0.5, pad=0.08)

# -----------------------------
# Panel 2: Contour + search trace
# -----------------------------
ax2 = fig.add_subplot(gs[0, 1])
ax2.set_facecolor('#0d1117')
ax2.contourf(X, Y, Z_log, levels=50, cmap='inferno')
positions_arr = np.array(positions_history)
if len(positions_arr) > 0:
    ax2.plot(positions_arr[:, 0], positions_arr[:, 1],
              color='cyan', linewidth=1.2, alpha=0.8, label='DE search trace')
ax2.scatter([3.0], [0.5], color='lime', s=160, marker='*',
            edgecolor='white', linewidth=0.5, label='Global Minimum', zorder=5)
ax2.set_xlabel('x', color='white')
ax2.set_ylabel('y', color='white')
ax2.set_title('Contour Map with Search Path', color='white', fontsize=13)
ax2.tick_params(colors='white')
ax2.legend(facecolor='#161b22', edgecolor='gray', labelcolor='white', fontsize=8, loc='upper left')

# -----------------------------
# Panel 3: Convergence curve
# -----------------------------
ax3 = fig.add_subplot(gs[1, 1])
ax3.set_facecolor('#0d1117')
conv_arr = np.array(convergence_history) + 1e-16
ax3.semilogy(range(len(conv_arr)), conv_arr, color='#ff6b6b', linewidth=2)
ax3.set_xlabel('Generation', color='white')
ax3.set_ylabel('f(x, y)  [log scale]', color='white')
ax3.set_title('Differential Evolution Convergence', color='white', fontsize=13)
ax3.grid(True, alpha=0.2)
ax3.tick_params(colors='white')

plt.tight_layout()
plt.show()

Interpreting the Visualization

The 3D surface plot (left panel) reveals why this function is such a good stress test: the terrain is dominated by a steep-walled basin curling from the upper-left toward the bottom-right, with the true minimum sitting at the bottom of a long, narrow trough. We plot $\log(1+f)$ instead of raw $f$ because the true function values span from $0$ to over $10^5$ across the domain — without the log transform, the entire interesting region near the minimum would be crushed flat and invisible.

The contour map with search trace (top right) shows the differential evolution population converging generation by generation. Notice how the cyan trace initially explores broadly across the domain before narrowing sharply into the valley containing the star-marked global minimum — this is the hallmark of a well-functioning global optimizer.

The convergence curve (bottom right) plots the best function value found at each generation on a log scale. The steep initial drop reflects DE quickly identifying the correct basin, while the long, near-flat tail shows the algorithm fine-tuning within that basin before we hand off to L-BFGS-B for the final high-precision polish.

Conclusion

The Beale function is a textbook example of why optimization algorithm choice matters as much as the algorithm’s raw speed. A purely local method risks getting misled by the function’s sharp curvature and vast scale differences, while a purely global method wastes time achieving precision it isn’t designed for. By pairing differential evolution’s broad exploration with L-BFGS-B’s precise convergence, we consistently land within machine-precision distance of the true minimum at $(3, 0.5)$ — a pattern that generalizes well beyond this single benchmark function to many real-world non-convex optimization problems.