基于投影重采样的多元非线性降维

On a Projective Resampling Method for Nonlinear Dimension Reduction with Multivariate Responses

摘要: 本文研究响应变量和预测因子均为向量的充分降维问题。投影重采样方法的核心思想是将多元响应投影到随机采样的方向上，以获取标量响应的样本，并反复应用单变量响应的降维方法来解决问题。该方法已被证明对多元线性降维有效。本文将投影重采样方法推广到非线性情境，并通过核映射提出了四种新的估计方法。研究结果表明，新方法具有优良的性质，并能在温和条件下完整恢复降维空间。最后，通过数值模拟和真实数据集分析验证了所提方法的有效性和可行性。

Abstract: This paper addresses the problem of sufficient dimension reduction where both the response and predictor are vectors. The core idea of projective resampling method is to project the multivariate responses along randomly sampled directions to obtain samples of scalar-valued responses. A univariate-response dimension reduction method is then applied repeatedly for solving the problem. This has proven effective for multivariate linear dimension reduction. In this paper, we extend the projective resampling method to nonlinear scenarios and use the mapping induced by kernels to develop four novel estimation methods. The research results suggest that the new methods exhibit excellent properties and ensure full recovery of the dimension reduction space under mild conditions. Finally, we validate the effectiveness and feasibility of the proposed methods through numerical simulations and real data analysis.

1. 引言

在经典情形，充分降维(Sufficient Dimension Reduction, SDR)考虑在给定$X\in {ℝ}^p$ 的条件下，$Y\in {ℝ}^1$ 的分布仅通过X的一组线性组合来依赖X [1]。如果存在一个矩阵$\beta \in {ℝ}^{p\times d}$ 且$d<p$ ，使得$Y\perp \text{}\text{}\text{}\text{}\perp X|{\beta }^{\text{T}}X$ ，其中$\perp \text{}\text{}\text{}\text{}\perp$ 表示独立性，那么由$\beta$ 的列向量张成的子空间表示中心子空间$S\left(Y|X\right)=span\left(\beta \right)$ ，d称为中心子空间的维数。参考切片逆回归(Sliced Inverse Regression, SIR) [1]，切片平均方差估计(Sliced Average Variance Estimate, SAVE) [2]，轮廓回归(Contour Regression, CR) [3]和方向回归(Directional Regression, DR) [4]。

近年来，一个重要的发展是将上述方法推广到多元响应变量情形。Setodji和Cook (2004) [5]提出了对多元Y进行高效分组的策略，并直接估计整个降维空间。Saracco (2005) [6]通过组合边际降维子空间来估计中心子空间，并且发展了估计的渐近分布。然而，这些方法受限于维度灾难，无法完全恢复中心子空间。Li等(2008) [7]提出的投影重采样方法(Projective Resampling Method, PR)克服了这两大难题，解决了响应变量是多维的问题。

随着数据生成和获取技术的快速发展，线性形式的${\beta }^{\text{T}}X$ 已经无法满足我们的需求，例如数据集中在流形时，线性方法通常无法寻找这种非线性结构。幸运的是，可以将数据嵌入希尔伯特空间的非线性映射，并且在这种空间中的线性方法可以捕捉低维结构。此时，可以通过核方法将数据映射到由核函数诱导的特征空间来解决问题。如果在希尔伯特空间中通过内积能够计算出低维结构的投影，就可以使用核技巧(kernel trick) [8]来获得简单且高效的估计方法。假设将p维数据通过特征映射$\varphi :X\to \mathscr{H}$ 投影到希尔伯特空间$\mathscr{H}$ 中，${\beta }^{\text{T}}X$ 可以表示为$\langle \beta ,\varphi \left(X\right)\rangle$ 。估计标量Y的降维方法可参考Wu (2008) [9]提出的核切片逆回归(Kernel Sliced Inverse Regression, KSIR)，Shotaro (2001) [10]、Bach和Jordan (2002) [11]、Fukumizu等(2007) [12]提出的核典型相关分析(Kernel Canonical Correlation Analysis, KCCA)，以及Lee等人(2013) [13]提出的广义切片平均方差估计(Generalized Sliced Average Variance Estimator, GSAVE)。

然而，当响应Y是一个向量时，还没有好的估计方法。Li等(2008) [7]提出了投影重采样方法来处理多元线性降维，并验证了其具有良好的性质。然而，该方法未考虑非线性情形，因此在处理多元非线性降维问题时仍然存在挑战。受此启发，我们将投影重采样方法推广至非线性情形，使用核方法来估计$\langle \beta ,\varphi \left(X\right)\rangle$ ，并且提出了四种新的降维方法。本文着重于研究这些方法的性质及适用条件。

本文的组织结构如下：第二节介绍了Mercer核，并推导了新的关系表达式。进一步，我们得到了四种新的非线性投影重采样方法。第三节给出了所提方法的理论性质。第四节通过模拟和实际数据集验证了这些方法的性质。最后，第五节对本文进行了总结。定理的证明见附录。

2. 非线性投影重采样方法

