C#数据结构(4) 稀疏矩阵与稀疏方阵

导言

线性代数是大学理工科学生的必修课。学过线性代数的同学一定对矩阵不陌生，因为线性代数就是一门关于矩阵的学科。

程序设计中有一种储存数据的方式是二维数组，而二维数组本质上就是矩阵。但是，假如我们想要用二维数组去储存一个大规模的矩阵并进行运算的话，会造成很大的资源浪费。举个例子，假如我们想要用二维数组去储存一个20W*20W的单位矩阵，事实上其中只有20W个数是1，其他数字都是0。所以，我们或许可以利用一种“忽略矩阵中的0项”的方式，来实现对矩阵的压缩储存，这种储存方式就叫做稀疏矩阵。对于大部分位置都是0，只有少部分位置有值的矩阵来说，使用稀疏矩阵可以让矩阵的储存密度大大提高。

我们可以使用一个二维单向链表来储存稀疏矩阵——因为链表在插入方面具有极低的时间复杂度，可以保证一个稀疏矩阵在经过各种运算后依然保持行与列间有序的组织。为了能够实现单向链表的删除，我们需要指向每一个有值结点前一个结点的引用。为此，我们还需要在第一行设置一个空行结点，每一行的第一列设置一个空列结点。形象地说，这种组织方式如下图所示：

使用链表的代价是寻址成本的提高，因此二维迭代器的封装对于基于链表的稀疏矩阵也是必不可少的。此外，比起二维数组矩阵，在各项操作方面，稀疏矩阵更难做到高效。因此如何高效地进行稀疏矩阵的各项操作也是本文需要探讨的话题。

本文中需要实现的稀疏矩阵基本操作如下：

增加与修改单个元素
从二维数组复制元素
矩阵复制
矩阵加减法
数乘
矩阵乘法
矩阵转置

同时稀疏矩阵还有派生类稀疏方阵，它实现的基本操作如下：

求行最简矩阵
求行列式的值
求逆矩阵

1.1行索引链表

基于链表的稀疏矩阵的基础是行索引。因为我们定位一个矩阵中的元素，是先定位其行位置，再定位其列位置的。

因此矩阵包含了一个作为基准的行索引链表结点。这个结点包含一个值，代表行号；包含一个横向指针指向一个列索引链表结点（在C#中就是对一个列索引链表结点的引用），也就是该行第一个有非0元素的列；包含一个纵向指针指向下一个有非0元素的行索引链表结点。

矩阵的基础构造与行索引链表结点如下所示：

namespace exp6lib
{public class mat{private class rowNode{}protected class lineNode{public int lineId;public rowNode firstRow;public lineNode nextLine;public lineNode(){firstRow = new rowNode();firstRow.rowId = -1;}}private lineNode baseLine;protected int rowCount;protected int lineCount;}
}

1.2列索引链表

当确定矩阵中元素所在的行时，再通过列索引链表确定元素所在的列，就可以将该元素定位了。因此，列索引链表中包含两个值，一个表示自己的列号，一个表示由行索引和列索引确定的元素的值。此外，还包含一个指向和自己同行的下一个非0元素的引用。

protected class rowNode
{public int rowId;public double value;public rowNode nextRow;
}

1.3二维迭代器

由于二维链表在定位元素时需要先在行索引链表中寻址，再在列索引链表中寻址，其时间复杂度为 $o(n)$ ，因此在实现诸如矩阵遍历等有序操作时必须依赖二维迭代器，否则会造将大量时间浪费在寻址上，其时间复杂度将是无法想象的。为矩阵维护一个包含一个行索引结点引用、一个列索引结点引用的对象，它将通过四个引用，分别指向一行、一列和该行的前一行、该行的前一列，并且通过所引用的结点内部包含的、指向其他结点的引用来实现自身所指元素的移动。

迭代器在初始化或者行复位时需要将指向前一行的引用放置到矩阵结构的“第0行”上，将指向当前行的引用放置到第一行。同理，每移动到下一行或者进行列复位，都要把指向前一列的引用放置到该行的“第0列”上，将指向当前列的引用放置到第一列。

