Flexible Bayesian Nonparametric Product Mixtures for Multi-scale Functional Clustering

Functional data analysis is perhaps one of the most prominent areas of research in statistics and machine learning (Ramsay & Silverman, 2005). This is not surprising given the vast amounts of functional data collected in practice, ranging from time-series data (Aghabozorgi et al., 2015), spatial data (Mateu & Romano, 2017), spatio-temporal data (Atluri et al., 2018), among others. In contrast to a rich literature for analyzing functional datasets (Wang et al., 2015), Bayesian literature on functional data clustering is (somewhat surprisingly) not very well-developed (Zhang & Parnell, 2023). Most Bayesian functional data clustering methods rely on parametric finite mixture models (Bouveyron et al., 2019; Chamroukhi & Nguyen, 2019) or more flexible non-parametric infinite mixtures (Frühwirth-Schnatter & Malsiner-Walli, 2019).

Bayesian functional clustering approaches typically assign a global membership, creating subgroups in which functions share the same cluster membership within each subgroup. These global clustering methods produce distinct groups of functions defined by overarching patterns; however, they may not capture the heterogeneity in localized features within the functions lying in the same cluster. Some examples of global clustering methods include (Ray & Mallick, 2006; Crandell & Dunson, 2011; Bigelow & Dunson, 2009; Rodriguez & Dunson, 2014). We refer the readers to Wade (2023) for a more comprehensive review of such approaches. The majority of Bayesian global clustering methods have focused their implementation on one-dimensional functions such as time-series curves or longitudinal trajectories. While extending these approaches to higher order functions (such as spatial functions or images) should be feasible in principle, it is not straightforward in practice. This is due to the fact that higher order functions involve high-dimensional parameters that are challenging to cluster using classical global clustering methods. In this context, Chandra et al. (2023) showed that classical Bayesian mixture modeling approaches produce global clustering that are susceptible to poor performance and degenerated clustering in high-dimensional settings. In particular, these approaches may potentially produce an unrealistic large number of clusters or collapse to just one cluster when the model dimension far exceeds the sample size. While these results in Chandra et al. (2023) were developed with respect to clustering of multivariate vectors, they are naturally generalizable to our functional data analysis settings of interest.

We conjecture that a potential reason for inadequate performance of global clustering approaches in high-dimensional settings may stem from their inability to account for: (i) heterogeneity of localized features (Di Iorio & Vantini, 2023); and (ii) scenarios where a small subset of important features separate the clusters (Floriello & Vitelli, 2017). To circumvent difficulties associated with global clustering, a limited number of local clustering methods have been proposed that enable flexible clustering in different resolutions or domains (e.g. Dunson, 2009; Petrone et al., 2009; Rodríguez et al., 2010). These approaches produce distinct clustering allocations across various subsets of model parameters characterizing the function and thereby relax the restrictions inherent in global clustering. This enables one to account for greater heterogeneity and flexibility when clustering high dimensional functions, which is desirable.

Local clustering approaches can be broadly categorized into two classes. The first class of methods often directly cluster the element-wise mean parameters corresponding to the different observed instances within the function domain, while incorporating spatial dependence via these mean parameters. Petrone et al. (2009) proposed a hybrid Dirichlet process prior that formulates the individual functions as hybrids of global functions drawn from a Gaussian process. A latent Gaussian copula is used to allow local surface selection. Related alternative approaches were proposed by Rodríguez et al. (2010), and Nguyen & Gelfand (2011). The above approaches are quite useful, but can be computationally intensive due to the need to update as many labeling indicator vectors as the number of locations within MCMC sampling. The second class of methods rely on basis expansions to smooth the mean function, while imposing suitable mixture priors on the basis coefficients. Dunson (2009) proposed local partition processes (LPP) that enabled independent clustering across basis coefficients, while also specifying an additional shared mixture component for global clustering. Suarez & Ghosal (2016) proposed a spike and slab prior independently for wavelet basis coefficients, and proposed a post-MCMC decision-theoretic criteria to consolidate the local clusters to inform global clustering. Recently, Fan & Sarkar (2023) proposed a hidden Markov model for local partitioning and clustering of multiple time-series data.

