c++ - 加快地面组中所有线对之间的 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，在特征 1 中的值为 4，在特征 -1 中的值为 0。 this 的值是一般实数( double )。

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

我已经使用 Eigen 库来执行我的大部分计算。为了计算 L2(欧几里得)，我使用以下观察:对于大小为 n 的两个 vector 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

In other words, we rewrite the squared norm using the dot product of the vector by themselves and subtract twice the dot product between them. From that, we just take the square and we are done. In time, I found this trick a long time ago and unfortunately I lost the reference for the author.

Anyway, this enable to write a fancy code using Eigen (in 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 进行快速计算吗？或者更好的形式来做到这一点，我只是没有弄清楚。

非常感谢您的帮助!

卡洛斯

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