基于組合優(yōu)化的Kriging參數(shù)估計算法

王 紅

(大連東軟信息學(xué)院計算機學(xué)院,遼寧 大連 116023)

0 引 言

數(shù)學(xué)模型是理解復(fù)雜系統(tǒng)的結(jié)構(gòu)、功能，并對系統(tǒng)行為進(jìn)行設(shè)計、優(yōu)化和控制的重要工具[1-2]。實際生產(chǎn)中的很多數(shù)學(xué)模型包含了大量的參數(shù)，由于實驗條件及經(jīng)費的限制，很多參數(shù)不能從實驗中直接獲取，因此必須依靠優(yōu)化方法根據(jù)已有的實驗數(shù)據(jù)對這些參數(shù)進(jìn)行估計來獲得。當(dāng)模型具有強非線性以及sloppiness屬性(即模型的行為只由少數(shù)幾個剛性參數(shù)的組合便可決定)時[3]，參數(shù)估計過程中的大部分迭代過程都發(fā)生在目標(biāo)函數(shù)的平原區(qū)域，參數(shù)的更新可能把參數(shù)推向其取值范圍的邊界值，此時耗費大量的優(yōu)化時間但目標(biāo)函數(shù)值的變化卻很不明顯[4]。當(dāng)模型以常微分方程組(ODEs)形式描述時，由于模型響應(yīng)與待估參數(shù)之間的關(guān)系無法以顯示函數(shù)的形式表達(dá)，所以需要使用數(shù)值積分計算，而優(yōu)化過程中頻繁的積分運算會消耗大量的計算時間[5]，特別是當(dāng)常微分方程組具有剛性(stiff)特征時，求解時間可能會占據(jù)優(yōu)化總時間的95%[6]。

為了減少優(yōu)化過程中進(jìn)行積分計算所耗費的時間，可以通過一定數(shù)目的訓(xùn)練樣本構(gòu)建待估參數(shù)和目標(biāo)函數(shù)之間的近似關(guān)系，即所謂的代理模型[7]。這樣在參數(shù)估計的優(yōu)化過程中，用這些代理模型函數(shù)替代原有的目標(biāo)函數(shù)進(jìn)行迭代優(yōu)化，就不再需要計算數(shù)值積分過程，從而大大地節(jié)省了優(yōu)化時間。也就是說，基于代理模型進(jìn)行參數(shù)估計的常規(guī)思路是分為2個階段：1)代理模型構(gòu)建；2)基于代理模型進(jìn)行參數(shù)估計。

但僅依賴初始訓(xùn)練樣本構(gòu)建的代理模型很難達(dá)到預(yù)期精度要求。因此，需要依據(jù)加點準(zhǔn)則借助優(yōu)化算法查找增加新的樣例點，通過新增樣例點更新模型以提高代理模型精度[8-9]。如果采用的加點準(zhǔn)則能充分考慮目標(biāo)函數(shù)局部與全局的平衡關(guān)系，則用于精化代理模型的新增樣例點有可能就是目標(biāo)函數(shù)的全局最優(yōu)解。也就是說，在增加新樣例點精化代理模型的過程中同步完成了參數(shù)的優(yōu)化過程，即將前面的2個階段合成了一個階段。

在所有的代理模型中，Kriging代理可以等價于任何階的多項式函數(shù)，因此非常適用于求解具有多個極值的非線性函數(shù)插值問題[10]。本文基于Kriging代理模型估計ODEs類型函數(shù)的參數(shù)，以解決ODEs積分時間長的問題，同時組合Adam及SGD with momentum這2個算法的優(yōu)勢，搜索新樣例點更新代理模型，以提高最優(yōu)查找速度及收斂時間。

本文主要工作：

1)介紹Kriging代理模型的一般構(gòu)成及影響模型優(yōu)化結(jié)果的3個重要因素即初始樣例集合、加點準(zhǔn)則及新增樣例點的優(yōu)化算法。

