在一些需要用单片机、DSP、ARM等系统中用C语言实现高维滤波器的场合，例如卡尔曼滤波器，我们会经常遇到用C语言求解方阵的逆的情况，当阶数为2阶或者3阶的时候我们可以用公式法直接求解，但是当阶数一旦达到4阶或者4阶以上时利用公式求解将会非常的麻烦，而且极易出错，这时我们需要寻求一种可以求解任意阶数的算法。由线性代数的知识可以知道求任意阶数矩阵逆矩阵有2种算法，一种是初等行（列）变换，一种是伴随矩阵法。其中初等行（列）变换由于需要进行行（列）的加减或交换，以及乘上（除以）非零的数这3种步骤，且顺序不固定，因此不适合用C程序进行实现。而伴随矩阵的方法求解步骤固定因此适合用C程序实现。

算法描述

根据线性代数的知识，我们可以知道对于N阶方阵 $A=\begin{bmatrix} a_{11} &a_{12} &... &a_{1N} \\ a_{21}&a_{22} &... &a_{2N} \\ ...&... &... &... \\ a_{N1}&a_{N2} &... &a_{NN} \end{bmatrix}$ ，它的伴随矩阵为 $A^{*}=\begin{bmatrix} A_{11} &A_{21} &... &A_{N1} \\ A_{12}&A_{22} &... &A_{N2} \\ ...&... &... &... \\ A_{1N}&A_{2N} &... &A_{NN} \end{bmatrix}$ ，其中 $A_{ij}=(-1)^{i+j}M_{ij}$ 称为元素 $a_{ij}$ 的代数余子式。则A的逆矩阵为 $A^{-1}=A^{*}/\left | A \right |$ ，其中 $\left | A \right |$ 为矩阵A的行列式。

C语言实现代码

/***************************************************************************************
 * 函数名称：Matrix_inv
 * 输入参数：x输入的NxN矩阵、y输出的NxN矩阵、N矩阵的阶数
 * 输出参数：1：表示输入矩阵有逆、0：表示输入矩阵的逆不存在
 * 实现功能：通过求伴随矩阵的方法求NxN矩阵的逆
 * 注意事项：需要的动态内存空间最少为：(2*2 + 3*3 + 4*4 + ... + (N-1)*(N-1))*4 = ((N-1)N(2N-1)/6 - 1)*4
 ****************************************************************************************/
int Matrix_inv(float *x, float *y, int N)
{
	float det_x;
	float recip_det_x;
	float coff;
	float *sub_x;
	int i, j, k, sub_N;

	det_x = Matrix_det(x, N);   //求行列式的值
	if(det_x == 0)
	{
		return 0;     //表示矩阵的逆不存在
	}
	else
	{
		recip_det_x = 1/det_x;
		coff = 1;
		sub_N = N-1;
		sub_x = (float *)malloc(4*sub_N*sub_N);   //初始值随机，动态分配N-1阶的数组
		for(i=0; i<N; i++)
		{
			for(j=0; j<N; j++)
			{
				for(k = 0; k< sub_N; k++)    //提取去除第i行，第j列的数组
				{
					if(i<=k)
					{
						if(j == 0)                      { memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k+1)*N+1]),      4*sub_N); }            //第0列时
						else if(j == sub_N) { memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k+1)*N]),            4*sub_N); }            //最后第N-1列时
						else                                {  memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k+1)*N]),            4*j);
															         memcpy((void *)(&sub_x[k*sub_N+j]),(void *)(&x[(k+1)*N+j+1]), 4*(sub_N-j)); }
					}
					else
					{
						if(j == 0)                      { memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k)*N+1]),      4*sub_N); }            //第0列时
						else if(j == sub_N) { memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k)*N]),            4*sub_N); }            //最后第N-1列时
						else                                {  memcpy((void *)(&sub_x[k*sub_N]),    (void *)(&x[(k)*N]),            4*j);
															         memcpy((void *)(&sub_x[k*sub_N+j]),(void *)(&x[(k)*N+j+1]), 4*(sub_N-j)); }
					}
				}
				y[j*N+i] = coff*Matrix_det(sub_x, sub_N)*recip_det_x;     //第i行j列的代数余子式，并除上行列式的值
				coff = -coff;
			}
		}
		free((void *)sub_x);      //使用完后释放内存
	}

	return 1;     //表示矩阵的逆存在
}

从上述代码里我们可以看到，程序的实现步骤主要是：1、首先利用函数Matrix_det（关于行列式的C语言实现请参考博客任意阶数实数方阵的行列式的值的C语言实现详解。里面有关于矩阵行列式求解的具体C语言实现过程。）计算矩阵A的行列式，2、然后判断行列式的值是否为0，若为0则返回0表示矩阵的逆不存在。3、若非零，则逐个求解各元素的代数余子式并除以A的行列式，一共3大步骤。在实现的过程中，使用到了malloc（包含在头文件 <stdlib.h>中）函数动态分配内存，一共4*(N-1)*(N-1)个字节，用于存储N-1阶余子式的各元素，在Matrix_det中也使用了malloc函数进行内存分配，其中首次Matrix_det用求解A的行列式的时候占用的内存大小为 $((N-1)N(2N-1)/6-1)\times 4$ ，但是用完后会被释放掉，与后面的不存在同时使用的情况。因此只需要考虑后面的Matrix_det使用的字节数量，它的大小为 $((N-2)(N-1)(2N-3)/6-1)\times 4$ 再加上与它同时使用用于存储N-1阶余子式元素的4*(N-1)*(N-1)个字节一共还是 $((N-1)N(2N-1)/6-1)\times 4$ 个字节，因此执行此函数需要的总的内存空间最低为 $((N-1)N(2N-1)/6-1)\times 4$ 。而这些内存都是使用的堆(heap)空间，因此要保证系统的稳定运行heap的大小不能低于该值。程序中使用了memcpy（包含在头文件 <string.h>中）函数进行数据搬移，主要是为了提高程序执行效率，它的效率比直接的单个数的赋值要高。

原文链接：https://blog.csdn.net/HJ199404182515/article/details/111881189