C#数据结构(4) 稀疏矩阵与稀疏方阵

Posted 2023-04-23

参考技术A

线性代数是大学理工科学生的必修课。学过线性代数的同学一定对矩阵不陌生，因为线性代数就是一门关于矩阵的学科。

程序设计中有一种储存数据的方式是二维数组，而二维数组本质上就是矩阵。但是，假如我们想要用二维数组去储存一个大规模的矩阵并进行运算的话，会造成很大的资源浪费。举个例子，假如我们想要用二维数组去储存一个20W*20W的单位矩阵，事实上其中只有20W个数是1，其他数字都是0。所以，我们或许可以利用一种“忽略矩阵中的0项”的方式，来实现对矩阵的压缩储存，这种储存方式就叫做稀疏矩阵。对于大部分位置都是0，只有少部分位置有值的矩阵来说，使用稀疏矩阵可以让矩阵的储存密度大大提高。

我们可以使用一个 二维单向链表 来储存稀疏矩阵——因为链表在插入方面具有极低的时间复杂度，可以保证一个稀疏矩阵在经过各种运算后依然保持行与列间有序的组织。为了能够实现单向链表的删除，我们需要指向每一个有值结点前一个结点的引用。为此，我们还需要在第一行设置一个空行结点，每一行的第一列设置一个空列结点。形象地说，这种组织方式如下图所示：

使用链表的代价是寻址成本的提高，因此二维迭代器的封装对于基于链表的稀疏矩阵也是必不可少的。此外，比起二维数组矩阵，在各项操作方面，稀疏矩阵更难做到高效。因此如何高效地进行稀疏矩阵的各项操作也是本文需要探讨的话题。

本文中需要实现的稀疏矩阵基本操作如下：

同时稀疏矩阵还有派生类稀疏方阵，它实现的基本操作如下：

基于链表的稀疏矩阵的基础是行索引。因为我们定位一个矩阵中的元素，是先定位其行位置，再定位其列位置的。

因此矩阵包含了 一个作为基准的行索引链表结点 。这个结点包含一个值，代表行号；包含一个横向指针指向一个列索引链表结点（在C#中就是对一个列索引链表结点的引用），也就是该行第一个有非0元素的列；包含一个纵向指针指向下一个有非0元素的行索引链表结点。

矩阵的基础构造与行索引链表结点如下所示：

当确定矩阵中元素所在的行时，再通过列索引链表确定元素所在的列，就可以将该元素定位了。因此，列索引链表中包含两个值，一个表示自己的列号，一个表示由行索引和列索引确定的元素的值。此外，还包含一个指向和自己同行的下一个非0元素的引用。

由于二维链表在定位元素时需要先在行索引链表中寻址，再在列索引链表中寻址，其时间复杂度为 $o(n)$ ，因此在实现诸如矩阵遍历等有序操作时必须依赖二维迭代器，否则会造将大量时间浪费在寻址上，其时间复杂度将是无法想象的。 为矩阵维护一个包含一个行索引结点引用、一个列索引结点引用的对象，它将通过四个引用，分别指向一行、一列和该行的前一行、该行的前一列 ，并且通过所引用的结点内部包含的、指向其他结点的引用来实现自身所指元素的移动。

迭代器在初始化或者行复位时需要 将指向前一行的引用放置到矩阵结构的“第0行”上，将指向当前行的引用放置到第一行 。同理，每移动到下一行或者进行列复位，都要 把指向前一列的引用放置到该行的“第0列”上，将指向当前列的引用放置到第一列 。

迭代器实现了五个方法，分别是 移动到下一列、移动到下一行并指向该行第一列、移动到该行第一列、在当前指向行和前一行之间插入行、在当前指向列和前一列之间增加列 。此外，在大类中还实现了一个私有方法reLine，用来 将迭代器恢复至整个稀疏矩阵第一 行。在大类中 声明一个迭代器对象 ，用来方便省时地访问其中的元素。

利用迭代器指向矩阵中的元素，可以单独设置这个元素的内容。而如果要求设置一个不超出矩阵大小，但原本不存在（也就是为0）的位置的元素，就需要在二维链表中插入行或列。这一判断需要依靠迭代器来作出。比如，若迭代器判定其指向位置的前一行小于目标行、而当前行大于目标行，就要在当前位置前插入目标行。同理，如果把一个原本不为0的位置的元素设为0，就要删除当前列，而列为空时，则要删除当前行。

而单个元素的加减和乘除本质也是相同的。记住要分成三种情况考虑需要实现的操作：

而获取单个元素使用的方法也运用同样的思想，只不过在获取不到元素时返回元素类型的默认值。

默认初始化时，用户需要指定该矩阵的行数和列数。此外，将迭代器置于第一行第一列，虽然此时第一行第一列没有声明，但是要记得我们矩阵的链表结构中存在一个空的“第0行”。将迭代器的preline引用指向该“第0行”，其他引用则设为null。

