CODEWITHSUNDEEP
pythongaussian
Abrar
Abrar

Reputation: 631

Kullback-Leibler divergence from Gaussian pm,pv to Gaussian qm,qv

I am trying to calculate the Kullback-Leibler divergence from Gaussian#1 to Gaussian#2 I have the mean and the standard deviation for both Gaussians I tried this code from http://www.cs.cmu.edu/~chanwook/MySoftware/rm1_Spk-by-Spk_MLLR/rm1_PNCC_MLLR_1/rm1/python/sphinx/divergence.py

def gau_kl(pm, pv, qm, qv):
    """
    Kullback-Leibler divergence from Gaussian pm,pv to Gaussian qm,qv.
    Also computes KL divergence from a single Gaussian pm,pv to a set
    of Gaussians qm,qv.
    Diagonal covariances are assumed.  Divergence is expressed in nats.
    """
    if (len(qm.shape) == 2):
        axis = 1
    else:
        axis = 0
    # Determinants of diagonal covariances pv, qv
    dpv = pv.prod()
    dqv = qv.prod(axis)
    # Inverse of diagonal covariance qv
    iqv = 1./qv
    # Difference between means pm, qm
    diff = qm - pm
    return (0.5 *
            (numpy.log(dqv / dpv)            # log |\Sigma_q| / |\Sigma_p|
             + (iqv * pv).sum(axis)          # + tr(\Sigma_q^{-1} * \Sigma_p)
             + (diff * iqv * diff).sum(axis) # + (\mu_q-\mu_p)^T\Sigma_q^{-1}(\mu_q-\mu_p)
             - len(pm)))                     # - N

I use the mean and the standard deviation as input, but the last line of the code (len(pm)) cause an error because the mean is one number and I don't understand the len function here.

Note. The two sets(i.e., Gaussians) are not equal that's why I couldn't use the scipy.stats.entropy

Upvotes: 4

Views: 6135

Answers (2)

honeybadger
honeybadger

Reputation: 1575

The following function computes the KL-Divergence between any two multivariate normal distributions (no need for the covariance matrices to be diagonal) (where numpy is imported as np)

def kl_mvn(m0, S0, m1, S1):
    """
    Kullback-Liebler divergence from Gaussian pm,pv to Gaussian qm,qv.
    Also computes KL divergence from a single Gaussian pm,pv to a set
    of Gaussians qm,qv.
    

    From wikipedia
    KL( (m0, S0) || (m1, S1))
         = .5 * ( tr(S1^{-1} S0) + log |S1|/|S0| + 
                  (m1 - m0)^T S1^{-1} (m1 - m0) - N )
    """
    # store inv diag covariance of S1 and diff between means
    N = m0.shape[0]
    iS1 = np.linalg.inv(S1)
    diff = m1 - m0

    # kl is made of three terms
    tr_term   = np.trace(iS1 @ S0)
    det_term  = np.log(np.linalg.det(S1)/np.linalg.det(S0)) #np.sum(np.log(S1)) - np.sum(np.log(S0))
    quad_term = diff.T @ np.linalg.inv(S1) @ diff #np.sum( (diff*diff) * iS1, axis=1)
    #print(tr_term,det_term,quad_term)
    return .5 * (tr_term + det_term + quad_term - N)

Upvotes: 4

mbsariyildiz
mbsariyildiz

Reputation: 43

If you are still interested ...

That function expects diagonal entries of covariance matrix of multivariate Gaussians, not standard deviations as you mention. If your inputs are univariate Gaussians, then both pv and qv are vectors of length 1 for variances of corresponding Gaussians.

Besides, len(pm) corresponds to dimension of mean vectors. It is indeed k in Multivariate normal distributions section here. For univariate Gaussians, k is 1, for bivariate ones k is 2, so on.

Upvotes: 1

Related Questions