自動微分轉換系統應用論文

導語：自動微分轉換系統應用論文一文來源于網友上傳，不代表本站觀點，若需要原創文章可咨詢客服老師，歡迎參考。

摘要自動微分轉換系統（DFT）由LASG和LSEC聯合研制開發，目前已擁有成熟的版本。本文對DFT系統的功能、特色及其基本應用作了全面的介紹，并給出了一些頗具說服力的數值試驗結果。同時，本文提出了統計準確率評價的概念，這對評價一類自動微分工具及其微分模式代碼的可靠性與有效性提供了一種客觀的尺度。最后，本文還詳細討論了運用切線性模式求解雅可比矩陣的問題，給出了求解初始輸入矩陣的有效算法。

關鍵詞自動微分切線性模式數據相關分析統計準確率

1.引言

計算微分大致經歷了從商微分，符號微分，手寫代碼到自動微分幾個階段。與其它幾種微分方法相比，自動微分具有代碼簡練、計算精度高及投入人力少等優點。自動微分實現的基本出發點是：一個數據相對獨立的程序對象（模式、過程、程序段、數值語句乃至數值表達式），無論多么復雜，總可以分解為一系列有限數目的基本函數（如sin、exp、log）和基本運算操作（加、減、乘、除、乘方）的有序復合；對所有這些基本函數及基本運算操作，重復使用鏈式求導法則，將得到的中間結果自上而下地做正向積分就可以建立起對應的切線性模式，而自下而上地做反向積分就可以建立起對應的伴隨模式[1]。基于自動微分方法得到的切線性模式和伴隨模式，在變分資料同化[2]、系統建模與參數辨識[3]、參數的敏感性分析[4]、非線性最優化以及數值模式的可預測性分析[5]等問題中有著十分廣泛的應用。

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

統計準確率被我們視為評價一類自動微分工具及其微分模式代碼可靠性與有效性的重要尺度。其基本假設是：如果對于定義域空間內隨機抽樣獲得的至多有限個n維初始場（或網格點），微分模式輸出的差分和微分逼近是成功的；那么對于定義域空間內所有可能初始場（或網格點），微分模式輸出的差分和微分逼近都是成功的。微分模式統計準確率評價的具體方法是：在所有隨機抽樣得到的初始場（或網格點）附近，當輸入擾動逐漸趨向于機器有效精度所能表示的最小正值時，模式輸出的差分和微分之間應該有足夠精度有效位數上的逼近。

DFT系統具有許多優點，它能夠完全接受用FORTRAN77語言編寫的源代碼，微分代碼結構清晰，其微分處理能力與問題和對象的規模及復雜性無關。它基于YACC實現，具有很強的可擴展性。DFT系統具有四個重要特色。它通過對象全局依賴相關分析，準確求解雅可比矩陣的稀疏結構，自動計算有效初始輸入矩陣，從而可以用較小的代價求得整個雅可比矩陣。同時，它可以自動生成客觀評價微分模式效率與可靠性的測試程序，對奇異函數做等價微分處理，并采用二元歸約的方法，在語句級層次上實現微分代碼優化。

2.系統概況

DFT系統主要由兩部分組成：微分代碼轉換和微分代碼評價，圖2.1。微分代碼轉換部分接受用戶輸入指令并自動分析對象模式，生成切線性模式代碼及其相關測試代碼，后者直接構成微分代碼評價系統的主體。微分代碼評價是DFT系統的一個重要特色。DFT系統的開發小組認為，一個微分模式如果在可靠性、時間和存儲效率上沒有得到充分的驗證，至少對實際應用而言，它將是毫無意義的。

原模式切線性模式

統計評價結果

2.1微分代碼轉換

DFT系統是基于YACC在UNIX環境下開發的，其結構圖2.2所示。通過DFT系統產生的切線性模式代碼成對出現，并在語句級程度上做了簡化，可讀性很強，如圖2.4。

切線性模式

評價函數集

圖2.2微分代碼轉換

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

2.2微分代碼評價

通常，評價一個編譯系統的性能有很多方面，如處理速度、結果代碼可靠性及質量、出錯診斷、可擴展和可維護性等。對于一類自動微分系統來說，由于軟件開發人力的局限以及對象模式的復雜多樣性,通過自動轉換得到的微分模式并非常常是有效而可靠的（即無論是在數學意義上還是在程序邏輯上應與期待的理想結果一致），因而在微分模式被投入實際應用前，往往需要投入一定的人力來對其做嚴格的分析測試。

