On Wednesday 29 June 2005 12:31, Henning Thielemann wrote:

On Wed, 29 Jun 2005, Alberto Ruiz wrote:

...

In my work I often need linear algebra algorithms and other numeric computations.

Nice coincidence: http://www.haskell.org//pipermail/libraries/2005-June/003936.html

Wow! It is exactly the same idea! I did not find the above message by Keean in my google searchs when I decided to work on this, it is very recent! After a quick look to the thread I wish I would have followed the discussions... A much more serious work was in progress. At least it was fun and I have learnt a lot of Haskell... I'm going to subscribe immediately to the libraries Haskell list. I stupidly thought that it was only for discussion about existing libraries :( And next I will study Keean's library. Thank you for your message, Alberto

On Wed, 29 Jun 2005, Henning Thielemann wrote:

On Wed, 29 Jun 2005, Alberto Ruiz wrote:

...

On Wednesday 29 June 2005 12:31, Henning Thielemann wrote:

...

On Wed, 29 Jun 2005, Alberto Ruiz wrote:

...

In my work I often need linear algebra algorithms and other numeric computations.

Nice coincidence: http://www.haskell.org//pipermail/libraries/2005-June/003936.html

Wow! It is exactly the same idea! I did not find the above message by Keean in my google searchs when I decided to work on this, it is very recent! After a quick look to the thread I wish I would have followed the discussions... A much more serious work was in progress. At least it was fun and I have learnt a lot of Haskell...

I'm going to subscribe immediately to the libraries Haskell list. I stupidly thought that it was only for discussion about existing libraries :(

From a first look at your approach I like it even more because you explicitly distinguish between Matrix and Vector.

I was wrong, the different names are synonymes for the same type. :-(

I would recommend that you look very closely at the design of the LinearAlgebra package, the Matrix and Vector constructors, and the underlying implementation data-structure rtable() for Maple's implementation of all these ideas. About 250 person-days were spent on just the high-level design, as we knew that the underlying implementation would take a lot longer (it ended up taking about 6 person-years). The main design points were: 1. complete separation of the user-visible types (Matrix, Vector, Array) from the implementation type (rtable) 2. complete separation of the usual functions for user use (matrix addition, multiplication, LUDecomposition, etc) from the actual underlying implementations 3. for efficiency purposes (only), allow access to the low-level details of both data-structures and implementations, but spend no effort on making these easy to use. Making them efficient was the only thing that mattered. Thus the calling sequences for the low-level routines are quite tedious. 4. the high-level commands should do what one expects, but they are expected to pay a (small) efficiency cost for this convenience 5. incorporate all of LAPACK (through ATLAS, raw LAPACK is now too slow), as well as *all* of NAG's LinearAlgebra routines in a completely seamless manner. 6. do #6 for double precision as well as arbitrary precision arithmetic, again seamlessly 7. use as much structure information to do algorithm dispatch as possible 8. Matrices and Arrays are different high-level types, with different defaults. A.B for Matrices is matrix multiplication, A.B for Arrays is component-wise. To serves the purpose of both those who want to do linear algebra and those who want to do signal processing. 9. There are row vectors and column vectors, and these are different types. You get type errors if you mix them incorrectly. Things that follow from that are numerous. For example, for the Matrix constructor there are: 1. 15 different abstract 'shapes' a matrix can have, with the band[n1,n2] 'shape' having two integer parameters 2. 14 different storage formats (element *layouts*) 3. each storage format can come in C or Fortran order 4. datatype which can be complex[8], float[8], integer[1], integer[2], integer[4], integer[8] or Maple 5. a variety of attributes (like positive_definite) which can be used by various routines for algorithm selection This implies that, in Maple (unlike in Matlab), the eigenvalues of a real symmetric matrix are guaranteed to be real, ie not contain spurious small imaginary parts. Solving Ax=b for A positive definite will use a Cholesky decomposition (automatically), etc. Computations on arbitrary length integers will also use different algorithms than for matrices of polynomials, for example. Where appropriate, fraction-free algorithms are used, etc. Some routines are even designed to be a bit 'lazy': one can use an LUDecomposition of a matrix A and return essentially only the information necessary to then do a large number of Ax=b solves for different b. Matrix access is also very convenient, at the element level, but also at the sub-matrix level, rather like Matlab's notation. If A is a 5x5 Matrix, then A[1..2,3..4] := A[2..3,-2..-1]; will assign the submatrix of A consisting of the elements from row 2,3 and columns 4,5 to the submatrix of A in rows 1,2 and columns 3,4. Notice the indexing from the right (positive numbers) and the left (negative). One can also say A[[2,1],3..4] := A[2..3,[-1,-2]]; As this is somewhat difficult to explain, let me show what this does:

A := Matrix(5,5,(i,j)->i*j);

[1 2 3 4 5] [ ] [2 4 6 8 10] [ ] A := [3 6 9 12 15] [ ] [4 8 12 16 20] [ ] [5 10 15 20 25]

A[[2,1],3..4] := A[2..3,[-1,-2]]: A;

[1 2 15 12 5] [ ] [2 4 10 8 10] [ ] [3 6 9 12 15] [ ] [4 8 12 16 20] [ ] [5 10 15 20 25] (hopefully the above will come out in a fixed-width font!). This makes writing block-algorithms extremely easy. Overall, this design lets one easily add new algorithms without having to worry about breaking user code. Jacques

On Wed, 29 Jun 2005, Jacques Carette wrote:

Distinction of row and column vectors is a misconcept

Row and column vectors are sometimes worth distinguishing because they can represent entirely different types of object. For example, if a column vector represents an element of a vector space V over a field F, then row vectors can be used to represent elements of the dual space, V* = {f:V->F, f linear}. Quite different objects and in some applications it makes sense to distinguish them.

You would never try to transpose a function.

Funny you say that, that's exactly what I've been writing Haskell code to do, at least for linear functions. Given any linear function f:V->W then the dual function, written f* maps W*->V*. In effect this is the transpose and it makes perfect sense.

So what is the operation of transposition? A type conversion?

It maps objects of type V to objects of type V* and objects of type V->W to objects of type W*->V*. -- Dan

On Wednesday 29 June 2005 22:54, Henning Thielemann wrote:

On Wed, 29 Jun 2005, Dan Piponi wrote:

...

...

On Wed, 29 Jun 2005, Jacques Carette wrote:

Distinction of row and column vectors is a misconcept

Row and column vectors are sometimes worth distinguishing because they can represent entirely different types of object. For example, if a column vector represents an element of a vector space V over a field F, then row vectors can be used to represent elements of the dual space, V* = {f:V->F, f linear}. Quite different objects and in some applications it makes sense to distinguish them.

If f is a function f :: a -> a then, of course, the dual space is again a vector space. But it contains functionals and they have very different type, namely f' :: (a -> a) -> a If dual spaces would represent the concept of transposition then the dual space of the dual space should be original space. It is not, it can even not be identified in many cases, only if the space is reflexive. f'' :: ((a -> a) -> a) -> a has certainly a type very different from (a -> a)!

IIRC, finite-dimensional spaces are always reflexive, as witnessed by the identification of elements of the dual space f with the vector y in f x = (real scalars, here), where the identification is, of course, with respect to a given basis. I bet you know all this very well! Ben

On row & column vectors, do you really want to think of them as {1,...,n)->R? They often represent linear maps from R^n to R or R to R^n, which are very different types. Similarly, instead of working with matrices, how about linear maps from R^n to R^m? In this view, column and row vectors, matrices, and often scalars are useful as representations of linear maps. I've played around some with this idea of working with linear maps instead of the common representations, especially in the context of derivatives (including higher-dimensional and higher-order), where it is the view of calculus on manifolds. It's a lovely, unifying approach and combines all of the various chain rules into a single one. I'd love to explore it more thoroughly with one or more collaborators. Cheers, - Conal -----Original Message----- Henning Thielemann wrote:

On Wed, 29 Jun 2005, Jacques Carette wrote:

...

9. There are row vectors and column vectors, and these are different types. You get type errors if you mix them incorrectly.

What do you mean with "row vectors and column vectors are different types"? Do you mean that in a well designed library they should be distinguished? I disagree to this point of view which is also represented by MatLab.

This was hotly debated during our design sessions too. It does appear to be a rather odd decision, I agree! It was arrived at after trying many alternatives, and finding all of them wanting. Mathematically, vectors are elements of R^n which does not have any inherent 'orientation'. They are just as simply abstracted as functions from {1,...,n) -> R, as you point out. But matrices are then just linear operators on R^n, and if you go there, then you have to first specify a basis for R^n before you can write down a representation for any matrix. Not useful. However, when you step away from the pure mathematics point of view and instead look at actual mathematical usage, things change significantly. Mathematics is full of abuse of notation, and to make usage of matrices and vectors both convenient yet 'safe' required us to re-examine these many de facto conventions that mathematicians routinely use. And try to adapt our design to find the middle road between safety and usefulness. One of the hacks in Matlab, because it only has matrices, is that vectors are also 1xN and Nx1 matrices. In Maple, these are different types. This helps a lot, because then a dot product is a function from R^n x R^n -> R. 1x1 matrices are not reals in Maple. Nor are Nx1 matrices Vectors.

If we instead distinguish row and column vectors because we treat them as matrices, then the quadratic form x^T * A * x denotes a 1x1 matrix, not a real.

But if you consider x to be a vector without orientation, writing down x^T is *completely meaningless*! If x is oriented, then x^T makes sense. Also, if x is oriented, then x^T * (A * x) = (x^T * A) * x. What is the meaning of (x * A) for a 'vector' x ? It gets much worse when the middle A is of the form B.C. To ensure that everything is as associative as can be, but no more, is very difficult. I don't have the time to go into all the details of why this design 'works' (and it does!), giving the 'expected' result to all meaningful linear algebra operations, but let me just say that this was the result of long and intense debate, and was the *only* design that actually allowed us to translate all of linear algebra idioms into convenient notation. Please believe me that this design was not arrived at lightly! Jacques _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

Let me a bit elaborate on what I wrote yesterday. On Wed, 29 Jun 2005, Henning Thielemann wrote:

I think matrices and derivatives are very different issues. I have often seen that the first derivative is considered as vector, and the second derivative is considered as matrix. In this spirit it is used like x^T * (D2 f)(x) * x

It should be h^T * D2 f x * h (if at all the transposition notation is used) if the curvature of f at x is meant.

but this is only abuse of the common multiplication definitions. A good interpretation and notation should seamless extend to higher derivatives. But the interpretation above does not work in higher dimensions.

I like the following type for derivation. derive :: ((i -> a) -> b) -> ((i -> a) -> (i -> b)) Here i is the index type, (i -> a) is the vector type, b is the type the vector function maps to. Its derivative has the same type of argument, but the result is a vector with indices of type i. You see that it is easy to repeat the application of 'derive', just replace b by say i->b. The second derivative yields vectors of type (i -> i -> b). This can be interpreted as matrix because it has two indices. But this is certainly not a matrix which represents a linear mapping as usual, but it is a matrix representing a bilinear form. The only thing we need is a multiplication to reduce one level of indices. mul :: (i -> c) -> (i -> b) -> b Though, what we still need is a general (overloaded?) definition of the scaling of b by c and a sum of b.

That is b must be a vector space with respect to c. 'derive' in this form is a bundled partial derivative, but it can be identified with the total derivative. That is derive f x i is the partial derivative of f at x with respect to the i-th component derive f x is the total derivative of f at x.

On Wed, 29 Jun 2005, Jacques Carette wrote:

...

If we instead distinguish row and column vectors because we treat them as matrices, then the quadratic form x^T * A * x denotes a 1x1 matrix, not a real.

But if you consider x to be a vector without orientation, writing down x^T is *completely meaningless*!

That's what I stated.

If x is oriented, then x^T makes sense. Also, if x is oriented, then x^T * (A * x) = (x^T * A) * x. What is the meaning of (x * A) for a 'vector' x ?

It has of course no meaning. Mathematical notation has the problem that it doesn't distinguish between things that are different but in turn discriminates things which are essentially the same. If your design goal is to keep as close as possible to common notational habits you have already lost! As I already pointed out in an earlier discussion I see it the other way round: Computer languages are the touchstones for mathematical notation because you can't tell a computer about an imprecise expression: "Don't be stupid, you _know_ what I mean!" More specific: You give two different things the same name. You write A*B and you mean a matrix multiplication. Matrix multiplication means finding a representation of the composition of the operators represented by A and B. But you also write A*x and you mean the matrix-vector multiplication. This corresponds to the application of the operator represented by A to the vector x. You see: Two different things, but the same sign (*). Why? You like this ambiguity because of its conciseness. You are used to it. What else? But you have to introduce an orientation of vectors, thus you discriminate two things (row and column vectors) which are essentially the same! What is the overall benefit? It seems to me like the effort of most Haskell XML libraries in trying to have as few as possible combinator functions (namely one: o) which forces you to not discriminate the types for the functions to be combined (the three essential types are unified to a -> [a]) and even more it forces you to put conversions from the natural type (like a->Bool, a->a) in every atomic function!

It gets much worse when the middle A is of the form B.C. To ensure that everything is as associative as can be, but no more, is very difficult.

I don't see the problem. There are three very different kinds of multiplication, they should also have their own signs: Scalar product, matrix-vector multiplication, matrix-matrix multiplication.

I don't have the time to go into all the details of why this design 'works' (and it does!),

I have worked with Maple and I have finally dropped it because of its design. I dropped MatLab, too, because the distinction of row and column vectors, because it makes no sense to distinguish between e.g. convolving row vectors or column vectors. Many routines have to be aware of this difference for which it is irrelevant and many routines work only with one of these kinds and you are often busy with transposing them.

giving the 'expected' result to all meaningful linear algebra operations, but let me just say that this was the result of long and intense debate, and was the *only* design that actually allowed us to translate all of linear algebra idioms into convenient notation.

If translating all of existing idioms is your goal, then this is certainly the only design. But adapting the sloppy (not really convenient) mathematical notation is not a good design guideline. I advise everyone who likes this kind of convenience to use Maple, MatLab, and friends instead!

On Wed, 29 Jun 2005, Jacques Carette wrote:

In fact, type classes in Haskell is a *great* way to do just that!

I agree. I'm also aware of that I mean different objects when I write uniformly '1'. But I know that they are somehow different. I'm also ok with not writing a conversion from say the integer 1 to the fraction 1/1, but I know that there had to be one. We should be aware of such abbreviations before adding them to a programming language.

This is called operator overloading. It is completely harmless because you can tell the two * apart from their type signatures. It is a complete and total waste of time to use two different symbols to mean the same thing.

But the multiplication in A*x already needs multi-parameter type classes. :-) Btw. A function in prefix notation for matrix vector multiplication would also stress the role of a matrix as linear mapping. This is because this function would just map the matrix to the linear mapping it represents. Matrix.mulVec :: Matrix -> (Vector -> Vector)

