# Sharp-SSL: Selective high-dimensional axis-aligned random projections for
semi-supervised learning

Tengyao Wang∗, Edgar Dobriban†, Milana Gataric‡
and Richard J. Samworth‡
∗Department of Statistics, London School of Economics
†Department of Statistics and Data Science, University of Pennsylvania
‡Statistical Laboratory, University of Cambridge

###### Abstract

We propose a new method for high-dimensional semi-supervised learning problems
based on the careful aggregation of the results of a low-dimensional procedure
applied to many axis-aligned random projections of the data. Our primary goal
is to identify important variables for distinguishing between the classes;
existing low-dimensional methods can then be applied for final class
assignment. Motivated by a generalized Rayleigh quotient, we score projections
according to the traces of the estimated whitened between-class covariance
matrices on the projected data. This enables us to assign an importance weight
to each variable for a given projection, and to select our signal variables by
aggregating these weights over high-scoring projections. Our theory shows that
the resulting Sharp-SSL algorithm is able to recover the signal coordinates
with high probability when we aggregate over sufficiently many random
projections and when the base procedure estimates the whitened between-class
covariance matrix sufficiently well. The Gaussian EM algorithm is a natural
choice as a base procedure, and we provide a new analysis of its performance
in semi-supervised settings that controls the parameter estimation error in
terms of the proportion of labeled data in the sample. Numerical results on
both simulated data and a real colon tumor dataset support the excellent
empirical performance of the method.

##  1 Introduction

Semi-supervised learning, where we attempt to assign observations to one of
finitely many groups based on partially-labeled training data, represents a
core modern statistical challenge. It is sufficiently general to incorporate,
at either extreme, the unsupervised case of no labeled training data
(clustering) and the supervised setting of fully-labeled training data
(classification). Such tasks abound in many application areas, including
genomics (e.g., Eisen et al., 1998), image processing (Jain and Flynn, 1996;
Cheplygina, de Bruijne and Pluim, 2019), natural language processing (Liang,
2005; Turian, Ratinov and Bengio, 2010) and anomaly detection (Akcay, Atapour-
Abarghouei and Breckon, 2019; Wang et al., 2019). Entry points to the
literature on semi-supervised learning include Zhu (2005), Zhu and Goldberg
(2009), Chapelle, Schölkopf and Zien (2006) and Van Engelen and Hoos (2020).
For introductions to clustering, see Xu and Wunsch (2005), Kaufman and
Rousseeuw (2009) and Xu and Tian (2015), and for classification, see Devroye,
Györfi and Lugosi (2013) and Hastie, Tibshirani and Friedman (2009).

A common feature of contemporary semi-supervised learning problems is high-
dimensionality, since we may record many covariates having a possible
association with the labels corresponding to different observations. This
represents a significant challenge, as can be seen by considering a simple
two-class problem with more covariates than observations. For any given
assignment of class labels, if no subset of $n_{0}$ observations lies in an
$(n_{0}-2)$-dimensional affine space, then we can find hyperplanes with
orthogonal normal vectors, each of which achieves zero training error (in
other words, they perfectly separate the classes). Nevertheless, even in the
simple setting where the true Bayes decision boundary is linear, many such
hyperplanes may be little better than a random guess on test data.

An appealing approach to tackling high-dimensionality is via random
projections into lower-dimensional spaces. Such projections may almost
preserve the pairwise distances between observations, as seen from the
Johnson–Lindenstrauss lemma (Johnson and Lindenstrauss, 1984; Dasgupta and
Gupta, 2003). Moreover, in cases where we have reason to believe that only a
relatively small proportion of the variables recorded are relevant for the
learning task, we can choose our random projections to be axis-aligned in
order to preserve this structure. A third benefit is the possibility of
aggregating results over multiple random projections, though this must be done
with care so as to avoid noise accumulation. These attractions have meant that
random projections have now been employed in many high-dimensional statistical
problems, including precision matrix estimation (Marzetta, Tucci and Simon,
2011), two-sample mean testing (Lopes, Jacob and Wainwright, 2011),
classification (Durrant and Kabán, 2015; Cannings and Samworth, 2017),
(sparse) principal component analysis (Yang et al., 2021; Gataric, Wang and
Samworth, 2020), linear regression (Thanei, Heinze and Meinshausen, 2017;
Slawski, 2018; Dobriban and Liu, 2019; Ahfock, Astle and Richardson, 2021),
clustering (Dasgupta, 1999; Fern and Brodley, 2003; Han and Boutin, 2015;
Yellamraju and Boutin, 2018; Anderlucci, Fortunato and Montanari, 2022) and
dimensionality reduction (Bingham and Mannila, 2001; Reeve, Kabán and
Bootkrajang, 2022). See Cannings (2021) for a review of recent developments in
the area.

In this paper, we propose a new method, called Sharp-SSL (short for Selective
high-dimensional, axis-aligned random projections for Semi-Supervised
Learning). Our primary goal is to identify a small subset of variables that
are particularly helpful for label assignment; existing low-dimensional
methods can then be used to complete the learning task. To this end, we
generate a large number of axis-aligned random projections, and apply a base
learning procedure such as a semi-supervised version of the Gaussian
Expectation–Maximization (EM) algorithm to our projected data. Motivated by
the notion of a generalized Rayleigh quotient (see (2) below for a formal
definition), and to avoid the noise accumulation issue mentioned above, we
score the projections by computing the trace of the corresponding estimated
whitened between-class covariance matrices. This enables us to assign an
importance weight to each variable for a given projection, and we select our
signal variables by aggregating these importance weights over the high-scoring
projections. See Section 2 for a more detailed description of our methodology.

Section 3 is devoted to a theoretical analysis of our Sharp-SSL algorithm. We
first show in Theorem 2 that provided the low-dimensional base learning
procedure satisfies a guarantee on the proximity of the estimated whitened
between-class covariance matrix to its population analogue, the corresponding
high-dimensional semi-supervised learning algorithm can recover the signal
coordinates with high probability when we aggregate over sufficiently many
random projections. It turns out that both Linear Discriminant Analysis and an
EM algorithm are examples of low-dimensional learning procedures that satisfy
this proximity guarantee, as we prove in Theorems 3 and 6 respectively. The
latter is particularly challenging, and one of the main novel contributions of
our analysis is to provide a guarantee on the performance of a $d$-dimensional
Gaussian EM algorithm in a semi-supervised setting. In particular, we control
the parameter estimation error in terms of the proportion of labeled data in
the sample, showing that with a sample size of $n$ it smoothly interpolates
between the $(d/n)^{1/4}$ rate for unsupervised learning and the $(d/n)^{1/2}$
rate for fully-labeled data, up to logarithmic factors. An advantage of the
modular approach to our analysis is that it illustrates the way in which the
Sharp-SSL algorithm can be combined with different base learning algorithms to
adapt to different problem settings and reflect the preferences of the
practitioner.

In Section 4, we study the numerical performance of the Sharp-SSL algorithm.
Our first goal, in Section 4.1, is to study the effect of the choices of input
parameters to our method, which allows us to recommend sensible default
choices for application in our subsequent comparisons. Section 4.2 presents
the results of a simulation study involving the Sharp-SSL method, as well as
five alternative approaches, on high-dimensional clustering tasks (since not
all of the competing methods are able to leverage partial label information).
We find that the Sharp-SSL algorithm is able to attain a misclustering rate
very close to that of the optimal Bayes classifier, even with only around 50
observations per cluster, in settings where these alternative techniques may
perform poorly. In Section 4.3, we investigate the extent to which the
different versions of the Sharp-SSL method are able to leverage partial label
information. The results here are consistent with the phase transition
phenomenon articulated by our theory. Finally, in Section 4.4, we apply the
Sharp-SSL algorithm, as well as the other methods from our simulation study,
on a colon tumor dataset, where we withhold the true labels from the
algorithms in order to assess performance. Our analysis supports the ability
of the Sharp-SSL algorithm to identify signal coordinates (genes) that are
useful for identifying patients with and without tumors.

In the broader literature on high-dimensional learning problems, a large
number of methods have been developed to leverage sparse low-dimensional
structures for both clustering (Witten and Tibshirani, 2010; Azizyan, Singh
and Wasserman, 2013; Wasserman, Azizyan and Singh, 2014; Azizyan, Singh and
Wasserman, 2015; Jin and Wang, 2016; Verzelen and Arias-Castro, 2017; Löffler,
Wein and Bandeira, 2022; Löffler, Zhang and Zhou, 2021) and classification
(Cai and Liu, 2011; Witten and Tibshirani, 2011; Mai, Zou and Yuan, 2012; Cai
and Zhang, 2019). These methods are not designed for partially-labeled (semi-
supervised) settings. Another common approach is to project the data into the
span of the top few principal components, and run a standard low-dimensional
method such as $k$-means clustering or the EM algorithm (Butler et al., 2018).
This approach can fail if the directions of largest variation in the data are
not aligned with the directions separating the clusters. Finally, recent
developments in other aspects of semi-supervised learning include self-
training (Oymak and Gulcu, 2020), mean estimation (Zhang, Brown and Cai,
2019), choice of $k$ in $k$-nearest neighbour classification (Cannings,
Berrett and Samworth, 2020) and linear regression (Chakrabortty and Cai,
2018).

Proofs of all of our results are provided in Section 5, and we conclude this
introduction with some notation used throughout the paper. We write
$\mathbb{S}^{d\times d}$ for the set of $d$-dimensional symmetric matrices,
write $\mathbb{S}_{+}^{d\times d}$ for the subset that are invertible, and
write $\mathbb{S}^{d\times d}_{K-1}$ the subset of matrices in
$\mathbb{S}^{d\times d}$ of rank at most $K-1$. We write $\mathbb{R}^{d\times
d}$ for the set of $d$-dimensional matrices. For $p\geq d$, let
$\mathbb{O}^{p\times d}$ denote the set of $p\times d$ matrices with
orthonormal columns. The Euclidean norm is denoted by $\|\cdot\|$, and the
operator norm of a matrix is denoted by $\|\cdot\|_{\mathrm{op}}$, so that
$\|A\|_{\mathrm{op}}:=\sup_{\\{x:\|x\|=1\\}}\|Ax\|$. Given two sequences
$(a_{n})$ and $(b_{n})$, we write $a_{n}\lesssim b_{n}$ when there exists a
universal constant $C>0$ such that $a_{n}\leq Cb_{n}$, and, given an
additional problem parameter $R$, we write $a_{n}\lesssim_{R}b_{n}$ when there
exists $C>0$, depending only on $R$, such that $a_{n}\leq Cb_{n}$.

For any set $S\subseteq\mathbb{R}^{d}$ and $d\leq|S|$, we write
$\binom{S}{d}:=\\{A\subseteq S:|A|=d\\}$. If $S\subseteq\mathbb{R}^{d}$, we
define $\operatorname*{sargmax}S$ to be the smallest element in the
$\operatorname*{argmax}$ in the lexicographic order. For a positive integer
$k$, we define $[k]:=\\{1,\ldots,k\\}$. For a vector
$v=(v_{1},\ldots,v_{k})^{\top}\in\mathbb{R}^{k}$, and $j\in[k]$, we define
$v_{-j}=(v_{1},\ldots,v_{j-1},v_{j+1},\ldots,v_{k})^{\top}\in\mathbb{R}^{k-1}$.

##  2 The Sharp-SSL algorithm

In this section, we describe in detail the Sharp-SSL algorithm for $K$-class
semi-supervised learning, with $K\geq 2$. We aim to provide a unified
treatment of clustering, semi-supervised learning and classification. To this
end, we assume that for $i\in[n]$, the observation $x_{i}\in\mathbb{R}^{p}$
has a true label $y_{i}^{*}\in[K]$, but it may be the case that we do not
observe $y_{i}^{*}$. Instead, we assume that our observed label $y_{i}$ takes
values in $[K]\cup\\{0\\}$, where $y_{i}:=y_{i}^{*}$ when the true class label
is observed, and $y_{i}:=0$ otherwise. Thus, our data can be regarded as
$(x_{1},y_{1}),\ldots,(x_{n},y_{n})\in\mathbb{R}^{p}\times([K]\cup\\{0\\})$,
and our goal is to construct a _data-dependent classifier_ 111It is convenient
to use the term ‘classifier’ here, even though some or all of the labels may
be unobserved., i.e. a Borel measurable function
$C:\mathbb{R}^{p}\times\bigl{(}\mathbb{R}^{p}\times([K]\cup\\{0\\})\bigr{)}^{n}\rightarrow[K]$,
with the interpretation that
$C\bigl{(}x;(x_{1},y_{1}),\ldots,(x_{n},y_{n})\bigr{)}$ is the predicted class
of $x\in\mathbb{R}^{p}$.

To motivate our Sharp-SSL algorithm, it is instructive first to consider a
canonical Gaussian classification problem, where our data can be regarded as
$n$ independent realizations of a pair $(X,Y)$ taking values in
$\mathbb{R}^{p}\times[K]$, with prior probability $\pi_{k}:=\mathbb{P}(Y=k)$
for the $k$th class and $X\mid
Y=k\sim\mathcal{N}_{p}(\nu_{k},\Sigma_{\mathrm{w}})$, for class means
$\nu_{1},\ldots,\nu_{K}\in\mathbb{R}^{p}$ and _within-class covariance matrix_
$\Sigma_{\mathrm{w}}\in\mathbb{S}_{+}^{p\times p}$. Let
$\nu:=\sum_{k=1}^{K}\pi_{k}\nu_{k}\in\mathbb{R}^{p}$ denote the grand
population mean, let

$\Sigma_{\mathrm{b}}:=\sum_{k=1}^{K}\pi_{k}(\nu_{k}-\nu)(\nu_{k}-\nu)^{\top}\in\mathbb{S}_{K-1}^{p\times
p}$ (1)

denote the _between-class covariance matrix_ , and consider
$D\in\mathbb{O}^{p\times(K-1)}$ with a column space spanned by
$\Sigma_{\mathrm{w}}^{-1}(\nu_{1}-\nu),\ldots,\Sigma_{\mathrm{w}}^{-1}(\nu_{K}-\nu)$.
Observe that for $k\neq\ell$,

$\log\biggl{\\{}\frac{\mathbb{P}(Y=k\mid X=x)}{\mathbb{P}(Y=\ell\mid
X=x)}\biggr{\\}}=\log\Bigl{(}\frac{\pi_{k}}{\pi_{\ell}}\Bigr{)}-\frac{1}{2}(\nu_{k}+\nu_{\ell})^{\top}\Sigma_{\mathrm{w}}^{-1}(\nu_{k}-\nu_{\ell})+x^{\top}\Sigma_{\mathrm{w}}^{-1}(\nu_{k}-\nu_{\ell}),$

from which we deduce that this likelihood ratio, and hence the Bayes
classifier $x\mapsto\operatorname*{argmax}_{k\in[K]}\mathbb{P}(Y=k\mid X=x)$,
only depends on $x$ through $D^{\top}x$. Thus, for the purposes of
classification, no signal would be lost (and the noise would be reduced) if
$X$ were replaced with $D^{\top}X$.

In high-dimensional settings with $p\gg n$, the matrix
$\Sigma_{\mathrm{w}}^{-1}$ is not consistently estimable in general, but we
can nevertheless make progress if the vectors
$\Sigma_{\mathrm{w}}^{-1}(\nu_{1}-\nu),\ldots,\Sigma_{\mathrm{w}}^{-1}(\nu_{K}-\nu)$
are sparse. In other words, writing $S_{0}$ for the union of the set of
coordinates for which these vectors are non-zero, we suppose that $|S_{0}|\ll
p$; this is a very common assumption in high-dimensional LDA (e.g. Cai and
Liu, 2011; Witten and Tibshirani, 2011; Mai, Zou and Yuan, 2012; Cai and
Zhang, 2019).

In such a setting, the column space of $D$ has a sparse basis, so it is
natural to consider projecting the data onto a small subset of its
coordinates. For $d\in[p]$, define the set of axis-aligned projection matrices
$\mathcal{P}_{d}:=\bigl{\\{}P\in\\{0,1\\}^{d\times
p}:PP^{\top}=I_{d}\bigr{\\}}$, i.e. the set of binary $d\times p$ matrices
with orthonormal rows. By the argument above, if $d\geq|S_{0}|$ then there
exists $P^{*}\in\mathcal{P}_{d}$ such that the error of the Bayes classifier
is unchanged by projecting the data along $P^{*}$. In practice, it would
typically be computationally too expensive to enumerate through all
$p(p-1)\cdots(p-d+1)$ axis-aligned projections. Instead, we consider a
randomly chosen subset of projections within $\mathcal{P}_{d}$. An axis-
aligned projection chosen uniformly at random is unlikely to capture all the
signal coordinates $S_{0}$, but by aggregating over a carefully-chosen subset
of these random projections, we can nevertheless recover the set of signal
coordinates under suitable conditions; see Theorem 2 below. To describe our
method for choosing good projections, for $V\in\mathbb{O}^{p\times d}$, we
define the _generalized Rayleigh quotient_ along $V$ by

$J(V;\Sigma_{\mathrm{b}},\Sigma_{\mathrm{w}}):=\mathrm{tr}\\{(V^{\top}\Sigma_{\mathrm{w}}V)^{-1}(V^{\top}\Sigma_{\mathrm{b}}V)\\}.$
(2)

Proposition 1 below motivates seeking to choose projections to maximize the
generalized Rayleigh quotient by showing that the column span of any maximizer
$J(V;\Sigma_{\mathrm{b}},\Sigma_{\mathrm{w}})$ over $V\in\mathbb{O}^{p\times
d}$ must contain the column space of $D$.

######  Proposition 1.

Let $K\geq 2$ and $d\geq K-1$. Assume that the convex hull of
$\nu_{1},\ldots,\nu_{K}$ is $(K-1)$-dimensional, and let
$V^{*}\in\operatorname*{argmax}_{V\in\mathbb{O}^{p\times
d}}J(V;\Sigma_{\mathrm{b}},\Sigma_{\mathrm{w}})$. Then the column space of
$V^{*}$ contains the eigenspace corresponding to the $K-1$ non-zero
eigenvalues222Even though $\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$ is not
guaranteed to be symmetric, it is similar (i.e. conjugate) to the symmetric
matrix
$\Sigma_{\mathrm{w}}^{-1/2}\Sigma_{\mathrm{b}}\Sigma_{\mathrm{w}}^{-1/2}$, so
has real eigenvalues and eigenvectors. of
$\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$, which is equal to the space
spanned by $\bigl{(}\Sigma_{\mathrm{w}}^{-1}(\nu_{k}-\nu):k\in[K]\bigr{)}$.

Based on Proposition 1, a natural conceptual approach to maximizing the
generalized Rayleigh quotient is to compute the leading $(K-1)$-dimensional
eigenspace of $\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$. This strategy,
however, runs into difficulties when we replace these population quantities
with their sample versions in the setting of the opening paragraph of this
section. More precisely, writing
$n_{k}:=\sum_{i=1}^{n}\mathbbm{1}_{\\{y_{i}=k\\}}$ for $k\in[K]$, as well as

$\displaystyle\tilde{\Sigma}_{\mathrm{w}}$
$\displaystyle:=\frac{1}{n}\sum_{k=1}^{K}\sum_{i=1}^{n}(x_{i}-\hat{\nu}_{k})(x_{i}-\hat{\nu}_{k})^{\top}\mathbbm{1}_{\\{y_{i}=k\\}}\in\mathbb{R}^{p\times
p}$ $\displaystyle\tilde{\Sigma}_{\mathrm{b}}$
$\displaystyle:=\sum_{k=1}^{K}\frac{n_{k}}{n}(\hat{\nu}_{k}-\hat{\nu})(\hat{\nu}_{k}-\hat{\nu})^{\top}\in\mathbb{R}^{p\times
p},$

for the sample versions of the within-class and between-class covariance
matrices respectively, the matrix $\tilde{\Sigma}_{\mathrm{w}}$ is not
invertible whenever $p>n$. Fortunately, though, this issue can be resolved by
working with the projected data, as long as we choose $d\leq n-K$: the
projected data $\\{PX_{i}:i\in[n]\\}$ has within-class covariance matrix
$P\Sigma_{\mathrm{w}}P^{\top}\in\mathbb{R}^{d\times d}$ and between-class
covariance matrix $P\Sigma_{\mathrm{b}}P^{\top}\in\mathbb{R}^{d\times d}$, so
with probability one, the sample version
$P\tilde{\Sigma}_{\mathrm{w}}P^{\top}$ is invertible.

Returning to the general setting of the opening paragraph of this section,
then, we seek projections $P$ with large
$J(P^{\top};\tilde{\Sigma}_{\mathrm{b}},\tilde{\Sigma}_{\mathrm{w}})=\mathrm{tr}\bigl{\\{}(P\tilde{\Sigma}_{\mathrm{w}}P^{\top})^{-1}(P\tilde{\Sigma}_{\mathrm{b}}P^{\top})\bigr{\\}}$.
To this end, for fixed $A,B\in\mathbb{N}$, we sample a set of projections
$\\{P^{a,b}:a\in[A],b\in[B]\\}$ uniformly at random from $\mathcal{P}_{d}$.
For each $a$ and $b$, we apply a low-dimensional base algorithm
$\psi:\bigl{(}\mathbb{R}^{d}\times([K]\cup\\{0\\})\bigr{)}^{n}\rightarrow\mathbb{S}_{+}^{d\times
d}$ to the projected data $(P^{a,b}x_{1},y_{1}),\ldots,(P^{a,b}x_{n},y_{n})$
to obtain an estimator $\hat{Q}^{a,b}$ of
$(P^{a,b}\Sigma_{\mathrm{w}}P^{a,b,\top})^{-1}(P^{a,b}\Sigma_{\mathrm{b}}P^{a,b,\top})$,
the _whitened between-class covariance matrix_ of the projected data. We
assume throughout for convenience that $\psi$ is _permutation equivariant_ in
the sense that $\psi\bigl{(}(\Pi z_{1},y_{1}),\ldots,(\Pi
z_{n},y_{n})\bigr{)}=\Pi\psi\bigl{(}(z_{1},y_{1}),\ldots,(z_{n},y_{n})\bigr{)}\Pi^{\top}$
for every permutation matrix $\Pi\in\mathbb{R}^{d\times d}$. One choice for
the base algorithm is to set
$\hat{Q}^{a,b}=(\hat{\Sigma}^{a,b}_{\mathrm{w}})^{-1}\hat{\Sigma}^{a,b}_{\mathrm{b}}$,
where $\hat{\Sigma}^{a,b}_{\mathrm{w}}$ and $\hat{\Sigma}^{a,b}_{\mathrm{b}}$
are estimated projected within- and between-class (or cluster) covariance
matrices.

