On Stein’s Identity and Near-Optimal Estimation in

arXiv:1709.08795v1 [] 26 Sep 2017

High-dimensional Index Models
Zhuoran Yang∗1 , Krishnakumar Balasubramanian†1 , and Han Liu‡1
1

Department of Operations Research and Financial Engineering, Princeton
University
September 27, 2017

Abstract
We consider estimating the parametric components of semi-parametric multiple index
models in a high-dimensional non-Gaussian setting. Our estimators leverage the score
function based second-order Stein’s lemma and do not require Gaussian or elliptical
symmetry assumptions made in the literature. Moreover, to handle score functions and
response variables that are heavy-tailed, our estimators are constructed via carefully
thresholding their empirical counterparts. We show that our estimator achieves nearoptimal statistical rate of convergence in several settings. We supplement our theoretical
results via simulation experiments that confirm the theory.

∗

zy6@princeton.edu
kb18@princeton.edu
‡
hanliu@princeton.edu
†

1

1

Introduction

Consider the semi-parametric index model relating the response (Y ) and the covariate (X) by
Y = f (hβ1∗ , Xi, . . . , hβk∗ , Xi) + ,

(1.1)

where X, {β`∗ }`∈[k] ∈ Rd and {β`∗ }`∈[k] are assumed to be structured and  is a zero-mean noise
that is independent of X. Here the vectors {β`∗ }`∈[k] are the parametric components and the
function f is the nonparametric component or the link function. Such a model is called
as sparse multiple index model (MIM) in the literature (Chen et al., 2010). In this work,
given n i.i.d samples {Xi , Yi }ni=1 from the above model, where n < d, we are concerned with
estimating the parametric components {β`∗ }`∈[k] when f is unknown. More importantly, we
do not impose the assumption that X is Gaussian, which is commonly made in the literature.
Two special cases of our model include phase retrieval, for which k = 1, and dimensionality
reduction, for which k ≥ 1. Motivated by these applications, we make a distinction between
the case of k = 1 (which is also called as single index model (SIM)) and k > 1 in the rest of
the paper.
Furthermore, note that when k = 1 and f is the identity function, we obtain the standard
linear regression model. Estimators for high-dimensional linear models have been developed
and analyzed extensively in the last decade (see for example (Bühlmann and van de Geer,
2011; Vershynin, 2015) for comprehensive overviews). While being a useful testbed for illustrating conceptual phenomenon, they often suffer from a lack of flexibility in modeling realworld situations. On the other hand, completely nonparametric models, although flexible,
suffer from the curse of dimensionality unless restrictive additive sparsity or smoothness assumptions are imposed (Yuan et al., 2016). An interesting compromise between the parametric and nonparametric models is provided by the semi-parametric index models (Horowitz,
2009) described in (1.1).
Estimating the parametric components {β`∗ }`∈[k] without the dependence on the exact
form of the link function appears naturally in several situations. For example, in one-bit
compressed sensing (Boufounos and Baraniuk, 2008) and sparse generalized linear models (Loh and Wainwright, 2015), we are interested in recovering the true signal vector based
on nonlinear measurements. In sufficient dimensionality reduction, where k is typically a
fixed number greater than one but much less than d, we would like to estimate the projection
2

onto the subspace spanned by the parametric components {β`∗ }`∈[k] without depending on
the specific form of the function f . Furthermore, in deep neural networks (DNN), which are
cascades of the MIM, the nonparametric component corresponds to the activation function
which is pre-specified and the task is to estimate the linear components, which are used for
prediction in the test stage. Hence, it is crucial to develop estimators for the linear component with both statistical accuracy and computational efficiency for a wide class of link
functions.

1.1

Intriguing Aspects of Index Models.

Several subtle issues arise when we consider optimal estimation in SIM and MIM. Specifically,
most existing results depend crucially on the assumption made on X or f and fail to hold
when those assumptions are relaxed. Such issues arise even in low-dimensional settings,
where n > d. Consider, for example, the case of k = 1 and a known link function f (u) = u2 .
This corresponds to phase retrieval, which is a challenging inverse problem that has regained
interest in the last few years along with the success of compressed sensing. A straightforward
way to estimate β ∗ is to do nonlinear least squares regression (Lecué and Mendelson, 2015),
which is a nonconvex optimization problem. Candès et al. (2013) propose an estimator based
on convex relaxations. Although their estimator is optimal when X is sub-Gaussian, they
are not agnostic to the link function., i.e., the same result does not hold if the link function
is changed. Direct optimization of the nonconvex phase retrieval problem was considered
by Candès et al. (2015) and Sun et al. (2016), which propose estimators based on iterative
algorithms that are statistical optimal. However, they rely on the assumption that X is
Gaussian. A careful look at their proofs reveal that extending them to a wider class of
distributions is significantly challenging – for example, they require sharp concentration
inequalities for polynomials of degree four of X, which would lead to suboptimal rate even
when X is sub-Gaussian. Furthermore, their results are not agnostic to the link function
as well. Similar observations could be made for both convex (Li and Voroninski, 2013) and
nonconvex estimators (Cai et al., 2015) for sparse phase retrieval in high dimensions.
In addition, a surprising result for SIM was established in Plan and Vershynin (2016).
They show that when X is Gaussian, even when the link function is unknown, one could
estimate β ∗ at the optimal statistical rate with the convex Lasso estimator. Unfortunately,
3

their assumptions on the link function is rather restrictive and rule out several interesting
models including phase retrieval. Furthermore, none of the above procedures are applicable
to the case of MIMs. A line of work pioneered by Ker-Chau Li, focused on estimation in
MIMs in the low-dimensional setting. We provide a discuss about them in the related work
section, but they again require restrictive assumption on either the link functions or on the
distribution of X. For example, X is required to be elliptically symmetric, which prevents
it from being widely applicable.

1.2

Motivation

Our work is primarily motivated by an interesting phenomenon illustrated in (Plan and
Vershynin, 2016) for a class of high-dimensional SIM. Below, we first briefly summarize the
result from (Plan and Vershynin, 2016) and then provide our alternative justification for the
same result via Stein’s identity. We mainly leverage this alternative justification and propose
our estimators for the more general setting we consider. Assuming, for simplicity, we work
in the one-dimensional setting and are given n i.i.d. samples from the SIM. Consider the
least-squares estimator
n

1X
(Yi − Xi β)2 .
βbLS = argmin
n
β∈R
i=1
Note that the above estimator is the standard least-squares estimator assuming a linear
model (i.e., identity link function). The surprising observation from (Plan and Vershynin,
2016) is that, under the crucial assumption that X is standard Gaussian, βbLS is a good
estimator of β ∗ (up to a scaling) even when the data is generated from a nonlinear SIM. The
same holds true for the high-dimensional setting when the minimization is performed in an
appropriately constrained norm-ball (for example, the `1 -ball). Hence the theory developed
for the linear setting could be leveraged to understand the performance in the SIM setting.
Below, we give an alternative justification for the above estimator as an implication of Stein’s
identity in the Gaussian case, which is summarized as follows.
Proposition 1.1 (Gaussian Stein’s Identity (Stein, 1972)). Let X ∼ N (0, 1) and g : R → R
be a continuos function such that E|g 0 (X)| ≤ ∞. Then we have E[g(X)X] = E[g 0 (X)].
Note that in our context, we have E[f 0 (X)] ∝ β ∗ and E[f (X)X] = E[Y ·X]. Now consider
the following estimator, which is based on performing least-squares on the sample version of
4

the above proposition:
n

1X
βbSL = argmin
(Yi Xi − β)2
n
β∈R
i=1
Note that βbLS and βbSL are the same estimators assuming X ∼ N (0, 1), as n → ∞. This
observation leads to an alternative interpretation of the estimator proposed by (Plan and
Vershynin, 2016) via Stein’s identity for Gaussian random variables. Thus it provides an
alternative justification for why the linear least-squares estimator should work in the SIM
setting. Interestingly, a similar procedure based on second-order Stein’s identity (see §2 for
precise definitions) was used in Candès et al. (2015) to provide a favorable initializer for their
gradient descent algorithm for phase retrieval. Our observation also provides an alternative
interpretation of the initialization method used in Candès et al. (2015) by appealing to
Stein’s identity. These observations also naturally leads to leveraging non-Gaussian versions
of Stein’s identity for dealing with non-Gaussian covariates. Our estimators based on this
motivation is described in detail in §3 and §4.

1.3

Related Work

There is a significant body of work on SIMs in the low-dimensional setting. They are based
on assumptions on either the distribution of the covariate or the link functions. Assuming a
monotonic link function, Han (1987); Sherman (1993) propose the maximum rank correlation
estimator exploiting the relationship between monotonic functions and rank-correlations.
Furthermore, Li and Duan (1989) propose an estimator for a wide class of unknown link
functions under the assumption that the covariate follows a symmetric elliptical distribution. This assumption is restrictive as often times the covariates are not from a symmetric
distribution. For example, in several economic applications where the covariates are usually
highly skewed and heavy-tailed (Horowitz, 2009).
The success of Lasso and related linear estimators in high-dimensions (Bühlmann and
van de Geer, 2011), also enabled the exploration of high-dimensional SIMs. Although, this
is very much work in progress. As mentioned previously, Plan and Vershynin (2016) show
that the Lasso estimator works for the SIMs in high dimensions when the data is Gaussian.
A more tighter albeit an asymptotic results under the same setting was proved in Thrampoulidis et al. (2015). Very recently Goldstein et al. (2016) extend the results of Li and Duan
5

(1989) to the high dimensional setting but it suffers from similar problems as mentioned in
the low-dimensional setting. For the case of monotone nonparametric component, Yang et al.
(2015) analyze a non-convex least squares approach under the assumption that the data is
sub-Gaussian. However, the success of their method hinges on the knowledge of the link
function. Furthermore, Jiang and Liu (2014); Lin et al. (2015); Zhu et al. (2006) analyze
the sliced inverse regression estimator in the high-dimensional setting concentrating mainly
on support recovery and consistency properties. Similar to the low-dimensional case, the
assumptions made on the covariate distribution restrict them from several real-world applications involving non-Gaussian or non-symmetric covariate, for example high-dimensional
problems in economics (Fan et al., 2011). Furthermore, several results are established on a
case-by-case basis for fixed link function. Specifically Boufounos and Baraniuk (2008); Ai
et al. (2014) and Davenport et al. (2014) consider 1-bit compressed sensing and matrix completion respectively, where the link is assumed to be the sign function. Also, Waldspurger
et al. (2015) and Cai et al. (2015) propose and analyze convex and non-convex estimators
for phase retrieval respectively, in which the link is the square function. All the above works,
except Ai et al. (2014) make Gaussian assumptions on the data and are specialized for the
specific link functions. The non-asymptotic result obtained in Ai et al. (2014) is under subGaussian assumptions, but the estimator is not consistent. Finally, there is a line of work
focussing on estimating both the parametric and the nonparametric component Kalai and
Sastry (2009); Kakade et al. (2011); Alquier and Biau (2013); Radchenko (2015). We do not
focus on this situation in this paper as mentioned before.
For multiple index models, relatively less work exist in the high-dimensional setting. In
the low-dimensional setting, a line of work for estimation in MIMs is proposed by KerChau Li, which include inverse regression (Li, 1991), principal Hessian directions (Li, 1992)
and regression under link violation (Li and Duan, 1989). The proposed estimators are
applicable for a class of unknown link functions under the assumption that the covariate
follows a Gaussian or symmetric elliptical distribution. Such an assumption is restrictive
as often times the covariates are heavy-tailed or skewed (Horowitz, 2009; Fan et al., 2011).
Furthermore, they concentrate only on the low-dimensional setting establishing asymptotic
statements. Estimation in high-dimensional MIM under the subspace sparsity assumption
was considered in Chen et al. (2010), where the results are asymptotic and the proposed

6

estimators are not computable in polynomial time.
To summarize, all the above works require restrictive assumption on either the data
distribution or on the link function. We propose and analyze an estimator for a class of
(unknown) link functions for the case when the covariates are drawn from a non-Gaussian
distribution – under the assumption that we know the distribution a priori. Note that in
several situations, one could fit specialized distributions, to real-world data that is often
times skewed and heavy-tailed, so that it provides a good generative model of the data.
Also, mixture of Gaussian distribution, with the number of components selected appropriately, approximates the set of all square integrable distributions to arbitrary accuracy (see
for example McLachlan and Peel (2004)). Furthermore, since this is a density estimation
problem it is unlabeled and there is no issue of label scarcity. Hence it is possible to get
accurate estimate of the distribution in most situations of interest. Thus our work is complementary to the existing literature and provides an estimator for a class of models that is
not addressed in the previous works.

1.4

Contributions

As we saw before, there are several subtleties based on the interplay between the assumptions made on X and f when dealing with estimation in SIM and MIM. Thus an interesting
question is, whether it is possible to estimate the linear components in SIMs and MIMs
with milder assumptions on both X and f in the high-dimensional setting. In this work,
we provide a partial answer to this question. We construct estimators that work for a wide
class of link functions, including the phase retrieval link function, and for a large family of
distributions of X, which is assumed to be known a priori. We particularly focus on the
case when X follows a non-Gaussian distribution that need not be elliptically symmetric or
sub-Gaussian, thus making our method applicable to several situations not possible before.
Our estimators are based on Stein’s identity for non-Gaussian distributions, which utilizes
the score function. Estimating with the score function is challenging due to their heavy tails.
In order to illustrate that, consider the univariate histograms provided in Figure-1. The dark
shaded, more concentrated one corresponds to the histogram of 10000 samples from Gamma
distribution with shape and scale parameters set to 5 and 0.2 respectively. The transparent
histogram corresponds to the distribution of the score function of the same Gamma dis7

1400
1200
1000
800
600
400
200
0
-10

0

10

20

30

40

50

Figure 1: Histogram of Score Function based on 10000 independent samples from the Gamma
distribution with shape 5 and scale 0.2. The dark histogram (we recommend the reader to
zoom in to notice it) concentrated around zero corresponds to the Gamma distribution and
the transparent histogram corresponds to the distribution of the score of the same Gamma
distribution.
tribution. Note that even when the actual Gamma distribution is well concentrated, the
distribution is the corresponding score function is well-spread and heavy-tailed. In the high
dimensional setting, in order to estimate with the score functions, we require certain vectors
or matrices based on the score functions to be well-concentrated in appropriate norms. In
order to achieve that, we construct robust estimators via careful truncation arguments to
balance the bias (due to thresholding)-variance (of the estimator) tradeoff and achieve the
required concentration.
• We construct estimators for the parametric component of a sparse SIM and MIM for
a class of unknown link function under the assumption that the covariate distribution
is non-Gaussian but known a priori.
• We establish near-optimal statistical rates for our estimators. Our results complement
the existing ones in the literature and hold in several case not possible before.
• We provide alternative justifications based on the Stein’s identity for the estimator
used in Plan and Vershynin (2016) for sparse SIM and the initializer used in Candès
et al. (2015) for phase retrieval.

8

• As a consequence of our results for SIM and MIM, we also obtain a near-optimal
estimator for sparse PCA with heavy-tailed data in the moderate sample size regime.
• We provide numerical simulations that confirm our theoretical results.

1.5

Notations

In this section, we introduce the notation and define the single index models. Throughout
this work, we use [n] to denote the set {1, . . . , n}. In addition, for a vector v ∈ Rd , we denote
by kvkp the `p -norm of v for any p ≥ 1. We use S d−1 to denote the unit sphere in Rd , which
is defined as S d−1 = {v ∈ Rd : kvk2 = 1}. In addition, we define the support of v ∈ Rd as
supp(v) = {j ∈ [d], vj 6= 0}. Moreover, we denote the nuclear norm, operator norm, elementwise max norm and Frobenius norm of a matrix A ∈ Rd1 ×d2 by k · k? , k · kop , k · k∞ and k · kfro ,
respectively. We denote by vec(A) the vectorization of matrix A, which is a vector in Rd1 ·d2 .
For two matrices A, B ∈ Rd1 ×d2 we define the trace inner product as hA, Bi = trace(A> B).
Note that it can be viewed as the standard inner product between vec(A) and vec(B). In
addition, for an univariate function g : R → R, we denote by g ◦ (v) and g ◦ (A) the output
of applying g to each element of a vector v and a matrix A, respectively. Finally, for a
random variable X ∈ R with density p, we use p⊗d : Rd → R to denote the joint density
of {X1 , · · · , Xd }, which are d identical copies of X. We also require some notations about
tensors. We concentrate on fourth-order tensors for simplicity. For any fourth-order tensor
Z ∈ Rd1 ×d2 ×d3 ×d4 , we denote its (j1 , j2 , j3 , j4 )-th entry by Z(j1 , j2 , j3 , j4 ). If d` = d for all
⊗

4

` ∈ [4], we denote the tensor as Z ∈ Rd 4 . Similar to the matrix case, we define vec(Z) ∈ Rd
⊗

as the vectorization of the tensor Z. For two tensors W, Z ∈ Rd 4 , we define their inner
inner product as
hZ, W i = vec(Z)> vec(W ) =

X

Z(j1 , j2 , j3 , j4 ) · W (j1 , j2 , j3 , j4 )

(1.2)

j1 ,j2 ,j3 ,j4 ∈[d]

The tensor Frobenius norm of is also denoted by k · kfro .

2

Index Models

Now we are ready to define the precise statistical models that we consider in this work. As
mentioned above, we consider the case of k = 1 (SIM) and k > 1 (MIM) separately. We
9

primarily distinguish our models based on the assumption made on the link functions. We
also require the following definition of score function of random variable. Let p : Rd → R be
a probability density function defined on Rd . The score function Sp : Rd → R associated to
p is defined as
Sp (x) = −∇x [log p(x)] = −∇x p(x)/p(x).
Note that in the above definition, the derivative is taken with respect to x. This is different
from the more traditional definition of the score function where the density belongs to a
parametrized family and the derivative is taken with respect to the parameters. In the rest
of the paper to simplify the notation, we omit the subscript x from ∇x . We also omit the
subscript p from Sp when the underlying density p is clear from the context.

2.1

First-order Link Functions

We first discuss a class of SIM that are based on a certain first-order link functions. We
discuss the motivation for our estimator, which automatically highlights the first-order assumption on the link function as well. Recall that our estimators are based on Stein’s
identity. To begin with, we present the first-order non-Gaussian Stein’s identity.
Proposition 2.1 (First-order Stein’s Identity (Stein et al., 2004)). Let X ∈ Rd be a realvalued random vector with density p. Assume that p : Rd → R is differentiable. In addition,
let g : Rd → R be a continous function such that E[∇g(X)] exists. Then it holds that
E[g(X) · S(X)] = E[∇g(X)],
where S(x) = −∇p(x)/p(x) is the score function of p.
One could apply the above Stein’s Identity to SIMs to obtain an estimate of β ∗ . To see
this, note that when X ∼ N (0, Id ) we have S(x) = x, ∀x ∈ Rd . In this case, as E() = 0, we
have
E(Y · X) = E[f (hX, β ∗ i) · X] = E[f 0 (hX, β ∗ i)] · β ∗ .
Hence one could estimate β ∗ based on estimating the moment E(Y · X). This observation
leads to the estimator proposed in Plan and Vershynin (2016). This motivates the following
definition of SIM with first order link functions.
10

Definition 2.2 (Vector SIM with First-order Links). Under this model, we assume that the
response variable Y ∈ R and the covariate X ∈ Rd are linked via
Y = f (hX, β ∗ i) + ,

(2.1)

where f : R → R is an unknown univariate function, β ∗ ∈ Rd is the parameter of interest,
and  ∈ R is the exogenous random noise such that E() = 0. In addition, we assume that
the entries of X are i.i.d. random variables with density p0 and that β ∗ is s∗ -sparse, that is,
β ∗ contains only s∗ nonzero entries such that s∗  n  d. Moreover, since the norm of β ∗
can be absorbed in f , we further let kβ ∗ k2 = 1 for identifiability. Finally, we assume f and
X are such that E[f 0 (hX, β ∗ i)] 6= 0.
Note that the SIM depends only on covariate only via inner products. Hence it is natural
to generalize it to the case of matrix and tensor valued covariates. To enable estimation in
a high-dimensional setting, we enforce low-rank constraints that we describe below.
Definition 2.3 (Matrix SIM with First-order Links). For the low-rank case SIM, we assume
that β ∗ ∈ Rd1 ×d2 has rank r∗  min{d1 , d2 }. In this scenario, X ∈ Rd1 ×d2 and the inner
product in (2.1) is hX, β ∗ i = trace(X > β ∗ ). For model identifiability, we further assume
that kβ ∗ kfro = 1, similar to the sparse case. Finally, we assume f and X are such that
E[f 0 (hX, β ∗ i)] 6= 0.
Before we lay out the first-order low-rank tensor single index model, we first introduce
⊗4

additional notation for tensors. Denote by u ⊗ v ⊗ s ⊗ t ∈ Rd

a rank-one tensor. The

minimum value of r such that the tensor Z could be written as a summation of r rank-one
P
tensors, i.e., Z = rj=1 uj ⊗ vj ⊗ sj ⊗ tj , is called as the CP-rank of the tensor, denoted by
rankCP (Z) = r. We now describe the low-rank tensor model that we consider in this work.
Definition 2.4 (Tensor SIM with First-order Links). For the low-rank Tensor SIM model,
⊗4

we assume that β ∗ ∈ Rd

⊗4

and has CP-rank, rankCP (β ∗ ) = r∗ . In this scenario, X ∈ Rd

and the inner product in (2.1) is understood as defined in (1.2). For model identifiability,
we further assume that kβ ∗ kfro = 1, similar to the matrix case. Finally, we assume f and X
are such that E[f 0 (hX, β ∗ i)] 6= 0.

11

2.2

Second-order Link Functions

In the above models, it is crucial that E[f 0 (hX, β ∗ i)] 6= 0, for it to work. Such a restriction
prevents it from being applicable to some widely used cases of SIM, for example, phase
retrieval where f is the quadratic function. This limitation of the first order Stein’s identity,
motivates us to examine the second order Stein’s identity which is summarized below.
Proposition 2.5 (Second-order Stein’s Identity (Janzamin et al., 2014)). Assume that the
density of X is twice differentiable. In addition, define the second-order score function
T : Rd → Rd×d as T (x) = ∇2 p(x)/p(x). Then for any twice differentiable function g : Rd →
R such that E[∇2 g(X)] exists, we have




E g(X) · T (X) = E ∇2 g(X) .

(2.2)

Back to the phase retrieval example, when X ∼ N (0, Id ), the second order score function
now becomes T (x) = xx> − Id , ∀x ∈ Rd . Setting g(x) = hx, β ∗ i2 in (2.2), we have
E[g(X) · T (X)] = E[g(X) · (XX > − I)] = E[hX, β ∗ i2 · (XX > − I)] = 2β ∗ β ∗ > .

(2.3)

Thus for phase retrieval, one could extract ±β ∗ based on second order Stein’s identity even
in the situation where the first order Stein’s identity fails. Indeed, (2.3) used in Candès
et al. (2015) implicitly to provide a spectral initialization for the Wirtinger flow algorithm
in the case of Gaussian phase retrieval. Here, we provided an alternative justification based
on Stein’s identity, for why such an initializer works. Motivated by the this observation, we
propose to use the second order Stein’s identity to estimate the parametric component of
SIMs and MIMs with a class of unknown link functions with non-Gaussian covariates. The
precise statistical models that we consider are defined as follows.
Definition 2.6 (Vector SIM with Second-order Links). Under this model, we assume that
the response variable Y ∈ R and the covariate X ∈ Rd are linked via
Y = f (hX, β ∗ i) + ,

(2.4)

where f : R → R is an unknown univariate function, β ∗ ∈ Rd is the parameter of interest,
and  ∈ R is the exogenous random noise such that E() = 0. In addition, we assume that
the entries of X are i.i.d. random variables with density p0 and that β ∗ is s∗ -sparse, that
12

is, β ∗ contains only s∗ nonzero entries. Moreover, since the norm of β ∗ can be absorbed in
f , we further let kβ ∗ k2 = 1 for identifiability. Finally, we assume f and X are such that
E[f 00 (hX, β ∗ i)] > 0.
Note that in the definition of the SIMs, we require that E[f 00 (hX, β ∗ i)] positive. Since if
E[f 00 (hX, β ∗ i)] is negative , we could always replace f by −f by flipping the sign of Y , we
essentially assume that E[f 00 (hX, β ∗ i)] is nonzero. Intuitively, such restriction on f implies
that the second order moments contains the information of β ∗ , thus we call such a function
the second order link. Similar to the first-order case, one could define matrix and tensor
versions of the second-order SIMs but we do not concentrate on such models in this work.
Thus far, we considered SIMs. We now define a class of MIMs with second order links.
For MIMs the notion of first order link functions is naturally not sufficient to estimate the
projector onto the subspace.
Definition 2.7 (MIM with Second-order Links). Under this model, we assume that the
response variable Y ∈ R and the covariate X ∈ Rd are linked via
Y = f (hX, β1∗ i, . . . , hX, βk∗ i) + ,

(2.5)

where f : Rk → R is an unknown function, {β`∗ }`∈[k] ⊆ Rd are the parameters of interest,
and  ∈ R is the exogenous random noise such that E() = 0. In addition, we assume
that the entries of X are i.i.d. random variables with density p0 and that {β`∗ }`∈[k] span a
k-dimensional subspace of Rd . Moreover, we denote B ∗ = (β1∗ . . . βk∗ ) ∈ Rd×k . Then the
model in (2.5) can be written as Y = f (XB ∗ ) + . By the QR-factorization, we can write
B ∗ as Q∗ R∗ , where Q∗ ∈ Rd×k is an orthonormal matrix and R∗ ∈ Rk×k is invertible.
Since f is unknown, R∗ can be absorbed into the link function. Thus, we assume that B ∗
is orthonormal for identifiability. Furthermore, we further assume that B ∗ is s∗ -row sparse,
that is, B ∗ contains only s∗ nonzero rows. We note that such a definition of sparsity for B ∗
does not depends on the choice of coordinate system. Finally, we assume f and X are such
that λmin (E[∇2 f (XB ∗ )]) > 0.
The assumption E[∇2 f (XB ∗ )] is positive definite, in MIM, is a multivariate generalization of the condition that E[f 00 (hX, β ∗ i)] > 0 in SIM. It essentially guarantees that estimation
of the projector onto the subspace spanned by the k components is well-defined. We now
13

introduce our estimators and provide theoretical results that are near-optimal in several
settings.

3

Theoretical Results for Index Models with First-order
Links

Recall that in the single index models introduced in §2.1, X in (2.1) has i.i.d. entries with
density p0 . To unify the vector, matrix and tensor settings, we identify X with vec(X) ∈ Rd
where d = d1 · d2 · d3 · d4 . In this case, X has density p = p⊗d
0 and the corresponding score
function S : Rd → Rd is given by
S(x) = −∇ log p(x) = −∇p(x)/p(x) = s0 ◦ (x),

(3.1)

where the univariate function s0 = p00 /p0 is applied to each entry of x. Thus S(X) has i.i.d.
entries. In addition, by Lemma 2.1, we have E[S(X)] = 0 by setting g to be a constant
function. Moreover, in the context of SIMs specified in (2.1), we have


E[Y · S(X)] = E f (hX, β ∗ i) · S(X) = E[f 0 (hX, β ∗ i)] · β ∗ ,
as long as the density and the link function satisfy the conditions stated in Lemma 2.1. This
implies that optimization problem

minimize hβ, βi − 2E[Y · hS(X), βi]
β∈Rd

(3.2)

has solution β = µ · β ∗ , where µ = E[f 0 (hX, β ∗ i)]. Hence the above program could be used to
obtain the unknown β ∗ as long as µ 6= 0. Before we proceed to describe the sample version
of the above program, we make the following brief remark. The requirement µ 6= 0 rules out
in particular the use of our approach for non-Gaussian phase retrieval (where f (u) = u2 ) as
in that case we have µ = 0 when X is centered. But we emphasize that the same holds true
in the Gaussian and elliptical setting as well, as noted in Plan and Vershynin (2016) and
Goldstein et al. (2016). Their methods also fail to recover the true β ∗ when the SIM model
corresponds to phase retrieval. We refer the reader to §4 for overcoming this limitation using
second-order Stein’s identity.
We use a sample version of the above program as an estimator for the unknown β ∗ . In
order to deal with the high-dimensional setting, we consider a regularized version of the
14

above formulation. More specifically, we use the `1 -norm and nuclear norm regularization
in the vector and matrix/tensor settings respectively. However, a major difficulty in the
sample setting for this procedure is that E[Y · S(X)] and its empirical counterpart may
not be close enough due to a lack of concentration. Recall our discussion from §1 that
even if the random variable X is light-tailed, its score-function S(x) might be arbitrarily
heavy-tailed. Furthermore, bounded-fourth moment assumption on the noise, Y too can be
heavy-tailed. Thus the naive method of using the sample version of (3.2) to estimate β ∗
leads to sub-optimal statistical rates of convergence.
To improve concentration and obtain optimal rates of convergence, we replace Y · S(X)
with a transformed random variable T (Y, X), which will be defined precisely later for the
sparse and low-rank cases. In particular, T (Y, X) is a carefully truncated version of Y ·S(X),
introduced and analyzed in Catoni et al. (2012); Fan et al. (2016) for related problems, that
enables us to obtain well-concentrated estimators. Thus our final estimator βb is defined as
the solution to the following regularized optimization problem
n

2X
hT (Yi , Xi ), β ,
minimize L(β) + λ · R(β), L(β) = hβ, βi −
n i=1
β∈Rd

(3.3)

where λ > 0 is the regularization parameter which will be specified later and R(·) is the `1 norm in the vector case and the nuclear norm in the matrix/tensor case. We now introduce
our main moment assumption for first-order SIM. This assumption is made apart from the
assumptions made on the noise and the link function. Recall that each entry of the score
function defined in (3.1) is equal to s0 (u) = −p00 (u)/p0 (u). We first state the assumption
and make a few remarks about it.
Assumption 3.1 (Moment Assumptions). There exists an absolute constant M > 0 such
that E(Y 4 ) ≤ M and Ep0 [s40 (U )] ≤ M , where random variable U ∈ R has density p0 .
Consider the assumption E(Y 4 ) ≤ M . By Cauchy-Schwarz inequality we have E(Y 4 ) ≤
4E(4 ) + 4E[f 4 (hX, β ∗ )]. Note that we assum  to be centered, independent of X and has
bounded fourth moment (see §2). If the covariate X has bounded fourth moment along the
direction of true parameter, since f (·) is continuously differentiable, f (hX, β ∗ i) has bounded
fourth moment as well if f (·) is defined on a compact subset of R. . Hence the condition
E(Y 4 ) ≤ M is relatively easy to satisfy and significantly milder than assuming that Y is
bounded or has lighter tails. Furthermore, Ep0 [s40 (U )] ≤ M is relatively mild and satisfied by
15

a wide class of random variables. Specifically random variables that are non-symmetric and
non-Gaussian satisfy this property thereby allowing our approach to work with covariates
not previously possible. We believe it is highly non-trivial to weaken this condition without
losing significantly in the rates of convergence that we discuss below.

3.1

Sparse Vector SIM

Under the above assumptions, we first state our theorem on the sparse SIM. As discussed
above, Y ·S(X) can by heavy-tailed and hence we apply truncation to achieve concentration.
Denote the j-th entry of the score function S in (3.1) as Sj : Rd → R, j ∈ [d]. We define the
truncated response and score function as


Ye = sign(Y ) · (|Y | ∧ τ ), Sej (x) = sign[Sj (x)] · |Sj (x)| ∧ τ ,

(3.4)

where τ > 0 is a predetermined threshold value. We define Yei similarly for all Yi , i ∈ [n].
Then we define the estimator βb as the solution to the optimization problem in (3.3) with
e i ) and R(β) = kβk1 . Here we apply elementwise truncation in T to ensure
T (Yi , Xi ) = Yei ·S(X
the sample average of T converges to E[Y ·S(X)] in the `∞ -norm for an appropriately chosen
τ . Note that the `∞ -norm is the dual norm of the `1 -norm. Such a convergence requirement
in the dual norm is standard in the analysis of regularized M -estimators (Negahban et al.,
2012) to achieve optimal rates. The following theorem characterizes the convergence rates
b
of β.
Theorem 3.2 (Signal Recovery for Sparse Vector SIM). For the sparse SIM defined in §2,
we assume that β ∗ ∈ Rd has s∗ nonzero entries. Under Assumption 4.1, we let
τ = 2(M · log d/n)1/4
in (3.4) and set the regularization parameter λ in (3.3) as
λ=C

p
M · log d/n,

where C > 0 is an absolute constant. Then with probability at least 1−d−2 , the `1 -regularized
estimator βb defined in (3.3) satisfies
kβb − µβ ∗ k2 ≤

√

s∗ · λ, kβb − µβ ∗ k1 ≤ 4s∗ · λ.
16

p
From this theorem, the `1 - and `2 -convergence rates of βb are kβb−µβ ∗ k1 = O(s∗ log d/n)
p
and kβb − µβ ∗ k2 = O( s∗ log d/n), respectively. These rates match the convergence rates of
sparse generalized linear models (Loh and Wainwright, 2015) and sparse single index models
with Gaussian and symmetric elliptical covariates (Plan and Vershynin, 2016; Goldstein
et al., 2016) which are known to be minimax-optimal for this problem via matching lower
bounds.

3.2

Low-rank Matrix SIM

We next state our theorem for the low-rank Matrix SIM. In this case, we apply the nuclear
norm regularization to promote low-rankness. Note that by definition, T is matrix-valued.
Since the dual norm of the nuclear norm is the operator norm, we need the sample average
of T to converge to E[Y · S(X)] in the operator norm rapidly to achieve optimal rates of
convergence. To achieve such a goal, we leverage the truncation argument from Catoni et al.
(2012); Minsker (2016) to construct T (Y, X).
Let φ : R → R be a non-decreasing function such that
− log(1 − x + x2 /2) ≤ φ(x) ≤ log(1 + x + x2 /2), ∀x ∈ R.
Based on φ, we define a linear mapping ψ : Rd1 ×d2 → Rd1 ×d2 as follows. For any A ∈ Rd1 ×d2 ,
let


0
e=
A
A>


A

0

e In addition, let B = Υψ ◦ (Λ)Υ> , where
and let ΥΛΥ> be the eigenvalue composition of A.
ψ is applied elementwisely on Λ. Then we write B in block from as


B11 B12

B=
B21 B22


and define ψ(A) = B12 . Finally, we define T (Y, X) = 1/κ · ψ κ · Y · S(X) , where κ > 0 will
be specified later. Therefore, our final estimator βb ∈ Rd1 ×d2 is defined as the solution to the
optimization problem in (3.3) with R(β) = kβk? . We note here the minimization in (3.3) is
taken over Rd1 ×d2 . The following theorem quantifies the convergence rates of the proposed
estimator.
17

Theorem 3.3 (Signal Recovery for Low-rank Matrix SIM). For the low-rank single index
model defined in §2, we assume that rank(β ∗ ) = r∗ . Under Assumption 4.1, we let
p
p
κ = 2 n · log(d1 + d2 )/ (d1 + d2 )M
in T (Y, X) and set λ in (3.3) as
λ=C

p
M · (d1 + d2 ) · log(d1 + d2 )/n,

where C > 0 is an absolute constant. Then with probability at least 1 − (d1 + d2 )−2 , the
nuclear norm regularized estimator βb satisfies
√
kβb − µβ ∗ kfro ≤ 3 r∗ · λ, kβb − µβ ∗ k? ≤ 12r∗ · λ.
p
By this theorem, we have kβb − µβ ∗ kfro = O( r∗ (d1 + d2 ) · log(d1 + d2 )/n) and kβb −
p
µβ ∗ k? = O(r∗ · (d1 + d2 ) · log(d1 + d2 )/n). Note that the rate obtained is minimax-optimal
up to a logarithmic factor. Furthermore, it matches the rates for low-rank single index
models with Gaussian and symmetric elliptical distributions up to a logarithmic factor Plan
and Vershynin (2016); Goldstein et al. (2016).

3.3

Low-rank Tensor SIM

We now state our result for low-rank tensor SIM. The notion of rank of a tensor is more
delicate compared to that of a matrix. Several generalizations of the matrix rank exist for
the case of tensors. Recall from Definition 2.4, that we assumed that the structure on β ∗
is that it has low CP-rank. Unfortunately, enforcing such a constraint via a direct tensor
nuclear norm relaxation (similar to that of the matrix nuclear norm) is NP-hard (Friedland
and Lim, 2014).
One way to overcome such a computational hurdle is to deal with tensors via appropriately matricized forms. In order to enable computable estimators, we specifically
leverage the results of Mu et al. (2014) and define the following square-unfolding of a
⊗4

tensor. Denote by Mat : Rd

2 ×d2

→ Rd

the operation of tensor square-unfolding, which

maps a fourth-order tensor to a square matrix. More specifically, the entries of Mat(Z)
are specified by [Mat(Z)]k1 ,k2 = Z(j1 , j2 , j3 , j4 ), where the indices satisfy the relationship
k1 = 1 + (j1 − 1) + (j2 − 1) · d and k2 = 1 + (j3 − 1) + (j4 − 1) · d. Intuitively, the matrix
18

obtained by the square-unfolding operation is as square as possible, i.e., it is d2 × d2 rather
than the rectangular d × d3 or d3 × d matrices. It is shown in Mu et al. (2014) such a square
matricization preserves the low CP-rank of the original tensor. Hence one could use the
matrix nuclear norm relaxation on the square-unfolded tensor. Furthermore, for the case
of tensor SIM as in Definition 2.4, note that we have hX, β ∗ i = hMat(X), Mat(β ∗ )i. Combining the above observations, the low CP-rank tensor SIM problem could be reduced to
that of low-rank matrix SIM problem, where matrix low-rank constraint, via nuclear norm,
is enforced on Mat(β ∗ ). Thus, we use the estimator in (3.3) with R(β) = k Mat(β)k∗ and
Mat(Xi ) for all i = 1, . . . , n along with the truncation operation T described in §3.2. We
now have the following theorem for the low-rank tensor SIM.
Theorem 3.4 (Signal Recovery for Low-rank Tensor SIM). For the low-rank tensor single
index model defined in §2, Definition 2.4, we assume that rankCP (β ∗ ) = r∗ . Under Assumption 4.1, we let
p
p
κ = 2 2n · log d/ (2d2 )M
in T (Y, X) and set λ in (3.3) as
λ=C

p
2M · (2d2 ) · log d/n,

where C > 0 is an absolute constant. Then with probability at least 1 − 2d, the nuclear-norm
regularized estimator βb satisfies
√
kβb − µ Mat(β ∗ )kfro ≤ 3 r∗ · λ
We omit the proof of the above theorem as it is follows the exact steps of Theorem 3.3
proved in Appendix A.2. From the above theorem, we see that as long as n = Ω(r∗ d2 ),
we achieve consistent estimation of β ∗ up to scaling. This improves upon recent results
established in Chen et al. (2016), that established similar results under restrictive Gaussian
covariate assumption and required knowledge of the link functions (i.e., generalized linear
models). Furthermore our results significantly generalizes the results of Mu et al. (2014)
that considered only linear link functions. Finally, although our structure on β ∗ was a
low CP-rank structure, the square matricization technique also applies for the case of low
Tucker-rank, which is yet another notion of rank for tensors with several applications. It is
straightforward to extend our results to this case of low Tucker-rank.
19

4

Theoretical Results for Index Models with SecondOrder Links

We now introduce our estimators and establish their statistical rates of convergence for the
case of index models with second-order link functions. Discussions on optimality of the
established rates and connection to sparse PCA problem is deferred to §4.3. Similar to the
first-order case, we focus on the case where X has i.i.d. entries with density p0 : R → R.
Qd
Thus the joint density of X is p(x) = p⊗d
0 (x) =
j=1 p0 (xj ). We define a univariate function
s0 : R → R by s0 (u) = p00 (u)/p0 (u). Then the first-order score function associated with p
is given by S(x) = s0 ◦ (x). Equivalently, the j-th entry of the first-order score function
associated with p is given by [S(x)]j = s0 (xj ). Moreover, the second order score function is
T (x) = S(x)S > (x) − ∇S(x) = S(x)S > (x) − diag[s00 ◦ (x)].

(4.1)

Before we present our estimator, we introduce the assumption on Y and s0 (·).
Assumption 4.1 (Moment Assumptions). We assume that there exists a constant M such
that Ep0 [s60 (U )] ≤ M and E(Y 6 ) ≤ M . We denote σ02 = Ep0 [s20 (U )] = Varp0 [s0 (U )].
The assumption that Ep0 [s60 (U )] ≤ M allows wide family of distributions of including
Gaussian and more heavy-tailed random variables. Furthermore, we do not require the
covariate to be elliptically symmetric as is commonly seen in prior work, which enables our
method to be applicable for skewed covariates. As for the assumption that E(Y 6 ) ≤ M , in the
case of SIMs, we have E(Y 6 ) ≤ C (E(6 ) + E[f 6 (hX, β ∗ i)]). Thus this assumption is satisfied
as long as both  and f (hX, β ∗ i) have bounded sixth moments. This is a significantly milder
assumption which allows for heavy-tailed response as opposed to bounded or light-tailed
response.

4.1

Sparse Vector SIM

Now we are ready to describe our estimator for the sparse SIMs in Definition 2.6. Note that
by Proposition 2.5 we have


E Y · T (X) = C0 · β ∗ β ∗ > ,

20

(4.2)

where C0 = 2E[f 00 (hX, β ∗ i)] > 0 as per Definition 2.6. Therefore, one way to estimator β ∗
is to obtain the leading eigenvector of E[Y · T (X)] from the samples. Since β ∗ is sparse, we
formulate our estimator as a sparsity constrained semi-definite program:
maximize hW, Σ∗ i + λkW k1
subject to 0  W  Id , trace(W ) = 1.

(4.3)

where Σ∗ = E[Y · T (X)]. Note that both the score T (X) and the response variable Y can
be heavy-tailed. In order to obtain near-optimal estimates in the sample setting, we apply
truncation to handle the heavy-tails. Specifically, for a positive parameter τ ∈ R, we define
the truncated random variables by

Yei = sign(Yi ) · min{|Yi |, τ } and Tejk (Xi ) = sign{Tjk (Xi )} · min |Tjk (Xi )|, τ 2 .

(4.4)

Then we define an robust estimator of Σ∗ as
n

X
e= 1
Yei · Te(Xi ).
Σ
n i=1

(4.5)

e let W
c be the solution of the following convex optimization problem
Given Σ,
e + λkW k1
maximize hW, Σi
subject to 0  W  Id , trace(W ) = 1.

(4.6)

Here λ is a regularization parameter to be specified later. The final estimator is defined
c . The following theorem quantifies the statistical rates of
as the leading eigenvector of W
convergence of the proposed estimator.
c be the solution of the optimization
Theorem 4.2 (Signal Recovery for Spare SIM). Let W
c . We set the regularization parameter λ in
in (4.6) and let βb be the leading eigenvector of W
p
(4.6) as λ = 10 M log d/n and set τ = (1.5M n/ log d)1/6 in (4.4). Under Assumption 4.1,
√
we have kβb − β ∗ k2 ≤ 4 2s∗ λ with probability at least 1 − d−2 .
By this theorem, the `2 -error of the proposed estimator is O(s∗
that consistent estimation requires n = Ω(s∗ 2 log d) samples.

21

p
log d/n), which implies

4.2

Subspace-Sparse MIM

Now we introduce the estimator for B ∗ of the sparse MIM in Definition 2.7. Proposition 2.5
implies that E[Y · T (X)] = B ∗ D0 B ∗ , where D0 = E[∇2 f (XB ∗ )] is positive definite. Similar
to (4.6), we recover the column space of B ∗ by solving
e + λkW k1 ,
maximize hW, Σi
subject to 0  W  Id , trace(W ) = k.

(4.7)

e is defined in (4.5), λ > 0 is a regularization parameter and k is the number of indices
where Σ
c be the solution of (4.7), the final estimator is the top
which is assumed to be known. Let W
c . For the above estimator, we have the following theorem quantifying
k eigenvectors of W
the statistical rate of convergence. Let ρ0 = λmin (E[∇2 f (XB ∗ )]).
c be the solution of the optimization
Theorem 4.3 (Signal Recovery for Sparse MIM). Let W
b be the leading eigenvector of W
c . We set the regularization
problem in (4.7) and let B
p
parameter in (4.6) as λ = 10 M log d/n and let the truncation parameter in (4.4) be τ =
(1.5M n/ log d)1/6 . Under Assumption 4.1, with probability at least 1 − d−2 , we have
√
b − B ∗ Ok2 ≤ 4 2/ρ0 · s∗ λ.
inf kB
O∈Ok

Minimax lower bounds for subspace estimation for MIM was established recently in Lin
et al. (2017). For a fixed k, the above theorem is near-optimal from a minimax point of view.
√
That is, the difference between the optimal rate and the above theorem is a factor of s.
We discuss more about this gap in Section 4.3. The proofs of Theorem 4.2 and Theorem 4.3
are in the supplementary material.
Remark 1. Recall that our discussion in §3 and §4 was under the assumption that the entries
in X are i.i.d. This could be relaxed to the case of weak dependence between the covariates
without any significant loss in the statistical rates we present in the theorems above. We
do not focus on such an extension in this paper as we wanted to clearly convey the main
message of the paper in a simpler setting.

4.3

Optimality and Relation to Sparse PCA

Now we discuss the optimality of the results presented in §4. Throughout the discussion we
assume that k is fixed and does not increase with n. Note that the estimator for SIM in
22

(4.6) and MIM in (4.7) are closely related to the semidefinite program based estimator for
Sparse PCA problem (Vu et al., 2013). Let X ∈ Rd be a random vector such that E(X) = 0
and covariance matrix Σ = E(XX > ) which is symmetric and positive definite. The problem
of sparse PCA is to estimate projector onto the subspace spanned by top k eigenvectors,
{v`∗ }`∈[k] of Σ under the subspace sparsity assumption as discussed in Definition 2.7. An
estimator based on semidefinite programing with sparsity constraints was analyzed in Vu
et al. (2013); Wang et al. (2016), which is based on solving the following program
b + λkW k1
maximize hW, Σi
subject to 0  W  Id , trace(W ) = k.

(4.8)

b = n−1 Pn Xi X > is the sample covariance matrix given n i.i.d copies {Xi }n of X.
Here Σ
i=1
i
i=1
Note that the main difference between the SIM estimator and the sparse PCA estimator is the
e in place of Σ.
b It is known that sparse PCA problem exhibits interesting statisticaluse of Σ
computational tradeoff (Krauthgamer et al., 2015; Wang et al., 2016) which naturally appears
in the context of SIM as well. Indeed while the minimax optimal statistical rate for sparse
p
p
PCA is O( s∗ log d/n), the SDP estimator achieves O(s∗ log d/n) under the assumption
that X is light-tailed. It is also known that when n = Ω(s∗ 2 log d), one can obtain the
p
optimal statistical rate of O( s∗ log d/n) either by nonconvex methods (Wang et al., 2014),
or refinements to the output of the SDP estimator (Wang et al., 2016). However their results
b to Σ in the restricted operator norm:
rely on the sharp concentration of Σ
p

b − Σ∗ kop,s = sup w> (Σ
b − Σ)w : kwk2 = 1, kwk0 ≤ s = O( s log d/n).
kΣ

(4.9)

When X has heavy-tailed entries, for example bounded fourth moment assumptions, its
highly unlikely that, (4.9) holds. Indeed the results in Wang et al. (2016) and Wang et al.
(2014) are applicable only to the case of Gaussian or light-tailed X.
4.3.1

Heavy-tailed Sparse PCA

Recall that our estimators utilize a data-driven truncation argument to handle heavy-tailed
distributions. Owing to the close relationship between our SIM/MIM estimators and the
sparse PCA estimator, it is natural to ask whether such a truncation argument could lead
to sparse PCA estimators for heavy tailed X. Below we show that it is indeed possible
23

to obtain a near-optimal estimator for sparse PCA with heavy-tailed data based on the
truncation argument. For a vector v ∈ Rd , let ϑ(v) be a truncation operation that operators
entry-wise as ϑj (v) = sign[vj ] · min {|vj |, τ } for j = 1, . . . d. Then, our estimator is defined
as follows.
maximize hW, Σi + λkW k1
subject to 0  W  Id , trace(W ) = k.
where Σ = n−1

Pn

i=1

(4.10)

X iX >
i and X i = ϑ(Xi ), for i = 1, . . . n. For the above estimator, we

have the following theorem under the assumption that X has heavy-tailed marginals. Let
V ∗ = (v1∗ . . . vk∗ ) ∈ Rd×k and assume that ρ0 = λk (Σ) − λk+1 (Σ) > 0.
c be the solution of the optimization in (4.10) and let Vb be the leading
Theorem 4.4. Let W
p
c . We set the regularization parameter in (4.10) as λ = C1 M log d/n and
eigenvector of W
set the truncation parameter by τ = (C2 M n/ log d)1/4 , where C1 and C2 are some positive
constants. Furthermore, assume that V ∗ contains only s∗ nonzero rows and that X satisfies
E|Xj |4 ≤ M and E|Xi · Xj |2 ≤ M . Then, with probability at least 1 − d−2 , we have
√
inf kVb − V ∗ Ok2 ≤ 4 2/ρ0 · s∗ λ.

O∈Ok

The proof of the above theorem is similar to that of Theorem 4.3 and hence we omit
it. The above theorem shows that with elementwise truncation, as long as X satisfies a
bounded fourth moment condition, the SDP estimator for sparse PCA achieves the nearp
optimal statistical rate of O(s∗ log d/n). We end this section with the following questions
based on the above discussions:
p
1. Can we obtain optimal statistical rate for sparse PCA problem (O( s log d/n)) when
X has only bounded fourth moment in the high sample size regime n = Ω(s∗ 2 log d) ?
p
2. Can we obtain optimal statistical rate (O( s∗ log d/n)) when n = Ω(s∗ 2 log d) and
when f, X and Y satisfies the heavy-tail condition in Assumption 4.1 for the MIM
problem ?
The answer to both questions lie in constructing truncation based estimators that concentrate sharply in restricted operator norm as in (4.9) or more realistically exhibit one-sided
24

concentration bounds (see Mendelson (2014) and Oliveira (2013) for related results and discussion). Obtaining such an estimator seems to be challenging for heavy-tailed sparse PCA
and it it not immediately clear if it is possible. We plan to report our findings for the above
problem in the near future.

5

Numerical Experiments

We now provide simulation experiments for the case of first-order and second-order SIMs.
For the first-order SIM, we concentrate on the sparse vector and low-rank matrix model.
Note that our tensor estimator is similar to the low-rank matrix estimator. Furthermore, for
the second-order case, we concentrate on the problem of robust sparse phase retrieval.
First-order SIM: We let  ∼ N (0, 1) and set the link function in (2.1) as one of
√
f1 (u) = 3u + 10 sin(u) and f2 (u) = 2u + 4 exp(−2u2 ), which are plotted in Figure 2. We
set p0 to be one of (i) Gamma distribution with shape parameter 5 and scale parameter
1, (ii) Student’s t-distribution with 5 degrees of freedom, and (iii) Rayleigh distribution
with scale parameter 2. To measure the estimation accuracy, we use the cosine distance
b β ∗ i|, where • stands for the Euclidean norm in the vector case
b β ∗ ) = 1 − kβk
b −1 |hβ,
cos θ(β,
•
and the Frobenius norm when β ∗ is a matrix. Here we report the cosine distance rather than
kβb − µβ ∗ k• to compare the performances for X having different distributions, where µ may
have different values.
For the vector case, we fix d = 2000, s∗ = 5 and vary n. The support of β ∗ is chosen
√
uniformly random among all subsets of {1, . . . , d}. For each j ∈ supp(β ∗ ), we set βj∗ = 1/ s∗ ·
γj , where each γj is an i.i.d. Rademacher random variable. In addition, the regularization
p
parameter λ is set to 4 log d/n. We plot the cosine distance against the signal strength
p
s∗ log d/n in Figure 5-(a) and (b) for f1 and f2 respectively, based on 200 independent
trials for each n. As shown in this figure, the estimation error grows sub-linearly as a function
of the signal strength.
As for the matrix case, we fix d1 = d2 = 20, r∗ = 3 and let n vary. The signal parameter β ∗ is equal to U SV > , where U, V ∈ Rd×d are random orthogonal matrices and S
is a diagonal matrix with r∗ nonzero entries. Moreover, we set the nonzero diagonal en√
tries of S as 1/ r∗ , which implies kβ ∗ kfro = 1. We set the regularization parameter as
p
λ = 2 (d1 + d2 ) log(d1 + d2 )/n. Furthermore, we use the proximal gradient descent algo25

rithm (with the learning rate fixed to 0.05) to solve the nuclear norm regularization problem in (3.3). To present the result, we plot the cosine distant against the signal strength
p
r∗ (d1 + d2 ) log(d1 + d2 )/n in Figure 5-(b) based on 200 independent trials. As shown in
this figure, the error is bounded by a linear function of the signal strength, which corroborates
Theorem 3.3.
40

15

2a + 4 exp(!2a2 )

20
10

-10

p

0

-20

f (u) =

f (u) = 3 " u + 10 sin(u)

30

-30
-40
-10

-5

0

5

10

10
5
0
-5
-10
-15
-10

u

-5

0

5

10

u

Figure 2: Plot of the link functions f1 (u) = 3u + 10 · sin(u) and f2 (u) =

√

2u + 4 exp(−2u2 ).

Second-order SIM: We now concentrate on the problem of sparse phase retrieval using
the SDP based estimators proposed based on second-order Stein’s identity. Recall that in
this case, the link function is known and existing convex and non-convex based estimators
are applicable predominantly for the case of Gaussian or light-tailed data. The question
of de-randomization or what are the necessary assumptions on the measurement vectors for
(sparse) phase retrieval to work is an intriguing one (Gross et al., 2015). Here we demonstrate
that using the proposed score-based estimators, one could deal with heavy-tailed and skewed
measurement as well, which significantly extend the class of measurement vectors applicable
for sparse phase retrieval.
Recall that the covariate X has i.i.d. entries with distribution p0 . We set p0 to be
one of Gamma distribution with shape parameter 5 and scale parameter 1 and Rayleigh
distribution with scale parameter 2. The random noise  is set to be standard Gaussian.
Moreover, we solve the optimization problems in (4.6) and (4.7) via the alternating direction
method of multipliers (ADMM) algorithm proposed in Vu et al. (2013), which introduces a
dual variable to encode the constrains and updates the primal and dual variables iteratively.
26

0.3

0.1
Gamma(5; 2)
t(5)
Rayleigh(2)

0.25

0.09
0.08

Gamma(5; 2)
t(5)
Rayleigh(2)

0.07

b -$)
cos 3(-;

b -$)
cos 3(-;

0.2

0.06

0.15

0.05
0.04

0.1

0.03
0.02

0.05

0.01
0
0.1

0.2

p 0.3
s$ log d=n

0.4

0
0.1

0.5

0.2

p 0.3
s$ log d=n

0.4

0.5

Figure 3: Cosine distances between the true parameter and the estimated parameter in the
sparse SIM with the link function in 2.1 set to one of f1 and f2 . Here we set d = 2000.
s∗ = 5 and vary n.
We set the link function to be one of f3 (u) = u2 , f4 = |u|, and f5 (u) = 4u2 + 3 cos(u).
Here f3 corresponds to the phase retrieval model and f4 and f5 can be viewed as its robust
extension. Throughout the experiment we fix d = 500, s∗ = 5 and vary n. The support of
β ∗ is chosen uniformly random among all subsets of {1, . . . , d}with cardinality s∗ . For each
√
j ∈ supp(β ∗ ), we set βj∗ = 1/ s∗ · γj , where γj ’s are i.i.d. Rademacher random variables.
p
Furthermore, we fix the regularization parameter λ = 4 log d/n and threshold parameter
b β ∗ ) = 1 − |hβ,
b β ∗ i|, to measure
τ = 20. In addition, we adopt the cosine distance cos θ(β,
the estimation error. We plot the cosine distance against the statistical rate of convergence
p
s∗ log d/n in Figure 5-(a)-(c) for each link function, respectively. The plot is based on 100
independent trials for each n, which shows that the estimation error is bounded by a linear
p
function of s∗ log d/n, which corroborate the theory.

6

Conclusion

In this work, we consider estimating the parametric components of single and multiple index
models in the high-dimensional setting, under fairly general assumptions on the link function
f and response Y . Furthermore, our estimators are applicable in the non-Gaussian setting
where X is not required to satisfy restrictive Gaussian or elliptical symmetry assumptions.
27

0.6

0.6
Gamma(5; 2)
t(5)
Rayleigh(2)

0.5

0.5

b -$)
cos 3(-;

0.4

b -$)
cos 3(-;

0.4

Gamma(5; 2)
t(5)
Rayleigh(2)

0.3

0.3

0.2

0.2

0.1

0.1

0
0.1

0.2 p

0.3

0.4

0.5

0.6

0
0.1

0.7

r$ (d1 + d2 ) log(d1 + d2 )=n

0.2 p

0.3

0.4

0.5

0.6

0.7

r$ (d1 + d2 ) log(d1 + d2 )=n

Figure 4: Cosine distances between the true parameter and the estimated parameter in the
low-rank SIM for with link function in 2.1 set to one of f1 and f2 . Here we set d1 = d2 = 20.
r∗ = 3 and vary n.
Our estimators are based on a data-driven truncation argument in combination with first
and second-order Stein’s identity. Furthermore, we show that proposed estimators are nearoptimal for several different settings.
Recently in the low-dimensional setting, for 2-layer neural networks Janzamin et al.
(2015) proposed a tensor-based method for estimating the parametric components. Their
estimators are sub-optimal even when we consider the low-dimensional Gaussian setting.
An immediate application of our truncation based estimators enables us to obtain optimal
results for a fairly general class of covariates in the low-dimensional setting. Obtaining
similar optimal or near-optimal results in the high-dimensional setting is of great interest
for 2-layer neural networks, albeit challenging. We plan to extend the result of this paper
for 2-layer neural networks in the high-dimensional setting and report our results in the near
future.

A

Proofs of the Main Results

In this section, we lay out the proofs of the theorems in §3 and §4, which establish the
statistical rates of convergence of our estimators.

28

0.3

0.3

0.3

0.25

0.25

0.25

0.2

0.2

0.2

0.15

0.15

0.15

0.1

0.1

0.1

0.05

0.05

0.05

0
0.1

0.15

0.2

0.25

0.3

0.35

0
0.1

0.4

f3 (u) = u2

0.15

0.2

0.25

0.3

f4 (u) = |u|,

0.35

0.4

0
0.1

0.15

0.2

0.25

0.3

0.35

0.4

f5 (u) = 4u2 + 3 cos(u)

Figure 5: Cosine distances between the true parameter β ∗ and the estimated parameter βb in
the sparse SIM with the link function in one of f1 , f2 , and f3 . Here we set d = 500. s∗ = 5
and vary n.

A.1

Proof of Theorem 3.2

Proof. Since βb is the solution of the optimization problem in (3.3), the first-order optimality
condition states that
b + λξ = 0, where ξ ∈ ∂kβk
b 1.
∇L(β)

(A.1)

Then the entries of ξ ∈ Rd are given by
b ξj ∈ [−1, 1], ∀j ∈
b
ξj = sign(βbj ), ∀j ∈ supp(β);
/ supp(β).
For any index set A ⊆ [d] and z ∈ Rd , we define the restriction of z to A, zA ∈ Rd , by letting
[zA ]j = zj if j ∈ A, [zA ]j = 0 otherwise.
Here [zA ]j is the j-th entry of zA . Let S = supp(β ∗ ), then we can write ξ = ξS + ξS c . For
notational simplicity, in the sequel, we define θ = βb − µ · β ∗ . Thus by (A.1) it holds that
b − ∇L(µβ ∗ ), θi = h−λ · ξ − ∇L(µβ ∗ ), θi
h∇L(β)
≤ h−λ · ξS − λ · ξS c , θi + k∇L(µβ ∗ )k∞ · kθk1 .

(A.2)

By the definition of ξ, we have
h−λ · ξS c , βb − µβ ∗ i = −λ · kβbS k1 .

29

(A.3)

Moreover, since kξk∞ ≤ 1, Hölder’s inequality implies that
h−λ · ξS , θi ≤ kθS k1 .

(A.4)

Note that ∇2 L(β) = 2Id . Combining (A.9), (A.3), and (A.4), we obtain
b − ∇L(µβ ∗ ), θi ≤ −λkθS c k1 + λkθS k1 + k∇L(µβ ∗ )k∞ · kθk1 .
2kθk22 = h∇L(β)

(A.5)

For an upper bound of the right-hand side of (A.5), we apply the following lemma to obtain
an upper bound on k∇L(µβ ∗ )k∞ .
Lemma 1 (Bound on k∇L(µβ ∗ )k∞ ). We set the truncation level in (3.4) as τ = 2(M ·
n/ log d)1/4 . Then we have
h
i
p
P k∇L(µβ ∗ )k∞ > 7 M · log d/n ≤ d−2 .
Proof. See §B.1 for a detailed proof.
Thus by Lemma 1 and the choice of λ, we have λ > 2k∇L(µβ ∗ )k∞ with probability at
least 1 − d−2 . This implies that
2kθk22 ≤ −λ/2 · kθS c k1 + 3λ/2 · kθS k1 ≤ 2λ · kθS k1 .

(A.6)

Since the leftmost term in (A.6) is nonnegative, we obtain kθS c k1 ≤ 3 · kθS k1 . In addition,
√
√
since S| = s∗ , kθS k1 ≤ s∗ · kθS k2 . Thus by (A.6) we have kθk2 ≤ s∗ · λ. Moreover, we
also have kθS k1 ≤ s∗ λ, which further implies that
kθk1 = kθS k1 + kθS c k1 ≤ 4 · kθS k1 ≤ 4s∗ λ.
Therefore, we conclude the proof.

A.2

Proof of Theorem 3.3

Proof. The proof of Theorem 3.3 is parallel to that of Theorem 3.2. Here the difference is
to handle the nuclear norm regularization, instead of the `1 -penalty. Since βb is the solution
of the optimization problem in (3.3), the first order optimality condition states that
b + λkβk
b ? ≤ L(µβ ∗ ) + λkµβ ∗ k? .
L(β)
30

(A.7)

To simplify the notation, we define Θ = βb − µ · β ∗ . Since L is quadratic,
b − L(µβ ∗ ) = h∇L(µβ ∗ ), Θi + 2kΘk2 ,
L(β)
fro

(A.8)

where ∇L takes values in Rd1 ×d2 . Then combining (A.7), (A.8), and Hölder’s inequality, we
have
b ? ≤ ∇L(µβ ∗ )
kΘk2fro ≤ −h∇L(µβ ∗ ), Θi + λkµβ ∗ k? − λkβk

op

b ?.
· kΘk? + λkµβ ∗ k? − λkβk
(A.9)

b ? in (A.9). Let U Λ∗ V > be the singular
In the following, we focus on the term kµβ ∗ k? − kβk
value decomposition of µβ ∗ , where U ∈ Rd1 ×d1 and V ∈ Rd2 ×d2 are orthogonal matrices, and
Λ∗ ∈ Rd1 ×d2 be formed by the singular values of µβ ∗ . Moreover, since rank(β ∗ ) = r∗ , Λ∗ can
be written in block form as


∗
Λ11 0
,
Λ∗ = 
0 0
where Λ∗11 ∈ Rr

∗ ×r ∗

(A.10)

is a diagonal matrix whose diagonal elements are the nonzero singular

values of µβ ∗ . We define Γ = U > ΘV , which can be written in block form as


Γ11 Γ12
,
Γ=
Γ21 Γ22
where Γ11 ∈ Rr

∗ ×r ∗

. In addition, we define matrices




Γ
Γ
0
0
 and Γ(2) =  11 12  .
Γ(1) = 
Γ21 0
0 Γ22

Then by (A.10) and triangle inequality of the nuclear norm, we have
b ? = kµβ ∗ + Θk? = kU (Λ∗ + Γ)V > k? = kΛ∗ + Γk?
kβk
≥ kΛ∗ + Γ(1) k? − kΓ(2) k? = kΛ∗ k? + kΓ(1) k? − kΓ(2) k? ,

(A.11)

where the last equality follows from the fact that Λ∗ + Γ(1) is block diagonal. Since kµβ ∗ k? =
kΛ∗ k? , by (A.11) we obtain
b ? ≤ kΓ(2) k? − kΓ(1) k? .
kµβ ∗ k? − kβk
31

(A.12)

In addition, triangle inequality implies that
kΘk? = kU ΓV > k? ≤ kΓ(1) k? + kΓ(2) k? .

(A.13)

Thus combining (A.11), (A.12), (A.13), we have
kΘk2fro ≤

∇L(µβ ∗ )


+
λ
· kΓ(2) k? +
op

∇L(µβ ∗ )


−
λ
· kΓ(1) k? ,
op

(A.14)

We utilize the following lemma to obtain an upper bound of k∇L(µβ ∗ )kop .
Lemma 2 (Upper bound of k∇L(µβ ∗ )kop ). Let loss function L : Rd1 ×d2 → R be defined in
p
p
(3.3) for the matrix setting. Setting κ = 2 n · log(d1 + d2 )/ (d1 + d2 )M , then it holds that
h
i
p
P k∇L(µβ ∗ )kop > 6 (d1 + d2 )/n ≤ (d1 + d2 )−2 .
Proof. See B.2 for a detailed proof.
By Lemma 2 and the choice of λ, we conclude that λ > 2 · k∇L(µβ ∗ )kop with probability
at least 1 − (d1 + d2 )−3 . Thus by (A.14) we have
kΘk2fro ≤ 3λ/2 · kΓ(2) k? − λ/2 · kΓ(1) k?

(A.15)

which implies that kΓ(1) k? ≤ 3 · kΓ(2) k? . Moreover, by the subadditivity of rank, we obtain








Γ
/2
Γ
Γ
/2
0
11
12
11
(2)
 + rank 
 = 2r∗ ,
rank(Γ ) ≤ rank 
0
0
Γ21 0
√
which implies that kΓ(2) k? ≤ 2r∗ · kΓ(2) kfro Then by (A.15) we obtain that kΘkfro ≤
√ √
3/ 2 · r∗ · λ. Finally, by triangle inequality for the nuclear norm,
√
kΘk? = kΓk? ≤ kΓ(1) k? + kΓ(2) k? ≤ 4 · kΓ(2) k? ≤ 4 2r∗ kΓ(2) kfro = 12r∗ λ.
Thus we conclude the proof of Theorem 3.3.

A.3

Proof of Theorem 4.2

c the solution of the optimization problem in (4.6). In addition, we
Proof. We denote by W
c − W ∗ kop .
let W ∗ = β ∗ β ∗ > . In the following, we establish an upper bound for kW
32

Since W ∗ is feasible for the optimization problem in (4.6), we have
c , Σi
e + λkW
c k1 ≥ hW ∗ , Σi
e + λkW ∗ k1 .
hW

(A.16)

We denote Σ∗ = E[Y · T (X)]. Note that β ∗ is the leading eigenvector of Σ∗ . Then (A.16) is
equivalent to
c − W ∗, Σ
e − Σ∗ i + λkW
c k1 − λkW ∗ k1 ≥ hΣ∗ , W ∗ − W
c i.
hW

(A.17)

The following Lemma in Vu et al. (2013) establishes an upper bound for the first term on
the left-hand side of (A.17).
Lemma 3. Let Ω ∈ Rd×d be a symmetric matrix and let λ1 ≥ λ2 ≥ . . . λd the eigenvalues
of Ω in descending order. For any ` ∈ [d − 1] such that λ` − λ`+1 > 0, let Π` ∈ Rd×d
be the projection matrix for the subspace spanned by the eigenvectors of Ω corresponding to
λ1 , . . . , λ` . Then for any Λ ∈ Rd×d satisfying 0  λ  Id and trace(Λ) = k, we have
(λ` − λ`+1 ) · kΠk − Λk2fro ≤ 2hΩ, Πk − Λi.
Note that W ∗ is the projection matrix for the subspace spanned by β ∗ . Applying Lemma
3 to Σ∗ with ` = 1, we have
c i ≥ C0 /2 · kW
c − W ∗ k2fro ,
hΣ∗ , W ∗ − W

(A.18)

where C0 > 0 is defined in (4.2). In addition, by Hölder’s inequality, we have
c − W ∗, Σ
e − Σ∗ i ≤ kΣ
e − Σ∗ k∞ · kW
c − W ∗ k1 .
hW

(A.19)

e − Σ∗ k∞ .
In what follows, we bound kΣ
e be defined in (4.5) and we define Σ∗ = E[Y · T (X)]. Under Assumption
Lemma 4. Let Σ
4.1, for any truncation level τ > 0 in (4.4), with probablity at least 1 − d−2 , we have
e − Σ∗ k∞ ≤ 9M · τ −3 + 2τ 3 · log d/n + 2
kΣ

p
5M · log d/n.

(A.20)

Proof. See §B.3 for a detailed proof.
By this lemma, if we set τ = (1.5M n/ log d)1/6 , then with probability at least 1 − d−1 ,
p
p
√
√
e − Σ∗ k∞ ≤ (2 5 + 2 6) · M log d/n ≤ 10 M log d/n.
kΣ
33

(A.21)

p
∗
e
Thus by setting λ = 10 M log d/n we have kΣ−Σ
k∞ ≤ λ with probability at least 1−d−2 .
Then combining (A.17), (A.18), and (A.19) we have


c − W ∗ k1 + kW
c k1 − kW ∗ k1 ≥ C0 /2 · kW
c − W ∗ k2 .
λ kW
fro

(A.22)

Note that W ∗ = β ∗ β ∗ > and that β ∗ is s∗ -sparse. We denote the support of W ∗ by J , which
is given by

J = (j, k) ∈ [d] × [d] : βj∗ · βk∗ 6= 0 .
Then by seperation of the `1 -norm, we have
cJ c k1 ,
c k1 = kW
cJ k1 + kW
cJ c k1 and kW
c − W ∗ k1 = kW
cJ − WJ∗ k1 + kW
kW
which implies that
c − W ∗ k1 + kW
c k1 − kW ∗ k1 = kW
cJ − WJ∗ k1 + kW
cJ k1 − kWJ∗ k1
kW
cJ − WJ∗ k1 ≤ 2s∗ kW
c − W ∗ kfro .
≤ 2kW

(A.23)

Here the last inequality in (A.23) follows from the fact that |J | = s∗ 2 . Combining (A.22)
and (A.23), we obtain
c − W ∗ kfro ≤ 4/C0 · s∗ λ.
kW
c , we have kβb−β ∗ k2 ≤
Since βb is the leading eigenvector of W

(A.24)
√

√
c −W ∗ kfro ≤ 4 2/C0 ·s∗ λ,
2kW

which concludes the proof.

A.4

Proof of Theorem 4.3

Proof. The proof is similar to that of Theorem 4.2. In the case of sparse MIM, we denote
c is the solution to the optimization problem in (4.7) and that B
b
W ∗ = B ∗ B ∗ > . Note that W
c . Then by Corollary 3.2 in Vu et al. (2013), we have
consists of the top-k eigenvectors of W
b − B ∗ Okfro ≤
inf kB

O∈Ok

√
c − W ∗ kfro .
2kW

(A.25)

c − W ∗ . Note that since B ∗ is orthonormal,
In what follows, we derive an upper bound for W
trace(W ∗ ) = k. Thus W ∗ is feasible for (4.7), which implies
c − W ∗, Σ
e − Σ∗ i + λkW
c k1 − λkW ∗ k1 ≥ hΣ∗ , W ∗ − W
c i.
hW
34

(A.26)

Here we define Σ∗ = E[Y · T (X)]. Note that W ∗ is the projection matrix for the subspace
spanned by the top-k leading eigenvectors of Σ∗ . By Lemma 3 with ` = k, we have
c i ≥ ρ0 /2 · kW
c − W ∗ k2 ,
hΣ∗ , W ∗ − W
fro
where ρ0 is the smallest eigenvalue of E[∇2 f (XB ∗ )]. Similar to the proof of Theorem 4.2,
by Hölder’s inequality and (A.26), we have
e − Σ∗ k∞ · kW
c − W ∗ k1 + λkW
c k1 − λkW ∗ k1 ≥ ρ0 /2 · kW
c − W ∗ k2 .
kΣ
fro

(A.27)

p
By Lemma ??, if we set λ = 10 M log d/n, with probability at least 1 − d−2 , we have
b − Σ∗ k∞ ≤ λ.
kΣ

(A.28)

Note that the support of W ∗ is

∗
J ⊆ (j, k) ∈ [d] × [d] : kBj·∗ k2 · kBk·
k2 6= 0 .
Since B ∗ is s∗ -row sparse, |J | ≤ s∗ 2 . Thus (A.23) also hold for the MIM. Combining (A.27),
(A.28), and (A.23), we obtain
c − W ∗ kfro ≤ 4/ρ0 · s∗ λ.
kW

(A.29)

Finally, combining (A.25) and (A.29), we conclude the proof.

B
B.1

Proof of Auxiliary Results
Proof of Lemma 1

Proof. By definition of the loss function L in (3.3), we have
n
n

 2X
2Xe e
e i ).
∇L(µβ ) = 2µβ −
Yi · S(Xi ) = E 2Yi · S(Xi ) −
Yei · S(X
n i=1
n i=1
∗

∗

By triangle inequality,




e
k∇L(µβ ∗ )k∞ ≤ E 2Y · S(X) − E 2Ye · S(X)

n

 2X
e
e i)
+ E 2Ye · S(X)
−
Yei · S(X
n i=1
∞

.
∞

(B.1)
35

e we
For any j ∈ [d], by the definition of the truncated response Ye and truncated score S,
obtain
n





o


E Ye · Sej (X) − E Y · Sj (X) ≤ E Ye · Sej (X) − Sj (X)
+ E (Ye − Y ) · Sj (X)




(B.2)
= E Ye · Sj (X) · 1{|Sj (X)| > τ } + E Y · Sj (X) · 1{|Y | > τ } .
{z
} |
{z
}
|
a1

a2

By Cauchy-Schwarz inequality, we have
 


a21 ≤ E Ye 2 Sj2 (X) · P |Sj (X)| ≥ τ
q


 
≤ E(Ye 4 ) · E Sj4 (X) · E Sj4 (X) · τ −4
= M 2 · τ −4 ,

(B.3)

where the second inequality follows from Chebyshev’s inequality. Similarly, for a2 we have



a22 ≤ E Y 2 Sj2 (X) · P |Y | ≥ τ
q


≤ E(Ye 4 ) · E S 4 (X) · E(Y 4 ) · τ −4
j

≤ M 2 · τ −4 .

(B.4)

Thus combining (B.2), (B.3), and (B.4), we conclude that




E Ye · Sej (X) − E Y · Sj (X) ≤ a1 + a2 ≤ 2M · τ −2
for all j ∈ [d]. Thus choosing τ = 2(M · n/ log d)1/4 , we have




E Ye · Sej (X) − E Y · Sj (X)

∞

≤ 1/2 ·

p
M · log d/n.

(B.5)

Furthermore, under Assumption 4.1, the variance of Ye · ·Sej (X) is bounded by
Var[Ye · Sej (X)] ≤ E[Ye 2 · Sej2 (X)] ≤ E[Y 2 · Sj2 (X)] ≤

q
E(Y 4 ) · E[Sj4 (X)] ≤ M.

Thus for the second term in (B.1), since |Ye · Sej (X)| ≤ τ 2 , by the Bernstein inequality in
Boucheron et al. (2013) (Theorem 2.10), for any j ∈ [d] and any t > 0, we have
r
 X

n


1
2M · t τ 2 · t
e
e
e
e
P
Yi · Sj (Xi ) − E Y · Sj (X) ≥
+
≤ exp(−t).
n i=1
n
3n

36

(B.6)

Taking union bound over j ∈ [t] in (B.6) yields
r


n
2


1Xe e
2M
·
t
τ
·
t
P
+
≤ exp(−t + log d).
Yi · Sj (Xi ) − E Ye · Sej (X)
≥
n i=1
n
3n
∞
Finally, we plug in τ = 2(M · n/ log d)1/4 and set t = 3 log d in (B.7) to obtain that
r
n
√


1Xe e
M · log d
Yi · Sj (Xi ) − E Ye · Sej (X)
≤ (4 + 6)
n i=1
n
∞

(B.7)

(B.8)

with probability at least 1 − d−2 . Finally, combining (B.1), (B.5), and (B.8), we conclude
the proof.

B.2

Proof of Lemma 2

Proof. For loss function L defined in (3.3) in the matrix setting, we have
n
n


2 X 
2 X 
ψ κ · Yi · S(Xi ) = 2E[Y · S(X)] −
ψ κ · Yi · S(Xi ) .
∇L(µβ ) = 2µβ −
κ · n i=1
κ · n i=1
∗

∗

(B.9)
Here the last equality follows from the generalized Stein’s identity. In the sequel, we apply
results in Minsker (2016) to bound k∇L(µβ ∗ )kop . To begin with, we first consider the
operator norm of E[Y 2 · S(X)S(X)> ] ∈ Rd! ×d2 and E[Y 2 · S(X)> S(X)] ∈ Rd2 ×d2 . For
notational simplicity, we denote by Sj,· (·) ∈ Rd2 S·,k (·) ∈ Rd1 the j-th row and k-the column
of the score function S(·), respectively. For any u ∈ S d1 −1 , by Cauchy-Schwarz inequality
we have
d2
q
X

 2
> 2
E[Y · u S(X)S(X) u] =
E [Y · S·,k (X) u] ≤ d2 · E(Y 4 ) · E [S·,1 (X)> u]4 ,
2

>

>

k=1

(B.10)
where we use the fact that the entries of S(X) are i.i.d. Since E[Sij (X)] = 0 and E[Sij4 (X)] ≤
M , by Cauchy-Schwarz inequality we obtain that
d X
d
X

> 4
E [S·,1 (X) u] =
E[Sj1 ,1 (X)2 · Sj22 ,1 (X)] · u2ji u2j2
j1 =1 j2 =1

≤

d X
d q
X

E[Sj41 ,1 (X)] · E[Sj42 ,1 (X)] · u2ji u2j2 ≤ M

j1 =1 j2 =1

d X
d
X
j1 =1 j2 =1

37

u2ji u2j2 = M.

(B.11)

Thus combining (B.10) and (B.11) we obtain that
E[Y 2 · u> S(X)S(X)> u] ≤ d2 · M,
which implies that kE[Y 2 ·S(X)S(X)> ]kop ≤ d2 ·M . Similarly, we obtain kE[Y 2 ·S(X)> S(X)]kop ≤
d1 · M. Thus by Corollary 3.1 in Minsker (2016), we have


n


1 X 
t
P
ψ κ · Yi · S(Xi ) − E[Y · S(X)]
≥√
κ · n i=1
n
op


√
≤ 2(d1 + d2 ) exp −κt n + κ2 (d1 + d2 )M/2



(B.12)

p
p
p
for any t > 0 and κ > 0. We set κ = 2 n · log(d1 + d2 )/ (d1 + d2 )M and t = (d1 + d2 )M · s
in (B.12), which implies that
r

n

1 X 
(d1 + d2 )M
P
ψ κ · Yi · S(Xi ) − E[Y · S(X)]
≥
·s
κ · n i=1
n
op

 p
≤ 2(d1 + d2 ) exp −2 log(d1 + d2 ) · s + 2 · log(d1 + d2 ) .


Now we set s = 3 ·

(B.13)

p
log(d1 + d2 ), which implies that the right-hand side of (B.13) is less

than


2(d1 + d2 ) exp −6 log(d1 + d2 ) + 2 · log(d1 + d2 )


≤ (d1 + d2 )2 · exp −4 · log(d1 + d2 ) = (d1 + d2 )−2 .
Therefore, combining (B.9) and (B.13) we conclude that
p
p
k∇L(µβ ∗ )kop ≤ 6 (d1 + d2 ) · M/n · log(d1 + d2 )
with probability at least 1 − (d1 + d2 )−2 , which concludes the proof.

B.3

Proof of Lemma 4

Proof. By triangle inequailty, we have
e − Σ∗ k∞ ≤ kΣ
e − EΣk
e ∞ + kEΣ
e − Σ∗ k ∞ .
kΣ

38

(B.14)

In the sequel, we bound the second term on the right-hand side of (B.14), which controls
the bias of truncation. For each j, k ∈ [d], we have




e jk − Σ∗ ≤ E Ye · Tejk (X) − E Y · Tjk (X)
EΣ
jk





≤ E Ye · Tejk (X) − Tjk (X)
+ E (Ye − Y ) · Tjk (X) .

(B.15)

For the first term in (B.15), note that Tejk (X) − Tjk (X) = Tjk (X) · 1{|Tjk (X)| ≥ τ 2 }. Then
by Cauchy-Schwarz inequality we have
n
h
io 2
h
i
2
e
e
e
E Y · Tjk (X) − Tjk (X)
= E Y · Tjk (X) · 1{|Tjk (X)| ≥ τ }
h
i 

2
≤ E Ye 2 · Tjk
(X) · P |Tjk (X)| ≥ τ 2 .
Furthermore, by Hölder’s inequality, we have
h
i h
i1/3  

2
E Ye 2 · Tjk
(X) ≤ E(Ye 6 )
· E |Tjk (X)|3

2/3


1/3   3

≤ E(Y 6 )
E |Tjk (X)|

2

(B.16)

2/3

. (B.17)

If j 6= k, by the definition of T (x) in (4.1), we have Tjk (x) = Sj (x) · Sk (x), ∀x ∈ Rd . Then
by Cauchy-Schwarz inequality, we have
 3


 q
E |Tjk
(X)| = E |Sj (X)|3 · |Sk (X)|3 ≤ E[Sj6 (X)] · E[Sk6 (X)] = E[Sj6 (X)].

(B.18)

In addition, if j = k, by (4.1), Tjj (x) = Sj2 (x) − s1 (xj ). Since (a + b)3 ≤ 4(a3 + b3 ) for any
a, b > 0, we have




E |Tjj3 (X)| ≤ 4E[Sj6 (X)] + 4E |s31 (Xj )| .

(B.19)

Moreover, by (B.16), (B.17), and the Markov’s inequality


 3

P |Tjk (X)| ≥ τ 2 ≤ E |Tjk
(X)| · τ −6 ,
we further have



E Ye · Tejk (X) − Tjk (X)

2


1/3   3

≤ E(Y 6 )
· E |Tjk (X)|

5/3

· τ −6 ≤ 32M 2 · τ −6 . (B.20)

Here the last inequality follows from combining Assumption 4.1, (B.18), and (B.19).
Similarly, for the second term in (B.15), by the Hölder’s inequality and the Markov’s
inequality we obtain
h
i2 
1/3   3
 2/3
e
E (Y − Y ) · Tjk (X) ≤ E(Y 6 )
· E |Tjk (X)|
· P(|Y | ≥ τ )

4/3   3
 2/3 −6
≤ E(Y 6 )
· E |Tjk (X)|
· τ ≤ 4M 2 · τ −6 .
39

(B.21)

Therefore, combining (B.15), (B.20), and (B.21) , we obtain
e − Σ∗ k∞ ≤ 9M · τ −3 .
kEΣ

(B.22)

e − EΣk
e ∞ using concentration
In what follows, we give a high-probability bound on kΣ
inequalities, which combined with B.22, concludes the proof.
For any j, k ∈ [d], note that |Ye · Tejk (X)| ≤ τ 3 . In addition, by assumption 4.1, its variance
is bounded by
h
i

 
1/3   3

2
e
e
Var Y · Tjk (X) ≤ E Y 2 · Tjk
(X) ≤ E(Y 6 )
· E |Tjk (X)|

2/3

≤ 2M.

Now we apply the Bernstein’s inequality (Boucheron et al., 2013) (Theorem 2.10) to {Yei ·
Tejk (Xi )}i∈[n] and obtain that
(
P

n

h
i
1Xe e
Yi · Tjk (Xi ) − E Ye · Tejk (X) ≥
n i=1

r

4M · t τ 3 · t
+
n
3n

)
≤ 2 exp(−t).

(B.23)

Taking a union bound over j, k ∈ [d] in (B.23), we obtain that
h

e − EΣk
e ∞≥
P kΣ

i
p
3
4M · t/n + τ · t/(3n) ≤ 2 exp(−t + 2 log d).

(B.24)

Choosing t = 5 log d in (B.24), we have
p
e − EΣk
e ∞ ≤ 2 5M log d/n + 2τ 3 · log d/n
kΣ

(B.25)

with probablity at least 1 − d−2 . Finally, combining (B.22) and (B.25), we conclude the
proof of Lemma 4.

References
Albert Ai, Alex Lapanowski, Yaniv Plan, and Roman Vershynin. One-bit compressed sensing
with non-gaussian measurements. Linear Algebra and its Applications, 441:222–239, 2014.
Pierre Alquier and Gérard Biau. Sparse single-index model. The Journal of Machine Learning
Research, 14(1):243–280, 2013.
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities: A
nonasymptotic theory of independence. Oxford university press, 2013.
40

Petros T Boufounos and Richard G Baraniuk. 1-bit compressive sensing. In Information
Sciences and Systems, 2008. CISS 2008. 42nd Annual Conference on, pages 16–21. IEEE,
2008.
Peter Bühlmann and Sara van de Geer. Statistics for high-dimensional data: methods, theory
and applications. Springer Science & Business Media, 2011.
T Tony Cai, Xiaodong Li, and Zongming Ma. Optimal rates of convergence for noisy sparse
phase retrieval via thresholded wirtinger flow. arXiv preprint arXiv:1506.03382, 2015.
Emmanuel J Candès, Thomas Strohmer, and Vladislav Voroninski. Phaselift: Exact and
stable signal recovery from magnitude measurements via convex programming. Communications on Pure and Applied Mathematics, 66(8):1241–1274, 2013.
Emmanuel J Candès, Yonina C Eldar, Thomas Strohmer, and Vladislav Voroninski. Phase
retrieval via matrix completion. SIAM Review, 57(2):225–251, 2015.
Olivier Catoni et al. Challenging the empirical mean and empirical variance: a deviation
study. Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 48(4):1148–1185,
2012.
Han Chen, Garvesh Raskutti, and Ming Yuan. Non-convex projected gradient descent for
generalized low-rank tensor regression. arXiv preprint arXiv:1611.10349, 2016.
Xin Chen, Changliang Zou, and R Dennis Cook. Coordinate-independent sparse sufficient
dimension reduction and variable selection. The Annals of Statistics, 38(6):3696–3723,
2010.
Mark A Davenport, Yaniv Plan, Ewout van den Berg, and Mary Wootters. 1-bit matrix
completion. Information and Inference, 3(3):189–223, 2014.
J. Fan, J. Lv, and L. Qi. Sparse high-dimensional models in economics. Annual review of
economics, 3(1):291–317, 2011.
Jianqing Fan, Weichen Wang, and Ziwei Zhu. Robust low-rank matrix recovery. arXiv
preprint arXiv:1603.08315, 2016.

41

Shmuel Friedland and Lek-Heng Lim. Nuclear norm of higher-order tensors. arXiv preprint
arXiv:1410.6072, 2014.
Larry Goldstein, Stanislav Minsker, and Xiaohan Wei. Structured signal recovery from nonlinear and heavy-tailed measurements. arXiv preprint arXiv:1609.01025, 2016.
David Gross, Felix Krahmer, and Richard Kueng. A partial derandomization of phaselift
using spherical designs. Journal of Fourier Analysis and Applications, 2015.
Aaron K Han. Non-parametric analysis of a generalized regression model: the maximum
rank correlation estimator. Journal of Econometrics, 35(2-3):303–316, 1987.
Joel L Horowitz. Semiparametric and nonparametric methods in econometrics, volume 12.
Springer, 2009.
Majid Janzamin, Hanie Sedghi, and Anima Anandkumar. Score function features for discriminative learning: Matrix and tensor framework. arXiv preprint arXiv:1412.2863, 2014.
Majid Janzamin, Hanie Sedghi, and Anima Anandkumar.

Beating the perils of non-

convexity: Guaranteed training of neural networks using tensor methods. arXiv preprint
arXiv:1506.08473, 2015.
B. Jiang and J. S. Liu. Variable selection for general index models via sliced inverse regression. The Annals of Statistics, 42(5):1751–1786, 2014.
Sham M Kakade, Varun Kanade, Ohad Shamir, and Adam Kalai. Efficient learning of
generalized linear and single index models with isotonic regression. In Advances in Neural
Information Processing Systems, pages 927–935, 2011.
Adam Tauman Kalai and Ravi Sastry. The isotron algorithm: High-dimensional isotonic
regression. In Conference on Learning Theory, 2009.
Robert Krauthgamer, Boaz Nadler, Dan Vilenchik, et al. Do semidefinite relaxations solve
sparse pca up to the information limit? The Annals of Statistics, 43(3):1300–1322, 2015.
Guillaume Lecué and Shahar Mendelson. Minimax rate of convergence and the performance
of empirical risk minimization in phase retrieval. Electron. J. Probab, 20(57):1–29, 2015.

42

Ker-Chau Li. Sliced inverse regression for dimension reduction. Journal of the American
Statistical Association, 86(414):316–327, 1991.
Ker-Chau Li. On principal Hessian directions for data visualization and dimension reduction:
Another application of Stein’s lemma. Journal of the American Statistical Association, 87
(420):1025–1039, 1992.
Ker-Chau Li and Naihua Duan. Regression analysis under link violation. The Annals of
Statistics, 17(3):1009–1052, 1989.
Xiaodong Li and Vladislav Voroninski. Sparse signal recovery from quadratic measurements
via convex programming. SIAM Journal on Mathematical Analysis, 45(5):3019–3033,
2013.
Q. Lin, Z. Zhao, and J. S. Liu. On consistency and sparsity for sliced inverse regression in
high dimensions. arXiv preprint arXiv:1507.03895, 2015.
Qian Lin, Xinran Li, Dongming Huang, and Jun S Liu. On the optimality of sliced inverse
regression in high dimensions. arXiv preprint arXiv:1701.06009, 2017.
Po-Ling Loh and Martin J Wainwright. Regularized m-estimators with nonconvexity: Statistical and algorithmic theory for local optima. Journal of Machine Learning Research,
16:559–616, 2015.
Geoffrey McLachlan and David Peel. Finite mixture models. John Wiley & Sons, 2004.
Shahar Mendelson. Learning without concentration. In Proceedings of The 27th Conference
on Learning Theory, pages 25–39, 2014.
Stanislav Minsker. Sub-gaussian estimators of the mean of a random matrix with heavytailed entries. arXiv preprint arXiv:1605.07129, 2016.
Cun Mu, Bo Huang, John Wright, and Donald Goldfarb. Square deal: Lower bounds and improved relaxations for tensor recovery. In Proceedings of The 31st International Conference
on Machine Learning, pages 73–81, 2014.

43

Sahand N. Negahban, Pradeep Ravikumar, Martin J. Wainwright, and Bin Yu. A unified
framework for high-dimensional analysis of M -estimators with decomposable regularizers.
Statistical Science, 27(4):538–557, 11 2012.
Roberto Imbuzeiro Oliveira.

The lower tail of random quadratic forms, with applica-

tions to ordinary least squares and restricted eigenvalue properties.

arXiv preprint

arXiv:1312.2903, 2013.
Yaniv Plan and Roman Vershynin. The generalized lasso with non-linear observations. IEEE
Transactions on information theory, 62(3):1528–1537, 2016.
Peter Radchenko. High dimensional single index models. Journal of Multivariate Analysis,
139:266–282, 2015.
Robert P Sherman. The limiting distribution of the maximum rank correlation estimator.
Econometrica: Journal of the Econometric Society, 61(1):123–137, 1993.
C. Stein. A bound for the error in the normal approximation to the distribution of a sum of
dependent random variables. In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory. The Regents of the
University of California, 1972.
Charles Stein, Persi Diaconis, Susan Holmes, Gesine Reinert, et al. Use of exchangeable pairs
in the analysis of simulations. In Stein’s Method. Institute of Mathematical Statistics, 2004.
Ju Sun, Qing Qu, and John Wright. A geometric analysis of phase retrieval. arXiv preprint
arXiv:1602.06664, 2016.
Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi. Lasso with non-linear measurements is equivalent to one with linear measurements. Advances in Neural Information
Processing Systems, 2015.
Roman Vershynin. Estimation in high dimensions: a geometric perspective. In Sampling
theory, a renaissance, pages 3–66. Springer, 2015.
Vincent Q Vu, Juhee Cho, Jing Lei, and Karl Rohe. Fantope projection and selection: A
near-optimal convex relaxation of sparse pca. In Advances in neural information processing
systems, pages 2670–2678, 2013.
44

Irène Waldspurger, Alexandre dAspremont, and Stéphane Mallat. Phase recovery, maxcut and complex semidefinite programming. Mathematical Programming, 149(1-2):47–81,
2015.
Tengyao Wang, Quentin Berthet, Richard J Samworth, et al. Statistical and computational
trade-offs in estimation of sparse principal components. The Annals of Statistics, 44(5):
1896–1930, 2016.
Zhaoran Wang, Huanran Lu, and Han Liu. Tighten after relax: Minimax-optimal sparse
pca in polynomial time. In Advances in neural information processing systems, pages
3383–3391, 2014.
Zhuoran Yang, Zhaoran Wang, Han Liu, Yonina C Eldar, and Tong Zhang. Sparse nonlinear
regression: Parameter estimation and asymptotic inference. International Conference on
Machine Learning, 2015.
Ming Yuan, Ding-Xuan Zhou, et al. Minimax optimal rates of estimation in high dimensional
additive models. The Annals of Statistics, 44(6):2564–2593, 2016.
Lixing Zhu, Baiqi Miao, and Heng Peng. On sliced inverse regression with high-dimensional
covariates. Journal of the American Statistical Association, 101(474):630–643, 2006.

45