迭代器实现了五个方法，分别是移动到下一列、移动到下一行并指向该行第一列、移动到该行第一列、在当前指向行和前一行之间插入行、在当前指向列和前一列之间增加列。此外，在大类中还实现了一个私有方法reLine，用来将迭代器恢复至整个稀疏矩阵第一行。在大类中声明一个迭代器对象，用来方便省时地访问其中的元素。

protected class iterator
{public rowNode row;public lineNode line;public rowNode prerow;public lineNode preline;public void nextRow(){prerow = row;row = row.nextRow;}public void nextLine(){preline = line;line = line.nextLine;if (line != null){prerow = line.firstRow;row = line.firstRow.nextRow;}}public void reRow(){row = line.firstRow.nextRow;prerow = line.firstRow;}public void addLine(int lineId){lineNode newLine = new lineNode();newLine.lineId = lineId;newLine.nextLine = preline.nextLine;preline.nextLine = newLine;line = newLine;prerow = newLine.firstRow;row = null;}public void addRow(int rowId, double value){rowNode newRow = new rowNode();newRow.rowId = rowId;newRow.nextRow = prerow.nextRow;prerow.nextRow = newRow;newRow.value = value;row = newRow;}
}protected iterator it;protected void reLine()
{it.preline = baseLine;it.line = baseLine.nextLine;if (it.line != null){it.prerow = it.line.firstRow;it.row = it.line.firstRow.nextRow;}else{it.prerow = null;it.row = null;}
}

1.4 单个元素的操作

利用迭代器指向矩阵中的元素，可以单独设置这个元素的内容。而如果要求设置一个不超出矩阵大小，但原本不存在（也就是为0）的位置的元素，就需要在二维链表中插入行或列。这一判断需要依靠迭代器来作出。比如，若迭代器判定其指向位置的前一行小于目标行、而当前行大于目标行，就要在当前位置前插入目标行。同理，如果把一个原本不为0的位置的元素设为0，就要删除当前列，而列为空时，则要删除当前行。

public void set(int line, int row, double value)
{if (line < 0 || row < 0 ||line>=lineCount||row>=rowCount)return;while (true){if (it.line == null){it.addLine(line);break;}if (it.line.lineId == line)break;if (line<it.line.lineId&&line>it.preline.lineId){it.nextLine();it.addLine(line);break;}if (line <= it.preline.lineId){reLine();continue;}it.nextLine();}while (true){if (it.row == null){if (value == 0)break;it.addRow(row, value);break;}if (it.row.rowId == row){it.row.value = value;if (it.row.value == 0){it.row = it.row.nextRow;it.prerow.nextRow = it.row;if (it.line.firstRow.nextRow == null){it.line = it.line.nextLine;it.preline.nextLine = it.line;it.row = it.line.firstRow.nextRow;it.prerow.nextRow = it.line.firstRow;}}break;}if (row < it.row.rowId && row > it.prerow.rowId){if (value == 0)break;it.addRow(row, value);break;}if (row <= it.prerow.rowId){it.reRow();continue;}it.nextRow();}
}

而单个元素的加减和乘除本质也是相同的。记住要分成三种情况考虑需要实现的操作：

对已有元素操作
把为0的元素变成非0
把非0的元素变成0

protected void add(int line, int row, double value)
{if (line < 0 || row < 0 || line >= lineCount || row >= rowCount)return;while (true){if (line <= it.preline.lineId){reLine();continue;}if (it.line == null){it.addLine(line);break;}if (it.line.lineId == line)break;if (line < it.row.rowId && line > it.preline.lineId){it.nextLine();it.addLine(line);break;}it.nextLine();}while (true){if (row <= it.prerow.rowId){it.reRow();continue;}if (it.row == null){if (value == 0)break;it.addRow(row, value);break;}if (it.row.rowId == row){it.row.value += value;if (it.row.value == 0){it.row = it.row.nextRow;it.prerow.nextRow = it.row;if (it.line.firstRow.nextRow == null){it.line = it.line.nextLine;it.preline.nextLine = it.line;it.row = it.line.firstRow.nextRow;it.prerow.nextRow = it.line.firstRow;}}break;}if (row < it.row.rowId && row > it.prerow.rowId){if (value == 0)break;it.addRow(row, value);break;}it.nextRow();}
}protected void times(int line, int row, double value)
{if (line < 0 || row < 0 || line >= lineCount || row >= rowCount)return;while (true){if (line <= it.preline.lineId){reLine();continue;}if (it.line == null)return;if (it.line.lineId == line)break;if (line < it.row.rowId && line > it.preline.lineId)return;it.nextLine();}while (true){if (row <= it.prerow.rowId){it.reRow();continue;}if (it.row == null)return;if (it.row.rowId == row){it.row.value *= value;if (it.row.value == 0){it.row = it.row.nextRow;it.prerow.nextRow = it.row;if (it.line.firstRow.nextRow == null){it.line = it.line.nextLine;it.preline.nextLine = it.line;it.row = it.line.firstRow.nextRow;it.prerow.nextRow = it.line.firstRow;}}return;}if (row < it.row.rowId && row > it.prerow.rowId)return;it.nextRow();}
}

