# How to vectorize/matricize multivariate Gaussian PDF for more efficient computation?

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)}$$

Notation:

• $$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.

• 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? – TheCoolDrop Mar 22 at 7:16
• 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. – charlieshades Mar 22 at 7:26