# Precise Asymptotics for Spectral Methods
in Mixed Generalized Linear Models

Yihan Zhang Institute of Science and Technology Austria. Email:
<EMAIL_ADDRESS>Marco Mondelli Institute of Science and Technology
Austria. Email<EMAIL_ADDRESS>Ramji Venkataramanan University of
Cambridge. Email<EMAIL_ADDRESS>
###### Abstract

In a mixed generalized linear model, the objective is to learn multiple
signals from unlabeled observations: each sample comes from exactly one
signal, but it is not known which one. We consider the prototypical problem of
estimating two statistically independent signals in a mixed generalized linear
model with Gaussian covariates. Spectral methods are a popular class of
estimators which output the top two eigenvectors of a suitable data-dependent
matrix. However, despite the wide applicability, their design is still
obtained via heuristic considerations, and the number of samples $n$ needed to
guarantee recovery is super-linear in the signal dimension $d$. In this paper,
we develop exact asymptotics on spectral methods in the challenging
proportional regime in which $n,d$ grow large and their ratio converges to a
finite constant. By doing so, we are able to optimize the design of the
spectral method, and combine it with a simple linear estimator, in order to
minimize the estimation error. Our characterization exploits a mix of tools
from random matrices, free probability and the theory of approximate message
passing algorithms. Numerical simulations for mixed linear regression and
phase retrieval display the advantage enabled by our analysis over existing
designs of spectral methods.

###### Contents

  1. 1 Introduction
    1. 1.1 Main contributions
    2. 1.2 Related work
  2. 2 Preliminaries
    1. 2.1 Notation
    2. 2.2 Model
    3. 2.3 Linear estimator
    4. 2.4 Spectral estimator
  3. 3 Main results
    1. 3.1 Linear estimator
    2. 3.2 Spectral estimator
    3. 3.3 Two illustrative examples
    4. 3.4 Technical tools
  4. 4 Numerical experiments
  5. 5 Eigenvalues via random matrix theory
    1. 5.1 Reduction to free additive convolution
    2. 5.2 Right edge of the bulk of $D$
    3. 5.3 Outlier eigenvalues of $D$
  6. 6 Joint distribution via Generalized Approximate Message Passing
    1. 6.1 GAMP as a method to compute the linear and spectral estimators
    2. 6.2 Proof of Theorem 3.1
      1. 6.2.1 Proof of Equation 3.7
      2. 6.2.2 Proof of Equation 3.8
  7. 7 Discussion
  8. A Bayes-optimal combination (proof of Corollary 3.2)
  9. B Additional proofs for linear estimator
    1. B.1 Optimization of overlap (proof of Proposition 3.4)
    2. B.2 Evaluation of overlap (proof of Corollary 3.7)
  10. C Additional proofs for spectral estimator (Proposition 3.6)
    1. C.1 Optimization of spectral threshold
    2. C.2 Optimization of spectral overlap
  11. D Universal lower bound on spectral threshold (missing proof in Remark 3.8)
  12. E Evaluation of spectral threshold and overlap
    1. E.1 Mixed linear regression (proof of Corollary 3.8)
    2. E.2 Mixed phase retrieval (proof of Corollary 3.9)
  13. F Random matrix and free probability theory background
  14. G Theorem 2.1 of [BBCF17]
    1. G.1 Assumptions
    2. G.2 Theorem statement
  15. H State evolution of GAMP for mixed GLMs (proof of Proposition 6.1)
  16. I Auxiliary lemmas

##  1 Introduction

We consider the problem of learning multiple $d$-dimensional vectors from $n$
unlabeled observations coming from a _mixed_ generalized linear model (GLM):

$\displaystyle y_{i}$ $\displaystyle=q\left(\left\langle
a_{i},x_{\eta_{i}}^{*}\right\rangle,\varepsilon_{i}\right),\qquad
i\in[n]=\\{1,\ldots,n\\}.$ (1.1)

Here, $x_{1}^{*},\ldots,x_{\ell}^{*}\in\mathbb{R}^{d}$ are the $\ell$ signals
(regression vectors) to be recovered from the observation vector
$y=(y_{1},\ldots,y_{n})\in\mathbb{R}^{n}$ and the known design matrix
$A=[a_{1},\ldots,a_{n}]^{\top}\in\mathbb{R}^{n\times d}$. For $i\in[n]$,
$\varepsilon_{i}$ is a noise variable, and $\eta_{i}$ is an $[\ell]$-valued
latent variable, i.e., it indicates which signal each observation comes from,
and is unknown to the statistician. The notation
$\left\langle\cdot,\cdot\right\rangle$ denotes the Euclidean inner product,
and $q:\mathbb{R}^{2}\to\mathbb{R}$ is a known link function. For $\ell=1$,
Equation 1.1 reduces to a generalized linear model [MN89], which covers many
widely studied problems in statistical estimation including linear regression,
logistic regression, phase retrieval [SEC+15, FS20], and 1-bit compressed
sensing [BB08]. The regression model with $\ell=1$ implicitly assumes a
homogeneous population, in which a single regression vector suffices to
capture the features of the entire sample. In practice, it is often the case
that the observations may come from multiple sub-populations. Mixed GLMs offer
a flexible solution in settings with unlabeled heterogeneous data, and have
found applications in a variety of fields including biology, physics, and
economics [MP04, GL07, LSL19, DGP20]. When $q(g,\varepsilon)=g+\varepsilon$,
Equation 1.1 reduces to the widely studied mixture of linear regressions
[VT02, FS10, SBvdG10, CL13, YCS14, ZJD16, SS19, ZMCL20, GK20].

A natural approach to estimate the vectors $x_{1}^{*},\ldots,x_{\ell}^{*}$
from $y$ and $A$ is via the maximum-likelihood estimator (assuming a
statistical model for $(\varepsilon_{i})_{i\in[n]}$ is available). However,
the corresponding optimization problem is non-convex and NP-hard [YCS14].
Thus, various low-complexity alternatives — mostly focusing on mixed linear
regression — have been proposed: examples include expectation-maximization
(EM) [KC07, FS10, SBvdG10], alternating minimization [YCS14, SS19, GK20],
convex relaxation [CYC14], moment descent methods [LL18, CLS20], and the use
of tractable non-convex objectives [ZJD16, BH22]. Many of these methods are
iterative in nature and require a “warm start” with an initial guess
correlated with the ground truth. Spectral methods represent a popular way to
provide such initialization [YCS14]. A variety of estimators based on the
spectral decomposition of data-dependent matrices or tensors has been proposed
for mixed GLMs [CL13, YCS14, SJA16]. In this paper, we focus on a spectral
method that estimates the $\ell$ signals via the top-$\ell$ principal
eigenvectors of the following data-dependent matrix:

$\displaystyle
D=\frac{1}{n}\sum_{i=1}^{n}{\mathcal{T}}(y_{i})a_{i}a_{i}^{\top}\in{\mathbb{R}}^{d\times
d},$ (1.2)

where $\mathcal{T}:\mathbb{R}\to\mathbb{R}$ is a suitably chosen preprocessing
function. This spectral estimator with the preprocessing function
${\mathcal{T}}(y)=y^{2}$ was studied for mixed linear regression by Yi et al.
[YCS14], who showed that the signals can be accurately recovered when the
number of observations $n$ is of order $d\log d$. Furthermore, existing
theoretical results for all estimators (including spectral, alternating
minimization and EM) require the number of observations $n$ to be of order at
least $d\log d$ to guarantee accurate recovery [CL13, YCS14, SJA16, LL18,
CLS20]. This leads to the following natural questions:

> _What is the optimal sample complexity of a spectral estimator based on
> Equation 1.2?
> Can we carry out a principled optimization of the preprocessing function
> ${\mathcal{T}}$?_

A simpler alternative to obtain an initial estimate is to employ the linear
estimator

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}{\mathcal{L}}(y_{i})a_{i}\
\in{\mathbb{R}}^{d},$ (1.3)

where ${\mathcal{L}}:{\mathbb{R}}\to{\mathbb{R}}$ is a suitable preprocessing
function. The performance analysis of this linear estimator for the mixed GLM
can be carried out similarly to that for the non-mixed case ($\ell=1$); the
analysis for the latter is given in [PVY17, Proposition 1] and in [MTV21,
Lemma 2.1]. Thus, a second natural question is:

> _What is the optimal way to combine a spectral estimator based on Equation
> 1.2 and the linear estimator in Equation 1.3?_

###  1.1 Main contributions

In this paper, we resolve the questions above for the prototypical setting of
the recovery of two independent signals $x_{1}^{*},x_{2}^{*}$ with a Gaussian
design matrix $A$. This is achieved by characterizing the high-dimensional
limit of the joint empirical distribution of _(i)_ the signals
$x_{1}^{*},x_{2}^{*}$, _(ii)_ the linear estimator in Equation 1.3, and
_(iii)_ spectral estimators based on the matrix in Equation 1.2. Our analysis
holds in the proportional setting where $n,d\to\infty$ with
$n/d\to\delta\in(0,\infty)$. That is, we consider the regime where the ratio
between sample size and signal dimension tends to a constant, as opposed to
most analyses of mixed GLMs in the literature which assume $n=\Omega(d\log
d)$. Our major findings are summarized as follows.

  * •

Our master theorem (Theorem 3.1) characterizes the joint distribution of the
linear estimator, the spectral estimator, and the signals in the high-
dimensional limit. This joint distribution characterization holds for
arbitrary preprocessing functions
${\mathcal{L}},{\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$ in Equations
1.2 and 1.3 (subject to certain mild regularity conditions). The limiting
joint distribution is expressed as the law of a set of jointly Gaussian random
variables whose covariance structure is explicitly derived in terms of the
model and the preprocessing functions.

  * •

As an immediate consequence of the above characterization, we derive the
normalized correlations (or ‘overlaps’) between the linear/spectral estimator
and the signals (Corollary 3.3/Corollary 3.5). The linear estimator achieves a
strictly positive overlap with each signal for any $\delta>0$, provided a
strictly positive overlap can be attained for some $\delta>0$. In contrast,
for the spectral estimator, we identify a threshold (depending on the
preprocessing function ${\mathcal{T}}$) such that strictly positive overlap is
attained as soon as $\delta$ exceeds this threshold.

  * •

Our master theorem also allows us to compute the limiting overlap of a class
of combinations of linear and spectral estimators. In particular, the Bayes-
optimal combination can be derived, which turns out to be a linear combination
of the two estimators due to the Gaussianity of their high-dimensional limits
(Corollary 3.2).

  * •

We determine the optimal preprocessing functions
${\mathcal{L}}^{*},{\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}\colon{\mathbb{R}}\to{\mathbb{R}}$
for the linear and spectral estimators that maximize the overlap between the
estimator and each signal (Proposition 3.4 and Proposition 3.6). The optimal
overlaps of linear and spectral estimators reveal intriguing behaviors of
mixed models. In particular, there is a _single_ function ${\mathcal{L}}^{*}$
that simultaneously maximizes the overlap between the linear estimator and
each signal. In contrast, for the spectral method, one needs to employ two
_different_ functions ${\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}$ in order
to achieve the maximal overlaps with $x_{1}^{*}$, $x_{2}^{*}$, respectively.
Furthermore, the optimal overlap of the spectral estimator with each signal
approaches $1$ — the best possible value — as the aspect ratio $\delta$ grows.
We remark that the same is not true for the linear estimator: the optimal
overlap with each signal remains strictly less than 1 even as
$\delta\to\infty$, as long as there is a strictly positive fraction of
observations corresponding to each signal.

  * •

Finally, we specialize our results to two canonical settings: mixed linear
regression (Corollaries 3.7 and 3.8) and mixed phase retrieval (Corollary
3.9). In the noiseless case, even though the two models are distinct111The
phase retrieval model, by definition, does not preserve the sign of the linear
observations, whereas the linear regression model preserves it., spectral
estimators have the same performance on both. In contrast, the linear
estimator exhibits non-trivial performance for linear regression, whereas it
always results in vanishing overlap for phase retrieval.

Our precise asymptotic analysis leads to a significant improvement over
previous designs of spectral methods, as showcased in Figure 1 for noiseless
mixed linear regression. The continuous lines correspond to our theoretical
predictions (“pred.”), which closely match the points coming from the
simulations (“sim.”). The following methods are compared: _(i)_ optimal
spectral method (black), obtained from Proposition 3.6; _(ii)_ optimal linear
method (blue), obtained from Proposition 3.4; _(iii)_ combined estimator
(“combo”) (red), obtained from Corollary 3.2; _(iv)_ spectral estimator for
mixed linear regression proposed in [YCS14] (yellow); _(v)_ spectral estimator
which optimizes the overlap in the non-mixed setting (green), proposed in
[LAL19]. The spectral methods resulting from our sharp analysis (red, black)
significantly outperform existing methods (green, yellow), especially for low
values of $\delta$. More details on the experimental setup and additional
simulation results can be found in Section 4.

(a) Recovery of $x_{1}^{*}$

(b) Recovery of $x_{2}^{*}$

Figure 1: Noiseless mixed linear regression with mixing parameter (i.e.,
probability that a sample corresponds to $x_{1}^{*}$) $\alpha=0.6$. Overlaps
with the first signal $x_{1}^{*}$ (left) and the second signal $x_{2}^{*}$
(right), computed via simulation (“sim.”) and the theoretical prediction
(“pred.”), are plotted as a function of the aspect ratio $\delta=n/d$. The
signal dimension is $d=2000$.

##### Proof techniques.

We exploit a combination of tools from free probability, random matrices and
the theory of approximate message passing. More specifically, generalized
approximate message passing (GAMP) refers to a family of iterative algorithms
[Ran11] with the following key feature: the joint distribution of the iterates
is accurately tracked by a simple deterministic recursion, called _state
evolution_. Our strategy to obtain the joint distribution of the
linear/spectral estimators and the signals in the master theorem (Theorem 3.1)
is to design a GAMP that _(i)_ outputs the linear estimator as the first
iterate and _(ii)_ then implements a power method, so that its fixed point
corresponds to the spectral estimator. One challenge in the implementation of
this strategy is that the state evolution of GAMP, in its original form for
vanilla (non-mixed) GLMs, only records the correlation of its iterates with a
single signal. To circumvent this issue, we equip GAMP with a state evolution
recursion involving both signals, and run a _pair_ of GAMP iterations
converging to the first and second top eigenvector of the spectral matrix in
Equation 1.2, respectively. A second – even more fundamental – challenge is
that, for the power method to converge to the desired eigenvector, a spectral
gap between the corresponding eigenvalue and the rest of the spectrum is
required. For non-mixed GLMs, the spectral analysis was carried out in earlier
work [LL20, MM19], and _both_ eigenvalues _and_ overlaps were characterized
using tools from random matrix theory. Here, the difficulty comes from the
mixed effect of the model, leading to additional matrix terms which appear to
be difficult to control. Our approach is to decompose $D$ into the sum of two
matrices, $D_{1}$ and $D_{2}$, each consisting of components only pertaining
to the first and second signal, respectively. Now, $D_{1},D_{2}$ can be
individually viewed as generated from a non-mixed GLM, hence their limiting
spectra are well understood. The key observation is then that, by assuming
both signals to be independent and uniformly distributed on the sphere,
$D_{1}$ and $D_{2}$ become _asymptotically free_ 222Asymptotic freeness can be
thought of as the random matrix analogue of _independence_ of random
variables.. Thus, we are able to characterize the sum of these two spiked
matrices by using techniques from free probability.

###  1.2 Related work

We review below the three lines of literature most closely related to our
work, which concern mixtures of generalized linear models, spectral methods
and the theory of approximate message passing algorithms.

##### Mixtures of generalized linear models.

Mixtures of generalized linear models have been studied in machine learning
under the name ‘hierarchical mixtures of experts’, see e.g., [JJ94]. Bayesian
methods for inference in this model were investigated in [PJT96] and [WMR95],
and Bayesian inference for the special case of mixed linear regression (MLR)
was analyzed in [VT02].

Khalili and Chen [KC07] proposed a penalized likelihood approach for variable
selection in mixed GLMs, and showed consistency of the procedure in the low-
dimensional setting (where the dimension $d$ is fixed as $n$ grows). Städler
et al. [SBvdG10] analyzed $\ell_{1}$-penalized estimators for high-dimensional
mixed linear regression. Recently Zhang et al. [ZMCL20] studied estimation and
inference for high-dimensional mixed linear regression with two sparse
components, in the setting where the mixing proportion and the covariance
structure of the covariates are unknown. The works [KC07, SBvdG10, ZMCL20] all
use variants of the EM algorithm for optimizing a suitable penalized
likelihood function.

Balakrishnan et al. [BWY17] and Klusowski et al. [KYB19] obtained statistical
guarantees on the performance of the EM algorithm for a class of problems,
including the special case of symmetric mixed linear regression (MLR) where
$x_{1}^{*}=-x_{2}^{*}$. Variants of the EM algorithm for symmetric MLR in the
high-dimensional setting (with sparse signals) were analyzed in [WGNL15, YC15,
ZWZG17]. Minimax lower bounds for a class of computationally feasible
algorithms for symmetric MLR were obtained in [FLWY18]. Kong et al. [KSS+20]
studied MLR as a canonical example of meta-learning: in the setting where the
number of signals ($\ell$) is large, they derived conditions under which a
large number of signals with a few observations can compensate for the lack of
signals with abundantly many observations. The prediction error of MLR in the
non-realizable setting, where no generative model is assumed for the data, was
studied in [PMSG22]. Chandrasekher et al. [CPT21] recently analyzed the
performance of a class of iterative algorithms (not including AMP) for
mixtures of GLMs. They provide a sharp characterization of the per-iteration
error with sample-splitting in the regime $n\asymp d\,\mathrm{polylog}(d)$,
assuming a Gaussian design and a random initialization.

##### Spectral methods.

We first review spectral methods based on Equation 1.2 for standard GLMs (non-
mixed, with $\ell=1$), which were introduced in [Li92]. For the special case
of phase retrieval, a series of works has provided increasingly refined bounds
on the number of samples needed to guarantee signal recovery via the spectral
method [NJS13, CSV13, CC15]. This type of analysis is based on matrix
concentration inequalities, a technique that typically does not return exact
values for the overlap between the signal and the estimate. More recently, an
exact high-dimensional analysis for generalized linear models was carried out
in [LL20, MM19]. These works focus on the regime of interest in this paper:
$n$ and $d$ growing at a proportional rate $\delta\coloneqq\lim\frac{n}{d}$.
This sharp analysis allows for the optimization of the preprocessing function:
the choice of ${\mathcal{T}}$ minimizing the value of $\delta$ (and, hence,
the amount of data) needed to achieve a strictly positive overlap was provided
in [MM19]; furthermore, the choice of ${\mathcal{T}}$ maximizing the overlap
was provided in [LAL19]. The analysis of the spectral methods in the works
above assumes a Gaussian design matrix. Going beyond this assumption, [DBMM20]
provides precise asymptotics for design matrices sampled from the Haar
distribution, and [MKLZ22] studies rotationally invariant designs.

Moving to the mixed regression setting ($\ell>1$), Yi et al. [YCS14] proposed
a spectral estimator based on Equation 1.2 with ${\mathcal{T}}(y)=y^{2}$. The
analysis is based on concentration inequalities and requires the number of
samples $n$ to be of order $d\log d$ for accurate recovery. Estimators based
on spectral decomposition of data-dependent tensors were proposed for MLR in
[CL13] and for mixed GLMs in [SJA16]. However, these methods require $n$ to be
of order at least $d^{3}$ for accurate recovery. Our work is the first to
establish exact asymptotics for a mixed GLM in the linear sample-size regime:
$n,d\to\infty$ with $n/d\to\delta\in(0,\infty)$. To achieve this goal, our
strategy differs from analyses of spectral methods in the non-mixed setting
[LL20, MM19] which reduces the study of the spectrum of $D$ to that of a
rank-1 perturbation. In contrast, our analysis is based on a combination of
techniques from free probability and approximate message passing.

##### Approximate message passing (AMP) algorithms.

AMP is a family of iterative algorithms that has been applied to several
problems in high-dimensional statistics, including estimation in linear models
[DMM09, BM11, KMS+12], generalized linear models [Ran11, SR14, SC19], and low-
rank matrix estimation [DM14, RFG09, LKZ17]. For a broad overview, we refer
the reader to [FVRS22]. A key feature of AMP algorithms is that under suitable
model assumptions, the empirical joint distribution of their iterates can be
exactly characterized in the high-dimensional limit, in terms of a simple
scalar recursion called _state evolution_. By taking advantage of this
characterization, AMP methods have been used to derive exact high-dimensional
asymptotics for convex penalized estimators such as LASSO [BM12], M-estimators
[DM16], logistic regression [SC19], and SLOPE [BKRS20]. AMP algorithms have
been initialized via spectral methods in the context of low-rank matrix
estimation [MV21c] and generalized linear models [MV21a]. Furthermore, they
have been used – in a non-mixed setting – to combine linear and spectral
estimators [MTV21]. A finite-sample analysis which allows the number of
iterations to grow roughly as $\log n$ ($n$ being the ambient dimension) was
put forward in [RV18], and the recent paper [LW22] improves this guarantee to
a linear (in $n$) number of iterations. This could potentially allow to study
settings in which $\delta=n/d$ approaches the spectral threshold. The works on
AMP discussed above all assume i.i.d. Gaussian matrices. A number of recent
papers has proposed generalizations of AMP for the much broader class of
rotationally invariant matrices, e.g., [OCW16, MP17, RSF19, Tak20, ZSF21,
Fan22, MV21b, VKM22].

##### Organization of the paper.

Formal definitions for the mixed GLM, linear estimators and spectral
estimators are presented in Section 2. Our main results are stated in Section
3: these include the master theorem (Theorem 3.1) and its various consequences
(whose proofs are deferred to Appendices A, B, C, D and E). Numerical
simulations are provided in Section 4. The proof of Theorem 3.1 is divided
across two sections. Section 5 contains a characterization of the top three
eigenvalues of the matrix $D$ which will be useful in the analysis of GAMP in
the following section. The limiting joint law of the signal, the linear and
the spectral estimators in Theorem 3.1 is then proved in Section 6 using a
GAMP algorithm and its characterization via state evolution. The proof of the
state evolution characterization is deferred to Appendix H. The main body of
the paper is concluded with discussions in Section 7. Some background on
random matrix and free probability theory is provided in Appendices F and G,
and several auxiliary lemmas used in our proofs are listed in Appendix I.

##  2 Preliminaries

###  2.1 Notation

The $i$-th element in a vector $a\in{\mathbb{R}}^{p}$ is denoted by
$a_{i}\in{\mathbb{R}}$. If a vector has multiple subscripts, the component
index is always the last one. For a symmetric matrix
$M\in{\mathbb{R}}^{p\times p}$, we denote by $\mu_{M}$ the empirical spectral
distribution of $M$. The (real) eigenvalues of $M$ are denoted by
$\lambda_{1}(M)\geq\lambda_{2}(M)\geq\cdots\geq\lambda_{p}(M)$, and the
corresponding eigenvectors of unit norms are denoted by
$v_{1}(M),v_{2}(M),\cdots,v_{p}(M)$. The $(i,j)$-th entry of $M$ is denoted by
$M_{i,j}$. For a random variable $X$, we use $\operatorname{supp}(X)$ to
denote the support of its density function. For a set
${\mathcal{S}}\subset{\mathbb{R}}$ of real numbers, denote by
$\sup({\mathcal{S}})\coloneqq\sup\\{x:x\in{\mathcal{S}}\\}$ its least upper
bound. The orthogonal group in dimension $p$ is denoted by
${\mathbb{O}}(p)\coloneqq\left\\{O\in{\mathbb{R}}^{p\times
p}:OO^{\top}=O^{\top}O=I_{p}\right\\}$. The unit sphere in dimension $p$ is
denoted by
${\mathbb{S}}^{p-1}\coloneqq\left\\{x\in{\mathbb{R}}^{p}:\left\|x\right\|_{2}=1\right\\}$.
For two distributions $P$ and $Q$, we use $P\otimes Q$ to denote the product
distribution with $P$ (resp. $Q$) being its first (resp. second) marginal. For
an integer $k\geq 2$, $P^{\otimes k}$ denotes the $k$-fold product
distribution of $P$.

###  2.2 Model

We consider a two-component mixed generalized linear model (mixed GLM) with
signal vectors $x_{1}^{*},x_{2}^{*}\in{\mathbb{S}}^{d-1}$, covariate vectors
$a_{1},\cdots,a_{n}\in{\mathbb{R}}^{d}$ , and a known link function
$q\colon{\mathbb{R}}^{2}\to{\mathbb{R}}$. Let $P_{\varepsilon}$ be a noise
distribution over ${\mathbb{R}}$. The $n$ observations
$y_{1},\cdots,y_{n}\in{\mathbb{R}}$ are generated as:

$\displaystyle y_{i}$ $\displaystyle=q\left(\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle,\,\varepsilon_{i}\right),\qquad
i\in[n].$ (2.1)

Here, the vector of latent variables
${\underline{\eta}}\coloneqq(\eta_{1},\cdots,\eta_{n})\sim\operatorname{Bern}(\alpha)^{\otimes
n}$ indicates which signal is selected by each observation, and is
_unobserved_. The vector ${\underline{\eta}}$ is independent of the signals
$x_{1}^{*},x_{2}^{*}$ and the covariate vectors $a_{1},\ldots,a_{n}$. The
noise vector
${\underline{\varepsilon}}\coloneqq(\varepsilon_{1},\cdots,\varepsilon_{n})\sim
P_{\varepsilon}^{\otimes n}$ is independent of
$x_{1}^{*},x_{2}^{*},{\underline{\eta}}$, in which case Equation 2.1 becomes

$y_{i}\,|\,\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle\sim
p(\cdot\,|\,\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle),$ (2.2)

where $p(\cdot|g)$ denotes the distribution of $q(g,\varepsilon)$ for a fixed
$g\in\mathbb{R}$ and $\varepsilon\sim P_{\varepsilon}$ independent of $g$. The
design matrix $A\in{\mathbb{R}}^{n\times d}$ is formed by collecting all
$(a_{i})_{i\in[n]}$ as rows:

$\displaystyle A$ $\displaystyle\coloneqq\begin{bmatrix}a_{1}^{\top}\\\
\vdots\\\ a_{n}^{\top}\end{bmatrix}\in{\mathbb{R}}^{n\times d}.$

Given $A$, upon observing $y=(y_{1},\cdots,y_{n})\in{\mathbb{R}}^{n}$, our
goal is to estimate $x_{1}^{*}$ and $x_{2}^{*}$. Given a pair of estimators
$\widehat{x}_{1}=\widehat{x}_{1}(y,A),\,\widehat{x}_{2}=\widehat{x}_{2}(y,A)\in{\mathbb{R}}^{d}$,
a performance measure of central interest is their _overlap_ with the
respective signals:

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle\widehat{x}_{1},x_{1}^{*}\right\rangle\right|}{\left\|\widehat{x}_{1}\right\|_{2}\left\|x_{1}^{*}\right\|_{2}},\quad\lim_{d\to\infty}\frac{\left|\left\langle\widehat{x}_{2},x_{2}^{*}\right\rangle\right|}{\left\|\widehat{x}_{2}\right\|_{2}\left\|x_{2}^{*}\right\|_{2}}.$

