Label Switching Problem in Bayesian Analysis for Gravitational Wave Astronomy

TL;DR

This paper addresses the label switching problem in Bayesian analysis, proposing a hypertriangle mapping solution that improves sampling efficiency for multi-modal posteriors in gravitational wave astronomy.

Contribution

It introduces a novel hypertriangle mapping method to effectively resolve label switching in Bayesian models with multiple parameters.

Findings

The method efficiently handles large parameter sets.

It is compatible with existing stochastic sampling algorithms.

Demonstrated effectiveness in gravitational wave astronomy examples.

Abstract

The label switching problem arises in the Bayesian analysis of models containing multiple indistinguishable parameters with arbitrary ordering. Any permutation of these parameters is equivalent, therefore models with many such parameters have extremely multi-modal posterior distributions. It is difficult to sample efficiently from such posteriors. This paper discusses a solution to this problem which involves carefully mapping the input parameter space to a high dimensional hypertriangle. It is demonstrated that this solution is efficient even for large numbers of parameters and can be easily applied alongside any stochastic sampling algorithm. This method is illustrated using two example problems from the field of gravitational wave astronomy.

Tables

Table 1. Table 1: Log-evidences for mixtures with different number of components K 𝐾 K . The variables Z 𝒯 subscript 𝑍 𝒯 Z_{\mathcal{T}} and Z 𝒞 subscript 𝑍 𝒞 Z_{\mathcal{C}} denote the evidences obtained by sampling on the hypertriangle parameter space in Eq.  14 and the (multimodal) hypercube parameter space in Eq.  15 respectively. Mathematically we have already proved that these parameter spaces are equivalent and therefore the two evidences are equal; these two columns serve to demonstrate this numerically. For the K = 1 𝐾 1 K=1 component case there is no distinction between the two parameter spaces (the map ϕ italic-ϕ \phi reduces to the identity in this case). The errors on the CPNest evidence integrals were estimated by a combination of the internal CPNest error estimate (as described in Skilling ( 2004 ) ) and examination of the spread of results from multiple runs. The Z 𝒯 subscript 𝑍 𝒯 Z_{\mathcal{T}} and Z 𝒞 subscript 𝑍 𝒞 Z_{\mathcal{C}} evidences are broadly consistent; however for large K 𝐾 K there is some tension. We think this is due to CPNest systematically underestimating the Z 𝒞 subscript 𝑍 𝒞 Z_{\mathcal{C}} evidence which comes from a high dimensional and highly multimodal posterior. Alternative nested samplers Feroz  et al. ( 2009 ); Handley  et al. ( 2015 ); Brewer and Foreman-Mackey ( 2016 ); Speagle ( 2019 ) have been shown to reliably estimate evidences for problems of similar complexity.

$K$	$\log Z_{𝒯}$	$\log Z_{𝒞}$
$1$	$-74.76\pm 0.09$
$2$	$-78.37\pm 0.05$	$-78.30\pm 0.05$
$3$	$-81.82\pm 0.09$	$-81.66\pm 0.08$
$4$	$-84.58\pm 0.06$	$-84.2\pm 0.1$

Table 2. Table 2: The parameters of the six injected GBs, with f ⋆ = 1 ​ mHz subscript 𝑓 ⋆ 1 mHz f_{\star}=1\,\mathrm{mHz} . The amplitudes were chosen such that the signal-to-noise ratio is 10 in each case.

$f-f_{\star }$	${\log }_{10}A$	$\iota$[rad]	$\lambda$ [rad]	$\beta$[rad]	$\psi$[rad]	$\varphi$[rad]	$\rho$
$0$	-22.15	0.246	-0.096	0.218	1.640	1.795	10
$2/\mathrm{yr}$	-22.13	0.403	0.091	0.294	1.066	4.249	10
$4/\mathrm{yr}$	-22.13	0.376	-0.055	0.359	0.794	4.760	10
$6/\mathrm{yr}$	-22.15	0.284	0.031	0.248	1.127	2.078	10
$8/\mathrm{yr}$	-22.13	0.390	0.006	0.223	0.775	4.537	10
$10/\mathrm{yr}$	-22.12	0.428	0.091	0.296	1.088	5.765	10

Table 3. Table S1: Priors on the BBH observed population parameters.

${\mu }_k^{\left(\log m_1\right)}$	${\mu }_k^{\left(\log m_1\right)}$	${\lambda }_k^1$	${\lambda }_k^2$	${\varphi }_k$ [rad]	$w_k$
(0, 2)	(0, 2)	(0.01, 4)	(0.01, 4)	(0, $\pi /2$)	(0, 1)

Table 4. Table S2: Priors on the GB parameters. The frequency prior spans twelve bins around f ⋆ = 1 ​ mHz subscript 𝑓 ⋆ 1 mHz f_{\star}=1\,\mathrm{mHz} .

$\left(f-f_{\star }\right)$ $\left[{\text{yr}}^{-1}\right]$	${\log }_{10}A$	$\lambda$ [rad]	$\sin \beta$	$\cos \iota$	$\psi$ [rad]	$\Phi$ [rad]
(-1, 11)	(-23.0, -21.8)	(0, 1)	(-0.75, 0.75)	(0,1)	(0, $\pi$)	(0, $2\pi$)

