<figure><img src="image_1.png" /></figure>

HAL
open science

# Bayesian nonparametric mixture of experts for inverse problems

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

<figure><img src="image_2.png" /></figure>

## To cite this version:

Trungtin Nguyen, Florence Forbes, Julyan Arbel, Hien Duy Nguyen. Bayesian nonparametric mixture of experts for inverse problems. 2023. hal-04015203v3

HAL Id: hal-04015203

https://hal.science/hal-04015203v3

Preprint submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Bayesian nonparametric mixture of experts for inverse problems

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

► To cite this version:

Trungtin Nguyen, Florence Forbes, Julyan Arbel, Hien Duy Nguyen. Bayesian nonparametric mixture of experts for inverse problems. Journal of Nonparametric Statistics, 2024, pp.1-60. 10.1080/10485252.2024.2426091 . hal-04015203v3

HAL Id: hal-04015203

https://hal.science/hal-04015203v3

Submitted on 17 Jun 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bayesian nonparametric mixture of experts for inverse problems

TrungTin Nguyen$^{a,d}$, Florence Forbes$^a$, Julyan Arbel$^a$, and Hien Duy Nguyen$^{b,c}$

$^a$Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France; $^b$School of Computing, Engineering and Mathematical Sciences, La Trobe University, Melbourne Australia; $^c$Institute of Mathematics for Industry, Kyushu University, Fukuoka Japan;

$^d$School of Mathematics and Physics, The University of Queensland, St. Lucia, Australia.

## ARTICLE HISTORY

Compiled June 17, 2024

## Word Limits

11578 words for the main text and appendices.

## 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

Bayesian 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
# Bayesian nonparametric mixture of experts for inverse problems

TrungTin Nguyen$^{a,d}$, Florence Forbes$^a$, Julyan Arbel$^a$, and Hien Duy Nguyen$^{b,c}$

$^a$Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France; $^b$School of Computing, Engineering and Mathematical Sciences, La Trobe University, Melbourne Australia; $^c$Institute of Mathematics for Industry, Kyushu University, Fukuoka Japan; $^d$School of Mathematics and Physics, The University of Queensland, St. Lucia, Australia.

## ARTICLE HISTORY

Compiled June 17, 2024

## Word Limits

11578 words for the main text and appendices.

## 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

Bayesian 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., 2024), 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. (2024b, 2023a, 2024a). 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 theories 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) 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 is performed using an information criterion, such as the Akaike information criterion (AIC; Akaike, 1974), the Bayesian information criterion (BIC; Schwarz, 1978, BIC-GLLIM; Forbes et al., 2022a,b) or the BIC-like approximation of integrated classification likelihood (ICL-BIC; Biernacki et al., 2000). However, an important limitation of these criteria is that they are only asymptotically valid (Westerhout et al., 2024). This means that there are no finite sample guarantees when using AIC, ICL-BIC or BIC 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,b, 2023d,c) recently established non-asymptotic risk bounds in the form of weak oracle inequalities, provided that lower bounds on the penalties hold, in
et al., 2021), planetary science (Kugler et al., 2022; Nguyen et al., 2024), 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. (2024b, 2023a, 2024a). 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 theories 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) 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 is performed using an information criterion, such as the Akaike information criterion (AIC; Akaike, 1974), the Bayesian information criterion (BIC; Schwarz, 1978, BIC-GLLIM; Forbes et al., 2022a,b) or the BIC-like approximation of integrated classification likelihood (ICL-BIC; Biernacki et al., 2000). However, an important limitation of these criteria is that they are only asymptotically valid (Westerhout et al., 2024). This means that there are no finite sample guarantees when using AIC, ICL-BIC or BIC 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,b, 2023d,c) 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, from Nguyen et al. (2022a), is based on the closed testing principle and leads to a sequential testing procedure that allows for confidence statements to be made regarding the order of a finite mixture model. These works lead to an optimal data-driven choice of the number of components in finite-sample settings. 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. 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
high-dimensional regression scenarios for a variety of MoE models, including GLLiM. Another approach, from Nguyen et al. (2022a), is based on the closed testing principle and leads to a sequential testing procedure that allows for confidence statements to be made regarding the order of a finite mixture model. These works lead to an optimal data-driven choice of the number of components in finite-sample settings. 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. 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, \dots, D$} is abbreviated as $[D]$ for $D \in \mathbb{N}^*$, where $\mathbb{N}^*$ denotes the positive natural numbers. The notation $\equiv$ 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, $[\cdot; \cdot]$ 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 to 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 $X = (X_l)_{l \in [L]}$ conditionally on a high-dimensional $Y = (Y_d)_{d \in [D]}$, where typically $D \gg L$. We follow an inverse regression framework as in e.g., Li (1991); Deleforge et al. (2015). Therefore, in training, the low-dimensional variable $X$ plays the role of the regressor, while the response $Y$ is a function of $X$, possibly corrupted by noise through inverse cPDF $p(Y | X; \psi)$, where $\psi$ is an inverse parameter. The low dimension of the regressor $X$ allows to drastically reduce the number of parameters to be estimated. In addition, the forward parameter $\psi^*$ and cPDF $p(X | Y; \psi^*)$ are tractable after estimating the inverse parameter $\psi$. Therefore, this density can be used to predict the low-dimensional response $x$ of a high-dimensional test point $y$. This inverse-then-forward regression strategy justifies the unconventional notation: $Y$ for the high-dimensional input and $X$ for the low-dimensional response. Here and subsequently, we refer to the low-dimensional data sample as $\mathcal{X} = \{x_n\}_{n \in [N]} \subset (\mathbb{R}^L)^N$, the high-dimensional data sample as $\mathcal{Y} = \{y_n\}_{n \in [N]} \subset (\mathbb{R}^D)^N$. We denote the observed values as $(x, y)$, which are independently
of MoE models for cPDFs.

**Notations.** Throughout this paper, {$1, \dots, D$} is abbreviated as $[D]$ for $D \in \mathbb{N}^*$, where $\mathbb{N}^*$ denotes the positive natural numbers. The notation $\equiv$ 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, $[\cdot; \cdot]$ 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 to 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 $X = (X_l)_{l \in [L]}$ conditionally on a high-dimensional $Y = (Y_d)_{d \in [D]}$, where typically $D \gg L$. We follow an inverse regression framework as in e.g., Li (1991); Deleforge et al. (2015). Therefore, in training, the low-dimensional variable $X$ plays the role of the regressor, while the response $Y$ is a function of $X$, possibly corrupted by noise through inverse cPDF $p(Y | X; \psi)$, where $\psi$ is an inverse parameter. The low dimension of the regressor $X$ allows to drastically reduce the number of parameters to be estimated. In addition, the forward parameter $\psi^*$ and cPDF $p(X | Y; \psi^*)$ are tractable after estimating the inverse parameter $\psi$. Therefore, this density can be used to predict the low-dimensional response $x$ of a high-dimensional test point $y$. This inverse-then-forward regression strategy justifies the unconventional notation: $Y$ for the high-dimensional input and $X$ for the low-dimensional response. Here and subsequently, we refer to the low-dimensional data sample as $\mathcal{X} = \{x_n\}_{n \in [N]} \subset (\mathbb{R}^L)^N$, the high-dimensional data sample as $\mathcal{Y} = \{y_n\}_{n \in [N]} \subset (\mathbb{R}^D)^N$. We denote the observed values as $(x, y)$, which are independently
and identically distributed (i.i.d.) samples from the random variables (X, Y).

## 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 X becomes the regressor as in Section 2.1. GLLiM then relies on a piecewise linear model in the following way. The high-dimensional response Y is approximated by the local affine mappings K:

$$ Y = \sum_{k=1}^{K} \mathbb{I}(Z = k) (A_k X + b_k + E_k). $$

Here, $\mathbb{I}$ is an indicator function and $Z$ is a latent variable that captures a cluster relationship, such that $Z = k$ if $Y$ comes from cluster $k \in [K]$. Matrices $A_k \in \mathbb{R}^{D \times L}$ and vectors $b_k \in \mathbb{R}^D$ define cluster-specific affine transformations. In addition, $E_k$ are error terms that capture both the reconstruction error due to the local affine approximations as well as the observation noise in $\mathbb{R}^D$.

Following the usual assumption that $E_k$ is a zero-mean Gaussian variable with a covariance matrix $\Sigma_k \in \mathbb{R}^{D \times D}$, it follows that

$$ p(y | \mathbf{x}, Z = k; \psi) = \mathcal{N}_D(y | A_k \mathbf{x} + \mathbf{b}_k, \Sigma_k), \quad (1) $$

where $\psi$ is the vector of model parameters, $\mathcal{N}_D(y; A_k \mathbf{x} + \mathbf{b}_k, \Sigma_k)$ is the Gaussian cPDF of dimension $D$. In order to enforce the affine transformations to be local, $\mathbf{X}$ is defined as a mixture of $K$ Gaussian components as follows:

$$ p(\mathbf{x} | Z = k; \psi) = \mathcal{N}_L(\mathbf{x} | \mathbf{c}_k, \Gamma_k), \quad p(Z = k; \psi) = \pi_k, $$

where $\mathbf{c}_k \in \mathbb{R}^L$, $\Gamma_k \in \mathbb{R}^{L \times L}$, and $\pi = (\pi_k)_{k \in [K]}$ belongs to a probability simplex, defined as $\{ (\pi_k)_{k \in [K]} \in (\mathbb{R}^+)^K, \sum_{k=1}^K \pi_k = 1 \}$. Then, via the conditional property of Gaussian variables and hierarchical decomposition, given any $\psi = (\pi_k, \mathbf{c}_k, \Gamma_k, A_k, \mathbf{b}_k, \Sigma_k)_{k \in [K]} \in \Psi$, we obtain the following inverse conditional density:

$$ p(y | \mathbf{x}; \psi) = \sum_{k=1}^{K} \frac{\pi_k \mathcal{N}_L(\mathbf{x} | \mathbf{c}_k, \Gamma_k)}{\sum_{l=1}^{K} \pi_l \mathcal{N}_L(\mathbf{x} | \mathbf{c}_l, \Gamma_l)} \mathcal{N}_D(y | A_k \mathbf{x} + \mathbf{b}_k, \Sigma_k). \quad (2) $$

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

$$ p(\mathbf{x} | y; \psi^*) = \sum_{k=1}^{K} \frac{\pi_k^* \mathcal{N}_D(y | \mathbf{c}_k^*, \Gamma_k^*)}{\sum_{l=1}^{K} \pi_l^* \mathcal{N}_D(y | \mathbf{c}_l^*, \Gamma_l^*)} \mathcal{N}_L(\mathbf{x} | A_k^* y + \mathbf{b}_k^*, \Sigma_k^*). \quad (3) $$
and identically distributed (i.i.d.) samples from the random variables (X, Y).

## 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 X becomes the regressor as in Section 2.1. GLLiM then relies on a piecewise linear model in the following way. The high-dimensional response Y is approximated by the local affine mappings K:

$$Y = \sum_{k=1}^{K} \mathbb{I}(Z = k) (A_k X + b_k + E_k).$$

Here, $\mathbb{I}$ is an indicator function and $Z$ is a latent variable that captures a cluster relationship, such that $Z = k$ if $Y$ comes from cluster $k \in [K]$. Matrices $A_k \in \mathbb{R}^{D \times L}$ and vectors $b_k \in \mathbb{R}^D$ define cluster-specific affine transformations. In addition, $E_k$ are error terms that capture both the reconstruction error due to the local affine approximations as well as the observation noise in $\mathbb{R}^D$.

Following the usual assumption that $E_k$ is a zero-mean Gaussian variable with a covariance matrix $\Sigma_k \in \mathbb{R}^{D \times D}$, it follows that

$$p(y | x, Z = k; \psi) = \mathcal{N}_D(y | A_k x + b_k, \Sigma_k), \quad (1)$$

where $\psi$ is the vector of model parameters, $\mathcal{N}_D(y; A_k x + b_k, \Sigma_k)$ is the Gaussian cPDF of dimension $D$. In order to enforce the affine transformations to be local, $X$ is defined as a mixture of $K$ Gaussian components as follows:

$$p(x | Z = k; \psi) = \mathcal{N}_L(x | c_k, \Gamma_k), \quad p(Z = k; \psi) = \pi_k,$$

where $c_k \in \mathbb{R}^L$, $\Gamma_k \in \mathbb{R}^{L \times L}$, and $\pi = (\pi_k)_{k \in [K]}$ belongs to a probability simplex, defined as $\{(\pi_k)_{k \in [K]} \in (\mathbb{R}^+)^K, \sum_{k=1}^K \pi_k = 1\}$. Then, via the conditional property of Gaussian variables and hierarchical decomposition, given any $\psi = (\pi_k, c_k, \Gamma_k, A_k, b_k, \Sigma_k)_{k \in [K]} \in \Psi$, we obtain the following inverse conditional density:

$$p(y | x; \psi) = \sum_{k=1}^{K} \frac{\pi_k \mathcal{N}_L(x | c_k, \Gamma_k)}{\sum_{l=1}^{K} \pi_l \mathcal{N}_L(x | c_l, \Gamma_l)} \mathcal{N}_D(y | A_k x + b_k, \Sigma_k). \quad (2)$$

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

$$p(x | y; \psi^*) = \sum_{k=1}^{K} \frac{\pi_k^* \mathcal{N}_D(y | c_k^*, \Gamma_k^*)}{\sum_{l=1}^{K} \pi_l^* \mathcal{N}_D(y | c_l^*, \Gamma_l^*)} \mathcal{N}_L(x | A_k^* y + b_k^*, \Sigma_k^*). \quad (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 $\psi^*$ in the forward cPDF, defined in (3), can then be deduced from $\psi$ in (2) via the following one-to-one correspondence:

$$ \begin{pmatrix} \pi_k \\ c_k \\ \Gamma_k \\ A_k \\ b_k \\ \Sigma_k \end{pmatrix}_{k \in [K]} \mapsto \begin{pmatrix} \pi_k^* \\ c_k^* \\ \Gamma_k^* \\ A_k^* \\ b_k^* \\ \Sigma_k^* \end{pmatrix}_{k \in [K]} = \begin{pmatrix} \pi_k \\ A_k c_k + b_k \\ \Sigma_k + A_k \Gamma_k A_k^\top \\ \Sigma_k^* A_k^\top \Sigma_k^{-1} \\ \Sigma_k^* (\Gamma_k^{-1} c_k - A_k^\top \Sigma_k^{-1} b_k) \\ (\Gamma_k^{-1} + A_k^\top \Sigma_k^{-1} A_k)^{-1} \end{pmatrix}_{k \in [K]} $$

**Remark 1 (GLLIM models are computationally efficient).** Without assuming anything about the structure of the parameters, the dimension of GLLiM is $\dim(\Psi) = K(1 + D(L+1) + D(D+1)/2 + L(L+1)/2 + L) - 1$. It is worth noting that $\dim(\Psi)$ 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 $D \gg L$. Furthermore, under the assumption that the $K$ transformations are affine, it is more realistic to make the assumption on the residual covariance matrices $\Sigma_k$ of the error vectors $E_k$ rather than on $\Gamma_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, $E_k$ can be modelled with equal isotropic Gaussian noise, so we have $\Sigma_k = \tilde{\sigma}^2 I_D, \forall k \in [K]$, with some positive $\tilde{\sigma}^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,030 $\gg$ 30,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 ($D \approx 1000 \gg L$). 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 $(y_n, x_n)$ within our BNP-GLLIM model:
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 $\psi^*$ in the forward cPDF, defined in (3), can then be deduced from $\psi$ in (2) via the following one-to-one correspondence:

$$ \begin{pmatrix} \pi_k \\ c_k \\ \Gamma_k \\ A_k \\ b_k \\ \Sigma_k \end{pmatrix}_{k \in [K]} \mapsto \begin{pmatrix} \pi_k^* \\ c_k^* \\ \Gamma_k^* \\ A_k^* \\ b_k^* \\ \Sigma_k^* \end{pmatrix}_{k \in [K]} = \begin{pmatrix} \pi_k \\ A_k c_k + b_k \\ \Sigma_k + A_k \Gamma_k A_k^\top \\ \Sigma_k^* A_k^\top \Sigma_k^{-1} \\ \Sigma_k^* (\Gamma_k^{-1} c_k - A_k^\top \Sigma_k^{-1} b_k) \\ (\Gamma_k^{-1} + A_k^\top \Sigma_k^{-1} A_k)^{-1} \end{pmatrix}_{k \in [K]} $$

**Remark 1 (GLLIM models are computationally efficient).** Without assuming anything about the structure of the parameters, the dimension of GLLiM is $\dim(\Psi) = K(1 + D(L + 1) + D(D + 1)/2 + L(L + 1)/2 + L) - 1$. It is worth noting that $\dim(\Psi)$ 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 $D \gg L$. Furthermore, under the assumption that the $K$ transformations are affine, it is more realistic to make the assumption on the residual covariance matrices $\Sigma_k$ of the error vectors $E_k$ rather than on $\Gamma_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, $E_k$ can be modelled with equal isotropic Gaussian noise, so we have $\Sigma_k = \tilde{\sigma}^2 I_D, \forall k \in [K]$, with some positive $\tilde{\sigma}^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,030 $\gg$ 30,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 ($D \approx 1000 \gg L$). 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 $(y_n, x_n)$ within our BNP-GLLIM model:
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

$$ (1) \quad \text{BNP prior: } G \sim \text{BNP}(\alpha, \sigma, G_0). $$

$$ (\alpha, \sigma) | s_1, s_2, a \sim \mathcal{SG}(\alpha | \sigma; s_1, s_2) p(\sigma | a) \equiv \text{Gam}(\alpha + \sigma | s_1, s_2) p(\sigma | a), \quad (4) $$

$$ \tau_k | \alpha, \sigma \stackrel{\text{ind}}{\sim} \text{Beta}(\tau_k | 1 - \sigma, \alpha + k\sigma), \quad k \in \mathbb{N}^*. \quad (5) $$

Define $\pi_k(\tau) = \tau_k \prod_{l=1}^{k-1} (1 - \tau_l)$, $k \in \mathbb{N}^*$, $\theta_k^* | G_0 \stackrel{\text{iid}}{\sim} G_0$, where $\theta_k^* \equiv (c_k, \Gamma_k, A_k, b_k, \Sigma_k)$, $k \in \mathbb{N}^*$, and finally define $G = \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^*}$.

$$ (2) \quad \text{BNP-GLLIM model: for each } n \in [N], y_n \sim \text{BNP-GLLIM}(s_1, s_2, a, G). $$

$$ \theta_n | G \stackrel{\text{iid}}{\sim} G, \quad \text{if } \theta_n = \theta_k^*, \text{ set } z_n = k; \quad (6) $$

$$ \mathbf{x}_n | z_n, c, \Gamma \stackrel{\text{iid}}{\sim} \mathcal{N}_L(\mathbf{x}_n | c_{z_n}, \Gamma_{z_n}), \quad (c, \Gamma) \equiv (c_k, \Gamma_k)_{k \in \mathbb{N}^*}; \quad (7) $$

$$ y_n | \mathbf{x}_n, z_n, \mathbf{A}, \mathbf{b}, \mathbf{\Sigma} \stackrel{\text{iid}}{\sim} \mathcal{N}_D(y_n | \mathbf{A}_{z_n} \mathbf{x}_n + \mathbf{b}_{z_n}, \mathbf{\Sigma}_{z_n}), \quad (\mathbf{A}, \mathbf{b}, \mathbf{\Sigma}) \equiv (\mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)_{k \in \mathbb{N}^*}. $$

