稀疏矩阵

概念

• 对于一个矩阵，我们非常自然的是将其存储在一个二维数组$a[m][n]$中，但对于一个矩阵，它的很多元素都为0，这样的矩阵我们叫做 “稀疏矩阵” ；
• 比如一个$1000\times 1000$的矩阵，它只有$1000$个非零元素，如果我们用$1000\times 1000$的数组存储的话，需要$1000000$个存储空间，同时在进行转置、乘法、加法时会花费大量时间，所以我们将进一步思考，如何更好的表示稀疏矩阵；

稀疏矩阵的表示

类似于”稀疏多项式“的表示方法（详细内容见：数据结构-数组-多项式），将矩阵中每个非零元素唯一表示为$<row,col,value>$形式的三元组（假设$value$为整数），非零元素按照由行到列的顺序排列，同时加上转置、乘法、加法等必要的操作，ADT定义如下：

//稀疏矩阵
class SparseMatrix;

//非零元素数据
class MatrixTerm
{
	friend class SparseMatrix;
private:
	int row, col, value;
};

class SparseMatrix
{
private:
	//行数、列数、非零数
	int rows, cols, terms;
	//非零元素数组
	MatrixTerm* smArray;

public:
	//构造函数1
	SparseMatrix(int r, int c, int t, int rs[],int cs[],int vs[]);

	//构造函数2
	SparseMatrix(int r, int c, int t);

	//默认构造函数
	SparseMatrix();

	//转置矩阵
	SparseMatrix Transpose();

	//快速转置
	SparseMatrix FastTranspose();

	//扩展容量
	void NewTerm(const int c, const int r, const int v);

	//矩阵相加
	SparseMatrix Add(SparseMatrix temp);

	//矩阵相乘
	SparseMatrix Multiply(SparseMatrix temp);

	//重载=号，后续传参
	const SparseMatrix& operator=(const SparseMatrix& temp);

	//矩阵打印
	void S_Show();
};

比如，对于下面的矩阵来说：

$rows=6,cols=6,terms=8,smArray数组如下表$：

矩阵转置

一般方法

对于转置一个二维数组$a[rows][cols]$,我们一般用下面的方法：

for(int i=0; i<cols; i++)
	for(int j=0; j<rows; j++)
		b[i][j]=a[j][i];

那我门接下来将要用稀疏矩阵表示法实现矩阵转置：

实现函数

• 在思考前我们要注意smArray数组是 **按行有序**排列，所以虽然转置只是交换行号与列号，但排列顺序发生了变化，如下图：
  
• 所以我们的目标就是将第$i$列的元素储存到第$i$行中，所以实现函数如下：

SparseMatrix SparseMatrix::Transpose()
{
	cout << "转置矩阵：" << endl;
	//含有非零元素
	if (terms > 0)
	{
		SparseMatrix result(cols, rows, terms);
		int bPos = 0;
		for (int i = 0; i < cols; i++)
		{
			for (int j = 0; j < terms; j++)
			{
				//找到smArray中col等于i
				if (smArray[j].col == i)
				{
					result.smArray[bPos].row = i;
					result.smArray[bPos].col = smArray[j].row;
					result.smArray[bPos++].value = smArray[j].value;
				}
			}
		}
		return result;
	}
	else return *this;
}

函数分析

• 从算法中看出，该函数的运行时间为$O(cols\times terms)$；
• 已知用二维数组表示是转置矩阵运行时间为$O(rows\times cols)$，但用稀疏矩阵表示时，当$terms$达到$rows\times cols$数量级时，此时转置的运行时间达到了$O(rows\times col{s}^2)$，可见我们在节省内存空间时，可能浪费了大量的运行时间；
• 所以我们可以适当的使用些许内存获取$result.smArray$的起始位置，将$\ast this$中元素逐个移到$result$的正确位置上，接下来介绍这种实现方法，“快速转置”（FastTranspose）。

快速转置(FastTranspose)

实现函数