While local clustering approaches are potentially more adept at modeling high-dimensional functions, there are still several unmet challenges. This is exemplified by an increase in the number of possible clusters to $\prod_{j=1}^{p}k_{j}$ when clustering a function with $p$ observed points, where $k_{j}$ denotes the number of mixture components for the $j$th coefficient ($j=1,\ldots,p$). Therefore, the increase in flexibility due to a greater number of mixture components under local clustering comes at a cost of greater computational complexity that must be controlled via careful modeling considerations. Moreover, the theoretical implications of Bayesian local clustering methods are not particularly well-understood, to our knowledge. Using wavelet expansions, Suarez & Ghosal (2016) established asymptotic consistency in recovering the true (finite dimensional) model parameters as well as the true partition structure in the ideal case when the observed measurement error associated with the function approaches zero. Unfortunately, the assumption of vanishing noise may not be supported in most practical settings. Nguyen & Gelfand (2011) used a finite mixture model and explored posterior consistency for density estimation in the limiting case when the number of components $k\to\infty$. While useful, they assumed the same form of the true density as the working model for their theoretical derivations that may be restrictive. Further, the working model based on finite-dimensional representations (via truncation) is deemed less flexible compared to the infinite dimensional mixture representation that is ultimately used for their theoretical analyses. Beyond these limited findings, posterior consistency guarantees pertaining to density estimation of random functions under Bayesian non-parametric local clustering approaches are limited or essentially absent to our knowledge.

Our focus is on developing practical non-parametric Bayesian product mixture models for local clustering of spatial functions that is scalable to high-dimensions and to investigate the theoretical properties such as posterior consistency. We leverage classical basis representations (such as wavelets) of functions to induce spatial smoothing for the observed noisy functional data, while simultaneously allowing for spatially correlated measurement errors for a more flexible characterization. To encourage local clustering among the high-dimensional basis coefficients, we apply independent Dirichlet process (DP) priors to coefficients corresponding to distinct wavelet resolutions. This structure introduces resolution-specific clustering and increases flexibility by delineating the clustering of high-level (or global) features and finer resolution (or local) features. By allowing the global features to drive the overall clustering process, the method diminishes the impact of local features on the clustering mechanism. This strategy is particularly useful when local features introduce additional noise and heterogeneity, or when the clusters are separated by a small subset of global features only. The proposed approach facilitates dimension reduction by dividing the basis coefficients into smaller subsets, which enhances model parsimony, reduces computational complexity, and improves scalability. Although our method is reminiscent of the local partitioning approach by Dunson (2009), there are important distinctions in terms of the local partitioning methodology and theoretical contributions, as discussed in Section 2.3.

The practical advantages of the proposed approach are complimented by desirable posterior consistency properties for density estimation, under mild assumptions. These results are established via careful sieve constructions and deriving corresponding entropy bounds and builds upon recent developments in non-parametric Bayesian density estimation literature (Canale & De Blasi, 2017; Kundu & Lukemire, 2024) focused on multivariate data. We develop an efficient MCMC algorithm for implementation of the proposed approach based on wavelet representations. Simulation studies illustrate strong advantages of the proposed approach assuming correlated or independent errors compared to global clustering approaches, in terms of clustering. The proposed method is applied to a spatial transcriptomics data focused on clustering the spatially varying gene expression profiles in a breast cancer study where it infers interpretable and reproducible gene clusters. While we focus on spatial functions for this article, the proposed approach is broadly applicable to a wide class of functions that admit suitable basis representations.

In summary, our approach based on non-parametric Bayesian product mixture models for high-dimensional functional data introduces four main contributions:

