Are new sequences really O(1)?

Adrian Hey

24 May 2005 24 May '05

3:29 p.m.

Hello, Just been trying a few simple benchmarks to compare the new sequences with AVL trees for simple deque operations and I'm getting some strange results. Code for both is attached at the end of this post, but basically the test in each case is: 1- Build a sequence of sz elements by pushing them one by one on the right. 2- Rotate the sequence right sz times, each rotation pops the rightmost element and pushes it on the left. 3- Empty the sequence by popping the leftmost element sz times. With ghc-6.4, results are (time in no particular units) sz AVL Sequence ---------------------- 2^10 0.278 0.177 2^11 0.666 0.425 2^12 1.464 0.958 2^13 3.112 2.214 2^14 6.664 5.525 2^15 14.18 12.18 2^16 30.16 27.31 2^17 62.14 60.6 2^18 130.4 142.6 2^19 269.3 365.1 2^20 564.4 ***** The AVL figures seem to show roughly the kind of n*(log n) growth I would expect. The Sequence figures start out promisingly enough, but seem to get progressively worse until it's actually slower than AVL. This isn't what I'd expect from an algorithm advertised as having O(1) asymptotic complexity for pushing/popping. It seems more like O(log n) or worse? Also there are no figures for 2^20 for Sequence because I get a stack overflow at this point. So any idea whether this is a bug in my understanding?, or a bug in the theory?, or a bug in the code perhaps? Maybe it's being excessively lazy somewhere (often the cause of stack overflows IME). Regards -- Adrian Hey Code follows: {-# OPTIONS -fno-cse -fno-full-laziness #-} module Main (main) where import Data.Sequence import System.CPUTime(getCPUTime,cpuTimePrecision) import System.Mem(performGC) result :: Int -> Seq () result sz = rep pop $ rep rot $ rep push empty where rep = rep' sz rep' 0 f x = x rep' n f x = let x' = f x in x' `seq` rep' (n-1) f x' push sq = sq |> () rot sq = case viewR sq of sq' :> x -> x <| sq' EmptyR -> empty pop sq = case viewL sq of _ :< sq' -> sq' EmptyL -> undefined test :: (Int,Int) -> IO () test (n,sz) = do performGC t0 <- getCPUTime rep n t1 <- getCPUTime putStr $ show sz ++ " : " putStrLn $ show $ (fromIntegral ((t1-t0) `div` cpuTimePrecision)) / (fromIntegral n) where rep 0 = return () rep m = (result sz) `seq` rep (m-1) main :: IO () main = mapM_ test [(10*2^(maxP-p), 2^p) | p <- [10..maxP]] where maxP = 20 ---------------------------------------------------------------- {-# OPTIONS -fno-cse -fno-full-laziness #-} module Main (main) where import Data.Tree.AVL import System.CPUTime(getCPUTime,cpuTimePrecision) import System.Mem(performGC) result :: Int -> AVL () result sz = rep pop $ rep rot $ rep push empty where rep = rep' sz rep' 0 f x = x rep' n f x = let x' = f x in x' `seq` rep' (n-1) f x' push sq = pushR sq () rot sq = case popR sq of (sq', x ) -> pushL x sq' pop sq = case popL sq of (_ , sq') -> sq' test :: (Int,Int) -> IO () test (n,sz) = do performGC t0 <- getCPUTime rep n t1 <- getCPUTime putStr $ show sz ++ " : " putStrLn $ show $ (fromIntegral ((t1-t0) `div` cpuTimePrecision)) / (fromIntegral n) where rep 0 = return () rep m = (result sz) `seq` rep (m-1) main :: IO () main = mapM_ test [(10*2^(maxP-p), 2^p) | p <- [10..maxP]] where maxP = 20

Show replies by date

Ross Paterson

26 May 26 May

9:35 a.m.

On Tue, May 24, 2005 at 08:29:57PM +0100, Adrian Hey wrote:

...

Just been trying a few simple benchmarks to compare the new sequences with AVL trees for simple deque operations and I'm getting some strange results.

The stack overflow is interesting, and relates to all representations based on lazy number systems. If you build a structure without accessing it, you get n/3 nested applications at the second level. That's no more space than they will evaluate to, but as soon as you access the second level, they all go on the stack. GCs that happen during this process will be more expensive, as they have to scan the stack. I suspect that GC costs are swamping everything else for large n. I've fixed the stack overflow; I believe this doesn't break the persistent bounds.

Adrian Hey

27 May 27 May

10:30 a.m.

On Thursday 26 May 2005 2:35 pm, Ross Paterson wrote:

...

On Tue, May 24, 2005 at 08:29:57PM +0100, Adrian Hey wrote:

...
Just been trying a few simple benchmarks to compare the new sequences with AVL trees for simple deque operations and I'm getting some strange results.

The stack overflow is interesting, and relates to all representations based on lazy number systems. If you build a structure without accessing it, you get n/3 nested applications at the second level. That's no more space than they will evaluate to, but as soon as you access the second level, they all go on the stack. GCs that happen during this process will be more expensive, as they have to scan the stack. I suspect that GC costs are swamping everything else for large n.

I've fixed the stack overflow; I believe this doesn't break the persistent bounds.

This seems to have improved things quite a bit. Hope your right about the persistent bounds (I wouldn't know :-). But it would be good to test this somehow. Anyway, latest test results from 2^8 to 2^22 are.. sz AVL Sequence 2^ 8 0.0448 0.0331 2^ 9 0.1015 0.0676 2^10 0.2789 0.1490 2^11 0.6707 0.3479 2^12 1.473 0.8010 2^13 3.142 1.676 2^14 6.733 3.461 2^15 14.39 7.181 2^16 30.41 15.36 2^17 62.73 30.82 2^18 130.8 62.89 2^19 270.9 126.6 2^20 559.2 254.8 2^21 1155 513.3 2^22 2411 1035 This is quite impressive. For moderate length sequences, it seems to be about twice as fast as AVL (which AFAIK is the fastest balanced tree implementation currently available for Haskell). The execution times are growing slower than AVL for longer sequences too (looks more like O(1) than it did before, but still not quite O(1) for some reason). Regards -- Adrian Hey

Einar Karttunen

11:20 a.m.

Adrian Hey writes:

...

This is quite impressive. For moderate length sequences, it seems to be about twice as fast as AVL (which AFAIK is the fastest balanced tree implementation currently available for Haskell). The execution times are growing slower than AVL for longer sequences too (looks more like O(1) than it did before, but still not quite O(1) for some reason).

Garbage collection is not O(1) and data moves out of cache, so there will be no true O(1), just something which looks like it would be O(1) in a perfect world. - Einar Karttunen

Keean Schupke

9 Jun 9 Jun

7:07 a.m.

New subject: Haskell Matrix Library...

Hi, I have a simple matrix library (styled after MatLab's use of matrices), and was wondering if any of the standard libraries include matrix functionality. The library as stands is not suitable for general use - but it could be adapted. I guess the questions are: - is there already a matrix library in the standard libraries? - is anyone interested in a matrix library? - would people want a library like MatLab where matrices are instances of Num/Fractional/Floating and can be used in normal math equations... - operations like cos/sin/log can be applied to a matrix (and apply to each element like in MatLab)... What do people think? Keean.

David Roundy

7:29 a.m.

New subject: Haskell Matrix Library...

[Turning questions around...] On Thu, Jun 09, 2005 at 12:07:54PM +0100, Keean Schupke wrote:

...

What do people think?

I think this would be cool. One of the things that saddens me about haskell is that I can't really use it for work (or at least, I don't think I could). One thing you haven't mentioned is the ability to create data structures containing matrices and use them nicely (i.e. with operator overloading, etc). This is something that I'd like, so I could create a data type WaveFunction that is an array, but also contains its quantum numbers and other relevant information. In matlab you're stuck with things being a pure matrix, which is highly inconvenient, and that's also the case for most

...

- is anyone interested in a matrix library?

Yes.

...

- would people want a library like MatLab where matrices are instances of Num/Fractional/Floating and can be used in normal math equations...

Absolutely. The only annoyance would be that due to haskell's strongly typed approach and lack of automatic type conversions, *everything* would have to be a matrix. That's no worse than matlab, but it's not one of the features of matlib I like. But at least haskell's way of handling rational numbers lends itself to this, since when you make a matrix an instance of Rational (which means defining fromRational), floating point constants can be automatically treated as 1x1 matrices.

...

- operations like cos/sin/log can be applied to a matrix (and apply to each element like in MatLab)...

Yeah, this does seem necesary. Especially since I envision that *everything* is likely to be a matrix (including scalars). If you could implement even a moderate fraction of octave's functionality (octave being a free matlab clone, for any listeners unfamiliar with it), it would be great. Matlab's got a horrible language, but it's just so darn convenient for dealing with linear algebra. I'd love to be able to recommend that new students learn haskell rather than octave (which I despise). I'm guessing that your simple matrix class doesn't call lapack or do blocked matrix multiplies? A blocked matrix multiply of the fftw/atlas variety would be an interesting project to code in template haskell... -- David Roundy http://www.darcs.net

Henning Thielemann

7:54 a.m.

New subject: Haskell Matrix Library...

On Thu, 9 Jun 2005, David Roundy wrote:

...

In matlab you're stuck with things being a pure matrix, which is highly inconvenient, and that's also the case for most

MatLab has also 'cell arrays' which are dynamically typed records. That's not a good solution because you never know what's actually in a record but at least the statement "you are stuck to matrices" is not true. :-)

Henning Thielemann

7:35 a.m.

New subject: Haskell Matrix Library...

On Thu, 9 Jun 2005, Keean Schupke wrote:

...

Hi, I have a simple matrix library (styled after MatLab's use of matrices), and was wondering if any of the standard libraries include matrix functionality.

Not more than arrays, as far as I know.

...

The library as stands is not suitable for general use - but it could be adapted. I guess the questions are:

- is there already a matrix library in the standard libraries?

There is a small linear algebra library of Jan Skibinski and a linear algebra part of the HaskellDSP library. I remember that this question already arose and resulted in some more interesting links.

...

- is anyone interested in a matrix library?

Yes.

...

- would people want a library like MatLab where matrices are instances of Num/Fractional/Floating and can be used in normal math equations...

I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components. What is the natural definition of matrix multiplication? The element-wise multiplication (which is commutative like the multiplication of other Num instances) or the matrix multiplication? I vote for separate functions or infix operators for matrix multiplications, linear equation system solvers, matrix-vector-multiplication, matrix and vector scaling. I wished to have an interface to LAPACK for the floating point types because you can hardly compete with algorithms you write in one afternoon. Btw. when considering the class hierarchy of the numeric prelude project http://cvs.haskell.org/darcs/numericprelude/ you would naturally put matrices and vectors in Additive, Module and VectorSpace class but not in Fractional.

...

- operations like cos/sin/log can be applied to a matrix (and apply to each element like in MatLab) ...

MatLab must provide these automatisms because it doesn't have proper higher functions. I think the Haskell way is to provide a 'map' for matrices. Otherwise the question is: What is the most natural 'exp' on matrices? The elementwise application of 'exp' to each element or the matrix exponentation?

Keean Schupke

8:08 a.m.

New subject: Haskell Matrix Library...

Henning Thielemann wrote:

...

...
- would people want a library like MatLab where matrices are instances of Num/Fractional/Floating and can be used in normal math equations...

I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components. What is the natural definition of matrix multiplication? The element-wise multiplication (which is commutative like the multiplication of other Num instances) or the matrix multiplication? I vote for separate functions or infix operators for matrix multiplications, linear equation system solvers, matrix-vector-multiplication, matrix and vector scaling. I wished to have an interface to LAPACK for the floating point types because you can hardly compete with algorithms you write in one afternoon. Btw. when considering the class hierarchy of the numeric prelude project http://cvs.haskell.org/darcs/numericprelude/ you would naturally put matrices and vectors in Additive, Module and VectorSpace class but not in Fractional.

I think we can have the best of both worlds...

...

...
- operations like cos/sin/log can be applied to a matrix (and apply to each element like in MatLab) ...

MatLab must provide these automatisms because it doesn't have proper higher functions. I think the Haskell way is to provide a 'map' for matrices. Otherwise the question is: What is the most natural 'exp' on matrices? The elementwise application of 'exp' to each element or the matrix exponentation?

MatLab has "exp" which works element wise and "expm" which does the matrix-exponentation. I have made the matrix an instance of Num/Fractional and Float. The float operations all apply element wise. I have constructors "scalar" and "vector" to produce a 1x1 and a row-matrix respectively. The '*' operator does matrix multiplication (as in MatLab) and I have a new class for matrices which defines ".*" for element wise multiplication. Like matlab, if you try to multiply a scalar (1x1 matrix) by another matrix you get element wise multiplication instead of matrix multiplication... (as a scalar can only really be multiplied by another scalar this is unabiguous). What about implementation? At the moment I am coding the vector operations in Haskell (a foreign interface to lapack or Intel MKL can come later)... and am using a UArray - is this portable enough - or should I parameterise by array type (using IArray) even though this might result in type ambiguities? Keean.

Henning Thielemann

8:16 a.m.

New subject: Haskell Matrix Library...

On Thu, 9 Jun 2005, Keean Schupke wrote:

...

Henning Thielemann wrote:

...
MatLab must provide these automatisms because it doesn't have proper higher functions. I think the Haskell way is to provide a 'map' for matrices. Otherwise the question is: What is the most natural 'exp' on matrices? The elementwise application of 'exp' to each element or the matrix exponentation?

MatLab has "exp" which works element wise and "expm" which does the matrix-exponentation.

This would be the worsed possible mix - just in the spirit of MatLab. :-) With (*) meaning a matrix multiplication and 'exp' meaning elementwise exponentation, the exponential law exp (x+y) = exp x * exp y does not hold. For elementwise (*) operation it would hold and for 'exp' as matrix exponentation it would hold at least for matrices x and y which are diagonalisable by the same similarity transform.

David Roundy

8:23 a.m.

New subject: Haskell Matrix Library...

On Thu, Jun 09, 2005 at 01:08:48PM +0100, Keean Schupke wrote:

...

What about implementation?

At the moment I am coding the vector operations in Haskell (a foreign interface to lapack or Intel MKL can come later)... and am using a UArray - is this portable enough - or should I parameterise by array type (using IArray) even though this might result in type ambiguities?

Personally, I'd lean towards a ForeignPtr for the back end from the beginning--but that's partly just because I'm very comfortable with ForeignPtrs. This is partly because I'm not sure how to pass IArray through the FFI, and if you can't figure out how to do that you'll have to change your storage mechanism (or do a copy... yuck) when you want to interface with lapack. Also, the FFI is pretty portable, and there's obviously going to be no boxing going on. It means using unsafePerformIO all over the place, but you can pretty easily encapsulate that in a completely safe manner. Another interesting question is what to do with complex numbers. Matlab just make everything complex, which is pretty easy, but sometimes leads to uglinesses (real matrices with tiny imaginary parts due to roundoff). But this question is one that can probably be deferred. On another interesting note, I wonder how hard it would be to write a compile-time optimized block matrix multiply in the spirit of ATLAS using template haskell? One could write a little program to generate blocking algorithms and test them all, and then generate the best one... :) I'm sure it'd be a lot easier to write in haskell than in C, and wonder how close we could come in performance. -- David Roundy http://www.darcs.net

Keean Schupke

12:52 p.m.

New subject: Haskell Matrix Library...

Please CC me in to replies, as I don't seem to get all the posts sent to the list for some reason... I will post a basic version of the Matrix library in the next couple of days for comments, bear in mind that my goals are: - keep it as small as possible. - only implement what I need now, but make sure it is extensible to cope with other functions that could be added in future. - make is as much like MatLab to use as possible within Haskell. - whilst a strongly types matrix library may be better, there are difficaulties to implementing this in haskell, but this does not remove the need for a matrix library. (sometimes the best is the enemy of the good!) Keean.

Bulat Ziganshin

10 Jun 10 Jun

1:19 a.m.

New subject: Re[2]: Haskell Matrix Library...

Hello Keean, Thursday, June 09, 2005, 8:52:55 PM, you wrote: KS> (sometimes the best is the enemy of the good!) i think that best in this case will be splitting library to "low-level" and "high-level" parts, where low-level part implements functions like "miltipleMatrixByScalar", "miltipleMatrixByMatrix" and so on and high-level part implements matlab-like or strong-typed access to this functions via appropriate overloaded operators and so on. so in this case low-level==semantics and high-level==syntax moreover, low-level interface must be general enough to allow different implementations - haskell, template haskell, ffi->lipack -- Best regards, Bulat mailto:bulatz@HotPOP.com

John Meacham

1:30 a.m.

New subject: Haskell Matrix Library...

On Fri, Jun 10, 2005 at 09:19:40AM +0400, Bulat Ziganshin wrote:

...

Hello Keean,

Thursday, June 09, 2005, 8:52:55 PM, you wrote:

KS> (sometimes the best is the enemy of the good!)

i think that best in this case will be splitting library to "low-level" and "high-level" parts, where low-level part implements functions like "miltipleMatrixByScalar", "miltipleMatrixByMatrix" and so on and high-level part implements matlab-like or strong-typed access to this functions via appropriate overloaded operators and so on. so in this case low-level==semantics and high-level==syntax

perhaps the low-level part could be the LAPACK/BLAS interface that was floating around? John -- John Meacham - ⑆repetae.net⑆john⑈

Keean Schupke

6:33 a.m.

New subject: Haskell Matrix Library...

Just looked at the matlab code... It doesn't in general do too much overloading - the only obvious example I can find is with matrix multiplication '*'. Here it multiplies matrices normally unless either of the arguments is a scalar, in which case the scalar is multipled by each element of the other matrix. The other multiply operator '.*' multiplies each element in one matix by the same element in the other matrix (as in addition, so matrices must be the same sizes)... If I were going to implement everything in Haskell I would probably use a data type like: data Matrix a e = Matrix [Int] (a Int e) type UMatrix e = Matrix UArray e Where [Int] is a list of dimensions, ie for a scalar [ ], a vector [3], 2D matrix [3,4], or a higher dimensional matrix [3,4,5,6]. It is possible to get static typing using type level naturals... something like: data Matrix d a e = Matrix d (a Int e) The use an HList (or equivalent) to encode the dimensionality of the matrix in the type. For example define [D0 x .. D9 x, DNil] as decimal type constructors, you would then get: type Matrix2x3 a e = Matrix (HCons (D2 DNil) (HCons (D3 DNil) HNil)) a Int e However as the built in operators +,-,*,/ etc have the type "a -> a -> a" you could not define '*' to mean matrix multiplication. At the moment I favour using dynamic size checking on the matrices, and overloading all the standard Haskell Numeric operators for matrices. I seems easier to use FFI and LAPACK for the more complex matrix operations... anyone know how to marshall data for Fortran? Keean. Bulat Ziganshin wrote:

...

Hello Keean,

Thursday, June 09, 2005, 8:52:55 PM, you wrote:

KS> (sometimes the best is the enemy of the good!)

i think that best in this case will be splitting library to "low-level" and "high-level" parts, where low-level part implements functions like "miltipleMatrixByScalar", "miltipleMatrixByMatrix" and so on and high-level part implements matlab-like or strong-typed access to this functions via appropriate overloaded operators and so on. so in this case low-level==semantics and high-level==syntax

moreover, low-level interface must be general enough to allow different implementations - haskell, template haskell, ffi->lipack

Keean Schupke

7:11 a.m.

New subject: Haskell Matrix Library...

I think people misunderstood this quote... The idea is that you may devote a lot of time to developing the 'Best' solution to a problem, but never get to use the solution (because it took so long to get right). On the other hand you can develop a 'good' solution much quicker, and have time to use it. I guess it is a rephrasing of the classic "95% is good enough" rule. In other words I think it is better to have a usable, but imperfect, matrix library soon, rather than a perfect one at some time in the future. Of course this library does not have to be included in the standard libraries, but if it is useful enough and is reasobably clean I don't see why it would not be wanted. LAPACK does not appear to support higher dimensional matrices... is this right? Anyone got any pointers to definitions for transpose and multiply for arbitrary higher dimensional matrices? Keean.

...

...
Thursday, June 09, 2005, 8:52:55 PM, you wrote:

KS> (sometimes the best is the enemy of the good!)

Henning Thielemann

7:36 a.m.

New subject: Haskell Matrix Library...

On Fri, 10 Jun 2005, Keean Schupke wrote:

...

I think people misunderstood this quote... The idea is that you may devote a lot of time to developing the 'Best' solution to a problem, but never get to use the solution (because it took so long to get right). On the other hand you can develop a 'good' solution much quicker, and have time to use it. I guess it is a rephrasing of the classic "95% is good enough" rule.

I gave you suggestions that can be implemented immediately. They resulted from my bad experiences with MatLab and better experiences I made with the suggested approach in Modula-3. Intentionally re-implementing a bad thing is different from implementing something you still don't know enough about.

Henning Thielemann

7:51 a.m.

New subject: Haskell Matrix Library...

On Fri, 10 Jun 2005, Henning Thielemann wrote:

...

On Fri, 10 Jun 2005, Keean Schupke wrote:

...
I think people misunderstood this quote... The idea is that you may devote a lot of time to developing the 'Best' solution to a problem, but never get to use the solution (because it took so long to get right). On the other hand you can develop a 'good' solution much quicker, and have time to use it. I guess it is a rephrasing of the classic "95% is good enough" rule.

I gave you suggestions that can be implemented immediately. They resulted from my bad experiences with MatLab and better experiences I made with the suggested approach in Modula-3. Intentionally re-implementing a bad thing is different from implementing something you still don't know enough about.

... as I also wrote, there are already at least two libraries providing some of the functionality you need. Do we really need a new library whose sloppy design is intended?

Henning Thielemann

7:57 a.m.

New subject: Haskell Matrix Library...

On Fri, 10 Jun 2005, Keean Schupke wrote:

...

LAPACK does not appear to support higher dimensional matrices... is this right?

Yes, why should it do so? Every finite dimensional operator can be represented by a matrix.

...

Anyone got any pointers to definitions for transpose and multiply for arbitrary higher dimensional matrices?

The generalised transposition should be able to arbitrarily permute the dimensions of a tensor. A generalised multiplication may fuse the last dimension of the first tensor with the first dimension with the second tensor. E.g. dimensions [a,b,c,d] and [d,e,f,g] result in [a,b,c,e,f,g].

Keean Schupke

10:56 a.m.

New subject: Haskell Matrix Library...

Henning Thielemann wrote:

...

I gave you suggestions that can be implemented immediately. They resulted from my bad experiences with MatLab and better experiences I made with the suggested approach in Modula-3. Intentionally re-implementing a bad thing is different from implementing something you still don't know enough about.

Jan Skibinski's page seems to have disappeared so there is no source I can find for the Orthoganal stuff - besides which I am interested in efficency - and that basically said the algorithms were experimental and use Haskell lists... The HaskellDSP defines simple matrix operations, but does not put the operations into the numeric, fractional, and floating classes. This makes the resulting expressions hard to read. I can appreciate what you are saying - but I remain to be persuaded.

...

... as I also wrote, there are already at least two libraries providing

...

some of the functionality you need. Do we really need a new library whose sloppy design is intended?

Well, I do as my requirements are different than those of the libraries above. I am merely offering to do a little extra work and put that into a generic library if people want it. If people don't want it, thats fine...

...

...
LAPACK does not appear to support higher dimensional matrices... is this right?

Yes, why should it do so? Every finite dimensional operator can be represented by a matrix.

I have found a couple of papers on generalising matrix algebra to higher dimensional matrices... So it seems to be possible... Keean.

Henning Thielemann

26 Jun 26 Jun

1:22 p.m.

New subject: Haskell Matrix Library...

On Fri, 10 Jun 2005, Keean Schupke wrote:

...

Jan Skibinski's page seems to have disappeared so there is no source I can find for the Orthoganal stuff - besides which I am interested in efficency - and that basically said the algorithms were experimental and use Haskell lists...

I have some private copy: http://www.math.uni-bremen.de/~thielema/Research/LinearAlgebra/

...

The HaskellDSP defines simple matrix operations, but does not put the operations into the numeric, fractional, and floating classes. This makes the resulting expressions hard to read.

Numeric classes do not exist for readability in the sense of infix operators. If infix operators were the main reason simple operators would suffice. The reason why these operators are methods of classes is that it shall be possible to write mathematical expressions which can be universally used for different types of numbers, such as Float, Double, Rational, Integer and so on. Prelude's class framework is oriented at the scalar types it provides. My concern is: Are matrices a natural extension of these scalar types? What is (fromInteger 1) for the matrix class? The identity matrix or the matrix which consists entirely of ones? Are there non-trivial algorithms written for numbers which can be re-used for matrices without re-writing? If not, why pressing matrices into this framework? You could also define (+), (-), (*) for Strings. But what is the benefit? Are there numerical algorithms which can be re-used for Strings? (Now when I write it, I remember the algorithm for searching a dictionary of strings by the method of false position. But the Prelude numeric class framework wouldn't help here, too, since you need a division of type String -> String -> Double, or so.) If readability in the sense of "writing less" is your main goal then just define infix operators for the functions found in HaskellDSP. (*>) is a good choice for scaling a vector or a matrix, this symbol is already used in numericprelude. You may define (.*) for element-wise multiplication which is already known from MatLab. It is more readable than everywhere (*) because it tells the user something about the types of the operands.

