Computationally Efficient Exact Calculation of Kernel Density Derivatives

Abstract

Machine learning research related to the derivatives of the kernel density estimator has received limited attention compared to the density estimator itself. This is despite of the general consensus that most of the important features of a data distribution, such as modes, curvature or even cluster structure, are characterized by its derivatives. In this paper we present a computationally efficient algorithm to calculate kernel density estimates and their derivatives for linearly separable kernels, with significant savings especially for high dimensional data and higher order derivatives. It significantly reduces the number of operations (multiplications and derivative evaluations) to calculate the estimates, while keeping results exact (i.e. no approximations are involved). The main idea is that the calculation of multivariate separable kernels and their derivatives, such as the gradient vector and the Hessian matrix involves significant number of redundant operations that can be eliminated using the chain rule. A tree-based algorithm that calculates exact kernel density estimate and derivatives in the most efficient fashion is presented with the particular focus being on optimizing kernel evaluations for individual data pairs. In contrast, most approaches in the literature resort to approximations of functions or downsampling. Overall computational savings of the presented method could be further increased by incorporating such approximations, which aim to reduce the number of pairs of data considered. The theoretical computational complexity of the tree-based and direct methods that perform all multiplications are compared. In experimental results, calculating separable kernels and their derivatives is considered, as well as a measure that evaluates how close a point is to the principal curve of a density, which employs first and second derivatives. These results indicate considerable improvement in computational complexity, hence time over the direct approach.

Appendices

Appendix A:

Proof of fact 1

Subject to the constraint $\sum _{l} \left | v_{lm} \right | \leq p$ on the growth of T, the number of nodes in the l-th layer is equal to $\binom {l+p}{p}$. The proof of this is based on mathematical induction for all l=1,…,d. Let the number of nodes in each layer be n _l .

Base case

When l=1, the first dimension contains all the derivative orders j=0,…,p which makes n ₁=p+1. The nodes with index j at the first layer are expanded by p+1−j child nodes in the second layer, and therefore $n_{2}=(p+1)+p+\dots +2+1=\binom {2+p}{2}$ which supports the statement for l=2.

Induction step

Let l ^′ be given as the layer number, and suppose the hypothesis is true for l=l ^′. Using Pascal’s rule about binomial coefficients, n _l ^′+1 can be written as

$$ \binom{l'+1+p}{p}=\binom{l'+p}{p}+\binom{l'+p}{p-1} $$

(19)

where the left term is the number of nodes at the layer l ^′+1. Looking back at the tree structure and Figure 1 as an example, it can be observed that the tree has a self-repeating pattern. At the new layer l ^′+1, the nodes are divided into two parts. The upper part is originally expanded from v ₁₁, and it exactly repeats the nodes of previous layer. The reason is that the starting node of this branch is the node containing j=0 at the first layer, which does not change the maximum derivative order constraint. As a result, the branch can be shifted one layer back to the root, and therefore this recursive nature propagates to the whole tree. In (19), this branch is represented by the term $\binom {l'+p}{p}$ which is equal to n _l ^′, according to the inductionhypothesis.

The lower part of the tree is expanded from all other nodes of the first layer expect v ₁₁. One can observe that these branches are identical to those of a tree with the same number of layers but with the maximum derivative order of p−1. Again the reasoning is based on the maximum derivative constraint. Since at the first layer the node value is nonzero, the nodes have one less derivative order for the valid expansion. Based on the induction hypothesis, the number of nodes in the layer l ^′+1 subject to the constraint $\sum _{l} \left | v_{lm} \right | \leq (p-1)$, will be $\binom {(l'+1)+(p-1)}{p-1}$ or $\binom {l'+p}{p-1}$, and equal to the second term of the the right side of (19). Thus the induction hypothesis is proved for l=l ^′+1, and the statement is generally true.

Appendix B:

In the case of multivariate separable Gaussian kernels, the exponent has the Mahalanobis distance form. Therefore, decomposing the inverse of covariance matrix (precision matrix) gives the matrix that scales and orients the input data relative to the desired kernel bandwidth.

$$ \mathbf{\Sigma}_{i}^{-1} = \mathbf{S}_{i}^{T} \mathbf{S}_{i}. $$

(20)

In order to extract S _i , the eigendecomposition of C _i is utilized

$$ \mathbf{\Sigma}_{i} = \mathbf{Q}_{i} \mathbf{\Lambda}_{i} \mathbf{Q}_{i}^{T}, $$

(21)

where the eigenvectors and eigenvalues of C _i are the columns and diagonal elements of orthogonal Q _i and diagonal Λ _i , respectively. From (20) and (21), the scale matrix S _i is written as follows

$$ \mathbf{S}_{i} = \mathbf{\Lambda}_{i}^{-1/2} \mathbf{Q}_{i}^{T}. $$

(22)

Using this formation, S _i is the bandwidth matrix affiliated with X _i . Similarly, a Cholesky factorization of $\mathbf {C}_{i} =\mathbf {R}^{T}_{i}\mathbf {R}_{i}$ would yield S _i =R _i as an option.