Throughout the paper, the following assumptions are imposed.

  1. (A1)

$x_{1}^{*},x_{2}^{*}$ are independent and uniformly distributed on unit
sphere,
$(x_{1}^{*},x_{2}^{*})\sim\operatorname{Unif}({\mathbb{S}}^{d-1})^{\otimes
2}$.

  2. (A2)

$\alpha\in(1/2,1)$.

  3. (A3)

The noise sequence ${\underline{\varepsilon}}\in{\mathbb{R}}^{n}$ is i.i.d.
according to ${\underline{\varepsilon}}\sim P_{\varepsilon}^{\otimes n}$, and
$P_{\varepsilon}$ has bounded second moment.

  4. (A4)

$a_{1},\cdots,a_{n}\in{\mathbb{R}}^{d}$ are i.i.d., each distributed according
to $a_{i}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0_{d},I_{d})$.

  5. (A5)

We consider the proportional regime where $n,d\to\infty$ and $n/d\to\delta$
for some constant $\delta\geq 0$ which we call _aspect ratio_.

As for Item (A1), choosing signals uniformly distributed on the sphere
corresponds to a setting in which no structural information about them is
available (and, therefore, the uniform prior is selected by the statistician).
We note that this requirement is natural, since our focus is on spectral
methods which are typically unable to exploit prior information about the
signal. Understanding the effect of correlation on the design of spectral
methods is an exciting avenue for future research.

Item (A2), which implies that a larger fraction of the $n$ observations come
from $x_{1}^{*}$ than from $x_{2}^{*}$, is without loss of generality.
Otherwise, if $0<\alpha<1/2$, one can simply interchange the roles of
$x_{1}^{*}$ and $x_{2}^{*}$. The case $\alpha=1/2$ is special. In this case,
as $n\to\infty$, the top two eigenvectors given by the spectral method
correspond to the same limiting eigenvalue. These eigenvectors provide an
estimate on the space spanned by $x_{1}^{*},x_{2}^{*}$ and, in order to
estimate the individual signals, an additional $1$-dimensional grid search has
to be performed. Provided this extra step is carried out, our results can be
shown to extend to the case $\alpha=1/2$. See Remarks 3.4, 5.2 and 6.2 for
more details.

Item (A4) is common in theoretical analyses of spectral estimators for
standard (non-mixed) GLMs [YCS14, MM19, LL20, LAL19]. The potential
universality of the predictions obtained with Gaussian design matrices is
discussed in Section 7.

###  2.3 Linear estimator

Let ${\mathcal{L}}\colon{\mathbb{R}}\to{\mathbb{R}}$ be a preprocessing
function. Then, the _linear estimator_ is defined as

$\displaystyle\widehat{x}^{\mathrm{lin}}$
$\displaystyle\coloneqq\frac{1}{n}A^{\top}{\mathcal{L}}(y)=\frac{1}{n}\sum_{i=1}^{n}{\mathcal{L}}(y_{i})a_{i}\in{\mathbb{R}}^{d},$
(2.3)

where ${\mathcal{L}}$ is applied component-wise, i.e.,
${\mathcal{L}}(y)=({\mathcal{L}}(y_{1}),\cdots,{\mathcal{L}}(y_{n}))$. Define
the random variable $Y$ as

$\displaystyle Y=q(G,\varepsilon),\quad\text{ where
}(G,\varepsilon)\sim{\mathcal{N}}(0,1)\otimes P_{\varepsilon}.$ (2.4)

We make the following assumption on the preprocessing function ${\mathcal{L}}$
used in the linear estimator.

  1. (A6)

${\mathcal{L}}\colon{\mathbb{R}}\to{\mathbb{R}}$ is Lipschitz and satisfies
the following conditions:

$\displaystyle\mathbb{E}\left[G{\mathcal{L}}(Y)\right]$ $\displaystyle\neq
0,\quad\mathbb{E}\left[|G{\mathcal{L}}(Y)|\right]<\infty.$

As shall be seen in our main result (Theorem 3.1), the first condition in Item
(A6) guarantees that the linear method w.r.t. ${\mathcal{L}}$ attains positive
overlaps with both signals. The second condition is rather mild and purely
technical.

###  2.4 Spectral estimator

Let ${\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$ be a preprocessing
function, and consider

$\displaystyle T$
$\displaystyle\coloneqq\mathrm{diag}({\mathcal{T}}(y))\in{\mathbb{R}}^{n\times
n},\quad
D\coloneqq\frac{1}{n}A^{\top}TA=\frac{1}{n}\sum_{i=1}^{n}{\mathcal{T}}(y_{i})a_{i}a_{i}^{\top}\in{\mathbb{R}}^{d\times
d},$ (2.5)

where we use again the notation
${\mathcal{T}}(y)=({\mathcal{T}}(y_{1}),\cdots,{\mathcal{T}}(y_{n}))$. Then,
the spectral method computes the top two eigenvectors $v_{1}(D),v_{2}(D)$ of
$D$ as estimates of $x_{1}^{*},x_{2}^{*}$, respectively. We make the following
assumption on the preprocessing function ${\mathcal{T}}$ used in the spectral
estimator.

  1. (A7)

${\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$ is Lipschitz and satisfies

$\displaystyle\inf_{y\in\operatorname{supp}(Y)}{\mathcal{T}}(y)>-\infty,\quad\text{and}\quad
0<\sup_{y\in\operatorname{supp}(Y)}{\mathcal{T}}(y)<\infty,$

where $Y$ is defined in Equation 2.4. We also assume that the random variable
${\mathcal{T}}(Y)$ is not almost surely zero, i.e.,
$\Pr\left[{\mathcal{T}}(Y)=0\right]<1$.

In words, Item (A7) requires ${\mathcal{T}}$ to be bounded, with the upper
edge of its range being strictly positive. Having a bounded preprocessing
function is necessary for the spectral method to be effective in the non-mixed
setting as well [MM19, LL20]. Furthermore, the requirement on the $\sup$ to be
strictly positive is purely technical, and is simply to rule out the trivial
cases in which the spectral matrix $D$ is all-zero with high probability. We
also note that Item (A7) is satisfied by the preprocessing function that
maximizes the overlap (cf. Proposition 3.6).

##  3 Main results

We start by defining a few auxiliary quantities. Let
$\delta_{1}=\alpha\delta,\delta_{2}=(1-\alpha)\delta$, and
$Z={\mathcal{T}}(Y)$, with $Y$ as defined in Equation 2.4. Define
$\varphi\colon(\sup\operatorname{supp}(Z),\infty)\to{\mathbb{R}}$ and
$\psi\colon(\sup\operatorname{supp}(Z),\infty)\times(0,\infty)\to{\mathbb{R}}$
as

$\displaystyle\varphi(\lambda)$
$\displaystyle\coloneqq\lambda\,\mathbb{E}\left[\frac{ZG^{2}}{\lambda-Z}\right],$
(3.1) $\displaystyle\psi(\lambda;\Delta)$
$\displaystyle\coloneqq\lambda\left(\frac{1}{\Delta}+\mathbb{E}\left[\frac{Z}{\lambda-Z}\right]\right).$
(3.2)

In what follows, we will set the second argument $\Delta$ of $\psi$ to
$\delta,\delta_{1}$ and $\delta_{2}$. For
$\Delta\in\\{\delta,\delta_{1},\delta_{2}\\}$, let
$\overline{\lambda}(\Delta)>\sup\operatorname{supp}(Z)$ be the minimum point
of $\psi(\cdot;\Delta)$, i.e.,

$\displaystyle\overline{\lambda}(\Delta)$
$\displaystyle\coloneqq\operatorname*{argmin}_{\lambda>\sup\operatorname{supp}(Z)}\psi(\lambda;\Delta).$
(3.3)

Since $\psi$ is convex in its first argument (see Lemma I.2), this minimum
point is readily obtained by setting the derivative to $0$. Furthermore,
define
$\zeta\colon(\sup\operatorname{supp}(Z),\infty)\times(0,\infty)\to{\mathbb{R}}$
as

$\displaystyle\zeta(\lambda;\Delta)$
$\displaystyle\coloneqq\psi(\max\\{\lambda,\overline{\lambda}(\Delta)\\};\Delta).$
(3.4)

Finally, for $i\in\\{1,2\\}$, by [MM19, Lemma 2], we have that the equation
$\zeta(\lambda;\delta_{i})=\varphi(\lambda)$ admits a unique solution in
$\lambda\in(\sup\operatorname{supp}(Z),\infty)$ which we call
$\lambda^{*}(\delta_{i})$. The functions
$\psi(\lambda;\Delta),\varphi(\lambda),\zeta(\lambda;\Delta)$ together with
the parameters $\lambda^{*}(\Delta),\overline{\lambda}(\Delta)$ are plotted in
Figure 2 for $\Delta\in\\{\delta,\delta_{1},\delta_{2}\\}$. Some convexity and
monotonicity properties (which will be useful later in the proofs) of these
functions can be found in Lemma I.2.

Figure 2: Plot of
$\psi(\lambda;\Delta),\varphi(\lambda),\zeta(\lambda;\Delta)$ as functions of
$\lambda$ with $\Delta\in\\{\delta,\delta_{1},\delta_{2}\\}$.

The empirical distribution of a vector $u\in{\mathbb{R}}^{d}$ is given by
$\frac{1}{d}\sum_{i=1}^{d}\delta_{u_{i}}$, where $\delta_{u_{i}}$ denotes a
Dirac delta mass on $u_{i}$. Similarly, the joint empirical distribution of
the rows of a matrix $(u^{1},u^{2},\ldots,u^{t})\in{\mathbb{R}}^{d\times t}$
is $\frac{1}{d}\sum_{i=1}^{d}\delta_{(u^{1}_{i},\ldots,u^{t}_{i})}$. Our
master theorem is an exact characterization in the high-dimensional limit of
the joint empirical distribution of the rows of the signals, the linear
estimator, and the spectral estimators. In particular, we show that this joint
empirical distribution converges to the law of a Gaussian random vector with a
specified covariance matrix. The result is stated in terms of the following
parameters: the asymptotic correlations
$\rho_{1}^{\mathrm{lin}},\rho_{2}^{\mathrm{lin}}$ between the linear estimator
and the two signals; the normalized length $n^{\mathrm{lin}}$ of the linear
estimator; and the asymptotic correlations
$\rho_{1}^{\mathrm{spec}},\rho_{2}^{\mathrm{spec}}$ between the spectral
estimators and the two signals. The formulas for these quantities are:

$\displaystyle
n^{\mathrm{lin}}\coloneqq\bigg{(}(\alpha^{2}+(1-\alpha)^{2})\mathbb{E}\left[G{\mathcal{L}}(Y)\right]^{2}+\frac{\mathbb{E}\left[{\mathcal{L}}(Y)^{2}\right]}{\delta}\bigg{)}^{\frac{1}{2}},$
(3.5a)
$\displaystyle\rho_{1}^{\mathrm{lin}}\coloneqq\frac{\alpha\mathbb{E}\left[G{\mathcal{L}}(Y)\right]}{n^{\mathrm{lin}}},\quad\rho_{2}^{\mathrm{lin}}\coloneqq\frac{(1-\alpha)\mathbb{E}\left[G{\mathcal{L}}(Y)\right]}{n^{\mathrm{lin}}},$
(3.5b)
$\displaystyle\rho_{1}^{\mathrm{spec}}\coloneqq\left(\frac{\frac{1}{\delta}-\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right)^{2}\right]}{\frac{1}{\delta}+\alpha\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right)^{2}(G^{2}-1)\right]}\right)^{\frac{1}{2}},\;\rho_{2}^{\mathrm{spec}}\coloneqq\left(\frac{\frac{1}{\delta}-\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{2})-Z}\right)^{2}\right]}{\frac{1}{\delta}+(1-\alpha)\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{2})-Z}\right)^{2}(G^{2}-1)\right]}\right)^{\frac{1}{2}}.$
(3.5c)

