[Haskell-cafe] matrix computations based on the GSL

[Henning Thielemann]

On Wed, 29 Jun 2005, Conal Elliott wrote:

> On row & column vectors, do you really want to think of them as
> {1,...,n)->R?  They often represent linear maps from R^n to R or R to
> R^n, which are very different types.  Similarly, instead of working with
> matrices, how about linear maps from R^n to R^m?  In this view, column
> and row vectors, matrices, and often scalars are useful as
> representations of linear maps.

We should not identify things which can be mapped bijectively. "1" and 1
are very different, although they may mean the same in a given context.
Yes it is possible to describe linear maps with vectors but vectors are
not linear maps. Namely, x |-> <x,y> is a linear map, so is y itself a
linear map? Certainly not!
 If you want to put an interpretation of row or column into a vector then
do it when you actually do the linear map but keep the vector itself free
of this information.

> I've played around some with this idea of working with linear maps
> instead of the common representations, especially in the context of
> derivatives (including higher-dimensional and higher-order), where it is
> the view of calculus on manifolds.  It's a lovely, unifying approach and
> combines all of the various chain rules into a single one.  I'd love to
> explore it more thoroughly with one or more collaborators.

I think matrices and derivatives are very different issues. I have often
seen that the first derivative is considered as vector, and the second
derivative is considered as matrix. In this spirit it is used like
  x^T * (D2 f)(x) * x
 but this is only abuse of the common multiplication definitions. A good
interpretation and notation should seamless extend to higher derivatives.
But the interpretation above does not work in higher dimensions.

I like the following type for derivation.
 derive :: ((i -> a) -> b) -> ((i -> a) -> (i -> b))
  Here i is the index type, (i -> a) is the vector type, b is the type the
vector function maps to. Its derivative has the same type of argument, but
the result is a vector with indices of type i. You see that it is easy to
repeat the application of 'derive', just replace b by say i->b. The second
derivative yields vectors of type (i -> i -> b). This can be interpreted
as matrix because it has two indices. But this is certainly not a matrix
which represents a linear mapping as usual, but it is a matrix
representing a bilinear form. The only thing we need is a multiplication
to reduce one level of indices.
 mul :: (i -> c) -> (i -> b) -> b
  Though, what we still need is a general (overloaded?) definition of the
scaling of b by c and a sum of b.