Improving Nonparametric Density Estimation with Tensor Decompositions

Nonparametric density estimation has been extensively studied with the histogram estimator and KDE being by far the most well known. There do exist, however, alternative methods for density estimation, e.g. the forest density estimator [19] and $k$-nearest neighbor density estimator [20]. The $L_{1},L_{2}$ and $L_{\infty}$ convergence of the histogram and KDE has been studied extensively [12, 8, 33, 14]. The KDE is generally accepted as being the superior density estimator with some mathematical justification [29]. Numerous modifications and extensions of the KDE have been proposed including utilizing variable bandwidth [32], robust KDEs [15, 34, 35], and methods for enforcing support boundary constraints [27]. Finally we mention one recent paper [16] that demonstrated that uniform convergence of a KDE to its population estimate suffered when the intrinsic dimension of the data was lower than the ambient dimension, a phenomenon seemingly at odds with the curse of dimensionality.

For our review of NNTF models we also include a general review of tensor/matrix factorizations since both can be viewed being low-rank models. In particular, for the multi-view model we have

$$\sum_{i=1}^{k}w_{i}f_{i,1}\left(x_{1}\right)\cdots f_{i,d}\left(x_{d}\right)\sim\sum_{i=1}^{k}\lambda_{i}\mathbf{v}_{i,1}\otimes\cdots\otimes\mathbf{v}_{i,d}.$$

(4)

We further note that for histogram estimation, once data has been assigned to bins, finding a good multi-view or Tucker histogram is analogous to estimating a probability tensor with a nonnegative factorization (we show this rigorously in Section 3).

A great deal of work has gone into leveraging low-rank assumptions to improve matrix estimation, particularly in the field of compressed sensing [10, 24]. In compressed sensing one has access to a collection of random linear measurements of an unknown low-rank matrix to be estimated. Fitting a matrix to these measurements with a nuclear norm regularized optimization problem achieves estimation bounds better than those possible without the low-rank assumption. These techniques have been effectively applied to problems such as matrix completion, multivariate regression, and estimating autoregressive models [21, 22]. Unfortunately such techniques are not extensible to histogram estimation because, in the density estimation setting, data are not linearly sampled from the target model. Furthermore how to extend compressed sensing techniques to general tensors is not clear.

General matrix/tensor factorization, including nonnegative matrix/tensor factorizations, have been extensively studied despite their inherent difficulty due to non-convexity. The works [9, 3] present potential theoretical grounds for avoiding the computational difficulties of nonnegative matrix factorization. Some algorithms for finding nonnegative tensor factorizations are mentioned in Section 3. One notable approach to tensor factorization is to assume, in the tensor representation in (4), that $d\geq 3$ and the collections of vectors $\mathbf{v}_{1,j},\ldots,\mathbf{v}_{k,j}$ are linearly independent for all $j$. Under this assumption we are guaranteed that the factorization (4) is unique [1]. In [2] the authors present a method for recovering this factorization efficiently and demonstrate its utility for a variety of tasks. This work was extended in [31] to recover a multi-view KDE satisfying an analogous linear independence assumption. This is the only work of its type of which we are aware. In this work the authors investigate the sample complexity of their estimator but do not demonstrate that their technique has potential for improving rates for nonparametric density estimation in general. Finally we note that nonparametric applications of the Tucker decomposition have been utilized in Bayesian statistics [26]. We are unaware of any literature describing the model we introduce in (3).

All norms in Section 2 are the $\ell^{1},L^{1},$ or total variation norm for vectors/tensors, densities, and measures respectively. Note that these norms are analogous e.g. the $L^{1}$ norm of a probability density function is the same as the total variation norm on the probability measure associated with the density. Let $\mathcal{D}_{d}$ be the set of all densities on $\left[0,1\right)^{d}$. By density we mean probability measures that are absolutely continuous with respect to the $d$-dimensional Lebesgue measure. We define a probability vector or probability tensor to simply mean a vector or tensor whose entires are nonnegative and sum to one. Let $\Delta_{b}$ denote the set of probability vectors in $\mathbb{R}^{b}$and $\mathcal{T}_{d,b}$ the set of probability tensors in $\mathbb{R}^{b^{\times d}}$. Let $\mathcal{T}_{d,b}^{k}$ be the set of tensors which are a convex combination of $k$ separable probability tensors, i.e.

