基于大偏差理論非高斯隨機動力系統(tǒng)離出行為研究

李 揚 趙 鋒 劉先斌 ,*

1 南京理工大學自動化學院, 南京 210094

2 南京航空航天大學機械結構力學及控制國家重點實驗室, 南京 210016

1 引 言

與非線性一樣, 隨機性也是客觀世界的本質屬性. 客觀世界豐富多彩的復雜現象, 其成因通常是源于系統(tǒng)內稟的非線性、隨機性及其之間的相互作用. 由于噪聲樣本的多樣性及不確定性,受隨機擾動的非線性系統(tǒng)所呈現的與原確定性系統(tǒng)截然不同的復雜且有趣的隨機動力學行為一直是人們關注的焦點! 而其中, 受弱噪聲擾動的系統(tǒng)動力學行為則更為有趣: 對于未受擾系統(tǒng)與受擾系統(tǒng), 一方面, 在噪聲充分弱的前提下, 人們會自然期待后者行為與前者充分接近 (依概率、概率1、均方意義等); 另一方面, 前者是確定性的, 其動力學行為由初始條件唯一確定, 而弱噪聲是否會導致后者長期的隨機動力學行為完全不同則更值得關注! 對此, 大偏差原理及其與之密切相關的離出問題研究給出了一個很好的回答. 與大數定律 (law of large numbers) 和中心極限定理 (central limit theorem) 不同, 這兩者只關心大概率事件及其極限行為, “大偏差理論 (large deviation theory)” (Freidlin & Wentzell 2012) 則試圖揭示“非線性隨機系統(tǒng)中最初的小概率事件經長時間演化可成為大概率事件” 的機理. 這是因為受微擾系統(tǒng)行為對未受擾系統(tǒng)行為的大偏差在有限的時間內是小概率事件. 大偏差現象產生的原因在于, 概率理論中, 所謂一個“弱”隨機擾動, 指的是其統(tǒng)計意義, 即: 擾動幅值很小的概率很大, 擾動幅值很大的概率很小, 且擾動的均值很小 (陳朕 2018).

數學上,粗略地說,大偏差理論源于一個事實:一族由實參數化的概率測度對于可測事件的估值呈現參數的指數變化率. 理論起源可追溯到北歐的精算師對保險行業(yè)風險估計研究——Cramér (1944) 的關于獨立隨機變量和一般的大偏差結果. Sanov (1957) 推導了獨立同分布隨機變量的經驗分布. 隨后, Donsker 和Varadhan (1975a, 1975b, 1976, 1983), G?rtner (1977) 將其推廣至一般的Markov 鏈和Markov 過程. 而動力系統(tǒng)受弱隨機擾動所產生的樣本路徑大偏差的結果最初由Varadhan (1967)、Ventsel 和Freidlin (1970)、Freidlin 和Wentzell (2012) 給出. 上述研究均基于一種廣義Cramér 變換, 即通過選取參考測度, 使其對于所考察的事件具有大的概率測度,然后對于所考察的概率族建立關于參考測度的Radon-Nikodym 導數, 并進行估值 (Freidlin &Wentzell 2012; Wentzell 2012; Varadhan 1984, 2016; Deuschel & Stroock 2001; Dembo & Zeitouni 1998; Feng & Kurtz 2006). 20 世紀90 年代, Puhalskii, O’Brien 和Vervaat, Acosta 等通過類似于“基于Prohorov 緊性證明一族概率測度弱收斂”的方法, 將大偏差理論推廣至樣本路徑左連續(xù)且具有右極限的隨機過程類 (Varadhan 1984, 2016; Deuschel & Stroock 2001; Dembo & Zeitouni 1998; Feng & Kurtz 2006).

在隨機動力系統(tǒng)中, 噪聲誘導的大偏差現象會導致系統(tǒng)遠離平衡態(tài)的復雜動力學行為, 即從穩(wěn)態(tài)的吸引域中的離出現象. 離出現象是噪聲誘導的非線性系統(tǒng)的全局動力學行為, 是非線性隨機系統(tǒng)所特有的復雜現象! 在弱隨機擾動下, 即使一個結構穩(wěn)定的動力系統(tǒng), 其受擾后的軌線可與未擾系統(tǒng)軌線之間產生大偏差并因此失去“結構穩(wěn)定性” ! 一個即使只有漸近穩(wěn)定平衡點的動力系統(tǒng), 在非常弱的隨機擾動下, 其長期作用的累積效應可使其樣本軌線具有從吸引域中概率1的離出概率且可以到達任何有界的區(qū)域 (朱金杰等 2020, 劉先斌等 1996). 數學嚴格意義上, 反映確定性動力系統(tǒng)長期行為的非游蕩點集、平衡點和極限環(huán)的穩(wěn)定性、全局結構等概念已不再適用于隨機系統(tǒng), 因為此時相空間結構可能已面目全非. 而對系統(tǒng)隨機動力學行為“如何偏差”以及“偏差發(fā)生的統(tǒng)計特征”等問題的探索便自然成為大偏差理論研究意義之所在 (Freidlin &Wentzell 2012, 朱金杰等 2020, 李揚 2021).

在離出問題中, 最優(yōu)離出路徑 (optimal exit path) 或稱最大可能離出路徑 (most probable exit path, MPEP)、平均首次離出時間 (mean first passage time, MFPT) 和離出位置分布 (exit location distribution, ELD) 是刻畫系統(tǒng)“怎樣離出” “何時離出”和“從哪離出”的三個重要的特征量,對它們的計算和分析一直是隨機動力學領域的核心問題. Freidlin 和Wentzell (2012) 建立的大偏差理論從數學上為這個問題的解答提供了理論基礎, 其核心思想在于, 定義了作用量泛函 (action functional) 的概念, 導出了軌道鄰域小管道的概率的指數估計式, 從而將稀有事件的概率計算問題轉化為泛函的變分問題. 基于經典力學的結果, 推導了最優(yōu)離出路徑滿足的歐拉-拉朗格朗日方程 (Euler-Lagrange equation), 通過定義作用量泛函最小值為擬勢 (quasi-potential), 同時證明了平均離出時間和離出位置分布與擬勢的指數關系.

對于Markov 過程, MFPT 是停時 (stopping time) (Freidlin & Wentzell 2012, Bernt 2010)的均值. 它在諸多學科領域中都有其實際的應用價值, 如動力系統(tǒng)中衡量吸引子之間的相對穩(wěn)定性(K?osek-Dygas et al. 1988)、神經元自發(fā)性放電現象的峰峰間期 (interspike interval) (朱金杰2018, Ryashko 2018)、工程領域中結構的可靠性分析 (Chen & Zhu 2009, 2010)、生物學中基因的轉錄與翻譯 (Bressloff 2017) 等. 針對MFPT 的研究由來已久, 最早可以追溯到20 世紀40 年代,Kramers (1940) 研究了朗之萬方程 (Langevin equations) 的離出問題, 發(fā)現MFPT 與噪聲強度的指數依賴關系. 針對擬可積哈密頓系統(tǒng), 朱位秋等 (1992, 2017) 利用隨機平均法研究了首次穿越時間和首次穿越失效的問題, 并以此刻畫系統(tǒng)的可靠性. 此外, 胞映射作為一種有效的全局動力學行為分析方法, 也是計算MFPT 的有效數值工具 (徐偉等 2013). 另外, 許多研究者發(fā)現, 當MFPT 與系統(tǒng)內稟的或外在的時間尺度匹配時, 系統(tǒng)會產生類似于共振的行為, 如隨機共振和自誘導隨機共振等現象, 在機械系統(tǒng)和醫(yī)學等領域都有重要的應用 (Kang et al. 2017, Zhang et al.2021, 王煒 2020, 顏志2021).

對于受高斯白噪聲擾動的非線性系統(tǒng), MFPT 方程v(x)=Exτ滿足如下Pontryagin 方程 (狄利克雷邊值問題) (Bernt 2010, Duan 2015)

其中A為生成元 (即Fokker-Planck 算子的伴隨算子),D為待離出區(qū)域. 利用奇異攝動法和特征線法, Matkowsky 和Schuss 等 (1982, 1983, ) 及Naeh 等(1990)解析近似計算了在特征邊界 (characteristic boundary) 和非特征邊界 (non-characteristic boundary) 條件下的MFPT 漸進展開式. 基于van der Pol 變換、奇異攝動法和隨機平均法, 孔琛和劉先斌 (2014) 計算了一類分段線性系統(tǒng)的MFPT, 發(fā)現分形邊界會造成理論與數值結果的偏差. 之后, Kong 和Liu (2017) 利用廣義圖胞映射方法進行全局動力學分析以及結合隨機平均法降維, 研究了受弱高斯白噪聲擾動的分段線性系統(tǒng)的噪聲誘導的混沌現象, 發(fā)現MFPT 可以作為“噪聲誘導混沌”的指示器.

離出位置分布刻畫了系統(tǒng)跨越邊界時在邊界上的離出點的分布情況, 同樣是一個重要的物理量.對任意的連續(xù)函數f(x),函數u(x)=f(y)pe(x,y)dSy(pe(x,y)為離出位置分布) 滿足方程

在弱高斯噪聲條件下, Matkowsky 和Schuss 等 (1982, 1983,) 及Naeh 等(1990)利用奇異攝動法推導了離出位置分布的解析近似表達式. 針對高斯白噪聲擾動的FitzHugh-Nagumo 系統(tǒng),Chen 等 (2017) 發(fā)現模擬軌道離出位置與邊界擬勢最低點的偏差, 即鞍點回避現象(saddle point avoidance), 本質上來源于有限噪聲強度的影響. 對于另一個可激系統(tǒng), Li 等 (2020b) 同樣揭示了鞍點回避現象, 并利用WKB 近似攝動計算了離出位置分布, 發(fā)現一階修正能夠彌補這個偏差.由于方程 (2) 的求解較為困難, 許多研究者設計了一系列計算離出位置分布的數值方法. 例如, 廣義胞映射方法作為分析非線性系統(tǒng)全局動力學結構的有效工具, 同樣能夠用來計算離出位置分布 (Han et al. 2016). Zhu 等 (2018) 綜合了正向通量采樣 (forward flux sampling) (Allen et al.2005) 和快速Monte Carlo 模擬方法的優(yōu)點, 提出了概率演化算法 (probability evolution method),能使用模擬方法快速得到離出位置分布.

