增量式約簡最小二乘孿生支持向量回歸機

曹 杰，顧斌杰，熊偉麗，潘 豐

江南大學 輕工過程先進控制教育部重點實驗室，江蘇 無錫214122

支持向量機（support vector machine,SVM）自20世紀90 年代由Vapnik[1]提出以來，受到廣泛的關注。SVM的主要思想是通過最大化間隔實現結構風險最小化，使得模型具有良好的學習性能和泛化性能[2-3]。目前，SVM 已經成功運用到特征選擇[4]、文本分類[5]、時間序列預測[6]、目標跟蹤[7]和行人檢測[8]等領域。之后，Collobert 等人[9]把支持向量機的思想用于解決回歸問題，提出了支持向量回歸機（support vector regression，SVR）。

近年來，SVM 又衍生出了新的算法。2007 年，Jayadeva等人[10]提出了孿生支持向量機（twin support vector machine，TSVM）。TSVM 通過求解兩個小規模的二次規劃問題，構建兩個非平行超平面，訓練時間縮短到原SVM的1/4左右，并且同樣具有很好的泛化性能。2010 年，Peng[11]基于TSVM 的思想，提出了孿生支持向量回歸機（twin support vector regression，TSVR），通過尋找兩個非平行的超平面，構建回歸模型。實驗結果表明，TSVR在訓練速度以及泛化性能方面均優于傳統SVR。因此，TSVR迅速成為了支持向量機領域的研究熱點。同年，Singh 等人[12]提出了約簡孿生支持向量回歸機，通過隨機選擇部分支持向量構成行列不等的核矩陣，獲得稀疏解，加快了訓練和預測的速度，但預測精度有所降低。2013 年，Shao 等人[13]提出了ε-孿生支持向量回歸機，不僅改進了損失函數，還通過添加正則化項實現了結構風險的最小化。2014年，盧振興等人[14]提出了最小二乘孿生支持向量回歸機，用等式約束代替不等約束，極大地加快了訓練的速度。2016 年，Rastogi 等人[15]提出了v型孿生支持向量回歸機，通過在目標函數中引入vε項，把不敏感參數ε作為變量，實現了自動控制ε的目的。

然而，以上方法都是關于離線訓練算法的研究。在實際應用中，由于數據的更新會使舊模型失效，無法滿足實際的需要。而增量學習算法能夠充分利用歷史的訓練結果，減少模型更新所需要的時間，能夠很好地解決時變問題。目前為止，在增量式SVM 以及增量式SVR 的研究上，學者們取得了很多研究成果，如Cauwenberghs 等人[16]提出的增量式與減量式SVM，何麗等人[17]提出的自適應SVM 增量算法，Ma等人[18]提出的精確在線支持向量回歸機，以及顧斌杰等人[19]提出的精確增量式在線ν型支持向量回歸機，都是在原有模型的基礎上通過少量迭代，更新模型，使得所有樣本滿足KKT條件。此外，張浩然等人[20]提出了最小二乘支持向量回歸機的增量學習算法，Zhao 等人[21]提出了在線不相關約簡最小二乘支持向量回歸機，避免了求解二次規劃問題，算法的實時性得到很大的提升，并且后者能夠獲得稀疏解，使得預測速度也得到了提升。

目前為止，關于增量式孿生支持向量回歸機的研究并不多見。2016年，郝運河等人[22]提出了增量式最小二乘孿生支持向量回歸機。該方法基于最小二乘孿生支持向量回歸機模型，在增量之前，隨機選取部分支持向量生成預訓練模型，在接下來的增量過程中，只增加核矩陣的行向量而不增加核矩陣的列向量，因此求解過程中逆矩陣的維數不變，加快了增量更新模型的速度，并且使得模型解的規模不變。但是這樣獲得的核矩陣可能會忽略掉原核矩陣中的線性無關的列向量，使得核矩陣無法很好地逼近原核矩陣，導致模型的泛化性能與離線算法相比大大下降。為此，本文提出一種增量式約簡最小二乘孿生支持向量回歸機算法（incremental reduced least squares twin support vector regression,IRLSTSVR）。IRLSTSVR

算法嘗試利用約簡方法，選取使核矩陣中列向量線性無關的樣本作為支持向量，構建行列不等的核矩陣，并且通過分塊矩陣求逆引理，充分利用已知計算結果高效更新逆矩陣，在保證模型解的稀疏性的同時，更好地逼近離線算法的泛化性能。最后通過實驗，驗證算法的可行性和有效性。

1 最小二乘孿生支持向量回歸機

