1 Introduction↩︎

Machine learning systems inevitably face data that deviate from their training distributions. Generally, this data is either sparsely labeled or wholly unsupervised. Faced with such a dynamic environment, an intelligent system must accurately detect outliers and respond appropriately. This capability is the subject of modern research on out-of-distribution (OOD) detection [@HendrycksD17], anomaly detection [@chandola2009anomaly], and open-set recognition [@scheirer12openset].

Outlier detection is one of the oldest problems in statistics [@anscombe1960rejection], and there are well-established methods for tackling this canonical problem. For example, Bayesian nonparametric methods offer a coherent, probabilistic framework for estimating the probability that a data point belongs to a new cluster [@ferguson1973bayesian; @antoniak74dpmm; @lo1984class; @Sethuraman94; @maceachern1994estimating; @neal00]. Here, we use Dirichlet Process Mixture Models (DPMMs) to fit a generative model to the training data. Under this model, OOD detection reduces to a straightforward computation of the probability that a data point belongs to a novel class.

Care must be taken with Bayesian nonparametric methods, however. Modern machine learning systems are often built on foundation models that have been trained on massive datasets [@dosovitskiy20vit; @caron21dino; @oquab23dinov2; @darcet23needreg; @chen20simclr; @chen20simclrv2]. These models yield feature embeddings of data points that can be used for several downstream tasks, but the embeddings are high-dimensional. When fitting a Gaussian DPMM, for example, we must implicitly estimate the covariance of embeddings within and across classes, which presents both computational and statistical challenges. We propose a hierarchical model that adaptively shares statistical strength across classes when estimating these high-dimensional covariance matrices.

Generative classifiers like these have been largely eschewed in favor of simpler distance metrics, like the relative Mahalanobis distance score [@ren21rmds]. Here, we show both theoretically and empirically that RMDS is a close approximation to the outlier probability under a Gaussian DPMM with tied covariance matrices, connecting this widely-used approach to inference in a Bayesian nonparametric model. From this perspective, we propose hierarchical models that generalize RMDS by relaxing the assumption of equal covariance matrices across classes.

We find that hierarchical Gaussian DPMMs offer a well-grounded and practically competitive approach to OOD detection. Section 2 reviews related work, and Section 3 covers important background on DPMMs. We make a theoretical connection between RMDS and DPMMs in Section 4. This connection motivates our use of hierarchical models for estimating the high-dimensional covariance matrices in DPMMs — including a novel “coupled diagonal covariance” model — which we describe in 5 and evaluate in 6. We compare these models to existing baselines on synthetic datasets as well as the OpenOOD benchmark to characterize the regimes where Bayesian nonparameteric yield improved OOD performance ¹.

2 Related Work↩︎

The OOD detection task for image classification is the binary decision that predicts if an image belongs to a class that is outside of the set of known classes that have been observed during training. This task has been widely studied and many solutions have been proposed. For example, some approaches alter the architecture or objective of a classifier [@TackJ20; @HuangR21; @wei2022mitigating; @linderman23], and others exploit auxiliary outlier datasets [@HendrycksD19; @zhang2021mixture]. While training and outlier exposure based methods are important research directions for OOD detection the focus of this work is post-hoc methods. We refer readers to the available high quality OOD detection surveys to learn more about these methods [@yang2022openood; @zhang23openood15]. Our approach is related to a class of post-hoc methods including max softmax probability [@HendrycksD17], temperature-scaled MSP [@guo17tempscale], ODIN [@LiangS18], energy-based OOD [@LiuW20], the Mahalanobis distance score [@lee18mds], and the Relative MDS [@ren21rmds], which derive OOD scores from embeddings or activations of a pre-trained network.

Recently, [@zhang23openood15] proposed a set of Near and Far OOD benchmarks, as well as a leaderboard named OpenOOD to facilitate comparison across methods. The OpenOOD benchmarks found (1) that post-hoc methods are more scalable to large datasets, (2) there is no method that is best on all datasets, and (3) methods are sensitive to which model was used for embedding. The best performing OpenOOD methods for vision transformer (ViT) feature embeddings are the MDS and RMDS. The relative Mahalanobis distance score was inspired by earlier work by [@ren2019likelihood] that addressed the poor performance of the OOD performance with density estimation methods. [@sun2022out] propose to relax some of the assumptions of Mahalanobis distance methods by using the negative \(k\)-th nearest neighbor distances instead. We will show that the relative Mahalanobis distance score (RMDS) is similar to scores derived from Bayesian nonparametric mixture models in 4.

Bayesian nonparametric methods have previously been proposed for outlier detection and used in several applications. [@shotwell2011] proposed to detect outliers within datasets by partitioning data via a DPMM and identifying clusters containing a small number of samples as outliers. [@varadarajan17] developed a method for detecting anomalous activity in video snippets by modeling object motion with DPMMs. Another line of work explored Dirichlet prior networks [@malinin2018predictive; @malinin2019reverse] that explicitly model distributional uncertainty arising from dataset shift as a Dirichlet distribution over the categorical class probabilities. More recently, [@Kim2024UnsupervisedOD] performed unsupervised anomaly detection through an ensemble of Gaussian DPMMs fit to random projections of a subset of datapoints. Our work focuses on connecting DPMMs to post-hoc confidence scores and developing hierarchical Gaussian DPMMs that share statistical strength across classes in order to estimate their high-dimensional covariance matrices.

3 Background↩︎

We start with background on Dirichlet process mixture models (DPMMs) and the special case of a Gaussian DPMM with tied covariance.

3.1 Dirichlet process mixture models↩︎

Dirichlet process mixture models [@lo1984class] are Bayesian nonparametric models for clustering and density estimation that allow for a countably infinite number of clusters. There is always some probability that a new data point could come from a cluster that has never been seen before — i.e., that the new point is an outlier.

Let \(\mathcal{D}= \{x_n, y_n\}_{n=1}^N\) denote a set of training data points \(x_n \in \mathbb{R}^D\) and labels \(y_n \in [K]\). Likewise, let \(\mathcal{D}_k = \{x_n: y_n=k\}\) denote the subset of points assigned to cluster \(k\), and let \(N_k = |\mathcal{D}_k|\) denote the number of such points. Now consider a new, unlabeled data point \(x\). Under a DPMM, its corresponding label \(y\) has probability, \[\begin{align} \label{eq:outlier95prob} p(y = k ~\vert~x, \mathcal{D}) &\propto \begin{cases} N_k \, p(x ~\vert~\mathcal{D}_k) &\text{if } k\in[K] \\ \alpha \, p(x) &\text{if } k = K+1, \end{cases} \end{align}\tag{1}\] where the hyperparameter \(\alpha \in \mathbb{R}_+\) specifies the concentration of the Dirichlet process prior.

The first case captures the probability that the new point belongs to one of the training clusters (that it is an inlier). That probability depends on two factors: 1) the number of training data points in that cluster since, intuitively, larger clusters are more likely; 2) the posterior predictive probability, which is obtained by integrating over the posterior distribution of cluster parameters, \(\theta_k\), \[\begin{align} p(x ~\vert~\mathcal{D}_k) &= \int p(x ~\vert~\theta_k) \, p(\theta_k ~\vert~\mathcal{D}_k) \mathop{}\!\mathrm{d}\theta_k \\ \nonumber &\propto \int p(x ~\vert~\theta_k) \, \bigg[\prod_{x_n \in \mathcal{D}_k} p(x_n ~\vert~\theta_k) \bigg] \, p(\theta_k) \mathop{}\!\mathrm{d}\theta_k. \end{align}\] The second case of eq. 1 captures the probability that the new point is an outlier. It depends on the concentration \(\alpha\) and the prior predictive probability obtained by integrating over the prior distribution of cluster parameters, \[\begin{align} p(x) &= \int p(x ~\vert~\theta_k) \, p(\theta_k) \mathop{}\!\mathrm{d}\theta_k. \end{align}\] For many models of interest, the posterior and prior predictive distributions have closed forms.

3.2 Gaussian DPMM with Tied Covariance↩︎

For example, consider a Gaussian DPMM in which each cluster is parameterized by a mean and covariance, \({\theta_k = (\mu_k, \Sigma_k)}\). Assume a conjugate prior for the mean, \(\mu_k \sim \mathcal{N}(\mu_0, \Sigma_0)\). For now, assume that all clusters share the same covariance matrix, which we express through an atomic prior, \(\Sigma_k \sim \delta_\Sigma\). The hyperparameters of the prior are \(\eta = (\mu_0, \Sigma_0, \Sigma)\).