• 首先我们先定义两个数组：
  • $\ast rowSize$：$rowSize[i]$表示 $i$列非零元素的数目 ；
  • $\ast rowStart$：$rowStart[j]$表示$this->smArray[i]$在$result.smArray[j]$位置上，$j=this->smArray[i].col$
    • *即rowStart[i]代表result中第i行（this中的第i）result.smArray[]的开始位置；
    • 可以推出$\ast rowStart$的计算公式为：
      $$rowStart[i]=\{\begin{array}{l}0,i=0\\ rowStart[i-1]+rowSize[i-1],1\le i\le cols\end{array}$$
• 比如下面的矩阵：

$rows=6,cols=6,terms=8,smArray数组如下表$：

两数组的取值如下表：

• 通过定义的两个数组确定函数定义，代码如下：

SparseMatrix SparseMatrix::FastTranspose()
{
	cout << "进行快转置：" << endl;

	if (terms > 0)
	{
		SparseMatrix result(cols, rows, terms);
		int* rowSize = new int[cols];
		int* rowStart = new int[cols];
		fill(rowSize, rowSize + cols, 0);//rowSize全部填充为0

		//统计*this每列(result每行)非0的数目
		for (int i = 0; i < terms; i++)
			rowSize[smArray[i].col]++;

		//rowStart[i]代表result中第i行result.smArray[]的开始位置
		rowStart[0] = 0;
		for (int i = 1; i < cols; i++)
			rowStart[i] = rowStart[i - 1] + rowSize[i - 1];

		for (int i = 0; i < terms; i++)
		{
			int j = rowStart[smArray[i].col];
			result.smArray[j].row = smArray[i].col;
			result.smArray[j].col = smArray[i].row;
			result.smArray[j].value = smArray[i].value;
			rowStart[smArray[i].col]++;
		}

		//删掉缓存
		delete[]rowSize;
		delete[]rowStart;

		return result;
	}

	else return *this;
}

函数分析：

• 通过分析算法可以得到整体算法的时间复杂度为$O(cols+terms)$；
• 当$terms$达到$rows\times cols$数量级时，$O(cols+terms)$变成了$O(cols\times rows)$，与二维数组一样，但算法的常数因子大于二维数组算法；
• 当$terms$比$rows\times cols$小的多时，算法既节省了空间和时间。

稀疏矩阵乘法

• 要实现矩阵a乘以矩阵b得到矩阵c前，要可以对矩阵实现smArray容量拓展，在不改变之前的数据，同时添加新的数据，代码如下：

void SparseMatrix::NewTerm(const int r, const int c, const int v)
{
	MatrixTerm* temp = new MatrixTerm[terms + 1];
	if (terms > 0)
	{
		//将旧smArray拷贝到temp中
		copy(smArray, smArray + terms, temp);

		//删除旧smArray
		delete[]smArray;
	}
	//smArray指针指向temp
	smArray = temp;

	smArray[terms].row = r;
	smArray[terms].col = c;
	smArray[terms++].value = v;
}

函数实现

• 因为$smArray$是按行排序，为了方便后矩阵的遍历 ，要将后矩阵进行转置，这样就可以在前矩阵的每一行遍历后矩阵的每一行（转置后），这样就可以实现矩阵乘法，具体代码如下：

SparseMatrix SparseMatrix::Multiply(SparseMatrix temp)
{
	cout << "进行矩阵乘法：" << endl;

	if (cols != temp.rows)throw"前矩阵列数不等于后矩阵行数，无法相乘！";

	SparseMatrix result(rows, temp.cols, 0);
	//将后矩阵进行转置，方便循环遍历
	SparseMatrix temp_t = temp.FastTranspose();

	//aPos：this->smArray索引；
	//aPosNext：this->smArray下一行索引；
	//rowIdex：result目前行数索引
	int aPos = 0, aPosNext = 0;
	int rowIdex = smArray[0].row;
	int sum = 0;

	while (aPos < terms)
	{
		//bPos：temp_t.smArray索引；
		///colIdex:result目前列数索引
		int bPos = 0;
		int colIdex = temp_t.smArray[0].row;

		while (bPos < temp_t.terms)
		{
			if (smArray[aPos].row != rowIdex)
			{
				result.NewTerm(rowIdex, colIdex, sum);
				sum = 0;
				aPos = aPosNext;
				//矩阵temp下一列
				while (temp_t.smArray[bPos].row == colIdex && bPos < temp_t.terms)
					bPos++;
				colIdex = temp_t.smArray[bPos].row;
			}
			else if (temp_t.smArray[bPos].row != colIdex)
			{
				result.NewTerm(rowIdex, colIdex, sum);
				sum = 0;
				aPos = aPosNext;
				//矩阵temp下一列
				colIdex = temp_t.smArray[bPos].row;
			}
			else
			{
				if (smArray[aPos].col < temp_t.smArray[bPos].col)
					aPos++;
				else if (smArray[aPos].col == temp_t.smArray[bPos].col)
				{
					sum += smArray[aPos].value * temp_t.smArray[bPos].value;
					aPos++;
					bPos++;
				}
				else
					bPos++;
			}
		}
		//矩阵this下一行
		while (smArray[aPos].row == rowIdex && aPos < terms)
			aPos++;
		aPosNext = aPos;
		rowIdex = smArray[aPos].row;
	}

	return result;
}

函数分析

• 通过分析算法，循环总耗时为$O({\sum }_r(b.cols\times t_r+b.terms))$，其中$t_r$为第$r$行的非零元素数，极端情况下，循环的总耗时$O({\sum }_r(b.cols\times t_r+b.terms))=O(b.cols\times a.terms+a.rows\times b.terms)$；
• 与数组表示矩阵进行比较，算法为：

for(int i = 0; i<a.rows; i++)
	for(int j=0; j<b.cols; j++)
	{
		sum=0;
		for(int k = 0; k<a.cols; k++)
			sum+=a[i][j] * b[k][j];
		c[i][j] = sum;
	{

• 该算法时间复杂度为$O(a.rows\times a.cols\times b.cols)$，因为$a.terms\le a.rows\times b.cols$，并且$b.terms\le a.cols\times b.rows$，所以还是稀疏矩阵表示时间开销更小。最坏情况下，算法较慢，但terms远小于最大值时，算法更加优越。

稀疏矩阵加法

• 加法并不是很复杂，循环遍历就好，代码如下：

SparseMatrix SparseMatrix::Add(SparseMatrix temp)
{
	cout << "稀疏矩阵相加：" << endl;
	if (rows != temp.rows || cols != temp.cols)throw"两矩阵行列数不相等，不能相加！！！";

	int idxa = 0, idxb = 0;
	SparseMatrix result(rows, cols, 0);

	while (idxa < terms && idxb < temp.terms)
	{
		if (smArray[idxa].row < temp.smArray[idxb].row)
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
			idxa++;
		}
		else if (smArray[idxa].row > temp.smArray[idxb].row)
		{
			result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);
			idxb++;
		}
		else if (smArray[idxa].col < temp.smArray[idxb].col)
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
			idxa++;
		}
		else if (smArray[idxa].col > temp.smArray[idxb].col)
		{
			result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);
			idxb++;
		}
		else
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value + temp.smArray[idxb].value);
			idxa++;
			idxb++;
		}
	}
	//未完全遍历
	for (; idxa < terms; idxa++)
		result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
	for (; idxb < temp.terms; idxb++)
		result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);

	return result;
}