To select projections, for each $a\in[A]$, we define

$b^{*}(a):=\operatorname*{sargmax}_{b\in[B]}\mathrm{tr}(\hat{Q}^{a,b}),$

and select $P^{a,b^{*}(a)}$. The main rationale for dividing the projections
into $A$ groups and selecting one within each group—as opposed to selecting
the $A$ projections with the largest values of $\mathrm{tr}(\hat{Q}^{a,b})$—is
that, conditional on the original data, the selected projections are
independent and identically distributed. This facilitates our theoretical
analysis by enabling the application of concentration inequalities in the
proof of Theorem 2.

The diagonal entries of $\\{\hat{Q}^{a,b^{*}(a)}:a\in[A]\\}$ measure the
importance of the projected variables for the semi-supervised learning task.
These can then be converted into importance scores for the original variables
by ‘back-projecting’ into the higher-dimensional space, i.e. by forming the
vector $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{p})^{\top}\in\mathbb{R}^{p}$
given by

$\hat{w}_{j}:=\frac{1}{A}\sum_{a=1}^{A}\left[P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\right]_{j,j},\quad
j\in[p].$

Finally, we rank the variables by their importance scores, and our estimate
$\hat{S}$ of the set of signal coordinates is given by the largest $\ell$
entries in $\hat{w}$, breaking ties arbitrarily if necessary, where
$\ell\in[p]$ is specified by the practitioner. Pseudocode for the Sharp-SSL
procedure is given in Algorithm 1.

Input: Data $(x_{1},y_{1}),\ldots,(x_{n},y_{n})$ $\in$ $\mathbb{R}^{p}\times$
$([K]\cup\\{0\\})$ (where $y_{i}=0$ denotes a missing label);
Projected dimension $d\in[\min(p,n-K)]$, number of selected signal coordinates
$\ell\in[p]$;
Number $A\in\mathbb{N}$ of groups of projections, number $B\in\mathbb{N}$ of
projections in each group;
Permutation equivariant base algorithm
$\psi:\bigl{(}\mathbb{R}^{d}\times([K]\cup\\{0\\})\bigr{)}^{n}$
$\to\mathbb{S}_{+}^{d\times d}$.

Generate axis-aligned random projections $\\{P^{a,b}:a\in[A],b\in[B]\\}$
independently and uniformly from $\mathcal{P}_{d}$.

for  _$a\in[A]$_ do

for  _$b\in[B]$_ do

Let
$\hat{Q}^{a,b}:=\psi\bigl{(}(P^{a,b}x_{1},y_{1}),\ldots,(P^{a,b}x_{n},y_{n})\bigr{)}$.

end for

Set $b^{*}(a):=\operatorname*{sargmax}_{b\in[B]}\mathrm{tr}(\hat{Q}^{a,b})$.

end for

Let $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{p})^{\top}$, where
$\hat{w}_{j}:=\frac{1}{A}\sum_{a=1}^{A}[P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}]_{j,j}$.

Output: $\hat{S}\subseteq[p]$, defined as the index set of the $\ell$ largest
components of $\hat{w}$, breaking ties randomly.

Algorithm 1 Sharp-SSL: Clustering and semi-supervised learning via ensembles
of axis-aligned random projections.

After applying Algorithm 1 to obtain an estimated set $\hat{S}$ of signal
variables, we can then apply any existing semi-supervised learning method for
low-dimensional data with input $(P_{\hat{S}}x_{i,},y_{i})_{i\in[n]}$, where
$P_{\hat{S}}$ denotes the projection onto the coordinates in $\hat{S}$.

###  2.1 Base learning methods

Algorithm 1 relies on a base learning method for low-dimensional data to
estimate the projected whitened between-class covariance matrix from the
projected data. When all or almost all of the input data are labeled, we can
use the procedure outlined in Algorithm 2, which ignores any unlabeled data,
for this purpose. On the other hand, when we have a substantial amount of
unlabeled data, Algorithm 2 may be inaccurate. In such circumstances, it may
be preferable to use Algorithm 3, which runs an _Expectation–Maximization_
(EM) procedure to predict the unobserved labels and subsequently estimate the
whitened between-class covariance matrix. More precisely, from $M$ random
initializations of the cluster means and the within-class covariance matrix,
Algorithm 3 uses the EM algorithm to update these quantities, and thereby
compute the whitened between-cluster sample covariance matrix estimators
$\bigl{\\{}\hat{Q}^{[m]}=(\hat{\Sigma}_{\mathrm{w}}^{[m]})^{-1}(\hat{\Sigma}_{\mathrm{b}}^{[m]}):m\in[M]\bigr{\\}}$.
We select $\hat{m}\in[M]$ such that $\hat{Q}^{[\hat{m}]}$ is in best agreement
with results from the other EM runs; this is made precise in (6).

The algorithm also allows the practitioner to incorporate prior knowledge
about the true cluster means and within-cluster covariance matrices, both
through optimizing over a restricted constraint set $\mathcal{C}$ in the M
step of the EM algorithm, and through the choice of a distribution supported
on $\mathcal{C}$ for the initialization of these quantities. An alternative to
the EM algorithm for unsupervised learning would be to apply $k$-means
clustering as a base procedure. Previous studies have suggested that these
approaches have comparable empirical performance (e.g., de Souto et al., 2008;
Rodriguez et al., 2019, and references therein), but the EM algorithm is more
amenable to theoretical analysis in our setting.

Input:
$(z_{1},y_{1}),\ldots,(z_{n},y_{n})\in\mathbb{R}^{d}\times([K]\cup\\{0\\})$,
where the convex hull of $(z_{i})_{i:y_{i}=k}$ is $d$-dimensional for at least
one $k\in[K]$.

Compute $\hat{\mu}:=n^{-1}\sum_{i=1}^{n}z_{i}$

for  _$k\in[K]$_ do

Set $n_{k}:=|\\{i:y_{i}=k\\}|$ and $n^{\prime}:=\sum_{k=1}^{K}n_{k}$

Compute $\hat{\mu}_{k}:=n_{k}^{-1}\sum_{i:y_{i}=k}z_{i}$ (with the convention
that $\hat{\mu}_{k}:=0$ if $n_{k}=0$).

end for

Compute the within-class and between-class covariance matrices as
$\hat{\Sigma}_{\mathrm{w}}:=\frac{1}{n^{\prime}}\sum_{i=1}^{n}(z_{i}-\hat{\mu}_{y_{i}})(z_{i}-\hat{\mu}_{y_{i}})^{\top}\quad\text{and}\quad\hat{\Sigma}_{\mathrm{b}}:=\sum_{k=1}^{K}\frac{n_{k}}{n^{\prime}}(\hat{\mu}_{k}-\hat{\mu})(\hat{\mu}_{k}-\hat{\mu})^{\top}.$
(3)

Output: $\hat{Q}:=\hat{\Sigma}_{\mathrm{w}}^{-1}\hat{\Sigma}_{\mathrm{b}}$.

Algorithm 2 Base learning using only labeled data

Input: Data
$(z_{1},y_{1}),\ldots,(z_{n},y_{n})\in\mathbb{R}^{d}\times([K]\cup\\{0\\})$. A
constraint set
$\mathcal{C}\subseteq(\mathbb{R}^{d})^{K}\times\mathbb{S}_{+}^{d\times d}$ and
a probability distribution $\pi_{\mathcal{C}}$ supported on $\mathcal{C}$.
Number of random initializations $M$. Number of iterations $T$.

for  _$m\in[M]$_ do

Randomly sample
$(\hat{\mu}_{1},\ldots,\hat{\mu}_{K},\hat{\Sigma}_{\mathrm{w}})\sim\pi_{\mathcal{C}}$.

for  _$t\in[T]$_ do

(E step) Compute the soft-label matrix $(L_{i,k})_{i\in[n],k\in[K]}$
$L_{i,k}:=\biggl{(}\frac{e^{-\frac{1}{2}(z_{i}-\hat{\mu}_{k})^{\top}\hat{\Sigma}_{\mathrm{w}}^{-1}(z_{i}-\hat{\mu}_{k})}}{\sum_{\ell=1}^{K}e^{-\frac{1}{2}(z_{i}-\hat{\mu}_{\ell})^{\top}\hat{\Sigma}_{\mathrm{w}}^{-1}(z_{i}-\hat{\mu}_{\ell})}}\biggr{)}\mathbbm{1}_{\\{y_{i}=0\\}}+\mathbbm{1}_{\\{y_{i}=k\\}}.$
(4)

(M step) Update parameter estimates by  $\displaystyle(\hat{\mu}_{1}$
$\displaystyle,\ldots,\hat{\mu}_{K},\hat{\Sigma}_{\mathrm{w}})$
$\displaystyle:=\operatorname*{argmin}_{(\mu_{1},\ldots,\mu_{K},\Sigma)\in\mathcal{C}}\biggl{\\{}\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{K}L_{i,k}(z_{i}-\mu_{k})^{\top}\Sigma^{-1}(z_{i}-\mu_{k})+\log\det\Sigma\biggr{\\}}.$
(5)

end for

Compute $(L_{i,k})_{i\in[n],k\in[K]}$ using the final values of
$(\hat{\mu}_{1},\ldots,\hat{\mu}_{K},\hat{\Sigma}_{\mathrm{w}})$ as in (4).

Compute
$\hat{\mu}_{\mathrm{tot}}:=n^{-1}\sum_{i=1}^{n}\sum_{k=1}^{K}L_{i,k}\hat{\mu}_{k}$
and the between-class covariance matrix
$\hat{\Sigma}_{\mathrm{b}}:=\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{K}L_{i,k}(\hat{\mu}_{k}-\hat{\mu}_{\mathrm{tot}})(\hat{\mu}_{k}-\hat{\mu}_{\mathrm{tot}})^{\top}.$
Set $\hat{Q}^{[m]}:=\hat{\Sigma}_{\mathrm{w}}^{-1}\hat{\Sigma}_{\mathrm{b}}$.

end for

Set
$\hat{m}\in\operatorname*{argmin}_{m\in[M]}\mathrm{median}\bigl{(}\|\hat{Q}^{[m]}-\hat{Q}^{[m^{\prime}]}\|_{\mathrm{op}}:m^{\prime}\in[M]\setminus\\{m\\}\bigr{)}.$
(6)

Output: $\hat{Q}:=\hat{Q}^{[\hat{m}]}$.

Algorithm 3 Base learning using partially labeled data via an EM algorithm

##  3 Theoretical guarantees

###  3.1 Results for the high-level algorithm

In this subsection, we consider independent triples
$(X_{1},Y_{1},Y_{1}^{*}),\ldots,(X_{n},Y_{n},Y_{n}^{*})$ taking values in
$\mathbb{R}^{p}\times([K]\cup\\{0\\})\times[K]$. We recall that $Y_{i}^{*}$
denotes the true label of the $i$th observation, and that $Y_{i}:=Y_{i}^{*}$
if the $i$th label is observed, and $Y_{i}:=0$ otherwise. For $k\in[K]$, let
$\pi_{k}:=\mathbb{P}(Y_{1}^{*}=k)$ and $\nu_{k}^{*}:=\mathbb{E}(X_{1}\mid
Y_{1}^{*}=k)$ denote the prior probability and the cluster mean of the $k$th
cluster respectively, let $\nu^{*}:=\sum_{k=1}^{K}\pi_{k}\nu_{k}^{*}$ denote
the weighted cluster mean and let $\Sigma_{\mathrm{w}}:=\mathrm{Cov}(X_{1}\mid
Y_{1}^{*}=k)$ denote the common within-cluster covariance matrix. With the
between-cluster covariance matrix $\Sigma_{\mathrm{b}}$ from (1), our goal is
to estimate the set of signal coordinates,

$S_{0}:=\bigl{\\{}j\in[p]:(\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}})_{j,j}\neq
0\bigr{\\}},$

and we write $s_{0}:=|S_{0}|$.

Our first main theoretical result shows that if the base algorithm is accurate
on each low-dimensional projection and $A$ is large, then with high
probability, all signal coordinates are selected.

######  Theorem 2.

Define $\gamma_{\min}:=\min_{j\in
S_{0}}(\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}})_{j,j}$ and
$\gamma_{\max}:=\max_{j\in
S_{0}}(\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}})_{j,j}$. Let $\hat{S}$ be
the output of Algorithm 1 with input $K$, $p$,
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, $A$, $B$, $d\geq s_{0}$, $\ell\geq
s_{0}$ and permutation equivariant base procedure $\psi$. Write

$\varepsilon:=\mathbb{P}\biggl{(}\max_{P\in\mathcal{P}_{d}}\bigl{\|}\psi\bigl{(}(PX_{i},Y_{i})_{i\in[n]}\bigr{)}-P\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}P^{\top}\bigr{\|}_{\mathrm{op}}\geq\frac{\gamma_{\min}}{4(K-1)}\biggr{)}.$
(7)

Then

$\mathbb{P}(S_{0}\subseteq\hat{S})\geq 1-\varepsilon-
pe^{-A\gamma_{\min}^{2}/(50p^{2}\gamma_{\max}^{2})}.$

In fact, we can see from the proof of Theorem 2 that the following stronger
conclusion holds: for any realization $(x_{i},y_{i})_{i\in[n]}$ of the data
satisfying

$\max_{P\in\mathcal{P}_{d}}\bigl{\|}\psi\bigl{(}(Px_{i},y_{i})_{i\in[n]}\bigr{)}-P\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}P^{\top}\bigr{\|}_{\mathrm{op}}<\frac{\gamma_{\min}}{4(K-1)},$
(8)

we have
$\mathbb{P}\bigl{(}S_{0}\subseteq\hat{S}\mid(X_{i},Y_{i})_{i\in[n]}=(x_{i},y_{i})_{i\in[n]}\bigr{)}\geq
1-pe^{-A\gamma_{\min}^{2}/(50p^{2}\gamma_{\max}^{2})}$. Note here that, after
conditioning on the data, the probability is taken over the randomness in the
projections. An attraction of Theorem 2 is its generality, and in particular
the fact that we do not impose strong distributional assumptions — we simply
require control of $\varepsilon$ in (7). The price we pay for this generality
is that the probability bound may be loose in particular cases; for example,
the bound holds even with $B=1$, though in practice we would expect it to
improve as $B$ increases.

###  3.2 Theory for base learning using labeled data

In this subsection, we demonstrate how the high-level result in Theorem 2 can
be used to derive performance guarantees for a high-dimensional classification
algorithm that uses the Sharp-SSL procedure in Algorithm 1 in conjunction with
the low-dimensional base method described in Algorithm 2 for estimating the
projected whitened between-class covariance matrix. The following theorem
provides uniform control of the output of Algorithm 2 for all axis-aligned
$d$-dimensional projected datasets.

######  Theorem 3.

Fix $\varepsilon\in(0,1]$, $K\in\\{2,3,\ldots,\\}$, $p,d\in\mathbb{N}$ with
$p\geq d$ and $n\in\mathbb{N}$ with $n\geq Kd+1$. Suppose that
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ are independent and identically
distributed pairs, with $\mathbb{P}(Y_{1}=k)=\pi_{k}$ and $X_{1}\mid
Y_{1}=k\sim\mathcal{N}_{p}(\nu_{k}^{*},\Sigma_{\mathrm{w}})$ for $k\in[K]$,
and let $\psi\bigl{(}(PX_{i},Y_{i})_{i\in[n]}\bigr{)}$ be the output of
Algorithm 2 with input $(X_{i},Y_{i})_{i\in[n]}$, for $P\in\mathcal{P}_{d}$.
Suppose that $\|\nu_{k}^{*}-\nu_{\ell}^{*}\|\leq R_{1}$ for all $k,\ell\in[K]$
and some $R_{1}>0$, and that $\Sigma_{\mathrm{w}}$ is diagonal and well-
conditioned in the sense that
$\max\\{\|\Sigma_{\mathrm{w}}\|_{\mathrm{op}},\|\Sigma_{\mathrm{w}}^{-1}\|_{\mathrm{op}}\\}\leq
R_{2}$ for some $R_{2}\geq 1$. If

$\frac{16R_{2}^{2}K}{n}\leq
1\quad\text{and}\quad\frac{32R_{2}^{2}}{n^{1/2}}\log^{1/2}\biggl{(}\frac{8\cdot
9^{d}\binom{p}{d}}{\varepsilon}\biggr{)}\leq 1,$ (9)

then with probability at least $1-\varepsilon$, we have

$\displaystyle\max_{P\in\mathcal{P}_{d}}\bigl{\|}\psi\bigl{(}(PX_{i},Y_{i})_{i\in[n]}\bigr{)}-(P\Sigma_{\mathrm{w}}P^{\top})^{-1}P\Sigma_{\mathrm{b}}P^{\top}$
$\displaystyle\bigr{\|}_{\mathrm{op}}$
$\displaystyle\lesssim_{R_{1},R_{2}}\frac{K}{n}+\sqrt{\frac{d\log(ep/d)+\log(1/\varepsilon)}{n}}.$

The sample size condition (9) can be restated as $n\gtrsim_{R_{1},R_{2}}d\log
p+\log(1/\varepsilon)+K$, so may be regarded as mild. Regarding $K$ as a
constant, Theorem 3 confirms that the uniform control of Algorithm 2 is at the
parametric rate, up to a logarithmic factor. The following corollary then
follows immediately by combining Theorems 2 and 3.

######  Corollary 4.

Fix $\varepsilon\in(0,1]$. Suppose that the conditions of Theorem 3 hold, and
moreover that $\lambda_{\min}(\Sigma_{\mathrm{b}})\geq 1/R_{3}$ for some
$R_{3}>0$. Then there exist $C_{1},C_{2}>0$, depending only on $R_{1},R_{2}$
and $R_{3}$, such that if

$C_{1}\biggl{(}\frac{K}{n}+\sqrt{\frac{d\log(ep/d)+\log(1/\varepsilon)}{n}}\biggr{)}\leq\frac{1}{K},$

then the output $\hat{S}$ of Algorithm 1 with input $K$, $p$, $d\geq s_{0}$,
$\ell\geq s_{0}$, $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, $A$, $B$, and base
procedure $\psi$ from Algorithm 2 satisfies

$\mathbb{P}(S_{0}\subseteq\hat{S})\geq
1-\varepsilon-p\exp\biggl{(}-\frac{A}{C_{2}p^{2}}\biggr{)}.$

Thus, under the conditions of Corollary 4, the Sharp-SSL algorithm can, with
high probability, select the signal variables in the top $s_{0}$ output
variables, provided that the number $A$ of groups of random projections is
large by comparison with $p^{2}$. In other words, the algorithm reduces the
problem to a low-dimensional one, for which standard learning techniques can
be applied. The guarantees for these methods (e.g. Anderson, 2003, Theorem
6.6.1) can then be combined on the high-probability event of Corollary 4 to
establish theoretical results for the full procedure.

###  3.3 Theory for semi-supervised based learning

When the proportion of labeled data is low, Algorithm 2 may be inaccurate when
used as the base procedure in Algorithm 1. The aim of this subsection,
therefore, is to study the base procedure of Algorithm 3, which is able to
leverage both the labeled and unlabeled data via an EM algorithm to estimate
the whitened between-class covariance matrix for each projected data set. Our
analysis builds on several recent breakthroughs in our understanding of the EM
algorithm. This line of work includes Balakrishnan, Wainwright and Yu (2017),
Daskalakis, Tzamos and Zampetakis (2017), Yan, Yin and Sarkar (2017), Dwivedi
et al. (2020a), Dwivedi et al. (2020b), Davis, Diaz and Wang (2021), Ho et al.
(2020), Ndaoud (2022), Wu and Zhou (2022) and Doss et al. (2023), all of which
focus on the unsupervised case.

For simplicity, we will focus on the setting where independent and identically
distributed $(X_{1},Y_{1}^{*}),\ldots,(X_{n},Y_{n}^{*})$ are generated from a
mixture of two Gaussians with opposite means and identity covariance matrix:

$Y_{i}^{*}\sim\mathrm{Unif}(\\{1,2\\}),\,X_{i}\mid
Y_{i}^{*}\sim\mathcal{N}_{p}\bigl{(}(-1)^{Y_{i}^{*}}\nu^{*},I_{p}\bigr{)},\text{
and }Y_{i}=Y_{i}^{*}\mathbbm{1}_{\\{i\leq n_{\mathrm{L}}\\}}\text{ for all
$i\in[n]$}.$ (10)

