naranjal naranjal

# MIXANDMIX: numerical techniques for the computation of empirical spectral
distributions of population mixtures

L Cordero-Grande

###### Abstract

The MIXANDMIX (mixtures by Anderson mixing) tool for the computation of the
empirical spectral distribution of random matrices generated by mixtures of
populations is described. Within the population mixture model the mapping
between the population distributions and the limiting spectral distribution
can be obtained by solving a set of systems of non-linear equations, for which
an efficient implementation is provided. The contributions include 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 approximation,
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.

Index terms— large dimensional statistics, random matrix theory, generalized
Marčenko-Pastur equations, asymptotic eigenvalue distribution, numerical
solutions.

##  1 Introduction

Random matrix theory is at the core of modern high dimensional statistical
inference (Yao et al., , 2015) with applications in physics, biology,
economics, communications, computer science 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 dimension $M$ are no longer
valid. A key contribution in this setting is the Marčenko-Pastur theorem
(Marčenko and Pastur, , 1967; Silverstein and Bai, , 1995), which
characterizes the limiting behaviour of the empirical spectral distribution
(ESD) for matrices with random entries when $N\to\infty$. Namely, there exists
a fixed point equation relating the eigenvalues of the empirical and
population distributions, which can be used for inference under the mediation
of appropriate numerical techniques. Using this characterization, subsequent
asymptotics-based inference can be performed, for instance, for dimensionality
reduction, hypothesis testing, signal retrieval, classification or covariance
estimation (Yao et al., , 2015). In addition, Silverstein and Choi,  (1995)
develops on ideas outlined in Marčenko and Pastur,  (1967) to analyze the
support of the ESD 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 limiting behaviour of the
sample eigenvalues is governed by a system of non-linear equations and,
although asymptotic sample eigenvalue confinement has also been proved
(Kammoun and Alouini, , 2016), a simple description of the support is no
longer at hand.

Despite the variety of potential applications, not many works have practically
confronted the numerical issues involved in computing the ESD 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 author shows that the proposed method compares favourably with
straightforward fixed point solvers, both in terms of accuracy and
computational efficiency. In addition, accuracy limitations of Monte Carlo
simulations (Jing et al., , 2010) are showcased and, while recognizing their
independent interest, arguments are provided about the practical limitations
(Rao and Edelman, , 2008) or limited applicability (Olver and Nadakuditi, ,
2013) of other approaches. In our experiments, this package has revealed an
exquisite accuracy and efficiency in computing the ESD and determining its
support. However, the applicability or efficiency of an ODE approach may be
compromised 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 observed eigenvalues 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. 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.

This technical note describes a series of tools that have been developed to
compute the ESD for the mixture of populations case in (Wagner et al., ,
2012). This model or certain analogous forms, has attracted interest in areas
such as telecommunications (Moustakas and Simon, , 2007; Couillet et al., ,
2011; Dupuy and Loubaton, , 2011), machine learning (Benaych-Georges and
Couillet, , 2016; Couillet and Benaych-Georges, , 2016), medical imaging
(Cordero-Grande et al., , 2018) or genetics (Fan and Johnstone, , 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 ESD 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 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 they
are reliable enough and their 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 they are also capable of
providing accurate results 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.1.0. This note is
organized as follows: in § 2 we review different random matrix models, in § 3
we describe the main features of our method, in § 4 we validate the proposed
technique, in § 5 we discuss the implications of the obtained results and in §
6 we end up with some conclusions.

##  2 Theory

Consider a complex random matrix $\mathbf{X}$ of size $N\times M$. We are
interested in the eigenvalue distribution of the sample covariance
$\mathbf{Y}=\displaystyle N^{-1}\mathbf{X}^{H}\mathbf{X}$ in the asymptotic
regime where both $M\to\infty$ and $N\to\infty$ but they maintain a fixed
ratio $\gamma=M/N$ with $\gamma>0$. 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
contemplated:

###  2.1 IID standard entries

Within this model, that we call _standard_ model, we can write
$\mathbf{X}^{H}\sim\boldsymbol{\mathcal{CN}}\displaystyle\left(\mathbf{0}_{MN},\mathbf{I}_{MN}\right)$,
with $\boldsymbol{\mathcal{CN}}$ denoting a circularly symmetric complex
Gaussian distribution, $\mathbf{0}_{MN}$ the $M\times N$ matrix with zero
entries, and $\mathbf{I}_{MN}$ the $MN\times MN$ identity matrix. Marčenko and
Pastur,  (1967) showed that in this case the distribution of eigenvalues of
the sample covariance asymptotically converges to

$f_{\gamma}^{\mathbf{I}}(x)=\begin{cases}\displaystyle\frac{\sqrt{(\gamma_{+}-x)(x-\gamma_{-})}}{2\pi\gamma
x}&\mbox{if }\gamma_{-}\leq x\leq\gamma_{+}\\\ 0&\mbox{otherwise,}\end{cases}$
(1)

with $\gamma_{-}=(1-\sqrt{\gamma})^{2}$ and $\gamma_{+}=(1+\sqrt{\gamma})^{2}$
defining the support of the distribution. Note that if $\gamma\geq 1$, the
distribution has a $1-\gamma^{-1}$ point mass ($\gamma>1$) or is locally
unbounded ($\gamma=1$) at $x=0$. Thus, in this case, (1) gives an explicit
characterization of the ESD.

###  2.2 IID rows

In this scenario, we can write
$\mathbf{X}^{H}\sim\boldsymbol{\mathcal{CN}}\displaystyle\left(\mathbf{0}_{MN},\boldsymbol{\Lambda}_{M}\otimes\mathbf{I}_{N}\right)$,
with $\boldsymbol{\Lambda}_{M}$ a given population covariance matrix. Marčenko
and Pastur,  (1967) and Silverstein and Bai,  (1995) respectively derived
expressions for probabilistic and almost sure limits of the sample covariance
spectrum using its Stieltjes transform

$m(z)=\int_{\mathbb{R}}\frac{f_{\gamma}^{\boldsymbol{\Lambda}_{M}}(x)}{x-z}dx,\quad
z\in\mathbb{C}\setminus\mathbb{R}$ (2)

for which the inversion formula would give back the ESD by

$f_{\gamma}^{\boldsymbol{\Lambda}_{M}}(x)=\frac{1}{\pi}\lim_{\epsilon\to
0^{+}}\Im\left\\{m(x+i\epsilon)\right\\}.$ (3)

If, in a discrete setting, we denote the increasingly sorted eigenvalues of
$\boldsymbol{\Lambda}_{M}$ by $\\{\lambda_{1},\ldots,\lambda_{M}\\}$, when
$N\to\infty$ the eigenvalue distribution of the matrix
$\tilde{\mathbf{Y}}=\displaystyle N^{-1}\mathbf{X}\mathbf{X}^{H}$ converges to
a function whose Stieltjes transform $\tilde{m}(z)$ is related to the
population distribution by the fixed point equation
$\tilde{m}(z)=\displaystyle\left(-z+\frac{\gamma}{M}\sum_{m=1}^{M}\frac{\lambda_{m}}{1+\lambda_{m}\tilde{m}(z)}\right)^{-1}$
and to the Stieltjes transform of the corresponding limiting function for the
spectral distribution of the sample covariance matrix $\mathbf{Y}$, $m(z)$, by
$\tilde{m}(z)=\gamma m(z)+(\gamma-1)/z$. In addition, Silverstein and Choi,
(1995), following the guidelines in Marčenko and Pastur,  (1967), showed that
the support of the ESD can be characterized by studying the zeros of the
derivative of the inverse map, $z^{\prime}(\tilde{m})$, in appropriate
domains. This property, together with the analyticity of the empirical 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)=-\displaystyle\frac{1}{\gamma z\tilde{m}(z)}-\frac{1}{\gamma}$ for an
equivalent fixed point equation,

$e(z)=\frac{1}{M}\sum_{m=1}^{M}\lambda_{m}\left(\frac{\lambda_{m}}{1+\gamma
e(z)}-z\right)^{-1},$ (4)

and an expression for the Stieltjes transform of the ESD,

$m(z)=\frac{1}{M}\sum_{m=1}^{M}\left(\frac{\lambda_{m}}{1+\gamma
e(z)}-z\right)^{-1}.$ (5)

We refer to this model as the _single population_ model.

###  2.3 Independent rows

