如何对存储为“压缩稀疏行”的矩阵进行稀疏矩阵索引？

Posted 2023-03-24

【中文标题】如何对存储为“压缩稀疏行”的矩阵进行稀疏矩阵索引？

【英文标题】：How sparse matrix-indexing can be done for a matrix stored as "Compressed Sparse Row"?

【发布时间】：2019-03-12 01:23:22

【问题描述】：

我使用Intel MKL 将我的大型稀疏对称矩阵存储为压缩稀疏行 (CSR)。例如，假设我的对称稀疏矩阵是5x5：

A =
    1    -1     0    -3     0
   -1     5     0     0     0
    0     0     4     6     4
   -3     0     6     7     0
    0     0     4     0    -5

values   = 1,  -1,  -3,   5,   4,   6,   4,   7,  -5; // symmetric sparse matrix
columns  = 0,   1,   3,   1,   2,   3,   4,   3,   4; // zero-based
rowIndex = 0,   3,   4,   7,   8,   9; // zero-based

我正在尝试在给定行和列的情况下找到A 的子矩阵，例如A(1:3, 2:4)：

A(1:3,2:4) =
   0     0     0
   4     6     4
   6     7     0

values   = 4,   6,   4,   6,   7; // General sparse matrix (sub-matrix is not necessarily symmetric)
columns  = 0,   1,   2,   0,   1; // zero-based
rowIndex = 0,   0,   3,   5; // zero-based

如果知道如何进行矩阵索引，我将不胜感激。我能想到的一种方法是将CSR 转换为坐标格式COO 并应用矩阵索引，然后将其转换回CSR，我认为这不是一种有效的方法。

有人可以告诉我一种有效或常用的稀疏矩阵索引方法吗？

【问题讨论】：

您似乎知道这一点，但重要的是要注意结果通常不是对称的（当然，当它保持对称时很容易注意到）。

如果你有一个完整的矩阵存储为 CSR，它应该很容易。在这里，您必须重建（显式或隐式）缺失的部分。我的第一个猜测是它应该类似于 CSR CSC 格式转换，我建议你看看那个算法（github.com/scipy/scipy/blob/…）。

【参考方案1】：

诀窍是通过输出列（即它们的行）在下三角形中查找值。您可以为每一行的数据保留一个索引，因为您在输出的行顺序中按列顺序访问条目。

带说明型

struct CSR   // sometimes implicitly symmetric
  std::vector<...> vals;
  std::vector<int> cols,rowStart;
;

我们有

// Return the [r0,r1) by [c0,c1) submatrix, never
// using any symmetry it might have.
CSR submatrix(const CSR &sym,int r0,int r1,int c0,int c1) 
  const int m=r1-r0,n=c1-c0;
  std::vector<int> finger(sym.rowStart.begin()+c0,sym.rowStart.begin()+c1);
  CSR ret;
  ret.rowStart.reserve(m+1);
  ret.rowStart.push_back(0);
  for(int r=0,rs=r0;r<m;++r,++rs) 
    // (Strictly) lower triangle:
    for(int cs=c0,c=0;cs<rs;++cs,++c)
      for(int &f=finger[c],f1=sym.rowStart[cs+1];f<f1;++f) 
        const int cf=sym.cols[f];
        if(cf>rs) break;
        if(cf==rs) 
          ret.vals.push_back(sym.vals[f]);
          ret.cols.push_back(c);
        
      
    // Copy the relevant subsequence of the upper triangle:
    for(int f=sym.rowStart[rs],f1=sym.rowStart[rs+1];f<f1;++f) 
      const int c=sym.cols[f]-c0;
      if(c<0) continue;
      if(c>=n) break;
      ret.vals.push_back(sym.vals[f]);
      ret.cols.push_back(c);
    
    ret.rowStart.push_back(ret.vals.size());
  
  return ret;

对于大型矩阵，上三角循环可以通过使用二分搜索来优化f的相关范围。

【讨论】：

确保使用正确的版本；第一个有错别字但仍然编译！