最近在开发的库涉及到对CSR矩阵删除指定行列的功能，在网上检索了很多都没有很好的答案，现成的稀疏矩阵库（python中的scipy、C++中的Eigen等）也都没有对应的方法，在这里分享一下我的解决方案。
何谓CSR矩阵
稀疏矩阵（Sparse Matrix）是指矩阵中含有大量元素的矩阵，在一般的计算密集型任务中，我们需要计算的矩阵维度高达数百万甚至上亿级，为了节省存储空间和提高计算效率，人们设计了种种稀疏矩阵的存储方式来进行存储，即不必将矩阵中的每一个元素都存储下来，我们只需找到那些非零元素（Non-zero, NNZ）并存储它们在矩阵中的位置和值即可。
CSR矩阵的全称是压缩行稀疏矩阵（Compressed Row Sparse Matrix），它是一种稀疏矩阵的存储格式，使用该数据结构存储矩阵能大大加快矩阵组装和计算的效率。除此之外，稀疏矩阵的存储格式还有COO、CSC、BSR等，COO格式（Coordinate Format）是一种更容易被理解的格式，因此我们以此来与CSR做对比来帮助理解。在高性能有限元求解的问题上，我们常用的是CSR矩阵。
假设有如下图所示的矩阵，这是一个的对称矩阵，尽管它的稀疏度不是那么高，但是作为一个案例来说明问题足以。
对于COO格式，它的存储是一个三元组列表，每一个非零元素都由一个三元组来存储，三元组中的三个元素分别是行索引、列索引和值，这样就可以知道每个非零元素的确切位置以及对应的值，在矩阵计算中可以快速访问这些元素。使用COO格式存储的好处有：
（1）它非常方便构建矩阵，我们只需给定行列位置和值的大小即可确定一个值list移除指定元素，这也是我们最常用的矩阵赋值方法。
（2）它不存在三元组之间的顺序的问题，我们可以任意插入或删除某个非零值。
（3）COO构建好的矩阵也很容易转换成其他类型矩阵。
对于矩阵，使用COO格式的存储结果为：
(0, 0, 10)
(0, 1, 20)
(1, 0, 20)
(1, 1, 30)
......
(5, 5, 100)
为了方便后面的对比，我们把行、列、值单独拿出来，得到下面的三个数组：
Rows = [0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5, 5]
Cols = [0, 1, 0, 1, 2, 3, 1, 2, 5, 1, 3, 4, 5, 2, 4, 5]
Vals = [10, 20, 20, 30, 35, 40, 35, 50, 60, 40, 70, 80, 90, 60, 90, 100]
对于CSR格式list移除指定元素，它是一种行优先的压缩存储格式，它不会再像COO格式一样把每一行非零值的行列坐标都记录下来，而是从0开始，只记录每一行的非零元素的数量，然后由列索引来找出每一个元素具体的位置。例如矩阵的第一行只有两个非零值、，所以Rows[1] = 2，第二行有4个非零值，因此Rows[2] = 2 + 4 = 6。
Rows = [0, 2, 6, 9, 11, 13, 16]
Cols = [0, 1, 0, 1, 2, 3, 1, 2, 5, 1, 3, 4, 5, 2, 4, 5]
Vals = [10, 20, 20, 30, 35, 40, 35, 50, 60, 40, 70, 80, 90, 60, 90, 100]
可以看出，对于一个矩阵，行数组的元素数量是N+1个，其中最后一个元素的值就是该稀疏矩阵的非零元素数量，即Rows[N] = nnz。与COO矩阵不同的是，为了确保每个元素的位置，每一行中的列索引和值要按照元素在矩阵中的位置顺序排列，即从左至右、从上至下。
因此，如果我们要构造一个CSR Matrix类，其构造的参考代码为：
public class CSRMatrix
{
///
/// Rows contains indices of the first non-zero elements in each row. The size of the rows array is N+1, and Rows[N] = NNZ.
///
public int[] Rows; // structure of the sparse matrix
///
/// The Cols array holds column indices of the elements, hence its size is also equal to NNZ.
///
public int[] Cols;
///
/// The Vals array holds only the non-zero elements listed left-to-right, top-to-bottom, and its size is NNZ.
///
public double[] Vals;
///
/// The size of CSR matrix.
///
public int N;
///
/// The number of non-zero entries
///
public int NNZ;
public CSRMatrix(int n, int nnz)
{
N = n;
NNZ = nnz;
Rows = new int[n + 1];
Rows[n] = nnz;
Cols = new int[NNZ];
Cols[n] = nnz;
Vals = new double[NNZ];
}
}
删除指定的行列
在组装矩阵准备求解大型稀疏方程组时，以有限元为例，我们需要在总体刚度矩阵中移除被约束的自由度（即指定的行列），这在matlab中可以直接通过局部指定行列来执行，对于python，可以通过numpy中的delete来执行，然而，python中并没有适用于稀疏矩阵的方法。对于C++而言，常用的线性代数库Eigen也没有现成的方法供我们使用，我曾在网上搜索过答案，Eigen开发者的意思是要通过置换矩阵把被约束的行列置换到最后，然后再取前面的行列构造的子矩阵来进行计算，我尝试了这种方法，但是结果看上去并不是我想要的（也可能是我没有仔细较对代码），因此该方法暂存疑，有了解的小伙伴欢迎私信我交流分享。由于笔者的代码主体部分是使用C#语言编写，仅求解过程使用C++语言，因此下面实现部分采用C#语言，完整代码请参考我的有限元库中的CSRMatrix.cs部分
下面来分享下我的方法，其实是个笨办法，只是直接计算出了删除指定行列后的三个数组，由于过程有点麻烦，自己也不想在下次遇到的时候再回过头来想，因此记录在此以便以后需要的时候再次查阅。
我们还是以上面的矩阵为例，矩阵及原本的CSR格式矩阵表示如下，假设我们要删除的是第2行和第2列（注意行列编号从0开始）。
Rows = [0, 2, 6, 9, 11, 13, 16]
Cols = [0, 1, 0, 1, 2, 3, 1, 2, 5, 1, 3, 4, 5, 2, 4, 5]
Vals = [10, 20, 20, 30, 35, 40, 35, 50, 60, 40, 70, 80, 90, 60, 90, 100]
假定我们要从行开始，先删除第2行，则删除后的矩阵为
Rows = [0, 2, 6, 8, 10, 13]
Cols = [0, 1, 0, 1, 2, 3, 1, 3, 4, 5, 2, 4, 5]
Vals = [10, 20, 20, 30, 35, 40, 70, 80, 90, 60, 90, 100]
接下来删除第2列，删除后的矩阵为
Rows = [0, 2, 5, 7, 9, 11]
Cols = [0, 1, 0, 1, 2, 1, 2, 3, 4, 3, 4]
Vals = [10, 20, 20, 30, 40, 70, 80, 90, 90, 100]
实现
仔细比对删除前后所做的操作，在上面的示例中，我们删除的是第2行2列，这里我们定义为第id行id列。
删除指定行
在删除第id行时，首先要计算出被删行的非零数量row_nnz。
int row_nnz = Rows[id+1] - Rows[id];
下面我们来移除三个数组中被删除的元素，对于Cols数组和Vals数组，只需要通过一个for循环，不断地移除第id个数据，移除row_nnz次即可得到更新后的数组。对于Rows数组，我们需要注意的是行数组是从0开始的，因此我们要移除的数据是Rows[id+1]。
for (int i = 0; i < row_nnz; i++){ new_Vals.RemoveAt(Rows[id]); new_Cols.RemoveAt(Rows[id]);}new_Rows.RemoveAt(id + 1);
接着，我们要改变更新后元素的值，对于Cols数组和Vals数组，我们无需修改，因为每一个数据都是独立的，删除某几个数据不会其他元素造成影响。但是对于Rows数组，被删行之后的每一行的值都要减去被删行的非零元素数目个row_nnz。
for (int i = id + 1; i < new_Rows.Count; i++){ new_Rows[i] -= row_nnz;}
删除指定列
删除指定列要比删除指定行要稍麻烦些，因为这是以行为主的稀疏矩阵表示方法。同样，我们是炫耀计算出被删列的非零元素数量col_nnz，并同时找出被删元素在Cols数组中的位置， 只需要遍历每一个列元素，找到列索引等于id的元素，每找到一个，col_nnz累加1，并收集所在位置的索引。
int col_nnz = 0;List<int> col_ids = new List<int>();for (int i = 0; i < new_Cols.Count; i++){ if (new_Cols[i] == id) { col_nnz++; col_ids.Add(i); }}
接下来，我们要移除被删掉的元素，由于移除指定列对Rows数组的数量不造成影响，因此这一步仅考虑Cols数组和Vals数组。移除时需要注意，我们每次在循环中移除一个数据，后续数据的索引都会向前移动，因此我们在每次循环中要移除的实际索引是col_ids[i]-i。
for (int i = 0; i < col_nnz; i++){ new_Vals.RemoveAt(col_ids[i] -i); new_Cols.RemoveAt(col_ids[i] -i);}
下一步是更新Cols中的元素，位于每个被删元素后的元素的值都要累减1。
for (int i = 0; i < new_Cols.Count; i++) if (new_Cols[i] >= id) new_Cols[i] -= 1;
然后就是更新Rows了，我们需要知道在被删列的非零元素的所在行，只要找到一行，那么该行后面的值都要累减1，这里我使用的方法比较简单粗暴，即遍历每一行，检查该行之前的非零元素数是否大于被删元素的位置，判断为真一次则该行的值累减一次。
for (int i = 0; i < new_Rows.Count; i++)
{
int num = 0;
for (int j = 0; j < col_ids.Count; j++)
if (new_Rows[i] > col_ids[j])
num++;
new_Rows[i] -= num;
}
最后更新矩阵即可。
Rows = new_Rows.ToArray();
Cols = new_Cols.ToArray();
Vals = new_Vals.ToArray();
NNZ = Rows.Last();
N -= 1;
结语
最后，我认为该方法并不是一个十分巧妙的方法，应该有很多办法来避免做这样麻烦的操作，例如在组装矩阵一开始就算好删除元素后的位置，本文仅用作记录，代码仅供参考。
———END———
限 时 特 惠：本站每日持续更新海量各大内部创业教程，一年会员只需128元，全站资源免费下载点击查看详情
站 长 微 信：jiumai99