最優(yōu)離出路徑的計算是最后也是最難的問題, 它旨在刻畫離出發(fā)生的內在機理, 具有結構性的難度, 是大偏差理論的核心問題. 基于作用量泛函的變分, 最優(yōu)離出路徑滿足歐拉-拉格朗日方程或對應的哈密頓系統(tǒng). 此時, 哈密頓相空間相對于原確定性相空間擴展一倍. 若考慮具有穩(wěn)定不動點的二維系統(tǒng), 則確定性相空間作為不動點的二維穩(wěn)定流形嵌入到四維拓展哈密頓相空間中. 由于噪聲的加入, 不動點多出兩個不穩(wěn)定方向, 張成了二維不穩(wěn)定流形, 即拉格朗日流形(Lagrangian manifold), 此流形到確定性相空間的投影構成離出路徑的圖案. 由于最優(yōu)離出路徑的解析計算比較困難, 眾多研究者發(fā)展了大量的數值方法, 如action plot (Beri et al. 2005)、幾何最小作用量方法 (geometric minimum action method, GMAM) (Heymann & Vanden-Eijnden 2008)、有序逆風法 (ordered upwind method, OUM) (Cameron 2012)、歷程概率分布 (prehistory probability distribution) (Dykman et al. 1992) 和機器學習 (Machine Learning) (Li et al. 2021a) 等方法.

一般情況下, 若終點坐標在相空間中連續(xù)變化, 人們自然認為連接始末端點的最優(yōu)路徑也隨之連續(xù)變化. 然而, 對于非平衡態(tài)系統(tǒng), 情況卻未必如此. Chen 和Liu (2016) 通過對周期激勵雙勢阱系統(tǒng)的離出問題的研究, 發(fā)現最優(yōu)路徑可能隨著終點的變化發(fā)生突變, 出現這種現象的根本原因在于拉格朗日流形出現折疊結構, 導致泛函極值的多值性和離出軌跡圖案的奇異性. 在折疊的邊緣——焦散線 (caustic) 處, 擬勢的二階導數發(fā)散, WKB 近似指數前因子也因此發(fā)散 (孔琛2018). 在折疊區(qū)域有一條明顯的切換線 (switching line), 切換線兩側即使相距很近的兩個點, 會有兩條拓撲性質完全不同的最優(yōu)路徑到達. 在切換線上, 擬勢不可微. Chen 等 (2019) 研究了Type-I 和Type-II 可激性條件下的Morris-Lecar 系統(tǒng), 揭示了擬勢的非可微性、最優(yōu)路徑的非光滑性和拉格朗日流形的奇異性的聯系, 同時發(fā)現最優(yōu)路徑滿足的非光滑向量場會出現閉軌, 本質上是非光滑系統(tǒng)的擦切分岔所產生的擦切環(huán), 或由滑移分岔產生的滑移環(huán).

另外, 離出問題中某些奇異現象的出現來源于動力系統(tǒng)本身的復雜結構, 如混沌現象和分形邊界等. 針對非雙曲混沌吸引子的離出問題, 最新研究 (Chen et al. 2016) 發(fā)現其離出過程是沿著一個異宿軌道的交叉層次結構進行的, 噪聲不斷地將離出軌跡從當前鞍型周期軌道的不穩(wěn)定流形擾至包含其上一層不穩(wěn)定流形, 而離出軌跡在不穩(wěn)定流形上的每一次躍遷都伴隨著噪聲的一次劇烈波動. Kraut 和Feudel (2003) 研究了吸引域中存在混沌鞍的Ikeda 映射的離出問題, 發(fā)現在混沌鞍上擬勢幾乎水平, 由此極大降低了離出所需的能量. 還有一些學者研究了具有局部不連通 (locally disconnected) 和局部連通 (locally connected) 分形邊界的連續(xù)和離散動力系統(tǒng)的離出問題, 發(fā)現離出總是發(fā)生于唯一的可達邊界點 (accessible boundary point) 上 (Silchenko et al.2005).

目前, 除了涉及混沌系統(tǒng)的復雜離出機理和離出路徑的拓撲奇異性問題之外, 有關高斯噪聲離出問題的結果相對較為完善. 而關于非高斯隨機系統(tǒng)的研究則較為匱乏, 近年來逐漸成為離出問題研究的焦點之一. 事實上, 相比高斯擾動, 自然界中的大量涉及跳躍、突變、間歇、爆破、躍遷等客觀現象, 如基因突變、火山爆發(fā)、能級躍遷等, 更適合利用具有非高斯噪聲的隨機動力系統(tǒng)建立模型 (Duan 2015, 李揚2021). 因此, 非高斯隨機動力系統(tǒng)是具有深刻科學意義的現象學模型, 相關的離出問題的研究不僅對于隨機動力系統(tǒng)理論具有推動作用, 更有著重要的實際應用價值.

非高斯擾動往往涉及噪聲樣本的跳躍特性, 目前研究較多和應用廣泛的非高斯隨機系統(tǒng)主要有三種. 第一種是隨機混合系統(tǒng) (stochastic hybrid system), 即連續(xù)變量與離散變量耦合的過程. 控制連續(xù)變量演化的方程的向量場依賴于離散變量, 離散變量是轉移率依賴于連續(xù)變量的跳躍Markov 過程, 因此兩者構成耦合隨機過程. 第二種是具有指數輕跳躍擾動 (exponentially light jump fluctuation) 的隨機動力系統(tǒng), 這種噪聲是Lévy 噪聲的一種, 其跳躍測度 (jump measure) 指數衰減. 第三種非高斯噪聲也是Lévy 噪聲的一種, 即α穩(wěn)定Lévy 噪聲, 其跳躍測度多項式衰減.本文第二節(jié)引入大偏差理論的基本思想, 第三、四、五節(jié)分別介紹三種常見非高斯隨機系統(tǒng)的離出問題的主要處理方法和近期研究進展, 并在最后一節(jié)提出一些開放性問題.

2 大偏差理論構建

2.1 大偏差理論基本思想

在非線性動力學中, 人們往往更加關注極限行為, 如平衡點、極限環(huán)和奇怪吸引子等不變集結構. 類似地, 眾所周知的大數定律和中心極限定理是概率理論中有關極限行為的兩個主要結果. 直觀上, 對于類型的隨機動力系統(tǒng),ξt為給定的隨機擾動,ε為小的噪聲強度,ε=0對應確定性系統(tǒng)t=b(xt),大數定律給出了擾動系統(tǒng)的解在有限時間內對于未受擾軌道xt的依概率收斂性結果, 中心極限定理則進一步斷言一階近似依分布收斂于正態(tài)分布.

與大數定律和中心極限定理不同, 這兩者只關心未受擾軌道附近的大概率事件, 而大偏差理論則反其道而行之, 它關注稀有事件 (rare events), 如噪聲誘導的從吸引域中的逃逸現象 (Maier &Stein 1992, Li & Liu 2019)、噪聲誘導的不同穩(wěn)態(tài)之間的轉移 (Li et al. 2020a, Chen & Liu 2017)、神經元系統(tǒng)跨閾值放電現象 (Khovanov et al. 2013, Li et al. 2020b, Chen et al. 2017)、噪聲誘導的混沌 (Tamás & Lai 2010, Kong &Liu 2017) 等, 均屬大偏差理論的應用范疇. 在有限時間內, 大偏差現象發(fā)生的概率很小, 這樣深度挖掘極度稀有的事件似乎與概率論的宗旨相悖, 因為似乎大概率事件對預測事物演化更具啟發(fā)式意義. 然而, 若將時間尺度延至無窮, 小概率事件長期的累積效應將可能逐漸演化為大概率事件. 因此, 在有限時間內精細區(qū)分出哪些事件是“更不可能的(more improbable)”, 哪些事件是“稍不可能的 (less improbable)”, 成為揭示無限時間尺度下概率1 的大偏差現象本質機理的一把鑰匙.

數學上, Freidlin-Wentzell (2012) 大偏差理論建立了隨機動力系統(tǒng)的解穿過一個給定路徑φt的鄰域的概率估計式

其中S0T[φ]=即為引言中提到的作用量泛函,L為拉格朗日量, 其形式取決于所研究的隨機過程 (擴散過程、隨機混合系統(tǒng)等), sup|·|記為上確界范數,δ是一個正的小參數, Px的下標為初值=φ0=x. 由于單條路徑的概率必然為零, 為了對比不同路徑實現的可能性, 式 (3)利用軌道附近小鄰域的概率代替原軌道的概率. 根據Laplace 方法, 諸多事件的概率均取決于此泛函的約束最小化, 例如B為Rn的某一Borel 子集, 則∈B的概率對數等價于exp( inf表示下確界). 因此, 大偏差理論的精辟之處在于, 它將概率的復雜計算問題轉化為作用量泛函的約束最小化問題, 從而可以利用分析力學和變分學的經典結果進行求解.

根據變分原理, 作用量泛函的駐值解-最優(yōu)離出路徑, 滿足歐拉-拉格朗日方程

或對應的輔助哈密頓系統(tǒng)

其中哈密頓量H是L的Legendre 變換. 盡管離出路徑的眾多樣本是隨機的, 然而它以壓倒性的概率落在最優(yōu)路徑附近, 且ε →0時, 它依概率收斂于最優(yōu)路徑. 大偏差理論的結果使得隨機動力學的某些復雜行為變得近乎可預測! 在離出問題中, 最優(yōu)路徑的研究始終是一個結構性的難題. 分析最優(yōu)路徑的結構對于理解噪聲在離出問題中的作用至關重要, 是隨機動力學的核心問題之一.

2.2 擬勢

如引言所述, 平均離出時間和離出位置分布均指數依賴于作用量泛函的最小值, Freidlin 和Wentzell (2012) 將其定義為擬勢. 具體來說, 擬勢具有如下形式

其中φ ∈(0,T)表示φ(0)=,φ(T)=x,記為絕對連續(xù)函數集. 換句話說, 擬勢是連接初始位置和終點位置x之間所有絕對連續(xù)函數的作用量的最小值, 且要求對軌跡和時間都取最小. 通常, 初始位置選為穩(wěn)定不變集 (如穩(wěn)定不動點O), 此時可將擬勢V(O,x)簡記為V(x). 對于多穩(wěn)態(tài)的系統(tǒng), 需要先計算出每個穩(wěn)態(tài)的吸引域擬勢, 然后將幾個區(qū)域的擬勢拼接起來, 得到全局范圍的擬勢.

針對簡單的加性高斯噪聲的情況, 擬勢的意義可以按下述方式理解. 首先, 假定向量場b(x)具有如下的分解

其中連續(xù)可微函數U(x)滿足U(O)=0(O為D中的穩(wěn)定平衡點), 對x/=0有U(x)>0且?U(x)/=0,(l(x),?U(x))=0.則系統(tǒng)關于點O的擬勢V(x)滿足V(x)=2U(x) (Freidlin & Wentzell 2012). 若U(x)是二次連續(xù)可微的, 則從O到x的最大可能路徑φs,-∞≤s≤T, 滿足下面的方程