In this _mixture of populations_ model, the matrix is drawn from
$\mathbf{X}^{H}\mathbin{\sim}\boldsymbol{\mathcal{CN}}\displaystyle\left(\mathbf{0}_{MN},\sum_{k=1}^{K}\boldsymbol{\Lambda}^{k}_{M}\otimes\mathbf{D}^{k}_{N}\right)$,
where $\mathbf{D}^{k}_{N}$ is a diagonal indicator matrix with ones in the
diagonal elements corresponding to those rows sampled according to the
population covariance $\boldsymbol{\Lambda}^{k}_{M}$ and zero otherwise. The
equations relating the empirical and population distributions for $N\to\infty$
have been derived in Wagner et al.,  (2012) also making use of the Stieltjes
transform of the limiting spectral distribution
$f_{\gamma}^{\prescript{}{K}{\boldsymbol{\Lambda}}_{M}}$, $m(z)$, and the
auxiliary functions $e_{j}(z)$, $1\leq j\leq K$. These functions are related
to the population covariances by a system of nonlinear equations

$e_{j}(z)=\frac{1}{M}\operatorname{tr}\left(\boldsymbol{\Lambda}_{j}\left(\sum_{k=1}^{K}\frac{\alpha_{k}\boldsymbol{\Lambda}_{k}}{1+\gamma
e_{k}(z)}-z\mathbf{I}_{M}\right)^{-1}\right),$ (6)

for which there is a unique solution in $\mathbb{C}\setminus\mathbb{R}^{+}$.
$m(z)$ is expressed in terms of these functions as

$m(z)=\frac{1}{M}\operatorname{tr}\left(\left(\sum_{k=1}^{K}\frac{\alpha_{k}\boldsymbol{\Lambda}_{k}}{1+\gamma
e_{k}(z)}-z\mathbf{I}_{M}\right)^{-1}\right),$ (7)

