如何在Matlab与Eigen(C++)间传递稀疏数组及实现稀疏数组乘法？

一、稀疏数组在Matlab与Eigen间的双向传递 #

1. 从Matlab传递稀疏数组到Eigen #

Matlab的稀疏矩阵底层用**压缩列存储(CSC)**格式，和Eigen的SparseMatrix默认存储格式完全匹配，对接起来很顺畅。你只需要从Matlab的mxArray里提取三个核心数据：

• 非零元素的值数组：通过mxGetPr()获取
• 行索引数组：通过mxGetIr()获取
• 列指针数组：通过mxGetJc()获取

具体实现代码片段如下：

#include <Eigen/Sparse>
#include "mex.h"

using Eigen::SparseMatrix;

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // 先检查输入是否为稀疏矩阵
    if (!mxIsSparse(prhs[0])) {
        mexErrMsgTxt("Input must be a sparse matrix!");
    }

    int rows = mxGetM(prhs[0]);
    int cols = mxGetN(prhs[0]);
    int nnz = mxGetNzmax(prhs[0]); // 非零元素总数

    // 提取Matlab稀疏矩阵的核心数据
    double* values = mxGetPr(prhs[0]);
    mwIndex* row_indices = mxGetIr(prhs[0]);
    mwIndex* col_ptrs = mxGetJc(prhs[0]);

    // 构造Eigen的SparseMatrix
    SparseMatrix<double> eigen_sparse(rows, cols);
    eigen_sparse.reserve(nnz);
    for (int j = 0; j < cols; ++j) {
        for (mwIndex i = col_ptrs[j]; i < col_ptrs[j+1]; ++i) {
            eigen_sparse.insert(row_indices[i], j) = values[i];
        }
    }
    eigen_sparse.makeCompressed(); // 压缩成标准CSC格式
}

2. 从Eigen返回稀疏数组到Matlab #

把Eigen的稀疏矩阵传回Matlab也很简单，步骤如下：

1. 用mxCreateSparse()创建一个空的Matlab稀疏矩阵
2. 把Eigen的values()、innerIndices()、outerIndices()数据直接拷贝到Matlab的稀疏矩阵中
3. 因为Eigen和Matlab都是列优先存储，无需额外转置操作

代码示例：

// 假设已经构造好Eigen稀疏矩阵eigen_sparse
int rows = eigen_sparse.rows();
int cols = eigen_sparse.cols();
int nnz = eigen_sparse.nonZeros();

// 创建Matlab稀疏矩阵
plhs[0] = mxCreateSparse(rows, cols, nnz, mxREAL);

// 获取Matlab稀疏矩阵的指针
double* mx_values = mxGetPr(plhs[0]);
mwIndex* mx_ir = mxGetIr(plhs[0]);
mwIndex* mx_jc = mxGetJc(plhs[0]);

// 拷贝数据
std::copy(eigen_sparse.valuePtr(), eigen_sparse.valuePtr() + nnz, mx_values);
std::copy(eigen_sparse.innerIndexPtr(), eigen_sparse.innerIndexPtr() + nnz, mx_ir);
std::copy(eigen_sparse.outerIndexPtr(), eigen_sparse.outerIndexPtr() + cols + 1, mx_jc);

二、修改稠密数组相乘的MEX代码适配稀疏g #

你的原代码是处理稠密矩阵相乘，现在要适配稀疏的g，核心是把MatrixXd换成SparseMatrix<double>，同时修改输入数据的提取逻辑。下面是完整的修改后代码：

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include "mex.h"

using Eigen::MatrixXd;
using Eigen::SparseMatrix;
using namespace Eigen;

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // 检查输入数量是否正确
    if (nrhs != 2) {
        mexErrMsgTxt("需要输入两个参数：稀疏矩阵g和稠密矩阵G");
    }
    // 确保第一个输入是稀疏矩阵
    if (!mxIsSparse(prhs[0])) {
        mexErrMsgTxt("第一个参数g必须是稀疏矩阵！");
    }
    // 这里暂时假设G是稠密矩阵，如需支持稀疏G可自行扩展
    if (mxIsSparse(prhs[1])) {
        mexErrMsgTxt("第二个参数G暂时按稠密矩阵处理，如需稀疏可自行扩展");
    }

    // 提取稀疏矩阵g的信息
    int nRows_g = mxGetM(prhs[0]);
    int nCols_g = mxGetN(prhs[0]);
    int nnz_g = mxGetNzmax(prhs[0]);
    double* g_vals = mxGetPr(prhs[0]);
    mwIndex* g_ir = mxGetIr(prhs[0]);
    mwIndex* g_jc = mxGetJc(prhs[0]);

    // 构造Eigen稀疏矩阵g
    SparseMatrix<double> g_sparse(nRows_g, nCols_g);
    g_sparse.reserve(nnz_g);
    for (int j = 0; j < nCols_g; ++j) {
        for (mwIndex i = g_jc[j]; i < g_jc[j+1]; ++i) {
            g_sparse.insert(g_ir[i], j) = g_vals[i];
        }
    }
    g_sparse.makeCompressed();

    // 提取稠密矩阵G的信息
    int nRows_G = mxGetM(prhs[1]);
    int nCols_G = mxGetN(prhs[1]);
    double* Gr = mxGetPr(prhs[1]);
    Map<MatrixXd> G_map(Gr, nRows_G, nCols_G);

    // 稀疏矩阵乘稠密矩阵，结果默认是稠密矩阵（如需返回稀疏可修改类型）
    MatrixXd result = g_sparse * G_map;

    // 准备输出到Matlab
    plhs[0] = mxCreateDoubleMatrix(result.rows(), result.cols(), mxREAL);
    double* out_ptr = mxGetPr(plhs[0]);
    Map<MatrixXd> out_map(out_ptr, result.rows(), result.cols());
    out_map = result;
}

关键修改点说明： #

• 新增Eigen/Sparse头文件，引入SparseMatrix类型
• 增加输入类型校验，确保g是稀疏矩阵
• 替换g的处理逻辑：从提取稠密数组改为提取稀疏矩阵的CSC三元组数据，构造Eigen稀疏矩阵
• 稀疏矩阵与稠密矩阵相乘时，Eigen会自动处理运算；如果需要返回稀疏结果，只需把MatrixXd result改成SparseMatrix<double> result，再按前面的方法返回稀疏矩阵到Matlab即可

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