Matrices and vectors are the fundamental building blocks of many of the areas of programming I find most interesting: machine learning, simulations, statistical analysis, 3D games, computer-generated art etc.

Unfortunately, using matrices presents some problems in the Clojure world. There are actually a number of really good libraries available, but they are all incompatible and designed / optimised for different use cases. Choosing the right one is hard, and a lot of effort gets wasted converting between them and duplicating code.

This post is about how we might fix this in the Clojure world

The problem

At a superficial level, matrices and vectors are pretty simple. They are just grids of numbers with one or more dimensions. And they have a few well-defined mathematical operations you can apply to them, like matrix multiplication, cross products, evaluation of determinants etc.

You might think therefore that it would be easy to write a library to do all this. Mathematical objects are, after all, generally quite amenable to representation and manipulation in code. But it actually turns out that it is much more complicated, in part because there are so many trade-offs to consider:

Specialised vs. general matrices? – some matrices are special and can be stored or processed much more efficiently. A good example is the diagonal matrix, which is a square matrix that has zeros everywhere except the leading (top left to bottom right) diagonal. A diagonal matrix can be stored and processed much more efficiently than a general purpose matrix – on the other hand distinguishing it as a special case requires extra code and complexity. And there are many more special cases……

Specialised vectors? For similar reasons to matrices, you might want to make specialised vectors. A good example is the classic 3D vector – by implementing this as a simple class with x,y,z primitive values, your code would be much faster for 3D processing.

For similar reasons to matrices, you might want to make specialised vectors. A good example is the classic 3D vector – by implementing this as a simple class with x,y,z primitive values, your code would be much faster for 3D processing. Higher dimensions? Mostly you are interested in matrices (2 dimensional arrays) and vectors (1 dimensional arrays). But in some cases you want higher dimensions. The famous NumPy library in Python allows arbitrary numbers of dimensions for its matrix objects

Mostly you are interested in matrices (2 dimensional arrays) and vectors (1 dimensional arrays). But in some cases you want higher dimensions. The famous NumPy library in Python allows arbitrary numbers of dimensions for its matrix objects API compatibility? Different libraries / APIs expect their vectors and matrices in different formats. What do you need to be compatible with? Do you need to write a conversion layer?

Different libraries / APIs expect their vectors and matrices in different formats. What do you need to be compatible with? Do you need to write a conversion layer? Precision? Do you want 32-bit or 64-bit floats? What about integer matrices? what about arbitrary precision numbers?

Do you want 32-bit or 64-bit floats? What about integer matrices? what about arbitrary precision numbers? Concurrency? If your matrices are very large, you might want to use concurrent algorithms so that the work can be broken up over many CPUs

If your matrices are very large, you might want to use concurrent algorithms so that the work can be broken up over many CPUs Distribution? If you are working with extremely large matrices, you might even want to split up the computations over multiple machines.

If you are working with extremely large matrices, you might even want to split up the computations over multiple machines. Native code? There are lots of very fast native code libraries for matrix maths (e.g. BLAS / ATLAS). Do you want to use these to get access to highly optimised native routines? Or do you need the portability of pure Java?

All of these are genuine trade-off, i.e. there is no “right” answer – it depends on your circumstances and problems. As a result, there are a host of different matrix libraries available for Java / Clojure that occupy different areas of the tradeoff space. A few of these:

JBLAS / Clatrix – wraps the native BLAS code, so not pure Java. Extremely good performance for big matrices because of years of optimisations in the native libraries.

– wraps the native BLAS code, so not pure Java. Extremely good performance for big matrices because of years of optimisations in the native libraries. javax.vecmath – pretty basic but needed for compatibility with Java3D

– pretty basic but needed for compatibility with Java3D Vectorz – optimised for small matrices (2D/3D) and aimed at games and simulations. Very high throughput (billions of operations/sec) but not thread safe. Has a nice Clojure wrapper and is pure Java.