Equations

$$\varphi :\mathcal{C}\to \mathcal{T}.$$

$$x_1^{\prime }$$

$$x_2^{\prime }$$

$$x_1^{\prime }$$

$$x_2^{\prime }$$

$$x_i^{\prime }=x_{i-1}^{\prime }+(1-x_{i-1}^{\prime })\left[1-(1-x_i{)}^{\frac{1}{K+1-i}}\right],$$

$$x_i^{\prime }=1-\prod _{j=1}^i(1-x_j{)}^{\frac{1}{K+1-j}}.$$

$$x_i^{\prime }$$

$$=\left[1-{\left(1-x_i\right)}^{\frac{1}{K+1-i}}\right]-x_{i-1}^{\prime }\left[1-{\left(1-x_i\right)}^{\frac{1}{K+1-i}}-1\right]$$

$$=\left[1-{\left(1-x_i\right)}^{\frac{1}{K+1-i}}\right]+{\left(1-x_i\right)}^{\frac{1}{K+1-i}}{x}_{i-1}^{\prime }$$

$$=1-{\left(1-x_i\right)}^{\frac{1}{K+1-i}}\left(1-x_{i-1}^{\prime }\right).$$

$$x_i^{\prime }$$

$$=\dots$$

$$=1-\prod _{j=1}^i{\left(1-x_j\right)}^{\frac{1}{K+1-j}}.$$

$$J$$

$$\times \frac{1}{\left(1-x_1\right)}\left[{\left(1-x_1\right)}^{\frac{1}{K+1-1}}\right]$$

$$\times \frac{1}{\left(1-x_2\right)}\left[{\left(1-x_1\right)}^{\frac{1}{K+1-1}}{\left(1-x_2\right)}^{\frac{1}{K+1-2}}\right]$$

$$\times \dots$$

$$\times \frac{1}{\left(1-x_K\right)}\left[{\left(1-x_1\right)}^{\frac{1}{K+1-1}}\dots {\left(1-x_K\right)}^{\frac{1}{K+1-K}}\right]$$

$$J=\frac{1}{K!}.$$

$$\Pi (x_1,\dots ,x_K)=\prod _{k=1}^K\pi (x_k)$$

$$F(x)={\int }_0^x\pi (s)\text{d}s.$$

$$\left[\begin{array}{c}\log m_1\\ \log m_2\end{array}\right]\sim \sum _{k=1}^K{w}_k\mathcal{N}\left(\left[\begin{array}{c}{\mu }_k^{(\log m_1)}\\ {\mu }_k^{(\log m_2)}\end{array}\right],{\mathbf{\Sigma }}_k\right).$$

$${\mathbf{\Lambda }}_k=\left({\mu }_k^{(\log m_1)},{\mu }_k^{(\log m_2)},{\lambda }_k^1,{\lambda }_k^2,{\varphi }_k,w_k\right).$$

$${\mathbf{\Lambda }}_k=\left({\chi }_k,{\mu }_k^{(\log m_2)},{\lambda }_k^1,{\lambda }_k^2,{\varphi }_k,w_k\right).$$

$$p\left(\{d\}\mid \mathbf{\Lambda }\right)=\prod _{i=1}^{N_{\text{obs}}}\frac{\int \mathrm{d}\theta p(d\mid \theta )p_{\text{pop}}(\theta \mid \mathbf{\Lambda })}{\int \mathrm{d}\theta p_{\text{pop}}(\theta \mid \mathbf{\Lambda })}.$$

$$p\left(\mathbf{\Lambda }\mid \{d\}\right)=\varpi (\mathbf{\Lambda })\prod _{i=1}^{N_{\text{obs}}}\frac{\frac{1}{N_i}{\sum }_{j=1}^{N_i}\frac{p_{\text{pop}}({\theta }_i^j\mid \mathbf{\Lambda })}{\pi ({\theta }_i^j)}}{\int \mathrm{d}\theta p_{\text{pop}}(\theta \mid \mathbf{\Lambda })},$$

$$h_{+}(t;\mathbf{\Lambda })$$

$$h_{\times }(t;\mathbf{\Lambda })$$

$${\mathbf{\Lambda }}_k=\{{\log }_{10}{A}_k,f_k,{\lambda }_k,\sin {\beta }_k,\cos {\iota }_k,{\psi }_k,{\Phi }_k\}.$$

$$\log \mathcal{L}({\mathbf{\Lambda }}_k)\propto -\frac{1}{2}\sum _{\alpha }{\left|s_{\alpha }-\sum _{k=1}^K{h}_{\alpha }({\mathbf{\Lambda }}_k)\right|}_{(\alpha )}^2,$$

$${⟨a|b⟩}_{(\alpha )}=4\Re \left\{{\int }_0^{\infty }\mathrm{d}f\frac{\tilde{a}(f)\tilde{b}(f)}{S_{\alpha }(f)}\right\},$$