Theorem 3.1 is stated in terms of _pseudo-Lipschitz_ test functions. A
function $\Psi:{\mathbb{R}}^{m}\to{\mathbb{R}}$ is pseudo-Lipschitz of order
$k\geq 1$, denoted $\Psi\in{\rm PL}(k)$, if there is a constant $C>0$ such
that

$\left\|\Psi(x)-\Psi(y)\right\|_{2}\leq
C(1+\|x\|_{2}^{k-1}+\|y\|_{2}^{k-1})\left\|x-y\right\|_{2},$ (3.6)

for all $x,y\in\mathbb{R}^{m}$. Examples of pseudo-Lipschitz functions of
order two are: $\Psi(u)=u^{2}$ and $\Psi(u,v)=\left|uv\right|$, for
$u,v\in{\mathbb{R}}$. For simplicity, we consider pseudo-Lipschitz test
functions of order two, as those suffice to compute the asymptotic overlaps
between the signals and the various estimators. We note that one could extend
Theorem 3.1 to test functions in ${\rm PL}(k)$ for $k>2$, at the cost of a
more involved argument and an additional assumption on the boundedness of the
moments of $P_{\varepsilon}$.

######  Theorem 3.1 (Master theorem on joint distribution).

Consider the setting of Section 2, and let Items (A1), (A2), (A3), (A4), (A5),
(A6) and (A7) hold. Define the following rescaled vectors of length
$\sqrt{d}$:
$x^{\mathrm{lin}}=\sqrt{d}\,\widehat{x}^{\mathrm{lin}}/\left\|\widehat{x}^{\mathrm{lin}}\right\|_{2}$,
$\overline{x}^{*}_{i}=\sqrt{d}x_{i}^{*}$,
$x^{\mathrm{spec}}_{i}=s_{i}\sqrt{d}\,v_{i}(D)$, where the sign
$s_{i}\in\\{-1,1\\}$ is chosen such that $\left\langle
s_{i}v_{i}(D),x_{i}^{*}\right\rangle\geq 0$ ($i\in\\{1,2\\}$). Then, the
following holds almost surely for any ${\rm PL}(2)$ function
$\Psi:{\mathbb{R}}^{3}\to{\mathbb{R}}$. If
$\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$, then

$\displaystyle\lim_{d\to\infty}\,\frac{1}{d}\sum_{i=1}^{d}\Psi(\overline{x}^{*}_{1,i},x^{\mathrm{lin}}_{i},x^{\mathrm{spec}}_{1,i})=\mathbb{E}\left[\Psi(X_{1},\,\rho_{1}^{\mathrm{lin}}X_{1}+\rho_{2}^{\mathrm{lin}}X_{2}+W^{\mathrm{lin}},\,\rho_{1}^{\mathrm{spec}}X_{1}+W_{1}^{\mathrm{spec}})\right].$
(3.7)

Similarly, if $\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$, then

$\displaystyle\lim_{d\to\infty}\,\frac{1}{d}\sum_{i=1}^{d}\Psi(\overline{x}^{*}_{2,i},x^{\mathrm{lin}}_{i},x^{\mathrm{spec}}_{2,i})=\mathbb{E}\left[\Psi(X_{2},\,\rho_{1}^{\mathrm{lin}}X_{1}+\rho_{2}^{\mathrm{lin}}X_{2}+W^{\mathrm{lin}},\,\rho_{2}^{\mathrm{spec}}X_{2}+W_{2}^{\mathrm{spec}})\right].$
(3.8)

Here $(X_{1},X_{2})\sim{\mathcal{N}}(0,1)^{\otimes 2}$, the pairs
$(W^{\mathrm{lin}},W_{1}^{\mathrm{spec}})$ and
$(W^{\mathrm{lin}},W_{2}^{\mathrm{spec}})$ are independent of $(X_{1},X_{2})$
and each pair is jointly Gaussian with zero mean and covariance given by

$\begin{split}&{\mathbb{E}}[(W^{\mathrm{lin}})^{2}]=1-(\rho_{1}^{\mathrm{lin}})^{2}-(\rho_{2}^{\mathrm{lin}})^{2},\quad{\mathbb{E}}[(W_{1}^{\mathrm{spec}})^{2}]=1-(\rho_{1}^{\mathrm{spec}})^{2},\quad{\mathbb{E}}[(W_{2}^{\mathrm{spec}})^{2}]=1-(\rho_{2}^{\mathrm{spec}})^{2},\\\
&{\mathbb{E}}[W^{\mathrm{lin}}W_{1}^{\mathrm{spec}}]=\frac{\alpha\rho_{1}^{\mathrm{spec}}}{n^{\mathrm{lin}}}\mathbb{E}\left[\frac{G{\mathcal{L}}(Y)Z}{\lambda^{*}(\delta_{1})-Z}\right],\quad{\mathbb{E}}[W^{\mathrm{lin}}W_{2}^{\mathrm{spec}}]=\frac{(1-\alpha)\rho_{2}^{\mathrm{spec}}}{n^{\mathrm{lin}}}\mathbb{E}\left[\frac{G{\mathcal{L}}(Y)Z}{\lambda^{*}(\delta_{2})-Z}\right].\end{split}$

The proof of this result, given in Section 6, relies on the characterization
of the eigenvalues of the spectral matrix $D$ carried out in Theorem 5.1,
which is stated and proved in Section 5.

######  Remark 3.1 (Equivalence to convergence of empirical distribution).

The result in Equation 3.7 is equivalent to the statement that the joint
empirical distribution of
$(\overline{x}^{*}_{1},x^{\mathrm{lin}},x^{\mathrm{spec}}_{1})$ converges in
Wasserstein-2 distance to the joint law of
$(X_{2},\,\rho_{1}^{\mathrm{lin}}X_{1}+\rho_{2}^{\mathrm{lin}}X_{2}+W^{\mathrm{lin}},\,\rho_{1}^{\mathrm{spec}}X_{1}+W_{1}^{\mathrm{spec}})$.
A proof of the equivalence between convergence of empirical distributions in
Wasserstein distance and convergence of empirical averages of pseudo-Lipschitz
functions can be found in [FVRS22, Corollary 7.21].

######  Remark 3.2 (What if either the linear or spectral estimator is
ineffective).

The validity of the description of the joint law of the first signal and the
linear/spectral estimators in Equation 3.7 relies on two assumptions:
$\mathbb{E}\left[G{\mathcal{L}}(Y)\right]\neq 0$ for the linear estimator, and
$\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$ for the spectral one.
They guarantee that both estimators achieve nonzero asymptotic overlaps with
$x_{1}^{*}$, namely, $\rho_{1}^{\mathrm{lin}}\neq 0$ and
$\rho_{1}^{\mathrm{spec}}>0$. If one of ${\mathcal{L}}$ or ${\mathcal{T}}$
fails to satisfy the respective condition, then a conclusion similar to
Equation 3.7 still holds with
$\Psi\colon{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}$ only taking
$x_{1}^{*}$ and the effective estimator as inputs. This follows from a
straightforward adaptation of the proof of Theorem 3.1 and the formal
justification is omitted. An analogous discussion can be done as concerns the
second signal.

######  Remark 3.3 (Sign calibration of spectral estimator).

As the eigenvectors of a matrix are insensitive to sign flip, the spectral
estimators $x^{\mathrm{spec}}_{1},x^{\mathrm{spec}}_{2}$ are defined up to a
change of sign. In Theorem 3.1, we pick the signs so that the resulting
overlaps $\rho_{1}^{\mathrm{spec}},\rho_{2}^{\mathrm{spec}}$ are positive. In
practice, there is a simple way to resolve the sign ambiguity: one can match
the sign of
$\mathbb{E}\left[(\rho_{1}^{\mathrm{lin}}X_{1}+\rho_{2}^{\mathrm{lin}}X_{2}+W^{\mathrm{lin}})\,(\rho_{i}^{\mathrm{spec}}X_{i}+W_{i}^{\mathrm{spec}})\right]$
with that of the scalar product $\left\langle
x^{\mathrm{lin}},x^{\mathrm{spec}}_{i}\right\rangle$, as the latter can be
computed empirically (without knowing $x_{1}^{*},x_{2}^{*}$).

######  Remark 3.4 (Master theorem for $\alpha=1/2$).

Theorem 3.1 assumes $\alpha\in(1/2,1)$ (see Item (A2)). Nonetheless, the
conclusion of Theorem 3.1 continues to hold for $\alpha=1/2$ with a twist in
the definition of the spectral estimators. In this case, as $n\to\infty$ the
top two eigenvectors given by the spectral method correspond to the same
limiting eigenvalue. These eigenvectors, $v_{1}(D)$ and $v_{2}(D)$, provide an
estimate on the subspace spanned by $x_{1}^{*},x_{2}^{*}$. To estimate each
individual signal, we search for a vector in
$\operatorname{span}\left\\{v_{1}(D),v_{2}(D)\right\\}$ whose correlation with
$x^{\mathrm{lin}}$ is closest to the theoretical prediction from Theorem 3.1.
Indeed, let $x_{1}^{\mathrm{spec}},x_{2}^{\mathrm{spec}}$ be defined as

$\displaystyle x_{i}^{\mathrm{spec}}$
$\displaystyle\coloneqq\operatorname*{argmin}_{v\,\in\,\operatorname{span}\left\\{v_{1}(D),v_{2}(D)\right\\}\cap\sqrt{d}\,{\mathbb{S}}^{d-1}}\left|\frac{\left\langle
v,x^{\mathrm{lin}}\right\rangle}{\sqrt{d}}-\left(\rho_{i}^{\mathrm{lin}}\rho_{i}^{\mathrm{spec}}+\mathbb{E}\left[W^{\mathrm{lin}}W_{i}^{\mathrm{spec}}\right]\right)\right|,\quad\text{for
}i\in\\{1,2\\}.$ (3.9)

Then, Equations 3.7 and 3.8 hold, provided
$\mathbb{E}\left[G{\mathcal{L}}(Y)\right]\neq 0$ (which guarantees that the
linear estimator attains nonzero overlaps; see Item (A6) and Equation 3.5b).
We stress that Equation 3.9 is practically computable since it only involves
$x^{\mathrm{lin}}$ and theoretical predictions. If $x^{\mathrm{lin}}$ is
ineffective (which is the case, for example, in mixed phase retrieval; see
Section B.2), a similar grid search can still be performed if the statistician
is given as side information a vector with known correlation with a signal.
The reader is referred to Remarks 5.2 and 6.2 for the adaptation of our proofs
to the case $\alpha=1/2$.

Equipped with Theorem 3.1, we can combine the linear and spectral estimators
to improve the performance in the recovery of $x_{1}^{*}$ and $x_{2}^{*}$.
Formally, consider the (rescaled) linear and spectral estimators
$x^{\mathrm{lin}}\in\sqrt{d}\,{\mathbb{S}}^{d-1}$ and
$x_{1}^{\mathrm{spec}},x_{2}^{\mathrm{spec}}\in\sqrt{d}\,{\mathbb{S}}^{d-1}$.
Define

$\displaystyle
X^{\mathrm{lin}}\coloneqq\rho_{1}^{\mathrm{lin}}X_{1}+\rho_{2}^{\mathrm{lin}}X_{2}+W^{\mathrm{lin}},\quad
X_{1}^{\mathrm{spec}}\coloneqq\rho_{1}^{\mathrm{spec}}X_{1}+W_{1}^{\mathrm{spec}},\quad
X_{2}^{\mathrm{spec}}\coloneqq\rho_{2}^{\mathrm{spec}}X_{2}+W_{2}^{\mathrm{spec}}.$
(3.10)

Theorem 3.1 tells us the joint empirical distribution of the estimators
$(x^{\mathrm{lin}},x_{1}^{\mathrm{spec}},x_{2}^{\mathrm{spec}})$ converges to
the law of $(X^{\mathrm{lin}},X_{1}^{\mathrm{spec}},X_{2}^{\mathrm{spec}})$.
For $i\in\\{1,2\\}$, define the set of functions

$\displaystyle{\mathcal{C}}_{i}$
$\displaystyle\coloneqq\left\\{C_{i}\colon{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}\mbox{
s.t.
}\mathbb{E}\left[C_{i}(X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})^{2}\right]\in(0,\infty)\right\\}.$
(3.11)

Then, for any $C_{i}\in{\mathcal{C}}_{i}$, the _combined estimator_
$x_{i}^{\mathrm{comb}}$ is defined as

$\displaystyle x^{\mathrm{comb}}_{i}\coloneqq
C_{i}(x^{\mathrm{lin}},x_{i}^{\mathrm{spec}}),$ (3.12)

where $C_{i}$ acts on its inputs component-wise, i.e.,
$x_{i,j}^{\mathrm{comb}}=C_{i}(x_{j}^{\mathrm{lin}},x_{i,j}^{\mathrm{spec}})$
for any $j\in[d]$. Now, Equation 3.7 says that we can reduce the vector
problem of estimating $x_{i}^{*}$ given
$(x^{\mathrm{lin}},x^{\mathrm{spec}}_{i})$ to an estimation problem over
scalar random variables, i.e., how to optimally estimate $X_{i}$ from
observations $X^{\mathrm{lin}}$ and $X_{i}^{\mathrm{spec}}$. The Bayes-optimal
combined estimator that minimizes the expected squared error for this scalar
problem is
$\mathbb{E}\left[X_{i}\,|\,X^{\mathrm{lin}},X_{i}^{\mathrm{spec}}\right]$.
Recalling from Theorem 3.1 that
$(X_{i},X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})$ are jointly Gaussian, the
Bayes-optimal combined estimator is a linear combination of
$(X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})$. The performance of this combined
estimator is formalized in the following corollary, whose proof is contained
in Appendix A.

######  Corollary 3.2 (Bayes-optimal linear-spectral combination).

Consider the setting of Theorem 3.1. For $i\in\\{1,2\\}$, define
$C_{i}^{*}\colon{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}$ as follows:

$\displaystyle C^{*}_{i}(X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})$
$\displaystyle\coloneqq\mathbb{E}\left[X_{i}\,|\,X^{\mathrm{lin}},\,X_{i}^{\mathrm{spec}}\right]=\frac{1}{1-\nu_{i}^{2}}\,\left(\xi_{i}X^{\mathrm{lin}}\,+\,\zeta_{i}X_{i}^{\mathrm{spec}}\right),$
(3.13)

where

$\displaystyle\nu_{i}\coloneqq\rho_{i}^{\mathrm{lin}}\rho_{i}^{\mathrm{spec}}+\mathbb{E}\left[W^{\mathrm{lin}}W_{i}^{\mathrm{spec}}\right],\quad\xi_{i}\coloneqq\rho_{i}^{\mathrm{lin}}-\rho_{i}^{\mathrm{spec}}\nu_{i},\quad\zeta_{i}\coloneqq\rho_{i}^{\mathrm{spec}}-\rho_{i}^{\mathrm{lin}}\nu_{i}.$

For $i\in\\{1,2\\}$, let $x_{i}^{\mathrm{comb}}$ be the combined estimators
defined in Equation 3.12 w.r.t. $C_{i}^{*}$, respectively. Then, almost surely
we have

$\displaystyle\lim_{d\to\infty}\frac{\left\langle
x_{i}^{\mathrm{comb}},x_{i}^{*}\right\rangle}{\left\|x_{i}^{\mathrm{comb}}\right\|_{2}\left\|x_{i}^{*}\right\|_{2}}$
$\displaystyle=\frac{1}{1-\nu_{i}^{2}}\left(\xi_{i}^{2}+\zeta_{i}^{2}+2\xi_{i}\zeta_{i}\left(\rho_{i}^{\mathrm{lin}}\rho_{i}^{\mathrm{spec}}+\mathbb{E}\left[W^{\mathrm{lin}}W_{i}^{\mathrm{spec}}\right]\right)\right)^{1/2}\eqqcolon\mathsf{OL}_{i}^{\mathrm{comb}}.$

Furthermore, for any
$(C_{1},C_{2})\in{\mathcal{C}}_{1}\times{\mathcal{C}}_{2}$, the corresponding
combined estimators
$\widetilde{x}_{1}^{\mathrm{comb}},\widetilde{x}_{2}^{\mathrm{comb}}$ defined
w.r.t. $C_{1},C_{2}$ through Equation 3.12, satisfy

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle\widetilde{x}_{i}^{\mathrm{comb}},x_{i}^{*}\right\rangle\right|}{\left\|\widetilde{x}_{i}^{\mathrm{comb}}\right\|_{2}\left\|x_{i}^{*}\right\|_{2}}$
$\displaystyle=\frac{\left|\mathbb{E}\left[X_{i}C_{i}(X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})\right]\right|}{\sqrt{\mathbb{E}\left[C_{i}(X^{\mathrm{lin}},X_{i}^{\mathrm{spec}})^{2}\right]}}\leq\mathsf{OL}_{i}^{\mathrm{comb}},$

for $i\in\\{1,2\\}$.

###  3.1 Linear estimator

Theorem 3.1 allows us to derive the asymptotic overlap of each signal with a
linear estimator defined via a given preprocessing function
${\mathcal{L}}\colon{\mathbb{R}}\to{\mathbb{R}}$.

