肖文鑫 張文文

高斯過程回歸(Gaussian process regression,GPR)廣泛應(yīng)用于機器人、控制系統(tǒng)和航空航天飛行器等領(lǐng)域,如機器人的逆動力學(xué)模型估計[1-3].它是僅利用簡單的線性代數(shù)處理非線性模型的有效工具,提供了一種簡單但有效的方法來表示數(shù)據(jù)的先驗分布,其特點是使用較少的參數(shù)得到準(zhǔn)確的學(xué)習(xí)結(jié)果[1].作為應(yīng)用范圍廣泛的非參數(shù)回歸,它在樣本外預(yù)測測試中脫穎而出,已成為機器學(xué)習(xí)的重要手段[4].然而,隨著現(xiàn)實世界中訓(xùn)練集和測試集的增長,高斯過程回歸模型在面對大數(shù)據(jù)量時暴露了缺點.當(dāng)數(shù)據(jù)量n很大時,標(biāo)準(zhǔn)高斯過程模型的立方復(fù)雜度 O(n3)會導(dǎo)致不可接受的計算量,使其不適用于高維系統(tǒng)或在線學(xué)習(xí)[5].

降低這種復(fù)雜度的解決方案主要分為兩種類型,其中一種方法是稀疏高斯過程回歸和高斯回歸混合模型[6-7].另一種是將數(shù)據(jù)分割并使用多個局部模型代替全局高斯過程模型,這種方法稱為局部高斯過程回歸(Local GPR,LGPR).在稀疏高斯過程中,按照設(shè)定的規(guī)則減少輸入點的數(shù)量,以降低計算復(fù)雜度[8].局部高斯過程回歸用多個局部子模型替換整體高斯過程模型,由于每個局部模型具有較少數(shù)量的數(shù)據(jù)點并且被單獨訓(xùn)練和更新,因此總體計算成本降低[4].對于局部高斯過程回歸,最關(guān)鍵的是如何劃分?jǐn)?shù)據(jù)并將其分配給不同的模型.

根據(jù)最近的文獻,夏嘉欣等[9]提出一種基于噪聲輸入的稀疏高斯算法,并將其應(yīng)用于解決人體姿態(tài)估計問題.喬少杰等[10]利用高斯過程回歸預(yù)測移動對象最可能運動軌跡.王傳云等[11]對圖像局部區(qū)域建立高斯混合模型解決圖像的背景模型構(gòu)建問題.Da 等[12]研究了一種基于單任務(wù)聚合方案的局部高斯過程回歸模型.Liu 等[13]提出了一種基于流形學(xué)習(xí)方法的局部高斯過程相似度維持算法.Binois 等[14]使用基于似然的高斯過程回歸方法,將異方差性問題表征為單目標(biāo)函數(shù),使用完整的閉合形式導(dǎo)數(shù),實現(xiàn)了基于庫的優(yōu)化.任志剛等[15]利用加權(quán)的優(yōu)秀樣本預(yù)估高斯均值改進了傳統(tǒng)高斯分布估計算法.Sarkar 等[16]使用高斯過程回歸構(gòu)建了一種魯棒算法,能夠處理不確定性和數(shù)據(jù)中的噪聲,并驗證了相比經(jīng)典諧波分析方法的優(yōu)勢.文獻[17]中的工作通過測量距離來選擇合適的局部模型并最終計算加權(quán)預(yù)測,但該方法在每次預(yù)測時仍然考慮所有模型,因此不完全遵循局部支持的概念.

為了進一步降低計算量,本文擬提出一種基于邊界平滑過渡的局部高斯過程回歸模型,即訓(xùn)練集中的數(shù)據(jù)不是確定性地插入到單個局部模型中,而是按概率分布分配給若干相鄰的局部模型.這意味著數(shù)據(jù)點屬于某一模型的概率分布在特定區(qū)域內(nèi)是連續(xù)的,而在超出該區(qū)域的概率分布則為零.這種分布將每個模型限制在一個定義清晰的區(qū)域,通過僅考慮實際相關(guān)且可能對所考慮的數(shù)據(jù)點有貢獻的模型來提高計算效率方面的性能,也就是說對于特定的數(shù)據(jù)點只需要考慮少量數(shù)據(jù)點所屬范圍內(nèi)的局部模型.另外,模型更新過程中使用矩陣塊的求逆方法進行矩陣求逆,大大減少了計算量.為了實現(xiàn)更新與預(yù)測階段局部模型的快速檢索,本文還引入并利用了來自計算幾何領(lǐng)域的k-d 樹最近鄰搜索算法[18].

1 局部高斯過程回歸模型

1.1 高斯過程回歸

高斯過程回歸是機器人和控制理論領(lǐng)域中用于監(jiān)督機器學(xué)習(xí)的常用技術(shù),它的優(yōu)點是通過少量的超參數(shù)實現(xiàn)高預(yù)測精度.但是,當(dāng)數(shù)據(jù)點的總數(shù)n過大時,高計算復(fù)雜度 O(n3)將帶來嚴(yán)重問題.對于由n個數(shù)據(jù)點組成的數(shù)據(jù)集 {X,y},每個數(shù)據(jù)點由d維輸入數(shù)據(jù)x和一維輸出數(shù)據(jù)y構(gòu)成.高斯過程回歸假設(shè)輸出數(shù)據(jù)由函數(shù)f(x)和附加噪聲?組成

