Beyond Smoothness: Incorporating Low-Rank Analysis into Nonparametric Density Estimation
Year: 2022

Abstract:
The construction and theoretical analysis of the most popular universally consistent nonparametric density estimators hinge on one functional property: smoothness. In this paper we investigate the theoretical implications of incorporating a multiview latent variable model, a type of low-rank model, into nonparametric density estimation. To do this we perform extensive analysis on histogram-style estimators that integrate a multi-view model. Our analysis culminates in showing that there exists a universally consistent histogram-style estimator that converges to any multi-view model with a finite number of Lipschitz continuous components at a rate of Õ(1/ 3 √ n) in L error. In contrast, the standard histogram estimator can converge at a rate slower than 1/ d √ n on the same class of densities. We also introduce a new nonparametric latent variable model based on the Tucker decomposition. A rudimentary implementation of our estimators experimentally demonstrates a considerable performance improvement over the standard histogram estimator. We also provide a thorough analysis of the sample complexity of our Tucker decomposition-based model and a variety of other results. Thus, our paper provides solid theoretical foundations for extending low-rank techniques to the nonparametric setting.

text: √ n) in L1 error. In contrast, the standard histogram estimator
can converge at a rate slower than 1/ d √ n on the same class of densities. We also introduce a new nonparametric latent variable model based on the Tucker decomposition. A rudimentary implementation of our estimators experimentally demonstrates a considerable performance improvement over the standard histogram estimator. We also provide a thorough analysis of the sample complexity of our Tucker decomposition-based model and a variety of other results. Thus, our paper provides solid theoretical foundations for extending low-rank techniques to the nonparametric setting.

heading: 1 Introduction
text: Nonparametric density estimators are density estimators capable of estimating a density p while making few to no assumptions on p. Two commonly used nonparametric density estimators include the histogram estimator and the kernel density estimator (KDE). A common characteristic of these two estimators is that they make estimation tractable via a hyperparameter that relates to smoothness, namely bin width and bandwidth. Large bin width or bandwidth allows for good control of estimation error* at the cost of increased approximation error via large bin volume for the histogram and smoothing for KDEs. Selection of this parameter is crucial for estimator performance. In fact a recent survey on bandwidth selection for KDEs found at least 30 methods for setting this value along with a few surveys dedicated to the topic [15].
While nonparametric density estimation has been shown to be effective for many tasks, it has been observed empirically that estimator performance typically declines as data dimensionality increases, a manifestation of the curse of dimensionality. For the histogram and KDE this phenomenon has concrete mathematical analogs. For example, these estimators are only universally consistent† if
*For an estimator V restricted space a of densities P , the estimation error refers to the difference between ‖V − p‖ and minq∈P ‖p− q‖, where p is the target density. This is similar to estimator variance.
†A density estimator is universally consistent if it asymptotically recovers any density.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
n → ∞ and h → 0, with nhd → ∞, where n is the number of samples, d is the data dimension, and h is the bin width for the histogram and bandwidth parameter for the KDE [14]. One may then wonder whether there exists some other way to constrain model capacity so as to alleviate this exponential dependence on dimension.
In this paper we theoretically analyze the advantages of including a constraint akin to matrix/tensor rank in a nonparametric density estimator. We do so by analyzing histogram-style estimators (estimators that output a density that are piecewise constant on bins defined by a grid) that are enforced to have a low-rank PARAFAC or Tucker decomposition. Enforcing this low-rank constraint in the histogram estimator allows for much faster rates on h→ 0, while still controlling estimation error. This can remove the exponential penalty on rate of convergence that dimensionality produces with the standard histogram estimator.
The bulk of this work focuses on analysis of hypothetical estimators that, while offering large statistical advantages over the standard histogram estimator, are not computationally tractable. At the end of this work we include experiments demonstrating that low-rank histograms, computed using an off-the-shelf nonnegative tensor factorization library [23], consistently outperform the standard histogram estimator.

