Eigen：比较两个可能具有不同稀疏模式的稀疏矩阵答案

【问题标题】：Eigen: Compare two sparse matrix with probably different sparsity patternEigen：比较两个可能具有不同稀疏模式的稀疏矩阵
【发布时间】：2019-06-11 14:02:39
【问题描述】：

我想比较两个Eigen::SparseMatrix

存在 res.isApprox(ans) 方法，但不幸的是，据我所知，在不同的稀疏模式的情况下，它会因断言而失败

AffineInvariantDeformerTest01: Eigen/src/SparseCore/SparseMatrix.h:934: void Eigen::internal::set_from_triplets(const InputIterator&, const InputIterator&, SparseMatrixType&, DupFunctor) [with InputIterator = __gnu_cxx::__normal_iterator<Eigen::Triplet<double, int>*, std::vector<Eigen::Triplet<double, int>, std::allocator<Eigen::Triplet<double, int> > > >; SparseMatrixType = Eigen::SparseMatrix<double, 0, int>; DupFunctor = Eigen::internal::scalar_sum_op<double, double>]: Assertion `it->row()>=0 && it->row()<mat.rows() && it->col()>=0 && it->col()<mat.cols()' failed.

我希望它用于单元测试，所以如果它不会那么快也没关系。可能我可以将稀疏矩阵转换为密集矩阵，但我希望存在更优雅的解决方案

编辑：只是一种可以比较两个矩阵的稀疏模式的方法对我来说也可以

【问题讨论】：

如果稀疏模式不同，它应该返回false，还是应该仍然是true，如果结构上不同的条目恰好/接近于零？
@chtz 感谢您的评论！对于我的问题，不同的模式意味着不同的矩阵。

标签： c++ sparse-matrix eigen eigen3

【解决方案1】：

首先，isApprox() 工作正常，即使使用不同的稀疏模式。看起来您的错误在其他地方（可能已经在您使用 setFromTriplets() 设置矩阵时）

当给定两个具有不同稀疏模式的矩阵时，如果不同的条目（几乎）为零，则它们将被视为近似相等。以下应计算两次true（增加1e-19 以查看差异）：

#include <Eigen/SparseCore>
#include <iostream>
#include <array>

int main() {
    Eigen::SparseMatrix<double> Mat1(2,2), Mat2(2,2);

    std::array<Eigen::Triplet<double,int>, 2> data1  {{{0,0,1.0}, {1,1, 1e-19}}};
    std::array<Eigen::Triplet<double,int>, 2> data2  {{{0,0,1.0}, {1,0, 1e-19}}};

    Mat1.setFromTriplets(data1.begin(), data1.end());
    Mat2.setFromTriplets(data2.begin(), data2.end());

    std::cout << "Mat1.isApprox(Mat1) == " << Mat1.isApprox(Mat1) << "\nMat1.isApprox(Mat2) == " << Mat1.isApprox(Mat2) << "\n";
}

如果要比较稀疏模式，可以检查从内部和外部索引指针开始的数据（两者都使用 Map）。如果两个矩阵具有相同的类型并且被压缩，则以下工作：

template<class Derived>
bool hasSamePattern(Eigen::SparseCompressedBase<Derived> const& A, Eigen::SparseCompressedBase<Derived> const& B)
{
    assert(A.isCompressed() && B.isCompressed());
    if(A.rows() != B.rows() || A.cols() != B.cols() || A.nonZeros() != B.nonZeros())
        return false;
    typedef Eigen::Matrix<typename Derived::StorageIndex, Eigen::Dynamic, 1> IndexVector;
    Eigen::Index outerSize = A.outerSize(), nnz = A.nonZeros();
    if(IndexVector::Map(A.outerIndexPtr(), outerSize) != IndexVector::Map(B.outerIndexPtr(), outerSize))
        return false;
    if(IndexVector::Map(A.innerIndexPtr(), nnz) != IndexVector::Map(B.innerIndexPtr(), nnz))
        return false;

    return true;
}

【讨论】：