Gaussian Mixture Models and Expectation-Maximization › Formula

This section derives the mixture density, introduces the latent-variable formulation, and develops the E- and M-step updates.

The derivation starts with the smallest object we need: a single Gaussian bump. In one dimension, with mean $\mu$, positive standard deviation $\sigma>0$ and variance $\sigma^2$, its density is

$$\mathcal{N}(x \mid \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).$$

The fraction in front is not decoration. It is the exact constant that makes the area under the curve equal one, so the result can be treated as a probability density. The exponent measures how far $x$ sits from $\mu$ in units of standard deviations, then squares that distance; the negative sign turns distance into a penalty, so points near the mean get high density and points far away get low density. A single Gaussian therefore encodes one centre and one spread. That is useful, but it cannot describe data with two separated modes.

A mixture removes that limitation by averaging several Gaussian densities together:

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

where:

• $\boldsymbol{\theta}$: the full parameter set collecting all means $\mu_k$, variances $\sigma_k^2$, and mixing weights $\pi_k$
• $\pi_k$: the mixing weight of component $k$, satisfying $\pi_k \ge 0$ and $\sum_{k=1}^K \pi_k = 1$, so the weights behave like probabilities
• $\mu_k$: the mean (centre) of component $k$
• $\sigma_k^2$: the variance of component $k$ (one-dimensional case; the multivariate case replaces it with a covariance matrix $\boldsymbol{\Sigma}_k$ below)
• $\mathcal{N}(x \mid \mu_k, \sigma_k^2)$: the Gaussian density of component $k$ evaluated at $x$

The mixing weights describe prior component membership within the model: before seeing an observation's features, $\pi_k$ is its chance of coming from component $k$. These are mixture probabilities, not a Bayesian prior on the fitted parameters. The sum says that the overall density at $x$ is a weighted average of the component densities at $x$: each component contributes in proportion to how much prior mass it owns.

This weighted sum is also a generative recipe. To produce a point, first choose a component $k$ with probability $\pi_k$, then draw the point from that component's Gaussian. The recipe contains a hidden step, and naming that step is what makes inference possible. Let $z_i \in \{1, \dots, K\}$ record which component produced observation $i$. We never observe $z_i$, which is why it is called a latent variable. Writing the model jointly over the observation and its hidden label gives

$$p(\mathbf{x}_i, z_i = k \mid \boldsymbol{\theta}) = \pi_k \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k),$$

and the law of total probability recovers the marginal density we use:

$$p(\mathbf{x}_i \mid \boldsymbol{\theta}) = \sum_{k=1}^{K} p(\mathbf{x}_i, z_i = k \mid \boldsymbol{\theta}).$$

That identity is the reason the latent variable earns its place. The marginal likelihood contains a sum over hidden states, and that sum is exactly what blocks a direct closed-form maximum. If we could condition on $z_i$, the sum would disappear; the latent variable lets us separate the difficult marginal problem into a sequence of easier complete-data problems.

Bayes' rule then turns the generative model into an inference procedure. Given the current parameters, the posterior probability that observation $i$ belongs to component $k$ is

$$\gamma_{ik} = p(z_i = k \mid \mathbf{x}_i, \boldsymbol{\theta}) = \frac{\pi_k \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}.$$

The denominator is the total density at $\mathbf{x}_i$ contributed by all components, and the numerator is the piece of that total contributed by component $k$. Their ratio is the share of the density that component $k$ explains, which is precisely a probability. These $\gamma_{ik}$ are the soft assignments. They replace the missing labels with graded beliefs instead of forcing each point into one cluster.

Before generalizing, it is worth seeing why the obvious optimization strategy fails. The log-likelihood of $n$ independent observations is

$$\ell(\boldsymbol{\theta}) = \sum_{i=1}^{n} \log\!\left( \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right).$$

We take logs because the likelihood itself is a product of $n$ densities, and products of many small numbers underflow quickly on a computer; logs turn products into sums without moving the maximum, since $\log$ is strictly increasing. The obstacle is the outer logarithm sitting on top of a sum. Differentiating with respect to $\boldsymbol{\mu}_k$ and setting the result to zero produces

$$\sum_{i=1}^{n} \gamma_{ik} \, \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0},$$