######  Corollary 3.3 (Overlaps, linear).

Consider the setting of Section 2, and let Items (A1), (A2), (A3), (A4), (A5)
and (A6) hold. Then, almost surely,

$\displaystyle\lim_{d\to\infty}\,\frac{\left\langle\widehat{x}^{\mathrm{lin}},x_{i}^{*}\right\rangle}{\left\|\widehat{x}^{\mathrm{lin}}\right\|_{2}\left\|x_{i}^{*}\right\|_{2}}$
$\displaystyle=\rho_{i}^{\mathrm{lin}}\,,\quad i\in\\{1,2\\}.$ (3.14)

###### Proof.

Choose $\Psi(a,b,c)=ab$, and note that $\Psi\in{\rm PL}(2)$. Then, recalling
that
$\left\|\widehat{x}^{\mathrm{lin}}\right\|_{2}=\left\|\overline{x}^{*}_{i}\right\|_{2}=\sqrt{d}$,
the left side of Equations 3.7 and 3.8 recover the overlaps in Equation 3.14
for $i=1,2$, respectively. The right sides of Equations 3.7 and 3.8 become
$\rho_{1}^{\mathrm{lin}},\rho_{2}^{\mathrm{lin}}$ (defined in Equation 3.5b).∎

######  Remark 3.5 (Overlap of linear estimator does not approach $1$).

From Equation 3.14 and the definitions of
$\rho_{1}^{\mathrm{lin}},\rho_{2}^{\mathrm{lin}}$ in Equation 3.5b, we observe
that the linear estimator achieves positive overlap with each signal for _any_
positive $\delta$, as long as $\mathbb{E}\left[G{\mathcal{L}}(Y)\right]>0$. As
$\delta\to\infty$, the limiting overlaps approach

$\displaystyle\sqrt{\frac{\alpha^{2}}{\alpha^{2}+(1-\alpha)^{2}}},\quad\sqrt{\frac{(1-\alpha)^{2}}{\alpha^{2}+(1-\alpha)^{2}}},$

respectively. These quantities are strictly less than $1$ for any
$\alpha\in(1/2,1)$. In contrast, the overlap of the spectral estimator becomes
positive only when $\delta$ exceeds a certain threshold (see Remark 3.8).
However, once this threshold is exceeded, the spectral estimator yields
overlaps approaching $1$ as $\delta$ grows (see Remark 3.9).

With the above characterization of the limiting overlap of a linear estimator,
we can optimize the performance over the choice of preprocessing function
${\mathcal{L}}\colon{\mathbb{R}}\to{\mathbb{R}}$ (subject to Item (A6)). Let

$\displaystyle{\mathcal{I}}$
$\displaystyle\coloneqq\left\\{{\mathcal{L}}\colon{\mathbb{R}}\to{\mathbb{R}}\,\text{
Lipschitz}\mbox{ s.t. }\mathbb{E}\left[G{\mathcal{L}}(Y)\right]\neq
0,\;\mathbb{E}\left[\left|G{\mathcal{L}}(Y)\right|\right]<\infty\right\\}$

be the set of functions ${\mathcal{L}}$ satisfying Item (A6). For
$i\in\\{1,2\\}$ and $\delta\in(0,\infty)$, define the _optimal overlaps_ among
linear estimators as

$\displaystyle\mathsf{OL}_{i}^{\mathrm{lin}}$
$\displaystyle\coloneqq\sup_{{\mathcal{L}}\in{\mathcal{I}}}\rho_{i}^{\mathrm{lin}}.$

Furthermore, if ${\mathcal{I}}=\emptyset$, we simply set
$\mathsf{OL}_{1}^{\mathrm{lin}}=\mathsf{OL}_{2}^{\mathrm{lin}}=0$. Recall that
$\rho_{1}^{\mathrm{lin}},\rho_{2}^{\mathrm{lin}}$ defined in Equation 3.5b
depend on $\delta$ and ${\mathcal{L}}$. Therefore,
$\mathsf{OL}_{1}^{\mathrm{lin}},\mathsf{OL}_{2}^{\mathrm{lin}}$ also depend on
$\delta$. In words, $\mathsf{OL}_{i}^{\mathrm{lin}}$ ($i\in\\{1,2\\}$) is the
largest overlap with the $i$-th signal that can be achieved by a linear
estimator. Then, we have the following characterization of the optimal
overlaps. The proof is contained in Section B.1.

######  Proposition 3.4 (Optimal linear estimator).

Consider the setting of Section 2, and let Items (A1), (A2), (A3), (A4) and
(A5) hold. Assume further that

$\displaystyle\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[Gp(y|G)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle\in(0,\infty),$ (3.15)

where $p(y|g)$ is the conditional law in Equation 2.2 and the expectation is
taken w.r.t. $G\sim{\mathcal{N}}(0,1)$. Then, for any $\delta\in(0,\infty)$,
we have

$\displaystyle\mathsf{OL}_{1}^{\mathrm{lin}}$
$\displaystyle=\left(\frac{\alpha^{2}+(1-\alpha)^{2}}{\alpha^{2}}+\frac{1}{\alpha^{2}\delta}\cdot\frac{1}{\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[Gp(y|G)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y}\right)^{-1/2},$
$\displaystyle\mathsf{OL}_{2}^{\mathrm{lin}}$
$\displaystyle=\left(\frac{\alpha^{2}+(1-\alpha)^{2}}{(1-\alpha)^{2}}+\frac{1}{(1-\alpha)^{2}\delta}\cdot\frac{1}{\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[Gp(y|G)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y}\right)^{-1/2}.$
(3.16)

Moreover, define ${\mathcal{L}}^{*}\colon{\mathbb{R}}\to{\mathbb{R}}$ as

$\displaystyle{\mathcal{L}}^{*}(y)$
$\displaystyle=\frac{\mathbb{E}\left[Gp(y|G)\right]}{\mathbb{E}\left[p(y|G)\right]}.$

Then, ${\mathcal{L}}^{*}\in{\mathcal{I}}$ and for any $\delta\in(0,\infty)$,
both $\mathsf{OL}_{1}^{\mathrm{lin}},\mathsf{OL}_{2}^{\mathrm{lin}}$ are
simultaneously achieved by ${\mathcal{L}}^{*}$.

######  Remark 3.6 (When linear estimator is ineffective).

The condition in Equation 3.15 ensures that the linear estimator
asymptotically achieves strictly positive overlap with the signals. In fact,
if

$\displaystyle\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[Gp(y|G)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=0,$

then, by inspecting the RHS of Equation 3.16, we readily obtain that
$\mathsf{OL}_{1}^{\mathrm{lin}}=\mathsf{OL}_{2}^{\mathrm{lin}}=0$ for any
$\delta\in(0,\infty)$. For example, this is the case for mixed phase retrieval
(see Section B.2). We remark that the condition in Equation 3.15 also appears
in the non-mixed setting (see Appendix C.1 of [MTV21]).

###  3.2 Spectral estimator

Theorem 3.1 allows us to identify a spectral threshold for an arbitrary
preprocessing function (subject to Item (A7)). Specifically, for any fixed
${\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$, it provides an explicit
sufficient condition for the spectral estimator defined via ${\mathcal{T}}$ to
have strictly positive overlap with the signals. The limiting value of the
overlaps can also be obtained similarly to Corollary 3.7.

######  Corollary 3.5 (Overlaps, spectral).

Consider the setting of Section 2, and let Items (A1), (A2), (A3), (A4), (A5)
and (A7) hold. Then, for $i\in\\{1,2\\}$, if
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$, we have that, almost
surely,

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{i}(D),x_{i}^{*}\right\rangle\right|}{\left\|v_{i}(D)\right\|_{2}\left\|x_{i}^{*}\right\|_{2}}$
$\displaystyle=\rho_{i}^{\mathrm{spec}}.$ (3.17)

######  Remark 3.7 (Condition for vanishing overlap).

We focus here on the recovery of the first signal, and an analogous discussion
is valid for the second one. Let us note that, as $\overline{\lambda}(\delta)$
approaches $\lambda^{*}(\delta_{1})$, the RHS of Equation 3.17 is $0$. Indeed,
as $\lambda^{*}(\delta_{1})\searrow\overline{\lambda}(\delta)$, we have
$\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right)^{2}\right]\nearrow\frac{1}{\delta}$
using Item 4 of Lemma I.4 and Equation I.2, and consequently the numerator of
$\rho_{1}^{\mathrm{spec}}$ (cf. Equation 3.5c) decreases to $0$. Furthermore,
in the non-mixed setting ($\alpha=1$), the analysis of [LL20, MM19] gives
that, when $\lambda^{*}(\delta_{1})<\overline{\lambda}(\delta)$, the
corresponding overlap vanishes. While we do not formally prove that the
condition $\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$ is _necessary_
for the spectral method to have non-vanishing overlap, these two observations
point strongly in that direction. A third piece of supporting evidence is
provided in Remark 5.1.

Equipped with Corollary 3.5, we can optimize both _(i)_ the spectral
threshold, namely, the minimum value of $\delta$ needed to satisfy the
condition $\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$ which gives a
strictly positive overlap, and _(ii)_ the limiting overlap given by the right
side of Equation 3.17. Formally, for $i\in\\{1,2\\}$ and
$\delta\in(0,\infty)$, let

$\displaystyle{\mathcal{H}}_{i}$
$\displaystyle\coloneqq\left\\{{\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}\text{
Lipschitz}\mbox{ s.t.
}\begin{array}[]{c}\displaystyle\inf_{y\in\operatorname{supp}(Y)}{\mathcal{T}}(y)>-\infty,\;\displaystyle
0<\sup_{y\in\operatorname{supp}(Y)}{\mathcal{T}}(y)<\infty,\\\
\Pr\left[{\mathcal{T}}(Y)=0\right]<1,\;\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)\end{array}\right\\}$
(3.20)

be the set of functions ${\mathcal{T}}$ satisfying Item (A7) such that
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$ holds. Here, we recall
that $\delta_{1}=\alpha\delta,\delta_{2}=(1-\alpha)\delta$ and
$\lambda^{*}(\cdot),\overline{\lambda}(\cdot)$ depend on the choice of the
preprocessing function. We stress that ${\mathcal{H}}_{i}$ depends on
$\delta$. Now, we can define the _spectral threshold_ for the $i$-th signal as

$\displaystyle\delta_{i}^{\mathrm{spec}}$
$\displaystyle\coloneqq\inf\left\\{\delta\in(0,\infty):{\mathcal{H}}_{i}\neq\emptyset\right\\}\,,\quad
i\in\\{1,2\\}.$

In words, this is the smallest $\delta$ such that there exists a preprocessing
function satisfying $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$ (and,
hence, leading to non-vanishing limiting overlap). Furthermore, for
$i\in\\{1,2\\}$ and $\delta>\delta_{i}^{\mathrm{spec}}$, define the _optimal
overlap_ as

$\displaystyle\mathsf{OL}_{i}^{\mathrm{spec}}$
$\displaystyle\coloneqq\sup_{{\mathcal{T}}\in{\mathcal{H}}_{i}}\rho_{i}^{\mathrm{spec}}.$

Recall that $\rho_{1}^{\mathrm{spec}},\rho_{2}^{\mathrm{spec}}$ defined in
Equation 3.5c depend on $\delta$ and $Z$ which in turn depends on
${\mathcal{T}}$. In words, $\mathsf{OL}_{i}^{\mathrm{spec}}$ is the largest
overlap with preprocessing functions that satisfy
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$. We note that the
supremum is guaranteed to be taken over a nonempty set as
$\delta>\delta_{i}^{\mathrm{spec}}$, and naturally
$\mathsf{OL}_{i}^{\mathrm{spec}}$ also depends on $\delta$. At this point, we
can state the following result whose proof is given in Appendix C.

######  Proposition 3.6 (Optimal spectral estimator).

Consider the setting of Section 2, and let Items (A1), (A2), (A3), (A4) and
(A5) hold. Then, we have

$\displaystyle\delta_{1}^{\mathrm{spec}}$
$\displaystyle=\frac{1}{\alpha^{2}\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y},\quad\delta_{2}^{\mathrm{spec}}=\frac{1}{(1-\alpha)^{2}\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y},$
(3.21)

and

$\displaystyle\mathsf{OL}_{1}^{\mathrm{spec}}$
$\displaystyle=\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha)+\alpha}},\quad\mathsf{OL}_{2}^{\mathrm{spec}}=\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha)+(1-\alpha)}},$
(3.22)

where $\beta_{1}^{*}(\delta,\alpha)\in(1-\alpha,\infty)$ and
$\beta_{2}^{*}(\delta,\alpha)\in(\alpha,\infty)$ are the unique solutions to
the following fixed point equations:

$\displaystyle(\beta_{1}^{*}(\delta,\alpha)-(1-\alpha))\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+\beta_{1}^{*}(\delta,\alpha)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\frac{1}{\alpha^{2}\delta},$ (3.23)
$\displaystyle(\beta_{2}^{*}(\delta,\alpha)-\alpha)\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{(1-\alpha)\mathbb{E}\left[p(y|G)G^{2}\right]+\beta_{2}^{*}(\delta,\alpha)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\frac{1}{(1-\alpha)^{2}\delta},$ (3.24)

respectively. Finally, define
${\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}\colon{\mathbb{R}}\to{\mathbb{R}}$
as follows

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$
$\displaystyle=1-\frac{1}{\alpha\cdot\frac{\mathbb{E}\left[p(y|G)G^{2}\right]}{\mathbb{E}\left[p(y|G)\right]}+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)\cdot\frac{\mathbb{E}\left[p(y|G)G^{2}\right]}{\mathbb{E}\left[p(y|G)\right]}+\alpha},$
(3.25)

where $G\sim{\mathcal{N}}(0,1)$. Then, for $i\in\\{1,2\\}$ and
$\delta>\delta_{i}^{\mathrm{spec}}$, we have that _(i)_
${\mathcal{T}}_{i}^{*}\in{\mathcal{H}}_{i}$, and _(ii)_ the value of
$\mathsf{OL}_{i}^{\mathrm{spec}}$ is achieved by ${\mathcal{T}}_{i}^{*}$.

######  Remark 3.8 (Universal lower bounds on spectral thresholds).

In Appendix D, we show that the spectral thresholds
$\delta_{1}^{\mathrm{spec}}$ and $\delta_{2}^{\mathrm{spec}}$ are always at
least $\delta_{1}^{*}\coloneqq\frac{1}{2\alpha^{2}}$ and
$\delta_{2}^{*}\coloneqq\frac{1}{2(1-\alpha)^{2}}$, respectively, regardless
of the choice of the conditional law $p(\cdot\,|\,g)$ in Equation 2.2 (i.e.,
regardless of the model). These lower bounds are met by both noiseless linear
regression and noiseless phase retrieval, see Remark 3.10. The bounds imply,
in particular, that unlike the linear estimator (cf. Section 3.1), the
spectral estimator (even the optimal one) does not achieve weak recovery for
all $\delta>0$; it starts producing positive overlaps only when the aspect
ratio $\delta$ exceeds a certain (strictly positive) spectral threshold.
Furthermore, in stark contrast with the non-mixed setting, having access to
the sign of the observations does not help spectral methods for noiseless
phase retrieval, since the optimal preprocessing function for mixed linear
regression effectively transforms the problem into mixed phase retrieval.

######  Remark 3.9 (Overlap of spectral estimator approaches $1$).

The optimal limiting overlaps in Equation 3.22 approach $1$ as
$\delta\to\infty$ provided

$\displaystyle\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+(1-\alpha)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle\,\in[0,\infty).$ (3.26)

To show this, consider the optimal limiting overlap between the spectral
estimator and the first signal, which by Equation 3.22 equals
$\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha)+\alpha}}$. To show the claim, it
suffices to show
$\beta_{1}^{*}(\delta,\alpha)\xrightarrow{\delta\to\infty}1-\alpha$. From
Equation 3.23, the fixed point equation defining
$\beta_{1}^{*}(\infty,\alpha)$ becomes

$\displaystyle(\beta_{1}^{*}(\infty,\alpha)-(1-\alpha))\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+\beta_{1}^{*}(\infty,\alpha)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=0,$ (3.27)

as $\delta\to\infty$. Since Equation 3.26 holds, the unique solution to
Equation 3.27 has to be $\beta_{1}^{*}(\infty,\alpha)=1-\alpha$. This proves
the claim. We note that the condition Equation 3.26 is satisfied by the mixed
linear regression model (see Section E.1).

###  3.3 Two illustrative examples

We specialize the results in Sections 3.2 and 3.1 regarding linear and
spectral estimators to two prototypical examples of mixed GLMs: the _mixed
linear regression_ model with link function given by

$\displaystyle q(g,\varepsilon)$
$\displaystyle=g+\varepsilon,\quad{\underline{\varepsilon}}\sim{\mathcal{N}}(0,\sigma^{2}I_{n}),$
(3.28)

and the _mixed phase retrieval_ model where

$\displaystyle q(g,\varepsilon)$
$\displaystyle=|g|+\varepsilon,\quad{\underline{\varepsilon}}\sim{\mathcal{N}}(0,\sigma^{2}I_{n}).$
(3.29)

The explicit formulas for the optimal preprocessing functions, the
corresponding optimal overlaps and the thresholds of aspect ratios (for
spectral estimators) are collected in the following corollaries.

Let us first consider linear estimators.

######  Corollary 3.7 (Mixed linear regression, linear estimator).

Consider the mixed linear regression model in Equation 3.28, and let Items
(A1), (A2), (A4) and (A5) hold. Then, the optimal preprocessing function
${\mathcal{L}}^{*}$ defined in Proposition 3.4 is given by

$\displaystyle{\mathcal{L}}^{*}(y)$ $\displaystyle=\frac{y}{1+\sigma^{2}}.$
(3.30)

With
$\widehat{x}^{\mathrm{lin}}\coloneqq\frac{1}{n}A^{\top}{\mathcal{L}}^{*}(y)$,
we have that, almost surely,

$\displaystyle\lim_{d\to\infty}\,\frac{\left\langle\widehat{x}^{\mathrm{lin}},x_{1}^{*}\right\rangle}{\left\|\widehat{x}^{\mathrm{lin}}\right\|_{2}\left\|x_{1}^{*}\right\|_{2}}$
$\displaystyle=\left(\frac{\alpha^{2}+(1-\alpha)^{2}}{\alpha^{2}}+\frac{1+\sigma^{2}}{\alpha^{2}\delta}\right)^{-1},$
$\displaystyle\lim_{d\to\infty}\,\frac{\left\langle\widehat{x}^{\mathrm{lin}},x_{2}^{*}\right\rangle}{\left\|\widehat{x}^{\mathrm{lin}}\right\|_{2}\left\|x_{2}^{*}\right\|_{2}}$
$\displaystyle=\left(\frac{\alpha^{2}+(1-\alpha)^{2}}{(1-\alpha)^{2}}+\frac{1+\sigma^{2}}{(1-\alpha)^{2}\delta}\right)^{-1}.$

