趙梓成，龍 潛，董小波，孟潤宇，鐘詩言，諶俊毅，向梓琨

(1. 中國科學(xué)院云南天文臺，云南 昆明 650216；2. 中國科學(xué)院大學(xué)，北京 100049；3. 中國科學(xué)技術(shù)大學(xué)物理學(xué)院，安徽 合肥 230026)

引力N體模擬是研究星系、暗物質(zhì)暈以及宇宙大尺度結(jié)構(gòu)等天文系統(tǒng)形成和演化的主要方法之一。由于引力是長程力，對于星系以及更大尺度或更多粒子數(shù)的引力系統(tǒng)，可以把成員恒星(或暗物質(zhì)粒子)的運(yùn)動看作是總體引力勢下的運(yùn)動，而不需要考慮恒星之間的相互作用，即無碰撞引力系統(tǒng)，也是一種平均場近似的思想。對這類系統(tǒng)隨時間演化的N體數(shù)值模擬，可以不考慮N體粒子兩兩之間的引力作用，僅計算共同的引力場，以及引力場對具體粒子的作用，即無碰撞引力N體數(shù)值模擬[1]。在近來的應(yīng)用中(Gadget-3, Gadget-4)，PM-Tree方法是無碰撞引力N體模擬的主要方法之一，而PM-Tree方法中求解勢能的步驟是主要耗時的部分[2]。所以求解泊松方程的速度是影響無碰撞引力N體數(shù)值模擬耗時的主要因素之一。

在具體的引力N體模擬中，N體粒子的數(shù)目不一定對應(yīng)，往往不可能達(dá)到恒星或暗物質(zhì)粒子的真實數(shù)量，數(shù)目N的設(shè)定，取決于模擬精度的需求以及計算能力。傳統(tǒng)無碰撞引力N體數(shù)值模擬方法的計算速度適用于粒子數(shù)規(guī)模在106～1013的宇宙系統(tǒng)[1]，模擬更大規(guī)模的宇宙系統(tǒng)的分布演化則需要計算速度更快的方法。深度神經(jīng)網(wǎng)絡(luò)在諸多解方程問題中都證明是有效的，并且還有速度快的優(yōu)點，所以我們提出在傳統(tǒng)無碰撞引力N體數(shù)值模擬方法PM-Tree中，使用深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法計算主要耗時的求解勢能部分。以粒子數(shù)規(guī)模為1013的系統(tǒng)為例，本文用深度神經(jīng)網(wǎng)絡(luò)的方法，單次勢能求解的速度理論上能提高40倍左右，循環(huán)疊加后對整體模擬的加速較為可觀，從各方面說明深度神經(jīng)網(wǎng)絡(luò)加速求解泊松方程的有效性，為后續(xù)加速無碰撞引力N體數(shù)值模擬的研究打下基礎(chǔ)。因為深度神經(jīng)網(wǎng)絡(luò)具有可擴(kuò)展的特點，所以在其他需要求解泊松方程的物理問題中，深度神經(jīng)網(wǎng)絡(luò)模型也能適用。

引力N體數(shù)值模擬方法主要有：(1)PP(Particle-to-Particle)方法，直接計算粒子間的牛頓引力，計算時間復(fù)雜度為O(N2)[3]。(2)當(dāng)系統(tǒng)中粒子數(shù)量較多時，可以使用自1980年延用至今的Tree代碼加快模擬速度[4-5]。Tree代碼使用遞歸方式計算子樹，省去了重復(fù)計算，將引力N體數(shù)值模擬的計算時間復(fù)雜度降低到O(NlogN)，例如GASOLINE[6]，iVINE[7]等。(3)PM(Particle Mesh)方法：首先使用CIC(Cloud-In-Cell)方法將粒子的位置信息轉(zhuǎn)換為均勻網(wǎng)格的密度信息，然后通過快速傅里葉變換求解基于這些網(wǎng)格的泊松方程。該方法的計算時間復(fù)雜度為O(NlogN)[8]，N代表網(wǎng)格單元數(shù)。(4)為了提高PM方法中網(wǎng)格內(nèi)模擬粒子的精確率，在網(wǎng)格內(nèi)使用Tree代碼計算粒子間的相互作用，稱為PM-Tree方法[9]。由于Tree代碼每次在PM方法的一個小步驟內(nèi)并行計算，而且一次Tree代碼計算時間小于或等于PM方法中一個小步驟計算時間的二分之一[10]，所以PM-Tree方法的計算時間復(fù)雜度和PM方法相同，為O(NlogN)。至今，PM-Tree方法為無碰撞引力N體數(shù)值模擬的主要方法之一，例如Gadget-2[11]。(5)近年來的N體程序(如Gadget-3, Gadget-4)[2]還實現(xiàn)了快速多極子方法(Fast Multipole Method, FMM)[12]。FMM可以將計算時間復(fù)雜度降到O(N)，在相同精度下性能優(yōu)于Tree方法，但算法更為復(fù)雜，特別是在動態(tài)步長、特殊邊界條件、高效并行化等方面實現(xiàn)較為困難。

