2017-02-17 6 views
9
I 하스켈 python3과의 exponentially weighted moving average (EWMA)를 구현

보다 하스켈 느리다 (컴파일). 그것은 거의 같은 시간이 걸립니다. 이 함수를 두 번인가 그러나 때 (파이썬 버전은 약 2 배 느린 반면,보다 1000 배)을 하스켈 버전 다운 예기치 느려.유사한 코드 파이썬

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 

하스켈 : 어레이와

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) 

하스켈 모든 경우

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 요소). "GHC -O2는"어떤 차이를하지 않았지만 하스켈 코드는 어떤 플래그없이 컴파일.

이 계산 된 절대 값과 절대 편차를 계산하기 위해 계산 된 값을 사용했습니다. 그런 다음이 편차에 ewma 함수를 적용했습니다.

Python3은 :

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

이 두 번 이상 EWMA에 비해 실행됩니다. 목록과

하스켈 :

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

하스켈 벡터와는 :

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 

모두 자신의 EWMA에 비해> 1000 느리게 실행 ewmd. 내가와 하스켈 코드 평가

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)) 

:

나는 함께 python3 코드를 평가 동등한 보일 수 있습니다

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 
+1

다음은 해당 기준 파일의 모양입니다. http://sprunge.us/VONE – Michael

답변

12

파이썬과 하스켈 알고리즘을하지만, 실제로 다른 점근 적 복잡성이 :

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 구현을 (n2)으로 만든다. 이는 받아 들일 수 없다. V.last . ewmaV 만 평가할 때 이것을 알지 못하는 이유는 마지막 요소 만 평가하기 때문에 실제로 전체 벡터를 처리 할 필요가 없으며 x에 하나의 순환 루프 만 가져옵니다. OTOH, ewmdV은 실제로 모든 요소를 ​​강제하므로 추가 비용이 발생합니다.

이 주위에 얻을 수있는 방법은 결과 memoise하는 하나의 간단한 (하지만 최적하지, 나는 아마 ... 일 것)

:

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 

지금 ewmdVewmaV만큼 ≈twice 취

[email protected]:/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는이 색인 생성 사업을 완전히 피할 것입니다. 우리는 포트란을 쓰지 않을 것입니까? EWMA와 같은 IIR은 순전히 "로컬"방식으로 더 잘 구현됩니다. 이것은 당신이 무엇을 컨테이너 데이터 선박의 독립적 인 것 때문에, 상태 모나드와 하스켈에서 잘 표현 될 수

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 

이렇게하면 원리에 저장된 모션 블러 링 된 비디오 데이터에 대해 동일한 필터를 사용할 수 느리게 저역 통과 필터링 박싱 Vector에 저장된 무선 신호 펄스에 관해서는, 화상 프레임의 무한리스트 스트리밍. (VU.Vector 실제로 더 Traversable 인스턴스가 없습니다, 당신은 다음 oforM을 대체해야합니다.)

+0

Monad.State 및 VectorSpace에 대한 훌륭한 답변과 포인터를 보내 주셔서 감사합니다. 결국 V.scanl 및 unboxed vector를 사용하여 ewma를 구현하는 합리적인 속도를 달성했습니다. –

4

다음은 두 개의 재귀 호출한다 :

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 

는 또한 그래서 모든 계산에 강제로 벤치 마크리스트의 합을 선택했다 그런데 n=20000

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

n=10000

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

와개

결과, 하나는이 특정 재귀 표준 라이브러리 루프를 사용할 수 있습니다. 여기 나는 mapAccumL에 의지했다 : 목록을 이중으로 뒤집을 필요가 없다.

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 

3

leftaroundabout's great answer을 추가하기 위해 (사실, 난 (nextAvg,nextAvg)를 사용하고 있다는 사실은 간단한 scanl도 ... 작업. 잘 아을 할 것이라는 점을 의미한다), 나는 결국 해결 한 내 문제는 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 2)로 내 초기 코드를 가지고있다. 박스가없는 벡터를 사용하면 꽤 괜찮은 성능을 얻을 수 있습니다.