[Haskell] ANNOUNCE: An index-aware linear algebra library in Haskell

[Frederik Eaton]

An index-aware linear algebra library in Haskell

I've been exploring the implementation of a library for linear
algebra, i.e. manipulating vectors and matrices and so forth, which
has as a fundamental design goal the exposure of index types and
ranges to the type system so that operand conformability can be
statically guaranteed (e.g. an attempt to add two matrices of
incompatible sizes, or multiply an NxK matrix by an MxL matrix where K
/= M, is flagged statically by the compiler).

A significant contribution of this work is a mechanism for supporting
interactive use. This is challenging because it requires instantiating
values with user-specified index types directly, rather than via CPS
rank-2 polymorphic "with" functions as in Kiselyov and Shan's
"Implicit Configurations". The solution is based on Template Haskell.


MOTIVATION:

Haskell has never been a serious language choice for numerics, partly
because of efficiency problems. However, a large portion of the demand
for number-processing in science is met by other interpreted languages
such as Octave and Matlab. These languages are just as slow, or
slower, than Haskell when it comes to executing a large number of
operations in succession. Their advantage is in making available to
the user a collection of highly-optimized linear algebra routines. 
Since many numerical tasks consist of simple, repetitive operations on
large amounts of data, and because these operations can often be
reformulated in terms of operations on vectors and matrices, such
matrix languages can give good performance on appropriately written
code, and have become standard, viable workbenches for numerical
analysis.

However, matrix languages are often good for little more than matrix
manipulations - their support for more complicated data structures is
poor. Furthermore, they are generally dynamically typed, which means
that type errors only become evident at runtime. In particular,
operand conformability errors, caused for instance by multiplying
matrices A*B instead of B*A or A'*B' (in Matlab notation), are runtime
errors, and such errors can be masked when dimensions happen to be
equal, for instance when A and B are square matrices as in the above
example. I believe that with a good linear algebra library which
solves this problem via static typing, Haskell could make a valuable
contribution in the area of technical computing.


DESIGN:

- The fundamental type is a "vector", which includes an element type
and an index type. Matrices are vectors indexed by pairs. 

Vectors are members of the Vector class. The element type is a class
parameter, to allow vectors over specialized types:

class Vector v e | v -> e

(note, the index type is not a class parameter)

  http://ofb.net/~frederik/futility/src/Vector/Base.hs

- The building-block index type wraps a type-level integer N, and
accepts a range of possible values 0..(N-1). I use a modified version
of Oleg Kiselyov and Chung-chieh Shan's "Implicit Configurations"
paper for the type-level integers.

Index types must implement the Dom ("domain") class. Superclasses are
Bounded, Enum, Ix, Eq, and Show. Instances are defined so that index
types may be combined with (,) and Either.

Additionally, there is a generics-style method domCastPair which
allows us to detect when a vector is also a matrix. The primary intent
of this is to make it possible for libraries such as Alberto Ruiz's
GSL library to implement our interface, using different data types for
matrices and vectors "under the hood". I also use domCastPair to
display vectors and matrices differently.

  http://ofb.net/~frederik/futility/src/Domain.hs
  http://ofb.net/~frederik/futility/src/Prepose.hs

- I have provided an example implementation of most of the operations
using the Array type. This is very slow for numerical operations
(about 200x slower than Octave), but it allows arbitrary (boxed)
element types, which faster implementations are unlikely to do. 

  http://ofb.net/~frederik/futility/src/Vector/Array.hs

Thus in the current state the library is a working prototype. For it
to be practically useful, it still needs a Vector instance which wraps
some fast linear algebra library such as the GSL (at the expense of
only supporting a restricted set of element types). A faster
Haskell-only implementation would also be possible using unboxed
arrays, but I think it is unlikely that this would be able to come
close to the performance of an external C or Fortran library such as
GSL.

