Matrix Mathematics

Mathematics can sometimes be a bit daunting, and one thing I've always found a bit difficult is matrix mathematics. So, I am writing down some notes regarding this particular subject, and the implementation of said subject in C/C++.

Mathematical Notation

Matricies are very abstract in the sense that we can write them down in many different ways. For example a 4×4 identity matrix written in standard mathematical notation:

When we consider standard mathematical notation, the origin is in the top left hand side, and we specify size in the format (rows, columns). Never confuse rows and columns with width and height, as things may get confusing.

Mathematically, here are the coordinates into the matrix, when we consider individual items, R designating the row count, and C designating the column count:

Memory layout

In C when we can access data in whatever order we like.

Row Major: Elements of a row are stored contiguously in memory.
Column Major: Elements of a column are stored contiguously in memory.

In C we also have multi-dimensional arrays to assist with accessing data. These are per the C specification laid out in memory in sequential arrays the size of the last dimension. For example:

float v[2][4];

What this is saying is "2 arrays of 4 float values". And that is the way they are arranged in memory also:

Memory Offset	0	1	2	3	4	5	6	7
Index[row][col]	v[0][0]	v[0][1]	v[0][2]	v[0][3]	v[1][0]	v[1][1]	v[1][2]	v[1][3]

Here is an alternative view (values are memory offsets):

	Col 0	Col 1	Col 2	Col 3
Row 0	0	1	2	3
Row 1	4	5	6	7

The function for row and column major order are not complex:

// If we increase row by 1, the offset will increase by sz (number of elements per row i.e. number of columns)
// If we increase col by 1, the offset will increase by 1
unsigned rowMajorOffset(unsigned row, unsigned col, unsigned sz)
{
  return col + row * sz;
}
 
// If we increase col by 1, the offset will increase by sz (number of elements per column i.e. number of rows)
// If we increase row by 1, the offset will increase by 1
unsigned columnMajorOffset(unsigned row, unsigned col, unsigned sz)
{
  return row + col * sz;
}

It is important you understand the following are equivalent terms:

number of elements per column	⇔	number of rows
number of elements per row	⇔	number of columns

Here we can use the rowMajorLookup function which is equivalent to the standard C array indexing:

// These values are equivalent (in function)
float a = v[0][3];
float b = ((float *)v)[rowMajorLookup(0, 3, 4)];

So, we know the standard mathematical representation of a matrix, and we know how they can be represented in memory, but where does that leave us? Right, we need to make some concrete implementation so we can really see how it all comes together.

Firstly, we have two options to consider, as mentioned earlier:

Row major order

A single row should be sequential in memory: float values[2][4]

	Col 0	Col 1	Col 2	Col 3
Row 0	0	1	2	3
Row 1	4	5	6	7

Column major order

A single column should be sequential in memory float values[4][2]

	Col 0	Col 1	Col 2	Col 3
Row 0	0	2	4	6
Row 1	1	3	5	7

As you can see, it is as simple as swapping (the technical term is transpose) around the dimensions for the array declaration. We can design a simple template class which will handle just about anything we can throw at it.

#include <iostream>
 
// If we increase row by 1, the offset will increase by sz (number of elements per row i.e. number of columns)
// If we increase col by 1, the offset will increase by 1
unsigned rowMajorOffset(unsigned row, unsigned col, unsigned sz)
{
  return col + row * sz;
}
 
// If we increase col by 1, the offset will increase by sz (number of elements per column i.e. number of rows)
// If we increase row by 1, the offset will increase by 1
unsigned columnMajorOffset(unsigned row, unsigned col, unsigned sz)
{
  return row + col * sz;
}
 
template <typename _ValueT, unsigned _R, unsigned _C, bool _ColumnMajor>
class Matrix {
protected:
	enum { ColumnMajor = _ColumnMajor };
	enum { R = _R };
	enum { C = _C };
 
	typedef _ValueT ValueT;
 
	ValueT m_values[C*R];
 
public:
	const ValueT & at (unsigned r, unsigned c) const
	{
		if (ColumnMajor)
			return m_values[columnMajorOffset(r, c, R)];
		else
			return m_values[rowMajorOffset(r, c, C)];
	}
 
	ValueT & at (unsigned r, unsigned c)
	{
		if (ColumnMajor)
			return m_values[columnMajorOffset(r, c, R)];
		else
			return m_values[rowMajorOffset(r, c, C)];
	}
 
	void loadTestPattern ()
	{
		for (unsigned r = 0; r < R; r += 1)
			for (unsigned c = 0; c < C; c += 1)
				at(r, c) = (r+1) * 1000 + (c+1);
	}
 