...

...
...
LAPACK does not appear to support higher dimensional matrices... is this right?

Yes, why should it do so? Every finite dimensional operator can be represented by a matrix.

I have found a couple of papers on generalising matrix algebra to higher dimensional matrices... So it seems to be possible...

Certainly, but LAPACK is very low-level. It's your task to break operations on various types of vectors down to basic linear algebra. The point is that every vector from a finite dimensional space can be represented by an array and each linear operator can be represented by a matrix. Everything which is structured differently can be converted to this but it is your task not that of LAPACK or any other linear algebra package.

Benjamin Franksen

2:44 p.m.

New subject: Haskell Matrix Library...

On Sunday 26 June 2005 19:22, Henning Thielemann wrote:

...

The reason why these operators are methods of classes is that it shall be possible to write mathematical expressions which can be universally used for different types of numbers, such as Float, Double, Rational, Integer and so on. Prelude's class framework is oriented at the scalar types it provides. My concern is: Are matrices a natural extension of these scalar types? What is (fromInteger 1) for the matrix class? The identity matrix or the matrix which consists entirely of ones?

The identity matrix, of course. For a given dimension, numeric literals of value x should be converted as (x .* m_id), where (.*) is scalar multiplication, and m_id is the identity matrix. This is BTW an embedding of the number field into the (quadratic) matrices (over teh numbers in question) as a subring. The real problem is the dimension, which we would like to have as an extra type argument and this is incompatible with the Prelude. Ben

