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.