其中,yi是輸出集y中的第i點,xi是輸入集X中的第i點,?i是均值為零、方差為的正態(tài)分布噪聲,n是數(shù)據(jù)點的數(shù)量.輸出集滿足下面的正態(tài)分布[19]

其中,K(X,X)是通過所有訓(xùn)練點計算得到的協(xié)方差矩陣,I表示單位矩陣.輸出訓(xùn)練集y和對問詢點x*的預(yù)測輸出f(x*)的聯(lián)合分布由下式給出

f(x*)的條件分布服從一個新的正態(tài)分布

其中,f(x*)的均值和方差的表達式為

因此,對問詢點x*的預(yù)測值可以通過下面表達式計算

對于內(nèi)核函數(shù)k(xp,xq)的選擇,本文采用充分靈活的平方指數(shù)協(xié)方差函數(shù),其定義為

使用文獻[20]中提出的共軛梯度算法最大化對數(shù)邊際似然

其中,C=K(X,X)+.由以下對數(shù)邊際似然導(dǎo)數(shù)的表達式可以推導(dǎo)出對數(shù)邊際似然取最大值時的超參數(shù)為

1.2 局部高斯過程回歸

由于高斯過程回歸在訓(xùn)練點數(shù)量大時表現(xiàn)不佳,為了在保持高斯過程回歸的優(yōu)勢的同時降低計算復(fù)雜度,提出了局部高斯過程回歸.局部高斯過程回歸可以實現(xiàn)模型的快速更新,因此更適合使用在需要通過實時產(chǎn)生的數(shù)據(jù)不斷更新模型的領(lǐng)域.局部高斯過程回歸方法背后的基本思想是引入多個回歸模型并減少每個模型的數(shù)據(jù)點數(shù)量,從而降低計算復(fù)雜度.該方法將數(shù)據(jù)點分配到多個局部模型并單獨訓(xùn)練[1].對應(yīng)于平方指數(shù)協(xié)方差函數(shù),點與模型中心之間的距離通過如下距離測量來計算,以將數(shù)據(jù)點分配給最近的局部模型

其中,ci表示的是第i個模型的中心,wi是模型i和給定x之間的距離.因此,為了正確分配新數(shù)據(jù)點,這個方法必須為每個數(shù)據(jù)點計算M次距離,其中M是局部模型的總數(shù).基于此度量,選擇具有最小距離的模型,并且將點添加到該最近的模型.一旦選擇了模型,該方法后續(xù)步驟與第1.1 節(jié)中的常規(guī)高斯過程回歸完全相同.如果沒有足夠接近新數(shù)據(jù)點的模型,則會創(chuàng)建一個僅包含此新點的新模型.亦即如果新點到所有中心的距離大于某個閾值,則創(chuàng)建一個新模型,其中心是新數(shù)據(jù)點.

式(12)中的距離度量也用于加權(quán)預(yù)測.使用式(7)對第i個局部模型的預(yù)測值進行計算,然后通過對這些模型的加權(quán)預(yù)測求和來獲得全局預(yù)測值

由于每個局部模型中的數(shù)據(jù)點的數(shù)量比高斯過程回歸中的低,用于更新模型的計算工作量得以減少,并且與高斯過程回歸相比可以實現(xiàn)整體性能改進.在第2.1 節(jié),本文提出的算法對以上方法進行了改進,采用概率數(shù)據(jù)關(guān)聯(lián)將數(shù)據(jù)點在其最近的 2d個局部模型中進行分配.

1.3 高效k 近鄰查詢樹

為了實現(xiàn)快速搜索算法,本文使用k-d 樹來存儲數(shù)據(jù).利用k-d 樹進行局部模型檢索,可以有效提高模型更新與對未知數(shù)據(jù)點預(yù)測的速度.k-d 樹是二叉搜索樹(Binary search tree,BST)的多維形式,而二叉搜索樹是用于快速可靠查詢的合適方法[18].文獻[21]中提出了這個k-d 樹的優(yōu)化版本.

1.3.1 k-d 樹的創(chuàng)建

對于樹的第零層,首先在第1 個維度分割數(shù)據(jù),并使用此維度中數(shù)據(jù)的中值作為分離的邊界.此步驟將數(shù)據(jù)集分為兩部分:一部分是由小于或等于中值的數(shù)據(jù)組成的樹的左側(cè);另一部分是剩余數(shù)據(jù)構(gòu)成的樹的右側(cè).在每個子樹上重復(fù)此過程,直到只剩下包含一個點的節(jié)點,將其作為樹的子節(jié)點.由于在每個節(jié)點中重復(fù)劃分并將剩余的邊處理為新的子樹,因此這是一個遞歸定義.分割的維度是根據(jù)樹中當(dāng)前節(jié)點的級別來決定的,即在層級i上,選擇維度i+1 進行分割.如果迭代遍歷了所有維度,則在樹的下一層過程將再次從第1 個維度開始分割,直到只剩下具有一個數(shù)據(jù)點的子集表示樹的子節(jié)點.