最小二乘孿生支持向量回歸機用等式約束替代不等約束，并加入正則項，把經驗風險與結構風險相結合，尋找下界回歸函數和上界回歸函數，其中，ω1,ω2∈Rm為權重向量，b1,b2∈R 為偏置，即得到兩個非平行的超平面。因此，線性回歸問題可以表示為如下約束最優化問題[14]：

其中，C1,C2＞0 為正則化常數，ε1,ε2＞0 為不敏感因子，ξ,η∈Rn為松弛向量，e為相應維數元素全為1的列向量。

其中，I為相應維數的單位陣。

因此，對于某一測試輸入x，可以通過式（7）預測其對應輸出：

在解決非線性回歸問題時，需要借助核函數K(·,·)，將樣本從原始空間映射到一個更高維度的特征空間，則此時的下界回歸函數為f1(x)=K(xT,AT)ω1+b1，上界回歸函數為f2(x)=K(xT,AT)ω2+b2，其中ω1,ω2∈Rn，b1,b2∈R。

令G=[e K(A,AT)]，則非線性回歸問題目標函數的解也可以表示為式（5）和式（6）的形式。對于某一測試輸入x，其對應輸出如下：

2 增量式約簡最小二乘孿生支持向量回歸機

對于線性回歸問題，模型解的規模只和樣本的特征數有關，即輸入量的維數。如果輸入量維數過高，可以預先進行降維，再通過文獻[22]中的方法進行增量訓練，本文不再贅述。然而，對于非線性回歸問題，隨著樣本個數的增加，支持向量的個數和解的維數會不斷增加，導致一次增量更新所需的時間和模型預測所需的時間快速增長，無法滿足實際應用中實時性的要求。此外，現有的增量式最小二乘孿生支持向量回歸機通過隨機選取支持向量構建預訓練模型的方法，可能會導致忽略了一些特征相異的樣本，保留了一些特征相似的樣本，使得核矩陣丟失了原核矩陣中線性無關的列向量，不能很好地逼近原核矩陣。

因此，本文針對非線性回歸問題，嘗試利用約簡方法，選擇特征相異的樣本，使得構建的核矩陣能夠保留原核矩陣中線性無關的列向量，以便更好地逼近原核矩陣。最后以約簡后的列線性無關的核矩陣為基礎，通過增量方法，構建約簡最小二乘孿生支持向量回歸機模型。

2.1 約簡方法

在約簡過程中，把訓練樣本分為三個集合，分別為約簡集S、約束集P和無關集O。其中，約簡集S為約簡后作為支持向量的樣本集合，該集合中樣本的特征差異性較大，它的樣本個數即為核矩陣的列數；約束集P為約簡集S中的樣本以及預測值誤差較大的樣本的集合，它的樣本個數即為核矩陣的行數；無關集O為預測值誤差較小的樣本集合。在增量過程中，忽略無關集O中的樣本。

假設在t時刻模型的S集合中有l1個m維樣本，P集合中有l2個m維樣本，并且P集合包含S集合，即l1≤l2，S集合的樣本矩陣為ASt，P集合的樣本矩陣為APt，核矩陣為。

在t+1 時刻，新增一個樣本(xn+1,yn+1)，如果把該樣本同時加入S集合和P集合中，則形成的新的核矩陣如下：

在選擇同時加入到S集合和P集合中的樣本時，需要考慮的就是核矩陣Kt+1中的列向量nt+1和矩陣Nt+1中的列向量是否線性相關，即式（10）是否有解：

其中，αt+1∈Rl1為線性方程組式（10）的解。

由于精確的線性相關要求會導致數值的不穩定性，因此通過求解如下最小化問題來近似判斷線性相關關系：

把式（13）代入式（12），可求得αt+1的值。

最后，將αt+1的值代入式（11）可求得目標函數的值。

若式（11）的值約等于0，則列向量nt+1和矩陣Nt+1中的列向量線性相關。在實際算法中，設置約簡常數λ∈(0,1)和δ(αt+1)進行比較。如果δ(αt+1)＞λ，則列向量nt+1和矩陣Nt+1中的列向量線性無關，新增樣本(xn+1,yn+1)同時添加到S集合和P集合，否則需要通過式（14）判斷樣本是否位于t時刻模型上界函數和下界函數之間：

如果不滿足式（14）則說明t時刻模型不滿足當前需求，新增樣本(xn+1,yn+1)需要添加到P集合中，否則新增樣本(xn+1,yn+1)添加到O集合中，無需更新模型。

