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

[Michael Chen]

Thanks for the explanation, now I get it. In case you are interested,
there is excellent article about monad style c++ template meta
programming by Bartosz Milewski which might be helpful for compile
time optimization for evaluation order.
Really looking forward to official release of the SciD. D is really
suitable for scientific computation. It will be great to have an
efficient and easy-to-use linear algebra library.


On Thu, Apr 5, 2012 at 7:42 AM, Cristi Cobzarenco
<cristi.cobzarenco at gmail.com> wrote:
> Thanks for the feedback!
>
> On 4 April 2012 10:21, Michael Chen <sth4nth at gmail.com> wrote:
>>
>> 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.
>
>
> I'll look into it, thanks.
>
>>
>> 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
>
>
> Well, as I said before, there's nothing stopping us from providing the free
> functions that call the member functionsl. However there is something
> stopping us from not providing a member function alternative. D's lack of
> ADL means we would not allow easy extensibility of possible operations. Say
> Alice invents a new kind of matrix type, the DiagonalMatrix type which
> stores its elements in a 1D array (this isn't exactly how it would work with
> the design we have, she would in fact have to define a storage type, but
> bear with me). If she wants to implement sum on her matrix, she can't simply
> add a specialisation to the sum( T ) function because of the lack of ADL. If
> instead we implemented sum(T) as:
>
> auto sum( T )( T matrix ) {
>      static if( is( typeof( T.init.sum() ) ) )
>          return matrix.sum;
>      else
>          return reduce!"a+b"( 0, matrix.elements() );
> }
>
> then Alice could simply define a DiagonalMatrix.sum() method that wouldn't
> need to go through all the zero elements, for example. The idea of rowReduce
> and columnReduce still works - we can provide this for when a user wants to
> use her own operation for reducing. We could also have a way of
> automatically optimising rowReduce!"a+b"( mat ) by calling
> mat.rowwise().sum() if the operation is available - but that wouldn't be
> exactly high priority.
>
>>
>> .
>> >>
>> >> 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
>> >>>
>
>