一种方法是使用
自动微分而不是符号微分;这是一种同时计算
f(
x)和
f′(
x)的方法。有一种非常酷的方法可以使用
双重数实现,我从{{link3:Dan "sigfpe" Piponi在自动微分方面的博客文章中学到了这种方法。你应该去阅读一下,但基本思想如下:不再使用实数(或者我们最喜欢的
Double
),而是定义一个新集合D,通过将一个新元素
ε附加到ℝ上来定义这个集合,使得
ε2=0。这就像我们通过添加一个新元素
i到ℝ上来定义复数ℂ,使得
i2=-1。(如果你喜欢代数,这就相当于说D=ℝ[x]/⟨x
2⟩。)因此,D的每个元素都是形如
a+
bε的形式,其中
a和
b是实数。双重数上的算术运算与你期望的一样:
- (a + bε) ± (c + dε) = (a + c) ± (b + d)ε; 并且
- (a + bε)(c + dε) = ac + bcε + adε + bdε2 = ac + (bc + ad)ε。
(由于 ε2 = 0,因此除法更加复杂,尽管你用于复数的共轭乘法技巧仍然有效;有关更多信息,请参见Wikipedia's explanation。)
现在,这些有什么用呢?直觉上,
ε 就像一个无穷小量,使您可以使用它来计算导数。事实上,如果我们使用不同的名称重写乘法规则,它将变成
(
f +
f′
ε)(
g +
g′
ε) =
fg + (
f′
g +
fg′)
ε
并且那里的
ε 的系数看起来很像
函数乘积求导法则!
那么,让我们来研究一大类函数的情况。由于我们忽略了除法,所以假设我们有一个定义在实数集上的函数 f: ℝ → ℝ,它由幂级数定义(可能是有限的,因此任何多项式都可以,例如 sin(x),cos(x) 和 e^x)。然后我们可以按照显然的方式定义一个新函数 f_D: D → D:我们不再添加实数,而是添加双线性数等等。然后我声称 f_D(x + ε) = f(x) + f′(x)ε。首先,我们可以通过归纳证明对于任何自然数 i,都有 (x + ε)^i = x^i + ix^(i-1)ε;这将为 f(x) = x^k 的情况建立我们的导数结果。在基本情况下,当 i = 0 时,这个等式显然成立。然后假设它对于 i 成立,我们就有
- (x + ε)i+1 = (x + ε)(x + ε)i 通过分解出一个(x + ε)的拷贝实现
- = (x + ε)(xi + ixi-1ε) 根据归纳假设得出
- = xi+1 + (xi + x(ixi-1))ε 根据双数乘法的定义
- = xi+1 + (i+1)xiε 简单代数可得。
我們希望得到這樣的結果。現在,考慮我們的冪級數f,我們知道
- f(x) = a0 + a1x + a2x2 + … + aixi + …
然後我們有
- fD(x + ε) = a0 + a1(x + ε) + a2(x + ε)2 + … + ai(x + ε)i + …
- = a0 + (a1x + a1ε) + (a2x2 + 2a2xε) + … + (aixi + iaixi-1ε) + … 根据引理得到
- = (a0 + a1x + a2x2 + … + aixi + …) + (a1ε + 2a2xε + … + iaixi-1ε + …) 通过交换律得到
- = (a0 + a1x + a2x2 + … + aixi + …) + (a1 + 2a2x + … + iaixi-1 + …)ε 通过提取出ε得到
- = f(x) + f′(x)ε 根据定义得到。
太好了!对于这种情况来说,双重数(至少是这种情况,但结果通常是正确的)可以为我们进行微分。我们所要做的就是将原始函数应用于不是实数
x,而是双重数
x+
ε,然后提取
ε的结果系数。我敢打赌你可以看到如何在Haskell中实现这一点:
data Dual a = !a :+? !a deriving (Eq, Read, Show)
infix 6 :+?
instance Num a => Num (Dual a) where
(a :+? b) + (c :+? d) = (a+c) :+? (b+d)
(a :+? b) - (c :+? d) = (a-c) :+? (b-d)
(a :+? b) * (c :+? d) = (a*c) :+? (b*c + a*d)
negate (a :+? b) = (-a) :+? (-b)
fromInteger n = fromInteger n :+? 0
abs _ = error "No abs for dual numbers."
signum _ = error "No signum for dual numbers."
differentiate :: Num a => (Dual a -> Dual a) -> (a -> a)
differentiate f x = case f (x :+? 1) of _ :+? f'x -> f'x
f :: Num a => a -> a
f x = 3*x^2 + 5*x + 9
f' :: Num a => a -> a
f' = differentiate f
然后,瞧:
*Main> f 42
5511
*Main> f' 42
257
正如Wolfram Alpha可以证实的那样,这就是正确答案。
关于这个东西的更多信息确实可以找到。我不是这方面的专家,我只是觉得这个想法很酷,所以我借此机会引用一下我读过的内容,并尝试证明一些简单的问题。Dan Piponi写了更多关于双线性数/自动微分的内容,包括其中一篇文章,他展示了一个更一般的构造,允许求偏导数。Conal Elliott在其中一篇文章中展示了如何类比地计算导数塔(f(x), f′(x), f″(x), …)。维基百科关于自动微分的文章详细介绍了一些其他方法。(这显然是一种“前向模式自动微分”,但是“反向模式”也存在,而且似乎更快。)
最后,这里有一篇关于自动微分的Haskell维基页面
(链接在此),其中包含一些文章和重要的
Hackage软件包!我从未使用过这些软件包,但似乎
Edward Kmett的ad
软件包是最完整的,可以处理多种不同的自动微分方式,并且他上传了该软件包
在撰写另一个Stack Overflow问题的答案之后。
我希望能补充一件事情。您说“然而,数据类型不应该表示函数(除了解析器)。”在这一点上我必须不同意 - 将函数转化为数据类型对于这种情况非常有用。(那么解析器有什么特别之处呢?)无论何时您想要内省一个函数,将其转化为数据类型都是一个很好的选择。例如,这里是一个符号微分的编码,与上面的自动微分编码类似:
data Symbolic a = Const a
| Var String
| Symbolic a :+: Symbolic a
| Symbolic a :-: Symbolic a
| Symbolic a :*: Symbolic a
deriving (Eq, Read, Show)
infixl 6 :+:
infixl 6 :-:
infixl 7 :*:
eval :: Num a => (String -> a) -> Symbolic a -> a
eval env = go
where go (Const a) = a
go (Var x) = env x
go (e :+: f) = go e + go f
go (e :-: f) = go e - go f
go (e :*: f) = go e * go f
instance Num a => Num (Symbolic a) where
(+) = (:+:)
(-) = (:-:)
(*) = (:*:)
negate = (0 -)
fromInteger = Const . fromInteger
abs = error "No abs for symbolic numbers."
signum = error "No signum for symbolic numbers."
differentiate :: Num a => Symbolic a -> String -> Symbolic a
differentiate f x = go f
where go (Const a) = 0
go (Var y) | x == y = 1
| otherwise = 0
go (e :+: f) = go e + go f
go (e :-: f) = go e - go f
go (e :*: f) = go e * f + e * go f
f :: Num a => a -> a
f x = 3*x^2 + 5*x + 9
f' :: Num a => a -> a
f' x = eval (const x) $ differentiate (f $ Var "x") "x"
再一次:
*Main> f 42
5511
*Main> f' 42
257
这两种解决方案的美妙之处(或其中一部分)在于,只要您的原始
f
是多态的(类型为
Num a => a -> a
或类似类型),您就永远不需要修改
f
!您唯一需要放置导数相关代码的地方是在定义新数据类型和差异化函数中;您可以免费获得现有函数的导数。