今年夏天，我有幸参加了首届 **Julia 暑期代码 (JSoC)**，该项目得到了 Gordon 和 Betty Moore 基金会的慷慨赞助。我的 JSoC 项目是探索 Julia 在自动微分 (AD) 中的应用，AD 是优化领域中广泛应用的一个主题。

在 Miles Lubin 和 Theodore Papamarkou 的指导下，我完成了 **ForwardDiff.jl** 的重大改版，这是一个用于计算原生 Julia 函数（或任何可调用的 Julia 类型）的导数、梯度、雅可比矩阵、海森矩阵和高阶导数的 Julia 包。

在本文结束时，您希望能够了解 ForwardDiff.jl 的工作原理、它为何有用以及为什么与其他语言相比，Julia 非常适合 AD。

什么是自动微分？

广义上讲，自动微分 描述了一类用于自动获取用户提供函数的精确导数的算法。除了产生更准确的结果之外，AD 方法通常也比其他常见的微分方法（例如 有限差分）更快。

AD 的两种主要类型称为前向模式和反向模式。正如您可能已经猜到的，本文只讨论前向模式，它是 ForwardDiff.jl 实现的 AD 类型。

在行动中看到 ForwardDiff.jl

在我们深入了解具体细节之前，先看一个简单的示例，它演示了 ForwardDiff.jl API 中的各种方法，可能会很有帮助。

下面的代码段是一个有点人为的示例，但足以作为对该包的介绍。首先，我们定义要微分的目标函数，然后使用 ForwardDiff.jl 计算该函数在给定输入处的某些导数

julia> using ForwardDiff

julia> f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);

julia> x = rand(5)
5-element Array{Float64,1}:
 0.986403
 0.140913
 0.294963
 0.837125
 0.650451

julia> g = ForwardDiff.gradient(f); # g = ∇f

julia> g(x)
5-element Array{Float64,1}:
 1.01358
 2.50014
 1.72574
 1.10139
 1.2445

julia> j = ForwardDiff.jacobian(g); # j = J(∇f)

julia> j(x)
5x5 Array{Float64,2}:
 0.585111  3.48083  1.7706    0.994057  1.03257
 3.48083   1.06079  5.79299   3.25245   3.37871
 1.7706    5.79299  0.423981  1.65416   1.71818
 0.994057  3.25245  1.65416   0.251396  0.964566
 1.03257   3.37871  1.71818   0.964566  0.140689

julia> ForwardDiff.hessian(f, x) # H(f)(x) == J(∇f)(x), as expected
5x5 Array{Float64,2}:
 0.585111  3.48083  1.7706    0.994057  1.03257
 3.48083   1.06079  5.79299   3.25245   3.37871
 1.7706    5.79299  0.423981  1.65416   1.71818
 0.994057  3.25245  1.65416   0.251396  0.964566
 1.03257   3.37871  1.71818   0.964566  0.140689

太棒了！

好吧，这并**不**太令人兴奋 - 我可以用 Calculus.jl 做同样的事情。我为什么要使用 ForwardDiff.jl 呢？

简单的答案是，在许多情况下，ForwardDiff.jl 的基于 AD 的方法比其他包中实现的有限差分方法性能要高得多。

ForwardDiff.jl 的工作原理 - 概述

ForwardDiff.jl 利用的关键技术是实现了几种不同的 `ForwardDiffNumber` 类型，每种类型都为正常值和导数值分配存储空间。然后，重载 `ForwardDiffNumber` 上的基本数值函数以评估原始函数及其导数，并将结果以新的 `ForwardDiffNumber` 的形式返回。

因此，我们可以将这些数字类型传递到一个通用函数 $f$（假设它由重载的基本函数组成）中，然后通过链式法则在计算的每个步骤中自然地传播导数信息。评估的最终结果（通常是 `ForwardDiffNumber` 或其数组）然后包含 $f(x)$ 和 $f'(x)$，其中 $x$ 是原始的评估点。

