矩阵的压缩存储

目录

       矩阵

    特殊矩阵

    稀疏矩阵

---

       矩阵是很多科学与工程计算问题中研究的数学对象。通常在编写程序时，用二维数组来表示矩阵，最简单，最直观，用于实现各种矩阵运算也最方便。

       矩阵

       这里首先给出一个以二维数组为存储结构的矩阵的类型实现，实现了矩阵的转置和乘法运算。

class matrix:
    def __init__(self, m):
        self.m = m
        self.rn = len(m)
        self.cn = len(m[0])

    def transpose(self):
        m = self.m
        r = [[0] * self.rn for i in range(self.cn)]
        for i in range(self.rn):
            for j in range(self.cn):
                r[j][i] = m[i][j]
        return self.__class__(r)

    def __str__(self):
        s = ""
        for i in self.m:
            s += " ".join([str(j) for j in i]) + "\n"
        return s

    def __mul__(self, r):
        n = [[0] * r.cn for i in range(self.rn)]
        m = self.m
        _m = r.m
        for i in range(self.rn):
            for j in range(r.cn):
                for k in range(self.cn):
                    n[i][j] += m[i][k] * _m[k][j]
        return self.__class__(n)
                
m = matrix([[3, 0, 0, 5], [0, -1, 0, 0], [2, 0, 0, 0]])
n = matrix([[0, 2], [1, 0], [-2, 4], [0, 0]])
t = m.transpose()
print(m, n, t, m * n)

        然而，在数值分析中经常出现一些阶数很高的矩阵，同时矩阵中有许多值相同的元素或者零元素。有时为了节省存储空间，对这类矩阵进行压缩存储。所谓压缩存储是指，为多个值相同的元只分配一个存储空间，对零元不分配空间。

       假若值相同的元素或者零元素在矩阵中的分布有一定规律，我们称此类矩阵为特殊矩阵，反之，称之为稀疏矩阵。

    特殊矩阵

       若n阶矩阵A中的元满足性质：aij = aji 对取0-n之间的i, j均成立，则称之为n阶对称矩阵。对于对称矩阵，我们可以为每一对对称元分配一个存储空间，则可将n^2个元压缩存储到n*(n+1)/2个元的空间中。

       不失一般性，假定对称矩阵是下三角矩阵，我们以行进行分割表示如下：a00;a10 a11;a20 a21 a22;……;an0 an1…ann。假设以一维数组sa作为对称矩阵的存储结构，并且是依行号由小到大逐行存储，那么对于任意的aij，它在数组sa中的index是多少？这是个数学问题就是求f((i, j)) = ?，这个问题并不复杂，求aij的位置，首先要求ai0在sa中的起始位置，它等于前i行的元个数的和：1+2+3+…+i = i*(i+1)/2，那么aij在sa中的位置为：f((i, j)) = i*(i+1)/2+j （i<=j)，若i > j时，f((i, j)) = j*(j+1)/2+i。

       其他类型的一些特殊矩阵，非零元的分布都有一个明显的规律，只要能找到一个{(i, j)|i, j in[1-n]}和自然数集N的一个首部子集之间的一一映射，就可以将矩阵存储在一维数组中。

    稀疏矩阵

        什么是稀疏矩阵？假设在一个m行n列的矩阵中，有t个非零元，非零元占比t/(m*n)小于等于5%，这样的矩阵被称之为稀疏矩阵。如何进行稀疏矩阵的压缩存储呢？仍然是存储到一个一维数组中，但它不是直接存储矩阵元的值了，而是存储一个三元组(i, j, v)，i、j表示矩阵元在矩阵中的行、列号，v则是矩阵元的值了。

      来看看三元组顺序表表示的稀疏矩阵的转置和乘法如何实现。

             3  0  0  5                       3  0   2

    M =  0 -1  0  0              M’ =  0  -1  0

            2  0  0  0                         0  0   0

                                                    5  0   0

       例如求M的转置矩阵M’，最简单的一种办法就是通过一轮对M的三元组顺序表遍历得到M’的一行，也就是以M’的行号（也是M的列号）去M中检索符合条件的元的三元组，把三元组的行列号对掉，再把它们存到M’的三元组顺序表中。以M’的行号由0到n的顺序依次去M的三元组顺序表中遍历检索符合条件的元的三元组，处理之后追加到M’的三元组顺序表中。

        这个算法的思路还是很清晰的，但是效率比较差。

        改进的一个算法是通过分析M得倒它的一些特征值来辅助M的转置，M的一列对应M’的一行，可以先分析M，得倒M的每一列包含的非零元的数目，假设用一个数组col_ele_counts来表示，进而根据它得到M的每一列，也就是M’的每一行的第一个非零元在M’的三元组顺序表中的起始index，假设用一个数组col_first_indexes来存这些index，这里存在一个递推关系：col_first_indexes[0] = 0，col_first_indexes[i] = col_first_indexes[i-1] + col_ele_counts[i-1]，有了col_first_indexes，完成M到M’的转换就简单了，遍历M就可以完成，比方说遍历得到一个元的三元组(i, j, v)，根据col_first_indexes[j]找到它在M’的三元组顺序表的填入index，然后把三元组(j, i, v)存到那个位置，然后col_first_indexes[j]++就可以了。

      下面再来看如何实现两个用三元组顺序表表示的稀疏矩阵的乘法。

             3   0  0  5                0   2                  0   6

    M =  0  -1  0  0        N =  1   0         Q =  -1   0

             2   0  0  0               -2   4                 0   4

                                             0   0

       M是m1行n1列的矩阵，N是m2行n2列的矩阵，M和N满足条件，M的列数等于N的行数时，也就是n1等于m2，M和N可以相乘得到一个结果矩阵Q，Q的行数为m1，列数为n2，如果拿M的第i行和N的第j列相乘得到cij有cij = ai0 * b0j + ai1 * b1j + … + ain1-1 * bn1-1j，要得到Q的一行：ci0 ci1 … cin2怎么做呢？拿M的第i行和N的所有列去相乘，要得到Q，怎么做呢？依次拿M的行去和N的所有列相乘。

       但在矩阵的三元组顺序表的表示里，非零元是按行来进行存储的，要得到矩阵的一列并不方面，所以直接实现M的行和N的列的相乘并不容易，那该怎么办呢？一种办法是把N转置为N’，然后对M和N’实行行和行的相乘，那不是很容易？这样确实可以，但需要对N转置，不转置行不行？它给了我们一个启发，那就是也许可以从行相乘的角度去考虑解决矩阵相乘的问题。在做M的行和N的列相乘时发现，M行中的元只和N列中的特定元结合，也就是M的第i行j列的元aij只和N的任意一列的第j行的元进行结合，既如此，我们可以让aij和N的所有列的第j行的元（也就是N的第j行）相乘得到一个序列：aij * bj0, aij * bj1, … , aij * bjn2，序列中的元素分别为ci0, ci1, …, cin2（Q的第i行）的值的一部分，可以预设一个初始值为0的数组larr来存放Q的一行的各个矩阵元的值，序列中的项aij * bjk加到larr[k]中（larr[k] += aij * bjk)，我们按照列号从小到大的顺序依次让M中第i行的元aij和N的第j行进行上面的相乘处理，最后就能得到Q的第i行，我们按照行号从小到大的顺序依次对M的行执行上述的操作，便能按照行号从小到大的顺序依次得到Q的各行，也就能得到完整的Q了。

