Haskell에서 이 수치 계산의 성능을 향상시키는 방법은 무엇입니까?

Haskell에서 이 수치 계산의 성능을 향상시키는 방법은 무엇입니까?

David Blei의 잠복 디리클레 할당의 원래 C 구현을 Haskell에 이식하는 중인데, C에 하위 수준의 것들을 남겨둘지 결정하려고 합니다.다음 함수는 한 가지 예입니다. 이것은 의 2차 도함수의 근사치입니다.lgamma:

double trigamma(double x)
{
    double p;
    int i;

    x=x+6;
    p=1/(x*x);
    p=(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)
         *p-0.033333333333333)*p+0.166666666666667)*p+1)/x+0.5*p;
    for (i=0; i<6 ;i++)
    {
        x=x-1;
        p=1/(x*x)+p;
    }
    return(p);
}

나는 이것을 다소 관용적인 Haskell로 다음과 같이 번역했다.

trigamma :: Double -> Double
trigamma x = snd $ last $ take 7 $ iterate next (x' - 1, p')
  where
    x' = x + 6
    p  = 1 / x' ^ 2
    p' = p / 2 + c / x'
    c  = foldr1 (\a b -> (a + b * p)) [1, 1/6, -1/30, 1/42, -1/30, 5/66]
    next (x, p) = (x - 1, 1 / x ^ 2 + p)

문제는 Criterion을 통해 둘 다 실행하면 Haskell 버전이 6~7배 느리다는 것입니다(컴파일 중).-O2(GHC 6.12.1)에 있습니다.유사한 기능은 더 나쁘다.

저는 Haskell의 성능에 대해 아는 것이 거의 없고, FFI를 통해 수학 집약적인 C 함수를 언제든지 호출할 수 있기 때문에 Core를 파헤치거나 하는 것에는 그다지 관심이 없습니다.

하지만 제가 놓치고 있는 작은 과일이 있는지 궁금합니다. 너무 추하게 만들지 않고 숫자를 빠르게 표시하기 위해 사용할 수 있는 확장 기능이나 라이브러리, 주석 같은 것이 있는지 궁금해요.

---

업데이트: Don Stewart와 Yitz 덕분에 두 가지 더 나은 솔루션이 있습니다.Yitz의 답변을 약간 수정해서Data.Vector.

invSq x = 1 / (x * x)
computeP x = (((((5/66*p-1/30)*p+1/42)*p-1/30)*p+1/6)*p+1)/x+0.5*p
  where p = invSq x

trigamma_d :: Double -> Double
trigamma_d x = go 0 (x + 5) $ computeP $ x + 6
  where
    go :: Int -> Double -> Double -> Double
    go !i !x !p
        | i >= 6    = p
        | otherwise = go (i+1) (x-1) (1 / (x*x) + p)

trigamma_y :: Double -> Double
trigamma_y x = V.foldl' (+) (computeP $ x + 6) $ V.map invSq $ V.enumFromN x 6

컴파일러 플래그에 따라 어느 한쪽이 퍼센트 포인트 또는2 포인트 차이로 승리하는 등, 양쪽의 퍼포먼스는 거의 같은 것 같습니다.

Reddit에서 camccann이 말했듯이 이 이야기의 교훈은 "최상의 결과를 얻으려면 GHC 백엔드 코드 생성기로 Don Stewart를 사용하라"입니다.이 솔루션을 제외한다면 가장 안전한 방법은 C 제어 구조를 Haskell로 직접 변환하는 것입니다.단, 루프 퓨전에서는 보다 관용적인 스타일로 유사한 퍼포먼스를 제공할 수 있습니다.

저는 아마 이 머신을 사용하게 될 겁니다.Data.Vector내 코드로 접근합니다.

동일한 제어 및 데이터 구조를 사용하여 다음을 산출합니다.

{-# LANGUAGE BangPatterns #-}
{-# OPTIONS_GHC -fvia-C -optc-O3 -fexcess-precision -optc-march=native #-}

{-# INLINE trigamma #-}
trigamma :: Double -> Double
trigamma x = go 0 (x' - 1) p'
    where
        x' = x + 6
        p  = 1 / (x' * x')

        p' =(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)
                  *p-0.033333333333333)*p+0.166666666666667)*p+1)/x'+0.5*p

        go :: Int -> Double -> Double -> Double
        go !i !x !p
            | i >= 6    = p
            | otherwise = go (i+1) (x-1) (1 / (x*x) + p)

테스트 슈트는 가지고 있지 않습니다만, 다음과 같은 결과가 됩니다.

A_zdwgo_info:
        cmpq    $5, %r14
        jg      .L3
        movsd   .LC0(%rip), %xmm7
        movapd  %xmm5, %xmm8
        movapd  %xmm7, %xmm9
        mulsd   %xmm5, %xmm8
        leaq    1(%r14), %r14
        divsd   %xmm8, %xmm9
        subsd   %xmm7, %xmm5
        addsd   %xmm9, %xmm6
        jmp     A_zdwgo_info

괜찮은 것 같네요.이건 코드 종류야-fllvm백엔드는 잘 작동합니다.

단, GCC는 루프를 언롤합니다.이러한 방법은 템플릿 해스켈 또는 수동 언롤 중 하나입니다.이 작업을 많이 수행하는 경우 (TH 매크로)를 고려할 수 있습니다.

실제로 GHC LLVM 백엔드는 루프를 언롤합니다:-)

마지막으로 원본 Haskell 버전이 마음에 드신다면 스트림 퓨전 콤비네이터를 사용하여 작성하시면 GHC가 다시 루프로 변환됩니다.(독자 연습)

최적화 작업 전에, 당신의 원번역이 Haskell에서 C코드가 무엇을 하고 있는지를 표현하는 가장 관용적인 방법이라고는 말할 수 없습니다.

대신 다음과 같이 시작한다면 최적화 프로세스가 어떻게 진행되었을까요?

trigamma :: Double -> Double
trigamma x = foldl' (+) p' . map invSq . take 6 . iterate (+ 1) $ x
where
  invSq y = 1 / (y * y)
  x' = x + 6
  p  = invSq x'
  p' =(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)
              *p-0.033333333333333)*p+0.166666666666667)*p+1)/x'+0.5*p

언급URL : https://stackoverflow.com/questions/2978979/how-to-improve-performance-of-this-numerical-computation-in-haskell