std.linalg

[John Colvin]

On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:
> There is "Matrices and linear algebra" module in wish list. 
> Let's discuss its design. D is complicated language so it's 
> difficult to choose the right way here. We need to find 
> compromise between efficiency and convenient interface. I'm 
> going to make some suggestions how this module should look like.
>
> First of all, it should provide two templates for matrices. 
> Let's call them StaticMatrix and DynamicMatrix. The first one 
> has "templated" size and therefore may use static arrays and 
> compile-time checks. It can be useful when the size is 
> determined by our needs, for example, in graphics. 
> DynamicMatrix has variable size, i.e. it should be created in 
> heap. It can be useful in all other math areas.
>
> Both templates should support all floating point types and 
> moreover user-defined (for example wrappers for GMP library and 
> others).
>
> For efficiency in both cases matrices should use 
> one-dimensional array for inner representation. But actually 
> I'm not sure if matrices should support other container types 
> besides standard D arrays. The good thing about one-dimensional 
> arrays is that they can be easily exposed to foreign functions, 
> for example, to C libraries and OpenGL. So we should take care 
> about memory layout - at least row-major and column-major. I 
> think it can be templated too.
>
> But another question arises - which "majority" should we use in 
> interface? Interface should not depend on inner representation. 
> All functions need unambiguity to avoid complication and 
> repetition of design. Well, actually we can deal with different 
> majority in interface - we can provide something like 
> "asTransposed" adapter, that will be applied by functions if 
> needed, but then we will force user to check majority of matrix 
> interface, it's not very good approach.
>
> Sometimes user takes data from some other source and wants to 
> avoid copying in Matrix construction, but she also wants to get 
> matrix functionality. So we should provide "arrayAsMatrix" 
> adapter, that can adopt one-dimensional and two-dimensional 
> arrays making them feel like matrices. It definitely should not 
> make copy of dynamic array, but I'm not sure about static.
>
> About operation overloading. It's quite clear about 'add' and 
> 'subtruct' operations, but what's about product? Here I think 
> all 'op'-functions should be 'element by element' operations. 
> So we can use all other operations too without ambiguity. For 
> actual matrix multiplication it can provide 'multiply' or 
> 'product' function. It's similar to Maxima approach, besides 
> Maxima uses dot notation for these needs.
>
> Transposition. I've already mentioned "asTransposed" adapter. 
> It should be useful to make matrix feel like transposed without 
> its copying. We also can implement 'transpose' and 'transposed' 
> functions. The first one transposes matrix in place. It's 
> actually not allowed for non-square StaticMatrix since we can't 
> change the size of this type of matrices at runtime. The second 
> one returns copy so it's applicable in all cases. Actually I'm 
> not sure should these functions be member functions or not.
>
> Invertible matrix. It must not be allowed for square 
> StaticMatrix. DynamicMatrix may make checks at runtime. Same as 
> for transposition we can implement 'invert' to invert in place 
> and 'inverted' to make copy. There is issue here - what should 
> we do when determinant is 0? I believe the best approach here 
> is to throw exception since if user needs invertible matrix 
> it's definitely exception when it can't be calculated.
>
> Please, make your suggestions too.

I'd like to echo some of the points made by CJS.

Rolling your own linear algebra (for anything other than very 
small matrices) is an exercise in futility unless you are 
prepared to become a world-class expert in what is an incredibly 
mature and detailed field.

Clever wrapping of industry standard libraries is the way to go.

I would love to see this https://github.com/cristicbz/scid taken 
further, which was a fork of scid that used clever magic to 
reduce matrix expressions to an semi-optimal set of blas calls. 
It doesn't appear to have got much love in recent times.