Gravitational wave detectors like the Einstein Telescope and LISA generate long multivariate time series which pose significant challenges in spectral density estimation due to a number of overlapping signals as well as the presence of correlated noise. Addressing both issues is crucial for accurately interpreting the signals detected by these instruments. This paper presents an application of a variational inference spectral density estimation method specifically tailored for dealing with correlated noise in the data. It is flexible in that it does not rely on any specific parametric form for the multivariate spectral density. The method employs a blocked Whittle likelihood approximation for stationary time series and utilizes the Cholesky decomposition of the inverse spectral density matrix to ensure a positive definite estimator. A discounted regularized horseshoe prior is applied to the spline coefficients of each Cholesky factor and the posterior distribution is computed using a stochastic gradient variational Bayes approach. This method is particularly effective in addressing correlated noise a significant challenge in the analysis of multivariate data from colocated detectors. The method is demonstrated by analyzing 2000 seconds of simulated Einstein Telescope noise which shows its ability to produce accurate spectral density estimates and quantify coherence between time series components. This makes it a powerful tool for analyzing correlated noise in gravitational wave data.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Statistics - Applications']
Variational inference for correlated gravitational wave detector network noise
{'Variational inference for correlated gravitational wave detector network noise': "Jianan Liu 1 , Avi Vajpeyi 1 , Renate Meyer 1 , Kamiel Janssens 2 , 3 , 4 , 5 , Jeung Eun Lee 1 , Patricio Maturana-Russel 1 , 6 , Nelson Christensen 3 , and Yixuan Liu 7 1 Department of Statistics, University of Auckland, 38 Princes St, Auckland, New Zealand 2 Universiteit Antwerpen, Prinsstraat 13, 2000 Antwerpen, Belgium 3 Universit'e Cˆote d'Azur, Observatoire Cˆote d'Azur, CNRS, Artemis, 06304 Nice, France 4 Department of Physics, The University of Adelaide, Adelaide, SA 5005, Australia 5 ARC Centre of Excellence for Dark Matter Particle Physics, Australia 6 Department of Mathematical Sciences, Auckland University of Technology, Auckland, New Zealand and 7 MRC Biostatistics Unit, University of Cambridge, Cambridge, UK  \n(Dated: 2024-09-23)  \nGravitational wave detectors like the Einstein Telescope and LISA generate long multivariate time series, which pose significant challenges in spectral density estimation due to a number of overlapping signals as well as the presence of correlated noise. Addressing both issues is crucial for accurately interpreting the signals detected by these instruments. This paper presents an application of a variational inference spectral density estimation method specifically tailored for dealing with correlated noise in the data. It is flexible in that it does not rely on any specific parametric form for the multivariate spectral density. The method employs a blocked Whittle likelihood approximation for stationary time series and utilizes the Cholesky decomposition of the inverse spectral density matrix to ensure a positive definite estimator. A discounted regularized horseshoe prior is applied to the spline coefficients of each Cholesky factor, and the posterior distribution is computed using a stochastic gradient variational Bayes approach. This method is particularly effective in addressing correlated noise, a significant challenge in the analysis of multivariate data from co-located detectors. The method is demonstrated by analyzing 2000 s of simulated Einstein Telescope noise, which shows its ability to produce accurate spectral density estimates and quantify coherence between time series components. This makes it a powerful tool for analyzing correlated noise in gravitational wave data.", 'I. INTRODUCTION': 'The next-generation gravitational-wave (GW) detectors, e.g. Einstein Telescope (ET) [1], Cosmic Explorer (CE) [2], and Laser Interferometer Space Antenna (LISA) [3], will usher in a transformative era of GW astronomy. The increased sensitivity to GW will lead to a much higher GW detection rate, improving our understanding of astrophysical and cosmic phenomena, ranging from the dynamics of black hole mergers to the nature of the early universe [2-6].  \nHowever, multivariate time series analysis of future GW detector data faces challenges in managing correlated multivariate noise [7, 8]. ET and LISA, which consist of co-located interferometers, may experience correlated noise in their data streams [9]. For ET, this correlated noise can stem from seismic, Newtonian, and magnetic sources [10-13], while LISA may encounter noise from temperature variations and micro-thrusters [14, 15].  \nIgnoring correlated noise in analyses of multivariate GWdata can lead to biased parameter estimates and incorrect astrophysical conclusions (e.g., in analyzing the stochastic GW background [16-18] and transient signals from events like binary black hole mergers [19]). Recent studies have shown the need for methods that can handle correlated noise [9], but did not include a demonstration of Bayesian methods capable of effectively doing so. Other work, [19] did demonstrate the capability to estimate the physical parameters of a compact binary  \ncoalescence event in a detector network with correlated noise. Their likelihood included the multivariate spectral density but they treated the spectral density as fixed with a focus on estimating the parameters of the binary chirp signal. While this is an important step forward in the field, the authors of [19] made a number of simplifying assumptions. First of all, they considered the correlated noise to be in phase. Secondly, they considered the detectors to be (almost) maximally correlated across the entire frequency spectrum. So far no Bayesian method has been suggested to estimate multivariate spectral densities of gravitational wave network data including correlated noise.  \nIn this work, we demonstrate a method which is able to account for correlated noise as well as uncorrelated noise in a simultaneous and consistent approach. Additionally we make no assumptions about the phase of the different noise terms and model both real and imaginary cross-terms. While we demonstrate our method on a simplified noise example, the method forms the basis for future work with more complex noise scenarios, but this work shows the effectiveness of the chosen approach. Additionally our proposed method allows for a flexible Bayesian approach for matrix-valued spectral density estimation, in a nonparametric way. Such a method is robust with respect to deviations from certain parametric shapes such as power laws in order to be able to capture all potential small- and large-scale noise characteristics. Furthermore, the approach must exhibit computational efficiency when applied to long time series of GW ob-  \nservations, as in the case of ET and LISA. Meier et al. [20] and Liu et al. [21] have developed non-parametric Bayesian approaches to multivariate spectral density estimation based on the Whittle likelihood and a Bernstein - Hermitian positive semidefinite matrix-Gamma process prior on the spectral density matrix and used an adaptive MCMC algorithm to sample from the posterior. Although theoretically attractive with proven consistency properties and contraction rates, MCMC methods require significant computation time due to the complexity of sampling from high-dimensional posterior distributions, making them impractical for long time series of GW observations, where the computational burden increases dramatically.  \nIn this paper, a stochastic gradient variational Bayes (SGVB) approach, as developed by Hu and Prado [22], is proposed to estimate the multivariate Power Spectral Density (PSD) for correlated ET noise. As a development within the broader framework of variational inference (VI), SGVB optimizes a surrogate posterior distribution by minimizing the Kullback-Leibler (KL) divergence from the true posterior distribution [23-25]. This method transforms complex posterior sampling into an optimization problem, enhancing sampling efficiency and reducing computational demands [26, 27]. Recent advances in VI, such as normalizing flows, have demonstrated its effectiveness for pulsar timing-array datasets [28].  \nIn this article, we demonstrate the usefulness and flexibility of a SGVB approach to estimate the spectral density matrix of the correlated instrumental noise from a network of GW detectors. In particular, the approach is illustrated by using simulated ET noise and is also readily applicable to estimate LISA instrumental noise or correlated LIGO-VIRGO-KAGRA (LVK) noise induced by magnetic noise such as from the Schumann resonances [16]. To ease computational efforts of handling large time series, a blocked Whittle likelihood approximation is introduced and utilized instead of the traditional Whittle likelihood. Advice is provided on choosing tuning parameters such as the learning rate and number of basis functions in the SGVB approach. Finally, as a by-product of this method, the squared coherence is computed as a function of frequency, quantifying the degree of inter-temporal correlation between components of the multivariate time series.  \nThe article is structured as follows. Section II introduces the stochastic gradient variational Bayes (SGVB) method and defines the blocked Whittle likelihood. The efficiency and accuracy of the SGVB approach are tested using simulations in Section III. Section IV presents results of the method applied to simulated ET datasets consisting of varying levels of correlated noise. Lastly, the article concludes with a summary of findings in Section V.', 'A. Likelihood': "Assume we are given a time series of gravitational wave observations of length n from p channels. To set notation, let Z = ( Z 1 , . . . , Z n ) ⊺ ∈ R n × p be a p -dimensional stationary, mean-zero time series, sampled at time intervals ∆ t = 1 / (2 f Ny ), where f Ny is the Nyquist frequency, for a total observation time of T with a total of n = T/ ∆ t sampled values for each of the p dimensions. The frequency resolution, ∆ f , is given by  \n∆ f = 1 n ∆ t = 1 T . (1)  \nLet d k be the Discrete Fourier Transform (DFT) of Z given by  \nd ( f k ) = √ ∆ t n n ∑ t =1 Z t exp ( -2 πi k n t ) , (2)  \nwhere f k = k ∆ f = k/ ( n ∆ t ) = k/T for k = 1 , . . . , N , where N = n/ 2 if n is even, N = ( n -1) / 2 if n is odd. In this analysis, we use the two-sided Fourier transform, which includes both positive and negative frequencies, allowing for a full representation of the frequency spectrum. However, for the purposes of this study, we only retain the positive frequency components.  \nDue to the normalizing and decorrelating properties of the Fourier transform, the discrete Fourier coefficients d ( f k ) are asymptotically independent and have a complex Gaussian distribution with covariance matrix given by the spectral density matrix S ( f k ), the Fourier transform of the autocovariance function. This asymptotic Gaussian distribution is the basis of the multivariate Whittle likelihood approximation in the frequency domain, expressed as  \nL ( d | S ) ∝ N ∏ k =1 det( S ( f k )) -1 × exp ( -d ( f k ) ∗ S ( f k ) -1 d ( f k ) ) (3)  \nwhere d ( f k ) ∗ is the conjugate transpose of the d ( f k ) and S ( f k ) is a p × p Hermitian positive semidefinite spectral density matrix at each f k . Therefore, the unknown quantity here is S , a matrix-valued function at each frequency with the additional restriction that its value at each frequency is Hermitian positive semidefinite. Note that when estimating the spectral density matrix, it is important that the estimate is again positive semidefinite at each f k so that the quadratic form in the exponent of the Whittle likelihood remains positive and thus defines a valid likelihood.  \nGiven d , one method to estimate the multivariate spectral density matrix S is Bayesian inference. In particular, having a flexible Bayesian method instead of simply using  \na frequentist Welch estimate is important when the ultimate task is to simultaneously estimate the parameters of a GW signal while properly taking the uncertainties of the instrumental noise estimation into account. The Bayes' theorem is given by  \np ( S | d ) = L ( d | S ) π ( S ) Z ( d ) ∝L ( d | S ) π ( S ) , (4)  \nwhere π ( S ) is the prior distribution and p ( S | d ) is the posterior density of unknown S given d , and Z ( d ) is the Bayesian evidence (see Thrane and Talbot [29] and Christensen and Meyer [30] for discussions on GW Bayesian inference).  \nIn GW parameter estimation with large datasets, such as those from ET and LISA, the computational burden of directly applying the Whittle likelihood can be significant. To mitigate this challenge, a 'blocked' Whittle likelihood is adopted [31]. The time series is divided into N b equal-sized blocks Z = ( Z (1) , . . . , Z ( N b ) ) ⊺ with each block a p -dimensional time series of each of length n/N b . The discrete Fourier transform of each block is denoted by d ( i ) , i = 1 , . . . , N b . The stationarity assumption implies that the statistical properties of each block, in particular, their spectral densities are the same. Assuming independence among different blocks, the likelihood then becomes the product of the individual blocked Whittle likelihoods,  \nL b ( d | S ) = N b ∏ i =1 L ( d ( i ) | S ) . (5)  \nNote that as the length of the blocked data is less than the original dataset, the frequency resolution of the spectral density will become coarser as the number of blocks N b increases. In practice, the number of blocks should be chosen to achieve a required frequency resolution while also remaining computationally tractable.", 'B. Parametrization of S': "The prior defined by Hu and Prado [22], Rosen and Stoffer [32] models the components of the Cholesky factorization of the inverse spectral density matrix via smoothing splines (for details, refer to Hu and Prado [22]). In their approach, the inverse of S ( f k ) can be represented as S ( f k ) -1 = T ∗ k D -1 k T k , where D k is a diagonal matrix with diagonal elements δ 2 1 k , δ 2 2 k , ..., δ 2 pk and  \nT k =         1 0 0 · · · 0 -θ ( k ) 21 1 0 · · · 0 -θ ( k ) 31 -θ ( k ) 32 1 . . . . . . . . . . . . . . . . . . 0 -θ ( k ) p 1 -θ ( k ) p 2 · · · -θ ( k ) p,p -1 1         (6)  \nis a p × p complex unit lower triangular matrix. Thus the Whittle likelihood can be rewritten as the product of its p components depending on θ j , δ j , for j = 1 , . . . , p ,  \nL ( d | S ) ∝ p ∏ j =1 L j ( d j | θ j , δ j ) (7)  \nwhere  \nL j ( d j | θ j , δ j ) ∝ N ∏ k =1 δ -2 jk exp    -∣ ∣ ∣ d j ( f k ) -∑ j -1 l =1 θ ( k ) jl d l ( f k ) ∣ ∣ ∣ 2 δ 2 jk    (8)  \nthe parameters θ j , δ j represent the set of θ k and δ k for the j th component of the multivariate time series, θ ( k ) jl represents the corresponding element in the matrix T k for any j > l , and θ ( k ) jl = 1 when j = l , d j ( f k ) denotes the Fourier coefficient of the j th component, and d j = ( d j ( f 1 ) , . . . , d j ( f N )) ⊺ . Then log( δ 2 jk ) and the real and imaginary parts of θ ( k ) jl are modeled by DemmlerReinsch basis functions in terms of M truncated smoothing splines, given by  \nℜ ( θ ( k ) jl ) = α jl, 0 + α jl, 1 f k + M -1 ∑ s =1 ψ ( f k ) α jl,s +1 , (9)  \nℑ ( θ ( k ) jl ) = β jl, 0 + β jl, 1 f k + M -1 ∑ s =1 ψ ( f k ) β jl,s +1 , (10)  \nlog δ 2 jk = γ j, 0 + γ j, 1 f k + M -1 ∑ s =1 ψ ( f k ) γ j,s +1 , (11)  \nwhere ψ ( x ) = √ 2 cos( sπx ), with s being an integer that ranges from 1 to M -1, represents the spline basis function. The flexibility of the model can be adjusted by choosing the number of basis splines M . Thus, a flexible model for the spectral densities S = S ( ν ) depending on a parameter vector ν that comprises all α, β and γ parameters is constructed. Following Hu and Prado [22], discounted regularized horseshoe prior[33] are used for the α, β and γ spline coefficients. This prior is adept at handling varying degrees of smoothness in the individual components of the spectral density matrix while avoiding overfitting (for details, refer to [22, 33]). Thus, in line with Equation 7, the subvector of the parameter vector ν that contains all parameters for component j is denoted as ν j . This allows us to decompose the posterior into the product of p posteriors for each individual component, given by  \np ( ν | d ) = p ∏ j =1 p j ( ν j | d j ) . (12)  \nThis prescription allows the application of the SGVB approach to each p j ( ν j | d j ) in parallel. The approach can be futher parallelised by computing each of the likelihood 'blocks' from Equation 5 independently.", 'C. Stochastic Gradient Variational Bayes': 'This section provides a brief review of SGVB, along with discussions on how to tune the learning rate and the number of basis splines.  \n- a. Variational Inference Review: The fundamental concept of variational inference (VI) is to approximate the posterior distribution p j ( ν j | d j ) by using a surrogate distribution from a family of variational distributions Q j = { q ϕ j ( ν j ); ϕ j ∈ Φ j } , which depends on a parameter vector ϕ j within a parameter space Φ j . Typically, the variational family is selected for its simplicity and computational tractability. A product of normal distributions with mean µ ji and variance σ 2 ji is utilized for each element of the parameter ν j . The goal of the variational approach is to identify the parameters ϕ j = ( µ ji , σ 2 ji ), for i = 1 , . . . , dim( ν j ), that minimize the reverse Kullback-Leibler (KL) divergence between the variational family and the true posterior distribution, denoted as d KL ( q ϕ j || p j )), i.e.,  \nϕ ∗ j = argmin ϕ j d KL ( q ϕ j || p j ) = argmin ϕ j ∫ log q ϕ j ( ν j ) p j ( ν j | d j ) q ϕ j ( ν j ) d ν j . (13)  \nThe optimization algorithm employed in this work uses a SGVB approach, as described by Kingma and Welling [27], Xu et al. [34], Domke [35].  \n- b. Choice of the Number of Basis Functions: The variation of parameters (log δ 2 jk , ℜ ( θ ( k ) j · ) , ℑ ( θ ( k ) j · )) across frequencies f k is modeled using M truncated smoothing splines. Choosing an appropriate M is crucial: a low M may cause underfitting, while a high M can lead to overfitting and increased computational complexity due to the higher dimensionality of ϕ j . To address overfitting concerns, a horseshoe prior on the coefficients acts as a regularization mechanism, effectively reducing the risk of overfitting [36].  \nThe likelihood function (Equation 8) measures the consistency of the data with the spectral density, parameterized by (log δ 2 jk , ℜ ( θ ( k ) j · ) , ℑ ( θ ( k ) j · )). We use the likelihood in the maximum likelihood estimate (MLE) to choose M . The likelihood at MLE tends to increase with M , reflecting an improvement in model fit as the model complexity increases. When M is sufficiently large, the likelihood tends to stabilize or increase steadily. We chose M to prevent underfitting by conservatively selecting a slightly larger M than the point where the likelihood stabilizes.  \n- c. Optimization of the Learning Rate: The KL minimizer ϕ ∗ j (Equation 13) is equivalent to the maximizer of the evidence lower bound (ELBO) between q ϕ j and p j ( ν j | d j ) i.e.,  \nELBO( q ϕ j , p j ( ν j | d j )) = E ν j ∼ q ϕ j ( ν j ) [log p j ( ν j | d j ) -log q ϕ j ( ν j )] . (14)  \nThe SGVB approach is used to estimate ϕ ∗ j , j = 1 , ..., p and it consists of the next two steps:  \n̸  \n- 1. Maximize log p ( ν j | d j ) with respect to ν j .\n- 2. Maximize the ELBO (Equation 14) with respect to ϕ j after substituting the maximized log p ( ν j | d j ) from the first step.  \nThe performance of each maximizer depends on the learning rate; if the learning rate is too small or too large, the optimization algorithm may get stuck at a local maximum or take too long to find the global maximum. Thus, choosing the optimal learning rate is crucial. Recently [37] proposed automated techniques for the blackbox VB in which the learning rate is adaptively decreased. The reliability of the stochastic optimization methods for VB due to substantial hand-tuning parameters to apply effectively were discussed in literature [37, 38].  \nWhile Hu and Prado [22] selected a specific but arbitrary value for the learning rate, this study proposes a method to select the optimal learning rate and automate the selection procedure.  \nLet τ 1 and τ 2 be the learning rates associated with the first and second steps, respectively. As τ 2 has a broad range of plausible values, it is unnecessary to deviate from the default settings proposed by Hu and Prado [22]. An appropriate τ 1 is obtained by maximizing the ELBO  \nτ ∗ 1 = argmax τ 1 p ∏ j =1 ELBO( q ϕ j , p j ( ν j | d j )) . (15)  \nThis optimization is carried out over a continuous parameter space using the Python package HyperOpt [39] and the tree-structured Parzen Estimator (TPE) algorithm [40]. From our preliminary study, the optimal learning rates obtained by fixing τ 2 coincide with those obtained when optimizing τ 1 and τ 2 jointly. The maximum number of iterations for the optimization procedure is set to 10 000.', 'D. Squared Coherence': 'To quantify the frequency-dependent relationship between the channels of the multivariate ET time series, we employ the squared coherence, widely used in many fields [e.g., 41, 42]. Squared coherence, C xy ( f k ), is a normalized measure of association between two time series at frequency f k , ranging from 0 (no correlation) to 1 (perfect correlation). For two components of a multivariate time series, it is defined as:  \nC xy ( f k ) = | S xy ( f k ) | 2 S xx ( f k ) S yy ( f k ) (16)  \n̸  \nwhere S xy ( f k ) is the cross-spectral density between components x and y , and S xx ( f k ), S yy ( f k ) are the spectral densities of components x and y , respectively, at frequency f k ( x, y = 1 , 2 , ..., p ; x = y ), for x, y = 1 , 2 , ..., p, with x = y .  \nFigure 3 displays the speed-up factor of the SGVB compared to the VNPC method. Additionally, Table I  \n<!-- image -->  \nM  \nFIG. 1: The relationship between the number of basis functions ( M ) and the log-normalized likelihood at the MLE for VAR(2) and VMA(1) models (solid and dotted lines respectively) with different data lengths n = 256 (blue), n = 512 (green), and n = 1024 (orange). A fixed learning rate of 0.002 was used across all cases.', 'III. SIMULATION STUDY': "In order to test the efficiency and accuracy of SGVB, 500 independent realisations of a bivariate time series are generated from a vector autoregressive model of order 2 (VAR(2)) and a vector moving-average model of order 1 (VMA(1)), using three different sample sizes n = 256; 512; 1024 (refer to Liu et al. [Section 4.2, 21] for the definitions of the VAR(2) and VMA(1) models). The spectral densities are estimated using both the SGVB method and Liu et al. [21]'s vectorized nonparametrically corrected method (VNPC), an MCMC-based approach that samples from the exact posterior distribution without a variational approximation 1 .  \nIn order to avoid underfitting in SGVB, an appropriate number of basis functions is determined for each dataset using the method discussed in Section II C. Figure 1 shows the log MLE (normalized for comparison between datasets) plotted against the number of basis functions. M = 30 is set as the log MLE enters a steady ascent phase and before any perturbations appear, thereby avoiding both underfitting across the three different dataset sizes. When n = 256, a rapid increase of the log-likelihood is observed as M > 50, indicating overfitting. As shown in equations (9), (10), and (11), the number of parameters is 4( M +2) for a bivariate time series, which exceeds the number of data points.  \nGiven an estimated spectral density matrix ˆ S and the true spectral density matrix S 0 , the accuracy of PSD estimates can be assessed using the L 2 error. The L 2 error  \nr  \no  \nr  \nr  \ne  \nL  \n<!-- image -->  \nn  \nFIG. 2: Violin plots of VNPC and SGVB L 2 errors from 500 realizations for VAR(2) and VMA(1) models, for different data lengths n . The SGVB method uses optimized values for M and τ 1 .  \nSpeed up  \n<!-- image -->  \nn  \nFIG. 3: Violin plots of the speed up gained by using SGVB compared to VNPC from 500 realizations for VAR(2) and VMA(1) models, for different data lengths n . The SGVB method uses optimized values for M and τ 1 .  \nprovides a quantitative measure of the discrepancy between ˆ S and S 0 , encapsulating the overall difference between the estimated and true spectral density matrices across the frequency range. A smaller L 2 error indicates that ˆ S closely approximates S 0 , while a larger L 2 error reflects greater deviation. The L 2 error can be approximated as  \n|| ˆ S -S 0 || L 2 ≈ ( 1 N N ∑ k =1 || ˆ S 0 ( f k ) -S 0 ( f k ) || 2 ) 1 / 2 , (17)  \nwhere ∥·∥ denotes the Frobenius norm, which for a complex p × p matrix A is defined by  \n∥ A ∥ = √ ∑ i,j | a ij | 2 . (18)  \nFigure 2 displays the L 2 error distributions from both the SGVB and VNPC methods across the different sample sizes. The L 2 errors are calculated by using the estimated spectral densities and the true spectral densities of the VAR(2) simulations (top panel) and the VMA(1) simulations (bottom panel), respectively.  \npresents the simulation study's median L 2 errors, pointwise coverage, and median width of pointwise 90% credible regions, as well as the median computation time (in seconds) 2 .  \nFigure 2 reveals that as the sample size increases, both methods achieve lower median L 2 errors, and the median L 2 errors are similar between the two methods. Additionally, Figure 3 shows the median speed-up factor for SGVB becomes increasingly pronounced, ultimately surpassing VNPC by more than 50 times, as n is increased. Notably, when n = 1024, SGVB's accuracy is only marginally lower than VNPC's, yet it reduces median computation time by a factor of approximately 50. The simulation study provides empirical evidence for the consistency of the SGVB approach, demonstrating that as sample size grows, the posterior distributions contracts around the true spectral density, following the same trend as the VNPC approach. However, it does so at a fraction of the computation time.  \nFinally, Table I illustrates that, even though VNPC method consistently outperforms SGVB across all sample sizes in terms of pointwise coverage, typically achieving 83-90% coverage compared to SGVB's 61-66%, SGVB generally produces narrower credible regions. The limited coverage and narrow credible intervals observed from the SGVB results can be attributed to the mean field approximation [25] which assumes independent components ignoring correlations between parameters, and using the KLdivergence as distance measure. The combined effects tend to underestimate the marginal variances of the target distribution [26, 43]. To address these limitations, researchers have explored more flexible surrogate distributions and a broader range of divergence measures in variational inference[25, 44-47].  \nAdditionally, studies have investigated model misspecification and developed assessment methods for approximation inference [48-51]. While more sophisticated frameworks yield improved posterior approximations, they often demand greater computational resources. In this work, we prioritize computational efficiency while maintaining an acceptable level of accuracy.", 'A. Data Generation': 'Three independent realizations of Gaussian noise were synthesized (one for each of the XYZ channels), each spanning 2000 s, spectrally shaped to match the design sensitivity of the ET xylophone configuration [52, 53] 3 .  \nThe noise was sampled at a frequency of 2048 Hz, resulting in a multivariate time series comprising 4096 k data points per channel. In the rest of the paper we refer to this colored Gaussian noise as ET noise .  \nFollowing Janssens et al. [9], non-identical correlated noise in the X, Y, and Z channels is simulated by injecting colored Gaussian noise characterized by frequencydomain Gaussian peaks. The power spectral density is given by  \nS GP n ( f, µ, A ) = ( A √ 2 π exp { -1 2 ( f -µ ) 2 }) 2 , (19)  \nwhere A represents the amplitude, and µ the frequency peak location.  \nGaussian noise is injected at specific frequencies in each channel as follows  \nX : S GP n ( µ =10Hz , A =4 × 10 -24 ) , S GP n ( µ =50Hz , A =2 × 10 -24 ) ,  \nY : S GP n ( µ =10Hz , A =4 × 10 -24 ) , S GP n ( µ =90Hz , A =1 . 5 × 10 -24 ) ,  \nZ : S GP n ( µ =50Hz , A =2 × 10 -24 ) , S GP n ( µ =90Hz , A =1 . 5 × 10 -24 ) .  \nTo introduce correlated noise between channel pairs, we utilize identical Gaussian peaks at matching frequencies across the paired channels. In contrast, for uncorrelated noise scenarios, independent Gaussian peaks are simulated for each channel, ensuring that there is no crosschannel correlation.  \nAs stated in Janssens et al. [9], this dataset is inherently simplified due to the distinct nature of the correlated noise terms, which are highly differentiable. Nevertheless, it serves as a valuable initial demonstration of our approach. Future work aims to investigate more realistic and complex noise scenarios, such as the presence of correlated magnetic Janssens et al. [12, 13] and/or Newtonian noise [10, 11] as well as (multiple) GW signals. However, this is considered beyond the scope of this work where we want to lay the foundation for the estimation framework.', 'B. Data Analysis': "Three analyses are performed to estimate the power spectral density (PSD) and the corresponding squared coherences under the following different noise conditions,  \nCase A: ET noise with correlated Gaussian peaks,  \nCase B: ET noise with uncorrelated Gaussian peaks,  \nCase C: ET noise with uncorrelated Gaussian peaks, assuming independence between channels.  \nM  \nTABLE I: Comparison of median L 2 errors, empirical pointwise coverage, median width of pointwise 90% credible intervals, and median computation time (in seconds) with their respective MAD (median absolute deviation) for 500 simulations using VNPC and SGVB methods across different sample sizes ( n = 256, 512, and 1024) for both VAR(2) and VMA(1) models.  \n4  \nFIG. 4: Relationship between the number of basis functions ( M ) and the normalized log maximum likelihood estimate (MLE) for different scenarios of Einstein Telescope noise analysis. Case A (blue line) represents correlated noise, while Cases B and C (green line) represent uncorrelated noise scenarios.  \n<!-- image -->  \nThe data is divided into 125 equal blocks, each consisting of 16 384 frequency points and a maximum frequency of 1024 Hz for the blocked Whittle likelihood estimation (Equation 5). Each Whittle likelihood block uses identical basis function expressions defined in Section II B for the spectral density matrix.  \nFigure 4 displays the log maximum likelihood estimates (MLE) across a range of M for each case. As the log MLE stabilizes around M ∼ 400, we set M = 450. Using the methodology discussed in Section II, the PSDs for each case are estimated and plotted in Figure 5 and Figure 6. For each case, the diagonal subplots represent the spectral densities for X, Y and Z channels. Since the off-diagonal elements of the PSD matrix at each frequency are complex numbers, with the elements above the diagonal being the conjugates of those below, we plot the real part of the off-diagonal elements in the upper subplot and the imaginary part in the lower subplot. Additionally, the squared coherences are plotted for the cases in Figure 7. The plot shows the squared coherences estimation between any pair of channels for case A and case B.  \nIn all three cases, prominent spectral peaks are observed at the specific frequencies in each channel, consistent with the injected Gaussian noise peaks. This demonstrates the ability of our proposed method to accurately capture and identify these peaks. We note that our PSD estimates exhibit oscillatory features around the Gaussian peaks, this is most likely linked to the relatively loud and sharp Gaussian peak features. This highlights one of the potential weaknesses of the SGVB method: its potential difficulty in fitting very sharp features. For more realistic data sets including correlated noise from magnetic and/or Newtonian noise origin, this will most likely  \nFIG. 5: Case A PSD. The plot shows the periodogram (gray), the PSD estimates (colored, with 90% credible interval shaded) and the true PSD (black curve in the diagonals). The upper triangle subplots show the real parts of the off-diagonal elements of the spectral density matrix, while the lower triangle subplots show the imaginary parts of the off-diagonal elements.  \n<!-- image -->  \nform less of an issue due to the smooth, power-law-like behavior of such noise sources Janssens et al. [9], Janssens et al. [10, 11, 13]. However, this could form an issue with actual detector data which is contaminated with a large amount of instrumental spectral noise lines. Future work will aim to investigate the capabilities of the SGVB method and its efficiency in deal with such complexity us-  \ning real GW data and whether pre-cleaning steps should be undertaken. However, the dramatic speedup of PSD estimation using the SGVB remains an extremely potent benefit for next generation GW detectors. Without faster novel noise estimates techniques, such as the one proposed in this work, the resulting science will very likely be limited by available computation time.  \nFIG. 6: Case B and C PSDs. The plot shows the periodogram (gray), the PSD estimates (colored, with 90% credible interval shaded) and the true PSD (black curve in the diagonals). The upper triangle subplots show the real parts of the off-diagonal elements of the spectral density matrix, while the lower triangle subplots show the imaginary parts of the off-diagonal elements.  \n<!-- image -->  \nIn case A, where correlated noise is intentionally introduced, substantial coherence between channel pairs is observed at their shared peak frequencies (see the top panel of Figure 7). Specifically, high coherence is evident between channels X and Y at 10 Hz, X and Z at 50 Hz, and Y and Z at 90 Hz. The median estimates for the squared coherence aligns with the true coherence, demonstrat-  \ning that our estimate accurately reflects the underlying correlated noise structure. This alignment validates the effectiveness of our spectral analysis in capturing interchannel relationships.  \nCase B, which considers ET noise data with nominally uncorrelated noise, reveals minor fluctuations in the off-diagonal subplots representing the real and imag-  \nC  \nFIG. 7: Squared coherence estimation for any pair of channels in cases A (top) and B (bottom). The solid curves represent the median estimated squared coherences, while the shaded areas indicate the corresponding 90% credible interval. The black solid line is the true squared coherence, and the black dashed line is the expected coherence from uncorrelated white noise data.  \n<!-- image -->  \ninary parts of the estimated spectral densities, despite the absence of deliberately injected correlations. These fluctuations manifest as small coherence values (median C xy < 0 . 008) that fluctuate around the true zero coherence in the corresponding squared coherence plots of Figure 7 (visible as they are plotted on a log-scale). Notably, the median squared coherence fluctuations are lower than those expected from uncorrelated white noise data, which scale as approximately the inverse of the number of data segments used to average over. In this case, i.e. 1 /N b = 0 . 008. The squared coherence estimates based on the SGVB method are below this estimate and hence are consistent with zero coherence.  \nIn case C, where the three channels are assumed independent, the model assumes zero cross spectrum (the XY, YZ, XZ components) in the PSD estimates. Consequently, the squared coherence between any pair of channels is uniformly zero across all frequencies, which is consistent with theoretical expectations for uncorrelated channels, providing a baseline to assess the degree of correlation in the other cases. Comparing case B and C highlights an important feature about the effectiveness of our SGVB model. As one would expect, a model with less complexity (case C) typically yields a more accurate estimate compared to a more complex model (case B). More concretely, the average root mean square deviation (RMSD) 4 over the frequency band 5 to 128 Hz is 9 . 773 × 10 -50 for case B and 9 . 725 × 10 -50 for case C. The  \ndifference between both different scenario's is relatively small, case C's RMSD is about 0.5% lower than that of case B, showcasing the huge potential and flexibility of the SGVB approach in estimating detector noise whether correlated noise is present or not.", 'V. CONCLUSIONS': "Recent work has pointed out the importance of taking into account correlated interferometer network noise when estimating the parameters of gravitational wave signals observed by ground-based detectors such as ET or space-based detectors such as LISA [7-9]. However, so far, no statistical methodology has been suggested to estimate the multivariate noise spectral density of a detector network. In this study, we propose a computationally efficient variational Bayesian approach for estimating the spectral density of multivariate time series. It is important to note that this technique does not assume any restrictive parametric form of the spectral density, allowing it to adapt to any shape of the PSD. To handle the large length of GW network measurements, we use a blocked Whittle likelihood. We provide a hyperparameter optimization method and guidelines to tune the settings of the method. Simulation studies presented for VAR(2) and VMA(1) demonstrate that our approach achieves accuracy comparable to MCMC while significantly reducing computation time, although with potentially narrower credible intervals.  \nThe application of our method to simulated ET noise data revealed its ability to accurately identify and characterize spectral features across different noise scenarios. In the correlated noise case, our approach successfully detected the injected Gaussian peaks and quantified the coherence between channels at specific frequencies.  \nThese findings have significant implications for GW data analysis. The ability to efficiently and accurately estimate spectral densities and inter-channel correlations is crucial for optimizing detection algorithms and understanding noise characteristics in GW detectors. Our method's computational efficiency makes it particularly suitable for analyzing large-scale datasets expected from next-generation detectors (e.g. ET and LISA). The flexibility of our method, combined with its computational efficiency, makes it particularly suitable for analyzing large-scale datasets expected from next-generation GW detectors such as ET and LISA.  \nFuture work could explore the application of this method to real GW detector data and its extension to simultaneously analyzing a stochastic GW background signal. Jointly analyzing the six channels of two detector networks such as TianQin and LISA could provide an more flexible alternative to cross-correlation detection [55]. Further investigations could examine the implications of these spectral characteristics on GW detection sensitivity and methods to mitigate correlated noise effects in ET or LISA. Extending the analysis to  \nlonger time series or different ET configurations could provide insights into the stability and generalizability of these spectral features. Additionally, investigating the method's performance on GPUs could further validate its scalability, robustness and versatility.  \nIn conclusion, the variational Bayes approach offers a promising tool for spectral density estimation in GW data analysis, combining accuracy with computational efficiency. This method has the potential to enhance our ability to characterize detector noise and ultimately improve GW detection capabilities.", 'DATA AND SOFTWARE AVAILABILITY': 'The software developed for this research is available on GitHub under the MIT License [56]. Documentation and examples of the software can be found at https://nz-gravity.github.io/sgvb\\_psd/ . The raw data and data products for this work have been archived on Zenodo [57].  \n- [1] M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, et al. , The einstein telescope: a thirdgeneration gravitational wave observatory, Classical and Quantum Gravity 27 , 194002 (2010).\n- [2] M. Evans, R. X. Adhikari, C. Afle, S. W. Ballmer, S. Biscoveanu, S. Borhanian, et al. , A horizon study for cosmic explorer: Science, observatories, and community (2021), arXiv:2109.09882 [astro-ph.IM].\n- [3] P. Amaro-Seoane, J. Andrews, M. Arca Sedda, A. Askar, Q. Baghi, R. Balasov, et al. , Astrophysics with the laser interferometer space antenna, Living Reviews in Relativity 26 , 10.1007/s41114-022-00041-y (2023).\n- [4] ET steering committee, Einstein telescope: Science case, design study and feasibility report (2020).\n- [5] M. Maggiore, C. V. D. Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, et al. , Science case for the einstein telescope, Journal of Cosmology and Astroparticle Physics 2020 (03), 050.\n- [6] M. Branchesi, M. Maggiore, D. Alonso, C. Badger, B. Banerjee, F. Beirnaert, et al. , Science with the einstein telescope: a comparison of different designs, Journal of Cosmology and Astroparticle Physics 2023 (07), 068.\n- [7] ET steering committee, ET Design Report Update 2020 (2020).\n- [8] M. Colpi, K. Danzmann, M. Hewitson, K. HolleyBockelmann, P. Jetzer, and et al., LISA Definition Study Report, arXiv e-prints , arXiv:2402.07571 (2024), arXiv:2402.07571 [astro-ph.CO].\n- [9] K. Janssens, G. Boileau, M.-A. Bizouard, N. Christensen, T. Regimbau, and N. v. Remortel, Formalism for power spectral density estimation for non-identical and correlated noise using the null channel in Einstein Telescope, European Physical Journal Plus 138 , 352 (2023), arXiv:2205.00416 [gr-qc].\n- [10] K. Janssens, G. Boileau, N. Christensen, F. Badaracco, and N. van Remortel, Impact of correlated seismic and', 'ACKNOWLEDGMENTS': "We thank Zhixiong Hu and Racquel Prado for making their code publicly available and for helpful discussions. AV, JEL, NC, PMR and RM gratefully acknowledge support by the Marsden Fund Council grant MFP-UOA2131 from New Zealand Government funding, managed by the Royal Society Te Ap¯arangi. KJ was supported by FWOVlaanderen via grant number 11C5720N during part of this work. We thank the New Zealand eScience Infrastructure (NeSI) https://www.nesi.org.nz for the use of their high performance computing facilities and the Centre for eResearch at the University of Auckland for their technical support.  \nSoftware used: Python [58], TensorFlow Probability [59, 60] numpy [61], scipy [62], pandas [63], matplotlib [64], Jupuyter Book [65].  \n- correlated newtonian noise on the einstein telescope, Phys. Rev. D 106 , 042008 (2022).\n- [11] K. Janssens, G. Boileau, N. Christensen, N. van Remortel, F. Badaracco, B. Canuel, A. Cardini, A. Contu, M. W. Coughlin, J.-B. Decitre, R. De Rosa, M. Di Giovanni, D. D'Urso, S. Gaffet, C. Giunchi, J. Harms, S. Koley, V. Mangano, L. Naticchioni, M. Olivieri, F. Paoletti, D. Rozza, D. O. Sabulsky, S. Shani-Kadmiel, and L. Trozzo, Correlated 0.01-40 hz seismic and newtonian noise and its impact on future gravitational-wave detectors, Phys. Rev. D 109 , 102002 (2024).\n- [12] K. Janssens, M. Ball, R. M. S. Schofield, N. Christensen, R. Frey, N. van Remortel, et al. , Correlated 1-1000 hz magnetic field fluctuations from lightning over earthscale distances and their impact on gravitational wave searches, Phys. Rev. D 107 , 022004 (2023).\n- [13] K. Janssens, K. Martinovic, N. Christensen, P. M. Meyers, and M. Sakellariadou, Impact of schumann resonances on the einstein telescope and projections for the magnetic coupling function, Phys. Rev. D 104 , 122006 (2021).\n- [14] G. Anderson, J. Anderson, M. Anderson, et al. , Experimental results from the ST7 mission on LISA Pathfinder, Phys. Rev. D 98 , 102005 (2018).\n- [15] M. Armano, H. Audley, J. Baird, et al. (LISA Pathfinder Collaboration), LISA Pathfinder micronewton cold gas thrusters: In-flight characterization, Phys. Rev. D 99 , 122003 (2019).\n- [16] E. Thrane, N. Christensen, and R. M. S. Schofield, Correlated magnetic noise in global networks of gravitationalwave detectors: Observations and implications, Phys. Rev. D 87 , 123009 (2013).\n- [17] N. Christensen, Stochastic gravitational wave backgrounds, Reports on Progress in Physics 82 , 016903 (2018).\n- [18] G. Boileau, N. Christensen, and R. Meyer, Figures of\n- merit for a stochastic gravitational-wave background measurement by lisa: Implications of lisa pathfinder noise correlations, Physical Review D 106 , 063025 (2022).\n- [19] F. Cireddu, M. Wils, I. C. F. Wong, P. T. H. Pang, T. G. F. Li, and W. Del Pozzo, Likelihood for a Network of Gravitational-Wave Detectors with Correlated Noise, arXiv e-prints , arXiv:2312.14614 (2023), arXiv:2312.14614 [gr-qc].\n- [20] A. Meier, C. Kirch, and R. Meyer, Bayesian nonparametric analysis of multivariate time series: A matrix gamma process approach, Journal of multivariate analysis 175 , 104560 (2020).\n- [21] Y. Liu, C. Kirch, J. E. Lee, and R. Meyer, A nonparametrically corrected likelihood for bayesian spectral analysis of multivariate time series (2024).\n- [22] Z. Hu and R. Prado, Fast bayesian inference on spectral analysis of multivariate stationary time series, Computational Statistics & Data Analysis 178 , 107596 (2023).\n- [23] M. I. Jordan, Z. Ghahramani, T. Jaakkola, and L. Saul, Introduction to variational methods for graphical models, Machine Learning , 183-233.\n- [24] M. J. Wainwright and M. I. Jordan, Graphical models, exponential families, and variational inference, Foundations and Trends in Machine Learning , 1.\n- [25] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, Variational inference: A review for statisticians, Journal of the American Statistical Association 112 , 859-877 (2017).\n- [26] D. M. Blei and M. I. Jordan, Variational inference for dirichlet process mixtures, Bayesian Analysis , 121-143 (2006).\n- [27] D. P. Kingma and M. Welling, Auto-encoding variational bayes (2022), arXiv:1312.6114 [stat.ML].\n- [28] M. Vallisneri, M. Crisostomi, A. D. Johnson, and P. M. Meyers, Rapid parameter estimation for pulsar-timingarray datasets with variational inference and normalizing flows (2024), arXiv:2405.08857 [gr-qc].\n- [29] E. Thrane and C. Talbot, An introduction to Bayesian inference in gravitational-wave astronomy: Parameter estimation, model selection, and hierarchical models, Publications of the Astronomical Society of Australia 36 , e010 (2019), arXiv:1809.02293 [astro-ph.IM].\n- [30] N. Christensen and R. Meyer, Parameter estimation with gravitational waves, Rev. Mod. Phys. 94 , 025001 (2022).\n- [31] B. A. Vu, D. Gunawan, and A. Zammit-Mangion, Recursive variational Gaussian approximation with the Whittle likelihood for linear non-Gaussian state space models, arXiv e-prints , arXiv:2406.15998 (2024), arXiv:2406.15998 [stat.CO].\n- [32] O. Rosen and D. S. Stoffer, Automatic estimation of multivariate spectra via smoothing splines, Biometrika 94 , 335 (2007).\n- [33] J. Piironen and A. Vehtari, Sparsity information and regularization in the horseshoe and other shrinkage priors, Electronic journal of statistics 11 , 5018 (2017).\n- [34] M. Xu, M. Quiroz, R. Kohn, and S. A. Sisson, Variance reduction properties of the reparameterization trick, in Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics , Proceedings of Machine Learning Research, Vol. 89, edited by K. Chaudhuri and M. Sugiyama (PMLR, 2019) pp. 27112720.\n- [35] J. Domke, Provable gradient variance guarantees for black-box variational inference, in Advances in Neural Information Processing Systems , Vol. 32, edited by H. Wal-  \n- lach, H. Larochelle, A. Beygelzimer, F. d ' Alch'e-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).\n- [36] J. Piironen and A. Vehtari, Sparsity information and regularization in the horseshoe and other shrinkage priors, Electronic Journal of Statistics 11 , 5018 (2017).\n- [37] M. Welandawe, M. R. Andersen, A. Vehtari, and J. H. Huggins, A framework for improving the reliability of Black-box variational inference, Journal of Machine Learning Research 25 , 1 (2024).\n- [38] A. Agrawal, D. Sheldon, and J. Domke, Advances in black-box vi: Normalizing flows, importance weighting, and optimization, in 34th Conference on Neural Information Processing Systems (NeurIPS 2020) (2020).\n- [39] J. Bergstra, D. Yamins, and D. Cox, Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures, in Proceedings of the 30th International Conference on Machine Learning , Proceedings of Machine Learning Research, Vol. 28, edited by S. Dasgupta and D. McAllester (PMLR, Atlanta, Georgia, USA, 2013) pp. 115-123.\n- [40] J. Bergstra, R. Bardenet, Y. Bengio, and B. K'egl, Algorithms for hyper-parameter optimization, in Advances in Neural Information Processing Systems , Vol. 24, edited by J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Weinberger (Curran Associates, Inc., 2011).\n- [41] V. Sakkalis, Review of advanced techniques for the estimation of brain connectivity measured with eeg/meg, Computers in Biology and Medicine 41 , 1110 (2011), special Issue on Techniques for Measuring Brain Connectivity.\n- [42] C. W. J. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica 37 , 424 (1969).\n- [43] B. Wang and D. M. Titterington, Inadequacy of interval estimates corresponding to variational bayesian approximations, in Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics , edited by R. G. Cowell and Z. Ghahramani (Society for Artificial Intelligence and Statistics, 2005) pp. 373-380.\n- [44] D. P. Kingma and M. Welling, Auto-Encoding Variational Bayes, arXiv e-prints , arXiv:1312.6114 (2013), arXiv:1312.6114 [stat.ML].\n- [45] R. Danilo and M. Shakir, Variational inference with normalizing flows, in Proceedings of the 32nd International Conference on Machine Learning, PMLR , Vol. 37 (2015) pp. 1530-1538.\n- [46] Y. Yang, D. Pati, and A. Bhattacharya, α -variational inference with statistical guarantees, The Annals of Statistics 48 , 886 (2020).\n- [47] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, beta-VAE: Learning basic visual concepts with a constrained variational framework, in International Conference on Learning Representations (2017).\n- [48] Y. Wang and D. Blei, Variational bayes under model misspecification, in Advances in Neural Information Processing Systems , Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d ' Alch'e-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).\n- [49] Y. Yao, A. Vehtari, D. Simpson, and A. Gelman, Yes, but did it work?: Evaluating variational inference, in Proceedings of the 35th International Conference on Machine Learning , Proceedings of Machine Learning Research, Vol. 80, edited by J. Dy and A. Krause (PMLR,  \n2018) pp. 5581-5590.  \n- [50] J. Lee, G. K. Nicholls, and R. J. Ryder ‡ , Calibration procedures for approximate Bayesian credible sets, Bayesian Analysis 14 , 1245 (2019).\n- [51] X. Yu, D. Nott, M.-N. Tranc, and N. Kleind, Assessment and adjustment of approximate inference algorithms using the law of total variance, JOURNAL OF COMPUTATIONAL AND GRAPHICAL STATISTICS 30 , 977 (2021).\n- [52] S. Hild, S. Chelkowski, A. Freise, J. Franc, N. Morgado, R. Flaminio, and R. DeSalvo, A xylophone configuration for a third-generation gravitational wave detector, Classical and Quantum Gravity 27 , 015003 (2009).\n- [53] S. Hild et al. , Sensitivity Studies for Third-Generation Gravitational Wave Observatories, Class. Quant. Grav. 28 , 094013 (2011), arXiv:1012.0908 [gr-qc].\n- [54] M. Branchesi, M. Maggiore, D. Alonso, C. Badger, B. Banerjee, F. Beirnaert, E. Belgacem, S. Bhagwat, G. Boileau, S. Borhanian, et al. , Science with the Einstein Telescope: a comparison of different designs, J. Cosmol. Astropart. Phys. 2023 (7), 068, arXiv:2303.15923 [gr-qc].\n- [55] Z.-C. Liang, Z.-Y. Li, E.-K. Li, J. dong Zhang, and Y.-M. Hu, Unveiling a multi-component stochastic gravitational-wave background with the tianqin + lisa network (2024), arXiv:2409.00778 [gr-qc].\n- [56] J. Liu, A. Vajpeyi, R. Meyer, K. Janssens, J. E. Lee, P. Maturana-Russel, N. Christensen, and Y. Liu, sgvb psd (2024).\n- [57] J. Liu, A. Vajpeyi, R. Meyer, K. Janssens, J. E. Lee, P. Maturana-Russel, N. Christensen, and Y. Liu, Dataset accompaning 'Variational inference for correlated gravitational wave detector network noise', 10.5281/zenodo.13789536 (2024).\n- [58] G. van Rossum and F. L. Drake Jr, Python programming  \nlanguage, Python Software Foundation (2020).  \n- [59] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Man'e, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Vi'egas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Largescale machine learning on heterogeneous systems (2015), software available from tensorflow.org.\n- [60] J. V. Dillon, I. Langmore, D. Tran, E. Brevdo, S. Vasudevan, D. Moore, B. Patton, A. Alemi, M. Hoffman, and R. A. Saurous, Tensorflow distributions, arXiv preprint arXiv:1711.10604 (2017).\n- [61] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R'ıo, M. Wiebe, P. Peterson, P. G'erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585 , 357 (2020).\n- [62] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al. , Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17 , 261 (2020).\n- [63] T. pandas development team, pandas-dev/pandas: Pandas (2020).\n- [64] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in science & engineering 9 , 90 (2007).\n- [65] The Executable Book Project, Jupyter book (2020)."}