4.7 矩阵的压缩存储
4.7 矩阵的压缩存储
早在线性代数中，我们就学习了关于矩阵的知识。矩阵是科学计算、工程数学，尤其是数值分析经常研究的对象。在计算机高级语言中，矩阵通常可以采用二维数组的形式来描述。
但是在有些高阶矩阵中，非零元素非常少（远小于m×n），此时若仍采用二维数组顺序存放就不合适了，因为很多存储空间存储的都是0，只有很少的一些空间存放的是有效数据，这将造成存储单元的很大浪费。另外，还有一些矩阵元素的分布有一定规律，我们可以利用这些规律，只存储部分元素，从而提高存储空间的利用率。上述矩阵叫做特殊矩阵。在实际应用中这类矩阵往往阶数很高，如2000×2000的矩阵。对于这种大容量的存储，在程序设计中，必须考虑对其进行有效的压缩存储。压缩原则是：对有规律的元素和值相同的元素只分配一个存储空间，对于零元素不分配空间。
下面介绍几种矩阵及对它们进行压缩存储的方式。
4.7.1 特殊矩阵
二维数组也称为矩阵，特殊矩阵是指非零元素或零元素的分布有一定规律的矩阵，为了节省存储空间，特别是在高阶矩阵的情况下，可以利用特殊矩阵的规律，对它们进行压缩存储，也就是说，使多个相同的非零元素共享同一个存储单元，对零元素不分配存储空间。特殊矩阵的主要形式有对称矩阵、对角矩阵等，它们都是方阵，即行数和列数相同。
1．对称矩阵的压缩存储
若一个n阶方阵A[n][n]中的元素满足ai,j=aj, i（0≤i，j≤n-1），则称其为n阶对称矩阵。
由于对称矩阵中的元素关于主对角线对称，因此在存储时可只存储对称矩阵中上三角或下三角中的元素，使得对称的元素共享一个存储空间。这样，就可以将n2个元素压缩存储到 
个元素的空间中。不失一般性，以行序为主序存储其下三角（包括对角线）的元素。
假设以一维数组B[0..
-1]作为n阶对称矩阵A的存储结构，在B中只存储对称矩阵A的下三角元素ai, j（i≥j），则上三角元素aj, i在B数组中对应的元素就是ai, j在B数组中对应的元素。对于元素ai,j，不包括当前行，它前面共有i行（行下标为0～i-1，行下标为0的行有1个元素，行下标为1的行有2个元素，……，行下标为i-1的行有i个元素），则第i行有1+2+…+i=
个元素，在当前行中，元素ai, j的前面亦有j个元素，则A中任一元素ai, j和bk之间存在着如下对应关系：
有些非对称的矩阵也可借用此方法存储，如n阶下（上）三角矩阵。所谓n阶下（上）三角矩阵，是指矩阵的上（下）三角（不包括对角线）中的元素均为常数c或0的n阶方阵。设以一维数组B[0..
]作为n阶三角矩阵A的存储结构，则A中任一元素ai, j和B中元素bk之间存在着如下对应关系。
上三角矩阵：
下三角矩阵：
其中，B的元素bn(n+1)/2中存放着常数c。
2．对角矩阵的压缩存储
|
图4-7 半带宽为b的对角矩阵
|
）的对角矩阵，其|i-j|≤b的元素ai, j不为零，其余元素为零。如图4-7是半带宽为b的对角矩阵示意图。
对于b=1的三对角矩阵，只存储其非零元素到一维数组B中，即将A的非零元素ai,j存储到B的元素bk中。A中第0行和第n-1行都只有两个非零元素，其余各行有3个非零元素。对于不在第0行的非零元素ai,j来说，在它前面存储了矩阵的前i行元素，这些元素的总数为2+3(i-1)。
· 若ai, j是本行中需要存储的第1个元素，则k=2+3(i-1)=3i-1，此时j=i-1，即k=2i+i-1=2i+j。
· 若ai, j是本行中需要存储的第2个元素，则k=2+3(i-1)+1=3i，此时i=j，即k=2i+i=2i+j。
· 若ai, j是本行中需要存储的第3个元素，则k=2+3(i-1)+2=3i+1，此时j=i+1，即k=2i+i+1=2i+j。
归纳起来有：k=2i+j。
以上讨论的对称矩阵、对角矩阵、三对角矩阵的压缩存储方法是将有一定分布规律的值相同的元素（包括0）压缩存储到一个存储空间中。这样的压缩存储只需在算法中按公式作一映射即可实现矩阵元素的随机存取。
4.7.2 稀疏矩阵
一个阶数较大的矩阵中的非零元素个数s相对于矩阵元素的总个数t非常小时，即s<<t时，称该矩阵为稀疏矩阵。例如，一个100×100的矩阵，若其中只有100个非零元素，就可以称其为稀疏矩阵。
稀疏矩阵与抽象数据类型d（d=2）维数组定义相似，这里不再介绍。
1．稀疏矩阵的三元组表示
不同于以上讨论的几种特殊矩阵的压缩存储方法，稀疏矩阵的压缩存储方法是只存储非零元素。由于稀疏矩阵中非零元素的分布没有任何规律，所以在存储非零元素时还必须同时存储该非零元素所对应的行下标和列下标。这样稀疏矩阵中的每一个非零元素需由一个三元组(i,j,ai,j)唯一确定，稀疏矩阵中的所有非零元素构成三元组线性表。
假设有一个6×7阶稀疏矩阵A（为图示方便，我们所取的行列数都很小），A中元素如图4-8所示。
图4-8 稀疏矩阵A
则对应的三元组线性表为：
((0,2,1)，(1,1,2)，(2,0,3)，(3,3,5)，(4,4,6)，(5,5,7)，(5,6,4))
若把稀疏矩阵的三元组线性表按顺序存储结构存储，则称为稀疏矩阵的三元组顺序表。则三元组顺序表的数据结构可定义如下：
#define MaxSize 100 /*矩阵中非零元素的最多个数*/
typedef struct
{
int r; /*行号*/
int c; /*列号*/
ElemType d; /*元素值*/
} TupNode; /*三元组定义*/
typedef struct
{
int rows; /*行数值*/
int cols; /*列数值*/
int nums; /*非零元素个数*/
TupNode data[MaxSize];
} TSMatrix; /*三元组顺序表定义*/
其中，data域中表示的非零元素通常以行序为主序顺序排列，它是一种下标按行有序的存储结构。这种有序存储结构可简化大多数矩阵运算算法。下面的讨论假设data域按行有序存储。
矩阵运算通常包括矩阵转置、矩阵加、矩阵减、矩阵乘等。这里仅讨论基本运算和矩阵转置运算算法。
（1）从一个二维矩阵创建其三元组表示的运算算法。
以行序方式扫描二维矩阵A，将其非零的元素插入到三元组t的后面。算法如下：
void CreatMat(TSMatrix &t,ElemType A[M][N])
{
int i,j;
t.rows=M;t.cols=N;t.nums=0;
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
if(A[i][j]!=0) /*只存储非零元素*/
{
t.data[t.nums].r=i;t.data[t.nums].c=j;
t.data[t.nums].d=A[i][j];t.nums++;
}
}
}
（2）三元组元素赋值运算算法。
先在三元组t中找到适当的位置k，将k～t.nums个元素后移一位，将指定元素x插入到t.data[k]处。算法如下：
int Value(TSMatrix &t,ElemType x,int rs,int cs)
{
int i,k=0;
if(rs>=t.rows || cs>=t.cols)
return 0;
while(k<t.nums && rs>t.data[k].r) k++; /*查找行*/
while(k<t.nums && cs>t.data[k].c) k++; /*查找列*/
if (t.data[k].r==rs && t.data[k].c==cs) /*存在这样的元素*/
t.data[k].d=x;
else /*不存在这样的元素时插入一个元素*/
{ for (i=k;i<t.nums;i++)
{
t.data[i+1].r=t.data[i].r;
t.data[i+1].c=t.data[i].c;
t.data[i+1].d=t.data[i].d;
}
t.data[k].r=rs;t.data[k].c=cs;t.data[k].d=x;
t.nums++;
}
return 1;
}
（3）将指定位置的元素值赋给变量的运算算法。
先在三元组t中找到指定的位置，将该处的元素值赋给x。算法如下：
int Assign(TSMatrix t,ElemType &x,int rs,int cs)
{
int k=0;
if(rs>=t.rows || cs>=t.cols)
return 0;
while(k<t.nums && rs>t.data[k].r) k++;
while(k<t.nums && cs>t.data[k].c) k++;
if(t.data[k].r==rs && t.data[k].c==cs)
{
x=t.data[k].d;
return 1;
}
else return 0;
}
（4）输出三元组运算算法。
从头到尾扫描三元组t，依次输出元素值。算法如下：
void DispMat(TSMatrix t)
{
int i;
if(t.nums<=0) return;
printf(\t%d\t%d\t%d\n",t.rows,t.cols,t.nums);
printf(" ------------------\n");
for(i=0;i<t.nums;i++)
printf("\t%d\t%d\t%d\n",t.data[i].r,t.data[i].c,t.data[i].d);
}
（5）矩阵转置运算算法。
对于一个m×n的矩阵Am×n，其转置矩阵是一个n×m的矩阵。设为Bn×m，满足ai, j=bj, i，其中1≤i≤m，1≤j≤n。其完整的转置算法如下：
void TranTat(TSMatrix t,TSMatrix &tb)
{
int p,q=0,v; /*q为tb.data的下标*/
tb.rows=t.cols;tb.cols=t.rows;tb.nums=t.nums;
if (t.nums!=0)
{
for (v=0;v<t.cols;v++) /*tb.data[q]中的记录以c域的次序排列*/
for(p=0;p<t.nums;p++) /*p为t.data的下标*/
if(t.data[p].c==v)
{
tb.data[q].r=t.data[p].c;
tb.data[q].c=t.data[p].r;
tb.data[q].d=t.data[p].d;
q++;
}
}
}
以上算法的时间复杂度为O(t.cols*t.nums)，而将二维数组存储在一个m行n列的矩阵中时，其转置算法的时间复杂度为O(m*n)。最坏情况是当稀疏矩阵中的非零元素个数t.nums和m*n同数量级时，上述转置算法的时间复杂度就为O(m*n2)。对其他几种矩阵运算也是如此。可见，常规的非稀疏矩阵应采用二维数组存储，只有当矩阵中非零元素个数s满足s<<m*n时，方可采用三元组顺序表存储结构。这个结论也同样适用于下面要讨论的三元组的十字链表。
【例4-2】采用三元组的存储结构，设计两个矩阵相加的运算算法。
解：用i和j作为两个三元组的指针，以行序对它们的当前值进行相加运算，并将结果存放在c中。本例算法如下。
int MatAdd(TSMatrix a,TSMatrix b,TSMatrix &c)
{
int i=0,j=0,k=0;
ElemType v;
if(a.rows!=b.rows || a.cols!=b.cols)
return 0; /*行数或列数不等时不能进行相加运算*/
c.rows=a.rows;c.cols=a.cols; /*c的行列数与a的相同*/
while(i<a.nums && j<b.nums) /*处理a和b中的每个元素*/
{
if (a.data[i].r==b.data[j].r) /*行号相等时*/
{ if(a.data[i].c<b.data[j].c)/*a元素的列号小于b元素的列号*/
{
c.data[k].r=a.data[i].r;/*将a元素添加到c中*/
c.data[k].c=a.data[i].c;
c.data[k].d=a.data[i].d;
k++;i++;
}
else if (a.data[i].c>b.data[j].c)/*a元素的列号大于b元素的列号*/
{
c.data[k].r=b.data[j].r; /*将b元素添加到c中*/
c.data[k].c=b.data[j].c;
c.data[k].d=b.data[j].d;
k++;j++;
}
else /*a元素的列号等于b元素的列号*/
{
v=a.data[i].d+b.data[j].d;
if(v!=0) /*只将不为0的结果添加到c中*/
{
c.data[k].r=a.data[i].r;
c.data[k].c=a.data[i].c;
c.data[k].d=v;
k++;
}
i++;j++;
}
}
else if(a.data[i].r<b.data[j].r) /*a元素的行号小于b元素的行号*/
{
c.data[k].r=a.data[i].r; /*将a元素添加到c中*/
c.data[k].c=a.data[i].c;
c.data[k].d=a.data[i].d;
k++;i++;
}
else /*a元素的行号大于b元素的行号*/
{
c.data[k].r=b.data[j].r; /*将b元素添加到c中*/
c.data[k].c=b.data[j].c;
c.data[k].d=b.data[j].d;
k++;j++;
}
c.nums=k;
}
return 1;
}
2．稀疏矩阵的十字链表表示
十字链表为稀疏矩阵的每一行设置一个单独链表，同时也为每一列设置一个单独链表。这样稀疏矩阵的每一个非零元素就同时包含在两个链表中，即每一个非零元素同时包含在所在行的行链表中和所在列的列链表中。这就大大降低了链表的长度，方便了算法中行方向和列方向的搜索，因而大大降低了算法的时间复杂度。
对于一个m×n的稀疏矩阵，每个非零元素用一个节点表示，节点结构可以设计成如 图4-9（a）所示的结构。其中i，j，value分别代表非零元素所在的行号、列号和相应的元素值；down和right分别称为向下指针和向右指针，分别用来链接同列中和同行中的下一个非零元素节点。也就是说，稀疏矩阵中同一列的所有非零元素通过down指针链接成一个列链表，同一行中的所有非零元素通过right指针链接成一个行链表。对稀疏矩阵的每个非零元素来说，它既是某个行链表中的一个节点，同时又是某个列链表中的一个节点。每个非零元素就好比在一个十字路口，由此称做十字链表。
十字链表中设置行头节点、列头节点和链表头节点。它们采用和非零元素节点类似的节点结构，具体如图4-9（b）所示。其中行头节点和列头节点的i，j域值均为0；行头节点的right指针指向该行链表的第一个节点，它的down指针为空；列头节点的down指针指向该列链表的第一个节点，它的right指针为空。行头节点和列头节点必须顺序链接，这样当需要逐行（列）搜索时，才能一行（列）搜索完后顺序搜索下一行（列），行头节点和列头节点均用link指针完成顺序链接。对比行头节点和列头节点可以看到，行头节点中未用down指针，列头节点中未用right指针，link指针完成行或列头节点的顺序链接，i域和j域未用。因此行和列的头节点可以合用，即第i行和第i列头节点共用一个头节点。这些合并后的头节点称为行列头节点，行列头节点数为矩阵行数m和列数n的最大值。
（a）节点结构 （b）头节点结构
图4-9 十字链表节点结构
十字链表头指针head指向链表头节点，链表头节点的i，j域分别存放稀疏矩阵的行数m和列数n，链表头节点的link指针指向行列头节点链表中的第一个行列头节点。由于矩阵运算中常常是一行（列）操作完后进行下一行（列）操作，所以十字链表中的所有单链表均链接成循环链表。这样就可方便地完成一行（列）操作后又回到该行列头节点，由link指针进入下一行列头节点，重新开始下一行（列）的相同操作。
设稀疏矩阵为：
则对应上述的十字链表如图4-10所示。为图示方便清楚，把每个行列头节点分别画成两个，而实际上，行头节点hi（0≤i≤4）与列头节点hi只存在一个这样的节点。
图4-10 对应稀疏矩阵的十字链表
十字链表节点结构和头节点的数据结构可定义如下：
#define M 3 /*矩阵行*/
#define N 4 /*矩阵列*/
#define Max((M)>(N)?(M)：(N)) /*矩阵行列较大者*/
typedef struct mtxn
{
int row; /*行号*/
int col; /*列号*/
struct mtxn *right,*down; /*向右和向下的指针*/
union
{
int value;
struct mtxn *link;
} tag;
} MatNode; /*十字链表类型定义*/
有关十字链表表示的运算算法与三元组表示类似，这里不再讨论。