Linear index for a diagonal run of an upper triangular matrix

Question

Given a NxN matrix, I would like to linearly index into its upper right triangle, following a diagonal by diagonal pattern, starting after the main diagonal.

For example, given a 4x4 matrix

X 0 3 5
X X 1 4
X X X 2
X X X X

I'm looking for a non recursive (closed form) function mapping linear indices from 0 to 5 to (x,y) achieving

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

Related for row by row runs:

• Linear index upper triangular matrix
• algorithm for index numbers of triangular matrix coefficients

Answer 1

This is fairly easy if we use numpy with python. I am writing the code if any Python users come by.

import numpy as np
import scipy.sparse

def linear_to_upper_triangle_index(N, idx):
    # Create a sparse upper triangular matrix to find the coordinates
    #The k=1 argument specifies that we want to exclude the main diagonal.
    upper_tri_indices = np.triu_indices(N, k=1)
    # Get the specific index from the upper triangular indices
    row_index = upper_tri_indices[0][idx]
    col_index = upper_tri_indices[1][idx]
    row_col_idx = (row_index,col_index)
    return (row_col_idx)

# Example usage:
N = 4
for idx in range(4):
    print('idx :',idx)
    print(f"f({idx}) = {linear_to_upper_triangle_index(N, idx)}")
'''
idx : 0
f(0) = (0, 1)
idx : 1
f(1) = (0, 2)
idx : 2
f(2) = (0, 3)
idx : 3
f(3) = (1, 2)

'''

Answer 2

So you want the inverse of the following function

• Zero-based indexing form of element [i,j] for a n×n upper triangular matrix including the diagonal

  index = i*n-i*(i+1)/2+j
  
  i=0..4, j=0..4, index=
  |  0 |  1 |  2 |  3 |  4 |
  |  X |  5 |  6 |  7 |  8 |
  |  X |  X |  9 | 10 | 11 |
  |  X |  X |  X | 12 | 13 |
  |  X |  X |  X |  X | 14 |
  

The easiest algorithm I can think of is to loop for all rows i and see if there is a match for the column j such that:

• i <= j
• j>=0
• j<n

Here is a sample code given index and n

for(i=0; i<n; i++)
{
    j = index - i*n + i*(i+1)/2
    if( j>=0 && j<n && j>= i)
    {
       break;
    }
}

And example with n=7 and [i,j]=[1,5] produces index=11. Now the coordinates of this index are

i	j	i<=j && j>=0 && j<7
0	11
1	5	valid
2	0
3	-4
4	-7
5	-9
6	-10

---

If you want strictly the upper triangular elements, excluding the diagonal then

• Zero-based indexing form of element [i,j] for a n×n upper triangular matrix excluding the diagonal

  index = i*n-i*(i+3)/2+j-1
  
  i=0..3, j=0..4, index=
  |  X |  0 |  1 |  2 |  3 |
  |  X |  X |  4 |  5 |  6 |
  |  X |  X |  X |  7 |  8 |
  |  X |  X |  X |  X |  9 |
  |  X |  X |  X |  X |  X |
  

The algorithm now is to loop for all rows i and see if there is a match for the column j such that:

• i < j
• j>0
• j<n

Here is a sample code given index and n

for(i=0; i<n; i++)
{
    j = index - i*n + i*(i+3)/2 + 1
    if( j>0 && j<n && j>i)
    {
       break;
    }
}

And example with n=7 and [i,j]=[1,5] produces index=9. Now the coordinates of this index are

i	j	i<j && j>0 && j<7
0	10
1	5	valid
2	1
3	-2
4	-4
5	-5

Answer 3

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!

Answer 4

Maybe someone can come up with a math formula that doesn't require a loop, but until then I've come up with a O(N) solution:

#include <utility>

constexpr std::pair<int, int> f(int n, int idx)
{
    int group_size = n - 1;
    int rest = idx + 1;

    while (rest > group_size)
    {
        rest = rest - group_size;
        --group_size;
    }
    return {(rest - 1) % group_size,
            n - group_size +  (rest - 1) % group_size};
}

/* 3x3
X 0 2
X X 1
X X X
*/
static_assert(f(3, 0) == std::pair{0, 1});
static_assert(f(3, 1) == std::pair{1, 2});
static_assert(f(3, 2) == std::pair{0, 2});

// 4x4
static_assert(f(4, 0) == std::pair{0, 1});
static_assert(f(4, 1) == std::pair{1, 2});
static_assert(f(4, 2) == std::pair{2, 3});
static_assert(f(4, 3) == std::pair{0, 2});
static_assert(f(4, 4) == std::pair{1, 3});
static_assert(f(4, 5) == std::pair{0, 3});

/* 5x5
X 0 4 7 9
X X 1 5 8
X X X 2 6
X X X X 3
X X X X X
*/
static_assert(f(5, 0) == std::pair{0, 1});
static_assert(f(5, 1) == std::pair{1, 2});
static_assert(f(5, 2) == std::pair{2, 3});
static_assert(f(5, 3) == std::pair{3, 4});
static_assert(f(5, 4) == std::pair{0, 2});
static_assert(f(5, 5) == std::pair{1, 3});
static_assert(f(5, 6) == std::pair{2, 4});
static_assert(f(5, 7) == std::pair{0, 3});
static_assert(f(5, 8) == std::pair{1, 4});
static_assert(f(5, 9) == std::pair{0, 4});

Answer 5

I created a custom method for the array and value you gave.

int a[4][4] ={{-1, 0, 3, 5}, {-1, -1, 1, 4}, {-1, -1, -1, 2}, {-1, -1, -1, -1}};

The code is exactly like this. You give the array And whatever you give to the second value in the Func method, the indexes of the value in the upper diagonal will reach you.

#include <iostream>

using namespace std;

int b[2] ={-1,-1};

int Func(int a[4][4],int n)
{
    
    for(int i =0;i<4;i++)
    {
        for(int j=0;j<4;j++)
        {
            if(a[i][j]==n)
            {
                if(i<j)
                {
                    b[0]=i;
                    b[1]=j;
                    return 0;
                }
            }
        }
    }
}
int main()
{
    int a[4][4] ={{-1, 0, 3, 5}, {-1, -1, 1, 4}, {-1, -1, -1, 2}, {-1, -1, -1, -1}};
    Func(a,5);
    for(int i=0;i<2;i++)
    {
        cout<<b[i]<<" ";
    }
    return 0;
}

thank you USEFUL for feedback if it worked for you