ANN: hmatrix-static: statically-sized linear algebra

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. Installation and requirements: Because of the use of type families and view patterns, hmatrix-static requires GHC 6.10 to build. The linear algebra routines are ultimately implemented by BLAS, GSL and LAPACK; see the hmatrix manual [3] for instructions on installing these, prior to installing hmatrix. Once this has been done, hmatrix-static can be installed by $ cabal install hmatrix-static I would be very interested in your feedback, so please try it out. Kind regards, Reiner Pope [1] http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix-static [2] http://www.hmatrix.googlepages.com/ [3] http://www.hmatrix.googlepages.com/installation

Reiner Pope
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/* <\ * (__ */\ <

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
Reiner Pope
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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

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
wrote: Reiner Pope
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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
participants (3)
-
Alberto Ruiz
-
Reiner Pope
-
Xiao-Yong Jin