Dimension Reduction for Fréchet Regression

Let $\mathfrak{F}$ be a family of measurable functions $f:\Omega_{Y}\to\mathbb{R}$, and for an $f\in\mathfrak{F}$, let $\mathcal{S}_{E[f(Y)|X]}$ be the central mean subspace of $f(Y)$ versus $X$. As mentioned in Section 1, we use two types of ensembles to recover the Fréchet central subspace. The first type is any $\mathfrak{F}$ that satisfies

This is the same ensemble as that in Yin & Li (2011), except that, here, the right-hand side is the Fréchet central subspace. Relation (2) allows us to recover the Fréchet central subspace $\mathcal{S}_{Y|X}$ by a collection of classical central mean subspaces. We call a class $\mathfrak{F}$ that satisfies (2) a CMS-ensemble. The second type of ensembles is any family $\mathfrak{F}$ that satisfies

which we call a CS-ensemble. Proposition 1 shows that a CMS ensemble is a CS-ensemble.

We next develop a sufficient condition for an $\mathfrak{F}$ to be a CMS-ensemble and hence also a CS-ensemble. Let $\mathfrak{B}=\{I_{B}:B\ \text{is Borel set in}\ \Omega_{Y}\}$ be the family of measurable indicator functions on $\Omega_{Y}$, and let $\mathrm{span}(\mathfrak{F})=\left\{\sum_{i=1}^{k}\alpha_{i}f_{i}:k\in\mathbb{N},\alpha_{1},\dots,\alpha_{k}\in\mathbb{R},f_{1},\dots,f_{k}\in\mathfrak{F}\right\}$ be the linear span of $\mathfrak{F}$, where $\mathbb{N}=\{1,2,\dots\}$. Yin & Li (2011) showed that if $\mathfrak{F}$ is a subset of $L_{2}(P_{Y})$ that is dense in $\mathfrak{B}$, then (2) holds for the classical $\mathcal{S}_{Y|X}$. Here, we generalize that result to our setting by requiring only $\mathrm{span}(\mathfrak{F})$ to be dense in $\mathfrak{B}$.

To construct a CMS-ensemble, we resort to the notion of the universal kernel. Let $C(\Omega_{Y})$ be the family of continuous real-valued functions on $\Omega_{Y}$. When $\Omega_{\scriptscriptstyle Y}$ is compact, Steinwart (2001) defined a continuous kernel $\kappa$ as universal (we refer to it as c-universal) if its associated RKHS ${\cal{H}}_{\scriptscriptstyle Y}$ is dense in $C(\Omega_{\scriptscriptstyle Y})$ under the uniform norm. To relax the compactness assumption, Micchelli et al. (2006) proposed the following notion of universality, which is referred to cc-univsersal in Sriperumbudur et al. (2011). For any compact set $K\subseteq\Omega_{\scriptscriptstyle Y}$, let ${\cal{H}}_{\scriptscriptstyle Y}(K)$ be the RKHS generated by $\{\kappa(\cdot,y):y\in K\}$. We should note that a member $f$ of ${\cal{H}}_{\scriptscriptstyle Y}(K)$ is supported on $\Omega_{\scriptscriptstyle Y}$, rather than $K$. Let $f|K$ denote the restriction of $f$ on $K$, and $C(K)$ the class of all continuous functions with respect to the topology in $(\Omega_{\scriptscriptstyle Y},d)$ restricted on $K$.

When $\Omega_{\scriptscriptstyle Y}$ is compact, Sriperumbudur et al. (2011) showed that two notions of universality are equivalent. In the following, we look into the conditions under which a metric space has a cc-universal kernel and how to construct such a kernel when it does.

Micchelli et al. (2006) showed that when $\Omega_{\scriptscriptstyle Y}=\mathbb{R}^{d}$, many standard kernels, including Laplacian kernels and Gaussian RBF kernels, are cc-universal. Unfortunately, when $\Omega_{Y}$ is a general metric space, direct extension of these types of kernels, for example, $k(y,y^{\prime})=\exp(-\gamma d(y,y^{\prime})^{2})$, are no longer guaranteed to be cc-universal. Christmann & Steinwart (2010) showed that for compact $\Omega_{\scriptscriptstyle Y}$, if there exists a separable Hilbert space $\mathcal{H}$ and a continuous injection $\rho:\Omega_{Y}\to\mathcal{H}$, then for any analytic function $F:\mathbb{R}\to\mathbb{R}$ whose Taylor series at zero has strictly positive coefficients, the function $\kappa(y,y^{\prime})=F(\langle\rho(y),\rho(y^{\prime})\rangle_{\mathcal{H}})$ defines a c-universal kernel on $\Omega_{Y}$. They also provide an analogous definition of the Gaussian-type kernel in the above case. We extend this result to construct cc-universal kernels on non-compact metric space. The proof is given in the Supplementary Material.

