如何处理Haskell中的表达式?

9

假设我有:

f :: Double -> Double
f x = 3*x^2 + 5*x + 9

我想要计算这个函数的导数并写下来。
derivate f

为了让...
derivate f == \x -> 6*x + 5

但是如何定义派生呢?
derivate :: (a -> a) -> (a -> a)
derivate f = f' -- how to compute f'?

我知道没有本地的方法可以做到这一点,但是有没有可以使用的库呢?

我们是否必须依靠“元”数据类型来实现此功能?

data Computation = Add Exp Expr | Mult Expr Expr | Power Expr Expr -- etc

那么,为每个函数制作相应的构造函数不是很繁琐吗?然而,数据类型不应该表示函数(解析器除外)。

Pure 是否是一个好的替代方案,因为它具有术语重写功能?它难道没有缺点吗?

列表是否负担得起?

f :: [Double]
f = [3, 5, 9]

derivate :: (a -> [a])
derivate f = (*) <$> f <*> (getNs f)

compute f x = sum $
    ((*) . (^) x) <$> (getNs f) <*> f

getNs f = (reverse (iterate (length f) [0..]))

现在Haskell看起来像依赖于LISP的语法不太合适。等待一起使用的函数和参数相当于存储在数据类型中。此外,它不是很自然。
它们似乎不够“灵活”,不能像有理函数那样使用我的derivate函数。
例如,现在我想在游戏中使用导数。角色在使用函数制作的地板上奔跑,如果地面足够陡峭,我希望他滑行。
我还需要解决各种问题的方程式。以下是一些示例:
我是一艘宇宙飞船,想要小睡一下。在我的睡眠期间,如果我没有仔细放置自己,由于重力可能会在行星上坠毁。我没有足够的燃料远离天体,也没有地图。因此,我必须将自己放置在这个区域的物体之间,以使它们对我产生的引力总和被取消。 xy是我的坐标。 gravity是一个函数,它接受两个对象并返回它们之间的引力矢量。如果有两个物体,比如地球和月球,除了我之外,我只需要解决:
gravity earth spaceship + gravity moon spaceship == (0, 0)

相比于从头开始创建一个新的函数equigravityPoint :: Object -> Object -> Object -> Point,使用已有的函数会更加简单快捷。

如果除了我还有3个对象,那仍然很简单。

gravity earth spaceship + gravity moon spaceship + gravity sun spaceship == (0, 0)

同样适用于4和n。使用这种方法处理对象列表比使用“equigravityPoint”更简单。
另一个例子。 我想编写一个敌人机器人来射击我。 如果他只是瞄准我的当前位置射击,我向他跑时他会击中我,但如果我跳起来落在他身上,他就会错过我。 一个更聪明的机器人会这样想:"嗯,他从墙上跳下来了。如果我瞄准他现在所在的位置开枪,子弹不会击中他,因为到那时他会移动。所以我要预测他几秒钟后会在哪里,然后瞄准那里开枪,让子弹和他同时到达这个点"。 基本上,我需要计算轨迹的能力。例如,对于这种情况,我需要解决方程“trajectoryBullet == trajectoryCharacter”,找到直线和抛物线相交的点。
一个类似但更简单的例子,不涉及速度。 我是一个消防机器人,有一栋大楼着火了。另一组消防员正在用水枪灭火。我和其他人正在建筑物周围。当我的朋友们在喷水时,我拿着蹦床。 我需要在人们掉下来之前去他们将要着陆的地方。所以我需要轨迹和方程求解能力。

1
你可能可以用模板Haskell来完成这个任务,但是需要付出很多的努力! - pat
1
听起来你想要一个符号代数系统,例如muPad、Axiom、Mathematica、Maple等。但是如果你只想计算地板的陡度,那么数值方法可能就足够了。 - Theodore Norvell
1
好的,你应该明天来纽约,我们可以一起学习关于一个据称(我在互联网上找不到任何信息)由Greg Wright开发的Haskell新符号代数库"Wheeler"。 - jberryman
2
我相信你要找的是“自动微分”。Conal Elliott的博客是一个不错的起点:http://conal.net/blog/posts/what-is-automatic-differentiation-and-why-does-it-work - John L
1
相关问题:https://dev59.com/rkbRa4cB1Zd3GeqP1Y_f - ShreevatsaR
显示剩余2条评论
3个回答

29