2.1. Mercer核和非线性降维

给定预测变量$X\in \mathcal{X}\subseteq {ℝ}^p$ ，Mercer核[14]是一个连续且半正定的函数$k\left(\cdot ,\cdot \right):\mathcal{X}\times \mathcal{X}\to ℝ$ ，其具有如下谱分解形式：

$k\left(x,z\right)={\displaystyle \sum _j{\lambda }_j{\varphi }_j\left(x\right){\varphi }_j\left(z\right)}$ ，

其中${\varphi }_j$ 为特征函数，${\lambda }_j$ 为相应的非负、非递增的特征值。Mercer核的一个重要性质是，每个核函数k唯一对应一个再生核希尔伯特空间(reproducing kernel Hilbert space, RKHS)，如下所示：

$\mathscr{H}=\left\{f|f\left(x\right)={\displaystyle \sum _{j\in \Lambda }{a}_j{\varphi }_j\left(x\right)\text{ with }{\displaystyle \sum _{j\in \Lambda }\frac{a_j^2}{{\lambda }_j}<\infty }}\right\}$ ，

其中集合$\Lambda :=\left\{j:{\lambda }_j>0\right\}$ 的基数即为再生核希尔伯特空间的维度，该维度可能是无穷的[14]。

给定一个Mercer核，存在一个唯一的映射$\varphi$ 从$\mathcal{X}$ 到由该核的特征值和特征函数定义的希尔伯特空间。该映射的形式如下：

$\varphi \left(x\right)=\left(\sqrt{{\lambda }_1}{\varphi }_1\left(x\right),\sqrt{{\lambda }_2}{\varphi }_2\left(x\right),\cdots ,\sqrt{{\lambda }_{\left|\Lambda \right|}}{\varphi }_{\left|\Lambda \right|}\left(x\right)\right)$ 。

$\varphi \left(Y\right)$ 由该映射所诱导的希尔伯特空间具有标准内积$k\left(x,z\right)=\langle \varphi \left(x\right),\varphi \left(z\right)\rangle$ ，并且该空间与再生核希尔伯特空间是同构的，我们将两个希尔伯特空间都记作$\mathscr{H}$ [15]。当核k是无限维的情况下，$\varphi :\mathcal{X}\to l_2$ 。此外，在核方法中，表示定理的关键作用在于，将可能在无限维希尔伯特空间$\mathscr{H}$ 上的优化问题简化为求解有限维的优化问题。

设Z是$\mathscr{H}$ 中的随机元素，满足$E\Vert Z\Vert <\infty$ ，其中$\Vert \cdot \Vert$ 表示由内积$\langle \cdot ,\cdot \rangle$ 诱导的$\mathscr{H}$ 中的范数。期望$E\left(Z\right)$ 被定义为$\mathscr{H}$ 中的一个元素，对于所有$a\in \mathscr{H}$ ，满足$\langle a,E\left(Z\right)\rangle =E\langle a,Z\rangle$ 。若$E\Vert Z\Vert <\infty$ ，则Z的协方差算子被定义为$E\left[\left(Z-EZ\right)\otimes \left(Z-EZ\right)\right]$ ，其中$\left(z\otimes b\right)f=\langle b,f\rangle a$ ，对于$\forall f\in \mathscr{H}$ 。随机变量$X\in \mathcal{X}$ 诱导了再生核希尔伯特空间中的一个随机元素$\varphi \left(X\right)$ 。不失一般性，假设$E\left(\varphi \left(X\right)\right)=0$ ，并且有界性也意味着协方差算子${\Sigma }_{XX}=E\left[\varphi \left(X\right)\otimes \varphi \left(X\right)\right]$ 是紧的，并且具有谱分解：

${\Sigma }_{XX}={\displaystyle \sum _{i=1}^{\infty }{w}_i{e}_i\otimes e_i}$ ，

其中$w_i$ 和$e_i\in \mathscr{H}$ 分别为特征值和特征函数。类似地，随机变量$Y\in \mathcal{X}$ 诱导了再生核希尔伯特空间中的一个随机元素。因此，有${\sum }_{YY}=E\left[\varphi \left(Y\right)\otimes \varphi \left(Y\right)\right]$ ，${\sum }_{YX}=E\left[\varphi \left(Y\right)\otimes \varphi \left(X\right)\right]$ 和${\sum }_{XY}=E\left[\varphi \left(X\right)\otimes \varphi \left(Y\right)\right]$ 。

我们假设Y和X之间的关系具有如下模型：

$Y=F\left(\langle {\beta }_1,\varphi \left(X\right)\rangle ,\cdots ,\langle {\beta }_d,\varphi \left(X\right)\rangle ,ϵ\right)$

其中${\beta }_j\in \mathscr{H}$ ，且误差项$ϵ$ 的分布独立于X。该模型意味着响应变量Y仅通过一个d维的总结统计量$\mathcal{S}\left(Y|X\right)=\left(\langle {\beta }_1,\varphi \left(X\right)\rangle ,\cdots ,\langle {\beta }_d,\varphi \left(X\right)\rangle \right)$ 来依赖于X。虽然$S\left(Y|X\right)$ 在$\mathscr{H}$ 中是一个线性总结统计量，但它提取了原始预测变量X空间中的非线性特征。