對比方程 (8) 和原始確定性方程后發(fā)現, 在噪聲的驅動下, 離出路徑保持旋轉分量l(x)不變,僅僅將向內的法向分量-?U(x)翻轉向外, 使得噪聲能夠消耗盡可能小的能量驅動系統(tǒng)離出. 若考慮有勢系統(tǒng), 即旋轉分量l(x)=0,則擬勢與勢函數U(x)僅僅相差一個常數因子. 換句話說, 在有勢系統(tǒng)中, 擬勢與勢函數完全等價. 即使在非有勢系統(tǒng)中, 擬勢仍然可以看成是廣義的勢函數.這也是擬勢這個名稱的來源.

通常, 對向量場進行線性化并計算特征值或最大李雅普諾夫指數 (largest Lyapunov exponent) 是最常見的系統(tǒng)穩(wěn)定性分析方法. 特征值取決于勢函數的二階導數, 刻畫了系統(tǒng)離開不動點鄰域的難易程度, 因此可以理解為系統(tǒng)離出時“爬得多陡”. 然而, 這種方式只能刻畫系統(tǒng)的局部穩(wěn)定性, 無法展現涉及全局結構的信息. 在隨機擾動下, 系統(tǒng)從不動點的整個吸引域中的離出行為是許多動力學問題的核心, 而局部穩(wěn)定性顯然不足以刻畫此類問題. 另一種穩(wěn)定性是根據吸引域半徑的長度來定義, 即考察系統(tǒng)“走得多遠”. 這種定義方式在某些情況下非常適用, 例如噪聲強度較大時的跳躍噪聲. 盡管這種穩(wěn)定性是比前一種更加全面的全局穩(wěn)定性, 在一定程度上能夠刻畫平衡點的穩(wěn)定范圍, 然而單純以吸引域大小定義穩(wěn)定性仍有一定的局限性. 例如, 定義在區(qū)域[-1,1]上的兩個分別具有勢函數U1(x)=x2與U2(x)=0.01x2的系統(tǒng), 盡管穩(wěn)定區(qū)域一致, 兩個系統(tǒng)的穩(wěn)定性明顯不同, 相應的離出難度也不一樣. 這就涉及第三種穩(wěn)定性的定義方式, 即從能量的角度考察離出需要越過的勢壘高度, 在物理化學系統(tǒng)中稱為活化能, 可以理解為“登得多高”(朱金杰等 2020). 在大偏差理論中, 擬勢扮演了這樣一個勢壘的角色, 在有勢系統(tǒng)中, 它也能與勢函數重合. 因此, 擬勢的概念從能量的角度為吸引子的全局穩(wěn)定性提供了一種定量判據, 從而噪聲誘導的不同穩(wěn)態(tài)之間的轉移現象轉化為不同擬勢勢阱的躍遷行為.

3 隨機混合系統(tǒng)研究進展

生物系統(tǒng)中的噪聲依據來源可以分為兩類, 以神經元系統(tǒng)為例, 可以分為固有噪聲 (intrinsic noise) 和外在噪聲 (extrinsic noise) (Bressloff & Newby 2013). 后者來源于環(huán)境因素, 主要由突觸輸入的持續(xù)性沖擊構成, 相關的系統(tǒng)一般模型化為基于郎之萬方程的連續(xù)Markov 過程. 而固有噪聲產生于有限個離子通道的隨機開閉, 具有明顯的離散特性, 一般模型化為跳躍Markov 過程.一般地, 若x(t)表示連續(xù)的突觸變量,n(t)∈{0,1,2,··· ,K}為離散的開通道數, 只考慮固有噪聲的情況下,x(t)在n(t)兩次跳躍之間滿足確定性方程

其中τ為突觸時間常數, 而n(t)的演化服從轉移率為Wnn′(x)/τn的跳躍Markov 過程n′ →n, 因此[x(t),n(t)]構成耦合Markov 過程, 稱為分段確定馬氏過程 (piecewise deterministic Markov process), 或隨機混合系統(tǒng) (stochastic hybrid system) (Bressloff & Newby 2014). 為了簡化, 這里考慮連續(xù)變量是一維的情況. 除了神經元模型外, 它在物理化學、分子動力學、基因轉錄等領域都有著重要應用 (Assaf et al. 2011, Newby 2014a).

對于這個耦合系統(tǒng), 引入對應的概率密度pn(x,t)

則概率密度的演化服從微分形式的Chapman-Kolmogorov (CK) 方程

注意, 在這個系統(tǒng)中引入了兩個時間尺度, 弛豫時間尺度τ和轉移時間尺度τn. 一方面, 若考慮極 限τ →0,根據式 (9) 得到vn(x)=0,則x為依賴于n的常數, 從而概率P(n,t)=P{n(t)=n}的演化控制方程轉化為

即跳躍Markov 過程. 若n(t)每次只跳躍一步, 即n →n±1, 則系統(tǒng)簡化為更為人們所熟知且應用更廣泛的生滅過程 (birth-and-death process) (Dykman et al. 1994), 即

其中邊界條件為P(-1,t)=0,P(K+1,t)=0.令Ω±(n/K)=ω±(n)/K,x=n/K, 結合方程 (13),當K →∞時得到確定性方程

另一方面, 在諸多實際的生物物理模型中, 跳躍過程的動力學遠快于連續(xù)變量的弛豫行為,即τn ?τ(Bressloff & Newby 2014). 因此, 令τ=1,τn=ε, 可以將方程 (11) 改寫成另一種更加緊湊的形式

在ε →0極限下, 式 (15) 約化為確定性方程或稱平均場方程 (mean field equation)

其中ρn(x)是滿足Anm(x)ρm(x)=0唯 一的平穩(wěn)分布. 這里假定, 對固定的x, 矩陣A(x)是不可約的.

注意到當前者K →∞或后者ε →0時, 這兩種系統(tǒng)分別收斂于對應的確定性系統(tǒng). 換句話說,1/K和ε分別刻畫了它們的噪聲強度. 在K較大或ε較小的情況下, 許多學者研究了這兩種系統(tǒng)中弱噪聲誘導的擾動系統(tǒng)與未受擾系統(tǒng)之間的大偏差現象及相關的復雜隨機動力學行為. 下面兩子節(jié)分別介紹這兩種系統(tǒng)的研究進展.

3.1 生滅過程

考慮上述生滅過程 (13), 寫成關于x的密度形式如下

此處KΩ+(x) (KΩ-(x))表示x跳躍至x+1/K(x-1/K) 的轉移率. 對于任意光滑函數f(x), 上述過程的無限小生成元為

給定任意初始分布, 上式唯一定義了一個跳躍Markov 過程. 一個非常有趣的問題是: 當K →∞時, 該Markov 過程具有怎樣的極限行為?

事實上, 當K →∞,O(1/K)階次的跳躍步長意味著上述過程逼近其連續(xù)極限, 而O(K)階次的轉移率說明系統(tǒng)關于時間快速變化, 系統(tǒng)將快速收斂至其唯一的平穩(wěn)不變分布. 與此同時, 若系統(tǒng)具有某種遍歷性特征, 則這種“快速震蕩”使得該Markov 過程的軌道呈現出大概率集中于其“平均系統(tǒng)”, 而小概率發(fā)生偏離的現象. 這一現象可由大偏差原理解釋.

令H(x,β)=Ω+(x)(eβ -1)+Ω-(x)(e-β -1),L(x,α)=(αβ-H(x,β))為H(x,β)的Legendre 變換. 記為上述過程以x0為起始點的概率分布, 則根據大偏差原理 (Freidlin & Wentzell 2012), 有

由于H(x,0)=0,可以推得L(x,α)≥0.同時L(x,α)=0當且僅當α=Ω+(x)-Ω-(x), 立即有推論

即, 當K →∞,在O(1)時間尺度, 系統(tǒng)大概率沿著確定性軌線運動. 這實際上是大數定律的結果,該確定性運動正是上述“快速震蕩”導致的系統(tǒng)平均動力學行為. 大偏差理論相比于傳統(tǒng)大數定律和中心極限定理的優(yōu)越性在于, 它使得系統(tǒng)相對于其確定性平均的有限偏差的指數小概率可以得到估計, 并可以應用于研究系統(tǒng)長期動力學的漸近行為. 如系統(tǒng)的轉移概率PK(t,x0,x)取決于V(t,x0,x)=,平穩(wěn)概率分布(x)由擬勢V(x)支配, 平均離出時間近似具有exp(KC)階次, 而常數C由擬勢確定. 這些結果與擴散情形一致, 此處不再闡述.

一個有趣的現象是, 當K →∞且t →∞時, 系統(tǒng)在不同的時間尺度呈現不同的動力學行為.如先有K →∞再有t →∞,此時xK(t)大概率趨于初值x0的ω極限集; 而若先有t →∞再有K →∞,系統(tǒng)呈現位于所有極限集上的集中分布; 而當t ～exp(KC), 如前所述, 系統(tǒng)呈現有界區(qū)域的離出行為 (Schuss 2009). 那么, 在其他時間尺度, 隨機系統(tǒng)又有怎樣的行為呢?

令x(t)為任意一條確定性平均系統(tǒng)的軌線, 可以證明標準偏差K1/2(xK(t)-x(t))在任意有限時間區(qū)間弱收斂于高斯擴散過程ξ(t), 其對應的無限小生成元為

這實際上是中心極限定理類型的結果 (Ge & Qian 2011). 取x(t)=為確定性系統(tǒng)的穩(wěn)定平衡點,由ξ離出O(1)區(qū)域的時間近似是O(1)階的可知,xK(t)相比于的O(K-1/2)偏差發(fā)生在O(1)時間尺度.

進一步, 考慮偏差ηK(t)=Kγ(xK(t)-x(t))=(xK(t)-x(t)),0＜γ ＜1/2, 有

同樣地, 取x(t)=為確定性系統(tǒng)的穩(wěn)定平衡點, 由ηK(t)發(fā)生O(1)偏差的時間近似是O(exp(K1-2γC))階次的可知,xK(t)相比于的O(K-γ)偏差發(fā)生在O(exp(K1-2γC))階時間尺度.

相似的結果可以推廣至高維. 考慮狀態(tài)空間為εZn(Z為整數集) 的連續(xù)時間跳躍Markov 過程, 其生成元為

利用與一維情形相似的推導過程, 上述結論可自然在高維進行推廣, 此處不再贅述.

3.2 隨機混合系統(tǒng)

