数据结构5.2—矩阵的压缩存储
一、整体介绍
矩阵是很多科学与工程计算问题中研究的数学对象。这里研究的不是矩阵本身，而是如何存储矩阵的元。从而使矩阵的各类运算能够有效展开。
用高级语言编制程序时，都是用二维数组来存储矩阵元。有些设计语言程序还提供了各种矩阵运算，用户使用非常方便。
存储矩阵的一般方法是采用二维数组，其优点是可以随机地访问每一个元素，因而能够较容易地实现矩阵的各种运算，如转置运算、加法运算、乘法运算等。
对于稀疏矩阵来说，采用二维数组的存储方法既浪费大量的存储单元用来存放零元素，又要在运算中花费大量的时间来进行零元素的无效计算。所以必须考虑对稀疏矩阵进行压缩存储。
尤其是在面对阶数很高的矩阵时，有时候为了节省存储空间，有必要对这类矩阵进行压缩。尤其是矩阵中存在许多值相同的元素或零元素。这就引申出一个课题：矩阵的压缩存储。
这类矩阵分两种：
特殊矩阵：假如值相同的元素或者零元素在矩阵中的分布有一定的规律
稀疏矩阵：元素分布没有规律
这里讨论的特殊矩阵都是方阵，主要有三种：
对称矩阵、上/下三角矩阵、对角矩阵；
稀疏矩阵：无法给出确切的定义，只是凭人们直觉来了解的一种概念。可以用稀疏因子去描述。
对于稀疏矩阵的压缩存储，只存储非零元的值。
二、稀疏矩阵的抽象数据类型定义及基本操作
ADT SparseMatrix{
数据对象：D={aij | i=1,2,...,m; j=1,2,...,n;
ai,j ∈ElemSet, m和n分别称为矩阵的行数和列数}
数据关系：R={Row, Col}
Row={<ai,j , ai,j+1>|1<=i<=m, 1<=j<=n-1}
Col={<ai,j , ai+1,j>|1<=i<=m-1, 1<=j<=n}
基本操作：
CreateSMatrix(&M); //创建稀疏矩阵M
DestroySMatrix(&M); //销毁稀疏矩阵M
PrintSMatrix(M); //输出稀疏矩阵M
CopySMatrix(M, &T) //由稀疏矩阵M复制得到T
AddSMatrix(M, N, &Q) //Q=M+N
SubSMatrix(M, N, &Q) //Q=M-N
MultSMatrix(M, N, &Q) //Q=M*N
TransposeMatrix(M, &T); //求稀疏矩阵M的转置矩阵T
}ADT SparseMatrix
三、稀疏矩阵压缩存储方法
这里总结一下稀疏矩阵压缩存储的几种实现方法：
COO：Coordinate 最简单的一种格式，每个非零元素需要用一个三元组来表示，分别表示行号、列号、值；
CSR ：Compressed Sparse Row ，也需要三种数据来表达：列号、值、行偏移（行偏移比较难理解）；--->一种改进型的三元组；
三元组顺序表：
一个表，其中的元素是一个三元组。像这样：（（1,2,12），（1,3,9），（3,1，-3），（3,6,14））；即COO 顺序表
三元组（i,j,a）->i行号，列号，a值；
---------稀疏矩阵的三元组顺序表存储表示-----------
1 #define MAXSIZE 12500 //假设非零元个数的最大值12500 2 typedef struct{ 3 int i,j; //非零元个数的行下标和列下标 4 ElemType e; //值 5 }Triple; //定义一个三元组 6 7 typedef struct{ 8 Triple data[MAXSIZE+1]; //非零元三元组表 data[0]未用 9 int mu,nu,tu; //行数、列数、非零元个数 10 }TSMatrix; //定义一个矩阵
问题分析：
假设a和b是TSMatrix型的变量，分别表示矩阵M和T。那么如何由a得到b呢？
分析a和b的差异即可知道：1）将矩阵的行列值相互交换；2）将每个三元组中的i和j相互交换；3）重排三元组之间的次序便可实现矩阵的转置；
前两条比较容易实现
如何实现第三条：
—————————————— ——————————————
i j v i j v
—————————————— ——————————————
1 2 12 1 3 -3
1 3 9 1 6 15
3 1 -3 2 1 12
3 6 14 2 5 18
4 3 24 3 1 9
5 2 18 3 4 24
6 1 15 4 6 -7
6 4 -7 6 3 14
—————————————— ——————————————
a.data b.data
如何使b.data中的三元组是以T的行（M的列）为主序排列。
关键就在于以T的行排列时，这就意味着以M的列为主序排列。
实现思路：
我们先定义a.data的下标为q，b.data的下标为p；
nu是TSMatrix的列数，mu是TSMatrix的行数，tu是非零元的个数，即三元组的长度；
外循环：遍历a.data的列号，即b.data的行号，即循环计数变量col，遍历的limit是M矩阵的列数M.nu；
内循环：在得到指定a.data的列号，即b.data的行号的情况下，遍历a.data三元组（a.data的下标p为循环计数变量，遍历的limit是非零元个数）
如果遍历到非零元的列号等于外循环中a.data指定的列号
将对应a.data的下标q的i，j，v赋值给对应b.data的下标p的i, j，v。（p=0不用赋值，所以p从1开始赋值）
代码实现：
1 Status TransposeSMatrix(TSMatrix M, TSMatrix &T){ 2 //采用三元组表存储表示，求稀疏矩阵M的转置矩阵T 3 T.mu=M.nu; 4 T.nu=M.mu; 5 T.tu=M.tu; 6 if(T.tu) 7 { 8 q=1; 9 for(col=1;col<=M.nu;col++) 10 for(p=1;p<=M.tu;p++) 11 if(M.data[p].j==col) 12 { 13 T.data[q].i=M.data[p].j; 14 T.data[q].j=M.data[p].i; 15 T.data[q].e=M.data[p].e; 16 ++q; 17 } 18 19 }
return OK; 20 }
算法分析：
该算法的主要工作量在p和col两重循环中完成。所以算法的时间复杂度为O(nu*tu)。即和M的列数及非零元的个数的乘积成正比。
当非零元的个数tu和mu*nu同数量级时。算法的时间复杂度为O(mu*nu*nu)了。
所以虽然节约了空间，但是时间复杂度提高了。
所以该算法仅适合于tu<<mu*nu的情况。适合非零元个数较少的情况。
接下来尝试改进上述算法，实现更快的转置。
问题分析：
矩阵M中的每一列，就是矩阵T中的每一行；
矩阵M中的每一列的第一个非零元在b.data中应有位置如果已知的话。那么在转置时，便可直接放到b.data中恰当的位置上去。
为了确定这些位置，在转置前，应先求得M的每一列中非零元的个数。
进而求得每一列的第一个非零元在b.data中应有的位置。
附设两个向量，num和cpot，num[col]代表矩阵M中第col列中非零元的个数，cpot[col]代表矩阵M中第col列的第一个非零元在b.data中的恰当位置。
cpot[1]=1;
cpot[col]=cpot[col-1]-num[col-1] 2<=col<=a.nu;
1 Status FastTransposeSMatrix(TSMatrix m, TSMatrix &T){ 2 T.mu=M.nu; T.nu=M.mu; T.tu=M.tu; 3 if(T.tu){ 4 q=1; 5 for(col=1;col<=M.nu;++col)num[col]=0; //设置num[]向量为0 6 for(t=1;t<M.tu;++t) ++num[M.data[t].j]; //初始化num[]向量 7 cpot[1]=1; 8 for(col=2;col<=M.nu;++col) cpot[col]=cpot[col-1]+num[col-1]; 9 for(p=1;p<=M.tu;++p); //遍历a.data的各个非零元 10 { 11 col=M.data[p].j; 12 q=cpot[col]; //M的某列第一个非零元数据告诉要插入到b.data的位置下标 13 T.data[q].i=M.data[p].j; 14 T.data[q].j=M.data[p].i; 15 T.data[q].e=M.data[p].e; 16 ++cpot[col];//该列的下一个非零元在b.data中的位置 17 } 18 } 19 return OK; 20 }
算法分析：
这个算法比之前的算法多用了两个辅助向量。从时间上看，算法有4个并列的单循环。
循环次数分别为nu和tu。因而总的时间复杂度为O(nu+tu)。
在M的非零元个数tu和mu*nu同等数量级时，其时间复杂度为O(mu*nu)
行逻辑链接的顺序表
为了便于随机存取任意一行的非零元，则需知道每一行的第一个非零元在三元组表中的位置。就是CSR表
这种带行链接信息的三元组表为行逻辑链接的顺序表，类型描述如下：
1 typedef struct{ 2 Triple data[]; //非零元三元组表 3 int rpos[]; //各行第一个非零元的位置表 4 int mu,nu,tu; //矩阵的行数、列数、和非零元个数 5 }RLSMatrix;
Q=M*N
M是m1*n1矩阵，N是m2*n2矩阵。当n1=m2时有：
1 for(i=1;i<=m1;++i) 2 for(j=1;j<=n2;++j) 3 { 4 Q[i][j]=0; 5 for(k=1;k<=n1;++k) Q[i][j] += M[i][k]+N[k][j]; 6 }
这个算法的时间复杂度是O(m1*n1*n2)
但是这种算法对于用三元组表作为存储结构的稀疏矩阵来说并不适用。
在上面的算法中，如论M或N的值是否为零，都要进行一次乘法运算，而实际上，这两者有一个值为零时，其乘积也为零。
因此对于稀疏矩阵来说应该免去该操作。换句话说就是在M和N中找到对应的各对元素相乘即可。找到配对的元素即可。
例如M.data[1]表示矩阵元(2,1,3)只要和N.data[1]表示的矩阵元(1,2,2)相乘。
rpos[row]指的是矩阵N的第row行中第一个非零元在N.data中的序号。
rpos[row+1]-1指的是矩阵N的第row行中最后一个非零元在N.data中的序号。
最后一行中最后一个非零元在N.data中的位置显然就是N.tu中。
乘积矩阵Q的元素是否为非零元，只有在求得其累加和后才能得知。
所以一般先将Q中第arrow行累计求和的中间结果求出来，存放到ctemp[]中，然后将ctemp[]中非零元压缩存储到Q.data；
1 Status MultiSMatrix(RLSMatrix M, RLSMatrix N, RLSMatrix &Q){ 2 //求矩阵乘积Q=M*N,采用行逻辑链接存储表示。 3 if(M.nu ！= N.mu) return ERROR; 4 Q.mu=M.mu;Q.nu=M.nu;Q.tu=0; //Q的初始化 5 if(M.tu*N.tu!-0) //Q是非零矩阵 6 { 7 for(arrow=1;arrow<=M.mu;++arrow) //处理M的每一行 8 { 9 ctemp[]=0; //当前各行元素累加清零 10 Q.rpos[arrow]=Q.tu+1; //新的一行，要确定该行非零元再Q.data中的序号 11 if(arrow<M.mu) tp=M.rpos[arrow+1]; 12 else{tp=M.tu+1;} 13 for(p=M.rpos[arrow];p<tp;++p){ //对当前行中每个非零元找到对应元在N中的行号 14 brow=M.data[p].j; 15 if(brow<N.mu) tq=N.rpos[brow+1]; 16 else{tq=N.tu+1;} 17 for(q=N.rpos[brow];q<tq;++q) 18 { 19 ccol=N.data[q].j; 20 ctemp[ccol]+=M.data[p].e * N.data[q].e; 21 } 22 } 23 for(ccol=1;ccol<=Q.nu;++ccol) //压缩存储该行非零元 24 if(ctemp[cool]){ 25 if(++Q.tu>MAXSIZE) return ERROR; 26 Q.data[Q.tu]=(arow,ccol,ctemp[ccol]); 27 } 28 } 29 } 30 return OK; 31 }
M,N,Q都是以行为主序；
求Q时，先得出一行的所有列的非零矩阵元；
十字链表
当矩阵的非零元个数和位置在操作过程中变化较大。就不宜采用顺序存储结构来表示三元组的线性表。
在链表中，每个非零元可用一个含5个域的结点表示。其中i、j和e这三个域分别表示该非零元所在的行、列和非零元的值。
向右域right用以链接同一行中下一个非零元；向下域用以链接同一列中下一个非零元；
同一行的非零元通过right域链接成一个线性链表；同一列的非零元通过down域链接成一个线性链表；
每个非零元既是某个行链表中的一个结点，有事某个列链表中的一个结点，整个矩阵就构成了一个十字交叉的链表。
故称这样的存储结构为十字链表。
可以用两个分别存储行链表的头指针和列链表的头指针的一维数组表示之。
1 typedef struct OLNode{ 2 int i,j; 3 ElemType e; 4 struct OLNode *right, *down; 5 }OLNode; *OLink; 6 7 typedef struct{ 8 OLink *rhead *chead; //行和列链表头指针向量基址由CreateSMatrix分配 9 int mu,nu,tu; 10 }CrossList;
1 Status CreateMatrix_OL (CrossList &M){ 2 //创建稀疏矩阵M。采用十字链表存储表示 3 if(M) free(M); 4 scanf(&m, &n, &t); 5 M.mu=m; N.nu=n; M.tu=t; 6 if(!(M.rhead=(OLink*)malloc((m+1)*sizeof(OLink)))) 7 exit(OVERFLOW); 8 if(!(M.chead=(OLink*)malloc((n+1)*sizeof(OLink)))) 9 exit(OVERFLOW); 10 M.rhead[]=M.chead[]=NULL; //初始化行列头指针向量；各行列链表为空链表 11 for(scanf(&i, &j, &e); i!=0; scanf(&i, &j, &e)) 12 {//输入非零元 13 if(!(p=(OLNode*)malloc(sizeof(OLNode)))) 14 exit(OVERFLOW); 15 p->i; p->j=j; p->e=e; //生成结点 16 if(M.rhead[i]=NULL||M.rhead[i]->j>j) 17 {p->right=M.rhead[i]; M.head[i]=p;} 18 else{//寻找在行表中插入的位置 19 for(q=M.rhead[i];(q->right)&&q->right->j<j;q=q->right); 20 p->right=q->right; 21 q->right=p; //完成行插入 22 } 23 if(M.chead[j]==NULL||M.chead[j]->i>i){p->down=M.chead[j];M.chead[j]=p;} 24 else{ //寻找在列表中插入的位置 25 for(q=M.chead[j];(q->down)&&q->down->i<i;q=q->down); 26 p->down=q->down; 27 q->down=p; //完成列插入 28 } 29 } 30 return OK; 31 }
现在开始讨论如何将矩阵B加到矩阵A上；
矩阵相加类似于一元多项式相加；所不同的是一元多项式中只有一个变元（指数项）。
而矩阵中每个非零元有两个变元（行值和列值）。
每个结点既在行表中又在列表中，致使插入和删除指针的修改稍为复杂，故需要更多的辅助指针。
假设两个矩阵相加后的结果为C，则矩阵C中的非零元cij只可能有三种情况；
1）它是aij+bij；
2）aij(bij=0)；
3）bij(aij=0)；
所以对于A矩阵的十字链表来说，可能是改变val域的值；或者不变（bij=0）；或者是加入一个新结点(aij=0)；
还有一种可能的情况是删除一个结点（aij+bij=0）；
由此整个运算过程可以从矩阵的第一行起逐行进行，对每一行都从行表头出发，分别找到A和B在该行的第一个非零元结点后开始比较。按照上述4种情况分别处理之。
假设非空指针pa和pb分别指向矩阵A和B中行值相同的两个结点，pa==NULL表明该行没有非零元
1）若pa==NULL或pa->j > pb->j，则需要在A矩阵的链表中插入一个值为bij的结点。此时，需改变同一行中前一个结点right域值，以及同一列中前一个结点的down域值；
2）若pa->j < pb->j，则只要将pa指针向右推进一步；
3）若pa->j==pb->j且pa->e+pb->e ！= 0，则需要将aij+bij的值送到pa所指结点的e域即可。
4）若pa->j==pb->j且pa->e+pb->e == 0，则需要在A矩阵的链表中删除pa所指的结点。此时，需改变同一行中前一个结点right域值，以及同一列中前一个结点的down域值；
为了便于插入和删除结点，还需要设立一些辅助指针。
其一是：在A的行链表上设pre指针，指示pa所指结点的前驱结点；
其二是：在A的每一列链表上设一个指针hl[j]，它的处置和列链表的头指针相同即hl[j]=chead[j]；
1 //十字链表实现矩阵相加A=A+B 2 Status AddSMatrix_OL(CrossList &A, CrossList B) 3 { 4 OLNode *pa *pb *pre 5 Olink *hl=(OLink*)malloc(A.nu*sizeof(OLink)); 6 pre=NULL; 7 for(j=1;j<=A.nu;++j) //对hl数组进行初始化，指向每一列的第一个非0元素 8 hl[j]=A.chead[j]; 9 for(int i=1;i<M.mu;i++) //按照行进行遍历 10 { 11 pa=A.rhead[i]; 12 pb=B.rhead[i]; 13 while(pb !=NULL) //当pb为NULL时，说明矩阵当前的非零元已经遍历完。 14 { 15 OLNode *p=(OLNode*)malloc(sizeof(OLNode)); //创建一个新结点，每次复制pb的结点； 16 p=pb; 17 p->down=NULL; 18 p->right=NULL; 19 if(pa==NULL||pa->j>pb->j)//第一种情况：要在pa所指结点之前插入pb所指的结点。 20 { 21 if(pre==NULL) 22 A.rhead[p->j]=p;//修改行表的头结点，如果pre为空； 23 else 24 pre->right=p;//如果pre不为空，修改pre的right域 25 p->right=pa; //插入结点的right域指向pa 26 pre=p;//更新pre，相当于pre右移了 27 if(!A.chead[p->j]||A.chead[p->j]->i>p->i) //列表down域的修改 28 { 29 p->down=M.chead[p->j]; 30 M.chead[p->j]=p; 31 } 32 else 33 { 34 p->down = hl[p->j]->down; 35 hl[p->j]->down = p; 36 } 37 hl[p->j] = p; 38 39 } 40 else 41 { 42 //第二种情况，只需要移动pa指向的位置，继续判断pa和pb；一定要有continue 43 if(pa->j < pb->j) 44 { 45 pre = pa; 46 pa = pa->right; 47 continue; 48 } 49 if(pa->j == pb->j) //第三、四种情况 50 { 51 pa->e += pa->e; 52 if(pa->e == 0) //如果为0，拆除当前结点，并释放所占空间 53 { 54 if(pre == NULL) 55 { 56 A.rhead[pa->i] = pa->right; 57 } 58 else 59 { 60 pre->right = pa->right; 61 } 62 p = pa; 63 pa = pa->right; //pa右移； 64 if(M.chead[p->j]==p) 65 { 66 M.chead[p->j] = hl[p->j] = p->down; 67 } 68 else 69 { 70 hl[p->j]->down = p->down; 71 } 72 free(p); 73 } 74 } 75 } 76 } 77 pb=pb->right; 78 } 79 return OK; 80 }
该算法的时间复杂度主要取决于A和B矩阵中非零元的个数，即ta和tb。由此算法的时间复杂度O(ta+tb)。
四、相关链接：按顺序阅读为宜
特殊矩阵的存储压缩：https://blog.csdn.net/kong_xz/article/details/79470663 //图示给得很足，适合入门了解概念，感性直观认识
数据结构-特殊矩阵的压缩存储：https://blog.csdn.net/qq_29721419/article/details/70482392 //介绍了三种特殊矩阵，及其压缩存储的代码实现。也把稀疏矩阵介绍了一下
稀疏矩阵的概念及简单实现：https://blog.csdn.net/yhb1047818384/article/details/78996906 //详细介绍了稀疏矩阵的概念，及提供了压缩方法的介绍，含图示COO和CSR，其中CSR的行偏移会比较难理解（理解成是对行坐标的一种压缩）
稀疏矩阵的三种存储方式：https://www.cnblogs.com/zhoug2020/p/7767085.html
矩阵加法基于十字链表：http://www.cnblogs.com/ciyeer/p/9039743.html