以上傳統(tǒng)引力N體數(shù)值模擬方法的計算速度較慢，在超大規(guī)模粒子數(shù)量問題中難以適用，所以需要速度更快的方法。目前深度神經(jīng)網(wǎng)絡(luò)在很多領(lǐng)域的成果顯著[13-14]，并且神經(jīng)網(wǎng)絡(luò)也常用于解方程[15-16]。文[17]證明神經(jīng)網(wǎng)絡(luò)能解常微分方程和偏微分方程。引力N體數(shù)值模擬問題中的高維數(shù)據(jù)所需計算量巨大，而卷積神經(jīng)網(wǎng)絡(luò)的核心能力是從高維輸入中提取多尺度特征[18-20]，有效應(yīng)對數(shù)值解偏微分方程中的維度詛咒問題[21]。所以我們采用深度卷積神經(jīng)網(wǎng)絡(luò)的方法，在無碰撞引力N體數(shù)值模擬方法PM-Tree中代替快速傅里葉變換法，更快求解勢能。

1 相關(guān)工作

近幾年出現(xiàn)了一些使用深度神經(jīng)網(wǎng)絡(luò)求解偏微分方程的方法。由無碰撞玻爾茲曼方程的速度矩定義的Jeans方程，與流體力學(xué)中(無粘的)納維-斯托克斯(Navier-Stokes)方程具有相同的形式，所以在這類流體領(lǐng)域解Navier-Stokes方程的相關(guān)研究結(jié)論也適用于引力N體數(shù)值模擬[22]。穩(wěn)定流近似是一個Navier-Stokes流體問題,是一種難以計算且非常耗時的迭代過程，文[23]使用卷積神經(jīng)網(wǎng)絡(luò)實現(xiàn)了快速準(zhǔn)確的擬合。文[24]在隨機(jī)多孔介質(zhì)流動問題中使用卷積神經(jīng)網(wǎng)絡(luò)，僅用很少的數(shù)據(jù)就得到了滿意的效果，優(yōu)于蒙特卡洛法。文[25]證明在流場預(yù)測問題上使用卷積神經(jīng)網(wǎng)絡(luò)方法的有效性，而且速度比雷諾斯平均納維-斯托克斯(Reynolds Averaged Navier-Stokes)法快。文[26]證明了神經(jīng)網(wǎng)絡(luò)可以簡潔地求解具有不確定性的偏微分方程問題。文[27]在不適定(ill-posed)逆問題上使用深度神經(jīng)網(wǎng)絡(luò)迭代求解，速度提高十分顯著。

不過這類用深度神經(jīng)網(wǎng)絡(luò)求解偏微分方程的方法是將求解算子參數(shù)化，并不直接，而且對于網(wǎng)格化尺寸不具有魯棒性，但對于引力N體粒子數(shù)值模擬來說，要處理的問題是不停迭代的，需要更簡單高效的求解方法。最近出現(xiàn)了一些針對特定問題直接求解偏微分方程的方法。文[28]在電磁領(lǐng)域求解二維圖像形式的泊松方程，得到了很好的效果。文[29]提出了深度瑞茲(Deep Ritz)方法，在使用瑞茲方法求解偏微分方程過程中使用深度神經(jīng)網(wǎng)絡(luò)來尋找Ritz法的解，證明深度神經(jīng)網(wǎng)絡(luò)方法優(yōu)于有限差分法。文[30]在無限維空間中證明使用深度神經(jīng)網(wǎng)絡(luò)參數(shù)化求解橢圓偏微分方程的有效性和收斂性。文[31]提出了一個艾克網(wǎng)絡(luò)(Eiko Net)，證明深度神經(jīng)網(wǎng)絡(luò)方法能夠穩(wěn)定求解艾克諾(Eikonal)方程，并能應(yīng)用到多種領(lǐng)域，具有非常顯著的計算速度優(yōu)勢。

2 本文工作及貢獻(xiàn)

