Linear index upper triangular matrix Linear index upper triangular matrix numpy numpy

Linear index upper triangular matrix


The equations going from linear index to (i,j) index are

i = n - 2 - floor(sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5)j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2

The inverse operation, from (i,j) index to linear index is

k = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1

Verify in Python with:

from numpy import triu_indices, sqrtn = 10for k in range(n*(n-1)/2):    i = n - 2 - int(sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5)    j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2    assert np.triu_indices(n, k=1)[0][k] == i    assert np.triu_indices(n, k=1)[1][k] == jfor i in range(n):    for j in range(i+1, n):        k = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1        assert triu_indices(n, k=1)[0][k] == i        assert triu_indices(n, k=1)[1][k] == j


First, let's renumber a[k] in opposite order. We'll get:

0  a9  a8  a7  a60   0  a5  a4  a30   0   0  a2  a10   0   0   0  a00   0   0   0   0

Then k2ij(k, n) will become k2ij(n - k, n).

Now, the question is, how to calculate k2ij(k, n) in this new matrix. The sequence 0, 2, 5, 9 (indices of diagonal elements) corresponds to triangular numbers (after subtracting 1): a[n - i, n + 1 - i] = Ti - 1. Ti = i * (i + 1)/2, so if we know Ti, it's easy to solve this equation and get i (see formula in the linked wiki article, section "Triangular roots and tests for triangular numbers"). If k + 1 is not exactly a triangular number, the formula will still give you the useful result: after rounding it down, you'll get the highest value of i, for which Ti <= k, this value of i corresponds to the row index (counting from bottom), in which a[k] is located. To get the column (counting from right), you should simply calculate the value of Ti and subtract it: j = k + 1 - Ti. To be clear, these are not exacly i and j from your problem, you need to "flip" them.

I didn't write the exact formula, but I hope that you got the idea, and it will now be trivial to find it after performing some boring but simple calculations.


The following is an implimentation in matlab, which can be easily transferred to another language, like C++. Here, we suppose the matrix has size m*m, ind is the index in the linear array. The only thing different is that here, we count the lower triangular part of the matrix column by column, which is analogus to your case (counting the upper triangular part row by row).

function z= ind2lTra (ind, m)  rvLinear = (m*(m-1))/2-ind;  k = floor( (sqrt(1+8*rvLinear)-1)/2 );  j= rvLinear - k*(k+1)/2;  z=[m-j, m-(k+1)];