Context: I was recently implementing (in Python) the Expectation-Maximization (EM) algorithm for Gaussian mixture models, and part of that process involves computing the Gaussian PDF for various points $x_i$ with various parameters. I implemented it somewhat naively, which was fast enough for what I needed, but I feel like it could be computed faster if it was written in a nice, clean matrix form. However, I don't know how you'd approach getting it in the form you'd need.

The multivariate Gaussian PDF I need is in the form:

$$ p_k(x_i \mid \mu_k, \Sigma_k) = \frac{1}{(2\pi)^\frac{d}{2} \vert \Sigma_k \vert^{\frac{1}{2}}} e^{-\frac{1}{2} (x_i - \mu_i)^T \Sigma_k^{-1} (x_i - \mu_k)} $$


  • $x_i$ is a $d$-dimensional vector of input data, where $1 \le i \le N$ ($N$ is total number of data points)
  • $k \in {1, 2, ..., K}$ is another parameter over which I have to calculate $p(x_i)$ for each value of $k$
  • $\mu_k$ is $d$-dimensional set of means for the data
  • $\Sigma_k$ is a $d \times d$ covariance matrix for the data

My main points of confusion:

The equation as I typed above gives a scalar probability for a single $i$ and a single $k$. I need to calculate it for all $i$, $k$, so in the end I'd want an $N \times K$ matrix of values.

The main portion of the exponential part, $(x_i - \mu_i)^T \Sigma_k^{-1} (x_i - \mu_k)$, effectively takes a $1 \times d$ vector, multiplies by a $d \times d$ matrix, and then multiplies again by a $d \times 1$ vector. With a single data point we get a $1 \times 1$ scalar as a result, but with more than one data point wouldn't we get an $N \times N$ matrix?

Can anyone help me figure out a more elegant way to write this? Again, what I'd ultimately like to end up with is some one-liner to get an $N \times K$ matrix that I supposed you'd write as: $p(X \mid \mu, \Sigma)$ where $X$ is an $N \times d$ matrix, $\mu$ is a $k \times d$ matrix, and $\Sigma$ is a $k \times d \times d$ matrix.

  • $\begingroup$ Here is how I would go about it. First use Numpy. Second the exponential part can be vectorized by computing the exponent first,and not inputing one by one value, but rather whole vectors/matrices, then you can exponentiate the obtained matrix elementwise ( which is also vectorized in numpy) and then you can multiply it by a matrix. But this is more of a programming than math question. Why not consider asking at stackoverflow? $\endgroup$ – TheCoolDrop Mar 22 at 7:16
  • $\begingroup$ In my implementation I am using numpy, but to compute the PDF I'm also using the multivariate_normal package from scipy.stats, and if I try to just give it everything at one time it complains about dimensions mismatching. I also don't understand how the exponent part can be computed element-wise given the $N \times d$ input, the $k \times d$ means vector, and the $k \times d \times d$ covariance matrix. Edit: Also I did ask on SO and nobody responded, but also I'm more curious about how you'd represent the equation I'm looking for mathematically. $\endgroup$ – charlieshades Mar 22 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.