heading: 1.1 Density Models
text: Here we introduce the low-rank models that will be used to constrain our estimators. The first is a multi-view model. A multi-view model p is a finite mixture p = ∑k i=1 wifi whose components are separable densities fi = pi,1 ⊗ pi,2 ⊗ · · · ⊗ pi,d ‡ [19, 1, 37, 36]. A multi-view model has the following form
p (x1, x2, . . . , xd) = k∑ i=1 wipi,1 (x1) pi,2 (x2) · · · pi,d (xd) . (1)
A multi-view model with one component is typically called a naive Bayes model. For the estimators we propose, the component marginals pi,j will have the form of one-dimensional histograms and the number of components k will be limited to restrict estimator capacity. When k = 1 the model is equivalent to a naive Bayes model, p (x1, x2, . . . , xd) = p1 (x1) p2 (x2) · · · pd (xd). When the component marginals pi,j are histograms increasing k expands the set of potential estimates from naive Bayes models when k = 1 to all possible histogram estimates.
The models in this paper are motivated by nonnegative tensor factorizations so we term them generally as nonparametric nonnegative tensor factorization (NNTF) models. The previous model was related to nonnegative PARAFAC [34]. Our second model is based on the nonnegative Tucker decomposition [22]. A density in this model utilizes d collections, F1, . . . ,Fd, of k one-dimensional densities, Fi = {pi,1, . . . , pi,k}, and some probability measure that randomly selects one density from each Fi. This measure a can be represented by a tensor W ∈ Rk ×d where the probability of selecting (p1,i1 , . . . , pd,id) from F1×· · ·×Fd isWi1,...,id . To sample from this model we first randomly select the marginal distributions p1,i1 , . . . , pd,id according to W , and an observation is sampled according to the d-dimensional distribution p1,i1 ⊗ p2,i2 ⊗ · · · ⊗ p2,i2 . The density of this model is
p (x1, x2, . . . , xd) = k∑ i1=1 · · · k∑ id=1 Wi1,...,idp1,i1(x1)p2,i2(x2) · · · pd,id(xd). (2)
We are unaware of previous works investigating this model for general probability distributions so we will simply term it the Tucker model. Again we will investigate estimators where the component marginals are one-dimensional histograms and small k corresponds to reduced model capacity. We remark that Tucker decompositions typically have a rank vector [k1, . . . , kd] which for our model would mean that each Fi contains ki vectors and W would lie in Rk1×···×kd . This sort of rank restriction could be used in our methods however, for simplicity, we just set k1 = k2 = · · · = kd , k.
‡For two functions f, g let f ⊗ g : (x, y) 7→ f(x)g(y). This is analogous to the tensor product of L2 functions.

heading: 1.2 Overview of Results
text: The principal contribution of this paper is to analyze the advantage of incorporating rank restriction into density estimation. In Section 2.1 we precisely introduce the multi-view histogram and a histogram based on the Tucker model. With these models we investigate how quickly we can let bin width h go to zero and rank k grow, in relation to the amount of data n, while still being able to control estimation error and select an estimator that is nearly optimal for the allowable set. For the multi-view histogram we show that one can control estimation error so long as k/h is asymptotically dominated by n and for the Tucker histogram we need k/h + kd to be asymptotically dominated by n (we are omitting logarithmic factors here for convenience). This stands in stark contrast to standard space of histogram estimators which requires 1/hd to be asymptotically dominated by n. We furthermore show that these estimators are universally consistent and that these rates cannot be significantly improved.
For a second style of analysis we provide finite-sample bounds on the convergence of NNTF histogram estimators. We then construct a class of universally consistent density estimators that converge at rate Õ(1/ 3 √ n)§ on all densities that are a multi-view model with Lipschitz continuous component marginals. Note that the NNTF histogram estimators we construct do not require the target density to be an NNTF model to function well. Our estimators will select a good estimator so long as there exists any NNTF histogram that approximates the target density well; so our hypothetical estimators “fail elegantly” in some sense. We further show that the standard histogram can converge at a rate ω(1/ d √ n)¶ on the same class of densities. In Section 3 we experimentally investigate the efficacy of using NNTF histograms on real-world data. In lieu of the computationally intractable methods we investigated in our theoretical analyses, we use an existing low-rank nonnegative Tucker factorization algorithm to fit an NNTF histogram to data. Surprisingly even this method outperforms the standard histogram estimator with very high statistical significance. Lastly we mention that this paper is an extension of [42] and contains a fair amount of overlap with that text.

