自動微分轉換系統及其應用
摘要 自動微分轉換系統（DFT）由LASG和LSEC聯合研製開發，目前已擁有成熟的版本。本文對DFT系統的功能、特色及其基本應用作了全面的介紹，並給出了一些頗具說服力的數值試驗結果。同時，本文提出了統計準確率評價的概念，這對評價一類自動微分工具及其微分模式程式碼的可靠性與有效性提供了一種客觀的尺度。最後，本文還詳細討論了運用切線性模式求解雅可比矩陣的問題，給出了求解初始輸入矩陣的有效演算法。
關鍵詞 自動微分 切線性模式 資料相關分析 統計準確率
1.引言
計算微分大致經歷了從商微分，符號微分，手寫程式碼到自動微分幾個階段。與其它幾種微分方法相比，自動微分具有程式碼簡練、計算精度高及投入人力少等優點。自動微分實現的基本出發點是：一個數據相對獨立的程式物件（模式、過程、程式段、數值語句乃至數值表示式），無論多麼複雜，總可以分解為一系列有限數目的基本函式（如sin、exp、log）和基本運算操作（加、減、乘、除、乘方）的有序複合；對所有這些基本函式及基本運算操作，重複使用鏈式求導法則，將得到的中間結果自上而下地做正向積分就可以建立起對應的切線性模式，而自下而上地做反向積分就可以建立起對應的伴隨模式[1]。基於自動微分方法得到的切線性模式和伴隨模式，在變分資料同化[2]、系統建模與引數辨識[3]、引數的敏感性分析
[4]、非線性最優化以及數值模式的可預測性分析[5]等問題中有著十分廣泛的應用。
統計準確率被我們視為評價一類自動微分工具及其微分模式程式碼可靠性與有效性的重要尺度。其基本假設是：如果對於定義域空間內隨機抽樣獲得的至多有限個n維初始場（或網格點），微分模式輸出的差分和微分逼近是成功的；那麼對於定義域空間內所有可能初始場（或網格點），微分模式輸出的差分和微分逼近都是成功的。微分模式統計準確率評價的具體方法是：在所有隨機抽樣得到的初始場（或網格點）附近，當輸入擾動逐漸趨向於機器有效精度所能表示的最小正值時，模式輸出的差分和微分之間應該有足夠精度有效位數上的逼近。
2.系統概況
DFT系統主要由兩部分組成：微分程式碼轉換和微分程式碼評價，圖2.1。微分程式碼轉換部分接受使用者輸入指令並自動分析物件模式，生成切線性模式程式碼及其相關測試程式碼，後者直接構成微分程式碼評價系統的主體。微分程式碼評價是DFT系統的一個重要特色。DFT系統的開發小組認為，一個微分模式如果在可靠性、時間和儲存效率上沒有得到充分的驗證，至少對實際應用而言，它將是毫無意義的。
原模式 切線性模式
統計評價結果
圖2.1 DFT系統結構簡圖
2.1 微分程式碼轉換
DFT系統是基於YACC在UNIX環境下開發的，其結構圖2.2所示。通過DFT系統產生的切線性模式程式碼成對出現，並在語句級程度上做了簡化，可讀性很強，如圖2.4。
切線性模式
評價函式集
圖2.2 微分程式碼轉換
微分程式碼轉換部分從功能上分為四個部分：詞法分析，語義分析，物件複雜性及資料相關分析和微分程式碼轉換。對於一組具有複雜資料相關的程式模式物件，通常需要系統執行兩遍才能得到有效而可靠的微分程式碼。這主要有兩方面的考慮：其一，根據物件的複雜性（如最大語句長度、最大變數維數、子過程或函式數目、子過程或函式內最大變數數目等物件特徵）選擇合適的系統引數以求最優的執行代價；其二，模式內各子過程或函式之間以及一個子過程或函式內往往具有很強的資料相關性，需要事先儲存物件的相關資訊並且在考慮當前物件的屬性之前必須做上下文相關分析。
圖2.3 PERIGEE源程式程式碼 圖2.4 DFT系統生成的切線性程式碼
2.2 微分程式碼評價
通常，評價一個編譯系統的效能有很多方面，如處理速度、結果程式碼可靠性及質量、出錯診斷、可擴充套件和可維護性等。對於一類自動微分系統來說，由於軟體開發人力的侷限以及物件模式的複雜多樣性,通過自動轉換得到的微分模式並非常常是有效而可靠的（即無論是在數學意義上還是在程式邏輯上應與期待的理想結果一致），因而在微分模式被投入實際應用前，往往需要投入一定的人力來對其做嚴格的分析測試。
對切線性模式做統計評價測試的主要內容可以簡單敘述為：在網格化的模式定義域空間內，選擇所有可能的網格點形成微分模式計算的初始場；在不同的網格點附近，隨機選取至少個線性無關的初始擾動，對每個擾動輸入分別進行網格點逼近，統計考察模式輸出差分和微分在有效位數上的逼近程度。圖2.5描述了整個測試過程，它包含網格點資料隨機取樣（1）和網格點資料逼近（2）兩級迴圈。
圖2.5 切線性模式程式碼的測試過程
3.系統主要特色
DFT系統並不是一個完整的FORTRAN編譯器，但它幾乎可以接受和處理所有FORTRAN 77編寫的源模式程式碼，並且可以很方便地擴充套件並接受FORTRAN 90編寫的源模式程式碼。本節將著重介紹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.1 GPS Rayshooting模式的相關樹結構片段
3.3 自動生成測試程式
基於IO相關分析的結果，DFT系統自動生成微分測試程式碼，分別對切線性模式的可靠性和執行代價做統計評價測試。特別地，DFT系統還可將任何模式引數都視為輸入輸出引數，生成在數學意義上等價的測試程式碼，這樣處理的不利之處在於往往需要極高的儲存開銷。
3.4 基於語句級的程式碼優化
目前，DFT系統僅僅具備局地優化能力。在語句級微分實現上採用二元歸約的方法對微分程式碼進行優化是DFT系統的一個重要特色。根據右端表示式的乘法複雜性及含變元數目的不同，DFT系統採取不同的分解策略。二元歸約的方法避免了微分計算中的許多冗餘計算，在一些複雜的非線性表示式的微分計算中具有最小的計算代價，同時也非常適合於微分系統的軟體實現。同時，對於某些特殊的運算操作（除法、乘方）和特殊函式（如sqrt、exp），DFT系統較好地利用了基態值計算得到的中間結果，避免了微分實現中的冗餘計算。
4.系統應用
運用自動微分工具得到的切線性模式，可以在無截斷誤差意義下求解函式的數值微分和導數、稀疏雅可比矩陣。同時這些結果在數值引數敏感性分析、非線性最優化以及其它數值理論分析中有著非常重要的應用。這裡簡單介紹切線性模式的幾個基本應用。
4.1 符號導數和微分
如果輸入為數學關係式，DFT系統可以自動生成對應的微分表示式和梯度，而與數學關係式的複雜程度無關。例如我們輸入關係式：
， （1）
DFT系統將自動生成其符號微分形式及其梯度形式分別為
， （2）
4.2 數值導數和微分
切線性模式最基本的應用就是在一定擾動輸入下求解輸出變數的擾動（響應）。表4.1給出了DFT系統在對IAP 9L模式、GPS Rayshooting模式和GPS Raytrace模式三個數值模式做切線性化的具體應用中，一些不同計算粒度、不同引用深度和不同程式風格的核心子過程，以及它們的切線性模式在SGI 2000上執行的統計評價測試結果，其中切線性模式的可靠性指標都準確到六個有效數字以上，在執行時間、儲存開銷和程式碼複雜性方面分別是原模式的兩倍左右，比較接近於理想的微分代價結果（1.5倍）。除了IAP 9L模式由於過於複雜僅做粗略統計外，其餘模式都用非註釋語句行數來表示各自的程式碼複雜性。
表4.1 DFT系統在三個數值模式中的統計評價測試結果
效能指標
物件模式 執行時間（10-3秒） 儲存開銷（位元組數） 程式碼複雜性
原模式 切線性
模式
原模式 切線性
模式
原模式 切線性
模式
Xyz2g 2.530 6.160 5524 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 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
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 13 27
Hlimits 0.830 1.880 242577 484254 37 74
Int3sL 26.90 51.20 820029 1639458 46 90
MAKE
NCEP 1340 3920 722925 1445850 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
×103 3000 5000* 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系統主要用於求解FORTRAN
77語言編寫的切線性模式，具有很強的全域性資料相關分析能力。此外，DFT系統還具有其它幾個重要特色，如結構化的微分實現、自動生成微分測試程式以及基於語句級的微分程式碼優化。本文簡單給出了DFT系統在求解數值和符號導數和微分、稀疏雅可比矩陣中的應用。為評價一類自動微分系統，本文初步提出了統計準確率的概念。
參考文獻
[1] Andreas Griewank. On Automatic Differentiation. In M. Iri and K. Tanabe, editors, Mathematical Programming:
Recent Developments and Applications. Kluwer Academic Publishers, 1989
[2] Le Dimet, F. X and O. Talagrand, Variational algorithms for analysis and assimilation of meteorological
observations: theoretical aspects, Tellus, 1986, 38A, 97-110
[3] P.Werbos, Applications of advances in nonlinear sensitivity analysis, In systems Modeling
and Optimization, New York, 1982, Springer Verlag, 762-777
[4 ] Christian Bischof, Gordon Pusch, and Ralf Knoesel. "Sensitivity Analysis of the MM5 Weather Model using
Automatic Differentiation," Computers in Physics, 0:605-612, 1996
[5] Mu Mu, et al，The predictability problem of weather and climate prediction, Progress in Nature Science, accepted.
[6] Giering R. et al. Recipes for Adjoint Code Construction. ACM Trans. On Math. Software. 1998，24（4）：
437-474.
[7] C. Bischof, A. Carle, P. Khademi, and G. Pusch. "Automatic Differentiation: Obtaining Fast and Reliable
Derivatives--Fast" in Control Problems in Industry, edited by I. Lasiecka and B. Morton, pages 1-16,Birkhauser,
Boston, 1995. CRPC version.
Rostaining N. et al. Automatic Differentiation in Odyssee, Tellus, 1993，45A：558-568