Under this Gaussian DPMM, the conditional distribution of a new data point’s label is, \[\begin{align} \label{eq:gaussian95dpmm95outlier95unnormalized} p(y = k ~\vert~x, \mathcal{D}) &\!\propto\! \begin{cases} N_k \, \mathcal{N}(x ~\vert~\mu_k', \Sigma_k' + \Sigma) &\!\!\!\text{if } k\in[K] \\ \alpha \, \mathcal{N}(x ~\vert~\mu_0, \Sigma_0 + \Sigma) &\!\!\!\text{if } k = K+1. \end{cases} \end{align}\tag{2}\] where \[\begin{align} \nonumber \mu_k' &= \Sigma_k' \left(\Sigma_0^{-1} \mu_0 + N_k \Sigma^{-1} \bar{x}_k\right), \\ \label{eq:gaussian95tied95postpred} \Sigma_k' &= \left( \Sigma_0^{-1} + N_k \Sigma^{-1} \right)^{-1}, \end{align}\tag{3}\] and \(\bar{x}_k = \tfrac{1}{N_k} \sum_{x_n \in \mathcal{D}_k} x_n\) is the mean of the data points assigned to cluster \(k\).

The relative probability of theses cases is an intuitive measure of how likely a point is to be an outlier. Indeed, the next section shows that the outlier probabilities from this Bayesian nonparametric model are closely related to another common outlier detection score.

4 Theory: Connecting Relative Mahalanobis Distance and DPMMs↩︎

Here we show that a widely used outlier detection method called the relative Mahalanobis distance score [@ren21rmds] is closely related to the outlier probabilities obtained using a Gaussian DPMM with tied covariances. RMDS outputs a score, \(C(x)\), where smaller values indicate that a data point \(x\) is more likely to be an outlier. The RMDS score of a new point \(x\) is defined as follows,² \[\begin{align} \nonumber \mathrm{MD}_0(x) &= (x - \hat{\mu}_0)^\top \hat{\Sigma}_0^{-1} (x - \hat{\mu}_0) \\ \nonumber \mathrm{MD}_k(x) &= (x - \hat{\mu}_k)^\top \hat{\Sigma}^{-1} (x - \hat{\mu}_k) \\ \nonumber \mathrm{RMD}_k(x) &= \mathrm{MD}_0(x) - \mathrm{MD}_k(x)\\ C(x) &= \max_k \; \mathrm{RMD}_k(x), \end{align}\] where \(\mathrm{MD}_0(x)\) and \(\mathrm{MD}_k(x)\) are squared Mahalanobis distances, \(\hat{\mu}_0\) and \(\hat{\Sigma}_0\) are the sample mean and covariance of the data, \(\hat{\mu}_k\) is the sample mean of cluster \(k\), and \(\hat{\Sigma} = \tfrac{1}{N} \sum_{n=1}^N (x_n - \hat{\mu}_{y_n})(x_n - \hat{\mu}_{y_n})^\top\) is the sample within-class covariance.

[@ren21rmds] motivated this score in terms of log density ratios between a Gaussian distribution for each cluster and a Gaussian “background” model. Specifically, \[\begin{align} \label{eq:rmdk} \mathrm{RMD}_k(x) &= 2 \log \frac{\mathcal{N}(x ~\vert~\hat{\mu}_k, \hat{\Sigma})}{\mathcal{N}(x ~\vert~\hat{\mu}_0, \hat{\Sigma}_0)} + d, \end{align}\tag{4}\] where \(d=\log |\hat{\Sigma}| - \log|\hat{\Sigma}_0|\) does not depend on \(x\) or \(k\). Larger values of \(\mathrm{RMD}_k(x)\) indicate that \(x\) is more likely under cluster \(k\) than under the background model.

The procedure for mapping \(\mathrm{RMD}_k(x)\) values to the score \(C(x)\) is inherited from the Mahalanobis distance score [@lee18mds]. If the log density ratio is negative for all \(k\), then the background model is more likely than all of the existing clusters, and hence \(x\) is likely to be an outlier. Propositions 1 and 2 show that a similar computation is at work in the outlier probabilities for DPMMs.

First, we show that the inlier probabilities under a general DPMM (not necessarily Gaussian) can be expressed in terms of a quantity analogous to \(C(x)\).

This proposition says that the log-odds of data point \(x\) belonging to an existing cluster is the difference of a DPMM score, \(\widetilde{C}(x)\), which is analogous to the relative Mahalanobis score, and a threshold, \(\log \frac{\alpha}{\overline{N}}\), which is tuned by the hyperparameter \(\alpha\).

Next, we show that in certain regimes, the DPMM score from a Gaussian DPMM with tied covariance is almost perfectly correlated with the RMDS. Below, we define the relative covariance matrix \(R = \hat{\Sigma}_0^{-1/2}\hat{\Sigma}\hat{\Sigma}_0^{-1/2}\) and note that as its operator norm \(\kappa = \|R\|_{op}\) goes to \(0\), then we intuitively have that \(\hat{\Sigma}_0\) is growing larger with respect to \(\hat{\Sigma}\).

This proposition establishes the close correspondence between the relative Mahalanobis distance score and the log-odds that a point is an inlier under a Gaussian DPMM with tied covariance. Note that the \(\log K\) factor in 2 is irreducible due to the difference between the \(\max\) used in RMDS and the smooth approximation used in DPMMs. We view this as a feature not a bug of the DPMMs: Since the \(\log K\) factor only appears when the scores for different clusters are close to identical, the gap arises when the DPMM aggregates evidence across multiple equally plausible components rather than arbitrarily selecting a single cluster with RMDS. In practice, we find that a close correspondence holds in the experiments below and in 16. This correspondence provides further support for using RMDS for outlier detection, beyond the original motivation in terms of log likelihood ratios. However, from this perspective, we also recognize several natural generalizations of RMDS that could improve outlier detection through richer DPMMs. We present three such generalizations below.

5 Hierarchical Gaussian DPMMs↩︎

RMDS has proven to be a highly effective outlier detection method, but it assumes that all clusters share the same covariance. This assumption helps avoid overfitting the covariance matrices for each class [@ren21rmds], but it is not always warranted. 1 shows a histogram of differences between empirical covariance matrices \(\hat{\Sigma}_k\) and \(\hat{\Sigma}_{k'}\) for all pairs of classes \((k,k')\) in the Imagenet-1K dataset, as measured by the Förstner-Moonen distance [@forstner_metric_2003].

These pairwise distances are systematically larger than what we would expect under a null distribution where the true covariance matrices are the same for all classes, and the empirical estimates differ solely due to sampling variability. Complete details of this analysis are provided in 9. This analysis suggests that the covariance matrices are significantly different across classes and motivates a more flexible approach.

The connection between RMDS and Gaussian DPMMs established above suggests a natural way of relaxing the tied-covariance assumption without sacrificing statistical power: Instead of estimating covariance matrices independently, we could infer them jointly under a hierarchical Bayesian model [@gelman1995bayesian]. With such a model, we can estimate separate covariance matrices for each cluster, while sharing information via a prior. By tuning the strength of the prior, we can obtain the tied covariance model in one limit and a fully independent model in the other. Finally, we can estimate these hierarchical prior parameters using a simple expectation-maximization algorithm that runs in a matter of minutes, even with large, high-dimensional datasets.

5.1 Full Covariance Model↩︎

First, we propose a hierarchical Gaussian DPMM with full covariance matrices and a conjugate prior. The cluster parameters, \(\theta_k = (\mu_k, \Sigma_k)\), are drawn from a conjugate, normal-inverse Wishart (NIW) prior, \[p(\theta_k) = \mathrm{IW}\big(\Sigma_k ~\vert~\nu_0, (\nu_0 - D - 1) \Sigma_0 \big) \times \mathcal{N}\big(\mu_k ~\vert~\mu_0, \kappa_0^{-1} \Sigma_k \big),\] where \(\mathrm{IW}\) denotes the inverse Wishart density. Under this parameterization, \(\mathbb{E}[\Sigma_k] = \Sigma_0\) for \(\nu_0 > D + 1\). The hyperparameters of the prior are \(\eta = (\nu_0, \kappa_0, \mu_0, \Sigma_0)\).

The most important hyperparameters are \(\nu_0\) and \(\Sigma_0\), as they specify the prior on covariance matrices. As \(\nu_0 \to \infty\), the prior concentrates around its mean and we recover a tied covariance model. For small values of \(\nu_0\), the hierarchical model shares little strength across clusters, and the covariance estimates are effectively independent.

We propose a simple approach to estimate these hyperparameters in 13. Briefly, we use empirical Bayes estimates for the prior mean and covariance, setting \(\mu_0 = \hat{\mu}_0\) and \(\Sigma_0 = \hat{\Sigma}\). We derive an expectation-maximization (EM) algorithm to optimize \(\nu_0\) and \(\kappa_0\). Thanks to the conjugacy of the model, the E-step and the M-step for \(\kappa_0\) can be computed in closed form. We leverage a generalized Newton method [@minka2000beyond] to update the concentration hyperparameter, \(\nu_0\), effectively learning the strength of the prior to maximize the marginal likelihood of the data.

Finally, the prior and posterior predictive distributions are multivariate Student’s t distributions with closed-form densities. The log density ratios derived from these predictive distributions form the basis of the DPMM scores, \(\widetilde{C}(x)\).

5.2 Diagonal Covariance Model↩︎

Even with the hierarchical prior, we find that the full covariance model can still overfit to high-dimensional embeddings. Thus, we also consider a simplified version of the hierarchical Gaussian DPMM with diagonal covariance matrices. Here, the cluser parameters are \(\theta_k = \{\mu_{k,d}, \sigma_{k,d}^2\}_{d=1}^D\), and the conjugate prior is, \[p(\theta_k) = \prod_{d=1}^D \chi^{-2}(\sigma_{k,d}^2 ~\vert~\nu_{0,d}, \sigma_{0,d}^2) \times \mathcal{N}(\mu_{k,d} ~\vert~\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2)\] where \(\chi^{-2}\) is the scaled inverse chi-squared density.

In addition to having fewer parameters per cluster, another advantage of this model is that it allows for different concentration hyperparameters for each dimension, \(\nu_{0,d}\). We estimate the hyperparameters using a procedure that closely parallels the full covariance model. Likewise, the prior and posterior predictive densities reduce to products of scalar Student’s t densities, which are even more efficient to compute. Complete details are in 14.

5.3 Coupled Diagonal Covariance Model↩︎

The diagonal covariance model dramatically reduces the number of parameters per cluster, but it also makes a strong assumption about the per-class covariance matrices. Specifically, it assumes the variances, \(\sigma_{k,d}^2\), are conditionally independent across dimensions. 2 suggests that this is not the case: the diagonals of the empirical covariance matrices, \(\hat{\Sigma}_k\), tend to be systematically larger or smaller than those of the average covariance matrix, \(\hat{\Sigma}\). This analysis suggests that \(\sigma_{k,d}^2\) are not independent; rather, if \(\sigma_{k,d}^2\) is larger than average, then \(\sigma_{k,d'}^2\) is likely to be larger as well.

We propose a novel, coupled diagonal covariance model to capture these effects. Specifically, we introduce a scale factor \(\gamma_k \in \mathbb{R}_+\) that scales the variances for class \(k\) compared to the average. In this model, the cluster parameters are \(\theta_k = (\gamma_k, \{\mu_{k,d}, \sigma_{k,d}^2\}_{d=1}^D)\), and the prior is, \[p(\theta_k) = \chi^2(\gamma_k ~\vert~\alpha_0) \prod_{d=1}^D \bigg[\chi^{-2}(\sigma_{k,d}^2 ~\vert~\nu_{0,d}, \gamma_k \sigma_{0,d}^2) \times \mathcal{N}(\mu_{k,d} ~\vert~\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2) \bigg]\] where \(\gamma_k\) scales the means of \(\sigma_{k,d}^2\) for all dimensions \(d\) to capture the correlations seen above.

Our procedure for hyperparameter estimation and computing DPMM scores is very similar to those described above. The only complication is that with the \(\gamma_k\), the posterior distribution no longer has a simple closed form. However, for any fixed value of \(\gamma_k\), the coupled model is a straightforward generalization of the diagonal model above. Since \(\gamma_k\) is a one-dimensional variable, we can use numerical quadrature to integrate over its possible values. Likewise, we can estimate the hyperparameter \(\alpha_0\) using a generalized Newton method, just like for the concentration parameter \(\nu_0\). See 15 for complete details of this model.

6 Experiments↩︎

We experimentally tested these hierarchical Gaussian DPMMs on real and synthetic datasets. First we used simulated datasets to build intuition for where hierarchical models improve performance. Then we compared hierarchical Gaussian DPMMs to other widely-used OOD metrics on the OpenOOD Benchmark, and we studied performance versus dimensionality of the embeddings.

6.1 Synthetically generated dataset experiments↩︎

To understand the regimes in which hierarchical models yield benefits, we simulated \(D=2\) dimensional data from full covariance models with varying \(\nu_0\) and \(N_k\). When \(\nu_0\) is small, covariances differ considerably across clusters, and the assumptions of the tied model (and of RMDS) are not well met. Conversely, as \(\nu_0 \to \infty\), the prior concentrates on \(\Sigma_0\), and the covariances are effectively tied. 3A and 3B show simulated datasets from these two regimes. As expected, 3C shows that hierarchical Gaussian DPMMs yielded considerable improvements when \(\nu_0\) was small, and the largest improvements came from the full covariance model, which matched the data generating process. Notably, the tied model matched the RMDS performance, as predicted in 4.

We then asked if the hierarchical model was strictly necessary or whether a simpler model would suffice. For example, we considered an “Independent RMDS” based on Mahalanobis distances to the per-class covariance estimates, \(\hat{\Sigma}_k\), instead of the average covariance \(\hat{\Sigma}\). Intuitively, we expected the hierarchical models to perform best when there were few data points per cluster relative to the dimensionality of the embeddings; i.e., when \(N_k/D\) is small. Indeed, 3D shows that DPMMs offered substantial improvements in this regime, with diminishing gains as \(N_k\) increased.

These analyses suggest that hierarchical Gaussian DPMMs should yield benefits in regimes where covariances differ across classes and the number of data points per class is relatively small.

6.2 OpenOOD Benchmark↩︎

Next, we compared hierarchical Gaussian DPMMs to other widely-used OOD detection methods on the OpenOOD benchmark datasets [@yang2022openood; @zhang23openood15]. The benchmark consists of 3 in-distribution (ID) datasets, CIFAR-10, CIFAR-100, and Imagenet-1K. For each ID dataset, several OOD datasets are grouped into Near and Far, where the Near OOD datasets are more similar to ID. For the Imagenet-1K experiment, we used \(D=768\) dimensional embeddings from the ViT-B-16 model trained according to the DeIT method [@touvron2021training], which are available in the Pytorch torchvision [@torchvision2016] package. The CIFAR experiments use OpenOOD’s pretrained ResNet18 [@HeK16] features, \(D=512\). We preprocessed the embeddings as described in 8. As baselines, we considered both MDS [@lee18mds] and RMDS [@ren21rmds]. We also compare to maximum softmax probability (MSP) [@HendrycksD17] and temperature scaled MSP with \(T=1000\) (Temp. Scale) [@guo17tempscale] . We trained a single linear layer with gradient descent and supervised cross-entropy loss for the MSP methods. For all models, we measured the accuracy of classifying which class an in-distribution test image came from, as well as the AUROC score for outlier detection across the OpenOOD datasets. For DPMMs, we computed AUROC scores using \(\widetilde{C}(x)\).

Table 1 shows that hierarchical models do offer improved performance, but the results are nuanced. We found the complexity and scale of the ID dataset determines which modeling assumptions are most appropriate. The full covariance model yields the highest performance on Near and Far OOD on the small scale CIFAR-10 task. Whereas the coupled diagonal model is the best performing model in Near OOD settings for the large-scale Imagenet-1K task. The full covariance model performed surprisingly poorly on the Imagenet-1K task, and we suspected it was due to the high-dimensional embeddings. We provide expanded results tables for each ID dataset that show performance on each OOD dataset in [app:ablation] [app:cifar10-ablation] [app:cifar100-ablation].

6.3 Performance vs. Dimensionality↩︎

Finally, we asked how these methods compare as we vary the dimensionality of the embeddings. We suspected that the full covariance model would perform better in lower dimensions for two reasons. First, it has \(\mathcal{O}(KD^2)\) parameters compared to only \(\mathcal{O}(KD)\) in the diagonal models and \(\mathcal{O}(D^2)\) in the tied model, so even with a hierarchical prior, the full model could still overfit. This problem is exacerbated for classes that have fewer data points than the feature dimension, in which case the prior has a strong effect on the conditional distribution of the per-class covariance matrix and the posterior predictive distributions. Second, we suspected that the inverse Wishart prior distribution, which has only a scalar concentration \(\nu_0\), may be a poor prior for high-dimensional covariance matrices.

To test this hypothesis, we swept the number of principal components retained in preprocessing (see 8). We found that out-of-distribution detection of the full-covariance hierarchical model plateaued for \(D \geq 128\) dimensional embeddings (4). By contrast, the diagonal and coupled diagonal models performed considerably better, especially on Near OOD benchmarks. The diagonal covariance model outperforms the tied model across all dimensions on Near OOD detection, as well as in lower dimensions for Far OOD. However, in-distribution classification accuracy plateaus around 256 dimensions.

Altogether, these analyses of synthetic and real datasets show that hierarchical Gaussian DPMMs are advantageous for OOD detection, especially in regimes where: i) covariance matrices differ across clusters; ii) the number of data points per cluster is small compared to the dimension; and iii) detection relies on fine-grained distinctions between training data and Near OOD test points.