我們結(jié)合目前深度學(xué)習(xí)速度快、有效解決高維度問題、能求解偏微分方程等特點，提出使用深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法，求解無碰撞引力N體數(shù)值模擬方法PM-Tree中的勢能。這種方法理論上可以更快地完成大規(guī)模的無碰撞引力N體數(shù)值模擬，并且可擴(kuò)展性好，易于實現(xiàn)。我們搭建的深度神經(jīng)網(wǎng)絡(luò)模型，以及對應(yīng)的訓(xùn)練和測試代碼都在Github中開源，項目地址：https://github.com/hi76/Poisson-slover-based-on-Encoder-Decoder-Neural-Network。

3 方 法

N體數(shù)值模擬最早的直接求解法耗時巨大。平均場近似思想將難解的N體數(shù)值模擬轉(zhuǎn)化為一個近似的統(tǒng)計問題，更易求解，但速度依然不快。后來陸續(xù)出現(xiàn)了一些優(yōu)化的方法如PM方法、PM-Tree方法、FMM方法等，其中PM-Tree方法相比其他大多數(shù)N體數(shù)值模擬方法有速度快、精度高的優(yōu)點，所以是目前主流的方法之一。在PM-Tree方法中主要耗時的步驟是求解勢能。我們用深度神經(jīng)網(wǎng)絡(luò)模型代替快速傅里葉變換法，加速PM-Tree方法中勢能的求解，以驗證深度學(xué)習(xí)方法在加速無碰撞引力N體數(shù)值模擬中的可行性。

3.1 N體數(shù)值模擬

物理學(xué)和天文學(xué)中，在一些物理作用力的影響下，對一個由粒子組成的系統(tǒng)的動力學(xué)模擬稱為N體模擬。引力N體模擬主要是對N粒子引力相互作用方程求數(shù)值解。引力N體模擬在天體物理學(xué)中是一個常用的工具，從少體系統(tǒng)到宇宙大尺度結(jié)構(gòu)，比如地球-月球-太陽系統(tǒng)到銀河系中都有應(yīng)用。N粒子引力相互作用方程為[32]

(1)

其中，G=6.673 00 × 10-11m3·kg-1·s-2；Φext為外部勢能。

N粒子引力相互作用系統(tǒng)定義了一個6N+ 1維的相空間，相空間中每個時刻每個粒子由N個位置和速度矢量確定。這個方程的解是相空間中的一個運(yùn)動軌跡。直接求解這個方程十分耗時，所以需要簡化的數(shù)值方法求解。

3.2 平均場法

平均場法是出現(xiàn)較早的一種N體數(shù)值模擬方法。用統(tǒng)計方法描述引力N體問題時，系統(tǒng)分布f(x1,v1, …,xn,vn,t)描述了粒子的全部信息。當(dāng)不考慮粒子間的關(guān)聯(lián)，即認(rèn)為f(2)(x1,v1,x2,v2,t)=f(1)(x1,v1,t)f(1)(x2,v2,t)，平均場近似，即無碰撞N體系統(tǒng)。平均場近似中，粒子受到的作用力取決于分布函數(shù)，不必計算粒子兩兩之間的引力相互作用，這時我們可以用單粒子分布函數(shù)f(x,v,t)構(gòu)造系統(tǒng)的平均場，即在d3xd3v的六維相空間中發(fā)現(xiàn)一個粒子的概率為f(x,v,t)d3xd3v。這個分布函數(shù)的演化可以由無碰撞玻爾茲曼方程描述[1]：

(2)

總潛在勢能為外部勢能和內(nèi)部勢能之和：

(3)

(4)

3.3 PM-Tree

PM-Tree方法也是一種N體數(shù)值模擬的方法。它在PM方法的網(wǎng)格內(nèi)加入Tree代碼，以實現(xiàn)更精細(xì)的模擬。因速度較快且精度較高，在最近的N體數(shù)值模擬應(yīng)用中，常以PM-Tree方法為主進(jìn)行模擬[2]。PM-Tree方法在每個時間步長(取決于模擬的時間區(qū)間以及模擬的精度要求)內(nèi)模擬的基本步驟可以概括為[10]

(1)首先確認(rèn)Tree代碼應(yīng)用的網(wǎng)格以及當(dāng)中的粒子，除此以外的粒子只由PM方法模擬計算；

(2)運(yùn)行PM方法的一半步驟；

(3)計算每個Tree代碼應(yīng)用的網(wǎng)格對應(yīng)的外部總勢能——其余網(wǎng)格對它的勢能和；

(4)對這些應(yīng)用的網(wǎng)格使用Tree代碼同時計算，并將計算結(jié)果加入正在運(yùn)行的PM方法步驟；

