eigen稀疏矩阵拼接（基于块操作的二维拼接）的思考

转载请说明出处：eigen稀疏矩阵拼接（块操作）

eigen稀疏矩阵拼接（块操作）

关于稀疏矩阵的块操作：参考官方链接

However, for performance reasons, writing to a sub-sparse-matrix is much more limited, and currently only contiguous sets of columns (resp. rows) of a column-major (resp. row-major) SparseMatrix are writable. Moreover, this information has to be known at compile-time, leaving out methods such as block(…) and corner*(…). The available API for write-access to a SparseMatrix are summarized below:

上面规则对应的示例：

SparseMatrix<double,ColMajor> sm1;
sm1.col(j) = ...;
sm1.leftCols(ncols) = ...;
sm1.middleCols(j,ncols) = ...;
sm1.rightCols(ncols) = ...;

SparseMatrix<double,RowMajor> sm2;
sm2.row(i) = ...;
sm2.topRows(nrows) = ...;
sm2.middleRows(i,nrows) = ...;
sm2.bottomRows(nrows) = ...;

总结为一句话：列优先的稀疏矩阵只能进行列向操作，行优先的稀疏矩阵只能进行横向操作

• 这个限制说明很多Dense矩阵的块操作对Sparse矩阵来说就不成立了：

  • 对正常的Dense矩阵来说，拼接比较方便，可以直接使用：

    M << m1, m2;
    
    M << m1,
         m2;

  • Dense矩阵可以对整个稀疏矩阵的某个block进行修改，比如只是对上述D部分进行修改，稀疏矩阵是不行的

  • Sparse矩阵使用的块操作更多的用读取权限，而没有写权限

• 在这个限制的前提下，当我们想进行两维稀疏矩阵的拼接的时候就会很矛盾，举个例子，我们现在有四个稀疏矩阵A B C D，如下：

  $M= \left[ \begin{matrix} A & B \\ C & D \\ \end{matrix} \right]$

  参考上面的规则代码，如果所有的稀疏矩阵都使用列优先，那么可以将A B进行列向拼接，即$\left[\begin{matrix}A & B \\\end{matrix}\right]$，C D也一样$\left[\begin{matrix}C & D \\\end{matrix}\right]$，但是之后我们就无法将结果的$\left[\begin{matrix}A & B \\\end{matrix}\right]$和$\left[\begin{matrix}C & D \\\end{matrix}\right]$的结果进行横向拼接为$\left[ \begin{matrix} A & B \\ C & D \\ \end{matrix} \right]$了。

