基于張量的正則化多線性回歸算法及其應(yīng)用*

路子祥，黃嘉爽，屠黎陽(yáng)，徐西嘉，張道強(qiáng)+

1.南京航空航天大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，南京 211106

2.南京醫(yī)科大學(xué)附屬南京腦科醫(yī)院 精神科，南京 210029

1 引言

線性回歸是統(tǒng)計(jì)學(xué)和機(jī)器學(xué)習(xí)領(lǐng)域內(nèi)重要的分析工具，主要目的是在因變量和自變量之間建立線性關(guān)系。近些年來(lái)，線性回歸被廣泛應(yīng)用于機(jī)器學(xué)習(xí)和模式識(shí)別等領(lǐng)域。線性回歸的公式可以表示為：

其中，J(·)表示損失函數(shù)；Ω(w)表示約束項(xiàng)；w是一維的回歸系數(shù)向量；·,·表示標(biāo)準(zhǔn)的歐拉內(nèi)積。常用的線性回歸算法LASSO（least absolute shrinkage and selection operator）[1]使用最小二乘作為損失項(xiàng)，一范數(shù)作為約束項(xiàng)，通?？梢缘玫浇^對(duì)值較小的回歸系數(shù)，其中某些回歸系數(shù)等于0，因此LASSO同時(shí)具有特征選擇和嶺回歸的優(yōu)點(diǎn)。

在做線性回歸時(shí)，數(shù)據(jù)通常需要先進(jìn)行向量化。然而，現(xiàn)實(shí)世界中的數(shù)據(jù)大多是以多維形式存在的，也就是張量的形式。例如，一幅灰度圖片本質(zhì)上是一個(gè)矩陣，或者說(shuō)是一個(gè)二階的張量，其行和列之間的像素值存在著某種程度的結(jié)構(gòu)關(guān)聯(lián)。然而，在模式識(shí)別的早期工作中，以視覺(jué)跟蹤算法為例[2-6]，大多是將輸入圖像變換為一個(gè)向量。在這個(gè)向量化的過(guò)程中不僅會(huì)破壞數(shù)據(jù)原始的結(jié)構(gòu)和內(nèi)在的相關(guān)性，而且忽略了數(shù)據(jù)存在的高階依賴性。向量化三維腦圖像數(shù)據(jù)通常會(huì)致使局部信息難以得到完整的反映，特別是當(dāng)病變區(qū)域集中在海馬等體積較小的區(qū)域時(shí)，使用向量化后的數(shù)據(jù)常常無(wú)法準(zhǔn)確定位這些區(qū)域。除此之外，數(shù)據(jù)向量化會(huì)導(dǎo)致數(shù)據(jù)維數(shù)過(guò)高，計(jì)算復(fù)雜和存儲(chǔ)困難。如一個(gè)61×73×61的三維腦影像數(shù)據(jù)形成的向量維數(shù)是271 633維，當(dāng)樣本過(guò)小的時(shí)候會(huì)導(dǎo)致維數(shù)災(zāi)難和過(guò)擬合問(wèn)題。因此，數(shù)據(jù)以原始的張量形式解決上述問(wèn)題更具有意義。

近年來(lái)，已經(jīng)有眾多研究人員提出了基于張量模式的學(xué)習(xí)算法，并被廣泛利用在機(jī)器學(xué)習(xí)的領(lǐng)域[7-9]。例如，Sanguansat等人[10]利用基于張量的2D-LDA算法從訓(xùn)練圖像中提取判別性特征進(jìn)行分析。Lu等人[11]提出了一種基于張量對(duì)象的多線性主成成分分析算法，在張量的各個(gè)模上做特征降維。然而，該方法并不能直接對(duì)多線性數(shù)據(jù)做回歸分析。Su等人[12]提出了基于張量的多線性多變量回歸算法解決張量數(shù)據(jù)的回歸問(wèn)題。Guo等人[13]和Zhou等人[14]通過(guò)加入張量分解的秩約束，得到了許多局部極小值。然而上述方法均未嵌入特征選擇算法，不具備特征選擇的功能。