As an example, Corollary 1 shows that the Gaussian-type kernel is cc-universal on $\Omega_{Y}$.

The second part of Corollary 1 is straightforward since an isometry is an injection. Similar results can be established for Laplacian-type kernel $\kappa_{\gamma}(y,y^{\prime})=\exp(-\gamma\|\rho(y)-\rho(y^{\prime})\|_{\mathcal{H}})$.

The continuous embedding condition in Proposition 2 covers several metric spaces often encountered in statistical applications. Section 3 employs it to construct cc-universal kernels on the space of univariate distributions endowed with Wasserstein-2 distance, correlation matrices endowed with Frobenius distance, and spheres endowed with geodesic distance.

By using the notion of regular probability measure, we connect the cc-universal kernel on $(\Omega_{Y},d)$ with the CMS-ensemble, which is the theoretical foundation of our method. Recall that a measure $P_{\scriptscriptstyle Y}$ on $(\Omega_{\scriptscriptstyle Y},d)$ is regular if, for any Borel subset $B\subseteq\Omega_{\scriptscriptstyle Y}$ and any $\varepsilon>0$, there is a compact set $K\subseteq B$ and an open set $G\supseteq B$, such that $P(G\backslash K)<\varepsilon$.

The proof of Theorem 1 is given in the Supplementary Material. Condition (2), which requires $P_{\scriptscriptstyle Y}$ to be regular, is quite mild: it is known that any Borel measure on a complete and separable metric space is regular (see Granirer (1970, Chapter 2: Theorem 1.2, Theorem 3.2)). Thus, a sufficient condition of Condition (2) is $(\Omega_{\scriptscriptstyle Y},d)$ being complete and separable, which is satisfied by all the metric spaces we consider. Specifically, note that if $M$ is separable and complete, then so is the Wasserstein-2 space $W_{\scriptscriptstyle 2}(M)$ (Panaretos & Zemel, 2020, Proposition 2.2.8, Theorem 2.2.7). Therefore, $W_{\scriptscriptstyle 2}(\mathbb{R})$ is complete and separable. Similarly, the SPD matrix space endowed with Frobenius distance and the sphere endowed with geodesic distance are both completely separable metric spaces. Furthermore, the Gaussian kernel and Laplacian kernel we considered satisfy Condition (1) in Theorem 1.

Thus, Proposition 2 and Theorem 1 provide a general mechanism to construct the CMS-ensemble over any separable and complete metric space without a linear structure, provided it can be continuously embedded in a separable Hilbert space. For the case where multiple cc-universal kernels exist, we design a cross-validation framework in Section 6 to choose the kernel type and the bandwidth $\gamma$.

Let $I$ be $\mathbb{R}$ or a closed interval of $\mathbb{R}$, $\mathcal{B}(I)$ the $\sigma$-field of Borel subsets of $I$, and $\mathcal{P}(I)$ the collection of all probability measures on $(I,\mathcal{B}(I))$. The Wasserstein space $\mathcal{W}_{2}(I)$ is defined as the subset of $\mathcal{P}(I)$ with finite second moment, that is, $\mathcal{W}_{2}(I)=\{\mu\in\mathcal{P}(I):\int_{I}t^{2}\,d\mu(t)<\infty\},$ endowed with the quadratic Wasserstein distance $d_{\mathrm{W}}(\mu_{1},\mu_{2})=\left(\int_{0}^{1}\left[F_{\mu_{1}}^{-1}(s)-F_{\mu_{2}}^{-1}(s)\right]^{2}\,ds\right)^{1/2},$ where $\mu_{1}$ and $\mu_{2}$ are members of $\mathcal{W}_{2}(I)$ and $F^{-1}_{\mu_{1}}$ and $F^{-1}_{\mu_{2}}$ are the quantile functions of $\mu_{1}$ and $\mu_{2}$, which we assume to be well defined. This distance can be equivalently written as $d_{\mathrm{W}}(\mu_{1},\mu_{2})=\left(\int_{I}\left[F_{\mu_{1}}^{-1}\circ F_{\mu_{2}}(t)-t\right]^{2}\,d\mu_{2}(t)\right)^{1/2}.$ The set $\mathcal{W}_{2}(I)$ endowed with $d_{\mathrm{W}}$ is a metric space with a formal Riemannian structure (Ambrosio et al., 2004).

