[Haskell-cafe] matrix computations based on the GSL

[David Roundy]

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

(where n, m and o are dimensions) but we can't do that.  So we're stuck
with dynamic size checking for all matrix and vector operations that
involve two or more operands.  This is just something we have to live with
for now.  If anyone knows differently, please speak up! Perhaps one could
do this with template haskell, but that sounds scary and a bit ugly
(although I'd love to be proven wrong).

You'd like to have certain special cases where a dimension 1 is treated
differently from other size dimensions.  I'd say that it'd be simpler to
start out with everything done in a general manner, and then you can always
add special cases later if you'd like.  Given that we're stuck with a
dynamically typed system, I'd prefer to have a simple dynamically typed
system.

You think of vectors as fundamental and matrices as a sort of derived
quantity.  I'd prefer to think the other way around, with vectors being
special cases of matrices, since it's simpler to implement and work with.
Since matrices form a complete self-contained algebra (and I'm sure I'm
massacring the terminology--sorry), it seems simplest to implement that
first, especially since you'd want to implement it anyways.

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.  Writing
functions of vectors instead as functions of 1xN matrices allows them to be
efficiently applied to multiple vectors simultaneously, provided they are
written carefully.  Alas, if we had static dimension checking, we wouldn't
have to worry about writing carefully, but alas that doesn't seem to be an
option.
-- 
David Roundy
http://www.darcs.net