$\displaystyle\mathcal{T}_{d,b}^{k}\triangleq\left\{\sum_{i=1}^{k}w_{i}\prod_{j=1}^{d}p_{i,j}\middle|w\in\Delta_{k},p_{i,j}\in\Delta_{b}\right\}.$

In this paper, the product symbol $\prod$ will always mean the standard outer product, e.g. set product or tensor product. For any natural number $b$ let $[b]=\left\{1,\ldots,b\right\}$. For a multi-index $A\in\left[b\right]^{d}$ we define $\mathbf{e}_{d,b,A}$ as the element of $\mathcal{T}_{d,b}$ where the $(A_{1},\ldots,A_{d})$-th entry is one and is zero elsewhere.

The following is the set of probability tensors constructed via a nonnegative Tucker factorization

$\displaystyle\widetilde{\mathcal{T}}_{d,b}^{k}\triangleq\left\{\sum_{S\in\left[k\right]^{\times d}}W_{S}\prod_{j=1}^{d}p_{i,S_{j}}\middle|W\in\mathcal{T}_{b,k},p_{i,j}\in\Delta_{b}\right\}.$

We will now construct the space of histograms on $\left[0,1\right)^{d}$. We begin with one-dimensional histograms. We define $h_{1,b,i}$ with $i\in\left[b\right]$ to be the one dimensional histogram where all weight is allocated to the $i$th bin. Formally we define this as

$\displaystyle h_{1,b,i}\left(x\right)\triangleq b\mathbbm{1}\left(\frac{i-1}{b}\leq x<\frac{i}{b}\right).$

Note that this is a valid density due to the leading $b$ coefficient. We use these to construct higher-dimensional histograms. For a multi-index $A\ \in\left[b\right]^{d}$, let

$\displaystyle h_{d,b,A}\triangleq\prod_{i=1}^{d}h_{1,b,A_{i}},$

the histogram whose entire density is allocated to the bin indexed by $A$. Finally we define $\Lambda_{d,b,A}$ to be the support of $h_{d,b,A}$, i.e. the “bins” of a histogram estimator,

$\displaystyle\Lambda_{d,b,A}\triangleq\prod_{i=1}^{d}\left[\frac{A_{i}-1}{b},\frac{A_{i}}{b}\right).$

For a sequence of points in $\left[0,1\right)^{d}$, $\mathcal{X}=\left(X_{1},\ldots,X_{n}\right)$, the standard histogram estimator is

$\displaystyle H_{d,b}\left(\mathcal{X}\right)\triangleq\frac{1}{n}\sum_{i=1}^{n}\sum_{A\in[b]^{d}}h_{d,b,A}\mathbbm{1}\left(X_{i}\in\Lambda_{d,b,A}\right).$

Let $\mathcal{H}_{d,b}\triangleq\operatorname{Conv}\left(\left\{h_{d,b,A}\middle|A\in\left[b\right]^{d}\right\}\right)$, the set of all $d$-dimensional histograms with $b$ bins per dimension. Let $\mathcal{H}_{d,b}^{k}$ be the set of histograms with at most $k$ separable components, i.e.

$$\mathcal{H}_{d,b}^{k}\triangleq\left\{\sum_{i=1}^{k}w_{i}\prod_{j=1}^{d}f_{i,j}\middle|w\in\Delta_{k},f_{i,j}\in\mathcal{H}_{1,b}\right\}.$$

(5)

We will refer elements in this space multi-view histograms. Analogously we define the space of Tucker histograms to be

$$\widetilde{\mathcal{H}}_{d,b}^{k}=\left\{\sum_{S\in\left[k\right]^{\times d}}W_{S}\prod_{i=1}^{d}f_{i,S_{i}}\middle|W\in\mathcal{T}_{d,k},f_{i,j}\in\mathcal{H}_{1,b}\right\}.$$

We emphasize that the collections of densities $\mathcal{H}_{d,b}^{k}$ and $\widetilde{\mathcal{H}}_{d,b}^{k}$ the primary objects of interest in this paper as they represent NNTF histograms. The theoretical results we present are concerned with finding good density estimators restricted to these sets as $k$ and $b$ vary.