with $\alpha_{k}=\operatorname{tr}(\mathbf{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). 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
penalty in computational efficiency. Thus, we have focused on developing a
reliable method to solve the system (6) not requiring the Jacobian.

##  3 Methods

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

###  3.1 Support detection

When $\gamma\to 0$ the ESD tends to the population distribution for the
standard and single population models, while for the mixture of populations it
is governed by an effective single population distribution
$f_{\gamma}^{\prescript{}{K}{\boldsymbol{\Lambda}}_{M}}=f_{0}^{\overline{\boldsymbol{\Lambda}}_{M}}$
with
$\overline{\boldsymbol{\Lambda}}_{M}=\displaystyle\sum_{k=1}^{K}\alpha_{k}\boldsymbol{\Lambda}_{M}^{k}$,
so that we have
$f_{0}^{\overline{\boldsymbol{\Lambda}}_{M}}(x)=\displaystyle\frac{1}{M}\sum_{m=1}^{M}\delta(x-\overline{\lambda}_{m})$,
with $\overline{\lambda}_{m}$ the $m$-th eigenvalue of
$\overline{\boldsymbol{\Lambda}}_{M}$. Thus, in this limiting case the ESD
support only includes the eigenvalues of an equivalent single population
distribution. On the other side, for $\gamma>0$,
$[a,b]=\left[t^{-1}(1-\sqrt{\gamma})^{2}\displaystyle\min_{k}\lambda_{1}^{k},t(1+\sqrt{\gamma})^{2}\max_{k}\lambda_{M}^{k}\right]$,
with $\lambda_{m}^{k}$ denoting the $m$-th eigenvalue of
$\boldsymbol{\Lambda}_{M}^{k}$ and $t>1$, provide lower and upper bounds on
the limiting support.

To consider these two extreme cases, we group and sort the set of $P=M(K+1)$
eigenvalues $\lambda_{p}=\\{\lambda_{m}^{k},\overline{\lambda_{m}}\\}$,
$1<k<K$, $1<m<M$. Our method attends to the overlap of the set of intervals

$[a_{p},b_{p}]=[t^{-1}(1-\sqrt{\gamma})^{2}\lambda_{p},t(1+\sqrt{\gamma})^{2}\lambda_{p}],$
(8)

i.e., the intervals generated by the maximum possible spectral dispersion of
each grouped eigenvalue, including a numerical safety margin $t$, which we
have set to $t=1.001$ (and constraining $|\gamma-1|\geq 10^{-9}$ to avoid the
potential singularity at $x=0$ when $\gamma=1$).

The support of the limiting distribution will be contained in $[a,b]$, as this
interval considers the maximum possible eigenvalue dispersion as produced by
minorizing ($a$) and majorizing ($b$) point mass population distributions
governed by (1). On the contrary, there is no guarantee for
$\displaystyle\bigcup_{p}[a_{p},b_{p}]$ to contain the whole support of the
limiting distribution. However, the intervals in (8) serve to cover the
$\gamma\to 0$ regime, which originates challenging highly localized support
intervals, as in this case our construction guarantees that the support will
be contained in the defined intervals. Although we have empirically observed
cases where part of the support lies on
$[a,b]\setminus\displaystyle\bigcup_{p}[a_{p},b_{p}]$, we have never observed
that these intervals contain any of the support edges. In this situation, the
adaptive regridding procedure in § 3.2 should generally be enough to
interpolate the distribution estimates on
$[a,b]\setminus\displaystyle\bigcup_{p}[a_{p},b_{p}]$.

Lack of overlap among adjacent intervals is detected by checking the condition
$b_{p}<a_{p+1}$. Assuming this is observed $I\leq P-1$ times, we can get a
partition into $I+1$ segments, each one induced by $P_{i}$ eigenvalues. Each
of these support segments is gridded using a total of
$\max(M^{\text{(o)}}P_{i},M^{\text{(i)}})$ points, with $M^{\text{(i)}}$ the
minimum number of points per segment and $M^{\text{(o)}}\leq M^{\text{(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 the defined grid locations.

The output of this first step is a set of estimates for the distribution in a
non-uniform grid $\mathbf{x}=x_{p}$ with $1\leq p\leq P_{0}$ and
$M^{\text{(o)}}P\leq P_{0}\leq M^{\text{(i)}}P$. The lower and upper bounds on
the number of grid points, $M^{\text{(o)}}P$ and $M^{\text{(i)}}P$, are
respectively attained for large and small enough $\gamma$. Fixing
$M^{\text{(o)}}=3$ and $M^{\text{(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 geometric
multiplicity of the discretized population covariances is given by
$M^{\text{(u)}}$, equivalent expressions of the problem can be obtained for
any $M=SM^{\text{(u)}}$ with $S\in\mathbb{N}_{>0}$. In our experiments $S$ has
been selected by defining a minimum number of grid points to approximate the
densities by $M^{\text{(m)}}$, and making
$M=\max(M^{\text{(u)}},M^{\text{(m)}})$ with $M^{\text{(m)}}=100$.

###  3.2 Adaptive regridding

Additional points are added to the grid by pursuing
$g^{\prime}(x)\sqrt{xf^{\prime\prime}(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^{\prime\prime}(x)$
and the grid value $x$. The first feature, $f^{\prime\prime}(x)$, favours the
allocation of grid points near the support edges, in accordance to the
$\sqrt{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(\mathbf{x})$ is called on the new set of points to allow for an
update of the $f^{\prime\prime}(x)$ values to be used at the next iterative
regridding step. This whole process is repeated $L$ times, so $R_{l}$, $1\leq
l\leq L$ control the final resolution of the distribution computations. The
parameters by default in our implementation are $R_{l}=1$ $\forall l$ and
$L=1$.

###  3.3 Homotopy continuation

The calculation of the ESD 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 $\mathbf{z}=\mathbf{x}+\boldsymbol{\xi}^{2}i$, with
$\boldsymbol{\xi}=\xi^{0}\mathbf{1}_{P_{l}}$ for big enough $\xi^{0}$ common
for all $1\leq p\leq P_{l}$. At each iteration $i$ we compute
$\varepsilon^{i}_{p}=\displaystyle\max_{k}|e^{i}_{k}(x_{p})-e^{i-1}_{k}(x_{p})|$,
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 $\epsilon$
for the distribution $f(x)$, we need to solve the system of equations in a
complex grid given by $x+i\epsilon^{2}$, we can perform the update
$\xi^{j+1}_{p}=\max(\xi^{j}_{p}/\beta,\epsilon)$ whenever
$\varepsilon^{i}_{p}\leq\varepsilon^{i-1}_{p}$. The iterations at the grid
location indexed by $p$ are terminated when the prescribed accuracy is
reached, namely when $\xi^{j}_{p}=\epsilon$ and
$\varepsilon^{i}_{p}<\epsilon$. In our experiments we have used $\xi^{0}=1$
and $\beta=10$, for which we have observed robust and efficient performance.

###  3.4 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\times$ quicker while simultaneously obtaining $1000\times$ 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
$\mathbf{g}(\mathbf{e})$ such as (6), Anderson iterations are computed as

$\mathbf{e}^{i+1}=\mathbf{g}(\mathbf{e}^{i})-\sum_{q=1}^{Q_{i}}(\mathbf{g}(\mathbf{e}^{i-Q_{i}+q})-\mathbf{g}(\mathbf{e}^{i-Q_{i}+q-1}))\nu^{i}_{q},$
(9)

with $Q_{i}$ denoting the number of iterations whose history is used for the
update at iteration $i$ and
$\boldsymbol{\nu}^{i}=(\nu^{i}_{1},\ldots,\nu^{i}_{Q_{i}})^{T}$ obtained by
solving a linear least squares problem involving the fixed point updates
$\mathbf{h}^{i}(\mathbf{e}^{i})=\mathbf{g}(\mathbf{e}^{i})-\mathbf{e}^{i}$ and
their differences
$\Delta\boldsymbol{h}^{i}=\boldsymbol{h}^{i}-\boldsymbol{h}^{i-1}$ arranged in
a $K\times Q_{i}$ matrix
$\Delta\boldsymbol{H}^{i}=[\Delta\boldsymbol{h}^{i-Q_{i}+1},\ldots,\Delta\boldsymbol{h}^{i}]$.
Due to the potential ill-posedness of this system, we have actually solved a
damped version (Scieur et al., , 2019; Henderson and Varadhan, , 2019)

$\boldsymbol{\nu}^{i}=\operatornamewithlimits{argmin}_{\boldsymbol{\nu}}\|\boldsymbol{h}^{i}-\Delta\boldsymbol{H}^{i}\boldsymbol{\nu}\|_{2}^{2}+\lambda^{i}\|\boldsymbol{\nu}\|_{2}^{2},$
(10)

with damping parameter given by $\lambda^{i}=\displaystyle
0.1\max_{k,q}|\Delta 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).

###  3.5 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 ESD computation in the QuEST method (Ledoit and Wolf, ,
2017) seems more plausible as it 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
ESDs, required, for instance, in patch-based image denoising applications
(Cordero-Grande et al., , 2018).

##  4 Results

In this Section we first validate the beneficial effects of Anderson
acceleration and homotopy continuation in ESD computations (§ 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.

###  4.1 Validation of introduced refinements

In Dobriban,  (2015) an experiment was performed illustrating the limitations
of fixed point iterations to obtain accurate results in ESD calculations.
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 $\gamma=0.5$ to
assess the accuracy. Here we repeat 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/\epsilon$
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.

Figure 1: a) Averaged accuracy of different methods $\overline{\Delta
f}=\displaystyle\sum_{p=1}^{P}|\hat{f}_{0.5}^{\mathbf{I}}(x_{p})-f_{0.5}^{\mathbf{I}}(x_{p})|/P$
with $\hat{f}$ denoting the computed density. b) Computation times $t$. c)
Accuracy $\Delta f(x)=|\hat{f}_{0.5}^{\mathbf{I}}(x)-f_{0.5}^{\mathbf{I}}(x)|$
throughout the support (case $\epsilon=10^{-5}$).

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 that
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
analyticity 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 results, $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+\delta
i)$ as we approach the real line with $\delta=\\{1,0.1,0.01,0.001\\}$. First,
for $z=x+1i$ there is a unique minimum in $\mathbb{C}^{+}$ whose basin of
attraction covers the whole of $\mathbb{C}^{+}$. As we decrease $\delta$ (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 $\delta$, the basin of
attraction of this minimum in $\mathbb{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.

Figure 2: a) Calculated ESDs when using 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$.

###  4.2 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 ESD. In Figs. 3a,b, we show
respectively the averaged accuracy and computation times for a set of aspect
ratios ranging from $\gamma=0.05$ to $\gamma=0.95$ in steps of $0.1$,
$\gamma=1$, and reciprocal $1/\gamma$ ranging from $\gamma=0.05$ to
$\gamma=0.95$ in steps of $0.1$ 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
$\gamma=0.5$ case. Analogous plots are provided in Figs. 3d-f for a two-delta
($\delta\delta$) distribution with equiprobable eigenvalues at $\lambda_{1}=1$
and $\lambda_{2}=8$, for which the ESD can be obtained by solving a third
order polynomial equation (Dobriban, , 2015; Rao and Edelman, , 2008). 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 $\epsilon=10^{-6}/\epsilon=10^{-5}$ (MP,
$\gamma<1/\gamma\geq 1$) and $\epsilon=10^{-5}/\epsilon=10^{-4}$
($\delta\delta$, $\gamma<1/\gamma\geq 1$) for SPECTRODE and
$\epsilon=10^{-5}$, $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 $\gamma\to 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 results near the lower edge for the MP case and
throughout the support for the $\delta\delta$ case. The grid sizes used by
each method for $\gamma=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 / $\delta\delta$ problems.

Figure 3: a,d) Averaged accuracy, b,e) computation times, and c,f) accuracy
throughout the domain ($\gamma=0.5$) for a-c) the MP problem and d-f) the
$\delta\delta$ problem with equiprobable $\boldsymbol{\lambda}=\\{1,8\\}$.

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 $\delta\delta$ problems with $\boldsymbol{\lambda}=\\{1,100\\}$ with
respective multiplicity ratios $\mathbf{w}=\\{0.99,0.01\\}$ (Figs. 4a-c), left
skewed (LS) problem, and $\mathbf{w}=\\{0.01,0.99\\}$ (Figs. 4d-f), right
skewed (RS) problem. Note that due to the scale differences in the population
eigenvalue locations, results on the accuracy throughout the domain are more
conveniently represented in a logarithmic scale. SPECTRODE parameters for
approximately matched average computation times have now been modified to
$\epsilon=10^{-4}/\epsilon=10^{-3}$ (LS, $\gamma<1/\gamma\geq 1$) and
$\epsilon=10^{-3}/\epsilon=10^{-2}$ (RS, $\gamma<1/\gamma\geq 1$) while
MIXANDMIX parameters could be kept the same as in previous experiments while
still matching QuEST computation times. Once again, results in Figs. 4a,d show
an improvement of averaged accuracy by several orders of magnitude when using
the MIXANDMIX method. Results on Figs. 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.

Figure 4: a,d) Averaged accuracy, b,e) computation times, and c,f) accuracy
throughout the domain ($\gamma=0.5$) for the $\delta\delta$ problem with
$\boldsymbol{\lambda}=\\{1,100\\}$ skewed towards a-c) the smallest
($\mathbf{w}=\\{0.99,0.01\\}$) and d-f) largest ($\mathbf{w}=\\{0.01,0.99\\}$)
eigenvalue.

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
$\gamma=0.025$, $\gamma=0.5$ and $\gamma=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 $\epsilon=10^{-3}$, the visual
impression is of limited performance. In contrast, MIXANDMIX appears to be
powerful, showing strengths not only in detecting all the support intervals,
but also in picking-up the density oscillations observed for $\gamma=0.025$
and the steeped lower edge for $\gamma=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, $\epsilon=10^{-6}$ for both $\gamma=0.025$
and $\gamma=0.5$, or the maximum code accuracy, $\epsilon=10^{-8}$, for
$\gamma=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 $\gamma=0.025$ and $\gamma=0.5$, and more
plausible functional shape of the former around the left edge for
$\gamma=0.975$, even though approximately $200\times$ less computational
resources were used. Finally, QuEST results show a remarkable ability to
correctly determine the support intervals, but limitations to accurately
approximate the spiked density areas or capture the density oscillations.

Figure 5: ESD for the Comb problem with $100$ equiprobable point masses evenly
distributed in the interval $[0.1,10]$. a) $\gamma=0.025$, b) $\gamma=0.5$ and
c) $\gamma=0.975$.

###  4.3 ESD of population mixtures

In this Section we provide an illustration of the MIXANDMIX results for the
mixture of populations model and the benefits of the GPU architecture. Fig. 6a
shows the calculated ESD for a mixture of populations with $K=6$ equiprobable
populations, so $\alpha_{k}=1/K$, drawn from population covariances
$\boldsymbol{\Lambda}_{M\text{(DIAG)}}^{k}=\Lambda_{mn\text{(DIAG)}}^{k}=((m+k)\bmod
K)+1)\delta[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
$\boldsymbol{\Lambda}_{M\text{(CORR)}}^{k}=\Lambda_{mn\text{(CORR)}}^{k}=\rho^{|m-n|^{l}}\sqrt{\Lambda_{mm\text{(DIAG)}}^{k}\Lambda_{nn\text{(DIAG)}}^{k}}$,
$\rho<1$, $l>0$, so this CORR problem reduces to the DIAG problem when
$l\to\infty$. Namely, Fig. 6b shows the results for $\rho=0.2$, $l=0.25$ and
$\gamma=0.5$, here with $M=120$ for convenience. We can see that calculations
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\geq 2$ the computation times grow with $K$ following a lower than $1$
ratio. 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.

Figure 6: Calculated ESDs and simulations for the a) DIAG and b) CORR
problems. c) computation times for solving the CORR problem for different
values of $\gamma$ and $K$.