Abstract Algebra is the branch of mathematics where you want to abstract out the *unimportant* details.

But typically they don't say 'identical' if they mean 'isomorphic'.

Note that the one giant try at a statically typed mathematics system (Axiom w/ programming language Aldor) never caught on,

You mean I should have a closer look on it?

...

If translating all of existing idioms is your goal, then this is certainly the only design. But adapting the sloppy (not really convenient) mathematical notation is not a good design guideline.

I should know better than to get into a discussion like this with someone who believes they singlehandedly know better than tens of thousands of mathematicians...

Millions of people believe that it makes no sense to work with infinite dimensional spaces ... The next argument, please. :-) I learned from "History of mathematical notation" by Florian Cajori that defending _the_ mathematical notation is useless because there is not the one notation. There are so much more notations than we know of. Notations were chosen because of typographical issues, because there was no notation for a specific operation before and so on. But there was never a mathematician who designed a consistent notation from scratch, there was no commitee to setup a standard, mathematical notation was always a patchwork from very different sources. It looks like too many cooks spoil the broth.

On Wed, 29 Jun 2005, Henning Thielemann wrote:

On Wed, 29 Jun 2005, Jacques Carette wrote:

...

This is called operator overloading. It is completely harmless because you can tell the two * apart from their type signatures. It is a complete and total waste of time to use two different symbols to mean the same thing.

But the multiplication in A*x already needs multi-parameter type classes. :-)

