• C语言矩阵求逆

    之前介绍了C语言用代数余子式求行列式

    本次开始介绍如何用公式法对矩阵求逆，并用C语言将其实现。

    之前程序有点小bug，已于2022年11月29日修改。

    更新：

            伴随法只适合求低阶矩阵的逆，对于相对高阶（20维以上）对矩阵求逆用高斯法求解效率更高，此外本文中使用了_msize函数用于判断内存维数，但该函数只适合winodows系统，Linux和Mac系统无法使用（笔者也是在用了Mac系统后才发现），对于上述两个问题，您应该可以在：C语言求矩阵的逆（高斯法）得到满意的答案。

            如果矩阵接近奇异值，求逆的数值将不稳定，那么使用C语言LU分解法求逆将会得到更好的效果。

    目录

    数学原理

    矩阵求逆的公式

    数乘矩阵 

    程序设计

    整体代码

    求行列式的值（Det、Cof）

    求伴随矩阵（FindCof）

    求逆矩阵的主函数（matrix_inver）

    测试

    数学原理

    矩阵求逆的方法有很多种，本次主要介绍如何利用公式法求解。

    矩阵求逆的公式

    A^{-1}=\frac{1}{\left | A \right |}A^{*}

    其中，\left | A \right |为A的行列式的值，A^{*}为矩阵A的伴随矩阵。

    伴随矩阵的表达方式为：

    A_{ij}^{*}=\begin{bmatrix} A_{11} &A_{12} &... &A_{1j} \\ A_{21}&A_{22} & ... & ...\\ ... &... & ...& ..\\ A_{i1} &A_{i2} &.. & A_{ij} \end{bmatrix}^{T}

    其中，A_{ij}为代数余子式，代数余子式A_{ij}与余子式M_{ij}的关系为：

    A_{ij}=\left ( -1 \right )^{i+j}M_{ij}

    M_{ij}是矩阵A去掉i行j列，剩下元素重新组成的新矩阵 。

    数乘矩阵 

    假设矩阵A为：

    A=\begin{bmatrix} a_{11} &a_{12} &a_{13} \\ a_{21}& a_{22} & a_{23}\\ a_{31}& a_{32} & a_{33} \end{bmatrix}

    如果一个数b乘以一个矩阵A，那么A中的每一个元素都要乘以b，其表示如下： 

    bA=\begin{bmatrix} ba_{11} &ba_{12} &ba_{13} \\ ba_{21}& ba_{22} & ba_{23}\\ ba_{31}& ba_{32} & ba_{33} \end{bmatrix}

    程序设计

    为了后面方便调试，先利用malloc函数创建一个矩阵，并为其赋值。关于malloc函数有不理解的可以参考之前写的C语言动态内存管理。

    malloc使用的基本方式：

    • void* malloc (size_t size);
    • void free (void* ptr);
    • size是指定的开辟内存的大小，单位是字节
    • size_t的无符号整型则限制程序员误操作开辟负字节的空间
    • 如果开辟成功，malloc会返回一个void*类型的指针
    • 如果开辟失败，则返回的是空指针，所以在malloc之后需要对指针进行检查
    • 当malloc的东西不再使用时，需要free对其进行释放，否则会造成内存泄漏
    • malloc和free均需要包含头文件
    1. //创建n维矩阵空间，并初始化
    2. double** test1(int n)
    3. {
    4. 	double** arr = (double**)malloc(sizeof(double*) * n);
    5. 	int i, j;
    6. 	if (arr != NULL)
    7. 	{
    8. 		for (i = 0; i < n; i++)
    9. 		{
    10. 			arr[i] = (double*)malloc(sizeof(double) * n);
    11. 		}
    12. 		//为矩阵赋值
    13. 		if (*arr != NULL)
    14. 		{
    15. 			for (i = 0; i < n; i++)
    16. 			{
    17. 				for (j = 0; j < n; j++)
    18. 				{
    19. 					arr[i][j] = pow(i, j);
    20. 				}
    21. 			}
    22. 		}
    23. 	}
    24. 	return arr;
    25. }

    其中：

    • double** arr = (double**)malloc(sizeof(double*) * n);表示开辟n个double*类型的空间，将其地址赋给arr
    • 然后再开辟n个double类型的空间，将其地址赋给double**中每个double*
    • 配合下图更好理解

    整体代码

    1. #include
    2. #include
    3. #include
    4. #include
    5. #include
    6. #include
    7.  
    8. #define MatMax 20
    9.  
    10. //函数声明
    11. double Det(const double arr[MatMax][MatMax], int n);
    12. double Cof(const double arr[MatMax][MatMax], int i, int n);
    13. void FindCof(double arr[MatMax][MatMax], double arr2[MatMax][MatMax], int i, int j, int n);
    14. double** matrix_inver(double** arr);
    15.  
    16. //计算行列式
    17. double Det(const double arr[MatMax][MatMax], int n)
    18. {
    19. 	assert(n > 0);
    20. 	double sum = 0;
    21. 	int i = 0;
    22. 	if (n == 1)//1阶行列式直接得出结果
    23. 	{
    24. 		sum = arr[0][0];
    25. 	}
    26. 	else if (n == 2)
    27. 	{
    28. 		sum = arr[0][0] * arr[1][1] - arr[0][1] * arr[1][0];//杀戮法求解
    29. 	}
    30. 	else if (n == 3)
    31. 	{
    32. 		sum = arr[0][0] * arr[1][1] * arr[2][2]
    33. 			+ arr[0][1] * arr[1][2] * arr[2][0]
    34. 			+ arr[1][0] * arr[2][1] * arr[0][2]
    35. 			- arr[0][2] * arr[1][1] * arr[2][0]
    36. 			- arr[0][1] * arr[1][0] * arr[2][2]
    37. 			- arr[1][2] * arr[2][1] * arr[0][0];//划线法求解
    38. 	}
    39. 	else
    40. 	{
    41. 		for (i = 0; i < n; i++)//按第一行展开
    42. 		{
    43. 			if (arr[0][i] != 0)//展开项不为0才计算
    44. 			{
    45. 				sum += ((int)pow(-1, i + 0)) * arr[0][i] * (Cof(arr, i, n));//2阶以上继续递归		
    46. 			}
    47. 			else
    48. 				sum += 0;//展开项为0
    49. 		}
    50. 	}
    51. 	return sum;
    52. }
    53. //找到余子式
    54. double Cof(const double arr[MatMax][MatMax], int i, int n)
    55. {
    56. 	assert(n > 0);
    57. 	int k = 0;
    58. 	int j = 0;
    59. 	double arr2[MatMax][MatMax] = { 0 };
    60. 	for (k = 0; k < n - 1; k++)//去除0行i列，剩下的组成新的矩阵
    61. 	{
    62. 		for (j = 0; j < n - 1; j++)
    63. 		{
    64. 			if (j < i)
    65. 			{
    66. 				arr2[k][j] = arr[k + 1][j];
    67. 			}
    68. 			else
    69. 			{
    70. 				arr2[k][j] = arr[k + 1][j + 1];
    71. 			}
    72. 		}
    73. 	}
    74. 	return Det(arr2, n - 1);
    75. }
    76. //找到去掉i行j列的余子式
    77. void FindCof(double arr[MatMax][MatMax], double arr2[MatMax][MatMax], int i, int j, int n)
    78. {
    79. 	int m = 0;
    80. 	int k = 0;
    81. 	for (m = 0; m < n - 1; m++)
    82. 	{
    83. 		for (k = 0; k < n - 1; k++)
    84. 		{
    85. 			if (k < j)
    86. 			{
    87. 				if (m < i)
    88. 				{
    89. 					arr2[m][k] = arr[m][k];
    90. 				}
    91. 				else
    92. 				{
    93. 					arr2[m][k] = arr[m + 1][k];
    94. 				}
    95. 			}
    96. 			else
    97. 			{
    98. 				if (m < i)
    99. 				{
    100. 					arr2[m][k] = arr[m][k + 1];
    101. 				}
    102. 				else
    103. 				{
    104. 					arr2[m][k] = arr[m + 1][k + 1];
    105. 				}
    106. 			}
    107. 		}
    108. 	}
    109. }
    110.  
    111.  
    112. //计算逆的主函数
    113. double** matrix_inver(double** arr)
    114. {
    115. 	int i, j, n;
    116. 	double** res=NULL;
    117. 	n = (int)_msize(arr) / (int)sizeof(double*);
    118. 	res = (double**)malloc(sizeof(double*) * n);
    119. 	if (res == NULL)exit(-1);
    120. 	for (i = 0; i < n; i++)
    121. 	{
    122. 		res[i] = (double*)malloc(sizeof(double) * n);
    123. 		memset(res[i], 0, sizeof(double) * n);
    124. 	}
    125. 	double tmp[MatMax][MatMax] = { 0 };
    126. 	//保护arr，将arr指向内存的数据拷贝到tmp二维数组中
    127. 	for (i = 0; i < n; i++)
    128. 	{
    129. 		memcpy(tmp[i],arr[i],sizeof(double) * n);
    130. 	}
    131. 	double a = 1.0 / (Det(tmp, n));
    132. 	for (i = 0; i < n; i++)
    133. 	{
    134. 		for (j = 0; j < n; j++)
    135. 		{
    136. 			double tmp2[MatMax][MatMax] = { 0 };
    137. 			FindCof(tmp, tmp2, j, i, n);//求转置后的伴随
    138. 			double b = pow(-1, i + j) * Det(tmp2, n - 1);
    139. 			res[i][j] = a * b;
    140. 		}
    141. 	}
    142. 	return res;
    143. }

    接下来针对上述代码进行讲解

    求行列式的值（Det、Cof）

    通过观察发现，利用公式法求矩阵的逆离不开求行列式的值。我们创建两个函数Det和Cof2

    利用代数余子式的方法求解行列式：

    Det：求解行列式的主函数，判断行列式的维数，如果维数大于3则进行代数余子式计算，按第一行展开，将其传给Cof函数；如果小于等于3则直接求解，避免过度递归。

    Cof：求余子式，将新的余子式传给Det函数。

    行列式的计算在这篇博客中C语言计算行列式

    求伴随矩阵（FindCof）

    该函数与cof函数类似，都是求余子式的，区别在于：

    • FindCoi是找到去掉第i行第j列的余子式；Cof函数固定从第1行第j列
    • FindCof函数不参与递归，只要找到余子式即可

    其中：

    • if部分是用来判断组成新的行列式的位置

    求逆矩阵的主函数（matrix_inver）

    对照之前的公式，对n维的方阵，其伴随矩阵需要求解n*n个代数余子式，利用两个for循环逐行逐列进行处理。

    主要思想：

    • 计算行列式的值a，计算1/a
    • 从i行j列开始，求代数余子式，并将求得的余子式部分传给Det函数，计算其行列式的值b
    • 计算a*b，并赋值

    主函数中：

    •  n = (int)_msize(arr) / (int)sizeof(double*);求矩阵的维数，_msize是库函数，需要包含头文件：求其指针指向地址的内存大小，单位是字节；对照res开辟内存的过程，即可明白其计算原理。
    • memcpy(tmp[i],arr[i],sizeof(double) * n);进行内存拷贝，将动态内存中数据拷贝到二维数组中。有两层意义：一：防止递归中造成内存泄漏；二：保护原矩阵

    测试

    1. int main()
    2. {
    3. 	int n = 5;
    4. 	double** arr = test1(n);
    5. 	printf("原矩阵:>\n");
    6. 	print(arr);
    7. 	double** res = matrix_inver(arr);
    8. 	printf("逆矩阵:>\n");
    9. 	print(res);
    10. 	return 0;
    11. }

    其中： 

    • test1是上述创建矩阵的函数
    • print是用来打印矩阵的函数，代码如下：
    1. //打印矩阵
    2. void print(double** arr)
    3. {
    4. 	putchar('\n');
    5. 	int i, j, row, col;
    6. 	row = (int)_msize(arr) / (int)sizeof(double*);//判断行数
    7. 	col = (int)_msize(*arr) / (int)sizeof(double);//判断列数
    8. 	for (i = 0; i < row; i++)
    9. 	{
    10. 		for (j = 0; j < col; j++)
    11. 		{
    12. 			printf("%10.5lf ", arr[i][j]);
    13. 		}
    14. 		putchar('\n');
    15. 	}
    16. 	putchar('\n');
    17. }

    测试结果如下： 

    与matlab计算的结果进行比对：

     如果矩阵治亏，则会输出：

    因为公式法涉及到递归，因此无法计算高维逆矩阵，想要保证精度和效率，可以尝试利用LU分解法进行矩阵求逆，这也是目前大多数计算机在处理高维矩阵时的策略。

  • 相关阅读:
    位、比特、字节、字、帧等概念关系的理解
    Apache Hudi技术与架构-1
    vue2和vue3响应式的理解及如何进行依赖收集
    CMake教程-第 3 步：添加库的使用要求
    回归预测 | MATLAB实现BP神经网络多输入多输出回归预测
    基于openwrt交叉编译opencv4.9.0版本
    SpringCloud学习（四）——Nacos注册中心
    扩散模型（Diffusion Model）简介
    （二十六）大数据实战——kafka集群之Kraft模式安装与部署
    HMM隐马尔可夫模型
  • 原文地址：https://blog.csdn.net/why1472587/article/details/127558213