• a.

  Local clustering to address the curse of dimensionality resulting in improved computation complexity: Using independent DPs on basis coefficients corresponding to distinct wavelet resolutions, the proposed approach partitions the parameter space and results in local clustering across resolutions and global clustering within a given resolution. The resolution-specific priors reduce the dimensionality and improve the computation complexity as illustrated in Lemma 1.

• b.

  Theoretical justifications based on posterior consistency: Under mild conditions, we establish posterior consistency for density estimation of the random function (Lemma 2 and Theorem 2) that is satisfied under practical prior specifications (Corollary 2).

• c.

  Inclusion of correlated noise: Our framework incorporates correlated residual noise terms (with independent errors as a special case) that allows an improved characterization of local perturbations that are unrelated to clustering and may accentuate heterogeneity. We introduce low-rank representations of the error covariance that improves scalability to high-dimensions while retaining desirable theoretical properties.

• d.

  Efficient Gibbs sampler : We implement an efficient Markov chain Monte Carlo (MCMC) algorithm (Section 4) that facilitates model implementation for high dimensional spatial functions.

Among the above notable contributions, the inclusion of correlated noise (c) is motivated by recent literature on measurement error models (Ma & Kundu, 2022). Existing Bayesian functional clustering literature typically ignores correlated error terms distributed over spatial locations, which may amplify local perturbations and ultimately hinder the ability of clustering methods to delineate clear patterns in the data. The proposed method addresses this gap by not only incorporating correlated measurement errors, but also allowing these correlations to vary across subjects to accommodate additional heterogeneity.

The article is structured as follows. We develop the proposed methodology for noisy functional data in Section 2. Section 3 shows the theoretical guarantees including posterior consistency, while Section 4 outlines posterior inference strategy. In Section 5, we present results from extensive simulation studies. Section 6 applies the approach to a spatial transcriptomics data for clustering gene expressions. The Discussion section comments on potential limitations of the proposed methods as well as extensions of this approach to other contexts. Supplementary materials contain proofs of theoretical results, posterior computation steps, and additional details.

The proposed method defines a prior on the observed functions $y_{i}(v),v\in\operatorname{\mathcal{V}}$ of Equation (6) that satisfies several desirable theoretical properties. In our theoretical analysis, particularly for deriving posterior consistency results, we omit the subscript $i$ when appropriate, following standard convention. All theoretical results are constructed under the conditions of $n\rightarrow\infty$ and fixed $L$. For our theoretical treatment, we assume that the number of observed support points for the spatial function is $2^{J+1}$ that produces a square orthonormal basis matrix $\Psi(L\times L)$. In general, the theoretical results hold for all orthogonal basis representations with the product of DP prior on basis coefficients, as long as the true density $f_{0}$ satisfies certain reasonable assumptions. Notably, our posterior consistency results generalized the rich theoretical results for Bayesian non-parametric multivariate density estimation to the case of Bayesian functional data analysis, for which there are limited results. In particular, our focus is on estimating densities of random functions admitting suitable basis expansions.

The detailed proof of Lemma 2 is available in Appendix. Notably, the assumptions on the true density in Lemma 2 are routinely assumed in multivariate density estimation literature (Wu & Ghosal, 2008). Under an orthonormal basis (resulting in a finite determinant and Frobenius norm of $\Psi_{j}$), the proof mainly follows from the results in Wu & Ghosal (2008). Following the argument of Schwartz (1965), a positive prior support for arbitrarily small KL neighborhoods of $f_{0}$ from Lemma 2 provides the weak consistency of fPDPM.

While weak consistency is useful, it is of interest to investigate whether fPDPM also results in a more salient strong consistency property. Strong consistency is desirable since it ensures that the posterior concentrates in arbitrarily small $L_{1}$ neighborhoods of $f_{0}$, thereby guaranteeing more reasonable numerical estimation. However, given that the space $\operatorname{\mathcal{F}}$ can be non-compact, additional conditions are required to establish strong consistency (Canale & De Blasi, 2017). The basic strategy follows the framework in Ghosal & van der Vaart (2017) by constructing sieves $\operatorname{\mathcal{F}}_{n}$ (a compact subset of $\operatorname{\mathcal{F}}$) that possesses certain required properties and that eventually grows to cover the whole space $\operatorname{\mathcal{F}}$ as $n\rightarrow\infty$. Next, we consider how to construct such sieves.