而获取单个元素使用的方法也运用同样的思想，只不过在获取不到元素时返回元素类型的默认值。

public double get(int line, int row)
{if (line < 0 || row < 0 || line >= lineCount || row >= rowCount)return 0;while (true){if (it.line == null)return 0;if (it.line.lineId == line)break;if (line < it.line.lineId && line > it.preline.lineId)return 0;if (line <= it.preline.lineId){reLine();continue;}it.nextLine();}while (true){if (it.row == null)return 0;if (it.row.rowId == row)return it.row.value;if (row < it.row.rowId && row > it.prerow.rowId)return 0;if (row <= it.prerow.rowId){it.reRow();continue;}it.nextRow();}
}

2.1 稀疏矩阵的初始化

默认初始化时，用户需要指定该矩阵的行数和列数。此外，将迭代器置于第一行第一列，虽然此时第一行第一列没有声明，但是要记得我们矩阵的链表结构中存在一个空的“第0行”。将迭代器的preline引用指向该“第0行”，其他引用则设为null。

public mat(int line,int row)
{baseLine = new lineNode();baseLine.lineId = -1;it = new iterator();it.preline = baseLine;it.line = null;it.row = null;it.prerow = null;lineCount = line;rowCount = row;
}

使用其他mat初始化时，稀疏矩阵会调用目标矩阵中的迭代器和自身的迭代器来实现矩阵遍历。遍历只需要不断调用nextRow方法，在指向的列为空时调用nextLine方法，直到指向的行为空为止即可。

在这里强调一下reLine的作用，它是迭代器的复位方法，也可以视为声明一个指向第一行第一列的迭代器。如果不进行reLine的话，假如目标矩阵中的迭代器已经指向了一个很后面的元素，那么利用迭代器nextRow和nextLine方法进行的遍历就是不成立的。

public mat(mat tar)
{baseLine = new lineNode();baseLine.lineId = -1;it = new iterator();it.preline = baseLine;it.line = null;it.row = null;it.prerow = null;tar.reLine();lineCount = tar.lineCount;rowCount = tar.rowCount;while (tar.it.line!=null){while (tar.it.row != null){set(tar.it.line.lineId, tar.it.row.rowId, tar.it.row.value);tar.it.nextRow();}tar.it.nextLine();}}

使用二维数组初始化通过遍历二维数组的方式，将二维数组中的数据一个一个set到矩阵当中。

public mat(int[][] tar)
{baseLine = new lineNode();baseLine.lineId = -1;it = new iterator();it.preline = baseLine;it.line = null;it.row = null;it.prerow = null;for (int i = 0; i < tar.Length; i++)for (int j = 0; j < tar[i].Length; j++)set(i, j, tar[i][j]);
}

2.2 矩阵加减和数乘

矩阵的加减用第一个矩阵初始化一个新矩阵，同时用新矩阵和第二个矩阵的迭代器遍历两个矩阵，再在新矩阵的每个位置调用加方法，加上第二个矩阵对应位置的元素值（或其相反数），最后把新矩阵返回。

static public mat operator+(mat obj,mat tar)
{if (tar.rowCount != obj.rowCount || tar.lineCount != obj.lineCount)return new mat(0, 0);mat result=new mat(obj);tar.reLine();while (tar.it.line != null){while (tar.it.row != null){result.add(tar.it.line.lineId, tar.it.row.rowId, tar.it.row.value);tar.it.nextRow();}tar.it.nextLine();} return result;
}static public mat operator -(mat obj, mat tar)
{if (tar.rowCount != obj.rowCount || tar.lineCount != obj.lineCount)return new mat(0, 0);mat result = new mat(obj);tar.reLine();while (tar.it.line != null){while (tar.it.row != null){result.add(tar.it.line.lineId, tar.it.row.rowId, -tar.it.row.value);tar.it.nextRow();}tar.it.nextLine();}return result;
}

