Jeffries Matusita distance for 14 variables

Question

I wish to perform Jeffries-Matusita distance on 14 spectral bands. Is there anyone who can help with how it is done in R?

Answer 1

As Nick Cox points out in a comment to the question, the Jeffries-Matusita distance should be called the Jeffreys-Matusita distance due to its origin in the work of Harold Jeffreys.

Whatever you choose to call it, distance between two continuous distributions $F$ and $G$ with PDFs $f$ and $g$ is a re-expression of their Bhattacharyya distance:

$$D_{JM}(F,G) = \sqrt{2(1 - \exp(-D_B(F,G))}$$

where

$$D_B(F,G) = \int_\mathbb{R}\sqrt{f(x)g(x)}dx.$$

For nondegenerate multivariate Normal distributions with means $\mu_i$ and covariance matrices $\Sigma_i$ this works out to

$$\eqalign{ &D_B(\mathcal{N}(\mu_F,\Sigma_F), \mathcal{N}(\mu_G,\Sigma_G)) = \\ &\frac{1}{8}D_\Sigma(\mu_F,\mu_G) + \frac{1}{2}\left(\log(|\Sigma|) - \log(|\Sigma_F|)/2 - \log(|\Sigma_G|)/2\right) }$$

where $\Sigma=(\Sigma_F+\Sigma_G)/2$ is the component-wise average of the two covariance matrices, $|\cdot|$ denotes the determinant, and $D_\Sigma$ is the Mahalanobis distance with respect to $\Sigma$,

$$D_\Sigma(\mu_F, \mu_G) = (\mu_F - \mu_G)^\prime \Sigma^{-1}(\mu_F - \mu_G).$$

I presume that the distance is desired between two "signatures." These are estimated multivariate normal distributions of the vectors $(x_1,x_2,\ldots,x_{14})$ associated with distinct groups of pixels in the image. (The components of these vectors are the "bands".)

---

Implementation considerations--all platforms

There are some important implementation details to note, especially for those who wish to generalize this to larger multispectral data, which can have hundreds or even thousands of bands. In these cases the covariance matrices become huge and their determinants can easily overflow double precision floating point numbers. In fact, it's usually not a good idea to compute determinants of large matrices in the first place, but if you have to, an efficient accurate way is to reduce each matrix to a triangular form. Once that is done, the determinant is found as the product of the diagonal elements. Equivalently, the log determinant (which is involved in $D_{JM}$) can be found by summing the logs of the diagonal elements. This avoids over- and underflow.

It also is not necessary to invert $\Sigma$, since all that is needed is the product $\nu=\Sigma^{-1}(\mu_F-\mu_G)$. Find this instead by solving the system of equations

$$\Sigma(\nu)= \mu_F-\mu_G.$$

That will be much faster and a little more precise.

It can be flexible and convenient to implement all three distances mentioned here, because this assists testing and provides additional functionality. Thus, using R to illustrate, we might code something like these functions:

#
# Compute the Mahalanobis distance between two vectors.
#
mahalanobis <- function(m1, m2, sigma) {m <- m1 - m2; m %*% solve(sigma, m)}
#
# Compute the Bhattacharyya distance between two multivariate normal distributions
# given by their means and covariance matrices.
#
bhattacharyya <- function(m1, s1, m2, s2) {
  d <- function(u) determinant(u, logarithm=TRUE)$modulus # Log determinant of matrix u
  s <- (s1 + s2)/2                                        # mean covariance matrix
  mahalanobis(m1, m2, s)/8 + (d(s) - d(s1)/2 - d(s2)/2)/2
}
#
# Re-express the Bhattacharyya distance as the Jeffries-Matusita distance.
#
jeffries.matusita <- function(...) sqrt(2*(1-exp(-bhattacharyya(...))))

That is simple, straightforward, fast, and readily portable to most computing platforms.

---

Illustration in R

An example of the use of these functions in R follows. It generates sample multi-spectral data for a given set of classes, representing them as an array in which each band's values appear in the columns together with a parallel array identifying the class of each row (aka "pixel") with an index $1, 2, \ldots$. It computes the class "signatures" (the estimated multivariate normal parameters). The signature data structure is a list consisting of the mean vector and covariance matrix; the class signatures are collected into one large list. It then applies $D_{JM}$ to all pairs of distinct distance signatures in that list to produce a distance matrix indexed by the classes.

#
# Generate sets of data for d bands aggregated into classes.
#
d <- 14                          # Number of bands
n <- rep(1000, 5)                # Class pixel counts
n.class <- length(n)             # Number of classes
require(MASS)                    # For generating multivariate normals
set.seed(17)                     # Allows reproducible results
Create random mean vectors for the classes
mu <- round(matrix(rnorm(d*n.class, 128, 1), ncol=n.class, byrow=TRUE), 0)
Initialize the data structure {x, classes}
x <- matrix(double(), ncol=d, nrow=0)
classes <- integer()
Generate the random data
for (i in 1:n.class) {
Create a random valid covariance matrix for this class
f <- svd(matrix(rnorm(d^2), ncol=d))
  sigma <- t(f$v) %% diag(rep(10, d)) %% f$v
Generate d-variate normals
x <- rbind(x, mvrnorm(n[i], mu[, i], sigma))
  classes <- c(classes, rep(i, n[i]))
}
#------------------------------------------------------------------------------------#

Given a set of bands (as the columns of x) and a grouping variable in class,
compute the class means and covariance matrices (the "signatures").

classes.stats <- by(x, classes, function(y) list(mean=apply(y, 2, mean), cov=cov(y)))
#------------------------------------------------------------------------------------#

Compute the J-M distances between the classes.

distances <- matrix(0.0, n.class, n.class)
for (i in 2:n.class) {
  m1 <- classes.stats[[i]]$mean; s1 <- classes.stats[[i]]$cov
  for (j in 1:(i-1)) {
    m2 <- classes.stats[[j]]$mean; s2 <- classes.stats[[j]]$cov
    distances[i,j] <- distances[j,i] <- jeffries.matusita(m1,s1,m2,s2)
  }
}
print(distances)