The second shoe drops
It's time to wrap up the series on matrices. I think I'm finally close enough to be able to give you a matrix class that's serviceable.
Because this series has been strung out so long (see the sidebar–the first matrix article I wrote was “Who needs matrices?” in the December 2007 issue of Embedded Systems Design) , I'll take the time to recap where we came from and where we've been.
For decades, I've been wanting a language that would let me do arithmetic with vectors and matrices, like this:
Vector x, y, u; Matrix A, B; ... y = A*x + B*u;
I first suggested a vector class in the August 2006 issue (“Loose ends,” page 14). In it, I showed a usable structure representing the mathematical entity called a vector. Over the next few issues, I developed the class until I now consider it complete. During the process, I showed you a few cases where use of vectors can make life bearable.
In the December 2007 issue (“Who Needs Matrices”–again I refer you to the sidebar for links), we wrapped up the Vector class and began the definition of a matrix class. Consider this the last issue in that series. At this point, I think we are close enough to a complete matrix class that we can wrap it up and declare victory. This is not to say that there is nothing left to do. The world of vectors and matrices is a rich, rich world, with a huge set of operations one can do with/on them. I could never in my lifetime give you code for every kind of manipulation one can do with them. Even then, there would still be the ones I've never even heard of.
A syllabus: matrix and vector math Jack Crenshaw has devoted many of his columns for the last 2 1/2 years to vector and matrix math. Here's a list of what he's written on vector and matrix math, with links to the online version. “Loose ends,” “Motivationally speaking,” October 2006, pg. 15, www.embedded.com/193402564 “Making the tough coding decisions,” December 2006, pg. 9, www.embedded.com/196601281 “SimpleVec wrap,” February 2007, pg. 11, www.embedded.com/197002378 “On to objects,” April 2007, pg. 13, www.embedded.com/198701651 “The vector class,” June 2007, pg. 13, www.embedded.com/199703501 “Completing the vector class,” August 2007, pg.13, www.embedded.com/201202113 “Vectors–the last word,” October 2007, pg. 9, http://www.embedded.com/202102663 “Who needs matrices,” December 2007, pg. 11, www.embedded.com/204300487 “Why multiply matrices?” February 2008, pg 9, www.embedded.com/205918951 “The matrix reprogrammed,” April 2008, pg. 9, “Fleshing out the matrix class,” June 2008, pg. 14, www.embedded.com/208400912 “Reusing code and inverting matrices,” August 2008, pg. 11, www.embedded.com/209600564 “A funny thing happened…” December 2008, pg. 9, www.embedded.com/212200180

