Dot and cross products in Clojure/JVM
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
: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
operations to do the rest of the vector math.
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
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
(and its factory function is called the same). To convert it back to Incanter's
matrix, I use
cross convert their arguments if
necessary. The former works best when it receives
latter still needs to deconstruct vector elements, and works best with
just 3-sequences as the first vector.
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
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
In : x=array([1,2,3],dtype=float64) In : %timeit dot(x,x) 100000 loops, best of 3: 3.4 us per loop ... In : y = array([4,5,6],dtype=float64) In : %timeit cross(x,y) 10000 loops, best of 3: 45.4 us per loop