title: Matrix Multiplication in Clojure vs Numpy
body: <p>I'm working on an application in Clojure that needs to multiply large matrices and am running into some large performance issues compared to an identical Numpy version. Numpy seems to be able to multiply a 1,000,000x23 matrix by its transpose in under a second, while the equivalent clojure code takes over six minutes. (I can print out the resulting matrix from Numpy, so it's definitely evaluating everything).</p>
  
  <p>Am I doing something terribly wrong in this Clojure code? Is there some trick of Numpy that I can try to mimic?</p>
  
  <p>Here's the python:</p>
  
  <pre><code>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):
  # &gt;&gt;&gt; results = test_my_mult(1000000)
  # 
  # 0.906938076019
  # (23, 23)
  </code></pre>
  
  <p>And the clojure:</p>
  
  <pre><code>(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) =&gt; "Elapsed time: 32.626 msecs"
  ;; (test-my-mult 10000 i/mmult) =&gt; "Elapsed time: 628.841 msecs"
  
  ;; (test-my-mult 1000 matrix-mult) =&gt; "Elapsed time: 14.748 msecs"
  ;; (test-my-mult 10000 matrix-mult) =&gt; "Elapsed time: 434.128 msecs"
  ;; (test-my-mult 1000000 matrix-mult) =&gt; "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&gt; (matrix-mult A B)
  ;; ((273 455) (243 235) (244 205) (102 160))
  </code></pre>
  
  <p>UPDATE: I implemented the same benchmark using the <a href="http://jblas.org/" rel="noreferrer">JBLAS</a> library and found massive, massive speed improvements. Thanks to everyone for their input! Time to wrap this sucker in Clojure. Here's the new code:</p>
  
  <pre><code>(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&gt; (test-mult 10000)
  ;; "Elapsed time: 6.99 msecs"
  ;; [23 23]
  
  ;; user&gt; (test-mult 100000)
  ;; "Elapsed time: 43.88 msecs"
  ;; [23 23]
  
  ;; user&gt; (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&gt; (square-benchmark 10)
  ;; "Elapsed time: 0.113 msecs"
  ;; nil
  ;; forma.matrix.jblas&gt; (square-benchmark 100)
  ;; "Elapsed time: 0.548 msecs"
  ;; nil
  ;; forma.matrix.jblas&gt; (square-benchmark 1000)
  ;; "Elapsed time: 107.555 msecs"
  ;; nil
  ;; forma.matrix.jblas&gt; (square-benchmark 2000)
  ;; "Elapsed time: 793.022 msecs"
  ;; nil
  </code></pre>