Linear index upper triangular matrix

Question

If I have the upper triangular portion of a matrix, offset above the diagonal, stored as a linear array, how can the (i,j) indices of a matrix element be extracted from the linear index of the array?

For example, the linear array [a0, a1, a2, a3, a4, a5, a6, a7, a8, a9 is storage for the matrix

0  a0  a1  a2  a3 0   0  a4  a5  a6 0   0   0  a7  a8 0   0   0   0  a9 0   0   0   0   0 

And we want to know the (i,j) index in the array corresponding to an offset in the linear matrix, without recursion.

A suitable result, k2ij(int k, int n) -> (int, int) would satisfy, for example

k2ij(k=0, n=5) = (0, 1) k2ij(k=1, n=5) = (0, 2) k2ij(k=2, n=5) = (0, 3) k2ij(k=3, n=5) = (0, 4) k2ij(k=4, n=5) = (1, 2) k2ij(k=5, n=5) = (1, 3)  [etc] 

Answer 1

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, sqrt n = 10 for 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] == j  for 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 

Answer 2

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

0  a9  a8  a7  a6 0   0  a5  a4  a3 0   0   0  a2  a1 0   0   0   0  a0 0   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.