加速接地组中所有对之间的L1距离

Speeding up the L1 distance between all pairs in a ground set

本文关键字：之间 L1 距离 加速      更新时间：2023-10-16

我有一个矩阵NxM(通常是10k X 10k元素)来描述一个接地集。每一行代表一个对象，每一列代表一个特定的特征。例如，在矩阵

中
   f1 f2 f3
x1 0  4  -1
x2 1  0  5
x3 4  0  0
x4 0  1  0

对象x1在特征1中值为0，在特征1中值为4，在特征-1中值为0。this的值是一般实数(双精度)。

我必须计算所有对对象(所有对线)之间的几个自定义距离/不相似性。为了比较，我想计算L1(曼哈顿)和L2(欧氏)距离

我已经使用特征库来执行我的大部分计算。为了计算L2(欧几里得)，我使用以下观察:对于大小n的两个向量a和b，我们有:<>之前a - b | | | | ^ 2 = (a_1 - b_1) ^ 2 + (a₂(b_2) ^ 2 +…+(a_n - b_n)^2= a_1 ^ 2 + b_1 ^ 2 - 2 a₁b_1 + a₂^ 2 + b_2 ^ 2 - 2 a₂b_2 +……+ a_n^2 + b_n^2 - 2a_n b_n= a。A + b。B - 2ab之前

换句话说，我们用向量本身的点积重写平方范数然后减去它们之间的点积的两倍。然后取平方，就搞定了。我很久以前就发现了这个技巧，不幸的是我丢失了作者的参考资料。

无论如何，这使得使用Eigen(在c++中)编写奇特的代码成为可能:

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix, XX, D;
// Load matrix here, for example
// matrix << 0, 4, -1,
//           1, 0,  5,
//           4, 0,  0,
//           0, 1,  0;
const auto N = matrix.rows();
XX.resize(N, 1);
D.resize(N, N);
XX = matrix.array().square().rowwise().sum();
D.noalias() = XX * Eigen::MatrixXd::Ones(1, N) +
              Eigen::MatrixXd::Ones(N, 1) * XX.transpose();
D -= 2 * matrix * matrix.transpose();
D = D.cwiseSqrt();

对于矩阵10k X 10k，我们能够在不到1分钟的时间内计算所有对对象/线的L2距离(2核/4线程)，我个人认为这对于我的目的来说是一个很好的性能。Eigen能够组合这些操作，并使用几个低级/高级优化来执行这些计算。在这种情况下，Eigen使用所有可用的内核(当然，我们可以对其进行调整)。

然而，我仍然需要计算L1距离，但我无法找到一个好的代数形式来与Eigen一起使用，并获得良好的性能。到目前为止，我有以下内容:

const auto N = matrix.rows();
for(long i = 0; i < N - 1; ++i) {
    const auto &row = matrix.row(i);
    #ifdef _OPENMP
    #pragma omp parallel for shared(row)
    #endif
    for(long j = i + 1; j < N; ++j) {
        distance(i, j) = (row - matrix.row(j)).lpNorm<1>();
    }
}

但是这需要更长的时间:对于相同的10k X 10k矩阵，此代码使用3.5分钟，考虑到L1和L2在其原始形式中非常接近，这要糟糕得多:

L1(a, b) = sum_i |a_i - b_i|
L2(a, b) = sqrt(sum_i |a_i - b_i|^2)

有什么想法，如何改变L1使用快速计算的特征?或者有更好的形式，我就是想不出来。

非常感谢你的帮助!

卡洛斯

让我们同时计算这两个距离。它们实际上只共享行差(虽然两者都可以是绝对差，但欧几里得距离使用的是平方，这并不等效)。现在我们只遍历n^2行

const auto N = matrix.rows();
for(long i = 0; i < N - 1; ++i) {
    const auto &row = matrix.row(i);
    #ifdef _OPENMP
    #pragma omp parallel for shared(row)
    #endif
    for(long j = i + 1; j < N; ++j) {
        const auto &rowDiff = row - matrix.row(j);
        distanceL1(i, j) = rowDiff.cwiseAbs().sum(); // or .lpNorm<1>(); if it's faster
        distanceL2(i, j) = rowDiff.norm()
    }
}

EDIT另一个内存更密集/未经测试的方法可能是每次迭代计算整个距离行(不知道这是否会是一个改进)

const auto N = matrix.rows();
#ifdef _OPENMP
#pragma omp parallel for shared(matrix)
#endif
for(long i = 0; i < N - 1; ++i) {
    const auto &row = matrix.row(i);
    // you could use matrix.block(i,j,k,l) to cut down on the number of unnecessary operations
    const auto &mat = matrix.rowwise() - row;
    distanceL1(i) = mat.cwiseAbs().sum().transpose();
    distanceL2(i) = mat.rowwise().norm().transpose();
}

这是图像处理中非常常见的两个操作。第一个是平方差和(SSD)，第二个是绝对差和(SAD)。

您已经正确地确定SSD只需要一个来计算两个系列之间的相互关系作为主要计算。然而，你可能想要考虑使用FFT来计算这些a.b项，它将大大减少L2情况下所需的操作次数(然而我不知道多少，这取决于什么矩阵-矩阵乘法算法Eigen使用)。如果你需要我解释一下，我可以，但我想你也可以把它作为fft的标准使用来查找。OpenCV有一个(非常糟糕/有bug的)模板匹配实现，这是你在使用CV_TM_SQDIFF时想要的。

L1的情况比较棘手。L1的情况不能很好地分解，但它也是您可以做的最简单的操作之一(只需添加&绝对价值。)因此，许多计算体系结构都将其并行实现为指令或硬件实现函数。其他架构的研究人员正在尝试计算这个的最佳方法。

您可能还想查看Intel Imaging Primitives的交叉相关和快速FFT库，如FFTW &CUFFT。如果您买不起Intel Imaging primitives，您可以使用SSE指令来大大加快处理速度，达到几乎相同的效果。

相关文章：