$${\mathbf{\Lambda }}_k=\{{\log }_{10}{A}_k,{\chi }_k,\cos {\iota }_k,{\lambda }_k,\sin {\beta }_k,{\psi }_k,{\varphi }_k\}.$$

Full text

I Introduction

It sometimes occurs in Bayesian inference problems that the target distribution depends on several parameters whose ordering is arbitrary. Three examples are immediately apparent from the field of gravitational wave (GW) astronomy alone. Firstly, when describing a compact binary with component masses $m_{1}$ and $m_{2}$, the likelihood is symmetric under exchange of the labels 1 and 2 (provided all other relevant parameters are suitably adjusted simultaneously). Secondly, when analysing GW time series data containing two or more overlapping sources of the same type, the likelihood is invariant under exchanging all of the parameters of any pair of sources. And thirdly, when analysing the parameters of a population of observed GW events, mixture models can be used to model the population and/or to infer the presence of distinct astrophysical populations. In this case the hyper-likelihood for the population parameters may be invariant under exchanging the parameters of the population components.

Sometimes a simple reparametrisation and restricting the parameter range is enough to remove the degeneracy arising from the arbitrary ordering. In the first case of the binary with two component masses, it is possible to define, say, the total mass $M=m_{1}+m_{2}$ and mass ratio $q=m_{2}/m_{1}$ and to sample these over the restricted ranges $M>0$ and $q\leq 1$. This covers only the restricted portion of the parameter space $m_{1}\geq m_{2}$, thereby removing the symmetry from the likelihood.

The second and third examples are more problematic as they are not restricted to just 2 degrees of freedom. In each case the target distribution has a high degree of symmetry and is invariant under permutations of some number of labels, $K$. A great deal of literature is devoted to this label switching problem in the context of mixture models Redner and Walker (1984); Diebolt and Robert (1994); Richardson and Green (1997); Stephens (2000); Jasra et al. (2005); Hurn et al. (2003); Fruhwirth-Schnatter (2001). The invariance of the target distribution under permutations means that if the posterior has a peak (or mode) at a particular point in parameter space it will necessarily have peaks at all $K!$ points related by symmetry. The extreme scaling of this multimodality poses a serious obstacle to any sampling algorithm in moderate or high dimensional problems.

The most natural way to solve the label switching problem is to impose an artificial identifiability constraint. Searching over the restricted region $m_{1}\geq m_{2}$ of the binary component mass space is an example of such a constraint in 2 dimensions. In the $K$ dimensional problem this can be generalised by demanding a certain ordering of the parameters; see, for example, Diebolt and Robert (1994); Stephens (2000); Jasra et al. (2005); Hurn et al. (2003). Restricting to this small region of parameter space avoids all symmetries and removes the excess multimodality. It is also obvious that if one can adequately explore the restricted parameter space satisfying the artificial identifiability constraint then, by symmetry, this is equivalent to exploring the full space.

It remains to implement a suitable artificial identifiability constraint in practical inference problems. This problem can be approached in several ways. For example, when using an Markov chain Monte Carlo (MCMC) to explore the target distribution the proposal can be augmented by composing with a sorting function; i.e. propose a point then reorder the parameters such that the constraint is satisfied Fruhwirth-Schnatter (2001). Alternatively, the log-prior distribution can be crudely modified so that it returns $-\infty$ for any point not satisfying the constraint. Either of these will ensure the chain never leaves the desired region of parameter space.

While undoubtedly simple, neither of these approaches are completely satisfactory. The former approach requires the user to modify their MCMC proposal distribution and it is difficult to apply when using other stochastic sampling algorithms which don’t have a user-accessible proposal distribution (such as nested sampling Skilling (2004)). For this reasons such an approach is not compatible with the modern approach of treating the sampler, as far as possible, as a black box to which the user must only provide a likelihood and a prior. The latter approach is easy to implement for all samplers, but has the significant drawback of being extremely inefficient in high numbers of dimensions. This is because the sampler only proposes useful points satisfying the identifiability constraint a tiny fraction $1/K!$ of the time.

This paper presents a solution to the label switching problem. Our approach follows that of Handley et al. (2015) (see, in particular, Eq. A13; however, this equation contains a typographical error as pointed out by Greif et al. (2019)). This solution to the label switching problem has been implemented in Handley et al. (2019a) and has been widely used extensively in the astronomical and cosmological literature Handley et al. (2019b); Greif et al. (2019); Higson et al. (2019); Planck Collaboration et al. (2018); Hall et al. (2018); Millea and Bouchet (2018); Olamaie et al. (2018); Finelli et al. (2018); Hee et al. (2017); Planck Collaboration et al. (2016); Hee et al. (2016). Here we describe the solution in detail, including proofs of certain important properties of the solution. This solution is mathematically elegant, efficient in high dimensions, and can be easily integrated with any sampling algorithm while treating it as a black box.

In Sec. II the label switching problem is described in detail and the idea behind our proposed solution is illustrated in 2 dimensions. Our solution, for an arbitrary number of dimensions, is presented in Sec. III. In Sec. IV the efficacy of our proposed solution is demonstrated by applying it to the second (Sec. IV.1) and third (Sec. IV.2) example problems described in the opening paragraph of this section. These example applications are drawn from the field of GW astronomy, but we stress that this method has been more generally applied to inference in astronomy already.