多维数组

• 一种常用的表示方法：行主序法，即将多维数组展开为一维数组，这里不做详细介绍。

代码总览

• 可自己用数据测试；

#include <iostream>
using namespace std;

//稀疏矩阵
class SparseMatrix;

//非零元素数据
class MatrixTerm
{
	friend class SparseMatrix;
private:
	int row, col, value;
};

class SparseMatrix
{
private:
	//行数、列数、非零数
	int rows, cols, terms;
	MatrixTerm* smArray;

public:
	//构造函数1
	SparseMatrix(int r, int c, int t, int rs[], int cs[], int vs[]);

	//构造函数2
	SparseMatrix(int r, int c, int t);

	//默认构造函数
	SparseMatrix();

	//转置矩阵
	SparseMatrix Transpose();

	//快速转置
	SparseMatrix FastTranspose();

	//扩展容量
	void NewTerm(const int c, const int r, const int v);

	//矩阵相加
	SparseMatrix Add(SparseMatrix temp);

	//矩阵相乘
	SparseMatrix Multiply(SparseMatrix temp);

	//重载=号，后续传参
	const SparseMatrix& operator=(const SparseMatrix& temp);

	//矩阵打印
	void S_Show();
};

SparseMatrix::SparseMatrix(int r, int c, int t, int rs[], int cs[], int vs[])
{
	rows = r;
	cols = c;
	terms = t;
	smArray = new MatrixTerm[terms];

	for (int i = 0; i < terms; i++)
	{
		smArray[i].row = rs[i];
		smArray[i].col = cs[i];
		smArray[i].value = vs[i];
	}
}

