蔣子超 江俊揚 姚清河 楊耿超

(中山大學航空航天學院力學系,廣州 510006)

引言

近年來,人工神經(jīng)網(wǎng)絡(luò)(artificial neural networks,ANN),尤其是深度神經(jīng)網(wǎng)絡(luò)(deep neural networks,DNN)因其在處理和預測高維復雜系統(tǒng)方面的優(yōu)異能力而被廣泛應(yīng)用于自然語言處理[1]、計算機視覺[2]與圖像處理[3]等領(lǐng)域.此外,深度神經(jīng)網(wǎng)絡(luò)憑借其在GPU 等異構(gòu)硬件架構(gòu)上的高并行效率,已經(jīng)成為當前科學計算領(lǐng)域的一個熱點研究方向,甚至為偏微分方程的實時求解提供了可行的實現(xiàn)途徑[4].目前國內(nèi)外學者已經(jīng)提出了多種基于深度神經(jīng)網(wǎng)絡(luò)的數(shù)值方法,如Ray 和Hesthaven[5]提出了基于深度神經(jīng)網(wǎng)絡(luò)的守恒定律輔助求解方法;Chan 和Wang 等[6-7]采用深度神經(jīng)網(wǎng)絡(luò)解決了多尺度問題.

根據(jù)萬能逼近定理[8-9],要使深度神經(jīng)網(wǎng)絡(luò)可以有效地逼近待擬合算子,該算子必須是分段連續(xù)的.然而偏微分方程,尤其是非線性偏微分方程的解通常正則性較低,深度神經(jīng)網(wǎng)絡(luò)對解的直接逼近往往會有精度不足的問題,且在外推預測時還可能受到可解釋性方面的質(zhì)疑.因此,目前已發(fā)表的研究主要集中在解的回歸分析和半解析方法上.例如,Owhadi[10]首先提出了物理信息神經(jīng)網(wǎng)絡(luò),基于訓練過程實現(xiàn)變分問題的優(yōu)化,是一種具有代表性的半解析方法.Raissi 等[11-12]在機器學習法的基礎(chǔ)上采用高斯過程擬合線性算子,并進一步將該方法推廣到非線性算子的回歸.此外,在復雜偏微分方程的輔助計算方面,深度神經(jīng)網(wǎng)絡(luò)也有諸多應(yīng)用,例如,謝晨月等[13]提出了基于深度神經(jīng)網(wǎng)絡(luò)的湍流大渦模擬方法.

在傳統(tǒng)數(shù)值模擬中,大型線性方程組通常采用Gauss-Seidel 等直接求解法或基于Krylov 子空間的BiCG 和ICPCG 等迭代方法[14-16]求解.這些方法雖保證了數(shù)值精度,但也在內(nèi)存消耗與計算時間方面帶來了巨大的挑戰(zhàn),成為目前大規(guī)模數(shù)值模擬領(lǐng)域的關(guān)鍵瓶頸之一[17].因此,通過深度神經(jīng)網(wǎng)絡(luò)加速求解線性方程組被認為是一種很有前景的研究思路[18-28].例如,Wang 等[18-19]將簡單結(jié)構(gòu)的神經(jīng)網(wǎng)絡(luò)直接應(yīng)用于求解線性方程組,并在后續(xù)研究中繼續(xù)完善了這一方法[20-23];多學者[24-27]也實現(xiàn)了將循環(huán)神經(jīng)網(wǎng)絡(luò)應(yīng)用于線性矩陣方程系統(tǒng)的求解中.

但是,在上述研究中,深度神經(jīng)網(wǎng)絡(luò)模型的精度成為了其在實際數(shù)值計算中的瓶頸[28-29].例如,Xiao等[30]提出了一種基于卷積神經(jīng)網(wǎng)絡(luò)的泊松方程求解器,并將其算法應(yīng)用于水和煙的實際流體模擬中.該算法相對經(jīng)典算法具有顯著的加速效果,在可視化效果上仿真結(jié)果也與經(jīng)典方法較為接近.但是,該方法的數(shù)值相對誤差大約在0.01 的量級,難以應(yīng)用于精度要求更高的數(shù)值仿真中.

