{'fraction_non_alphanumeric': 0.034543715767815224, 'fraction_numerical': 0.017764511548381486, 'mean_word_length': 4.634879163945134, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical techniques for the computation of sample spectral distributions of population mixtures

This note describes some techniques developed for the computation of the sample eigenvalue distribution of random matrices generated by mixtures of populations. Within this model the mapping between the population distributions and the asymptotic sample distribution can be obtained by solving a set of systems of non-linear equations, for which we provide an efficient implementation. This work contributes by describing a method for accelerated fixed point convergence, a homotopy continuation strategy to prevent convergence to non-admissible solutions, a blind non-uniform grid construction for effective distribution support detection and improved approximation at its edges, and a parallel computing architecture. Comparisons are performed with available packages for the single population case and with results obtained by simulation for the more general model implemented here. Results show competitive performance and improved flexibility.


Introduction
Random matrix theory is at the core of modern high dimensional statistical inference with applications in physics, biology, economics, communications or imaging (Couillet and Debbah, 2013;Paul and Aue, 2014;Bun et al., 2017). In the large dimensional scenario, classical asymptotics where the available number of samples of a given population N is much larger than the population size M are no longer valid. A key contribution in this setting is the Marčenko-Pastur theorem (Marčenko and Pastur, 1967), which provides an explicit characterization of the limiting behaviour of the sample or empirical eigenvalue distribution for matrices with independent and identically distributed (IID) entries in terms of the matrix aspect ratio γ = M/N . Using this characterization, subsequent asymptotics-based inference can be performed, for instance, for dimensionality reduction, hypothesis testing, signal retrieval, or covariance estimation. A different version of the Marčenko-Pastur theorem was proved in Silverstein and Bai (1995) for matrices with rows drawn from a distribution that converges asymptotically to a given population distribution, not necessarily IID. In this case, they showed that there exist a fixed point equation relating the eigenvalues of the sample and population distributions, which can be used for inference under the mediation of appropriate numerical techniques. In addition, in Silverstein and Choi (1995), they provided a simple characterization of the support of the sample distribution based on the monotonicity of the inverse of the function involved in the fixed point equation. Emerging from statistical models in array processing, the extension in Wagner et al. (2012) focuses on generalizing previous results to the case where the matrix rows are independent but drawn from a collection of population distributions. Here, the relation between the population covariances and the sample eigenvalues is governed by a system of non-linear equations and, although sample eigenvalue confinement has also been proved (Kammoun and Alouini, 2016), a simple characterization of the support is no longer at hand.
Despite the number of potential applications of this theory, not many works have practically addressed the numerical issues involved in computing the sample distribution from a given population distribution. The most flexible package that we have identified has been described in Dobriban (2015) (SPECTRODE), where the fixed point equation in Silverstein and Bai (1995) is transformed into an ordinary differential equation (ODE) with starting point obtained by the solution of the fixed point equation on a single point within the support. In this work, the authors show that their method compares favourably with straightforward fixed point solvers, both in terms of accuracy and computational efficiency. In addition, they showcase the accuracy limitations of Monte Carlo simulations (Jing et al., 2010). Finally, they provide arguments about the limited applicability or lack of efficiency of other approaches. In our experiments, this package has revealed an exquisite accuracy and efficiency in computing the sample distribution and estimating its support. However, the applicability of an ODE approach may be limited for more general models such as Wagner et al. (2012). This is mainly due to the increased computational complexity of evaluating the Jacobian of the system of fixed point equations, which no longer depends on the eigenvalues of the populations but on a combination of their covariances. Another interesting tool, conceived to solve the more general problem of recovering the population distribution from the sample distribution by means of the so called QuEST function, has been described in Ledoit and Wolf (2017). Solving this problem is required, for instance, for covariance estimation. However, for other applications, such as denoising, the interest may lie in accessing the sample distribution of noise from the knowledge of its statistical properties. The literature review in this work points to a systematic limitation of most precedent methods, as they are only capable to obtain estimates for very particular forms of the population distribution. In addition, it introduces an interesting feature not present in Dobriban (2015), the use of a non-uniform grid with increased resolution near the support edges, which allows more efficient approximations. However, once again, this technique is designed for cases where a single nonlinear equation is to be solved while generalizations to systems of equations do not seem straightforward.
In this technical note we describe a series of tools that we have developed to obtain sample distribution estimates for the mixture of populations case in (Wagner et al., 2012). This model or certain analogous forms, has attracted interest in works such as Moustakas and Simon (2007) (2019). Our method is based on directly solving the system of nonlinear equations. However, the results in Dobriban (2015) and our own analysis, have identified certain limitations in commonly reported algorithms based on fixed point iteration solvers, so a set of technical refinements are proposed in this note. These include the use of Anderson mixing to accelerate the fixed point iterations, an homotopy continuation method to prevent non-admissible solutions, a set of heuristics to detect the support of the distribution and to adapt the approximation grid to the sample distribution shape, and a formulation that allows for efficient computations in graphical processing units (GPU). We validate our method by comparison with Dobriban (2015) and Ledoit and Wolf (2017), both in terms of computational efficiency and accuracy. As our methods are envisaged to operate on more general models than those contemplated in Dobriban (2015); Ledoit and Wolf (2017), they do not make use of any precomputed information about the distribution support. Nevertheless, we show that our approach is reliable enough and its efficiency and accuracy is comparable or superior to that in Dobriban (2015); Ledoit and Wolf (2017). In addition, comparisons with Monte Carlo simulations show that our technique is also capable of providing accurate estimates for more general models. A MATLAB implementation of our approach, that we will refer to as the MIXANDMIX (Mixtures by Anderson Mixing) method, including the scripts required to replicate the experiments in this note, has been made publicly available at https: //github.com/mriphysics/MixAndMix/releases/tag/1.0.0. This note is organized as follows: in § 2 we review different random matrix models, in § 3 we describe the main functionalities of our method, in § 4 we validate our technique, in § 5 we discuss the implications of the obtained results and in § 6 we end up with some conclusions.

Theory
Consider a complex random matrix X of size N × M . We are interested in the eigenvalue distribution of the sample covariance Y = 1 N X H X in the asymptotic regime where both M → ∞ and N → ∞ but they maintain a fixed ratio γ = M/N with 0 < γ < 1. For simplicity we assume the entries of the matrix are zero mean Gaussian distributed, but keeping in mind that the literature contemplates different generalizations. Three main scenarios are considered:

IID standard entries
Within this model, that we call standard model, we can write X H ∼ CN (0 M N , I M N ), with CN denoting a circularly symmetric complex Gaussian distribution, 0 M N the M × N matrix with zero entries, and I M N the M N × M N identity matrix. Marčenko and Pastur (1967) showed that in this case the eigenvalues of the sample covariance follow with γ − = (1 − √ γ) 2 and γ + = (1 + √ γ) 2 defining the support of the distribution. Thus, in this case, (1) gives us an explicit characterization of the sample spectrum.

IID rows
In this scenario, we can write X H ∼ CN (0 M N , Λ M ⊗ I N ), with Λ M a given population covariance matrix. Silverstein and Bai (1995) characterized the sample distribution using its Stieltjes transform for which the inversion formula would give back the sample distribution by If, in a discrete setting, we denote the increasingly sorted eigenvalues of Λ M by {λ 1 , . . . , λ M }, Silverstein and Bai (1995) showed that the Stieltjes transform of the sample distribution of the matrixỸ = 1 N XX H , m(z), is related to the population distribution by the fixed point equationm(z) = ) −1 and to the Stieltjes transform of the spectral distribution of the sample covariance matrix Y, m(z), bỹ m(z) = γm(z) + (γ − 1)/z. In addition, Silverstein and Choi (1995) showed that the support of the distribution can be characterized by studying the zeros of the derivative of the inverse map, z ′ (m), in appropriate domains. This property, together with the analyticity of the sample density within its support are the keys to the SPECTRODE approach in Dobriban (2015).
To clarify the connections between the different models, it is more convenient to express the previous relations in terms of auxiliary functions e(z) = − 1 γzm(z) − 1 γ for an equivalent fixed point equation, and an expression for the Stieltjes transform of the sample distribution, We refer to this model as the single population model.