Henning Thielemann

3:25 p.m.

New subject: Haskell Matrix Library...

On Sun, 26 Jun 2005, Benjamin Franksen wrote:

...

On Sunday 26 June 2005 19:22, Henning Thielemann wrote:

...
The reason why these operators are methods of classes is that it shall be possible to write mathematical expressions which can be universally used for different types of numbers, such as Float, Double, Rational, Integer and so on. Prelude's class framework is oriented at the scalar types it provides. My concern is: Are matrices a natural extension of these scalar types? What is (fromInteger 1) for the matrix class? The identity matrix or the matrix which consists entirely of ones?

The identity matrix, of course.

Certainly not "of course", because e.g. 'exp' and 'sin' are "of course" applied element-wise in the library proposal. Thus e.g. exp (fromScalar 0) == fromScalar (exp 0) is no longer true. (fromInteger does not work for this example, but you see the problem in principle.)

...

The real problem is the dimension, which we would like to have as an extra type argument and this is incompatible with the Prelude.

Btw. it is not good to code a fixed matrix size into the type because there are many applications where the size is not known at compile time. The only thing we should assert is that the sizes of matrices in an operation fit together. A general type signature like matrix_mul :: Num a => Matrix i j a -> Matrix j k a -> Matrix i k a is a step towards this goal. Making it really type-safe sounds like the implicit configuration problem. http://www.eecs.harvard.edu/~ccshan/prepose/prepose.pdf