because the ratio $\pi_k \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) / \sum_j \pi_j \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)$ is precisely $\gamma_{ik}$. Solving for the mean gives

$$\boldsymbol{\mu}_k = \frac{\sum_{i=1}^{n} \gamma_{ik}\, \mathbf{x}_i}{\sum_{i=1}^{n} \gamma_{ik}}.$$

This looks like a closed-form answer, but it is not: the responsibilities $\gamma_{ik}$ on the right-hand side are themselves functions of the whole parameter set $\boldsymbol{\theta}$, including $\boldsymbol{\mu}_k$. We have a fixed-point equation, not a formula. EM's insight is to freeze the right-hand side at its current value, recompute the left-hand side, and repeat. The E-step supplies the current beliefs about the hidden labels; the M-step turns those beliefs into new parameters.

Now the multivariate case, which is what we will implement. For an observation $\mathbf{x} \in \mathbb{R}^d$ the component density is

$$\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = (2\pi)^{-d/2} \lvert \boldsymbol{\Sigma} \rvert^{-1/2} \exp\!\left(-\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\right),$$

where $\boldsymbol{\mu}$ is a $d$-vector, $\boldsymbol{\Sigma}$ is a $d \times d$ symmetric positive definite covariance matrix, the superscript $\top$ means transpose (turn the column vector into a row vector so the matrix product is defined), $\boldsymbol{\Sigma}^{-1}$ is the matrix inverse, and $\lvert \boldsymbol{\Sigma} \rvert$ is the determinant. The quadratic form $(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ is the squared Mahalanobis distance: it rescales and de-rotates space so that distance is measured in units of the component's own spread. The factor $\lvert \boldsymbol{\Sigma} \rvert^{-1/2}$ does the same job the $1/\sigma$ factor did in one dimension: it divides by the volume the Gaussian occupies, so the density integrates to one.

EM alternates between two problems, each of which is easy on its own. The first problem is inference: given parameters, what do we believe about the hidden labels? The second is estimation: given beliefs about the labels, what parameters are most likely?

Both problems are phrased in terms of the complete-data log-likelihood, the quantity we would maximize if the labels were observed. With $z_{ik}$ a one-hot indicator equal to one when observation $i$ comes from component $k$,

$$\ell_c(\boldsymbol{\theta}) = \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \Big[ \log \pi_k + \log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \Big].$$

where:

• $\ell_c(\boldsymbol{\theta})$: the complete-data log-likelihood, the log-likelihood we would maximize if the latent labels were observed
• $z_{ik}$: a one-hot indicator equal to $1$ when observation $i$ came from component $k$, and $0$ otherwise
• $\pi_k$: the mixing weight of component $k$
• $\mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$: the Gaussian density of component $k$ evaluated at $\mathbf{x}_i$
• $n$: the number of observations, $K$: the number of components

Because the indicators select exactly one term per observation, each component's mean and covariance can be fitted from its assigned observations. The mixing weights still share their sum-to-one constraint and are estimated from the component counts. That is the structure the marginal likelihood destroyed.

Since $z_{ik}$ is unknown, we replace it with its conditional expectation, which is the responsibility: $\mathbb{E}[z_{ik} \mid \mathbf{x}_i, \boldsymbol{\theta}] = \gamma_{ik}$. This produces the auxiliary function

$$Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{t}) = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} \Big[ \log \pi_k + \log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \Big],$$

where $\boldsymbol{\theta}^{t}$ denotes the parameters at iteration $t$. The E-step consists of computing the $\gamma_{ik}$ from the current parameters using Bayes' rule above. Each $\gamma_{ik}$ depends on all components through the denominator. Holding responsibilities fixed separates the mean and covariance estimation by component, while the mixing weights retain their shared normalization.

The M-step maximizes $Q$ over $\boldsymbol{\theta}$. For the means, only the quadratic term depends on $\boldsymbol{\mu}_k$:

$$\log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) = -\tfrac{d}{2}\log(2\pi) - \tfrac{1}{2}\log\lvert \boldsymbol{\Sigma}_k \rvert - \tfrac{1}{2}(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k).$$

Differentiating the quadratic form with respect to the vector $\boldsymbol{\mu}_k$ uses the chain rule on the squared Mahalanobis distance,

