Clojure/Incanter中的快速矢量数学

Question

我目前正在研究Clojure和Incanter作为R的替代品.(不是我不喜欢R,但尝试使用新语言会很有趣.)我喜欢Incanter并且发现语法很吸引人,但是矢量化操作比较慢例如,R或Python.

作为一个例子,我想使用Incanter向量运算,Clojure map和R来获得向量的第一阶差分.以下是所有版本的代码和时间.如你所见,R显然更快.

Incanter和Clojure:

(use '(incanter core stats)) 
(def x (doall (sample-normal 1e7))) 
(time (def y (doall (minus (rest x) (butlast x))))) 
"Elapsed time: 16481.337 msecs" 
(time (def y (doall (map - (rest x) (butlast x))))) 
"Elapsed time: 16457.850 msecs"

R:

rdiff <- function(x){ 
   n = length(x) 
   x[2:n] - x[1:(n-1)]} 
x = rnorm(1e7) 
system.time(rdiff(x)) 
   user  system elapsed 
  1.504   0.900   2.561

所以我想知道有没有办法加速Incanter/Clojure中的矢量操作？还欢迎涉及使用来自Clojure的循环,Java数组和/或库的解决方案.

我还向Incanter Google小组发布了此问题,目前尚无回复.

更新:我已将Jouni的答案标记为已接受,请参阅下面的我自己的答案,我已经清理了他的代码并添加了一些基准测试.

Answer 1

我的最终解决方案

在所有测试之后,我发现两种略有不同的方式以足够的速度进行计算.

首先我使用了diff具有不同类型返回值的函数,下面是返回向量的代码,但我还定时返回一个双数组的版本(用y替换(vec y))和Incanter.matrix(替换( vec y)与矩阵y).此函数仅基于java数组.这是基于Jouni的代码,删除了一些额外的类型提示.

另一种方法是使用Java数组进行计算,并将值存储在瞬态向量中.从时间上看,如果你不想返回和数组的功能,这比接近1快一点.这是在功能上实现的difft.

因此,选择实际上取决于您不想对数据做什么.我想一个好的选择是重载函数,以便它返回调用中使用的相同类型.实际上将java数组传递给diff而不是向量会使~1s更快.

不同功能的时间:

差异返回矢量:

(time (def y (diff x)))
"Elapsed time: 4733.259 msecs"

diff returns Incanter.matrix:

(time (def y (diff x)))
"Elapsed time: 2599.728 msecs"

diff返回双数组:

(time (def y (diff x)))
"Elapsed time: 1638.548 msecs"

difft:

(time (def y (difft x)))
"Elapsed time: 3683.237 msecs"

功能

(use 'incanter.stats)
(def x (vec (sample-normal 1e7)))

(defn diff [x]
  (let [y (double-array (dec (count x)))
        x (double-array x)] 
   (dotimes [i (dec (count x))]
     (aset y i
       (- (aget x (inc i))
                   (aget x i))))
   (vec y)))


(defn difft [x]
  (let [y (vector (range n))
        y (transient y)
        x (double-array x)]
   (dotimes [i (dec (count x))]
     (assoc! y i
       (- (aget x (inc i))
                   (aget x i))))
   (persistent! y))) 

Answer 2

这是我的系统上的Java数组实现比R代码(YMMV)更快.注意启用反射警告,这在优化性能时是必不可少的,并且y上的重复类型提示(def上的那个似乎对aset没有帮助)并将所有内容都转换为原始double值(dotimes确保我是一个原始的int).

(set! *warn-on-reflection* true)
(use 'incanter.stats)
(def ^"[D" x (double-array (sample-normal 1e7)))

(time
 (do
   (def ^"[D" y (double-array (dec (count x))))
   (dotimes [i (dec (count x))]
     (aset ^"[D" y
       i
       (double (- (double (aget x (inc i)))
                  (double (aget x i))))))))