heading: 1.3 Previous Work
text: Nonparametric density estimation has been extensively studied with the histogram estimator and KDE being some of the most well-known methods. There do exist, however, alternative methods for density estimation, e.g. the forest density estimator [24] and k-nearest neighbor density estimator [25]. The L1, L2, and L∞ convergence of the histogram and KDE has been studied extensively [14, 9, 40, 16]. The KDE is generally regarded as the superior density estimator, with some mathematical justification [14, 35]. Numerous modifications and extensions of the KDE have been proposed including utilizing variable bandwidth [39], robust KDEs [20, 41, 44], methods for enforcing support boundary constraints [33], and a supervised variant [43]. The work [46] investigated using KDEs for nonparametric mixture modeling. Finally we mention [21] 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 the following analogy
k∑ i=1 wipi,1 (x1) pi,2 (x2) · · · pi,d (xd) ∼ k∑ i=1 λivi,1 ⊗ vi,2 ⊗ · · · ⊗ vi,d. (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, 30]. The most basic version of compressed sensing is concerned with estimating a “tall” vector x from a “short” vector y , Ax where A is a known “short and fat” matrix. One can recover y if it is sparse and A satisfies a property known as the restricted isometry property (RIP) [47, 7]. These methods can be extended to the estimation of matrices when x is a low-rank matrix [30, 26, 27]. In this extension A is an order-3 tensor which acts
§fn ∈ Õ(gn) ⇐⇒ ∃k s.t. fn ∈ O ( gn log k gn )
¶fn ∈ ω (gn) ⇐⇒ |fn/gn| → ∞ in probability. In d-dimensional space the standard histogram can converge slower than 1/ d √ n.
as a linear operator on x and satisfies an adjusted form of RIP. RIP commonly arises from matrices and tensors whose entries are random. Because of this compressed sensing techniques are useful in settings where one wants to estimate x from random linear transforms of x. For example, in matrix completion observing the (i, j)-th entry of x can be represented as an inner product of x with an indicator matrix ei,j , i.e. 〈x, ei,j〉F . Thus the random observed indices (i1, j1) , (i2, j2) , . . . can be represented as random matrices ei1,j1 , ei2,j2 , . . . which are then stacked into A. Now A is a random linear operator that is applied to x to represent the observation of random entries of x and the methods of compressed sensing can be used to recover x. Compressed sensing has also proven useful for multivariate regression and autoregressive models [26, 27]. Such techniques don’t appear to be extensible to histogram estimation due to the lack of a linear sampling scheme.
General matrix/tensor factorization, including nonnegative matrix/tensor factorizations, has been extensively studied despite being inherently difficult due to non-convexity. The identifiability and recovery of mixtures of product distributions was studied in [18, 38, 13]. The special case of mixtures of power distributions has also been investigated [29, 45]. The works [11, 5] present potential theoretical grounds for avoiding the computational difficulties of nonnegative matrix factorization. One notable approach to tensor factorization is to assume, in the tensor representation in (3), that d ≥ 3 and the collections of vectors v1,j , . . . ,vk,j are linearly independent for all j. Under this assumption we are guaranteed that the factorization (3) is unique [1, 18]. In [4] the authors present a method for recovering this factorization efficiently and demonstrate its utility for a variety of tasks. This work was extended in [36] to recover a multi-view KDE satisfying an analogous linear independence assumption and theoretically analyze the estimator’s convergence to the true low-rank components. In [36] 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. In [37] it was observed that using low-rank embeddings can improve density estimation. A multi-view histogram was investigated in [17] where the authors present an identifiability result and algorithm for recovering latent factors of the distribution. The related works [2, 3] consider a low-rank characteristic function as an approach to improving nonparametric density estimation. Though earlier works have observed that a low-rank approach improves nonparametric density estimation [37, 36, 17, 2, 3], we are the first to demonstrate this through theoretical analysis of sample complexity. Finally we note that the Tucker decomposition has been utilized in Bayesian statistics [32]. We are unaware of any literature on factoring functions in Rd → R in a Tucker-inspired as we do in (2).

heading: 2 Theoretical Results
text: In this section we mathematically demonstrate that histogram estimators can achieve greater performance by restricting to NNTF models and using a proper procedure to select a representative from these using data. To simplify analysis we will only consider densities on the unit cube [0, 1)d and analyze the number of bins per dimension b which is the inverse of the bin width, i.e. b = 1/h. To state our results precisely we must introduce a fair amount of notation.

