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  

上一页  [1] [2]