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

What is the fastest way to compute the Euclidean distances of a very large matrix with complex numbers?

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

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

我想计算每对行的欧式距离。我发现 scipy.spatial.distance.pdist 在使用实数矩阵时计算欧氏距离最快(例如,计算 data 上的距离需要约 15 秒)。但是,此函数不适用于复数。

我已经尝试了 this SO post, but this gave me serious memory issues (i.e. "Unable to allocate 191. GiB for an array with shape (50000, 50000, 9) and data type complex128"). I have also tried using the EDM defined in this Medium article 中提供的解决方案,但这也给了我类似的内存问题。

最初,我能够通过使用定义 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

当涉及复数时,是否有更快的方法来获得这些距离?

您可以使用 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]

我不明白你对 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即可。