Publications · Journal

Bayesian nonparametric mixture of experts for inverse problems

TrungTin Nguyen, Florence Forbes, Julyan Arbel, Hien Duy Nguyen

† Corresponding author.

J. Nonparametric Statistics · Journal Journal of Nonparametric Statistics. Journal article (2024).

Abstract

Large classes of problems can be formulated as inverse problems, where the goal is to find parameter values that best explain some observed measures. The relationship between parameters and observations is typically highly non-linear, with relatively high dimensional observations and correlated multidimensional parameters. To deal with these constraints via inverse regression strategies, we consider the Gaussian Local Linear Mapping (GLLiM) model, a special instance of mixture of expert models. We propose a general scheme to design a Bayesian nonparametric GLLiM model to avoid any commitment to an arbitrary number of experts. A tractable estimation algorithm is designed using variational Bayesian expectation-maximisation. We establish posterior consistency for the number of mixture components after the merge-truncate-merge algorithm post-processing. Illustrations on simulated data show good results in terms of recovering the true number of experts and the regression function.

Keywords: B

ayesian nonparametrics; mixture of experts; inverse problems; Gaussian locally-linear mapping models; linear cluster-weighted models; variational inference; clustering; regression; model selection.

1 Introduction

Inverse problems.

Many problems can be formulated as inverse problems, where the goal is to find parameter values that best explain an observed phenomenon. Typical constraints in practice are that the relationships between parameters and observations are highly non-linear, with relatively high dimensional observations and correlated multivariate parameters. To handle these constraints, we consider probabilistic mixtures of locally linear models, namely the Gaussian locally linear mapping (GLLiM) approach of Deleforge et al. (2015), which includes affine instances of the mixture of experts (MoE) model (Xu et al., 1995) and the classical inverse linear regression model (Hoadley, 1970; Li, 1991), as special cases. We propose to address inverse problems in a Bayesian framework, making use of the availability of simulations from forward models of interest. The GLLiM approach has been used in many applications, e.g., in medical imaging (Boux et al., 2021), planetary science (Kugler et al., 2022; Nguyen et al., 2024a), head-pose estimation problem in computer vision (Lathuilière et al., 2017) and quantitative trait prediction from biological data (Blein-Nicolas et al., 2024).

Mixture of experts models (MoE)

are generalizations of neural network architectures proposed by Jacobs et al. (1991). Further, these flexible models also generalize the classical mixture models (MMs) and mixture of regression models (McLachlan and Peel, 2000). Their flexibility comes from the fact that they allow the mixture weights (or the gating functions) to depend on the explanatory variables, together with the component densities (or the experts). In regression, MoE models with Gaussian experts and softmax or normalized Gaussian gating functions (as in GLLiM) are the most popular choices. These models are powerful tools for modelling complex nonlinear relationships between outputs (responses) and inputs (predictors) that arise from different subpopulations. The popularity of these conditional mixture density models is largely due to their universal approximation properties (Norets, 2010; Nguyen et al., 2016, 2019, 2021a) as well as their good convergence rate, see, e.g., for mixture of regression in Ho et al. (2022) and for MoE in Jiang and Tanner (1999); Nguyen et al. (2024c, 2023a, b). It is worth noting that these results improve the approximation capabilities and convergence rates of unconditional MMs, as discussed in Genovese and Wasserman (2000); Rakhlin et al. (2005); Nguyen (2013); Shen et al. (2013); Ho and Nguyen (2016a, b); Nguyen et al. (2020, 2023b). At a high level, universal approximation theorems state that given a large enough number of components, MMs and MoE models can approximate a large class of unconditional and conditional probability density functions (cPDF), respectively, to any degree of accuracy. See, e.g., Yuksel et al. (2012); Nguyen and Chamroukhi (2018); Do et al. (2023); Nguyen (2021); Chen et al. (2022), for further detailed reviews of practical and theoretical aspects of MoE models in statistics and in diverse domains, e.g.,  natural language processing and computer vision.

Model selection in MoE models.

Although universal approximation allows us to conclude that, given a sufficient number of components, a finite MoE can approximate any other cPDF to an arbitrary degree of accuracy, it is not clear how to choose a large enough number of components for realistic problems. This motivates a careful study of interesting and important model selection problems for MMs and MoE models, which have attracted much attention in statistics and machine learning over the last 50 years, see, e.g., Celeux et al. (2019); Gormley and Frühwirth-Schnatter (2019) for a recent comprehensive review.

When selecting the best data-driven number of components for MoE models, there are several approaches to controlling and accounting for model complexity. Typically, model selection in MoE models is performed using an information criterion, such as the Akaike information criterion (AIC; Akaike, 1974; Frühwirth-Schnatter et al., 2018), the Bayesian information criterion (BIC; Schwarz, 1978; Berrettini et al., 2024, BIC-GLLiM; Forbes et al., 2022a, b), the BIC-like approximation of integrated classification likelihood (ICL-BIC; Biernacki et al., 2000; Frühwirth-Schnatter et al., 2012), extended BIC for mixture Gaussian graphical models (eBIC; Nguyen and Li, 2024; Foygel and Drton, 2010), or Sin and White information criterion (SWIC; Sin and White, 1996; Westerhout et al., 2024) for stationary sequences that are not independently and identically distributed. However, an important limitation of these criteria is that they are only asymptotically valid. This means that there are no finite sample guarantees when using AIC, BIC, ICL-BIC or SWIC to choose between different levels of complexity. Therefore, their use in small sample settings is ad hoc.

To overcome such difficulties, and to partially support the so-called slope heuristic approach (Birgé and Massart, 2007, also see Arlot, 2019 for a recent review), Nguyen et al. (2021b, 2022c, 2022b, 2023d, 2023c) recently established non-asymptotic risk bounds in the form of weak oracle inequalities, provided that lower bounds on the penalties hold, in high-dimensional regression scenarios for a variety of MoE models, including GLLiM. Another approach, as proposed by Nguyen et al. (2022a), is based on the closed testing principle, which leads to a sequential testing procedure that enables confidence statements regarding the order of a finite mixture model. These methods facilitate an optimal data-driven selection of the number of components in finite-sample settings. Additionally, Bayesian model selection, which employs marginal likelihoods as a criterion, provides another strategy for determining the number of components in MoE models (Frühwirth-Schnatter, 2019; Zens, 2019). However, all previous approaches require that a range of models with different values be trained and compared, which can be a computational bottleneck in a high-dimensional framework.

Recently, Kock et al. (2022) proposed computationally efficient variational inference approaches for architecture selection in high-dimensional deep Gaussian mixture models using overfitted mixtures (see, e.g., Rousseau and Mengersen, 2011; Forbes et al., 2019), where unnecessary components are dropped in the estimation. The same idea is applied in the context of overfitting sparse Bayesian factor analysis in Frühwirth-Schnatter et al. (2024), where the authors employed the Markov Chain Monte Carlo (MCMC) procedure, followed by post-processing of posterior draws, to facilitate practical Bayesian inference. However, we are interested here in the more general context of the Bayesian nonparametric (BNP) approach, (see, e.g., Hjort et al., 2010; Ghosal and Van der Vaart, 2017), where it is not necessary to know an upper bound of the true number of components as in Bayesian overfitted MM. This is one motivation for the BNP priors that are considered here for GLLiM.

Dirichlet process mixture models (DP-MMs) and Pitman-Yor process mixture models (PYP-MMs) are among the most popular BNP models, particularly suitable for density estimation and probabilistic clustering. However, the posterior of the DP-MMs or PYP-MMs are inconsistent for the number of components if the true number of components is finite and the concentration parameter is fixed (see, e.g., Miller and Harrison, 2014, and Alamichel et al., 2024 for a review). This is because a BNP prior such as DP or PYP places zero probability on mixing measures with a finite number of supporting atoms. An interesting recent result in Ascolani et al. (2022) is that consistency for the number of components can be achieved if a prior is placed on the concentration parameter of the DP-MM, under some assumptions on this prior.

It appears that BNP-GLLiM tends to produce many small extraneous components around the true clusters in our numerical experiment Section 6. This makes it difficult to use them to infer the true number of components when this becomes a quantity of interest (Maceachern and Müller, 1998; Green and Richardson, 2001). This encourages the use of a novel, simple post-processing algorithm in the spirit of the Merge-Truncate-Merge (MTM) introduced by Guha et al. (2021). This post-processing procedure consistently estimates the number of components for any general Bayesian prior, even without knowing its exact structure, as long as the posterior for that prior contracts to the true mixing distribution at a known rate.

Contributions.

To address the challenges of highly nonlinear inverse problems with relatively high dimensional observations and correlated parameters, we propose a novel BNP-GLLiM model and inference procedure, which is computationally efficient and avoids any commitment to an arbitrary number of components. Although BNP mixture models (BNP-MM), which are special cases of the BNP-GLLiM model, have been extensively studied in the literature (Escobar and West, 1995; Maceachern and Müller, 1998; Neal, 2000; Arbel et al., 2021; Li et al., 2022; Durand et al., 2022), the extension to MoE models in an inverse regression framework has not been covered. In addition, we establish theoretical properties such as posterior consistency for recovering the true number of components in BNP-GLLiM using the post-processing merge-truncate-merge (MTM) algorithm. Finally, our illustrations on simulated data show good results in terms of recovering the true number of components and mean regression functions. It is worth emphasising that, for the first time, we provide evidence that MTM consistency holds not only for the MMs results of Guha et al. (2021), but also in the more general context of MoE models for cPDFs.

Notations.

Throughout this paper, {1,,D} is abbreviated as [D] for D, where denotes the positive natural numbers. The notation refers to a definition. It is used to simplify the notation or expression. For a parametric model S, dim(S) refers to its dimension, i.e.,  the total number of parameters to be estimated. Furthermore, [;] denotes vertical vector concatenation. Throughout, we use the following colour rule for observations and parameters: observations are represented in green; latent, random or unknown parameters in red; and (fixed) hyperparameters in blue.

Outline.

The paper is organized as follows. In Section 2, we first discuss how to construct the BNP-GLLiM model. A VBEM algorithm and the corresponding ELBO are described in Section 3, and predictive cPDFs in Section 4. Next, Section 5 shows how we can integrate the Merge-Truncate-Merge (MTM) post-processing procedure and prove consistency for the MTM output. This is useful to perform regression, clustering and model selection, simultaneously. We experimentally evaluate our new results on simulated datasets in Section 6. Some perspectives are provided in Section 7. We recall the standard Bayesian nonparametric priors and variational Bayesian expectation-maximisation principle in Appendices A and B, respectively. All details of VBEM for the BNP-GLLiM model, evidence lower-bound, and technical proofs not included in the main paper are relegated to Appendices C, D and E, respectively. Appendix F presents a more general model with a hyperprior on the gating parameters, called BNP-GLLiM2.

2 BNP-GLLiM model

In Sections 2.1 and 2.2, we present the advantage of adopting a GLLiM model, which uses an inverse regression approach to estimate nonlinear high-to-low dimensional mappings. Such a strategy allows to greatly reduce the number of required parameters. To avoid any commitment to an arbitrary number of components, we then construct the BNP-GLLiM model in Section 2.3.

2.1 Inverse regression framework

We are interested in estimating the law of a low-dimensional random variable 𝐗=(𝐗l)l[L] conditionally on a high-dimensional 𝐘=(𝐘d)d[D], where typically DL. We follow an inverse regression framework as in e.g.,  Li (1991); Deleforge et al. (2015). Therefore, in training, the low-dimensional variable 𝐗 plays the role of the regressor, while the response 𝐘 is a function of 𝐗, possibly corrupted by noise through inverse cPDF p(𝐘𝐗;𝝍), where 𝝍 is an inverse parameter. The low dimension of the regressor 𝐗 allows to drastically reduce the number of parameters to be estimated. In addition, the forward parameter 𝝍 and cPDF p(𝐗𝐘;𝝍) are tractable after estimating the inverse parameter 𝝍. Therefore, this density can be used to predict the low-dimensional response 𝐱 of a high-dimensional test point 𝐲. This inverse-then-forward regression strategy justifies the unconventional notation: 𝐘 for the high-dimensional input and 𝐗 for the low-dimensional response. Here and subsequently, we refer to the low-dimensional data sample as 𝒳{𝐱n}n[N](L)N, the high-dimensional data sample as 𝒴{𝐲n}n[N](D)N. We denote the observed values as (𝐱,𝐲), which are independently and identically distributed (i.i.d.) samples from the random variables (𝐗,𝐘).

2.2 Nonlinear high-to-low dimensional mapping via GLLiM model

The GLLiM models, as originally introduced in Deleforge et al. (2015), are used to capture the non-linear relationship between the response and the set of covariates from a high-to-low heterogeneous data. Specifically, Deleforge et al. (2015) overcame the difficulty of high-to-low regression by tackling the problem the other way round, i.e.,  low-to-high. This means that the roles of input and response variables are swapped so that the low-dimensional variable 𝐗 becomes the regressor as in Section 2.1. GLLiM then relies on a piecewise linear model in the following way. The high-dimensional response 𝐘 is approximated by the local affine mappings K:

𝐘=k=1K𝕀(Z=k)(𝐀k𝐗+𝐛k+𝐄k).

Here, 𝕀 is an indicator function and Z is a latent variable that captures a cluster relationship, such that Z=k if 𝐘 comes from cluster k[K]. Matrices 𝐀kD×L and vectors 𝐛kD define cluster-specific affine transformations. In addition, 𝐄k are error terms that capture both the reconstruction error due to the local affine approximations as well as the observation noise in D.

Following the usual assumption that 𝐄k is a zero-mean Gaussian variable with a covariance matrix 𝚺kD×D, it follows that

p(𝐲𝐱,Z=k;𝝍)=𝒩D(𝐲𝐀k𝐱+𝐛k,𝚺k), (1)

where 𝝍 is the vector of model parameters, 𝒩D(𝐲;𝐀k𝐱+𝐛k,𝚺k) is the Gaussian cPDF of dimension D. In order to enforce the affine transformations to be local, 𝐗 is defined as a mixture of K Gaussian components as follows:

p(𝐱Z=k;𝝍)=𝒩L(𝐱𝐜k,𝚪k),p(Z=k;𝝍)=πk,

where 𝐜kL, 𝚪kL×L, and 𝝅=(πk)k[K] belongs to a probability simplex, defined as {(πk)k[K](+)K,k=1Kπk=1}. Then, via the conditional property of Gaussian variables and hierarchical decomposition, given any 𝝍=(πk,𝐜k,𝚪k,𝐀k,𝐛k,𝚺k)k[K]𝚿, we obtain the following inverse conditional density:

p(𝐲𝐱;𝝍)=k=1Kπk𝒩L(𝐱𝐜k,𝚪k)l=1Kπl𝒩L(𝐱𝐜l,𝚪l)𝒩D(𝐲𝐀k𝐱+𝐛k,𝚺k). (2)

Using the inverse regression framework, in (1), the roles of input and response variables should be reversed so that 𝐘 becomes the covariate and 𝐗 plays the role of the multivariate response. Therefore, based on a similar previous hierarchical one in (2), its corresponding forward conditional density from D to L is defined by

p(𝐱𝐲;𝝍) =k=1Kπk𝒩D(𝐲𝐜k,𝚪k)l=1Kπl𝒩D(𝐲𝐜l,𝚪l)𝒩L(𝐱𝐀k𝐲+𝐛k,𝚺k). (3)

A useful feature of GLLiM models is described in the following Lemma 2.1, established for multivariate Gaussian and Student components in Deleforge et al. (2015); Perthame et al. (2018) and which can be straightforwardly extended to Gaussian scale mixtures and elliptical distributions (Nguyen et al., 2022c; Ingrassia et al., 2012).

Lemma 2.1.

The parameter 𝛙 in the forward cPDF, defined in (3), can then be deduced from 𝛙 in (2) via the following one-to-one correspondence:

(πk𝐜k𝚪k𝐀k𝐛k𝚺k)k[K](πk𝐜k𝚪k𝐀k𝐛k𝚺k)k[K]=(πk𝐀k𝐜k+𝐛k𝚺k+𝐀k𝚪k𝐀k𝚺k𝐀k𝚺k1𝚺k(𝚪k1𝐜k𝐀k𝚺k1𝐛k)(𝚪k1+𝐀k𝚺k1𝐀k)1)k[K].
Remark 1 (GLLiM models are computationally efficient).

Without assuming anything about the structure of the parameters, the dimension of GLLiM is dim(𝚿)=K(1+D(L+1)+D(D+1)/2+L(L+1)/2+L)1. It is worth noting that dim(𝚿) can be very large compared to the sample size (see, e.g., Deleforge et al. 2015 for real data sets) whenever D is large and DL. Furthermore, under the assumption that the K transformations are affine, it is more realistic to make the assumption on the residual covariance matrices 𝚺k of the error vectors 𝐄k rather than on 𝚪k (cf.,  Section 3 Deleforge et al., 2015). This justifies using the inverse regression trick from Deleforge et al. (2015), drastically reducing the number of parameters to be estimated. For instance, 𝐄k can be modelled with equal isotropic Gaussian noise, so we have 𝚺k=σ~2ID,k[K], with some positive σ~2. The number of parameters to be estimated is then K+K(L+L(L+1)/2+DL+D). For example, it is 30,060 if K=10, L=2, D=1000. However, if a high-to-low regression is estimated directly instead, the size of the parameter vector will be K+K(D+LD+D(D+1)/2+L), which is 5,035,03030,060 in the previous example.

A notable recent illustration of the GLLiM good features is described in Blein-Nicolas et al. (2024), tackling the challenge of nonlinear quantitative trait prediction using biological data. In this work, the authors focused on predicting a small set of continuous quantitative traits (L=2) from a large set of biomarkers (D1000L). Their inverse regression approach is not only computationally efficient, but also has the advantage of preserving all covariates while assuming a specific structure for the covariance matrix, which aids dimensionality reduction and improves prediction accuracy. This method operates by estimating an inverse model through several low-dimensional regressions and then inverting these estimators to solve the initial high-dimensional regression problem using Lemma 2.1.

2.3 Construction of BNP-GLLiM model

We propose the following hierarchical representation of BNP-GLLiM model to generate a data point (𝐲n,𝐱n) within our BNP-GLLiM model:

  1. 1.

    BNP prior: GBNP(α,σ,G0).

    (α,σ)s1,s2,a𝒮𝒢(ασ;s1,s2)p(σa)Gam(α+σs1,s2)p(σa), (4)
    τkα,σindBeta(τk1σ,α+kσ),k. (5)

    Define πk(𝝉)=τkl=1k1(1τl),k, 𝜽kG0iidG0, where 𝜽k(𝐜k,𝚪k,𝐀k,𝐛k,𝚺k),k, and finally define G=k=1πk(𝝉)δ𝜽k.

  2. 2.

    BNP-GLLiM model: for each n[N], 𝐲nBNP-GLLiM(s1,s2,a,G).

    𝜽nGiidG,if 𝜽n=𝜽k,set zn=k; (6)
    𝐱nzn,𝐜,𝚪iid𝒩L(𝐱n𝐜zn,𝚪zn),(𝐜,𝚪)(𝐜k,𝚪k)k; (7)
    𝐲n𝐱n,zn,𝐀,𝐛,𝚺iid𝒩D(𝐲n𝐀zn𝐱n+𝐛zn,𝚺zn),(𝐀,𝐛,𝚺)(𝐀k,𝐛k,𝚺k)k.

3 Variational inference for BNP-GLLiM

For a brief summary of variational Bayesian expectation maximisation (VBEM) and its notation, see Appendix B. The task of conditional density estimation and clustering using the BNP-GLLiM model is mainly to estimate the unknown labels 𝒵=(𝐳n)n[N] from the observed data (𝒴,𝒳)=(𝐲n,𝐱n)n[N], whose joint distribution p(𝒴,𝒳,𝒵,𝚯;ϕ) is determined by a set of BNP prior parameters 𝚯=(𝝉,α,σ,𝜽), namely the stick-breaking construction of Pitman–Yor process (PYP) (Pitman and Yor, 1997) in Appendix A; and by additional hyperparameters ϕ=(s1,s2,a). Then the desired joint distribution is given by:

p(𝒴,𝒳,𝒵,𝚯;ϕ) =n=1Np(𝐲n𝐱n,zn;𝐀,𝐛,𝚺)p(𝐱nzn,𝐜,𝚪)p(𝒵𝝉)
×kp(τkα,σ)p(α,σs1,s2,a)p(𝜽|G0). (8)

In most variational approximations, the posterior for the stick-breaking variables is approximated in a factorized form (mean-field approximation). Following the same approach, by factorizing the latent variables and the parameters, we choose the following variational distribution: q(𝒵,𝚯)=q𝒵(𝒵)q𝚯(𝚯).

𝝉zn𝐱n𝐲nN𝐜,𝚪αs1,s2σa𝐀,𝐛,𝚺
Figure 1: Graphical representation of BNP-GLLiM: the plate denotes N i.i.d. observations, white-filled circles correspond to latent variables and random or unknown parameters represented in red, while grey-filled circles correspond to observed variables represented in green. Hyperparameters are represented in blue.

In particular, the intractable posterior on 𝒵 is approximated as q𝒵(𝒵) that factorizes so as to handle intractability, namely

q𝒵(𝒵)=n=1Nqzn(zn). (9)

Then the infinite state space for each zj is dealt with by choosing a truncation of the state space to a maximum label K, see, e.g., Blei and Jordan (2006); Wang et al. (2011). In practice, this consists of assuming that the variational distributions qzn for n[N], satisfy qzn(k)=0 for k>K and that the variational distribution on 𝝉 also factorizes as q𝝉(𝝉)=k=1K1qτk(τk) with the additional condition that τK=1. Thus the truncated variational posterior of parameters 𝚯 is given by

q𝚯(𝚯)=qα,σ(α,σ)k=1K1qτk(τk)k=1Kq𝜽k(𝜽k).

