路子祥，黄嘉爽，屠黎阳，徐西嘉，张道强+

1.南京航空航天大学 计算机科学与技术学院，南京 211106

2.南京医科大学附属南京脑科医院 精神科，南京 210029

1 引言

线性回归是统计学和机器学习领域内重要的分析工具，主要目的是在因变量和自变量之间建立线性关系。近些年来，线性回归被广泛应用于机器学习和模式识别等领域。线性回归的公式可以表示为：

其中，J(·)表示损失函数；Ω(w)表示约束项；w是一维的回归系数向量；·,·表示标准的欧拉内积。常用的线性回归算法LASSO（least absolute shrinkage and selection operator）[1]使用最小二乘作为损失项，一范数作为约束项，通常可以得到绝对值较小的回归系数，其中某些回归系数等于0，因此LASSO同时具有特征选择和岭回归的优点。

在做线性回归时，数据通常需要先进行向量化。然而，现实世界中的数据大多是以多维形式存在的，也就是张量的形式。例如，一幅灰度图片本质上是一个矩阵，或者说是一个二阶的张量，其行和列之间的像素值存在着某种程度的结构关联。然而，在模式识别的早期工作中，以视觉跟踪算法为例[2-6]，大多是将输入图像变换为一个向量。在这个向量化的过程中不仅会破坏数据原始的结构和内在的相关性，而且忽略了数据存在的高阶依赖性。向量化三维脑图像数据通常会致使局部信息难以得到完整的反映，特别是当病变区域集中在海马等体积较小的区域时，使用向量化后的数据常常无法准确定位这些区域。除此之外，数据向量化会导致数据维数过高，计算复杂和存储困难。如一个61×73×61的三维脑影像数据形成的向量维数是271 633维，当样本过小的时候会导致维数灾难和过拟合问题。因此，数据以原始的张量形式解决上述问题更具有意义。

近年来，已经有众多研究人员提出了基于张量模式的学习算法，并被广泛利用在机器学习的领域[7-9]。例如，Sanguansat等人[10]利用基于张量的2D-LDA算法从训练图像中提取判别性特征进行分析。Lu等人[11]提出了一种基于张量对象的多线性主成成分分析算法，在张量的各个模上做特征降维。然而，该方法并不能直接对多线性数据做回归分析。Su等人[12]提出了基于张量的多线性多变量回归算法解决张量数据的回归问题。Guo等人[13]和Zhou等人[14]通过加入张量分解的秩约束，得到了许多局部极小值。然而上述方法均未嵌入特征选择算法，不具备特征选择的功能。

本文旨在解决张量模式的多线性回归问题，提出了一种新的正则化多线性回归算法（multilinear LASSO，mLASSO）。该算法是LASSO（least absolute shrinkage and selection operator）算法在张量空间上的一个扩展，与LASSO不同的是，mLASSO直接在原始的张量数据上分析计算，在张量数据的每一个模式上都加入一范数对权重向量进行稀疏约束，从而在做回归分析的同时具备特征选择的功能。通过基于三维脑影像数据回归临床变量值这一实验来验证本文算法的有效性，实验证明该算法比基于向量的算法表现出了更加良好的回归性能。本文的主要贡献包括以下几点：

（1）算法处理的数据是张量模式，这就避免了将张量数据转化为向量形式的预处理，从而有效地保存了数据信息的同时简化了算法流程。

（2）嵌入特征选择算法，避免无用特征对回归分析的干扰，提高了模型的泛化能力。

本文组织结构如下：第2章介绍了多线性代数基础、多线性回归和LASSO模型；第3章提出了本文模型，并对该模型应用于三阶张量这一特殊情况进行分析求解；第4章介绍了实验所应用的数据集；第5章详细介绍了实验的各种方法之间的比较，并且对实验结果进行了仔细的分析；最后对本文进行总结，并指出进一步的工作。

2 相关工作

本章主要介绍预备知识，包括多线性代数基础、多线性回归模型以及LASSO模型。

2.1 多线性代数基础

为了描述张量问题，首先给出有关张量的常用记号和基本运算[15]。

