Designing a matrix library for D

[H. S. Teoh via Digitalmars-d]

On Mon, Jun 23, 2014 at 04:45:27PM +0000, bearophile via Digitalmars-d wrote:
> I'd like a matrix library for D, able to handle sparse matrices
> too. But first I'd like to discuss its API and general design a
> bit.

When it comes to libraries that need to handle diverse concrete types
(e.g. compact matrices vs. sparse matrices), my take is always to write
the matrix algorithms in a way that's independent of the concrete matrix
type(s), and to provide concrete types as a separate submodule (or
logically distinct section of the API).


[...]
> Those collections of X10 are based on "points" and "regions".
> 
> A Point is a vector of N keys, that are indexes of a vector (like
> the row,col for a matrix, or the key for an associative array).
> If N is a compile-time constant we lose some flexibility but we
> gain efficiency. In theory we want both.
> 
> A Region defines part of a collection. In the simpler case it's a
> slice of a 1D array, or it can be a 2D rectangular region of a 2D
> array, or more. So it's an array of M points. Usually you want M
> to be a compile-time constant, but in some cases you want more
> flexibility. In some cases you want to intersect regions, or you
> want non-convex regions, so the situation gets a little more
> complex.

I think supporting non-convex (or non-rectangular) regions is a bit too
ambitious. Basic support should focus on rectangular/cuboidal/etc.
regions, which have the advantage that they are closed under
intersections (the intersection of two rectangular regions is always
another rectangular region). Non-rectangular regions would destroy a lot
of simplifying assumptions that can be made to streamline the design,
and yet most applications I can think of need only rectangular regions.
I don't think it's worth the price to support such unusual usages.

(Besides, if we design the API correctly, they should still be
supportable via wrapped matrix types -- this is why matric algorithms
should not depend on concrete types, so that for special needs the user
can simply define a wrapper that remaps points in such a way as to
achieve the desired non-convex regions.)


> The paper lists some Region operations, that I think can be
> implemented with no changes to the D language:
> 
> R.rank ::= # dimensions in region;
> R.size() ::= # points in region
> R.contains(P) ::= predicate if region R contains point P
> R.contains(S) ::= predicate if region R contains region S
> R.equal(S) ::= true if region R and S contain same set of points
> R.rank(i) ::= projection of region R on dimension i (a
> one-dimensional region)
> R.rank(i).low() ::= lower bound of i-th dimension of region R
> R.rank(i).high() ::= upper bound of i-th dimension of region R
> R.ordinal(P) ::= ordinal value of point P in region R
> R.coord(N) ::= point in region R with ordinal value = N
> R1 && R2 ::= region intersection (will be rectangular if R1 and
> R2 are rectangular)
> R1 || R2 ::= union of regions R1 and R2 (may or may not be
> rectangular,compact)
> R1 - R2 ::= region difference (may or may not be
> rectangular,compact)

My current thoughts on this, is that concrete matrix/array types should
only support array element access and size measurements. Operations like
intersections, subarrays, etc., should be done via wrapper types that
remap element access and size measurements, but forward element accesses
to the underlying concrete type. This way, we avoid excessive copying,
which get very expensive with large or high-dimensional arrays.


> The paper also shows some ideas for array operations:
> 
> A.rank ::= # dimensions in array
> A.region ::= index region (domain) of array
> A.distribution ::= distribution of array A

What's a "distribution"?


> A[P] ::= element at point P, where P belongs to A.region
> A | R ::= restriction of array onto region R (returns copy of
> subarray)

IMO, we should avoid all implicit copying. Let the user explicitly ask
for a copy if a copy is desired. Implicit copying is the mortal enemy of
high performance.


> A.sum(), A.max() ::= sum/max of elements in array
> A1 <op> A2 ::= returns result of applying a point-wise op on
> array elements,
> when A1.region = A2. region
> (<op> can include +, -, *, and / )
> A1 || A2 ::= disjoint union of arrays A1 and A2
> (A1.region and A2.region must be disjoint)
> A1.overlay(A2) ::= array with region, A1.region || A2.region,
>    with element value A2[P] for all points P in A2.region
>    and A1[P] otherwise.

All of these should be generic functions / wrappers, I think.


[...]
> I've found a second source of ideas in this page:
> http://dis.um.es/~alberto/hmatrix/static.html
[...]
> >The size of the result of certain computations can only be known at
> >run time. For instance, the dimension of the nullspace of matrix
> >depends on its rank, which is a nontrivial property of its elements:<
> 
> The library solves this problem using a GHC compiler language
> extension:
> 
> >By hiding the unknown dimension in an existential type we can still
> >compute safely with a strongly typed nullspace.<
[...]
> With this strategy in most cases the compiler is able to enforce
> the arrays are of right types at compile-time, and allows
> run-time sizes for the few remaining cases.
> 
> In D this idea should avoid template bloat caused by all the
> different array dimensions, but allow optimizations based on the
> compile-time knowledge of sizes in the points where the
> programmer asks for more efficiency (like unrolling on request in
> certain zones if the size is just 4).
> 
> Is all this possible with the type system of D? Perhaps it's
> possible with the enum preconditions or some of the stuff
> discussed in the "Value range propagation for if-else" thread.
> Improving the D type system could be an option, if useful.
[...]

IMO, if we separate matrix algorithms from concrete matrix types, we've
already won half the battle.

Some applications already know the exact sizes of the matrices they will
use (e.g., 3D graphics apps / games will almost certainly use only 4x4
homogeneous matrices), so it doesn't make sense to force them to pay for
the runtime penalty of dynamically-sized m*n matrices. So in this case,
they would use matrices whose size is compile-time determined, probably
as compile-time size parameters to a StaticMatrix template.

But other applications, e.g., computer algebra systems, cannot know the
size of matrices they will need until runtime. For these applications,
compile-time matrix sizes can't be used, so they will need to work with
dynamically-sized matrices.

Furthermore, most applications only use matrices of a fixed rank
(usually rank 2, rectangular/square matrices), so it makes sense to take
rank as a compile-time parameter. But some applications (e.g., programs
that perform tensor calculations) need to support runtime-determined
ranks. But if we make rank a runtime parameter, the previous two
categories of applications will suffer unacceptable performance
degradation.

But if we design a proper generic matrix concept, where *any* concrete
type can be used by a given matrix algorithm provided it satisfies some
required properties (analogous to input ranges satisfying the range
API), then we can solve this problem by having multiple matrix modules
that provide different concrete types (StaticMatrix, Matrix (of fixed
rank), Tensor (dynamic rank), SparseMatrix, etc.), but we only need to
implement a single set of matrix operations / algorithms that can work
with any combination of concrete matrix types. (And of course, these
algorithms can specialize on certain concrete types if it helps improve
performance, but at least they should be able to handle anything that
looks like a matrix.)

The key to success, then, lies in how we design the generic matrix API.
If we do it right, then it should be easy to implement features like
compile-time size verifications, etc.. Otherwise it can be very painful.


T

-- 
Life is complex. It consists of real and imaginary parts. -- YHL