import copy
class seqmatrix:
    def __init__(self, m, rn, cn, recs = None, rtab=None):
        self.m = m
        self.rn =rn
        self.cn = cn
        self.en = len(m)
        if recs is None:
            recs = [0] * rn
            for i in m:
                recs[i[0]] += 1
        if rtab is None:
            rtab = [0] * rn
            for i in range(1, len(rtab)):
                rtab[i] = rtab[i - 1] + recs[i - 1]
        self.recs = recs
        self.rtab = rtab

    def transpose(self):
        cecs = [0] * self.cn
        ctab = [0] * self.cn
        m = [()] * self.en
        for e in self.m:
           cecs[e[1]] += 1
        for i in range(1, len(ctab)):
            ctab[i] = ctab[i - 1] + cecs[i - 1]
        _ctab = copy.copy(ctab)
        for e in self.m:
            m[ctab[e[1]]] = (e[1], e[0], e[2])
            ctab[e[1]] += 1
        return self.__class__(m, self.cn, self.rn, cecs, _ctab)

    def __str__(self):
        m = [[0] * self.cn for i in range(self.rn)]
        for e in self.m:
            m[e[0]][e[1]] = e[2]
        return matrix(m).__str__()

    def __mul__(self, r):
        if not self and r:
            return self.__class__([], self.rn, r.cn)
        m = []
        for i in range(self.rn):
            lcache = [0] * r.cn
            for j in range(self.rtab[i], self.rtab[i + 1] if i < self.rn - 1 else self.en):
                l = self.m[j][1]
                for k in range(r.rtab[l], r.rtab[l + 1] if l < r.rn - 1 else r.en):
                    lcache[r.m[k][1]] += self.m[j][2] * r.m[k][2]

            for _j, v in enumerate(lcache):
                if v:
                    m.append((i, _j, v))
        return self.__class__(m, self.rn, r.cn)

m = seqmatrix([(0, 0, 3), (0, 3, 5), (1, 1, -1), (2, 0, 2)], 3, 4)
t = m.transpose()
print(m, t)
r = seqmatrix([(0, 1, 2), (1, 0, 1), (2, 0, -2), (2, 1, 4)], 4, 2)
print(r, m * r)