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()

Answer 3

The explanation of ShreevatsaR is excellent and helped me solve my problem. However, the explanation and code provided for the column index give a slightly different answer than what the question asks for.

To reiterate, there are j' = i' - K(K+1)/2 elements in the row after i. But the row, like every other row, has M elements. Hence the (zero-based) column index is y = M - 1 - j'.

The corresponding pseudo code is:

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

The answer given by ShreevatsaR, K - j', is the column index when starting to count (with zero) from the diagonal of the matrix. Hence, his calculation gives column_index(7,4) = 0 and not, as specified in the question, column_index(7,4) = 2.

Answer 4

There may be a clever one liner for these, but (minus any error checking):

unsigned int row_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    for( unsigned int x = delta; x < i; x += delta-- ){
        row++;
    }
    return row;
}

unsigned int column_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    unsigned int x;
    for( x = delta; x < i; x += delta-- ){
        row++;
    }
    return M + i - x - 1;
}

Answer 5

It has to be that

i == col + row*(M-1)-row*(row-1)/2

So one way to find col and row is to iterate over possible values of row:

for(row = 0; row < M; row++){
  col = i - row*(M-1)-row*(row-1)/2
  if (row <= col < M) return (row,column);
}

This is at least non-recursive, I don't know if you can do this without iteration.

As can be seen from this and other answers, you pretty much have to calculate row in order to calculate column, so it could be smart to do both in one function.

Answer 6

I thought a little bit and I got the following result. Note that you get both row and columns on one shot.

Assumptions: Rows start a 0. Columns start at 0. Index start at 0

Notation

N = matrix size (was M in the original problem)

m = Index of the element

The psuedo code is

function ind2subTriu(m,N)
{
  d = 0;
  i = -1;
  while d < m
  {
    i = i + 1
    d = i*(N-1) - i*(i-1)/2
  }
  i0 = i-1;
  j0 = m - i0*(N-1) + i0*(i0-1)/2 + i0 + 1;
  return i0,j0
}

And some octave/matlab code

function [i0 j0]= ind2subTriu(m,N)
 I = 0:N-2;
 d = I*(N-1)-I.*(I-1)/2;
 i0 = I(find (d < m,1,'last'));
 j0 = m - d(i0+1) + i0 + 1;

What do you think?

As of December 2011 really good code to do this was added to GNU/Octave. Potentially they will extend sub2ind and ind2sub. The code can be found for the moment as private functions ind2sub_tril and sub2ind_tril

Answer 7

Took me some time to understand what you needed! :)

unsigned int row_index(int i, int m)
{
    int iCurrentRow = 0;
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return iCurrentRow;

        iCurrentRow ++;
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}

Sorry forgot the other function.....

unsigned int column_index(int i, int m)
{
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return m - (iTotalItems - i);
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}