隨機混合系統(tǒng)在神經元模型、分子動力學、基因轉錄等領域都有大量的應用 (Assaf et al.2011, Newby 2014a, Li & Liu 2019). 其連續(xù)變量和離散變量耦合的特殊結構, 既能模型化一大類實際物理系統(tǒng), 同時也帶來了高度的復雜性. 盡管擾動一般較弱, 即ε通常較小, 長期作用下離出行為仍十分普遍. 為此, 眾多研究者分析了隨機混合系統(tǒng)的大偏差現象, 設計了多種解析與數值方法, 并得到了一系列引人注目的結果. 例如, 一些研究者利用攝動法近似計算了基因調控和離子通道等模型的平均首次離出時間 (Bressloff & Lawley 2017). 基于WKB 近似和路徑積分公式,Bressloff 和Newby (2013, 2014) 推導了隨機混合系統(tǒng)的作用量泛函和最優(yōu)離出路徑滿足的輔助系統(tǒng), 計算了隨機混合神經網絡系統(tǒng)的離出路徑和離出時間, 分析了隨機混合Morris-Lecar 系統(tǒng)離出路徑的奇異性結構. 分析發(fā)現, 類似于擴散過程的情況, 隨機混合系統(tǒng)的最優(yōu)路徑仍然滿足哈密頓系統(tǒng). 然而, 區(qū)別在于, 哈密頓量由一個符號矩陣的特征值提供, 而作用量泛函為哈密頓量的Legendre 變換 (即拉格朗日量) 的積分, 因此作用量泛函與哈密頓量一般無法給出顯式形式,這為問題的求解帶來了極大的困難與挑戰(zhàn). 為此, 許多研究者提出了諸多近似方法, 如擬穩(wěn)態(tài)擴散近似 (quasi-steady-state diffusion approximation, QSS), 絕熱近似 (adiabatic approximation), 系統(tǒng)尺寸擴展 (system-size expansion) 等方法 (Bressloff & Newby 2011, Brooks & Bressloff 2015,Keener & Newby 2011), 試圖利用擴散過程近似隨機混合系統(tǒng). 這些方法對于捕捉系統(tǒng)的高斯型漲落以及揭露擬環(huán) (quasi-cycle) (Brooks & Bressloff 2015) 的出現等局部性質具有突出的優(yōu)越性.然而, 在考慮亞穩(wěn)態(tài) (metastability) 現象時, 擬勢和最優(yōu)路徑會出現不可接受的誤差 (Li & Liu 2019).

下面引入一些典型的隨機混合系統(tǒng)的理論和數值方法, 如擴散近似方法中最常用的QSS, 處理大偏差問題的WKB 近似, 并介紹近期發(fā)展的重要結果.

3.2.1 擬穩(wěn)態(tài)擴散近似

假定ε很小, 則內部變量n的 轉移率充分大, 使得n快 速跳躍的同時x的變化很小. 這個結果表明, 系統(tǒng)位于狀態(tài)n的概率快速收斂于ρn(x). 這一事實引出了QSS 方法的基本思想, 即, 概率密度pn(x,t)可以近似分解為給定x的狀態(tài)n的平穩(wěn)分布ρn(x)與x的邊際分布的乘積 (Bressloff & Newby 2014). 具體來說, 假設概率密度為

代入微分CK 方程 (15) 中, 并利用攝動法得到C(x,t)的演化滿足的FPK 方程

其中漂移系數和擴散系數分別為

Zn(x)是滿足條件Zn(x)=0, 且是方程(28)的唯一解

注意在這里只保留了FPK 方程 (26) 的漂移與擴散的最低階項. 函數C(x,t)表示近似的擴散過程的概率密度函數. 由此可以利用擴散過程 (26) 近似原先的隨機混合系統(tǒng), 從而使系統(tǒng)簡化.

實際上, QSS 擴散近似方法的優(yōu)勢主要在于二維情況下能夠捕捉到不動點附近的高斯型漲落. 具體來說, 考察平均場方程處于接近超臨界Hopf 分岔點的情況, 系統(tǒng)存在一個穩(wěn)定的焦點,雅可比矩陣具有一對復共軛特征值. 此時確定性流表現為暫態(tài)振蕩模式, 衰減率為特征值實部.而當存在固有噪聲時, 即對于隨機混合系統(tǒng), 暫態(tài)振蕩轉化為長期的大振幅振蕩, 對應的功率譜出現一個峰值, 稱為擬環(huán) (quasi-cycle) 現象. 通過QSS 方法計算出近似的It?隨機微分方程, 在不動點處線性化后經過Fourier 變換可以得到系統(tǒng)的功率譜. 近似結果得到的功率譜會出現一個峰值, 與原隨機混合系統(tǒng)的模擬結果相當吻合 (Bressloff 2010, 2017; Brooks & Bressloff 2015). 因此,QSS 方法對于擬環(huán)現象的分析具有顯著的優(yōu)勢.

以單集群隨機混合神經網絡模型 (Li & Liu 2019) 為例, 討論QSS 擴散近似方法的應用. 對于此模型, 連續(xù)變量x代表總的突觸電流, 離散的內部變量n表示激活的神經元個數.n(t)的演化服從一個轉移率為Ω±(x,n)/ε的單步跳躍Markov 過程

其中Ω+=n(t),Ω-=F(x).F是Sigmoid 激活率, 定義為

其中F0,Γ和θ分別為速率常數, 增益和閾值. 突觸電流x(t)滿足分段確定方程

則系統(tǒng)的概率密度pn(x,t)滿足的CK 方程形式為

其中的矩陣A(x)滿足

通過計算, 平穩(wěn)分布ρn(x)滿足參數為F(x)的泊松分布, 則可得到平均場方程

平均場方程關于參數θ的分岔圖如圖1 所示,θ=0.95時,存在兩個穩(wěn)定不動點x1,x2和一個不穩(wěn)定不動點x0.

圖 1 平均場方程關于 θ從單穩(wěn)態(tài)到雙穩(wěn)態(tài)的分岔圖.θ =0.95時,存在兩個穩(wěn)定不動點 x1, x2和一個不穩(wěn)定不動點x0

其中ξt是一個高斯白噪聲, 漂移系數f(x)是由式 (34) 定義的平均場方程的向量場, 擴散系數滿足D(x)=w2F(x).

若要考察系統(tǒng)的離出現象, 則根據人們所熟知的擴散過程的大偏差理論的結果, 易得擬勢VFP(x)滿足哈密頓雅可比方程 (李揚2021)

根據經典力學的結果, 這個方程可以沿著下面的參數化軌道進行積分求解

事實上, 對于此一維模型, 擬勢可以直接通過式 (36) 求解得到

當然, 后面會看到, 擴散近似后的大偏差結果相對于原混合系統(tǒng), 在遠離不動點處會出現較大的誤差.

3.2.2 WKB 近似

如前文所述, 盡管擴散近似在捕捉不動點鄰域內的高斯型漲落等局部性質時具有突出的優(yōu)越性, 然而, 當考察離出問題等全局現象時會出現不可接受的誤差. 此時, 直接對原隨機混合系統(tǒng)應用WKB 近似方法是十分必要的.

一般地, 假設CK 方程 (15) 的平穩(wěn)概率密度滿足下面的WKB 形式

其中V(x)為擬勢,Rn(x)是給定x的狀態(tài)n的條件分布.pn(x)的值主要取決于指數項,N為前因子函數, 由ε的更高階的齊次方程決定, 這里不做討論. 將方程 (39) 代入CK 方程的時齊版本并合并ε的最低階項得到

其矩陣形式為

眾所周知, 在弱噪聲極限下, 擴散過程對應的擬勢滿足哈密頓雅可比方程. 類似地, Bressloff和Newby (2014) 證明了隨機混合系統(tǒng)的擬勢滿足下面的哈密頓雅可比方程

其中,λ0是矩陣Q的主特征值, 稱為Perron 特征值. 它是實的單特征值, 并且大于所有余下的特征值的實部. Perron-Frobenius 定理保證了它對應的特征向量是嚴格正的. 因此, 最優(yōu)離出路徑滿足輔助哈密頓系統(tǒng)

沿著這些極值路徑, 擬勢的演化滿足

可以看出, 隨機混合系統(tǒng)的哈密頓量一般沒有解析表達式, 因為它是一個符號矩陣的主特征值. 這正是隨機混合系統(tǒng)復雜性的根源所在. 因此, 用隨機微分方程逼近隨機混合系統(tǒng)的QSS 等方法就顯示出它們的重要意義.

事實上, 由方程 (41) 的解決定的哈密頓量的選擇并不唯一. 例如, 它也可以定義為(x,p)=det(Q)(Newby 2014b). 盡管這個定義有一定的優(yōu)勢, 即哈密頓量總能寫出解析表達式,但它有一個缺點, 即這樣定義的哈密頓量不是動量p的凸函數, 而凸函數性質是處理更高維問題的幾何最小作用量方法 (Heymann & Vanden-Eijnden 2008)、有序逆風法 (Cameron 2012) 等數值方法所需的本質特性. 盡管如此, 可以證明H和具有相同的極值路徑, 只是時間參數化不同(Newby 2014b). 一般地, 將系統(tǒng)的哈密頓量取為Perron 特征值.

上一子節(jié)的隨機混合神經網絡模型是少有的能夠解析計算出哈密頓量的例子. 應用WKB近似方法, 通過將Rn(x)假設成Poisson 分布的形式, 根據方程 (41), 可推導出哈密頓量 (Li &Liu 2019)

最優(yōu)離出路徑和擬勢可以通過積分對應的常微分方程組 (43) 和 (44) 得到.

根據式 (45), 零值哈密頓量分別導出兩條零能量軌道, 即對應于p=0的弛豫軌道和對應于p=1/w-F(x)/x的活化軌道, 如圖2 所示. 容易看出, 每條最優(yōu)路徑都是連接不動點和鞍點的異宿軌. 如果p=0, 可以發(fā)現輔助哈密頓系統(tǒng) (43) 約化為平均場方程, 因此系統(tǒng)在弛豫軌道上運動不需要能量. 這表明動量實際上度量了最優(yōu)路徑和確定性路徑之間的距離. 因此, 類似于擴散過程, 動量p的值可以看成是反映噪聲影響的指示器. 然而, 系統(tǒng)在活化軌道上的運動需要克服穩(wěn)定不動點的吸引性, 能量當然是需要的, 大小為

圖 2 哈密頓系統(tǒng)的相圖和離出路徑.粉色和亮藍曲線分別代表從 x1到x2和x2到 x1的離出路徑.綠色和棕色虛線表示擴散近似的離出路徑

