Linear Index Upper Triangular Matrix

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, 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

Linear index for a diagonal run of an upper triangular matrix

Thanks to @loopy-walt's observation, we have an answer!
Using the result from Linear index upper triangular matrix, a transformation of the result

(i, j) |-> (j-i-1, j)

Gives the expected outcome.

Here is a C++ implementation.

#include<tuple>
#include<cmath>

// Linear indexing of the upper triangle, row by row
std::tuple<size_t, size_t> k2ij(size_t n, size_t k){
  size_t i = n - 2 - (size_t)std::floor(std::sqrt(4*n*(n-1) - (8*k) -7)/2.0 - 0.5);
  size_t j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2;
  return {i,j};
}

// Linear indexing of the upper triangle, diagonal by diagonal
std::tuple<size_t, size_t> d2ij(size_t n, size_t d){
  const auto [i, j] = k2ij(n, d);
  return {j-i-1, j}; // Conversion from row by row to diag by diag
}

#include<iostream>
#include<set>
int main(int argc, char** argv) {

  size_t n = 4;
  size_t top = n*(n-1)/2;

  for(size_t d=0; d<top; ++d){
    const auto [i,j] = d2ij(n, d);
    std::cout << "d2ij(" << n << ", " << d << ") = (" << i << ", " << j << ")" << std::endl;
  }

  return 0;
}

Producing

d2ij(4, 0) = (0, 1)
d2ij(4, 1) = (1, 2)
d2ij(4, 2) = (2, 3)
d2ij(4, 3) = (0, 2)
d2ij(4, 4) = (1, 3)
d2ij(4, 5) = (0, 3)

Note: if someone wishes the form f(d) instead, a lambda can be used to capture the dimension 'n'

auto f = [n](size_t d){return d2ij(n, d);};
const auto [i,j] = f(5);

Thanks to everybody that took the time to read and help!

algorithm for index numbers of triangular matrix coefficients

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

Giving an element from lower/upper triangular matrix

My solution might be equivalent to your’s, I haven’t checked:

index = N * i - ((i - 1) * i) / 2 + (j - i)

Here’s a complete Python test for it. I used Python because Numpy has triu_indices, which gives the upper-triangular indexes.

import numpy as np

def mksquare(N):
    """Make a square N by N matrix containing 0 .. N*N-1"""
    return np.arange(N * N).reshape(N, N)

def mkinds(N):
    """Return all triu indexes for N by N matrix"""
    return [(i,j) for i in range(N) for j in range(N) if i <= j]

def ij2linear(i, j, N):
    """Convert (i,j) 2D index to linear triu index for N by N array"""
    return N * i - ((i - 1) * i) // 2 + (j - i)

def test(N):
    """Make sure my `mkinds` works for given N"""
    arr = mksquare(N)
    vec = arr[np.triu_indices(N)]

    inds = mkinds(N)
    expected = [arr[i, j] for (i, j) in inds]

    actual = [vec[ij2linear(i, j, N)] for (i, j) in inds]

    return np.all(np.equal(actual, expected))

"""Run `test` for a bunch of `N`s and make sure they're all right"""
print(all(map(test, range(2, 20))))
# prints True br>

Worth a blog post explaining how to arrive at this conclusion, but this’ll do for now .

How to convert triangular matrix indexes in to row, column coordinates?

Not optimized at all :

int j = idx;
int i = 1;

while(j > i) {
    j -= i++;
}

Optimized :

int i = std::ceil(std::sqrt(2 * idx + 0.25) - 0.5);
int j = idx - (i-1) * i / 2;

And here is the demonstration:

You're looking for i such that :

sumRange(1, i-1) < idx && idx <= sumRange(1, i)

when sumRange(min, max) sum integers between min and max, both inxluded.
But since you know that :

sumRange(1, i) = i * (i + 1) / 2

Then you have :

idx <= i * (i+1) / 2
=> 2 * idx <= i * (i+1)
=> 2 * idx <= i² + i + 1/4 - 1/4
=> 2 * idx + 1/4 <= (i + 1/2)²
=> sqrt(2 * idx + 1/4) - 1/2 <= i

Map upper triangular matrix on vector skipping the diagonal

Be careful. The formula you found from that link, index = x + (y+1)*y/2, includes the diagonal entries, and is for a lower triangular matrix, which I gather is not what you want. The exact formula you are looking for is actually index = x + ((y-1)y)/2
(see: https://math.stackexchange.com/questions/646117/how-to-find-a-function-mapping-matrix-indices).

Again, watch out, this formula I gave you assumes your indices: x,y, are zero-based. Your MiniZinc code is using indices i,j that start from 1 (1 <= i <= m), 1 <= j <= m)). For indices that start from 1, the formula is T(i,j) = i + ((j-2)(j-1))/2. So your code should look like:

