Computationally efficient estimation of multivariate mode

Question

Short version: What's the most computationally efficient method of estimating the mode of a multidimensional data set, sampled from a continuous distribution?

Long version: I've got a data set that I need to estimate the mode of. The mode does not coincide with the mean or median. A sample is shown below, this is a 2D example, but an N-D solution would be better:

Currently, my method is

1. Calculate kernel density estimate on a grid equal to the desired resolution of the mode
2. Look for the greatest calculated point

Obviously, this calculates the KDE at a lot of non-plausible points, which is especially bad if there are a lot of data points of high dimensions or I expect good resolution on the mode.

An alternate would be to use a simulated annealing, genetic algorithm, etc to find the global peak in the KDE.

The question is whether there's a smarter method of performing this calculation?

Answer 1

The method that would fit the bill for what you want to do is the mean-shift algorithm. Essentially, mean-shift relies on moving along the direction of the gradient, which is estimated non-parametrically with the "shadow", $K'$ of a given kernel $K$. To wit, if the density $f(x)$ is estimated by $K$, then $\nabla f(x)$ is estimated by $K'$. Details of estimating the gradient of a kernel density are described in Fukunaga and Hostetler (1975), which also happened to introduce the mean-shift algorithm.

A very detailed exposition on the algorithm is also given in this blog entry.

REFERENCES:

• K. Fukunaga and L. Hostetler, "The estimation of the gradient of a density function, with applications in pattern recognition, " IEEE Transactions on Information Theory 21(1), January 1975.

Answer 2

If your main interest is 2-Dimensional problems, I would say that the kernel density estimation is a good choice because it has nice asymptotical properties (note that I am not saying that it is the best). See for example

Parzen, E. (1962). On estimation of a probability density function and mode. Annals of Mathematical Statistics 33: 1065–1076.

de Valpine, P. (2004). Monte Carlo state space likelihoods by weighted posterior kernel density estimation. Journal of the American Statistical Association 99: 523-536.

For higher dimensions (4+) this method is really slow due to the well-known difficulty in estimating the optimal bandwidth matrix, see.

Now, the problem with the command ks in the package KDE is, as you mentioned, that it evaluates the density in a specific grid which can be very limiting. This issue can be solved if you use the package KDE for estimating the bandwidth matrix, using for example Hscv, implement the Kernel density estimator and then optimise this function using the command optim. This is shown below using simulated data and a Gaussian kernel in R.

rm(list=ls())

# Required packages
library(mvtnorm)
library(ks)

# simulated data
set.seed(1)
dat = rmvnorm(1000,c(0,0),diag(2))

# Bandwidth matrix
H.scv=Hlscv(dat)

# [Implementation of the KDE](http://en.wikipedia.org/wiki/Kernel_density_estimation)
H.eig = eigen(H.scv)
H.sqrt = H.eig$vectors %*% diag(sqrt(H.eig$values)) %*% solve(H.eig$vectors)
H = solve(H.sqrt)
dH = det(H.scv)

Gkde = function(par){
return( -log(mean(dmvnorm(t(H%*%t(par-dat)),rep(0,2),diag(2),log=FALSE)/sqrt(dH))))
}

# Optimisation
Max = optim(c(0,0),Gkde)$par
Max

Shape-restricted estimators tend to be faster, for example

Cule, M. L., Samworth, R. J. and Stewart, M. I. (2010). Maximum likelihood estimation of a multi-dimensional log-concave density. Journal Royal Statistical Society B 72: 545–600.

But they are too peaked for this purpose.

The problem in high dimensions is difficult to attack independently of the method used due to the nature of the question itself. For example, the method propopsed in another answer (mean-shift) is nice but it is known that estimating the derivative of a density is even more difficult that estimating the density itself in terms of the errors (I am not criticising this, just pointing out how difficult this problem is). Then you will probably need thousands of observations for accuarately estimating the mode in dimensions higher than $4$ in non-toy problems.

Other methods that you may consider using are: fitting a multivariate finite mixture of normals (or other flexible distributions) or

Abraham, C., Biau, G. and Cadre, B. (2003). Simple estimation of the mode of a multivariate density. The Canadian Journal of Statistics 31: 23–34.

I hope this helps.

Answer 3

Recently we have published a paper suggesting a fast consistent mode estimator.

P.S. Ruzankin and A.V. Logachov (2019). A fast mode estimator in multidimensional space. Statistics & Probability Letters

Our estimator has time complexity $O(dn)$, where $d$ is the dimensionality and $n$ is the number of the observed points. Though our method may be not as precise as the other ones already mentioned here, we write out complete proofs for consistency and strong consistency.

I would also suggest the new minimal variance mode estimators from my recent paper

P.S. Ruzankin (2020). A class of nonparametric mode estimators. Communications in Statistics - Simulation and Computation

Those estimators have time complexity $O(dn^2)$ for $n$ points in ${\mathbb R}^d$. Please see Section 2.3 there. The estimators have precision similar to that of the known algorithms.