std.linalg

[FreeSlave]

On Friday, 11 October 2013 at 17:49:32 UTC, H. S. Teoh wrote:
> On Fri, Oct 11, 2013 at 06:10:19PM +0200, 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.
>
> I think we need to differentiate between multidimensional 
> arrays (as a
> data storage type) and linear algebra (operations performed on 
> 2D
> arrays). These two intersect, but they also have areas that are 
> not
> compatible with each other (e.g. matrix product vs. 
> element-by-element
> product). Ideally, we should support both in a clean way.
>
> As far as the former is concerned, Denis has implemented a
> multidimensional array library, and I've independently done the 
> same,
> with a slightly different interface. I think one or two others 
> have
> implemented similar libraries as well. It would be good if we
> standardize the API so that our code can become interoperable.

Can you please give links to both libraries?

> As far as the latter is concerned, I've been meaning to 
> implement a
> double-description convex hull algorithm, but have been too 
> busy to
> actually work on it. This particular algorithm is interesting, 
> because
> it stress-tests (1) performance of D algorithms, and (2) 
> challenges the
> design of matrix APIs because while the input vertices (resp.
> hyperplanes) can be interpreted as a matrix, the algorithm 
> itself also
> needs to permute rows, which means it is most efficient when 
> given an
> array-of-pointers representation, contrary to the usual 
> flattened
> representations (as proposed below). I think there's a place 
> for both,
> which is why we need to distinguish between data representation 
> and the
> algorithms that work on them.
>
>
>> 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.
>
> I like this idea. Ideally, we should have many possible 
> representations,
> but all conforming to a single API understood by all 
> algorithms, so that
> you only have to write algorithms once, and they will work with 
> any data
> structure. That's one key advantage of D, and we should make 
> good use of
> it.

The problem is that algorithms still should know matrix template 
to provide compile-time checks if possible or throw exceptions at 
runtime if something gone wrong.

> I do not want to see a repetition of the C++ situation where 
> there are
> so many different matrix/multidimensional array libraries, and 
> all of
> them use incompatible representations and you cannot freely 
> pass data
> from one to algorithms in the other. Then when no library meets
> precisely what you need, you're forced to reinvent yet another 
> matrix
> class, which is a waste of time.
>
>
>
>> Both templates should support all floating point types and 
>> moreover
>> user-defined (for example wrappers for GMP library and others).
>
> Definitely, yes.
>
>
>> 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.
>
> We should not tie algorithms to specific data representations 
> (concrete
> types). One key advantage of D is that you can write algorithms
> generically, such that they can work with *any* type as long as 
> it
> conforms to a standard API.  One excellent example is the range 
> API:
> *anything* that conforms to the range API can be used with
> std.algorithm, not just a specific representation. In fact, 
> std.range
> provides a whole bunch of different ranges and range wrappers, 
> and all
> of them automatically can be used with std.algorithm, because 
> the code
> in std.algorithm uses only the range API and never (at least in 
> theory
> :P) depends on concrete types. We should take advantage of this 
> feature.
>
> It would be good, of course, to provide some standard, 
> commonly-used
> representations, for example row-major (or column-major) matrix 
> classes
> / structs, etc.. But the algorithms should not directly depend 
> on these
> concrete types. An algorithm that works with a matrix stored as 
> a 1D
> array should also work with a matrix stored as a nested array 
> of arrays,
> as well as a sparse matrix representation that uses some other 
> kind of
> storage mechanism. As long as a type conforms to some standard 
> matrix
> API, it should Just Work(tm) with any std.linalg algorithm.
>
>
>> 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.
>
> Algorithms shouldn't even care what majority the data 
> representation is
> in. It should only access data via the standardized matrix API 
> (whatever
> it is we decide on). The input type should be templated so that 
> *any*
> type that conforms to this API will work.
>
> Of course, for performance-sensitive code, the user should be 
> aware of
> which representations are best-performing, and make sure to 
> pass in the
> appropriate type of representations; but we should not 
> prematurely
> optimize here. Any linear algebra algorithms should be able to 
> work with
> *any* type that conforms to a standard matrix API.
>

I'm not sure if you understand idea of differences between inner 
implementation majority and interface majority. I agree that 
inner majority should be defined by inner type. Interface 
majority is just choice between

matrix[rowIndex, columnIndex]

and

matrix[columnIndex, rowIndex]

In case of interface majority we just must choose the appropriate 
one and use it all over the library. It does not relate to 
performance.