Li (1991) [1]使用矩阵${\left[var\left(X\right)\right]}^{-1}var\left[E\left(X|Y\right)\right]$ 的特征向量来估计线性充分降维中的$\beta$ ；而非线性特征则通过算子$T={\left[var\left(\varphi \left(X\right)\right)\right]}^{-1}var\left[E\left(\varphi \left(X\right)|Y\right)\right]$ 的特征函数来恢复。例如，Wu (2008) [9]提出的核切片逆回归(KSIR)方法首先将Y的样本分成若干个切片，然后取$\varphi \left(X\right)$ 的切片均值来形成候选矩阵，并在此矩阵上进行谱分解。Lee 等(2016) [16]则考虑回归算子来处理给定Y时X的函数的条件期望。即，对于$\forall f\in \mathscr{H},{\Sigma }_{YY}^{-1}{\Sigma }_{YX}f=E\left[f\left(X\right)|Y\right]+常数$ 。因此，Li (2018) [17]给出了T的两种形式：

$T_1={\sum }_{XX}^{-1}{\sum }_{XY}{\sum }_{YX}{\sum }_{XX}^{-1},T_2={\sum }_{XX}^{-1}{\sum }_{XY}{\sum }_{YY}^{-1}{\sum }_{YX}{\sum }_{XX}^{-1},$

这些方法被称为广义切片逆回归(Generalized Sliced Inverse Regression, GSIR)。虽然KSIR在数值上与GSIR相似，但KSIR背后的理论基础却有很大不同：它假设$E\left[\varphi \left(X\right)|{\beta }^{\text{T}}\varphi \left(X\right)\right]$ 对于$\varphi \left(X\right)$ 是线性的，而这一假设类似于线性充分降维(SDR)的假设。而GSIR并不需要这一假设。需要注意的是，中心化的Gram矩阵K是由核函数$k\left(\cdot ,\cdot \right)$ 定义的，其中$K_{ij}=\langle \varphi \left(x_i\right),\varphi \left(x_j\right)\rangle =k\left(x_i,x_j\right)$ 。此外，我们可以为协方差算子中添加一个正则化项[17]，以防止其不可逆。那么，

$T_{n,1}={\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}{K}_X{K}_Y{K}_X{\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1},$ (1)

$T_{n,2}={\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}{K}_X{K}_Y{\left(K_Y+{\epsilon }_Y{\lambda }_{\max }\left(K_Y\right)I_n\right)}^{-1}{K}_X{\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}$ .(2)

其中${\epsilon }_X$ 和${\epsilon }_Y$ 是通过交叉验证算法[13]获得的参数，${\lambda }_{\max }(\cdots )$ 表示一个矩阵最大的特征值。

根据矩阵的选择，GSIR可以分为两种形式，称为GSIR1和GSIR2，分别对应于表达式(1)和(2)的基矩阵。类似地，其他非线性方法也是在寻找不同的矩阵T来获得解。下一节将PR与这些核方法结合起来，发展出新的降维方法。

2.2. 投影重采样方法

回到多元非线性问题中，设X和Y的维度分别为p和q。对于$\forall u\in {ℝ}^q$ ，考虑n个投影样本 $\left(X_1,u^{\text{T}}{Y}_1\right),\cdots ,\left(X_n,u^{\text{T}}{Y}_n\right)$ 。假设给定一个用于单响应变量的基本方法，例如前文描述的GSIR1。则总体水平的矩阵T由$T=var\left\{E\left[\varphi \left(X\right)|u^{\text{T}}Y\right]\right\}$ 给出。那么，对于每一个$u\in {ℝ}^q$ ，可以获得一个半正定矩阵$T_n\left(u\right)$ ，它用于估计$\mathcal{S}\left(u^{\text{T}}Y|X\right)$ 。最后，这些矩阵取平均值以估计$S\left(Y|X\right)$ 。

当然，实际操作中不可能对${ℝ}^q$ 中每一个u都执行$u^{\text{T}}Y$ 对X的降维。Li等(2008) [7]提出，我们可以通过适当的分布生成一个有限集u。研究表明，只要这个集合相对于样本量足够大，基于该有限集u所得到的估计与对${ℝ}^q$ 中每一个u进行降维所得到的估计是渐近等价的。接下来我们给出对于给定d (其估计将在后续讨论)和给定的基估计矩阵$T_n\left(u\right)$ 的具体算法：

算法1：PR-GSIR1算法