(5)計算PM方法中的勢能以及對應(yīng)每個粒子的加速度；

(6)將每個Tree整合到PM方法的最后步驟；

(7)更新粒子的速度和位置信息，從步驟(1)開始循環(huán)計算。

因為在PM-Tree方法步驟中，Tree代碼是和PM方法并行計算的，并且每次Tree代碼計算所需時間不大于PM方法每次計算時間的一半，所以PM-Tree方法耗費的時間主要為PM步驟耗費的時間。在PM方法中，求解勢能是主要耗時的計算步驟，而勢能需要求解泊松方程得到。所以加速PM-Tree方法的關(guān)鍵在于加速求解泊松方程。PM-Tree方法的應(yīng)用中，常用快速傅里葉變換法求解泊松方程[33]。本文用深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法求解PM-Tree方法中的所有泊松方程，以加速PM-Tree方法，即加速N體數(shù)值模擬。

3.4 泊松方程數(shù)值求解

傳統(tǒng)求解泊松方程的方法有有限元法、有限差分法、有限體積法、快速傅里葉變換法等。利用計算機(jī)數(shù)值求解泊松方程，首先將泊松方程轉(zhuǎn)換為代數(shù)方程，即差分方程。二維泊松方程的差分方程表示為

(5)

分子是對φ做卷積，DT為離散化參數(shù)。(5)式在計算機(jī)中的操作可以簡化表示為生成的勢能φ和一個二階差分算子?的卷積：

(6)

其中，?為卷積操作符，卷積操作符的左邊項是二階差分算子?。有限差分法求解泊松方程是將(5)式全部展開為一個有解的線性方程組，再用各類迭代法求解。有限差分法求解泊松方程的計算時間復(fù)雜度為O(N2)[34]。

快速傅里葉變換法是利用傅里葉變換將實數(shù)域中的卷積轉(zhuǎn)換為復(fù)數(shù)域的乘法，將(6)式轉(zhuǎn)換為一個復(fù)數(shù)域中易解的方程，求解后再做傅里葉逆變換，得到泊松方程的數(shù)值解[35]?？焖俑道锶~變換法的計算時間復(fù)雜度為O(NlogN)。因速度較快，在近來的引力N體數(shù)值模擬應(yīng)用中，涉及求解泊松方程的步驟，主要用快速傅里葉變換法進(jìn)行求解[2]。

三維泊松方程的二階差分算子?為

(7)

快速傅里葉變換求解三維泊松方程的原理和過程與二維相同，只是二階差分算子?不同，計算量更大。我們的目標(biāo)是訓(xùn)練一個深度神經(jīng)網(wǎng)絡(luò)，以O(shè)(N)的計算時間復(fù)雜度，實現(xiàn)求解(6)式的反卷積過程，驗證深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法求解泊松方程，實現(xiàn)加速無碰撞引力N體數(shù)值模擬的可行性。因為我們關(guān)注的問題是深度神經(jīng)網(wǎng)絡(luò)是否可以實現(xiàn)求解(6)式的反卷積過程，所以選用消耗計算資源更少的二維泊松方程進(jìn)行驗證。

目前N體數(shù)值模擬的粒子數(shù)規(guī)模主要在106～1013之間。隨著計算機(jī)計算能力的提升，能夠模擬的粒子數(shù)規(guī)模以每年兩倍左右的速度增長，符合摩爾定律。在無碰撞引力N體數(shù)值模擬中，影響模擬速度的主要因素為解泊松方程的速度。傳統(tǒng)方法中常用的快速傅里葉變換求解的計算時間復(fù)雜度為O(NlogN)[1]，而深度神經(jīng)網(wǎng)絡(luò)模型的計算時間復(fù)雜度為O(N)，以粒子數(shù)規(guī)模1013為例，每一次泊松方程求解的速度理論上能提高40倍左右。在PM-Tree方法中，求解泊松方程的步驟需要循環(huán)執(zhí)行很多次，所以循環(huán)疊加以后，計算時間復(fù)雜度O(N)的方法對無碰撞引力N體數(shù)值模擬的加速較為可觀。而且隨著將來模擬的粒子數(shù)規(guī)模不斷增大，計算時間復(fù)雜度O(N)的方法能節(jié)省的時間也更加明顯。

3.5 深度神經(jīng)網(wǎng)絡(luò)模型

圖1 神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)。數(shù)字是卷積層通道數(shù)，左側(cè)向下代表特征提取部分，右側(cè)向上代表上采樣部分。左側(cè)綠色為卷積層，紅色為激活函數(shù)；右側(cè)藍(lán)色為反卷積層，黃色為激活函數(shù)