$$\frac{\partial}{\partial \boldsymbol{\mu}_k} (\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = -2\,\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_i - \boldsymbol{\mu}_k),$$

where the minus sign appears because $\boldsymbol{\mu}_k$ enters with a negative sign inside the parenthesis. Setting $\partial Q / \partial \boldsymbol{\mu}_k = \mathbf{0}$ gives

$$\sum_{i=1}^{n} \gamma_{ik} \, \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0}.$$

Because $\boldsymbol{\Sigma}_k^{-1}$ is invertible, we may multiply both sides on the left by $\boldsymbol{\Sigma}_k$ to remove it, leaving $\sum_i \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0}$. Splitting the sum and solving,

$$\boldsymbol{\mu}_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik} \, \mathbf{x}_i, \qquad N_k = \sum_{i=1}^{n} \gamma_{ik}.$$

where:

• $N_k = \sum_{i=1}^{n} \gamma_{ik}$: the effective number of observations assigned to component $k$, formed by summing the responsibilities
• $\gamma_{ik}$: the responsibility of component $k$ for observation $i$, its posterior probability of membership
• $\mathbf{x}_i$: the $i$-th observation

The mean update requires $N_k>0$. When a component has zero effective count, this weighted-mean formula is undefined. We can stop with an error, reinitialize or remove the component, or retain its parameters as an inactive zero-weight component and skip its update. Our implementation below stops with an error. The effective count need not be an integer, because responsibilities are fractional: a point with $\gamma_{ik} = 0.3$ contributes $0.3\mathbf{x}_i$ to the weighted feature sum, then division by $N_k$ converts that sum to the mean. This is the soft-count idea in its simplest form.

For the covariance, we need two matrix-calculus facts. First, $\partial \log\lvert \boldsymbol{\Sigma} \rvert / \partial \boldsymbol{\Sigma} = \boldsymbol{\Sigma}^{-1}$ for a symmetric positive definite matrix, which is the matrix analogue of $\partial \log \sigma^2 / \partial \sigma^2 = 1/\sigma^2$. Second, $\partial (\mathbf{a}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{a})/\partial \boldsymbol{\Sigma} = -\boldsymbol{\Sigma}^{-1}\mathbf{a}\mathbf{a}^{\top}\boldsymbol{\Sigma}^{-1}$, which follows from differentiating the inverse. Applying both,

$$\frac{\partial Q}{\partial \boldsymbol{\Sigma}_k} = \sum_{i=1}^{n} \gamma_{ik} \left[ -\tfrac{1}{2}\boldsymbol{\Sigma}_k^{-1} + \tfrac{1}{2}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}\boldsymbol{\Sigma}_k^{-1} \right] = \mathbf{0}.$$

Multiply on the left and right by $\boldsymbol{\Sigma}_k$ (which is legal since it is invertible) and use $N_k = \sum_i \gamma_{ik}$:

$$-N_k \boldsymbol{\Sigma}_k + \sum_{i=1}^{n} \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} = \mathbf{0} \quad\Longrightarrow\quad \boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}.$$

where:

• $\boldsymbol{\Sigma}_k$: the $d \times d$ covariance matrix of component $k$
• $\gamma_{ik}$: the responsibility of component $k$ for observation $i$
• $\mathbf{x}_i - \boldsymbol{\mu}_k$: the centered observation, its deviation from component $k$'s mean
• $N_k$: the effective count for component $k$

For this to be a valid nondegenerate Gaussian covariance we need $N_k>0$ and a positive-definite scatter. This is the weighted scatter matrix of component $k$, with each observation's outer product $(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}$ counted $\gamma_{ik}$ times. It is the same estimator we would use with known labels, except that the counts are fractional.

The mixing weights carry a constraint, $\sum_k \pi_k = 1$, so we cannot differentiate and set to zero naively; we add the constraint times a Lagrange multiplier $\lambda$ and differentiate the augmented expression:

$$\mathcal{L} = \sum_{i=1}^{n}\sum_{k=1}^{K} \gamma_{ik} \log \pi_k + \lambda\left(\sum_{k=1}^{K}\pi_k - 1\right), \qquad \frac{\partial \mathcal{L}}{\partial \pi_k} = \frac{N_k}{\pi_k} + \lambda = 0.$$