Benjamin Franksen

8:39 p.m.

New subject: Haskell Matrix Library...

On Sunday 26 June 2005 21:25, Henning Thielemann wrote:

...

On Sun, 26 Jun 2005, Benjamin Franksen wrote:

...
On Sunday 26 June 2005 19:22, Henning Thielemann wrote:

...
The reason why these operators are methods of classes is that it shall be possible to write mathematical expressions which can be universally used for different types of numbers, such as Float, Double, Rational, Integer and so on. Prelude's class framework is oriented at the scalar types it provides. My concern is: Are matrices a natural extension of these scalar types? What is (fromInteger 1) for the matrix class? The identity matrix or the matrix which consists entirely of ones?

The identity matrix, of course.

Certainly not "of course", because e.g. 'exp' and 'sin' are "of course" applied element-wise in the library proposal.

I wasn't refering to the library proposal. And element-wise definitions of 'exp' and 'sin' for matrices may be natural in certain contexts, but certainly not as /the/ general definition.

...

Thus e.g. exp (fromScalar 0) == fromScalar (exp 0) is no longer true. (fromInteger does not work for this example, but you see the problem in principle.)

Yea, this is ugly and inconsistent. I'd rather use explicit lifting to express that a numeric function should be applied element-wise. BTW, isn't Matrix a functor over Num so we can use fmap for lifting? Ah, doesn't work because of the extra Num constraint, grrr...

...

...
The real problem is the dimension, which we would like to have as an extra type argument and this is incompatible with the Prelude.

Btw. it is not good to code a fixed matrix size into the type because there are many applications where the size is not known at compile time.

Right.

...

The only thing we should assert is that the sizes of matrices in an operation fit together. A general type signature like

matrix_mul :: Num a => Matrix i j a -> Matrix j k a -> Matrix i k a

is a step towards this goal. Making it really type-safe sounds like the implicit configuration problem. http://www.eecs.harvard.edu/~ccshan/prepose/prepose.pdf

Yes, Matrix dimensions is a good example, one where the author's solution isn't too involved, since dimension is a plain integer. Cheers, Ben

Ketil Malde

27 Jun 27 Jun

