基于线性索引在Eigen中实现矩阵映射的问题排查与解决

现有三个二维矩阵：

• shiftedData：尺寸为[tSize*TXSize,xSize]的复数双精度矩阵
• Data：尺寸为[tSize*TXSize,xSize]的复数双精度矩阵
• indices：尺寸为[tSize,xSize]的整数矩阵

需求是根据indices中的线性索引，将Data的每个子块映射到shiftedData中。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函数的实现方式及调用循环：

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) {

    // 方法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);
        }
    }

    // 方法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);
    }


    // 方法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);
    }

    // 方法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无法得到正确结果，推测原因是向extractData传递了shiftedData和Data的非连续块。请问如何正确处理该问题，实现与Matlab示例类似的功能？

错误原因分析 #

方法3、4出错的核心问题是Eigen通过seq切片得到的子块引用是非连续存储的，而Eigen::Map要求内存必须连续。当直接用Data.data()获取子块首地址并映射为向量时，会按内存物理顺序读取数据，而非子块的逻辑行列顺序，最终导致数据读取错位。

修正后的实现 #

要实现Matlab中shiftedData(indices(:)) = Data(:)的向量化效果，需先将非连续子块转换为连续存储的临时矩阵，再进行映射或索引操作。以下是几种可行方案：

方案1：拷贝为连续矩阵后执行向量化赋值

修改extractData函数，先将非连续的子块拷贝为连续矩阵，再使用向量化逻辑：

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) {
    // 将非连续子块拷贝为连续存储的临时矩阵
    Eigen::MatrixXcd Data_contiguous = Data;
    Eigen::MatrixXcd shiftedData_contiguous = shiftedData;

    // 映射为向量执行赋值
    Eigen::Map<const Eigen::VectorXcd> linData(Data_contiguous.data(), xSize*tSize);
    Eigen::Map<const Eigen::VectorXi> linIDX(indices.data(), xSize*tSize);
    Eigen::Map<Eigen::VectorXcd> linView(shiftedData_contiguous.data(), xSize*tSize);
    // 注意Matlab是1索引，Eigen是0索引，需做转换
    linView(linIDX.array() - 1) = linData;

    // 将结果写回原非连续子块
    shiftedData = shiftedData_contiguous;
}

方案2：使用Eigen 3.4+的Reshape功能

如果你的Eigen版本≥3.4，可以用reshape直接将矩阵转为向量视图，无需手动映射：

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) {
    // 先确保数据连续
    Eigen::MatrixXcd Data_contiguous = Data;
    Eigen::MatrixXcd shiftedData_contiguous = shiftedData;

    // 用reshape转为向量视图（按列优先匹配Matlab线性索引逻辑）
    auto linData = Data_contiguous.reshaped();
    auto linIDX = indices.reshaped();
    auto linView = shiftedData_contiguous.reshaped();

    linView(linIDX.array() - 1) = linData;
    shiftedData = shiftedData_contiguous;
}

方案3：优化循环（无拷贝，兼顾效率）

如果不想拷贝数据，可优化方法1的循环逻辑，利用Eigen的表达式模板减少冗余开销：

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) {
    for (int nx = 0; nx < xSize; ++nx) {
        for (int nt = 0; nt < tSize; ++nt) {
            int idx = indices(nt, nx);
            // Matlab 1索引转Eigen 0索引
            int row = (idx - 1) % tSize;
            int col = (idx - 1) / tSize;
            shiftedData(row, col) = Data(nt, nx);
        }
    }
}

关键注意事项 #

1. 索引转换：Matlab采用1-based索引，Eigen采用0-based索引，必须将indices中的值减1后再使用（对应Matlab代码中的indices+1）。
2. 内存连续性：Eigen的切片子块默认非连续，不能直接用Map或reshape，必须先拷贝为连续矩阵。
3. 存储顺序：Matlab线性索引按列优先，Eigen默认也是列优先，若需匹配行优先逻辑，可在reshaped中指定Eigen::RowMajor。

内容的提问来源于stack exchange，提问作者drakon101