凝聚矩阵函数找到对

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

凝聚基质的指数的公式是

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 - 2d)*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 = math.floor((-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 2

您可能会发现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)