- Due to the need to specify index types at some point, input of
vectors is difficult. I have provided two functions in Fu.Vector.Base
which should cover most of the cases:

listVec :: Vector v e => [e] -> (forall n . (ReflectNum n) => v (L n) -> w) -> w 
listMat :: Vector v e => [[e]] -> 
    (forall n m . (ReflectNum n, ReflectNum m) => v (L m, L n) -> w) -> w 

However, these aren't useful in interactive situations. So I have also
provided some template-haskell routines

  http://ofb.net/~frederik/futility/src/Vector/Template.hs

which can be used to instantiate vectors directly. For example:

 -- 'ident' is an identity matrix with polymorphic dimensions
 -- $(cm 3 3) specifies a 3x3 matrix
 -- aV specifies use of the Array implementation
> let x = aV $ $(cm 3 3) ident
> x
<# 1, 0, 0; 0, 1, 0; 0, 0, 1 #>
> let v = x + ones
> v
<# 2, 1, 1; 1, 2, 1; 1, 1, 2 #>

See also

  http://ofb.net/~frederik/futility/src/Vector/examples.hs

I haven't finalized the names in Vector/Template.hs. While it seems
desirable that they be short ('cm', 'aV'), I imagine the names I have
chosen, and the interface in general, could be improved upon.


PROBLEMS stemming from shortcomings in the language / GHC:

- Any code using template haskell cannot be profiled.

- When vector index types are known explicitly, for instance with
interactive use, errors can be difficult to parse:

> ($(cv 5) ones) + ($(cv 4) ones)

<interactive>:1:18:
    Couldn't match `Succ (Twice (Twice (Succ Zero)))'
           against `Twice (Twice (Succ Zero))'
      Expected type: v (L (Succ (Twice (Twice (Succ Zero)))))
      Inferred type: v (L (Twice (Twice (Succ Zero))))
    In the application `$[splice](cv 4) ones'
    In the second argument of `(+)', namely `($[splice](cv 4) ones)'

A base-10 implementation of the type-level integers might help, but
even better would be language support for dependent types.

- Class annotations are repetitive and tiresome. The code starts to
look like Prolog.

  vzipWith_ :: (Dom a, Vector v k, Vector u l, Vector w m) => (k -> l -> m) -> v a -> u a -> w a 

A "functional" syntax would be cleaner:

  vzipWith_ :: (k -> l -> m) -> V k a -> V l a -> V m a

where "v = V k a" implies a constraint (Vector v k, Dom a). Oleg
suggests this might be accomplished with better support for type
aliases:

  type V k a = (Vector v k, Dom a) => v a

"Partial signatures", i.e. being able to specify types without class
constraints, would also be nice.


PROBLEMS stemming from shortcomings in the approach:

- In Octave, it is easy to create a new vector from a subset of the
indices of an existing vector, by specifying a boolean relation on
those indices:

> a=[1,0,1,1];
> b=[1,2,3,4];
> b(a==1)
ans =
  1  3  4

Doing this sort of thing in an environment where vector indices are
strongly-typed seems difficult. Probably the best option in this
particular case would be returning a list of values.

- Lazy evaluation is likely to be a source of performance headaches,
although I haven't actually investigated this.

- Many expressions can be better optimized if one has access to the
structure of the whole expression. For instance, given matrices A, B,
and C, a specialized routine might calculate A*B+C in one step, and it
may be that A*(B*C) is more efficient to calculate than (A*B)*C. A
special-purpose language such as Matlab can perform these kinds of
optimizations easily. It is more difficult to do them in a library,
although it is not impossible.


FUTURE work:

The set of standard Vector operations needs to be fleshed out. A
wrapper for Alberto's GSL library, or another fast matrix library,
should be written.


ACKNOWLEDGMENTS:

Thanks to Oleg Kiselyov for many useful discussions, and to Ralf
Lammel for help understanding generics.

-- 
http://ofb.net/~frederik/