We assume that we observe
$(X_{1},Y_{1}),\ldots,(X_{n_{\mathrm{L}}},Y_{n_{\mathrm{L}}}),X_{n_{\mathrm{L}}+1},\ldots,X_{n}$
for some $n_{\mathrm{L}}\in\\{0,\ldots,n\\}$. In other words, we are given
$n_{\mathrm{L}}$ labeled observations and $n_{\mathrm{U}}:=n-n_{\mathrm{L}}$
unlabeled ones. Thus, $n_{\mathrm{L}}=0$ corresponds to the fully unsupervised
case, i.e., clustering, while $n_{\mathrm{L}}=n$ corresponds to the supervised
case, i.e., classification. We define $Y_{i}=Y_{i}^{*}$ for
$i\in[n_{\mathrm{L}}]$, and $Y_{i}=0$ for
$i\in\\{n_{\mathrm{L}}+1,\ldots,n\\}$.

We first study the performance of the EM procedure after the covariates have
been projected into a lower-dimensional space. In other words, for some fixed
$P\in\mathcal{P}_{d}$, define $Z_{i}:=PX_{i}$ for $i\in[n]$ and
$\mu^{*}:=P\nu^{*}\in\mathbb{R}^{d}$, so that $Z_{i}\mid
Y_{i}^{*}\sim\mathcal{N}_{d}((-1)^{Y_{i}^{*}}\mu^{*},I_{d})$. In this setting,
we have a single unknown parameter $\mu^{*}$ to estimate, and this can be
achieved by applying Algorithm 3 to $(Z_{i},Y_{i})_{i\in[n]}$ with $K=2$ and
the constraint set

$\mathcal{C}:=\bigl{\\{}(-\mu,\mu,I_{d}):\mu\in\mathbb{R}^{d}\bigr{\\}}.$ (11)

After initializing the EM algorithm at some fixed
$(-\hat{\mu}^{(0)},\hat{\mu}^{(0)},I_{d})\in\mathcal{C}$, for
$t\in\mathbb{N}$, the $t$th iterate of the EM iteration described in (4) and
(3) is $(-\hat{\mu}^{(t)},\hat{\mu}^{(t)},I_{d})$, where

$\hat{\mu}^{(t)}:=\frac{1}{n}\biggl{\\{}\sum_{i:Y_{i}\neq
0}(-1)^{Y_{i}}Z_{i}+\sum_{i:Y_{i}=0}Z_{i}\tanh\big{\langle}Z_{i},\hat{\mu}^{(t-1)}\big{\rangle}\biggr{\\}};$
(12)

see Lemma 15. Since we allow $n_{\mathrm{L}}=0$, where $\mu$ is only
identifiable up to sign, and since the between-class sample covariance matrix
$\hat{\Sigma}_{\mathrm{b}}$ computed in Algorithm 3 is equal to
$\hat{\Sigma}_{\mathrm{b}}=\hat{\mu}_{1}\hat{\mu}_{1}^{\top}-\hat{\mu}_{\mathrm{tot}}\hat{\mu}_{\mathrm{tot}}^{\top}$,
which is invariant to flipping the signs of $\hat{\mu}_{1}$ and
$\hat{\mu}_{2}$ simultaneously, it is natural to consider the loss function
$L:\mathbb{R}^{d}\times\mathbb{R}^{d}\to[0,\infty)$ given by

$L(\mu,\mu^{\prime}):=\|\mu-\mu^{\prime}\|\wedge\|\mu+\mu^{\prime}\|.$

Proposition 5 below provides a theoretical guarantee for this semi-supervised
EM algorithm. For notational simplicity, we define $\gamma:=n_{\mathrm{L}}/n$,
$\omega_{0}:=\sqrt{\\{d\log n+\log(1/\delta)\\}/{n_{\mathrm{U}}}}$ and
$\zeta_{0}:=\min\\{\omega_{0}\gamma^{-1/2},\omega_{0}^{1/2}\\}$ throughout
this section. Thus, treating $d$ as a constant and ignoring polylogarithmic
terms, $\omega_{0}$ is of order $n_{\mathrm{U}}^{-1/2}$ and $\zeta_{0}$ is of
order $\min\\{n_{\mathrm{L}}^{-1/2},n_{\mathrm{U}}^{-1/4}\\}$ when
$\gamma<1/2$. We remark that $n_{\mathrm{L}}^{-1/2}$ is the critical
$\ell_{2}$-testing radius for distinguishing the means of two labeled Gaussian
distributions with identity covariance using $n_{\mathrm{L}}$ observations. On
the other hand, as we show in Lemma 16, no test of the null hypothesis
$H_{0}:\mathcal{N}_{d}(0,I_{d})$ against the two-component mixture alternative
$H_{1}:\frac{1}{2}\mathcal{N}_{d}(\mu^{*},I_{d})+\frac{1}{2}\mathcal{N}_{d}(-\mu^{*},I_{d})$
based on $n_{\mathrm{U}}$ observations can have large power unless the signal
strength $\|\mu^{*}\|$ is at least of order $n_{\mathrm{U}}^{-1/4}$.

######  Proposition 5.

Fix $\delta\in(2e^{-n},1]$ and $r\geq 1$, and suppose that $\|\mu^{*}\|\leq r$
and $\gamma<1/2$. There exists $c>0$, depending only on $r$, such that if
$\omega_{0}\leq c$ and $n\geq 3$, then the following statements hold:

  1. (i)

For any $\hat{\mu}^{(0)}\in\mathbb{R}^{d}$ with $\|\hat{\mu}^{(0)}\|\leq r+3$,
we have with probability at least $1-2\delta$ that

$\limsup_{t\to\infty}L(\hat{\mu}^{(t)},\mu^{*})\lesssim_{r}\zeta_{0}\vee\|\mu^{*}\|.$

  2. (ii)

There exists $C>0$, depending only on $r$, such that if $\|\mu^{*}\|\geq
C\zeta_{0}\sqrt{d\log n}$ and $\hat{\mu}^{(0)}=(\zeta_{0}\vee
r\omega_{0})\eta_{0}$ with $\eta_{0}\sim\mathrm{Unif}(\mathbb{S}^{d-1})$, then
with probability at least $1-2\delta-\sqrt{2/(\pi\log n_{\mathrm{U}})}$, we
have

$\limsup_{t\to\infty}L(\hat{\mu}^{(t)},\mu^{*})\lesssim_{r}\frac{\omega_{0}}{\|\mu^{*}\|}\wedge\frac{\omega_{0}}{\gamma^{1/2}}.$

In order to interpret Proposition 5 _(i)_ , consider the regime where
$\|\mu^{*}\|\leq\zeta_{0}$. In this case, as discussed above, the two mixture
components are essentially indistinguishable, and the bound reveals that the
EM algorithm performs no worse than the trivial zero estimator, up to constant
factors. On the other hand, part  _(ii)_ studies the more interesting regime
where the two mixture components are distinguishable, and we establish a
faster convergence rate for the EM algorithm in this strong signal regime.

The following theorem combines the two convergence regimes in Proposition 5 to
derive a convergence guarantee for the estimated whitened between-class
covariance matrix output by Algorithm 3. To state the result, recall the
definition of $\mathcal{C}$ from (11). For any $\zeta>0$, we write $U(\zeta)$
for the pushforward measure on $\mathcal{C}$ induced by
$\mathrm{Unif}(\zeta\mathbb{S}^{d-1})$ under the map
$\mu\mapsto(-\mu,\mu,I_{d})$.

######  Theorem 6.

Fix $\delta\in(2e^{-n},1]$, and $r\geq 1$ and suppose that $\|\mu^{*}\|\leq r$
and $\gamma<1/2$. There exists $c>0$, depending only on $r$, such that if
$\omega_{0}\leq\min\\{c,(d\log n)^{-3}\\}$ and $n\geq 108$, then the sequence
of outputs $(\hat{Q}^{(T)})_{T\in\mathbb{N}}$ of Algorithm 3 with inputs
$(Z_{1},Y_{1}),\ldots,(Z_{n},Y_{n})$, $\mathcal{C}$,
$\pi_{\mathcal{C}}=U(\zeta_{0}\vee r\omega_{0})$, $M\in\mathbb{N}$ and
$T\in\mathbb{N}$ satisfies with probability at least $1-3\delta-e^{-M/50}$
that

$\limsup_{T\to\infty}\|\hat{Q}^{(T)}-\mu^{*}\mu^{*\top}\|_{\mathrm{op}}\lesssim_{r}\frac{\omega_{0}}{\|\mu^{*}\|}\wedge\zeta_{0}.$

Finally in this section, we study the implications of Theorem 6 for the
recovery of the signal coordinates in the semi-supervised learning setting. We
write $\psi^{(M,T)}$ for the base procedure that takes
$(z_{i},y_{i})_{i\in[n]}\in\mathbb{R}^{d}\times([K]\cup\\{0\\})$ as input and
returns the output of Algorithm 3 when run with these inputs together with
$\mathcal{C}$, $\pi_{\mathcal{C}}$, $M$ and $T$. Let $S_{0}$ denote the set of
coordinates where $\nu^{*}\in\mathbb{R}^{p}$ is non-zero, and let
$s_{0}:=|S_{0}|$.

######  Corollary 7.

Fix $\varepsilon\in(8e^{-n/2},1]$, $r\geq 1$, and suppose that
$\|\mu^{*}\|\leq r$, $M\geq 50\log(4/\varepsilon)+50d\log p$ and $\gamma<1/2$.
Let $\nu^{*}_{\max}:=\|\nu^{*}\|_{\infty}$ and let $\nu_{\min}^{*}$ denote the
minimum absolute value of a non-zero component of $\nu^{*}$. There exist
$C_{1},C_{2}>0$, depending only on $r$, such that if $n\geq C_{1}(d\log
p)^{6}\\{d\log p+\log(1/\varepsilon)\\}$, and

$C_{2}\min\biggl{[}\biggl{\\{}\frac{d\log(p\vee
n)+\log(1/\varepsilon)}{n}\biggr{\\}}^{1/4},\sqrt{\frac{d\log(p\vee
n)+\log(1/\varepsilon)}{n_{\mathrm{L}}}}\biggr{]}\leq\frac{(\nu^{*}_{\min})^{2}}{4},$

then the sequence of outputs $(\hat{S}^{(T)})_{T\geq 1}$ of Algorithm 1 with
inputs $K=2$, $p$, $d\geq s_{0}$, $\ell\geq s_{0}$, $(X_{i},Y_{i})_{i\in[n]}$,
$A$, $B$ and base procedure $\psi^{(M,T)}$ satisfies

$\liminf_{T\to\infty}\mathbb{P}(S_{0}\subseteq\hat{S}^{(T)})\geq
1-\varepsilon-pe^{-A(\nu^{*}_{\min})^{4}/(50p^{2}(\nu^{*}_{\max})^{4})}.$

Corollary 7 reveals in particular that, treating $\nu^{*}_{\max}$ and
$\nu^{*}_{\min}$ as constants and under the stated sample size conditions, we
again recover all of the signal coordinates in the top $s_{0}$ output entries,
provided that $A$ is large by comparison with $p^{2}$. Thus, in this sense, we
can achieve a similar guarantee to that provided by Corollary 4, though the
number of groups of projections required for a high probability guarantee in
Corollary 7 may be significantly larger in settings where the ratio
$\nu_{\max}^{*}/\nu_{\min}^{*}$ is large.

##  4 Numerical studies

Throughout this section, unless otherwise stated, data
$(X_{i},Y_{i},Y_{i}^{*})_{i\in[n]}$ are sampled from an equal-probability
normal mixture as follows: $\mathbb{P}(Y_{i}^{*}=k)=1/K$ for $k\in[K]$,
$\mathbb{P}(Y_{i}=Y_{i}^{*})=1-\mathbb{P}(Y_{i}=0)=\gamma$ and $X_{i}\mid
Y_{i}^{*}\sim N_{p}(\mu_{Y_{i}^{*}},\Sigma_{\mathrm{w}})$. The cluster means
$(\mu_{k})_{k\in[K]}$ are chosen to be $s_{0}$-sparse and we define the
signal-to-noise ratio of the problem to be333In some of our simulations,
$\Sigma_{\mathrm{w}}$ was generated randomly for convenience. In such
settings, we replaced $\mathrm{tr}(\Sigma_{\mathrm{w}})/p$ in the denominator
of (13) with $\mathbb{E}\\{\mathrm{tr}(\Sigma_{\mathrm{w}})\\}/p$.

$\mathrm{SNR}:=\frac{\min_{k,k^{\prime}\in[K],k\neq
k^{\prime}}\|\mu_{k}-\mu_{k^{\prime}}\|}{\sqrt{\mathrm{tr}(\Sigma_{\mathrm{w}})/p}}.$
(13)

In our numerical studies, we slightly modify Algorithm 3 so that instead of
randomly initializing the cluster means and the covariance matrix, we use the
output of hierarchical clustering to initialize the EM algorithm as
implemented in the mclust R package (Fraley and Raftery, 1998). This allow us
to run Algorithm 3 with $M=1$.

###  4.1 Choice of tuning parameters

The purpose of this subsection is to investigate the effect of the various
input parameters $A$, $B$, $d$ and $\ell$ in Algorithm 1, and to recommend
sensible default choices. In Figure 1, we plot the misclustering rate with
Algorithm 3 as a base procedure in our Gaussian semi-supervised learning
setting as each of these parameters varies, for four different $\mathrm{SNR}$
levels. After applying Algorithm 1, we obtain our final estimated cluster
labels by using Algorithm 3 again on the data projected onto the selected
coordinates with a single hierarchical clustering initialization. We then
output the predicted labels, computed as
$\hat{y}_{i}:=\operatorname*{sargmax}_{k\in[K]}L_{i,k}$ for $i\in[n]$, instead
of $\hat{Q}$.

The panels in Figure 1 reveal that the misclustering rate is quite robust to
the choices of $A$, $B$ and $d$, and that it is less serious (and may even
help) to choose $\ell$ larger—rather than smaller—than $s_{0}$. In particular,
it seems that $A=150$ suffices for almost optimal performance (though there
appears to be some penalty for choosing it to be as small as $50$), and $B=75$
appears adequate. There is no clear trend on performance with the choice of
$d$, so for simplicity we took $d=s_{0}$ in our remaining simulations below.
Finally, the misclustering rate appears to decrease as $\ell$ increases, with
an elbow in the curve visible at the highest value of the $\mathrm{SNR}$ when
$\ell$ is set to the true sparsity level $s_{0}$. Of course, if $\ell$ is
chosen to be very large, then we will include many noise variables, and the
misclustering rate will eventually deteriorate. Nevertheless, the bottom-right
panel of Figure 1 indicates that the gain in increasing the probability of
including all signal variables may outweigh the penalty of also including more
noise variables—as expected, this effect is larger when the $\mathrm{SNR}$ is
larger. For simplicity we choose $\ell=s_{0}$ in our remaining simulations,
though we recommend practitioners err on the side of choosing larger $\ell$.

|
---|---
|
Figure 1: Average misclustering rate over 200 repetitions in our anisotropic
Gaussian semi-supervised learning problem with $n=250$, $p=600$, $s=4$, $K=3$,
$\gamma=0.05$, $\mathrm{SNR}\in\\{2.5,3,3.5,4\\}$ and
$\Sigma_{\mathrm{w}}=V\Lambda V^{\top}$, where $\Lambda\in\mathbb{R}^{p\times
p}$ is diagonal with independent $\mathrm{Unif}[0,2]$ diagonal entries and $V$
is independent of $\Lambda$, and generated according to the Haar measure on
$\mathbb{O}^{p\times p}$. For each of the four panels, we fix three of $d=4$,
$\ell=4$, $A=150$, $B=75$, and vary the remaining one. The shaded regions
represent interpolated 95% confidence intervals at each of the points.

###  4.2 Comparison with existing methods

In this subsection, we compare the empirical performance of the Sharp-SSL
algorithm in high-dimensional clustering tasks with several existing
approaches. We apply the Sharp-SSL algorithm using the EM algorithm of
Algorithm 3 as a base procedure, with input parameters $A=150$, $B=75$,
$d=\ell=s_{0}$ as discussed in Section 4.1, and our final estimated cluster
labels are then obtained as described there.

$n=250$, $p=200$, isotropic |  $n=250$, $p=600$, anisotropic
---|---
|
$\mathrm{SNR}=3$, $p=200$, isotropic |  $\mathrm{SNR}=3$, $p=600$, anisotropic
|
Figure 2:  Average misclustering rate over 100 repetitions using Sharp-SSL
followed by the EM algorithm, as well as using the other methods from Section
4.2. Data are generated from the normal mixture distribution described at the
beginning of Section 4 with $K=3$ and $p=200$ (left) as well as $p=600$
(right). The three cluster means are given by $\mu_{1}=a(1,1,0,\bm{0}_{p-3})$,
$\mu_{2}=a(-1,0,1,\bm{0}_{p-3})$ and $\mu_{3}=a(0,-1,-1,\bm{0}_{p-3})$, where
the scale $a$ is chosen such that their pairwise distances are all equal to
$\mathrm{SNR}$. For isotropic settings (left), $\Sigma_{\mathrm{w}}=I_{p}$;
for anisotropic settings (right), $\Sigma_{\mathrm{w}}=V\Lambda V^{\top}$,
where $\Lambda\in\mathbb{R}^{p\times p}$ is diagonal with independent
$\mathrm{Unif}[0,2]$ diagonal entries and $V$ is independent of $\Lambda$, and
sampled from the Haar measure on $\mathbb{O}^{p\times p}$. The Bayes risk is
shown as the gray dashed line. In the top panels, $n=250$ and the
$\mathrm{SNR}$ varies; in the bottom panels, $\mathrm{SNR}=3$ and $n$ varies.
The shaded regions represent interpolated 95% confidence intervals at each of
the points.

We compare the Sharp-SSL algorithm with five alternative high-dimensional
clustering methods: spectral clustering (e.g. von Luxburg, 2007), the
$\ell_{1}$-penalized approach of Witten and Tibshirani (2010) and the RPEClus
algorithm of Anderlucci, Fortunato and Montanari (2022) as well as a pair of
methods that, like Sharp-SSL, apply dimension reduction prior to a low-
dimensional clustering algorithm.

In more detail, the spectral clustering approach first constructs a
$J$-nearest neighbour graph adjacency matrix
$A=(A_{i,i^{\prime}})_{i,i^{\prime}\in[n]}\in\\{0,1\\}^{n\times n}$, where
$A_{i,i^{\prime}}:=1$ if either $X_{i}$ is one of the $J=10$ nearest
neighbours of $X_{i^{\prime}}$ in Euclidean distance or vice versa, and
$A_{i,i^{\prime}}:=0$ otherwise. It then computes an $n\times K$ matrix of
eigenvectors associated with the $K$ smallest nonzero eigenvalues of the
Laplacian matrix $L:=D-A$, where $D\in\mathbb{R}^{n\times n}$ is a diagonal
matrix with diagonal entries
$D_{i,i}:=\sum_{i^{\prime}\in[n]}A_{i,i^{\prime}}$. The final step is to apply
the $K$-means clustering algorithm (Lloyd, 1982), as implemented in the kmeans
base R function with 100 random initializations, to the rows of $L$ with the
oracle choice of $K$.

The Witten and Tibshirani (2010) method, which is implemented in the sparcl R
package, determines the estimated cluster memberships by maximizing a
coordinatewise-weighted between-cluster sum of squares criterion, subject to
an $\ell_{1}$ constraint on the weights. A permutation approach is used to
select the $\ell_{1}$ tuning parameter.

In the RPEClus algorithm of Anderlucci, Fortunato and Montanari (2022), we
generate $B$ random orthogonal projections and incorporate the $d$-dimensional
projected data as covariates for a linear regression with the orthogonal
complement of the projected data as the response. We then use the Bayesian
Information Criteria (BIC) from both an application of the EM algorithm to the
projected data and the aforementioned regression to identify good projections,
and aggregate using the consensus clustering technique of Dimitriadou,
Weingessel and Hornik (2002) over the best $B^{*}$ projections chosen
according to the sum of the BIC scores. Following the recommendation of
Anderlucci, Fortunato and Montanari (2022), we took $B=1000$ and $B^{*}=100$
as well as $d=s_{0}$. It turned out that this approach had a misclustering
rate almost identical to that of a random guess, likely because it did not
leverage the sparsity of the signal. We therefore modified this method by
generating random axis-aligned projections instead of orthogonal ones, and
report this version in our comparison.

The first of the two-stage approaches applies principal component analysis
(PCA) to project the data into the oracle choice of $K-1$ dimensions (the
dimension of the space spanned by the $K$ cluster means); the second uses
sparse principal component analysis (SPCA), as implemented in the SPCAvRP
algorithm (Gataric, Wang and Samworth, 2020) with inputs $A=600$, $B=200$, and
the oracle choices $d=\ell=s_{0}$, to project into $s_{0}$ dimensions.
Thereafter, both algorithms apply $K$-means to the projected data as above. We
also explored the option of replacing the $K$-means steps in these latter
algorithms with the EM algorithm, but observed very little difference, so do
not report these results here.

Given true labels $y_{1},\ldots,y_{n}\in[K]$ and estimated labels
$\hat{y}_{1},\ldots,\hat{y}_{n}\in[K]$ from a clustering algorithm, we measure
the performance of the algorithm via its _misclustering rate_ , defined
as444Here, the minimum over permutations is taken because it is only the
cluster groupings, and not the labels themselves, that are important.

$L(\\{y_{1},\ldots,y_{n}\\},\\{\hat{y}_{1},\ldots,\hat{y}_{n}\\}):=\min_{\sigma\in\mathcal{S}_{K}}\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}_{\\{\sigma(\hat{y}_{i})\neq
y_{i}\\}},$

where $S_{K}$ is the group of all permutations of $[K]$. In particular, Figure
2 presents the average misclustering rates over 100 Monte Carlo repetitions of
the different high-dimensional clustering algorithms described above. Across
two different dimensions $p\in\\{200,600\\}$, isotropic and anisotropic
settings, and for different values of $n$ and $\mathrm{SNR}$, we see a
consistent picture of the Sharp-SSL algorithm combined with EM producing the
lowest misclustering rates, often by a large margin. Indeed, for all but the
smallest sample sizes or values of $\mathrm{SNR}$, the Sharp-SSL+EM algorithm
nearly attains the Bayes risk in all of the problems considered here.

