Seurat 与 Cellranger 之间互通的二三事中,我们遇到了 dgTMatrixdgCMatrix 这两个稀疏矩阵的不同表示。先前不清楚的时候,在必应中搜索稀疏矩阵中,出现最多的文章就是诸如 《理解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 compressedSparse(matrix, type="csc"):
""" 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]

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

参考

Comments