Dot and cross products in Clojure/JVM
20120531
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 builtin vector products in its 3D geometry package)
First option. Using just Incanter.Core, [incanter/incantercore "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 asmatrix
. dot
and cross
convert their arguments if
necessary. The former works best when it receives mvector
s, the
latter still needs to deconstruct vector elements, and works best with
just 3sequences 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/commonsmath3 "3.0"]
in lein
:dependencies
. The entire commonsmath33.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. (doublearray [1 2 3])))
#'user/v3d
user=> (def u3d (Vector3D. (doublearray [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