2:58 a.m.

New subject: Haskell Matrix Library...

Benjamin Franksen writes:

...

...
Certainly not "of course", because e.g. 'exp' and 'sin' are "of course" applied element-wise in the library proposal.

...

I wasn't refering to the library proposal. And element-wise definitions of 'exp' and 'sin' for matrices may be natural in certain contexts, but certainly not as /the/ general definition.

Perhaps I haven't followed this discussion closely enough, but I'm not sure this kind of element-wise functions should be provided at all. Isn't a Functor¹ instance (and the consequent "fmap sin") so much clearer that it is worth the extra typing? -k ¹ Or something similar, if you really must preserve the matrix' type. -- If I haven't seen further, it is by standing in the footprints of giants

Henning Thielemann

7:05 a.m.

New subject: Haskell Matrix Library...

On Mon, 27 Jun 2005, Ketil Malde wrote:

...

Benjamin Franksen writes:

...
I wasn't refering to the library proposal. And element-wise definitions of 'exp' and 'sin' for matrices may be natural in certain contexts, but certainly not as /the/ general definition.

Perhaps I haven't followed this discussion closely enough, but I'm not sure this kind of element-wise functions should be provided at all. Isn't a Functor� instance (and the consequent "fmap sin") so much clearer that it is worth the extra typing?

I also plead for a 'map' function instead of a Floating instance.

Marcin 'Qrczak' Kowalczyk

7:10 a.m.

New subject: Haskell Matrix Library...

Benjamin Franksen writes:

...

BTW, isn't Matrix a functor over Num so we can use fmap for lifting? Ah, doesn't work because of the extra Num constraint, grrr...