In practice, for tractability reasons, we have to restrict to p(𝜽k|G0)1 and q𝜽k(𝜽k) to Dirac distributions, q𝜽k=δ𝜽k, which is equivalent to treating the 𝜽k as fixed unknown hyperparameters, as illustrated graphically in Figure 1. These forms of q𝒵 and q𝚯 lead to our three VB-E steps and four VB-M steps, summarized below and with more detail in Appendix C. Set the initial value of ϕ to ϕ(0). Then repeat the following steps, iteratively. The iteration index is omitted in the update formulas for simplicity. Note that a more complex version with a normal-inverse-Wishart (NIW) distribution on the gating parameters (𝐜k,𝚪k), referred to as BNP-GLLiM2, is presented in Appendix F.

3.1 VB-E steps

VB-E-𝝉 step.

The VB-E-𝝉 step corresponds to a variational approximation in the exponential family case and results in a posterior from the same family as the prior. More precisely, to achieve this, we use (5), (6), (8), and are only interested in the functional dependence of (22) on the variable τk. Thus, any terms that do not depend on τk can be included in the additive normalization constant. Then, given for k[K1] that Nk=n=1Nqzn(k) corresponds to the weight of the cluster k, see more details in Section C.1, it holds that

qτk(τk) =Beta(τkγ^k,1,γ^k,2), where,
γ^k,1 =1𝔼qα,σ[σ]+Nk,γ^k,2=𝔼qα,σ[α]+k𝔼qα,σ[σ]+l=k+1KNl.
VB-E-(α,σ) step.

The (α,σ) variational posterior is more complex, but has a simple form in the DP case (σ=0). Specifically, we have to compute

s^1 =s1+K1,s^2=s2k=1K1ψ(γ^k,2)ψ(γ^k,1+γ^k,2), (10)

where ψ() is the digamma function defined by ψ(γ)=ddγlogΓ(γ)=Γ(γ)Γ(γ).

When σ=0 then qα,σqα,0 is a gamma distribution Gam(s^1,s^2) with 𝔼qα,σ[α]=s^1/s^2. Otherwise (PYP case), qα,σ is only identified up to a normalizing constant but the required 𝔼qα,σ[α] and 𝔼qα,σ[σ] can be computed by importance sampling, see Section C.2 for more details.

We next consider the derivation of the update equation for the factor q𝒵(𝒵).

VB-E-𝒵 step.

By using the mean-field approximation (9) and the truncation (see Section C.3), for all n[N] and for all k[K], this step consists in computing

qzn(k)=ρnkl=1Kρnl, where we define logρnk as follows: (11)
logρnk= 12{Dlog(2π)+log|𝚺^k|+(𝐲n𝐀^k𝐱n𝐛^k)𝚺^k1(𝐲n𝐀^k𝐱n𝐛^k)
+Llog(2π)+log|𝚪^k|+(𝐱n𝐜^k)𝚪^1k(𝐱n𝐜^k)}
+ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1[ψ(γ^l,2)ψ(γ^l,1+γ^l,2)].

Note that in the above formula, symbols (𝐜^k,𝚪^k,𝐀^k,𝐛^k,𝚺^k) are the hyperparameters Specifically defined in the following Sections 3.2 and 3.2. It is important to note that (11) provides assignment probabilities qzn(k) rather than intermediate commitments to hard assignments of zn. However, the hard assignments can be postponed to the end if desired to obtain a segmentation by the following MAP estimation z^n=argmaxk[K]qzn(k).

3.2 VB-M steps

The maximisation step consists of updating the hyperparameters ϕ=(s1,s2,a,𝜽[K]), where 𝜽[K]=(𝐜k,𝚪k,𝐀k,𝐛k,𝚺k)k[K], by maximizing the free energy as follows:

ϕ(r)=argmaxϕ𝔼q𝒵(r)q𝝉(r)qα,σ(r)[logp(𝒴,𝒳,𝒵,𝝉,α,σ;ϕ)].

The VB-M-step can therefore be divided into 4 independent sub-steps, as listed below. From the conditional independence of (s1,s2,a) and (𝒴,𝒳,𝒵) given (𝝉,α,σ), the solution for the VB-M-(s1,s2) (in the DP case) step is straightforward. Only the M-(s1,s2,a) (in the PYP case) and M-𝜽[K] steps are more involved.

VB-M-(s1,s2,a) step.

This step is straightforward in the DP case (σ=0). It can be expressed easily using the fact that both the prior and the variational posterior are Gamma distributions, and using the cross-entropy properties,

(s1,s2)(r)=argmax(s1,s2)𝔼qα,0(r)[logp(αs1,s2)]=(s^(r)1,s^(r)2),

where (s^(r)1,s^(r)2) is given in (10). We can also solve this step numerically using importance sampling in the more general case of PYP (σ0). For more details, see Appendix A.7 in et al. (2020).

VB-M-(𝐜,𝚪) step.

This step is divided into solving K optimisation problems:

(𝐜^k,𝚪^k)(𝐜^k(r),𝚪^k(r)) =argmax(𝐜k,𝚪k)𝔼q𝒵(r)[logp(𝒳𝒵;𝐜k,𝚪k)].

We can then update the Gaussian gating parameters as follows:

𝐜^k=1Nkn=1Nqzn(k)𝐱n,𝚪^k=1Nkn=1Nqzn(k)(𝐱n𝐜^k)(𝐱n𝐜^k).

The technical details will be left to the Section C.4.

VB-M-(𝐀,𝐛,𝚺) step.

Using the same idea, this step is divided into K sub-steps, which include the following optimisation problems

(𝐀^k,𝐛^k,𝚺^k)(𝐀^k(r),𝐛^k(r),𝚺^k(r)) =argmax(𝐀k,𝐛k,𝚺k)𝔼q𝒵(r)[logp(𝒴𝒳,𝒵;𝐀k,𝐛k,𝚺k)].

Given the following quantities:

𝐱¯k =1Nkn=1Nqzn(k)𝐱n,𝑿k=1Nk(qz1(k)(𝐱1𝐱¯k),,qzN(k)(𝐱N𝐱¯k)),
𝐲¯k =1Nkn=1Nqzn(k)𝐲n,𝒀k=1Nk(qz1(k)(𝐲1𝐲¯k),,qzN(k)(𝐲N𝐲¯k)),

We can update the parameters for the Gaussian experts as follows (cf. Section C.5):

𝐀^k =𝒀k𝑿k(𝑿k𝑿k)1,𝐛^k=1Nkn=1Nqzn(k)(𝐲n𝐀^k𝐱n),
𝚺^k =1Nkn=1Nqzn(k)(𝐲n𝐀^k𝐱n𝐛^k)(𝐲n𝐀^k𝐱n𝐛^k).

3.3 Evidence lower-bound (ELBO)

Evaluating the ELBO in (20) allows us to not only monitor the bound during the re-estimation to test for convergence but also to check both the mathematical expressions for the solutions and their software implementation. Indeed, the value of this bound (20) at each step of the iterative re-estimation procedure should not decrease (Svensén and Bishop, 2005), in particular, see recent results for Bayesian nonparametric mixture models in Appendix A of Durand et al. (2022). Recall that ϕ^=(s^1,s^2,a^,(𝐜^k,𝚪^k,𝐀^k,𝐛^k,𝚺^k)k). Here, in order to keep the notation uncluttered, we will sometimes omit the subscripts on the expectation operators because each expectation is taken with respect to all of the random variables in its argument, and the hat superscript ^ on the hyperparameters ϕ^ of q distribution. If σ0 and there are enough training data, the ELBO can be evaluated via the fact that the integral reduces to a point evaluation at the posterior mean of each parameter, see, e.g., Yuan and Neubauer (2008); Luo and Sun (2017); Wu and Ma (2018); Nguyen and Bonilla (2014). When σ=0, we can analytically compute the ELBO in the BNP-GLLiM via Proposition 3.1 which is proved in Section E.2.

Proposition 3.1.

If σ=0, the ELBO in the BNP-GLLiM is decomposed as follows:

(q𝒵,q𝚯,ϕ^) =𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]+𝔼[logp(𝒳𝒵,𝚯;ϕ^)]+𝔼[logp(𝒵𝚯;ϕ^)]
+𝔼[logp(𝚯;ϕ^)]𝔼[logq(𝒵)]𝔼[logq(𝚯)], (12)

where all the terms have a closed-form expression.

The closed-form expressions for the terms of the right-hand side of Equation 12 are provided in Appendix D. Note that if the free energy is computed at the end of each VBEM iteration, as in Section 3.2, we have 𝔼[logqα,0(α)]=𝔼[logp(αs^1,s^2)].

4 Predictive conditional density

The most popular uses of BNP-GLLiM with discrete random probability measures, such as the one displayed in (7), relate to conditional density estimation and data clustering. Specifically, we are interested in the predicted conditional density for a new value (𝐲^,𝐱^) of the observed variables. Note that there will be a corresponding latent variable 𝐳^ associated with these observations. If σ0, we can use the previous remark in Section 3.3, where the integral reduces to a point evaluation at the posterior mean of each parameter. When σ=0, we can analytically approximate such densities via several following theorems. In Theorems 4.1, 4.2 and 4.3, the notation “” means that we approximate the desired densities of the BNP-GLLiM by a mixture of Gaussians using factorized variational approximation posteriors and a truncation of K.

4.1 Joint density

We first compute the joint density via Theorem 4.1, which is proved in Section E.3.

Theorem 4.1.

With 𝐰^[𝐱^;𝐲^], we have

𝝁k𝔼[𝐰^] =(𝐜^k𝐀^k𝐜^k+𝐛^k),𝐕k=cov[𝐰^]=(𝚪^k𝚪^k𝐀^k𝐀^k𝚪^k𝚺^k+𝐀^k𝚪^k𝐀^k), (13)
p(𝐲^,𝐱^,𝒳,𝒴) k=1K𝔼qτk[τk]l=1k1𝔼qτl[1τl]𝒩L+D(𝐰^𝝁k,𝐕k), (14)
𝔼qτk[τk] =γ^k,1γ^k,1+γ^k,2,𝔼qτk[1τk]=γ^k,2γ^k,1+γ^k,2.

4.2 Inverse conditional density

We then show how to approximate the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴). This predictive density in BNP-GLLiM is approximated by a GLLiM via Theorem 4.2 with the proof in Section E.4.

Theorem 4.2.

We approximate the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴) and its conditional expectation for prediction by:

p(𝐲^𝐱^,𝒳,𝒴) k=1KgLk(𝐱^𝚯^,ϕ^,𝒳,𝒴)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k), (15)
𝔼[𝐲^𝐱^,𝒳,𝒴] k=1KgLk(𝐱^𝚯^,ϕ^,𝒳,𝒴)[𝐀^k𝐱^+𝐛^k], where
gLk(𝐱^𝚯^,ϕ^,𝒳,𝒴) =𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)l=1K𝔼q𝝉[πl(𝝉)]𝒩L(𝐱^𝐜^l,𝚪^l),k[K].

4.3 Forward conditional density

Given the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴), we approximate the following forward conditional density via Theorem 4.3, whose proof is provided in Section E.5.

Theorem 4.3.

We approximate the forward conditional density and its corresponding conditional expectation and variance for prediction and uncertainty estimation by

p(𝐱^𝐲^,𝒳,𝒴) k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k),
𝔼[𝐱^𝐲^,𝒳,𝒴] k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)(𝐀^k𝐲^+𝐛^k),
var[𝐱^𝐲^,𝒳,𝒴] k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)[𝚺^k+(𝐀^k𝐲^+𝐛^k)(𝐀^k𝐲^+𝐛^k)]
𝔼(𝐱^𝐲^,𝒳,𝒴)𝔼(𝐱^𝐲^,𝒳,𝒴), where
gDk(𝐲^𝚯^,ϕ^,𝒳,𝒴) =𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k)l=1K𝔼q𝝉[πl(𝝉)]𝒩D(𝐲^𝐜^l,𝚪^l),
𝚺^k =(𝚪^1k+𝐀^k𝚺^1k𝐀^k)1,𝐛^k=𝚺^k[𝚪^1k𝐜^k𝐀^k𝚺^1k𝐛^k],
𝐀^k =𝚺^k𝐀^k𝚺^1k,𝐜^k=𝐀^k𝐜^k+𝐛^k,𝚪^k=𝚺^k+𝐀^k𝚪^k𝐀^k.

5 Bayesian nonparametric model selection

Notations. A coupling between 𝝅(πk)k[K] and 𝝅0(π0l)l[K0] is a joint distribution 𝐐 on [K]×[K0], which is expressed as a matrix 𝐐=(qkl)k[K],l[K0][0,1]K×K0 with marginal probabilities k[K]qkl=π0l and l[K0]qkl=πk, for any k[K] and l[K0]. We use 𝒬(𝝅,𝝅0) to denote the space of all such couplings. Regarding the space of mixing measures, let KK(𝚯) and 𝒪K𝒪K(𝚯) respectively denote the space of all mixing measures with exactly and at most K support points, all in some parameter space 𝚯. Additionally, with 𝒢𝒢(𝚯)=K+K to denote the set of all discrete measures with finite supports on 𝚯. Moreover, 𝒢¯(𝚯) denotes the space of all discrete measures (including those with countably infinite supports) on 𝚯. Finally, 𝒫(𝚯) stands for the space of all probability measures on 𝚯.

5.1 Posterior contraction rate in Bayesian infinite mixtures

Problem setup.

We first recall the GMM where we have i.i.d. samples (𝐖n)n[N]𝒲 coming from a true but unknown distribution PG0 with given PDF

pG0𝒩(𝐰𝜽)dG0(𝜽)=k[K0]πk0𝒩(𝐰𝜽k0),𝜽0k(𝝁0k,𝐕0k), (16)

where G0=k[K0]πk0δ𝜽k0 is a true but unknown mixing distribution with exactly K0 number of support points, where K0 is also unknown. Furthermore, 𝚯 is a chosen parameter space, to which we believe that the true parameters belong. In a well-specified setting, all support points of G0 reside in 𝚯, but this may not be the case in a misspecified setting. In this section, we assume that the GMM is well-specified, i.e.,  the data are i.i.d. samples from the mixture density pG0, where mixing measure G0 has K0 support atoms in compact parameter space 𝚯.

A Bayesian modeller places a prior distribution Π on a suitable subspace of 𝒢¯(𝚯). Then, the posterior distribution over G is given by:

Π(GB𝒲)Bn=1NpG(𝐖n)dΠ(G)𝒢¯(𝚯)n=1NpG(𝐖n)dΠ(G).

Here, the GMM pG is defined in (16) with K unknown number of support points. We are interested in the posterior contraction behaviour of G toward G0, in addition to recovering the true number of components K0.

We next recall the notion of Wasserstein distance for mixing measures that prove useful in the next sections.

Wasserstein distance for MM.

It is useful to analyze the identifiability and convergence of parameter estimation in mixture models using the notion of Wasserstein distance, as in Nguyen (2013); Ho and Nguyen (2016b). This distance can be defined as the optimal cost of moving masses in the transformation from one probability measure to another (Villani, 2003, 2009).

Definition 5.1.

Suppose 𝚯 is equipped with a metric d. The Wasserstein distance Wr between two discrete measures G=k[K]πkδ𝜽k and G0=l[K0]π0lδ𝜽0l is given by

Wr(G,G0)=inf𝐐𝒬(𝝅,𝝅0)[k[K],l[K0]qkl[d(𝜽k,𝜽0l)]r]1/r,

where couplings 𝐐 are defined at the beginning of this section. See Delon and Desolneux (2020) for more details.

It should be emphasized that if a sequence of probability measures GN𝒪K0 converges to G0K0 under the Wr metric at a rate ωN=o(1) for some r1, then there exists a subsequence of GN such that the set of atoms of GN converges to the K0 atoms of G0, up to a permutation of the atoms, at the same rate ωN.

Posterior contraction rate in infinite mixtures.

With a similar idea as in Guha et al. (2021), our starting point is the availability of a mixing measure sample G that is drawn from the posterior distribution Π(G𝒲), where 𝒲 are i.i.d. samples of the mixing density pG0. Under certain conditions on the kernel density, it can be established that for some Wasserstein metric Wr,

Π(G𝒢¯(𝚯):Wr(G,G0)δωN𝒲)N1 in PG0-probability, (17)

for all constant δ>0, while ωN=o(1) is a vanishing rate. Thus ωN can be assumed to be (slightly) slower than the actual rate of posterior contraction of the mixture measure. We can also write that ωN is a rate such that, under the posterior distribution Π(G𝒲), Wr(G,G0)=oPG0(ωN). See Nguyen (2013); Gao and Vaart (2016); Ho and Nguyen (2016b) for concrete examples of posterior contraction rates in infinite and (overfitted) finite mixtures.

5.2 Merge-Truncate-Merge (MTM) algorithm for BNP-GLLiM

Link between GLLiM and joint GMM.

We start by noting that a GLLiM model on (𝐗,𝐘), see (2), with unconstrained parameters 𝝍=(πk,𝐜k,𝚪k,𝐀k,𝐛k,𝚺k)k[K], is equivalent to a GMM on the joint variable [𝐗;𝐘] with unrestricted parameters, via  Lemma 5.2, which is briefly proved in Section E.1.

Lemma 5.2.

A GLLiM model on (𝐗,𝐘) with unconstrained parameters 𝛙=(πk,𝐜k,𝚪k,𝐀k,𝐛k,𝚺k)k[K], defined in (2), is equivalent to a GMM on the joint variable [𝐗;𝐘]𝐖 with unconstrained parameters 𝛎={𝛍k,𝐕k,ρk}k=1K, i.e., 

p(𝐰𝝍)=k=1Kρk𝒩L+D(𝐰𝝁k,𝐕k).

The parameter 𝛙 can be expressed as a function of 𝛎 by:

πk =ρk,𝐜k=𝝁𝐱k,𝚪k=𝐕𝐱𝐱k,𝐀k=𝐕𝐱𝐲k(𝐕𝐱𝐱k)1,
𝐛k =𝝁𝐲k𝐕𝐱𝐲k(𝐕𝐱𝐱k)1𝝁𝐱k,𝚺k=𝐕𝐲𝐲k𝐕𝐱𝐲k(𝐕𝐱𝐱k)1𝐕𝐱𝐲k. (18)

Here, we have defined

𝝁k=[𝝁𝐱k𝝁𝐲k],𝐕k=[𝐕𝐱𝐱k𝐕𝐱𝐲k𝐕𝐲𝐱k𝐕𝐲𝐲k].

Note that the symmetry 𝐕k=𝐕k of the covariance matrix implies that 𝐕𝐱𝐱k and 𝐕𝐲𝐲k are symmetric, while 𝐕𝐱𝐲k=𝐕𝐲𝐱k. The parameter vector 𝛎 can be expressed as a function of 𝛙 by:

ρk =πk,𝝁k=[𝐜k𝐀k𝐜k+𝐛k],𝐕k=[𝚪k(𝐀k𝚪k)𝐀k𝚪k𝚺k+𝐀k𝚪k𝐀k]. (19)
Merge-Truncate-Merge (MTM) algorithm

is a post-processing procedure applied to a posterior sample of the mixing measure G in BNP-MM, essential for achieving posterior contraction rates under the Wasserstein metric (Guha et al., 2021). We propose in Algorithm 1 an algorithm for BNP joint GMM which follows the same steps as the original MTM algorithm for BNP-MM from Guha et al. (2021). This algorithm involves two main stages: the Merge procedure and the Truncate-Merge procedure. In the Merge stage, atoms are reordered by simple random sampling without replacement, ensuring random permutation. Sequentially, atoms are merged based on a distance threshold ωN, updating weights and removing merged atoms from the set, resulting in a new measure G with reordered weights. In the Truncate-Merge stage, atoms are divided into two sets based on a weight threshold (cωN)r. For each atom in the significant weight set, if another atom within the threshold distance exists, it is moved to the negligible weight set. Atoms in the negligible set are then merged with the nearest significant atom, resulting in the final measure G~ and the number of its supporting atoms K~. See Algorithm 1 for a pseudo-code description.

Algorithm 1 MTM Algorithm for BNP joint GMM
1:Posterior sample G=k=1Kπkδ𝜽k, posterior contraction rate ωN from from (17), and a tunning parameter c.
2:Discrete measure G~ and its number of supporting atoms K~. Stage 1: Merge procedure:
3:Reorder atoms {𝜽k}k[K] by simple random sampling without replacement with corresponding weights {π1,π2,}. Let τ1,τ2, denote the new indices, and set ={τj}j as the existing set of atoms.
4:Sequentially for each index τj, if there exists an index τi<τj such that d(𝜽τi,𝜽τj)ωN, then update πτi=πτi+πτj, and remove τj from .
5:Collect G=j:τjπτjδ𝜽τj. Then write G as k=1Kqkδϕk so that q1q2. Stage 2: Truncate-Merge procedure:
6:Set 𝒜={i:qi>(cωN)r}, 𝒩={i:qi(cωN)r}.
7:For each index i𝒜, if there is j𝒜 such that j<i and qiϕiϕjr(cωN)r, then remove i from 𝒜 and add it to 𝒩.
8:For each i𝒩, find atom ϕj among j𝒜 that is nearest to ϕi, then update qj=qj+qi.
9:Return G~=j𝒜qjδϕj and K~=|𝒜|.

As a consequence of our MTM algorithm, we obtain the theoretical guarantee of Theorem 5.3 for the outcome of Algorithm 1.

Theorem 5.3 (MTM consistency for BNP joint GMM).

Let G be a posterior sample from the posterior distribution of any Bayesian procedure, namely, Π(G𝒲) according to which the upper bound (17) holds for all δ>0. Let G~ and K~ be the outcome of Algorithm 1 applied to G, for an arbitrary constant c>0. Then the following hold

  1. (a)

    Π(K~=K0𝒲)N1 in PG0-probability.

  2. (b)

    For all δ>0, Π(G𝒢¯(𝚯):Wr(G~,G0)δωN𝒲)N1 in PG0-probability.

