std.linalg

[H. S. Teoh]

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.