	void debug ()
	{
		using namespace std;
 
		if (ColumnMajor)
			cout << "Column-Major Matrix " << "(" << R << "," << C << ")" << " @ " << this << endl;
		else
			cout << "Row-Major Matrix " << "(" << R << "," << C << ")" << " @ " << this << endl;
 
		cout << "Memory Offset: ";
		for (unsigned i = 0; i < (R*C); i += 1)
			cout << i << "    ";
		cout << endl;
 
		cout << "       Values: ";	
		for (unsigned i = 0; i < (R*C); i += 1)
			cout << m_values[i] << " ";
		cout << endl;
 
		cout << "Standard Mathematical Notation:" << endl;
		cout << "      ";
		for (unsigned c = 0; c < C; c += 1)
			cout << "Col " << c << " ";
		cout << endl;
 
		for (unsigned r = 0; r < R; r += 1) {
			cout << "Row " << r << " ";
			for (unsigned c = 0; c < C; c += 1)
				cout << at(r, c) << "  ";
			cout << endl;
		}
		cout << endl;
	}
 
	Matrix<ValueT, R, C, !ColumnMajor> transposeStorage () const
	{
		Matrix<ValueT, R, C, !ColumnMajor> result;
 
		for (unsigned r = 0; r < R; r += 1)
			for (unsigned c = 0; c < C; c += 1)
				result.at(r, c) = at(r, c);
 
		return result;
	}
 
	Matrix<ValueT, C, R, !ColumnMajor> transposeMatrix () const
	{
		Matrix<ValueT, C, R, !ColumnMajor> result;
 
		memcpy(&result.at(0,0), m_values, sizeof(m_values));
 
		return result;
	}
};
 
int main (int argc, char * const argv[]) {
	Matrix<float, 4, 2, false> rowMajorMatrix;
	Matrix<float, 4, 2, true> columnMajorMatrix;
 
	rowMajorMatrix.loadTestPattern();
	rowMajorMatrix.debug();
 
	columnMajorMatrix.loadTestPattern();
	columnMajorMatrix.debug();
 
	rowMajorMatrix = columnMajorMatrix.transposeStorage();
	rowMajorMatrix.debug();
 
	Matrix<float, 2, 4, false> transposedMatrix = columnMajorMatrix.transposeMatrix();
	transposedMatrix.debug();
 
	return 0;
}

The code is available as an Xcode project.

Here is the result after running this code

Row-Major Matrix (4,2) @ 0xbffff6a0
Memory Offset: 0    1    2    3    4    5    6    7    
       Values: 1001 1002 2001 2002 3001 3002 4001 4002 
Standard Mathematical Notation:
      Col 0 Col 1 
Row 0 1001  1002  
Row 1 2001  2002  
Row 2 3001  3002  
Row 3 4001  4002  

Column-Major Matrix (4,2) @ 0xbffff680
Memory Offset: 0    1    2    3    4    5    6    7    
       Values: 1001 2001 3001 4001 1002 2002 3002 4002 
Standard Mathematical Notation:
      Col 0 Col 1 
Row 0 1001  1002  
Row 1 2001  2002  
Row 2 3001  3002  
Row 3 4001  4002  

Row-Major Matrix (4,2) @ 0xbffff6a0
Memory Offset: 0    1    2    3    4    5    6    7    
       Values: 1001 1002 2001 2002 3001 3002 4001 4002 
Standard Mathematical Notation:
      Col 0 Col 1 
Row 0 1001  1002  
Row 1 2001  2002  
Row 2 3001  3002  
Row 3 4001  4002  

Row-Major Matrix (2,4) @ 0xbffff660
Memory Offset: 0    1    2    3    4    5    6    7    
       Values: 1001 2001 3001 4001 1002 2002 3002 4002 
Standard Mathematical Notation:
      Col 0 Col 1 Col 2 Col 3 
Row 0 1001  2001  3001  4001  
Row 1 1002  2002  3002  4002

As you can see, it doesn't matter how we store the data in memory (depending on what you are doing, it may be more or less convenient to use either way), we can still retrieve the data in standard mathematical form.

Conclusion

Ultimately, depending on what other libraries you want to integrate your code with, you may need to choose a specific format. Other reasons, such as performance, may dictate the best layout for your data - i.e. vectors of data are best laid out sequentially for quick and easy access. SSE/3DNow/MMX may also have strict requirements on data layout and alignment for parallelization and optimization of specific operations.

It is also worthwhile to point out, that at least in my case, it was confusing to think of matrices having "width" and "height". In order to consider a matrix with width and height, we actually need to rotate the matrix to get a co-ordinate system that makes sense, i.e. x, y in the origin.

I personally find the transpose functions fascinating. For example, simply by reinterpreting the data from row major to column major we can transpose the data.