Proof of Theorem 5.3.

Lemma 5.2 implies that BNP joint GMM and BNP-GLLiM are considered equivalent with respect to the number of components in the model selection problem. Therefore, using Theorem 3.2 from Guha et al. (2021) and  Lemma 5.2, it follows that the result of the MTM Algorithm 1 for BNP joint GMM is a consistent estimate of both the number of components and the mixing measure. The latter also admits the upper bound of the posterior contraction rate ωN, which leads to the desired Theorem 5.3. ∎

Remark 2.

Regarding the above theorem, we provide the following comments on posterior consistency for the number of components in BNP-GLLiM after the MTM algorithm post-processing.

  1. (i)

    As a complementary result to Guha et al. (2021), the aim of this paper is to study the practical viability of MTM Algorithm 1 and Theorem 5.3 in the context of high-to-low dimensional inverse regression via BNP-GLLiM model. In order to do this, we first need to specify the metric d in 𝚯, e.g., d(𝛉τi,𝛉τj)=𝛍τi𝛍τj+𝐕τi𝐕τj. Here, denotes either the l2-norm elements in L+D or the entrywise l2-norm for matrices in (L+D)×(L+D).

  2. (ii)

    In practice, one may not have a mixing measure G sampled from the posterior Π(𝒲), but rather a sample of G itself. In particular, to deal with large data sets, we need to use VBEM. Therefore, in BNP-GLLiM, instead, we only obtain a sample FN from the variational posterior GV=k=1K𝔼q𝝉[πk(𝝉)]δ𝜽k. Here, 𝔼q𝝉[πk(𝝉)] and 𝜽k=(𝝁k,𝐕k) are defined in (14) and (13), respectively. However, as long as FN is sufficiently close to G in the sense that Wr(FN,G)Wr(G,G0), we can still apply the MTM algorithm to FN, instead. This requires an extension of the above Theorem 5.3 to cover this scenario and verify this approximation condition, which we leave for future work.

6 Numerical experiments

The code to reproduce our simulation study is publicly available111https://github.com/Trung-TinNGUYEN/BNP-GLLiM and all simulations below were performed in Python 3.9.13 on a standard Unix machine. For proof-of-concept numerical experiments, we consider only the simple data generating mechanism with D=1 and L=1 and demonstrate that BNP-GLLiM performs well in model selection, clustering and cPDF estimation with the MTM procedure. Real world examples are postponed to future work.

6.1 Data generation

We illustrate our theoretical results on simulated datasets in a more general setting for the BNP approach compared to those considered by Chamroukhi et al. (2010); Montuelle and Le Pennec (2014); Nguyen et al. (2022c). Specifically, we consider the following true inverse cPDF from GLLiM model as follows:

s0(𝐲𝐱;𝝍0)=k=1K0π0k𝒩L(𝐱𝐜0k,𝚪0k)l=1Kπ0l𝒩L(𝐱𝐜0l,𝚪0l)𝒩D(𝐲𝐀0k𝐱+𝐛0k,𝚺0k).

Here K0=3, L=D=1, and 𝝍0=(𝝅0,𝐜0,𝚪0,𝐀0,𝐛0,𝚺0), where

𝝅0 =(0.3,0.4,0.3),𝐜0=(1,0.05,1),𝚪0=(0.1,0.2,0.1),
𝐀0 =(15,3,15),𝐛0=(2,1,2),𝚺0=(0.5,0.3,0.5).

Figure 2a shows typical N=1,000 realisations of the true inverse cPDF from GLLiM, representing a π-shape simulation with three clusters without labels.

6.2 Model selection, clustering and regression tasks via MTM-BNP-GLLiM

Our goal is to evaluate the inverse and forward cPDF, as well as the conditional means, to investigate the empirical performance of our MTM-BNP-GLLiM in the previous simulation. In Figure 2 it is clear that with the help of MTM Algorithm 1, MTM-BNP-GLLiM can simultaneously perform regression, clustering and model selection well. Without the MTM procedure, BNP-GLLiM performs poorly in model selection, clustering and cPDF estimation, except for conditional expectations as shown in Figure 3.

Refer to caption
(a) Typical realisations without labels.
Refer to caption
(b) True conditional expectations with labels.
Refer to caption
(c) True and estimated inverse means.
Refer to caption
(d) True and estimated forward means.
Refer to caption
(e) Contour of estimated inverse cPDF.
Refer to caption
(f) Contour of estimated forward cPDF.
Refer to caption
(g) Contour of true inverse cPDF.
Refer to caption
(h) Contour of true forward cPDF.
Figure 2: Top row: Typical 1, 000 realisations of GLLiM’s true inverse cPDF with its true conditional expectations. 2nd row and below: True and estimated inverse and forward cPDF of GLLiM with the true number of components (KMTM=3) using MTM algorithm for post-processing in MTM-BNP-GLLiM.
Refer to caption
(a) True and estimated inverse means.
Refer to caption
(b) True and estimated forward means.
Refer to caption
(c) Contour of estimated inverse cPDF.
Refer to caption
(d) Contour of estimated forward cPDF.
Figure 3: True and estimated inverse and forward means and CPDFs of GLLiM without MTM algorithm for post-processing in BNP-GLLiM with truncated number of components K=20.

Next, we illustrate the performance of the MTM algorithm when applied to the variational posterior from BNP-GLLiM. Specifically, the samples in our 100 trials are drawn from GV=k=1K𝔼q𝝉[πk(𝝉)]δ𝜽k, wherere 𝔼q𝝉[πk(𝝉)] and 𝜽k=(𝝁k,𝐕k) are defined in (14) and (13), respectively. We know that for some constant C~, which depends on the covariance matrix 𝐕k0, the location parameters 𝝁0k and the weights πk0, the contraction rate of mixing measures under the location Gaussian DP-MM is C~(log(N))1/2 with respect to the W2-norm. Similar to Guha et al. (2021), our first attempt to choose wN to satisfy (17) is wN=(log(log(N))log(N))1/2. In fact, we can choose any wN, as long as wNlog(N)1/2, in order for wN to satisfy (17).

Since we only work with finite sample N, it is not expected that the posterior probability for KMTM=K0 is close to 1 and the input c to Algorithm 1 should be chosen so that C~(log(log(N)))1/2c. Furthermore, based on Equation (26) from Guha et al. (2021), with a useful lower bound on the posterior mass the mode, for any 1>ϵ>0, (1ϵ)(1k=13c0r/2πk0)>12, we hope to identify K0 via the posterior mode with a reasonable estimate. To guarantee K=K0 consistently using the posterior mode safety, we have to choose c<c0, with c0 satisfying

(1ϵ)(1k=13c0r/2πk0)>1212ϵ2(1ϵ)(k=131πk0)1>c0r/2>cr/2.

Therefore, we can choose

[12ϵ4(1ϵ)(k=131πk0)1]2/r=c0>c, for all 12>ϵ>0.

In particular, it is unrealistic to obtain the exact computation of the upper bound c0 and the lower bound C~(log(log(N)))1/2. However, a reasonable estimate may be possible by considering a large range of c, and show that there is a range where we can robustly identify the true number of components via the posterior mode. Guha et al. (2021) also used the same setting in their experiments. Figure 4 indicates that c=0.45 leads to a quite good posterior mode in our experiments.

Although we do not have a theoretical result for the convergence rate of the variational posterior of BNP-GLLiM to the true data generating process, Figure 4 seems to suggest that MTM-BNP-GLLiM gives a comparable good result to the location Gaussian DP-MM in the simulation studies in Guha et al. (2021).

Refer to caption
Figure 4: Histogram of KMTM using variational posterior sample with 100 trials and c=0.45.

7 Perspectives

To address the dimensionality issue when N is less than D, we could incorporate sparsity penalty terms as described in the GLLiM context in Chamroukhi et al. (2019). Alternatively, we could impose additional structural restrictions, such as block-diagonal covariance matrices as in Blein-Nicolas et al. (2024)), which extends the GLLiM method to account for hidden module-structured regulatory networks of predictors through block-diagonal structures. In particular, Blein-Nicolas et al. (2024) applied their method to the prediction of drought-related traits (L=2) from protein abundance (D=973) in N=233 maize genotypes. We leave for future research the intriguing but challenging questions of how to integrate BNP priors into the GLLiM model and how to establish posterior convergence theory in high-dimensional settings where DN.

As indicated in Remark 2, there is a crucial need to formally establish general conditions on the prior, the likelihood and the variational class to characterise the convergence rate of the variational posterior of BNP-GLLiM to the true data generating process. Using the similar “prior mass and testing” conditions as in Ghosal et al. (2000), we believe that an interesting but challenging extension of the work on variational posterior unconditional distributions for MMs (Zhang and Gao, 2020) and on adaptive Bayesian estimation for MMs and MoE models but for true posterior distribution (Kruijer et al., 2010; Shen et al., 2013; Norets and Pati, 2017) can help shed some light and answer this important question. Furthermore, it is important to establish an extensional convergence property of our VBEM algorithm for BNP-GLLiM. This property is only known for GMM from Titterington and Wang (2006). A potential improvement of the VBEM algorithm developed for BNP-GLLiM can be achieved by combining it with MCMC, taking advantage of both inference approaches as in Ruiz and Titsias (2019). Finally, as mentioned in Section 6.2, the selection of a good data-driven tuning parameter c as the same idea from the slope heuristic of Birgé and Massart (2007) is crucial for the success of the MTM procedure for any BNP model. We leave these interesting but challenging questions for future research.

Acknowledgments

We would like to express our sincere gratitude to the Reviewers and Editors for their valuable feedback, which has greatly enhanced the quality of this paper. All authors acknowledge funding from the Australian Research Council grant DP230100905, and from Inria Project WOMBAT.

