我正在移植 David Blei 的原版 C implementation将 Latent Dirichlet Allocation 分配给 Haskell,我正在尝试决定是否在 C 中保留一些低级内容。以下函数是一个示例 - 它是 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 版本慢了六到七倍(我在 GHC 6.12.1 上用 -O2
编译)。一些类似的功能更糟糕。
我对 Haskell 的性能几乎一无所知,而且我对 digging through Core 也不是很感兴趣或类似的东西,因为我总是可以通过 FFI 调用少数数学密集型 C 函数。
但我很好奇是否有我遗漏的唾手可得的成果——某种扩展、库或注释,我可以用来加速这些数字的东西,而不会使它变得太难看。
更新:这里有两个更好的解决方案,感谢 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
两者的性能似乎几乎完全相同,根据编译器标志,一个或另一个获胜一个百分点或两个百分点。
作为camccann说over at Reddit ,故事的寓意是“为了获得最佳结果,请使用 Don Stewart 作为您的 GHC 后端代码生成器。”除了该解决方案,最安全的选择似乎只是将 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 会展开循环,但唯一的方法是通过 Template Haskell 或手动展开。如果经常这样做,您可能会考虑(一个 TH 宏)。
实际上,GHC LLVM 后端确实展开了循环:-)
最后,如果你真的喜欢原始的 Haskell 版本,请使用 stream fusion combinators, 编写它GHC 会将其转换回循环。 (读者练习)。
关于c - 如何提高 Haskell 中这种数值计算的性能?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/2978979/