7 Discussion↩︎

We developed a theoretical connection between the relative Mahalanobis distance score for outlier detection and the outlier probability under a Gaussian DPMM with tied covariance. This Bayesian nonparametric perpsective led us to propose hierarchical Gaussian DPMMs that allow each cluster a different covariance matrix, while still sharing statistical strength across classes via the prior. We developed efficient EM algorithms to estimate the hyperparameters of the hierarchical models, and we studied their performance on synthetic data as well as the OpenOOD benchmarks. We found that these models — especially the coupled diagonal covariance model — yielded improved performance on some benchmarks, especially the Near OOD benchmarks.

7.0.0.1 Limitations and Future Work

Despite the hierarchical prior, we found that the full covariance Gaussian DPMMs were prone to overfitting. Further work could explore low-rank plus diagonal covariance matrices, which would interpolate between the diagonal and full covariance models. More generally, like RMDS and MDS, we assume that features are Gaussian distributed within each class. The competitive performance on the OpenOOD Imagenet benchmark using ViT features suggests that this assumption is reasonable [@yang2022openood; @zhang23openood15], but there is no guarantee. Future work could consider fine-tuning the features or learning a nonlinear transformation to address this potential source of model misspecification, as in prototype networks [@snell2017prototypical].

Bayesian nonparametric approaches naturally extend to other closely related problems, like generalized category discovery [@vaze22] and continual learning [@van2022three]. For example, given a collection of data points, a DPMM may have sufficient evidence to allocate new classes for the out-of-distribution data. More generally, casting OOD as inference in a generative model brings modeling choices to the fore. Here, we focused on the challenge of modeling high-dimensional covariance matrices that may vary across classes, but there are several other ways in which the simple Gaussian DPMM could be improved. For example, we could attempt to capture the nonstationarity inherent in the OOD setting by allowing the prior predictive distribution to drift from what was inferred based on the training data. Such a model could afford greater robustness on OOD detection tasks.

7.0.0.2 Conclusion

In summary, we find that Bayesian nonparametric methods with hierarchical priors are a promising approach for OOD detection. If the features extracted from foundation models are reasonably well approximated as realizations of Gaussian DPMMs, the posterior inferences under such models can provide accurate estimates of outlier probability. This probabilistic perspective not only casts widely used methods in a new light, it also leads to practical model improvements and enables several lines of future inquiry.

Acknowledgments↩︎

We thank Jingyang Zhang for his feedback and detailed knowledge of the OpenOOD codebase. RWL was supported in part by the National Science Foundation (NSF) under Grant No. 2112562 to YC, the U.S. Army Research Laboratory and the U.S. Army Research Office (ARL/ARO) under grant number ARO-W911NF-23-2-0224 to YC, and the Department of Defense (DoD) through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program. SWL was supported in part by fellowships from the Simons Collaboration on the Global Brain, the Alfred P. Sloan Foundation, and the McKnight Foundation.

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF, the ARL/ARO, or the DoD.

8 Preprocessing↩︎

Before computing OOD scores, we first preprocessed the embeddings using PCA to whiten and sort the dimensions in order of decreasing variance. We discarded dimensions with near zero variance to ensure the empirical covariance matrices were full rank. We scaled each dimension by the inverse square root of the eigenvalues so that the transformed embeddings had identity covariance. Finally, we rotated the embeddings using the eigenvectors of the average covariance matrix, so that the average within-class covariance matrix is diagonal.

More precisely, the preprocessing steps are as follows:

1. Let \(\{x_i\}_{i=1}^N\) denote the mean-centered embeddings.

2. Let \(\hat{\Sigma}_0 = U \Lambda U^\top\) denote the covariance of the centered embeddings and its eigendecomposition. Discard any dimensions with eigenvalues less than a threshold of approximately \(10^{-7}\), and then we project and scale the embeddings by, \[\begin{align} x_i' &\leftarrow \Lambda^{-\frac{1}{2}} U^T x_i, \end{align}\] so that the empirical covariance of \(\{x_i'\}_{i=1}^n\) is the identity matrix.