References

  • H. Akaike (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control 19 (6), pp. 716–723. External Links: Document Cited by: §1. a
  • L. Alamichel, D. Bystrova, J. Arbel, and G. Kon Kam King (2024) Bayesian mixture models (in)consistency for the number of clusters. Scandinavian Journal of Statistics. Cited by: §1. a
  • J. Arbel, G. Kon Kam King, A. Lijoi, L. Nieto-Barajas, and I. Prünster (2021) BNPdensity: Bayesian nonparametric mixture modelling in R. Australian & New Zealand Journal of Statistics 63 (3), pp. 542–564. External Links: Link, Document Cited by: §1. a
  • S. Arlot (2019) Minimal penalties and the slope heuristics: a survey. Journal de la Société Française de Statistique 160 (3), pp. 1–106. Cited by: §1. a
  • F. Ascolani, A. Lijoi, G. Rebaudo, and G. Zanella (2022) Clustering consistency with Dirichlet process mixtures. Biometrika 110 (2), pp. 551–558. External Links: ISSN 1464-3510, Link, Document Cited by: §1. a
  • M. J. Beal and Z. Ghahramani (2003) The variational Bayesian EM algorithm for incomplete data: with application to scoring graphical model structures. Bayesian statistics 7 (453-464), pp. 210. Cited by: Appendix B. a
  • M. J. Beal (2003) Variational algorithms for approximate Bayesian inference. Ph.D. Thesis, University of London, University College London (United Kingdom). Cited by: Appendix B, Appendix B. a b
  • M. Berrettini, G. Galimberti, S. Ranciati, and T. B. Murphy (2024) Identifying Brexit voting patterns in the British house of commons: an analysis based on Bayesian mixture models with flexible concomitant covariate effects. Journal of the Royal Statistical Society Series C: Applied Statistics 73 (3), pp. 621–638. External Links: ISSN 0035-9254, Link, Document Cited by: §1. a
  • C. Biernacki, G. Celeux, and G. Govaert (2000) Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22 (7), pp. 719–725. External Links: Document Cited by: §1. a
  • L. Birgé and P. Massart (2007) Minimal penalties for Gaussian model selection. Probability Theory and Related Fields 138 (1), pp. 33–73. External Links: Link, Document Cited by: §1, §7. a b
  • C. M. Bishop (2006) Pattern Recognition and Machine Learning (Information Science and Statistics). Springer-Verlag, Berlin, Heidelberg. External Links: ISBN 0-387-31073-8 Cited by: Appendix B, §E.3, §E.5, §F.1.2, §F.4, §F.6. a b c d e f
  • D. M. Blei and M. I. Jordan (2006) Variational inference for Dirichlet process mixtures. Bayesian analysis 1 (1), pp. 121–143. Cited by: Appendix A, §E.3, §E.3, §F.4, §3. a b c d e
  • M. Blein-Nicolas, E. Devijver, M. Gallopin, and E. Perthame (2024) Nonlinear network-based quantitative trait prediction from biological data. Journal of the Royal Statistical Society Series C: Applied Statistics, pp. qlae012. Cited by: §1, §7, Remark 1. a b c d
  • F. Boux, F. Forbes, J. Arbel, B. Lemasson, and E. L. Barbier (2021) Bayesian inverse regression for vascular magnetic resonance fingerprinting. IEEE Transactions on Medical Imaging 40 (7), pp. 1827–1837. Cited by: §1. a
  • G. Celeux, S. Frühwirth-Schnatter, and C. P. Robert (2019) Model selection for mixture models–perspectives and strategies. In Handbook of mixture analysis, pp. 117–154. Cited by: §1. a
  • L. Chaari, T. Vincent, F. Forbes, M. Dojat, and P. Ciuciu (2013) Fast joint detection-estimation of evoked brain activity in event-related fMRI using a variational approach. IEEE transactions on Medical Imaging 32 (5), pp. 821–837. Cited by: Appendix B. a
  • F. Chamroukhi, F. Lecocq, and H. D. Nguyen (2019) Regularized Estimation and Feature Selection in Mixtures of Gaussian-Gated Experts Models. In Statistics and Data Science, H. Nguyen (Ed.), Singapore, pp. 42–56. External Links: ISBN 978-981-15-1960-4 Cited by: §7. a
  • F. Chamroukhi, A. Samé, G. Govaert, and P. Aknin (2010) A hidden process regression model for functional data description. Application to curve discrimination. Neurocomputing 73 (7-9), pp. 1210–1221. Cited by: §6.1. a
  • Z. Chen, Y. Deng, Y. Wu, Q. Gu, and Y. Li (2022) Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In Advances in Neural Information Processing Systems, A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho (Eds.), External Links: Link Cited by: §1. a
  • P. De Blasi, S. Favaro, A. Lijoi, R. H. Mena, I. Prünster, and M. Ruggiero (2015) Are Gibbs-type priors the most natural generalization of the Dirichlet process?. IEEE transactions on pattern analysis and machine intelligence 37 (2), pp. 212–229. Cited by: Appendix A. a
  • A. Deleforge, F. Forbes, and R. Horaud (2015) High-dimensional regression with gaussian mixtures and partially-latent response variables. Statistics and Computing 25 (5), pp. 893–911. External Links: ISSN 1573-1375, Link, Document Cited by: §E.1, §1, §2.1, §2.2, §2.2, Remark 1. a b c d e f g h i
  • J. Delon and A. Desolneux (2020) A Wasserstein-type distance in the space of Gaussian mixture models. SIAM Journal on Imaging Sciences 13 (2), pp. 936–970. Cited by: Definition 5.1. a
  • A. P. Dempster, N. M. Laird, and D. B. Rubin (1977) Maximum Likelihood from Incomplete Data Via the EM Algorithm. Journal of the Royal Statistical Society: Series B (Methodological) 39 (1), pp. 1–22. External Links: ISSN 0035-9246, Link, Document Cited by: Appendix B. a
  • T. G. Do, H. K. Le, T. Nguyen, Q. Pham, B. T. Nguyen, T. Doan, C. Liu, S. Ramasamy, X. Li, and S. HOI (2023) HyperRouter: Towards Efficient Training and Inference of Sparse Mixture of Experts. In Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing, Singapore. Cited by: §1. a
  • J.-B. Durand, F. Forbes, C.D. Phan, L. Truong, H.D. Nguyen, and F. Dama (2022) Bayesian non-parametric spatial prior for traffic crash risk mapping: A case study of Victoria, Australia. Australian & New Zealand Journal of Statistics 64 (2), pp. 171–204. External Links: Link, Document Cited by: Appendix A, §1, §3.3. a b c
  • M. D. Escobar and M. West (1995) Bayesian Density Estimation and Inference Using Mixtures. Journal of the American Statistical Association 90 (430), pp. 577–588. External Links: Link, Document Cited by: §1. a
  • S. Favaro, A. Lijoi, C. Nava, B. Nipoti, I. Pruenster, and Y. W. Teh (2016) On the stick-breaking representation for homogeneous NRMIs. Bayesian Analysis 11 (3), pp. 697–724. Cited by: Appendix A. a
  • T. S. Ferguson (1973) A Bayesian Analysis of Some Nonparametric Problems. The Annals of Statistics 1 (2), pp. 209 – 230. External Links: Link, Document Cited by: Appendix A. a
  • F. Forbes, A. Arnaud, B. Lemasson, and E. Barbier (2019) Component elimination strategies to fit mixtures of multiple scale distributions. In RSSDS 2019 - Research School on Statistics and Data Science, Communications in Computer and Information Science, Vol. 1150, Melbourne, Australia, pp. 81–95. External Links: Link, Document Cited by: §1. a
  • F. Forbes, H. D. Nguyen, T. Nguyen, and J. Arbel (2022a) Mixture of expert posterior surrogates for approximate Bayesian computation. In JDS 2022 - 53èmes Journées de Statistique de la Société Française de Statistique (SFdS), Lyon, France. Cited by: §1. a
  • F. Forbes, H. D. Nguyen, T. Nguyen, and J. Arbel (2022b) Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. Statistics and Computing 32 (5), pp. 85. External Links: ISSN 1573-1375, Link, Document Cited by: §1. a
  • R. Foygel and M. Drton (2010) Extended Bayesian Information Criteria for Gaussian Graphical Models. In Advances in Neural Information Processing Systems, J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta (Eds.), Vol. 23. External Links: Link Cited by: §1. a
  • S. Frühwirth-Schnatter, D. Hosszejni, and H. F. Lopes (2024) Sparse Bayesian Factor Analysis When the Number of Factors Is Unknown. Bayesian Analysis, pp. 1 – 48. Note: Publisher: International Society for Bayesian Analysis External Links: Link, Document Cited by: §1. a
  • S. Frühwirth-Schnatter, C. Pamminger, A. Weber, and R. Winter-Ebmer (2012) Labor market entry and earnings dynamics: Bayesian inference using mixtures-of-experts Markov chain clustering. Journal of Applied Econometrics 27 (7), pp. 1116–1137. External Links: Link, Document Cited by: §1. a
  • S. Frühwirth-Schnatter, S. Pittner, A. Weber, and R. Winter-Ebmer (2018) Analysing plant closure effects using time-varying mixture-of-experts Markov chain clustering. The Annals of Applied Statistics 12 (3), pp. 1796 – 1830. Note: Publisher: Institute of Mathematical Statistics External Links: Link, Document Cited by: §1. a
  • S. Frühwirth-Schnatter (2019) Keeping the balance—Bridge sampling for marginal likelihood estimation in finite mixture, mixture of experts and Markov mixture models. Brazilian Journal of Probability and Statistics 33 (4), pp. 706 – 733. External Links: Link, Document Cited by: §1. a
  • F. Gao and A. v. d. Vaart (2016) Posterior contraction rates for deconvolution of Dirichlet-Laplace mixtures. Electronic Journal of Statistics 10 (1), pp. 608 – 627. External Links: Link, Document Cited by: §5.1. a
  • C. R. Genovese and L. Wasserman (2000) Rates of convergence for the Gaussian mixture sieve. The Annals of Statistics 28 (4), pp. 1105 – 1127. External Links: Link, Document Cited by: §1. a
  • S. Ghosal, J. K. Ghosh, and A. W. v. d. Vaart (2000) Convergence rates of posterior distributions. The Annals of Statistics 28 (2), pp. 500 – 531. External Links: Link, Document Cited by: §7. a
  • S. Ghosal and A. Van der Vaart (2017) Fundamentals of nonparametric Bayesian inference. Vol. 44, Cambridge University Press. External Links: ISBN 0-521-87826-8 Cited by: Appendix A, §E.3, §1. a b c
  • I. C. Gormley and S. Frühwirth-Schnatter (2019) Mixture of experts models. In Handbook of mixture analysis, pp. 271–307. Cited by: §1. a
  • P. J. Green and S. Richardson (2001) Modelling Heterogeneity With and Without the Dirichlet Process. Scandinavian Journal of Statistics 28 (2), pp. 355–375. External Links: Link, Document Cited by: §1. a
  • A. Guha, N. Ho, and X. Nguyen (2021) On posterior contraction of parameters and interpretability in Bayesian mixture modeling. Bernoulli 27 (4), pp. 2159 – 2188. External Links: Link, Document Cited by: §1, §1, item i, §5.1, §5.2, §5.2, §6.2, §6.2, §6.2. a b c d e f g h i j k
  • N. L. Hjort, C. Holmes, P. Müller, and S. G. Walker (Eds.) (2010) Bayesian Nonparametrics. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press. External Links: ISBN 978-0-521-51346-3, Link, Document Cited by: §1. a
  • N. Ho and X. Nguyen (2016a) Convergence rates of parameter estimation for some weakly identifiable finite mixtures. The Annals of Statistics 44 (6), pp. 2726 – 2755. External Links: Link, Document Cited by: §1. a
  • N. Ho and X. Nguyen (2016b) On strong identifiability and convergence rates of parameter estimation in finite mixtures. Electronic Journal of Statistics 10 (1), pp. 271–307. Cited by: §1, §5.1, §5.1. a b c
  • N. Ho, C. Yang, and M. I. Jordan (2022) Convergence Rates for Gaussian Mixtures of Experts. Journal of Machine Learning Research 23 (323), pp. 1–81. External Links: Link Cited by: §1. a
  • B. Hoadley (1970) A Bayesian Look at Inverse Linear Regression. Journal of the American Statistical Association 65 (329), pp. 356–369. External Links: Link, Document Cited by: §1. a
  • S. Ingrassia, S. C. Minotti, and G. Vittadini (2012) Local Statistical Modeling via a Cluster-Weighted Approach with Elliptical Distributions. Journal of Classification 29 (3), pp. 363–401. External Links: ISSN 1432-1343, Link, Document Cited by: §2.2. a
  • H. Ishwaran and L. F. James (2001) Gibbs Sampling Methods for Stick-Breaking Priors. Journal of the American Statistical Association 96 (453), pp. 161–173. External Links: ISSN 0162-1459, Link, Document Cited by: Appendix A. a
  • R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton (1991) Adaptive mixtures of local experts. Neural computation 3 (1), pp. 79–87. Cited by: §1. a
  • W. Jiang and M. A. Tanner (1999) Hierarchical mixtures-of-experts for exponential family regression models: approximation and maximum likelihood estimation. Annals of Statistics, pp. 987–1011. Cited by: §1. a
  • L. Kock, N. Klein, and D. J. Nott (2022) Variational inference and sparsity in high-dimensional deep Gaussian mixture models. Statistics and Computing 32 (5), pp. 70. External Links: ISSN 1573-1375, Link, Document Cited by: §1. a
  • W. Kruijer, J. Rousseau, and A. v. d. Vaart (2010) Adaptive Bayesian density estimation with location-scale mixtures. Electronic Journal of Statistics 4 (none), pp. 1225 – 1257. External Links: Link, Document Cited by: §7. a
  • B. Kugler, F. Forbes, and S. Douté (2022) Fast Bayesian Inversion for high dimensional inverse problems. Statistics and Computing 32 (2). Cited by: §1. a
  • S. Lathuilière, R. Juge, P. Mesejo, R. Muñoz-Salinas, and R. Horaud (2017) Deep mixture of linear inverse regressions applied to head-pose estimation. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 4817–4825. Cited by: §1. a
  • K. Li (1991) Sliced Inverse Regression for Dimension Reduction. Journal of the American Statistical Association 86 (414), pp. 316–327. External Links: ISSN 01621459, Link, Document Cited by: §1, §2.1. a b
  • N. Li, W. Li, Y. Jiang, and S. Xia (2022) Deep Dirichlet process mixture models. In Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence, J. Cussens and K. Zhang (Eds.), Proceedings of Machine Learning Research, Vol. 180, pp. 1138–1147. External Links: Link Cited by: §1. a
  • H. Lü, J. Arbel, and F. Forbes (2020) Bayesian nonparametric priors for hidden Markov random fields. Statistics and Computing 30 (4), pp. 1015–1035. External Links: ISSN 1573-1375, Link, Document Cited by: Appendix A, Appendix A, §3.2. a b c
  • C. Luo and S. Sun (2017) Variational Mixtures of Gaussian Processes for Classification. In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17, pp. 4603–4609. External Links: Link, Document Cited by: §3.3. a
  • S. N. Maceachern and P. Müller (1998) Estimating Mixture of Dirichlet Process Models. Journal of Computational and Graphical Statistics 7 (2), pp. 223–238. External Links: Link, Document Cited by: §1, §1. a b
  • G. J. McLachlan and D. Peel (2000) Finite Mixture Models. John Wiley & Sons. Cited by: §1. a
  • G. J. McLachlan and T. Krishnan (1997) The EM algorithm and extensions. Wiley, New York, USA. Cited by: Appendix B. a
  • J. W. Miller and M. T. Harrison (2014) Inconsistency of Pitman-Yor Process Mixtures for the Number of Components. Journal of Machine Learning Research 15 (96), pp. 3333–3370. External Links: Link Cited by: §1. a
  • L. Montuelle and E. Le Pennec (2014) Mixture of Gaussian regressions model with logistic weights, a penalized maximum likelihood approach. Electronic Journal of Statistics 8 (1), pp. 1661–1695. Cited by: §6.1. a
  • R. M. Neal and G. E. Hinton (1998) A View of the EM Algorithm that Justifies Incremental, Sparse, and other Variants. In Learning in Graphical Models, M. I. Jordan (Ed.), pp. 355–368. External Links: ISBN 978-94-011-5014-9, Link, Document Cited by: Appendix B. a
  • R. M. Neal (2000) Markov Chain Sampling Methods for Dirichlet Process Mixture Models. Journal of Computational and Graphical Statistics 9 (2), pp. 249–265. External Links: Link, Document Cited by: §1. a
  • D. T. Nguyen, S. Jacquemoud, A. Lucas, S. Douté, C. Ferrari, S. Coustance, S. Marcq, and A. Meygret (2024a) Unveiling the Characteristics of the Lunar Surface by Massive Inversion of a Photometric Model. In LPI Contributions, LPI Contributions, Vol. 3040, pp. 1998. Cited by: §1. a
  • D. N. Nguyen and Z. Li (2024) Joint learning of Gaussian graphical models in heterogeneous dependencies of high-dimensional transcriptomic data. In The 16th Asian Conference on Machine Learning (Conference Track), External Links: Link Cited by: §1. a
  • H. D. Nguyen, F. Chamroukhi, and F. Forbes (2019) Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. Neurocomputing 366, pp. 208–214. External Links: ISSN 0925-2312, Link, Document Cited by: §1. a
  • H. D. Nguyen and F. Chamroukhi (2018) Practical and theoretical aspects of mixture-of-experts modeling: An overview. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 8 (4), pp. e1246. Cited by: §1. a
  • H. D. Nguyen, L. R. Lloyd-Jones, and G. J. McLachlan (2016) A universal approximation theorem for mixture-of-experts models. Neural computation 28 (12), pp. 2585–2593. Cited by: §1. a
  • H. D. Nguyen, D. Fryer, and G. J. McLachlan (2022a) Order selection with confidence for finite mixture models. Journal of the Korean Statistical Society. External Links: ISSN 2005-2863, Link, Document Cited by: §1. a
  • H. D. Nguyen, T. Nguyen, F. Chamroukhi, and G. J. McLachlan (2021a) Approximations of conditional probability density functions in Lebesgue spaces via mixture of experts models. Journal of Statistical Distributions and Applications 8 (1), pp. 13. External Links: ISSN 2195-5832, Link, Document Cited by: §1. a
  • H. Nguyen, P. Akbarian, T. Nguyen, and N. Ho (2024b) A General Theory for Softmax Gating Multinomial Logistic Mixture of Experts. In Proceedings of The 41st International Conference on Machine Learning, Cited by: §1. a
  • H. Nguyen, T. Nguyen, and N. Ho (2023a) Demystifying Softmax Gating in Gaussian Mixture of Experts. In Advances in Neural Information Processing Systems, Cited by: §1. a
  • H. Nguyen, T. Nguyen, K. Nguyen, and N. Ho (2024c) Towards Convergence Rates for Parameter Estimation in Gaussian-gated Mixture of Experts. In Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, S. Dasgupta, S. Mandt, and Y. Li (Eds.), Proceedings of Machine Learning Research, Vol. 238, pp. 2683–2691. External Links: Link Cited by: §1. a
  • T. Nguyen and E. Bonilla (2014) Fast Allocation of Gaussian Process Experts. In Proceedings of the 31st International Conference on Machine Learning, E. P. Xing and T. Jebara (Eds.), Proceedings of Machine Learning Research, Vol. 32, Bejing, China, pp. 145–153. External Links: Link Cited by: §3.3. a
  • T. Nguyen, F. Chamroukhi, H. D. Nguyen, and G. J. McLachlan (2023b) Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. Communications in Statistics - Theory and Methods 52 (14), pp. 5048–5059. External Links: Link, Document Cited by: §1. a
  • T. Nguyen, F. Chamroukhi, H. D. Nguyen, and F. Forbes (2021b) Non-asymptotic model selection in block-diagonal mixture of polynomial experts models. Preprint. arXiv:2104.08959. External Links: Link Cited by: §1. a
  • T. Nguyen, F. Chamroukhi, H. D. Nguyen, and F. Forbes (2022b) Model selection by penalization in mixture of experts models with a non-asymptotic approach. In JDS 2022 - 53èmes Journées de Statistique de la Société Française de Statistique (SFdS), Lyon, France. Cited by: §1. a
  • T. Nguyen, D. N. Nguyen, H. D. Nguyen, and F. Chamroukhi (2023c) A non-asymptotic theory for model selection in high-dimensional mixture of experts via joint rank and variable selection. In AJCAI Australasian Joint Conference on Artificial Intelligence 2023, External Links: Link Cited by: §1. a
  • T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan (2023d) Non-asymptotic oracle inequalities for the Lasso in high-dimensional mixture of experts. arXiv:2009.10622. External Links: Link Cited by: §1. a
  • T. Nguyen, H. D. Nguyen, F. Chamroukhi, and F. Forbes (2022c) A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models. Electronic Journal of Statistics 16 (2), pp. 4742 – 4822. External Links: Link, Document Cited by: §1, §2.2, §6.1. a b c
  • T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan (2020) Approximation by finite mixtures of continuous density functions that vanish at infinity. Cogent Mathematics & Statistics 7 (1), pp. 1750861. External Links: ISSN null, Link, Document Cited by: §1. a
  • T. Nguyen (2021) Model Selection and Approximation in High-dimensional Mixtures of Experts Models: from Theory to Practice. PhD Thesis, Normandie Université, (en). External Links: Link Cited by: §1. a
  • X. Nguyen (2013) Convergence of latent mixing measures in finite and infinite mixture models. The Annals of Statistics 41 (1), pp. 370–400. External Links: Link, Document Cited by: §1, §5.1, §5.1. a b c
  • A. Norets and D. Pati (2017) Adaptive Bayesian estimation of conditional densities. Econometric Theory 33 (4), pp. 980–1012. External Links: Document Cited by: §7. a
  • A. Norets (2010) Approximation of conditional densities by smooth mixtures of regressions. The Annals of Statistics 38 (3), pp. 1733 – 1766. External Links: Link, Document Cited by: §1. a
  • E. Perthame, F. Forbes, and A. Deleforge (2018) Inverse regression approach to robust nonlinear high-to-low dimensional mapping. Journal of Multivariate Analysis 163 (C), pp. 1–14. Cited by: §2.2. a
  • J. Pitman and M. Yor (1997) The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator. The Annals of Probability 25 (2), pp. 855–900. External Links: Document Cited by: Appendix A, Appendix A, §3. a b c
  • A. Rakhlin, D. Panchenko, and S. Mukherjee (2005) Risk bounds for mixture density estimation. ESAIM: PS 9, pp. 220–229. External Links: Link, Document Cited by: §1. a
  • J. Rousseau and K. Mengersen (2011) Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (5), pp. 689–710. External Links: Link, Document Cited by: §1. a
  • F. Ruiz and M. Titsias (2019) A Contrastive Divergence for Combining Variational Inference and MCMC. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, pp. 5537–5545. External Links: Link Cited by: §7. a
  • G. Schwarz (1978) Estimating the dimension of a model. The Annals of Statistics 6 (2), pp. 461–464. Cited by: §1. a
  • J. Sethuraman (1994) A constructive definition of Dirichlet priors. Statistica Sinica 4, pp. 639–650. Cited by: Appendix A. a
  • W. Shen, S. T. Tokdar, and S. Ghosal (2013) Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. Biometrika 100 (3), pp. 623–640. Cited by: §1, §7. a b
  • C. Sin and H. White (1996) Information criteria for selecting possibly misspecified parametric models. Journal of Econometrics 71 (1), pp. 207–225. External Links: ISSN 0304-4076, Link, Document Cited by: §1. a
  • M. Svensén and C. M. Bishop (2005) Robust Bayesian mixture modelling. Neurocomputing 64, pp. 235–252. External Links: ISSN 0925-2312, Link, Document Cited by: §3.3. a
  • D. M. Titterington and B. Wang (2006) Convergence properties of a general algorithm for calculating variational Bayesian estimates for a normal mixture model. Bayesian Analysis 1 (3), pp. 625 – 650. External Links: Link, Document Cited by: §7. a
  • C. Villani (2003) Topics in optimal transportation. Graduate Studies in Mathematics, Vol. 58, American Mathematical Society. External Links: ISBN 0-8218-3312-X Cited by: §5.1. a
  • C. Villani (2009) Optimal transport: old and new. Vol. 338, Springer. Cited by: §5.1. a
  • C. Wang, J. Paisley, and D. M. Blei (2011) Online Variational Inference for the Hierarchical Dirichlet Process. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, G. Gordon, D. Dunson, and M. Dudík (Eds.), Proceedings of Machine Learning Research, Vol. 15, Fort Lauderdale, FL, USA, pp. 752–760. External Links: Link Cited by: §3. a
  • J. Westerhout, T. Nguyen, X. Guo, and H. D. Nguyen (2024) On the Asymptotic Distribution of the Minimum Empirical Risk. In Forty-first International Conference on Machine Learning, External Links: Link Cited by: §1. a
  • D. Wu and J. Ma (2018) A Two-Layer Mixture Model of Gaussian Process Functional Regressions and Its MCMC EM Algorithm. IEEE Transactions on Neural Networks and Learning Systems 29 (10), pp. 4894–4904. External Links: Document Cited by: §3.3. a
  • L. Xu, M. Jordan, and G. E. Hinton (1995) An Alternative Model for Mixtures of Experts. In Advances in Neural Information Processing Systems, G. Tesauro, D. Touretzky, and T. Leen (Eds.), Vol. 7. External Links: Link Cited by: §1. a
  • C. Yuan and C. Neubauer (2008) Variational Mixture of Gaussian Process Experts. In Advances in Neural Information Processing Systems, D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou (Eds.), Vol. 21. External Links: Link Cited by: §3.3. a
  • S. E. Yuksel, J. N. Wilson, and P. D. Gader (2012) Twenty Years of Mixture of Experts. IEEE Transactions on Neural Networks and Learning Systems 23 (8), pp. 1177–1193. External Links: ISSN 2162-2388 VO - 23, Document Cited by: §1. a
  • G. Zens (2019) Bayesian shrinkage in mixture-of-experts models: identifying robust determinants of class membership. Advances in Data Analysis and Classification 13 (4), pp. 1019–1051. External Links: ISSN 1862-5355, Link, Document Cited by: §1. a
  • F. Zhang and C. Gao (2020) Convergence rates of variational posterior distributions. The Annals of Statistics 48 (4), pp. 2180 – 2207. External Links: Link, Document Cited by: §7. a

Supplementary Materials for
“Bayesian nonparametric mixture of experts for inverse problems”

In this supplementary material, we first recall the standard Bayesian nonparametric priors and variational Bayesian expectation-maximisation principle in Appendices A and B, respectively. Then, all specifications of the VBEM for the BNP-GLLiM model, the evidence lower bound and the technical proofs that are not included in the main paper are placed in Appendices C, D and E. . Finally, Appendix F proposes a more general model with a hyper prior on the gating parameters, referred to as BNP-GLLiM2.

Appendix A Bayesian nonparametric priors

Stick-breaking construction of Dirichlet process.

Note that the Dirichlet process (DP) (Ferguson, 1973) is a central BNP prior and is the infinite-dimensional generalization of the Dirichlet distribution. Therefore, for the sake of completeness, let us first recall the definition of the DP. A DP on the space 𝒢 is defined as a random process characterized by a concentration parameter α and a base distribution G0, denoted by GDP(α,G0), such that for any finite partition {A1,,Ap} of 𝒢, the random vector (G(A1),,G(Ap)) is Dirichlet distributed:

(G(A1),,G(Ap))Dir(αG0(A1),,αG0(Ap)).

We use the stick-breaking construction of the DP G (SBDP), due to Sethuraman (1994):

𝜽0kG0iidG0,k,
τkαiidBeta(τk1,α),k,
πk(𝝉)=τkl=1k1(1τl),k,
G=k=1πk(𝝉)δ𝜽0kDP(α,G0).
Pitman–Yor process.

As a generalized version of the Dirichlet process, in the Pitman–Yor process (PYP) (Pitman and Yor, 1997), the τk’s are independent (ind) but not identically distributed. Specifically,

𝜽0kG0iidG0,k,
τkα,σindBeta(1σ,α+kσ),k,
πk(𝝉)=τkl=1k1(1τl),k,
G=k=1πk(𝝉)δ𝜽0kPYP(α,σ,G0).

Here σ(0,1) is a discount parameter and α is a concentration parameter α>σ. The PYP is a two-parameter generalisation of the DP that allows one to control the tail behaviour when modelling data with either exponential or power-law tails (Ishwaran and James, 2001; Pitman and Yor, 1997). The PYP reduces to a DP when σ=0. More general stick-breaking representations are possible, e.g., Gibbs-type priors (De Blasi et al., 2015; Ghosal and Van der Vaart, 2017) or homogeneous normalised random measures with independent increments (Favaro et al., 2016). The PYP has a power-law behaviour for the number of clusters. This can make it more suitable for a number of applications. In other words, the number of clusters grows as 𝒪(Nσ) for PYP, while growing more slowly as 𝒪(logN) for DP.

Since the hyperparameters α and σ can have a significant effect on the growth of the number of clusters with data sample size, it is possible to specify priors for them. For the DP case obtained with σ=0, it is suggested in Blei and Jordan (2006) to use a gamma prior, αGam(s1,s2), where the hyperparameters s1 and s2 can be estimated or fixed. A natural question is whether one can also find a tractable prior for the discount parameter σ. Following the work of et al. (2020), we use the following prior that satisfies the constraints σ(0,1) and α>σ,

p(α,σs1,s2,a)=p(ασ;s1,s2)p(σa),

where p(ασ;s1,s2) is a shifted gamma distribution 𝒮𝒢(ασ;s1,s2) and p(σ,a) is a distribution depending on some parameter a which is not specified at the moment but which can typically be assumed to be a uniform distribution on the interval (0,1). Such a shifted gamma distribution is the distribution of a variable Uσ, where σ is considered fixed and U follows a gamma distribution Gam(s1,s2). The PDF of this shifted gamma distribution is obtained from the standard gamma distribution as 𝒮𝒢(ασ;s1,s2)=Gam(α+σs1,s2).

Hierarchical representation of BNP-MM.

BNP-MMs, including DP-MMs and PYP-MMs, see, e.g.,  et al. (2020); Durand et al. (2022), have the following hierarchical representation to generate a data point 𝐱n as a special case of BNP-GLLiMs:

1. BNP prior: aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
G=k=1πk(𝝉)δ𝜽kBNP(α,σ,G0),𝜽k=(𝝁k,𝐕k),
2. BNP-MM: for each n[N],
𝜽nGiidG,
If 𝜽n=𝜽k,setzn=k,
𝐱nzn,𝜽iid𝒩L(𝐱n𝜽zn),𝜽(𝜽k)k[K].

Appendix B Variational Bayesian expectation-maximisation principle

The clustering task consists mainly of estimating the unknown labels 𝒵=(𝐳n)n[N] from the observed data (𝒴,𝒳)=(𝐲n,𝐱n)n[N], whose joint distribution p(𝒴,𝒳,𝒵,𝚯;ϕ) is determined by a set of parameters denoted by 𝚯 and often by additional hyperparameters ϕ.

The expectation-maximisation (EM) algorithm (Dempster et al., 1977; McLachlan and Krishnan, 1997) is a generative technique for maximum likelihood estimation (MLE) in the presence of unobserved latent variables or missing data. An EM iteration consists of two steps usually referred to as the E-step in which the expectation of the so-called complete log-likelihood is computed and the M-step in which this expectation is maximized over 𝚯. An equivalent way to define EM is the following. As discussed in Neal and Hinton (1998), EM can be viewed as an alternating maximisation procedure of a function 0 defined, for any probability distribution q𝒵 over labels 𝒵, by

0(q𝒵,𝚯,ϕ) =𝒵q𝒵(𝒵)logp(𝒴,𝒳,𝒵𝚯,ϕ)𝔼q𝒵[logq𝒵(𝒵)]
=𝔼q𝒵[logp(𝒴,𝒳,𝒵𝚯,ϕ)q𝒵(𝒵)],

where 𝔼q𝒵[logq𝒵(𝒵)] is the entropy of q𝒵 and 𝔼q[] is the expectation with respect to q. The function 0 depends on the observations (𝒴,𝒳), which are fixed throughout and are therefore omitted from the notation.

When prior knowledge on the parameters is available, an alternative approach consists of replacing the MLE by a maximum a posteriori (MAP) estimation of 𝚯 using the prior knowledge encoded in a distribution p(𝚯). More precisely, the MLE of 𝚯 is then replaced by a point estimation 𝚯^=argmax𝚯𝚯p(𝚯𝒴,𝒳). In this paper, instead of considering only point estimation of 𝚯, we carry out a fully Bayesian approach. That is, we integrate out 𝚯 as follows

p(𝒵𝒴,𝒳)=𝚯p(𝒵𝒴,𝒳,𝚯)p(𝚯𝒴,𝒳)d𝚯.

This integration requires the computation of the density p(𝚯𝒴,𝒳), which is usually not available in closed-form. As an alternative to costly simulation-based methods (e.g., Markov chain Monte Carlo (MCMC)), an EM-like procedure using variational approximation can provide approximations of the marginal posterior distributions p(𝚯𝒴,𝒳) and p(𝒵𝒴,𝒳). This approach is referred to as VBEM for variational Bayesian EM, as introduced by Beal and Ghahramani (2003).

To deal with the BNP-GLLiM model, we need to use the VBEM with hyperparameter optimisation of Beal (2003, Figure 2.5 and Algorithm 5.3). Let q𝒵 and q𝚯 denote the distributions over 𝒵 and 𝚯, respectively, which will serve as approximations to the true posteriors. Specifically, in the Bayesian setting, the intractable posterior p(𝒵,𝚯𝒴,𝒳;ϕ) is approximated by the variational posterior q(𝒵,𝚯)=q𝒵(𝒵)q𝚯(𝚯).

Similar to standard EM, VBEM maximizes the following evidence lower bound (often abbreviated ELBO, and sometimes called the variational lower bound or negative variational free energy), defined for arbitrary q𝒵 and q𝚯 distributions by

(q𝒵,q𝚯,ϕ) =𝔼q𝒵q𝚯[logp(𝒴,𝒳,𝒵,𝚯;ϕ)q𝒵(𝒵)q𝚯(𝚯)] (20)
=logp(𝒴,𝒳ϕ)KL(q𝒵q𝚯p(𝒵,𝚯𝒴,𝒳,ϕ))logp(𝒴,𝒳ϕ),

alternatively over q𝒵,q𝚯 and ϕ. Here, KL stands for Kullback-Leibler divergence. It is worth noting that adding a prior on 𝚯 is formally equivalent to considering 𝚯 as missing variables, while the hyperparameters ϕ play the role of the parameters of interest in MLE.

The alternate maximisation over leads to the VBEM algorithm, which can be decomposed into three steps. It is easy to show, using the KL divergence properties, that the maximisation over q𝒵 and q𝚯 leads to the following E-steps, see, e.g., Chaari et al. (2013, Appendix A), Beal (2003, Theorem 2.1) and Bishop (2006, Section 10.1.1), which is essentially coordinate ascent in the function space of variational distributions. Furthermore, the following update rules for E-steps converge to a local maximum of (q𝒵,q𝚯,ϕ). At the rth iteration, using current values ϕ(r1) and q(r1)Θ, we get the following updating,

VB-E-𝒵q(r)𝒵(𝒵)exp𝔼q(r1)Θ[logp(𝒴,𝒳,𝒵,𝚯;ϕ(r1))], (21)
VB-E-𝚯q(r)𝚯(𝚯)exp𝔼q(r)𝒵[logp(𝒴,𝒳,𝒵,𝚯;ϕ(r1))], (22)
VB-M-ϕϕ(r)argmaxϕ𝔼q(r)𝒵q(r)𝚯[logp(𝒴,𝒳,𝒵,𝚯;ϕ)].

In practice, we can decide which parameters to treat as genuine parameters 𝚯 or as hyperparameters ϕ, depending on whether some prior knowledge is available for only a subset of the parameters, or whether the model has hyperparameters ϕ for which no prior information is available. Furthermore, for complex models, q𝚯 and q𝒵 may need to be further restricted to simpler forms, such as factorised forms, to ensure tractable VB-E steps. This is illustrated in the next Section F.1 for the BNP-GLLiM inference.

Appendix C Details of VBEM for BNP-GLLiM model

C.1 VB-E-𝝉 step from Section 3.1

To achieve results from Section 3.1, we make use of (5), (6), (40), and are only interested in the functional dependence of the right-hand side of (22) on the variable τk. Thus, any terms that do not depend on τk can be absorbed into the additive normalization constant, giving

qτk(τk) =exp{𝔼qα,σ[logp(τkα,σ)]+n=1N𝔼qznqτ{k}[logπzn(𝝉)]}+constant
exp{𝔼qα,σ[σ]logτk+[𝔼qα,σ[α]+k𝔼qα,σ[σ]1]log(1τk)
+n=1Nqzn(k)logτk+n=1Nj=k+1Kqzn(j)log(1τk)}
=Beta(τkγ^k,1,γ^k,2).

Here,

γ^k,1 =1𝔼qα,σ[σ]+n=1Nqzn(k)=1𝔼qσ[σ]+n=1Nqzn(k),
γ^k,2 =𝔼qα,σ[α]+k𝔼qα,σ[σ]+n=1Nj=k+1Kqzn(j)=𝔼qα,σ[α]+k𝔼qα,σ[σ]+j=k+1Kn=1Nqzn(j).

Furthermore, we used the fact that

logp(τkα,σ) =log[Beta(τk1σ,α+kσ)]=log[Γ(1σ+kσ)Γ(1σ)Γ(α+kσ)τk1σ1(1τk)α+kσ1],
logπzn(𝝉) =log[τznl=1zn1(1τl)]=logτzn+l=1zn1log(1τl),
qτ{k}(τl) =i=1,ikK1qτi(τl),dτ{k}=i=1,ikK1dτi.

Finally, we have for k[K], let Nk=n=1Nqzn(k) correspond to the weight of cluster k, then

qτk(τk) =Beta(τkγ^k,1,γ^k,2),
γ^k,1 =1𝔼qα,σ[σ]+Nk,γ^k,2=𝔼qα,σ[α]+k𝔼qα,σ[σ]+l=k+1KNl.

C.2 VB-E-(α,σ) step from Section 3.1

In the PY case, to achieve results from Section 3.1, we make use of (5), (4), (40), (41), and are only interested in the functional dependence of the right-hand side of (22) to the variables (α,σ). Thus, any terms that do not depend on (α,σ) can be included in the additive normalization constant. This results in qα,σ(α,σ) being proportional to

q~α,σ(α,σ)
=p(α,σs1,s2,a)exp{𝔼k=1K1qτk[logk=1K1p(τkα,σ)]}
=p(α,σs1,s2,a)k=1K1Γ(1σ+α+kσ)Γ(1σ)Γ(α+kσ)
×exp{k=1K1𝔼qτk[σlogτk+(α1+kσ)log(1τk)]}
=p(α,σs1,s2,a)1Γ(1σ)K1k=1K1[α+(k1)σ]Γ(α+(k1)σ)Γ(α+kσ)
×exp{k=1K1𝔼qτk[σ(logτkklog(1τk))+(α1)log(1τk)]}
=p(α,σs1,s2,a)1Γ(1σ)K1k=1K1[α+(k1)σ]l=0K2Γ(α+lσ)k=1K1Γ(α+kσ)
×exp{k=1K1𝔼qτk[σ(logτkklog(1τk))+(α1)log(1τk)]}
=p(α,σs1,s2,a)Γ(α)Γ(1σ)K1Γ(α+(K1)σ)k=1K1[α+(k1)σ]
×exp{σ[k=1K1𝔼qτk[logτk]k=1K1k𝔼qτk[log(1τk)]]+(α1)k=1K1𝔼qτk[log(1τk)]},

where we used the fact that Γ(x+1)=xΓ(x). Except in the DP-GLLiM case, i.e.,  σ=0, the normalizing constant, (q~α,σ(α,σ)d(α,σ))1, for q~α,σ is not tractable. However, to perform VBEM in Section F.1, we do not need the full qα,σ distribution, but only the means 𝔼qα,σ[α] and 𝔼qα,σ[σ]. One solution, therefore, is to use importance sampling or MCMC to compute these expectations by means of Monte Carlo sums. Via the prior on (α,σ) given in (4), it holds that

q~α,σ(α,σ)
=1Γ(s1)s2s1(α+σ)s11exp{s2(α+σ)}Γ(α)p(σa)Γ(1σ)K1Γ(α+(K1)σ)k=1K1[α+(k1)σ]
×exp{σ[k=1K1𝔼qτk[logτk]k=1K1k𝔼qτk[log(1τk)]]+(α1)k=1K1𝔼qτk[log(1τk)]}
=1Γ(s1)s2s1(α+σ)s11exp{[s2k=1K1𝔼qτk[log(1τk)]](α+σ)}
×eσξexp{k=1K1𝔼qτk[log(1τk)]}Γ(α)p(σa)Γ(1σ)K1Γ(α+(K1)σ)k=1K1[α+(k1)σ]
=1Γ(s1)(α+σ)s11exp{[s2k=1K1𝔼qτk[log(1τk)]](α+σ)}[s2k=1K1𝔼qτk[log(1τk)]]s1s1
×eσξexp{k=1K1𝔼qτk[log(1τk)]}Γ(α)p(σa)Γ(1σ)K1Γ(α+(K1)σ)k=1K1[α+(k1)σ]s2s1
Gam(α+σs^1,s^2)eσξΓ(α)p(σa)Γ(1σ)K1Γ(α+(K1)σ)k=1K1α+(k1)σα+σ. (23)

Here, given that ψ() is the digamma function defined by ψ(z)=ddzlogΓ(z)=Γ(z)Γ(z), we have

ξ =k=1K1𝔼qτk[logτk]k=1K1(k1)𝔼qτk[log(1τk)],
𝔼qτk[logτk] =ψ(γ^k,1)ψ(γ^k,1+γ^k,2),
𝔼qτk[log(1τk)] =ψ(γ^k,2)ψ(γ^k,1+γ^k,2),
s^1 =s1+K1,s^2=s2k=1K1[ψ(γ^k,2)ψ(γ^k,1+γ^k,2)].

We propose to use the following important distribution ν(α,σ)=Gam(α+σs^1,s^2)p(σa) where p(σa) is the uniform distribution on [0,1], denoted as 𝒰[0,1](σ). Then we obtain an expression for the importance weights,

W(α,σ) =q~α,σ(α,σ)ν(α,σ)=eσξΓ(α)Γ(1σ)K1Γ(α+(K1)σ)k=1K1α+(k1)σα+σ.

The importance sampling scheme then consists of the following steps

  • For m[M], first simulate independently σm from 𝒰[0,1](σ) and then simulate conditionally αm with the σm-shifted gamma 𝒢(σm,s^1,s^2). This later simulation is easily obtained by simulating a standard γ(αms^1,s^2) and then subtracting σm from the result.

  • Compute the importance weights wm=W(αm,αm).

  • Approximate the means

    𝔼qα,σ[α] m=1Mwmi=1Mwiαm,𝔼qα,σ[σ]m=1Mwmi=1Mwiσm.

Note that this complication is due to the PY. In the DP-GLLiM case, by substituting σ=0 in (23), the E step is much simpler, as it reduces to computing the approximate posterior expectation of α, namely,

𝔼qα,0[α] =s^1s^2,qα,0=Gam(αs^1,s^2).

C.3 VB-E-𝐙 step from Section 3.1

In some situations, it is useful to use a 1-of-K binary vector 𝐳n to represent the latent variable zn for each observation (𝐲n,𝐱n). To be more precise, we introduce a K-dimensional binary random variable 𝐳n=(znk)k[K],K, with a 1-of-K representation in which a particular element znk is equal to 1, i.e.,  zn=k, and all other elements are equal to 0. The values of znk thus satisfy znk{0,1} and k[K]znk=1,n. If there is no confusion, we also denote 𝒵 as the latent matrix 𝒵=(znk)n[N],k[K]. It is worth to mentioning that when using a 1-of-K representation of 𝐳n, we can also write down marginal the conditional distributions of 𝒳 and 𝒴𝒳, corresponding to (7), in the form

p(𝒳𝒵,𝒄,𝚪) =n=1Nk=1K𝒩L(𝐱n𝒄k,𝚪k)znk, (24)
p(𝐱n𝒄,𝚪) =k=1Kpzn(k)𝒩L(𝐱n𝒄k,𝚪k),pzn(k)p(zn=k)=πk(𝝉), (25)
p(𝒴𝒳,𝒵;𝐀,𝐛,𝚺) =n=1Nk=1K𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k)znk, (26)
p(𝐲n𝐱n;𝐀,𝐛,𝚺) =k=1Kpzn(k)𝒩L(𝐱n𝒄k,𝚪k)l=1Kpzn(l)𝒩L(𝐱n𝒄l,𝚪l)𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k).

