remarkRemark hypothesisHypothesis claimClaim factFact Smoothed Gradient Outer
ProductG. Yuan, M. Xu, S. Kpotufe and D. Hsu supplementary

# Efficient Estimation of the Central Mean Subspace via Smoothed Gradient
Outer Products

Gan Yuan Department of Statistics, Columbia University, New York, NY, 10027 (,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Mingyue Xu Department of
Computer Science, Purdue University, West Lafayette, IN, 47907 ().
<EMAIL_ADDRESS>Samory Kpotufe 22footnotemark: 2 Daniel Hsu Department of
Computer Science and Data Science Institute, Columbia University, New York,
NY, 10027<EMAIL_ADDRESS>
###### Abstract

We consider the problem of _sufficient dimension reduction_ (SDR) for multi-
index models. The estimators of the central mean subspace in prior works
either have slow (non-parametric) convergence rates, or rely on stringent
distributional conditions (e.g., the covariate distribution
$P_{\boldsymbol{X}}$ being elliptical symmetric). In this paper, we show that
a fast parametric convergence rate of form $C_{d}\cdot n^{-1/2}$ is achievable
via estimating the _expected smoothed gradient outer product_ , for a general
class of distribution $P_{\boldsymbol{X}}$ admitting Gaussian or heavier
distributions. When the link function is a polynomial with a degree of at most
$r$ and $P_{\boldsymbol{X}}$ is the standard Gaussian, we show that the
prefactor depends on the ambient dimension $d$ as $C_{d}\propto d^{r}$..

###### keywords:

central mean subspace, multi-index model, smoothed gradient outer product,
sufficient dimension reduction

62B05, 62G08

##  1 Introduction

We consider a regression model for data
$(\boldsymbol{X},Y)\in\mathbb{R}^{d}\times\mathbb{R}$ in which the regression
function
$g(\boldsymbol{x})=\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}=\boldsymbol{x}]$
has the form

$g(\boldsymbol{x})=f(\boldsymbol{U}^{\scriptscriptstyle{\mathsf{T}}}\boldsymbol{x})$

for some unknown matrix $\boldsymbol{U}\in\mathbb{R}^{d\times k}$ of rank $k$
and some unknown function $f\colon\mathbb{R}^{k}\to\mathbb{R}$. The goal of
this work is to provide computationally tractable methods for estimating the
column space of $\boldsymbol{U}$, denoted
$\operatorname{col}(\boldsymbol{U})$, from a finite sample of
$(\boldsymbol{X},Y)$.

The model described above is called a _multi-index model_ in the literature,
and $\operatorname{col}(\boldsymbol{U})$ is referred to as the _index space_.
The index space is also known as the _central mean subspace_ in the context of
sufficient dimension reduction [8, 10] under mild identifiability conditions
on $\operatorname{supp}(\boldsymbol{X})$111See Section 2.2 for a discussion of
these conditions.. An accurate estimate of the index space has many potential
applications, including the following:

  * •

The index space has, in the context of regression, an explanatory role
analogous to that of the leading principal components of the data. In
particular, it reveals which (linear combinations of) covariates are able to
best predict the response $Y$.

  * •

If $k\ll d$, then the achievable rates for estimating the regression function
$g$ may be greatly improved by first projecting the covariates
$\boldsymbol{x}$ to the index space before performing estimation. For
instance, if $g$ is only known to be Lipschitz, then we may hope to improve
the rate of estimation from $n^{-1/(d+2)}$ to $n^{-1/(k+2)}$.

The problem of estimating the index space has received much attention, but
important gaps remain in terms of what approximation error rates are
achievable by computationally tractable estimators. Results may be cataloged
into two main directions, one emphasizing fast convergence rates222Under
suitable subspace approximation measures, e.g., angles between spaces (see
Section 2.2). but requiring strict distributional conditions, the other
preferring general distributional settings but admitting significantly worse
convergence rates. For example, on one hand, so-called _inverse regression_
approaches—which estimate the index space via functionals of
$\operatorname*{\mathbb{E}}[\boldsymbol{X}\mid Y]$—often admit fast
convergence rates of the form $O(n^{-1/2})$ however at the cost of restrictive
conditions tying the marginal distribution $P_{\boldsymbol{X}}$ to
$\operatorname{col}(\boldsymbol{U})$; furthermore, outside of these conditions
on $P_{X}$, only a strict subspace of $\operatorname{col}(\boldsymbol{U})$ may
be recovered. On the other hand, _forward regression_ methods—which leverage
functionals of the regression function
$g:=\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$, e.g., its
derivatives—can ensure full recovery of the index space, importantly without
un-naturally forcing a relation between $P_{\boldsymbol{X}}$ and
$\operatorname{col}(\boldsymbol{U})$; however such milder distributional
conditions often come at the cost of significantly slower estimation rates,
often of a nonparametric form $n^{-1/\Omega(k)}\gg n^{-1/2}$.

Our main contribution is to elucidate the following middle-ground: we show
that fast rates of the form $C_{d}\cdot n^{-1/2}$ remain possible for a
reweighted variant of forward methods, for a prefactor $C_{d}$ depending on
data distribution. While our distributional conditions—knowledge of the
marginal distribution $P_{\boldsymbol{X}}$ and certain moment assumptions—are
stricter than for some forward regression methods, our rate guarantees are of
a parametric form even for nonparametric link $g$, similar to usual
convergence rates for inverse regression methods. Yet, unlike inverse
regression methods, we guarantee full recovery of the index space
$\operatorname{col}(\boldsymbol{U})$ under qualitatively milder distributional
conditions, where in particular, we impose no condition on the relation
between $P_{\boldsymbol{X}}$ and the unknown index space
$\operatorname{col}(\boldsymbol{U})$ as is common in that literature.

Our result relies on estimating outer products of a _smoothed_ version of
$\nabla f$, where the smoothing itself—induced by careful reweighting of the
training sample—leverages a first-order Stein’s Lemma to guarantee full
recovery of $\operatorname{col}(\boldsymbol{U})$. Similar to other existing
analyses, we use knowledge of $P_{\boldsymbol{X}}$ to properly reweighting of
the training sample.

To better situate our results, we expand next on the state-of-the-art on the
subject, followed by a more detailed overview of the results.

###  1.1 Prior Theoretical Works

The problem of index space estimation has received much attention over the
last few decades from both statistics and machine learning communities [21,
22, 9, 28, 15, 16, 11, 1, 19, 14, 5, 23]. For ease of comparison, we first
note that a prevalent approach in the literature, also adopted here, is to
proceed by defining and estimating a population object
$\boldsymbol{M}\in\mathbb{R}^{d\times d}$, whose column space estimates (part
of) $\operatorname{col}(\boldsymbol{U})$; the object $\boldsymbol{M}$ is
typically a functional of either the inverse or forward conditional
distribution $\boldsymbol{X}\mid Y$ or $Y\mid\boldsymbol{X}$.

##### Inverse regression methods

As the name suggests, the class of inverse regression methods leverage
functionals of the inverse relation $\boldsymbol{X}\mid Y$ that may reveal the
index space or a subspace thereof. Although these methods are initially
designed to estimate the _central subspace_ , they can also serve the purpose
of estimating the central mean subspace when the covariates $\boldsymbol{X}$
and the error $\varepsilon=Y-g(\boldsymbol{X})$ are independent, under which
case these two subspaces are identical. For instance, Sliced Inverse
Regression (SIR) by [21] estimates
$\boldsymbol{M}=\text{Cov}[\operatorname*{\mathbb{E}}[\boldsymbol{X}\mid Y]]$
and Sliced Average Variance Estimator (SAVE) by [9] estimates
$\boldsymbol{M}=\operatorname*{\mathbb{E}}[(\boldsymbol{I}_{d}-\text{Cov}[\boldsymbol{X}\mid
Y])^{2}]$. By carefully estimating means or covariances over _slices_ , i.e.,
partitions of the $Y$ space, parametric convergence rates of the form
$O(n^{-1/2})$ may be guaranteed, albeit under restrictive distributional
conditions. For instance, SIR imposes a so-called linear conditional mean
condition (where
$\operatorname*{\mathbb{E}}[\boldsymbol{X}\mid\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X}]$
is linear in
$\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X}$) in order to
ensure
$\operatorname{col}(\boldsymbol{M})\subset\operatorname{col}(\boldsymbol{U})$,
i.e., in order to guarantee _partial_ recovery of
$\operatorname{col}(\boldsymbol{U})$; if in addition we impose a so-called
random conditional projection condition (requiring for all
$v\in\mathbb{R}^{d}$ that
$\operatorname*{\mathbb{E}}[\boldsymbol{v}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X}\mid
Y]$ is not constant as a r.v.) then one can ensure
$\operatorname{col}(\boldsymbol{M})=\operatorname{col}(\boldsymbol{U})$. These
distributional conditions, by restricting the relation between
$P_{\boldsymbol{X}}$ and the index space $\operatorname{col}(\boldsymbol{U})$,
rule out some rather benign settings: for instance, for $k=1$, it is well
known that under Gaussian $P_{\boldsymbol{X}}$, these conditions cannot hold
for _even_ link functions such as $f(z)=z^{2}$ (also referred as phase
retrieval problem). Similar examples of failure modes are discussed for SAVE
in Appendix A. Additional discussions and examples for non-exhaustiveness
issues can be found in [14] and [19].

##### Forward regression methods

In contrast to inverse regression methods, the class of forward regression
methods leverage functionals of the regression function
$\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$. An important line of work
relies on Stein’s Lemma in constructing $\boldsymbol{M}$. The earliest work
along this line by [3] considered a single-index model ($k=1$) estimates the
index space direction, say $\boldsymbol{u}\in\mathbb{R}^{d}$, via
$\boldsymbol{M}=\operatorname*{\mathbb{E}}[Y\cdot\boldsymbol{X}]$: suppose
that $P_{\boldsymbol{X}}$ were standard Gaussian, and assuming an additive
noise model, then Stein’s Lemma implies that
$\operatorname*{\mathbb{E}}[Y\cdot\boldsymbol{X}]=\operatorname*{\mathbb{E}}[f^{\prime}(\boldsymbol{u}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})]\cdot\boldsymbol{u}$.
For the general case $k>1$, [29] proposed an estimator based on a second-order
Stein’s Lemma that is $\sqrt{n}$-consistent. Similar to our work, they assume
the setting on known $P_{X}$, and admit a broad class of distributions
including Gaussian and more heavy-tailed distributions. However, their
theoretical results relies on the non-singularity of the expected Hessian
$\operatorname*{\mathbb{E}}[\nabla^{2}g(\boldsymbol{X})]$, again ruling out
some apparently benign settings: for instance for $k=1$, and
$\operatorname*{\mathbb{E}}X=0$, simple link functions such as $f(z)=z^{3}$
are not admitted. In contrast, the approach analyzed here relies on a first-
order Stein’s Lemma (Lemma 3.6) and thus imposes no condition on the Hessian
or any other higher order derivatives of $g$.

Another class of forward regression methods are those that rely on functionals
of the gradient of $g$. In the simplest case of single index models ($k=1$),
we may take $\boldsymbol{M}=\operatorname*{\mathbb{E}}[\nabla
g(\boldsymbol{X})]=\operatorname*{\mathbb{E}}[f^{\prime}(\boldsymbol{u}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})]\cdot\boldsymbol{u}$,
for an index direction $\boldsymbol{u}$; thus, provided
$\operatorname*{\mathbb{E}}[f^{\prime}(\boldsymbol{u}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})]$
does not vanish (similar to [3] discussed above), [16] designs an estimator of
$\boldsymbol{u}$ with convergence rate $n^{-1/2}$. In the more general case
$k\geq 1$, one may consider a so-called _expected gradient outer product_
(EGOP), defined as $\boldsymbol{M}=\operatorname*{\mathbb{E}}[\nabla
g(\boldsymbol{X})\nabla
g(\boldsymbol{X})^{{\scriptscriptstyle{\mathsf{T}}}}]=\boldsymbol{U}\operatorname*{\mathbb{E}}[\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})^{{\scriptscriptstyle{\mathsf{T}}}}]\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}$.
This approach is quite appealing as it holds that, under quite mild conditions
(see Section 2.2), one can show that
$\operatorname{col}(\boldsymbol{M})=\operatorname{col}(\boldsymbol{U})$. Thus,
a number of estimators of the EGOP $\boldsymbol{M}$ have been proposed over
the years, e.g., [25, 28, 27, 24]. However, the convergence guarantees
available for these methods are, at best, of the nonparametric form
$n^{-1/\Omega(k)}$, without additional assumptions on $P_{\boldsymbol{X}}$.
For the restricted case of $k\leq 3$, [15, 11] shows that the parametric rate
$n^{-1/2}$ may be achieved under known Gaussian noise and in a fixed design
setting (under some conditions on the observed $\boldsymbol{X}$ data itself).
We have no such restriction on the index space dimension $k$, and operate in
the more general random design setting, however while assuming known
$P_{\boldsymbol{X}}$.

