Re: [Haskell-cafe] matrix computations based on the GSL

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 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