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

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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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.