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

数值,计算,反向,自动,微分,python,实现 · 浏览次数 : 466

小编点评

value: 3.694528049465325 deriv: 9.236320123663312 value: 3.694528049465325 deriv: 9.236320123663312 value: 3.694528049465325 deriv: 9.236320123663312 value: 3.694528049465325 deriv: 9.236320123663312 value: 3.694528049465325 deriv: 9.236320123663312

正文

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实现)相似的内容:

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

自动微分技术(称为“automatic differentiation, autodiff”)是介于符号微分和数值微分的一种技术,它是在计算效率和计算精度之间的一种折衷。自动微分不受任何离散化算法误差的约束,它充分利用了微分的链式法则和其他关于导数的性质来准确地计算它们。我们可以选择定义一种新的数据类型来存储[u,u′]和[v,v′]这类数对。我们可以在这些数对上定义一种代数运算,这些代数运算编码了一些经典的操作。

[转帖]9.2 TiFlash 架构与原理

9.2 TiFlash 架构与原理 相比于行存,TiFlash 根据强 Schema 按列式存储结构化数据,借助 ClickHouse 的向量化计算引擎,带来读取和计算双重性能优势。相较于普通列存,TiFlash 则具有实时更新、分布式自动扩展、SI(Snapshot Isolation)隔离级别读

【RcoketMQ】RcoketMQ 5.0新特性(一)- Proxy

为了向云原生演进,提高资源利用和弹性能力,RcoketMQ在5.0进行了架构的调整与升级,先来看新特性之一,增加了Proxy层。 增加Proxy代理层 计算存储分离 计算存储分离是一种分层架构,将计算层与存储层分开。 计算层指的是一些消耗计算资源的功能模块比如协议解析、消费管理等,存储指的是数据存储

记一次 .NET某网络边缘计算系统 卡死分析

一:背景 1. 讲故事 早就听说过有什么 网络边缘计算,这次还真给遇到了,有点意思,问了下 chatgpt 这是干嘛的 ? 网络边缘计算是一种计算模型,它将计算能力和数据存储位置从传统的集中式数据中心向网络边缘的用户设备、传感器和其他物联网设备移动。这种模型的目的是在接近数据生成源头的地方提供更快速

WebGPU缓冲区更新最佳实践

介绍 在WebGPU中,GPUBuffer是您将要操作的主要对象之一。它与GPUTextures一同代表了您的应用程序向GPU传递用于渲染的大部分数据。在WebGPU中,缓冲区用于顶点和索引数据、uniforms、计算和片段着色器的通用存储,以及作为纹理数据的临时存储区域。 本文档专注于找到将数据有

【numpy基础】--通用计算

`numpy`提供了简单灵活的接口,用于优化数据数组的计算。 通用计算最大的优势在于通过向量化操作,将循环推送至`numpy`之下的编译层,从而取得更快的执行效率。 `numpy`的通用计算让我们计算数组时就像计算单独一个变量一样, 不用写循环去遍历数组中的各个元素。 比如,对于一般的`python

半同态计算芯片

半同态计算芯片 学习该文章:华控清交推出业界首款半同态计算芯片 赋能隐匿查询实用化 摘要 隐匿查询是指在不向数据提供方暴露查询方的查询意图,同时又能在保护数据提供方数据库中其他数据的情况下让查询方获得相关查询结果。实际使用的场景大多是跨广域网环境下基于关键字的查询。当前的常用方法要么需要传输大量的数

数据分析---numpy模块

前戏 NumPy(Numerical Python) 是 Python 语言中做科学计算的基础库。重在于数值计算,也是大部分Python科学计算库的基础,多用于在大型、多维数组上执行的数值运算。 快捷键的使用: 添加cell:a或者b 删除:x 修改cell的模式: m:修改成markdown模式

2023年最具威胁的25种安全漏洞(CWE TOP 25)

摘要: CWE Top 25 是通过分析美国国家漏洞数据库(NVD)中的公共漏洞数据来计算的,以获取前两个日历年 CWE 弱点的根本原因映射。 本文分享自华为云社区《2023年最具威胁的25种安全漏洞(CWE TOP 25)》,作者: Uncle_Tom 。 CWE Top 25 是通过分析美国国家

基于阿里云服务实现短信验证码功能

## 前言: 阿里云短信服务是一项基于云计算和大数据技术的企业级短信平台服务。它能够为企业和开发者提供高可用、高性能、高稳定性的短信发送服务,可以快速地将各类业务通知、验证码、营销推广等信息发送给用户。在我们经常登录一些系统或者APP时候,经常会遇到其他登录登录方式——短信验证码登录。这也是我前一段