Clojure 中惰性卷积 fn 的问题

Question

我正在编写一些信号处理软件，我首先编写一个离散卷积函数< /a>。

这对于前一万个左右的值列表来说效果很好，但随着它们变得更大（比如 100k），我当然开始收到 StackOverflow 错误。

不幸的是，我在将命令式卷积算法转换为递归和递归算法时遇到了很多麻烦。惰性版本实际上足够快，可以使用（至少有一点优雅也很好）。

我也不能 100% 确定我的这个功能完全正确，但如果我遗漏了什么/做错了什么，请告诉我。我认为这是正确的。

(defn convolve
  "
    Convolves xs with is.

    This is a discrete convolution.

    'xs  :: list of numbers
    'is  :: list of numbers
  "
  [xs is]
  (loop [xs xs finalacc () acc ()]
    (if (empty? xs)
      (concat finalacc acc)
      (recur (rest xs)
             (if (empty? acc)
               ()
               (concat finalacc [(first acc)]))
             (if (empty? acc)
               (map #(* (first xs) %) is)
               (vec-add
                (map #(* (first xs) %) is)
                (rest acc)))))))

我非常感谢任何类型的帮助：我仍然在 Clojure 中了解自己的方位，并且使这种优雅、懒惰和/或递归会很棒。

我有点惊讶的是，在 Lisp 中用命令式语言表达相当的算法是多么困难。但也许我做错了！

编辑： 只是为了展示用命令式语言表达是多么容易，并为人们提供运行良好且易于阅读的算法，这里是 Python 版本。除了更短、更简洁、更容易推理之外，它的执行速度比 Clojure 代码快几个数量级：甚至是我使用 Java 数组的命令式 Clojure 代码。

from itertools import repeat

def convolve(ns, ms):
    y = [i for i in repeat(0, len(ns)+len(ms)-1)]
    for n in range(len(ns)):
        for m in range(len(ms)):
            y[n+m] = y[n+m] + ns[n]*ms[m]
    return y

另一方面，这里是命令式 Clojure 代码。它还会从卷积中删除最后一个非完全沉浸的值。因此，除了缓慢且丑陋之外，它并不是 100% 实用。也不实用。

(defn imp-convolve-1
  [xs is]
  (let [ys (into-array Double (repeat (dec (+ (count xs) (count is))) 0.0))
        xs (vec xs)
        is (vec is)]
     (map #(first %)
          (for [i (range (count xs))]
            (for [j (range (count is))]
              (aset ys (+ i j)
                    (+ (* (nth xs i) (nth is j))
                       (nth ys (+ i j)))))))))

这太令人沮丧了。请有人告诉我我刚刚错过了一些明显的事情。

编辑3：

这是我昨天想到的另一个版本，展示了我希望如何表达它（尽管其他解决方案非常优雅；我只是把另一个解决方案放在那里！）

(defn convolve-2
  [xs is]
  (reduce #(vec-add %1 (pad-l %2 (inc (count %1))))
         (for [x xs]
           (for [i is]
             (* x i)))))

它使用这个实用函数vec-add：

(defn vec-add
  ([xs] (vec-add xs []))
  ([xs ys]
     (let [lxs (count xs)
           lys (count ys)
           xs (pad-r xs lys)
           ys (pad-r ys lxs)]
       (vec (map #(+ %1 %2) xs ys))))
  ([xs ys & more]
     (vec (reduce vec-add (vec-add xs ys) more))))
     (vec (reduce vec-add (vec-add xs ys) more))))

Answer 1

(defn ^{:static true} convolve ^doubles [^doubles xs ^doubles is]
  (let [xlen (count xs)
        ilen (count is)
        ys   (double-array (dec (+ xlen ilen)))]
    (dotimes [p xlen]
      (dotimes [q ilen]
        (let [n (+ p q), x (aget xs p), i (aget is q), y (aget ys n)]
          (aset ys n (+ (* x i) y)))))
    ys))

如果我在 Clojure equiv 分支中执行此操作，则重复 jg-faustus 的版本。对我有用。在 i7 Mackbook Pro 上，1,000,000 个点约为 400 毫秒，100,000 个点约为 25 毫秒。

Answer 2

堆栈溢出错误的可能原因是惰性重击变得太深。 （concat 和 map 是惰性的）。尝试将这些调用包装在 doall 中以强制评估其返回值。

至于更实用的解决方案，请尝试这样的方法：

(defn circular-convolve
  "Perform a circular convolution of vectors f and g"
  [f g]
  (letfn [(point-mul [m n]
        (* (f m) (g (mod (- n m) (count g)))))
      (value-at [n]
        (reduce + (map #(point-mul % n) (range (count g)))))]
    (map value-at (range (count g)))))

使用可以使用reduce轻松执行求和，并且由于map产生惰性序列，因此该函数也是惰性的。

Answer 3

无法帮助实现高性能函数式版本，但通过放弃惰性并添加类型提示，您可以使命令式版本获得 100 倍的加速：

(defn imp-convolve-2 [xs is]
  (let [^doubles xs (into-array Double/TYPE xs)
        ^doubles is (into-array Double/TYPE is)
        ys (double-array (dec (+ (count xs) (count is)))) ]
    (dotimes [i (alength xs)]
      (dotimes [j (alength is)]
        (aset ys (+ i j)
          (+ (* (aget xs i) (aget is j))
            (aget ys (+ i j))))))
    ys))

使用 xs 大小 100k 且 is< /code> size 2，当你的 imp-convolve-1 包裹在 doall 中时，你的 imp-convolve-1 在我的机器上需要大约 6,000 毫秒，而这个需要大约 35 毫秒。

更新

这是一个惰性函数版本：

(defn convolve 
  ([xs is] (convolve xs is []))
  ([xs is parts]
    (cond
      (and (empty? xs) (empty? parts)) nil
      (empty? xs) (cons
                    (reduce + (map first parts))
                    (convolve xs is
                      (remove empty? (map rest parts))))
      :else (cons
              (+ (* (first xs) (first is))
                (reduce + (map first parts)))
              (lazy-seq
                (convolve (rest xs) is
                  (cons 
                    (map (partial * (first xs)) (rest is))
                    (remove empty? (map rest parts)))))))))

在大小 100k 和 2 上，它的时钟时间约为 600 毫秒（450-750 毫秒不等），而 imp-convolve-1 约为 6,000 毫秒，imp-convolve 约为 35 毫秒-2。

所以它是实用的、惰性的并且具有可接受的性能。不过，它的代码量是命令式版本的两倍，而且我又花了 1-2 个小时才找到，所以我不确定我是否明白这一点。

当纯函数使代码更短或更简单，或者比命令式版本有其他一些好处时，我完全支持纯函数。当他们不这样做时，我不反对切换到命令式模式。

这是我认为 Clojure 很棒的原因之一，因为您可以根据需要使用任何一种方法。

更新2：

我会修改我的“功能上这样做的意义是什么”，说我喜欢此功能实现（第二个，页面下方），作者：David Cabana。

它简短、可读，时间约为 140 毫秒，输入大小与上面相同（100,000 和 2），使其成为迄今为止我尝试过的性能最佳的功能实现。

考虑到它是函数式的（但不是惰性的），不使用类型提示并且适用于所有数字类型（不仅仅是双精度数），这是相当令人印象深刻的。

Answer 4

(defn convolve [xs is]
  (if (> (count xs) (count is))
    (convolve is xs)
    (let [os (dec (+ (count xs) (count is)))
          lxs (count xs)
          lis (count is)]
      (for [i (range os)]
        (let [[start-x end-x] [(- lxs (min lxs (- os i))) (min i (dec lxs))]
              [start-i end-i] [(- lis (min lis (- os i))) (min i (dec lis))]]
          (reduce + (map *
                         (rseq (subvec xs start-x (inc end-x)))
                         (subvec is start-i (inc end-i)))))))))

可以用简洁的术语表达一个惰性的、函数式的解决方案。唉，> 的表现2k不切实际。我很想知道是否有方法可以在不牺牲可读性的情况下加快速度。

编辑：
阅读 drcabana 关于该主题的信息帖子（http： //erl.nfshost.com/2010/07/17/discrete-convolution-of-finite-vectors/），我已经更新了我的代码以接受不同大小的向量。他的实现表现更好：
对于 xs 尺寸 3，尺寸为 1000000，~2s 与 ~3s

编辑 2：
采用 drcabana 简单反转 xs 和 padding 的想法，我得出的结论是：

(defn convolve [xs is]
  (if (> (count xs) (count is))
    (convolve is xs)
    (let [is (concat (repeat (dec (count xs)) 0) is)]
      (for [s (take-while not-empty (iterate rest is))]
         (reduce + (map * (rseq xs) s))))))

这可能是尽可能简洁的，但总体来说仍然较慢，可能是由于需要一段时间。感谢博客作者经过深思熟虑的方法。这里唯一的优点是上面的代码确实很懒，因为如果我询问 (nth res 10000)，它只需要前 10k 次计算即可得出结果。

Answer 5

并不是对您提出的众多问题中的任何一个的真正答案，但我对您没有提出的问题有一些评论。

1. 您可能不应该在向量上使用nth。是的，它的时间复杂度为 O(1)，但由于 nth 在 O(n) 中适用于其他序列，因此 (a) 并没有明确表明您期望输入是向量，并且 ( b) 意味着如果你犯了一个错误，你的代码会神秘地变得非常慢，而不是立即失败。

2. for 和 map 都是惰性的，而 aset 仅具有副作用。这种组合会导致灾难：对于类似 for 的副作用，请使用 doseq。

3. for 和 doseq 允许多个绑定，因此您不需要像在 Python 中那样（显然）堆积它们的负载。

   (doseq [i（范围（计数 cs））
           j(范围(计数为))] 
     ...)
   

   会做你想做的事。

4. #(first %) 更简洁地写为 first；同样，#(+ %1 %2) 是 +。

5. 对一堆不需要向量的中间结果调用 vec 会降低速度。具体来说，在 vec-add 中，仅在创建最终返回值时调用 vec 就足够了：在 (vec (reduce foo bar)) 中，有如果 foo 从未将其中间结果用于随机访问，则它没有理由将其转换为向量。