heading: 2.1 Notation
text: We will denote the L1 and L2 Lebesgue space norms via a 1 or 2 subscript. Let Dd be the set of all densities on [0, 1)d. By density we mean probability measures that are absolutely continuous with respect to the d-dimensional Lebesgue measure on [0, 1)d. We define a probability vector or probability tensor to simply mean a vector or tensor whose entries are nonnegative and sum to one. Let ∆b denote the set of probability vectors in Rb and Td,b the set of probability tensors in Rb×d ||. The product symbol ∏ will always mean the standard outer product, e.g. set product** or tensor product, when the multiplicands are not real numbers††. The natural numbers N will always denote positive integers. For any natural number b let [b] = {1, . . . , b}. We will let 1 be the indicator function and Conv be the convex hull. Later we will use projection operator where ProjS x , arg mins∈S ‖x− s‖2; for every instance in this work this projection is unique.
||Rb ×d is the set of b× · · · × b︸ ︷︷ ︸ d times tensors. For example Rb ×2 is the set of b× b matrices.
**For sets S1, . . . , Sd we have ∏d
i=1 Si = S1 × · · · × Sd = {(s1, . . . , sd) : si ∈ Si∀i}. ††For functions f1, . . . , fd then ∏d i=1 fi = f1 ⊗ · · · ⊗ fd : (x1, . . . , xd) 7→ ∏d i=1 fi (xi).
We will now construct the space of histograms on [0, 1)d. We begin with one-dimensional histograms, which will serve as the pi,j terms in (1) or (2). We define h1,b,i with i ∈ [b] to be the one-dimensional histogram where all weight is allocated to the ith bin. Formally we define this as
h1,b,i (x) , b1 ( i− 1 b ≤ x < i b ) .
Note that this is a valid density due to the leading b coefficient. We use these to construct higherdimensional histograms. For a multi-index A ∈ [b]d, let
hd,b,A , d∏ i=1 h1,b,Ai ,
i.e. the d-dimensional histogram with b bins per dimension whose entire density is allocated to the bin indexed by A. Finally we define Λd,b,A to be the support of hd,b,A, i.e. the “bins” of a histogram estimator,
Λd,b,A , d∏ i=1 [ Ai − 1 b , Ai b ) .
For a dataset X = (X1, . . . , Xn) in [0, 1)d, the standard histogram estimator is
Hd,b (X ) , 1
n n∑ i=1 ∑ A∈[b]d hd,b,A1 (Xi ∈ Λd,b,A) .
Let Hd,b , Conv ({ hd,b,A ∣∣∣A ∈ [b]d}), the set of all d-dimensional histograms with b bins per dimension. LetHkd,b be the set of histograms with at most k separable components, i.e.
Hkd,b ,  k∑ i=1 wi d∏ j=1 pi,j ∣∣∣∣∣∣w ∈ ∆k, pi,j ∈ H1,b  . (4)
We will refer to elements in this space as multi-view histograms. Elements in this space have the same form as (1) in Section 1.1. Similarly we define the space of Tucker histograms to be
H̃kd,b =  ∑ S∈[k]d WS d∏ i=1 pi,Si ∣∣∣∣∣∣W ∈ Td,k, pi,j ∈ H1,b  .
These have the same form as (2) in Section 1.1.
We emphasize that the collections of densitiesHkd,b and H̃kd,b are primary objects of interest in this paper. The results we present are concerned with finding good density estimators restricted to these sets as k and b vary.

heading: 2.2 Estimator Theoretical Results
text: We present two approaches to the analysis of NNTF histogram estimators. All proofs and additional results are contained in the appendix.
First we provide an asymptotic analysis of NNTF histogram estimators in terms of estimation error control: how fast can we let b and k grow, with respect to n, while still controlling estimation error over all densities? An advantage of this analysis is that we can demonstrate that these rates are approximately optimal (up to logarithmic terms).
For our second approach we present finite-sample bound analysis. In this analysis we first present a distribution-dependent bound that, for an estimator restricted to Hkd,b, depends on n, b, k and minq∈Hkd,b ‖p− q‖1 where p is the data generating “target” distribution. We follow this up with distribution-free bounds.
Distribution-free bounds for nonparametric density estimation require that the target distribution belong to a well-behaved class of distributions (such as Sobolev or Hölder classes) to enable bounding of the approximation error [40]. Our distribution-free finite-sample analysis assumes that the target density is a multi-view model whose component marginals are Lipschitz continuous. We construct an estimator that converges at a rate of approximately 1/ 3 √ n on this class of densities, independent of d. For comparison we show that the standard histogram estimator can converge at a rate worse than 1/ d √ n on this same class of densities. We mention again that this estimator shouldn’t fail catastrophically when these distributional assumptions aren’t exactly met so long as its approximation error, minp∈Hkd,b ‖p− q‖1, isn’t large. For brevity, and because the results are virtually direct analogues of their multi-view histogram counterparts, we reserve all finite-sample results for Tucker histogram for the appendix.
Main Technical Tools Our results rely on finding L1 ε-coverings of the spaces of NNTF histograms and using Theorem 3.7 from [6] to select a good representative from that collection. The aforementioned theorem is a slight extension of Theorem 6.3 in [9] and is essentially proven in [48]. As was mentioned in [6], the application of these results typically does not yield a computationally practical algorithm. Likewise our results are simply meant to highlight the potential of NNTF models and are not practically implementable as is.

