Clojure 与 Numpy 中的矩阵乘法
Posted
技术标签:
【中文标题】Clojure 与 Numpy 中的矩阵乘法【英文标题】：Matrix Multiplication in Clojure vs Numpy 【发布时间】：2012-02-12 13:27:30 【问题描述】：我正在使用 Clojure 开发一个应用程序，该应用程序需要将大型矩阵相乘，与相同的 Numpy 版本相比，我遇到了一些较大的性能问题。 Numpy 似乎能够在一秒钟内将一个 1,000,000x23 矩阵乘以它的转置，而等效的 clojure 代码则需要六分钟以上。 （我可以从 Numpy 打印出结果矩阵，所以它肯定会评估所有内容）。
我在这个 Clojure 代码中做错了什么吗？是否有一些我可以尝试模仿的 Numpy 技巧？
这是蟒蛇：
import numpy as np
def test_my_mult(n):
A = np.random.rand(n*23).reshape(n,23)
At = A.T
t0 = time.time()
res = np.dot(A.T, A)
print time.time() - t0
print np.shape(res)
return res
# Example (returns a 23x23 matrix):
# >>> results = test_my_mult(1000000)
#
# 0.906938076019
# (23, 23)
还有clojure：
(defn feature-vec [n]
(map (partial cons 1)
(for [x (range n)]
(take 22 (repeatedly rand)))))
(defn dot-product [x y]
(reduce + (map * x y)))
(defn transpose
"returns the transposition of a `coll` of vectors"
[coll]
(apply map vector coll))
(defn matrix-mult
[mat1 mat2]
(let [row-mult (fn [mat row]
(map (partial dot-product row)
(transpose mat)))]
(map (partial row-mult mat2)
mat1)))
(defn test-my-mult
[n afn]
(let [xs (feature-vec n)
xst (transpose xs)]
(time (dorun (afn xst xs)))))
;; Example (yields a 23x23 matrix):
;; (test-my-mult 1000 i/mmult) => "Elapsed time: 32.626 msecs"
;; (test-my-mult 10000 i/mmult) => "Elapsed time: 628.841 msecs"
;; (test-my-mult 1000 matrix-mult) => "Elapsed time: 14.748 msecs"
;; (test-my-mult 10000 matrix-mult) => "Elapsed time: 434.128 msecs"
;; (test-my-mult 1000000 matrix-mult) => "Elapsed time: 375751.999 msecs"
;; Test from wikipedia
;; (def A [[14 9 3] [2 11 15] [0 12 17] [5 2 3]])
;; (def B [[12 25] [9 10] [8 5]])
;; user> (matrix-mult A B)
;; ((273 455) (243 235) (244 205) (102 160))
更新：我使用JBLAS 库实现了相同的基准测试，并发现了巨大的速度改进。感谢大家的意见！是时候用 Clojure 包装这个傻瓜了。这是新代码：
(import '[org.jblas FloatMatrix])
(defn feature-vec [n]
(FloatMatrix.
(into-array (for [x (range n)]
(float-array (cons 1 (take 22 (repeatedly rand))))))))
(defn test-mult [n]
(let [xs (feature-vec n)
xst (.transpose xs)]
(time (let [result (.mmul xst xs)]
[(.rows result)
(.columns result)]))))
;; user> (test-mult 10000)
;; "Elapsed time: 6.99 msecs"
;; [23 23]
;; user> (test-mult 100000)
;; "Elapsed time: 43.88 msecs"
;; [23 23]
;; user> (test-mult 1000000)
;; "Elapsed time: 383.439 msecs"
;; [23 23]
(defn matrix-stream [rows cols]
(repeatedly #(FloatMatrix/randn rows cols)))
(defn square-benchmark
"Times the multiplication of a square matrix."
[n]
(let [[a b c] (matrix-stream n n)]
(time (.mmuli a b c))
nil))
;; forma.matrix.jblas> (square-benchmark 10)
;; "Elapsed time: 0.113 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 100)
;; "Elapsed time: 0.548 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 1000)
;; "Elapsed time: 107.555 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 2000)
;; "Elapsed time: 793.022 msecs"
;; nil
【问题讨论】：
【参考方案1】：Python 版本正在编译为 C 中的循环，而 Clojure 版本正在为此代码中的每个映射调用构建一个新的中间序列。您看到的性能差异很可能来自数据结构的差异。
为了比这更好，您可以使用像 Incanter 这样的库，或者按照 this SO question 中的说明编写自己的版本。另见this one、neanderthal 或nd4j。如果您真的想保留序列以保持惰性评估属性等，那么您可以通过查看transients 进行内部矩阵计算来获得真正的提升
编辑：忘记添加调优clojure的第一步，打开“反射警告”
【讨论】：
啊，这些链接真的很有帮助。我会试着看看我能从一些矩阵包中得到什么。 不是 C；大部分 BLAS 和 LAPACK 都是 Fortran。 Python 版本也没有“编译到”C 中的循环；它只是使用用 C/Fortran 编写的库。最后，使用的算法完全不同。 Numpy 方法的描述不太可能与我们原始海报显示的代码相关。 @pavpanchekha：另外，BLAS/LAPACK 矩阵乘法的底层算法当然不是上面闭包中实现的简单算法（它具有 O(n³) 复杂度（对于方阵）。参见@987654327 @ 作为一个很好的答案，我冒昧地添加了几个额外的现代选项。【参考方案2】：Numpy 链接到 BLAS/Lapack 例程，这些例程已在机器架构级别优化了数十年，而 Clojure 则以最直接和幼稚的方式实现乘法。
任何时候您需要执行重要的矩阵/向量运算，您可能应该链接到 BLAS/LAPACK。
唯一不会更快的情况是对于来自语言的小矩阵，其中在语言运行时和 LAPACK 之间转换数据表示的开销超过了计算所花费的时间。
【讨论】：
【参考方案3】：我刚刚在Incanter 1.3 和jBLAS 1.2.1 之间进行了一场小型枪战。代码如下：
(ns ml-class.experiments.mmult
[:use [incanter core]]
[:import [org.jblas DoubleMatrix]])
(defn -main [m]
(let [n 23 m (Integer/parseInt m)
ai (matrix (vec (double-array (* m n) (repeatedly rand))) n)
ab (DoubleMatrix/rand m n)
ti (copy (trans ai))
tb (.transpose ab)]
(dotimes [i 20]
(print "Incanter: ") (time (mmult ti ai))
(print " jBLAS: ") (time (.mmul tb ab)))))
在我的测试中，Incanter 在普通矩阵乘法中始终比 jBLAS 慢 45%。但是，Incanter trans 函数不会创建矩阵的新副本，因此 jBLAS 中的 (.mmul (.transpose ab) ab) 占用两倍的内存，并且仅比 Incanter 中的 (mmult (trans ai) ai) 快 15%。
鉴于 Incanter 丰富的功能集（尤其是绘图库），我认为我不会很快切换到 jBLAS。尽管如此，我还是希望看到 jBLAS 和 Parallel Colt 之间的另一场枪战，也许值得考虑在 Incanter 中用 jBLAS 替换 Parallel Colt？ :-)
编辑：这是我在（相当慢的）PC 上获得的绝对数字（以毫秒为单位）：
Incanter: 665.362452
jBLAS: 459.311598
numpy: 353.777885
对于每个库，我从 20 次运行中选出了最佳时间，矩阵大小为 23x400000。
PS。 Haskell hmatrix 结果接近 numpy，但我不确定如何正确对其进行基准测试。
【讨论】：
【参考方案4】：Numpy 代码使用内置库，这些库是在过去几十年里用 Fortran 编写的，并由作者、您的 CPU 供应商和您的操作系统分销商（以及 Numpy 人员）进行了优化，以实现最佳性能。您刚刚完成了矩阵乘法的完全直接、明显的方法。确实，性能不同并不奇怪。
但是，如果您坚持在 Clojure 中执行此操作，请考虑查找 better algorithms，使用直接循环而不是像 reduce 这样的高阶函数，或者为 Java 找到合适的矩阵代数库（我对此表示怀疑）是 Clojure 中的好东西，但我真的不知道）由一位称职的数学家编写。
最后，看看如何正确编写快速的 Clojure。使用类型提示，在你的代码上运行一个分析器（惊喜！你的点积函数用尽了大部分时间），然后把高级特性放到你的紧密循环中。
【讨论】：
一个好的 Java 矩阵代数库就是一个好的矩阵代数库，供 Clojure 等待包装器使用。 你说的。我已经找到了 JBLAS，并且正在着手整理它。【参考方案5】：正如@littleidea 和其他人所指出的，您的 numpy 版本正在使用 LAPACK/BLAS/ATLAS，这将比您在 clojure 中所做的任何事情都快得多，因为它已经经过多年的微调。 :)
也就是说，Clojure 代码的最大问题是它使用了 Doubles，就像在盒装的 doubles 中一样。我称之为“懒惰的 Double”问题，我在工作中遇到过很多次。到目前为止，即使是 1.3，clojure 的集合也不是原始友好的。 （您可以创建一个原始向量，但它对您没有任何帮助，因为所有 seq. 函数最终都会对它们进行装箱！我还应该说 1.3 中的原始改进非常好并且最终会有所帮助..我们只是集合中不是 100% 支持 WRT 原语。）
在 clojure 中进行任何类型的矩阵数学运算时，您确实需要使用 java 数组，或者更好的是矩阵库。 Incanter 确实使用了 parrelcolt，但您需要注意您使用的 incanter 函数......因为其中很多使矩阵可序列化，最终将双打装箱，从而为您提供与您当前看到的类似的性能。 （顺便说一句，我有自己设置的 parrelcolt 包装器，如果您认为它们会有所帮助，我可以发布它们。）
为了使用 BLAS 库，您在 java-land 中有几个选项。使用所有这些选项，您必须支付 JNA 税……您的所有数据都必须在被处理之前被复制。当您执行诸如矩阵分解之类的 CPU 绑定操作并且其处理时间比复制数据所需的时间长时，这种税是有意义的。对于小矩阵的简单操作，那么留在 java-land 中可能会更快。您只需要像上面所做的那样进行一些测试，看看哪种方法最适合您。
以下是从 java 中使用 BLAS 的选项：
http://jblas.org/
http://code.google.com/p/netlib-java/
以上API：http://code.google.com/p/matrix-toolkits-java/我应该指出 parrelcolt 使用 netlib-java 项目。这意味着，我相信，如果您正确设置它，它将使用 BLAS。但是，我还没有验证这一点。有关 jblas 和 netlib-java 之间差异的解释，请参阅我在 jblas 的邮件列表上开始的这个线程：
http://groups.google.com/group/jblas-users/browse_thread/thread/c9b3867572331aa5
我还应该指出通用 Java 矩阵包库：
http://sourceforge.net/projects/ujmp/
它包含了我提到的所有库，还有一些！尽管我没有过多地研究 API，但我知道它们的抽象有多么泄漏。这似乎是一个不错的项目。我最终使用了我自己的 parrelcolt clojure 包装器，因为它们足够快，而且我实际上非常喜欢 colt API。 （Colt 使用函数对象，这意味着我可以轻松地传递 clojure 函数！）
【讨论】：
【参考方案6】：如果您想在 Clojure 中处理数字，我强烈建议您使用 Incanter 而不是尝试滚动您自己的矩阵函数等。
Incanter 在后台使用Parallel Colt，速度非常快。
编辑：
截至 2013 年初，如果您想在 Clojure 中处理数字，我强烈建议您查看 core.matrix
【讨论】：
【参考方案7】：Numpy 针对线性代数进行了高度优化。当然对于大型矩阵，大部分处理都在本机 C 代码中。
为了匹配这种性能（假设它在 Java 中是可能的），您必须去除 Clojure 的大部分抽象：在迭代大型矩阵时不要将 map 与匿名函数一起使用，添加类型提示以启用原始 Java数组等
可能最好的选择就是使用针对数值计算优化的现成 Java 库（http://math.nist.gov/javanumerics/ 或类似的）。
【讨论】：
实际上，numpy 甚至没有必要在 C 中进行计算。它只是调用BLAS 例程。如果您安装了MKL 或ATLAS，（并且如果numpy 配置为使用它们）它将调用他们的BLAS 例程进行矩阵乘法，其中一部分基本上是手动调整的汇编。无论如何，我的观点是 numpy （如果配置正确）将击败相当于用 C 编写的 OP 的 clojure 示例，所以最好的选择可能更像是 Incanter 等（正如你所建议的，我只是漫无目的） .
Incanter 使用 parallel-colt，它也可以利用 MKL / ATLAS【参考方案8】：
我没有任何具体的答案给你；只是一些建议。
-
使用分析器找出时间花在哪里
设置 warn-on-reflection 并在需要时使用类型提示
您可能不得不放弃一些高级构造并使用循环递归来挤出最后一盎司的性能
IME，Clojure 代码的性能应该非常接近 Java（2 或 3X）。但你必须努力。
【讨论】：
【参考方案9】：只有在有意义的情况下才使用 map()。这意味着：如果你有一个特定的问题，比如将两个矩阵相乘，不要尝试 map() 它，只需将矩阵相乘即可。
我倾向于只在它具有语言意义时才使用 map()（即，如果程序真的比没有它更具可读性）。矩阵相乘是如此明显的循环，以至于映射它毫无意义。
你的。
佩德罗·福图尼。
【讨论】：
以上是关于Clojure 与 Numpy 中的矩阵乘法的主要内容，如果未能解决你的问题，请参考以下文章