1) 选择蒙特卡洛样本大小$m_n$ 满足$m_n$ 比n趋于无穷的速度快，比如$m_n=n\log \left(n\right)$ 或$m_n=n^{3/2}$ 。取$G_1,\cdots ,G_{m_n}$ 独立同分布于$N\left(0,I_q\right)$ ，那么得到$m_n$ 个投影$u_j=G_j/\Vert G_j\Vert ,j=1,\cdots ,m_n$ 。 2) 对每个$u_j,j=1,\cdots ,m_n$ ，对$\left(X_1,u_j^{\text{T}}{Y}_1\right),\cdots ,\left(X_n,u_j^{\text{T}}{Y}_n\right)$ 使用GSIR1来计算$T_n\left(u_j\right)$ ，即： $T_n\left(u_j\right)={\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}{K}_X{K}_{u_j^{T}Y}{K}_X{\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}$ ， 这个过程得到了$m_n$ 个矩阵。 3) 计算$T_{n,m_n}=m_n^{-1}{\displaystyle \sum _{j=1}^{m_n}{T}_n}\left(u_j\right)$ ，对$T_{n,m_n}$ 执行谱分解得到d个主特征向量${\widehat{w}}_1,\cdots ,{\widehat{w}}_d$ ，那么总结统计量$\mathcal{S}\left(Y|X\right)$ 的估计为${\left(K_X+{\epsilon }_X{\lambda }_{\text{max}}\left(K_X\right)I_n\right)}^{-1}{\widehat{w}}_i,i=1,\cdots ,d$。

上述方法是投影重采样与广义切片逆回归的第一种矩阵结合起来得到的新方法，我们把它命名为PR-GSIR1。需要注意的是，参数${\epsilon }_X$ 是通过Lee等(2013) [13]提出的交叉验证算法得到的，在最后一步对求逆过程进行类似于岭回归的正则化是为了提高精确度。类似地，算法PR-GSIR2可以通过将算法1的第2步替换为$T_n\left(u_j\right)$ 来得到：

$T_n\left(u_j\right)={\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}{K}_X{K}_{u_j^{\text{T}}Y}{\left(K_{u_j^{\text{T}}Y}+{\epsilon }_{u_j^{\text{T}}Y}{\lambda }_{\max }\left(K_{u_j^{\text{T}}Y}\right)I_n\right)}^{-1}{K}_X{\left(K_X+{\epsilon }_X{\lambda }_{\max }\left(K_X\right)I_n\right)}^{-1}$ 。

算法与第一种方法类似，我们不做过多解释。为了表明投影重采样适用于非线性充分降维，我们也将投影重采样与Akaho (2001) [10]，Bach和Jordan (2002) [11]和Fukumizu等(2007) [12]发展的核典型相关分析结合起来，可以得到第三种方法(PR-KCCA)通过使用$T_n\left(u_j\right)$ ：

$T_n\left(u_j\right)={\left(K_X+{\epsilon }_X{\lambda }_{\text{max}}\left(K_X\right)I_n\right)}^{-1}{K}_X{K}_{u_j^{\text{T}}Y}{\left(K_{u_j^{\text{T}}Y}+{\epsilon }_{u_j^{\text{T}}Y}{\lambda }_{\text{max}}\left(K_{u_j^{\text{T}}Y}\right)I_n\right)}^{-2}{K}_{u_j^{\text{T}}Y}{K}_X{\left(K_X+{\epsilon }_X{\lambda }_{\text{max}}\left(K_X\right)I_n\right)}^{-1}$ 。

但是，上述提到的三种方法都只适合于提取响应变量的条件均值，因此我们也引入使用广义切片平均方差估计得到的新方法，其对于提取响应的条件方差效果更好，我们会在第4节模拟部分验证。为了引入广义切片平均方差估计的算法，还需要建立一些关系。$1_n$ 是元素全是1的n维向量，矩阵$Q=I_n-1_n{1}_n^{\text{T}}/n$ ，$L_X={\left(1_n,K_X\right)}^{\text{T}}$ ，$L_Y={\left(1_n,K_Y\right)}^{\text{T}}$ 。那么可以计算：$C=L_Y^{\text{T}}{\left(L_Y{L}_Y^{\text{T}}+{\epsilon }_Y{I}_{n+1}\right)}^{-1/2}{L}_Y$ ，令$C_i$ 为C的第i列，计算${\Lambda }_i=\text{diag}\left(C_i\right)-C_i{C}_i^{\text{T}}$ ，以及${\Gamma }_i=Q/n-{\Lambda }_i$ ，对于$i=1,\cdots ,n$ 。第四种降维方法(PR-GSAVE)的算法如下：

算法2：PR-GSAVE算法