### 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 $\mathcal{Z} = (z_n)_{n \in [N]}$ from the observed data $(\mathcal{Y}, \mathcal{X}) = (y_n, \mathbf{x}_n)_{n \in [N]}$, whose joint distribution $p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)$ is determined by a set of BNP prior parameters $\Theta = (\tau, \alpha, \sigma, \theta^*)$, namely the stick-breaking construction of Pitman-Yor process (PYP) (Pitman and Yor, 1997) in Appendix A; and by additional hyperparameters $\phi = (s_1, s_2, a)$. Then the desired joint distribution is given by:

$$ p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi) = \prod_{n=1}^{N} p(y_n | \mathbf{x}_n, z_n; \mathbf{A}, \mathbf{b}, \mathbf{\Sigma}) p(\mathbf{x}_n | z_n, \mathbf{c}, \Gamma) p(\mathcal{Z} | \tau) \\ \times \prod_{k \in \mathbb{N}^*} p(\tau_k | \alpha, \sigma) p(\alpha, \sigma | s_1, s_2, a) p(\theta^* | G_0). \quad (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(\mathcal{Z}, \Theta) = q_{\mathcal{Z}}(\mathcal{Z}) q_{\Theta}(\Theta)$. In particular, the intractable posterior
<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

$$ (1) \quad \text{BNP prior: } G \sim \text{BNP}(\alpha, \sigma, G_0). $$

$$ (\alpha, \sigma) | s_1, s_2, a \sim \mathcal{SG}(\alpha | \sigma; s_1, s_2) p(\sigma | a) \equiv \text{Gam}(\alpha + \sigma | s_1, s_2) p(\sigma | a), \quad (4) $$

$$ \tau_k | \alpha, \sigma \stackrel{\text{ind}}{\sim} \text{Beta}(\tau_k | 1 - \sigma, \alpha + k\sigma), \quad k \in \mathbb{N}^*. \quad (5) $$

Define $\pi_k(\tau) = \tau_k \prod_{l=1}^{k-1} (1 - \tau_l)$, $k \in \mathbb{N}^*$, $\theta_k^* | G_0 \stackrel{\text{iid}}{\sim} G_0$, where $\theta_k^* \equiv (c_k, \Gamma_k, A_k, b_k, \Sigma_k)$, $k \in \mathbb{N}^*$, and finally define $G = \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^*}$.

$$ (2) \quad \text{BNP-GLLIM model: for each } n \in [N], y_n \sim \text{BNP-GLLIM}(s_1, s_2, a, G). $$

$$ \theta_n | G \stackrel{\text{iid}}{\sim} G, \quad \text{if } \theta_n = \theta_k^*, \text{ set } z_n = k; \quad (6) $$

$$ x_n | z_n, c, \Gamma \stackrel{\text{iid}}{\sim} \mathcal{N}_L(x_n | c_{z_n}, \Gamma_{z_n}), \quad (c, \Gamma) \equiv (c_k, \Gamma_k)_{k \in \mathbb{N}^*}; \quad (7) $$

$$ y_n | x_n, z_n, A, b, \Sigma \stackrel{\text{iid}}{\sim} \mathcal{N}_D(y_n | A_{z_n}x_n + b_{z_n}, \Sigma_{z_n}), \quad (A, b, \Sigma) \equiv (A_k, b_k, \Sigma_k)_{k \in \mathbb{N}^*}. $$

### 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 $\mathcal{Z} = (z_n)_{n \in [N]}$ from the observed data $(\mathcal{Y}, \mathcal{X}) = (y_n, x_n)_{n \in [N]}$, whose joint distribution $p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)$ is determined by a set of BNP prior parameters $\Theta = (\tau, \alpha, \sigma, \theta^*)$, namely the stick-breaking construction of Pitman-Yor process (PYP) (Pitman and Yor, 1997) in Appendix A; and by additional hyperparameters $\phi = (s_1, s_2, a)$. Then the desired joint distribution is given by:

$$ p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi) = \prod_{n=1}^{N} p(y_n | x_n, z_n; A, b, \Sigma) p(x_n | z_n, c, \Gamma) p(\mathcal{Z} | \tau) \\ \times \prod_{k \in \mathbb{N}^*} p(\tau_k | \alpha, \sigma) p(\alpha, \sigma | s_1, s_2, a) p(\theta^* | G_0). \quad (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(\mathcal{Z}, \Theta) = q_{\mathcal{Z}}(\mathcal{Z})q_{\Theta}(\Theta)$. In particular, the intractable posterior
on $\mathcal{Z}$ is approximated as $q_{\mathcal{Z}}(\mathcal{Z})$ that factorizes so as to handle intractability, namely

$$ q_{\mathcal{Z}}(\mathcal{Z}) = \prod_{n=1}^{N} q_{z_n}(z_n). \qquad (9) $$

Then the infinite state space for each $z_j$ is dealt with by choosing a truncation of the state space to a maximum label $K \in \mathbb{N}^*$, see, e.g., Blei and Jordan (2006); Wang et al. (2011). In practice, this consists of assuming that the variational distributions $q_{z_n}$ for $n \in [N]$, satisfy $q_{z_n}(k) = 0$ for $k > K$ and that the variational distribution on $\tau$ also factorizes as $q_{\tau}(\tau) = \prod_{k=1}^{K-1} q_{\tau_k}(\tau_k)$ with the additional condition that $\tau_K = 1$. Thus the truncated variational posterior of parameters $\Theta$ is given by

$$ q_{\Theta}(\Theta) = q_{\alpha, \sigma}(\alpha, \sigma) \prod_{k=1}^{K-1} q_{\tau_k}(\tau_k) \prod_{k=1}^{K} q_{\theta_k^*}(\theta_k^*). $$

In practice, for tractability reasons, we have to restrict to $p(\theta_k^*|G_0) \propto 1$ and $q_{\theta_k^*}(\theta_k^*)$ to Dirac distributions, $q_{\theta_k^*} = \delta_{\theta_k^*}$, which is equivalent to treating the $\theta_k^*$ as fixed unknown hyperparameters, as illustrated graphically in Figure 1. These forms of $q_{\mathcal{Z}}$ and $q_{\Theta}$ 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 $\phi$ to $\phi^{(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 ($c_k, \Gamma_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 (B3) on the variable $\tau_k$. Thus, any terms that do not depend on $\tau_k$ can be included in the additive normalization constant. Then, given for $k \in [K-1]$ that $N_k = \sum_{n=1}^N q_{z_n}(k)$ corresponds to the weight of the cluster $k$, see more details in Appendix C.1, it holds that

$$ q_{\tau_k}(\tau_k) = \text{Beta}(\tau_k | \hat{\gamma}_{k,1}, \hat{\gamma}_{k,2}), \quad \text{where,} \\ \hat{\gamma}_{k,1} = 1 - \mathbb{E}_{q_{\alpha,\sigma}}[\sigma] + N_k, \quad \hat{\gamma}_{k,2} = \mathbb{E}_{q_{\alpha,\sigma}}[\alpha] + k\mathbb{E}_{q_{\alpha,\sigma}}[\sigma] + \sum_{l=k+1}^{K} N_l. $$

**VB-E-(α, σ) step.** The $(\alpha, \sigma)$ variational posterior is more complex, but has a simple form in the DP case ($\sigma = 0$). Specifically, we have to compute

$$ \hat{s}_1 = s_1 + K - 1, \quad \hat{s}_2 = s_2 - \sum_{k=1}^{K-1} \psi(\hat{\gamma}_{k,2}) - \psi(\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}), \qquad (10) $$

where $\psi(\cdot)$ is the digamma function defined by $\psi(\gamma) = \frac{d}{d\gamma} \log \Gamma(\gamma) = \frac{\Gamma'(\gamma)}{\Gamma(\gamma)}$.

When $\sigma = 0$ then $q_{\alpha,\sigma} = q_{\alpha,0}$ is a gamma distribution $\text{Gam}(\hat{s}_1, \hat{s}_2)$ with $\mathbb{E}_{q_{\alpha,\sigma}}[\alpha] = \hat{s}_1/\hat{s}_2$. Otherwise (PYP case), $q_{\alpha,\sigma}$ is only identified up to a normalizing constant but
on $\mathcal{Z}$ is approximated as $q_{\mathcal{Z}}(\mathcal{Z})$ that factorizes so as to handle intractability, namely

$$ q_{\mathcal{Z}}(\mathcal{Z}) = \prod_{n=1}^{N} q_{z_n}(z_n). \qquad (9) $$

Then the infinite state space for each $z_j$ is dealt with by choosing a truncation of the state space to a maximum label $K \in \mathbb{N}^*$, see, e.g., Blei and Jordan (2006); Wang et al. (2011). In practice, this consists of assuming that the variational distributions $q_{z_n}$ for $n \in [N]$, satisfy $q_{z_n}(k) = 0$ for $k > K$ and that the variational distribution on $\tau$ also factorizes as $q_{\tau}(\tau) = \prod_{k=1}^{K-1} q_{\tau_k}(\tau_k)$ with the additional condition that $\tau_K = 1$. Thus the truncated variational posterior of parameters $\Theta$ is given by

$$ q_{\Theta}(\Theta) = q_{\alpha, \sigma}(\alpha, \sigma) \prod_{k=1}^{K-1} q_{\tau_k}(\tau_k) \prod_{k=1}^{K} q_{\theta_k^*}(\theta_k^*). $$

In practice, for tractability reasons, we have to restrict to $p(\theta_k^*|G_0) \propto 1$ and $q_{\theta_k^*}(\theta_k^*)$ to Dirac distributions, $q_{\theta_k^*} = \delta_{\theta_k^*}$, which is equivalent to treating the $\theta_k^*$ as fixed unknown hyperparameters, as illustrated graphically in Figure 1. These forms of $q_{\mathcal{Z}}$ and $q_{\Theta}$ 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 $\phi$ to $\phi^{(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 ($c_k, \Gamma_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 (B3) on the variable $\tau_k$. Thus, any terms that do not depend on $\tau_k$ can be included in the additive normalization constant. Then, given for $k \in [K-1]$ that $N_k = \sum_{n=1}^{N} q_{z_n}(k)$ corresponds to the weight of the cluster $k$, see more details in Appendix C.1, it holds that

$$ q_{\tau_k}(\tau_k) = \text{Beta}(\tau_k | \hat{\gamma}_{k,1}, \hat{\gamma}_{k,2}), \quad \text{where,} \\ \hat{\gamma}_{k,1} = 1 - \mathbb{E}_{q_{\alpha,\sigma}}[\sigma] + N_k, \quad \hat{\gamma}_{k,2} = \mathbb{E}_{q_{\alpha,\sigma}}[\alpha] + k\mathbb{E}_{q_{\alpha,\sigma}}[\sigma] + \sum_{l=k+1}^{K} N_l. $$

**VB-E-(α, σ) step.** The $(\alpha, \sigma)$ variational posterior is more complex, but has a simple form in the DP case ($\sigma = 0$). Specifically, we have to compute

$$ \hat{s}_1 = s_1 + K - 1, \quad \hat{s}_2 = s_2 - \sum_{k=1}^{K-1} \psi(\hat{\gamma}_{k,2}) - \psi(\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}), \qquad (10) $$

where $\psi(\cdot)$ is the digamma function defined by $\psi(\gamma) = \frac{d}{d\gamma} \log \Gamma(\gamma) = \frac{\Gamma'(\gamma)}{\Gamma(\gamma)}$.

When $\sigma = 0$ then $q_{\alpha,\sigma} = q_{\alpha,0}$ is a gamma distribution $\text{Gam}(\hat{s}_1, \hat{s}_2)$ with $\mathbb{E}_{q_{\alpha,\sigma}}[\alpha] = \hat{s}_1/\hat{s}_2$. Otherwise (PYP case), $q_{\alpha,\sigma}$ is only identified up to a normalizing constant but
the required $\mathbb{E}_{q_{\alpha,\sigma}}[\alpha]$ and $\mathbb{E}_{q_{\alpha,\sigma}}[\sigma]$ can be computed by importance sampling, see Appendix C.2 for more details.

We next consider the derivation of the update equation for the factor $q_Z(\mathcal{Z})$.

**VB-E-**$\mathcal{Z}$ step. By using the mean-field approximation (9) and the truncation (see Appendix C.3), for all $n \in [N]$ and for all $k \in [K]$, this step consists in computing

$$q_{z_n}(k) = \frac{\rho_{nk}}{\sum_{l=1}^{K} \rho_{nl}}, \quad \text{where we define } \log \rho_{nk} \text{ as follows:} \qquad (11)$$

$$\log \rho_{nk} = -\frac{1}{2} \left\{ D \log(2\pi) + \log |\hat{\Sigma}_k| + (y_n - \hat{\mathbf{A}}_k x_n - \hat{\mathbf{b}}_k)^{\top} \hat{\Sigma}_k^{-1} (y_n - \hat{\mathbf{A}}_k x_n - \hat{\mathbf{b}}_k) \right. \\ + L \log(2\pi) + \log |\hat{\Gamma}_k| + (x_n - \hat{c}_k)^{\top} \hat{\Gamma}_k^{-1} (x_n - \hat{c}_k) \\ \left. + \psi(\hat{\gamma}_{k,1}) - \psi(\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}) + \sum_{l=1}^{k-1} [\psi(\hat{\gamma}_{l,2}) - \psi(\hat{\gamma}_{l,1} + \hat{\gamma}_{l,2})]. \right\}$$

Note that in the above formula, symbols $(\hat{c}_k, \hat{\Gamma}_k, \hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\Sigma}_k)$ are the hyperparameters. Specifically defined in the following Section 3.2. It is important to note that (11) provides assignment probabilities $q_{z_n}(k)$ rather than intermediate commitments to hard assignments of $z_n$. However, the hard assignments can be postponed to the end if desired to obtain a segmentation by the following MAP estimation $\hat{z}_n = \underset{k \in [K]}{\operatorname{argmax}} q_{z_n}(k)$.

## 3.2. VB-M steps

The maximisation step consists of updating the hyperparameters $\phi = (s_1, s_2, a, \theta_{[K]}^*)$, where $\theta_{[K]}^* = (\mathbf{c}_k, \Gamma_k, \mathbf{A}_k, \mathbf{b}_k, \Sigma_k)_{k \in [K]}$, by maximizing the free energy as follows:

$$\phi^{(r)} = \underset{\phi}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)} q_{\tau}^{(r)} q_{\alpha,\sigma}^{(r)}} [\log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \tau, \alpha, \sigma; \phi)].$$

The VB-M-step can therefore be divided into 4 independent sub-steps, as listed below. From the conditional independence of $(s_1, s_2, a)$ and $(\mathcal{Y}, \mathcal{X}, \mathcal{Z})$ given $(\tau, \alpha, \sigma)$, the solution for the VB-M-$(s_1, s_2)$ (in the DP case) step is straightforward. Only the M-$(s_1, s_2, a)$ (in the PYP case) and M-$\theta_{[K]}^*$ steps are more involved.

**VB-M-(s₁, s₂, a) step.** This step is straightforward in the DP case ($\sigma = 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,

$$(s_1, s_2)^{(r)} = \underset{(s_1, s_2)}{\operatorname{argmax}} \mathbb{E}_{q_{\alpha,0}^{(r)}} [\log p(\alpha | s_1, s_2)] = (\hat{s}_1^{(r)}, \hat{s}_2^{(r)}),$$

where $(\hat{s}_1^{(r)}, \hat{s}_2^{(r)})$ is given in (10). We can also solve this step numerically using importance sampling in the more general case of PYP ($\sigma \neq 0$). For more details, see Appendix A.7 in Lü et al. (2020).
the required $\mathbb{E}_{q_{\alpha,\sigma}}[\alpha]$ and $\mathbb{E}_{q_{\alpha,\sigma}}[\sigma]$ can be computed by importance sampling, see Appendix C.2 for more details.

We next consider the derivation of the update equation for the factor $q_Z(\mathcal{Z})$.

**VB-E-\mathcal{Z} step.** By using the mean-field approximation (9) and the truncation (see Appendix C.3), for all $n \in [N]$ and for all $k \in [K]$, this step consists in computing

$$q_{z_n}(k) = \frac{\rho_{nk}}{\sum_{l=1}^{K} \rho_{nl}}, \quad \text{where we define } \log \rho_{nk} \text{ as follows:} \quad (11)$$

$$\log \rho_{nk} = -\frac{1}{2} \left\{ D \log(2\pi) + \log |\hat{\Sigma}_k| + (y_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k)^{\top} \hat{\Sigma}_k^{-1} (y_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k) \right. \\ + L \log(2\pi) + \log |\hat{\Gamma}_k| + (\mathbf{x}_n - \hat{\mathbf{c}}_k)^{\top} \hat{\Gamma}_k^{-1} (\mathbf{x}_n - \hat{\mathbf{c}}_k) \\ \left. + \psi(\hat{\gamma}_{k,1}) - \psi(\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}) + \sum_{l=1}^{k-1} [\psi(\hat{\gamma}_{l,2}) - \psi(\hat{\gamma}_{l,1} + \hat{\gamma}_{l,2})]. \right\}$$

Note that in the above formula, symbols $(\hat{\mathbf{c}}_k, \hat{\mathbf{\Gamma}}_k, \hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)$ are the hyperparameters. Specifically defined in the following Section 3.2. It is important to note that (11) provides assignment probabilities $q_{z_n}(k)$ rather than intermediate commitments to hard assignments of $z_n$. However, the hard assignments can be postponed to the end if desired to obtain a segmentation by the following MAP estimation $\hat{z}_n = \operatorname*{argmax}_{k \in [K]} q_{z_n}(k)$.

## 3.2. VB-M steps

The maximisation step consists of updating the hyperparameters $\phi = (s_1, s_2, a, \theta^*_{[K]})$, where $\theta^*_{[K]} = (\mathbf{c}_k, \mathbf{\Gamma}_k, \mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)_{k \in [K]}$, by maximizing the free energy as follows:

$$\phi^{(r)} = \underset{\phi}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)} q_{\tau}^{(r)} q_{\alpha, \sigma}^{(r)}} [\log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \tau, \alpha, \sigma; \phi)].$$

The VB-M-step can therefore be divided into 4 independent sub-steps, as listed below. From the conditional independence of $(s_1, s_2, a)$ and $(\mathcal{Y}, \mathcal{X}, \mathcal{Z})$ given $(\tau, \alpha, \sigma)$, the solution for the VB-M-$(s_1, s_2)$ (in the DP case) step is straightforward. Only the M-$(s_1, s_2, a)$ (in the PYP case) and M-$\theta^*_{[K]}$ steps are more involved.

**VB-M-(s₁, s₂, a) step.** This step is straightforward in the DP case ($\sigma = 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,

$$(s_1, s_2)^{(r)} = \underset{(s_1, s_2)}{\operatorname{argmax}} \mathbb{E}_{q_{\alpha,0}^{(r)}} [\log p(\alpha | s_1, s_2)] = (\hat{s}_1^{(r)}, \hat{s}_2^{(r)}),$$

where $(\hat{s}_1^{(r)}, \hat{s}_2^{(r)})$ is given in (10). We can also solve this step numerically using importance sampling in the more general case of PYP ($\sigma \neq 0$). For more details, see Appendix A.7 in Lü et al. (2020).
**VB-M-(c, Γ) step.** This step is divided into solving K optimisation problems:

$$ (\hat{\mathbf{c}}_k, \hat{\Gamma}_k) \equiv (\hat{\mathbf{c}}_k^{(r)}, \hat{\Gamma}_k^{(r)}) = \underset{(\mathbf{c}_k, \mathbf{\Gamma}_k)}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} [\log p(\mathcal{X} | \mathcal{Z}; \mathbf{c}_k, \mathbf{\Gamma}_k)]. $$

We can then update the Gaussian gating parameters as follows:

