草庐IT

数值计算:前向和反向自动微分(Python实现)

Orion's Blog 2023-03-28 原文

1 自动微分

我们在《数值分析》课程中已经学过许多经典的数值微分方法。许多经典的数值微分算法非常快,因为它们只需要计算差商。然而,他们的主要缺点在于他们是数值的,这意味着有限的算术精度和不精确的函数求值,而这些都从根本上限制了求解结果的质量。因此。充满噪声的、复杂多变的函数很难得到精准的数值微分。

自动微分技术(称为“automatic differentiation, autodiff”)是介于符号微分和数值微分的一种技术,它是在计算效率和计算精度之间的一种折衷。自动微分不受任何离散化算法误差的约束,它充分利用了微分的链式法则和其他关于导数的性质来准确地计算它们。

2 前向自动微分

我们先来计算简单的前向自动微分。假设我们有两个变量\(u\)\(v\),使用浮点数存储。我们将变量\(u′=du/dt\)\(v′=dv/dt\)和这些变量一起存储,这里\(t\)是独立的变量。在一些程序设计语言(如Python)中,我们可以选择定义一种新的数据类型来存储\([u,u′]\)\([v,v′]\)这类数对。我们可以在这些数对上定义一种代数运算,这些代数运算编码了一些经典的操作:

\[\begin{gathered} {\left[u, u^{\prime}\right]+\left[v, v^{\prime}\right] \equiv\left[u+v^{\prime}, u^{\prime}+v^{\prime}\right]} \\ c\left[u, u^{\prime}\right] \equiv\left[c u, c u^{\prime}\right] \\ {\left[u, u^{\prime}\right] \cdot\left[v, v^{\prime}\right] \equiv\left[u v, u v^{\prime}+u^{\prime} v\right]} \\ {\left[u, u^{\prime}\right] /\left[v, v^{\prime}\right] \equiv\left[u / v,\left(v u^{\prime}-u v^{\prime}\right) / v^2\right]} \\ \exp (\left[u, u^{\prime}\right]) \equiv\left[e^u, u^{\prime} e^u\right] \\ \ln (\left[u, u^{\prime}\right]) \equiv\left[\ln u_{,} u^{\prime} / u\right] \\ \cos (\left[u, u^{\prime}\right]) \equiv\left[\cos u,-u^{\prime} \sin u^{\prime}\right] \\ \vdots \quad\vdots \end{gathered} \]

在进行前向自动微分之前,我们需要先将计算\(f(t)\)所产生的操作序列表示为计算图。接着,采用自底向上的递推算法的思想,从做为递推起点的数对\(t≡[t_0,1]\)(因为\(dt/dt= 1\))开始,我们能够按照我们上述编码规则同时对函数\(f(t)\)和它的导数\(f′(t)\)进行求值。我们在编程语言中可以选择令数对重载运算符,这样额外的求导数运算就可以对用户透明地执行了。

例1 比如,对于函数\(f(x) = \exp(x^2 - x)/{x}\),想要依次计算\({dy}_i/dx\)(这里\(y_i\)为所有计算中间项)。则我们先从\(x\)开始将表达式分解为计算图:

\[\begin{aligned} & x \\ & y_1= x^2\\ & y_2=y_1 - x\\ & y_3 = \exp(y_2)\\ & y_4 = y_3/x \end{aligned} \]

然后前向递推地按照我们之前所述的编码规则来进行求导

\[\begin{aligned} & \frac{dy_1}{dx} = 2x\\ &\frac{dy_2}{dx} = \frac{dy_1}{dx} - \frac{dx}{dx} = 2x-1\\ & \frac{dy_3}{dx} = \exp(y_2)\cdot \frac{dy_2}{dx} \\ & \frac{dy_4}{dx} = \frac{\frac{dy_3}{dx}x - y_3}{x^2} \end{aligned} \]

注意链式法则(chain rule)告诉我们:

\[(f(g(x)))' = f'(g(x))\cdot g'(x) \]

所以我们对

\[y_k = g(y_i) \]

\[y'_k = g'(y_i)\cdot y_i' \]

事实上,我们也能够处理有多个输入的函数\(g\)

\[y_k = g(y_i,\cdots, y_j) \]

多元微分链式法则如下:

\[\begin{aligned} \frac{d}{dx} y_k(x) &= \frac{d}{dx} g(y_i(x),\cdots, y_j(x))\\ &= \sum_{h=i}^j\frac{\partial g}{\partial y_h} \frac{d y_h}{dx} \end{aligned} \]

比如,对于

\[\begin{aligned} & y_1 = x\\ & y_2 = x \\ & y_3 = y_2 \\ & y_4 = y_1\cdot y_2\cdot y_3 \\ \end{aligned} \]

我们有

\[\begin{aligned} \frac{dy_1}{dx} &=1 \\ \frac{dy_2}{dx} &= 1\\ \frac{dy_3}{dx} &= 1\cdot \frac{dy_2}{dx} = 1 \\ \frac{dy_4}{dx} &= y_2 y_3\cdot \frac{dy_1}{dx} + y_1 y_3\frac{dy_2}{dx} + y_1 y_2 \frac{dy_3}{dx}\\ &= y_2 y_3 + y_1 y_3 + y_1y_2 \\ &= 3x^2 \end{aligned} \]

下面展示了一个对二元函数模拟前向自动微分的过程。

例2\(f(x_1, x_2) = x_1\cdot \exp(x_2) - x_1\),模拟前向微分过程。

\[\begin{aligned} y_1 = \exp(x_2)\\ y_2 = x_1 \cdot y_1\\ y_3 = y_2 - x_1 \end{aligned} \]

\[\begin{aligned} & \frac{d y_1}{ d x_2} = \exp(x_2)\\ & \frac{d y_2}{d x_1} = y_1=\exp(x_2) \quad \frac{d y_2}{dx_2} = x_1 \cdot \frac{dy_1}{dx_2} = x_1\cdot \exp(x_2) \\ & \frac{d y_3}{d x_1} = \frac{dy_2}{dx_1} - \frac{dx_1}{dx_1} =\exp(x_2) -1 \quad \frac{dy_3}{dx_2} = \frac{dy_2}{dx_2} = x_1\cdot \exp(x_2) \\ \end{aligned} \]

接下来我们看如何用Python代码来实现单变量函数的前向自动微分过程。为了简便起见,我们下面只编码了几个常用的求导规则。

import math

class Var:
    def __init__(self, val, deriv=1.0):
        self.val = val
        self.deriv = deriv
    
    def __add__(self, other):
        if isinstance(other, Var):
            val = self.val + other.val
            deriv = self.deriv + other.deriv
        else:
            val = self.val + other
            deriv = self.deriv
        return Var(val, deriv)
    
    def __radd__(self, other):
        return self + other

    def __sub__(self, other):
        if isinstance(other, Var):
            val = self.val - other.val
            deriv = self.deriv - other.deriv
        else:
            val = self.val - other
            deriv = self.deriv
        return Var(val, deriv)
    
    def __rsub__(self, other):
        val = other - self.val
        deriv = - self.deriv
        return Var(val, deriv)

    def __mul__(self, other):
        if isinstance(other, Var):
            val = self.val * other.val
            deriv = self.val * other.deriv + self.deriv * other.val
        else:
            val = self.val * other
            deriv = self.deriv * other
        return Var(val, deriv)
    
    def __rmul__(self, other):
        return self * other

    def __truediv__(self, other):
        if isinstance(other, Var):
            val = self.val / other.val
            deriv = (self.deriv * other.val - self.val * other.deriv)/other.val**2
        else:
            val = self.val / other
            deriv = self.deriv / other
        return Var(val, deriv)

    def __rtruediv__(self, other):
        val = other / self.val
        deriv = other * 1/self.val**2
        return Var(val, deriv)
    
    def __repr__(self):
        return "value: {}\t gradient: {}".format(self.val, self.deriv)
        

def exp(f: Var):
    return Var(math.exp(f.val), math.exp(f.val) * f.deriv)

例如,我们若尝试计算函数\(f(x) = \exp(x^2 - x)/{x}\)\(x=2.0\)处的导数\(f'(2.0)\)如下:

fx = lambda x: exp(x*x - x)/x
df = fx(Var(2.0))
print(df) 

打印输出:

value: 3.694528049465325         deriv: 9.236320123663312

可见,前向过程完成计算得到\(f(2.0)\approx 3.69\), \(f'(2.0)\approx 9.24\)

3 反向自动微分

我们前面介绍的前向自动微分方法在计算\(y = f(t)\)的时候并行地计算\(f'(t)\)。接下来我们介绍一种“反向”自动微分方法,相比上一种的方法它仅需要更少的函数求值,不过需要以更多的内存消耗和更复杂的实现做为代价。

同样,这个技术需要先将计算\(f(t)\)所产生的操作序列表示为计算图。不过,与之前的从\(dt/dt = 1\)开始,然后往\(dy/dt\)方向计算不同,反向自动求导算法从\(dy/dy = 1\)开始并且按与之前同样的规则往反方向计算,一步步地将分母替换为\(dt\)。反向自动微分可以避免不必要的计算,特别是当\(y\)是一个多元函数的时候。例如,对\(f(t_1, t_2) = f_1(t_1) + f_2(t_2)\),反向自动微分并不需要计算\(f_1\)关于\(t_2\)的微分或\(f_2\)关于\(t_1\)的微分。

例3\(f(x_1, x_2) = x_1\cdot \exp(x_2) - x_1\),模拟反向自动微分过程。

\[\begin{aligned} y_1 = \exp(x_2)\\ y_2 = x_1 \cdot y_1\\ y_3 = y_2 - x_1 \end{aligned} \]

\[\begin{aligned} & \frac{\partial f}{\partial y_3} = 1\\ & \frac{\partial f}{\partial y_2} = \frac{\partial f}{\partial y_3}\frac{\partial y_3}{\partial y_2} = 1 \cdot 1 = 1\\ & \frac{\partial f}{\partial y_1} = \frac{\partial f}{\partial y_2} \frac{\partial y_2}{\partial y_1} = 1 \cdot x_1 = x_1\\ & \frac{\partial f}{\partial x_2} = \frac{\partial f}{\partial y_1} \frac{\partial y_1}{\partial x_2} = x_1 \cdot \exp(x_2)\\ & \frac{\partial f}{\partial x_1} = \frac{\partial f}{\partial y_2}\frac{\partial y_2}{x_1} + \frac{\partial f}{\partial y_3}\frac{\partial y_3}{\partial x_1} = 1\cdot y_1 + 1\cdot (-1) = \exp(x_2) - 1 \end{aligned} \]

可见若采用反向自动微分,我们需要存储计算过程中的所有东西,故内存的使用量会和时间成正比。不过,在现有的深度学习框架中,对反向自动微分的实现进行了进一步优化,我们会在深度学习专题文章中再进行详述。

4 总结

自动微分被广泛认为是一种未被充分重视的数值技术, 它可以以尽量小的执行代价来产生函数的精确导数。它在软件需要计算导数或Hessian来运行优化算法时显得格外有价值,从而避免每次目标函数改变时都去重新手动计算导数。当然,做为其便捷性的代价,自动微分也会带来计算的效率问题,因为在实际工作中自动微分方法并不会去化简表达式,而是直接应用最显式的编码规则。

参考

有关数值计算:前向和反向自动微分(Python实现)的更多相关文章

  1. ruby-on-rails - 使用 Ruby on Rails 进行自动化测试 - 最佳实践 - 2

    很好奇,就使用ruby​​onrails自动化单元测试而言,你们正在做什么?您是否创建了一个脚本来在cron中运行rake作业并将结果邮寄给您?git中的预提交Hook?只是手动调用?我完全理解测试,但想知道在错误发生之前捕获错误的最佳实践是什么。让我们理所当然地认为测试本身是完美无缺的,并且可以正常工作。下一步是什么以确保他们在正确的时间将可能有害的结果传达给您? 最佳答案 不确定您到底想听什么,但是有几个级别的自动代码库控制:在处理某项功能时,您可以使用类似autotest的内容获得关于哪些有效,哪些无效的即时反馈。要确保您的提

  2. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  3. ruby-on-rails - 如果为空或不验证数值,则使属性默认为 0 - 2

    我希望我的UserPrice模型的属性在它们为空或不验证数值时默认为0。这些属性是tax_rate、shipping_cost和price。classCreateUserPrices8,:scale=>2t.decimal:tax_rate,:precision=>8,:scale=>2t.decimal:shipping_cost,:precision=>8,:scale=>2endendend起初,我将所有3列的:default=>0放在表格中,但我不想要这样,因为它已经填充了字段,我想使用占位符。这是我的UserPrice模型:classUserPrice回答before_val

  4. ruby-on-rails - 使用一系列等级计算字母等级 - 2

    这里是Ruby新手。完成一些练习后碰壁了。练习:计算一系列成绩的字母等级创建一个方法get_grade来接受测试分数数组。数组中的每个分数应介于0和100之间,其中100是最大分数。计算平均分并将字母等级作为字符串返回,即“A”、“B”、“C”、“D”、“E”或“F”。我一直返回错误:avg.rb:1:syntaxerror,unexpectedtLBRACK,expecting')'defget_grade([100,90,80])^avg.rb:1:syntaxerror,unexpected')',expecting$end这是我目前所拥有的。我想坚持使用下面的方法或.join,

  5. ruby - RuntimeError(自动加载常量 Apps 多线程时检测到循环依赖 - 2

    我收到这个错误:RuntimeError(自动加载常量Apps时检测到循环依赖当我使用多线程时。下面是我的代码。为什么会这样?我尝试多线程的原因是因为我正在编写一个HTML抓取应用程序。对Nokogiri::HTML(open())的调用是一个同步阻塞调用,需要1秒才能返回,我有100,000多个页面要访问,所以我试图运行多个线程来解决这个问题。有更好的方法吗?classToolsController0)app.website=array.join(',')putsapp.websiteelseapp.website="NONE"endapp.saveapps=Apps.order("

  6. ruby - 如何根据特征实现 FactoryGirl 的条件行为 - 2

    我有一个用户工厂。我希望默认情况下确认用户。但是鉴于unconfirmed特征,我不希望它们被确认。虽然我有一个基于实现细节而不是抽象的工作实现,但我想知道如何正确地做到这一点。factory:userdoafter(:create)do|user,evaluator|#unwantedimplementationdetailshereunlessFactoryGirl.factories[:user].defined_traits.map(&:name).include?(:unconfirmed)user.confirm!endendtrait:unconfirmeddoenden

  7. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  8. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  9. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  10. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

随机推荐