[Haskell-cafe] matrix computations based on the GSL

[David Roundy]

On Thu, Jul 07, 2005 at 03:08:50PM +0200, Henning Thielemann wrote:
> 
> On Thu, 7 Jul 2005, David Roundy wrote:
> 
> > On Tue, Jul 05, 2005 at 08:17:58PM +0200, Henning Thielemann wrote:
> > > The example, again: If you write some common expression like
> > >
> > > transpose x * a * x
> > >
> > > then both the human reader and the compiler don't know whether x is a
> > > "true" matrix or if it simulates a column or a row vector. It may be that
> > > 'x' is a row vector and 'a' a 1x1 matrix, then the expression denotes a
> > > square matrix of the size of the vector simulated by 'x'. It may be that
> > > 'x' is a column vector and 'a' a square matrix. Certainly in most cases I
> > > want the latter one and I want to have a scalar as a result. But if
> > > everything is a matrix then I have to check at run-time if the result is a
> > > 1x1 matrix and then I have to extract the only element from this matrix.
> > > If I omit the 1x1 test mistakes can remain undiscovered. I blame the
> > > common notation x^T * A * x for this trouble since the alternative
> > > notation <x, A*x> can be immediately transcribed into computer language
> > > code while retaining strong distinction between a vector and a matrix
> > > type.
> >
> > The issue is that Haskell (as far as I understand, and noone has suggested
> > anything to the contrary) doesn't have a sufficiently powerful type system
> > to represent matrices or vectors in a statically typed way.  It would be
> > wonderful if we could represent matrix multiplication as
> >
> > matrix_mul :: Matrix n m -> Matrix m o -> Matrix n o
> 
> Even if we had that I would vote for distinction of Vector and Matrix! I
> wanted to show with my example that vectors and matrices are different
> enough that it makes sense to give them different types. In my opinion
> multiplying a matrix with a so-called column vector is a more fundamental
> bug than multiplying a matrix with a vector of non-compatible size. That
> is I like static check of matrix-vector combinations and a dynamic check
> of their size. The latter also because in many application you don't know
> the vector size at compile time.

You wouldn't need to know the vector size at compile time in order to get
static type checking of vector sizes--if the type system is sufficiently
powerful.

> > Also note that if you have several vectors x for which you want to compute
> > the dot product with metric A, and if you want to do this efficiently,
> > you'll have to convert your list of vectors into a matrix anyways.
> 
> If you bundle some vectors as columns in matrix B and want to compute
> norms with respect to matrix A writing
>  B^T * A * B
>   you will not only get the norms of the vectors in B but also many mixed
> scalar products. This example shows to me that matrices are not simply
> collections of vectors.

Good point, and that does illustrate the pain of doing this (treating sets
of vectors as matrices).

> Btw. we should not mess up the design with decisions on efficiency
> concerns. If you want efficiency then you can abuse matrices for that but
> I consider a 'map' over a list of vectors as the cleaner way (and more
> type safe, because you can't make transposition errors).

Mapping is cleaner, but something like 10 times slower... not for my
example above, but for a simple y = A*x computation--you then could do the
vector-vector inner product better than transpose y * x.

Efficiency concerns shouldn't be separated from design decisions.  An API
shouldn't force an implementation to be inefficient.  On the other hand,
this is sort of a silly debate, since the API *I* want is a subset of the
API *you* want.  (Except that I'd like matrices to be an instance of
Floating, but that's easy enough to do once we've got all the matrix
operations...)
-- 
David Roundy