###  4.3 Effect of observed fraction on misclustering rate

One of the key attractions of our procedure is that it offers a unified
framework to perform classification or clustering with an arbitrary fraction
of labeled observations. In this subsection, we explore the performance of the
algorithm as we vary the proportion of observed labels.

$\mathrm{SNR}=3.5$, $p=200$, $s=3$ |  $\mathrm{SNR}=3.5$, $p=600$, $s=3$
---|---
|
$\mathrm{SNR}=1.5$, $p=s=3$ |  $\mathrm{SNR}=2$, $p=s=3$
|
Figure 3:  Effect of label fraction on performance of supervised, semi-
supervised and unsupervised Sharp-SSL learning methods. Data are generated
from the normal mixture distribution described at the beginning of Section 4
with $K=2$ and $\Sigma_{\mathrm{w}}=I_{p}$,
$\mu_{1}=-\mu_{2}=a(\bm{1}_{s},\bm{0}_{p-s})^{\top}\in\mathbb{R}^{p}$, where
$a$ is chosen such that $\|\mu_{1}-\mu_{2}\|=\mathrm{SNR}$. Bottom: average
Frobenius loss of estimating the $(\mu_{1},\mu_{2})\in\mathbb{R}^{p\times 2}$
over 100 repetitions via the semi-supervised approach (Algorithm 3),
supervised approach (Algorithm 2) and unsupervised approach (Algorithm 3
without using the labels). Top: average misclustering rate over 100
repetitions from applying the above three methods as base algorithms in
Algorithm 1. The shaded regions represent interpolated 95% confidence
intervals at each of the points.

Recall that we have two different options for the way in which we implement
the Sharp-SSL algorithm to estimate the set of signal coordinates: we can
either use only the labeled data, as in the supervised learning approach of
Algorithm 2, or we can try to leverage in addition the unlabeled data via the
semi-supervised EM approach of Algorithm 3. In the extreme case of this latter
version, we have no labeled data, so the algorithm is unsupervised. In Figure
3 we compare the performance of these three methods in both high- and low-
dimensional versions of the normal mixture distribution data generation
mechanism described at the beginning of Section 4 as the proportion $\gamma$
of observed labels varies.

More precisely, for the semi-supervised and unsupervised algorithms, we adopt
the same implementation of Sharp-SSL as described at the beginning of Section
4.2. The supervised algorithm is very similar, but applies Algorithm 2 in
place of Algorithm 3 to select coordinates, and obtains final predicted labels
by applying LDA again on the projected labeled data. In cases where the
proportion of labeled data was so small that the convex hull of the projected
labeled data was less than full-dimensional for every class, we forced
Algorithm 2 to return a zero matrix (this only happened when $\gamma$ was very
small).

The top panels of Figure 3 present the results in high-dimensional settings
with $p\in\\{200,600\\}$. Since the unsupervised approach has no access to the
labels, it has constant misclustering rate. The performance of the semi-
supervised approach is always at least as good as that of the unsupervised
algorithm, and improves as $\gamma$ increases. In other words, it effectively
leverages the additional information provided by the class labels. When
$\gamma$ is very small, the supervised algorithm—which ignores the unlabeled
data—is inaccurate, as it has very little data to work with. On the other
hand, its performance also improves as $\gamma$ increases, and once around
5$\%$ of our data are labeled, it outperforms the unsupervised algorithm.
Further, it essentially matches the semi-supervised approach when about a
third of the data are labeled. We truncate the plot at $\gamma=1/2$ to ensure
that we have enough test data on which to compute the misclustering rate.

In the bottom panels of Figure 3, we explore the performance of the three
algorithms above in two low-dimensional settings with different values of SNR,
in order to provide further insight into the phenomena described in the
previous paragraph. Here, we take $K=2$ and report the average Frobenius norm
loss

$\mathcal{L}\bigl{(}(\hat{\mu}_{1},\hat{\mu}_{2}),(\mu_{1},\mu_{2})\bigr{)}:=\min\bigl{\\{}\|(\hat{\mu}_{1},\hat{\mu}_{2})-(\mu_{1},\mu_{2})\|_{\mathrm{F}},\|(\hat{\mu}_{2},\hat{\mu}_{1})-(\mu_{1},\mu_{2})\|_{\mathrm{F}}\bigr{\\}}$

of the estimated means, over 100 repetitions. If there are insufficient
labeled data to run Algorithm 2, then we output
$\hat{\mu}_{1}=\hat{\mu}_{2}=\bm{0}_{p}$. We see that, already in these low-
dimensional problems, a similar picture emerges: if the proportion of labeled
data is small, then the unsupervised algorithm outperforms the supervised one,
but this situation may be reversed when $\gamma$ is larger. The semi-
supervised algorithm is able to leverage both the unlabeled and labeled data
to obtain the best of both worlds. These empirical observations agree with our
theory from Section 3, in particular in the way in which Theorem 6 bounds the
accuracy of mean estimation for the semi-supervised algorithm by a minimum of
a term that does not depend on $\gamma$ and one that decreases as $\gamma$
increases. It appears that the switch in the minimum occurs around
$\gamma=0.02$ in these examples.

###  4.4 Empirical data analysis

In this subsection we apply the Sharp-SSL algorithm, as well as several
competing methods, to the gene expression data set from Alon et al. (1999),
which contains observations on 62 patients. A preprocessed version of the data
can be downloaded from the R package ‘datamicroarray’ (Ramey, 2016), with a
total of $2000$ features (genes) measured on $40$ patients with colon tumors
and $22$ without tumors. We first exclude 9 genes to remove perfect
collinearity and then standardize each of the remaining $p=1991$ columns of
the dataset to have unit variance.

We apply the Sharp-SSL algorithm using EM (Algorithm 3) as the base procedure,
with input parameters $A=150$, $B=75$, $d=\ell=5$. In addition to our approach
(Sharp-SSL+EM), we also compare the performance of the spectral clustering
(SC) method, the Witten and Tibshirani (2010) method (WT2010, as well as four
two-stage methods (PCA+Kmeans, PCA+EM, SPCA+Kmeans, SPCA+EM), where we first
reduce dimension of the data to a $5$-dimensional subspace using either PCA or
SPCA and then apply either the EM algorithm or $K$-means clustering on the
low-dimensional data. For SPCA, we use the SPCAvRP algorithm (Gataric, Wang
and Samworth, 2020) with inputs $A=600$, $B=200$ and $d=\ell=5$. The true
labels are hidden to all algorithms and are only used to evaluate the final
misclustering rate.

Over 100 Monte Carlo repetitions of the randomized algorithms, the Sharp-
SSL+EM method had an average misclustering rate of 28.8%, whereas all other
competitors had a misclustering rate above 40%, as can be seen from the right-
hand data points in Figure 4. To investigate this performance further, we
applied each method to a subset of the features. These were constructed from
the top $\ell=5$ genes identified through Sharp-SSL, together with
$m=0,10,50,200$ and $600$ randomly chosen genes from the remaining $1986$. The
results are presented as the other data points in Figure 4. We see that the
improved performance of the Sharp-SSL+EM relative to the other methods
persists, even when only a small number of potentially non-discriminative
covariates are present. When $m=0$, Sharp-SSL+EM has a slight disadvantage as
other algorithms benefit from the ensemble effect of combining two different
learning methods; nevertheless it remains competitive. This reinforces the
point that the primary contribution of the Sharp-SSL algorithm is to identify
signal coordinates that are helpful for semi-supervised learning, and once
this task has been accomplished, a variety of low-dimensional procedures are
available to the practitioner.

Figure 4: Average misclustering rate (over 100 repetitions for randomized
algorithms) for the colon tumor data, using Sharp-SSL followed by the EM
algorithm, as well as the other methods described in Section 4.4. The right-
hand data points plot the average misclustering rate on the full data set. The
other points were obtained by applying each method to a subset of genes formed
from the top five genes identified by Sharp-SSL together with randomly sampled
genes. The shaded regions represent interpolated 95% confidence intervals at
each of the points.

##  5 Proof of the main results

###  5.1 Proof of Proposition 1

The assumption that the convex hull of $\nu_{1},\ldots,\nu_{K}$ is
$(K-1)$-dimensional implies that $\Sigma_{\mathrm{b}}$ is of rank $K-1$.
Define
$A:=\Sigma_{\mathrm{w}}^{-1/2}\Sigma_{\mathrm{b}}\Sigma_{\mathrm{w}}^{-1/2}\in\mathbb{R}^{p\times
p}$, which has rank $K-1$. Given $V\in\mathbb{O}^{p\times d}$, we can find
$Q\in\mathbb{O}^{p\times d}$ with the same column span as that of
$\Sigma_{\mathrm{w}}^{1/2}V$ and let
$R:=Q^{\top}\Sigma_{\mathrm{w}}^{1/2}V\in\mathbb{R}^{d\times d}$, so that $R$
is invertible, and $\Sigma_{\mathrm{w}}^{1/2}V=QR$. We observe that

$\mathrm{tr}\\{(V^{\top}\Sigma_{\mathrm{w}}V)^{-1}(V^{\top}\Sigma_{\mathrm{b}}V)\\}=\mathrm{tr}\\{(R^{\top}R)^{-1}(R^{\top}Q^{\top}AQR)\\}=\mathrm{tr}(Q^{\top}AQ).$

Thus, $J(V;\Sigma_{\mathrm{b}},\Sigma_{\mathrm{w}})$ depends on $V$ only
through the column space of $\Sigma_{\mathrm{w}}^{1/2}V$. Moreover,
$\mathrm{tr}(Q^{\top}AQ)$ is maximized when $Q$, or equivalently
$\Sigma_{\mathrm{w}}^{1/2}V$, spans a $d$-dimensional space that contains the
$(K-1)$-dimensional eigenspace corresponding to the non-zero eigenvalues of
$A$. Note that if for some $v\in\mathbb{R}^{p}\setminus\\{0\\}$ and
$\lambda\geq 0$, we have $Av=\lambda v$, then
$\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}\Sigma_{\mathrm{w}}^{-1/2}v=\Sigma_{\mathrm{w}}^{-1/2}Av=\lambda\Sigma_{\mathrm{w}}^{-1/2}v$,
so $\Sigma_{\mathrm{w}}^{-1/2}v$ is an eigenvector of
$\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$ with eigenvalue $\lambda$. Hence
$V$ maximizes $J(V;\Sigma_{\mathrm{b}},\Sigma_{\mathrm{w}})$ over
$\mathbb{O}^{p\times d}$ if and only if $V$ spans a $d$-dimensional space that
contains the $(K-1)$-dimensional eigenspace corresponding to the $K-1$ non-
zero eigenvalues of $\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$. Finally,
for any $v\in\mathbb{R}^{p}\setminus\\{0\\}$,

$v^{\top}\Sigma_{\mathrm{w}}^{-1/2}\Sigma_{\mathrm{b}}\Sigma_{\mathrm{w}}^{-1/2}v=\sum_{k=1}^{K}\pi_{k}v^{\top}\Sigma_{\mathrm{w}}^{-1/2}(\nu_{k}-\nu)(\nu_{k}-\nu)^{\top}\Sigma_{\mathrm{w}}^{-1/2}v\neq
0$

if and only if $v^{\top}\Sigma_{\mathrm{w}}^{-1/2}(\nu_{k}-\nu)\neq 0$ for
some $k\in[K]$. Thus, the eigenspace corresponding to the non-zero eigenvalues
of $A$ is spanned by
$\bigl{(}\Sigma_{\mathrm{w}}^{-1/2}(\nu_{k}-\nu):k\in[K]\bigr{)}$, and so the
eigenspace corresponding to non-zero eigenvalues of
$\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}}$ is spanned by
$\bigl{(}\Sigma_{\mathrm{w}}^{-1}(\nu_{k}-\nu):k\in[K]\bigr{)}$.

###  5.2 Proof of Theorem 2

We write $S^{a,b}:=\\{j\in[p]:(P^{a,b,\top}P^{a,b})_{j,j}=1\\}$. For any
$S=\\{j_{1},\ldots,j_{d}\\}^{\top}\in\binom{[p]}{d}$, we identify the set $S$
with the sequence $j_{i_{1}}<\ldots<j_{i_{d}}$ sorted in increasing order; and
with a slight abuse of notation, we will use $S$ to refer to either object,
which will always be clear depending on the context. We define
$P^{S}\in\mathcal{P}_{d}$ by
$(P^{S})_{\ell,j}:=\mathbbm{1}_{\\{j=j_{\ell}\\}}$, so that
$P^{S,\top}P^{S}=\mathrm{diag}\bigl{(}(\mathbbm{1}_{\\{j\in
S\\}})_{j\in[p]}\bigr{)}$. Define
$Q^{S}:=(P^{S}\Sigma_{\mathrm{w}}P^{S,\top})^{-1}P^{S}\Sigma_{\mathrm{b}}P^{S,\top}\in\mathbb{R}^{d\times
d}$ and
$\hat{Q}^{S}:=\psi\bigl{(}(P^{S}X_{i},Y_{i})_{i\in[n]}\bigr{)}\in\mathbb{R}^{d\times
d}$. Note that $\hat{Q}^{a,b}=\hat{Q}^{S^{a,b}}$ in this notation, and we will
similarly denote $Q^{a,b}:=Q^{S^{a,b}}$ for simplicity. Recalling the
definition of $\Omega$ from (8), in our new notation, and recalling that
$\psi$ is permutation-equivariant, we can write

$\displaystyle\Omega=\biggl{\\{}\max_{S\in\binom{[p]}{d}}\bigl{\|}\hat{Q}^{S}-Q^{S}\bigr{\|}_{\mathrm{op}}<\frac{\gamma_{\min}}{4(K-1)}\biggr{\\}},$

and have $\mathbb{P}(\Omega)\geq 1-\varepsilon$ by (7). We will work on the
event $\Omega$ throughout the remainder of the proof, and assume also that
$\gamma_{\min}>0$, because otherwise the conclusion is trivial.

By Weyl’s inequality (Weyl, 1912; Stewart and Sun, 1990, Corollary IV.4.9), we
have on $\Omega$ that

$\displaystyle|\mathrm{tr}(\hat{Q}^{S})-\mathrm{tr}(Q^{S})|\leq(K-1)\|\hat{Q}^{S}-Q^{S}\|_{\mathrm{op}}\leq\frac{\gamma_{\min}}{4}.$

On the other hand, we note that $\mathrm{tr}(Q^{S})=\sum_{j\in S\cap
S_{0}}(\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}})_{j,j}$. Therefore, by the
triangle inequality, for any $S,S^{\prime}\in\binom{[p]}{d}$ such that $S\cap
S_{0}$ is a proper subset of $S^{\prime}\cap S_{0}$, we have on $\Omega$ that

$\mathrm{tr}(\hat{Q}^{S})-\mathrm{tr}(\hat{Q}^{S^{\prime}})\leq\frac{\gamma_{\min}}{2}-\\!\\!\sum_{j\in(S^{\prime}\setminus
S)\cap S_{0}}(\Sigma_{\mathrm{w}}^{-1}\Sigma_{\mathrm{b}})_{j,j}<0.$ (14)

Fix $a\in[A]$, and for any $\tilde{j}\in[p]$, define
$q_{\tilde{j}}:=\mathbb{P}\bigl{(}\tilde{j}\in
S^{a,b^{*}(a)}\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}$. Now fix some $j\in S_{0}$
and $j^{\prime}\in[p]\setminus S_{0}$. We claim that

$(q_{j}-q_{j^{\prime}})\mathbbm{1}_{\Omega}\geq 0.$ (15)

To verify this claim, define for $\tilde{j}\in\\{j,j^{\prime}\\}$ and
$b\in[B]$ the sets

$\mathcal{S}_{b,\tilde{j}}:=\bigl{\\{}(S^{a,1},\ldots,S^{a,B}):b^{*}(a)=b,\tilde{j}\in
S^{a,b}\bigr{\\}}\quad\text{and}\quad\mathcal{S}_{b}:=\bigl{\\{}(S^{a,1},\ldots,S^{a,B}):b^{*}(a)=b\bigr{\\}}.$

Let $f:\binom{[p]}{d}\to\binom{[p]}{d}$ be a map defined by

$f(S):=\begin{cases}(S\setminus\\{j^{\prime}\\})\cup\\{j\\}&\text{if $j\notin
S$ and $j^{\prime}\in S$}\\\ S&\text{otherwise}.\end{cases}$

If $j\notin S^{a,b}$ and $j^{\prime}\in S^{a,b}$, then $f(S^{a,b})\cap
S_{0}=(S^{a,b}\cup\\{j\\})\cap S_{0}=(S^{a,b}\cap S_{0})\cup\\{j\\}$, so
$S^{a,b}\cap S_{0}$ is a proper subset of $f(S^{a,b})\cap S_{0}$; on the other
hand, if either $j\in S^{a,b}$ or $j,j^{\prime}\notin S^{a,b}$, then
$f(S^{a,b})=S^{a,b}$. It follows by (14) that on $\Omega$ we have

$\mathrm{tr}(\hat{Q}^{S^{a,b}})\leq\mathrm{tr}(\hat{Q}^{f(S^{a,b})}).$ (16)

Now let $F:\mathcal{S}_{b,j^{\prime}}\to\mathcal{S}_{b,j}$ be defined as

$F(S^{a,1},\ldots,S^{a,B}):=\bigl{(}S^{a,1},\ldots,S^{a,b-1},f(S^{a,b}),S^{a,b+1},\ldots,S^{a,B}\bigr{)}.$

We claim that $F$ is both well-defined and injective on $\Omega$. For the
first of these claims, we note that since $j^{\prime}\in S^{a,b}$, we must
have $j\in f(S^{a,b})$. Moreover, if
$(S^{a,1},\ldots,S^{a,B})\in\mathcal{S}_{b,j^{\prime}}$, then $b^{*}(a)=b$.
But (16) holds on $\Omega$, so
$\bigl{(}S^{a,1},\ldots,S^{a,b-1},f(S^{a,b}),S^{a,b+1},\ldots,S^{a,B}\bigr{)}\in\mathcal{S}_{b,j}$.
Hence $F$ is well-defined. For the second claim, suppose that
$S_{1},S_{2}\in\binom{[p]}{d}$ are such that $j^{\prime}\in S_{1}\cap S_{2}$
and $f(S_{1})=f(S_{2})$. If $j\in S_{1}\cap S_{2}$, then
$S_{1}=f(S_{1})=f(S_{2})=S_{2}$; if $j\in S_{1}$ but $j\notin S_{2}$, then
$j\in f(S_{2})\setminus f(S_{1})$, a contradiction. Similarly, we cannot have
$j\notin S_{1}$ but $j\in S_{2}$. Finally, if $j\in S_{1}^{\mathrm{c}}\cap
S_{2}^{\mathrm{c}}$, then
$S_{1}=\bigl{(}f(S_{1})\setminus\\{j\\}\bigr{)}\cup\\{j^{\prime}\\}=\bigl{(}f(S_{2})\setminus\\{j\\}\bigr{)}\cup\\{j^{\prime}\\}=S_{2}$.
We deduce that $f$ is injective on $\\{S:j^{\prime}\in S\\}$. Since
$j^{\prime}\in S^{a,b}$ for
$(S^{a,1},\ldots,S^{a,B})\in\mathcal{S}_{b,j^{\prime}}$, this establishes the
injectivity of $F$. In particular,
$|\mathcal{S}_{b,j^{\prime}}|\leq|\mathcal{S}_{b,j}|$. Consequently, on
$\Omega$, we have for all $b\in[B]$ that

$\displaystyle\mathbb{P}\bigl{(}j\in
S^{a,b^{*}(a)}\mid(X_{i},Y_{i})_{i\in[n]},b^{*}(a)=b\bigr{)}$
$\displaystyle=\frac{\mathbb{P}\bigl{(}j\in
S^{a,b^{*}(a)},b^{*}(a)=b\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}}{\mathbb{P}\bigl{(}b^{*}(a)=b\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}}=\frac{|\mathcal{S}_{b,j}|}{|\mathcal{S}_{b}|}$
$\displaystyle\geq\frac{|\mathcal{S}_{b,j^{\prime}}|}{|\mathcal{S}_{b}|}=\frac{\mathbb{P}\bigl{(}j^{\prime}\in
S^{a,b^{*}(a)},b^{*}(a)=b\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}}{\mathbb{P}\bigl{(}b^{*}(a)=b\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}}$
$\displaystyle=\mathbb{P}\bigl{(}j^{\prime}\in
S^{a,b^{*}(a)}\mid(X_{i},Y_{i})_{i\in[n]},b^{*}(a)=b\bigr{)},$

which implies Claim (15). We remark that one consequence of (15) is that,
since $d\geq s_{0}$, we have on $\Omega$ that

$q_{j}\geq\frac{\sum_{\tilde{j}\in([p]\setminus
S_{0})\cup\\{j\\}}q_{\tilde{j}}}{p-s_{0}+1}=\frac{d-\sum_{\tilde{j}\in
S_{0}\setminus\\{j\\}}q_{\tilde{j}}}{p-s_{0}+1}\geq\frac{d-s_{0}+1}{p-s_{0}+1}\geq\frac{1}{p}.$
(17)

Again fixing $j\in S_{0}$ and $j^{\prime}\notin S_{0}$, we observe on
$\Omega\cap\\{j\in S^{a,b^{*}(a)}\\}$ that

