（运筹学与控制论专业论文）自动微分广义系统与最优化方法.pdf

摘要 随着科学技术的发展，非线性最优化方法在科学计算和工程分析中起着越来越 重要的作用。而在菲线性最优化的计算中，大多依赖目标函数、约束函数的一阶或 高阶导数及其相关项( 如j a c o b i a n 矩阵与向量的乘积等) 的求解。自动微分是计算这 些导数项的有效工具，与传统的微分方法相比具有计算成本低、计算精度高等优点。 本文从基本系统和广义系统多角度对比和研究了求解一阶导数的自动微分方法，其 中包括逆向模式、正向模式等，并将其应用到增广l a g r a n g e 乘子法中。在求解高阶 导数时讨论了三阶h a l l e y 类方法中的c h e b y s h e v 方法，并对原有算法进行了改进， 最后通过数值试验实现了改进算法并验证了算法的高效性。 关键词自动微分；正向模式；逆向模式；广义系统；增广l a g r a n g e 乘子法；c h e b y s h e v 方法 北京工业大学理学硕士学位论文 a bs t r a c t w i t ht h ed e v e l o p m e n to fs c i e n c ea n dt e c h n o l o g y , n o n l i n e a ro p t i m i z a t i o nm e t h o d p l a y sam o r ea n dm o r es i g n i f i c a n tr o l ei ns c i e n t i f i cc o m p u t a t i o n a n de n g i n e e r i n ga n a l y s i s f i e l d s a sf a ra sn o n l i n e a ro p t i m i z a t i o ni sc o n c e m e d ，t h ec o m p u t a t i o nm o s t l yd e p e n d so n t h ef i r s ta n dh i g ho r d e rd e r i v a t i v e sa sw e l la so t h e ri t e m s ( e g t h em u l t i p l i c a t i o nb e t w e e n j a c o b i a nm a t r i xa n dv e c t o r ) o ft h eo b j e c t i v ef u n c t i o na n dc o n s t r a i n tf u n c t i o n a u t o m a t i c d i f f e r e n t i a t i o n ，w h i c hi sc h a r a c t e r i s t i co fl o wc o s ta n dh i g ha c c u r a c yc o m p a r e dw i t ht h e t r a d i t i o n a ld i f f e r e n t i a t i o nm e t h o d s ，i sa ne f f i c i e n tt o o lt oc a l c u l a t et h ed e r i v a t i v e s i nt h i s p a p e r , f r o md i f f e r e n ts t a n d p o i n t so ft h eb a s i ca n de x t e n d e ds y s t e m ，s o m er e s e a r c hw a s c a r r i e do nt oc o m p a r ea n da n a l y z et h ea u t o m a t i cd i f f e r e n t i a t i o nm e t h o df o rt h ef i r s to r - d e rd e r i v a t i v e ，i n c l u d i n gt h er e v e r s em o d ea n df o r w a r dm o d ea n ds oo n a tt h es a m e t i m e ，t h e s em o d e sw e r ea p p l i e dt ot h ea u g m e n t e dl a g r a n g em u l t i p l i e r w h e na t t e m p t i n gt os o l v et h eh i g ho r d e rd e r i v a t i v e s ，t h i sp a p e rd i s c u s s e d t h ec h e b y s h e vm e t h o dw h i c h b e l o n g st ot h et h i r do r d e rd e r i v a t i v eh a l l e yc l a s sm e t h o d ，a n di m p r o v e dt h eo r i g i n a la l g o r i t h m a tl a s t ，s o m en u m e r i c a le x p e r i m e n t sw e r ed o n et or e a l i z et h ea l g o r i t h ma n d d e m o n s t r a t ei t sh i g he f f i c i e n c y k e y w o r d sa u t o m a t i cd i f f e r e n t i a t i o n ；f o r w a r dm o d e ；r e v e r s em o d e ；e x t e n d e ds y s t e r n ；a u g m e n t e dl a g r a n g em u l t i p l i e r ；c h e b y s h e vm e t h o d i i 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研究成 果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经 发表或撰写过的研究成果，也不包含为获得北京工业大学或其它教育机构的学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 签名：垒! 至1 日期：翌：2 ：笸：曼 关于论文使用授权的说明 本人完全了解北京工业大学有关保留、使用学位论文的规定，即：学校有权保 留送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内 容，可以采用影印、缩印或其他复制手段保存论文 ( 保密的论文在解密后应遵守此规定) 签名：嗵导师签名：雪鬈蝗驺年日期：丛鲨善盥 。、f 第1 章前言 1 1 历史背景 第1 章前言 在科学计算和工程分析中，常常会计算函数一阶乃至高阶导数，如函数 的梯度、j a c o b i a n 矩阵、h e s s i a n 矩阵以及三阶张量目前，主要有四种不同的 计算导数的方法：手写代码、有限差分、符号微分以及自动微分( 也称算法微 分，a u t o m a t i cd i f f e r e n t i a t i o n ) 。比起传统意义上的微分方法，自动微分( a d ) 具 有代码简练、计算精度高及投入人力少等优点如果我们直接编写导数代码， 不仅费时，而且极易出错，当问题的规模非常巨大时，单靠人力几乎是不可 能完成的有限差分可以用差商来近似导数，但求解函数的梯度是相对于计 算函数本身具有极高的计算复杂性( 其计算代价依赖于独立变量的数目) ，而 且由于受机器有效精度的限制不可避免地存在截断误差。符号微分一般基于 一些成熟的数学通用软件，如m a t h e m a t i c a ，m a p l e ，r e d u c e 等，准确求得复杂表 达式或函数的导数，但对具有数据相关性的复杂对象模式来说，则具有很大 的局限性基于原模式代码，自动微分依籁机器自动生成数值函数的导数或 微分代码，通过整体分析对象模式的结构特征和数据相关性，可以得到一今 合理的与问题本身规模无关的计算时间和空间存储代价。同时，自动微分针 对程序代码分析求解数值函数的导数，因而它的一个突出优点是无截断误差， 可以准确到机器精度。另外，基于自动微分得到的微分模式代码在变分资料 同化、最优化方法、奇异向量与奇异值问题、系统参数识别，多学科综合设计 以及敏感性分析等领域中有着非常广泛和成功的应用 3 5 3 6 1 1 3 7 。 自动微分方法最早是由w e n g e r tr e ( 1 9 6 4 ) 、b e d al m ( 1 9 5 9 ) 等人正式提 出，不过那时人们一般普遍认为分析求解函数的导数具有极高的计算代价。 随后，l i n n a i n m as ( 1 9 7 9 ) 、w a e n e rd d ( 1 9 7 5 ) 、k e d e mg ( 1 9 8 0 ) 、b a u rw ( 1 9 8 2 ) 等人 先后在求解j a c o b i a n 矩阵的系数结构、正向和逆向模式的最小代价实现、导数 求解的计算复杂性以及求解h e s s i a n 矩阵等方面做过许多重要的工作直到八 十年代，随着计算机技术与软件技术的飞速发展及自动微分基本理论和方法 北京工业大学理学硕士学位论文 的逐渐成熟，r a l ll b ( 1 9 8 1 ) 、g r i e w a n ka ( 1 9 8 9 ) 、b i s c h o f c ( 1 9 9 1 ) 和c h r i s t i a n s o n b ( 2 0 0 0 ) 等人在a d 基本方法、高阶导数分析求解，特别是在揭示计算导数的 真实计算代价方面做了系统的工作，并得到了极为广泛的应用。在这方面比较 有代表性的文献有 6 7 8 3 1 9 。u w en a u m a n n ( 1 9 9 9 ) 在最近完成的博士论文 1 0 】 中采用计算图分析方法，系统地讨论了基于非增量积分法的导数计算，并提 出了基于节点消去法来寻找最优消去次序等价于一个n p 难解问题。g r i e w a n k a ( 2 0 0 0 ) 在最近出版的著作e v a l u a t i n gd e r i c a t i c e s 一书中【11 ，对目前a d 研究 的基本现状以及a d 基本理论框架结构作了一个非常好的概括和论述，为今 后a d 的发展及应用提供了非常有价值的思路和文献资料，这也标志着a d 已 经成为一门独立的学科分支，其主要参考文献见 2 3 1 3 1 4 。 通常，用来计算函数的导数有两种不同的基本模式，即正向模式( f o r w a r d m o d e ) 和逆向模式( r e v e r s em o d e ) 正向模式就是按照程序运行的自然顺序自 上而下的计算函数的梯度，而逆向模式则是按照与程序运行相反的顺序自下 而上的计算函数的梯度作为计算科学领域的一个独立分支，目前自动微分 研究包括；稀疏j a c o b i a n 矩阵及稀疏h e s s i a n 矩阵，高阶导数求解，奇异微分， 并行计算微分，隐式方程求解，微分代码的评价和优化等系列问题 1 】 2 。 自动微分在实际的计算中有着极其广泛的应用，如数学物理中的许多数 值模拟问题都与非线性最优化方法有关，或者最终可以归结为一个非线性最 优化问题来求解。运用自动微分方法得到的微分模式代码，统称以合理的代 价求解目标函数的梯度和h e s s i a n 矩阵，以及向量值函数的j a c o b i a n 矩阵和各 种形式的方向导数这些问题包括：系统参数辨别、变分资料同化、多学科最 优设计、模式的可预测性分析等，关于这方面的文献如 3 】 4 5 等。 最优化理论与算法是一个重要的数学分支，它所研究的问题是讨论在众 多方案中什么样的方案最优以及怎样找出最优方案长期以来，人们对最优 化问题一直进行着不懈地探讨和研究早在1 7 世纪，英国科学家n e w t o n 发明 微积分的时代，就已提出极值问题，后来又出现l a g r a n g e 乘子法。1 8 4 7 年法 国数学家c a u c h y 研究了函数值沿什么方向下降最快的问题，提出了最速下降 2 第1 章前言 法。1 9 3 9 年前苏联数学家j b k a h t o p o b h y 提出了解决下料问题和运输问题这两 种线性规划问题的求解方法关于最优化问题的研究工作，随着历史的发展 不断深入。2 0 世纪4 0 年代以来，由于生产和科学技术突飞猛进的发展，特别 是电子计算机应用日益广泛，使最优化问题的研究不仅成为一种迫切需要， 而且有了有力的求解工具。与此同时最优化理论和算法迅速发展起来，形成 了一个新的学科 在众多最优化方法中，非线性最优化方法在科学计算和工程分析中起着 越来越重要的作用，带有等式约束的非线性最优化方法是科研人员和工程技 术人员所关注的热点。通过求解一系列无约束最优化问题来得到约束问题的 最优解，称这类方法为序列无约束极小化方法本文所涉及到的增广l a g r a n g e 乘子法就属于这种方法，另外还讨论了求解无约束c h e b y s h e v 算法的改进算 法 1 2本文工作概述 本文第二章主要分析了自动微分方法原理及其实现。第三章结合广义系 统相关结论，研究从广义系统中得到j a c o b i a n 矩阵的方法第四章主要研究分 析了将自动微分应用到最优化理论中的增广l a g r a n g e 乘子法和三阶c h e b y s h e v 方法并对其中的算法进行改进，通过数值试验验证了算法的可行性与优越性。 最后本文作了总结和展望。 3 第2 章自动微分算法及其实现 第2 章自动微分算法及其实现 自动微分方法讨论如何以合理的代价求解数值函数的梯度、j a c o b i a n 矩阵、 更高阶导数和微分以及多种矩阵一向量乘积形式的微分模式其基本出发点 是：任何程序意义的数值模式，无论多么复杂，总可以分解为一系列有序的 有限数目的基本函数和基本运算操作的运算复合。按照链式求导法则，可以 形式地对程序代码做微分线性化处理，从而达到对整个模式对象做微分处理 的目的 基于不同的数值需要，自动微分方法依赖机器自动形式化地对程序模式 代码做一阶或高阶微分线性化处理其处理的对象通常是一组具有数据相关 性的数值程序模式或单个独立子过程或函数，也可以是单个的数值表达式、 程序语句或程序段。自动微分主要基于这样的思路实现：一个相对独立的程 序对象、数值模式、独立子过程或函数，按照其输入输出关系可被视为一个 分段可微的算子，而该算子又可以按不同粒度、过程或函数划分为一系列最 小算子的有序复合重复使用链式求导法则，可以形式地对程序对象做一阶 或高阶微分线性化处理 2 1 自动微分( a d ) 的两种基本模式 自动微分包括正向和逆向两种基本模式，对于一个给定数值模式，如果 从独立变量到中间变量按照程序执行的自然路径求导，就得到正向模式；相 反地，如果从依赖变量到中间变量按照与程序执行相反的计算路径反向求导， 就得到逆向模式。通过正向模式和逆向模式来计算j a c o b i a n 矩阵具有不同的 计算代价，当依赖变量的数目小于独立变量的数目时，运用逆向模式往往只 需较小的计算时间代价；而当独立变量的数目小于依赖变量的数目时，运用 正向模式往往只需较小的计算时间和空间存储代价。 下面以求解j a c o b i a n 矩阵为例来具体说明正向模式和逆向模式的实现形 式 北京工业大学理学硕士学位论文 对于如下的数值模式y = r ( x ) ，其中y r ”为函数变量集，x r “为 自变量集。除有限点外，f 具有一阶连续导数 求解原函数f ( z ) 的顺序执行的w e n g e r t 表为： 砺一。= x i ， l 珏 v i = 咿j ( ) j i 0 其中f 为中间变量的个数。 求解j a c o b i a n 矩阵的正向模式的基本结构为： * h i n = 蕊，i = 1 ，住 啦筹江n “n + 2 ，n + f j ，1 5 j 曼 j 雪m t = 心z i ，i = m 一1 ，0 求解j a c o b i a n 矩阵的逆向模式的基本结构为： 玩= 0 ，i = 1 一n ：，z 奶一t = 哥。一 ，i = 0 ，m 一1 吗+ = 砚+ 型o v j ，歹 i ，i = f ，l 牙t = 移i n ，i = n ，1 ( 2 1 ) ( 2 2 ) 2 3 ) ( 2 4 ) ( 2 - 5 ) ( 2 6 ) ( 2 7 ) ( 2 - 8 ) ( 2 9 ) ( 2 1 0 ) 对于具体的问题，就是把函数划分为一些基本的程序段，运用中间变量把 每个程序段联系起来，然后运用正向模式或逆向模式的链式法则求出j a c o b i a n 矩阵。相对于计算函数本身的代价而言运用正向模式求解雅可比矩阵，其计 算量不超过0 ( n ) ；而运用逆向模式来求解，其计算量不超过d ( m ) ( 本文所有的 计算量均指乘除法的次数) 3 8 。下面用具体实例来说明自动微分算法的实现 过程。 第2 章自动微分算法及其实现 2 2 自动微分的实现 考虑向量值函数f ，z 。，z ：，z 。为基本变量，求解v f = 箬，篷，筹】_ t ，其 中f = e x l 。2 + z l z 2s i n x 3 引入中间变量v i , i = l ，2 ，3 ，4 ，作如下代换： 1 ) 利用正向模式求解： 由链式法则整理得： u l2x l 丰3 ：2 y 22s i nx 3 v 32e ”1 0 42v l 求? 2 2 f = 7 3 3 + u 4 d v l2d z l 术x 2 + z 1 术d x 2 d v 2 = c o s x 3 木d x 3 d v 32y 3 术d r l d v 4 。d v l 木v 2 十v l 木d 2 d f = d v 3 牟d v 4 d f = 向2 e 。1 2 2 + z 2s i n ：r z ) d z l + ( z l e 。1 。2 + z 1s i nx a ) d z 2 + x l x 2 c o sx 3 d x 3 分别令 ( d z l ，d z 2 ，d x 3 ) t ：( 1 ，o ，o ) 丁，( d x l ，d x 2 ，d x 3 ) t = ( o ，1 ，o ) 丁，( d x l ，d x 2 ，d x 3 ) t = ( o ，0 ，1 ) 下 从而得到： v f = 署，鸶) 丁= 3 z i x 2s i n 2 2 x 3 z 。s i n x 2 c o s z 2 + z s i n 2 x 2 - t 7 。 北京工业大学理学硕士学位论文 2 ) 利用逆向模式求解： d x l = 0 ，d x 2 = 0 ，d x a = 0 ，d v i = 0 ，i = 1 ，2 ，4 ，5 d u 3 + = d f , d 口4 + = d f d 钉1 十= v 2 d v 4 ，d 移2 + = v l d v 4 d v l + = v 3 d v 3 d z 3 + = c o s x a d v 2 d z 2 + = z l d v l ，d x l + = x 2 d v l 由链式法则整理得： v f = 石d f ，荔d f ，坚d x 3 ) 丁= x 2 e z z 2 + x 2 s i n x 3 , x 1 e z l x 2 j _ x 1s i n x 3 , x l x 2c o s x a 7 由此可见，正向模式结果与逆向模式结果相同。下面我们给出具体的自 动微分算法 2 3自动微分算法 考虑函数厂( z ) ：眇_ r ，其梯度和h e s s i a n 矩阵分别用v f ( x ) ，v 2 f ( x ) 表示， 种子向量用圣( 逆向模式的输入值) 表示 给定初始点z 和种子向量圣，下列步骤可以得到，( z ) ：v f ( x ) ，v 2 ，( z ) 圣。 第一步：使用链式规则得到，( z ) ； 第二步：利用逆向规则计算函数梯度v 歹( 。) ； 第三步：令f ( z ) = ( ，( z ) ，v ，p ) ) ：r ”_ 咒叶1 ，对f ( z ) 再次使用正向规则 便得到v 2 ，( z ) 士 说明：可以在第二步与第三步之间设置两个开关，当仅需要v f ( x ) 时，输 出v f ( x ) 而不继续第三步；当需计算v 2 f ( x ) 时，可以分别取圣= e ，e 2 ， 并输出v 2 ，( z ) ，其中e l ，e 2 ，为单位矩阵厶的行向量 8 第2 章自动微分算法及其实现 2 4 自动微分两种模式的计算量 以上详细介绍了自动微分的两种基本模式，在这一节中以计算函数梯度 为例，分别讨论两种模式的计算量首先看一下用正向模式计算函数梯度所 需的计算量。 一、正向模式自动微分的计算量 首先给出几个定义： 定义2 4 1 ： ( 1 ) w o n k t a s k ( f ) 表示计算与函数f 导数相关量的复杂性 ( 2 ) w o r k e v a l ( f ) 表示计算函数f 的复杂性 假设2 4 1 ； ， f w o r k t a s 南( f ) ) w o n n ( t a s k f ) l ( 2 - 11 ) i = 1 w o r k e v a l ( p ) = ：w o r k e v a l ( q o i ) ( 2 1 2 ) i = l 其中t a s k 7 ( 妒) 不同于t a s k ( 妒) ，因为t a s k ( ) 仅表示函数妒中的计算量，而没有 包含这些元素的复合计算量，所以这里用t a s k 7 ( 妒) 来表示妒的总计算量 在自动微分中，我们规定把复杂函数分解为若干初等函数，因此w o r n t a s k 中只包含四种初等运算，即位移、加法、乘法、非线性操作因为除法操作是 乘法的逆运算，幂运算和三角运算都为非线性操作，这里认为它们的计算量 相当现举一个例子来说明，例如乘法算子：u 兰妒( u ，u ) = u 车u ，这个运算中 包括3 个存取运算和1 个乘法运算，于是有w o r k e v a l ( q o ) = ( 3 ，o ，1 ，o ) r ，同 样当矽( _ l z ) = s i n ( u ) ，此时d r k e 懈f ( 砂) ) = ( 2 ，0 ，0 ，1 ) 7 表2 一l 给出这四种运 算各包含的操作数： 北京工业大学理学硕士学位论文 表2 1 四种运算的操作数 t a b l e2 - 1t h eo p e r a n d so ft h e s ef o u ro p e r a t i o n s e u n fc土 木 砂 位移 l332 加法 0 l00 乘法 0010 非线性操作 00o1 根据假设2 4 1 ，有： o r k e u n f ( f ) ) 三眠。l f l 三 1332 0100 0010 00o1 ( 2 1 3 ) 其中i f l = ( 1 1 ，1 2 ，f 3 ，1 4 ) 丁称为f 的操作频数向量我们用矩阵m 。与矩阵1 哌删 的比值g 。七来表示f 的计算复杂度，即 c t a s k 篇 p 面斫面面丽 毕。1 叫 由假设2 4 1 还可得出 z t i m e 伽f ( f ) = t i m e e u a f ( 妒i ) ) ( 2 1 5 ) i = 1 f t i m e t a s 七( f ) ) t i m e t a s 后 ) ) ( 2 1 6 ) i = 1 ，l o 第2 章自动微分算法及其实现 于是可用下式将函数计算时间与函数计算量联系起来： t i m e t a s k ( f ) = 。t 彬0 冗k t q s ( f ) ) ( 2 1 7 ) 其中是d i m ( w o r k ) 的权重向量，设 u t = ( p ，1 ，丌，) ，p2m a x 1 ，7 r 2 ) ，7 r 1 ，2 z r ( 2 1 8 ) 也就是说一个非线性操作至少是乘法计算量的两倍，而存取操作至少相当于 乘法计算量或为乘法计算量的一半 有了上述的定义，现在来分析一下正向模式自动微分计算函数梯度所需 的计算量由正向模式可知，计算函数梯度所需的计算量如表2 - 2 所示： 表2 2 正向模式计算函数梯度计算量 t a b l e2 - 2t h eo p e r a t i o nf o rc o m p u t i n gg r a d su s i n gf o r w a r dm o d e t a n g c土 木 移 位移 l + l 3 + 3 3 十32 + 2 加法 ol + 10 + 1o + 0 乘法 o0l + 2o + 1 非线性操作 0001 + 1 表2 2 中加号左边的数字表示计算函数值所需各种操作的计算量，右边 的数字表示运用正向模式求阶导数所需的计算量。例如：，( z ) = z - 乖z 2 ，计 算函数本身需要3 个位移操作，1 个乘法操作现在对f ( x ) 求一阶导数，有 ，7 ( z ) = 2 1 丰x 2 + 窖2 木z l ，这个过程需要3 个位移操作，1 个加法操作，2 个乘法 操作，其它操作亦按照这种方法计算 北京工业大学理学硕士学位论文 此外从表2 2 中得到正向模式计算量矩阵为： 彤8 n 9 三 w ：。l 三 基于假设( 2 1 8 ) 和( 2 1 4 ) 式，可以得到正向模式计算量的估计值： = m a x 警，筹， 6 “+ 1 + 3 7 r4 3 “十7 r “+ 7 r + 2 工， 2 “+ ( 2 1 9 ) ) 2 ：5 2 ( 2 2 0 ) 可见利用正向模式自动微分计算函数梯度的计算量至多是函数计算量的2 5 倍。 其次再分析一下用逆向模式自动微分计算函数梯度所需的计算量。 二、逆向模式自动微分的计算量 由逆向模式计算函数梯度所需的计算量如表2 3 所示。 表2 3 逆向模式计算函数梯度计算量 t a b l e2 3t h eo p e r a t i o nf o rc o m p u t i n gt h eg r a d su s i n gr e v e r s em o d e g r a d c士 毒 移 位移 l + 13 + 63 + 82 + 5 加法 0l + 2o + 20 + 1 乘法 o 01 + 2 0 + 1 非线性操作 0oo1 + l 表2 3 中加号左边的数字表示计算函数值所需各种操作的计算量，右边 的数字表示运用逆向模式求一阶导数所需的计算量例如：妒( z ) = e x p ( z ) ，计 1 2 2 o o l 3 0 l q 3 1 0 o 1 o 0 0 4 o 1 2 6 1 3 0 6 2 0 o 2 o 0 0 第2 章自动微分算法及其实现 算函数本身需要2 个位移操作，1 个非线性操作现在对v ( x ) 求一阶导数，有 妒面) + = e x p ( z ) 车牙，这个过程需要5 个位移操作，1 个加法操作，1 个乘法操 作，1 个非线性操作，其它操作也是按照这种方法计算于是，从表2 3 中得 到逆向模式计算量矩阵： 。d 三 291 17 o321 0031 0002 w ：z 兰 13 32 o1 0o o 010 000 1 基于上述假设( 2 1 8 ) ，便可以得到逆向模式计算量的估计值： ( 2 2 1 ) 啪= 酬罟，鬻，紫，等等竽) 3 州( z - z 2 ) u 舭d2m a x t 了，石玎石石_ ，1 再丁一 h4 j 可见利用逆向模式计算函数梯度的计算量大概是函数计算量的3 至4 倍 2 5 本章小结 本章主要介绍了自动微分的两种基本形式：正向模式和逆向模式及两种 基本模式的执行过程和计算量。文中用具体的实例详细说明了自动微分的应 用过程，并且给出算法的执行步骤。在计算量的讨论中，主要得出以下结论： 利用正向模式计算函数梯度计算量至多是函数计算量的2 5 倍，而逆向模式计 算函数梯度计算量至多是函数计算量的4 倍这一结果很好地说明了不管采 用自动微分正向模式还是逆向模式计算函数梯度，都具有计算量小的优点， 为多种情况的问题求解作了概括和总结。 1 3 第3 章基于自动微分的广义系统 第3 章基于自动微分的广义系统 增广l a g r a n g e 乘子法在执行过程中需要用到增广l a g r a n g e 函数的梯度、j a c o b i a n 矩阵的信息，怎样才能提高增广l a g r a n g e 函数梯度、j a c o b i a n 矩阵的计算效率 成为本章研究的重点问题本章介绍和讨论广义系统，并且从广义系统中得 到j a c o b i a n 矩阵，从而提高计算j a c o b i a n 矩阵的效率 广义系统e ( z ，) = 0 这里将第2 章中求解原函数f ( z ) 顺序执行的w e n g e r t 表压缩成非线性等 式系统 1 1 ： 0 = e ( x ，l , ) 兰( i ( u i ) 一a ) i ；1 飞，f( 3 一1 ) 其中e 的前礼个分量定义为初始化函数0 个输入变量) ： ( u t ) 三x i + 。，i = 1 一祝，0 ( 3 - 2 ) m 个相互独立的输出变量： y m - = i ? i = m l ，0( 3 - 3 ) 定义 = 筹( 3 - 4 ， 2 瓦 j 当i f m 有c i j 三0 ，则 三( 乱t ) 三o j “_ 2 1 一n i ，j f( 3 - 5 ) e 的j a c o b i a n 矩阵有n4 - 1 个变量吩，j = 1 一n ，l ，且该j a c o b i a n 矩阵为下 三角矩阵： u ( x ，) = ( c 妇一文j ) 蒿二：0 = c 一，= 1 5 北京工业大学理学硕士学位论文 。 1 l 一_ ，0l ( 3 - 6 ) i ti1 由初始化函数定义，当i k 2 0 ( 3 - 2 4 ) ( 3 2 5 ) 第3 章基于自动微分的广义系统 v l v 2 可l y 2 x lx 2x 3 v l忱u 3魄 0 - - 1000 00 0 - - 100 00 0o 210 10 00011 00 000 oo 0 o003 第一步：消去地所在的行、列 2 1 7 1x , 2z 31 1 忱v 3 魄 0 10 00 0 一1 0 0 10 002 1 00 00 000oo 3 2 1 ( 3 2 6 ) ( 3 2 7 ) 魄 收 i | 、 一l 丁 b r ， n 沈 地 玑 驰 北京工业大学理学硕士学位论文 第二步：消去1 2 所在的行、列 扩1 1 2 v 3 v 4 y l y 2 z l2 7 2x 3v l 忱怕蚝 0 - - 1 q l0 o一1 00 0 00 00 3 第三步：消去1 所在的行、列 1 1 垤 v 4 y t x lx 22 7 3 1 屹1 3 吮 木0 1 00 000 3 第四步：消去地所在的行、列 2 2 ( 3 2 8 ) ( 3 2 9 ) 第3 章基于自动微分的广义系统 z 1z 2x 31 11 2 的v 4 o ( 3 3 0 ) 此时便得到了所求的原函数f ( x ) 的j a c o b i a n 矩阵，即上述过程最后一个 矩阵的左下角矩阵。现在将两组数字作一下比较： 表3 1 消去法与自动微分计算量比较 t a b l e3 - 1t h eo p e r a t i o nc o m p a r i s o nb e t w e e ne l i m i n a t i o nm e t h o da n da u t o m a t i o n d i f f e r e n t i a t i o n 求j a c o b i a n 矩阵方法加法数量乘法数量 自动微分法( a d ) 9 2 1 消去法 17 从表3 1 中可以看出，消去法显然比正向模式自动微分计算j a c o b i a n 矩阵 所需的计算量小由图论知识可知，顶点消去法和边消去法的本质是变换图 形，使其相互的对应关系保持不变。从而保证j a c o b i a n 矩阵不变因为这种操 作是将线性方程组以某种线性变换进行转换，而没有改变输入变量与输出变 2 3 n 忱 地 魄 玑 驰 北京工业大学理学硕士学位论文 量之间的关系( 证明见引理3 3 ) 由广义系统( 3 6 ) 知 ，、 l 000 l il u2i 口l0 i ll r t o ( 3 - 3 1 ) g = ( 三三) c 3 3 2 ， 引理3 3 阵汐、y ：= ( 二；) ，y = ( 二三) ，其中w 和z 都是非奇异阵。有 2 4 第3 章基于自动微分的广义系统 所给矩阵不是稀疏矩阵，而是满的矩阵，那么两种方法所用到的计算量是一 样的但在实际应用中，这个矩阵往往都是稀疏形式的，所以用消去法要比 用自动微分方法节省计算量。 综上可知j a c o b i a n 矩阵有以下两种计算形式： f 7 ( z ) = r t ( l 一，) 一1 b = 兄一t ( 一，) 1 b 】= r 一 t ( l 一，) 一1 b( 3 3 3 ) 现在若要计算j a c o b i a n 矩阵乘向量的形式，即f 7 ( z ) 圣，可以用下面两种思路 来实现首先令圣为输入向量： f 7 ( z ) = r 圣一t ( l 一，) 一1 b j c 此式中涉及求矩阵的逆，避免求逆计算，令： ( l t ) z = b j c 因为l j 是严格下三角矩阵，所以z 很容易求，此时： ( 3 3 4 ) ( 3 3 5 ) f 7 ( z ) 圣= r d ：一t ( l 一，) 一1 b 窑= j r 孟一丁2 ( 3 3 6 ) 其次还可以令雪为输入向量： 雪丁f 7 ( z ) = 雪t r 圣一9 t t ( l 一，) 一1 b ( 3 3 7 ) 同样，为避免求逆计算，令： 护( l 一7 ) = 雪丁t ( 3 3 8 ) 即( l 一，) 丁2 = t 丁妒 因为( l ，) t 是严格上三角矩阵，所以2 很容易求，此时： 雪t f 7 ( z ) = 9 t r 一9 r t ( l j ) 一1 b = 哥t r 一2 r b ( 3 3 9 ) 2 5 北京工业大学理学硕士学位论文 3 3本章小结 本章主要对广义系统作了详细的阐述。广义系统的j a c o b i a n 矩阵和原函数 f ( z ) 的j a c o b i a n 矩阵有着紧密的联系。通过引入广义系统，本章分析了一些 求解原函数j a c o b i a n 矩阵的有效方法，并且通过理论和例子分别说明了消去 法和自动微分方法的优缺点 2 6 第4 章基于自动微分的最优化方法 第4 章基于自动微分的最优化方法 由第2 章可知自动微分在计算函数梯度、j a c o b i a n 矩阵以及高阶导数方面 有其独特的优点，本章把自动微分技术应用到最优化算法中，以此来提高原 有算法的效率 4 1应用自动微分的约束最优化方法 罚函数法最早是由c o u r a n t 在1 9 4 3 年提出，其基本思想是：对不满足约束 条件的点进行惩罚，通过求解多个罚函数的极小值得到约束问题的最优解 考虑等式约束问题 m 砌i n ，( z ) ( 4 - i ) s t c ( 2 7 、= 0 其中f ：r “一r ，c ：r “_ 胛都是平滑的非凸函数。现将问题( 4 1 ) 转化为无 约束问题，即希望构造相应的无约束问题，而该无约束问题的解恰好是约束 问题的最优解，要做到这一点就必须增大在非可行域处的目标函数值，即对 非可行域处的目标函数加以惩罚由此( 4 1 ) 的罚函数可定义为： p ( z ，仃) = f ( x ) + i o c t ( z ) c ( z ) ( 4 2 ) 厶 可以期望，当盯充分大时，无约束问题m i ne ( x ，盯) = f ( z ) + 暑c t ( 2 7 ) c ( 2 7 ) 的最 优解牙= 牙( 仃) 接近约束问题的最优解z 罚函数的优点是将约束问题化成一系列无约束问题，可用求解无约束问 题的算法得到约束问题的最优解。但是当惩罚因子吼充分大后，函数e ( 2 7 ，o k ) 通常是一个病态函数，即其h e s s i a n 矩阵v ：e ( 2 7 ( 扪，o r k ) 的条件数非常大，这给 无约束问题的求解带来了困难为克服这一缺点，h e s t e n e s 和p o w e l l 于1 9 6 4 年 各自独立地提出了乘子法。对于问题( 4 1 ) 设2 7 + 为其全局解，考察矿是否为 2 7 北京工业大学理学硕士学位论文 无约束问题m i np ( z ，盯) = ，( z ) + 号c t ( z ) c ( z ) 的全局解。p ( z ，盯) 在3 7 4 处的梯度为 v 。p ( x ，盯) = v ( x + ) + 口vc ( 矿) t c ( z + ) = v f ( z + ) 。通常意义下，v f ( z + ) 0 ，也 就是说矿不是m i np ( z ，盯) = 厂( z ) + 号c 7 ( z ) c ( z ) 的全局解，因此就不能期望选择 一个固定的口求解一个无约束问题，得到约束问题的最优解 现在就提出一个问题：是否存在这样的函数( z ，o r ) 恰好是无约束问题 r n i n p ，g r ) 的全局解或局部解，以期望通过求解一个无约束问题而得到约束 问题的最优解，这样可大大提高算法效率。于是对l a g r a n g e 函数加以惩罚， 这样得到函数( z ，盯) = ，( z ) + a + c ( z ) + 量c t ( z ) c ( z ) ，可以期望，当仃充分大 时，有v ( ，盯) 正定。当构造出函数( z ，伊) 后，可以通过求解一个无约束 问题得到约束问题的最优解但事实上，并不能做到这一点。因为( z ，矿) 中 的砖未知，克服上述困难的方法是用参数入代替 ，得到增广l a g r a n g e 函数 a ，盯) = f ( z ) + a c ( z ) + 署c t ( z ) c ( z ) 考虑相应的无约束问题r a i n ( z ，a ，盯) ，其 最优解为主= 孟( a ，口) 。这样做可以保证只要当盯大到一定量后，调整参数a 使 之趋于a + ，就能得到约束问题的最优解。上述是关于最早的罚函数到乘子罚 函数的发展概述，现在有仍许多学者在这个领域继续研究 2 1 2 2 】 2 9 【3 9 4 0 】 在这一节中，介绍其中的伴随值作为增广l a g r a n g e 乘子法中的乘子，计算 增广l a g r a n g e 函数的梯度在自动微分中，对应的将一阶逆向模式中的输入 和输出变量分别替换成输入扰动d y + 和微分输出扰动d x + ，或简单置换线性 模式中输入扰动和输出扰动向量的位置，就得到如下形式的伴随模式： d x = v z t 。f d y + 其中d x 和d y + 分别为n 维和m 维向量。 显然，伴随模式只是线性模式的一种简单“转置”形式从输入输出关系 来看，后者反映了输入扰动对每个输出扰动的集中影响程度，而前者则反映 了输出扰动对每个输入扰动的依赖程度( 敏感性) 线性模式具有相对固定的 程序结构和计算代价，而伴随模式由于实现过程中在处理基态值重复计算和 数据存取上的差异，其形式则具有多样性且计算代价各异伴随模式常用来 计算函数梯度，具有与独立变元数目无关的理想计算代价 2 8 第4 章基于自动微分的最优化方法 4 1 1 伴随值作为增广l a g r a n g e 乘子 在第二章中介绍了自动微分逆向模式的基本形式，在这里将中间变量仇 视为具有等式约束忱( u i ) 一= 0 、目标函数为鼢这样一个等式约束问题的 l a g r a n g e 乘子，即为： m i n9 y ( 4 - 3 ) s t 9 i ( u i l 一玩= 0 设吻，j = 1 一礼，f 为相互独立的最优解，并且满足约束 e ( z ，) = ( 妒i ( u i ) 一耽) 仁1 一，l = 0 ，妒i ( u i ) 兰z 件。，i = 1 1 7 , ，0 ( 4 - 4 ) 这个广义系统是亩( 3 6 ) 给出的。现在将 m l ，y 三勋= 孰埘一l ( 4 5 ) i = 0 作为目标函数对于固定的z 和雪，k k t 条件中要求j a c o b i a n 矩阵e 7 = 石o e 的