$$ \hat{\mathbf{c}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{x}_n, \quad \hat{\mathbf{\Gamma}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{x}_n - \hat{\mathbf{c}}_k) (\mathbf{x}_n - \hat{\mathbf{c}}_k)^{\top}. $$

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

**VB-M-(A, b, Σ) step.** Using the same idea, this step is divided into K sub-steps, which include the following optimisation problems

$$ (\hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \equiv (\hat{\mathbf{A}}_k^{(r)}, \hat{\mathbf{b}}_k^{(r)}, \hat{\mathbf{\Sigma}}_k^{(r)}) = \underset{(\mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} [\log p(\mathcal{Y} | \mathcal{X}, \mathcal{Z}; \mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)]. $$

Given the following quantities:

$$ \bar{\mathbf{x}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{x}_n, \quad \mathbf{X}_k = \frac{1}{\sqrt{N_k}} \left( \sqrt{q_{z_1}(k)}(\mathbf{x}_1 - \bar{\mathbf{x}}_k), \dots, \sqrt{q_{z_N}(k)}(\mathbf{x}_N - \bar{\mathbf{x}}_k) \right), $$

$$ \bar{\mathbf{y}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{y}_n, \quad \mathbf{Y}_k = \frac{1}{\sqrt{N_k}} \left( \sqrt{q_{z_1}(k)}(\mathbf{y}_1 - \bar{\mathbf{y}}_k), \dots, \sqrt{q_{z_N}(k)}(\mathbf{y}_N - \bar{\mathbf{y}}_k) \right), $$

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

$$ \hat{\mathbf{A}}_k = \mathbf{Y}_k \mathbf{X}_k^\top (\mathbf{X}_k \mathbf{X}_k^\top)^{-1}, \quad \hat{\mathbf{b}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n), $$

$$ \hat{\mathbf{\Sigma}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k)(\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k)^{\top}. $$

### 3.3. Evidence lower-bound (ELBO)

Evaluating the ELBO in (B1) 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 (B1) 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 $\hat{\phi} = (\hat{s}_1, \hat{s}_2, \hat{a}, (\hat{\mathbf{c}}_k, \hat{\Gamma}_k, \hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)_{k \in \mathbb{N}^*})$. 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 $\hat{}$ on the hyperparameters $\hat{\phi}$ of q distribution. If $\sigma \neq 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
**VB-M-(c, Γ) step.** This step is divided into solving K optimisation problems:

$$ (\hat{\mathbf{c}}_k, \hat{\Gamma}_k) \equiv (\hat{\mathbf{c}}_k^{(r)}, \hat{\Gamma}_k^{(r)}) = \underset{(\mathbf{c}_k, \Gamma_k)}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} [\log p(\mathcal{X} | \mathcal{Z}; \mathbf{c}_k, \Gamma_k)]. $$

We can then update the Gaussian gating parameters as follows:

$$ \hat{\mathbf{c}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{x}_n, \quad \hat{\Gamma}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{x}_n - \hat{\mathbf{c}}_k) (\mathbf{x}_n - \hat{\mathbf{c}}_k)^{\top}. $$

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

**VB-M-(A, b, Σ) step.** Using the same idea, this step is divided into K sub-steps, which include the following optimisation problems

$$ (\hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \equiv (\hat{\mathbf{A}}_k^{(r)}, \hat{\mathbf{b}}_k^{(r)}, \hat{\mathbf{\Sigma}}_k^{(r)}) = \underset{(\mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} [\log p(\mathcal{Y} | \mathcal{X}, \mathcal{Z}; \mathbf{A}_k, \mathbf{b}_k, \mathbf{\Sigma}_k)]. $$

Given the following quantities:

$$ \bar{\mathbf{x}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{x}_n, \quad \mathbf{X}_k = \frac{1}{\sqrt{N_k}} \left( \sqrt{q_{z_1}(k)}(\mathbf{x}_1 - \bar{\mathbf{x}}_k), \dots, \sqrt{q_{z_N}(k)}(\mathbf{x}_N - \bar{\mathbf{x}}_k) \right), $$

$$ \bar{\mathbf{y}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) \mathbf{y}_n, \quad \mathbf{Y}_k = \frac{1}{\sqrt{N_k}} \left( \sqrt{q_{z_1}(k)}(\mathbf{y}_1 - \bar{\mathbf{y}}_k), \dots, \sqrt{q_{z_N}(k)}(\mathbf{y}_N - \bar{\mathbf{y}}_k) \right), $$

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

$$ \hat{\mathbf{A}}_k = \mathbf{Y}_k \mathbf{X}_k^\top (\mathbf{X}_k \mathbf{X}_k^\top)^{-1}, \quad \hat{\mathbf{b}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n), $$

$$ \hat{\mathbf{\Sigma}}_k = \frac{1}{N_k} \sum_{n=1}^{N} q_{z_n}(k) (\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k)(\mathbf{y}_n - \hat{\mathbf{A}}_k \mathbf{x}_n - \hat{\mathbf{b}}_k)^{\top}. $$

### 3.3. Evidence lower-bound (ELBO)

Evaluating the ELBO in (B1) 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 (B1) 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 $\hat{\phi} = (\hat{s}_1, \hat{s}_2, \hat{a}, (\hat{\mathbf{c}}_k, \hat{\Gamma}_k, \hat{\mathbf{A}}_k, \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)_{k \in \mathbb{N}^*})$. 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 $\hat{}$ on the hyperparameters $\hat{\phi}$ of q distribution. If $\sigma \neq 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 $\sigma = 0$, we can analytically compute the ELBO in the BNP-GLLIM via Proposition 3.1 which is proved in Appendix E.2.

**Proposition 3.1.** If $\sigma = 0$, the ELBO in the BNP-GLLIM is decomposed as follows:

$$ \mathcal{F}(q_Z, q_\Theta, \hat{\phi}) = \mathbb{E}[\log p(\mathcal{Y} | \mathcal{X}, \mathcal{Z}, \Theta; \hat{\phi})] + \mathbb{E}[\log p(\mathcal{X} | \mathcal{Z}, \Theta; \hat{\phi})] + \mathbb{E}[\log p(\mathcal{Z} | \Theta; \hat{\phi})] \\ + \mathbb{E}[\log p(\Theta; \hat{\phi})] - \mathbb{E}[\log q(\mathcal{Z})] - \mathbb{E}[\log q(\Theta)], \quad (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 $\mathbb{E}[\log q_{\alpha,0}(\alpha)] = \mathbb{E}[\log p(\alpha | \hat{s}_1, \hat{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 $(\hat{y}, \hat{x})$ of the observed variables. Note that there will be a corresponding latent variable $\hat{z}$ associated with these observations. If $\sigma \neq 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 $\sigma = 0$, we can analytically approximate such densities via several following theorems. In Theorems 4.1 to 4.3, the notation “$\approx$” 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 Appendix E.3.

**Theorem 4.1.** With $\hat{w} = [\hat{x}; \hat{y}]$, we have

$$ \mu_k \equiv \mathbb{E}[\hat{w}] = \begin{pmatrix} \hat{c}_k \\ \hat{A}_k \hat{c}_k + \hat{b}_k \end{pmatrix}, \quad V_k = \operatorname{cov}[\hat{w}] = \begin{pmatrix} \hat{\Gamma}_k & \hat{\Gamma}_k \hat{A}_k^\top \\ \hat{A}_k \hat{\Gamma}_k & \hat{\Sigma}_k + \hat{A}_k \hat{\Gamma}_k \hat{A}_k^\top \end{pmatrix}, \quad (13) $$

$$ p(\hat{y}, \hat{x}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} \mathbb{E}_{q_{\tau_k}}[\tau_k] \prod_{l=1}^{k-1} \mathbb{E}_{q_{\tau_l}}[1 - \tau_l] \mathcal{N}_{L+D}(\hat{w} | \mu_k, V_k), \quad (14) $$

$$ \mathbb{E}_{q_{\tau_k}}[\tau_k] = \frac{\hat{\gamma}_{k,1}}{\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}}, \quad \mathbb{E}_{q_{\tau_k}}[1 - \tau_k] = \frac{\hat{\gamma}_{k,2}}{\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}}. $$

### 4.2. Inverse conditional density

We then show how to approximate the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$. This predictive density in BNP-GLLIM is approximated by a GLLIM via Theorem 4.2 with the proof in Appendix E.4.
(2008); Luo and Sun (2017); Wu and Ma (2018); Nguyen and Bonilla (2014). When $\sigma = 0$, we can analytically compute the ELBO in the BNP-GLLIM via Proposition 3.1 which is proved in Appendix E.2.

**Proposition 3.1.** If $\sigma = 0$, the ELBO in the BNP-GLLIM is decomposed as follows:

$$ \mathcal{F}(q_Z, q_\Theta, \hat{\phi}) = \mathbb{E}[\log p(\mathcal{Y} | \mathcal{X}, \mathcal{Z}, \Theta; \hat{\phi})] + \mathbb{E}[\log p(\mathcal{X} | \mathcal{Z}, \Theta; \hat{\phi})] + \mathbb{E}[\log p(\mathcal{Z} | \Theta; \hat{\phi})] \\ + \mathbb{E}[\log p(\Theta; \hat{\phi})] - \mathbb{E}[\log q(\mathcal{Z})] - \mathbb{E}[\log q(\Theta)], \quad (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 $\mathbb{E}[\log q_{\alpha,0}(\alpha)] = \mathbb{E}[\log p(\alpha | \hat{s}_1, \hat{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 $(\hat{y}, \hat{x})$ of the observed variables. Note that there will be a corresponding latent variable $\hat{z}$ associated with these observations. If $\sigma \neq 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 $\sigma = 0$, we can analytically approximate such densities via several following theorems. In Theorems 4.1 to 4.3, the notation “$\approx$” 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 Appendix E.3.

**Theorem 4.1.** With $\hat{w} \equiv [\hat{x}; \hat{y}]$, we have

$$ \boldsymbol{\mu}_k \equiv \mathbb{E}[\hat{w}] = \begin{pmatrix} \hat{c}_k \\ \hat{A}_k \hat{c}_k + \hat{b}_k \end{pmatrix}, \quad \mathbf{V}_k = \operatorname{cov}[\hat{w}] = \begin{pmatrix} \hat{\Gamma}_k & \hat{\Gamma}_k \hat{A}_k^\top \\ \hat{A}_k \hat{\Gamma}_k & \hat{\Sigma}_k + \hat{A}_k \hat{\Gamma}_k \hat{A}_k^\top \end{pmatrix}, \quad (13) $$

$$ p(\hat{y}, \hat{x}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} \mathbb{E}_{q_{\tau_k}}[\tau_k] \prod_{l=1}^{k-1} \mathbb{E}_{q_{\tau_l}}[1 - \tau_l] \mathcal{N}_{L+D}(\hat{w} | \boldsymbol{\mu}_k, \mathbf{V}_k), \quad (14) $$

$$ \mathbb{E}_{q_{\tau_k}}[\tau_k] = \frac{\hat{\gamma}_{k,1}}{\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}}, \quad \mathbb{E}_{q_{\tau_k}}[1 - \tau_k] = \frac{\hat{\gamma}_{k,2}}{\hat{\gamma}_{k,1} + \hat{\gamma}_{k,2}}. $$

### 4.2. Inverse conditional density

We then show how to approximate the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$. This predictive density in BNP-GLLIM is approximated by a GLLIM via Theorem 4.2 with the proof in Appendix E.4.
**Theorem 4.2.** We approximate the inverse conditional density $p(\hat{y} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y})$ and its conditional expectation for prediction by:

$$p(\hat{y} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} g_{Lk} (\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k), \quad (15)$$

$$\mathbb{E}[\hat{y} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Lk} (\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) [\hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k], \text{ where}$$

$$g_{Lk}(\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_L(\hat{\mathbf{x}} | \hat{\mathbf{c}}_k, \hat{\mathbf{\Gamma}}_k)}{\sum_{l=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_l(\tau)] \mathcal{N}_L(\hat{\mathbf{x}} | \hat{\mathbf{c}}_l, \hat{\mathbf{\Gamma}}_l)}, \quad k \in [K].$$

### 4.3. Forward conditional density

Given the inverse conditional density $p(\hat{y} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y})$, we approximate the following forward conditional density via Theorem 4.3, whose proof is provided in Appendix 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(\hat{\mathbf{x}} | \hat{y}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) \mathcal{N}_L (\hat{\mathbf{x}} | \hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*, \hat{\mathbf{\Sigma}}_k^*),$$

$$\mathbb{E}[\hat{\mathbf{x}} | \hat{y}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*),$$

$$\mathrm{var}[\hat{\mathbf{x}} | \hat{y}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) [\hat{\mathbf{\Sigma}}_k^* + (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*) (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*)^{\top}] \\ - \mathbb{E}(\hat{\mathbf{x}} | \hat{y}, \mathcal{X}, \mathcal{Y}) \mathbb{E}(\hat{\mathbf{x}} | \hat{y}, \mathcal{X}, \mathcal{Y})^{\top}, \text{ where}$$

$$g_{Dk}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*)}{\sum_{l=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_l(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_l^*, \hat{\mathbf{\Gamma}}_l^*)}, \\ \hat{\mathbf{\Sigma}}_k^* = (\hat{\mathbf{\Gamma}}_k^{-1} + \hat{\mathbf{A}}_k^{\top} \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{A}}_k)^{-1}, \quad \hat{\mathbf{b}}_k^* = \hat{\mathbf{\Sigma}}_k^* [\hat{\mathbf{\Gamma}}_k^{-1} \hat{\mathbf{c}}_k - \hat{\mathbf{A}}_k^{\top} \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{b}}_k], \\ \hat{\mathbf{A}}_k^* = \hat{\mathbf{\Sigma}}_k^* \hat{\mathbf{A}}_k^{\top} \hat{\mathbf{\Sigma}}_k^{-1}, \quad \hat{\mathbf{c}}_k^* = \hat{\mathbf{A}}_k \hat{\mathbf{c}}_k + \hat{\mathbf{b}}_k, \quad \hat{\mathbf{\Gamma}}_k^* = \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k \hat{\mathbf{\Gamma}}_k \hat{\mathbf{A}}_k^{\top}.$$

## 5. Bayesian nonparametric model selection

**Notations.** A coupling between $\pi \equiv (\pi_k)_{k \in [K]}$ and $\pi^0 \equiv (\pi_l^0)_{l \in [K_0]}$ is a joint distribution $\mathbf{Q}$ on $[K] \times [K_0]$, which is expressed as a matrix $\mathbf{Q} = (q_{kl})_{k \in [K], l \in [K_0]} \in [0, 1]^{K \times K_0}$ with marginal probabilities $\sum_{k \in [K]} q_{kl} = \pi^0$ and $\sum_{l \in [K_0]} q_{kl} = \pi_k$, for any $k \in [K]$ and $l \in [K_0]$. We use $\mathcal{Q}(\pi, \pi^0)$ to denote the space of all such couplings. Regarding the space of mixing measures, let $\mathcal{E}_K \equiv \mathcal{E}_K(\Theta)$ and $\mathcal{O}_K \equiv \mathcal{O}_K(\Theta)$ respectively denote the space of all mixing measures with exactly and at most $K$ support points, all in some parameter space $\Theta$. Additionally, with $\mathcal{G} \equiv \mathcal{G}(\Theta) = \bigcup_{K \in \mathbb{N}_+} \mathcal{E}_K$ to denote the set of all discrete measures with finite supports on $\Theta$. Moreover, $\bar{\mathcal{G}}(\Theta)$ denotes the space of
**Theorem 4.2.** We approximate the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$ and its conditional expectation for prediction by:

$$p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} g_{Lk} (\hat{x} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k), \quad (15)$$

$$\mathbb{E}[\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Lk} (\hat{x} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) [\hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k], \text{ where}$$

$$g_{Lk}(\hat{x} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_L(\hat{x} | \hat{\mathbf{c}}_k, \hat{\Gamma}_k)}{\sum_{l=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_l(\tau)] \mathcal{N}_L(\hat{x} | \hat{\mathbf{c}}_l, \hat{\Gamma}_l)}, \quad k \in [K].$$

### 4.3. Forward conditional density

Given the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$, we approximate the following forward conditional density via Theorem 4.3, whose proof is provided in Appendix 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(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}) \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) \mathcal{N}_L (\hat{x} | \hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*, \hat{\mathbf{\Sigma}}_k^*),$$

$$\mathbb{E}[\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*),$$

$$\mathrm{var}[\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}] \approx \sum_{k=1}^{K} g_{Dk} (\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) [\hat{\mathbf{\Sigma}}_k^* + (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*) (\hat{\mathbf{A}}_k^* \hat{y} + \hat{\mathbf{b}}_k^*)^\top] \\ - \mathbb{E}(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}) \mathbb{E}(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y})^\top, \text{ where}$$

$$g_{Dk}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\Gamma}_k^*)}{\sum_{l=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_l(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_l^*, \hat{\Gamma}_l^*)}, \\ \hat{\mathbf{\Sigma}}_k^* = (\hat{\Gamma}_k^{-1} + \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{A}}_k)^{-1}, \quad \hat{\mathbf{b}}_k^* = \hat{\mathbf{\Sigma}}_k^* [\hat{\Gamma}_k^{-1} \hat{\mathbf{c}}_k - \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{b}}_k], \\ \hat{\mathbf{A}}_k^* = \hat{\mathbf{\Sigma}}_k^* \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1}, \quad \hat{\mathbf{c}}_k^* = \hat{\mathbf{A}}_k \hat{\mathbf{c}}_k + \hat{\mathbf{b}}_k, \quad \hat{\mathbf{\Gamma}}_k^* = \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k \hat{\mathbf{\Gamma}}_k \hat{\mathbf{A}}_k^\top.$$

## 5. Bayesian nonparametric model selection

**Notations.** A coupling between $\pi \equiv (\pi_k)_{k \in [K]}$ and $\pi^0 \equiv (\pi_l^0)_{l \in [K_0]}$ is a joint distribution $\mathbf{Q}$ on $[K] \times [K_0]$, which is expressed as a matrix $\mathbf{Q} = (q_{kl})_{k \in [K], l \in [K_0]} \in [0, 1]^{K \times K_0}$ with marginal probabilities $\sum_{k \in [K]} q_{kl} = \pi^0$ and $\sum_{l \in [K_0]} q_{kl} = \pi_k$, for any $k \in [K]$ and $l \in [K_0]$. We use $\mathcal{Q}(\pi, \pi^0)$ to denote the space of all such couplings. Regarding the space of mixing measures, let $\mathcal{E}_K \equiv \mathcal{E}_K(\Theta)$ and $\mathcal{O}_K \equiv \mathcal{O}_K(\Theta)$ respectively denote the space of all mixing measures with exactly and at most $K$ support points, all in some parameter space $\Theta$. Additionally, with $\mathcal{G} \equiv \mathcal{G}(\Theta) = \bigcup_{K \in \mathbb{N}_+} \mathcal{E}_K$ to denote the set of all discrete measures with finite supports on $\Theta$. Moreover, $\bar{\mathcal{G}}(\Theta)$ denotes the space of
all discrete measures (including those with countably infinite supports) on $\Theta$. Finally, $\mathcal{P}(\Theta)$ stands for the space of all probability measures on $\Theta$.

## 5.1. Posterior contraction rate in Bayesian infinite mixtures

**Problem setup.** We first recall the GMM where we have i.i.d. samples ($\mathbf{W}_n$)$_{n \in [N]} \equiv \mathcal{W}$ coming from a true but unknown distribution $P_{G_0}$ with given PDF

$$ p_{G_0} \equiv \int \mathcal{N}(\mathbf{w} | \boldsymbol{\theta}) dG_0(\boldsymbol{\theta}) = \sum_{k \in [K_0]} \pi_k^0 \mathcal{N}(\mathbf{w} | \boldsymbol{\theta}_k^0), \quad \boldsymbol{\theta}_k^0 \equiv (\boldsymbol{\mu}_k^0, \mathbf{V}_k^0), \quad (16) $$

