C语言求矩阵的逆（高斯法）

C语言，利用高斯消去法求矩阵的逆

文章共4,721字 · 阅读需要大约16分钟

一键AI生成摘要，助你高效阅读

问答

是王久久阿

7723人浏览 · 2022-11-30 23:13:21

是王久久阿 · 2022-11-30 23:13:21 发布

初等变换法是常用的矩阵求逆方法之一

相对于伴随法，初等行变换法有着较低的时间复杂度，可以进行相对高维的矩阵运算，但同时也会损失一点点精度。

伴随法可参考之前的博客：C语言求矩阵的逆（伴随法）

LU矩阵分解法求逆可参考：LU分解法求矩阵逆

更新：

原代码在进行增广矩阵进行内存拷贝的时，使用的是n（2*col），应改为col(原矩阵的列)，此处已做更改，并在原代码中指出。
最后释放内存时需要循环释放，已做修改。
_msize函数在Mac/Linux中无法使用，这里给出了另一种解决思路，见方法二。
更改选主元时的bug。

数学原理

矩阵的初等行变换

矩阵的初等变换又分为矩阵的初等行变换和矩阵的初等列变换。矩阵的初等行变换和初等列变换统称为初等变换。另外：分块矩阵也可以定义初等变换。

定义：如果B可以由A经过一系列初等变换得到，则称矩阵A与B称为等价。

数域V上的矩阵具有以下三种行变换：

1）以V中一个非零的数，乘以矩阵的某一行

2）把矩阵的某一行的a倍加到另一行(a∈V)

3）互换矩阵中两行的位置

矩阵A通过以上三种行变换变成矩阵B，称为A→B

（注意：不要与行列式的行变换性质混淆）

利用增广矩阵求逆

如果想求矩阵A的逆，可以通过拼上一个单位E矩阵，形成A|E的增广矩阵，通过初等行变换使其成为E|B，B就是A的逆矩阵。

例如：

$gif.latex?A%3D%5Cbegin%7Bbmatrix%7D%201%20%260%20%260%20%5C%5C%201%26%201%20%260%20%5C%5C%201%261%20%261%20%5Cend%7Bbmatrix%7D$

(1)形成A|E

(2)进行初等行变换，A|E→E|B

$gif.latex?C%3D%5Cbegin%7Bbmatrix%7D%201%20%26%200%26%200%20%261%20%26%200%260%20%5C%5C%201%26%201%26%200%26%200%20%26%201%260%20%5C%5C%201%26%201%20%261%20%26%200%20%260%20%26%201%20%5Cend%7Bbmatrix%7D%5Crightarrow%20%5Cbegin%7Bbmatrix%7D%201%20%26%200%26%200%20%261%20%26%200%260%20%5C%5C%200%26%201%26%200%26%20-1%20%26%201%260%20%5C%5C%201%26%201%20%261%20%26%200%20%260%20%26%201%20%5Cend%7Bbmatrix%7D$

$gif.latex?%5Crightarrow%20%5Cbegin%7Bbmatrix%7D%201%20%26%200%26%200%20%261%20%26%200%260%20%5C%5C%200%26%201%26%200%26%20-1%20%26%201%260%20%5C%5C%200%26%201%20%261%20%26%20-1%20%260%20%26%201%20%5Cend%7Bbmatrix%7D%5Crightarrow%20%5Cbegin%7Bbmatrix%7D%201%20%26%200%26%200%20%261%20%26%200%260%20%5C%5C%200%26%201%26%200%26%20-1%20%26%201%260%20%5C%5C%200%26%200%20%261%20%26%200%20%26-1%20%26%201%20%5Cend%7Bbmatrix%7D$

逆矩阵B则是

$gif.latex?B%3D%5Cbegin%7Bbmatrix%7D%201%26%200%20%260%20%5C%5C%20-1%26%201%20%260%20%5C%5C%200%26-1%20%261%20%5Cend%7Bbmatrix%7D$

选择主元

主元就是在矩阵消去过程中，每列的要保留的非零元素，用它可以把该列其他消去。在阶梯型矩阵中，主元就是每个非零行第一个非零元素就是主元。

高斯消去法在消元过程中可能出现零主元，即 $a_{kk}^{n}=0$ ，这时消元过程将无法进行；也可能主元绝对值非常小，用它做除法将会导致舍入误差的扩散，使数值解不可靠。解决该问题的办法是避免使用绝对值过小的元素作主元。

选择主元的方法：

1）找到主对角线以下每列最大的数Max所在的行数k

2）利用初等行变换——换行变换，将k行与当前主元行互换（记录总共换行次数n）

3）以当前主元行为基，利用初等行变换——消法变换，将主对角线下方消0

4）行列式每次换行需变号，行列式最后的符号为 $-1^{n}$

5）每次进行高斯消去前都必须选择主元，计算n维的行列式，则需要进行n-1次主元选择

程序设计

法一

整体流程：

1）判断传入指针是否为空，判断矩阵是否为方针

2）为增广矩阵、输出的逆矩阵开辟空间，并初始化为0

3）将原矩阵中的数据拷贝到增广矩阵中，并将增广矩阵右侧化为单位阵