1) 与算法1中第1步一致，得到$m_n$ 个投影。 2) 对每个$u_j,j=1,\cdots ,m_n$ ，对$\left(X_1,u_j^{\text{T}}{Y}_1\right),\cdots ,\left(X_n,u_j^{\text{T}}{Y}_n\right)$ 使用GSAVE来计算$T_n\left(u_j\right)$ ，即： $T_n\left(u_j\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(L_{X}Q{L}_X^{\text{T}}+{\epsilon }_X{\lambda }_{\max }\left(L_{X}Q{L}_X^{\text{T}}\right)I_{n+1}\right)}^{-1/2}}\cdot$ $L_{X}Q{\Gamma }_{i}Q{\Gamma }_{i}Q{L}_X^{\text{T}}{\left(L_{X}Q{L}_X^{\text{T}}+{\epsilon }_X{\lambda }_{\max }\left(L_{X}Q{L}_X^{\text{T}}\right)I_{n+1}\right)}^{-1/2}$ 这个过程得到了$m_n$ 个矩阵。 3) 计算$T_{n,m_n}=m_n^{-1}{\displaystyle \sum _{j=1}^{m_n}{T}_n}\left(u_j\right)$ ，对$T_{n,m_n}$ 执行谱分解得到d个主特征向量${\widehat{w}}_1,\cdots ,{\widehat{w}}_d$ ，那么总结统计量$\mathcal{S}\left(Y|X\right)$ 的估计为${\left(L_{X}Q{L}_X^{\text{T}}+{\epsilon }_X{\lambda }_{\max }\left(L_{X}Q{L}_X^{\text{T}}\right)I_{n+1}\right)}^{-1/2}{\widehat{w}}_i,i=1,\cdots ,d$。

对上述PR-GSAVE中的投影体现在从$K_Y$ 到$K_{u_j^{\text{T}}Y}$ 的改变。这4种方法的性质会在第4节的模拟部分验证。

在算法过程中，我们使用Zhu等(2006) [18]介绍的贝叶斯信息准则来确定中心子空间的维数d。其他贝叶斯信息准则可参考，在Wang和Yin (2008) [19]，Li等(2010) [20]以及Guo等(2015) [21]。下一节研究非线性投影重采样方法的性质。

3. 非线性投影重采样的性质

3.1. 费希尔一致性

对于每个$u\in {ℝ}^q$ ，令$T\left(u\right)$ 是一个$p\times p$ 的半正定矩阵满足$span\left[T\left(u\right)\right]=\mathcal{S}\left(u^{\text{T}}Y|X\right)$ ，其中$T\left(u\right)$ 是非线性降维方法的总体水平矩阵。例如，对于广义切片逆回归，$T\left(u\right)=var\left\{E\left[\varphi \left(X\right)|u^{T}Y\right]\right\}$ 。那么可以得到以下定理：

定理3.1：假设u是均匀分布在${\mathbb{S}}^q$ 上的随机向量，对于每个固定的$u\in {ℝ}^q$ ，$T\left(u\right)\ge 0$ 都有 $\text{span}\left\{T\left(u\right)\right\}=\mathcal{S}\left(u^{\text{T}}Y|X\right)$ 。那么$span\left\{E\left[T\left(u\right)\right]\right\}=\mathcal{S}\left(Y|X\right)$ 。

对于定理3.1需要说明的是，虽然使用均匀分布得到投影，但是分布的选取对定理证明并没有起到关键的作用。所以我们可以使用其他分布，例如算法中使用的正态分布。并且定理3.1表明，如果我们想估计$S\left(Y|X\right)$ ，就可以用对$E\left[T\left(u\right)\right]$ 的估计来代替。这表明了算法的可行性及合理性。而期望$E\left[T\left(u\right)\right]$ 可以用蒙特卡洛近似来求解积分

${\displaystyle {\int }_{{\mathbb{S}}^q}{T}_n}\left(u\right)\lambda \left(\text{d}u\right)$ ，

来近似估计，其中$\lambda$ 是在${ℝ}^q$ 上的均匀分布。对于$T\left(u\right)\ge 0$ 的原因在于确保$E\left[T\left(u\right)\right]$ 的贡献是可测的，即：对${ℝ}^q$ 中的两个不同的点$u_1$ 和$u_2$ ，那么$T\left(u_1\right)$ 和$T\left(u_2\right)$ 对积分$E\left[T\left(u\right)\right]$ 的贡献不会被抵消。

3.2. $\sqrt{n}$ 一致性和蒙特卡洛效应

很明显，只要$T_n\left(u\right)$ 对每个u都是$\sqrt{n}$ 一致的并且满足一些均匀分布条件，那么执行蒙特卡洛过程无穷多次，可以得到估计量$T_{n,m_n}$ 将会收敛于$T_{n,\infty }\equiv E\left(T_n\left(u\right)|X_1,Y_1,\cdots ,X_n,Y_n\right)$ ，那么$T_{n,m_n}$ 也是$\sqrt{n}$ 一致的。假设$T_n\left(u\right)$ 有以下展开式：

$T_n\left(u\right)=T\left(u\right)+E_n\Phi \left(X,Y,u\right)+R_n\left(u\right)$ ，(3)

其中，对每个u，$\Phi \left(X,Y,u\right)$ 这一项都是均值为0且有限二阶距。余项$R_n\left(u\right)$ 满足

$\sup \Vert R_n\left(u\right)\Vert =o_p\left(n^{-1/2}\right)$ ，(4)

其中$\Vert \cdot \Vert$ 表示矩阵的F范数。对于条件(3)和(4)的细节和应用可参考Fernholz (2012) [22]。下一个定理表明$T_{n,m_n}$ 的$\sqrt{n}$ 一致性。