where $G_0 = \sum_{k \in [K_0]} \pi_k^0 \delta_{\theta_k^0}$ is a true but unknown mixing distribution with exactly $K_0$ number of support points, where $K_0$ is also unknown. Furthermore, $\Theta$ is a chosen parameter space, to which we believe that the true parameters belong. In a well-specified setting, all support points of $G_0$ reside in $\Theta$, but this may not be the case in a mis-specified 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 $p_{G_0}$, where mixing measure $G_0$ has $K_0$ support atoms in compact parameter space $\Theta$.

A Bayesian modeller places a prior distribution $\Pi$ on a suitable subspace of $\bar{\mathcal{G}}(\Theta)$. Then, the posterior distribution over $G$ is given by:

$$ \Pi(G \in B | \mathcal{W}) \equiv \frac{\int_B \prod_{n=1}^N p_G(\mathbf{W}_n) d\Pi(G)}{\int_{\bar{\mathcal{G}}(\Theta)} \prod_{n=1}^N p_G(\mathbf{W}_n) d\Pi(G)}. $$

Here, the GMM $p_G$ is defined in (16) with $K \le \infty$ unknown number of support points. We are interested in the posterior contraction behaviour of $G$ toward $G_0$, in addition to recovering the true number of components $K_0$.

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 $\Theta$ is equipped with a metric $d$. The Wasserstein distance $W_r$ between two discrete measures $G = \sum_{k \in [K]} \pi_k \delta_{\theta_k}$ and $G_0 = \sum_{l \in [K_0]} \pi_l^0 \delta_{\theta_l^0}$ is given by

$$ W_r(G, G_0) = \inf_{\mathbf{Q} \in \mathcal{Q}(\pi, \pi^0)} \left[ \sum_{k \in [K], l \in [K_0]} q_{kl} \left[ d(\boldsymbol{\theta}_k, \boldsymbol{\theta}_l^0) \right]^r \right]^{1/r}, $$

where couplings $Q$ 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 $G_N \in \mathcal{O}_{K_0}$ converges to $G_0 \in \mathcal{E}_{K_0}$ under the $W_r$ metric at a rate $\omega_N = o(1)$ for some $r \ge 1$, then there
all discrete measures (including those with countably infinite supports) on $\Theta$. Finally, $\mathcal{P}(\Theta)$ stands for the space of all probability measures on $\Theta$.

## 5.1. Posterior contraction rate in Bayesian infinite mixtures

**Problem setup.** We first recall the GMM where we have i.i.d. samples ($\mathbf{W}_n$)$_{n \in [N]} \equiv \mathcal{W}$ coming from a true but unknown distribution $P_{G_0}$ with given PDF

$$ p_{G_0} \equiv \int \mathcal{N}(\mathbf{w} | \boldsymbol{\theta}) dG_0(\boldsymbol{\theta}) = \sum_{k \in [K_0]} \pi_k^0 \mathcal{N}(\mathbf{w} | \boldsymbol{\theta}_k^0), \quad \boldsymbol{\theta}_k^0 \equiv (\boldsymbol{\mu}_k^0, \mathbf{V}_k^0), \qquad (16) $$

where $G_0 = \sum_{k \in [K_0]} \pi_k^0 \delta_{\theta_k^0}$ is a true but unknown mixing distribution with exactly $K_0$ number of support points, where $K_0$ is also unknown. Furthermore, $\Theta$ is a chosen parameter space, to which we believe that the true parameters belong. In a well-specified setting, all support points of $G_0$ reside in $\Theta$, but this may not be the case in a mis-specified 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 $p_{G_0}$, where mixing measure $G_0$ has $K_0$ support atoms in compact parameter space $\Theta$.

A Bayesian modeller places a prior distribution $\Pi$ on a suitable subspace of $\bar{\mathcal{G}}(\Theta)$. Then, the posterior distribution over $G$ is given by:

$$ \Pi(G \in B | \mathcal{W}) \equiv \frac{\int_B \prod_{n=1}^N p_G(\mathbf{W}_n) d\Pi(G)}{\int_{\bar{\mathcal{G}}(\Theta)} \prod_{n=1}^N p_G(\mathbf{W}_n) d\Pi(G)}. $$

Here, the GMM $p_G$ is defined in (16) with $K \le \infty$ unknown number of support points. We are interested in the posterior contraction behaviour of $G$ toward $G_0$, in addition to recovering the true number of components $K_0$.

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 $\Theta$ is equipped with a metric $d$. The Wasserstein distance $W_r$ between two discrete measures $G = \sum_{k \in [K]} \pi_k \delta_{\theta_k}$ and $G_0 = \sum_{l \in [K_0]} \pi_l^0 \delta_{\theta_l^0}$ is given by

$$ W_r(G, G_0) = \inf_{\mathbf{Q} \in \mathcal{Q}(\pi, \pi^0)} \left[ \sum_{k \in [K], l \in [K_0]} q_{kl} \left[ d(\boldsymbol{\theta}_k, \boldsymbol{\theta}_l^0) \right]^r \right]^{1/r}, $$

where couplings $Q$ 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 $G_N \in \mathcal{O}_{K_0}$ converges to $G_0 \in \mathcal{E}_{K_0}$ under the $W_r$ metric at a rate $\omega_N = o(1)$ for some $r \ge 1$, then there
exists a subsequence of $G_N$ such that the set of atoms of $G_N$ converges to the $K_0$ atoms of $G_0$, up to a permutation of the atoms, at the same rate $\omega_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 $\Pi(G | W)$, where $W$ are i.i.d. samples of the mixing density $p_{G_0}$. Under certain conditions on the kernel density, it can be established that for some Wasserstein metric $W_r$,

$$ \Pi(G \in \bar{G}(\Theta) : W_r(G, G_0) \le \delta\omega_N | W) \xrightarrow{N \to \infty} 1 \text{ in } P_{G_0}\text{-probability,} \quad (17) $$

for all constant $\delta > 0$, while $\omega_N = o(1)$ is a vanishing rate. Thus $\omega_N$ can be assumed to be (slightly) slower than the actual rate of posterior contraction of the mixture measure. We can also write that $\omega_N$ is a rate such that, under the posterior distribution $\Pi(G | W)$, $W_r(G, G_0) = o_{P_{G_0}}(\omega_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 $(X, Y)$, see (2), with unconstrained parameters $\psi = (\pi_k, c_k, \Gamma_k, A_k, b_k, \Sigma_k)_{k \in [K]}$, is equivalent to a GMM on the joint variable $[X; Y]$ with unrestricted parameters, via Lemma 5.2, which is briefly proved in Appendix E.1.

**Lemma 5.2.** A GLLiM model on $(X, Y)$ with unconstrained parameters $\psi = (\pi_k, c_k, \Gamma_k, A_k, b_k, \Sigma_k)_{k \in [K]}$, defined in (2), is equivalent to a GMM on the joint variable $[X; Y] \equiv W$ with unconstrained parameters $\nu = \{\mu_k, V_k, \rho_k\}_{k=1}^K$, i.e.,

$$ p(w | \psi) = \sum_{k=1}^{K} \rho_k \mathcal{N}_{L+D}(w | \mu_k, V_k). $$

The parameter $\psi$ can be expressed as a function of $\nu$ by:

$$ \begin{aligned} \pi_k &= \rho_k, & c_k &= \mu_k^x, & \Gamma_k &= V_k^{xx}, & A_k &= V_k^{xy\top} (V_k^{xx})^{-1}, \\ b_k &= \mu_k^y - V_k^{xy\top} (V_k^{xx})^{-1} \mu_k^x, & \Sigma_k &= V_k^{yy} - V_k^{xy\top} (V_k^{xx})^{-1} V_k^{xy}. \end{aligned} \quad (18) $$

Here, we have defined

$$ \mu_k = \begin{bmatrix} \mu_k^x \\ \mu_k^y \end{bmatrix}, \quad V_k = \begin{bmatrix} V_k^{xx} & V_k^{xy} \\ V_k^{yx} & V_k^{yy} \end{bmatrix}. $$

Note that the symmetry $V_k^{\top} = V_k$ of the covariance matrix implies that $V_k^{xx}$ and $V_k^{yy}$ are symmetric, while $V_k^{xy\top} = V_k^{yx}$. The parameter vector $\nu$ can be expressed as a
exists a subsequence of $G_N$ such that the set of atoms of $G_N$ converges to the $K_0$ atoms of $G_0$, up to a permutation of the atoms, at the same rate $\omega_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 $\Pi(G | W)$, where $W$ are i.i.d. samples of the mixing density $p_{G_0}$. Under certain conditions on the kernel density, it can be established that for some Wasserstein metric $W_r$,

$$ \Pi(G \in \bar{G}(\Theta) : W_r(G, G_0) \le \delta\omega_N | W) \xrightarrow{N \to \infty} 1 \text{ in } P_{G_0}\text{-probability,} \quad (17) $$

for all constant $\delta > 0$, while $\omega_N = o(1)$ is a vanishing rate. Thus $\omega_N$ can be assumed to be (slightly) slower than the actual rate of posterior contraction of the mixture measure. We can also write that $\omega_N$ is a rate such that, under the posterior distribution $\Pi(G | W)$, $W_r(G, G_0) = o_{P_{G_0}}(\omega_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 $(X, Y)$, see (2), with unconstrained parameters $\psi = (\pi_k, c_k, \Gamma_k, A_k, b_k, \Sigma_k)_{k \in [K]}$, is equivalent to a GMM on the joint variable $[X; Y]$ with unrestricted parameters, via Lemma 5.2, which is briefly proved in Appendix E.1.

**Lemma 5.2.** A GLLiM model on $(X, Y)$ with unconstrained parameters $\psi = (\pi_k, c_k, \Gamma_k, A_k, b_k, \Sigma_k)_{k \in [K]}$, defined in (2), is equivalent to a GMM on the joint variable $[X; Y] \equiv W$ with unconstrained parameters $\nu = \{\mu_k, V_k, \rho_k\}_{k=1}^K$, i.e.,

$$ p(\mathbf{w} | \psi) = \sum_{k=1}^{K} \rho_k \mathcal{N}_{L+D}(\mathbf{w} | \mu_k, \mathbf{V}_k). $$

The parameter $\psi$ can be expressed as a function of $\nu$ by:

$$ \begin{aligned} \pi_k &= \rho_k, & c_k &= \mu_k^x, & \Gamma_k &= V_k^{xx}, & A_k &= V_k^{xy\top} (V_k^{xx})^{-1}, \\ b_k &= \mu_k^y - V_k^{xy\top} (V_k^{xx})^{-1} \mu_k^x, & \Sigma_k &= V_k^{yy} - V_k^{xy\top} (V_k^{xx})^{-1} V_k^{xy}. \end{aligned} \quad (18) $$

Here, we have defined

$$ \mu_k = \begin{bmatrix} \mu_k^x \\ \mu_k^y \end{bmatrix}, \quad V_k = \begin{bmatrix} V_k^{xx} & V_k^{xy} \\ V_k^{yx} & V_k^{yy} \end{bmatrix}. $$

Note that the symmetry $V_k^{\top} = V_k$ of the covariance matrix implies that $V_k^{xx}$ and $V_k^{yy}$ are symmetric, while $V_k^{xy\top} = V_k^{yx}$. The parameter vector $\nu$ can be expressed as a
function of $\psi$ by:

$$ \rho_k = \pi_k, \quad \mu_k = \begin{bmatrix} \mathbf{c}_k \\ \mathbf{A}_k \mathbf{c}_k + \mathbf{b}_k \end{bmatrix}, \quad \mathbf{V}_k = \begin{bmatrix} \mathbf{\Gamma}_k & (\mathbf{A}_k \mathbf{\Gamma}_k)^{\top} \\ \mathbf{A}_k \mathbf{\Gamma}_k & \mathbf{\Sigma}_k + \mathbf{A}_k \mathbf{\Gamma}_k \mathbf{A}_k^{\top} \end{bmatrix}. \quad (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 $\omega_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\omega_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 $\tilde{G}$ and the number of its supporting atoms $\tilde{K}$. See Algorithm 1 for a pseudo-code description.

**Algorithm 1 MTM Algorithm for BNP joint GMM**

**Require:** Posterior sample $G = \sum_{k=1}^{K} \pi_k \delta_{\theta_k}$, posterior contraction rate $\omega_N$ from (17), and a tuning parameter $c$.
**Ensure:** Discrete measure $\tilde{G}$ and its number of supporting atoms $\tilde{K}$.

**Stage 1: Merge procedure:**

1. Reorder atoms $\{\theta_k\}_{k \in [K]}$ by simple random sampling without replacement with corresponding weights $\{\pi_1, \pi_2, \dots\}$. Let $\tau_1, \tau_2, \dots$ denote the new indices, and set $\mathcal{E} = \{\tau_j\}_j$ as the existing set of atoms.

2. Sequentially for each index $\tau_j \in \mathcal{E}$, if there exists an index $\tau_i < \tau_j$ such that $d(\theta_{\tau_i}, \theta_{\tau_j}) \le \omega_N$, then update $\pi_{\tau_i} = \pi_{\tau_i} + \pi_{\tau_j}$, and remove $\tau_j$ from $\mathcal{E}$.

3. Collect $G' = \sum_{j: \tau_j \in \mathcal{E}} \pi_{\tau_j} \delta_{\theta_{\tau_j}}$. Then write $G'$ as $\sum_{k=1}^{K} q_k \delta_{\phi_k}$ so that $q_1 \ge q_2 \ge \dots$.

**Stage 2: Truncate-Merge procedure:**

4. Set $\mathcal{A} = \{i: q_i > (c\omega_N)^r\}$, $\mathcal{N} = \{i: q_i \le (c\omega_N)^r\}$.

5. For each index $i \in \mathcal{A}$, if there is $j \in \mathcal{A}$ such that $j < i$ and $q_i \|\phi_i - \phi_j\|^r \le (c\omega_N)^r$, then remove $i$ from $\mathcal{A}$ and add it to $\mathcal{N}$.

6. For each $i \in \mathcal{N}$, find atom $\phi_j$ among $j \in \mathcal{A}$ that is nearest to $\phi_i$, then update $q_j = q_j + q_i$.

7. Return $\tilde{G} = \sum_{j \in \mathcal{A}} q_j \delta_{\phi_j}$ and $\tilde{K} = |\mathcal{A}|$.

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, $\Pi(G | W)$ according to which the upper bound (17) holds for all $\delta > 0$. Let $\tilde{G}$ and $\tilde{K}$ be the outcome of Algorithm 1 applied to $G$, for an arbitrary constant $c > 0$. Then the following hold
function of $\psi$ by:

$$ \rho_k = \pi_k, \quad \mu_k = \begin{bmatrix} c_k \\ A_k c_k + b_k \end{bmatrix}, \quad V_k = \begin{bmatrix} \Gamma_k & (A_k \Gamma_k)^\top \\ A_k \Gamma_k & \Sigma_k + A_k \Gamma_k A_k^\top \end{bmatrix}. \quad (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 $\omega_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\omega_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 $\tilde{G}$ and the number of its supporting atoms $\tilde{K}$. See Algorithm 1 for a pseudo-code description.

<caption>**Algorithm 1 MTM Algorithm for BNP joint GMM**</caption>

**Require:** Posterior sample $G = \sum_{k=1}^K \pi_k \delta_{\theta_k}$, posterior contraction rate $\omega_N$ from (17), and a tuning parameter $c$.
**Ensure:** Discrete measure $\tilde{G}$ and its number of supporting atoms $\tilde{K}$.

**Stage 1: Merge procedure:**
1: Reorder atoms $\{\theta_k\}_{k \in [K]}$ by simple random sampling without replacement with corresponding weights $\{\pi_1, \pi_2, \dots\}$. Let $\tau_1, \tau_2, \dots$ denote the new indices, and set $\mathcal{E} = \{\tau_j\}_j$ as the existing set of atoms.
2: Sequentially for each index $\tau_j \in \mathcal{E}$, if there exists an index $\tau_i < \tau_j$ such that $d(\theta_{\tau_i}, \theta_{\tau_j}) \le \omega_N$, then update $\pi_{\tau_i} = \pi_{\tau_i} + \pi_{\tau_j}$, and remove $\tau_j$ from $\mathcal{E}$.
3: Collect $G' = \sum_{j: \tau_j \in \mathcal{E}} \pi_{\tau_j} \delta_{\theta_{\tau_j}}$. Then write $G'$ as $\sum_{k=1}^K q_k \delta_{\phi_k}$ so that $q_1 \ge q_2 \ge \dots$.

**Stage 2: Truncate-Merge procedure:**
4: Set $\mathcal{A} = \{i: q_i > (c\omega_N)^r\}$, $\mathcal{N} = \{i: q_i \le (c\omega_N)^r\}$.
5: For each index $i \in \mathcal{A}$, if there is $j \in \mathcal{A}$ such that $j < i$ and $q_i \|\phi_i - \phi_j\|^r \le (c\omega_N)^r$, then remove $i$ from $\mathcal{A}$ and add it to $\mathcal{N}$.
6: For each $i \in \mathcal{N}$, find atom $\phi_j$ among $j \in \mathcal{A}$ that is nearest to $\phi_i$, then update $q_j = q_j + q_i$.
7: Return $\tilde{G} = \sum_{j \in \mathcal{A}} q_j \delta_{\phi_j}$ and $\tilde{K} = |\mathcal{A}|$.

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, $\Pi(G | W)$ according to which the upper bound (17) holds for all $\delta > 0$. Let $\tilde{G}$ and $\tilde{K}$ be the outcome of Algorithm 1 applied to $G$, for an arbitrary constant $c > 0$. Then the following hold
(a) $\Pi(\tilde{K} = K_0 | W) \xrightarrow{N \to \infty} 1$ in $P_{G_0}$-probability.

(b) For all $\delta > 0$, $\Pi(G \in \bar{G}(\Theta) : W_r(\tilde{G}, G_0) \le \delta\omega_N | W) \xrightarrow{N \to \infty} 1$ in $P_{G_0}$-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 $\omega_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.

(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 $\Theta$, e.g., $d(\theta_{\tau_i}, \theta_{\tau_j}) = \|\mu_{\tau_i} - \mu_{\tau_j}\| + \|V_{\tau_i} - V_{\tau_j}\|$. Here, $\|\cdot\|$ denotes either the $l_2$-norm elements in $\mathbb{R}^{L+D}$ or the entrywise $l_2$-norm for matrices in $\mathbb{R}^{(L+D)\times(L+D)}$.

(ii) In practice, one may not have a mixing measure $G$ sampled from the posterior $\Pi(\cdot | W)$, 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 $F_N$ from the variational posterior $G_V = \sum_{k=1}^K \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \delta_{\theta_k}$. Here, $\mathbb{E}_{q_{\tau}}[\pi_k(\tau)]$ and $\theta_k = (\mu_k, V_k)$ are defined in (14) and (13), respectively. However, as long as $F_N$ is sufficiently close to $G$ in the sense that $W_r(F_N, G) \lesssim W_r(G, G_0)$, we can still apply the MTM algorithm to $F_N$, 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 available<fnref n="1" /> 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

<footnote n="1">¹https://github.com/Trung-TinNGUYEN/BNP-GLLiM</footnote>
(a) $\Pi(\tilde{K} = K_0 | W) \xrightarrow{N \to \infty} 1$ in $P_{G_0}$-probability.

(b) For all $\delta > 0$, $\Pi(G \in \bar{G}(\Theta) : W_r(\tilde{G}, G_0) \le \delta\omega_N | W) \xrightarrow{N \to \infty} 1$ in $P_{G_0}$-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 $\omega_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.

(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 $\Theta$, e.g., $d(\theta_{\tau_i}, \theta_{\tau_j}) = \|\mu_{\tau_i} - \mu_{\tau_j}\| + \|V_{\tau_i} - V_{\tau_j}\|$. Here, $\|\cdot\|$ denotes either the $l_2$-norm elements in $\mathbb{R}^{L+D}$ or the entrywise $l_2$-norm for matrices in $\mathbb{R}^{(L+D)\times(L+D)}$.

(ii) In practice, one may not have a mixing measure $G$ sampled from the posterior $\Pi(\cdot | W)$, 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 $F_N$ from the variational posterior $G_V = \sum_{k=1}^K \mathbb{E}_{q_\tau} [\pi_k(\tau)] \delta_{\theta_k}$. Here, $\mathbb{E}_{q_\tau} [\pi_k(\tau)]$ and $\theta_k = (\mu_k, V_k)$ are defined in (14) and (13), respectively. However, as long as $F_N$ is sufficiently close to $G$ in the sense that $W_r(F_N, G) \lesssim W_r(G, G_0)$, we can still apply the MTM algorithm to $F_N$, 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 available<fnref n="1" /> 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

<footnote n="1">¹https://github.com/Trung-TinNGUYEN/BNP-GLLiM</footnote>
following true inverse cPDF from GLLiM model as follows:

$$s_0(\mathbf{y} | \mathbf{x}; \psi^0) = \sum_{k=1}^{K_0} \frac{\pi_k^0 \mathcal{N}_L(\mathbf{x} | \mathbf{c}_k^0, \mathbf{\Gamma}_k^0)}{\sum_{l=1}^{K} \pi_l^0 \mathcal{N}_L(\mathbf{x} | \mathbf{c}_l^0, \mathbf{\Gamma}_l^0)} \mathcal{N}_D(\mathbf{y} | \mathbf{A}_k^0 \mathbf{x} + \mathbf{b}_k^0, \mathbf{\Sigma}_k^0).$$

Here $K_0 = 3$, $L = D = 1$, and $\psi^0 = (\pi^0, \mathbf{c}^0, \mathbf{\Gamma}^0, \mathbf{A}^0, \mathbf{b}^0, \mathbf{\Sigma}^0)$, where

$$\begin{align*}
\pi^0 &= (0.3, 0.4, 0.3), & \mathbf{c}^0 &= (1, 0.05, -1), & \mathbf{\Gamma}^0 &= (0.1, 0.2, 0.1), \\
\mathbf{A}^0 &= (-15, 3, 15), & \mathbf{b}^0 &= (-2, 1, -2), & \mathbf{\Sigma}^0 &= (0.5, 0.3, 0.5).
\end{align*}$$

Figure 2a shows typical $N = 1,000$ realisations of the true inverse cPDF from GLLiM, representing a $\pi$-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.

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 $G_V = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \delta_{\theta_k}$, where $\mathbb{E}_{q_{\tau}}[\pi_k(\tau)]$ and $\theta_k = (\mu_k, V_k)$ are defined in (14) and (13), respectively. We know that for some constant $\tilde{C}$, which depends on the covariance matrix $V_k^0$, the location parameters $\mu_k^0$ and the weights $\pi_k^0$, the contraction rate of mixing measures under the location Gaussian DP-MM is $\tilde{C}(\log(N))^{-1/2}$ with respect to the $W_2$-norm. Similar to Guha et al. (2021), our first attempt to choose $w_N$ to satisfy (17) is $w_N = (\frac{\log(\log(N))}{\log(N)})^{1/2}$. In fact, we can choose any $w_N$, as long as $\frac{w_N}{\log(N)^{-1/2}} \to \infty$, in order for $w_N$ to satisfy (17).

Since we only work with finite sample $N$, it is not expected that the posterior probability for $K_{MTM} = K_0$ is close to 1 and the input $c$ to Algorithm 1 should be chosen so that $\frac{\tilde{C}}{(\log(\log(N)))^{1/2}} \le c$. Furthermore, based on Equation (26) from Guha et al. (2021), with a useful lower bound on the posterior mass the mode, for any $1 > \epsilon > 0$, $(1 - \epsilon) \left(1 - \sum_{k=1}^3 \frac{c_0^{r/2}}{\pi_k^0}\right) > \frac{1}{2}$, we hope to identify $K_0$ via the posterior mode with a reasonable estimate. To guarantee $K = K_0$ consistently using the posterior mode safety, we have to choose $c < c_0$, with $c_0$ satisfying

$$(1 - \epsilon) \left(1 - \sum_{k=1}^{3} \frac{c_0^{r/2}}{\pi_k^0}\right) > \frac{1}{2} \Leftrightarrow \frac{1 - 2\epsilon}{2(1 - \epsilon)} \left(\sum_{k=1}^{3} \frac{1}{\pi_k^0}\right)^{-1} > c_0^{r/2} > c^{r/2}.$$
following true inverse cPDF from GLLiM model as follows:

$$s_0(\mathbf{y} | \mathbf{x}; \psi^0) = \sum_{k=1}^{K_0} \frac{\pi_k^0 \mathcal{N}_L(\mathbf{x} | \mathbf{c}_k^0, \mathbf{\Gamma}_k^0)}{\sum_{l=1}^{K} \pi_l^0 \mathcal{N}_L(\mathbf{x} | \mathbf{c}_l^0, \mathbf{\Gamma}_l^0)} \mathcal{N}_D(\mathbf{y} | \mathbf{A}_k^0 \mathbf{x} + \mathbf{b}_k^0, \mathbf{\Sigma}_k^0).$$

Here $K_0 = 3$, $L = D = 1$, and $\psi^0 = (\pi^0, \mathbf{c}^0, \mathbf{\Gamma}^0, \mathbf{A}^0, \mathbf{b}^0, \mathbf{\Sigma}^0)$, where

$$\begin{align*}
\pi^0 &= (0.3, 0.4, 0.3), & \mathbf{c}^0 &= (1, 0.05, -1), & \mathbf{\Gamma}^0 &= (0.1, 0.2, 0.1), \\
\mathbf{A}^0 &= (-15, 3, 15), & \mathbf{b}^0 &= (-2, 1, -2), & \mathbf{\Sigma}^0 &= (0.5, 0.3, 0.5).
\end{align*}$$

Figure 2a shows typical $N = 1,000$ realisations of the true inverse cPDF from GLLiM, representing a $\pi$-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.

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 $G_V = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \delta_{\theta_k}$, where $\mathbb{E}_{q_{\tau}}[\pi_k(\tau)]$ and $\theta_k = (\mu_k, \mathbf{V}_k)$ are defined in (14) and (13), respectively. We know that for some constant $\tilde{C}$, which depends on the covariance matrix $\mathbf{V}_k^0$, the location parameters $\mu_k^0$ and the weights $\pi_k^0$, the contraction rate of mixing measures under the location Gaussian DP-MM is $\tilde{C}(\log(N))^{-1/2}$ with respect to the $W_2$-norm. Similar to Guha et al. (2021), our first attempt to choose $w_N$ to satisfy (17) is $w_N = (\frac{\log(\log(N))}{\log(N)})^{1/2}$. In fact, we can choose any $w_N$, as long as $\frac{w_N}{\log(N)^{-1/2}} \to \infty$, in order for $w_N$ to satisfy (17).

Since we only work with finite sample $N$, it is not expected that the posterior probability for $K_{MTM} = K_0$ is close to 1 and the input $c$ to Algorithm 1 should be chosen so that $\frac{\tilde{C}}{(\log(\log(N)))^{1/2}} \le c$. Furthermore, based on Equation (26) from Guha et al. (2021), with a useful lower bound on the posterior mass the mode, for any $1 > \epsilon > 0$, $(1 - \epsilon) (1 - \sum_{k=1}^3 \frac{c_0^{r/2}}{\pi_k^0}) > \frac{1}{2}$, we hope to identify $K_0$ via the posterior mode with a reasonable estimate. To guarantee $K = K_0$ consistently using the posterior mode safety, we have to choose $c < c_0$, with $c_0$ satisfying

$$(1 - \epsilon) \left( 1 - \sum_{k=1}^{3} \frac{c_0^{r/2}}{\pi_k^0} \right) > \frac{1}{2} \Leftrightarrow \frac{1 - 2\epsilon}{2(1 - \epsilon)} \left( \sum_{k=1}^{3} \frac{1}{\pi_k^0} \right)^{-1} > c_0^{r/2} > c^{r/2}.$$
<figure><img src="image_4.png" /><figcaption>(a) Typical realisations without labels.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) True conditional expectations with labels.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) True and estimated inverse means.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) True and estimated forward means.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) Contour of estimated inverse cPDF.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) Contour of estimated forward cPDF.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(g) Contour of true inverse cPDF.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(h) Contour of true forward cPDF.</figcaption></figure>

**Figure 2.** Top row: Typical 1,000 realisations of GLLiM's true inverse cPDF with its true conditional expectations. 2^nd^ row and below: True and estimated inverse and forward cPDF of GLLiM with the true number of components ($K_{MTM} = 3$) using MTM algorithm for post-processing in MTM-BNP-GLLIM.
<figure><img src="image_3.png" /><figcaption>(a) Typical realisations without labels.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) True conditional expectations with labels.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) True and estimated inverse means.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) True and estimated forward means.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) Contour of estimated inverse cPDF.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) Contour of estimated forward cPDF.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(g) Contour of true inverse cPDF.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(h) Contour of true forward cPDF.</figcaption></figure>

