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

Wed Apr 4 16:42:07 PDT 2012

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
> >>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.puremagic.com/pipermail/digitalmars-d/attachments/20120405/2095e9fd/attachment-0001.html>