2

I am working on a computer program (in CUDA C++) that solves issues related to triangular matrices. I can not find a way to correctly index the triangular matrix entries.

Let's assume I have a $N$x$N$ lower triangular matrix.

The maximum number of nonzero entries in such a matrix is $Q = N ( N + 1 ) / 2$.

The program generates a set $W = \{ 0, 1, 2, 3, ..., Q - 1 \}$ to work with.

My responsibility is to map set $W$ into a set of pairs { matrix row, matrix column }

I have those two expressions:
$n = \frac {i}{N}$
$m = i - n * N$

where:

  • $i \in W$
  • $n$ is a row index, $n \in N^0$
  • $m$ is a column index. $m \in N^0$

(I know that $m$ should be equal to $0$ but we are working with natural numbers. In this case, it works as a "round down to the nearest whole number" function.)

let's assume $N = 4$
Those two expressions generate a set of pairs $\{ n, m\}$ from the set $W$:

$$ \begin{matrix} \{0,0\} & \{0,1\} & \{2,0\} & -\\ \{1,0\} & \{1,1\} & \{2,1\} & -\\ \{2,0\} & \{2,1\} & - & -\\ \{3,0\} & \{3,1\} & - & -\\ \end{matrix} $$

(Shown as a matrix to help visualize the issue.)

As you can see this is not a triangular matrix.

Question:

Is there a way to generate a set of pairs $\{ n, m\}$ from the set $W$, which would look like this (for $N = 4$):

$$ \begin{matrix} \{0,0\} & - & - & - \\ \{1,0\} & \{1,1\} & - & - \\ \{2,0\} & \{2,1\} & \{2,2\} & - \\ \{3,0\} & \{3,1\} & \{3,2\} & \{3,3\} \\ \end{matrix} $$

PatrykB
  • 239

2 Answers2

3

Here's one way to do it . . .

For $i \in W$, let $r(i),c(i)$ denote the row and column indices, respectively, for the matrix entry corresponding to the index $i$ (with row and column indices starting at $0$).

Then $r(i)$ and $c(i)$ are given by \begin{align*} r(i)&=\left\lfloor\frac{-1+\sqrt{1+8i}}{2}\right\rfloor\\[4pt] c(i)&=i-\left(\frac{r(i)^2 + r(i)}{2}\right)\\[4pt] \end{align*}

For example, based on the above formulas, $$(r(7),c(7)) = (3,1)$$ which matches your specification.

quasi
  • 58,772
  • The expression $\sqrt{1+8i}$ seems weird, at first sight. :-) – José Carlos Santos Aug 11 '17 at 12:05
  • @quasi is there a way that does not involve a square root? They are expensive to calculate by a computer. – PatrykB Aug 12 '17 at 16:35
  • @cukier9a7b5: One alternative is to do a search for $r(i)$ (e.g., a binary search, successively halving the search interval), based on the fact that $r(i)$ is the largest nonnegative integer such that $$\frac{r(i)^2 + r(i)}{2} \le i$$ – quasi Aug 13 '17 at 00:00
  • @quasi thank you. Can I ask you one last thing? How did You derive the equation for a row index? I understand that you started with something like: $i(i + 1)/2 + j$. Right now I need a similar equation but for a lower triangular matrix without a diagonal. And I would like to exercise it myself. – PatrykB Aug 13 '17 at 12:23
  • @quasi Ok now I understand it... You started with $ x(x + 1) / 2 = i$. In order to get a similar solution for a lower triangular matrix without a diagonal i need to start with $ x(x - 1) / 2 = i$. P.S. I have tested your solution in terms of performance. it's very efficient. One more time thank you. – PatrykB Aug 13 '17 at 12:46
  • @quasi Is there a formula for a matrix of arbitrary dimensions rather than just n=4? – Marko Grdinić Jul 08 '18 at 11:55
  • @ Marko Grdinic: The formula I gave works for any positive integer $n$. – quasi Jul 08 '18 at 13:46
1

For those wishing for an answer without requiring a square-root, here is a variant based upon the approach described in the other answer with Python:

n=5

import math
[(i, (r, i-(rr+r)//2)) for i in range(n  (n+1) // 2) for r in ((-1+math.isqrt(1+8*i))//2,)]

Performing the binary search amounts to:

def triangular_index_to_coordinates(i, n):
    lbound, ubound = 0, n-1
    while lbound <= ubound:
        mid = (lbound + ubound) // 2
        val = (mid * mid + mid) // 2
        if val < i: lbound = mid + 1
        elif val > i: ubound = mid - 1
        else: return (mid, i - val)
    return (ubound, i - (ubound*ubound+ubound)//2)
[(i, triangular_index_to_coordinates(i, n)) for i in range(n * (n+1) // 2)]

The output in both cases is:

[(0, (0, 0)), (1, (1, 0)), (2, (1, 1)), (3, (2, 0)), (4, (2, 1)), (5, (2, 2)), (6, (3, 0)), (7, (3, 1)), (8, (3, 2)), (9, (3, 3)), (10, (4, 0)), (11, (4, 1)), (12, (4, 2)), (13, (4, 3)), (14, (4, 4))]

If a variant without the diagonal included is desired, then the lower bound is adjusted from $0$ to $1$ along with changing $\frac{r(i)^2+r(i)}{2}$ to $\frac{r(i)^2-r(i)}{2}$. Or as code:

def triangular_index_to_coordinates_nodiag(i, n):
    lbound, ubound = 1, n-1
    while lbound <= ubound:
        mid = (lbound + ubound) // 2
        val = (mid * mid - mid) // 2
        if val < i: lbound = mid + 1
        elif val > i: ubound = mid - 1
        else: return (mid, i - val)
    return (ubound, i - (ubound*ubound-ubound)//2)

[(i, triangular_index_to_coordinates_nodiag(i, n)) for i in range(n * (n-1) // 2)]
```
Gregory Morse
  • 117
  • I'm confused. You write "For those wishing for an answer without requiring a square-root", then present a line of code containing the evaluation of a square-root, "for r in ((-1+int(math.sqrt(1+8*i)))//2,)]"... – Gonçalo Nov 11 '23 at 20:55
  • 1
    @Gonçalo, I understand that part to simply be an implementation of the square-root part of the other answer, for comparison. The variant comes afterwards. The wording is a little unclear though indeed. Gregory, a small pointer: from Python 3.8 onwards there is a function math.isqrt which is generally a bit safer than int(math.sqrt(...)). In particular if $N$ is on the order of $2^{53}$ then your implementation might not give the right answer :) – Izaak van Dongen Nov 12 '23 at 00:59
  • Yes it was just an oracle to compare against to show correctness. The code after is the binary search version. As for isqrt, this cannot be used. As that's a floor square root. And here the floor is after subtracting 1 and dividing by 2. I thought about integer Square root as a solution. In fact in Cuda, it might be better if the API supports it otherwise its probably comparable to binary search if having to manually implement it. But to use isqrt you must prove it is an equivalent formula which it does seem to be as the fractional part won't change the outer floor. – Gregory Morse Nov 12 '23 at 03:16
  • Cuda has no builtin integer square root. This algorithm is already a binary search integer square root like algorithm if inspecting closely. See the typical binary search square root algo on Wikipedia which is quite common. Digit by digit variants aren't great either. I'm pretty sure this prototype is quite elegant for Cuda though it is best to bound the while loop iterations so it unrolls based on value of n. – Gregory Morse Nov 12 '23 at 03:24
  • Ah, OK, you may be right. I came here from a review queue so I hadn't fully read the other answer. The reason for my suggestion is that int(math.sqrt(...)) is also a floor square root. So it does exactly the same as what math.isqrt would do. So perhaps your implementation there is not quite right? – Izaak van Dongen Nov 12 '23 at 03:25
  • Yea so you are right, we can prove floor((-1+x)/2)==floor((-1+floor(x))/2) and therefore my implementation is correct but suboptimal to isqrt. Fixed :). Thanks. – Gregory Morse Nov 12 '23 at 03:28