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

關(guān)鍵詞 自動微分 切線性模式 數(shù)據(jù)相關(guān)分析 統(tǒng)計準確率

1.引言

計算微分大致經(jīng)歷了從商微分，符號微分，手寫代碼到自動微分幾個階段。與其它幾種微分方法相比，自動微分具有代碼簡練、計算精度高及投入人力少等優(yōu)點。自動微分實現(xiàn)的基本出發(fā)點是：一個數(shù)據(jù)相對獨立的程序?qū)ο螅Ｊ健⑦^程、程序段、數(shù)值語句乃至數(shù)值表達式），無論多么復雜，總可以分解為一系列有限數(shù)目的基本函數(shù)（如sin、exp、log）和基本運算操作（加、減、乘、除、乘方）的有序復合；對所有這些基本函數(shù)及基本運算操作，重復使用鏈式求導法則，將得到的中間結(jié)果自上而下地做正向積分就可以建立起對應的切線性模式，而自下而上地做反向積分就可以建立起對應的伴隨模式[1]。基于自動微分方法得到的切線性模式和伴隨模式，在變分資料同化[2]、系統(tǒng)建模與參數(shù)辨識[3]、參數(shù)的敏感性分析 [4]、非線性最優(yōu)化以及數(shù)值模式的可預測性分析[5]等問題中有著十分廣泛的應用。
迄今為止，已有數(shù)十所大學和研究所各自開發(fā)了能夠用于求解切線性模式的自動微分系統(tǒng)，比較典型的有TAMC系統(tǒng)[6]、ADJIFOR系統(tǒng)[7] 和ODYSSEE系統(tǒng)[8]。在一些特定的運用中，它們都是比較成功的，但在通用性和復雜問題的處理效率上還存在許多不足。通常，自動生成切線性模式的關(guān)鍵難題在于對象自身的強相關(guān)性，這給系統(tǒng)全局分析（如數(shù)據(jù)IO相關(guān)分析和數(shù)據(jù)依賴相關(guān)分析）和微分代碼的整體優(yōu)化都帶來了很多困難。同時，對于程序?qū)ο蟛豢蓪幍臏蚀_識別和微分處理，至今仍還沒有一個統(tǒng)一而有效的算法。另外，最優(yōu)或有效求解稀疏雅可比矩陣一直是衡量一個自動微分系統(tǒng)有效性的重要尺度。
統(tǒng)計準確率被我們視為評價一類自動微分工具及其微分模式代碼可靠性與有效性的重要尺度。其基本假設(shè)是：如果對于定義域空間內(nèi)隨機抽樣獲得的至多有限個n維初始場（或網(wǎng)格點），微分模式輸出的差分和微分逼近是成功的；那么對于定義域空間內(nèi)所有可能初始場（或網(wǎng)格點），微分模式輸出的差分和微分逼近都是成功的。微分模式統(tǒng)計準確率評價的具體方法是：在所有隨機抽樣得到的初始場（或網(wǎng)格點）附近，當輸入擾動逐漸趨向于機器有效精度所能表示的最小正值時，模式輸出的差分和微分之間應該有足夠精度有效位數(shù)上的逼近。
DFT系統(tǒng)具有許多優(yōu)點，它能夠完全接受用FORTRAN 77語言編寫的源代碼，微分代碼結(jié)構(gòu)清晰，其微分處理能力與問題和對象的規(guī)模及復雜性無關(guān)。它基于YACC實現(xiàn)，具有很強的可擴展性。DFT系統(tǒng)具有四個重要特色。它通過對象全局依賴相關(guān)分析，準確求解雅可比矩陣的稀疏結(jié)構(gòu)，自動計算有效初始輸入矩陣，從而可以用較小的代價求得整個雅可比矩陣。同時，它可以自動生成客觀評價微分模式效率與可靠性的測試程序，對奇異函數(shù)做等價微分處理，并采用二元歸約的方法，在語句級層次上實現(xiàn)微分代碼優(yōu)化。

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

原模式 切線性模式

統(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)，并在語句級程度上做了簡化，可讀性很強，如圖2.4。

切線性模式

評價函數(shù)集

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

圖2.3 PERIGEE源程序代碼 圖2.4 DFT系統(tǒng)生成的切線性代碼

2.2 微分代碼評價

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

圖2.5 切線性模式代碼的測試過程

3.系統(tǒng)主要特色

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

3.1 結(jié)構(gòu)化的微分實現(xiàn)

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

3.2 全局數(shù)據(jù)相關(guān)分析

DFT系統(tǒng)具有較強的數(shù)據(jù)相關(guān)分析能力，它包括全局數(shù)據(jù)IO相關(guān)分析、全局數(shù)據(jù)依賴相關(guān)分析、全局過程相關(guān)分析以及數(shù)據(jù)迭代相關(guān)分析幾個不同方面。數(shù)據(jù)依賴相關(guān)與數(shù)據(jù)IO相關(guān)關(guān)系密切，但又存在根本不同。前者強調(diào)每個變量在數(shù)學關(guān)系上的依賴性；而后者描述了一個對象的輸入輸出特性，且具有相對性，即任何一個變量參數(shù)，無論它是獨立變量還是依賴變量，在數(shù)學意義上都可等價為一個既是輸入又是輸出的參數(shù)來處理。
DFT系統(tǒng)記錄所有過程參數(shù)的IO屬性表，通過深度遞歸相關(guān)計算，準確計算每個過程參數(shù)的最終IO屬性。DFT系統(tǒng)通過對數(shù)據(jù)相關(guān)矩陣做模二和及自乘迭代計算（An+1= An⊕An2）來完成數(shù)據(jù)的依賴相關(guān)分析，這種算法具有很好的對數(shù)收斂特性。DFT系統(tǒng)通過全局過程相關(guān)分析的結(jié)果，自動生成模式的局部或整體相關(guān)引用樹結(jié)構(gòu)（如圖3.1），這對用戶分析復雜數(shù)值模式和微分評價測試都具有很好的指導作用。DFT系統(tǒng)還具有分析局部數(shù)據(jù)迭代相關(guān)和函數(shù)迭代相關(guān)的能力，這兩種形式的數(shù)據(jù)迭代相關(guān)是自動微分實現(xiàn)頗具挑戰(zhàn)的難題之一。

