前言

在许多实际问题中，所产生的矩阵往往都是对称矩阵，比如我们耳熟能详的实对称矩阵也是重要的研究对象。以下就从实对称矩阵的角度出发，利用特征值的极小极大原理，从普通特征值问题$Ax=\lambda x$衍生到广义特征值问题$Ax=\lambda Bx$逐步讨论其特征值的性质。

【广义特征值问题】设$A=(a_{ij})\in {\mathbb{R}}^{n\times n}$是$n$阶实对称矩阵，$B=(b_{ij})\in {\mathbb{R}}^{n\times n}$是$n$阶实对称正定矩阵，使下式 $$\mathbf{A}\mathbf{x}=\lambda \mathbf{B}\mathbf{x}$$ 有非零解向量$x\in {\mathbb{R}}^n$，则称$\lambda$是矩阵$A$相对于矩阵$B$的特征值，且$x$是属于$\lambda$的特征向量。该问题常见于振动理论。

我们可以发现

• 当$B\ne I$时，该问题是广义特征值问题
• 当$B=I$时，该问题是普通特征值问题

思路：如何利用极小极大原理求第$k$个特征值及奇异值？

利用极大极小原理，我们先确定$n$阶实对称阵的最大最小特征值，然后逐步求第2大和第2小特征值进而归纳到求第$k$大和第$k$小特征值。

本文就对称矩阵特征值的极性与直积做以梳理，完整定理证明请参考西工大的《矩阵论》[1]。

一、实对称矩阵的瑞利商与广义瑞利商性质

我们在讨论实对称矩阵的特征值时，往往会通过实对称阵的瑞利商来研究，因为瑞利商是由如下特征值问题推导出来的，它可以直接求出矩阵的特征值。
$$Ax=\lambda x\Rightarrow x^{T}Ax=\lambda x^{T}x\Rightarrow \lambda =\frac{x^{T}Ax}{x^{T}x}=R(x)$$

【瑞利商定义】设$A=(a_{ij})\in {\mathbb{R}}^{n\times n}$是$n$阶实对称矩阵，$x\in {\mathbb{R}}^n$，则称下式为矩阵$A$的瑞利商($\text{Rayleigh}$商) $$\mathbf{R}(\mathbf{x})=\frac{{\mathbf{x}}^{\mathbf{T}}\mathbf{A}\mathbf{x}}{{\mathbf{x}}^{\mathbf{T}}\mathbf{x}}(x\ne 0)$$

【广义瑞利商定义】设$A=(a_{ij})\in {\mathbb{R}}^{n\times n},B=(b_{ij})\in {\mathbb{R}}^{n\times n}$均是$n$阶实对称矩阵，且$B$正定，$x\in {\mathbb{R}}^n$，则称下式为矩阵$A$相对于矩阵$B$的广义瑞利商$$\mathbf{R}(\mathbf{x})=\frac{{\mathbf{x}}^{\mathbf{T}}\mathbf{A}\mathbf{x}}{{\mathbf{x}}^{\mathbf{T}}\mathbf{B}\mathbf{x}}(x\ne 0)$$

• 【性质1】：$R(x)$是$x$的连续函数
• 【性质2】：$R(x)$是$x$的零次齐次函数（齐次性$R(kx)=R(x)$）
  事实上，对于任意实数$\lambda \ne 0$有下式分别满足齐次性和零次
  $$R(\lambda x)=R(x)={\lambda }^{0}R(x)$$
• 【性质3】：当$x$是由$x_0\ne 0$张成的空间时，$R(x)$是一常数
• 【性质4】：$R(x)$的最大最小值存在，且能够在单位球面 S = { x ∣ x ∈ R n , ∥ x ∥ 2 = 1 } S=\{x|x\in \mathbb{R}^n,\|x\|_2=1\} S={x∣x∈Rn,∥x∥2​=1}上达到
• 【性质5】：非零向量$x_0$是$R(x)$的驻点$\iff x_0$是$Ax=\lambda Bx$的特征向量，当$B=I$时对应于瑞利商问题同理，通过矩阵求导可得

一般情况下，我们令实对称矩阵$A$的特征值按从小到大顺序排列如下
$${\lambda }_1\le {\lambda }_2\le ...\le {\lambda }_n$$
对应标准正交特征向量系为$p_1,p_2,...,p_n$。

