[Haskell-cafe] ANN: hmatrix-static: statically-sized linear algebra

[Reiner Pope]

There should be essentially no performance penalty in going from
hmatrix to hmatrix-static, for most functions. The Vector and Matrix
types I define are just newtypes of hmatrix's Vector and Matrix types,
so they will have the same runtime representation.

There are a few cases where performance could potentially differ,
described below.

Reflecting and reifying Ints (i.e. converting between types and values):
(The ideas for this come from the Implicit Configurations paper[1].)
Some Int arguments in hmatrix have been written as types in
hmatrix-static, such as the function "constant".

hmatrix type:
constant :: Element a => a -> Int -> Vector a
hmatrix-static type:
constant :: (Element a, PositiveT n) => a -> Vector n a

The PositiveT constraint is essentially a way of passing the Int
parameter as an implicit parameter. To demonstrate this, we use two
library functions which allow us to convert from Ints to types with
PositiveT constraints, and back:

value :: PositiveT n => Proxy n -> Int    -- type -> value
reifyPositive :: Int -> (forall n. PositiveT n => n -> r) -> r   --
value -> type.

The use of these functions is nice in some cases, such as "constant"
above, because it allows us to pass parameters implicitly. On the
other hand, the conversion between types and values imposes an O(log
n) cost, where n is the size of the number we are converting. In my
experience, this cost has not been significant (although it was
previously, when I used a naive O(n) reifyIntegral implementation!).

Newtype conversion costs:
Occasionally, unwrapping newtypes can actually have a runtime cost.
For instance, the hmatrix-static function

joinU :: Storable t => [Vector Unknown t] -> Vector Unknown t

needs to do a conversion :: [Vector Unknown t] -> [HMatrix.Vector t].
Now, the conversion "unwrap :: Vector Unknown t -> HMatrix.Vector t"
is a noop, as it unwraps a newtype. However, to get the list
conversion, we need to do "map unwrap", which requires mapping a noop
over a list. However, because of map's recursion, GHC may not be able
to recognise that "map unwrap" is a noop, and traverse the list
anyway, causing a performance loss.

However, there aren't many recursive data structures used in
hmatrix-static, so this problem mostly doesn't exist.

Cheers,
Reiner

[1] http://okmij.org/ftp/Haskell/types.html#Prepose

On Sun, Apr 12, 2009 at 12:18 PM, Xiao-Yong Jin <xj2106 at columbia.edu> wrote:
> Reiner Pope <reiner.pope at gmail.com> writes:
>
>> Hi everyone,
>>
>> I am pleased to announce hmatrix-static[1], a thin wrapper over
>> Alberto Ruiz's excellent hmatrix library[2] for linear algebra.
>>
>> The main additions of hmatrix-static over hmatrix are:
>>  - vectors and matrices have their length encoded in their types
>>  - vectors and matrices may be constructed and destructed using view
>> patterns, affording a clean, safe syntax.
>>
>> hmatrix-static includes statically-sized versions of hmatrix's linear
>> algebra routines, including:
>>  - simple arithmetic (addition, multiplication, of matrices and vectors)
>>  - matrix inversion and least squares solutions
>>  - determinants / rank / condition number
>>  - computing eigensystems
>>  - factorisations: svd, qr, cholesky, hessenberg, schur, lu
>>  - exponents and sqrts of matrices
>>  - norms
>>
>> See http://code.haskell.org/hmatrix-static/examples/ for example code.
>>
>
> This is quite interesting.  I would like to know the
> performance penalty of such a wrapper.  I'll test it by
> myself when I've got time.  But can you give me some idea of
> how such a wrapper can be optimized by ghc in terms of space
> and time performance?
> --
>    c/*    __o/*
>    <\     * (__
>    */\      <
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>