机器如何求导?编程是一门艺术,它在与熵增和人类的愚蠢作斗争

许多导数其实很容易求,尤其是多项式函数的导数。对于这些情况,我们可以直接用幂次法则(power rule)。这是一个方便的简写方法,可以手算导数。在物理计算中,导数经常出现,因此让计算机能计算导数非常有用。但是,计算机是如何确定一个函数的导数的呢?导数的计算方式取决于函数的表示形式。如果是数值表示(即一组数据点序列),所用方法会与显式表示(数学符号)不同。显式表示本质上就是把数学函数改写成编程中的函数,例如将

要在某个特定点上求导数,我们需要创建一个新的函数来表示这个导数。这很直接:只需进行符号微分,然后把结果转换成所需编程语言中的函数即可。然而,当涉及多个函数时,这种方法就会变得笨拙,因为每次要求导都要实现一个新的函数。如果我们需要对几十甚至几百个函数求导,这会很快变得令人厌烦。更糟糕的是,如果函数需要用到链式法则(chain rule)或其他复杂法则,计算过程会迅速变得复杂。例如,假设我们要求以下函数的导数:

它的导数表达式相对简单:

但如果我们要求两个函数的乘积的导数,就会得到下面的结果:

函数从一个单一的乘积项变成了两个乘积项。如果被相乘的项数增加,导数中的项数也会随之增加:

可以看到,导数比原函数更复杂。而这仅仅是在简单多项式中使用乘法和乘积法则时的情况。同样的事情在使用链式法则时也会发生:

函数中嵌套的项越多,导数就越复杂。这意味着,在特定数值下计算导数所需的运算更多。这种现象被称为 表达式膨胀(expression swell):在计算过程中,表达式中的数字和项的数量不断增加。正如前面的例子所示,计算 f(5) 要比计算 f′(5) 容易得多。那么,我们能不能换一种方式来求导呢?幸运的是,我们有一个替代方案:对偶数(dual numbers)自动微分(automatic differentiation)

什么是对偶数?

在数学上,对偶数是一种超复数系统(hypercomplex number system)。在实际意义上,它们类似于复数(虚数),但性质不同。复数的形式是 a + b*i,其中 i^2 = -1,而对偶数的形式是 a + bε,其中 ε^2 = 0 且 ε ≠ 0。我们可以把 ε 看作是加在实数上的一个无穷小值(infinitesimal)。它的代数运算几乎与复数相同,只是当相乘时,我们将 ε^2 替换为 0。

对偶数有一个很巧妙的性质:只要这个数是 x + ε 的形式,我们就可以把它代入任意函数,从而在一次运算中同时得到该函数的值和它的导数。怎么做到的呢?例如,我们想要求当 x = 5 时的斜率(slope)。我们只需用对偶数 5 + ε 代入表达式,就会得到另一个对偶数:

结果中的对偶部分(ε 的系数)就是 x = 5 时导数的值,在这个例子中是 13。这一点可以通过手动求导并代入 x = 5 来验证。而结果中的实部则是函数在 x = 5 时的值。也就是说,在一次计算中,我们同时得到了函数值和导数值。

更有趣的是,这对任何多项式乃至一般函数都成立。我们可以用以下证明来说明原因:

所有多项式表达式都可以写成以下形式:

将对偶数 x + ε 代入,得到:

我们对等式右边应用二项式定理:

结果是:

稍加简化:第二项的系数是 n,因为 “n 取 1” 总是 n。而任何 ε 的指数大于 1 的项都会消失,因为它们包含 ε^2,而 ε^2 = 0。这样,上式就简化成:

这个第二项看起来已经很眼熟了,但我们继续。将它代回多项式表达式的原公式:

可以看到,把 x + ε 形式的对偶数代入任意多项式,得到的对偶数的实部就是函数值,而 ε 的系数就是该点处的导数值。这不仅对多项式成立,对任意函数也成立。

假设我们想在某个特定 x 处求一个三角函数的导数:

将 x + ε 代入,我们得到:

虽然这个表达式看起来比多项式的情况更难化简,但我们依然可以用泰勒展开(Taylor series)来证明。泰勒展开是关于某个特定点对函数进行近似的无限项和。它的公式为:

令 h = ε,我们得到:

由于 ε 的高次项全部消失(因为它们等于 0),一次项 ε 的系数仍然给出 f(x) 的导数。有趣的是,如果我们稍作整理,会得到:

它看起来和导数的极限定义非常相似,只是这里没有取极限(严格来说,我们不能除以 ε,因为这会导致矛盾:0 = 0/ε = ε^2/ε = ε,而根据定义 ε ≠ 0,但思想是类似的)。

在之前的例子中,我们得到了一个等式,证明了 sin(x) 的导数是 cos(x)。

既然我们已经知道什么是对偶数,以及它们如何用来计算导数,那么它们与自动微分(automatic differentiation)以及计算机中导数的计算又有什么关系呢?

自动微分,尤其是正向模式自动微分(forward-mode automatic differentiation),其实就是在算法/计算层面上实现对偶数的方法。它的做法是使用一种对象或其他数据类型,让它的行为与对偶数相同。实现方式之一,就是创建一个 Dual 类,用来存储对偶数的实部与对偶部,并在其上实现各种数学运算,就像处理复数一样。我们可以在 Python 中看到一个示例:每个数学运算都通过双下划线方法(dunder method)来实现。

class Dual:

 def __init__(self, real, dual):
  self.real = real
  self.dual = dual

 def __str__(self) -> str:
  sign = "+" if self.dual > 0 else ""
  return f"{self.real}{sign}{self.dual}{chr(949)}"

 def __neg__(self):
  return Dual(-self.real, -self.dual)
  
 def __add__(self, number):
  if isinstance(number, (int, float)):
   return Dual(self.real + number, self.dual)
  if isinstance(number, Dual):
   return Dual(self.real + number.real, self.dual + number.dual)
  return NotImplemented
  
 def __radd__(self, number):
  return self + number

 def __sub__(self, number):
  if isinstance(number, (int, float)):
   return Dual(self.real - number, self.dual)
  if isinstance(number, Dual):
   return Dual(self.real - number.real, self.dual - number.dual)
  return NotImplemented

 def __rsub__(self, number):
  return -self + number

 def __mul__(self, number):
  if isinstance(number, (int, float)):
   return Dual(self.real*number, self.dual*number)
  if isinstance(number, Dual):
   coeff1 = self.real*number.real
   coeff2 = self.real*number.dual + self.dual*number.real
   return Dual(coeff1, coeff2)
  return NotImplemented

 def __rmul__(self, number):
  return self*number

 def __truediv__(self, number):
  if isinstance(number, (int, float)):
   return Dual(self.real/number, self.dual/number)
  if isinstance(number, Dual):
   return self*number.conjugate()/number.real**2
  return NotImplemented

 def __rtruediv__(self, number):
  return number * self.conjugate()/self.real**2

 def __pow__(self, exponent):
   if isinstance(exponent, Dual):
     raise ValueError("Powers between two dual numbers is not defined.")
   new_real = self.real ** exponent
   new_dual = exponent * (self.real ** (exponent - 1)) * self.dual
   return Dual(new_real, new_dual) 

 def conjugate(self):
  return Dual(self.real, -self.dual)

有了上面的类,实现我们最初那个函数的求导就很简单:只需要创建一个 Dual 对象,把它传进函数里,然后读取结果即可。例如,如果我们要在 x = 5 处求导,只需要运行以下代码:

a = Dual(5, 1)
result = f(a)
print(result)

如果你是从算法的角度思考,可能会问:这个算法的计算复杂度是多少?或者,根本的算法到底是什么?与我们之前探讨过的其他算法不同,这里似乎没有明显的“算法步骤”——它更像是一种软件设计选择(创建一个类)并用它来做一些计算。

但实际上,这里面确实有一个算法。当我们探讨霍纳法则(Horner’s method)来在特定值处计算多项式时,发现,通过一些巧妙的代数变形,这可以在 O(n) 步内完成。这一次也是同样的原理,只不过我们传入的不是单个数值,而是一个包含两个值的对象。两个值会被并行计算,从而让我们在 O(n) 步内同时得到函数值和导数值。对偶数就是实现这种并行计算的工具。

这种微分方法的一个优点是:它给出的导数值是精确的,不像数值算法那样可能存在近似误差,同时还能避免我们在符号微分时遇到的表达式膨胀问题。

需要注意的是,对偶数只有在我们有实际函数表达式时才有效。如果只有一组数值序列(数值采样点),就必须使用有限差分(finite differences)或其他数值算法。