类似的代码（对于指数加权偏差）在 Haskell 中比在 Python 中慢

Question

我在 python3 和 Haskell 中实现了 exponentially weighted moving average (ewma)（已编译）。它需要大约相同的时间。然而，当这个函数被应用两次时，haskell 版本会意外地变慢（超过 1000 次，而 python 版本只慢大约 2 倍）。

Python3版本：

import numpy as np
def ewma_f(y, tau):
    a = 1/tau
    avg = np.zeros_like(y)
    for i in range(1, len(y)):
        avg[i] = a*y[i-1]+(1-a)*avg[i-1]
    return avg

Haskell 列表：

ewmaL :: [Double] -> Double -> [Double]
ewmaL ys tau = reverse $ e (reverse ys) (1.0/tau)
    where e [x] a    = [a*x]
          e (x:xs) a = (a*x + (1-a)*(head $ e xs a) : e xs a)

Haskell 数组：

import qualified Data.Vector as V
ewmaV :: V.Vector Double -> Double -> V.Vector Double
ewmaV x tau = V.map f $ V.enumFromN 0 (V.length x)
    where
        f (-1) = 0
        f n    = (x V.! n)*a + (1-a)*(f (n-1))
        a = 1/tau

在所有情况下，计算时间大致相同（针对包含 10000 个元素的数组进行测试）。 Haskell 代码编译时没有任何标志，尽管 "ghc -O2" 没有任何区别。

我使用计算出的 ewma 来计算与这个 ewma 的绝对偏差；然后我将 ewma 函数应用于此偏差。

Python3:

def ewmd_f(y, tau):
    ewma = ewma_f(y, tau)
    return ewma_f(np.abs(y-ewma), tau)

它运行是 ewma 的两倍。

Haskell 列表：

ewmdL :: [Double] -> Double -> [Double]
ewmdL xs tau = ewmaL devs tau
    where devs = zipWith (\ x y -> abs $ x-y) xs avg
          avg  = (ewmaL xs tau)

Haskell 向量：

ewmdV :: V.Vector Double -> Double -> V.Vector Double
ewmdV xs tau = ewmaV devs tau
    where devs = V.zipWith (\ x y -> abs $ x-y) xs avg
          avg  = ewmaV xs tau

两个 ewmd 运行 > 1000 比他们的 ewma 同行慢。

我评估了 python3 代码：

from time import time
x = np.sin(np.arange(10000))
tau = 100.0

t1 = time()
ewma = ewma_f(x, tau)
t2 = time()
ewmd = ewmd_f(x, tau)
t3 = time()

print("EWMA took {} s".format(t2-t1))
print("EWMD took {} s".format(t3-t2))

我评估了 Haskell 代码：

import System.CPUTime

timeIt f = do
    start <- getCPUTime
    end   <- seq f getCPUTime
    let d = (fromIntegral (end - start)) / (10^12) in
        return (show d)

main = do
    let n = 10000 :: Int
    let tau = 100.0
    let l = map sin [0.0..(fromIntegral $ n-1)]
    let x = V.map sin $ V.enumFromN 0 n

    putStrLn "Vectors"
    aV <- timeIt $ V.last $ ewmaV x tau
    putStrLn $ "EWMA (vector) took "++aV
    dV <- timeIt $ V.last $ ewmdV x tau
    putStrLn $ "EWMD (vector) took "++dV
    putStrLn ""
    putStrLn "Lists"
    lV <- timeIt $ last $ ewmaL l tau
    putStrLn $ "EWMA (list) took "++lV
    lD <- timeIt $ last $ ewmdL l tau
    putStrLn $ "EWMD (list) took "++lD

Answer 1

您的 Python 和 Haskell 算法可能看起来相同，但它们实际上具有不同的渐近复杂度：

ewmaV x tau = V.map f $ V.enumFromN 0 (V.length x)
    where
        f (-1) = 0
        f n    = (x V.! n)*a + (1-a)
                     *(f (n-1)) -- Recursion!
        a = 1/tau

这使得 Haskell 实现 O (n²)，这是不可接受的。仅评估 V.last . ewmaV 时您没有注意到这一点的原因是懒惰：仅评估最后一个元素，您实际上并不需要处理整个向量，而是只得到一个 recursion-loop x。 OTOH，ewmdV 实际上会强制所有元素，因此会产生额外费用。

解决这个问题的一个简单（但我敢说不是最佳）方法是记住结果：

