首页 > 编程语言 >数值计算:前向和反向自动微分(Python实现)

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

时间:2022-12-28 16:12:19浏览次数:63  
标签:frac deriv val Python self 微分 other 反向 partial

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来运行优化算法时显得格外有价值,从而避免每次目标函数改变时都去重新手动计算导数。当然,做为其便捷性的代价,自动微分也会带来计算的效率问题,因为在实际工作中自动微分方法并不会去化简表达式,而是直接应用最显式的编码规则。

参考

标签:frac,deriv,val,Python,self,微分,other,反向,partial
From: https://www.cnblogs.com/orion-orion/p/17010353.html

相关文章

  • python中的集合推导式
    集合推导式可用来去重需求:将列表list1=[2,2,2,3,4,4,4]中的偶数进行筛选,并且去重list1=[2,2,2,3,4,4,4]set1={iforiinlist1ifi%2==0}print(set......
  • python中的列表推导式
    1.单列表,单条件求1-20之间的偶数list1=[]foriinrange(1,21):ifi%2==0:list1.append(i)print(list1)列表推导式list2=[iforiinrange(1,21)if......
  • Python爬虫(一)热身
    基础操作一importurllibimportchardet#字符集检测url="http://www.163.com/"html=urllib.urlopen(url)printhtml.headers#头部信息printhtml.getcode()#状态......
  • 当我把用Python做的课堂点名系统献给各科老师后,再也没挂过科
    刚上大学的表弟问我,大学准备好好玩玩,问我有没有什么不挂科的秘诀。哎,这可就问对人了,要想不挂科,先把老师贿赂好,当然,咱们说的贿赂不是送钱啥的,这不是侮辱老师吗?于是我......
  • python中resp.json()与json.loads(str)的区别
    resp=resquests.get(url)print(type(resp))#<class'requests.models.Response'>第一行代码使用requests库发送get请求,得到响应数据resp。第二行代码的输......
  • 读python代码-学到的python
    1.withopen(data_path,'r')asf:withopen()是python用来打开本地文件的,他会在使用完毕后,自动关闭文件,无需手动书写close().三种打开模式:r:只读 用read()w:只写用w......
  • python模块之psutil详解
     一、psutil模块:1.psutil是一个跨平台库(​​http://pythonhosted.org/psutil/​​)能够轻松实现获取系统运行的进程和系统利用率(包括CPU、内存、磁盘、网络等)信息。它主......
  • 【python】抽象基类 from abc import ABC, abstractmethod
    abc模块作用Python本身不提供抽象类和接口机制,要想实现抽象类,可以借助abc模块。ABC是Abstract BaseClass的缩写。假设我们定义一些抽象方法,然后子类继承的时候必须要重......
  • 逆向工程 Python 逆向
    逆向工程Python逆向Salarypython逆向https://github.com/SKPrimin/HomeWork/ReverseEngineering/lab1_python(选做)运行Salary.pyc,要求输出flag代表成功。直接运行......
  • python的list的用法
    #ReadMe#本工具是根据用户选择的条目来打印该列表下的内容#例如选择“北京”就会打印北京下面的“海淀”“昌平”“朝阳”,选择“海淀”然后会打印海淀下面的“清华大学”和......