Hi, I have written a simple 2D matrix template where the extents of the matrix are template parameters. The multi_array container and the uBLAS matrix class both provide matrices of fixed dimension (template parameter for multi_array and 2 for uBLAS) but variable extents. I think my matrix definition is complementary to these. There are several advantages to having an array of fixed extents; but the implementation is, of course, unusable when the extents are not known at compile-time. If the extents of the matrix are part of the type, the compiler can enforce correct matrix extents for operations, avoiding the cost of runtime checks. For example: double b[3][2] = { ... }; matrix<3,3> A(2.0); // 2 * identity matrix<3,2> B(b); matrix<2,3> C; matrix<3,2> D; D = A * B; // OK C = A * B; // compile error - mismatched types C = transpose(B) * A; // OK Optimizing compilers can trivially unroll the loops used to define the operations. Use of this class does not require heap allocation. The storage required for the matrix is the minimal possible (assuming no assumptions about the values of the elements.) An array of matrices can be "serialized" to an array of doubles with a single reinterpret_cast. The proposed implementation is unsuitable for very large matrices: defining the element storage as a contiguous array may be problematic, most of the performance gains become insignificant, and the implicit copying will become significant. I've attached my current implementation. In it's current form, it is unsuitable for inclusion in Boost. I will fix any required issues if there is any interest in including it in Boost. The things that need to be changed (that I'm aware of) are: o license - GPL is not acceptable for Boost o naming (namespace and class name) o change explict 'double' for element type to a template parameter Some other changes that should probably be done before inclusion: o implement inverse for other than 2x2 and 3x3 matrices o more efficient determinant implementation for larger matrices o LU factorization o QR factorization In timing tests using G++ 4.1 on a Linux/Pentium4 platform, the performance of the arithmetic operations was equivalent to those of a special 3x3 matrix implementation with all operations explicitly unrolled (the 'old' implementation this one was to replace.) Is there any interest in including such a matrix definition in Boost? Is there already an effort to include something equivalent? I searched back a few years in the mailing list archive and didn't notice anything. Any suggested changes (other than those I listed above) for the attached implementation?
I did something very much like this a while back, which was submitted as a potential jumping off point for a boost geometry library. It is in the boost vault under Math - Geometry (http://tinyurl.com/ymjsoz), along with two other submissions. (Mine is geometryJH_01.zip). -Jason On 11/7/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
Hi,
I have written a simple 2D matrix template where the extents of the matrix are template parameters. The multi_array container and the uBLAS matrix class both provide matrices of fixed dimension (template parameter for multi_array and 2 for uBLAS) but variable extents. I think my matrix definition is complementary to these. There are several advantages to having an array of fixed extents; but the implementation is, of course, unusable when the extents are not known at compile-time.
If the extents of the matrix are part of the type, the compiler can enforce correct matrix extents for operations, avoiding the cost of runtime checks. For example: double b[3][2] = { ... }; matrix<3,3> A(2.0); // 2 * identity matrix<3,2> B(b); matrix<2,3> C; matrix<3,2> D; D = A * B; // OK C = A * B; // compile error - mismatched types C = transpose(B) * A; // OK Optimizing compilers can trivially unroll the loops used to define the operations. Use of this class does not require heap allocation. The storage required for the matrix is the minimal possible (assuming no assumptions about the values of the elements.) An array of matrices can be "serialized" to an array of doubles with a single reinterpret_cast.
The proposed implementation is unsuitable for very large matrices: defining the element storage as a contiguous array may be problematic, most of the performance gains become insignificant, and the implicit copying will become significant.
I've attached my current implementation. In it's current form, it is unsuitable for inclusion in Boost. I will fix any required issues if there is any interest in including it in Boost. The things that need to be changed (that I'm aware of) are: o license - GPL is not acceptable for Boost o naming (namespace and class name) o change explict 'double' for element type to a template parameter Some other changes that should probably be done before inclusion: o implement inverse for other than 2x2 and 3x3 matrices o more efficient determinant implementation for larger matrices o LU factorization o QR factorization
In timing tests using G++ 4.1 on a Linux/Pentium4 platform, the performance of the arithmetic operations was equivalent to those of a special 3x3 matrix implementation with all operations explicitly unrolled (the 'old' implementation this one was to replace.)
Is there any interest in including such a matrix definition in Boost?
Is there already an effort to include something equivalent? I searched back a few years in the mailing list archive and didn't notice anything.
Any suggested changes (other than those I listed above) for the attached implementation?
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
Jason Hise wrote:
I did something very much like this a while back, which was submitted as a potential jumping off point for a boost geometry library. It is in the boost vault under Math - Geometry (http://tinyurl.com/ymjsoz), along with two other submissions. (Mine is geometryJH_01.zip).
I looked at your matrix implementation briefly, and have a few suggestions: o specialized implementations of inverse and determinant for common sizes (2x2, 3x3, maybe 4x4) will probably give better performance. o I got substantially better performance for matrix multiplies by explicitly unrolling the innermost loop. but I only tested on one platform, so that my not be a universal improvement. I think that if Boost includes such a matrix implementation, it should be in some general math primitives library, because it is applicable to more than just geometry. I think that my design of the matrix class is a little better because a vector is defined as a Nx1 matrix, rather than the matrix being defined as an array of vectors. This avoids the need for special operators for vector-matrix multiplication, is closer to the underlying mathematical concept, and makes the equivalency of a vector and column matrix explicit.
On 11/8/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
Jason Hise wrote:
I did something very much like this a while back, which was submitted as a potential jumping off point for a boost geometry library. It is in the boost vault under Math - Geometry (http://tinyurl.com/ymjsoz), along with two other submissions. (Mine is geometryJH_01.zip).
I think that my design of the matrix class is a little better because a vector is defined as a Nx1 matrix, rather than the matrix being defined as an array of vectors. This avoids the need for special operators for vector-matrix multiplication, is closer to the underlying mathematical concept, and makes the equivalency of a vector and column matrix explicit.
I suppose this makes sense from an efficiency standpoint. I prefered the matrix as a vector of vectors because it made defining algorithms in terms of elementary row operations cleaner, and because it allowed the more natural syntax of indexing the matrix as though it were a 2D array. There are two usability features that my implementation provides which yours does not. The first is being able to index a vector by name (_x, _y, _z, _w). This is easily implemented using a properly scoped enum. I used leading underscores for my component names as a convention to indicate that they were property names, though admittedly this naming convention could be debated. vector <float, 3> up; up[_x] = 0; up[_y] = 1; up[_z] = 0; Another thing that my code allows that yours does not is explicit construction of a vector from the number of parameters that matches the size (at least from 1 through 4), using the SFINAE techniques of boost.enable_if. It can be much nicer for client code to be able to write: vector <float, 3> up (0, 1, 0); than to have to write the code segment I used for my previous example. Just a few thoughts. -Jason
Jason Hise wrote:
On 11/8/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
Jason Hise wrote:
I did something very much like this a while back, which was submitted as a potential jumping off point for a boost geometry library. It is in the boost vault under Math - Geometry (http://tinyurl.com/ymjsoz), along with two other submissions. (Mine is geometryJH_01.zip).
I think that my design of the matrix class is a little better because a vector is defined as a Nx1 matrix, rather than the matrix being defined as an array of vectors. This avoids the need for special operators for vector-matrix multiplication, is closer to the underlying mathematical concept, and makes the equivalency of a vector and column matrix explicit.
I suppose this makes sense from an efficiency standpoint. I prefered the matrix as a vector of vectors because it made defining algorithms in terms of elementary row operations cleaner, and because it allowed the more natural syntax of indexing the matrix as though it were a 2D array.
Indexing the matrix as if it were a 2D array can be accomplished other ways. For example: private: Type data[R*C]; public: Type operator[](unsigned i) { return data+i*C; } const Type* operator[](unsigned i) const { return data+i*C; } To index it as if it were a 2D array, the operator[] must return something that also has an operator[]. An array of works as well as a vector. And I fail to see how row operations make the implementation cleaner or easier. Most operators can be implemented with a single loop in both cases. For example: private: Type data[R*C]; public: matrix<R,C,Type>& operator+=( const matrix<R,C,Type>& other ) { for (unsigned i = 0; i < R*C; ++i) data[i] += other.data[i]; return *this; }
There are two usability features that my implementation provides which yours does not. The first is being able to index a vector by name (_x, _y, _z, _w). This is easily implemented using a properly scoped enum. I used leading underscores for my component names as a convention to indicate that they were property names, though admittedly this naming convention could be debated.
vector <float, 3> up; up[_x] = 0; up[_y] = 1; up[_z] = 0;
I've never found such a convention very useful. It is a trivial thing to add if others do.
Another thing that my code allows that yours does not is explicit construction of a vector from the number of parameters that matches the size (at least from 1 through 4), using the SFINAE techniques of boost.enable_if. It can be much nicer for client code to be able to write:
vector <float, 3> up (0, 1, 0);
than to have to write the code segment I used for my previous example.
That is a useful feature. I had considered several ways of implementing it. However, I had doubts about all of them, so I didn't provide any. 1) varargs o matrix( Type ii, ... ); o no type checking on arguments 2) multiple constructors with different numbers of args o what you have implemented o seems conceptually wrong to have a four-arg constructor for a 2-element vector, or visa versa. o doesn't work for larger matrices 3) one constructor with unrolled assignment o most of the same limitations as 2) o one function o safe (will never write off the end of the array) o matrix( Type a, Type b , Type c = 0, Type d = 0, Type e = 0, Type f = 0, Type g = 0, Type h = 0, Type i = 0, Type j = 0, Type k = 0, Type l = 0, Type m = 0, Type n = 0, Type o = 0, Type p = 0 ) { switch (R*C) { case 16: data[15] = p; case 15: data[14] = o; ... case 1: data[ 0] = a; } o excessively large argument list for e.g. a 2x1 matrix - jason
On 11/9/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
I fail to see how row operations make the implementation cleaner or easier. Most operators can be implemented with a single loop in both cases. For example: private: Type data[R*C]; public: matrix<R,C,Type>& operator+=( const matrix<R,C,Type>& other ) { for (unsigned i = 0; i < R*C; ++i) data[i] += other.data[i]; return *this; }
I was thinking of the inverse algorithm specifically, and other algorithms in general. I'd like to distinguish here between a function which is a necessary part of the interface to the type, and algorithms which are separate extensions that operate on that type. Being able to perform elementary row operations on matricees is important (I just had to implement the simplex algorithm for a class I am taking). However, it may be more important to provide for matricees which do not have fixed sizes at compile time than for those that do (simplex required me to be able to insert rows and columns at run time).
There are two usability features that my implementation provides which yours does not. The first is being able to index a vector by name (_x, _y, _z, _w). This is easily implemented using a properly scoped enum. I used leading underscores for my component names as a convention to indicate that they were property names, though admittedly this naming convention could be debated.
vector <float, 3> up; up[_x] = 0; up[_y] = 1; up[_z] = 0;
I've never found such a convention very useful. It is a trivial thing to add if others do.
I think it falls into the same category as not using magic numbers. To give another example, if you use a vector to store a color, being able to refer to the components as r, g, b, and a channels instead of by number improves readability. But perhaps enumerating these constants is better left to client code.
Another thing that my code allows that yours does not is explicit construction of a vector from the number of parameters that matches the size (at least from 1 through 4), using the SFINAE techniques of boost.enable_if. It can be much nicer for client code to be able to write:
vector <float, 3> up (0, 1, 0);
than to have to write the code segment I used for my previous example.
That is a useful feature. I had considered several ways of implementing it. However, I had doubts about all of them, so I didn't provide any.
1) varargs o matrix( Type ii, ... ); o no type checking on arguments 2) multiple constructors with different numbers of args o what you have implemented o seems conceptually wrong to have a four-arg constructor for a 2-element vector, or visa versa. o doesn't work for larger matrices
This is why I used boost::enable_if... it makes it impossible to call a 4 arg constructor on a two component vector and vice versa. The functions can never be instantiated if the dimensions don't match. The code could probably be further improved using BOOST_PP to automatically generate an arbitrary number of constructors to match larger numbers of initializers.
3) one constructor with unrolled assignment o most of the same limitations as 2) o one function o safe (will never write off the end of the array) o matrix( Type a, Type b , Type c = 0, Type d = 0, Type e = 0, Type f = 0, Type g = 0, Type h = 0, Type i = 0, Type j = 0, Type k = 0, Type l = 0, Type m = 0, Type n = 0, Type o = 0, Type p = 0 ) { switch (R*C) { case 16: data[15] = p; case 15: data[14] = o; ... case 1: data[ 0] = a; } o excessively large argument list for e.g. a 2x1 matrix
- jason
Another alternative is to use creative operator overloading: vector<float, 3> up ((vec_builder<float>, 0, 1, 0)); or: vector<float, 3> up (vec_builder<float>(0)(1)(0)); or maybe: vector<float, 3> up; up << 0 << 1 << 0; though I do happen to prefer my former solution, which looks and acts like a constructor call. -Jason
Jason Hise wrote:
There are two usability features that my implementation provides which yours does not. The first is being able to index a vector by name (_x, _y, _z, _w). This is easily implemented using a properly scoped enum. I used leading underscores for my component names as a convention to indicate that they were property names, though admittedly this naming convention could be debated.
vector <float, 3> up; up[_x] = 0; up[_y] = 1; up[_z] = 0;
I've never found such a convention very useful. It is a trivial thing to add if others do.
I think it falls into the same category as not using magic numbers. To give another example, if you use a vector to store a color, being able to refer to the components as r, g, b, and a channels instead of by number improves readability. But perhaps enumerating these constants is better left to client code.
I also think that this is best left to the client code. The concept of a mathematical vector should be kept separate from any domain-specific meaning (Cartesian vector, color space, etc.) Also, I think its best not to pollute the namespace with constants that are so likely conflict with other names. - jason
I think that the two important issues where are: 1) Boost should have these types defined outside of any domain-specific library. I think they equivalent to other types such a complex numbers, and should exist as primitives for use in many different libraries. There seems to be some support for this from others: http://lists.boost.org/MailArchives/ublas/2006/11/1519.php 2) Should a vector be defined as a column matrix or a separate template. This is the principal differentiator between your implementation and mine. The other issues are just simple features that can be added to one or the other. My argument for a column matrix definition of a vector is that it is closer to the mathematical concept we are trying to model. All matrix operations (multiplication between vectors and matrices, transpose, etc) work for vectors as well as matrices w/out needing special vector forms of the operations. My understanding of your argument for vectors as a separate class is that it allows the matrix class to be defined as an array of vectors, where each vector is a row of the matrix. This, in turn, allows for simpler coding of matrix row operations for algorithms such as Gauss-Jordan elimination. My implementation also provides row access to matrices as 1xC matrices, and provides vector specific operations (vector product, dot product) for both Rx1 and 1xC matrices. The only drawback to my approach is that a reinterpret_cast is required to return a reference to a row as a 1xC matrix. It is portable to do so, as the underlying cast is just operating on array data (no alignment issues for individual members, no vtable, etc.) But it is ugly. Overall, 1) what I'd most like to see. I'm willing to contribute to whichever implementation ends up being used for 1). - jason Jason Hise wrote:
On 11/9/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
I fail to see how row operations make the implementation cleaner or easier. Most operators can be implemented with a single loop in both cases. For example: private: Type data[R*C]; public: matrix<R,C,Type>& operator+=( const matrix<R,C,Type>& other ) { for (unsigned i = 0; i < R*C; ++i) data[i] += other.data[i]; return *this; }
I was thinking of the inverse algorithm specifically, and other algorithms in general. I'd like to distinguish here between a function which is a necessary part of the interface to the type, and algorithms which are separate extensions that operate on that type. Being able to perform elementary row operations on matricees is important (I just had to implement the simplex algorithm for a class I am taking). However, it may be more important to provide for matricees which do not have fixed sizes at compile time than for those that do (simplex required me to be able to insert rows and columns at run time).
There are two usability features that my implementation provides which yours does not. The first is being able to index a vector by name (_x, _y, _z, _w). This is easily implemented using a properly scoped enum. I used leading underscores for my component names as a convention to indicate that they were property names, though admittedly this naming convention could be debated.
vector <float, 3> up; up[_x] = 0; up[_y] = 1; up[_z] = 0;
I've never found such a convention very useful. It is a trivial thing to add if others do.
I think it falls into the same category as not using magic numbers. To give another example, if you use a vector to store a color, being able to refer to the components as r, g, b, and a channels instead of by number improves readability. But perhaps enumerating these constants is better left to client code.
Another thing that my code allows that yours does not is explicit construction of a vector from the number of parameters that matches the size (at least from 1 through 4), using the SFINAE techniques of boost.enable_if. It can be much nicer for client code to be able to write:
vector <float, 3> up (0, 1, 0);
than to have to write the code segment I used for my previous example. That is a useful feature. I had considered several ways of implementing it. However, I had doubts about all of them, so I didn't provide any.
1) varargs o matrix( Type ii, ... ); o no type checking on arguments 2) multiple constructors with different numbers of args o what you have implemented o seems conceptually wrong to have a four-arg constructor for a 2-element vector, or visa versa. o doesn't work for larger matrices
This is why I used boost::enable_if... it makes it impossible to call a 4 arg constructor on a two component vector and vice versa. The functions can never be instantiated if the dimensions don't match. The code could probably be further improved using BOOST_PP to automatically generate an arbitrary number of constructors to match larger numbers of initializers.
3) one constructor with unrolled assignment o most of the same limitations as 2) o one function o safe (will never write off the end of the array) o matrix( Type a, Type b , Type c = 0, Type d = 0, Type e = 0, Type f = 0, Type g = 0, Type h = 0, Type i = 0, Type j = 0, Type k = 0, Type l = 0, Type m = 0, Type n = 0, Type o = 0, Type p = 0 ) { switch (R*C) { case 16: data[15] = p; case 15: data[14] = o; ... case 1: data[ 0] = a; } o excessively large argument list for e.g. a 2x1 matrix
- jason
Another alternative is to use creative operator overloading:
vector<float, 3> up ((vec_builder<float>, 0, 1, 0));
or:
vector<float, 3> up (vec_builder<float>(0)(1)(0));
or maybe:
vector<float, 3> up; up << 0 << 1 << 0;
though I do happen to prefer my former solution, which looks and acts like a constructor call.
-Jason
Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
2) Should a vector be defined as a column matrix or a separate template.
If you are ever going to do objects with rank three or four (e.g. tensors), then it really does not make sense to make vectors just be column matrices. So my feeling is that vectors==column matrix is unnecessarily limiting. Cheers, Walter Landry wlandry@ucsd.edu
On Fri, 10 Nov 2006, Jason Kraftcheck wrote:
I think that the two important issues where are:
1) Boost should have these types defined outside of any domain-specific library.
I think they equivalent to other types such a complex numbers, and should exist as primitives for use in many different libraries. There seems to be some support for this from others: http://lists.boost.org/MailArchives/ublas/2006/11/1519.php
I personally agree that Boost should have this type of matrices and vectors. But for one thing, instead to designing a new interface or library, I would seriously think about integrating static sized matrices and vectors into uBLAS. I don't see any good reasons why we shouldn't mix the two kinds. -- François Duranleau LIGUM, Université de Montréal
François Duranleau wrote:
On Fri, 10 Nov 2006, Jason Kraftcheck wrote:
I think that the two important issues where are:
1) Boost should have these types defined outside of any domain-specific library.
I think they equivalent to other types such a complex numbers, and should exist as primitives for use in many different libraries. There seems to be some support for this from others: http://lists.boost.org/MailArchives/ublas/2006/11/1519.php
I personally agree that Boost should have this type of matrices and vectors. But for one thing, instead to designing a new interface or library, I would seriously think about integrating static sized matrices and vectors into uBLAS. I don't see any good reasons why we shouldn't mix the two kinds.
It would, of course, be best if any static-size matrix implementation worked with uBLAS. But it should be possible to have the matrix definition to exist outside uBLAS. It is a fairly fundamental type. I think it'd be nice to be able to use it without depending on the rest of uBLAS. - jason
On 11/10/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote:
I think that the two important issues here are:
1) Boost should have these types defined outside of any domain-specific library.
I think they are equivalent to other types such as complex numbers, and should exist as primitives for use in many different libraries. There seems to be some support for this from others: http://lists.boost.org/MailArchives/ublas/2006/11/1519.php
Vectors do seem to be important as a primitive foundation for many things, including complex numbers, quaternions, colors, and matricees ;)
2) Should a vector be defined as a column matrix or a separate template.
This is the principal differentiator between your implementation and mine. The other issues are just simple features that can be added to one or the other.
My argument for a column matrix definition of a vector is that it is closer to the mathematical concept we are trying to model. All matrix operations (multiplication between vectors and matrices, transpose, etc) work for vectors as well as matrices w/out needing special vector forms of the operations.
Although I understand the desire to merge the matrix and vector into a single type in order to shrink the interface that must be implemented, and though it may make sense from a mathematical and/or efficiency standpoint, I do not think that it makes sense from a design or extensibility standpoint. Good designs tend to compose more complex concepts from simpler ones, rather than defining simpler concepts as special cases of more complex ones. Generally, programmers do not tend to view floats as complex numbers that happen to have a zero imaginary component, or view ints as floats that have nothing following the decimal point, even though that is the mathematic concept being modeled. If I were defining a color class based on a vector, I don't think that it would occur to me that I was defining a specialized matrix that only had one row. Also, consider the library from a dependency standpoint. People that need vectors may well not need matricees. But people that are using matricees are very likely to need vectors as well.
My understanding of your argument for vectors as a separate class is that it allows the matrix class to be defined as an array of vectors, where each vector is a row of the matrix. This, in turn, allows for simpler coding of matrix row operations for algorithms such as Gauss-Jordan elimination.
Yes, but it is more than that. It seems to me that a vector is not just a specialized matrix, but is a type in its own right which supports its own operations. Matrix multiplication is not defined between two vectors of the same orientation, and vector multiplication does not apply to matricees.
My implementation also provides row access to matrices as 1xC matrices, and provides vector specific operations (vector product, dot product) for both Rx1 and 1xC matrices.
This means that you are supporting two different types of vectors. So which type should be used as the foundation for vector based types, like complex numbers and colors?
The only drawback to my approach is that a reinterpret_cast is required to return a reference to a row as a 1xC matrix. It is portable to do so, as the underlying cast is just operating on array data (no alignment issues for individual members, no vtable, etc.) But it is ugly.
This issue arises directly because of the fact that a matrix conceptually is a compound type, whose components should be acessible to client code. But a matrix is not just a compound type composed of elements. A matrix is composed of rows of elements. -Jason
On Fri, 10 Nov 2006, Jason Hise wrote: [snip]
This issue arises directly because of the fact that a matrix conceptually is a compound type, whose components should be acessible to client code. But a matrix is not just a compound type composed of elements. A matrix is composed of rows of elements.
I agree on your stand about matrix vs vector types as being distinct, but I don't agree that a matrix is just composed of rows of elements. It could also be columns of elements. Which is it? A matrix should be left as a two dimensional array (conceptually), either row-major or column-major or the exact layout could be a parameter, as for boost::multi_array, except this time maybe as a template parameter. -- François Duranleau LIGUM, Université de Montréal
On 11/10/06, François Duranleau <duranlef@iro.umontreal.ca> wrote:
On Fri, 10 Nov 2006, Jason Hise wrote:
[snip]
This issue arises directly because of the fact that a matrix conceptually is a compound type, whose components should be accessible to client code. But a matrix is not just a compound type composed of elements. A matrix is composed of rows of elements.
I agree on your stand about matrix vs vector types as being distinct, but I don't agree that a matrix is just composed of rows of elements. It could also be columns of elements. Which is it? A matrix should be left as a two dimensional array (conceptually), either row-major or column-major or the exact layout could be a parameter, as for boost::multi_array, except this time maybe as a template parameter.
Maybe. I have to wonder though if this flexibility adds functionality, or just room for confusion. I've heard of doing elementary row opeerations, but I have not heard of doing elementary column operations. I would think that a standardized means of storage and accesss would make code clearer, but I will concede this point if you can show me use cases where this flexibility comes in handy for client code. -Jason
On Fri, 10 Nov 2006, Jason Hise wrote:
On 11/10/06, François Duranleau <duranlef@iro.umontreal.ca> wrote:
On Fri, 10 Nov 2006, Jason Hise wrote:
I agree on your stand about matrix vs vector types as being distinct, but I don't agree that a matrix is just composed of rows of elements. It could also be columns of elements. Which is it? A matrix should be left as a two dimensional array (conceptually), either row-major or column-major or the exact layout could be a parameter, as for boost::multi_array, except this time maybe as a template parameter.
Maybe. I have to wonder though if this flexibility adds functionality, or just room for confusion. I've heard of doing elementary row opeerations, but I have not heard of doing elementary column operations. I would think that a standardized means of storage and accesss would make code clearer, but I will concede this point if you can show me use cases where this flexibility comes in handy for client code.
I can think of a few. One is when you do graphics with OpenGL (I don't know how it is with other APIs). The latter uses column-major format for matrices. Thus if you need to often interface with OpenGL, it's better (faster) if you're internal representation is in the same order (as a result, I often see matrix implementations in column-major order instead of standard C ordering in the graphics community). Another one, still in graphics (or actually, any applications dealing with geometry), concerns rotation matrices. If you have a basis of three orthonormal vectors (assuming 3D here), you can get a rotation matrix by assign two vectors as the columns or rows of a 3x3 matrix, depending on the direction of the transformation (local coordinates to world or vice-versa). It is also often useful to extract back those vectors (columns or rows). Also, when you do some linear solving, if you want to implement a scheme with full pivoting (rows and columns), you have more than row operations. Concerning the cost of that flexibility, well, after posting my previous reply, I went back to my matrix implementation and see what I had to change to had a choice of ordering as a template parameter. I was surprised that I had only a few changes to do. The only thing that poses some problem is when applying an operator on two matrices of different ordering: what should be the returned ordering type? -- François Duranleau LIGUM, Université de Montréal
On 11/10/06, Jason Hise <0xchaos@gmail.com> wrote:
On Fri, 10 Nov 2006, Jason Hise wrote:
[snip]
This issue arises directly because of the fact that a matrix conceptually is a compound type, whose components should be accessible to client code. But a matrix is not just a compound type composed of elements. A matrix is composed of rows of elements.
I agree on your stand about matrix vs vector types as being distinct, but I don't agree that a matrix is just composed of rows of elements. It could also be columns of elements. Which is it? A matrix should be left as a two dimensional array (conceptually), either row-major or column-major or
exact layout could be a parameter, as for boost::multi_array, except
On 11/10/06, François Duranleau <duranlef@iro.umontreal.ca> wrote: the this
time maybe as a template parameter.
Maybe. I have to wonder though if this flexibility adds functionality, or just room for confusion. I've heard of doing elementary row opeerations, but I have not heard of doing elementary column operations. I would think that a standardized means of storage and accesss would make code clearer, but I will concede this point if you can show me use cases where this flexibility comes in handy for client code.
-Jason _______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
above certain rings sometimes it is necessary to perform elementary column operations rather then row operations, something that i run into when i attempted to build my own matrix implementation of fixed dimensions (sadly before i was aware of boost) -- -- lee Lee Elenbaas lee.elenbaas@gmail.com
Jason Hise wrote:
On 11/10/06, Jason Kraftcheck <kraftche@cae.wisc.edu> wrote: ...
2) Should a vector be defined as a column matrix or a separate template.
This is the principal differentiator between your implementation and mine. The other issues are just simple features that can be added to one or the other.
My argument for a column matrix definition of a vector is that it is closer to the mathematical concept we are trying to model. All matrix operations (multiplication between vectors and matrices, transpose, etc) work for vectors as well as matrices w/out needing special vector forms of the operations.
Although I understand the desire to merge the matrix and vector into a single type in order to shrink the interface that must be implemented, and though it may make sense from a mathematical and/or efficiency standpoint, I do not think that it makes sense from a design or extensibility standpoint. Good designs tend to compose more complex concepts from simpler ones, rather than defining simpler concepts as special cases of more complex ones. Generally, programmers do not tend to view floats as complex numbers that happen to have a zero imaginary component, or view ints as floats that have nothing following the decimal point, even though that is the mathematic concept being modeled. If I were defining a color class based on a vector, I don't think that it would occur to me that I was defining a specialized matrix that only had one row.
There are obvious performance, storage, and implementation issues for the analogies you give. If these issues didn't exist, I don't think we'd have the distinct types. Defining a vector as 1-column matrix doesn't have any such issues. If the vector-as-matrix is implemented correctly, you should never need to know that the vector is a matrix unless you're using it with matrices.
Also, consider the library from a dependency standpoint. People that need vectors may well not need matricees. But people that are using matricees are very likely to need vectors as well.
If all the linear-algebra type algorithms such as matrix factorizations, inversion, etc. were moved to a separate file, the number of operations that apply to matrices other than vectors is trivially small. Addition, subtraction, scalar operators, etc. all share a single implementation for vectors and matrices. The 'extra' dependency on matrices is only the multiplication implementation and some methods for accessing data.
My understanding of your argument for vectors as a separate class is that it allows the matrix class to be defined as an array of vectors, where each vector is a row of the matrix. This, in turn, allows for simpler coding of matrix row operations for algorithms such as Gauss-Jordan elimination.
Yes, but it is more than that. It seems to me that a vector is not just a specialized matrix, but is a type in its own right which supports its own operations. Matrix multiplication is not defined between two vectors of the same orientation, and vector multiplication does not apply to matricees.
Matrix multiplication is still applicable for vectors. matrix<double,3,1> v1, v2; matrix<double,3,3> M; ... v2 = M * v1; // matrix multiplication v2 = transpose(v1) * M; // same as above M = v1 * transpose(v2); // outer product There are a few special cases for vector multiplication: - vector (cross) products - inner (dot) product returns a scalar rather than a 1x1 matrix: double d = v1 % v1; is much more convenient and readable than double d = (transpose(v1) * v2)(0,0); But this implies that the definition of vectors requires additional operations beyond those of a matrix, not the other way around. I think that defining matrices in terms of vectors introduces a lot more junk into the matrix definition than defining vectors as matrices does for the vector definition. Anyway, multiplication is the only special case. All the other operations apply to both (addition, scalar multiply, etc.)
My implementation also provides row access to matrices as 1xC matrices, and provides vector specific operations (vector product, dot product) for both Rx1 and 1xC matrices.
This means that you are supporting two different types of vectors. So which type should be used as the foundation for vector based types, like complex numbers and colors?
Column matrices: template<Typename T, unsigned l> class svector : public smatrix<T,l,1> { ... T& operator()( unsigned i } { return smatrix<T,l,1>::operator()(i,0); } ... } I provide the other because sometimes it is useful to have "row" vectors rather than "column" vectors. Why impose artificial limitations? Anyway, for the specific case of row operations for elimination, why does it matter if the rows are vectors or matrices. Operations like scalar multiply and addition are the same for either.
The only drawback to my approach is that a reinterpret_cast is required to return a reference to a row as a 1xC matrix. It is portable to do so, as the underlying cast is just operating on array data (no alignment issues for individual members, no vtable, etc.) But it is ugly.
This issue arises directly because of the fact that a matrix conceptually is a compound type, whose components should be acessible to client code. But a matrix is not just a compound type composed of elements. A matrix is composed of rows of elements.
Or columns of elements, or submatrices, or whatever... I don't think there's anything special about rows of elements. I had a bit of trouble deciding if I should define my matrix in row-major or column-major order. The former corresponds to the normal C array layout and allows for easier row operations, but for my work I often need to work with the columns of matrices. - jason
Hi, we recently investigated available C++ matrix libraries in the context of a new implementation of our Finite Element library eXperiment FeLyX. Within the element formulations (e.g. evaluation of stiffness matrices) of such FE libraries there are lots of fixed-size matrix operations where the sizes are known at compile time. Usually the matrix sizes are given through the actual element type used. We looked at MTL, UBLAS and TVMET (Tiny Vector Matrix library using Expression Templates, http://tvmet.sourceforge.net/), and finally decided to use TVMET. The reasons were: - As far as I understand your 2D-matrix template, TVMET offers exactly this type of fixed-size matrix library you are describing, providing compile-time checks for matrix sizes etc. - TVMET is fast. For matrices from size 3x3 up to app. 60x60 we did a very basic performance comparison (Matrix-Matrix / Matrix-Vector products) with MTL, UBLAS and TVMET (Tiny Vector Matrix library using Expression Templates, http://tvmet.sourceforge.net/). We were really impressed by the performance of TVMET, in our basic test it was significantly faster than the other two. - The operator overlads of TVMET (implemented using expression templates) provide a very convenient and easy-to-use interface (especially compared to MTL which we used before). I think TVMET could represent an interesting base for such a library within BOOST? Sincerely ok ----------- Oliver Koenig EVEN - Evolutionary Engineering AG www.even-ag.ch Jason Kraftcheck wrote:
Hi,
I have written a simple 2D matrix template where the extents of the matrix are template parameters. The multi_array container and the uBLAS matrix class both provide matrices of fixed dimension (template parameter for multi_array and 2 for uBLAS) but variable extents. I think my matrix definition is complementary to these. There are several advantages to having an array of fixed extents; but the implementation is, of course, unusable when the extents are not known at compile-time.
If the extents of the matrix are part of the type, the compiler can enforce correct matrix extents for operations, avoiding the cost of runtime checks. For example: double b[3][2] = { ... }; matrix<3,3> A(2.0); // 2 * identity matrix<3,2> B(b); matrix<2,3> C; matrix<3,2> D; D = A * B; // OK C = A * B; // compile error - mismatched types C = transpose(B) * A; // OK Optimizing compilers can trivially unroll the loops used to define the operations. Use of this class does not require heap allocation. The storage required for the matrix is the minimal possible (assuming no assumptions about the values of the elements.) An array of matrices can be "serialized" to an array of doubles with a single reinterpret_cast.
The proposed implementation is unsuitable for very large matrices: defining the element storage as a contiguous array may be problematic, most of the performance gains become insignificant, and the implicit copying will become significant.
I've attached my current implementation. In it's current form, it is unsuitable for inclusion in Boost. I will fix any required issues if there is any interest in including it in Boost. The things that need to be changed (that I'm aware of) are: o license - GPL is not acceptable for Boost o naming (namespace and class name) o change explict 'double' for element type to a template parameter Some other changes that should probably be done before inclusion: o implement inverse for other than 2x2 and 3x3 matrices o more efficient determinant implementation for larger matrices o LU factorization o QR factorization
In timing tests using G++ 4.1 on a Linux/Pentium4 platform, the performance of the arithmetic operations was equivalent to those of a special 3x3 matrix implementation with all operations explicitly unrolled (the 'old' implementation this one was to replace.)
Is there any interest in including such a matrix definition in Boost?
Is there already an effort to include something equivalent? I searched back a few years in the mailing list archive and didn't notice anything.
Any suggested changes (other than those I listed above) for the attached implementation?
Oliver Koenig wrote:
Hi, we recently investigated available C++ matrix libraries in the context of a new implementation of our Finite Element library eXperiment FeLyX. Within the element formulations (e.g. evaluation of stiffness matrices) of such FE libraries there are lots of fixed-size matrix operations where the sizes are known at compile time. Usually the matrix sizes are given through the actual element type used.
We looked at MTL, UBLAS and TVMET (Tiny Vector Matrix library using Expression Templates, http://tvmet.sourceforge.net/), and finally decided to use TVMET. The reasons were: - As far as I understand your 2D-matrix template, TVMET offers exactly this type of fixed-size matrix library you are describing, providing compile-time checks for matrix sizes etc. - TVMET is fast. For matrices from size 3x3 up to app. 60x60 we did a very basic performance comparison (Matrix-Matrix / Matrix-Vector products) with MTL, UBLAS and TVMET (Tiny Vector Matrix library using Expression Templates, http://tvmet.sourceforge.net/). We were really impressed by the performance of TVMET, in our basic test it was significantly faster than the other two. - The operator overlads of TVMET (implemented using expression templates) provide a very convenient and easy-to-use interface (especially compared to MTL which we used before).
Thanks for the info. I looked at this implementation briefly, and I think my implementation is better for the reasons I outline below. But as a said, I only looked at it briefly, so it is possible that I missed or mis-understood some portion of it. I'd greatly appreciate any corrections to my statements about TVMET below. I think that most of my original operator implementations are equivalent to those of TVMET. I later tuned the multiply operation. I'm sure that for the specific case of multiplying 3x3 matrices on an Intel Netburst-architecture CPU with G++ >= 4.0, my implementation will be faster than TVMET's. I think it will be faster for many other cases also, but that's just a guess. The TVMET library seems to define implementations of every libm function for matrices (sin, pow, etc.). It implements them by applying the equivalent scalar function to each matrix element. This kind of stuff results in a lot of noise in the library, is unlikely to be very useful, and the theoretical meaning of these operations is dubious. For example, mathematically, pow( M, 3 ) for a matrix M means M * M * M, which is not the same as cubing each element of M. I didn't notice any common square matrix operations such as inverse, determinant, etc. Did I miss these? Also, see my other message about a matrix as an array of vectors vs. a vector as a Nx1 matrix. - jason
participants (6)
-
François Duranleau -
Jason Hise -
Jason Kraftcheck -
Lee Elenbaas -
Oliver Koenig -
Walter Landry