【定理】设$A=(a_{ij})\in {\mathbb{R}}^{n\times n}$是$n$阶实对称矩阵，则有 $$\underset{\mathbf{x}\ne 0}{\min }\mathbf{R}(\mathbf{x})={\lambda }_1,\underset{\mathbf{x}\ne 0}{\max }\mathbf{R}(\mathbf{x})={\lambda }_{\mathbf{n}},{\lambda }_1\le \mathbf{R}(\mathbf{x})\le {\lambda }_{\mathbf{n}}$$

【证明】任取$0\ne x\in {\mathbb{R}}^n$，则有
$$x=c_1{p}_1+c_2{p}_2+...+c_n{p}_n(c_1^2+c_2^2+...+c_n^2\ne 0)$$
由于$p_1,p_2,...,p_n$是正交特征向量系，所以有$x_i=c_i{p}_i$
于是有
$$\begin{array}{rl}Ax& =\lambda x={\lambda }_1{c}_1{p}_1+{\lambda }_2{c}_2{p}_2+...+{\lambda }_n{c}_n{p}_n\\ x^{T}Ax& =c_1^2{\lambda }_1+c_2^2{\lambda }_2+...+c_n^2{\lambda }_n\\ x^{T}x& =c_1^2+c_2^2+...+c_n^2\end{array}$$
令$k_i=\frac{c_i^2}{c_1^2+c_2^2+...+c_n^2}$，其中$k_1+k_2+...+k_n=1$，则有
$$R(x)=\frac{x^{T}Ax}{x^{T}x}=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n$$
简单起见，假设$A$是$2$阶实对称阵，即仅有两个特征值${\lambda }_1,{\lambda }_2$满足$R(x)=k_1{\lambda }_1+k_2{\lambda }_2(k_1+k_2=1)$，则如下图所示

从上图，我们可以清晰的看出$R(x)$是$x$的连续函数，该集合也被称为凸包，由此可得
$${\lambda }_1\le R(x)\le {\lambda }_n$$
可以通过如下式子验证$R(p_1)={\lambda }_1$
$$R(p_i)=\frac{p_i^{T}A{p}_i}{p_i^T{p}_i}={\lambda }_i$$
有了$p_k$或$x_k$，我们可以直接求得第$k$小特征值${\lambda }_k$。但问题来了，如果我们不知道$p_k$或者不想依赖于$x_k$，我们如何求得第$k$小特征值${\lambda }_k$呢？这就需要下面一章的极小极大原理了。

【重要推论】若${\lambda }_1=...={\lambda }_k(1\le k\le n)$，则在$\Vert x{\Vert }_2=1$上，$R(x)$的所有极小点为 $${\mathbf{l}}_1{\mathbf{p}}_1+{\mathbf{l}}_2{\mathbf{p}}_2+...+{\mathbf{l}}_{\mathbf{k}}{\mathbf{p}}_{\mathbf{k}}$$ 其中，$l_i\in R(i=1,...,k)$，且满足$l_1^2+l_1^2+..+l_k^2=1$.

二、普通与广义特征值的极小极大原理

