優化簡單Common Lisp gibbs採樣器程序

Question

作爲練習，我重寫了Darren Wilkinson在博客文章Gibbs sampler in various languages (revisited)中的示例程序。

代碼如下所示。此代碼我（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)) 

然後我生成的可執行​​與調用sh gibbs.sh與gibbs.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 

我得到6只編譯器注意到SBCL 1.0.56，其下面重現編譯時。我不知道該怎麼做，但會對任何提示感激。

; 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's answer指出，SQRT的說法可能是負的，這 是模糊的編譯器的筆記我看到的來源，即

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

的編譯器抱怨因爲它不知道參數的值是否爲正， 結果可能是一個複數。因爲在這個例子中，值x是伽瑪分佈的一個樣本變量，它總是大於0.它在SBCL用戶郵件列表 （見線程"Optimizing a simple Common Lisp Gibbs sampler program"中的第二條消息，由Stephan幫助指出）這可以通過聲明x到大於或零以下來解決，

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

見Common Lisp的Hyperspec的相關文件約 FLOAT types 和 Interval Designators。

氏這似乎加快了代碼的速度。現在它可靠地低於52秒，但仍然沒有太大的收益。 這也留下要注意的事項

note: doing float to pointer coercion (cost 13) 

如果這是不是可以解決由於某些原因，我想知道這是爲什麼。此外，無論如何，說明該筆記意味着什麼會很有趣。 具體來說，這個詞pointer這個詞是什麼意思？這與C函數被調用的事實有關嗎？另外，費用13似乎不是 非常有用。什麼是測量？

另外，Rainer建議可以減少consing，這可能會減少運行時間。 我不知道是否可能減少收縮，或者是否會減少運行時間，但我會對意見和方法感興趣。總的來說，似乎沒有什麼可以改進這個功能的性能。 也許它太小，簡單。

Answer 1

請注意，Common Lisp有一個THE特殊運算符。它允許您爲表達式結果聲明類型。例如，這可以讓你儘可能縮小類型。

例如(SQRT somefloat)的結果是什麼？它可以是一個浮點數，但如果somefloat爲負數，它可能是一個複數。如果你知道somefloat總是正面的（並且只有這樣），那麼你可以寫(the double-float (sqrt somefloat))。編譯器可能會生成更高效的代碼。

另請注意，Common Lisp有OPTIMIZE聲明。如果你想要最快的代碼，你需要確保你相應地設置它們。可能僅限於個人功能。通常情況下，它比改變全球優化要好得多。

Common Lisp有一個功能DISASSEMBLE它可以讓你看看編譯的代碼。

然後是宏TIME。您從中獲得的有趣信息包括它的作用。使用雙浮點算術可能會引起大量的考慮。在SBCL郵件列表上尋求幫助會很有幫助。也許有人可以告訴你如何避免這種顧慮。

Answer 2

這個工作對我來說：

(sqrt (the (double-float 0d0) (+ (* 2d0 x) 2d0)))