本文旨在解決張量模式的多線性回歸問(wèn)題，提出了一種新的正則化多線性回歸算法（multilinear LASSO，mLASSO）。該算法是LASSO（least absolute shrinkage and selection operator）算法在張量空間上的一個(gè)擴(kuò)展，與LASSO不同的是，mLASSO直接在原始的張量數(shù)據(jù)上分析計(jì)算，在張量數(shù)據(jù)的每一個(gè)模式上都加入一范數(shù)對(duì)權(quán)重向量進(jìn)行稀疏約束，從而在做回歸分析的同時(shí)具備特征選擇的功能。通過(guò)基于三維腦影像數(shù)據(jù)回歸臨床變量值這一實(shí)驗(yàn)來(lái)驗(yàn)證本文算法的有效性，實(shí)驗(yàn)證明該算法比基于向量的算法表現(xiàn)出了更加良好的回歸性能。本文的主要貢獻(xiàn)包括以下幾點(diǎn)：

（1）算法處理的數(shù)據(jù)是張量模式，這就避免了將張量數(shù)據(jù)轉(zhuǎn)化為向量形式的預(yù)處理，從而有效地保存了數(shù)據(jù)信息的同時(shí)簡(jiǎn)化了算法流程。

（2）嵌入特征選擇算法，避免無(wú)用特征對(duì)回歸分析的干擾，提高了模型的泛化能力。

本文組織結(jié)構(gòu)如下：第2章介紹了多線性代數(shù)基礎(chǔ)、多線性回歸和LASSO模型；第3章提出了本文模型，并對(duì)該模型應(yīng)用于三階張量這一特殊情況進(jìn)行分析求解；第4章介紹了實(shí)驗(yàn)所應(yīng)用的數(shù)據(jù)集；第5章詳細(xì)介紹了實(shí)驗(yàn)的各種方法之間的比較，并且對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了仔細(xì)的分析；最后對(duì)本文進(jìn)行總結(jié)，并指出進(jìn)一步的工作。

2 相關(guān)工作

本章主要介紹預(yù)備知識(shí)，包括多線性代數(shù)基礎(chǔ)、多線性回歸模型以及LASSO模型。

2.1 多線性代數(shù)基礎(chǔ)

為了描述張量問(wèn)題，首先給出有關(guān)張量的常用記號(hào)和基本運(yùn)算[15]。

一個(gè) N 階張量定義為 A∈RI1×I2×…×IN,并稱 in∈{1,2,…,In}為張量A的第n個(gè)指標(biāo)，每個(gè)指標(biāo)對(duì)應(yīng)A的一個(gè)模式。若A的第n個(gè)指標(biāo)in變動(dòng)而其他指標(biāo)固定，則所得的In維向量稱為A的n模向量。將張量A的所有n模向量“展開(kāi)”所得到的矩陣稱為A的n模展開(kāi)矩陣，記為 A(n)。張量 A ∈ RI1×I2×…×IN與矩陣U∈的n模乘記為 A×nU，定義為(A×nU)(i1,i2,…,張量 A,B ∈的內(nèi)積定義為iN)·B(i1,i2,…,iN)。張量 A 的 Frobenius范數(shù)為 ||A||F= A,A ，并且||A-B||F=||A(n)-B(n)||F=||vec(A)-vec(B)||2，其中vec(A)表示張量A的向量化。張量A的“n模切片”是一個(gè)N-1階張量，通過(guò)固定 A的模n值為in:A(:,…,:,in,:,…,:)。如果張量A能表示成N個(gè)向量外積 ，即 A=u(1)°u(2)°…°u(N)，其 中 u(n)∈ RIn，n=1,2,…,N，則稱A為秩-1張量。

2.2 多線性回歸模型

在統(tǒng)計(jì)學(xué)中，多線性回歸是指在觀測(cè)數(shù)據(jù)集為張量的情況下，通過(guò)回歸模型學(xué)習(xí)一組系數(shù)，從而擬合出預(yù)測(cè)值。N階張量數(shù)據(jù)的多線性回歸模型如下：

