稀疏矩阵的存储和乘法操作

一 稀疏矩阵的存储

1.三元组顺序表

三元组表示法就是在存储非零元的同时，存储该元素所对应的行下标和列下标。稀疏矩阵中的每一个非零元素由一个三元组(i,j,a_ij)唯一确定。矩阵中所有非零元素存放在由三元组组成的顺序表中(通常用数组)。所以三元组的逻辑结构如下：

//————稀疏矩阵的三元组表示法————//
#define MAX_SIZE 1500  //表示稀疏矩阵的非零元素的最大个数
class Triple
{
    int i,j;//表示非零元素的行下表和列下标
    int val;//非零元素的值，此处以int类型为例
};
class TSMatrix
{
    Triple data[MAX_SIZE];
    int row_num,col_num,cnt;//稀疏矩阵的行数、列数以及非零元素的个数
};

注意，此处的非零元素的三元组是以行序为主序顺序排列的。

2.行逻辑链接顺序表

行逻辑链接顺序表的实质就是在三元组顺序表的基础上加了一个数组，这个数组用于存储稀疏矩阵中每行的第一个非零元素的在三元组顺序表中的位置(此处一定要理解对，是在三元组顺序表中的位置)。所以其逻辑结构如下：

//————稀疏矩阵的行逻辑链接表示法————//
#define MAX_SIZE 1500  //表示稀疏矩阵的非零元素的最大个数
#define MAX_ROW 1500  //表示稀疏矩阵的行数的最大个数
class Triple
{
    int i,j;//表示非零元素的行下表和列下标
    int val;//非零元素的值，此处以int类型为例
};
class RLSMatrix
{
    Triple data[MAX_SIZE]; //非零元三元组表
    int rpos[MAX_ROW];//每行第一个非零元素的位置
    int row_num,col_num,cnt;//稀疏矩阵的行数、列数以及非零元素的个数
};

 3.十字链表

当稀疏矩阵的非零元个数和位置在操作过程中变化较大时，就不易采用顺序存储结构来表示三元组的线性表。对于这种类型得矩阵，采用链式存储结构表示三元组的线性表更为恰当。

在链表中，每个非零元可用一个含5个域的结点表示，其中$i$、$j$、$val$这三个域分别表示该非零元所在的行、列和非零元的值(就是三元组中的那三个域)，向右域right用以链接同一行中下一个非零元，向下域down用以链接同一行中的下一个非零元。

所以，同一行的非零元通过right域链接成一个线性链表，同一列的非零元通过down域链接成一个线性表，每个非零元既是某个行链表中的一个结点，又是某个列链表中的一个结点，整个矩阵构成了一个十字交叉的链表，所以才称为十字链表。十字链表中，可以使用两个分别存储行链表的头指针和列链表的头指针的一维数组表示。所以，其逻辑结构如下：

//————稀疏矩阵的十字链表表示法————//
class CrossNode
{
    int i,j;//表示非零元素的行下表和列下标
    int val;//非零元素的值，此处以int类型为例
    CrossNode *right,*down; //非零元所在行表和列表的后继指针域
};
class CrossList
{
    CrossNode *row_head,*col_head; //行链表和列链表的头指针
    int row_num,col_num,cnt;//稀疏矩阵的行数、列数以及非零元素的个数
};

二 稀疏矩阵的乘法　

1.一般矩阵的乘法操作　

对于一般矩阵而言，通常用二维数组表示矩阵。对于矩阵$M∈R^{m imes n}$和矩阵$N∈R^{n imes p}$而言，它们的乘法操作表示如下：

#include<iostream>
#include<unordered_map>
#include<queue>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<sstream>
#include<set>
#include<map>
using namespace std;
#define MAX_NUM 2
void MatrixMultiply(int M[][MAX_NUM],int N[][MAX_NUM],int rs[][MAX_NUM],int m,int n,int p)
{
    for(int i = 0; i < m; i++)
    {
        for(int j = 0; j < p; j++)
        {
            rs[i][j] = 0;//初始化第行第j列的元素
            for(int k = 0;k < n;k++)
            {
                rs[i][j] += (M[i][k] * N[k][j]);//第i行与第j列相乘求和
            }
        }
    }
}

int main()
{
    
    int M[][2] = {{1,1},{2,3}};
    int N[][2] = {{2,2},{4,5}};
    int rs[][2] ={{0,0},{0,0}};
    MatrixMultiply(M,N,rs,2,2,2);
    for(int i = 0 ;i < 2;i++)
    {

        for(int j = 0 ; j < 2; j++)
        cout<<rs[i][j]<<" ";
        cout<<endl;
    }
}

