压缩矩阵函数以查找配对

Question

对于一组观察数据：

[a1,a2,a3,a4,a5]

它们之间的成对距离

d=[[0,a12,a13,a14,a15]
   [a21,0,a23,a24,a25]
   [a31,a32,0,a34,a35]
   [a41,a42,a43,0,a45]
   [a51,a52,a53,a54,0]]

以压缩矩阵的形式给出（上三角部分，来自 scipy.spatial.distance.pdist 的计算）：

c=[a12,a13,a14,a15,a23,a24,a25,a34,a35,a45]

问题是，既然我已经有了压缩矩阵中的索引，是否有一个函数（最好是用python写的）f 可以快速告诉我是哪两个观察数据被用来计算这些距离的？

f(c,0)=(1,2)
f(c,5)=(2,4)
f(c,9)=(4,5)
...

我尝试了一些解决方案，但没有一个值得一提的 :(

Answer 1

很明显，你要找的这个函数 f 需要一个第二个参数：矩阵的维度——在你的例子中，就是 5。

第一次尝试：

def f(dim,i): 
  d = dim-1 ; s = d
  while i<s: 
    s+=d ; d-=1
  return (dim-d, i-s+d)

Answer 2

压缩矩阵的索引公式是

index = d * (d - 1) / 2 - (d - i) * (d - i - 1) / 2 + j - i - 1

这里的 i 是行索引，j 是列索引，而 d 是原始上三角矩阵的行长度（d X d）。

假设这个索引指的是原始矩阵某一行中最左边的非零元素。对于所有最左边的索引，

j == i + 1

所以

index = d * (d - 1) / 2 - (d - i) * (d - i - 1) / 2 + i + 1 - i - 1
index = d * (d - 1) / 2 - (d - i) * (d - i - 1) / 2

通过一些代数运算，我们可以把这个公式改写为

i ** 2 + (1 - (2 * d)) * i + 2 * index == 0

接着，我们可以用二次公式来找出这个方程的根，而我们只关心正根。

如果这个索引确实对应于最左边的非零单元格，那么我们会得到一个正整数作为解，这个解对应于行号。然后，找到列号就只是简单的算术运算了。

j = index - d * (d - 1) / 2 + (d - i) * (d - i - 1)/ 2 + i + 1

如果这个索引不对应于最左边的非零单元格，那么我们就找不到整数根，但可以取正根的下限作为行号。

def row_col_from_condensed_index(d,index):
    b = 1 - (2 * d) 
    i = (-b - math.sqrt(b ** 2 - 8 * index)) // 2
    j = index + i * (b + i + 2) // 2 + 1
    return (i,j)  

如果你不知道 d 的值，可以通过压缩矩阵的长度来计算出来。

((d - 1) * d) / 2 == len(condensed_matrix)
d = (1 + math.sqrt(1 + 8 * len(condensed_matrix))) // 2 

Answer 3

你可能会觉得 triu_indices 这个函数很有用。比如说，

In []: ti= triu_indices(5, 1)
In []: r, c= ti[0][5], ti[1][5]
In []: r, c
Out[]: (1, 3)

要注意的是，索引是从0开始的。你可以根据自己的需要进行调整，比如：

In []: def f(n, c):
   ..:     n= ceil(sqrt(2* n))
   ..:     ti= triu_indices(n, 1)
   ..:     return ti[0][c]+ 1, ti[1][c]+ 1
   ..:
In []: f(len(c), 5)
Out[]: (2, 4)