变量循环重新标号法求对称正定矩阵逆矩阵

一 算法原理

1. 对称矩阵特征值算法

   雅可比方法用于求解实对称矩阵的特征值和特征向量,对于实对称矩阵$A$,必有正交矩阵$U$,使得$U^{T}AU=D$.$D$是一个对角阵,主对角线的元素是矩阵$A$的特征值,正交矩阵$U$的每一列对应于属于矩阵$D$的主对角线对应元素的特征向量.
   雅可比方法用平面旋转对矩阵$A$做相似变换,化$A$为对角阵,从而求解出特征值和特征向量.旋转矩阵$U_p{}_q$,是一个单位阵在第$p$行,第$p$列,第$q$行,第$q$列,元素为$cos\phi$,第$p$行第$q$列为$-sin\phi$,第$q$行第$p$列为$sin\phi$.对于这样的平面旋转矩阵,不难验证其是一种正交矩阵.因此对于向量$x$,$U_p{}_{q}x$等同于把第$p$个坐标轴和第$q$个坐标轴共同所确定的平面旋转了$\phi$度.记矩阵$A_1=U_p{{}_q}^{T}A{U}_p{}_q$.因为旋转矩阵是正交阵,因此实际上矩阵$A_1$与矩阵$A$是相似的,因此其特征值是相同的.
   设矩阵$A_1$第$i$行,第$j$列的元素为$b_i{}_j$,矩阵$A$的第$i$行,第$j$列的元素为$a_i{}_j$($i=0,1,2,...,n-1,j=0,1,2,...,n-1$).式(1-1-1)给出了两矩阵元素之间的运算关系.
   $$\{\begin{array}{l}{b}_p{}_p=a_p{}_{p}co{s}^2\phi +a_q{}_{q}si{n}^2\phi +2a_p{}_{q}cos\phi sin\phi \\ b_q{}_q=a_p{}_{p}si{n}^2\phi +a_q{}_{q}co{s}^2\phi -2a_p{}_{q}cos\phi sin\phi \\ b_p{}_q=b_q{}_p=\frac{1}{2}(a_q{}_q-a_p{}_p)sin2\phi +a_p{}_{q}cos2\phi \\ b_p{}_i=a_p{}_{i}cos\phi +a_q{}_{i}sin\phi ,(i\ne p,q)\\ b_q{}_i=-a_p{}_{i}sin\phi +a_q{}_{i}cos\phi ,(i\ne p,q)\\ b_j{}_p=a_j{}_{p}cos\phi +a_j{}_{q}sin\phi ,(j\ne p,q)\\ b_j{}_q=-a_j{}_{q}sin\phi +a_j{}_{q}cos\phi ,(j\ne p,q)\\ b_i{}_j=b_j{}_i=a_i{}_j,i\ne p,q;j\ne p,q\end{array}\text{\hspace{1em}(1-1-1)}$$
   其中有两点需要说明:(1) $p$,$q$分别是前一次的迭代矩阵的非主对角线上绝对值最大元素的行列号
   (2)$\phi$是旋转角度,可以由式(1-1-2)确定
   $$tan2\phi =\frac{-2a_p{}_q}{a_q{}_q-a_p{}_p}\text{\hspace{1em}(1-1-2)}$$

   归纳得到雅可比方法求解矩阵特征值和特征向量的具体步骤如下:
   (1).初始化特征向量为对角阵V,主对角线元素为1,其他元素为0.
   (2).在$A$的非主对角线的元素中,找到绝对值最大元素$a_p{}_q$.
   (3).用式(1-1-2)计算出旋转矩阵
   (4).计算矩阵$A1$,用当前的矩阵$V$乘旋转矩阵得到当前的特征矩阵$V$.
   (5).若当前迭代的矩阵$A$的非主对角线元素最大值小于给定阈值,停止计算,否则执行上述过程.停止计算时,特征值为矩阵$A$的主对角线元素,特征矩阵为矩阵$V$.

2. 对称正定矩阵求逆算法

   首先,我们先考虑这样的一个线性方程组(式1-2-1),我们现在想做的事情是把每一个$x$写成是$y$的线性组合,这样就等同于将方程组的系数矩阵进行了求逆运算.现在介绍一种求解对称的正定矩阵逆矩阵方法,考虑该方程组的第一个方程,将第一个方程改写成为式(1-2-2)的形式,即将$x_1$表达为其余项的线性组合.现在我们将$x_1$带入到式(1-2-1)中,得到了式(1-2-3).如果我们对每一个$x$都这样计算的话,这和普通求逆的方法别无它异,这里我们使用一种新的方法,变量循环重新标号法,这种方法的核心就是让带入变量的操作只发生在$x_1$中.在式(1-2-3)中,我们对所有的变量执行下述操作(1-2-4),称执行(1-2-4)的操作为一次变量更替,容易发现,连续做$n$次变量更替后方程组会回到本身.算法执行的步骤是这样的:我们先对原方程(1-2-1)进行(1-2-2)的那样操作,这一步操作记作操作$1$,然后将$x_1$代入到其他方程中记作操作$2$.然后进行变量更替,此时继续上述操作,共执行$n$(系数矩阵维度)次,可以得到系数矩阵对应的逆矩阵.
   在下一节中的函数$fun1$实现的就是**操作$1$的功能,而函数$n2$实现的是操作$2$**的功能.
   $$\{\begin{array}{l}{y}_1=C_1{}_1{x}_1+C_1{}_2{x}_2+C_1{}_3{x}_3+\cdots +C_1{}_n{x}_n\\ y_2=C_2{}_1{x}_1+C_2{}_2{x}_2+C_2{}_3{x}_3+\cdots +C_2{}_n{x}_n\\ y_3=C_3{}_1{x}_1+C_3{}_2{x}_2+C_3{}_3{x}_3+\cdots +C_3{}_n{x}_n\\ ⋮\\ y_n=C_4{}_1{x}_1+C_4{}_2{x}_2+C_4{}_3{x}_3+\cdots +C_4{}_n{x}_n\end{array}\text{\hspace{1em}(1-2-1)}$$

$$x_1=C_1{{}_2}^{\prime }{x}_2+C_1{{}_3}^{\prime }{x}_3+C_1{{}_4}^{\prime }{x}_4+\cdots +C_1{{}_n}^{\prime }{y}_1\text{\hspace{1em}(1-2-2)}$$

---

二 C语言实现

//本程序的功能是输入一个矩阵,如果它是对称正定的,则用变量循环标号法计算其逆矩阵,若不是,则退出.
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void Show_Matrix(float **array,int n);//输出矩阵
int Check(float **array,int n);//返回1代表是实对称矩阵，返回其他代表不是实对称矩阵
void fun1(float **array,int n,int k);//x1=x1(x2,x3,.....,y1)
void fun2(float **array,int n,int k);//将x1带入到其他方程中
int Matrix_Free(float **tmp, int m, int n);//释放申请的内存空间
float** Matrix_Jac_Eig(float **array, int n, float *eig);
int main(void)
{
	int n;
	printf("请输入矩阵维度:\n");
	scanf("%d", &n);
	float **array = (float **)malloc(n * sizeof(float *));
	for (int i = 0; i < n; i++)
	{
		array[i] = (float *)malloc(n * sizeof(float));
	}
	printf("请输入矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			scanf("%f", &array[i][j]);
		}
	}
	if (Check(array,n) == 1)//判断是否为对称正定矩阵
	{
		for (int i = 0; i < n; i++)
		{
			fun1(array,n,i);
			fun2(array,n,i);
		}
		Show_Matrix(array,n);
	}
	else
	{
		printf("您输入的矩阵不是对称正定矩阵\n");
		return -1;
	}
	return 0;
}
void fun1(float **array,int n,int k)//操作1
{
	float temp = array[k][0];
	for (int j = 0; j < 2; j++)
	{
		array[k][j] = -array[k][j + 1] / temp;
	}
	array[k][n - 1] = 1 / temp;
}
void fun2(float **array,int n,int k)//操作2
{
	int i,j;
	float temp;
	for (i = 0; i < n; i++)
	{
		if (i == k)
			continue;
		else
		{
			temp = array[i][0];
			for (j = 0; j < 2; j++)
			{
				array[i][j] = array[i][j + 1] + array[k][j] * temp;
			}
			array[i][2] = temp*array[k][j];
		}
	}
}
void Show_Matrix(float **array,int n)
{
	int i, j;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			printf("%f ", array[i][j]);
		}
		printf("\n");
	}
}
int Check(float **array,int n)
{
	//检查是否为对称正定矩阵,如果返回的是1,则说明是对称正定矩阵
	//返回其他则说明不是对称正定矩阵.通过判断其特征值来判断是否为正定矩阵.
	//如果有一个特征值是负数,则离开返回.
	int flag = 1;
	int i = 0;
	float *eig = (float *)malloc(n * sizeof(float));;
	float **special_vector = Matrix_Jac_Eig(array, n, eig);
	for (i = 0; i < n; i++)
	{
		if (eig[i] < 0)
		{
			flag = 0;
			break;
		}
	}
	Matrix_Free(special_vector, n, n);
	free(eig);
	return flag;
}
float** Matrix_Jac_Eig(float **array, int n, float *eig)
{
	//先copy一份array在temp_mat中，因为我实在堆区申请的空间,在对其进行处理
	//的过程中会修改原矩阵的值,因此要存储起来,到最后函数返回的
	//时候再重新赋值
	int i, j, flag, k;
	flag = 0;
	k = 0;
	float sum = 0;
	float **temp_mat = (float **)malloc(n * sizeof(float *));
	for (i = 0; i < n; i++)
	{
		temp_mat[i] = (float *)malloc(n * sizeof(float));
	}
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			temp_mat[i][j] = array[i][j];
		}
	}
	//判断是否为对称矩阵
	for (i = 0; i < n; i++)
	{
		for (j = i; j < n; j++)
		{
			if (array[i][j] != array[j][i])
			{
				flag = 1;
				break;
			}
		}
	}
	if (flag == 1)
	{
		printf("error in Matrix_Eig: 输入并非是对称矩阵:\n");
		return NULL;
	}
	else
	{
		//开始执行算法
        //result存放特征向量,result_temp为result的临时存储
        //rot是A1矩阵,array是A矩阵
        //mat旋转矩阵
		int p, q;
		float thresh = 0.0000000001;
		float max = array[0][1];
		float tan_angle, sin_angle, cos_angle;
		float **result = (float **)malloc(n * sizeof(float *));
		if (result == NULL)
		{
			printf("error in result:申请空间失败\n");
			return NULL;
		}
        
		float **result_temp = (float **)malloc(n * sizeof(float *));
		if (result_temp == NULL)
		{
			printf("error in result_temp:申请空间失败\n");
			return NULL;
		}
		float **rot = (float **)malloc(n * sizeof(float *));
		if (rot == NULL)
		{
			printf("error in rot:申请空间失败\n");
			return NULL;
		}
		float **mat = (float **)malloc(n * sizeof(float *));
		if (mat == NULL)
		{
			printf("error in mat:申请空间失败\n");
			return NULL;
		}
		for (i = 0; i < n; i++)
		{
			result[i] = (float *)malloc(n * sizeof(float));
			if (result[i] == NULL)
			{
				printf("error in result:申请子空间失败\n");
				return NULL;
			}
			result_temp[i] = (float *)malloc(n * sizeof(float));
			if (result_temp[i] == NULL)
			{
				printf("error in result_temp:申请子空间失败\n");
				return NULL;
			}
			rot[i] = (float *)malloc(n * sizeof(float));
			if (rot[i] == NULL)
			{
				printf("error in rot:申请子空间失败\n");
				return NULL;
			}
			mat[i] = (float *)malloc(n * sizeof(float));
			if (mat[i] == NULL)
			{
				printf("error in mat:申请子空间失败\n");
				return NULL;
			}
		}
        //初始化特征矩阵result
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					result[i][j] = 1;
				}
				else
				{
					result[i][j] = 0;
				}
			}
		}
        //初始旋转矩阵mat
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					mat[i][j] = 1;
				}
				else
				{
					mat[i][j] = 0;
				}
			}
		}
        //寻找p,q
		max = array[0][1];
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					continue;
				}
				else
				{
					if (fabs(array[i][j]) >= fabs(max))
					{
						max = array[i][j];
						p = i;
						q = j;
					}
					else
					{
						continue;
					}
				}
			}
		}
		while (fabs(max) > thresh)
		{
			if (fabs(max) < thresh)
			{
				break;
			}
            //式(1-1-2)
			tan_angle = -2 * array[p][q] / (array[q][q] - array[p][p]);
			sin_angle = sin(0.5*atan(tan_angle));
			cos_angle = cos(0.5*atan(tan_angle));
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					
					if (i == j)
					{
						mat[i][j] = 1;
					}
					else
					{
						mat[i][j] = 0;
					}
				}
			}
			mat[p][p] = cos_angle;
			mat[q][q] = cos_angle;
			mat[q][p] = sin_angle;
			mat[p][q] = -sin_angle;
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					rot[i][j] = array[i][j];
				}
			}
			for (j = 0; j < n; j++)
			{
				rot[p][j] = cos_angle*array[p][j] + sin_angle*array[q][j];
				rot[q][j] = -sin_angle*array[p][j] + cos_angle*array[q][j];
				rot[j][p] = cos_angle*array[j][p] + sin_angle*array[j][q];
				rot[j][q] = -sin_angle*array[j][p] + cos_angle*array[j][q];
			}
			rot[p][p] = array[p][p] * cos_angle*cos_angle +
				array[q][q] * sin_angle*sin_angle +
				2 * array[p][q] * cos_angle*sin_angle;
			rot[q][q] = array[q][q] * cos_angle*cos_angle +
				array[p][p] * sin_angle*sin_angle -
				2 * array[p][q] * cos_angle*sin_angle;
			rot[p][q] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
				array[p][q] * (2 * cos_angle*cos_angle - 1);
			rot[q][p] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
				array[p][q] * (2 * cos_angle*cos_angle - 1);
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					array[i][j] = rot[i][j];
				}
			}
			max = array[0][1];
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					if (i == j)
					{
						continue;
					}
					else
					{
						if (fabs(array[i][j]) >= fabs(max))
						{
							max = array[i][j];
							p = i;
							q = j;
						}
						else
						{
							continue;
						}
					}
				}
			}
			for (i = 0; i < n; i++)
			{
				eig[i] = array[i][i];
			}
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					sum = 0;
					for (k = 0; k < n; k++)
					{
						sum = sum + result[i][k] * mat[k][j];
					}
					result_temp[i][j] = sum;
				}
			}
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					result[i][j] = result_temp[i][j];
				}
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				array[i][j] = temp_mat[i][j];
			}
		}
		Matrix_Free(result_temp, n, n);
		Matrix_Free(rot, n, n);
		Matrix_Free(mat, n, n);
		Matrix_Free(temp_mat, n, n);
		return result;
	}
}
int Matrix_Free(float **tmp, int m, int n)
{
	//释放动态开辟的内存
	int i, j;
	if (tmp == NULL)
	{
		return(1);
	}
	for (i = 0; i < m; i++)
	{
		if (tmp[i] != NULL)
		{
			free(tmp[i]);
			tmp[i] = NULL;
		}
	}
	if (tmp != NULL)
	{
		free(tmp);
		tmp = NULL;
	}
	return(0);
}

---

三 程序运行结果