Note that there exists a $\ell^{1}\to L^{1}$ linear isometry $U_{d,b}:\mathcal{T}_{d,b}\to\mathcal{H}_{d,b}$ with $U_{d,b}$ defined as

$$U_{d,b}(\mathbf{e}_{d,b,A})=h_{d,b,A}.$$

The inverse function, $U^{-1}_{d,b}$, simply transforms a histogram to the tensor representing its bin weights and $U_{d,b}$ performs the reverse transformation. Note that $U_{d,b}$ is also a bijection between $\mathcal{T}_{d,b}^{k}\to\mathcal{H}_{d,b}^{k}$ and $\widetilde{\mathcal{T}}_{d,b}^{k}\to\widetilde{\mathcal{H}}_{d,b}^{k}$. Much of our analysis on histograms will be performed on the space of probability tensors with the analysis being translated to histograms via this operator.

For a set of vectors $\mathcal{V}$ we define $k\operatorname{-mix}\left(\mathcal{V}\right)\triangleq\left\{\sum_{i=1}^{k}w_{i}v_{i}\middle|w\in\Delta_{k},v_{i}\in\mathcal{V}\right\}$, i.e. the set of convex combinations of collections of $k$ vectors from $\mathcal{V}$. We define $N\left(\mathcal{V},\varepsilon\right)$ to be the minimum cardinality for a subset of of $\mathcal{V}$ which $\varepsilon$-covers $\mathcal{V}$ (with closed balls) with respect to the $\left\|\cdot\right\|$ metric. It will be clear from context whether $\left\|\cdot\right\|$ represents the $\ell^{1}$, $L^{1}$ or total variation norm.

For brevity the main text only contains a full proof of Lemma 2.7. The remaining full proofs can be found in the appendix. We include general descriptions of the proof techniques we use for multi-view histogram results. These are similar to the techniques we use for Tucker histograms. Our general proof technique is to find good covers of spaces of densities, i.e. $\mathcal{H}_{d,b}^{k}$ and $\widetilde{\mathcal{H}}_{d,b}^{k}$, and then apply an existing algorithm for selecting a good estimator from finite collections of densities given data. We begin by establishing a covering number bound on the space of probability vectors via an adaptation of a standard result presented in [8].

We can extend this to a covering number for the space of separable probability tensors.

Now we establish the following lemma for covering numbers of mixtures of densities.

By combining Lemma 2.3 with Lemma 2.2 we arrive at covering numbers for the space of multi-view probability tensors.

Through application of the $U_{d,b}$ operator we now have a characterization of the complexity of the space $\mathcal{H}_{d,b}^{k}$.

The following are analogous results for Tucker histograms.

The following lemma from [4] provides us with a way to choose good estimators from finite collections of densities. It can be proven by applying a Chernoff bound to [8], Theorem 6.3.

We present the following asymptotic version of the previous lemma and include its full proof. We highlight the use of finding sufficiently slow rates on parameters in order to establish asymptotic results, a technique which we will use in later proofs.

We can now state the central results of this paper. The following theorem states that one can control the estimation error of multi-view histograms with $k$ components and $b$ bins per dimension so long as $n\sim bk$ (omitting logarithmic factors). Recall that the standard histogram estimator requires $n\sim b^{d}$, so we have removed the exponential dependence of bin rate on dimensionality. Here and elsewhere the $\sim$ symbol is not a precise mathematical statement but rather describes that the two values should be of the same order. In the following $b$ and $k$ are functions of $n$; the space of histograms which we are fitting changes as we acquire more data.

The sample complexity for the multi-view histogram is perhaps more accurately approximated as being on the order of $dbk$ however the $d$ disappears in the asymptotic analysis.

The following theorem states that we can control the error of Tucker histogram estimates so long as $n\sim bk+k^{d}$ (omitting logarithmic factors).

Allowing $b$ to grow as aggressively as possible we achieve consistent estimation, using either the multi-view or Tucker histograms, so long as $n\sim b\log b$ regardless of dimensionality.

Replacing $b:=1/h$ allows us to arrive at the rates mentioned in Section 1. The following result shows that the bias of the estimators in Theorems 2.1 and 2.2 go to zero for all densities. Thus these estimators are universally consistent even when the NNTF model assumption is not satisfied.

