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.

Good idea since there is so much implementations of "fixed-width vectors/matrices" since the beginning of times. There has been such efforts in the past to makes it more standardized.

There will be

On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:
> Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others).

And integers for eg. video processing. Thankfully it's not that much harder.

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

About 3D graphics: I don't think choosing one or the other is that big a deal. I tend to prefer row major order but eg. gl3n uses the reverse.

It's up to you. If you templated against order, you will be able to implement you "asTransposed" idea.

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

Possible with pointer cast? (strict aliasing rule aside)

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

The product to me is a special case.
Element-wise matrix multiplication is seldom used, it's a lot more common to want normal matrix multiplication. I and others use "*" for that purpose.
https://github.com/Dav1dde/gl3n/blob/master/gl3n/linalg.d#L1734
https://github.com/p0nce/gfm/blob/master/gfm/math/matrix.d#L150
https://github.com/gecko0307/dlib/blob/master/dlib/math/matrix.d#L171

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

But it's dissimilar to shader languages and existing D libraries.

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

Use enforce to crash. Dividing by zero is not recoverable.

On Friday, 11 October 2013 at 17:08:26 UTC, ponce wrote:
> Good idea since there is so much implementations of "fixed-width vectors/matrices" since the beginning of times. There has been such efforts in the past to makes it more standardized.
>
> There will be
>
> On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:
>> Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others).
>
> And integers for eg. video processing. Thankfully it's not that much harder.
>
>> 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.
>
> About 3D graphics: I don't think choosing one or the other is that big a deal. I tend to prefer row major order but eg. gl3n uses the reverse.
>
> It's up to you. If you templated against order, you will be able to implement you "asTransposed" idea.
>
>> 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.
>
> Possible with pointer cast? (strict aliasing rule aside)

Well, pointer cast would work for one-dimensional arrays, but not for dynamic two-dimensional matrices. Anyway it's just implementation specific, not design.

>> 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.
>
> The product to me is a special case.
> Element-wise matrix multiplication is seldom used, it's a lot more common to want normal matrix multiplication. I and others use "*" for that purpose.
> https://github.com/Dav1dde/gl3n/blob/master/gl3n/linalg.d#L1734
> https://github.com/p0nce/gfm/blob/master/gfm/math/matrix.d#L150
> https://github.com/gecko0307/dlib/blob/master/dlib/math/matrix.d#L171
>> 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.
>
> But it's dissimilar to shader languages and existing D libraries.

That's actual a common issue to use operator overloading for multiplication or not. Some people does not like it, others prefer this way. That's why I want to discuss it. "*" is more practical, but user always expects "/" as opposite. Matrix has no division operation (except of dividing by scalar), only multiplication by invertible matrix.

> Invertible matrix. It must not be allowed for square StaticMatrix.
I meant non-square of course. Sorry for this mistake. I don't know how to edit my posts here.



The another issue is getting of minors. It can use some operator overload magic and look like

matrix(start, middle1, middle2, end)(start, middle1, middle2, end)

where the first () stands for rows and the second () stands for columns. Gap is between middle1 and middle2. There is efficiency issue - should it make copy or work like some kind of two-dimensional range?

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.

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.

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.


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


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

Here I'd just like to say that inversion is specific to linear algebra, and doesn't apply to multidimensional arrays in general. So it should be implemented in the Matrix wrapper, not tied to any specific data representation.

As for trying to invert non-invertible matrices: it's just another form of "division by zero". I'd say throwing an exception is the right thing to do here.

Of course, you might want to provide an alternative API that checks for invertibility when the caller is not sure whether it can be inverted: since this check pretty much already computes the inverse if it exists, it's probably more efficient to have a single method of doing both. Since it's probably cheaper to do this inversion / check in-place, you could just have a method that takes a Matrix by ref:

	struct Matrix(T) {
		...
		// Returns true if m is invertible, in which case m will
		// be replaced with its inverse; otherwise returns
		// false, and the contents of m will be scrambled. Make
		// a copy of m if you wish the preserve its original
		// value.
		static bool invert(Matrix m) { ... }
		...
	}

This is a sort of low-level operation for performance-sensitive code, so the "nicer" invert operation can be implemented in terms of this one:

	struct Matrix(T) {
		...
		static bool invert(Matrix m) { ... }
		...
		@property Matrix inverse() {
			auto result = this.dup; // make a copy first
			if (!invert(result))
				throw new Exception("Tried to invert non-invertible matrix");
			return result;
		}
	}

	Matrix!(Array!(2,float)) A = ...;
	auto B = A.inverse; // this returns a separate inverse matrix
	A.invert(A);	// this modifies A in-place to become its inverse

Of course, these method names are tentative; I'm sure you can think of better names for them.


T

-- 
Computers shouldn't beep through the keyhole.

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

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

First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?

On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:
> First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?

I picked it up immediately. It's an obvious abbreviation of "linear algebra".

I'm not sure I understand this resentment against abbreviated names.  By your argument, we should rename std.xml to std.extensibleMarkupLanguage instead. Bring on the rainbow, the bikeshed is up for painting.


T

-- 
You are only young once, but you can stay immature indefinitely. -- azephrahel

On Fri, 11 Oct 2013 15:42:17 -0700, H. S. Teoh wrote:

> On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:
>> First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?
> 
> I picked it up immediately. It's an obvious abbreviation of "linear algebra".
> 
> I'm not sure I understand this resentment against abbreviated names.  By your argument, we should rename std.xml to std.extensibleMarkupLanguage instead. Bring on the rainbow, the bikeshed is up for painting.
> 
> 
> T

Don't be silly.  It'd obviously have to be "standardLibrary.extensibleMarkupLanguage"

On 10/11/13 3:42 PM, H. S. Teoh wrote:
> On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:
>> First thing. I see std.linarg I have no clue whatsoever what it even
>> may be about. Do we really want to have a standard lib with names
>> that look like base64 encoded ?
>
> I picked it up immediately. It's an obvious abbreviation of "linear
> algebra".

Honestly I first thought it was "linear-time algorithms".

Andrei