优化简单的Common Lisp gibbs采样程序

Question

作为练习,我在Darren Wilkinson的各种语言(重新访问)的博客文章Gibbs采样器中重写了示例程序.

代码如下所示.这个代码在我的(5岁)机器上运行大约53秒,使用SBCL 1.0.56,使用buildapp创建核心映像,然后运行它

time ./gibbs > gibbs.dat

由于这是如何计算帖子中其他语言的时间,我想我会做一些类似的事情.帖子中的C代码在大约25秒内运行.如果可能的话,我想尝试加速Lisp代码.

##############################
gibbs.lisp
##############################
(eval-when (:compile-toplevel :load-toplevel :execute)
  (require :cl-rmath) (setf *read-default-float-format* 'double-float))

(defun gibbs (N thin)
  (declare (fixnum N thin))
  (declare (optimize (speed 3) (safety 1)))
  (let ((x 0.0) (y 0.0))
    (declare (double-float x y))
    (print "Iter x y")
    (dotimes (i N)
      (dotimes (j thin)
    (declare (fixnum i j))
    (setf x (cl-rmath::rgamma 3.0 (/ 1.0 (+ (* y y) 4))))
    (setf y (cl-rmath::rnorm (/ 1.0 (+ x 1.0)) (/ 1.0 (sqrt (+ (* 2 x) 2))))))
      (format t "~a ~a ~a~%" i x y))))

(defun main (argv)
  (declare (ignore argv))
  (gibbs 50000 1000))

然后我用as gibbs调用构建了可执行文件sh gibbs.shgibbs.sh

##################
gibbs.sh
##################
buildapp --output gibbs --asdf-tree /usr/share/common-lisp/source/ --asdf-tree /usr/local/share/common-lisp/source/ --load-system cl-rmath --load gibbs.lisp --entry main

在使用SBCL 1.0.56进行编译时,我得到了6个编译器注释,它在下面重现.我不知道该怎么办,但对任何提示都会感激不尽.

; compiling file "/home/faheem/lisp/gibbs.lisp" (written 30 MAY 2012 02:00:55 PM):

; file: /home/faheem/lisp/gibbs.lisp
; in: DEFUN GIBBS
;     (SQRT (+ (* 2 X) 2))
; 
; note: unable to
;   optimize
; due to type uncertainty:
;   The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT))
;                           &OPTIONAL), not a (VALUES FLOAT &REST T).

;     (/ 1.0d0 (SQRT (+ (* 2 X) 2)))
; 
; note: unable to
;   optimize
; due to type uncertainty:
;   The second argument is a (OR (DOUBLE-FLOAT 0.0)
;                                (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
;                                                                DOUBLE-FLOAT).
; 
; note: forced to do static-fun Two-arg-/ (cost 53)
;       unable to do inline float arithmetic (cost 12) because:
;       The second argument is a (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
;       The result is a (VALUES (OR (COMPLEX DOUBLE-FLOAT) (DOUBLE-FLOAT 0.0))
;                               &OPTIONAL), not a (VALUES DOUBLE-FLOAT &REST T).

;     (CL-RMATH:RGAMMA 3.0d0 (/ 1.0d0 (+ (* Y Y) 4)))
; 
; note: doing float to pointer coercion (cost 13)

;     (SQRT (+ (* 2 X) 2))
; 
; note: doing float to pointer coercion (cost 13)

;     (CL-RMATH:RNORM (/ 1.0d0 (+ X 1.0d0)) (/ 1.0d0 (SQRT (+ (* 2 X) 2))))
; 
; note: doing float to pointer coercion (cost 13)
; 
; compilation unit finished
;   printed 6 notes

; /home/faheem/lisp/gibbs.fasl written
; compilation finished in 0:00:00.073

更新1:Rainer Joswig的回答指出SQRT的论点可能是否定的,这是我看到的模糊编译器笔记的来源,即

The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT))
    ;                           &OPTIONAL), not a (VALUES FLOAT &REST T).

编译器抱怨说,由于它不知道参数的值是否为正,因此结果可能是一个复数.由于在该示例中,该值x是与gamma分布的样本变量,因此总是大于0. Stephan在SBCL用户邮件列表中有人指出了这一点(请参阅"优化简单公共Lisp"主题中的第二条消息)Gibbs采样器程序",这可以通过声明x大于或等于如下来解决,

(declare (type (double-float 0.0 *) x))

有关FLOAT类型 和 间隔指示符的相关文档,请参阅Common Lisp Hyperspec .

这似乎加快了代码的速度.它现在可靠地低于52秒,但仍然没有太大的收益.这也留下了笔记

note: doing float to pointer coercion (cost 13)

如果出于某种原因这是不可修复的,我想知道原因.此外,无论如何,对说明意味着什么的解释都会很有趣.具体来说,这个词pointer在这里意味着什么？这与C函数被调用的事实有关吗？此外,成本13似乎没有用.测量什么？

此外,Rainer建议可以减少耗尽,这可能会减少运行时间.我不知道是否可以减少减少,或者是否会减少运行时间,但我对意见和方法感兴趣.总的来说,似乎没有太多可以改善这个功能的性能.也许它太小而且简单.

Answer 1

请注意，Common Lisp 有一个特殊的运算符。它允许您声明表达式结果的类型。例如，如果可能的话，这可以让您缩小类型范围。

例如，结果是什么(SQRT somefloat)？它可以是浮点数，但如果为负数，则可以是复数somefloat。如果您知道 somefloat 始终为正（并且只有在那时），那么您可以编写(the double-float (sqrt somefloat)). 然后编译器可能能够生成更有效的代码。

另请注意，Common Lisp 有OPTIMIZE声明。如果您想要最快的代码，您需要确保相应地设置它们。可能仅用于个人功能。通常，这比非常激进地全局更改优化要好。

Common Lisp 有一个函数DISASSEMBLE可以让你查看编译后的代码。

然后是宏TIME。您从中获得的有趣信息包括它的作用有多大。对于双浮点算术，可能需要大量的操作。在 SBCL 邮件列表上寻求帮助会很有用。也许有人可以告诉你如何避免这种欺骗。