一个 N 阶张量定义为 A∈RI1×I2×…×IN,并称 in∈{1,2,…,In}为张量A的第n个指标，每个指标对应A的一个模式。若A的第n个指标in变动而其他指标固定，则所得的In维向量称为A的n模向量。将张量A的所有n模向量“展开”所得到的矩阵称为A的n模展开矩阵，记为 A(n)。张量 A ∈ RI1×I2×…×IN与矩阵U∈的n模乘记为 A×nU，定义为(A×nU)(i1,i2,…,张量 A,B ∈的内积定义为iN)·B(i1,i2,…,iN)。张量 A 的 Frobenius范数为 ||A||F= A,A ，并且||A-B||F=||A(n)-B(n)||F=||vec(A)-vec(B)||2，其中vec(A)表示张量A的向量化。张量A的“n模切片”是一个N-1阶张量，通过固定 A的模n值为in:A(:,…,:,in,:,…,:)。如果张量A能表示成N个向量外积 ，即 A=u(1)°u(2)°…°u(N)，其 中 u(n)∈ RIn，n=1,2,…,N，则称A为秩-1张量。

2.2 多线性回归模型

在统计学中，多线性回归是指在观测数据集为张量的情况下，通过回归模型学习一组系数，从而拟合出预测值。N阶张量数据的多线性回归模型如下：

其中，{Xm,m=1,2,…,M}表示M个张量数据的样本，{ym,m=1,2,…,M}表示样本的标签，{u(n)∈RIn,n=1,2,…,N}是学习出的权重向量。

2.3 LASSO模型

LASSO[1]是一种基于一范式的特征选择方法和线性回归算法，与已有的特征选择方法相比较，LASSO不仅可以准确地选择出与类相关的变量，同时还具有特征选择的稳定性。LASSO的基本思想是在回归系数的绝对值之和小于一个常数的约束条件下，使得残差平方和最小化，从而可以得到绝对值较小的回归系数，因此得到解释力较强的回归模型，即LASSO回归模型。现有的研究表明如果样本较少，同时具有大量不相关的特征时，LASSO方法效果显著。LASSO回归估计稀疏表示系数w可描述如下：

其中，X表示样本的数据集；y表示样本的标签；w表示权值向量；λ表示正则化参数。

3 方法

3.1 mLASSO模型

为了解决多线性回归的问题，本节提出了一种基于多线性的mLASSO算法，首先给定如下符号表示：

令{Xm,m=1,2,…,M}表示M个张量数据的样本，其中每个样本维度为2,…,M}表示样本的标签。{λn,n=1,2,…,N}表示N阶张量在N个方向上的正则化参数。

mLASSO算法所要解决的问题是，寻找一组向量{u(n)∈RIn,n=1,2,…,N}可以在张量的各个方向上稀疏约束，并在各个方向上做张量n模乘法运算，使得计算出的结果和样本标签尽可能得相近，即平方差最小。模型公式化描述如下：

3.2 模型求解

目前没有一个解决方案能够同时优化N个向量值。本文算法首先使用张量模乘运算将数据从张量空间变换到向量空间中，然后张量空间的mLASSO问题转化为向量空间内的LASSO模型。接着，将这个问题转换为多个LASSO子问题并采用交替迭代的方法求解各个方向的最优加权向量。最后，使用各个方向的最优加权向量和张量数据做模乘运算得到预测变量值。对于一个三阶张量，算法流程如下：

算法1基于三阶张量的正则化多线性回归算法

注意：式（5）、（6）、（7）与式（3）的优化方法相同，因此任何用于求解LASSO问题的算法都可以用在此处。

3.3 收敛证明

本节提出了一种关于上面迭代计算的收敛证明方法。有以下定理：

定理1解决优化问题（5）、（6）、（7）将单调地减小目标函数值（4），因此mLASSO算法收敛。

证明对于三阶张量数据，定义函数 f如下：