如前所述, 若考慮噪聲誘導的離出現象, 則隨機混合系統(tǒng)與其擴散近似的結果相比有顯著的偏差, 如圖3(a)所示的擬勢. 對于亞穩(wěn)態(tài)問題, 這種偏差即使很小, 也會導致QSS 方法的失效. 其原因是, 擬勢出現在MFPT 的指數項中, 因此擬勢的微小差異將導致MFPT 的較大誤差. 通過Monte Carlo 模擬 (Li & Liu 2019), 針對不同的噪聲強度, 計算出從右邊不動點穿過鞍點的MFPT,如圖3(b)所示. 利用MFPT 的對數與擬勢、噪聲強度的倒數成正比的事實, 可以很容易地得到MFPT 的對數與噪聲強度倒數的函數關系. 在圖3(b)中, 擬合直線的斜率代表了擬勢的大小, 約為0.097 8, 這與原系統(tǒng)的WKB 近似結果0.0964 一致, 而與QSS 給出的結果有很大差異.

Li & Liu (2019) 給出了QSS 方法失效的原因. 通過對比擴散近似前后的兩個哈密頓量H(x,p)=-xp+wp·F(x)/(1-wp)和(x,p)=[-x+wF(x)]p+w2F(x)p2,容易發(fā)現(x,p)恰好是H(x,p)關于動量p的二階近似. 因此, 在不動點或鞍點鄰域內,H(x,p)和(x,p), 以及它們對應的離出軌道, 吻合地很好. 在圖2 中, 動量p關于坐標x的曲線的斜率代表擬勢的二階導數. 基于此, 可以發(fā)現, 在不動點鄰域, 有效FPK 和原系統(tǒng)的結果幾乎是重合的. 這一特殊事實意味著,QSS 前后計算的擬勢關于x的二階以內是一致的. 然而, 在相空間的其他區(qū)域, 由于動量很大, 哈密頓量的二階近似是不夠精確的. 事實上, 正是遠離不動點處的較大的動量推動系統(tǒng)克服了不動點的吸引力, 使稀有的離出事件得以實現. 這就是QSS 失效的原因. 值得注意的是, FPK 方程的解所提供的哈密頓量必然是動量p的 二次函數, 而矩陣Q的Perron 特征值一般不會正好是這種形式. 鑒于此, 當考慮離出事件時, 擴散近似的失效是自然的.

圖 3 (a)擴散近似和原始系統(tǒng)的擬勢對比,(b)平均離出時間的對數和噪聲強度的倒數之間的關系

圖 4 細致平衡條件分析示意圖

3.2.3 細致平衡條件

細致平衡 (detailed balance) 原是物理學中的一個概念, 由van Kampen (1957)、Graham 和Haken (1971) 獨立地應用于FPK 方程. 粗略地說, 一個馬爾可夫過程, 如果在平穩(wěn)情形下, 每一個可能的狀態(tài)轉移都與其逆轉移平衡, 那它就處于細致平衡. 細致平衡現象與大偏差問題具有密切聯系, 若系統(tǒng)處于細致平衡, 則從穩(wěn)態(tài)漲落到某一點的最優(yōu)離出路徑與該點在無噪聲情況下弛豫到穩(wěn)態(tài)的確定性軌跡完全重合, 只需將時間反向. 因此, 離出路徑構成的拉格朗日流形的奇異性也不會出現, 問題得到了極大地簡化.

表1 中列出了幾種典型隨機過程的細致平衡條件. 眾所周知, 對于加性高斯噪聲的具有勢函數的隨機微分方程, 細致平衡自動成立. 對于其他情況, 朱位秋 (1992) 給出了擴散過程處于細致平衡時漂移和擴散系數需要滿足的條件, 并以此可得到系統(tǒng)的平穩(wěn)勢, 從而能夠求解平穩(wěn)FPK方程. 對于一般的齊次Markov 過程, 在滿足若干假設下, Gardiner (1985) 推導了此類過程微分形式的CK 方程, 并利用其系數導出了此過程滿足細致平衡的充要條件. 基于此, Dykman等 (1994) 證明了跳躍Markov 過程的細致平衡條件為WjiPi=WijPj.

表 1 幾種典型隨機過程的細致平衡條件

對于隨機混合系統(tǒng), Li 和Liu (2019) 證明了此類過程必然不滿足細致平衡. 數學上, Gardiner(1985) 推導的齊次Markov 過程處于細致平衡的充要條件在隨機混合系統(tǒng)的情況下約化為

根據第一個方程, 得到vn(x)=0.存在兩種可能: 一方面, 若vn(x)與x無關, 則可解出狀態(tài)n為某個常數, 即狀態(tài)n不存在轉移, 此時離散變量已不再是隨機的, 這時系統(tǒng)已不再是隨機混合系統(tǒng). 另一方面, 若vn(x)與x相關, 則可得到x=f(n), 意味著一個連續(xù)變量等于一個離散變量. 顯然, 這個事件概率為零. 綜上, 隨機混合系統(tǒng)必然不滿足細致平衡.

前面的解釋是由數學語言給出的, 然而怎樣才能直觀地理解這種奇特的現象呢? 如圖4所示, 假設系統(tǒng)在時刻t跳到狀態(tài)n, 然后它會在狀態(tài)n上停留一個正隨機時間τ(τ=0的概率為零), 然后它將跳到另一個狀態(tài). 在這個時間間隔內, 系統(tǒng)沿著由x˙=vn(x)定義的確定軌跡移動.假設系統(tǒng)起始于(n,x),在充分小的時間區(qū)間 Δt內,系統(tǒng)從狀態(tài)n跳躍到其他狀態(tài)的概率為CΔt, 留在狀態(tài)n的概率為1-CΔt.因此, 若vn(x)不等于零, 則轉移概率

然而, 若系統(tǒng)起始于狀態(tài)x+vn(x)Δt,在n保持不變的情況下, 系統(tǒng)在 Δt時間之后不可能到達狀態(tài)(n,x),這是因為x˙=vn(x)是確定性方程, 而相軌跡不可能與自身相交. 因此,n必須在 Δt時間區(qū)間內跳躍數次才能抵達狀態(tài)(n,x). 鑒于隨機軌道精確到達某一點的概率必然為零, 因此轉移概率

于是, 正向與反向的運動概率不可能相等, 除非vn(x)=0. 這與前面的分析一致, 即對于隨機混合系統(tǒng), 時間不可逆性總是存在的.

4 指數輕跳躍過程研究進展

作為高斯白噪聲的推廣, Lévy 噪聲是另一種近年來備受關注的隨機過程. 類似于高斯白噪聲, Lévy 噪聲可以看作是Lévy 過程的形式導數. 布朗運動是具有平穩(wěn)獨立增量性和概率1 連續(xù)樣本路徑的高斯過程, 而Lévy 過程去掉了高斯分布的要求, 并將概率1 連續(xù)退化至隨機連續(xù). 由此, Lévy 過程及Lévy 噪聲驅動的隨機微分方程的解一般是不連續(xù)的, 甚至跳躍點可能是稠密的(Duan 2015). 這一方面使得此類問題的處理變得相當困難, 但另一方面在研究地球氣候的tipping現象 (Serdukova et al. 2017, Zheng et al. 2020)、基因轉錄過程中的突變 (Cheng et al. 2019, Chen X et al. 2019) 與股票市場的“笑狀波幅 (implied volatility smile)” (Poirot & Tankov 2006) 等現象時, 人們發(fā)現Lévy 過程更加符合實際模型.

具體來說, 令Lt為Rn上的Lévy 運動, 定義跳躍過程ΔL(t)為

其中Lt-是Lt在t時刻的左極限.對于一個Borel 集S ∈B(Rn{0})和t>0,定義[0, t)時間區(qū)間內跳躍幅度落在S內的次數為

其中ω為樣本路徑.

定義一個在B(Rn{0})上的Borel 測度ν為

稱為Poisson 隨機測度, 或跳躍測度 (jump measure). 進一步, 定義補償Poisson 隨機測度為

則對任意t≥0,Lévy 運動Lt具有如下Lévy-It?分解 (Duan 2015)

其中b ∈Rn,為 具有某個協方差矩陣Q的n維布朗運動. 數字1 將大跳躍與小跳躍分開, 實際上可以用任何正數替代. (b, Q, ν)稱為Lévy 運動的生成三元組 (generating triplet).

根據Lévy-Khintchine 公式 (Duan 2015), Lévy 運動的特征函數為

其中I為示性函數. Lévy 運動的無限小生成元為

其中H為函數的Hessian 矩陣, Tr表示矩陣取跡. 實際上, 對于受Lévy 噪聲激勵的隨機動力系統(tǒng),式 (54) 中的前兩項可以分別并入漂移與擴散項中. 因此, 人們通?？紤]純跳的Lévy 運動(0,0, ν).關于Lévy 運動及具有Lévy 噪聲的隨機微分方程的更多性質, 見文獻 (Duan 2015, Applebaum 2009).

在滿足某些假定下, Freidlin 和Wentzell (2012) 證明了Lévy 噪聲和高斯白噪聲驅動的隨機微分方程滿足大偏差理論, 并隱式的給出了其作用量泛函. 然而, 這些條件很難驗證. 為此, Lipster和Puhalskii (1992) 利用隨機微分方程的系數給出了這些假定的部分充分條件. Budhiraja 等(2011) 對具有Lévy 噪聲的隨機系統(tǒng)建立了Poisson 隨機測度的泛函的變分公式, 并且推導了另一個滿足大偏差原理的充分條件. 隨后, 他們成功地將這一研究結果應用于滿足合適條件的有限維 (Budhiraja et al. 2013) 和無限維 (Budhiraja & Fischer 2012) 隨機動力系統(tǒng)的研究之中. 對于一維梯度系統(tǒng), Imkeller 等 (2009) 研究了指數輕跳躍過程 (exponentially light jump process) 擾動的首次離出時間問題.

指數輕跳躍過程是一類特殊的Lévy 過程, 其跳躍測度滿足ν(dy)=e-|y|αdy. 考慮隨機動力系統(tǒng)

4.1 大偏差原理

為了簡化, 考慮一維的情況. 當α≥1時, 隨機系統(tǒng) (57) 的解服從大偏差原理 (de Oliveira Gomes 2018, Yuan & Duan 2019)

其作用量泛函為

對任意的x,z ∈R, 定義擬勢為 (de Oliveira Gomes 2018)

一般地, 若起始點z為系統(tǒng)的不動點, 則擬勢簡記為V(x).若D為包含不動點的待離出區(qū)域, 記和首次離出時間

則對任意的δ >0和x ∈D, 有 (de Oliveira Gomes 2018)

然而, 作用量泛函 (59) 的形式并不適合最優(yōu)離出路徑的計算. Yuan 和 Duan (2019) 給出了等價的另一種形式的作用量泛函