Finally, the popular MAVE approach of [28], while still a forward regression
method, differs considerably from all the above in that it does not explicitly
define a surrogate object $\boldsymbol{M}$ for $\boldsymbol{U}$; instead, it
works by optimizing a suitable objective over all possible projections
$\boldsymbol{U}$. MAVE is appealing in that, similar to EGOP, it can guarantee
full recovery of the index space under mild conditions. However, this also
comes at the cost of a nonparametric rate of the form $n^{-1/\Omega(k)}$.

###  1.2 Detailed Contributions

We introduce an object called the _Expected Smoothed Gradient Outer Product
(ESGOP)_ $\overline{\boldsymbol{M}}$, which can be viewed as a surrogate for
the EGOP. The ESGOP is based on a _smoothed_ version of the gradient $\nabla
g(\boldsymbol{X})$, as well as a subtle change of distribution (from
$P_{\boldsymbol{X}}$ to a Gaussian) that enables application of Stein’s
inequality.

We first show that $\overline{\boldsymbol{M}}$ is _exhaustive_ in that
$\operatorname{col}(\overline{\boldsymbol{M}})=\operatorname{col}(\boldsymbol{U})$
(Proposition 3.2). The surrogate $\overline{\boldsymbol{M}}$ can then be
estimated from a finite sample from $P_{\boldsymbol{X},Y}$ via a reweighting
scheme; in particular, the estimation of the _smoothed_ gradients may be
interpreted as local linear estimator with a data-dependent kernel (Remark
3.9).

Our main result states that, for a general class of distributions
$P_{\boldsymbol{X}}$ admitting Gaussians and heavier tail distributions, our
estimator converges at a parametric rate of the form $C_{d}\cdot n^{-1/2}$.
Finally, to elucidate the leading constant $C_{d}$, we consider the special
case of polynomial links $g$—which has received much recent interest—and show
that it behaves as $d^{O(r)}$ for polynomials of degree at most $r$ (Corollary
3.19), in line with recent lower-bounds of [12] for certain classes of
computationally tractable procedures.

###  1.3 Outline of the Paper

The remainder of the paper is organized as follows. In Section 2, we introduce
the basic setups and some notations. In Section 3, we propose a novel
estimator for central mean subspace based on the expected smoothed gradient
outer product, and provide upper-bounds for the estimation errors for general
link functions and polynomial link functions. In Section 4, we outline the
analysis of the main results presented in Section 3. In Section 5, we
demonstrate the choice of hyperparameters in the main algorithm through a
simulation study. Omitted proofs can be found in the Supplementary Material.

##  2 Preliminaries

###  2.1 Notations

Throughout, we use bold font letters (e.g. $\boldsymbol{x}$, $\boldsymbol{X}$)
to refer to vectors and matrices. Let $\langle\cdot,\cdot\rangle$ and
$\lVert\cdot\rVert$ denote Euclidean inner product and norm on
$\mathbb{R}^{d}$. For a probability measure $\nu$ on $\mathbb{R}^{d}$, we let
$\langle
f,g\rangle_{\nu}:=\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim\nu}[f(\boldsymbol{X})\cdot
g(\boldsymbol{X})]$ and $\lVert f\rVert_{\nu}:=\sqrt{\langle
f,f\rangle_{\nu}}$ denote the inner-product and norm on the function space
$L_{2}(\nu)$. The operator norm of a matrix $\boldsymbol{M}$ (with respect to
Euclidean norm) is denoted by
$\lVert\boldsymbol{M}\rVert_{\operatorname{op}}$. For a positive integer $d$,
we use $[d]$ to denote the set $\\{1,\ldots,d\\}$.

###  2.2 Central Mean Subspace

Let $(\boldsymbol{X},Y)\sim P_{\boldsymbol{X},Y}$, with $P_{\boldsymbol{X}}$
the marginal distribution assumed to be known, and let
$g(\boldsymbol{X}):=\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$ denote
the regression function.

######  Definition 2.1 (Mean dimension-reduction subspaces and central mean
subspace [7]).

A _mean dimension-reduction subspace_ for
$\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$ is a subspace
$W\subseteq\mathbb{R}^{d}$ such that
$\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$ is a function only of
$\boldsymbol{V}^{\scriptscriptstyle{\mathsf{T}}}\boldsymbol{X}$ for some
matrix $\boldsymbol{V}$ whose columns form a basis for $W$. Define

$\mathcal{U}:=\bigcap\\{W\subseteq\mathbb{R}^{d}:\text{$W$ is a mean
dimension-reduction subspace for
$\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$}\\}.$

We say $\mathcal{U}$ is the _central mean subspace (CMS)_ for
$\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$ if $\mathcal{U}$ is a mean
dimension-reduction subspace. (If $\mathcal{U}$ is not a mean dimension-
reduction subspace, then the CMS does not exist.)

The CMS is known to exist under mild conditions on
$\operatorname{supp}(\boldsymbol{X})$, as discussed by [7]; an example
sufficient condition is that $\operatorname{supp}(\boldsymbol{X})$ is open and
convex. Henceforth, we shall assume that $\operatorname{supp}(\boldsymbol{X})$
satisfies a sufficient condition that guarantees the existence of the CMS.
(Such an assumption will be explicit given later in 1.)

Let $k=\dim(\mathcal{U})$, and let the columns of
$\boldsymbol{U}:=[\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{k}]\in\mathbb{R}^{d\times
k}$ be an orthonormal basis for $\mathcal{U}$. Because $\mathcal{U}$ is the
CMS for $\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}]$, the regression
function has the form of a _$k$ -index model_, specifically

(1) $\displaystyle
g(\boldsymbol{X})=f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})$

for some function $f\colon\mathbb{R}^{k}\to\mathbb{R}$, which we call the
_link function_. Our goal is to estimate $\mathcal{U}$ (or equivalently,
$\operatorname{col}(\boldsymbol{U})$) from an i.i.d. sample from
$P_{\boldsymbol{X},Y}$. We evaluate our estimate using the following measure
of error.

######  Definition 2.2 (Distance with Optimal Rotation, Adapted from [6]).

The distance with optimal rotation between $\boldsymbol{U}$ and its estimates
$\hat{}\boldsymbol{U}$ is defined as:

$d(\hat{}\boldsymbol{U},\boldsymbol{U}):=\min_{\boldsymbol{R}\in\mathcal{O}^{k\times
k}}\lVert\hat{}\boldsymbol{U}\boldsymbol{R}-\boldsymbol{U}\rVert_{\operatorname{op}}\,,$

where $\mathcal{O}^{k\times k}$ is the set of all $k\times k$ orthogonal
matrices.

The distance with optimal rotation in Definition 2.2 turns out to be
equivalent to the so-called _$\mathbf{\sin}_{\Theta}$ distance_, another
popular notion of distance between subspaces given as follows. Let
$\boldsymbol{U}_{\perp}\in\mathbb{R}^{d\times(d-k)}$ denote a matrix whose
columns form an orthonormal basis for the orthogonal complement
$\mathcal{U}^{\perp}$ of $\mathcal{U}$. Then, the $\sin_{\Theta}$ distance
between $\boldsymbol{U}$ and $\widehat{\boldsymbol{U}}$ is defined as
$\sin_{\Theta}(\widehat{\boldsymbol{U}},\boldsymbol{U}):=\lVert\widehat{\boldsymbol{U}}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{U}_{\perp}\rVert_{\operatorname{op}}.$

It can be shown (Lemma 2.6 from [6]) that there exists $c,C>0$ such that

$c\cdot
d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq\sin_{\Theta}(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq
C\cdot d(\widehat{\boldsymbol{U}},\boldsymbol{U})_{\operatorname{op}}.$

The distance with optimal rotation is appealing in that it readily yields
control on regression error when data is first projected to
$\operatorname{col}(\widehat{\boldsymbol{U}})$ rather than $\mathcal{U}$. For
instance, suppose $f$ is $L$-Lipschitz (e.g, as in linear regression, or as in
polynomial regression with a bounded domain), and let
$\tilde{f}:=f\circ\boldsymbol{R}^{{\scriptscriptstyle{\mathsf{T}}}}$, where
$\boldsymbol{R}\in\mathcal{O}^{k\times k}$ is such that
$d(\widehat{\boldsymbol{U}},\boldsymbol{U})=\lVert\widehat{\boldsymbol{U}}\boldsymbol{R}-\boldsymbol{U}\rVert_{\operatorname{op}}$,
then:

$\lVert\tilde{f}(\widehat{\boldsymbol{U}}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})-f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})\rVert_{P_{\boldsymbol{X}}}:=\lVert
f((\widehat{\boldsymbol{U}}\boldsymbol{R})^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})-f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})\rVert_{P_{\boldsymbol{X}}}\leq
L\cdot\lVert\boldsymbol{X}\rVert_{P_{\boldsymbol{X}}}\cdot
d(\widehat{\boldsymbol{U}},{\boldsymbol{U}})\,.$

To estimate $\tilde{f}$, one can apply the classical Nadaraya-Watson estimator
on the projected sample of
$(\widehat{\boldsymbol{U}}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X},Y)$.
If $d(\widehat{\boldsymbol{U}},\boldsymbol{U})\lesssim n^{-1/(k+2)}$, the
regression error rate is the same as if $\boldsymbol{U}$ is known.

##  3 Main Results

We first introduce a surrogate object $\overline{\boldsymbol{M}}$ in Section
3.1 below, and show that it indeed is _exhaustive_ , i.e., it recovers all of
$\operatorname{col}({\bf U})$. Following this, we derive a separate surrogate
object $\widetilde{\boldsymbol{M}}$ in Section 3.3 which approximates
$\overline{\boldsymbol{M}}$ and turns out to be simpler to estimate and
retains the same exhaustiveness properties. Our main algorithm then estimates
$\widetilde{\boldsymbol{M}}$, and returns the column space $\hat{U}$ of this
estimate. Our main result of Theorem 3.13 which bounds the estimation error is
presented in Section 3.4, followed by an instantiation for polynomial link
functions in Section 3.5.

We will need the following few assumptions.

######  Assumption 1.

The joint distribution $P_{\boldsymbol{X},Y}$ satisfies the following basic
conditions.

  1. (1)

The link function $f$ is differentiable, with gradient $\nabla
f(\boldsymbol{z}),\forall\boldsymbol{z}\in\mathbb{R}^{k}$;

  2. (2)

The marginal $P_{\boldsymbol{X}}$ admits a Lebesgue density $p$ supported on
$\mathbb{R}^{d}$, i.e., $p>0$ everywhere.

The assumption of differentiability above implies that of the regression
function $g$, and is needed as our approach relies on gradient estimates. The
second assumption above will be required in designing the estimation
procedure, and also ensures the existence of $\mathcal{U}$ as previously
remarked.

###  3.1 Expected Smoothed Gradient Outer Product

We now present a smoothed version of the expected gradient outer product
$\boldsymbol{M}:=\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P_{\boldsymbol{X}}}[\nabla g(\boldsymbol{X})\nabla
g(\boldsymbol{X})^{{\scriptscriptstyle{\mathsf{T}}}}]$, which was discussed in
the introduction.

######  Definition 3.1.

The Expected Smoothed Gradient Outer Product (ESGOP)
$\overline{\boldsymbol{M}}$ is defined as follows. First, let $h>0$ be given,
and define the _smoothed gradient_ of $g$ at
$\boldsymbol{\theta}\in\mathbb{R}^{d}$ as:

$\displaystyle\boldsymbol{\beta}_{h}(\boldsymbol{\theta})$
$\displaystyle:=\boldsymbol{U}\cdot\overline{\nabla}_{h}f(\boldsymbol{\theta}),\quad\text{where
}\overline{\nabla}_{h}f(\boldsymbol{\theta}):=\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0}_{d},h^{2}\boldsymbol{I}_{d})}[\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}(\boldsymbol{Z}+\boldsymbol{\theta}))].$

Let $\sigma_{\theta}>0$. We then have:

(2)
$\overline{\boldsymbol{M}}:=\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})}\left[\boldsymbol{\beta}_{h}(\boldsymbol{\theta})\boldsymbol{\beta}_{h}(\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}\right]=\boldsymbol{U}\left(\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})}\left[\overline{\nabla}_{h}f(\boldsymbol{\theta})\overline{\nabla}_{h}f(\boldsymbol{\theta})^{\scriptscriptstyle{\mathsf{T}}}\right]\right)\boldsymbol{U}^{\scriptscriptstyle{\mathsf{T}}}.$

