Efficient metric for to measure whether two clouds of points overlap?

Viewed 170

Assume I have two sets of points A and B. A is a matrix of size N-by-D and B is a matrix of size M-by-D. Both sets have the same dimensions D but may be composed of different numbers of samples N,M > 1 (assume that number of samples is sufficiently large).

I am looking for an efficient metric to determine to what degree set A overlaps with set B. This metric should have the following properties:

  1. It should return a large value when set A is contained in set B
  2. It should return an intermediate value when set A partially overlaps with set B
  3. It should return a small value when set A and set B do not overlap

I have thought of some ways to achieve this, but none of them quite make the cut:

  • Determine the convex hull of B, then calculate what percentage of A lies within this convex hull. This is more or less reliable (assuming B is sufficiently convex), but calculating the convex hull becomes prohibitively expensive for large D.
  • Estimate the mean and covariance of A and B and calculate a Kullback-Leibler divergence between the two resulting multivariate Gaussian distributions. This is reasonably efficient, but does not distinguish the case when A is completely embedded in B but has significantly lower spread, and when A and B have similar spread but only overlap partially.

Do you have any other ideas on how I could tackle this issue? Below I have provided an example code illustrating the problem with using the Kullback-Leibler divergence:

import numpy as np
import scipy.stats

N   = 1000
M   = 1000
D   = 3

def metric_KLD(A,B):
    
    # Get the dimension of A and B
    D       = A.shape[-1]
    
    # Estimate mean and cov of A
    A_mean  = np.mean(A,axis=0)
    A_cov   = np.cov(A.T)

    # Estimate mean and cov of B
    B_mean  = np.mean(B,axis=0)
    B_cov   = np.cov(B.T)
    
    # Calculate the KLD
    KLD     = 0.5*(np.log(np.linalg.det(B_cov)/np.linalg.det(A_cov)) - \
        D   + np.trace(np.dot(np.linalg.inv(B_cov),A_cov)) + \
        np.linalg.multi_dot((
            (B_mean - A_mean)[np.newaxis,:],
            np.linalg.inv(B_cov),
            (B_mean - A_mean)[:,np.newaxis])))
        
    return KLD



# Case 1: Both distributions overlap perfectly
A1  = scipy.stats.multivariate_normal.rvs(
    mean    = np.zeros(D),
    cov     = np.identity(D)*10**2,
    size    = N)
B1  = scipy.stats.multivariate_normal.rvs(
    mean    = np.zeros(D),
    cov     = np.identity(D)*10**2,
    size    = M)

print('KLD case 1: '+str(metric_KLD(A1,B1)))



# Case 2: Both distributions overlap partially
A2  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([0,0,0]),
    cov     = np.identity(D)*10**2,
    size    = N)
B2  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([10,0,0]),
    cov     = np.identity(D)*10**2,
    size    = M)

print('KLD case 2: '+str(metric_KLD(A2,B2)))


# Case 3: Both distributions don't overlap at all
A3  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([0,0,0]),
    cov     = np.identity(D)*10**2,
    size    = N)
B3  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([30,30,30]),
    cov     = np.identity(D)*10**2,
    size    = M)

print('KLD case 3: '+str(metric_KLD(A3,B3)))


# Case 4: A is included in B, but A has significantly smaller spread
A4  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([0,0,0]),
    cov     = np.identity(D)*1**2,
    size    = N)
B4  = scipy.stats.multivariate_normal.rvs(
    mean    = np.asarray([0,0,0]),
    cov     = np.identity(D)*10**2,
    size    = M)

# This is problematic: Should be better than case 2 but isn't
print('KLD case 4: '+str(metric_KLD(A4,B4)))
1 Answers

Another idea is to compare the standard deviation within each set to the overall standard deviation in the pooled set. The more the sets overlap, the more similar these numbers will be.

To handle the case of one set with a small spread within the other one, you could just use the larger individual standard deviation for comparison with the pooled one.

This could be done for each dimension, and the final metric would be the average over all dimensions:

def metric_std(A, B):
        
    # standard deviations 
    A_std  = np.std(A, axis=0)
    B_std  = np.std(B, axis=0)
    pooled_std = np.std(np.vstack([A, B]), axis=0)
    
    # metric: ratio of higher individual to pooled std    
    std_ratio = np.max([A_std, B_std], axis=0) / pooled_std
        
    return np.mean(std_ratio)


print('case 1: '+str(metric_std(A1, B1)))
print('case 2: '+str(metric_std(A2, B2)))
print('case 3: '+str(metric_std(A3, B3)))
print('case 4: '+str(metric_std(A4, B4)))
case 1: 1.012685295955191
case 2: 0.973728554241719
case 3: 0.5569612281519413
case 4: 1.4071848466046386
Related