使用Eigen加速稀疏矩阵乘法的方法及与MATLAB性能差异原因

我需要使用Eigen进行大规模稀疏矩阵的乘法、加法或减法运算，但部分场景下其运算速度远慢于MATLAB。例如以下代码：

Eigen::SparseMatrix<double> JPSlice = JP.block(0, 3, JP.rows(), JP.cols()-3);
Eigen::SparseMatrix<double> NonZeroJD = RemoveSpMatrixZeroCols(JD);
JPSlice.makeCompressed();
NonZeroJD.makeCompressed();
Eigen::SparseMatrix<double> U = JPSlice.transpose() * JPSlice;
Eigen::SparseMatrix<double> V = NonZeroJD.transpose() * NonZeroJD;
Eigen::SparseMatrix<double> InvV = InverseDiagonalSpMatrix(V);
InvV.makeCompressed();
Eigen::SparseMatrix<double> VM1 = JPSlice.transpose() * NonZeroJD;
Eigen::SparseMatrix<double> Block = VM1 * InvV;
Eigen::SparseMatrix<double> VM = Block * VM1.transpose();

JPSlice和NonZeroJD的维度分别为26400000×801和26400000×529530。

使用Eigen计算VM1和VM分别耗时5.8秒和25秒，而MATLAB仅需约0.8秒和0.7秒，但计算V的耗时两者相近（约3秒）。

我在Mac上通过Clang和CMake编译代码，编译器设置如下：

set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -O3 -march=native -mavx -mfma -fopenmp-simd -fopenmp=libomp -DNDEBUG -DEIGEN_USE_MKL_ALL")

请问如何加速这些稀疏矩阵运算？为何Eigen与MATLAB的运算速度存在如此大的差异？

速度差异的核心原因 #

• 后端实现差异：MATLAB底层稀疏矩阵乘法依赖高度优化的Intel MKL稀疏BLAS或自有优化库，针对大规模稀疏矩阵的转置乘场景做了深度适配；Eigen默认的稀疏乘法虽支持MKL，但部分场景的调度逻辑、内存布局适配不如MATLAB紧密。
• 维度特性利用不足：你的场景中JPSlice是超大型行稀疏矩阵（2640万行×801列），转置后变成列稀疏的小宽矩阵（801×2640万），和NonZeroJD相乘时，MATLAB会自动识别这种特殊结构，采用缓存友好的分块计算策略；Eigen通用稀疏乘法逻辑未针对性优化这类极端维度组合。
• 并行调度效率差异：Mac上Clang+OpenMP组合在稀疏矩阵运算的并行负载均衡上，不如MATLAB针对MKL做的线程调度优化，非均匀稀疏结构下线程负载差异会拉低整体效率。

加速Eigen运算的具体方案 #

1. 强制启用MKL稀疏后端并优化调用 #

虽然定义了DEIGEN_USE_MKL_ALL，但Eigen部分场景可能未自动切换到MKL后端，显式调用MKL稀疏乘法接口可大幅提速：

#include <mkl_spblas.h>

// 转换Eigen矩阵到MKL CSR格式（MKL对CSR支持更高效）
mkl_sparse_matrix mkl_JPSlice_T, mkl_NonZeroJD, mkl_VM1;
struct matrix_descr descr = {SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT};

// 转置JPSlice并转为行优先存储
Eigen::SparseMatrix<double, Eigen::RowMajor> JPSlice_T = JPSlice.transpose();
mkl_sparse_d_create_csr(&mkl_JPSlice_T, SPARSE_INDEX_BASE_ZERO, JPSlice_T.rows(), JPSlice_T.cols(),
                        JPSlice_T.outerIndexPtr(), JPSlice_T.outerIndexPtr()+1,
                        JPSlice_T.innerIndexPtr(), JPSlice_T.valuePtr());

mkl_sparse_d_create_csr(&mkl_NonZeroJD, SPARSE_INDEX_BASE_ZERO, NonZeroJD.rows(), NonZeroJD.cols(),
                        NonZeroJD.outerIndexPtr(), NonZeroJD.outerIndexPtr()+1,
                        NonZeroJD.innerIndexPtr(), NonZeroJD.valuePtr());

// 执行乘法：VM1 = JPSlice_T * NonZeroJD
mkl_sparse_d_mm(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, mkl_JPSlice_T, descr,
                SPARSE_OPERATION_NON_TRANSPOSE, mkl_NonZeroJD, 0.0, &mkl_VM1);

// 转换回Eigen格式
Eigen::SparseMatrix<double> VM1;
VM1.resize(JPSlice_T.rows(), NonZeroJD.cols());
VM1.makeCompressed();
mkl_sparse_d_export_csr(mkl_VM1, &descr, &VM1.rows(), &VM1.cols(),
                        &VM1.outerIndexPtr(), &VM1.innerIndexPtr(), &VM1.valuePtr());

// 释放MKL资源
mkl_sparse_destroy(mkl_JPSlice_T);
mkl_sparse_destroy(mkl_NonZeroJD);
mkl_sparse_destroy(mkl_VM1);

2. 利用对角矩阵特性重构计算 #

InvV是对角矩阵，Block = VM1 * InvV无需完整稀疏乘法，直接遍历VM1非零元素做逐列缩放即可：

Eigen::SparseMatrix<double> Block = VM1;
// 提前提取InvV的对角元素到向量中
Eigen::VectorXd inv_v_diag(InvV.rows());
for (int i=0; i<InvV.rows(); ++i) {
    inv_v_diag[i] = InvV.coeff(i,i);
}

// 对VM1的每一列元素乘以对应对角值
for (int k=0; k<Block.outerSize(); ++k) {
    for (Eigen::SparseMatrix<double>::InnerIterator it(Block, k); it; ++it) {
        it.valueRef() *= inv_v_diag[it.row()];
    }
}

此操作时间复杂度为O(VM1非零元素数)，远低于稀疏矩阵乘法，能大幅节省时间。

3. 优化编译与并行设置 #

• 确保OpenMP正确链接，在CMake中添加：

target_link_libraries(your_target_name PRIVATE omp)

• 替换-fopenmp-simd为-ffast-math，配合-march=native最大化硬件浮点优化能力。
• 尝试使用Homebrew安装的GCC替代Clang，GCC对OpenMP和MKL的兼容性在稀疏运算场景下表现更好。

4. 调整Eigen稀疏存储格式 #

将涉及转置乘法的矩阵改为**行优先（RowMajor）**存储，适配MKL的CSR格式偏好，减少格式转换开销。

内容的提问来源于stack exchange，提问作者Yingyu.Wang