The observations 𝒳 and 𝒴 are therefore i.i.d. and generated from the same GMM (25) and infinite GLLiM (26), respectively. Similarly, (6) can be written down in the form

p(𝒵𝝉)=n=1Nk=1Kπk(𝝉)znk.

By using the decomposition (8), the representation (24), (26) and absorbing any terms that are independent on 𝒵 into the additive normalization constant, we obtain

logq𝒵(𝒵)logq(r)𝒵(𝒵)
=𝔼q𝚯[logp(𝒴𝒳,𝒵;𝐀^,𝐛^,𝚺^)p(𝒳𝒵;𝐜^,𝚪^)p(𝒵𝝉)]+constant
𝔼q𝚯[logp(𝒴𝒳,𝒵;𝐀^,𝐛^,𝚺^)]+𝔼q𝚯[logp(𝒳𝒵;𝐜^,𝚪^)]+𝔼q𝚯[logp(𝒵𝝉)]
=n=1Nk=1Kznk𝔼q𝚯[log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)]+n=1Nk=1Kznk𝔼q𝚯[log𝒩L(𝐱n𝐜^k,𝚪^k)]
+𝔼q𝚯[n=1Nlog(πzn(𝝉))]
=n=1Nk=1Kznklog𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)+n=1Nk=1Kznklog𝒩L(𝐱n𝐜^k,𝚪^k)
+n=1Nk=1Kznk𝔼q𝝉[log(πk(𝝉))]=n=1Nk=1Kznklogρnk.

Here, we used the fact that

logρnk =log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)+log𝒩L(𝐱n𝐜^k,𝚪^k)+𝔼q𝝉[log(πk(𝝉))]
=D2log(2π)12log|𝚺^k|12(𝐲n𝐀^k𝐱n𝐛^k)𝚺^1k(𝐲n𝐀^k𝐱n𝐛^k)
L2log(2π)12log|𝚪^k|12(𝐱n𝐜^k)𝚪^1k(𝐱n𝐜^k)
+ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1ψ(γ^l,2)ψ(γ^l,1+γ^l,2).

By taking exponential of both sides and taking into account the normalized constant, it holds that

q𝒵(𝒵) =1l=1Kρnln=1Nk=1Kρnkznk=n=1Nk=1Krnkznk,rnk=ρnkl=1Kρnl.

Note also that zn=k if and only if the latent matrix 𝒵 reduces to a sparse matrix 𝒵nk which has only one position different from 0, namely znk=1. This leads to the following simplified notation:

logqzn(k) logqzn(zn=k)logq𝒵(𝒵nk)=rnk.

C.3.1 Updating 𝚺k

By using matrix derivatives, the derivative of the log likelihood with respect to 𝚺1k is given by

𝚺1kf1(𝐀^k,𝐛k,𝚺1k)
=12n=1Nqzn(k)𝚺1k[log|𝚺1k|+Tr[𝚺1k(𝐲n𝐀k𝐱n𝐛k)(𝐲n𝐀k𝐱n𝐛k)]]
=12n=1Nqzn(k)[𝚺k+(𝐲n𝐀k𝐱n𝐛k)(𝐲n𝐀k𝐱n𝐛k)]
=Nk2𝚺k12n=1Nqzn(k)(𝐲n𝐀k𝐱n𝐛k)(𝐲n𝐀k𝐱n𝐛k).

Finally, setting to zero yields

𝚺^k =1Nkn=1Nqzn(k)(𝐲n𝐀^k𝐱n𝐛^k)(𝐲n𝐀^k𝐱n𝐛^k).

C.4 VB-M-(𝐜,𝚪) step from Section 3.2

This step divides into K sub-steps that involve the following optimisations

(𝐜^k,𝚪^k)(𝐜^k(r),𝚪^k(r)) =argmax(𝐜k,𝚪k)𝔼q𝒵(r)[logp(𝒳𝒵;𝐜k,𝚪k)].

By definition, we have

𝔼q𝒵(r)[log(p(𝒳𝒵;𝐜k,𝚪k))]
=𝔼q𝒵(r)[logn=1N𝒩L(𝐱n𝒄k,𝚪k)znk]
=n=1N𝔼q𝒵(r)[znklog𝒩L(𝐱n𝒄k,𝚪k)]
=n=1N𝔼q𝒵(r)[znk]log𝒩L(𝐱n𝒄k,𝚪k)
=n=1Nqzn(k)log𝒩L(𝐱n𝒄k,𝚪k)
=n=1Nqzn(k)[L2log(2π)12log|𝚪k|12(𝐱n𝐜k)𝚪1k(𝐱n𝐜k)]
f2(𝐜k,𝚪k).

We aim to solve the following optimisation

(𝐜^k,𝚪^k) =argmax(𝐜k,𝚪k)f2(𝐜k,𝚪k).

Similarly with Sections C.5.1 and C.3.1, we obtain the following update:

𝐜^k =1Nkn=1Nqzn(k)𝐱n,
𝚪^k =1Nkn=1Nqzn(k)(𝐱n𝐜^k)(𝐱n𝐜^k).

C.5 VB-M-(𝐀,𝐛,𝚺) step from Section 3.2

By definition, we have

𝔼q𝒵(r)[logp(𝒴𝒳,𝒵;𝐀k,𝐛k,𝚺k)]
=𝔼q𝒵(r)[logn=1N𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k)znk]
=n=1N𝔼q𝒵(r)[znklog𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k)]
=n=1N𝔼q𝒵(r)[znk]log𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k)
=n=1Nqzn(k)log𝒩D(𝐲n𝐀k𝐱n+𝐛k,𝚺k)
=n=1Nqzn(k)[D2log(2π)12log|𝚺k|12(𝐲n𝐀k𝐱n𝐛k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)]
f1(𝐀k,𝐛k,𝚺k).

We aim to solve the following optimisation

(𝐀^k,𝐛^k,𝚺^k) =argmax(𝐀k,𝐛k,𝚺k)f1(𝐀k,𝐛k,𝚺k).

C.5.1 Updating 𝐛k

The derivative of the log likelihood with respect to 𝐛k is given by

𝐛kf1(𝐀k,𝐛k,𝚺k) =12n=1Nqzn(k)𝐛k[(𝐲n𝐀k𝐱n𝐛k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)]
=n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)𝐛k(𝐲n𝐀k𝐱n𝐛k)
=n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k).

Setting this derivative to zero, we obtain the solution for VB-M-𝐛 step given by

n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)=0
n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n)n=1Nqzn(k)𝚺1k𝐛k=0
n=1Nqzn(k)(𝐲n𝐀k𝐱n)n=1Nqzn(k)𝐛k=0(left multiplying by 𝚺k)
𝐛k=1n=1Nqzn(k)n=1Nqzn(k)(𝐲n𝐀k𝐱n)1Nkn=1Nqzn(k)(𝐲n𝐀k𝐱n). (27)

C.5.2 Updating 𝐀k

The derivative of the log likelihood with respect to 𝐀k is given by

𝐀kf1(𝐀k,𝐛k,𝚺k) =12n=1Nqzn(k)𝐀k[(𝐲n𝐀k𝐱n𝐛k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)]
=n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)𝐀k(𝐲n𝐀k𝐱n𝐛k)
=n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)𝐱n.

Then, we set this derivative w.r.t. 𝐀k equal to zero, giving

n=1Nqzn(k)𝚺1k(𝐲n𝐀k𝐱n𝐛k)𝐱n=0
n=1Nqzn(k)(𝐲n𝐀k𝐱n𝐛k)𝐱n=0(left multiplying by 𝚺k)
n=1Nqzn(k)𝐲n𝐱nn=1Nqzn(k)𝐀k𝐱n𝐱nn=1Nqzn(k)𝐛k𝐱n=0
n=1Nqzn(k)𝐲n𝐱n𝐀kn=1Nqzn(k)𝐱n𝐱n𝐛kn=1Nqzn(k)𝐱n=0
𝐀kn=1Nqzn(k)𝐱n𝐱n
=n=1Nqzn(k)𝐲n𝐱n1Nkn=1Nqzn(k)(𝐲n𝐀k𝐱n)n=1Nqzn(k)𝐱n(using (27) for 𝐛k)
𝐀kn=1Nqzn(k)𝐱n(𝐱n1Nkn=1Nqzn(k)𝐱n)=n=1Nqzn(k)(𝐲n1Nkn=1Nqzn(k)𝐲n)𝐱n
Nk𝐀k𝑿k𝑿k=Nk𝒀k𝑿k𝐀k=𝒀k𝑿k(𝑿k𝑿k)1.

Here, the last equality is obtained by firstly define the following quantities,

𝐱¯k =1Nkn=1Nqzn(k)𝐱n,
𝐲¯k =1Nkn=1Nqzn(k)𝐲n,
𝑿k =1Nk(qz1(k)(𝐱1𝐱¯k),,q𝐳n(k)(𝐱N𝐱¯k))L×N,
𝒀k =1Nk(qz1(k)(𝐲1𝐲¯k),,q𝐳n(k)(𝐲N𝐲¯k))D×N.

Then, we used the fact that

n=1Nqzn(k)𝐱n(𝐱n1Nkn=1Nqzn(k)𝐱n)
=n=1Nqzn(k)𝐱n(𝐱n𝐱¯k)
=n=1Nqzn(k)(𝐱n𝐱¯k)(𝐱n𝐱¯k)+n=1Nqzn(k)𝐱¯k(𝐱n𝐱¯k)
=n=1Nqzn(k)(𝐱n𝐱¯k)(𝐱n𝐱¯k)
=Nk𝑿k𝑿k,

and

n=1Nqzn(k)(𝐲n1Nkn=1Nqzn(k)𝐲n)𝐱n
=n=1Nqzn(k)(𝐲n𝐲¯k)𝐱n
=n=1Nqzn(k)(𝐲n𝐲¯k)(𝐱n𝐱¯k)+n=1Nqzn(k)(𝐲n𝐲¯k)𝐱¯k
=n=1Nqzn(k)(𝐲n𝐲¯k)(𝐱n𝐱¯k)
=Nk𝒀k𝑿k.

Here, we also use the equalities

n=1Nqzn(k)𝐱¯k(𝐱n𝐱¯k) =n=1Nqzn(k)𝐱¯k𝐱n𝐱¯kNk𝐱¯k=n=1Nqzn(k)𝐱¯k𝐱nn=1Nqzn(k)𝐱¯k𝐱n=0,
n=1Nqzn(k)(𝐲n𝐲¯k)𝐱¯k =n=1Nqzn(k)𝐲n𝐱¯kNk𝐲¯k𝐱¯k=n=1Nqzn(k)𝐲n𝐱¯kn=1Nqzn(k)𝐲n𝐱¯k=0,

and for each i,j[L], it holds that