1.3.2 k-d 樹最近鄰搜索

k-d 樹最近鄰搜索是一種快速且可靠的搜索算法,其基于預(yù)定義的相似性或距離度量來判斷距查詢點最近的k個局部模型.對于如何找到最佳匹配和k個最近鄰項存在多種方法[22-23],這項工作將搜索算法制定為遞歸策略,其中算法的輸入是每次遞歸對應(yīng)的節(jié)點和其子節(jié)點.首先,包含所有數(shù)據(jù)的根節(jié)點作為首次的輸入傳遞到算法.當(dāng)前最近鄰項的候選項被存儲在優(yōu)先隊列中,其中查詢點和近鄰項之間的相似性是判斷的依據(jù).接下來會將每個節(jié)點與優(yōu)先隊列中的條目進行比較,如果該節(jié)點比隊列中的一個條目更接近查詢點,則相應(yīng)地更新優(yōu)先隊列使該節(jié)點成為新的候選項.如果節(jié)點不是終端,即不是葉子節(jié)點,則查詢點所在的子樹再次遞歸地傳遞給算法.如果節(jié)點是葉子節(jié)點,則此子樹上的遞歸結(jié)束.k-d 樹最近鄰搜索的計算復(fù)雜度是O(k×n×logn),其中k是數(shù)據(jù)的維度,n是數(shù)據(jù)點的個數(shù).

2 基于概率關(guān)聯(lián)進行模型選擇的局部高斯過程回歸

本研究提出一種無硬邊界局部模型概率關(guān)聯(lián)方法,這種方法確保了在僅有局部支持的條件下實現(xiàn)模型之間的平滑過渡.模型建立階段,通過定義模型相鄰部分的概率分布,數(shù)據(jù)點被根據(jù)其坐標(biāo)分配到局部模型中.對未知數(shù)據(jù)預(yù)測階段,由距離查詢點最近的若干模型計算加權(quán)平均值,從而求和得到全局預(yù)測值.因此,該算法包括3 個階段:1)初始模型的建立;2)局部模型的更新;3)測試集的加權(quán)預(yù)測.第2.1 節(jié)將討論局部模型的概率型數(shù)據(jù)關(guān)聯(lián),第2.2 節(jié)和第2.3 節(jié)描述了模型的建立和模型的更新,第2.4 節(jié)提出了對未知數(shù)據(jù)的預(yù)測方法.

2.1 邊界約束的局部高斯過程回歸模型

緊致支持的邊界約束局部模型是本在線學(xué)習(xí)模型的基本組成部分,本節(jié)介紹如何構(gòu)建本地模型.

根據(jù)文獻[1,17]局部高斯過程回歸方法通常使用模型內(nèi)數(shù)據(jù)點的平均值作為每個模型的中心.因此,對于每個新輸入的數(shù)據(jù)點,這些中心都需要更新,即局部模型需要持續(xù)地改變.為了減少每次更新造成的計算量,本文提出一種固定的模型網(wǎng)絡(luò),同時為了實現(xiàn)局部模型與數(shù)據(jù)之間的概率關(guān)聯(lián),我們提出了一種激活函數(shù).激活函數(shù)φ(ci|x)表示數(shù)據(jù)點x屬于局部模型ci的概率,滿足概率分布的條件.距數(shù)據(jù)點最近的 2d個局部模型對x的激活函數(shù)可以組成一個離散概率密度函數(shù)f(c),其中,d是輸入數(shù)據(jù)點x的維度,c=c1,···,ci,···,,ci是局部模型i的中心,本文以模型中心代指具體的模型.最后,根據(jù)f(c)中的概率分布產(chǎn)生一個偽隨機數(shù),以決定數(shù)據(jù)點 (x,y)分配到哪個局部模型中.

所有局部模型組成一個覆蓋數(shù)據(jù)集空間的模型網(wǎng)格,在維度為一維時,該網(wǎng)格由每個中心之間的距離w定義.當(dāng)維度k增加時,在所有維度中應(yīng)用這種設(shè)置來延伸網(wǎng)格,從而生成由各局部模型組成k維立方體.

相鄰的局部模型在共享邊界處做出不同的預(yù)測,導(dǎo)致邊界處的預(yù)測是不連續(xù)的[24].實際上,對于幾個模型的交匯處的點,屬于每個模型的概率可能都大于 0,這意味著它不是絕對屬于最近的模型.然而,當(dāng)局部模型被設(shè)計為具有硬邊界時,預(yù)測值僅由最近的模型確定,這可能導(dǎo)致大的誤差.為了實現(xiàn)相鄰局部模型之間的平滑過渡,在創(chuàng)建局部高斯過程回歸模型時,本文中使用緊湊支持的邊界約束來構(gòu)建局部模型.這里采用激活函數(shù)定義模型的邊界約束,以表示某一坐標(biāo)的數(shù)據(jù)點屬于這個模型的概率.在模型中心位置,定義只有這個模型是“激活”的,這意味著在這個位置模型ci的激活函數(shù)φ(ci|x)等于 1,其余模型的激活函數(shù)等于 0,位于這個位置的數(shù)據(jù)點只可能屬于這個模型.通過這樣的設(shè)置,局部模型與其相鄰的各模型共享邊界區(qū)域.由于本算法面向的是工程中在線更新的場景,主要處理的是多維數(shù)據(jù),因此激活函數(shù)應(yīng)當(dāng)盡可能考慮對稱的形式,以避免多維空間中模型構(gòu)造過于復(fù)雜.基于通過簡潔的結(jié)構(gòu)實現(xiàn)更好的性能的理念,將共享區(qū)域的激活函數(shù)定義為范圍從 0 到 1 的斜坡函數(shù).這樣,所有局部模型使用相同形式的激活函數(shù).如圖1 所示,激活函數(shù)可表示為