The proof of the above corollary is in Section B.2. In contrast, for mixed
phase retrieval, it is easy to check (see again Section B.2) that if the
linear estimator is applied, the overlaps with both signals are always
vanishing regardless of the choice of the preprocessing function.

We then turn to spectral estimators. The proofs of the following two
corollaries (Corollaries 3.8 and 3.9) can be found in Sections E.1 and E.2,
respectively.

######  Corollary 3.8 (Mixed linear regression, spectral estimator).

Consider the mixed linear regression model and let Items (A1), (A2), (A4) and
(A5) hold. Then, the optimal preprocessing functions
${\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}$ defined in Proposition 3.6 are
given by:

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$
$\displaystyle=1-\frac{1}{\alpha\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\alpha}.$
(3.31)

Let $T_{i}^{*}=\mathrm{diag}({\mathcal{T}}_{i}^{*}(y))$ and
$D_{i}^{*}=\frac{1}{n}A^{\top}T_{i}^{*}A$ for $i\in\\{1,2\\}$. Denote by
$v_{1}(D_{i}^{*}),v_{2}(D_{i}^{*})$ the eigenvectors of $D_{i}^{*}$
corresponding to the two largest eigenvalues, respectively. Then for
$\delta>\frac{(1+\sigma^{2})^{2}}{2\alpha^{2}}$, we have

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{1}(D_{1}^{*}),x_{1}^{*}\right\rangle\right|}{\left\|v_{1}(D_{1}^{*})\right\|_{2}\left\|x_{1}^{*}\right\|_{2}}$
$\displaystyle=\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha,\sigma)+\alpha}},$

almost surely, where $\beta_{1}^{*}(\delta,\alpha,\sigma)$ is the unique
solution in $(1-\alpha,\infty)$ to the following fixed point equation:

$(\beta_{1}^{*}(\delta,\alpha,\sigma)-(1-\alpha))\left[-\frac{\alpha+\beta_{1}^{*}(\delta,\alpha,\sigma)}{\alpha^{2}}+\left(\frac{\alpha+\beta_{1}^{*}(\delta,\alpha,\sigma)}{\alpha}\right)^{2}\sqrt{\frac{\pi(1+\sigma^{2})^{2}}{2\alpha(\sigma^{2}\alpha+(1+\sigma^{2})\beta_{1}^{*}(\delta,\alpha,\sigma))}}\right.\\\
\left.\times\exp\left(\frac{\sigma^{2}\alpha+(1+\sigma^{2})\beta_{1}^{*}(\delta,\alpha,\sigma)}{2\alpha}\right)\mathrm{erfc}\left(\sqrt{\frac{\sigma^{2}\alpha+(1+\sigma^{2})\beta_{1}^{*}(\delta,\alpha,\sigma)}{2\alpha}}\right)\right]=\frac{1}{\alpha^{2}\delta}.$
(3.32)

For $\delta>\frac{(1+\sigma^{2})^{2}}{2(1-\alpha)^{2}}$, we have

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{2}(D_{2}^{*}),x_{2}^{*}\right\rangle\right|}{\left\|v_{2}(D_{2}^{*})\right\|_{2}\left\|x_{2}^{*}\right\|_{2}}$
$\displaystyle=\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha,\sigma)+(1-\alpha)}},$

almost surely, where $\beta_{2}^{*}(\delta,\alpha,\sigma)$ is the unique
solution in $(\alpha,\infty)$ to the following fixed point equation:

$(\beta_{2}^{*}(\delta,\alpha,\sigma)-\alpha)\left[-\frac{1-\alpha+\beta_{2}^{*}(\delta,\alpha,\sigma)}{(1-\alpha)^{2}}+\left(\frac{1-\alpha+\beta_{2}^{*}(\delta,\alpha,\sigma)}{1-\alpha}\right)^{2}\right.\\\
\times\sqrt{\frac{\pi(1+\sigma^{2})^{2}}{2(1-\alpha)(\sigma^{2}(1-\alpha)+(1+\sigma^{2})\beta_{2}^{*}(\delta,\alpha,\sigma))}}\exp\left(\frac{\sigma^{2}(1-\alpha)+(1+\sigma^{2})\beta_{2}^{*}(\delta,\alpha,\sigma)}{2(1-\alpha)}\right)\\\
\left.\times\mathrm{erfc}\left(\sqrt{\frac{\sigma^{2}(1-\alpha)+(1+\sigma^{2})\beta_{2}^{*}(\delta,\alpha,\sigma)}{2(1-\alpha)}}\right)\right]=\frac{1}{(1-\alpha)^{2}\delta}.$
(3.33)

######  Corollary 3.9 (Mixed phase retrieval, spectral).

Consider the mixed phase retrieval model in Equation 3.29, and let Items (A1),
(A2), (A4) and (A5) hold. Then, the optimal preprocessing functions
${\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}$ defined in Proposition 3.6 are
given by:

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$
$\displaystyle=1-\frac{1}{\alpha\Delta(y)+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)\Delta(y)+\alpha},$
(3.34)

where the auxiliary function $\Delta\colon{\mathbb{R}}\to{\mathbb{R}}$ is
defined as

$\displaystyle\Delta(y)$
$\displaystyle\coloneqq\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\sqrt{\frac{2}{\pi}}\cdot\frac{\sigma
y\exp\left(-\frac{y^{2}}{2\sigma^{2}(1+\sigma^{2})}\right)}{(1+\sigma^{2})^{3/2}\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]}.$

Let $T_{i}^{*}=\mathrm{diag}({\mathcal{T}}_{i}^{*}(y))\in{\mathbb{R}}^{n\times
n}$ and $D_{i}^{*}=\frac{1}{n}A^{\top}T_{i}^{*}A\in{\mathbb{R}}^{d\times d}$
for $i\in\\{1,2\\}$. Denote by $v_{1}(D_{i}^{*}),v_{2}(D_{i}^{*})$ the
eigenvectors of $D_{i}^{*}$ corresponding to the two largest eigenvalues,
respectively. Let

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{1}{\alpha^{2}}\left(\frac{2}{(1+\sigma^{2})^{2}}+\frac{4\sigma^{5}h(\sigma^{2})}{\pi^{3/2}(1+\sigma^{2})^{2}}\right)^{-1},\quad\delta^{*}_{2}=\frac{1}{(1-\alpha)^{2}}\left(\frac{2}{(1+\sigma^{2})^{2}}+\frac{4\sigma^{5}h(\sigma^{2})}{\pi^{3/2}(1+\sigma^{2})^{2}}\right)^{-1},$

where
$h(\sigma^{2})\coloneqq\int_{\mathbb{R}}\frac{\exp\left(-(2+\sigma^{2})z^{2}\right)z^{2}}{1+\mathrm{erf}(z)}\mathrm{d}z$.
Define the functions $m_{0},m_{1}\colon{\mathbb{R}}\to{\mathbb{R}}$ and
$I\colon[1/2,1]\times(0,\infty)\to{\mathbb{R}}$ as

$\displaystyle m_{0}(y)$
$\displaystyle\coloneqq\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right],$
$\displaystyle m_{1}(y)$ $\displaystyle\coloneqq
m_{0}(y)\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right),$
$\displaystyle I(\alpha,\beta)$
$\displaystyle\coloneqq\int_{\operatorname{supp}(Y)}\frac{(m_{2}(y)-m_{0}(y))^{2}}{\alpha
m_{2}(y)+\beta m_{0}(y)}\mathrm{d}y.$

Then for $\delta>\delta^{*}_{1}$, we have

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{1}(D_{1}^{*}),x_{1}^{*}\right\rangle\right|}{\left\|v_{1}(D_{1}^{*})\right\|_{2}\left\|x_{1}^{*}\right\|_{2}}$
$\displaystyle=\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha,\sigma)+\alpha}},$

almost surely, where $\beta_{1}^{*}(\delta,\alpha,\sigma)$ is the unique
solution in $(1-\alpha,\infty)$ to the following fixed point equation:

$\displaystyle(\beta_{1}^{*}(\delta,\alpha,\sigma)-(1-\alpha))I(\alpha,\beta_{1}^{*}(\delta,\alpha,\sigma))$
$\displaystyle=\frac{1}{\alpha^{2}\delta}.$

For $\delta>\delta^{*}_{2}$, we have

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{2}(D_{2}^{*}),x_{2}^{*}\right\rangle\right|}{\left\|v_{2}(D_{2}^{*})\right\|_{2}\left\|x_{2}^{*}\right\|_{2}}$
$\displaystyle=\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha,\sigma)+(1-\alpha)}},$

almost surely, where $\beta_{2}^{*}(\delta,\alpha,\sigma)$ is the unique
solution in $(\alpha,\infty)$ to the following fixed point equation:

$\displaystyle(\beta_{2}^{*}(\delta,\alpha,\sigma)-\alpha)I(1-\alpha,\beta_{2}^{*}(\delta,\alpha,\sigma))$
$\displaystyle=\frac{1}{(1-\alpha)^{2}\delta}.$

######  Remark 3.10 (Linear regression vs. phase retrieval).

We note that the performance of the optimal spectral estimators given in
Corollaries 3.8 and 3.9 coincides for mixed _noiseless_ linear regression and
mixed _noiseless_ phase retrieval. Specifically, for both models, when
$\sigma^{2}=0$, the spectral thresholds are:

$\displaystyle\delta_{1}^{*}$
$\displaystyle=\frac{1}{2\alpha^{2}},\quad\delta_{2}^{*}=\frac{1}{2(1-\alpha)^{2}},$

the optimal preprocessing functions are:

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$ $\displaystyle=1-\frac{1}{\alpha
y^{2}+1-\alpha},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)y^{2}+\alpha},$
(3.35)

and the corresponding overlaps are:

$\displaystyle\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha)+\alpha}},\quad\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha)+1-\alpha}},$

where
$\beta_{1}^{*}(\delta,\alpha)\in(1-\alpha,\infty),\,\beta_{2}^{*}(\delta,\alpha)\in(\alpha,\infty)$
are the unique solutions to the following fixed point equations:

$(\beta_{1}^{*}(\delta,\alpha)-(1-\alpha))\left[-\frac{\alpha+\beta_{1}^{*}(\delta,\alpha)}{\alpha^{2}}+\sqrt{\frac{\pi}{2\alpha\beta_{1}^{*}(\delta,\alpha)}}\left(\frac{\alpha+\beta_{1}^{*}(\delta,\alpha)}{\alpha}\right)^{2}\right.\\\
\left.\exp\left(\frac{\beta_{1}^{*}(\delta,\alpha)}{2\alpha}\right)\mathrm{erfc}\left(\sqrt{\frac{\beta_{1}^{*}(\delta,\alpha)}{2\alpha}}\right)\right]=\frac{1}{\alpha^{2}\delta},$
(3.36)

and

$(\beta_{2}^{*}(\delta,\alpha)-\alpha)\left[-\frac{(1-\alpha)+\beta_{2}^{*}(\delta,\alpha)}{(1-\alpha)^{2}}+\sqrt{\frac{\pi}{2(1-\alpha)\beta_{2}^{*}(\delta,\alpha)}}\left(\frac{(1-\alpha)+\beta_{2}^{*}(\delta,\alpha)}{1-\alpha}\right)^{2}\right.\\\
\left.\exp\left(\frac{\beta_{2}^{*}(\delta,\alpha)}{2(1-\alpha)}\right)\mathrm{erfc}\left(\sqrt{\frac{\beta_{2}^{*}(\delta,\alpha)}{2(1-\alpha)}}\right)\right]=\frac{1}{(1-\alpha)^{2}\delta},$
(3.37)

respectively. In fact, even the first-order dependence of the spectral
threshold on the noise variance $\sigma$ coincides for the noisy versions of
these two problems (cf. Equations E.8 and E.12):

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{1+2\sigma^{2}}{2\alpha^{2}}+{\mathcal{O}}(\sigma^{4}),\quad\delta^{*}_{2}=\frac{1+2\sigma^{2}}{2(1-\alpha)^{2}}+{\mathcal{O}}(\sigma^{4}).$

This phenomenon is due to the fact that the optimal preprocessing functions
${\mathcal{T}}_{1}^{*}(y),{\mathcal{T}}_{2}^{*}(y)$ in both cases depend only
on $y^{2}$, and are therefore invariant to the signs of the channel outputs
$(y_{i})_{i=1}^{n}$.

###  3.4 Technical tools

The proof of the master theorem (Theorem 3.1) uses a combination of two tools:
approximate message passing (AMP) and random matrix theory (RMT). We now
outline the high-level ideas in the analysis.

##### Joint distribution via GAMP.

The convergence results in Equations 3.7 and 3.8 are obtained using AMP of a
particular form known as generalized approximate message passing (GAMP)
[Ran11]. This is a family of iterative algorithms defined by two sequences of
_denoising functions_ $f_{t+1},g_{t}$ for each iteration $t\geq 0$. Under the
assumption that the design matrix is Gaussian, the joint empirical
distribution of the first $t$ iterates of GAMP converges to the law of $t$
jointly Gaussian random variables. (Here, $t$ is held fixed and the
convergence is with respect to the limit $n,d\to\infty$ with $n/d\to\delta$.)
The covariance structure of these jointly Gaussian random variables is
described by a set of recursions called _state evolution_.

The linear estimator is readily obtained via the iterate of GAMP run for one
step ($t=0$). For $t\geq 1$, we tailor the denoisers $(f_{t+1},g_{t})_{t\geq
1}$ so that the iterates of GAMP resemble a power method, which for large
enough $t$, gives the first (resp. second) eigenvector of the spectral matrix
$D$ (defined in Equation 2.5). We highlight a few challenges and our solutions
in implementing the above idea.

  * •

In a mixed GLM, since the observations $(y_{i})_{i\in[n]}$ are unlabeled
(i.e., it is unknown to the estimator whether each $y_{i}$ is generated from
the first or the second signal), estimating _both_ signals is more challenging
than estimating each one from an individual non-mixed GLM. However, the
existing state evolution result for GAMP [Ran11], [FVRS22, Sec. 4] is derived
for a non-mixed model, and only keeps track of the effect of a single signal.
We generalize the GAMP state evolution result to mixed GLMs, so that the state
evolution recursion tracks the effect of both signals. This result is derived
by reducing GAMP to an abstract AMP with matrix-valued iterates for which a
state evolution result has been established [JM13, FVRS22].

  * •

To study the limiting joint distribution of each signal and the corresponding
spectral estimator, we analyze a _pair_ of GAMPs with two different choices
for the denoisers $(f_{t+1},g_{t})_{t\geq 0}$. These choices ensure that the
GAMP equation for large $t$ becomes essentially an eigen-equation for the
first (resp. second) eigenvalue of the matrix $D$ defined in Equation 2.5. In
other words, the GAMP iteration effectively implements a power method.
However, power methods crucially require a spectral gap to converge to the
desired eigenvector. We show the existence of this spectral gap using tools
from random matrix theory, discussed below.

We stress that GAMP in this paper is used only as a tool for analysis and is
not part of the estimators. The actual estimators (spectral and linear) can be
computed by a combination of the following simple operations: _(i)_ applying a
component-wise nonlinearity, _(ii)_ matrix-vector/-matrix multiplication,
_(iii)_ computation of eigenvectors.

##### Eigenvalues via random matrix theory.

With the goal of proving the spectral gap needed by GAMP to approach its top
eigenvectors, in Theorem 5.1 we characterize the bulk and the outliers of the
limiting spectrum of $D$. The spectral analysis involves the following
challenges:

  * •

The matrix $D$ can be thought of as an instance of _spiked matrix model_. Its
structure is, however, more sophisticated than the canonical “signal plus
noise” model. Indeed, the potential spikes of $D$ result from two signals
through the composition of the link function $q$ and the spectral
preprocessing function ${\mathcal{T}}$.

  * •

The analysis of the limiting spectrum for non-mixed GLMs is provided in [LL20,
MM19]. In our mixed setting, applying the strategy of [LL20, MM19] to analyze
the spectrum of $D$ results in additional matrix terms which are hard to
bound.

The key idea is to decompose $D$ into the sum of two random matrices, $D_{1}$
and $D_{2}$, each consisting of the observations corresponding to the first
and second signal, respectively. When considered in isolation, $D_{1}$ and
$D_{2}$ are obtained from a non-mixed GLM with aspect ratio discounted by a
factor $\alpha$ and $1-\alpha$, respectively. Thanks to [LL20, MM19], their
limiting spectra are well understood. Now, the crucial observation is that, by
rotational invariance of the Gaussian design matrix and the fact that the
signals are independent and uniform on the sphere, $D_{1}$ and $D_{2}$ are
asymptotically free. This allows us to characterize their (free) sum using the
tools developed in [BBCF17]. Background on random matrix and free probability
theory is provided in Appendices F and G.

##### Optimal linear and spectral estimators.

