# Sparse Equisigned PCA:
Algorithms and Performance Bounds in the Noisy Rank-1 Setting111This work was
partially supported by the NSF and the ONR.

Arvind Prasadan222University of Michigan, Dept of EECS<EMAIL_ADDRESS>Raj Rao Nadakuditi 333University of Michigan, Dept of EECS<EMAIL_ADDRESS>and Debashis Paul444University of California, Davis, Department of Statistics,
<EMAIL_ADDRESS>
###### Abstract

Singular value decomposition (SVD) based principal component analysis (PCA)
breaks down in the high-dimensional and limited sample size regime below a
certain critical eigen-SNR that depends on the dimensionality of the system
and the number of samples. Below this critical eigen-SNR, the estimates
returned by the SVD are asymptotically uncorrelated with the latent principal
components. We consider a setting where the left singular vector of the
underlying rank one signal matrix is assumed to be sparse and the right
singular vector is assumed to be equisigned, that is, having either only
nonnegative or only nonpositive entries. We consider six different algorithms
for estimating the sparse principal component based on different statistical
criteria and prove that by exploiting sparsity, we recover consistent
estimates in the low eigen-SNR regime where the SVD fails. Our analysis
reveals conditions under which a coordinate selection scheme based on a sum-
type decision statistic outperforms schemes that utilize the $\ell_{1}$ and
$\ell_{2}$ norm-based statistics. We derive lower bounds on the size of
detectable coordinates of the principal left singular vector and utilize these
lower bounds to derive lower bounds on the worst-case risk. Finally, we verify
our findings with numerical simulations and a illustrate the performance with
a video data where the interest is in identifying objects.

##  1 Introduction

It is well-understood that singular value decomposition (SVD) based principal
component analysis (PCA) breaks down in the high-dimensional and limited
sample size regime below a certain critical eigen-SNR (eigenvalue signal-to-
noise ratio) that depends on the dimensionality of the system and the number
of samples [21, 5]. Several sparse PCA algorithms have been proposed in the
literature (see [21, 5, 11, 25, 41, 4]) and have been shown to successfully
estimate the principal components in the low eigen-SNR regime where the SVD
fails.

Prior work in this area primarily considers the Gaussian signal-plus-noise
model with random effects, where the signal matrix is assumed to have sparse
left singular vectors, normally distributed right singular vectors, and the
noise matrix is assumed to have normally distributed _i.i.d._ entries. Here,
we consider the setting where the left singular vector of the rank one signal
matrix is sparse _and_ the right singular vector is assumed to be equisigned.
We say that a vector is equisigned if its entries are all non-negative or all
non-positive. This is motivated by applications such as diffusion imaging in
MRI where the right singular vector represents a physical quantity (e.g.
intensity as the diffusion agent is absorbed by a tissue) that is non-
negative, by imaging problems such as foreground-background separation in
video data [32, 39] and object detection in astronomy [35], where the data are
naturally non-negative, and by problems in bioinformatics where the data are
(non-negative) counts of genes [38]. When analyzing data that are non-
negative, it is logical to take advantage of this property, and investigate
how we may use this knowledge to do better than the (generic) alternatives.
Alternatively, a practitioner may seek to use techniques that constrain or
impose non-negativity to preserve interpretability of the results, e.g., non-
negative matrix factorization. Additionally, we motivate the rank-$1$
assumption by noting that for a video with a static background, the foreground
is a perturbation of a rank-$1$ background [27, 16]. Finally, even though we
do not pursue this angle here, our framework can be extended to deal with the
scenario where the signal can be viewed of a rank 1 tensor with all but one of
the representors in the Kroneker product representation of the tensor is an
equisigned vector.

There is precedent for and prior work on non-negative PCA, including the
sparse biased PCA in [9], the sparse PCA with non-negativity priors in [31],
and the work in [26]. These works differ from our work in that they impose
non-negativity on the factors or the left singular vectors. In this work, we
study sparse factors with non-negative loadings; i.e., we are solving a
different problem in this work.

A natural question at this juncture is the following: _how does our problem
differ from that solved by Non-Negative Matrix Factorization (NNMF)?_ NNMF
takes a given matrix $\operatorname{\mathbf{X}}$ and looks for non-negative
matrices $\operatorname{\mathbf{F}}$ and $\operatorname{\mathbf{G}}$ such that
$\operatorname{\mathbf{X}}=\operatorname{\mathbf{F}}\operatorname{\mathbf{G}}^{T}$
[19, 40]. Ordinary NNMF has no sparsity constraints. We might impose such
constraints, as is done in [18] and [24], but except in special cases, these
solutions have no known theoretical guarantee of statistical performance. This
problem partly stems from the fact that solutions to the corresponding
optimization problems may not be unique. In contrast, our problem only
constrains the right singular vectors, while the left singular vectors are
free to take any sign. The work in [12] extends the NNMF framework to one
wherein only one of the factors is non-negative; nevertheless, the rest of the
constraints we impose are not included. The work in [42] seeks factors (left
singular vectors) with disjoint supports and non-negative loadings, but this
definition of sparsity does not match that from the sparse PCA literature.
Hence, NNMF is not an answer to the problem we consider herein.

The main contribution of this paper is a rigorous sparsistency analysis of the
various algorithms that brings into focus the various very-low eigen-SNR
regimes where the new algorithms work and the SVD based methods provably fail.
Additionally, a major novelty of this work is the integration of FDR-
controlling (False Discovery Rate) hypothesis testing to the Sparse PCA
problem.

Our analysis illustrates the situations where the sum based coordinate
selection scheme dramatically outperforms the $\ell_{1}$ and $\ell_{2}$ [21,
5] based sparse PCA schemes. Additionally, our proposed algorithms are non-
iterative, do not require the computation of the sample covariance matrix, and
do not require knowledge of the sparsity level. We separate our algorithms
into two groups: one where the Family-Wise Error Rate (FWER) is controlled,
and another where the False Discovery Rate (FDR) is controlled. We utilize
sharp tail probability bounds for relevent statistics to derive our FWER-
controlling estimators [7]. For the FDR controlling estimators, we relate the
problem at hand to that of the sparse normal means problem [13].

This paper is organized as follows. In Section 3, we describe three algorithms
for estimating the sparse principal component that utilize a coordinate
selection scheme based on the sum, $\ell_{1}$, and $\ell_{2}$ norm-based
statistics respectively. We call our family of algorithms _SEPCA_ , an
abbreviation for Sparse Equisigned PCA. Section 4 proposes three FDR-
controlling refinements of the sum- and $\ell_{2}$-based algorithms in Section
3 by relating coordinate detection to the sparse normal means estimation
problem. In Section 5 we show how the estimation performance is governed by
the size of the smallest detectable coordinate, which we analyze in Section 6
and validate using numerical simulations in Section 7. In Section 8, we
provide some geometric intuitions about the relative performance of three of
our algorithms. We show that the sum statistic is potentially the most
powerful, while the $\ell_{1}$ is the least powerful. We provide some
concluding remarks in Section 9.

##  2 Problem Formulation

Let $\operatorname{\mathbf{X}}\in\mathbb{R}^{p\times n}$ be a real-valued
signal-plus-noise data matrix of the form

$\operatorname{\mathbf{X}}=\theta\operatorname{\mathbf{u}}\operatorname{\mathbf{v}}^{T}+\sigma\operatorname{\mathbf{G}}.$
(1)

The columns of the $p\times n$ data matrix $\operatorname{\mathbf{X}}$
represent $p$-dimensional observations. In (1), $\operatorname{\mathbf{u}}$
and $\operatorname{\mathbf{v}}$ are the left and right singular vectors of the
rank-one latent signal matrix, and have entries $u_{i}$ and $v_{j}$,
respectively. The entries of $\operatorname{\mathbf{G}}$, the noise matrix,
are assumed to be _i.i.d._ Gaussian random variables with mean $0$ and
variance ${1}/{n}$. We assume that
$\operatorname{\mathbf{u}}\in\mathbb{R}^{p}$ has unit norm and is sparse in
the sense of small $\ell_{0}$ norm, with $s\ll p$ non-zero entries, where
$s/n\rightarrow 0$. That is, for a set
$I=\\{i_{1},\cdots,i_{s}\\}\subset\\{1,\cdots,p\\}$,

$\begin{array}[]{ll}u_{i}\neq 0&\text{for $i\in I$,}\\\ u_{i}=0&\text{for
$i\in I^{C}$,}\end{array}$ (2)

where $I^{C}$ denotes the complement of $I$. We further assume
$\operatorname{\mathbf{v}}\in\mathbb{R}^{n}$ to be of unit norm,
deterministic, and equisigned. Given $\operatorname{\mathbf{X}}$, our goal is
to recover $\operatorname{\mathbf{u}}$ and $\operatorname{\mathbf{v}}$.

Note that the $(i,k)$ entry of $\operatorname{\mathbf{X}}$, $X_{ik}$, is a
Gaussian random variable with mean $[\theta u_{i}]v_{k}$ and variance
${\sigma^{2}}/{n}$. Moreover, it follows that

$\mathbb{E}(\operatorname{\mathbf{X}}\operatorname{\mathbf{X}}^{T})=\theta^{2}\operatorname{\mathbf{u}}\operatorname{\mathbf{u}}^{T}+{\sigma^{2}}\operatorname{\mathcal{I}}_{p},$

where $\operatorname{\mathcal{I}}_{p}$ denotes the $p\times p$ identity
matrix. The quantity $\left(\theta/\sigma\right)^{2}$ is, for this model, the
eigen-SNR (signal-to-noise ratio).

###  2.1 Motivation: Breakdown of PCA / SVD

From [3], we have the following result: let
$\widehat{\operatorname{\mathbf{u}}}$ be the estimate of
$\operatorname{\mathbf{u}}$ given by the Singular Value Decomposition (SVD) of
$\operatorname{\mathbf{X}}$, and let ${p(n)}/{n}$ have limit $c\in[0,\infty]$
as $n$ grows, with $\theta$ fixed and $\sigma=1$. Then, with probability $1$,

$\left|\langle\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}}\rangle\right|^{2}\rightarrow\left\\{\begin{array}[]{ll}1-\frac{c\left(1+\theta^{2}\right)}{\theta^{2}\left(c+\theta^{2}\right)}&\text{if
$\theta\geq c^{1/4}$,}\\\ 0&\text{otherwise}.\end{array}\right.$ (3)

For general $\sigma$, we replace $\theta$ by $\theta/\sigma$ in (3). Hence,
SVD based PCA leads to inconsistent estimates of $\operatorname{\mathbf{u}}$
(and also for $\operatorname{\mathbf{v}}$, which can be deduced from (3)) when
the dimension $p$ is comparable to or larger than the sample size $n$.
Moreover, in the low eigen-SNR regime, the estimates break down completely.
SVD does not exploit any assumed structure in $\operatorname{\mathbf{u}}$ and
$\operatorname{\mathbf{v}}$. Consequently, (3) holds for arbitrary
$\operatorname{\mathbf{u}}$ and $\operatorname{\mathbf{v}}$, including our
setting where $\operatorname{\mathbf{u}}$ is sparse and/or
$\operatorname{\mathbf{v}}$ is equisigned. Our goal, in what follows, is to
derive consistent estimators for $\operatorname{\mathbf{u}}$ and
$\operatorname{\mathbf{v}}$ that outperform the SVD by exploiting the sparsity
of $\operatorname{\mathbf{u}}$ and the equisigned nature of
$\operatorname{\mathbf{v}}$.

###  2.2 Problem Statement

Note that we have assumed that $\operatorname{\mathbf{u}}$ is sparse and that
the sparsity $s$ is such that $s/n$ has limit zero. Hence, if we had oracle
knowledge of the sparsity pattern $I$ (the indices of
$\operatorname{\mathbf{u}}$ that have non-zero coordinates), restricting the
matrix $X$ to those rows indexed by $I$ and performing the SVD on the smaller
matrix would yield a consistent estimator for the non-zero elements of
$\operatorname{\mathbf{u}}$ and the vector $\operatorname{\mathbf{v}}$. This
conclusion follows from (3), since the value $c$ is replaced with $s/n$, which
has limit zero. Thus, if we derived consistent estimators of the support of
$\operatorname{\mathbf{u}}$, we have a consistent two-stage estimation
procedure of the vectors $\operatorname{\mathbf{u}}$ and
$\operatorname{\mathbf{v}}$.

Formally, we are interested in finding a procedure that estimates $I$ by
$\widehat{I}$ such that in the limit $n\rightarrow\infty$,

$\left\\{\begin{array}[]{ll}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)\rightarrow
1&\textrm{ if $i\in I$},\\\
\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)\rightarrow 0&\textrm{ if
$i\notin I$}.\end{array}\right.$ (4)

Equivalently, noting that the Hamming distance of $I$ and $\widehat{I}$,
denoted by $d_{H}\left(I,\widehat{I}\right)$, is given by the cardinality of
their symmetric set difference,

$d_{H}\left(I,\widehat{I}\right)=\left|\left(I\cup\widehat{I}\right)\backslash\left(I\cap\widehat{I}\right)\right|,$

we want the expected Hamming distance
$\operatorname{\mathbb{E}}d_{H}\left(I,\widehat{I}\right)$ to have limit $0$,
which is stronger than requiring consistency in recovering the support (or
sparsity pattern) of $\operatorname{\mathbf{u}}$. However, as the work in [10,
34] indicates, this limit will not in general be zero, and will depend on the
noise level, signal strength, and sparsity.

##  3 Proposed Algorithms

We propose six different two-stage algorithms for estimating
$\operatorname{\mathbf{u}}$. The first three algorithms are designed to
control the family-wise error rate (FWER), or, the probability of obtaining a
false positive in the coordinate selection. The last three algorithms aim to
control the false discovery rate (FDR), or, the proportion of false
discoveries (coordinate detections) among all discoveries. We defer discussion
of the FDR-based algorithms to Section 4.

