新型魯棒孿生支持向量回歸機

陳素根，石婷

1.安慶師范大學 數理學院，安徽 安慶246133

2.安徽省大別山區域復雜生態系統建模、仿真與控制重點實驗室，安徽 安慶246133

3.安徽省皖江流域種群生態模擬與控制國際聯合研究中心，安徽 安慶246133

回歸問題一直是機器學習領域備受關注的問題之一，支持向量回歸機（support vector regression,SVR）是解決回歸問題的重要工具[1]。SVR是通過求解一個較大的凸二次規劃問題（quadratic programming problem，QPP）來構造一個回歸函數，一方面使更多的訓練樣本點落在f(x)-ε和f(x)+ε之間的ε-不敏感帶內，另一方面通過在目標函數中引入正則項以實現結構風險最小化原則，使f(x)盡可能平滑。但是，SVR 訓練速度較慢，為了降低計算成本，受孿生支持向量機（twin support vector machine，TWSVM）[2]的啟發，Peng提出了孿生支持向量回歸機（twin support vector regression，TSVR）[3]，TSVR 是通過求解兩個較小規模的二次規劃問題（QPP）構造兩個非平行ε-不敏感邊界函數，它的訓練時間大約是SVR的1/4。為了克服TSVR存在的問題，Shao等人[4]提出了ε-孿生支持向量回歸機（epsilon twin support vector regression，ε-TSVR），通過在目標函數中引入正則項，實現結構風險最小化，解決了TSVR 的Hessian 矩陣可能出現的奇異性問題。從此，SVR 和TSVR逐漸成為機器學習的研究熱點之一，涌現出了大量相關研究成果[5-7]。雖然取得了豐富的成果，但是SVR、TSVR和ε-TSVR算法依然存在著對噪聲敏感等問題。為了克服噪聲對模型性能的影響，近年來，諸多學者通過在目標函數中調整損失函數來加強模型對噪聲的魯棒性，如：Ye 等人[8]提出了基于二次非凸ε-不敏感損失函數的魯棒支持向量回歸機，該算法可以減少噪聲對模型性能的影響。之后，研究人員通過構造新型損失函數，提出了一系列魯棒的回歸模型[9-15]，如：Dong等人[14]基于Huber損失函數與Pinball損失函數提出了對非對稱噪聲魯棒的支持向量回歸機；馬夢萍等人[15]通過引入pinball損失函數提出了非對稱v-無核二次曲面支持向量回歸機；同時，文獻[11]已經證明了支持向量機原始空間優化的近似解優于對偶空間優化的近似解。受此啟發，在原始空間中求解優化問題的支持向量回歸機被廣泛研究，如：Peng 提出了原空間孿生支持向量回歸機（primal twin support vector regression，PTSVR）[16]；Chen 等人[17]提出了原始空間中基于修剪的Huber 損失函數的魯棒支持向量回歸機。

綜上分析，本文基于ε-不敏感損失函數和Huber 損失函數構造了新型混合Hε損失函數，提出了一種基于原始空間求解的新型魯棒孿生支持向量回歸機，簡稱為Hε-TSVR。該算法有以下優點：（1）混合Hε損失函數是可微的凸函數，繼承了ε-不敏感損失函數和Huber損失函數的優點，使模型具有更好的魯棒性；（2）基于結構風險最小化，在目標函數中加入正則項，可以有效防止過擬合，控制模型的復雜程度；（3）基于牛頓算法在原始空間中求解，無需轉化為對偶問題來求解。最后分別在有噪聲和無噪聲人工數據集、UCI 數據集上進行實驗，實驗驗證了所提算法的有效性。

1 相關工作

給定訓練數據集X={(xi,yi)|i=1,2,…,n}，其中xi∈Rm為輸入樣本，yi∈R 為輸入樣本的輸出值，令A=[x1,x2,…,xn]T∈Rn×m，Y=(y1,y2,…,yn)T∈Rn。下面簡單介紹SVR與TSVR。

1.1 SVR

