如何在Haskell中实现卡尔曼矩阵？

Question

我正在尝试使用 Numeric.AD 库在 Haskell 中实现可微函数的卡尔曼矩阵。我正在使用 https://en.wikipedia.org/wiki/Carleman_matrix 作为参考。

到目前为止我有以下内容，

{-# LANGUAGE ParallelListComp #-}

module Sqrt where
import Data.List
import Numeric.AD

-- Find the kth order derivative
kdiff :: (a -> a) -> Int -> a
kdiff f 0 = diff f 0
kdiff f k = kdiff (diff f) (k-1)

-- Takes a function f and generates the Carleman matrix of f
carleman :: Floating a => (a -> a) -> Int ->  [[Float]]
carleman f size = [[ aux j k | k <- [0..size-1]] | j <- [0..size-1]]
    where
        aux :: Int -> Int -> Float
        aux j k = fromIntegral(kdiff (\x -> f x ^ j) k) / fromIntegral(product [1..k])

我不断收到此错误，我不明白如何使

f

与库兼容。我可以在 REPL 中执行类似

diff (\x -> x^2 + 1) 0

的操作，但不能在函数中执行。

[1 of 1] Compiling Sqrt             ( Sqrt.hs, interpreted )

Sqrt.hs:9:18: error:
    • Couldn't match type ‘a’
                     with ‘AD s (Numeric.AD.Internal.Forward.Forward a)’
      Expected: AD s (Numeric.AD.Internal.Forward.Forward a)
                -> AD s (Numeric.AD.Internal.Forward.Forward a)
        Actual: a -> a
      ‘a’ is a rigid type variable bound by
        the type signature for:
          kdiff :: forall a. (a -> a) -> Int -> a
        at Sqrt.hs:8:1-29
    • In the first argument of ‘diff’, namely ‘f’
      In the expression: diff f 0
      In an equation for ‘kdiff’: kdiff f 0 = diff f 0
    • Relevant bindings include
        f :: a -> a (bound at Sqrt.hs:9:7)
        kdiff :: (a -> a) -> Int -> a (bound at Sqrt.hs:9:1)
  |
9 | kdiff f 0 = diff f 0
  |                  ^

Sqrt.hs:10:25: error:
    • Couldn't match type ‘a’
                     with ‘AD s (Numeric.AD.Internal.Forward.Forward a)’
      Expected: AD s (Numeric.AD.Internal.Forward.Forward a)
                -> AD s (Numeric.AD.Internal.Forward.Forward a)
        Actual: a -> a
      ‘a’ is a rigid type variable bound by
        the type signature for:
          kdiff :: forall a. (a -> a) -> Int -> a
        at Sqrt.hs:8:1-29
    • In the first argument of ‘diff’, namely ‘f’
      In the first argument of ‘kdiff’, namely ‘(diff f)’
      In the expression: kdiff (diff f) (k - 1)
    • Relevant bindings include
        f :: a -> a (bound at Sqrt.hs:10:7)
        kdiff :: (a -> a) -> Int -> a (bound at Sqrt.hs:9:1)
   |
10 | kdiff f k = kdiff (diff f) (k-1)
   |                         ^
Failed, no modules loaded.

另一种方法会产生同样的错误，

-- Bell polynomial
bellY :: (Floating a, Enum a) => a -> a -> a
bellY x k = sum [stirling2 k i * x^i | i <- [0..k]]

-- Stirling numbers of the second kind
stirling2 :: (Eq a, Num a) => a -> a -> a
stirling2 0 0 = 1
stirling2 n k
  | n == 0 || k == 0 = 0
  | otherwise = k * stirling2 (n - 1) k + stirling2 (n - 1) (k - 1)

-- CarlemanMatrix function
carlemanMatrix :: (Floating a, Enum a) => (a -> a) -> a -> (Integer, Integer) -> [[a]]
carlemanMatrix f x0 (m, n) =
  map (\j -> map (\k -> if k == 0
                        then f x0 ^ j
                        else bellY (diff f x0) k / product [1..k])
                 [0..n])
      [0..m]

参考。 https://mathematica.stackexchange.com/a/836/72682

Answer 1

签名

carleman :: Floating a => (a -> a) -> Int ->  [[Float]]

...这是

的简写

carleman :: ∀ a . Floating a => (a -> a) -> Int ->  [[Float]]

表示

carleman

应该能够使用任何类型的函数。例如，我应该能够将它与

一起使用

f :: Double -> Double
f x = fromIntegral (denominator (toRational x))

这是非常不连续的，更不用说可微分了。

要能够区分函数，不能简单地是

Double

或

Float

。在这种情况下，正如错误消息所述，它需要对

AD s (Numeric.AD.Internal.Forward.Forward a)

值进行操作，这些值负责自动微分中涉及的额外簿记。实际上这还不够，因为你需要更高的导数；那里应该是

AD s (Tower a)

。并且

s

类型变量需要在函数内进行通用量化才能微分。 （这本质上确保您不会混淆不同函数的微分的无穷小值。）

所有这些繁琐的第 k 导数创建不仅是一个问题，因为每个导数都需要有不同的类型（每个订单需要不同的评估），而且效率非常低，因为你要分别区分每个订单而不是重复使用工作，而且 完全没有必要，因为您可以简单地一次性获取整个衍生品列表并使用它：

carleman :: ∀ a . Floating a => (∀ s . AD s (Tower a) -> AD s (Tower a)) -> Int -> [[a]]
carleman f size = [take size $ diffs (\x -> f x ^ j) 0] | j <- [0..size-1]]

此外，我想指出，根本没有充分的理由预先限制大小。您可以简单地创建一个无限矩阵，然后将其修剪为您可能需要的任何有限大小：

carleman :: ∀ a . Floating a => (∀ s . AD s (Tower a) -> AD s (Tower a)) -> [[a]]
carleman f = [diffs (\x -> f x ^ j) 0] | j <- [0..]]

如果您发现在签名中使用像

AD s (Tower a)

这样晦涩的类型很尴尬，另一种方法是将输入设为通用数值函数，即多态

Floating

函数，并让它在内部专门化为特定可微类型：

carleman :: ∀ a . Floating a => (∀ b . Floating b => b -> b) -> [[a]]

让函数与结果矩阵条目具有相同的类型是根本不可能的。这是不可能的，因为需要通过函数传递额外的信息才能计算这些条目。