All of the algorithms have the same basic form given in Algorithm 1. Given
$\operatorname{\mathbf{X}}$, we associate a test statistic $T_{i}$ to each row
of $\operatorname{\mathbf{X}}$. The sparsity of $\operatorname{\mathbf{u}}$
implies that the majority of the rows of $\operatorname{\mathbf{X}}$ are
purely noise, so that the majority of the $T_{i}$ come from the null, noise-
only distribution. Hence, based on the statistics $\\{T_{i}\\}$, we perform a
form of multiple hypotheses testing procedure, and select the set
$\widehat{I}$ of indices that are non-null. In this way, we can estimate the
_support_ of $\operatorname{\mathbf{u}}$, thereby isolating the the rows of
$\operatorname{\mathbf{X}}$ that contain the signal. Then, taking the SVD of
this submatrix (comprised of only the selected rows of
$\operatorname{\mathbf{X}}$) yields a better estimate of the non-zero
coordinates in $\operatorname{\mathbf{u}}$, as well as
$\operatorname{\mathbf{v}}$.

Algorithm 1 Variable Selection and Estimation Algorithm

0: Threshold $\tau_{n,p}$ and form of Test Statistic $T_{i}$ from Table 1

Let $\widehat{I}$ be an empty list

for all Rows $i$ of $\operatorname{\mathbf{X}}$, $1\leq i\leq p$ do

Form test statistic $T_{i}$ from row $i$ of $\operatorname{\mathbf{X}}$

if $T_{i}\geq\tau_{n,p}$ then

Add $i$ to $\widehat{I}$

end if

end for

Let
$[\widetilde{\operatorname{\mathbf{u}}},\widetilde{\theta},\widehat{\operatorname{\mathbf{v}}}]=\text{SVD}(\operatorname{\mathbf{X}}_{\widehat{I},:})$
be the rank-$1$ SVD of $\operatorname{\mathbf{X}}$ restricted to rows in
$\widehat{I}=[i_{1},\cdots,i_{|\widehat{I}|}]$

For $i_{k}\in\widehat{I}$, let $\widehat{u}_{i_{k}}=\widetilde{u}_{k}$; the
other entries of $\widehat{\operatorname{\mathbf{u}}}$ are set to $0$.

We begin by discussing the FWER-controlling algorithms. The work in [21]
proposed a covariance thresholding method for Sparse PCA called DT-SPCA; this
is equivalent to a coordinate selection scheme based on the $\ell_{2}$ norm-
based statistic. In our terminology and with our choice of thresholds, we
label it as $\ell_{2}$-SEPCA. We label the coordinate selection scheme based
on the $\ell_{1}$ norm-based statistic $\ell_{1}$-SEPCA. Finally, the sum-
SEPCA algorithm utilizes row sums of the data matrix.

###  3.1 Computational Complexity

Note that the variable selection part of our procedures has a computational
complexity that is $O\left(pn\right)$: the formation of the test statistic is
linear in the number of columns, and the formation is repeated once per row.
Noting that for a $p\times n$ matrix, the complexity of the rank-$1$ SVD is
$O\left(1\times pn\right)$, we find that if $\left|\widehat{I}\right|$
coordinates are selected, we have an overall complexity of
$O\left(pn+\left|\widehat{I}\right|n\right)=O\left(pn\right)$ [1].

Computation of the covariance matrix has a (naive) complexity of
$O\left(p^{2}n\right)$, and in practice is somewhere between
$O\left(p^{2}\right)$ and $O\left(p^{3}\right)$ [15]. Immediately, our methods
here are faster than those requiring explicit formation of the covariance
matrix [5, 21, 25, 41, 4]. Additionally, there is no iteration or convergence
of any optimization problems required. Note that a semi-definite programming-
based formulation is at least polynomial in the problem size:
$O\left(p^{4}\right)$ [11] or $O\left(p^{3}\right)$ [4]. The ITSPCA method
applied to our rank-$1$ setting would have a cost of $O\left(ps\right)$ per
iteration [25, Sec. 4]. TPower has a similar complexity of
$O\left(sp+p\right)$ per iteration [41].

###  3.2 The DT-SPCA Algorithm and Two-Stage Procedures

The DT-SPCA algorithm was proposed in [21] and later used as the first stage
of the ASPCA algorithm given in [5]. The algorithm thresholds the diagonals of
the matrix $XX^{T}$ to perform variable selection: note that in our setting,
these values are

${\theta^{2}}u_{i}^{2}+{\sigma^{2}}\sum_{j=1}^{n}G_{ij}^{2},$

with expectation $\left(\theta^{2}u_{i}^{2}+\sigma^{2}\right)$. The DT-SPCA
algorithm thresholds these diagonal values at
$\sigma^{2}\left(1+\gamma\sqrt{\frac{\log p}{n}}\right)$, where $\gamma>0$,
and then performs PCA on the reduced matrix formed from the selected
variables. Noting that the diagonals of $XX^{T}$ are the same as the row sum-
of-squares of $X$, we see that $\ell_{2}$-SEPCA is essentially the same (up to
choice of threshold) as DT-SPCA.

However, the innovation of [21, 5] that we carry forward is the two-stage
procedure. That is, we perform some sort of testing to estimate the support of
the sparse singular vector $\operatorname{\mathbf{u}}$, and then perform an
SVD on the reduced matrix. As we will see in what follows, there is
flexibility in the choice of testing or support estimation method.

###  3.3 Statement of Thresholds

We shall choose the thresholds $\tau_{n,p}$ for the coordinate selection
scheme so that in the noise-only case,

$\mathbb{P}\left(\max_{1\leq i\leq
p}T_{i}\geq\tau_{n,p}\right)\leq\frac{1}{ep}\rightarrow 0,$ (5)

where $e$ is Euler’s number, or the base of the natural logarithm. This choice
ensures that the probability of a false positive tends to zero as
$p\to\infty$. That is, the FWER is asymptotically zero and is bounded by
$1/ep$ in the finite-dimensional case. Note that the constraint used to
control the FWER is simply that the distribution of the noise is log-concave.
In the Gaussian case, we obtain the specific expressions given summarized in
Table 1; however, with knowledge of the moments $\mathbb{E}T_{i}$ and
$\text{Var }T_{i}$, we can repeat our analysis and find thresholds for the
$\ell_{1}$ and $\ell_{2}$-SEPCA algorithms with _any_ log-concave noise
distribution. The thresholds are summarized in Table 1.

Table 1: Test Statistics and Thresholds for Algorithm (1)

Algorithm | Statistic $T_{i}$ | Threshold $\tau_{n,p}$
---|---|---
$\ell_{1}$-SEPCA | $\frac{1}{\sqrt{n}}\sum_{k=1}^{n}|X_{i,k}|$ | $\sigma\left(\sqrt{\frac{2}{\pi}}+C_{1}\frac{\log ep}{\sqrt{n}}\right)$
$\ell_{2}$-SEPCA | $\sum_{k=1}^{n}X_{i,k}^{2}$ | $\sigma^{2}\left(1+C_{2}\frac{\log ep}{\sqrt{n}}\right)$
sum-SEPCA | $\frac{1}{\sqrt{n}}|\sum_{k=1}^{n}X_{i,k}|$ | $\sigma C_{U}\sqrt{\frac{\log p}{n}}$

*See (6) and (10) for definitions of the constants $C_{2}$, $C_{1}$, and $C_{U}$.

  * Remark

Note that we impose strong control over the FWER and seek to reject individual
null hypotheses, instead of weak control and considering the global null
hypothesis as in [4].

###  3.4 FWER Thresholds

####  3.4.1 $\ell_{2}$\- and $\ell_{1}$-SEPCA

In the noise-only cases, the statistics for $\ell_{2}$\- and $\ell_{1}$-SEPCA
are distributed as scaled $\chi_{n}^{2}$ and sums of half-normal,
respectively. Both of these quantities are log-concave random variables, so we
may apply the result in [23] to set the threshold $\tau_{n,p}$ in both cases.

Defining $K$ to be some absolute constant (we may use $K=e$, as in [6]), we
define the constants

$C_{2}=\sqrt{2}K\text{ and }C_{1}=K\sqrt{\left(1-{2}/{\pi}\right)}.$ (6)

####  3.4.2 sum-SEPCA

From Proposition $4.4$ of [8], we obtain that the threshold for sum-SEPCA is
given by

$\displaystyle\tau_{n,p}=\frac{\sigma}{\sqrt{n}}\left(\sqrt{2\log
p}+\frac{1}{U(p)}\left(\frac{1}{3}\log{ep}+\sqrt{\log{ep}}\right)+\delta_{p}\right).$
(7)

In (7), we have that

$U(p)=\sqrt{2}\text{ Erf}^{-1}\left(1-\frac{1}{p}\right)\text{ and
}\delta_{p}\asymp\frac{\pi^{2}}{12}\left(\log p\right)^{-3/2},$ (8)

where Erf denotes the _error function_ , or alternatively, the cumulative
distribution function of a standard Gaussian random variable is given by

$\Phi(x)=\frac{1}{2}\left(1+\text{Erf}\left(\frac{x}{\sqrt{2}}\right)\right).$
(9)

Moreover, $\tau_{n,p}\leq\sigma C_{U}\sqrt{\frac{\log p}{n}}$ for some
constant $C_{U}$. For a fixed value of $p$, choosing

$\kappa_{U}\geq\frac{\sqrt{2}}{U(p)}\left(3+\sqrt{\log{p}}\right)>1\text{ and
}C_{U}=\sqrt{2}+\frac{\kappa_{U}}{3\sqrt{2}}$ (10)

is sufficient. The choice of $1/ep$ is the largest bound justified by
Proposition 4.4 of [8], so we have calibrated all of our algorithms to the
same constant factor times $1/p$. The thresholds are summarized in Table 1.

###  3.5 Estimation of the Noise Variance, $\sigma^{2}$

In this work, we assume that the noise variance $\sigma^{2}$ is known;
however, in general, estimation of $\sigma^{2}$ may not be straightforward
[29]. Recently proposed procedures such as those proposed in [29, 30, 37]
could be employed to estimate the noise variance, and we point the interested
reader to these references for more theoretical background on the problem. We
note that in most applications, including the video example we consider, one
can obtain a relatively sparse representation of the object in a multiscale
basis such as a wavelet basis [20, Sec. 7.5]. Under such circumstances, under
the assumed additive, isotropic noise model, we can easily obtain a consistent
estimate of $\sigma^{2}$ by utilizing the inherent sparsity of the signal,
especially in finer scales. This can be done, for example, by computing the
variance of the wavelet coefficients in the finest scale [20, Sec. 7.5]. One
can obtain a more robust estimate by taking the median absolute deviation of
the coefficients about their median and then by multiplying its square with a
known scale factor (assuming normality) [30, 20].

##  4 Controlling the False Discovery Rate

So far, we have controlled the probability of a false alarms when detecting
coordinates. However, there are two relevant observations to make. First,
under the Gaussian noise, rank-1, and equisigned assumptions, the vector of
test statistics $\\{T_{i}\\}$ in the sum-SEPCA algorithm looks like a sparse
vector plus Gaussian noise (or a vector of $\chi_{n}^{2}$-variates with
varying non-centralities, in the $\ell_{2}$-SEPCA algorithm). Secondly,
controlling the false discovery rate, that is, the proportion of rejected
nulls that are false positives, can lead to increased detection power relative
to controlling the false positive rate. We hence look at FDR-controlling tests
for the _Sparse Normal Means_ problem.

That is, given a vector of test statistics (as before), we replace the
thresholding and selection in Algorithm 1 with an FDR-controlling selection
procedure. We summarize this change in Algorithm 2. There are three procedures
we consider. The first two are known as Higher Criticism, and directly extend
the sum- and $\ell_{2}$-SEPCA algorithms [13, 14]. The third is a method for
detection in the sparse normal means problem that comes out of complexity-
penalized estimation theory for linear inverse problems [22].

Algorithm 2 FDR-Controlling Variable Selection and Estimation Algorithm

0: Test Statistic $T_{i}$ from Table 1 and Selection Procedure

Let $\widehat{I}$ be an empty list

for all Rows $i$ of $\operatorname{\mathbf{X}}$, $1\leq i\leq p$ do

Form test statistic $T_{i}$ from row $i$ of $\operatorname{\mathbf{X}}$

end for

Perform an FDR-Controlling selection procedure, and add the selected indices
to $\widehat{I}$

Let
$[\widetilde{\operatorname{\mathbf{u}}},\widetilde{\theta},\widehat{\operatorname{\mathbf{v}}}]=\text{SVD}(\operatorname{\mathbf{X}}_{\widehat{I},:})$
be the rank-$1$ SVD of $\operatorname{\mathbf{X}}$ restricted to rows in
$\widehat{I}=[i_{1},\cdots,i_{|\widehat{I}|}]$

For $i_{k}\in\widehat{I}$, let $\widehat{u}_{i_{k}}=\widetilde{u}_{k}$; the
other entries of $\widehat{\operatorname{\mathbf{u}}}$ are set to $0$.

###  4.1 Higher Criticism

####  4.1.1 Formulation

Assume we have $p$ independent tests of the form

$\begin{array}[]{ll}H_{o,i}:&W_{i}\sim\mathcal{N}\left(0,1\right),\\\
H_{1,i}:&W_{i}\sim\mathcal{N}\left(\mu_{i},1\right),\end{array}$ (11)

and assume that at most $p^{1-\beta}$ of the $p$ hypotheses are truly non-
null, for some $\beta\in(1/2,1)$. Further assume that the non-null means have
magnitude

$\mu_{i}=\mu_{p}=\sqrt{2r\log p},$

for $r\in(0,1)$. Here, the means will correspond to the coordinate size. Note
that the expected maximum of $p$ standard Gaussian random variables is upper
bounded by $\sqrt{2\log p}$, with the bound being asymptotically sharp.

If we let $p_{(1)}\leq p_{(2)}\leq\cdots\leq p_{(p)}$ be the sorted p-values
of the individual tests, we may define the Higher Criticism statistic:

$HC_{p}=\max_{i:1/p\leq p_{(i)}\leq
1/2}\frac{\sqrt{p}\left(i/p-p_{(i)}\right)}{\sqrt{p_{(i)}(1-p_{(i)})}}.$ (12)

Rejecting the global null hypothesis (that there are no non-null coordinates)
when $HC_{p}>\sqrt{2\log\log p}(1+o(1))$ leads to asymptotically full power
when $r$ is greater than some decision boundary $\rho$, and that under the
global null,

$\frac{HC_{p}}{\sqrt{2\log\log p}}\rightarrow 1$ (13)

