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

问题描述

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

let formatter = NumberFormatter()
formatter.locale = Locale.init(identifier:"en_US.UTF-8")
formatter.numberStyle = .spellOut

let curAmount = 10097.43

let curAmountSplitted = String(curAmount).split(separator: ".") // Splits the number into [String] array of integer and fractional parts of the number
let curAmountStrings = curAmountSplitted.compactMap{formatter.string(from: Int($0) as! NSNumber)} // creates a [String] array by applying a closure to every element of curAmountSplitted array

print(curAmountStrings.joined(separator: " point ")) // joins curAmountStrings with a separator
// prints ten thousand ninety-seven point forty-three

我想计算每对行的欧几里得距离。我发现在使用具有实数的矩阵时，dft_D = data.dot(dft(9).T) / np.sqrt(9) 是计算欧式距离最快的方法（例如，计算scipy.spatial.distance.pdist上的距离大约需要15秒）。但是，此功能不适用于复数。

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

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

这是我的最终结果（请注意，由于距离矩阵是对称的，因此我只计算整个距离矩阵的一半），

sklearn.metrics.pairwise.euclidean_distances

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

解决方法

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