自動微分轉(zhuǎn)換系統(tǒng)及其應(yīng)用

1、自動微分轉(zhuǎn)換系統(tǒng)及其應(yīng)用    摘要    自動微分轉(zhuǎn)換系統(tǒng)（DFT）由LASG和LSEC聯(lián)合研制開發(fā)，目前已擁有成熟的版本。本文對DFT系統(tǒng)的功能、特色及其基本應(yīng)用作了全面的介紹，并給出了一些頗具說服力的數(shù)值試驗(yàn)結(jié)果。同時，本文提出了統(tǒng)計準(zhǔn)確率評價的概念，這對評價一類自動微分工具及其微分模式代碼的可靠性與有效性提供了一種客觀的尺度。最后，本文還詳細(xì)討論了運(yùn)用切線性模式求解雅可比矩陣的問題，給出了求解初始輸入矩陣的有效算法。關(guān)鍵詞 自動微分  切線性模式  數(shù)據(jù)相關(guān)分析  統(tǒng)計準(zhǔn)確率1.引言計算

2、微分大致經(jīng)歷了從商微分，符號微分，手寫代碼到自動微分幾個階段。與其它幾種微分方法相比，自動微分具有代碼簡練、計算精度高及投入人力少等優(yōu)點(diǎn)。自動微分實(shí)現(xiàn)的基本出發(fā)點(diǎn)是：一個數(shù)據(jù)相對獨(dú)立的程序?qū)ο螅Ｊ健⑦^程、程序段、數(shù)值語句乃至數(shù)值表達(dá)式），無論多么復(fù)雜，總可以分解為一系列有限數(shù)目的基本函數(shù)（如sin、exp、log）和基本運(yùn)算操作（加、減、乘、除、乘方）的有序復(fù)合；對所有這些基本函數(shù)及基本運(yùn)算操作，重復(fù)使用鏈?zhǔn)角髮?dǎo)法則，將得到的中間結(jié)果自上而下地做正向積分就可以建立起對應(yīng)的切線性模式，而自下而上地做反向積分就可以建立起對應(yīng)的伴隨模式1?；谧詣游⒎址椒ǖ玫降那芯€性模式和伴隨模式，在變分資料同化

3、2、系統(tǒng)建模與參數(shù)辨識3、參數(shù)的敏感性分析4、非線性最優(yōu)化以及數(shù)值模式的可預(yù)測性分析5等問題中有著十分廣泛的應(yīng)用。迄今為止，已有數(shù)十所大學(xué)和研究所各自開發(fā)了能夠用于求解切線性模式的自動微分系統(tǒng)，比較典型的有TAMC系統(tǒng)6、ADJIFOR系統(tǒng)7 和ODYSSEE系統(tǒng)8。在一些特定的運(yùn)用中，它們都是比較成功的，但在通用性和復(fù)雜問題的處理效率上還存在許多不足。通常，自動生成切線性模式的關(guān)鍵難題在于對象自身的強(qiáng)相關(guān)性，這給系統(tǒng)全局分析（如數(shù)據(jù)IO相關(guān)分析和數(shù)據(jù)依賴相關(guān)分析）和微分代碼的整體優(yōu)化都帶來了很多困難。同時，對于程序?qū)ο蟛豢蓪?dǎo)處的準(zhǔn)確識別和微分處理，至今仍還沒有一個統(tǒng)一而有效的算法。另外，最優(yōu)

4、或有效求解稀疏雅可比矩陣一直是衡量一個自動微分系統(tǒng)有效性的重要尺度。統(tǒng)計準(zhǔn)確率被我們視為評價一類自動微分工具及其微分模式代碼可靠性與有效性的重要尺度。其基本假設(shè)是：如果對于定義域空間內(nèi)隨機(jī)抽樣獲得的至多有限個n維初始場（或網(wǎng)格點(diǎn)），微分模式輸出的差分和微分逼近是成功的；那么對于定義域空間內(nèi)所有可能初始場（或網(wǎng)格點(diǎn)），微分模式輸出的差分和微分逼近都是成功的。微分模式統(tǒng)計準(zhǔn)確率評價的具體方法是：在所有隨機(jī)抽樣得到的初始場（或網(wǎng)格點(diǎn)）附近，當(dāng)輸入擾動逐漸趨向于機(jī)器有效精度所能表示的最小正值時，模式輸出的差分和微分之間應(yīng)該有足夠精度有效位數(shù)上的逼近。DFT系統(tǒng)具有許多優(yōu)點(diǎn)，它能夠完全接受用FORTRA

