algorithm for index numbers of triangular matrix coefficients

Question

I think this must be simple but I can't get it right...

I have an MxM triangular matrix, the coefficients of which are stored in a vector, row by row. For example:

M =   [ m00 m01 m02 m03 ]        [     m11 m12 m13 ]       [         m22 m23 ]       [             m33 ] 

is stored as

coef[ m00 m01 m02 m03 m11 m12 m13 m22 m23 m33 ] 

Now I'm looking for a non-recursive algorithm that gives me for matrix size M and coefficient array index i

unsigned int row_index(i,M) 

and

unsigned int column_index(i,M) 

of the matrix element that it refers to. So, row_index(9,4) == 3, column_index(7,4) == 2 etc. if the index counting is zero-based.

EDIT: Several replies using an iteration have been given. Does anyone know of algebraic expressions?

Answer 1

One-liners at the end of this reply, explanation follows :-)

The coefficient array has: M elements for the first row (row 0, in your indexing), (M-1) for the second (row 1), and so on, for a total of M+(M-1)+…+1=M(M+1)/2 elements.

It's slightly easier to work from the end, because the coefficient array always has 1 element for the last row (row M-1), 2 elements for the second-last row (row M-2), 3 elements for row M-3, and so on. The last K rows take up the last K(K+1)/2 positions of the coefficient array.

Now suppose you are given an index i in the coefficient array. There are M(M+1)/2-1-i elements at positions greater than i. Call this number i'; you want to find the largest k such that k(k+1)/2 ≤ i'. The form of this problem is the very source of your misery -- as far as I can see, you cannot avoid taking square roots :-)

Well let's do it anyway: k(k+1) ≤ 2i' means (k+1/2)² - 1/4 ≤ 2i', or equivalently k ≤ (sqrt(8i'+1)-1)/2. Let me call the largest such k as K, then there are K rows that are later (out of a total of M rows), so the row_index(i,M) is M-1-K. As for the column index, K(K+1)/2 elements out of the i' are in later rows, so there are j'=i'-K(K+1)/2 later elements in this row (which has K+1 elements in total), so the column index is K-j'. [Or equivalently, this row starts at (K+1)(K+2)/2 from the end, so we just need to subtract the starting position of this row from i: i-[M(M+1)/2-(K+1)(K+2)/2]. It is heartening that both expressions give the same answer!]

(Pseudo-)Code [using ii instead of i' as some languages don't allow variables with names like i' ;-)]:

row_index(i, M):     ii = M(M+1)/2-1-i     K = floor((sqrt(8ii+1)-1)/2)     return M-1-K  column_index(i, M):     ii = M(M+1)/2-1-i     K = floor((sqrt(8ii+1)-1)/2)     return i - M(M+1)/2 + (K+1)(K+2)/2 

Of course you can turn them into one-liners by substituting the expressions for ii and K. You may have to be careful about precision errors, but there are ways of finding the integer square root which do not require floating-point computation. Also, to quote Knuth: "Beware of bugs in the above code; I have only proved it correct, not tried it."

If I may venture to make a further remark: simply keeping all the values in an M×M array will take at most twice the memory, and depending on your problem, the factor of 2 might be insignificant compared to algorithmic improvements, and might be worth trading the possibly expensive computation of the square root for the simpler expressions you'll have.

[Edit: BTW, you can prove that floor((sqrt(8ii+1)-1)/2) = (isqrt(8ii+1)-1)/2 where isqrt(x)=floor(sqrt(x)) is integer square root, and the division is integer division (truncation; the default in C/C++/Java etc.) -- so if you're worried about precision issues you only need to worry about implementing a correct integer square root.]

Answer 2

Here's an algebraic (mostly) solution:

unsigned int row_index( unsigned int i, unsigned int M ){     double m = M;     double row = (-2*m - 1 + sqrt( (4*m*(m+1) - 8*(double)i - 7) )) / -2;     if( row == (double)(int) row ) row -= 1;     return (unsigned int) row; }   unsigned int column_index( unsigned int i, unsigned int M ){     unsigned int row = row_index( i, M);     return  i - M * row + row*(row+1) / 2; } 

EDIT: fixed row_index()