时间戳

首页 > 解决方案 > 你将如何合并稀疏矩阵来创建一个新的稀疏矩阵？

algorithm - 你将如何合并稀疏矩阵来创建一个新的稀疏矩阵？

问题描述

我正在努力实现以下目标：

Given A(m,n), B(m,q), C(p,n), D(p,q), sparse matrices.
Create E(m+p,n+q), a sparse matrix, like :
E = | A B |
    | C D |

我尝试了以下方法：

这些库的问题在于 E 的构造就像任何稀疏矩阵一样，而没有使用 E 和 A、B、C、D 的部分之间存在冗余的事实。我想应该可以通过重新索引 A、B、C 和 D 的内部结构中存储的内容来构造 E。

问题如下：您将如何使用稀疏矩阵实现这种合并操作？你会用什么软件？你会使用哪种算法？

理想的解决方案不会使用块存储方案，因此稀疏性将基于零值，而不是零块。

编程语言无关紧要。

提前致谢。

标签： algorithmmatrixsparse-matrixeigen3intel-mkl

解决方案

老问题，新解决方案！似乎每个人都找到了一种方法来通过创建一个新的三元组容器并相应地添加来做到这一点。这无疑是安全和直观的，但即使您相应地保留非零值，也不会那么快。

在这种仅适用于列优先矩阵的方法中（您可以很容易地制作行优先版本），我直接使用 CSC 数据结构。我提出了 3 种方法：1）“垂直堆叠”两个具有相同列数的稀疏矩阵，2）“水平堆叠”两个具有相同行数的稀疏矩阵，以及 3）“就地水平堆叠”两个稀疏矩阵具有相同的行数。3）本质上更有效，因为左矩阵的存储不变。

如您所见，即使直接使用原始指针，水平堆叠两个矩阵也几乎是微不足道的。垂直堆叠需要一些思考，但并不多。如果有人能找出一种避免一些复制（或任何改进）的垂直堆叠的就地方法，请告诉我！

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

template<typename Scalar, typename StorageIndex>
void sparse_stack_v(
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& top,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& bottom,
    SparseMatrix<Scalar, ColMajor, StorageIndex>& stacked)
{
    assert(top.cols() == bottom.cols());
    stacked.resize(top.rows() + bottom.rows(), top.cols());
    stacked.resizeNonZeros(top.nonZeros() + bottom.nonZeros());

    StorageIndex i = 0;

    for (StorageIndex col = 0; col < top.cols(); col++)
    {
        stacked.outerIndexPtr()[col] = i;

        for (StorageIndex j = top.outerIndexPtr()[col]; j < top.outerIndexPtr()[col + 1]; j++, i++)
        {
            stacked.innerIndexPtr()[i] = top.innerIndexPtr()[j];
            stacked.valuePtr()[i] = top.valuePtr()[j];
        }

        for (StorageIndex j = bottom.outerIndexPtr()[col]; j < bottom.outerIndexPtr()[col + 1]; j++, i++)
        {
            stacked.innerIndexPtr()[i] = (StorageIndex)top.rows() + bottom.innerIndexPtr()[j];
            stacked.valuePtr()[i] = bottom.valuePtr()[j];
        }
    }
    stacked.outerIndexPtr()[top.cols()] = i;
}

template<typename Scalar, typename StorageIndex>
void sparse_stack_h(
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& right,
    SparseMatrix<Scalar, ColMajor, StorageIndex>& stacked)
{
    assert(left.rows() == right.rows());

    stacked.resize(left.rows(), left.cols() + right.cols());
    stacked.resizeNonZeros(left.nonZeros() + right.nonZeros());

    std::copy(left.innerIndexPtr(), left.innerIndexPtr() + left.nonZeros(), stacked.innerIndexPtr());
    std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), stacked.innerIndexPtr() + left.nonZeros());

    std::copy(left.valuePtr(), left.valuePtr() + left.nonZeros(), stacked.valuePtr());
    std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), stacked.valuePtr() + left.nonZeros());

    std::copy(left.outerIndexPtr(), left.outerIndexPtr() + left.cols(), stacked.outerIndexPtr());//dont need the last entry of A.outerIndexPtr() -- total length is AB.cols() + 1 = A.cols() + B.cols() + 1
    std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, stacked.outerIndexPtr() + left.cols(), [&](StorageIndex i) { return i + left.nonZeros(); });
}

template<typename Scalar, typename StorageIndex>
void sparse_stack_h_inplace(
    SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& right)
{
    assert(left.rows() == right.rows());

    const StorageIndex leftcol = (StorageIndex)left.cols();
    const StorageIndex leftnz = (StorageIndex)left.nonZeros();

    left.conservativeResize(left.rows(), left.cols() + right.cols());
    left.resizeNonZeros(left.nonZeros() + right.nonZeros());

    std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), left.innerIndexPtr() + leftnz);
    std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), left.valuePtr() + leftnz);
    std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, left.outerIndexPtr() + leftcol, [&](StorageIndex i) { return i + leftnz; });
}

推荐阅读