5、N 77語言編寫的源代碼，微分代碼結(jié)構(gòu)清晰，其微分處理能力與問題和對象的規(guī)模及復(fù)雜性無關(guān)。它基于YACC實(shí)現(xiàn)，具有很強(qiáng)的可擴(kuò)展性。DFT系統(tǒng)具有四個重要特色。它通過對象全局依賴相關(guān)分析，準(zhǔn)確求解雅可比矩陣的稀疏結(jié)構(gòu)，自動計算有效初始輸入矩陣，從而可以用較小的代價求得整個雅可比矩陣。同時，它可以自動生成客觀評價微分模式效率與可靠性的測試程序，對奇異函數(shù)做等價微分處理，并采用二元?dú)w約的方法，在語句級層次上實(shí)現(xiàn)微分代碼優(yōu)化。2.系統(tǒng)概況DFT系統(tǒng)主要由兩部分組成：微分代碼轉(zhuǎn)換和微分代碼評價，圖2.1。微分代碼轉(zhuǎn)換部分接受用戶輸入指令并自動分析對象模式，生成切線性模式代碼及其相關(guān)測試代碼，后者直接構(gòu)成

6、微分代碼評價系統(tǒng)的主體。微分代碼評價是DFT系統(tǒng)的一個重要特色。DFT系統(tǒng)的開發(fā)小組認(rèn)為，一個微分模式如果在可靠性、時間和存儲效率上沒有得到充分的驗(yàn)證，至少對實(shí)際應(yīng)用而言，它將是毫無意義的。                             原模式     

7、60;                                                 

8、60;                    切線性模式                              &#

9、160;                                                  &

10、#160;                                                 &

11、#160;                   統(tǒng)計評價結(jié)果     圖2.1  DFT系統(tǒng)結(jié)構(gòu)簡圖2.1 微分代碼轉(zhuǎn)換DFT系統(tǒng)是基于YACC在UNIX環(huán)境下開發(fā)的，其結(jié)構(gòu)圖2.2所示。通過DFT系統(tǒng)產(chǎn)生的切線性模式代碼成對出現(xiàn)，并在語句級程度上做了簡化，可讀性很強(qiáng)，如圖2.4。       &#

12、160;                                                   

13、                                                  

14、                                                   

15、;                                                  

16、;                         切線性模式                         

17、;                                                  

18、;                                          評價函數(shù)集圖2.2   微分代碼轉(zhuǎn)換     

19、;    微分代碼轉(zhuǎn)換部分從功能上分為四個部分：詞法分析，語義分析，對象復(fù)雜性及數(shù)據(jù)相關(guān)分析和微分代碼轉(zhuǎn)換。對于一組具有復(fù)雜數(shù)據(jù)相關(guān)的程序模式對象，通常需要系統(tǒng)運(yùn)行兩遍才能得到有效而可靠的微分代碼。這主要有兩方面的考慮：其一，根據(jù)對象的復(fù)雜性（如最大語句長度、最大變量維數(shù)、子過程或函數(shù)數(shù)目、子過程或函數(shù)內(nèi)最大變量數(shù)目等對象特征）選擇合適的系統(tǒng)參數(shù)以求最優(yōu)的運(yùn)行代價；其二，模式內(nèi)各子過程或函數(shù)之間以及一個子過程或函數(shù)內(nèi)往往具有很強(qiáng)的數(shù)據(jù)相關(guān)性，需要事先保存對象的相關(guān)信息并且在考慮當(dāng)前對象的屬性之前必須做上下文相關(guān)分析。    

20、                     圖2.3   PERIGEE源程序代碼                         

