如何将稀疏矩阵和密集向量相乘？答案

【问题标题】：How to multiply a sparse matrix and a dense vector?如何将稀疏矩阵和密集向量相乘？
【发布时间】：2021-12-18 05:52:57
【问题描述】：

我正在尝试以下方法：

Eigen::SparseMatrix<double> bijection(2 * face_count, 2 * vert_count);
/* initialization */
Eigen::VectorXd toggles(2 * vert_count);
toggles.setOnes();
Eigen::SparseMatrix<double> deformed;
deformed = bijection * toggles;

Eigen 正在返回一个错误声明：

 error: static assertion failed: THE_EVAL_EVALTO_FUNCTION_SHOULD_NEVER_BE_CALLED_FOR_DENSE_OBJECTS
  586 |       EIGEN_STATIC_ASSERT((internal::is_same<Dest,void>::value),THE_EVAL_EVALTO_FUNCTION_SHOULD_NEVER_BE_CALLED_FOR_DENSE_OBJECTS);

根据eigen documentaion

允许使用稀疏矩阵和向量积。我做错了什么？

【问题讨论】：

标签： c++ math matrix eigen

【解决方案1】：

问题是您的产品输出类型错误。

Eigen 文档指出定义了以下乘法类型：

dv2 = sm1 * dv1;

稀疏矩阵乘以密集向量等于密集向量。

如果您确实需要稀疏表示，我认为没有比执行上述乘法然后使用sparseView 成员函数将乘积转换为稀疏矩阵更好的方法了。例如

Eigen::SparseMatrix<double> bijection(2 * face_count, 2 * vert_count);
/* initialization */
Eigen::VectorXd toggles(2 * vert_count);
toggles.setOnes();
Eigen::VectorXd deformedDense = bijection * toggles;
Eigen::SparseMatrix<double> deformedSparse = deformedDense.sparseView();

【讨论】：

【解决方案2】：

如果它非常稀疏，这可能比输出到密集向量更快。否则，比传统产品快 99/100 倍。

void sparsem_densev_sparsev(const SparseMatrix<double>& A, const VectorX<double>& x, SparseVector<double>& Ax)
{
    Ax.resize(x.size());

    for (int j = 0; j < A.outerSize(); ++j)
    {
        if (A.outerIndexPtr()[j + 1] - A.outerIndexPtr()[j] > 0)
        {
            Ax.insertBack(j) = 0;
        }
    }

    for (int j_idx = 0; j_idx < Ax.nonZeros(); j_idx++)
    {
        int j = Ax.innerIndexPtr()[j_idx];

        for (int k = A.outerIndexPtr()[j]; k < A.outerIndexPtr()[j + 1]; ++k)
        {
            int i = A.innerIndexPtr()[k];
            Ax.valuePtr()[j_idx] += A.valuePtr()[k] * x.coeff(i);
        }
    }
}

对于（可能不是最佳的）自伴随版本（下三角形），将j_idx 循环更改为：


for (int j_idx = 0; j_idx < Ax.nonZeros(); j_idx++)
    {
        int j = Ax.innerIndexPtr()[j_idx];
        int i_idx = j_idx;//i>= j, trick to improve binary search

        for (int k = A.outerIndexPtr()[j]; k < A.outerIndexPtr()[j + 1]; ++k)
        {
            int i = A.innerIndexPtr()[k];
            Ax.valuePtr()[j_idx] += A.valuePtr()[k] * x.coeff(i);
            if (i != j)
            {
                i_idx = std::distance(Ax.innerIndexPtr(), std::lower_bound(Ax.innerIndexPtr() + i_idx, Ax.innerIndexPtr() + Ax.nonZeros(), i));
                Ax.valuePtr()[i_idx] += A.valuePtr()[k] * x.coeff(j);
            }
        }
    }

【讨论】：