##  5 Discussion

We have presented a set of numerical techniques to aid in the computation of
the ESD 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 contemplated even more general models such as the Kronecker model
(Zhang et al., , 2013), with a recent contribution for the separable case in
Leeb,  (2019), 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 potential to reuse or adapt our tools to tackle
them. The population model discussed in this paper admits a free probability
(Mingo and Speicher, , 2017) based formulation given by

$f_{\gamma}^{\prescript{}{K}{\boldsymbol{\Lambda}}_{M}}=\displaystyle\mathop{\vphantom{\bigoplus}\mathchoice{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}}\displaylimits_{k=1}^{K}\alpha_{k}(f_{0}^{\boldsymbol{\Lambda}_{M}^{k}}\boxtimes
f^{\mathbf{I}}_{\gamma/\alpha_{k}}),$ (11)

where $\boxtimes$ represents the free multiplicative convolution and
$\mathop{\vphantom{\bigoplus}\mathchoice{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}{\vbox{\hbox{\leavevmode\resizebox{0.0pt}{}{$\boxplus$}}}}}\displaylimits$
the free additive convolution over the summands. This implies that the
machinery described in Belinschi et al.,  (2017) could be used to obtain the
system of equations in § 2.3. However, the models in Belinschi et al.,  (2017)
are more general, including self-adjoint polynomials not necessarily built
from a combination of Marčenko-Pastur and atomic distributions, but simply
from asymptotically freely independent matrices. Thus, investigations are
required to discern under what conditions our numerical tools can be extended
to cope with these models.

The proposed 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 (although recent results
such as those in Bao et al.,  (2019); Ji,  (2019) may pave the way for them).
Anyhow, the gridding procedure has shown a robust behaviour in challenging
practical scenarios, even when compared with methods that exploit the
properties of the support. 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. Although this grid
refinement criterion has proven effective for all the tested cases, generally
providing finer spectral resolvability than previously proposed methods, other
criteria may be more appropriate for different applications. In this regard,
we should mention that our software is built on top of a generic grid
refinement method that allows to test other possibilities by simply defining
different criteria for grid cell subdivision, with some exemplary alternatives
already included in the code. In addition, support detection correctness can
generally be inferred from the numerical integration of the estimated
measures, which is provided as an output.

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 the tool should be readily usable for many applications.
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^{\text{(i)}}$, robustness of approximation by increasing
$\xi^{0}$ and/or decreasing $\beta$ and accuracy by decreasing $\epsilon$
and/or increasing $L$. However, a line for future research could involve the
incorporation of refined techniques for parameter selection. A possibility
could be to investigate more efficient and robust interplays between homotopy
continuation and nonlinear acceleration, for instance based on adaptive
regularization schemes for nonlinear acceleration (Scieur et al., , 2019).

##  6 Conclusions

This work has introduced a set of techniques to compute the ESD in a mixture
of populations model. A generic procedure using only the functional form of
the fixed point equations relating the population and limiting empirical
distributions has been proposed. Efficient convergence is achieved by Anderson
acceleration and homotopy continuation and novel strategies for grid
construction have been provided. The 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 tool to be of practical interest for
statisticians working in the field.

## Acknowledgments

This work received funding from the European Research Council under the
European Union’s Seventh Framework Programme (FP7/20072013/ERC grant agreement
no. [319456] dHCP project). The research was supported by the Wellcome/EPSRC
Centre for Medical Engineering at King’s College London [WT 203148/Z/16/Z];
the Medical Research Council [MR/K006355/1]; and the National Institute for
Health Research (NIHR) Biomedical Research Centre based at Guy’s and St
Thomas’ NHS Foundation Trust and King’s College London. The views expressed
are those of the author and not necessarily those of the NHS, the NIHR or the
Department of Health. The author also acknowledges the Department of Perinatal
Imaging & Health at King’s College London and advice and support from Edgar
Dobriban, Jo V. Hajnal and Daan Christiaens.

## References

  * Anderson, (1965) Anderson, D. G. (1965).  Iterative procedures for nonlinear integral equations.  J. Assoc. Comput. Mach., 12(4):547–560.
  * Bao et al., (2019) Bao, Z., Erdős, L., and Schnelli, K. (2019).  On the support of the free additive convolution.  J. Anal. Math., in press. https://arxiv.org/abs/1804.11199:1–19.
  * Belinschi et al., (2017) Belinschi, S. T., Mai, T., and Speicher, R. (2017).  Analytic subordination theory of operator-valued free additive convolution and the solution of a general random matrix problem.  J. Reine Ange. Math., 732:21–53.
  * Benaych-Georges and Couillet, (2016) Benaych-Georges, F. and Couillet, R. (2016).  Spectral analysis of the Gram matrix of mixture models.  ESAIM: Probability and Statistics, 20:217–237.
  * Bun et al., (2017) Bun, J., Bouchaud, J.-P., and Potters, M. (2017).  Cleaning large correlation matrices: Tools from random matrix theory.  Phys. Rep., 666:1–109.
  * Cordero-Grande et al., (2018) Cordero-Grande, L., Christiaens, D., Hutter, J., Price, A. N., and Hajnal, J. V. (2018).  Complex diffusion-weighted image estimation via matrix recovery under general noise models.  NeuroImage, in press. https://doi.org/10.1016/j.neuroimage.2019.06.039:1–21.
  * Couillet and Benaych-Georges, (2016) Couillet, R. and Benaych-Georges, F. (2016).  Kernel spectral clustering of large dimensional data.  El. J. Stat., 10:1393–1454.
  * Couillet and Debbah, (2013) Couillet, R. and Debbah, M. (2013).  Signal processing in large systems: A new paradigm.  IEEE Signal Process. Mag., 30(1):24–39.
  * Couillet et al., (2011) Couillet, R., Debbah, M., and Silverstein, J. W. (2011).  A deterministic equivalent for the analysis of correlated MIMO multiple access channels.  IEEE Trans. Inf. Theory, 57(6):3493–3514.
  * Dobriban, (2015) Dobriban, E. (2015).  Efficient computation of limit spectra of sample covariance matrices.  Rand. Matr. Th. Appl., 4(4):1550019:1–36.
  * Dobriban and Owen, (2019) Dobriban, E. and Owen, A. B. (2019).  Deterministic parallel analysis: an improved method for selecting factors and principal components.  J. R. Statist. Soc. B, 81(1):163–183.
  * Dupuy and Loubaton, (2011) Dupuy, F. and Loubaton, P. (2011).  On the capacity achieving covariance matrix for frequency selective MIMO channels using the asymptotic approach.  IEEE Trans. Inf. Theory, 57(9):5737–5753.
  * Fan and Johnstone, (2019) Fan, Z. and Johnstone, I. (2019).  Eigenvalue distributions of variance components estimators in high-dimensional random effects models.  Ann. Stat., in press. https://arxiv.org/abs/1607.02201:1–34.
  * Henderson and Varadhan, (2019) Henderson, N. C. and Varadhan, R. (2019).  Damped Anderson acceleration with restarts and monotonicity control for accelerating EM and EM-like algorithms.  J. Comput. Graph. Stat., in press. https://doi.org/10.1080/10618600.2019.1594835:1–42.
  * Ji, (2019) Ji, H. C. (2019).  Properties of free multiplicative convolution.  arXiv:1903.02326v1, https://arxiv.org/abs/1903.02326:1–31.
  * Jing et al., (2010) Jing, B.-Y., Pan, G., Shao, Q.-M., and Zhou, W. (2010).  Nonparametric estimate of spectral density functions of sample covariance matrices: a first step.  Ann. Stat., 38(6):3724–3750.
  * Kammoun and Alouini, (2016) Kammoun, A. and Alouini, M.-S. (2016).  No eigenvalues outside the limiting support of generally correlated Gaussian matrices.  IEEE Trans. Inf. Theory, 62(7):4312–4326.
  * Ledoit and Wolf, (2017) Ledoit, O. and Wolf, M. (2017).  Numerical implementation of the QuEST function.  Comput. Stat. Data Anal., 115:199–223.
  * Leeb, (2019) Leeb, W. (2019).  Evaluating the boundary and Stieltjes transform of limiting spectral distributions for random matrices with a separable variance profile.  arXiv:1904.11665v1, https://arxiv.org/abs/1904.11665:1–13.
  * Lu et al., (2016) Lu, A.-A., Gao, X., and Xiao, C. (2016).  Free deterministic equivalents for the analysis of MIMO multiple access channel.  IEEE Trans. Inf. Theory, 62(8):4604–4629.
  * Marčenko and Pastur, (1967) Marčenko, V. A. and Pastur, L. A. (1967).  Distribution of eigenvalues for some sets of random matrices.  Math. USSR - Sbornik, 1(4):457–483.
  * Mingo and Speicher, (2017) Mingo, J. A. and Speicher, R. (2017).  Free probability and random matrices.  Springer, New York, USA.
  * Moustakas and Simon, (2007) Moustakas, A. L. and Simon, S. H. (2007).  On the outage capacity of correlated multiple-path MIMO channels.  IEEE Trans. Inf. Theory, 53(11):3887–3903.
  * Nadakuditi, (2014) Nadakuditi, R. R. (2014).  OptShrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage.  IEEE Trans. Inf. Theory, 60(5):3002–3018.
  * Olver and Nadakuditi, (2013) Olver, S. and Nadakuditi, R. R. (2013).  Numerical computation of convolutions in free probability theory.  arXiv:1203.1958v2, https://arxiv.org/abs/1203.1958:1–31.
  * Paul and Aue, (2014) Paul, D. and Aue, A. (2014).  Random matrix theory in statistics: A review.  J. Stat. Plan. Inf., 150:1–29.
  * Ramière and Helfer, (2015) Ramière, I. and Helfer, T. (2015).  Iterative residual-based vector methods to accelerate fixed point iterations.  Comput. Math. Appl., 70(9):2210–2226.
  * Rao and Edelman, (2008) Rao, N. R. and Edelman, A. (2008).  The polynomial method for random matrices.  Found. Comput. Math., 8(6):649–702.
  * Scieur et al., (2019) Scieur, D., d’Aspremont, A., and Bach, F. (2019).  Regularized nonlinear acceleration.  Math. Programmming, in press. https://doi.org/10.1007/s10107-018-1319-8:1–37.
  * Silverstein and Bai, (1995) Silverstein, J. W. and Bai, Z. D. (1995).  On the empirical distribution of eigenvalues of a class of large dimensional random matrices.  J. Multivar. Anal., 54(2):175–192.
  * Silverstein and Choi, (1995) Silverstein, J. W. and Choi, S.-I. (1995).  Analysis of the limiting spectral distribution of large dimensional random matrices.  J. Multivar. Anal., 54(2):295–309.
  * Wagner et al., (2012) Wagner, S., Couillet, R., Debbaha, M., and Slock, D. T. M. (2012).  Large system analysis of linear precoding in correlated MISO broadcast channels under limited feedback.  IEEE Trans. Inf. Theory, 58(7):4509–4537.
  * Wen et al., (2011) Wen, C.-K., Jin, S., and Wong, K.-K. (2011).  On the sum-rate of multiuser MIMO uplink channels with jointly-correlated Rician fading.  IEEE Trans. Comm., 59(10):2883–2895.
  * Yao et al., (2015) Yao, J., Zheng, S., and Bai, Z. (2015).  Large sample covariance matrices and high-dimensional data analysis.  Cambridge University Press, New York, USA.
  * Zhang et al., (2013) Zhang, J., Wen, C.-K., Jin, S., Gao, X., and Wong, K.-K. (2013).  On capacity of large-scale MIMO multiple access channels with distributed sets of correlated antennas.  IEEE J. Sel. Area. Comm., 31(2):133–148.


