Comment améliorer les performances de ce calcul numérique en Haskell?

Question

Je suis dans le milieu de portage de David Blei d'origine C mise en œuvre de Latent Dirichlet Allocation de Haskell, et je suis en train de décider de quitter le faible niveau de stuff en C. La fonction suivante est un exemple-c'est une approximation de la dérivée seconde de l' 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);
}

J'ai traduit cela en plus ou moins idiomatiques Haskell comme suit:

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)

Le problème est que quand je lance la fois par Critère, mon Haskell version est de six ou sept fois plus lent (je suis de la compilation avec -O2 sur le GHC 6.12.1). Certaines des fonctions similaires sont encore pire.

Je connais pratiquement rien au sujet de Haskell performance, et je ne suis pas vraiment intéressée à creuser par le biais de Cœur ou quelque chose comme ça, car je peux toujours appeler la poignée de math-intensif C fonctions par le biais de FFI.

Mais je suis curieux de savoir si il y a des branches basses de fruits que je suis en manque, une sorte d'extension ou de la bibliothèque ou de l'annotation que je pourrais utiliser pour accélérer le numérique trucs sans en faire trop laid.

---

Mise à JOUR: Voici les deux meilleures solutions, grâce à Don Stewart et Yitz. J'ai modifié Yitz de réponse légèrement à utiliser 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

Les performances des deux semble être presque exactement la même, avec l'un ou l'autre de gagner par un point de pourcentage ou deux selon les drapeaux du compilateur.

Comme camccann dit plus sur Reddit, la morale de l'histoire est "Pour de meilleurs résultats, utilisez Don Stewart comme votre GHC backend générateur de code." Sauf que la solution, le pari le plus sûr semble être simplement de traduire la C structures de contrôle directement dans Haskell, bien que la boucle de fusion peut donner des performances similaires à un plus idiomatique.

Je vais probablement jusqu'à la fin à l'aide de l' Data.Vector approche dans mon code.

Answer 1

Utiliser la même commande et les structures de données de rendement:

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

Je n'ai pas votre suite de tests, mais cela aboutit à la suite de l'asm:

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

Qui semble ok. C'est le genre de code de l' -fllvm backend fait un bon travail.

GCC se déroule la boucle si, et la seule façon de le faire est soit via le Modèle Haskell ou manuel dérouler. Vous pourriez envisager qu' (TH macro) si vous faites beaucoup de cela.

En fait, le GHC LLVM backend n'dérouler la boucle :-)

Enfin, si vous aimez vraiment l'original Haskell version, écrire à l'aide de flux de fusion combinators, et GHC va la convertir en boucles. (Exercice pour le lecteur).

Answer 2

Avant le travail d’optimisation, je ne dirais pas que votre traduction originale est le moyen le plus idiomatique d’exprimer en Haskell ce que fait le code C.

Comment se serait déroulé le processus d'optimisation si nous avions plutôt commencé avec les éléments suivants:

 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