In my experience, it's always a mistake to try to anticipate every possible use one should expect, or every possible feature one might need, in a piece of software. I tried that approach once before, long ago. As I was building a huge simulation of a tankmounted rocket launcher, I'd ask my boss “Do you think we should include the effect of (mass imbalance, thrust misalignment, rocket asymmetry, etc., etc.)?” His answer was always yes. I put all those effects in, with added parameters to define the additional effects. During the life of the program, not one simulation run was ever made with anything but zeroes in those parameters.
I learned my lesson then. Today, I'm a fervent followed of the P.J. Plauger dictum: “Never put off until tomorrow what you can put off forever.” These days, though I try to allow enough flexibility and modularity in my software to meet expected needs, I don't add features that no one has asked for.
The Matrix class I'm going to leave with you this month is a working, serviceable class, that you can use with confidence for quite a number of problems. Over time, I'm sure that both you and I will come across features we wish it had. We can add them then. The structure of the class–the relationship between the mathematical concept of a matrix and the computer implementation of it–is in place. The rest is simply a collection of member functions. If we've done our job correctly, we can always add more member functions without breaking the class itself.
In the remainder of this column, We won't have room for a lot of my usual philosophical ramblings, or even detailed explanations. Instead, I'm going to focus on describing the changes/additions to the files SimpleMat.cpp and Matrices.cpp . The code for the whole collection, including the regression test drivers, will be posted on the web site.
SimpleMat and friends
When I was developing class Vector , I worried about the issue of efficiency vs. ease of use. I wanted to be able to manipulate the vectors using overloaded operators, but sometimes I also wanted efficiency. My solution was to split the functionality into two files and two representations of a vector. The first, which I called SimpleVec , used the same approach we used back in my Fortran days: Vectors are treated as simple arrays of doubles. The user must pass the dimensions of the vector to the functions, and face the consequences if he gets them wrong. I eased the pain where I could by using default values for the dimensions. The only inefficiency in the SimpleVec functions is the overhead of the function call itself–not much of a penalty, these days.
In class Vector , on the other hand, I wanted the convenience and robustness one gets by carrying the vector size along as part of the class. You can't give the wrong dimensions because they're part of the object. In class Vector , I made heavy use of overloaded operators to let me write expressions like
y = (a+2*b)/(a*b); // last term is a dot product
The end result is, I get the efficiency of Fortranlike functions when I need it, but the ease of use of a class with overloaded operators, when I prefer that.
There's a lot of similarity between vector and matrix arithmetic. In fact, from one point of view a matrix can be considered to be an array of vectors. The big difference is that the elements of a vector map nicely onto the computer implementation we call an array. The matrix doesn't map so easily. At least, it doesn't if we want to have variable array dimensions. If you want to do matrix arithmetic in C/C++, it's necessary to store the data in a onedimensional array, just as for vectors. Then you use simple rules for indexing into the array. To the computer, the “meat” of both data types look like simple onedimensional arrays.
When I first began the class Matrix , I hoped to be able to leverage this similarity by using calls to the SimpleVec functions. For example, to add two vectors, you simply add their elements, two by two. The same is true of matrices. It seemed a shame to repeat the same operation twice, when a call to vAdd( ) works for both cases.
Other operations, however, are completely different for the two cases. There are two multiplication operators for vectors, only one for matrices, and the algorithms for multiplying them are different. Vector math has no function of division; matrix math does (sort of). In the end, I found it much cleaner to separate the primitive operations into two files, SimpleVec and SimpleMat. So what if a couple of functions look alike? The code we're talking about is tiny either way.
This month, I've made a few additions to SimpleMat . First, I noticed that I'd left out comparison functions–functions that were in SimpleVec . These are the functions mIsEqual( ) and mIsUnequal( ) . They work in the way you might expect: two matrices are equal if and only if all their elements are equal.
The usual caveat concerning comparison of floating point numbers apply. For those who don't know (yes, there are plenty), many numbers can't be represented exactly in floating point, so when you do arithmetic with them, you may be surprised that the result isn't quite what you expect. Thus:
(1)
because 1/3 cannot be represented exactly in floating point format. Indeed, no fraction can, except negative powers of two. Comparing real numbers for equality is only reliable if there is no intervening arithmetic to cause roundoff error.
In the case of the function mIsEqual , I still lean back on the vector function. The code is simply:
// Test for equalitybool mIsEqual(const double a[ ], const double b[ ], const size_t m, const size_t n){ return vIsEqual(a, b, m*n);}
The vector class also had tests for >, < , etc. They were based on the absolute value (norm) of the vector. We don't attempt similar tests for matrices, because the concept of one matrix being numerically (not dimensionally) larger than another is so fuzzy. I suppose we could use the determinant as a measure of size, but I can't think of a good reason to do so. Here, I'm going to invoke the Plauger doctrine.
I also noticed that I'd left out the transpose operator. This is important, because we use the operation a lot. To transpose a matrix, you simply exchange rows with columns, and vice versa. Mathematically, b _{ji } =a _{ij } , where i and j are the row and column indices, respectively. The code reflects this relationship, while taking into account the conversion from matrix indices to simple array index. The key line in the code is:
b[j*m+i] = a[i*n+j];
where m and n are the row and column dimensions of A .
Be warned that you can't transpose a matrix back on itself. That is, you can't write:
mTrans(a, a, m, n);
Doing so scrambles the matrix that you're working on.
Consider, however, what happens if the matrix is square. The diagonal elements don't move during the tranposition, and the offdiagonals only move twoby two, exchanging them symmetrically across the diagonal. It's definitely possible to write a function that will transpose a square matrix in place, though the indexing is a bit tricky.
Everything goes all to heck, though, if the matrix isn't square. The process goes from exchanging elements, two by two, to something more akin to rotating a Rubik's cube. It can be done, but the logic is just horrible, and not at all recommended. In the end, I chose not to provide a function that only works for square matrices. I may revisit this decision in the future; we'll see.
Transpose multiplies
On the other hand, operations like A ^{T} B are extremely common in matrix math. It seems a shame to create a new temporary every time you need to perform the operation. Furthermore, it's not really necessary to transpose either matrix, as long as you know what you're doing. After all, the matrix elements are stored in some order, and you know what that order is. So instead of moving the elements around, we only need to access them in a different order.
The authors of the old IBM Scientific Subroutine Package thought so, too. They supplied a subroutine (I think it was TMATX) to do the job. I've included that function in SimpleMat.cpp , under the name mTmult( ) .
For about 40 years I managed to get by without feeling the need for the corresponding operation, AB ^{T} . Since memory is cheap these days, I went ahead and wrote mMultT( ) as well, (thereby violating Plauger's rule).
Be warned: the ground rules for these functions are tricky, tricky, tricky. When you multiply two matrices, the inner dimensions have to match. That is, the number of columns in the first matrix (say A ) must match the number of rows in the second (B ). The result has the rows of A , and the columns of B . In short:
(2)
The calling list for function mMult( ) has the form:
mMult(A, B, C, L, m, n);
Since we know that B must have m rows, we saw/see no reason to repeat it. The dimensions of the product, C will be L Ã— n .
That's all well and good for the simple matrix product. But what if I'm going to be transposing A before the multiply? The call to mTmult( ) looks the same:
mTmult(A, B, C, L, m, n);
where, again, L and m represent the number of rows and columns in A . But this time, A is going to get conceptually tranposed to an (m Ã— L ) matrix before the multiply. So the inner dimension that B must match–that is, the row dimension of B –must be L , not m . And it's a dimension not even mentioned in the calling list. The number of rows in the transposed A is m , not L , and so the dimensions of C must come out to be (m Ã— n ).
It gets even worse for the other case. In:
mMultT(A, B, C, L, m, n);
It's the B matrix that's transposed. The dimensions L and m still refer to the dimensions of A , and the number of rows in B ^{T} must still be equal to m. But that's the column dimension of B , not the row dimension. Therefore n must now be the row dimension.
If you find all this too confusing for words, not to worry. In many cases we're dealing with square matrices anyhow, so the dimensions are all equal. If they're not, you can always choose the simpler approach: first transpose one of the matrices, then multiply them.
When we get to the matrix class, it will be simpler yet. First of all, there is no T* or *T operator in C++, so there is nothing we can overload. The only operators I can think of that can appear on either side of an argument are ++ and — , but using them would be just too bizarre for words. So if you want to use the overloaded operators in class Matrix , you're going to have to settle for:
C=T(A)*B;
or:
C=A*T(B);
For the purist, I'm including functions in the matrix class called Tmult( ) and multT( ) . They work the same as the functions in SimpleMat , except that you don't have to deal with the dimensions. Those dimensions are carried along as part of the objects, and the functions will check them and declare an error if they don't match up.
The test driver
To test SimpleMat.cpp , I wrote a test driver, imaginatively called TestSimpleMat . This driver is called, in turn, by the overall regression test program, TestMatrices .
The philosophy behind these test programs bears repeating.
Each function in the classes Vectors and Matrices is, itself, relatively simple. In the old days, I would have written a test driver for each function, which asked me to input some test values of the inputs. Then it would output the results, and I'd have to verify that these results were correct. To do that requires a hand check. If I had to test the same function again, I'd probably have to repeat the hand check because I couldn't remember the values I'd used the first time. It does you absolutely no good to run the test without verifying the answer. To do so would violate the criterion voiced by my pal, Jim Adams: “When testing, never ask a question unless you already know the answer.” Just saying, “Yep, that looks right” is not an option.
These days, I incorporate both the test inputs and expected outputs directly in the test driver. I use asserts to verify that the answers are correct. This results in what I call “Silent Testing.” If all the tests pass, no asserts are violated, no exceptions are thrown, and the test program returns with no output at all. In this case, no news is good news.
When I was developing SimpleMat , I deviated from this path for the sake of expedience, reverting instead to the “Give me a vector to test,” highly verbose sort of exchange. That's why I didn't post the code on Embedded.com; I didn't want you to see just how badly I can write code when time is tight.
Since then, I've spent a lot of time cleaning up TestSimpleMat and making it conform to the Silent Test sort of behavior. That code is posted, along with all the other files, on our Embedded.com source code library (Editor's Note: you now have to sign in to download code) .
I'll only mention in passing that, some months back, I promised (threatened?) to reduce the number of files involved by combining SimpleVec , SimpleMat , Vectors , and Matrices into a single file and its corresponding header file. Fortunately, I took a cold shower and forgot that idea. As with everyone else, not all my ideas are good ones.
Matrix math functions
In the end, adding the math operations to the matrix class is the easiest part of the job. I had similar functions in Vectors.cpp , so setting up the ones in the matrix class was pretty easy. You'll find the code on our web site. In Listing 1 , I show the typical usage of each function.
Understand, these functions by no means exhaust the possible functions we can provide in this class. As I mentioned earlier, the math for matrices is a rich one indeed, and even moreso when you combine matrices with vectors. We could go on adding functions for the rest of eternity.
Listing 1 Matrix math functions.
B = T(A); // transpose a matrixA.Negate( ); // replace A by its negativeB = A; // unary minusA += B; // add in placeA += 1.0; // add scalar matrix in placeA = B; // subtract in placeA = 1.0; // subtract scalar matrix in placeC = A+B; // matrix addC = A+s; // add matrix and scalar (square only)C = AB; // matrix subtractC = As; // subtract scalar from matrix (square only)C = A*B; // matrix productC *= 2.0; // scale matrix in placeA = 2.0*A; // multiply matrix by scalarA = A*2.0; // same resultA.Invert( ); // matrix inverseC = B/A; // multiply matrix by inverse: C=A^{1} By = A*x; // multiply vector by matrixy = x/A; // multiply vector by inverse: y=A^{1} x
I can think of many more operations I wish I had. I'd like to be able to construct a matrix from a set of vectors. I'd like to be able to extract rows and columns from a matrix. I'd like to be able to build a matrix from a set of smaller ones, and I'd like to be able to subdivide them again. I'd like to create a diagonal matrix whose diagonal elements are defined by a vector. I'd like to be able to extract the diagonal into a vector.
I want at least two nontrivial functions that I've not written so far. One is to compute the determinant of a matrix. That determinant is computed within the matrix inverter, mInv( ) , and returned from that function. But I have no “compute determinant function” otherwise.
The other function that's often essential for work in control theory, is computing the eigenvalues and eigenvectors of a square matrix (if you don't know what these are, you probably don't need the function). The problem with that one is that the results are often complex numbers. So far, I've not included complex arithmetic in these classes, so that one is going to be a special case.
So there's no doubt that the class Matrix needs more functions. I'll be adding them to the class over time and will show them to you as I do so. In the meantime, however, the class as I've presented it here, and its related files, are nice, serviceable classes in their own right. So I'm going to invoke Plauger's dictum and leave the rest until someone asks for them. Use the libraries in good health.
Jack Crenshaw is a systems engineer and the author of Math Toolkit for RealTime Programming. He holds a PhD in physics from Auburn University. Email him at jcrens@earthlink.net. For more information about Jack click here