深度神經(jīng)網(wǎng)絡(luò)模型無法實現(xiàn)更高的預測精度通常有兩方面的原因:(1)訓練過程停滯于局部最優(yōu)或梯度消失等問題;(2)模型深度(層數(shù))不足等問題導致的模型參數(shù)過少[29].因此,為了提升深度神經(jīng)網(wǎng)絡(luò)求解的精度,通常需要增加網(wǎng)絡(luò)模型的深度與參數(shù)數(shù)量.但隨著層數(shù)的加深,網(wǎng)絡(luò)退化問題成為了一個新的挑戰(zhàn).因此,He 等[31]首先提出了殘差網(wǎng)絡(luò)(Res-Net)結(jié)構(gòu):Res-Net 將層間短接結(jié)構(gòu)引入了深度神經(jīng)網(wǎng)絡(luò)網(wǎng)絡(luò)中,實現(xiàn)了一個天然的恒等映射,以保證誤差信息可以在深度神經(jīng)網(wǎng)絡(luò)中各隱藏層間的傳播.目前,Res-Net 已經(jīng)被大量的實際應(yīng)用證明了其解決網(wǎng)絡(luò)退化與梯度消失問題時的有效性[32-35].在數(shù)值計算領(lǐng)域,Qin 等[36]將Res-Net 引入了非線性的自恰系統(tǒng)的模擬中.

雖然萬能逼近定理[8-9]在理論上證明了當神經(jīng)元數(shù)量充分大時,深度神經(jīng)網(wǎng)絡(luò)模型可以實現(xiàn)任意精度直接預測線性方程組的解.但實際上,受制于深度神經(jīng)網(wǎng)絡(luò)模型的泛化能力等因素,對于一般的線性方程組,通過單次深度神經(jīng)網(wǎng)絡(luò)模型預測而直接獲得足夠精確(誤差量級小于離散誤差)的解往往是不切實際的.受傳統(tǒng)迭代算法的啟發(fā),本文提出一種迭代校正方法來提升深度神經(jīng)網(wǎng)絡(luò)模型預測解的精度,且該方法不需要過大的深度神經(jīng)網(wǎng)絡(luò)模型.該方法通過計算每一次深度神經(jīng)網(wǎng)絡(luò)預測解的殘差,將殘差作為右端項輸入深度神經(jīng)網(wǎng)絡(luò)模型計算校正量,如此循環(huán)直至達到收斂準則,由此以校正迭代的方式逐步減少數(shù)值的誤差.

深度神經(jīng)網(wǎng)絡(luò)模型中,輸入變量通常為向量或向量組,因此,為兼顧線性方程組的向量化表示與算法的適用范圍,本文所提出的基于深度神經(jīng)網(wǎng)絡(luò)的求解算法主要針對多對角矩陣.在數(shù)值計算中,多對角矩陣被廣泛應(yīng)用于有限差分法與譜方法中,此外,多對角矩陣具有天然的向量化表示方法,可直接作為深度神經(jīng)網(wǎng)絡(luò)模型的輸入層.

本工作將采用基于Res-Net 的深度神經(jīng)網(wǎng)絡(luò)模型與校正迭代方法,開發(fā)一種兼具高計算效率與高數(shù)值精度的線性方程組求解算法.此外,本文將研究該算法與有限差分法的結(jié)合,以實現(xiàn)偏微分方程數(shù)值的加速求解,并驗證該求解算法的計算效率與數(shù)值精度,以期為偏微分方程數(shù)值計算提供新型的高效求解方法.

1 求解算法

1.1 神經(jīng)網(wǎng)絡(luò)求解方法

考慮到多對角線性系統(tǒng)Ax=b,其左側(cè)矩陣A可以式(1)的形式分解

