# Filling up a numpy array using data from a list

I am looking to optimize the nested for-loops shown in the code below.

I have a list, mixtures, that contains points that represent 1 Gaussian Mixture (GM) on each line. My goal is to run a function kl_divergence between all combinations of GM's (with i≠j) and store the results in an array PxP.

Note: The following mixtures list is a simplified version of the actual list. In my original code, I am working with a lot more data.

import numpy as np
from sklearn.mixture import GaussianMixture

def kl_divergence(gmm_i, gmm_j, n_samples=10**5):
X = gmm_i.sample(n_samples)
log_p_i = gmm_i.score_samples(X[0])
log_p_j = gmm_j.score_samples(X[0])
return log_p_i.mean() - log_p_j.mean()

mixtures = [[[1, 3, 4, 7], [3, 5, 9, 2], [4, 3, 6, 1], [4, 3, 6, 3]],
[[3, 4, 8, 2], [3, 6, 3, 7], [2, 6, 8, 4]],
[[4, 8, 9, 3], [2, 6, 5, 8], [2, 5, 3, 5]],
[[4, 9, 0, 2], [2, 4, 8, 3], [9, 8, 2, 3], [2, 6, 8, 3]]]
#Note: Each line has a different number of points

P_comp = len(mixtures)
PxP = np.zeros((P_comp, P_comp), dtype=int)

for row_i, GM_i in enumerate(mixtures):    #Loop through mixtures using i
gmm_i = GaussianMixture(n_components=1).fit(GM_i)
for row_j, GM_j in enumerate(mixtures): #Loop through mixtures using j
if GM_i is not GM_j: #Skip if same mixture
gmm_j = GaussianMixture(n_components=1).fit(GM_j)
PxP[row_i, row_j] = kl_divergence(gmm_i, gmm_j) #Store result at index (row_i, row_j)

print("PxP\n", PxP)

• Welcome to CodeReview! Looks like you instantiate a GaussianMixture over and again for any given element of mixtures - is that a cheap operation? Did you try and find out which part of the consumes most time with a lot more data? – greybeard Jul 22 '19 at 3:43
• Are you working with so much data that space might start to become a problem, or "just" so much data that performance might be significant? – Ninetails Jul 22 '19 at 7:35
• @greybeard Here are some performance figures: kl_divergence(gmm_i, gmm_j) runs in 45ms on average. Preparing a GM gmm_i or gmm_j takes anywhere between 10ms and 1000ms, and depends on how many points are in row_i or row_j. – Honey Gourami Jul 23 '19 at 1:02
• @Ninetails Memory hasn't been an issue yet, so it's purely a performance issue. – Honey Gourami Jul 23 '19 at 1:09

The code to establish PxP should be in a function of its own, featuring a doc string and comments where necessary (why not use GaussianMixture(covariance_type ='spherical')?)
Same goes for kl_divergence(): when(/why?) is n_samples=10**5 appropriate?
From eyeballing your code, it would seem that you instantiate a GaussianMixture() not just once for each mixture, but (len(mixtures))²? times.
Try to turn mixtures into GaussianMixture upfront: mixtures = [GaussianMixture().fit(mix) for mix in mixtures] (n_components=1 is default - do you need covariance_type='full'?).
• The improvement I got in computation time after I turned mixtures into GaussianMixture upfront is unbelievable. I was expecting a 1.5-2x boost in performance but it turned out to be much more: 30x. Regarding the doc string and comments, I believe you are right. It will help explain/clarify the parameters I'm using. – Honey Gourami Jul 25 '19 at 4:12