圖1 一維局部模型激活函數(shù)示意圖Fig.1 Schematic diagram of one-dimensional local model activation function

其中,x是數(shù)據(jù)點在當(dāng)前維度的值,ci是局部模型的中心在當(dāng)前維度的坐標(biāo),w是相鄰局部模型的間距.由圖中可以看出,每個局部模型在當(dāng)前維度的寬度為 2w.對于處于相鄰模型中心ci和ci+1之間的x,滿足

由此該激活函數(shù)實現(xiàn)了數(shù)據(jù)點與鄰近局部模型的概率關(guān)聯(lián).

圖1 為一維情況下局部模型的初始化示意圖.基于這種一維方法,本文開發(fā)了具有更高維度的數(shù)據(jù)的解決方案.對于模型布局和激活功能,基本原理保持不變.但是,每個數(shù)據(jù)點可能被分配進的模型最多為 2d個.將一維網(wǎng)格布局自然延伸至d維,得到多維模型網(wǎng)絡(luò).

圖2 展通過在二維空間內(nèi)構(gòu)建c1,c2,c3,c4為中心的4 個局部模型,展示了訓(xùn)練集數(shù)據(jù)為二維時的局部模型布局和可能激活的模型的情況.示意圖中的c1,c2,c3,c4分別為4 個局部模型的中心.圖中的區(qū)域A由于不與其他局部模型相交,此區(qū)域內(nèi)的點只可能屬于模型c1,區(qū)域B與一維時的情況相同,最多只有兩個模型處于激活狀態(tài),而這些激活中心c1和c2共享相同的第2 維分量.對于區(qū)域A中的點,可能的激活局部模型僅是c1模型;對于區(qū)域B,可能的激活模型是c1和c2模型.C區(qū)域則有4 個可能激活的模型.另外,該圖還顯示了各局部模型在二維立方體的等距網(wǎng)格中的位置.

圖2 二維局部模型分布示意圖Fig.2 Schematic diagram of two-dimensional local model distribution

對于多維數(shù)據(jù),使用輸入數(shù)據(jù)點在每個維度上的最大值和最小值來計算各局部模型的中心.局部模型的數(shù)量和中心可表示為

其中,numi是局部模型在維度i上的個數(shù),和分別是輸入數(shù)據(jù)在維度i上的最大值和最小值,是相鄰模型在維度i上的間距,ci,n表示的是第i維第n個模型中心的位置.維度i上模型中心的最小值設(shè)置為,以保證所有數(shù)據(jù)點在各維度屬于各局部模型的概率和始終為1.因此,圖1 中的最外側(cè)的兩個局部模型均只利用了靠內(nèi)的一半.擴展到圖2 的情況,數(shù)據(jù)點范圍是集中在C區(qū)域中的,其他區(qū)域由于處在整個模型的邊緣部分,未加利用.當(dāng)模型網(wǎng)格擴大時,始終未使用整個模型的邊緣區(qū)域.計算出局部模型網(wǎng)絡(luò)中所有模型中心的坐標(biāo)后,用這些坐標(biāo)值建立k-d 樹,每個模型都用一個葉節(jié)點表示.通過這個數(shù)據(jù)結(jié)構(gòu),可以使用文獻[22]中的快速搜索算法來執(zhí)行有效的最近鄰居查詢,以便后續(xù)更新模型和預(yù)測未知輸出時搜索模型時使用.

初步確定模型的位置后,根據(jù)式(14)的激活函數(shù)將訓(xùn)練集中的所有數(shù)據(jù)點分配到各個局部模型中.在每個維度上模型中心之間的距離都是w,使得局部模型的范圍覆蓋邊長為w的d維立方體.這些立方體的中心即為局部模型的中心.在整個模型的非邊緣區(qū)域,每個位置均由 2d個立方體覆蓋,也就是說可以有 2d個激活模型.為了得到d維情況下的概率分布,首先使用式(14)在每個維度上計算局部模型對應(yīng)的激活函數(shù)φt(ci|x),它們僅取決于該特定維度中的位置,使得它們符合獨立分布的概率.然后通過這些一維值的乘積來計算整體激活函數(shù)φ(ci|x),局部模型ci的全局激活函數(shù)由其所有子概率的乘積得到,即

2.2 局部模型的初始化

本文提出的在線學(xué)習(xí)模型由初始化和更新階段兩個階段建立,本節(jié)解釋了如何使用訓(xùn)練集構(gòu)建初始局部模型.