>> 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.
>
> If a function expects a 1xN matrix, we should be able to pass 
> in an
> array and it should Just Work. Manually using adapters should 
> not be
> needed. Of course, standard concrete matrix types provided by 
> the
> library should have ctors / factory methods for initializing a 
> matrix
> object that uses some input array as initial data -- if we 
> design this
> correctly, it should be a cheap operation (the matrix type 
> itself should
> just be a thin wrapper over the array to provide methods that 
> conform to
> the standard matrix API). Then if some function F requires a 
> matrix
> object, we should be able to just create a Matrix instance with 
> our
> input array as initial data, and pass it to F.
>
>
>> 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.
>
> Here is where we see the advantage of separating representation 
> from
> algorithm. Technically, a matrix is not the same thing as a 2D 
> array,
> because a matrix has a specific interpretation in linear 
> algebra,
> whereas a 2D array is just a 2D container of some elements. My
> suggestion would be to write a Matrix struct that wraps around 
> a 2D
> array, and provides / overrides the overloaded operators to 
> have a
> linear algebra interpretation.
>
> So, a 2D array type should have per-element operations, but 
> once wrapped
> in a Matrix struct, it will acquire special matrix algebra 
> operations
> like matrix products, inversion, etc.. In the most general 
> case, a 2D
> array should be a specific instance of a multidimensional 
> array, and a
> Matrix struct should be able to use any underlying 
> representation that
> conforms to a 2D array API. For example:
>
> 	// Example of a generic multidimensional array type
> 	struct Array(int dimension, ElemType) {
> 		...
> 		Array opBinary(string op)(Array x)
> 		{
> 			// implement per-element operations here
> 		}
> 	}
>
> 	// A matrix wrapper around a 2D array type.
> 	struct Matrix(T)
> 		if (is2DArray!T)
> 	{
> 		T representation;
> 		Matrix opBinary(string op)(Matrix x)
> 			if (op == "*")
> 		{
> 			// implement matrix multiplication here
> 		}
>
> 		Matrix opBinary(string op)(Matrix x)
> 			if (op != "*")
> 		{
> 			// forward to representation.opBinary to default
> 			// to per-element operations
> 		}
>
> 		// Provide operations specific to matrices that don't
> 		// exist in general multidimensional arrays.
> 		Matrix invert() {
> 			...
> 		}
> 	}
>
> 	Array!(2,float) myArray, myOtherArray;
> 	auto arrayProd = myArray * myOtherArray; // per-element 
> multiplication
>
> 	auto A = Matrix(myArray);	// wrap array in Matrix wrapper
> 	auto B = Matrix(myOtherArray);
> 	auto C = A * B;			// matrix product
>
> The idea of the Matrix struct here is that the user should be 
> free to
> choose any underlying matrix representation: a 1D array in 
> row-major or
> column-major representation, or a nested array of arrays, or a 
> sparse
> array with some other kind of representation. As long as they 
> provide a
> standard way of accessing array elements, Matrix should be able 
> to
> accept them, and provide matrix algebra semantics for them.
>
>
>> 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.
>
> The most generic approach to transposition is simply a 
> reordering of
> indices. This difference is important once you get to 3D arrays 
> and
> beyond, because then there is no unique transpose, but any 
> permutation
> of array indices should be permissible. Denis' multidimensional 
> arrays
> have a method that does O(1) reordering of array indices: 
> basically, you
> create a "view" of the original array that has its indices 
> swapped
> around. So there is no data copying; it's just a different 
> "view" into
> the same underlying data.
>
> This approach of using "views" rather than copying data allows 
> for O(1)
> submatrix extraction: if you have a 50x50 matrix, then you can 
> take
> arbitrary 10x10 submatrices of it without needing to copy any 
> of the
> data, which would be very expensive. Avoiding unnecessary 
> copying
> becomes very important when the dimension of the array 
> increases; if you
> have a 3D or 5D array, copying subarrays become extremely 
> expensive very
> quickly.
>
> A .dup method should be provided in the cases where you 
> actually *want*
> to copy the data, of course.
>
> Basically, subarrays / transpositions / index reordering should 
> be
> regarded as generalizations of D's array slices. No data should 
> be
> copied until necessary.
>
>
>> Invertible matrix. It must not be allowed for square 
>> StaticMatrix.
>
> You mean for non-square StaticMatrix?

Yes, non-square. My bad.



Well, ok. We want to abstract from inner representation to 
provide freedom for users. We fall in metaprogramming and generic 
programming here, so we need to define concepts just like 
Boost/STL/std.range do. The good thing is that in D types with 
different interfaces and syntax constraints can satisfy same 
concept that would be impossible or very difficult in C++. Thanks 
to static if and is(typeof()). For example inner representation 
type can provide [][] operator or [,] operator and Matrix type 
will understand both cases.

Suppose:

template canBeMatrixRepresentation(T)
{
     enum bool canBeMatrixRepresentation = is(typeof(
         {
             T t; //default constructable
             const(T) ct;
             alias ElementType!T E; //has ElementType
             E e; //element type is default constructable
             static if (/*has [,] operator*/)
             {
                 t[0,0] = e; //can be assigned
                 e = ct[0,0]; //can retrive element value from 
const(T)
             }
             else static if (/*has [][] operator*/)
             {
                 t[0][0] = e; //can be assigned
                 e = ct[0][0]; //can retrive element value from 
const(T)
             }
             else
             {
                 static assert(false);
             }

             size_t rows = ct.rowNum; //has row number
             size_t cols = ct.columnNum; //has column number

             t.rowNum = size_t.init;
             t.columnNum = size_t.init;
         }));
}

We see that two-dimensional D array does not satisfy this concept 
because it has no rowNum and columnNum so it should be handled 
separately. This concept is not ideal since not all types may 
provide variable rowNum and columnNum. Also concept should expose 
information whether is "static" type or not, so algorithms will 
know can they use compile-time checks or not.
Types also can provide copy constructor. If they do then Matrix 
will use it, if they don't then Matrix will do element-by-element 
copy. It also can try .dup property.

It's just example how it should work, but I hope the point is 
clear. We also need Matrix concept (or separate concepts for 
StaticMatrix and DynamicMatrix).