Triangular Matrix
published: Tue, 12-Dec-2006 | updated: Tue, 12-Dec-2006
Every time I have to do this, I forget and have to look it up again. The last time, I thought I'd remembered how to do it and then wasted the next 10 minutes because I'd actually remembered it the wrong way round (well, I was tired, and an article had to be finished). The "this" I'm talking about? Implementing triangular matrices.
A lot of times when we need to implement a square matrix, the cells in the matrix are symmetric about the main diagonal, that is, the data in cell (x, y) is the same as the data in the cell (y, x) for all x and y. In my recent example I had a matrix for the distances between a set of cities — it was a Symmetric Travelling Salesman Problem — and the distance between city A and city B is assumed to be the same as the distance between B and A.
In this kind of situation it's a waste of space to store the duplicate information. We should just store the data for the main diagonal and for the upper (or lower) triangle and then encapsulate the data and its access into a normal-looking matrix class.
If we store the data in the lower triangle then all row numbers are greater than or equal to the column numbers. So if we're asked for the data in cell (1,2) all we need to do is reverse the coordinates and return the value for (2,1).
But how do we efficiently store the values for the lower triangle? (And my asking this question means that it is the lower triangle we are going to use.) Take a look at this triangular matrix where I've numbered the cells starting at the top left and counted down the rows (the "unused" cells in the matrix have periods):
0 . . . . 1 2 . . . 3 4 5 . . 6 7 8 9 . 10 11 12 13 14
This ladies and gentlemen is known as Pascal's Triangle (well, at least the zero-based version of it). There is one filled cell on the first row, three in the first two rows, six in the first three, 10 in the first four rows. If you know about Pascal's Triangle, you'll know the formula for calculating the number of items in an n-row triangle, but let's derive it for fun.
The number of cells is 1 + 2 +...+ n-1 + n Or, in reverse n + n-1 +...+ 2 + 1 Add these two up n+1 + n+1 +...+ n+1 + n+1 This is twice what we want, so divide by two (n+1 + n+1 +...+ n+1 + n+1) / 2 There are n items being added up in the brackets, so we can get rid of the ellipsis n * (n+1) / 2
So there are n*(n+1)/2 items in an n-row Pascal's Triangle. (Notice something: this is always a whole number since either n or n+1 is going to be even and so their product will be divisible by 2.)
(Aside: the great mathematician Gauss, as a little boy, was in class one day when the teacher decided he wanted some peace and quiet for a while. So he set the class a problem: sum all the numbers from 1 to 100, and figured he'd have a good 10 minutes rest. After a few seconds, young Gauss put up his hand to say that he'd finished and the answer was 5050. Gauss had worked out the process I'd gone through above and calculated that the answer must be 100*101/2.)
Where were we? Oh yes. From the figure above you can see pretty easily that we can simulate a 5 row triangular matrix with a 15 element array, and the numbers of the cells indicate the element in the array. So element 7 of the array would marry up to cell (3,1), and of course (1,3). But how do we calculate the element number given the cell address? Using the formula we just derived of course.
Suppose we want to find the element number for cell (x,y) where x >= y (all I'm saying here is the row number is greater than or equal to the column number, if not just swap them over). How many elements cover all the rows less than x? Simple: x*(x+1)/2. (Try it out.) We can now easily count along the row until we reach the column we need, or simpler, just add the column number.
So, the element in the array for cell (x,y), where x >= y, is (x*(x+1)/2 + y). (Test this with cell (3,1) => element 7.)
After that minor bit of algebra, the code is pretty boring since all that needs to be done is encapsulate this formula. Here it is in Delphi (this had be be first in honor of Pascal):
type
TSymmetricMatrix = class
strict private
FMatrix : array of double;
FLength : integer;
function getIndex(row, column : integer) : integer;
function getValue(row, column : integer) : double;
procedure setValue(row, column : integer; aValue : double);
public
constructor Create(aRowCount : integer);
property Value[row, column : integer] : double
read getValue write setValue; default;
end;
constructor TSymmetricMatrix.Create(aRowCount : integer);
begin
inherited Create;
FLength := aRowCount;
SetLength(FMatrix, (aLength * (aLength + 1)) div 2);
end;
function TSymmetricMatrix.getIndex(row, column : integer) : integer;
var
temp : integer;
begin
if not ((0 <= row) and (row < FLength)) then
raise Exception.Create('Symmetric matrix row index out of bounds');
if not ((0 <= column) and (column < FLength)) then
raise Exception.Create('Symmetric matrix column index out of bounds');
if (row < column) then begin
temp := row;
row := column;
column := temp;
end;
result := ((row * (row + 1)) div 2) + column;
end;
function TSymmetricMatrix.getValue(row, column : integer) : double;
begin
result := FMatrix[getIndex(row, column)];
end;
procedure TSymmetricMatrix.setValue(row, column : integer; aValue : double);
begin
FMatrix[getIndex(row, column)] := aValue;
end;
And here it is in C# using generics.
public class SymmetricMatrix<T> {
T[] matrix;
int length;
public SymmetricMatrix(int rowCount) {
this.matrix = new T[(rowCount * (rowCount + 1) / 2)];
this.length = rowCount;
}
public int getIndex(int row, int column) {
if (row < 0 || row >= length)
throw new ArgumentOutOfRangeException("Symmetric matrix row index out of bounds");
if (column < 0 || column >= length)
throw new ArgumentOutOfRangeException("Symmetric matrix column index out of bounds");
if (row < column) {
int temp = row;
row = column;
column = temp;
}
return ((row * (row + 1)) / 2) + column;
}
public T this[int row, int column] {
get {
return matrix[getIndex(row, column)];
}
set {
matrix[getIndex(row, column)] = value;
}
}
}
There, now I have it online the next time...
