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

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