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

[Alberto Ruiz]

Hi Reiner,

Fantastic work! User-friendly static dimension checking is an essential 
feature for any decent linear algebra library. Your interface using 
quasiquotation and view patterns is very elegant and practical. I am 
happy that hmatrix is useful, but I'm afraid that its primitive dynamic 
interface will no longer be used :)

Many thanks for your excellent work!

Alberto

Reiner Pope wrote:
> 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
>>
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>