矩阵的数乘用第一个矩阵初始化一个新矩阵，同时用新矩阵的迭代器遍历新矩阵，再在新矩阵的每个位置调用乘方法，乘上目标值，最后把新矩阵返回。

static public mat operator *(mat obj,double value)
{mat result = new mat(obj);result.reLine();while (result.it.line!= null){while (result.it.row != null){result.times(result.it.line.lineId, result.it.row.rowId, value);result.it.nextRow();}result.it.nextLine();}return result;
}static public mat operator *(double value,mat obj)
{mat result = new mat(obj);result.reLine();while (result.it.line != null){while (result.it.row != null){result.times(result.it.line.lineId, result.it.row.rowId, value);result.it.nextRow();}result.it.nextLine();}return result;
}

2.3 矩阵乘法

依然通过迭代器来选取两个矩阵里的元素，通过矩阵乘法的特殊运算顺序来运算，将结果放进新矩阵中。返回的矩阵是通过new创建的新矩阵，因此不会对原来两个矩阵产生影响。

static public mat operator *(mat tar,mat obj)
{mat result = new mat(tar.lineCount,obj.rowCount);tar.reLine();obj.reLine();while (tar.it.line!=null){for (int i = 0; i < obj.rowCount; i++){double sum = 0;while (tar.it.row != null){sum += tar.it.row.value * obj.get(tar.it.row.rowId, i);tar.it.nextRow();}result.set(tar.it.line.lineId, i,sum);tar.it.reRow();}tar.it.nextLine();}return result;
}

2.4 矩阵的转置

矩阵转置创建一个新的矩阵，其行数为原矩阵的列数，列数为原矩阵的行数，并将原矩阵中的元素所在行列对换后放进新矩阵中。

public mat turn()
{mat result = new mat(rowCount,lineCount);reLine();while (it.line!= null){while (it.row!= null){result.set(it.row.rowId, it.line.lineId,it.row.value);it.nextRow();}it.nextLine();}return result;
}

2.5 打印矩阵

public void print()
{reLine();int lineNum = 0;while(it.line!=null){int rowNum = 0;while(it.row!=null){for (; rowNum < it.row.rowId; rowNum++)Console.Write("0\t");rowNum++;Console.Write(it.row.value);Console.Write("\t");it.nextRow();}for(;rowNum<rowCount;rowNum++)Console.Write("0\t");Console.Write("\n");lineNum++;it.nextLine();}for(;lineNum<lineCount;lineNum++){for (int rowNum=0; rowNum < rowCount; rowNum++)Console.Write("0\t");Console.Write("\n");}
}

3.1 稀疏方阵的初始化

基本上就是从矩阵当中继承了初始化方法。注意构造函数是要继承自base(父类构造函数参数)的。以及行和列的数目应该相同。

public class smat : mat
{public smat(int n):base(n,n){}public smat(smat tar) : base(tar){}
}

3.2 稀疏方阵的基本运算

都是对稀疏矩阵基本运算的重写。由于行列数目相同，实际实现会容易得多。

