加快地面组中所有线对之间的 L1 距离答案

【问题标题】：Speeding up the L1 distance between all pairs in a ground set加快地面组中所有线对之间的 L1 距离
【发布时间】：2015-08-07 10:19:11
【问题描述】：

我有一个矩阵 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（欧几里得）距离。

我已经使用 Eigen 库来执行我的大部分计算。为了计算 L2（欧几里得），我使用以下观察：对于大小为 n 的两个向量 a 和 b，我们有：

||a - b||^2 = (a_1 - b_1)^2 + (a_2 - b_2)^2 + ... +(a_n - b_n)^2 = a_1^2 + b_1^2 - 2 a_1 b_1 + a_2^2 + b_2^2 - 2 a_2 b_2 + ... + a_n^2 + b_n^2 - 2 a_n b_n = 一个。一个 + 乙。 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 分钟（2 核 / 4 线程）内计算所有对象/线对的 L2 距离，我个人认为这对我的目的来说是一个很好的性能。 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 以使用 Eigen 进行快速计算吗？或者更好的形式来做到这一点，我只是没有弄清楚。

非常感谢您的帮助！

卡洛斯

【问题讨论】：

这不能回答您的问题，但请注意，如果您只有 2 个物理内核，那么您应该只启用 2 个线程，因为超线程会减慢 CPU 密集型操作。您还可以使用复制初始化 D：D = XX.replicate(1,n) + XX.transpose().replicate(n,1);
我要在这里冒个险...请注意，您正在操作行。但是，默认情况下，特征矩阵按列优先顺序排列 (eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html)。这意味着每当您调用 row() 时，Eigen 都必须从大量不连续的内存区域中读取。如果切换到行优先顺序，您是否会获得更好的性能/更少的缓存未命中？请注意，L2 范数使用的矩阵乘法不受此影响，因为基础操作通过 dgemm 中的“T”参数针对两个订单进行了优化
@PatrickMineault 是的，你是对的。确实，我已经更改了矩阵顺序以加快一点速度。它有很大的不同，但不是我正在寻找的那个。无论如何，谢谢你的通知。
如果可以使用 8 位值，您可以使用指令或内在 _mm_mpsadbw_epu8。这样，就可以在 9 个时钟周期内做 8 个 8 字节的绝对差之和。 software.intel.com/sites/default/files/m/a/9/b/7/b/1000-SSE.pdf.
我无法在数学上为您提供帮助，但我可以向您保证，有更快的方法可以以编程方式进行。对于 10k x 10k 矩阵，GPU 可能值得考虑。此外，我自己的经验表明，使用 SIMD 矢量指令比仅使用 openmp 进行并行化要快得多。因此，无论是否使用 openmp，您都应该编写代码以使用 SIMD

标签： c++ algorithm matrix parallel-processing eigen

【解决方案1】：

让我们同时计算两个距离。他们只真正共享行差异（虽然两者都可能是绝对差异，但欧几里德距离使用平方，这并不是真正等效的）。所以现在我们只循环 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()
    }
}

编辑另一种更占用内存/未经测试的方法可能是每次迭代计算一个完整的距离行（不知道这是否会有所改进）

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();
}

【讨论】：

【解决方案2】：

这是图像处理中两个非常常见的操作。第一个是Sum of Squared Differences (SSD)，第二个是Sum of Absolute Differences (SAD)。

您已正确确定 SSD 只需要一个来计算两个系列之间的 cross-correlation 作为主要计算。但是，您可能需要考虑使用 FFT 来计算这些 a.b 项，它将显着减少 L2 情况所需的操作数量（但是我不知道多少，这取决于什么矩阵- Eigen 使用的矩阵乘法算法。）如果您需要我解释这一点，我可以，but I figure you can also look it up as its a standard use of FFTs。 OpenCV 有一个（相当糟糕/错误的）模板匹配实现，这是您在使用 CV_TM_SQDIFF 时想要的。

L1 案例比较棘手。 L1 案例不能很好地分解，但它也是您可以执行的最简单的操作之一（只是加法和绝对值。）因此，a lot of computation architectures have parallelized implementations 将其作为指令或硬件实现的功能。其他架构有researchers experimenting with the best way to compute this.

您可能还想查看 Intel Imaging Primitives 的互相关，以及快速 FFT 库，例如 FFTW 和 CUFFT。如果您买不起 Intel Imaging Primitves，您可以使用SSE instructions 大大加快您的处理速度，达到几乎相同的效果。

【讨论】：