加速计算一个给定集合内所有元素之间的L1距离

Question

我有一个矩阵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，在特征 2 中的值为 4，在特征 -1 中的值为 0。这些值是一般实数（双精度）。

我需要计算所有对象（所有行对）之间的多个自定义距离/不相似度。为了比较，我想计算 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
            = 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，使用Eigen进行快速计算？或者有更好的方法来做到这一点，我只是没有想出来。

非常感谢您的帮助！

卡洛斯

- an_drade

---

这并不回答你的问题，但请注意，如果你只有两个物理核心，那么你应该只启用两个线程，因为超线程会减慢CPU密集型操作。你也可以使用replicate来初始化D: D = XX.replicate(1,n) + XX.transpose().replicate(n,1); - ggael

我要冒个险...注意到你正在操作行。然而，默认情况下，Eigen矩阵是按列主序排列的（http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html）。这意味着每当你调用row()时，Eigen必须从许多不连续的内存区域中读取。如果你切换到行主序，你是否会获得更好的性能/更少的缓存未命中次数？请注意，由L2范数使用的矩阵乘法不会受到太大影响，因为底层操作通过dgemm中的'T'参数进行了优化。 - Patrick Mineault

@PatrickMineault 是的，你说得对。我确实改变了矩阵顺序以加快速度。这确实有所改善，但并不是我想要的那种程度。无论如何，感谢你的提醒。 - an_drade

如果可以使用8位值，您可以使用指令或内置函数_mm_mpsadbw_epu8。通过这种方式，您可以在9个时钟周期内执行8个8字节的绝对差总和。https://software.intel.com/sites/default/files/m/a/9/b/7/b/1000-SSE.pdf。 - Jens Munk

我无法在数学上帮助您，但我可以向您保证，在编程方面有更快的方法。对于一个10k乘10k的矩阵，考虑使用GPU可能是值得的。此外，我的经验表明，使用SIMD向量指令比仅使用openmp进行并行化要快得多。因此，无论是否使用openmp，您都应该编写代码以使用SIMD。 - Sebastian Cabot

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

- Louis Ricci

Answer 2

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

你已经正确地确定了 SSD 只需要计算两个序列之间的交叉相关作为主要计算。 然而，你可能想考虑使用 FFT 来计算这些a.b项，它将显著减少 L2 情况下所需的操作次数（但我不知道具体减少多少，这取决于 Eigen 使用的矩阵乘法算法）。如果你需要我解释一下，我可以，但我认为你也可以查一下，因为这是 FFT 的标准用法。OpenCV有一个（相当糟糕/有 bug）的模板匹配实现，这是在使用 CV_TM_SQDIFF 时所需的。

L1情况更加棘手。L1情况无法像其他情况那样被很好地分解，但它也是您可以执行的最简单的操作之一（只需进行加法和绝对值运算）。因此，许多计算架构都有并行化实现这个操作，作为指令或硬件实现函数。其他架构则有研究人员尝试找到计算此操作的最佳方法。

您可能还想了解 Intel Imaging Primitives 进行互相关，以及快速的 FFT 库，例如 FFTW 和 CUFFT。如果无法购买 Intel Imaging Primitves，您可以使用 SSE 指令 大大加速处理速度，实现几乎相同的效果。

- ldog