其中,d igi(ai,1,ai,2,···,ai,(n?i))表示將向量(ai,1,ai,2,···,ai,(n?i))的元素排列在第i條對角線上,其他位置均為0 的矩陣,特別地,d ig0(·) 表示主對角線;當i>0 時,digi(·)表示主對角線上方的對角線,當i<0時,digi(·)表示主對角線下方的對角線.該表示方法即將任意一個線性系統(tǒng)用數(shù)個向量進行表示,這一表示方法對于多對角矩陣與類多對角矩陣有天然的兼容性.實際上,本文提出的算法可適用于多種具有類似結(jié)構(gòu)的矩陣,如五對角矩陣等,在本文中,由于三對角矩陣具有更廣泛的應(yīng)用場景,且為了簡化表述,故選用了三對角矩陣作為示例以說明.

因為任意的三對角矩陣可以轉(zhuǎn)變成主對角線處的所有元素均為1 的三對角矩陣,所以任意三對角線性方程組均可用3 個向量來表示:左對角線向量l,右對角線向量r與右側(cè)向量b,而這3 個向量構(gòu)成的向量組即可直接作為神經(jīng)網(wǎng)絡(luò)模型的輸入層.

實際上,本文提出的方法可以看作是對算子S(A,b)=A?1b的回歸,相應(yīng)地,選擇均方誤差作為網(wǎng)絡(luò)的損失函數(shù)L(x?;x),以實現(xiàn)深度神經(jīng)網(wǎng)絡(luò)模型對算子S的擬合

其中M為訓練集的大小.

本文中,采用的訓練集由服從高斯分布的隨機向量構(gòu)成,即訓練集輸入為,其中,ρ為訓練集的標準差,D為待求解方程組的維數(shù),均為生成訓練集時需給定的參數(shù).

1.2 殘差網(wǎng)絡(luò)

Res-Net 通過在隱藏層之間建立“短路”連接(稱作殘差塊,如圖1 所示),從而使得網(wǎng)絡(luò)可以直接近似恒等算子,即將原本的待擬合算子轉(zhuǎn)換為對輸入與輸出層間的殘差的擬合[36-38].

圖1 殘差網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig.1 Schematic of the Res-Net

在圖1 中,傳入隱藏層的輸入向量在完成激活函數(shù)(本文選用ReLU 激活函數(shù))的計算后,會再與輸入向量自身相加而得到殘差塊的輸出值.由此可見,殘差塊結(jié)構(gòu)使網(wǎng)絡(luò)可以直接近似恒等算子,即將原本的待擬合算子轉(zhuǎn)換為對輸入與輸出層間的殘差的擬合[36-38].目前,Res-Net 已經(jīng)被證明是一種可靠的避免網(wǎng)絡(luò)退化問題的手段,且被廣泛應(yīng)用于多種成熟的深度神經(jīng)網(wǎng)絡(luò)應(yīng)用中[32-35].

除殘差塊外,還引入到另一個覆蓋整個網(wǎng)絡(luò)的短連接,以進一步逼近線性系統(tǒng)的解,即講深度神經(jīng)網(wǎng)絡(luò)模型的輸出替換為 (I?A?1)b.這一變換將使輸出向量的均值為0,相當于通過正則化的方法優(yōu)化了采用ReLU 激活函數(shù)的深度神經(jīng)網(wǎng)絡(luò)模型的訓練效果.

為了評估通過Res-Net 的實際優(yōu)化效果,測試了一個數(shù)值實驗進行比較,圖2 展示了無Res-Net 時深度神經(jīng)網(wǎng)絡(luò)模型訓練損失函數(shù)與有Res-Net 模型的損失函數(shù)的比值.

圖2 求解不同大小線性方程時Res-Net 與非Res-Net 的損失函數(shù)對比Fig.2 Comparison of the loss function for different sizes of the linear equations with Res-Net and no Res-Net