A straightforward consequence of this is that the Tucker histogram bias also goes to zero.

Finally we have that rate on $bk$ in Theorem 2.1 cannot be made significantly faster.

Likewise the rate on $bk+k^{d}$ can also not be significantly improved in Theorem 2.2.

Naturally real world data likely never exactly satisfies the NNTF model assumption. Our results are meant to highlight a trade-off between model assumptions of smoothness (low $b$) and simple dependence between features (low $k$). Here we will explore this trade-off for multi-view histogram. Letting $k=b^{d}$ gives $\mathcal{H}_{d,b}^{k}=\mathcal{H}_{d,b}$ since we can allocate one component to each bin. Using the estimator from Theorem 2.1 gives a sample complexity of approximately $b^{d+1}$. Thus setting $k=b^{d}$ in the multi-view histogram gives us something which behaves similarly to the standard histogram estimator with a similar sample complexity. On the other hand setting $k=1$ gives a naive Bayes model with a sample complexity of approximately $b$. The Tucker histogram can be similarly analyzed with $k=b$ corresponding to the standard histogram. Thus we have a span of $k$ yielding different estimators with maximal $k$ corresponding to the standard histogram and minimal $k$ corresponding to a naive Bayes assumption. We observe in Section 3 that this trade-off is useful in practice: we virtually never want $k$ to be maximized.

Our experiments were performed on the Scikit-learn “toy” datasets MNIST and Diabetes [23], with labels removed. We used the expression inside the minimization in (10) to evaluate performance. Our experiments considered estimating histograms in $d=2,3,4,5$ dimensional space. We consider two forms of dimensionality reduction. First we consider projecting the dataset onto its top $d$ principle components. As an alternative we consider projecting our dataset onto a random subspace of dimension $d$. We have constructed our random subspace dimensionality reduction so that each additional dimension adds a new index without affecting the others. For each dataset we randomly select an orthonormal basis that remains unchanged for all experiments $v_{1},v_{2},\ldots$. To transform a point $X$ to dimension $d$ we perform the following transform

$\displaystyle X_{\text{reduced dim.}}=\left[v_{1}\cdots v_{d}\right]^{T}X.$

We consider both transforms since PCA may select dimensions where the features tend to be independent, as is the case when the distribution is a multivariate Gaussian. After dimensionality reduction we scale and translate the data to fit in the unit cube.

With our preprocessed dataset, each experiment consisted of randomly selecting 200 samples for training and using the rest to evaluate performance. For the estimators we tested all combinations using 1 to $b_{\text{max}}$ bins per dimension and $k$ from 1 to $k_{\text{max}}$. As $d$ increased the best cross validated $b$ and $k$ value decreased, so we reduced $b_{\text{max}}$ and $k_{\text{max}}$ for larger $d$ to reduce computational time, while still leaving a sizable gap between the best cross validated $b$ and $k$ across all runs of all experiment. For $d=2,3$ we have $b_{\text{max}}=15$ and $k_{\text{max}}=10$; for $d=4$ we have $b_{\text{max}}=12$ and $k_{\text{max}}=8$; for $d=5$ we have $b_{\text{max}}=8$ and $k_{\text{max}}=6$. For parameter fitting we used random subset cross validation repeated 80 times using 40 of the 200 samples to cross validate. Performing 80 folds of cross validation was necessary because of the variance of the histogram estimated loss. This high variance is likely due to the noncontinuous nature of the histogram estimator itself and the noncontinuity of the histogram as a function of the data, i.e. slightly moving one training sample can potentially the change histogram bin in which it lies. Each experiment was run 32 times and we report the mean and standard deviations of estimator performance as well as the best parameters found from cross validation. We additionally apply the Wilcoxon signed rank test to the 32 pairs of performance results to statistically determine if the mean performance between the standard histogram and our algorithm are different and report the corresponding $p$-value. Our results are in Table 1 where the Tucker histogram dominates. We also observe that the Tucker histogram can estimate more bins per dimension than the standard histogram and is able to estimate more bins per dimension when the number of components of components is lower. This corroborates the number of components versus the bins per dimension trade-off from Section 2.4.