其中，{Xm,m=1,2,…,M}表示M個(gè)張量數(shù)據(jù)的樣本，{ym,m=1,2,…,M}表示樣本的標(biāo)簽，{u(n)∈RIn,n=1,2,…,N}是學(xué)習(xí)出的權(quán)重向量。

2.3 LASSO模型

LASSO[1]是一種基于一范式的特征選擇方法和線性回歸算法，與已有的特征選擇方法相比較，LASSO不僅可以準(zhǔn)確地選擇出與類相關(guān)的變量，同時(shí)還具有特征選擇的穩(wěn)定性。LASSO的基本思想是在回歸系數(shù)的絕對(duì)值之和小于一個(gè)常數(shù)的約束條件下，使得殘差平方和最小化，從而可以得到絕對(duì)值較小的回歸系數(shù)，因此得到解釋力較強(qiáng)的回歸模型，即LASSO回歸模型?，F(xiàn)有的研究表明如果樣本較少，同時(shí)具有大量不相關(guān)的特征時(shí)，LASSO方法效果顯著。LASSO回歸估計(jì)稀疏表示系數(shù)w可描述如下：

其中，X表示樣本的數(shù)據(jù)集；y表示樣本的標(biāo)簽；w表示權(quán)值向量；λ表示正則化參數(shù)。

3 方法

3.1 mLASSO模型

為了解決多線性回歸的問(wèn)題，本節(jié)提出了一種基于多線性的mLASSO算法，首先給定如下符號(hào)表示：

令{Xm,m=1,2,…,M}表示M個(gè)張量數(shù)據(jù)的樣本，其中每個(gè)樣本維度為2,…,M}表示樣本的標(biāo)簽。{λn,n=1,2,…,N}表示N階張量在N個(gè)方向上的正則化參數(shù)。

mLASSO算法所要解決的問(wèn)題是，尋找一組向量{u(n)∈RIn,n=1,2,…,N}可以在張量的各個(gè)方向上稀疏約束，并在各個(gè)方向上做張量n模乘法運(yùn)算，使得計(jì)算出的結(jié)果和樣本標(biāo)簽盡可能得相近，即平方差最小。模型公式化描述如下：

3.2 模型求解

目前沒(méi)有一個(gè)解決方案能夠同時(shí)優(yōu)化N個(gè)向量值。本文算法首先使用張量模乘運(yùn)算將數(shù)據(jù)從張量空間變換到向量空間中，然后張量空間的mLASSO問(wèn)題轉(zhuǎn)化為向量空間內(nèi)的LASSO模型。接著，將這個(gè)問(wèn)題轉(zhuǎn)換為多個(gè)LASSO子問(wèn)題并采用交替迭代的方法求解各個(gè)方向的最優(yōu)加權(quán)向量。最后，使用各個(gè)方向的最優(yōu)加權(quán)向量和張量數(shù)據(jù)做模乘運(yùn)算得到預(yù)測(cè)變量值。對(duì)于一個(gè)三階張量，算法流程如下：

算法1基于三階張量的正則化多線性回歸算法

注意：式（5）、（6）、（7）與式（3）的優(yōu)化方法相同，因此任何用于求解LASSO問(wèn)題的算法都可以用在此處。

3.3 收斂證明

本節(jié)提出了一種關(guān)于上面迭代計(jì)算的收斂證明方法。有以下定理：

定理1解決優(yōu)化問(wèn)題（5）、（6）、（7）將單調(diào)地減小目標(biāo)函數(shù)值（4），因此mLASSO算法收斂。

證明對(duì)于三階張量數(shù)據(jù)，定義函數(shù) f如下：

令 v0、w0為初始化值。固定 v0、w0，可以通過(guò)式（5）優(yōu)化求解 u0。同理，固定 u0、w0，可以通過(guò)式（6）優(yōu)化求解 v1；固定 u0、v1，可以通過(guò)式（7）優(yōu)化求解 w1。注意，此處優(yōu)化求解的LASSO問(wèn)題是一個(gè)凸函數(shù)，因此求解LASSO問(wèn)題將得到一個(gè)全局最優(yōu)解[16-17]。特別的，求解問(wèn)題（5）、（6）、（7）也將得到全局最優(yōu)解。因此，有如下不等式：

