詮釋自動微分轉換系統(tǒng)及其應用

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

2、及投入人力少等優(yōu)點。自動微分實現(xiàn)的基本出發(fā)點是：一個數(shù)據(jù)相對獨立的程序對象（模式、過程、程序段、數(shù)值語句乃至數(shù)值表達式），無論多么復雜，總可以分解為一系列有限數(shù)目的基本函數(shù)（如sin、exp、log）和基本運算操作（加、減、乘、除、乘方）的有序復合；對所有這些基本函數(shù)及基本運算操作，重復使用鏈式求導法則，將得到的中間結果自上而下地做正向積分就可以建立起對應的切線性模式，而自下而上地做反向積分就可以建立起對應的伴隨模式1?；谧詣游⒎址椒ǖ玫降那芯€性模式和伴隨模式，在變分資料同化2、系統(tǒng)建模與參數(shù)辨識3、參數(shù)的敏感性分析4、非線性最優(yōu)化以及數(shù)值模式的可預測性分析5等問題中有著十分廣泛的應用。 迄

3、今為止，已有數(shù)十所大學和研究所各自開發(fā)了能夠用于求解切線性模式的自動微分系統(tǒng)，比較典型的有TAMC系統(tǒng)6、ADJIFOR系統(tǒng)7和ODYSSEE系統(tǒng)8。在一些特定的運用中，它們都是比較成功的，但在通用性和復雜問題的處理效率上還存在許多不足。通常，自動生成切線性模式的關鍵難題在于對象自身的強相關性，這給系統(tǒng)全局分析（如數(shù)據(jù)IO相關分析和數(shù)據(jù)依賴相關分析）和微分代碼的整體優(yōu)化都帶來了很多困難。同時，對于程序對象不可導處的準確識別和微分處理，至今仍還沒有一個統(tǒng)一而有效的算法。另外，最優(yōu)或有效求解稀疏雅可比矩陣一直是衡量一個自動微分系統(tǒng)有效性的重要尺度。 統(tǒng)計準確率被我們視為評價一類自動微分工具及其微分

4、模式代碼可靠性與有效性的重要尺度。其基本假設是：如果對于定義域空間內(nèi)隨機抽樣獲得的至多有限個n維初始場（或網(wǎng)格點），微分模式輸出的差分和微分逼近是成功的；那么對于定義域空間內(nèi)所有可能初始場（或網(wǎng)格點），微分模式輸出的差分和微分逼近都是成功的。微分模式統(tǒng)計準確率評價的具體方法是：在所有隨機抽樣得到的初始場（或網(wǎng)格點）附近，當輸入擾動逐漸趨向于機器有效精度所能表示的最小正值時，模式輸出的差分和微分之間應該有足夠精度有效位數(shù)上的逼近。 DFT系統(tǒng)具有許多優(yōu)點，它能夠完全接受用FORTRAN77語言編寫的源代碼，微分代碼結構清晰，其微分處理能力與問題和對象的規(guī)模及復雜性無關。它基于YACC實現(xiàn)，具有很

5、強的可擴展性。DFT系統(tǒng)具有四個重要特色。它通過對象全局依賴相關分析，準確求解雅可比矩陣的稀疏結構，自動計算有效初始輸入矩陣，從而可以用較小的代價求得整個雅可比矩陣。同時，它可以自動生成客觀評價微分模式效率與可靠性的測試程序，對奇異函數(shù)做等價微分處理，并采用二元歸約的方法，在語句級層次上實現(xiàn)微分代碼優(yōu)化。 2.系統(tǒng)概況 DFT系統(tǒng)主要由兩部分組成：微分代碼轉換和微分代碼評價，圖2.1。微分代碼轉換部分接受用戶輸入指令并自動分析對象模式，生成切線性模式代碼及其相關測試代碼，后者直接構成微分代碼評價系統(tǒng)的主體。微分代碼評價是DFT系統(tǒng)的一個重要特色。DFT系統(tǒng)的開發(fā)小組認為，一個微分模式如果在可靠

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

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

8、計考察模式輸出差分和微分在有效位數(shù)上的逼近程度。圖2.5描述了整個測試過程，它包含網(wǎng)格點數(shù)據(jù)隨機采樣（1）和網(wǎng)格點數(shù)據(jù)逼近（2）兩級循環(huán)。 3.系統(tǒng)主要特色 DFT系統(tǒng)并不是一個完整的FORTRAN編譯器，但它幾乎可以接受和處理所有FORTRAN77編寫的源模式代碼，并且可以很方便地擴展并接受FORTRAN90編寫的源模式代碼。本節(jié)將著重介紹DFT系統(tǒng)（版本3.0）的以下幾個重要特色。 3.1結構化的微分實現(xiàn) DFT系統(tǒng)采用標準化的代碼實現(xiàn)，切線性模式的擾動變量和基態(tài)值變量、微分計算語句和基態(tài)值計算語句總是成對出現(xiàn)，并具有清晰的程序結構。微分代碼保持了原模式本身的結構和風格（如并行和向量特性、