**Figure 2.** Top row: Typical 1,000 realisations of GLLiM's true inverse cPDF with its true conditional expectations. 2^nd^ row and below: True and estimated inverse and forward cPDF of GLLiM with the true number of components ($K_{MTM} = 3$) using MTM algorithm for post-processing in MTM-BNP-GLLIM.
<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

Therefore, we can choose

$$ \left[ \frac{1 - 2\epsilon}{4(1 - \epsilon)} \left( \sum_{k=1}^{3} \frac{1}{\pi_k^0} \right)^{-1} \right]^{2/r} = c_0 > c, \text{ for all } \frac{1}{2} > \epsilon > 0. $$

In particular, it is unrealistic to obtain the exact computation of the upper bound $c_0$ and the lower bound $\frac{\tilde{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).

## 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
<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

Therefore, we can choose

$$ \left[ \frac{1 - 2\epsilon}{4(1 - \epsilon)} \left( \sum_{k=1}^{3} \frac{1}{\pi_k^0} \right)^{-1} \right]^{2/r} = c_0 > c, \text{ for all } \frac{1}{2} > \epsilon > 0. $$

In particular, it is unrealistic to obtain the exact computation of the upper bound $c_0$ and the lower bound $\frac{\tilde{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).

## 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
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Histogram of $K_{MTM}$ using variational posterior sample with 100 trials and $c = 0.45$.</figcaption></figure>

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 $D \gg N$.

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.

## References

* Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19(6):716-723. 2

* Alamichel, L., Bystrova, D., Arbel, J., and Kon Kam King, G. (2024). Bayesian mixture models (in)consistency for the number of clusters. *Scandinavian Journal of Statistics*. 3

* Arbel, J., Kon Kam King, G., Lijoi, A., Nieto-Barajas, L., and Prünster, I. (2021). BNPdensity: Bayesian nonparametric mixture modelling in R. *Australian & New Zealand Journal of Statistics*, 63(3):542-564. 3

* Arlot, S. (2019). Minimal penalties and the slope heuristics: a survey. *Journal de la Société Française de Statistique*, 160(3):1-106. 2

* Ascolani, F., Lijoi, A., Rebaudo, G., and Zanella, G. (2022). Clustering consistency with Dirichlet process mixtures. *Biometrika*, 110(2):551-558. 3

* Beal, M. J. (2003). *Variational algorithms for approximate Bayesian inference*. PhD thesis, University of London, University College London (United Kingdom). 30
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Histogram of $K_{MTM}$ using variational posterior sample with 100 trials and $c = 0.45$.</figcaption></figure>

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 $D \gg N$.

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.

## References

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19(6):716-723. 2

Alamichel, L., Bystrova, D., Arbel, J., and Kon Kam King, G. (2024). Bayesian mixture models (in)consistency for the number of clusters. *Scandinavian Journal of Statistics*. 3

Arbel, J., Kon Kam King, G., Lijoi, A., Nieto-Barajas, L., and Prünster, I. (2021). BNPdensity: Bayesian nonparametric mixture modelling in R. *Australian & New Zealand Journal of Statistics*, 63(3):542-564. 3

Arlot, S. (2019). Minimal penalties and the slope heuristics: a survey. *Journal de la Société Française de Statistique*, 160(3):1-106. 2

Ascolani, F., Lijoi, A., Rebaudo, G., and Zanella, G. (2022). Clustering consistency with Dirichlet process mixtures. *Biometrika*, 110(2):551-558. 3

Beal, M. J. (2003). *Variational algorithms for approximate Bayesian inference*. PhD thesis, University of London, University College London (United Kingdom). 30
Beal, M. J. and Ghahramani, Z. (2003). The variational Bayesian EM algorithm for incomplete data: with application to scoring graphical model structures. *Bayesian statistics*, 7(453-464):210. 30

Biernacki, C., Celeux, G., and Govaert, G. (2000). Assessing a mixture model for clustering with the integrated completed likelihood. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 22(7):719-725. 2

Birgé, L. and Massart, P. (2007). Minimal penalties for Gaussian model selection. *Probability Theory and Related Fields*, 138(1):33-73. 2, 20

Bishop, C. M. (2006). *Pattern Recognition and Machine Learning (Information Science and Statistics)*. Springer-Verlag, Berlin, Heidelberg. 30, 48, 51, 55, 63, 67

Blei, D. M. and Jordan, M. I. (2006). Variational inference for Dirichlet process mixtures. *Bayesian analysis*, 1(1):121-143. 8, 28, 47, 49, 61

Blein-Nicolas, M., Devijver, E., Gallopin, M., and Perthame, E. (2024). Nonlinear network-based quantitative trait prediction from biological data. *Journal of the Royal Statistical Society Series C: Applied Statistics*, page qlae012. 2, 6, 19

Boux, F., Forbes, F., Arbel, J., Lemasson, B., and Barbier, E. L. (2021). Bayesian inverse regression for vascular magnetic resonance fingerprinting. *IEEE Transactions on Medical Imaging*, 40(7):1827-1837. 1

Celeux, G., Frühwirth-Schnatter, S., and Robert, C. P. (2019). Model selection for mixture models–perspectives and strategies. In *Handbook of mixture analysis*, pages 117–154. Chapman and Hall/CRC. 2

Chaari, L., Vincent, T., Forbes, F., Dojat, M., and Ciuciu, P. (2013). Fast joint detection-estimation of evoked brain activity in event-related fMRI using a variational approach. *IEEE transactions on Medical Imaging*, 32(5):821-837. 30

Chamroukhi, F., Lecocq, F., and Nguyen, H. D. (2019). Regularized Estimation and Feature Selection in Mixtures of Gaussian-Gated Experts Models. In Nguyen, H., editor, *Statistics and Data Science*, pages 42–56, Singapore. Springer Singapore. 19

Chamroukhi, F., Samé, A., Govaert, G., and Aknin, P. (2010). A hidden process regression model for functional data description. Application to curve discrimination. *Neurocomputing*, 73(7-9):1210-1221. 16

Chen, Z., Deng, Y., Wu, Y., Gu, Q., and Li, Y. (2022). Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K., editors, *Advances in Neural Information Processing Systems*. 2

De Blasi, P., Favaro, S., Lijoi, A., Mena, R. H., Prünster, I., and Ruggiero, M. (2015). Are Gibbs-type priors the most natural generalization of the Dirichlet process? *IEEE transactions on pattern analysis and machine intelligence*, 37(2):212-229. 28

Deleforge, A., Forbes, F., and Horaud, R. (2015). High-dimensional regression with gaussian mixtures and partially-latent response variables. *Statistics and Computing*, 25(5):893-911. 1, 4, 5, 6, 43

Delon, J. and Desolneux, A. (2020). A Wasserstein-type distance in the space of Gaussian mixture models. *SIAM Journal on Imaging Sciences*, 13(2):936-970. 13

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data Via the EM Algorithm. *Journal of the Royal Statistical Society: Series B (Methodological)*, 39(1):1-22. 29

Do, T. G., Le, H. K., Nguyen, T., Pham, Q., Nguyen, B. T., Doan, T.-N., Liu, C., Ramasamy, S., Li, X., and HOI, S. (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. Association for Computational Linguistics. 2

Durand, J.-B., Forbes, F., Phan, C., Truong, L., Nguyen, H., and Dama, F. (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):171-204. 3, 10, 29

Escobar, M. D. and West, M. (1995). Bayesian Density Estimation and Inference Using Mixtures. *Journal of the American Statistical Association*, 90(430):577-588. 3

Favaro, S., Lijoi, A., Nava, C., Nipoti, B., Pruenster, I., and Teh, Y. W. (2016). On the
Beal, M. J. and Ghahramani, Z. (2003). The variational Bayesian EM algorithm for incomplete data: with application to scoring graphical model structures. *Bayesian statistics*, 7(453-464):210. 30

Biernacki, C., Celeux, G., and Govaert, G. (2000). Assessing a mixture model for clustering with the integrated completed likelihood. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 22(7):719-725. 2

Birgé, L. and Massart, P. (2007). Minimal penalties for Gaussian model selection. *Probability Theory and Related Fields*, 138(1):33-73. 2, 20

Bishop, C. M. (2006). *Pattern Recognition and Machine Learning (Information Science and Statistics)*. Springer-Verlag, Berlin, Heidelberg. 30, 48, 51, 55, 63, 67

Blei, D. M. and Jordan, M. I. (2006). Variational inference for Dirichlet process mixtures. *Bayesian analysis*, 1(1):121-143. 8, 28, 47, 49, 61

Blein-Nicolas, M., Devijver, E., Gallopin, M., and Perthame, E. (2024). Nonlinear network-based quantitative trait prediction from biological data. *Journal of the Royal Statistical Society Series C: Applied Statistics*, page qlae012. 2, 6, 19

Boux, F., Forbes, F., Arbel, J., Lemasson, B., and Barbier, E. L. (2021). Bayesian inverse regression for vascular magnetic resonance fingerprinting. *IEEE Transactions on Medical Imaging*, 40(7):1827-1837. 1

Celeux, G., Frühwirth-Schnatter, S., and Robert, C. P. (2019). Model selection for mixture models–perspectives and strategies. In *Handbook of mixture analysis*, pages 117–154. Chapman and Hall/CRC. 2

Chaari, L., Vincent, T., Forbes, F., Dojat, M., and Ciuciu, P. (2013). Fast joint detection-estimation of evoked brain activity in event-related fMRI using a variational approach. *IEEE transactions on Medical Imaging*, 32(5):821-837. 30

Chamroukhi, F., Lecocq, F., and Nguyen, H. D. (2019). Regularized Estimation and Feature Selection in Mixtures of Gaussian-Gated Experts Models. In Nguyen, H., editor, *Statistics and Data Science*, pages 42–56, Singapore. Springer Singapore. 19

Chamroukhi, F., Samé, A., Govaert, G., and Aknin, P. (2010). A hidden process regression model for functional data description. Application to curve discrimination. *Neurocomputing*, 73(7-9):1210-1221. 16

Chen, Z., Deng, Y., Wu, Y., Gu, Q., and Li, Y. (2022). Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K., editors, *Advances in Neural Information Processing Systems*. 2

De Blasi, P., Favaro, S., Lijoi, A., Mena, R. H., Prünster, I., and Ruggiero, M. (2015). Are Gibbs-type priors the most natural generalization of the Dirichlet process? *IEEE transactions on pattern analysis and machine intelligence*, 37(2):212-229. 28

Deleforge, A., Forbes, F., and Horaud, R. (2015). High-dimensional regression with gaussian mixtures and partially-latent response variables. *Statistics and Computing*, 25(5):893-911. 1, 4, 5, 6, 43

Delon, J. and Desolneux, A. (2020). A Wasserstein-type distance in the space of Gaussian mixture models. *SIAM Journal on Imaging Sciences*, 13(2):936-970. 13

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data Via the EM Algorithm. *Journal of the Royal Statistical Society: Series B (Methodological)*, 39(1):1-22. 29

Do, T. G., Le, H. K., Nguyen, T., Pham, Q., Nguyen, B. T., Doan, T.-N., Liu, C., Ramasamy, S., Li, X., and HOI, S. (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. Association for Computational Linguistics. 2

Durand, J.-B., Forbes, F., Phan, C., Truong, L., Nguyen, H., and Dama, F. (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):171-204. 3, 10, 29

Escobar, M. D. and West, M. (1995). Bayesian Density Estimation and Inference Using Mixtures. *Journal of the American Statistical Association*, 90(430):577-588. 3

Favaro, S., Lijoi, A., Nava, C., Nipoti, B., Pruenster, I., and Teh, Y. W. (2016). On the
stick-breaking representation for homogeneous NRMIs. *Bayesian Analysis*, 11(3):697-724.
28

Ferguson, T. S. (1973). A Bayesian Analysis of Some Nonparametric Problems. *The Annals of Statistics*, 1(2):209 – 230. 27

Forbes, F., Arnaud, A., Lemasson, B., and Barbier, E. (2019). Component elimination strategies to fit mixtures of multiple scale distributions. In *RSSDS 2019 - Research School on Statistics and Data Science*, volume 1150 of *Communications in Computer and Information Science*, pages 81–95, Melbourne, Australia. Springer. 3

Forbes, F., Nguyen, H. D., Nguyen, T., and Arbel, J. (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. 2

Forbes, F., Nguyen, H. D., Nguyen, T., and Arbel, J. (2022b). Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. *Statistics and Computing*, 32(5):85. 2

Gao, F. and Vaart, A. v. d. (2016). Posterior contraction rates for deconvolution of Dirichlet-Laplace mixtures. *Electronic Journal of Statistics*, 10(1):608 – 627. 14

Genovese, C. R. and Wasserman, L. (2000). Rates of convergence for the Gaussian mixture sieve. *The Annals of Statistics*, 28(4):1105 – 1127. 2

Ghosal, S., Ghosh, J. K., and Vaart, A. W. v. d. (2000). Convergence rates of posterior distributions. *The Annals of Statistics*, 28(2):500 – 531. 20

Ghosal, S. and Van der Vaart, A. (2017). *Fundamentals of nonparametric Bayesian inference*, volume 44. Cambridge University Press. 3, 28, 49

Green, P. J. and Richardson, S. (2001). Modelling Heterogeneity With and Without the Dirichlet Process. *Scandinavian Journal of Statistics*, 28(2):355–375. 3

Guha, A., Ho, N., and Nguyen, X. (2021). On posterior contraction of parameters and interpretability in Bayesian mixture modeling. *Bernoulli*, 27(4):2159 – 2188. 3, 14, 15, 16, 17, 19

Hjort, N. L., Holmes, C., Müller, P., and Walker, S. G., editors (2010). *Bayesian Nonparametries*. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press. 3

Ho, N. and Nguyen, X. (2016a). Convergence rates of parameter estimation for some weakly identifiable finite mixtures. *The Annals of Statistics*, 44(6):2726 – 2755. 2

Ho, N. and Nguyen, X. (2016b). On strong identifiability and convergence rates of parameter estimation in finite mixtures. *Electronic Journal of Statistics*, 10(1):271–307. 2, 13, 14

Ho, N., Yang, C.-Y., and Jordan, M. I. (2022). Convergence Rates for Gaussian Mixtures of Experts. *Journal of Machine Learning Research*, 23(323):1–81. 2

Hoadley, B. (1970). A Bayesian Look at Inverse Linear Regression. *Journal of the American Statistical Association*, 65(329):356–369. 1

Ingrassia, S., Minotti, S. C., and Vittadini, G. (2012). Local Statistical Modeling via a Cluster-Weighted Approach with Elliptical Distributions. *Journal of Classification*, 29(3):363–401. 6

Ishwaran, H. and James, L. F. (2001). Gibbs Sampling Methods for Stick-Breaking Priors. *Journal of the American Statistical Association*, 96(453):161–173. 28

Jacobs, R. A., Jordan, M. I., Nowlan, S. J., and Hinton, G. E. (1991). Adaptive mixtures of local experts. *Neural computation*, 3(1):79–87. 2

Jiang, W. and Tanner, M. A. (1999). Hierarchical mixtures-of-experts for exponential family regression models: approximation and maximum likelihood estimation. *Annals of Statistics*, pages 987–1011. 2

Kock, L., Klein, N., and Nott, D. J. (2022). Variational inference and sparsity in high-dimensional deep Gaussian mixture models. *Statistics and Computing*, 32(5):70. 3

Kruijer, W., Rousseau, J., and Vaart, A. v. d. (2010). Adaptive Bayesian density estimation with location-scale mixtures. *Electronic Journal of Statistics*, 4(none):1225 – 1257. 20

Kugler, B., Forbes, F., and Douté, S. (2022). Fast Bayesian Inversion for high dimensional inverse problems. *Statistics and Computing*, 32(2). 2
stick-breaking representation for homogeneous NRMIs. *Bayesian Analysis*, 11(3):697-724.
28

Ferguson, T. S. (1973). A Bayesian Analysis of Some Nonparametric Problems. *The Annals of Statistics*, 1(2):209 – 230. 27

Forbes, F., Arnaud, A., Lemasson, B., and Barbier, E. (2019). Component elimination strategies to fit mixtures of multiple scale distributions. In *RSSDS 2019 - Research School on Statistics and Data Science*, volume 1150 of *Communications in Computer and Information Science*, pages 81–95, Melbourne, Australia. Springer. 3

Forbes, F., Nguyen, H. D., Nguyen, T., and Arbel, J. (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. 2

Forbes, F., Nguyen, H. D., Nguyen, T., and Arbel, J. (2022b). Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. *Statistics and Computing*, 32(5):85. 2

Gao, F. and Vaart, A. v. d. (2016). Posterior contraction rates for deconvolution of Dirichlet-Laplace mixtures. *Electronic Journal of Statistics*, 10(1):608 – 627. 14

Genovese, C. R. and Wasserman, L. (2000). Rates of convergence for the Gaussian mixture sieve. *The Annals of Statistics*, 28(4):1105 – 1127. 2

Ghosal, S., Ghosh, J. K., and Vaart, A. W. v. d. (2000). Convergence rates of posterior distributions. *The Annals of Statistics*, 28(2):500 – 531. 20

Ghosal, S. and Van der Vaart, A. (2017). *Fundamentals of nonparametric Bayesian inference*, volume 44. Cambridge University Press. 3, 28, 49

Green, P. J. and Richardson, S. (2001). Modelling Heterogeneity With and Without the Dirichlet Process. *Scandinavian Journal of Statistics*, 28(2):355–375. 3

Guha, A., Ho, N., and Nguyen, X. (2021). On posterior contraction of parameters and interpretability in Bayesian mixture modeling. *Bernoulli*, 27(4):2159 – 2188. 3, 14, 15, 16, 17, 19

Hjort, N. L., Holmes, C., Müller, P., and Walker, S. G., editors (2010). *Bayesian Nonparametrics*. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press. 3

Ho, N. and Nguyen, X. (2016a). Convergence rates of parameter estimation for some weakly identifiable finite mixtures. *The Annals of Statistics*, 44(6):2726 – 2755. 2

Ho, N. and Nguyen, X. (2016b). On strong identifiability and convergence rates of parameter estimation in finite mixtures. *Electronic Journal of Statistics*, 10(1):271–307. 2, 13, 14

Ho, N., Yang, C.-Y., and Jordan, M. I. (2022). Convergence Rates for Gaussian Mixtures of Experts. *Journal of Machine Learning Research*, 23(323):1–81. 2

Hoadley, B. (1970). A Bayesian Look at Inverse Linear Regression. *Journal of the American Statistical Association*, 65(329):356–369. 1

Ingrassia, S., Minotti, S. C., and Vittadini, G. (2012). Local Statistical Modeling via a Cluster-Weighted Approach with Elliptical Distributions. *Journal of Classification*, 29(3):363–401. 6

Ishwaran, H. and James, L. F. (2001). Gibbs Sampling Methods for Stick-Breaking Priors. *Journal of the American Statistical Association*, 96(453):161–173. 28

Jacobs, R. A., Jordan, M. I., Nowlan, S. J., and Hinton, G. E. (1991). Adaptive mixtures of local experts. *Neural computation*, 3(1):79–87. 2

Jiang, W. and Tanner, M. A. (1999). Hierarchical mixtures-of-experts for exponential family regression models: approximation and maximum likelihood estimation. *Annals of Statistics*, pages 987–1011. 2

Kock, L., Klein, N., and Nott, D. J. (2022). Variational inference and sparsity in high-dimensional deep Gaussian mixture models. *Statistics and Computing*, 32(5):70. 3

Kruijer, W., Rousseau, J., and Vaart, A. v. d. (2010). Adaptive Bayesian density estimation with location-scale mixtures. *Electronic Journal of Statistics*, 4(none):1225 – 1257. 20

Kugler, B., Forbes, F., and Douté, S. (2022). Fast Bayesian Inversion for high dimensional inverse problems. *Statistics and Computing*, 32(2). 2
Lathuilière, S., Juge, R., Mesejo, P., Muñoz-Salinas, R., and Horaud, R. (2017). Deep mixture of linear inverse regressions applied to head-pose estimation. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition*, pages 4817–4825. 2

Li, K.-C. (1991). Sliced Inverse Regression for Dimension Reduction. *Journal of the American Statistical Association*, 86(414):316–327. 1, 4

Li, N., Li, W., Jiang, Y., and Xia, S.-T. (2022). Deep Dirichlet process mixture models. In Cussens, J. and Zhang, K., editors, *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, volume 180 of *Proceedings of Machine Learning Research*, pages 1138–1147. PMLR. 3

Luo, C. and Sun, S. (2017). Variational Mixtures of Gaussian Processes for Classification. In *Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17*, pages 4603–4609. 11

Lü, H., Arbel, J., and Forbes, F. (2020). Bayesian nonparametric priors for hidden Markov random fields. *Statistics and Computing*, 30(4):1015–1035. 9, 28, 29

Maceachern, S. N. and Müller, P. (1998). Estimating Mixture of Dirichlet Process Models. *Journal of Computational and Graphical Statistics*, 7(2):223–238. 3

McLachlan, G. J. and Krishnan, T. (1997). *The EM algorithm and extensions*. Wiley, New York, USA. 29

McLachlan, G. J. and Peel, D. (2000). *Finite Mixture Models*. John Wiley & Sons. 2

Miller, J. W. and Harrison, M. T. (2014). Inconsistency of Pitman-Yor Process Mixtures for the Number of Components. *Journal of Machine Learning Research*, 15(96):3333–3370. 3

Montuelle, L. and Le Pennec, E. (2014). Mixture of Gaussian regressions model with logistic weights, a penalized maximum likelihood approach. *Electronic Journal of Statistics*, 8(1):1661–1695. 16

Neal, R. M. (2000). Markov Chain Sampling Methods for Dirichlet Process Mixture Models. *Journal of Computational and Graphical Statistics*, 9(2):249–265. 3

Neal, R. M. and Hinton, G. E. (1998). A View of the EM Algorithm that Justifies Incremental, Sparse, and other Variants. In Jordan, M. I., editor, *Learning in Graphical Models*, pages 355–368. Springer Netherlands, Dordrecht. 29

Nguyen, D. T., Jacquemoud, S., Lucas, A., Douté, S., Ferrari, C., Coustance, S., Marcq, S., and Meygret, A. (2024). Unveiling the Characteristics of the Lunar Surface by Massive Inversion of a Photometric Model. In *LPI Contributions*, volume 3040 of *LPI Contributions*, page 1998. 2

Nguyen, H., Akbarian, P., Nguyen, T., and Ho, N. (2024a). A General Theory for Softmax Gating Multinomial Logistic Mixture of Experts. In *Proceedings of The 41st International Conference on Machine Learning*. 2

Nguyen, H., Nguyen, T., and Ho, N. (2023a). Demystifying Softmax Gating in Gaussian Mixture of Experts. In *Advances in Neural Information Processing Systems*. 2

Nguyen, H., Nguyen, T., Nguyen, K., and Ho, N. (2024b). Towards Convergence Rates for Parameter Estimation in Gaussian-gated Mixture of Experts. In Dasgupta, S., Mandt, S., and Li, Y., editors, *Proceedings of The 27th International Conference on Artificial Intelligence and Statistics*, volume 238 of *Proceedings of Machine Learning Research*, pages 2683–2691. PMLR. 2

Nguyen, H. D. and Chamroukhi, F. (2018). Practical and theoretical aspects of mixture-of-experts modeling: An overview. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 8(4):e1246. 2

Nguyen, H. D., Chamroukhi, F., and Forbes, F. (2019). Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. *Neurocomputing*, 366:208–214. 2

Nguyen, H. D., Fryer, D., and McLachlan, G. J. (2022a). Order selection with confidence for finite mixture models. *Journal of the Korean Statistical Society*. 3

Nguyen, H. D., Lloyd-Jones, L. R., and McLachlan, G. J. (2016). A universal approximation theorem for mixture-of-experts models. *Neural computation*, 28(12):2585–2593. 2

Nguyen, H. D., Nguyen, T., Chamroukhi, F., and McLachlan, G. J. (2021a). Approximations of
Lathuilière, S., Juge, R., Mesejo, P., Muñoz-Salinas, R., and Horaud, R. (2017). Deep mixture of linear inverse regressions applied to head-pose estimation. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition*, pages 4817–4825. 2

Li, K.-C. (1991). Sliced Inverse Regression for Dimension Reduction. *Journal of the American Statistical Association*, 86(414):316–327. 1, 4

Li, N., Li, W., Jiang, Y., and Xia, S.-T. (2022). Deep Dirichlet process mixture models. In Cussens, J. and Zhang, K., editors, *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, volume 180 of *Proceedings of Machine Learning Research*, pages 1138–1147. PMLR. 3

Luo, C. and Sun, S. (2017). Variational Mixtures of Gaussian Processes for Classification. In *Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17*, pages 4603–4609. 11

Lü, H., Arbel, J., and Forbes, F. (2020). Bayesian nonparametric priors for hidden Markov random fields. *Statistics and Computing*, 30(4):1015–1035. 9, 28, 29

Maceachern, S. N. and Müller, P. (1998). Estimating Mixture of Dirichlet Process Models. *Journal of Computational and Graphical Statistics*, 7(2):223–238. 3

McLachlan, G. J. and Krishnan, T. (1997). *The EM algorithm and extensions*. Wiley, New York, USA. 29

McLachlan, G. J. and Peel, D. (2000). *Finite Mixture Models*. John Wiley & Sons. 2

Miller, J. W. and Harrison, M. T. (2014). Inconsistency of Pitman-Yor Process Mixtures for the Number of Components. *Journal of Machine Learning Research*, 15(96):3333–3370. 3

Montuelle, L. and Le Pennec, E. (2014). Mixture of Gaussian regressions model with logistic weights, a penalized maximum likelihood approach. *Electronic Journal of Statistics*, 8(1):1661–1695. 16

Neal, R. M. (2000). Markov Chain Sampling Methods for Dirichlet Process Mixture Models. *Journal of Computational and Graphical Statistics*, 9(2):249–265. 3

Neal, R. M. and Hinton, G. E. (1998). A View of the EM Algorithm that Justifies Incremental, Sparse, and other Variants. In Jordan, M. I., editor, *Learning in Graphical Models*, pages 355–368. Springer Netherlands, Dordrecht. 29

Nguyen, D. T., Jacquemoud, S., Lucas, A., Douté, S., Ferrari, C., Coustance, S., Marcq, S., and Meygret, A. (2024). Unveiling the Characteristics of the Lunar Surface by Massive Inversion of a Photometric Model. In *LPI Contributions*, volume 3040 of *LPI Contributions*, page 1998. 2

Nguyen, H., Akbarian, P., Nguyen, T., and Ho, N. (2024a). A General Theory for Softmax Gating Multinomial Logistic Mixture of Experts. In *Proceedings of The 41st International Conference on Machine Learning*. 2

Nguyen, H., Nguyen, T., and Ho, N. (2023a). Demystifying Softmax Gating in Gaussian Mixture of Experts. In *Advances in Neural Information Processing Systems*. 2

Nguyen, H., Nguyen, T., Nguyen, K., and Ho, N. (2024b). Towards Convergence Rates for Parameter Estimation in Gaussian-gated Mixture of Experts. In Dasgupta, S., Mandt, S., and Li, Y., editors, *Proceedings of The 27th International Conference on Artificial Intelligence and Statistics*, volume 238 of *Proceedings of Machine Learning Research*, pages 2683–2691. PMLR. 2

Nguyen, H. D. and Chamroukhi, F. (2018). Practical and theoretical aspects of mixture-of-experts modeling: An overview. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 8(4):e1246. 2

Nguyen, H. D., Chamroukhi, F., and Forbes, F. (2019). Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. *Neurocomputing*, 366:208–214. 2

Nguyen, H. D., Fryer, D., and McLachlan, G. J. (2022a). Order selection with confidence for finite mixture models. *Journal of the Korean Statistical Society*. 3

Nguyen, H. D., Lloyd-Jones, L. R., and McLachlan, G. J. (2016). A universal approximation theorem for mixture-of-experts models. *Neural computation*, 28(12):2585–2593. 2

Nguyen, H. D., Nguyen, T., Chamroukhi, F., and McLachlan, G. J. (2021a). Approximations of
conditional probability density functions in Lebesgue spaces via mixture of experts models. *Journal of Statistical Distributions and Applications*, 8(1):13. 2

Nguyen, T. (2021). *Model Selection and Approximation in High-dimensional Mixtures of Experts Models: from Theory to Practice*. PhD Thesis, Normandie Université. 2

Nguyen, T. and Bonilla, E. (2014). Fast Allocation of Gaussian Process Experts. In Xing, E. P. and Jebara, T., editors, *Proceedings of the 31st International Conference on Machine Learning*, volume 32 of *Proceedings of Machine Learning Research*, pages 145–153, Beijing, China. PMLR. 11

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and Forbes, F. (2021b). Non-asymptotic model selection in block-diagonal mixture of polynomial experts models. *Preprint*. arXiv:2104.08959. 2

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and Forbes, F. (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. 2

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and McLachlan, G. J. (2023b). Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. *Communications in Statistics - Theory and Methods*, 52(14):5048–5059. 2

Nguyen, T., Nguyen, D. N., Nguyen, H. D., and Chamroukhi, F. (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*. 2

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and Forbes, F. (2022c). A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models. *Electronic Journal of Statistics*, 16(2):4742 – 4822. 2, 6, 16

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and McLachlan, G. J. (2020). Approximation by finite mixtures of continuous density functions that vanish at infinity. *Cogent Mathematics & Statistics*, 7(1):1750861. 2

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and McLachlan, G. J. (2023d). Non-asymptotic oracle inequalities for the Lasso in high-dimensional mixture of experts. *arXiv:2009.10622*. 2

Nguyen, X. (2013). Convergence of latent mixing measures in finite and infinite mixture models. *The Annals of Statistics*, 41(1):370–400. 2, 13, 14

Norets, A. (2010). Approximation of conditional densities by smooth mixtures of regressions. *The Annals of Statistics*, 38(3):1733 – 1766. 2

Norets, A. and Pati, D. (2017). Adaptive Bayesian estimation of conditional densities. *Econometric Theory*, 33(4):980–1012. 20

Perthame, E., Forbes, F., and Deleforge, A. (2018). Inverse regression approach to robust nonlinear high-to-low dimensional mapping. *Journal of Multivariate Analysis*, 163(C):1–14. 6

Pitman, J. and Yor, M. (1997). The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator. *The Annals of Probability*, 25(2):855–900. 7, 28

Rakhlin, A., Panchenko, D., and Mukherjee, S. (2005). Risk bounds for mixture density estimation. *ESAIM: PS*, 9:220–229. 2

Rousseau, J. and Mengersen, K. (2011). Asymptotic behaviour of the posterior distribution in overfitted mixture models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 73(5):689–710. 3

Ruiz, F. and Titsias, M. (2019). A Contrastive Divergence for Combining Variational Inference and MCMC. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 5537–5545. PMLR. 20

Schwarz, G. (1978). Estimating the dimension of a model. *The Annals of Statistics*, 6(2):461–464. 2

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. *Statistica Sinica*, 4:639–650. 28
conditional probability density functions in Lebesgue spaces via mixture of experts models. *Journal of Statistical Distributions and Applications*, 8(1):13. 2

Nguyen, T. (2021). *Model Selection and Approximation in High-dimensional Mixtures of Experts Models: from Theory to Practice*. PhD Thesis, Normandie Université. 2

Nguyen, T. and Bonilla, E. (2014). Fast Allocation of Gaussian Process Experts. In Xing, E. P. and Jebara, T., editors, *Proceedings of the 31st International Conference on Machine Learning*, volume 32 of *Proceedings of Machine Learning Research*, pages 145–153, Beijing, China. PMLR. 11

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and Forbes, F. (2021b). Non-asymptotic model selection in block-diagonal mixture of polynomial experts models. *Preprint*. arXiv:2104.08959. 2

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and Forbes, F. (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. 2

Nguyen, T., Chamroukhi, F., Nguyen, H. D., and McLachlan, G. J. (2023b). Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. *Communications in Statistics - Theory and Methods*, 52(14):5048–5059. 2

Nguyen, T., Nguyen, D. N., Nguyen, H. D., and Chamroukhi, F. (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*. 2

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and Forbes, F. (2022c). A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models. *Electronic Journal of Statistics*, 16(2):4742 – 4822. 2, 6, 16

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and McLachlan, G. J. (2020). Approximation by finite mixtures of continuous density functions that vanish at infinity. *Cogent Mathematics & Statistics*, 7(1):1750861. 2

Nguyen, T., Nguyen, H. D., Chamroukhi, F., and McLachlan, G. J. (2023d). Non-asymptotic oracle inequalities for the Lasso in high-dimensional mixture of experts. *arXiv:2009.10622*. 2

Nguyen, X. (2013). Convergence of latent mixing measures in finite and infinite mixture models. *The Annals of Statistics*, 41(1):370–400. 2, 13, 14

Norets, A. (2010). Approximation of conditional densities by smooth mixtures of regressions. *The Annals of Statistics*, 38(3):1733 – 1766. 2

Norets, A. and Pati, D. (2017). Adaptive Bayesian estimation of conditional densities. *Econometric Theory*, 33(4):980–1012. 20

Perthame, E., Forbes, F., and Deleforge, A. (2018). Inverse regression approach to robust nonlinear high-to-low dimensional mapping. *Journal of Multivariate Analysis*, 163(C):1–14. 6

Pitman, J. and Yor, M. (1997). The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator. *The Annals of Probability*, 25(2):855–900. 7, 28

Rakhlin, A., Panchenko, D., and Mukherjee, S. (2005). Risk bounds for mixture density estimation. *ESAIM: PS*, 9:220–229. 2

Rousseau, J. and Mengersen, K. (2011). Asymptotic behaviour of the posterior distribution in overfitted mixture models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 73(5):689–710. 3

Ruiz, F. and Titsias, M. (2019). A Contrastive Divergence for Combining Variational Inference and MCMC. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 5537–5545. PMLR. 20

Schwarz, G. (1978). Estimating the dimension of a model. *The Annals of Statistics*, 6(2):461–464. 2

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. *Statistica Sinica*, 4:639–650. 28
Shen, W., Tokdar, S. T., and Ghosal, S. (2013). Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. *Biometrika*, 100(3):623–640. 2, 20

Svensén, M. and Bishop, C. M. (2005). Robust Bayesian mixture modelling. *Neurocomputing*, 64:235–252. 10

Titterington, D. M. and Wang, B. (2006). Convergence properties of a general algorithm for calculating variational Bayesian estimates for a normal mixture model. *Bayesian Analysis*, 1(3):625–650. 20

Villani, C. (2003). *Topics in optimal transportation*, volume 58 of *Graduate Studies in Mathematics*. American Mathematical Society. 13

Villani, C. (2009). *Optimal transport: old and new*, volume 338. Springer. 13

Wang, C., Paisley, J., and Blei, D. M. (2011). Online Variational Inference for the Hierarchical Dirichlet Process. In Gordon, G., Dunson, D., and Dudík, M., editors, *Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics*, volume 15 of *Proceedings of Machine Learning Research*, pages 752–760, Fort Lauderdale, FL, USA. PMLR. 8

Westerhout, J., Nguyen, T., Guo, X., and Nguyen, H. D. (2024). On the Asymptotic Distribution of the Minimum Empirical Risk. In *Forty-first International Conference on Machine Learning*. 2

Wu, D. and Ma, J. (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):4894–4904. 11

Xu, L., Jordan, M., and Hinton, G. E. (1995). An Alternative Model for Mixtures of Experts. In Tesauro, G., Touretzky, D., and Leen, T., editors, *Advances in Neural Information Processing Systems*, volume 7. MIT Press. 1

Yuan, C. and Neubauer, C. (2008). Variational Mixture of Gaussian Process Experts. In Koller, D., Schuurmans, D., Bengio, Y., and Bottou, L., editors, *Advances in Neural Information Processing Systems*, volume 21. Curran Associates, Inc. 10

Yuksel, S. E., Wilson, J. N., and Gader, P. D. (2012). Twenty Years of Mixture of Experts. *IEEE Transactions on Neural Networks and Learning Systems*, 23(8):1177–1193. 2

Zhang, F. and Gao, C. (2020). Convergence rates of variational posterior distributions. *The Annals of Statistics*, 48(4):2180–2207. 20
Shen, W., Tokdar, S. T., and Ghosal, S. (2013). Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. *Biometrika*, 100(3):623–640. 2, 20

Svensén, M. and Bishop, C. M. (2005). Robust Bayesian mixture modelling. *Neurocomputing*, 64:235–252. 10

Titterington, D. M. and Wang, B. (2006). Convergence properties of a general algorithm for calculating variational Bayesian estimates for a normal mixture model. *Bayesian Analysis*, 1(3):625–650. 20

Villani, C. (2003). *Topics in optimal transportation*, volume 58 of *Graduate Studies in Mathematics*. American Mathematical Society. 13

Villani, C. (2009). *Optimal transport: old and new*, volume 338. Springer. 13

Wang, C., Paisley, J., and Blei, D. M. (2011). Online Variational Inference for the Hierarchical Dirichlet Process. In Gordon, G., Dunson, D., and Dudík, M., editors, *Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics*, volume 15 of *Proceedings of Machine Learning Research*, pages 752–760, Fort Lauderdale, FL, USA. PMLR. 8

Westerhout, J., Nguyen, T., Guo, X., and Nguyen, H. D. (2024). On the Asymptotic Distribution of the Minimum Empirical Risk. In *Forty-first International Conference on Machine Learning*. 2

Wu, D. and Ma, J. (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):4894–4904. 11

Xu, L., Jordan, M., and Hinton, G. E. (1995). An Alternative Model for Mixtures of Experts. In Tesauro, G., Touretzky, D., and Leen, T., editors, *Advances in Neural Information Processing Systems*, volume 7. MIT Press. 1

Yuan, C. and Neubauer, C. (2008). Variational Mixture of Gaussian Process Experts. In Koller, D., Schuurmans, D., Bengio, Y., and Bottou, L., editors, *Advances in Neural Information Processing Systems*, volume 21. Curran Associates, Inc. 10

Yuksel, S. E., Wilson, J. N., and Gader, P. D. (2012). Twenty Years of Mixture of Experts. *IEEE Transactions on Neural Networks and Learning Systems*, 23(8):1177–1193. 2

Zhang, F. and Gao, C. (2020). Convergence rates of variational posterior distributions. *The Annals of Statistics*, 48(4):2180–2207. 20
# 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 to E. Finally, Appendix F proposes a more general model with a hyper prior on the gating parameters, referred to as BNP-GLLIM2.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>BNP-GLLIM model</td><td>4</td></tr><tr><td>2.1</td><td>Inverse regression framework</td><td>4</td></tr><tr><td>2.2</td><td>Nonlinear high-to-low dimensional mapping via GLLIM model</td><td>5</td></tr><tr><td>2.3</td><td>Construction of BNP-GLLIM model</td><td>6</td></tr><tr><td>3</td><td>Variational inference for BNP-GLLIM</td><td>7</td></tr><tr><td>3.1</td><td>VB-E steps</td><td>8</td></tr><tr><td>3.2</td><td>VB-M steps</td><td>9</td></tr><tr><td>3.3</td><td>Evidence lower-bound (ELBO)</td><td>10</td></tr><tr><td>4</td><td>Predictive conditional density</td><td>11</td></tr><tr><td>4.1</td><td>Joint density</td><td>11</td></tr><tr><td>4.2</td><td>Inverse conditional density</td><td>11</td></tr><tr><td>4.3</td><td>Forward conditional density</td><td>12</td></tr><tr><td>5</td><td>Bayesian nonparametric model selection</td><td>12</td></tr><tr><td>5.1</td><td>Posterior contraction rate in Bayesian infinite mixtures</td><td>13</td></tr><tr><td>5.2</td><td>Merge-Truncate-Merge (MTM) algorithm for BNP-GLLIM</td><td>14</td></tr><tr><td>6</td><td>Numerical experiments</td><td>16</td></tr><tr><td>6.1</td><td>Data generation</td><td>16</td></tr><tr><td>6.2</td><td>Model selection, clustering and regression tasks via MTM-BNP-GLLIM</td><td>17</td></tr><tr><td>7</td><td>Perspectives</td><td>19</td></tr><tr><td>A</td><td>Bayesian nonparametric priors</td><td>27</td></tr><tr><td>B</td><td>Variational Bayesian expectation-maximisation principle</td><td>29</td></tr><tr><td>C</td><td>Details of VBEM for BNP-GLLIM model</td><td>31</td></tr><tr><td>C.1</td><td>VB-E-τ step from Section 3.1</td><td>31</td></tr><tr><td>C.2</td><td>VB-E-(α, σ) step from Section 3.1</td><td>32</td></tr><tr><td>C.3</td><td>VB-E-Z step from Section 3.1</td><td>34</td></tr><tr><td>C.3.1</td><td>Updating Σ<sub>k</sub></td><td>36</td></tr><tr><td>C.4</td><td>VB-M-(c, Γ) step from Section 3.2</td><td>36</td></tr></table>
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 to E. Finally, Appendix F proposes a more general model with a hyper prior on the gating parameters, referred to as BNP-GLLIM2.

Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>BNP-GLLIM model</td><td>4</td></tr><tr><td>2.1</td><td>Inverse regression framework</td><td>4</td></tr><tr><td>2.2</td><td>Nonlinear high-to-low dimensional mapping via GLLIM model</td><td>5</td></tr><tr><td>2.3</td><td>Construction of BNP-GLLIM model</td><td>6</td></tr><tr><td>3</td><td>Variational inference for BNP-GLLIM</td><td>7</td></tr><tr><td>3.1</td><td>VB-E steps</td><td>8</td></tr><tr><td>3.2</td><td>VB-M steps</td><td>9</td></tr><tr><td>3.3</td><td>Evidence lower-bound (ELBO)</td><td>10</td></tr><tr><td>4</td><td>Predictive conditional density</td><td>11</td></tr><tr><td>4.1</td><td>Joint density</td><td>11</td></tr><tr><td>4.2</td><td>Inverse conditional density</td><td>11</td></tr><tr><td>4.3</td><td>Forward conditional density</td><td>12</td></tr><tr><td>5</td><td>Bayesian nonparametric model selection</td><td>12</td></tr><tr><td>5.1</td><td>Posterior contraction rate in Bayesian infinite mixtures</td><td>13</td></tr><tr><td>5.2</td><td>Merge-Truncate-Merge (MTM) algorithm for BNP-GLLIM</td><td>14</td></tr><tr><td>6</td><td>Numerical experiments</td><td>16</td></tr><tr><td>6.1</td><td>Data generation</td><td>16</td></tr><tr><td>6.2</td><td>Model selection, clustering and regression tasks via MTM-BNP-GLLIM</td><td>17</td></tr><tr><td>7</td><td>Perspectives</td><td>19</td></tr><tr><td>A</td><td>Bayesian nonparametric priors</td><td>27</td></tr><tr><td>B</td><td>Variational Bayesian expectation-maximisation principle</td><td>29</td></tr><tr><td>C</td><td>Details of VBEM for BNP-GLLIM model</td><td>31</td></tr><tr><td>C.1</td><td>VB-E-τ step from Section 3.1</td><td>31</td></tr><tr><td>C.2</td><td>VB-E-(α, σ) step from Section 3.1</td><td>32</td></tr><tr><td>C.3</td><td>VB-E-Z step from Section 3.1</td><td>34</td></tr><tr><td>C.3.1</td><td>Updating Σ<sub>k</sub></td><td>36</td></tr><tr><td>C.4</td><td>VB-M-(c, Γ) step from Section 3.2</td><td>36</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>C.5</td>
      <td>VB-M-(<b>A</b>, <b>b</b>, <b>Σ</b>) step from Section 3.2</td>
      <td>37</td>
    </tr>
    <tr>
      <td>C.5.1</td>
      <td>Updating <b>b</b><sub>k</sub></td>
      <td>38</td>
    </tr>
    <tr>
      <td>C.5.2</td>
      <td>Updating <b>A</b><sub>k</sub></td>
      <td>38</td>
    </tr>
    <tr>
      <td>D</td>
      <td>Details on the ELBO</td>
      <td>41</td>
    </tr>
    <tr>
      <td>E</td>
      <td>Technical proofs</td>
      <td>42</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>Proof of Lemma 5.2</td>
      <td>42</td>
    </tr>
    <tr>
      <td>E.2</td>
      <td>Proof of Proposition 3.1</td>
      <td>43</td>
    </tr>
    <tr>
      <td>E.3</td>
      <td>Proof of Theorem 4.1</td>
      <td>47</td>
    </tr>
    <tr>
      <td>E.4</td>
      <td>Proof of Theorem 4.2</td>
      <td>50</td>
    </tr>
    <tr>
      <td>E.5</td>
      <td>Proof of Theorem 4.3</td>
      <td>50</td>
    </tr>
    <tr>
      <td>F</td>
      <td>BNP-GLLIM2: a model with an hyperprior on the gating parameters</td>
      <td>52</td>
    </tr>
    <tr>
      <td>F.1</td>
      <td>VBEM for BNP-GLLIM2</td>
      <td>52</td>
    </tr>
    <tr>
      <td>F.1.1</td>
      <td>VB-E-<i>Z</i> step</td>
      <td>54</td>
    </tr>
    <tr>
      <td>F.1.2</td>
      <td>VB-E-<i>θ</i>* step</td>
      <td>55</td>
    </tr>
    <tr>
      <td>F.1.3</td>
      <td>VB-M-<i>ρ</i> step</td>
      <td>55</td>
    </tr>
    <tr>
      <td>F.1.4</td>
      <td>ELBO for BNP-GLLIM2</td>
      <td>56</td>
    </tr>
    <tr>
      <td>F.2</td>
      <td>Predictive conditional density for BNP-GLLIM2</td>
      <td>57</td>
    </tr>
    <tr>
      <td>F.2.1</td>
      <td>Joint density</td>
      <td>57</td>
    </tr>
    <tr>
      <td>F.2.2</td>
      <td>Inverse conditional density</td>
      <td>57</td>
    </tr>
    <tr>
      <td>F.2.3</td>
      <td>Forward conditional density</td>
      <td>58</td>
    </tr>
    <tr>
      <td>F.3</td>
      <td>Proof of Proposition F.1</td>
      <td>58</td>
    </tr>
    <tr>
      <td>F.4</td>
      <td>Proof of Theorem F.2</td>
      <td>61</td>
    </tr>
    <tr>
      <td>F.5</td>
      <td>Proof of Theorem F.3</td>
      <td>66</td>
    </tr>
    <tr>
      <td>F.6</td>
      <td>Proof of Theorem F.4</td>
      <td>67</td>
    </tr>
  </tbody>
</table>

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 $\mathcal{G}$ is defined as a random process characterized by a concentration parameter $\alpha$ and a base distribution $G_0$, denoted by $G \sim \text{DP}(\alpha, G_0)$, such that for any finite partition $\{A_1, \dots, A_p\}$ of $\mathcal{G}$, the random vector $(G(A_1), \dots, G(A_p))$ is Dirichlet distributed:

$$ (G(A_1), \dots, G(A_p)) \sim \text{Dir}(\alpha G_0(A_1), \dots, \alpha G_0(A_p)). $$
<table><tbody><tr><td>C.5</td><td>VB-M-(A, b, Σ) step from Section 3.2</td><td>37</td></tr><tr><td>C.5.1</td><td>Updating b<sub>k</sub></td><td>38</td></tr><tr><td>C.5.2</td><td>Updating A<sub>k</sub></td><td>38</td></tr><tr><td>D</td><td>Details on the ELBO</td><td>41</td></tr><tr><td>E</td><td>Technical proofs</td><td>42</td></tr><tr><td>E.1</td><td>Proof of Lemma 5.2</td><td>42</td></tr><tr><td>E.2</td><td>Proof of Proposition 3.1</td><td>43</td></tr><tr><td>E.3</td><td>Proof of Theorem 4.1</td><td>47</td></tr><tr><td>E.4</td><td>Proof of Theorem 4.2</td><td>50</td></tr><tr><td>E.5</td><td>Proof of Theorem 4.3</td><td>50</td></tr><tr><td>F</td><td>BNP-GLLIM2: a model with an hyperprior on the gating parameters</td><td>52</td></tr><tr><td>F.1</td><td>VBEM for BNP-GLLIM2</td><td>52</td></tr><tr><td>F.1.1</td><td>VB-E-ζ step</td><td>54</td></tr><tr><td>F.1.2</td><td>VB-E-θ* step</td><td>55</td></tr><tr><td>F.1.3</td><td>VB-M-ρ step</td><td>55</td></tr><tr><td>F.1.4</td><td>ELBO for BNP-GLLIM2</td><td>56</td></tr><tr><td>F.2</td><td>Predictive conditional density for BNP-GLLIM2</td><td>57</td></tr><tr><td>F.2.1</td><td>Joint density</td><td>57</td></tr><tr><td>F.2.2</td><td>Inverse conditional density</td><td>57</td></tr><tr><td>F.2.3</td><td>Forward conditional density</td><td>58</td></tr><tr><td>F.3</td><td>Proof of Proposition F.1</td><td>58</td></tr><tr><td>F.4</td><td>Proof of Theorem F.2</td><td>61</td></tr><tr><td>F.5</td><td>Proof of Theorem F.3</td><td>66</td></tr><tr><td>F.6</td><td>Proof of Theorem F.4</td><td>67</td></tr></tbody></table>

## 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 $\mathcal{G}$ is defined as a random process characterized by a concentration parameter $\alpha$ and a base distribution $G_0$, denoted by $G \sim \text{DP}(\alpha, G_0)$, such that for any finite partition $\{A_1, \dots, A_p\}$ of $\mathcal{G}$, the random vector $(G(A_1), \dots, G(A_p))$ is Dirichlet distributed:

$$ (G(A_1), \dots, G(A_p)) \sim \text{Dir}(\alpha G_0(A_1), \dots, \alpha G_0(A_p)). $$
We use the stick-breaking construction of the DP G (SBDP), due to Sethuraman (1994):

$$
\begin{align*}
\theta_k^0 | G_0 &\stackrel{\text{iid}}{\sim} G_0, \quad k \in \mathbb{N}^*, \\
\tau_k | \alpha &\stackrel{\text{iid}}{\sim} \text{Beta}(\tau_k | 1, \alpha), \quad k \in \mathbb{N}^*, \\
\pi_k(\tau) &= \tau_k \prod_{l=1}^{k-1} (1 - \tau_l), \quad k \in \mathbb{N}^*, \\
G &= \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^0} \sim \text{DP}(\alpha, G_0).
\end{align*}
$$

**Pitman–Yor process.** As a generalized version of the Dirichlet process, in the Pitman–Yor process (PYP) (Pitman and Yor, 1997), the $\tau_k$'s are independent ($\stackrel{\text{iid}}{\sim}$) but not identically distributed. Specifically,

$$
\begin{align*}
\theta_k^0 | G_0 &\stackrel{\text{iid}}{\sim} G_0, \quad k \in \mathbb{N}^*, \\
\tau_k | \alpha, \sigma &\stackrel{\text{ind}}{\sim} \text{Beta}(1 - \sigma, \alpha + k\sigma), \quad k \in \mathbb{N}^*, \\
\pi_k(\tau) &= \tau_k \prod_{l=1}^{k-1} (1 - \tau_l), \quad k \in \mathbb{N}^*, \\
G &= \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^0} \sim \text{PYP}(\alpha, \sigma, G_0).
\end{align*}
$$

Here $\sigma \in (0, 1)$ is a discount parameter and $\alpha$ is a concentration parameter $\alpha > \sigma$. 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 $\sigma = 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 $O(N^\sigma)$ for PYP, while growing more slowly as $O(\log N)$ for DP.

Since the hyperparameters $\alpha$ and $\sigma$ 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 $\sigma = 0$, it is suggested in Blei and Jordan (2006) to use a gamma prior, $\alpha \sim \text{Gam}(s_1, s_2)$, where the hyperparameters $s_1$ and $s_2$ can be estimated or fixed. A natural question is whether one can also find a tractable prior for the discount parameter $\sigma$. Following the work of Lü et al. (2020), we use the following prior that satisfies the constraints $\sigma \in (0, 1)$ and $\alpha > -\sigma$,

$$p(\alpha, \sigma | s_1, s_2, a) = p(\alpha | \sigma; s_1, s_2) p(\sigma | a),$$

where $p(\alpha | \sigma; s_1, s_2)$ is a shifted gamma distribution $SG(\alpha | \sigma; s_1, s_2)$ and $p(\sigma, 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 - \sigma$, where $\sigma$ is considered fixed and $U$ follows a gamma distribution $\text{Gam}(s_1, s_2)$. The PDF of this shifted gamma
We use the stick-breaking construction of the DP G (SBDP), due to Sethuraman (1994):

$$
\begin{align*}
\theta_k^0 | G_0 &\stackrel{\text{iid}}{\sim} G_0, \quad k \in \mathbb{N}^*, \\
\tau_k | \alpha &\stackrel{\text{iid}}{\sim} \text{Beta}(\tau_k | 1, \alpha), \quad k \in \mathbb{N}^*, \\
\pi_k(\tau) &= \tau_k \prod_{l=1}^{k-1} (1 - \tau_l), \quad k \in \mathbb{N}^*, \\
G &= \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^0} \sim \text{DP}(\alpha, G_0).
\end{align*}
$$

**Pitman-Yor process.** As a generalized version of the Dirichlet process, in the Pitman-Yor process (PYP) (Pitman and Yor, 1997), the $\tau_k$'s are independent ($\stackrel{\text{iid}}{\sim}$) but not identically distributed. Specifically,

$$
\begin{align*}
\theta_k^0 | G_0 &\stackrel{\text{iid}}{\sim} G_0, \quad k \in \mathbb{N}^*, \\
\tau_k | \alpha, \sigma &\stackrel{\text{ind}}{\sim} \text{Beta}(1 - \sigma, \alpha + k\sigma), \quad k \in \mathbb{N}^*, \\
\pi_k(\tau) &= \tau_k \prod_{l=1}^{k-1} (1 - \tau_l), \quad k \in \mathbb{N}^*, \\
G &= \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^0} \sim \text{PYP}(\alpha, \sigma, G_0).
\end{align*}
$$

Here $\sigma \in (0, 1)$ is a discount parameter and $\alpha$ is a concentration parameter $\alpha > \sigma$. 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 $\sigma = 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 $O(N^\sigma)$ for PYP, while growing more slowly as $O(\log N)$ for DP.

Since the hyperparameters $\alpha$ and $\sigma$ 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 $\sigma = 0$, it is suggested in Blei and Jordan (2006) to use a gamma prior, $\alpha \sim \text{Gam}(s_1, s_2)$, where the hyperparameters $s_1$ and $s_2$ can be estimated or fixed. A natural question is whether one can also find a tractable prior for the discount parameter $\sigma$. Following the work of Lü et al. (2020), we use the following prior that satisfies the constraints $\sigma \in (0, 1)$ and $\alpha > -\sigma$,

$$p(\alpha, \sigma | s_1, s_2, a) = p(\alpha | \sigma; s_1, s_2) p(\sigma | a),$$

where $p(\alpha | \sigma; s_1, s_2)$ is a shifted gamma distribution $SG(\alpha | \sigma; s_1, s_2)$ and $p(\sigma, 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 - \sigma$, where $\sigma$ is considered fixed and $U$ follows a gamma distribution $\text{Gam}(s_1, s_2)$. The PDF of this shifted gamma
distribution is obtained from the standard gamma distribution as $\mathcal{SG}(\alpha | \sigma; s_1, s_2) = \text{Gam}(\alpha + \sigma | s_1, s_2)$.

**Hierarchical representation of BNP-MM.** BNP-MMs, including DP-MMs and PYP-MMs, see, e.g., Lü et al. (2020); Durand et al. (2022), have the following hierarchical representation to generate a data point $\mathbf{x}_n$ as a special case of BNP-GLLIMs:

**1. BNP prior:**

$$G = \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^*} \sim \text{BNP}(\alpha, \sigma, G_0), \quad \theta_k^* = (\mu_k, \mathbf{V}_k),$$

**2. BNP-MM:** for each $n \in [N]$,

$$\theta_n | G \stackrel{\text{iid}}{\sim} G,$$

If $\theta_n = \theta_k^*$, set $z_n = k$,

$$\mathbf{x}_n | z_n, \theta^* \stackrel{\text{iid}}{\sim} \mathcal{N}_L(\mathbf{x}_n | \theta_{z_n}^*), \quad \theta^* = (\theta_k^*)_{k \in [K]}.$$

## Appendix B. Variational Bayesian expectation-maximisation principle

The clustering task consists mainly of estimating the unknown labels $\mathcal{Z} = (z_n)_{n \in [N]}$ from the observed data $(\mathcal{Y}, \mathcal{X}) = (y_n, \mathbf{x}_n)_{n \in [N]}$, whose joint distribution $p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)$ is determined by a set of parameters denoted by $\Theta$ and often by additional hyperparameters $\phi$.

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 $\Theta$. 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 $\mathcal{F}_0$ defined, for any probability distribution $q_{\mathcal{Z}}$ over labels $\mathcal{Z}$, by

$$\begin{aligned} \mathcal{F}_0(q_{\mathcal{Z}}, \Theta, \phi) &= \sum_{\mathcal{Z}} q_{\mathcal{Z}}(\mathcal{Z}) \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z} | \Theta, \phi) - \mathbb{E}_{q_{\mathcal{Z}}}[\log q_{\mathcal{Z}}(\mathcal{Z})] \\ &= \mathbb{E}_{q_{\mathcal{Z}}}\left[\log \frac{p(\mathcal{Y}, \mathcal{X}, \mathcal{Z} | \Theta, \phi)}{q_{\mathcal{Z}}(\mathcal{Z})}\right], \end{aligned}$$

where $-\mathbb{E}_{q_{\mathcal{Z}}}[\log q_{\mathcal{Z}}(\mathcal{Z})]$ is the entropy of $q_{\mathcal{Z}}$ and $\mathbb{E}_q[\cdot]$ is the expectation with respect to $q$. The function $\mathcal{F}_0$ depends on the observations $(\mathcal{Y}, \mathcal{X})$, 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 $\Theta$ using the prior knowledge encoded in a distribution $p(\Theta)$. More precisely, the MLE of $\Theta$ is then replaced by a point estimation $\hat{\Theta} = \operatorname*{argmax}_{\Theta \in \Theta} p(\Theta | \mathcal{Y}, \mathcal{X})$. In this paper, instead of considering only point estimation of $\Theta$, we carry out a fully Bayesian approach. That
distribution is obtained from the standard gamma distribution as $\mathcal{SG}(\alpha | \sigma; s_1, s_2) = \text{Gam}(\alpha + \sigma | s_1, s_2)$.

**Hierarchical representation of BNP-MM.** BNP-MMs, including DP-MMs and PYP-MMs, see, e.g., Lü et al. (2020); Durand et al. (2022), have the following hierarchical representation to generate a data point $\mathbf{x}_n$ as a special case of BNP-GLLIMs:

**1. BNP prior:**

$$G = \sum_{k=1}^{\infty} \pi_k(\tau) \delta_{\theta_k^*} \sim \text{BNP}(\alpha, \sigma, G_0), \quad \theta_k^* = (\mu_k, V_k),$$

**2. BNP-MM:** for each $n \in [N]$,

$$\theta_n | G \stackrel{\text{iid}}{\sim} G,$$

If $\theta_n = \theta_k^*$, set $z_n = k$,

$$\mathbf{x}_n | z_n, \theta^* \stackrel{\text{iid}}{\sim} \mathcal{N}_L(\mathbf{x}_n | \theta_{z_n}^*), \quad \theta^* = (\theta_k^*)_{k \in [K]}.$$

## Appendix B. Variational Bayesian expectation-maximisation principle

The clustering task consists mainly of estimating the unknown labels $\mathcal{Z} = (z_n)_{n \in [N]}$ from the observed data $(\mathcal{Y}, \mathcal{X}) = (y_n, \mathbf{x}_n)_{n \in [N]}$, whose joint distribution $p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)$ is determined by a set of parameters denoted by $\Theta$ and often by additional hyperparameters $\phi$.

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 $\Theta$. 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 $\mathcal{F}_0$ defined, for any probability distribution $q_{\mathcal{Z}}$ over labels $\mathcal{Z}$, by

$$\begin{aligned} \mathcal{F}_0(q_{\mathcal{Z}}, \Theta, \phi) &= \sum_{\mathcal{Z}} q_{\mathcal{Z}}(\mathcal{Z}) \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z} | \Theta, \phi) - \mathbb{E}_{q_{\mathcal{Z}}}[\log q_{\mathcal{Z}}(\mathcal{Z})] \\ &= \mathbb{E}_{q_{\mathcal{Z}}}\left[\log \frac{p(\mathcal{Y}, \mathcal{X}, \mathcal{Z} | \Theta, \phi)}{q_{\mathcal{Z}}(\mathcal{Z})}\right], \end{aligned}$$

where $-\mathbb{E}_{q_{\mathcal{Z}}}[\log q_{\mathcal{Z}}(\mathcal{Z})]$ is the entropy of $q_{\mathcal{Z}}$ and $\mathbb{E}_q[\cdot]$ is the expectation with respect to $q$. The function $\mathcal{F}_0$ depends on the observations $(\mathcal{Y}, \mathcal{X})$, 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 $\Theta$ using the prior knowledge encoded in a distribution $p(\Theta)$. More precisely, the MLE of $\Theta$ is then replaced by a point estimation $\hat{\Theta} = \underset{\Theta \in \Theta}{\operatorname{argmax}} p(\Theta | \mathcal{Y}, \mathcal{X})$. In this paper, instead of considering only point estimation of $\Theta$, we carry out a fully Bayesian approach. That
is, we integrate out $\Theta$ as follows

$$p(\mathcal{Z} | \mathcal{Y}, \mathcal{X}) = \int_{\Theta} p(\mathcal{Z} | \mathcal{Y}, \mathcal{X}, \Theta)p(\Theta | \mathcal{Y}, \mathcal{X})d\Theta.$$

This integration requires the computation of the density $p(\Theta | \mathcal{Y}, \mathcal{X})$, 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(\Theta | \mathcal{Y}, \mathcal{X})$ and $p(\mathcal{Z} | \mathcal{Y}, \mathcal{X})$. 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_{\mathcal{Z}}$ and $q_{\Theta}$ denote the distributions over $\mathcal{Z}$ and $\Theta$, respectively, which will serve as approximations to the true posteriors. Specifically, in the Bayesian setting, the intractable posterior $p(\mathcal{Z}, \Theta | \mathcal{Y}, \mathcal{X}; \phi)$ is approximated by the variational posterior $q(\mathcal{Z}, \Theta) = q_{\mathcal{Z}}(\mathcal{Z})q_{\Theta}(\Theta)$.

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_{\mathcal{Z}}$ and $q_{\Theta}$ distributions by

$$ \begin{aligned} \mathcal{F}(q_{\mathcal{Z}}, q_{\Theta}, \phi) &= \mathbb{E}_{q_{\mathcal{Z}}q_{\Theta}} \left[ \log \frac{p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)}{q_{\mathcal{Z}}(\mathcal{Z})q_{\Theta}(\Theta)} \right] \quad (B1) \\ &= \log p(\mathcal{Y}, \mathcal{X} | \phi) - \text{KL}(q_{\mathcal{Z}}q_{\Theta} \| p(\mathcal{Z}, \Theta | \mathcal{Y}, \mathcal{X}, \phi)) \le \log p(\mathcal{Y}, \mathcal{X} | \phi), \end{aligned} $$