The first term's derivative is $N_k/\pi_k$ because $\partial \log \pi_k/\partial \pi_k = 1/\pi_k$, and the Lagrange multiplier $\lambda$ is a constant; note that although $\gamma_{ik}$ itself depends on $\pi_k$ through the Bayes formula, in $Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t)$ the responsibilities are held fixed at their current E-step values, so within this maximization they act as constants on the $\pi_k$ derivative. Rearranging gives $\pi_k = N_k / (-\lambda)$. Summing both sides over $k$ and using $\sum_k \pi_k = 1$ and $\sum_k N_k = \sum_i \sum_k \gamma_{ik} = n$ (each observation's responsibilities sum to one),

$$1 = \sum_{k=1}^{K}\pi_k = \frac{\sum_{k=1}^{K} N_k}{-\lambda} = \frac{n}{-\lambda} \quad\Longrightarrow\quad -\lambda = n \quad\Longrightarrow\quad \pi_k = \frac{N_k}{n}.$$

where:

• $\pi_k$: the mixing weight of component $k$
• $N_k$: the effective count for component $k$
• $n$: the number of observations
• $\lambda$: the Lagrange multiplier enforcing the constraint $\sum_k \pi_k = 1$

So each mixing weight is just the share of the effective count that component $k$ holds.

Collecting everything, one full EM iteration is:

$$\gamma_{ik} = \frac{\pi_k \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} \qquad\text{(E-step)}$$ $$N_k = \sum_i \gamma_{ik}, \qquad \boldsymbol{\mu}_k = \frac{1}{N_k}\sum_i \gamma_{ik}\mathbf{x}_i, \qquad \boldsymbol{\Sigma}_k = \frac{1}{N_k}\sum_i \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}, \qquad \pi_k = \frac{N_k}{n} \qquad\text{(M-step)}$$

The matrix form of the M-step makes the computation obvious. Gather the responsibilities into an $n \times K$ matrix $\mathbf{R}$ whose entry $(i,k)$ is $\gamma_{ik}$, and the data into an $n \times d$ matrix $\mathbf{X}$ whose rows are the observations. Let $\mathbf{1}_n$ be a column vector of ones and $\mathbf{n} = \mathbf{R}^{\top}\mathbf{1}_n$ the column vector of effective counts. Then

$$\mathbf{M} = \operatorname{diag}(\mathbf{n})^{-1} \mathbf{R}^{\top}\mathbf{X}, \qquad \boldsymbol{\Sigma}_k = \frac{1}{n_k} (\mathbf{X} - \mathbf{1}_n \boldsymbol{\mu}_k^{\top})^{\top} \operatorname{diag}(\mathbf{r}_k) (\mathbf{X} - \mathbf{1}_n \boldsymbol{\mu}_k^{\top}),$$

where:

• $\mathbf{R}$: the $n \times K$ responsibility matrix with entry $(i,k) = \gamma_{ik}$
• $\mathbf{X}$: the $n \times d$ data matrix whose rows are the observations
• $\mathbf{1}_n$: a column vector of ones of length $n$
• $\mathbf{n} = \mathbf{R}^{\top}\mathbf{1}_n$: the column vector of effective counts whose $k$-th entry is $N_k$
• $n_k$: the $k$-th entry of $\mathbf{n}$ (same as $N_k$), the effective count of component $k$
• $\mathbf{r}_k$: the $k$-th column of $\mathbf{R}$ (the responsibilities for component $k$)
• $\mathbf{M}$: the $K \times d$ matrix whose $k$-th row is $\boldsymbol{\mu}_k^{\top}$
• $\boldsymbol{\mu}_k$: the mean of component $k$
• $\operatorname{diag}(\cdot)$: turns a vector into a diagonal matrix

Multiplying $\mathbf{R}^{\top}$ by $\mathbf{X}$ sums the rows of $\mathbf{X}$ weighted by the responsibilities (each column accumulates a weighted sum), and dividing by $\mathbf{n}$ turns those weighted sums into weighted averages, giving the $K \times d$ matrix $\mathbf{M}$ whose $k$-th row is $\boldsymbol{\mu}_k^{\top}$. The covariance expression is the same operation applied to the centered data: subtracting $\mathbf{1}_n\boldsymbol{\mu}_k^{\top}$ centers every row, the diagonal matrix reweights the rows, and the product $\mathbf{A}^{\top}\mathbf{W}\mathbf{A}$ accumulates the weighted outer products. In code this is a handful of matrix multiplications rather than $n$ nested loops, which matters because the E-step is repeated until convergence.