線性SVR要尋找一個線性回歸函數：

其中，ω∈Rm為法向量，b∈R 為偏置。線性SVR 的原始問題如下：

其中，c>0 為懲罰參數，ξ,ξ*∈Rn為松弛變量，e=(1,1,…,1)T∈Rn。

引入拉格朗日函數，得到式（2）的對偶問題：

其中，α,β∈Rn是拉格朗日乘子。

求解式（2）和式（3）可得ω和b，進而求得式（1）。非線性SVR的相關內容可參見文獻[1]。

1.2 TSVR

線性TSVR要尋找下邊界函數和上邊界函數：

線性TSVR的原始問題如下：

其中，c1,c2>0 為懲罰參數，ε1,ε2≥0 為不敏感參數，ξ,η∈Rn為松弛變量，e=(1,1,…,1)T∈Rn。

引入拉格朗日函數，式（5）的對偶問題如下：

其中，G=[A e]，f=Y-eε1。

類似地，引入拉格朗日函數，可得式（6）的對偶問題如下：

其中，h=Y+eε2。

求解對偶問題式（7）和式（8）得α和γ，代入：

進而得到最終的回歸函數f(x)=(f1(x)+f2(x))/2。對于非線性TSVR的相關內容可參見文獻[3]。

2 Hε-TSVR

2.1 Hε 損失函數

近年來，許多文獻對損失函數進行了研究，常用的損失函數有ε-不敏感損失函數Lε(u)和Huber損失函數Lh(u)。

ε-不敏感損失函數定義如下：

其中，ε為不敏感參數。

該損失函數引入了ε-不敏感區域，落入ε-帶中的樣本點對應損失函數值為0，這使得算法具有稀疏性。

Huber損失函數定義如下：

其中，δ為從二次函數轉換到線性函數的轉換參數。二次函數部分使得該損失函數更適合于具有高斯分布的誤差，而線性函數部分使得該損失函數更適合于具有拉普拉斯分布的誤差，因此Huber損失函數可以對不同誤差進行調節。

受此啟發，本文基于ε-不敏感損失函數Lε(u)和Huber損失函數Lh(u)構造混合Hε損失函數如下：

由式（11）易知，Lhε(u)由三部分組成：當|u|＜ε時，該損失函數值為0，不懲罰小于ε的誤差，使算法更加稀疏；當ε≤|u|≤δ+ε時，損失函數為二次函數，類似于Huber損失函數中的二次函數部分，能有效抑制服從高斯分布的噪聲；當 ||u>δ+ε時，損失函數為一次函數，類似于Huber 損失函數中的線性函數部分，能有效降低較大噪聲和異常數據的影響。圖1 給出了Lε(u)、Lh(u)和Lhε(u)三種損失函數的示意圖。

圖1 Lε(u)、Lh(u)和Lhε(u)的示意圖Fig.1 Illustration of Lε(u), Lh(u) and Lhε(u)

下面給出混合Hε損失函數Lhε(u)的兩個性質：

性質1當δ趨近于0時，Lhε(u)退化為Lε(u)。

證明當δ趨近于0時，Lhε(u)的表達式轉化為：

顯然有Lhε(u)=Lε(u)，即證。

性質2當ε趨近于0時，Lhε(u)退化為Lh(u)。

證明當ε趨近于0時，Lhε(u)的表達式轉化為：

顯然有Lhε(u)=Lh(u)，即證。

2.2 線性Hε-TSVR

下面介紹線性Hε-TSVR 模型及其求解方法，尋求兩個非平行函數：

其中，Gi表示G的第i行。

在目標函數式（13）中，第一項表示正則項，最小化這一項是為了防止過擬合訓練樣本；第二項表示平方誤差損失函數，它嘗試最小化樣本點與下邊界函數之間的平方距離；第三項表示混合Hε損失函數，它使得位于ε帶內的、上方的和下方的樣本點具有不同的懲罰，使整個目標函數更加魯棒。對于目標函數式（14）有類似的解釋，不同的是它嘗試最小化樣本點與上邊界函數之間的平方距離。