3. Compute the average within-class covariance \(\hat{\Sigma} = \frac{1}{N} \sum_{i=1}^M (x_i' - \hat{\mu}_{y_i})(x_i' - \hat{\mu}_{y_i})^\top\), where \(\hat{\mu}_k = \frac{1}{N_k} \sum_{i:y_i=k} x_i'\).

4. Compute the eigendecomposition \(\hat{\Sigma} = V S V^\top\), with eigenvalues \(S = \mathrm{diag}(\sigma_1^2, \ldots, \sigma_D^2)\) sorted in increasing order of magnitude so that the first dimension has the smallest within-class covariance. The embeddings \(x_i'\) have unit variance in all dimensions, but along the dimension of the first eigenvector in \(V\), the average within-class covariance is smallest.

5. Project the embeddings into this eigenbasis, \[\begin{align} z_i &\leftarrow V^\top x_i'. \end{align}\]

After these preprocessing steps, the resulting embeddings \(\{z_i\}_{i=1}^N\) are zero mean (\(\hat{\mu}_0=0\)), their empirical covariance is the identity (\(\hat{\Sigma}_0=I\)), and the average within-class covariance is diagonal (\(\hat{\Sigma}=\mathrm{diag}(\sigma_1^2, \ldots, \sigma_D^2)\)). The empirical within-class covariance matrix \(\hat{\Sigma}_k\) for class \(k\) will not generally be diagonal, but this sequence of preprocessing steps is intended to make them closer to diagonal on average.

Note that the relative Mahalanobis distance score is invariant to these linear transformations. They simply render the embeddings more amenable to our hierarchical models with diagonal covariance. We further test the effect of these preprocessing steps via the ablation experiment described in 17.

9 Further Details of Exploratory Analyses↩︎

First, we investigated the degree to which the sample covariance matrices differ between the 1000 classes in the Imagenet-1K dataset. In accordance with the OpenOOD benchmark, we used embeddings from the ViT-B-16 model trained according to the DeIT method [@touvron2021training], which are \(D=768\) dimensional. The embeddings are available in the Pytorch torchvision [@torchvision2016] package. We preprocessed the embeddings as described in 8. For this analysis, we only kept the top 128 PCs to speed computation.

To measure the distance between covariance matrices for two clusters, we used the Förstner-Moonen (FM) metric [@forstner_metric_2003], \[\begin{align} d(\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2) = \bigg[\sum_{i=1}^n (\log \lambda_i (\boldsymbol{\Sigma}_1^{-1} \boldsymbol{\Sigma}_2))^{2}\bigg]^{\frac{1}{2}}, \end{align}\] where \(\lambda_i(\boldsymbol{\Sigma}_1^{-1} \boldsymbol{\Sigma}_2)\) is the \(i\)-th eigenvalue of \(\boldsymbol{\Sigma}_1^{-1} \boldsymbol{\Sigma}_2\). We computed the FM metric for all pair of empirical covariance matrices \((\hat{\Sigma}_k, \hat{\Sigma}_{k'})\). We compared the distribution of FM distances under the real data to distances between covariance matrices sampled from the null model, in which \(\Sigma_k\) truly equals \(\hat{\Sigma}\) for all \(k\), and the differences in the estimates \(\hat{\Sigma}_k\) arise solely from sampling error. The corresponding null distribution is a Wishart distribution, \(\hat{\Sigma}_k \sim \mathrm{W}(\overline{N}, \hat{\Sigma}/\overline{N})\), where \(\overline{N}\) is the average number of data points per class.

10 Compute Resources↩︎

The OOD experiments were performed on a cluster consisting of compute nodes with 8 NVIDIA RTX A5000 GPUs. The OpenOOD [@yang2022openood; @zhang23openood15] experiments on the Imagenet-1K dataset [@RussakovskyO15] utilized weights available from Pytorch’s [@PYTORCH] torchvision [@torchvision2016] models. The compute across experiments was reduced by storing summary statistics of the embeddings for the Gaussian models. The DPMM fitting and prediction was performed on the CPU except for calculating the posterior and prior predictive distributions for each sample.

11 Computational Complexity↩︎

The complexity of computing outlier scores is not that much greater than computing the relative Mahalanobis distance score (RMDS). Our preprocessing procedure (8) involves whitening and rotating the embeddings, which requires \(\mathcal{O}(D^2)\) memory and \(\mathcal{O}(D^3)\) time, just like RMDS.

In order to account for different covariances across classes, we incur additional costs. The first costs are those of estimating the hyperparameters via expectation-maximization (EM). In the EM algorithm, we need to store expected covariance matrices (or rather, their inverses) for each class. This costs \(O(KD^2)\) space for the full covariance model and \(\mathcal{O}(KD)\) space for the (coupled) diagonal model. Each M-step of the EM algorithm also involves inverting covariance matrices, so for the full covariance model the time complexity is \(\mathcal{O}(TKD^3)\), and for the diagonal and coupled-diagonal covariance model it is only \(\mathcal{O}(TKD)\), where \(T\) is the number of EM iterations. Typically, we find that \(T \approx 10\) iterations is enough to converge.

Finally, there is the cost of computing the final outlier scores, which involves computing Mahalanobis distances for RMDS and a marginal likelihood ratio the DPMM. For both RMDS and full covariance DPMM scores, these can be computed in \(\mathcal{O}(ND^2)\) time using the Cholesky factors of the covariance matrices from the preprocessing or EM algorithm; for the (coupled) diagonal covariance model, this cost is only \(\mathcal{O}(ND)\).

In sum, the additional expressivity of the Bayesian nonparametric mixture models incurs some additional cost. For the full covariance model, these costs can become burdensome (especially in memory), but for the diagonal and coupled diagonal model they are quite reasonable and, as we show in our experiments, they can lend a significant boost.

12 Proof of Proposition 2↩︎

Using the definitions, we can write \(\lambda_k\) exactly as a difference of log Gaussian densities:

\[\lambda_k \;=\; \log \frac{\mathcal{N}(x ~\vert~\mu_k',\,\Sigma_k'+\hat{\Sigma})}{\mathcal{N}(x ~\vert~\hat{\mu}_0,\,\hat{\Sigma}_0+\hat{\Sigma})} \;=\; \Big[\log \mathcal{N}(x ~\vert~\mu_k',\,\Sigma_k'+\hat{\Sigma})\;-\;\log \mathcal{N}(x ~\vert~\hat{\mu}_0,\,\hat{\Sigma}_0+\hat{\Sigma})\Big]\,.\]

Similarly, the RMDS for cluster \(k\) can be written as a log-density ratio between a cluster Gaussian and the “background” Gaussian model:

\[\mathrm{RMD}_k(x) \;=\; 2\log\frac{\mathcal{N}(x ~\vert~\hat{\mu}_k,\,\hat{\Sigma})}{\mathcal{N}(x ~\vert~\hat{\mu}_0,\,\hat{\Sigma}_0)} + d,\]

where \(d=\log|\hat{\Sigma}|-\log|\hat{\Sigma}_0|\). For clarity, define the ideal log-ratio using true cluster parameters as

\[\rho_k(x) \;:=\; \log\frac{\mathcal{N}(x ~\vert~\mu_k,\,\Sigma)}{\mathcal{N}(x ~\vert~\mu_0,\,\Sigma_0)} \,.\]

Notice that \(\frac{1}{2}[\mathrm{RMD}_k(x)-d] = \log\frac{\mathcal{N}(x ~\vert~\hat{\mu}_k,\,\hat{\Sigma})}{\mathcal{N}(x ~\vert~\hat{\mu}_0,\,\hat{\Sigma}_0)} =: \hat{\rho}_k(x)\) is the same ratio but with empirical estimates \((\hat{\mu}_k,\hat{\Sigma})\) in place of \((\mu_k,\Sigma)\), and \(\rho_k(x)\) uses the true cluster mean and covariance. We wish to bound \(2|\lambda_k - \hat{\rho}_k|\).
We now can use the triangle inequality and decompose this into four separate terms: \[\begin{align} 2|\lambda_k - \hat{\rho}_k| &\leq \left|\log|\hat{\Sigma}_0 + \hat{\Sigma}| - \log{|\hat{\Sigma}_0|} \right|\bigg\} & & & (\mathsf{A}) \\ &\quad \quad + \left|(x - \hat{\mu}_0)^T\left[(\hat{\Sigma}_0 + \hat{\Sigma})^{-1} - \hat{\Sigma}_0^{-1}\right](x - \hat{\mu}_0)\right|\bigg\} & & & (\mathsf B) \\ & \quad \quad + \left|\log{|\Sigma_k' + \hat{\Sigma}|} - \log{|\hat{\Sigma}|}\right|\bigg\}& & & (\mathsf C) \\ &\quad \quad +\left|(x - \hat{\mu}_k)^T\hat{\Sigma}^{-1}(x - \hat{\mu}_k)- (x - \mu_k')^T(\Sigma_k + \hat{\Sigma})^{-1}(x - \mu_k')\right|\bigg\}& & & (\mathsf D) \end{align}\] To bound \(\mathsf{A}\), we first note that \[\log|\hat{\Sigma}_0 + \hat{\Sigma}| - \log{|\hat{\Sigma}_0|} = \log|I + \hat{\Sigma}_0^{-1}\hat{\Sigma}| = \log |I + \hat{\Sigma}_0^{-1/2}\hat{\Sigma}\hat{\Sigma}_0^{-1/2}| = \log |I + R|\] using the fact that the determinant behaves multiplicatively with respect to the matrix product. Then given that the determinant is the product of the eigenvalues we finish our bound on \(\mathsf{A}\) \[\mathsf{A} \leq \sum_{i=1}^D\log(1 + \lambda_i^{(R)}) \leq D\log(1 + \kappa) \leq D\kappa\] thus this term is less than \(\delta_1\) when \(\kappa < \kappa_1 \leq \frac{\delta_1}{D}\). For term \(\mathsf{B}\), recall that \[(A + B)^{-1} = A^{-1} - A^{-1}B(A + B)^{-1}\] (which can be verified by right multiplying by \(A + B\)) combined with the fact that \[\begin{align} \hat{\Sigma}_0^{1/2}\hat{\Sigma}_0^{-1}\hat{\Sigma}(\hat{\Sigma}_0 + \hat{\Sigma})^{-1}\hat{\Sigma}_0^{1/2} &= \hat{\Sigma}_0^{-1/2}\hat{\Sigma}(\hat{\Sigma}_0 + \hat{\Sigma})^{-1}\hat{\Sigma}_0^{1/2} \\ &= \hat{\Sigma}_0^{-1/2}\hat{\Sigma}(\hat{\Sigma}_0(I + \hat{\Sigma}_0^{-1}\hat{\Sigma}))^{-1}\hat{\Sigma}_0^{1/2} \\ &= \hat{\Sigma}_0^{-1/2}\hat{\Sigma}(I + \hat{\Sigma}_0^{-1}\hat{\Sigma})^{-1}\hat{\Sigma}_0^{-1}\hat{\Sigma}_0^{1/2} \\ &= \hat{\Sigma}_0^{-1/2}\hat{\Sigma}(I + \hat{\Sigma}_0^{-1}\hat{\Sigma})^{-1}\hat{\Sigma}_0^{-1/2} \\ &= \hat{\Sigma}_0^{-1/2}\hat{\Sigma}[\hat{\Sigma}_0^{-1/2}(I + R)\hat{\Sigma}_0^{1/2}]^{-1}\hat{\Sigma}_0^{-1/2}\\ &= R[1 + R]^{-1} \end{align}\] Means that we can say \[\begin{align} \hat{\Sigma}_0^{-1}\hat{\Sigma}(\hat{\Sigma}_0 + \hat{\Sigma})^{-1} &\preceq \frac{\kappa}{1 + \kappa}\hat{\Sigma}_0^{-1} \end{align}\] By the standard fact that for \(B\) positive definite, \(A^{-1/2}BA^{-1/2} \preceq \alpha I\) if and only if \(B \preceq \alpha A\). So combining all of these identities we have \[\mathsf{B} \leq \frac{\kappa}{1 + \kappa}(x - \hat{\mu}_0)^\top \hat{\Sigma}_0^{-1}(x - \hat{\mu}_0)\] Which we then bound with the standard Hanson Wright inequality [@vershynin2018high] to get that with probability \(1 - \epsilon_1\) \[\begin{align} \mathsf{B} \leq \kappa \sigma^2 (D + \sqrt{D\log\frac{1}{\epsilon_1}} + 2\log \frac{1}{\epsilon_1}) \end{align}\] Which is less than \(\delta_2\) when \(\kappa < \kappa_2 = \left[\sigma^2 (D + \sqrt{D\log\frac{1}{\epsilon_1}} + 2\log \frac{1}{\epsilon_1})\right]^{-1}\delta_2\)

Now to deal with term \(\mathsf{C}\), we again note the multiplicativity of the determinant with respect to matrix products and see that \[\log |\Sigma_k' + \hat{\Sigma}| - \log |\hat{\Sigma}| = \log | I + \hat{\Sigma}^{-1}\Sigma_k'|\]

But recall that \(\Sigma_k' = (\hat{\Sigma}_0^{-1} + N\hat{\Sigma}^{-1})^{-1}\), which immediately gives a useful result when combined with the fact that \(\hat{\Sigma}_0^{-1}\hat{\Sigma} = \hat{\Sigma}_0^{-1/2}R\hat{\Sigma}_0^{1/2}\) \[\begin{align} \log | I + \hat{\Sigma}^{-1}\Sigma_k'| &= \log | I + \hat{\Sigma}^{-1}(\hat{\Sigma}_0^{-1} + N\hat{\Sigma}^{-1})^{-1}| \\ &= \log |I + (NI + \hat{\Sigma}_0^{-1}\hat{\Sigma})^{-1}| \\ &= \log |I + (NI + \hat{\Sigma}_0^{-1/2}R\hat{\Sigma}_0^{1/2})^{-1}| \\ &= \log |\hat{\Sigma}_0^{1/2}\hat{\Sigma}_0^{-1/2} + (N\hat{\Sigma}_0^{-1/2}\hat{\Sigma}_0^{1/2} + \hat{\Sigma}_0^{-1/2}R\hat{\Sigma}_0^{1/2})^{-1}| \\ &= \log |\hat{\Sigma}_0^{1/2}\hat{\Sigma}_0^{-1/2} + \hat{\Sigma}_0^{1/2}(NI + R)^{-1}\hat{\Sigma}_0^{-1/2}| \\ &= \log |\hat{\Sigma}_0^{1/2}||I + (NI + R)^{-1}||\hat{\Sigma}_0^{-1/2}| \\ &= \log |I + (NI + R)^{-1}| \end{align}\] Now we are functionally done, since we can bound this easily using known bounds of the logarithm (easy to check by concavity) to say \[\begin{align} \log |I + (NI + R)^{-1}| &= \sum_{i=1}^D\log \left(1 + \frac{1}{N + \lambda_i}\right) \end{align}\] \[\mathsf{C} \leq D \log\left(1 + \frac{1}{N}\right) \leq \frac{D}{N}\] So when \(N \geq N_1 > \frac{D}{\delta_3}\) then \(\mathsf{C} \leq \delta_3\).

Finally, we turn to \(\mathsf{D}\). We note that we can re-write it as \[(x - \mu_k' + \mu_k' - \hat{\mu}_k)^T\hat{\Sigma}^{-1}(x - \mu_k' + \mu_k' - \hat{\mu}_k) - (x - \mu_k')^T(\Sigma_k + \hat{\Sigma})^{-1}(x - \mu_k')\] which decomposes into three parts. Trying to ease notation, let \(v = x - \mu_k'\) and \(u = \mu_k' - \hat{\mu}_k\). So then \(\mathsf{D}\) is just \[\underbrace{v^\top[\hat{\Sigma}^{-1} - (\Sigma_k + \hat{\Sigma})^{-1}]v}_{(i)} + \underbrace{u^\top\hat{\Sigma}^{-1}u}_{(ii)} + \underbrace{2 u^\top\hat{\Sigma}^{-1}v}_{(iii)}\] For the first piece, we again use the identity \[(A + B)^{-1} = A^{-1} - A^{-1}B(A + B)^{-1}\] to get that \((i)\) is equal to \[v^\top \hat{\Sigma}^{-1}\Sigma_k'(\hat{\Sigma} + \Sigma_k')^{-1}v\] And recall the Loewner-Heinz inequality [@Loewner34; @Carlen10] tells us that when \(A + B \succ A\), then \((A + B)^{-1} \prec A^{-1}\). Then note that we can show with not too much difficulty (using the fact that since \(A\) and \(B\) are positive definite then \(B^{-1/2}AB^{-1/2} \prec I \iff A \prec B\)) that this implies \[\hat{\Sigma}^{-1}\Sigma_k'(\hat{\Sigma} + \Sigma_k')^{-1} \preceq \hat{\Sigma}^{-1}\Sigma_k'\hat{\Sigma}^{-1}\] But since \(\Sigma_k' \preceq \frac{1}{N}\hat{\Sigma}\) we get that \((i)\) is bounded by \[\frac{1}{N}v^\top \hat{\Sigma}^{-1}v\] now here we again use Hanson Wright to say with probability at least \(1-\epsilon_2\), as long as \(N > N_2 = \frac{1}{\delta_4}\sigma^2\left(D + \sqrt{D \log \frac{1}{\epsilon_2}} + 2\log \frac{1}{\epsilon_2}\right)\) it is bounded. For term \((ii)\), we have rapid convergence of \(\mu_k'\) to \(\hat{\mu}_k\) and note that \(\hat{\Sigma}^{-1}\) is bounded above by \(\frac{1}{N}\|\hat{\Sigma}\|_{op}\). While the difference between the means is bounded by \(\frac{1}{N}\left(2\|\hat{\Sigma}\|_{op}\lambda_{\min}(\hat{\Sigma}_0)\right)\) so if \(N \geq N_3 = \frac{1}{\delta_5}\left(2\|\hat{\Sigma}\|_{op}\lambda_{\min}(\hat{\Sigma}_0)\right)\), then we have convergence. Finally, we can just use Cauchy-Schwarz to argue that if \(N\) is set as above, then the bounds for both of the above hold for the cross term. Combining all terms together, we have that if \(\kappa < \min\{\kappa_1, \kappa_2\}\) and \(N > \max\{N_1, N_2, N_3\}\), then the difference between \(\lambda_k\) and \(\rho_k\) is less than \(\delta = \delta_1 + \delta_2 + \delta_3 + \delta_4 + \delta_5\) with probability \(1 - \epsilon\) as long as \(\epsilon_1 + \epsilon_2 \leq \epsilon\).
Finally, \[\tilde{C}(x) = \log \sum_{k} \exp\bigl(\lambda_k(x)\bigr).\] The map \((\lambda_1,\dots,\lambda_K)\mapsto \log\bigl(\sum_k e^{\lambda_k}\bigr)\) is 1-Lipschitz continuous in the \(\ell_\infty\)-norm in the sense that if \(|\lambda_k - \lambda_k'|\le \delta\) for all \(k\), then \(\bigl|\mathrm{LogSumExp}\{\lambda_k\} \;-\; \mathrm{LogSumExp}\{\lambda_k'\}\bigr|\le \,\epsilon\). Hence if \(\lambda_k(x)\) is within \(\delta\) of \(\tfrac{1}{2}\,[\mathrm{RMD}_k(x)-d]\) uniformly in \(k\), then \(\tilde{C}(x)\) is within \(\delta\) of \(\mathrm{LogSumExp}_k \{\tfrac12\,\mathrm{RMD}_k(x)-d\}\). Noting that \[\begin{align} |\mathrm{LogSumExp}\{\lambda_k\} - \max_k \{C_k(x)\}| &\leq |\mathrm{LogSumExp}\{\lambda_k\} - \mathrm{LogSumExp}\{C_k(x)\}| \\ &\quad\quad+ |\mathrm{LogSumExp}\{C_k(x)\} - \max_kC_k(x)| \\ &\leq \delta + \log K \end{align}\] Thus we conclude \[\bigl|\, \tilde{C}(x) \;-\; \bigl[\tfrac12\,C(x) - d\bigr] \bigr| \;\le\; \, \delta + \log K,\] with probability at least \(1-\epsilon\). This completes the proof.

13 EM Algorithm for the Full Covariance Model↩︎

Here we describe an expectation-maximization (EM) algorithm for estimating the hyperparameters of the hierarchical covariance model. Recall that under this model, \[\begin{align} \Sigma_k &\sim \mathrm{IW}(\nu_0, (\nu_0 - D - 1) \Sigma_0) \\ \mu_k &\sim \mathrm{N}(\mu_0, \kappa_0^{-1} \Sigma_k) \end{align}\] so that the prior hyperparameter \(\Sigma_0\) specifies the mean of the per-class covariances, \(\mathbb{E}[\Sigma_k] = \Sigma_0\). This hyperparameter should not be confused with \(\hat{\Sigma}_0\) defined in the main text, which denoted the empirical estimate of the marginal covariance. Also note that this prior formulation requires \(\nu_0 > D+1\).

First, we set the hyperparameters \(\mu_0\) and \(\Sigma_0\) using empirical Bayes estimates, \[\begin{align} \mu_0 &= \hat{\mu}_0 = \frac{1}{N} \sum_{n=1}^N x_n \\ \Sigma_0 &= \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x_n - \hat{\mu}_{y_n})(x_n - \hat{\mu}_{y_n})^\top \end{align}\] where \(\hat{\mu}_{y_n} = \frac{1}{N_k} \sum_{n: y_n=k} x_n\).

To set the remaining hyperparameters, \(\kappa_0\) and \(\nu_0\), we use EM. The expected log likelihood is separable over these two hyperparameters, \[\begin{align} \mathcal{L}(\nu_0, \kappa_0) &= \mathcal{L}(\nu_0) + \mathcal{L}(\kappa_0) \\ \mathcal{L}(\nu_0) &= \mathbb{E}\left[\sum_{k=1}^K \log \mathrm{IW}(\Sigma_k ~\vert~\nu_0, (\nu_0 - D - 1)\Sigma_0) \right] \\ \mathcal{L}(\kappa_0) &= \mathbb{E}\left[ \sum_{k=1}^K \log \mathrm{N}(\mu_k ~\vert~\mu_0, \kappa_0^{-1} \Sigma_k) \right] \end{align}\] where the expectations are taken with respect to the posterior distribution over \(\{\mu_k, \Sigma_k\}_{k=1}^K\).

13.1 M-step for \(\nu_0\)↩︎

Expanding the first objective yields, \[\begin{align} \mathcal{L}(\nu_0) &= \sum_{k=1}^K \mathbb{E}\left[ \log \left\{ \frac{|(\nu_0 - D - 1)\Sigma_0|^{\frac{\nu_0}{2}}}{2^{\frac{\nu_0 D}{2}} \Gamma_D(\frac{\nu_0}{2}) } |\Sigma_k|^{-\frac{\nu_0 + D + 1}{2}} e^{-\mathrm{Tr}(\tfrac{\nu_0 - D - 1}{2} \Sigma_0 \Sigma_k^{-1})} \right\} \right] \\ &= \sum_{k=1}^K \tfrac{\nu_0 D}{2} \log \left(\tfrac{\nu_0 - D - 1}{2}\right) + \tfrac{\nu_0}{2} \log |\Sigma_0| - \log \Gamma_D(\tfrac{\nu_0}{2}) -\tfrac{\nu_0 + D + 1}{2} \mathbb{E}[\log |\Sigma_k|] - \tfrac{\nu_0 - D - 1}{2} \mathrm{Tr}(\Sigma_0 \mathbb{E}[\Sigma_k^{-1}]) \\ &= \sum_{k=1}^K \tfrac{\nu_0 D}{2} \log \left(\tfrac{\nu_0 - D - 1}{2}\right) + \tfrac{\nu_0}{2} \left[\log |\Sigma_0| - \mathbb{E}[\log |\Sigma_k|] - \mathrm{Tr}(\Sigma_0 \mathbb{E}[\Sigma_k^{-1}]) \right] - \log \Gamma_D(\tfrac{\nu_0}{2}). \end{align}\] We can maximize this objective using a generalized Newton’s method [@minka2000beyond]. We need the first and second derivatives of the objective, \[\begin{align} \mathcal{L}'(\nu_0) &= \sum_{k=1}^K \tfrac{D}{2} \left[ \log \left(\tfrac{\nu_0 - D - 1}{2}\right) + \tfrac{\nu_0}{\nu_0 - D - 1} \right] + \tfrac{1}{2} \left[\log |\Sigma_0| - \mathbb{E}[\log |\Sigma_k|] - \mathrm{Tr}(\Sigma_0 \mathbb{E}[\Sigma_k^{-1}]) \right] - \tfrac{1}{2} \psi_D(\tfrac{\nu_0}{2}) \\ \mathcal{L}''(\nu_0) &= \sum_{k=1}^K \tfrac{D}{2} \left[\tfrac{1}{\nu_0 - D - 1} - \tfrac{D+1}{(\nu_0 - D - 1)^2}\right] - \tfrac{1}{4} \psi_D^{(2)}(\tfrac{\nu_0}{2}). \end{align}\] The idea is to lower bound the objective with a concave function of the form, \[\begin{align} g(\nu_0) = k + a \log \nu_0 + b \nu_0 \end{align}\] which has derivatives \(g'(\nu_0) = \tfrac{a}{\nu_0} + b\) and \(g''(\nu_0) = -\tfrac{a}{\nu_0^2}\). Matching derivatives implies, \[\begin{align} a &= -\nu_0^2 \mathcal{L}''(\nu_0) \\ b &= \mathcal{L}'(\nu_0) - \tfrac{a}{\nu_0} \\ k &= \mathcal{L}(\nu_0) - a\log \nu_0 -b \nu_0. \end{align}\] For \(a > 0\) and \(b<0\), the maximizer of the lower bound is obtained at \[\begin{align} \nu_0^\star &= -\frac{a}{b} = \frac{\nu_0^2 \mathcal{L}''(\nu_0)}{\mathcal{L}'(\nu_0) + \nu_0 \mathcal{L}''(\nu_0)} \end{align}\]

13.2 M-step for \(\kappa_0\)↩︎

Expanding the second objective, \[\begin{align} \mathcal{L}(\kappa_0) &= \sum_{k=1}^K \tfrac{D}{2} \log \kappa_0 - \tfrac{\kappa_0}{2} \mathbb{E}\left[(\mu_k - \mu_0)^\top \Sigma_k^{-1} (\mu_k - \mu_0) \right] + c. \end{align}\] The maximum is obtained at, \[\begin{align} \kappa_0^\star &= \left(\frac{1}{KD} \sum_{k=1}^K \mathbb{E}\left[(\mu_k - \mu_0)^\top \Sigma_k^{-1} (\mu_k - \mu_0) \right]\right)^{-1}. \end{align}\]

13.3 Computing the posterior expectations↩︎

Under the conjugate prior, those posteriors are normal inverse Wishart distributions, \[\begin{align} \mu_k, \Sigma_k ~\vert~\{x_n: y_n = k\} &\sim \mathrm{NIW}(\nu_k', \Sigma_k', \kappa_k', \mu_k') \\ \nu_k' &= \nu_0 + N_k \\ \kappa_k' &= \kappa_0 + N_k \\ \mu_k' &= \frac{1}{\kappa_k'} \left(\kappa_0 \mu_0 + \sum_{n: y_n=k} x_n \right) \\ \Sigma_k' &= (\nu_0 - D - 1) \Sigma_0 + \kappa_0 \mu_0 \mu_0^\top + \sum_{n:y_n=k} x_n x_n^\top - \kappa_k' \mu_k' \mu_k'^\top \end{align}\] To evaluate the objectives above, we need the following expected sufficient statistics of the normal inverse Wishart distribution, \[\begin{align} \mathbb{E}[\Sigma_k^{-1}] &= \nu_k' \Sigma_k'^{-1} \\ \mathbb{E}[\log |\Sigma_k|] &= \log |\Sigma_k'| - \psi_D(\tfrac{\nu_k'}{2}) - D \log 2 \\ \mathbb{E}[(\mu_k - \mu_0)^\top \Sigma_k^{-1} (\mu_k - \mu_0)] &= \frac{1}{\kappa_k'} + (\mu_k' - \mu_0)^\top \mathbb{E}[\Sigma_k^{-1}] (\mu_k' - \mu_0) \end{align}\]

Since the tied covariance approach in RMDS already works quite well, we recommend initializing the EM iterations by setting \(\nu_0 \approx \bar{N}_k\) and \(\kappa_0 \approx 0\). That way, the covariances are strongly coupled across clusters and the means have an uninformative prior.

13.4 Marginal Likelihood↩︎

This EM algorithm maximizes the marginal likelihood, \[\begin{align} \log p(\{x_n, y_n\}_{n=1}^N) &= \sum_{k=1}^K \log p(\{x_n: y_n =k\}) \\ &= \sum_{k=1}^K \log \int p(\{x_n: y_n=k\} ~\vert~\mu_k, \Sigma_k) \, p(\mu_k, \Sigma_k) \mathop{}\!\mathrm{d}\mu_k \mathop{}\!\mathrm{d}\Sigma_k \\ &= \sum_{k=1}^K \log \int \left[\prod_{n:y_n=k} \mathrm{N}(x_n ~\vert~\mu_k, \Sigma_k) \right] \mathrm{NIW}(\mu_k, \Sigma_k ~\vert~\nu_0, \Sigma_0, \kappa_0, \mu_0) \mathop{}\!\mathrm{d}\mu_k \mathop{}\!\mathrm{d}\Sigma_k \\ &= \sum_{k=1}^K \log Z(\nu_k', \Sigma_k', \kappa_k', \mu_k') - \log Z(\nu_0, (\nu_0 - D - 1) \Sigma_0, \kappa_0, \mu_0) + c \end{align}\] where \[\begin{align} \log Z(\nu, \Sigma, \kappa, \mu) = -\tfrac{D}{2} \log \kappa + \log \Gamma_D(\tfrac{\nu}{2}) + \tfrac{\nu D}{2}\log 2 - \tfrac{\nu}{2} \log |\Sigma| \end{align}\] is the log normalizer of the normal inverse Wishart distribution, and \(c\) is constant with respect to the hyperparameters being optimized (but it is data dependent).

14 EM Algorithm for the Diagonal Covariance Model↩︎

Here we describe an expectation-maximization (EM) algorithm for estimating the hyperparameters of the hierarchical diagonal covariance model. Recall that under this model, \[\begin{align} x_{n,d} ~\vert~y_n=k &\sim \mathrm{N}(\mu_{k,d}, \sigma_{k,d}^2) \end{align}\] where \[\begin{align} \sigma_{k,d}^2 &\sim \chi^{-2}(\nu_{0,d}, \sigma_{0,d}^2) \\ \mu_{k,d} &\sim \mathrm{N}(\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2), \end{align}\] for each dimension \(d=1,\ldots,D\) independently. Under the prior \(\mathbb{E}[\sigma_{k,d}^2] = \frac{\nu_{0,d}}{\nu_{0,d} - 2} \sigma_{0,d}^2\), which is approximately \(\sigma_{0,d}^2\) for large degrees of freedom \(\nu_{0,d}\).

First, we set the hyperparameters \(\mu_{0,d}\) and \(\sigma_{0,d}^2\) using empirical Bayes estimates, \[\begin{align} \mu_{0,d} &= \hat{\mu}_{0,d} = \frac{1}{N} \sum_{n=1}^N x_{n,d} \\ \sigma_{0,d}^2 &= \hat{\sigma}_{d}^2 = \frac{1}{N} \sum_{n=1}^N (x_{n,d} - \hat{\mu}_{y_n,d})(x_{n,d} - \hat{\mu}_{y_n,d})^\top \end{align}\] where \(\hat{\mu}_{y_n,d} = \frac{1}{N_k} \sum_{n: y_n=k} x_{n,d}\).

To set the remaining hyperparameters, \(\kappa_{0,d}\) and \(\nu_{0,d}\), we use EM. The expected log likelihood is separable over these two hyperparameters, \[\begin{align} \mathcal{L}(\nu_{0,d}, \kappa_{0,d}) &= \mathcal{L}(\nu_{0,d}) + \mathcal{L}(\kappa_{0,d}) \\ \mathcal{L}(\nu_{0,d}) &= \mathbb{E}\left[\sum_{k=1}^K \log \chi^{-2}(\sigma_{k,d}^2 ~\vert~\nu_{0,d}, \sigma_{0,d}^2) \right] \\ \mathcal{L}(\kappa_{0,d}) &= \mathbb{E}\left[ \sum_{k=1}^K \log \mathrm{N}(\mu_{k,d} ~\vert~\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2) \right] \end{align}\] where the expectations are taken with respect to the posterior distribution over \(\{\mu_{k,d}, \sigma_{k,d}\}_{k=1}^K\).

14.1 M-step for \(\nu_{0,d}\)↩︎

Expanding the first objective yields, \[\begin{align} \mathcal{L}(\nu_{0,d}) &= \sum_{k=1}^K \mathbb{E}\left[ \log \left\{ \frac{(\tfrac{\nu_{0,d}}{2} \sigma_{0,d}^2)^{\tfrac{\nu_{0,d}}{2}}}{\Gamma(\frac{\nu_{0,d}}{2}) } (\sigma_{k,d}^2)^{-\frac{\nu_{0,d} + 2}{2}} e^{-\frac{\nu_{0,d} \sigma_{0,d}^2}{2 \sigma_{k,d}^2}} \right\} \right] \\ &= \sum_{k=1}^K \tfrac{\nu_{0,d}}{2} \log \left(\tfrac{\nu_{0,d}}{2}\right) + \tfrac{\nu_{0,d}}{2} \log \sigma_{0,d}^2 - \log \Gamma(\tfrac{\nu_{0,d}}{2}) -\tfrac{\nu_{0,d} + 2}{2} \mathbb{E}[\log \sigma_{k,d}^2] - \tfrac{\nu_{0,d}}{2} \sigma_{0,d}^2 \mathbb{E}[\sigma_{k,d}^{-2}] \\ &= \sum_{k=1}^K \tfrac{\nu_{0,d}}{2} \log \left(\tfrac{\nu_{0,d}}{2}\right) + \tfrac{\nu_{0,d}}{2} \left[\log \sigma_{0,d}^2 - \mathbb{E}[\log \sigma_{k,d}^2] - \sigma_{0,d}^2 \mathbb{E}[\sigma_{k,d}^{-2}]) \right] - \log \Gamma(\tfrac{\nu_{0,d}}{2}) + c. \end{align}\] We can maximize this objective using a generalized Newton’s method [@minka2000beyond]. We need the first and second derivatives of the objective, \[\begin{align} \mathcal{L}'(\nu_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2} \left[ \log \left(\tfrac{\nu_{0,d}}{2}\right) + 1 \right] + \tfrac{1}{2} \left[\log \sigma_{0,d}^2 - \mathbb{E}[\log \sigma_{k,d}^2] - \sigma_{0,d}^2 \mathbb{E}[\sigma_{k,d}^{-2}]) \right] - \tfrac{1}{2} \psi(\tfrac{\nu_{0,d}}{2}) \\ \mathcal{L}''(\nu_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2 \nu_{0,d}} - \tfrac{1}{4} \psi'(\tfrac{\nu_{0,d}}{2}). \end{align}\] The idea is to lower bound the objective with a concave function of the form, \[\begin{align} g(\nu_{0,d}) = k + a \log \nu_{0,d} + b \nu_{0,d} \end{align}\] which has derivatives \(g'(\nu_{0,d}) = \tfrac{a}{\nu_{0,d}} + b\) and \(g''(\nu_{0,d}) = -\tfrac{a}{\nu_{0,d}^2}\). Matching derivatives implies, \[\begin{align} a &= -\nu_{0,d}^2 \mathcal{L}''(\nu_{0,d}) \\ b &= \mathcal{L}'(\nu_{0,d}) - \tfrac{a}{\nu_{0,d}} \\ k &= \mathcal{L}(\nu_{0,d}) - a\log \nu_{0,d} -b \nu_{0,d}. \end{align}\] For \(a > 0\) and \(b<0\), the maximizer of the lower bound is obtained at \[\begin{align} \nu_{0,d}^\star &= -\frac{a}{b} = \frac{\nu_{0,d}^2 \mathcal{L}''(\nu_{0,d})}{\mathcal{L}'(\nu_{0,d}) + \nu_{0,d} \mathcal{L}''(\nu_{0,d})} \end{align}\]

14.2 M-step for \(\kappa_{0,d}\)↩︎

Expanding the second objective, \[\begin{align} \mathcal{L}(\kappa_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2} \log \kappa_{0,d} - \tfrac{\kappa_{0,d}}{2} \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right] + c. \end{align}\] The maximum is obtained at, \[\begin{align} \kappa_{0,d}^\star &= \left(\frac{1}{K} \sum_{k=1}^K \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right]\right)^{-1}. \end{align}\]