在明確當前新增樣本的所屬集合之后，為了判斷下一次新增樣本的所屬集合，以下兩種情況需要更新Kt+1和Mt+1：

情況1新增樣本添加到S集合和P集合中，Kt+1更新同式（9），Mt+1的更新公式如下：

引理2設A是n×n的可逆矩陣，b是n×1 的向量，d是標量，且d-bTA-1b≠0，則有：

引理2的證明見文獻[24]，此處省略。

令Z=V-1v，J=q-vTZ，則由引理2可得：

式（16）中的V-1由式（13）通過引理1 求得，而J為一個標量，因此整個更新無需矩陣求逆。

情況2新增樣本只添加到P集合中，更新公式如下：

式（18）由式（13）通過引理1 獲得，同樣，以上更新過程無需矩陣求逆。

上述約簡方法主要是針對核矩陣為半正定矩陣，存在線性相關的列向量的情況，對核矩陣進行約簡。由于添加了線性相關的列向量，導致核矩陣中含有值為0的特征值，對應核矩陣中的冗余信息。因此剔除線性相關的列向量，在減小核矩陣規模的同時，可以獲得一個近似原核矩陣的約簡核矩陣。判斷線性相關的一般方法為求解相應的線性方程組，但是線性方程組的等式要求過于嚴格，不便于求解。約簡方法采用最小化問題替代線性方程組，度量新增列向量與核矩陣中已有列向量之間的相關程度，并且為了便于求解采用二次損失函數。如果線性相關，那么最小化問題會收斂到0 值附近，否則就會大于0值。

在實際操作中，采用非精確的線性相關性要求。當最小化問題目標函數的值小于一個接近于0的約簡常數時，滿足線性相關要求，該新增列向量不添加到約簡核矩陣中，使得原核矩陣中約等于0的特征值被剔除。但是由于剔除的特征值并非嚴格為0，因此有效樣本信息會有一定的損失。并且隨著約簡常數的增大，較大的特征值也被剔除，約簡核矩陣的有效樣本信息損失增大，列向量個數也會越少，相應解的規模也會越小。對于約簡常數的設置，往往需要權衡核矩陣的樣本信息損失和解的稀疏程度。

2.2 非線性增量式約簡最小二乘孿生支持向量回歸機

2.1 節詳細闡述了劃分S集合、P集合和O集合的方法，而非線性增量式約簡最小二乘孿生支持向量回歸機就是基于S集合和P集合構建的凸優化問題：

2.3 增量式約簡最小二乘孿生支持向量回歸機算法步驟

增量式約簡最小二乘孿生支持向量回歸機的步驟歸納如下，其中，步驟1～步驟3為模型初始化，步驟4～步驟7為增量更新模型：

步驟1設置參數C1、C2、ε1、ε2、λ。

步驟2把已有樣本劃分為S集合、P集合和O集合，并且計算Kt和Mt（詳細方法見2.1節）。

步驟3由式（23）和式（24）計算t時刻模型的解u1t和u2t，同時保存模型更新時所需的量。

步驟4在t+1時刻，讀入一個新的樣本(xn+1,yn+1)。

步驟5由式（11）～式（13）可求得δ(αt+1)。

步驟6如果δ(αt+1)＞λ，則先由式（9）、式（15）和式（16）更新Kt+1和Mt+1，再由式（25）～式（30）求得u1(t+1)以及下一次增量更新所需的值，同理可得u2(t+1)及下一次增量更新所需的值；否則，判斷式（14）是否成立，如果成立，則直接舍棄樣本(xn+1,yn+1)；不成立則先由式（17）、式（18）更新Kt+1和Mt+1，再由式（31）～式（34）求得u1(t+1)以及更新下一次增量所需的值，同理可得u2(t+1)及下一次增量更新所需的值。

步驟7如果繼續讀入一個新樣本，則重復步驟4～步驟6；否則，算法結束。

2.4 運算復雜度分析

針對2.1 節和2.2 節中闡述的增量過程，下面詳細分析本文算法新增一個樣本所需的時間復雜度。考慮到乘法的時間復雜度遠高于加法，因此以下時間復雜度的分析只考慮算法所需乘法次數。此外，由于構造核矩陣的時間消耗是所有算法都不可避免的部分，本文算法并沒有對此部分進行改進，故以下分析中不考慮構造核矩陣的時間復雜度。

綜上可知，使用矩陣求逆引理進行增量求逆的時間復雜度為平方階，而直接求解逆矩陣的時間復雜度為立方階，因此使用矩陣求逆引理降低了算法的時間復雜度。此外，本文算法需要存儲三個逆矩陣，逆矩陣的大小與S集合的大小有關。并且增量求逆方法對于近似奇異矩陣的求逆精度會有所下降，需要保證矩陣的正定性。