[n=1Nqzn(k)(𝐱n𝐱¯k)(𝐱n𝐱¯k)]ij=n=1Nqzn(k)[(𝐱n𝐱¯k)(𝐱n𝐱¯k)]ij
=n=1Nqzn(k)[(𝐱n𝐱¯k)]i1[(𝐱n𝐱¯k)]1j
n=1Nqzn(k)[(𝐱n𝐱¯k)]i[(𝐱n𝐱¯k)]j
=[(qz1(k)(𝐱1𝐱¯k),,q𝐳n(k)(𝐱N𝐱¯k))]i[(qz1(k)(𝐱1𝐱¯k),,q𝐳n(k)(𝐱N𝐱¯k))]j
=Nk[𝑿k𝑿k]ij.

Appendix D Details on the ELBO

In this section, we provide the closed-form expressions for the ELBO stated in Proposition 3.1. Let us recal that when σ=0, the ELBO in the BNP-GLLiM is derived as follows:

(q𝒵,q𝚯,ϕ^) =𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]+𝔼[logp(𝒳𝒵,𝚯;ϕ^)]+𝔼[logp(𝒵𝚯;ϕ^)]
+𝔼[logp(𝚯;ϕ^)]𝔼[logq(𝒵)]𝔼[logq(𝚯)],

The terms of the right-hand side of the above equation have the following closed-form expressions:

𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)] =n=1Nk=1Kqzn(k)log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k) (28)
𝔼[logp(𝒳𝒵,𝚯;ϕ^)] =n=1Nk=1Kqzn(k)log𝒩L(𝐱n𝒄^k,𝚪^k), (29)
𝔼[logp(𝒵𝚯;ϕ^)] =k=1KNk[ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1[ψ(γ^l,2)ψ(γ^l,1+γ^l,2)]], (30)
𝔼[logp(𝚯;ϕ^)] =k=1K1𝔼[logp(τkα)]+𝔼[logp(αs^1,s^2)], (31)
𝔼[logp(τkα)] =s^1s^2s^2[ψ(γ^k,2)ψ(γ^k,1+γ^k,2)]+ψ(s^1)log(s^2),
𝔼[logp(αs^1,s^2)] =logΓ(s^1)+(s^11)ψ(s^1)+log(s^2)s^1,
𝔼[logq(𝒵)] =n=1Nk=1Kqzn(k)logqzn(k), (32)
𝔼[logq(𝚯)] =𝔼[logqα,0(α)]+k=1K1𝔼[logqτk(τk)], (33)
𝔼[logqα,0(α)] =logΓ(s^1)+(s^11)ψ(s^1)+log(s^2)s^1,
𝔼[logqτk(τk)] =l=12(γ^k,l1){ψ(γ^k,l)ψ(γ^k,1+γ^k,2)}+logΓ(γ^k,1+γ^k,2)Γ(γ^k,1)Γ(γ^k,2).

Appendix E Technical proofs

E.1 Proof of Lemma 5.2

We first want to prove (5.2). Using the partition of a joint Gaussian with 𝐱b𝐱,𝝁b=𝝁𝐱k,𝚺bb𝐕𝐱𝐱k,𝐱a𝐲,𝝁a𝝁𝐲k,𝚺aa𝐕𝐲𝐲k, we obtain

p(𝐱a𝐱b) =𝒩(𝐱a𝝁ab,𝚪1aa),𝝁ab=𝝁a𝚪1aa𝚪ab(𝒙b𝝁b)=𝝁a𝚺ab𝚺1bb(𝒙b𝝁b),
p(𝒙b) =𝒩(𝒙b𝝁b,𝚺bb). (34)

Recall that

p(𝐲𝐱,Z=k;𝝍) =𝒩D(𝐲𝐀k𝐱+𝐛k,𝚺k),
p(𝐱Z=k;𝝍) =𝒩L(𝐱𝐜k,𝚪k),p(Z=k;𝝍)=πk. (35)

By identifying the parameters of (E.1) and (E.1), it holds that

πk =ρk
𝐜k =𝝁𝐱k,
𝚪k =𝐕𝐱𝐱k,
𝐀k =𝚪1aa𝚪ab=𝐕𝐱𝐲k(𝐕𝐱𝐱k)1,
𝐛k =𝝁a+𝚪1aa𝚪ab𝝁b=𝝁𝐲k𝐕𝐱𝐲k(𝐕𝐱𝐱k)1𝝁𝐱k,
𝚺k =𝚪1aa=𝚺aa𝚺ab𝚺bb1𝚺ba=𝐕𝐲𝐲k𝐕𝐱𝐲k(𝐕𝐱𝐱k)1𝐕𝐱𝐲k.

The following decomposition of the joint probability distribution will be used:

p(𝐰𝝍) =k=1Kp(𝐲𝐱,Z=k;𝝍)p(𝐗=𝐱Z=k;𝝍)p(Z=k;𝝍)
=k=1Kπk𝒩D(𝐲𝐀k𝐱+𝐛k,𝚺k)𝒩L(𝐱𝐜k,𝚪k)
k=1Kρk𝒩L+D(𝐰𝝁k,𝐕k).

By using result for the joint Gaussian, see e.g., (39), we obtain the desired result (19).

Finally, Lemma 5.2 is proved via using the following two statements (Deleforge et al., 2015, Lemmas 1 and 2):

  1. (i)

    For any ρk,𝝁kL+D, and 𝐕k𝒮L+D+, there is a set of parameters 𝐜kL,λk𝒮L+,πk,𝐀kD×L,𝐛kD,𝚺k𝒮D+ such that (5.2) holds.

  2. (ii)

    Reciprocally, for any 𝐜kL,𝚲k𝒮L+,πk,𝐀kD×L,𝐛kD,𝚺k𝒮D+, there is a set of parameters ρk,𝝁kL+D and 𝐕k𝒮L+D+ such that (19) holds.

E.2 Proof of Proposition 3.1

Using the sum and product rules for both discrete and continuous variables, the ELBO in BNP-GLLiM (20) is given by

(q𝒵,q𝚯,ϕ^) =𝔼q𝒵q𝚯[logp(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q(𝚯)]𝔼[logp(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q(𝚯)]
=𝒵q(𝒵)q(𝚯)log[p(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q(𝚯)]d𝒵d𝚯
=𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]+𝔼[logp(𝒳𝒵,𝚯;ϕ^)]+𝔼[logp(𝒵𝚯;ϕ^)]
+𝔼[logp(𝚯;ϕ^)]𝔼[logq(𝒵)]𝔼[logq(𝚯)]. (36)

Next, we evaluate the various terms in the ELBO (36).

Proof of (28)

Via the mean field approximation and the truncation, we have the following computations:

𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)] =𝔼[logn=1Np(𝐲n𝐱n,zn,𝚯;ϕ^)]
=𝔼[logn=1Nk=1K𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)znk]
=n=1Nk=1K𝔼[znklog𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)]
=n=1Nk=1K𝔼q𝒵[znk]log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k)
=n=1Nk=1Kqzn(k)log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k),

where

log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k) =D2log(2π)12log|𝚺^k|
(𝐲n𝐀^k𝐱n𝐛^k)𝚺^k1(𝐲n𝐀^k𝐱n𝐛^k).

Proof of (29)

Similarly to the previous proof, we obtain

𝔼[logp(𝒳𝒵,𝚯;ϕ^)] =𝔼[logn=1Np(𝐱nzn,𝚯;ϕ^)]
=𝔼[logn=1Nk=1K𝒩L(𝐱n𝒄^k,𝚪^k)znk]
=n=1Nk=1K𝔼[znklog𝒩L(𝐱n𝒄^k,𝚪^k)]
=n=1Nk=1K𝔼q𝒵[znk]log𝒩L(𝐱n𝒄^k,𝚪^k)
=n=1Nk=1Kqzn(k)log𝒩L(𝐱n𝒄^k,𝚪^k),

where

log𝒩L(𝐱n𝒄^k,𝚪^k)=L2log(2π)12log|𝚪^k|12(𝐱n𝒄^k)𝚪^k1(𝐱n𝒄^k).

Proof of (30)

Via calculation, it follows the expressions of the following quantities,

𝔼qτk[log(τk)] =ψ(γ^k,1)ψ(γ^k,1+γ^k,2),
𝔼qτk[log(1τk)] =ψ(γ^k,2)ψ(γ^k,1+γ^k,2). (37)

Via (37), it holds that

𝔼[logp(𝒵𝚯;ϕ^)] =𝔼[logn=1Nk=1K[πk(𝝉)]znk]
=n=1Nk=1K𝔼q𝒵[znk]𝔼q𝚯[log[τkl=1k1(1τl)]]
=k=1Kn=1Nqznk[𝔼qτk[logτk]+l=1k1𝔼qτl[log(1τl)]]
=k=1KNk[𝔼qτk[logτk]+l=1k1𝔼qτl[log(1τl)]]
=k=1KNk[ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1[ψ(γ^l,2)ψ(γ^l,1+γ^l,2)]].

Proof of (31)

Given a chosen truncated value K, it holds that

𝔼q𝚯[logp(𝚯;ϕ^)] =k=1K1𝔼q𝚯[logp(τkα,σ)]+𝔼q𝚯[logp(α,σs^1,s^2,a^)].

Here, we have

𝔼q𝚯[logp(τkα,σ)] =𝔼q𝚯[logBeta(τk1σ,α+kσ)]
=𝔼q𝚯[logτkσ(1τk)α+kσ1+logC(α,σ)],
=𝔼qα,σ[σ]𝔼qτk[logτk]+𝔼qα,σ[α+kσ1]𝔼qτk[log(1τk)]
+𝔼qα,σ[logC(α,σ)],

where we have defined

C(α,σ)=Γ(1σ+α+kσ)Γ(1σ)Γ(α+kσ).

Next, for the sake of simplicity, for σ, we use a uniform prior 𝒰[0,1](σ) so that parameter a does not have to be taken into account. Then it holds that

𝔼q𝚯[logp(α,σs^1,s^2)] =𝔼qα,σ[logGam(α+σs^1,s^2)]+𝔼qα,σ[log𝒰[0,1](σ)]
=log[1Γ(s^1)s^2s^1]+(s^11)𝔼qα,σ[log(α+σ)]s^2𝔼qα,σ[α+σ]
=log[1Γ(s^1)s^2s^1]+(s^11)𝔼qα,σ[log(α+σ)]s^2𝔼qα,σ[α+σ].

When σ0, the normalizing constant for qα,σ(α,σ) is not tractable. Nevertheless, to compute the ELBO, we do not need the full qα,σ distribution but only the means 𝔼qα,σ[σ], 𝔼qα,σ[α+kσ1], 𝔼qα,σ[logC(α,σ)], 𝔼qα,σ[log(α+σ)] and 𝔼qα,σ[α+σ]. One solution is therefore to use importance sampling or MCMC to compute these expectations via Monte Carlo sums.

When σ=0, using integration by parts, it holds that Γ(α+1)=αΓ(α) and hence C(α,σ)C(α)=α. Furthermore, the posterior qα,σqα is again a gamma distribution Gam(αs^1,s^2) with 𝔼qα,σ[α]𝔼qα[α]=s^1s^2 and 𝔼qα,σ[logα]𝔼qα[logα]=ψ(s^1)log(s^2). Therefore, we have the following tractable formulas:

𝔼q𝚯[logp(τkα,σ)] 𝔼q𝚯[logp(τkα)]
=[𝔼qα,0[α]1]𝔼qτk[log(1τk)]+𝔼qα,0[logα],
=s^1s^2s^2[ψ(γ^k,2)ψ(γ^k,1+γ^k,2)]+ψ(s^1)log(s^2),
𝔼q𝚯[logp(α,σs^1,s^2)] 𝔼q𝚯[logp(αs^1,s^2)]
=log[1Γ(s^1)s^2s^1]+(s^11)[ψ(s^1)log(s^2)]s^1.

Proof of (32)

Due to the mean-field approximation (9) and truncation, this step is analytically computed as follows:

𝔼q𝒵[logq(𝒵)] =𝔼q𝒵[logn=1Nqzn(zn)]=𝔼q𝒵[logn=1Nk=1Kqzn(k)znk]
=n=1Nk=1Klogqzn(k)𝔼q𝒵[znk]=n=1Nk=1Kqzn(k)logqzn(k).

Proof of (33)

We have

𝔼[logq(𝚯)] =𝔼[logqα,σ(α,σ)]+k=1K1𝔼[logqτk(τk)].

Note that these terms involving expectations of the logs of the q distributions simply represent the negative entropies of those distributions.

Since qα,σ(α,σ) is not tractable, when σ0, we cannot calculate analytically 𝔼[logqα,σ(α,σ)]. Furthermore, it is also difficult to approximate it using MCMC or importance sampling.

When σ=0, the posterior qα,σqα is again a gamma distribution Gam(αs^1,s^2) with

𝔼[logqα,0(α)] 𝔼[logGam(αs^1,s^2)]
=H[Gam(αs^1,s^2)]
=logΓ(s^1)+(s^11)ψ(s^1)+log(s^2)s^1.

Since we had qτk(τk)=Beta(τkγ^k,1,γ^k,2), its differential entropy is given by

𝔼[logqτk(τk)] =H[Beta(τkγ^k,1,γ^k,2)]
=l=12(γ^k,l1){ψ(γ^k,l)ψ(γ^k,1+γ^k,2)}+logΓ(γ^k,1+γ^k,2)Γ(γ^k,1)Γ(γ^k,2).

E.3 Proof of Theorem 4.1

Recall that 𝚯=(𝝉,α,σ). Then,

p(𝐲^,𝐱^,𝒳,𝒴) =𝐳^p(𝐲^𝐱^,𝐳^,𝚯,𝒳,𝒴)p(𝐱^𝐳^,𝚯,𝒳,𝒴)p(𝐳^𝚯,𝒳,𝒴)p(𝚯𝒳,𝒴)d𝚯
=𝐳^p(𝐲^𝐱^,𝐳^;𝑨,𝒃,𝚺)p(𝐱^𝐳^,𝒄,𝚪)p(𝐳^𝝉;β)p(𝚯𝒳,𝒴)d𝚯D1. (38)

Note that in (E.3), p(𝚯𝒳,𝒴) is in fact the (unknown) true posterior distribution of the parameters given a sample (𝒳,𝒴). Because the integrations w.r.t. true posterior distribution are intractable, we approximate the predictive conditional density by replacing the true posterior distribution p(𝚯𝒳,𝒴) with its truncated variational posterior of parameters 𝚯 given by

qΘ(𝚯𝒳,𝒴)=qα,σ(α,σ𝒳,𝒴)k=1K1qτk(τk𝒳,𝒴).

Recall that the infinite state space for each zj is dealt with by choosing a truncation of the state space to a maximum label K (Blei and Jordan, 2006). In practice, this consists of assuming that the variational distributions qzn for n[N], satisfy qzn(k)=0 for k>K and that the variational distribution on τ also factorizes as q𝝉(𝝉)=k=1K1qτk(τk) with the additional condition that τK=1. In particular, here we choose K=K^, where K^ is estimated from some suitable procedures.

For simplicity, we consider the case when β=0,σ=0. Then we have

D1 𝐳^p(𝐲^𝐱^,𝐳^;𝑨^,𝒃^,𝚺^)p(𝐱^𝐳^,𝒄^,𝚪^)p(𝐳^𝝉)qΘ(𝚯𝒳,𝒴)d𝚯
=k=1𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜^k,𝚪^k)πk(𝝉)qΘ(𝚯𝒳,𝒴)d𝚯
k=1K𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)πk(𝝉)qΘ(𝚯𝒳,𝒴)d𝚯
=k=1K𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)πk(𝝉)q𝝉(𝝉𝒳,𝒴)d𝝉qα,0(α𝒳,𝒴)dα=1
=k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)
k=1K𝔼q𝝉[πk(𝝉)]𝒩L+D(𝐰^𝔼[𝐰],cov[𝐰])
=k=1K𝔼qτk[τk]l=1k1𝔼qτl[1τl]𝒩L+D(𝐰^𝔼[𝐰],cov[𝐰]).

Here, by defining 𝐰^[𝐱^;𝐲^], we used the fact that

𝔼[𝐰] =(𝐜^k𝐀^k𝐜^k+𝐛^k),cov[𝐰]=(𝚪^k𝚪^k𝐀^k𝐀^k𝚪^k𝚺^k+𝐀^k𝚪^k𝐀^k).

Indeed, we made use of the following result for the joint Gaussian, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for 𝐱 and a conditional Gaussian distribution for 𝐲 given 𝐱 in the form

p(𝐱) =𝒩(𝐱𝝁,𝚪1),
p(𝐲𝐱) =𝒩(𝐲𝐀𝐱+𝐛,𝐋1),

then the joint distribution of 𝐰[𝐱;𝐲] is given by

p(𝐰) =𝒩(𝐰𝔼[𝐰],cov[𝐰]), where
cov[𝐰] =(𝚪1𝚪1𝐀𝐀𝚪1𝐋1+𝐀𝚪1𝐀),𝔼[𝐰]=(𝝁𝐀𝝁+𝐛). (39)

In our situation, the desired result is obtained via using 𝐲𝐲^,𝐀𝐀^k,𝐛𝐛^k,𝐋1=𝚺^k, 𝐱𝐱^,𝝁𝐜^k,𝚪1𝚪^k.

Furthermore, we also used the fact that

𝔼q𝝉[πk(𝝉)] =τkqτk(τk𝒳,𝒴)dτkl=1k1(1τl)j=1,jkK1qτj(τj𝒳,𝒴)j=1,jkKdτj
=𝔼qτk[τk]l=1k1(1τl)j=1k1qτj(τj𝒳,𝒴)j=k+1K1qτj(τj𝒳,𝒴)j=k+1K1dτj=1j=1k1dτj
=𝔼qτk[τk]l=1k1(1τl)qτl(τl𝒳,𝒴)dτl
=𝔼qτk[τk]l=1k1𝔼qτl[1τl].

Next, we aim to prove that

k=1K𝔼q𝝉[πk(𝝉)]=1.

Indeed, recall that we have defined

τkα,σindBeta(τk1σ,α+kσ),k,
πk(𝝉)=τkl=1k1(1τl),k,
p(𝒵𝝉)n=1Nπzn(𝝉),

and to deal with the infinite state space for each zj, we considered a truncation of the state space to a maximum label KKmax,Kmax (Blei and Jordan, 2006). In practice, this consists of assuming that the variational distributions qzn for n[N], satisfy qzn(k)=0 for k>K and that the variational distribution on τ also factorizes as q𝝉(𝝉)=k=1K1qτk(τk) with the additional condition that τK=1. Based on the proof from Ghosal and Van der Vaart (2017, Lemma 3.4), it holds that a necessary and sufficient condition to guarantee that these πk’s sum to 1 almost surely, i.e.,

k=1πk(𝝉)=k=1τkl=1k1(1τl)=1,

is that the expectation 𝔼[l=1k1(1τl)] tends to 0 as k tends to . In particular, if τ1,τ2, are i.i.d., e.g., when σ=0, it suffices that p(τ1>0)>0. Then

1 =𝔼q𝝉[k=1πk(𝝉)]=k=1𝔼q𝝉[πk(𝝉)]=k=1K𝔼q𝝉[πk(𝝉)].

E.4 Proof of Theorem 4.2

From the product rule of probability, we see that this conditional distribution can be evaluated from the joint and marginal distributions. Furthermore, by integrating out 𝐳^ and 𝚯, the predictive conditional density is then given by

p(𝐲^𝐱^,𝒳,𝒴) =p(𝐲^,𝐱^𝒳,𝒴)p(𝐱^𝒳,𝒴)=𝐳^p(𝐲^,𝐱^,𝐳^,𝚯𝒳,𝒴)d𝐳^d𝚯𝐳^p(𝐱^,𝐳^,𝚯𝒳,𝒴)d𝐳^d𝚯D1D2.

Next, with a similar step as in the proof of Theorem 4.1, we also obtain

D2k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k).

Therefore, we obtain

p(𝐲^𝐱^,𝒳,𝒴)
k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)
=k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)
k=1Kgk(𝐱^𝚯^,ϕ^,𝒳,𝒴)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k),

which is a mixture of Gaussian experts since we have

gLk(𝐱^𝚯^,ϕ^,𝒳,𝒴)=𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k),k[K],

belongs to a K1 dimensional probability simplex.

E.5 Proof of Theorem 4.3

To deal with high-dimensional regression data, namely high-to-low regression, given the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴), we want to compute the following forward conditional density

p(𝐱^𝐲^,𝒳,𝒴) =p(𝐱^,𝐲^𝒳,𝒴)p(𝐲^𝒳,𝒴)=p(𝐱^,𝐲^𝒳,𝒴)𝐱^p(𝐱^,𝐲^𝒳,𝒴)d𝐱^=D1𝐱^D1(𝐱^)d𝐱^D1D3.

Then, we have to compute or approximate D3. Using Theorem 4.1, we obtain

D3 k=1K𝔼q𝝉[πk(𝝉)]𝒩L(𝐱^𝐜^k,𝚪^k)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)d𝐱^
=k=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐀^k𝐜^k+𝐛^k,𝚺^k+𝐀^k𝚪^k𝐀^k).

Indeed, we made use of the following results for marginal and conditional Gaussians, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for 𝐱 and a conditional Gaussian distribution for 𝐲 given 𝐱 in the form

p(𝐱) =𝒩(𝐱𝝁,𝚪1),
p(𝐲𝐱) =𝒩(𝐲𝐀𝐱+𝐛,𝐋1),

then the marginal distribution of 𝐲 and the conditional distribution of 𝐱 given 𝐲 are given by

p(𝐲) =p(𝐲𝐱)p(𝐱)d𝐱=𝒩(𝐲𝐀𝝁+𝐛,𝐋1+𝐀𝚪1𝐀),
p(𝐱𝐲) =𝒩(𝐱𝚺[𝐀𝐋(𝐲𝐛)+𝚪𝝁],𝚺),𝚺=(𝚪+𝐀𝐋𝐀)1.