The reason EM works at all is that it is coordinate ascent on a lower bound. For a probability distribution $q(z)$ that is positive on the joint model's positive support, multiplying and dividing inside the log by $q(z)$ and applying Jensen's inequality (which says that for a concave function like $\log$, the average of the function is at most the function of the average) gives

$$\log p(\mathbf{x}\mid\boldsymbol{\theta}) = \log \sum_{z} q(z)\frac{p(\mathbf{x}, z \mid \boldsymbol{\theta})}{q(z)} \;\ge\; \sum_{z} q(z) \log \frac{p(\mathbf{x}, z \mid \boldsymbol{\theta})}{q(z)}.$$

The gap between the two sides is the Kullback-Leibler divergence from $q$ to the posterior $p(z \mid \mathbf{x}, \boldsymbol{\theta})$, which is non-negative and equals zero exactly when $q$ is the posterior. The E-step sets $q$ to that posterior, closing the gap so the bound touches the log-likelihood at the current parameters. The M-step then pushes the bound up by maximizing it over $\boldsymbol{\theta}$. Since the log-likelihood is at least the bound, and the bound has increased, the log-likelihood cannot decrease. This proves likelihood nondecrease for the exact bound-increasing updates. It does not prove parameter convergence, a finite maximum, or a global optimum, and it does not apply automatically after an arbitrary covariance modification.

A few properties follow directly from the derivation and shape how the method is used.

• Monotone under exact updates. With valid parameters and an exact E-step followed by a bound-increasing M-step, $\ell(\boldsymbol{\theta}^{t+1}) \ge \ell(\boldsymbol{\theta}^{t})$. The likelihood is non-concave; finite interior limits under suitable regularity are stationary, and they need not be global maxima. Singular boundary behaviour remains possible.
• Unbounded above for unrestricted mixtures with $K\ge2$. A component can collapse onto one observation while another covers the remaining points, driving the sample likelihood to $+\infty$. A single Gaussian whose empirical centered covariance is positive definite does have a finite maximum-likelihood fit; the collapse claim is not a statement about every Gaussian model.
• Label switching. Any permutation of the component indices leaves the likelihood unchanged, so its component labels are arbitrary. Comparing parameter estimates across runs requires matching components first.
• Affine transformation of full-covariance models. For $\mathbf{x}'=A\mathbf{x}+\mathbf{b}$ with invertible $A$, transform each mean to $A\boldsymbol{\mu}_k+\mathbf{b}$ and covariance to $A\boldsymbol{\Sigma}_k A^{\top}$. The density divides by $|\det A|$, so responsibilities remain unchanged and total log-likelihood shifts by $-n\log|\det A|$. This relates corresponding unrestricted fits, not arbitrary numerical runs with fixed regularization. Diagonal models are closed under axis scaling but not general rotations; spherical models permit common uniform scaling and rotations, not arbitrary axis scaling.
• Asymptotic relation to k-means. For fixed observations and centroids, suppose all components have equal fixed weights and common spherical covariance $\sigma^2 \mathbf{I}$. As the externally specified common variance tends to zero, responsibilities at observations with a unique nearest centroid approach nearest-centroid indicators. This varies a prescribed variance across models rather than estimating it in the M-step. For each component with a positive limiting hard-assignment count, its mean update approaches the k-means update. Empty limiting clusters require a separate handling rule. This is the small-variance connection between the mixture and k-means.
• Parameter counts. A full-covariance mixture has $K\big(d + d(d+1)/2\big) + (K-1)$ free parameters: $Kd$ means, $K$ covariance matrices with $d(d+1)/2$ distinct entries each, and $K-1$ free mixing weights. Diagonal covariances cost $2Kd + (K-1)$, spherical ones cost $K(d+1) + (K-1)$, and a tied covariance uses $Kd+d(d+1)/2+(K-1)$ parameters by sharing one covariance across components.