21、 圖2.4  DFT系統(tǒng)生成的切線性代碼2.2 微分代碼評價通常，評價一個編譯系統(tǒng)的性能有很多方面，如處理速度、結(jié)果代碼可靠性及質(zhì)量、出錯診斷、可擴(kuò)展和可維護(hù)性等。對于一類自動微分系統(tǒng)來說，由于軟件開發(fā)人力的局限以及對象模式的復(fù)雜多樣性,通過自動轉(zhuǎn)換得到的微分模式并非常常是有效而可靠的（即無論是在數(shù)學(xué)意義上還是在程序邏輯上應(yīng)與期待的理想結(jié)果一致），因而在微分模式被投入實(shí)際應(yīng)用前，往往需要投入一定的人力來對其做嚴(yán)格的分析測試。  對切線性模式做統(tǒng)計評價測試的主要內(nèi)容可以簡單敘述為：在網(wǎng)格化的模式定義域空間內(nèi)，選擇所有可能的網(wǎng)格點(diǎn)形成微分模式計算的初始場；在不同的網(wǎng)格點(diǎn)附近，隨

22、機(jī)選取至少 個線性無關(guān)的初始擾動，對每個擾動輸入分別進(jìn)行網(wǎng)格點(diǎn)逼近，統(tǒng)計考察模式輸出差分和微分在有效位數(shù)上的逼近程度。圖2.5描述了整個測試過程，它包含網(wǎng)格點(diǎn)數(shù)據(jù)隨機(jī)采樣（1）和網(wǎng)格點(diǎn)數(shù)據(jù)逼近（2）兩級循環(huán)。     圖2.5   切線性模式代碼的測試過程3.系統(tǒng)主要特色DFT系統(tǒng)并不是一個完整的FORTRAN編譯器，但它幾乎可以接受和處理所有FORTRAN 77編寫的源模式代碼，并且可以很方便地擴(kuò)展并接受FORTRAN 90編寫的源模式代碼。本節(jié)將著重介紹DFT系統(tǒng)（版本3.0）的以下幾個重要特色。3.1 結(jié)構(gòu)化的微分實(shí)現(xiàn)

23、 DFT系統(tǒng)采用標(biāo)準(zhǔn)化的代碼實(shí)現(xiàn)，切線性模式的擾動變量和基態(tài)值變量、微分計算語句和基態(tài)值計算語句總是成對出現(xiàn)，并具有清晰的程序結(jié)構(gòu)。微分代碼保持了原模式本身的結(jié)構(gòu)和風(fēng)格（如并行和向量特性、數(shù)據(jù)精度等），即語句到語句、結(jié)構(gòu)到結(jié)構(gòu)的微分實(shí)現(xiàn)。在奇異點(diǎn)或不可導(dǎo)處，DFT系統(tǒng)對微分?jǐn)_動采取簡單的清零處理，實(shí)踐證明這對抑制擾動計算溢出具有重要意義，但并不影響評價測試結(jié)果。3.2 全局?jǐn)?shù)據(jù)相關(guān)分析        DFT系統(tǒng)具有較強(qiáng)的數(shù)據(jù)相關(guān)分析能力，它包括全局?jǐn)?shù)據(jù)IO相關(guān)分析、全局?jǐn)?shù)據(jù)依賴相關(guān)分析、全局過程相關(guān)分析以及數(shù)據(jù)迭代相關(guān)分析幾

24、個不同方面。數(shù)據(jù)依賴相關(guān)與數(shù)據(jù)IO相關(guān)關(guān)系密切，但又存在根本不同。前者強(qiáng)調(diào)每個變量在數(shù)學(xué)關(guān)系上的依賴性；而后者描述了一個對象的輸入輸出特性，且具有相對性，即任何一個變量參數(shù)，無論它是獨(dú)立變量還是依賴變量，在數(shù)學(xué)意義上都可等價為一個既是輸入又是輸出的參數(shù)來處理。        DFT系統(tǒng)記錄所有過程參數(shù)的IO屬性表，通過深度遞歸相關(guān)計算，準(zhǔn)確計算每個過程參數(shù)的最終IO屬性。DFT系統(tǒng)通過對數(shù)據(jù)相關(guān)矩陣做模二和及自乘迭代計算（An+1= AnAn2）來完成數(shù)據(jù)的依賴相關(guān)分析，這種算法具有很好的對數(shù)收斂特性。DFT系統(tǒng)通過全局過程