The ESGOP $\overline{\boldsymbol{M}}$ differs from the EGOP $\boldsymbol{M}$
in two important ways. First, the gradient $\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})$ is
replaced by a smoothed gradient $\overline{\nabla}_{h}f(\boldsymbol{\theta})$.
Second, the expectation in the definition of $\overline{\boldsymbol{M}}$ is
taken over a Gaussian distribution (i.e.,
$\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})$),
instead of the data distribution $P_{\boldsymbol{X}}$ as in the EGOP. As we
will see later in Section 3.3, such a design will allow us to invoke Stein’s
Lemma (Lemma 3.6) for Gaussians.

####  3.1.1 Exhaustiveness of the ESGOP $\overline{\boldsymbol{M}}$

From Eq. 2 it is clear that $\operatorname{col}(\overline{\boldsymbol{M}})$ is
contained in $\mathcal{U}:=\operatorname{col}(\boldsymbol{U})$. We argue here
that the two column spaces are in fact equal. From the same equation, this
boils down to arguing that

$\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})}\left[\overline{\nabla}_{h}f(\boldsymbol{\theta})\overline{\nabla}_{h}f(\boldsymbol{\theta})^{\scriptscriptstyle{\mathsf{T}}}\right]$

is full-rank. To gain some intuition, let’s first revisit a similar argument
for the EGOP
$\boldsymbol{M}=\boldsymbol{U}\left(\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P_{\boldsymbol{X}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})^{{\scriptscriptstyle{\mathsf{T}}}}\right)\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}.$

Suppose $\boldsymbol{M}$ were not exhaustive, then one can find a direction
$\boldsymbol{v}\in\mathbb{R}^{k}$ with $\lVert\boldsymbol{v}\rVert=1$, such
that $\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P_{\boldsymbol{X}}}[(\boldsymbol{v}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X}))^{2}]=0$.
Equivalently, $\boldsymbol{v}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})=0$ for
$P_{\boldsymbol{X}}$-almost every $\boldsymbol{x}$. Then,
$\boldsymbol{U}\boldsymbol{v}\in\mathbb{R}^{d}$ is a direction along which the
directional derivative of $g$ is zero almost surely. In other words,
$\mathcal{U}$ admits an irrelevant direction for the regression function $g$,
which contradicts the assumption that $\mathcal{U}$ is the CMS and hence
cannot be further reduced.

The exhaustiveness of $\overline{\boldsymbol{M}}$ is established in a similar
way as above, with the added difficulty that we now have to instead argue
about directional _smoothed_ derivatives
$\boldsymbol{v}^{\scriptscriptstyle{\mathsf{T}}}\overline{\nabla}_{h}f$ and
how they might similarly lead to a contradiction on the minimality of
$\mathcal{U}$. This is established by relating
$\boldsymbol{v}^{{\scriptscriptstyle{\mathsf{T}}}}\overline{\nabla}_{h}f$ to
$\boldsymbol{v}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla f$ by a Hermite
expansion of $f$. Note that, now the contradiction on minimality of
$\mathcal{U}$ can only be established $\mathcal{N}$-a.s., which fortunately
implies $P_{X}$-a.s. since $P_{X}$ has a density and is therefore
$\mathcal{N}$ dominated. Such arguments lead to the following proposition
whose proof is given in Section 4.1. The result requires an additional minor
condition on square integrability of $f$ so that Hermite expansions are well
defined.

######  Proposition 3.2 (Exhaustiveness of $\overline{\boldsymbol{M}}$).

Suppose that 1 holds, and that the link function $f$ satisfies
$\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0}_{k},h^{2}\boldsymbol{I}_{k})}[f(\boldsymbol{Z})^{2}]<\infty$.
Then, for $h,\sigma_{\theta}>0$, the ESGOP $\overline{\boldsymbol{M}}$ is
exhaustive, i.e., we have
$\operatorname{col}({\overline{\boldsymbol{M}}})=\mathcal{U}.$

###  3.2 Average Smoothed Gradient Outer Product

We now derive a first approximation to $\overline{\boldsymbol{M}}$ by
approximating the outer-expectation in Eq. 2 from a random sample of
$\boldsymbol{\theta}$ locations.

######  Definition 3.3.

Let $m,h>0$, and let
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{m}\overset{i.i.d.}{\sim}\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})$.
The Average Smoothed Gradient Outer Product (ASGOP) is obtained as:

(3)
$\displaystyle\widetilde{\boldsymbol{M}}:=\frac{1}{m}\sum_{j=1}^{m}\boldsymbol{\beta}_{h}(\boldsymbol{\theta}_{j})\boldsymbol{\beta}_{h}(\boldsymbol{\theta}_{j})^{{\scriptscriptstyle{\mathsf{T}}}}=\boldsymbol{U}\left(\frac{1}{m}\sum_{j=1}^{m}\left[\overline{\nabla}_{h}f(\boldsymbol{\theta}_{j})\overline{\nabla}_{h}f(\boldsymbol{\theta}_{j})^{\scriptscriptstyle{\mathsf{T}}}\right]\right)\boldsymbol{U}^{\scriptscriptstyle{\mathsf{T}}}.$

We will require that $m$ is sufficiently large so that
$\widetilde{\boldsymbol{M}}$ inherits the exhaustiveness of
$\overline{\boldsymbol{M}}$. This is established in the following proposition
via a concentration argument on eigenvalues (see Section 4.1). We require the
following moment condition for matrix concentration.

######  Assumption 2.

The smooth gradient $\overline{\nabla}_{h}f$ satisfies:

$\mu_{\nabla}:=\left(\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})}\left[\left\lVert\overline{\nabla}_{h}f\left({\boldsymbol{\theta}}\right)\right\rVert^{4}\right]\right)^{1/4}<\infty.$

Examples of link function $f$ for which the condition $\mu_{\nabla}<\infty$ is
satisfied include Lipschitz functions (in which case $\nabla f$ is bounded and
hence so is $\overline{\nabla}_{h}f$) and polynomials (see Lemma 4.6 (ii)).

######  Corollary 3.4 (Exhaustiveness of $\widetilde{\boldsymbol{M}}$).

Suppose that 2 and the condition of Proposition 3.2 holds, and Let
$0<\delta<1/2$, if
$m\geq{8\mu_{\nabla}^{4}}{\lambda_{k}(\overline{\boldsymbol{M}})^{-2}}\delta^{-1}(\log\left({4d}/{\delta}\right))^{2}$,
then with probability at least $1-\delta$, the ASGOP
$\widetilde{\boldsymbol{M}}$ satisfies
$\operatorname{col}(\widetilde{\boldsymbol{M}})=\mathcal{U}.$

######  Remark 3.5.

The above corollary is important in that we can now focus on estimating the
simpler object $\widetilde{\boldsymbol{M}}$, which turns out to be exhaustive
given just $m=O(1)$ locations $\boldsymbol{\theta}_{j}$’s. Such intuition
leads to the estimator of the next section, which our main theorem relies upon
(Theorem 3.13). Note that, estimating $\overline{\boldsymbol{M}}$ itself
appears to require more resources, namely, $m=O(n)$ (see Remark 3.15).

###  3.3 Estimation of CMS

We now present a simple estimator $\widehat{}\boldsymbol{M}$ of the
$\widetilde{\boldsymbol{M}}$ from i.i.d. samples
$\mathcal{D}:=\\{\boldsymbol{X}_{i},Y_{i}\\}_{i=1}^{n}\sim
P_{\boldsymbol{X},Y}$; let $\widehat{}\boldsymbol{U}$ denote the top $k$
eigenvectors of $\widehat{}\boldsymbol{M}$, then
$\operatorname{col}(\widehat{\boldsymbol{U}})$ estimates $\mathcal{U}$. We
assume access to the density $p$ of $P_{\boldsymbol{X}}$.

To this end, we first propose an unbiased estimator of the smoothed gradient
$\boldsymbol{\beta}_{h}(\boldsymbol{\theta})$ of the regression function $g$,
by relying on the following Stein’s Lemma:

######  Lemma 3.6 (Stein’s Lemma, [4]).

Let
$\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta},h^{2}\boldsymbol{I}_{d})$
and suppose $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is differentiable with
$\operatorname*{\mathbb{E}}\lVert\nabla g(\boldsymbol{Z})\rVert<\infty$. We
then have

$\displaystyle\operatorname*{\mathbb{E}}[g(\boldsymbol{Z})(\boldsymbol{Z}-\boldsymbol{\theta})]=h^{2}\operatorname*{\mathbb{E}}[\nabla
g(\boldsymbol{Z})].$

Algorithm 1 Estimation of the CMS $\mathcal{U}$

1:Input: $h>0,\sigma_{\theta}>0,m\in\mathbb{Z}$, dataset
$\mathcal{D}\subset\mathbb{R}^{d}\times\mathbb{R}$

2:Randomly sample $\\{\boldsymbol{\theta}_{j}\\}_{j=1}^{m}$ from
$\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})$

3:Split $\mathcal{D}$ into $\mathcal{D}_{1}$, $\mathcal{D}_{2}$, …, and
$\mathcal{D}_{m}$ with equal sizes

4:Split each $\mathcal{D}_{j}$ into $\mathcal{D}_{j,1}$ and
$\mathcal{D}_{j,2}$ with equal sizes

5:$\forall j\in[m]$, let
$\hat{\boldsymbol{\beta}}_{j,1}\leftarrow\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j},\mathcal{D}_{j,1})$,
and
$\hat{\boldsymbol{\beta}}_{j,2}\leftarrow\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j},\mathcal{D}_{j,2})$
$\triangleright$ As defined in Equation 5.

6:$\widehat{}\boldsymbol{M}\leftarrow\frac{1}{2m}\sum_{j=1}^{m}(\hat{\boldsymbol{\beta}}_{j,1}\hat{\boldsymbol{\beta}}^{{\scriptscriptstyle{\mathsf{T}}}}_{j,2}+\hat{\boldsymbol{\beta}}_{j,2}\hat{\boldsymbol{\beta}}^{{\scriptscriptstyle{\mathsf{T}}}}_{j,1})$

7:Return: $\widehat{\boldsymbol{U}}\in\mathbb{R}^{d\times k}$, the top-$k$
eigenvectors of $\widehat{}\boldsymbol{M}$

For intuition, notice that we have
$\boldsymbol{\beta}_{h}(\boldsymbol{\theta})=\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},h^{2}\boldsymbol{I}_{d})}\nabla
g(\boldsymbol{Z}+\boldsymbol{\theta})=\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta},h^{2}\boldsymbol{I}_{d})}\nabla
g(\boldsymbol{Z})$. Thus, suppose for a moment we were to sample
$\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\theta},h^{2}\boldsymbol{I}_{d})$,
then we would have

$\operatorname*{\mathbb{E}}Y\cdot(\boldsymbol{X}-\boldsymbol{\theta})=\operatorname*{\mathbb{E}}g(\boldsymbol{X})\cdot(\boldsymbol{X}-\boldsymbol{\theta})=h^{2}\boldsymbol{\beta}_{h}(\boldsymbol{\theta}),$

that is $h^{-2}Y\cdot(\boldsymbol{X}-\boldsymbol{\theta})$ would be unbiased
for $\beta_{h}(\boldsymbol{\theta})$. Our main insight therefore is to proceed
by importance-weighting, that is, reweighting the estimator as
$\rho_{h}(\boldsymbol{X};\boldsymbol{\theta})\cdot
Y\cdot(\boldsymbol{X}-\boldsymbol{\theta})$ so as to recover expectation under
$\mathcal{N}(\boldsymbol{\theta},h^{2}\boldsymbol{I}_{d})$. In particular,
recalling $p$ the density of $P_{\boldsymbol{X}}$, we let

(4)
$\displaystyle\rho_{h}(\boldsymbol{x};{\boldsymbol{\theta}}):=\frac{\varphi_{h}(\boldsymbol{x}-\boldsymbol{\theta})}{p(\boldsymbol{x})}\,,\quad\text{where
}\varphi_{h}(\cdot)\text{ denotes the density of
$\mathcal{N}(\boldsymbol{0}_{d},h^{2}\boldsymbol{I}_{d})$ }.$

##### Estimating $\boldsymbol{\beta}_{h}(\boldsymbol{\theta}_{j})$’s

Consider the $m$ values $\\{\boldsymbol{\theta}_{j}\\}_{j\in[m]}$ drawn from
$\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})$, upon
which $\widetilde{\boldsymbol{M}}$ is defined (Definition 3.3). In order to
evaluate $\boldsymbol{\beta}_{h}(\boldsymbol{\theta}_{j})$’s, assume (without
loss of generality) that the sample size $n$ is divisible by $2m$, and
partition the dataset $\mathcal{D}$ into independent
$\\{\mathcal{D}_{j,\ell}\\}_{j\in[m],\ell=1,2}$ of equal size. We then define
the following estimators:

