Probabilistic Generative Models

A latency dashboard reports a mean of 60 milliseconds across production traffic. The operations dashboard displays a green status, but when you inspect the raw request logs [28, 30, 32, 86, 90, 94], not a single request actually took 60 ms. Three queries finished in roughly 30 ms (in-memory cache hits), while the remaining three stalled near 90 ms (unindexed disk queries).

Fitting a single standard Gaussian to this sample produces a bell curve centered at 60 ms: the exact region where the probability of observing an actual request is near zero. A hard clustering algorithm like K-means assigns each point to a nearest centroid, but it can't quantify how confident it feels about boundary observations. A Gaussian mixture model resolves this by modeling the data as a combination of multiple hidden sub-populations, keeping uncertain membership visible through probability distributions.^{[1]Reference 1Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

Generative versus discriminative modeling

Machine learning systems answer two distinct modeling questions depending on whether they estimate conditional or joint distributions.

Discriminative models estimate the conditional distribution $P(Y \mid X)$. Given observed features $X$ (like request latency or query size), a discriminative classifier determines the probability of class label $Y$ (such as cache hit versus database fallback). Discriminative models focus strictly on the decision boundary separating classes. If you feed them corrupted inputs or impossible latencies like -500 ms, they still output confident class probabilities because they never modeled what valid inputs look like. They can't generate new data points, evaluate input plausibility, or identify missing features.

Generative models estimate how data came to exist in the world. They learn either the joint distribution $P(X, Y) = P(Y)P(X \mid Y)$ or, in the unsupervised setting, the marginal data distribution $P(X)$. By capturing the data-generating mechanism, a generative model evaluates whether a new observation looks probable under normal traffic, generates synthetic samples by sampling from the learned distribution, and accounts for unobserved mechanisms.

In real-world request logs, we don't have labeled class tags $Y$ telling us which infrastructure component served each query. We only observe raw latencies $X$. To explain multi-modal observations without human labels, generative models introduce unobserved hidden variables: latent variables.

Latent variables and the intractable log-sum

Start with the six observed request latencies: [28, 30, 32, 86, 90, 94]. Their arithmetic mean is 60 ms, but the numbers naturally cluster into two distinct groups.

A probabilistic generative model formulates a two-stage story for each observation:

First, nature selects an unobserved latent component $z_i = k \in \{1, \dots, K\}$ with prior probability $P(z_i = k) = \pi_k$, where $\sum_{k=1}^K \pi_k = 1$ and $\pi_k \ge 0$. Second, conditioned on that choice, nature samples the observable latency $x_i$ from that component's density $p(x_i \mid z_i = k, \theta) = \mathcal{N}(x_i \mid \mu_k, \sigma_k^2)$.

Because the request log only records the observed values $X = \{x_1, \dots, x_n\}$ while the latent choices $Z = \{z_1, \dots, z_n\}$ remain hidden, evaluating the probability of an observed point requires marginalizing over all possible latent assignments:

$$p(x_i \mid \theta) = \sum_{k=1}^{K} p(x_i, z_i = k \mid \theta) = \sum_{k=1}^{K} \pi_k\,\mathcal{N}(x_i \mid \mu_k, \sigma_k^2).$$

Assuming requests arrive independently, maximum likelihood estimation aims to find parameters $\theta = \{\pi_k, \mu_k, \sigma_k^2\}_{k=1}^K$ that maximize the joint probability of all observed data. Taking the logarithm converts the product of independent probabilities into a sum:

$$\ell(\theta) = \log p(X \mid \theta) = \sum_{i=1}^{n} \log \left( \sum_{k=1}^{K} \pi_k\,\mathcal{N}(x_i \mid \mu_k, \sigma_k^2) \right).$$

Notice the mathematical barrier: the summation across components sits inside the logarithm.