定理3.2：如果条件(3)和(4)均成立，$n=O\left(m_n\right)$ ，那么$T_{n,m_n}=E\left[T\left(u\right)\right]+O_p\left(n^{-1/2}\right)$ 。

该定理(参见Li等(2008) [7]——虽然仅展示了线性情况，但将其扩展到流形上的推导是直接的)告诉我们，为了实现对联合统计量$S\left(Y|X\right)$ 的$\sqrt{n}$ 一致性，只需对$O\left(n\right)$ 个随机采样u的边际统计量$\mathcal{S}\left(u^{\text{T}}Y|X\right)$ 进行估计即可。接下来，我们研究$T_{n,m_n}$ 和$T_{n,\infty }$ 渐进相等时，$m_n\to \infty$ 的速度，这意味着蒙特卡洛效应会消失。令$\text{vec}\left(A\right)$ 表示对矩阵A的向量化，令B表示矩阵$B=\mathrm{var}\left\{E\left[vec\left(\Phi \left(X,Y,u\right)\right)|X,Y\right]\right\}$ ，$T=E\left[T\left(u\right)\right]$ 。那么可以得到如下定理：

定理3.3：如果条件(3)和(4)均成立，那么当$n\to \infty$ 时，

$\sqrt{n}\left[\text{vec}\left(T_{n,\infty }\right)-\text{vec}\left(T\right)\right]\stackrel{D}{\to }N\left(0,B\right)$。

并且也满足$\underset{n\to \infty }{\lim }\left(\frac{n}{m_n}\right)=0$ ，那么当$n\to \infty$ 时，

$\sqrt{n}\left[\text{vec}\left(T_{n,m_n}\right)-\text{vec}\left(T\right)\right]\stackrel{D}{\to }N\left(0,B\right)$。

这个定理(参见Li等(2008) [7]——虽然仅展示了线性情况，但将其扩展到流形上的推导是直接的)很明显地说明，如果$m_n$ 趋于无穷的速度比n快，那么$\text{vec}\left(T_{n,m_n}\right)$ 和$\text{vec}\left(T_{n,\infty }\right)$ 有相同的渐进正态分布。这就意味着使用$\text{vec}\left(T_{n,m_n}\right)$ 和$\text{vec}\left(T_{n,\infty }\right)$ 来估计T是渐进一致的，也验证了3.1节中当$n\to \infty$ 时用蒙特卡洛近似求积分的误差可以忽略。下一节我们通过模拟来检验新方法的性质及适用条件。

4. 模拟比较

4.1. 模拟数据对比

由于2.2节所述的原因，我们在充分预测因子出现在条件均值时比较PR-GSIR1与PR-GSIR2和PR-KCCA，同时在充分预测因子出现在条件方差比较PR-GSAVE与PR-GSIR1、PR-GSIR2和PR-KCCA。使用如下的可加高斯核：

$k\left(x,z\right)={\displaystyle \sum _{j=1}^d\text{exp}\left(-\frac{{\left(x_j-z_j\right)}^2}{2{\sigma }^2}\right)}$ ，

其中$\sigma =2$ 。为了验证我们的方法确实适用于非线性降维问题，我们还采用了Li等(2008) [7]提出的PR方法进行比较。由于该方法以SIR作为其基础算法，将其记作PR-SIR。

首先选择充分预测因子为条件均值，考虑以下三个模型：

模型I：$Y_1=\sqrt{X_1^2+X_2^2}\log \left(\sqrt{X_1^2+X_2^2}\right)+{ϵ}_1;Y_2={ϵ}_2$

模型II：$Y_1=\frac{X_1}{1+{\text{e}}^{X_2}}+{ϵ}_1;Y_2={ϵ}_2$

模型III：$Y_1=\sin \left(\frac{\pi \left(X_1+X_2\right)}{10}\right)+{ϵ}_1;Y_2={ϵ}_2$

其中$X~N\left(0,I_p\right),X\perp \text{}\text{}\text{}\text{}\perp ϵ,ϵ~N\left(0,\Delta \right)$ ，Δ满足

$\Delta =\left(\begin{array}{cc}1& 0.5\\ 0.5& 1\end{array}\right)$ 。

X的维数是p，Y的维数q固定为3，生成的样本数为n，选择蒙特卡洛样本大小$m_n=2000$ 。通过估计的充分预测因子与多维响应的接近程度来评估算法结果，由于只关注预测因子的单调函数，我们使用Spearman相关性作为接近程度的衡量标准。对于样本数$n=100,200$ 和X的维数$p=10,20$ 的四种结合，分别产生模型I、模型II和模型III的100个样本。对于PR-SIR，切片的数量是10。然后求出相应的Spearman相关系数的平均值。结果见表1。

Table 1. Average absolute Spearman correlation coefficients between estimated predictors and each Y_i for PR-SIR, PR-GSIR1, PR-GSIR2, and PR-KCCA across different sample sizes and predictor dimensions for Models I, II, and III (100 Repetitions)