II The Label Switching Problem

We wish to treat problems containing multiple indistinguishable components. Each of the $K$ components is modeled by some parameters $\mathbf{\Lambda}_{k}\in U$, where the parameter space $U$ is an open set of $\mathbb{R}^{n}$ and $k\in\{1,2,\ldots,K\}$. We will choose to distinguish components based on the values of one of these parameters, $x_{k}\equiv\mathbf{\Lambda}_{k}^{1}\in I$ where $I$ is an set of $\mathbb{R}$. For simplicity, in this section we will consider $x_{k}\in(0,1)$ and use a flat prior on each $x_{k}$, although these restrictions can be relaxed later.

In the case where there are two components, $K=2$, our full parameter space is $U\times U$. However we will mainly be interested in the subspace spanned by $\vec{x}=(x_{1},x_{2})$, which covers the unit square $I\times I$ (in the general $K$-dimensional case this will be a hypercube which will be denoted $\mathcal{C}$). For the remainder of this section we suppress the other components of $\mathbf{\Lambda}_{k}$ from our notation for clarity.

Since the two components are indistinguishable, the points $(x_{1},x_{2})$ and $(x_{2},x_{1})$ are equivalent; both the likelihood, $\mathcal{L}(\vec{x})$, and prior distributions are symmetric under interchange of the labels 1 and 2 (provided we also remember to relabel all the other components of $\mathbf{\Lambda}_{1}$ and $\mathbf{\Lambda}_{2}$ simultaneously). As a result, the parameter space is twice as large as it needs to be. Evaluating $\mathcal{L}(\vec{x})$ over the square will typically lead to a distribution with two global maxima (an exception occurs when the true maximum is on the boundary $x_{1}=x_{2}$); secondary peaks, ridges and other structures in the likelihood are also duplicated. In higher dimensions this duplication and multimodality increases in proportion to $K!$ and becomes a serious obstacle to sampling the target distribution.

To avoid sampling multiple identical copies of the same likelihood modes we will enforce the identifiability constraint $x_{2}\geq x_{1}$. This amounts to labelling the component with the smallest $x$ as #1, the component with the next largest $x$ as #2, and so on in higher dimensions. In two dimensions, this restricts the parameter space to the triangle $x_{2}\geq x_{1}$ (see the off principal diagonal panels in Fig. 1). In higher dimensions, the parameter space is restricted to the region $x_{K}\geq x_{K-1}\geq\cdots\geq x_{1}$, which is hereafter referred to as the hypertriangle and denoted $\mathcal{T}$.

Samplers naturally propose points in a hypercube. To avoid modifying the sampler itself, we wish to map points in the hypercube to points in the hypertriangle (following the strategy first introduced in the astronomy literature by Feroz et al. (2009)):

[TABLE]

Naively, we might try to choose $\phi$ to be the sorting function. Unfortunately, although it does map into the hypertriangle, it doesn’t solve the multimodality problem, since sorting is a many-to-one map. If the sampler proposes a point $\vec{x}=(x_{1},x_{2})$ in the hypercube and then the user applies the sorting function $\vec{x}^{\prime}=\mathrm{sort}(\vec{x})$ before evaluating the likelihood $\mathcal{L}(\vec{x}^{\prime})$, nothing restricts the sampler from searching over the full hypercube. In fact, this procedure is identical to sampling the original hypercube with no sorting.

This is to be distinguished from the procedure of sorting inside the proposal distribution, as referenced in the introduction, which does restrict sampling to the hypertriangle. This is because the newly sorted points are kept by the sampler and used for generating the next set of proposed points. However, this approach violates our desire to treat samplers as black boxes.

To solve the problem we seek a function, $\phi$, which is one-to-one. One possibility, in 2 dimensions, is to leave the $x_{1}$ coordinate invariant and shift/rescale the $x_{2}$ coordinate such that it lies in the desired range:

[TABLE]

To see that points are indeed mapped to $\mathcal{T}$ it is sufficient to note that the correct ordering is enforced by adding a positive quantity to $x^{\prime}_{1}$ to get a larger value for $x^{\prime}_{2}$. The range $x^{\prime}_{2}\in(0,1)$ is in turn ensured by scaling $x_{2}$ with the factor $1-x^{\prime}_{1}$. This map is a indeed one-to-one map from the square to the triangle, thereby removing the problem of multiple modes. However, this map has the unfortunate property that it distorts the prior on the $x_{2}$ component, favoring larger values (see the red distribution in Fig. 1).

The map in Eq. 2 can be “fixed” by revising the $x_{1}$ coordinate downwards, before shifting/rescaling $x_{2}$:

[TABLE]