令 v0、w0为初始化值。固定 v0、w0，可以通过式（5）优化求解 u0。同理，固定 u0、w0，可以通过式（6）优化求解 v1；固定 u0、v1，可以通过式（7）优化求解 w1。注意，此处优化求解的LASSO问题是一个凸函数，因此求解LASSO问题将得到一个全局最优解[16-17]。特别的，求解问题（5）、（6）、（7）也将得到全局最优解。因此，有如下不等式：

最终，可以得到：

因此，函数 f是收敛的。同理，可以将三阶张量的证明方法推广到更高阶。 □

4 数据集

4.1 数据说明

本文所分析的数据集来自南京医科大学附属南京脑科医院精神科，包括20个精神分裂症住院患者，表1列出了数据集病人的具体统计信息。

Table 1 Demographic information of subjects表1 被试者信息统计表

数据样本纳入标准包括：

（1）符合难治性精神分裂症的诊断标准；

（2）年龄在20岁至45岁；

（3）均征得患者及其监护人同意，由其法定监护人签署知情同意书。

数据样本排除标准包括：

（1）患有脑器质性疾病、感染性疾病或其他慢性躯体疾病、精神活性物质滥用史以及明显的药物不良反应史；

（2）具有磁共振检查禁忌症。

目前，精神分裂症症状评估的经典量表包括简明精神病评定量表（brief psychiatric rating scale，BPRS）和阳性与阴性症状量表（positive and negative symptom scale，PANSS）。PANSS作为测量精神分裂症的一个标准化测量工具，已经获得了世界精神病学领域的认可和广泛应用。本文数据集中符合标准的难治性精神分裂症患者由两名专业的精神科医师使用PANSS确诊。

4.2 数据采集

所有被试者接受功能磁共振成像扫描，使用Siemens 3.0T成像系统，在标准的头像圈内完成扫描。图像扫描参数如下：重复时间（TR）=3 000 ms，回波时间（TE）=30 ms，层厚为3 mm，翻转角度为90°，视野（FOV）=220 mm×220 mm。获得无缺口的跨轴切片。为了涵盖整个大脑体积，静息态fMRI的扫描时间为6 min。

4.3 数据预处理

对于每一个被试者，图像处理步骤如下：剔除前10个时间点；时间校正；估计扫描期间的头部参数，进行头部校正；将扫描所得结构图像标准化到SPM5的蒙特利尔神经学研究所模板[18]，并将体素重采样为3 mm×3 mm×3 mm。经过头动检测后，所有在各个方向的平动均小于2 mm，转动的角度均小于1°的研究对象入选进行进一步分析。值得注意的是，在这个过程中没有实现时间过滤。这保证了不同的频段可以被用于接下来的分析。经过图像预处理后[19]，由REST软件可以计算出低频振幅[20]。经过预处理操作，可以将原始的四维数据变化成三维数据，如图1所示。

5 实验结果与分析

5.1 实验设置

实验在真实的三阶fMRI张量数据的4个不同频段上进行，三阶分别对应人脑的横轴面、矢状面、冠状面成像数据，在坐标中分别以XYZ坐标轴表示。其中频段包括Slow-5（0.010～0.027 Hz）、Slow-4（0.027～0.073 Hz）、Slow-3（0.073～0.198 Hz）和Slow-2（0.198～0.250 Hz）。每个数据维度均为 61×73×61（271 633个体素值）。为了评估所有比较方法的性能，计算预测PANSS评分和真实PANSS评分之间的相关系数（correlation coefficient，CC）和均方根误差（root mean squared error，RMSE）。相关系数越大说明预测值和真实值这两组数据之间相关性越强。均方根误差越小说明预测模型越能更好地描述实验数据。

本文比较的所有方法均采用留一法进行交叉验证。具体而言，在每次实验中，选取一个样本作为测试样本，其余样本作为训练样本。

Fig.1 Preprocessing changes original four-dimensional data into three-dimensional data图1 预处理将原始的四维原始数据变化成三维数据

5.2 比较方法

本文选取如下的方法进行比较实验：

LASSO模型，将三阶张量数据拉成一维向量，使用LASSO模型做回归预测。

LASSO+SVR模型，将三阶张量拉成一维向量，使用LASSO模型进行特征选择后用SVR做回归预测。