然而, 不同于高斯噪聲的情況, 式中的拉格朗日量不能顯式解析地表達. 幸運的是, 對應的哈密頓量, 即拉格朗日量的Legendre 變換, 具有下面的解析表達式

這里,p=?L/?x˙為動量. 積分項ypχ|y|≤1實際上為零, 因為測度ν關于y是對稱的.

基于此, Li 等 (2020a) 推廣了高斯情況下的計算最優(yōu)路徑和擬勢的哈密頓形式 (Hamiltonian formalism) 應用于此類系統(tǒng). 如前所述, 作用量泛函實際上支配著軌跡實現的可能性, 其全局最小值對應于概率最大的路徑. 于是, 通過作用量泛函變分得到的輔助哈密頓系統(tǒng)的解提供了最優(yōu)離出路徑

于是, 這個哈密頓系統(tǒng)的積分解(x(t), p(t))到坐標空間的投影x(t)提供了最優(yōu)離出路徑.

從式 (64) 中注意到H(x,0)=0,則H(x1,0)=0 (x1為確定性系統(tǒng)的穩(wěn)定不動點). 由于這樣兩個事實, 即最優(yōu)路徑連接不動點且哈密頓量在連續(xù)的相連軌道上為常數, 因此在整個最優(yōu)路徑上H(x,p)=0.注意到起始于x1的最優(yōu)路徑在x →x1時t0→-∞, 根據式 (60)、式(63) 及Legendre 變換, 可以得到, 在最優(yōu)路徑上

或等價地

注意, 方程 (65) 和 (67) 構成一組計算離出路徑和擬勢的完備微分-積分方程組. 然而, 仍然需要初始條件來進行數值積分. 從理論上講, 當初始時間趨于負無窮時, 最優(yōu)路徑應該起始于不動點, 而在實際中這是不可能的. 于是, 初始點應選在擴展相空間的不動點的不穩(wěn)定流形上, 并接近不動點. 因此, 需要在不動點附近進行局部線性化. 將B記為方程 (65) 的雅可比矩陣

可以看到B有兩個特征值C ＜0,-C >0.對應于特征值-C的不動點的不穩(wěn)定特征向量為(1,-2C/D)T. 因此, 最優(yōu)路徑的初始條件可以選為

其中δ是一個小參數.

Li 等 (2020a) 給出了類似高維情況下具有指數輕跳躍噪聲的隨機動力系統(tǒng)的哈密頓形式, 以能量平衡模型和Maier-Stein 系統(tǒng)為例, 數值計算了系統(tǒng)的最優(yōu)離出路徑和擬勢. 研究發(fā)現, 在一維情況下, 擬勢的形狀與勢函數完全相同, 表明這種噪聲定性類似于高斯噪聲. 研究結果同時給出了擬勢與跳躍測度的參數之間的解析近似關系, 并分析了跳躍測度對最優(yōu)離出路徑的影響. 在高維情況下, 發(fā)現了跳躍測度誘導的拉格朗日流形的奇異性.

4.2 中度偏差原理

中度偏差原理與大偏差原理的形式是類似的, 只是指數衰減率的階數更低一些. 一般地, 中度偏差原理建立了隨機過程xt穿過給定函數φt的小鄰域的概率估計式

其中,S0T[φ]為作用量泛函,ε →0時,a(ε)→0且ε/a2(ε)→0,因此階數ε/a2(ε)低于大偏差原理的ε.

由于理論上的難度, 目前關于中度偏差原理的研究還相對較少. Hu 和Shi (2004) 對一類具有“布朗勢”的擴散過程建立了長時行為的中度偏差估計. Budhiraja 等 (2012, 2013) 對滿足充分條件的有限維和無限維Lévy 型隨機系統(tǒng)建立了中度偏差原理. 在此基礎上, de Oliveira Gomes(2018) 證明了0＜α ＜1時的指數輕跳躍過程滿足中度偏差原理, 并估計了平均離出時間的大小.

當0＜α ＜1時,可以證明方程 (57) 的解服從速率為εα的中度偏差原理, 即a(ε)=ε(1-α)/2, 其作用量泛函為

為此, 定義一個重整化偏差過程(renormalized deviation process)

可以證明, 此過程同樣服從速率為εα的中度偏差原理 (de Oliveira Gomes 2018), 作用量泛函為

類似地, 定義擬勢為

此時的擬勢在有界區(qū)域內始終是有界的. 類似地, 若起始點z為系統(tǒng)的不動點, 則擬勢簡記為V(x).若D為 包含不動點的待離出區(qū)域, 記

給定y0=x ∈D, 考察離出時間

則對任意的δ >0和x ∈D, 有

根據過程yt的定義, 這個隨機時間τε實際上是xt從區(qū)域+a(ε)D中的首次離出時間. 若考察未受擾系統(tǒng)的不動點, 即則中度偏差所關心的問題的本質實際上是系統(tǒng)從不動點的a(ε)階小鄰域中的離出現象.

5α 穩(wěn)定Lévy 噪聲研究進展

非高斯α穩(wěn)定Lévy 噪聲是最常見的一種Lévy 噪聲, 具有廣泛的實際應用, 其原因主要有兩點. 一方面, 穩(wěn)定隨機變量X的定義為:X是一個尺度化序列(Sn-bn)/an的依分布極限, 其中SnX1+X2+···+Xn,{Xi}是一列獨立同分布的隨機變量序列,an >0,bn是某個實序列, 這里并不要求{Xi}具有有限的均值或方差 (Duan 2015). 也就是說, 穩(wěn)定隨機變量相當于高斯隨機變量的推廣, 其定義可以看成是廣義的中心極限定理. 另一方面, 以旋轉對稱α穩(wěn)定Lévy 運動為例,其跳躍測度為

其中的常數為

其中Gamma 函數Γ定義為. 無限小生成元為

其中, 右端項理解為柯西主值積分 (Cauchy principal value integral). 基于Fourier 變換發(fā)現, 這個積分算子本質上是一個分數階拉普拉斯算子. 這兩點使得α穩(wěn)定Lévy 噪聲不僅具有理論上的研究價值, 對實際物理系統(tǒng)也有重要的應用.

例如, 基于真實的格陵蘭冰芯測量數據, Ditlevsen (1999) 發(fā)現氣候變化系統(tǒng)可以建模為具有布朗運動和α穩(wěn)定Lévy 過程的隨機微分方程. 隨后Zheng 等 (2020) 發(fā)展了一個概率框架來研究在溫室效應和非高斯Lévy 噪聲共同作用下的能量平衡系統(tǒng)的最大可能氣候變化. 一些研究者發(fā)現在可激神經元模型中α穩(wěn)定Lévy 噪聲能夠誘導隨機共振現象 (Patel & Kosko 2008, Cognata et al. 2010). 此外, Liu 等 (2018) 研究了具有暫時免疫和Lévy 跳躍的延遲接種SIR 流行病模型的持久性和滅絕性. Guarcello 等 (2016) 通過將電流偏置的長約瑟夫森結建模為受振蕩場驅動并受外部非高斯噪聲影響的sine-Gordon 方程, 對其動力學進行了數值研究. 理論上, 具有布朗運動和Lévy 過程的隨機微分方程可以用來模型化一大類重要的馬氏過程, 稱為Feller 過程 (B?ttcher 2010).

這些工作表明, 具有α穩(wěn)定Lévy 噪聲的隨機動力系統(tǒng)是具有深刻科學意義的現象學模型. 本節(jié)針對此類系統(tǒng), 介紹離出現象及其相關動力學行為的主要處理方法和近期研究進展.

5.1 平均首次離出時間和離出概率

考慮一維隨機動力系統(tǒng)

其中Bt為布朗運動,為對稱α穩(wěn) 定Lévy 過程, 跳躍測度滿足να(dy)=c(1,α)|y|-1-αdy,Λ和σ分別是對應的噪聲強度. 則系統(tǒng)的生成元為

引言中方程 (1) 和 (2) 分別給出了高斯白噪聲條件下平均首次離出時間和離出位置分布滿足的狄利克雷邊值問題. 根據Dynkin 公式, 針對形如式 (80) 的系統(tǒng), 平均離出時間v(x)仍滿足類似的方程 (Duan 2015)

只是其中的算子改為生成元 (81), 邊界條件由D的 邊界上的值改成D的補集Dc上的值.

另一方面, 由于Lévy 過程的跳躍特性, 離出位置分布的概念已不再適用, 因為系統(tǒng)離出時未必與邊界相交. 取而代之的是離出概率 (escape probability) 的概念, 即系統(tǒng)起始于區(qū)域D中的x,首次離出區(qū)域D時落在Dc的某個子集E內的概率, 記為PE(x), 滿足下面的方程 (Duan 2015)

根據生成元的形式 (81), 可以發(fā)現方程 (82) 和 (83) 是微分-積分方程, 并且積分算子是奇異積分, 因此求解的難度遠大于僅有高斯噪聲的情況. 由于復雜性較高, 目前處理此類問題的解析方法較少, 只有在特殊情況下能近似求解. 例如Λ=0,σ為小量, 即無高斯噪聲且Lévy 噪聲較弱,此時基于α穩(wěn)定Lévy 噪聲尾部的多項式衰減特性, 可利用多項式級數攝動近似平均離出時間和離出概率 (Qiao & Duan 2015).

然而, 多數情況下, 得到這兩個方程的理論解是較為困難的, 為此發(fā)展相應的數值方法是十分必要的. 通過將積分分割為并數值離散化, 且利用差分代替微分,Gao 等 (2014) 設計了一個有效的數值方法用來計算平均離出時間和離出概率. 基于此方法, Wu等 (2018) 計算了一個基因表達模型的平均離出時間、離出概率和隨機吸引域 (stochastic basins of attraction), 采用非高斯與高斯噪聲強度之比來識別轉移過程中的最優(yōu)選擇. Cai 等 (2017) 通過計算離出時間與離出概率來刻畫神經元系統(tǒng)中α穩(wěn)定Lévy 噪聲誘導的興奮事件. Xu 等 (2016)通過計算穩(wěn)態(tài)概率密度和平均離出時間, 研究了非高斯噪聲對相干轉換和開/關轉換的影響. 反過來, 基于對平均離出時間的觀測, Zhang 等 (2020) 設計了一個數據驅動方法, 能夠從離出時間數據中提取具有高斯布朗噪聲或非高斯Lévy 噪聲的隨機微分方程.

5.2 Onsager-Machlup 理論