I see no need of putting the Num constraint on the type definition. fmap will work. -- __("< Marcin Kowalczyk \__/ qrczak@knm.org.pl ^^ http://qrnik.knm.org.pl/~qrczak/

Keean Schupke

28 Jun 28 Jun

6:22 a.m.

New subject: Haskell Matrix Library...

Marcin 'Qrczak' Kowalczyk wrote:

...

Benjamin Franksen writes:

...
BTW, isn't Matrix a functor over Num so we can use fmap for lifting? Ah, doesn't work because of the extra Num constraint, grrr...

I see no need of putting the Num constraint on the type definition. fmap will work.

Hmm... It appears I was being dense (or a perfectionist) as I couldn't define a generic instance of Functor for all Matrices generalised over Array type. Of course the solution is just to define an instance for each array type... Here's the latest version ... now smaller with the Floating instance removed, but an instance of Functor added, and a better Show instance (that generates valid Haskell using the constructor functions provided)... Have also removed the Matlab style "scalar" by "matrix" operations, so to multiply a scalar by each element use fmap (or mmap) fmap (+1) matrix I guess in the end it's just more Haskell like - even if it is harder to read than the Matlab code. Keean.

Keean Schupke

6:39 a.m.

New subject: Haskell Matrix Library...

Further to my last email, It appears that the instances of Functor for UArray and IOUArray don't work... Also I think sin and cos could be defined for a matrix using the matrix version of 'exp' this means that Floating could be defined for matrices - but using proper matrix operations... I might do this eventually. Keean.

Jan-Willem Maessen

10 Jun 10 Jun

10:50 a.m.

New subject: Haskell Matrix Library...

On Jun 10, 2005, at 6:33 AM, Keean Schupke wrote:

...

[Making Matrices work in Haskell]...

I seems easier to use FFI and LAPACK for the more complex matrix operations...

And it'll be tough to impossible to match, say, FFTW and Atlas in Haskell. FFTW at least uses the C compiler like an assembler.

...

anyone know how to marshall data for Fortran?

Find something with a C interface and use that. Most LAPACK and BLAS libraries should have a C binding. My understanding is that Fortran dope vector formats are not standard, though the latest versions of the Fortran standard aim to fix that non-portability problem. -Jan-Willem Maessen

...

Keean.

PS - I had the following code kicking around, written with the hope of specializing the Matrices for maximal unboxing (thus the use of constructors with explicit type signatures, which may not be strictly necessary anymore). Is this useful to anyone? No FFI here, it's all just Haskell code, but using an FFI binding instead wouldn't be hard. I'd been hoping to push more towards specialized matrix representations (sparse, banded, upper/lower triangular), to see how ugly the code would get. And of course there should be "solve" and "decompose" routines of various sorts.

David Roundy

11 Jun 11 Jun

10:30 a.m.

New subject: Haskell Matrix Library...

On Fri, Jun 10, 2005 at 10:50:01AM -0400, Jan-Willem Maessen wrote:

...

On Jun 10, 2005, at 6:33 AM, Keean Schupke wrote:

...
anyone know how to marshall data for Fortran?

Find something with a C interface and use that. Most LAPACK and BLAS libraries should have a C binding. My understanding is that Fortran dope vector formats are not standard, though the latest versions of the Fortran standard aim to fix that non-portability problem.

I've never heard of a problem with the Fortran matrix format. We call the fortran LAPACK/BLAS routines straight from C, using an autoconf macro to figure out how the fortran symbols are mangled. It's not so hard. Basically, a fortran matrix is just a C matrix (2D array) transposed. It's either M_ij = M[i+N*j] or M[j+N*i]. I can never remember which. Fortran also always passes by reference, so when calling from C you need to make every argument a pointer, but other than that you can pretty much read the LAPACK docs and code in C straight from that. Oddly enough, I think calling the fortran routines from C is more portable than trying to find a C binding. P.S. for examples of calling lapack from C/C++ (although it also does a lot more) you can run see the dft++ source code. darcs get http://dft.physics.cornell.edu/dft You can grep for F77_FUNC to find fortran-calling code. The fortran-calling and lapack-finding autoconf macros are cribbed from mpb, which is written by the Steven Johnson, who is one of the authors of fftw. So I suppose you could also look at the mpb source... -- David Roundy http://www.darcs.net

Keean Schupke

12 Jun 12 Jun

3 p.m.

New subject: Haskell Matrix Library...

Well, here for your delight is my first attempt at what some will undoutably think is a broken matrix library. However I would argue that my principal aim is to allow functions to be polymorphic over matrices as well as scalars. Features are: - Scalars and Vectors supported as [1x1] and [Nx1] matrices. - Standard arithmetic operators do standard matrix operations (except "/" is unimplemented - until I write a Gaussian-elimination routine) - 'dot' operators '.*' './' do elementwise operations on matrices, and are composed into a matrix level analogue of the standard Numeric class hierachy. - Any operation involving a scalar applies that operation to each element in the matrix. - operations in Floating that have no obvious matrix equivalents are applied elementwise (sin/cos/tan) - operations in Floating that have matrix equivalents (exp) are applied elementwise for consistancy with the other functions in the Floating class. So functions which generalise over scalars, vectors and matrices can be written... for example element-wise generalisation: f :: FloatingMatrix a i e => Matrix a i e -> Matrix a i e f x = 1./ (1 + exp (-x)) can be used with: f (scalar 2) f (vector [1,2,3]) f (matrix [[1,2,3],[3,4,5],[6,7,8]]) On reflection, am not that fussy over the naming of functions. 'exp' could be matrix exponentiation rather than the element-wise one. It could be that just defining a map for matrices is enough, in which case 'f' could be written: f :: FloatingMatrix a i e => Matrix a i e -> Matrix a i e f z = mmap (\x -> 1 / (1 + exp (-x))) z I will add features as I require them, but if anyone actually wants to use it, I will consider patches and requests for features... Idea's and improvements greatly appreciated... The archive includes the module, and a modified copy of lambdaTeX which uses smaller fonts, more appropriate to modern high resolution printers and displays. Keean.

David Roundy

13 Jun 13 Jun

7:11 p.m.

New subject: Haskell Matrix Library...

On Sun, Jun 12, 2005 at 08:00:02PM +0100, Keean Schupke wrote:

...

- Standard arithmetic operators do standard matrix operations (except "/" is unimplemented - until I write a Gaussian-elimination routine)

It would be nice to at least implement "/" for division by scalars... yes one could use ./ for that, but it's less nice.

...

- 'dot' operators '.*' './' do elementwise operations on matrices, and are composed into a matrix level analogue of the standard Numeric class hierachy. - Any operation involving a scalar applies that operation to each element in the matrix.

Looks nice.

...

- operations in Floating that have no obvious matrix equivalents are applied elementwise (sin/cos/tan)

At least sin and cos do have matrix equivalents (imag(mexp(M)), etc)...

...

- operations in Floating that have matrix equivalents (exp) are applied elementwise for consistancy with the other functions in the Floating class.

Okay... I don't really expect to use mexp anyways (although it *is* handy in quantum mechanics, if you can store the entire hamiltonian). Shouldn't this be an instance of Functor, and have the mmap be called fmap? I'd would lean towards making 'exp' be matrix exponentiation.

...

I will add features as I require them, but if anyone actually wants to use it, I will consider patches and requests for features...

I don't suppose there's a chance you could set up a public darcs repository for this?

...

Idea's and improvements greatly appreciated...

One thing that would really help would be a nicer show instance. As a corollary, a read instance would also be nice. -- David Roundy

Keean Schupke

14 Jun 14 Jun

7:54 a.m.

New subject: Haskell Matrix Library...

David Roundy wrote:

...

On Sun, Jun 12, 2005 at 08:00:02PM +0100, Keean Schupke wrote:

...
- Standard arithmetic operators do standard matrix operations (except "/" is unimplemented - until I write a Gaussian-elimination routine)

It would be nice to at least implement "/" for division by scalars... yes one could use ./ for that, but it's less nice.

Can be done... matrix division also not that hard...

...

...
- 'dot' operators '.*' './' do elementwise operations on matrices, and are composed into a matrix level analogue of the standard Numeric class hierachy. - Any operation involving a scalar applies that operation to each element in the matrix.

Looks nice.

...
- operations in Floating that have no obvious matrix equivalents are applied elementwise (sin/cos/tan)

At least sin and cos do have matrix equivalents (imag(mexp(M)), etc)...

You learn something everyday... I guess I knew this already as I know about the power series definitions of sin and cos... guess I just didn't think about it too much. Would we rather have: sin and sinm (elementwise and matrix) - same as matlab OR sine and sin (elementwise and matrix) - applies matrix form as 'normal' OR some other.

...

...
- operations in Floating that have matrix equivalents (exp) are applied elementwise for consistancy with the other functions in the Floating class.

Okay... I don't really expect to use mexp anyways (although it *is* handy in quantum mechanics, if you can store the entire hamiltonian).

Shouldn't this be an instance of Functor, and have the mmap be called fmap? I'd would lean towards making 'exp' be matrix exponentiation.

Ah... No. class Functor m where fmap :: a -> b -> m a -> m b But the type of "mmap" is mmap :: e -> e -> Matrix a i e -> Matrix a i e in other words we map values to the same type, not a different type.

...

...
I will add features as I require them, but if anyone actually wants to use it, I will consider patches and requests for features...

I don't suppose there's a chance you could set up a public darcs repository for this?

Yes ... will do.

...

...
Idea's and improvements greatly appreciated...

One thing that would really help would be a nicer show instance. As a corollary, a read instance would also be nice.

Definitely, had planned to do this. Keean.

Bulat Ziganshin

8 a.m.

New subject: Re[2]: Haskell Matrix Library...

Hello Keean, Tuesday, June 14, 2005, 3:54:08 PM, you wrote: KS> class Functor m where KS> fmap :: a -> b -> m a -> m b KS> But the type of "mmap" is KS> mmap :: e -> e -> Matrix a i e -> Matrix a i e KS> in other words we map values to the same type, not a different type. mmap should be mmap :: (e -> e1) -> Matrix a i e -> Matrix a i e1 for example, we can map matrix of strings to matrix of its lengths -- Best regards, Bulat mailto:bulatz@HotPOP.com

Keean Schupke

9:58 a.m.

New subject: Haskell Matrix Library...

Yes, it should be - I have corrected: mmap :: (IArray a e, Num e, IArray a e', Num e') => (e -> e') -> Matrix a e -> Matrix a e' It still cannot be an instance of functor though as: instance Functor (Matrix a) where fmap f = mmap x We cannot constrain the instance correctly... We need "IArray a e", but fmap is (expanding) fmap :: b -> c -> Matrix a b -> Matrix a c But 'b' and 'c' cannot be constrained by the instance... because they do not occur in the class, only the function definition. Keean. Bulat Ziganshin wrote:

...

Hello Keean,

Tuesday, June 14, 2005, 3:54:08 PM, you wrote:

KS> class Functor m where KS> fmap :: a -> b -> m a -> m b

KS> But the type of "mmap" is

KS> mmap :: e -> e -> Matrix a i e -> Matrix a i e

KS> in other words we map values to the same type, not a different type.

mmap should be

mmap :: (e -> e1) -> Matrix a i e -> Matrix a i e1

for example, we can map matrix of strings to matrix of its lengths

Frederik Eaton

24 Jun 24 Jun

10:49 p.m.

New subject: Haskell Matrix Library...

I haven't read this whole thread in detail, but the matlab/octave thing has also been something of a concern of mine since I will probably be doing a master's in machine learning next year. However, I was just planning to spend a few weeks writing a embedded domain specific language wrapper in Haskell for matlab/octave. Maybe that's harder than I'd imagined it would be - but the problem with writing a Haskell matrix library from scratch is that not only do you lose a lot of the functionality which has already been implemented in these programs, but you also lose all the work which has gone into tuning them for speed. They support a wide variety of optimized internal matrix representations which are transparent to the user, special fast linear algebra operations on those matrices, etc. And for the kind of applications that people use matlab and octave for, speed and scalability tend to be important. So an arrangement where you keep the optimized matlab kernel as a backend and just send it commands would seem to be preferable. Plus, with an EDSL you could improve the type system, index matrices over arbitrary Enum-like types (which are translated to and from integers), graft rudimentary tensor support into Octave, etc. Maybe there would also be uses for a native Haskell matrix library but my intuition is that something like what I described above would be more useful to more people. Am I wrong? Frederik -- http://ofb.net/~frederik/

Daan Leijen

11:35 p.m.

New subject: Haskell Matrix Library...

Frederik Eaton wrote:

...

However, I was just planning to spend a few weeks writing a embedded domain specific language wrapper in Haskell for matlab/octave. Maybe that's harder than I'd imagined it would be

By all means: go for it. Writing an EDSL is a lot of fun and normally it leads to nice design -- often nicer than the original API. I am very interested to see the end result if you go through with it. -- Daan Leijen ps. I like your mounting proposal too :-)

David Roundy

25 Jun 25 Jun

7:09 a.m.

New subject: Haskell Matrix Library...

On Fri, Jun 24, 2005 at 07:49:47PM -0700, Frederik Eaton wrote:

...

However, I was just planning to spend a few weeks writing a embedded domain specific language wrapper in Haskell for matlab/octave. Maybe that's harder than I'd imagined it would be - but the problem with writing a Haskell matrix library from scratch is that not only do you lose a lot of the functionality which has already been implemented in these programs, but you also lose all the work which has gone into tuning them for speed. They support a wide variety of optimized internal matrix representations which are transparent to the user, special fast linear algebra operations on those matrices, etc. And for the kind of applications that people use matlab and octave for, speed and scalability tend to be important. So an arrangement where you keep the optimized matlab kernel as a backend and just send it commands would seem to be preferable.

I'd expect that simply calling LAPACK you'd get as good performance as octave for anything that octave does efficiently, and you'd get better performance on operations such as for loops, which are extremely slow in octave--at least, assuming you're not running interactively... but maybe even then. I may be wrong, but I think that LAPACK *is* the optimized matlab kernel. I guess there's a bit of smarts to decide whether octave can get by with real arrays, but if you're using haskell you might prefer to decide by hand whether matrices are real or complex.

...

Plus, with an EDSL you could improve the type system, index matrices over arbitrary Enum-like types (which are translated to and from integers), graft rudimentary tensor support into Octave, etc.

I'd think that all these thing could be done more easily with a native Haskell matrix library, since in that case you'd actually have control over your data. With an octave backend, I don't see how you could even conveniently separate real and complex matrices... maybe that isn't one of your goals, but it would seem to me like an improvement in the type system.

...

Maybe there would also be uses for a native Haskell matrix library but my intuition is that something like what I described above would be more useful to more people. Am I wrong?

Either would be useful, but I'd guess that a simple Haskell matrix library calling LAPACK would be more useful. To seriously replace octave, you'd need plotting functions, but I certainly wouldn't want to call octave, which calls gnuplot--it'd be much more flexible to call gnuplot directly (or use some other plotting plan). -- David Roundy http://www.darcs.net

Henning Thielemann

26 Jun 26 Jun

12:41 p.m.

New subject: Haskell Matrix Library...

On Fri, 24 Jun 2005, Frederik Eaton wrote:

...

Maybe there would also be uses for a native Haskell matrix library but my intuition is that something like what I described above would be more useful to more people. Am I wrong?

It would be nice to not restrict the element type of the matrix. There might be matrices of Integers or polynomials and you want to compute a determinant. LAPACK can't do that. So there should be some framework of type classes where the floating point number instances use LAPACK and the others use Haskell code.

Henning Thielemann

12:39 p.m.

New subject: Haskell Matrix Library...

On Tue, 14 Jun 2005, Bulat Ziganshin wrote:

...

Hello Keean,

Tuesday, June 14, 2005, 3:54:08 PM, you wrote:

KS> class Functor m where KS> fmap :: a -> b -> m a -> m b

KS> But the type of "mmap" is

KS> mmap :: e -> e -> Matrix a i e -> Matrix a i e

KS> in other words we map values to the same type, not a different type.

mmap should be

mmap :: (e -> e1) -> Matrix a i e -> Matrix a i e1

for example, we can map matrix of strings to matrix of its lengths

... or a real valued matrix to a complex valued one.

David Roundy

9 Jun 9 Jun

8:09 a.m.

New subject: Haskell Matrix Library...

On Thu, Jun 09, 2005 at 01:35:26PM +0200, Henning Thielemann wrote:

...

On Thu, 9 Jun 2005, Keean Schupke wrote:

...
- would people want a library like MatLab where matrices are instances of Num/Fractional/Floating and can be used in normal math equations...

I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components.

The thing is that all these things are what make matlab so easy to use. Would you really like to have separate multiplication operators for every combination of scalars, vectors and matrices? It may be possible to come up with a better scheme, but personally I'd be very happy with a library that just does what matlab does, i.e. define (*), (/), (+), (-), (.*), (./), and treat every object as a matrix.

...

What is the natural definition of matrix multiplication? The element-wise multiplication (which is commutative like the multiplication of other Num instances) or the matrix multiplication? I vote for separate functions or infix operators for matrix multiplications, linear equation system solvers, matrix-vector-multiplication, matrix and vector scaling. I wished to have an interface to LAPACK for the floating point types because you can hardly compete with algorithms you write in one afternoon.

I'd tend to say that matlab has a working system that has a number of satisfied users. Inventing a new interface for linear algebra is an interesting problem, but I'd be satisfied with a "matlab with a nice programming language"... a tool for those who just want to use matrices rather than for studying matrices.

...

...
- operations like cos/sin/log can be applied to a matrix (and apply to each element like in MatLab) ...

MatLab must provide these automatisms because it doesn't have proper higher functions. I think the Haskell way is to provide a 'map' for matrices. Otherwise the question is: What is the most natural 'exp' on matrices? The elementwise application of 'exp' to each element or the matrix exponentation?

The utility of these functions would depend on whether one makes matrices be instances of floating and *doesn't* implement a separate operator for the scalar-matrix operations. If that is the case, I think we'd want these functions defined you could do something like let a = sqrt(2.0)*b On the other hand, this argument only says that I would like these to be defined for 1x1 matrices, and doesn't specify what they should do to bigger matrices. Of course, the trancendental functions are only defined (except in an element-wise sense) for square matrices... But it would be an interesting trick to allow exponentiation of a matrix (presumably by diagonalizing, exponentiating the eigenvalues and then applying the unitary transformation back again). And we could afford that, since element-wise exp could always be done with (map exp a)... :) -- David Roundy http://www.darcs.net

Henning Thielemann

8:53 a.m.

New subject: Haskell Matrix Library...

On Thu, 9 Jun 2005, David Roundy wrote:

...

On Thu, Jun 09, 2005 at 01:35:26PM +0200, Henning Thielemann wrote:

...
I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components.

The thing is that all these things are what make matlab so easy to use.

Isn't it easier to actually use MatLab if you like it? I don't see the benefit of a MatLab library embedded in Haskell.

...

Would you really like to have separate multiplication operators for every combination of scalars, vectors and matrices?

There are more functions but I think they are worth the distinction. In the case of (+) and (-) I agree to overloaded operators. In the case of multiplication we have these classes of functions: - multiply a scalar with a vector or a matrix or a general tensor (the (*>) method in class Module in NumericPrelude) - multiply matrix with vector (the matrix is considered as a canonical representation of a linear operator, which is hereby applied to the vector) - multiply matrix with matrix (representing the composition of linear operators) - scalar product (could be part of a Hilbert Space type class) - a general product between tensors were arbitrary ranks are fused

...

It may be possible to come up with a better scheme, but personally I'd be very happy with a library that just does what matlab does, i.e. define (*), (/), (+), (-), (.*), (./), and treat every object as a matrix.

Treating everything as a matrix is a big bug of MatLab: Matrices are a special case of the tensors MatLab also support. So tensors of at most rank 2 are treated as matrices and above they are treated differently. This break is quite arbitrary because you often switch from vectors to arrays of vectors (= matrices) and from matrices to arrays of matrices (= higher rank tensors). In MatLab you must distinguish between row and column vectors, which is again quite arbitrary (although it is adapted from common usage in linear algebra). Further you never know if a 1x1 matrix is actually meant as a scalar, a one-element row vector, a one-element column vector or a one by one matrix. In my applications I can always decide statically which of these four (essentially three) alternatives is appropriate and I want to express this by types and I want that the compiler checks the proper usage.

...

I'd tend to say that matlab has a working system that has a number of satisfied users.

I know also a lot of unsatisfied people - no longer being "users". :-)