對切線性模式做統計評價測試的主要內容可以簡單敘述為：在網格化的模式定義域空間內，選擇所有可能的網格點形成微分模式計算的初始場；在不同的網格點附近，隨機選取至少個線性無關的初始擾動，對每個擾動輸入分別進行網格點逼近，統計考察模式輸出差分和微分在有效位數上的逼近程度。圖2.5描述了整個測試過程，它包含網格點數據隨機采樣（1）和網格點數據逼近（2）兩級循環。

3.系統主要特色

DFT系統并不是一個完整的FORTRAN編譯器，但它幾乎可以接受和處理所有FORTRAN77編寫的源模式代碼，并且可以很方便地擴展并接受FORTRAN90編寫的源模式代碼。本節將著重介紹DFT系統（版本3.0）的以下幾個重要特色。

3.1結構化的微分實現

DFT系統采用標準化的代碼實現，切線性模式的擾動變量和基態值變量、微分計算語句和基態值計算語句總是成對出現，并具有清晰的程序結構。微分代碼保持了原模式本身的結構和風格（如并行和向量特性、數據精度等），即語句到語句、結構到結構的微分實現。在奇異點或不可導處，DFT系統對微分擾動采取簡單的清零處理，實踐證明這對抑制擾動計算溢出具有重要意義，但并不影響評價測試結果。

3.2全局數據相關分析

DFT系統具有較強的數據相關分析能力，它包括全局數據IO相關分析、全局數據依賴相關分析、全局過程相關分析以及數據迭代相關分析幾個不同方面。數據依賴相關與數據IO相關關系密切，但又存在根本不同。前者強調每個變量在數學關系上的依賴性；而后者描述了一個對象的輸入輸出特性，且具有相對性，即任何一個變量參數，無論它是獨立變量還是依賴變量，在數學意義上都可等價為一個既是輸入又是輸出的參數來處理。

DFT系統記錄所有過程參數的IO屬性表，通過深度遞歸相關計算，準確計算每個過程參數的最終IO屬性。DFT系統通過對數據相關矩陣做模二和及自乘迭代計算（An+1=An⊕An2）來完成數據的依賴相關分析，這種算法具有很好的對數收斂特性。DFT系統通過全局過程相關分析的結果，自動生成模式的局部或整體相關引用樹結構（如圖3.1），這對用戶分析復雜數值模式和微分評價測試都具有很好的指導作用。DFT系統還具有分析局部數據迭代相關和函數迭代相關的能力，這兩種形式的數據迭代相關是自動微分實現頗具挑戰的難題之一。

3.3自動生成測試程序

基于IO相關分析的結果，DFT系統自動生成微分測試代碼，分別對切線性模式的可靠性和運行代價做統計評價測試。特別地，DFT系統還可將任何模式參數都視為輸入輸出參數，生成在數學意義上等價的測試代碼，這樣處理的不利之處在于往往需要極高的存儲開銷。

3.4基于語句級的代碼優化

目前，DFT系統僅僅具備局地優化能力。在語句級微分實現上采用二元歸約的方法對微分代碼進行優化是DFT系統的一個重要特色。根據右端表達式的乘法復雜性及含變元數目的不同，DFT系統采取不同的分解策略。二元歸約的方法避免了微分計算中的許多冗余計算，在一些復雜的非線性表達式的微分計算中具有最小的計算代價，同時也非常適合于微分系統的軟件實現。同時，對于某些特殊的運算操作（除法、乘方）和特殊函數（如sqrt、exp），DFT系統較好地利用了基態值計算得到的中間結果，避免了微分實現中的冗余計算。

4.系統應用

運用自動微分工具得到的切線性模式，可以在無截斷誤差意義下求解函數的數值微分和導數、稀疏雅可比矩陣。同時這些結果在數值參數敏感性分析、非線性最優化以及其它數值理論分析中有著非常重要的應用。這里簡單介紹切線性模式的幾個基本應用。

4.1符號導數和微分

如果輸入為數學關系式，DFT系統可以自動生成對應的微分表達式和梯度，而與數學關系式的復雜程度無關。例如我們輸入關系式：

，（1）

DFT系統將自動生成其符號微分形式及其梯度形式分別為

，（2）

.2數值導數和微分

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

表4.1DFT系統在三個數值模式中的統計評價測試結果

性能指標

對象模式運行時間（10-3秒）存儲開銷（字節數）代碼復雜性

原模式切線性

模式

原模式切線性

模式

原模式切線性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注釋行)2826*