static public smat operator +(smat obj, smat tar)
{if (tar.rowCount != obj.rowCount || tar.lineCount != obj.lineCount)return new smat(0);smat result = new smat(obj);tar.reLine();while (tar.it.line != null){while (tar.it.row != null){result.add(tar.it.line.lineId, tar.it.row.rowId, tar.it.row.value);tar.it.nextRow();}tar.it.nextLine();}return result;
}static public smat operator -(smat obj, smat tar)
{if (tar.rowCount != obj.rowCount || tar.lineCount != obj.lineCount)return new smat(0);smat result = new smat(obj);tar.reLine();while (tar.it.line != null){while (tar.it.row != null){result.add(tar.it.line.lineId, tar.it.row.rowId, -tar.it.row.value);tar.it.nextRow();}tar.it.nextLine();}return result;
}static public smat operator *(smat obj, double value)
{smat result = new smat(obj);result.reLine();while (result.it.line != null){while (result.it.row != null){result.times(result.it.line.lineId, result.it.row.rowId, value);result.it.nextRow();}result.it.nextLine();}return result;
}static public smat operator *(double value, smat obj)
{smat result = new smat(obj);result.reLine();while (result.it.line != null){while (result.it.row != null){result.times(result.it.line.lineId, result.it.row.rowId, value);result.it.nextRow();}result.it.nextLine();}return result;
}static public smat operator *(smat tar, smat obj)
{smat result = new smat(tar.lineCount);tar.reLine();obj.reLine();while (tar.it.line != null){for (int i = 0; i < obj.rowCount; i++){double sum = 0;while (tar.it.row != null){sum += tar.it.row.value * obj.get(tar.it.row.rowId, i);tar.it.nextRow();}result.set(tar.it.line.lineId, i,sum);tar.it.reRow();}tar.it.nextLine();}return result;
}new public smat turn()
{smat result = new smat(rowCount);reLine();while (it.line != null){while (it.row != null){result.set(it.row.rowId, it.line.lineId, it.row.value);it.nextRow();}it.nextLine();}return result;
}

3.3 求行最简矩阵

求行最简矩阵的目的就是把原矩阵通过初等行变换化为相抵的上三角矩阵。其具体做法是：

从原矩阵中寻找第一个非0元素最靠左（行坐标最小）的一行
从原矩阵中寻找第一个非0元素次靠左（行坐标次小）的一行
如果1、2步找出的行的第一个非0元素行坐标一致，说明需要通过高斯消元法相消。将第一行第一个非0元素设为a，第二行第一个非0元素设为b，第二行每个位置减去第一行每个位置对应元素乘以b/a后的结果，就可以把第二行第一个非0元素变成0，然后重新开始第1步
如果1、2步找出的行的第一个非0元素行坐标不一致，则将第1步找出的行从原矩阵中删除，放入新矩阵，然后重新开始第1步
如果找不到行，说明原矩阵中所有行以及被移入新矩阵，循环结束

public smat tri()
{smat tmp = new smat(this);smat result = new smat(lineCount);lineNode maxline = null;lineNode max2line = null;lineNode premaxline = null;lineNode premax2line = null;tmp.reLine();int i = 0;while (tmp.it.line!= null){while (tmp.it.line!= null){if (maxline == null || tmp.it.row.rowId < maxline.firstRow.nextRow.rowId){max2line = maxline;premax2line = premaxline;maxline = tmp.it.line;premaxline = tmp.it.preline;}else if ((max2line == null || tmp.it.row.rowId < max2line.firstRow.nextRow.rowId) && tmp.it.line.lineId != maxline.lineId){max2line = tmp.it.line;premax2line = tmp.it.preline;}tmp.it.nextLine();}if (max2line!=null&&maxline.firstRow.nextRow.rowId == max2line.firstRow.nextRow.rowId){double basis = max2line.firstRow.nextRow.value / maxline.firstRow.nextRow.value;rowNode rowit = maxline.firstRow.nextRow;while (rowit != null){tmp.add(max2line.lineId, rowit.rowId, -(rowit.value * basis));rowit = rowit.nextRow;}}else{premaxline.nextLine = maxline.nextLine;maxline.nextLine = null;if (premax2line == maxline)premax2line = premaxline;rowNode j = maxline.firstRow.nextRow;while(j!=null){result.set(i, j.rowId, j.value);j = j.nextRow;}i++;result.it.nextLine();maxline = null;premaxline = null;max2line = null;premax2line = null;}tmp.reLine();}return result;
}

3.4 求行列式的值

求完行最简矩阵之后，求行列式的值也变成了非常简单的工作。只需把对角线上所有元素相乘即可。

public double det()
{smat result = tri();int i = 0;double sum = 1;result.reLine();while (result.it.line != null){if (result.it.row.rowId != i)return 0;i++;sum *= result.it.row.value;result.it.nextLine();}return sum;
}

3.5 求逆矩阵

逆矩阵同样通过初等行变换法来解出，方法是把一个行列数与原矩阵相同的单位矩阵进行和原矩阵同步的行变换。当把原矩阵变成单位矩阵时，单位矩阵也变成了原矩阵的逆矩阵。

