"core.matrix" proposal

[Mikera]

Hello all,

I've been experimenting with a common API / abstraction for matrix and vector maths in Clojure:

  https://github.com/mikera/matrix-api

The idea is:

 - Provide a clear, consistent API for matrix and vector operations

 - Support multiple different underlying implementations (e.g. native code via JBLAS vs pure Java like vectorz-clj)

 - Provide a base which other libraries that need to consume matrices can build upon (e.g. Incanter)

 - Offer good performance while still presenting a reasonably flexible high level API

I think this could be very useful for the Clojure community, especially given the interest in big data, simulations, machine learning, 3D graphics etc. If it goes well and there is enough interest, I guess that this could form the basis for a future "core.matrix" library.

Comments / ideas / patches welcome.

  Mike.

[Konrad Hinsen]

Mikera writes:

> I think this could be very useful for the Clojure community, especially given the
> interest in big data, simulations, machine learning, 3D graphics etc. If it goes well
> and there is enough interest, I guess that this could form the basis for a future
> "core.matrix" library.
>
> Comments / ideas / patches welcome.

I haven't looked at your library yet, but my immediate comment is that
I'd much prefer to have an API for N-dimensional arrays rather than
just for matrices (2d) and vectors (1d).

Konrad.

[Rob Lachlan]

I really like this idea -- I think there's a need for a dedicated matrix computation library in clojure.  I really like the idea of having matrix operations implemented in clojure (I think that you have this in persistent_vector.clj) but also being able to call on java libraries.

[James Sofra]

+1 on Konrad's comment about N-dimensional arrays for me.

I think the issue of using immutable arrays (which can be important for performance) in Clojure is interesting to think about.

I have had some ideas about how to use arrays with restricted or limited mutability but haven't solidified anything yet.

[Mikera]

Initial tests seem to suggest that supporting N-dimensional arrays should be pretty easy in terms of the API itself. We could also provide fast-paths for 1D and 2D arrays.

Of course, actual support for N-dimensional arrays would depend on the underlying implementation. 

[Mikera]

Yep, the idea is to be flexible enough to support many different implementations.

The pure Clojure version should be very easy to use and flexible since it uses regular Clojure persistent vectors. The trade-off is of less performance compared to the Java/native implementations.

As an added bonus, writing a pure Clojure version is useful for testing / validating the design of the API before we extend it to more complex implementations.

[Konrad Hinsen]

Mikera writes:

> Initial tests seem to suggest that supporting N-dimensional arrays
> should be pretty easy in terms of the API itself. We could also
> provide fast-paths for 1D and 2D arrays.
>
> Of course, actual support for N-dimensional arrays would depend on
> the underlying implementation.

There are a few array implementations that support N-d arrays, for example
the ones in JHDF5 (https://wiki-bsse.ethz.ch/display/JHDF5) and in netCDF
(http://www.unidata.ucar.edu/software/netcdf-java/).

> Yep, the idea is to be flexible enough to support many different implementations.

The real problem I see in making all this practically usable is
efficient conversion between implementations. The fragmentation of the
scientific computing ecosystem in Java makes it nearly impossible to
get any real work done without some conversion between array
libraries. If I want to read a matrix from an HDF5 file and then
calculate its eigenvalues, I need to convert from HDF5 arrays to Colt
(or whatever) arrays. Some libraries provide a way to construct an N-d
array from a 1-d data storage area plus a shape vector, without
copying the data. I think it's important to leverage such
functionality in a Clojure API in order to let conversion happen as
much as possible behind the scenes.

BTW, I started my own attempt at something similar a while ago, but
never found the time to get it to a usable state:

https://code.google.com/p/clj-multiarray/

Konrad.

[Brandon Bloom]

There's a bunch of useful operations that games & animation systems perform frequently that are less common in other uses of linear algebra. For example: linear and spline interpolations.

The DirectXMath library is worth studying: http://msdn.microsoft.com/en-us/library/windows/desktop/hh437833(v=vs.85).aspx

Similarly, the XNA math libraries: http://msdn.microsoft.com/en-us/library/bb196942(v=xnagamestudio.35).aspx

[Mikera]

On Monday, 7 January 2013 15:53:20 UTC+8, Konrad Hinsen wrote:

Mikera writes:

 > Initial tests seem to suggest that supporting N-dimensional arrays
 > should be pretty easy in terms of the API itself. We could also
 > provide fast-paths for 1D and 2D arrays.
 >
 > Of course, actual support for N-dimensional arrays would depend on
 > the underlying implementation.

There are a few array implementations that support N-d arrays, for example
the ones in JHDF5 (https://wiki-bsse.ethz.ch/display/JHDF5) and in netCDF
(http://www.unidata.ucar.edu/software/netcdf-java/).

 > Yep, the idea is to be flexible enough to support many different implementations.

The real problem I see in making all this practically usable is
efficient conversion between implementations. The fragmentation of the
scientific computing ecosystem in Java makes it nearly impossible to
get any real work done without some conversion between array
libraries. If I want to read a matrix from an HDF5 file and then
calculate its eigenvalues, I need to convert from HDF5 arrays to Colt
(or whatever) arrays. Some libraries provide a way to construct an N-d
array from a 1-d data storage area plus a shape vector, without
copying the data. I think it's important to leverage such
functionality in a Clojure API in order to let conversion happen as
much as possible behind the scenes.

I've actually got this working to an extent via a "coercion" protocol that can be extended to different implementations. It's not perfect yet but it works, e.g. you can do

(+ m [[2 0] [0 2]])

And it is smart enough to figure out how to coerce the nested persistent vectors into something that the implementation for matrix m understands.

There's a lot of scope to make this much more efficient of course....

 

BTW, I started my own attempt at something similar a while ago, but
never found the time to get it to a usable state:

   https://code.google.com/p/clj-multiarray/

Konrad.

Very interesting, thanks! I like the use of protocols to construct the specialised array types in particular.

Would you be willing and interested to work on merging this into the "core.matrix" code base? I had started on a simple NDArray implementation myself but your code is clearly much further along than mine.....

[Konrad Hinsen]

Mikera writes:

> BTW, I started my own attempt at something similar a while ago, but
> never found the time to get it to a usable state:
>
> https://code.google.com/p/clj-multiarray/
>
> Konrad.
>
> Very interesting, thanks! I like the use of protocols to construct the specialised array
> types in particular.
>
> Would you be willing and interested to work on merging this into the "core.matrix" code
> base? I had started on a simple NDArray implementation myself but your code is clearly
> much further along than mine.....

Sure! I won't be able to work much on this over the next few weeks, but feel free
to take over any code that looks good to you.

See you on Numerical Clojure...

Konrad.