Independent rows
In this mixture of populations model, the matrix is drawn from , m(z), and the auxiliary functions e j (z), 1 ≤ j ≤ K. These functions are related to the population covariances by a system of nonlinear equations for which there is a unique solution in C \ R + . These functions are used to obtain an expression for m(z), with α k = tr(D k )/N . Note that (6) and (7) reduce to (4) and (5) when K = 1. There are two main limitations to extend the SPECTRODE method to this setting. First, we are unaware of studies characterizing the support of the measures inducing e j (z) by means of some analogue to Silverstein and Choi (1995) results. Second, both potential extensions of support characterizations or usage of ODE solvers would require the Jacobian of (6), which involves additional matrix multiplications, with a negative impact in computational efficiency. Thus, we have focused on developing a reliable method to solve the system (6) not requiring the Jacobian.

Methods
In this Section we describe our numerical solver for the system of equations in (6).

Support detection
When γ → 0 the sample distribution tends to the population distribution for the standard and single population models, while for the mixture of populations, the spectral distribution will be governed by an Thus, in this limiting case the sample distribution support only includes the eigenvalues of an equivalent single population distribution. On the other side, , with λ k m denoting the m-th eigenvalue of Λ k M and t > 1, provide lower and upper bounds on the sample spectral distribution support. To consider these two extreme cases, we group and sort the set of P = M (K + 1) eigenvalues λ p = {λ k m , λ m }, 1 < k < K, 1 < m < M . Our method attends to the overlap of the set of intervals i.e., the intervals generated by the maximum expected spectral dispersion of each grouped eigenvalue, as given by the worst case scenario in (1), including a numerical safety margin t, which we have set to t = 1.001. Lack of overlap among adjacent eigenvalues is detected by checking the condition b p < a p+1 . Assuming this is observed I ≤ P − 1 times, we can get a partition into I + 1 segments, each one induced by P i eigenvalues, 1 ≤ i ≤ I +1. Each of these support segments is gridded using a total of max(M (o) P i , M (i) ) points, with M (i) the minimum number of points per segment and M (o) ≤ M (i) the minimum ratio of points per number of grouped eigenvalues. Importantly, to improve detectability, by noting the multiplicative dependence of the spectral dispersion width with the spectral location in (8), gridding is performed uniformly in logarithmic units. Then, we can call the solver of (6) (to be described in § 3.3), and compute (7) and (3) over these grid locations. The output of this first step is a set of estimates for the distribution in a non-uniform grid The lower bound on the number of grid points M (o) P will be attained for λ → 1 and clustered eigenvalues, and the upper bound M (i) P when λ → 0. Fixing M (o) = 3 and M (i) = 15 has detected at least a single point within all the segments conforming the support of the distribution for the range of problems studied in § 4. For problems where the total geomeric multiplicity of the discretized population covariance is given by M (u) , equivalent expressions of the problem can be obtained for any M = SM (u) with S ∈ N >0 . In our experiments S has been selected by defining a minimum number of grid points to approximate the densities, M (m) , and making M = max(M (u) , M (m) ) with M (m) = 100.