2)介紹混合Adam&SGD的優(yōu)化算法。

3)介紹將混合Adam&SGD的優(yōu)化算法和Kriging代理模型結(jié)合的具體過程。

4)在計算結(jié)果部分首先介紹驗證本文算法的算例血管中性粒細(xì)胞模型，然后給出了本文提出的優(yōu)化算法對該模型中參數(shù)進(jìn)行估計的結(jié)果，并將該結(jié)果和其他優(yōu)化算法進(jìn)行比較。

5)總結(jié)概括了本文的算法及下一步的工作方向。

1 算 法

1.1 Kriging代理模型

Kriging代理模型[11]是基于回歸分析技術(shù)提出的，整個函數(shù)由代表全局趨勢函數(shù)和表征局部特點函數(shù)2個部分構(gòu)成，其中FT(x)βf為全局趨勢函數(shù)，Z(x)為局部特征函數(shù)。全局趨勢可以用常數(shù)或低階多項式表示，而局部特征函數(shù)Z(x)是一個動態(tài)隨機過程，服從式(2)所示的高斯分布。

f(x)=FT(x)βf+Z(x),x∈Rm

(1)

(2)

式(1)中，F(xiàn)(x)是不同取樣點x所對應(yīng)的目標(biāo)值所構(gòu)成的集合，βf是和F(x)對應(yīng)的因子向量。式(2)中，R是相關(guān)函數(shù),根據(jù)式(2)所示，在樣例集合S=[x1,x2,…,xn]及對應(yīng)的目標(biāo)值集合F=[f1,f2,…,fn]已知的情況下，任給數(shù)據(jù)點xnew，根據(jù)已構(gòu)建的代理模型，可知該點的目標(biāo)預(yù)測值為式(3):

fp=CF

(3)

其中，C是權(quán)重系數(shù)，代表樣例集合中每個點對新點xnew的貢獻(xiàn)度。根據(jù)Kriging的估計具有無偏性及最優(yōu)性2個特征，可推導(dǎo)出期望為0時，預(yù)測值和真實值之間的差平方和最小。據(jù)此可將式(3)細(xì)化為式(4)，其中σ2是預(yù)測值與真實值之間的方差估計。

fp=σ2-λT(FTC-f)

(4)

通用Kriging的優(yōu)化算法流程如圖1所示。在圖1所示的流程里，對最終優(yōu)化結(jié)果起到影響作用的因素包括3個：

圖1 常用Kriging優(yōu)化算法流程

1)構(gòu)建初始代理模型的樣例點集。其中初始樣例點集合的好壞不僅決定了模型的優(yōu)劣還直接影響了后續(xù)模型尋找最優(yōu)值的優(yōu)化速度，因此要求初始集合點應(yīng)該能提供足夠多的信息量，即在對參數(shù)空間沒有任何先驗知識的前提下盡量均勻取樣。

2)加點準(zhǔn)則。加點準(zhǔn)則是選擇新樣例點的基本標(biāo)準(zhǔn)，該標(biāo)準(zhǔn)就是優(yōu)化算法的目標(biāo)函數(shù)。作為一種基于高斯過程的代理模型，Kriging代理可以在優(yōu)化過程中確定模型的不確定性，因此能預(yù)知在參數(shù)空間的哪些區(qū)域增加新的樣例點可以提高構(gòu)建模型的精度。常見的加點準(zhǔn)則有可能提高加點準(zhǔn)則(PI)[12]、期望提高(EI)[13]和廣義期望提高(GEI)[14]等。

3)新增數(shù)據(jù)點的優(yōu)化算法。新增數(shù)據(jù)點的優(yōu)化算法是基于Kriging代理模型優(yōu)化中討論最少的議題，針對非線性問題，通?？梢圆捎靡呀?jīng)存在的任何優(yōu)化算法。

