I have the upper triangular portion of a matrix, with the main diagonal stored as a linear array, how can the (i,j) indices of a matrix element be extracted from the linear index of the array?
For example the linear array :[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10]
is storage for the matrix
a0 a1 a2 a3
0 a4 a5 a6
0 0 a7 a8
0 0 0 a10
I've found solutions for this problem but without the main diagonal which is:
index = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1
And solution for the same problem but for a lower triangular matrix with the diagonal:
index = ((i + 1) * i / 2 + i).
Regards,
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 👌
Worth a blog post explaining how to arrive at this conclusion, but this’ll do for now 😇.