On Mon, 26 Apr 2010 03:46:47 -0300, Gareth Charnock <gareth.tpc@gmail.com> wrote:

> Thanks. To quickly answer this:
>
>  > - The version I'm seeing on github doesn't seem to have all the features
>  > you're referencing (i.e. v*v). Why are some ops limited to */ and other +-?
>
> It was quite late (3am) when I typed up that email. I'm sorry if I got the ops wrong. v*v was actually rejected in favour of dot(v,v). Unless I've made another mistake the implemented operators are:
>
> v*s    v/s    v*=s
> v/=s   v+v    v-v
> v+=v   v-=v   m+m
> m-m    m*m    m*v
>
> (where v is a vector, m a matrix, s a scalar)
>
> The reason the feature set is very limited is I wanted to get feedback before implementing a large set of features and then having to change things. The idea was to implement just enough so that it would be clear where I was going. For example, I wasn't aware that compile time foreach was a bad way to unroll loops (perhaps I should build a mixin string?).

- Element wise *, etc are important operators that you need to cleanly support. Design-wise, since we can't do Matlab style operators, I prefer */+- to be element wise, and then use dot/cross/solve for matrix ops. This seems to avoid the confusion of as having some ops be matrix style and some ops be array style. Since D is unicode, you could also add the real inner product operator as an alias. Also, don't forget the outer product, aka broadcasting, i.e. v.dot(v.T).
- static foreach seems to add some overhead, although I don't know why. What you should really do is write some benchmarks and test static foreach, array ops and hand-unrolling for several vector sizes.

Robert Jacques:
> - Element wise *, etc are important operators that you need to cleanly support. Design-wise, since we can't do Matlab style operators, I prefer */+- to be element wise,

But keep in mind in D the built-in array ops always require the [].

Bye,
bearophile

On Mon, 26 Apr 2010 13:16:02 -0300, bearophile <bearophileHUGS@lycos.com> wrote:

> Robert Jacques:
>> - Element wise *, etc are important operators that you need to cleanly
>> support. Design-wise, since we can't do Matlab style operators, I prefer
>> */+- to be element wise,
>
> But keep in mind in D the built-in array ops always require the [].
>
> Bye,
> bearophile

Yes, but emulating that behavior is difficult. And array ops are still buggy and broken to the degree of not being usable. In theory you shouldn't need to define a small vector or matrix in D: fixed size arrays are value types, array op support and free functions as member function support. In practice, though, we still need it.

Robert Jacques:
> Yes, but emulating that behavior is difficult. And array ops are still buggy and broken to the degree of not being usable. In theory you shouldn't need to define a small vector or matrix in D: fixed size arrays are value types, array op support and free functions as member function support. In practice, though, we still need it.

You are right, in practice we need it. But practice is not enough, we have to keep in mind what's right too. And the right thing is to improve the built-in array ops :-)

And array ops don't replace everything you can do with a vector lib, you can define dot and cross product, small matrices, matrix operations, sparse arrays, truly rectangular arrays made of a single block of memory, etc.

Bye,
bearophile

> - Element wise *, etc are important operators that you need to cleanly support. Design-wise, since we can't do Matlab style operators, I prefer */+- to be element wise, and then use dot/cross/solve for matrix ops. This seems to avoid the confusion of as having some ops be matrix style and some ops be array style. Since D is unicode, you could also add the real inner product operator as an alias. Also, don't forget the outer product, aka broadcasting, i.e. v.dot(v.T).
> - static foreach seems to add some overhead, although I don't know why. What you should really do is write some benchmarks and test static foreach, array ops and hand-unrolling for several vector sizes.

I'm afraid you've just set me off a rant. I apologise, but you just told me that your web browser was Google.

I think we may be heading for a culture clash. I approach vectors and matrices from a math/physics background. Here vector's primary and only purpose is to model a directed arrow in some space (3D space, function space (lots and lots of things can be shoehorned into the idea of a space) but always in some space.). Consequently, the only things that matter about a vector are physical quantities. What are physical quantities? Anything that is invariant under a coordinate transform. For example, if you have a unit stick sitting due north and and you measure it with a ruler sitting due north and one due east you conclude this:

(1,0)

But if you measure with two rulers due NW and NE you conclude:

(0.707,0.707)

These are a same physical vector. You can see this by using the rotation  matrix between the coordinate systems:

(0.707 -0.707)(1)=(0.707)
(0.707  0.707)(0) (0.707)

So the vector coordinates (a coincidence of the placement of the rulers) are not physical but the vector is. I don't think I have never seen a physical equation requiring a element by element multiplication but I have seen lots of dot products. There is a reason for this. The dot product of a vector with another vector gives you a scalar. A scalar is not just a number (vector components are numbers but not scalars), a scalar is a number that does not change under a coordinate transform. The height of a hill does not change when you draw a map in different ways so it that is a scalar. The element by element product does not yield anything other than N numbers which are an accident of the coordinate system.

Matrices are either a system of linear equations or a representation of a tensor. Again, the most natural way to combine them is the matrix product (although there are other natural ways: direct product, kronchecker product and index contraction). The best way to think about a matrix is not to think about the matrix components, but as an indivisible whole. Element by element multiplication seems equally arbarant.

These things have been so deeply ingrained that I didn't even consider an element by element product in the design and I had to do a double take before I got what was being asked for. When the dot product was asked for as a "dot" I thought it would replace

I'd consider any use of element by element multiplication has probably a sign of a hack, and possibly a sign that what you are trying to use a vector for is actually an array of numbers. Hacks are fine but I don't think they're worth an operator. I'd want a big, explicit elementWiseProduct() to tell me something was going on.

Good user interfaces are modelled on the principle of least surprise. I think there are two different models of least surprise here. I want transcribe my matrix-vector equations directly into D and:

u=Ap+Bq   p,q,u vectors, A matrix, B matrix

Looks better as:
u=A*p+B*q;

Than:
u=matrix_vector_product(A,p) + matrix_vector_product(B,q)

But I think programmers are more used to using vector and array interchangeably. Possibly this is a reason why a vector struct is needed rather than defining free operators over just any old array. You want the operators to have their vector-meanings for vectors while retaining their field-meanings when dealing with fields like float, double, real etc.

Here's another one from the C++ standard library. std::complex is a pair of numbers. Why isn't operator* a element wise multiplication with an separate complex_product(.,.)? Because if you use std::complex you probably want to model complex numbers.

tl;dr
Google is not really your web browser, even though I know what you mean. An array is not really a vector, especially when  3D (or any other D) space .


....

PS: Okay so I just had a looked at the matrix and vector classes in Ogre3D and irrlicht. Looks like they both define v*v as element wise multiplication but m*m is matrix multiplication. That just seems even more inconsistent.

bearophile wrote:
> Gareth Charnock:
> 
>> http://github.com/gcharnock/phoboslinalgebra
> 
> I think that a single module is better. If seen fitting, some functionality can be moved inside other already present modules of Phobos.
> The module will need a good amount of unit tests.
> 
> In the code I see no point in calling functions like:
> CTFENextToken
> Just call them:
> nextToken
> Note: functions, template functions, and ctfe functions, start with a lower case letter (and templates generally with with an upper case).
> 
I should probably make this private as I don't see them being generally useful.

> Don't call this module and its contents matrix or vector or Vector, etc. Call them SmallVector, or SmallVec, SmallMat, ecc, because this lib is not designed to be a general purpose matrix or vector lib, it's designed for small number of items.
> You can call this module std.tinyarray :-)
> 
Good idea although I'm not sure tinyarray is a good name and I wouldn't look there for vectors and matrices. std.tinylinalg comes off as a mess of letters.

> 
>> alias Matrix!(float,3) matrix_t;
> 
> To define the sizes of a matrix you need two numbers:
> alias Matrix!(float, 3, 5) Matrix35;
> 
Fair enough.
> 
>> //Very natural initialisation
>> auto m=matrix_t(1,2,0,0,1,0,0,3,1);
> 
> A 1D initialization of a matrix is not so natural :-)