从原矩阵中寻找第一个非0元素最靠左（行坐标最小）的一行，并选择单位矩阵中行id对应的一行
从原矩阵中寻找第一个非0元素次靠左（行坐标次小）的一行，并选择单位矩阵中行id对应的一行
如果1、2步找出的行的第一个非0元素行坐标一致，说明需要通过高斯消元法相消。将第一行第一个非0元素设为a，第二行第一个非0元素设为b，第二行每个位置减去第一行每个位置对应元素乘以b/a后的结果，就可以把第二行第一个非0元素变成0，单位矩阵重复同样的内容，然后重新开始第1步
如果1、2步找出的行的第一个非0元素行坐标不一致，则将第1步找出的行从原矩阵中删除，放入新矩阵，然后重新开始第1步，单位矩阵重复同样的内容
如果找不到行，说明原矩阵中所有行以及被移入新矩阵，进入第6步
从最后一行开始，将原矩阵对应新矩阵每一行乘以第一个非0元素的倒数，再用行id比当前行更大的每一行的第一个元素来和当前行每一个非0元素相消，同时单位矩阵对应新矩阵也作同步变换，即可得到逆矩阵

public smat re()
{smat tmp = new smat(this);smat tmp2 = new smat(lineCount);smat resultTmp = new smat(lineCount);smat result = new smat(lineCount);for (int l = 0; l < lineCount; l++)resultTmp.set(l, l, 1);tmp.reLine();resultTmp.reLine();lineNode maxline = null;lineNode max2line = null;lineNode premaxline = null;lineNode premax2line = null;lineNode resMaxline = null;lineNode resMax2line = null;lineNode preresMaxline = null;lineNode preresMax2line = null;int i = 0;while (tmp.it.line != null){while (tmp.it.line != null){if (maxline == null || tmp.it.row.rowId < maxline.firstRow.nextRow.rowId){max2line = maxline;premax2line = premaxline;resMax2line = resMaxline;preresMax2line = preresMaxline;maxline = tmp.it.line;premaxline = tmp.it.preline;resMaxline = resultTmp.it.line;preresMaxline = resultTmp.it.preline;}else if ((max2line == null || tmp.it.row.rowId < max2line.firstRow.nextRow.rowId) && tmp.it.line.lineId != maxline.lineId){max2line = tmp.it.line;premax2line = tmp.it.preline;resMax2line = resultTmp.it.line;preresMax2line = resultTmp.it.preline;}tmp.it.nextLine();resultTmp.it.nextLine();}if (max2line != null && maxline.firstRow.nextRow.rowId == max2line.firstRow.nextRow.rowId){double basis = max2line.firstRow.nextRow.value / maxline.firstRow.nextRow.value;rowNode rowit = maxline.firstRow.nextRow;rowNode resRowit = resMaxline.firstRow.nextRow;while (rowit != null){add(max2line.lineId, rowit.rowId, -(rowit.value * basis));add(resMax2line.lineId, resRowit.rowId, -(resRowit.value * basis));rowit = rowit.nextRow;resRowit = resRowit.nextRow;}}else{premaxline.nextLine = maxline.nextLine;preresMaxline.nextLine = resMaxline.nextLine;maxline.nextLine = null;resMaxline.nextLine = null;if (premax2line == maxline){premax2line = premaxline;preresMax2line = preresMaxline;}rowNode j = maxline.firstRow.nextRow;rowNode k = resMaxline.firstRow.nextRow;while (j != null){tmp2.set(i, j.rowId, j.value/ maxline.firstRow.nextRow.value);j = j.nextRow;}while (k != null){result.set(i, k.rowId, k.value/ maxline.firstRow.nextRow.value);k = k.nextRow;}i++;tmp2.it.nextLine();result.it.nextLine();maxline = null;resMaxline = null;premaxline = null;preresMaxline = null;max2line = null;resMax2line = null;premax2line = null;preresMax2line = null;}tmp.reLine();resultTmp.reLine();}if (i != lineCount)return new smat(0);tmp2.reLine();result.reLine();for (int l = lineCount - 1; l >= 0; l--)for (int j = lineCount - 1; j > l; j--){for(int k=j;k<lineCount;k++)result.add(l, k, -result.get(j, k) * tmp2.get(l, j));tmp2.it.reRow();}return result;
}