另一方面，本文算法的時間復雜度和l1、l2有關，而l1、l2大小能夠分別由約簡方法中的參數λ和ε1、ε2進行控制。本文算法能夠在盡可能地降低精度損失的前提下，通過調節這三個參數進一步降低算法的時間復雜度。

3 數值實驗與分析

3.1 實驗設計

為了驗證算法的可行性和有效性，選取離線最小二乘支持向量回歸機（least squares support vector regression，LSSVR）[25]、離線最小二乘孿生支持向量回歸機（least squares twin support vector regression，LSTSVR）[14]、增量式最小二乘孿生支持向量回歸機（incremental least squares twin support vector regression，ILSTSVR）[22]和本文算法IRLSTSVR，在基準測試數據集上進行對比。所有實驗在Intel i5-8300H（@2.30 GHz）處理器，8 GB 內存的計算機，Matlab2014a 軟件平臺上完成。

實驗中使用的基準測試數據集如表1所示，其中最小規模為167，最大規模為9 568。Servo、Boston Housing、Concrete CS、Airfoil Self Noise、Abalone、Triazines 和CCPP 來 自 于UCI 數 據 庫，PM10 和Space_ga來自于StatLib數據庫。

Table 1 9 benchmark datasets used in experiment表1 實驗中使用的9個基準測試數據集

首先把每組數據集的輸入數據歸一化為[0,1]，然后劃分訓練集和測試集，在訓練集上采用5次五折交叉驗證的方式，共25 次實驗的平均值進行參數尋優，最終以訓練集上的最優模型在測試集上的表現來評判該模型的好壞。采用如下性能評價指標：

其中，RMSE為均方根誤差；MAE為平均絕對誤差；ρ和ψ分別表示RMSE和MAE的損失程度；φ為解的稀疏率；為第i個樣本的預測值；yi為第i個樣本的實際輸出值；N為訓練樣本的個數；Ron為增量算法的RMSE；Roff為相應離線算法的RMSE；Mon為增量算法的MAE；Moff為相應離線算法的MAE；NS為S集合中樣本的個數。

此外，還統計了S集合的樣本個數NS和P集合的樣本個數NP，平均一次增量所需的訓練時間和最終模型對測試集進行預測所需的預測時間作為性能評價指標。

3.2 參數設置

為克服離線算法參數尋優的困難，在本文實驗中，離線LSSVR 的參數C=2i，在i∈[-10,10]的范圍內尋找最優值。離線LSTSVR的參數C1=C2=2i，在i∈[-10,10]的范圍內尋找最優值，考慮到不敏感因子ε1、ε2對離線LSTSVR的影響較小[26-28]以及為了與離線LSSVR進行公平對比，取ε1=ε2=0.01。

對于本文算法，參數λ=10i，i∈[-5,-1]，在保證解的稀疏率的前提下，選擇具有較優RMSE和MAE指標的情況。由于ε1、ε2的大小會對P集合的規模產生影響，并且ε1、ε2的設置不能夠太大，實驗中在相同λ下，選擇ε1=ε2=0 的情況和ε1=ε2≠0 且使P集合的規模產生顯著變化的情況，其他參數使用相應離線算法的尋優參數。對于ILSTSVR算法，隨機選取和本文算法S集合大小相同的樣本訓練預訓練模型，其他參數使用相應離線算法的尋優參數。核函數統一選取高斯徑向基函數K(xi,xj)=exp(-||xi-xj||2/2σ2)，核參數σ=2i在i∈[-5,5]的范圍內尋優。

3.3 實驗結果與分析

表2 所示為選取的四種不同算法在基準測試數據集上的實驗結果，其中IRLSTSVR 列舉了兩種情況，一種是ε=ε1=ε2=0 的情況，另一種是ε=ε1=ε2≠0 的情況，“—”表示該處指標無意義，粗體表示最優指標。

