Dot and cross products in Clojure/JVM

2012-05-31

Dot and cross vector products are essential for geometric calculations in 3D. Incanter which is the best known "scientific" library for Clojure wraps 2D matrices from the Parallel Colt, and can do some linear algebra too, but it doesn't wrap neither 1D vectors, nor vector operations.

I tried several options for vector math in Clojure:

  • use Incanter matrix to store both matrices and vectors, and write vector operations by hand (it did work, and was not that bad)

  • use Incanter matrix for matrices and Parallel Colt's 1D matrices for vectors (with inevitable conversion costs between 1D and 2D representations)

  • use Apache Commons Math library (it has built-in vector products in its 3D geometry package)

First option. Using just Incanter.Core, [incanter/incanter-core "1.3.0"] in lein :dependencies. It adds 5.3 MB to the lib folder.

The naive implementation for dot product is just (apply + (map * xs ys)), and cross product uses matrix multiplication. We can use Incanter.Core's elementwise plus, minus, and mult operations to do the rest of the vector math.

Implementation: naive.clj

Benchmark:

user> (use 'incanter.core '(incanter.contrib.vecmath (naive :as naive)))
nil
user> (def x (matrix [1 2 3]))
#'user/x
user> (def y (matrix [4 5 6]))
#'user/y
user> (bench (naive/dot x y))
...
             Execution time mean : 14.427771 us ...
...
user> (bench (naive/cross x y))
...
             Execution time mean : 41.505528 us ...
...

Things worked surprisingly well for small (3D) vectors, but dot deteriorates very quickly if you multiply longer vectors.

Second option. Using Colt's DenseDoubleMatrix1D for vectors, BLAS' ddot for dot product, and Incanter.Core for everything else. To distinguish Colt vector from the Clojure vectors, I called it mvector (and its factory function is called the same). To convert it back to Incanter's matrix, I use as-matrix. dot and cross convert their arguments if necessary. The former works best when it receives mvectors, the latter still needs to deconstruct vector elements, and works best with just 3-sequences as the first vector.

Implementation: mvector.clj

Benchmark:

user> (require '(incanter.contrib.vecmath (mvector :as mvector)))
nil
user> (def xv (mvector/mvector [1 2 3]))
#'user/xv
user> (def yv (mvector/mvector [4 5 6]))
#'user/yv
user> (bench (mvector/dot xv yv))
...
             Execution time mean : 15.618065 ...
...
user> (bench (mvector/cross xv yv))
...
             Execution time mean : 62.234059 us ...
...

Oops. Not much improvement for dot, and significant conversion costs in cross. The only benefit is that this dot is fast even on long vectors (who cares).

The last option I tried was Commons Math. Just [org.apache.commons/commons-math3 "3.0"] in lein :dependencies. The entire commons-math3-3.0.jar is only 1.3 MB. No Clojure wrappers, just plain Java API. There is already a cross product implementation in 3D Euclidean geometry package.

user=> (import '(org.apache.commons.math3.geometry.euclidean.threed Vector3D))
org.apache.commons.math3.geometry.euclidean.threed.Vector3D
user=> (def v3d (Vector3D. (double-array [1 2 3])))
#'user/v3d
user=> (def u3d (Vector3D. (double-array [4 5 6])))
#'user/u3d
user=> (bench (.dotProduct v3d u3d))
...
             Execution time mean : 15.088191 us  ...
...
user=> (bench (.crossProduct v3d u3d))
...
             Execution time mean : 17.698480 us  ...
...

Of all three options I now prefer Commons Math. It is small and reasonably fast. Its Java API is easy to use directly. And it offers many other useful algorithms beyond just simple vector math.

For reference, in Python 2.7.3 & SciPy 0.10.1 on the same machine I have a much faster dot, but slower cross product:

In [2]: x=array([1,2,3],dtype=float64)

In [3]: %timeit dot(x,x)
100000 loops, best of 3: 3.4 us per loop
...
In [4]: y = array([4,5,6],dtype=float64)

In [5]: %timeit cross(x,y)
10000 loops, best of 3: 45.4 us per loop