C++ Eigen转换for循环矩阵乘法计算结果错误问题

原始参考代码 #

以下是基于手动内存管理、双层循环实现的矩阵计算代码，运行结果符合预期：

static const int ny = 10; 
std::vector<double> ygl(ny+1);

double *dVGL; 
dVGL = (double*) fftw_malloc(((ny+1)*(ny+1))*sizeof(double));
memset(dVGL, 42, ((ny+1)*(ny+1))* sizeof(double));
double *dummy1; 
dummy1 = (double*) fftw_malloc(((ny+1)*(ny+1))*sizeof(double));
memset(dummy1, 42, ((ny+1)*(ny+1))* sizeof(double));

double *dummy2; 
dummy2 = (double*) fftw_malloc(((ny+1)*(ny+1))*sizeof(double));
memset(dummy2, 42, ((ny+1)*(ny+1))* sizeof(double));

for (int i = 0; i < ny+1; i++){
    for (int j = 0; j < ny+1; j++){
        ygl[j] = -1. * cos(((j) * EIGEN_PI )/ny);
        dummy1[j + ny*i] =  1./(sqrt(1-pow(ygl[j],2))); 
        dummy2[j + ny*i] = 1. * (i);
        dVGL[j + ny*i] =  dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i];
    }
}

问题Eigen代码 #

转换为Eigen实现后计算结果存在偏差，问题代码如下：

Eigen::Matrix< double, 1, ny+1> v1; 
v1.setZero();
std::iota(v1.begin(), v1.end(), 0);

Eigen::Matrix< double, ny+1, ny+1> dummy1; 
dummy1.setZero();

Eigen::Matrix< double, ny+1, ny+1> dummy2; 
dummy2.setZero();

for (int j = 0; j < ny+1; j++){
    v[j] = 1./(sqrt(1-pow(ygl[j],2)));
}
dummy1 = v.array().matrix().asDiagonal(); 
dummy2 = v1.array().matrix().asDiagonal(); 

Eigen::Matrix< double, ny+1, ny+1> dVGL; 
dVGL.setZero();
for (int i = 0; i < ny+1; i++){
    for (int j = 0; j < ny+1; j++){
        ygl[j] = -1. * cos(((j) * EIGEN_PI )/ny);
        dVGL(j + ny*i) =   sin(acos(ygl[j]) * (i)); 
    }
}

dv1 = (dummy1) * (dVGL) * (dummy2); // 结果错误

问题根因 #

核心问题共3处，和你猜测的索引步长问题直接相关：

• 扁平化索引步长错误：矩阵维度为(ny+1) × (ny+1)，行优先存储时每行步长应为ny+1，代码中用ny作为步长计算偏移j + ny*i，会导致不同行的元素内存位置重叠，本身存在逻辑隐患，仅因原始代码申请的内存空间足够大未触发崩溃。
• Eigen存储顺序不匹配：Eigen默认矩阵为列优先存储，使用单参数dVGL(index)访问元素时，按列优先顺序计算内存偏移，和你手动按行优先计算偏移的逻辑完全不一致，直接导致元素存储位置错位。你在循环内用相同的单索引写入后立刻读取，会误以为赋值正确，但实际矩阵对应(i,j)坐标位置的值完全错误。
• 对角矩阵乘法顺序错误：对角矩阵左乘是对矩阵做行缩放（每行乘对应对角元），右乘是对矩阵做列缩放（每列乘对应对角元）。原始逻辑中dummy1对应列缩放因子（和j相关），应该右乘；dummy2对应行缩放因子（和i相关），应该左乘，代码中的乘法顺序完全颠倒。此外原始代码中重复在i循环内计算固定值ygl[j]，属于无意义性能损耗。

修复方案 #

• 所有矩阵元素访问统一使用Eigen双索引(i,j)，明确指定行号i、列号j，完全规避扁平化索引的步长、存储顺序问题。
• 将和j相关的ygl、缩放因子v的计算移到i循环外层，避免重复计算。
• 调整对角矩阵乘法顺序：行缩放对角阵左乘，列缩放对角阵右乘。
• 如果不需要保留对角矩阵乘法的写法，也可以直接用Eigen广播机制构造全量缩放矩阵，通过逐元素乘积得到结果，逻辑更直观。

修正后代码（对角矩阵乘法版本，和原始计算逻辑完全对齐） #

static const int ny = 10;
std::vector<double> ygl(ny+1);
// 提前计算固定值
Eigen::Matrix<double, ny+1, 1> col_scale(ny+1);
for (int j = 0; j < ny+1; j++){
    ygl[j] = -1. * cos((j * EIGEN_PI)/ny);
    col_scale[j] = 1./(sqrt(1 - pow(ygl[j],2)));
}

// 构造行缩放对角矩阵：对角元为行号i，左乘实现行缩放
Eigen::Matrix<double, ny+1, 1> row_scale(ny+1);
std::iota(row_scale.begin(), row_scale.end(), 0);
auto D_row = row_scale.asDiagonal();
// 构造列缩放对角矩阵：对角元为列对应缩放因子，右乘实现列缩放
auto D_col = col_scale.asDiagonal();

// 构造正弦值矩阵，用双索引赋值
Eigen::Matrix<double, ny+1, ny+1> sin_mat;
for (int i = 0; i < ny+1; i++){
    for (int j = 0; j < ny+1; j++){
        sin_mat(i,j) = sin(acos(ygl[j]) * i);
    }
}

// 按正确顺序计算最终结果
Eigen::Matrix<double, ny+1, ny+1> dv1 = D_row * sin_mat * D_col;

注：如果需要和原始带步长错误的代码结果逐比特对齐，只需要调整sin_mat的赋值偏移即可，但原始步长写法属于未定义行为，存在内存越界风险，不建议保留。

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