為了從理論上確定最優(yōu)離出路徑, 本節(jié)引入Onsager-Machlup (OM) 理論. OM 理論是研究隨機動力系統(tǒng)離出現象的常用方法, 其基本思想是利用圍繞某條路徑的固定厚度的管道的概率來代替這條路徑實現的可能性, 可以看成是Freidlin-Wentzell 大偏差理論在有限噪聲強度下的推廣. 在高斯噪聲情況下, OM 作用量泛函比Freidlin-Wentzell 作用量泛函多了一項漂移系數的散度項, 在弱噪聲極限下OM 作用量泛函退化為Freidlin-Wentzell 作用量泛函的結果.

早在1953 年, Onsager 和Machlup (1953) 首先導出了具有線性漂移系數和常擴散系數的擴散過程的OM 泛函. Tisza 和Manning (1957) 隨后對非線性系統(tǒng)進行了擴展. 此外, Dürr 和Bach(1978) 提出了另一種推導OM 泛函的方法, 將Girsanov 變換應用于擴散過程誘導的測度. Chao和Duan (2019) 將此方法推廣到求解 (非高斯) Lévy 噪聲和 (高斯) 布朗噪聲下隨機動力系統(tǒng)的更復雜情況. Tang 等 (2014, 2017) 從另一個方面進一步推導了帶乘性高斯噪聲的過阻尼朗之萬方程和一般隨機解釋的OM 泛函. 綜合以上研究成果, 本節(jié)介紹較為一般的乘性高斯噪聲和加性非高斯Lévy 噪聲驅動下的隨機動力系統(tǒng)的OM 理論.

考慮下面的n維隨機動力系統(tǒng)

其中Bt=[B1,t, B2,t ··· , Bn,t]T是n維布朗運動,Lt=[L1,t, L2,t ··· , Ln,t]T是非高斯Lévy 過程, 其跳躍測度ν滿足. Lévy 過程的大幅值跳躍可以通過交錯 (interlacing) (Applebaum 2009) 方法來處理, 因此這里只考慮小跳躍是第4 節(jié)引入的補償Poisson 隨機測度. 向量f(x)=[f1(x), f2(x), ··· , fn(x)]T為漂移系數,ΛΛT是擴散矩陣.

對于具有乘性高斯噪聲的方程 (84), 選擇積分方法時的模糊性會導致不同的隨機解釋, 一般表示法為α-解釋 (Shi et al. 2012). 為了避免與α穩(wěn)定Lévy 噪聲的參數混淆, 采用κ解釋代替α.κ=0,κ=1/2和κ=1分別對應于It?, Stratonovich 和anti- It?的解釋. 通過修正漂移項, 隨機微分方程 (84) 可以轉化為統(tǒng)一的Stratonovich 形式

其中

記為修正的漂移系數. 使用這種統(tǒng)一的Stratonovich 形式的優(yōu)點是可以簡單地應用普通的微積分法則, 由此為路徑積分法推導后續(xù)的OM 函數 (Tang et al. 2014) 帶來了一些便利. 事實上,統(tǒng)一隨機解釋形式的選擇只影響推導過程, 而不影響OM 函數的結果.

在隨機擾動下, 亞穩(wěn)態(tài)之間可能發(fā)生躍遷. 在這個過程中, 最感興趣的是確定最大可能的轉移路徑. 由于單條路徑的概率必然為零, 需要研究隨機軌跡通過某條路徑的鄰域或管道的概率.在給定管道厚度的條件下, 軌道留在管內的概率實際上描述了這一特定路徑實現的可能性. 更準確地說, 考慮一個圍繞參考路徑φ(t),t ∈[0, T]的管道. 如果δ足夠小, 則解過程x(t)落在管道內的概率可按以下形式估算

因此, 類似于Freidlin-Wentzell 大偏差理論的結果, OM 泛函的全局最小值對應于具有最大概率的路徑, 即最優(yōu)路徑. 于是, 概率的計算轉化為OM 泛函的變分問題. 基于分析力學的經典結果, 連接點x0和點xf的最優(yōu)路徑滿足歐拉-拉格朗日方程

其邊界條件為x(0)=x0和x(T)=xf. 這里將符號φ替換為了x. 因為這是一個二階微分方程, 將它轉化為一個等價的哈密頓系統(tǒng)對于數值積分會更為方便

這里,H(x,p)是OM 函數的Legendre 變換,p=OM(x˙,x)/x˙為動量. 解(x(t),p(t))到坐標空間的投影x(t)給出了最優(yōu)路徑.

根據文獻 (Chao & Duan 2019, Tang et al. 2014), 方程 (87) 的OM 函數在不考慮加性常數的情況下由下式給出

其中Tr[A]=. 為了方便, 記yνdy. 則Legendre 變換計算出的哈密頓量為

因此, 哈密頓系統(tǒng)具有下面的形式

若高斯噪聲為加性的, 方程約化為

對式 (90) 中OM 函數的形式作兩點說明. 一方面, 如果考慮α穩(wěn)定的Lévy 運動, 則條件＜∞要求0＜α ＜1.另一方面, 當令式 (90) 中的ν=0時, 它恢復為擴散過程的OM函數. 換句話說, Lévy 噪聲對躍遷現象的影響反映在式 (90) 中OM 函數的第三項中. 特別地, 如果Lévy 運動是對稱的, 則結果與擴散過程的結果是一致的.

5.3 計算最大可能路徑的機器學習方法

注意5.2 節(jié)中的方程 (93) 是一個常微分方程的兩點邊值問題, 通?？梢杂胹hooting 方法來處理. 也就是說, 可以調整動量的初始值并對方程進行數值積分, 直到終點x(T)抵達xf. 然而, 這種方法在實踐中仍存在兩個不足之處. 第一, 通常選擇一個稱為亞穩(wěn)態(tài)的穩(wěn)定不動點作為初始點x0來考慮它的躍遷. 由于共軛動量方程 (即式 (93) 中的第二個方程) 包含-(?f)Tp, 因此在時間上向前的數值積分在數值上是不穩(wěn)定的, 甚至是不適定的. 此外, 在高維系統(tǒng)中會出現另一個問題,即很難決定應該調整初始動量的哪一分量. Li 等 (2021) 設計了一個機器學習方法來解決這兩個問題, 算法分為哈密頓系統(tǒng)的重公式化和神經網絡算法兩部分.

首先, 為了處理動量的發(fā)散問題, 記ST[x(t)]=dt并且

其中C1={x ∈C[0,T]|x(0)=x0, x(T)=xf}. 另外, 定義

其中C0={x ∈C[0,T]|x(0)=x0}. 注意這個極小化并不在終點上施加約束. 換句話說, 函數空間C0刻畫了起始于 x0無所謂終點的連續(xù)軌道的集合. 事實上,I?(λ)和I(xf)構成一組Fenchel-Legendre 變換對.

I?(λ)和I(xf)的關系事實上允許處理最小化問題 (95) 來代替 (94). 式(95) 的變分結果導出下面的哈密頓系統(tǒng)

從式 (89) 可以看出, 邊界條件約束在初始坐標和最終坐標上, 但在動量上不存在. 而在式(96) 中, 它們被轉換成一個坐標的邊界條件和另一個動量的邊界條件. 這種操作的好處是顯而易見的. 式 (96) 的計算是時間前向積分x(t),然后時間后向積分動量p(t). 這種方法成功地解決了前面提到的shooting 方法的第一個缺點, 因為兩個方向的積分都是收斂的.

這個方法的算法總結如下:

算法1:

Step 1 給定λ的一個值和一條初始軌跡xk(t)(初始k=1);

Step 2 從p(T)=λ時間后向積分方程(xk,p)得到pk;

Step 3 從x(0)=x0時間前向積分方程(x,pk)得到xk+1;

Step 4 迭代Step 2 和Step 3 直到收斂.

接下來, 如何根據x的 終點邊界條件確定λ的值可以通過一個深度神經網絡來完成, 該神經網絡的結構如圖5 所示. 根據研究目的, 輸入和輸出狀態(tài)分別固定為xf和λ,每個狀態(tài)都有n個分量.假設在輸入和輸出之間有L層隱藏層, 第l層有nl個神經元,l=1,2, ··· , L. 因此, 一個復雜的非線性函數可以用一系列簡單函數的復合來近似

圖 5具有 L層隱藏層的神經網絡的結構. xfi和λi,i=1, 2, ··· , n分別是神經網絡的輸入和輸出.表示第 l層第 j個神經元的值,其中j =1, 2, ··· , nl,l=1, 2, ··· , L

每層的函數g(l)(·)定義為

為方便起見, 記a(0)=xf和a(L+1)=λ.此外,W(l)和b(l)分別稱為權重矩陣和偏置向量.σ(l)表示一個非線性激活函數, 應用于變量的每個分量上, 常用的函數有sigmoid, tanh, ReLu 等. 接下來, 對所有隱藏層使用ReLu 函數ReLu(x)=max{0,x}.而輸出層選擇恒同函數, 以覆蓋λ的空間.

令θ為這個神經網絡參數的集合, 即θ={W(l), b(l):l=1,2, ··· , L+1}. 神經網絡的訓練是通過對參數θ進行優(yōu)化, 使其輸出狀態(tài)能夠最優(yōu)地逼近λ(xf).具體來說,假設有數據集D=, 并引入網絡預測和目標之間的L2距離作為損失函數(loss function). 然后通過搜索最優(yōu)參數對神經網絡進行訓練, 求解損失函數的非線性回歸問題

λNN表示輸入數據時網絡的輸出. 一般地, 最小化參數通過梯度下降法進行迭代

其中學習率 (learning rate)η是一個較小的數.

神經網絡的訓練過程總結如下:

算法2:

Step 1 給定數據集和初始參數集θ0;

Step 2 應用前向傳播計算網絡的輸出和代價函數;

Step 3 應用后向傳播計算代價函數對參數的梯度;

Step 4 利用方程 (100) 迭代參數;

Step 5 迭代直到收斂.

作為總結, 算法框架如下:

算法3:

Step 1 生成數據: 具體來說, 在λ的空間的一個合適的區(qū)域隨機選擇M個點. 將初始點固定在系統(tǒng)的穩(wěn)態(tài), 利用算法1, 對λ的每個點計算最大可能路徑x(t);

Step 2 訓練神經網絡: 基于算法2, Step 1 生成的數據可用于訓練神經網絡, 其中x(T)為輸入,λ為輸出;

Step 3 測試: 給定某個終點xtest,利用訓練好的神經網絡來計算對應的輸出λtest;

Step 4 計算最大可能路徑: 再次利用算法1 對λtest計算最大可能路徑, 將終點與xtest做對比.

此機器學習方法在兩個典型例子中的成功應用(Li et al. 2021a), 證實了它對于有或沒有Lévy 噪聲、加性或弱乘性高斯噪聲、不同隨機積分解釋 (It?意義、Stratonovich 意義和anti-It?意義) 的系統(tǒng)以及各種維數的系統(tǒng)的有效性. 結果表明, 該方法相比非高斯噪聲更適用于高斯噪聲, 相比乘性噪聲更適用于加性噪聲, 該算法對Stratonovich 意義是最有效的, 然后anti- It?的,再之后It?的.

