最近花了一段时间去温习矩阵的基础算法,要去敲代码的时候才发现原来以前的线性代码全都白学了,没有去用没有去实现的结果就等于什么都没用。
自己为了练习实现了一些矩阵的基础算法,算是一个温习,一条条罗列下来。其中我的矩阵数据为浮点数格式,0这样的精度判断上(+/-0.000001f之间)可能还有些局限,仅供学习参考,如果发现跑不通的用例也欢迎及时反馈给我。 首先是矩阵基本格式的定义,比较简单。
<pre name="code" class="csharp">public class Matrix
{
public const float MAX = 10000000f;
public const float ZERO = 0.000001f; //极小值,认为是0
//空对象
private Matrix NULL = null;
public Dictionary<int, List<float>> data = new Dictionary<int, List<float>>();
public Matrix()
{
}
/// <summary>
/// 带参构造:float二维数组构造一个矩阵
/// ex:float[,] data = new float[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
/// </summary>
/// <param name="data"></param>
public Matrix(float[,] dataArray)
{
for (int index = 0; index < dataArray.GetLength(0); index++)
{
List<float> list = new List<float>();
for (int subindex = 0; subindex < dataArray.GetLength(0); subindex++)
{
list.Add(dataArray[index, subindex]);
}
this.data.Add(index, list);
}
}
}
如上是矩阵的默认默认构造和基本构造函数,传入的数据格式为一个二维的浮点矩阵。先来介绍具体的实现:
(1)矩阵相乘矩阵A和B相乘的基本思路(要满足A行等于B列的条件):结果矩阵的i行等于A的i行所有元素和B的i列所有元素依次乘积的和。结果矩阵的行列就为A行B列。开始写的是个矩阵类的成员函数,后来觉得这样很不好用,还是重载了*运算,可以使矩阵直接相乘运算matrest=mata*matb这样:
/// <summary>
/// 重载矩阵乘法*,
/// </summary>
/// <param name="matrixA"></param>
/// <param name="matrixB"></param>
/// <returns></returns>
public static Matrix operator *(Matrix matrixA, Matrix matrixB)
{
Matrix resMatrix = new Matrix();
List<int> keyList = new List<int>(matrixA.data.Keys);
Dictionary<int, List<float>> dic = matrixA.data;//[keyList[inIndex]];
//现在是A矩阵的inIndex行,要分别乘以B的inIndex列
//遍历没一列
for (int lineIndex = 0; lineIndex < matrixA.LieCount; lineIndex++)
{
List<float> resList = new List<float>();
//即将遍历每一行
for (int rowIndex = 0; rowIndex < matrixA.HangCount; rowIndex++)
{
float rowNum = 0;
//将A行rowIndex行subIndex列的数字和B中rowIndex列subIndex行的数字相乘
List<float> aRow = matrixA.GetHang(lineIndex); //A行的数据
List<float> bLine = matrixB.GetLie(rowIndex);//B列的数据
for (int subIndex = 0; subIndex < aRow.Count; subIndex++)
{
rowNum += aRow[subIndex] * bLine[subIndex];
}
resList.Add(rowNum);
}
resMatrix.data.Add(lineIndex, resList);
}
return resMatrix;
}
(2)矩阵相加:
矩阵A+B相对简单,直接每个元素相加的结果即可,需要保证A与B的行列都相同:
/// <summary>
/// 重载+号,计算矩阵相加
/// </summary>
/// <param name="lMat"></param>
/// <param name="rmat"></param>
/// <returns></returns>
public static Matrix operator +(Matrix lMat, Matrix rMat)
{
if (lMat.HangCount != rMat.HangCount)
{
return null;
}
if (lMat.LieCount != rMat.LieCount)
{
return null;
}
//计算 矩阵相加
Matrix resMat = new Matrix();
List<int> keyList = new List<int>(lMat.data.Keys);
for (int index = 0; index < keyList.Count; index++)
{
List<float> lines = lMat.data[keyList[index]];
List<float> linesB = rMat.data[keyList[index]];
List<float> resList = new List<float>();
for (int subIndex = 0; subIndex < lines.Count; subIndex++)
{
resList.Add(lines[subIndex] + linesB[subIndex]);
}
resMat.data.Add(index, resList);
}
return resMat;
}
(3)求矩阵的模,矩阵的模仅针对方针讨论,同时2阶矩阵的处理有别于多阶矩阵,注意点都注明了,同时overflow字段的结果为true的话也就表明模也就失效,表示不是矩阵:
/// <summary>
/// 计算矩阵的模
/// </summary>
/// <returns></returns>
public float GetModule(out bool overflow)
{
overflow = false;
//行数不等于列数,
if (LieCount != HangCount)
{
overflow = true;
return 0;
}
//2阶矩阵的模为0 啊
///a b
///c d
///模= a*d - c*b
if (LieCount == 2)
{
return getMultioflist(getRightDownList(0)) - getMultioflist(getLeftUpList(0));
}
//测试通过
List<float> rightDown = new List<float>();//右下各个列的积
for (int index = 0; index < LieCount; index++)
{
rightDown.Add(getMultioflist(getRightDownList(index)));
}
List<float> leftUp = new List<float>();
for (int index = 0; index < LieCount; index++)
{
leftUp.Add(getMultioflist(getLeftUpList(index)));
}
float result = getSumofList(rightDown) - getSumofList(leftUp);
return result;
}
(4)求矩阵的转置:
矩阵的转置的结果就是原来的i行变成了i列:
/// <summary>
/// 获得转置矩阵,不改变原矩阵
/// </summary>
/// <returns></returns>
public Matrix GetTranspose()
{
//非方阵
if (HangCount != LieCount)
{
return null;
}
Matrix transpose = new Matrix();
List<int> keyList = new List<int>(data.Keys);
for (int index = 0; index < keyList.Count; index++)
{
List<float> hang = new List<float>(data[keyList[index]]);
for (int subIndex = 0; subIndex < hang.Count; subIndex++)
{
if (transpose.data.ContainsKey(subIndex))
{
transpose.data[subIndex].Add(hang[subIndex]);
}
else
{
transpose.data.Add(subIndex, new List<float> { hang[subIndex] });
}
}
}
return transpose;
}
(5)求矩阵的标准伴随矩阵:
矩阵A的标准伴随矩阵,就是A的代数余子式矩阵(所有行列的代数余子式的模组成的矩阵)转置得到的,得到标准伴随矩阵的方法为 GetStanderedAdjoint(),最重要的还是求标准伴随矩阵的方法。
其中i行j列的代数余子式的矩阵其实是原矩阵A去掉i行j列的结果矩阵,而代数余子式矩阵就是由这些余子式的模*符号sign 组成的矩阵。具体符号sign的值为1要求i+j为偶数,否则sign为-1。
//A的标准伴随矩阵
public Matrix GetStanderedAdjoint()
{
return GetConfactorMat().GetTranspose();
}
//获取代数余子式组成的矩阵
public Matrix GetConfactorMat()
{
Matrix resMat = new Matrix();
for (int hangInd = 0; hangInd < HangCount; hangInd++)
{
List<float> list = GetHang(hangInd);
List<float> resList = new List<float>();
for (int lieInd = 0; lieInd < list.Count; lieInd++)
{
//获取index行,subindex列的代数余子式的模啊
float sign = (hangInd + lieInd) % 2 == 0 ? 1f : -1f;//结果行列的符号
float mode = GetConfactor(hangInd, lieInd).Mode; //.Mode为矩阵的模
resList.Add(sign * mode);
}
resMat.data.Add(hangInd, resList);
}
return resMat;
}
(6)求矩阵的逆矩阵,都不用多说,代码中都已经有详细的算法注明。除了MultiScalar(float)方法是没有粘出来的矩阵乘以浮点数。
/// <summary>
/// 获得一个矩阵的逆矩阵:如果是,返回的逆矩阵不为空
/// 注:如果一个矩阵有逆矩阵,那矩阵和逆矩阵相乘结果等于单位矩阵
/// </summary>
/// <returns></returns>
public Matrix GetInverseMatrix()
{
//step1:先筛选:某行或某列全为0的矩阵没有逆矩阵(因为乘自己得到的结果是0矩阵)
if (CheckZeroList() == false)
return null;
//几个简化原则:
//1.单位矩阵的逆矩阵是本身
if (GetUnitMatrixN(this.HangCount) == this)
return this;
///其他性质
///2.A是奇异矩阵的话,他的逆矩阵的逆等于自己
///3.矩阵转置的逆等于矩阵的逆的转置
///4.(AB)-1 = (B-1)*(A-1) -1表示-1次方,即矩阵乘积的逆等于反顺序的矩阵逆的乘积
//step2:先得到:代数余子式的矩阵
//GetConfactorMat();
//step3:余子式的矩阵要转置:得到A的标准伴随矩阵
Matrix resMat = GetConfactorMat().GetTranspose();
//step3:矩阵除以原矩阵的模
float moderate = 1f / this.Mode;
//Log.MyDebug("原矩阵的1/模=" + moderate);
return resMat.MultiScalar(moderate);
}
(7)判断一个矩阵是否正交:
一个矩阵如果是正交的,那它的转置矩阵和自身相乘的结果(矩阵为N行)应该等于N阶单位矩阵。其中GetUnitMatrixN(int n)方法为获取n阶的单位矩阵:
//检查一个矩阵是否正交?
public bool CheckOrthogonal()
{
Matrix lMat = this;
if (lMat == NULL)
return false;
if (lMat.HangCount == 1 || lMat.LieCount == 1)
return false;
Matrix rMat = this.GetTranspose();
if (rMat == NULL)
return false;
Matrix res = lMat * rMat;
if (res == NULL)
return false;
Matrix resOther = GetUnitMatrixN(res.HangCount);
return (res == resOther);
}
最后附上工程的完整地址,我的github地址:
https://github.com/yuhezhangyanru/UnityStudy/tree/master/UnityStudy/Study1