$\displaystyle\frac{3}{4}\gamma_{\min}$
$\displaystyle\leq\bigl{[}P^{a,b^{*}(a),\top}Q^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j,j}-\bigl{\|}\hat{Q}^{a,b^{*}(a)}\\!-\\!Q^{a,b^{*}(a)}\bigr{\|}_{\mathrm{op}}\leq\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j,j}$
$\displaystyle\leq[P^{a,b^{*}(a),\top}Q^{a,b^{*}(a)}P^{a,b^{*}(a)}]_{j,j}+\bigl{\|}\hat{Q}^{a,b^{*}(a)}-Q^{a,b^{*}(a)}\bigr{\|}_{\mathrm{op}}\leq\frac{5}{4}\gamma_{\max},$

and similarly on $\Omega\cap\\{j^{\prime}\in S^{a,b^{*}(a)}\\}$ that
$\bigl{|}\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j^{\prime},j^{\prime}}\bigr{|}\leq\gamma_{\min}/4$.
Recall also that
$\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{{\tilde{j},\tilde{j}}}=0$
for all $\tilde{j}\notin S^{a,b^{*}(a)}$. Combining the above bounds on the
diagonal entries of $\hat{Q}^{a,b^{*}(a)}$ with (15) and (17), we have on
$\Omega$ that

$\displaystyle\mathbb{E}\bigl{(}$
$\displaystyle\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j,j}-\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j^{\prime},j^{\prime}}\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}$
$\displaystyle=\mathbb{E}\bigl{(}\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j,j}\mathbbm{1}_{\\{j\in
S^{a,b^{*}(a)}\\}}$ $\displaystyle\hskip
56.9055pt-\bigl{[}P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}\bigr{]}_{j^{\prime},j^{\prime}}\mathbbm{1}_{\\{j^{\prime}\in
S^{a,b^{*}(a)}\\}}\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}$
$\displaystyle\geq\frac{q_{j}\gamma_{\min}}{2}\geq\frac{\gamma_{\min}}{2p}.$
(18)

Now, let $a,j,j^{\prime}$ be freely varying again. Since
$\hat{w}_{j}=A^{-1}\sum_{a\in[A]}[P^{a,b^{*}(a),\top}\hat{Q}^{a,b^{*}(a)}P^{a,b^{*}(a)}]_{j,j}$,
on $\Omega$ we have for any $j\in S_{0}$ and $j^{\prime}\notin S_{0}$ that
$\mathbb{E}\bigl{(}\hat{w}_{j}-\hat{w}_{j^{\prime}}\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}\geq\gamma_{\min}/(2p)$
from (18). Since $\ell\geq s_{0}$, we have by Hoeffding’s inequality that on
$\Omega$,

$\displaystyle\mathbb{P}\bigl{(}S_{0}\not\subseteq\hat{S}\mid(X_{i},Y_{i})_{i\in[n]}\bigr{)}$
$\displaystyle\leq\mathbb{P}\Bigl{(}\min_{j\in
S_{0}}\hat{w}_{j}\leq\max_{j^{\prime}\notin
S_{0}}\hat{w}_{j^{\prime}}\Bigm{|}(X_{i},Y_{i})_{i\in[n]}\Bigr{)}$
$\displaystyle\leq\sum_{j\in
S_{0}}\mathbb{P}\biggl{\\{}\hat{w}_{j}-\mathbb{E}(\hat{w}_{j}\mid(X_{i},Y_{i})_{i\in[n]})\leq-\frac{\gamma_{\min}}{4p}\biggm{|}(X_{i},Y_{i})_{i\in[n]}\biggr{\\}}$
$\displaystyle\quad+\sum_{j\notin
S_{0}}\mathbb{P}\biggl{\\{}\hat{w}_{j}-\mathbb{E}(\hat{w}_{j}\mid(X_{i},Y_{i})_{i\in[n]})\geq\frac{\gamma_{\min}}{4p}\biggm{|}(X_{i},Y_{i})_{i\in[n]}\biggr{\\}}$
$\displaystyle\leq
p\exp\biggl{\\{}-\frac{A}{2}\biggl{(}\frac{\gamma_{\min}}{4p}\biggr{)}^{2}\biggm{/}\biggl{(}\frac{5\gamma_{\max}}{4}\biggr{)}^{2}\biggr{\\}}\leq
pe^{-A\gamma_{\min}^{2}/(50p^{2}\gamma_{\max}^{2})},$

as desired.

###  5.3 Proof of Theorem 3

The main ingredient of the proof of Theorem 3 is the following proposition,
which controls the rate of convergence of the sample between- and within-class
covariance matrices to their respective population versions in a
classification problem.

######  Proposition 8 (Rate of convergence for LDA).

Suppose that $(Z_{1},Y_{1}),\ldots,(Z_{n},Y_{n})\in\mathbb{R}^{d}\times[K]$
are independent and identically distributed data-label pairs, such that
$\mathbb{P}(Y_{1}=k)=\pi_{k}$ and $Z_{1}\mid
Y_{1}=k\sim\mathcal{N}_{d}(\mu_{k},\Sigma_{\mathrm{w}})$ for $k\in[K]$. Write
$\mu:=\sum_{k\in[K]}\pi_{k}\mu_{k}$ and
$\Sigma_{\mathrm{b}}:=\sum_{k\in[K]}\pi_{k}(\mu_{k}-\mu)(\mu_{k}-\mu)^{\top}$
and let $\hat{\Sigma}_{\mathrm{w}}$ and $\hat{\Sigma}_{\mathrm{b}}$ be
computed as in (3) applied to $(Z_{1},Y_{1}),\ldots,(Z_{n},Y_{n})$. If
$\|\mu_{k}-\mu\|\leq R_{1}$ for all $k\in[K]$ and
$\|\Sigma_{\mathrm{w}}\|_{\mathrm{op}}\leq R_{2}$ for some $R_{1},R_{2}>0$,
then for every $\delta\in(0,1/4]$, we have with probability at least
$1-\delta$ that

$\displaystyle\|\hat{\Sigma}_{\mathrm{b}}-\Sigma_{\mathrm{b}}\|_{\mathrm{op}}$
$\displaystyle\leq\frac{12R_{2}\\{K+\log(8\cdot
9^{d}/\delta)\\}}{n}+\bigl{(}4R_{1}\sqrt{R_{2}}+R_{1}^{2}+R_{1}\bigr{)}\sqrt{\frac{2\log(8\cdot
9^{d}/\delta)}{n}},$
$\displaystyle\|\hat{\Sigma}_{\mathrm{w}}-\Sigma_{\mathrm{w}}\|_{\mathrm{op}}$
$\displaystyle\leq\frac{4R_{2}\\{K+\log(8\cdot
9^{d}/\delta)\\}}{n}+4R_{2}\sqrt{\frac{\log(8\cdot 9^{d}/\delta)}{n}}.$

###### Proof.

We first control the rate of convergence of $\hat{\Sigma}_{\mathrm{b}}$. For
$n_{k}:=\sum_{i=1}^{n}\mathbbm{1}_{\\{Y_{i}=k\\}}$ and
$\tilde{\mu}:=\sum_{k\in[K]}(n_{k}/n)\mu_{k}$, we define
$\tilde{\Sigma}_{\mathrm{b}}^{(1)}:=\sum_{k=1}^{K}(n_{k}/n)(\mu_{k}-\mu)(\mu_{k}-\mu)^{\top}$
and
$\tilde{\Sigma}_{\mathrm{b}}^{(2)}:=\sum_{k=1}^{K}(n_{k}/n)(\mu_{k}-\tilde{\mu})(\mu_{k}-\tilde{\mu})^{\top}$.
We have the following decomposition

$\|\hat{\Sigma}_{\mathrm{b}}-\Sigma_{\mathrm{b}}\|_{\mathrm{op}}\leq\|\tilde{\Sigma}_{\mathrm{b}}^{(1)}-\Sigma_{\mathrm{b}}\|_{\mathrm{op}}+\|\tilde{\Sigma}_{\mathrm{b}}^{(2)}-\tilde{\Sigma}_{\mathrm{b}}^{(1)}\|_{\mathrm{op}}+\|\hat{\Sigma}_{\mathrm{b}}-\tilde{\Sigma}_{\mathrm{b}}^{(2)}\|_{\mathrm{op}}.$
(19)

We will control the three terms on the right-hand side above separately. For
the first term, we note that $\|\mu_{k}-\mu\|\leq R_{1}$ for all $k\in[K]$.
Since $n_{1},\ldots,n_{K}$ are functions of $Y_{1},\ldots,Y_{n}$, we have by
McDiarmid’s inequality (see, e.g. Boucheron, Lugosi and Massart, 2013, Theorem
6.2) that with probability at least $1-\delta/4$,

$\|\tilde{\Sigma}_{\mathrm{b}}^{(1)}-\Sigma_{\mathrm{b}}\|_{\mathrm{op}}\leq
R_{1}^{2}\sum_{k=1}^{K}\biggl{|}\frac{n_{k}}{n}-\pi_{k}\biggr{|}\leq
R_{1}^{2}\sqrt{\frac{2\log(4/\delta)}{n}}.$ (20)

For the second term, we first apply McDiarmid’s inequality again to see that
with probability at least $1-\delta/4$, we have

$\|\tilde{\mu}-\mu\|\leq\sum_{k=1}^{K}\biggl{|}\frac{n_{k}}{n}-\pi_{k}\biggr{|}\|\mu_{k}-\mu\|\leq
R_{1}\sqrt{\frac{2\log(4/\delta)}{n}}.$

Thus, we have with probability at least $1-\delta/4$ that

$\displaystyle\|\tilde{\Sigma}_{\mathrm{b}}^{(2)}-\tilde{\Sigma}_{\mathrm{b}}^{(1)}\|_{\mathrm{op}}$
$\displaystyle\leq
2\sum_{k=1}^{K}\frac{n_{k}}{n}\|(\mu_{k}-\mu)(\mu-\tilde{\mu})^{\top}\|_{\mathrm{op}}+\sum_{k=1}^{K}\frac{n_{k}}{n}\|(\mu-\tilde{\mu})(\mu-\tilde{\mu})^{\top}\|_{\mathrm{op}}$
$\displaystyle\leq 2R_{1}\|\mu-\tilde{\mu}\|+\|\mu-\tilde{\mu}\|^{2}\leq
3R_{1}\|\mu-\tilde{\mu}\|\leq 3R_{1}^{2}\sqrt{\frac{2\log(4/\delta)}{n}}.$
(21)

Finally, for the third term, we write $\hat{\mu}_{k}:=\sum_{i:Y_{i}=k}Z_{i}$
and note that $V_{k}:=n_{k}^{1/2}\hat{\mu}_{k}$ satisfies $V_{k}\mid
Y_{1},\ldots,Y_{n}\sim\mathcal{N}_{d}(n_{k}^{1/2}\mu_{k},\Sigma_{\mathrm{w}})$.
Defining $N:=(n_{1}^{1/2},\ldots,n_{K}^{1/2})^{\top}$ and
$P:=NN^{\top}/n\in\mathbb{R}^{K\times K}$, we may write

$n\hat{\Sigma}_{\mathrm{b}}=\sum_{k=1}^{K}n_{k}\hat{\mu}_{k}\hat{\mu}_{k}^{\top}-n\hat{\mu}\hat{\mu}^{\top}=V^{\top}(I_{K}-P)V,$

where $V:=(V_{1},\ldots,V_{K})^{\top}$, and where
$\hat{\mu}:=n^{-1}\sum_{i=1}^{n}Z_{i}=\sum_{k=1}^{K}(n_{k}/n)\hat{\mu}_{k}$.
By Lemma 17, we deduce that $n\hat{\Sigma}_{\mathrm{b}}$ conditional on
$Y_{1},\ldots,Y_{n}$ has a $d$-dimensional non-central Wishart distribution
with $K-1$ degrees of freedom, covariance matrix $\Sigma_{\mathrm{w}}$ and
non-centrality matrix $n\tilde{\Sigma}_{\mathrm{b}}^{(2)}$, which we denote as

$n\hat{\Sigma}_{\mathrm{b}}\mid
Y_{1},\ldots,Y_{n}\sim\mathcal{W}_{d}(K-1,\Sigma_{\mathrm{w}};n\tilde{\Sigma}_{\mathrm{b}}^{(2)});$

a formal definition is given just before Lemma 17. For any fixed
$u\in\mathbb{S}^{d-1}$, we have by Muirhead (2009, Theorem 10.3.6) that

$u^{\top}\hat{\Sigma}_{\mathrm{b}}u\mid
Y_{1},\ldots,Y_{n}\sim\frac{u^{\top}\Sigma_{\mathrm{w}}u}{n}\chi^{2}_{K-1}\biggl{(}\frac{nu^{\top}\tilde{\Sigma}_{\mathrm{b}}^{(2)}u}{u^{\top}\Sigma_{\mathrm{w}}u}\biggr{)},$

where $\chi^{2}_{r}(\lambda)$ denotes a non-central chi-squared distribution
with $r$ degrees of freedom and non-centrality parameter $\lambda$. By Birgé
(2001, Lemma 8.1), for every $\delta^{\prime}\in(0,1/2]$, we have with
probability at least $1-2\delta^{\prime}$ conditional on $Y_{1},\ldots,Y_{n}$,
that

$\displaystyle|u^{\top}(\hat{\Sigma}_{\mathrm{b}}-\tilde{\Sigma}_{\mathrm{b}}^{(2)})u|$
$\displaystyle\leq\frac{u^{\top}\Sigma_{\mathrm{w}}u}{n}\biggl{\\{}K+2\sqrt{\biggl{(}K+\frac{2nu^{\top}\tilde{\Sigma}_{\mathrm{b}}^{(2)}u}{u^{\top}\Sigma_{\mathrm{w}}u}\biggr{)}\log(1/\delta^{\prime})}+2\log(1/\delta^{\prime})\biggr{\\}}$
$\displaystyle\leq\frac{u^{\top}\Sigma_{\mathrm{w}}u}{n}\bigl{\\{}2K+3\log(1/\delta^{\prime})\bigr{\\}}+\sqrt{\frac{8u^{\top}\Sigma_{\mathrm{w}}uu^{\top}\tilde{\Sigma}_{\mathrm{b}}^{(2)}u\log(1/\delta^{\prime})}{n}}$
$\displaystyle\leq\frac{3\bigl{(}K+\log(1/\delta^{\prime})\bigr{)}\|\Sigma_{\mathrm{w}}\|_{\mathrm{op}}}{n}+\sqrt{\frac{8\|\Sigma_{\mathrm{w}}\|_{\mathrm{op}}\|\tilde{\Sigma}_{\mathrm{b}}^{(2)}\|_{\mathrm{op}}\log(1/\delta^{\prime})}{n}}.$

Let $\mathcal{N}$ be a 1/4-net of the sphere $\mathbb{S}^{d-1}$, which can be
chosen to have cardinality at most $9^{d}$ (Vershynin, 2012, Lemma 5.2).
Hence, through a union bound, and taking $\delta^{\prime}:=\delta/(8\cdot
9^{d})$, we have with probability at least $1-\delta/4$ conditional on
$Y_{1},\ldots,Y_{n}$ that

$\displaystyle\|\hat{\Sigma}_{\mathrm{b}}-\tilde{\Sigma}_{\mathrm{b}}^{(2)}\|_{\mathrm{op}}$
$\displaystyle\leq
2\max_{u\in\mathcal{N}}|u^{\top}(\hat{\Sigma}_{\mathrm{b}}-\tilde{\Sigma}_{\mathrm{b}}^{(2)})u|$
$\displaystyle\leq\frac{6R_{2}\bigl{(}K+\log(8\cdot
9^{d}/\delta)\bigr{)}}{n}+\sqrt{\frac{32R_{1}^{2}R_{2}\log(8\cdot
9^{d}/\delta)}{n}}.$ (22)

Combining (19), (20), (21) and (22), we have that the desired bound on
$\|\hat{\Sigma}_{\mathrm{b}}-\Sigma_{\mathrm{b}}\|_{\mathrm{op}}$ occurs on an
event with probability at least $1-3\delta/4$.

We now turn to control
$\|\hat{\Sigma}_{\mathrm{w}}-\Sigma_{\mathrm{w}}\|_{\mathrm{op}}$. Let
$Q:=\sum_{k\in[K]}(n_{k}^{-1}\mathbbm{1}_{\\{Y_{i}=k,Y_{i^{\prime}}=k\\}})_{i,i^{\prime}=1}^{n}\in\mathbb{R}^{n\times
n}$, so that

$n\hat{\Sigma}_{\mathrm{w}}=\sum_{i=1}^{n}Z_{i}Z_{i}^{\top}-\sum_{k=1}^{K}n_{k}\hat{\mu}_{k}\hat{\mu}_{k}^{\top}=Z^{\top}(I-Q)Z,$

where $Z:=(Z_{1},\ldots,Z_{n})^{\top}$. It therefore follows again by Lemma 17
that $\hat{\Sigma}_{\mathrm{w}}\mid Y_{1},\ldots,Y_{n}\sim
n^{-1}\mathcal{W}_{d}(n-K,\Sigma_{\mathrm{w}})$. Another application of
Muirhead (2009, Theorem 10.3.6) then yields for any $u\in\mathbb{S}^{d-1}$
that

$u^{\top}\hat{\Sigma}_{\mathrm{w}}u\mid
Y_{1},\ldots,Y_{n}\sim\frac{u^{\top}\Sigma_{\mathrm{w}}u}{n}\chi^{2}_{n-K}.$

By Laurent and Massart (2000, Lemma 1), we have with probability at least
$1-2\delta^{\prime}=1-\delta/(4\cdot 9^{d})$ that

$\displaystyle|u^{\top}(\hat{\Sigma}_{\mathrm{w}}-\Sigma_{\mathrm{w}})u|$
$\displaystyle\leq\frac{R_{2}}{n}\bigl{\\{}K+2\sqrt{n\log(1/\delta^{\prime})}+2\log(1/\delta^{\prime})\bigr{\\}}.$

Again, taking a union bound over the $1/4$-net $\mathcal{N}$, we conclude that
with probability at least $1-\delta/4$, we have

$\|\hat{\Sigma}_{\mathrm{w}}-\Sigma_{\mathrm{w}}\|_{\mathrm{op}}\leq
2\max_{u\in\mathcal{N}}|u^{\top}(\hat{\Sigma}_{\mathrm{w}}-\Sigma_{\mathrm{w}})u|\leq\frac{4R_{2}}{n}\Bigl{\\{}K+\sqrt{n\log(8\cdot
9^{d}/\delta)}+\log(8\cdot 9^{d}/\delta)\Bigr{\\}},$

as desired. ∎

###### Proof of Theorem 3.

Define $\delta:=\binom{p}{d}^{-1}\varepsilon$. Since Algorithm 2 is
permutation-equivariant, by a union bound, it suffices to show that for every
$P\in\mathcal{P}_{d}$, with probability at least $1-\delta$, the desired upper
bound holds for
$\bigl{\|}\psi\bigl{(}(PX_{i},Y_{i})_{i\in[n]}\bigr{)}-(P\Sigma_{\mathrm{w}}P^{\top})^{-1}P\Sigma_{\mathrm{b}}P^{\top}\bigr{\|}_{\mathrm{op}}$.
Recall that $n_{k}:=\sum_{i=1}^{n}\mathbbm{1}_{\\{Y_{i}=k\\}}$. Write
$\Sigma_{\mathrm{w},P}:=P\Sigma_{\mathrm{w}}P^{\top}$ and
$\Sigma_{\mathrm{b},P}:=P\Sigma_{\mathrm{b}}P^{\top}$,
$\hat{\mu}_{k,P}:=n_{k}^{-1}\sum_{i:y_{i}=k}PX_{i}$,
$\hat{\mu}_{P}:=n^{-1}\sum_{i=1}^{n}PX_{i}$ and

$\hat{\Sigma}_{\mathrm{w},P}:=\frac{1}{n}\sum_{i=1}^{n}(PX_{i}-\hat{\mu}_{Y_{i},P})(PX_{i}-\hat{\mu}_{Y_{i},P})^{\top}\quad\text{and}\quad\hat{\Sigma}_{\mathrm{b},P}:=\sum_{k=1}^{K}\frac{n_{k}}{n}(\hat{\mu}_{k,P}-\hat{\mu}_{P})(\hat{\mu}_{k,P}-\hat{\mu}_{P})^{\top}.$

Observe that since $n\geq Kd+1$, we have $\max_{k\in[K]}n_{k}\geq d+1$, so
$\hat{\Sigma}_{\mathrm{w},P}$ is positive definite with probability 1. Thus,
by the triangle inequality, for each $P\in\mathcal{P}_{d}$, we have

$\displaystyle\bigl{\|}\psi\bigl{(}(PX_{i},Y_{i})_{i\in[n]}\bigr{)}$
$\displaystyle-(P\Sigma_{\mathrm{w}}P^{\top})^{-1}P\Sigma_{\mathrm{b}}P^{\top}\bigr{\|}_{\mathrm{op}}=\bigl{\|}\hat{\Sigma}_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{w},P}^{-1}\Sigma_{\mathrm{b},P}\bigr{\|}_{\mathrm{op}}$
$\displaystyle\leq\|\hat{\Sigma}_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}\|_{\mathrm{op}}+\|\Sigma_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{w},P}^{-1}\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}.$
(23)

By Proposition 8 and our hypothesis, there is an event $\Omega_{P}$ with
probability at least $1-\delta$, on which

$\displaystyle\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}$
$\displaystyle\leq\frac{4R_{2}\\{K+\log(8\cdot
9^{d}/\delta)\\}}{n}+4R_{2}\sqrt{\frac{\log(8\cdot
9^{d}/\delta)}{n}}\leq\frac{1}{2R_{2}}.$ (24)
$\displaystyle\|\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}$
$\displaystyle\lesssim_{R_{1},R_{2}}\frac{K+d+\log(1/\delta)}{n}+\sqrt{\frac{d+\log(1/\delta)}{n}}\lesssim_{R_{2}}1.$

