介绍矩阵的编程表示，矩阵的初等变换，化阶梯矩阵及求方阵的行列式。
摘要
介绍矩阵的编程表示，矩阵的初等变换，化阶梯矩阵及求方阵的行列式。
矩阵介绍
矩阵就是一个M×N平面的表格，用一个两维数组就可以表示,为了输入方便，我们用一个特殊格式的字符串对矩阵进行初始化，用"|"分割每一行，用","分割每一列，并增加一个Show的方法打印出矩阵，为了测试及调试的需要，还重写了ToString()方法。
代码
public class Matrix {
public int M { get; private set; }
public int N { get; private set; }
private readonly double[,] num = null;
#region 构造函数/Show
public Matrix(int m, int n, string input) {
M = m;
N = n;
num = new double[m, n];
if (!string.IsNullOrEmpty(input))
parseInput(input);
}
private void parseInput(string input) {
string[] rows = input.Split(new[] { '|' });
if (rows.Length != M)
throw new ArgumentException("row count err");
for (int i = 0; i < M; i++) {
string row = rows[i];
string[] cells = row.Split(new[] { ',' });
if (cells.Length != N)
throw new ArgumentException(string.Format("cells counte err:{0}", row));
for (int j = 0; j < N; j++) {
int cellValue;
if (!int.TryParse(cells[j], out cellValue))
throw new ArgumentException(string.Format("cell error:{0}", cells[j]));
num[i, j] = cellValue;
}
}
}
public void Show() {
for (int i = 0; i < M; i++) {
for (int j = 0; j < N; j++)
Console.Write("{0}\t", num[i, j]);
Console.WriteLine();
}
}
public override string ToString() {
StringBuilder sb = new StringBuilder();
for (int i = 0; i < M; i++) {
for (int j = 0; j < N; j++) {
sb.Append(num[i, j]);
if (j != N - 1)
sb.Append(',');
}
if (i != M - 1)
sb.Append('|');
}
return sb.ToString();
}
}
先对这些代码进行单元测试，为了让代码块覆盖率达到100%,对构造函数的测试要故意模拟错误的输入，以便覆盖抛出异常的语句。
代码
[TestMethod()]
public void MatrixConstructorTest() {
Matrix m;
try
{
m = new Matrix(2,2,"1,2|3,4|5,6");
Assert.Fail();
}
catch{}
try {
m = new Matrix(2, 2, "1,2|3,4,5");
Assert.Fail();
}
catch { }
try {
m = new Matrix(2, 2, "1,2|3,a");
Assert.Fail();
}
catch { }
}
[TestMethod]
public void toStringTest()
{
Matrix matrix = new Matrix(3, 3, "1,2,3|4,5,6|7,8,9");
Assert.IsTrue(matrix.ToString() == "1,2,3|4,5,6|7,8,9");
}
[TestMethod()]
public void ShowTest() {
Matrix matrix = new Matrix(2, 2, "1,2|3,4");
matrix.Show();
//这玩意还真没法儿测，不出错就算成功
}
初等行变化
矩阵的初等行变换虽然是很简单的变换，但用处非常大，包括行交换（两行交换位置），行倍乘（给某行乘以一个非零常数），行倍（某行乘以一个非零常数加到另一行上）加三种变换，因为行倍乘在求阶梯矩阵时用不到，就不写了。
代码
//r1→r2,r2→r1
public void swapRow(int r1, int r2) {
double[] temp = new double[N];
for (int i = 0; i < N; i++)
temp[i] = num[r1, i];
for (int i = 0; i < N; i++)
num[r1, i] = num[r2, i];
for (int i = 0; i < N; i++)
num[r2, i] = temp[i];
}
//c×r1+r2
public void double_add(int r1, int r2, double c) {
for (int i = 0; i < N; i++) {
num[r2, i] += num[r1, i] * c;
}
}
单元测试如下，能达到预期就可以。
代码
[TestMethod()]
public void swapRowTest() {
Matrix matrix = new Matrix(2,2,"1,2|3,4");
matrix.swapRow(0, 1);
Assert.IsTrue(matrix.ToString() == "3,4|1,2");
}
[TestMethod()]
public void double_addTest() {
Matrix matrix = new Matrix(2,2,"1,2|3,4");
matrix.double_add(0,1,2);
Assert.IsTrue(matrix.ToString() == "1,2|5,8");
}
化阶梯矩阵
阶梯矩阵在矩阵算法中有很重要的作用，比如在解线性方程时，先化为阶梯阵，再从未知数最少的方程逐步回代求解，或者求矩阵的秩的时候先化成阶梯矩阵，再去数非零的行数。
任何一个矩阵都可以转换成阶梯矩阵，方法是一列一列去转换。
注：A表示矩阵，_标识下标，第一个下标表示行,第二个下标表示列
1、从A_1_1开始，假设A_1_1非0（如果为0，则和其他首非零行交换），把A_1_1记做a
2、如果A_2_1为非0，记做b，那么第一行乘以-(b/a)加到第二行上，这时第二行首为零。
3、依次类推，把从A_2_1,A_3_1,...,A_m_1都化为0.
4、转向A_2_2，假设假设A_2_2非0（如果为0，则和其他第二列非零行交换），把A_2_2记做a
5、如果A_3_2为非0，记做b，那么第一行乘以-(b/a)加到第三行上，这时第三行第二为零。
6、依次类推，把从A_3_2,A_4_2,...,A_m_2都化为0.
7、依次类推，没重复一次1-3步的操作，都会把一列的指定行下面均化为0，最终可得到阶梯矩阵。
代码
public void toStrassen() {
int start_row = 0, start_col = 0;
for (int col = 0; col < N; col++) {//一列一列去处理
var no_zero_col = get_no_zero_col(start_row, start_col);
if (no_zero_col == -1) return; //从该行该列向下向右全为0,直接返回
make_col_zero(no_zero_col, start_row);
start_row++;
if (no_zero_col == N - 1)
return;
start_col = ++no_zero_col;
}
}
//通过倍加运算使某一列从某个起始行开始下面的数为0
private void make_col_zero(int col, int start_row) {
for (int row = start_row + 1; row < M; row++) {
var z = num[row, col] / num[start_row, col];
double_add(start_row, row, -z);
}
}
// 获取从指定起始行和起始列开始的非零列。
//如果在起始列上从起始行开始找到非零行，但不是给定起始行，那么利用初等变化交换这两行。
// 如果在起始列上未找到非零行，那么起始行不变，起始列右移重新计算
private int get_no_zero_col(int start_row, int start_col) {
for (int row = start_row; row < M; row++) {
if (num[row, start_col] == 0) continue;
swapRow(row, start_row);
return start_col;
}
for (int col = start_col + 1; col < N; col++) {
for (int row = start_row; row < M; row++) {
if (num[row, col] == 0) continue;
swapRow(start_row, row);
return col;
}
}
return -1;
}
这个用例的单元测试稍微复杂，因为有很多分支，所以要考虑多种情况，比如矩阵本来就是阶梯矩阵，或者有0行，或者有多列只有前N行为0，对角线均为0等情况。
代码
[TestMethod()]
public void toStrassenTest() {
Matrix matrix = new Matrix(4, 4, "0,1,-1,1|2,-2,1,1|2,3,-5,7|2,16,-14,24");
matrix.toStrassen();
Assert.IsTrue(matrix.ToString() == "2,-2,1,1|0,1,-1,1|0,0,-1,1|0,0,0,8");
matrix = new Matrix(3, 5, "0,0,0,5,6|1,2,3,4,5|0,0,0,10,8");
matrix.toStrassen();
Assert.IsTrue(matrix.ToString() == "1,2,3,4,5|0,0,0,5,6|0,0,0,0,-4");
matrix = new Matrix(3, 3, "1,2,3|4,0,0|0,0,0");
matrix.toStrassen();
Assert.IsTrue(matrix.ToString() == "1,2,3|0,-8,-12|0,0,0");
}
求方阵的行列式
只有方阵才可以求行列式的值，求行列式的值是一个多项式想加算法，公式如下
其中a是n阶方阵对应的行列式，τ表示对一个数列求逆序数，j1,j2,...jn是列下标，大概意思是行下标不变，是1,2,3,...,n，列下标是一个1,2,3,...,n的一个全排列（N级排列），然后每个多项式是n个分量相乘，这个分量的行下标是1-n顺序排列，列下标是n级排列的其中一个。
这里用到了N级排列，所以先写个辅助类来求N级排列
代码
//计算一个N级排列
public static IEnumerable<int[]> permute(int n) {
if (n < 0) throw new ArgumentException();
List<int[]> ret = new List<int[]>();
int[] used = new int[n];
int[] p = new int[n];
int[] data = new int[n];
for (int i = 0; i < n; i++)
data[i] = i;
permute(ret, n, 0, used, p, data);
return ret;
}
//从一个数的集合里选注N个数的排列组合，使用回溯算法，复杂度应该是O(n!)
//http://blog.chinaunix.net/u/10638/showart.php?id=52814
private static void permute(ICollection<int[]> list, int n, int pos, int[] used, int[] p, int[] data) {
if (pos == n) {
int[] temp = new int[n];
p.CopyTo(temp, 0);
list.Add(temp);
return;
}
for (int i = 0; i < n; i++) {
if (used[i] != 0) continue;
used[i]++;
p[pos] = data[i];
permute(list, n, pos + 1, used, p, data);
used[i]--;
}
}
核心算法是从网上找的，先对其进行单元测试
代码
[TestMethod()]
public void permuteTest() {
try
{
MathHelper.permute(-1);
}
catch (Exception ex)
{
Assert.IsInstanceOfType(ex, typeof(ArgumentException));
}
const int n = 3;
IEnumerable<int[]> actual = MathHelper.permute(n);
int len = 0;
List<string> expected = new List<string>() { "012", "021", "102", "120", "201", "210" };
List<string> actualstrs = new List<string>();
foreach (int[] ints in actual) {
string temp = string.Empty;
for (int i = 0; i < ints.Length; i++)
temp += ints[i].ToString();
actualstrs.Add(temp);
len++;
}
Assert.IsTrue(len == 6);
foreach (string s in expected)
Assert.IsTrue(actualstrs.Contains(s));
}
求一个数列中的逆序数，直接用遍历数一遍就行了，不过也得做单元测试。
代码
//计算一个排列中的逆序数，O(n^2)复杂度，没做优化
public static int ReverseNumber(int[] cols) {
if (cols == null) throw new ArgumentNullException();
int ret = 0;
int len = cols.Length;
for (int i = 0; i < len; i++) {
for (int j = i + 1; j < len; j++)
if (cols[i] > cols[j])
ret++;
}
return ret;
}
2, 4, 3, 1, 5这个排列中，2后面有1，4后面有3和1，3后面有1，共4个逆序，用例做单元测试的数据。
代码
[TestMethod()]
public void ReverseNumberTest() {
try {
MathHelper.ReverseNumber(null);
}
catch (Exception ex) {
Assert.IsInstanceOfType(ex, typeof(ArgumentNullException));
}
int[] cols = new[] { 2, 4, 3, 1, 5 };
const int expected = 4;
int actual = MathHelper.ReverseNumber(cols);
Assert.AreEqual(expected, actual);
}
有了两个辅助方法，就是用代码把公式写出来就行了
代码
//∑_（j1,j2,…,jn）?〖(-1)^τ(j1,j2,…,jn) a_(1j_1 ) a_(2j_2 ) 〗…a_(nj_n )
public double detA() {
if (M != N)
throw new NotSupportedException();
double sum = 0.0D;
var permute = MathHelper.permute(M);
foreach (int[] cols in permute) {
int reverse_number = MathHelper.ReverseNumber(cols);
double temp = 1.0D;
for (int i = 0; i < M; i++) {
int j = cols[i];
temp *= num[i, j];
}
if (reverse_number % 2 == 0)
sum += temp;
else
sum -= temp;
}
return sum;
}
也做下单元测试
代码
[TestMethod()]
public void detATest() {
const int m = 4;
const int n = 4;
const string input = "4,-5,10,3|1,-1,3,1|2,-4,5,2|-3,2,-7,-1";
Matrix target = new Matrix(m, n, input);
const double expected = 1D;
double actual = target.detA();
Assert.AreEqual(expected, actual);
Matrix matrix = new Matrix(3, 2, "1,2|3,4|5,6");
try
{
matrix.detA();
Assert.Fail();
}
catch{}
}
小节
这里演示了写一个简单的类库，并用单元测试来保证其质量的过程，矩阵还涉及到其它算法，数乘，乘法，求逆，求特征值，特征向量等，慢慢再实现。
其中求特征值是一个解一元N次方程的问题，周末用二分法弄了好几个小时也没弄出来，这里求一下代码，呵呵。