為了初始化模型,本研究使用了一種對訓(xùn)練集數(shù)據(jù)點進行聚類從而創(chuàng)建局部模型的方法.首先,通過訓(xùn)練集中的數(shù)據(jù)計算各個局部模型中心的坐標(biāo),構(gòu)建第2.1 節(jié)中所述的模型網(wǎng)絡(luò).然后將所有數(shù)據(jù)點分配到相應(yīng)的局部模型中.為了確定最終分配數(shù)據(jù)點的模型ci,必須在具有非零概率的激活模型中確定一個模型.如果直接選擇最有可能的模型,會在模型相鄰處形成硬邊界,導(dǎo)致預(yù)測在這個區(qū)域的不連續(xù)性.為了解決這個問題,對于訓(xùn)練集中的數(shù)據(jù)點 (x,y),通過第1.3.2 節(jié)介紹的的k-d 樹實現(xiàn)k-最近鄰搜索,找到距離數(shù)據(jù)點最近的 2d個局部模型的編號.通過式(18)計算數(shù)據(jù)點對應(yīng)這些局部模型的激活函數(shù)φ(ci|x)的值,可得到至多 2d個非零概率,構(gòu)造由這些概率值組成的離散概率密度函數(shù)f(c),其中的概率值分別對應(yīng)一個局部模型.根據(jù)各模型對應(yīng)的概率從f(c)中抽取隨機數(shù)i,確定激活的局部模型ci,將數(shù)據(jù)點分配到這個模型.每個數(shù)據(jù)點僅分配給一個模型,但與其他局部高斯過程回歸方法不同,這種關(guān)聯(lián)是用概率方法完成的.完成訓(xùn)練集中所有數(shù)據(jù)點的分配后,通過第2.1 節(jié)中的方法計算每個局部模型的預(yù)測矢量α,并創(chuàng)建相應(yīng)的LGPR 模型.

2.3 局部模型的更新

本文提出的緊致支持的邊界約束局部模型最突出的特點是支持在線學(xué)習(xí),這意味著局部模型能夠通過實時數(shù)據(jù)不斷進行更新.模型的在線更新可以應(yīng)對一些現(xiàn)實世界的工程問題,例如機械特性和機器人交互隨時間的變化[25].在更新階段,依據(jù)第2.2節(jié)中分配數(shù)據(jù)點的方法,新輸入點 (x,y)分配到一個局部模型,需要對這個模型的參數(shù)進行更新,供預(yù)測階段調(diào)用.對未知點的預(yù)測是通過式(7)計算得到,其中協(xié)方差函數(shù)K(x*,X)只能在查詢點x*已知的情況下計算,而預(yù)測矢量α=(K(X,X)+)-1y可以在模型更新階段利用已知數(shù)據(jù)更新,供預(yù)測階段直接調(diào)用.

為了更新預(yù)測矢量,首先,將局部模型的數(shù)據(jù)集更新為

其次,預(yù)測矢量α中更新后的協(xié)方差矩陣Kn(Xn,Xn)+項可以表示為

其中,K(X,X)+是更新前的協(xié)方差矩陣,K(X,x)是由X和x的協(xié)方差組成的矢量,K(x,x)+由x的自協(xié)方差和模型的超參數(shù)相加得到.使用式(20)對協(xié)方差矩陣求逆,得到

其中

根據(jù)以上方法得到 (Kn(Xn,Xn)+))-1和yn后,可以通過下式計算更新的預(yù)測矢量,即

由式(21)可知,每次在計算新的協(xié)方差矩陣的逆時,都利用了它在更新前的值以減少計算量.因此,在模型的更新階段需要更新并儲存數(shù)據(jù)集、協(xié)方差矩陣的逆和預(yù)測矢量.

如果在使用新數(shù)據(jù)點更新模型階段檢測到數(shù)據(jù)點不在任何已有局部模型范圍內(nèi),即在儲存局部模型的k-d 樹中沒有搜索到足夠接近數(shù)據(jù)點的模型,則會創(chuàng)建一個僅包含此新點的新模型,其中心是新數(shù)據(jù)點,并補齊模型網(wǎng)格的其余部分.最后使用第1.1 節(jié)中的方法計算這個模型的超參數(shù)θ和預(yù)測矢量α,并創(chuàng)建LGPR 模型,計算新建模型的中心坐標(biāo)和間距參數(shù).

2.4 局部模型對未知數(shù)據(jù)的預(yù)測

基于輸入數(shù)據(jù)預(yù)測輸出數(shù)據(jù)是本局部高斯過程回歸模型的一個重要功能.在多個模型處于激活狀態(tài)的區(qū)域,可以通過激活函數(shù)計算查詢點屬于各局部模型的概率,再通過各模型的預(yù)測結(jié)果計算加權(quán)預(yù)測.首先,通過k-d 樹搜索查找到距查詢點x*最近的 2d個局部模型,使用激活函數(shù)φ(ci|x)算出各模型對應(yīng)的概率,這些概率值組成離散概率密度函數(shù)f(c).基于這些概率值和這 2d個局部模型各自的預(yù)測值,可以計算得到全局加權(quán)預(yù)測值.如第1.1 節(jié)所述,局部模型ci的預(yù)測遵循高斯分布,其平均值為