Thus, for the first term in (5.3), by Weyl’s inequality, on $\Omega_{P}$, we
have

$\displaystyle\|\hat{\Sigma}_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}\|_{\mathrm{op}}$
$\displaystyle\leq\|\Sigma_{\mathrm{w},P}^{-1}-\hat{\Sigma}_{\mathrm{w},P}^{-1}\|_{\mathrm{op}}\|\hat{\Sigma}_{\mathrm{b},P}\|_{\mathrm{op}}$
$\displaystyle\leq\|\Sigma_{\mathrm{w},P}^{-1}\|_{\mathrm{op}}\|\hat{\Sigma}_{\mathrm{w},P}^{-1}\|_{\mathrm{op}}\|\hat{\Sigma}_{\mathrm{b},P}\|_{\mathrm{op}}\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}$
$\displaystyle\leq\frac{\bigl{(}\|\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}+\|\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}\bigr{)}\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}}{\lambda_{\min}(\Sigma_{\mathrm{w},P})\bigl{(}\lambda_{\min}(\Sigma_{\mathrm{w},P})-\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}\bigr{)}}$
$\displaystyle\leq\frac{R_{2}\bigl{(}R_{1}^{2}+\|\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}\bigr{)}\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}}{\bigl{(}1/R_{2}-\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}\bigr{)}}$
$\displaystyle\lesssim_{R_{1},R_{2}}\|\hat{\Sigma}_{\mathrm{w},P}-\Sigma_{\mathrm{w},P}\|_{\mathrm{op}}\lesssim_{R_{2}}\frac{K}{n}+\sqrt{\frac{d+\log(1/\delta)}{n}},$
(25)

where we used (24) in the penultimate inequality. For the second term in
(5.3), we also have on $\Omega_{P}$ that

$\bigl{\|}\Sigma_{\mathrm{w},P}^{-1}\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{w},P}^{-1}\Sigma_{\mathrm{b},P}\bigr{\|}_{\mathrm{op}}\leq\|\Sigma_{\mathrm{w},P}^{-1}\|_{\mathrm{op}}\|\hat{\Sigma}_{\mathrm{b},P}-\Sigma_{\mathrm{b},P}\|_{\mathrm{op}}\lesssim_{R_{1},R_{2}}\frac{K}{n}+\sqrt{\frac{d+\log(1/\delta)}{n}}.$
(26)

The desired result follows by combining (25) and (26), and using the fact that
$\log(1/\delta)\leq d\log(ep/d)+\log(1/\varepsilon)$. ∎

###  5.4 Proofs of Proposition 5 and Theorem 6

In the proof of Proposition 5, we show the convergence of the EM iterates
$\hat{\mu}^{(t)}$ by analyzing their components parallel and orthogonal to
$\mu^{*}$ separately. Writing $\eta:=\mu^{*}/\|\mu^{*}\|$, let
$\alpha_{t}\in\mathbb{R}$, $\beta_{t}\geq 0$ be defined by

$\hat{\mu}^{(t)}=\alpha_{t}\eta+\beta_{t}\xi_{t},$ (27)

where $\xi_{t}\in\mathbb{S}^{d-1}$ is orthogonal to $\eta$. Our proof will
combine several propositions that control $\alpha_{t}$ and $\beta_{t}$ under
different conditions. We begin by laying some groundwork and defining some
quantities that will be used throughout this subsection.

First, it will be convenient to relabel the two classes as $\\{-1,1\\}$
instead of $\\{1,2\\}$. By the rotational symmetry of the problem, we may
assume without loss of generality that
$\mu^{*}=(s,0,\ldots,0)^{\top}\in\mathbb{R}^{d}$ for some $s\geq 0$, and that
the first $n_{\mathrm{L}}$ observations are labeled (i.e., $Y_{i}\neq 0$ for
$i\in[n_{\mathrm{L}}]$). We assume throughout this section that $s\leq r$ and
$r\geq 1$. Let
$\hat{\mu}_{n_{\mathrm{L}}}:=n_{\mathrm{L}}^{-1}\sum_{i=1}^{n_{\mathrm{L}}}Z_{i}Y_{i}$,
with the convention that $\hat{\mu}_{n_{\mathrm{L}}}:=0$ if
$n_{\mathrm{L}}=0$, and define the function
$f_{n_{\mathrm{U}}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ by

$f_{n_{\mathrm{U}}}(v):=\frac{1}{n_{\mathrm{U}}}\sum_{i=n_{\mathrm{L}}+1}^{n}Z_{i}\tanh\langle
Z_{i},v\rangle,$ (28)

with $f_{n_{\mathrm{U}}}:=0$ if $n_{\mathrm{U}}=0$. Throughout, and without
further comment, we assume that $n=n_{\mathrm{L}}+n_{\mathrm{U}}\geq 2$. In
this notation, the EM update (12) can be rewritten, defining the function
$g_{n}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, as

$\hat{\mu}^{(t)}=g_{n}(\hat{\mu}^{(t-1)}):=\gamma\hat{\mu}_{n_{\mathrm{L}}}+(1-\gamma)f_{n_{\mathrm{U}}}(\hat{\mu}^{(t-1)}).$

The corresponding population quantities are

$f(v):=\mathbb{E}Z_{1}\tanh\langle v,Z_{1}\rangle\quad\text{and}\quad
g(v):=\gamma\mu^{*}+(1-\gamma)f(v).$

Writing $\Delta_{n_{\mathrm{U}}}:=f_{n_{\mathrm{U}}}-f$, we have

$g_{n}(v)=g(v)+(1-\gamma)\Delta_{n_{\mathrm{U}}}(v)+\gamma(\hat{\mu}_{n_{\mathrm{L}}}-\mu^{*}).$
(29)

For $\omega,\phi>0$ and $r\geq 1$, we define the following two events that
control the terms in the EM iteration involving the unlabeled and labeled data
respectively:

$\displaystyle\Omega_{1}(\omega)$
$\displaystyle:=\biggl{\\{}\sup_{v\in\mathbb{R}^{d}}\|g_{n}(v)\|\leq
2(r+\sqrt{d})\biggr{\\}}\cap\biggl{\\{}\sup_{\begin{subarray}{c}\|v\|\leq
2(r+\sqrt{d})\\\ v\neq
0\end{subarray}}\frac{\|\Delta_{n_{\mathrm{U}}}(v)\|}{\|v\|}\leq\omega\biggr{\\}}$
(30) $\displaystyle\Omega_{2}(\phi)$
$\displaystyle:=\bigl{\\{}\|\hat{\mu}_{n_{\mathrm{L}}}-\mu^{*}\|\leq\phi\bigr{\\}}.$

######  Proposition 9.

There exists $C_{r}>0$, depending only on $r$, such that for any
$\delta\in(2e^{-n},1]$ and $\omega=C_{r}\sqrt{\frac{d\log
n+\log(1/\delta)}{n_{\mathrm{U}}}}$, we have
$\mathbb{P}\bigl{(}\Omega_{1}(\omega)^{\mathrm{c}}\bigr{)}\leq\delta$.
Moreover, for any $\delta\in(0,1]$ and for
$\phi=\sqrt{\frac{2d+3\log(1/\delta)}{n_{\mathrm{L}}}}$, we have
$\mathbb{P}\bigl{(}\Omega_{2}(\phi)^{\mathrm{c}}\bigr{)}\leq\delta$.

###### Proof.

For any $v\in\mathbb{R}^{d}$,

$\displaystyle\|g_{n}(v)\|=\|(1-\gamma)f_{n_{\mathrm{U}}}(v)+\gamma\hat{\mu}_{n_{\mathrm{L}}}\|$
$\displaystyle\leq(1-\gamma)\cdot\frac{1}{n_{\mathrm{U}}}\sum_{i=n_{\mathrm{L}}+1}^{n}\|Z_{i}\|+\frac{\gamma}{n_{\mathrm{L}}}\sum_{i=1}^{n_{\mathrm{L}}}\|Z_{i}\|$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\|Z_{i}\|\leq\biggl{(}\frac{1}{n}\sum_{i=1}^{n}\|Z_{i}\|^{2}\biggr{)}^{1/2}.$

Since $\sum_{i=1}^{n}\|Z_{i}\|^{2}\sim\chi^{2}_{nd}(ns^{2})$, by Birgé (2001,
Lemma 8.1), we have with probability at least $1-\delta/2$ that

$\displaystyle\sup_{v\in\mathbb{R}^{d}}\|g_{n}(v)\|^{2}$ $\displaystyle\leq
d+s^{2}+2\sqrt{\frac{(d+2s^{2})\log(2/\delta)}{n}}+\frac{2\log(2/\delta)}{n}$
$\displaystyle\leq 2d+3s^{2}+\frac{3\log(2/\delta)}{n}\leq 4(r+\sqrt{d})^{2}.$
(31)

Also, by a very similar argument as in the proof of Wu and Zhou (2022, Theorem
4), we have with probability at least $1-\delta/2$ that

$\sup_{\begin{subarray}{c}\|v\|\leq 2(r+\sqrt{d})\\\ v\neq
0\end{subarray}}\frac{\|\Delta_{n_{\mathrm{U}}}(v)\|}{\|v\|}\leq
C_{r}\sqrt{\frac{d\log n+\log(1/\delta)}{n_{\mathrm{U}}}},$ (32)

for some $C_{r}>0$ depending only on $r$. The first claim follows by combining
(5.4) and (32).

For the second claim, we have $\hat{\mu}_{n_{\mathrm{L}}}\sim
N_{d}(\mu^{*},n_{\mathrm{L}}^{-1}I_{d})$. Hence, by Laurent and Massart (2000,
Lemma 1), we have

$\displaystyle\mathbb{P}\bigl{(}\Omega_{2}(\phi)^{\mathrm{c}}\bigr{)}$
$\displaystyle=\mathbb{P}\bigl{(}n_{\mathrm{L}}\|\hat{\mu}_{n_{\mathrm{L}}}-\mu^{*}\|^{2}>n_{\mathrm{L}}\phi^{2}\bigr{)}$
$\displaystyle\leq\mathbb{P}\bigl{(}n_{\mathrm{L}}\|\hat{\mu}_{n_{\mathrm{L}}}-\mu^{*}\|^{2}>d+2\sqrt{d\log(1/\delta)}+2\log(1/\delta)\bigr{)}\leq\delta,$

as required. ∎

For any $a\in\mathbb{R}$, $b\in[0,\infty)$ and $\xi\in\mathbb{S}^{d-1}$ that
is orthogonal to $\eta$, we define $F(a,b):=\eta^{\top}f(a\eta+b\xi)$ and
$G(a,b):=\|(I_{d}-\eta\eta^{\top})f(a\eta+b\xi)\|$. Note that the distribution
of $Z_{1}$ is orthogonally invariant along the axis $\mu^{*}$; in other words,
if $P\in\mathbb{R}^{d\times d}$ is orthogonal and has $\mu^{*}$ as an
eigenvector with eigenvalue 1, then
$PZ_{1}\stackrel{{\scriptstyle\mathrm{d}}}{{=}}Z_{1}$. It follows that
$f(a\eta+b\xi)$, and hence $F(a,b)$ and $G(a,b)$, do not depend on $\xi$. We
remark that

$f(\alpha_{t}\eta+\beta_{t}\xi_{t})=F(\alpha_{t},\beta_{t})\eta+G(\alpha_{t},\beta_{t})\xi_{t+1}^{\prime}$

for some $\xi_{t+1}^{\prime}\in\mathbb{S}^{d-1}$ that is orthogonal to $\eta$.

Proposition 10 controls the magnitude of the component $\beta_{t}$ of the EM
algorithm iterates that is orthogonal to the signal direction $\eta$. We
define $\zeta:=\omega\gamma^{-1/2}\wedge\omega^{1/2}$.

######  Proposition 10.

Assume that $\phi\gamma^{1/2}\leq\omega\leq\min\\{1/12,1/(r+3)\\}$ and that
$\|\hat{\mu}^{(0)}\|\leq r+3$. On the event
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have

$\limsup_{t\to\infty}\beta_{t}\leq 60(\zeta\vee r\omega).$

Moreover, on the same event, if $\beta_{t_{0}}\leq 60(\zeta\vee r\omega)$ for
some $t_{0}\in\mathbb{N}_{0}$, then $\beta_{t}\leq 60(\zeta\vee r\omega)$ for
all $t\geq t_{0}$.

###### Proof.

We first claim that on the event $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we
have $\|\hat{\mu}^{(t)}\|\leq r+3$ for all $t\in\mathbb{N}_{0}$. The case
$t=0$ is true by the assumption on the initializer $\hat{\mu}^{(0)}$, and if
the claim holds for $t\in\mathbb{N}_{0}$, then since $2(r+\sqrt{d})\geq
2(r+1)\geq r+3$, we have on $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$ that

$\displaystyle\|\hat{\mu}^{(t+1)}\|$
$\displaystyle\leq(1-\gamma)\bigl{\\{}|F(\alpha_{t},\beta_{t})|+G(\alpha_{t},\beta_{t})+\|\Delta_{n_{\mathrm{U}}}(\hat{\mu}^{(t)})\|\bigr{\\}}+\gamma(s+\phi)$
$\displaystyle\leq s+2\sqrt{2/\pi}+\omega\|\hat{\mu}^{(t)}\|+\gamma\phi\leq
r+2+\frac{\|\hat{\mu}^{(t)}\|}{r+3}\leq r+3,$

where the second inequality uses Wu and Zhou (2022, Lemma 5(5)). Moreover,
from (29), we have on $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$ that for
$t\in\mathbb{N}$,

$\displaystyle\beta_{t+1}$
$\displaystyle=\bigl{\|}(I_{d}-\eta\eta^{\top})\bigl{\\{}(1-\gamma)\bigl{(}f(\hat{\mu}^{(t)})+\Delta_{n_{\mathrm{U}}}(\hat{\mu}^{(t)})\bigr{)}+\gamma\hat{\mu}_{n_{\mathrm{L}}}\bigr{\\}}\bigr{\|}$
$\displaystyle\leq(1-\gamma)\bigl{\\{}G(\alpha_{t},\beta_{t})+\omega(|\alpha_{t}|+\beta_{t})\bigr{\\}}+\gamma\phi$
$\displaystyle\leq\beta_{t}(1-\gamma)\biggl{\\{}1+\omega-\frac{(\alpha_{t}^{2}+\beta_{t}^{2})\wedge
1}{6}\biggr{\\}}+\gamma\phi+\omega|\alpha_{t}|,$ (33)

where the final bound uses Wu and Zhou (2022, Lemma 5(8)). If
$\alpha_{t}^{2}+\beta_{t}^{2}>1$ or $\gamma>1/2$, then using the fact that
$\omega\leq 1/12$, we have from (33) that

$\beta_{t+1}\leq\frac{11}{12}\beta_{t}+\gamma\phi+(r+3)\omega\leq\frac{11}{12}\beta_{t}+(r+4)\omega.$
(34)

On the other hand, if $\alpha_{t}^{2}+\beta_{t}^{2}\leq 1$ and $\gamma\leq
1/2$, then

$\displaystyle\beta_{t+1}$
$\displaystyle\leq\beta_{t}\biggl{(}1+\omega-\gamma-\frac{\alpha_{t}^{2}+\beta_{t}^{2}}{12}\biggr{)}+\gamma\phi+\omega|\alpha_{t}|.$
(35)

Note that the right-hand side of (34) is increasing in $\beta_{t}$ and the
right-hand side of (35) is increasing in $\beta_{t}$ for
$\alpha_{t}^{2}+\beta_{t}^{2}\leq 1$ and $\gamma\leq 1/2$. Combining (34) and
(35), denoting $\beta_{\infty}:=\limsup_{t\to\infty}\beta_{t}$ and using the
fact that
$0\leq\frac{3}{\beta_{\infty}}(\omega-|\alpha_{t}|\beta_{\infty}/6)^{2}=3\omega^{2}/\beta_{\infty}-\omega|\alpha_{t}|+\alpha_{t}^{2}\beta_{\infty}/12$,
we have

$\beta_{\infty}\leq\max\biggl{\\{}\frac{11}{12}\beta_{\infty}+(r+4)\omega,\,\beta_{\infty}\biggl{(}1+\omega-\gamma-\frac{\beta_{\infty}^{2}}{12}\biggr{)}+\gamma\phi+\frac{3\omega^{2}}{\beta_{\infty}}\biggr{\\}}.$
(36)

From the first term in the maximum in (36), we obtain

$\beta_{\infty}\leq(1r+38)\omega\leq 60r\omega.$ (37)

From the second term in the maximum in (36), we obtain

$\beta_{\infty}\biggl{(}\gamma-\omega+\frac{\beta_{\infty}^{2}}{12}\biggr{)}\leq\gamma\phi+\frac{3\omega^{2}}{\beta_{\infty}}.$
(38)

If $\gamma<2\omega$, then from (38),

$\beta_{\infty}\leq 5\omega^{1/2}\leq 5\sqrt{2}\zeta,$ (39)

since otherwise we would have that the left-hand side would be at least
$(-5+125/12)\omega^{3/2}$ and the right-hand side would at most
$(\sqrt{2}+3/5)\omega^{3/2}$, contradicting the inequality. On the other hand,
if $\gamma\geq 2\omega$, then we derive from (38) that

$\beta_{\infty}^{2}-2\phi\beta_{\infty}-\frac{6\omega^{2}}{\gamma}\leq 0.$

Solving this inequality, we find that

$\beta_{\infty}\leq\phi+\sqrt{\phi^{2}+6\omega^{2}/\gamma}\leq
4\omega\gamma^{-1/2}=4\zeta.$ (40)

The first claim of the proposition follows by combining (36), (37), (39) and
(40). We now prove the second claim by induction on $t$. The base case
$t=t_{0}$ is true by assumption, so we assume that $\beta_{t}\leq 60(\zeta\vee
r\omega)$ for some $t\geq t_{0}$. Again we consider two cases. If
$\beta_{t}\leq 36(\zeta\vee r\omega)$, then from (33) and using that
$|\alpha_{t}|\leq\|\hat{\mu}^{(t)}\|\leq r+3$ for $t\geq 2$,

$\beta_{t+1}\leq\beta_{t}(1+\omega)+\gamma\phi+\omega|\alpha_{t}|\leq
44(\zeta\vee r\omega),$

as desired. On the other hand, if $\beta_{t}>36(\zeta\vee r\omega)$, then
combining (34) and (35), we obtain that

$\displaystyle\beta_{t+1}\leq\max\biggl{\\{}60(\zeta\vee
r\omega),\,\beta_{t}\biggl{(}1+\omega-\gamma-\frac{\beta_{t}^{2}}{12}\biggr{)}+\gamma\phi+\frac{3\omega^{2}}{\beta_{t}}\biggr{\\}}.$

It suffices to show the second term in the maximum is no larger than
$\beta_{t}$. To this end, if $\gamma\leq 2\omega$, then
$\zeta^{2}\leq\omega\leq 2\zeta^{2}$, and so

$\displaystyle\beta_{t}\biggl{(}1$
$\displaystyle+\omega-\gamma-\frac{\beta_{t}^{2}}{12}\biggr{)}+\gamma\phi+\frac{3\omega^{2}}{\beta_{t}}\leq\beta_{t}+2\zeta^{2}\beta_{t}-\frac{\beta_{t}^{3}}{12}+\gamma^{1/2}\omega+\frac{3\omega^{2}}{\beta_{t}}$
$\displaystyle\leq\beta_{t}+120(\zeta\vee
r\omega)\zeta^{2}-\frac{36^{3}(\zeta\vee
r\omega)^{3}}{12}+4\zeta^{3}+\frac{\zeta^{3}}{3}\leq\beta_{t}.$

On the other hand, if $\gamma>2\omega$, then
$\zeta=\omega\gamma^{-1/2}\geq\phi$, and so

$\displaystyle\beta_{t}\biggl{(}1+\omega-\gamma-\frac{\beta_{t}^{2}}{12}\biggr{)}+\gamma\phi+\frac{3\omega^{2}}{\beta_{t}}$
$\displaystyle\leq\beta_{t}-\frac{\gamma\beta_{t}}{2}+\gamma\phi+\frac{3\omega^{2}}{\beta_{t}}$
$\displaystyle\leq\beta_{t}-18\gamma\zeta+\gamma\zeta+\frac{\omega^{2}}{12\zeta}\leq\beta_{t},$

as desired, which completes the induction. ∎

The following result bounds the magnitude of the signal component,
$\alpha_{t}$, of the EM iterates.

######  Proposition 11.

Assume that $\phi\gamma^{1/2}\leq\omega\leq\min\\{1/12,1/(r+3)\\}$ and that
$\|\hat{\mu}^{(0)}\|\leq r+3$. Then there exists $C_{r}>0$, depending only on
$r$, such that on the event $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have

$\limsup_{t\to\infty}|\alpha_{t}|\leq C_{r}(\zeta\vee s).$

###### Proof.

By definition of $\alpha_{t+1}$ and (29), we have for every
$t\in\mathbb{N}_{0}$ that

$\alpha_{t+1}=\eta^{\top}\bigl{\\{}(1-\gamma)\bigl{(}f(\hat{\mu}^{(t)})+\Delta_{n_{\mathrm{U}}}(\hat{\mu}^{(t)})\bigr{)}+\gamma\hat{\mu}_{n_{\mathrm{L}}}\bigr{\\}}.$

Thus, by the first claim in the proof of Proposition 10, we have on the event
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$ that

$\bigl{|}\alpha_{t+1}-(1-\gamma)F(\alpha_{t},\beta_{t})-\gamma
s\bigr{|}\leq(1-\gamma)\omega(|\alpha_{t}|+\beta_{t})+\gamma\phi$ (41)

for every $t\in\mathbb{N}_{0}$. From Wu and Zhou (2022, Lemma 5(1) and Lemma
5(7)), $\alpha\mapsto F(\alpha,\beta)$ is an increasing and odd function
satisfying $|F(\alpha,\beta)-F(\alpha,0)|\leq(1+s^{2})|\alpha|\beta^{2}$ for
every $\alpha,\beta\in\mathbb{R}$. Hence, by (41), we have on
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$ that

$|\alpha_{t+1}|\leq(1-\gamma)\Bigl{\\{}F(|\alpha_{t}|,0)+(1+s^{2})|\alpha_{t}|\beta_{t}^{2}+\omega(|\alpha_{t}|+\beta_{t})\Bigr{\\}}+\gamma(s+\phi).$
(42)

Note the right-hand side of (42) is increasing in $|\alpha_{t}|$. Define
$\alpha_{\infty}:=\limsup_{t\to\infty}|\alpha_{t}|$, so that
$\alpha_{\infty}\leq r+3$ on $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, again
by the first claim in the proof of Proposition 10. We may also assume that
$\alpha_{\infty}>s$, because otherwise the result is clear. Since
$\alpha\mapsto F(\alpha,0)$ is continuous, we have from (42) that on
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$,

$\alpha_{\infty}\leq(1-\gamma)\bigl{\\{}F(\alpha_{\infty},0)+(1+r^{2})\alpha_{\infty}\beta_{\infty}^{2}+\omega(\alpha_{\infty}+\beta_{\infty})\bigr{\\}}+\gamma(s+\phi),$
(43)

where we recall that $\beta_{\infty}:=\limsup_{t\rightarrow\infty}\beta_{t}$.
Define $q:[0,\infty)\rightarrow\mathbb{R}$ by

$q(\alpha):=\left\\{\begin{array}[]{ll}F(\alpha,0)/\alpha&\mbox{if $\alpha\neq
0$}\\\ 1+s^{2}&\mbox{if $\alpha=0$.}\end{array}\right.$ (44)

By Lemma 20, we have $q(s)=1$ (which confirms that $\mu_{*}$ is a fixed point
of the population EM iteration), and that $q^{\prime}(\alpha)\leq-c_{r}\alpha$
for all $\alpha\in(0,r]$, where $c_{r}\in(0,1]$ depends only on $r$. Thus,
dividing both sides of (43) by $\alpha_{\infty}$, we have

$\displaystyle 1$
$\displaystyle\leq(1-\gamma)\biggl{\\{}q(s)+\int_{s}^{\alpha_{\infty}}q^{\prime}(\alpha)\,d\alpha+(1+r^{2})\beta_{\infty}^{2}+\omega\biggl{(}1+\frac{\beta_{\infty}}{\alpha_{\infty}}\biggr{)}\biggr{\\}}+\frac{\gamma(s+\phi)}{\alpha_{\infty}}$
$\displaystyle\leq(1-\gamma)\biggl{\\{}1-\frac{c_{r}}{2}(\alpha_{\infty}^{2}-s^{2})+(1+r^{2})\beta_{\infty}^{2}+\omega\biggl{(}1+\frac{\beta_{\infty}}{\alpha_{\infty}}\biggr{)}\biggr{\\}}+\frac{\gamma(s+\phi)}{\alpha_{\infty}}.$
(45)

Now $\beta_{\infty}\leq 60(1+r)\zeta$ by Proposition 10. We now claim that
$\alpha_{\infty}\leq 4s+120c_{r}^{-1/2}(1+r)^{2}\zeta$. Indeed, assuming the
contrary, we would have
$c_{r}\alpha_{\infty}^{2}/4>c_{r}s^{2}/2+(1+r^{2})\beta_{\infty}^{2}$ and
$\beta_{\infty}/\alpha_{\infty}<1$. Hence from (5.4), we have

$\displaystyle 1$
$\displaystyle\leq(1-\gamma)\biggl{(}1-\frac{c_{r}\alpha_{\infty}^{2}}{4}\biggr{)}+\frac{\gamma
s}{\alpha_{\infty}}+\biggl{(}2+\frac{\gamma^{1/2}}{\alpha_{\infty}}\biggr{)}\omega.$

We consider two cases. First, if $\gamma\leq 4\omega$, then
$2\zeta\geq\omega^{1/2}\geq\gamma^{1/2}/2$ and hence

$\displaystyle(1-\gamma)\biggl{(}1-\frac{c_{r}\alpha_{\infty}^{2}}{4}\biggr{)}+\frac{\gamma
s}{\alpha_{\infty}}+\biggl{(}2+\frac{\gamma^{1/2}}{\alpha_{\infty}}\biggr{)}\omega$
$\displaystyle\leq\max\biggl{\\{}1-\frac{c_{r}\alpha_{\infty}^{2}}{4},0\biggr{\\}}+\frac{\gamma}{4}+3\omega\leq\max\biggl{\\{}1-\zeta^{2},\frac{\gamma}{4}+3\omega\biggr{\\}}<1,$

a contradiction. Second, if $\gamma>4\omega$, then $\zeta=\omega\gamma^{-1/2}$
and

$(1-\gamma)\biggl{(}1-\frac{c_{r}\alpha_{\infty}^{2}}{4}\biggr{)}+\frac{\gamma
s}{\alpha_{\infty}}+\biggl{(}2+\frac{\gamma^{1/2}}{\alpha_{\infty}}\biggr{)}\omega\leq
1-\gamma+\frac{\gamma}{4}+2\omega+\frac{\gamma\zeta}{\alpha_{\infty}}<1,$

again a contradiction. This establishes the claimed upper bound on
$\alpha_{\infty}$. ∎

Recall the definitions $L(\mu,\mu^{*})=\|\mu-\mu^{*}\|\wedge\|\mu+\mu^{*}\|$.
Our next result shows that if $\alpha_{t}$ ever becomes sufficiently large,
then improved bounds can be derived on the limiting behaviour of $\alpha_{t}$,
$\beta_{t}$ and $L(\hat{\mu}^{(t)},\mu^{*})$.

######  Proposition 12.

Assume that $\gamma\in[0,1/2)$. Given any $c>0$, there exists $C,c_{1}>0$,
depending only on $r$ and $c$, such that if $|\alpha_{t_{0}}|\geq cs$,
$\beta_{t_{0}}\leq 60(\zeta\vee r\omega)$ for some iteration $t_{0}$, and
$\phi\gamma^{1/2}\leq\omega\leq c_{1}$ and $s\geq C\zeta$, then on the event
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have

$\displaystyle\limsup_{t\to\infty}|\alpha_{t}-s|$
$\displaystyle\lesssim_{r,c}\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}},$
(46) $\displaystyle\limsup_{t\to\infty}\beta_{t}$
$\displaystyle\lesssim_{r,c}\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}},$
(47) $\displaystyle\limsup_{t\to\infty}L(\hat{\mu}^{(t)},\mu^{*})$
$\displaystyle\lesssim_{r,c}\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}.$
(48)