我們的訓(xùn)練目標(biāo)是最小化預(yù)測值和勢能真值的平均平方誤差

(8)

使用最大差值絕對值來評估誤差

(9)

卷積神經(jīng)網(wǎng)絡(luò)模型的計算時間復(fù)雜度為[38]

(10)

4 訓(xùn)練模型

我們的目標(biāo)是用深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法，求解PM-Tree方法中的勢能，以驗證深度學(xué)習(xí)加速無碰撞引力N體數(shù)值模擬的可行性。其中關(guān)鍵任務(wù)是驗證深度神經(jīng)網(wǎng)絡(luò)是否可以實現(xiàn)求解(6)式中的反卷積過程，所以需要訓(xùn)練得到一個可以作為泊松方程求解器的深度神經(jīng)網(wǎng)絡(luò)模型。本文使用隨機(jī)生成的矩陣模擬網(wǎng)格化后需要預(yù)測的勢能φ，即泊松方程的數(shù)值解，對生成圖像做二階差分得到模擬的輸入分布ρ；建立一個深度神經(jīng)網(wǎng)絡(luò)模型，以多個卷積層提取特征，多個反卷積層表達(dá)特征輸出預(yù)測勢能；生成100萬幅數(shù)據(jù)圖，用平均平方誤差對這個深度神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練；最后對訓(xùn)練好的模型做精度測試，用于簡單等時勢問題的評估，并與快速傅里葉變換法做精度對比，與快速傅里葉變換法和有限差分法做時間對比，用多種網(wǎng)格化方式進(jìn)行測試。

4.1 數(shù)據(jù)集

我們生成100萬幅圖作為訓(xùn)練集，1萬幅圖作為測試集。這里我們假設(shè)將求解區(qū)域采樣為N×N的網(wǎng)格，則使用的數(shù)據(jù)中模擬勢能φ真值的構(gòu)成為隨機(jī)生成的N×N矩陣，然后將周圍一圈值設(shè)置為0作為狄利克雷(Dirichlet)邊界條件。采用Dirichlet邊界條件的原因是在引力N體數(shù)值模擬中，通常我們關(guān)注的物質(zhì)團(tuán)周圍很大區(qū)域是真空，可以近似為0，并且為了避免網(wǎng)絡(luò)訓(xùn)練過于受到隨機(jī)函數(shù)本身影響(隨機(jī)函數(shù)得到的模式是不規(guī)律跳變值，相同的或平滑的情況很少)，我們?nèi)藶榧尤胍恍┨厥饽Ｊ阶鳛閿?shù)據(jù)增強(qiáng)，對中間部分區(qū)域做高斯平滑，將中間部分區(qū)域數(shù)值設(shè)置為0，將四角部分區(qū)域數(shù)值設(shè)置為0.99等。這些圖像是模擬的泊松方程的解，也就是勢能分布φ，是深度神經(jīng)網(wǎng)絡(luò)模型的預(yù)測目標(biāo)，如圖2～圖5。

圖2 全隨機(jī)Fig.2 Rand

圖3 中間平滑F(xiàn)ig.3 Middle smooth

圖4 中間置0Fig.4 Middle setting 0

圖5 四角置0.99Fig.5 Quadrant setting 0.99

需要模擬求解的分布ρ圖，可以由以上模擬的勢能φ通過(6)式得到，如圖6～圖9，這些分布ρ圖就是深度神經(jīng)網(wǎng)絡(luò)模型的輸入。

每次訓(xùn)練給深度神經(jīng)網(wǎng)絡(luò)模型輸入ρ，用模型輸出與φ的平均平方誤差進(jìn)行訓(xùn)練。因為數(shù)據(jù)是隨機(jī)生成的，所以有不限數(shù)量的訓(xùn)練數(shù)據(jù)的優(yōu)點。每次生成100萬對數(shù)據(jù)用于訓(xùn)練。

圖6 對圖2的二階差分

圖7 對圖3的二階差分

圖8 對圖4的二階差分

圖9 對圖5的二階差分

4.2 訓(xùn)練參數(shù)

我們在Pytorch中實現(xiàn)深度神經(jīng)網(wǎng)絡(luò)模型。每次訓(xùn)練使用100萬對數(shù)據(jù)，批大小設(shè)置為100，優(yōu)化器選用Adam[39]，學(xué)習(xí)率從0.01視情況調(diào)整到0.000 001。每組數(shù)據(jù)訓(xùn)練100次，訓(xùn)練使用的圖形處理器是Tesla P100。

