Filling up a numpy array using data from a list
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'?).