In our situation, the desired result is obtained via using 𝐲𝐲^,𝐀𝐀^k,𝐛𝐛^k,𝐋1=𝚺^k, 𝐱𝐱^,𝝁𝐜^k,𝚪1𝚪^k.

Finally, we obtain

p(𝐱^𝐲^,𝒳,𝒴)
k=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜^k,𝚪^k)k=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐀^k𝐜^k+𝐛^k,𝚺^k+𝐀^k𝚪^k𝐀^k)
=k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k),

where

gDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)=𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k)k=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k).

Here, we used the fact that p(𝐲^,𝐱^z^=k)=p(𝐱^𝐲^,z^=k)p(𝐲^z^=k), namely,

𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜^k,𝚪^k)
=𝒩L(𝐱^𝚺^k[𝐀^k𝚺^1k(𝐲^𝐛^k)+𝚪^1k𝐜^k],𝚺^k)𝒩D(𝐲^𝐀^k𝐜^k+𝐛^k,𝚺^k+𝐀^k𝚪^k𝐀^k)
=𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k)𝒩D(𝐲^𝐜^k,𝚪^k),

with

𝚺^k =(𝚪^1k+𝐀^k𝚺^1k𝐀^k)1,
𝐀^k =𝚺^k𝐀^k𝚺^1k,
𝐛^k =𝚺^k[𝚪^1k𝐜^k𝐀^k𝚺^1k𝐛^k],
𝐜^k =𝐀^k𝐜^k+𝐛^k,
𝚪^k =𝚺^k+𝐀^k𝚪^k𝐀^k.

When required, it is straightforward to approximate the expectation and covariance matrix of 𝐱^𝐲^,𝒳,𝒴 as follows:

𝔼[𝐱^𝐲^,𝒳,𝒴] (𝐱^𝐲^,𝒳,𝒴)k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k)d𝐱^
=k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)(𝐱^𝐲^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k)d𝐱^
=k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)(𝐀^k𝐲^+𝐛^k),
var[𝐱^𝐲^,𝒳,𝒴] =𝔼[(𝐱^𝐲^,𝒳,𝒴)(𝐱^𝐲^,𝒳,𝒴)]𝔼(𝐱^𝐲^,𝒳,𝒴)𝔼(𝐱^𝐲^,𝒳,𝒴)
(𝐱^𝐲^,𝒳,𝒴)(𝐱^𝐲^,𝒳,𝒴)k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k)d𝐱^
𝔼(𝐱^𝐲^,𝒳,𝒴)𝔼(𝐱^𝐲^,𝒳,𝒴)
k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)(𝐱^𝐲^,𝒳,𝒴)(𝐱^𝐲^,𝒳,𝒴)𝒩L(𝐱^𝐀^k𝐲^+𝐛^k,𝚺^k)d𝐱^
𝔼(𝐱^𝐲^,𝒳,𝒴)𝔼(𝐱^𝐲^,𝒳,𝒴)
k=1KgDk(𝐲^𝚯^,ϕ^,𝒳,𝒴)[𝚺^k+(𝐀^k𝐲^+𝐛^k)(𝐀^k𝐲^+𝐛^k)]
𝔼(𝐱^𝐲^,𝒳,𝒴)𝔼(𝐱^𝐲^,𝒳,𝒴),

where we used the following definitions

cov(𝐗,𝐘) =𝔼(𝐗𝐘)𝔼(𝐗)𝔼(𝐘),var(𝐗)=cov(𝐗,𝐗).

Appendix F BNP-GLLiM2: a model with an hyperprior on the gating parameters

F.1 VBEM for BNP-GLLiM2

A more general BNP-GLLiM model, referred to as BNP-GLLiM2 can be considered by specifying a prior on the gating parameters (𝐜k,𝚪k) as a normal-inverse-Wishart (NIW) distribution parameterized by 𝝆k=(𝐦k,λk,𝚿k,νk) with a PDF

p(𝐜k,𝚪k𝝆k)𝒩𝒲(𝐜k,𝚪k𝝆k)=𝒩(𝐜k𝐦k,λk1𝚪k)𝒲(𝚪k𝚿k,νk).

The assumptions on the other parameters are not changed, so that hyperparameters and parameters are now as follows:

ϕ=(s1,s2,a,(𝝆k,𝐀k,𝐛k,𝚺k)k), while 𝚯=(𝝉,α,σ,𝜽),𝜽=(𝜽k)k(𝐜k,𝚪k)k.

BNP-GLLiM2 can be represented graphically as in Figure 5.

𝝉zn𝐱n𝐲nN𝐜𝚪𝐦,λ𝚿,ναs1,s2σa𝐀,𝐛,𝚺
Figure 5: Graphical representation of BNP-GLLiM2: the plate denotes N i.i.d. observations, white-filled circles correspond to unobserved (latent) variables and random or unknown parameters represented in red, while grey-filled circles correspond to observed variables represented in green. Hyperparameters are represented in blue.

The joint distribution of the observed data 𝒳,𝒴 and all latent variables can be expressed hierarchically as

p(𝒴,𝒳,𝒵,𝚯;ϕ) =n=1Np(𝐲n𝐱n,zn,𝚯;ϕ)p(𝐱nzn,𝚯;ϕ)p(𝒵𝚯;ϕ)p(𝚯;ϕ)
=n=1Np(𝐲n𝐱n,zn;𝐀,𝐛,𝚺)p(𝐱nzn,𝐜,𝚪)p(𝒵𝝉)
kp(τkα,σ)p(α,σs1,s2,a)kp(𝐜k,𝚪k;𝝆k). (40)

Following the same idea as in Section 3, we only consider the truncated variational posterior of parameters 𝚯 as follows

q𝚯(𝚯)=qα,σ(α,σ)k=1K1qτk(τk)k=1Kq𝜽k(𝜽k). (41)

These forms of q𝒵 and q𝚯 lead to our four VB-E steps and three VB-M steps, summarized below with details in Appendix C. Set the initial value of ϕ to ϕ(0). Then, repeat iteratively the following steps. The iteration index is omitted in the update formulas for simplicity.

VB-E steps

Note that the VB-E-𝝉, VB-E-(α,σ) steps are the same as in Section 3. We only highlight the modified steps as follows.

We first consider the derivation of the update equation for the factor q𝒵(𝒵).

F.1.1 VB-E-𝒵 step

By using the mean-field approximation (9) and the truncation, see Section C.3 for more details, for all n[N] and for k[K], this step consists of computing

qzn(k)=ρnkl=1Kρnl. (42)

Here, given 𝒩n represents the neighbors of n, we define logρnk by

12{log|𝚺^k|+(𝐲n𝐀^k𝐱n𝐛^k)𝚺^k1(𝐲n𝐀^k𝐱n𝐛^k)+log|𝚿^k2|l=1Lψ(ν^k+(1l)2)+ν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k)+Lλ^k}+ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1[ψ(γ^l,2)ψ(γ^l,1+γ^l,2)]. (43)

Note that in the above formula, symbols (𝐦^k,λ^k,𝚿^k,ν^k) and (𝐀^k,𝐛^k,𝚺^k) are the hyperparameters Specifically defined in the following Sections F.1.2 and 3.2.

Proof of (43).

With respect to the VBEM for BNP-GLLiM2 model from Section F.1.1, it is almost similar to the previous step in Section C.3, except that we have to take into account the randomness of 𝐜 and 𝚪. Namely, we have

qzn(zn)
exp𝔼q𝚯[log(p(𝐲n𝐱n,zn;𝐀^zn,𝐛^zn,𝚺^zn)p(𝐱nzn,𝐜zn,𝚪zn)p(𝐳𝝉))]
=exp{logp(𝐲n𝐱n,zn;𝐀^zn,𝐛^zn,𝚺^zn)+𝔼q𝜽zn[logp(𝐱nzn,𝐜zn,𝚪zn)]+𝔼q𝝉[logπzn(𝝉)]}. (44)

Here, for zn=k, it holds that

𝔼q𝜽zn[logp(𝐱nzn,𝐜^zn,𝚪^zn)]=𝔼q𝜽zn[log𝒩L(𝐱n𝐜^k,𝚪^k)]
=L2log(2π)12𝔼q𝚪k[log|𝚪^k|]12𝔼q𝜽zn[(𝐱n𝐜^k)𝚪^1k(𝐱n𝐜^k)],

where we used the fact that

𝔼q𝚪k[log|𝚪^k|] =log|𝚿^k2|l=1Lψ(ν^k+(1l)2),
𝔼q𝜽zn[(𝐱n𝐜^k)𝚪^1k(𝐱n𝐜^k)] =ν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k)+Lλ^k.

Plugging in all of the above expression back into (44) yields the desired results in (43). ∎

F.1.2 VB-E-𝜽 step

This step is divided into K parts where the computation is similar to that in standard Bayesian GMM with a choice of conjugate prior. Hence, for each kK, the variational posterior is a Normal-inverse-Wishart density defined as

q𝜽k(𝐜k,𝚪k)=𝒩𝒲(𝐜k,𝚪k𝐦^k,λ^k,𝚿^k,ν^k). (45)

Here, the hyperparameters are updated as follows (see, e.g., Bishop (2006, Section 10.2.1)):

λ^k =λk+Nk,ν^k=νk+Nk,Nk=n=1Nqzn(k)
𝚿^k =𝚿k+Nk𝑺k+λkNkλk+Nk(𝐦k𝒄¯k)(𝐦k𝒄¯k),
𝐦^k =λk𝐦k+Nk𝒄¯kλk+Nk=λk𝐦k+Nk𝒄¯kλ^k, (46)
𝒄¯k =1Nkn=1Nqzn(k)𝐱n,
𝑺k =1Nkn=1Nqzn(k)(𝐱n𝒄¯k)(𝐱n𝒄¯k).

VB-M steps

The maximisation step consists of updating the hyperparameters ϕ=(s1,s2,a,(𝝆k,𝐀k,𝐛k,𝚺k)k[K]), where 𝝆k=(𝐦k,λk,𝚿k,νk),k[K], by maximizing the free energy, if they are not set heuristically:

ϕ(r)=argmaxϕ𝔼q𝒵(r)q𝝉(r)qα,σ(r)q𝜽(r)[logp(𝒴,𝒳,𝒵,𝝉,α,σ,𝜽;ϕ)]. (47)

The VB-M-step can therefore be divided into four independent sub-steps as listed below. From the conditional independence of (s1,s2,a,𝝆) and (𝒴,𝒳,𝒵) given (𝝉,α,σ,𝜽), the solutions for the VB-M-(s1,s2) (in the DP case) and VB-M-𝝆 steps are straightforward. Only the M-(s1,s2,a) step (in the PYP case) and (𝐀k,𝐛k,𝚺k)k[K] are more involved.

Note that the VB-M-(s1,s2,a), VB-M-(𝐀,𝐛,𝚺) steps are the same as in Section 3. We only highlight the modified step below.

F.1.3 VB-M-𝝆 step

This step divides into K sub-steps that involve again cross-entropies,

ρk(r)=argmaxρ𝔼qθk(r)[logp(𝐜k,𝚪k;ρk)]=ρ^k(r),

where ρ^k(r)=(λ^k(r),ν^k(r),𝚿^k(r),𝐦^k(r)) is given in (F.1.2).

F.1.4 ELBO for BNP-GLLiM2

Proposition F.1.

When σ=0, the ELBO in BNP-GLLiM2 is determined analytically as follows:

(q𝒵,q𝚯,ϕ^) =𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]+𝔼[logp(𝒳𝒵,𝚯;ϕ^)]+𝔼[logp(𝒵𝚯;ϕ^)]
+𝔼[logp(𝚯;ϕ^)]𝔼[logq(𝒵)]𝔼[logq(𝚯)]. (48)

Here, we have the following update formulas:

𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]=n=1Nk=1Kqzn(k)log𝒩D(𝐲n𝐀^k𝐱n+𝐛^k,𝚺^k), (49)
𝔼[logp(𝒳𝒵,𝚯;ϕ^)]
=12k=1KNk[log𝚪~kLlog(2π)Lλ^1kν^kTr(𝐒k𝚿^k1)ν^k(𝐱¯k𝐦^k)𝚿^k1(𝐱¯k𝐦^k)], (50)
𝔼[logp(𝒵𝚯;ϕ^)]=k=1KNk[ψ(γ^k,1)ψ(γ^k,1+γ^k,2)+l=1k1[ψ(γ^l,2)ψ(γ^l,1+γ^l,2)]], (51)
𝔼[logp(𝚯;ϕ^)]=k=1K1𝔼[logp(τkα)]+𝔼[logp(αs^1,s^2)]+k=1K𝔼[logp(𝐜k,𝚪k;𝝆^k)], (52)
𝔼[logp(τkα)]=s^1s^2s^2[ψ(γ^k,2)ψ(γ^k,1+γ^k,2)]+ψ(s^1)log(s^2),
𝔼[logp(αs^1,s^2)]=logΓ(s^1)+(s^11)ψ(s^1)+log(s^2)s^1,
𝔼[logp(𝐜k,𝚪k;𝝆^k)]=12Llog(λ^k2π)L2L2ν^k+logB(𝚿^k1,ν^k)ν^kL2log𝚪~k,
log𝚪~k=l=1Lψ(ν^k+1l2)+Llog2+log|𝚿^k|,
𝔼[logq(𝒵)]=n=1Nk=1Kqzn(k)logqzn(k), (53)
𝔼[logq(𝚯)]=𝔼[logqα,0(α)]+k=1K1𝔼[logqτk(τk)]+k=1K𝔼[logq𝐜k,𝚪k(𝐜k,𝚪k)], (54)
𝔼[logqα,0(α)]=logΓ(s^1)+(s^11)ψ(s^1)+log(s^2)s^1,
𝔼[logqτk(τk)]=l=12(γ^k,l1){ψ(γ^k,l)ψ(γ^k,1+γ^k,2)}+logΓ(γ^k,1+γ^k,2)Γ(γ^k,1)Γ(γ^k,2),
𝔼[logq𝐜k,𝚪k(𝐜k,𝚪k)]=L2logλ^k2πL2+logB(𝚿^k1,ν^k)ν^kL2log𝚪~kν^kL2.

F.2 Predictive conditional density for BNP-GLLiM2

F.2.1 Joint density

We first show how to compute the joint density p(𝐲^,𝐱^,𝒳,𝒴) via Theorem F.2, which is proved in Section F.4.

Theorem F.2.

We approximate the joint density of BNP-GLLiM2 by a mixture of product between Gaussian and Student’s t-distributions as follows:

p(𝐲^,𝐱^,𝒳,𝒴) k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k). (55)

Here, the positive semidefinite shape matrices of Student’s t-distributions are given by

𝐋k =(ν^k+1L)λ^k1+λ^k𝚿^k. (56)

F.2.2 Inverse conditional density

We then show how to approximate the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴). This predictive density in BNP-GLLiM2 is approximated by a MoE via Theorem F.3 with the proof in Section F.5.

Theorem F.3.

We approximate the inverse conditional density of BNP-GLLiM2 by a MoE as follows:

p(𝐲^𝐱^,𝒳,𝒴) k=1Kgk(𝐱^𝚯^,ϕ^,𝒳,𝒴)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k).

Here, the gating posteriors are defined as

gk(𝐱^𝚯^,ϕ^,𝒳,𝒴)=𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)l=1K𝔼q𝝉[πl(𝝉)]St(𝐱^𝐦^l,𝐋l,ν^l+1L),k[K].

Furthermore, for any k[K], it holds that

𝔼q𝝉[πk(𝝉)] =𝔼qτk[τk]l=1k1𝔼qτl[1τl],
𝔼qτk[τk] =γ^k,1γ^k,1+γ^k,2,𝔼qτk[1τk]=1𝔼qτk[τk]=γ^k,2γ^k,1+γ^k,2,
γ^k,1 =1𝔼qα,σ[σ]+Nk,γ^k,2=𝔼qα,σ[α]+k𝔼qα,σ[σ]+l=k+1KNl,Nk=n=1nqzn(k).

The prediction task is carried out via the following approximation

𝔼[𝐲^𝐱^,𝒳,𝒴] k=1Kgk(𝐱^𝚯^,ϕ^,𝒳,𝒴)[𝐀^k𝐱^+𝐛^k].

F.2.3 Forward conditional density

To deal with high-dimensional regression data, namely high-to-low regression, given the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴), we want to approximate the following forward conditional density via Theorem F.4, whose proof is provided in Section F.6.

Theorem F.4.

It holds that

p(𝐱^𝐲^,𝒳,𝒴) i=1Ik=1Kgki(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k(ηi)𝐲^+𝐛^k(ηi),𝚺^k(ηi)),

which is a mixture of Gaussian experts, where, for all k[K],i[I],

gki(𝐲^𝚯^,ϕ^,𝒳,𝒴) =𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k(ηi))Gam(ηiν^k+1L2,ν^k+1L2)i=1Il=1K𝔼q𝝉[πl(𝝉)]𝒩D(𝐲^𝐜^l,𝚪^l(ηi))Gam(ηiν^l+1L2,ν^l+1L2),
𝚺^k(ηi) =(ηi𝐋k+𝐀^k𝚺^1k𝐀^k)1,
𝐀^k(ηi) =𝚺^k(ηi)𝐀^k𝚺^1k,
𝐛^k(ηi) =𝚺^k(ηi)[η𝐋k𝐦^k𝐀^k𝚺^1k𝐛^k],
𝐜^k =𝐀^k𝐦^k+𝐛^k,
𝚪^k(ηi) =𝚺^k+𝐀^k(ηi𝐋k)1𝐀^k.

Here, ηi,i[I], are chosen via discretizing η-space, [0,Uη], into a grid, e.g., uniform. Note that for simplicity, we evaluate the integrand as a Riemann integral with a truncated value 0<Uη< and a number of point I for approximating the integration but we can use any scheme to approximate such 1-dimensional integration.

F.3 Proof of Proposition F.1

Using the sum and product rules for both discrete and continuous variables, the ELBO in BNP-GLLiM (20) is given by

(q𝒵,q𝚯,ϕ^) =𝔼q𝒵q𝚯[logp(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q𝚯(𝚯)]𝔼[logp(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q(𝚯)]
=𝒵q(𝒵)q(𝚯)log[p(𝒴,𝒳,𝒵,𝚯;ϕ^)q(𝒵)q(𝚯)]d𝒵d𝚯
=𝔼[logp(𝒴𝒳,𝒵,𝚯;ϕ^)]+𝔼[logp(𝒳𝒵,𝚯;ϕ^)] (57)
+𝔼[logp(𝒵𝚯;ϕ^)]+𝔼[logp(𝚯;ϕ^)]𝔼[logq(𝒵)]𝔼[logq(𝚯)].

Note that the proofs of (49), (51), (53) are the same as in the proof of Proposition 3.1.

Proof of (50)

𝔼[logp(𝒳𝒵,𝚯;ϕ^)]=𝔼[logn=1Np(𝐱nzn,𝚯;ϕ^)]=𝔼[logn=1Nk=1K𝒩L(𝐱n𝒄k,𝚪k)znk]
=n=1Nk=1K𝔼[znklog𝒩L(𝐱n𝒄k,𝚪k)]
=n=1Nk=1K𝔼q𝒵[znk]𝔼q𝒄k,𝚪k[log𝒩L(𝐱n𝒄k,𝚪k)]
=n=1Nk=1Kqzn(k)[L2log(2π)12𝔼[log|𝚪k|]12𝔼[(𝐱n𝒄k)𝚪k1(𝐱n𝒄k)]]
=k=1Kn=1Nqzn(k)[L2log(2π)12𝔼q𝚪k[log|𝚪k|]12𝔼q𝒄k,𝚪k[(𝐱n𝒄k)𝚪k1(𝐱n𝒄k)]](Lemma F.5)
=12k=1Kn=1Nqzn(k)[log𝚪~kLlog(2π)Lλ^1kν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k)]
=12k=1KNk[log𝚪~kLlog(2π)Lλ^1k]12k=1Kn=1Nqzn(k)[ν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k)]
=12k=1KNk[log𝚪~kLlog(2π)Lλ^1kν^kTr(𝐒k𝚿^k1)ν^k(𝐱¯k𝐦^k)𝚿^k1(𝐱¯k𝐦^k)]
(using (60) from Lemma F.5). (58)

To obtain (58), we have to use the following Lemma F.5.

Lemma F.5.

We can compute the expectations w.r.t. the variational distributions of the parameters as follows:

log𝚪~k𝔼q𝚪k[log|𝚪k|] =l=1Lψ(ν^k+1l2)+Llog2+log|𝚿^k|,
𝔼q𝒄k,𝚪k[(𝐱n𝒄k)𝚪k1(𝐱n𝒄k)] =Lλ^1k+ν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k). (59)

Furthermore, for each k[K], it holds that

n=1Nqzn(k)[ν^k(𝐱n𝐦^k)𝚿^k1(𝐱n𝐦^k)] =Nk[ν^kTr(𝐒k𝚿^k1)+ν^k(𝐱¯k𝐦^k)𝚿^k1(𝐱¯k𝐦^k)]. (60)

Proof of (52)

Given a chosen truncated value K, it holds that

𝔼q𝚯[logp(𝚯;ϕ^)] =k=1K1𝔼q𝚯[logp(τkα,σ)]+𝔼q𝚯[logp(α,σs^1,s^2,a^)]
+k=1K𝔼q𝚯[logp(𝐜k,𝚪k;𝝆^k)].

Note that 𝔼q𝚯[logp(τkα,σ)] and 𝔼q𝚯[logp(α,σs^1,s^2,a^)] are calculated in the same way as in Proposition 3.1.

Finally, we have to compute the remaining term

