第38卷第2期
2021年2月统计研究Statistical Ｒesearch Vol．38，No.2Feb．2021
大数据下张量充分降维方法
及其应用研究*
马少沛孙庆慧武雅萱田茂再
内容提要：在大数据时代，金融学、基因组学和图像处理等领域产生了大量的张量数据。Zhong 等
（2015）提出了张量充分降维方法，并给出了处理二阶张量的序列迭代算法。鉴于高阶张量在实际生活
中的广泛应用，本文将Zhong 等（2015）的算法推广到高阶，以三阶张量为例，提出了两种不同的算法：结
构转换算法和结构保持算法。两种算法都能够在不同程度上保持张量原有结构信息，同时有效降低变
量维度和计算复杂度，避免协方差矩阵奇异的问题。将两种算法应用于人像彩图的分类识别，以二维和
三维点图等形式直观展现了算法分类结果。将本文的结构保持算法与K-
means 聚类方法、t-SNE 非线性降维方法、
多维主成分分析、多维判别分析和张量切片逆回归共五种方法进行对比，结果表明本文所提方法在分类精度方面有明显优势，因此在图像识别及相关应用领域具有广阔的发展前景。
关键词：张量；充分降维；迭代算法；图像识别
DOI ：10.19343/j．cnki．11－1302/c．2021.02.009
中图分类号：O212文献标识码：A 文章编号：1002－4565（2021）02－0114－21
*基金项目：中国人民大学科学研究基金（中央高校基本科研业务费专项资金资助）项目“高维复杂数据充分降维的稳健方法”
（20XHN106）。Ｒesearch on Tensor Sufficient Dimension Ｒeduction
Method and Its Application in Big Data
Ma Shaopei Sun Qinghui Wu Yaxuan Tian Maozai
Abstract ：In the era of big data ，a large amount of tensor data are generated in many fields such as
finance ，genomics and image processing．Zhong et al．（2015）propose tensor sufficient dimension reduction and
present a sequential iterative algorithm for second-order tensors．In view of the wide application of high-order
tensors in real life ，this paper extends the iterative algorithm of Zhong et al．（2015）to higher-
order tensors．Taking third-order tensors as an example ，two different algorithms are proposed ：structural transformation
algorithm and structural maintenance algorithm．The two algorithms can effectively reduce the variable dimension
and computational complexity while maintaining the original structural information of tensors ，avoiding the
singular problem of covariance matrix．The two algorithms are applied to the classification and recognition of
color facial images．Then 2-D and 3-D scatter plots are used to display the classification results of the
algorithms．Comparing the structural maintenance algorithm with the K-means clustering method ，t-
SNE nonlinear dimension reduction method ，multilinear principal component analysis ，multilinear discriminant
analysis and tensor SIＲ，the paper shows that the proposed method has obvious advantages in classification
accuracy ，and has broad development prospect in image recognition and related fields of application．
Key words ：Tensor ；Sufficient Dimension Ｒeduction ；Iterative Algorithm ；Image Ｒecognition
第38卷第2期马少沛等：大数据下张量充分降维方法及其应用研究·115·
一、引言
科学技术的飞速发展使人们拥有海量的数据来挖掘变量之间的结构性关系，并由此进行深入的数据分析和趋势预测。田茂再（2015）对大数据时代下统计学理论研究中的热点问题做了前瞻性的介绍。其中，高维数据下的维数灾难问题受到了社会的广泛关注。当预测变量的维度远远超过样本量时，往往会面临
函数逼近、参数估计、信息提取等方面的困难。因此，在保留重要信息的前提下，对预测变量进行降维处理非常有必要。
常用的降维方法有两类：变量选择和低维投影。变量选择的基本思想是从所有待选择的变量中，挑选出对响应变量影响最大的一个变量子集。低维投影的目的是到一组投影方向，使得原始的高维变量在该方向上的投影能够包含与响应变量有关的绝大部分信息，主要方法有主成分回归法、投影寻踪法和充分降维法。陈道铨等（2010）通过蒙特卡洛模拟对不同的投影子空间估计方法进行了比较和推广。现有的文献对基于向量的变量选择和降维问题已有较深入的研究。
然而，在大数据时代，仅仅依靠向量这种表示形式难以全面考虑原始数据各个维度的变化情况，并且容易丢失有效信息。作为向量的扩展，张量是一个多方向的数组，能够最大化地保留数据的原始结构信息，更好地适应了实际生产对高通量数据的储存和处理要求。随着社会科学、成像技术等领域的快速发展，张量的处理和应用也变得越来越重要。
以往文献对张量型数据的主要处理方法是对其进行向量化处理，然后应用基于向量的传统方法进行降维分析。但是简单的向量化操作极易破坏原始的数据结构，并使算法陷入维数灾难，带来复杂的计算问题。因此许多学者开始考虑直接对张量结构进行降维，Cheng等（2020）对现有的张量研究方法做了总结，主要分为无监督和有监督两类。
在无监督学习方法中，主成分分析（PCA）具有极高的应用价值，因此很多文献对主成分分析方法进行了不同形式的张量型推广。例如，对于二阶张量的主成分分析，Yang等（2004）提出了2－DPCA算法，对矩阵形式的输入变量进行线性变换以达到降维目的。为了处理更高阶的张量数据，Lu等（2008）提出了多维主成分分析（MPCA），在满足投影数据方差最大的要求下，可以将任意阶的张量投影到低维空间加以分析。在理论方面，Hung等（2012）研究了二阶张量的主成分所具有的渐近分布和渐近效率等理论性质。在实际数据应用中，为了进一步降低张量维度，Chen等（2014）利用两阶段PCA方法对低温电子显微镜（cryo-EM）图像数据进行了聚类分析。
作为一种有监督的降维方法，线性判别分析（LDA）方法在张量数据中也得到了极大的推广和发展。Ye等（2005）较早地提出了2DLDA方法，将矩阵的行列结构同时进行线性变换，从而得到低维投影矩阵，但存在迭代算法不收敛的问题。Tao等（2007）提出了一般张量判别分析（GTDA）方法，在散点差值准则下首次提出了具有收敛性质的算法，但其表现效果依赖于准则中调节参数的选择。Li等（2014）在此基础上改进后得到了DGTDA方法，避免了参数选择的问题，同时能够提高时间效率和分类精度。Tao等（2017）从张量秩保持的角度对判别分析进行了研究，同时指出LDA方法通常需要较大的样本量才能获得良好的模型表现，因此经常面临小样本问题。
而张量充分降维是监督学习框架下的一种新兴处理张量数据的方法，近年来得到学者们的广泛关注。该方法旨在从原有的张量结构中到一组低维解释变量，使其能够包含类别标签或回归响应中的绝大部分
信息。Li等（2010）率先提出了一种维数折叠的思想，并进一步提出了Folded-SIＲ，Folded-SAVE等方法，可以对矩阵或数组形式的自变量进行多方向降维；但这些方法在处理高阶张量时，存在计算效率较低的问题。Ding和Cook（2015）把传统的SIＲ方法拓展到一般的m阶
·116·统计研究2021年2月
张量预测变量中，得到了一种高阶张量的充分降维方法。Zhong等（2015）进一步探寻更小的降维子空间，通过假定投影向量的m阶可分解结构，得到了比上述Folded-SIＲ等方法更简练的降维结构。但是Zhong等（2015）在求解投影方向时，仅针对二阶张量（矩阵）形式的数据提出了序列迭代算法，并没有对算法和相关估计量的理论性质进行证明。
因此，本文在Zhong等（2015）的研究基础上，对高阶张量充分降维的迭代算法进行了研究，主要针对在实际应用中最广泛的三阶张量数据，提出两种不同的降维算法：张量结构转换算法和张量结构保持算法。首先在理论上证明了算法的收敛性和一定条件下输出的投影向量的相合性，补充了目前已有的部分文献在理论方面的研究欠缺；继而将其应用于当前模式识别的研究热点———人脸彩图像的识别问题；并将分类结果与K-means聚类方法、t-SNE非线性降维方法、多维主成分分析、多维判别分析和张量切片逆回归共五种方法进行对比，通过二维和三维散点图的形式进行了直观的可视化比较，展现了本文方法的实际应用价值。
文章内容的结构安排如下：第二节介绍了张量充分降维模型；第三节在充分降维框架下，提出了基于三阶张量的序列迭代算法：张量结构转换算法和张量结构保持算法；第四节以Georgia Tech 人脸数据为样本，对算法的分类性能进行了检验，并比较了转换算法和保持算法的差异；最后将结构保持算法与文献中已有的其他降维方法进行对比讨论，发现使用本文提出的张量充分降维方法进行人脸识别效果突出，在实际数据分析中具有良好的实证表现。
二、张量充分降维模型
（一）充分降维方法（SDＲ）
假设Y∈瓗为响应变量，X=（x1，…，x p）T∈瓗p为预测变量，满足E（X）=0，cov（X）= X。对于变量维数p非常大的情况，Li（1991）提出了充分降维模型，Chen和Li（1998）、Cook等（2004）、Zhong等（2012）对此模型进行了不断的丰富和发展。模型形式如下：
Y=f（βT
1X，…，βT
K
X，ε）（1）
其中，f（·）是定义在瓗K+1上的未知连接函数，β1，…，βK是p维向量，ε是与X独立的随机误差项。当模型（1）成立时，p维变量X被投影到由β1，…，βK张成的K维子空间S=span｛β1，…，βK｝中。该
子空间包含了Y的全部信息，即
Y⊥X|P
S
X（2）其中，P（·）是基于标准内积的投影算子，S是预测变量空间的子空间，“⊥”代表“相互独立”。
SDＲ模型假设Y和X在给定P
S
X的条件下是相互独立的，这意味着高维预测变量可以投影到低维子空间而不会丢失信息。
未知系数向量β1，…，βK被称为SDＲ方向，由这些方向张成的空间称为SDＲ子空间。SDＲ方法意味着X中包含的所有关于Y的信息都含在了K个投影βT1X，…，βT K X中。
目前已有很多降维方法可以用来估计SDＲ方向，例如切片逆回归（Sliced InverseＲegression）、切片平均方差估计（Sliced Average Variance Estimation）、主海森方向（Principal Hessian Direction）等。其中Li（1991）提出的切片逆回归（SIＲ）凭借其良好的性质成为目前广泛使用的降维方法。
Li（1991）证明了当自变量向量满足线性设计条件，即对任意的b∈瓗p，存在常数c
0，c
1
，…，c
K
，
使得E（b T X|βT1X，…，βT K X）=c0+c1βT1X+…+c KβT K X时，中心化的逆回归曲线E（X|Y）－E（X）包含在由ΣXβk（k=1，…，K）张成的线性子空间中。因此，对SDＲ方向β1，…，βK的估计可以转化为以下的特征值分解问题：
第38卷第2期马少沛等：大数据下张量充分降维方法及其应用研究·117·ΣX|Y u i=λiΣX u i
u T
i
ΣX u i=1，i=1，…，p
λ1≥λ2≥…≥λp
{
其中，ΣX|Y=cov［E（X|Y）］。与Σ－1XΣX|Y的前K个最大特征值对应的特征向量u k（k=1，…，K）即为待求的充分降维方向，将其称为SIＲ方向。
实质上，切片逆回归与以下寻最优投影方向的方法具有等价性：即SIＲ方向实际上就是使得因变量与投影自变量之间的线性相关程度达到最大的投影方向。
从变量转换和投影追踪的角度出发，寻SDＲ方向β1，…，βK也可以通过对P2（η）进行最大化得到：
P2（η）=max
T
corr2（T（Y），ηT X）（3）其中，T（Y）代表对因变量Y进行的某种变换，投影方向η∈瓗p，corr（·，·）计算两个变量间的相关系数。P2（η）旨在寻某种变换形式T，使得变换后的因变量T（Y）与投影变量ηT X之间的相关系数平方达到最大值。Chen和Li（1998）证明了最优变换为T（Y）=E（ηT X|Y），此时P2（η）有以下解析形式：
P2（η）=ηT cov［E（X|Y）］η
ηT cov（X）η
ηTΣX|Yη
ηTΣXη
（4）
因此，求得的第一投影方向η1∈瓗p，使得P2（η1）达到最大；之后，在与η1正交的方向上，寻η2∈瓗p，使得P2（η2）达到次大；以此类推，重复以上过程，直至到一组相互正交的基向量η1，…，ηK，从而张成K维的SDＲ子空间。
由上述分析可知，ηk（k=1，…，K）刚好是对应Σ－1XΣX|Y的第k个最大特征值的特征向量。因此，切片逆回归方法和上述相关系数最大化方法得到的充分降维方向是完全一致的，且有特征根λk=P2（ηk），k=1，…，K。关于ΣX|Y的计算，Li（1991）用切片的思想构造了ΣX|Y的加权样本协方差矩阵估计。
（二）张量充分降维（TSDＲ）
当观测到的自变量是二阶或更高阶的张量时，本文需要考虑将传统的充分降维方法推广到张量型的预测变量。TSDＲ模型假设响应变量依赖于一些张量的低维投影，但连接函数的形式未知。Zhong等（2015）通过假设投影方向满足某种张量结构，大大减少了参数个数，能够重建一个包含预测变量重要回归信息的低维子空间。
假设X∈瓗p1ˑ…ˑp m是一个m阶张量，vec（X）表示向量化的X。设γj=β（1）j … β（m）j是向量β（1）j，…，β（m）j的Kronecker乘积，其中β（i）j∈瓗p i。给定张量预测变量X∈瓗p1ˑ…ˑp m和响应变量Y∈瓗，TSDＲ模型的形式如下：
Y=f（γT
1vec（X），…，γT
K
vec（X），ε）（5）
其中，f（·）为未知函数，γj，j=1，…，K为回归指标，ε为与vec（X）独立的随机误差项。
因此，TSDＲ模型假设在给定（γT1vec（X），…，γT K vec（X））的情况下，Y与vec（X）是相互独立的。为了保持张量结构和模型（5）中指标的可解释性，现引入可分解张量的概念。称向量γ∈瓗p1 … 瓗p m是m阶－可分解的，如果它能写成ν1 … νm的形式，其中νi∈瓗p i，i=1，…，m。记Ｒ是瓗p1 … 瓗p m中所有可分解向量的集合。
回顾Li（1991）定理3.1可知，经典的切片逆回归方法要求自变量向量X的分布满足一定的线
·118·统计研究2021年2月
性设计条件。因此，当把张量X向量化为vec（X）后，为了估计投影方向γ1，…，γK，同样需要考虑vec（X）的分布条件。记cov（vec（X））=Σ
vec（X）
Σv，Γ=（γ1，…，γK）。受到Li等（2010a）的启发，本文发现只要vec（X）的分布仍然满足线性条件，即E（vec（X）|ΓT vec（X））是ΓT vec（X）的线性组合，则中心化的逆回归曲线E（vec（X）|Y）－E（vec（X））将包含在由Σvγk（k=1，…，K）张成的线性子空间中。根据Li（2018），相对于协方差矩阵Σv，到中心子空间S Y|vec（X）=span｛γ1，…，γK｝的
投影矩阵定义为P
Γ（Σ
v
）=Γ（ΓTΣ
v
Γ）－1ΓTΣv。则有以下引理成立：
引理：假设ΓTΣvΓ是正定的，当vec（X）的分布满足线性假设条件时，有
E［vec（X）－E（vec（X））|ΓT vec（X）］=P T
Γ（Σ
v
）［vec（X）－E（vec（X））］
下面将从与Li（1991）定理3.1不同的角度出发，对vec（X）在满足线性设计条件下的性质进行证明。
证明：不失一般性，假设E（vec（X））=0。则有
E（vec（X）|Y）=E［E（vec（X）|ΓT vec（X），Y）|Y］
=E［E（vec（X）|ΓT vec（X））|Y］
=E［P T
Γ（Σ
v
）vec（X）|Y］张雅玫
=P T
Γ（Σ
v
）E［vec（X）|Y］
其中，第二个等式成立是因为Y⊥vec（X）|ΓT vec（X）；第三个等式根据线性假设条件和上述引理
即可得出。因为P T
Γ（Σ
v
）=Σ
v
P
Γ
（Σ
v
）Σ－1
v
，所以第四个等式右边可以写成Σ
v
P
Γ
（Σ
v
）
Σ－1
v
E［vec（X）|
Y］的形式，因此有：
E［vec（X）|Y］=Σ
v P
Γ
（Σ
v
）Σ－1
v
E［vec（X）|Y］
所以可得，Σ－1v E［vec（X）|Y］∈span｛P
Γ（Σ
v
）｝=S
Y|vec（X）
。当E（vec（X））≠0时，只需将中心化后
的逆回归曲线代入上述过程即可得证。证毕。
因此，当vec（X）的分布满足线性设计条件时，从相关系数最大化的角度可以得到式（3）的推广形式：
P2（η）=max
T
corr2（T（Y），ηT vec（X））
其中，对应的投影向量η∈Ｒ。与式（4）类似，根据Chen和Li（1998）的证明过程，P2（η）有以下等价的表示形式：
P2（η）=ηT cov［E（vec（X）|Y）］η
ηT cov（vec（X））η
（6）
通过对P2（η）最大化可以得到第一投影方向（充分降维方向）γ1为：
γ1=arg max
η∈Ｒ
P2（η）
为了求解出全部的充分降维方向，在得到γj（j=1，…，K－1）以后，需要在与γj正交的方向上继续寻γj+1使得P2（γj+1）达到次大。式（6）与式（4）最大的不同之处在于，由张量展开得到的vec（X）通
常有较高的维数，不能通过简单的谱分解技术直接计算。
针对二阶张量（m=2）的情况，Zhong等（2015）提出了一种迭代算法来寻充分降维方向γj。假设X∈瓗p2ˑp1，vec（X）是X的向量化，对于任意的η=β（1）1 β（2）1∈Ｒ，有ηT vec（X）=β（2）1T Xβ（1）1，其中β（i）1∈瓗p i，i=1，2。因此，式（6）的最大化等价于最大化关于β（1）1，β（2）1的下式：
P2（β（1）
1，β（2）
1
）=
cov［E（β（2）
1
T Xβ（1）
1
|Y）］
cov（β（2）
1
T Xβ（1）
1
）
（7）
为了求解式（7），采用下面的迭代算法：
发布评论