I'm uncertain about how who want to put the different kinds of multiplication into one method, even with multi-parameter type classes. You need instances (*) :: Matrix -> Matrix -> Matrix (*) :: RowVector -> Matrix -> RowVector (*) :: Matrix -> ColumnVector -> ColumnVector (*) :: RowVector -> ColumnVector -> Scalar (*) :: ColumnVector -> RowVector -> Matrix (*) :: Scalar -> RowVector -> RowVector (*) :: RowVector -> Scalar -> RowVector (*) :: Scalar -> ColumnVector -> ColumnVector (*) :: ColumnVector -> Scalar -> ColumnVector but you have to make sure that it is not possible to write an expression which needs (*) :: Matrix -> RowVector -> RowVector Further you need transpose :: RowVector -> ColumnVector transpose :: ColumnVector -> RowVector transpose :: Matrix -> Matrix and you must forbid, say transpose :: RowVector -> RowVector Not to mention many special (linear) operators which need to be prepared for both column and row vectors, like Fourier transform, component permutation, convolution, which essentially do the same on both types. I assume that you mean with "common linear algebra idioms" transforms like the following one. (x and y are column vectors) (x * y^T)^2 = x * y^T * x * y^T = x * (y^T * x) * y^T | because y^T * x is "scalar" = (y^T * x) * (x * y^T) (Analogously this can be done with bra-ket notation, i.e. , then transposition must be adjoint) Now your concerns are that this becomes too complicated if the distinction between row and column vectors is dropped and instead different operators are used?

...

I should know better than to get into a discussion like this with someone who believes they singlehandedly know better than tens of thousands of mathematicians...

Is the design discussion related to linear algebra for Maple documented somewhere? For Modula-3 such a discussion was recorded on video tape (then got lost :-( ) and was later edited in "Systems programming with Modula-3" by Greg Nelson. It's fun to read and very informative and I haven't seen such a discussion for other languages.

On Thu, 30 Jun 2005, Jacques Carette wrote:

Henning Thielemann wrote:

...

I'm uncertain about how who want to put the different kinds of multiplication into one method, even with multi-parameter type classes. You need instances

(*) :: Matrix -> Matrix -> Matrix (*) :: RowVector -> Matrix -> RowVector

[many other instances removed.]

Definitely not. You could do: Data Orientation = Row | Column Data Vector a = Vector Orientation [a]

In the first mail you wrote "9. There are row vectors and column vectors, and these are different types. You get type errors if you mix them incorrectly." I interpreted that you want to encode the information Row or Column into the type. This sounds reasonable to me because multiplying e.g. a column vector by a matrix is so obviously wrong that it should be detected statically.

On Thu, 30 Jun 2005, Jacques Carette wrote:

Henning Thielemann wrote:

...

...

Data Orientation = Row | Column Data Vector a = Vector Orientation [a]

In the first mail you wrote "9. There are row vectors and column vectors, and these are different types. You get type errors if you mix them incorrectly."

I interpreted that you want to encode the information Row or Column into the type. This sounds reasonable to me because multiplying e.g. a column vector by a matrix is so obviously wrong that it should be detected statically.

Sorry, I should have been more precise, I used too straightforward an encoding from Maple (which is dynamically typed) into Haskell. I should have used the usual type-trickery to encode the Orientation as 2 different types, to have ``types reflect values''.

You mean phantom types? data Vector orient a = Vector [a] I'm excited to add lots of new orientations then! :-]

On Tue, 5 Jul 2005, Keean Schupke wrote:

Henning Thielemann wrote:

...

I'm uncertain about how who want to put the different kinds of multiplication into one method, even with multi-parameter type classes. You need instances

(*) :: Matrix -> Matrix -> Matrix (*) :: RowVector -> Matrix -> RowVector (*) :: Matrix -> ColumnVector -> ColumnVector (*) :: RowVector -> ColumnVector -> Scalar (*) :: ColumnVector -> RowVector -> Matrix (*) :: Scalar -> RowVector -> RowVector (*) :: RowVector -> Scalar -> RowVector (*) :: Scalar -> ColumnVector -> ColumnVector (*) :: ColumnVector -> Scalar -> ColumnVector

but you have to make sure that it is not possible to write an expression which needs (*) :: Matrix -> RowVector -> RowVector

This was my reply to Jacques Carette who suggested to distinguish row and column vectors by their _type_. His revised suggestion was to use phantom types. Then the types changes to (*) :: Matrix -> Matrix -> Matrix (*) :: Vector Row -> Matrix -> Vector Row (*) :: Matrix -> Vector Column -> Vector Column (*) :: Vector Row -> Vector Column -> Scalar ... but the problem remains ...

...

Further you need transpose :: RowVector -> ColumnVector transpose :: ColumnVector -> RowVector transpose :: Matrix -> Matrix and you must forbid, say transpose :: RowVector -> RowVector

Of course if they are all of type Matrix this problem disappears.

All type problems disappear when we switch to exclusive String representation. 8-]

What is the difference between a 1xN matrix and a vector? Please explain...

My objections to making everything a matrix were the objections I sketched for MatLab. The example, again: If you write some common expression like transpose x * a * x then both the human reader and the compiler don't know whether x is a "true" matrix or if it simulates a column or a row vector. It may be that 'x' is a row vector and 'a' a 1x1 matrix, then the expression denotes a square matrix of the size of the vector simulated by 'x'. It may be that 'x' is a column vector and 'a' a square matrix. Certainly in most cases I want the latter one and I want to have a scalar as a result. But if everything is a matrix then I have to check at run-time if the result is a 1x1 matrix and then I have to extract the only element from this matrix. If I omit the 1x1 test mistakes can remain undiscovered. I blame the common notation x^T * A * x for this trouble since the alternative notation can be immediately transcribed into computer language code while retaining strong distinction between a vector and a matrix type. More examples? More theory?

On Tue, 5 Jul 2005, Michael Walter wrote:

On 7/5/05, Henning Thielemann wrote:

...

The example, again: If you write some common expression like

transpose x * a * x

then both the human reader and the compiler don't know whether x is a "true" matrix or if it simulates a column or a row vector.

But since a, by definition (question), is a 1xN matrix, what's the ambiguity exactly?

If you want this definition then you must also interpret any 1x1 matrix as a real. That's what I wanted to show with my example, that's the way MatLab works and why it sucks. Multiplication of reals is commutative, reals are naturally totally ordered and so on, matrices (including 1x1 matrices) don't have these properties. Since it is sensible to work with one dimensional vectors it is also sensible to work with 1x1 matrices. But 1x1 matrices of this kind are certainly different from 1x1 matrices produced by transpose x * a * x. A 1x1 matrix in MatLab can thus mean a scalar, a row 1-vector, a column 1-vector or a 1x1-matrix. (If you accept the differences between these terms.) Alternatively you could convert the expected 1x1-matrix into a real which must be checked at run time. Theoretically vectors are objects which can be scaled and added and matrices represent linear operators on vectors. (Operators including linear operators may build a vector space itself, but that's a different issue.) Why should we convert each vector into the representation of some linear operator before doing linear algebra and why should we convert this representation of a linear operator back to the vector after linear algebra has happened? Would you load an audio signal into a 1xN-matrix or into a N-vector? Loading it into a 1xN-matrix forces you to check dynamically and repeatedly if the matrix has really only row. Btw. I would also load an image into a vector, but a vector with a two-dimensional index set.

David Roundy wrote:

The issue is that Haskell (as far as I understand, and noone has suggested anything to the contrary) doesn't have a sufficiently powerful type system to represent matrices or vectors in a statically typed way. It would be wonderful if we could represent matrix multiplication as

matrix_mul :: Matrix n m -> Matrix m o -> Matrix n o

Actually, Haskell does allow you to do that. But the syntax of the types gets pretty horrendous. I'm sure Oleg will show you how any moment now. :) -- Lennart

On Thu, 7 Jul 2005, Conal Elliott wrote:

Maybe we could solve this problem in a simple and general way by working with a more abstract notion of linear maps, rather than the matrices commonly used to represent linear maps. Instead of "Matrix n m", where n and m are either integers (requiring something like dependent types) or type encodings of integers, use "LMap u v", where u and v are vector spaces, e.g., R^n and R^m. This change would eliminate the need for dependent types or integer encodings.

I would also dislike integer encodings in the matrix type because it is very sensible to have other types of indices. I would just do it how it is solved for Arrays: Two index types for a matrix but the actual matrix size is not coded in the type. E.g. I would like to multiply a matrix of type Matrix (Bool, Char) Double with a vector of type Vector Char Double obtaining a result of type Vector Bool Double.

Also, it's clear that "LMap R R^n" and "LMap R^n R" are different types (for n/=1), so we can distinguish between "row" and "column" vectors without creating more types. Are there other reasons for orienting vectors?

My point was that vectors naturally do _not_ represent linear maps at all, but they are the objects linear maps act on. If I process an audio signal or an image I can consider it well as vector but why should I consider it as linear map?

I can't call these comments a concrete proposal, as it's not clear to me how to define LMap and what operations it supports.