这种最通用的矩阵乘法操作算法的时间复杂度为O(mnp)。在这个算法中，最关键之处在于求解rs[i][j]的值，但是仔细观察可以知道，上述算法中不论M[i][k]和N[k][j]是否为零，都要进行一次乘法运算，而实际上，这两者有一个值为零，其乘积必然为零。所以，在进行矩阵乘法时，应免去这种操作(如果用这种方法来进行稀疏矩阵的乘法运算的话，更应该免去这种无效操作)，相当于在第三个循环中加入判断条件条件：

if(M[i][k] != 0 && N[k][j] != 0)

　　

2.稀疏矩阵的乘法操作　

对于稀疏矩阵而言，我们可以通过上述1中的通用的矩阵乘法操作进行运算，但是这样往往会有很多不必要的运算。所以我们可以进一步优化，最常见的是用行逻辑链接的顺序表的表示形式进行乘法操作，本文也只讲解这种表示形式的矩阵乘法的操作。

实现代码如下(参考严蔚敏《数据结构》P102)：

#define MAX_SIZE 1500  //表示稀疏矩阵的非零元素的最大个数
#define MAX_ROW 1500  //表示稀疏矩阵的行数的最大个数
class Triple
{
    public:
    int i,j;//表示非零元素的行下表和列下标
    int val;//非零元素的值，此处以int类型为例
};
class RLSMatrix
{
    public:
    Triple data[MAX_SIZE]; //非零元三元组表
    int rpos[MAX_ROW];//每行第一个非零元素的位置
    int row_num,col_num,cnt;//稀疏矩阵的行数、列数以及非零元素的个数
};

void MultRLSMatrix(RLSMatrix M,RLSMatrix N,RLSMatrix &rs){
	int arow,brow,p,q,ccol,ctemp[MAX_ROW + 1],t,tp;
	if(M.col_num != N.row_num){//不能相乘
		return;
	}
	if(0 == M.cnt * N.cnt ){//有一个是零矩阵
		return;
	}
	//rs初始化
	rs.row_num = M.row_num;
	rs.col_num = N.col_num;
	rs.cnt = 0;
	//从M的第一行开始到最后一行，arow是M的当前行
	for(arow = 1;arow <= M.row_num;arow++){
		for(ccol=1;ccol <= rs.col_num;ccol++){
			ctemp[ccol] = 0;//rs的当前行的各列元素清零
		}
		rs.rpos[arow] = rs.cnt + 1;//开始时从第一个存储位置开始存，后面是基于前面的
		if(arow < M.row_num){
			tp = M.rpos[arow+1];//下一行的起始位置
		}else{
			tp = M.cnt + 1;//最后一行的边界
		}
		for(p = M.rpos[arow];p < tp;p++){
			//对M当前行的每一个非零元
			//找到对应元素在N中的行号，即M中当前元的列号
			brow = M.data[p].j;
			//原理同上
			if(brow < N.row_num){
				t = N.rpos[brow + 1];
			}else{
				t = N.cnt + 1;
			}
			for(q = N.rpos[brow];q < t;q++){
				ccol = N.data[q].j;//乘积元素在rs中列的位置
				ctemp[ccol] += M.data[p].val * N.data[q].val;
			}//for_q
		}//for_p
		//该压缩存储该行非零元了
		for(ccol = 1;ccol <= rs.col_num;ccol++){
			if(0 != ctemp[ccol]){
				if(++rs.cnt > MAX_SIZE){//注意这里有个++
					return;
				}
				rs.data[rs.cnt].i = arow;
				rs.data[rs.cnt].j = ccol;
				rs.data[rs.cnt].val = ctemp[ccol];
			}
		}
	}//for_arow
}

　算法的主要思想如下：

　实际上它的求解过程与1中的算法本质是一样的，只不过充分利用了三元组的性质和rpos数组的信息。为了得到非零元的乘积，只要对M.data[1,2...]中的每个元素即$(i,k,M[i]k])$，找到N.data中所有相应的元素$(k,j,N[k][j])$相乘即可，为此

只需在N.data中找到矩阵N中第k行的所有非零元；而rpos数组正好提供了相关的信息。rpos[row]表示矩阵N中第row行，第一个非零元在N.data中的序号；而rpos[row+1]-1表示矩阵N中第row行最后一个非零元素在N.data中的序号(因为data数组中存储的全都是非零元素，而rpos[row+1]是第row+1行中第一个非零元素在data中的位置，所以它的前一个位置一定是第row行中最后一个非零元素在data中的位置)。

此外，有几点需要注意：

(1)由于两个稀疏矩阵相乘的结果矩阵不一定是稀疏矩阵，所以结果矩阵的表示方式不一定要用行逻辑链接顺序表，也可以用二维数组。

(2)该算法的时间复杂度具体分析见《数据结构》p103，总之，矩阵越稀疏，该算法时间复杂度越低，能达到O(mp)数量级。