Denote $F_{0}^{n}$ as the $n$-product of the measure $F_{0}$ that is associated with the true density $f_{0}$ and $\log N(\epsilon,\operatorname{\mathcal{G}},d)$ as the entropy of the space of densities $\operatorname{\mathcal{G}}\subset\operatorname{\mathcal{F}}$, where $N(\epsilon,\operatorname{\mathcal{G}},d)$ is the minimum integer for which there exists densities $f_{1},\ldots,f_{N}\in\operatorname{\mathcal{F}}$ such that $\operatorname{\mathcal{G}}\subset\cup_{j=1}^{N}\{f:d(f,f_{j})<\epsilon\}$. For the metric $d$ used here, we consider the Hellinger distance, $d^{2}(f,g)=\int(\surd{f}-\surd{g})^{2}$, and $L_{1}$ distance, $\lVert f-g\rVert_{1}=\int\lvert f-g\rvert$. We state below the sufficient conditions for the sieves that are required for strong consistency and was originally proposed in Canale & De Blasi (2017). This result will be used in our theoretical calculations when deriving strong consistency.

Conditions (1A) and (1B) in Lemma 3 regulate the growth of sieve in $\operatorname{\mathcal{F}}_{n}$ when sample size $n$ increases. Specifically, (1A) suggests the prior probability assigned to $\operatorname{\mathcal{F}}_{n}^{c}$ needs to shrink exponentially in $n$, implying that the tails with only exponentially small probability are excluded. Condition (1B) can be considered as a summability condition that stipulates that the weighted sum of the entropy number (weighted by the square root of the prior probability) is required to go toward zero when $n$ increases. To achieve both conditions, we carefully construct sieves (compact sub-spaces) as outlined below.

Under the above construction, the entropy of the sieves as given by the following Theorem.

Detailed proof is available in the Appendix. Having established the entropy bounds in Theorem 1, we are now in a position to establish strong consistency results via the Theorem below.

The proof is a direct result of Corollary 1 of Canale & De Blasi (2017). In this paper, we generalize these priors to Laplace priors (Park & Casella, 2008) on the base measures $P^{*}(\beta)$ that achieve better shrinkage and therefore results in improved model fitting. In particular, we express the Laplace priors via scale mixtures of normals as

In addition for the subject-specific covariance matrices, we impose the low-rank factor model to further improve the scalability. In particular, we specify

where $\Lambda_{i}$ is the factor loading matrix of dimension $L\times K_{i}$ with $K_{i}\ll L$ as the number of factors for $i$th unit, and $P^{*}_{\Sigma}=P^{*}_{\Lambda}\times P^{*}_{\sigma}$. Denote $h_{is}$ as the membership of covariance for $i$th unit. We specify the base measure as

The multiplicative variance terms involves local shrinkages $\phi$ that control the shrinkage of individual factor loadings, and global terms $\zeta,e$ that control the number of factors and the cluster-specific shrinkage, respectively. A key advantage of the model specified in (8) is that the number of factors $K_{i}$ is decided by the data instead of pre-specification (Murphy et al., 2020) and can vary among different clusters. The above prior construction is borrows from the infinite mixture of factor analyzers considered in Murphy et al. (2020), and results in base measures on the covariance the satisfy conditions (2B)-(2D) in Theorem 2 as shown in the following Corollary. Ultimately, this implies that the resulting base measures leads to strong consistency.

Detailed proof is available in Appendix. Note that the case of independent covariance of $\Sigma_{i}=\sigma^{2}_{i}I_{L}$ can be easily obtain as a special case Corollary 2 with $\Lambda_{i}=0$.