The new map in Eq. 3 solves the problem in 2 dimensions (see the blue distribution in Fig. 1). The sampler can propose points $\vec{x}=(x_{1},x_{2})$ uniformly in the square. These points are mapped to the triangle $\vec{x}^{\prime}=\phi(\vec{x})$. Finally, the likelihood is evaluated at the mapped points, $\mathcal{L}(\vec{x}^{\prime})$. This procedure correctly covers the parameter space just once with the desired flat prior. To prove that the proposed map does indeed maintain the desired flat prior on the individual components one can evaluate the Jacobian of the transformation $\vec{x}^{\prime}=\phi(\vec{x})$ and show that is constant. This is done in the next section for the $K$-dimensional case. Because the Jacobian is constant, this transformation will correctly preserve the flat prior that is imposed on the original $x_{k}$. 111For an extension of our solution to the wider class of separable priors, see Sec. III.3.

To state the problem formally: to solve the label-switching problem, we seek a bijection (a “one-to-one” and “onto” map) ${\phi:\mathcal{C}\rightarrow\mathcal{T}}$, for an arbitrary number of dimensions, with components ${x^{\prime}_{\kappa}=\phi_{\kappa}(x_{k})}$, such that the Jacobian $J\equiv\left|\partial x^{\prime}_{\kappa}/\partial x_{k}\right|=\mathrm{const}$.

III The hypertriangle map in arbitrary dimensions

Our proposed generalization of the 2-dimensional map in Eq.3, $x^{\prime}=\phi(x)$, is defined recursively as

[TABLE]

where $i\in{1,\ldots,K}$ and $x^{\prime}_{0}=0$ by definition. This closely resembles Equation (A13) of Handley et al. (2015), although here we have corrected a typographical error. Eq. 4 can be expressed non-recursively as:

[TABLE]

If the inputs are in the correct range $x_{j}\in(0,1)$, i.e. $x_{j}\in\mathcal{C}$, it can be shown that the output falls in $\mathcal{T}$ (the logic as outlined in Sec. II for Eq. 2 still applies). It can also be shown that this map is a bijection by inverting Eq. 5.

In the remainder of the section, we will first prove that Eq.5 is equivalent to Eq.4, and then that the Jacobian of Eq.5 is constant.

III.1 Equivalent Representations of $\phi$

Starting with the recursive version of the map given in Eq.4, we rearrange it as follows:

[TABLE]

This procedure can be repeated for the $x_{i-1}^{\prime}$ term inside the final set of parentheses, and then for $x_{i-2}^{\prime}$ and so on down to $x^{\prime}_{1}$. This gives the equivalent representation to Eq. 5;

[TABLE]

III.2 The Jacobian of $\phi$

As discussed in Sec. II, to maintain the correct prior on the hypertriangle, it is necessary that the map $\phi$ has a constant Jacobian. To prove that our proposed hypertriangle map has this property, we start with the form of the map in Eq.5. The Jacobian matrix for this specific transformation is lower-triangular because the component $x^{\prime}_{i}$ depends only on $x_{j}$ with $j\leq i$. Its determinant is therefore equal to the product of the diagonal terms:

[TABLE]

where in the final step a factor has been moved inside of the second product and the upper limit of the product has been changed accordingly. Writing out the products explicitly gives

[TABLE]

Careful counting of all the terms reveals that everything cancels and we are left with

[TABLE]

The Jacobian is equal to one over the number of times the original parameter space was covered by the hypercube.

III.3 Extension to separable priors

The above derivation considered only flat priors on the $x_{k}$. Here we consider the applicability of our hypertriangulation map to separable priors of the form

[TABLE]

In such cases it is first necessary to transform to new coordinates such that the prior is flat before proceeding to apply the hypertriangulation map as before.

In order to find the new coordinates with flat priors, first evaluate the cumulative distribution function

[TABLE]

Then define new coordinates $y_{k}=F(x_{k})$ which lie in the range $[0,1]$. The prior on these new coordinates is now flat and the hypertriangulation map may now be applied to the $y_{k}$.

III.4 Implementation of $\phi$

For concreteness, we provide here a pseudo-Python implementation of Eq. 5. The input x (in $\mathcal{C}$) and output X (in $\mathcal{T}$) are numpy arrays where all values are in the prior range $(0,1)$. The values of x may be in any order whilst the values of X are, by construction, in ascending order. If a different prior range is needed then the input and output must be shifted and rescaled as appropriate. A full Python implementation (including the shifting and rescaling) is provided at the GitHub repository Buscicchio et al. (2019a).

IV Example GW Applications

In this section we present two applications of our hypertriangle method to two rather different Bayesian inference problems drawn from the field of GW astronomy.

The first example in Sec. IV.1 is a Gaussian mixture model; models of this type have been studied extensively in the context of the label switching problem Hurn et al. (2003); Richardson and Green (1997); Jasra et al. (2005); Fruhwirth-Schnatter (2001); Stephens (2000); Diebolt and Robert (1994); Redner and Walker (1984).

The second example in Sec. IV.2 involves the identification of multiple overlapping signals in time series data. The label switching problem has not often been explicitly considered in this context. However Gillard and Kvasov (2017) discuss it when fitting multiple damped sinusoids to time series data.

IV.1 The Observed Mass Function of LIGO/Virgo Binary Black Holes