使用其他mat初始化时，稀疏矩阵会调用目标矩阵中的迭代器和自身的迭代器来实现矩阵遍历。遍历只需要 不断调用nextRow方法，在指向的列为空时调用nextLine方法，直到指向的行为空为止即可 。

在这里强调一下reLine的作用，它是 迭代器的复位方法 ，也可以视为声明一个指向第一行第一列的迭代器。如果不进行reLine的话，假如目标矩阵中的迭代器已经指向了一个很后面的元素，那么利用迭代器nextRow和nextLine方法进行的遍历就是不成立的。

使用二维数组初始化通过遍历二维数组的方式，将二维数组中的数据一个一个set到矩阵当中。

矩阵的加减用第一个矩阵初始化一个新矩阵，同时用新矩阵和第二个矩阵的迭代器遍历两个矩阵，再在新矩阵的每个位置调用加方法，加上第二个矩阵对应位置的元素值（或其相反数），最后把新矩阵返回。

矩阵的数乘用第一个矩阵初始化一个新矩阵，同时用新矩阵的迭代器遍历新矩阵，再在新矩阵的每个位置调用乘方法，乘上目标值，最后把新矩阵返回。

依然通过迭代器来选取两个矩阵里的元素，通过矩阵乘法的特殊运算顺序来运算，将结果放进新矩阵中。返回的矩阵是通过new创建的新矩阵，因此不会对原来两个矩阵产生影响。

矩阵转置创建一个新的矩阵，其行数为原矩阵的列数，列数为原矩阵的行数，并将原矩阵中的元素所在行列对换后放进新矩阵中。

基本上就是从矩阵当中继承了初始化方法。注意构造函数是要继承自base(父类构造函数参数)的。以及行和列的数目应该相同。

都是对稀疏矩阵基本运算的重写。由于行列数目相同，实际实现会容易得多。

求行最简矩阵的目的就是把原矩阵通过初等行变换化为相抵的上三角矩阵。其具体做法是：

求完行最简矩阵之后，求行列式的值也变成了非常简单的工作。只需把对角线上所有元素相乘即可。

逆矩阵同样通过初等行变换法来解出，方法是把一个行列数与原矩阵相同的单位矩阵进行和原矩阵同步的行变换。当把原矩阵变成单位矩阵时，单位矩阵也变成了原矩阵的逆矩阵。

在 cython 中快速访问稀疏矩阵：memoryview 与字典向量

【中文标题】在 cython 中快速访问稀疏矩阵：memoryview 与字典向量

【英文标题】：fast access of sparse matrix in cython: memoryview vs vector of dictionaries

【发布时间】：2021-09-24 05:46:49

【问题描述】：

我使用 cython 来加速我在 python 中的瓶颈。任务是计算稀疏矩阵的选择性逆（低于 S），该稀疏矩阵由其以 csc 格式（数据、indptr、索引）提供的cholesky 分解给出。但任务并不是很重要，最后它是一个 3 次嵌套的 for 循环，我必须快速访问 S 的元素。

当我使用完整/巨大矩阵的内存视图时

double[:,:] Sfull

然后访问条目然后算法非常快并且符合我的期望。但很明显，这只有在矩阵 Sfull 适合内存时才有可能。

我的方法是使用字典/地图的列表/向量，这样我也可以相对快速地访问元素。

cdef vector[map[int, double]] S

事实证明，使用这种数据结构访问循环内的元素要慢大约 20 倍。这是预期的还是有其他问题？你看到其他数据结构了吗？

非常感谢您的任何 cmets 或帮助！ 最好的， 曼努埃尔

下面是 cython 代码，其中带有完整内存视图的版本被注释掉了。

cdef int invTakC12( double[:] id_diag, double[:] data, int len_i, int[:] indptr, int[:] indices, double[:, :] Sfull):


cdef vector[map[int, double]] S = testDictC(len_i-1) #list of empty dicts
cdef int i, j, j_i, lc
cdef double q 




for i in range(len_i-2, -1, -1):

    for j_i in range(indptr[i+1]-1, indptr[i]-1, -1):
        j = indices[j_i]

        q = 0

        for lc in range(indptr[i+1] -1, indptr[i], -1):

            q += data[lc] * S[j][ indices[lc] ]
            #q += data[lc] * Sfull[ indices[lc], j ]
  
        S[j][i] = -q
        #Sfull[i,j] = -q
        

        if i==j:
            S[j][i] +=  id_diag[i]
            #Sfull[i,j] += id_diag[i]
        else:
            S[i][j] -=  q
            #Sfull[j,i] -= id_diag[i]
    

return 0

【问题讨论】：

使用unordered_map 可能比使用map 更好。很多操作都更快

【参考方案1】：

您可以独立访问数组 - 例如：

cdef double[:] S_data = S.data
cdef np.int32_t[:] S_ind = S.indices
cdef np.int32_t[:] S_indptr = S.indptr

如果太不方便，可以将它们作为指针放在 C 结构中：

cdef struct mapped_csc:
    double *data
    np.int32_t *indices
    np.int32_t *indptr