The problem will be that there is no canonical representation for linear maps in general vector spaces. I consider a function of type (Double -> Double) as an approximation to a real function but we don't have a general representation for linear operators on real functions. The remaining problem is that Vector and Linear Map are relative terms. That is, as you point out, a Linear Map can be a Vector with respect to a higher order Linear Map. I think in contrast to real mathematics we have to distinct due to computational restrictions. On the one hand we can work with a vector space type class which only specifies what basic operations (namely add and scale) shall be supported by types that want to be called vectors. This type class has the full flexibility but no canonical representation for linear maps. On the other hand we define finite dimensional vectors (array like) and matrices as representations for linear maps. For the case of higher order linear maps there should be a conversion between Matrix and Vector.

On Thu, 7 Jul 2005, David Roundy wrote:

On Tue, Jul 05, 2005 at 08:17:58PM +0200, Henning Thielemann wrote:

...

The example, again: If you write some common expression like

transpose x * a * x

then both the human reader and the compiler don't know whether x is a "true" matrix or if it simulates a column or a row vector. It may be that 'x' is a row vector and 'a' a 1x1 matrix, then the expression denotes a square matrix of the size of the vector simulated by 'x'. It may be that 'x' is a column vector and 'a' a square matrix. Certainly in most cases I want the latter one and I want to have a scalar as a result. But if everything is a matrix then I have to check at run-time if the result is a 1x1 matrix and then I have to extract the only element from this matrix. If I omit the 1x1 test mistakes can remain undiscovered. I blame the common notation x^T * A * x for this trouble since the alternative notation can be immediately transcribed into computer language code while retaining strong distinction between a vector and a matrix type.

The issue is that Haskell (as far as I understand, and noone has suggested anything to the contrary) doesn't have a sufficiently powerful type system to represent matrices or vectors in a statically typed way. It would be wonderful if we could represent matrix multiplication as

matrix_mul :: Matrix n m -> Matrix m o -> Matrix n o

Even if we had that I would vote for distinction of Vector and Matrix! I wanted to show with my example that vectors and matrices are different enough that it makes sense to give them different types. In my opinion multiplying a matrix with a so-called column vector is a more fundamental bug than multiplying a matrix with a vector of non-compatible size. That is I like static check of matrix-vector combinations and a dynamic check of their size. The latter also because in many application you don't know the vector size at compile time.

You think of vectors as fundamental and matrices as a sort of derived quantity. I'd prefer to think the other way around, with vectors being special cases of matrices, since it's simpler to implement and work with.

It's less type-safe and thus cumbersome to work with. That's my experience with MatLab.

Also note that if you have several vectors x for which you want to compute the dot product with metric A, and if you want to do this efficiently, you'll have to convert your list of vectors into a matrix anyways.

If you bundle some vectors as columns in matrix B and want to compute norms with respect to matrix A writing B^T * A * B you will not only get the norms of the vectors in B but also many mixed scalar products. This example shows to me that matrices are not simply collections of vectors. Btw. we should not mess up the design with decisions on efficiency concerns. If you want efficiency then you can abuse matrices for that but I consider a 'map' over a list of vectors as the cleaner way (and more type safe, because you can't make transposition errors).

On Thu, 7 Jul 2005, Henning Thielemann wrote:

...

...

fft([1,0;0,0])

ans =

1 0 1 0

Also funny:

...

conv([1;1],[1,1])

ans = 1 2 1

...

conv([1;1;1],[1,1])

ans = 1 2 2 1

On Fri, 8 Jul 2005, Keean Schupke wrote:

Okay, this approach is starting to make sense to me... I can see now that Vectors are a different type of object to Matrices. Vectors represent points in N-Space and matrices represent operations on those points

That's what I wanted to express.

(say rotations or translations)... But it seems we can represent translations as adding vectors or matrix operations (although we need to introduce the 'extra' dimension W... and have an extra field in vectors that contains the value '1').

(3D translation)

[x,y,z,1] * [[0,0,0,0],[0,0,0,0],[0,0,0,0],[dx,dy,dz,dw]] = [x+dx,y+dy,z+dz,1+dw]

Do you mean [x,y,z,1] * [[1,0,0,0],[0,1,0,0],[0,0,1,0],[dx,dy,dz,dw+1]] ?

but how is this different from adding vectors? If we allow vector addition then we no longer have the nice separation between values and linear operators, as a value can also be a linear operator (a translation)?

???

On Fri, 8 Jul 2005, Keean Schupke wrote:

Henning Thielemann wrote:

...

Do you mean [x,y,z,1] * [[1,0,0,0],[0,1,0,0],[0,0,1,0],[dx,dy,dz,dw+1]] ?

Erm, yes thats what I meant ... but you obviously got the point.

...

...

but how is this different from adding vectors? If we allow vector addition then we no longer have the nice separation between values and linear operators, as a value can also be a linear operator (a translation)?

???

Well if a vector can be a linear-operator, then surely it _is_ a matrix!

In general a vector need not to be a linear operator. You talked about vector translation, translation is not a linear operator. You gave some process to map the problem to somewhere, where it becomes a linear operator. Other people said that the scalar product with a fixed vector is a linear operator. That's true. Now what is a natural interpretation of a vector as linear operator? The scalar product or the translation? Vectors can be used and abused for many things but an object which can be called a vector (because of its ability of to be added and to be scaled) is not a linear operator itself and does not naturally represent one.

On Fri, 8 Jul 2005, Keean Schupke wrote:

So the linear operator is translation (ie: + v)... effectively 'plus' could be viewed as a function which takes a vector and returns a matrix (operator)

(+) :: Vector -> Matrix

Since a matrix _is_ not a linear map but only its representation, this would not make sense. As I said (v+) is not a linear map thus there is no matrix which represents it. A linear map f must fulfill f 0 == 0 But since v+0 == v the function (v+) is only a linear map if 'v' is zero. I can't see how to fit in your vector extension by the 1-component.

On Fri, 8 Jul 2005, Keean Schupke wrote:

Henning Thielemann wrote:

...

On Fri, 8 Jul 2005, Keean Schupke wrote:

...

So the linear operator is translation (ie: + v)... effectively 'plus' could be viewed as a function which takes a vector and returns a matrix (operator)

(+) :: Vector -> Matrix

Since a matrix _is_ not a linear map but only its representation, this would not make sense. As I said (v+) is not a linear map thus there is no matrix which represents it. A linear map f must fulfill f 0 == 0

But since v+0 == v the function (v+) is only a linear map if 'v' is zero.

I can't see how to fit in your vector extension by the 1-component.

Eh?

Today: +----------------+ | The Linear Map | +----------------+ If F is a field, V a set of vectors, and (F,V,+,*) a vectorspace then a function f of V -> V is linear if: for all x and y from V f (x+y) == f x + f y for all k from F and x from V f (k*x) == k * f x Lemma If f is a linear map, then f 0 = 0. Proof: For any v from V it is v+(-1)*v = 0. f 0 = f (v+(-1)*v) = f v + f ((-1)*v) = f v + (-1) * f v = 0 Theorem If (v+) is a linear map then v == 0. Proof: (v+0) == 0 (conclusion from the lemma) => v == 0 8-]

(or the second by the first - the two are isomorphic)... A translation can be represented by the matrix:

1 0 0 0 0 1 0 0 0 0 1 0 dx dy dz 1

So the result of "v+" is this matrix.

But the vectors I can add to v have one component less than necessary for multiplication with this matrix. Indeed you can map all v's with three components to an affine sub-space of the 4-dimensional space, namely to those vectors with a 1 in the last component. But you are no longer working with three dimensional vectors but with four-dimensional ones. Again: isomorphy is not identity.

On Friday 08 July 2005 17:46, Henning Thielemann wrote:

Vectors can be used and abused for many things but an object which can be called a vector (because of its ability of to be added and to be scaled) is not a linear operator itself and does not naturally represent one.

At least for finite dimensional spaces (and these are the only ones under consideration here, right?) scalar multiplication is a very nice and natural way to view a vector as a linear operation (into the scalar field). I know, in linear algebra the told us all that this corespondence is not a 'natural' or 'canonic' one because it depends on a chosen basis. Well, well. For practical purposes of programming, we always use the 'canonic base' right? So if the base is canonic, then so is the correspondence between vector space and its dual. On a different not, one could argue that a 1xn matrix M is indeed a vector of dimension n, an then M' = M^T is a nx1 matrix, that is also a vector of dimension n, but this is the same vector as the non-transposed version. Now, the two things (M and M') are the same, if viewed as a vector, but not the same if viewed as a matrix. Can we express this in Haskell? Ben

