More fun with NumPy, CuPy, Clojure and GPU acceleration. Hold my Cider 2!

The following two tables contain the summary of the micro-benchmarks that I run on my CPU and GPU.

Since NumPy and CuPy forcefully do float64 computation regardless of data type, there is one column for them. Neanderthal () follows what it's told, so I include both float64 (double in Javaspeak) for direct comparison, and float32 (float) as a reference.

NumPy and CuPy code is run from the interactive Python console, while Clojure uses interactive REPL.

Note that my code is not a port of NumPy, nor a drop-in replacement; this is simply a head-to-head comparison of two different implementations. Python has its strengths and weaknesses, but JVM has its own constraints, too, so don't blame me for either of them, nor expect that I care about fixing them.

To prevent folks bugging me about why I called NumPy from Clojure, I wrote a Python script this time. Micro benchmarks in Python code that I see on the Internet are clunky at best, and more sophisticated tools are equally ridiculous. Having to call my functions in a string- embedded DSL? Thanks, but no.

So I did the lazy thing, and used a simple time_ns() based measurement of a loop that calls the corrcoef function enough times to mitigate the imprecision of the timestamp that Python gets. I believe that's a good approximation of how would a user experience this function when calling it from the code. It's good enough considering that we are calling exactly one function that delegates to native code anyway.

On the CPU, it looks like this:

def numpy_corrcoef(m, n, reps):
    a = numpy.random.random(m * n).reshape(m,n)
    s = time.time_ns()
    for x in range(reps):
        numpy.corrcoef(a)
    e = time.time_ns()
    print((e-s)/reps/(10**6), " milliseconds")

numpy_corrcoef(1000, 100000, 1)
numpy_corrcoef(10, 10000, 1000)
# etc.

On the GPU, it takes care of synchronization:

def cupy_corrcoef(m, n, reps):
    a = cupy.random.random(m * n).reshape(m,n)
    cupy.cuda.Stream.null.synchronize()
    s = time.time_ns()
    for x in range(reps):
        cupy.corrcoef(a)
    cupy.cuda.Stream.null.synchronize()
    e = time.time_ns()
    print((e-s)/reps/(10**6), " milliseconds")

cupy_corrcoef(1000, 100000, 1)
# etc.

In Clojure, we have sophisticate benchmarking capabilities right there in the REPL through the Criterium library. Although would like to use its quick-bench, as it gives reliable and precise measurements with one simple line of code, I opted to follow the Python way, and use time and the appropriate loop, to make it as similar as the Python code as possible.

(defn bench-neanderthal-cpu [m n factory reps]
  (with-release [x (rand-uniform! (ge factory m n))]
    (time
     (dotimes [i reps]
       (matrix-correlation! x)))))

(defn bench-neanderthal-gpu [m n factory reps]
  (with-default
    (cuda/with-engine factory default-stream
      (with-release [x (rand-uniform! (cuge m n))]
        (synchronize!)
        (time
         (do
           (dotimes [i reps]
             (release (matrix-correlation! x)))
           (synchronize!)))))))

In the previous article you can find the code for matrix-correlation!. That implementation did a full general matrix multiplication, although the correlation matrix is symmetric, since for most matrix sizes it is faster to do twice as many flops as necessary instead of using logic to avoid duplicate work. Especially on the GPU, the matrix should have to be unrealistically large for the symmetric matrix multiply to pay off.

Also, consider what corrcoef computes: correlation between different variables. There will typically be much fewer variables than observations, so the matrix will be a long rectangle. There, my old implementation will be optimal, even for large matrices.

However, I confirmed by measurement and comparison that NumPy uses symmetric multiply for large matrices. This cannot be seen directly from the code, but is 100% certain since it takes NumPy 1.5 second to do the full matrix multiply of this matrix, and around 1 sec for the whole corrcoef.

That's why I added this option to Neanderthal (), and made a simple heuristic to use symmetric multiply for huge fat matrices (the mmt! function, and GEMM for everything else. That made the code less elegant, and several lines longer, but, hey, everything for performance! :)

(defn matrix-correlation! [a!]
  (with-release [ones (entry! (vctr a! (ncols a!)) 1.0)
                 work (mv a! ones)]
    (let [a-mean (rk! (/ -1.0 (dim ones)) work ones a!)]
      (let-release [large (and (< 512 (mrows a!)) (< 8196 (ncols a!)))
                    sigma (sy a! (mrows a!))]
        (if large
          (mmt! a-mean sigma)
          (mm! 1.0 a-mean (trans a-mean) (view-ge sigma)))
        (sqrt! (copy! (dia sigma) work))
        (with-release [sigma-x*y (rk work)]
          (div! (view-ge sigma) (view-ge sigma-x*y)))
        sigma))))

I guess that this is enough info. If you feel inspired, you can try to write the Python equivalent, and see whether it is faster than the NumPy's and CuPy's implementation.

The books that I write are written for all programmers, not just Clojure programmers. I use Clojure because it is very enjoyable due to its elegance and power, but almost all the lessons in the book are applicable in other environments. Clojure helped me in writing the whole stack of these libraries that support CPU, GPU, CUDA, OpenCL, etc. in surprisingly small numbers of lines of code.

It can help you, too, to be more productive, but even if Clojure is not your cup of tea, the books can help you learn this stuff and apply it in other programming languages, since the books make clear connection from theory to implementation, without skipping any steps in between!

The key difference compared to other books, is that here you can try each small building block, which is usually one or a few lines of code directly in the interactive REPL, see the results right away, experiment, and understand the theoretic concept that is demonstrated. Equally important is that the performance of that code is top notch right away!

• Numerical Linear Algebra for Programmers
• Deep Learning for Programmers

There's already 400 pages written. By subscribing to these books you get the access to the drafts right away, and the final book this summer. There's no middle man; 100% of the proceeds goes towards my work on these open source libraries!