9、數(shù)據(jù)精度等），即語句到語句、結構到結構的微分實現(xiàn)。在奇異點或不可導處，DFT系統(tǒng)對微分擾動采取簡單的清零處理，實踐證明這對抑制擾動計算溢出具有重要意義，但并不影響評價測試結果。 3.2全局數(shù)據(jù)相關分析 DFT系統(tǒng)具有較強的數(shù)據(jù)相關分析能力，它包括全局數(shù)據(jù)IO相關分析、全局數(shù)據(jù)依賴相關分析、全局過程相關分析以及數(shù)據(jù)迭代相關分析幾個不同方面。數(shù)據(jù)依賴相關與數(shù)據(jù)IO相關關系密切，但又存在根本不同。前者強調(diào)每個變量在數(shù)學關系上的依賴性；而后者描述了一個對象的輸入輸出特性，且具有相對性，即任何一個變量參數(shù)，無論它是獨立變量還是依賴變量，在數(shù)學意義上都可等價為一個既是輸入又是輸出的參數(shù)來處理。 DFT系統(tǒng)

10、記錄所有過程參數(shù)的IO屬性表，通過深度遞歸相關計算，準確計算每個過程參數(shù)的最終IO屬性。DFT系統(tǒng)通過對數(shù)據(jù)相關矩陣做模二和及自乘迭代計算（An+1=AnAn2）來完成數(shù)據(jù)的依賴相關分析，這種算法具有很好的對數(shù)收斂特性。DFT系統(tǒng)通過全局過程相關分析的結果，自動生成模式的局部或整體相關引用樹結構（如圖3.1），這對用戶分析復雜數(shù)值模式和微分評價測試都具有很好的指導作用。DFT系統(tǒng)還具有分析局部數(shù)據(jù)迭代相關和函數(shù)迭代相關的能力，這兩種形式的數(shù)據(jù)迭代相關是自動微分實現(xiàn)頗具挑戰(zhàn)的難題之一。 3.3自動生成測試程序 基于IO相關分析的結果，DFT系統(tǒng)自動生成微分測試代碼，分別對切線性模式的可靠性和運行

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

12、系統(tǒng)較好地利用了基態(tài)值計算得到的中間結果，避免了微分實現(xiàn)中的冗余計算。 4.系統(tǒng)應用 運用自動微分工具得到的切線性模式，可以在無截斷誤差意義下求解函數(shù)的數(shù)值微分和導數(shù)、稀疏雅可比矩陣。同時這些結果在數(shù)值參數(shù)敏感性分析、非線性最優(yōu)化以及其它數(shù)值理論分析中有著非常重要的應用。這里簡單介紹切線性模式的幾個基本應用。 4.1符號導數(shù)和微分 如果輸入為數(shù)學關系式，DFT系統(tǒng)可以自動生成對應的微分表達式和梯度，而與數(shù)學關系式的復雜程度無關。例如我們輸入關系式：，（1） DFT系統(tǒng)將自動生成其符號微分形式及其梯度形式分別為，（2） 4.2數(shù)值導數(shù)和微分 切線性模式最基本的應用就是在一定擾動輸入下求解輸出變量

13、的擾動（響應）。表4.1給出了DFT系統(tǒng)在對IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三個數(shù)值模式做切線性化的具體應用中，一些不同計算粒度、不同引用深度和不同程序風格的核心子過程，以及它們的切線性模式在SGI2000上運行的統(tǒng)計評價測試結果，其中切線性模式的可靠性指標都準確到六個有效數(shù)字以上，在運行時間、存儲開銷和代碼復雜性方面分別是原模式的兩倍左右，比較接近于理想的微分代價結果（1.5倍）。除了IAP9L模式由于過于復雜僅做粗略統(tǒng)計外，其余模式都用非注釋語句行數(shù)來表示各自的代碼復雜性。 適當設置輸入擾動的初值，運用切線性模式可以簡單求解輸出變量對輸入的偏導數(shù)

14、。例如，對于一個含有個輸入?yún)?shù)的實型函數(shù)（3） 這里設，。運用DFT系統(tǒng)，可以得到對應的切線性模式（4） 其中，為切線性模式的擾動輸入?yún)?shù)。可以通過以下辦法來求得偏導數(shù)：（5） 其中。如果對于某個既是輸入?yún)?shù)又是輸出參數(shù)，可以類似以下過程引用的辦法來處理。對于過程引用的情形，例如一個含有個輸入?yún)?shù)的子過程（6） 其中，為輸入?yún)?shù)；，為輸出參數(shù)；，既為輸入?yún)?shù)又為輸出參數(shù)。運用DFT系統(tǒng)，可以得到對應的切線性模式為（7） 其中，分別為切線性模式的微分擾動輸入、輸出和輸入輸出參數(shù)?？梢酝ㄟ^以下輸入擾動設置并引用切線性模式（7）來求得偏導數(shù)：a)設置；（，）；（）可以同時求得（）和（），其中。 b)