14.3 Computing the posterior expectations↩︎

Under the conjugate prior, those posteriors are normal inverse chi-squared distributions, \[\begin{align} \mu_{k,d}, \sigma_{k,d}^2 ~\vert~\{x_n: y_n = k\} &\sim \mathrm{NIX}(\nu_{k,d}', \sigma_{k,d}'^2, \kappa_{k,d}', \mu_{k,d}') \\ \nu_{k,d}' &= \nu_{0,d} + N_{k} \\ \kappa_{k,d}' &= \kappa_{0,d} + N_{k} \\ \mu_{k,d}' &= \frac{1}{\kappa_{k,d}'} \left(\kappa_0 \mu_{0,d} + \sum_{n: y_n=k} x_n \right) \\ \sigma_{k,d}'^2 &= \frac{1}{\nu_{k,d}'} \left[\nu_{0,d} \sigma_{0,d}^2 + \kappa_{0,d} \mu_{0,d}^2 + \sum_{n:y_n=k} x_{n,d}^2 - \kappa_{k,d}' \mu_{k,d}'^2 \right] \end{align}\] To evaluate the objectives above, we need the following expected sufficient statistics of the normal inverse chi-squared distribution, \[\begin{align} \mathbb{E}[\sigma_{k,d}^{-2}] &= \sigma_{k,d}'^{-2} \\ \mathbb{E}[\log \sigma_{k,d}^2] &= \log \tfrac{\nu_{k,d}' \sigma_{k,d}'^2}{2} - \psi(\tfrac{\nu_{k,d}'}{2}) \\ \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right] &= \frac{1}{\kappa_{k,d}'} + \frac{(\mu_{k,d}' - \mu_{0,d})^2}{\sigma_{k,d}'^2} \end{align}\]

Since the tied covariance approach in RMDS already works quite well, we recommend initializing the EM iterations by setting \(\nu_{0,d} \approx \bar{N}_k\) and \(\kappa_{0,d} \approx 0\). That way, the covariances are strongly coupled across clusters and the means have an uninformative prior.

14.4 Marginal Likelihood↩︎

This EM algorithm maximizes the marginal likelihood, \[\begin{align} \log p(\{x_n, y_n\}_{n=1}^N) &= \sum_{k=1}^K \log p(\{x_n: y_n =k\}) \\ &= \sum_{d=1}^D \sum_{k=1}^K \log \int p(\{x_{n,d}: y_n=k\} ~\vert~\mu_{k,d}, \sigma_{k,d}^2) \, p(\mu_{k,d} \sigma_{k,d}^2) \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \\ &= \sum_{d=1}^D \sum_{k=1}^K \log \int \left[\prod_{n:y_n=k} \mathrm{N}(x_{n,d} ~\vert~\mu_{k,d}, \sigma_{k,d}^2) \right] \mathrm{NIX}(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\nu_{0,d}, \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d}) \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \\ &= \sum_{d=1}^D \left[ \sum_{k=1}^K \left(\log Z(\nu_{k,d}', \sigma_{k,d}'^2, \kappa_{k,d}', \mu_{k,d}') - \log Z(\nu_{0,d}, \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d})\right) + c \right] \end{align}\] where \[\begin{align} \log Z(\nu, \sigma^2, \kappa, \mu) = -\tfrac{1}{2} \log \kappa + \log \Gamma(\tfrac{\nu}{2}) - \tfrac{\nu}{2}\log \tfrac{\nu \sigma^2}{2} \end{align}\] is the log normalizer of the normal inverse chi-squared distribution, and \(c=\log (2\pi)^{-N/2}\) is constant with respect to the hyperparameters being optimized.