alternatively over $q_{\mathcal{Z}}$, $q_{\Theta}$ and $\phi$. Here, KL stands for Kullback-Leibler divergence. It is worth noting that adding a prior on $\Theta$ is formally equivalent to considering $\Theta$ as missing variables, while the hyperparameters $\phi$ play the role of the parameters of interest in MLE.

The alternate maximisation over $\mathcal{F}$ 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_{\mathcal{Z}}$ and $q_{\Theta}$ 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 $\mathcal{F}(q_{\mathcal{Z}}, q_{\Theta}, \phi)$. At the $r$th iteration, using current values $\phi^{(r-1)}$ and $q_{\Theta}^{(r-1)}$, we get the following updating,

$$ \text{VB-E-}\mathcal{Z}: q_{\mathcal{Z}}^{(r)}(\mathcal{Z}) \propto \exp \mathbb{E}_{q_{\Theta}^{(r-1)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi^{(r-1)}) \right], \quad (B2) $$

$$ \text{VB-E-}\Theta: q_{\Theta}^{(r)}(\Theta) \propto \exp \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi^{(r-1)}) \right], \quad (B3) $$

$$ \text{VB-M-}\phi: \phi^{(r)} \propto \underset{\phi}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}q_{\Theta}^{(r)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi) \right]. $$

