矩阵变换
在一些算法中需要用到矩阵,自然就需要用到矩阵的一些操作,比如行变换、列变换、最简式、求矩阵的秩等,下面是实现的代码
public class Matrix { #region 属性 /// <summary> /// 行数 /// </summary> public int RowCount { get { return _rowCount; } } /// <summary> /// 列数 /// </summary> public int ColumnCount { get { return _columnCount; } } #endregion #region 内部变量 /// <summary> /// 行数 /// </summary> /// <remarks>Using this instead of the RowCount property to speed up calculating /// a matrix index in the data array.</remarks> private readonly int _rowCount; /// <summary> /// 列数 /// </summary> /// <remarks>Using this instead of the ColumnCount property to speed up calculating /// a matrix index in the data array.</remarks> private readonly int _columnCount; /// <summary> /// 矩阵数据 /// </summary> private readonly double[] _data; /// <summary> /// 矩阵的元素总数 /// </summary> private readonly int _count; #endregion #region 构造函数 /// <summary> /// 矩阵构造函数 /// </summary> /// <param name="rowCount">行数</param> /// <param name="columnCount">列数</param> public Matrix(int rowCount, int columnCount) { if (rowCount <= 0 || columnCount <= 0) { throw new Exception("矩阵的行数和列数必须大于0"); } _data = new double[rowCount * columnCount]; _rowCount = rowCount; _columnCount = columnCount; _count = _rowCount * _columnCount; } #endregion #region 矩阵元素的取值和赋值 /// <summary> /// 得到指定位置的矩阵元素(索引从零开始) /// </summary> public double At(int row, int column) { if (!IsLegalIndex(row, column)) { return 0xFFFFFFFF; } return _data[(column * _rowCount) + row]; } /// <summary> ///给指定位置的矩阵元素赋值(索引从零开始) /// </summary> public void At(int row, int column, double value) { if (!IsLegalIndex(row, column)) { return; } _data[(column * _rowCount) + row] = value; } #region AtRow /// <summary> /// 给指定的行赋值 /// </summary> /// <param name="row"></param> /// <param name="values"></param> public void AtRow(int row, params double[] values) { //for (int column = 0; column < _columnCount; column++) //{ // At(row, column, values[column]); //} AtRow(row, values, 0, ColumnCount); } /// <summary> /// 给指定的行赋值 /// </summary> /// <param name="row"></param> /// <param name="values"></param> /// <param name="valuesStartIndex">values数组的起始索引,通常是values的长度大于矩阵的列数时需要指定.</param> /// <param name="valuesCount">values从起始索引算起需要赋值的个数.</param> public void AtRow(int row, double[] values, int valuesStartIndex, int valuesCount) { for (int column = 0; column < valuesCount; column++) { At(row, column, values[valuesStartIndex + column]); } } /// <summary> /// 给指定的行赋值 /// </summary> /// <param name="row"></param> /// <param name="matrix">将指定矩阵的对应行的数据赋值到当前矩阵的行</param> public void AtRow(int row, Matrix matrix) { AtRow(row, matrix, row); } /// <summary> /// 给指定的行赋值 /// </summary> /// <param name="row">目标矩阵的行</param> /// <param name="rowSource">来源矩阵的行</param> /// <param name="matrixSource">将指定矩阵的对应行的数据赋值到当前矩阵的行</param> public void AtRow(int row, Matrix matrixSource, int rowSource) { if (this.ColumnCount != matrixSource.ColumnCount) { throw new Exception("指定矩阵的列数与当前矩阵的列数不相等"); } for (int column = 0; column < ColumnCount; column++) { At(row, column, matrixSource.At(rowSource, column)); } } #endregion /// <summary> /// 是否为合法的索引 /// </summary> /// <param name="row"></param> /// <param name="column"></param> /// <returns></returns> private bool IsLegalIndex(int row, int column) { if (row < 0 || column < 0) { throw new Exception("行列索引必须大于或等于0"); } if ((column * _rowCount) + row > _count - 1) { throw new Exception("行列索引超出范围"); } return true; } #endregion #region Copy Function public Matrix Copy() { Matrix target = new Matrix(this._rowCount, this._columnCount); Array.Copy(this._data, target._data, _count); return target; } #endregion #region Clear /// <summary> /// 清空矩阵.所有元素都置为0. /// </summary> public void Clear() { Array.Clear(_data, 0, _data.Length); } #endregion #region IsMostSimpleRow /// <summary> /// 是否为最简行. /// 即除了指定列和最后一列为非零外,其他列都为零. /// </summary> /// <param name="row"></param> /// <param name="notZeroColumn"></param> /// <returns></returns> public bool IsMostSimpleRow(int row, int notZeroColumn) { if (At(row, notZeroColumn) == 0 || At(row, ColumnCount - 1) == 0) {//指定列和最后一列为零 return false; } for (int column = 0; column < ColumnCount; column++) { if (column == notZeroColumn || column == ColumnCount - 1) { continue; } if (At(row, column) != 0) {//除了指定列和最后一列为非零外,其他列都为零 return false; } } return true; } #endregion #region Rank Function /// <summary> /// 计算矩阵的秩 /// </summary> /// <returns></returns> public int Rank() { //matrix为空则直接默认已经是最简形式 if (_count == 0) { return 0; } //复制一个matrix到martrix_copy,之后因计算需要改动矩阵时并不改动matrix本身 Matrix matrix_copy = this.Copy(); matrix_copy.TransformMostSimple(); //行最简矩阵的秩即为所求 int rank = matrix_copy.RankOfMostSimpleMatrix(); return rank; } #region SortByLeftNotZeroPosition /// <summary> /// 排序(按左侧最前非零位位置自上而下升序排列) /// </summary> /// <param name="matrix">矩阵</param> private void SortByLeftNotZeroPosition() { //统计每行第一个非零元素的出现位置 int[] counter = new int[_rowCount]; for (int r = 0; r < _rowCount; r++) { for (int c = 0; c < _columnCount; c++) { if (At(r, c) == 0) { counter[r]++; } else { break; } } } //按每行非零元素的出现位置升序排列 for (int r = 0; r < _rowCount; r++) { for (int j = r; j < _rowCount; j++) { if (counter[r] > counter[j]) { ExchangeRow(r, j); } } } } #endregion #region IsMostSimple /// <summary> /// 判断矩阵是否变换到最简形式(非零行数达到最少) /// </summary> /// <param name="matrix"></param> /// <returns>true:</returns> private bool IsMostSimple() { //统计每行第一个非零元素的出现位置 int[] counter = new int[_rowCount]; for (int r = 0; r < _rowCount; r++) { for (int c = 0; c < _columnCount; c++) { if (At(r, c) == 0) { counter[r]++; } else break; } } //后面行的非零元素出现位置必须在前面行的后面,全零行除外 for (int i = 1; i < counter.Length; i++) { if (counter[i] <= counter[i - 1] && counter[i] != _columnCount) { return false; } } return true; } #endregion #region ElementaryTrasform Function /// <summary> /// 行初等变换(左侧最前非零位位置最靠前的行,只保留一个) /// </summary> /// <param name="matrix">矩阵</param> private void ElementaryTrasform() { //统计每行第一个非零元素的出现位置,数值从1开始. int[] firstNotZeroPosArr = new int[_rowCount]; for (int r = 0; r < _rowCount; r++) { for (int c = 0; c < _columnCount; c++) { if (At(r, c) == 0) { firstNotZeroPosArr[r]++; } else { break; } } } for (int row = 1; row < firstNotZeroPosArr.Length; row++) { if (firstNotZeroPosArr[row] == firstNotZeroPosArr[row - 1] && firstNotZeroPosArr[row] != _columnCount) {//该行的非零位置与上面一行的非零位置相同,并且不是最后一列. //上面一行非零位置的值 double upRowNotZeroValue = At(row - 1, firstNotZeroPosArr[row - 1]); //当前行非零位置的值 double currentRowNotZeroValue = At(row, firstNotZeroPosArr[row]); At(row, firstNotZeroPosArr[row], 0);//将当前位置的值设为0. for (int j = firstNotZeroPosArr[row] + 1; j < _columnCount; j++) { //上面一行非零位置的下一个位置的值 double upRowNotZeroNextPosValue = At(row - 1, j); double value = At(row, j) - (upRowNotZeroNextPosValue * currentRowNotZeroValue / upRowNotZeroValue); At(row, j, value); } break; } } } #endregion #region AssignZeroAlmost /// <summary> /// 将和0非常接近的数字视为0 /// </summary> /// <param name="matrix"></param> private void AssignZeroAlmost() { for (int i = 0; i < _count; i++) { if (Math.Abs(_data[i]) <= 0.00001) { _data[i] = 0; } } } #endregion #region RankOfMostSimpleMatrix Function /// <summary> /// 计算行最简矩阵的秩 /// </summary> /// <param name="matrix"></param> /// <returns></returns> public int RankOfMostSimpleMatrix() { int rank = -1; bool isAllZero = true; for (int r = 0; r < _rowCount; r++) { isAllZero = true; //查看当前行有没有0 for (int j = 0; j < _columnCount; j++) { if (At(r, j) != 0) { isAllZero = false; break; } } //若第i行全为0,则矩阵的秩为i if (isAllZero) { rank = r; break; } } //满秩矩阵的情况 if (rank == -1) { rank = _rowCount; } return rank; } #endregion #endregion #region TransformMostSimple Function /// <summary> /// 变换为最简矩阵 /// </summary> public void TransformMostSimple() { //先以最左侧非零项的位置进行行排序 SortByLeftNotZeroPosition(); //循环化简矩阵 while (!IsMostSimple()) { ElementaryTrasform(); SortByLeftNotZeroPosition(); } //过于趋近0的项,视作0,减小误差 AssignZeroAlmost(); } #endregion #region ExchangeRow Function /// <summary> /// 交换矩阵的行 /// </summary> /// <param name="sourceRow">源行</param> /// <param name="targetRow">目标行</param> public void ExchangeRow(int sourceRow, int targetRow) { double sourceTemp = 0; double targetTemp = 0; for (int c = 0; c < _columnCount; c++) { sourceTemp = At(sourceRow, c); targetTemp = At(targetRow, c); At(sourceRow, c, targetTemp); At(targetRow, c, sourceTemp); } } #endregion #region ExchangeColumn Function /// <summary> /// 交换矩阵的列 /// </summary> /// <param name="sourceColumn">源行</param> /// <param name="targetColumn">目标行</param> public void ExchangeColumn(int sourceColumn, int targetColumn) { double sourceTemp = 0; double targetTemp = 0; for (int row = 0; row < RowCount; row++) { sourceTemp = At(row, sourceColumn); targetTemp = At(row, targetColumn); At(row, sourceColumn, targetTemp); At(row, targetColumn, sourceTemp); } } #endregion #region OfRowArray Function /// <summary> /// 得到行数组 /// </summary> /// <param name="rowIndex"></param> /// <returns></returns> public double[] OfRowArray(int rowIndex) { double[] rowArray = new double[_columnCount]; for (int c = 0; c < _columnCount; c++) { rowArray[c] = At(rowIndex, c); } return rowArray; } #endregion #region OfColumnArray Function /// <summary> /// 得到列数组 /// </summary> /// <param name="rowIndex"></param> /// <returns></returns> public double[] OfColumnArray(int columnIndex) { double[] columnArray = new double[_rowCount]; for (int r = 0; r < _rowCount; r++) { columnArray[r] = At(r, columnIndex); } return columnArray; } #endregion #region ToString public override string ToString() { //return base.ToString(); StringBuilder sb = new StringBuilder(); for (int r = 0; r < _rowCount; r++) { for (int c = 0; c < _columnCount; c++) { //sb.Append(String.Format("{0:000.000},", At(r, c))); //sb.Append(String.Format("{0:00.0},", At(r, c))); sb.Append(At(r, c) + ","); } sb.Append(";" + Environment.NewLine); } return sb.ToString(); } #endregion #region TransformRow /// <summary> /// 行变换(row1=row1*factor1-row2*factor2) /// </summary> /// <param name="row1"></param> /// <param name="row2"></param> /// <param name="factor1">行1乘以的系数</param> /// <param name="factor2">行2乘以的系数</param> public void TransformRow(int row1, int row2, double factor2) { double value1 = 0; double value2 = 0; for (int column = 0; column < ColumnCount; column++) { value1 = At(row1, column); value2 = At(row2, column) * factor2; At(row1, column, value1 - value2); } } #endregion #region TransformColumn /// <summary> /// 列变换(column1=column1*factor1-column2*factor2) /// </summary> /// <param name="column1"></param> /// <param name="column2"></param> /// <param name="factor1">列1乘以的系数</param> /// <param name="factor2">列2乘以的系数</param> public void TransformColumn(int column1, int column2, double factor2) { double value1 = 0; double value2 = 0; for (int row = 0; row < RowCount; row++) { value1 = At(row, column1); value2 = At(row, column2) * factor2; At(row, column1, value1 - value2); } } #endregion #region Assign /// <summary> /// 将源矩阵的值赋入当前矩阵中 /// </summary> /// <param name="matrixSource"></param> public void Assign(Matrix matrixSource) { for (int row = 0; row < matrixSource.RowCount; row++) { AtRow(row, matrixSource); } } /// <summary> /// 将源矩阵的值赋入当前矩阵中 /// </summary> /// <param name="matrixSource"></param> /// <param name="skipRowOfSource">源矩阵中需要跳过的行</param> public void Assign(Matrix matrixSource, int skipRowOfSource) { int rowNew = 0; for (int row = 0; row < matrixSource.RowCount; row++) { if (row == skipRowOfSource) { continue; } AtRow(rowNew, matrixSource, row); rowNew++; } } #endregion #region EqualeValue public bool EqualeValue(Matrix compareMatrix) { if (compareMatrix == null || compareMatrix.RowCount != this.RowCount || compareMatrix.ColumnCount != this.ColumnCount) { return false; } for (int row = 0; row < RowCount; row++) { for (int column = 0; column < ColumnCount; column++) { if (At(row, column) != compareMatrix.At(row, column)) { return false; } } } return true; } #endregion #region 重载运行符 public static Matrix operator *(Matrix source, double factor) { for (int r = 0; r < source.RowCount; r++) { for (int c = 0; c < source.ColumnCount; c++) { source.At(r, c, source.At(r, c) * factor); } } return source; } #endregion #region RowMulti /// <summary> /// 行乘以一个系数 /// </summary> /// <param name="row"></param> /// <param name="factor"></param> public void RowMulti(int row, double factor) { for (int c = 0; c < ColumnCount; c++) { At(row, c, At(row, c) * factor); } } #endregion #region ColumnMulti /// <summary> /// 列乘以一个系数 /// </summary> /// <param name="column"></param> /// <param name="factor"></param> public void ColumnMulti(int column, double factor) { for (int row = 0; row < RowCount; row++) { At(row, column, At(row, column) * factor); } } #endregion }转载请注明出处