14.5 Predictive Distributions↩︎

Under this model, the predictive distribution of a data point given its cluster assignment is, \[\begin{align} p(x ~\vert~y=k, X_{\mathsf{tr}}, y_{\mathsf{tr}}) &= \prod_{d=1}^D \int p(x_d ~\vert~\mu_{k,d}, \sigma_{k,d}^2) \, p(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~X_{\mathsf{tr}}, y_{\mathsf{tr}}) \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \\ &= \prod_{d=1}^D \int \mathrm{N}(x_d ~\vert~\mu_{k,d}, \sigma_{k,d}^2) \, \mathrm{NIX}(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\nu_{k,d}', \sigma_{k,d}'^2, \kappa_{k,d}', \mu_{k,d}') \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \\ &= \prod_{d=1}^D \mathrm{St}(x_d ~\vert~\nu_{k,d}', \mu_{k,d}', \tfrac{\kappa_{k,d}' + 1}{\kappa_{k,d}'} \sigma_{k,d}'^2), \end{align}\] where \(\mathrm{St}(x ~\vert~\nu, \mu, \sigma^2)\) denotes the univariate Student’s t distribution with \(\nu\) degrees of freedom, location \(\mu\), and scale \(\sigma\).

Under this model, the prior predictive distribution is, \[\begin{align} p(x ~\vert~y=K+1, X_{\mathsf{tr}}, y_{\mathsf{tr}}) &= \prod_{d=1}^D \mathrm{St}(x_d ~\vert~\nu_{0,d}, \mu_{0,d}, \tfrac{\kappa_{0,d} + 1}{\kappa_{0,d}} \sigma_{0,d}^2), \end{align}\] which is approximately Gaussian, \(\mathrm{N}(x_d ~\vert~\mu_{0,d}, \sigma_{0,d}^2)\) when \(\nu_{0,d}, \kappa_{0,d} \gg 1\).

15 EM Algorithm for the Coupled Diagonal Covariance Model↩︎

This model introduces a scale factor \(\gamma_k \in \mathbb{R}_+\) that is shared by all dimensions. The model is, \[\begin{align} \gamma_k &\sim \chi^2(\alpha_0) \\ \sigma_{k,d} &\sim \chi^{-2}(\nu_{0,d}, \gamma_k \sigma_{0,d}^2) & \text{for } d &=1,\ldots D \\ \mu_{k,d} &\sim \mathrm{N}(\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2) & \text{for } d &=1,\ldots D \\ \end{align}\] Since \(\mathbb{E}[\gamma_k] = 1\), under the prior \(\mathbb{E}[\sigma_{k,d}^2] = \tfrac{\nu_{0,d}}{\nu_{0,d}-2} \sigma_{0,d}^2\), which is approximately \(\sigma_{0,d}^2\) for large \(\nu_{0,d}\).

The hyperparameters of the model are \(\eta = (\alpha_0, \{\nu_{0,d}, \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d}\})\). We set the hyperparameters \(\mu_{0,d}\) and \(\sigma_{0,d}^2\) using empirical Bayes estimates, \[\begin{align} \mu_{0,d} &= \hat{\mu}_{0,d} = \frac{1}{N} \sum_{n=1}^N x_{n,d} \\ \sigma_{0,d}^2 &= \hat{\sigma}_{d}^2 = \frac{1}{N} \sum_{n=1}^N (x_{n,d} - \hat{\mu}_{y_n,d})^2 \end{align}\] where \(\hat{\mu}_{k} = \frac{1}{N_k} \sum_{n: y_n=k} x_n\) and \(N_k = \sum_n \mathbb{I}[y_n = k]\).

To set the remaining hyperparameters, we use EM.

15.1 E-step↩︎

Note that the posterior distribution factors as, \[\begin{align} p(\gamma_k, \{\mu_{k,d}, \sigma_{k,d}^2\}_{d=1}^D ~\vert~X_k) &= p(\gamma_k ~\vert~X_k) \, p(\{\mu_{k,d}, \sigma_{k,d}^2\}_{d=1}^D ~\vert~\gamma_k, X_k) \\ &= p(\gamma_k ~\vert~X_k) \prod_{d=1}^D p(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\gamma_k, X_k). \end{align}\]

The posterior distribution over \(\gamma_k\) doesn’t have a simple closed form, but since it’s only one-dimensional, we can approximate it on a dense grid of points, \(\{\gamma^{(p)}\}_{p=1}^P\). Conditioned on \(\gamma_k = \gamma^{(p)}\), the distribution of \(\mu_{k,d}\) and \(\sigma_{k,d}^2\) is a normal inverse chi-squared. For each point, \[\begin{align} p(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\gamma_k = \gamma^{(p)}, X_k) &=\mathrm{NIX}(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\nu_{k,d}', \sigma_{k,p,d}'^2, \kappa_{k,d}', \mu_{k,d}') \\ \nu_{k,d}' &= \nu_{0,d} + N_{k} \\ \kappa_{k,d}' &= \kappa_{0,d} + N_{k} \\ \mu_{k,d}' &= \frac{1}{\kappa_{k,d}'} \left(\kappa_0 \mu_{0,d} + \sum_{n: y_n=k} x_n \right) \\ \sigma_{k,p,d}'^2 &= \frac{1}{\nu_{k,d}'} \left[\nu_{0,d} \gamma^{(p)} \sigma_{0,d}^2 + \kappa_{0,d} \mu_{0,d}^2 + \sum_{n:y_n=k} x_{n,d}^2 - \kappa_{k,d}' \mu_{k,d}'^2 \right] \end{align}\] Note that this is practically the same as above, but with \(\gamma_k\) scaling the prior for \(\sigma_{k,d}^2\). For any value of \(\gamma_k\), the posterior probability is, \[\begin{align} p(\gamma_k = \gamma^{(p)} ~\vert~X_k) &\propto p(\gamma_k = \gamma^{(p)}) \, p(X_k ~\vert~\gamma_k = \gamma^{(p)}) \\ &= p(\gamma_k = \gamma^{(p)}) \prod_{d=1}^D p(X_{k,d} ~\vert~\gamma_k = \gamma^{(p)}) \\ &= p(\gamma_k = \gamma^{(p)}) \prod_{d=1}^D \frac{Z(\nu_{k,d}', \sigma_{k,p,d}'^2, \kappa_{k,d}', \mu_{k,d}')}{Z(\nu_{0,d}, \gamma^{(p)} \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d})} \\ &\triangleq \tilde{w}_{k,p}'. \end{align}\] where we reused the marginal likelihood calculation from the hierarchical diagaonal DPMM above. Finally, denote the normalized posterior probabilities as, \[\begin{align} w_{k,p}' &= \frac{\tilde{w}_{k,p}'}{\sum_r \tilde{w}_{k,r}'}. \end{align}\]

