39
我正在研究Clojure中的一個應用程序，該應用程序需要乘以大型矩陣，並且遇到了與相同的Numpy版本相比的一些大型性能問題。 Numpy似乎能夠在一秒之內通過轉置乘以一個1,000,000x23的矩陣，而相同的clojure代碼則需要六分鐘。 （我可以打印Numpy生成的矩陣，所以它絕對評估一切）。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
啊，這些鏈接真的很有幫助。我會試着看看我能從一些矩陣包中得到什麼。 – 2012-01-17 19:18:48
不是C;大部分BLAS和LAPACK都是Fortran。 Python版本也不是「編譯」成C中的循環;它只是使用C/Fortran編寫的庫。最後，使用的算法是完全不同的。對Numpy的方法的描述完全與我們的原始海報顯示的代碼相關的可能性不大。 – pavpanchekha 2012-01-17 21:28:26
@pavpanchekha：另外，BLAS/LAPACK矩陣乘法的基本算法當然不是在上面的閉包中實現的天真的算法（其具有O（n 3）複雜性（對於方陣）。參見http://en.wikipedia.org /維基/ Matrix_multiplication＃Algorithms_for_efficient_matrix_multiplication – 2012-01-28 18:04:51