25、相關(guān)分析的結(jié)果，自動生成模式的局部或整體相關(guān)引用樹結(jié)構(gòu)（如圖3.1），這對用戶分析復(fù)雜數(shù)值模式和微分評價測試都具有很好的指導(dǎo)作用。DFT系統(tǒng)還具有分析局部數(shù)據(jù)迭代相關(guān)和函數(shù)迭代相關(guān)的能力，這兩種形式的數(shù)據(jù)迭代相關(guān)是自動微分實(shí)現(xiàn)頗具挑戰(zhàn)的難題之一。       圖3.1  GPS Rayshooting模式的相關(guān)樹結(jié)構(gòu)片段3.3 自動生成測試程序        基于IO相關(guān)分析的結(jié)果，DFT系統(tǒng)自動生成微分測試代碼，分別對切線性模式的可靠性和運(yùn)行代

26、價做統(tǒng)計評價測試。特別地，DFT系統(tǒng)還可將任何模式參數(shù)都視為輸入輸出參數(shù)，生成在數(shù)學(xué)意義上等價的測試代碼，這樣處理的不利之處在于往往需要極高的存儲開銷。3.4 基于語句級的代碼優(yōu)化       目前，DFT系統(tǒng)僅僅具備局地優(yōu)化能力。在語句級微分實(shí)現(xiàn)上采用二元?dú)w約的方法對微分代碼進(jìn)行優(yōu)化是DFT系統(tǒng)的一個重要特色。根據(jù)右端表達(dá)式的乘法復(fù)雜性及含變元數(shù)目的不同，DFT系統(tǒng)采取不同的分解策略。二元?dú)w約的方法避免了微分計算中的許多冗余計算，在一些復(fù)雜的非線性表達(dá)式的微分計算中具有最小的計算代價，同時也非常適合于微分系統(tǒng)的軟件實(shí)現(xiàn)。同時，對于某

27、些特殊的運(yùn)算操作（除法、乘方）和特殊函數(shù)（如sqrt、exp），DFT系統(tǒng)較好地利用了基態(tài)值計算得到的中間結(jié)果，避免了微分實(shí)現(xiàn)中的冗余計算。4.系統(tǒng)應(yīng)用運(yùn)用自動微分工具得到的切線性模式，可以在無截斷誤差意義下求解函數(shù)的數(shù)值微分和導(dǎo)數(shù)、稀疏雅可比矩陣。同時這些結(jié)果在數(shù)值參數(shù)敏感性分析、非線性最優(yōu)化以及其它數(shù)值理論分析中有著非常重要的應(yīng)用。這里簡單介紹切線性模式的幾個基本應(yīng)用。4.1 符號導(dǎo)數(shù)和微分 如果輸入為數(shù)學(xué)關(guān)系式，DFT系統(tǒng)可以自動生成對應(yīng)的微分表達(dá)式和梯度，而與數(shù)學(xué)關(guān)系式的復(fù)雜程度無關(guān)。例如我們輸入關(guān)系式：       

28、;                    ，                               （1

29、）DFT系統(tǒng)將自動生成其符號微分形式及其梯度形式分別為                                ，              &

30、#160;   （2）4.2 數(shù)值導(dǎo)數(shù)和微分 切線性模式最基本的應(yīng)用就是在一定擾動輸入下求解輸出變量的擾動（響應(yīng)）。表4.1給出了DFT系統(tǒng)在對IAP 9L模式、GPS Rayshooting模式和GPS Raytrace模式三個數(shù)值模式做切線性化的具體應(yīng)用中，一些不同計算粒度、不同引用深度和不同程序風(fēng)格的核心子過程，以及它們的切線性模式在SGI 2000上運(yùn)行的統(tǒng)計評價測試結(jié)果，其中切線性模式的可靠性指標(biāo)都準(zhǔn)確到六個有效數(shù)字以上，在運(yùn)行時間、存儲開銷和代碼復(fù)雜性方面分別是原模式的兩倍左右，比較接近于理想的微分代價結(jié)果（1.5倍）。除了IAP 9L模式由于過于復(fù)雜僅做粗略統(tǒng)

31、計外，其余模式都用非注釋語句行數(shù)來表示各自的代碼復(fù)雜性。表4.1   DFT系統(tǒng)在三個數(shù)值模式中的統(tǒng)計評價測試結(jié)果 性能指標(biāo)對象模式  運(yùn)行時間（10-3秒） 存儲開銷（字節(jié)數(shù)） 代碼復(fù)雜性      原模式 切線性模式 原模式 切線性模式 原模式 切線性模式          Xyz2g 2.530 6.160 55