4）逐列开始，选择主元，将其他位置化为0，每列阶梯位置化为1

6）将增广矩阵右侧的逆矩阵拷贝到输入矩阵中，并释放增广矩阵内存

整体代码

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>

double** Matrix_inver(double** src)
{
	//step 1
	//判断指针是否为空
	if (src == NULL)exit(-1);
	int i, j, k, row, col, n，principal;
	double** res, ** res2,tmp;//res为增广矩阵，res为输出的逆矩阵
    double Max;
	//判断矩阵维数
	row = (double)_msize(src) / (double)sizeof(double*);
	col = (double)_msize(*src) / (double)sizeof(double);
	if (row != col)exit(-1);
	//step 2
	res = (double**)malloc(sizeof(double*) * row);
	res2 = (double**)malloc(sizeof(double*) * row);
	n = 2 * row;
	for (i = 0; i < row; i++)
	{
		res[i] = (double*)malloc(sizeof(double) * n);
		res2[i] = (double*)malloc(sizeof(double) * col);
		memset(res[i], 0, sizeof(res[0][0]) * n);//初始化
		memset(res2[i], 0, sizeof(res2[0][0]) * col);
	}
	//step 3
	//进行数据拷贝
	for (i = 0; i < row; i++)
	{
    //此处源代码中的n已改为col，当然方阵的row和col相等，这里用col为了整体逻辑更加清晰
		memcpy(res[i], src[i], sizeof(res[0][0]) * col);
	}
	//将增广矩阵右侧变为单位阵
	for (i = 0; i < row; i++)
	{
		for (j = col; j < n; j++)
		{
			if (i == (j - row))
				res[i][j] = 1.0;
		}
	}
	
	for (j = 0; j < col; j++)
	{
       //step 4
	   //整理增广矩阵，选主元
		principal = j;
        Max = fabs(res.[principal][j]); // 用绝对值比较
        // 默认第一行的数最大
        // 主元只选主对角线下方
        for (i = j; i < row; i++)
        {
            if (fabs(res.[i][j]) > Max)
            {
                principal = i;
                Max = fabs(res.[i][j]);
            }
        }
        if (j != principal)
        {
            for (k = 0; k < n; k++)
            {
                tmp = res.[principal][k];
                res.[principal][k] = res.[j][k];
                res.[j][k] = tmp;
            }
        }
        //step 5
		//将每列其他元素化0
		for (i = 0; i < row; i++)
		{
			if (i == j || res[i][j] == 0)continue;
			double b = res[i][j] / res[j][j];
			for (k = 0; k < n; k++)
			{
				res[i][k] += b * res[j][k] * (-1);
			}
		}
		//阶梯处化成1
		double a = 1.0 / res[j][j];
		for (i = 0; i < n; i++)
		{
			res[j][i] *= a;
		}
	}
	//step 6
	//将逆矩阵部分拷贝到res2中
	for (i = 0; i < row; i++)
	{
		memcpy(res2[i], res[i] + row, sizeof(res[0][0]) * row);
	}
	//必须释放res内存！
    for(i = 0; i < row; i++)
    {
        free(res[i]);
    }
    free(res);
	return res2;
}

上述代码中：

函数传入的参数需是以malloc开辟的动态矩阵，如固定为二维数组，需自行更改

step 1：

row = (double)_msize(src) / (double)sizeof(double*);
col = (double)_msize(*src) / (double)sizeof(double);
为判断矩阵的维数；
其中_msize为库函数，需要包含头文件<malloc.h>
返回指针指向内存的小心，单位为字节。

step 2：

增广矩阵的行数不变，列数是原来的2倍；
对照该方法开辟矩阵空间，便可以理解step1中如何计算矩阵维数。

step 3：

memset和memcpy都需要包含头文件<string.h>
需要注意memset是以字节对内存进行赋值的，对非字符（char类型）的类型进行赋值时，只使用赋0或-1。

step 4：

选择主元：该步骤具有两层意义：
1）相对程度上提高计算精度
2）保证阶梯处的数都非0

step 5：

j代表从每列开始进行消0，i和k代表行和列；
不要混淆j和k：j代表消零的大循环，k代表赋值的小循环。

step 6：

res[i] + row代表只拷贝增广矩阵中逆矩阵部分，函数结束前必须释放增广矩阵的内存，否则将造成内存泄漏。

关于上述函数不太理解的地方可以参考之前发的判断矩阵维数和伴随法求逆矩阵

C语言动态内存部分可参考C语言动态内存管理

测试

相关测试函数（创建矩阵，初始化矩阵，打印矩阵）

double** MakeMat(int n)
{
	int i = 0;
	if (n <= 0)exit(-1);
	double** res = (double**)malloc(sizeof(double*) * n);
	if (res == NULL)exit(-1);
	for (i = 0; i < n; i++)
	{
		res[i] = (double*)malloc(sizeof(double) * n);
	}
	return res;
}

void InitMat(double** src)
{
	if (src == NULL)exit(-1);
	int i, j, n;
	n = (double)_msize(src) / (double)sizeof(double*);
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			src[i][j] = pow(i,j);
		}
	}
}