(5)
$\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j};\mathcal{D}_{j,\ell}):=\frac{h^{-2}}{|\mathcal{D}_{j,\ell}|}\sum_{(\boldsymbol{X},Y)\in\mathcal{D}_{j,\ell}}\rho_{h}(\boldsymbol{X};\boldsymbol{\theta}_{j})\cdot
Y\cdot(\boldsymbol{X}-\boldsymbol{\theta}_{j}).$

Thus, for each $\boldsymbol{\theta}_{j}$, we have two independent estimates
$\hat{}\boldsymbol{\beta}_{j,\ell}:=\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j};\mathcal{D}_{j,\ell})$,
$\ell=1,2$. Each such estimate is unbiased, following the above intuition, as
formalized below.

######  Proposition 3.7.

We have
$\operatorname*{\mathbb{E}}\hat{\boldsymbol{\beta}}_{j,\ell}={\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j})$
for each $j\in[m]$ and $\ell=1,2.$

######  Remark 3.8.

Note that $P_{\boldsymbol{X}}$ being fully supported is a necessary condition
for importance reweighting to recover exactly the expectation under Gaussian
measure by sampling from $P_{\boldsymbol{X}}$. If $P_{\boldsymbol{X}}$ has a
bounded support, one can still approximate the expectation by choosing the
variance parameter $h$ small, see Appendix C for a detailed discussion.

Notice that, by the independence of the two estimates of
$\boldsymbol{\beta}_{h}(\boldsymbol{\theta}_{j})$, we also immediately have
that
$\hat{\boldsymbol{\beta}}_{j,1}\hat{\boldsymbol{\beta}}_{j,2}^{{\scriptscriptstyle{\mathsf{T}}}}={\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j}){\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j})^{{\scriptscriptstyle{\mathsf{T}}}}$.
This leads to the estimator of $\widetilde{\boldsymbol{M}}$ given below:

$\displaystyle\widehat{}\boldsymbol{M}:=\frac{1}{2m}\sum_{j=1}^{m}\left[\hat{\boldsymbol{\beta}}_{j,1}\hat{\boldsymbol{\beta}}_{j,2}^{{\scriptscriptstyle{\mathsf{T}}}}+\hat{\boldsymbol{\beta}}_{j,2}\hat{\boldsymbol{\beta}}_{j,1}^{{\scriptscriptstyle{\mathsf{T}}}}\right].$

The form of $\widehat{}\boldsymbol{M}$ ensures that it is symmetric. The full
procedure is detailed in Algorithm 1.

######  Remark 3.9 (Connection to Local Linear Regression).

Finally, we remark that the estimator
$\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta})$ is tightly related to a
_local linear regression estimator_ as they can be shown to match
asymptotically (see Appendix B).

###  3.4 Upper-bound on Estimation Error

We start with some further assumptions, some of which imply the various
conditions laid out so far in Proposition 3.2 and Corollary 3.4, as explained
below. These various conditions are required for concentration arguments on
$\hat{\boldsymbol{\beta}}_{h}$ and $\widehat{}\boldsymbol{M}$.

######  Assumption 3 (Tail Conditions for Noise).

The noise has a subgaussian tail, i.e., $\exists C,\sigma_{Y}^{2}>0$ such that
$\forall\boldsymbol{x}\in\mathbb{R}^{d}$,
$\operatorname*{\mathbb{P}}\\{Y-\operatorname*{\mathbb{E}}[Y\mid\boldsymbol{X}=\boldsymbol{x}]>t\mid\boldsymbol{X}=\boldsymbol{x}\\}\leq
C\cdot\exp(-t^{2}/(2\sigma_{Y}^{2}))$.

######  Assumption 4 (Moment Condition on Density Ratio).

The density ratio $\rho_{h}(\boldsymbol{X};\boldsymbol{0}_{d})$, as defined in
Eq. 4, satisfies:

$\displaystyle\mu_{\rho}:=\left(\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P_{\boldsymbol{X}}}\left[\rho_{h}(\boldsymbol{X};\boldsymbol{0}_{d})^{5}\right]\right)^{1/5}<\infty\,.$

######  Remark 3.10.

4 holds when $P_{\boldsymbol{X}}$ has a heavier tail than
$\mathcal{N}(\boldsymbol{0}_{d},h^{2}\boldsymbol{I}_{d})$, since
$\rho_{h}(\boldsymbol{X};\boldsymbol{0}_{d})$ is then bounded. Thus,
intuitively, lighter tail $P_{\boldsymbol{X}}$ are admitted by smaller choices
of $h$. This is illustrated in the following proposition.

######  Proposition 3.11 (Example of $\mu_{\rho}$ and choice of $h$).

Suppose
$P_{\boldsymbol{X}}=\mathcal{N}(\boldsymbol{0}_{d},\sigma^{2}\boldsymbol{I}_{d})$.
Then,

$\mu_{\rho}=\begin{cases}\left(\frac{5h^{8}}{\sigma^{8}}-\frac{4h^{10}}{\sigma^{10}}\right)^{-d/10}\,,&\text{
when }0<h<\sqrt{5}\sigma/2\\\ \infty\,,&\text{ when
}h\geq\sqrt{5}\sigma/2\end{cases}.$

Next, we require some additional moment conditions on the link function $f$.

######  Assumption 5 (Moment Condition on Link Function).

The link function $f$ satisfies:

$\displaystyle\mu_{f}:=\left(\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P_{\boldsymbol{X}}}[f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})^{6}]\right)^{1/6}<\infty\,.$

######  Remark 3.12.

We note that the condition of Proposition 3.2 (exhaustiveness of
$\overline{\boldsymbol{M}}$), i.e., $f$ is square integrable with respect to
$\mathcal{N}(\boldsymbol{0}_{k},h^{2}I_{k})$, holds under 4 and 5.

The following result establishes an upper-bound of order $O(n^{-1/2})$ on the
estimation error for the CMS $\mathcal{U}$, and holds for general link
functions $f$ satisfying the minor moment conditions discussed above. In
particular, this main result allows nonparametric, i.e., highly complex link
and regression functions $f$ and $g$, and yet admits fast convergence to
$\mathcal{U}$.

######  Theorem 3.13 (Main Result).

Let $\mathcal{D}$ be a random sample of size $n$ from $P_{\boldsymbol{X},Y}$,
$h,\delta>0$,
${16\mu_{\nabla}^{4}}{\lambda_{k}(\overline{\boldsymbol{M}})^{-2}}\delta^{-1}(\log\left({8d}/{\delta}\right))^{2}\leq
m\leq n/(2d)$ and $0<\sigma_{\theta}<{h}/\sqrt{20}$. Let
$\widehat{\boldsymbol{U}}$ be the output from Algorithm 1 with inputs
$(h,\sigma_{\theta},m,\mathcal{D})$. Under 1-5, there exists some absolute
constant $C_{1}>0$, such that with probability at least $1-\delta$,

(6) $\displaystyle d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq
C_{1}\Bigg{[}{C_{h,\sigma_{\theta},1}\cdot{d(\mu_{f}+\sigma_{Y})^{2}}}\cdot\frac{\sqrt{m}}{n}+{C_{h,\sigma_{\theta},2}\cdot{d^{1/2}(\mu_{f}+\sigma_{Y})}}\cdot\frac{1}{\sqrt{n}}\
\Bigg{]}\cdot\frac{\log({8d}/{{\delta}})}{\sqrt{\delta}}\,,$

where $C_{h,\sigma_{\theta},1}$ and $C_{h,\sigma_{\theta},2}$ are given by:

$\displaystyle C_{h,\sigma_{\theta},1}$
$\displaystyle:=(1-20\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot{\mu_{\rho}^{5/3}}\cdot(\lambda_{k}(\overline{\boldsymbol{M}}))^{-1},\text{
and }$ (7) $\displaystyle C_{h,\sigma_{\theta},2}$
$\displaystyle:=(1-10\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot\mu_{\nabla}\cdot\mu_{\rho}^{5/6}\cdot(\lambda_{k}(\overline{\boldsymbol{M}}))^{-1}.$

######  Remark 3.14 (Constants in Rate).

Since the term $\sqrt{m}/n\leq 1/\sqrt{n}$, Theorem 3.13 establishes that
$\mathcal{U}$ can be estimated at a parametric rate $C_{d}\cdot n^{-1/2}$, for
a constant

$C_{d}\leq C\cdot\max\\{C_{h,\sigma_{\theta},1}\cdot d,\
C_{h,\sigma_{\theta},2}\cdot d^{1/2}\\}.$

This prefactor $C_{d}$ depends indirectly on the condition number
$\lambda_{1}/\lambda_{k}$ of $\overline{\boldsymbol{M}}$ since the largest
eigenvalue $\lambda_{1}(\overline{\boldsymbol{M}})$ could be bounded by some
of the other moment parameters appearing in $C_{d}$; for instance it is easy
to see that
$\lambda_{1}(\overline{\boldsymbol{M}})\leq\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{0}_{d},\sigma_{\theta}^{2}\boldsymbol{I}_{d})}\left\lVert\overline{\nabla}_{h}f\left({\boldsymbol{\theta}}\right)\right\rVert^{2}\leq\mu_{\nabla}^{2}$.
For a more explicit dependence, see LABEL:rmk:largest_eigen relating to the
proof of LABEL:lem:beta-error.

Clearly, the various moment parameters, i.e.,
$\mu_{f},\mu_{\nabla},\mu_{\rho}$ appearing in $C_{d}$ depend tightly on
$P_{X,Y}$, and most importantly on the behavior of the link function $f$. In
Section 3.5 below we consider explicit polynomial conditions on $f$ along with
the example of Gaussian $P_{\boldsymbol{X}}$, and illustrate $C_{d}$ via
explicit bounds on moment conditions.

In the above result, while the lower-bound on $m$ was required for the
exhaustiveness of $\widetilde{\boldsymbol{M}}$ (Corollary 3.4), the upper-
bound $m\leq n/(2d)$ ensures that the datasets used to estimate smoothed
gradients $\overline{\nabla}_{h}f(\boldsymbol{\theta}_{j})$ are not empty.

######  Remark 3.15 (Estimating $\overline{\boldsymbol{M}}$).

As it appears, estimation of $\overline{\boldsymbol{M}}$ may require larger
values of $m$ than required in the above main theorem for estimation of
$\mathcal{U}$. For instance consider upper-bounding

$\lVert\widehat{}\boldsymbol{M}-\overline{\boldsymbol{M}}\rVert_{\operatorname{op}}\leq\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}+\lVert\widetilde{\boldsymbol{M}}-\overline{\boldsymbol{M}}\rVert_{\operatorname{op}}.$

The second term
$\lVert\widetilde{\boldsymbol{M}}-\overline{\boldsymbol{M}}\rVert_{\operatorname{op}}$
would necessarily introduce an $O(1/\sqrt{m})$ term which is of order
$O(n^{-1/2})$ only for large $m=\Omega(n)$.

######  Remark 3.16 (Dependence on $\delta$).

The dependence on the failure probability $\delta$ is of order
$\tilde{O}(\delta^{-1/2})$ which may appear large in light of usual high-
probability results. This is due to the fact that we make minimal
distributional assumptions, which in particular result in lower order moment
conditions on the random matrices
$\hat{\boldsymbol{\beta}}_{j,1}\hat{\boldsymbol{\beta}}_{j,2}^{{\scriptscriptstyle{\mathsf{T}}}}$,
preventing exponential concentration guarantees. Nonetheless, heavy tail
approaches such as _median-of-means_ (see e.g. [17]) may be employed to
improve such dependence on $\delta$.

###  3.5 Instantiation for Polynomial Links

In this section, we instantiate our main Theorem 3.13 under the assumptions of
a polynomial link $f$ and standard Gaussian $P_{X}$. Such assumptions have
been considered extensively in the recent literature on CMS estimation (see
e.g. [14, 5, 12]) and are formalized below.

######  Assumption 6.

We consider the following conditions on $P_{X,Y}$:

  1. (1)

The marginal density $P_{X}$ is
$\mathcal{N}(\boldsymbol{0}_{d},\boldsymbol{I}_{d})$ ;

  2. (2)

The link function $f$ is a polynomial of degree at most $r$. Furthermore, we
assume w.l.o.g that $f$ is normalized, i.e.,
$\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0}_{k},\boldsymbol{I}_{k})}[f(\boldsymbol{Z})^{2}]=1$
;

  3. (3)

3 holds, namely the noise is subgaussian with parameter $\sigma_{Y}^{2}>0$ .

6 above supersedes all previous 1-5 as will be evident in the proof of
Corollary 3.19 (see Section 4.4). In particular, as per Remark 3.14, we are
interested in understanding the constant $C_{d}$ in the rate of Theorem 3.13,
by explicitly bounding all the various moment parameters
$\mu_{f},\mu_{\rho},\mu_{\nabla}$ involved.