從表2 可以看出，對于本文算法，當ε=0 時，ρ最大為4.4%，Ψ最大為4.1%，φ最小為3.8%，最大為57%。當ε≠0 時，ρ最大為5.6%，Ψ最大為13%，φ最小為3.7%，最大為54%。對于ILSTSVR算法，ρ最大為135%，Ψ最大為97%，φ最小為3.8%，最大為57%。與離線LSSVR 和離線LSTSVR 相比，本文算法具有相近的泛化性能，解的稀疏性大為提高，預測時間減少以及平均一次增量更新所需的時間大為減少。與ILSTSVR相比，本文算法的泛化性能更好，更加接近離線算法，但是平均一次增量更新所需的時間略大。原因是本文算法在添加新樣本時，進行了約簡，考慮了核矩陣列向量之間的線性相關性。由于特征相似樣本的添加會導致核矩陣中含有線性相關的列向量，利用凸優化問題來近似判斷列向量之間的線性相關性，就能夠排除特征相似的樣本作為支持向量，獲得一個與原核矩陣近似的約簡核矩陣，有效減小模型解的規模獲得稀疏解，減少預測時間。同時該約簡核矩陣列向量對應的特征值都遠大于0，含有原核矩陣絕大部分的信息，能獲得和基于原核矩陣的離線算法相近的泛化性能；并且充分利用前一次更新時的計算結果，對矩陣求逆進行高效更新；此外和ILSTSVR不同，本文算法并不是只增加核矩陣的行數，當新添加的列向量線性無關時，需要計算添加了列向量與行向量后的逆矩陣，此時計算復雜度要高于只添加行向量的情況，而且多了約簡部分的計算，因此本文算法更新時間會較長。

Table 2 Experimental results of four different algorithms on benchmark datasets表2 四種不同算法在基準數據集上的實驗結果

對于本文算法，與當ε=0 時的情況相比，當ε≠0 時，P集合的個數減少，也就意味著減少了樣本的更新次數，導致平均一次增量更新的時間會減少，但由于ε1、ε2的設置使樣本偏離模型有了上限，因此依然能獲得相近的泛化性能。

為了更好地說明增量過程中RMSE、MAE、總訓練時間三種指標的變化，圖1 給出了在Abalone 數據集上，隨著訓練樣本個數的增加，三種指標的變化過程，其中離線算法只給出最終模型的RMSE和MAE，作為增量算法的對比基準。從圖1（a）和圖1（b）中可以看出，本文算法在RMSE和MAE的下降趨勢方面要明顯優于ILSTSVR，這是因為與ILSTSVR 算法中隨機選取支持向量相比，本文算法通過約簡方法，在樣本增加過程中，能夠更有效地保留原核矩陣中的大部分信息，獲得更加近似于原核矩陣的約簡核矩陣；并且當ε=0 時的模型指標要略勝于當ε≠0 時的模型指標，這是由于當ε≠0 時模型允許少量的誤差存在，使得參與模型更新的樣本減少，模型的性能通常也會有一定的損失。

Fig.1 Performance comparison of different algorithms on Abalone dataset圖1 Abalone數據集上不同算法的性能對比

從圖1（c）中可以看出，ILSTSVR 算法的總訓練時間要小于IRLSTSVR 算法，而且大致呈線性增長。這是由于ILSTSVR 算法每次增加樣本時，僅添加核矩陣的行數，每次更新的時間復雜度是固定的。相當于每次更新僅需做2.4節所述情況2中模型解的更新部分，具體的時間復雜度為4(l1+1)2+3(l1+1)，而l1是由預訓練模型提前確定的，在增量過程中不會增長，因此ILSTSVR 算法的時間復雜度要小于IRLSTSVR 算法,且總訓練時間呈線性增長。而IRLSTSVR算法增加了約簡過程，并且添加核矩陣列向量和行向量的時間復雜度要高于只添加行向量的時間復雜度，即2.4節中所述，情況1時的時間復雜度要高于情況2時的時間復雜度，因此本文算法的總訓練時間要大于ILSTSVR 算法。而當ε≠0 時模型的總訓練時間要大于當ε=0 時，原因是模型允許誤差存在使得參與模型更新的樣本減少，由2.4 節可知，本文算法的時間復雜度由l1、l2決定，參與模型更新的樣本減少，即l2減小，總的時間復雜度也會相應地減小，也就是說參數ε能夠起到控制總訓練時間的作用。

4 結束語

本文結合約簡孿生支持向量回歸機和最小二乘孿生支持向量回歸機的思想，提出了一種增量式約簡最小二乘孿生支持向量回歸機算法。該算法通過約簡支持向量以及矩陣求逆引理，在不進行矩陣求逆的前提下，獲得約簡最小二乘孿生支持向量回歸機模型，實現了解的稀疏化以及高效增量更新。實驗結果表明，本文算法獲得的模型和離線訓練模型具有相近的回歸精度并且能夠獲得稀疏解，與增量式最小二乘孿生支持向量回歸機相比，泛化性能更加接近離線算法，因此本文算法具備可行性和有效性。