15.2 M-Step↩︎

To set the hyperparameters, \(\kappa_{0,d}\), \(\nu_{0,d}\), and \(\alpha_0\), we use EM. The expected log likelihood is separable over these two hyperparameters, \[\begin{align} \mathcal{L}(\nu_{0,d}, \kappa_{0,d}, \alpha_0) &= \mathcal{L}(\nu_{0,d}) + \mathcal{L}(\kappa_{0,d}) + \mathcal{L}(\alpha_0) \\ \mathcal{L}(\nu_{0,d}) &= \mathbb{E}\left[\sum_{k=1}^K \log \chi^{-2}(\sigma_{k,d}^2 ~\vert~\nu_{0,d}, \gamma_k \sigma_{0,d}^2) \right] \\ \mathcal{L}(\kappa_{0,d}) &= \mathbb{E}\left[ \sum_{k=1}^K \log \mathrm{N}(\mu_{k,d} ~\vert~\mu_{0,d}, \kappa_{0,d}^{-1} \sigma_{k,d}^2) \right] \\ \mathcal{L}(\alpha_0) &= \mathbb{E}\left[ \log \mathrm{Ga}(\gamma_k ~\vert~\alpha_0, \alpha_0) \right] \end{align}\] where the expectations are taken with respect to the posterior distribution over \(\{\gamma_k, \{\mu_{k,d}, \sigma_{k,d}\}_{d=1}^D \}_{k=1}^K\).

15.3 M-step for \(\nu_{0,d}\)↩︎

Expanding the first objective yields, \[\begin{align} \mathcal{L}(\nu_{0,d}) &= \sum_{k=1}^K \mathbb{E}\left[ \log \left\{ \frac{(\tfrac{\nu_{0,d}}{2} \gamma_k \sigma_{0,d}^2)^{\tfrac{\nu_{0,d}}{2}}}{\Gamma(\frac{\nu_{0,d}}{2}) } (\sigma_{k,d}^2)^{-\frac{\nu_{0,d} + 2}{2}} e^{-\frac{\nu_{0,d} \gamma_k \sigma_{0,d}^2}{2 \sigma_{k,d}^2}} \right\} \right] \\ &= \sum_{k=1}^K \tfrac{\nu_{0,d}}{2} \log \left(\tfrac{\nu_{0,d}}{2}\right) + \tfrac{\nu_{0,d}}{2} \mathbb{E}[\log \gamma_k] + \tfrac{\nu_{0,d}}{2} \log \sigma_{0,d}^2 - \log \Gamma(\tfrac{\nu_{0,d}}{2}) -\tfrac{\nu_{0,d} + 2}{2} \mathbb{E}[\log \sigma_{k,d}^2] - \tfrac{\nu_{0,d}}{2} \sigma_{0,d}^2 \mathbb{E}[\gamma_k \sigma_{k,d}^{-2}] \\ &= \sum_{k=1}^K \tfrac{\nu_{0,d}}{2} \log \left(\tfrac{\nu_{0,d}}{2}\right) + \tfrac{\nu_{0,d}}{2} \left[\mathbb{E}[\log \gamma_k] + \log \sigma_{0,d}^2 - \mathbb{E}[\log \sigma_{k,d}^2] - \sigma_{0,d}^2 \mathbb{E}[\gamma_k \sigma_{k,d}^{-2}]) \right] - \log \Gamma(\tfrac{\nu_{0,d}}{2}) + c. \end{align}\] We can maximize this objective using a generalized Newton’s method [@minka2000beyond]. We need the first and second derivatives of the objective, \[\begin{align} \mathcal{L}'(\nu_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2} \left[\log \left(\tfrac{\nu_{0,d}}{2}\right) + 1 \right] + \tfrac{1}{2} \left[\mathbb{E}[\log \gamma_k] + \log \sigma_{0,d}^2 - \mathbb{E}[\log \sigma_{k,d}^2] - \sigma_{0,d}^2 \mathbb{E}[\gamma_k \sigma_{k,d}^{-2}]) \right] - \tfrac{1}{2} \psi(\tfrac{\nu_{0,d}}{2}) \\ \mathcal{L}''(\nu_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2 \nu_{0,d}} - \tfrac{1}{4} \psi'(\tfrac{\nu_{0,d}}{2}). \end{align}\] The idea is to lower bound the objective with a concave function of the form, \[\begin{align} g(\nu_{0,d}) = k + a \log \nu_{0,d} + b \nu_{0,d} \end{align}\] which has derivatives \(g'(\nu_{0,d}) = \tfrac{a}{\nu_{0,d}} + b\) and \(g''(\nu_{0,d}) = -\tfrac{a}{\nu_{0,d}^2}\). Matching derivatives implies, \[\begin{align} a &= -\nu_{0,d}^2 \mathcal{L}''(\nu_{0,d}) \\ b &= \mathcal{L}'(\nu_{0,d}) - \tfrac{a}{\nu_{0,d}} \\ k &= \mathcal{L}(\nu_{0,d}) - a\log \nu_{0,d} -b \nu_{0,d}. \end{align}\] For \(a > 0\) and \(b<0\), the maximizer of the lower bound is obtained at \[\begin{align} \nu_{0,d}^\star &= -\frac{a}{b} = \frac{\nu_{0,d}^2 \mathcal{L}''(\nu_{0,d})}{\mathcal{L}'(\nu_{0,d}) + \nu_{0,d} \mathcal{L}''(\nu_{0,d})} \end{align}\]

15.4 M-step for \(\kappa_{0,d}\)↩︎

Expanding the second objective, \[\begin{align} \mathcal{L}(\kappa_{0,d}) &= \sum_{k=1}^K \tfrac{1}{2} \log \kappa_{0,d} - \tfrac{\kappa_{0,d}}{2} \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right] + c. \end{align}\] The maximum is obtained at, \[\begin{align} \kappa_{0,d}^\star &= \left(\frac{1}{K} \sum_{k=1}^K \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right]\right)^{-1}. \end{align}\]