這種方法還可以推廣到其他類似的問題. 例如, 除了OM 理論外, 在計算Freidlin-Wentzell 的大偏差理論下的最大可能路徑時仍然可行. 更廣泛地說, 它也可以應用于其他領域, 如最優(yōu)控制問題等.

最后, 值得注意的是, 該算法在實際應用中還存在一個挑戰(zhàn). 如果積分結果的數據xf分布相對均勻, 則可以更有效、更準確地訓練神經網絡. 然而, 如果這些點集中在幾個特定的區(qū)域中, 由于損失函數在各個區(qū)域的權重差異過大, 則效果不佳. 因此, 該算法的成功應用要求系統(tǒng)的向量場具有某種均勻性.

5.4 計算最大可能路徑的最大似然法

盡管具有Lévy 噪聲的隨機動力系統(tǒng)的OM 作用量泛函已有結果, 然而它只對某一類Lévy噪聲成立, 應用具有局限性. 另一方面, 理論結果也需得到數值方法的驗證. 因此, 發(fā)展其他的計算最大可能路徑的方法是十分必要的. 基于轉移概率密度滿足的Fokker-Planck 方程, Zheng 等(2020) 開發(fā)了一種計算最大可能路徑的最大似然法.

為了方便表示, 以一維系統(tǒng)為例

p(x(t)=u)代表隨機微分方程 (101) 的解位于x(t)=u的概率密度, p(u,t|z,s)表示給定x(s)=z的條件下到x(t)=u的轉移概率密度, 其中0 ≤s ＜t≤tf, 即

根據 (Duan 2015), 轉移概率密度p(u,t|z,s)滿足下面的Fokker-Planck 方程

右邊的積分部分實際上是非局部或分數階拉普拉斯算子, 對應非高斯α穩(wěn)定Lévy 擾動. 轉移概率密度的初始條件滿足

對任意的t ∈[0,tf],x, x0, xf ∈Rn,假設條件概率密度PA(x,t)用于刻畫在條件A下x(t)在時刻t位于點x的概率密度. 這里, 條件A指初始條件x(0)=x0和終止條件x(tf)=xf. 根據馬氏性, 它可以表示為

在條件A下, 對任一時刻t ∈[0,tf],密度函數PA(x,t)出現的峰值位置對應狀態(tài)xm(t). 這意味著, 在給定時刻t, 條件概率密度PA(x,t)的峰值位置xm(t)記為這些隨機軌道的最大可能位置, 即

將所有這些位置連接起來即得到從初始位置x(0)=x0轉移到終點位置x(tf)=xf的最大可能軌跡.

式 (105) 中的轉移概率密度p是 Fokker-Planck 方程 (103) 的解, 理論求解相當困難. 為此, 基于時間-空間離散化的Toeplitz 矩陣結構, Gao 等 (2016) 提出了一種快速精確的數值算法用于計算具有吸收邊界或無窮區(qū)域的非局部Fokker-Planck 方程, 并證明了其收斂性. 目前, 這個方法已經廣泛應用于離出現象 (Chen X 2019)、隨機分岔 (Wang 2020) 和數據驅動方法 (Zhang 2017) 等相關問題中. 另外, 機器學習方法也可以應用于求解Fokker-Planck 方程, 例如Xu 等 (2020) 利用深度神經網絡計算高斯噪聲下的概率密度, Chen 等 (2020) 基于樣本路徑數據求解同時具有高斯布朗噪聲和非高斯Lévy 噪聲的隨機系統(tǒng)Fokker-Planck 方程.

5.5 非局部Kramers-Moyal 公式

對于實際工程結構和自然科學等問題的研究, 人們通常利用牛頓第二定律或拉格朗日力學等理論建立系統(tǒng)的動力學方程. 通過求解此控制方程, 希望得到系統(tǒng)的響應, 或得到系統(tǒng)的一些動力學性質, 例如計算不同穩(wěn)態(tài)的轉移路徑. 然而, 由于人們對許多復雜現象的內部機制缺乏完備的理解, 如生物醫(yī)藥、腦科學、大氣科學等領域, 很難直接建立其顯式的運動方程. 幸運的是,隨著科研工具和仿真能力的進步, 人們能夠獲得越來越多的復雜系統(tǒng)數據集. 因此, 從噪聲數據中提取控制規(guī)律或動力學特性在各個科學領域都起著至關重要的作用.

針對具有高斯布朗噪聲的隨機動力系統(tǒng), Dai 等 (2020) 提出了一個數據驅動框架, 用以從樣本數據中提取最大可能轉移路徑. 此方法的基本思想是, 先利用Kramers-Moyal 公式從樣本路徑數據中提取隨機微分方程的漂移和擴散系數, 然后利用前一子節(jié)的最大似然法計算最大可能路徑. 然而, 這種方法無法直接推廣到具有非高斯Lévy 噪聲的情況. 其原因在于,α穩(wěn)定分布的尾部多項式衰減, 具有重尾 (heavy-tailed) 特性, 其二階矩及以上必然發(fā)散, 一階矩是否發(fā)散依賴于α的取值. 因此, Kramers-Moyal 公式不再適用Lévy 噪聲的情況. 另一個難點在于此時的Fokker-Planck 方程是微分-積分方程, 這也使得問題變得更加復雜. 幸運地是, Li 和Duan (2021) 推導的非局部 (nonlocal) Kramers-Moyal 公式有效解決了這個問題, 下面對此公式做簡要介紹.

考慮一個n維 隨機動力系統(tǒng)

其中b(x)=[b1(x), b2(x), ··· , bn(x)]T是 Rn中的漂移系數,Bt=[B1,t, B2,t, ··· , Bn,t]T是n維布朗運動,Λ(x)是一個n×n的矩陣,Lt=[L1,t, L2,t, ··· , Ln,t]T是具有正常數噪聲強度σ的對稱Lévy運動. 假定初始條件為x(0)=z,a(x)=ΛΛT為擴散矩陣,Lt的跳躍測度滿足ν(dy)=W(y)dy,y ∈Rn{0}.由于Lévy 運動的對稱性,W(-y)=W(y).

基于系統(tǒng) (107) 的Fokker-Planck 方程, 可以推導出下面的定理和推論.

定理1對每個ε>0,概率密度函數p(x,t|z,0)和跳躍測度、漂移、擴散具有如下的關系:

(1) 對每個滿足|x-z|>ε的x和z

對x和z一致成立;

(2) 對i=1,2, ··· , n

(3) 對i,j=1,2, ··· , n

這個定理的思想可以按照下述方式理解. 在不含Lévy 噪聲的情況下, It?隨機微分方程的樣本軌道是概率1 連續(xù)的. 也就是說, 漂移系數和擴散系數貢獻的是連續(xù)部分, 而Lévy 噪聲貢獻的是跳躍部分. 在很短的一段時間內, 大的跳躍幾乎都來自于Lévy 噪聲. 依據這個思想, 這個定理將相空間分割為球內和球外兩部分. 球外部分 (即定理的第一條) 用來表達Lévy 跳躍項, 球內部分 (即定理的第二條和第三條) 用來表達漂移和擴散項. 由于將積分區(qū)域限制在球內, 從而有效避開了重尾特性導致積分發(fā)散的問題, 成功地將漂移、擴散與Lévy 跳躍項利用轉移概率密度來表達. 然而, 這個定理對于數值算法的設計并不友好. 為此, 可以將其改寫成下面的推論形式, 從而將漂移、擴散與Lévy 跳躍項利用樣本路徑數據來表達:

推論1對每個ε>0, 隨機微分方程 (107) 的解和跳躍測度、漂移、擴散具有如下的關系:

(1) 對每個m>1

這個推論第二條和第三條與Kramers-Moyal 公式具有類似的形式, 因此這個推論可以稱為非局部Kramers-Moyal 公式. 對于推論的第一條, 考慮多個不同的區(qū)間, 可以近似得到Lévy 跳躍測度和噪聲強度. 關于第二條和第三條, 利用基函數的線性組合近似, 并結合最小二乘法, 可以估計出漂移與擴散項的近似表達式. 具體的數據驅動算法的設計參考文獻 (Li & Duan 2021), 這里不再詳細介紹. 因此, 給定樣本路徑數據, 利用非局部Kramers-Moyal 公式, 可以識別出具有高斯布朗噪聲和非高斯Lévy 噪聲的隨機動力系統(tǒng), 結合上一節(jié)介紹的最大似然法, 可以計算出系統(tǒng)亞穩(wěn)態(tài)間轉移的最大可能路徑, 從而從樣本路徑數據中提取出了最大可能路徑.

最后, 可以發(fā)現, 將這種方法擴展到具有乘性Lévy 噪聲的隨機動力系統(tǒng)的數據集是一個挑戰(zhàn). 在這種情況下, 乘性Lévy 噪聲破壞了定理1 中第一個斷言的“空間齊次性”, 導致函數W同時依賴于x和z(而不僅僅依賴于空間平移x-z), 從而給數據驅動算法的設計帶來了困難 (Li &Duan 2021).

6 展 望

嚴格意義上, 在自然科學、社會科學和工程科學等領域中, 非高斯模型比高斯過程更具有普遍性, 但同時非線性非高斯噪聲系統(tǒng)的隨機動力學行為亦更為復雜. 隨著人們對自然現象的描述越發(fā)精細, 新的問題層出不窮, 對大偏差理論在非高斯隨機動力系統(tǒng)的應用提出了更多的挑戰(zhàn).以發(fā)展的眼光來看, 從以下幾個方面提出一些潛在性問題:

(1) 如前所述, 用于計算隨機混合系統(tǒng)的離出問題的哈密頓量是一個符號矩陣的主特征值.然而, 對于高維的問題, 目前只有個別特殊系統(tǒng)能夠解析求解. 對于一般的隨機混合系統(tǒng), 如何解析近似或數值計算其哈密頓量仍是一個巨大的挑戰(zhàn);

(2) 近年來隨著計算機技術的蓬勃發(fā)展, 機器學習方法以其強大的功能被廣泛應用于解決動力系統(tǒng)中的難題, 如利用機器學習方法求解哈密頓-雅可比方程、Fokker-Planck 方程和變分問題等. 由于計算大偏差問題中的最大可能路徑和擬勢的難度較大, 在大量觀測和模擬數據的基礎上, 利用機器學習方法來求解大偏差問題具有廣闊的應用前景.

致 謝國家自然科學基金資助項目 (11772149).