LIGO and Virgo LIGO Scientific Collaboration et al. (2015); Acernese et al. (2015) are ground-based GW detectors operating in the $(10^{1}-10^{4})\,\mathrm{Hz}$ frequency range. The network has been operating since September 2015 and has so far confidently detected 10 binary black hole (BBH) mergers and 1 binary neutron star merger The LIGO Scientific Collaboration et al. (2018a). The third observation run is ongoing and low latency pipelines Sachdev et al. (2019); Usman et al. (2015); Hooper et al. (2011); Adams et al. (2015); Klimenko et al. (2008) have produced a number of public alerts associated with event candidates Gra (oorg); Cooper et al. (acuk). It is likely that by the end of the current run dozens more detections will be available Abbott et al. (2018) for further investigation.

Detailed waveform models for BBH signal calibrated against numerical relativity are now available Bohé et al. (2017); Husa et al. (2016); Khan et al. (2016); Hannam et al. (2014); Dal Canton and Harry (2017). These are used in the LALInference Bayesian analysis software package Veitch et al. (2015) to construct posterior distributions on the parameters of each event. These include both intrinsic (component masses, spins, angular momentum, etc) and extrinsic (sky position, distance, inclination) parameters.

Of these parameters, the best measured and most astrophysically interesting are the individual black hole masses. Parameter reconstruction is crucial from an astrophysical perspective, because it allows both for in-depth studies of individual objects The LIGO Scientific Collaboration et al. (2018a, b); Abbott et al. (2019); The LIGO Scientific Collaboration et al. (2019); Isi et al. (2019); Abbott et al. (2017) and of populations masses The LIGO Scientific Collaboration et al. (2018c); Kapadia et al. (2019); Talbot and Thrane (2018).

From a statistical point of view, Bayesian inference on a population of events with imperfect measurements has a well established formalism Mandel et al. (2019); Hogg et al. (2010). A residual freedom remains in the choice of parameterization for the population. Previous studies have used astrophysically motivated functional dependencies The LIGO Scientific Collaboration et al. (2018c); Talbot and Thrane (2018); Kapadia et al. (2019); Wysocki (2017); Wysocki et al. (2018); Stevenson et al. (2017). For example, one parameter in such models might be the location of a mass gap in the black hole population Roulet and Zaldarriaga (2019); Bai et al. (2018). Other studies have used a broader family of somewhat non-parametric models Del Pozzo et al. (2018); Powell et al. (2019); Keitel (2019); Hjort et al. (2010); Orbanz and Teh (2010).

Within the latter formalism, greater flexibility can be achieved by fitting the observed data with an unknown number of sub-components. No a priori physical meaning is necessarily associated with these components, and they are usually sampled from a common hyper-parameter space. The lack of any hierarchy among these components naturally introduces a symmetry under permutations and leads to the label switching problem.

Here we apply our hypertriangle approach to inference on the population of observed BBH component masses, $m_{1}\geq m_{2}$. We model the observed distribution of source frame Abbott et al. (2016) component black hole masses (in solar mass units) as a mixture $p_{\mathrm{pop}}(\log m_{1},\log m_{2})$ of $K$ bivariate Gaussians;

[TABLE]

Each component has a pair of means, $\mu^{(\log m_{1})}_{k}$ and $\mu^{(\log m_{2})}_{k}$, a symmetric $2\times 2$ covariance matrix, $\mathbf{\Sigma}_{k}$, and a weight, $w_{k}$. The covariance matrix is described by its two eigenvalues, $\lambda^{1}_{k}$ and $\lambda^{2}_{k}$, and a rotation angle $\phi_{k}$. Overall, each component is fully described by the parameter vector

[TABLE]

We choose to enforce the artificial identifiability constraint $\mu_{k+1}^{(\log m_{1})}\geq\mu_{k}^{(\log m_{1})}$. This is done by applying our map $\phi$ from Eq. 5 to the vector of components $\mu_{k}^{(m_{1})}$ with $k=1,2,\ldots,K$. We can sample on the modified parameter space covered by

[TABLE]

where $\mu^{(m_{1})}_{k}=\phi(\chi_{k})$. In the language of Sec. II, sampling on the parameter space in Eq. 14 covers $\mathcal{C}$ (with multimodality) while sampling on Eq. 15 covers $\mathcal{T}$.

The priors are taken to be flat on all of the components in Eqs. 14 and 15, except for the $\lambda^{1}_{k},\lambda^{2}_{k}$ which we take log-uniformly distributed within their ranges.

The ranges for $\chi_{k}$, $\mu_{k}^{(\log m_{1})}$, $\mu_{k}^{(\log m_{2})}$ are $(0,2)$, with the additional constraint of $\mu_{k}^{(\log m_{1})}>\mu_{k}^{(\log m_{2})}$. The range on the angle $\phi_{k}$ is $(0,\pi/2)$ and the ranges on $\lambda_{k}^{1}$ and $\lambda_{k}^{2}$ are $(0.01,4)$. Finally, the weights $w_{k}$ were sampled in the range $(0,1)$ and then normalized such that $\sum_{k}w_{k}=1$.