其中,x*是輸入的查詢點,Xi是分配到模型ci中的輸入數(shù)據(jù)點集,αi是第i個局部模型的預(yù)測矢量.

通過每個局部模型對應(yīng)的概率對其預(yù)測值進行加權(quán)求和,計算對查詢點x*的預(yù)測輸出y*.如前所述,激活模型ci的概率P(ci|x)通過激活函數(shù)φ(ci|x)表示,可以導(dǎo)出y*的期望值的表達式為

因此,最終的全局預(yù)測值表示為

3 實驗結(jié)果及分析

本節(jié)通過實驗,在更新時間、預(yù)測時間和誤差三個方面評估邊界約束的局部高斯過程回歸算法的性能.第3.1 節(jié)通過人工生成的低維數(shù)據(jù),將本文算法與全局高斯過程回歸方法和具有硬邊界的局部高斯回歸過程方法進行比較,驗證了低維數(shù)據(jù)時本算法的優(yōu)勢.第3.2 節(jié)使用來自F16 飛機的真實多維工程數(shù)據(jù)進行實驗并將本算法應(yīng)用于機械手臂上,分析了本算法在多維數(shù)據(jù)和大數(shù)據(jù)量的場景下,實現(xiàn)在線學(xué)習(xí)的可行性.

3.1 低維數(shù)據(jù)實驗

為了驗證該基于邊界約束的局部高斯過程回歸算法的有效性,本文通過模擬帶有噪聲的高斯過程對該算法進行測試,并與全局高斯過程回歸方法和具有硬邊界的局部高斯回歸過程方法進行比較[17,19].對于訓(xùn)練和測試,數(shù)據(jù)集由正弦函數(shù)的組合產(chǎn)生,超參數(shù)通過第1.1 節(jié)中的方法計算.本算法的性能取決于局部模型間距參數(shù)的選擇,本文首先評估它們對預(yù)測結(jié)果的影響.通過下面正弦函數(shù)組合生成二維數(shù)據(jù)集

并為輸出y附加σn=0.1 的噪聲.首先,將各維度的輸入值都限定在 0 到 40 之間,取 6 400 個數(shù)據(jù)點作為總的數(shù)據(jù)集.從中隨機選取 80% 作為訓(xùn)練集,剩余部分用于測試集.使用這組數(shù)據(jù)重復(fù)模擬邊界約束的高斯局部模型 50 次,得到了不同參數(shù)情況下的平均預(yù)測誤差與運行時間,如圖3 所示.

圖3 局部模型參數(shù)對邊界約束模型性能的影響Fig.3 Influence of local model parameters on performance of boundary constraint model

圖中展示了局部模型間距以8 為間隔由4 增加到44 時,用單個數(shù)據(jù)點更新模型和預(yù)測未知輸出所用時間及預(yù)測誤差值的變化情況.小的局部模型寬度意味著總的模型數(shù)量多,即每個局部模型分配的數(shù)據(jù)點較少,對應(yīng)著較大的誤差和較短的模型更新和預(yù)測時間.當(dāng)局部模型寬度增加時,每個局部模型分配的數(shù)據(jù)點也隨之增加,使得預(yù)測結(jié)果更為準(zhǔn)確.但同時由于求協(xié)方差逆矩陣的計算工作量隨著更多數(shù)據(jù)點而上升,消耗的時間也隨模型寬度增加而延長.因此,模型設(shè)計過程中需要對預(yù)測精度和計算時間做出權(quán)衡.由圖中可見,當(dāng)模型寬度增加到20 以后,模型預(yù)測的準(zhǔn)確度變化不大,但更新和預(yù)測時間仍然隨模型寬度增加而大幅增加,因此20 是這個數(shù)據(jù)集范圍內(nèi)較好的局部模型寬度參數(shù).

為了直觀地比較3 種方法的預(yù)測結(jié)果,本文使用以下等式生成一維輸入與輸出數(shù)據(jù)集

圖4 展示了3 種方法的預(yù)測情況,其中本文提出的算法的預(yù)測誤差為0.03292,全局回歸方法的誤差為0.03239,硬邊界的局部回歸方法的誤差為0.50499.在這個數(shù)據(jù)集中,本算法實現(xiàn)了與全局高斯過程回歸同一量級的精度,比硬邊界局部高斯過程回歸方法的精度高一個數(shù)量級.

圖4 對一維測試集的預(yù)測結(jié)果Fig.4 Prediction results for one-dimensional test sets

為了進一步比較3 種算法的性能,本文使用式(24)生成的二維數(shù)據(jù),將3 種方法的模型各重復(fù)50 次,分別得到它們的平均預(yù)測誤差、單個數(shù)據(jù)點的更新和預(yù)測時間,如表1 所示.其中,邊界約束的局部高斯過程模型的模型寬度參數(shù)設(shè)置為20.

表1 3 種方法的性能對比Table 1 Performance comparison of three methods