– optimised for small matrices (2D/3D) and aimed at games and simulations. Very high throughput (billions of operations/sec) but not thread safe. Has a nice Clojure wrapper and is pure Java. Parallel Colt – designed for multi-threaded concurrency, large matrices. Scientific focus (e.g. includes Complex matrices)

designed for multi-threaded concurrency, large matrices. Scientific focus (e.g. includes Complex matrices) EJML – fast all-round Java matrix library

– fast all-round Java matrix library …… the list goes on……

Thus, anyone trying to use a few simple matrices in Clojure faces a bewildering array of options.

What would be better?

There is obvious appeal in having a de-facto standard in a core matrix library that enables other code to be build on top of it without worrying about all the trade-offs.

In fact, this is the case in the Python world, where the NumPy library forms the basis for an impressive range of libraries and tools. As a result, Python is a pretty attractive choice for research and scientific computing.

So how can we get this kind of goodness in the Clojure world?

In my view, the ideal situation would be to create a common high level abstraction / API for matrix maths in Clojure with the following objectives:

Provide a shared, flexible high level matrix API in Clojure comparable to NumPy in Python.

Support multiple underlying matrix implementations . The existing matrix libraries mentioned above are all great for their specific use cases – so it would be great to allow users to select the one best suited to their needs but still access them through a common API. Ideally switching between implementations should be realtively painless / transparent.

. The existing matrix libraries mentioned above are all great for their specific use cases – so it would be great to allow users to select the one best suited to their needs but still access them through a common API. Ideally switching between implementations should be realtively painless / transparent. Offer good performance – an important requirement for applications that need matrix computations as these situations are often CPU-bound

an important requirement for applications that need matrix computations as these situations are often CPU-bound Provide a platform for building a library of common algorithms that will be useful across a wide range of domains (e.g. linear equation solvers)

If we could achieve all of these objectives, I think Clojure would have a much more compelling story for matrix computation, which would in turn make it much more attractive as a language of choice in many more domains (big data, machine learning, scientific computing, simulations, etc.)

The solution – core.matrix?

From a number of discussions online, it seems like there is general appetite and demand in the broader Clojure community for something like this. So as an experiment, I’ve decided to create a proof of concept matrix API to achieve these objectives.

In many languages building this kind of API / abstraction would be difficult. In Java, for example, to build this kind of API you would typically need to define a common interface and then develop a whole layer of wrapper classes to wrap the underlying implementations in a way that conforms to the common interface.

In Clojure however we have protocols – which are a powerful tool for building APIs and shared abstractions in a way that solves the expression problem. This makes it possible to create an API that can be extended to deal with different implementations *without* having to either modify the original implementation or create a wrapper. For example suppose we create a protocol for matrix addition:

(defprotocol PMatrixAdd "Protocol to support matrix addition on an matrices of same size" (matrix-add [m a]))

We can then extend this protocol to handle any different underlying implementation:

(extend-protocol PMatrixAdd my.special.MatrixClass (matrix-add [m a] (.specialMatrixAdd m (.someMethod a))))

Note that we are free to define how matrix addition is achieved for a specific implementation in any way we like: perhaps even deferring the computation or delegating it to a remote machine.

Finally, we can build a regular, easy to use Clojure API on top of the protocols. Perhaps you might want to define vector addition, with something like the following:

(defn + "Matrix addition operator" ([a] a) ([a b] (if (and (number? a) (number? b)) (clojure.core/+ a b) (matrix-add a b))) ([a b & more] (reduce + (+ a b) more)))

With this kind of approach, we can create a common API that can handle all different kinds of underlying matrix implementations. Furthermore, protocols are very efficient because they make use of the highly optimised method dispatch on the JVM.

Overall, I think this is a very promising approach to creating a common matrix API / abstraction in Clojure. If anyone is interested, here is the experimental repo on GitHub where I am developing this idea:

Thoughts / suggestions / patches welcome as always!