What is the status of symbolic matrices?

[rbuchmann]

I saw you have a dedicated repository, but since its documentation still is a little sparse, I am not sure if that's still the preferred way or even how to add the implementation? I would be willing to help with the integration & documentation if there is something I can do.

[mschuene]

Hi,
I am still very interested in symbolic matrices using expresso.
See this old thread on the clojure numerics mailing list for reference
https://groups.google.com/forum/#!msg/numerical-clojure/uIfIAcIufvE/n7CYfVJn0mEJ

There are two ways to support symbolic matrices (or complex ones etc)
1.) make a wrapper type around a matrix implementation and implement the core.matrix
protocols to use the symbolic operations, so you could write
(clojure.core.matrix/add (symbolic [[1 2][3 4]]) [[5 6][7 8]])) and it returns a symbolic matrix.
The other direction is more problematic
(clojure.core.matrix/add [[5 6][7 8]](symbolic [[1 2][3 4]])) wouldn't work as expected.

the symbolic wrapper could also wrap an arbitrary core.matrix implementation that allows for Object elements and this way for the cost of another indirection all existing core.matrix implementation can be leveraged.

2.) the symbolic matrices repository tested another aproach:
To support an arbitrary scalar type it should be sufficiend (and seems natural to do in functional languages) to provide the corresponding scalar-functions. One could then define a generic api that accepts for each function such a specialisation as first argument and uses the functions specified there in it's operation. This quite nicely solves the order issue from the approach above but I didn't fully test the performance implications. The current(old) status is just to give a proof of concept. It supports symbolic matrices with the following specialization
(def symbolic (map->Specialisation
{:add (fn [& args](symb/simplify %28con/cev '+ args%29%29%29
:mul %28fn [& args] %28symb/simplify %28con/cev '* args)))
:sub (fn [& args](symb/simplify %28con/cev '- args%29%29%29
:div %28fn [& args] %28symb/simplify %28con/cev 'clojure.core// args)))
:abs (fn [x](if %28number? x%29 %28Math/abs %28double x%29%29
%28symb/ex' %28abs x%29%29))
:scalar? symb/expression?
:zero 0
:one 1
:= (fn [a b](if %28and %28number? a%29 %28number? b%29%29
%28== a b%29 %28= a b%29))
:supports-inequality? false
:sqrt (fn [& args] (symb/simplify (con/cev 'sqrt args)))
}))

I don't have enough time currently to give serious implementing symbolic matrices a try myself, but I would gladly support you (or anyone willing to do this) with guidance on using expresso / fixing bugs that pop up :)

I haven't looked in detail at the recent try to support complex numbers in core.matrix,but it looks promising.
https://github.com/mikera/core.matrix.complex

and the discussion on the clojure numerics mailing list
https://groups.google.com/forum/#!topic/numerical-clojure/gocQu7XTaNo

[rbuchmann]

I see, thanks for the detailed overview and the links. As far as I can tell the thread didn't really conclude with a definite answer, but at least it mentioned some things to consider.

The implementation for complex matrices is a good reference. I find it interesting that in the thread you linked @mikera seems to argue in favor of your approach 2., but implemented complex matrices as a separate class anyway. I can see why this approach is compelling in this case though, since splitting the matrices into real and complex part allows reusing all other implementations.

I would very much like to try the second direction you mentioned, since I also believe that it is the most natural way to do it and other systems seem to fare quite well with it, e.g. sage's symbolic rings.

I would define protocols for rings or fields, and then define an abstract matrix class parameterized by its underlying ring. I hope that is what you had in mind :)

This might also provide an elegant (albeit partial and potentially slow) solution for the AB vs. BA problem: You could just define one function per pair that exhibits the underlying ring of A as a subring of the one underlying B and vice versa, very much limiting duplication.

[mikera]

FWIW I think that different use cases require different implementations. There is a big tradeoff between making things very generic through abstraction and making them efficient.

https://github.com/mikera/core.matrix.complex is meant to be fairly lightweight and fast. In order to do that, the implementation is deliberately specialised to complex values only.

I'd fully support experiments to make more generic implementations, I just think that these would have quite different design goals. They would probably be of more interest in academic context, whereas the more specialised implementations are likely to appeal to industry practitioners.

The nice thing about core.matrix is that these implementations can inter-operate.... so it isn't really an either/or decision, we can have both.

[rbuchmann]

Sure, there is always a tradeoff, even with macros. Like I said, I see why this particular implementation is very appealing. If the parameterized arrays work as expected, I can probably provide the other end of the spectrum pretty easily.

[rbuchmann]

Ok, so I played around with it a bit, but I ran into a problem that I could use your help with: When I implement the arrays on top of the most general object container arrays, having something like (* 5 a) in them messes up the built in dimensionality checkers. Do I need to put every entry in a wrapper type, or is there a smarter solution?

I like how it is handled in core.matrix.complex, where the explicit element type is only constructed when needed, as for example in the get-FOO-d functions.

[mschuene]

Hmm, I think that is a current limitation of core.matrix because it regads everything seqable as an array and not a scalar. Is this correct @mikera ?
One would need a is-scalar? protocol for implementations, which for symbolic matrices would be expression?. This is what I did for the symbolic matrices, where it could be specified in the specialization map.

If this is the case, then I see no better solution as using a wrapper type.

using the wrapper type one still could extract the underlying expression from the wrapper type in the get-FOO-d functions.

[mikera]

Yeah, we probably shouldn't interpret arbitrary sequences as arrays. I've been playing around with this a bit, it gets tricky because of the number of permutations. But I think we will revert to considering plain sequences as scalar values (except in a few limited cases like array construction, perhaps...)

We need to consider vectors as arrays though - we have to make them count as arrays or else the :persistent-vector implementation can't work.

My suggested workaround is that if you want to put sequences / vectors as scalar elements then you should use an implementation that explicitly stores the array shape (e.g. NDArray). Then you avoid any ambiguity.

[rbuchmann]

Hm, I would like to keep the implementation as general as possible, but I think using NDArray and specifying the dimensions is the most elegant way to do it here (Apart from core.matrix dropping the support for sequences as arrays of course). Since my approach is supposed to support arbitrary rings, and since construct-matrix should be implemented independently of the structure of the underlying ring, it would feel pretty weird to have to do something like

(matrix [[(Fraction. 1/2) (Fraction. 0/1)]...])

Thanks for the suggestions, I'll try using NDArray now.