# On pseudo-absence generation and machine learning for locust breeding ground
prediction in Africa

Ibrahim Salihu Yusuf
InstaDeep &Kale-ab Tessera
InstaDeep &Thomas Tumiel
InstaDeep &Zohra Slim
InstaDeep &Amine Kerkeni
InstaDeep &Sella Nevo
Google Research &Arnu Pretorius
InstaDeep Correspondence<EMAIL_ADDRESS>
###### Abstract

Desert locust outbreaks threaten the food security of a large part of Africa
and have affected the livelihoods of millions of people over the years.
Furthermore, these outbreaks could potentially become more severe and frequent
as a result of global climate change. Machine learning (ML) has been
demonstrated as an effective approach to locust distribution modelling which
could assist in early warning. However, ML requires a significant amount of
labelled data to train. Most publicly available labelled data on locusts are
presence-only data, where only the sightings of locusts being present at a
particular location are recorded. Therefore, prior work using ML have resorted
to pseudo-absence generation methods as a way to circumvent this issue and
build balanced datasets for training. The most commonly used approach is to
randomly sample points in a region of interest while ensuring these sampled
pseudo-absence points are at least a specific distance away from true presence
points. In this paper, we compare this random sampling approach to more
advanced pseudo-absence generation methods, such as environmental profiling
and optimal background extent limitation, specifically for predicting desert
locust breeding grounds in Africa. Interestingly, we find that for the
algorithms we tested, namely logistic regression, gradient boosting, random
forests and MaxEnt, all popular in prior work, the linear logistic model
performed significantly better than the more sophisticated ensemble methods,
both in terms of prediction accuracy and F1 score. Although background extent
limitation combined with random sampling seemed to boost performance for
ensemble methods, no statistically significant differences were detected
between the pseudo-absence generation methods used to train the logistic
model. In light of this, we conclude that simpler approaches such as random
sampling for pseudo-absence generation combined with linear classifiers such
as logistic regression are sensible and effective for predicting locust
breeding grounds across Africa.111Data preprocessing and modelling code:
https://github.com/instadeepai/locust-predict

##  1 Introduction

Climatic conditions leading to cyclones and monsoons, causing heavy rains,
prompted the recent 2019-2021 upsurge in desert locusts (FAO, 2021a). These
upsurges pose a significant threat to food security in affected areas,
especially in the Northern parts of the African continent. Furthermore, the
occurrence and severity of such upsurges could potentially be exacerbated by
global climate change (Vallebona et al., 2008; FAO, 2016; Zhang et al., 2019;
Salih et al., 2020).

The Food and Agriculture Organization (FAO) of the United Nations operate a
sophisticated monitoring and early warning system for locust outbreaks (FAO,
2020). The system relies on a range of technologies serving field survey
operators, control centres and researchers (Cressman, 2008). In particular,
remote sensing has become an invaluable component for early warning because of
its usefulness in locust distribution modelling (Latchininsky and Sivanpillai,
2010; Cressman, 2013; Latchininsky, 2013; Klein et al., 2021). Female locusts
typically lay their eggs in wet, warm soil and wingless nymph locusts,
referred to as hoppers, require specific vegetation nearby to sustain them
before their wings develop (Symmons and Cressman, 2001; FAO, 2021b). This
connection between certain environmental variables and locust behaviour makes
it possible to attempt to model locust distribution through remote sensing
combined with survey data (Piou et al., 2013; Escorihuela et al., 2018; Piou
et al., 2019; Chen et al., 2020; Ellenburg et al., 2021).

Machine learning (ML) has been shown to be a valuable tool for species
distribution modeling (Beery et al., 2021). In particular, many recent papers
have looked at using ML specifically for modelling locusts (Gómez et al.,
2018, 2019; Kimathi et al., 2020; Gómez et al., 2020, 2021). However, even
when remote sensing is capable of providing useful features for such models,
ML still heavily relies on large amounts of labelled data for training.
Currently, the FAO provides many years worth of labelled data on locusts
hosted through their Locust Hub.222https://locust-hub-hqfao.hub.arcgis.com/
This is an extremely useful resource and contains recorded sightings of
locusts in various phases and stages of their lifecycle. That said, survey
teams in general only record the presence of locusts and rarely their absence.
This is typical of many ecological surveys and data of this kind are referred
to as presence-only data. To overcome the lack of negative labels when
training ML models, past work have made use of pseudo-absence generation
(Gómez et al., 2018, 2020, 2021). A commonly used approach is to randomly
sample points in a region of interest while ensuring that pseudo-absence
points are sampled a minimum distance away from any true presence points
(Barbet-Massin et al., 2012). More advanced pseudo-absence generation methods
also exist, such as environmental profiling and background extent limitation
techniques (Iturbide et al., 2015).