1.2 混合Adam&SGD的優(yōu)化算法

所謂參數(shù)估計就是借助于某種優(yōu)化算法，通過最小化目標(biāo)函數(shù)找到一組參數(shù)值p，以和實驗數(shù)據(jù)最好地擬合。以常微分方程組形式呈現(xiàn)的參數(shù)估計問題可以表示為：

(5)

其中，x代表ODEs，p是待估計參數(shù)，f是ODEs右側(cè)函數(shù)。xei(tj)代表tj時刻實驗值，xi(tj，p)代表tj時刻計算值，N是樣例點數(shù)目，n是ODEs的數(shù)目，PU和pL代表參數(shù)取值范圍的上下限。由于f的非線性特性，因此x很難有解析解。

機器學(xué)習(xí)也是一種代理模型，其中一系列的優(yōu)化算法各具特色。動態(tài)梯度下降法(SGD)[15]是其中最簡單的優(yōu)化算法，又稱小批量梯度下降法。通過計算目標(biāo)函數(shù)對已有參數(shù)的梯度方向，并將其作為參數(shù)更新下降的最快移動方向，公式表示為：

(6)

這里的α是尺度系數(shù)，決定了在最大下降方向上的步長。mt和Vt分別代表一階動量和二階動量，而gt是目標(biāo)函數(shù)相對于參數(shù)的梯度。該增量通用公式適用于機器學(xué)習(xí)中的所有優(yōu)化算法[16]。

對目標(biāo)函數(shù)比較簡單的優(yōu)化問題，梯度下降法收斂效果比較好，而對于參數(shù)以千計算的優(yōu)化問題來說，梯度下降法可能是唯一的選擇[17]。但該方法在遇到V型搜索區(qū)域時，會在波谷兩側(cè)遞增持續(xù)震蕩，浪費大量的迭代時間，最后可能還會停留在一個局部最優(yōu)點。針對SGD在峰谷震蕩的情況，引入Nesterov提出的梯度加速法[18],通過修改梯度的計算公式為式(7)，將未來即將走到的點梯度方向與歷史累積動量相結(jié)合，以達(dá)到預(yù)看多一步的目的。

(7)

同時為提高下坡時的優(yōu)化速度，可在梯度下降過程中加入慣性[19]，此時一階動量修改為式(8)所示，使得當(dāng)前點的梯度下降方向由累積下降方向與當(dāng)前下降方向共同決定。

mt=β1·mt-1+(1-β1)·gt

(8)

以上SGD及其變形都是針對原有算法中存在的問題，基于一階動量提出的改進(jìn)。為提高算法的收斂速度以及使其具有自適應(yīng)的特點，Adam算法[2]在步長更新中引入二階動量即梯度方差元素，修改如式(9)：

(9)

該算法適用于處理大型數(shù)據(jù)及稀疏數(shù)據(jù)，但缺點是二階動量Vt不是單調(diào)變化，因此最終可能不收斂，而且可能會錯過全局最優(yōu)解[20]。

考慮一階動量和二階動量的各自優(yōu)勢，Keskar等人[21]提出可以將二者組合，同時享有Adam優(yōu)化速度快和SGD最終收斂的雙重優(yōu)點。本文利用該思路，將其作為尋找新增點的優(yōu)化算法。算法流程如算法1所示。

算法1 本文尋找新增點的優(yōu)化算法的流程

輸入：設(shè)定優(yōu)化過程所需初始參數(shù)

輸出：待增加點集合

過程:

1.計算目標(biāo)函數(shù)相對于各參數(shù)的梯度。

2.每循環(huán)10次按照特定規(guī)律減小學(xué)習(xí)率。

3.判斷是否切換到SGD算法，否就跳轉(zhuǎn)到步驟4，是就計算SGD學(xué)習(xí)率后，跳轉(zhuǎn)到步驟5。

4.執(zhí)行Adam算法，得到參數(shù)變化量。

