How to estimate the mode using non-parametric methods of a 4-variate random vector drawn from a continuous multivariate distribution? I have a sample of size 10,000 of a 4-variate random vector coming from a (unknown) continuous multivariate distribution.  
How can I estimate the mode of this density using nonparametric methods?
I am currently estimating the corresponding density using the sample and a kernel estimator.  Then I take the value in the grid that maximises this estimator (there are several theoretical results that justifies this procedure).  The packages ks and np are very slow in this case,  particularly in the calculation of the bandwidth matrix. 
Do you have any suggestion?
 A: This is not an answer. I am just posting a code that @PeterEllis suggested but the edit was not accepted for some reason.
library(mvtnorm)
library(misc3d)
library(ks)
library(MASS)


# simulated data, toy example, still hard
dat =  mvrnorm(n=10000, rep(0, 6),diag(6))

#Bandwidth matrix,  takes a looooooooong time
H.scv=Hscv(dat)

# Kernel density estimator, how to specify eval.points?
den = kde(dat, H=H.scv, eval.points=???????)

#Point that maximises the density estimator
ind=which(den$estimate==max(den$estimate), arr.ind=TRUE)

A: A quick glance suggests that this is over my head... but here I go anyway.
In terms of options the docs make it seem like if you had a prior for Hstart things might go a little faster.  Things also look like they might go faster if binned, but I get the sense that binning is forced off for if ncol(dat) > 4.  However at four (as in your text, but not the purported author's code) it looks like you could turn on binning.  The consequences of this are entirely unknown to me.
Manually debugging the code... I saw that most of the time seems to be lost in Hscv (at least for k < 6 (where my computer will finish this month)) during calls in that function to Hpi, gamse.scv, and ultimately nlm. nlm is being passed around as an optimization function.  No nasty loops or easily paralleled apply statements jumped out at me (which didn't/doesn't mean much).  Both nlm and optim (the other choice mentioned in the docs) are .Internal, so I don't think we'll speed them up by much.  
Trying Rprof:  Poking way beyond my level of familiarity, I tried using Rprof in utils to profile the code.  But, I got nowhere.  Most of the crunching that is happening is related to matrix multiplication.  I haven't used it myself yet, but I understand that R 2.14 has 'parallel' built in with parApply, parCapply, and parRapply functions.  If you can find a high enough order non-serial loop or apply statement, perhaps replacing the existing calls to be parallel would help and not just load you down with overhead.
Perhaps someone with a beefer computer than mine, experience in R 2.14, and/or more experience profiling can tell you more.
A: Recently we have suggested a fast consistent mode estimator:

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

Besides, there is a theme on this site devoted to fast mode estimators: Computationally efficient estimation of multivariate mode