Adaptive regridding
Additional points are added to the grid by pursuing g ′ (x) √ xf ′′ (x) = c, with g(x) the grid mapping function and c a constant. This non-uniform grid construction criterion is based on both the second order derivative of the density f ′′ (x) and the grid value x. The first feature, f ′′ (x), favours the allocation of grid points near the support edges, in accordance to the √ x − x 0 behaviour of the distribution at the boundaries (Silverstein and Choi, 1995), as well as in those areas where linear interpolation results in larger approximation errors. This is similar in spirit to the arcsine criterion in Ledoit and Wolf (2017) but does not use any prior information about the support edges as it is not available in the mixture of populations model. The second feature, x, favours the allocation of grid points near the upper edge of the support, which could be important for applications related to signal detection (Nadakuditi, 2014;Dobriban and Owen, 2019). After P l = R l P 0 points are added to the grid, the solver for f (x) is called on the new set of points to allow for an update of the f ′′ (x) values to be used at the next iterative regridding step. This whole process is repeated L times, so R l , 1 ≤ l ≤ L control the final resolution of the distribution computations. The parameters by default in our implementation are R l = 1 ∀l and L = 1.

Homotopy continuation
The calculation of the sample distribution involves a pass to the limit in (3) as the Stieltjes transform does not converge in the real line. In addition, the solution of (6) is not unique on the real line. Numerically, this may provoke spurious fixed point convergence when the current solution is far away from the optimum and the computations are being performed in locations that are close to the real line. To prevent these situations, we have emulated (3) by homotopy continuation. We start by obtaining an approximate solution for (6) in a grid given by z = x + ξi, with ξ = ξ 0 1 P l for big enough ξ 0 common for all 1 ≤ p ≤ P l . At each iteration i we compute ε i p = max where the updates on e are to be described in § 3.4. Considering that, as discussed in Dobriban (2015), to obtain an accuracy of at least ϵ for the distribution f (x), we need to solve the system of equations in a complex grid given by x + iϵ 2 , we can perform the update ξ j+1 p = max(ξ j p /β 2 , ϵ 2 ) whenever ε i p ≤ ε i−1 p . The iterations at the grid location indexed by p are terminated when the prescribed accuracy is reached, namely when ξ j p = ϵ 2 and ε i p < ϵ. In our experiments we have used ξ 0 = 1 and β = 10, for which we have observed robust and efficient performance.