5.執(zhí)行SGD算法，得到參數(shù)變化量。

6.參數(shù)更新。

7.優(yōu)化算法是否收斂，如收斂跳轉(zhuǎn)到步驟8，否則跳轉(zhuǎn)到步驟1。

8.得到待增加點集合。

其中算法Adam和SGD切換條件需滿足式(10)：

(10)

即Adam算法當(dāng)前的學(xué)習(xí)率和平均學(xué)習(xí)率特別接近時，表示應(yīng)該由Adam算法轉(zhuǎn)化為SGD算法。轉(zhuǎn)換之后SGD學(xué)習(xí)率采用式(11)計算，其中g(shù)t是梯度。

(11)

1.3 混合Adam&SGD優(yōu)化算法的Kriging代理模型

使用混合Adam&SGD優(yōu)化算法的Kriging代理模型流程如圖2所示。

圖2 混合Adam&SGD優(yōu)化算法的Kriging流程

(12)

為提高建模速度，這里采用多點加點原則，每次增加3個點，即GEI值最大的點、實際值與預(yù)測值相差最大的點、目標(biāo)函數(shù)值最小的點用來充實模型構(gòu)建樣例點集。

2 計算結(jié)果

2.1 算例介紹

本文使用血管中性粒細(xì)胞(PMN)[22]代謝模型作為檢測算法有效性的算例。該模型使用13個常微分方程組表示細(xì)胞中的酶促反應(yīng)動力學(xué)特征即米氏方程，該方程具體如式(13)和式(14)：

(13)

其中,[S]是底物濃度，[Et]是酶的整個濃度，Kcat是酶的轉(zhuǎn)換數(shù)，Km是米氏常數(shù)。而加入可逆抑制劑后式(13)可轉(zhuǎn)化為式(14)。其中[I]是抑制劑的濃度，KI是抑制常數(shù)。

(14)

在13個方程中共含有41個參數(shù)。根據(jù)這些參數(shù)在酶促反應(yīng)中承擔(dān)的角色不同將其分為8個類別，不同類別的取值范圍從跨越2個數(shù)量級到跨越6個數(shù)量級變化不等，具體信息如表1所示。

表1 待估計參數(shù)信息

因為這41個參數(shù)的值無法通過實驗獲取到，只能通過優(yōu)化算法對其進(jìn)行數(shù)值估計。本文使用基于混合Adam&SGD優(yōu)化算法的Kriging代理模型對這8個類別中的41個參數(shù)值進(jìn)行估計。

因為沒有實驗測量得到的參數(shù)真實數(shù)值作為比對，因此參數(shù)估計結(jié)果的優(yōu)劣只能通過目標(biāo)函數(shù)值的大小進(jìn)行衡量，后續(xù)的結(jié)果說明就是基于這樣的前提完成的。

2.2 計算條件設(shè)定及結(jié)果說明

初始樣例采用確定性均勻分布與隨機性均勻分布相結(jié)合方式來生成，考慮PMN模型中有的參數(shù)變化跨越幾個指數(shù)級別，因此取對數(shù)后均勻取樣。抽取4組樣例，樣例信息如表2所示。

表2 4組樣例信息說明

算法優(yōu)化終止條件有3個：1)滿足式(15)；2)迭代次數(shù)≥30；3)連續(xù)3次無新樣例點加入。其中終止條件中的第2個是為了體現(xiàn)資源受限的特征而設(shè)置的。

(15)

算法開發(fā)平臺是Matlab2009b，CPU是Intel(R) Core(TM) 4.2 GHz。

因為PMN模型的高度非線性所以很難有解析解，在求解雅可比矩陣時需要使用數(shù)值形式計算多次積分，可能會增加積分時間。為節(jié)省計算時間，本文利用Broyden[23]提出的擬牛頓根發(fā)現(xiàn)法來估計雅可比矩陣。該方法最大的特點就是依賴前次矩陣值可以計算出本次矩陣值，具體如式(16)：