一种方法是使用自动微分而不是符号微分;这是一种同时计算f(x)和f′(x)的方法。有一种非常酷的方法可以使用双重数实现,我从{{link3:Dan "sigfpe" Piponi在自动微分方面的博客文章中学到了这种方法。你应该去阅读一下,但基本思想如下:不再使用实数(或者我们最喜欢的Double),而是定义一个新集合D,通过将一个新元素ε附加到ℝ上来定义这个集合,使得ε2=0。这就像我们通过添加一个新元素i到ℝ上来定义复数ℂ,使得i2=-1。(如果你喜欢代数,这就相当于说D=ℝ[x]/⟨x2⟩。)因此,D的每个元素都是形如a+的形式,其中ab是实数。双重数上的算术运算与你期望的一样:
  • (a + ) ± (c + ) = (a + c) ± (b + d)ε; 并且
  • (a + )(c + ) = ac + bcε + adε + bdε2 = ac + (bc + ad)ε

(由于 ε2 = 0,因此除法更加复杂,尽管你用于复数的共轭乘法技巧仍然有效;有关更多信息,请参见Wikipedia's explanation。)

现在,这些有什么用呢?直觉上,ε 就像一个无穷小量,使您可以使用它来计算导数。事实上,如果我们使用不同的名称重写乘法规则,它将变成
(f + fε)(g + gε) = fg + (fg + 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 + 2a2) + … + (aixi + iaixi-1ε) + … 根据引理得到
  • = (a0 + a1x + a2x2 + … + aixi + …) + (a1ε + 2a2 + … + 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 and signum might actually exist, but I'm not sure what they are.
  abs    _              = error "No abs for dual numbers."
  signum _              = error "No signum for dual numbers."

-- Instances for Fractional, Floating, etc., are all possible too.

differentiate :: Num a => (Dual a -> Dual a) -> (a -> a)
differentiate f x = case f (x :+? 1) of _ :+? f'x -> f'x

-- Your original f, but with a more general type signature.  This polymorphism is
-- essential!  Otherwise, we can't pass f to differentiate.
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
  -- Ignoring abs and signum again
  abs         = error "No abs for symbolic numbers."
  signum      = error "No signum for symbolic numbers."

-- Instances for Fractional, Floating, etc., are all possible too.

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!您唯一需要放置导数相关代码的地方是在定义新数据类型和差异化函数中;您可以免费获得现有函数的导数。

这太酷了,而且你实际上包含了一个工作的(而且非常简短!)实现,使得这个答案真的很棒。谢谢。 - Niklas B.
很棒的答案!我想指出:仅计算导数值时,自动微分要好得多-这样做符号化可能会导致一些效率问题,特别是涉及“let”语句时。这并不是说我反对更多的符号方法,但通常最好为您要计算的内容使用专门的数据类型。正如您所说,由于多态性,您可以毫不费力地以多种不同的方式评估相同的表达式 - leftaroundabout
非常感谢自动微分。我可能会将其与符号计算相结合,用于其他目的,例如方程求解(我更新了问题以提供此类用途的示例)。 - L01man
我不喜欢符号计算和“常规”计算分开的方式,因为这会导致我认为可以避免的样板文件。这就是为什么我会看看Maxima(Mathematica不是开源的)。我想知道Haskell是否可以成为这些CAS语言之一...而解析器的特殊之处在于,您无法避免使用数据类型来表示解析的语言的语法。 - L01man

4

数值求导可以很容易地完成:

derive f x = (f (x + dx) - f (x - dx)) / (2 * dx) where dx = 0.00001

然而,对于符号导数,您需要创建一个AST,然后通过匹配和重写AST来实现求导规则。


3
要是你能在那里加入当dx趋近于零时的极限就好了! - pat

2
我不理解您对使用自定义数据类型的问题所在。
data Expr = Plus Expr Expr 
           | Times Expr Expr 
           | Negate Expr 
           | Exp Expr Expr 
           | Abs Expr
           | Signum Expr
           | FromInteger Integer
           | Var

instance Num Expr where
   fromInteger = FromInteger
   (+) = Plus
   (*) = Times
   negate = Negate
   abs = Abs
   signum = Signum

toNumF :: Num a => Expr -> a -> a
toNumF e x = go e where
  go Var = x
  go (FromInteger i) = fromInteger i
  go (Plus a b) = (go a) + (go b)
  ...

你可以像使用 IntDouble 一样使用此方法,所有操作都会正常运行!你可以定义一个函数。

deriveExpr :: Expr -> Expr

然后,您就可以定义以下(RankN)函数。
derivate :: Num b => (forall a. Num a => a -> a) -> b -> b
derivate f = toNumF $ deriveExpr (f Var)

你可以将此扩展到与数值层次结构的其他部分配合使用。

它可以正常工作,但我不喜欢例如Plus作为(+)的重复。也许这是Haskell或库的问题,但我认为您应该能够随意使用表达式1 + 3x。由于引用透明性,Haskell计算它并且不允许您操作其术语。您需要一个特殊的结构,将许多现有函数重新定义为“值”!也许最好为函数设置两种模式:正常模式和代数模式,在代数模式下,您可以处理起始表达式而不是其结果。 - L01man
但你只需要编写这段代码一次。一旦你完成后,就可以使用任何类型的数字,编写 derivate (\x -> 2x*x) 并得到答案。 - Philip JF

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接