We adopt the following notion of minimum signal strength from [14], which will
serve to characterize $\lambda_{k}(\overline{\boldsymbol{M}})$.

######  Definition 3.17 (Minimum Signal Strength).

Let $\tau>0$ be define as:

$\tau:=\min_{\boldsymbol{\eta}\in\mathbb{R}^{k}}\
\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0}_{k},\boldsymbol{I}_{k})}\left[\left(\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{Z})\right)^{2}\right]\,.$

######  Remark 3.18.

The fact that $\tau\neq 0$ is guaranteed by the minimality of $\mathcal{U}$,
since otherwise, for some $\boldsymbol{\eta}\in\mathbb{R}^{k}$,
$\boldsymbol{U}\boldsymbol{\eta}\in\mathcal{U}$ would be an irrelevant
direction for the regression function $g$.

The following corollary establishes an upper-bound on subspace estimation
error, with an explicit dependence on $d,h$ and $\sigma_{\theta}$.

######  Corollary 3.19.

Suppose that 6 holds. Let $\widehat{\boldsymbol{U}}$ be the output from
Algorithm 1 with inputs $0<h<\sqrt{5}/2$, $0<\sigma_{\theta}<h/\sqrt{20}$ and
$m=n/(2d)$. Then, there exists $C_{2}>0$, independent of $n$, $d$, $\delta$,
$h$ and $\sigma_{\theta}$, such that with probability at least $1-\delta$,

$d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq C_{2}\cdot
A_{\sigma_{\theta}^{2}/h^{2},d,r}\cdot
B_{h,d,r}\cdot\frac{d}{\tau\sqrt{n}}\cdot\frac{\log({8d}/{{\delta}})}{\sqrt{\delta}}\,.$

where $A_{h,d,r},B_{\sigma_{\theta}^{2}/h^{2},d,r}$ are given as:

$\displaystyle A_{\sigma_{\theta^{2}/h^{2}},d,r}$
$\displaystyle:={\left(\sigma_{\theta}^{2}/h^{2}\right)^{-(r-1)}\cdot(1-20\sigma_{\theta}^{2}/h^{2})^{-d/2}}\,;$
$\displaystyle B_{h,d,r}$
$\displaystyle:=\max\\{h^{-2r+1},h^{r-1}\\}\cdot{(5h^{8}-4h^{10})^{-d/6}}\,.$

######  Remark 3.20 (Choice of $h,\sigma_{\theta}$).

In order to minimize the prefactor $C_{2}\cdot
A_{\sigma_{\theta}^{2}/h^{2},d,r}\cdot B_{h,d,r}$, we can minimize
$A_{\sigma_{\theta}^{2}/h^{2},d,r}$ (with respect to
$\sigma_{\theta}^{2}/h^{2}$) and $B_{h,d,r}$ (with respect to $h$) separately.
These terms are minimized at $\sigma^{2}_{\theta}/h^{2}=(r-1)/(20(r-1)+10d)$
and $h=1$. With such optimal choices, we have
$A_{\sigma_{\theta}^{2}/h^{2},d,r}\leq(20+30d/(r-1))^{r-1}$, and
$B_{h,d,r}=1$. Hence $d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq
C_{d}\cdot n^{-1/2}$ where $C_{d}\propto d^{r}$, ignoring the logarithmic term
and $\delta$.

######  Remark 3.21 (Bound Optimality).

The dependence on $d$ in the bound of Corollary 3.19 is hard to improve since
our estimator falls in the class of so-called _Correlational Statistical Query
(CSQ)_ learners, i.e., learners that rely solely on statistics of form
$\sum_{i}Y_{i}\cdot\Phi(\boldsymbol{X}_{i})$ for some function $\Phi$ (see [2,
18, 13] for an exact definition). Namely, such CSQ learners have been shown in
the recent work of [12] to require a sample complexity $n=\Omega(d^{r/2})$ to
achieve nontrivial error in subspace estimation. In comparison, our method
requires $n=O(d^{2r})$, which corresponds to a constant factor mismatch in the
exponent of $d$. It remains unclear how such a gap may be closed for CSQ
learners as the result of [12] also leaves this open. We note however that,
outside of CSQ procedures, much more benign sample size requirements are
possible, for instance $n=O(d)$ was shown for a recent method of [5].

##  4 Analysis Overview

In this section, we outline the proofs of results in Section 3. Some further
details can be found in LABEL:sec:omitted_proof.

###  4.1 Proof of Proposition 3.2 (Exhaustiveness of
$\overline{\boldsymbol{M}}$)

######  Proof 4.1.

We first consider the case $h=1$. Using the Hermite expansion of $f$:
$f\left(\boldsymbol{z}\right)=\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}w_{\boldsymbol{\alpha}}H_{\boldsymbol{\alpha}}\left(\boldsymbol{z}\right)$,
we have

$\displaystyle\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{d})}[\partial_{i}f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}(\boldsymbol{Z}+\boldsymbol{\theta}))]$
$\displaystyle=\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}:\alpha_{i}\geq
1}w_{\boldsymbol{\alpha}}\left(\prod_{j\neq
i}\frac{1}{\sqrt{\alpha_{j}!}}(\boldsymbol{U}_{j}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{\theta})^{\alpha_{j}}\right)\frac{\sqrt{\alpha_{i}}}{\sqrt{\left(\alpha_{i}-1\right)!}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}_{i}\boldsymbol{\theta})^{\alpha_{i}-1}$
$\displaystyle{=}\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}\sqrt{\alpha_{i}+1}\cdot
w_{\boldsymbol{\alpha}{(i)}}\prod_{j\in[k]}\left(({\boldsymbol{U}_{j}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{\theta})^{\alpha_{j}}}/{\sqrt{\alpha_{j}!}}\right),$
(8)
$\displaystyle=:\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}\lambda_{\boldsymbol{\alpha},i}S_{\boldsymbol{\alpha}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{\theta})\,,$

where
$\boldsymbol{\alpha}{(i)}:=(\alpha_{1},\ldots,\alpha_{i-1},\alpha_{i}+1,\alpha_{i+1},\ldots,\alpha_{k})\in\mathbb{R}^{k}$,
$\lambda_{\boldsymbol{\alpha},i}:=\sqrt{\alpha_{i}+1}\cdot
w_{\boldsymbol{\alpha}{(i)}}\in\mathbb{R}$ and
$S_{\boldsymbol{\alpha}}(\boldsymbol{z})=\prod_{j\in[k]}({z_{j}^{\alpha_{j}}}/{\sqrt{\alpha_{j}!}})$
is a scaled monomial basis indexed by $\boldsymbol{\alpha}$.

According to the definition,
$\boldsymbol{M}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}$,
where
$\boldsymbol{\Lambda}:=\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}}[\overline{\nabla}_{h}f(\boldsymbol{\theta})\overline{\nabla}_{h}f(\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}]\in\mathbb{R}^{k\times
k}$. To show
$\operatorname{col}({\boldsymbol{M}})=\operatorname{col}(\boldsymbol{U})$, we
only need to show that the matrix $\boldsymbol{\Lambda}$ is full-rank. Suppose
not, then there exists a unit vector $\boldsymbol{\eta}\in\mathbb{R}^{k}$ such
that
$\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{\Lambda}\boldsymbol{\eta}=0$,
which indicates
$\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\overline{\nabla}_{h}f(\boldsymbol{\theta})=0$,
for almost every $\boldsymbol{\theta}$ with respect to Gaussian measure.
Therefore,

(9)
$\displaystyle\forall\boldsymbol{\alpha}\in\mathbb{N}^{k},\quad\sum_{i\in[k]}\eta_{i}\lambda_{\boldsymbol{\alpha},i}=0\,.$

Next, we show that Eq. 9 implies that
$\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})=0$ for
almost every $\boldsymbol{x}\in\mathbb{R}^{d}$. For each $i\in[k]$,

$\displaystyle\partial_{i}f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}w_{\boldsymbol{\alpha}}\cdot{\partial_{i}H_{\boldsymbol{\alpha}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})}$
$\displaystyle=\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}:\alpha_{i}\geq
1}w_{\boldsymbol{\alpha}}\cdot(\sqrt{\alpha_{i}}H_{\alpha_{i}-1}(\boldsymbol{U}_{i}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x}))\prod_{j\neq
i}H_{\alpha_{j}}(\boldsymbol{U}_{j}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})$
$\displaystyle{=}\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}\sqrt{\alpha_{i}+1}\cdot
w_{\boldsymbol{\alpha}{(i)}}\cdot
H_{\boldsymbol{\alpha}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})$
(10)
$\displaystyle=\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}\lambda_{\boldsymbol{\alpha},i}\cdot
H_{\boldsymbol{\alpha}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})\,.$

Therefore, we have for almost every $\boldsymbol{x}\in\mathbb{R}^{d}$ with
respect to the Gaussian measures, and hence $P_{\boldsymbol{X}}$-a.s. (since
$P_{\boldsymbol{X}}$ is dominated by Gaussian),

$\displaystyle\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})\overset{\lx@cref{creftype~refnum}{eq:partial}}{=}\sum_{\boldsymbol{\alpha}\in\mathbb{N}^{k}}\left(\sum_{i\in[k]}\eta_{i}\lambda_{\boldsymbol{\alpha},i}\right)H_{\boldsymbol{\alpha}}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})\overset{\lx@cref{creftype~refnum}{eqn:zero_coeff}}{=}0\,.$

Then, it follows that

$(\boldsymbol{U}\boldsymbol{\eta})^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
g(\boldsymbol{x})=\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{U}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})=\boldsymbol{\eta}^{{\scriptscriptstyle{\mathsf{T}}}}\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{x})=0$

In other words, $\boldsymbol{U}\boldsymbol{\eta}\in\mathcal{U}$ is a direction
along which the directional derivative of $g$ is zero almost surely. In other
words, $\mathcal{U}$ admits an irrelevant direction for the regression
function $g$, which contradicts the assumption that $\mathcal{U}$ is the CMS
and hence cannot be further reduced.

For the case $h\neq 1$, one can let
${f}_{h}(\boldsymbol{z}):=f(h\boldsymbol{z})$, then
$\operatorname*{\mathbb{E}}_{{\boldsymbol{Z}}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{d})}[f_{h}(\boldsymbol{Z})^{2}]=1$,
and

$\overline{\nabla}_{h}f(\boldsymbol{\theta}):=\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},h^{2}\boldsymbol{I}_{d})}[\nabla
f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}(\boldsymbol{Z}+\boldsymbol{\theta}))]=h^{-1}\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{d})}[\nabla
f_{h}(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}(\boldsymbol{Z}+\boldsymbol{\theta})))]=h^{-1}\overline{\nabla}_{1}f_{h}(\boldsymbol{\theta}).\,$

Clearly,
$\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}}[\overline{\nabla}_{1}\tilde{f}(\boldsymbol{\theta})\overline{\nabla}_{1}\tilde{f}(\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}]$
is full-rank by the reasoning for the case $h=1$, and it follows that
$\overline{\boldsymbol{M}}=h^{-2}\boldsymbol{U}\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}}[\overline{\nabla}_{1}\tilde{f}(\boldsymbol{\theta})\overline{\nabla}_{1}\tilde{f}(\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}]\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}$
is exhaustive.

###  4.2 Proof of Corollary 3.4 (Exhaustiveness of
$\widetilde{\boldsymbol{M}}$)

The following lemma (see LABEL:sec:proof_ASGOP for the proof) establishes the
concentration of the ASGOP $\widetilde{\boldsymbol{M}}$ around ESGOP
$\overline{\boldsymbol{M}}$, and is the key ingredient for proving the
exhaustiveness of $\widetilde{\boldsymbol{M}}$.

######  Lemma 4.2.

Under the setting of Theorem 3.13, with probability $1-\delta$,

$\lVert\widetilde{\boldsymbol{M}}-\overline{\boldsymbol{M}}\rVert_{\operatorname{op}}\leq\left({\sqrt{2}\mu_{\nabla}^{2}}\cdot\frac{1}{\sqrt{m}}\right)\cdot\frac{1}{\sqrt{{\delta}}}\cdot\log\left(\frac{4d}{{\delta}}\right)\,,$

for all $0<\delta<1$.

Now we proceed to the proof of the Corollary 3.4.

######  Proof 4.3 (Proof of Corollary 3.4 (Exhaustiveness of
$\widetilde{\boldsymbol{M}}$)).

By Lemma 4.2, we have with probability at least $1-\delta$,

$\lVert\widetilde{\boldsymbol{M}}-\overline{\boldsymbol{M}}\rVert_{\operatorname{op}}\leq{\lambda_{k}(\overline{\boldsymbol{M}})}/{2}\,,$

