使用线性索引列表在特征值中的输入和输出矩阵之间进行映射

Question

假设我有 3 个 2D 矩阵：

shiftedData - 复数双倍大小 [tSize*TXSize,xSize]
数据 - 大小为复数双精度 [tSize*TXSize,xSize]
indices - 大小 [tSize,xSize] 的整数

我想根据矩阵索引中定义的线性索引将数据的每个子块映射到移位数据中。在 Matlab 中，我可能会执行以下操作：

shiftedData = zeros(tSize*TXSize,xSize);
for nT = 0:TXSize-1
shiftedData( (nT*tSize+1):((nT+1)*tSize),:) = ...
            extractData(Data((nT*tSize+1):((nT+1)*tSize),:), ...
            indices(:,(nR*xSize+1):((nR+1)*xSize))+1,xSize,tSize);
end

function [shiftedData] = extractData(Data,indices,xSize,tSize)

    shiftedData(indices(:)) = Data(:);

end

我正在尝试在 Eigen 中实现类似的东西。我已经为类似的 extractData 函数编写了 4 种不同的方法，并在下面的伪代码中展示了它们以及从 for 循环调用它们的行：

Eigen::MatrixXcd shiftedData = Eigen::MatrixXcd::Zero(tSize*TXSize,xSize);
for (int nT = 0; nT < TXSize; nT++){
   extractData(shiftedData(Eigen::seq(nT*tSize,(nT+1)*tSize-1),Eigen::all), 
   Data(Eigen::seq(nT*tSize,(nT+1)*tSize-1),Eigen::all), indices, xSize, tSize);
}


void extractData(Eigen::Ref<Eigen::MatrixXcd> shiftedData,
            const Eigen::Ref<const Eigen::MatrixXcd>& Data, 
            const Eigen::Ref<const Eigen::MatrixXi>& indices, 
            const long int& xSize, const long int& tSize) {

    // Approach 1
    for (int nt = 0; nt < tSize; nt++){
        for (int nx = 0; nx < xSize; nx++){
            int row = indices(nt,nx) % tSize;
            int col = indices(nt,nx) / tSize;
            shiftedData(row,col) = Data(nt,nx);
        }
    }

    // Approach 2
    for (int i = 0; i < tSize*xSize; i++) {
        int nt = i % tSize;
        int nx = i / tSize;
        int row = indices(nt,nx) % tSize;
        int col = indices(nt,nx) / tSize;
        shiftedData(row,col) = Data(nt,nx);
    }


    // Approach 3
    Eigen::Map<const Eigen::VectorXcd> linData(Data.data(), xSize*tSize);
    Eigen::Map<const Eigen::VectorXi> linIDX(indices.data(), xSize*tSize);
    for (int i = 0; i < xSize*tSize; i++){
        int col = linIDX(i) / tSize;
        int row = linIDX(i) % tSize;
        shiftedData(row,col) = linData(i);
    }

    // Approach 4
    Eigen::Map<const Eigen::VectorXcd> linData(Data.data(), xSize*tSize);
    Eigen::Map<const Eigen::VectorXi> linIDX(indices.data(), xSize*tSize);
    Eigen::Map<Eigen::VectorXcd> linView(shiftedData.data(), xSize*tSize);
    linView(linIDX) = linData;

}

方法 1 和 2 生成正确的输出。然而，方法 3 和 4 却没有。我相信这个问题与我将不连续的

shiftedData

和

Data

块传递给

extractData

这一事实有关。我怎样才能正确地解释这一点并执行与 Eigen 中的 Matlab 示例更类似的操作？

Answer 1

我将重点改进您的第一种方法，因为它似乎是最简单的。正如您所怀疑的，

Map

无法工作。您可以使用

Reshape

但这只会隐藏 Eigen 方法内的除法和模数。让我们尝试一些不同的方法。

首先，进行一些小的清理。我将完全忽略

TXSize

上的外循环，并专注于

extractData

函数。界面有两个小问题：

您将两个维度作为
```
const long int&
```
传递，这是毫无意义的，创建对简单整数的引用只会带来性能损失的风险（由于别名而导致冗余重新加载的潜在风险）。只需按值传递它们即可
```
long int
```
本身并不是最好的类型。从您链接的 CodeReview 中我知道您在 Windows 上操作，这意味着
```
long int
```
与
```
int
```
基本相同。最好使用在所有平台上尺寸正确的类型。当您使用 Eigen 时，请使用
```
Eigen::Index
```
，默认情况下为
```
std::ptrdiff_t
```
，并且可能是您输入
```
long int
```
时所期望的结果。

回到主菜。为了进行测试，我使用 CodeReview 页面中的尺寸，即

tSize = 2688; xSize = 192

迭代顺序

对于方法 1，你的循环是

    for (int nt = 0; nt < tSize; nt++){
        for (int nx = 0; nx < xSize; nx++){
            ...
            ... = indices(nt,nx) ...;
            ... = Data(nt,nx);
        }
    }

