Nonparametric estimation of a multidimensional probability density

A sample of size $n$ from a $k$-dimensional absolutely continuous distribution being available, the function
$$ f_n(x_1,\dots,x_k)=\frac1n\sum_{i=1}^n\prod_{l=1}^k\frac1{h_l(n)}K_l\biggl(\frac{x_l-x_l^{(i)}}{h_l(n)}\biggr) $$
is taken as a density function estimator, where $K_l(y)$'s are given real-valued functions symmetric with respect to $y=0$ and having bounded moments. $f_n(x_1,\dots,x_k)$ is shown to be asymptotically unbiased and consistent estimate of the probability density at each point $(x_1,\dots,x_k)$ provided that $\lim\limits_{n\to\infty}h_l(n)=0$, $\lim\limits_{n\to\infty}\prod_{l=1}^kh_l(n)\to\infty$. Optimal functions $K_l(y)$ are found which reduce the asymptotic relative total mean-square error to the minimum.