In practice, we can decide which parameters to treat as genuine parameters $\Theta$ or as hyperparameters $\phi$, depending on whether some prior knowledge is available for only a subset of the parameters, or whether the model has hyperparameters $\phi$ for which no prior information is available. Furthermore, for complex models, $q_{\Theta}$ and $q_{\mathcal{Z}}$ 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 Appendix F.1 for the BNP-GLLIM inference.
is, we integrate out $\Theta$ as follows

$$p(\mathcal{Z} | \mathcal{Y}, \mathcal{X}) = \int_{\Theta} p(\mathcal{Z} | \mathcal{Y}, \mathcal{X}, \Theta)p(\Theta | \mathcal{Y}, \mathcal{X})d\Theta.$$

This integration requires the computation of the density $p(\Theta | \mathcal{Y}, \mathcal{X})$, 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(\Theta | \mathcal{Y}, \mathcal{X})$ and $p(\mathcal{Z} | \mathcal{Y}, \mathcal{X})$. 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_{\mathcal{Z}}$ and $q_{\Theta}$ denote the distributions over $\mathcal{Z}$ and $\Theta$, respectively, which will serve as approximations to the true posteriors. Specifically, in the Bayesian setting, the intractable posterior $p(\mathcal{Z}, \Theta | \mathcal{Y}, \mathcal{X}; \phi)$ is approximated by the variational posterior $q(\mathcal{Z}, \Theta) = q_{\mathcal{Z}}(\mathcal{Z})q_{\Theta}(\Theta)$.

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_{\mathcal{Z}}$ and $q_{\Theta}$ distributions by

