基于深度學(xué)習(xí)的無(wú)碰撞引力N體數(shù)值模擬的可行性研究 *

趙梓成，龍 潛，董小波，孟潤(rùn)宇，鐘詩(shī)言，諶俊毅，向梓琨

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

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

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

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

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

1 相關(guān)工作

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

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

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

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

3 方 法

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

3.1 N體數(shù)值模擬

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

(1)

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

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

3.2 平均場(chǎng)法

平均場(chǎng)法是出現(xiàn)較早的一種N體數(shù)值模擬方法。用統(tǒng)計(jì)方法描述引力N體問(wèn)題時(shí)，系統(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)，平均場(chǎng)近似，即無(wú)碰撞N體系統(tǒng)。平均場(chǎng)近似中，粒子受到的作用力取決于分布函數(shù)，不必計(jì)算粒子兩兩之間的引力相互作用，這時(shí)我們可以用單粒子分布函數(shù)f(x,v,t)構(gòu)造系統(tǒng)的平均場(chǎng)，即在d3xd3v的六維相空間中發(fā)現(xiàn)一個(gè)粒子的概率為f(x,v,t)d3xd3v。這個(gè)分布函數(shù)的演化可以由無(wú)碰撞玻爾茲曼方程描述[1]：

(2)

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

(3)

(4)

3.3 PM-Tree

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

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

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

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

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

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

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

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

因?yàn)樵赑M-Tree方法步驟中，Tree代碼是和PM方法并行計(jì)算的，并且每次Tree代碼計(jì)算所需時(shí)間不大于PM方法每次計(jì)算時(shí)間的一半，所以PM-Tree方法耗費(fèi)的時(shí)間主要為PM步驟耗費(fèi)的時(shí)間。在PM方法中，求解勢(shì)能是主要耗時(shí)的計(jì)算步驟，而勢(shì)能需要求解泊松方程得到。所以加速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ì)算機(jī)數(shù)值求解泊松方程，首先將泊松方程轉(zhuǎn)換為代數(shù)方程，即差分方程。二維泊松方程的差分方程表示為

(5)

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

(6)

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

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

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

(7)

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

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

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ù)測(cè)值和勢(shì)能真值的平均平方誤差

(8)

使用最大差值絕對(duì)值來(lái)評(píng)估誤差

(9)

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

(10)

4 訓(xùn)練模型

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

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

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

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

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

圖6 對(duì)圖2的二階差分

圖7 對(duì)圖3的二階差分

圖8 對(duì)圖4的二階差分

圖9 對(duì)圖5的二階差分

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

我們?cè)赑ytorch中實(shí)現(xiàn)深度神經(jīng)網(wǎng)絡(luò)模型。每次訓(xùn)練使用100萬(wàn)對(duì)數(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)練成果及測(cè)試

5.1 訓(xùn)練成果

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

×10-4

圖10 模擬ρFig.10 Simulation ρ

圖11 模擬φFig.11 Simulation φ

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

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

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

5.2 測(cè)試結(jié)果

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

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

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

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

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

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

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

由圖15～圖16可以看出，最大誤差趨勢(shì)漸平，不隨網(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ù)測(cè)值

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

圖21 模擬1 024 × 1 024 ρ

圖22 模擬1 024 × 1 024 φ

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

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

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

5.4 等時(shí)勢(shì)問(wèn)題測(cè)試

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

(11)

(12)

×10-3

圖25 模擬ρFig.25 Simulation ρ

圖26 模擬φFig.26 Simulation φ

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

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

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

×10-2

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

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

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

5.5 與快速傅里葉變換法及有限差分法的計(jì)算時(shí)間對(duì)比

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

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

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

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

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

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

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

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

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

5.8 實(shí)驗(yàn)總結(jié)

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

圖33 模擬32 × 32 ρ

圖34 模擬32 × 32 φ

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

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

圖37 模擬64 × 64 ρ

圖38 模擬64 × 64 φ

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

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

6 結(jié) 論

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

7 未來(lái)與展望

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

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