Anderson acceleration
The experiments in Dobriban (2015) showed that when solving the IID rows problem in § 2.2 by a straightforward fixed point algorithm, in our context when performing the updates on e(z) directly using (4), convergence is often very slow, so they reported situations where their SPECTRODE code could be 1000× quicker while simultaneously obtaining 1000× higher accuracy. In this note we show that this apparent limitation of the fixed point iterates can be overcome by using techniques to accelerate their convergence. As discussed in § 2.3, the accelerated convergence achievable by methods requiring the Jacobian of the fixed point identities may not compensate for the increased cost per iteration involved in computing the Jacobian. Thus, we have resorted to Anderson mixing (Anderson, 1965), a technique not requiring explicit Jacobian calculations that has demonstrated good practical performance, in occasions providing competitive results when compared to gradient-based approaches (Ramière and Helfer, 2015).
Considering a given multidimensional fixed point mapping g(e) such as (6), Anderson iterations are computed as with Q i denoting the number of iterations whose history is used to compute the update at iteration i and γ i = (γ i 1 , . . . , γ i Q i ) T obtained by solving a linear least squares problem involving the fixed point update h i (e i ) = g(e i )−e i and its differences ∆h Due to the potential ill-posedness of this system, we have actually solved the damped version (Henderson and Varadhan, 2018) γ with damping parameter given by λ i = 0.1 max k,q |∆H i k,q |. We have set Q i = min(2, i − 1) on the basis of our empirical testing and in agreement with the experimental results in Ramière and Helfer (2015).

GPU acceleration
GPU-based implementations of the SPECTRODE method (Dobriban, 2015) appear involved due to the sequential nature of ODE solvers. In contrast, acceleration of the sample spectral distribution computation in the QuEST method (Ledoit and Wolf, 2017) seems more plausible as the method solves for a zero of a function independently for the different grid locations, but the authors have not discussed this aspect. Our code has been architectured so that most demanding routines support both CPU and GPU based parallel computations. This includes the parallel computation of the solutions of the system of equations in (6) for the different grid locations but also the parallel computation of different sample spectra, required, for instance, in patch-based image denoising applications (Cordero-Grande et al., 2018).

Results
In this Section we first justify the beneficial effects of Anderson acceleration and homotopy continuation in sample spectral distribution estimation ( § 4.1), then compare our method to the SPECTRODE and QUEST proposals in those regimes in which they can operate ( § 4.2) and finally provide some results on the application of our technique to the mixture of populations model ( § 4.3). Unless otherwise stated experiments are performed using CPU computations.