On Thu, 7 Jul 2005, David Roundy wrote:

On the other hand, this is sort of a silly debate, since the API *I* want is a subset of the API *you* want.

The API you want is certainly not just a subset of what I want. E.g. the singular value decomposition in your favorite API will probably return a 1xN matrix or a MxN diagonal matrix containing the singular values, while with my favorite API it will return a vector or a list. Right?

On Fri, 8 Jul 2005, David Roundy wrote:

I don't particularly care what API you use for svd, since it's trivial to convert from one API to the other. It's matrix arithmetic I care about, since that's the complicated part of the API.

Of course I want to use the results of more complicated routines with basic matrix arithmetic and I like to reduce the number of conversions. The other reason for the debate is that if you don't like the extra vector type then you will not use it. If I want to apply a routine of you to my, say, audio data I will have to decide whether I shall store it as column or as row vector/matrix although this decision seems to me rather irrelevant and annoying.

On the other hand, the most natural return value for svd would be a diagonal matrix, since that is what the objects are, right?

Hm, since SVD means Singular Value Decomposition I like to have the singular values as they are. I don't want to search them in a sparse matrix.

svd returns three matrices, which when multiplied together give the original matrix ...

This would be a nice property, though. I could do it by converting the singular value list to a diagonal matrix. So I need a conversion, hm.

or at least that's how I think of it. But I'll grant that a diagonal matrix isn't the most convenient representation, and certain is far from the most efficient, unless we introduce a diagonal matrix constructor (which would certainly be nice).

I would not like to obtain a value of general Matrix type since it is statically guaranted that it is always diagonal.

I guess you'd prefer that svd returns a list of doubles and two lists of vectors? Or a list of triplets of a double and two vectors?

Either there is a function to scale the columns of a matrix separately, then two matrices and a list/vector of doubles are ok. Or there is a function to multiply two vectors to obtain a matrix with rank 1 (Outer product or tensor product might be the right name. Btw. the (<>) operator has the problem, that it is always interpreted as scalar product. But it could also mean this kind of multiplication.) Then adding such products of left and right singular vectors scaled by the corresponding singular value let me reconstruct the matrix. The triplet approach has the advantage that it is statically sure that the number of singular values and singular vectors match. The matrix approach has the advantage that it is statically guaranteed that the dimension of all singular vectors match. With the (orthogonal) singular vector matrices we map vectors from one basis to the representation in another basis. So these matrices represents naturally linear maps and it makes sense to use them.

Hello! Thank you very much for all your suggestions. A new version of the library can be found at: http://dis.um.es/~alberto/hmatrix/matrix.html Here are the main changes: - Vector and Matrix are now different types with different functions operating on them. They cannot be interchanged and we must use explicit conversion functions. Things like (v::Vector) + (m::Matrix), or svd (v::Vector) are statically rejected. You must use functions like matrix_x_matrix, matrix_x_vector, dot :: Vector->Vector->Double, etc. On the other hand, there are functions to build a matrix from a list of vectors, to concatenate the columns of a matrix as a vector, etc. If I am not mistaken, now the ambiguities explained by Henning do not occur. - Only Eq, Show, and Read instances are now defined for Vector and Matrix. No numeric instances are defined, although it is very easy to do it using the functions provided. I have included a tentative "Interface" module including some hopefully user friendly operators. But it can be replaced or deleted without any problem... :) If this approach makes sense, we can proceed to include more standard linear algebra functions based on GSL, LAPACK, or on any other available resource. We can use a darcs repository to allow contributions from all interested people. And a few comments: David Roundy wrote:

It would be wonderful if we could represent matrix multiplication as

matrix_mul :: Matrix n m -> Matrix m o -> Matrix n o

(where n, m and o are dimensions) but we can't do that. So we're stuck with dynamic size checking for all matrix and vector operations that involve two or more operands.

Lennart Augustsson wrote:

Actually, Haskell does allow you to do that. But the syntax of the types gets pretty horrendous.

We must also "add" types. To build a matrix from blocks we need something like: matrix_block_row :: Matrix n m -> Matrix n o -> Matrix n (o+m) I understand from some papers and messages to the list that it can be done, but then we may also want lists of matrices with different dimensions: matrix_block :: [[Matrix ? ?]] -> Matrix The compile-time consistency check is not trivial. On the other hand, many useful functions have a type which depends on a value: vector :: [Double] -> Vector n We should statically reject this: vector [1,2,3] `vector_add` vector [1,2] (or perhaps we can imitate the zip behavior? :) The sizes of some vectors and matrices depend on run time computations, even without taking into account IO: nullspace :: Matrix m n -> Matrix k n k depends on the particular matrix. We can improve this to: nullspace :: Matrix m n -> [Vector n] but if we need the nullspace in another matrix operation we must build a matrix from a list whose size is not known at compile time... Static size checking is an advanced topic! Alberto

I' sorry, our web server is temporarily down :-( On Thursday 07 July 2005 20:37, Alberto Ruiz wrote:

Hello! Thank you very much for all your suggestions. A new version of the library can be found at:

http://dis.um.es/~alberto/hmatrix/matrix.html

On Fri, 8 Jul 2005, Alberto Ruiz wrote:

The server is working again.

On Thursday 07 July 2005 20:58, Alberto Ruiz wrote:

...

I' sorry, our web server is temporarily down :-(

...

http://dis.um.es/~alberto/hmatrix/matrix.html

I would remove the 'matrix' portions of the function names, since this information is already contained in the module name. E.g. after importing LinearAlgebra.Matrix as Matrix, Matrix.matrix look strange, but Matrix.fromList says everything. I suggest matrix -> fromList matrixElements -> toList displayMatrix -> display readMatrix -> fromFile :: FilePath -> IO Matrix writeMatrix -> toFile :: FilePath -> Matrix -> IO () blockMatrix -> fromBlocks or fromMatrices subMatrix -> clip or cut or take matrixScale -> scale matrixOffset -> offset matrixAdd -> add matrixSub -> sub matrixMul -> mul (element-wise multiplication? then better zipMul or elementwiseMul, elMul what know I :-) matrixDiv -> div matrixSigns -> signs matrixAbs -> abs matrix_x_matrix ... hm difficult matrix_x_vector vector_x_matrix msin :: Matrix -> Matrix mcos :: Matrix -> Matrix mtan :: Matrix -> Matrix masin :: Matrix -> Matrix macos :: Matrix -> Matrix matan :: Matrix -> Matrix matan2 :: Matrix -> Matrix -> Matrix msinh :: Matrix -> Matrix mcosh :: Matrix -> Matrix mtanh :: Matrix -> Matrix masinh :: Matrix -> Matrix macosh :: Matrix -> Matrix matanh :: Matrix -> Matrix mexp :: Matrix -> Matrix mlog :: Matrix -> Matrix If there are no efficiency concerns, I would drop element-wise operations and prefer a matrix-map and a matrix-zipWith. If these operations shall remain I would somehow point to their element-wise operation in the name.

Hello! I have a few doubts concerning the LinearAlgebra library... On Friday 08 July 2005 11:29, Henning Thielemann wrote:

I would remove the 'matrix' portions of the function names, since this information is already contained in the module name. E.g. after importing LinearAlgebra.Matrix as Matrix, Matrix.matrix look strange, but Matrix.fromList says everything.

I have changed the function names as suggested. This new style is clearly better, allowing Vector.add, Matrix.add, Vector.Complex.add, Matrix.Complex.add, etc.

If the Matrix type would be parametrised (...)

Now we can have Vector.T a and Matrix.T a for any storable a (although at this point most functions are only defined for Double). For example: eig :: Matrix.T Double -> ([Double], Matrix.T Double)

I would also prefer a vector of complex numbers for the FFT function instead of a nx2 matrix.