(16)

其中，ΔOi是目標(biāo)函數(shù)變化，Δpi是參數(shù)值變化。為避免累計誤差，算法每5次迭代執(zhí)行后重新計算該矩陣值。

為驗證本文優(yōu)化算法對Kriging代理模型的有效性，本文選擇SSMGO(SS)[24]、Matlab工具箱自帶的fmincon (FM)、SGD with Momentum(SM)[19]及Adam(AD)[2]這4種算法與本文算法進(jìn)行比較。針對這5種算法，在4組集合上各執(zhí)行3次。對每種算法的運行結(jié)果，挑選優(yōu)化中目標(biāo)函數(shù)最小的那組作比較，結(jié)果如表3所示(其中初值指初始目標(biāo)值，終值指優(yōu)化后目標(biāo)值)。

表3 各優(yōu)化算法結(jié)果比較

從表3給出的5種算法比較結(jié)果可以得到如下結(jié)論：

1)從最優(yōu)解角度看:建模數(shù)據(jù)集中樣例多的優(yōu)化效果整體好于樣例少的；初始數(shù)據(jù)集中的最小值對某些優(yōu)化結(jié)果有一定影響，但并未看出直接影響關(guān)系。大部分情況下本文算法和SSMGO能得到相對好的優(yōu)化效果。

2)從優(yōu)化時間看:在同等條件下，數(shù)據(jù)集中樣例數(shù)目越少，優(yōu)化時間也相對短；對5種算法而言，相同數(shù)據(jù)集條件下，所用時間在同一數(shù)量級。本文算法和SSMGO算法所用時間略長于其他優(yōu)化算法。

以表3中數(shù)據(jù)集4為例，給出本文算法估計出的41個參數(shù)結(jié)果值如表4所示。

表4 數(shù)據(jù)集4的本文算法估計得到的參數(shù)值

以表3中數(shù)據(jù)集4為例，本文展示5種優(yōu)化算法在優(yōu)化過程中目標(biāo)值的變化過程，如圖3所示。

圖3 數(shù)據(jù)集4優(yōu)化過程所對應(yīng)目標(biāo)值變化過程

從圖3的變化過程可見，本文算法和SW在規(guī)定的優(yōu)化次數(shù)內(nèi)有收斂趨勢。所有5種優(yōu)化算法中，F(xiàn)M和AD最先達(dá)到該算法的優(yōu)化最小值，而SW后達(dá)到優(yōu)化最小值。所有算法都有震蕩過程，即目標(biāo)值升高后有一個下降過程。

仍以數(shù)據(jù)集4為例，本文展現(xiàn)不同初始學(xué)習(xí)率對本文算法優(yōu)化的影響，過程如圖4所示。

圖4 不同初始學(xué)習(xí)率對優(yōu)化過程的影響

由圖4可見2種初始學(xué)習(xí)率，在優(yōu)化過程中均有抖動，但在后續(xù)均趨于平穩(wěn)。低的初始學(xué)習(xí)率比高的初始學(xué)習(xí)率在優(yōu)化過程中的下降過程要緩慢一些。

3 結(jié)束語

本文針對Kriging代理方法在解決非線性且具有sloppiness屬性模型參數(shù)估計問題中的新增點查找算法進(jìn)行改進(jìn)，將Adam算法和帶有動量的SGD算法的優(yōu)勢相結(jié)合，以期達(dá)到提高搜索質(zhì)量以及最后收斂的雙重目的。將該算法應(yīng)用于PMN代謝模型，并與其他4種算法進(jìn)行了比較。結(jié)果顯示，在可比擬的時間內(nèi)該算法能得到較好的優(yōu)化效果。下一步的工作是探索學(xué)習(xí)率的修改規(guī)則以及進(jìn)入局部最優(yōu)的調(diào)整規(guī)則。