Validation of introduced refinements
In Dobriban (2015) an experiment was performed illustrating the limitations of fixed point iterations to obtain accurate estimates of the sample spectral distribution. Their method is compared with a fixed point iteration solver for the standard model in § 2.1, using the closed form density in (1) with γ = 0.5 to assess the accuracy. Here we replicate that experiment adding our Anderson acceleration technique to the fixed point solver. The results are presented in Fig. 1. First, we have been able to replicate the results in Dobriban (2015); the fixed point algorithm, grossly equivalent to the MIXANDMIX implementation with Q = 0, i.e., without Anderson mixing, is only able to provide very moderate accuracies, despite being run for 1/ϵ iterations. However, when introducing the Anderson acceleration scheme, the results are dramatically better, with MIXANDMIX and SPECTRODE demonstrating comparable performance. In this experiment the curves show a slightly better accuracy (Fig. 1a) and worse computational efficiency ( Fig. 1b) for MIXANDMIX, but this should be taken with caution as these tests have been conducted without considering the influence of grid sizes on the approximation, which will be taken into account in the experiments in § 4.2. In addition, we show (Fig. 1c) that despite the accuracy obtained by a straightforward fixed point algorithm is poor everywhere within the support, the accuracy curve after Anderson mixing remains below the SPECTRODE curve almost everywhere.
In Fig 2a we compare the results of our method without and with homotopy continuation to those of SPECTRODE. The SPECTRODE method and ours with homotopy continuation are observed to overlap at the scale of the plot. However, when running MIXANDMIX without homotopy continuation, we can see there exist some grid points for which the computations spuriously converge to the zero solution. We know this solution is infeasible because it provokes discontinuities in the distribution, which contradicts its expected analiticity properties (Silverstein and Choi, 1995). To illustrate the reasons for these numerical issues, we have taken a grid point corresponding to one of these infeasible estimates, x = 2.2. For this point, Figs. 2b-e show the squared magnitude of the residuals of the fixed point maps, |h(e)| 2 , at different complex plane values of the auxiliary function e(z) = e(x + δi) as we are approaching the real line with δ = {1, 0.1, 0.01, 0.001}. First, for z = x + 1i there is a unique minimum in C + whose basin of attraction covers the whole of C + . As we decrease δ (see for instance z = x + 0.1i) we can track this minimum in a neighborhood of its previous location and check that it is still the only one in the upper half of the complex plane. However, a new local minimum has emerged in the lower half of the plane but so close to the real line that its basin of attraction extends to the upper half. As we keep decreasing δ, the basin of attraction of this minimum in C + gets bigger; however, by analyticity there has to be an area around the global optimum for which the method should still converge to the global optimum, which can be ensured by homotopy continuation. This explains the problem we are observing in the left hand side: the fixed point algorithm has entered the basin of attraction of the minimum in the lower half and it has not been able to escape from this area. In addition, the location of the attractor explains why the spurious distribution value obtained by the fixed point algorithm is generally pushed to 0 when failing to converge to the right optimum.

