在scipy中计算含缺失值的成对距离

Question

我对 scipy.spatial.distance.pdist 处理缺失值（nan）的方式有点困惑。

为了确保我的矩阵维度没有搞错，先把这个问题解决掉。从文档中可以看到：

这些点在矩阵 X 中被安排成 m 个 n 维的行向量。

接下来，我们生成三个在 10 维空间中的点，并且这些点有缺失值：

numpy.random.seed(123456789)
data = numpy.random.rand(3, 10) * 5
data[data < 1.0] = numpy.nan

如果我计算这三个观察值之间的欧几里得距离：

pdist(data, "euclidean")

我得到的结果是：

array([ nan,  nan,  nan])

不过，如果我过滤掉所有有缺失值的列，我就能得到正确的距离值：

valid = [i for (i, col) in enumerate(data.T) if ~numpy.isnan(col).any()]
pdist(data[:, valid], "euclidean")

我得到的结果是：

array([ 3.35518662,  2.35481185,  3.10323893])

这样做的话，我丢掉的数据比我想要的还要多，因为我其实不需要过滤整个矩阵，只需要在比较的向量对之间进行过滤。有没有办法让 pdist 或类似的函数进行成对的掩码处理呢？

---

编辑：

由于我的完整矩阵相当大，我在这里提供的小数据集上做了一些时间测试。

1.) 使用 scipy 函数。

%timeit pdist(data, "euclidean")

10000 loops, best of 3: 24.4 µs per loop

2.) 不幸的是，目前提供的解决方案大约慢了 10 倍。

%timeit numpy.array([pdist(data[s][:, ~numpy.isnan(data[s]).any(axis=0)], "euclidean") for s in map(list, itertools.combinations(range(data.shape[0]), 2))]).ravel()

1000 loops, best of 3: 231 µs per loop

3.) 然后我测试了一下“纯” Python，结果让我感到惊喜：

from scipy.linalg import norm

%%timeit
m = data.shape[0]
dm = numpy.zeros(m * (m - 1) // 2, dtype=float)
mask = numpy.isfinite(data)
k = 0
for i in range(m - 1):
    for j in range(i + 1, m):
        curr = numpy.logical_and(mask[i], mask[j])
        u = data[i][curr]
        v = data[j][curr]
        dm[k] = norm(u - v)
        k += 1

10000 loops, best of 3: 98.9 µs per loop

所以我认为接下来的方向是把上面的代码用 Cython 优化成一个函数。

Answer 1

如果我理解得没错，你是想要计算两个向量在所有有效维度上的距离。

可惜的是，pdist 这个函数并不支持这种带有掩码的数组，所以我对你的半成品做了一些修改，以避免丢失信息。不过，这个方法并不是最有效的，也不是最容易理解的：

np.array([pdist(data[s][:, ~numpy.isnan(data[s]).any(axis=0)], "euclidean") for s in map(list, itertools.combinations(range(data.shape[0]), 2))]).ravel()

外层的数组和 ravel 只是为了让数据的形状符合你的预期。

itertools.combinations 会生成 data 数组中所有可能的成对索引。

然后我就根据这些索引来切片数据（这里必须用 list 而不是 tuple，这样才能正确切片），并像你的代码那样进行成对的 nan 过滤。

Answer 2

其实你可以考虑使用这个现成的解决方案：https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.nan_euclidean_distances.html

不过，它的一个缺点是，当数据中有缺失值的时候，添加权重会比较麻烦。