如何用CWiseNullaryOp构造CSR稀疏单位矩阵？

核心思路 #

你不需要编写三个仿函数，也不用把整个矩阵写成仿函数。仅需一个简单的仿函数判断对角元素，结合Eigen的NullaryExpr（它会生成CWiseNullaryOp表达式），就能快速得到CSR格式的单位矩阵。另外，如果你预先知晓CSR的三个数组结构，直接构造数组初始化稀疏矩阵会更高效。

方法一：用CWiseNullaryOp（NullaryExpr）生成 #

这种方法通过仿函数定义元素值的规则，自动生成稀疏矩阵的CSR结构：

#include <Eigen/Sparse>
#include <iostream>

// 定义仿函数：仅对角元素返回1，其余返回0
struct IdentityFunctor {
    template<typename Index>
    double operator()(Index row, Index col) const {
        return (row == col) ? 1.0 : 0.0;
    }
};

int main() {
    const int matrix_size = 5; // 5阶单位矩阵

    // 创建NullaryExpr表达式，转成RowMajor的SparseMatrix（即CSR格式）
    Eigen::SparseMatrix<double, Eigen::RowMajor> identity_matrix = 
        Eigen::SparseMatrix<double, Eigen::RowMajor>(
            Eigen::NullaryExpr(matrix_size, matrix_size, IdentityFunctor())
        );

    // 验证矩阵内容
    std::cout << "CSR单位矩阵元素：\n";
    for (int outer_idx = 0; outer_idx < identity_matrix.outerSize(); ++outer_idx) {
        for (Eigen::SparseMatrix<double>::InnerIterator it(identity_matrix, outer_idx); it; ++it) {
            std::cout << "(" << it.row() << "," << it.col() << ") = " << it.value() << "\n";
        }
    }

    // 查看CSR的三个核心数组
    std::cout << "\nrow_ptr数组：";
    for (int i = 0; i <= matrix_size; ++i) {
        std::cout << identity_matrix.outerIndexPtr()[i] << " ";
    }
    std::cout << "\ncol_indices数组：";
    for (int i = 0; i < identity_matrix.nonZeros(); ++i) {
        std::cout << identity_matrix.innerIndexPtr()[i] << " ";
    }
    std::cout << "\nvalues数组：";
    for (int i = 0; i < identity_matrix.nonZeros(); ++i) {
        std::cout << identity_matrix.valuePtr()[i] << " ";
    }

    return 0;
}

代码说明

• IdentityFunctor：仿函数接收行、列索引，判断是否为对角元素，返回对应值。
• NullaryExpr：生成一个CWiseNullaryOp表达式，描述矩阵每个元素的计算规则。
• Eigen::RowMajor：指定稀疏矩阵按行优先存储，对应CSR格式（outerIndexPtr是row_ptr，innerIndexPtr是col_indices，valuePtr是values）。

方法二：直接构造CSR数组初始化 #

如果你已经明确CSR三个数组的结构，直接构造数组初始化会更高效（避免遍历所有矩阵元素）：

#include <Eigen/Sparse>
#include <iostream>
#include <vector>
#include <algorithm>

int main() {
    const int matrix_size = 5;
    int nnz = matrix_size; // 单位矩阵的非零元素数等于阶数

    // 构造CSR的三个数组
    std::vector<int> row_ptr(matrix_size + 1);
    std::vector<int> col_indices(matrix_size);
    std::vector<double> values(matrix_size, 1.0);

    // 生成row_ptr：[0,1,2,...,matrix_size]
    for (int i = 0; i <= matrix_size; ++i) {
        row_ptr[i] = i;
    }
    // 生成col_indices：[0,1,...,matrix_size-1]
    for (int i = 0; i < matrix_size; ++i) {
        col_indices[i] = i;
    }

    // 初始化CSR稀疏矩阵
    Eigen::SparseMatrix<double, Eigen::RowMajor> identity_matrix(matrix_size, matrix_size);
    identity_matrix.resizeNonZeros(nnz);
    std::copy(row_ptr.begin(), row_ptr.end(), identity_matrix.outerIndexPtr());
    std::copy(col_indices.begin(), col_indices.end(), identity_matrix.innerIndexPtr());
    std::copy(values.begin(), values.end(), identity_matrix.valuePtr());

    // 验证矩阵内容
    std::cout << "直接构造的CSR单位矩阵元素：\n";
    for (int outer_idx = 0; outer_idx < identity_matrix.outerSize(); ++outer_idx) {
        for (Eigen::SparseMatrix<double>::InnerIterator it(identity_matrix, outer_idx); it; ++it) {
            std::cout << "(" << it.row() << "," << it.col() << ") = " << it.value() << "\n";
        }
    }

    return 0;
}

代码说明

• 直接生成row_ptr、col_indices、values三个数组，完全符合单位矩阵的CSR结构。
• 通过resizeNonZeros指定非零元素数，再用std::copy将数组内容复制到稀疏矩阵的内部指针中。

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