$$ \begin{aligned} \mathcal{F}(q_{\mathcal{Z}}, q_{\Theta}, \phi) &= \mathbb{E}_{q_{\mathcal{Z}}q_{\Theta}} \left[ \log \frac{p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi)}{q_{\mathcal{Z}}(\mathcal{Z})q_{\Theta}(\Theta)} \right] \quad (B1) \\ &= \log p(\mathcal{Y}, \mathcal{X} | \phi) - \text{KL}(q_{\mathcal{Z}}q_{\Theta} \| p(\mathcal{Z}, \Theta | \mathcal{Y}, \mathcal{X}, \phi)) \le \log p(\mathcal{Y}, \mathcal{X} | \phi), \end{aligned} $$

alternatively over $q_{\mathcal{Z}}$, $q_{\Theta}$ and $\phi$. Here, KL stands for Kullback-Leibler divergence. It is worth noting that adding a prior on $\Theta$ is formally equivalent to considering $\Theta$ as missing variables, while the hyperparameters $\phi$ play the role of the parameters of interest in MLE.

The alternate maximisation over $\mathcal{F}$ 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_{\mathcal{Z}}$ and $q_{\Theta}$ 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 $\mathcal{F}(q_{\mathcal{Z}}, q_{\Theta}, \phi)$. At the $r$th iteration, using current values $\phi^{(r-1)}$ and $q_{\Theta}^{(r-1)}$, we get the following updating,

$$ \text{VB-E-} \mathcal{Z}: q_{\mathcal{Z}}^{(r)}(\mathcal{Z}) \propto \exp \mathbb{E}_{q_{\Theta}^{(r-1)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi^{(r-1)}) \right], \quad (B2) $$

$$ \text{VB-E-} \Theta: q_{\Theta}^{(r)}(\Theta) \propto \exp \mathbb{E}_{q_{\mathcal{Z}}^{(r)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi^{(r-1)}) \right], \quad (B3) $$

$$ \text{VB-M-} \phi: \phi^{(r)} \propto \underset{\phi}{\operatorname{argmax}} \mathbb{E}_{q_{\mathcal{Z}}^{(r)}q_{\Theta}^{(r)}} \left[ \log p(\mathcal{Y}, \mathcal{X}, \mathcal{Z}, \Theta; \phi) \right]. $$

In practice, we can decide which parameters to treat as genuine parameters $\Theta$ or as hyperparameters $\phi$, depending on whether some prior knowledge is available for only a subset of the parameters, or whether the model has hyperparameters $\phi$ for which no prior information is available. Furthermore, for complex models, $q_{\Theta}$ and $q_{\mathcal{Z}}$ 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 Appendix F.1 for the BNP-GLLIM inference.