In this paper, we compare different pseudo-absence generation methods used in
conjunction with ML, including methods such as random sampling, environmental
profiling and background extent limitation, specifically when modelling the
desert locust species in Africa. We focus on ML algorithms commonly used in
prior work: logistic regression (LR), gradient boosting (XGBoost), random
forests (RF) and maximum entropy (MaxEnt – a presence-background modelling
approach333Presence-background approaches make use of only labelled presence
data and points sampled across the entire study area of interest, referred to
as background data, without using any absence or pseudo-absence data. Even
though MaxEnt does not rely on pseudo-absence generation, we include it in
this work as a strong baseline for comparison.). We train each algorithm on
specific environmental variables combined with presence labels from the FAO’s
Locust Hub. For LR, XGBoost and RF we generate pseudo-absence labels using
each of the above-mentioned pseudo-absence generation methods. Our results
show LR performing significantly better in terms of prediction accuracy and F1
score compared to the ensemble methods XGBoost and RF as well as MaxEnt. For
LR training, there are no statistically significant differences between the
pseudo-absence generation methods, whereas for the ensemble methods, random
sampling combined with background extent limitation significantly improved
performance. We therefore conclude, as well as simultaneously validate, the
approach of random sampling for pseudo-absence generation, used and studied by
Gómez et al. (2018, 2019) and Barbet-Massin et al. (2012) to be effective.
However, in contrast to these works, we find the simpler linear model LR, as
opposed to the RF, to be more suitable when used for breeding ground
prediction over a large region of Africa.

##  2 Methodology

We are interested in testing the effectiveness of different pseudo-absence
generation methods combined with ML for modelling locusts. Here we discuss our
methodology concerning our data, including our choice of environmental
variables and preprocessing. We explain the different pseudo-absence
generation methods and ML algorithms in more detail and formally state our
research hypothesis.

Data. We focus on modelling desert locusts over the entire affected region of
the African continent (we provide the full list of countries in the
supplementary material (SM)). We use the FAO’s Locust Hub observation data in
this study, as it contains the geo-locations of areas where locusts were
observed, type of locust observed and some environmental conditions. It has a
temporal range of 1985 to 2021. The observations were enriched with
environmental data from
NASA444https://disc.gsfc.nasa.gov/datasets/GLDAS_NOAH025_3H_2.1/summary and
ISRIC
SoilGrids555https://data.isric.org/geonetwork/srv/eng/catalog.search#/home
respectively (names and descriptions of each variable are provided in SM). Our
variables include soil characteristics such as moisture, profile (type) and
temperature as well as air pressure, humidity and surface level temperature.
The environmental data from NASA have a temporal range of 2000 to 2021, while
the soil profile information is non-temporal. The three datasets were combined
by selecting a region of temporal overlap from 2000 to 2021. As in Gómez et
al. (2018), we use hopper presence/absence as a proxy for locust breeding
grounds. Given that the maximum time period between the start of egg laying to
the end of the hopper phase is approximately 95 days and that hoppers are not
able to fly, they remain close to the breeding ground and therefore act as a
good proxy. Therefore, we use a 95-day history for all temporal variables and
predict hopper presence 7-days into the future as a proxy for predicting the
presence of potential future breeding grounds. The resulting dataset was split
into train and test sets based on time; the period from 2000-2014 was used for
training, while the period from 2015-2021 was used for testing. Pseudo-absence
generation was also performed separately on each split ensuring a balanced
distribution of hopper presence/absence. After generating pseudo-absences the
train set was further split into a smaller train set and a validation set for
parameter tuning in the ratio of 80:20. To ensure that different algorithm and
pseudo-absence pairs could be tested fairly and on the same test set, we
constructed pseudo-absences for the test set by randomly sampling a balanced
mixture of test pseudo-absence points generated by the different methods
operating on the test set presence points. For more details regarding the
dataset, we refer the reader to the SM.

Figure 1:  Illustrations of pseudo-absence and background data generation.
Blue square is a presence point and red squares represent generated pseudo-
absence points. (a) Random sampling (RS). Orange circle is the exclusion
buffer. (b) Random sampling with environment profiling (RSEP). Dark blue line
shows the boundary between environmentally suitable and unsuitable regions.
(c) Random sampling with background extent limitation (RS+). Light blue lines
show excluded background regions. (d) Random sampling with environment
profiling and background extent limitation (RSEP+). (e) Background sampling
(BS) for presence-background modelling. It is important to note that for BS
the red points do not necessarily correspond to absence but could represent
both presence and absence.

Pseudo-absence generation. We consider four different approaches for pseudo-
absence generation detailed in Iturbide et al. (2015) and depicted in Figure 1
(a)-(d):

  1. 1.

Random sampling (RS): Pseudo-absences are sampled at random across all points
in the study area that are not within some minimum selected distance to any
presence point. The minimum distance between a presence and any absence point
is referred to as the exclusion buffer and is set to 30km for all methods.

  2. 2.

Random sampling with environmental profiling (RSEP): The RSEP method is aimed
at defining the environmental range of the background from which pseudo-
absences are sampled. Environmentally unsuitable areas for locusts are defined
using a presence-only profiling algorithm (one-class SVM) trained on soil
moisture. Once these unsuitable regions have been established, pseudo-absence
points are randomly sampled from within them.

  3. 3.

RS with background extent limitation (RS+): Pseudo-absences are sampled at
random within a limited background extent not including the full study region.
This optimum limited background is determined by a multi-step process as
outlined in (Iturbide et al., 2015).666We perform all pseudo-absence
generation using the mopa R package (Iturbide et al., 2015).

  4. 4.

RSEP with background extent limitation (RSEP+): This method is similar to the
above RS+, but instead of using unconditioned random sampling within the
limited background extent, samples are only within environmentally unsuitable
regions identified through profiling.

As a strong baseline we compare these pseudo-absence generation methods to
only using presence-background data (see Figure 1 (e)) modelling where
background samples are generated over the entire study region of interest
without any constraint on their location with respect to presence points.

ML algorithms. We consider the following algorithms for comparison: logistic
regression (LR), gradient boosting (XGBoost) (Freund and Schapire, 1997;
Friedman, 2001), random forests (RF) (Breiman, 2001) and maximum entropy
(MaxEnt) (Phillips et al., 2006). We provide hyperparameter details in SM.

Hypothesis. Our null hypothesis, $H_{0}$, is that of no difference between the
mean performances of the different pseudo-absence generation methods across
all the algorithms tested including the mean performance of MaxEnt using only
presence-background data. Specifically, let $G=\\{rs,rsep,rs+,rsep+\\}$,
$A=\\{lr,xgboost,rf\\}$ and let $\mu^{g}_{a}$ represent the mean performance
for pseudo-absence generation method $g\in G$, used when training algorithm
$a\in A$. Then the null hypothesis is given by

$\displaystyle
H_{0}:\mu^{g}_{a}=\mu^{g^{\prime}}_{a^{\prime}}=\mu^{BS}_{MaxEnt}\hskip
28.45274pt\forall g,g^{\prime}\in G,a,a^{\prime}\in A|g\neq g^{\prime}\land
a\neq a^{\prime}$ (1)

We are required to reject $H_{0}$ to have evidence that there are indeed
differences between the pseudo-absence generation methods used for training
the different algorithms as well as the presence-background MaxEnt approach.
Our significance level, i.e. the point at which the probability of equal mean
performances under the null hypothesis is low enough to reject the null
hypothesis, is set to $\alpha=0.05$. If $H_{0}$ is rejected, we can continue
with more specific pairwise tests, with appropriate $p$-value adjustments to
control for the family-wise error (Demšar, 2006).

##  3 Experiments

As an example of the different generation methods, in Figure 2, we show the
pseudo-absence points generated by each method in the countries Niger,
Mauritania, Mali, Algeria, Western Sahara and Morocco for November 2003.

Figure 2:  Pseudo-absence and background data generation: an example on a
subset of African countries Niger, Mauritania, Mali, Algeria, Western Sahara
and Morocco for November 2003. (a) Random sampling (RS). (b) Random sampling
with environment profiling (RSEP). White regions indicate environmentally
suitable regions as identified through environmental profiling, i.e. these are
the regions where pseudo-absence points should not be sampled. (c) Random
sampling with background extent limitation (RS+). (d) Random sampling with
environment profiling and background extent limitation (RSEP+). (e) Background
sampling (BS). Note that for BS the red points do not necessarily correspond
to absence but could represent both presence and absence.

Results. The mean accuracy and F1 score over 100 runs for each algorithm and
generation method is shown in Table 1. In general, the performances are
respectable. However, interestingly, the linear model LR is seen to outperform
the more sophisticated ensemble methods, XGBoost and RF, across all generation
methods, as well as MaxEnt, both in terms of absolute mean accuracy and F1
score. Next, we ascertain the statistical significance of these results.

Table 1: Performance comparison between generation methods and ML algorithms.
Bold values indicate top performance across generation methods for a specific
algorithm.