heading: 2.2.1 Asymptotic Error Control
text: 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 ∼ bk (omitting logarithmic factors). Recall that the standard histogram estimator requires n ∼ bd, so we have removed the exponential dependence of bin rate on dimensionality. Here and elsewhere the ∼ symbol is not a precise mathematical statement but rather signifies that the two values should be of the same order in a general sense. In the following b and k are functions of n so the space of histograms changes as one acquires more data. Theorem 2.1. For any pairs of sequences b→∞ and k →∞ satisfying
n/(bk log(b) + k log(k))→∞, there exists an estimator Vn ∈ Hkd,b such that, for all ε > 0
sup p∈Dd P
( ‖Vn − p‖1 > 3 min
q∈Hkd,b ‖p− q‖1 + ε
) → 0,
where Vn is a function of X1, . . . , Xn iid∼ p.
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 one can control the error of Tucker histogram estimates so long as n ∼ bk + kd (omitting logarithmic factors). Theorem 2.2. For any pairs of sequences b→∞ and k →∞ satisfying
n/ ( bk log(b) + kd log ( kd )) →∞,
there exists an estimator Vn ∈ H̃kd,b such that, for all ε > 0
sup p∈Dd P
( ‖Vn − p‖1 > 3 min
q∈H̃kd,b ‖p− q‖1 + ε
) → 0,
where Vn is a function of X1, . . . , Xn iid∼ p.
Allowing b to grow as aggressively as possible we achieve consistent estimation so long as n ∼ b log b and k grows sufficiently slowly, regardless of dimensionality.
Corollary 2.1. For all d, b, k fix Rkd,b to be either Hkd,b or H̃kd,b‡‡. For any sequence b → ∞ with n/ (b log b)→∞, there exists a sequence k →∞ and estimator Vn ∈ Rkd,b such that, for all ε > 0
sup p∈Dd P
( ‖Vn − p‖1 > 3 min
q∈Rkd,b ‖p− q‖1 + ε
) → 0,
‡‡R is fixed toH or H̃ and doesn’t change as n, b, k vary.
where Vn is a function of X1, . . . , Xn iid∼ p.
The following result shows that the approximation error of the estimators in Theorem 2.1, Theorem 2.2, and Corollary 2.1 go to zero for all densities. Thus these estimators are universally consistent even when the NNTF model assumption is not satisfied.
Lemma 2.1. Let p ∈ Dd. If k →∞ and b→∞ then minq∈Hkd,b ‖p− q‖1 → 0.
A straightforward consequence of this is that the Tucker histogram approximation error also goes to zero.
Lemma 2.2. Let p ∈ Dd. If k →∞ and b→∞ then minq∈H̃kd,b ‖p− q‖1 → 0.
The next theorem shows that the rate on bk in Theorem 2.1 cannot be made significantly faster.
Theorem 2.3. Let d ≥ 2, b → ∞, and k → ∞ with b ≥ k and n/ (bk) → 0. There exists no estimator Vn ∈ Hkd,b such that, for all ε > 0, the following limit holds
sup p∈Dd P
( ‖Vn − p‖1 > 3 min
q∈Hkd,b ‖p− q‖1 + ε
) → 0,
where Vn is a function of X1, . . . , Xn iid∼ p.
Likewise the rate on bk + kd can also not be significantly improved in Theorem 2.2. Theorem 2.4. Let d ≥ 2, b→∞, and k →∞ with b ≥ k and n/ ( bk + kd ) → 0. There exists no
estimator Vn ∈ H̃kd,b such that, for all ε > 0, the following limit holds
sup p∈Dd P
( ‖Vn − p‖1 > 3 min
q∈H̃kd,b ‖p− q‖1 + ε
) → 0,
where Vn is a function of X1, . . . , Xn iid∼ p.