This comes from the old version in which vectors and matrices were the same type :-) Now we have: fft1D :: Vector.T (Complex Double) -> Vector.T (Complex Double) ifft1D :: Vector.T (Complex Double) -> Vector.T (Complex Double) fft2D :: Matrix (Complex Double) -> Matrix (Complex Double) etc. But now I have two problems: 1) If I define foo :: Vector.T a -> Matrix.T a Haddock (version 0.6) shows just this: foo :: T a -> T a I don't know how to tell haddock that I want the full names in the signatures. 2) When setting up the "main" module which reexports all the modules exposed to the user I cannot write this: --------------------------------- module LinearAlgebra ( module Vector, module Matrix, module Matrix.Complex, module LinearAlgebra.Algorithms --etc. ) where import LinearAlgebra.Vector as Vector import LinearAlgebra.Matrix as Matrix import LinearAlgebra.Matrix.Complex as Matrix.Complex import LinearAlgebra.Algorithms --etc. ------------------------------------- since we get lots of conflicting exports for T, add, scale, etc. I we use qualified imports: ----------------------------------------------------------- module LinearAlgebra ( module Vector, module Matrix, module Matrix.Complex, module LinearAlgebra.Algorithms --etc. ) where import qualified LinearAlgebra.Vector as Vector import qualified LinearAlgebra.Matrix as Matrix import qualified LinearAlgebra.Matrix.Complex as Matrix.Complex import LinearAlgebra.Algorithms ----------------------------------------------------------- then all the names in Vector, Matrix and Matrix.Complex must be always qualified, which is not very practical. I am afraid that I am doing something wrong... Another possibility is using type classes, since there is a lot of types addmiting add, sub, (ew)mul, scale, etc. (e.g., Vector Int, Vector Double, Matrix (Complex Float), etc.) We could also define type classes for higher level linear algebra algorithms working with different base types in the Matrices. Any suggestion? And some final comments.

If there are no efficiency concerns, I would drop element-wise operations and prefer a matrix-map and a matrix-zipWith. If these operations shall remain I would somehow point to their element-wise operation in the name.

There is about 5x speed gain if we map in the C side. The "optimized" floating map functions could be moved to a separate module.

If the Matrix type would be parametrised then Matrix.fromBlocks could use a more natural indexing.

Matrix.fromBlocks :: [[Matrix a b]] -> Matrix (Int,a) (Int,b)

The Int of the index pair would store the block number and the type a index would store the index within the block. Hm, but it would not be possible to store the sizes of the sub-matrices. It would be possible only if the sub-matrices have the same size. Maybe it is better to allow matrices of matrices and to be able to apply a matrix-matrix multiplication which in turn employs a matrix-matrix multiplication of its elements. But the matrix decompositions will fail on this structure - but possibly those which fail aren't appropriate for block matrices.

We can have a Matrix (Matrix Double) if we define the instance Storable (Matrix a). Although a "block matrix" is perhaps better represented at a higher level (storing only one copy of repeated blocks, etc.).

It would be nice to have functions which don't simply write the formatted vectors and matrices to stdout but which return them as strings. E.g.

Matrix.format :: Int -> Matrix.T -> String Vector.format :: Int -> Vector.T -> String

(or Matrix.toString)

The results can be reused in other contexts and the 'display' routines can use their results.

Done. --Alberto

On Wed, Jul 13, 2005 at 06:13:48PM +0200, Alberto Ruiz wrote:

I have changed the function names as suggested. This new style is clearly better, allowing Vector.add, Matrix.add, Vector.Complex.add, Matrix.Complex.add, etc. ... Now we can have Vector.T a and Matrix.T a for any storable a (although at this point most functions are only defined for Double). For example:

I would like to bristle mildly against the style of using Vector.T to represent the vector type. The reasons are 1) it is cryptic to those not used to the convention; 2) enshrining one-type-per-module in the naming convention is not IMO justified, and may prove limiting; 3) it doesn't work out well when you import a module that reexports Vector. I would say leave it Vector.Vector. Then the user may import the module unqualified, or qualified with an abbreviation like V, define his own synonym ("type Vector = Vector.Vector"), or just put up with the light annoyance of writing "Vector.Vector".

1) If I define

foo :: Vector.T a -> Matrix.T a

Haddock (version 0.6) shows just this:

foo :: T a -> T a

I think this is evidence that trying to treat the "primary type" as a special case in naming creates more annoyances than it eliminates. Andrew

On Fri, Jul 15, 2005 at 10:48:04AM +0200, Henning Thielemann wrote:

On Thu, 14 Jul 2005, Andrew Pimlott wrote:

...

2) enshrining one-type-per-module in the naming convention is not IMO justified, and may prove limiting;

Other languages like Modula-3 and Oberon do it with great success. The limit in Haskell is that most compilers don't conform to the Haskell 98 report which allows mutually recursive modules. But I think the compilers should allow them instead of forcing users to put many type and class definitions into one module.

Even given an ideal implementation (I would add that it should allow multiple modules in one file), I don't find one type or class per module preferable. I think it's usually a false division. For one, there will be functions that operate on multiple types, so where do you put those? For another, the user will probably have to import several modules, and remember which which module contains a function that operates on multiple types. That said, my point is not that this style of choosing module boundaries is always wrong, or that you can't answer the above objections in many cases. It's that you shouldn't build a general naming convention around it, especially one that some will find confusing or awkward, and when it has only small practical advantages (and some disadvantages). But I'm not going to push this any further: If people still like it, it's not that big a deal to me. :-)

Vector.T is just the consequence of removing type-related prefixes from all identifiers.

I interpreted the style as removing anything that is redundant in the context of the module. In math, when you're talking about vector spaces, you say "add" and "multiply", but you don't say "thing" or "one"--you still say "vector".

If you stick to qualified naming you only need to import the module once. Otherwise you must import a module twice.

import LinearAlgebra.Vector (Vector) import qualified LinearAlgebra.Vector as V

I don't think I understand what you're saying with this example. If you export the type as "Vector", the user could "import LinearAlgebra.Vector as Vector" and then refer to both "Vector" and "Vector.add". Or "import LinearAlgebra.Vector as V" and then refer to "Vector" (or "V.Vector") and "V.Add". Andrew

On Fri, 15 Jul 2005, Andrew Pimlott wrote:

Even given an ideal implementation (I would add that it should allow multiple modules in one file),

Why?

I don't find one type or class per module preferable. I think it's usually a false division.

It helped me to decide for divisions early. What do you use as division criterion?

For one, there will be functions that operate on multiple types, so where do you put those?

I try to find out which is the more complex one, which one depends on the other. In this example I consider a matrix as a representation of a linear map (see the long preceding discussion about how matrices can be interpreted :-), where the linear maps operate on vectors. Thus matrices depend on vectors but not vice versa, so I put (basic) matrix-vector operations to the matrix module. More complicated functions can be moved in a separate module.

For another, the user will probably have to import several modules, and remember which which module contains a function that operates on multiple types.

Yes, that's the way modularization work. That's not a problem with the T name but general to modularization, isn't it?

When it has only small practical advantages (and some disadvantages).

The big advantage is that you can keep the interfaces of all modules consistent. Ideally you can switch an 'import qualified ... as' statement to a different module which serves the same interface with different implementation. For instance you could switch between matrix implementations based on different back-ends. Though, this is probably better solved with type classes.