從表中可以看出,全局高斯過程回歸方法的預(yù)測精度最佳,而本算法的預(yù)測誤差與它處于同一數(shù)量級,硬邊界局部模型方法的預(yù)測誤差約為全局高斯過程回歸方法的80 倍.考慮到全局高斯過程回歸使用全部數(shù)據(jù)構(gòu)建單個模型,本算法的精度稍低于它是可以接受的.本算法通過與硬邊界局部模型方法接近的局部模型數(shù)量(各約100 個),實現(xiàn)了更好的預(yù)測精度.比較更新模型的時間時,全局高斯過程回歸方法的弊端顯現(xiàn)出來,其使用單個數(shù)據(jù)點更新模型的時間為132.753 ms,是局部模型方法的100 倍以上.當(dāng)數(shù)據(jù)量繼續(xù)增大時,計算量將迅速增大直至無法被接受,而本算法和硬邊界局部模型方法的模型更新時間處于同一數(shù)量級.另外,3 種方法對單個數(shù)據(jù)點的預(yù)測時間為同一數(shù)量級,其中本文算法的預(yù)測時間最短,為1.342 ms.綜上可見,本文提出的算法在保證高精度預(yù)測的同時,能大量減少局部模型的更新時間.

3.2 多維工程數(shù)據(jù)實驗

在機器人控制和航空航天等工程領(lǐng)域,輸入數(shù)據(jù)通常是多維的,離線學(xué)習(xí)的數(shù)據(jù)集很難覆蓋所有情景.即時能夠使用大量數(shù)據(jù)集提前進行離線學(xué)習(xí),計算量也會是巨大的.針對這些工程場景,本文提出的在線學(xué)習(xí)模型能夠使用實時數(shù)據(jù)更新模型,實現(xiàn)預(yù)測功能.

為了分析該模型處理真實多維工程數(shù)據(jù)時的性能,本文采用葡萄牙波爾圖大學(xué)公開的來自F16 型號飛機升降電梯的數(shù)據(jù)集進行實驗[26].F16 型號飛機升降電梯的數(shù)據(jù)集共包括9 517 個數(shù)據(jù)點.輸入數(shù)據(jù)有6 個屬性:爬升率、高度、轉(zhuǎn)輪角速度、轉(zhuǎn)輪角度、爬升加速度和爬升加速度的變化率,這6 個屬性的數(shù)值組成了6 維輸入數(shù)據(jù).輸出是控制信號Se,它是一個處于-0.014 到0.013 間的值.為了分析在線學(xué)習(xí)過程中,邊界約束回歸模型的性能是如何隨著數(shù)據(jù)量的增長變化的,將9 000 個數(shù)據(jù)點作為訓(xùn)練集,以數(shù)據(jù)流的形式不斷輸入模型進行更新.每增加1 000 個數(shù)據(jù)點,記錄當(dāng)前數(shù)據(jù)量對應(yīng)的更新時間、預(yù)測時間和誤差.這3 項是判斷模型是否滿足在線學(xué)習(xí)要求的關(guān)鍵指標(biāo),下面對這3 項指標(biāo)進行分析.

首先,用一個新數(shù)據(jù)點更新模型所需的時間Tu是衡量在線學(xué)習(xí)模型性能的重要指標(biāo),只有模型能夠隨數(shù)據(jù)流快速更新,才能實現(xiàn)在線學(xué)習(xí).圖5 展示了Tu隨數(shù)據(jù)量增長的變化趨勢.

圖5 更新時間隨數(shù)據(jù)量增長的變化趨勢Fig.5 The trend of update time as data increases

圖5 中N是用于訓(xùn)練模型的數(shù)據(jù)點的數(shù)量,m是每個維度的局部模型數(shù)量,通過m和每個維度輸入數(shù)據(jù)的范圍可以確定模型寬度.可以看到,當(dāng)m=3 時,Tu增長的趨勢過快,而m=5 時Tu由N=1 000 到N=9 000 只增加了0.75 倍,m=7 時Tu僅小幅度波動.因此m=5 或 7 對于此數(shù)據(jù)集是更為合適的參數(shù).

圖6 展示了對一個未知數(shù)據(jù)點進行預(yù)測所需的時間Tp隨數(shù)據(jù)量增長的變化趨勢.同樣顯示了m=5 或 7 時預(yù)測時間增長平緩,分別只增加了0.75 倍和0.28 倍,相比m=3 是更合適的參數(shù).

圖6 預(yù)測時間隨數(shù)據(jù)量增長的變化趨勢Fig.6 The trend of prediction time as the data increases

圖7 顯示了平均預(yù)測誤差隨數(shù)據(jù)量增長的變化趨勢.該圖的趨勢表明,隨著訓(xùn)練集中數(shù)據(jù)點的增加,模型的預(yù)測誤差減少.也就是說隨著數(shù)據(jù)流的輸入,在線學(xué)習(xí)模型能夠越來越精確.

圖7 預(yù)測誤差隨數(shù)據(jù)量增長的變化趨勢Fig.7 The trend of prediction error as data increases