將式（11）代入式（13）中：

其中，W1(u1)=diag(s1(u1))，W2(u1)=diag(s2(u1))，W3(u1)=diag(s3(u1))，e=(1,1,…,1)T∈Rn。

根據式（16）計算一階梯度和Hessian矩陣，得：

根據上述分析，構造牛頓算法的迭代公式如下：

通過式（19）可求出u1=[ω1Tb1]T，再由式（12）進一步求出f1(x)=ω1Tx+b1。

類似地，定義向量si(u2)∈Rn(i=5,6,7,8)，它們的第j個分量定義為(u2)(j=1,2,…,n)：

將式（11）代入式（14），得到相應的目標函數，并計算一階梯度和Hessian矩陣，得：

其中，W2=2δ2W5(u2)+c2W6(u2)+c2W7(u2)，S2=δ2c2s5(u2)-c2s6(u2)ε2-c2s7(u2)ε2，W5(u2)=diag(s5(u2))，W6(u2)=diag(s6(u2))，W7(u2)=diag(s7(u2))，I是n×n的單位矩陣。

因此，構造牛頓算法的迭代公式如下：

通過式（23）可求出u2=[ω2Tb2]T，再由式（12）進一步求出f2(x)=ω2Tx+b2。當求出f1(x)和f2(x)后，進而得到最終的回歸函數f(x)=(f1(x)+f2(x))/2 。綜上所述，給出本文線性Hε-TSVR算法如下：

算法1線性Hε-TSVR

2.3 非線性Hε-TSVR

為了將線性Hε-TSVR擴展到非線性Hε-TSVR，先將數據映射到高維特征空間中，再通過引入核函數K(x,xT)來構造模型。尋求兩個非平行函數：

令ui=[ωiTbi]T，H=[K(A,AT)e]，r(u1)=Y-Hu1，q(u2)=Hu2-Y。類似地，構造非線性Hε-TSVR 原始問題如下：

其中，Hi表示H的第i行。

類似于線性Hε-TSVR的推導過程，對于優化問題（25）和優化問題（26），可以分別構造類似于式（19）和式（23）的牛頓算法的迭代公式。通過迭代算法可分別求出u1=[ω1Tb1]T和u2=[ω2Tb2]T，進而計算求出f1(x)=K(xT,AT)ω1+b1和f2(x)=K(xT,AT)ω2+b2。對于非線性Hε-TSVR 算法，可類似于線性Hε-TSVR 算法步驟構造，唯一不同的是需要先選擇合適的核函數K將原始數據映射到高維特征空間，在此不再贅述。

3 實驗結果與分析

為了驗證本文所提Hε-TSVR 算法的有效性，本文分別在無噪聲、有噪聲人工數據集和UCI數據集上進行實驗，并與SVR、TSVR和ε-TSVR 算法進行比較。所有實驗都是在硬件平臺為macOS Mojave，2.4 GHz 的處理器，8 GB 內存的MATLAB R2017b 上進行實驗，其中SVR 算法是直接調用Libsvm 工具箱（https://www.csie.ntu.edu.tw/～cjlin/libsvm/）中的函數來求解，TSVR算法和ε-TSVR 算法采用超松弛技術（successive over relaxation，SOR）求解最優化問題。為了對算法性能進行評價，表1給出了常用的回歸評價指標。

對于非線性算法，都采用高斯核函數進行實驗，即K(u,v)=exp(-||u-v||2/2σ2)，其中σ是核參數。對于參數選取問題，本文采取網格尋優與十折交叉驗證的方式在{2i|i=-5,-4,…,5}范圍內選取最優參數。簡單起見，實驗中設置c1=c2,ε1=ε2,δ1=δ2。

3.1 人工數據集

首先在無噪聲和有噪聲的人工數據集上實驗，驗證算法的有效性。定義Sinc函數如下：

