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

路子祥，黃嘉爽，屠黎陽，徐西嘉，張道強+

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回歸模型?，F有的研究表明如果樣本較少，同時具有大量不相關的特征時，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算法自帶的特征選擇功能，使得計算出的權重向量具有稀疏性，從而模型更具有泛化能力。在求解中采用交替迭代算法使得模型可求解，并可求得最優權值向量。在真實的三維腦影像數據集上進行實驗，實驗結果也表明了本文方法的有效性。