Sparse Equisigned PCA: Algorithms and Performance Bounds in the Noisy Rank-1 Setting

TL;DR

This paper introduces algorithms for sparse equisigned PCA that outperform traditional SVD in noisy, high-dimensional settings, providing theoretical bounds and practical validation for detecting sparse principal components.

Contribution

It proposes new algorithms exploiting sparsity and equisignedness, with theoretical performance bounds and empirical validation in noisy, high-dimensional regimes.

Findings

Sparse PCA algorithms recover signals below SVD's eigen-SNR threshold.

Coordinate selection based on sum-type statistic outperforms norm-based schemes.

Lower bounds on detectable coordinate size and worst-case risk are established.

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 \textit{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 illustrate the performance with a video data example, where the interest is in identifying objects.

Tables

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

Equations

$$X=\theta u{v}^T+\sigma G.$$

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

$$\mathbb{E}(X{X}^T)={\theta }^{2}u{u}^T+{\sigma }^2{I}_p,$$

$$\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.$$

$$\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.$$

$$d_H\left(I,\hat{I}\right)=\left|\left(I\cup \hat{I}\right)\backslash \left(I\cap \hat{I}\right)\right|,$$

$${\theta }^2{u}_i^2+{\sigma }^2\sum _{j=1}^n{G}_{ij}^2,$$

$$\mathbb{P}\left(\underset{1\le i\le p}{\max }{T}_i\ge {\tau }_{n,p}\right)\le \frac{1}{ep}\to 0,$$

$$C_2=\sqrt{2}K\text{ and }{C}_1=K\sqrt{\left(1-2/\pi \right)}.$$

$${\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).$$

$$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},$$

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

$${\kappa }_U\ge \frac{\sqrt{2}}{U(p)}\left(3+\sqrt{\log p}\right)>1\text{ and }{C}_U=\sqrt{2}+\frac{{\kappa }_U}{3\sqrt{2}}$$

$$\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}$$

$${\mu }_i={\mu }_p=\sqrt{2r\log p},$$

$$H{C}_p=\underset{i:1/p\le p_{(i)}\le 1/2}{\max }\frac{\sqrt{p}\left(i/p-p_{(i)}\right)}{\sqrt{p_{(i)}(1-p_{(i)})}}.$$

$$\frac{H{C}_p}{\sqrt{2\log \log p}}\to 1$$

$$\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.$$

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

$$\hat{\mu }=\arg \underset{\mu }{\min }\Vert y-\mu {\Vert }_2^2+{\sigma }^2\text{pen}\left(\Vert \mu {\Vert }_0\right),$$

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

$$\hat{k}=\arg \underset{k\ge 0}{\min }\sum _{i>k}|y{|}_{(i)}^2+{\sigma }^2\text{pen}(k),$$

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

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

$$L(\hat{u},u)={\Vert u-\text{sign}(⟨u,\hat{u}⟩)\hat{u}\Vert }_2^2.$$

$$\underset{u\in {\mathbb{S}}^{p-1}:\Vert u{\Vert }_0\le s}{\sup }\mathbb{E}L\left(\hat{u},u\right),$$

$$\Vert u-\hat{u}{\Vert }_2^2=\underset{\begin{array}{c}\text{Estimation Error from }\\ \text{detected coordinates}\end{array}}{\underbrace{\Vert u_{\hat{I}}-\hat{u}{\Vert }_2^2}}+\underset{\begin{array}{c}\text{Error from }\\ \text{missed coordinates}\end{array}}{\underbrace{\Vert u_{{\hat{I}}^c}{\Vert }_2^2}}\ge \Vert u_{{\hat{I}}^c}{\Vert }_2^2.$$

$$\begin{array}{rl}\mathbb{E}L(u,\hat{u})& \ge \sum _{k=1}^p|u_k{|}^2\mathbb{P}\left(\text{Not Selecting Coordinate k}\right)\\ & \ge \sum _{k=2}^{m+1}|u_k{|}^2\mathbb{P}\left(\text{Not Selecting Coordinate k}\right),\end{array}$$

$$\mathbb{E}L(u,\hat{u})\ge r^2\mathbb{P}\left(\text{Not Selecting Coordinate k}\right).$$

$$\underset{i\in I^c}{\max }\mathbb{P}\left(i\in \hat{I}\right)\to 0,$$

Full text

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

[TABLE]

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 [math] 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\}$,

[TABLE]

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

[TABLE]

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$,

[TABLE]

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$,

[TABLE]

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,

[TABLE]

we want the expected Hamming distance $\operatorname{\mathbb{E}}d_{H}\left(I,\widehat{I}\right)$ to have limit [math], 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}}$.

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

[TABLE]

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,

[TABLE]

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]

*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

[TABLE]

3.4.2 sum-SEPCA

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

[TABLE]

In (7), we have that

[TABLE]

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

[TABLE]

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

[TABLE]

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].

4.1 Higher Criticism

4.1.1 Formulation

Assume we have $p$ independent tests of the form

[TABLE]

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

[TABLE]

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:

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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$:

[TABLE]

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

[TABLE]

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

[TABLE]

defining

[TABLE]

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

In this set-up, we have that

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math]. Then, the expected loss has the lower bound

[TABLE]

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

[TABLE]

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}$,

[TABLE] 2. b.

For $i\in I$,

[TABLE]

[TABLE]

Here

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$,

[TABLE] 2. b.

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

[TABLE]

[TABLE] 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

[TABLE]

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

[TABLE]

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 [math] 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).

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.

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

[TABLE]

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).

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 [math] 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.

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.

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:

[TABLE]

It is sufficient to guarantee that

[TABLE]

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

[TABLE]

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

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]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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 [math]. 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

[TABLE]

and

[TABLE]

Since $\sigma$ and $\theta$ are fixed, the variance always decays to [math].

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

[TABLE]

which is less than or equal to

[TABLE]

Since $\|\operatorname{\mathbf{v}}\|_{2}=1$, the variance of $T_{i}$ has limit [math]. 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

[TABLE]

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

[TABLE]

and write

[TABLE]

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])

[TABLE]

so that

[TABLE]

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

[TABLE]

In particular, from (10), we have that

[TABLE]

where for $p>1$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Note that sum-SEPCA can detect coordinates of size

[TABLE]

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

[TABLE]

Assuming oracular knowledge of $\sigma$, the statistic

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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],

[TABLE]

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}$,

[TABLE]

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

[TABLE]

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}$,

[TABLE] 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.,

[TABLE]

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

Theorem 3**.**

Let

[TABLE]

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

[TABLE] 2. b.

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

[TABLE] 3. c.

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

[TABLE]

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 [math], 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:

[TABLE]

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

[TABLE]

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

[TABLE]

so that,

[TABLE]

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:

[TABLE]

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,

[TABLE]

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

[TABLE]

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

With these choices of parameters,

[TABLE]

and

[TABLE]

so that

[TABLE]

Noting that

[TABLE]

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

In summary:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] and covariance $\operatorname{\bm{\Sigma}}$, and $\operatorname{\bm{\Sigma}}$ satisfies

[TABLE]

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}}$:

[TABLE]

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

[TABLE]

with $\zeta>1$ and

[TABLE]

The parameter $\beta$ may be set to [math] 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

[TABLE]

defining

[TABLE]

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

In this set-up, we have that

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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],

[TABLE]

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

[TABLE]

the ideal risk is bounded as

[TABLE]

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 [math] and $C^{q}$ with $s$ in the case of hard sparsity with $s$ non-zero coordinates. Doing so leads to the bound:

[TABLE]

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