對于無噪聲數據，分別給出252個不帶噪聲的訓練點和500 個不帶噪聲的測試點進行實驗。對于有噪聲數據，將訓練樣本加上零均值高斯噪聲和均勻分布噪聲。特別地，本文選取以下四種類型的人工數據集：

其中，U[a,b]為[a,b]內的均勻隨機變量，N(c,σ2)為服從均值c和方差σ2的高斯隨機變量。分別生成252個帶有噪聲的訓練點和500 個不帶噪聲的測試點進行實驗。每種類型的噪聲數據都做10次獨立實驗，最后取10次實驗結果的平均值與方差。表2和表3分別給出了有噪聲和無噪聲情形下的實驗結果，表中加粗數字表示效果最好的評價指標和最少的訓練時間。

表2 不同算法在無噪聲人工數據集上的結果比較Table 2 Result comparison of different algorithms on artificial datasets without noise

表3 不同算法在有噪聲人工數據集上的結果比較Table 3 Result comparison of different algorithms on artificial datasets with noise

圖2和圖3分別給出了SVR、TSVR、ε-TSVR 和Hε-TSVR 在有噪聲和無噪聲人工數據集上運行一次的擬合效果圖。從表2 和圖2 的實驗結果可以看出，在無噪聲情形下，Hε-TSVR 與其他算法相比都有較好的效果。從表3 和圖3 中的實驗結果可以看出，在有噪聲情形下，Hε-TSVR 在大部分數據集上都取得了較好的效果，各類評價指標均優于SVR、TSVR和ε-TSVR，唯一不足的是Hε-TSVR 的訓練時間相對較長。

圖2 無噪聲人工數據集上的結果Fig.2 Results on artificial datasets without noise

圖3 不同算法在有噪聲人工數據集上的結果Fig.3 Results of different algorithms on artificial datasets with noise

3.2 UCI數據集

為了進一步驗證Hε-TSVR 算法的有效性，本文分別在10 個UCI 數據集Servo、Auto Price、Machine CPU、Auto-Mpg、Wisconsin B.C.、Bodyfat、Boston housing、Concrete CS、Diabetes和Pyrim（http://kdd.ics.uci.edu/）上進行實驗，表4列出了上述數據集的具體信息。

表4 UCI數據集的基本信息Table 4 Basic information of UCI datasets

實驗只考慮有噪聲的情形，先隨機將數據集劃分為70%部分和30%部分，30%部分作為測試集，對70%部分的數據隨機選取10%的數據加入均值為0、標準差為0.2的噪聲數據構成訓練集。由于Diabetes和Pyrim兩個數據集樣本較少，均采用留一法交叉驗證進行實驗，其余UCI數據集上采取十折交叉驗證進行實驗，每個數據集上獨立重復實驗5次，將這5 次結果的平均值作為最終的實驗結果。每個數據集在訓練之前都被標準化到[0,1]區間。

表5 給出了有10%噪聲情形下UCI 數據集上的實驗結果，實驗結果表明了所提Hε-TSVR 算法具有較好的性能。從表5 可看出Hε-TSVR 大部分評價指標要優于SVR、TSVR和ε-TSVR，如在Auto Price數據集上Hε-TSVR 與SVR、TSVR和ε-TSVR 相比有較小的SSE/SST 值和MAE值，較大的SSR/SST值。然而，所提Hε-TSVR 的訓練時間比其他算法相對較長，時間較長的原因可能是所構造的牛頓算法在每一次迭代中都需要計算相應的梯度和Hessian矩陣，并需要計算Hessian矩陣的逆矩陣，當迭代次數較多時，算法運行的時間較長。這個問題需要深入研究，構建快速求解算法，進一步提高模型效率。

表5 有10%噪聲UCI數據集上不同算法的結果比較Table 5 Result comparison of different algorithms on UCI datasets with 10%noise

為了進一步驗證所提算法的有效性，本文利用統計分析中的Friedman test 和Nemenyi test方法來分析各算法之間的差異性[18]。Friedman test定義如下：

服從k-1 自由度的分布，其中vj代表各個算法的平均等級值。

