王 紅

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

0 引 言

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

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

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

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

本文主要工作：

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

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

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

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

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

1 算 法

1.1 Kriging代理模型

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

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

(1)

(2)

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

fp=CF

(3)

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

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

(4)

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

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

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

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

3)新增數據點的優(yōu)化算法。新增數據點的優(yōu)化算法是基于Kriging代理模型優(yōu)化中討論最少的議題，針對非線性問題，通?？梢圆捎靡呀洿嬖诘娜魏蝺?yōu)化算法。

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

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

(5)

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

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

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

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

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

輸出：待增加點集合

過程:

1.計算目標函數相對于各參數的梯度。

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

3.判斷是否切換到SGD算法，否就跳轉到步驟4，是就計算SGD學習率后，跳轉到步驟5。

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

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

6.參數更新。

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

8.得到待增加點集合。

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

(10)

即Adam算法當前的學習率和平均學習率特別接近時，表示應該由Adam算法轉化為SGD算法。轉換之后SGD學習率采用式(11)計算，其中gt是梯度。

(11)

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

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

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

(12)

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

2 計算結果

2.1 算例介紹

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

(13)

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

(14)

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

表1 待估計參數信息

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

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

2.2 計算條件設定及結果說明

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

表2 4組樣例信息說明

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

(15)

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

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

(16)

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

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

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

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

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

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

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

表4 數據集4的本文算法估計得到的參數值

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

圖3 數據集4優(yōu)化過程所對應目標值變化過程

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

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

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

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

3 結束語

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