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

Je suis en train de porter l’ implémentation C originale de Latent Dirichlet Allocation de David Blei sur Haskell, et j’essaie de décider s’il faut laisser une partie des choses de bas niveau dans C. La fonction suivante est un exemple – c’est une approximation de la dérivée seconde de lgamma :

 double sortinggamma(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); } 

Je l’ai traduit en Haskell plus ou moins idiomatique comme suit:

 sortinggamma :: Double -> Double sortinggamma 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 (\ab -> (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 lorsque je cours à la fois via Criterion , ma version Haskell est six ou sept fois plus lente (je comstack avec -O2 sur GHC 6.12.1). Certaines fonctions similaires sont encore pires.

Je ne connais pratiquement rien de la performance de Haskell, et je ne suis pas très intéressé par le fait de creuser dans Core ou quoi que ce soit de ce genre, car je peux toujours appeler la poignée de fonctions C intensives en mathématiques via FFI.

Mais je suis curieux de savoir s’il me manque des éléments faciles, comme une extension, une bibliothèque ou une annotation que je pourrais utiliser pour accélérer cette tâche numérique sans la rendre trop laide.


MISE À JOUR: Voici deux meilleures solutions, grâce à Don Stewart et Yitz . J’ai légèrement modifié la réponse de Data.Vector pour 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 sortinggamma_d :: Double -> Double sortinggamma_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) sortinggamma_y :: Double -> Double sortinggamma_y x = V.foldl' (+) (computeP $ x + 6) $ V.map invSq $ V.enumFromN x 6 

La performance des deux semble être presque identique, l’un ou l’autre gagnant un ou deux points de pourcentage selon les indicateurs du compilateur.

Comme le disait camccann chez Reddit , la morale de l’histoire est la suivante: «Pour de meilleurs résultats, utilisez Don Stewart comme générateur de code backend du GHC». Sauf cette solution, le pari le plus sûr semble être de traduire les structures de contrôle C directement dans Haskell, bien que la fusion en boucle puisse donner des performances similaires dans un style plus idiomatique.

Je vais probablement finir par utiliser l’approche Data.Vector dans mon code.

Utilisez les mêmes structures de contrôle et de données pour obtenir:

 {-# LANGUAGE BangPatterns #-} {-# OPTIONS_GHC -fvia-C -optc-O3 -fexcess-precision -optc-march=native #-} {-# INLINE sortinggamma #-} sortinggamma :: Double -> Double sortinggamma 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 donne les résultats suivants:

 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 

Ce qui a l’air correct C’est le genre de code que le backend -fllvm fait bien.

GCC déroule la boucle, et la seule façon de le faire est soit via Template Haskell, soit en mode manuel. Vous pourriez envisager cela (une macro TH) si vous en faites beaucoup.

En fait, le backend GHC LLVM déroule la boucle 🙂

Enfin, si vous aimez vraiment la version originale de Haskell, écrivez-la en utilisant des combinateurs de fusion de stream, et GHC le reconvertira en boucles. (Exercice pour le lecteur).

Avant le travail d’optimisation, je ne dirais pas que votre traduction originale est la manière la plus idiomatique d’exprimer dans Haskell ce que fait le code C.

Comment le processus d’optimisation se serait-il déroulé si nous avions commencé avec les éléments suivants:

 sortinggamma :: Double -> Double sortinggamma 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