in probability as $n,p\rightarrow\infty$. The function $\rho$ depends on the
sparsity index $\beta$, and as [14] indicate:

$\rho(\beta)=\left\\{\begin{array}[]{ll}\beta-1/2&\text{ when
}\beta\in(1/2,3/4),\\\ \left(1-\sqrt{1-\beta}\right)^{2}&\text{ when
}\beta\in(3/4,1).\end{array}\right.$ (14)

If we replace the normal distribution with a $\chi_{n}^{2}$ distribution, the
same results hold for tests of the form

$\begin{array}[]{ll}H_{o,i}:&W_{i}\sim{\chi_{n}^{2}},\\\
H_{1,i}:&W_{i}\sim{\chi_{n}^{2}}(\delta),\end{array}$ (15)

where $\delta$ is a non-centrality parameter and we consider $r\in(0,1)$ such
that $\delta=2r\log p$.

  * Remark

While Higher Criticism is typically formulated for the case of identical non-
null means or parameters (all of the non-zero $\mu_{i}$ are identical), this
constraint is not mandatory [2, 17]. Indeed, the results hold without
modification for the Gaussian model with non-null means of size
$\mu_{i}=\alpha_{i}\sqrt{2\log p}$, where $\alpha_{i}$ is a non-negative
random variable with the property that $\mathbb{P}(\alpha_{i}\leq\sqrt{r})=1$
and $\mathbb{P}(\alpha_{i}>\sqrt{r}-\epsilon)>0$ for all $\epsilon>0$ [17].
The case of a $\chi_{n}^{2}$ distribution is similar.

As a point of interest, the test in (11) can be extended to (and potentially
strengthened in) the case where the $p$ tests are correlated, i.e., when the
additive Gaussian noise has a non-identity covariance [17].

####  4.1.2 Application to our Problem

Recall that for the sum-SEPCA algorithm, we formed a vector of row-sums. That
is, in the equisigned setting, taking sums across the rows of
$\operatorname{\mathbf{X}}$, we obtain a vector $\operatorname{\mathbf{y}}$
where $y_{i}=\mu_{i}+\sigma z_{i}$, with $\mu_{i}=(\theta
u_{i})\|\operatorname{\mathbf{v}}\|_{1}$: this situation is exactly that of a
sparse mean vector embedded in Gaussian noise. Similarly, taking sums of
squares across the rows of $\operatorname{\mathbf{X}}$ (as in the
$\ell_{2}$-SEPCA algorithm) yields scaled $\chi^{2}_{n}$ distributed random
variables, of which only a few have non-zero non-centrality parameters.

With knowledge of the noise distribution, we may compute the p-values of each
row statistic: these p-values are used to form the Higher Criticism statistic
(12). As in [14], we may adapt the original global testing problem to a
selection problem. For each p-value $p_{(i)}$, we have a value $HC_{p,i}$ of
the higher criticism statistic (the value that is maximized in (12)).
Rejecting each null hypothesis (that the coordinate of the corresponding row
is zero) when $HC_{p,i}$ is larger than the threshold $\sqrt{2\log\log p}$ is
a variable selection procedure. We refer to the procedure based on the sum
statistic as HC-sum-SEPCA and that based on the sum of squares statistic as
HC-$\ell_{2}$-SEPCA. Importantly, we note that the form of the decision
boundary $\rho$ is identical to the global testing case, and that applying
Higher Criticism to our row statistics is a viable global testing procedure
[14].

###  4.2 FDR-SEPCA

In this section, we give an summary of the algorithm for uncorrelated noise
and defer the general case and details to Appendix D. We continue in the same
vein as in the previous section on Higher Criticism.

We note that in the equisigned, rank-$1$ setting, coordinate selection is
equivalent to the estimation of a sparse mean vector. Let
$y_{i}=\mu_{i}+\sigma z_{i}$, where $i\in\\{1,\cdots,p\\}$ and the vector
$\operatorname{\mathbf{z}}$ of the $z_{i}$ is normally distributed with mean
$0$ and covariance $\mathcal{I}_{p}$. The mean vector
$\operatorname{\bm{\mu}}$ of the $\mu_{i}$ is assumed to be sparse; the goal
is to estimate $\operatorname{\bm{\mu}}$. Taking sums across the rows of
$\operatorname{\mathbf{X}}$, we obtain a vector $\operatorname{\mathbf{y}}$
where $y_{i}=\mu_{i}+\sigma z_{i}$, with $\mu_{i}=(\theta
u_{i})\|\operatorname{\mathbf{v}}\|_{1}$. Hence, we are in the same setting as
in the previous section.

The following penalized least squares formulation, taken from [22], yields an
estimator for $\mu$:

$\widehat{\operatorname{\bm{\mu}}}=\arg\min_{\operatorname{\bm{\mu}}}\|\operatorname{\mathbf{y}}-\operatorname{\bm{\mu}}\|_{2}^{2}+\sigma^{2}\text{pen}\left(\|\operatorname{\bm{\mu}}\|_{0}\right),$
(16)

where $\text{pen}(k)$ is defined as

$\text{pen}(k)=\zeta k\left(1+\sqrt{2\log(\nu p/k)}\right)^{2},$ (17)

with $\zeta>1$; we may take $\zeta=1+o(1)$. The parameter $\nu$ is no smaller
than $e$. We define $\|\operatorname{\bm{\mu}}\|_{0}$ to be the number of non-
zero coordinates of $\mu$.

The solution to (16) is given by hard-thresholding. Let $|y|_{(i)}$ be the
$i^{th}$ order statistic of $|y_{i}|$, namely
$|y|_{(1)}\geq\cdots\geq|y|_{(p)}$. Then if

$\widehat{k}=\arg\min_{k\geq
0}\sum_{i>k}|y|_{(i)}^{2}+\sigma^{2}\text{pen}(k),$ (18)

defining

$t_{k}^{2}=\text{pen}(k)-\text{pen}(k-1),$ (19)

the solution is to hard threshold at $t_{\widehat{k}}$.

In this set-up, we have that

$t_{k}\approx\sqrt{\zeta}(1+\sqrt{2\log(\nu p/k)}).$

We provide a precise quantification of $t_{k}$ in Appendix D.

Hence, by computing $t_{k}$ and performing hard thresholding of the row sums,
we can perform coordinate selection. Once again, this procedure replaces the
test statistic/thresholding in Algorithm 1.

##  5 Estimation Error and Smallest Detectable Coordinate

As we will see, our theorems discuss the “detectability” of the coordinates
$u_{i}$ of $\operatorname{\mathbf{u}}$. However, it is common in the sparse
PCA literature to discuss lower bounds for the risk (estimation error) [21, 5,
25]. In what follows, we will show that these two notions are equivalent.

We define the $L^{2}$ estimation error for a principal component estimator as

$L(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}})=\left\|\operatorname{\mathbf{u}}-\text{sign}(\langle\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}}\rangle)\widehat{\operatorname{\mathbf{u}}}\right\|_{2}^{2}.$
(20)

The quantity in (20) is upper bounded by $2$; this bound is attained when
$\operatorname{\mathbf{u}}$ and $\widehat{\operatorname{\mathbf{u}}}$ are unit
norm and mutually orthogonal. Following [5], we want to compute a lower bound
for the maximum expected loss for the $s$-sparse vectors
$\operatorname{\mathbf{u}}$ (in the sense of $\ell_{0}$ sparsity) defined as

$\sup_{\operatorname{\mathbf{u}}\in\mathbb{S}^{p-1}:\|\operatorname{\mathbf{u}}\|_{0}\leq
s}\mathbb{E}L\left(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}}\right),$
(21)

where $\mathbb{S}^{p-1}$ denotes the unit sphere in $\mathbb{R}^{p}$. Let
$\widehat{I}$ be some index set of coordinates selected by an algorithm of the
form given in Algorithm (1). We may take
$\langle\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}}\rangle$
to be non-negative, and decompose the loss as

$\|\operatorname{\mathbf{u}}-\widehat{\operatorname{\mathbf{u}}}\|_{2}^{2}=\underbrace{\|\operatorname{\mathbf{u}}_{\widehat{I}}-\widehat{\operatorname{\mathbf{u}}}\|_{2}^{2}}_{\begin{subarray}{c}\text{Estimation
Error from }\\\ \text{detected
coordinates}\end{subarray}}+\underbrace{\|\operatorname{\mathbf{u}}_{\widehat{I}^{c}}\|_{2}^{2}}_{\begin{subarray}{c}\text{Error
from }\\\ \text{missed
coordinates}\end{subarray}}\geq\|\operatorname{\mathbf{u}}_{\widehat{I}^{c}}\|_{2}^{2}.$
(22)

Equation (22) shows that the loss is lower-bounded by the squared sum of the
missed coordinates. Indeed, it is a natural consequence of the result in [3]
that if the sparsity $s$ grows slower than does $n$, and we have a consistent
estimate of the support of $\operatorname{\mathbf{u}}$, the estimation error
will asymptotically be small. Essentially, we are estimating the singular
vectors of an $s\times n$ matrix instead of a $p\times n$ matrix, so that if
the ratio $s/n$ has limit zero, our estimates will be consistent (see (3) and
[3]). This suggests the following strategy for lower-bounding (21): we want to
construct a non-trivial ‘worst-case’ sparse vector. That is, we want a vector
$\operatorname{\mathbf{u}}$ that has a non-trivial loss (less than $2$), is
sparse (fewer than $s$ non-zero coordinates), and has maximal error from
missed coordinates. To ensure a non-trivial loss, we set the first coordinate
$u_{1}$ to be large, i.e., $u_{1}=\sqrt{1-r^{2}}$, where $r=o(1)$. To ensure
sparsity, we set $u_{2},\cdots,u_{m+1}$ to be non-zero for some $m\leq s-1$,
with the subsequent coordinates of $u$ set to $0$. Then, the expected loss has
the lower bound

$\begin{split}\mathbb{E}L(\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}})&\geq\sum_{k=1}^{p}|u_{k}|^{2}\mathbb{P}\left(\text{Not
Selecting Coordinate k}\right)\\\
&\geq\sum_{k=2}^{m+1}|u_{k}|^{2}\mathbb{P}\left(\text{Not Selecting Coordinate
k}\right),\end{split}$ (23)

since $u_{1}$ is detected with probability approaching $1$ and $u_{k}$ is zero
for $k>m+1$. Now, let $u_{2}$ through $u_{m+1}$ all have value ${r}/\sqrt{m}$,
so that we may simplify the lower bound to

$\mathbb{E}L(\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}})\geq
r^{2}\mathbb{P}\left(\text{Not Selecting Coordinate k}\right).$ (24)

If coordinates of size $r/\sqrt{m}$ are not detected with a probability
approaching $1$, $r^{2}$ is a lower-bound on the risk. This construction shows
that specifying the sizes of coordinates that are not detected with
probability approaching $1$ is equivalent to specifying a worst-case risk
lower bound. Note that the value of $r^{2}$ depends on the specific algorithm
and estimator, and that this is not a general or universal bound. Rather, the
purpose of this construction is to show the equivalence between the two
perspectives (a lower bound and detectable coordinate size).

Consequently, in what follows we focus on the smallest detectable and largest
undetectable coordinates because they directly shed light on the attainable
estimation error. The details of the risk calculations and extensions to
approximate sparsity are deferred to Appendix C, where we summarize our
findings in Theorem 3.

##  6 Main Results

The following theorem characterizes consistent support recovery conditions.
These results are the analogue of the ‘sparsistency’ guarantees found in the
LASSO and $\ell_{1}$-norm minimization literature [33]. Throughout,
$\widehat{I}$ denotes the set of coordinates selected by the coordinate
selection scheme.

######  Theorem 1.

For the model specified in (1) and (2) and the algorithms specified in Table
1, assume that $p(n),n\rightarrow\infty$, $s(n)/n\rightarrow 0$, and $\log
p(n)=o(n)$. Let $\epsilon\in(0,1)$. We have that

  1. a.

For $i\in I^{c}$,