(含注釋行)

適當設置輸入擾動的初值，運用切線性模式可以簡單求解輸出變量對輸入的偏導數。例如，對于一個含有個輸入參數的實型函數

（3）

這里設，。運用DFT系統，可以得到對應的切線性模式

（4）

其中，為切線性模式的擾動輸入參數。可以通過以下辦法來求得偏導數：

（5）

其中。如果對于某個既是輸入參數又是輸出參數，可以類似以下過程引用的辦法來處理。對于過程引用的情形，例如一個含有個輸入參數的子過程

（6）

其中，為輸入參數；，為輸出參數；，既為輸入參數又為輸出參數。運用DFT系統，可以得到對應的切線性模式為

（7）

其中，，，分別為切線性模式的微分擾動輸入、輸出和輸入輸出參數。可以通過以下輸入擾動設置并引用切線性模式（7）來求得偏導數：

a)設置；（，）；（）可以同時求得（）和（），其中。

b)設置（）；；（，）可以同時求得（）和（），其中。

4.3稀疏雅可比矩陣

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

正向積分當時，通常采用切線性模式來計算雅可比矩陣。根據雅可比矩陣的稀疏結構，適當選擇右乘初始輸入矩陣，可以獲得接近的計算時間代價。DFT系統采用一種逐列（行）求解的方法，來有效求解右（左）乘初始輸入矩陣。其基本思路是：按照某種列次序考察雅可比矩陣的各列；考察當前列中所有非零元素，并對這些非零元素所在行的行向量做類似模二和累加運算（即將非零元素視為邏輯“1”，零元素視為邏輯“0”），從而得到一個描述當前列與各行存在“某種”相關的標志向量（其元素都是“1”或“0”）；依據此標志向量，就很容易得到一個與之正交的列初始向量，其中與當前列序號對應的元素設置為“1”，而與標志向量中非零元素序號對應的元素設置為“0”，與標志向量中非零元素序號對應的元素設置為“-1”，顯然，該列初始向量是唯一的，并且對應著當前右乘初始輸入矩陣的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量與當前求解得到的列初始向量按下面定義的乘法（見過程4）正交，那么這兩列就可以合并，即將當前列初始向量中非“-1”的元素按照對應關系分別賦值給該初始向量，并從記錄中刪除當前列初始向量；重復以上過程，繼續按照給定列次序考察雅可比矩陣的“下一列”。不難說明，按照不同列次序求解得到的右乘初始輸入矩陣可能不同。其中逐列求解右乘初始輸入矩陣的過程可以簡單敘述為：

1）將右乘初始輸入矩陣所有元素的初值均設置為，，。。

2）如果，轉6）。否則，如果雅可比矩陣的第列中的所有元素均為，，重復2）的判斷。否則轉3）。

3）計算標志向量。令，做如下計算：

，;

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

,;

，6);

2);

5）令，并做如下計算：

，;

令,。如果，轉6）；否則，重復2）的判斷。

6）對，，如果，則。取的前列，這樣，我們就得到了一個維右乘初始輸入矩陣。

這里需要說明的是，運用上面的方法求得的右乘初始輸入矩陣不僅與求解雅可比矩陣的列序有關，而且與過程4）中的合并順序也有關系。至于如何最優求解右乘初始輸入矩陣，目前還很難討論清楚。但是，大量模擬試驗結果表明，運用上面自然次序求得的右乘初始輸入矩陣寬度已經非常接近于其下界值。

反向積分當和時，通常采用伴隨模式來計算雅可比矩陣。根據雅可比矩陣的稀疏結構，適當選擇左乘初始輸入矩陣，可以獲得接近的計算時間代價。其中左乘初始輸入矩陣的求解過程完全可以按照上面的方法進行，但是在處理前必須先將雅可比矩陣轉置，最后還需將得到的初始輸入矩陣轉置才能最終得到左乘初始輸入矩陣。同時，其行寬度也已經非常接近于其下界值。

混合積分如果將切線性模式和伴隨模式相結合，往往可以避免梯度向量運算中的諸多冗余計算。例如，ADJIFOR系統在求解雅可比矩陣時，在語句級微分實現中首先用伴隨方法求得所有偏導數，然后做梯度向量積分；其計算時間代價與和模式的語句數目有關，而其存儲代價為。具體討論可參考文獻[7]。

5．結論

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

參考文獻

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

Boston,1995.CRPCversion.

RostainingN.etal.AutomaticDifferentiationinOdyssee,Tellus,1993，45A：558-568