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

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

It's less type-safe and thus cumbersome to work with. That's my experience with MatLab.

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