最終，可以得到：

因此，函數(shù) f是收斂的。同理，可以將三階張量的證明方法推廣到更高階。 □

4 數(shù)據(jù)集

4.1 數(shù)據(jù)說(shuō)明

本文所分析的數(shù)據(jù)集來(lái)自南京醫(yī)科大學(xué)附屬南京腦科醫(yī)院精神科，包括20個(gè)精神分裂癥住院患者，表1列出了數(shù)據(jù)集病人的具體統(tǒng)計(jì)信息。

Table 1 Demographic information of subjects表1 被試者信息統(tǒng)計(jì)表

數(shù)據(jù)樣本納入標(biāo)準(zhǔn)包括：

（1）符合難治性精神分裂癥的診斷標(biāo)準(zhǔn)；

（2）年齡在20歲至45歲；

（3）均征得患者及其監(jiān)護(hù)人同意，由其法定監(jiān)護(hù)人簽署知情同意書(shū)。

數(shù)據(jù)樣本排除標(biāo)準(zhǔn)包括：

（1）患有腦器質(zhì)性疾病、感染性疾病或其他慢性軀體疾病、精神活性物質(zhì)濫用史以及明顯的藥物不良反應(yīng)史；

（2）具有磁共振檢查禁忌癥。

目前，精神分裂癥癥狀評(píng)估的經(jīng)典量表包括簡(jiǎn)明精神病評(píng)定量表（brief psychiatric rating scale，BPRS）和陽(yáng)性與陰性癥狀量表（positive and negative symptom scale，PANSS）。PANSS作為測(cè)量精神分裂癥的一個(gè)標(biāo)準(zhǔn)化測(cè)量工具，已經(jīng)獲得了世界精神病學(xué)領(lǐng)域的認(rèn)可和廣泛應(yīng)用。本文數(shù)據(jù)集中符合標(biāo)準(zhǔn)的難治性精神分裂癥患者由兩名專業(yè)的精神科醫(yī)師使用PANSS確診。

4.2 數(shù)據(jù)采集

所有被試者接受功能磁共振成像掃描，使用Siemens 3.0T成像系統(tǒng)，在標(biāo)準(zhǔn)的頭像圈內(nèi)完成掃描。圖像掃描參數(shù)如下：重復(fù)時(shí)間（TR）=3 000 ms，回波時(shí)間（TE）=30 ms，層厚為3 mm，翻轉(zhuǎn)角度為90°，視野（FOV）=220 mm×220 mm。獲得無(wú)缺口的跨軸切片。為了涵蓋整個(gè)大腦體積，靜息態(tài)fMRI的掃描時(shí)間為6 min。

4.3 數(shù)據(jù)預(yù)處理

對(duì)于每一個(gè)被試者，圖像處理步驟如下：剔除前10個(gè)時(shí)間點(diǎn)；時(shí)間校正；估計(jì)掃描期間的頭部參數(shù)，進(jìn)行頭部校正；將掃描所得結(jié)構(gòu)圖像標(biāo)準(zhǔn)化到SPM5的蒙特利爾神經(jīng)學(xué)研究所模板[18]，并將體素重采樣為3 mm×3 mm×3 mm。經(jīng)過(guò)頭動(dòng)檢測(cè)后，所有在各個(gè)方向的平動(dòng)均小于2 mm，轉(zhuǎn)動(dòng)的角度均小于1°的研究對(duì)象入選進(jìn)行進(jìn)一步分析。值得注意的是，在這個(gè)過(guò)程中沒(méi)有實(shí)現(xiàn)時(shí)間過(guò)濾。這保證了不同的頻段可以被用于接下來(lái)的分析。經(jīng)過(guò)圖像預(yù)處理后[19]，由REST軟件可以計(jì)算出低頻振幅[20]。經(jīng)過(guò)預(yù)處理操作，可以將原始的四維數(shù)據(jù)變化成三維數(shù)據(jù)，如圖1所示。

