线性索引上三角矩阵

Question

如果我有一个矩阵的上三角形部分,在对角线上方偏移,存储为线性数组,那么如何(i,j)从数组的线性索引中提取矩阵元素的索引？

例如,线性阵列[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9是矩阵的存储

0  a0  a1  a2  a3
0   0  a4  a5  a6
0   0   0  a7  a8
0   0   0   0  a9
0   0   0   0   0

并且我们想要知道数组中的(i,j)索引,该索引对应于线性矩阵中的偏移,而没有递归.

例如,合适的结果k2ij(int k, int n) -> (int, int)将满足

k2ij(k=0, n=5) = (0, 1)
k2ij(k=1, n=5) = (0, 2)
k2ij(k=2, n=5) = (0, 3)
k2ij(k=3, n=5) = (0, 4)
k2ij(k=4, n=5) = (1, 2)
k2ij(k=5, n=5) = (1, 3)
 [etc]

Answer 1

从线性指数到(i,j)指数的方程是

i = n - 2 - floor(sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5)
j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2

从(i,j)索引到线性索引的逆操作是

k = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1

在Python中验证:

from numpy import triu_indices, sqrt
n = 10
for k in range(n*(n-1)/2):
    i = n - 2 - int(sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5)
    j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2
    assert np.triu_indices(n, k=1)[0][k] == i
    assert np.triu_indices(n, k=1)[1][k] == j

for i in range(n):
    for j in range(i+1, n):
        k = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1
        assert triu_indices(n, k=1)[0][k] == i
        assert triu_indices(n, k=1)[1][k] == j

Answer 2

首先，让我们以相反的顺序对a [k]重新编号。我们会得到：

0  a9  a8  a7  a6
0   0  a5  a4  a3
0   0   0  a2  a1
0   0   0   0  a0
0   0   0   0   0

然后k2ij（k，n）将变为k2ij（n-k，n）。

现在的问题是，如何在这个新矩阵中计算k2ij（k，n）。序列0、2、5、9（对角线元素的索引）对应于三角数（减去1后）：a [n-i，n + 1-i] = Ti-1。Ti = i *（i + 1）/ 2，因此，如果我们知道Ti，则很容易求解该方程并得到i （请参阅链接的Wiki文章中的公式，“三角形根和三角形数的测试”部分）。如果k + 1不完全是一个三角数，该公式仍将为您提供有用的结果：将其四舍五入后，将获得i的最大值，对于Ti <= k，i的该值对应于行索引（从底部开始计数），其中a [k]位于其中。要获取该列（从右边开始计数），您应该简单地计算Ti的值并将其减去：j = k + 1-Ti。需要明确的是，这些并不是您问题中的i和j，您需要“翻转”它们。

我没有写确切的公式，但是希望您能想到，现在在执行一些无聊但简单的计算后，找到它就变得微不足道了。

Answer 3

下面是matlab中的实现，可以很容易地转移到另一种语言，比如C++。这里，我们假设矩阵的大小为 m*m，ind 是线性数组中的索引。唯一不同的是，在这里，我们逐列计算矩阵的下三角部分，这与您的情况类似（逐行计算上三角部分）。

function z= ind2lTra (ind, m)
  rvLinear = (m*(m-1))/2-ind;
  k = floor( (sqrt(1+8*rvLinear)-1)/2 );

  j= rvLinear - k*(k+1)/2;

  z=[m-j, m-(k+1)];