konovod / linalg

Linear algebra library based on LAPACK
MIT License
50 stars 4 forks source link

`sum` requires a `zero` method if no initial is given #6

Closed KCErb closed 5 years ago

KCErb commented 5 years ago

Example program

require "linalg"

m = LA::GMat32.new(3, 3)
m1 = LA::GMat32.new(3, 3)
[m, m1].sum

throws

Error in line 5: instantiating 'Array(LA::GeneralMatrix(Float32))#sum()'

in /usr/local/Cellar/crystal/0.27.0/src/enumerable.cr:1112: undefined method 'zero' for LA::GeneralMatrix(Float32).class

    sum Reflect(T).first.zero

The idea is that if we don't pass an initial, then the sum has to start at zero (https://crystal-lang.org/api/0.27.2/Enumerable.html#sum-instance-method).

This issue is to discuss a possible solution. I'm guessing that we can overwrite sum to create a zeros to match the size of the first element (and even assert_same_size) when an initial is missing.

If you like that idea, or have a better one, let me know and I can probably put together a PR (unless you'd rather).

konovod commented 5 years ago

Huh, it's sad we can't use sum (and product) from stdlib Enumerable. I was sure it uses first element and adds other to it (like reduce do when initial not provided). I think it's better way than using zeros. There is still a question - what to return if a collection is empty? Return empty matrix? Raise an exception? I'm not sure what is better.

Another idea is to create a special "zero" matrix without size that can be added to any matrix, then return it as LA::Matrix(T).zero implementation. In that case sum doesn't need to be reimplemented, but this looks pretty hacky.

KCErb commented 5 years ago

Yes, I thought of that too and felt the same way, it would be too special for such a use case.

But the more I think about it the more it "feels" ok from a group theoretical perspective. The ability to support

matrix * Sizeless(1) = matrix

and

matrix + Sizeless(0) = matrix

could be nice for the same reason that [].sum == 0 is nice. On the one hand it's a silent failure, but on the other it makes types easier to reason about which is sometimes a good motivator for working with arrays (at least I can say in my experience that if I don't know the cardinality of something, I like putting it into an array and I like that if the array is empty nothing happens).

Then again, this is starting to feel like the equivalent of just supporting numeric + matrix (i.e. element wise addition). How do you feel about that?

....

Ok one more thought here, do you support a type that fixes matrix size (sorry still new to the library)? Like

a = [] of LA::GMat(3,3)

Because sum and product should only be defined for an array where the matrices are all the same size right? In that case it would make sense to me to have our own sum and product and throw an error (at compile time) when calling sum or product on an array of matrices of arb dimension.

konovod commented 5 years ago

Yes, i think the best way is just to add support for numeric + matrix operation (and make Matrix(T).zero returning T.new(0)). After all, this works in Matlab and Numpy. You can do PR if you want or I'll do it tomorrow.

As for matrices of fixed size - I decided to don't support them. This library is built around LAPACK and oriented to large matrix processing while fixed size is mostly useful for small matrices (used in, say, game engines). Library for small matrices should be pretty different - no need in LAPACK, should allocate them on stack, use inline routines and SIMD to achieve maximum speed.

KCErb commented 5 years ago

OK cool, makes sense. I probably wouldn't have time to work on this as quickly as tomorrow. So you can go ahead :)

konovod commented 5 years ago

fixed in https://github.com/konovod/linalg/commit/79c2eee8984c079970bb34cdee871fd42d59aefd