Here, we present some basic results that characterize $\mathcal{W}_{2}(I)$, whose proofs can be found, for example, in Ambrosio et al. (2004) and Bigot et al. (2017). For $\mu_{1},\,\mu_{2}\in\mathcal{W}_{2}(I)$, we say that a $\mathcal{B}(I)$-measurable map $r:I\to I$ transports $\mu_{1}$ to $\mu_{2}$ if $\mu_{2}=\mu_{1}\circ r^{-1}$. This relation is often written as $\mu_{2}={r}_{\#}\mu_{1}$. Let $\mu_{0}\in\mathcal{W}_{2}(I)$ be a reference measure with a continuous $F_{\mu_{0}}$. The tangent space at $\mu_{0}$ is $T_{\mu_{0}}=\mathrm{cl}_{L_{2}(\mu_{0})}{\{\lambda(F_{\mu}^{-1}\circ F_{\mu_{0}}-\mathrm{id}):\mu\in\mathcal{W}_{2}(I),\lambda>0\}},$ where, for a set $A\subseteq L_{2}(\mu_{0})$, $\mathrm{cl}_{L_{2}(\mu_{0})}(A)$ denotes the $L_{2}(\mu_{0})$-closure of $A$. The exponential map $\exp_{\mu_{0}}$ from $T_{\mu_{0}}$ to $\mathcal{W}_{2}(I)$, defined by $\exp_{\mu_{0}}({r})=({r}+\mathrm{id})_{\#}\mu_{0}$, is surjective. Therefore its inverse, $\log_{\mu_{0}}:\mathcal{W}_{2}(I)\to T_{\mu_{0}}$, defined by $\log_{\mu_{0}}(\mu)=F_{\mu}^{-1}\circ F_{\mu_{0}}-\mathrm{id}$, is well defined on $\mathcal{W}_{2}(I)$. It is well known that the exponential map restricted to the image of $\log$ map, denoted as $\exp_{\mu_{0}}|_{\log_{\mu_{0}}(\mu)(\mathcal{W}_{2}(I))}$, is an isometric homeomorphism (Bigot et al., 2017). Therefore, $\log_{\mu_{0}}$ is a continuous injection from $\mathcal{W}_{2}(I)$ to $L_{2}(\mu_{0})$. We can then construct CMS-ensembles using the general constructive method provided by Theorem 1 and Proposition 2. The next proposition gives two such constructions, where the subscripts “G” and “L” for the two kernels refer to “Gaussian” and “Laplacian”, respectively.

We first introduce some notations. Let $\mathrm{Sym}(r)$ be the set of $r\times r$ invertible symmetric matrices with real entries and $\mathrm{Sym}^{+}(r)$ the set of $r\times r$ symmetric positive definite (SPD) matrices. For any $Y\in\mathbb{R}^{\scriptscriptstyle r\times r}$, the matrix exponential of $Y$ is defined as the infinite power series $\exp(Y)=\sum_{k=0}^{\infty}Y^{k}/k!$. For any $X\in\mathrm{Sym}^{+}(r)$, the matrix logarithm of $X$ is defined as any $r\times r$ matrix $Y$ such that $\exp(Y)=X$ and denoted by $\log(X)$.

Let $d_{\mathrm{F}}$ be the Frobenius metric. Then $(\mathrm{Sym}(r),d_{\mathrm{F}})$ is a metric space continuously embedded by identity mapping in ${\mathrm{Sym}}(r)$, which is a Hilbert space with the Frobenius inner product $\langle A,B\rangle=\mathrm{tr}(A^{\mbox{\tiny{\sf T}}}B)$. Also, the identity mapping $\mathrm{id}:\mathrm{Sym}^{+}(r)\to\mathrm{Sym}(r)$ is obviously isometric. Therefore, by Corollary 1, the two types of radial basis function kernels for Wasserstein space can be similarly extended to $\mathrm{Sym}^{+}(r)$. That is, let $\kappa_{G}(y,y^{\prime})=\exp(-\gamma d_{\mathrm{F}}(y,y^{\prime})^{2})$ and $\kappa_{L}(y,y^{\prime})=\exp(-\gamma d_{\mathrm{F}}(y,y^{\prime})),$ then $\mathfrak{F}_{G}=\{\kappa_{G}(y,y^{\prime}),y^{\prime}\in\mathrm{Sym}^{+}(r)\}$ and $\mathfrak{F}_{L}=\{\kappa_{L}(y,y^{\prime}),y^{\prime}\in\mathrm{Sym}^{+}(r)\}$ are CMS-ensembles.

Another widely used metric over $\mathrm{Sym}^{+}(r)$ is the log-Euclidean distance that is defined as $d_{\log}(Y_{1},Y_{2})=\|\log(Y_{1})-\log(Y_{2})\|_{\mathrm{F}}.$ Basically, it pulls the Frobenius metric on $\mathrm{Sym}(r)$ back to $\mathrm{Sym}^{+}(r)$ by the matrix logarithm map. The matrix logarithm $\log(\cdot)$ is a continuous injection to Hilbert $\mathrm{Sym}(r)$. By Corollary 1, the two types of radial basis function kernels $\kappa_{G,\log}(y,y^{\prime})=\exp(-\gamma d_{\log}(y,y^{\prime})^{2})$ and $\kappa_{L,\log}(y,y^{\prime})=\exp(-\gamma d_{\log}(y,y^{\prime}))$ are cc-universal. Then, $\mathfrak{F}_{G,\log}=\{\kappa_{G,\log}(y,y^{\prime}),y^{\prime}\in\mathcal{W}_{2}(I)\}$ and $\mathfrak{F}_{L,\log}=\{\kappa_{L,\log}(y,y^{\prime}),y^{\prime}\in\mathcal{W}_{2}(I)\}$ are CMS-ensembles.

Consider the random vector taking values in the sphere $\mathbb{S}^{n}=\{x\in\mathbb{R}^{n}:\|x\|=1\}$. To respect the nonzero curvature of $\mathbb{S}^{n}$, the geodesic distance $d_{g}(Y_{1},Y_{2})=\arccos(Y_{1}^{\mbox{\tiny{\sf T}}}Y_{2})$, which is derived from its Riemannian geometry, is often used rather than the Euclidean distance. However, the popular Gaussian-type RBF kernel $\kappa_{G}(y,y^{\prime})=\exp(-\gamma d_{g}(y,y^{\prime})^{2})$ is not positive definite on $\mathbb{S}^{n}$ (Jayasumana et al., 2013). In fact, Feragen et al. (2015) proved that for complete Riemannian manifold $M$ with its associated geodesic distance $d_{g}$, $\kappa_{G}(y,y^{\prime})=\exp(-\gamma d_{g}(y,y^{\prime})^{2})$ is positive semidefinite only if $M$ is isometric to a Euclidean space. Honeine & Richard (2010) and Jayasumana et al. (2013) proved that the Laplacian-type kernel $\kappa_{L}(y,y^{\prime})=\exp(-\gamma d_{g}(y,y^{\prime}))$ is positive definite on the sphere $\mathbb{S}^{n}$. We show in the following proposition that $\kappa_{L}(y,y^{\prime})$ is cc-universal.

We first develop a general class of Fréchet SDR estimators based on the ensembled moment estimators of the CMS, such as the OLS, PHD, and IHT. Let ${\cal{P}}_{\scriptscriptstyle XY}$ be the collection of all distributions of $(X,Y)$, and let $M:{\cal{P}}_{\scriptscriptstyle XY}\to\mathbb{R}^{\scriptscriptstyle p\times p}$ be a measurable function to be used as an estimator of the Fréchet central subspace ${\cal{S}}_{\scriptscriptstyle Y|X}$. A function defined on ${\cal{P}}_{\scriptscriptstyle XY}$ is called statistical functional; see, for example, Chapter 9 of Li (2018). In the SDR literature, such a function is also called a candidate matrix (Ye & Weiss, 2003). Let $F_{\scriptscriptstyle XY}$ be a generic member of ${\cal{P}}_{\scriptscriptstyle XY}$, $F_{\scriptscriptstyle XY}^{\scriptscriptstyle(0)}$ the true distribution of $(X,Y)$, and $\hat{F}_{\scriptscriptstyle XY}^{\scriptscriptstyle(n)}$ the empirical distribution of $(X,Y)$ based on an i.i.d. sample $(X_{\scriptscriptstyle 1},Y_{\scriptscriptstyle 1}),\ldots,(X_{\scriptscriptstyle n},Y_{\scriptscriptstyle n})$. Extending the terminology of classical SDR (see, for example, Li 2018, Chapter 2), we say that the estimate $M(\hat{F}_{\scriptscriptstyle XY}^{\scriptscriptstyle(n)})$ is unbiased if $M(F_{\scriptscriptstyle XY}^{\scriptscriptstyle(0)})\subseteq{\cal{S}}_{\scriptscriptstyle Y|X}$, exhaustive if $M(F_{\scriptscriptstyle XY}^{\scriptscriptstyle(0)})\supseteq{\cal{S}}_{\scriptscriptstyle Y|X}$, and Fisher consistent if $M(F_{\scriptscriptstyle XY}^{\scriptscriptstyle(0)})={\cal{S}}_{\scriptscriptstyle Y|X}$. We refer to $M$ as the Fréchet candidate matrix.

Suppose we are given a CMS-ensemble $\mathfrak{F}$. Let $M_{\scriptscriptstyle 0}:{\cal{P}}_{\scriptscriptstyle XY}\times\mathfrak{F}\to\mathbb{R}^{\scriptscriptstyle p\times p}$ be a function to be used as an estimator of ${\cal{S}}_{\scriptscriptstyle E[f(Y)|X]}$ for each $f$. This is not a statistical functional in the classical sense, as it involves an additional set $\mathfrak{F}$. So, we redefine unbiasedness, exhaustiveness, and Fisher consistency for this type of augmented statistical functional.

Note that $M_{\scriptscriptstyle 0}(\cdot,f)$ is an estimator of the classical central mean subspace ${\cal{S}}_{\scriptscriptstyle E[f(Y)|X]}$, as $f(Y)$ is a random number rather than a random object. We refer to $M_{\scriptscriptstyle 0}$ as the ensemble candidate matrix or, when confusion is possible, the CMS-ensemble candidate matrix. Our goal is to construct a Fréchet candidate matrix $M:{\cal{P}}_{\scriptscriptstyle XY}\to\mathbb{R}^{\scriptscriptstyle p\times p}$ from the ensemble candidate matrix $M_{\scriptscriptstyle 0}:{\cal{P}}_{\scriptscriptstyle XY}\times\mathfrak{F}\to\mathbb{R}^{\scriptscriptstyle p\times p}$. To do so, we assume $\mathfrak{F}$ is of the form $\{\kappa(\cdot,y):y\in\Omega_{\scriptscriptstyle Y}\}$, where $\kappa:\Omega_{\scriptscriptstyle Y}\times\Omega_{\scriptscriptstyle Y}\to\mathbb{R}$ is a cc-universal kernel. Given such an $\mathfrak{F}$ and $M_{\scriptscriptstyle 0}$, we define $M$ as follows

where $F_{\scriptscriptstyle Y}$ is the distribution of $Y$ derived from $F_{\scriptscriptstyle XY}$.

We now adapt several estimates for the classical central mean subspace to the estimation of Fréchet SDR: the ordinary least squares (OLS; Li & Duan 1989), the principal Hessian directions (PHD; Li 1992), and the Iterative Hessian Transformation (IHT; Cook & Li 2002). These estimates are based on sample moments and require additional conditions on the predictor $X$ for their unbiasedness. Specifically, we make the following assumptions :

Under the first assumption, the ensemble OLS and IHT are unbiased for estimating the Fréchet central subspace; under both assumptions, the ensemble PHD is unbiased for estimating ${\cal{S}}_{\scriptscriptstyle Y|X}$. More detailed discussions on the unbiasedness and fisher consistency of ensemble estimators are presented in Section 4.4. In practice, the two assumptions above cannot be checked directly since we do not know $\beta$. However, as was shown by Eaton (1986), if Assumption 1 holds for all $\beta$, then the distribution of $X$ is elliptical, and vice versa. If further $X$ is multivariate normal, then Assumption 2 is satisfied. Currently, the scatter plot matrix is the most commonly used empirical method to check the elliptical distribution assumption. If non-elliptical features are observed, one can use marginal transformations of the predictors, such as the Box-Cox transformation, to mitigate the non-ellipticity problem. Furthermore, in practice, the SDR methods that require ellipticity usually still work reasonably well even when the elliptical distribution assumption is violated. This occurs particularly when the dimension $p$ of $X$ is high. See Hall & Li (1993) and Li & Yin (2007) for the theoretical supports. Our simulation results in Section 6 support this phenomenon.

It is most convenient to construct these ensemble estimators using standardized predictors. The theoretical basis for doing so is an equivariant property of the Fréchet central subspace, as stated in the next proposition.

The proof is essentially the same as that for the classical central subspace (see, for example, Li, 2018, page 24), and is omitted. Using this property, we first transform $X$ to $Z=\mathrm{var}(X)^{\scriptscriptstyle-1/2}(X-EX)$, estimate the Fréchet central subspace ${\cal{S}}_{\scriptscriptstyle Y|Z}$, and then transform it by $\mathrm{var}(X)^{\scriptscriptstyle-1/2}{\cal{S}}_{\scriptscriptstyle Y|Z}$, which is the same as ${\cal{S}}_{\scriptscriptstyle Y|X}$. The candidate matrices $M_{\scriptscriptstyle 0}$ and $M$ for ensemble OLS, PHD, and IHT are formulated in Remark 1. Detailed motivation for each can be found in Li (2018, Chapter 8). The sample estimates can then be constructed by replacing the expectations in $M_{\scriptscriptstyle 0}$ and $M$ with sample moments whenever possible. Algorithm 1 summarize the steps to implement an ensembled moment estimator, where $\kappa_{\scriptscriptstyle c}(y,y^{\prime})$ stands for the centered kernel $\kappa(y,y^{\prime})-E_{\scriptscriptstyle n}\kappa(Y,y^{\prime})$.

In this subsection, we adapt the OPG (Xia et al. 2002), a popular method for estimating the classical CMS based on nonparametric forward regression, to the estimation of the Fréchet central subspace, which do not require LCM and CCV conditions. The adaption of another forward regression method MAVE is similar and presented in Section S.3.2 of the Supplementary Material. The framework of the statistical functional $M_{\scriptscriptstyle 0}(F_{\scriptscriptstyle XY},f)$ is no longer sufficient to cover this case because we now have a tuning parameter here. So, we adopt the notion of tuned statistical functional in Section 11.2 of Li (2018) to accommodate a tuning parameter.

Let ${\cal{P}}_{\scriptscriptstyle XY}$, $F_{\scriptscriptstyle XY}$, $F_{\scriptscriptstyle XY}^{\scriptscriptstyle(0)}$ and $\hat{F}_{\scriptscriptstyle XY}^{\scriptscriptstyle(n)}$ be as defined in Section 4.1. For simplicity, we assume the tuning parameter $h$ to be a scalar, but it could also be a vector. Given a CMS-ensemble $\mathfrak{F}$, let $T_{0}:{\cal{P}}_{\scriptscriptstyle XY}\times\mathfrak{F}\times\mathbb{R}\to\mathbb{R}^{\scriptscriptstyle p\times p}$ be a tuned functional to be used as an estimator of ${\cal{S}}_{\scriptscriptstyle E[f(Y)|X]}$ for each $f$. We refer to $T_{0}$ as the ensemble-tuned candidate matrix. The unbiasedness, exhaustiveness, and Fisher consistency of $T_{0}$ are defined as follows.

Given $\mathfrak{F}=\{\kappa(\cdot,y):y\in\Omega_{Y}\}$ and $T_{\scriptscriptstyle 0}$, we define the tuned Fréchet candidate matrix $T:{\cal{P}}_{\scriptscriptstyle XY}\times\mathbb{R}\to\mathbb{R}^{p\times p}$ as $T(F_{\scriptscriptstyle XY},h)=\int_{\Omega_{Y}}T_{0}(F_{\scriptscriptstyle XY},\kappa(\cdot,y),h)dF_{\scriptscriptstyle Y}(y).$ We say that the estimate $T(T^{\scriptscriptstyle(n)}_{\scriptscriptstyle XY},h)$ is unbiased if $\mathrm{span}(\lim_{\scriptscriptstyle h\to 0}T(F^{\scriptscriptstyle(0)}_{\scriptscriptstyle XY},h))\subseteq{\cal{S}}_{\scriptscriptstyle Y|X}$, exhaustive if $\mathrm{span}(\lim_{\scriptscriptstyle h\to 0}T(F^{\scriptscriptstyle(0)}_{\scriptscriptstyle XY},h))\supseteq{\cal{S}}_{\scriptscriptstyle Y|X}$, and Fisher consistent if $\mathrm{span}(\lim_{\scriptscriptstyle h\to 0}T(F^{\scriptscriptstyle(0)}_{\scriptscriptstyle XY},h))={\cal{S}}_{\scriptscriptstyle Y|X}$.

For a function $h(x)$, we use $\partial h(X)/\partial X$ to denote $\partial h(x)/\partial x$ evaluated at $x=X$. The OPG aims to estimate central mean subspace ${\cal{S}}_{\scriptscriptstyle E[\kappa(Y,y)|X]}$ by $E\left[\frac{\partial E(\kappa(Y,y)|X)}{\partial X}\frac{\partial E(\kappa(Y,y)|X)}{\partial X^{\mbox{\tiny{\sf T}}}}\right]$ where the gradient $\partial E(\kappa(Y,y)|X)/{\partial X}$ is estimated by local linear approximation as follows. Let $K_{0}:\mathbb{R}\to[0,\infty)$ be a kernel function as used in kernel estimation. For any $v\in\mathbb{R}^{p}$ and bandwidth $h>0$, let $K_{h}(v)=h^{-p}K_{0}(\|v\|/h)$. At the population level, for fixed $x\in\Omega_{X}$ and $y\in\Omega_{Y}$, we minimize the objective function

over all $a\in\mathbb{R}$ and $b\in\mathbb{R}^{d_{0}}$. The minimizer depends on $x,y$ and we write it as $(a_{h}(x,y)$, $b_{h}(x,y))$. The ensemble tuned candidate matrix for estimating the central mean subspace ${\cal{S}}_{\scriptscriptstyle E[\kappa(Y,y)|X]}$ is $T_{\scriptscriptstyle 0}(F_{\scriptscriptstyle XY},\kappa(\cdot,y),h)=E[b_{\scriptscriptstyle h}(X,y)b_{\scriptscriptstyle h}(X,y)^{\mbox{\tiny{\sf T}}}]$ and the tuned Fréchet candidate matrix is $T(F_{\scriptscriptstyle XY},h)=E[b_{\scriptscriptstyle h}(X,Y)b_{\scriptscriptstyle h}(X,Y)^{\mbox{\tiny{\sf T}}}]$.

At the sample level, we minimize, for each $j,k=1,\dots,n$, the empirical objective function

over $a_{jk}\in\mathbb{R}$ and $b_{jk}\in\mathbb{R}^{p}$, where $w_{h}(X_{i},X_{j})=K_{h}(X_{i}-X_{j})/\sum_{l=1}^{n}K_{h}(X_{l}-X_{j}).$ Following Xia et al. (2002), we take the bandwidth to be $h=c_{0}n^{1/(p_{0}+6)}$ where $p_{0}=\max\{p,3\}$ and $c_{0}=2.34$, which is slightly larger than the optimal $n^{-1/(p+4)}$ in terms of the mean integrated squared errors. As proposed in Li (2018, Lemma 11.6), instead of solving $b_{jk}$ from (5) $n^{2}$ times, we solve multivariate weighted least squares to obtain $b_{j1},\dots,b_{jn}$ simultaneously. Computation details are given in Section S.3 of the Supplementary Material.

We can further enhance the performance of FOPG by projecting the original predictors onto the directions produced by the FOPG to re-estimate $\mathcal{S}_{Y|X}$. Specifically, after the first round of FOPG, we form the matrix $\hat{B}=(\hat{v}_{1},\dots,\hat{v}_{d})$ and replace the kernel $K_{h}(X_{j}-X_{i})$ by $K_{h}(\hat{B}^{\mbox{\tiny{\sf T}}}(X_{j}-X_{i}))$ with an updated bandwidth $h$, and complete the next round of iteration, which leads to an updated $\hat{B}$. We then iterate this process until convergence. In this way, we reduce the dimension of the kernel from $p$ to $d_{0}$ and mitigate the “curse of dimensionality”. To avoid confusion, We call this refined version of Fréchet OPG as FOPG in the following. The algorithms for FOPG are summarized as Algorithm 2.