從圖2 中可以看出,采用Res-Net 后,對于任意大小的方程組(深度神經(jīng)網(wǎng)絡(luò)模型),模型的損失函數(shù)可降低至不采用Res-Net 時模型損失函數(shù)的1/5000左右,這說明Res-Net 可以顯著提升訓練后的模型預測精度.

1.3 迭代校正方法

記深度神經(jīng)網(wǎng)絡(luò)模型的預測解為S(A,b),則迭代算法如圖3 所示.

圖3 迭代算法示意圖Fig.3 Schematic of the iteration algorithm

在圖3 中,參數(shù) ε 為迭代收斂條件;運算符 ? 表示兩個向量的對應(yīng)元素相乘.每次迭代中的向量m表示當前迭代步中Ari的值.在m的計算中,三對角矩陣和向量的乘積可以拆分成3 個向量內(nèi)積運算的求和,使計算復雜度和內(nèi)存使用量最小.

深度神經(jīng)網(wǎng)絡(luò)模型擬合的目標算子S′(A,b)=A?1b對于常數(shù)k應(yīng)具備性質(zhì)S′(A,kb)=kS′(A,b),受此啟發(fā),本文引入的迭代算法額外引入了一個放大因子(圖3 中的 α)用于進一步提升求解的精度.

為了論證該迭代算法可以有效地提升深度神經(jīng)網(wǎng)絡(luò)的求解精度,圖4 展示了對于給定方程組在設(shè)置不同的放大因子時迭代過程中殘差的變化情況以及對不同的放大因子,迭代最終收斂時的殘差大小.

圖4 (a)不同放大因子時殘差隨迭代的變化情況與(b)放大因子對收斂時殘差的影響Fig.4 (a) Variation of residuals with iterations for different amplification factors and (b) effect of amplification factors on convergence residuals

由圖4 可見,迭代算法可以有效地降低求解的誤差,對給定的線性方程組,殘差可被縮小至迭代初始時的 10?5倍.此外,求解精度隨放大因子 α 增大而提高,但是存在一個確定的閾值,當超過該值時,增大放大因子不會使精度進一步提升.

綜上所述,本文提出的求解器的體系結(jié)構(gòu)可以由圖5 來演示,它包括帶Res-Net 的深度神經(jīng)網(wǎng)絡(luò)和修正迭代.

圖5 算法求解器架構(gòu)示意圖Fig.5 The schematic of the linear solver

在圖5 中,DiagMatMult 子程序表示圖3 中提到的多對角矩陣與向量的乘積,Concatenate 子程序表示對輸入向量進行拼接均,帶上下標的b與w分別為神經(jīng)元的偏置與權(quán)重.

圖5 的下半部分為包含了Res-Net 的深度神經(jīng)網(wǎng)絡(luò)模型,上半部分為校正迭代過程:在每個迭代步中,先由深度神經(jīng)網(wǎng)絡(luò)模型得到預測解x?,再計算x?的殘差,再作為新的右端向量傳入深度神經(jīng)網(wǎng)絡(luò)進行下一輪的校正迭代(如圖3 所示).

由算法架構(gòu)可見,該算法除表示線性系統(tǒng)的向量組(即深度神經(jīng)網(wǎng)絡(luò)的輸入向量)外,不需要引入額外的矩陣結(jié)構(gòu),可完全依托向量化計算實現(xiàn),相對傳統(tǒng)方法可顯著減少空間消耗.此外,該方法可同時求解多個具有相同結(jié)構(gòu)的線性方程組,這對于基于交替方向隱式等方法的數(shù)值方法具有顯著的效率優(yōu)勢.

2 數(shù)值實驗和驗證

2.1 求解效率測試

為了驗證本文所提出的基于深度神經(jīng)網(wǎng)絡(luò)的求解算法在效率方面的優(yōu)勢以及數(shù)值計算中的應(yīng)用前景,采用本文提出的求解算法和傳統(tǒng)方法(LU 分解法等)分別對線性方程組進行求解以測試其求解效率.

