稀疏矩阵的表示法

在Seurat 与 Cellranger 之间互通的二三事中，我们遇到了 dgTMatrix 和 dgCMatrix 这两个稀疏矩阵的不同表示。先前不清楚的时候，在必应中搜索稀疏矩阵中，出现最多的文章就是诸如 《理解Compressed Sparse Column Format (CSC)》这一类文章，我就不吐槽 CSDN 了，唉。这也就说明了作为写代码的人，你为什么要去直接看英文的原始资料，因为你永远不知道，翻译资源的那个人英文水平怎么样，更别提他的技术水平了。我这里是为了记录自己的学习内容，目的是为了自己观看。如果有某位同学进来了，还是看英文资料比较好一点，放在 参考 中了，你可以自己查阅。

稀疏矩阵的表示法

所谓稀疏矩阵，就是矩阵包含的值有太多是 0 了，而关键信息的值却占比很少。如果把这些 0 也进行存储的话，无疑是浪费了太多空间。基于此，就有了稀疏矩阵的各种表示格式，以摘要的信息表示那些非零值在矩阵中的位置。总体上讲，稀疏矩阵的表示格式可以分为两大类：支持高效修改的、支持高效访问与矩阵操作的。

• Efficient modification
  • Dictonary of keys
  • List of lists
  • Coordinate list
• Efficient access & matrix operations
  • compressed sparse row
  • compressed sparse column

这里我们主要看一下后面的两种：CSR 和 CSC。其实它们两个是相似的，应该是转置的关系。所以只要了解 CSR 是怎么得到的，CSC 也能轻易了解；反之亦然。

Compressed sparse row format, CSR

1 2 3 4

0 0 0 0 5 8 0 0 0 0 3 0 0 6 0 0

对于上面的这个矩阵，我们先从左到右，从上到下获得这几个信息：

• 非零值 A: [5, 8, 3, 6]
• 每行的非零值个数 NNZ：[0, 2, 1, 1]
• 每个非零值的所在列 cols：[0, 1, 2, 1] ( 基于 0 开始的索引）
• rows: []

因为 CSR 是针对 row 的信息进行压缩的，所以我们要使用 NNZ 这个数据来获得每个非零值的所在行信息。首先我们看 0 + 2 + 1 + 1 = 4 是 A 的长度，而且是按每行进行排的，熟悉数组的同学应该就可以理解下面这个操作：例如第二行的元素可以通过 第一行的非零元素个数(这里是NNZ[0]=0) + 1作为切片的 start , start + 第二行的非零元素个数(NNZ[1]=2) 作为切片的 end，对非零值 A 集合进行切片得到。也就是说对于第 i 行，它的非零值应该为 A[ NNZ[i-1]+1, NNZ[i-1]+1+NNZ[i] ], 值对应的所在列应该为 A[ NNZ[i-1]+1, NNZ[i-1]+1+NNZ[i] ]。这样通过三个一维数组就可以得到稀疏矩阵的表示。但是，如果像前面那样做，我们每次又得在 end 加上 start，有点麻烦，所以对于 NNZ，我们可以在生成的时候直接获得，因而 NNZ = [0, 2, 3, 4] 。 到了这里，上述矩阵的稀疏矩阵表示法就是:

1 2 3

non_zero = [5, 8, 3, 6] # 非零值 col_name = [0, 1, 2, 1] # 非零值所在的列 row_indices = [0, 2, 3, 4] # 每行对应的非零值指针；对于第 0 行，取的切片应该为 0:row_indices[0]

基于上面还有额外添加 0 的信息，则可以直接把 0 加入到 row_indices 中，在进行切片的时候，变成了 row_indices[i]: row_indices[i+1].

1 2 3

non_zero = [5, 8, 3, 6] # 非零值 col_name = [0, 1, 2, 1] # 非零值所在的列 row_indices = [0, 0, 2, 3, 4] # 各个行的非零值指针

下面我们就使用 Python 来编写函数实现这个逻辑（不用 R 的原因在于还要处理切片的问题，1-based index）的问题。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

def (matrix, type="csc"): & 大专栏  稀疏矩阵的表示法#34;"" create compressed sparse column format (csc) or row format (csr)""" non_zero = [] indices = [] pointer = [0] count = 0 if type == "csr": for i in range(len(matrix)): for j in range(len(matrix[i])): if matrix[i][j] != 0: non_zero.append(matrix[i][j]) indices.append(j) # 添加所在列 count += 1 # 计数 pointer.append(count) # 每行都要添加非零值指针 else: for j in range(len(matrix[0])): for i in range(len(matrix)): if matrix[i][j] != 0: non_zero.append(matrix[i][j]) indices.append(j) # 添加所在列 count += 1 # 计数 pointer.append(count) # 每行都要添加非零值指针 description = "compressed sparse {type} format, x is non zeros, indices is the {type} pointer, line_num is the {type} line number in the matrix".format(type="row" if type == "csr" else "column") return { "x": non_zero, "indices": indices, "pointer": pointer, "type": type, "description": description "dim": [len(matrix), len(matrix[0])] } }

基于上面的函数得到的 Compressed sparse column format 为：

1 2 3

x [5, 8, 6, 3] indices [1, 1, 3, 2] pointer [0, 1, 3, 4, 4]

而 R 中得到的 Compressed sparse column format 为：

1 2 3 4 5 6 7 8 9

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:4] 1 1 3 2 ..@ p : int [1:5] 0 1 3 4 4 ..@ Dim : int [1:2] 4 4 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:4] 5 8 6 3 ..@ factors : list()

同样的，我们看一下 csr 是否同上面分析的一致。

1 2 3

x [5, 8, 3, 6] indices [0, 1, 2, 1] pointer [0, 0, 2, 3, 4]

都是一致，这表明，我们的结果是对的。

参考

• sparse matrix compression matrix
• 维基-sparse matrix