Label Switching Problem in Bayesian Analysis for Gravitational Wave Astronomy
Riccardo Buscicchio Elinore Roebber Janna M. Goldstein Christopher J. Moore
astro-ph.IM


Label Switching Problem]The Label Switching Problem in Bayesian Analysis for Gravitational Wave Astronomy


riccardo@star.sr.bham.ac.uk















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.

[
    Christopher J. Moore
    December 30, 2023
========================



§ 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≤ 1.
This covers only the restricted portion of the parameter space m_1 ≥ 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 <cit.>.
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 ≥ 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, <cit.>.
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 <cit.>.
Alternatively, the log-prior distribution can be crudely modified so that it returns -∞ 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 <cit.>).
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 <cit.> (see, in particular, Eq. A13; however, this equation contains a typographical error as pointed out by <cit.>).
This solution to the label switching problem has been implemented in <cit.> and has been widely used extensively in the astronomical and cosmological literature <cit.>.
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. <ref> 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. <ref>.
In Sec. <ref> the efficacy of our proposed solution is demonstrated by applying it to the second (Sec. <ref>) and third (Sec. <ref>) 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.



§ THE LABEL SWITCHING PROBLEM


We wish to treat problems containing multiple indistinguishable components.
Each of the K components is modeled by some parameters Λ_k∈ U, where the parameter space U is an open set of ℝ^n and k∈{1,2,…,K}. 
We will choose to distinguish components based on the values of one of these parameters, x_k≡Λ_k^1∈ I where I is an set of ℝ.
For simplicity, in this section we will consider x_k ∈(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 × U.
However we will mainly be interested in the subspace spanned by x⃗=(x_1, x_2), which covers the unit square I× I (in the general K-dimensional case this will be a hypercube which will be denoted 𝒞).
For the remainder of this section we suppress the other components of Λ_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, ℒ(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 Λ_1 and Λ_2 simultaneously). 
As a result, the parameter space is twice as large as it needs to be.
Evaluating ℒ(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≥ 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≥ x_1 (see the off principal diagonal panels in Fig. <ref>).  In higher dimensions, the parameter space is restricted to the region x_K≥ x_K-1≥⋯≥ x_1, which is hereafter referred to as the hypertriangle and denoted 𝒯.

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 <cit.>):


    ϕ:𝒞→𝒯 .


Naively, we might try to choose ϕ 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 x⃗=(x_1,x_2) in the hypercube and then the user applies the sorting function x⃗'=sort(x⃗) before evaluating the likelihood ℒ(x⃗'), 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, ϕ, 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: 

    x'_1   = x_1
     
        x'_2   = x'_1+(1-x'_1)x_2 .

To see that points are indeed mapped to 𝒯 it is sufficient to note that the correct ordering is enforced by adding a positive quantity to x'_1 to get a larger value for x'_2. 
The range x'_2∈(0, 1) is in turn ensured by scaling x_2 with the factor 1 - x'_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. <ref>).



The map in Eq. <ref> can be “fixed” by revising the x_1 coordinate downwards, before shifting/rescaling x_2:

    x'_1   = 1-√(1-x_1)
    
        x'_2   = x'_1+(1-x'_1)x_2 .

The new map in Eq. <ref> solves the problem in 2 dimensions (see the blue distribution in Fig. <ref>). 
The sampler can propose points x⃗=(x_1,x_2) uniformly in the square. These points are mapped to the triangle x⃗'=ϕ(x⃗). Finally, the likelihood is evaluated at the mapped points, ℒ(x⃗'). 
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 x⃗'=ϕ(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. [For an extension of our solution to the wider class of separable priors, see Sec. <ref>.]


To state the problem formally: to solve the label-switching problem, we seek a bijection (a “one-to-one” and “onto” map) ϕ:𝒞→𝒯, for an arbitrary number of dimensions, with components x'_κ=ϕ_κ(x_k), such that the Jacobian J≡|∂ x'_κ / ∂ x_k|=const.


§ THE HYPERTRIANGLE MAP IN ARBITRARY DIMENSIONS