本文使用Google[39]的TensorFlow 2.0 平臺來開發(fā)本工作中的所有計算程序.圖6 展示了本文新算法與NumPy 與SciPy 中線性求解器之間的計算效率比較.NumPy 默認調(diào)用LAPACK 庫中內(nèi)置的LU分解求解線性方程組;SciPy 基于BLAS 和LAPACK的標準線性代數(shù)算法.NumPy 與SciPy 兩者均能在多核平臺上實現(xiàn)最大CPU 核數(shù)(本例為28 核)的并行求解,而TensorFlow 平臺天然支持GPU 的加速.圖中測試采用的平臺是Intel Xeon-W3175X CPU(28 核,3.1 GHz)和Nvidia RTX2080Ti GPU.在測試中,不同求解器求解精度均為 10?9,分別測試了多種具有代表性的迭代法,測試的矩陣大小范圍為20×20～4000×4000之間,每次測試均求解5000 個重新隨機生成的矩陣,其高斯分布參數(shù)與訓練集一致,并計算其平均求解時間,結(jié)果如圖6 所示.

圖6 求解不同大小線性方程各種算法的求解時間Fig.6 Solving times of DNN and conventional algorithms for different size of linear equations

圖6 中,LU 算法由于雖然在小矩陣的求解中效率比其他方法高,但是由于其計算復雜度較高,隨著矩陣大小的增加其相對效率也逐漸變低.而本文所提出的基于深度神經(jīng)網(wǎng)絡(luò)的算法隨著矩陣大小增加,求解時間增長速度明顯較慢,其求解時間對矩陣大小的敏感度明顯低于BiCG 和GMRES 等傳統(tǒng)方法.由此,該算法計算復雜度較經(jīng)典算法低,因此在大規(guī)模線性系統(tǒng)中更具有優(yōu)勢.

2.2 線性偏微分方程算例

本節(jié)將以二維熱傳導方程為例,將本文所提出的基于深度神經(jīng)網(wǎng)絡(luò)的求解方法結(jié)合交替方向隱式方法實現(xiàn)對實際偏微分方程進行數(shù)值求解.交替方向隱式法在基于有限差分法的數(shù)值算法中得到了廣泛的應(yīng)用,它可以將多對角矩陣轉(zhuǎn)換為三對角矩陣,這與上面介紹的網(wǎng)絡(luò)具有良好相容性.