表1. PR-SIR、PR-GSIR1、PR-GSIR2和PR-KCCA在不同样本数和不同预测因子维数，对于模型I、II和III，估计预测因子与每个Y_i的Spearman相关系数绝对值在100重复得到的平均值

(n, p)	模型	PR-SIR		PR-GSIR1		PR-GSIR2		PR-KCCA
		Y1	Y2	Y1	Y2	Y1	Y2	Y1	Y2
(100, 10)	I	0.24	0.16	0.92	0.71	0.82	0.39	0.83	0.38
	II	0.47	0.10	0.87	0.71	0.82	0.54	0.79	0.57
	III	0.42	0.13	0.83	0.74	0.79	0.60	0.79	0.61
(100, 20)	I	0.26	0.16	0.94	0.74	0.91	0.19	0.91	0.57
	II	0.47	0.18	0.91	0.76	0.87	0.72	0.88	0.70
	III	0.44	0.22	0.88	0.77	0.86	0.74	0.85	0.74

续表

(200, 10)	I	0.22	0.12	0.89	0.64	0.80	0.29	0.80	0.32
	II	0.45	0.07	0.85	0.67	0.74	0.44	0.76	0.48
	III	0.39	0.07	0.81	0.70	0.73	0.51	0.71	0.52
(200, 20)	I	0.22	0.15	0.93	0.75	0.88	0.42	0.90	0.48
	II	0.49	0.10	0.89	0.78	0.87	0.65	0.87	0.62
	III	0.44	0.11	0.86	0.76	0.85	0.71	0.85	0.71

从上表可以看出三种非线性方法都有着满意的结果，且PR-GSIR1的结果比另两种方法更佳。然而，PR-SIR效果较差，这验证了本文提出的方法确实适合于处理非线性降维问题。并且可以看出，这三种非线性方法均不受维数灾难的影响。

接下来，我们比较当预测因子只影响方差时，PR-GSAVE、PR-GSIR1、PR-GSIR2和PR-KCCA这四种方法的性质。前三个模型中，我们计算多元响应的每个分量与估计的充分预测因子之间的相关性，接下来我们提出一个新的标准来综合考虑多元响应与预测因子的相关程度。因为$Y\perp \text{}\text{}\text{}\text{}\perp X|f\left(X\right)$ ，对于$f\left(X\right)$ 生成的可积函数类空间的正交补中的元素$g\left(X\right)$ ，则有$Y\perp \text{}\text{}\text{}\text{}\perp g\left(X\right)$ ，从而对于$\forall u\in {ℝ}^q$ ，$u^{\text{T}}Y\perp \text{}\text{}\text{}\text{}\perp g\left(X\right)$ ，所以$u^{\text{T}}Y$ 与$g\left(X\right)$ 的相关系数为0，因此一个好的估计应该使得$u^{\text{T}}Y$ 与$g\left(X\right)$ 的相关系数尽可能小。

Figure 1. Boxplots of the absolute values of the correlation coefficients between $u^{\text{T}}Y$ and $g\left(X\right)$ obtained by the four methods in Model IV

图1. 四种方法在模型IV得到的$u^{\text{T}}Y$ 与$g\left(X\right)$ 的相关系数绝对值的箱线图

Figure 2. Boxplots of the absolute values of the correlation coefficients between $u^{\text{T}}Y$ and $g\left(X\right)$ obtained by the four methods in Model V

图2. 四种方法在模型V得到的$u^{\text{T}}Y$ 与$g\left(X\right)$ 的相关系数绝对值的箱线图

考虑以下两个模型：

模型IV：$Y_1=X_1{ϵ}_1;Y_2=X_2{ϵ}_2;Y_3={ϵ}_3$

模型V：$Y_1=\frac{X_1}{1+{\text{e}}^{X_2}}{ϵ}_1;Y_2=\frac{X_3}{1+{\text{e}}^{X_4}}{ϵ}_2;Y_3={ϵ}_3$

其中$X~N\left(0,I_{10}\right),X\perp \text{}\text{}\text{}\text{}\perp ϵ,ϵ~N\left(0,\Delta \right)$ ，Δ满足

$\Delta =\left(\begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 0\\ 0& 0& 0.5\end{array}\right)$ 。

生成的样本数$n=100$ ，选择蒙特卡洛样本大小$m_n=2000$ 。结果见图1和图2。

图1，图2分别是模型IV，模型V产生来自$N\left(0,I_3\right)$ 三维正态随机数$u_1,\cdots ,u_{200}$ ，使用四种方法计算$u^{\text{T}}Y$ 与$g\left(X\right)$ 的复相关系数的箱线图，从图中可以看出，PR-GSAVE的效果是最好的，这验证了广义切片平均方差估计处理响应的条件方差效果更好。

4.2. 真实数据分析

