7

我正在编写一些信号处理软件,并从写出离散卷积函数开始。

这适用于前一万左右的值列表,但是随着它们变大(例如,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))))
4

5 回答 5

4
(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 毫秒。

于 2010-07-17T16:38:07.207 回答
4

堆栈溢出错误的可能原因是惰性 thunk 变得太深。(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产生一个惰性序列,这个函数也是惰性的。

于 2010-07-16T01:20:47.583 回答
4

高性能的函数式版本无能为力,但您可以通过消除惰性并添加类型提示来获得命令式版本的 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大小为 2 的情况下,你的 imp-convolve-1 在我的机器上包裹在一个 doall 中需要大约 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-2 约为 35 毫秒。

所以它是功能性的、惰性的并且具有可容忍的性能。尽管如此,它的代码量还是命令式版本的两倍,而且我花了 1-2 个小时才找到,所以我不确定我是否明白这一点。

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

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

更新 2:

我将通过说我喜欢David Cabana 的这个功能实现(第二个,在页面下方)来修改我的“在功能上做这个有什么意义”。

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

考虑到它是功能性的(但不是惰性的),不使用类型提示并且适用于所有数字类型(不仅仅是双精度数),这令人印象深刻。

于 2010-07-16T20:30:30.283 回答
3
(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 vs ~3s

编辑 2: 采用 drcabana 的简单反转 xs 和填充的想法,我得出:

(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 次计算就可以得出结果。

于 2010-07-17T19:02:24.510 回答
3

不是你问的许多问题的真正答案,但我对你没有问的问题有几条评论。

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

  2. for并且map都很懒惰,而且aset只有副作用。这种组合是灾难的根源:对于类似副作用的for行为,使用doseq.

  3. fordoseq允许多个绑定,因此您不需要像(显然)在 Python 中那样堆积它们的负载。

    (doseq [i (range (count cs))
            j (range (count is))] 
      ...)
    

    会做你想做的。

  4. #(first %)更简洁地写为first; 同样#(+ %1 %2)+

  5. 调用vec一堆不需要是向量的中间结果会减慢你的速度。具体来说,仅在您生成最终返回值时vec-add调用就足够了:如果它从不将它们用于随机访问,则没有理由将其中间结果转换为向量。vec(vec (reduce foo bar))foo

于 2011-03-26T03:42:28.213 回答