In standard maximum likelihood without latent variables, taking the logarithm moves directly onto exponential family densities. The log cancels the exponential term, creating simple linear or quadratic expressions whose derivatives set cleanly to zero. With latent variables, the sum inside the log prevents the logarithm from reaching individual Gaussian terms. All parameters across all components become nonlinearly coupled. Setting $\nabla_\theta \ell(\theta) = 0$ yields no closed-form analytical solution. If $Z$ were continuous rather than discrete, that inner sum would become an intractable high-dimensional integral.

Soft responsibilities in Gaussian mixtures

If we knew which latent component generated each request, fitting each Gaussian would be straightforward: compute the sample mean and variance for each group separately. If we already had the true parameters $\theta$, assigning requests to components would also be straightforward Bayes' rule.

Because we have neither, we start by evaluating provisional parameters. Consider two provisional components: a fast component at $\mu_1 = 30$ ms, a slow component at $\mu_2 = 90$ ms, equal standard deviations $\sigma_1 = \sigma_2 = 20$ ms, and equal prior weights $\pi_1 = \pi_2 = 0.5$.

Now evaluate a boundary request that took 50 ms. The fast Gaussian's density at 50 ms is roughly 0.01210 per ms; the slow Gaussian's density is roughly 0.00270 per ms. Weighting each contribution by prior probability 0.5 gives 0.00605 and 0.00135.

Normalizing these two numbers yields the posterior probability that component $k$ generated observation $x_i$, known as the component's responsibility $\gamma_{ik}$:^{[1]Reference 1Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

$$\gamma_{ik} = P(z_i = k \mid x_i, \theta) = \frac{\pi_k\,\mathcal{N}(x_i \mid \mu_k, \sigma_k^2)}{\sum_{j=1}^{K} \pi_j\,\mathcal{N}(x_i \mid \mu_j, \sigma_j^2)}.$$

For our 50 ms request, the fast responsibility is $0.00605 / (0.00605 + 0.00135) \approx 0.818$, and the slow responsibility is $0.00135 / 0.00740 \approx 0.182$. The next code snippet calculates these values without intermediate rounding:

The fast component receives greater responsibility here because its mean is closer and variances are equal. If components had unequal spreads or priors, the closer mean wouldn't automatically win.

Contrast this soft probabilistic assignment with hard K-means clustering from Clustering, PCA, and Representation Learning. K-means forces an absolute 0 or 1 choice based on Euclidean distance, discarding ambiguity. Mathematically, K-means is the limiting case of a Gaussian mixture model where all component covariances are identical and spherical ($\Sigma_k = \sigma^2 I$) as variance shrinks to zero ($\sigma^2 \to 0$). In that zero-variance limit, the softmax posterior collapses into a hard argmax. By retaining non-zero variance, a Gaussian mixture preserves the real uncertainty: a 50 ms latency leans fast, but still carries an 18.2% probability of being an unusually speedy database query.

Responsibilities are discrete probabilities summing to one for each observation ($\sum_k \gamma_{ik} = 1$). In contrast, $p(x \mid \theta)$ is a continuous probability density, whose values can exceed 1.0 for narrow spreads and whose integral over an interval gives probability.

A 50-millisecond request receives fast-component responsibility 0.818. What happens to the remaining 0.182, and does the larger number prove that an in-memory cache handled the request?

Answer

The remaining 0.182 belongs to the slow component, so both responsibilities sum to exactly one. The 0.818 score is a conditional probability under the assumed statistical model; no hardware traces or queue states were measured, so it doesn't establish physical causation.

Expectation-Maximization as lower-bound ascent