綜合圖5～7,對于這一特定數(shù)據(jù)集m=5 是最合適的參數(shù).在m=5 時,數(shù)據(jù)點數(shù)量N=9 000雖然是N=1 000 時的9 倍,但它的更新時間0.541 ms只是N=1 000 時的0.309 ms 的1.75 倍,能夠?qū)崿F(xiàn)快速更新模型,滿足在線學(xué)習(xí)的要求.預(yù)測時間由N=1 000 時的2 ms 增加到N=9 000 時的3.5 ms,只增加了0.75 倍,同樣滿足在線學(xué)習(xí)的實時性要求.本模型的預(yù)測誤差由N=1 000 時的0.0000026989降低到N=9 000 時的0.0000020318,隨著數(shù)據(jù)量的增加實現(xiàn)了更高的精度.

為了進一步驗證本模型在機器人控制場景下的可用性,本文將模型應(yīng)用到具有兩個自由度的機器人的機械手臂上進行實驗,測試其在機器人控制方面應(yīng)用的可行性.該機器人手臂由文獻[27]中提出的基于PD 反饋的計算扭矩控制方法進行控制.我們?yōu)闄C器人手臂設(shè)置期望的二維平面活動軌跡,通過模型的在線學(xué)習(xí),預(yù)測控制機械臂所需的關(guān)節(jié)扭矩,并以此控制機器人手臂使得它能按照期望的活動軌跡運動.本文通過下式為機器人手臂設(shè)定期望的活動軌跡

其中,θ1和θ2分別是兩段機器人手臂與水平線的夾角,t是從0 開始計時的時間.將θ1和θ2對t求導(dǎo),得到兩段機器人手臂的角速度ω1和ω2.由θ1,θ2,ω1和ω2組成四維輸入矢量x=[θ1θ2ω1ω2].我們的任務(wù)是使用這些輸入數(shù)據(jù)通過本文模型預(yù)測機器人手臂的關(guān)節(jié)扭矩,并使用該關(guān)節(jié)扭矩控制機械臂,使得機械臂能夠跟隨預(yù)設(shè)軌跡,同時將四維輸入數(shù)據(jù)和對應(yīng)的真實關(guān)節(jié)扭矩輸入模型進行在線學(xué)習(xí).為了更清晰地展現(xiàn)控制效果,本文截取了機械臂10 s內(nèi)的運動軌跡參數(shù),共計2 000 個預(yù)測點,即每0.005 s 預(yù)測一次,對模型有實時性要求.

圖8 顯示了預(yù)設(shè)和實際運動軌跡分別對應(yīng)的參數(shù)隨時間的變化情況,其中,4 條由星號繪制的粗曲線分別對應(yīng)的是實際運動軌跡的角度和角速度,4條虛線分別對應(yīng)的是預(yù)設(shè)運動軌跡的角度和角速度.從圖中可以看出,在開始階段,實際運動軌跡和預(yù)設(shè)運動軌跡的曲線相差較大,但隨著模型的運行,兩者的曲線趨于重合.

圖8 預(yù)設(shè)和實際運動軌跡的參數(shù)Fig.8 Preset and actual motion trajectory parameters

為了更直觀地展示圖8 中預(yù)設(shè)和實際運動軌跡間的誤差變化情況,使用軌跡誤差矢量的范數(shù)||xr－x||表示實際運動軌跡和預(yù)設(shè)運動軌跡之間的誤差

其中,θr,1,θr,2,ωr,1,ωr,2是實際運動軌跡的參數(shù),θ1,θ2,ω1,ω2是預(yù)設(shè)運動軌跡的參數(shù).由此得到圖9,它顯示了運動軌跡誤差隨著時間振蕩下降的趨勢.圖8 和圖9 驗證了本文模型在高實時性要求的場景下的可用性.

圖9 運動軌跡誤差Fig.9 Motion track error

以上各實驗從不同角度驗證了本文模型能夠適用于工程中的大數(shù)據(jù)量場景,并且滿足在線學(xué)習(xí)的實時性需求.

4 結(jié)束語

為了在保持全局高斯過程回歸預(yù)測精度高等優(yōu)點的情況下,減小其立方復(fù)雜度帶來的過大的計算量,本文提出了一種基于邊界約束的概率相關(guān)局部高斯過程回歸模型.本文在該模型中使用一種基于貝葉斯原理的數(shù)據(jù)關(guān)聯(lián)方法,定義了用于數(shù)據(jù)點的分配和預(yù)測的激活函數(shù).同時,本文將局部模型放置在空間網(wǎng)格中,實現(xiàn)空間重疊的局部模型設(shè)置,提高了模型的精度.在算法方面,利用對計算過程的優(yōu)化和k-d 樹最近鄰搜索,以及更新模型預(yù)測矢量時采用已儲存數(shù)據(jù)求矩陣的逆,降低了模型更新階段的計算量.評估結(jié)果表明,與其他局部高斯過程回歸方法相比,本文方法可以顯著縮短計算時間,同時實現(xiàn)與全局高斯過程回歸方法相似的預(yù)測精度,能夠更好地處理大量且持續(xù)更新的數(shù)據(jù).在真實工程場景下,能夠很好地完成在線學(xué)習(xí)功能.該模型具有能夠快速更新和預(yù)測的特點,滿足在線學(xué)習(xí)的需求,對于機器人控制、航空航天等實時產(chǎn)生大量多維數(shù)據(jù)的領(lǐng)域,具有一定的工程意義.