$\phantom{\qquad}\max_{i\in
I^{c}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow 0,$

  2. b.

For $i\in I$,

$\min_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|>\beta_{crit}(1+\epsilon)\end{subarray}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow
1,$ $\max_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|<\beta_{crit}(1-\epsilon)\end{subarray}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow
0.$

Here

$\beta_{crit}=\left\\{\begin{array}[]{ll}\sigma C_{U}\frac{\sqrt{\log
p}}{\left|\sum_{k}v_{k}\right|}&\text{for sum-SEPCA,}\\\
\sigma\sqrt{C_{2}}\sqrt{\frac{\log ep}{\sqrt{n}}}&\text{for
$\ell_{2}$-SEPCA,}\\\ \sigma t_{\ell_{1}}&\text{for
$\ell_{1}$-SEPCA,}\\\\[4.26773pt] \end{array}\right.$ (25)

and $t_{\ell_{1}}$ satisfies the relation

$\displaystyle\left(\sqrt{\frac{2}{\pi}}+C_{1}\frac{\log ep}{\sqrt{n}}\right)$
$\displaystyle=\frac{1}{n}\sqrt{\frac{2}{\pi}}[\sum_{k}\exp\left(-\left(\sqrt{n}\frac{(t_{\ell_{1}})v_{k}}{\sqrt{2}}\right)^{2}\right)+$
$\displaystyle\sqrt{\pi}\sum_{k}\left(\sqrt{n}\frac{(t_{\ell_{1}})v_{k}}{\sqrt{2}}\right)\text{Erf}\left(\sqrt{n}\frac{(t_{\ell_{1}})v_{k}}{\sqrt{2}}\right)].$

We defer the proof to Appendix A.

Theorem 1 identifies a phase transition in the ability of the algorithms to
accurately estimate the support of $\operatorname{\mathbf{u}}$. Note that the
analysis brings into sharp focus the dependence of $\beta_{crit}$ on
$\operatorname{\mathbf{v}}$ for the $\ell_{1}$\- and sum-SEPCA algorithms, but
not the $\ell_{2}$-SEPCA algorithm. Consequently, we can expect the algorithms
to perform differently depending on the structure of the underlying
$\operatorname{\mathbf{v}}$. It is important to note that the sparsity $s$ of
$\operatorname{\mathbf{u}}$ is not a parameter in the thresholds and results.

It is also important to note that $\ell_{2}$-SEPCA and $\ell_{1}$-SEPCA do not
rely on the equisigned character of $\operatorname{\mathbf{v}}$. However, it
is clear that the sum-SEPCA algorithm explicitly depends on the equisigned
assumption.

###  6.1 Hamming Loss

It is also possible to state the above results in terms of the Hamming loss
for the support of $\operatorname{\mathbf{u}}$, and prove consistency of the
coordinate selection scheme by assuming that all the nonzero coordinates of
$\operatorname{\mathbf{u}}$ lie above a critical threshold. A detailed
decision-theoretic analysis of variable selection under a sequence model with
_i.i.d._ noise and with respect to the Hamming loss has recently been carried
out by [10]. Recall that the Hamming loss measures the number of elements in
two sets that are different, so that here the loss between the true support
$I$ and the estimated support $\widehat{I}$ would be the size of the symmetric
set difference of $I$ and $\widehat{I}$. Let $d_{H}\left(I,\widehat{I}\right)$
denote the Hamming loss and assume that whatever algorithm we are using has a
threshold $\beta_{crit}$. Then, for any $\epsilon\in(0,1)$, we may write

$\displaystyle\operatorname{\mathbb{E}}d_{H}\left(I,\widehat{I}\right)$
$\displaystyle=$ $\displaystyle\sum_{i\in I:|\theta
u_{i}|>\beta_{crit}(1+\epsilon)}\operatorname{\mathbb{P}}\left(i\notin\widehat{I}\right)+\sum_{i\in
I:|\theta
u_{i}|<\beta_{crit}(1-\epsilon)}\operatorname{\mathbb{P}}\left(i\notin\widehat{I}\right)$
(26) $\displaystyle+\sum_{i\in I:(1-\epsilon)\leq|\theta
u_{i}|/\beta_{crit}\leq(1+\epsilon)}\operatorname{\mathbb{P}}\left(i\notin\widehat{I}\right)+\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right).$

We can then restate the results on coordinate selection in terms of the
Hamming loss under a more restricted setting that assumes an exact form of
sparsity of the vector $\operatorname{\mathbf{u}}$.

######  Corollary 1.

For the model specified in (1) and (2), and an algorithm specified in Table 1,
assume that the conditions of Theorem 1 hold. Let the support $I$ of
$\operatorname{\mathbf{u}}$ be estimated by $\widehat{I}$. Moreover, assume
that the algorithm has a threshold $\beta_{crit}$ (given in (25)) such that
for a small, fixed $\epsilon_{0}>0$, the set

$I_{0}:=\left\\{j:|\theta u_{j}|>\beta_{crit}(1+\epsilon_{0})\right\\}$ (27)

equals the set $I$. Then the expected Hamming loss satisfies

$\operatorname{\mathbb{E}}d_{H}\left(I,\widehat{I}\right)/|I|\rightarrow 0.$
(28)

The proof of the corollary, given in Appendix A.1, follows from applying
Theorem 1, with a more detailed enumeration of the sets and the inclusion
probabilities, to each term of (26).

###  6.2 FDR-Based Algorithms

We may summarize the coordinate selection properties of the FDR refinements as
follows:

######  Theorem 2.

For the model specified in (1) and (2) and the three FDR-controlling
algorithms summarized in Algorithm 2, assume that $p(n),n\rightarrow\infty$,
$s(n)/n\rightarrow 0$, and $\log p(n)=o(n)$. Let $\epsilon\in(0,1)$. We have
that

  1. a.

For all three algorithms and $i\in I^{c}$,

$\phantom{\qquad}\max_{i\in
I^{c}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow 0,$

  2. b.

For the Higher Criticism-based algorithms and $i\in I$,

$\min_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|>\beta_{crit}(1+\epsilon)\end{subarray}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow
1,$ $\max_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|<\beta_{crit}(1-\epsilon)\end{subarray}}~{}\mathbb{P}\left(i\in\widehat{I}\right)~{}\rightarrow
0.$

  3. c.

For the FDR-SEPCA algorithm, uniformly over $i\in I$,
if $|\theta u_{i}|>\beta_{crit}(1+\epsilon)$, coordinate $i$ is selected;
if $|\theta u_{i}|<\beta_{crit}(1-\epsilon)$, coordinate $i$ is not selected
with probability tending to 1.

Here

$\displaystyle\beta_{crit}=\left\\{\begin{array}[]{ll}\sigma\sqrt{\rho(\beta)}\frac{\sqrt{2\log
p}}{\|\operatorname{\mathbf{v}}\|_{1}}&\text{for HC-sum-SEPCA,}\\\
\sigma\rho(\beta)\frac{2\log p}{\sqrt{n}}&\text{for HC-$\ell_{2}$-SEPCA,}\\\
\sigma\left(1-o(1)\right)\sqrt{\zeta}\frac{1+\sqrt{2\log(\nu
p/\widehat{k})}}{\|\operatorname{\mathbf{v}}\|_{1}}&\text{for FDR-SEPCA,}\\\
\end{array}\right.$ (32)

where $\zeta>1$, $\nu>e$, and the FDR-SEPCA algorithm detects $\widehat{k}$
coordinates.

We defer the proof to Appendix B.

Once again, we see that the structure of the underlying
$\operatorname{\mathbf{v}}$ plays a role in the performance of the sum-based
algorithms, but not for the $\ell_{2}$-based HC-$\ell_{2}$-SEPCA algorithm.
Unlike in the FWER-controlling cases, the sparsity of
$\operatorname{\mathbf{u}}$ plays a (small) role here, via the constant
$\rho\left(\beta\right)$ for the Higher Criticism-based methods and via
$\widehat{k}$ for FDR-SEPCA. Moreover, $\ell_{2}$-HC-SEPCA, like
$\ell_{2}$-SEPCA, does not make use of the equisigned nature of
$\operatorname{\mathbf{v}}$.

###  6.3 Higher Ranks

In this work, we restrict our focus to the rank-$1$, equisigned setting. A
natural question is are our results extensible to the higher rank setting?

The first point is concerned with the right singular vectors. To preserve
orthogonality, we would need equisigned right-singular vectors
$\operatorname{\mathbf{v}}_{i}$ with disjoint supports. The second point is
concerned with the left singular vectors. Our algorithms are based on
thresholding row-statistics: it is possible that the union of supports of
several sparse vectors is a relatively large set. The FWER-controlling
algorithms (by design) are not sensitive to the increased supports, but the
FDR-controlling algorithms are sensitive to this. Indeed, the decision
boundaries for the FDR algorithms explicitly depend on the sparsity levels.
Third point, once again, is concerned with the left singular vectors. It is
possible that a sum-based statistic suffers from cancellations that decrease
the size of the row-statistic. For example, in a rank-$2$ setting, if
$\operatorname{\mathbf{u}}_{1}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1&0&\cdots&0\end{bmatrix}^{T}$
and
$\operatorname{\mathbf{u}}_{2}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-1&0&\cdots&0\end{bmatrix}^{T}$
and $\left\|\operatorname{\mathbf{v}}_{1}\right\|_{1}$ and
$\left\|\operatorname{\mathbf{v}}_{2}\right\|_{1}$ have similar values and are
both non-negative, the row-sum of the second row will be small. Note, however,
that the $\ell_{2}$-norm based methods do not suffer from this issue.

##  7 Simulations

To illustrate the relative powers of the six algorithms, we compute the
theoretical limits on the sizes of detectable coordinates as a function of
$n$. We use a unit-norm, equisigned $\operatorname{\mathbf{v}}$ such that

$v_{k}\propto\exp\left(-5\frac{k}{n}\right)\left|\sin\left(4\frac{k}{n}\right)\right|\text{
for }1\leq k\leq n.$ (33)

This choice of $\operatorname{\mathbf{v}}$ has a ‘rise and fall’ sort of
behavior, and is motivated by physical signals, e.g., chemical reactions or
nerve signals in the brain. The value of $\beta_{crit}$ is shown in Figure 1;
for this choice of $\operatorname{\mathbf{v}}$, it is clear that the sum-SEPCA
dramaticaly outperforms the other SEPCA variants in terms of size of the
smallest detectable component. The FDR-SEPCA algorithm has similar performance
to sum-SEPCA, and the HC-sum-SEPCA algorithm has the strongest performance.

In Figure 2, we plot the estimation error as a function of $n$ and $\theta$
for all six algorithms. We also include results for the SVD and competing
algorithms TPower [41] and ITSPCA [25]. In the simulations, we fix $p=1000$
and vary $n$, since the dependence in $p$ in the thresholds is logarithmic,
whereas that in $n$ is not. The left singular vector
$\operatorname{\mathbf{u}}$ is chosen to be the vector with $1$ in the first
coordinate and $0$ elsewhere. We fix the noise variance $\sigma^{2}$ at $1$,
so that $\theta^{2}$ is the eigen-SNR. The results should be interpreted as
follows. For the particular $\operatorname{\mathbf{v}}$ chosen here, we expect
HC-sum-SEPCA to have the lowest detectable limit, and $\ell_{1}$-SEPCA to have
the largest. This behavior is confirmed. Moreover, the sum-based algorithms
offer a slight strengthening of both ITSPCA and TPower. Importantly, note that
the sum-based algorithms explicitly take advantage of the equisigned nature of
$\operatorname{\mathbf{v}}$: that is, algorithms that explicitly use the
equisigned property outperform algorithms that do not (the $\ell_{2}$ and
$\ell_{1}$ algorithms, as well as ITSPCA and TPower).

Figure 1: This plot shows $\beta_{crit}$ for all six algorithms for the
$\operatorname{\mathbf{v}}$ described in (33).

Figure 2: The plots show the empirical estimation error for all six algorithms
for the $\operatorname{\mathbf{u}}$ with one non-zero coordinate and the
$\operatorname{\mathbf{v}}$ described in (33). We include results from TPower,
ITSPCA and the SVD for comparison.

We repeat our simulations for a
$\operatorname{\mathbf{u}}\in\operatorname{\mathbb{R}}^{p}$ with $\sqrt{p}$
non-zero coordinates (of equal size) and the same $\operatorname{\mathbf{v}}$,
as seen in Figure 3. We find similar conclusions as in Figure 2, where the
sum-based algorithms offer a strengthening over ITSPCA and TPower; the
$\ell_{2}$-norm based algorithms do not perform well. Note that a sparsity of
$\sqrt{p}$ is at the limit/valid edge for the higher criticism-based methods,
but that these methods still perform well.

Figure 3: The plots show the empirical estimation error for all six algorithms
for the $\operatorname{\mathbf{u}}$ with $\sqrt{p}$ non-zero coordinates and
the $\operatorname{\mathbf{v}}$ described in (33). We include results from
TPower, ITSPCA and the SVD for comparison.

###  7.1 Comments on the FDR-controlling procedures

The Higher Criticism for the $\chi_{n}^{2}$-variates ‘pushes back’ the phase
transition between detecting nothing and something to a lower value of
$\theta$ relative to the $\ell_{2}$-SEPCA algorithm, but is still less
powerful than any of the sum-based algorithms. Moreover, even above the phase
transition, the $\ell_{2}$-SEPCA algorithm may be preferable, as the error is
increased by unacceptably many false positives.

The Higher Criticism procedure for the sum statistic has the lowest phase
transition point and hence the highest power. Its transition is more gradual
than the penalized FDR thresholding procedure and sum-SEPCA, which have
roughly the same performance in this simulation.

###  7.2 An example where $\ell_{2}$-based algorithms outperform sum-based
algorithms

Sum-SEPCA has a $\beta_{crit}$ that depends on $\operatorname{\mathbf{v}}$.
Looking at the form in (25), if $\|\operatorname{\mathbf{v}}\|_{1}$ is smaller
than $n^{1/4}$, we would expect $\ell_{2}$-SEPCA to detect a smaller
coordinate size. Vectors with smaller coordinates have a smaller
$\ell_{1}$-norm, i.e., one that is closer to their $\ell_{2}$-norm. Hence, if
we choose

$v_{k}\propto\frac{1}{k^{2}}\text{ for }1\leq k\leq n,$ (34)

we expect sum-SEPCA to have worse performance relative to $\ell_{2}$-SEPCA.
Figures 4 and 5 confirm this expectation. The FDR refinements perform poorly.
It should be noted, however, that TPower and ITSPCA retain their performance.
This choice of $\operatorname{\mathbf{v}}$ effectively corresponds to a very
small value of $n$: the majority of coordinates are tiny in size and buried
beneath noise regardless of the value of $\theta$. If we ‘corrected’ the
scenario and used a smaller $n$ and a subset of $\operatorname{\mathbf{v}}$,
we would be in a situation closer to that given in (33).

Figure 4: This plot shows $\beta_{crit}$ for all six algorithms for the
$\operatorname{\mathbf{v}}$ described in (34).

Figure 5: The plots show the empirical estimation error for all six algorithms
for the $\operatorname{\mathbf{u}}$ and $\operatorname{\mathbf{v}}$ described
in (34). We include results from TPower, ITSPCA and the SVD for comparison.

###  7.3 A video data example

We conclude our sequence of examples with a real data study. This example is
motivated by the problem of foreground-background separation in videos.
Consider a grayscale video of stars twinkling against a black background [36].
Our goal is to estimate the locations of the stars: by reshaping the video, we
may treat each frame as a vector and hence treat the video as a sparse matrix.
Only a few locations have a star and are hence non-zero. The scale of the
video pixels is between $0$ and $255$. We examine the top-left $72\times 64$
pixels for $89$ frames, as shown in Figure 6(a). In Figure 6(b), we plot the
singular values of the video matrix. The first singular value stands out
strongly against the rest, and at most two more singular values are well-
separated from the bulk. This structure suggests that our rank-$1$ based
approach is well suited to this problem.

|
---|---

(a) The image shows the mean intensity of pixels from the top-left $72\times
64$ pixels for $89$ frames. White indicates the presence of a star.

(b) The plots shows the singular values of the video data. The spacing suggest
a low-rank-plus-noise structure.

Table 2: Video Example Figures

We add Gaussian noise of variance $\sigma^{2}$ and study the True Positive
Rates (TPR) and False Discovery Rates (FDR) across all algorithms and across
different values of $\sigma$. In Figure (6), we show the results of our
simulations. In terms of the TPR, everything other than the SVD has a similar
performance, while the test-statistic SEPCA-based algorithms enjoy the best
performance in terms of the FDR. In Figure 7 we zoom in on the top-right three
stars and show how the algorithms perform as noise increases. Here, we see
that the behavior alluded to in the TPR/FDR results actually occurs in the
video.

Figure 6: The left plot shows the True Positive Rate of the various algorithms
as a function of the noise level $\sigma$. The right plot shows the False
Discovery Rates.

(a) $\sigma=0$

(b) $\sigma=0.3$

(c) $\sigma=15$

(d) $\sigma=200$

Figure 7: A zoomed-in view of the three top-right stars in the video example.
White indicates a false negative (missed star), Red a false positive (a
guessed pixel where there was nothing), and Blue a true positive (correctly
identified pixel).

##  8 A geometric view: which algorithm to use?

We have stated detectability results for each algorithm in Section 6 and
provided a numerical verification and comparison in Section 7. In this
section, we wish to analytically compare the algorithms. In particular, we
have seen that the right singular vector $\operatorname{\mathbf{v}}$ plays a
critical role in the detectability and estimability of
$\operatorname{\mathbf{u}}$, and we will characterize this behavior carefully.

In this section, will use the following notational convenience: we absorb
$(\theta u_{i})$ into $\operatorname{\mathbf{v}}\in\mathbb{R}^{n}$, and write
the detectability of coordinates in terms of $\operatorname{\mathbf{v}}$. That
is, if $\operatorname{\mathbf{v}}^{T}$ is a row of
$\operatorname{\mathbf{X}}$, we specify when that row is selected. Moreover,
we take $\sigma=1$ for simplicity.

There are two ‘classes’ of detectability: in terms of
$\|\operatorname{\mathbf{v}}\|_{1}$ and in terms of
$\|\operatorname{\mathbf{v}}\|_{2}$. The sum-, HC-sum, and FDR-SEPCA
algorithms select a coordinate if
$\left|\sum_{k}v_{k}\right|=\|\operatorname{\mathbf{v}}\|_{1}$ is large enough
for a $\operatorname{\mathbf{v}}$ in the orthant with all non-negative or all
non-positive coordinates. Geometrically, the vector
$\operatorname{\mathbf{v}}$ is selected if it is ‘outside’ a hyperplane with a
normal vector proportional to the vector of all $1$s. The $\ell_{1}$-SEPCA
algorithm is similar, as it selects a coordinate when
$\|\operatorname{\mathbf{v}}\|_{1}$ is large enough, or if
$\operatorname{\mathbf{v}}$ lies outside an $\ell_{1}$-ball of some radius.
The connection between the previous three algorithms and $\ell_{1}$-SEPCA
comes from noting that the faces of an $\ell_{1}$-ball are sections of
hyperplanes with normal vectors proportional to a vector of $\pm 1$s. Finally,
the $\ell_{2}$\- and HC-$\ell_{2}$-SEPCA algorithms select a coordinate when
$\|\operatorname{\mathbf{v}}\|_{2}$ is large enough. I.e., when
$\operatorname{\mathbf{v}}$ lies outside some $\ell_{2}$-ball.

Our goal in this section is to derive comparisons between the six algorithms.
Specifically, for a given vector $\operatorname{\mathbf{v}}$, which algorithm
will have the greatest detection ability (we are, for the moment, only
concerned with maximizing power)? Note that when $\operatorname{\mathbf{v}}$
has a large norm, it does not matter which algorithm is used. Questions only
arise when $\|\operatorname{\mathbf{v}}\|_{1}$ or
$\|\operatorname{\mathbf{v}}\|_{2}$ are relatively small and are close to the
thresholds.

###  8.1 Intersection of a hyperplane and a hypersphere

We may think of the $\ell_{1}$ ball as a hyperplane when restricted to a
single orthant. If a hypersphere of radius $r$ intersects a hyperplane with a
normal vector proportional to the vector of all $\pm 1$s and minimum distance
to the origin of $r-h$, a hyperspherical cap of height $h$ is formed: see
Figure 8 for a simple illustration. Geometrically, a right triangle is formed,
with hypotenuse $r$ and leg $r-h$. Hence, the angle between the center of the
cap and the edge is:

$\theta_{lim}=\cos^{-1}\frac{r-h}{r}.$ (35)

It is sufficient to guarantee that

$0\leq\frac{r-h}{r}\leq 1$

for the hyperspherical cap to exist. Moreover, a vector
$\operatorname{\mathbf{v}}$ has a direction contained inside the cap when the
angle between $\operatorname{\mathbf{v}}$ and the vector of $\pm 1$ in the
orthant containing $\operatorname{\mathbf{v}}$ is smaller than $\theta_{lim}$.
In other words, defining the angle for a vector $\operatorname{\mathbf{v}}$ as

$\theta(\operatorname{\mathbf{v}})=\cos^{-1}\frac{\|\operatorname{\mathbf{v}}\|_{1}}{\|\operatorname{\mathbf{v}}\|_{2}\sqrt{n}},$
(36)

we need $\theta(\operatorname{\mathbf{v}})\leq\theta_{lim}$.

Figure 8: A spherical cap in $\mathbb{R}^{2}$

###  8.2 Comparison: $\ell_{2}$-based versus sum-based statistics

We begin with a summary of the performance of each individual algorithm in
Table 3. We first compare $\ell_{2}$-SEPCA and then compare
HC-$\ell_{2}$-SEPCA with sum-, HC-sum, FDR-SEPCA in Tables 4 and 5. In our
comparisons, we consider when the hyperspherical cap exists and give the angle
of the cap. These are routine calculations, so we omit the details. We also
omit $\ell_{1}$-SEPCA from our comparisons, as we lack a closed-form
expression for $t_{\ell_{1}}$.

Note that the existence of this cap is a proxy for the equivalent statement
that there exist vectors for which the sum-based algorithms are more powerful
than the $\ell_{2}$-based algorithm. While this existence is not the same as
attributing uniformly greater power to the sum-based algorithms relative to
the $\ell_{2}$-based algorithm, the cap _not_ existing is equivalent to the
$\ell_{2}$-based algorithm having uniformly greater power.

Essentially, we observe that for $n$ and $p$ sufficiently large, the cap will
exist. Moreover, for $\operatorname{\mathbf{v}}$ that is sufficiently dense
($\|\operatorname{\mathbf{v}}\|_{1}$ is sufficiently large), $\theta_{lim}$
will lie inside the cap. Hence, in these situations we would prefer a sum-
based algorithm over an $\ell_{2}$-based algorithm.

Table 3: A summary of the six algorithms.

Algorithm | Threshold | Geometric Quantity
---|---|---
sum | $\|\operatorname{\mathbf{v}}\|_{1}\geq C_{U}\sqrt{\log p}$ | $r-h=C_{U}\sqrt{\frac{\log p}{n}}$
HC-sum | $\|\operatorname{\mathbf{v}}\|_{1}\geq\sqrt{2\rho(\beta)\log p}$ | $r-h=\sqrt{2\rho(\beta)\frac{\log p}{n}}$
FDR | $\begin{aligned} &\|\operatorname{\mathbf{v}}\|_{1}\geq\left(1-o(1)\right)\sqrt{\zeta}\\\ &\left(1+\sqrt{2\log\left(\nu p/\widehat{k}\right)}\right)\end{aligned}$ | $\begin{aligned} r&-h=\left(1-o(1)\right)\sqrt{\zeta}\\\ &\frac{1+\sqrt{2\log\left(\nu p/\widehat{k}\right)}}{\sqrt{n}}\end{aligned}$
$\ell_{1}$ | $\|\operatorname{\mathbf{v}}\|_{1}\geq t_{\ell_{1}}$ | $r-h=t_{\ell_{1}}$
$\ell_{2}$ | $\|\operatorname{\mathbf{v}}\|_{2}\geq\sqrt{C_{2}}\sqrt{\frac{\log ep}{\sqrt{n}}}$ | $r=\sqrt{C_{2}}\sqrt{\frac{\log ep}{\sqrt{n}}}$
HC-$\ell_{2}$ | $\|\operatorname{\mathbf{v}}\|_{2}\geq 2\rho(\beta)\frac{2\log p}{\sqrt{n}}$ | $r=2\rho(\beta)\frac{2\log p}{\sqrt{n}}$

Table 4: The relative performance of $\ell_{2}$-SEPCA.

Algorithm | $\cos\theta_{lim}$ | Cap exists if
---|---|---
sum | $C_{U}/\sqrt{C_{2}}\sqrt{\frac{\log p}{(1+\log p)\sqrt{n}}}$ | $n\geq C_{U}^{4}/C_{2}^{2},p\geq 1$
HC-sum | $\sqrt{2\rho(\beta)/C_{2}}\sqrt{\frac{\log p}{(1+\log p)\sqrt{n}}}$ | $n\geq 1,p\geq 1$
FDR-sum | $\frac{1+\sqrt{2\log\nu p/\widehat{k}}}{\sqrt{C_{2}\sqrt{n}\log\nu p}}$ | $p\geq 11,n\geq 1$

Table 5: The relative performance of HC-$\ell_{2}$-SEPCA.

Algorithm | $\cos\theta_{lim}$ | Cap exists if
---|---|---
sum | $\frac{C_{U}}{2\rho(\beta)\sqrt{\log p}}$ | $p\geq\exp\left(\frac{C_{U}^{2}}{4\rho(\beta)^{2}}\right)$
HC-sum | $\left[2\rho(\beta)\log p\right]^{-1}$ | $p\geq\exp\left(\frac{1}{2\rho(\beta)}\right)$
FDR-sum | $\frac{1+\sqrt{2\log\nu p/\widehat{k}}}{2\rho(\beta)\log p}$ | $\begin{aligned} &\log p\geq\\\ &\frac{1}{4\rho(\beta)^{2}}\bigl{(}1+2\rho(\beta)\\\ &+\sqrt{8\rho(\beta)^{2}+4\rho(\beta)+1}\bigr{)}\end{aligned}$

###  8.3 HC-$\ell_{2}$-SEPCA versus $\ell_{2}$-SEPCA

Now, we consider when HC-$\ell_{2}$-SEPCA is more powerful than
$\ell_{2}$-SEPCA. The ratio of the radii is given by

$\frac{2\rho(\beta)}{\sqrt{C_{2}}n^{1/4}}\frac{\log p}{\sqrt{1+\log p}}.$ (37)

If this ratio is smaller than $1$, HC-$\ell_{2}$-SEPCA is more powerful than
$\ell_{2}$-SEPCA. Note that the quantity

$\frac{2\sqrt{2}}{\sqrt{C_{2}}}\sqrt{\frac{\log p}{\sqrt{n}}}$

is an upper bound for (37), so that if

$\frac{\log p}{\sqrt{n}}<\frac{e}{4\sqrt{2}},$

the original ratio is smaller than $1$ and HC-$\ell_{2}$-SEPCA is preferable
to $\ell_{2}$-SEPCA.

###  8.4 Comparing the sum-based algorithms

Finally, we compare sum-, HC-sum-, and FDR-SEPCA. First, the ratio of the
thresholds for HC-sum- and sum-SEPCA is

$\frac{\sqrt{2\rho(\beta)}}{C_{U}}.$ (38)

Noting that $\rho(\beta)\leq 1$ and that $C_{U}\geq\sqrt{2}+1/3\sqrt{2}$, it
is clear that this ratio is always smaller than $1$ so that HC-sum-SEPCA is a
strict improvement on sum-SEPCA.

Next, we compute the ratio of the thresholds for FDR- and sum-SEPCA:

$\frac{1+\sqrt{2\log\nu p/\widehat{k}}}{C_{U}\sqrt{\log p}}.$ (39)

Using the lower bound on $C_{U}$, we find that if $\widehat{k}\geq 11$ (and
$p\geq\widehat{k}$, naturally), FDR-SEPCA is always more powerful than sum-
SEPCA. For smaller values of $\widehat{k}$, for sufficiently large values of
$p$, the ratio will be smaller than $1$.

Lastly, we compare FDR-SEPCA to HC-sum-SEPCA, wherein the ratio of the
thresholds is (FDR to HC-sum):

$\frac{1+\sqrt{2\log\nu p/\widehat{k}}}{\sqrt{2\rho(\beta)}\sqrt{\log p}}.$
(40)

Because of involvement of $\rho(\beta)$, this quantity is hard to analyze. If
in an oracle manner, FDR-SEPCA obtained $\widehat{k}$ correctly as
$p^{1-\beta}$, we would find that this ratio is always larger than $1$ for
$p>1$. That is if $\widehat{k}$ assumes the the correct value, HC-sum-SEPCA is
more powerful than FDR-SEPCA. Alternatively, we can note that
$\rho(\beta)\in(0,1]$ and ask when the ratio is larger than $1$. Based on the
ratio above, we can see that in the following scenarios

$\displaystyle\left\\{\begin{array}[]{lcl}\widehat{k}=1&\text{and}&p>1,\\\
2\leq\widehat{k}\leq 18&\text{and}&p\geq\widehat{k}\,\text{(always)},\\\
\widehat{k}\geq 19&\text{and}&\log
p\geq\frac{1}{8}\left(4(\log\widehat{k})^{2}-4\log\widehat{k}+1\right),\end{array}\right.$
(44)

HC-sum-SEPCA is more powerful than FDR-SEPCA.

To summarize, we prefer the FDR-controlling alternatives to sum-SEPCA, but
depending on the output of FDR-SEPCA, HC-sum-SEPCA may be more powerful.
However, as the simulations in Section 7 revealed (see Figure 2), the number
of false positives with HC-sum-SEPCA may be higher than with FDR-SEPCA.

###  8.5 Overall Message

We have seen that for $n$ and $p$ sufficiently large and
$\operatorname{\mathbf{v}}$ that is sufficiently dense (in the sense of
$\|\operatorname{\mathbf{v}}\|_{1}$ being large), a sum-based statistic and
algorithm leads to better performance. This is expected behavior, as by using
a sum-based method, we are taking advantage of the equisigned nature of
$\operatorname{\mathbf{v}}$. Moreover, within the class of sum-based
algorithms, controlling the FDR leads to greater power, as expected. It is
difficult to clearly identity which of HC-sum- and FDR-SEPCA will have the
greatest power, and the end result may come down to a practitioner’s tolerance
for false discoveries.

##  9 Conclusions

We have considered the setting where the left singular vector of the
underlying rank one signal matrix plus noise data matrix is assumed to be
sparse and the right singular vector is assumed to be equisigned. We have
proposed six different SEPCA algorithms for estimating the sparse principal
component based on different decision statistics and provided sparsistency
conditions for the same. Our analysis reveals conditions where a coordinate
selection scheme based on a sum-based decision statistic outperforms schemes
that utilize the $\ell_{1}$ and $\ell_{2}$ decision statistics. Thereby, the
proposed algorithm outperforms known schemes such as diagonal thresholded PCA
[21] in terms of estimation of the singular vectors associated with the rank-1
component. We have derived lower bounds on the size of detectable coordinates
of the principal left singular vector, utilized these lower bounds to derive
lower bounds on the worst-case risk and verified our findings with numerical
simulations. Finally, we have discussed the results of our simulations
analytically, by providing a geometric interpretation of the differences in
power among the algorithms.

We note that while we have stated our results for Gaussian noise with identity
covariance, we can extend the FWER-controlling results to any log-concave
noise distribution, and the FDR-controlling procedures to Gaussian noise with
certain non-identity covariances. Additionally, another way to view this work
is that it proposes a two-stage procedure/framework for sparse PCA based
around hypothesis testing of statistics associated to each row. Some natural
extensions would be the inclusion or consideration of other testing
frameworks, e.g., that in [28], where knowledge of the size of coordinates is
taken into account.

##  Appendix A Proof of Theorem 1

  1. a.

Note that

$\mathbb{P}\left(T_{i}\geq\tau\right)\leq\mathbb{P}\left(\max_{j\in
I^{c}}T_{j}\geq\tau\right)$

for $i\in I^{c}$. Taking the maximum over the left-hand side and noting that
the right-hand side has limit zero yields the result. This follows from (5). ∎

  2. b.

We consider when true positives occur with probability approaching $1$. We
want to find the smallest coordinate $(\theta u_{i})$ such that the following
probability approaches $1$:

$\mathbb{P}\left(T_{i}>\tau_{n,p}\right)=\mathbb{P}\left(\frac{T_{i}-\mathbb{E}T_{i}}{\sqrt{\text{Var
}T_{i}}}>\frac{\tau_{n,p}-\mathbb{E}T_{i}}{\sqrt{\text{Var }T_{i}}}\right).$
(45)

Note that if $(\tau_{n,p}-\mathbb{E}T_{i})$ is negative and not tending to
zero as $n$ grows, and if the variance of $T_{i}$ decays to zero as $n$ grows,
the quantity

$\frac{\tau_{n,p}-\mathbb{E}T_{i}}{\sqrt{\text{Var }T_{i}}}$ (46)

tends toward negative infinity. Hence, we will specify conditions so that
$\text{Var }T_{i}$ decays to zero as $n$ grows and then compute when a
coordinate is detectable by considering when $\tau_{n,p}$ is strictly less
than $\mathbb{E}T_{i}$. For brevity, we omit the computations in solving
$\tau_{n,p}<\mathbb{E}T_{i}$ for $|\theta u_{i}|$ and present verifications
that the variance of $T_{i}$ has limit $0$. These results show that above the
_decision boundary_ , we have uniform detection.

In sum-SEPCA, $T_{i}$ is a Gaussian random variable with mean $\frac{(\theta
u_{i})}{\sqrt{n}}\sum_{k}v_{k}$ and variance $\frac{\sigma^{2}}{n}$. Since
$\sigma$ does not grow with $n$, $\text{Var }T_{i}$ always decays to zero.

In $\ell_{2}$-SEPCA, $T_{i}$ has

$\mathbb{E}T_{i}=\left(\theta u_{i}\right)^{2}+\sigma^{2}$

and

$\text{Var }T_{i}=\frac{2\sigma^{2}}{n}\left(\sigma^{2}+2\left(\theta
u_{i}\right)^{2}\right).$

Since $\sigma$ and $\theta$ are fixed, the variance always decays to $0$.

Let $x_{i,k}=\left(\sqrt{n}\frac{v_{k}\left(\theta
u_{i}\right)}{\sigma}\right)$. In $\ell_{1}$-SEPCA, $T_{i}$ has

$\displaystyle\text{Var }T_{i}$
$\displaystyle=\frac{\sigma^{2}}{n^{2}}\sum_{k}x_{i,k}^{2}\left(1-\left(\text{Erf
}\left(\frac{x_{i,k}}{\sqrt{2}}\right)\right)^{2}\right)$
$\displaystyle+\frac{\sigma^{2}}{n}\left(1-\frac{2}{n\pi}\sum_{k}\exp\left(-x_{i,k}^{2}\right)\right)$
$\displaystyle-2\sqrt{\frac{2}{\pi}}\frac{\sigma^{2}}{n^{2}}\sum_{k}x_{i,k}\exp\left(-\frac{x_{i,k}^{2}}{2}\right)\text{Erf
}\left(\frac{x_{i,k}}{\sqrt{2}}\right),$

which is less than or equal to

$\frac{\left(\theta
u_{i}\right)^{2}}{n}\sum_{k}v_{k}^{2}+\frac{\sigma^{2}}{n}+2\sqrt{\frac{2}{\pi}}\frac{\sigma}{n\sqrt{n}}|\left(\theta
u_{i}\right)|\sum_{k}|v_{k}|.$ (47)

Since $\|\operatorname{\mathbf{v}}\|_{2}=1$, the variance of $T_{i}$ has limit
$0$. Because we cannot solve the inequality $\tau_{n,p}<\mathbb{E}T_{i}$
analytically, we leave the bound in the form given previously. ∎

In the proof above, note that if $(\tau_{n,p}-\mathbb{E}T_{i})$ is positive
and not tending to zero as $n$ grows, the quantity in (46) tends to positive
infinity when the variance decays to zero. Hence, modifying the proof by
solving $\tau_{n,p}>\mathbb{E}T_{i}$ for $|\theta u_{i}|$ yields when a
coordinate is not detectable with probability approaching $1$: i.e., when
$|\theta u_{i}|$ is smaller than the values given in (25). ∎

###  A.1 Proof of Corollary 1

The first three terms of (26) are characterized by Theorem 1(b) and by noting
that $I_{0}=I$ in the corollary. The last term can be characterized as
follows. In particular, for an algorithm with row test statistics $T_{i}$ and
a threshold $\tau$ from Table 1, we may write

$\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)=\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(T_{i}\geq\tau\right).$ (48)

For the $\ell_{2}$\- and $\ell_{1}$-SEPCA algorithms, we may reindex the
$T_{i}$ according to their order statistics $T_{(i)}$, where

$\left|T_{(1)}\right|\geq\left|T_{(2)}\right|\geq\cdots\geq\left|T_{(p-s)}\right|,$

and write

$\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)=\sum_{i=1}^{p-s}\operatorname{\mathbb{P}}\left(T_{(i)}\geq\tau\right).$
(49)

Note that there are $p-s$ null entries. For the $\ell_{2}$\- and
$\ell_{1}$-SEPCA algorithms, we have that (as a consequence of [23, Thm. 3])

$\operatorname{\mathbb{P}}\left(T_{(i)}\geq\tau\right)\leq\left(\frac{1}{ep}\right)^{\sqrt{i}},$

so that

$\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)=\sum_{i=1}^{p-s}\operatorname{\mathbb{P}}\left(T_{(i)}\geq\tau\right)\leq\sum_{i=1}^{p-s}\left(\frac{1}{ep}\right)^{\sqrt{i}}.$
(50)