We adopt a fully Bayesian hierarchical approach. At the lowest level there are the short segments of time series data $\{d\}$ surrounding each of the $N_{\text{obs}}$ events. Each event is described by some parameters $\theta$ (e.g. masses, spins, etc). The likelihood that we wish to sample from is the probability of all the observed data given a certain value of the population parameters $\boldsymbol{\Lambda}=\{\boldsymbol{\Lambda}_{k}|k=1,2,\ldots,K\}$:

[TABLE]

Using Bayes theorem, the above likelihood can be turned into a posterior on the population parameters $\mathbf{\Lambda}$. This in turn can be expressed in terms of the $N_{i}$ posterior samples on $\theta$ from each individual event Mandel et al. (2019):

[TABLE]

where the posterior samples for each event are denoted $\theta_{i}^{j}$ ($i$ labels the event and $j$ labels the sample in the posterior chain), and $\varpi(\boldsymbol{\Lambda})$ and $\pi(\theta)$ respectively denote the priors on the population and individual event parameters. We will consider only the component masses as event parameters, $\theta=(m_{1},m_{2})$. Note that the normalization integral in the denominator of Eq. 17 is evaluated over the constrained prior range $\log m_{1}>\log m_{2}$. We use the publicly available posterior samples Vallisneri et al. (2015) for the 10 BBH events described in The LIGO Scientific Collaboration et al. (2018a).

As our focus here is on the label switching problem, and its solution using the hypertriangle map, for simplicity we do not consider selection effects Foreman-Mackey et al. (2014); Farr et al. (2014). Rather, we model the distribution of observed black hole masses. We defer a full treatment, including selection effects, to future work.

We model the observed mass distribution using $K=1,\dots,4$ Gaussian components. We sample the distribution in Eq. 17 using the nested sampling algorithm Skilling (2004) as implemented in CPNest Veitch et al. (2017). The primary output of the algorithm is the model evidence, which we use to determine which $K$ is favored; we find that the data favors a description using a 1-component Gaussian mixture. Additionally, the algorithm produces samples from the posterior in Eq. 17. The log-evidences for different $K$ are presented in Table 1, while the posterior samples for $K=2$ ($K=3$) on the $\vec{\mu}^{(\log m_{1})}$ parameters are shown in Fig. 2 (Fig. S1 in Suppl. S1). The full posterior chain on all of the parameters is provided at Buscicchio et al. (2019b).

Because this is a relatively low-dimensional problem (we consider $K\leq 4$ ) the analysis can be performed both with and without the hypertriangle map. If the map is not used then the posterior has $K!$ degenerate modes. If the map is used then there is just a single mode and, importantly, no information is lost. The elimination of the excess multimodality is shown for two dimensions in Fig.2. More impressive demonstrations of the elimination of the excess multimodality are possible in higher numbers of dimensions; a plot in 3 dimensions for $K=3$ component mixture is shown in the supplementary material. The preservation of information is demonstrated by the fact that the evidence in unchanged. That the evidence is unchanged can be shown analytically and is a consequence of the Jacobian for our transformation in Eq. 10; it is also demonstrated numerically for this specific problem in Table 1.

We can now use the posteriors on $\mathbf{\Lambda}$ to plot the observed black hole mass distribution. This can be done using the posterior on the $\mathbf{\Lambda}$ from either of Eqs. 14 or 15 with identical results. Although, the single mode posterior from Eq. 15 is naturally easier to sample from. The marginalised mass distributions on $m_{1}$ and $m_{2}$ are plotted in Fig. 3. As shown by the evidences in Table 1 a one component mass distribution is favoured. We stress again that we have not included selection effects; including these is expected to suppress the high mass tail (this is because high mass BBHs can be seen out to greater distances than lower mass systems) and therefore our results are not incompatible with the presence of a mass gap.

The hypertriangle map has demonstrated its utility. It eliminated the excess multimodality in the description of the observed BBH mass distribution. This renders the target posterior easier to sample. There is no loss in information incurred by sampling this remapped parameter space compared to sampling the full original space.

IV.2 Overlapping Galactic White-Dwarf

Binaries in LISA

LISA Amaro-Seoane et al. (2017) is a planned space-based mission which will observe GWs in the (0.1–100) mHz frequency range. The LISA band is source-rich, with many signals overlapping in both time and frequency. In particular, galactic white dwarf binaries (GBs) Cornish and Robson (2017) are so numerous at low frequencies that they form a confusion noise foreground for LISA. Several GBs have already been identified electromagnetically and will serve as verification sources for LISA Stroeer and Vecchio (2006).

The label-switching problem arises in the analysis of multiple sources, since the parameters of any pair of sources are interchangeable. In this section we will show how the application of the hypertriangle map allows for efficient Bayesian recovery of multiple GB signals without ambiguity arising from label switching.

The GWs emitted by a distant source are observed in the solar system as plane waves. There are two GW polarization components denoted $+$ and $\times$. Under the assumption that each source is monochromatic, these components are given by

[TABLE]

where $f$ is the GW frequency, $\iota$ is the inclination angle between the binary’s orbital angular momentum and the line of sight, and $\Phi$ is a phase offset.