5 訓(xùn)練成果及測試

5.1 訓(xùn)練成果

輸入模擬分布函數(shù)ρ圖，通過神經(jīng)網(wǎng)絡(luò)模型得到預(yù)測圖，與模擬勢能φ真值圖比較并作差，訓(xùn)練成果如圖10～圖13。

×10-4

圖10 模擬ρFig.10 Simulation ρ

圖11 模擬φFig.11 Simulation φ

圖12 預(yù)測值Fig.12 Prediction

圖13 預(yù)測值與模擬φ的相對誤差

由圖10～圖13可以看出，深度神經(jīng)網(wǎng)絡(luò)模型的預(yù)測圖和模擬勢能φ真值圖肉眼幾乎看不出區(qū)別，圖像特征基本吻合，它們的相對誤差最大為2 × 10-4，精度達(dá)到泊松方程的求解目標(biāo)。

5.2 測試結(jié)果

我們使用1萬對新生成數(shù)據(jù)用于測試，結(jié)果如圖14，統(tǒng)計得到不同最大誤差精度的準(zhǔn)確率如表1。

表1 不同精度下的準(zhǔn)確率Table 1 Accuracy at different precision

圖14 10 000組數(shù)據(jù)測試結(jié)果Fig.14 Test results for 10 000 sets of data

5.3 不同網(wǎng)格化數(shù)據(jù)測試

為了驗證模型的普適性，我們用16 × 16數(shù)據(jù)訓(xùn)練的神經(jīng)網(wǎng)絡(luò)模型用于漸增網(wǎng)格化尺寸(16 × 16-1 024 × 1 024)的圖像進(jìn)行測試，如圖15～圖16。

圖15 平均最大誤差趨勢圖Fig.15 Maximum error trend graph

圖16 平均平方誤差趨勢圖Fig.16 Mean squared error trend graph

由圖15～圖16可以看出，最大誤差趨勢漸平，不隨網(wǎng)格化尺寸增大而線性增大。從平均平方誤差圖也能看出，隨著網(wǎng)格化尺寸增大，平均以后的平方誤差反而減小，當(dāng)網(wǎng)格化尺寸足夠大以后，誤差基本不變。其中兩組模擬數(shù)據(jù)(256 × 256, 1 024 × 1 024)如圖17～圖24。

圖17 模擬256 × 256 ρ

圖18 模擬256 × 256 φ

圖19 256 × 256預(yù)測值

圖20 預(yù)測值與模擬φ的相對誤差

圖21 模擬1 024 × 1 024 ρ

圖22 模擬1 024 × 1 024 φ

圖23 1 024 × 1 024預(yù)測值

圖24 預(yù)測值與模擬φ的相對誤差

由圖17～圖24可以看出，用圖10這類16 × 16尺寸數(shù)據(jù)訓(xùn)練的網(wǎng)絡(luò)，在更大尺度上依然具有普適性。經(jīng)過各種尺度圖像的對比分析，最大誤差變大的原因是更大尺寸的數(shù)據(jù)包含的局部特殊模式更多。這些特殊模式是在測試集數(shù)據(jù)測試時出現(xiàn)的那些誤差較大的點。但總體可以說明，小尺寸數(shù)據(jù)訓(xùn)練的模型可以用于大尺寸數(shù)據(jù)。這意味著模型訓(xùn)練的時間成本可以非常低。

5.4 等時勢問題測試

我們將訓(xùn)練好的模型用于一個實際的物理問題——等時勢問題。按生成等時勢問題的φ和ρ圖像。將ρ輸入訓(xùn)練好的深度神經(jīng)網(wǎng)絡(luò)模型進(jìn)行測試，得到預(yù)測值，將預(yù)測值與模擬φ作差對比，如圖25～圖28。由圖25～圖28可以看出，在實際的物理問題中，深度神經(jīng)網(wǎng)絡(luò)模型依然表現(xiàn)良好。

(11)

(12)

×10-3

圖25 模擬ρFig.25 Simulation ρ

圖26 模擬φFig.26 Simulation φ

圖27 網(wǎng)絡(luò)預(yù)測值Fig.27 Network Prediction

圖28 網(wǎng)絡(luò)預(yù)測值與模擬φ的誤差

我們將同樣的數(shù)據(jù)用于快速傅里葉變換法求解，以做精度對比。如圖29～圖30。

×10-2

圖29 快速傅里葉變換求解值Fig.29 FFT Prediction