The right-hand side of (50) has limit zero, as needed.

For the sum-SEPCA algorithm, from (10), it follows that there exists a non-
zero constant $\epsilon>0$ such that the threshold $\tau$ satisfies

$\tau=\sigma\left(\sqrt{2}+\epsilon\right)\sqrt{\frac{\log p}{n}}.$

In particular, from (10), we have that

$\epsilon=\frac{1+\sqrt{\log
p}/3}{\sqrt{2}\textrm{Erf}^{-1}\left(1-1/p\right)},$

where for $p>1$,

$\frac{1}{3\sqrt{2}}<\epsilon<2.$

Hence, for any $T_{i}$ such that $i\notin I$,

$\operatorname{\mathbb{P}}\left(T_{i}\geq\tau\right)\leq\exp\left(-\frac{\left(\sqrt{2}+\epsilon\right)^{2}}{2}\log
p\right)=p^{-\left(1+\sqrt{2}\epsilon+\epsilon^{2}/2\right)},$

where we have used a Gaussian tail bound [7, Sec. 2.3]. Then,

$\begin{split}\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(i\in\widehat{I}\right)&=\sum_{i\notin
I}\operatorname{\mathbb{P}}\left(T_{i}\geq\tau\right)\\\
&\leq\left(p-s\right)p^{-\left(1+\sqrt{2}\epsilon+\epsilon^{2}/2\right)}\\\
&\leq p^{-\left(\sqrt{2}\epsilon+\epsilon^{2}/2\right)}.\end{split}$ (51)