The LISA detector response additionally depends on the ecliptic longitude and latitude $\{\lambda,\beta\}$ of the source and a polarization angle $\psi$. The GW amplitude $A$ can be further expressed in terms of physical quantities of the GB system (e.g. the component masses and the luminosity distance); however, these quantities are highly degenerate and are therefore not considered.

Each of the $K$ sources is described by seven parameters:

[TABLE]

We use flat priors on all parameters with ranges given in Sec. S2. The log-likelihood is given by

[TABLE]

where $k$ labels the various GBs, and where $s_{\alpha}$ denotes two approximately independent LISA output channels, with $\alpha\in\{\mathrm{A},\mathrm{E}\}$ (see, for example, Tinto and Dhurandhar (2014)). The model $h_{\alpha}$ is the LISA response to sinusoidal signals of the form in Eq. 18. The line brackets indicate a norm with respect to the usual signal inner product

[TABLE]

where $\tilde{a}(f)$ is the Fourier transform of $a(t)$. Each output channel is assumed to contain additive stationary Gaussian noise with a one-sided power spectral density $S_{\alpha}(f)$.

We simulate one year of mock LISA noise using LISA code Petiteau et al. (2008). For simplicity, we estimate the power spectral densities from these signal-free noise realizations using the Welch periodogram D. Welch (1967); Bartlett (1950).

We inject $K=6$ sources, each with a signal-to-noise ratio $\rho_{k}=10$ defined with respect to the inner product in Eq. 21. The six sources were chosen to have regularly spaced frequencies; other source parameters were chosen randomly and are given in Table 2. For simplicity, we perform a noise-free analysis.

The simulated data has a cadence of $5\,\mathrm{s}$ and a total duration of $1\,\mathrm{yr}$, resulting in arrays of length $6.3\times 10^{6}$. This data was heterodyned, filtered, and downsampled to isolate a narrow range of frequencies $f\in(f_{\star}-1/\mathrm{yr},f_{\star}+11/\mathrm{yr})$, where $f_{\star}=1$ mHz. For a one-year observation period, the expected frequency resolution of LISA is $\sim 1/1\,\text{yr}$, so this frequency range covers 12 bins.

We assume the number of sources $K$ is already known by other means; we do not address the problem of searching for an unknown number of sources (see, for example, Stroeer and Veitch (2009); Crowder and Cornish (2007); Cornish and Larson (2003)).

This is a $6\times 7=42$-dimensional Bayesian inference problem. The likelihood in Eq. (20) is invariant under permutations of the index $k$ (i.e. relabelling the GBs numbered 1 to 6). Naively sampling this distribution in the specified prior ranges will return a posterior distribution with (at least) $6!=720$ peaks. To remove this problem we enforce the artificial identifiability constraint $f_{k+1}\geq f_{k}$ by sampling on the parameters

[TABLE]

Here $f_{k}=\phi(\chi_{k})$ (see Eq. 5), and the prior on $\chi_{k}$ is the same as the prior on $f_{k}$. The sampler explores the space of $\chi_{k}\in\mathcal{C}$ which is mapped to the physical frequencies $f_{k}\in\mathcal{T}$.

The resultant distribution has a single global maximum and is therefore relatively easy to sample from (albeit in 42 dimensions). We use CPNest Veitch et al. (2017) to sample the distribution and correctly recover all sources. We note that without applying our hypertriangle map, it would be excessively difficult to sample from this 720-fold degenerate distribution.

The 1D marginalized posteriors on the physical frequencies $f_{k}$ are plotted in Fig. 4. Additional plots and the full posterior file are provided in the supplementary material.

V Discussion

We discussed a general solution to the label switching problem which allows the sampler to be treated as a black box, and is therefore widely applicable. To enforce the identifiability constraint, we map the sampled points from a hypercube with the desired prior to a hypertriangle, taking care to preserve the prior. We have successfully used this for two real-world problems from gravitational wave astrophysics. The hypertriangle transformation has the potential to greatly simplify a wide class of highly-degenerate Bayesian inference problems, with no loss of information.

Acknowledgements.

The authors gratefully acknowledge the contributions of Alberto Vecchio, Sebastian Gaebel, Antoine Klein, Davide Gerosa, Alberto Sesana, Christopher Berry and Siyuan Chen in the development of the Birmingham LISA data analysis code which was used for the GB problem; full details of this code will be published elsewhere. The authors also thank Alvin Chua, John Veitch and Walter Del Pozzo for discussions about the method and RB acknowledges Gianmarco Brocchi for useful discussions. The computations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service.

S1 Additional LIGO/Virgo population Data and Plots

In this section, we show prior and posteriors for parameters of the observed populations described in Sec. IV.1.

Priors on mixture parameters are given in Table S1. Posteriors on mixture primary log-mass means are shown in Fig. S1. Median a posteriori values of $p(\log m_{1},\log m_{2})$ are shown in Fig. S2 for one and two mixture components.

S2 Additional GB Data and Plots

The full posterior parameters are in Buscicchio et al. (2019b). A selection of these parameters are plotted in Fig. S3