###### Proof.

By flipping the sign of $\mu^{*}$ if necessary, we may assume without loss of
generality that $\alpha_{t_{0}}\geq 0$ and that
$c_{1}\leq\min\\{1/12,1/(r+3)\\}$. From (41) and the argument immediately
below it, we have

$\displaystyle(1-\gamma)\Bigl{\\{}F(\alpha_{t},0)$
$\displaystyle-(1+s^{2})\alpha_{t}\beta_{t}^{2}-\omega(\alpha_{t}+\beta_{t})\Bigr{\\}}+\gamma(s-\phi)\leq\alpha_{t+1}$
$\displaystyle\leq(1-\gamma)\Bigl{\\{}F(\alpha_{t},0)+(1+s^{2})\alpha_{t}\beta_{t}^{2}+\omega(\alpha_{t}+\beta_{t})\Bigr{\\}}+\gamma(s+\phi).$
(49)

For any $t$ such that $\alpha_{t}\geq cs$, since $\beta_{t}\leq
60r\omega^{1/2}$ by Proposition 10, we have that

$(1+s^{2})\alpha_{t}\beta_{t}^{2}+\omega(\alpha_{t}+\beta_{t})\leq\biggl{\\{}(1+r^{2})60^{2}r^{2}\omega+\omega\biggl{(}1+\frac{60s\zeta}{cs}\biggr{)}\biggr{\\}}\alpha_{t}\leq
c^{\prime\prime}\omega\alpha$ (50)

where $c^{\prime\prime}:=60^{2}(1+r^{2})r^{2}+(1+60rc^{-1}C^{-1})$. Moreover,
if $\alpha_{t}\geq cs$, then

$\gamma\phi\leq\omega\gamma^{1/2}\leq\begin{cases}\gamma\zeta\leq C^{-1}\gamma
s&\text{if $\gamma\geq\omega$}\\\ \omega^{3/2}=\omega\zeta\leq
c^{-1}C^{-1}\omega\alpha_{t}&\text{otherwise}.\end{cases}$ (51)

Let $c^{\prime}:=c^{\prime\prime}+2c^{-1}C^{-1}$ and define functions
$H,L:[0,\infty)\to\mathbb{R}$ by

$\displaystyle H(\alpha)$
$\displaystyle:=(1-\gamma)\bigl{\\{}F(\alpha,0)+c^{\prime}\omega\alpha\bigr{\\}}+\gamma(s+\phi),$
$\displaystyle L(\alpha)$
$\displaystyle:=(1-\gamma)\bigl{\\{}F(\alpha,0)-c^{\prime}\omega\alpha\bigr{\\}}+\gamma\max(s-\phi,s/2),$

From (5.4), (50) and (51), we obtain that for $\alpha_{t}\geq cs$ and $C\geq
2$,

$L(\alpha_{t})\leq\alpha_{t+1}\leq H(\alpha_{t}).$ (52)

Define auxiliary sequences $(\alpha_{t}^{+})_{t\geq t_{0}}$ and
$(\alpha_{t}^{-})_{t\geq 0}$ by
$\alpha_{t_{0}}^{+}:=\alpha_{t_{0}}=:\alpha_{t_{0}}^{-}$ and for $t\geq
t_{0}$,

$\alpha_{t+1}^{+}:=H(\alpha_{t}^{+})\quad\text{and}\quad\alpha_{t+1}^{-}:=L(\alpha_{t}^{-}).$

We first derive some properties of the two recursion maps $H$ and $L$. For the
former, we have by Wu and Zhou (2022, Lemma 3) that $F$, and hence $H$, is
increasing and concave on $[0,\infty)$ with $H(0)>0$ when $\gamma>0$ and
$H^{\prime}(0)>\partial_{1}F(0,0)>1$ when $\gamma=0$. Moreover, since $F$ is
bounded, we can choose $c_{1}>0$, depending only on $r$ and $c$, such that
$\lim_{\alpha\to\infty}H^{\prime}(\alpha)=(1-\gamma)c^{\prime}\omega\leq(1-\gamma)c^{\prime}c_{1}<1/2$.
On the other hand, we have $L(0)>0$ when $\gamma>0$. When $\gamma=0$, we have
$\omega^{1/2}=\zeta\leq s/C$, which means that after increasing $C\equiv
C(r,c)>0$ if necessary, $L^{\prime}(0)=\partial_{1}F(0,0)-c^{\prime}\omega\geq
1+s^{2}-c^{\prime}s^{2}/C^{2}>1$. By Wu and Zhou (2022, Lemma 3),
$\alpha\mapsto F(\alpha,0)$ is differentiable, increasing and concave for
$\alpha\in[0,\infty)$. Reducing $c_{1}\equiv c_{1}(r,c)>0$ if necessary to
ensure that $c_{1}\leq\partial_{1}F(r+3,0)/c^{\prime}$, we have for
$\alpha\in[0,r+3]$ that

$L^{\prime}(\alpha)=\partial_{1}F(\alpha,0)-c^{\prime}\omega\geq\partial_{1}F(r+3,0)-c^{\prime}c_{1}\geq
0.$

In other words, $L$ is increasing on $[0,r+3]$, and moreover, similarly to
$H$, it is also concave on this interval. Finally, we claim that for
$\tilde{c}:=\min\bigl{\\{}c,32(3+r^{4})/3\\}^{-1/2}\bigr{\\}}$, and
$\tilde{\alpha}:=\tilde{c}s$, we have $L(\tilde{\alpha})\geq\tilde{\alpha}$.
To verify this, we note by Wu and Zhou (2022, Lemma 3(3) and Equation (98)),
we have

$L(\tilde{\alpha})\geq(1-\gamma)\tilde{\alpha}\biggl{\\{}1+s^{2}-\frac{8}{3}(3+r^{4})\tilde{\alpha}^{2}-c^{\prime}\omega\biggr{\\}}+2\gamma\tilde{\alpha}\geq(1-\gamma)\tilde{\alpha}\biggl{(}1+\frac{3s^{2}}{4}-c^{\prime}\omega\biggr{)}+2\gamma\tilde{\alpha}.$
(53)

To control the right-hand side of (53), if $\gamma\leq c^{\prime}\omega$, we
have
$\zeta=\omega^{1/2}\wedge\omega\gamma^{-1/2}\geq(\omega/c^{\prime})^{1/2}$.
Hence, from the condition $s\geq C\zeta$, if we choose $C>2c^{\prime}$ (which
is possible because $c^{\prime}$ is a decreasing function of $C$), we have
$c^{\prime}\omega\leq(c^{\prime}s/C)^{2}\leq s^{2}/4$ and consequently the
right-hand side of (53) is at least $\tilde{\alpha}$. If
$\gamma>c^{\prime}\omega$, then we have
$L(\tilde{\alpha})\geq(1-\gamma)^{2}\tilde{\alpha}+2\gamma\tilde{\alpha}\geq\tilde{\alpha}$
as desired. This establishes the claim.

We now show by induction that for all $t\geq t_{0}$,

$\tilde{\alpha}\leq\alpha_{t}^{-}\leq\alpha_{t}\leq\alpha_{t}^{+}.$ (54)

The base case is clear by the definition of $\alpha_{t_{0}}^{-}$ and
$\alpha_{t_{0}}^{+}$ above. Now suppose that (54) holds for some iteration
$t\geq 0$, so in particular, (52) applies.

Using the monotonicity of $H$ and (52), we have $\alpha_{t+1}\leq
H(\alpha_{t})\leq H(\alpha_{t}^{+})=\alpha_{t+1}^{+}$. Observe that
$\alpha_{t}\leq\|\hat{\mu}^{(t)}\|\leq r+3$ by the proof of Proposition 10.
Using the monotonicity of $L$ on $[0,r+3]$, we find that $\alpha_{t+1}\geq
L(\alpha_{t})\geq L(\alpha_{t}^{-})=\alpha_{t+1}^{-}$. Moreover,

$\alpha_{t+1}^{-}=L(\alpha_{t}^{-})\geq L(\tilde{\alpha})\geq\tilde{\alpha},$

which completes the induction.

To prove (46), we will analyze the sequences $(\alpha_{t}^{+})_{t\geq t_{0}}$
and $(\alpha_{t}^{-})_{t\geq t_{0}}$, which sandwich $(\alpha_{t})_{t\geq
t_{0}}$. We start by considering the behaviour of $(\alpha_{t}^{+})_{t\geq
t_{0}}$. The properties of $H$ derived above mean that we can apply Lemma 19
to obtain that $\alpha_{t}^{+}$ converges to a limit, denoted $\alpha^{+}$,
satisfying $\alpha^{+}=H(\alpha^{+})$. By Lemma 20, we have $F(s,0)=s$ and so
$H(s)=(1-\gamma)(s+c^{\prime}\omega s)+\gamma(s+\phi)>s$. Hence from Lemma 19
again, we have $\alpha^{+}>s$. On the other hand, since
$F(\alpha,0)\leq\mathbb{E}|Z_{1,1}|\leq(\mathbb{E}Z_{1,1}^{2})^{1/2}\leq(1+s^{2})^{1/2}\leq
1+r$, we have $\alpha^{+}=H(\alpha^{+})\leq(1-\gamma)(1+r)+\alpha^{+}/2+\gamma
r+1\leq r+2+\alpha^{+}/2$, so $\alpha^{+}\leq 2r+4$. Recalling the definition
of $q$ from (44), by Lemma 20 again, we have

$q(\alpha^{+})=q(s)+\int_{s}^{\alpha^{+}}q^{\prime}(\alpha)\,d\alpha\leq
1-c_{2}\bigl{(}(\alpha^{+})^{2}-s^{2}\bigr{)}$

for some $c_{2}>0$ depending only on $r$. Consequently,

$\displaystyle\alpha^{+}=H(\alpha^{+})$
$\displaystyle=(1-\gamma)\alpha^{+}\bigl{\\{}q(\alpha^{+})+c^{\prime}\omega\bigr{\\}}+\gamma(s+\phi)$
$\displaystyle\leq(1-\gamma)\alpha^{+}\bigl{\\{}1-c_{2}\bigl{(}(\alpha^{+})^{2}-s^{2}\bigr{)}+c^{\prime}\omega\bigr{\\}}+\gamma(s+\omega\gamma^{-1/2}),$

so

$(\alpha^{+})^{2}-s^{2}\leq\frac{c^{\prime}\omega}{c_{2}}-\frac{\gamma}{(1-\gamma)c_{2}}\frac{\alpha^{+}-s-\omega\gamma^{-1/2}}{\alpha^{+}}.$
(55)

We now prove that

$\alpha^{+}-s\lesssim_{r,c}\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}$
(56)

by considering two cases. If $\alpha^{+}\leq 2s$, then from (55), we have

$\alpha^{+}-s\leq\frac{(1-\gamma)c^{\prime}\omega+\gamma^{1/2}\omega/\alpha^{+}}{(1-\gamma)c_{2}(\alpha^{+}+s)+\gamma/\alpha^{+}}\lesssim_{r,c}\frac{\omega}{s}\cdot\frac{1+\gamma^{1/2}/s}{1+\gamma/s^{2}}\lesssim\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}.$

On the other hand, if $\alpha^{+}>2s$, then we have from (55) again that

$\frac{3(\alpha^{+})^{2}}{4}+\frac{\gamma}{2c_{2}}\leq(\alpha^{+})^{2}-s^{2}+\frac{\gamma(\alpha^{+}-s)}{(1-\gamma)c_{2}\alpha^{+}}\leq\frac{(c^{\prime}+2\gamma^{1/2}/\alpha^{+})\omega}{c_{2}}.$
(57)

In particular, $\gamma/2\leq(c^{\prime}+\gamma^{1/2}/s)\omega\leq
c^{\prime}\omega+C^{-1}(\gamma^{1/2}\omega^{1/2}\vee\gamma)$, so
$\gamma\lesssim_{r,c}\omega$ and
$\alpha^{+}\gtrsim_{r,c}\zeta\gtrsim_{r,c}\gamma^{1/2}$. Consequently from
(57),

$\alpha^{+}-s\leq\alpha^{+}\lesssim_{r,c}\omega^{1/2}\lesssim_{r,c}\frac{\omega}{\alpha^{+}}\wedge\frac{\omega}{\gamma^{1/2}}\lesssim\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}},$

which establishes (56). We now consider $(\alpha_{t}^{-})_{t\geq t_{0}}$.
Define $\tilde{L}:[0,\infty)\to[0,\infty)$ by
$\tilde{L}(\alpha):=L\bigl{(}\alpha\wedge(r+3)\bigr{)}$. Since
$\alpha_{t}^{-}\leq\alpha_{t}\leq r+3$ for all $t\geq 0$, we have
$\alpha_{t+1}^{-}=\tilde{L}(\alpha_{t}^{-})$ for all $t\geq t_{0}$. From the
properties of $L$ derived above, we see that $\tilde{L}$ satisfies the
conditions of Lemma 19, and hence $\alpha_{t}^{-}$ converges to a limit,
denoted $\alpha^{-}$, satisfying
$\alpha^{-}=\tilde{L}(\alpha^{-})=L(\alpha^{-})$. By Lemma 20, $F(s,0)=s$, so
we have $\tilde{L}(s)=L(s)\leq(1-\gamma)(s-c^{\prime}\omega s)+\gamma s<s$, so
by Lemma 19, we must have $\alpha^{-}<s$. By Lemma 20 again, we have

$q(\alpha^{-})=q(s)-\int^{s}_{\alpha^{-}}q^{\prime}(\alpha)\,d\alpha\geq
1+c_{2}^{\prime}\bigl{(}s^{2}-(\alpha^{-})^{2}\bigr{)},$

where $c_{2}^{\prime}>0$ depends only on $r$. Consequently, we have

$\alpha^{-}=L(\alpha^{-})\geq(1-\gamma)\alpha^{-}\bigl{\\{}1+c_{2}^{\prime}\bigl{(}s^{2}-(\alpha^{-})^{2}\bigr{)}-c^{\prime}\omega\bigr{\\}}+\gamma(s-\omega\gamma^{-1/2}),$

which after rearranging and using the fact that
$\alpha^{-}\geq\tilde{\alpha}\gtrsim_{r}s$ leads to

$s-\alpha^{-}\leq\frac{(1-\gamma)c^{\prime}\omega+\gamma^{1/2}\omega/\alpha^{-}}{(1-\gamma)c_{2}^{\prime}(s+\alpha^{-})+\gamma/\alpha^{-}}\lesssim_{r,c}\frac{\omega}{s}\cdot\frac{1+\gamma^{1/2}/s}{1+\gamma/s^{2}}\lesssim\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}.$
(58)

Combining (54), (56) and (58), we have established (46).

We now turn to prove (47). By increasing $C$ if necessary, we have for all
sufficiently large $t$ that $s/2\leq\alpha_{t}\leq 2s$. Consequently, we have
by (33) that for all large $t$,