注意如何迭代内循环中的列和外循环中的行。 Eigen 默认使用列主格式，这使得此顺序在数据局部性较差的情况下不是最优的。

假设

indices

中的所有条目都是唯一的，我们可以交换循环头的顺序。在我的系统上，这给我带来了 50% 的改进。

void extractData2(Eigen::Ref<Eigen::MatrixXcd> shiftedData,
            const Eigen::Ref<const Eigen::MatrixXcd>& Data,
            const Eigen::Ref<const Eigen::MatrixXi>& indices,
            Eigen::Index xSize, Eigen::Index tSize) {

    for (Eigen::Index nx = 0; nx < xSize; nx++){
        for (Eigen::Index nt = 0; nt < tSize; nt++){
            Eigen::Index row = indices(nt,nx) % tSize;
            Eigen::Index col = indices(nt,nx) / tSize;
            shiftedData(row,col) = Data(nt,nx);
        }
    }
}

避免冗余划分

您可以对多个矩阵重复使用相同的索引。这意味着您要重复进行整数除法。在许多 CPU 上，整数除法非常慢，每条指令的吞吐量为 6-14 个周期。也许值得缓存它。 void extractData(Eigen::Ref<Eigen::MatrixXcd> shiftedData, const Eigen::Ref<const Eigen::MatrixXcd>& Data, const Eigen::Ref<const Eigen::Array2Xi>& indices, Eigen::Index xSize, Eigen::Index tSize) { for (Eigen::Index nx = 0; nx < xSize; nx++){ for (Eigen::Index nt = 0; nt < tSize; nt++){ Eigen::Index flatIdx = nx * tSize + nt; const auto outIndex = indices.col(flatIdx); shiftedData(outIndex.x(), outIndex.y()) = Data(nt, nx); } } } ... Eigen::Array2Xi splitIndices(2, xSize * tSize); for(Eigen::Index i = 0; i < xSize; ++i) { for(Eigen::Index j = 0; j < tSize; ++j) { auto out = splitIndices.col(i * tSize + j); const int flatIdx = indices(j, i); out.x() = flatIdx % tSize; out.y() = flatIdx / tSize; } } for(int i = 0; i < TXSize; ++i) extractData(ShiftedData, Data, splitIndices, xSize, tSize);

但是，在我的系统上我没有看到任何好处。我假设

shiftedData

上的随机内存访问成本占主导地位。

比起商店更喜欢随机加载

一般来说，CPU 比存储更擅长处理随机内存负载。原因是随机位置上的存储仍然必须将整个缓存行拉入 L1 缓存，然后可能再次将其从缓存中逐出，而无需写入整个缓存行。

我们可以反转索引并将其与上面的缓存索引计算相结合。

void extractData(Eigen::Ref<Eigen::MatrixXcd> shiftedData, const Eigen::Ref<const Eigen::MatrixXcd>& Data, const Eigen::Ref<const Eigen::Array2Xi>& indices, Eigen::Index xSize, Eigen::Index tSize) { for (Eigen::Index nx = 0; nx < xSize; nx++){ for (Eigen::Index nt = 0; nt < tSize; nt++){ Eigen::Index flatIdx = nx * tSize + nt; const auto inIndex = indices.col(flatIdx); shiftedData(nt, nx) = Data(inIndex.x(), inIndex.y()); } } } ... Eigen::Array2Xi inverseIndices(2, xSize * tSize); for(Eigen::Index i = 0; i < xSize; ++i) { for(Eigen::Index j = 0; j < tSize; ++j) { int flatIdx = indices(j, i); Eigen::Index row = flatIdx % rows; Eigen::Index col = flatIdx / rows; auto out = inverseIndices.col(col * tSize + row); out.x() = j; out.y() = i; } } for(int i = 0; i < TXSize; ++i) extractData(ShiftedData, Data, inverseIndices, cols, rows);

在我的特定平台（TigerLake）上，我只看到了最小的积极影响。然而，当我将矩阵的大小减小到适合 L2 缓存 (128 x 192) 的大小时，与上面的缓存索引相比，性能提高了 18%，并且缓存索引比不缓存除法有 7% 的优势。

使用线性索引列表在特征值中的输入和输出矩阵之间进行映射

问题描述投票：0回答：1

1个回答

迭代顺序

避免冗余划分

一般来说，CPU 比存储更擅长处理随机内存负载。原因是随机位置上的存储仍然必须将整个缓存行拉入 L1 缓存，然后可能再次将其从缓存中逐出，而无需写入整个缓存行。

最新问题

使用线性索引列表在特征值中的输入和输出矩阵之间进行映射

问题描述 投票：0回答：1

1个回答

迭代顺序

避免冗余划分

一般来说，CPU 比存储更擅长处理随机内存负载。原因是随机位置上的存储仍然必须将整个缓存行拉入 L1 缓存，然后可能再次将其从缓存中逐出，而无需写入整个缓存行。

最新问题

问题描述投票：0回答：1