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

[Michael Chen]

another btw, there is also another great c++ linear algebra library
besides Eigen: Amardillo which has very simple matlab like interface
 and performs on a par with Eigen.
On Wed, Apr 4, 2012 at 5:14 PM, Michael Chen <sth4nth at gmail.com> wrote:
> 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
>>>