生成稀疏的成对距离矩阵避免内存错误

Question

我有一些地理空间数据，我试图在点之间产生成对的哈弗线距离。不幸的是，大约有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”的数据帧。当我用大约50k点运行这段代码时，我得到一个内存错误：

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

是否有办法确保以（下三角或上三角）稀疏矩阵作为预期输出执行计算

Answer 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)