32、24 11048 55 89      IntCIRA 1.560 2.750 1334 2661 41 65      Dabel 0.035 0.072 60 120 27 49      LSS 8.300 17.50 669 1338 79&

33、#160;143      RP 42.40 85.10 3605 7210 22 38      Vgrad1 0.100 0.212 18564 36828 24 54      RefGr 43.00 86.00 718654 1437308 35 78 

34、;     LL2JK 0.626 1.350 2622 5244 22 32      RayFind4 62.70×103 125.4×103 9856 18212 111 179        EPSIMP 1.760 11.50 4455 8910&#

35、160;13 27      Hlimits 0.830 1.880 242577 484254 37 74      Int3sL 26.90 51.20 820029 1639458 46 90      MAKENCEP 1340 3920 722925 1445850&#

36、160;45 84      Curvcent 0.013 0.0385 27 54 27 54       DYFRM 3.800×103 7.250×103 5000* 9500* 161 279      PHYSIC 2.750×103 5.385&#

37、215;103 3000 5000* 1399*(含注釋行) 2826*(含注釋行)         適當(dāng)設(shè)置輸入擾動的初值，運(yùn)用切線性模式可以簡單求解輸出變量對輸入的偏導(dǎo)數(shù)。例如，對于一個含有 個輸入?yún)?shù)的實(shí)型函數(shù)                       

38、;                                                  

39、;         （3）這里設(shè) ， 。運(yùn)用DFT系統(tǒng)，可以得到對應(yīng)的切線性模式                       （4）其中 ， 為切線性模式的擾動輸入?yún)?shù)。可以通過以下辦法來求得偏導(dǎo)數(shù)：        &#

40、160;          （5）其中 。如果對于某個 既是輸入?yún)?shù)又是輸出參數(shù)，可以類似以下過程引用的辦法來處理。對于過程引用的情形，例如一個含有 個輸入?yún)?shù)的子過程                            

41、60;                           （6）其中 ， 為輸入?yún)?shù)； ， 為輸出參數(shù)； ， 既為輸入?yún)?shù)又為輸出參數(shù)。運(yùn)用DFT系統(tǒng)，可以得到對應(yīng)的切線性模式為              &

42、#160;                  （7）其中 ， ， ， 分別為切線性模式的微分?jǐn)_動輸入、輸出和輸入輸出參數(shù)?？梢酝ㄟ^以下輸入擾動設(shè)置并引用切線性模式（7）來求得偏導(dǎo)數(shù)：a) 設(shè)置 ； （ ， ）； （ ）可以同時求得 （ ）和 （ ），其中 。b) 設(shè)置 （ ）； ； （ ， ）可以同時求得 （ ）和 （ ），其中 。4.3 稀疏雅可比矩陣 運(yùn)用上節(jié)討論的方法來求解稀疏雅可比矩陣,具有極高的計算代價。例如，一個含

43、 個獨(dú)立和 個依賴參數(shù)的子過程，為求解整個雅可比矩陣就需要反復(fù)調(diào)用 次切線性模式，當(dāng) 相當(dāng)大時，這對許多實(shí)際的數(shù)值計算問題是不能接受的。事實(shí)上，如果雅可比矩陣的任意兩列（行）相互正交，那么可以通過適當(dāng)設(shè)置擾動輸入值，這兩列（行）的元素就可以通過一次引用切線性模式（伴隨模式）完全得到。設(shè) 和 分別為雅可比矩陣的行寬度和列寬度，即各行和各列非零元素數(shù)目的最大值，顯然有 ， 。這里介紹幾種常用的求解方法。正向積分  當(dāng) 時，通常采用切線性模式來計算雅可比矩陣。根據(jù)雅可比矩陣的稀疏結(jié)構(gòu)，適當(dāng)選擇右乘初始輸入矩陣，可以獲得接近 的計算時間代價。DFT系統(tǒng)采用一種逐列（行）求解的方法，來有效求解

