自动微分转换系统及其应用
摘要 自动微分转换系统（DFT）由LASG和LSEC联合研制开发，目前已拥有成熟的版本。本文对DFT系统的功能、特色及其基本应用作了全面的介绍，并给出了一些颇具说服力的数值试验结果。同时，本文提出了统计准确率评价的概念，这对评价一类自动微分工具及其微分模式代码的可靠性与有效性提供了一种客观的尺度。最后，本文还详细讨论了运用切线性模式求解雅可比矩阵的问题，给出了求解初始输入矩阵的有效算法。
关键词 自动微分 切线性模式 数据相关分析 统计准确率
1.引言
计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析 [4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。
迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7] 和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。
统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。
DFT系统具有许多优点，它能够完全接受用FORTRAN 77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。
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秒）存储开销（字节数）代码复杂性
原模式切线性
模式
原模式切线性
模式
原模式切线性
模式
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系统主要用于求解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