Comparison with the literature
In Fig. 3 we compare the accuracy and computational efficiency of MIXANDMIX with the SPECTRODE and QuEST approaches for two population distributions that admit a closed form expression for the sample Figure 2: a) Sample spectral distribution estimates for the SPECTRODE method, MIXANDMIX including homotopy continuation and MIXANDMIX without homotopy continuation. b-e) Base 10 logarithm of the squared magnitude of the fixed point update, i.e., log 10 (|h(e)| 2 ), together with corresponding isolines at b) z = 2.2 + 1i, c) z = 2.2 + 0.1i, d) z = 2.2 + 0.01i and e) z = 2.2 + 0.001i. spectral distribution. In Figs. 3a,b, we show respectively the averaged accuracy and computation times for a set of aspect ratios ranging from γ = 0.025 to 0.975 in steps of 0.05 for the standard distribution (MP) in (1). Corresponding accuracies throughout the support together with the gold standard density (in a logarithmic scale) are shown in Fig. 3c for the λ = 0.5 case. Analogous plots are provided in Figs. 3d-f for a two-delta (δδ) distribution with equiprobable eigenvalues at λ 1 = 1 and λ 2 = 8, for which the sample spectral distribution can be obtained by solving a third order polynomial equation (Dobriban, 2015). The SPECTRODE and MIXANDMIX approaches have been run with an accuracy parameter providing similar computation times than those of the QuEST method using 100 grid points, which corresponds to ϵ = 10 −6 (MP) and ϵ = 10 −5 (δδ) for SPECTRODE and ϵ = 10 −4 , L = 3 (both) for MIXANDMIX. To account for the relative grid complexities of different methods, linearly interpolated densities are compared with close-form solutions in a uniform grid comprised of 10000 evenly distributed points along the support. MIXANDMIX is roughly 2 and 1 orders of magnitude more accurate than QuEST and SPECTRODE respectively. We observe that the computation times of SPECTRODE largely depend on the aspect ratio, with increments of several orders of magnitude as γ → 1, while they are much more uniform for MIXANDMIX and QuEST. As for the accuracy distributions, they are generally satisfactory for all methods, but MIXANDMIX seems to provide improved estimates near the lower edge for the MP case and throughout the support for the δδ case. The grid sizes used by each method for γ = 0.5 have been 102/104 by QuEST (100 plus some additional points to localize the support limits), 2930/5734 by SPECTRODE and 1200/1200 by MIXANDMIX, for the MP / δδ problems. MIXANDMIX is potentially limited by risks of failures at detecting all the segments comprising the distribution support. Thereby, we have conducted some experiments to test the support detection reliability in challenging scenarios. In Fig. 4 we cover analogous experiments to those in Fig. 3 for skewed δδ . 4c,f show that MIXANDMIX has been capable to approximate the support for both problems and to provide more accurate results than the other two methods almost everywhere within the support. In this case the grid sizes have been 104/104 for QuEST, 6963/9130 for SPECTRODE and 1248/1248 for MIXANDMIX for the LS / RS problems.
In Fig. 5 we apply the three methods to a challenging comb-like (Dobriban, 2015) problem. In this case the population eigenvalues come from 100 equiprobable point masses uniformly distributed in the interval [0.1, 10], thus with a large number of components that differ by as much as two orders of magnitude. We show the results at selected areas of the support for γ = 0.025, γ = 0.5 and γ = 0.975 respectively in Figs. 5a-c, again using logarithmic scaling for improved visualization. When tuning the SPECTRODE parameter for similar computation times to those of QuEST and MIXANDMIX, which turned out to happen at ϵ = 10 −3 , the visual impression is of very limited performance. In contrast, MIXANDMIX appears to be very powerful, showing strengths not only in detecting all the support intervals, but also in picking-up the density oscillations observed for γ = 0.025 and the steeped lower edge for γ = 0.975. We have quantitatively assessed this perception by running the final experiment including also the SPECTRODE results at the accuracy level where we could not visually detect any differences with MIXANDMIX results, ϵ = 10 −6 for both γ = 0.025 and γ = 0.5, or the maximum code accuracy, ϵ = 10 −8 , for γ = 0.975. Grid sizes of the different methods were 1950/1113/1216 for SPECTRODE with matched computation times, 140/106/102 for QuEST, 1248/1200/1200 for MIXANDMIX, and 11678/19156/244296 for SPECTRODE with improved accuracy. Results show no visual differences between MIXANDMIX and SPECTRODE with improved accuracy for γ = 0.025 and γ = 0.5, and more plausible functional shape of the former around the left edge for γ = 0.975, even though approximately 200× less computational resources were used. Finally, QuEST results show a remarkable ability to correctly estimate the support intervals, but strong limitations to accurately approximate the spiked density areas or capture the density oscillations.