The master theorem (Theorem 3.1) holds for arbitrary linear and spectral
preprocessing functions
${\mathcal{L}},{\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$ satisfying the
stated assumptions. Specializing Theorem 3.1 to linear and spectral estimators
alone and using the explicit formulas for their limiting overlaps (given in
Corollaries 3.3 and 3.5), we can find the _optimal_ preprocessing functions
${\mathcal{L}}^{*},{\mathcal{T}}_{1}^{*},{\mathcal{T}}_{2}^{*}\colon{\mathbb{R}}\to{\mathbb{R}}$
that maximize the limiting overlaps. This is done in Propositions 3.4 and 3.6
by casting the optimization problem as a variational problem and solving it
explicitly.

##  4 Numerical experiments

The experimental results in Figure 1 and Figures 3, 4 and 5 show that the
empirical performance of the various estimators (linear, spectral and
combined) closely match the asymptotic predictions. In all plots, the signal
dimension is $d=2000$, and the vertical and horizontal axes represent the
overlap and the aspect ratio $\delta$, respectively. The solid curves
correspond to the theoretical predictions whose analytic expressions can be
found in Section 3.3. Discrete points (little squares, triangles, asterisks,
etc.) are computed using synthetic data. Each of these points represents the
mean of $10$ i.i.d. trials together with error bars at $1$ standard deviation.

  * •

Figures 1(a) and 1(b) show numerical results for the recovery of the first and
second signal, respectively, from a noiseless linear regression model (i.e.,
the model in Equation 3.28 with $\sigma=0$) with mixing parameter
$\alpha=0.6$. We plot overlaps obtained via _(i)_ the optimal spectral
estimator in Equation 3.35, _(ii)_ the optimal linear estimator in Equation
3.30, _(iii)_ the Bayes-optimal linear combination of the estimators in _(i)_
and _(ii)_ (as per Corollary 3.2), _(iv)_ the spectral estimator proposed in
[YCS14] whose preprocessing
${\mathcal{T}}^{\mathrm{YCS}}\colon{\mathbb{R}}\to{\mathbb{R}}$ is:

$\displaystyle{\mathcal{T}}^{\mathrm{YCS}}(y)=\min\left\\{y^{2},10\right\\},$
(4.1)

and _(v)_ the spectral estimator proposed in [LAL19] whose preprocessing
function ${\mathcal{T}}^{\mathrm{LAL}}\colon{\mathbb{R}}\to{\mathbb{R}}$ is
given by

$\displaystyle{\mathcal{T}}^{\mathrm{LAL}}(y)=\max\left\\{1-\frac{1}{y^{2}},-10\right\\}.$
(4.2)

The estimators in both Equations 4.1 and 4.2 are truncated at $+10$ and $-10$,
respectively, in order to be able to compute our theoretical predictions. We
remark that choosing a larger value in magnitude for the truncation does _not_
experimentally lead to improved performance. This is justified by the fact
that such choices are not optimal for estimation from mixed models. Our
combined estimator in _(iii)_ and our optimal design of the spectral method in
_(i)_ offer a remarkable performance improvement with respect to existing
heuristic choices, such as _(iv)_ or _(v)_.

  * •

In Figure 3, we consider the recovery of both signals for noiseless mixed
linear regression (link function given by Equation 3.28 with $\sigma=0$),
using the spectral estimator with optimal preprocessing functions given by
Equation 3.35. Overlaps are plotted for two values of the mixing parameter
$\alpha\in\\{0.6,0.8\\}$. Precisely the same results (for both the simulation
at $d=2000$ and the asymptotic prediction) hold for noiseless phase retrieval,
as noted in Remark 3.10.

  * •

In Figure 4, we compare mixed linear regression and mixed phase retrieval
(link functions given by Equations 3.28 and 3.29, respectively), under their
respective optimal spectral estimators (see Equations 3.31 and 3.34). For each
model, we plot the overlap with the first signal for two different values of
the standard deviation of the noise $\sigma\in\\{0.8,1.5\\}$. In all four
cases, the mixing parameter is fixed to be $\alpha=0.8$. Though for $\sigma=0$
the curves for both models coincide, the gap between phase retrieval and
linear regression grows with $\sigma$, with the former model leading to
increasingly better performance.

  * •

In Figure 5, we consider the recovery of both signals for mixed linear
regression and mixed phase retrieval with mixing parameter $\alpha=0.6$ and
standard deviation of the noise $\sigma=1.5$. We observe that the overlaps for
linear regression are noticeably lower than that for phase retrieval, which
shows how the model noise makes the latter problem easier than the former for
spectral estimation.

Figure 3: Spectral estimators for noiseless mixed linear regression, with
mixing parameter $\alpha\in\\{0.6,0.8\\}$. Optimal spectral estimators given
by Equation 3.35 are used. Overlaps with both signals $x_{1}^{*},x_{2}^{*}$,
computed from simulation (“sim.”) and prediction (“pred.”), are plotted as a
function of the aspect ratio $\delta$. Same numerics apply to noiseless phase
retrieval (see Remark 3.10). Figure 4: Spectral estimators for mixed linear
regression and mixed phase retrieval, with mixing parameter $\alpha=0.8$ and
standard deviation of the noise $\sigma\in\\{0.8,1.5\\}$. Optimal spectral
estimators (Equations 3.31 and 3.34) are used. Overlap with the first signal
$x_{1}^{*}$, computed from simulation (“sim.”) and prediction (“pred.”), is
plotted as a function of the aspect ratio $\delta$. Figure 5: Spectral
estimators for mixed linear regression and mixed phase retrieval, with fixed
mixing parameter $\alpha=0.6$ and fixed standard deviation of the noise
$\sigma=1.5$. Optimal spectral estimators (Equations 3.31 and 3.34) are used.
Overlaps with both signals $x_{1}^{*},x_{2}^{*}$, computed from simulation
(“sim.”) and prediction (“pred.”), are plotted as a function of the aspect
ratio $\delta$.

##  5 Eigenvalues via random matrix theory

The characterization of the limiting joint law of spectral and linear
estimators in Theorem 3.1 is based on the analysis of a Generalized
Approximate Message Passing (GAMP) algorithm. The proof of convergence of GAMP
to the desired high-dimensional limit, whenever the conditions
$\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$ and/or
$\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$ are satisfied, crucially
relies on the existence of an eigengap. In this section, we derive the limits
of the top three eigenvalues of the matrix $D$ defined in Equation 2.5. This
result, stated as Theorem 5.1 below, is then used in Section 6 to prove
Theorem 3.1.

######  Theorem 5.1 (Eigenvalues).

Consider the setting of Section 2.2 and let Items (A1), (A2), (A3), (A4), (A5)
and (A7) hold. Then we have

$\displaystyle\lim_{d\to\infty}\lambda_{1}(D)$
$\displaystyle=\zeta(\lambda^{*}(\delta_{1});\delta),\quad\lim_{d\to\infty}\lambda_{2}(D)=\zeta(\lambda^{*}(\delta_{2});\delta),\quad\lim_{d\to\infty}\lambda_{3}(D)=\zeta(\overline{\lambda}(\delta);\delta),$

almost surely. Furthermore,

  1. 1.

If
$\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$,
then

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)>\zeta(\lambda^{*}(\delta_{2});\delta)>\zeta(\overline{\lambda}(\delta);\delta);$

  2. 2.

If
$\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{2})$,
then

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)>\zeta(\lambda^{*}(\delta_{2});\delta)=\zeta(\overline{\lambda}(\delta);\delta);$

  3. 3.

If
$\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})$,
then

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)=\zeta(\lambda^{*}(\delta_{2});\delta)=\zeta(\overline{\lambda}(\delta);\delta).$

######  Remark 5.1 (Phase transition for eigenvalues).

Theorem 5.1 shows a _phase transition_ phenomenon for the top three
eigenvalues of $D$: _(i)_ the top two eigenvalues escape from the bulk of $D$
if
$\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$;
_(ii)_ only the largest eigenvalue escapes from the bulk if
$\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{2})$;
_(iii)_ no outlier eigenvalue exists if
$\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})$.
See Figure 2. In words, the condition
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$ is _necessary and
sufficient_ for the $i$-th eigenvalue to escape the bulk of the spectrum. This
provides an additional piece of evidence (see also Remark 3.7) suggesting that
such condition is also necessary and sufficient for the corresponding
eigenvector to have non-vanishing overlap with the signal. In fact, phase
transitions in the behavior of eigenvalues typically correspond to phase
transitions in the behavior of the related eigenvectors, see e.g. [BGN11,
BGN12, MM19, LL20].

######  Remark 5.2 (Eigenvalues for $\alpha=1/2$).

Similar results to Theorem 5.1 hold for $\alpha=1/2$. In this case, the limits
of the first and second eigenvalues of $D$ coincide and equal
$\zeta(\lambda^{*}(\delta/2);\delta)$. The limit of the right edge of the bulk
of $D$ does not depend on $\alpha$ and remains the same
($\overline{\lambda}(\delta)$) as in Theorem 5.1. Therefore, we get two cases:
_(i)_ if $\lambda^{*}(\delta/2)>\overline{\lambda}(\delta)$, the top two
eigenvalues of $D$ are repeated outliers; otherwise, _(ii)_
$\overline{\lambda}(\delta)\geq\lambda^{*}(\delta/2)$ and there is no outlier
eigenvalue in the limiting spectrum of $D$.

######  Remark 5.3 (Explicit formulas).

By the definition of $\zeta(\lambda;\delta)$ (cf. Equation 3.4), we can write
the limits of the eigenvalues in the following more explicit form, which will
be convenient in Section 6:

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)$
$\displaystyle=\begin{cases}\lambda^{*}(\delta_{1})\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right]\right),&\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)\\\
\overline{\lambda}(\delta)\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\overline{\lambda}(\delta)-Z}\right]\right),&\lambda^{*}(\delta_{1})\leq\overline{\lambda}(\delta)\end{cases},$
$\displaystyle\zeta(\lambda^{*}(\delta_{2});\delta)$
$\displaystyle=\begin{cases}\lambda^{*}(\delta_{2})\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\lambda^{*}(\delta_{2})-Z}\right]\right),&\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)\\\
\overline{\lambda}(\delta)\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\overline{\lambda}(\delta)-Z}\right]\right),&\lambda^{*}(\delta_{2})\leq\overline{\lambda}(\delta)\end{cases},$
$\displaystyle\zeta(\overline{\lambda}(\delta);\delta)$
$\displaystyle=\overline{\lambda}(\delta)\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\overline{\lambda}(\delta)-Z}\right]\right).$

###### Proof of Theorem 5.1..

The proof is divided into three steps. Specifically, we first write $D$ as the
sum of two asymptotically free spiked random matrices in Section 5.1. Then,
the limit of $\lambda_{3}(D)$ is determined in Lemma 5.2. Finally, the limits
of $\lambda_{1}(D),\lambda_{2}(D)$ and the monotonicity properties of the
limiting eigenvalues in Items 1, 2 and 3 are given by Lemma 5.3. ∎

###  5.1 Reduction to free additive convolution

To begin with, assume for notational convenience that $\eta_{i}=1$ for $1\leq
i\leq n_{1}$ and $\eta_{i}=0$ for $n_{1}+1\leq i\leq n$, for some $0\leq
n_{1}\leq n$. Let $n_{2}=n-n_{1}$. We have almost surely
$n_{1}/d\to\delta_{1},n_{2}/d\to\delta_{2}$ since the mixing variables are
independent and Bernoulli-distributed with mean $\alpha$. We can write the
matrices of interest in block forms

$\displaystyle A$ $\displaystyle=\begin{bmatrix}A_{1}\\\
A_{2}\end{bmatrix},\quad T=\begin{bmatrix}T_{1}&0_{n_{1}\times n_{2}}\\\
0_{n_{2}\times n_{1}}&T_{2}\end{bmatrix},$ (5.1)

where $A_{1}\in{\mathbb{R}}^{n_{1}\times d},A_{2}\in{\mathbb{R}}^{n_{2}\times
d}$ and $T_{1}\in{\mathbb{R}}^{n_{1}\times
n_{1}},T_{2}\in{\mathbb{R}}^{n_{2}\times n_{2}}$. We also let
${\underline{\varepsilon}}_{1}=(\varepsilon_{1},\cdots,\varepsilon_{n_{1}})$
and
${\underline{\varepsilon}}_{2}=(\varepsilon_{n_{1}+1},\cdots,\varepsilon_{n})$.
Then,

$\displaystyle A^{\top}TA$
$\displaystyle=\begin{bmatrix}A_{1}^{\top}&A_{2}^{\top}\end{bmatrix}\begin{bmatrix}T_{1}&0_{n_{1}\times
n_{2}}\\\ 0_{n_{2}\times n_{1}}&T_{2}\end{bmatrix}\begin{bmatrix}A_{1}\\\
A_{2}\end{bmatrix}=A_{1}^{\top}T_{1}A_{1}+A_{2}^{\top}T_{2}A_{2}.$ (5.2)

Note that, for $i\in\\{1,2\\}$,

$\displaystyle A_{i}^{\top}T_{i}A_{i}$
$\displaystyle=A_{i}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{i}x_{i}^{*},{\underline{\varepsilon}}_{i})))A_{i}.$

Since $A_{1},x_{1}^{*},{\underline{\varepsilon}}_{1}$ and
$A_{2},x_{2}^{*},{\underline{\varepsilon}}_{2}$ are mutually independent, we
know that $A_{1}^{\top}T_{1}A_{1}$ is independent of $A_{2}^{\top}T_{2}A_{2}$.
However, $A_{1}$ and $T_{1}$ are _not_ independent, neither are $A_{2}$ and
$T_{2}$.

Let $O\sim\operatorname{Haar}({\mathbb{O}}(d))$ be a matrix sampled uniformly
from the orthogonal group ${\mathbb{O}}(d)$ and independent of everything
else. Then, we have

$\displaystyle A_{1}^{\top}T_{1}A_{1}+A_{2}^{\top}T_{2}A_{2}$
$\displaystyle=A_{1}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{1}x_{1}^{*},{\underline{\varepsilon}}_{1})))A_{1}+A_{2}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{2}x_{2}^{*},{\underline{\varepsilon}}_{2})))A_{2}$
$\displaystyle\overset{\mathrm{d}}{=}A_{1}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{1}x_{1}^{*},{\underline{\varepsilon}}_{1})))A_{1}+(A_{2}O)^{\top}\mathrm{diag}({\mathcal{T}}(q((A_{2}O)x_{2}^{*},{\underline{\varepsilon}}_{2})))(A_{2}O)$
(5.3)
$\displaystyle=A_{1}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{1}x_{1}^{*},{\underline{\varepsilon}}_{1})))A_{1}+O^{\top}A_{2}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{2}(Ox_{2}^{*}),{\underline{\varepsilon}}_{2})))A_{2}O$
$\displaystyle\overset{\mathrm{d}}{=}A_{1}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{1}x_{1}^{*},{\underline{\varepsilon}}_{1})))A_{1}+O^{\top}A_{2}^{\top}\mathrm{diag}({\mathcal{T}}(q(A_{2}x_{2}^{*},{\underline{\varepsilon}}_{2})))A_{2}O$
(5.4) $\displaystyle=A_{1}^{\top}T_{1}A_{1}+O^{\top}A_{2}^{\top}T_{2}A_{2}O.$

Equation 5.3 follows from the independence of $A_{1},A_{2}$, and from the
rotational invariance of isotropic Gaussians. Equation 5.4 follows since $O$
and $Ox_{2}^{*}$ are independent if
$O\sim\operatorname{Haar}({\mathbb{O}}(d))$ and
$x_{2}^{*}\sim\operatorname{Unif}({\mathbb{S}}^{d-1})$. In this step, we
crucially use the assumption that $x_{1}^{*}$ and $x_{2}^{*}$ are independent
and each uniformly distributed over ${\mathbb{S}}^{d-1}$.

Now $A_{1}^{\top}T_{1}A_{1}$ and $O^{\top}A_{2}^{\top}T_{2}A_{2}O$ are
_asymptotically free_. (See Appendix F for a definition of asymptotic freeness
and a primer on free probability theory.) Therefore, we can study
$A_{1}^{\top}T_{1}A_{1}$ and $A_{2}^{\top}T_{2}A_{2}$ separately by using
existing results, and then apply tools from free probability theory to
characterize their sum. In particular, to understand $A_{i}^{\top}T_{i}A_{i}$
for $i\in\\{1,2\\}$, we use a theorem from [MM19] (transcribed in Theorem
I.6); to understand $A_{1}^{\top}T_{1}A_{1}+A_{2}^{\top}T_{2}A_{2}$, we use a
theorem from [BBCF17] (transcribed in Theorem G.1).

###  5.2 Right edge of the bulk of $D$

Before proceeding to the analysis, let us introduce some more notation. Let

$\displaystyle D_{1}=\frac{1}{n}A_{1}^{\top}T_{1}A_{1},\quad
D_{2}=\frac{1}{n}A_{2}^{\top}T_{2}A_{2}.$ (5.5)

Therefore $D=D_{1}+D_{2}$ according to Equation 5.2. We first calculate the
limiting value of the right edge of the bulk of the spectrum of $D$.

######  Lemma 5.2.

Consider the setting of Section 2.2 and let Items (A1), (A2), (A3), (A4), (A5)
and (A7) hold. Denote by $\mu_{D}$ the empirical spectral distribution of $D$.
Then

$\displaystyle\lim_{d\to\infty}\sup\operatorname{supp}(\mu_{D})$
$\displaystyle=\frac{1}{\delta}\cdot
s^{-1}_{\mu_{1}\boxplus\mu_{2}}(-1/\overline{\lambda}(\delta)),$ (5.6)

almost surely, where $\overline{\lambda}(\delta)$ is the solution to

$\displaystyle\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\delta)-Z}\right)^{2}\right]$
$\displaystyle=\frac{1}{\delta}$ (5.7)

and the function $s_{\mu_{1}\boxplus\mu_{2}}^{-1}$ is defined as

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}^{-1}(z)$
$\displaystyle=-\frac{1}{z}+\delta\,\mathbb{E}\left[\frac{Z}{1+zZ}\right].$
(5.8)

######  Remark 5.4.

The function $s_{\mu_{1}\boxplus\mu_{2}}^{-1}$ is in fact the inverse
_Stieltjes transform_ of the _free additive convolution_ of the limiting
spectral distributions $\mu_{1}$ of $\frac{n}{d}D_{1}$ and $\mu_{2}$ of
$\frac{n}{d}D_{2}$. The function $s_{\mu_{1}\boxplus\mu_{2}}^{-1}(\lambda)$ is
tightly connected to the function $\psi(\lambda;\delta)$ defined in Equation
3.2. Indeed, $\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$ is precisely
$\psi(\lambda;\delta)$. We note also that the parameter
$\overline{\lambda}(\delta)$ defined in Equation 5.7 is the same as that
defined through Equation 3.3. (See Lemma I.3.) The connection shall become
more transparent in the proof below.

###### Proof of Lemma 5.2.

First note that the scaling factor $\frac{1}{d}$ in Theorem I.7 (which will be
used momentarily) is different from our scaling $\frac{1}{n}$ in the
definition of $D$ (cf. Equation 2.5). We therefore consider
$\widetilde{D}=\frac{1}{d}A^{\top}TA$ for the convenience of applying Theorem
I.7. All results regarding the matrix $\widetilde{D}$ can be translated to $D$
by inserting a factor $\frac{d}{n}\to\frac{1}{\delta}$ at proper places, since
$D=\frac{d}{n}\widetilde{D}$.

Let

$\displaystyle\widetilde{D}_{1}=\frac{1}{d}A_{1}^{\top}T_{1}A_{1},\quad\widetilde{D}_{2}=\frac{1}{d}A_{2}^{\top}T_{2}A_{2}.$
(5.9)

By Equation 5.2, $\widetilde{D}=\widetilde{D}_{1}+\widetilde{D}_{2}$. Let
$\mu_{1}$ and $\mu_{2}$ be the limiting spectral distributions of
$\widetilde{D}_{1}$ and $\widetilde{D}_{2}$, respectively, as
$n_{1},n_{2},d\to\infty$ with $n_{1}/d\to\delta_{1}=\alpha\delta$ and
$n_{2}/d\to\delta_{2}=(1-\alpha)\delta$. As argued in Section 5.1,
$\widetilde{D}_{1}$ and $\widetilde{D}_{2}$ are asymptotically free. Hence,
the limiting spectral distribution of $\widetilde{D}$ is given by the free
additive convolution of $\mu_{1},\mu_{2}$, denoted by $\mu_{1}\boxplus\mu_{2}$
[Voi91, Spe93]. It remains to compute
$\sup\operatorname{supp}(\mu_{1}\boxplus\mu_{2})$.