...

Inventing a new interface for linear algebra is an interesting problem, but I'd be satisfied with a "matlab with a nice programming language"...

Haskell with weak matrix typing will be at least as broken as MatLab.

Keean Schupke

10:11 a.m.

New subject: Haskell Matrix Library...

Henning Thielemann wrote:

...

On Thu, 9 Jun 2005, David Roundy wrote:

...
On Thu, Jun 09, 2005 at 01:35:26PM +0200, Henning Thielemann wrote:

...
I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components.

The thing is that all these things are what make matlab so easy to use.

Isn't it easier to actually use MatLab if you like it? I don't see the benefit of a MatLab library embedded in Haskell.

Because Haskell is a compiled language - and because I program in Haskell - however every time I try and do something a little complicated I am frustrated by lack of matrix support. (and Haskell is free/open-source)

...

...
Would you really like to have separate multiplication operators for every combination of scalars, vectors and matrices?

There are more functions but I think they are worth the distinction. In the case of (+) and (-) I agree to overloaded operators. In the case of multiplication we have these classes of functions:

- multiply a scalar with a vector or a matrix or a general tensor (the (*>) method in class Module in NumericPrelude) - multiply matrix with vector (the matrix is considered as a canonical representation of a linear operator, which is hereby applied to the vector) - multiply matrix with matrix (representing the composition of linear operators) - scalar product (could be part of a Hilbert Space type class) - a general product between tensors were arbitrary ranks are fused

