Ускорение расстояния L1 между всеми парами в наземном наборе

У меня есть матрица NxM (обычно 10k X 10k элементов), описывающая набор оснований. Каждая строка представляет объект и каждый столбец конкретную функцию. Например, в матрице

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

объект x1 имеет значение 0 в функции 1, значение 4 в функции 1 и значение 0 в элементе -1. Значения этого представляют собой общие действительные числа (double).

Мне нужно вычислить несколько пользовательских расстояний/различий между всеми парами объектов (все пары линий). Для сравнения, я хочу вычислить L1 (manhattan) и L2 (евклидовых) расстояний.

У меня есть библиотека Eigen для выполнения большинства моих вычислений. Для вычисления L2 (евклидова) я использую следующее наблюдение: для двух векторов a и b размера n мы имеем:

||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 (в С++):

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 мы можем вычислить расстояние L2 для всех пар объектов/линий менее чем за 1 минуту (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? Или лучшая форма для этого, и я просто не понял.

Большое спасибо за вашу помощь!

Карлос

Ответы

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