void print(double** src)
{
	if (src == NULL)exit(-1);
	putchar('\n');
	int i, j, row,col;
	row = (double)_msize(src) / (double)sizeof(double*);
	col = (double)_msize(*src) / (double)sizeof(double);
	for (i = 0; i < row; i++)
	{
		for (j = 0; j < col; j++)
		{
			printf("%9.4lf", src[i][j]);
		}
		putchar('\n');
	}
}

主函数测试：

int main()
{
	int n = 5;
	double** arr = MakeMat(n);
	InitMat(arr);
	double** res = Matrix_inver(arr);
	printf("原矩阵:>");
	print(arr);
	printf("逆矩阵:>");
	print(res);
	return 0;
}

打印结果：

如果矩阵接近奇异值，计算结果会不准确；如果矩阵秩亏，会输出无效值（nan）。

法二

整体代码

由于无法利用_msize函数直接调取内存，又不想在计算时手动输入矩阵的维数，不妨转换一下思路：在创建矩阵时就把矩阵的行和列的信息保存下来，调用矩阵时直接读取矩阵信息即可。

于是需要创建一个矩阵的结构体：

typedef struct Matrix
{
    int row;
    int col;
    double **data;
} Matrix, *pMatrix;

其中row、col保存矩阵的行、列，data指针与原来的方法一样，通过row和col为矩阵开辟内存，因此开辟矩阵的函数如下：

Matrix MakeMatrix(int row, int col)
{
	int i = 0;
	Matrix arr = {0};
	arr.row = row;
	arr.col = col;
	arr.data = (double **)malloc(sizeof(double *) * arr.row);
	if (arr.data == NULL)
		exit(-1);
	for (i = 0; i < arr.row; i++)
	{
		arr.data[i] = (double *)malloc(sizeof(double) * arr.col);
		memset(arr.data[i], 0, sizeof(double) * arr.col);
	}
	return arr;
}

此外，由于生成矩阵时是循环开辟内存的，因此在销毁矩阵时依然需要循环销毁，销毁矩阵内存的函数如下：

void free_Matrix(Matrix arr)
{
	int i;
	for (i = 0; i < arr.row; i++)
	{
		free(arr.data[i]);
	}
    free(arr.data);
}

接下来只需要对求逆函数略作更改，将参数与返回值改为结构体形式即可：

Matrix Matrix_Guass_inver(Matrix src)
{
    // step 1
    // 判断指针是否为空,判断矩阵维数
    assert(src.data);
    assert(src.row == src.col);
    // step 2
    int i, j, k, n, principal;
    double Max, tmp;
    n = 2 * src.row;
    // res为增广矩阵，res为输出的逆矩阵
    Matrix res = MakeMatrix(src.row, n);
    Matrix res2 = MakeMatrix(src.row, src.col);
    // step 3
    // 进行数据拷贝
    for (i = 0; i < src.row; i++)
    {
        memcpy(res.data[i], src.data[i], sizeof(res.data[0][0]) * src.col);
    }
    // 将增广矩阵右侧变为单位阵
    for (i = 0; i < src.row; i++)
    {
        for (j = src.col; j < n; j++)
        {
            if (i == (j - src.row))
                res.data[i][j] = 1.0;
        }
    }
    for (j = 0; j < src.col; j++)
    {
        // step 4
        // 整理增广矩阵，选主元
        principal = j;
        Max = fabs(res.data[principal][j]); // 用绝对值比较
        // 默认第一行的数最大
        // 主元只选主对角线下方
        for (i = j; i < src.row; i++)
        {
            if (fabs(res.data[i][j]) > Max)
            {
                principal = i;
                Max = fabs(res.data[i][j]);
            }
        }
        if (j != principal)
        {
            for (k = 0; k < n; k++)
            {
                tmp = res.data[principal][k];
                res.data[principal][k] = res.data[j][k];
                res.data[j][k] = tmp;
            }
        }
        // step 5
        // 将每列其他元素化0
        for (i = 0; i < src.row; i++)
        {
            if (i == j || res.data[i][j] == 0)
                continue;
            double b = res.data[i][j] / res.data[j][j];
            for (k = 0; k < n; k++)
            {
                res.data[i][k] += b * res.data[j][k] * (-1);
            }
        }
        // 阶梯处化成1
        double a = 1.0 / res.data[j][j];
        for (i = 0; i < n; i++)
        {
            res.data[j][i] *= a;
        }
    }
    // step 6
    // 将逆矩阵部分拷贝到res2中
    for (i = 0; i < src.row; i++)
    {
        memcpy(res2.data[i], res.data[i] + src.row, sizeof(res.data[0][0]) * src.row);
    }
    // 释放增广矩阵内存
    free_Matrix(res);
    return res2;
}

测试

void test()
{
    Matrix a = MakeMatrix(3, 3);
    a.data[0][0] = 1;
    a.data[1][1] = 2;
    a.data[0][1] = 1;
    a.data[2][0] = 1;
    a.data[2][2] = 2;
    Matrix res = Matrix_Guass_inver(a);
}