Since $\epsilon$ is larger than $1/3\sqrt{2}$, the right-hand side of (51) is
upper bounded by $p^{-13/36}$, which has limit zero, as desired.

##  Appendix B Proof of Theorem 2

###  B.1 Size of Detectable Coordinates

####  B.1.1 Sum: HC-SEPCA

If $\operatorname{\mathbf{v}}$ is equisigned, summing across the rows of
$\operatorname{\mathbf{X}}$ yields a normally distributed quantity with mean
$(\theta u_{i})\|\operatorname{\mathbf{v}}\|_{1}$ and variance $\sigma^{2}$.
Dividing by $\sigma$ and adopting the notation of $HC$, we have that under the
alternative hypothesis, $\mu_{i}=\sqrt{2r\log p}$, so that

$r=\left(\frac{|\theta u_{i}|\|\operatorname{\mathbf{v}}\|_{1}}{\sqrt{2\log
p}}\right)^{2}.$

Rearranging the inequality $r>\rho(\beta)$ yields

$|\theta u_{i}|>\sigma\sqrt{\rho(\beta)}\frac{\sqrt{2\log
p}}{\|\operatorname{\mathbf{v}}\|_{1}}.$ (52)

Note that sum-SEPCA can detect coordinates of size

$|\theta u_{i}|>\sigma C_{U}\frac{\sqrt{\log
p}}{\|\operatorname{\mathbf{v}}\|_{1}}.$ (53)

However, $C_{U}$ is strictly larger than $\sqrt{2}+1/(3\sqrt{2})$. Thus, using
HC yields a threshold of the same order, but with a strictly smaller scaling.

####  B.1.2 Sum of squares: HC-$\ell_{2}$-SEPCA

If we sum the squares of the entries of rows of $\operatorname{\mathbf{X}}$,
abusing notation slightly and using $\mathcal{N}(\mu,\sigma^{2})$ to indicate
a Gaussian random variable with mean $\mu$ and variance $\sigma^{2}$, the
statistic for the $i^{th}$ coordinate is of the form

$\sum_{k=1}^{n}\left(\frac{\sigma}{\sqrt{n}}\mathcal{N}\left(\frac{\theta
u_{i}}{\sigma}v_{k}\sqrt{n},1\right)\right)^{2}.$

Assuming oracular knowledge of $\sigma$, the statistic

$\frac{n}{\sigma^{2}}\sum_{k}X_{ik}^{2}$

places us in the setting of (15). The non-centrality parameter $\delta$ is
given by

$\delta=\sqrt{\sum_{k=1}^{n}\left(\frac{\theta
u_{i}}{\sigma}v_{k}\sqrt{n}\right)^{2}}=\left|\frac{\theta
u_{i}}{\sigma}\right|\sqrt{n}.$

Setting $\delta=2r\log p$ and solving $r>\rho(\beta)$ yields

$|\theta u_{i}|>\sigma{\rho(\beta)}{\frac{2\log p}{\sqrt{n}}}.$ (54)

We have that $\ell_{2}$-SEPCA can detect coordinates with

$|\theta u_{i}|>\sigma\sqrt{e\sqrt{2}}\sqrt{\frac{1+\log p}{\sqrt{n}}}.$ (55)

Using HC offers a significant improvement over $\ell_{2}$-SEPCA. However, we
also expect HC with the $\chi_{n}^{2}$ statistic to have a smaller detectable
coordinate: $\|\operatorname{\mathbf{v}}\|_{1}\leq\sqrt{n}$, so that for fixed
$\beta$ and $p$, the threshold in (52) is asymptotically larger than that in
(54) (but potentially of the same order). This result is strange in context of
the non-FDR results. In any case, HC improves on $\ell_{2}$-SEPCA.

####  B.1.3 FDR-SEPCA

Recall that taking sums across the rows of $\operatorname{\mathbf{X}}$, we
obtain a vector $\operatorname{\mathbf{y}}$ where $y_{i}=\mu_{i}+\sigma
z_{i}$, with $\mu_{i}=(\theta u_{i})\|\operatorname{\mathbf{v}}\|_{1}$.
Moreover, we have noted that

$t_{k}\approx\sqrt{\zeta}(1+\sqrt{2\log(\nu p/k)}),$

where $t_{k}$ is the level at which $\operatorname{\mathbf{y}}$ is
thresholded. It follows that, entries of $\operatorname{\mathbf{y}}$ that are
of size at least

$y_{i}>(1-o(1))\sqrt{\zeta}\sigma\left(1+\sqrt{2\log(\nu
p/\widehat{k})}\right)$

are selected, or, since $\mu_{i}=(\theta
u_{i})\|\operatorname{\mathbf{v}}\|_{1}$ (when $v$ is equisigned), if we
select $\widehat{k}$ coordinates, we expect to detect

$\begin{split}|\theta
u_{i}|&>(1-o(1))\sqrt{\zeta}\sigma\frac{\left(1+\sqrt{2\log(\nu
p/\widehat{k})}\right)}{\|v\|_{1}}\\\
&=O\left(\sigma\frac{\sqrt{2\log\left(\nu
p/\widehat{k}\right)}}{\|\operatorname{\mathbf{v}}\|_{1}}\right).\end{split}$
(56)

