产生稀疏成对距离矩阵python避免内存错误答案

【问题标题】：Produce sparse pairwise distance matrix python avoiding memory error产生稀疏成对距离矩阵python避免内存错误
【发布时间】：2020-10-13 12:19:47
【问题描述】：

我有一些地理空间数据，我正在尝试生成点之间的成对半正弦距离。不幸的是，大约有 50k 个数据点，这将产生一个 50000x50000 的距离矩阵。由于这个距离矩阵是对称的，我认为这个矩阵可能是稀疏的。我希望能够直接从距离计算中生成一个稀疏矩阵，这样我就不会出现内存错误。到目前为止，这是我的代码：

def convert_to_arrays(df1, df2):
    d1 = np.array(df1[['x','y']].values.tolist())
    d2 = np.array(df2[['x','y']].values.tolist())
    return d1,d2

def haversine(data1, data2):
    data1 = np.deg2rad(data1)                     
    data2 = np.deg2rad(data2)                     

    lat1 = data1[:,0]                     
    lng1 = data1[:,1]         

    lat2 = data2[:,0]                     
    lng2 = data2[:,1]         

    diff_lat = lat1[:,None] - lat2
    diff_lng = lng1[:,None] - lng2
    d = np.sin(diff_lat/2)**2 + np.cos(lat1[:,None])*np.cos(lat2) * np.sin(diff_lng/2)**2
    return 2 * 6371 * np.arcsin(np.sqrt(d))

dist = haversine(*convert_to_arrays(df, df))

其中df 是带有“x”和“y”列的 pandas 数据框。当我以大约 50k 点运行此代码时，出现内存错误：

MemoryError: Unable to allocate 19.3 GiB for an array with shape (50893, 50893) and data type float64

有没有办法确保使用（下三角或上三角）稀疏矩阵作为预期输出来执行计算？

【问题讨论】：

数组大小太大。
您可以使用 float32 精度创建这种大小的稀疏矩阵。例如，A = scipy.sparse.csr_matrix((50893,50893), dtype=np.float32).toarray()。我希望能够从一开始就计算一个稀疏的下三角距离矩阵，而不必先计算一个密集的矩阵。
当稀疏度为 50%（三角距离矩阵将是）时，我不知道任何比标准浮点数组更有效的稀疏数据结构。
@CJR, scipy.distance 为此使用压缩距离矩阵表示。基本上，您只是存储一维距离数组，其中条目对应于[dist(x[0], x[1]), dist(x[0], x[2]), ..., dist(x[0], x[n]), dist(x[1], x[2]), dist(x[1], x[3]), ..., dist(x[n-1], x[n])]。 @CopyOfA，我不相信 scipy 有任何稀疏矩阵格式可以满足您的需求。

标签： python scipy sparse-matrix distance-matrix

【解决方案1】：

我发现实现这一点的方法相当蛮力；如果您有更优雅的解决方案，请发布。本质上，我从第一个数据帧中取出点，并计算与第二个数据帧中所有点的成对距离，产生一个 Nx50893 数组。然后我将此数组转换为稀疏的下三角矩阵。最后，我计算下一个块，转换为下三角形并将这些数组堆叠起来。

def haversine(data1, data2):
    data1 = np.deg2rad(data1)                     
    data2 = np.deg2rad(data2)                     

    lat1 = data1[:,0]                     
    lng1 = data1[:,1]         

    lat2 = data2[:,0]                     
    lng2 = data2[:,1]         

    #chunk it out 
    N = 2000 #chunk size
    start = 0
    end = N
    #First iteration
    diff_lat = lat1[:N,None] - lat2
    diff_lng = lng1[start:end,None] - lng2()
    d = np.sin(diff_lat/2)**2 + np.cos(lat1[start:end,None])*np.cos(lat2) * np.sin(diff_lng/2)**2
    D = scipy.sparse.tril((2 * 6378137 * np.arcsin(np.sqrt(d))).astype(np.float32))
    start = end
    end += N
    L = len(lat1) - end
    #A bunch more iterations
    while L>0:
        diff_lat = lat1[start:end,None] - lat2
        diff_lng = lng1[start:end,None] - lng2
        d = np.sin(diff_lat/2)**2 + np.cos(lat1[start:end,None])*np.cos(lat2) * np.sin(diff_lng/2)**2
        d = scipy.sparse.tril((2 * 6378137 * np.arcsin(np.sqrt(d))).astype(np.float32))
        D = scipy.sparse.vstack([D,d])
        start = end
        end += N
        L = len(lat1) - end
    #Last iteration    
    diff_lat = lat1[start:,None] - lat2
    diff_lng = lng1[start:,None] - lng2
    d = np.sin(diff_lat/2)**2 + np.cos(lat1[start:,None])*np.cos(lat2) * np.sin(diff_lng/2)**2
    d = scipy.sparse.tril((2 * 6378137 * np.arcsin(np.sqrt(d))).astype(np.float32))
    D = scipy.sparse.vstack([D,d])
        
    return D

dist = haversine(*convert_to_arrays(df, df)

【讨论】：