考慮在矩形區(qū)域 Ω={(x,y)|0

其中,常系數(shù) β 取值為0.1.本文所采用的交替方向隱式格式如式(4)所示,其空間精度為2 階,時間精度為1 階

在該方法中,每個時間步長被分割成兩個半步長,如式(4)所示,一個三對角線性方程可以代表每個半步,因此,本文所提出的求解方法可以直接求解該三對角線性方程.

選取時間步長為0.001,兩個方向上的空間步長均為0.01,方程(3)的求解結(jié)果如圖7 所示.圖7 分別展示了在第50,第100 與第200 時間步時中心截面(y=0.5)處經(jīng)典方法(LU 算法)和本文提出的算法的計算結(jié)果(圖7(a))及其差值分布(圖7(b)).

圖7 中,隨計算時間逐漸推進,由本文提出的基于深度神經(jīng)網(wǎng)絡(luò)的求解算法所得結(jié)果在整個計算區(qū)域內(nèi)的誤差均保持在 10?5量級,可以證明其量級不高于與式(4)中格式的截斷誤差.由此,證明了本文所提出的算法對于二維線性方程的求解具有較高的數(shù)值精度.

圖7 (a)二維熱傳導方程的數(shù)值求解結(jié)果與(b)相對經(jīng)典方法的求解誤差Fig.7 (a) Numerical solution of the 2D heat conduction equation and (b) the error distribution

圖7 (a)二維熱傳導方程的數(shù)值求解結(jié)果與(b)相對經(jīng)典方法的求解誤差 (續(xù))Fig.7 (a) Numerical solution of the 2D heat conduction equation and (b) the error distribution (continued)

2.3 非線性偏微分方程方程算例

Burgers 方程是非線性聲學、氣體動力學、流體力學等各種應(yīng)用場景中的基本偏微分方程之一[40],最早由Bateman[41]提出,在該算例中,本文將利用所提出算法對一維Burgers 方程進行數(shù)值求解.

考慮如式(5)所示的非線性Burgers 方程

其中,β 為擴散系數(shù),本算例中,取 β=0.025.

該問題存在如式(6)所示的解析解.

選擇如式(7)所示的二階精度的差分格式

其中,常數(shù) α=Δt/(4Δx),β ?=βΔt/(2Δx2).本算例中,時間步長Δt=0.001,空間步長Δx=0.01.

與2.2 中算例類似,圖8(a)分別展示了第50、第100、第200 時間步時的計算結(jié)果.特別地,圖8(b)中給出了經(jīng)典線性求解器與本文提出的方法相對解析解的誤差對比.

圖8 不同時刻下(a)一維Burges 方程的數(shù)值求解結(jié)果與(b)求解誤差Fig.8 Numerical solution of (a) the 1D Burgers equation and (b) the error distribution

圖8 不同時刻下(a)一維Burges 方程的數(shù)值求解結(jié)果與(b)求解誤差 (續(xù))Fig.8 Numerical solution of (a) the 1D Burgers equation and (b) the error distribution (continued)

由圖8 可見,隨計算時間推進,由本文提出的求解算法所得結(jié)果在整個計算區(qū)域內(nèi)的誤差均保持在10?5量級,與2.2 節(jié)中的線性方程組算例保持一致,由此說明本算法可適應(yīng)多種具有類似結(jié)構(gòu)的差分方程組求解.

此外,由圖8(b)中結(jié)果可見,本文所提出的算法求解精度與經(jīng)典求解方法的精度幾乎一致,均為10?5量級,即不會引入高于截斷誤差量級的計算誤差.與之相對地,其他大多數(shù)基于神經(jīng)網(wǎng)絡(luò)的求解器的誤差大小范圍在 10?1到10?3之間[28,42].因此,作為一種基于深度神經(jīng)網(wǎng)絡(luò)的求解方法,本文所提出的求解算法在顯著提升計算效率的情況下,仍保證了高計算精度.

3 結(jié)語

本文結(jié)合深度神經(jīng)網(wǎng)絡(luò)與Res-Net 架構(gòu),引入了修正迭代過程,提出了一種針對多對角線性系統(tǒng)的新型求解算法,并將其應(yīng)用于實際的偏微分方程的有限差分法求解中.

(1) 基于深度神經(jīng)網(wǎng)絡(luò),本文提出的新算法在GPU 等混合平臺上具有較高的計算效率和原生硬件兼容性.與經(jīng)典的求解方法相比,本文的算法可同時求解多個線性方程組而無需引入額外的數(shù)據(jù)處理過程,因此具有更高的求解效率.

(2)通過在算法中加入了Res-Net 架構(gòu),本文的算法避免了深層網(wǎng)絡(luò)的退化問題,從而提高了新算法的精度.在此基礎(chǔ)上,本文引入了修正迭代與放大因子,進一步降低誤差,使深度神經(jīng)網(wǎng)絡(luò)求解誤差量級低于一般差分格式的離散誤差.

(3)相對傳統(tǒng)方法,本文的深度神經(jīng)網(wǎng)絡(luò)求解線性方程組的另一優(yōu)勢在于,該方法可以同時求解多個獨立的方程組,省略了矩陣構(gòu)建與循環(huán)過程,可顯著提升交替方向隱式等方法的求解效率.

(4)由于多對角線性方程組廣泛存在于多種差分方法中,本文的方法可以被應(yīng)用于多種偏微分方程的求解中,本文給出了一維非線性方程與二維線性方程的算例,論證了該方法應(yīng)用的廣泛性.