有效地计算只有一个对称外积矩阵的独特元素

Question

我有一个大的集合，我通过取自外积起重Ñd维向量（驻留在一个矩阵）（即，每个矢量I倍本身）的。对于每一个矢量，这产生与对称矩阵（d + 1）选择两个独特的条目。针对整个数据，这是一个N×d X d张量。我想仅计算唯一的（d + 1）从下部对角张量的每个切片的选择项2，并将它们存储在一个载体中。我想用尽可能少的内存占用和尽可能快地在Python这样做 - 包括使用C绑定。

如果你做到这一点使用标准numpy的方法，它分配每个矩阵的全部。这是关于双什么是实际需要的存储复杂。

对于规模感这里，考虑如下情况，其中N = 20K和d = 20K。然后N * d ^ 2 *〜每元素=（2 * 10 ^ 4）^ 3 * 8个字节= 64个TB的8个字节。

如果只计算了编码的唯一条目的向量，我们有（20001选择2）* 20K×8 = 200010000 * 20000 * 8个字节= 32兆兆字节。

是否有一个快速的方法来做到这一点，而不诉诸慢的方法（如编码我自己外积在Python）？

编辑：我注意到，类似的问题被问Create array of outer products in numpy

我已经知道如何计算这个使用einsum（如上面的问题）。然而，没有答案是达到了约这样做没有额外的（D选择2）计算和分配

编辑2：此线程How to exploit symmetry in outer product in Numpy (or other Python solutions)?问一个相关的问题，但并没有解决内存的复杂性。顶端答案将仍然分配用于每个外产物d X d阵列。

此线程Numpy Performance - Outer Product of a vector with its transpose还讨论了自外积的计算考虑，但没有达到内存有效的解决方案。

编辑3：如果一个人想要分配整个数组，然后提取元素，np.tril_indices或scipy.spatial.distance.squareform会做的伎俩。

Answer 1

不知道究竟如何你想你的输出，但总是有使用Numba的选项：

import numpy as np
import numba as nb

# Computes unique pairwise products
@nb.njit(parallel=True)
def self_outer_unique(a):
    n, d = a.shape
    out = np.empty((n, (d * d + d) // 2), dtype=a.dtype)
    for i in nb.prange(n):
        for j1 in range(d):
            for j2 in range(j1, d):
                idx = j1 * (2 * d - j1 + 1) // 2 + j2 - j1
                out[i, idx] = a[i, j1] * a[i, j2]
    return out

这将使您对每行（即，满输出的扁平的上部三角形）中的所有独特的产品阵列。

import numpy as np

a = np.arange(12).reshape(4, 3)
print(a)
# [[ 0  1  2]
#  [ 3  4  5]
#  [ 6  7  8]
#  [ 9 10 11]]
print(self_outer_unique(a))
# [[  0   0   0   1   2   4]
#  [  9  12  15  16  20  25]
#  [ 36  42  48  49  56  64]
#  [ 81  90  99 100 110 121]]

性能的角度来看，它是比计算与NumPy的全外产物更快，虽然重新创建全阵列从该需要更长的时间。

import numpy as np

def np_self_outer(a):
    return a[:, :, np.newaxis] * a[:, np.newaxis, :]

def my_self_outer(a):
    b = self_outer_unique(a)
    n, d = a.shape
    b_full = np.zeros((n, d, d), dtype=a.dtype)
    idx0 = np.arange(n)[:, np.newaxis]
    idx1, idx2 = np.triu_indices(d)
    b_full[idx0, idx1, idx2] = b
    b_full += np.triu(b_full, 1).transpose(0, 2, 1)
    return b_full

n, d = 1000, 100
a = np.arange(n * d).reshape(n, d)
print(np.all(np_self_outer(a) == my_self_outer(a)))
# True

%timeit np_self_outer(a)
# 24.6 ms ± 248 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit self_outer_unique(a)
# 6.32 ms ± 69.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit my_self_outer(a)
# 124 ms ± 770 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)