44、右（左）乘初始輸入矩陣。其基本思路是：按照某種列次序考察雅可比矩陣的各列；考察當(dāng)前列中所有非零元素，并對這些非零元素所在行的行向量做類似模二和累加運(yùn)算（即將非零元素視為邏輯“1”， 零元素視為邏輯“0”），從而得到一個描述當(dāng)前列與各行存在“某種”相關(guān)的標(biāo)志向量（其元素都是“1”或“0”）；依據(jù)此標(biāo)志向量，就很容易得到一個與之正交的列初始向量 ，其中與當(dāng)前列序號對應(yīng)的元素設(shè)置為“1”，而與標(biāo)志向量中非零元素序號對應(yīng)的元素設(shè)置為“0”， 與標(biāo)志向量中非零元素序號對應(yīng)的元素設(shè)置為“-1”，顯然，該列初始向量是唯一的，并且對應(yīng)著當(dāng)前右乘初始輸入矩陣的最后一列；逐一考察已求解得到的列初始向量，如果某列初

45、始向量與當(dāng)前求解得到的列初始向量按下面定義的乘法（見過程4）正交，那么這兩列就可以合并，即將當(dāng)前列初始向量中非“-1”的元素按照對應(yīng)關(guān)系分別賦值給該初始向量，并從記錄中刪除當(dāng)前列初始向量；重復(fù)以上過程，繼續(xù)按照給定列次序考察雅可比矩陣的“下一列”。不難說明，按照不同列次序求解得到的右乘初始輸入矩陣可能不同。其中逐列求解右乘初始輸入矩陣的過程可以簡單敘述為：1）將右乘初始輸入矩陣 所有元素的初值均設(shè)置為 ， ， 。 。2）如果 ，轉(zhuǎn)6）。否則，如果雅可比矩陣 的第 列中的所有元素均為 ， ，重復(fù)2）的判斷。否則轉(zhuǎn)3）。3）計算標(biāo)志向量 。令 ，做如下計算：   &#

46、160;                                                  &

47、#160;                                                   

48、;                             ， ;4）設(shè) 為 的列向量。在 上定義乘法 ，對任意的 ，我們有：a） ；b)如果 ，必有 和 。然后，做如下計算：            

49、                                                   

50、0;                                                  

51、60;                                                  &#

52、160;        ,  ;                                         

53、60;                           ，  6);                     

54、60;               2);5）令 ，并做如下計算：                                 

55、0;                             ， ;令 , 。如果 ，轉(zhuǎn)6）；否則，重復(fù)2）的判斷。6）對 ， ，如果 ，則 。取 的前 列，這樣，我們就得到了一個 維右乘初始輸入矩陣。    這里需要說明的是，運(yùn)用上面的方法求得的右乘初始輸入矩陣不僅與求解雅可比矩陣的列序

56、有關(guān)，而且與過程4）中的合并順序也有關(guān)系。至于如何最優(yōu)求解右乘初始輸入矩陣，目前還很難討論清楚。但是，大量模擬試驗(yàn)結(jié)果表明，運(yùn)用上面自然次序求得的右乘初始輸入矩陣寬度 已經(jīng)非常接近于其下界值 。反向積分 當(dāng) 和 時，通常采用伴隨模式來計算雅可比矩陣。根據(jù)雅可比矩陣的稀疏結(jié)構(gòu)，適當(dāng)選擇左乘初始輸入矩陣，可以獲得接近 的計算時間代價。其中左乘初始輸入矩陣的求解過程完全可以按照上面的方法進(jìn)行，但是在處理前必須先將雅可比矩陣轉(zhuǎn)置，最后還需將得到的初始輸入矩陣轉(zhuǎn)置才能最終得到左乘初始輸入矩陣。同時，其行寬度 也已經(jīng)非常接近于其下界值 ?；旌戏e分  如果將切線性模式和伴隨模式相結(jié)合，往往可以避免梯度向量運(yùn)算中的諸多冗余計算。例如，ADJIFOR系統(tǒng)在求解雅可比矩陣時，在語句級微分實(shí)現(xiàn)中首先用伴隨方法求得所有偏導(dǎo)數(shù)，然后做梯度向量積分；其計算時間代價與 和模式的語句數(shù)目有關(guān)，而其存儲代價為 。具體討論可參考文獻(xiàn)7。5結(jié)論