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]

Solution 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

Solution 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.

Solution 3

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)];

Solution 4

For the records, this is the same function, but with one-based indexing, and in Julia:

function iuppert(k::Integer,n::Integer)
  i = n - 1 - floor(Int,sqrt(-8*k + 4*n*(n-1) + 1)/2 - 0.5)
  j = k + i + ( (n-i+1)*(n-i) - n*(n-1) )÷2
  return i, j
end

Solution 5

In python 2:

def k2ij(k, n):
    rows = 0
    for t, cols in enumerate(xrange(n - 1, -1, -1)):
        rows += cols
        if k in xrange(rows):
            return (t, n - (rows - k))
    return None

View more solutions

Comments

[n. m.]

Write a fomula for elements in the last column. To make it easier write a formula that computes the linear index from a row number (the column number is fixed), then invert it. Proceed to a general formula from there.

[Robert Crovella]

Note that the solution methods presented here can also be used to list the combinations of N things taken 2 at a time (without repetition), without the need for any iteration/recursion.

[Hope]

Same question without the offset was answered here

[linello]

Perfect! This helped me to reduce the number of variables in a linear program!

[Squidly]

There's an extra bracket in i = ...

[Imran]

can you explain formula of K / how it drives?

[Hugo Maxwell]

You need a couple more parenthesis to make sure to avoid rounding a temporary value when diving by 2: k = ((n*(n-1))/2) - ((n-i)*((n-i)-1))/2 + j - i - 1

[jay s]

thanks, this really helped me understand the solution

[Elarion]

Note that the equations are for zero-indexed variables...

[cmo]

you need import numpy as np in there for the assert np.triu_indices lines.

[seralouk]

what about an ij2k function that does the inverse ? @smac89

[seralouk]

xrange needs to be replaced by range

[smac89]

@seralouk in python2, xrange is more efficient than range

[seralouk]

thanks. is there any python function from numpy or scipy that does this as well?