问题描述

设

(; 表示纵向连接) 和

, 计算矩阵

中每一个行向量和矩阵

中每一个行向量的平方欧氏距离 (pairwise squared Euclidean distance), 即计算:

(这是一个

矩阵).

这个计算在度量学习, 图像检索, 行人重识别等算法的性能评估中有着广泛的应用.

公式转化

在 NumPy 中直接利用上述原式来计算两个矩阵的成对平方欧氏距离, 要显式地使用二重循环, 而在 Python 中循环的效率是相当低下的. 如果想提高计算效率, 最好是利用 NumPy 的特性将原式转化为数组/矩阵运算. 下面就尝试进行这种转化.

先将原式展开为:

下面逐项地化简或转化为数组/矩阵运算的形式:

式中,

表示按元素积 (element-wise product), 又称为 Hadamard 积;

表示维的全1向量 (all-ones vector), 余者类推. 上式中

的作用是计算

每行元素的和, 返回一个列向量;

的作用类似于 NumPy 中的广播机制, 在这里是将一个列向量扩展为一个矩阵, 矩阵的每一列都是相同的.

所以:

上述转化式中出现了

(矩阵乘) , 矩阵乘在 NumPy 等很多库中都有高效的实现, 对代码的优化是有好处的.

特别地, 当

时, 原式等于

, 注意到第一项和第二项互为转置. 当

且

(即

和

的每一个行向量的范数均为1时), 原式等于

,

是

全1矩阵.

代码实现

sklearn 中已经包含了用 NumPy 实现的计算 "两个矩阵的成对平方欧氏距离" 的函数 (sklearn.metrics.euclidean_distances

import numpy as np

def euclidean_distances(x, y, squared=True):

"""Compute pairwise (squared) Euclidean distances."""

assert isinstance(x, np.ndarray) and x.ndim == 2

assert isinstance(y, np.ndarray) and y.ndim == 2

assert x.shape[1] == y.shape[1]

x_square = np.sum(x*x, axis=1, keepdims=True)

if x is y:

y_square = x_square.T

else:

y_square = np.sum(y*y, axis=1, keepdims=True).T

distances = np.dot(x, y.T)

# use inplace operation to accelerate

distances *= -2

distances += x_square

distances += y_square

# result maybe less than 0 due to floating point rounding errors.

np.maximum(distances, 0, distances)

if x is y:

# Ensure that distances between vectors and themselves are set to 0.0.

# This may not be the case due to floating point rounding errors.

distances.flat[::distances.shape[0] + 1] = 0.0

if not squared:

np.sqrt(distances, distances)

return distances

如果想进一步加速, 可以将

x_square = np.sum(x*x, axis=1, keepdims=True)

替换为

x_square = np.expand_dims(np.einsum('ij,ij->i', x, x), axis=1)

以及将

y_square = np.sum(y*y, axis=1, keepdims=True).T

替换为

y_square = np.expand_dims(np.einsum('ij,ij->i', y, y), axis=0)

使用 np.einsum 的好处是不会产生一个和 x 或 y 同样形状的临时数组 (x*x 或 y*y 会产生一个和 x 或 y 同样形状的临时数组).

PyTorch 中也包含了计算 "两个矩阵的成对平方欧氏距离" 的函数

另外上述的转化公式也可以用在其他 Python 框架 (如 TensorFlow) 或其他语言中, 这里就不展开叙述了.

版权声明

版权声明：自由分享，保持署名-非商业用途-非衍生，知识共享3.0协议。

如果你对本文有疑问或建议，欢迎留言！转载请保留版权声明！

参考