[Haskell-cafe] matrix computations based on the GSL

[David Roundy]

On Wed, Jun 29, 2005 at 10:23:23PM +0200, Henning Thielemann wrote:
> More specific:
>  You give two different things the same name. You write
>   A*B
>    and you mean a matrix multiplication. Matrix multiplication means
> finding a representation of the composition of the operators represented
> by A and B.
>  But you also write
>   A*x
>  and you mean the matrix-vector multiplication. This corresponds to the
> application of the operator represented by A to the vector x.
>  You see: Two different things, but the same sign (*). Why? You like this
> ambiguity because of its conciseness. You are used to it. What else?
>  But you have to introduce an orientation of vectors, thus you
> discriminate two things (row and column vectors) which are essentially the
> same!
>  What is the overall benefit?

I just wanted to insert a few thoughts on the subject of row and column
vectors.  (And I'm not going to address most of the content of this
discussion.)

If we support matrix-matrix multiplication, we already automatically
support matrix-column-vector and row-vector-matrix multiplication, whether
or not we actually intend to, unless you want to forbid the use of 1xn or
nx1 matrices.  So (provided we *do* want to support matrix-matrix
multiplication, and *I* certainly would like that) there is no question
whether we'll have octavish/matlabish functionality in terms of row and
column vectors--we already have this behavior with a single multiplication
representing a single operation.  There is indeed a special case where one
or more of the dimensions is 1, but it's just that, a special case, not a
separate function.

If you want to introduce a more general set of tensor datatypes, you could
do that, but matrix arithmetic is all *I* really need.  I agree that when
you add other data types you'll need other operators (or will need to do
something more tricky), but starting with the simple case seems like a good
idea, especially since the simple case is the one for which there exist
fast algorithms (i.e. level 3 BLAS).

There are good reasons when if I've got a set of vectors that I'd like to
multiply by a matrix that I should group those vectors into a separate
matrix--it's far faster (I'm not sure how much, but usually around 10 times
faster, sometimes more).  In fact, I'd think it's more likely that a
general tensor library will be implemented using matrix operations than the
other way around, since it's the matrix-matrix operations that have
optimized implementations.

In short, especially since the folks doing the work (not me) seem to want
plain old octave-style matrix operations, it makes sense to actually do
that.  *Then* someone can implement an ultra-uber-tensor library on top of
that, if they like.  And I would be interested in a nice tensor
library... it's just that matrices need to be the starting point, since
they're where most of the interesting algorithms are (that is, the ones
that are interesting to me, such as diagonalization, svd, matrix
multiplication, etc).
-- 
David Roundy
http://www.darcs.net