when
$m\geq{8\mu_{\nabla}^{4}}{\lambda_{k}(\overline{\boldsymbol{M}})^{-2}}\delta^{-1}(\log\left({4d}/{\delta}\right))^{2}$.
Consequently,
$\lambda_{k}(\widetilde{\boldsymbol{M}})\geq\lambda_{k}(\overline{\boldsymbol{M}})/2>0$
by Weyl’s inequality. Therefore,
$\operatorname{col}(\widetilde{\boldsymbol{M}})$ is of rank $k$, and hence
equals to $\mathcal{U}$.

###  4.3 Proof of Theorem 3.13 (Main Theorem)

The following Lemma 4.4 (see LABEL:sec:proof_estimator for the proof)
establishes a probabilistic bound on
$\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}$.
It is a key step for the proof of the main theorem.

######  Lemma 4.4.

Under the setting of Theorem 3.13, for all $0<\delta<1$ and some absolute
constant $C>0$, with probability $1-\delta$,

(11)
$\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}\leq
C\left(\tilde{C}_{h,\sigma_{\theta},1}\cdot
d(\mu_{f}+\sigma_{Y})^{2}\cdot\frac{\sqrt{m}}{n}+\tilde{C}_{h,\sigma_{\theta},2}\cdot
d^{1/2}(\mu_{f}+\sigma_{Y})\cdot\frac{1}{\sqrt{n}}\right)\frac{\log({4d}/{{\delta}})}{\sqrt{{\delta}}}\,.$

where $\tilde{C}_{h,\sigma_{\theta},\ell}$, $\tilde{C}_{h,\sigma_{\theta},2}$
are as defined as:

$\displaystyle\tilde{C}_{h,\sigma_{\theta},1}:=(1-20\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot\mu_{\rho}^{5/3};\quad\tilde{C}_{h,\sigma_{\theta},2}:=(1-10\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot\mu_{\nabla}\cdot\mu_{\rho}^{5/6}.$

Now we proceed to the proof of the main theorem.

######  Proof 4.5 (Proof of Theorem 3.13 (Main Theorem)).

Let $\tilde{\delta}=\delta/2\in(0,1/2)$. We define two good events with
desired properties for the concentration of the subspace estimator: Let
$\xi_{1}$ denote the event where $\mathcal{U}$ can be fully recovered from
$\widetilde{\boldsymbol{M}}$, i.e.,
$\operatorname{col}(\widetilde{\boldsymbol{M}})=\mathcal{U}$. Since
$m\geq{8\mu_{\nabla}^{4}}{\lambda_{k}(\overline{\boldsymbol{M}})^{-2}}\tilde{\delta}^{-1}(\log({4d}/{\tilde{\delta}}))^{2}$,
we have the guarantee that $\operatorname*{\mathbb{P}}(\xi_{1})\geq
1-\tilde{\delta}$ by Corollary 3.4. Let $\xi_{2}$ denote the good event in
Lemma 4.4 with $\delta=\tilde{\delta}$. Next, we show that Eq. 6 holds on
$\xi_{1}\cap\xi_{2}$, whose probability is at least $1-\delta$. This follows
from a standard application of the _Davis-Kahan Theorem_ (Theorem V.3.6 from
[26]). In particular, the eigengap can be calculated as:

$\displaystyle\lambda_{k}(\widetilde{\boldsymbol{M}})-\lambda_{k+1}(\widehat{}\boldsymbol{M})$
$\displaystyle=\lambda_{k}(\widetilde{\boldsymbol{M}})-(\lambda_{k+1}(\widehat{}\boldsymbol{M})-\lambda_{k+1}(\widetilde{\boldsymbol{M}}))$
$\displaystyle\geq\lambda_{k}(\widetilde{\boldsymbol{M}})-\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}$
$\displaystyle\geq\lambda_{k}(\widetilde{\boldsymbol{M}})/2\,,$

for large $n$ that ensures
$\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}\leq\lambda_{k}(\widetilde{\boldsymbol{M}})/2$.
Therefore,

$d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq\frac{\sqrt{2}\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}}{\lambda_{k}(\widetilde{\boldsymbol{M}})-\lambda_{k+1}(\widehat{}\boldsymbol{M})}\leq\frac{2\sqrt{2}\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}}{\lambda_{k}(\widetilde{\boldsymbol{M}})}\leq\frac{4\sqrt{2}\lVert\widehat{}\boldsymbol{M}-\widetilde{\boldsymbol{M}}\rVert_{\operatorname{op}}}{\lambda_{k}(\overline{\boldsymbol{M}})}.$

By substituting Eq. 11 into the last inequality, we conclude the proof.

###  4.4 Proof of Corollary 3.19

To invoke the main Theorem 3.13 for the special case of polynomial links, we
rely on the following lemma (see LABEL:sec:proof:lemma_poly for proof) that
quantifies the various moment parameters and
$\lambda_{k}(\overline{\boldsymbol{M}})$ in terms of $d,r,\tau,h$ and
$\sigma_{\theta}$.

######  Lemma 4.6.

Under the conditions for Corollary 3.19, we have:

  1. (i)

$\mu_{\rho}=(5h^{8}-4h^{10})^{-d/10}\,;$

  2. (ii)

$\mu_{f}\leq 5^{r/2}\,;$

  3. (iii)

$\mu_{\nabla}\leq(r+1)^{(6k+r+2)/4}\cdot\max\\{h^{r-1},h^{-1}\\}\,;$

  4. (iv)

$\lambda_{k}\left(\overline{\boldsymbol{M}}\right)\geq r^{-r-4k}\cdot
2^{-3(r-1)-k}\cdot((r-1)!)^{-1}\cdot\left(\sigma_{\theta}^{2}/h^{2}\right)^{r-1}\cdot\min\\{h^{2(r-1)},1\\}\cdot\tau\,.$

Next, we proceed to prove Corollary 3.19, which merely substitutes the bounds
from Lemma 4.6 into Theorem 3.13.

######  Proof 4.7 (Proof of Corollary 3.19).

By Theorem 3.13 and Lemma 4.6, we have:

$\displaystyle C_{h,\sigma_{\theta},1}$ $\displaystyle\leq
C_{k,r}\cdot\max\\{h^{-2r+2},1\\}{(5h^{8}-4h^{10})^{-d/6}}\left(\sigma_{\theta}^{2}/h^{2}\right)^{-(r-1)}(1-20\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot\tau^{-1}\,;$
$\displaystyle C_{h,\sigma_{\theta},2}$ $\displaystyle\leq
C_{k,r}\cdot\max\\{h^{-2r+1},h^{r-1}\\}{(5h^{8}-4h^{10})^{-d/12}}\left(\sigma_{\theta}^{2}/h^{2}\right)^{-(r-1)}(1-10\sigma_{\theta}^{2}/h^{2})^{-d/2}\cdot\tau^{-1}\,,$

where $C_{k,r}$ is a constant depends only on $k$ and $r$. Thus, we have
$\max\\{d\cdot C_{h,\sigma_{\theta},1},d^{1/2}\cdot
C_{h,\sigma_{\theta},2}\\}\leq C_{k,r}\cdot
A_{\sigma_{\theta}^{2}/h^{2},d,r}\cdot B_{h,d,r}\cdot d$. Finally, we invoke
Theorem 3.13 and get:

$\displaystyle d(\widehat{\boldsymbol{U}},\boldsymbol{U})\leq C_{2}\cdot
A_{\sigma_{\theta}^{2}/h^{2},d,r}\cdot
B_{h,d,r}\cdot\frac{d}{\tau\sqrt{n}}\cdot\frac{1}{\tau}\cdot\frac{\log({8d}/{{\delta}})}{\sqrt{\delta}}\,,$

where $C_{2}:=C_{1}\cdot
C_{k,r}\cdot\max\\{\mu_{f}+\sigma_{Y},(\mu_{f}+\sigma_{Y})^{2}\\}$.

##  5 Simulations

In this section, we demonstrate how to choose the hyperparameters $h$ (the
“radius” of smoothing) and $m$ (the number of $\boldsymbol{\theta}_{j}$
sampled in ASGOP) properly through a simple example with polynomial link
function and standard Gaussian design. The related theoretical discussions can
be found in Remark 3.5, 3.10, 3.15 and 3.20.

##### Experiment Settings

Throughout the experiments, we consider a multi-index model Eq. 1 with
$d=10,k=3$, and
$\boldsymbol{U}=[\boldsymbol{e}_{1},\ldots,\boldsymbol{e}_{k}]\in\mathbb{R}^{d\times
k},$ where $\boldsymbol{e}_{i}$ denotes the $i$-th column of the $d\times d$
identity matrix. That is, only the first three coordinates of $\boldsymbol{X}$
are relevant directions for $Y$. The link function
$f:\mathbb{R}^{3}\to\mathbb{R}$ is a polynomial:
$f(\boldsymbol{z})=z_{1}^{2}+z_{2}z_{3}$.

##### Choices of $h$ under Different $P_{\boldsymbol{X}}$

Recall from Remark 3.10 that our method work for any $P_{\boldsymbol{X}}$ with
heavier tails than Gaussian. The choice of $h$ is critical when
$P_{\boldsymbol{X}}$ is a Gaussian distribution (see Proposition 3.11). We
verify the statement by applying Algorithm 1 with different choices of $h>0$
and $\sigma_{\theta}=h/\sqrt{20+10d}$ (see Remark 3.20 for the reasoning of
such choice) to datasets generated from $P_{\boldsymbol{X}}$ being either
standard Gaussian or standard Cauchy. Under appropriate choices of $h$, our
method works reasonably well for both cases.

Figure 1: The subspace estimation error
$d(\widehat{\boldsymbol{U}},\boldsymbol{U})$ v.s. sampling budget $n$. Here,
we fixed the number of partitions $m=15$ and
$\sigma_{\theta}=h/\sqrt{20+10d}$. We replicate 10 times for each pair of
$(n,h)$ and plot the mean (the dots) and the standard error (the error bars)
of the estimation error. When $P_{\boldsymbol{X}}=$ standard Gaussian (left),
the performance of Algorithm 1 is quite sensitive to the choice of $h$. The
optimal choice of $h$ is around the data variance 1. When $h$ gets smaller
(e.g., $h=0.5$) or larger (e.g., $h=1.2,1.5$), we observe larger errors under
the same budget level. This actually coincides approximately with the
minimizer of $\mu_{\rho}$ in terms of $h$ (c.f. Proposition 3.11). When
$P_{\boldsymbol{X}}=$ standard Cauchy, the method is more robust to the choice
of $h$.

##### Choice of the Number of Partitions $m$ under Gaussian design

The number of partition $m$ affects the subspace estimation error rates in two
ways. Firstly, when $m$ is not large enough, there is no guarantee that
$\widetilde{\boldsymbol{M}}$ is exhaustive (Corollary 3.4). With a positive
probability, one can only recover a proper subspace of the CMS $\mathcal{U}$,
hence the subspace estimation error should be large. When $m$ gets larger, we
expect the estimation error to grow in the order of $\sqrt{m}$ for a fixed $n$
according to the upper-bound in Theorem 3.13. We evaluate the performance of
Algorithm 1 under different choices of $m$, with a fixed sampling budget
$n=100,000$ under a standard Gaussian design. The related results are
presented in Figure 2.

Figure 2: The subspace estimation error
$d(\widehat{\boldsymbol{U}},\boldsymbol{U})$ v.s. the choice of $m$. We
replicate 10 times for each $m$, and plot the mean (the dots) and the standard
error (the error bars) of the subspace estimation errors. When $m$ is small,
the ASGOP $\widetilde{\boldsymbol{M}}$ is not guaranteed to be exhaustive, and
only a proper subspace of $\mathcal{U}$ can be recovered. This results in a
large subspace estimation error. When $m$ is larger than a certain threshold
(c.f., Corollary 3.4), the ASGOP $\widetilde{\boldsymbol{M}}$ is exhaustive
with high probability, and the subspace error has an upper-bound that grows at
the rate of $O(\sqrt{m})$. The result in the figure matches the reasoning
above, as the subspace estimation error drops sharply atk the regime where $m$
is small, and grows gradually for large $m$.

## References

  * [1] D. Babichev and F. Bach, Slice inverse regression with score functions, Electronic Journal of Statistics, 12 (2018), pp. 1507 – 1543, https://doi.org/10.1214/18-EJS1428, https://doi.org/10.1214/18-EJS1428.
  * [2] S. Bendavid, A. Itai, and E. Kushilevitz, Learning by distances, Information and Computation, 117 (1995), pp. 240–250, https://doi.org/https://doi.org/10.1006/inco.1995.1042, https://www.sciencedirect.com/science/article/pii/S0890540185710425.
  * [3] D. R. Brillinger, A generalized linear model with ’gaussian’ regressor variables, A Festschrift for Erich L. Lehmann, K. D. P.J. Bickel and J. Hodges, eds., (1982).
  * [4] L. H. Chen, L. Goldstein, and Q.-M. Shao, Normal Approximation by Stein’s Method, Springer, Berlin, Heidelberg, 2011, https://doi.org/10.1007/978-3-642-15007-4.
  * [5] S. Chen and R. Meka, Learning polynomials in few relevant dimensions, in Conference on Learning Theory, PMLR, 2020, pp. 1161–1227.
  * [6] Y. Chen, Y. Chi, J. Fan, and C. Ma, Spectral methods for data science: A statistical perspective, Foundations and Trends® in Machine Learning, 14 (2021), pp. 566–806, https://doi.org/10.1561/2200000079, https://doi.org/10.1561%2F2200000079.
  * [7] R. Cook and B. Li, Dimension reduction for conditional mean in regression, The Annals of Statistics, 30 (2002), pp. 455 – 474, https://doi.org/10.1214/aos/1021379861, https://doi.org/10.1214/aos/1021379861.
  * [8] R. D. Cook, On the interpretation of regression plots, Journal of the American Statistical Association, 89 (1994), pp. 177–189.
  * [9] R. D. Cook, SAVE: a method for dimension reduction and graphics in regression, Communications in Statistics - Theory and Methods, 29 (2000), pp. 2109–2121, https://doi.org/10.1080/03610920008832598, https://doi.org/10.1080/03610920008832598, https://arxiv.org/abs/https://doi.org/10.1080/03610920008832598.
  * [10] R. D. Cook, Regression graphics: Ideas for studying regressions through graphics, John Wiley & Sons, 2009.
  * [11] A. S. Dalalyan, A. Juditsky, and V. Spokoiny, A new algorithm for estimating the effective dimension-reduction subspace, Journal of Machine Learning Research, 9 (2008), pp. 1647–1678, http://jmlr.org/papers/v9/dalalyan08a.html.
  * [12] A. Damian, J. Lee, and M. Soltanolkotabi, Neural networks can learn representations with gradient descent, in Conference on Learning Theory, PMLR, 2022, pp. 5413–5452.
  * [13] I. Diakonikolas, D. M. Kane, V. Kontonis, and N. Zarifis, Algorithms and sq lower bounds for pac learning one-hidden-layer relu networks, in Proceedings of Thirty Third Conference on Learning Theory, J. Abernethy and S. Agarwal, eds., vol. 125 of Proceedings of Machine Learning Research, PMLR, 09–12 Jul 2020, pp. 1514–1539, https://proceedings.mlr.press/v125/diakonikolas20d.html.
  * [14] R. Dudeja and D. Hsu, Learning single-index models in gaussian space, in Conference On Learning Theory, PMLR, 2018, pp. 1887–1930.
  * [15] M. Hristache, A. Juditsky, J. Polzehl, and V. Spokoiny, Structure adaptive approach for dimension reduction, The Annals of Statistics, 29 (2001), pp. 1537–1566, http://www.jstor.org/stable/2699943 (accessed 2023-04-14).
  * [16] M. Hristache, A. Juditsky, and V. Spokoiny, Direct estimation of the index coefficient in a single-index model, Annals of Statistics, (2001), pp. 595–623.
  * [17] D. Hsu and S. Sabato, Loss minimization and parameter estimation with heavy tails, Journal of Machine Learning Research, 17 (2016), pp. 1–40, http://jmlr.org/papers/v17/14-273.html.
  * [18] M. Kearns, Efficient noise-tolerant learning from statistical queries, J. ACM, 45 (1998), p. 983–1006, https://doi.org/10.1145/293347.293351, https://doi.org/10.1145/293347.293351.
  * [19] T. Klock, A. Lanteri, and S. Vigogna, Estimating multi-index models with response-conditional least squares, Electronic Journal of Statistics, 15 (2021), pp. 589 – 629, https://doi.org/10.1214/20-EJS1785, https://doi.org/10.1214/20-EJS1785.
  * [20] B. Laurent and P. Massart, Adaptive estimation of a quadratic functional by model selection, The Annals of Statistics, 28 (2000), pp. 1302 – 1338, https://doi.org/10.1214/aos/1015957395, https://doi.org/10.1214/aos/1015957395.
  * [21] K.-C. Li, Sliced inverse regression for dimension reduction, Journal of the American Statistical Association, 86 (1991), pp. 316–327, http://www.jstor.org/stable/2290563.
  * [22] K.-C. Li, On principal hessian directions for data visualization and dimension reduction: Another application of stein’s lemma, Journal of the American Statistical Association, 87 (1992), pp. 1025–1039, http://www.jstor.org/stable/2290640 (accessed 2023-04-14).
  * [23] A. Mousavi-Hosseini, S. Park, M. Girotti, I. Mitliagkas, and M. A. Erdogdu, Neural networks efficiently learn low-dimensional representations with SGD, in The Eleventh International Conference on Learning Representations, 2023, https://openreview.net/forum?id=6taykzqcPD.
  * [24] A. Radhakrishnan, D. Beaglehole, P. Pandit, and M. Belkin, Feature learning in neural networks and kernel machines that recursively learn features, arXiv preprint arXiv:2212.13881, (2022).
  * [25] A. M. Samarov, Exploring regression structure using nonparametric functional estimation, Journal of the American Statistical Association, 88 (1993), pp. 836–847, http://www.jstor.org/stable/2290772 (accessed 2023-12-07).
  * [26] G. Stewart and J.-g. Sun, Matrix Perturbation Theory, Academic Press, 1990.
  * [27] S. Trivedi, J. Wang, S. Kpotufe, and G. Shakhnarovich, A consistent estimator of the expected gradient outerproduct, in Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI’14, Arlington, Virginia, USA, 2014, AUAI Press, p. 819–828.
  * [28] Y. Xia, H. Tong, W. K. Li, and L.-X. Zhu, An adaptive estimation of dimension reduction space, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64 (2002), pp. 363–410, https://doi.org/https://doi.org/10.1111/1467-9868.03411, https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.03411, https://arxiv.org/abs/https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/1467-9868.03411.
  * [29] Z. Yang, K. Balasubramanian, Z. Wang, and H. Liu, Estimating high-dimensional non-gaussian multiple index models via stein’s lemma, in Advances in Neural Information Processing Systems, I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds., vol. 30, Curran Associates, Inc., 2017, https://proceedings.neurips.cc/paper_files/paper/2017/file/4db0f8b0fc895da263fd77fc8aecabe4-Paper.pdf.

##  Appendix A SAVE Is Not Exhaustive

In this section, we provide an example where the Sliced Average Variance
Estimator (SAVE) proposed by [9] cannot extract all the relevant directions,
and hence cannot recover the whole central mean subspace (CMS), let alone the
central subspace (CS). In the setting of the multi-index model Eq. 1, SAVE
extracts the index space by estimating the SAVE matrix:

$\boldsymbol{V}=\operatorname*{\mathbb{E}}\left[\left(\operatorname*{\mathbb{E}}\left[\boldsymbol{X}\boldsymbol{X}^{T}\mid\
Y\right]-\boldsymbol{I}\right)^{2}\right].$

One critical condition needed for exhaustiveness here is that
$\operatorname{col}(\boldsymbol{V})$ should contain all the relevant
directions in the index space. This condition is not always satisfied.

Consider a simple noiseless single-index model:

(12) $\displaystyle
Y=f(\boldsymbol{u}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{X})\,,$

where
$\boldsymbol{X}=(X_{1},X_{2})^{{\scriptscriptstyle{\mathsf{T}}}}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{2})$,
and $\boldsymbol{u}=(1,0)^{{\scriptscriptstyle{\mathsf{T}}}}$. Next, we show
that there is a class of link functions with which $\boldsymbol{V}$ is a zero
matrix. Let
$\xi(z_{1},z_{2}):=\operatorname*{\mathbb{E}}_{Z\sim\mathcal{N}(0,1)}[Z^{2}|Z\in[z_{1},z_{2}]]$
for $0\leq z_{1}\leq 1\leq z_{2}$. Note that $\xi$ is continuous, and
monotonically increasing in each of $z_{1}$ and $z_{2}$,
$\xi(1,1)=\operatorname*{\mathbb{E}}[Z^{2}|Z=1]=1$ and
$\xi(0,\infty)=\operatorname*{\mathbb{E}}[Z^{2}|Z\geq 0]=1$. Therefore, there
exists a continuous decreasing function $\nu:[1,\infty)\to(0,1]$, such that
$\nu(1)=1,\lim_{z\to\infty}\nu(z)=0$ and

$\displaystyle\xi(\nu(z),z)=1\text{ for any }z\in[1,\infty)\,.$

Let $f_{0}\colon[0,1]\to\mathbb{R}_{+}$ be any continuous, non-negative, and
monotone increasing function defined on $[0,1]$ with $f_{0}=0$. Then, we
define $f\colon\mathbb{R}\to\mathbb{R}_{+}$ in terms of $f_{0}$ as follows:

(16) $\displaystyle f(z):=\left\\{\begin{array}[]{lll}0,&&z<0;\\\
f_{0}(z),&&0\leq z\leq 1;\\\ f_{0}(\nu(z)),&&1<z<\infty.\end{array}\right.\,$

The following proposition shows that the SAVE matrix $\boldsymbol{V}$ is a
zero matrix with link function $f$.

######  Proposition A.1.

Suppose $(\boldsymbol{X},Y)$ is an example generated from Model (12) with link
function as defined in (16). Then, we have $\boldsymbol{V}=0$.

######  Proof A.2.

It suffices to show that for any $y\geq 0$,
$\operatorname*{\mathbb{E}}[\boldsymbol{X}\boldsymbol{X}^{{\scriptscriptstyle{\mathsf{T}}}}\mathbbm{1}(Y\geq
y)]=\boldsymbol{I}_{2}\operatorname*{\mathbb{P}}(Y\geq y)$. For $y>f_{0}(1)$
or $y\leq 0$, this is vacuously true. For $0<y\leq f_{0}(1)$, we can show from
the definition of $f$ that:

(17) $\displaystyle\\{Y\geq y\\}=\\{X_{1}:y\leq f(X_{1})\leq
f_{0}(1)\\}=\\{f_{0}^{-1}(y)\leq X_{1}\leq\nu^{-1}(f_{0}^{-1}(y))\\}.$

See also Figure 3 for a graphical illustration. Therefore,

$\displaystyle\operatorname*{\mathbb{E}}\left[\boldsymbol{X}\boldsymbol{X}^{{\scriptscriptstyle{\mathsf{T}}}}\Big{|}Y\geq
y\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\begin{bmatrix}X_{1}^{2}&X_{1}X_{2}\\\
X_{1}X_{2}&X_{2}^{2}\end{bmatrix}\Bigg{|}Y\geq y\right]$
$\displaystyle=\begin{bmatrix}\operatorname*{\mathbb{E}}\left[X_{1}^{2}\Big{|}Y\geq
y\right]&\operatorname*{\mathbb{E}}\left[X_{1}\Big{|}Y\geq
y\right]\operatorname*{\mathbb{E}}[X_{2}]\\\
\operatorname*{\mathbb{E}}\left[X_{1}\Big{|}Y\geq
y\right]\operatorname*{\mathbb{E}}[X_{2}]&\operatorname*{\mathbb{E}}[X_{2}^{2}]\end{bmatrix}$
(17)=[ξ(f0-1(y), ν-1(f0-1(y))) 0 0 1 ] = $\boldsymbol{I}$ ,

Finally, we have:

$\displaystyle\operatorname*{\mathbb{E}}\left[\boldsymbol{X}\boldsymbol{X}^{{\scriptscriptstyle{\mathsf{T}}}}\mathbbm{1}(Y\geq
y)\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\boldsymbol{X}\boldsymbol{X}^{{\scriptscriptstyle{\mathsf{T}}}}\mid
Y\geq y\right]\operatorname*{\mathbb{P}}(Y\geq
y)=\boldsymbol{I}\operatorname*{\mathbb{P}}(Y\geq y)\,.$

The proof is completed.

Figure 3: An example plot of the link function $f$ as defined in (16). Here,
we have $\\{Y>y\\}=[z_{1},z_{2}]$ from the plot $\\{x:y\leq f(x)\\}$, where
$z_{1}=f_{0}^{-1}(y)$ and $z_{2}=\nu^{-1}(z_{1})=\nu^{-1}(f_{0}^{-1}(y))$.

##  Appendix B Connection to Local Linear Regression

In this section, we will show that our estimator
$\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}):=\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta};\mathcal{D})$,
which is defined similarly as in Eq. 5 but with a full sample
$\mathcal{D}=\\{(\boldsymbol{X}_{i},Y_{i})\\}_{i=1}^{n}$, is asymptotically
equivalent to the following local linear estimator (LLE) of $\nabla
g(\boldsymbol{\theta})$ with a data-dependent kernel:

$\hat{\boldsymbol{\beta}}^{\operatorname{LLE}}_{h}(\boldsymbol{\theta}):=\operatorname*{arg\,min}_{(\alpha,\boldsymbol{\beta})\in\mathbb{R}\times\mathbb{R}^{d}}\sum_{i=1}^{n}\left(Y_{i}-\alpha-\boldsymbol{\beta}^{{\scriptscriptstyle{\mathsf{T}}}}(\boldsymbol{X}_{i}-\boldsymbol{\theta})\right)^{2}K^{\boldsymbol{\theta}}_{h}(\boldsymbol{X}_{i})\,,$

where
$K^{\boldsymbol{\theta}}_{h}(\boldsymbol{x}):=\rho_{h}(\boldsymbol{x};\boldsymbol{\theta})$.

######  Proposition B.1.

For each $\boldsymbol{\theta}\in\mathbb{R}^{d}$,

$\left\lVert\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta})-\hat{\boldsymbol{\beta}}^{\operatorname{LLE}}_{h}(\boldsymbol{\theta})\right\rVert\overset{\operatorname*{\mathbb{P}}}{\rightarrow}{0}\,,$

here, the notation $\overset{\operatorname*{\mathbb{P}}}{\rightarrow}$
represents convergence in probability.

######  Proof B.2.

By standard solution for the weighted least square error problems gives:

(18)
$\displaystyle\hat{\boldsymbol{\beta}}^{\operatorname{LLE}}_{h}(\boldsymbol{\theta})=\boldsymbol{I}_{d}^{\backslash
1}(\boldsymbol{X}_{\boldsymbol{\theta}}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{W}^{\boldsymbol{\theta}}_{h}\boldsymbol{X}_{\boldsymbol{\theta}})^{-1}\boldsymbol{X}_{\boldsymbol{\theta}}^{{\scriptscriptstyle{\mathsf{T}}}}\boldsymbol{W}^{\boldsymbol{\theta}}_{h}\boldsymbol{Y}$

where $\boldsymbol{I}_{d+1}^{\backslash 1}\in\mathbb{R}^{d\times(d+1)}$ is
obtained by removing the first row of a $(d+1)\times(d+1)$ identity matrix
$\boldsymbol{I}_{d+1}$,

$\boldsymbol{X}_{\boldsymbol{\theta}}=\begin{bmatrix}1&(\boldsymbol{X}_{1}-\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}\\\
\vdots&\vdots\\\
1&(\boldsymbol{X}_{n}-\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}}\end{bmatrix},\quad\boldsymbol{W}_{\boldsymbol{\theta}}=\operatorname{diag}(K^{\boldsymbol{\theta}}_{h}(\boldsymbol{X}_{1}),\ldots,K^{\boldsymbol{\theta}}_{h}(\boldsymbol{X}_{n})),\quad\boldsymbol{Y}=\begin{bmatrix}Y_{1}\\\
\vdots\\\ Y_{n}\end{bmatrix}.$

One can show from Eq. 18 and the block matrix inversion formula
(LABEL:lemma:matrix) that

$\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta})=h^{-2}\left(\hat{\Sigma}(\boldsymbol{\theta})\,\hat{\boldsymbol{\beta}}^{\operatorname{LLE}}_{h}(\boldsymbol{\theta})+\hat{\mu}_{y}(\boldsymbol{\theta})\hat{\boldsymbol{\mu}}_{\boldsymbol{x}}(\boldsymbol{\theta})\right),$

where

$\displaystyle\hat{\Sigma}(\boldsymbol{\theta})$
$\displaystyle:=n^{-1}\sum_{i=1}^{n}\rho_{h}(\boldsymbol{X}_{i};\boldsymbol{\theta})(\boldsymbol{X}_{i}-\boldsymbol{\theta})(\boldsymbol{X}_{i}-\boldsymbol{\theta})^{{\scriptscriptstyle{\mathsf{T}}}};$
$\displaystyle\hat{\boldsymbol{\mu}}_{\boldsymbol{x}}(\boldsymbol{\theta})$
$\displaystyle:=n^{-1}\sum_{i=1}^{n}\rho_{h}(\boldsymbol{X}_{i};\boldsymbol{\theta})(\boldsymbol{X}_{i}-\boldsymbol{\theta});$
$\displaystyle\hat{\mu}_{y}(\boldsymbol{\theta})$
$\displaystyle:=n^{-1}\sum_{i=1}^{n}\rho_{h}(\boldsymbol{X}_{i};\boldsymbol{\theta})Y_{i}.$

By the weak law of large numbers,
$\hat{\Sigma}\overset{\operatorname*{\mathbb{P}}}{\rightarrow}h^{2}\boldsymbol{I}_{d},$
$\hat{\boldsymbol{\mu}}_{\boldsymbol{x}}\overset{\operatorname*{\mathbb{P}}}{\rightarrow}0,$
$\hat{\mu}_{y}\overset{\operatorname*{\mathbb{P}}}{\rightarrow}\operatorname*{\mathbb{E}}_{\tilde{\boldsymbol{X}}\sim\mathcal{N}(\boldsymbol{0}_{d},\boldsymbol{I}_{d})}[f(\boldsymbol{U}^{{\scriptscriptstyle{\mathsf{T}}}}\tilde{\boldsymbol{X}})].$
The proposition follows and the proof is completed.

##  Appendix C Estimation Bias Under $P_{\boldsymbol{X}}$ with a Bounded
Support

In this section, we discuss the performance of our estimator when 1(2) is not
satisfied and the Lebesgue density $p$ of $P_{\boldsymbol{X}}$ has a bounded
support $\\{\boldsymbol{x}\in\mathbb{R}^{d}:\lVert\boldsymbol{x}\rVert\leq
K\\}$ for some $K<\infty$. In this case, we note that the smoothed gradient
estimator
$\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j};\mathcal{D}_{j,l})$ as
defined in Eq. 5 is biased for
${\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j})$. In particular, for a
fixed $\boldsymbol{\theta}_{j}$,

$\displaystyle\operatorname*{\mathbb{E}}[\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j};\mathcal{D}_{j,l})]$
$\displaystyle={h^{-2}}\operatorname*{\mathbb{E}}[\rho_{h}(\boldsymbol{X};\boldsymbol{\theta})\cdot
Y\cdot(\boldsymbol{X}-\boldsymbol{\theta}_{j})]$
$\displaystyle={h^{-2}}\int_{\lVert\boldsymbol{X}\rVert\leq
K}\rho_{h}(\boldsymbol{X};\boldsymbol{\theta}_{j})\cdot
g(\boldsymbol{X})\cdot(\boldsymbol{X}-\boldsymbol{\theta}_{j})dP_{\boldsymbol{X}}$
$\displaystyle={h^{-2}}\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[\mathbbm{1}(\lVert\boldsymbol{Z}\rVert\leq
K)\cdot
g(\boldsymbol{Z})\cdot(\boldsymbol{Z}-\boldsymbol{\theta}_{j})\right]\,.$

Hence, the bias is given as:

$\displaystyle B_{h,j}$
$\displaystyle:=\operatorname*{\mathbb{E}}[\hat{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j};\mathcal{D}_{j,l})]-{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j})$
(19)
$\displaystyle=-{h^{-2}}\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[\mathbbm{1}(\lVert\boldsymbol{Z}\rVert>K)\cdot
g(\boldsymbol{Z})\cdot(\boldsymbol{Z}-\boldsymbol{\theta}_{j})\right]$

Intuitively, $\lVert B_{h,j}\rVert$ should shrink to zero as $h\rightarrow 0$
when $K\gg\lVert\boldsymbol{\theta}_{j}\rVert$, which holds with high
probability if we choose the variance $\sigma_{\theta}^{2}$ of
$\boldsymbol{\theta}_{j}$’s is small compared to $K$. On the other hand, it is
clear that if $\sup_{j\in[m]}\lVert B_{h,j}\rVert=o(n^{-1/2})$, then the bias
will not affect the rate in the upper-bound for subspace estimation error. The
following proposition shows how small $h$ needs to be to guarantee such small
biases.

######  Proposition C.1.

Suppose that $\lVert\boldsymbol{\theta}_{j}\rVert\leq h$. Let
$h=K(\sqrt{d}+\sqrt{\log n})^{-1}/4$, then we have the bias
$\lVert\operatorname*{\mathbb{E}}\hat{\boldsymbol{\beta}}_{j,\ell}-{\boldsymbol{\beta}}_{h}(\boldsymbol{\theta}_{j})\rVert=O(\log
n\cdot n^{-2/3})$.

######  Proof C.2.

By Cauchy-Schwarz Inequality and Appendix C, we have:

$\displaystyle\lVert B_{h,j}\rVert$
$\displaystyle\leq{h^{-2}}\Bigg{(}\operatorname*{\mathbb{E}}_{\widetilde{\boldsymbol{Z}}\sim\mathcal{N}(h^{-2}\boldsymbol{\theta}_{j},\boldsymbol{I}_{d})}\left[\mathbbm{1}(\lVert\widetilde{\boldsymbol{Z}}\rVert>Kh^{-1})\right]\cdot\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[g(\boldsymbol{Z})^{3}\right]$
$\displaystyle\qquad\qquad\qquad\cdot\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[\lVert\boldsymbol{Z}-\boldsymbol{\theta}_{j}\rVert^{3}\right]\Bigg{)}^{1/3}\,.$

We will upper-bound each of the terms in the parentheses. By 4 and 5, it can
be verified that there exists some $C>0$, independent of $h$ and $n$, such
that:

(20)
$\displaystyle\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[g(\boldsymbol{Z})^{3}\right]\leq
C\cdot\mu_{\rho}\cdot\mu_{f}\,.$

The condition $\lVert\boldsymbol{\theta}_{j}\rVert\leq h$ implies that:

(21)
$\displaystyle\operatorname*{\mathbb{E}}_{\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{\theta}_{j},h^{2}\boldsymbol{I}_{d})}\left[\lVert\boldsymbol{Z}-\boldsymbol{\theta}_{j}\rVert^{3}\right]\leq
C\cdot h^{3}\,.$

Finally, we have under the assumption
$\lVert\boldsymbol{\theta}_{j}\rVert<h<K/2$ that:

(22)
$\displaystyle\operatorname*{\mathbb{E}}_{\widetilde{\boldsymbol{Z}}\sim\mathcal{N}(h^{-2}\boldsymbol{\theta}_{j},\boldsymbol{I}_{d})}\left[\mathbbm{1}(\lVert\widetilde{\boldsymbol{Z}}\rVert>Kh^{-1})\right]$
$\displaystyle\leq\operatorname*{\mathbb{P}}\left(\lVert\mathcal{N}(\boldsymbol{0}_{d},\boldsymbol{I}_{d})\rVert>\frac{Kh^{-1}}{2}\right)=\operatorname*{\mathbb{P}}\left(\chi^{2}_{d}>\frac{K^{2}h^{-2}}{4}\right).$

Consider the following concentration inequality for the chi-squared
distribution:

######  Lemma C.3 (Concentration for $\chi^{2}_{d}$, adapted from Lemma 1 of
[20]).

Let $W$ be a random variable follows the centered chi-squared distribution
with degree of freedom $d$, i.e., $W\sim\chi^{2}_{d}$, then for any $t>0$:

$\operatorname*{\mathbb{P}}(W>d+2\sqrt{dt}+2t)\leq\exp(-t)\,.$

By Lemma C.3, when $h=K(\sqrt{d}+\sqrt{\log n})^{-1}/4$, we have:

(23)
$\displaystyle\operatorname*{\mathbb{P}}\left(\chi^{2}_{d}>\frac{K^{2}h^{-2}}{4}\right)\leq
n^{-2}\,.$

From Eq. 20 to Eq. 23, we have:

$\displaystyle\lVert B_{h,j}\rVert\leq
Ch^{-2}(n^{-2}\cdot\mu_{\rho}\cdot\mu_{f}\cdot
h^{3})^{1/3}\leq\tilde{C}\cdot\log n\cdot n^{-2/3}\,.$

where $\tilde{C}$ is independent of $n$.

######  Remark C.4.

The condition $\lVert\boldsymbol{\theta}_{j}\rVert\leq h$ in Proposition C.1
holds uniformly for $j\in[m]$, with probability at least $1-\delta$ with the
choice $\sigma_{\theta}\lesssim h\cdot\log(m/\delta)$. This is compatible with
our choice of $\sigma_{\theta}$ in the main theorem.

See pages - of supplementary.pdf