Subject: ANNOUNCE: An index-aware linear algebra library in Haskell Newsgroups: gmane.comp.lang.haskell.general Date: Friday 14th April 2006 15:02:01 UTC (over 12 years ago) 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) |
|||