Julia 拥有一个完整的生态系统，用于生成稀疏模式并以与科学计算和机器学习（或科学机器学习）相结合的方式进行稀疏自动微分。SparseDiffTools.jl、ModelingToolkit.jl和SparsityDetection.jl等工具将执行以下操作：

• 从代码中自动查找稀疏模式
• 生成解析雅可比表达式并自动并行化稀疏非分配雅可比表达式
• 为正向和反向模式 AD 进行矩阵着色

可以在此处找到一个集成示例，该示例使用稀疏性自动加速 ODE 求解以实现 55 倍加速。

为了看看实际效果如何，让我们对拉普拉斯方程进行简单的离散化：

fcalls = 0
function f(dx,x) # in-place
  global fcalls += 1
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end

我在那里放了一个小函数计数器来演示它是如何工作的。我们可以使用 SparsityDetection.jl 生成稀疏模式：

using SparsityDetection, SparseArrays
input = rand(10)
output = similar(input)
sparsity_pattern = jacobian_sparsity(f,output,input)
jac = Float64.(sparse(sparsity_pattern))

我们得到了我们都知道和喜爱的三对角矩阵。从这里我们执行矩阵着色：

using SparseDiffTools
colors = matrix_colors(jac)

因为maximum(colors)是 3，这意味着有限差分只需要 4 次函数评估来计算完整的雅可比行列式（要了解这一切是如何工作的，请参阅MIT 18.337 Parallel Computing and Scientific Machine Learning 讲座笔记，特别是关于前向模式 AD 的部分和求解刚性 ODE）。因此，我们可以通过以下方式快速计算整个稀疏雅可比行列式：

using FiniteDiff
FiniteDiff.finite_difference_jacobian!(jac, f, rand(30), colorvec=colors)
@show fcalls # 5

请注意，完整的函数调用为 5，因为自动稀疏检测f通过抽象解释使用虚假调用来生成稀疏模式。

然后，我们可以通过以下方式将前向模式 AD 用于稀疏模式：

forwarddiff_color_jacobian!(jac, f, x, colorvec = colors)

只需要总共 3 次f调用即可生成完整的雅可比行列式。包FiniteDiff.jl和SparseDiffTools.jl允许预缓存所有计算组件，因此您可以使这比演示更快，并使整个内部循环完全不分配。

请注意，反向模式 AD 的矩阵着色是 via matrix_colors(jac')，然后可以使用Zygote.jl、ReverseDiff.jl等用于稀疏反向模式。

但正如@chennaK 提到的，稀疏自动微分仍然会有一些开销。为了获得完全最优的结果，我们可以使用ModelingToolkit.jl生成完整漂亮的稀疏（和并行化代码）。我们可以通过抽象解释从我们的代码生成符号数学模型：

using ModelingToolkit
@variables u[1:10] du[1:10]
f(du,u)
du

10-element Array{Operation,1}:
        -2u₁ + u₂
  (u₁ - 2u₂) + u₃
  (u₂ - 2u₃) + u₄
  (u₃ - 2u₄) + u₅
  (u₄ - 2u₅) + u₆
  (u₅ - 2u₆) + u₇
  (u₆ - 2u₇) + u₈
  (u₇ - 2u₈) + u₉
 (u₈ - 2u₉) + u₁₀
        u₉ - 2u₁₀

现在我们可以使用sparsejacobian生成稀疏雅可比行列式的符号表达式：

sparsejac = ModelingToolkit.sparsejacobian(du,u)

然后我们可以告诉它生成一个快速的、非分配的、多线程的 Julia 代码：

build_function(sparsejac,u,parallel=ModelingToolkit.MultithreadedForm())[2]

它在此处生成代码，您可以eval在需要的任何其他代码中使用这些代码。这可以扩展到至少几百万个输入，因此我们在AutoOptimize.jl中使用它来执行用户代码的自动优化。

在 Julia 中完成这一切的好处是 Julia 将能够从所有这些调用中生成非常高效的机器代码，这意味着它更符合 C++ 而不是 Python。这方面的一个证明是纯 Julia 中的刚性 ODE 求解器比 CVODE 等 C++ 方法的性能高出 5x，因此在某种意义上，尽管 Julia 是一种高级语言，这都是一个有趣、快速和交互式的稀疏 AD 代码生成示例，只是因为它很简单不代表不快！

我还想指出支持稀疏雅可比矩阵的MatlabAutoDiff。自己尝试过：可以在很短的时间内计算大型雅可比行列式（尝试使用 N=1e5）。