矩阵求逆-高斯消元法介绍及其实现

数学基础之矩阵系列

1. 矩阵求逆-高斯消元法介绍及其实现
2. 矩阵求行列式-高斯消元法实现

---

矩阵求逆在实际问题中经常遇到，根据定义，对于任一个矩阵$A_{n\times n}$，其逆矩阵和其本身满足$A^{-1}A=E$，其中矩阵$E$ 是单位矩阵。

高斯消元法是一种求矩阵逆比较高效的方法，其方法是对于一个矩阵$A_{n\times n}$，先构造一个增广矩阵$W=[A|E]$，其中$E$ 是一个$n\times n$的单位矩阵，这样$W$就成了一个$n\times 2n$的矩阵。接下来对$W$进行行变换，使之变成$[E|B]$的形式，这样就可以确定$A^{-1}=B$。但为什么矩阵$B$就为矩阵$A$的逆呢，现证明如下：

因为是仅对增广矩阵$W$做行变换，因此，其做变换的过程可看为对$W$左乘一个可逆矩阵$P$，使之满足$PW=P[A|E]=[PA|P]=[E|B]$，因此，有$P=B$且$PA=E$，显然，就有$BA=E$，则矩阵$B$就为矩阵$A$的逆。

下面介绍其实现过程：

1. 从上往下做行变换，使增广矩阵$W$的前一部分的方阵变为一个上三角矩阵
2. 从下往上做行变换，使增广矩阵$W$的前一部分变成一个对角矩阵
3. 每一行乘以一个系数使增广矩阵的前一部分变为单位矩阵
4. 经过变换后的增广矩阵的后一部分即为所求矩阵的逆矩阵

当然，第3步可在第1步时执行。

下面以C#语言来实现该过程，如下：

/**
	matrix为所求矩阵
	dim为矩阵的维度
*/
   public static double[, ] InverseMatrix(double[,] matrix, int dim)
   {
        readonly double eps = 1e-6;
        double[,] mat = new double[dim * 2, dim * 2];
        //构造增广矩阵W
        for(int i = 0;i < dim; i++)
        {
           for(int j = 0;j < 2 * dim; j++)
           {
                if(j < dim)
                {
                    mat[i, j] = matrix[i, j];
                }
                else
                {
                    mat[i, j] = j - dim == i ? 1 : 0;
                }
            }
        }
        for(int i = 0;i < dim; i++)
        {
            if(Math.Abs(mat[i,i]) < eps)
            //若mat[i,i]为0，则往下找一个不为0的行，加到第i行
            {
                int j;
                for ( j = i + 1; j < dim; j++)
                {
                    if (Math.Abs(mat[j, i]) > eps) break;
                }
                if (j == dim) return null;
                for(int r = i; r < 2 * dim; r++)
                {
                    mat[i, r] += mat[j, r]; 
                }
            }
            double ep = mat[i, i];
            //将mat[i,i]变为1
            for (int r = i; r < 2 * dim; r++)
            {
                mat[i, r] /= ep;
            }

            for(int j = i + 1; j < dim; j++)
            {
                double e = -1 * (mat[j, i] / mat[i, i]);
                for(int r = i; r < 2 * dim; r++)
                {
                    mat[j, r] += e * mat[i, r];
                }
            }
        }

        for(int i = dim - 1; i >= 0; i--)
        {
            for (int j = i - 1; j >= 0; j--)
            {
                double e = -1 * (mat[j, i] / mat[i, i]);
                for (int r = i; r < 2 * dim; r++)
                {
                    mat[j, r] += e * mat[i, r];
                }
            }
        }

        double[,] result = new double[dim, dim];
        for(int i = 0;i < dim; i++)
        {
            for(int r = dim; r < 2 * dim; r++)
            {
                result[i, r - dim] = mat[i, r];
            }
        }
        return result;
   }

为了验证该算法的效率，使用VS2017做了实验。针对多个维度，随机产生方阵计算其逆矩阵，并执行30次求逆过程，平均所耗时间如下表所示：

矩阵维度	10	15	30	50	80	100	200	500	1000
时间(ms)	0.10	0.13	0.76	2.74	11.96	22.76	119.58	1657.01	13712.42