Vectors and Matrices

The topics covered in this section include:

Overview

SimTK has two different sets of classes for vector and matrix objects. You have seen types Vec3 (a 3-vector) and Vector (a variable-length vector) in OpenSim examples. Here we will discuss those in more detail. If you want to know more about the design goals and implementation of these classes, see the SimTK Simmatrix document here: https://simtk.org/home/simbody, Documents tab.

First, there are classes to represent small, fixed size vectors and matrices with zero runtime overhead: Vec for column vectors, and Mat for matrices There is also a Row type that does not normally appear in user programs.. These classes are templatized based on size and element type. Synonyms (typedefs) are defined for common combinations; for example, Vec3 is a synonym for Vec<3,double>, while Mat22 is a synonym for Mat<2,2,double>. You can also create other combinations, such as Mat<2,10,double> or Vec<4,std::complex<double>>. However, the size must always be determinable at compile time. The in-memory representation of these small objects is minimal: only the data elements are stored.

Second, there are classes to represent large vectors and matrices whose sizes are determined at runtime: Vector_ for column vectors and Matrix_ for matrices Again with a normally-hidden RowVector_ type.. These classes are templatized based on element type. Types Vector and Matrix are synonyms for Vector_<double> and Matrix_<double>. Again, you can use other element types. In fact, the element type can even be one of the fixed-size vector or matrix objects. For example, Vector_<Vec3> is a vector, where each element is itself a three-component vector. The type Vec<2,Vec3>, called a spatial vector, is useful for combining rotational and translational quantities into a single object representing a spatial velocity or spatial force, for example. However, it is not permissible to use the variable-size Vector_ or Matrix_ objects as element types. The in-memory representation of these objects includes, in addition to the data, an opaque descriptor containing the length and information on how the data is laid out; the declared objects actually consist only of a pointer (essentially a void*) to the descriptors. This has many advantages for implementation and binary compatibility, but makes it difficult to look through these objects in a debugger as you can with the small Vec and Mat classes.

Here are some sample declarations:

Vec<3> v; // a 3-vector of RealsVec3
w; // same thing, using abbr.
Vector b,x; // vectors of doubles
Matrix M; // mxn matrix of doubles
Matrix_<Complex> C; // mxn matrix of complex<double>
Vector_<Vec3> v3; // big vector of 3-vecs
 
// This type is a 2-element vector whose elements// are 3-vectors. Memory layout and computational// efficiency are identical to Vec<6>.
typedef Vec<2,Vec3> SpatialVec; 

Operators

All of these classes support standard mathematical operators like +, -, *, / and C-style assignment operators like +=, -=, *=, /=. In addition, SimTK overloads the ~ unary operator to indicate transpose (or more precisely, Hermitian conjugate). That is, for any vector or matrix x, SimTK's "~x" has the same meaning as Matlab's "x'". For Vec3 there is also a cross product operator % available so that you can write compact expressions like

Vec3 w, r; // defined somewhere
Vec3 a = w % (w % r); // a=w x (w x r)

Like Matlab, SimTK requires strict shape conformance for vector and matrix arguments. So for Vec3 v and w, ~v*w is a dot product (scalar=row*vector) while v*~w is an outer product (matrix=vector*row), while v*w fails to compile because the arguments are not conforming. Scalar multiplication acts as expected. Global functions dot(), cross(), and outer() are available for those who prefer them to using the operators.

There also are versions of many standard math functions that operate on vectors and matrices: sin(), exp(), sqrt(), etc. and additional functions abs(), min(), max(), sort(), mean(), median(). These allow many calculations to be written in a very concise way.

Construction and assignment

All vector and matrix types define a default constructor, that is, a constructor with no arguments. In Debug mode, the default constructor initializes all elements to NaN. In Release mode, all elements are left uninitialized.

Constructors are also available for initializing data elements from individual element values, or by copying compatible objects. Initialization values can be provided in the constructor or via a pointer (or C array) to values of the appropriate type. Assignment operators are available for copying whole objects, and for setting single elements, subvectors and submatrices.

One convention followed by SimTK, which is different from that of most similar systems, is the treatment of scalar assignment. We follow this convention: (1) when a scalar s is assigned to a vector, every element of the vector is set to s (this is the typical convention), and (2) when a scalar s is assigned to a matrix, the diagonal elements of the matrix are set to s while the off-diagonals are set to zero. Examples:

Vec3 v;
Mat22 m;
Vec3 v(0); // v=(0,0,0)
v=0; // "
v=3; // v=(3,3,3)
Mat33 m(0); // m=( 0,0 ) zero
m=0 // ( 0,0 )
Mat33 m(1); // m=( 1,0 ) identity
m=1; // ( 0,1 ) 
Vector b(10); // initial size 10 doubles
Matrix M(20,10); // initial size 20x10
b=0; // b=10 zeroes
M=1; // M=0, except M(i,i)=1, 0<=i<10

This convention is especially apt for matrices, because the matrix resulting from such a scalar assignment "acts like" that scalar. That is, if you multiply by this matrix the result is identical to a scalar multiply by the original scalar. Two important special cases enabled by this convention are: (1) setting a matrix to the scalar "1" results in the multiplicative identity matrix of that shape, and (2) setting a matrix to the scalar "0" results in the additive identity matrix of that shape. In general, in any operation involving a scalar s and a Matrix or Mat, the scalar is treated as if it were a conforming matrix whose main diagonal consists of all s's with all other elements zero. So Matrix m += s will result in s being added to m's diagonal, which is what would happen if s were replaced by diag(s) of the same dimension as m. m =1 thus subtracts an identity matrix from m, without touching any of the off-diagonal elements. Note that for multiply and divide this convention yields the ordinary scalar multiply and divide operations: ms (=m*diag(s)) multiplies every element of m by s, while m/s (=m (1/s)) divides every element of m by s.

Getting the Size of a Vector, Matrix, or Array

The example below shows how to get the size of a vector, matrix, or array:

Vector v(2, SimTK::Pi);
int length;
length = v.size();

Matrix m(Mat22(1,2,
			   3,4));
int numRows = m.nrow();
int numCols = m.ncol();

Array_<bool> a;
a.push_back(true);
int numElements = a.size();

Indexing

SimTK provides 0-based indexing using the [] operator. If a Matrix is modifiable (non-const), then the indexed element can be modified and that change affects the contents of the object. The [] operator applied to a Matrix returns a row, which may in turn be indexed to obtain an element in C style. SimTK also permits indexing using round brackets () yielding identical results to [] for Vector but selecting a column rather than a row when applied to a Matrix. A two-argument round bracket operator accesses a Matrix element.

Matrix m; Vector v; …
v[i] // ref to ith element of v, 0-based 
v(i) // same 
m[i][j] // ref to i,jth element of m, 0-based 
m(i,j) // same, but faster 
m[i] // ref to ith row of m, 0-based 
m(j) // ref to jth column of m, 0-based

For the variable-size Matrix and Vector classes only, there are also operators for selecting sub-vectors and sub-matrices. Like the indexing operators, these return references into the original object, not copies. Sub-matrices are thus "lvalues" (in C terminology), meaning that they can appear on the left hand side of an assignment.

Matrix m; Vector v; …
v(i,m) // ref to m-element subvector whose 0th element is v's ith element 
 m(i,j,m,n) // ref to mxn submatrix whose (0,0) element is m's (i,j) element

Output

The C++ operator << is overloaded for all the matrix and vector types so you can look at a human-readable version of their contents via statements like

Matrix m;
std::cout << "m=" << m;