Relative to HC and sum-SEPCA, the gain here is found when there are many
smaller coordinates of $\operatorname{\mathbf{u}}$ and $\widehat{k}$ is large.

###  B.2 Proofs for the Higher Criticism-Based Methods

  1. a.

From (2.8) in [17],

$\mathbb{P}\left(T_{i}\geq\tau\right)\leq\mathbb{P}\left(\max_{j\in
I^{c}}T_{j}\geq\tau\right)$

has limit zero. ∎

  2. b.

Let $I_{1}\subseteq I$ be the set of coordinates with signal larger than the
detection limit ($i\in I$ such that $|\theta
u_{i}|>\beta_{crit}(1+\epsilon)$), and let $I_{2}\subseteq I$ contain the rest
of the coordinates ($i\in I$ such that $|\theta
u_{i}|<\beta_{crit}(1-\epsilon)$). By Theorem 1 in [2], the asymptotic power
for detecting signals below the detection limit is one, and that for signals
below the limit is zero. Hence, for $i\in I_{1}$,

$\min_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|>\beta_{crit}(1+\epsilon)\end{subarray}}\mathbb{P}\left(\text{$i$
selected}\right)\rightarrow 1,$

and for $i\in I_{2}$,

$\max_{\begin{subarray}{c}i\in I~{}:~{}|\theta
u_{i}|<\beta_{crit}(1-\epsilon)\end{subarray}}\mathbb{P}\left(\text{$i$
selected}\right)\rightarrow 0.$

As with Theorem 1, we omit the computation of $\beta_{crit}$, as it follows
from the discussion in Section 4.1. ∎

###  B.3 FDR-SEPCA

The details of these computations are in Appendix D.1, so we will summarize
the properties here.

  1. a.

The choice of $\nu=2^{1/\omega}$ controls the FDR at level $\omega$ [22].
Choosing $\omega=\omega(p)\rightarrow 0$ as $p\rightarrow\infty$ leads to an
asymptotic FDR of zero. I.e., for $i\in I^{c}$,

$\max_{i\in I^{c}}\qquad\mathbb{P}\left(\text{$i$ selected}\right)\rightarrow
0.\qed$

  2. b.

Noting that the consistency of estimating the mean vector
$\operatorname{\bm{\mu}}=\left(\theta\|v\|_{1}\right)\operatorname{\mathbf{u}}$
encompasses the estimation of the support of $\operatorname{\mathbf{u}}$, risk
bounds for the estimation of $\operatorname{\bm{\mu}}$ yield the result. To be
precise, if the expected risk
$\mathbb{E}\|\operatorname{\bm{\mu}}-\widehat{\operatorname{\bm{\mu}}}\|_{2}^{2}\leq
B$ for some bound $B$, we expect to detect coordinates of size larger than $B$
and to not detect those smaller than $B$. ∎

##  Appendix C Risk bounds under $\ell_{q}$ sparsity

In this section, we simultaneously generalize our setting to approximate
sparsity, and specify the risk lower-bounds. We omit the $\ell_{1}$-SEPCA
algorithm from consideration.

Let $\operatorname{\mathbf{u}}\in\mathbb{R}^{p}$ have unit $\ell_{2}$-norm and
belong to an $\ell_{q}$ ball with radius $C$ for $q\in(0,2]$. I.e.,

$\sum_{i=1}^{p}|u_{i}|^{q}\leq C^{q}.$ (57)

When $q=0$, we replace $C^{q}$ with $s$, the level of ‘hard’ sparsity. We have
the following result:

######  Theorem 3.

Let

$L(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}})=\left\|\operatorname{\mathbf{u}}-\text{sign}(\langle\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}}\rangle)\widehat{\operatorname{\mathbf{u}}}\right\|_{2}^{2}$
(58)

be the risk of the estimator $\widehat{\operatorname{\mathbf{u}}}$ of
$\operatorname{\mathbf{u}}$, where $\operatorname{\mathbf{u}}$ is as specified
in (1) and the estimators are the six algorithms that we have previously
described. Then,

  1. a.

sum-, HC-sum, and FDR-SEPCA have expected risks lower-bounded by

$\mathbb{E}L(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}})\geq
O\left([C^{q}-1]\|\operatorname{\mathbf{v}}\|_{1}^{-(2-q)}\right).$ (59)

  2. b.

$\ell_{2}$-SEPCA has a risk lower-bounded by

$\mathbb{E}L(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}})\geq
O\left([C^{q}-1]n^{-\frac{1}{2}(1-q/2)}\right).$ (60)

  3. c.

HC-$\ell_{2}$-SEPCA has a risk lower-bounded by

$\mathbb{E}L(\widehat{\operatorname{\mathbf{u}}},\operatorname{\mathbf{u}})\geq
O\left([C^{q}-1]n^{-(1-q/2)}\right).$ (61)

The rest of this section contains the proof of Theorem 3.

###  C.1 Proof of Theorem 3

We construct a ‘worst-case’ sparse $\operatorname{\mathbf{u}}$. Note that
$C^{q}\geq 1$ necessarily, and that if $C\geq p^{1-q/2}$, every unit norm
vector is in the $\ell_{q}$ ball. Hence, we take $C\in[1,p^{1-q/2})$.

Let $\theta$ and $\sigma$ be fixed. We want a sparse vector with several
coordinates guaranteed to be missed (the probability of not detecting them is
asymptotically $1$). For this vector $\operatorname{\mathbf{u}}$ to be sparse
and for the loss to not be $1$, set $u_{1}$ to be $\sqrt{1-r_{n}^{2}}$, where
$r_{n}^{2}=o(1)$, and take $u_{2},\cdots,u_{m_{n}+1}$ to be
${r_{n}}/{\sqrt{m_{n}}}$. The other coordinates of $u$ are $0$, so that
$\operatorname{\mathbf{u}}$ has unit $\ell_{2}$-norm.

We assume that $u_{1}$ is detected with probability $1$ as
$n\rightarrow\infty$, and want to set $u_{2},\cdots,u_{m_{n}+1}$ so that the
expected loss is lower bounded by:

$\begin{split}\mathbb{E}L(\operatorname{\mathbf{u}},\widehat{\operatorname{\mathbf{u}}})&\geq\sum_{k=1}^{p}|u_{k}|^{2}\mathbb{P}(\text{Not
Selecting Coordinate k})\\\
&\geq\sum_{k=2}^{m_{n}+1}|u_{k}|^{2}\mathbb{P}(\text{Not Selecting Coordinate
k}).\end{split}$ (62)

If coordinates of size $\frac{r_{n}}{\sqrt{m_{n}}}$ are not detected, the
expected loss is lower bounded by $r_{n}^{2}$.

Let $m_{n}=\lfloor m\rfloor$ where

$m=\delta n^{\phi}r_{n}^{\psi}\|v\|^{\eta}.$

Note that we have not specified the norm used in $\|v\|$: we will choose the
norm at the very end of the calculation. Let

$r_{n}=[C^{q}-1]^{\alpha}n^{\beta+\gamma q}\|v\|^{\kappa},$

so that,

$\frac{r_{n}}{\sqrt{m_{n}}}\approx\frac{r_{n}}{\sqrt{m}}=\frac{1}{\delta}n^{-\phi/2}\|v\|^{\kappa-\eta/2}r_{n}^{1-\psi/2}.$

We will choose $\delta,\phi,\eta,\alpha,\beta,\gamma,\kappa,\psi$ so that the
$\ell_{q}$ sparsity constraint is met and the lower bound $r_{n}^{2}$ is
maximized. The sparsity constraint requires:

$\displaystyle\begin{split}\sum_{i=1}^{p}|u_{i}|^{q}=(1-r_{n})^{q/2}+m_{n}^{1-q/2}r_{n}^{q}\leq
1+m^{1-q/2}r_{n}^{q}\leq C^{q}.\end{split}$ (63)

First, we will assume (for now) that $r_{n}=o(1)$ and that via other
parameters we may control the scaling of the coordinate sizes; hence, we set
$\psi=2$. Then,

$\begin{split}r_{n}^{q}m^{1-q/2}=\delta^{1-q/2}n^{2(\beta+\gamma
q)+(1-q/2)\phi}[C^{q}-1]^{2\alpha}\|v\|^{2\kappa+\eta(1-q/2)}.\end{split}$
(64)

We need this quantity to be smaller than $C^{q}-1$. To eliminate the $n$
dependence, we set $\beta=\frac{-\phi}{2}$ and $\gamma=\frac{\phi}{4}$. We
choose $\alpha=\frac{1}{2}$ to match powers of $[C^{q}-1]$ on both sides of
the inequality. Defining another parameter $\rho$, let
$\delta=\rho\|\operatorname{\mathbf{v}}\|^{-\eta}$. Then, the inequality is

$\rho^{1-q/2}\|\operatorname{\mathbf{v}}\|^{2\kappa}[C^{q}-1]\leq[C^{q}-1].$

Choosing $\rho\leq\|\operatorname{\mathbf{v}}\|^{-2\kappa/(1-q/2)}$ is enough.

With these choices of parameters,

$r_{n}=\sqrt{[C^{q}-1]}n^{-\frac{1}{2}\phi(1-q/2)}\|\operatorname{\mathbf{v}}\|^{\kappa},$

and

$m=\rho n^{\phi}r_{n}^{2},$

so that

$\frac{r_{n}}{\sqrt{m}}=\frac{1}{\sqrt{\rho}}n^{-\phi/2}.$

Noting that

$\frac{1}{\sqrt{\rho}}\geq\|\operatorname{\mathbf{v}}\|^{\kappa/(1-q/2)},$

choosing $\rho=\|\operatorname{\mathbf{v}}\|^{-2\kappa/(1-q/2)}$ leads to the
smallest possible choice of coordinate.

In summary:

$r_{n}=\sqrt{[C^{q}-1]}n^{-\frac{1}{2}\phi(1-q/2)}\|\operatorname{\mathbf{v}}\|^{\kappa},$
(65)
$r_{n}^{2}=[{C^{q}-1}]n^{-\phi(1-q/2)}\|\operatorname{\mathbf{v}}\|^{2\kappa},$
(66) $m=\|\operatorname{\mathbf{v}}\|^{-2\kappa/(1-q/2)}n^{\phi}r_{n}^{2},$
(67)

and

$\frac{r_{n}}{\sqrt{m}}=\|\operatorname{\mathbf{v}}\|^{\kappa/(1-q/2)}n^{-\phi/2}.$
(68)

So, for a given algorithm, it remains to choose $\phi$ and $\kappa$ so that
the worst-case risk is lower-bounded by $r_{n}^{2}$. Sum-SEPCA misses
coordinates of size $O\left(\frac{\sqrt{\log
p}}{\|\operatorname{\mathbf{v}}\|_{1}}\right)$ and $\ell_{2}$-SEPCA misses
coordinates of size $O\left(\frac{\sqrt{\log
p}}{n^{1/4}\|\operatorname{\mathbf{v}}\|_{2}}\right)$. For sum-SEPCA,
$\kappa=\frac{q-2}{2}$, and for $\ell_{2}$-SEPCA, $\kappa$ is irrelevant, as
$\|\operatorname{\mathbf{v}}\|_{2}=1$. Sum-SEPCA uses $\phi=0$ and
$\ell_{2}$-SEPCA uses $\phi=\frac{1}{2}$. Hence, sum-SEPCA has a risk lower-
bounded by

$O\left([C^{q}-1]\|\operatorname{\mathbf{v}}\|_{1}^{-(2-q)}\right).$ (69)

Noting that $\|\operatorname{\mathbf{v}}\|_{2}=1$, $\ell_{2}$-SEPCA has

$\begin{split}&O\left([C^{q}-1]n^{-\frac{1}{2}(1-q/2)}\|\operatorname{\mathbf{v}}\|_{2}^{-(1-q/2)}\right)=\\\
&O\left([C^{q}-1]n^{-\frac{1}{2}(1-q/2)}\right).\end{split}$ (70)

In the $\ell_{0}$ case, i.e., when $\operatorname{\mathbf{u}}$ has no more
than $s$ non-zero entries, the preceding analysis goes through with $C^{q}$
replaced by $s$ and $q$ set to zero.

####  C.1.1 FDR Algorithms

For HC-sum-SEPCA, the $\beta_{crit}$ is of the same order as that for sum-
SEPCA. Similarly, for FDR-SEPCA, if $\widehat{k}$ is much smaller than $p$,
$\beta_{crit}$ is of roughly the same order. Hence, these two algorithms have
the same risk bound as sum-SEPCA. For HC-$\ell_{2}$-SEPCA, $\kappa=0$ and
$\phi=1$. The risk is therefore lower-bounded by

$O\left([C^{q}-1]n^{-(1-q/2)}\right).$ (71)

##  Appendix D FDR-SEPCA: Further Details

Let $y_{i}=\mu_{i}+\sigma z_{i}$, where $i\in\\{1,\cdots,p\\}$, the vector
$\operatorname{\mathbf{z}}$ of the $z_{i}$ is normally distributed with mean
$0$ and covariance $\operatorname{\bm{\Sigma}}$, and
$\operatorname{\bm{\Sigma}}$ satisfies

$\xi_{o}\mathcal{I}_{p}\leq\operatorname{\bm{\Sigma}}\leq\xi_{1}\mathcal{I}_{p}.$

Here, $\xi_{0}$ is the smallest eigenvalue of $\operatorname{\bm{\Sigma}}$ and
$\xi_{1}$ is the largest. The mean vector $\operatorname{\bm{\mu}}$ of the
$\mu_{i}$ is assumed to be sparse; the goal is to estimate
$\operatorname{\bm{\mu}}$. The following penalized least squares formulation
yields an estimator for $\operatorname{\bm{\mu}}$:

$\widehat{\operatorname{\bm{\mu}}}=\arg\min_{\operatorname{\bm{\mu}}}\|\operatorname{\mathbf{y}}-\operatorname{\bm{\mu}}\|_{2}^{2}+\sigma^{2}\text{pen}\left(\|\operatorname{\bm{\mu}}\|_{0}\right),$
(72)

where $\text{pen}(k)$ is defined as

$\text{pen}(k)=\xi_{1}\zeta k\left(1+\sqrt{2L_{p,k}}\right)^{2},$ (73)

with $\zeta>1$ and

$L_{p,k}=(1+2\beta)\log(\nu p/k).$ (74)