圖30 快速傅里葉變換求解值與模擬φ的誤差Fig.30 Error between FFT prediction and simulation φ

圖28和圖30對比可以看出，在同樣采樣率的情況下，深度神經(jīng)網(wǎng)絡(luò)模型的精度高于快速傅里葉變換法，滿足精度要求。

5.5 與快速傅里葉變換法及有限差分法的計算時間對比

實驗生成網(wǎng)格化尺寸逐漸增長的數(shù)據(jù)，分別對比測試訓(xùn)練好的深度神經(jīng)網(wǎng)絡(luò)模型和快速傅里葉變換法以及有限差分法的計算時間。測試結(jié)果如圖31～圖32。

實驗結(jié)果可以看出，深度神經(jīng)網(wǎng)絡(luò)方法求解泊松方程，隨著數(shù)據(jù)網(wǎng)格化尺寸增大，計算速度快于快速傅里葉變換法以及有限差分法。這歸功于兩點：(1)深度神經(jīng)網(wǎng)絡(luò)可以利用圖形處理器進(jìn)行快速計算；(2)從實驗結(jié)果圖可以看出，深度神經(jīng)網(wǎng)絡(luò)的計算時間趨勢為線性增長，而快速傅里葉變換法為NlogN倍增長，有限差分法更是N2倍增長。這個結(jié)論與前文卷積神經(jīng)網(wǎng)絡(luò)計算時間復(fù)雜度理論分析相互印證。

圖31 神經(jīng)網(wǎng)絡(luò)模型和快速傅里葉變換法計算時間對比。藍(lán)色代表神經(jīng)網(wǎng)絡(luò)模型，紅色代表快速傅里葉變換法

圖32 神經(jīng)網(wǎng)絡(luò)模型和有限差分法計算時間對比。藍(lán)色代表神經(jīng)網(wǎng)絡(luò)模型，黃色代表有限差分法

5.6 失敗神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)概述

這里描述一些試驗結(jié)果差的網(wǎng)絡(luò)結(jié)構(gòu)，數(shù)據(jù)生成方式和上述相同。第1種是純卷積網(wǎng)絡(luò)結(jié)構(gòu)，通道數(shù)為1-128-256-512-256-128，最后平均平方誤差收斂到0.000 3，誤差較大。在第1種網(wǎng)絡(luò)結(jié)構(gòu)的中間每個卷積層后插入歸一化層，效果變差，平均平方誤差只能收斂到0.14。加入一種網(wǎng)絡(luò)基礎(chǔ)模塊，對每次輸入進(jìn)行一次卷積，再將卷積后的輸出與輸入疊加作為下一層的輸入，下一層再通過卷積將維度變回初始輸入的維度，效果也不如純卷積網(wǎng)絡(luò)，平均平方誤差只能收斂到0.005。在純卷積網(wǎng)絡(luò)末端加入兩層全連接層，不僅參數(shù)量驟增，平均平方誤差也僅能達(dá)到0.01。將全連接層加到前端以及后端，平均平方誤差能到0.000 8，但參數(shù)量大，訓(xùn)練很慢，且效果不如簡單的純卷積網(wǎng)絡(luò)。將純卷積網(wǎng)絡(luò)中間兩層的步長改為2，提取更精簡的特征，并在后端的兩層將卷積改為反卷積進(jìn)行上采樣，此時平均平方誤差降到10-6，比純卷積網(wǎng)絡(luò)效果更好，但依然達(dá)不到精度要求。采用一次卷積一次反卷積迭代的網(wǎng)絡(luò)結(jié)構(gòu)，平均平方誤差只能到0.3。在純卷積網(wǎng)絡(luò)結(jié)構(gòu)上，加入模仿ResNet的網(wǎng)絡(luò)基礎(chǔ)模塊，對某次輸入進(jìn)行一次卷積，然后將輸入和輸出相加作為下一層的輸入，一次卷積，一次殘差模塊，連續(xù)進(jìn)行6次，平均平方誤差能達(dá)到10-5?；谝陨暇W(wǎng)絡(luò)結(jié)構(gòu)的試驗，最后綜合各種結(jié)構(gòu)的優(yōu)劣勢，我們設(shè)計了最終的深度神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)。

5.7 有趣實驗現(xiàn)象