|  | LR | XGBoost | RF | MaxEnt
---|---|---|---|---|---
| RS | $0.8499\pm 0.0021$ | $0.7320\pm 0.0016$ | $0.6623\pm 0.0112$ |
Accuracy | RSEP | $\mathbf{0.8541\pm 0.0020}$ | $0.7020\pm 0.0015$ | $0.6706\pm 0.0133$ | $0.7418\pm 0.0019$
| RS+ | $0.8417\pm 0.0020$ | $\mathbf{0.7580\pm 0.0021}$ | $\mathbf{0.7855\pm 0.0176}$ |
| RSEP+ | $0.8530\pm 0.0022$ | $0.7423\pm 0.0021$ | $0.7712\pm 0.0138$ |
| RS | $0.9064\pm 0.0012$ | $0.8157\pm 0.0009$ | $0.7495\pm 0.0104$ |
F1 | RSEP | $\mathbf{0.9098\pm 0.0011}$ | $0.7895\pm 0.0008$ | $0.7576\pm 0.0121$ | $0.5437\pm 0.0035$
| RS+ | $0.9008\pm 0.0011$ | $\mathbf{0.8448\pm 0.0012}$ | $\mathbf{0.8660\pm 0.0121}$ |
| RSEP+ | $0.9093\pm 0.0013$ | $0.8293\pm 0.0012$ | $0.8527\pm 0.0101$ |

Statistical analysis. We test $H_{0}$ using the Friedman aligned ranks test
(Friedman, 1937) and find that we can easily reject the null hypothesis of no
difference with $p$-value $<2\times 10^{-16}$. Having rejected $H_{0}$, we
perform additional pairwise tests using the $p$-value adjustment procedure for
multiple testing from Holm (1979). For XGBoost and RF significant differences
are detected between the generation methods (all $p$-values $<2\times
10^{-16}$)777We conducted our statistical testing in R using the stats library
where $2\times 10^{-16}$ is considered the minimum $p$-value. Values smaller
than this are indicated with the ‘$<$’ symbol in printed output., providing
evidence that background extent limitation can be of benefit to these ensemble
methods. For LR, the differences between the generation methods are not
significant ($p$-values provided in SM). Given the seemingly absolute
superiority in performance from LR with no difference between generation
methods, we conduct a pairwise test between only the top performing algorithms
and generation pairs and LR using the simplest random sampling. More
specifically, we test the null hypothesis
$H^{a}_{0}:\mu^{rs}_{lr}=\mu^{rs+}_{rf}=\mu^{rs+}_{xgboost}=\mu_{maxent}$. We
find that the difference between LR and the other algorithms to be significant
($p<2\times 10^{-16}$).

Figure 3:  SHAP analysis for logistic regression. An interpretation of the
logistic regression model on the different pseudo-absence generation methods
using Shapley additive explanations (Lundberg and Lee, 2017). (a) Random
sampling (RS). (b) Random sampling with environment profiling (RSEP).

Interpretation. In Figure 3, we provide a SHAP analysis (Lundberg and Lee,
2017) for variable importance of the LR model between using RS or RSEP. We
find that only a few of the total 174 features are important for prediction.
For both methods features such as Albedo_inst_bucket_14 and clay_0.5cm_mean
roughly amount to the same predictive effect as the bottom 165 variables
combined. This is in contrast to the ensemble methods where feature importance
is far more spread out. Therefore, we suspect the LR is performant due to its
ability to avoid overfitting to noisy features.

##  4 Discussion

Our study provides some evidence for the suitability of using random sampling
for pseudo-absence and logistic regression for predicting desert locust
breeding grounds in Africa. We note that a detailed analysis of pseudo-absence
generation appeared in Barbet-Massin et al. (2012), but was performed on
synthetic data, whereas we are primarily interested in locust modeling and
considered more modern methods such as background extent limitation. Finally,
although we find linear models to perform well in our study, recent papers
have started to use deep learning for locust modeling (Samil et al., 2020;
Tabar et al., 2021), moving towards more sophisticated ML model architectures
and leveraging promising new sources of data such as those generated by the
eLocust3m app (PlantVillage, 2020) to good effect. We hope to compare with and
potentially contribute to these approaches in future work.

## Acknowledgments and Disclosure of Funding

This research was supported in part through computational resources provided
by Google.

## References

  * FAO (2021a) FAO, “Desert locust upsurge in 2019–2021,” _http://www.fao.org/ag/locusts/en/info/2094/web18/index.html_ , 2021.
  * Vallebona et al. (2008) C. Vallebona, L. Genesio, A. Crisci, M. Pasqui, A. Di Vecchia, and G. Maracchi, “Large-scale climatic patterns forcing desert locust upsurges in west africa,” _Climate Research_ , vol. 37, no. 1, pp. 35–41, 2008.
  * FAO (2016) FAO, “Weather and desert locusts,” _http://www.fao.org/3/i6152en/i6152en.pdf_ , 2016.
  * Zhang et al. (2019) L. Zhang, M. Lecoq, A. Latchininsky, and D. Hunter, “Locust and grasshopper management,” _Annual review of entomology_ , vol. 64, pp. 15–34, 2019.
  * Salih et al. (2020) A. A. Salih, M. Baraibar, K. K. Mwangi, and G. Artan, “Climate change and locust outbreak in east africa,” _Nature Climate Change_ , vol. 10, no. 7, pp. 584–585, 2020.
  * FAO (2020) FAO, “Early warning system,” _http://www.fao.org/ag/locusts/en/activ/DLIS/earlywarning/index.html_ , 2020\.
  * Cressman (2008) K. Cressman, “The use of new technologies in desert locust early warning,” _Outlooks on Pest Management_ , vol. 19, no. 2, pp. 55–59, 2008.
  * Latchininsky and Sivanpillai (2010) A. V. Latchininsky and R. Sivanpillai, “Locust habitat monitoring and risk assessment using remote sensing and gis technologies,” in _Integrated management of arthropod pests and insect borne diseases_. Springer, 2010, pp. 163–188.
  * Cressman (2013) K. Cressman, “Role of remote sensing in desert locust early warning,” _Journal of Applied Remote Sensing_ , vol. 7, no. 1, p. 075098, 2013.
  * Latchininsky (2013) A. V. Latchininsky, “Locusts and remote sensing: a review,” _Journal of Applied Remote Sensing_ , vol. 7, no. 1, p. 075099, 2013.
  * Klein et al. (2021) I. Klein, N. Oppelt, and C. Kuenzer, “Application of remote sensing data for locust research and management—a review,” _Insects_ , vol. 12, no. 3, p. 233, 2021.
  * Symmons and Cressman (2001) P. Symmons and K. Cressman, “Desert locust guidelines: biology and behaviour,” _FAO, Rome_ , 2001.
  * FAO (2021b) FAO, “Standard operating procedues (sop) desert locust biology and behaviour,” _http://www.fao.org/ag/locusts/common/ecg/359/en/SOPBiologyEv2021.pdf_ , 2021\.
  * Piou et al. (2013) C. Piou, V. Lebourgeois, A. S. Benahi, V. Bonnal, M. el Hacen Jaavar, M. Lecoq, and J.-M. Vassal, “Coupling historical prospection data and a remotely-sensed vegetation index for the preventative control of desert locusts,” _Basic and Applied Ecology_ , vol. 14, no. 7, pp. 593–604, 2013\.
  * Escorihuela et al. (2018) M. J. Escorihuela, O. Merlin, V. Stefan, G. Moyano, O. A. Eweys, M. Zribi, S. Kamara, A. S. Benahi, M. A. B. Ebbe, J. Chihrane _et al._ , “Smos based high resolution soil moisture estimates for desert locust preventive management,” _Remote Sensing Applications: Society and Environment_ , vol. 11, pp. 140–150, 2018.
  * Piou et al. (2019) C. Piou, P.-E. Gay, A. S. Benahi, M. A. O. Babah Ebbe, J. Chihrane, S. Ghaout, S. Cisse, F. Diakite, M. Lazar, K. Cressman _et al._ , “Soil moisture from remote sensing to forecast desert locust presence,” _Journal of Applied Ecology_ , vol. 56, no. 4, pp. 966–975, 2019.
  * Chen et al. (2020) C. Chen, J. Qian, X. Chen, Z. Hu, J. Sun, S. Wei, and K. Xu, “Geographic distribution of desert locusts in africa, asia and europe using multiple sources of remote-sensing data,” _Remote Sensing_ , vol. 12, no. 21, p. 3593, 2020.
  * Ellenburg et al. (2021) W. L. Ellenburg, V. Mishra, J. B. Roberts, A. S. Limaye, J. L. Case, C. B. Blankenship, and K. Cressman, “Detecting desert locust breeding grounds: A satellite-assisted modeling approach,” _Remote Sensing_ , vol. 13, no. 7, p. 1276, 2021.
  * Beery et al. (2021) S. Beery, E. Cole, J. Parker, P. Perona, and K. Winner, “Species distribution modeling for machine learning practitioners: A review,” _arXiv preprint arXiv:2107.10400_ , 2021.
  * Gómez et al. (2018) D. Gómez, P. Salvador, J. Sanz, C. Casanova, D. Taratiel, and J. L. Casanova, “Machine learning approach to locate desert locust breeding areas based on esa cci soil moisture,” _Journal of Applied Remote Sensing_ , vol. 12, no. 3, p. 036011, 2018.
  * Gómez et al. (2019) D. Gómez, P. Salvador, J. Sanz, C. Casanova, D. Taratiel, and J. Casanova, “Desert locust detection using earth observation satellite data in mauritania,” _Journal of Arid Environments_ , vol. 164, pp. 29–37, 2019\.
  * Kimathi et al. (2020) E. Kimathi, H. E. Tonnang, S. Subramanian, K. Cressman, E. M. Abdel-Rahman, M. Tesfayohannes, S. Niassy, B. Torto, T. Dubois, C. M. Tanga _et al._ , “Prediction of breeding regions for the desert locust schistocerca gregaria in east africa,” _Scientific Reports_ , vol. 10, no. 1, pp. 1–10, 2020.
  * Gómez et al. (2020) D. Gómez, P. Salvador, J. Sanz, and J. L. Casanova, “Modelling desert locust presences using 32-year soil moisture data on a large-scale,” _Ecological Indicators_ , vol. 117, p. 106655, 2020.
  * Gómez et al. (2021) D. Gómez, P. Salvador, J. Sanz, J. F. Rodrigo, J. Gil, and J. L. Casanova, “Prediction of desert locust breeding areas using machine learning methods and smos (mir_smnrt2) near real time product,” _Journal of Arid Environments_ , vol. 194, p. 104599, 2021.
  * Barbet-Massin et al. (2012) M. Barbet-Massin, F. Jiguet, C. H. Albert, and W. Thuiller, “Selecting pseudo-absences for species distribution models: how, where and how many?” _Methods in ecology and evolution_ , vol. 3, no. 2, pp. 327–338, 2012.
  * Iturbide et al. (2015) M. Iturbide, J. Bedia, S. Herrera, O. del Hierro, M. Pinto, and J. M. Gutiérrez, “A framework for species distribution modelling with improved pseudo-absence generation,” _Ecological Modelling_ , vol. 312, pp. 166–174, 2015.
  * Freund and Schapire (1997) Y. Freund and R. E. Schapire, “A decision-theoretic generalization of on-line learning and an application to boosting,” _Journal of computer and system sciences_ , vol. 55, no. 1, pp. 119–139, 1997.
  * Friedman (2001) J. H. Friedman, “Greedy function approximation: a gradient boosting machine,” _Annals of statistics_ , pp. 1189–1232, 2001.
  * Breiman (2001) L. Breiman, “Random forests,” _Machine learning_ , vol. 45, no. 1, pp. 5–32, 2001.
  * Phillips et al. (2006) S. J. Phillips, R. P. Anderson, and R. E. Schapire, “Maximum entropy modeling of species geographic distributions,” _Ecological modelling_ , vol. 190, no. 3-4, pp. 231–259, 2006.
  * Demšar (2006) J. Demšar, “Statistical comparisons of classifiers over multiple data sets,” _The Journal of Machine Learning Research_ , vol. 7, pp. 1–30, 2006\.
  * Friedman (1937) M. Friedman, “The use of ranks to avoid the assumption of normality implicit in the analysis of variance,” _Journal of the american statistical association_ , vol. 32, no. 200, pp. 675–701, 1937.
  * Holm (1979) S. Holm, “A simple sequentially rejective multiple test procedure,” _Scandinavian journal of statistics_ , pp. 65–70, 1979.
  * Lundberg and Lee (2017) S. M. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions,” in _Advances in Neural Information Processing Systems 30_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds. Curran Associates, Inc., 2017, pp. 4765–4774. [Online]. Available: http://papers.nips.cc/paper/7062-a-unified-approach-to-interpreting-model-predictions.pdf
  * Samil et al. (2020) H. M. O. A. Samil, A. Martin, A. K. Jain, S. Amin, and S. E. Kahou, “Predicting regional locust swarm distribution with recurrent neural networks,” _arXiv preprint arXiv:2011.14371_ , 2020.
  * Tabar et al. (2021) M. Tabar, J. Gluck, A. Goyal, F. Jiang, D. Morr, A. Kehs, D. Lee, D. P. Hughes, and A. Yadav, “A plan for tackling the locust crisis in east africa: Harnessing spatiotemporal deep models for locust movement forecasting,” in _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, 2021, pp. 3595–3604.
  * PlantVillage (2020) PlantVillage, “elocust3m,” _https://play.google.com/store/apps/details?id=plantvillage.locustsurvey &hl=en_GB&gl=US_, 2020\.
  * Phillips et al. (2017) S. J. Phillips, R. P. Anderson, M. Dudík, R. E. Schapire, and M. E. Blair, “Opening the black box: An open-source release of maxent,” _Ecography_ , vol. 40, no. 7, pp. 887–893, 2017.

