c++arrayslinear-algebratriangular

Giving an element from lower/upper triangular matrix


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,


Solution

  • 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 😇.