我們所給的數(shù)據(jù)都設(shè)置了Dirichlet邊界條件的數(shù)據(jù)，嘗試不設(shè)置邊界條件會發(fā)生什么。用神經(jīng)網(wǎng)絡(luò)模型對不設(shè)置邊界的16 × 16數(shù)據(jù)進(jìn)行訓(xùn)練，訓(xùn)練過程和結(jié)果的精度與給定邊界條件時的情況相似。在同網(wǎng)格化尺寸下訓(xùn)練結(jié)果的相對誤差精度也可以達(dá)到10-4，但將訓(xùn)練的神經(jīng)網(wǎng)絡(luò)模型用于更大網(wǎng)格化尺寸的數(shù)據(jù)上時，就發(fā)生了有趣的現(xiàn)象，它只能預(yù)測到邊緣16寬度的范圍，中間部分幾乎為0，如圖33～圖40。

結(jié)合神經(jīng)網(wǎng)絡(luò)模型具有不同網(wǎng)格化數(shù)據(jù)的普適性分析，這樣的現(xiàn)象可能是因為訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)模型是局部迭代操作，同時又是從邊界開始向內(nèi)求解泊松方程造成的。

5.8 實驗總結(jié)

綜合以上實驗結(jié)果可以得出：深度神經(jīng)網(wǎng)絡(luò)模型具有快速求解PM-Tree方法中勢能的能力,精度不低于快速傅里葉變換法，速度優(yōu)于快速傅里葉變換法和有限差分法，特別是隨著無碰撞引力N體數(shù)值模擬的粒子數(shù)規(guī)模增加，深度神經(jīng)網(wǎng)絡(luò)模型的速度優(yōu)勢更加明顯；在未經(jīng)訓(xùn)練的更大網(wǎng)格化尺寸數(shù)據(jù)上也能使用，具有可擴(kuò)展性。所以深度學(xué)習(xí)的方法在加速無碰撞引力N體數(shù)值模擬中是可行的。

圖33 模擬32 × 32 ρ

圖34 模擬32 × 32 φ

圖35 32 × 32預(yù)測值

圖36 預(yù)測值與模擬φ的相對誤差

圖37 模擬64 × 64 ρ

圖38 模擬64 × 64 φ

圖39 64 × 64預(yù)測值

圖40 預(yù)測值與模擬φ的相對誤差

6 結(jié) 論

本文在神經(jīng)網(wǎng)絡(luò)解各類方程的研究和幾種無碰撞引力N體數(shù)值模擬方法的基礎(chǔ)上，提出了用一個輔以殘差局部結(jié)構(gòu)的Encoder-Decoder深度神經(jīng)網(wǎng)絡(luò)模型求解PM-Tree方法中的勢能，代替常用的快速傅里葉變換求解；驗證了深度神經(jīng)網(wǎng)絡(luò)模型可以快速求解泊松方程，理論上深度神經(jīng)網(wǎng)絡(luò)的計算時間復(fù)雜度為O(N)；驗證了深度神經(jīng)網(wǎng)絡(luò)模型計算耗時低于快速傅里葉變換法以及有限差分法；驗證了在相同離散化采樣率的情況下，深度神經(jīng)網(wǎng)絡(luò)模型求解精度不低于快速傅里葉變換法；驗證了深度神經(jīng)網(wǎng)絡(luò)具有可擴(kuò)展性，在小網(wǎng)格化尺寸數(shù)據(jù)集上訓(xùn)練的模型可以應(yīng)用于大網(wǎng)格化尺寸數(shù)據(jù)集。并且深度神經(jīng)網(wǎng)絡(luò)模型可以更好地利用圖形處理器的計算資源進(jìn)行并行處理，在粒子數(shù)規(guī)模較大時，速度優(yōu)勢更加明顯。綜合以上驗證結(jié)果，我們提出使用深度神經(jīng)網(wǎng)絡(luò)代替快速傅里葉變換法求解PM-Tree方法中的勢能，可以加速PM-Tree方法中主要的耗時步驟，即可以加速無碰撞引力N體數(shù)值模擬，驗證了深度學(xué)習(xí)方法在無碰撞引力N體數(shù)值模擬中的可行性，為將來實際應(yīng)用大規(guī)模無碰撞引力N體數(shù)值模擬打下基礎(chǔ)。

7 未來與展望

未來工作主要有3方面。(1)拓展解泊松方程由二維到三維(只需簡單地將現(xiàn)有神經(jīng)網(wǎng)絡(luò)模型中的二維卷積改為三維卷積即可)；(2)加入迭代求解過程，實現(xiàn)時變的模擬；(3)將模型嵌入無碰撞引力N體數(shù)值模擬的PM-Tree方法中實際應(yīng)用。

致謝：感謝劉冬明、翟順誠在工作中提出的意見；感謝顧建宇對文章部分內(nèi)容的修正。