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

数学基础之矩阵系列

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

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

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

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

下面介绍其实现过程:

  1. 从上往下做行变换,使增广矩阵 W W W的前一部分的方阵变为一个上三角矩阵
  2. 从下往上做行变换,使增广矩阵 W W 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次求逆过程,平均所耗时间如下表所示:

矩阵维度10153050801002005001000
时间(ms)0.100.130.762.7411.9622.76119.581657.0113712.42
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值