SparseMatrix::SparseMatrix(int r, int c, int t)
{
	rows = r;
	cols = c;
	terms = t;
	if (terms > 0)
		smArray = new MatrixTerm[terms];
	else smArray = NULL;
}

SparseMatrix::SparseMatrix()
{
	rows = 0;
	cols = 0;
	terms = 0;
	smArray = NULL;
}

SparseMatrix SparseMatrix::Transpose()
{
	cout << "转置矩阵：" << endl;

	if (terms > 0)
	{
		SparseMatrix result(cols, rows, terms);
		int bPos = 0;

		for (int i = 0; i < cols; i++)
		{
			for (int j = 0; j < terms; j++)
			{
				//找到smArray中col等于i
				if (smArray[j].col == i)
				{
					result.smArray[bPos].row = i;
					result.smArray[bPos].col = smArray[j].row;
					result.smArray[bPos++].value = smArray[j].value;
				}
			}
		}

		return result;
	}
	else return *this;
}

SparseMatrix SparseMatrix::FastTranspose()
{
	cout << "进行快转置：" << endl;

	if (terms > 0)
	{
		SparseMatrix result(cols, rows, terms);
		int* rowSize = new int[cols];
		int* rowStart = new int[cols];
		fill(rowSize, rowSize + cols, 0);//rowSize全部填充为0

		//统计*this每列(result每行)非0的数目
		for (int i = 0; i < terms; i++)
			rowSize[smArray[i].col]++;

		//rowStart[i]代表result中第i行result.smArray[]的开始位置
		rowStart[0] = 0;
		for (int i = 1; i < cols; i++)
			rowStart[i] = rowStart[i - 1] + rowSize[i - 1];

		for (int i = 0; i < terms; i++)
		{
			int j = rowStart[smArray[i].col];
			result.smArray[j].row = smArray[i].col;
			result.smArray[j].col = smArray[i].row;
			result.smArray[j].value = smArray[i].value;
			rowStart[smArray[i].col]++;
		}

		//删掉缓存
		delete[]rowSize;
		delete[]rowStart;

		return result;
	}

	else return *this;
}

void SparseMatrix::NewTerm(const int r, const int c, const int v)
{
	MatrixTerm* temp = new MatrixTerm[terms + 1];
	if (terms > 0)
	{
		//将旧smArray拷贝到temp中
		copy(smArray, smArray + terms, temp);

		//删除旧smArray
		delete[]smArray;
	}
	//smArray指针指向temp
	smArray = temp;

	smArray[terms].row = r;
	smArray[terms].col = c;
	smArray[terms++].value = v;
}

SparseMatrix SparseMatrix::Add(SparseMatrix temp)
{
	cout << "稀疏矩阵相加：" << endl;
	if (rows != temp.rows || cols != temp.cols)throw"两矩阵行列数不相等，不能相加！！！";

	int idxa = 0, idxb = 0;
	SparseMatrix result(rows, cols, 0);

	while (idxa < terms && idxb < temp.terms)
	{
		if (smArray[idxa].row < temp.smArray[idxb].row)
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
			idxa++;
		}
		else if (smArray[idxa].row > temp.smArray[idxb].row)
		{
			result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);
			idxb++;
		}
		else if (smArray[idxa].col < temp.smArray[idxb].col)
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
			idxa++;
		}
		else if (smArray[idxa].col > temp.smArray[idxb].col)
		{
			result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);
			idxb++;
		}
		else
		{
			result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value + temp.smArray[idxb].value);
			idxa++;
			idxb++;
		}
	}
	//未完全遍历
	for (; idxa < terms; idxa++)
		result.NewTerm(smArray[idxa].row, smArray[idxa].col, smArray[idxa].value);
	for (; idxb < temp.terms; idxb++)
		result.NewTerm(temp.smArray[idxb].row, temp.smArray[idxb].col, temp.smArray[idxb].value);

	return result;
}