We need responsibilities to learn parameters, but we need parameters to calculate responsibilities. The Expectation-Maximization (EM) algorithm, introduced by Dempster, Laird, and Rubin (1977)^{[2]Reference 2Maximum Likelihood from Incomplete Data via the EM Algorithmhttps://doi.org/10.1111/j.2517-6161.1977.tb01600.x}, resolves this mutual dependency by alternating two coordinated steps.

Rather than attacking the non-convex log-sum directly, EM introduces an arbitrary probability distribution $q(Z)$ over the latent variables and constructs a tractable lower bound: the Evidence Lower Bound (ELBO).

Let's derive the lower bound through Jensen's inequality. For any concave function $f$, Jensen's inequality guarantees that the function of an expectation is at least the expectation of the function: $f(\mathbb{E}[Y]) \ge \mathbb{E}[f(Y)]$. Because the logarithm function is strictly concave, multiplying and dividing by $q(Z)$ yields:

$$\log p(X \mid \theta) = \log \sum_{Z} q(Z) \frac{p(X, Z \mid \theta)}{q(Z)} = \log \mathbb{E}_q \left[ \frac{p(X, Z \mid \theta)}{q(Z)} \right] \ge \mathbb{E}_q \left[ \log \frac{p(X, Z \mid \theta)}{q(Z)} \right].$$

Expanding that expectation defines the Evidence Lower Bound $\mathcal{L}(q, \theta)$:^{[3]Reference 3CS229 Machine Learning Notes: The EM Algorithmhttps://cs229.stanford.edu/notes2022fall/main_notes.pdf}

$$\mathcal{L}(q, \theta) \equiv \sum_{Z} q(Z) \log \frac{p(X, Z \mid \theta)}{q(Z)}.$$

We can also express the relationship as an exact algebraic decomposition:

$$\log p(X \mid \theta) = \mathcal{L}(q, \theta) + D_{\mathrm{KL}}(q(Z) \parallel p(Z \mid X, \theta)),$$

where $D_{\mathrm{KL}}(q \parallel p)$ is the Kullback-Leibler divergence measuring the divergence between our distribution $q(Z)$ and the true posterior $p(Z \mid X, \theta)$:

$$D_{\mathrm{KL}}(q(Z) \parallel p(Z \mid X, \theta)) = \sum_{Z} q(Z) \log \frac{q(Z)}{p(Z \mid X, \theta)}.$$

Because KL divergence is always non-negative ($D_{\mathrm{KL}} \ge 0$) and equals zero if and only if $q(Z) = p(Z \mid X, \theta)$, $\mathcal{L}(q, \theta)$ forms a valid lower bound on the true marginal log likelihood $\log p(X \mid \theta)$.

The algorithm climbs this geometry through coordinate ascent:

1. E-step (Expectation / Bound Tightening): Hold model parameters fixed at $\theta^{(t)}$. Set $q(Z) = p(Z \mid X, \theta^{(t)})$, computing the conditional responsibilities $\gamma_{ik}$. This forces $D_{\mathrm{KL}}(q \parallel p) = 0$, causing the lower bound $\mathcal{L}(q, \theta)$ to touch the true log likelihood curve tangentially at $\theta^{(t)}$:

   $$\mathcal{L}(q, \theta^{(t)}) = \log p(X \mid \theta^{(t)}).$$

2. M-step (Maximization / Parameter Optimization): Freeze $q(Z)$ and maximize $\mathcal{L}(q, \theta)$ with respect to $\theta$:

   $$\theta^{(t+1)} = \arg\max_\theta \mathcal{L}(q, \theta) = \arg\max_\theta \sum_{Z} q(Z) \log p(X, Z \mid \theta).$$

   Because the logarithm now operates directly on the complete-data joint probability $p(X, Z \mid \theta)$, the sum inside the log is gone! For Gaussian components, this yields exact closed-form updates.

3. Monotonicity guarantee: Combining both steps confirms that likelihood never decreases:^{[2]Reference 2Maximum Likelihood from Incomplete Data via the EM Algorithmhttps://doi.org/10.1111/j.2517-6161.1977.tb01600.x}

   $$\log p(X \mid \theta^{(t+1)}) \ge \mathcal{L}(q, \theta^{(t+1)}) \ge \mathcal{L}(q, \theta^{(t)}) = \log p(X \mid \theta^{(t)}).$$

Jensen's inequality ensures the ELBO never exceeds true likelihood. Maximizing that bound over $\theta$ in the M-step guarantees progress, while the E-step closes the divergence gap at the current parameter. Together, these properties ensure monotonic non-decrease toward a stationary point, though they don't guarantee reaching the global optimum.

Closed-form updates and numerical stability

For Gaussian mixture models, maximizing $\sum_Z q(Z) \log p(X, Z \mid \theta)$ produces intuitive weighted updates for each component $k$:

$$N_k = \sum_{i=1}^{n} \gamma_{ik}, \qquad \pi_k^{\text{new}} = \frac{N_k}{n}, \qquad \mu_k^{\text{new}} = \frac{\sum_{i=1}^{n} \gamma_{ik} x_i}{N_k}, \qquad (\sigma_k^2)^{\text{new}} = \frac{\sum_{i=1}^{n} \gamma_{ik} (x_i - \mu_k^{\text{new}})^2}{N_k}.$$

Here $N_k$ represents the effective number of observations assigned to component $k$. The denominator for variance is $N_k$ rather than $N_k - 1$ because this is maximum likelihood estimation, not an unbiased sample variance calculation.

When implementing these updates in software, direct calculation of tiny Gaussian densities leads to floating-point underflow. For points far from a mean, raw evaluation of $\exp(-0.5 (x - \mu)^2 / \sigma^2)$ produces zero. When every component underflows, computing $\gamma_{ik} = 0 / 0$ results in NaN.

The log-sum-exp trick protects calculations by operating in log space:

$$\log \sum_{k=1}^{K} \exp(s_k) = s_{\max} + \log \left( \sum_{k=1}^{K} \exp(s_k - s_{\max}) \right),$$

where $s_{\max} = \max_k s_k$. Subtracting the maximum log score ensures that the largest term evaluates to $\exp(0) = 1.0$, preventing underflow during normalization.

The next snippet fits our six observed latencies [28, 30, 32, 86, 90, 94] starting from provisional means at 25 and 95 ms with initial variances of 900 ms²:

The model settles on component means at 30.0 and 90.0 ms, with variances of 8/3 (2.667 ms²) and 32/3 (10.667 ms²). The slow component's variance is four times larger because observations at [86, 90, 94] spread twice as far from their center as [28, 30, 32] do from theirs.

The variance_floor safeguard prevents an optimization pathology known as covariance collapse. If a component places its mean on a single point with zero variance, its likelihood spikes toward infinity. Imposing a floor keeps the optimization bounded.^{[4]Reference 4Gaussian mixture modelshttps://scikit-learn.org/stable/modules/mixture.html}

The log likelihood stops changing after four iterations. Does that prove the fitted parameters reached the global maximum?

Answer

No. EM guarantees monotonic improvement to a local stationary point, not the global optimum. Another initialization might converge to a distinct mode with higher likelihood.

Diagnosing failure modes: singularities and symmetric traps

Production systems encounter two classic failure modes when fitting probabilistic mixture models.

The first failure mode is covariance collapse (singularities). In standard single-Gaussian maximum likelihood, the likelihood function is strictly bounded. In a Gaussian mixture, however, the likelihood surface contains singularities. If component $k$ assigns its center $\mu_k$ directly to a single observation $x_j$ while its variance shrinks toward zero ($\sigma_k^2 \to 0$), the Gaussian density $\mathcal{N}(x_j \mid x_j, \sigma_k^2) = (2\pi\sigma_k^2)^{-1/2} \to \infty$. The overall log likelihood shoots toward positive infinity. This doesn't reflect a good fit; it's a degenerate singularity where one component collapses into an infinitely sharp spike. In production code, always regularize covariances by adding a diagonal ridge term (such as scikit-learn's reg_covar) or clamping variances with a minimum floor.^{[5]Reference 5GaussianMixturehttps://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html}

The second failure mode is the symmetric initialization trap. If both components start with identical means, variances, and weights (such as $\mu_1 = \mu_2 = 60.0$ ms), every single observation receives identical responsibilities: $\gamma_{i1} = \gamma_{i2} = 0.5$. The subsequent M-step updates both components to the identical pooled sample mean (60.0 ms) and pooled sample variance (878.7 ms²). Because the responsibilities never differentiate, exact EM steps can't break symmetry.

The next snippet verifies this behavior by testing both starting points against an illustrative validation set [29, 31, 88, 92]:

The separated initialization reaches a training log likelihood of -17.694 and a validation score of -2.637 per point. The symmetric initialization remains trapped at the pooled mean 60.0 ms, scoring an inferior -28.943 on training data and -4.822 on validation traffic. Always initialize GMM components with dispersed seeds, such as K-means++ centers or multiple random restarts (n_init=10).

Comparing models with different numbers of components ($K$) requires held-out validation likelihood or information criteria (AIC/BIC). Adding more components always increases training likelihood, but risks overfitting noise or creating degenerate single-point components.

To evaluate model fidelity beyond scalar likelihood scores, draw synthetic samples from the fitted mixture: choose component $k$ with probability $\pi_k$, then draw $x \sim \mathcal{N}(\mu_k, \sigma_k^2)$. Compare histograms of synthetic points against holdout traffic. If the simulated data misses the bimodal separation or generates negative latencies, the Gaussian assumption itself may need revision (for instance, adopting log-normal or Gamma component densities).^{[6]Reference 6Posterior Predictive Checkshttps://mc-stan.org/docs/stan-users-guide/posterior-predictive-checks.html}

Bridge to modern deep learning: Variational Autoencoders as amortized EM

In a Gaussian mixture model, the latent variable $Z$ is a discrete index $k \in \{1, \dots, K\}$. Because $K$ is small, calculating the exact posterior responsibility $\gamma_{ik} = P(Z = k \mid X)$ in the E-step is fast and exact.

Modern generative AI takes this formulation into high dimensions. What happens when:

1. The latent representation $z \in \mathbb{R}^d$ is a continuous, multi-dimensional vector (such as a 128-dimensional embedding representing an image or audio clip)?
2. The generative process $p_\theta(x \mid z)$ is parameterized by a deep neural network (the decoder)?

Now marginal likelihood faces double intractability:

• Marginalizing over continuous latent space requires an intractable high-dimensional integral: $p_\theta(x) = \int p_\theta(x, z)\,dz$.
• The exact posterior $p_\theta(z \mid x) = p_\theta(x, z) / p_\theta(x)$ is impossible to compute analytically, ruling out an exact classical E-step.

Kingma and Welling (2014)^{[7]Reference 7Auto-Encoding Variational Bayes.https://arxiv.org/abs/1312.6114} solved this in the Variational Autoencoder (VAE) framework through amortized variational EM:

Rather than optimizing a separate variational distribution for every data point, an encoder network with weights $\phi$ amortizes inference by predicting posterior parameters $\mu_\phi(x)$ and $\sigma_\phi(x)$ in a single forward pass.

The training objective is the exact same Evidence Lower Bound:

$$\mathcal{L}(\phi, \theta; x) = \mathbb{E}_{q_\phi(z \mid x)} [\log p_\theta(x \mid z)] - D_{\mathrm{KL}}(q_\phi(z \mid x) \parallel p(z)).$$

The first term acts as the M-step reconstruction objective, training the decoder to reconstruct input $x$ from latent code $z$. The second term acts as the E-step regularizer, penalizing the divergence between the approximate posterior $q_\phi(z \mid x)$ and standard Gaussian prior $p(z) = \mathcal{N}(0, I)$.

To backpropagate through stochastic latent variables, Kingma and Welling introduced the reparameterization trick: sample external noise $\epsilon \sim \mathcal{N}(0, I)$ and compute $z = \mu_\phi(x) + \sigma_\phi(x) \odot \epsilon$. Because randomness enters as an independent input, gradients flow smoothly through $\mu_\phi$ and $\sigma_\phi$ to train encoder and decoder end-to-end. The 1977 principle of optimizing marginal likelihood by tightening and climbing a lower bound directly powers modern latent variable architectures, from VAEs to latent diffusion models.