• 解决方案

  • 笨方法，用insert进行逐个插值

    此时借助迭代器是一个不错的选择，给出一段示例程序：

    Eigen::SparseMatrix<double> M_sparse(2*npixels, 2*npixels);
    M_sparse.reserve(FU_sparse.nonZeros() + FV_sparse.nonZeros());
    for(int c=0; c<FU_sparse.cols(); ++c)
    {
        for(Eigen::SparseMatrix<double>::InnerIterator itU(FU_sparse, c); itU; ++itU)
        {
            if (itU.value())
            {
                M_sparse.insert(itU.row(), c) = -itU.value();
            }
        }
    }

    问题是，当非零元素的数目过w的时候，上述hash-map的insert弊端就会凸显出来，也是很耗时的。

  • 在上述限制的前提下利用现有接口实现上述过程

    $A= \left[ \begin{matrix} 1 & 2 \\ 3 & 4 \\ \end{matrix} \right]$

    $B= \left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]$

    $C= \left[ \begin{matrix} 11 & 12 \\ 13 & 14 \\ \end{matrix} \right]$

    $D= \left[ \begin{matrix} A & B \\ C & D \\ \end{matrix} \right]$

    我们想要的结果是：
    $M= \left[ \begin{matrix} 1 & 2 & a & b \\ 3 & 4 & c & d\\ 11 & 12 & A & B \\ 13 & 14 & C & D\\ \end{matrix} \right]$

    1. 对A、B、C、D进行转置得到
       $A= \left[ \begin{matrix} 1 & 3 \\ 2 & 4 \\ \end{matrix} \right]$

       $B= \left[ \begin{matrix} a & c \\ b & d \\ \end{matrix} \right]$

       $C= \left[ \begin{matrix} 11 & 13 \\ 12 & 14 \\ \end{matrix} \right]$

       $D= \left[ \begin{matrix} A & C \\ B & D \\ \end{matrix} \right]$

    2. 将A和C、B和D进行列向拼接，得到：
       $m1= \left[ \begin{matrix} 1 & 3 & 11 & 13 \\ 2 & 4 & 12 & 14 \\ \end{matrix} \right]$

       $m2= \left[ \begin{matrix} a & c & A & C \\ b & d & B & D \\ \end{matrix} \right]$

    3. 将m1和m2进行转置，得到：
       $m1= \left[ \begin{matrix} 1 & 2 \\ 3 & 4 \\ 11 & 12 \\ 13 & 14 \\ \end{matrix} \right]$

       $m2= \left[ \begin{matrix} a & b \\ c & d \\ A & b \\ C & d \\ \end{matrix} \right]$

    4. 将m1和m2进行列项拼接，得到：

    $M= \left[ \begin{matrix} 1 & 2 & a & b \\ 3 & 4 & c & d\\ 11 & 12 & A & B \\ 13 & 14 & C & D\\ \end{matrix} \right]$

    这就是我们想要的结果。

    给出一段简单的示例程序（B和D为全零矩阵）：

    Eigen::SparseMatrix<double> m1(3,3);
    m1.insert(0,0)=1;
    m1.insert(0,1)=2;
    m1.insert(0,2)=3;
    m1.insert(1,0)=4;
    m1.insert(1,1)=5;
    m1.insert(1,2)=6;
    m1.insert(2,0)=7;
    m1.insert(2,1)=8;
    m1.insert(2,2)=9;
    Eigen::SparseMatrix<double> m1_tp  = m1.transpose();
    
    Eigen::SparseMatrix<double> m2(3,3);
    m2.insert(0,0)=11;
    m2.insert(0,1)=12;
    m2.insert(0,2)=13;
    m2.insert(1,0)=14;
    m2.insert(1,1)=15;
    m2.insert(1,2)=16;
    m2.insert(2,0)=17;
    m2.insert(2,1)=18;
    m2.insert(2,2)=19;
    Eigen::SparseMatrix<double> m2_tp  = m2.transpose();
    
    Eigen::SparseMatrix<double> zero(3,3);
    
    Eigen::SparseMatrix<double> m_top_re(3,6);
    Eigen::SparseMatrix<double> m_bottom_re(3,6);
    
    m_top_re.leftCols(3) = m1_tp;
    m_top_re.rightCols(3) = zero;
    
    m_bottom_re.leftCols(3) = zero;
    m_bottom_re.rightCols(3) = m2_tp;
    
    Eigen::SparseMatrix<double> m_leftcol_re  = m_top_re.transpose();
    Eigen::SparseMatrix<double> m_rightcol_re  = m_bottom_re.transpose();
    
    Eigen::SparseMatrix<double> m3(6,6);
    m3.leftCols(3)  = m_leftcol_re;
    m3.rightCols(3) = m_rightcol_re;
    
    cout <<m3<<endl;

  ---

  • 接下来思考一个问题

    • 为什么列优先的稀疏矩阵不支持横向拼接？

      即$\left[ \begin{matrix} A \\ B \\ \end{matrix} \right]$

      理论上hash表是有顺序的，比如某一个key后面对应多个value，它们是有先后顺序的，这个顺序也代表了他们在矩阵中的先后位置关系。因为如果是列优先的矩阵，那么hash表在进行存储的时候是列向连续读取并存入的，如果我们在列向方向继续加入新的矩阵就等同于hash表继续读取并存入。但加入我们在横向方向加入新的矩阵（上面的A B），就相当于每一列都加入了新的元素，得到的新矩阵也应该符合列有限的规则，那么就相当于我们要重新读取每一列的数据，进行hash，这个过程会涉及大量的数据拷贝操作，会非常耗时。