We implement the fPDPM approach via an efficient Gibbs sampler that proceeds via a slice sampling technique. The full details are provided in Supplementary Material Section S4. The Gibbs sampler results in posterior samples of coefficients and memberships across all resolution levels and units. The posterior samples are summarized for the following analysis.

Clustering.  The fPDPM may result in varying cluster memberships for different resolution levels. One can consolidate these local clusters to inform global clustering patterns as follows. The global clustering coefficients are summarized based on pairwise agreement among the subject-specific coefficients by adapting the method in Suarez & Ghosal (2016). Denote $h^{(r)}_{i}=[h^{(r)}_{i0},\ldots,h^{(r)}_{iJ}]^{\mathrm{\scriptscriptstyle T}}$ as a vector of members with $h_{ij}$ as the membership at resolution level $j$ of unit $i$ for the $r$th posterior sample with $r=1,\ldots,R$. Define the pairwise distance between units $i$ and $i^{\prime}$ as a weighted sum as

where $w_{j}$ is the weight for the agreement of coefficients at the resolution level $j$. To focus on the global pattern and reduce the impact of the local features on the overall clustering, we assign a higher weight on the global coefficients and downweight the coefficients corresponding to the remaining resolution levels. We achieve this by setting $w_{j}=\frac{1}{j}$ when $2^{j}<n$ and $w_{j}=\frac{1}{2j}$ if $2^{j}\geq n$. Based on the pairwise distance $\bar{d}(h_{i},h_{i^{\prime}})$, we build an agglomerative clustering with the complete linkage and decide the memberships by cutting the tree at different levels.

Function estimation.  We report the accuracy of the estimated functions via the mean squared error (MSE) of the posterior mean. Denote $R$ posterior samples of coefficients across all resolution levels for unit $i$ as $\{\beta^{(r)}_{ij},j=0,\ldots,J,r=1,\ldots,R\}$. Let $\theta^{(r)}_{i}=\sum_{j=1}^{J}\Psi_{j}\beta^{(r)}_{ij}$ be the estimated function from the $r$th posterior sample of coefficients. We calculate the posterior mean function as $\hat{\theta}_{i}=\sum_{r=1}^{R}\theta^{(r)}_{i}/R=[\hat{\theta}_{i}(v_{1}),\ldots,\hat{\theta}_{i}(v_{L})]^{\mathrm{\scriptscriptstyle T}}$ and obtain the MSE by $\sum_{i=1}^{n}(\hat{\theta}_{i}-\theta_{i})^{2}/nL$.

In this paper, we present a flexible Bayesian nonparametric framework for local clustering, termed fPDPM, designed to cluster high-dimensional functional data across multiple resolution levels. By employing a wavelet basis expansion combined with a product of Dirichlet process priors on the basis coefficients, fPDPM enables independent clustering of functional data at different resolutions. This approach achieves a form of local clustering that reduces computational complexity compared to traditional local clustering methods, while effectively mitigating the curse of dimensionality often encountered in global clustering techniques. Additionally, we propose a low rank decomposition for the residual covariance matrix that provides an improved characterization and simultaneously improves computational complexity. Under mild conditions, we show the salient strong consistency of fPDPM and exemplify some common priors that satisfies the assumed conditions. We also propose an efficient Gibbs sampler for posterior inference. Our simulations demonstrate that fPDPM outperforms existing methods that focus on global clustering under a variety of settings. We use fPDPM on spatial transcriptomics data of breast cancer and identify replicable clusters comprising genes that corroborate existing biological literature.

Currently, fPDPM considers cross-sectional functional data under the normality error assumption. One generalization is to consider the longitudinal functional data such as spatio-temporal data (Atluri et al., 2018) and model the temporal correlation by allowing the coefficients to vary on the temporal relationships. Another possible extension is to extend the fPDPM approach for clustering based on images from multiple modalities. All these directions are left for future investigations.