实际数据集使用R软件程序包dr (下载网址：http://mirrors.xmu.edu.cn/CRAN/)中自带数据集ais，为澳大利亚体育研究所收集202名男女运动员身体各个指标的数据，选取瘦体重(lean body mass, LBM)、脂肪率(body fat percentage, Bfat)为2维响应变量，身高(height, Ht)、体重(weight, Wt)、红细胞数量(red blood cell count, RCC)、白细胞数量(white blood cell count, WCC)、红细胞比容(hematocrit, Hc)、血色素(hemoglobin, Hg)、血浆铁蛋白浓度(plasma ferritin concentration, Ferr)、体质指数(body mass index, BMI)、皮肤褶皱(skinfold thickness, SSF)为9维协变量，结果见图3。图3为来自$N\left(0,I_2\right)$ 三维正态随机数$u_1,\cdots ,u_{200}$ ，使用四种方法计算$u^{\text{T}}Y$ 与$g\left(X\right)$ 的复相关系数的箱线图，可以看出PR-GSAVE的效果相对其他三种方法较差，而PR-GSIR1的效果最好。

Figure 3. Boxplots of the absolute values of the correlation coefficients between $u^{\text{T}}Y$ and $g\left(X\right)$ obtained by the four methods on the dataset of the athlete physical indicators

图3. 四种方法在运动员身体指标数据集得到的$u^{\text{T}}Y$ 与$g\left(X\right)$ 的相关系数绝对值的箱线图

综上分析，可以得出初步结论。对于一般的多元非线性降维问题，PR-GSIR1的效果应该是最好的，并且这种方法也是最简单的。这与线性降维中情况一致，即切片逆回归是简单且使用最广泛的方法之一。对于预测因子出现在条件方差的情形中，PR-GSAVE的效果就相对较好了，这也与线性中情况一致。我们可以将PR与其他降维方法结合，寻找更优的降维方法。

5. 结论

投影重采样基于这样一个事实：一个随机向量在所有方向上的投影唯一确定了该随机向量的联合分布。这一方法能够简单而有效地将多元响应降维转换为单变量响应降维，从而避免维数灾难，提高精确度。我们将投影重采样这种用于处理多元响应的降维技术与一般非线性降维方法相结合，得到了四种新的降维方法。令人振奋的是，我们的方法只需要在蒙特卡洛样本大小比n趋于无穷的速度快，就能保证良好的性质。在模拟实验中，我们验证了新方法的有效性，结果显示这些新方法仍然保留了原降维方法的特点。例如，PR-GSAVE能够更好地处理预测因子为条件方差的情况，这类似于广义切片平均方差估计和切片平均方差估计。我们也期望将投影重采样与其他核方法相结合，以获得更多具有良好性质的降维方法，这将成为未来的研究方向。

基金项目

本文得到山西省自然科学基金(项目编号：RD2200002021)的支持。

附录：证明

定理3.1的证明：

由于T表示$E\left[T\left(u\right)\right]$ ，令P表示在$\mathcal{S}\left(Y|X\right)$ 上的投影函数，$P_T$ 是在$span\left(T\right)$ 上的投影函数。

我们首先表明$span\left(T\right)\subseteq \mathcal{S}\left(Y|X\right)$ 。根据统计量$\mathcal{S}\left(Y|X\right)$ 的定义，$X\perp \text{}\text{}\text{}\text{}\perp Y|\langle \beta ,\varphi \left(X\right)\rangle$ 。因此，对$\forall u\in {\mathbb{S}}^q$ ，有$X\perp \text{}\text{}\text{}\text{}\perp u^{\text{T}}Y|\langle \beta ,\varphi \left(X\right)\rangle$ 。这表明，对$\forall u\in {\mathbb{S}}^q$ ，都可以得到$\mathcal{S}\left(u^{\text{T}}Y|X\right)\subseteq \mathcal{S}\left(Y|X\right)$ ，或者等价地$span\left[T\left(u\right)\right]\subseteq \mathcal{S}\left(Y|X\right)$ 。令$w\perp \mathcal{S}\left(Y|X\right)$ ，那么对$\forall u\in {\mathbb{S}}^q$ 都有：$w\perp T\left(u\right)$ ，因此$\langle T,w\rangle =\langle E\left[T\left(u\right)\right],w\rangle =0$ ，从而$w\perp span\left(T\right)$ 。

接下来我们证明$\mathcal{S}\left(Y|X\right)\subseteq span\left(T\right)$ 。只需证明对$\forall u\in {ℝ}^q$

$E\left({\text{e}}^{i{u}^{\text{T}}Y}|\varphi \left(X\right)\right)=E\left({\text{e}}^{i{u}^{\text{T}}Y}|\langle P_T,\varphi \left(X\right)\rangle \right)$ 。

因为$X\perp \text{}\text{}\text{}\text{}\perp Y|\langle P_T,\varphi \left(X\right)\rangle$ 等价于$X\perp \text{}\text{}\text{}\text{}\perp |\left|u\right|\left|Y\right|\langle P_T,\varphi \left(X\right)\rangle$ ，对所有$\forall u\in {\mathbb{S}}^q$ 都足以证明上述表达式。并且由于特征函数是连续的，这足以证明对${\mathbb{S}}^q$ 的一些密子集中的u上式均成立。这就证明了定理3.1。

NOTES

^*第一作者。

^#通讯作者。