Mapping Elements in 2D Upper Triangle and Lower Triangle to Linear Structure

Mapping elements in 3D lower triangle to linear structure

The 3D version of the equation is

index := (z * (z+1) * (z+2)) / 6 + (y * (y+1))/2 + x

efficient way to represent a lower/upper triangular matrix

Really, you're best off just using a regular two dimensional matrix. RAM is pretty cheap. If you really don't want to do that, then you can build a one-dimensional array with the right number of elements and then figure out how to access each element. For example, if the array is structured like this:

    j
    1234
i 1 A
  2 BC
  3 DEF
  4 GHIJ

and you have it stored as a one dimensional array, left to right, you'd access element C (2, 2) with array[3]. You can work out a function to go from [i][j] to [n] but I won't spoil your fun. But you don't have to do this unless the triangular array in question is really huge or you're very concerned about space.

Mapping a triangular matrix to vector

@reve_etrange I believe N should be the length of d, not the number of rows of your points matrix.

K=5; % Number of points
d = pdist(rand(K,3));
N = length(d);
m = ceil(sqrt(2*N));
I = 1:N;
CI = m - round( sqrt( 2*(1 + N - I) ) );
RI = mod(I + CI.*(CI+1)/2 - 1, m) + 1;

This results in

I =    1     2     3     4     5     6     7     8     9    10
RI =   2     3     4     5     3     4     5     4     5     5
CI =   1     1     1     1     2     2     2     3     3     4

Mapping a triangular matrix to vector

@reve_etrange I believe N should be the length of d, not the number of rows of your points matrix.

K=5; % Number of points
d = pdist(rand(K,3));
N = length(d);
m = ceil(sqrt(2*N));
I = 1:N;
CI = m - round( sqrt( 2*(1 + N - I) ) );
RI = mod(I + CI.*(CI+1)/2 - 1, m) + 1;

This results in

I =    1     2     3     4     5     6     7     8     9    10
RI =   2     3     4     5     3     4     5     4     5     5
CI =   1     1     1     1     2     2     2     3     3     4

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

CUDA loop over lower triangular matrix

The problem is that we are indexing a 1D array so in order to map it we need to multiply the row index with the number of columns, therefore following the example:

__global__ void Kernel(int N) {
        int row = blockIdx.x * blockDim.x + threadIdx.x;
        int col = blockIdx.y * blockDim.y + threadIdx.y;
        if((row < N) && (col<=row) )
           printf("%d\n", row*N + col);
 }

Neighbor index computation for diagonally flattened matrix

I actually already had the elements to solve it somewhere else in my code. As BLUEPIXY's solution hinted, I am using scatter/gather operations, which I had already implemented for layout transformation.

This solution basically rebuilds the original (x,y) index of the given element in the matrix, applies the index offset and translates the result back to the transformed layout. It splits the square in 2 triangles and adjust the computation depending on which triangle it belongs to.

It is an almost entirely algebraic transformation: it uses no loop and no table lookup, has a small memory footprint and little branching. The code can probably be optimized further.

Here is the draft of the code:

#include <stdio.h>
#include <math.h>

// size of the matrix
#define SIZE 4

// triangle number of X
#define TRIG(X) (((X) * ((X) + 1)) >> 1)
// triangle root of X
#define TRIROOT(X) ((int)(sqrt(8*(X)+1)-1)>>1);

// return the index of a neighbor given an offset
int getNGonalNeighbor(const size_t index,
                      const int x_offset,
                      const int y_offset){
    // compute largest upper triangle index
    const size_t upper_triangle = TRIG(SIZE);

    // position of the actual element of index
    unsigned int x = 0,y = 0;

    // adjust the index depending of upper/lower triangle.
    const size_t adjusted_index = index < upper_triangle ?
                index :
                SIZE * SIZE - index - 1;

    // compute triangular root
    const size_t triroot = TRIROOT(adjusted_index);
    const size_t trig = TRIG(triroot);
    const size_t offset = adjusted_index - trig;

    // upper triangle
    if(index < upper_triangle){
        x = offset;
        y = triroot-offset;
    }
    // lower triangle
    else {
        x = SIZE - offset - 1;
        y = SIZE - (trig + triroot + 1 - adjusted_index);
    }

    // adjust the offset
    x += x_offset;
    y += y_offset;

    // manhattan distance
    const size_t man_dist = x+y;

    // calculate index using triangular number
    return TRIG(man_dist) +
            (man_dist >= SIZE ? x - (man_dist - SIZE + 1) : x) -
            (man_dist > SIZE ? 2* TRIG(man_dist - SIZE) : 0);
}

int main(){
    printf("%d\n", getNGonalNeighbor(15,-1,-1)); // should return 11
    printf("%d\n", getNGonalNeighbor(15, 0,-1)); // should return 14
    printf("%d\n", getNGonalNeighbor(15,-1, 0)); // should return 13
    printf("%d\n", getNGonalNeighbor(11,-2,-1)); // should return 1
}

And the output is indeed:

11
14
13
1

If you think this solution looks over complicated and inefficient, I remind you that the target here is GPU, where computation costs virtually nothing compared to memory accesses, and all index computations are computed at the same time using massively parallel architectures.

Multiplying upper triangular matrix stored in 1D array alogrithm

Let's say the full matrix is of size N x N.

In the first row, row = 0, there are N elements.

In the second row, row = 1, there are N - 1 elements.

...

In the k-th row, row = k, there are N - k elements.

Here's an ASCII diagram:

|<-               N                    ->|

+--+--+--+       ---            +--+--+--+
|  |  |  |                      |  |  |  |  row = 0, N elements
+--+--+--+       ---            +--+--+--+
   |  |  |                      |  |  |  |  row = 1, N-1 elements
   +--+--+       ---            +--+--+--+
      |  |                      |  |  |  |  row = 2, N-2 elements
      +--+       ---            +--+--+--+

                 ---            +--+--+--+
                                |  |  |  |  row = k, N-k elements
                 ---            +--+--+--+

                                +--+--+--+
                                |  |  |  |  row = N-3, 3 elements
                                +--+--+--+
                                   |  |  |  row = N-2, 2 elements
                                   +--+--+
                                      |  |  row = N-1, 1 element
                                      +--+

The number of elements needed in the 1D array to store such a matrix is:

1 + 2 + ... + N = N*(N+1)/2

To access an element of the matrix, we need to map a row index and column index to an index in the 1D array.

To get the j-th column of the first row (row index = 0, column index = j),

index = j

To get the j-th column of the second row (row index = 1, column index = j),

index = N + (j-1)

To get the j-th column of the third row (row index = 2, column index = j),

index = N + (N-1) + (j-2)

...

To get the j-th column of the i-th row (row index = i, column index = j),

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

With the above information, you can write functions that:

1. Allocate the right amount of memory for the upper triangular matrix
2. Set the value of the matrix given a row index (i) and column index (j).
3. Get the value of the matrix given a row index (i) and column index (j).

That should be adequate to be able to multiply an upper triangular matrix with another upper triangular matrix as well as a regular matrix.

---

Related Topics