Estimating Distributions with Low-dimensional Structures Using Mixtures of Generative Models

Modelling and estimating complicated high-dimensional distributions with low-dimensional structures remains one of the major challenges in modern statistical learning. Suppose we observe $n$ i.i.d. samples $\{X_{1},X_{2},\ldots,X_{n}\}$ living in an ambient Euclidean space $\mathbb{R}^{D}$ according to some unknown distribution $\mu^{*}$. We wish to estimate $\mu^{*}$ based on the samples for conducting statistical inference and generating new samples. One of the most popular nonparametric methods for distribution estimation is kernel density estimation (KDE). It has been shown that when $\mu^{*}$ admits a density function relative to the Lebesgue measure of $\mathbb{R}^{D}$, and the density function is $\beta$-smooth, then KDE can achieve the optimal rate $n^{-\frac{\beta}{2\beta+D}}$ for recovering the density value at any point in $\mathbb{R}^{D}$  (Silverman, 2018; Tsybakov, 2009). However, the non-parametric rate $n^{-\frac{\beta}{2\beta+D}}$ suffers from the curse of dimensionality as the ambient dimension $D$ appears in the rate exponent and can be enormous in machine learning applications involving images and texts (Brock et al., 2018; Oord et al., 2016). In order to avoid this exponential blow-up of the dimension, a common practice is to assume some additional structure in the data so that the effective dimension of the data space is relatively low.

One such structure that has attracted much attention recently is the so-called manifold hypothesis, which assumes the date to live on a $d$-dimensional submanifold $\mathcal{M}$ embedded in the possibly high-dimensional ambient space $\mathbb{R}^{D}$. Although submanifolds have more complicated geometry than the conventional Euclidean spaces, the manifold hypothesis is a natural assumption to make in a number of areas of science and technology. For example, in computer vision and medical imaging, data are usually images represented as vectorized pixel intensities. Although images may contain millions of pixels, it is usually determined by a comparatively smaller set of global characteristics such as camera projection, lighting condition, texture, object position and orientation. Other examples of high-dimensional complex data with low-dimensional manifold structures appear in natural language processing (Luo et al., 2020; Ling et al., 2017), protein-protein interaction detection (You et al., 2010; Terradot et al., 2004), and astronomy and shape analysis (Mardia, 1999; Jupp & Mardia, 2009).

Statistical theory and methodology for modeling manifold valued data have been developed in various contexts (Lin et al., 2020; Lin et al., 2017; Zhang et al., 2022; Lan et al., 2021; Divol, 2022; Tang & Yang, 2022; Berenfeld et al., 2022). Specifically, the problem of estimating a probability measure lying on an unknown low-dimensional Riemannian submanifold has been studied in a number of recent works. For example, Divol (2022) consider a kernel density type estimator based on a preliminary step of estimating the volume measure of the submanifold using local polynomial estimation techniques. They prove that the developed estimator can achieve the minimax-optimal error bound under the Wasserstein loss. Tang & Yang (2022) construct a two-step estimator: the first step estimates the data supporting submanifold; and the second step recovers the distribution on the estimated submanifold based on wavelet type estimators. They also show that such an estimation is minimax-optimal with respect to certain adversarial loss functions.  Berenfeld et al. (2022) develop a Bayesian procedure based on location-scale mixtures of Gaussians for estimating the density of data living close to an unknown submanifold with theoretical guarantees. However, although these existing methods are theoretically appealing, they usually have poor computational scalability with the ambient dimensionality and the sample size, making them costly to implement for modeling massive and high-dimensional real data, such as images and texts.