由上章，我们得到几个工具，令 V n = span { x 1 , x 2 , . . . , x n }    ( λ 1 ≤ λ 2 ≤ . . . ≤ λ n ) V_n=\text{span}\{x_1,x_2,...,x_n\}\;(\lambda_1 \le \lambda_2 \le... \le \lambda_n ) Vn​=span{x1​,x2​,...,xn​}(λ1​≤λ2​≤...≤λn​)则有
$$R(x)=\frac{x^{T}Ax}{x^{T}x}=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n$$
$${\lambda }_1\le R(x)\le {\lambda }_n\Rightarrow \{\begin{array}{l}{\min }_{x\ne 0,x\in V_n}R(x)={\lambda }_1\\ {\max }_{x\ne 0,x\in V_n}R(x)={\lambda }_n\end{array}$$
当我们想求${\lambda }_2,{\lambda }_{n-1}$时，可以通过缩小张成的子空间得到
$$\begin{gathered} \begin{array}{rl}{\lambda }_2=\underset{x\ne 0}{\min }& R(x)=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n\\ s.t.& k_1=0\end{array} \\ ⋮ \\ \begin{array}{rl}{\lambda }_i=\underset{x\ne 0}{\min }& R(x)=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n\\ s.t.& k_1=k_2=...=k_{i-1}=0\end{array} \end{gathered}$$
同理得
$$\begin{gathered} \begin{array}{rl}{\lambda }_{n-1}=\underset{x\ne 0}{\max }& R(x)=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n\\ s.t.& k_n=0\end{array} \\ ⋮ \\ \begin{array}{rl}{\lambda }_{n-i-1}=\underset{x\ne 0}{\min }& R(x)=k_1{\lambda }_1+k_2{\lambda }_2+...+k_n{\lambda }_n\\ s.t.& k_n=k_{n-1}=...=k_{n-i}=0\end{array} \end{gathered}$$
因此，我们可以归纳出如下定理

【定理】设$x\in L(p_r,p_{r+1},...,p_s),1\le r\le s\le n$，则有 $$\underset{\mathbf{x}\ne 0}{\min }\mathbf{R}(\mathbf{x})={\lambda }_{\mathbf{r}}\underset{\mathbf{x}\ne 0}{\max }\mathbf{R}(\mathbf{x})={\lambda }_{\mathbf{s}}$$

2.1 引出问题：由于$V_k$不唯一导致得到多个特征值

但以上定理在$p_r,p_s$未知下无法使用，因此我们不再指定让某个系数$k_i=0$，而是选取$k$维子空间$V_k$来求，由于$V_k$是不唯一的，因此可能会得到多个特征值，例如我们想要得到${\lambda }_2$，则选取$V_{n-1}$，有如下两种情况

$$\underset{x\ne 0}{\min }R(x)=\{\begin{array}{l}{\lambda }_1\text{if}{x}_1\in V_{n-1}\\ {\lambda }_2\text{if}{x}_1\notin V_{n-1}\end{array}$$
$$\underset{x\ne 0}{\max }R(x)=\{\begin{array}{l}{\lambda }_n\text{if}{x}_n\in V_{n-1}\\ {\lambda }_{n-1}\text{if}{x}_n\notin V_{n-1}\end{array}$$

2.2 解决问题：使用极大极小原理固定特征向量

对于上述子空间$V_k$不唯一情况，得到
$$\underset{0\ne x\in V_{n-1}}{\min }R(x)\le {\lambda }_2\underset{0\ne x\in V_{n-1}}{\max }R(x)\ge {\lambda }_{n-1}$$
为解决此问题，我们使用极小极大原理得到
$${\lambda }_2=\underset{V_{n-1}}{\max }\left[\underset{0\ne x\in V_{n-1}}{\min }R(x)\right],{\lambda }_{n-1}=\underset{V_{n-1}}{\min }\left[\underset{0\ne x\in V_{n-1}}{\max }R(x)\right]$$
为此，我们归纳出一般的式子，我们

【定理】设$V_k$是${\mathbb{R}}^n$中的任意一个$k$维子空间，则普通特征值问题与广义特征值问题从小到大的第$k$个特征值和$n-(k-1)$个特征值具有如下极小极大性质
$${\lambda }_{\mathbf{n}-(\mathbf{k}-1)}=\underset{{\mathbf{V}}_{\mathbf{k}}}{\max }\left[\underset{0\ne \mathbf{x}\in {\mathbf{V}}_{\mathbf{k}}}{\min }\mathbf{R}(\mathbf{x})\right],{\lambda }_{\mathbf{k}}=\underset{{\mathbf{V}}_{\mathbf{k}}}{\min }\left[\underset{0\ne \mathbf{x}\in {\mathbf{V}}_{\mathbf{k}}}{\max }\mathbf{R}(\mathbf{x})\right]$$

• 左式被称为特征值的极大极小原理
• 右式被称为特征值的极小极大原理

三、矩阵奇异值的极小极大性质

我们通过矩阵瑞利商的极小极大原理，可以衍生到解决奇异值问题，我们将矩阵$A\in {\mathbb{R}}_r^{m\times n}$的奇异值排列如下 [其中，${\sigma }_i=\sqrt{{\lambda }_i(A^{T}A)}$]
$$0={\sigma }_1={\sigma }_2=...={\sigma }_{n-r}\le {\sigma }_{n-r+1}\le ...\le {\sigma }_n$$

我们令$B=A^{T}A$，则实对称矩阵$B$的瑞利商如下
$$R(x)=\frac{x^{T}Bx}{x^{T}x}=\frac{x^T(A^{T}A)x}{x^{T}x}=\frac{(Ax{)}^{T}Ax}{x^{T}x}=\frac{\Vert Ax{\Vert }_2^2}{\Vert x{\Vert }_2^2}=\lambda =\sqrt{\sigma }$$
则矩阵$A$的第$k$个奇异值和第$n-k+1$个奇异值具有如下极小极大性质
$${\sigma }_{n-(k-1)}=\underset{V_k}{\max }\left[\underset{0\ne x\in V_k}{\min }\frac{\Vert Ax{\Vert }_2}{\Vert x{\Vert }_2}\right],{\sigma }_k=\underset{V_k}{\min }\left[\underset{0\ne x\in V_k}{\max }\frac{\Vert Ax{\Vert }_2}{\Vert x{\Vert }_2}\right]$$
其中，$V_k$是${\mathbb{R}}^n$中的任意一个$k$维子空间。

附录：矩阵直积($\text{Kronecker}$积)的概念

运用矩阵的直积运算，能够将线性矩阵方程转换为线性代数方程组进行求解

【定义】设$A=(a_{ij})\in {\mathbb{C}}^{m\times n},B=(b_{ij})\in {\mathbb{C}}^{p\times q}$，则称如下分块矩阵为$A$与$B$的直积($\text{Kronecker}$积)

参考文献

程云鹏, 凯院, 仲. 矩阵论[M]. 西北工业大学出版社, 2006.

一 算法原理

雅可比方法用于求解实对称矩阵的特征值和特征向量,对于实对称矩阵$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$.

二 C语言实现

#include<stdio.h>
#include<stdlib.h>
float** Matrix_Jac_Eig(float **array, int n, float *eig);
int Matrix_Free(float **tmp, int m, int n);
int main(void)
{
	int n;
	printf("请输入矩阵维度:\n");
	scanf("%d", &n);
	float **array = (float **)malloc(n * sizeof(float *));
	if (array == NULL)
	{
		printf("error :申请数组内存空间失败\n");
		return -1;
	}
	for (int i = 0; i < n; i++)
	{
		array[i] = (float *)malloc(n * sizeof(float));
		if (array[i] == NULL)
		{
			printf("error :申请数组子内存空间失败\n");
			return -1;
		}
	}
	printf("请输入矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			scanf("%f", &array[i][j]);
		}
	}
	float *eig = (float *)malloc(n * sizeof(float));
	float **Result = Matrix_Jac_Eig(array, n, eig);
	printf("特征矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%f ", Result[i][j]);
		}
		printf("\n");
	}
	printf("特征矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		printf("%f \n", eig[i]);
	}
	Matrix_Free(Result, n, n);
	free(eig);
	eig = NULL;
	return 0;
}
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
	{
		//开始执行算法
		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 Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **result_temp = (float **)malloc(n * sizeof(float *));
		if (result_temp == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **rot = (float **)malloc(n * sizeof(float *));
		if (rot == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **mat = (float **)malloc(n * sizeof(float *));
		if (mat == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		for (i = 0; i < n; i++)
		{
			result[i] = (float *)malloc(n * sizeof(float));
			if (result[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			result_temp[i] = (float *)malloc(n * sizeof(float));
			if (result_temp[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			rot[i] = (float *)malloc(n * sizeof(float));
			if (rot[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			mat[i] = (float *)malloc(n * sizeof(float));
			if (mat[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					result[i][j] = 1;
				}
				else
				{
					result[i][j] = 0;
				}
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					mat[i][j] = 1;
				}
				else
				{
					mat[i][j] = 0;
				}
			}
		}
		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;
			}
			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);
}

三 结果

1. 2阶实对称矩阵特性

定理：2阶实对称矩阵H的特征值是实数

H=[a,b;b,c]    

    a,b,c是实数，λ 是特征值
A=[a-λ,b;b,c-λ]    

    特征值求解方法为：(a- λ )(c- λ) - b2  = 0

    求解方程得到两个根为：
λ=（a+c）±（a+c）2-4（ac-b2）2  

                   （a+c）2-4ac-b2=a-c2+4b2≥0

    所以，在a、b、c为实数时，特征值也是实数。

2、特征向量

根据特征值和特征向量的定义：HX=λX，（H-λE）X = 0；因此方程若有解，则

det（H-λE）=0；

设
A=[a-λ,b;b,c-λ]    

则有-b/(a-λ) = (c-λ)/b,   线性齐次方程组AX=0有非零解，其中之一解向量 [1,-b/(a-λ)],归一化后得到标准解。