GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)

[Michael Chen]

Btw, I really don't like the matrix api to be member functions. It is
hard for user to extend the library in an unified way. It is also ugly
when you want to chain the function calls.

On Wed, Apr 4, 2012 at 5:09 PM, Michael Chen <sth4nth at gmail.com> wrote:
> For the Point 4, I really like to have high order functions like
> reduceRow and reduceCol. then the function argument is simply the
> reduceRow!foo(0,mat), here the foo is not a function operating on the
> whole column but simply a function of two elements (e.g.
> reduceRow!("a+b")(0,mat)). Or even better we could have a general
> reduce function with row and col as template parameter so that we can
> do reduce!(foo,row)(0,mat). I dont know whether we can optimize such
> reduce function for different matrix type, but such a function will be
> extremely useful from a user perspective.
>
> On Tue, Apr 3, 2012 at 7:20 PM, Cristi Cobzarenco
> <cristi.cobzarenco at gmail.com> wrote:
>>
>>
>> On 3 April 2012 02:37, Caligo <iteronvexor at gmail.com> wrote:
>>>
>>> I've read **Proposed Changes and Additions**, and I would like to
>>> comment and ask a few questions if that's okay.  BTW, I've used Eigen
>>> a lot and I see some similarities here, but a direct rewrite may not
>>> be the best thing because D > C++.
>>>
>>> 2.  Change the matrix & vector types, adding fixed-sized matrix
>>> support in the process.
>>>
>>> This is a step in the right direction I think, and by that I'm talking
>>> about the decision to remove the difference between a Vector and a
>>> Matrix.  Also, fixed-size matrices are also a must.  There is
>>> compile-time optimization that you won't be able to do for
>>> dynamic-size matrices.
>>>
>>>
>>> 3. Add value arrays (or numeric arrays, we can come up with a good name).
>>>
>>> I really don't see the point for these.  We have the built-in arrays
>>> and the one in Phobos (which will get even better soon).
>>
>>
>> The point of these is to have light-weight element wise operation support.
>> It's true that in theory the built-in arrays do this. However, this library
>> is built on top BLAS/LAPACK, which means operations on large matrices will
>> be faster than on D arrays. Also, as far as I know, D doesn't support
>> allocating dynamic 2-D arrays (as in not arrays of arrays), not to mention
>> 2-D slicing (keeping track of leading dimension).
>> Also I'm not sure how a case like this will be compiled, it may or may not
>> allocate a temporary:
>>
>> a[] = b[] * c[] + d[] * 2.0;
>>
>> The expression templates in SciD mean there will be no temporary allocation
>> in this call.
>>
>>>
>>>
>>> 4. Add reductions, partial reductions, and broadcasting for matrices and
>>> arrays.
>>>
>>> This one is similar to what we have in Eigen, but I don't understand
>>> why the operations are member functions (even in Eigen).  I much
>>> rather have something like this:
>>>
>>> rowwise!sum(mat);
>>>
>>> Also, that way the users can use their own custom functions with much
>>> ease.
>>
>>
>> There is a problem with this design. You want each matrix type (be it
>> general, triangular, sparse or even an expression node)  do be able to
>> define its own implementation of sum: calling the right BLAS function and
>> making whatever specific optimisations they can. Since D doesn't have
>> argument dependent look-up (ADL), users can't provide specialisations for
>> their own types. The same arguments apply to rowwise() and columnwise()
>> which will return proxies specific to the matrix type. You could do
>> something like this, in principle:
>>
>> auto sum( T )( T mat ) {
>>   return mat.sum();
>> }
>>
>> And if we want that we can add it, but this will provide no addition in
>> extensibility. By the way, you can use std.algorithm with matrices since
>> they offer range functionality, but it will be much slower to use
>> reduce!mySumFunction(mat) than mat.sum() which uses a BLAS backend.
>>
>>
>>>
>>>
>>>
>>> 6. Add support for interoperation with D built-in arrays (or pointers).
>>>
>>> So I take that Matrix is not a sub-type?  why? If we have something like
>>> this:
>>>
>>> struct Matrix(Real, size_t row, size_t col) {
>>>
>>>  Real[row*col] data;
>>>  alias data this;
>>> }
>>>
>>> then we wouldn't need any kind of interoperation with built-in arrays,
>>> would we?  I think this would save us a lot of headache.
>>>
>>> That's just me and I could be wrong.
>>
>>
>> Inter-operation referred more to having a matrix object wrapping a pointer
>> to an already available piece of memory - maybe allocated through a region
>> allocator, maybe resulting from some other library. This means we need to
>> take care of different strides and different storage orders which cannot be
>> handled by built-in arrays. Right now, matrices wrap ref-counted
>> copy-on-write array types (ArrayData in the current code) - we decided last
>> year that we don't want to use the garbage collector, because of its current
>> issues. Also I would prefer not using the same Matrix type for
>> pointer-wrappers and normal matrices because the former must have reference
>> semantics while the latter have value semantics. I think it would be
>> confusing if some matrices would copy their data and some would share
>> memory.
>>
>>>
>>>
>>> I've got to tell you though, I'm very excited about this project and
>>> I'll be watching it closely.
>>>
>>> cheers.
>>
>> ---
>> Cristi Cobzarenco
>> BSc in Artificial Intelligence and Computer Science
>> University of Edinburgh
>> Profile: http://www.google.com/profiles/cristi.cobzarenco
>>