Re: matrix library ?

From: lums_at_[hidden]
Date: 2001-03-24 07:50:35

--- In boost_at_y..., hankel_o_fung_at_y... wrote:
> --- In boost_at_y..., "David Abrahams" <abrahams_at_m...> wrote:
> > > IMHO, it is natural to use both index styles in a single 
problem,
> > 
> > (!) How does /that/ situation come up?
> A typical scenario is that you want to index n+1 increasing time
> points as well as the n time intervals (or their lengths) bracketed
> by these points. In some problem domains (e.g., financial option
> pricing), it is natural to use 1-based indices for time periods
> and 0-based indices for the time points.

Our experience with developing (and using) MTL, interestingly enough, 
is that once you start writing your algorithms in a generic 
programming style, the fact that there is any base index at all 
becomes surprisingly irrelevant.  There are cases where one might 
have to expose the internal representations of a matrix (say, in a 
compressed sparse matrix) when dealing with external code (importing 
or exporting the matrix data), but to a very large extent, you just 
use iterators and the notion of an index having a value vanishes.

Here is an illustrative example.

template <class Matrix, class VecX, class VecZ>
void mult(const Matrix& A, const VecX& x, VecZ z)
{
  typename Matrix::const_iterator i;
  typename Matrix::OneD::const_iterator j;

  for (i = A.begin(); i != A.end(); ++i)
    for (j = (*i).begin(); j != (*i).end(); ++j)
      z[row(j)] += *j * x[column(j)];
}

This a matrix-vector product written in a generic style.  Note what 
is and what is not going on.  There is no assumption about the 
storage format of the matrix -- this will work for dense or sparse, 
row major or column major, triangular, banded, etc.  There is also no 
assumption of indexing.  All that is required is that the vector 
classes have a consistent domain/range relationship via the linear 
transformation represented by the matrix.  The vectors (at our 
mathematical layer of abstraction) are essentially tuples and we use 
some information from the linear operator to get to the element of 
interest within the tuple.  That we used the terms 'row' and 'column' 
to denote this is a throwback to the old way of thinking about this 
and somewhat a red herring.

Disclaimer: We are also somewhat mixing layers of abstraction in this 
example.  The iterator i for instance is iterating over a container 
of containers.  This really belongs more at the storage 
representation layer.  A more pedagogically correct implementation 
might have been to first get a handle on the storage and then jump 
into the loops:

  for (i = A.storage.begin(); i != A.storage.end(); ++i)
    for (j = (*i).begin(); j != (*i).end(); ++j)
      z[row(j)] += *j * x[column(j)];

(Just a random thought)

In a more mature version of such a routine there are also concept 
checks that one might want to add to assert, e.g., that the Matrix is 
a transformation from the one Vector Space to the other (and, one 
would also want checks that the vector quantities being passed in 
were elements of a Vector Space).

 
> > > and a good array library should support both styles.
> > > A matrix library, however, can only and should at least support
> > > the 1-based index, because matrix indices are 1-based in math.
> > 
> > Are you using some subtle semantic distinction between "array"
> > and "matrix" here?
> In terms of arrangement of quantities or symbols, a matrix is
> always an array. 

Conceptually, yes.  Since all finite dimensional vector spaces are 
isomorphic to linear arrangements of numbers (tuples), all linear 
operators mapping from one finite dimensional vector space to another 
are isomorphic to a collection of numbers arranged in tabular 
fashion -- "array" form.  Representing a matrix in a computer program 
as an actual array -- meaning there is space allocated for all 
possible entries in such a tabular form -- is only one way of 
representing a matrix (usually called "dense") and there are even 
variations within that (row-major, column-major, e.g.).

cf also a previous exchange I had with Dave Abrahams on this list 
about "matrix-free" linear operators (probably better-named "array-
free" operators or "storage-free" operators -- but "matrix-free" is 
the name used in the literature).

> > > Yet, even if the index convention problem is fixed, we still 
need
> > > to solve the row-major/column-major problem ...
> > 
> > What problem is that? I thought row-major/column-major was just a
> > distinction describing underlying storage format, and thus an
> > implementation detail.
> Sorry for my poor presentation. I simply meant that one has to
> make a design decision on whether the i in A(i,j), for instance,
> means the i-th row, the i-th column, or both (depending on the
> storage format). This is not just an implementation detail because
> it concerns whether expressions like
>    (A*u)(1) = A(1,1)*u(1) + ... + A(1,n)*u(n)
> are promised to be identities.

This example is not quite the same thing and is mixing two different 
levels of concepts.  To the extent that a matrix A is modeled as a 
tabular arrangement of numbers A(i,j) is the element at row i and 
column j.  This is the mathematical convention.  And, mathematically, 
the expression you wrote would hold (modulo the normal disclaimers 
about floating point arithmetic).  That is, mathematically, the 
relationship  v(i) = sum(A(i,j)*u(j)) would hold when multiplying A*u.

However, as shown above, you would never want to write the matrix-
vector product as two nested loops over integers i and j.  In between 
our mathematical abstraction of a matrix -- with indexing convention 
as above -- is a layer mapping from underlying storage (possibly an 
array, possibly a sparse format of some kind, possibly something 
else), to the mathematical.  

The issue of row major versus column major comes up in this layer 
when you have to decide how you actually want to stuff your elements 
into storage.  In the dense case this is an issue of mapping between 
the A(i,j) abstraction of the matrix to some A_rep[i][j] or A_rep[j]
[i] or A_rep[i*lda+j] or A_rep[j*lda+i] array representation.

> 
> > > I hope the Boost folks can draw a line between vectors, matrices
> > > and arrays in their matrix & array library.
> > 
> > What sort of line did you have in mind?
> Uhh, I don't have one in mind, because I don't know how other people
> use arrays. At any rate, it is essential to identify different
> orthogonal concepts, and not to let the base array class (if any)
> assume too much. If I just want to interface my old legacy array
> architecture with the library functions (exchange two rows, slice,
> ... etc.), I don't want to define array element multiplication
> just because the library assumes that my array is a matrix.

Based on previous exchanges across this list, I think it is becoming 
fairly certain that there will be a very careful distinction drawn 
between matrix and array.

Regards,
Andrew Lumsdaine