服從k-1 和(k-1)(N-1)自由度的F分布。其中k表示算法的數量，N表示數據集的數量。如果原假設被拒絕，本文將進行Nemenyi test來比較每兩個算法的差異性，如果每兩個算法平均等級的差大于臨界差異CD，認為這兩個算法的性能有顯著差異。臨界差異定義如下：

臨界值qα的計算參見文獻[18]。

表6給出了在有10%噪聲的UCI數據集上SVR、TSVR、ε-TSVR 和Hε-TSVR 算法關于評價指標MAE值的平均等級。從表6可看出，Hε-TSVR 算法的MAE 值的平均等級是最小的。根據式（27）和式（28），可求得?15.72，FF?9.907 6 。由于FF=9.907 6>2.960 0，拒絕原假設。本文進行Nemenyi test，在p=0.1 時的臨界值qα=2.291，根據式（29）計算可得相應的CD值為1.322 7。Hε-TSVR 算法與SVR、TSVR 算法之間的差分別是3.6-1.5=2.1，2.9-1.5=1.4，它們的值都大于CD值1.322 7，結果也表明Hε-TSVR 算法的性能明顯優于SVR 和TSVR算法。

表6 有10%噪聲UCI數據集上不同算法關于MAE值的平均等級Table 6 Average ranks of different algorithms on MAE on UCI datasets with 10%noise

為了驗證算法的魯棒性，進一步分析噪聲強度對所提算法性能的影響，選取上述UCI 數據集中的Servo、Machine CPU、Auto-Mpg、Wisconsin B.C.和Boston housing 這5 個數據集。與前面實驗設置類似，先隨機將數據集劃分為70%部分和30%部分，30%部分作為測試集，再對70%部分的數據分別隨機選取20%和30%的數據加入均值為0、標準差為0.2的噪聲數據構成訓練集，實驗結果分別見表7 和表8。從表7和表8中的結果可以看出本文算法在大多數數據集上依然具有較好的性能，充分說明了本文方法的魯棒性。

表7 有20%噪聲UCI數據集上不同算法的結果比較Table 7 Result comparison of different algorithms on UCI datasets with 20%noise

表8 有30%噪聲UCI數據集上不同算法的結果比較Table 8 Result comparison of different algorithms on UCI datasets with 30%noise

進一步討論核參數σ對Hε-TSVR 性能的影響，以有噪聲Servo 和Boston housing 數據集上實驗為例。方便起見，設置核參數σ1=σ2=σ，核參數的取值范圍為{2-5,2-4,…,24,25}，其他參數c1=c2,δ1=δ2,ε1=ε2都設置為2-1。圖4 給出了核參數對Hε-TSVR 算法評價指標MSE 值的影響，圖4（a）給出了在Servo 數據集上的結果，圖4（b）給出了在Boston housing 數據集上的結果。從圖4（a）中可看出，當核參數σ為2-3時MSE 值最小，在2-5到2-3之間的MSE 值在逐漸減少，在2-3到25之間的MSE 值在逐漸增加。從圖4（b）中可看出，當核參數σ為25時MSE值最小，在2-5到21之間的MSE值在逐漸增加，在21到25之間的MSE值在逐漸減少。

圖4 核參數對Hε-TSVR算法MSE值的影響Fig.4 Influence of kernel parameters on MSE of Hε-TSVR algorithm

4 結束語

本文基于ε-不敏感損失函數與Huber損失函數構造了混合Hε損失函數，并提出了一種新型的魯棒孿生支持向量回歸機（Hε-TSVR）。該算法可以降低噪聲對回歸算法性能的影響。同時，該算法考慮結構風險最小化原則，在目標函數中加入正則項，能有效防止過擬合。在有噪聲和無噪聲的人工數據集、UCI 數據集上進行實驗，實驗結果驗證了本文算法的有效性。然而，本文算法存在訓練時間較長和統計學習理論基礎不清晰等問題，構建快速算法提升效率和探究算法的統計學習理論基礎將是下一步的研究工作。