$\beta_{t+1}\leq\beta_{t}(1-\gamma)\biggl{(}1+\omega-\frac{s^{2}/4\wedge
1}{6}\biggr{)}+\gamma\phi+2s\omega\leq\beta_{t}(1-\gamma)(1+\omega-
c_{3}s^{2})+(\gamma^{1/2}+2s)\omega,$ (59)

with $c_{3}:=1/(6r^{2}+24)$. Denote
$\beta_{\infty}:=\limsup_{t\to\infty}\beta_{t}$. If $\gamma\leq 2\omega$, then
$\zeta\in[(\omega/2)^{1/2},\omega^{1/2}]$ so $\omega\leq 2s^{2}/C^{2}\leq
c_{3}s^{2}/2$ for $C$ sufficiently large. Hence, from (59), and the fact that
$\gamma^{1/2}\leq 2\zeta\leq s$,

$\beta_{\infty}\leq\frac{(\gamma^{1/2}+2s)\omega}{c_{3}s^{2}/2}\lesssim_{r}\frac{\omega}{s}=\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}.$

On the other hand, if $\gamma>2\omega$, then $(1-\gamma)(1+\omega-
c_{3}s^{2})\leq 1-\gamma/2-c_{3}s^{2}/2$ and from (59), we obtain

$\beta_{\infty}\leq\frac{(\gamma^{1/2}+2s)\omega}{(\gamma+c_{3}s^{2})/2}\lesssim_{r}\frac{\omega}{s}\cdot\frac{1+\gamma^{1/2}/s}{1+\gamma/s^{2}}\lesssim\frac{\omega}{s}\wedge\frac{\omega}{\gamma^{1/2}}.$

Combining the above two cases establishes (47).

Finally, recalling the decomposition of $\hat{\mu}^{(t)}$ in (27), we see that
(48) follows immediately from (46) and (47). ∎

Next, we show that provided the initialization is not too uncorrelated with
the true parameter, $|\alpha_{t}|$ reaches a level that makes Proposition 12
applicable after a sufficient number of iterations.

######  Proposition 13.

Assume that $n\geq 3$, that
$\phi\gamma^{1/2}\leq\omega\leq\min\\{1/12,1/(r+3)\\}$ and that
$\gamma\in[0,1/2)$. Suppose that $\hat{\mu}^{(0)}$ is chosen such that
$c^{\prime}(\zeta\vee r\omega)\leq\|\hat{\mu}^{(0)}\|\leq 60(\zeta\vee
r\omega)$ for some $c^{\prime}\in(0,1)$ and that
$|\langle\hat{\mu}^{(0)}/\|\hat{\mu}^{(0)}\|,\eta\rangle|\geq\sqrt{1/(d\log
n_{\mathrm{U}})}$. Then there exist $c_{4},c_{5}>0$, depending only on $r$ and
$c^{\prime}$, such that if $s\geq c_{4}\zeta\sqrt{d\log n_{\mathrm{U}}}$, then
on $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have

$\displaystyle|\alpha_{t}|\geq c_{5}s$

for some $t>0$.

###### Proof.

By flipping the sign of $\mu^{*}$ if necessary, we may assume without loss of
generality that $\alpha_{0}\geq 0$. Assuming that the desired result is not
true, we will prove by induction that on
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, (a) $\alpha_{t}/\beta_{t}\geq
c^{\prime}/(60\sqrt{d\log n_{\mathrm{U}}})$ and (b)
$\alpha_{t+1}\geq(1+\omega\sqrt{d\log n_{\mathrm{U}}})\alpha_{t}$ for all
$t\geq 0$. We show this by first verifying the base case of (a), then proving
that (a) implies (b) for each $t$, and finally proving that
$\alpha_{t+1}/\beta_{t+1}\geq 1/(60\sqrt{d\log n_{\mathrm{U}}})$ once (b)
holds for a given $t$.

For the base case, from the assumption on $\hat{\mu}^{(0)}$, we have

$\frac{\alpha_{0}}{\beta_{0}}\geq\frac{\alpha_{0}}{\|\hat{\mu}^{(0)}\|}=\bigg{\langle}\frac{\hat{\mu}^{(0)}}{\|\hat{\mu}^{(0)}\|},\eta\biggr{\rangle}\geq\frac{1}{\sqrt{d\log
n_{\mathrm{U}}}}\geq\frac{c^{\prime}}{60\sqrt{d\log n_{\mathrm{U}}}}$

since $c^{\prime}\leq 60$.

Now assume that $\alpha_{t}/\beta_{t}\geq c^{\prime}/(60\sqrt{d\log
n_{\mathrm{U}}})$ and that $\alpha_{0}\leq\alpha_{t}<c_{5}s$ for some $t\geq
0$. We aim to show that (b) holds for the same $t$, and start by controlling
$e_{1}^{\top}f(\hat{\mu}^{(t)})$. Let $W=(W_{1},\ldots,W_{d})^{\top}$ be an
independent copy of $Z_{1}$, independent of all other randomness in the
problem, and define
$u_{t}:=\bigl{\\{}\tanh(W_{1}\alpha_{t}+W_{-1}^{\top}\hat{\mu}^{(t)}_{-1})-\tanh(W_{1}\alpha_{t}-W_{-1}^{\top}\hat{\mu}^{(t)}_{-1})\bigr{\\}}/2$.
Then by applying the second part of Lemma 18 with $a=W_{1}\alpha_{t}$ and
$b=W_{-1}^{\top}\hat{\mu}^{(t)}_{-1}$ (so that $a+b=W^{\top}\hat{\mu}^{(t)}$),
we have

$\displaystyle e_{1}^{\top}f(\hat{\mu}^{(t)})$
$\displaystyle=\mathbb{E}\bigl{\\{}W_{1}\tanh(W_{1}\alpha_{t}+W_{-1}^{\top}\hat{\mu}^{(t)}_{-1})\bigr{\\}}$
$\displaystyle\geq\mathbb{E}\biggl{\\{}\alpha_{t}W_{1}^{2}-\frac{\alpha_{t}^{3}W_{1}^{4}}{3}-\alpha_{t}W_{1}^{2}(W_{-1}^{\top}\hat{\mu}^{(t)}_{-1})^{2}\biggr{\\}}+\mathbb{E}(W_{1}u_{t})$
$\displaystyle=\alpha_{t}(1+s^{2})(1-\beta_{t}^{2})-\alpha_{t}^{3}(1+2s^{2}+s^{4}/3),$

where in the final step we have used the fact that $u_{t}$ is an odd function
of $W_{-1}^{\top}\hat{\mu}^{(t)}_{-1}$, which has a symmetric distribution
about 0, conditional on $(W_{1},\hat{\mu}^{(t)}_{-1})$, and hence
$\mathbb{E}(W_{1}u_{t})=\mathbb{E}\bigl{\\{}\mathbb{E}(W_{1}u_{t}\mid
W_{1},\hat{\mu}_{-1}^{(t)})\bigr{\\}}=0$. From the assumption $s\geq
c_{4}\zeta$, and using $\beta_{t}\leq 60(1+r)\zeta$ from Proposition 10, we
have for sufficiently large $c_{4}$ that
$\beta_{t}^{2}(1+s^{2})\leq\\{60(1+r)\\}^{2}(1+r^{2})s^{2}/c_{4}^{2}\leq
s^{2}/4$. By choosing $c_{5}>0$, depending only on $r$, sufficiently small, we
may assume that
$\alpha_{t}^{2}(1+2s^{2}+s^{4}/3)<c_{5}^{2}(1+2r^{2}+r^{4}/3)s^{2}\leq
s^{2}/4$. Recall the definition of $f_{n_{\mathrm{U}}}$ from (28). Since on
the event $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have
$\|\hat{\mu}^{(t)}\|\leq r+3\leq 2(r+\sqrt{d})$ as in the first line of the
proof of Proposition 10, we have on the event
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$ that

$\displaystyle\alpha_{t+1}$
$\displaystyle=(1-\gamma)e_{1}^{\top}f_{n_{\mathrm{U}}}(\hat{\mu}^{(t)})+\gamma
e_{1}^{\top}\hat{\mu}_{n_{\mathrm{L}}}$
$\displaystyle\geq(1-\gamma)\bigl{\\{}e_{1}^{\top}f(\hat{\mu}^{(t)})-\omega\|\hat{\mu}^{(t)}\|\bigr{\\}}+\gamma(s-\phi)$
$\displaystyle\geq(1-\gamma)\bigl{\\{}\alpha_{t}(1+s^{2}/2)-\omega(\alpha_{t}+\beta_{t})\bigr{\\}}+\gamma(s-\phi).$

If $\gamma\geq\omega$, then $\phi\leq\omega\gamma^{-1/2}=\zeta\leq s/2$
(assuming $c_{4}\geq 2$). If $\gamma<\omega$, then

$\gamma\phi\leq\omega\gamma^{1/2}<\omega^{3/2}=\omega\zeta\leq\frac{\alpha_{0}\omega}{c^{\prime}}\sqrt{d\log
n_{\mathrm{U}}}\leq\frac{2(1-\gamma)\alpha_{t}\omega}{c^{\prime}}\sqrt{d\log
n_{\mathrm{U}}}.$

Hence, in either case, we have on $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$
that

$\displaystyle\alpha_{t+1}$
$\displaystyle\geq(1-\gamma)\alpha_{t}\bigl{\\{}1+s^{2}/2-\omega(1+62/c^{\prime})\sqrt{d\log
n_{\mathrm{U}}})\bigr{\\}}+\frac{\gamma s}{2}$
$\displaystyle\geq(1-\gamma)\alpha_{t}\bigl{\\{}1+s^{2}/2+(1+62/c^{\prime})(2\gamma-\omega\sqrt{d\log
n_{\mathrm{U}}})\bigr{\\}},$

where the final bound holds provided we reduce $c_{5}$ to be at most
$1/(2+124/c^{\prime})$ if necessary.

Now, when $\gamma\leq\omega\sqrt{d\log n_{\mathrm{U}}}$, we have
$\zeta=\omega^{1/2}\wedge\omega\gamma^{-1/2}\geq\omega^{1/2}(d\log
n_{\mathrm{U}})^{-1/4}$ and hence by the condition on $s$ in the proposition,
we have $s^{2}\geq c_{4}^{2}\zeta^{2}d\log n_{\mathrm{U}}\geq
c_{4}^{2}\omega\sqrt{d\log n_{\mathrm{U}}}$. Thus, by increasing $c_{4}$ to be
at least $\sqrt{4+248/c^{\prime}}$ if necessary, we have
$(1+62/c^{\prime})\omega\sqrt{d\log n_{\mathrm{U}}}\leq s^{2}/4$. Hence, in
this case, and on the event $\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$,

$\alpha_{t+1}\geq(1-\gamma)(1+s^{2}/4+(2+124/c^{\prime})\gamma)\alpha_{t}\geq\biggl{(}1+\frac{s^{2}}{8}\biggr{)}\alpha_{t}\geq(1+\omega\sqrt{d\log
n_{\mathrm{U}}})\alpha_{t}.$

On the other hand, when $\gamma>\omega\sqrt{d\log n_{\mathrm{U}}}$, we have

$\alpha_{t+1}\geq(1-\gamma)(1+(1+62/c^{\prime})\gamma)\alpha_{t}\geq(1+\gamma)\alpha_{t}\geq(1+\omega\sqrt{d\log
n_{\mathrm{U}}})\alpha_{t}.$

Combining the two bounds above proves (b) for this given $t$.

It remains to verify (a) for $t+1$, assuming that (a) and (b) hold up to and
including $t$. Since $\beta_{0}\leq\|\hat{\mu}^{(0)}\|\leq 60(\zeta\vee
r\omega)$, we have by Proposition 10 that $\beta_{t+1}\leq 60(\zeta\vee
r\omega)\leq 60(c^{\prime})^{-1}\|\hat{\mu}^{(0)}\|$. Thus,

$\frac{\alpha_{t+1}}{\beta_{t+1}}\geq\frac{\alpha_{0}}{60\|\hat{\mu}^{(0)}\|/c^{\prime}}\geq\frac{c^{\prime}}{60\sqrt{d\log{n_{\mathrm{U}}}}},$

which completes the induction. In particular, the geometric growth of
$\alpha_{t}$ implied by (b) means that $\alpha_{t}$ will exceed $c_{5}s$ for
sufficiently large $t>0$. This establishes our desired contradiction, and
hence proves the result. ∎

###### Proof of Proposition 5.

Define $\phi_{0}:=\omega_{0}\gamma^{-1/2}$, and recall the definitions of
$\Omega_{1}(\omega)$ and $\Omega_{2}(\phi)$ from (30). By Proposition 9, there
exists $C_{r}\geq 1$, depending only on $r$, such that for
$\omega=C_{r}\omega_{0}$ and $\phi=C_{r}\phi_{0}$, we have
$\mathbb{P}\bigl{(}\Omega_{1}(\omega)\cap\Omega_{2}(\phi)\bigr{)}\geq
1-2\delta$.

_(i)_ By the definition of $\omega$ and $\phi$, we have
$\omega=\phi\gamma^{1/2}$. If we choose $c$ such that $c\leq
C_{r}^{-1}\min\\{1/12,1/(r+3)\\}$, then $\omega\leq\min\\{1/12,1/(r+3)\\}$.
Thus, we may apply Propositions 10 and 11 to obtain that on
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)$, we have

$\limsup_{t\to\infty}\|\hat{\mu}^{(t)}-\mu^{*}\|\leq\limsup_{t\to\infty}(|\alpha_{t}|+\beta_{t}+\|\mu^{*}\|)\lesssim_{r}\zeta\vee\|\mu^{*}\|.$

The first claim follows.

_(ii)_ From Lemma 21 and by considering the case $d=1$ separately, for the
chosen $\eta_{0}$, we have

$\mathbb{P}\bigl{(}|e_{1}^{\top}\eta_{0}|\leq 1/\sqrt{d\log
n_{\mathrm{U}}}\bigr{)}\leq\sqrt{\frac{2}{\pi\log n_{\mathrm{U}}}}.$

Again, if we choose $c\leq C_{r}^{-1}\min\\{1/12,1/(r+3)\\}$, then
$\phi\gamma^{1/2}=\omega\leq\min\\{1/12,1/(r+3)\\}$. Also,
$\|\hat{\mu}^{(0)}\|=\zeta_{0}\vee r\omega_{0}\in[C_{r}^{-1}(\zeta\vee
r\omega),\zeta\vee r\omega]$. Thus, applying Proposition 13 with
$c^{\prime}=1/C_{r}$, there exists $c>0$, depending only on $r$, and
$t_{0}\in\mathbb{N}$ such that on
$\Omega_{1}(\omega)\cap\Omega_{2}(\phi)\cap\\{|e_{1}^{\top}\eta_{0}|>1/\sqrt{d\log
n_{\mathrm{U}}}\\}$, we have $|\alpha_{t_{0}}|\geq cs$.

Since $\beta_{0}\leq\|\hat{\mu}^{(0)}\|\leq\zeta\vee r\omega$, we can apply
Proposition 10 to obtain that $\beta_{t}\leq 60(\zeta\vee r\omega)$ for all
$t\geq 0$. Hence all conditions of Proposition 12 are satisfied, and the
desired result then follows from (48). ∎

To prove Theorem 6, we need the following proposition, which relates the loss
of estimating $\mu^{*}$ to the operator norm loss of estimating
$\mu^{*}\mu^{*\top}$.

######  Proposition 14.

Assume that $n\geq 3$ and that
$(Z_{1},Y_{1},Y_{1}^{*}),\ldots,(Z_{n},Y_{n},Y_{n}^{*})$ are independent with

$Y_{i}^{*}\sim\mathrm{Unif}\\{-1,1\\},\;Z_{i}\mid
Y_{i}^{*}\sim\mathcal{N}_{d}(Y_{i}^{*}\mu^{*},I_{d}),\;Y_{i}=Y_{i}^{*}\mathbbm{1}_{\\{i\leq
n_{\mathrm{L}}\\}}\quad\text{for $i\in[n]$.}$

For $\mu\in\mathbb{R}^{d}$ and $i\in[n]$, let
$L_{i}(\mu):=Y_{i}\mathbbm{1}_{\\{Y_{i}\neq 0\\}}+\tanh\langle
Z_{i},\mu\rangle\mathbbm{1}_{\\{Y_{i}=0\\}}$, $\mu_{\mathrm{tot}}(\mu):=\mu
n^{-1}\sum_{i=1}^{n}L_{i}(\mu)$ and
$\Sigma_{\mathrm{b}}(\mu):=\mu\mu^{\top}-\mu_{\mathrm{tot}}(\mu)\mu_{\mathrm{tot}}(\mu)^{\top}$.
For any $\delta\in(0,1)$ and $B>0$, we have with probability at least
$1-\delta$ that

$\sup_{\mu:\|\mu\|\leq
B}\bigl{\\{}\|\Sigma_{\mathrm{b}}(\mu)-\mu^{*}\mu^{*\top}\|_{\mathrm{op}}-(B+s)L(\mu,\mu^{*})\bigr{\\}}\lesssim
B^{2}(s^{2}\vee 1)\biggl{(}\frac{d\log(2Bn+e)+\log(1/\delta)}{n}\biggr{)}.$

###### Proof.

For any $\mu\in\mathbb{R}^{d}$ with $\|\mu\|\leq B$, we have

$\displaystyle\|\Sigma_{\mathrm{b}}(\mu)-\mu^{*}\mu^{*\top}\|_{\mathrm{op}}$
$\displaystyle\leq\|\mu\mu^{\top}-\mu^{*}\mu^{*\top}\|_{\mathrm{op}}+\|\mu_{\mathrm{tot}}(\mu)\mu_{\mathrm{tot}}(\mu)^{\top}\|_{\mathrm{op}}$
$\displaystyle\leq(B+s)L(\mu,\mu^{*})+B^{2}\biggl{(}\frac{1}{n}\sum_{i=1}^{n}L_{i}(\mu)\biggr{)}^{2}.$

Thus, it is enough to show that $\sup_{\mu:\|\mu\|\leq B}$
$n^{-1}\sum_{i=1}^{n}$ $L_{i}(\mu)$ $\lesssim\sqrt{\frac{(s^{2}\vee
1)\\{d\log(2Bn+e)+\log(1/\delta)\\}}{n}}$ with probability at least
$1-\delta$. To this end, we have

$\displaystyle\sup_{\mu:\|\mu\|\leq B}\frac{1}{n}\sum_{i=1}^{n}L_{i}(\mu)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n_{\mathrm{L}}}Y_{i}+\frac{1}{n}\sup_{\mu:\|\mu\|\leq
B}\sum_{i=n_{\mathrm{L}}+1}^{n}\tanh\langle Z_{i},\mu\rangle.$ (60)

For the first term on the right-hand side of (60), by Hoeffding’s inequality,
we have

$\mathbb{P}\biggl{(}\frac{1}{n}\sum_{i=1}^{n_{\mathrm{L}}}Y_{i}>\frac{\sqrt{2n_{\mathrm{L}}\log(3/\delta)}}{n}\biggr{)}\leq\frac{\delta}{3}.$
(61)

For the second term on the right-hand side of (60), let $\mathcal{N}$ be a
$\varepsilon$-net of $\\{v:\|v\|\leq B\\}$ with respect to the Euclidean
distance, for some $\varepsilon\in(0,1/2]$ to be specified later. Since a
maximal $\varepsilon$-packing set is an $\varepsilon$-net, we may assume that
$|\mathcal{N}|\leq(B+\varepsilon/2)^{d}/(\varepsilon/2)^{d}=(2B/\varepsilon+1)^{d}$.
Using the fact that $x\mapsto\tanh x$ is $1$-Lipschitz, together with the
Cauchy–Schwarz inequality, we have

$\displaystyle\sup_{v:\|v\|\leq B}\sum_{i=n_{\mathrm{L}}+1}^{n}$
$\displaystyle\tanh\langle Z_{i},v\rangle$
$\displaystyle\leq\sup_{v\in\mathcal{N}}\sum_{i=n_{\mathrm{L}}+1}^{n}\tanh\langle
Z_{i},v\rangle+\sup_{u,v:\|u-v\|\leq\varepsilon}\sum_{i=n_{\mathrm{L}}+1}^{n}(\tanh\langle
Z_{i},u\rangle-\tanh\langle Z_{i},v\rangle)$
$\displaystyle\leq\sup_{v\in\mathcal{N}}\sum_{i=n_{\mathrm{L}}+1}^{n}\tanh\langle
Z_{i},v\rangle+\varepsilon\sum_{i=n_{\mathrm{L}}+1}^{n}\|Z_{i}\|.$

Hence taking $\varepsilon=1/n$, and defining
$\tau:=\frac{\log(3/\delta)}{d\log(2Bn+e)}>0$, we have

$\displaystyle\mathbb{P}\biggl{(}\frac{1}{n}\sup_{v:\|v\|\leq
B}\sum_{i=n_{\mathrm{L}}+1}^{n}\tanh\langle
Z_{i},v\rangle>2\sqrt{\frac{2(s^{2}\vee 1)(1+\tau)d\log(2Bn+e)}{n}}\biggr{)}$
$\displaystyle\leq\mathbb{P}\biggl{(}\frac{1}{n}\sup_{v\in\mathcal{N}}\sum_{i=n_{\mathrm{L}}+1}^{n}\tanh\langle
Z_{i},v\rangle>\sqrt{\frac{2(s^{2}\vee 1)(1+\tau)d\log(2Bn+e)}{n}}\biggr{)}$
$\displaystyle\hskip
85.35826pt+\mathbb{P}\biggl{(}\frac{1}{n}\sum_{i=n_{\mathrm{L}}+1}^{n}\|Z_{i}\|\geq\sqrt{2(s^{2}\vee
1)(1+\tau)nd\log(2Bn+e)}\biggr{)}$