Numpy：将方阵的索引转换为其上三角索引的有效方法[关闭]答案

【问题标题】：Numpy: Efficient way to convert indices of a square matrix to its upper triangular indices [closed]Numpy：将方阵的索引转换为其上三角索引的有效方法[关闭]
【发布时间】：2018-11-09 21:52:43
【问题描述】：

问题：给定一个索引元组，返回它在上三角索引中的顺序。这是一个例子：

假设我们有一个形状为 (3, 3) 的方阵 A。

A有6个上三角索引，分别是(0, 0), (0, 1), (0, 2), (1, 1), (1, 2), (2, 2)。

现在我知道索引 (1, 2) 处的一个元素，该索引属于 A 的上三角部分。我想返回 4（这意味着它是所有上三角索引中的第 5 个元素。）

关于如何做到这一点的任何想法？

最好，志豪

【问题讨论】：

更一般地说，如果我有一个索引列表，我们是否有一个函数可以将它们一起转换为 triu 索引 [返回转换结果的列表/数组]？

标签： python numpy matrix

【解决方案1】：

可以写出显式公式：

def utr_idx(N, i, j):
    return (2*N+1-i)*i//2 + j-i

演示：

>>> N = 127
>>> X = np.transpose(np.triu_indices(N))
>>> utr_idx(N, *X[2123])
2123

【讨论】：

我才意识到哈哈
反过来呢？从向量返回到矩阵 (i,j) 索引？
Willem Van Onsem 在他的回答@makis 中写下了逆。

【解决方案2】：

对于一个n×n矩阵，上三角的第(i, j)项是i×(2×n-i +1)/2+j-i-矩阵的第一个元素。

我们也可以反过来做数学运算，得到第 k 个元素的 (i, j) 元素：

i = ⌊(-√((2n+1)²-8k)+2n+1)/2⌋ 和 j = k+i-i ×(2×n-i+1)/2

例如：

from math import floor, sqrt

def coor_to_idx(n, i, j):
    return i*(2*n-i+1)//2+j-i

def idx_to_coor(n, k):
    i = floor((-sqrt((2*n+1)*(2*n+1)-8*k)+2*n+1)/2)
    j = k + i - i*(2*n-i+1)//2
    return i, j

例如：

>>> [idx_to_coor(4, i) for i in range(10)]
[(0, 0), (0, 1), (0, 2), (0, 3), (1, 1), (1, 2), (1, 3), (2, 2), (2, 3), (3, 3)]
>>> [coor_to_idx(4, i, j) for i in range(4) for j in range(i, 4)]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

鉴于数字并不大（如果这些数字很大，则计算不再在恒定时间内完成），因此我们可以计算 k-th 坐标 O(1) ，例如：

>>> idx_to_coor(1234567, 123456789)
(100, 5139)

相当于通过枚举获得：

>>> next(islice(((i, j) for i in range(1234567) for j in range(i, 1234567)), 123456789, None))
(100, 5139)

此处将索引转换为坐标也可能会由于浮点不精确而产生一些舍入误差。

【讨论】：

【解决方案3】：

IIUC，你可以使用itertools组合with replacement获取索引

>>> ind = tuple(itertools.combinations_with_replacement(range(3),2))
((0, 0), (0, 1), (0, 2), (1, 1), (1, 2), (2, 2))

要检索索引，只需使用index 方法

>>> ind.index((1,2))
4

【讨论】：

【解决方案4】：

您可以使用np.triu_indices 和dictionary：

import numpy as np

iu1 = np.triu_indices(3)
table = {(i, j): c for c, (i, j) in enumerate(zip(*iu1))}
print(table[(1, 2)])

输出

【讨论】：

【解决方案5】：

与@DanielMesejo 类似，您可以将np.triu_indices 与argwhere 或nonzero 一起使用：

my_index = (1,2)

>>> np.nonzero((np.stack(np.triu_indices(3), axis=1) == my_index).all(1))
(array([4]),)
>>> np.argwhere((np.stack(np.triu_indices(3), axis=1) == my_index).all(1))
array([[4]])

解释：

np.stack(np.triu_indices(3), axis=1) 按顺序为您提供上三角形的索引：

array([[0, 0],
       [0, 1],
       [0, 2],
       [1, 1],
       [1, 2],
       [2, 2]])

因此，您所要做的就是找到与 [1,2] 匹配的位置（您可以使用 == 运算符和 all 来完成）

【讨论】：

【解决方案6】：

构建上层索引的成本很高。我们可以像这样直接获取对应的索引-

def triu_index(N, x, y):
    # Get index corresponding to (x,y) in upper triangular list
    idx = np.r_[0,np.arange(N,1,-1).cumsum()]
    return idx[x]+y-x

示例运行 -

In [271]: triu_index(N=3, x=1, y=2)
Out[271]: 4

【讨论】：