std.linalg

H. S. Teoh hsteoh at quickfur.ath.cx
Fri Oct 11 15:27:36 PDT 2013


On Fri, Oct 11, 2013 at 09:49:22PM +0200, FreeSlave wrote:
> 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?

I seem to have lost the link to Denis' code, but you can see the docs
here:

	http://denis-sh.bitbucket.org/unstandard/unstd.multidimensionalarray.html

I haven't put my code up in any public repo (yet), because it's not
quite ready for public consumption. But you can see the docs (somewhat
older version) here:

	http://eusebeia.dyndns.org/~hsteoh/tmp/doc/mda

I'm not 100% satisfied with its API design (spanningIndices probably
introduces an unacceptable bottleneck in many cases), though it does
employ the principle of generalizing D array slices, so all instances of
Arrays are "views" into some underlying data. Assignment, subarrays,
etc., are therefore cheap and have reference semantics. Subarrays can be
employed to allow block-initialization of arrays, etc.. A .dup method is
provided to make actual copies of data when needed.

There are some ugly hacks (IndexRange) in order to work around language
limitations until Kenji's pull for implementing multidimensional slicing
is committed. My code is also not as flexible as Denis' when it comes to
index reordering, but it does work in the 2D case.

Some of the unittests successfully demonstrate that the code is generic
enough to handle arrays of non-numeric types, such as arrays of strings,
in which case per-element unary/binary operations work as expected
(e.g., A ~ "a" appends "a" to all elements, when A is an array of
strings).


[...]
> >>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.

We should simply have standard APIs for checking for compile-time
available properties, similar to std.range.hasLength. Then algorithms
that care about the distinction can make use of it, and other algorithms
that don't care will automatically work with both variants.


[...]
> >>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.

Sorry, I misunderstood what you meant. I thought you were talking about
storage order, but apparently you were talking about index order in the
API, which is a different issue.

My personal preference is columnIndex first, then rowIndex, but that is
opposite of standard matrix notation in math (which is row-first then
column). I think we should just choose one way or the other, and just
stick with it consistently. The important thing is to be consistent to
prevent bugs; exactly which order is the "right" way is just
bikeshedding, I think.

There's another related point here, related to declaring array
dimensions. In my multi-dim array code, I chose to have the index order
match the order of the specified array dimensions; e.g., if you declare:

	auto A = Array!(2,int)([2, 5]);

then you'd index A from A[0,0] to A[1,4]. However, this is opposite of
what nested arrays in D does:

	int[2][5] B;

Here, the index range for B goes from B[0][0] to B[4][1], because B is a
5-element array of 2-element subarrays.

Again, I'm not sure what the "correct" way is; I chose my way because
it's easier to remember (the order you specify dimensions is the same
order the indices appear in, you don't have to keep remembering when to
reverse the order). But either way, a single standard approach should be
chosen, and we should stick with it consistently.


[...]
> 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)

<bikeshed> What about a nicer, more concise name, like is2DArray? :)
</bikeshed>


> {
>     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

I'm uneasy about this one. This is too specific to 2D arrays, and would
not work with general multidimensional arrays (that just happen to be
2D). I think it's better to use opDollar:

	size_t rows = ct.opDollar!0; // assume row-major indexing
	size_t cols = ct.opDollar!1;

Requiring user types to implement opDollar has the advantage that then
we can provide nice $ notation inside the slicing operator:

	auto topLeftElem = myMatrix[0,0];
	auto bottomRightElem = myMatrix[$-1, $-1];


[...]
> 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.

I think we should use opDollar instead of arbitrarily-named custom
fields like rowNum or columnNum. Now, opDollar doesn't quite support
built-in nested arrays, so I'm not sure how to fix that; perhaps provide
some wrappers in std.array? It's kinda annoying, though, because nested
arrays may be jagged, in which case it cannot be treated as a matrix.
But rejecting all nested arrays outright seems a bit too limited. Maybe
provide a wrapper to convert nested arrays into "real" 2D arrays (with
runtime checking against jaggedness)?

That also lets us check if a particular dimension is compile-time known
or not:

	template isStaticDim(Array, int dim) {
		enum isStaticDim = is(typeof({
			// If it's possible to put opDollar!dim into an
			// enum, then it's statically-known, otherwise
			// it's dynamic.
			enum length = Array.init.opDollar!dim;
		}));
	}
	StaticMatrix!(4,4) sm; // both dimensions known at compile-time
	DynamicMatrix dm(4,4); // both dimensions specified at runtime
	RowAppendableMatrix!4 am; // one dimension specified at runtime

	assert(isStaticDim!(typeof(sm), 0));
	assert(isStaticDim!(typeof(sm), 1));

	assert(!isStaticDim!(typeof(dm), 0));
	assert(!isStaticDim!(typeof(dm), 1));

	assert(isStaticDim!(typeof(am), 0));
	assert(!isStaticDim!(typeof(am), 1));


> 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.

I think it should be OK to standardize on ".dup", by analogy with
built-in arrays. If no such method is provided, then do
element-by-element copy.


> 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).

Looks like a hierarchy of concepts:

	           Matrix
		  /      \
	StaticMatrix     DynamicMatrix

The Matrix concept provides the common core APIs for accessing matrices;
the more specific concepts StaticMatrix and DynamicMatrix build on top
of Matrix by adding more API methods/properties specific to static or
dynamic matrices.

If an algorithm doesn't care which kind of matrix is passed to it, it
can simply take Matrix (via an isMatrix signature constraint, ala
std.algorithm.* specifying isInputRange, etc.). If it wants a specific
kind of matrix, it specifies that instead (e.g. isStaticMatrix). If it
can handle both, it can use static if to treat each case separately:

	auto myAlgo(M)(M matrix)
		if (isMatrix!M) // accept most general form of matrix
	{
		... // code
		static if (isStaticMatrix!M)
			... // stuff specific to static matrices
		else static if (isDynamicMatrix!M)
			... // stuff specific to dynamic matrices
		else static assert(0); // oops
		... // more code
	}


T

-- 
VI = Visual Irritation


More information about the Digitalmars-d mailing list