A careful inspection of the proof of [MM19, Lemma 2] shows that the bulk of
the spectrum of $\widetilde{D}_{i}$ (for $i\in\\{1,2\\}$ respectively), i.e.,
$\lambda_{2}(\widetilde{D}_{i})\geq\cdots\geq\lambda_{d}(\widetilde{D}_{i})$,
interlaces the spectrum of
$E_{i}\coloneqq\frac{1}{d}\widetilde{A}_{i}^{\top}T_{i}\widetilde{A}_{i}\in{\mathbb{R}}^{(d-1)\times(d-1)}$
for $i\in\\{1,2\\}$, respectively. Specifically,

$\displaystyle\lambda_{1}(E_{i})\geq\lambda_{2}(\widetilde{D}_{i})\geq\lambda_{2}(E_{i})\geq\lambda_{3}(\widetilde{D}_{i})\geq\cdots\geq\lambda_{d-2}(E_{i})\geq\lambda_{d-1}(\widetilde{D}_{i})\geq\lambda_{d-1}(E_{i})\geq\lambda_{d}(\widetilde{D}_{i}).$
(5.10)

Here,
$T_{i}=\mathrm{diag}({\mathcal{T}}(q(A_{i}x_{i}^{*},{\underline{\varepsilon}}_{i})))$
(recall Equation 5.1) and
$\widetilde{A}_{i}\in{\mathbb{R}}^{n_{i}\times(d-1)}$ is an _independent_
matrix with i.i.d. ${\mathcal{N}}(0,1)$ entries. In particular, $T_{i}$ and
$\widetilde{A}_{i}$ are independent. For the convenience of applying Theorem
I.7, we also define, for $i\in\\{1,2\\}$,
$\widetilde{E}_{i}\coloneqq\frac{1}{d-1}\widetilde{A}_{i}^{\top}T_{i}\widetilde{A}_{i}$.
Note that $E_{i}=\frac{d-1}{d}\widetilde{E}_{i}$. Also
$n_{1}/(d-1)\to\delta_{1}$ and $n_{2}/(d-1)\to\delta_{2}$.

Since each $y_{i}$ (for $1\leq i\leq n$) is i.i.d., the limiting spectral
distributions of $T_{1}$ and $T_{2}$ are in fact the same and both equal the
law of $Z$. Thus, Theorem I.7 provides us with a characterization of the
limiting spectral distribution of $\widetilde{E}_{i}$:

$\displaystyle\mu_{\widetilde{E}_{1}}$
$\displaystyle\to\widetilde{\mu}_{1},\quad\mu_{\widetilde{E}_{2}}\to\widetilde{\mu}_{2},$

weakly as $n_{1},n_{2},d\to\infty$ with
$n_{1}/(d-1)\to\delta_{1},n_{2}/(d-1)\to\delta_{2}$. Furthermore, the limiting
spectral distributions admit the following explicit description through the
inverse Stieltjes transform (see Equation F.1 in Appendix F for the definition
of _Stieltjes transform_ of a distribution):

$\displaystyle s_{\widetilde{\mu}_{1}}^{-1}(z)$
$\displaystyle=-\frac{1}{z}+\delta_{1}\mathbb{E}\left[\frac{Z}{1+zZ}\right],\quad
s_{\widetilde{\mu}_{2}}^{-1}(z)=-\frac{1}{z}+\delta_{2}\mathbb{E}\left[\frac{Z}{1+zZ}\right].$
(5.11)

In view of the scaling factor $\frac{d-1}{d}\to 1$, the limiting spectral
distributions of $E_{1},E_{2}$ are also given by
$\widetilde{\mu}_{1},\widetilde{\mu}_{2}$, respectively. Recall that the bulks
of the spectra of $\widetilde{D}_{1},\widetilde{D}_{2}$ interlace the spectra
of $E_{1},E_{2}$, respectively (cf. Equation 5.10). Since $\widetilde{D}_{1}$
and $\widetilde{D}_{2}$ can each have at most one outlier eigenvalue (cf.
Theorem I.6), the limiting spectral distributions $\mu_{1},\mu_{2}$ of
$\widetilde{D}_{1},\widetilde{D}_{2}$, respectively, are the same as
$\widetilde{\mu}_{1},\widetilde{\mu}_{2}$ whose inverse Stieltjes transforms
are shown in Equation 5.11.

The following two well-known facts in free probability theory (cf. Equations
F.2 and F.3 in Appendix F)

$\displaystyle R_{\mu_{1}\boxplus\mu_{2}}(z)$
$\displaystyle=R_{\mu_{1}}(z)+R_{\mu_{2}}(z),$

and

$\displaystyle R_{\mu}(z)$ $\displaystyle=s_{\mu}^{-1}(-z)-\frac{1}{z},$

then allows us to compute $s_{\mu_{1}\boxplus\mu_{2}}^{-1}$. Indeed

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}^{-1}(z)$
$\displaystyle=s_{\mu_{1}}^{-1}(z)+s_{\mu_{2}}^{-1}(z)+\frac{1}{z}=s_{\widetilde{\mu}_{1}}^{-1}(z)+s_{\widetilde{\mu}_{2}}^{-1}(z)+\frac{1}{z}=-\frac{1}{z}+\delta\mathbb{E}\left[\frac{Z}{1+zZ}\right].$
(5.12)

Given $s_{\mu_{1}\boxplus\mu_{2}}^{-1}$, one can calculate
$\sup\operatorname{supp}(\mu_{1}\boxplus\mu_{2})$ which is in turn the
limiting value of $\sup\operatorname{supp}(\mu_{\widetilde{D}})$, where
$\mu_{\widetilde{D}}$ denotes the empirical spectral distribution of
$\widetilde{D}$. This can be accomplished thanks to the results in [BY12,
Lemma 3.1] (see also [SC95, Sec. 4]):

$\displaystyle\lim_{d\to\infty}\sup\operatorname{supp}(\mu_{\widetilde{D}})$
$\displaystyle=\sup\operatorname{supp}(\mu_{1}\boxplus\mu_{2})$ (5.13)
$\displaystyle=\min_{\lambda>\sup\operatorname{supp}(Z)}s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$
$\displaystyle=\min_{\lambda>\sup\operatorname{supp}(Z)}\lambda+\delta\mathbb{E}\left[\frac{Z}{1-Z/\lambda}\right].$
(5.14)

The convergence in Equation 5.13 holds almost surely since

$\displaystyle\mu_{\widetilde{D}}$
$\displaystyle=\mu_{\widetilde{D}_{1}+\widetilde{D}_{2}}\to\mu_{1}\boxplus\mu_{2}$

weakly [Voi91, Spe93]. To solve the minimization problem in Equation 5.14, we
observe that the function $s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$ is
intimately related to $\psi(\lambda;\delta)$ defined in Equation 3.2:

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$
$\displaystyle=\lambda+\delta\mathbb{E}\left[\frac{Z}{1-Z/\lambda}\right]=\delta\lambda\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\lambda-Z}\right]\right)=\delta\cdot\psi(\lambda;\delta).$

Since $\psi(\lambda;\delta)$ is convex in the first argument (cf. Lemma I.2),
so is $s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$ as a function of
$\lambda$. As a result, the minimizer $\overline{\lambda}(\delta)$ of the
above minimization problem (Equation 5.14) is the critical point of
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$. That is,

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}\lambda}s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)\right|_{\lambda=\overline{\lambda}(\delta)}$
$\displaystyle=1-\delta\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\delta)-Z}\right)^{2}\right]=0,$

i.e., $\overline{\lambda}(\delta)$ is the solution to the following equation

$\displaystyle\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\delta)-Z}\right)^{2}\right]$
$\displaystyle=\frac{1}{\delta}.$ (5.15)

The minimum value of the above minimization problem (Equation 5.14) is
therefore

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\overline{\lambda}(\delta))$
$\displaystyle=\overline{\lambda}(\delta)\left(1+\delta\mathbb{E}\left[\frac{Z}{\overline{\lambda}(\delta)-Z}\right]\right).$
(5.16)

At this point, we have successfully computed the limiting value of
$\sup\operatorname{supp}(\mu_{\widetilde{D}})$. However, recall that the
original matrix we are interested in is
$D=\frac{d}{n}(\widetilde{D}_{1}+\widetilde{D}_{2})$. Therefore,

$\displaystyle\lim_{d\to\infty}\sup\operatorname{supp}(\mu_{D})$
$\displaystyle=\lim_{d\to\infty}\frac{d}{n}\sup\operatorname{supp}(\mu_{\widetilde{D}})=\overline{\lambda}(\delta)\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\overline{\lambda}(\delta)-Z}\right]\right),$

where $\overline{\lambda}(\delta)$ satisfies Equation 5.15. This concludes the
proof. ∎

###  5.3 Outlier eigenvalues of $D$

Finally, we need to understand the outliers in the spectrum of $D$.

######  Lemma 5.3.

Consider the setting of Section 2.2 and let Items (A1), (A2), (A3), (A4), (A5)
and (A7) hold. Let the function $s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$
be given by Equation 5.8. Then, the following statements hold.

  1. 1.

$\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})$;

  2. 2.

If
$\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$,
then
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{1}))>s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{2}))$;

  3. 3.

For $i\in\\{1,2\\}$, if $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$,
then
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{i}))>s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\overline{\lambda}(\delta))$;

  4. 4.

We have that, almost surely,

$\displaystyle\lim_{d\to\infty}\lambda_{1}(D)$
$\displaystyle=\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\max\\{\lambda^{*}(\delta_{1}),\overline{\lambda}(\delta)\\}),$
(5.17) $\displaystyle\lim_{d\to\infty}\lambda_{2}(D)$
$\displaystyle=\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\max\\{\lambda^{*}(\delta_{2}),\overline{\lambda}(\delta)\\}),$
(5.18) $\displaystyle\lim_{d\to\infty}\lambda_{3}(D)$
$\displaystyle=\frac{1}{\delta}\cdot\sup\operatorname{supp}(\mu_{1}\boxplus\mu_{2})=\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\overline{\lambda}(\delta)).$ (5.19)

######  Remark 5.5.

Recalling the definition
$\zeta(\lambda;\delta)=\psi(\max\\{\lambda,\overline{\lambda}(\delta)\\};\delta)$
(cf. Equation 3.4) and the relation $\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)=\psi(\lambda;\delta)$ (cf. Remark
5.4), we can write the limiting values of
$\lambda_{1}(D),\lambda_{2}(D),\lambda_{3}(D)$ in Equations 5.17, 5.18 and
5.19 as

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)\geq\zeta(\lambda^{*}(\delta_{2});\delta)\geq\zeta(\overline{\lambda}(\delta);\delta),$
(5.20)

respectively. To see why the above chain of inequalities holds, note that by
Item 3 of Lemma 5.3,
$\zeta(\lambda^{*}(\delta_{i});\delta)>\zeta(\overline{\lambda}(\delta);\delta)$
if $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$ and
$\zeta(\lambda^{*}(\delta_{i});\delta)=\zeta(\overline{\lambda}(\delta);\delta)$
otherwise. So

$\displaystyle\zeta(\lambda^{*}(\delta_{i});\delta)\geq\zeta(\overline{\lambda}(\delta);\delta)$
(5.21)

is always true for $i\in\\{1,2\\}$. Also, by Item 2 of Lemma 5.3,
$\zeta(\lambda^{*}(\delta_{1});\delta)\geq\zeta(\lambda^{*}(\delta_{2});\delta)>\zeta(\overline{\lambda}(\delta);\delta)$
if
$\lambda^{*}(\delta_{1})\geq\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$.
If
$\lambda^{*}(\delta_{1})\geq\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{2})$,
$\zeta(\lambda^{*}(\delta_{2});\delta)=\zeta(\overline{\lambda}(\delta);\delta)\leq\zeta(\lambda^{*}(\delta_{1});\delta)$
by Equation 5.21. If
$\overline{\lambda}(\delta)\geq\lambda^{*}(\delta_{1})\geq\lambda^{*}(\delta_{2})$,
$\zeta(\lambda^{*}(\delta_{1});\delta)=\zeta(\lambda^{*}(\delta_{2});\delta)=\zeta(\overline{\lambda}(\delta);\delta)$.
All cases have been exhausted in light of Item 1 of Lemma 5.3. In any case,

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta)\geq\zeta(\lambda^{*}(\delta_{2});\delta)$
(5.22)

holds. Equation 5.20 then follows from Equations 5.21 and 5.22.

###### Proof of Lemma 5.3..

The proof is divided into three parts. We first explicitly evaluate the
theoretical predictions of the limiting values of the top three eigenvalues of
$D$ given by Equation G.1 of Theorem G.1. The convergence of the outlier
eigenvalues and the right edge of the bulk to the respective predictions is
then formally justified in the second part. Finally, several properties
concerning the spectral threshold and the limiting eigenvalues are proved in
the third part.

##### Limiting eigenvalues.

To understand the outlier eigenvalues of $D=D_{1}+D_{2}$, we need to first
understand the outlier eigenvalues of $D_{1}$ and $D_{2}$ individually. The
latter quantities have been characterized by [MM19, Lemma 2] (transcribed in
Theorem I.6). To calibrate the scaling, let us define

$\displaystyle
D_{1}^{\prime}\coloneqq\frac{1}{n_{1}}A_{1}^{\top}T_{1}A_{1},\quad
D_{2}^{\prime}\coloneqq\frac{1}{n_{2}}A_{2}^{\top}T_{2}A_{2}.$

Theorem I.6 then applies to the above matrices $D_{1}^{\prime},D_{2}^{\prime}$
and implies that each of $D_{1}^{\prime}$ and $D_{2}^{\prime}$ has a potential
outlier eigenvalue $\lambda_{1}(D_{1}^{\prime})$ and
$\lambda_{1}(D_{2}^{\prime})$, respectively. As $n_{1},n_{2},d\to\infty$ with
$n_{1}/d\to\delta_{1}$ and $n_{2}/d\to\delta_{2}$, they converge almost surely
to the following limiting values:

$\displaystyle\lim_{d\to\infty}\lambda_{1}(D_{1}^{\prime})$
$\displaystyle=\zeta(\lambda^{*}(\delta_{1});\delta_{1}),\quad\lim_{d\to\infty}\lambda_{1}(D_{2}^{\prime})=\zeta(\lambda^{*}(\delta_{2});\delta_{2}),$

where $\lambda^{*}(\delta_{1})$ and $\lambda^{*}(\delta_{2})$ are the
solutions to

$\displaystyle\zeta(\lambda^{*}(\delta_{1});\delta_{1})$
$\displaystyle=\varphi(\lambda^{*}(\delta_{1})),\quad\zeta(\lambda^{*}(\delta_{2});\delta_{2})=\varphi(\lambda^{*}(\delta_{2})),$

respectively. For $i\in\\{1,2\\}$, let us assume that
$\lambda_{1}(D_{i}^{\prime})$ is indeed an outlier eigenvalue of
$D_{i}^{\prime}$, that is, its limiting value
$\zeta(\lambda^{*}(\delta_{i});\delta_{i})$ lies outside the bulk of the
limiting spectrum of $D_{i}^{\prime}$. According to Theorem I.6, this happens
if and only if $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta_{i})$. In
this case, the limiting value of the outlier eigenvalue can be written more
explicitly as

$\displaystyle\zeta(\lambda^{*}(\delta_{i});\delta_{i})=\psi(\lambda^{*}(\delta_{i});\delta_{i})=\lambda^{*}(\delta_{i})\left(\frac{1}{\delta_{i}}+\mathbb{E}\left[\frac{Z}{\lambda^{*}(\delta_{i})-Z}\right]\right),$
(5.23)

where $\lambda^{*}(\delta_{i})$ is the solution to

$\displaystyle\mathbb{E}\left[\frac{Z(G^{2}-1)}{\lambda^{*}(\delta_{i})-Z}\right]$
$\displaystyle=\frac{1}{\delta_{i}}.$ (5.24)

Let us first translate the above result (i.e., Equations 5.23 and 5.24)
regarding $D_{1}^{\prime},D_{2}^{\prime}$ to
$\widetilde{D}_{1},\widetilde{D}_{2}$ defined in Equation 5.9. Since
$n_{1}/d\to\delta_{1},n_{2}/d\to\delta_{2}$ and
$\widetilde{D}_{1}=\frac{n_{1}}{d}D_{1}^{\prime},\widetilde{D}_{2}=\frac{n_{2}}{d}D_{2}^{\prime}$,
we have that, almost surely,

$\displaystyle\lim_{d\to\infty}\lambda_{1}(\widetilde{D}_{i})$
$\displaystyle=\delta_{i}\lambda^{*}(\delta_{i})\left(\frac{1}{\delta_{i}}+\mathbb{E}\left[\frac{Z}{\lambda^{*}(\delta_{i})-Z}\right]\right)=\lambda^{*}(\delta_{i})\left(1+\delta_{i}\mathbb{E}\left[\frac{Z}{\lambda^{*}(\delta_{i})-Z}\right]\right)\eqqcolon\theta_{i},$

where we have denoted the limiting value of $\lambda_{1}(\widetilde{D}_{i})$
by $\theta_{i}$. In view of the definition of $s_{\mu_{i}}^{-1}$ in Equation
5.11, we recognize that

$\displaystyle\theta_{i}=s_{\mu_{i}}^{-1}(-1/\lambda^{*}(\delta_{i})).$ (5.25)

Provided with the individual outlier of $\widetilde{D}_{i}$ (cf. Equation
5.25), we now invoke [BBCF17, Theorem 2.1] (transcribed in Theorem G.1) to
determine how an outlier of $\widetilde{D}_{i}$ is mapped to the spectrum of
$\widetilde{D}$ by the free additive convolution. According to Theorem G.1,
the limiting value, denoted by $\rho_{i}$, of the potential outlier of
$\widetilde{D}=\widetilde{D}_{1}+\widetilde{D}_{2}$ resulting from
$\theta_{i}$ is given by

$\displaystyle\rho_{i}\coloneqq w_{i}^{-1}(\theta_{i}),$ (5.26)

where $w_{1},w_{2}$ are the pair of _subordination functions_ associated with
the pair of distributions $\mu_{1},\mu_{2}$.

As the name suggests, $w_{1},w_{2}$ enjoy the following subordination property
(cf. [BBCF17, Sec. 3.4.1] or Equation F.4 in Appendix F):

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}(z)$
$\displaystyle=s_{\mu_{1}}(w_{1}(z))=s_{\mu_{2}}(w_{2}(z))=\frac{1}{z-(w_{1}(z)+w_{2}(z))}.$
(5.27)

To understand the value of $\rho_{i}=w_{i}^{-1}(\theta_{i})$ (cf. Equation
5.26), let us compute

