I have two 2D arrays of size N x 3968 and M x 3968, respectively. They represent the N and M feature vectors of two sets of molecules (N and M), each of which consists of 3968 features (floating point numbers). I want to measure the divergence between these two multivariate distributions, which don't follow any norm (they are not multivariate normal or multinomial, etc.). Is there an implementation in Python to compute the Jensen-Shannon divergence or Kullback-Leibler divergence between two such multi-variate distributions?
I have found a function in dit package but it seems to take only 1D distributions: https://dit.readthedocs.io/en/latest/measures/divergences/jensen_shannon_divergence.html
Also scipy has an implementation of Jensen-Shannon distance but it also computes it between 1D distributions: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.jensenshannon.html?highlight=jensenshannon
The latter has an argument that controls along which axis the distances are computed and at the end returns a 1D array of distances. Given my two 2D arrays (N x 3968 and M x 3968), would computing the Jensen-Shannon distance along axis=0 and subsequently averaging the resulting 3968 distances to a single number, give a reasonable estimate of the divergence of the two multivariate distributions? No, because the mean doesn't capture variable (feature) dependences. What about their product?
What about a dimensionality reduction technique to remove nonlinear dependence between features? A technique that estimates the percentage of the variance "pv" that each lower space dimension explains (as PCA does for linear dependence), could facilitate calculation of the weighted mean of the JS distances on the low dimensions, where the weight will be the "pv" of each low dimension.