线性索引上三角矩阵答案

【问题标题】：Linear index upper triangular matrix线性索引上三角矩阵
【发布时间】：2015-01-21 01:33:22
【问题描述】：

如果我有一个矩阵的上三角部分，在对角线上偏移，存储为一个线性数组，如何从数组的线性索引中提取矩阵元素的(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]

【问题讨论】：

为最后一列的元素写一个公式。为了更容易编写一个从行号（列号是固定的）计算线性索引的公式，然后将其反转。从那里开始一个通用公式。
请注意，这里介绍的求解方法也可用于列出一次取 2 件的 N 个事物的组合（无需重复），无需任何迭代/递归。
回答了没有偏移量的相同问题here

标签： c++ arrays numpy linear-algebra triangular

【解决方案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

【讨论】：

完美！这帮助我减少了线性程序中的变量数量！
i = ...多了一个括号
你能解释一下K的公式/它是如何驱动的吗？
你需要更多的括号来确保避免在 2 下潜入时将临时值四舍五入：k = ((n*(n-1))/2) - ((ni)*( (ni)-1))/2 + j - i - 1
请注意，这些方程适用于零索引变量...

【解决方案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（对角线元素的索引）对应于triangular numbers（减去 1 后）：a[n - i, n + 1 - i] = Ti - 1. Ti = i * (i + 1)/2，所以如果我们知道 Ti，就很容易求解这个方程并得到 i（参见链接的 wiki 文章中的公式，“三角根和三角数测试”部分）。如果 k + 1 不完全是一个三角数，这个公式仍然会给你有用的结果：四舍五入后，你会得到 i 的最大值，对于其中 Ti

我没有写出确切的公式，但我希望你能明白这一点，在进行一些无聊但简单的计算之后，现在找到它是微不足道的。

【讨论】：

谢谢，这真的帮助我理解了解决方案

【解决方案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)];

【讨论】：

【解决方案4】：

对于记录，这是相同的函数，但使用从一开始的索引，并且在 Julia 中：

function iuppert(k::Integer,n::Integer)
  i = n - 1 - floor(Int,sqrt(-8*k + 4*n*(n-1) + 1)/2 - 0.5)
  j = k + i + ( (n-i+1)*(n-i) - n*(n-1) )÷2
  return i, j
end

【讨论】：

【解决方案5】：

在 python 2 中：

def k2ij(k, n):
    rows = 0
    for t, cols in enumerate(xrange(n - 1, -1, -1)):
        rows += cols
        if k in xrange(rows):
            return (t, n - (rows - k))
    return None

【讨论】：

反之的ij2k 函数呢？ @smac89