I have mixed feelings about Module.T, but consider it useful. basically, I use it only when the module naturally only exports a single type which is the purpose of the module (modules which incidentally only export a single type but have a different focus shouldn't use the type synonym) and the modules functions are named in such a way that it is intended to be used in a qualified way. (like the new Data.Map and Data.Set). However, I don't try to force any modules or types into this model if they don't naturally fall into them. I also make sure that the T is a type synonym for the actual name. as in module Vector where data Vector = ... type T = Vector This ensures that tools that strip the module name still give reasonable names to things since most expand type synonyms away. So, in conclusion, I think Vector and Matrix do naturally fall into these conditions so should have T but as type synonyms and not the base name. John -- John Meacham - ⑆repetae.net⑆john⑈

Hello Alberto, Wednesday, July 13, 2005, 8:13:48 PM, you wrote:

...

If there are no efficiency concerns, I would drop element-wise operations and prefer a matrix-map and a matrix-zipWith. If these operations shall remain I would somehow point to their element-wise operation in the name.

AR> There is about 5x speed gain if we map in the C side. The "optimized" floating AR> map functions could be moved to a separate module. GHC also have a RULES pragma which can be used to automatically convert, for example, "mmap (*)" to "multipleElementWise". below is examples of using this pragma in the standard GHC modules: {-# RULES "foldr/id" foldr (:) [] = \x->x "foldr/single" forall k z x. foldr k z [x] = k x z "foldr/nil" forall k z. foldr k z [] = z #-} -- Best regards, Bulat mailto:bulatz@HotPOP.com

Hello Bulat, thanks a lot for your message, the RULES pragma is just what we need! However, in some initial experiments I have observed some strange behavior. For instance, in the following program: ------------------------------------------ {-# OPTIONS_GHC -fglasgow-exts #-} apply :: (Int -> Int) -> Int -> Int apply f n = f n sqr :: Int -> Int sqr n = n * n optimized_sqr :: Int -> Int optimized_sqr n = n*n+1 -- to check that the rule works :-) {-# RULES "apply/sqr" apply sqr = optimized_sqr #-} main = do --print $ apply sqr 3 print $ apply sqr 5 ----------------------------------------- The rule is not applied. 1 RuleFired 1 *# if we uncomment the first line in the main function main = do print $ apply sqr 3 print $ apply sqr 5 then the rule is correctly applied: 6 RuleFired 2 *# 2 +# 2 apply/sqr Solution: include at the beginning of the file module Main where and then the rule works in both cases. I have a similar problem in the LinearAlgebra library but there, curiously, the rule only works if it is applied once: module Main where import (...) (...) main = do (...) print $ Vector.map cos v --print $ Vector.map cos v ==================== Grand total simplifier statistics ==================== Total ticks: 2584 461 PreInlineUnconditionally 230 PostInlineUnconditionally 387 UnfoldingDone 91 RuleFired 2 *# 16 *## 5 +## 8 ++ 5 -## 2 SPEC $fLinearArray1 2 SPEC $fLinearArray2 1 SPEC $fNumComplex 1 SPEC $fShowComplex 1 Vector.map/cos <--------------- OK 20 int2Double# 4 map 4 mapList 2 plusDouble 0.0 x 4 plusDouble x 0.0 2 timesDouble x 0.0 2 timesDouble x 1.0 3 unpack 3 unpack-list 2 zipWith 2 zipWithList 47 LetFloatFromLet 9 EtaReduction 1136 BetaReduction 6 CaseOfCase 217 KnownBranch 14 SimplifierDone But: (...) main = do (...) print $ Vector.map cos v print $ Vector.map cos v ==================== Grand total simplifier statistics ==================== Total ticks: 2664 470 PreInlineUnconditionally 240 PostInlineUnconditionally 402 UnfoldingDone 90 RuleFired 2 *# 16 *## 5 +## 8 ++ 5 -## 2 SPEC $fLinearArray1 2 SPEC $fLinearArray2 1 SPEC $fNumComplex 1 SPEC $fShowComplex 20 int2Double# 4 map 4 mapList 2 plusDouble 0.0 x 4 plusDouble x 0.0 2 timesDouble x 0.0 2 timesDouble x 1.0 3 unpack 3 unpack-list 2 zipWith 2 zipWithList 49 LetFloatFromLet 9 EtaReduction 1181 BetaReduction 5 CaseOfCase 218 KnownBranch 17 SimplifierDone I have tried several ideas, without any luck. Alberto On Monday 18 July 2005 10:14, Bulat Ziganshin wrote:

Hello Alberto,

Wednesday, July 13, 2005, 8:13:48 PM, you wrote:

...

...

If there are no efficiency concerns, I would drop element-wise operations and prefer a matrix-map and a matrix-zipWith. If these operations shall remain I would somehow point to their element-wise operation in the name.

AR> There is about 5x speed gain if we map in the C side. The "optimized" floating AR> map functions could be moved to a separate module.

GHC also have a RULES pragma which can be used to automatically convert, for example, "mmap (*)" to "multipleElementWise". below is examples of using this pragma in the standard GHC modules:

{-# RULES "foldr/id" foldr (:) [] = \x->x "foldr/single" forall k z x. foldr k z [x] = k x z "foldr/nil" forall k z. foldr k z [] = z #-}

Hi Alberto, I'm sorry if this has been discussed before... I'm reading your paper, at one point it says (re. the PCA example): "Octave achieves the same result, slightly faster. (In this experiment we have not used optimized BLAS libraries which can improve efficiency of the GSL)" That seems to imply that there is a way to use optimized BLAS libraries? How can I do that? Also, in my experiments (with matrix inversion) it seems, subjectively, that Octave is about 5 or so times faster for operations on large matrices. Presumably you've tested this as well, do you have any comparison results? Thanks, Frederik -- http://ofb.net/~frederik/

On Thursday 16 March 2006 18:13, Frederik Eaton wrote:

Also, in my experiments (with matrix inversion) it seems, subjectively, that Octave is about 5 or so times faster for operations on large matrices. Presumably you've tested this as well, do you have any comparison results?

Frederik, you are right, in my machine Octave is at least 3x faster than gsl (v1.5). Too much, specially in a simple matrix multiplication, and I don't know why. See below the times measured in the C side for the PCA example. I will look into this... Alberto gsl 1.5 GHC> main Loading package haskell98-1.0 ... linking ... done. GSL Wrapper submatrix: 0 s GSL Wrapper constant: 0 s GSL Wrapper multiplyR (gsl_blas_dgemm): 0 s GSL Wrapper vector_scale: 0 s GSL Wrapper vector_scale: 1 s GSL Wrapper gsl_vector_add: 0 s GSL Wrapper trans: 0 s GSL Wrapper multiplyR (gsl_blas_dgemm): 36 s <---- (784x5000) x (5000 x 784) GSL Wrapper vector_scale: 0 s GSL Wrapper trans: 0 s GSL Wrapper vector_scale: 0 s GSL Wrapper gsl_vector_add: 0 s GSL Wrapper toScalar: 0 s GSL Wrapper eigensystem: 11 s <---------- (784x784) [337829.65625394206,253504.7867502207, etc. (97.95 secs, 13096315340 bytes) (includes file loading, to be optimized) GNU Octave, version 2.1.57 (i686-pc-linux-gnu). octave:6> testmnist x - repmat(mean(x),rows(x),1) 0.46144 (xc'*xc)/rows(x) 11.623 <-------------- eig 4.3873 <-------------- 3.3776e+05 2.5345e+05 2.0975e+05 etc. ----------octave------------------ load mnist.txt x = mnist(:,1:784); d = mnist(:,785); t0=time(); xc = x - repmat(mean(x),rows(x),1); disp("x - repmat(mean(x),rows(x),1)"); disp(time()-t0) t0=time(); mc = (xc'*xc)/rows(x); disp("(xc'*xc)/rows(x)"); disp(time()-t0) t0=time(); [v,l]=eig(mc); disp("eig"); disp(time()-t0) disp(flipud(diag(l))(1:10)); ------------- GSL + Haskell --------- module Main where import GSL import GSL.Util mean x = sumCols x ./. fromIntegral (rows x) center x = x |-| fromRows (replicate (rows x) (mean x)) cov x = (trans xc <> xc) ./. n where n = fromIntegral (rows x -1) xc = center x m ./. r = m <> (1/r ::Double) test :: M -> [Double] test x = take 10 (toList lambdas) where mc = cov x (lambdas, _) = eig mc main = do m <- fromFile "mnist.txt" let x = subMatrix 0 (rows m-1) 0 (cols m -2) m print (test x)

it appears possible to define matrix operations like LU-decomposition, SVD, inverse etc. in a recursive fashion (see transpose in the prelude). is this possible? has anybody code in haskell which does this (not asking for high performance here, more instructive value). thank you! andrew

On Wed, 19 Apr 2006, Andrew U. Frank wrote:

it appears possible to define matrix operations like LU-decomposition, SVD, inverse etc. in a recursive fashion (see transpose in the prelude).

is this possible? has anybody code in haskell which does this (not asking for high performance here, more instructive value).

There is some code by Jan Skibinski called 'indexless linear algebra': http://www.haskell.org/haskellwiki/Libraries_and_tools/Mathematics

On Thu, 7 Jul 2005, Alberto Ruiz wrote:

Hello! Thank you very much for all your suggestions. A new version of the library can be found at:

http://dis.um.es/~alberto/hmatrix/matrix.html

If the Matrix type would be parametrised then Matrix.fromBlocks could use a more natural indexing. Matrix.fromBlocks :: [[Matrix a b]] -> Matrix (Int,a) (Int,b) The Int of the index pair would store the block number and the type a index would store the index within the block. Hm, but it would not be possible to store the sizes of the sub-matrices. It would be possible only if the sub-matrices have the same size. Maybe it is better to allow matrices of matrices and to be able to apply a matrix-matrix multiplication which in turn employs a matrix-matrix multiplication of its elements. But the matrix decompositions will fail on this structure - but possibly those which fail aren't appropriate for block matrices.

On Fri, 8 Jul 2005, Keean Schupke wrote:

Henning Thielemann wrote:

...

My objections to making everything a matrix were the objections I sketched for MatLab.

The example, again: If you write some common expression like

transpose x * a * x

Which just goes to show why haskell limits the '*' operator to multiplying the same types. Keep this to Matrix-times-Matrix operators.

Btw. the interface file by Alberto gives you uniform usage of an operator symbol (<>) while retaining full type safety by functional dependencies. http://dis.um.es/~alberto/hmatrix/doc/LinearAlgebra.Interface.html This is a good solution, I think.

I feel using separate types for vectors and scalars just overcomplicates things...

I'm excited if your code gets swamped by conversions between Double and Matrix then. I really plead for representing everything with strings, this is the most simple and most flexible solution! :-]

On Wed, 29 Jun 2005, Jacques Carette wrote:

<sarcasm>Next thing you know, you'll want a different 'application' symbol for every arity of function, because they are ``different''. </sarcasm>

Btw. there is less sarcasm in it as may you think. There was already a proposal to extend function application: http://www.haskell.org/pipermail/haskell/2002-October/010629.html and guess - I would be very unhappy about such an extension because it mixes the representation of a map with the map itself.

Jacques Carette writes:

Henning Thielemann wrote:

...

I don't see the problem. There are three very different kinds of multiplication, they should also have their own signs: Scalar product, matrix-vector multiplication, matrix-matrix multiplication.

You see 3 concepts, I see one: multiplication. Abstract Algebra is the branch of mathematics where you want to abstract out the *unimportant* details. Much progress was made in the late 1800's when this was discovered by mathematicians ;-).

One of the things I appreciate and I hate simultaneously in your postings is that you are so categorical. This time you will *not* convince me that there is "one concept: multipli- cation", moreover "abstracted over unimportant details". If matrices represent operators, their multiplication is a *group* operation, the op. composition. Acting of a matrix on a vector is not. "Multiplication" of two vectors giving a scalar (their contration) is yet another beast. I believe that some progress has been done in math, when people discovered that mixing-up things is not necessarily a good thing, and different entities should be treated differently. Jerzy Karczmarczuk

On Wed, Jun 29, 2005 at 10:23:23PM +0200, Henning Thielemann wrote:

More specific: You give two different things the same name. You write A*B and you mean a matrix multiplication. Matrix multiplication means finding a representation of the composition of the operators represented by A and B. But you also write A*x and you mean the matrix-vector multiplication. This corresponds to the application of the operator represented by A to the vector x. You see: Two different things, but the same sign (*). Why? You like this ambiguity because of its conciseness. You are used to it. What else? But you have to introduce an orientation of vectors, thus you discriminate two things (row and column vectors) which are essentially the same! What is the overall benefit?

I just wanted to insert a few thoughts on the subject of row and column vectors. (And I'm not going to address most of the content of this discussion.) If we support matrix-matrix multiplication, we already automatically support matrix-column-vector and row-vector-matrix multiplication, whether or not we actually intend to, unless you want to forbid the use of 1xn or nx1 matrices. So (provided we *do* want to support matrix-matrix multiplication, and *I* certainly would like that) there is no question whether we'll have octavish/matlabish functionality in terms of row and column vectors--we already have this behavior with a single multiplication representing a single operation. There is indeed a special case where one or more of the dimensions is 1, but it's just that, a special case, not a separate function. If you want to introduce a more general set of tensor datatypes, you could do that, but matrix arithmetic is all *I* really need. I agree that when you add other data types you'll need other operators (or will need to do something more tricky), but starting with the simple case seems like a good idea, especially since the simple case is the one for which there exist fast algorithms (i.e. level 3 BLAS). There are good reasons when if I've got a set of vectors that I'd like to multiply by a matrix that I should group those vectors into a separate matrix--it's far faster (I'm not sure how much, but usually around 10 times faster, sometimes more). In fact, I'd think it's more likely that a general tensor library will be implemented using matrix operations than the other way around, since it's the matrix-matrix operations that have optimized implementations. In short, especially since the folks doing the work (not me) seem to want plain old octave-style matrix operations, it makes sense to actually do that. *Then* someone can implement an ultra-uber-tensor library on top of that, if they like. And I would be interested in a nice tensor library... it's just that matrices need to be the starting point, since they're where most of the interesting algorithms are (that is, the ones that are interesting to me, such as diagonalization, svd, matrix multiplication, etc). -- David Roundy http://www.darcs.net

On Thu, Jun 30, 2005 at 02:20:16PM +0200, Henning Thielemann wrote:

On Thu, 30 Jun 2005, David Roundy wrote:

...

If we support matrix-matrix multiplication, we already automatically support matrix-column-vector and row-vector-matrix multiplication, whether or not we actually intend to, unless you want to forbid the use of 1xn or nx1 matrices. So (provided we *do* want to support matrix-matrix multiplication, and *I* certainly would like that) there is no question whether we'll have octavish/matlabish functionality in terms of row and column vectors--we already have this behavior with a single multiplication representing a single operation.

Of course, I only wanted a separate vector type (which also means a separate matrix-vector multiplication) and I argued against a further discrimination of row and column vectors.

...

If you want to introduce a more general set of tensor datatypes,

Did someone requested tensor support? At least not me. I used the tensor example to show that MatLab throws them all together with matrices and vectors and I wanted to give an idea how to make it better by separating them.

I disagree. Vectors *are* tensors... and once you've added rank-1 tensors, you'll also be wanting to add support for rank-0 tensors. These are nice objects (and would be nice to support eventually), but they add considerable complexity. There's only one matrix-matrix multiplication, but there are two vector-vector multiplications, two vector-matrix multiplications. All four of these multiplications can be expressed in terms of the single matrix-matrix multiplication. I agree that if we introduce a vector type there's no reason to introduce separate row and column vectors, but I think we can get by quite well for many purposes without introducing the vector type, which adds so much complexity. I guess I've left out the ".*" pointwise multiply, probably because I'm a physicist and am tensor-biased, and it's not a tensor operation in the physics sense (it doesn't behave right under unitary transformations of its arguments).

...

In short, especially since the folks doing the work (not me) seem to want plain old octave-style matrix operations, it makes sense to actually do that. *Then* someone can implement an ultra-uber-tensor library on top of that, if they like. And I would be interested in a nice tensor library... it's just that matrices need to be the starting point,

Matrices _and_ vectors! Because matrices represent operators on vectors and it is certainly not sensible to support only the operators but not the objects they act on ... Adding a vector type by a library that is build on top of a matrix library seems to me like making the first step after the second one.

No, matrices operate on matrices and return matrices. This is the wonderful thing about matrix arithmetic, why it's unique, and why I'd like to have a library that supports matrix arithemetic. Tensors are also nice, but are much more complicated, since most tensor multiplications change the rank of the tensor--and that complication is precisely what you want. You're taking an interperetation of what you mean by the math and wanting to encode it in the type system. Matrix arithmetic is entirely self-contained, and I'd like to have *that* reflected in the type system. -- David Roundy http://www.darcs.net

David Roundy wrote:

In short, especially since the folks doing the work (not me) seem to want plain old octave-style matrix operations, it makes sense to actually do that. *Then* someone can implement an ultra-uber-tensor library on top of that, if they like. And I would be interested in a nice tensor library... it's just that matrices need to be the starting point, since they're where most of the interesting algorithms are (that is, the ones that are interesting to me, such as diagonalization, svd, matrix multiplication, etc).

This is a really good idea. I would like a Matrix library soon, not in 6 years time. Slice it up into managable pieces and keep it simple! Keean.

On Tue, Jul 05, 2005 at 07:49:56PM +0200, Henning Thielemann wrote:

On Tue, 5 Jul 2005, Keean Schupke wrote:

...

David Roundy wrote:

...

In short, especially since the folks doing the work (not me) seem to want plain old octave-style matrix operations, it makes sense to actually do that. *Then* someone can implement an ultra-uber-tensor library on top of that, if they like. And I would be interested in a nice tensor library... it's just that matrices need to be the starting point, since they're where most of the interesting algorithms are (that is, the ones that are interesting to me, such as diagonalization, svd, matrix multiplication, etc).

This is a really good idea. I would like a Matrix library soon, not in 6 years time. Slice it up into managable pieces and keep it simple!

As I said, _that_ already exists: MatLab, HaskellDSP (with some simple new definitions for infix operators) ...

Except that matlab isn't a Matrix library--it's a horribly language--and HaskellDSP doesn't seem to implement linear algebra and its interface is such that it is probably very difficult to implement efficiently (since efficient implementation means calling lapack). In my opinion a decent matrix API needs to have an opaque data type. In other words, using HaskellDSP to create a decent haskell Matrix library would be about as much work as using lapack to do so. -- David Roundy http://www.darcs.net