15.5 M-step for \(\alpha_0\)↩︎

Expanding the final objective, \[\begin{align} \mathcal{L}(\alpha_0) &= \sum_{k=1}^K \alpha_0 \log \alpha_0 - \log \Gamma(\alpha_0) + \alpha_0 \mathbb{E}\left[ \log \gamma_k \right]- \alpha_0 \mathbb{E}\left[\gamma_k\right] + c. \end{align}\] Its derivatives are, \[\begin{align} \mathcal{L}'(\alpha_0) &= K\log \alpha_0 + K - K\psi(\alpha_0) + \sum_{k=1}^K \mathbb{E}\left[\log \gamma_k\right] - \mathbb{E}[\gamma_k] \\ \mathcal{L}''(\alpha_0) &= \frac{K}{\alpha_0} - K\psi'(\alpha_0) \end{align}\] We can optimize \(\alpha_0\) using the generalized Newton’s method described in the M-step for \(\nu_0\).

15.6 Computing the posterior expectations↩︎

To evaluate the objectives above, we need the following expected sufficient statistics of the normal inverse chi-squared distribution, \[\begin{align} \mathbb{E}[\gamma_k] &= \sum_p w_{k,p}' \gamma^{(p)} \\ \mathbb{E}[\log \gamma_k] &= \sum_p w_{k,p}' \log \gamma^{(p)} \\ \mathbb{E}[\gamma_k \sigma_{k,d}^{-2}] &= \mathbb{E}_{\gamma_k} [ \mathbb{E}_{\sigma_{k,d}^2 | \gamma_k} [\gamma_k \sigma_{k,d}^{-2}]] = \sum_p w_{k,p}' \gamma^{(p)} \sigma_{k,p,d}'^{-2} \\ \mathbb{E}[\log \sigma_{k,d}^2] &= \mathbb{E}_{\gamma_k} [\mathbb{E}_{\sigma_{k,d}^2 | \gamma_k} [\log \sigma_{k,d}^2]] = \sum_p w_{k,p}' [\log \tfrac{\nu_{k,d}' \sigma_{k,p,d}'^2}{2} - \psi(\tfrac{\nu_{k,d}'}{2})] \\ \mathbb{E}\left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right] &= \mathbb{E}_{\gamma_k} \left[\mathbb{E}_{\mu_{k,d}, \sigma_{k,d}^2 | \gamma_k} \left[\frac{(\mu_{k,d} - \mu_{0,d})^2}{\sigma_{k,d}^2}\right] \right] = \sum_p w_{k,p}' \left[\frac{1}{\kappa_{k,d}'} + \frac{(\mu_{k,d}' - \mu_{0,d})^2}{\sigma_{k,p,d}'^2} \right] \end{align}\] where \(\gamma^{(p)}\) are the centers of discretized posterior on \(\gamma_k\) and \(w_{k,p}\) are the corresponding weights.

Since the tied covariance approach in RMDS already works quite well, we recommend initializing the EM iterations by setting \(\nu_{0,d} \approx \bar{N}_k\) and \(\kappa_{0,d} \approx 0\). That way, the covariances are strongly coupled across clusters and the means have an uninformative prior.

15.7 Computing the predictive distributions↩︎

The prior predictive is, \[\begin{gather} p(x^\star; \eta_0) = \\ \int \left[\prod_{d=1}^D \int \mathrm{N}(x^\star_d ~\vert~\mu_d, \sigma^2_d) \mathrm{NIX}(\mu_d, \sigma^2_d ~\vert~\nu_{0,d}, \kappa_{0,d}, \mu_{0,d}, \gamma \sigma^2_{0,d}, \{x_n, y_n\}_{n=1}^N) \mathop{}\!\mathrm{d}\mu_d \mathop{}\!\mathrm{d}\sigma^2_d \right] \mathrm{Ga}(\gamma ~\vert~\alpha_0, \alpha_0) \mathop{}\!\mathrm{d}\gamma. \end{gather}\] where \(\eta_0 = \alpha_0, \{\mu_{0,d}, \sigma^2_{0,d}, \kappa_{0,d}, \nu_{0,d}\}_{d=1}^D\) are the model hyperparamters.

The \(\gamma\) integral can be estimated by numerical integration over a dense grid of points, \[\begin{align} p(x^\star; \eta_0) &\approx \sum_{p=1}^P w_{0,p} \left[ \prod_{d=1}^D \int \mathrm{N}(x^\star_d ~\vert~\mu_d, \sigma^2_d) \mathrm{NIX}(\mu_d, \sigma^2_d ~\vert~\nu_{0,d}, \kappa_{0,d}, \mu_{0,d}, \gamma^{(p)} \sigma^2_{0,d}, \{x_n, y_n\}_{n=1}^N) \mathop{}\!\mathrm{d}\mu_d \mathop{}\!\mathrm{d}\sigma^2_d \right] \\ &= \sum_{p=1}^P w_{0,p} \, \left[ \prod_{d=1}^D \mathrm{St}(x^\star_d ~\vert~\nu_{0,d}, \mu_{0,d}, \tfrac{\kappa_{0,d}+1}{\kappa_{0,d}} \gamma^{(p)}\sigma^2_{0,d}) \right] \end{align}\] where, \[\begin{align} w_{0, p} = \frac{\mathrm{Ga}(\gamma^{(p)} ~\vert~\alpha_0, \alpha_0) \Delta\gamma^{(p)}}{\sum_r \mathrm{Ga}(\gamma^{(r)} ~\vert~\alpha_0, \alpha_0) \Delta\gamma^{(r)}}. \end{align}\] We renormalize the weights to ensure that the numerical integration satisfies that \(\mathbb{E}_\gamma[1] = 1\).

In practice, we evaluate the prior log predictive probability using a log-sum-exp, \[\begin{align} \log p(x^\star; \eta_0) \approx \mathrm{logsumexp}_p \left[ \log w_{0,p} + \sum_{d=1}^D \log \mathrm{St}(x_d^\star ~\vert~\nu_{0,d}, \mu_{0,d}, \tfrac{\kappa_{0,d} + 1}{\kappa_{0,d}} \gamma^{(p)} \sigma_{0,d}^2) \right] \end{align}\]

By the same logic, the posterior log predictive is, \[\begin{align} \log p(x^\star ~\vert~y^\star=k, \{x_n, y_n\}_{n=1}^N; \eta_0) \approx \mathrm{logsumexp}_p \left[ \log w_{k,p}' + \sum_{d=1}^D \log \mathrm{St}(x_d^\star ~\vert~\nu_{k,d}', \mu_{k,d}', \tfrac{\kappa_{k,d}'+ 1}{\kappa_{k,d}'} \gamma^{(p)} \sigma_{k,p,d}'^2) \right] \end{align}\]

15.8 Marginal Likelihood↩︎

This EM algorithm maximizes the marginal likelihood, \[\begin{align} \log p(\{x_n, y_n\}_{n=1}^N) &= \sum_{k=1}^K \log p(\{x_n: y_n =k\}) \\ &= \sum_{k=1}^K \log \int p(\{x_n: y_n=k\} ~\vert~\{\mu_{k,d}, \sigma_{k,d}^2\}) \;p(\{\mu_{k,d} \sigma_{k,d}^2\} ~\vert~\gamma_k) \, p(\gamma_k) \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \mathop{}\!\mathrm{d}\gamma_k \\ &= \sum_k \log \int \left[ \prod_{d=1}^D \int \left[ \prod_{n:y_n=k} p(x_{n,d} ~\vert~\mu_{k,d}, \sigma_{k,d}^2) p(\mu_{k,d}, \sigma_{k,d}^2 ~\vert~\gamma_k) \mathop{}\!\mathrm{d}\mu_{k,d} \mathop{}\!\mathrm{d}\sigma_{k,d}^2 \right] p(\gamma_k) \mathop{}\!\mathrm{d}\gamma_k \right] \\ &= \sum_k \log \int \prod_{d=1}^D \left[\frac{Z(\nu_{k,d}', \sigma_{k,p,d}'^2, \kappa_{k,d}', \mu_{k,d}')}{Z(\nu_{0,d}, \gamma_k \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d})} \right] p(\gamma_k) \mathop{}\!\mathrm{d}\gamma_k + c \\ &\approx \sum_k \log \sum_p w_{0,p} \prod_{d=1}^D \frac{Z(\nu_{k,d}', \sigma_{k,p,d}'^2, \kappa_{k,d}', \mu_{k,d}')}{Z(\nu_{0,d}, \gamma^{(p)} \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d})} + c \\ &= \sum_k \mathrm{logsumexp}_p \left[ \ell_{k,p} \right] + c \end{align}\] where \[\begin{align} \ell_{k,p} &\triangleq \log w_{0,p} + \sum_{d=1}^D \left(\log Z(\nu_{k,d}', \sigma_{k,p,d}'^2, \kappa_{k,d}', \mu_{k,d}') - \log Z(\nu_{0,d}, \gamma^{(p)} \sigma_{0,d}^2, \kappa_{0,d}, \mu_{0,d}) \right) \end{align}\] and \(c = \tfrac{ND}{2} \log 2\pi\).

16 Score Correlation between RMDS and the Tied Covariance Model↩︎

We observe that the DPMM model with shared covariance and the RMDS [@ren21rmds] are highly correlated, as illustrated in 6. Here, we plot the RMDS vs the Tied Covariance Gaussian DPMM for all the real datasets (differentiated by color) in the Imagenet-1K OpenOOD task and note the tight agreement between the two. This empirical result supports the theoretical relationship derived in 2.

17 Imagenet-1K Experiment↩︎

18 CIFAR-10 Experiment↩︎

19 CIFAR-100 Experiment↩︎

20 Feature Dimension Study with Independent RMDS↩︎

---

1. The implementation for all models and experiments is available at https://github.com/rwl93/bnp4ood.↩︎

2. We flip the sign of \(\mathrm{RMD}_k(x)\) compared to [@ren21rmds], but account for it in the definition of \(C(x)\) so that the resulting score is unchanged. Our presentation is in keeping with the definition of the MDS score [@lee18mds].↩︎