## Supplementary Material

Countries in study region. Mauritania, Mali, Egypt, Morocco, Algeria, Sudan,
South Sudan, Niger, Eritrea, Senegal, Libya, Western Sahara, Uganda, Tunisia,
Cape Verde, Chad, Ethiopia, Kenya, Somalia, and Djibouti.

Dataset details – environmental variables, splits, sizes and features.
Temporal variables were extracted from NASA Global Land Data Assimilation
System Version 2.1 (GLDAS-2.1) Noah Land Surface Model with a temporal
resolution of 3 hours and spatial resolution of 0.25 x 0.25 degrees.
https://disc.gsfc.nasa.gov/datasets/GLDAS_NOAH025_3H_2.1/summary

Table 2: Environmental variables and their descriptions Name | Description
---|---
AvgSurfT_inst | Instantaneous average surface skin temperature (K)
Albedo_inst | Instantaneous albedo (%)
SoilMoi0_10cm_inst | Instantaneous soil moisture 0-10cm (kg m-2)
SoilMoi10_40cm_inst | Instantaneous soil moisture 10-40cm (kg m-2)
SoilTMP0_10cm_inst | Instantaneous soil temperature 0-10cm (K)
SoilTMP10_40cm_inst | Instantaneous soil temperature 0-10cm (kg m-2)
Tveg_tavg | 3-hour averaged Transpiration (W m-2)
Wind_f_inst | Instantaneous wind speed (m s-1)
Rainf_f_tavg | 3-hour averaged total precipitation rate (kg m-2 s-1)
Tair_f_inst | Instantaneous air temperature (K)
Qair_f_inst | Instantaneous specific humidity (kg kg-1)
Psurf_f_inst | Instantaneous surface pressure (Pa)

Soil profile variables were extracted from International Soil Reference and
Information Centre (ISRIC) SoilGrids250m 2.0 data product.
https://data.isric.org/geonetwork/srv/eng/catalog.search#/home

Table 3: Environmental variables and their descriptions Name | Description
---|---
sand_0.5cm_mean | Average sand content between 0-5cm (g/kg)
sand_5.15cm_mean | Average sand content between 5-15cm (g/kg)
clay_0.5cm_mean | Average clay content between 0-5cm (g/kg)
clay_5.15cm_mean | Average clay content between 5-15cm (g/kg)
silt_0.5cm_mean | Average silt content between 0-5cm (g/kg)
silt_5.15cm_mean | Average silt content between 5-15cm (g/kg)
Table 4: Dataset splits, sizes and number of features Split | Feature type(s) | Number of features | Presence | Pseudo-Absence | Total
---|---|---|---|---|---
Train | Numeric | 174 | 17007 | 9251 | 26258
Val | Numeric | 174 | 4206 | 2359 | 6565
Test | Numeric | 174 | 5842 | 1238 | 7080

We used a total of 12 temporal and 6 non-temporal variables. For each temporal
variable, we retrieved a 95-day history and removed the last 7 days (including
the observation day), to ensure that we predict 7 days ahead. Futhermore, we
took each 89-day window and bucketised over time by computing the mean value
for every 6-day window. We dropped the last window, which had less than 6
days. After doing this we arrived at 14 windows for each variable, so that in
total we had 168 (14*12) temporal features. In addition to 6 non-temporal
features the total number of features was 174.

ML model hyperparameters. We performed manual hyperparameter tuning on the
validation set resulting in the follow values: a maximum tree depth of 4 for
XGBoost and 15 for RF. The number of random variables used for node splitting
in RF was $\sqrt{p}$, where $p$ is the total number of variables. For our
MaxEnt model, we used a linear feature class and a regularization factor of 1.
These hyperparameters were found using a grid search on the training data over
the available feature classes (linear, quadratic, product, threshold and
hinge) and regularization factor in the range [0.1,1], with a step size of
0.1. We used the MaxNet library from Phillips et al. (2017) for modeling.