Generalising to arbitrary scalar types

[Mike Anderson]

Hi Guys,

Longish post to summarize a perspective on the thorny design topic of how to generalise core.matrix to arbitrary scalar types, motivated by the various discussions on symbolic and complex matrices etc. Thanks to Kovas and Maik in particular for the detailed and intelligent exploration of this space so far.

Here are the "requirements" as I see them:

- Allow matrix operations with arbitrary user / implementation defined scalar types (Complex numbers, symbolic expressions etc.)

- Allow generalised matrix / array algorithms using such user defined scalar types

- Do not adversely affect performance for regular numeric matrix operations (99% case)

- Do not create significant additional problems / complexity for core.matrix implementers

Would love thoughts / constructive critique. Once we've thrashed these out, I'll put them into a design doc on the Wiki.

Design ideas follow:

(0) Operations are separated from matrix types (with the important special exception detailed in points 2-4, see below)

core.matrix fundamentally provides an abstraction over array types, in the sense of n-dimensional arrays of values.

In a generic sense, arrays can contain anything. So we don't know what a user might want to do with an Object array for example. Also we don't want to force the creation of a new matrix type just to support a new generic operation (this is hard for implementers, and potentially creates a lot of overhead in terms of additional classes and protocol implementations to compile dynamically etc.)

So we allow a user to mix and match operations with relevant array types. The mathematical / logical domain used is a property of the *function*, not the array. As an example, consider two completely different operations that might be applied to the same arrays containg java.lang.Double values:

(complex-add matrixA matrixB) ;; should work for any two matrices that contain Double values and produce complex results

(symbolic-add matrixA matrixB) ;; should work for any two matrices that contain Double values and produce s-expression results

(1) Generic versions of core.matrix operations in new namespace "clojure.core.matrix.generic"

We create a new API for higher order generic operations that can be parameterised with generic operators. This will enable implementation of algorithms that work on arbitrary scalar types and operations.

Higher order functions can be used to construct generic operations, e.g.

(def my-inner-product (clojure.core.matrix.generic/construct-inner-product my+, my*, my0)

It is anticipated that users wanting specialised operations (e.g. symbolic addition with Clojure s-expressions :-) ) would either:

a) use this API directly, providing their own symbolic operations each time

b) use the higher order functions to build a namespace containing specialised operations (e.g. "expresso.symbolic.matrix" )

Facilities may be provided to support compilation of efficient generic operations for specific matrix types.

(2) Each matrix type declares a "default" numerical scalar type and corresponding operations

(numerical-type m)

=> java.lang.Double

Each matrix can define its own ==, +, *, 0, 1 etc. In 99% of cases (including every implementation in core.matrix itself) these will be the regular numerical definitions. 

Other implementations happy with clojure.core.+, clojure.core.== don't need to do anything special as these will be the defaults.

(3) the functions in the main "core.matrix" namespace will always use the default operations.

I.e. c.c.m/add will always use the default + operation for the given matrix type. Usually Double for numeric stuff, but it *could* be a complex matrix. The complex matrix type would need to re-implement the add protocols using complex addition if it wanted to perform these efficiently.

This is necessary to ensure no performance overhead for the 99% numeric case. It also means that implementations that are happy to use regular numeric operations don't need to do anything differently from the current approach: they just need to implement the core.matrix protocols in the obvious way.

This is logically consistent if you take the API contract of the "clojure.core.matrix" functions to be "perform computation using the array's default numerical operations".

(4) It's the user's responsibility to exercise caution when mixing matrix types with different default operations.

The primary matrix used in core.matrix API functions will determine the underlying operation used. We're pretty consistent with this design principle across core.matrix, and it's the most pragmatic option given how protocol dispatch works.

This means that we need to be *extra careful* about deviating from the regular default numerical operations: whenever we do so this will cause extra work for users.

(add complex-matrix vectorz-matrix) ;; should work: assuming double will get cast into complex

(add vectorz-matrix complex-matrix) ;; error! can't cast complex into double

We can provide cast operations to make explicit conversions easier:

(add vectorz-matrix (numerical complex-matrix)) ;; should work if complex-matrix contains only real values

(add (array :complex vectorz-matrix) complex-matrix) ;; should work by getting everything into the complex domain

[Maik Schünemann]

Gret post, Mike!
I am perfectly happy with the design you posted above

because this api would be - coincidentally :) - a perfect fit for
expresso I would
like to help contributing to the new higher order generic operation
API and build expressos
symbolic matrix manipulations on top of that.

I think this is good for now it is what the user expects from protocol dispatch.
Provided this will be used much by core.matrix users and it is possible
to get julia like method dispatch on both arguments without too much overhead
this may be something to look out for in the long term perspective
> --
> You received this message because you are subscribed to the Google Groups
> "Numerical Clojure" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to numerical-cloj...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Mike Anderson]

Cool, thanks Maik!

Your contributions in this space would be welcome if you want to take it on: the generic operation system is something that clearly needs a lead developer / maintainer (and it probably won't be me: I've got my hands full with other projects right now :-) ).