constraint
   forall ( i in 1..m-1 , j in i+1..m )
      ((distances[(i + ((j-2)*(j-1)) div 2]) >= ...

Note that (j-2)(j-1) will always be a multiple of 2, so we can just use integer division with divisor 2 (no need to worry about converting to/from floats).

---

The above assumes you are using a square m*m matrix.

To generalise to a M*N rectangular matrix, one formula could be:

where 0 <= i < M, 0<= j < N [If you again, need your indices to start from 1, replace i with i-1 and j with j-1 in the above formula]. This touches all of cells of an upper triangular matrix as well as the 'extra block on the side' of the square that occurs when N > M. That is, it touches all cells (i,j) such that i < j for 0 <= i < M, 0 <= j < N.

---

Full code:

% original: https://github.com/MiniZinc/minizinc-benchmarks/blob/master/golomb/golomb.mzn

include "alldifferent.mzn";

int: m;
int: n = m*m;
array[1..m] of var 0..n: mark;
array[1..(m*(m-1)) div 2] of var 0..n: differences;

constraint mark[1] = 0;
constraint forall ( i in 1..m-1 ) ( mark[i] < mark[i+1] );
constraint alldifferent(differences);
constraint forall (i,j in 1..m where j > i)
    (differences[i + ((j-1)*(j-2)) div 2] = mark[j] - mark[i]);
constraint forall (i,j in 1..m where j > i)
    (differences[i + ((j-1)*(j-2)) div 2] >= (floor(int2float(( j-i )*( j-i+1 )) / int2float(2))));
constraint differences[1] < differences[(m*(m-1)) div 2];

solve :: int_search(mark, input_order, indomain, complete)
    minimize mark[m];

output ["golomb ", show(mark), "\n"];

Lower triangular version (take previous code and swap i and j where necessary):

% original: https://github.com/MiniZinc/minizinc-benchmarks/blob/master/golomb/golomb.mzn

include "alldifferent.mzn";

int: m;
int: n = m*m;
array[1..m] of var 0..n: mark;
array[1..(m*(m-1)) div 2] of var 0..n: differences;

constraint mark[1] = 0;
constraint forall ( i in 1..m-1 ) ( mark[i] < mark[i+1] );
constraint alldifferent(differences);
constraint forall (i,j in 1..m where i > j)
    (differences[j + ((i-1)*(i-2)) div 2] = mark[i] - mark[j]);
constraint forall (i,j in 1..m where i > j)
    (differences[j + ((i-1)*(i-2)) div 2] >= (floor(int2float(( i-j )*( i-j+1 )) / int2float(2))));
constraint differences[1] < differences[(m*(m-1)) div 2];

solve :: int_search(mark, input_order, indomain, complete)
    minimize mark[m];

output ["golomb ", show(mark), "\n"];

Getting the Row and Column of a Triangular Matrix, Given the Index

New Answer

You can find the row and column using the following formulas:

int row = floor(-0.5 + sqrt(0.25 + 2 * index));
int triangularNumber = row * (row + 1) / 2;
int column = index - triangularNumber;

This works because the first item in each row is a triangular number (0, 1, 3, 6, 10, 15, ...). So the biggest triangular number that is lower than index gives us the row. Then column is simply the difference between index and that triangular number.

Also, note that you don't need the parameter M.

---

Old Answer

This code will give you both the row and column of index.

int triangularNumber = 0;
int row = 0;
while (triangularNumber + row < index) {
    row ++;
    triangularNumber += row;
}
int column = index - triangularNumber;

Upper-triangle matrix looping

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)
jj = ii - K(K+1)/2
return M-1-jj

where M is the size of the matrix
Here's a link algorithm for index numbers of triangular matrix coefficients
which gives more details about the answer algebraically. credit to https://stackoverflow.com/users/4958/shreevatsar

Along with the comment from https://stackoverflow.com/users/2019794/michael-bauer

How to merge an upper and lower triangular, also adding a diagonal of 1s

You can start from the identity matrix and then fill the upper/lower triangular parts.

To access the upper/lower triangular parts of both input matrices and output matrix you can use np.triu_indices and np.tril_indices.

The following code should create the expected out array.

n = len(upper) + 1
out = np.eye(n)
out[np.triu_indices(n, 1)] = upper[np.triu_indices(n-1)]
out[np.tril_indices(n, -1)] = lower[np.tril_indices(n-1)]

---

Related Topics