Efficient kernel estimation of density over grid of points I was wondering whether it is possible to calculate density estimates for $n \times n$ matrix in $o(n^2)$ computational time. Imagine we have a bivariate random vector $(X,Y)$ with a joint density $f(x,y)$ and realizations denoted by $\{x_i,y_i\}_{i \in 1..n}$, where $n$ is the number of realizations. A standard approach would be for each pair $i,j$ calculate
$$
 \hat{f}(x_i,y_j) = n^{-1}h^{-2} \sum_{k=1}^n K\left(\frac{x_i-x_k}{h},\frac{y_j-y_k}{h}\right),
$$
where $h$ is the bandwidth parameter and $K()$ is the multivariate kernel. For each pair $i,j$ this approach is of $o(n^3)$ computational time. Do you know a faster approach? Is it possible to get $o(n^2)$ computational efficiency?
 A: It seems that it is possible to improve efficiency of the method through matrix multiplication, assuming that $K()$ is a multiplicative kernel. Nevertheless, since multiplication of $n \times n$ matrices is not of $o(n^2)$, the method would be of slightly lower efficiency. 
The following approach would work:
1) calculate two matrices $MX$ and $MY$ where each entry $i,j$ is equal
$$
MX_{i,j} = K\left( \frac{x_i-x_j}{h}\right),\\
MY_{i,j} = K\left( \frac{y_i-y_j}{h}\right).
$$
Both matrices are $n\times n$ and are calculated in $o(n^2)$ computational time. 
2) The estimates $\hat{f}(x_i,y_j)$ are to be taken as $i,j$ elements of matrix $F$ calculated as
$$
F=n^{-1}h^{-2}MX\times MY^T,
$$
where '$\times$' denotes matrix multiplications and superscript $T$ is the matrix transposition. This follows directly from matrix multiplication https://en.wikipedia.org/wiki/Matrix_multiplication. The efficient matrix multiplication algorithm is $o(n^{2.3737})$. Not sure if the above problem can be calculated faster.  