5 實(shí)驗(yàn)結(jié)果與分析

5.1 實(shí)驗(yàn)設(shè)置

實(shí)驗(yàn)在真實(shí)的三階fMRI張量數(shù)據(jù)的4個(gè)不同頻段上進(jìn)行，三階分別對(duì)應(yīng)人腦的橫軸面、矢狀面、冠狀面成像數(shù)據(jù)，在坐標(biāo)中分別以XYZ坐標(biāo)軸表示。其中頻段包括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）。每個(gè)數(shù)據(jù)維度均為 61×73×61（271 633個(gè)體素值）。為了評(píng)估所有比較方法的性能，計(jì)算預(yù)測(cè)PANSS評(píng)分和真實(shí)PANSS評(píng)分之間的相關(guān)系數(shù)（correlation coefficient，CC）和均方根誤差（root mean squared error，RMSE）。相關(guān)系數(shù)越大說(shuō)明預(yù)測(cè)值和真實(shí)值這兩組數(shù)據(jù)之間相關(guān)性越強(qiáng)。均方根誤差越小說(shuō)明預(yù)測(cè)模型越能更好地描述實(shí)驗(yàn)數(shù)據(jù)。

本文比較的所有方法均采用留一法進(jìn)行交叉驗(yàn)證。具體而言，在每次實(shí)驗(yàn)中，選取一個(gè)樣本作為測(cè)試樣本，其余樣本作為訓(xùn)練樣本。

Fig.1 Preprocessing changes original four-dimensional data into three-dimensional data圖1 預(yù)處理將原始的四維原始數(shù)據(jù)變化成三維數(shù)據(jù)

5.2 比較方法

本文選取如下的方法進(jìn)行比較實(shí)驗(yàn)：

LASSO模型，將三階張量數(shù)據(jù)拉成一維向量，使用LASSO模型做回歸預(yù)測(cè)。

LASSO+SVR模型，將三階張量拉成一維向量，使用LASSO模型進(jìn)行特征選擇后用SVR做回歸預(yù)測(cè)。

5.3 實(shí)驗(yàn)結(jié)果

本文通過(guò)對(duì)陽(yáng)性和陰性癥狀量表（PANSS）總分的估計(jì)來(lái)評(píng)估回歸性能。表2列出了所有方法在Slow-5、Slow-4、Slow-3、Slow-2共4個(gè)頻段數(shù)據(jù)上得到的相關(guān)系數(shù)和均方誤差，圖2畫(huà)出了在Slow-5、Slow-4頻段數(shù)據(jù)上真實(shí)值和預(yù)測(cè)值的散點(diǎn)圖。從表2中可以看到，提出的基于張量模式的回歸算法性能優(yōu)于基于向量模式的回歸算法性能，即mLASSO算法在三階張量數(shù)據(jù)上表現(xiàn)出了良好的回歸性能。比如，在Slow-5頻段上，本文所提出的模型獲得了0.460 6的相關(guān)系數(shù)值，9.627 9的均方根誤差值。相比較本文提出的mLASSO方法，直接使用LASSO模型做回歸分析，得到的相關(guān)系數(shù)只有0.029 5，均方根誤差值為16.395 3，這說(shuō)明數(shù)據(jù)維度過(guò)高導(dǎo)致了LASSO模型失效。在第二個(gè)比較方法上，本文用LASSO進(jìn)行特征選擇，對(duì)數(shù)據(jù)首先做一個(gè)降維的預(yù)處理，預(yù)處理后，用支持向量回歸機(jī)做回歸分析。經(jīng)過(guò)多次實(shí)驗(yàn)，數(shù)據(jù)從原有的上萬(wàn)維度降到4 000維度以下時(shí)取得最好的回歸效果，相關(guān)系數(shù)值為0.183 3，均方根誤差值為10.507 4，相比直接利用LASSO進(jìn)行回歸的結(jié)果更好，但效果還是比不過(guò)本文提出的方法。因此可以看出mLASSO算法在多線性分析中的效果優(yōu)于LASSO算法。

Table 2 Regression ability of different frequencies in experiments表2 回歸實(shí)驗(yàn)中不同頻段的回歸能力記錄表

