两个矩阵的 C++ 高效元素匹配

Question

我正在使用 C++ 和 Eigen 库来完成匹配任务。但是，结果非常低效。下面将描述任务和我当前的解决方案。

任务

该任务涉及两个相同大小的大地图，mapL 和 mapR。匹配过程按行进行。例如，mapL 的r-th 行中的元素将与mapR 的r-th 行中的元素匹配。两个匹配元素的值可能不完全相同，所以我需要在mapR中寻找差异最小的元素。

为了计算差异，我执行以下步骤：

1. 在mapL 的r-th 行中提取一个元素。
2. 提取mapR 的r-th 行。
3. mapR 的行减去mapL 的元素）。
4. 对mapL 中的每个元素重复上述步骤。

示例

一个例子将有助于理解上述过程。假设我们有

mapL = [ 1  2  3  4, 5; 
         6  7  8  9 10;
        11 12 13 14 15];
mapR = [ 5  4  3  2, 1; 
        10  9  8  7  6;
        15 14 13 12 11];

对于mapL中的第一个元素，我想获取一个向量difference用于后续的一些计算

difference = [5 4 3 2 1] - [1] = [4 3 2 1 0];

当前解决方案

在构建 C++ 代码时，我尝试使用 .replicate 和 .colwise 对过程进行矢量化。下面是一个示例代码：

int nRow = 1800;
int nCol = 2000;

// Setup two random maps for illustration purposes
MatrixXf mapL = MatrixXf::Random(nRow, nCol);
MatrixXf mapR = MatrixXf::Random(nRow, nCol);

// Initialize some parameters
VectorXf rowL;
MatrixXf rowR(1, nCol);
MatrixXf repeatR;
MatrixXf difference;

// Loop through every row
for (int r = 0; r < nRow; r++) {

    // Extract a row from mapL and mapR
    rowL = mapL.row(r);
    rowR = mapR.row(r);

    // Repeat rowR along the col direction
    repeatR = rowR.replicate(nCol, 1);

    // Compute the difference by .colwise
    difference = repeatR.colwise()-=rowL;

    /* 
    Some other codes after the difference is computed
    */ 
}

问题

上面的代码在我的机器（VS 2017，发布模式，x64）上需要大约 61 秒来遍历每一行。我在 MATLAB 中实现了类似的代码，只需要大约 6 秒。有什么方法可以提高 C++ 代码的效率吗？我对 C++ 很陌生，因此如果我遗漏了任何重要的概念，请告诉我。非常感谢！

编辑

经过更多的测试，我找到了一种更有效的方法来完成任务。分析我的原始代码显示.replicate 和.colwise 操作消耗了大量的计算能力（~90%）。使用嵌套循环需要

int nRow = 1800;
int nCol = 2000;

MatrixXf mapL = MatrixXf::Random(nRow, nCol);
MatrixXf mapR = MatrixXf::Random(nRow, nCol);

VectorXf rowL;
ArrayXf rowR(1, nCol);
ArrayXf difference;

for (int r = 0; r < nRow; r++) {
    rowL = mapL.row(r);
    rowR = mapR.row(r).array();
    for (int c = 0; c < nCol; c++) {
        difference = rowR - rowL(c);
    }
}

Answer 1

两个建议：

• 默认情况下，Eigen 以列优先顺序存储矩阵。通过定义

  来代替使用行主要顺序可能会有所帮助

  using RowMajorMatrixXf = Matrix<float, Dynamic, Dynamic, RowMajor>;
  

  并使用它而不是 MatrixXf 以便行在内存中是连续的。

• .row(r) 方法返回一个自身不复制内存的表达式对象。但是将其分配给 VectorXf 或 ArrayXf 会生成一个副本。您可以通过保留 decltype(mapL)::RowXpr rowL = mapL.row(r) 之类的表达式对象来消除副本。