计算具有复数的非常大矩阵的欧几里得距离的最快方法是什么？

Question

我有一个非常大的输入数据集，包含 50,000 个 9 维样本（即 50000x9 矩阵）。此数据已使用 DFT 进行转换：

dft_D = data.dot(dft(9).T) / np.sqrt(9)

我想计算每对行的欧几里得距离。当使用实数矩阵时，我发现scipy.spatial.distance.pdist 在计算欧几里得距离方面是最快的（例如，计算data 上的距离需要~`15 秒）。但是，此函数不适用于复数。

我尝试了this SO post 中提出的解决方案，但这给我带来了严重的内存问题（即“无法为形状为 (50000, 50000, 9) 且数据类型为 complex128 的数组分配 191.GiB”）。我也尝试过使用 this Medium article 中定义的 EDM，但这也给了我类似的内存问题。

最初，我能够通过使用定义np.sqrt(np.sum(np.square(np.abs(data[i,:] - data[j,:])))) 遍历行和列来计算这些欧几里得距离。这非常慢。然后我将docs 中描述的定义用于sklearn.metrics.pairwise.euclidean_distances（它也不适用于复数），它稍微快一些，但仍然很慢（运行超过2 小时）。

这是我的最终结果（注意，由于距离矩阵是对称的，所以我只计算了完整距离矩阵的一半），

import numpy as np
def calculate_euclidean_distance(arr, num_rows):
    dist_matrix = np.empty(int((num_rows*(num_rows - 1))/2))
    idx = 0
    dot_dict = {}
    # get the 0th row out of the way
    dot_dict[0] = arr[0,:].dot(arr[0,:])
    
    for i in range(1,num_rows):
        # Save the value of dot(X,X) in dict to not recompute it every time when needed
        if i not in dot_dict:
            dot_dict[i] = arr[i,:].dot(arr[i,:])
        i_dot = dot_dict[i]
        for j in range(0,i):
            j_dot = dot_dict[j]
            dist_matrix[idx] = np.sqrt(i_dot - 2*arr[i,:].dot(arr[j,:]) + j_dot)
            idx+=1
    return dist_matrix

当涉及复数时，有没有更快的方法来获得这些距离？

Answer 1

您可以使用 numpy.roll() 以循环方式移动输入数组的行。它重复了很多计算，但尽管如此，它还是要快得多。下面的代码填充了距离矩阵的下半部分

dist_matrix = np.empty(shape = [inp_arr.shape[0], inp_arr.shape[0]])
for i in range(inp_arr.shape[0]):
    shifted_arr = np.roll(inp_arr, i, axis = 0)
    curr_dist = np.sqrt(np.sum(np.square(np.abs(inp_arr - shifted_arr)), axis = 1))
    for j in range(i, inp_arr.shape[0]):
        dist_matrix[j, j - i] = curr_dist[j]

Answer 2

我不明白你对dft_D 的定义。但是，如果您尝试计算原始数据的 DFT 行之间的距离，这将与原始数据的行之间的距离相同。

根据Parseval's theorem，向量的大小及其变换是相同的。并且通过线性，两个向量的差的变换等于它们变换的差。由于欧几里得距离是差值大小的平方根，因此使用哪个域来计算度量并不重要。我们可以用一个小样本来演示：

import numpy as np
import scipy.spatial

x = np.random.random((500,9)) #Use a smaller data set for the demo
Sx = np.fft.fft(x)/np.sqrt(x.shape[1]) #numpy fft doesn't normalize by default
xd = scipy.spatial.distance.pdist(x,metric='euclidean')
Sxd = np.array([np.sqrt(np.sum(np.square(np.abs(Sx[i,:] - Sx[j,:])))) for i in range(Sx.shape[0]) for j in range(Sx.shape[0])]).reshape((Sx.shape[0],Sx.shape[0])) #calculate the full square of pairwise distances
Sxd = scipy.spatial.distance.squareform(Sxd) #use scipy helper function to get back the same format as pdist
np.all(np.isclose(xd,Sxd)) # Should print True

因此，只需在原始数据上使用pdist。