圖3.1 GPS Rayshooting模式的相關(guān)樹結(jié)構(gòu)片段

3.3 自動生成測試程序

基于IO相關(guān)分析的結(jié)果，DFT系統(tǒng)自動生成微分測試代碼，分別對切線性模式的可靠性和運行代價做統(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系統(tǒng)較好地利用了基態(tài)值計算得到的中間結(jié)果，避免了微分實現(xiàn)中的冗余計算。

4.系統(tǒng)應用

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

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

， （1）

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

， （2）

4.2 數(shù)值導數(shù)和微分

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

表4.1 DFT系統(tǒng)在三個數(shù)值模式中的統(tǒng)計評價測試結(jié)果

性能指標
對象模式 運行時間（10-3秒）存儲開銷（字節(jié)數(shù)）代碼復雜性
原模式切線性
模式
原模式切線性
模式
原模式切線性
模式

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*
(含注釋行)
適當設(shè)置輸入擾動的初值，運用切線性模式可以簡單求解輸出變量對輸入的偏導數(shù)。例如，對于一個含有 個輸入?yún)?shù)的實型函數(shù)

（3）

這里設(shè) ， 。運用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ù)。可以通過以下輸入擾動設(shè)置并引用切線性模式（7）來求得偏導數(shù)：

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

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

4.3 稀疏雅可比矩陣

運用上節(jié)討論的方法來求解稀疏雅可比矩陣,具有極高的計算代價。例如，一個含 個獨立和 個依賴參數(shù)的子過程，為求解整個雅可比矩陣就需要反復調(diào)用 次切線性模式，當 相當大時，這對許多實際的數(shù)值計算問題是不能接受的。事實上，如果雅可比矩陣的任意兩列（行）相互正交，那么可以通過適當設(shè)置擾動輸入值，這兩列（行）的元素就可以通過一次引用切線性模式（伴隨模式）完全得到。設(shè) 和 分別為雅可比矩陣的行寬度和列寬度，即各行和各列非零元素數(shù)目的最大值，顯然有 ， 。這里介紹幾種常用的求解方法。
正向積分 當 時，通常采用切線性模式來計算雅可比矩陣。根據(jù)雅可比矩陣的稀疏結(jié)構(gòu)，適當選擇右乘初始輸入矩陣，可以獲得接近的計算時間代價。DFT系統(tǒng)采用一種逐列（行）求解的方法，來有效求解右（左）乘初始輸入矩陣。其基本思路是：按照某種列次序考察雅可比矩陣的各列；考察當前列中所有非零元素，并對這些非零元素所在行的行向量做類似模二和累加運算（即將非零元素視為邏輯“1”，零元素視為邏輯“0”），從而得到一個描述當前列與各行存在“某種”相關(guān)的標志向量（其元素都是“1”或“0”）；依據(jù)此標志向量，就很容易得到一個與之正交的列初始向量 ，其中與當前列序號對應的元素設(shè)置為“1”，而與標志向量中非零元素序號對應的元素設(shè)置為“0”，與標志向量中非零元素序號對應的元素設(shè)置為“-1”，顯然，該列初始向量是唯一的，并且對應著當前右乘初始輸入矩陣的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量與當前求解得到的列初始向量按下面定義的乘法（見過程4）正交，那么這兩列就可以合并，即將當前列初始向量中非“-1”的元素按照對應關(guān)系分別賦值給該初始向量，并從記錄中刪除當前列初始向量；重復以上過程，繼續(xù)按照給定列次序考察雅可比矩陣的“下一列”。不難說明，按照不同列次序求解得到的右乘初始輸入矩陣可能不同。其中逐列求解右乘初始輸入矩陣的過程可以簡單敘述為：
1）將右乘初始輸入矩陣 所有元素的初值均設(shè)置為 ， ， 。 。
2）如果 ，轉(zhuǎn)6）。否則，如果雅可比矩陣 的第 列中的所有元素均為 ， ，重復2）的判斷。否則轉(zhuǎn)3）。
3）計算標志向量 。令 ，做如下計算：
， ;
4）設(shè) 為 的列向量。在 上定義乘法 ，對任意的 ，我們有：a） ；b)如果 ，必有 和 。然后，做如下計算：
, ;
， 6);
2);
5）令 ，并做如下計算：
， ;
令 , 。如果 ，轉(zhuǎn)6）；否則，重復2）的判斷。
6）對 ， ，如果 ，則 。取 的前 列，這樣，我們就得到了一個 維右乘初始輸入矩陣。
這里需要說明的是，運用上面的方法求得的右乘初始輸入矩陣不僅與求解雅可比矩陣的列序有關(guān)，而且與過程4）中的合并順序也有關(guān)系。至于如何最優(yōu)求解右乘初始輸入矩陣，目前還很難討論清楚。但是，大量模擬試驗結(jié)果表明，運用上面自然次序求得的右乘初始輸入矩陣寬度 已經(jīng)非常接近于其下界值 。

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

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

5．結(jié)論

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

[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