The parameter $\beta$ may be set to $0$ here, and $\nu$ is chosen to be no
smaller than $e^{1/(1+2\beta)}$. We define $\|\operatorname{\bm{\mu}}\|_{0}$
to be number of non-zero coordinates of $\operatorname{\bm{\mu}}$.

The solution to (72) is given by hard-thresholding. Let $|y|_{(i)}$ be the
$i^{th}$ order statistic of $|y_{i}|$, namely
$|y|_{(1)}\geq\cdots\geq|y|_{(p)}$. Then if

$\widehat{k}=\arg\min_{k\geq
0}\sum_{i>k}|y|_{(i)}^{2}+\sigma^{2}\text{pen}(k),$ (75)

defining

$t_{k}^{2}=\text{pen}(k)-\text{pen}(k-1),$ (76)

the solution is to hard threshold at $t_{\widehat{k}}$.

In this set-up, we have that

$t_{k}\approx\lambda_{p,k}=\sqrt{\xi_{1}\zeta}(1+\sqrt{2L_{p,k}}),$

with $|t_{k}-\lambda_{p,k}|<c/\lambda_{k}$. More precisely, Lemma 11.7 of [20]
says that

$\lambda_{p,k}-\frac{4\zeta b}{\lambda_{p,k}}\leq t_{k}\leq\lambda_{p,k}.$

When $\nu\geq e^{2}$, we may take $b=(1+2\beta)$. In any case, if $k=o(n)$,
$\lambda_{p,k}\asymp\sqrt{\log p}$. Hence, entries of
$\operatorname{\mathbf{y}}$ that are of size at least

$y_{i}>(1-o(1))\sqrt{\xi_{1}\zeta}\sigma\left(1+\sqrt{2\log(\nu
p/\widehat{k})}\right)$

are selected, or, since $\mu_{i}=(\theta
u_{i})\|\operatorname{\mathbf{v}}\|_{1}$ (when $\operatorname{\mathbf{v}}$ is
equisigned), if we select $\widehat{k}$ coordinates, we expect to detect

$\begin{split}|\theta
u_{i}|&>(1-o(1))\sqrt{\xi_{1}\zeta}\sigma\frac{\left(1+\sqrt{2\log(\nu
p/\widehat{k})}\right)}{\|\operatorname{\mathbf{v}}\|_{1}}\\\
&=O\left(\sigma\sqrt{\xi_{1}}\frac{\sqrt{2\log\left(\nu
p/\widehat{k}\right)}}{\|\operatorname{\mathbf{v}}\|_{1}}\right).\end{split}$
(77)

###  D.1 Risk Behavior

Recalling that (72) solves a penalized least squares problem for
$\widehat{\operatorname{\mathbf{y}}}$ close to $\operatorname{\mathbf{y}}$, we
may discuss the statistical behavior of this estimator. The following
discussion follows and reproduces that in [22]

First, note that for $\beta=0$, the parameter $\nu$ directly controls the FDR
(where a false positive corresponds to selecting a zero coordinate in
$\operatorname{\mathbf{y}}$): a choice of $\nu=2^{1/\omega}$ for
$\omega\in(0,1)$ bounds the FDR at a level $\omega$.

Second, the expected risk,
$\mathbb{E}\|\operatorname{\mathbf{y}}-\widehat{\operatorname{\mathbf{y}}}\|_{2}^{2}$,
is bounded as follows. By Proposition 4.1 in [22],

$\mathbb{E}\|\operatorname{\mathbf{y}}-\widehat{\operatorname{\mathbf{y}}}\|_{2}^{2}\leq
D\left[2M_{p}^{\prime}\xi_{1}\sigma^{2}+\mathcal{R}(\operatorname{\mathbf{y}},\sigma)\right],$
(78)

where $D$ is a constant $2\zeta(\zeta+1)^{3}(\zeta-1)^{-3}=\Theta(1)$, we
assume that $\xi_{1}=1$, and $0\leq M_{p}^{\prime}\leq
C_{\beta}p^{-2\beta}\nu^{-1}$, for some $C_{\beta}>0$. Since $\beta=0$,
$M_{p}^{\prime}=O(1/\nu)=O(\omega)$, if we control the FDR at level $\omega$.

The second term in (78) is the ideal risk, or, the infimum of the penalized
least squares objective. If $\operatorname{\mathbf{y}}$ belongs to an
$\ell_{q}$ ball with radius $C$ and $0<q<2$, and we define

$r_{p,q}(C)=\left\\{\begin{array}[]{ll}C^{2}&\text{ if }C\leq\sqrt{1+\log
p},\\\ C^{q}[1+\log(p/C^{q})]^{1-q/2}&\text{ if }\sqrt{1+\log p}\leq C\leq
p^{1/q},\\\ p&\text{ if }C\geq p^{1/q},\end{array}\right.$ (79)

the ideal risk is bounded as

$\sup_{\operatorname{\mathbf{y}}\in\mathbb{R}^{p}:\sum_{i}|y_{i}|^{q}\leq
C^{q}}\mathcal{R}(\operatorname{\mathbf{y}},\sigma)\leq
c(\log\nu)\sigma^{2}r_{p,q}(C/\sigma),$ (80)

for some $c>0$. The supplementary results in [22] yield that
$\mathcal{R}(\operatorname{\mathbf{y}},\sigma)$ is bounded by $C^{2}\log\nu$,
and by $C^{2}$ when $C\leq\sqrt{1+\log p}$.

As in Appendix C, we may replace $q$ with $0$ and $C^{q}$ with $s$ in the case
of hard sparsity with $s$ non-zero coordinates. Doing so leads to the bound:

$\mathbb{E}\|\operatorname{\mathbf{y}}-\widehat{\operatorname{\mathbf{y}}}\|_{2}^{2}\leq
s\sigma\log\nu\log\frac{\sigma p\nu}{s}+\sigma_{1}\sigma^{2}\frac{2}{\nu}.$
(81)

Note that we have recovered the factor of $\log\nu p/s$ in $\beta_{crit}$.

## References

  * [1] Zeyuan Allen-Zhu and Yuanzhi Li.  LazySVD: even faster SVD decomposition yet without agonizing pain.  In Advances in Neural Information Processing Systems, pages 974–982, 2016.
  * [2] Ery Arias-Castro, Emmanuel J Candès, and Yaniv Plan.  Global testing under sparse alternatives: ANOVA, multiple comparisons and the higher criticism.  Annals of Statistics, pages 2533–2556, 2011.
  * [3] Florent Benaych-Georges and Raj Rao Nadakuditi.  The singular values and vectors of low rank perturbations of large rectangular random matrices.  Journal of Multivariate Analysis, 111:120–135, October 2012.
  * [4] Quentin Berthet, Philippe Rigollet, et al.  Optimal detection of sparse principal components in high dimension.  Annals of Statistics, 41(4):1780–1815, 2013.
  * [5] Aharon Birnbaum, Iain M Johnstone, Boaz Nadler, and Debashis Paul.  Minimax bounds for sparse PCA with noisy high-dimensional data.  Annals of Statistics, 41(3):1055, 2013.
  * [6] Sergey G Bobkov and Fedor L Nazarov.  On convex bodies and log-concave probability measures with unconditional basis.  In Geometric aspects of functional analysis, pages 53–69. Springer, 2003.
  * [7] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration inequalities: A nonasymptotic theory of independence.  Oxford University Press, 2013.
  * [8] Stéphane Boucheron and Maud Thomas.  Concentration inequalities for order statistics.  Electronic Communications in Probability, 17:no. 51, 1–12, 2012\.
  * [9] Matthew Brennan, Guy Bresler, and Wasim Huleihel.  Reducibility and computational lower bounds for problems with planted sparse structure.  arXiv preprint arXiv:1806.07508, 2018.
  * [10] Cristina Butucea, Mohamed Ndaoud, Natalia A Stepanova, Alexandre B Tsybakov, et al.  Variable selection with Hamming loss.  The Annals of Statistics, 46(5):1837–1875, 2018.
  * [11] Alexandre d’Aspremont, Laurent El Ghaoui, Michael I Jordan, and Gert RG Lanckriet.  A direct formulation for sparse PCA using semidefinite programming.  SIAM Review, 49(3):434–448, 2007.
  * [12] Chris HQ Ding, Tao Li, and Michael I Jordan.  Convex and semi-nonnegative matrix factorizations.  IEEE Transactions on Pattern analysis and Machine Intelligence, 32(1):45–55, 2010.
  * [13] David Donoho and Jiashun Jin.  Higher criticism for detecting sparse heterogeneous mixtures.  Annals of Statistics, pages 962–994, 2004.
  * [14] David Donoho, Jiashun Jin, et al.  Higher criticism for large-scale inference, especially for rare and weak effects.  Statistical Science, 30(1):1–25, 2015.
  * [15] François Le Gall and Florent Urrutia.  Improved rectangular matrix multiplication using powers of the Coppersmith-Winograd tensor.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1029–1046. SIAM, 2018.
  * [16] Chen Gao, Brian E Moore, and Raj Rao Nadakuditi.  Augmented robust PCA for foreground-background separation on noisy, moving camera video.  In 2017 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pages 1240–1244. IEEE, 2017.
  * [17] Peter Hall, Jiashun Jin, et al.  Innovated higher criticism for detecting sparse signals in correlated noise.  Annals of Statistics, 38(3):1686–1732, 2010.
  * [18] Patrik O Hoyer.  Non-negative matrix factorization with sparseness constraints.  Journal of Machine Learning Research, 5(Nov):1457–1469, 2004.
  * [19] Kejun Huang, Nicholas D Sidiropoulos, and Ananthram Swami.  Non-negative matrix factorization revisited: Uniqueness and algorithm for symmetric decomposition.  IEEE Transactions on Signal Processing, 62(1):211–224, 2014.
  * [20] Iain M Johnstone.  Gaussian estimation: Sequence and wavelet models.  Unpublished manuscript, 2017.  http://statweb.stanford.edu/~imj/GE_08_09_17.pdf.
  * [21] Iain M Johnstone and Arthur Yu Lu.  On consistency and sparsity for principal components analysis in high dimensions.  Journal of the American Statistical Association, 104(486):682, 2009\.
  * [22] Iain M Johnstone and Debashis Paul.  Adaptation in some linear inverse problems.  Stat, 3(1):187–199, 2014.
  * [23] Rafał Latała.  Order statistics and concentration of norms for log-concave vectors.  Journal of Functional Analysis, 261(3):681 – 696, 2011.
  * [24] Haifeng Liu, Zhaohui Wu, Xuelong Li, Deng Cai, and Thomas S Huang.  Constrained nonnegative matrix factorization for image representation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(7):1299–1311, 2012.
  * [25] Zongming Ma et al.  Sparse principal component analysis and iterative thresholding.  Annals of Statistics, 41(2):772–801, 2013.
  * [26] Andrea Montanari and Emile Richard.  Non-negative principal component analysis: Message passing algorithms and sharp asymptotics.  IEEE Transactions on Information Theory, 62(3):1458–1484, 2015\.
  * [27] Brian Moore, Chen Gao, and Raj Rao Nadakuditi.  Panoramic robust PCA for foreground-background separation on noisy, free-motion camera video.  IEEE Transactions on Computational Imaging, 2019.
  * [28] Mohamed Ndaoud.  Interplay of minimax estimation and minimax support recovery under sparsity.  arXiv preprint arXiv:1810.05478, 2018.
  * [29] Damien Passemier, Zhaoyuan Li, and Jianfeng Yao.  On estimation of the noise variance in high dimensional probabilistic principal component analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(1):51–67, 2017.
  * [30] Dominique Pastor and Francois-Xavier Socheleau.  Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences.  IEEE transactions on Signal Processing, 60(4):1545–1555, 2012.
  * [31] Amelia Perry, Alexander S Wein, Afonso S Bandeira, Ankur Moitra, et al.  Optimality and sub-optimality of PCA I: Spiked random matrix models.  The Annals of Statistics, 46(5):2416–2451, 2018.
  * [32] Massimo Piccardi.  Background subtraction techniques: A review.  In 2004 IEEE International Conference on Systems, Man and Cybernetics (IEEE Cat. No. 04CH37583), volume 4, pages 3099–3104. IEEE, 2004\.
  * [33] Pradeep Ravikumar, Martin J Wainwright, and et al. Lafferty, John D.  High-dimensional Ising model selection using $\ell_{1}$-regularized logistic regression.  Annals of Statistics, 38(3):1287–1319, 2010.
  * [34] Galen Reeves and Michael Gastpar.  Sampling bounds for sparse support recovery in the presence of noise.  In 2008 IEEE International Symposium on Information Theory, pages 2187–2191. IEEE, 2008.
  * [35] Bin Ren, Laurent Pueyo, Guangtun Ben Zhu, John Debes, and Gaspard Duchêne.  Non-negative matrix factorization: robust extraction of extended structures.  The Astrophysical Journal, 852(2):104, 2018.
  * [36] Paul Ross.  Stars.  https://archive.org/details/Stars_2D.  Accessed: 2016 November 16.
  * [37] Francois-Xavier Socheleau and Dominique Pastor.  Testing the energy of random signals in a known subspace: An optimal invariant approach.  IEEE Signal Processing Letters, 21(10):1182–1186, 2014.
  * [38] Leo Taslaman and Björn Nilsson.  A framework for regularized non-negative matrix factorization, with application to the analysis of gene expression data.  PloS one, 7(11):e46331, 2012.
  * [39] Namrata Vaswani, Thierry Bouwmans, Sajid Javed, and Praneeth Narayanamurthy.  Robust subspace learning: Robust PCA, robust subspace tracking, and robust subspace recovery.  IEEE signal processing magazine, 35(4):32–55, 2018.
  * [40] Yu-Xiong Wang and Yu-Jin Zhang.  Nonnegative matrix factorization: A comprehensive review.  IEEE Transactions on Knowledge and Data Engineering, 25(6):1336–1353, 2013.
  * [41] Xiao-Tong Yuan and Tong Zhang.  Truncated power method for sparse eigenvalue problems.  Journal of Machine Learning Research, 14(Apr):899–925, 2013.
  * [42] Ron Zass and Amnon Shashua.  Nonnegative sparse PCA.  In Advances in neural information processing systems, pages 1561–1568, 2007.