如何基于Eigen实现类CWiseNullaryOp的稀疏矩阵函子式映射？

为什么原代码无法运行 #

Eigen的Map<SparseMatrix>要求传入实际内存中的数组指针作为行指针、列索引和值的数据源，它不支持直接通过函子生成数据，因此直接传递函子会触发编译错误。

方案一：手动生成数组后构建稀疏矩阵 #

先通过函子生成行指针、列索引和值的数组，再用这些数组初始化稀疏矩阵，这是最直接的实现方式：

#include <iostream>
#include <Eigen/Sparse>
#include <vector>
using namespace std;
using namespace Eigen;

class index_class {
public:
    int operator() (int index) {
        return index;
    }
};

class one_class {
public:
    double operator() (int index) { // 改为double类型匹配矩阵元素类型
        return 1.0;
    }
};

int main(){
    int size = 4;
    int nnz = 4;

    one_class nnz_functor;
    index_class col_functor;
    index_class row_ptr_functor;

    // 生成行指针数组（RowMajor下长度为size+1）
    vector<int> row_ptr(size + 1);
    for(int i = 0; i <= size; ++i){
        row_ptr[i] = row_ptr_functor(i);
    }

    // 生成列索引数组
    vector<int> col_indices(nnz);
    for(int i = 0; i < nnz; ++i){
        col_indices[i] = col_functor(i);
    }

    // 生成值数组
    vector<double> values(nnz);
    for(int i = 0; i < nnz; ++i){
        values[i] = nnz_functor(i);
    }

    // 初始化稀疏矩阵
    SparseMatrix<double, RowMajor, int> mat(size, size);
    mat.makeCompressed();
    // 直接赋值内部数组指针
    mat.outerIndexPtr() = row_ptr.data();
    mat.innerIndexPtr() = col_indices.data();
    mat.valuePtr() = values.data();
    mat.nonZeros() = nnz;

    // 验证输出
    cout << "稀疏矩阵：" << endl << mat << endl;

    return 0;
}

方案二：使用Triplet列表构建（更灵活） #

如果函子生成行/列的逻辑更复杂，推荐用Triplet列表来构建稀疏矩阵，无需手动处理行指针的细节：

#include <iostream>
#include <Eigen/Sparse>
#include <vector>
using namespace std;
using namespace Eigen;

class index_class {
public:
    int operator() (int index) {
        return index;
    }
};

class one_class {
public:
    double operator() (int index) {
        return 1.0;
    }
};

int main(){
    int size = 4;
    int nnz = 4;

    one_class nnz_functor;
    index_class col_functor;
    index_class row_functor; // 这里直接生成行索引

    vector<Triplet<double>> triplets;
    triplets.reserve(nnz);

    // 用函子生成每个非零元素的(row, col, value)三元组
    for(int i = 0; i < nnz; ++i){
        int row = row_functor(i);
        int col = col_functor(i);
        double val = nnz_functor(i);
        triplets.emplace_back(row, col, val);
    }

    // 构建稀疏矩阵
    SparseMatrix<double, RowMajor, int> mat(size, size);
    mat.setFromTriplets(triplets.begin(), triplets.end());

    // 验证输出
    cout << "稀疏矩阵：" << endl << mat << endl;

    return 0;
}

关键注意事项 #

• 确保函子的返回值类型与稀疏矩阵的元素类型匹配（比如矩阵是double，值函子要返回double）
• RowMajor格式的稀疏矩阵，行指针数组长度为行数+1，最后一个元素等于非零元素总数
• 如果需要频繁复用这种逻辑，可以把数组生成和矩阵封装成一个工具函数

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