...
It may be possible to come up with a better scheme, but personally I'd be very happy with a library that just does what matlab does, i.e. define (*), (/), (+), (-), (.*), (./), and treat every object as a matrix.

Treating everything as a matrix is a big bug of MatLab: Matrices are a special case of the tensors MatLab also support. So tensors of at most rank 2 are treated as matrices and above they are treated differently. This break is quite arbitrary because you often switch from vectors to arrays of vectors (= matrices) and from matrices to arrays of matrices (= higher rank tensors). In MatLab you must distinguish between row and column vectors, which is again quite arbitrary (although it is adapted from common usage in linear algebra). Further you never know if a 1x1 matrix is actually meant as a scalar, a one-element row vector, a one-element column vector or a one by one matrix. In my applications I can always decide statically which of these four (essentially three) alternatives is appropriate and I want to express this by types and I want that the compiler checks the proper usage.

Yes, I haven't thought about general tensor support - it would be better - however the current project I am working on does not need tensors (yet!)

...

...
I'd tend to say that matlab has a working system that has a number of satisfied users.

I know also a lot of unsatisfied people - no longer being "users". :-)

I work in the electrical and electronic engineering dept of Imperial College London, MatLab is probably the most used engineering application (both for teaching and research)... Nearly everone used MatLab actively whether its the signal processing toolbox, the neural-network toolbox, or the control toolbox. People want to write equations like they look in maths text-books ... which as it happens use the same symbol for matrix and scalar multiplication...

...

...
Inventing a new interface for linear algebra is an interesting problem, but I'd be satisfied with a "matlab with a nice programming language"...

Haskell with weak matrix typing will be at least as broken as MatLab.

You could also say: Haskell with weak matrix typing will be at least as good as MatLab. Infact as 'broken' and 'good' merely represent opinions what you are really saying is: Haskell with weak matrix typing will be like MatLab... which is a tortology. Haskell cannot do strong matrix typing easily... you need to get into type level natural representations (to encode the dimensions of the matrix into the type)... At this point it is turning into a computer-science project rather than a useful matrix library. Keean.

David Roundy

11 Jun 11 Jun

10:15 a.m.

New subject: Haskell Matrix Library...

On Thu, Jun 09, 2005 at 02:53:27PM +0200, Henning Thielemann wrote:

...

On Thu, 9 Jun 2005, David Roundy wrote:

...
On Thu, Jun 09, 2005 at 01:35:26PM +0200, Henning Thielemann wrote:

...
I'm pretty unhappy with all these automatisms in MatLab which prevent fast detection of mistakes, e.g. treating vectors as column or row matrices, random collapses of singleton dimensions, automatic extension of singletons to vectors with equal components.

The thing is that all these things are what make matlab so easy to use.

Isn't it easier to actually use MatLab if you like it? I don't see the benefit of a MatLab library embedded in Haskell.

I don't like it. But I also know that the rest of my research group can do many things more easily in octave than I can without it. And as a programming tool it's just abominable.

...

...
Would you really like to have separate multiplication operators for every combination of scalars, vectors and matrices?

There are more functions but I think they are worth the distinction. In the case of (+) and (-) I agree to overloaded operators. In the case of multiplication we have these classes of functions:

How can we overload (+) and (-) without also overloading (*)? They're all in the same class.

...

...
It may be possible to come up with a better scheme, but personally I'd be very happy with a library that just does what matlab does, i.e. define (*), (/), (+), (-), (.*), (./), and treat every object as a matrix.

Treating everything as a matrix is a big bug of MatLab: Matrices are a special case of the tensors MatLab also support.

I must admit that my familiarity with matlab is mostly through octave, which doesn't even have as nice tensor support as matlab does (from one I hear). But the point of a simple matrix module would simpley a have matrix support, not to have tensor support. You can put full tensor support into your tensor module. Matrices are indeed a special case of tensors, but they are a special case that has simpler operators than those that exist for tensors, and on which there are more functions defined. The only reasonable notations that I'm aware of for tensor multiplication are either based on repeated index summation (or explicit summation) or treating a tensor as a matrix (e.g. treating a fourth-rank tensor as a 9x9 matrix). If you know of a nice representation for tensor multiplication, I'd be very interested to hear about it. If you don't, then I'm not sure how the creation of O(N^2) multiplication operators (where N is the highest rank tensor supported) can be viewed as an improvement over treating everything as a matrix.

...

...
Inventing a new interface for linear algebra is an interesting problem, but I'd be satisfied with a "matlab with a nice programming language"...

Haskell with weak matrix typing will be at least as broken as MatLab.

Fortunately, the brokenness of matlab that bothers me isn't its matrix operators, but rather it's features as a programming language. -- David Roundy http://www.darcs.net

Henning Thielemann

26 Jun 26 Jun

12:54 p.m.

New subject: Haskell Matrix Library...

On Sat, 11 Jun 2005, David Roundy wrote:

...

How can we overload (+) and (-) without also overloading (*)? They're all in the same class.

Unfortunately yes. This is the reason for the type class system of http://cvs.haskell.org/darcs/numericprelude/ Since this framework is under progress and it may annoy people to use non-standard classes you may to stick to Prelude's Num class which forces you to define (*) whenever you only want (+) and (-). So in this case you could only define (*) = undefined which is not nice because it is only checked at runtime, but the numeric type classes are also not nice. :-)

...

If you know of a nice representation for tensor multiplication, I'd be very interested to hear about it. If you don't, then I'm not sure how the creation of O(N^2) multiplication operators (where N is the highest rank tensor supported) can be viewed as an improvement over treating everything as a matrix.

I would define the types Vector Matrix Tensor (for any rank, including 1 and 2) A rank 1 tensor would be different from a vector and a rank 2 tensor would be different from a matrix. This is seems to be ok, since vectors and matrices belong to one theory and tensors are somehow different.

...

...
...
Inventing a new interface for linear algebra is an interesting problem, but I'd be satisfied with a "matlab with a nice programming language"...

Haskell with weak matrix typing will be at least as broken as MatLab.

Fortunately, the brokenness of matlab that bothers me isn't its matrix operators, but rather it's features as a programming language.

I hope you never had to search for a bug which is caused whenever two automatisms of MatLab collide. E.g. polyder([0,0,0,0,1,0]) is not [0,0,0,0,1] as you might expect but [1]. Thus polyder([0,0,0,0,1,0]) + [1,0,0,0,0] is [2,1,1,1,1] ! This is so because adding a scalar to a vector is interpreted as expanding the scalar to a vector by replication and add it to the vector.

7266

Age (days ago)

7301

Last active (days ago)

List overview

Download

45 comments

14 participants

participants (14)

Adrian Hey
Benjamin Franksen
Bulat Ziganshin
Daan Leijen
David Roundy
Einar Karttunen
Frederik Eaton
Henning Thielemann
Jan-Willem Maessen
John Meacham
Keean Schupke
Ketil Malde
Marcin 'Qrczak' Kowalczyk
Ross Paterson