$\displaystyle s_{\mu_{1}\boxplus\mu_{2}}(w_{i}^{-1}(\theta_{i}))$
$\displaystyle=s_{\mu_{i}}(\theta_{i})=-1/\lambda^{*}(\delta_{i}).$ (5.28)

The first equality is by the subordination property (Equation 5.27) and the
second one by the observation in Equation 5.25. Equation 5.28 then gives

$\displaystyle\rho_{i}$
$\displaystyle=w_{i}^{-1}(\theta_{i})=s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{i})).$
(5.29)

To translate the result in Equation 5.29 regarding $\widetilde{D}$ to $D$
defined in Equation 5.5, we simply note that $D=\frac{d}{n}\widetilde{D}$ and
$d/n\to 1/\delta$. Therefore, the limiting eigenvalue of $D$ resulting from
the outlier eigenvalue of $D_{i}$ is given by

$\displaystyle\frac{1}{\delta}\cdot\rho_{i}=\frac{1}{\delta}\cdot
s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{i})),$ (5.30)

almost surely.

##### Convergence of eigenvalues.

We then formally justify that the right edge of the bulk and the outlier
eigenvalues of $D$ indeed converge to the theoretical predictions in Equations
5.6 and 5.30, respectively, as $d\to\infty$, therefore confirming the validity
of the latter formulas. Let
${\mathcal{K}}_{0}\coloneqq\operatorname{supp}(\mu_{1}\boxplus\mu_{2})$. For
$i\in\\{1,2\\}$, let ${\mathcal{K}}_{i}$ be the singleton set $\\{\rho_{i}\\}$
if $\theta_{i}\notin\operatorname{supp}(\mu_{i})$ and $\emptyset$ otherwise.
Let
${\mathcal{K}}\coloneqq{\mathcal{K}}_{0}\cup{\mathcal{K}}_{1}\cup{\mathcal{K}}_{2}$.
Then the first statement of [BBCF17, Theorem 2.1] (transcribed in Theorem G.1)
guarantees that for any $\varepsilon>0$,

$\displaystyle\Pr\left[\exists d_{0},\,\forall
d>d_{0},\,\\{\lambda_{i}(\widetilde{D})\\}_{i=1}^{d}\subset{\mathcal{K}}_{\varepsilon}\right]$
$\displaystyle=1,$ (5.31)

where ${\mathcal{K}}_{\varepsilon}$ denotes the _$\varepsilon$ -enlargement_
of ${\mathcal{K}}$, i.e.,

$\displaystyle{\mathcal{K}}_{\varepsilon}\coloneqq\left\\{\rho\in{\mathbb{R}}:\inf_{\rho^{\prime}\in{\mathcal{K}}}|\rho-\rho^{\prime}|\leq\varepsilon\right\\}.$

In words, Equation 5.31 says that almost surely for every sufficiently large
dimension $d$, the spectrum of $\widetilde{D}$ is contained in an arbitrarily
small neighbourhood of ${\mathcal{K}}$. Furthermore, suppose
$\rho\in{\mathcal{K}}_{1}\cup{\mathcal{K}}_{2}$ and
$\rho\notin{\mathcal{K}}_{0}$, that is, $\rho$ is an outlier in the limiting
spectrum of $\widetilde{D}$. Assume also that $\varepsilon>0$ is sufficiently
small so that
$(\rho-2\varepsilon,\rho+2\varepsilon)\cap{\mathcal{K}}=\\{\rho\\}$. Then

$\displaystyle\Pr\left[\exists d_{0},\,\forall
d>d_{0},\,\left|\\{\lambda_{i}(\widetilde{D})\\}_{i=1}^{d}\cap(\rho-2\varepsilon,\rho+2\varepsilon)\right|=\mathds{1}{\left\\{w_{1}(\rho)=\theta_{1}\right\\}}+\mathds{1}{\left\\{w_{2}(\rho)=\theta_{2}\right\\}}\right]$
$\displaystyle=1.$ (5.32)

In words, Equation 5.32 says that almost surely for every sufficiently large
dimension $d$, the outlier $\theta_{1}$ (resp. $\theta_{2}$) in the limiting
spectrum of $\widetilde{D}_{1}$ (resp. $\widetilde{D}_{2}$) is mapped to
$w_{1}^{-1}(\theta_{1})$ (resp. $w_{2}^{-1}(\theta_{2})$) in the limiting
spectrum of $\widetilde{D}$. Since $D,D_{1},D_{2}$ and
$\widetilde{D},\widetilde{D}_{1},\widetilde{D}_{2}$ only differ by a $\delta$
factor, similar statements hold true for $D,D_{1},D_{2}$ as well.

Combining Equations 5.30, 5.31 and 5.32 yields Equations 5.17, 5.18 and 5.19
in Item 4 of Lemma 5.3.

##### Properties of spectral threshold and limiting eigenvalues.

We identify under what condition $\rho_{i}=w_{i}^{-1}(\theta_{i})$ is an
outlier in the limiting spectrum of $\widetilde{D}$. For this to be the case,
$\theta_{i}$ is necessarily an outlier in the limiting spectrum of
$\widetilde{D}_{i}$, which is assumed in the preceding derivations. As Theorem
I.6 guaranteed, a sufficient and necessary condition for this event is
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta_{i})$. Under the free
additive convolution, the outlier $\theta_{i}$ of $\widetilde{D}_{i}$ is then
mapped to $w_{i}^{-1}(\theta_{i})$. Let us compare $w_{i}^{-1}(\theta_{i})$
with $\sup\operatorname{supp}(\mu_{1}\boxplus\mu_{2})$, i.e., the right edge
of the bulk of the limiting spectral distribution of
$\widetilde{D}=\widetilde{D}_{1}+\widetilde{D}_{2}$. The former quantity
equals $s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{i}))$ (as
derived in Equation 5.29) and the latter one equals
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\overline{\lambda}(\delta))$ (see Equation
5.16 in the proof of Lemma 5.2). Recall the following two facts:

  1. 1.

$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)=\delta\cdot\psi(\lambda;\delta)$
(as observed in Remark 5.4);

  2. 2.

$\psi(\lambda;\delta)$ is convex in $\lambda$ and increasing for
$\lambda\in[\overline{\lambda}(\delta),\infty)$ (proved in Lemma I.2).

We therefore conclude that
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda^{*}(\delta_{i}))>s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\overline{\lambda}(\delta))$
if $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$. This establishes Item
3 of Lemma 5.3. This condition is more stringent than the previous one
$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta_{i})$. This can be seen by
inspecting the definitions (see, e.g., Equation I.3 in Lemma I.3) of
$\overline{\lambda}(\delta_{i})$ and $\overline{\lambda}(\delta)$:

$\displaystyle\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\delta)-Z}\right)^{2}\right]$
$\displaystyle=\frac{1}{\delta},\quad\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\delta_{i})-Z}\right)^{2}\right]=\frac{1}{\delta_{i}},$
(5.33)

respectively, and realizing that
$\overline{\lambda}(\delta)>\overline{\lambda}(\delta_{i})$ since
$\delta>\delta_{i}$.

We pause and make the following remark regarding the effect of the free
additive convolution on the outliers in the spectra of the addends. Comparing
Equation 5.29 with the limiting value of the right edge of the bulk (cf.
Equation 5.16), we note the following: $\lambda_{1}(\widetilde{D}_{i})$ being
an outlier eigenvalue of $\widetilde{D}_{i}$ does _not_ imply that its image
$\rho_{i}$ under the free additive convolution is also an outlier eigenvalue
of $\widetilde{D}=\widetilde{D}_{1}+\widetilde{D}_{2}$. In fact, it can be
buried strictly inside the bulk, which happens if
$\overline{\lambda}(\delta_{i})<\lambda^{*}(\delta_{i})<\overline{\lambda}(\delta)$.

We then show $\lambda^{*}(\delta_{1})>\lambda^{*}(\delta_{2})$ in Item 1 of
Lemma 5.3. Recall that $\lambda^{*}(\delta_{1})$ and $\lambda^{*}(\delta_{2})$
are the unique solutions to
$\zeta(\lambda^{*}(\delta_{1});\delta_{1})=\varphi(\lambda^{*}(\delta_{1}))$
and
$\zeta(\lambda^{*}(\delta_{2});\delta_{2})=\varphi(\lambda^{*}(\delta_{2}))$,
respectively. Since $\zeta(\cdot;\delta_{1}),\zeta(\cdot;\delta_{2})$ are non-
decreasing and $\varphi(\cdot)$ is strictly decreasing, it suffices to show

$\displaystyle\zeta(\lambda;\delta_{1})<\zeta(\lambda;\delta_{2})$ (5.34)

for any $\lambda>\sup\operatorname{supp}(Z)$. We do so in four steps. (The
following arguments are best understood with Figure 2 in mind.)

  1. 1.

First we claim that
$\overline{\lambda}(\delta_{1})>\overline{\lambda}(\delta_{2})$. This follows
from a similar observation as in Equation 5.33 and the assumption $\alpha>1/2$
(cf. Item (A2)) which implies $\delta_{1}>\delta_{2}$.

  2. 2.

Second we claim that
$\psi(\overline{\lambda}(\delta_{1});\delta_{1})<\psi(\overline{\lambda}(\delta_{2});\delta_{2})$.
Indeed,

$\displaystyle\psi(\overline{\lambda}(\delta_{1});\delta_{1})$
$\displaystyle<\psi(\overline{\lambda}(\delta_{2});\delta_{1})<\psi(\overline{\lambda}(\delta_{2});\delta_{2}).$

The first inequality follows since $\psi(\lambda;\delta_{1})$ is strictly
decreasing for $\lambda\leq\overline{\lambda}(\delta_{1})$ (see Item 2 of
Lemma I.2) and $\overline{\lambda}(\delta_{1})>\overline{\lambda}(\delta_{2})$
as shown in Item 1 above. The second inequality follows since

$\displaystyle\psi(\cdot;\delta_{1})<\psi(\cdot;\delta_{2})$ (5.35)

for any $\lambda>\sup\operatorname{supp}(Z)$ (see the definition of $\psi$ in
Equation 3.2 and also Item 3 of Lemma I.2). Note that in this step we use
$\sup\operatorname{supp}(Z)>0$ in Item (A7). This shows that Equation 5.34
holds for any $\lambda\leq\overline{\lambda}(\delta_{2})$.

  3. 3.

We then claim that Equation 5.34 holds for any
$\lambda\geq\overline{\lambda}(\delta_{1})$. This is because, in this regime,
we have

$\displaystyle\zeta(\lambda;\delta_{1})=\psi(\lambda;\delta_{1})<\psi(\lambda;\delta_{2})=\zeta(\lambda;\delta_{2})$

using the definition of $\zeta(\cdot;\delta_{i})$ (cf. Equation 3.4) and
Equation 5.35.

  4. 4.

Finally, it remains to verify that Equation 5.34 holds for
$\overline{\lambda}(\delta_{2})\leq\lambda\leq\overline{\lambda}(\delta_{1})$.
Indeed, we have

$\displaystyle\zeta(\lambda;\delta_{1})$
$\displaystyle=\psi(\overline{\lambda}(\delta_{1});\delta_{1})<\psi(\overline{\lambda}(\delta_{2});\delta_{2})<\psi(\lambda;\delta_{2}).$

The equality is by definition of $\zeta(\cdot;\delta_{1})$. The first
inequality is by Item 2 above. The second inequality follows since
$\psi(\cdot;\delta_{2})$ is strictly increasing for
$\lambda\geq\overline{\lambda}(\delta_{2})$ (see Item 2 of Lemma I.2).

Combining Items 1, 2, 3 and 4 above then proves Equation 5.34 which implies
Item 1 of Lemma 5.3.

Since $\overline{\lambda}(\delta)$ is the (unique) critical point of
$s_{\mu_{1}\boxplus\mu_{2}}^{-1}(-1/\lambda)$ which is increasing for
$\lambda\geq\overline{\lambda}(\delta)$, Item 2 of Lemma 5.3 then follows.

The proof of the whole lemma is complete. ∎

##  6 Joint distribution via Generalized Approximate Message Passing

The limiting joint distribution in Theorem 3.1 is obtained via a _generalized
approximate message passing_ (GAMP) algorithm whose iterates converge to the
top two eigenvectors of $D=A^{\top}TA$. Within this section, we adopt the
following rescaling for the convenience of applying the GAMP machinery:

$\displaystyle\overline{A}$
$\displaystyle\coloneqq\frac{1}{\sqrt{d}}\,A,\quad\overline{x}_{1}^{*}\coloneqq\sqrt{d}\,x_{1}^{*},\quad\overline{x}_{2}^{*}\coloneqq\sqrt{d}\,x_{2}^{*},\quad\overline{D}\coloneqq\overline{A}^{\top}\,T\,\overline{A}=\frac{n}{d}A^{\top}TA.$
(6.1)

According to Items (A1) and (A4), we have
$\overline{A}_{i,j}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,1/d)$ and
$\overline{x}_{1}^{*},\overline{x}_{2}^{*}\overset{\mathrm{i.i.d.}}{\sim}\operatorname{Unif}(\sqrt{d}\,{\mathbb{S}}^{d-1})$.
Let the $i$-th row of $\overline{A}$ be denoted by
$\overline{a}_{i}^{\top}\in{\mathbb{R}}^{d}$. Then, we have

$\displaystyle y_{i}$ $\displaystyle=q\left(\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle,\varepsilon_{i}\right)=q\left(\left\langle\overline{a}_{i},\eta_{i}\overline{x}_{1}^{*}+(1-\eta_{i})\overline{x}_{2}^{*}\right\rangle,\varepsilon_{i}\right).$

Therefore, $y\in{\mathbb{R}}^{n}$ and related quantities such as
$T\in{\mathbb{R}}^{n\times n}$ (defined in Equation 2.5) do not have to be
rescaled. The overlaps are invariant under rescaling of $D$. Furthermore,
since $n/d\to\delta$, the limiting eigenvalues of $\overline{D}$ are equal to
those of $D$ multiplied by $\delta$ in view of Equation 6.1. The results to be
proved in this section will therefore be consistent with those stated in
Theorem 3.1.

We first extend the GAMP algorithm for the non-mixed GLM [Ran11] and its
associated state evolution analysis to the mixed GLM model. The GAMP algorithm
is defined in terms of a sequence of Lipschitz functions
$g_{t}:\mathbb{R}^{2}\to\mathbb{R}$ and $f_{t+1}:\mathbb{R}^{3}\to\mathbb{R}$,
for $t\geq 0$. For $t\geq 0$, the algorithm iteratively computes
$u^{t},{\tilde{u}}^{t}\in{\mathbb{R}}^{n}$ and
$v^{t+1},{\tilde{v}}^{t+1}\in{\mathbb{R}}^{d}$ as follows:

$\displaystyle\begin{split}u^{t}&=\frac{1}{\sqrt{\delta}}\overline{A}{\tilde{v}}^{t}-{\sf
b}_{t}{\tilde{u}}^{t-1},\quad{\tilde{u}}^{t}=g_{t}(u^{t};y),\\\
v^{t+1}&=\frac{1}{\sqrt{\delta}}\overline{A}^{\top}{\tilde{u}}^{t}-{\sf
c}_{t}{\tilde{v}}^{t},\quad{\tilde{v}}^{t+1}=f_{t+1}(v^{t+1};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*}).\end{split}$
(6.2)

The iteration is initialized with a given ${\tilde{v}}^{0}\in{\mathbb{R}}^{d}$
and ${\tilde{u}}^{-1}=0_{n}$. The functions $f_{t}$ and $g_{t}$ are applied
component-wise, i.e.,
$f_{t}(v^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*})=(f_{t}(v^{t}_{1};\,\overline{x}_{1,1}^{*},\overline{x}_{2,1}^{*})$,
$\ldots,f_{t}(v^{t}_{d};\,\overline{x}_{1,d}^{*},\overline{x}_{2,d}^{*}))$ and
$g_{t}(u^{t};y)=(g_{t}(u^{t}_{1};y_{1}),\ldots,g_{t}(u^{t}_{n};y_{n}))$. The
scalars ${\sf b}_{t},{\sf c}_{t}$ are defined as

${\sf
b}_{t}=\frac{1}{n}\sum_{i=1}^{d}f_{t}^{\prime}(v_{i}^{t};\,\overline{x}_{1,i}^{*},\overline{x}_{2,i}^{*}),\qquad{\sf
c}_{t}=\frac{1}{n}\sum_{i=1}^{n}g_{t}^{\prime}(u_{i}^{t};y_{i}),$ (6.3)

where $f_{t}^{\prime}$ and $g_{t}^{\prime}$ each denote the derivative with
respect to the first argument.

An important feature of the GAMP algorithm is that as $d\to\infty$, the
empirical distributions of the iterates $u^{t}$ and $v^{t+1}$ converge to the
laws of well-defined scalar random variables $U_{t}$ and $V_{t+1}$,
respectively. Specifically, for $t\geq 0$, let

$U_{t}\coloneqq\mu_{1,t}G_{1}+\mu_{2,t}G_{2}+W_{U,t},\quad
V_{t+1}\coloneqq\chi_{1,t+1}X_{1}+\chi_{2,t+1}X_{2}+W_{V,t+1},$ (6.4)

where
$(G_{1},G_{2},W_{U,t})\sim{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,\sigma_{U,t}^{2})$,
and
$(X_{1},X_{2},W_{V,t+1})\sim{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,\sigma_{V,t+1}^{2})$.
We remark that the independence of $X_{1}$ and $X_{2}$ (and analogously of
$G_{1},G_{2}$) follows from the independence of the two signals
$\overline{x}_{1}^{*},\overline{x}_{2}^{*}$. In addition,
$X_{1},X_{2}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,1)$ as the signals
$\overline{x}_{1}^{*},\overline{x}_{2}^{*}$ are uniformly distributed on the
sphere and, hence, their empirical distribution converges to a standard
Gaussian. The deterministic coefficients
$(\mu_{1,t},\mu_{2,t},\sigma_{U,t},\chi_{1,t+1},\chi_{2,t+1},\sigma_{V,t+1})$
are computed using the following state evolution recursion:

$\displaystyle\mu_{1,t}=\frac{1}{\sqrt{\delta}}{\mathbb{E}}[X_{1}f_{t}(V_{t};\,X_{1},X_{2})],\quad\mu_{2,t}=\frac{1}{\sqrt{\delta}}{\mathbb{E}}[X_{2}f_{t}(V_{t};\,X_{1},X_{2})],$
(6.5)
$\displaystyle\sigma_{U,t}^{2}=\frac{1}{\delta}{\mathbb{E}}[f_{t}(V_{t};\,X_{1},X_{2})^{2}]-\mu_{1,t}^{2}-\mu_{2,t}^{2}\,,$