15、設置（）；（，）可以同時求得（）和（），其中。 4.3稀疏雅可比矩陣 運用上節(jié)討論的方法來求解稀疏雅可比矩陣,具有極高的計算代價。例如，一個含個獨立和個依賴參數(shù)的子過程，為求解整個雅可比矩陣就需要反復調(diào)用次切線性模式，當相當大時，這對許多實際的數(shù)值計算問題是不能接受的。事實上，如果雅可比矩陣的任意兩列（行）相互正交，那么可以通過適當設置擾動輸入值，這兩列（行）的元素就可以通過一次引用切線性模式（伴隨模式）完全得到。設和分別為雅可比矩陣的行寬度和列寬度，即各行和各列非零元素數(shù)目的最大值，顯然有，。這里介紹幾種常用的求解方法。 正向積分當時，通常采用切線性模式來計算雅可比矩陣。根據(jù)雅可比矩陣的稀疏

16、結構，適當選擇右乘初始輸入矩陣，可以獲得接近的計算時間代價。DFT系統(tǒng)采用一種逐列（行）求解的方法，來有效求解右（左）乘初始輸入矩陣。其基本思路是：按照某種列次序考察雅可比矩陣的各列；考察當前列中所有非零元素，并對這些非零元素所在行的行向量做類似模二和累加運算（即將非零元素視為邏輯“1”，零元素視為邏輯“0”），從而得到一個描述當前列與各行存在“某種”相關的標志向量（其元素都是“1”或“0”）；依據(jù)此標志向量，就很容易得到一個與之正交的列初始向量，其中與當前列序號對應的元素設置為“1”，而與標志向量中非零元素序號對應的元素設置為“0”，與標志向量中非零元素序號對應的元素設置為“-1”，顯然，該

17、列初始向量是唯一的，并且對應著當前右乘初始輸入矩陣的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量與當前求解得到的列初始向量按下面定義的乘法（見過程4）正交，那么這兩列就可以合并，即將當前列初始向量中非“-1”的元素按照對應關系分別賦值給該初始向量，并從記錄中刪除當前列初始向量；重復以上過程，繼續(xù)按照給定列次序考察雅可比矩陣的“下一列”。不難說明，按照不同列次序求解得到的右乘初始輸入矩陣可能不同。其中逐列求解右乘初始輸入矩陣的過程可以簡單敘述為： 1）將右乘初始輸入矩陣所有元素的初值均設置為，。 2）如果，轉6）。否則，如果雅可比矩陣的第列中的所有元素均為，重復2）的判斷。否則轉3

18、）。 3）計算標志向量。令，做如下計算：，; 4）設為的列向量。在上定義乘法，對任意的，我們有：a）；b)如果，必有和。然后，做如下計算：,;，6);2); 5）令，并做如下計算：，;令,。如果，轉6）；否則，重復2）的判斷。 6）對，如果，則。取的前列，這樣，我們就得到了一個維右乘初始輸入矩陣。 這里需要說明的是，運用上面的方法求得的右乘初始輸入矩陣不僅與求解雅可比矩陣的列序有關，而且與過程4）中的合并順序也有關系。至于如何最優(yōu)求解右乘初始輸入矩陣，目前還很難討論清楚。但是，大量模擬試驗結果表明，運用上面自然次序求得的右乘初始輸入矩陣寬度已經(jīng)非常接近于其下界值。 反向積分當和時，通常采用伴隨

19、模式來計算雅可比矩陣。根據(jù)雅可比矩陣的稀疏結構，適當選擇左乘初始輸入矩陣，可以獲得接近的計算時間代價。其中左乘初始輸入矩陣的求解過程完全可以按照上面的方法進行，但是在處理前必須先將雅可比矩陣轉置，最后還需將得到的初始輸入矩陣轉置才能最終得到左乘初始輸入矩陣。同時，其行寬度也已經(jīng)非常接近于其下界值。 混合積分如果將切線性模式和伴隨模式相結合，往往可以避免梯度向量運算中的諸多冗余計算。例如，ADJIFOR系統(tǒng)在求解雅可比矩陣時，在語句級微分實現(xiàn)中首先用伴隨方法求得所有偏導數(shù)，然后做梯度向量積分；其計算時間代價與和模式的語句數(shù)目有關，而其存儲代價為。具體討論可參考文獻7。 5結論 切線性模式在無截斷

20、誤差意義上計算函數(shù)的方向導數(shù)、梯度或雅可比矩陣，以及在模式的可預測性及參數(shù)敏感性分析、伴隨模式構造等相關問題中有著廣泛應用。DFT系統(tǒng)主要用于求解FORTRAN77語言編寫的切線性模式，具有很強的全局數(shù)據(jù)相關分析能力。此外，DFT系統(tǒng)還具有其它幾個重要特色，如結構化的微分實現(xiàn)、自動生成微分測試程序以及基于語句級的微分代碼優(yōu)化。本文簡單給出了DFT系統(tǒng)在求解數(shù)值和符號導數(shù)和微分、稀疏雅可比矩陣中的應用。為評價一類自動微分系統(tǒng)，本文初步提出了統(tǒng)計準確率的概念。 參考文獻 1AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,e

21、ditors,MathematicalProgramming: RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989 2LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological observations:theoreticalaspects,Tellus,1986,38A,97-110 3P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling andOptimization,NewYork,1982,SpringerVerlag,762-777