SparseMatrix SparseMatrix::Multiply(SparseMatrix temp)
{
	cout << "进行矩阵乘法：" << endl;

	if (cols != temp.rows)throw"前矩阵列数不等于后矩阵行数，无法相乘！";

	SparseMatrix result(rows, temp.cols, 0);
	//将后矩阵进行转置，方便循环遍历
	SparseMatrix temp_t = temp.FastTranspose();

	//aPos：this->smArray索引；
	//aPosNext：this->smArray下一行索引；
	//rowIdex：result目前行数索引
	int aPos = 0, aPosNext = 0;
	int rowIdex = smArray[0].row;
	int sum = 0;

	while (aPos < terms)
	{
		//bPos：temp_t.smArray索引；
		///colIdex:result目前列数索引
		int bPos = 0;
		int colIdex = temp_t.smArray[0].row;

		while (bPos < temp_t.terms)
		{
			if (smArray[aPos].row != rowIdex)
			{
				result.NewTerm(rowIdex, colIdex, sum);
				sum = 0;
				aPos = aPosNext;
				//矩阵temp一列
				while (temp_t.smArray[bPos].row == colIdex && bPos < temp_t.terms)
					bPos++;
				colIdex = temp_t.smArray[bPos].row;
			}
			else if (temp_t.smArray[bPos].row != colIdex)
			{
				result.NewTerm(rowIdex, colIdex, sum);
				sum = 0;
				aPos = aPosNext;
				//矩阵temp下一列
				colIdex = temp_t.smArray[bPos].row;
			}
			else
			{
				if (smArray[aPos].col < temp_t.smArray[bPos].col)
					aPos++;
				else if (smArray[aPos].col == temp_t.smArray[bPos].col)
				{
					sum += smArray[aPos].value * temp_t.smArray[bPos].value;
					aPos++;
					bPos++;
				}
				else
					bPos++;
			}
		}
		//矩阵this下一行
		while (smArray[aPos].row == rowIdex && aPos < terms)
			aPos++;
		aPosNext = aPos;
		rowIdex = smArray[aPos].row;
	}

	return result;
}

const SparseMatrix& SparseMatrix::operator=(const SparseMatrix& temp)
{
	cols = temp.cols;
	rows = temp.rows;
	terms = temp.terms;

	if (smArray != NULL)
		delete[]smArray;

	if (terms > 0)
	{
		smArray = new MatrixTerm[terms];
		copy(temp.smArray, temp.smArray + terms, smArray);
	}
	else smArray = NULL;

	return *this;
}

void SparseMatrix::S_Show()
{
	if (rows > 0 && cols > 0)
	{
		for (int i = 0; i < rows; i++)
		{
			cout << "[";
			for (int j = 0; j < cols; j++)
			{
				int out = 0;
				for (int t = 0; t < terms; t++)
				{
					if (smArray[t].row == i && smArray[t].col == j)
					{
						out = smArray[t].value;
					}
				}
				cout << out << "\t";
			}
			cout << "]\n";
		}
	}
	else cout << "[ ]\n";
	cout << "稀疏矩阵信息：\n行数：" << rows << "\t列数：" << cols << "\t非零项数：" << terms << endl << endl;
}

int main()
{
	int rs1[7] = { 0,0,1,1,3,5,7 }, rs2[10] = { 0,0,0,1,3,4,6,6,6,7 };
	int cs1[7] = { 1,2,4,7,5,3,4 }, cs2[10] = { 0,3,8,4,2,5,2,4,7,6 };
	int vs1[7] = { 20,11,23,45,13,6,16 }, vs2[10] = { 23,12,55,13,56,13,15,46,24,66 };

	SparseMatrix a(8, 9, 7, rs1, cs1, vs1), b(8, 9, 10, rs2, cs2, vs2), c;
	a.S_Show();
	b.S_Show();

	c = a.Transpose();
	c.S_Show();

	c = a.FastTranspose();
	c.S_Show();

	c = a.Multiply(c);
	c.S_Show();

	c = a.Add(b);
	c.S_Show();

	return 0;
}

上一篇：数据结构-数组-多项式