Our proposed generalization of the 2-dimensional map in Eq.<ref>,  x'=ϕ(x),  is defined recursively as

    x'_i = x'_i-1+(1-x'_i-1)[1-(1-x_i)^1/K+1-i]  ,

where i∈1,…,K and x'_0= 0 by definition. This closely resembles Equation (A13) of <cit.>, although here we have corrected a typographical error. Eq. <ref> can be expressed non-recursively as:

    x'_i = 1 - ∏_j=1^i (1 - x_j)^1/K+1-j .

If the inputs are in the correct range x_j∈(0,1), i.e. x_j∈𝒞, it can be shown that the output falls in 𝒯 (the logic as outlined in Sec. <ref> for Eq. <ref> still applies). It can also be shown that this map is a bijection by inverting Eq. <ref>.


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



 §.§ Equivalent Representations of ϕ


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

    x_i^'    =  x_i-1^'+(1-x_i-1^')[1-(1-x_i)^1/K+1-i]
        =  [1-(1-x_i)^1/K+1-i]-x_i-1^'[1-(1-x_i)^1/K+1-i-1]
        =  [1-(1-x_i)^1/K+1-i]+(1-x_i)^1/K+1-ix_i-1^'
        =  1-(1-x_i)^1/K+1-i(1-x_i-1^') .

This procedure can be repeated for the x_i-1^' term inside the final set of parentheses, and then for x_i-2^' and so on down to x^'_1.
This gives the equivalent representation to Eq. <ref>;

    x_i^'    =  1-(1-x_i)^1/K+1-i[(1-x_i-1)^1/K+1-(i-1)(1-x_i-2^')]
        = …
        =  1-∏_j=1^i(1-x_j)^1/K+1-j .

 



 §.§ The Jacobian of ϕ


As discussed in Sec. <ref>, to maintain the correct prior on the hypertriangle, it is necessary that the map ϕ has a constant Jacobian.  
To prove that our proposed hypertriangle map has this property, we start with the form of the map in Eq.<ref>.
The Jacobian matrix for this specific transformation is lower-triangular because the component x'_i depends only on x_j with j≤ i. Its determinant is therefore equal to the product of the diagonal terms:

    J     =    ∏_i=1^K∂ x_i^'/∂ x_i
        =    ∏_i=1^K1/K+1-i(1-x_i)^1/K+1-i-1∏_j=1^i-1(1-x_j)^1/K+1-j
        =    1/K!∏_i=1^K1/(1-x_i)(1-x_i)^1/K+1-i∏_j=1^i-1(1-x_j)^1/K+1-j
        =    1/K!∏_i=1^K1/(1-x_i)∏_j=1^i(1-x_j)^1/K+1-j ,

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 

    J   = 1/K!
        ×1/(1-x_1)[(1-x_1)^1/K+1-1]
        ×1/(1-x_2)[(1-x_1)^1/K+1-1(1-x_2)^1/K+1-2]
        ×…
        ×1/(1-x_K)[(1-x_1)^1/K+1-1…(1-x_K)^1/K+1-K]

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

    J = 1/K! .

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


 §.§ 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 

    Π(x_1,…,x_K)=∏_k=1^Kπ(x_k)

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

    F(x)=∫_0^xπ(s)ds .

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.



 §.§ Implementation of ϕ



For concreteness, we provide here a pseudo-Python implementation of Eq. <ref>.
The input  (in 𝒞) and output  (in 𝒯) are numpy arrays where all values are in the prior range (0,1). 
The values of  may be in any order whilst the values of  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 <cit.>. 



§ 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. <ref> is a Gaussian mixture model; models of this type have been studied extensively in the context of the label switching problem <cit.>. 

The second example in Sec. <ref> 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 <cit.> discuss it when fitting multiple damped sinusoids to time series data. 



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


LIGO and Virgo <cit.>
are ground-based GW detectors operating in the (10^1-10^4) 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 <cit.>. The third observation run is ongoing and low latency pipelines <cit.> have produced a number of public alerts associated with event candidates <cit.>.
It is likely that by the end of the current run dozens more detections will be available <cit.>
for further investigation.

Detailed waveform models for BBH signal calibrated against numerical relativity are now available  <cit.>.
These are used in the LALInference Bayesian analysis software package <cit.> 
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 <cit.> and of populations masses <cit.>.

From a statistical point of view, Bayesian inference on a population of events with imperfect measurements has a well established formalism <cit.>.
A residual freedom remains in the choice of parameterization for the population. 
Previous studies have used astrophysically motivated functional dependencies <cit.>.
For example, one parameter in such models might be the location of a mass gap in the black hole population <cit.>.
Other studies have used a broader family of somewhat non-parametric models <cit.>. 

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≥ m_2. 
We model the observed distribution of source frame <cit.> component black hole masses (in solar mass units) as a mixture p_pop(log m_1,log m_2) of K bivariate Gaussians; 

    [ log m_1; log m_2 ]∼∑^K_k=1 w_k 𝒩(
               [ μ_k^(log m_1); μ_k^(log m_2) ],
               Σ_k
              )  .

Each component has a pair of means, μ^(log m_1)_k and μ^(log m_2)_k, a symmetric 2× 2 covariance matrix, Σ_k, and a weight, w_k. 
The covariance matrix is described by its two eigenvalues, λ^1_k and λ^2_k, and a rotation angle ϕ_k. 
Overall, each component is fully described by the parameter vector

    Λ_k=(μ^(log m_1)_k,μ^(log m_2)_k, λ^1_k, λ^2_k, ϕ_k, w_k ) .


We choose to enforce the artificial identifiability constraint μ_k+1^(log m_1)≥μ_k^(log m_1). 
This is done by applying our map ϕ from Eq. <ref> to the vector of components μ_k^(m_1) with k=1,2,…,K.
We can sample on the modified parameter space covered by

    Λ_k=(χ_k,μ^(log m_2)_k, λ^1_k, λ^2_k, ϕ_k, w_k ) .

where μ^(m_1)_k=ϕ(χ_k). 
In the language of Sec. <ref>, sampling on the parameter space in Eq. <ref> covers 𝒞 (with multimodality) while sampling on Eq. <ref> covers 𝒯.

The priors are taken to be flat on all of the components in Eqs. <ref> and <ref>, except for the λ^1_k,λ^2_k which we take log-uniformly distributed within their ranges.

The ranges for χ_k, μ_k^(log m_1), μ_k^(log m_2) are (0,2), with the additional constraint of μ_k^(log m_1)>μ_k^(log m_2).
The range on the angle ϕ_k is (0,π/2) and the ranges on λ_k^1 and λ_k^2 are (0.01, 4).
Finally, the weights w_k were sampled in the range (0,1) and then normalized such that ∑_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_obs events.
Each event is described by some parameters θ (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 Λ={Λ_k|k=1,2,…,K}:

    p({d}|Λ)= ∏_i=1^N_obs∫dθ  p(d |θ) p_pop(θ|Λ)/∫dθ  p_pop(θ|Λ) .

Using Bayes theorem, the above likelihood can be turned into a posterior on the population parameters Λ.
This in turn can be expressed in terms of the N_i posterior samples on θ from each individual event <cit.>:

    p(Λ|{d}) = ϖ(Λ)∏_i=1^N_obs1/N_i∑_j=1^N_ip_pop(θ_i^j |Λ)π(θ_i^j)∫dθ  p_pop(θ|Λ),

where the posterior samples for each event are denoted θ_i^j (i labels the event and j labels the sample in the posterior chain), and ϖ(Λ) and π(θ) respectively denote the priors on the population and individual event parameters. 
We will consider only the component masses as event parameters, θ=(m_1,m_2). 
Note that the normalization integral in the denominator of Eq. <ref> is evaluated over the constrained prior range log m_1 > log m_2.
We use the publicly available posterior samples <cit.> for the 10 BBH events described in <cit.>. 



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 <cit.>.
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,…,4 Gaussian components.
We sample the distribution in Eq. <ref> using the nested sampling algorithm <cit.> 
as implemented in CPNest <cit.>. 
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. <ref>.
The log-evidences for different K are presented in  Table <ref>, while the posterior samples for K=2 (K=3) on the μ⃗^(log m_1) parameters are shown in Fig. <ref> (Fig. <ref> in Suppl. <ref>).
The full posterior chain on all of the parameters is provided at <cit.>.

Because this is a relatively low-dimensional problem (we consider K≤ 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.<ref>.
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. <ref>; it is also demonstrated numerically for this specific problem in Table <ref>.

We can now use the posteriors on Λ to plot the observed black hole mass distribution. 
This can be done using the posterior on the Λ from either of Eqs. <ref> or <ref> with identical results. Although, the single mode posterior from Eq. <ref> is naturally easier to sample from.
The marginalised mass distributions on m_1 and m_2 are plotted in Fig. <ref>.
As shown by the evidences in Table <ref> 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.


    


 §.§ Overlapping Galactic White-Dwarf 
 Binaries in LISA



LISA <cit.> 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) <cit.> 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 <cit.>.

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 ×.
Under the assumption that each source is monochromatic, these components are given by

    h_+(t; Λ)         = A (1+cos^2ι) cos(2π ft-Φ)  ,
    
        h_×(t; Λ)    = -2A cosι sin(2π ft-Φ)  ,

where f is the GW frequency, ι is the inclination angle between the binary's orbital angular momentum and the line of sight, and Φ is a phase offset. 

The LISA detector response additionally depends on the ecliptic longitude and latitude {λ, β} of the source and a polarization angle ψ.
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:

    Λ_k = {log_10 A_k, f_k, λ_k, sinβ_k, cosι_k, ψ_k, Φ_k} .

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

    logℒ(Λ_k) ∝ -1/2∑_α| s_α - ∑_k=1^K h_α(Λ_k)|_(α)^2 ,

where k labels the various GBs, and where s_α denotes two approximately independent LISA output channels, with α∈{A,E} (see, for example, <cit.>).
The model h_α is the LISA response to sinusoidal signals of the form in Eq. <ref>.
The line brackets indicate a norm with respect to the usual signal inner product

    <a|b>_(α) = 4{∫_0^∞df ã(f)b̃(f)/S_α(f)}  ,

where ã(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_α(f).

    

We simulate one year of mock LISA noise using LISA code <cit.>. For simplicity, we estimate the power spectral densities from these signal-free noise realizations using the Welch periodogram <cit.>.



We inject K=6 sources, each with a signal-to-noise ratio ρ_k=10 defined with respect to the inner product in Eq. <ref>.
The six sources were chosen to have regularly spaced frequencies;
other source parameters were chosen randomly and are given in Table <ref>.
For simplicity, we perform a noise-free analysis.

The simulated data has a cadence of 5 s and a total duration of 1 yr, resulting in arrays of length 6.3× 10^6.
This data was heterodyned, filtered, and downsampled to isolate a narrow range of frequencies f∈(f_⋆ -1/yr, f_⋆ +11/yr), where f_⋆  = 1 mHz. 
For a one-year observation period, the expected frequency resolution of LISA is ∼ 1/1 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, <cit.>).


This is a 6×7=42-dimensional Bayesian inference problem.
The likelihood in Eq. (<ref>) 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≥ f_k by sampling on the parameters 

    Λ_k = {log_10 A_k, χ_k, cosι_k, λ_k, sinβ_k, ψ_k, ϕ_k}  .

Here f_k = ϕ(χ_k) (see Eq. <ref>), and the prior on χ_k is the same as the prior on f_k.
The sampler explores the space of χ_k ∈𝒞 which is mapped to the physical frequencies f_k ∈𝒯.

The resultant distribution has a single global maximum and is therefore relatively easy to sample from (albeit in 42 dimensions).  
We use CPNest <cit.> 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. <ref>. Additional plots and the full posterior file are provided in the supplementary material.  


§ 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.





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.











  
  
   
   










SUPPLEMENTARY MATERIAL
 





§ ADDITIONAL LIGO/VIRGO POPULATION DATA AND PLOTS
 

In this section, we show prior and posteriors for parameters of the observed populations described in Sec. <ref>. 

Priors on mixture parameters are given in Table <ref>. Posteriors on mixture primary log-mass means are shown in Fig. <ref>. Median a posteriori values of p(log m_1,log m_2) are shown in Fig. <ref> for one and two mixture components.








§ ADDITIONAL GB DATA AND PLOTS


The full posterior parameters are in <cit.>. A selection of these parameters are plotted in Fig. <ref>