Fig.2 Scatter plots and respective correlation coefficients obtained by 3 methods圖2 3種方法的預(yù)測(cè)值和真實(shí)值之間的散點(diǎn)圖和相關(guān)系數(shù)值

圖3顯示了在訓(xùn)練集上做回歸分析時(shí)，目標(biāo)值的變化情況。從圖3中可以看出，目標(biāo)值在一定的迭代次數(shù)后，呈現(xiàn)單調(diào)遞減的趨勢(shì)，證明本文方法具有收斂性。

Fig.3 Convergence curve on training set圖3 訓(xùn)練集上的收斂曲線圖

5.4 最具判別性的腦區(qū)

除了預(yù)測(cè)腦疾病的得分，本文使用算法自帶的特征選擇功能進(jìn)一步挖掘出與腦疾病密切相關(guān)的腦區(qū)。在每次進(jìn)行交叉驗(yàn)證時(shí)，所用的訓(xùn)練集不同導(dǎo)致訓(xùn)練出的權(quán)重向量不同，因而所選出的特征體素不同。本文統(tǒng)計(jì)每一個(gè)特征體素出現(xiàn)在交叉驗(yàn)證中的次數(shù)，其中出現(xiàn)次數(shù)越多，說(shuō)明特征體素的判別性越強(qiáng)。根據(jù)劃分腦區(qū)的AAL模板（anatomical automatic labeling）[21]，本文找到每個(gè)判別性體素所對(duì)應(yīng)的腦區(qū)，作為判別性腦區(qū)。表3給出了在Slow-5頻段的實(shí)驗(yàn)中使用mLASSO算法挑選出的出現(xiàn)頻率最高的10個(gè)腦區(qū)。其中，表3第一列是腦區(qū)對(duì)應(yīng)的AAL模板編號(hào)，第二列是腦區(qū)對(duì)應(yīng)的AAL模板名稱，第三列是腦區(qū)對(duì)應(yīng)的解剖學(xué)名稱。圖4在一個(gè)大腦模板空間中畫(huà)出表3所示的10個(gè)最具有判別性的腦區(qū)。如表3所示，其中楔前葉、楔葉、舌回屬于默認(rèn)腦區(qū)，該區(qū)域被認(rèn)為是負(fù)責(zé)人類內(nèi)在精神活動(dòng)的核心區(qū)域，與精神分裂癥存在密切聯(lián)系[22]。此外，現(xiàn)有的研究證明本文發(fā)現(xiàn)的補(bǔ)充運(yùn)動(dòng)區(qū)與精神分裂癥患者所出現(xiàn)的幻覺(jué)癥狀相關(guān)[23]。因此，本文的結(jié)果和前人的研究結(jié)果是一致的，這也表明本文方法能夠發(fā)現(xiàn)與疾病相關(guān)的大腦區(qū)域，從而可以輔助醫(yī)生進(jìn)一步研究相關(guān)的腦疾病。

Table 3 Top 10 significant biomarkers of Slow-5表3 Slow-5頻段實(shí)驗(yàn)中出現(xiàn)頻率最高的10個(gè)腦區(qū)

Fig.4 Top 10 significant biomarkers of Slow-5圖4 Slow-5頻段實(shí)驗(yàn)中出現(xiàn)頻率最高的10個(gè)腦區(qū)

6 總結(jié)

本文提出了一種新的張量數(shù)據(jù)分析模型mLASSO，用于預(yù)測(cè)精神病患者的PANSS值。mLASSO使用張量的模乘運(yùn)算將張量空間變換為向量空間，在向量空間上使用LASSO算法做回歸分析，分別計(jì)算張量數(shù)據(jù)各個(gè)方向上的權(quán)重向量。由于LASSO算法自帶的特征選擇功能，使得計(jì)算出的權(quán)重向量具有稀疏性，從而模型更具有泛化能力。在求解中采用交替迭代算法使得模型可求解，并可求得最優(yōu)權(quán)值向量。在真實(shí)的三維腦影像數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果也表明了本文方法的有效性。