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