数学：《线性代数》方阵求逆 之 代码实现

背景

给定坐标，我可以使用变换方阵对其进行右乘进行坐标转换，如果使用一个堆栈记录下执行过的变换方阵，然后执行一系列出栈操作，对出栈的方阵求逆，然后使用逆方阵右乘坐标就可以得到之前的坐标。

方阵求逆

理论

求逆

对 A:E 进行矩阵的最简算法得出 E:A-¹。

证明：

D1D2DkA = E

D1D2DkAA-¹ = EA-¹

D1D2DkE = A-¹

上面的三个公式说明，对 A 做有限的初等变换可以得到 E，同样对 E 做同样的初等变换可得到 A-¹，因此 我们求 A:E 的最简矩阵就可以得出  A-¹。

代码实现

核心代码

            public Matrix Inverse()
            {
                if (_m != _n)
                {
                    throw new InvalidOperationException("只有方阵可以求逆");
                }

                var combinedMatrix = this.Combine(this.E());
                combinedMatrix = combinedMatrix.Simplify();

                var resultMatrix = new Matrix(_m, _m);
                for (var i = 0; i < resultMatrix._m; i++)
                {
                    Array.Copy(combinedMatrix._data[i], resultMatrix._m, resultMatrix._data[i], 0, resultMatrix._m);
                }

                return resultMatrix;
            }

全部代码

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace DataStuctureStudy.Matrixes
{
    class MatrixTest
    {
        public static void Test()
        {
            var matrixA =
                new Matrix(3, 3)
                .SetRowData(0, 1, 2, 3)
                .SetRowData(1, 1, 2, 7)
                .SetRowData(2, 4, 9, 2);

            var matrixB = matrixA.Simplify();
            var matrixC = matrixA.Inverse();

            matrixA.Dispaly();
            matrixB.Dispaly();
            matrixC.Dispaly();
            matrixA.Multiply(matrixC).Dispaly();

        }
        class Matrix
        {
            private readonly int _m;
            private readonly int _n;
            private readonly double[][] _data;

            public Matrix(int m, int n)
            {
                _m = m;
                _n = n;
                _data = new double[m][];
                for (var i = 0; i < m; i++)
                {
                    _data[i] = new double[n];
                }
            }

            public Matrix SetRowData(int row, params double[] values)
            {
                Array.Copy(values, _data[row], _n);

                return this;
            }

            public Matrix SwapRow(int rowX, int rowY)
            {
                var resultMatrix = this.ConeSelf();

                var temp = resultMatrix._data[rowX];
                resultMatrix._data[rowX] = resultMatrix._data[rowY];
                resultMatrix._data[rowY] = temp;

                return resultMatrix;
            }

            public Matrix MultiplyRow(int row, double operand)
            {
                var resultMatrix = this.ConeSelf();

                for (var i = 0; i < resultMatrix._n; i++)
                {
                    resultMatrix._data[row][i] *= operand;
                }

                return resultMatrix;
            }


            public Matrix AddSourceRowToTargetRow(int sourceRow, double operand, int targetRow)
            {
                var resultMatrix = this.ConeSelf();

                for (var i = 0; i < resultMatrix._n; i++)
                {
                    resultMatrix._data[targetRow][i] += resultMatrix._data[sourceRow][i] * operand;
                }

                return resultMatrix;
            }

            public Matrix Simplify()
            {
                var resultMatrix = this.ConeSelf();

                for (var col = 0; col < resultMatrix._n; col++)
                {
                    if (col >= _m)
                    {
                        break;
                    }

                    if (resultMatrix._data[col][col] == 0)
                    {
                        var nonZeroRowIndex = resultMatrix.FindNonZeroRowIndex(col, col + 1);
                        if (nonZeroRowIndex == -1)
                        {
                            break;
                        }
                        resultMatrix = resultMatrix.SwapRow(col, nonZeroRowIndex);
                    }

                    resultMatrix = resultMatrix.MultiplyRow(col, 1 / resultMatrix._data[col][col]);

                    for (var row = 0; row < _m; row++)
                    {
                        if (row == col)
                        {
                            continue;
                        }
                        resultMatrix = resultMatrix.AddSourceRowToTargetRow(col, -1 * resultMatrix._data[row][col], row);
                    }
                }

                return resultMatrix;
            }

            public Matrix Multiply(Matrix matrix)
            {
                if (this._n != matrix._m)
                {
                    throw new InvalidOperationException("左侧矩阵的列不等于右侧矩阵的行");
                }

                var resultMatrix = new Matrix(this._m, matrix._n);

                for (var i = 0; i < resultMatrix._m; i++)
                {
                    for (var j = 0; j < resultMatrix._n; j++)
                    {
                        resultMatrix._data[i][j] = Enumerable.Range(0, this._n).Sum(k => this._data[i][k] * matrix._data[k][j]);
                    }
                }

                return resultMatrix;
            }

            public Matrix Inverse()
            {
                if (_m != _n)
                {
                    throw new InvalidOperationException("只有方阵可以求逆");
                }

                var combinedMatrix = this.Combine(this.E());
                combinedMatrix = combinedMatrix.Simplify();

                var resultMatrix = new Matrix(_m, _m);
                for (var i = 0; i < resultMatrix._m; i++)
                {
                    Array.Copy(combinedMatrix._data[i], resultMatrix._m, resultMatrix._data[i], 0, resultMatrix._m);
                }

                return resultMatrix;
            }

            public Matrix Combine(Matrix matrix)
            {
                if (this._m != matrix._m)
                {
                    throw new InvalidOperationException("左侧矩阵的行不等于右侧矩阵的行");
                }

                var resultMatrix = new Matrix(this._m, this._n + matrix._n);
                for (var i = 0; i < resultMatrix._m; i++)
                {
                    Array.Copy(this._data[i], 0, resultMatrix._data[i], 0, this._data[i].Length);
                    Array.Copy(matrix._data[i], 0, resultMatrix._data[i], this._data[i].Length, matrix._data[i].Length);
                }

                return resultMatrix;
            }

            public Matrix E()
            {
                return E(_m);
            }

            public static Matrix E(int m)
            {
                var matrix = new Matrix(m, m);
                for (var i = 0; i < m; i++)
                {
                    matrix._data[i][i] = 1;
                }

                return matrix;
            }

            public Matrix ConeSelf()
            {
                var resultMatrix = new Matrix(this._m, this._n);
                for (var i = 0; i < this._m; i++)
                {
                    Array.Copy(this._data[i], resultMatrix._data[i], this._n);
                }

                return resultMatrix;
            }

            public void Dispaly()
            {
                Console.WriteLine();
                foreach (var row in _data)
                {
                    Console.WriteLine("[ " + String.Join(" , ", row) + " ]");
                }
                Console.WriteLine();
            }

            private int FindNonZeroRowIndex(int col, int startRow)
            {
                if (startRow >= _m)
                {
                    return -1;
                }

                for (var i = startRow; i < _m; i++)
                {
                    if (_data[i][col] != 0)
                    {
                        return i;
                    }
                }

                return -1;
            }
        }
    }
}