Solved – FFT for speeding up kernel density estimation

Question

One of the method for speeding up Kernel Density Estimation calculations is based on using FFT. This is implemented for example in kde{ks} R function. But unfortunately, a serious problem exists while trying to use this kde{ks} for unconstrained (full) bandwidth matrices. Below I try to explain it.

First, load a sample dataset and set some inputs:

library(ks)
library(mvtnorm)
data(unicef)
x <- cbind(unicef[,1], unicef[,2])
# grid sizes
gs1 <- 40
gs2 <- 40
# for nlevels in contour
nl <- 16
par(mfrow = c(2,2), pty="s")

Now execute kde with binned=F and plot a contour.
H is the unconstrained bandwidth matrix. The density seems quite good (Fig. 1).

H <- Hpi(x)
fhat1 <- kde(x=x, H=H, supp=99, gridsize=c(gs1,gs2), xmin=c(min(x[,1]),min(x[,2])), xmax=c(max(x[,1]),max(x[,2])), binned=F)
contour(fhat1$eval.points[[1]], fhat1$eval.points[[2]], fhat1$estimate, main="H is unconstrained, binned=F", nlevels=nl, xlab="Fig. 1")
points(fhat1$x, fhat1$y, col='red')

Now I do the same, but now with binned=T.
I get an alert: Binned estimation for non-diagonal bandwidth matrix H can be inaccurate.
fhat2$H shows that the unconstrained bandwidth was used as an input. The density seems rather strange (Fig. 2).

fhat2 <- kde(x=x, H=H, supp=99, bgridsize=c(gs1,gs2), xmin=c(min(x[,1]),min(x[,2])), xmax=c(max(x[,1]),max(x[,2])), binned=T)
contour(fhat2$eval.points[[1]], fhat2$eval.points[[2]], fhat2$estimate, main="H is unconstrained, binned=T", nlevels=nl, xlab="Fig. 2")

Now I change manually H into a diagonal form and execute kde again.
Fig. 3 is almost identical (or at least similar) to Fig. 2. This suggests that kde uses (internally) a form of the diagonal H, instead of the unconstrained H.

H <- Hpi.diag(x)
fhat3 <- kde(x=x, H=H, supp=99, bgridsize=c(gs1,gs2), xmin=c(min(x[,1]),min(x[,2])), xmax=c(max(x[,1]),max(x[,2])), binned=T)
contour(fhat3$eval.points[[1]], fhat3$eval.points[[2]], fhat3$estimate, main="H diagonal, binned=T", nlevels = 20, xlab="Fig. 3")

Now I become suspicious that turning on the binning makes some difficulties while using kde. After a short inspection of ks source codes I realized that FFT was used to speed up the calculations and this FFT procedure is not "ready" for supporting unconstrained H bandwidths.
Is it true? To check my suspicions I have create a simple code for calculating fhat WITH binning and WITHOUT using FFT (that is with radially symmetric kernels). Below is the code (not optimized for speed, literal implementation of the definition in fact):

bkde.no.fft.radial <-  function(x, bandwidth, gridsize) {
  range.x <- list(c(min(x[,1]),max(x[,1])),c(min(x[,2]),max(x[,2])))
  n <- nrow(x)
  M <- gridsize
  H <- bandwidth
  d <- ncol(x)

  a <- c(range.x[[1]][1], range.x[[2]][2])
  b <- c(range.x[[1]][3], range.x[[2]][4])

  binned <- binning(x=x, bgridsize=c(M), xmin=c(a[1],a[2]), xmax=c(b[1],b[2]))
  gcounts <- binned$counts
  gpoints1 <- binned$eval.points[[1]]
  gpoints2 <- binned$eval.points[[2]]

  temp <- 0
  row <- 0
  col <- 0
  result <- matrix(nrow=M[1], ncol=M[2])
  for(i1 in 1:M[1]) {
    row <- row + 1
    for(i2 in 1:M[2]) {
      # for observing the progress
      cat(i1, "/", M[1], "---", i2, "/", M[2], "\r")
      flush.console()
      col <- col + 1
      for(j1 in 1:M[1]) {
        for(j2 in 1:M[2]) {
          if (gcounts[j1,j2] != 0) {
            temp <- temp + dmvnorm(c(gpoints1[i1],gpoints2[i2]) - c(gpoints1[j1],gpoints2[j2]), sigma=H) * gcounts[j1,j2]
          }
        }
      }
      result[row, col] <- temp / n
      temp <- 0
    }
    col <- 0
  }
  list(x1=gpoints1, x2=gpoints2, estimate=result)
}

Execute the above function (works very slow):

H <- Hpi(x)
fhat4 <- bkde.no.fft.radial(x=x, bandwidth=H, gridsize=c(gs1,gs2))
contour(fhat4$x1, fhat4$x2, fhat4$estimate, main="H is unconstrained, binned=T, radially symetric kernels", nlevels = 20, xlab="Fig. 4")

It can be observed that Fig. 1 is identical with Fig. 4 and it proves that the unconstrained H was used correctly. This experiment shows that kde with binned=T and the unconstrained bandwidth H used CAN NOT make use of radially symmetric kernels if this is implemented with FFT!

And now my the most interesting and fascinating question:
Is it possible to utilize FFT for implementing fast kde with the unconstrained bandwidth H support?

Below I enclose the image generated after executing all the code fragments given above.
![figures][5]

The problem described above seems similar to the one given here, as in both posts FFT-like-nature seems to be responsible for "noise effects".

Answer 1

Recently, the problem was successfully solved. Please, see https://arxiv.org/abs/1508.02766 for details. Similar problem, relating the unconstrained bandwidth selection, was presented here https://arxiv.org/abs/1511.07482