heading: 2.2.2 Finite-Sample and Distribution-Independent Bounds
text: In this section we investigate the convergence of multi-view histogram estimators to densities that satisfy the multi-view assumption. This will be done via a standard bias/variance decomposition style of argument. We begin with the following distribution-dependent finite-sample bound.
Proposition 2.1. Let d, b, k, n ∈ N and 0 < δ ≤ 1. There exists an estimator Vn ∈ Hkd,b such that
sup p∈Dd P ‖p− Vn‖1 > 3 min q∈Hkd,b ‖p− q‖1 + 7 √ 2bdk log(4bdkn) n + 7 √ log( 3δ ) 2n  < δ, where Vn is a function of X1, . . . , Xn iid∼ p.
To analyze the approximation error, ‖p − q‖1, we first consider the case where p has a single component with L-Lipschitz marginals, so p = ∏d i=1 pi. Using the indicator function over the unit cube with Hölder’s Inequality we have that ‖p− q‖1 = ‖ (p− q)1‖1 ≤ ‖p− q‖2 ‖1‖2 = ‖p− q‖2. It is possible to show that the L2 projection of p onto H1d,b is achieved by simply projecting each marginal to its best approximating one-dimensional marginal i.e. (see the appendix)
arg min q∈H1d,b ∥∥∥∥∥ d∏ i=1 pi − q ∥∥∥∥∥ 2 = d∏ i=1 ProjH1,b pi.
Let LipL be the set of L-Lipschitz probability density functions [0, 1] and let mL , supf∈LipL ‖f‖2. The following theorem characterizes the approximation error for separable densities with Lipschitz continuous marginals.
Theorem 2.5. Let b2 ≥ L2/12 and f1, . . . , fd be elements of LipL then∥∥∥∥∥ d∏ i=1 fi − ProjH1d,b d∏ i=1 fi ∥∥∥∥∥ 1 ≤ √ m2dL − ( m2L − L2 12b2 )d .
We show in the appendix that this decays at a rate ofO ( b−1 )
and that this rate is tight. To characterize the behavior of mL we have the following.
Proposition 2.2. Let mL = supf∈LipL ‖f‖2, then
m2L =
{ L2/12 + 1 0 ≤ L ≤ 2√
8L 9 L ≥ 2
.
The following theorem is a finite-sample bound for a well-constructed multi-view histogram estimator.
Theorem 2.6. Fix n, d, k ∈ N, L ≥ 2, and 1 ≥ δ > 0. There exists b and an estimator Vn ∈ Hkd,b such that
sup p∈Q P ‖Vn − p‖1 > 21dk1/3L d+312 n 1 3 √ log(3Ldkn) + 7 √ log( 3δ ) 2n  < δ, where Q is the set of densities of the form ∑k i=1 wi ∏d j=1 pi,j with pi,j ∈ LipL, w ∈ ∆k, and Vn is a function of X1, . . . , Xn iid∼ p.
Here we analyze the asymptotic rate at which the estimator in Theorem 2.6 converges to a large class of multi-view models. To this end let
Qd =  k′∑ i=1 wipi,1 ⊗ · · · ⊗ pi,d | k′ ∈ N, w ∈ ∆k′ , L′ ≥ 0, pi,j ∈ LipL′  , i.e. Qd is the space of all multi-view models whose component marginals are all Lipschitz continuous. Consider letting L→∞, k →∞, and δ → 0 in Theorem 2.6 arbitrarily slowly as n→∞. For some element p in Qd its respective maximum Lipschitz constant and component number, L′ and k′, are fixed, so for sufficiently large n we have L > L′ and k > k′ and the bound from Theorem 2.6 applies. From this it follows that ‖Vn − p‖1 ∈ Õ (1/ 3 √ n). So we can construct an estimator that converges at rate Õ (1/ 3 √ n) for any multi-view model in Qd, independent of dimension! This rate appears to be approximately optimal since “for smooth densities, the average L1 error for the histogram estimate must vary at least as n−1/3” ([14], p. 99). For comparison the histogram estimator’s convergence rate is hindered exponentially in dimension.
Proposition 2.3. Let Vn be the standard histogram estimator with n/bd →∞. There exists p ∈ Qd such that ‖Vn − p‖1 ∈ ω(1/ d √ n).

heading: 2.3 Discussion
text: While Theorem 2.6 gives an estimator with good convergence, the class of densities Q is somewhat restrictive and likely not realistic for many situations where one would want to apply a nonparametric density estimator. Proposition 2.1, on the other hand, is not so restrictive since it depends on minq∈Hkd,b ‖p− q‖1, and more clearly conveys the message of this paper. Typical works on multiview nonparametric density estimation assume that the target density p is a multi-view model and are interested in recovering the model components pi,j and w from (1). Similarly to how the standard histogram estimator doesn’t assume p ∈ Hd,b our work isn’t meant to assume that p is a multi-view model, but is instead meant to explore the benefits of including a hyperparameter k, in addition to b, to restrict rank of the estimator. Just as the histogram can approximate any density as b → ∞, Lemmas 2.1 and 2.2 show that the inclusion of k does not hinder the approximation power of the estimator.
To explore the trade-off between b and k first observe that setting k = bd gives Hkd,b = Hd,b since one can allocate one component to each bin. Theorem 2.1 and Proposition 2.1 with k = bd gives a sample complexity of approximately n ∼ bd+1 which coincides with standard histogram estimator. Alternatively, setting k = 1 restricts the estimator to separable histograms H1d,b, with a sample complexity of approximately n ∼ b. 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.

heading: 3 Experiments
text: The previous section proved the existence of estimators that select NNTF histograms that offer advantages over the standard histogram. Unfortunately these estimators are not computationally tractable and only demonstrate the potential benefit incorporating an NNTF model in nonparametric density estimation. It would be nonetheless interesting to observe the behavior of an NNTF histogram estimator, even if it lacks the theoretical guarantees developed in the previous section. To this end we consider an L2-minimizing NNTF histogram estimator. For all d, b, k fixRkd,b to be eitherHkd,b or H̃kd,b. We consider an estimator Un, that attempts to minimize
Un = argminp̂∈Rkd,b ‖p̂− p‖ 2 2 .
Note that ‖p̂− p‖22 = ‖p̂‖ 2 2 − 2 〈p, p̂〉+ ‖p‖ 2 2 . (5)
The ‖p‖22 term in (5) is not relevant when minimizing over p̂. Additionally for data X1, . . . , Xn iid∼ p the law of large numbers gives us
〈p, p̂〉 = ∫ [0,1)d p(x)p̂(x)dx = EX∼p [p̂ (X)] ≈ 1 n n∑ i=1 p̂(Xi)
so we may consider
Un , argminp̂∈Rd,b ‖p̂‖ 2 2 −
2
n n∑ i=1 p̂(Xi) (6)
as a practical version of (5) which is conveniently equivalent to nonnegative tensor factorization (see the appendix).
Risk expressions for nonparametric density estimation based on L2-minimization similar to (6) have appeared in previous works related to kernel density estimation [41, 8, 31]. Proving optimal rates on bandwidth in these settings seems challenging, however.
For our experiments we used the Tensorly library [23] to perform the nonnegative Tucker decomposition [22] with Tucker rank [k, k, . . . , k] which was then projected to the simplex of probability
tensors using [12]. We also performed experiments with nonnegative PARAFAC decompositions using [34, 23]. These decompositions performed poorly. This is potentially because the PARAFAC optimization is more difficult or the additional flexibility of the Tucker decomposition was more appropriate for the experimental datasets.

heading: 3.1 Experimental Setup
text: Our experiments were performed on the Scikit-learn [28] datasets MNIST and Diabetes [28], with labels removed. We use the estimated risk from (6) to evaluate estimator performance (which may cause negative performance values, however lower values always indicate smaller estimated L2distance). 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. We also performed experiments projecting each dataset onto a random subspace of dimension d. These random subspaces were constructed so that each additional dimension adds a new index without affecting the others. To do this, we first randomly select an orthonormal basis for each dataset that remains unchanged for all experiments v1, v2, . . .. Then to transform a point X to dimension d we perform the following transform: Xreduced dim. = [v1 · · · vd]T X . We consider both transforms since PCA may select dimensions where the features tend to be independent, e.g. 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 (again using (6)). For the estimators we tested all combinations using 1 to bmax bins per dimension and k from 1 to kmax. As d increased the best cross validated b and k value decreased, so we reduced bmax and kmax for larger d to reduce computational time, while still leaving a sizable gap between the best cross validated b and k and bmax and kmax across all runs of all experiment. For d = 2, 3 we have bmax = 15 and kmax = 10; for d = 4 we have bmax = 12 and kmax = 8; for d = 5 we have bmax = 8 and kmax = 6. For parameter fitting we used random subset cross validation repeated 80 times using 40 of the 200 samples to evaluate the performance of the estimator fit using the other 160 samples. Performing 80 folds of cross validation was necessary because of the high variance of the histogram’s estimated risk. 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 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 two-tailed 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.

heading: 3.2 Results
text: Our results are presented in Table 1. Apart from the density estimators’ performance (“Hist. Perf.” and “Tucker Perf.”) this table contains the mean and standard deviation over the 32 trials for the optimal cross validated parameters for the estimators. This includes the number of bins, corresponding to b from the earlier sections, and the Tucker rank k. We see that the Tucker histogram usually outperforms the standard histogram (“Tucker Perf.” is less than “Hist. Perf.”) with high statistical significance (small p-val). As one would expect, the Tucker histogram cross validates for more bins, presumably reducing the k in exchange for more bins. Note that the MNIST PCA experiments always cross validated for the largest k and yielded less statistically significant performance increases (outside of the unusual Diabetes PCA d = 5 experiment). Presumably this experiment fit the NNTF assumption the least well (hence the large k) and thus received the least benefit.

heading: 4 Conclusion
text: In this paper we have theoretically and experimentally demonstrated the advantages of including rank restriction into nonparametric density estimation. In particular, rank restriction may be a way to overcome the curse of dimensionality since it reduces the sample complexity penalty incurred from dimensionality from exponential to linear. This paper is an initial theoretical foray for demonstrating the potential of combining rank restriction with nonparametric methods.
Acknowledgments and Disclosure of Funding
Robert A. Vandermeulen acknowledges support by the Berlin Institute for the Foundations of Learning and Data (BIFOLD) sponsored by the German Federal Ministry of Education and Research (BMBF). Antoine Ledent acknowledges support by the German Research Foundation (DFG) award KL 2698/5- 1 and the Federal Ministry of Science and Education (BMBF) award 01IS18051A. The authors thank Marius Kloft and Clayton Scott for useful discussions related to this work. The authors thank Nicholas D. Sidiropoulos for bringing to our attention some previous works on low-rank nonparametric density estimation.

Review summaries:
The authors propose and analyse a novel method for nonparametric density estimation. They provide a method for providing a density estimate with low rank structure which is a mixture of densities that factorizes across dimensions, with each factor defined as a one dimensional histogram. They provide conditions under which their estimator is consistent and characterize upper and lower bounds on its error.

The authors consider the density estimation problem in multiple dimensions. They study the convergence properties of a low-rank estimator: the estimator is the sum of a few tensor products of univariate histograms. They give asymptotic and finite sample bounds for convergence. When the true density is a sum of k tensor products of Lipschitz densities, the authors show an $L_1$ convergence rate of $n^{-⅓}.$   

The paper shows a number of results on non-parametric density estimation using mixtures of product histograms with two proposed methods called \textit{multi-view histograms} and \textit{Tucker histograms}. A simple experiment is also performed to show these methods work better than the usual histogram estimation in some specific cases. There are two major theoretical contributions as follows.  Firstly, the authors give some asymptotic error control for the proposed methods for any true density functions based on sample sizes, number of histograms' bins, and number of mixtures.  Secondly, they prove that when the true density is a mixture of product distributions, the proposed methods are more favorable compared to the traditional histograms method in terms of the convergence rate.

This work is concerned with the (high-dimensional) density estimation problem.  While most minimax rates  for  the  problem  are  expressed  in  terms  of  some  measure  of  smoothness  that  the  unknown  distribution is assumed to verify, the current approach adds more structure, by requiring that the joint density be tensor-decomposable (in a PARAFAC or in a Tucker sense), with Lipschitz components.  The authors show(theoretical analysis) that this subverts the smoothness-only minimax rates together with their curse of dimensionality. In order to illustrate their method, the authors further perform experiments in which –perhaps controversially– they relax the estimator search to an L2 problem.