It is in emacs, which will align the statement like this:
auto m=matrix_t(1,2,0,
		0,1,0,
		0,3,1)

But point taken, there needs to be other ways of initialising.

> 
>> //No matter how big N, elements always have a canonical name (probably more important for matrices where the letters of the alphabet run out faster although I've only done this for vectors so far)
> 
> What's the point of giving names to matrix elements?

Ease of access, just in case you need it. Also stuff like:
vector3 translate = m.a30a31a32;
Most of the metaprograming is already done so it might as well be reused.
(Actually that's so common it could do with its own function.)

> 
>> 2) Set a bool. Switch between two element iteration schemes for nearly every operation (with suitable use of templates this might not be as painful as it sounds.)
> 
> The point of a transposition is sometimes to actually have data in a different position. With cache effects the position of data can be quite important (a matrix multiplication can be 2 or 3 times faster if you transpose. This is true for larger matrices).

Which is why using this design will increase the complexity of the code a fair bit. I'm not a fan of 2 because I'd prefer a value type. However I want to check I'm not alone (as I appear to be on the definition of * for vectors).

> 
>> Transposed:
>> 1) makes a new matrix
>> 2) creates a transposed image which is forever linked to the original matrix: A=B.Transpose() => Changes to A affect B
> 
> The semantics of transpose/transposed can be like in Python: the transpose works in-place, transposed creates a new matrix.

Sorry, should have read A=B.Transposed() => A becomes linked to B, changes to A affect B.

> Methods and properties like transpose have to start with a lowercase:
> auto A = B.transposed;

I'll try to be more consistent.

> Bye,
> bearophile

On Mon, 26 Apr 2010 18:02:53 -0300, Gareth Charnock <gareth.tpc@gmail.com> wrote:
>> - Element wise *, etc are important operators that you need to cleanly support. Design-wise, since we can't do Matlab style operators, I prefer */+- to be element wise, and then use dot/cross/solve for matrix ops. This seems to avoid the confusion of as having some ops be matrix style and some ops be array style. Since D is unicode, you could also add the real inner product operator as an alias. Also, don't forget the outer product, aka broadcasting, i.e. v.dot(v.T).
>> - static foreach seems to add some overhead, although I don't know why. What you should really do is write some benchmarks and test static foreach, array ops and hand-unrolling for several vector sizes.
>
> I'm afraid you've just set me off a rant. I apologise, but you just told me that your web browser was Google.
>
> I think we may be heading for a culture clash. I approach vectors and matrices from a math/physics background. Here vector's primary and only purpose is to model a directed arrow in some space (3D space, function space (lots and lots of things can be shoehorned into the idea of a space) but always in some space.). Consequently, the only things that matter about a vector are physical quantities. What are physical quantities? Anything that is invariant under a coordinate transform. For example, if you have a unit stick sitting due north and and you measure it with a ruler sitting due north and one due east you conclude this:
>
> (1,0)
>
> But if you measure with two rulers due NW and NE you conclude:
>
> (0.707,0.707)
>
> These are a same physical vector. You can see this by using the rotation   matrix between the coordinate systems:
>
> (0.707 -0.707)(1)=(0.707)
> (0.707  0.707)(0) (0.707)
>
> So the vector coordinates (a coincidence of the placement of the rulers) are not physical but the vector is. I don't think I have never seen a physical equation requiring a element by element multiplication but I have seen lots of dot products. There is a reason for this. The dot product of a vector with another vector gives you a scalar. A scalar is not just a number (vector components are numbers but not scalars), a scalar is a number that does not change under a coordinate transform. The height of a hill does not change when you draw a map in different ways so it that is a scalar. The element by element product does not yield anything other than N numbers which are an accident of the coordinate system.
>
> Matrices are either a system of linear equations or a representation of a tensor. Again, the most natural way to combine them is the matrix product (although there are other natural ways: direct product, kronchecker product and index contraction). The best way to think about a matrix is not to think about the matrix components, but as an indivisible whole. Element by element multiplication seems equally arbarant.
>
> These things have been so deeply ingrained that I didn't even consider an element by element product in the design and I had to do a double take before I got what was being asked for. When the dot product was asked for as a "dot" I thought it would replace
>
> I'd consider any use of element by element multiplication has probably a sign of a hack, and possibly a sign that what you are trying to use a vector for is actually an array of numbers. Hacks are fine but I don't think they're worth an operator. I'd want a big, explicit elementWiseProduct() to tell me something was going on.
>
> Good user interfaces are modelled on the principle of least surprise. I think there are two different models of least surprise here. I want transcribe my matrix-vector equations directly into D and:
>
> u=Ap+Bq   p,q,u vectors, A matrix, B matrix
>
> Looks better as:
> u=A*p+B*q;
>
> Than:
> u=matrix_vector_product(A,p) + matrix_vector_product(B,q)
>
> But I think programmers are more used to using vector and array interchangeably. Possibly this is a reason why a vector struct is needed rather than defining free operators over just any old array. You want the operators to have their vector-meanings for vectors while retaining their field-meanings when dealing with fields like float, double, real etc.
>
> Here's another one from the C++ standard library. std::complex is a pair of numbers. Why isn't operator* a element wise multiplication with an separate complex_product(.,.)? Because if you use std::complex you probably want to model complex numbers.
>
> tl;dr
> Google is not really your web browser, even though I know what you mean. An array is not really a vector, especially when  3D (or any other D) space .
>
>
> ....
>
> PS: Okay so I just had a looked at the matrix and vector classes in Ogre3D and irrlicht. Looks like they both define v*v as element wise multiplication but m*m is matrix multiplication. That just seems even more inconsistent.

