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)
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'?).