5.3 实验结果

本文通过对阳性和阴性症状量表（PANSS）总分的估计来评估回归性能。表2列出了所有方法在Slow-5、Slow-4、Slow-3、Slow-2共4个频段数据上得到的相关系数和均方误差，图2画出了在Slow-5、Slow-4频段数据上真实值和预测值的散点图。从表2中可以看到，提出的基于张量模式的回归算法性能优于基于向量模式的回归算法性能，即mLASSO算法在三阶张量数据上表现出了良好的回归性能。比如，在Slow-5频段上，本文所提出的模型获得了0.460 6的相关系数值，9.627 9的均方根误差值。相比较本文提出的mLASSO方法，直接使用LASSO模型做回归分析，得到的相关系数只有0.029 5，均方根误差值为16.395 3，这说明数据维度过高导致了LASSO模型失效。在第二个比较方法上，本文用LASSO进行特征选择，对数据首先做一个降维的预处理，预处理后，用支持向量回归机做回归分析。经过多次实验，数据从原有的上万维度降到4 000维度以下时取得最好的回归效果，相关系数值为0.183 3，均方根误差值为10.507 4，相比直接利用LASSO进行回归的结果更好，但效果还是比不过本文提出的方法。因此可以看出mLASSO算法在多线性分析中的效果优于LASSO算法。

Table 2 Regression ability of different frequencies in experiments表2 回归实验中不同频段的回归能力记录表

Fig.2 Scatter plots and respective correlation coefficients obtained by 3 methods图2 3种方法的预测值和真实值之间的散点图和相关系数值

图3显示了在训练集上做回归分析时，目标值的变化情况。从图3中可以看出，目标值在一定的迭代次数后，呈现单调递减的趋势，证明本文方法具有收敛性。

Fig.3 Convergence curve on training set图3 训练集上的收敛曲线图

5.4 最具判别性的脑区

除了预测脑疾病的得分，本文使用算法自带的特征选择功能进一步挖掘出与脑疾病密切相关的脑区。在每次进行交叉验证时，所用的训练集不同导致训练出的权重向量不同，因而所选出的特征体素不同。本文统计每一个特征体素出现在交叉验证中的次数，其中出现次数越多，说明特征体素的判别性越强。根据划分脑区的AAL模板（anatomical automatic labeling）[21]，本文找到每个判别性体素所对应的脑区，作为判别性脑区。表3给出了在Slow-5频段的实验中使用mLASSO算法挑选出的出现频率最高的10个脑区。其中，表3第一列是脑区对应的AAL模板编号，第二列是脑区对应的AAL模板名称，第三列是脑区对应的解剖学名称。图4在一个大脑模板空间中画出表3所示的10个最具有判别性的脑区。如表3所示，其中楔前叶、楔叶、舌回属于默认脑区，该区域被认为是负责人类内在精神活动的核心区域，与精神分裂症存在密切联系[22]。此外，现有的研究证明本文发现的补充运动区与精神分裂症患者所出现的幻觉症状相关[23]。因此，本文的结果和前人的研究结果是一致的，这也表明本文方法能够发现与疾病相关的大脑区域，从而可以辅助医生进一步研究相关的脑疾病。

Table 3 Top 10 significant biomarkers of Slow-5表3 Slow-5频段实验中出现频率最高的10个脑区

Fig.4 Top 10 significant biomarkers of Slow-5图4 Slow-5频段实验中出现频率最高的10个脑区

6 总结

本文提出了一种新的张量数据分析模型mLASSO，用于预测精神病患者的PANSS值。mLASSO使用张量的模乘运算将张量空间变换为向量空间，在向量空间上使用LASSO算法做回归分析，分别计算张量数据各个方向上的权重向量。由于LASSO算法自带的特征选择功能，使得计算出的权重向量具有稀疏性，从而模型更具有泛化能力。在求解中采用交替迭代算法使得模型可求解，并可求得最优权值向量。在真实的三维脑影像数据集上进行实验，实验结果也表明了本文方法的有效性。