Julia 中的简单前向模式 AD

编写演示此技术的实际 Julia 代码的最简单方法是实现一个简单的 对偶数 类型。请注意，已经有一个专门用于此类实现的 Julia 包，但出于教学目的，我们将在这里自己编写。

以下是我们定义 `DualNumber` 类型的方式

immutable DualNumber{T} <: Number
    value::T
    deriv::T
end

value(d::DualNumber) = d.value
deriv(d::DualNumber) = d.deriv

接下来，我们可以开始定义 `DualNumber` 上的函数。以下是一些示例，让您了解此过程

function Base.sqrt(d::DualNumber)
    new_value = sqrt(value(d))
    new_deriv = 0.5 / new_value
    return DualNumber(new_value, new_deriv*deriv(d))
end

function Base.sin(d::DualNumber)
    new_value = sin(value(d))
    new_deriv = cos(value(d))
    return DualNumber(new_value, new_deriv*deriv(d))
end

function Base.(:+)(a::DualNumber, b::DualNumber)
    new_value = value(a) + value(b)
    new_deriv = deriv(a) + deriv(b)
    return DualNumber(new_value, new_deriv)
end

function Base.(:*)(a::DualNumber, b::DualNumber)
    val_a, val_b = value(a), value(b)
    new_value = val_a * val_b
    new_deriv = val_b * deriv(a) + val_a * deriv(b)
    return DualNumber(new_value, new_deriv)
end

现在，我们可以评估由上述基本函数组成的任何标量函数的导数。为此，我们只需将 `DualNumber` 类型的实例传递到函数中，然后从结果中提取导数。例如

julia> f(x) = sqrt(sin(x * x)) + x
f (generic function with 1 method)

julia> f(1.0)
1.8414709848078965

julia> d = f(DualNumber(1.0, 1.0))
DualNumber{Float64}(1.8414709848078965,1.5403023058681398)

julia> deriv1 = deriv(d)
1.589002649374538

julia> using Calculus; deriv2 = Calculus.derivative(f, 1.0)
1.5890026493377403

julia> deriv1 - deriv2
3.679767601738604e-11

请注意，我们的对偶数结果接近从 Calculus.jl 获得的结果，但实际上略有不同。这种细微的差异是由于 Calculus.jl 采用的有限差分方法固有的近似误差造成的。

实际上，ForwardDiff.jl 提供的数字类型比 `DualNumber` 复杂得多。而不是简单的对偶数，各种 `ForwardDiffNumber` 类型表现得像对偶数和 超对偶数（对偶数的高阶类似物）的集合。这种基于集合的方法允许在一次评估目标函数时同时计算多个高阶偏导数。

性能比较：Ackley 函数

说明使用 ForwardDiff.jl 可以实现的性能提升的最佳方法是进行一些基准测试。让我们比较使用 ForwardDiff.jl、Calculus.jl 和基于 Python 的 AD 工具 AlgoPy 计算函数梯度所需的时间。

我们将在测试中使用的函数是 Ackley 函数，它在数学上定义为

function ackley(x)
    a, b, c = 20.0, -0.2, 2.0*π
    len_recip = inv(length(x))
    sum_sqrs = zero(eltype(x))
    sum_cos = sum_sqrs
    for i in x
        sum_cos += cos(c*i)
        sum_sqrs += i^2
    end
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + e)
end

def ackley(x):
    a, b, c = 20.0, -0.2, 2.0*numpy.pi
    len_recip = 1.0/len(x)
    sum_sqrs, sum_cos = 0.0, 0.0
    for i in x:
        sum_cos += algopy.cos(c*i)
        sum_sqrs += i*i
    return (-a * algopy.exp(b*algopy.sqrt(len_recip*sum_sqrs)) -
            algopy.exp(len_recip*sum_cos) + a + numpy.e)