Personally, I'd expect a pure mathematician to argue that the * operator isn't equivalent to the • operator and might even argue for the more formal <,> operator. :) Remember, unlike complex numbers, * isn't defined between vector types, so the library needs to decide what operator it's mapping to it. As for principal of least surprise, lots of libraries op for using array-wise * and v.dot(v), including (IIRC) all of the various GPU shading/compute languages.

P.S. I use Opera for my web browsing.

Gareth Charnock wrote:
> PS: Okay so I just had a looked at the matrix and vector classes in Ogre3D and irrlicht. Looks like they both define v*v as element wise multiplication but m*m is matrix multiplication. That just seems even more inconsistent.

	Eigen (http://eigen.tuxfamily.org/ ) uses '*' for the matrix
multiplication. v*v is an error (incompatible shapes). Element wise
operations can be done like this: v.cwise()*v

		Jerome
-- 
mailto:jeberger@free.fr
http://jeberger.free.fr
Jabber: jeberger@jabber.fr

> As for principal of least surprise, lots of libraries op for using array-wise * and v.dot(v), including (IIRC) all of the various GPU shading/compute languages.


In GLSL, vector * vector is component-wise
      but  matrix * vector is not
             matrix * matrix is not,

and it feels very natural.

In HLSL, * is always component-wise and you have to use mul(a, b) to use the product... I don't buy it.

Gareth is right, element-wise multiplication for matrices is near to useless, whereas it's handy for vectors.

Jérôme M. Berger wrote:
> Gareth Charnock wrote:
>> PS: Okay so I just had a looked at the matrix and vector classes in
>> Ogre3D and irrlicht. Looks like they both define v*v as element wise
>> multiplication but m*m is matrix multiplication. That just seems even
>> more inconsistent.
> 
> 	Eigen (http://eigen.tuxfamily.org/ ) uses '*' for the matrix
> multiplication. v*v is an error (incompatible shapes). Element wise
> operations can be done like this: v.cwise()*v
> 
> 		Jerome
That's an very impressive looking library if even half of what they claim is true. I can't believe I've missed this for so long (found the boost matrix library, blitz++, the matrix template library, FLENS, something called armadillo but not eigen. I do believe all my C++ linear algebra woes are over.

Eigen seems to treat vectors as 1 by n matrices and if you do this you get the matrix-vector product and the dot product for free as these are all the same operation. Probably v^T*v would be the dot product. Expression templates should make these operations efficient.