ewmaV :: V.Vector Double -> Double -> V.Vector Double
ewmaV x tau = result
  where result = V.map f $ V.enumFromN 0 (V.length x)
        f 0 = V.head x * a
        f n    = (x V.! n)*a + (1-a)*(result V.! (n-1))
        a = 1/tau

现在 ewmdV 是 ewmaV:

的两倍

sagemuej@sagemuej-X302LA:/tmp$ ghc wtmpf-file6122.hs -O2 && ./wtmpf-file6122
[1 of 1] Compiling Main             ( wtmpf-file6122.hs, wtmpf-file6122.o )
Linking wtmpf-file6122 ...
Vectors
EWMA (vector) took 4.932e-3
EWMD (vector) took 7.758e-3

（这些时间不是很可靠；要进行准确的性能测试，请使用 criterion。）

IMO 更好的解决方案是完全避免这种索引业务——我们不是在编写 Fortran，是吗？像 EWMA 这样的 IIR 最好以纯粹的“本地”方式实施；这可以在 Haskell 中用状态 monad 很好地表达，所以你独立于数据发送的容器。

import Data.Traversable
import Control.Monad (forM)
import Control.Monad.State

ewma :: Traversable t => t Double -> Double -> t Double
ewma x tau = (`evalState`0) . forM x $
         \xi -> state $ \carry
            -> let yi = a*xi + (1-a)*carry
               in (yi, yi)
 where a = 1/tau

虽然我们在概括：没有理由将其限制为仅适用于 Double 数据；您可以过滤任何 kind of variable that can be scaled and interpolated.

{-# LANGUAGE FlexibleContexts #-}
import Data.VectorSpace

ewma :: (Traversable t, VectorSpace v, Fractional (Scalar v))
             => t v -> Scalar v -> t v
ewma x tau = (`evalState`zeroV) . forM x $
         \xi -> state $ \carry
            -> let yi = a*^xi ^+^ (1-a)*^carry
               in (yi, yi)
 where a = 1/tau

这样，原则上您可以对 motion-blurring 存储在惰性流式传输的无限相框列表中的视频数据使用相同的过滤器，对于 lowpass-filtering 存储在未装箱中的无线电信号脉冲Vector。（VU.Vector 实际上没有 Traversable 实例；您需要替换 oforM。）

Answer 2

下面进行了两次递归调用：

ewmaL ys tau = reverse $ e (reverse ys) (1.0/tau)
    where e [x] a    = [a*x]
          e (x:xs) a = (a*x + (1-a)*(head $ e xs a) : e xs a)

我们可以进行一次递归调用，并将结果用于两种情况：

ewmaLcse :: [Double] -> Double -> [Double]
ewmaLcse ys tau = reverse $ e (reverse ys) (1.0/tau)
    where e [x] a    = [a*x]
          e (x:xs) a = (a*x + (1-a)*(head zs) : zs)
                       where zs = e xs a

我还选择了对列表的总和进行基准测试，以便强制计算所有列表的总和：

lV <- timeIt $ sum $ ewmaL l tau
putStrLn $ "EWMA (list) took "++lV
lVcse <- timeIt $ sum $ ewmaLcse l tau
putStrLn $ "EWMAcse (list) took "++lVcse

结果，n=10000

Lists
EWMA (list) took 2.384
EWMAcse (list) took 0.0

和n=20000

Lists
EWMA (list) took 16.472
EWMAcse (list) took 4.0e-3

顺便说一下，对于这个特定的递归，也可以使用标准库循环。这里我使出了mapAccumL：不需要double-reverse列表。

ewmaL2 :: [Double] -> Double -> [Double]
ewmaL2 ys tau = snd $ mapAccumL e 0 ys
    where a = 1/tau
          e !prevAvg !y = (nextAvg,nextAvg)
             where !nextAvg = a*y+(1-a)*prevAvg

（实际上，我使用 (nextAvg,nextAvg) 意味着更简单的 scanl 也可以完成这项工作。哦，好吧...）

Answer 3

总计，我最终使用 V.scanl:

解决了我的问题

ewma :: V.Vector Double -> Double -> V.Vector Double
ewma x tau = V.tail $ V.scanl ma 0 x
    where
        ma avg x = x*a + (1-a)*avg
        a = 1/tau

我的猜测是 scanl 具有 O(n) 复杂性，而不是 O(n²) 作为我的初始代码。对于未装箱的向量，它提供了相当不错的性能。

类似的代码（对于指数加权偏差）在 Haskell 中比在 Python 中慢

Similar code (for exponentially weighted deviation) is slower in Haskell than in Python

python

performance

haskell

weighted-average