𝔼q𝚯[logp(𝐜k,𝚪k;𝝆^k)]
=𝔼q𝒄k,𝚪k[L2log(2π)12logλ^kL|𝚪k|12(𝒄k𝐦^k)(λ^k1𝚪k)1(𝒄k𝐦^k)]
+𝔼q𝚪k[log𝒲(𝚪k1𝚿^k1,ν^k)]
=12Llog(2π)+12Llogλ^k12𝔼q𝚪k[log|𝚪k|]12λ^k𝔼q𝒄k,𝚪k[(𝒄k𝐦^k)𝚪k1(𝒄k𝐦^k)]
+𝔼q𝚪k[logB(𝚿^k1,ν^k)+ν^kL12log|𝚪k1|12Tr(𝚿^k𝚪k1)]
=12Llog(λ^k2π)12𝔼q𝚪k[log|𝚪k|]12λ^k𝔼q𝒄k,𝚪k[(𝒄k𝐦^k)𝚪k1(𝒄k𝐦^k)]
+logB(𝚿^k1,ν^k)ν^kL12𝔼q𝚪k[log|𝚪k|]12Tr(𝚿^k𝔼q𝚪k[𝚪k1)]
=12Llog(λ^k2π)12λ^k[Lλ^1k+ν^k(𝐦^k𝐦^k)𝚿^k1(𝐦^k𝐦^k)]
+logB(𝚿^k1,ν^k)ν^kL2log𝚪~k12ν^kTr(𝚿^k𝚿^k1)(using Lemma F.5)
=12Llog(λ^k2π)L2L2ν^k+logB(𝚿^k1,ν^k)ν^kL2log𝚪~k,

where

log𝚪~k𝔼q𝚪k[log|𝚪k|] =l=1Lψ(ν^k+1l2)+Llog2+log|𝚿^k|.

Proof of (54)

We have

𝔼[logq(𝚯)] =𝔼[logqα,σ(α,σ)]+k=1K1𝔼[logqτk(τk)]+k=1K𝔼[logq𝐜k,𝚪k(𝐜k,𝚪k)].

Note that these terms involving expectations of the logs of the q distributions simply represent the negative entropies of those distributions. In particular, the first two terms are calculated in the same way as in Proposition 3.1.

Similarly, we obtain

𝔼[logq𝐜k,𝚪k(𝐜k,𝚪k)]
=𝔼[log𝒩L(𝐜k𝐦^k,λ^k1𝚪k)]+𝔼[log𝒲(𝚪k1𝚿^k1,ν^k)]
=H[𝒩L(𝐜k𝐦^k,λ^k1𝚪k)]H[𝒲(𝚪k1𝚿^k1,ν^k)]
=L2logλ^k2π+12𝔼[log|𝚪k|]L2+logB(𝚿^k1,ν^k)ν^kL12𝔼[log|𝚪k|]ν^kL2
=L2logλ^k2πL2+logB(𝚿^k1,ν^k)ν^kL2log𝚪~kν^kL2.

F.4 Proof of Theorem F.2

Recall that we defined 𝚯=(𝝉,α,σ,𝜽), 𝜽=(𝜽k)k(𝐜k,𝚪k)k. Then,

p(𝐲^,𝐱^,𝒳,𝒴) =𝐳^p(𝐲^𝐱^,𝐳^,𝚯,𝒳,𝒴)p(𝐱^𝐳^,𝚯,𝒳,𝒴)p(𝐳^𝚯,𝒳,𝒴)p(𝚯𝒳,𝒴)d𝚯
=𝐳^p(𝐲^𝐱^,𝐳^;𝑨,𝒃,𝚺)p(𝐱^𝐳^,𝒄,𝚪)p(𝐳^𝝉;β)p(𝚯𝒳,𝒴)d𝚯T1. (61)

Note that in (F.4), p(𝚯𝒳,𝒴) is in fact the (unknown) true posterior distribution of the parameters given a sample (𝒳,𝒴). Because the integrations w.r.t. true posterior distribution are intractable, we approximate the predictive conditional density by replacing the true posterior distribution p(𝚯𝒳,𝒴) with its truncated variational posterior of parameters 𝚯 given by

qΘ(𝚯𝒳,𝒴)=qα,σ(α,σ𝒳,𝒴)k=1K1qτk(τk𝒳,𝒴)k=1Kq𝜽k(𝜽k𝒳,𝒴).

Recall that the infinite state space for each zj is dealt with by choosing a truncation of the state space to a maximum label K (Blei and Jordan, 2006). In practice, this consists of assuming that the variational distributions qzn for n[N], satisfy qzn(k)=0 for k>K and that the variational distribution on τ also factorizes as q𝝉(𝝉)=k=1K1qτk(τk) with the additional condition that τK=1. In particular, here we choose K=K^, where K^ is estimated from some suitable procedures.

For simplicity, we consider the case when β=0,σ=0. Then, we obtain

T1 𝐳^p(𝐲^𝐱^,𝐳^;𝑨,𝒃,𝚺)p(𝐱^𝐳^,𝒄,𝚪)p(𝐳^𝝉)qΘ(𝚯𝒳,𝒴)d𝚯
=k=1𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜k,𝚪k)πk(𝝉)qΘ(𝚯𝒳,𝒴)d𝚯
k=1K𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜k,𝚪k)πk(𝝉)q𝝉(𝝉𝒳,𝒴)d𝝉
×qα,0(α𝒳,𝒴)dα=1k=1Kq𝜽k(𝒄k,𝚪k𝒳,𝒴)d𝒄d𝚪
=k=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜k,𝚪k)q𝜽k(𝒄k,𝚪k𝒳,𝒴)d𝒄kd𝚪k
×j=1,jkKq𝜽j(𝒄j,𝚪j𝒳,𝒴)j=1,jkKd𝒄jd𝚪j=1
=k=1K𝔼qτk[πk(𝝉)]𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐜k,𝚪k)q𝜽k(𝒄k,𝚪k𝒳,𝒴)d𝒄kd𝚪k=St(𝐱^𝐦^k,𝐋k,ν^k+1L)(Lemma F.6)
=k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k).

Here we used the following Lemma F.6

Lemma F.6.

For each k[K], it holds that

𝒩L(𝐱^𝐜k,𝚪k)q𝜽k(𝒄k,𝚪k𝒳,𝒴)d𝒄kd𝚪k=St(𝐱^𝐦^k,𝐋k,ν^k+1L).

Proof of Lemma F.6

By definition, we obtain

𝒩L(𝐱^𝒄k,𝚪k1)q(𝝅𝒳)q(𝐜k,𝚪k𝒳)d𝒄kd𝚪k
=𝒩L(𝐱^𝒄k,𝚪k1)𝒩L(𝐜k𝐦^k,(λ^k𝚪k)1,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)d𝒄kd𝚪k
=𝒩L(𝐱^𝒄k,𝚪k1)𝒩L(𝐜k𝐦^k,(λ^k𝚪k)1,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)d𝒄kd𝚪k
=[𝒩L(𝐱^𝒄k,𝚪k1)𝒩L(𝐜k𝐦^k,(λ^k𝚪k)1,𝒳)d𝒄k]𝒲(𝚪k𝚿^k,ν^k,𝒳)d𝚪k
=𝒩L(𝐱^𝐦^k,(1+λ^1k)𝚪1k,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)d𝚪k. (62)

When the size of the data set is large, i.e.,  N , this predictive distribution (F.4) becomes a mixture of Gaussians with component means 𝐦^k and precisions 𝐋k. In particular, we made use of the following results for marginal and conditional Gaussians, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for 𝐱 and a conditional Gaussian distribution for 𝐲 given 𝐱 in the form

p(𝐱) =𝒩(𝐱𝝁,𝚪1),
p(𝐲𝐱) =𝒩(𝐲𝐀𝐱+𝐛,𝐋1),

then the marginal distribution of 𝐲 and the conditional distribution of 𝐱 given 𝐲 are given by

p(𝐲) =𝒩(𝐲𝐀𝝁+𝐛,𝐋1+𝐀𝚪1𝐀),
p(𝐱𝐲) =𝒩(𝐱𝚺{𝐀𝐋(𝐲𝐛)+𝚪𝝁},𝚺),

where

𝚺=(𝚪+𝐀𝐋𝐀)1.

In our situation, via using 𝐲𝐱^,𝐱𝐜k,𝐀𝐈,𝐛𝟎,𝐋1=𝚪1k,𝝁𝐦^k,𝚪1(λ^k𝚪k)1, we obtain

p(𝐱^|𝚪1k,𝒳) =𝒩L(𝐱^𝒄k,𝚪k1)𝒩L(𝐜k𝐦^k,(λ^k𝚪k)1,𝒳)d𝒄k
=𝒩L(𝐱^𝐦^k,𝚪1k+(λ^k𝚪k)1,𝒳)
=𝒩L(𝐱^𝐦^k,(1+λ^kλ^k)𝚪1k,𝒳).

Notice that the Wishart distribution is a conjugate prior for the Gaussian distribution with known mean and unknown precision. Therefore, it holds that the product of

𝒩L(𝐱^𝐦^k,(1+λ^1k)𝚪1k,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)

is again a Wishart distribution without normalized. This can be verified by focusing on the dependency on 𝚪k. More precisely, by using the trace trick of quadratic form, (𝐱^𝐦^k)𝚪k(𝐱^𝐦^k)=Tr((𝐱^𝐦^k)(𝐱^𝐦^k)𝚪k), we obtain

𝒩L(𝐱^𝐦^k,(1+λ^1k)𝚪1k,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)
=B(𝚿^k,ν^k)(2π(1+λ^1k))L/2C(𝚿^k,ν^k,λ^k)|𝚪k|1/2+(ν^kL1)/2
×exp{12(1+λ^1k)(𝐱^𝐦^k)𝚪k(𝐱^𝐦^k)12Tr(𝚿^1k𝚪k)}
=C(𝚿^k,ν^k,λ^k)|𝚪k|(ν^k+1L1)/2exp{12Tr((1+λ^1k)1(𝐱^𝐦^k)(𝐱^𝐦^k)𝚪k+𝚿^1k𝚪k)}
=C(𝚿^k,ν^k,λ^k)|𝚪k|(ν^k+1L1)/2exp{12Tr{[(1+λ^1k)1(𝐱^𝐦^k)(𝐱^𝐦^k)+𝚿^1k]𝚪k}}
=C(𝚿^k,ν^k,λ^k)B(𝚿^k,ν^k)𝒲(𝚪k𝚿k,ν^k).

Here, ν^k=ν^k+1, and

𝚿k =[(1+λ^1k)1(𝐱^𝐦^k)(𝐱^𝐦^k)+𝚿^1k]1,
|𝚿^k|(ν^k+1)/2 =|(1+λ^1k)1(𝐱^𝐦^k)(𝐱^𝐦^k)+𝚿^1k|(ν^k+1)/2
=|𝚿^1k[(1+λ^1k)1𝚿^k(𝐱^𝐦^k)(𝐱^𝐦^k)+𝐈]|(ν^k+1)/2
=|𝚿^k|(ν^k+1)/2|(1+λ^1k)1𝚿^k(𝐱^𝐦^k)(𝐱^𝐦^k)+𝐈|(ν^k+1)/2
=|𝚿^k|(ν^k+1)/2[1+(1+λ^1k)1(𝐱^𝐦^k)𝚿^k(𝐱^𝐦^k)](ν^k+1)/2.

Via the normalization constant we have

𝒩L(𝐱^𝐦^k,(1+λ^1k)𝚪1k,𝒳)𝒲(𝚪k𝚿^k,ν^k,𝒳)d𝚪k
=C(𝚿^k,ν^k,λ^k)B(𝚿^k,ν^k)𝒲(𝚪k𝚿k,ν^k)d𝚪k=1=C(𝚿^k,ν^k,λ^k)B(𝚿^k,ν^k)=B(𝚿^k,ν^k)(2π(1+λ^1k))L/21B(𝚿^k,ν^k)
=1(2π(1+λ^1k))L/2|𝚿^k|ν^k/2(2ν^kL/2πL(L1)/4l=1LΓ(ν^k+1l2))1|𝚿^k|(ν^k+1)/2(2(ν^k+1)L/2πL(L1)/4l=1LΓ((ν^k+1)+1l2))1
=1(π(1+λ^1k))L/2|𝚿^k|ν^k/2Γ(ν^k+12)Γ(ν^k2)Γ(ν^k+2L2)|𝚿^k|(ν^k+1)/2Γ(ν^k2)Γ(ν^k12)Γ(ν^k+2L2)Γ(ν^k+1L2)
=Γ(ν^k+12)Γ(ν^k+1L2)πL/2|𝚿^k|ν^k/2(1+λ^1k)L/2|𝚿^k|(ν^k+1)/2[1+(1+λ^1k)1(𝐱^𝐦^k)𝚿^k(𝐱^𝐦^k)](ν^k+1)/2
=Γ(ν^k+12)Γ(ν^k+1L2)πL/2|𝚿^k|1/2(1+λ^1k)L/2[1+(1+λ^1k)1(𝐱^𝐦^k)𝚿^k(𝐱^𝐦^k)](ν^k+1)/2
=St(𝐱^𝐦^k,𝐋k,ν^k+1L).

Here,

𝐋k =(ν^k+1L)λ^k1+λ^k𝚿^k,

and Δ2 is the squared Mahalanobis distance defined by

Δ2 =(𝐱^𝐦^k)𝐋k(𝐱^𝐦^k).

Then, the last equality holds since we have

St(𝐱^𝐦^k,𝐋k,ν^k+1L)=Γ(ν^k+1L2+L2)Γ(ν^k+1L2)πL/2|𝐋k|1/2(ν^k+1L)L/2[1+Δ2ν^k+1L](ν^k+1L)/2L/2
=Γ(ν^k+12)Γ(ν^k+1L2)πL/2(ν^k+1L)L/2|𝚿^k|1/2(ν^k+1L)L/2(1+λ^1k)L/2[1+Δ2ν^k+1L](ν^k+1L)/2L/2
=Γ(ν^k+12)Γ(ν^k+1L2)πL/2|𝚿^k|1/2(1+λ^1k)L/2[1+(ν^k+1L)λ^k(1+λ^k)(𝐱^𝐦^k)𝚿^k(𝐱^𝐦^k)ν^k+1L](ν^k+1)/2
=Γ(ν^k+12)Γ(ν^k+1L2)πL/2|𝚿^k|1/2(1+λ^1k)L/2[1+(𝐱^𝐦^k)𝚿^k(𝐱^𝐦^k)(1+λ^1k)](ν^k+1)/2.

F.5 Proof of Theorem F.3

From the product rule of probability, we see that this conditional distribution can be evaluated from the joint and marginal distributions. Furthermore, by integrating out 𝐳^ and 𝚯, the predictive conditional density is then given by

p(𝐲^𝐱^,𝒳,𝒴) =p(𝐲^,𝐱^𝒳,𝒴)p(𝐱^𝒳,𝒴)=𝐳^p(𝐲^,𝐱^,𝐳^,𝚯𝒳,𝒴)d𝐳^d𝚯𝐳^p(𝐱^,𝐳^,𝚯𝒳,𝒴)d𝐳^d𝚯T1T2.

Next, with a similar step as in the proof of Theorem F.2, we also obtain

T2=k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L).

Therefore

p(𝐲^𝐱^,𝒳,𝒴)
k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)
=k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)
k=1Kgk(𝐱^𝚯^,ϕ^,𝒳,𝒴)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k),

which is a mixture of Gaussian experts since we have

gk(𝐱^𝚯^,ϕ^,𝒳,𝒴)=𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L),k[K],

belongs to a K1 dimensional probability simplex.

F.6 Proof of Theorem F.4

To deal with high-dimensional regression data, namely high-to-low regression, given the inverse conditional density p(𝐲^𝐱^,𝒳,𝒴), we want to compute the following forward conditional density

p(𝐱^𝐲^,𝒳,𝒴) =p(𝐱^,𝐲^𝒳,𝒴)p(𝐲^𝒳,𝒴)=p(𝐱^,𝐲^𝒳,𝒴)𝐱^p(𝐱^,𝐲^𝒳,𝒴)d𝐱^=T1𝐱^T1(𝐱^)d𝐱^T1T3.

Then, we have to compute or numerically approximate D3. Using Theorem F.2 and definition of Student’s t-distribution, we obtain

T3 =k=1K𝔼q𝝉[πk(𝝉)]Dk.

Then, by definition of Student’s t-distribution, it holds that

Dk =St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)d𝐱^
=0𝒩L(𝐱^𝐦^k,(η𝐋k)1)Gam(ην^k+1L2,ν^k+1L2)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)dηd𝐱^
=0𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐦^k,(η𝐋k)1)d𝐱^Gam(ην^k+1L2,ν^k+1L2)dη
=0𝒩D(𝐲^𝐀^k𝐦^k+𝐛^k,𝚺^k+η1𝐀^k𝐋k1𝐀^k)Gam(ην^k+1L2,ν^k+1L2)dη.

Furthermore, we used the fact that

𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐦^k,(η𝐋k)1)d𝐱^=𝒩D(𝐲^𝐀^k𝐦^k+𝐛^k,𝚺^k+𝐀^k(η𝐋k)1𝐀^k).

Indeed, we made use of the following results for marginal and conditional Gaussians, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for 𝐱 and a conditional Gaussian distribution for 𝐲 given 𝐱 in the form

p(𝐱) =𝒩(𝐱𝝁,𝚪1),
p(𝐲𝐱) =𝒩(𝐲𝐀𝐱+𝐛,𝐋1),

then the marginal distribution of 𝐲 and the conditional distribution of 𝐱 given 𝐲 are given by

p(𝐲) =p(𝐲𝐱)p(𝐱)d𝐱=𝒩(𝐲𝐀𝝁+𝐛,𝐋1+𝐀𝚪1𝐀),
p(𝐱𝐲) =𝒩(𝐱𝚺{𝐀𝐋(𝐲𝐛)+𝚪𝝁},𝚺),𝚺=(𝚪+𝐀𝐋𝐀)1.

In our situation, the desired result is obtained via using 𝐲𝐲^,𝐀𝐀^k,𝐛𝐛^k,𝐋1=𝚺^k, 𝐱𝐱^,𝝁𝐦^k,𝚪1(η𝐋k)1.

Therefore, we obtain

p(𝐱^𝐲^,𝒳,𝒴)
k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)k=1K𝔼q𝝉[πk(𝝉)]St(𝐱^𝐦^k,𝐋k,ν^k+1L)𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)d𝐱^
=k=1K𝔼q𝝉[πk(𝝉)]0𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐦^k,(η𝐋k)1)Gam(ην^k+1L2,ν^k+1L2)dηk=1K𝔼q𝝉[πk(𝝉)]0𝒩D(𝐲^𝐀^k𝐦^k+𝐛^k,𝚺^k+η1𝐀^k𝐋k1𝐀^k)Gam(ην^k+1L2,ν^k+1L2)dη
=k=1K𝔼q𝝉[πk(𝝉)]0𝒩D(𝐲^𝐜^k,𝚪^k(η))𝒩L(𝐱^𝐀^k(η)𝐲^+𝐛^k(η),𝚺^k(η))Gam(ην^k+1L2,ν^k+1L2)dηk=1K𝔼q𝝉[πk(𝝉)]0𝒩D(𝐲^𝐜^k,𝚪^k(η))Gam(ην^k+1L2,ν^k+1L2)dη
i=1Ik=1Kgki(𝐲^𝚯^,ϕ^,𝒳,𝒴)𝒩L(𝐱^𝐀^k(ηi)𝐲^+𝐛^k(ηi),𝚺^k(ηi)),

where, for all k[K],i[I],

gki(𝐲^𝚯^,ϕ^,𝒳,𝒴)=𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k(ηi))Gam(ηiν^k+1L2,ν^k+1L2)i=1Ik=1K𝔼q𝝉[πk(𝝉)]𝒩D(𝐲^𝐜^k,𝚪^k(ηi))Gam(ηiν^k+1L2,ν^k+1L2).

Here, we used the fact that p(𝐲^,𝐱^z^=k)=p(𝐱^𝐲^,z^=k)p(𝐲^z^=k), namely,

𝒩D(𝐲^𝐀^k𝐱^+𝐛^k,𝚺^k)𝒩L(𝐱^𝐦^k,(η𝐋k)1)
=𝒩L(𝐱^𝚺^k[𝐀^k𝚺^1k(𝐲^𝐛^k)+η𝐋k𝐦^k],𝚺^k)𝒩D(𝐲^𝐀^k𝐦^k+𝐛^k,𝚺^k+𝐀^k(η𝐋k)1𝐀^k)
=𝒩L(𝐱^𝐀^k(η)𝐲^+𝐛^k(η),𝚺^k(η))𝒩D(𝐲^𝐜^k,𝚪^k(η)),

where

𝚺^k(η) =(η𝐋k+𝐀^k𝚺^1k𝐀^k)1,
𝐀^k(η) =𝚺^k(η)𝐀^k𝚺^1k,
𝐛^k(η) =𝚺^k(η)[η𝐋k𝐦^k𝐀^k𝚺^1k𝐛^k],
𝐜^k =𝐀^k𝐦^k+𝐛^k,
𝚪^k(η) =𝚺^k+𝐀^k(η𝐋k)1𝐀^k.

The last approximation is deduced by using the fact that one simplistic strategy for evaluating integration would be to discretize η-space (1-dimensional) into a uniform grid and to evaluate the integrand as a Riemann integral with a truncated value 0<Uη< and a number of point I for approximating the integration.

Cite this paper

Please cite the published version. Venue: Journal of Nonparametric Statistics, Journal article (2024). DOI: 10.1080/10485252.2024.2426091. Official record: Taylor & Francis.

BibTeX
@article{nguyen2024bayesian,
  title     = {Bayesian nonparametric mixture of experts for inverse problems},
  author    = {Nguyen, TrungTin and Forbes, Florence and Arbel, Julyan and Nguyen, Hien Duy},
  journal   = {Journal of Nonparametric Statistics},
  pages     = {1--60},
  year      = {2024}, publisher = {Taylor \& Francis},
  doi       = {10.1080/10485252.2024.2426091},
}