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.**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**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?**Precision?**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**Distribution?**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.**javax.vecmath**– 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.**Parallel Colt -**designed for multi-threaded concurrency, large matrices. Scientific focus (e.g. includes Complex matrices)**EJML**– 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. - Offer
**good performance -**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!

Looks interesting. I’m interested in having matrices in clojure as well. After my exams I’ll probably take a look at your project. In the mean time, have fun!

I do think that using a protocol to create a matrix type is a good idea. ie. “PIndexedAccess” is a good idea. However I am not appealed by the fact that you want to make protocols for every single operation. eg: PMatrixMultiply

The reason you do this (as far as my understanding goes) is so that we can make special cases of matrix multiply for special types of matrices.

Two problems with this approach:

1. If you want to multiply two different “type” of matrices. Then which of the two multiply methods are you going to call? You will have to make priorities, how?

2. for every new kind of operation you will have to extend existing types.

I personally appeal more to the following approach:

you define functions like

‘(multiply m1 m2)’

‘(solve mA mx)’

‘(solve-sparse mA mx)’

Summarized

I like PIndexedAccess. This can be used to very efficiently store certain matrices. eg: Diagonal matrices

I think it is the programmers/scientist job to now what type of structures he is working with. So I would create separate functions for equivalent operations/algorithms but optimized for specific types of matrices.

Good thoughts.

On 1., my current idea is that the first matrix takes charge initially, since this determines what matrix implementation we are using (e.g. JBLAS). The implementation can then decide to either a) operate on the second matrix directly if it knows how or b) call a coercion methods to convert the second matrix to a type that it understands.

On 2, I actually think this is pretty much inevitable. There will always be a need to extend each protocol to different types of matrix implementation. Again, it will be up to the implementation to decide how to perform the operation, e.g. if the implementation uses sparse matrices it can call the optimised sparse version of the algorithm

I still don’t know how it will all work out of course – it’s pretty experimental at the moment so anything could change. In my experience, if you are trying to define a common API you need to do 2-3 different implementations before you work out the API quirks.

I feel that if you really want maximum flexibility and while still having a knob for performance, then you must build mathematical expressions and then compile it down to a target. The reach is greater, since you can extend to any implementation, including outside the language since you are building emitters. You can manipulate and analyze the expressions at compile time, allowing simplifications, proofs, etc. Since you usually you compile once, but run many times, you can use slower but more flexible dispatching methods such as multimethods or predicate dispatch on the expressions. I am currently working on my own implementation but it is highly experimental but I hope someone can steal my ideas and make a better implementation. https://github.com/bmillare/dj.math/blob/master/project.clj

(Note if you don’t want to use dj, you’ll need to use git to download its dependencies).

Totally agree! I think doing manipulation at the symbolic / logic level is critical.

I’m hoping that the core.matrix stuff should work as a reasonable target for such compilation, i.e. it should be at the right level for a symbolic computation to produce an expression, where the individual operations can then be handled by an arbitrary implementation, e.g. JBLAS.

Perhaps dj.math could try using the core.matrix API as a target? I added an issue in your project to consider this.

This is very interesting. I’ve not had time off late to stay abreast with Clojure, but I remember wanting something like this. I had used Parallel Colt to do some simple benchmarking, and found its performance to be good.

There is Incanter for Clojure, although that is a different use scenario, and doesn’t have (when I looked last) sparse matrices.

I’d be glad to get back into Clojure and help out with making a Numpy for the Clojure world.

Would be great to have your help! I agree with you that having a NumPy for the Clojure world would be very powerful.

I’m actually thinking that Incanter would become a consumer of the library – I raised this on the Incanter mailing list and there seemed to be a lot of interest.

Array computation is important. I’m glad that players within the Clojure community are engaging it. I want to stress that there is a large body of accrued wisdom on this subject.

An anecdote from recent history: Several years ago Python had an issue with two competing array implementations, Numeric and NumArray. This caused significant havoc before NumPy arrived. It is unlikely that Python would be a tenth as effective without NumPy’s unifying effect. This is probably obvious but I’ll say it anyway, care should go in to the design of any API expected to be used by the entire community; APIs need to anticipate future use cases.

I suggest looking deeply at the principles behind NumPy; they extend beyond broadcasted operators and indexing. I also suggest looking beyond NumPy; the scientific Python community is largely moving on to better things; NumPy is seen as somewhat antiquated.

For example the growth of parallel and heterogeneous computing caused the Python community to value the separation of the matrix expression from the numeric computation; matrix-expressions-as-data is a big deal now. The rising use of computation in statistical communities has highlighted the importance of heterogeneous datatypes (like pandas or R’s dataframe.)

I’m excited to see what comes of this project.

Great insights, thanks! I’m certainly looking at the Python world for inspiration and learnings.

To some extent I think that “matrix expressions as data” will come quite naturally in the Clojure world. Once we have a clean API for the matrix computations themselves, I think it should be possible to build higher level functionality on top (e.g. using core.logic to optimise matrix expressions etc.)

Perhaps traditional Common Lisp mechanisms, such as keywords and default values might come in handy.

In any case, you’ve certainly identified a need, and it’s good of you to begin investigating “core.matrix”.

Good idea – there’s definitely a lot of syntactical niceness that can come from the CL techniques.

If you have any good ideas in particular, do drop them into the project! I’m keen to get as many eyeballs on the design as possible at this early stage.

Hi, I’m the main developer of jblas and I’m all for making jblas better at multilanguage support. Originally, I wrote jblas as a layer to use in some private jruby matlab like clone. I’ve been thinking a lot about how to make jblas better, for example by extracting the actual computational routines in compact classes which are better suited to build the language specific syntactic sugar on top. I’m personally working a lot with Scala these days, so having support for Scala, Java, and Clojure looks like a nice set of languages.

If you have some ideas or want to discuss with me, head over to the jblas-users group: http://groups.google.com/group/jblas-users

-Mikio

Thanks Mikio! That’s great to hear, good JBLAS support is definitely one of the things we are very keen to interface with. Happy to give you some input from the Clojure perspective.

Thanks for the info and the libraries. I’m converting a program from Perl to clojure (for direct Neo4j access, long story…) and matrix math is exactly what I was missing. I’ve been using PDL in Perl, which is very fast, and was needed an analog in clojure.

Cool, let me know how you get on with it! FWIW, vectorz-clj is probably the most mature implementation so far:

https://github.com/mikera/vectorz-clj