Auto-encoder based deep generative modeling approaches in the machine learning literature, such as variational auto-encoder (VAE) (Kingma & Welling, 2013; Rezende et al., 2014; Kingma et al., 2016), Wasserstein auto-encoder (WAE) (Tolstikhin et al., 2019), InfoVAE (Zhao et al., 2019) and inferential Wasserstein generative adversarial networks (iWGAN) (Chen et al., 2022), have achieved great successes in generating synthetic realistic-looking images and texts, and are usually very efficient to implement. However, despite their empirical successes, a general theoretical framework explaining whether and how these generative modelling approaches benefit from the low-dimensional manifold structure is lacking, and it is also not clear whether these existing methods are theoretically optimal in the minimax sense. For example, the key step in the auto-encoder is the extraction of $d$ ($d\ll D$) latent features (via an encoder $Q:\,\mathbb{R}^{D}\to\mathbb{R}^{d}$) that can be used for accurately reconstructing the original data (via a decoder $G:\,\mathbb{R}^{d}\to\mathbb{R}^{D}$). In other words, these auto-encoder based methods implicitly assume data $X$ to have a low-dimensional structure so that they can be accurately reconstructed in the sense that $X\approx G(Q(X))$. Moreover, it is often the case that real-world data falls on a manifold that does not admit a global parametrization. For example, when the data space is a boundaryless manifold such as a sphere, or disconnected (Khayatkhoei et al., 2018). This lack of global parametrization makes conventional auto-encoder methods equipped with a single encoder/decoder pair incapable of recovering the entire data space without incurring distortions. Our empirical results (c.f. Fig. 1) also suggest that conventional auto-encoder based generative modelling approaches tend to generate off real-manifold samples with unrealistic appearances.

In this article, we propose a new generative modelling approach for learning manifold-supported distributions that is theoretically minimax-optimal, computationally efficient, and empirically promising in generating complicated yet realistic-looking data. Unlike most existing generative modelling procedures that rely on the strong assumption of the existence of a global parametrization of the data space, we employ multiple encoder/decoder pairs, where each pair corresponds to the parametrization of a local patch of the data supporting manifold. Moreover, we utilize the partition of unity technique for gluing local probability measures estimated in the patches to form a global estimation of the probability measure on the manifold. In addition, most existing methods simply plug in the data empirical distribution in constructing the objective function for defining a GAN (Goodfellow et al., 2014) type estimator, which may lead to theoretical deficiency due to the failure of taking the smoothness of the target distribution into account. We instead propose to plug in a smoothness-regularized version that provably improves the estimation accuracy. Concretely, we show that when the target distribution is $\alpha$-smooth and lies in a $\beta$-smooth $d$-dimensional submanifold in $\mathbb{R}^{D}$, then the corresponding estimator $\widehat{\mu}$ based on $n$ data points achieves a non-asymptotic error bound of order $O\big(\frac{\log n}{\sqrt{n}}\vee n^{-\frac{\alpha\wedge(\beta-1)+1}{2(\alpha\wedge(\beta-1))+d}}\big)$ (here $a\wedge b$ denotes $\min\{a,b\})$ under the $1$-Wasserstein distance, which corresponds to the minimax rate modulo a logarithmic factor when $\alpha\leq\beta-1$. The implied rate of convergence does not suffer from the “curse of dimensionality” and only depends on the intrinsic dimensionality $d$ of the data. Our numerical results also show that the proposed method tends to be more accurate than conventional auto-encoder based generative modelling approaches and classic kernel density estimators for learning target distributions with low-intrinsic dimensional structures.

From discussions in Section 2, we see that conventional auto-encoder based generative modelling approaches may suffer from low representation power when the target distribution to be estimated lies on a general submanifold without global parametrization. However, the property that any manifold-supported distribution can be expressed in the form of a mixture of generative models (c.f. decomposition (3)) motivates us to employ multiple encoder/decoder pairs, and to use the partition of unity to glue them together with proper weights.

Recall that we have a set of $n$ i.i.d observations $X^{(n)}=\{X_{1},\ldots,X_{n}\}$ sampled from the target distribution $\mu^{*}$ lying on a $d$-dimensional submanifold $\mathcal{M}$ embedded in $\mathbb{R}^{D}$ with $d\leq D$. Let $\{S_{k}\}_{k\in[K]}$ be a suitably chosen open cover to $\mathcal{M}={\rm supp}(\mu^{\ast})\subset\mathbb{R}^{D}$, fix a partition of unity $\{\rho_{k}\}_{k\in[K]}$ subordinate to $\{S_{k}\}_{k\in[K]}$,³³ 3 Here we may consider any function $\rho_{k}:\mathbb{R}^{D}\to[0,1]$ so that $\text{supp}(\rho_{k})\subset S_{k}\subset\mathbb{R}^{D}$ and $\sum_{k\in[K]}\rho_{k}(x)=1$ for any $x\in\cup_{k\in[K]}S_{k}$. Note that $\{\rho_{k}|_{S_{k}\cap\mathcal{M}}\}_{k\in[K]}$ would form a partition of unity to $\mathcal{M}$. which can be chosen without the knowledge of $\mathcal{M}$ (c.f. Section 4). For any generic approximation family $\mathcal{G}$ consists of $(\mathbb{G},\mathbb{Q},\mathbb{v})$, where $\mathbb{G}=(G_{1},G_{2},\cdots,G_{K})$ with $G_{{k}}:\mathbb{R}^{d}\to\mathbb{R}^{D}$, $\mathbb{Q}=(Q_{1},Q_{2},\cdots,Q_{K})$ with $Q_{{k}}:\mathbb{R}^{D}\to\mathbb{R}^{d}$, and $\mathbb{v}=(\nu_{1},\nu_{2},\cdots,\nu_{K})$ with $\nu_{{k}}\in\mathcal{P}(\mathbb{R}^{d})$, we define the following estimator, which we call mixture of latent distribution matched auto-encoder (MLDMAE) estimator:

where recall that $c(\cdot,\cdot)$ is the cost function, $\widetilde{\nu}_{k,Q_{k}}$ is a (smoothness-regularized) estimator to the density of $(Q_{k})_{\#}(\frac{\mu^{*}\cdot\rho_{k}}{\mathbb{E}_{\mu^{\ast}}[\rho_{k}]})$ (the precise definition is available in Appendix C), and $\mathcal{D}(\cdot,\cdot)$ is a generic discrepancy measure between distributions on the latent space. Different from conventional LDMAE estimators, MLDMAE can also use empirical Bayes method to select data-dependent prior distributions for local latent variables, which adds extra flexibility in the modeling and may potentially reduce the approximation error. We show in Theorem 1 that for some carefully chosen approximation family $\mathcal{G}$, cost function $c$, discrepancy metric $\mathcal{D}$, and smoothness-regularized estimator $\widetilde{\nu}_{k,Q_{k}}$, the resulting estimator $\widehat{\mu}$ attains the minimax rate of convergence under the $1$-Wasserstein distance as $W_{1}(\widehat{\mu},\mu^{\ast})\leq C\,n^{-\frac{\alpha\wedge(\beta-1)+1}{2(\alpha\wedge(\beta-1))+d}}\vee\frac{\log n}{\sqrt{n}}$ when $\mu^{\ast}$ is an $\alpha$-smooth distribution on an unknown $\beta$-smooth $d$-dimensional submanifold.

The objective function of MLDMAE can be decomposed into two parts: the reconstruction cost $n^{-1}\sum_{k=1}^{K}\sum_{i=1}^{n}c\big(X_{i},G_{k}(Q_{k}(X_{i}))\big)\cdot\rho_{k}(X_{i})$ and the penalty $\sum_{k=1}^{K}\lambda_{k}\cdot\mathcal{D}\big(\widetilde{\nu}_{k,Q_{k}},\nu_{k}\big)$. We may also allow the latent dimension $d$ to be different across encoder/decoder pairs over $k\in[K]$. The reconstruction cost aims to learn local parametrizations of the supporting manifold of $\mu^{\ast}$ by enforcing the encoder/decoder pair $(\widehat{Q}_{k},\widehat{G}_{k})$ to represent some coordinate system $(\varphi_{k},\varphi_{k}^{-1})$ of local patch $U_{k}=\mathcal{M}\cap S_{k}$ of $\mathcal{M}$ in decomposition (3). Therefore, it corresponds to the support recovery of $\mu^{\ast}$. By employing multiple encoder/decoder pairs, the MLDMAE estimator avoids the restrictive global-parametrization assumption that is implicitly assumed in conventional LDMAE estimators. As a result, MLDMAE is suitable for a wider range of problems (see Fig. 1 for an illustration). On the other hand, the penalty term aims to enforce the reweighted local latent distribution $(\widehat{Q}_{k})_{\#}\big(\frac{\mu^{\ast}\cdot\rho_{k}}{\mathbb{E}_{\mu^{*}}[\rho_{k}]}\big)$ to match some member $\widehat{\nu}_{k}$ in the pre-specified prior family, so that $\widehat{\nu}_{k}$ is close to the $\nu_{k}^{\ast}$ in decomposition (3).

In practice, instead of selecting the best data-dependent priors, we can also fix the prior as a simple distribution $\nu_{0}$ such as an isotropic Gaussian. Moreover, $\mathcal{D}(\cdot,\cdot)$ can be chosen as certain squared maximum mean discrepancy (MMD) loss⁴⁴ 4 For a positive-definite reproducing kernel $k$, the MMD loss is defined as ${\rm MMD}^{2}(\mu_{1},\mu_{2})=\mathbb{E}_{X,X^{\prime}\in\mu_{1}}[k(X,X^{\prime})]+\mathbb{E}_{Y,Y^{\prime}\in\mu_{2}}[k(Y,Y^{\prime})]-2\mathbb{E}_{X\in\mu_{1},Y\in\mu_{2}}[k(X,Y)]$ that can be efficiently computed in a closed-form formula. The $k$-th smoothness-regularized distribution $\widetilde{\nu}_{k,Q_{k}}$ in (4) can be constructed by applying kernel-smoothing to its (weighted) empirical counterpart $(Q_{{k}})_{\#}\widehat{\mu}_{n}^{{k}}$ with $\widehat{\mu}_{n}^{{k}}=(n\widehat{p}_{k})^{-1}\sum_{i=1}^{n}\rho_{k}(X_{i})\,\delta_{X_{i}}$, leading to $\widetilde{\nu}_{k,Q_{k}}(z)=(n\widehat{p}_{k})^{-1}\sum_{i=1}^{n}\widetilde{k}(z,Q_{k}(X_{i}))\rho_{k}(X_{i})$ for a suitable kernel $\widetilde{k}$. Note that when kernel $\widetilde{k}$ is the Gaussian kernel $\widetilde{k}(x,y)=(2\pi h)^{-\frac{d}{2}}\exp(-\frac{\|x-y\|^{2}}{2h})$ with bandwidth parameter $h$, then $\widetilde{\nu}_{k,Q_{k}}$ corresponds to the Gaussian-smoothed distribution $(\widetilde{Q}_{{k}})_{\#}\widehat{\mu}_{n}^{{k}}$, where $\widetilde{Q}_{{k}}$ is the randomly perturbed encoder defined by $\widetilde{Q}_{k}(X)=Q_{k}(X)+\sqrt{h}\cdot\mathcal{N}(0,I_{d})$.⁵⁵ 5 Here for randomized map $Q$, the push forward measure $Q_{\#}\mu$ is defined as the measure so that for any measureable function $f$, $\int f(x)\,{\rm d}[Q_{\#}\mu]=\mathbb{E}\big[\int f(Q(x))\,{\rm d}\mu\big]$ where the expectation is with respect to the randomness of $Q$. Employing such a smoothness-regularized distribution can be viewed as applying a randomized data augmentation to increase the variability of the encoded training samples, which mitigates potential overfitting to data and improves the generalization ability of the resulting estimator.

Introducing the encoder-decoder structure as in our estimator brings several benefits. Computationally, the encoders turn the high-dimensional data into low-dimensional latent variables so that we only need to compute a penalty term over low-dimensional distributions. Therefore, the MLDMAE framework brings less computational burden compared with generative modelling approaches (e.g. iWGAN) that directly deal with distributions in the ambient space. Theoretically, when $d\ll D$, the data distribution $\mu^{*}$ becomes a singular measure in $\mathbb{R}^{D}$. As a consequence, with the information about the supporting manifold of $\mu^{\ast}$, which is explicitly induced by the encoder-decoder pairs, it is possible to utilize classical techniques of nonparametric density estimation, such as wavelet truncation, to construct a minimax-optimal estimator. Specifically, the underlying true latent variable distribution $(Q_{k})_{\#}(\frac{\mu^{*}\cdot\rho_{k}}{\mathbb{E}_{\mu^{\ast}}[\rho_{k}]})$ defined in the mixture of generative models (3) is, with high probability, absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}^{d}$ (c.f. Lemma 4 in Appendix C), which enables us to develop smoothness-regularized estimators by borrowing techniques from Liang (2020) and Singh et al. (2018). Indeed, as suggested by Liang (2020) and Tang & Yang (2022), the rate $O\big(n^{-\frac{\alpha+1}{2\alpha+d}}\vee\frac{\log n}{\sqrt{n}}\big)$ achieved by the MLDMAE estimator when $\alpha\leq\beta-1$ is minimax-optimal up to logarithmic factor relative to the $1$-Wasserstein distance.

Based on the above discussion, we can now develop the algorithm as described in Algorithm 1 for implementing the MLDMAE estimation, where we use $\mathbb{G}_{\theta}=\{G_{\theta,1},G_{\theta,2},\cdots,G_{\theta,M}\}$ and $\mathbb{Q}_{\phi}=\{Q_{\phi,1},Q_{\phi,2},\cdots,Q_{\phi,M}\}$ to denote the decoders and encoders parametrized by $\theta$ and $\phi$ respectively.

In this section, we derive the finite sample error of the MLDMAE estimator. We first state the following assumptions on the target distribution $\mu^{\ast}$ and the approximation family used in defining the estimator (4).

In this section, we present some visual results of the MLDMAE approach when apply to common manifolds: $2$D-spiral and torus. The precise data generating distributions are given in Appendix A. The penalty term $\mathcal{D}(\cdot,\cdot)$ is chosen to be the $1$-Wasserstein distance computed by the “POT” package (Flamary et al., 2021) and the cost function $c(\cdot,\cdot)$ is the squared $\ell_{2}$ loss. As benchmarks, we also consider (1) the classic kernel density estimator (KDE) with Gaussian kernel that is commonly employed in statistics literature for density estimation; (2) the LDMAE estimator with the same kind of cost function and discrepancy metric as MLDMAE. The generated samples from the generators learned by different approaches are given in Fig. 3. The $1$-Wasserstein distance between the estimated distribution and the true distribution are given in Table 1. We can see that employing multiple encoders and decoders can lead to much better performance than employing a single pair of encoder and decoder. In particular, we can see even though we increase the complexity of the encoder/decoder family, LDMAE still can not capture the correct shape of these standard manifolds in statistics. On the other hand, by allowing multiple encoders/decoders in MLDMAE, the manifold structure can be correctly learned with a relatively simple encoder/decoder structure and a smaller number of training parameters. Moreover, MLDMAE can beat the classic KDE in both examples.

In this section, we empirically evaluate the proposed MLDMAE approach using three real-world datasets: MNIST handwritten digit (LeCun et al., 1995), Fashion-MNIST (Xiao et al., 2017) and CelebA $64\times 64$ (Krizhevsky et al., 2009). For comparison, we consider LDMAE approach and variational auto-enocoder (VAE) (Kingma & Welling, 2013; Rezende et al., 2014; Kingma et al., 2016), which are commonly used auto-encoder based generative modelling approaches. In all reported experiments, we set the reconstruction cost to be the squared $\ell_{2}$ loss, and fix the priors $\nu_{{k}}$ to be a standard Gaussian $\mathcal{N}(0,I_{d})$. The partition of unity is chosen to be the indicator functions described in Section 4. The encoders and decoders are modelled by convolutional deep neural networks with parameter sharing as described in Section 4, and further details are available in Appendix A. We consider two kinds of discrepancy metric $\mathcal{D}(\cdot,\cdot)$ in the penalty terms of LDMAE and MLDMAE, one is the $1$-Wasserstein distance computed through “POT” package, and the other one is MMD with the inverse multiquadratics (IMQ) kernel $k(x,y)=2d/(2d+\|x-y\|^{2}_{2})$. As described in Tolstikhin et al. (2019), the inverse multiquadratics kernel has a much heavier tail than the conventional RBF kernel $k(x,y)=\exp(-\|x-y\|^{2}/C)$, so it can provide more meaningful gradients for outliers. The numbers $K$ of clusters for MLDMAE are selected so that the number of free parameters is around $\frac{1}{5}\sim\frac{1}{2}$ of the number of sharing parameters, and it turns out $K=5$ works well for all the examples. Another important factor is the latent dimension, which is not explicit for the real dataset. Selecting a too small latent dimension would lead to large reconstruction errors and thus result in noisy generated samples. On the contrary, selecting a too large latent dimension would lead to the singularity of the encoded distribution and thus result in numerical instabilities. We use $d=4$ for Fashion-MNIST, $d=8$ for MNIST handwritten digit and $d=64$ for CelebA, which seems to work reasonably well, and trainings of MLDMAE are stable and robust to initialization.