Sample spectral distributions of mixtures of populations
In this Section we provide an illustration of the MIXANDMIX results for the mixture of populations model and the benefits of our GPU architecture. Fig. 6a shows the sample spectral distribution estimations corresponding to a mixture of populations with K = 6 equiprobable populations, so α k = 1/K, drawn from population covariances Λ k M (DIAG) = Λ k mn(DIAG) = ((m + k) mod K) + 1)δ[m, n], i.e., a diagonal matrix whose elements in the diagonal grow from 1 to K with period K with this pattern being shifted across the populations. Fig. 6b extends this DIAG problem to non-diagonal population covariances using Λ k M (CORR) = Λ k mn(CORR) = ρ |m−n| l √ Λ k mm(DIAG) Λ k nn(DIAG) , ρ < 1, l > 0, so this CORR problem reduces to the DIAG problem when l → ∞. Namely, Fig. 6b shows the results for ρ = 0.2, l = 0.25 and γ = 0.5, here with M = 120 for convenience. We can see that estimates are in agreement with simulations (with these being obtained using the biggest matrix sizes that we would fit in our GPU memory) for both the DIAG and CORR problems, with the CORR results showing a larger spectral dispersion due to the larger condition number of the non-diagonal matrix. Finally, Fig. 6c shows the GPU computation times for the CORR problem for a number of populations ranging from K = 1 to K = 6. First, we can appreciate a significant penalty when moving from K = 1 to K = 2, as that switches the problem from solving a single equation based only on the eigenvalues to a system of equations based on the whole structure of the covariance matrices. Second, we observe that for K ≥ 2, despite the number of equations to be solved and thus the dimensionality of the parameter space of our problem grows linearly with K, the computation times grow in a sublinear manner. This is to be attributed to an increased degree of parallelization of the GPU implementation for bigger problems and a stable fixed point Anderson acceleration in the multidimensional case.

Discussion
We have presented a set of numerical techniques to aid in the computation of the sample spectral distribution in a mixture of populations model. These include the use of Anderson mixing to accelerate the fixed point iterations, homotopy continuation for robust convergence to the right optimum, adaptive grid construction to efficiently detect the support and approximate the distribution, and a parallel architecture to tackle the increased computational demands in this setting. Results have shown that our method offers favorable practical efficiency-accuracy tradeoffs when compared with related approaches while being able to address more general models than available tools.
Our tool has focused on the model described in Wagner et al. (2012) but it could straightforwardly be adapted and optimized to address several analogous models in the literature, such as those mentioned in § 1. Nevertheless, we have not considered even more general models such as the Kronecker model (Zhang et al., 2013) or more general couplings between the matrix elements (Wen et al., 2011;Lu et al., 2016). However, these models generally involve the solution of more intricate systems of nonlinear equations with more auxiliary functions, but with strong functional resemblances to the system we have studied here, so there is also potential to reuse or adapt our tools to tackle them.
While our technique does not fully exploit the existing descriptions used for support detection in the single population model because we are not aware of any such descriptions for mixtures of populations, the gridding procedure has shown a robust behaviour in challenging practical scenarios, even when compared with methods that exploit those descriptions. The main reason is the introduction of a logarithmic grid, that enables efficient support searches at different spectral scales. This has been synergistically combined with an adaptive grid refinement making use of the second order derivative of the density (with a bias towards the upper edge information) generally with more efficient approximations than provided by previous methods, particularly for the challenging γ → 1 case. Although this grid refinement criterion has shown to be effective for the all the cases we have tested, other criteria may be more appropriate for different applications. In this regard, we should mention that our software provides with a generic grid refinement method that allows to test other possibilities by simply defining different criteria for grid cell subdivision, with some alternatives already included in the code.
Another difference with previous approaches is the dependence of the method on more parameters. Although this may add an extra degree of complexity for users, we have observed good behaviour for all test cases without resorting to parameter tuning, so for many applications it should be ready to use with the default parameters. From a different perspective, the combined inspection of simulations and manipulation of these parameters may allow to fine tune the method in most challenging scenarios, some of them, as shown in some of the experiments in § 4.2, not being adequately covered by the reduced flexibility of related approaches. In summary, grid detection robustness can be improved by increasing M (i) , robustness of approximation by increasing ξ 0 and/or decreasing β and accuracy by decreasing ϵ and/or increasing L.

Conclusions
This work has introduced a set of techniques to estimate sample spectral distributions in a mixture of populations model. A generic procedure using only the functional form of the fixed point equations relating the population and spectral distributions has been proposed. Efficient convergence is achieved by Anderson acceleration and homotopy continuation and novel strategies for grid construction have been provided. This method has compared well with related proposals in the literature which, to our knowledge, are only capable to address more restricted models. By providing this detailed description of our solution, we expect our distributed tools to be of practical interest for statisticians working in this field.