如何高效获取大尺寸NumPy坐标数组的对称差（含分块处理）

我有两个超大尺寸的NumPy坐标数组，格式为[[x1,y1,z1],[x2,y2,z2]...]，包含约10^9个元素。两个数组尺寸不同，坐标存在重叠（部分坐标仅存在于其中一个数组），需要快速获取所有仅存在于单个数组中的坐标。

示例代码 #

import numpy as np
array1 = np.array([[1,2,3],[2,3,4],[3,4,5],[4,5,6]])
array2 = np.array([[2,3,4],[3,4,5],[4,5,6],[5,6,7]])

array_diff = some_function_to_get_diff(array1,array2)

期望输出 #

array_diff: np.array([[1,2,3],[5,6,7]])

已知一维数组可用setdiff1d处理，但不清楚如何适配二维数组。文件大小达10GB级别，希望方案支持快速处理、分块并行。

已测试两种实现方案，性能对比如下：

方案1：基于集合的实现 #

def some_faster_function_to_get_diff(array1, array2):
    unique1 = np.unique(array1, axis=0)  # 获取array1中的唯一坐标
    unique2 = np.unique(array2, axis=0)  # 获取array2中的唯一坐标

    set1 = set(map(tuple, unique1))
    set2 = set(map(tuple, unique2))

    array1_only = np.array(list(set1 - set2))
    array2_only = np.array(list(set2 - set1))

    if array1_only.size == 0:
        return array2_only
    if array2_only.size == 0:
        return array1_only

    return np.vstack((array1_only, array2_only))

方案2：基于结构化数组的实现（由@Andrés提供） #

def symmetric_difference(arr1, arr2):        
    dtype = [('f{}'.format(i), arr1.dtype) for i in range(arr1.shape[1])]
    struct_arr1 = np.ascontiguousarray(arr1).view(dtype)
    struct_arr2 = np.ascontiguousarray(arr2).view(dtype)
    unique_to_arr1 = np.setdiff1d(struct_arr1, struct_arr2, assume_unique=True)
    unique_to_arr2 = np.setdiff1d(struct_arr2, struct_arr1,assume_unique=True)
    result = np.concatenate([unique_to_arr1, unique_to_arr2]).view(arr1.dtype).reshape(-1, arr1.shape[1])

    return result

性能测试代码及结果 #

from sys import getsizeof
import time

for i in range(2,9):
 
    no_ele=3*10**i
    arr1 = np.random.rand(no_ele, 3) * 100
    arr2 = np.random.rand(no_ele, 3) * 100

    print(f"i={i}: 内存占用 {getsizeof(arr2)/1024**3:2f}GB")
    arr2[50:60] = arr1[50:60] 
    t0=time.time()
    diff1=symmetric_difference(arr1,arr2)
    t1=time.time()
    print(f"{no_ele}: 结构化数组方案耗时 {t1-t0}")
    t0=time.time()
    diff2=some_faster_function_to_get_diff(arr1,arr2)
    t1=time.time()
    print(f"{no_ele}: 集合方案耗时 {t1-t0}")

测试输出：

i=2: 内存占用 0.000007GB
300: 结构化数组方案耗时 0.0004420280456542969
300: 集合方案耗时 0.0007300376892089844
i=3: 内存占用 0.000067GB
3000: 结构化数组方案耗时 0.0038480758666992188
3000: 集合方案耗时 0.0068819522857666016
i=4: 内存占用 0.000671GB
30000: 结构化数组方案耗时 0.052056074142456055
30000: 集合方案耗时 0.07555103302001953
i=5: 内存占用 0.006706GB
300000: 结构化数组方案耗时 0.6442217826843262
300000: 集合方案耗时 0.9859399795532227
i=6: 内存占用 0.067055GB
3000000: 结构化数组方案耗时 9.833515882492065
3000000: 集合方案耗时 12.196370840072632
i=7: 内存占用 0.670552GB
30000000: 结构化数组方案耗时 142.16755604743958
30000000: 集合方案耗时 138.87937593460083

测试结论：两种方案均可行，结构化数组方案在数据量较小时更快；当数据量达到千万级后，集合方案性能反超，整体性能相近。

当前数组无法放入内存，已知laspy包支持分块读取数据（示例代码如下），寻求基于分块读取的解决方案：

with laspy.open(las_file1) as f1, laspy.open(las_file2) as f2:
    for points1 in f1.chunk_iterator(chunk_size):
        # 处理逻辑

内容的提问来源于stack exchange，提问作者abinitio