Map upper triangular matrix on vector skipping the diagonal

Question

I have a problem that could be boiled down to finding a way of mapping a triangular matrix to a vector skipping the diagonal.

Basically I need to translate this C++ code using the Gecode libraries

// implied constraints
for (int k=0, i=0; i<n-1; i++)
  for (int j=i+1; j<n; j++, k++)
    rel(*this, d[k], IRT_GQ, (j-i)*(j-i+1)/2);

Into this MiniZinc (functional) code

constraint
   forall ( i in 1..m-1 , j in i+1..m )
      ( (differences[?]) >=  (floor(int2float(( j-i )*( j-i+1 )) / int2float(2)) ));

And I need to figure out the index in differences[?].

MiniZinc is a functional/mathematical language with no proper for loops. So I have to map those indexes i and j that are touching all and only the cells of an upper triangular matrix, skipping its diagonal, to a k that numbers those cells from 0 to whatever.

If this was a regular triangular matrix (it's not), a solution like this would do

index = x + (y+1)*y/2

The matrix I'm handling is a square n*n matrix with indexes going from 0 to n-1, but it would be nice to provide a more general solution for an n*m matrix.

Here's the full Minizinc code

% modified version of the file found at 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[int] of var 0..n: differences = [mark[j] - mark[i] | i in 1..m, j in i+1..m];

constraint mark[1] = 0;

constraint forall ( i in 1..m-1 ) ( mark[i] < mark[i+1] );

% this version of the constraint works
constraint forall ( i in 1..m-1 , j in i+1..m )
    ( (mark[j] - mark[i]) >= (floor(int2float(( j-i )*( j-i+1 )) / int2float(2))) );
    
%this version does not
%constraint forall ( i in 1..m-1, j in i+1..m )
%    ( (differences[(i-1) + ((j-2)*(j-1)) div 2]) >= (floor(int2float(( j-i )*( j-i+1 )) / int2float(2))) );

constraint alldifferent(differences);

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

Thanks.

Sansen

Answer 1

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

Deirdredeism

I'm still quite interested in the answer. I'm going to use this very formula again for another project. Could you please have a look at the code? Thanks. – Sansen

Sorry for the long reply. Take a look at my edit. Key point: use T(i,j) = i + (j-1)*(j-2))/2 to bijectively map {(i,j) | 1 <= i < j <= N} to {1,2,...,(N*(N-1))/2}. The original code from minizinc uses tuples i > j (the lower triangle), but my code uses tuples j > i (the upper triangle). Either way would work but you asked for upper triangle in question so I gave the code for that. If you want the code for lower triangular too, feel free to ask. – Deirdredeism

Remove the list comprehension. Instead of array[int] of var 0..n: differences = [mark[j] - mark[i] | i in 1..m, j in i+1..m]; you just need array[1..(m*(m-1)) div 2] of var 0..n: differences;. – Deirdredeism

Thanks. Add this additional constraint 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)))); and now it's faster than the original. Thanks! I updated the question with the correct code. If you want, tomorrow I can remove the solution code and you can post it in your answer. If you'd like to add the lower matrix formula too, that would be great. – Sansen

Thanks for the advice. I modified this post but I can't delete the other. I get an error message: 'You cannot delete this accepted answer' when I try. Perhaps if you accept this answer SO will then allow me to delete the other? – Deirdredeism

Okay, I deleted the other post and hopefully cleaned this one up enough. Thanks again :) – Deirdredeism