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/