矩阵的压缩存储

§5.3 矩阵的压缩存储 矩阵是很多科学与工程计算问题中研究的数学对象，在此，我们讨论如何存储矩阵的元，从而使矩阵的各种运算能有效第进行。对于一个矩阵结构显然用一个二维数组来表示是非常恰当的，但在有些情况下，比如常见的一些特殊矩阵，如三角矩阵、对称矩阵、带状矩阵、稀疏矩阵等，从节约存储空间的角度考虑，对这类矩阵进行压缩存储，即：为多个值相同的元素或只分配一个存储单元，对零元素不分配空间。下面从这一角度来考虑这些特殊矩阵的存储方法。 值相同的元素或者零元素在矩阵中的分布有一定的规律，则称此类矩阵为特殊矩阵；反之为稀疏矩阵。 5.3.1 特殊矩阵 对称矩阵: 对称矩阵的特点是：在一个n阶方阵中，有aij=aji ，其中1≤i , j≤n，如图5.5所示是一个5阶对称矩阵。对称矩阵关于主对角线对称，因此只需存储上三角或下三角部分即可，比如，我们只存储下三角中的元素aij，其特点是中j≤i且1≤i≤n，对于上三角中的元素aij ，它和对应的aji相等，因此当访问的元素在上三角时，直接去访问和它对应的下三角元素即可，这样，原来需要n*n个存储单元，现在只需要n(n+1)/2个存储单元了，节约了n(n-1)/2个存储单元，当n较大时，这是可观的一部分存储资源。 如何只存储下三角部分呢？对下三角部分以行为主序顺序存储到一个向量中去，在下三角中共有n*(n+1)/2个元素，因此，不失一般性，设存储到向量SA[n(n+1)/2]中，存储顺序可用图5.6示意，这样，原矩阵下三角中的某一个元素aij则具体对应一个sak，下面的问题是要找到k与i、j之间的关系。 对于下三角中的元素aij，其特点是：i≥j且1≤i≤n，存储到SA中后，根据存储原则，它前面有i-1行，共有1+2+…+i-1=i*(i-1)/2个元素，而aij又是它所在的行中的第j个，所以在上面的排列顺序中，aij是第i*(i-1)/2+j个元素，因此它在SA中的下标k与i、j的关系为： k=i*(i-1)/2+j-1 (0≤k<n*(n+1)/2 ) 若i<j，则aij是上三角中的元素，因为aij=aji ，这样，访问上三角中的元素aij时则去访问和它对应的下三角中的aji即可，因此将上式中的行列下标交换就是上三角中的元素在SA中的对应关系： k=j*(j-1)/2+i-1 (0≤k<n*(n+1)/2 ) 综上所述，对于对称矩阵中的任意元素aij，若令I=max(i,j)，J=min(i,j)，则将上面两个式子综合起来得到：　k=I*(I-1)/2+J-1。称SA[n(n+1)/2]为n阶对称矩阵A的压缩存储。 三角矩阵 形如图5.7的矩阵称为三角矩阵，其中c为某个常数。其中5.7(a)为下三角矩阵：主队角线以上均为同一个常数；(b)为上三角矩阵，主队角线以下均为同一个常数；下面讨论它们的压缩存储方法。 下三角矩阵 与对称矩阵类似，不同之处在于存完下三角中的元素之后，紧接着存储对角线上方的常量，因为是同一个常数，所以存一个即可，这样一共存储了n*(n+1)+1个元素，设存入向量：SA[n*(n+1)+1]中，这种的存储方式可节约n*(n-1)-1个存储单元，sak 与aji 的对应关系为： 2. 上三角矩阵 对于上三角矩阵，存储思想与下三角类似，以行为主序顺序存储上三角部分，最后存储对角线下方的常量。对于第1行，存储n个元素，第2行存储n-1个元素，…，第p行存储(n-p+1)个元素，aij的前面有i-1行，共存储： n+(n-1) +…+(n-i+1)= =(i-1)*(2n-i+2)/2个元素，而aij 是它所在的行中要存储的第（j-i+1）个；所以，它是上三角存储顺序中的第 (i-1)*(2n-i+2)/2+(j-i+1)个，因此它在SA中的下标为：k=(i-1)*(2n-i+2)/2+j-i。 综上， sak 与aji 的对应关系为： 对角矩阵 在这种矩阵中，所有的非零元都集中在以主对角线为中心的带状区域中，即除了主对角线上和直接在对角线上方、下方若干条对角线上的元外，所有的其它元皆为零。对于这种矩阵，我们也可按某个原则（以行为主，或以对角线为主的顺序）将其压缩存储道一维数组中。见P96 图5.4 5.3.2 稀疏矩阵 一、定义： 设m*n矩阵中有t个非零元素且t<<m*n，这样的矩阵称为稀疏矩阵。令,称为矩阵的稀疏因子，通常认为≤0.05时称为稀疏矩阵。很多科学管理及工程计算中，常会遇到阶数很高的大型稀疏矩阵。 二、抽象数据类型稀疏矩阵的定义 ADT SparseMatrix{ 数据对象： 数据关系： 基本操作: CreateSMatrix(&M) 创建稀疏矩阵 DestroySMatrix(&M) 销毁稀疏矩阵 PrintSMatrix(&M) 输出稀疏矩阵 CopySMatrix(M，&T) 复制稀疏矩阵 AddSMatrix(M，N，&Q) 求稀疏矩阵的和 SubSMatrix(&M) 求稀疏矩阵的差 MultSMatrix(&M) 求稀疏矩阵的积 TransposeSMatrix(&M) 求稀疏矩阵的转置矩阵 } ADT SparseMatrix 三、稀疏矩阵的三元组表存储 如果按常规分配方法，顺序分配在计算机内，那将是相当浪费内存的。为此提出另外一种存储方法，仅仅存放非零元素。但对于这类矩阵，通常零元素分布没有规律，为了能找到相应的元素，所以仅存储非零元素的值是不够的，还要记下它所在的行和列。于是采取如下方法：将非零元素所在的行、列以及它的值构成一个三元组（i,j,v），然后再按某种规律存储这些三元组，这种方法可以节约存储空间。下面讨论稀疏矩阵的压缩存储方法。 将三元组按行优先的顺序，同一行中列号从小到大的规律排列成一个线性表，称为三元组表（有序的双下标法），采用顺序存储方法存储该表。如图5.11稀疏矩阵对应的三元组表为图 5.12。 显然，要唯一的表示一个稀疏矩阵，还需要存储三元组表的同时存储该矩阵的行、列，为了运算方便，矩阵的非零元素的个数也同时存储。这种存储的思想实现如下 define MAXSIZE 12500 /*一个足够大的数*/ typedef struct { int i,j; /*非零元素的行、列*/ elemtype e; /*非零元素值*/ }Triple; /*三元组类型*/ typedef struct { int mu,nu,tu; /*矩阵的行、列及非零元素的个数*/ Triple data[MAXSIZE]; /*三元组表*/ } TSMatrix; /*三元组表的存储类型*/ 这样的存储方法确实节约了存储空间，但矩阵的运算从算法上可能变的复杂些。下面我们讨论这种存储方式下的稀疏矩阵的两种运算：转置和相乘。 稀疏矩阵的转置 设TSMatrix A; 表示一m*n的稀疏矩阵，其转置B则是一个n*m的稀疏矩阵，因此也有 TSMatrix B; 由A求B需要： A的行、列转化成B的列、行； 将A.data中每一三元组的行列交换后转化到B.data中； 看上去以上两点完成之后，似乎完成了B，没有。因为我们前面规定三元组的是按一行一行且每行中的元素是按列号从小到大的规律顺序存放的，因此B也必须按此规律实现，A的转置B如图5.13所示，图5.14是它对应的三元组存储，就是说，在A的三元组存储基础上得到B的三元组表存储（为了运算方便，矩阵的行列都从1算起，三元组表data也从1单元用起）。 转置实现算法思路： ①将矩阵的行列值相互交换； ②将每个三元组中的I和j相互调换； ③重排三元组之间的次序便可实现矩阵的转置。 第三步的实现： A、按照b.data中三元组的次序依次在a.data中找到相应的三元组进行转置。具体算法描述见P98 算法5.1。其仅适于tu《mu×nu的情况。 B、按照a.data中的三元组的次序进行转置，并将转置后的三元组置入b中恰当的位置。为求出每列中b.data的恰当位置，在转置前，应先求得M的每一列中非零元的个数。 需要附设num和cpot两个向量。Num[col]表示矩阵M中第col列中非零元的个数，cpot[col]指示M中第col列的第一个非零元在b.data中的恰当位置。显然有：copt[1]=1 copt[col]= copt[col-1]+ Num[col-1] 2≤col≤a.nu 算法见P100 算法5.2 将这种转置称为快速转置。 总结：非零元在表中按行序有序存储，因此便于进行以行顺序处理的矩阵运算，然而若序按行号存取某一行的非零元，则需从头开始进行查找。 2．稀疏矩阵的乘积 为了实现矩阵相乘，应便于随机存取任意一行的非零元，则需知道每一行的第一个非零元在三元组表中的位置，为此，将快速转置矩阵算法中创建的，指示“行”信息的辅助数组copt固定在稀疏矩阵的存储结构中，称这种“带行链接信息”的三元组表称为逻辑连接的顺序表，其类型描述如下： Typedef struct{ Triple data[MAXSIZE+1]; Int rpos[MAXRC+1]; Int mu,nu,tu; }RLSMatrix; 下面讨论两个矩阵相乘已知稀疏矩阵A(m1× n1)和B(m2× n2)，求乘积C(m1× n2)。 稀疏矩阵A、B、C及它们对应的三元组表A.data、B.data、C.data如图5.16所示。 由矩阵乘法规则知： C（i，j）=A(i,1)×B(1,j)+A(i,2)×B(2,j)+…+A(i,n)×B(n,j) = 这就是说只有A(i,k)与B(k,p)（即A元素的列与B元素的行相等的两项）才有相乘的机会，且当两项都不为零时，乘积中的这一项才不为零。 矩阵用二维数组表示时，传统的矩阵乘法是：A的第一行与B的第一列对应相乘累加后得到c11，A的第一行再与B的第二列对应相乘累加后得到c12，…，因为现在按三元组表存储，三元组表是按行为主序存储的，在B.data中，同一行的非零元素其三元组是相邻存放的，同一列的非零元素其三元组并未相邻存放，因此在B.data中反复搜索某一列的元素是很费时的，因此改变一下求值的顺序，以求c11和c12为例，因为 即a11只有可能和B中第1行的非零元素相乘，a12只有可能和B中第2行的非零元素相乘，…，而同一行的非零元是相邻存放的，所以求c11和c12同时进行：求a11*b11累加到c11，求a11*b12累加到c12，再求a12*b21累加到c11，再求a12*b22累加到c22.，…，当然只有aik和bkj(列号与行号相等)且均不为零（三元组c11= c12= 解释 a11*b11+ a11*b12+ a11只与B中1行元素相乘 a12*b21+ a12*b22+ a12只与B中2行元素相乘 a13*b31+ a13*b32+ a13只与B中3行元素相乘 a14*b41 a14*b42 a14只与B中4行元素相乘 存在）时才相乘，并且累加到cij当中去。 为了运算方便，设一个累加器：datatype temp[n+1]；用来存放当前行中cij的值，当前行中所有元素全部算出之后，再存放到C.data中去。 为了便于B.data中寻找B中的第k行第一个非零元素，与前面类似，在此需引入num和rpot两个向量。num[k]表示矩阵B中第k行的非零元素的个数；rpot[k]表示第k行的第一个非零元素在B.data中的位置。 于是有 rpot[1]=1 rpot[k]=rpot[k-1]+num[k-1] 2≤k≤n 例如 ，对于矩阵B的num和rpot如图5.17所示。 col 1 2 3 4 num[col] 2 0 2 0 cpot[col] 1 3 3 5 图5.17 矩阵B的num与rpot值 根据以上分析，稀疏矩阵的乘法运算的粗略步骤如下： ⑴初始化。清理一些单元，准备按行顺序存放乘积矩阵; ⑵求B的num，rpot; ⑶做矩阵乘法。将A.data中三元组的列值与B.data中三元组的行值相等的非零元素相乘，并将具有相同下标的乘积元素相加。见P102 算法5.3 分析上述算法的时间性能如下：（1）求num的时间复杂度为O(B->nu+B->tu)；（2）求rpot 时间复杂度为O(B->mu)；（3）求temp时间复杂度为O(A->mu*B->nu)；（4）求C的所有非零元素的时间复杂度为O(A->tu*B->tu/B->mu)；（5）压缩存储时间复杂度为O(A->mu*B->nu)；所以总的时间复杂度为O(A->mu*B->nu+(A->tu*B->tu)/B->nu)。 四、稀疏矩阵的十字链表存储 三元组表可以看作稀疏矩阵顺序存储，但是在做一些操作（如加法、乘法）时，非零项数目及非零元素的位置会发生变化，这时这种表示就十分不便。在这节中，我们介绍稀疏矩阵的一种链式存储结构——十字链表，它同样具备链式存储的特点，因此，在某些情况下，采用十字链表表示稀疏矩阵是很方便的。 图5.18是一个稀疏矩阵的十字链表。 1、▲：用十字链表表示稀疏矩阵的基本思想是：对每个非零元素存储为一个结点，结点由5个域组成，其结构如图5.19表示，其中：row域存储非零元素的行号，col域存储非零元素的列号， 域存储本元素的非零元值，right，down是两个指针域。 图5.19 十字链表的结点结构 图5.19 十字链表的结点结构 稀疏矩阵中每一行的非零元素结点按其列号从小到大顺序由right域链成一个带表头结点的循环行链表，同样每一列中的非零元素按其行号从小到大顺序由down域也链成一个带表头结点的循环列链表。即每个非零元素aij既是第i行循环链表中的一个结点，又是第j列循环链表中的一个结点。整个矩阵构成一个十字交叉的链表，顾称这样的存储结构为十字链表。可用两个存储行链表的头指针和列链表的头指针的一维数组之。见P104 图5.6 2、稀疏矩阵的十字链表存储表示：见P104 算法5.4 3、矩阵运算： ①．建立稀疏矩阵A的十字链表 首先输入的信息是：m（A的行数），n（A的列数），r（非零项的数目），紧跟着输入的是r个形如（i,j,aij）的三元组。 算法的设计思想是：首先建立每行（每列）只有头结点的空链表，并建立起这些头结点拉成的循环链表；然后每输入一个三元组（i，j，aij），则将其结点按其列号的大小插入到第i个行链表中去，同时也按其行号的大小将该结点插入到第j个列链表中去。在算法中将利用一个辅助数组MNode *hd[s+1]; 其中 s=max(m , n) , hd [i]指向第i行(第i列)链表的头结点。这样做可以在建立链表时随机的访问任何一行（列），为建表带来方便。 见P 104 算法5.4 建立稀疏矩阵的十字链表 上述算法中，建立头结点循环链表时间复杂度为O(S)，插入每个结点到相应的行表和列表的时间复杂度是O(t*S)，这是因为每个结点插入时都要在链表中寻找插入位置，所以总的时间复杂度为O(t*S)。该算法对三元组的输入顺序没有要求。如果我们输入三元组时是按以行为主序（或列）输入的，则每次将新结点插入到链表的尾部的，改进算法后，时间复杂度为O(S+t)。 ②．两个十字链表表示的稀疏矩阵的加法 v已知两个稀疏矩阵A和B，分别采用十字链表存储，计算C=A+B，C也采用十字链表方式存储，并且在A的基础上形成C。 由矩阵的加法规则知，只有A和B行列对应相等，二者才能相加。C中的非零元素cij只可能有3种情况：或者是aij+bij，或者是aij　(bij=0)，或者是bij　(aij=0)，因此当B加到A上时，对A十字链表的当前结点来说，对应下列四种情况：或者改变结点的值（aij+bij≠0），或者不变（bij＝0），或者插入一个新结点（aij＝0），还可能是删除一个结点（aij+bij＝0）。整个运算从矩阵的第一行起逐行进行。对每一行都从行表的头结点出发，分别找到A和B在该行中的第一个非零元素结点后开始比较，然后按4种不同情况分别处理。 设pa和pb分别指向A和B的十字链表中行号相同的两个结点，4种情况如下： (1) 若pa->col=pb->col且pa->v+pb->v≠0，则只要用aij+bij的值改写pa所指结点的值域即可。 (2) 若pa->col=pb->col且pa->v+pb->v=0，则需要在矩阵A的十字链表中删除pa所指结点，此时需改变该行链表中前趋结点的right域，以及该列链表中前趋结点的down域。 (3) 若pa->col < pb->col且pa->col≠0（即不是表头结点），则只需要将pa指针向右推进一步，并继续进行比较。 (4) 若pa->col > pb->col或pa->col＝0（即是表头结点），则需要在矩阵A的十字链表中插入一个pb所指结点。 为了便于运算，在A的行链表上设pre指针，指示pa所指结点的前驱结点；在A的每一列的链表上设一个指针hl[j]，它的初值和列链表的头指针相同，即hl[j]=chead[j]。矩阵相加操作的概要描述，见P105。